# Continuous Urban Change Detection from Satellite Image Time Series with
Temporal Feature Refinement and Multi-Task Integration

Sebastian Hafner, Heng Fang, Hossein Azizpour and Yifang Ban*, Senior Member,
IEEE The authors are with the Division of Geoinformatics (Sebastian Hafner and
Yifang Ban) and the Division of Robotics, Perception and Learning (Heng Fang
and Hossein Azizpour) at KTH Royal Institute of Technology, 114 28 Stockholm,
Sweden.*Corresponding author: Yifang Ban (e-mail<EMAIL_ADDRESS>
###### Abstract

Urbanization advances at unprecedented rates, resulting in negative effects on
the environment and human well-being. Remote sensing has the potential to
mitigate these effects by supporting sustainable development strategies with
accurate information on urban growth. Deep learning-based methods have
achieved promising urban change detection results from optical satellite image
pairs using convolutional neural networks (ConvNets), transformers, and a
multi-task learning setup. However, transformers have not been leveraged for
urban change detection with multi-temporal data, i.e., $>2$ images, and multi-
task learning methods lack integration approaches that combine change and
segmentation outputs. To fill this research gap, we propose a continuous urban
change detection method that identifies changes in each consecutive image pair
of a satellite image time series. Specifically, we propose a temporal feature
refinement (TFR) module that utilizes self-attention to improve ConvNet-based
multi-temporal building representations. Furthermore, we propose a multi-task
integration (MTI) module that utilizes Markov networks to find an optimal
building map time series based on segmentation and dense change outputs. The
proposed method effectively identifies urban changes based on high-resolution
satellite image time series acquired by the PlanetScope constellation (F1
score 0.551) and Gaofen-2 (F1 score 0.440). Moreover, our experiments on two
challenging datasets demonstrate the effectiveness of the proposed method
compared to bi-temporal and multi-temporal urban change detection and
segmentation methods. Code is available on GitHub:
https://github.com/SebastianHafner/ContUrbanCD.

###### Index Terms:

Earth observation, Remote sensing, Multi-temporal, Multi-task learning,
Transformers

##  I Introduction

Urbanization is progressing at unprecedented rates [1]. Thus, the global
amount of urban land is projected to increase by a factor of 2–6 over the 21st
century [2]. The rapid expansion of urban land, i.e., urban sprawl, is
associated with multiple negative effects on the environment and human well-
being [3, 4]. To mitigate urban sprawl, informed and sustainable urban
development strategies are crucial [5]. However, these strategies are
currently hampered by a lack of timely information on the extent of urban
land.

Remote sensing is an efficient tool to monitor the Earth’s surface [6]. Urban
changes are commonly detected from two satellite images acquired at different
times over the same geographical area. Traditional change detection methods
use arithmetic operations to derive change features from bi-temporal image
pairs. For example, various arithmetic methods have been developed to derive
change features from optical images, such as image differencing, image
regression, and change vector analysis [7]. These features are then classified
into changed/unchanged pixels or objects using different classification
algorithms, including machine learning algorithms [7, 6].

In recent years, deep learning has been continuously replacing traditional
change detection methods [8, 9]. Specifically, deep convolutional neural
networks (ConvNets) have been used extensively for change detection in bi-
temporal optical satellite image pairs (see Figure 1(a)). The simplest way of
adapting common ConvNets such as U-Net [10] or HRNet [11] for change detection
is with an input-level fusion (or early fusion [12]) strategy, referring to
the concatenation of image pairs before passing them to a ConvNet. Contrarily,
late fusion strategies typically process images separately in a Siamese
network consisting of two ConvNets with shared weights. Extracted bi-temporal
features are then fused using concatenation or absolute differencing [12, 13].
Since Siamese networks are generally considered preferable to input-level
fusion strategies, multiple studies developed modules that are incorporated
into Siamese network architectures to improve feature representations [14, 15,
16]. For very high-resolution (VHR) imagery, Chen et al. [17, 18] proposed to
refine features extracted by ConvNets using a transformer-based module that
alleviates the limited long-range context modeling capability of convolutions
with self-attention. Since then, self-attention has become a popular mechanism
for capturing long-range contextual information in VHR change detection [19,
16, 20, 21].

((a)) Bi-temporal urban change detection

((b)) Multi-temporal urban change detection

((c)) Continuous urban change detection using a bi-temporal model

((d)) Continuous urban change detection using the proposed model

Figure 1: Overview of standard urban change detection frameworks and the
proposed method. (a) Bi-temporal urban change detection is typically performed
on image pairs acquired multiple years apart. (b) Multi-temporal change
detection methods leverage image time series but only predict changes between
the first and last image. (c) A bi-temporal model is applied to consecutive
image pairs of a time series to perform continuous urban change detection;
however, this method fails to incorporate multi-temporal (i.e., $>2$ images)
information (c). On the other hand, we propose (d) a continuous urban change
detection method that incorporates multi-temporal information.

In recent years, high-resolution (i.e., 1 – 10 m) satellite image time series
have become increasingly available [22]. Those data have enabled a shift from
detecting land cover changes in image pairs acquired years apart to continuous
annual and sub-annual change detection [23]. In contrast, urban change
detection methods are predominantly designed for bi-temporal change detection
from image pairs acquired multiple years apart [24]. However, considering the
unprecedented rate of global urbanization [25], it is essential to develop a
new suite of methods that detect urban changes continuously. While continuous
change detection can be achieved by applying a bi-temporal model to
consecutive image pairs in satellite image time series (see Figure 1(c)), this
approach fails to exploit multi-temporal, i.e., $>$ 2 images, information.
Furthermore, recent studies [26, 27] have demonstrated the effectiveness of
multi-temporal change detection models which predict changes between the first
and last images of a satellite image time series (see Figure 1(b)). For
example, multi-temporal information can help to reduce commission errors from
registration errors, illumination differences, or other types of change
unrelated to the problem of interest [27]. In addition, it can mitigate the
effect of cloud artifacts in single images [26]. Existing multi-temporal
change detection methods employ either recurrent layers, such as long short-
term memory (LSTM) layers [28], or 3D convolutional layers to model temporal
information [26, 27, 29, 30]. While these layers effectively model short-range
temporal dependencies in time series data, the self-attention mechanism can
explicitly model temporal dependencies across all timestamps of a time series
[31]. However, self-attention is predominately used in urban change detection
to model long-range dependencies in space [17, 18, 19, 16]. In contrast, this
research aims to use self-attention to model multi-temporal information
extracted from satellite image time series for urban change detection.

Another promising avenue of research for change detection is multi-task
learning [32], where a related semantic segmentation task is trained parallel
to the change detection task using a shared feature representation. The change
detection task is commonly combined with building segmentation [33, 27, 34,
35, 36] or land cover segmentation [37, 38, 39]. To that end, Siamese networks
are extended with an additional decoder for the semantic segmentation task.
The feature maps extracted by the encoder are then shared between the change
decoder and the segmentation decoder. However, despite the attention multi-
task learning has attracted in change detection, effective methods to
integrate segmentation and change outputs have been largely unaddressed. For
example, multi-task binary change detection studies consider building and
change predictions independent outputs of the network [33, 27, 34, 35, 36]. On
the other hand, semantic change detection studies, concerned with identifying
not only whether changes occurred (i.e., binary change detection) but also
what types of changes occurred, apply a simple masking approach for multi-task
integration [37, 38, 39]. Specifically, land cover predictions and binary
change predictions are integrated by masking land cover changes derived from
the bi-temporal semantic predictions with the binary change prediction.
Therefore, existing multi-task change detection studies do not account for
inconsistencies between the segmentation and change predictions. Moreover,
these approaches fail to exploit the complementary information produced by
multi-task predictions.

In this paper, we propose a continuous urban change detection method (see
Figure 1(d)) and explore two research gaps in the current literature, namely
(1) the use of the self-attention mechanism to model multi-temporal
information in urban change detection and (2) the integration of segmentation
and change predictions in multi-task learning setups. Specifically, we propose
a new network architecture that relies on convolutions to extract multi-
temporal building representations and employs self-attention to model temporal
dependencies in feature space. We also propose a novel integration approach
that determines the most likely state for each timestamp in a time series by
considering multi-task predictions using Markov networks. The effectiveness of
the proposed architecture and integration approach is demonstrated on two
urban change detection datasets featuring high-resolution optical satellite
image time series, namely the SpaceNet 7 dataset [40] and the Wuhan Urban
Semantic Understanding (WUSU) dataset [41].

The following summarizes the main contributions of this paper.

  * •

We present a transformer-based feature refinement module that effectively
models temporal information in satellite image time series, leading to
improved building representations.

  * •

We introduce a new multi-task integration approach that represents
segmentation and change outputs in Markov networks to find the most likely
building map for each image in a time series using maximum a posteriori (MAP)
estimation.

  * •

Experiments on two datasets, namely SpaceNet 7 and WUSU, show that the
proposed continuous urban change detection method is more effective than
related methods.

##  II Related Work

###  II-A Bi-temporal change detection

In recent years, a plethora of deep learning-based bi-temporal change
detection methods have been proposed. Most of these works focus on developing
new Siamese network architectures and/or training strategies. Initially, Daudt
et al. [12, 13] proposed two Siamese ConvNet architectures for change
detection. The Siam-Diff and Siam-Conc architectures employ encoders with
shared weights for feature extraction from bi-temporal high-resolution image
pairs and combine the corresponding feature maps using a subtraction and
concatenation strategy, respectively. While encoders and decoders in these
models follow the U-Net architecture [10], Fang et al. [14] incorporated a
nested U-Net (i.e., UNet++ [42]) into a Siamese network to maintain high-
resolution fine-grained representations through dense skip connections. Many
works also improved Siamese networks by incorporating different modules into
the architecture. For example, an ensemble channel attention module was
proposed for feature refinement in [14], and a new spatial pyramid pooling
block was utilized in [15] to preserve shapes of change areas.

However, most recent methods are developed for bi-temporal change detection
from VHR image pairs. Consequently, many methods employ the self-attention
mechanism to improve the modeling of long-range dependencies in VHR imagery
[17, 18, 16, 19]. Both [17] and [18] extract image features with ConvNets and
employ self-attention modules to learn more discriminative features. Other
works combined ConvNets and transformers with attention modules and multi-
scale processing [16, 43]. Bandara et al. [19], on the other hand, proposed a
fully transformer-based change detection method. Specifically, ChangeFormer
combines two hierarchically structured transformer encoders with shared
weights and a multi-layer perception decoder in a Siamese network
architecture. Since transformer-based methods strongly rely on pretraining,
Noman et al. [21] recently proposed ScratchFormer which is a transformer-based
change detection method that is trained from scratch but achieves SOTA
performance. The ScrachFormer architecture utilizes shuffled sparse attention
layers that enable faster convergence due to their sparse structure. Although
these transformer-based methods are considered SOTA for urban change
detection, their effectiveness has been predominately demonstrated on bi-
temporal VHR datasets such as LEVIR-CD [17] and WHU-CD [44]. In comparison,
high-resolution imagery is acquired much more frequently by satellite
constellations such as PlanetScope, making it an invaluable data source for
change detection applications. Therefore, developing methods that effectively
leverage transformers for change detection from high-resolution imagery is
crucial.

###  II-B Change detection from satellite image time series

Few studies have developed deep learning methods for urban change detection
from high-resolution satellite image time series. For example, Papadomanolaki
et al. [26] proposed to incorporate LSTM networks into a U-Net model to
leverage optical satellite image time series for change detection. Their
L-UNet outperformed bi-temporal ConvNet-based methods on a bi-temporal dataset
enriched with intermediate satellite images [26]. Others proposed an encoder-
decoder LSTM model that is trained to rearrange temporally shuffled time
series [30]. The core assumption of this unsupervised method is that the model
fails to correctly rearrange shuffled data for changed pixels. Finally,
Meshkini et al. [29] recently proposed a weakly supervised change detection
method that employs 3D convolutional layers to capture spatial-temporal
information in satellite image time series. Despite these promising results
underlying the value of time series information, multi-temporal urban change
detection methods may benefit from layers with better temporal modeling
capabilities. Furthermore, most of these methods fail to detect urban changes
continuously and are limited to the detection of changes between the first and
the last image (i.e., multi-temporal change detection in Figure1(b)).

###  II-C Multi-task learning

Multi-task learning has been investigated by several studies for urban change
detection over the past years. Liu et al. [33] proposed a dual-task Siamese
ConvNet to learn more discriminative feature representations for building
change detection from bi-temporal image pairs. The proposed dual-task
constrained deep Siamese convolutional network (DTCDSCN) consists of three
main components: a shared ResNet-based encoder, a shared decoder for building
segmentation, and a separate decoder for change detection. On the other hand,
Papadomanolaki et al. [27] proposed a multi-task learning framework for urban
change detection from image time series by adding building segmentation tasks
for the first and last images of a time series to the urban change detection
task. While L-UNet [26] is employed to extract changes, the segmentation is
performed with a separate decoder that directly uses the feature maps
extracted for the image pair by the shared encoder.

Some urban change detection studies also combined multi-task learning with
semi-supervised learning [34, 35]. In [34], the Siam-Diff network [13] was
extended with an additional shared decoder for building segmentation, and an
unsupervised term was introduced to encourage consistency between the changes
derived from the building predictions and those predicted by the change
decoder. Shu et al. [35], on the other hand, proposed to learn consistency
between two building predictions corresponding to the pre-change image. The
first prediction is obtained by segmenting the pre-change image and the second
one by combining segmentation features of the post-change image with changes
features.

Apart from binary change detection, multi-task learning has been employed for
semantic change detection [37, 38, 39] and 3D change detection [20]. Daudt et
al. [37] proposed a learning strategy that integrates change detection and
land cover mapping from bi-temporal image pairs. Specifically, they propose a
network consisting of separate encoder-decoder sub-networks for land cover
mapping and binary change detection, where the change detection decoder also
receives features from the mapping encoder in the form of difference skip
connections. ChangeMask [38] uses a similar setup but introduces a transformer
applied to the bi-temporal feature maps to learn more robust representations
for the change task. Furthermore, [39] optimized semantic features in
ChangeMask by adding a multi-task join loss (i.e., ChangeMJ). Finally,
Marsocci et al. [20] extended the bi-temporal image transformer (BIT)
presented in [18] with a regression head to solve 2D and 3D change detection
simultaneously using shared feature representations.

##  III The proposed method

###  III-A Overview

Figure 2: Overview of the proposed method. First, an encoder extracts multi-
scale feature maps from a satellite image time series. Next, transformer-based
temporal feature refinement (TFR) modules enrich feature maps at each scale
with multi-temporal information, and the change feature (CF) module generates
bi-temporal difference feature maps from the temporally refined feature maps.
Then, two separate decoders are used to obtain segmentation and change
predictions from the respective feature maps. Finally, predictions for the two
tasks are combined using a multi-task integration (MTI) module. For brevity,
only three out of the five scales of the feature maps are shown.

As illustrated in Figure 2, the main components of the proposed method are a
ConvNet-based encoder, transformer-based temporal feature refinement (TFR)
modules, a change feature (CF) module, two task-specific ConvNet-based
decoders, and a Markovian multi-task (MTI) integration module. The following
summarizes the urbanization monitoring process of the proposed method for a
time series of satellite images:

  1. 1.

First, for each image in the time series, multi-scale feature maps are
extracted using an encoder with shared weights.

  2. 2.

Next, the above time series of feature maps are grouped by scale and fed to
separate TFR modules consisting of transformer encoder layers. The temporally
refined feature maps are regrouped according to their timestamp.

  3. 3.

Then, the CF module obtains change feature maps from the temporally refined
segmentation feature maps. The module considers changes between all possible
combinations of temporal pairs.

  4. 4.

Two task-specific decoders are deployed to obtain building segmentation
outputs for each image in the time series from the temporally refined
segmentation features maps and change outputs from the change feature maps.

  5. 5.

Finally, the building and urban change outputs are combined using the MTI
module. The module uses pixel-wise Markov networks to obtain optimal building
states for the satellite image time series.

Detailed descriptions of the components comprising the proposed method, as
well as the loss function, are given in the following sections.

###  III-B Satellite image time series encoding

We consider a time series of $T$ satellite images, represented as
$\mathbf{X}\in\mathbb{R}^{T\times C\times H\times W}$, where $C$, $H$, and $W$
denote the channel, height, and width dimensions, respectively. A ConvNet
encoder with shared weights is utilized to separately extract feature maps
$\mathbf{F}_{\rm seg}$ from each image in the time series, as follows:

$\mathbf{F}_{\rm seg}=e(\mathbf{X}),$ (1)

where $e(.)$ represents the encoder, and subscript ${\rm seg}$ indicates that
the feature maps contain representations for building segmentation.

The architecture of the encoder is based on the U-Net encoder [10].
Specifically, after an initial convolution block, the combination of a max-
pooling layer and a consecutive convolution block is applied four times. Each
of these four steps halves the spatial dimensions $H$ and $W$ due to the
pooling operation, whereas the number of features $D$ is doubled with the
convolution block. Importantly, U-Net achieves precise localization by
leveraging skip connections that forward the feature maps before each pooling
operation to the decoder. Therefore, the output of the encoder consists of
five feature map time series with different scales. To denote the scale of
these time series, we introduce superscript $s$ in the notation:
$\mathbf{F}_{\rm seg}^{s}$, where $s\in\\{0,1,2,3,4\\}$. For a given feature
map time series $\mathbf{F}_{\rm seg}^{s}$, the sizes of the height and width
dimensions, as well as the feature dimension, are dependent on $s$, as
follows:

$H^{s}=\frac{H}{2^{s}},\thinspace W^{s}=\frac{W}{2^{s}},\thinspace
D^{s}=64\cdot 2^{s}.$ (2)

It should be noted that for brevity, Figure 2 illustrates the proposed method
for $s\in\\{0,1,2\\}$.

###  III-C Temporal feature refinement

The TFR module, illustrated in Figure 3, creates temporally refined feature
maps using the self-attention mechanism along the temporal dimension [31].
Specifically, the module takes as input a time series of feature maps at the
same scale $s$ and reshapes this 4D tensor to a 3D tensor
$\mathbf{T}^{s}\in\mathbb{R}^{T\times D^{s}\times P^{s}}$ by flattening the
spatial dimensions $H^{s}$ and $W^{s}$. Consequently, $T$, $D$, and $P$
represent the temporal, feature embedding, and spatial dimensions,
respectively.

After reshaping the feature map time series, self-attention is applied along
the temporal dimension $T$ for each cell in the spatial dimension $P$.
However, since the self-attention mechanism contains no recurrence, it is
necessary to first inject information about the temporal position of the
feature vectors in the time series. Specifically, temporal encodings having
the same dimension as the feature vectors are generated based on sine and
cosine functions of different frequencies [31]. These relative temporal
encodings are then added to the feature vectors.

The tensor, enriched with relative temporal position information, is passed
through two transformer encoder layers (see Figure 4). The key component of
the transformer encoder layer is the multi-head attention block, which
performs self-attention defined as follows:

$\text{Att}(\mathbf{Q},\mathbf{K},\mathbf{V})=\text{softmax}\left(\frac{\mathbf{Q}\mathbf{K}^{T}}{\sqrt{D}}\right)\mathbf{V},$
(3)

where $\mathbf{Q}$, $\mathbf{K}$, and $\mathbf{V}$ are referred to as query,
key, and value, respectively. The query-key-value triplet is computed with
three linear projection layers with parameter matrices
$\mathbf{W}^{Q},\mathbf{W}^{K},\mathbf{W}^{V}\in\mathbb{R}^{D\times D}$ that
are separately applied to a given cell of the 3D tensor $\mathbf{T}^{s}_{p}$,
where $p$ denotes the cell index in $\mathbf{T}^{s}$.

The core idea of multi-head attention is, however, that self-attention is
performed multiple times in parallel using $h$ attention heads, as follows:

$\begin{split}\text{MultiHead}(\mathbf{Q},\mathbf{K},\mathbf{V})&=\text{Concat}\left(\text{head}_{1},...,\text{head}_{h}\right)\mathbf{W}^{O},\\\
\text{where
}\text{head}_{i}&=\text{Att}(\mathbf{Q}\mathbf{W}_{i}^{Q},\mathbf{K}\mathbf{W}_{i}^{K},\mathbf{V}\mathbf{W}_{i}^{V}).\end{split}$
(4)

Each $\text{head}_{i}$ performs self-attention on different projections of the
keys, values, and queries obtained from linear layers with parameter matrices
$\mathbf{W}_{i}^{Q},\mathbf{W}_{i}^{K},\mathbf{W}_{i}^{V}\in\mathbb{R}^{D\times
D_{\rm head}}$. Finally, the concatenated outputs of the heads are reprojected
using parameter matrix $\mathbf{W}^{Q}\in\mathbb{R}^{hD_{\rm head}\times D}$.
We employ $h=2$ attention heads, where the head dimension is given by $D_{\rm
head}=D/h$.

After applying self-attention to each multi-temporal feature vector, we obtain
a 3D tensor containing temporally refined building representations. In
practice, however, all cells of tensor $\mathbf{T}^{s}$ are processed in
parallel by incorporating them into the batch dimension which is omitted for
clarity. Finally, the 3D tensor is reshaped to the dimensions of the feature
map time series by unflattering dimension $P$. We denote this temporally
refined feature map time series with $\mathbf{\hat{F}}^{s}_{\rm seg}$.

Figure 3: Illustration of the temporal feature refinement (TFR) module. Figure
4: Illustration of the transformer encoder layer applying multi-headed self-
attention.

###  III-D Change feature extraction

The CF module is used to convert the temporally refined segmentation feature
maps to change feature maps. Specifically, we consider the temporally refined
feature maps $\mathbf{\hat{F}}_{\rm seg}^{t}$ and $\mathbf{\hat{F}}_{\rm
seg}^{k}$, corresponding to two images acquired at time $t$ and time $k$,
where $1\leq t<k\leq T$. Then, change feature map $\mathbf{F}_{\rm ch}^{n}$
corresponding to the urban changes between the bi-temporal image pair is
constructed as follows:

$\mathbf{F}_{\rm ch}^{n}=\mathbf{\hat{F}}_{\rm seg}^{k}-\mathbf{\hat{F}}_{\rm
seg}^{t},$ (5)

where $n$ denotes a change edge between timestamps $t$ and $k$. It should be
noted that this is done for each scale $s$ of the feature maps.

This operation is identical to the change feature computation in the Siam-Diff
method [13]. However, instead of only considering changes between the first
and the last images of a time series, the CF module computes change feature
maps for all possible combinations of image pairs. The total number of
combinations $N_{\rm dense}$ is given by the length of the time series $T$,
defined as follows:

$\displaystyle N_{\rm dense}=$ $\displaystyle\frac{T(T-1)}{2}.$ (6)

We refer to the scenario when all possible combinations of image pairs are
considered as dense. However, we also investigate sparser settings such as
adjacent, considering only changes between temporally adjacent images ($N_{\rm
adjacent}=T-1$), and cyclic, adding the changes between the first and last
image to the adjacent setting ($N_{\rm cyclic}=T$) (see Figures 5(a) – 5(c)
for visualizations of these settings in Markov networks with $T=5$). We
provide ablation results on the different edge settings in Table III. Finally,
it should be noted that the CF module does not have any trainable parameters.

###  III-E Multi-task decoding

Two separate decoders are deployed to convert the temporally enriched
segmentation feature maps and the change feature maps to building outputs and
urban change outputs, respectively. Formally, we obtain $T$ built-up area
segmentation outputs $\mathbf{O}_{\rm seg}\in(0,1)^{T\times H\times W}$ with
the segmentation decoder $d_{\rm seg}(.)$ as follows:

$\mathbf{O}_{\rm seg}=d_{\rm seg}(\mathbf{\hat{F}}_{\rm seg}).$ (7)

Furthermore, we obtain $N$ change detection outputs $\mathbf{O}_{\rm
ch}\in(0,1)^{N\times H\times W}$ with the change decoder $d_{\rm ch}(.)$ as
follows:

$\mathbf{O}_{\rm ch}=d_{\rm ch}(\mathbf{F}_{\rm ch}).$ (8)

Both decoders follow the architecture of the U-Net expansive path consisting
of four upsampling blocks followed by a $1\times 1$ convolution layer and a
sigmoid activation function. Upsampling blocks double the height and width of
feature maps via a transpose conv $2\times 2$ layer. Upsampled feature maps
are then concatenated along the channel dimension with the temporally refined
feature maps matching their scale (skip connection). Subsequently, the layer
triplet $3\times 3$ convolution, batch normalization, and ReLu activation is
applied twice.

###  III-F Loss function

The network is trained using a loss function composed of two terms, namely for
the urban change detection task ($\mathcal{L}_{\rm ch}$) and the building
segmentation task ($\mathcal{L}_{\rm seg}$). For both loss terms, a Jaccard
metric measuring the similarity between continuous network outputs
$\mathbf{O}\in(0,1)$ and binary labels $\mathbf{Y}\in\\{0,1\\}$ is used [45].
We denote the Jaccard metric by $J(,)$ and define the loss function as
follows:

$\mathcal{L}=\sum_{t=1}^{T}J(\mathbf{O}^{t}_{\rm seg},\mathbf{Y}^{t}_{\rm
seg})+\sum_{n=1}^{N}J(\mathbf{O}_{\rm ch}^{n},\mathbf{Y}_{\rm ch}^{n}),\\\ $
(9)

where $T$ denotes the length of the time series and $N$ denotes the number of
edges (i.e., combinations of bi-temporal image pairs) considered. Segmentation
and change labels are denoted by $\mathbf{Y}_{\rm seg}\in\\{0,1\\}^{T\times
H\times W}$ and $\mathbf{Y}_{\rm ch}\in\\{0,1\\}^{N\times H\times W}$,
respectively. Specifically, we assume that pixel-wise building annotations,
$\mathbf{Y}_{\rm seg}$, are available and derive pixel-wise built-up changes,
$\mathbf{Y}_{\rm ch}$, according to the considered edges.

###  III-G Multi-task integration

To combine segmentation and change predictions, we propose the MTI module
which determines the optimal building segmentation output for each image in a
time series. Since this is a pixel-based approach, we represent the location
of a specific pixel in the segmentation and change output notations by
introducing superscript coordinates $i$ and $j$. Following that,
$\mathbf{O}_{\rm seg}^{(i,j),t}$ denotes the segmentation output for a
specific pixel $i,j$ at timestamp $t$, where $i\in\\{1,\ldots,H\\}$,
$j\in\\{1,\ldots,W\\}$, and $t\in\\{1,\ldots,T\\}$. Likewise, the change
output for a specific pixel is denoted by $\mathbf{O}_{\rm ch}^{(i,j),n}$,
where $n$ denotes the change edge connecting timestamps $t$ and $k$.

The core idea of the MTI module is to represent segmentation and change
outputs in a pairwise Markov network. This subclass of Markov networks is
associated with an undirected graph $G=(\mathcal{N},\mathcal{E})$ in which the
nodes $\mathcal{N}$ correspond to random variables and the edges $\mathcal{E}$
represent pairwise relationships between the nodes (see [46]). For a given
pixel, we construct a Markov network with $T$ nodes corresponding to the
timestamps in an image time series. Specifically, all nodes in the network
correspond to a binary variable representing the presence of buildings (i.e.,
$\mathcal{N}^{t}\in\\{true,false\\}$). We use state 1 to denote $true$,
representing the presence of a building, and state 0 to denote $false$,
representing the absence of a building. Adjacent nodes in the network are
connected with $N-1$ edges, where we use $\mathcal{E}^{t\rightleftharpoons k}$
to denote an edge connecting timestamps $t$ and $k$. We refer to this Markov
network structure as an adjacent network (see Figure 5(a)).

To represent the segmentation and change outputs, the graph structure needs to
be associated with a set of parameters that capture the relationships between
nodes. The parameterization in a pairwise Markov network is achieved by
assigning factors over nodes or edges, where a factor $\phi$, also referred to
as potential, is a function from value assignments of random variables to real
positive numbers $\mathbb{R}^{+}$. Thus, a pairwise Markov network is
associated with a set of node potentials
$\\{\phi(\mathcal{N}_{t}):t=1,...,T\\}$ and a set of edge potentials
$\\{\phi(\mathcal{N}_{t},\mathcal{N}_{k}):(\mathcal{N}_{t},\mathcal{N}_{k})\in
G\\}$. The overall distribution represented by the network is then the
normalized product of all the node and edge potentials.

The segmentation outputs for a specific pixel are incorporated into the Markov
network by assigning a factor $\phi_{t}$ over each node $\mathcal{N}_{t}$.
Then, the segmentation outputs are encoded as node potentials, as follows:

$\begin{split}\phi_{t}(\mathcal{N}^{t}=1)&=\mathbf{O}_{\rm seg}^{(i,j),t}\\\
\phi_{t}(\mathcal{N}^{t}=0)&=1-\mathbf{O}_{\rm seg}^{(i,j),t}.\end{split}$
(10)

These node potentials characterize the bias of nodes towards state 1 or 0,
representing the presence or absence of a building, respectively. We refer to
this Markov network as degenenrate network, characterized by the absence of
functions that capture the interactions between nodes.

To incorporate the change outputs for a specific pixel, we first include
additional edges for the edge settings cyclic and dense. Specifically, for the
cyclic case (Figures 5(b)), edge $\mathcal{E}^{t1\rightleftharpoons T}$,
connecting the first node $\mathcal{N}^{t1}$ and the last node
$\mathcal{N}^{T}$, is added. On the other hand, for the dense case (Figure
5(c)), all possible non-adjacent edges are added to the graph. Then, we define
factors over the edges in the Markov network to add pairwise interactions of
connected nodes. Specifically, we define pairwise potentials $\phi_{n}$ for
each edge $\mathcal{E}^{n}$, connecting two nodes $\mathcal{N}^{t}$ and
$\mathcal{N}^{k}$. Since all variables in the network are binary, each factor
over an edge has four parameters. The change outputs are encoded as edge
potentials for the four combinations of states, as follows:

$\begin{split}\phi_{n}(\mathcal{N}^{t}=0,\mathcal{N}^{k}=1)&=\mathbf{O}_{\rm
ch}^{(i,j),n}\\\
\phi_{n}(\mathcal{N}^{t}=1,\mathcal{N}^{k}=0)&=\mathbf{O}_{\rm
ch}^{(i,j),n}\\\
\phi_{n}(\mathcal{N}^{t}=0,\mathcal{N}^{k}=0)&=1-\mathbf{O}_{\rm
ch}^{(i,j),n}\\\
\phi_{n}(\mathcal{N}^{t}=1,\mathcal{N}^{k}=1)&=1-\mathbf{O}_{\rm
ch}^{(i,j),n},\end{split}$ (11)

where edge $n$ is connecting timestamps $t$ and $k$.

The value associated with each particular assignment of states denotes the
affinity between the two states. Consequently, the higher the value assigned
to the edge potential for a particular combination of states, the more
compatible these two states are.

To define a global model from the local interactions defined in the
parameterization of the Markov network, we take the product of the local
factors and then normalize it. Once the distribution is defined, we perform a
MAP query to find the most likely state assignment for each node in the graph.
To that end, inference is performed using the belief propagation algorithm
(see Algorithm 10.4 in [46]). Due to the absence of loops in the graph, belief
propagation provides an exact solution. Finally, it should be noted that the
MTI module does not contain any trainable parameters and is only deployed
during inference.

((a)) Adjacent

((b)) Cyclic

((c)) Dense

Figure 5: Markov networks with different edge settings (exemplified for a time
series with length $T=5$). Nodes and Edges are denoted by $\mathcal{N}$ and
$\mathcal{E}$, respectively.

##  IV Experimental setting

###  IV-A Datasets

SpaceNet 7: The SpaceNet 7 dataset features time series of satellite images
acquired by the PlanetScope constellation between 2017 and 2020 for 60 sites
spread across the globe [41]. Each time series consists of about 24 monthly
mosaics with a spatial resolution of 4 m (approximately 1024 x 1024 pixels).
Furthermore, the SpaceNet 7 dataset provides manually annotated building
footprints, whereas annotations are missing for image parts affected by
clouds. While the task of the original SpaceNet 7 challenge was to track these
building footprints (i.e., vector data), the SpaceNet 7 dataset has been
leveraged for a diverse set of tasks such as urban mapping [47], building
counting [48], and urban change forecasting [49]. To split the SpaceNet 7
dataset into training, validation, and test areas, we apply the within-scene
splits recommended in [48]. Specifically, approximately 80 % of a scene (top
part) is used as a source of training patches, and the remaining 20 % (bottom
part) is divided evenly into validation (bottom-left part) and test (bottom-
right part) areas. Within-scene splits minimize the occurrence of out-of-
distribution data during testing while simultaneously avoiding data leakage
between the training and test set by utilizing spatially disjoint areas for
the different dataset splits. During training, samples from the training areas
are generated by randomly selecting $T$ timestamps from the time series of a
site. The rasterized building labels (see [34]) for these timestamps, were
used to compute the change label. We draw 100 samples from each site during an
epoch to reach an adequate number of steps before model evaluation. For model
evaluation (validation and testing), the first and the last cloud-free images
of a time series, in addition to evenly spaced intermediate images, were
selected.

WUSU: The Wuhan Urban Semantic Understanding (WUSU) dataset features tri-
temporal high-resolution Gaofen-2 images covering Wuhan’s Hongshan District
and Jiang’an District [39]. The preprocessing workflow of the satellite images
includes orthographic correction and multi-temporal image registration.
Furthermore, the four multi-spectral bands acquired at 4 m spatial resolution
are pansharpened to a spatial resolution of 1 m, resulting in images of size
6358 x 6382 and 7025 x 5500 pixels for Hongshan District and Jiang’an
District, respectively. In addition to the Gaofen-2 images, the WUSU dataset
provides corresponding land-use/land-cover (LULC) labels, including manually
refined building annotations (Class 2 Low building and Class 3 High building).
Since the proposed method requires binary building labels, Class 2 and Class 3
were remapped to the foreground class, whereas all other classes were remapped
to the background class. We follow the within-scene split recommended by the
authors, using the top halves of the six images for the test set and the
bottom halves for the training set that was further divided into training (90
%) and validation (10 %) tiles.

###  IV-B Baseline and benchmark methods

We selected a comprehensive set of baseline and benchmark methods for
quantitative and qualitative comparisons with the proposed methods. These
selected methods are listed under the categories bi-temporal change detection,
multi-temporal change detection, and segmentation methods below.

Bi-temporal urban change detection methods

  1. 1.

Siam-Diff [13] employs a Siamese encoder to extract feature maps from bi-
temporal images. A decoder produces the change prediction from the subtracted
feature maps. The encoders and decoder follow the U-Net architecture.

  2. 2.

SNUNet [14] replaces the architecture in the Siam-Diff network with a Nested
U-Net (UNet++ [42]). In addition, a channel attention module is incorporated
into the architecture.

  3. 3.

DTCDSCN [33] combines a typical Siamese ConvNet for bi-temporal change
detection with a dual attention module and two additional decoders with shared
weights for building segmentation.

  4. 4.

BIT [18] employs a bi-temporal image transformer module that operates in a
compact token space to refine features extracted by a Siamese ConvNet.

  5. 5.

AMTNet [16] also extracts features using a Siamese ConvNet, and combines
attention mechanisms and multi-scale processing techniques to model contextual
information in bi-temporal images.

  6. 6.

ScratchFormer [21] introduces shuffled sparse attention layers in a Siamese
ConvNet encoder to effectively capture semantic changes when training from
scratch.

Multi-temporal urban change detection methods

  1. 7.

L-UNet [26] combines a U-Net for feature extraction in image time series with
LSTM modules [28] that obtain rich temporal information from the extracted
feature maps.

  2. 8.

Multi-Task L-UNet [27] adds a semantic decoder to L-UNet to segment buildings
in the first and last images of a time series.

Building segmentation methods

  1. 9.

U-Net [10] is an encoder-decoder ConvNet widely used for semantic segmentation
tasks, including satellite image segmentation.

  2. 10.

HRNet [50], or High-Resolution Net, is a ConvNet that maintains high-
resolution representations through the whole network. For semantic
segmentation, the high-resolution feature map is combined with the upsampled
lower-resolution maps [11].

The Siamese ConvNets Siam-Diff, SNUNet, and DTCDSCN are commonly used as
change detection baselines, whereas BIT, AMTNet, and ScratchFormer represent
recent methods combining Siamese ConvNets with transformers. On the other
hand, L-UNet and Multi-Task L-UNet are benchmark methods for multi-temporal
change detection. Finally, HRNet was employed by the winning approach of the
SpaceNet 7 challenge for building segmentation, whereas U-Net was used as the
baseline method [40]. Both methods also served as a LULC segmentation baseline
for the WUSU dataset [39]. Finally, it should be noted that Multi-Task L-UNet
and DTCDSCN are multi-task methods that perform change detection and building
segmentation.

###  IV-C Model evaluation

Two accuracy metrics were used for the quantitative assessment of predicted
changes: F1 score and intersection over union (IoU). The combination of F1
score and IoU is commonly used for performance assessments in CD studies,
e.g., [35, 36]. Formulas for the metrics are given in Eq. 12 and 13, where TP,
FP, and FN represent the number of true positive, false positive, and false
negative pixels, respectively.

$\text{F1 score}=\frac{TP}{TP+\frac{1}{2}(FP+FN)}$ (12)

$\text{IoU}=\frac{TP}{TP+FP+FN},$ (13)

Using these two accuracy metrics, we assessed model performance across three
tasks to accommodate the large variety of baseline and benchmark methods.
These tasks are described in detail in the following:

  * •

Bi-temporal change detection measures the accuracy of the predicted changes
between the first and last image of a time series.

  * •

Continuous change detection measures the average accuracy of the predicted
changes between consecutive image pairs in a time series.

  * •

Segmentation measures the accuracy of the building predictions corresponding
to the first and last image of a time series.

The first task focuses on urban change detection from image pairs acquired
multiple years apart. This task is considered by most urban change detection
methods. For bi-temporal change detection methods, changes were directly
predicted based on the first-last image pair, ignoring intermediate images in
a time series. On the other hand, the second task focuses on assessing change
predictions between consecutive image pairs in a time series. Consequently,
the continuous urban change detection task focuses on image pairs with periods
between acquisition dates that are considerably shorter (i.e., annual and sub-
annual). The last task assesses the auxiliary segmentation task of multi-task
methods and segmentation models. It should be noted that the change detection
performance of segmentation models is not assessed because post-classification
comparison suffers from the accumulation of classification errors [51].

###  IV-D Implementation details

We implement the proposed method using the deep learning framework PyTorch
[52]. To efficiently reshape feature maps, we use the rearrange
operation111See https://einops.rocks/api/rearrange/. For the MTI module, we
used the pgmpy222See https://pgmpy.org/index.html package to implement the
Markov network and perform belief propagation. Models were trained for a
maximum duration of 100 epochs on NVIDIA GeForce RTX 3080 graphics cards,
using early stopping with patience 10 to prevent models from overfitting to
the training set. AdamW was used as optimizer [53] with a linear learning rate
scheduler. The remainder of this section describes the training setup in
detail.

Augmentations: To enhance the training dataset, we applied four data
augmentation operations, namely rotations ($k\cdot 90^{\circ}$, where $k$ is
randomly selected from $\\{0,1,2,3\\}$), flips (horizontal and vertical with a
probability of 50 %), Gaussian blur, and random color jittering. The
parameters that determine how much to jitter the brightness, contrast,
saturation, and hue of an image were set to 0.3 [19]. For validation and
testing, on the other hand, no data augmentation was applied.

Oversampling: To account for the fact that the occurrence of change is usually
considerably less frequent than no change [54], change areas were oversampled
during network training. For a given site, twenty patches of size $64\times
64$ pixels were randomly cropped from the change label, before assigning each
patch a probability according to its change pixel percentage, including a base
probability for patches with no change pixels. A single patch was chosen based
on those probabilities. For transformer-based methods (BIT, AMTNet, and
Scratchformer), the patch size was increased to $128\times 128$ pixels to
include more long-range spatial context.

Hyper-parameter tuning: For each model, hyper-parameters were tuned
empirically on the validation set using grid search. Specifically, an
exhaustive search with three learning rates ($1\cdot 10^{-5}$, $5\cdot
10^{-5}$, $1\cdot 10^{-4}$) and two batch sizes (8, 16) was performed to
determine the optimum values of hyper-parameters. Then, five models were
trained with the best hyper-parameters but different seeds for weight
initialization and data shuffling. Consequently, reported values correspond to
the average of five runs.

Multi-Task L-UNet requires additional hyper-parameters to balance the
contribution of the segmentation and change loss terms, which we adopted from
the paper [27]. All bi-temporal urban change detection methods were trained on
the cyclic edges setting (see Figure 5(b)).

##  V Experiment results

In this section, we present the quantitative and qualitative results on the
SpaceNet 7 and WUSU datasets, and the ablation study results. It should be
noted that all accuracy values are reported on the respective test sets and
correspond to the mean values obtained from five models. These were trained
with different seeds using the best hyper-parameters determined with a grid
search (see Section IV-D). On the other hand, the median model is used for the
qualitative results.

###  V-A SpaceNet 7

The quantitative results for the SpaceNet 7 dataset are listed in Table I. The
proposed method achieved the highest F1 scores and IoU values for both urban
change detection tasks (i.e., bi-temporal and continuous). For bi-temporal
change detection, methods using multi-temporal information (L-UNet and Multi-
Task L-UNet) clearly outperform bi-temporal methods. Among the latter methods,
the highest accuracy values were achieved by the ConvNet-based methods Siam-
Diff and SNUNet, in addition to the transformer-based method ScratchFromer.
For building segmentation, the best performance was also achieved by the
proposed method, followed by the segmentation methods U-Net and HRNet and the
multi-task methods Multi-Task L-UNet and DTCDSCN.

TABLE I: Quantitative results on the SpaceNet 7 test areas. The best and second-best performance is highlighted in red and blue, respectively. ”-” denotes that the accuracy metric does not apply to a specific method since the corresponding variable is not predicted. Method | Change detection | Segmentation
---|---|---
Bi-temporal | Continuous |  |
F1 | IoU | F1 | IoU | F1 | IoU
Siam-Diff | 0.453 | 0.293 | 0.273 | 0.158 | - | -
SNUNet | 0.454 | 0.294 | 0.300 | 0.177 | - | -
DTCDSCN | 0.413 | 0.260 | 0.250 | 0.143 | 0.488 | 0.323
BIT | 0.386 | 0.239 | 0.275 | 0.160 | - | -
AMTNet | 0.424 | 0.269 | 0.282 | 0.164 | - | -
ScratchF. | 0.469 | 0.306 | 0.326 | 0.195 | - | -
L-UNet | 0.519 | 0.350 | - | - | - | -
MT L-UNet | 0.515 | 0.347 | - | - | 0.512 | 0.344
U-Net | - | - | - | - | 0.527 | 0.359
HRNet | - | - | - | - | 0.527 | 0.358
Proposed | 0.551 | 0.381 | 0.414 | 0.261 | 0.596 | 0.424

Qualitative results for the SpaceNet 7 dataset are shown for urban change
detection and building segmentation in Figures 6 and 7, respectively. The
change detection results are shown for test sites located in the United
States, China, the Philippines, and Australia (top to bottom). Across these
sites, the proposed method detects urban changes more accurately than other
methods. For example, our method produces the fewest FP pixels (green) for the
site located in the Philippines (third row). In addition, it reduces the
number of FN pixels (magenta) compared to the other methods for the sites
located in China (second row) and Australia (fourth row). The proposed method
also accurately maps buildings, as demonstrated in Figure 7 for test sites
located (from top to bottom) in the United States (first and second row), the
Netherlands (third row), and Australia (fourth row).

Figure 6: Qualitative comparison of the predicted urban changes between the
first (first column) and last (second column) image of a time series for a
selection of the SpaceNet 7 test areas. Columns three to seven show the
benchmarks and the last column shows the proposed method (White: TPs, Black:
TNs, Green: FPs, and Magenta: FNs). Figure 7: Qualitative comparison of the
built-up area predictions corresponding to the last image (first column) of a
time series for a selection of the SpaceNet 7 test areas. Columns two to four
show the benchmarks and the last column shows the proposed method (White: TPs,
Black: TNs, Green: FPs, and Magenta: FNs).

###  V-B WUSU

The quantitative results for the WUSU datasets are listed in Table II. In
general, change detection performance on the WUSU dataset is lower than on the
SpaceNet 7 dataset. For bi-temporal change detection, only ScratchFormer,
L-UNet, and the multi-task methods exceed an F1 score of 0.275 and an IoU
value of 0.160. In comparison, accuracy values for the continuous change
detection task are slightly higher, except for the proposed method. Overall,
our method outperforms all other methods on both change detection tasks. For
building segmentation, F1 scores of 0.500 and IoU values of 0.300 are exceeded
by all methods with a segmentation output except Multi-Task L-UNet. U-Net and
HRNet outperform the multi-task methods; however, the best results were
achieved by the proposed method with an F1 score of 0.663 and an IoU value of
0.496.

TABLE II: Quantitative results on the WUSU test areas. The best and second-best performance is highlighted in red and blue, respectively. ”-” denotes that the accuracy metric does not apply to a specific method since the corresponding variable is not predicted. Method | Change detection | Segmentation
---|---|---
Bi-temporal | Continuous |  |
F1 | IoU | F1 | IoU | F1 | IoU
Siam-Diff | 0.175 | 0.096 | 0.236 | 0.134 | - | -
SNUNet | 0.188 | 0.104 | 0.211 | 0.118 | - | -
DTCDSCN | 0.278 | 0.162 | 0.318 | 0.189 | 0.539 | 0.369
BIT | 0.213 | 0.120 | 0.314 | 0.187 | - | -
AMTNet | 0.187 | 0.104 | 0.264 | 0.152 | - | -
ScratchF. | 0.324 | 0.193 | 0.353 | 0.214 | - | -
L-UNet | 0.279 | 0.162 | - | - | - | -
MT L-UNet | 0.276 | 0.161 | - | - | 0.479 | 0.315
U-Net | - | - | - | - | 0.576 | 0.406
HRNet | - | - | - | - | 0.611 | 0.440
Proposed | 0.440 | 0.282 | 0.389 | 0.242 | 0.663 | 0.496

Qualitative results for the WUSU dataset are shown for urban change detection
and building segmentation in Figures 8 and 9, respectively. All change
detection methods considered in this visual comparison suffer from a large
number of FN pixels (purple). In contrast, the proposed method is considerably
more effective in detecting urban changes. In addition, the segmentation
results for the last image in the time series (see Figure 9) demonstrate that
the selected methods accurately map built-up area, even under challenging
conditions with considerable off-nadir angles and shadows. Most incorrect
predictions (FNs and FPs) are present at the edges of buildings, but a few
buildings were also completely omitted. In general, the proposed method
produces fewer FPs and FNs than the selected methods in this comparison.

Figure 8: Qualitative comparison of the predicted urban changes between the
first (first column) and last (second column) image of a time series for a
selection of the WUSU test areas. Columns three to seven show the benchmarks
and the last column shows the proposed method (White: TPs, Black: TNs, Green:
FPs, and Magenta: FNs). Figure 9: Qualitative comparison of the built-up area
predictions corresponding to the last image (first column) of a time series
for a selection of the WUSU test areas. Columns two to four show the
benchmarks and the last column shows the proposed method (White: TPs, Black:
TNs, Green: FPs, and Magenta: FNs).

###  V-C Ablation study

Loss function and TFR module: To measure the contribution of applying the loss
function to additional change edges in image time series, we perform an
ablation experiment that tests different change loss settings on the SpaceNet
7 and WUSU datasets. Furthermore, we investigate the contribution of the TFR
module by removing it from the model. The results listed in Table III show
that adding change edges in the loss function generally improves performances
of the segmentation and change detection tasks. Regarding the TFR module, the
results in Table III demonstrate a clear performance gain for both tasks from
introducing transformer-based feature refinement. Specifically, adding the TFR
module improves segmentation and change detection performance on both datasets
considerably.

TABLE III: Ablation results for the loss function and the TFR module on the SpaceNet 7 (SN7) and WUSU datasets. The change loss settings cyclic and dense are equivalent for the WUSU dataset since its time series consists of three images. Dataset | Change loss | TFR module | Change detection | Segmentation
---|---|---|---|---
Bi-temporal | Continuous |  |
F1 | IoU | F1 | IoU | F1 | IoU
SN7 | adjacent | ✓ | 0.520 | 0.352 | 0.384 | 0.238 | 0.584 | 0.413
cyclic | ✓ | 0.532 | 0.363 | 0.377 | 0.233 | 0.581 | 0.409
dense | ✓ | 0.537 | 0.367 | 0.397 | 0.248 | 0.593 | 0.422
dense | ✗ | 0.490 | 0.324 | 0.314 | 0.186 | 0.549 | 0.379
WUSU | adjacent | ✓ | 0.356 | 0.217 | 0.373 | 0.229 | 0.650 | 0.482
cyclic | ✓ | 0.420 | 0.266 | 0.391 | 0.243 | 0.660 | 0.493
cyclic | ✗ | 0.342 | 0.208 | 0.364 | 0.223 | 0.583 | 0.412

MTI module: We perform another ablation experiment investigating performance
gains from leveraging additional edges in the MTI module. A degenerate Markov
network that uses no change information (i.e., only the segmentation
information represented as nodes) is added as a baseline. Table V-C lists the
ablation results for the SpaceNet 7 and WUSU datasets. Introducing change
information in the MTI module results in considerable change detection
performance gains compared to the degenerate setting while also providing
minor segmentation performance gains. For example, using the adjacent setting
improves change detection performance (first-last) by 44.8 % and 61.6 % in
terms of F1 score and IoU, respectively, compared to the degenerate setting on
the SpaceNet 7 dataset. On the WUSU dataset, the corresponding performance
improvements are 57.5 % and 73.6 % for the F1 score and IoU, respectively
(cyclic loss scenario). Table V-C also shows that introducing change
information beyond adjacent edges results in further change detection
performance gains, albeit to a much lesser extent. For example, compared to
only using adjacent change information, dense information improved the F1
score and IoU values by 2.0 % and 3.0 %, respectively. In the context of the
WUSU dataset, additional change information (i.e., adjacent vs. cyclic)
improved the F1 score and IoU values by 1.6 % and 2.2 %, respectively.
Therefore, the ablation study demonstrates that the proposed MTI module
effectively integrates the outputs of the segmentation and change detection
tasks at inference time.

TABLE IV: Ablation results for the MTI module on the SpaceNet 7 and WUSU
datasets. The change loss settings cyclic and dense are equivalent for the
WUSU dataset since its time series consists of three images.

Dataset | Change loss | Edges | Change detection | Segmentation
---|---|---|---|---
Bi-temporal | Continuous |  |
F1 | IoU | F1 | IoU | F1 | IoU
SpaceNet 7 | adjacent | degenerate | 0.369 | 0.226 | 0.170 | 0.093 | 0.584 | 0.413
adjacent | 0.533 | 0.363 | 0.398 | 0.248 | 0.586 | 0.415
2-9 | cyclic | degenerate | 0.362 | 0.221 | 0.164 | 0.089 | 0.581 | 0.409
| adjacent | 0.526 | 0.357 | 0.391 | 0.243 | 0.583 | 0.411
| cyclic | 0.535 | 0.366 | 0.394 | 0.246 | 0.583 | 0.412
2-9 | dense | degenerate | 0.373 | 0.229 | 0.173 | 0.095 | 0.593 | 0.422
| adjacent | 0.540 | 0.370 | 0.410 | 0.258 | 0.595 | 0.424
| cyclic | 0.547 | 0.377 | 0.412 | 0.260 | 0.596 | 0.424
| dense | 0.551 | 0.381 | 0.414 | 0.261 | 0.596 | 0.424
WUSU | adjacent | degenerate | 0.244 | 0.139 | 0.186 | 0.103 | 0.650 | 0.482
adjacent | 0.409 | 0.257 | 0.368 | 0.226 | 0.656 | 0.488
2-9 | cyclic | degenerate | 0.275 | 0.159 | 0.203 | 0.113 | 0.660 | 0.493
| adjacent | 0.433 | 0.276 | 0.389 | 0.242 | 0.663 | 0.497
| cyclic | 0.440 | 0.282 | 0.389 | 0.242 | 0.663 | 0.496

Time series length: To investigate the effect of satellite image time series
length on performance, we tested the proposed network with different settings
for $T$. For the change loss and MTI module, the maximum number of edges
(i.e., dense edge setting) for a given $T$ was used (see Equation 6). However,
it should be noted that for lengths $T=3$ and $T=2$, the edge settings dense
are equivalent to cyclic and adjacent, respectively. The results of this
experiment (Figure 10) show that adding intermediate images has a positive
impact on bi-temporal change detection performance. Specifically, while the
proposed method performed worse for $T=2$ compared to SOTA methods (Table I),
our model outperforms these methods, including multi-temporal change detection
methods, at longer time series lengths. Therefore, it effectively uses multi-
temporal information to improve change detection performance.

Figure 10: The effect of time series length on bi-temporal change detection
performance for the SpaceNet 7 dataset. The maximum number of edges was used
for each time series length in the change loss and MTI module. Values
represent the mean $\pm$ 1 standard deviation of 5 runs.

##  VI Discussion

Our results highlight the challenging nature of continuously detecting urban
changes from high-resolution satellite image time series. Most methods are
developed for bi-temporal urban change detection from VHR imagery, and SOTA
methods have achieved remarkable results with F1 scores exceeding 0.9 on
datasets such as LEVIR-CD [17] and WHU-CD [44]. In contrast, studies focusing
on urban change detection from high-resolution imagery (e.g., PlanetScope or
Sentinel-2) often fall short of F1 scores of 0.6, even when employing multi-
task learning and multi-temporal data [27]. The suite of methods that achieved
these impressive performances typically employs transformers to capture long-
range contextual information in bi-temporal VHR imagery [17, 18, 19, 16]. Our
experiments indicate, however, that improving long-range context is less
relevant in high-resolution imagery (Tables I and II). In addition,
transformer-based methods rely heavily on pretraining on large datasets, which
ScratchFormer successfully mitigates [21]. Furthermore, a recent study
suggests that a simple segmentation ConvNet such as U-Net is still a top
performer for change detection tasks [55]. We attribute those findings to the
capability of modern ConvNets to capture sufficient spatial context in
satellite imagery, especially when ground sampling distances are above 1 m. In
comparison to the SOTA methods for bi-temporal change detection, our method
leverages the self-attention mechanism to model multi-temporal information in
satellite image time series. We showed that the proposed TFR module
contributes to the network’s representation learning capability, resulting in
improved change detection and semantic segmentation performance (Table III).
Therefore, we deem self-attention an effective mechanism for feature
refinement in continuous urban change detection networks.

Our work also highlights the need for effective integration approaches in
multi-task learning schemes. Specifically, although multi-task learning is
commonly applied for change detection [33, 27], existing multi-task studies do
not address the integration of the semantic segmentation and change detection
outputs. To fill this research gap, we proposed the MTI module that represents
segmentation and change predictions using Markov networks to find the most
likely built-up area state for each timestamp in a pixel time series. Our
results demonstrate that the proposed integration approach improves both
tasks, namely the change detection and building segmentation task (Table V-C).
Furthermore, we demonstrate that the proposed approach benefits from
integrating dense change information, obtained from predicting changes between
all possible combinations of satellite image pairs in a time series, compared
to using only adjacent change information (Table V-C).

Despite the improvements our method achieves over existing methods, we also
identified several limitations related to our work. First of all, our
integration approach relies on meaningful potentials extracted from the multi-
task network outputs. However, the outputs of deep networks may not be well-
calibrated [56]. Furthermore, we assumed that our networks do not encounter
out-of-distribution data during deployment due to the within-scene splits. In
practice, however, urban mapping and change detection methods may encounter
domain shifts when deployed to unseen geographic areas [47, 36]. Therefore,
future work will test the susceptibility of our multi-task integration
approach to out-of-distribution data. For example, the effectiveness of the
MTI module could be improved by explicitly calibrating the segmentation and
change outputs of the model, using calibration techniques such as temperature
scaling [56].

Another limitation of the proposed method is that it requires continuous
building labels for training. Most popular urban change detection datasets are
bi-temporal and feature VHR imagery. On the other hand, few urban change
detection datasets featuring satellite image time series and corresponding
building labels for each image are available. Therefore, weakly supervised
methods, using partial annotation or less accurate labeling, should be
investigated for continuous urban change detection from satellite image time
series (e.g., [29]).

##  VII Conclusion

In this study, we propose a continuous urban change detection method using
optical satellite image time series. The proposed method leverages a
transformer-based module to temporally refine feature representations
extracted from image time series using a ConvNet. Furthermore, a new multi-
task integration approach based on pairwise Markov networks is proposed to
combine building segmentation and dense urban change information. We
demonstrated the effectiveness of our method on the SpaceNet 7 and WUSU
datasets, where it outperformed existing bi-temporal and multi-temporal change
detection methods and segmentation methods. The ablation study attributes
these accuracy gains to the proposed TFR and MTI modules. Therefore, this
research demonstrates the potential of high-resolution satellite image time
series for the continuous detection of urban changes.

## Acknowledgement

This work was supported by the EU Horizon 2020 project HARMONIA under
agreement No. 101003517, Digital Futures under the grant for the EO-
AI4GlobalChange project, and the EO-AI4Urban project within the ESA-MOST
Dragon 5 Programme.

## References

  * [1] X. Liu, Y. Huang, X. Xu, X. Li, X. Li, P. Ciais, P. Lin, K. Gong, A. D. Ziegler, A. Chen _et al._ , “High-spatiotemporal-resolution mapping of global urban change from 1985 to 2015,” _Nature Sustainability_ , pp. 1–7, 2020.
  * [2] J. Gao and B. C. O’Neill, “Mapping global urban land for the 21st century with data-driven simulations and shared socioeconomic pathways,” _Nature communications_ , vol. 11, no. 1, p. 2302, 2020.
  * [3] M. P. Johnson, “Environmental impacts of urban sprawl: a survey of the literature and proposed research agenda,” _Environment and planning A_ , vol. 33, no. 4, pp. 717–735, 2001.
  * [4] S. A. Sarkodie, P. A. Owusu, and T. Leirvik, “Global effect of urban sprawl, industrialization, trade and economic development on carbon dioxide emissions,” _Environmental Research Letters_ , vol. 15, no. 3, p. 034049, 2020.
  * [5] S. Arshad, S. R. Ahmad, S. Abbas, A. Asharf, N. A. Siddiqui, and Z. ul Islam, “Quantifying the contribution of diminishing green spaces and urban sprawl to urban heat island effect in a rapidly urbanizing metropolitan city of pakistan,” _Land Use Policy_ , vol. 113, p. 105874, 2022.
  * [6] Y. Ban and O. Yousif, “Change detection techniques: A review,” _Multitemporal Remote Sensing_ , pp. 19–43, 2016.
  * [7] D. Lu, P. Mausel, E. Brondizio, and E. Moran, “Change detection techniques,” _International journal of remote sensing_ , vol. 25, no. 12, pp. 2365–2401, 2004.
  * [8] X. X. Zhu, D. Tuia, L. Mou, G.-S. Xia, L. Zhang, F. Xu, and F. Fraundorfer, “Deep learning in remote sensing: A comprehensive review and list of resources,” _IEEE Geoscience and Remote Sensing Magazine_ , vol. 5, no. 4, pp. 8–36, 2017.
  * [9] H. Jiang, M. Peng, Y. Zhong, H. Xie, Z. Hao, J. Lin, X. Ma, and X. Hu, “A survey on deep learning-based change detection from high-resolution remote sensing images,” _Remote Sensing_ , vol. 14, no. 7, p. 1552, 2022.
  * [10] O. Ronneberger, P. Fischer, and T. Brox, “U-net: Convolutional networks for biomedical image segmentation,” in _International Conference on Medical image computing and computer-assisted intervention_. Springer, 2015, pp. 234–241.
  * [11] K. Sun, Y. Zhao, B. Jiang, T. Cheng, B. Xiao, D. Liu, Y. Mu, X. Wang, W. Liu, and J. Wang, “High-resolution representations for labeling pixels and regions,” _arXiv preprint arXiv:1904.04514_ , 2019.
  * [12] R. C. Daudt, B. Le Saux, A. Boulch, and Y. Gousseau, “Urban change detection for multispectral earth observation using convolutional neural networks,” in _IGARSS 2018-2018 IEEE International Geoscience and Remote Sensing Symposium_. IEEE, 2018, pp. 2115–2118.
  * [13] R. C. Daudt, B. Le Saux, and A. Boulch, “Fully convolutional siamese networks for change detection,” in _2018 25th IEEE International Conference on Image Processing (ICIP)_. IEEE, 2018, pp. 4063–4067.
  * [14] S. Fang, K. Li, J. Shao, and Z. Li, “Snunet-cd: A densely connected siamese network for change detection of vhr images,” _IEEE Geoscience and Remote Sensing Letters_ , vol. 19, pp. 1–5, 2021.
  * [15] K. Basavaraju, N. Sravya, S. Lal, J. Nalini, C. S. Reddy, and F. Dell’Acqua, “Ucdnet: A deep learning model for urban change detection from bi-temporal multispectral Sentinel-2 satellite images,” _IEEE Transactions on Geoscience and Remote Sensing_ , 2022.
  * [16] W. Liu, Y. Lin, W. Liu, Y. Yu, and J. Li, “An attention-based multiscale transformer network for remote sensing image change detection,” _ISPRS Journal of Photogrammetry and Remote Sensing_ , vol. 202, pp. 599–609, 2023.
  * [17] H. Chen and Z. Shi, “A spatial-temporal attention-based method and a new dataset for remote sensing image change detection,” _Remote Sensing_ , vol. 12, no. 10, p. 1662, 2020.
  * [18] H. Chen, Z. Qi, and Z. Shi, “Remote sensing image change detection with transformers,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 60, pp. 1–14, 2021.
  * [19] W. G. C. Bandara and V. M. Patel, “A transformer-based siamese network for change detection,” in _IGARSS 2022-2022 IEEE International Geoscience and Remote Sensing Symposium_. IEEE, 2022, pp. 207–210.
  * [20] V. Marsocci, V. Coletta, R. Ravanelli, S. Scardapane, and M. Crespi, “Inferring 3d change detection from bitemporal optical images,” _ISPRS Journal of Photogrammetry and Remote Sensing_ , vol. 196, pp. 325–339, 2023.
  * [21] M. Noman, M. Fiaz, H. Cholakkal, S. Narayan, R. M. Anwer, S. Khan, and F. S. Khan, “Remote sensing change detection with transformers trained from scratch,” _IEEE Transactions on Geoscience and Remote Sensing_ , 2024.
  * [22] D. P. Roy, H. Huang, R. Houborg, and V. S. Martins, “A global analysis of the temporal availability of planetscope high spatial resolution multi-spectral imagery,” _Remote Sensing of Environment_ , vol. 264, p. 112586, 2021.
  * [23] Z. Zhu, S. Qiu, and S. Ye, “Remote sensing of land change: A multifaceted perspective,” _Remote Sensing of Environment_ , vol. 282, p. 113266, 2022.
  * [24] A. Shafique, G. Cao, Z. Khan, M. Asad, and M. Aslam, “Deep learning-based change detection in remote sensing images: A review,” _Remote Sensing_ , vol. 14, no. 4, p. 871, 2022.
  * [25] X. Liu, G. Hu, Y. Chen, X. Li, X. Xu, S. Li, F. Pei, and S. Wang, “High-resolution multi-temporal mapping of global urban land using Landsat images based on the Google Earth Engine platform,” _Remote sensing of environment_ , vol. 209, pp. 227–239, 2018.
  * [26] M. Papadomanolaki, S. Verma, M. Vakalopoulou, S. Gupta, and K. Karantzalos, “Detecting urban changes with recurrent neural networks from multitemporal Sentinel-2 data,” in _IGARSS 2019-2019 IEEE International Geoscience and Remote Sensing Symposium_. IEEE, 2019, pp. 214–217.
  * [27] M. Papadomanolaki, M. Vakalopoulou, and K. Karantzalos, “A deep multitask learning framework coupling semantic segmentation and fully convolutional LSTM networks for urban change detection,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 59, no. 9, pp. 7651–7668, 2021.
  * [28] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” _Neural computation_ , vol. 9, no. 8, pp. 1735–1780, 1997.
  * [29] K. Meshkini, F. Bovolo, and L. Bruzzone, “Multi-annual change detection using a weakly supervised 3d cnn in hr sits,” _IEEE Geoscience and Remote Sensing Letters_ , pp. 1–1, 2024.
  * [30] S. Saha, F. Bovolo, and L. Bruzzone, “Change detection in image time-series using unsupervised LSTM,” _IEEE Geoscience and Remote Sensing Letters_ , 2020.
  * [31] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin, “Attention is all you need,” _Advances in neural information processing systems_ , vol. 30, 2017.
  * [32] R. Caruana, “Multitask learning,” _Machine learning_ , vol. 28, pp. 41–75, 1997.
  * [33] Y. Liu, C. Pang, Z. Zhan, X. Zhang, and X. Yang, “Building change detection for remote sensing images using a dual-task constrained deep siamese convolutional network model,” _IEEE Geoscience and Remote Sensing Letters_ , vol. 18, no. 5, pp. 811–815, 2020.
  * [34] S. Hafner, Y. Ban, and A. Nascetti, “Urban change detection using a dual-task siamese network and semi-supervised learning,” in _IGARSS 2022-2022 IEEE International Geoscience and Remote Sensing Symposium_. IEEE, 2022, pp. 1071–1074.
  * [35] Q. Shu, J. Pan, Z. Zhang, and M. Wang, “Mtcnet: Multitask consistency network with single temporal supervision for semi-supervised building change detection,” _International Journal of Applied Earth Observation and Geoinformation_ , vol. 115, p. 103110, 2022.
  * [36] S. Hafner, Y. Ban, and A. Nascetti, “Semi-supervised urban change detection using multi-modal sentinel-1 sar and sentinel-2 msi data,” _Remote Sensing_ , vol. 15, no. 21, p. 5135, 2023.
  * [37] R. C. Daudt, B. Le Saux, A. Boulch, and Y. Gousseau, “Multitask learning for large-scale semantic change detection,” _Computer Vision and Image Understanding_ , vol. 187, p. 102783, 2019.
  * [38] Z. Zheng, Y. Zhong, S. Tian, A. Ma, and L. Zhang, “Changemask: Deep multi-task encoder-transformer-decoder architecture for semantic change detection,” _ISPRS Journal of Photogrammetry and Remote Sensing_ , vol. 183, pp. 228–239, 2022.
  * [39] S. Shi, Y. Zhong, Y. Liu, J. Wang, Y. Wan, J. Zhao, P. Lv, L. Zhang, and D. Li, “Multi-temporal urban semantic understanding based on gf-2 remote sensing imagery: from tri-temporal datasets to multi-task mapping,” _International Journal of Digital Earth_ , vol. 16, pp. 3321–3347, 2023.
  * [40] A. Van Etten and D. Hogan, “The spacenet multi-temporal urban development challenge,” _arXiv preprint arXiv:2102.11958_ , 2021.
  * [41] A. Van Etten, D. Hogan, J. Martinez-Manso, J. Shermeyer, N. Weir, and R. Lewis, “The multi-temporal urban development SpaceNet dataset,” _arXiv preprint arXiv:2102.04420_ , 2021.
  * [42] Z. Zhou, M. M. Rahman Siddiquee, N. Tajbakhsh, and J. Liang, “Unet++: A nested u-net architecture for medical image segmentation,” in _Deep Learning in Medical Image Analysis and Multimodal Learning for Clinical Decision Support: 4th International Workshop, DLMIA 2018, and 8th International Workshop, ML-CDS 2018, Held in Conjunction with MICCAI 2018, Granada, Spain, September 20, 2018, Proceedings 4_. Springer, 2018, pp. 3–11.
  * [43] L. Wang, J. Zhang, Q. Guo, and D. Chen, “Iftsdnet: An interact-feature transformer network with spatial detail enhancement module for change detection,” _IEEE Geoscience and Remote Sensing Letters_ , 2023.
  * [44] S. Ji, S. Wei, and M. Lu, “Fully convolutional networks for multisource building extraction from an open aerial and satellite imagery data set,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 57, no. 1, pp. 574–586, 2018.
  * [45] D. Duque-Arias, S. Velasco-Forero, J.-E. Deschaud, F. Goulette, A. Serna, E. Decencière, and B. Marcotegui, “On power Jaccard losses for semantic segmentation,” in _VISAPP 2021: 16th International Conference on Computer Vision Theory and Applications_ , 2021.
  * [46] D. Koller and N. Friedman, _Probabilistic graphical models: principles and techniques_. MIT press, 2009.
  * [47] S. Hafner, Y. Ban, and A. Nascetti, “Unsupervised domain adaptation for global urban extraction using sentinel-1 sar and sentinel-2 msi data,” _Remote Sensing of Environment_ , vol. 280, p. 113192, 2022.
  * [48] M. T. Razzak, G. Mateo-García, G. Lecuyer, L. Gómez-Chova, Y. Gal, and F. Kalaitzis, “Multi-spectral multi-image super-resolution of sentinel-2 with radiometric consistency losses and its effect on building delineation,” _ISPRS Journal of Photogrammetry and Remote Sensing_ , vol. 195, pp. 1–13, 2023.
  * [49] N. Metzger, M. Ö. Türkoglu, R. C. Daudt, J. D. Wegner, and K. Schindler, “Urban change forecasting from satellite images,” _PFG–Journal of Photogrammetry, Remote Sensing and Geoinformation Science_ , pp. 1–10, 2023.
  * [50] J. Wang, K. Sun, T. Cheng, B. Jiang, C. Deng, Y. Zhao, D. Liu, Y. Mu, M. Tan, X. Wang _et al._ , “Deep high-resolution representation learning for visual recognition,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 43, no. 10, pp. 3349–3364, 2020.
  * [51] S. Liu, H. Su, G. Cao, S. Wang, and Q. Guan, “Learning from data: A post classification method for annual land cover analysis in urban areas,” _ISPRS Journal of Photogrammetry and Remote Sensing_ , vol. 154, pp. 202–215, 2019.
  * [52] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga _et al._ , “Pytorch: An imperative style, high-performance deep learning library,” _arXiv preprint arXiv:1912.01703_ , 2019.
  * [53] I. Loshchilov and F. Hutter, “Fixing weight decay regularization in adam,” 2018.
  * [54] F. Bovolo and L. Bruzzone, “The time variable in data fusion: A change detection perspective,” _IEEE Geoscience and Remote Sensing Magazine_ , vol. 3, no. 3, pp. 8–26, 2015.
  * [55] I. Corley, C. Robinson, and A. Ortiz, “A change detection reality check,” _arXiv preprint arXiv:2402.06994_ , 2024.
  * [56] C. Guo, G. Pleiss, Y. Sun, and K. Q. Weinberger, “On calibration of modern neural networks,” in _International conference on machine learning_. PMLR, 2017, pp. 1321–1330.