Spatial Context Awareness for Unsupervised Change Detection in Optical Satellite Images

Detecting changes on the ground in multitemporal Earth observation data is one of the key problems in remote sensing. In this paper, we introduce Sibling Regression for Optical Change detection (SiROC), an unsupervised method for change detection in optical satellite images with medium and high resolution. SiROC is a spatial context-based method that models a pixel as a linear combination of its distant neighbors. It uses this model to analyze differences in the pixel and its spatial context-based predictions in subsequent time periods for change detection. We combine this spatial context-based change detection with ensembling over mutually exclusive neighborhoods and transitioning from pixel to object-level changes with morphological operations. SiROC achieves competitive performance for change detection with medium-resolution Sentinel-2 and high-resolution Planetscope imagery on four datasets. Besides accurate predictions without the need for training, SiROC also provides a well-calibrated uncertainty of its predictions.

imagery, the possibilities of multi-temporal analysis have increased significantly [8].Combined with the open data policy of the Copernicus program, it is, for example, possible to acquire a Sentinel-2 image with 10m resolution per pixel of any region of interest on any continent every 5 days [9] free of charge.Commercial providers of satellite imagery can even offer almost daily coverage with high-resolution imagery for large parts of the planet [10].These trends emphasize the increasing opportunities in monitoring Earth from space and the relevance of CD as a field within remote sensing.Obtaining labeled data for CD, however, is costly in terms of time and effort, especially at scale.Therefore, a large focus of attention in the design of CD algorithms is unsupervised methods that do not require ground truth [11].
The applicability of unsupervised CD methods in multispectral satellite images varies depending on the spatial resolution of input images.For very-high-resolution (VHR) imagery with a spatial resolution up to 0.5m, deep learning-based methods tend to be in general preferable because of their elaborate capacity to model spatial context [11] although most of the work in this area focuses on supervised methods [12]- [17].Since for VHR imagery an object such as a building consists of a number of pixels, modeling spatial context is essential to provide accurate unsupervised change segmentations.Saha et al. (2019) introduce Deep Change Vector Analysis (DCVA) [11], a VHR CD framework that combines ideas from image differencing with feature extraction based on pre-trained neural networks.DCVA has also been combined with selfsupervised pre-training of the feature extractor specifically for remote sensing images [18].MSDRL [19] is a scale-driven unsupervised method that uses deep feature extraction to obtain a pseudo-classification of change superpixels.Superpixels with high certainty pseudo-labels are then taken as input to train a support vector machine which eventually classifies the uncertain superpixels.Such pre-classification schemes where pseudo-labels are obtained based on another method have also been presented in conjunction with methods for unsupervised change detection in Synthetic Aperture Radar images [4], [20], [21] such as PCANet [22], [23].Gong et al. [24] introduced modeling the difference image with a generative adversarial network (GAN).While the deep learning methods above were primarily designed for high-resolution imagery, some of them can be applied to medium-resolution imagery as well.In the case of DCVA, there also exists a variant adjusted to the spatial and spectral scale of Sentinel-2 [25].
For medium-resolution CD, non-deep learning methods based and improved on Change Vector Analysis (CVA) can arXiv:2110.02068v1[cs.CV] 5 Oct 2021 still compete.CVA takes the difference of radiometric values or features derived from it over time [11] and applies a threshold to this difference image.Examples of features that have been derived from radiometric values as input for image differencing are vegetation indices [26] or tasseled cap transformation features [27].Otsu thresholding [28] has been shown to be effective for thresholding the difference image [29] although a variety of approaches exist [30]- [32].Beyond, binary change detection, the signal in the CVA difference image can also be used to uncover the type of change [33], [34].
CVA-based methods can still be insightful especially with medium-resolution because the size of objects in these images is typically assumed to be similar to the spatial resolution of a pixel.However, extensions of CVA still fall short of the deep learning-based DCVA for unsupervised CD on the OSCD benchmark.Still, the relative performance of traditional methods based on CVA improves with medium-resolution imagery compared to higher resolutions.More recent versions of CVA that can also be applied to higher resolution imagery tend to include close spatial context of pixels to some extent.Parcel CVA (PCVA) includes surrounding information of pixels by independent hierarchical segmentation at several scales [35].Robust CVA (RCVA) improves on potential co-registration errors in the CVA framework by replacing a point in the difference image with the difference to a neighboring pixel if the difference to this neighbor is smaller [29].Object CVA (OCVA) computes histograms of object sizes in an image and incorporates this information into a CVA framework [36].Image differencing methods have also been successfully combined with morphological operations which allow transitioning from the pixel to the object level [37].
Although neighboring pixels are somewhat included in the change analysis of a pixel in these extensions of CVA, the spatial extent of incorporated information is small compared to the effective window of sequential convolutional operations in neural networks.Neighborhood in this context is defined not only as the immediate neighbors to a pixel of interest, but also its larger spatial context up to a distance measure.The distant neighborhood of a pixel may help to identify changes because it is also affected by local trends in the image but unaffected by the change itself.For example, consider an explosion of a building between pre and post-image such as in the Beirut dataset used below.Analyzing the distant neighborhood allows to separate the actual change (destroyed building) from local trends such as dust and dirt remains on surrounding buildings stirred by the explosion.However, the use of distant neighborhood context has only found limited application in change detection thus far.This is particularly surprising since applications of image differencing in other domains such as astronomy emphasize the importance of the relation of a pixel to its neighborhood [38].Wang et al. (2016) present the causal pixel model (CPM) for the study of multi-temporal Kepler data which is used to spot transiting exoplanets in front of distant stars observed by the space telescope.The method is also more generally known as half-sibling regression (HSR) [39] Their task is conceptually related to a CD problem in remote sensing since it is also centered around spotting changes in multitemporal reflection intensities which should be unrelated to the acquisition conditions.In their case, these deviations hint towards a transient object in front of a distant star rather than a change on the ground but the fundamental principle is similar.They solve this task by modeling pixels as a function of their distant neighbors.With this model, it is possible to obtain a prediction for pixels in subsequent time steps based on their distant neighbors and compare the prediction to the actual value of the pixel.The size of the difference between predictions of pixels and their actual values is interpreted as the strength of the change signal.
HSR is related to the application of local binary patterns [40] for change detection in more traditional image recognition problems.Bilodeau et al. (2013) [41] design a method based on local binary similarity patterns (LBSP) to separate the image background from changes in multitemporal images.In their method, a binary similarity measure is computed between a pixel of interest and its closest neighbors within an image.If the binary similarity pattern updates notably between images, this is considered to be a change signal.A version of LBSPs has also been applied to change detection in remote sensing where multi-temporal images are split into overlapping blocks and LBSPs of these blocks are compared across time [42].Similarly, the graph structure of image patches across time has been used for homogenous and heterogenous change detection [43].The shared principle between LBSP and HSR is the approach to compare a pixel to its neighborhood and inspect how this relationship changes over time to discover potential changes.However, HSR relies mostly on distant neighborhood information rather than a small set of close neighbors and models this relationship explicitly to obtain a prediction for subsequent time periods.
One key property of HSR is the fact that it is by design comparably robust to registration errors and varying acquisition conditions for a given sensor [38].This is because variations in the acquisition conditions can also affect distant neighbors in an image whereas actual changes at the pixel level should be independent of distant context.Changing acquisition conditions and registration errors, however, are two of the primary sources of false positives in change detection [44].Since HSR deals comparably well with these issues it may work especially well for change detection in remote sensing time-series data.In this paper, we therefore apply HSR for change detection in remote sensing.When we know from astronomy that distant spatial context can improve resilience against varying acquisition conditions, this might be particularly helpful in remote sensing CD.
We modify HSR in two major ways to apply it as SiROC for change detection in remote sensing: At first, we design an ensemble version of HSR based on mutually exclusive neighborhoods.Second, we make use of morphological operations to transition from pixel-level changes to object-level changes.SiROC is tested for urban CD in medium-resolution images on the Onera Satellite Change Detection Dataset (OSCD) and high-resolution images from the Beirut Harbor Explosion of 2020 (BHED).Outside of the urban context, we test SiROC on the Barrax Agriculture dataset and the Lamar Alpine dataset.
Our main contributions are three-fold: 1) We introduce SiROC, a robust method for unsupervised change detection in optical remote sensing that combines ideas from HSR with ensembles over mutually exclusive neighborhoods and morphological operations.2) SiROC achieves competitive performance for mediumas well as high-resolution unsupervised change detection with optical images.3) SiROC also returns a built-in, well-calibrated uncertainty score with its change segmentation which makes it well suited for applications in conjunction with deep learning such as pseudo-labeling in self-or unsupervised learning.

A. Half-Sibling Regression (HSR) Image Differencing
The foundation of our method is Half-Sibling Regression which was originally developed for time-series analysis of Kepler data in astronomy [38], [39], [45], [46].Figure 1 displays the intuition of HSR and how it is applied to obtain signals of changing pixels across time in three steps.
At first, HSR models the pixel value of a star at time t as a linear combination of the pixel values of many other stars from the distant neighborhood of the pixel (Figure 1a).The result of this first step is a linear coefficient for every included neighborhood pixel.In the second step (Figure 1b), predictions for the pixel at time t + 1 are obtained with HSR based on the neighboring pixels at t + 1 and the respective linear coefficients from step 1.If steps 1 and 2 are executed for the whole image, there is a prediction for any pixel at t + 1 as well as its actual value.The predicted image at t + 1 is then subtracted from the actual image value to obtain a change signal for every pixel (Figure 1c).Intuitively, one expects a change in the pixels where the predictions based on the pixel's relation to its neighborhood in the past divert from the actual realization of the pixel.In the following, we elaborate on the formal definition of the three steps described.We restrict our description of HSR to the case of two time periods for simplicity since this is how we apply the method to remote sensing as well.
Step 1 Let I x,y,t be a pixel in a single-channel, twodimensional image time-series (I) at time t with coordinates (x, y).HSR models the pixel I x,y,t as a linear combination of a set of distant neighbors N from I. The neighborhood set has the points I i,j,t as elements such that (i, j) ∈ N x,y : Neighbors are chosen from the distant neighborhood because they might be subject to the same noise as the pixel of interest when they are selected to close to it.Wang et al. (2016) [38] require that an eligible neighbor has a distance of at least 20 pixels from the pixel of interest to be considered.This ensures that the pixel of interest and the chosen neighbors have practically no overlap in stellar illumination.The number of neighboring pixels considered is generally large and Wang et al. (2016) [38] select 4000 neighboring pixels in their original proposal of HSR to model one pixel of interest for Kepler data.Given the high temporal density of observations for each pixel (every 30 minutes) in Wang et al. (2016) [38], this is still solvable because the number of observed time periods exceeds the number of neighboring pixels used.
However, in the bitemporal case, where only one period is used for fitting, there are many potential combinations of β which solve equation (1).We derive β as the closed form solution of the least squares problem.It is a function of the pixel of interest I x,y , the respective neighbor I i,j and the quadratic sum of all neighbors I i ,j : Step 2: With the coefficients obtained in step 1, I x,y,t+1 can be predicted as Îx,y,t+1 = (i,j)∈Nx,y With the expression for β from equation (2), equation ( 3), can be rearranged to: Îx,y,t+1 = (i,j)∈Nx,y I i,j,t+1 I i,j,t where g t+1 resembles a growth rate of the sum of pixel values in the selected neighbors from t to t + 1.In essence, the assumption is that if the pixel values around I x,y,t increase by a factor g t+1 and no changes occurred at this location, I x,y,t+1 should be close to I x,y,t g t+1 .We can circumvent the explicit calculation of beta and directly obtain Îx,y,t+1 based on equation (4) which is computationally efficient.
Step 3: The difference between I x,y,t+1 and Îx,y,t+1 is taken as the change signal for pixel I x,y between t and t + 1.
After obtaining I x,y,t+1 for all (x, y) ∈ I t+1 , the residual is given as the difference of the image matrices: Note that this is slightly different from the standard application of image differencing in change vector analysis in multitemporal remote sensing.We do not directly take the difference of the image vectors at t and t + 1.Instead, we predict how the image would have looked like in t + 1 if the local neighborhood relations persisted.Then, we use this predicted image as input for image differencing with the actual image in t + 1.The extension of HSR to images with several channels is straightforward as one can directly sum the absolute values of for each channel to incorporate HSR information from all channels.Let t+1,c be the residual of channel c of a multispectral image with C channels.Then, the aggregated chance signal can be computed as:  for (channel in channels) do for (pixel in I t ) do 6: Apply HSR(n,e) to get Ît+1

B. HSR for Earth Observation Data (SiROC)
We improve and adapt the standard HSR Image Differencing model to apply it effectively for CD in remote sensing as Sibling Regression for Optical Change detection (SiROC).Algorithm ?? outlines SiROC in pseudocode.In summary, there are two major differences between SiROC and HSR Image Differencing.At first, we redesign the notion of included and excluded pixels in the neighborhood selection to create an ensemble version of HSR over mutually exclusive neighborhoods.This does not only improve performance but also allows us to obtain an uncertainty along with the prediction.The rationale of splitting neighborhoods by distance is to inspect trends at different distances separately instead of pooling the trends together.For example, two trends at different spatial scales might offset each other when pooling them although both may be a signal for change.Second, we combine HSR with morphological profiles to move from pixel-level to object-level changes since changes in remote sensing typically occur at the object level.
Ensembling: The starting point for SiROC is applying HSR to I t to obtain Ît+1 based on a set of neighboring pixels.We use all pixels that have a distance of at least e but at most n rows or columns from the pixel of interest.Graphically, this corresponds to all points in a square with width 2n and I x,y,t in its center that are not in the smaller square with width 2e around I x,y,t .Formally, a pixel I x ,y ,t is included in the set of neighbors for With Ît+1 , a channel-level difference image is obtained by taking the difference Ît+1 -I t+1 .The absolute value of the change signal is summed across the channels.We apply Otsuthresholding [28] to the resulting difference image which has been successfully used for thresholding difference images in CD before [29].Further, the evaluation of competing methods is also based on this thresholding approach.This allows for comparing relevant methods in a fixed setting.Nevertheless, Otsu-thresholding is a design choice here with a variety of alternatives that can also be used in conjunction with SiROC including the T-point method [47], the Rosin method [48] or the expectation-maximization algorithm [49], [50].
The result of the thresholding step is a binary segmentation of the difference image on the pixel level.However, in remote sensing applications, changes such as the construction of roads or buildings tend to occur at the object level.This is why object-based methods often tend to be superior for these applications [51].We rely on morphological profiles which are an established tool to bridge the gap between pixel-level change segmentations and the object level [52].
Morphological Profile: A morphological profile is the sequential application of morphological opening and closing to an image [52].We employ morphological opening and closing at one spatial filter size p.Intuitively, morphological closing helps to fill in missed pixels in detected change objects as changed.On the other hand, morphological opening removes spurious false positives when there are no other changes around them.After obtaining an object-level change segmentation for a given neighborhood size n and exclusion window e, we repeat the procedure and use new neighbors that are further away than the current set.Both, n and e increase by the same added factor s. In the next iteration, the previous neighborhood window becomes the exclusion window and a new binary segmentation based on more distant neighbors is obtained.This procedure is repeated until a maximum neighborhood size is reached.The number of models F is given as F = ((n max − n min)//s) + 1.Every model in the ensemble classifies each pixel either as change (1) or nochange (0) which can be interpreted as a voting mechanism among models.Voting mechanisms across spatial scales [53] or different bands [25] are a common aggregation mechanism in change detection.The number of votes per pixel ranges between 0 and F.
Majority Voting: The matrix of votes per pixel can be visualized as a heatmap of agreement between different sets of neighbors if a change occurred.This also directly transports a measure of uncertainty embedded in SiROC.If a pixel has no or the maximum number of votes, the agreement is high and the method is confident in its prediction.If the number of votes is split the model shows low confidence in its prediction for this point.We threshold these votes with a pre-defined voting share 0 ≤ v ≤ 1 that is required to classify a pixel as changed.v is the sensitivity of our model towards change.The choice of v contains a trade-off between objectives.With a higher v, the number of false negatives rises but false positives decline (and vice versa).Since all models are equally weighted in the voting process, the importance of a single neighbor is decreasing in its distance to the pixel of interest.This is because the number of neighbors used per model is increasing in n.The underlying assumption is that pixels closer to the point of interest carry more information about its potential change.This assumption is domain-specific to Earth observation and stands in contrast to the idea of HSR in astronomy where there is no weighting based on distance.The application of the voting threshold is the last step of SiROC to obtain the final change segmentation.The voting matrix is normalized by the number of models before the percentage threshold is applied.
To summarize, SiROC has the following hyperparameters: 1) Maximum neighborhood size: n max 2) Initial exclusion window: e start 3) Step size of ensemble: s 4) Filter Size of Morphological Operations: p 5) Voting Threshold: The initial size of the neighborhood window n start is given as e start + s.

III. EXPERIMENTS AND RESULTS
Section III-A describes the datasets used to assess the performance of SiROC and competing methods.The competing methods used as a benchmark and the evaluation criteria are described in more detail in Section III-B.The results on OSCD, BHED, the Agriculture Dataset and the Alpine Dataset are presented in depth in Sections III-C, III-D III-E, and III-F respectively.

A. Description of Datasets
Onera Change Detection Dataset (OSCD): OSCD is a benchmark for bitemporal urban CD based on multispectral Sentinel-2 images [54].It contains manual annotations of binary changes for 24 cities across the globe where 14 are used for training and 10 for testing.The labels focus on urban changes such as newly constructed buildings and natural changes such as sea-level rise or differences in vegetation are not annotated.The two images per city are selected to be cloud-free and are generally taken about 1-3 years apart.While there are 13 bands available in Sentinel-2 images, we restrict our focus to the RGB channels here.Although SiROC is able to handle channels outside of the visible spectrum as well our experiments show that the inclusion of the NIR band does not add value in the urban applications considered here.This may be different in vegetation monitoring where NIR bands tend to be more insightful.Spatial bands beyond RGB and NIR do not have a spatial resolution of 10m and are therefore excluded as well.
Beirut Harbor Explosion Dataset (BHED): On 4 th August 2020, a devastating explosion of large amounts of ammonium nitrate occurred in the port of Beirut in Lebanon.It led to over 200 deaths and left more than 300,000 people homeless because of heavy damages to buildings in the city. 1 We collect a pair of cloud-free Planetscope images with 3m per pixel resolution on 1 st and 5 st August before and after the explosion.We combine these images with ground truth on destroyed buildings provided by the Center for Satellite Based Crisis Information (ZKI) of the German Aerospace Center. 2 The building destruction map is based on manual annotation of very-high resolution images and field reports on the ground.Note that the annotations contain building destruction rather than building damage.Therefore, partial damages to buildings that withstood the explosion are not included.With this dataset, we aim to test the applicability of SiROC not only in medium but also in higher resolution images in problems where fast and accurate annotations are essential.
Agriculture Dataset: To test SiROC also outside the urban domain, we include two other test datasets from Saha et al. (2019) [11] as reference points.The first one, the Agricultural dataset, is a scene with bitemporal Sentinel-2 images from July 2015 over Barrax, Spain with 600 × 600 pixels in size.Between the two images is a time period of 10 days between which agricultural field activity changed notably.The reference map was manually annotated by the authors of Saha et al. (2019) [11].
Alpine Dataset: The second dataset consists of pre and post Sentinel-2 images of a fire in an alpine region close to Trento, Italy in spring 2019.A variety of other seasonal vegetation trends like ice and snow complicate this dataset.The scene has a size of 350 × 350 pixels with ground truth annotated manually by the authors of Saha et al. (2019) [11].

B. Competing Methods and Criteria
We compare our results to a variety of state-of-the-art unsupervised methods for optical CD in remote sensing.Since SiROC needs no training and does not rely on pre-trained neural networks its primary group of comparison consists of other image differencing-based methods.This makes SiROC fast compared to deep learning methods with comparable speed to traditional methods.We include several frameworks that improve on classical CVA.RCVA [29] incorporates close neighborhood information to make CVA more robust against misregistration.PCVA [35] uses change vector analysis of multi-level parcels to improve on CVA.DCVA is based on deep feature extraction with a deep neural network pretrained on imagenet [11].While DCVA was originally developed for high-resolution images, it is a resolution-agnostic framework relying on deep feature extraction from RGB channels.We also include a version of this method which we call DCVAMR specifically adjusted for medium-resolution, multispectral Sentinel-2 imagery [25] for the OSCD dataset.For BHED, we include DCVA, RCVA and PCVA as baselines as DCVAMR is not capable of handling Planetscope input channels.The most recent advancement in unsupervised CD for high-resolution imagery is Saha et al. (2020) [18] who employ self-supervised pre-training on remote sensing images in combination with a DCVA framework.We call this refined version of DCVA 'SSDCVA' and include it as the primary baseline besides general DCVA for BHED.
In line with previous evaluations on OSCD [25], [54], we analyze the performance of SiROC against the state-of-the-art in binary change segmentation based on specificity and sensitivity.Specificity is defined as the number of true positives (TP) over the sum of TP and false positives (FP): Specificity = T N T N +F P .Sensitivity is the number of TP over the sum of TP and false negatives: Sensitivity = T P T P +F N .This criterion is also known as recall.A method that is sensitive towards changes has a high sensitivity but a low specificity (and vice versa).A superior method should balance these objectives and evaluate better in both criteria.To further elaborate on the balance of change and no change class, we also report Precision = T P T P +F P and F1-Score = 2 * P recision * Recall P recision+Recall .

C. Results on OSCD
Parameters: We tune the parameters of SiROC on the OSCD training set resulting in the following parameter specifications: 1) Maximum neighborhood size: n max=200 2) Initial exclusion window: e start=0 3) Step size of ensemble: s=8 4) Filter Size of Morphological Operations: p=5 The maximum neighborhood size is 200 with stepsize 8. Contrary to the original idea of HSR in astronomy, we do not find it to be optimal to exclude direct neighbors of the pixel of interest from the analysis resulting in an initial exclusion window of zero.While neighboring pixels may be subject to the same kind of object-level change on the ground, they still can contribute important information if their weight is moderate.We find the best results with morphological opening and closing with a filter size of 5. We do not tune the To understand the origin of the performance difference to the previous state-of-the-art in more detail, we provide two ablation scores of SiROC.At first, we remove the morphological operations in SiROC.While morphological profiles help to transition to an object-level change mask, SiROC still exceeds previous unsupervised performance without them.No MP performance improves by 2 p.p. on specificity and 5 p.p. sensitivity vs. DCVAMR and by 4 p.p. specificity and 1 p.p. on sensitivity vs. DCVA.The resulting F1 score is 3-4 p.p. higher than deep learning-based methods and 5-10 p.p. higher than traditional methods here.To evaluate the effectiveness of ensembling, we also provide a score for a vanilla HSR with the same neighborhood size and no exclusion window.The Vanilla HSR performs slightly better but in the range of DCVA and DCVAMR with a specificity of 79.45% and a sensitivity of 70.24%.The F1 score is about 1 p.p. lower without ensemble voting.
Therefore, the majority voting mechanism is an effective tool to extract a more granular signal from the general HSR predictions.Further, the use of wide spatial context pioneered in astronomy is advantageous for change detection in remote sensing as well.In summary of Table I, SiROC sets a new state-of-the-art of unsupervised CD in medium-resolution images on OSCD.Even without morphological filters, the method still notably outperforms previous scores which points to a strong signal for change information in the original HSR method and the effectiveness of the majority voting mechanism.Combined with ensembling over different neighborhoods and morphological profiles, this exceeds previous quantitative results on the OSCD dataset.
Qualitative Results: The edge of SiROC compared to other unsupervised methods in the quality of change annotations for medium-resolution imagery is also visible when inspecting the predictions for specific scenes.Figure 2 displays exemplary change masks for Las Vegas.For SiROC, a threshold of v = 1/2 was used to obtain the images with a specificity of 95.28%, sensitivity of 78.75% and precision of 58.14% for this scene.Figure 2a visualizes the confidence of SiROC in the change propensity of a pixel as a heatmap from dark purple (0% votes) to yellow (100% votes).When comparing this to the ground truth on the bottom right, one can see that change confidence is strongly associated with the occurrence of a change.Figure 2b shows the binary change map after applying the threshold to the uncertainty map.Not only does SiROC pick up on the changed areas in the image, but it also fits the shapes of changing buildings fairly well.The visual similarities between Figure 2b and 2h are striking, especially compared to the other segmentations of competing methods.Also before applying the morphological operations, SiROC identifies the areas of interest in the image well although the predicted mask is naturally slightly more spurious.The morphological operations help to remove these spurious changes but the change signal in the predictions is in line with the ground truth (2c).DCVAMR is generally able to discover the changing regions of an image but struggles to identify the shapes of changing objects and rather fits round blobs (Figure 2d).DCVA tends to discover large changes and overestimate their size whereas smaller changes go undetected (Figure 2e).This might be related to the fact that DCVA was originally designed for high-resolution optical imagery in which building changes are larger in terms of pixel size.This is in line with the fact that DCVAMR, which is explicitly adjusted for Sentinel 2, tends to fit the size of changes better even though it also struggles with change shapes.PCVA and RCVA seem to extract building footprints rather than building changes here which leads to overcrowding of the segmentation mask.
A similar picture emerges when inspecting results for Dubai, which is a slightly more complex scene since the shapes of changes differ widely.SiROC detects changing regions again well but seems to struggle with the shape of changes in the upper part of the image.The newly constructed road is identified well.Consequently, the quantitative scores on this scene are slightly lower compared to the Las Vegas Scene with 86.87 % specificity, 76.61% sensitivity and 39.14% precision.The struggles of the competing methods are similar to the Las Vegas Scene: DCVAMR fits round shapes to any kind of change (Figure 3d), DCVA overestimates the size of large changes (Figure 3e), and PCVA extracts a spurious change map that rather looks like building footprints (Figure 3f).Therefore, the quality inspection of visual results confirms that SiROC obtains superior results on OSCD.
Uncertainty Estimation: To properly analyze if the confidence of SiROC also corresponds to well-calibrated uncertainties, we test this with calibration curves.For this, we split pixels into subsets based on the SiROC confidence and Average Specificity and Sensitivity on the OSCD training set when varying SiROC hyperparameters.The average scores are compared to the selected optimal selection of parameters and their training performance in the first line.
For single parameters (rows 2-5), we vary only the mentioned parameter on the grid given in the column range and leave the others at the selected optimum.
The column Evaluations gives the number of runs that were executed to obtain the average scores.Finally, in the last row "Joint" we vary all parameters on the given intervals in the rows above simultaneously.
analyze the respective performance for a level of confidence.
If the performance of SiROC is in principle increasing with the confidence, the uncertainty levels in fact correspond to the certainty of the prediction the model has. Figure 4 plots these confidence-performance curves for four cities in the OSCD test set.For all four cities, we see that model precision is non-decreasing in the confidence of the SiROC.Most of the time, the prediction increases notably in the confidence which means that SiROC not only performs well for this task but also returns well-calibrated uncertainties as part of its prediction.
Sensitivity to Hyperparamenters: To allow effective use of SiROC in practice, we offer a sensitivity analysis of the hyperparameter choice on OSCD along with recommendations for this choice in other applications.This sensitivity analysis is executed on the training set to avoid multiple evaluations on the test set.The results are shown in Table II.While the performance of the method naturally varies with the choice of hyperparameters, SiROC looks fairly robust against its hyperparameter choices.The first row gives the training set performance based on the selected parameters described in section C as a comparison point.Varying only the maximum neighborhood N max, the number of rows excluded e start and the stepsize s at the selected parameter specification influences the training performance marginally, at most.For all three parameters, the average specificity decreases while average sensitivity increases slightly.These three parameters essentially navigate how to group and prioritize neighborhoods.Excluding close context (e start), including more distant context (N max), and aggregating neighborhoods into larger groups (s) therefore does not seem to matter notably in practice to achieve good performance.The performance is slightly more sensitive towards the size of the morphological profile (p) where average specificity increases marginally and sensitivity drops by 9 p.p. if this is varied leaving other parameters untouched.Similarly, when varying all four parameters simultaneously in 75 random draws, performance drops with a difference of about 8 p.p. in sensitivity with similar specificity.In terms of magnitude, this performance drop is only a fraction of the difference between SiROC and its closest competitors  For each of these buckets, the performance is estimated separately.Since performance is generally non-decreasing in confidence, the uncertainty measure is considered well calibrated.on the OSCD test set.This implies that SiROC would likely outperform competing methods on this dataset for a variety of hyperparameter choices.
For potential other applications of SiROC in the future we suggest using the obtained parameter combination initially.This provides a starting point for further analysis in different contexts.Since the performance seems to be comparably susceptible to the size of the morphological profile, this parameter may deserve special attention during tuning.In the following, results on the remaining three datasets are obtained with this parameter combination which was the result of tuning on OSCD.Even though this may not necessarily give the best possible performance, we aim to validate that SiROC achieves convincing results in other applications without finetuning on single scenes.

D. Results on BHED
Quantitative Results: Table III displays specificity, sensitivity, precision and F1 scores on the scene.Generally, scores on BHED are higher than on OSCD since the changes are centered around the same area and have similar shapes.SiROC with default parameters achieves a specificity of 92.01% and a sensitivity of 83.38%.DCVA achieves a similar specificity with 91.87% but falls short in terms of sensitivity by about 4 p.p. with a score of 79.85%.SSDCVA places slightly When we adjust the scale of morphological operations to 10, SiROC performs significantly better which suggests that there may be notable tuning potential for higher resolution inputs.Still, the baseline parameters perform well on this scene.Therefore, SiROC demonstrates its usefulness beyond medium-resolution images and can also be used in conjunction with high-resolution images for CD.The other ablation scores again point towards the most important steps within SiROC to achieve this performance.Without morphological profiles, the scores of SiROC drop about 4-9 p.p. in all four categories.Still, it achieves slightly superior precision and F1 scores but falls short to DCVA with a difference of about 3 p.p. in specificity and similar sensitivity.This is a notable difference to medium-resolution imagery on OSCD where the exclusion of morphological filters decreased the performance of SiROC but it was still superior to DCVAbased methods.This is not necessarily surprising since deep learning-based methods tend to relatively improve their CD performance compared to traditional methods with increasing spatial resolution.Without majority voting over different neighborhoods, the Vanilla HSR version performs better but in the range of RCVA and PCVA.Again, it is the combination of HSR, ensembling over different neighborhoods and transitioning to the object level with morphological operations that all contribute significantly to the overall performance of SiROC.
Qualitative Results: Figure 5 shows visual comparisons of the discussed methods on BHED.The first row of images presents the pre-explosion image (Figure 5a), the post image (5b) and the ground truth (5c).The heart of the explosion in the port can be found in the middle of the image with almost the entirety of buildings completely destroyed around it.Panel 5d presents the binary SiROC segmentation with baseline parameters obtained on OSCD.While SiROC is missing some granularity in its segmentation of destroyed building footprints, the changing areas are well identified with few false positives outside of the port.For a larger morphological filter size (p), the main area is identified more densely with better quantitative results but the shapes of buildings vanish (5e).Without morphological operations, the core change is still well-segmented although the number of false positives in the outer regions of the image increases (5f).SSDCVA shows similar tendencies to summarize the port area as one large change with a number of spurious false positives (5g).DCVA shows fewer salt and pepper noise than SSDCVA here and generally segments the exploded buildings similar to SiROC, however, with a slightly more perforated shape (5h).The segmentation by RCVA is not really competitive here since the maps are spurious and changes are not well identified (5i).Results for PCVA are similar to RCVA and hence omitted.

E. Results on Agriculture Dataset
Quantitative Results: Table IV displays the results of SiROC and competing methods on the agriculture scene.SiROC is applied to the dataset with the parameters obtained on Onera without further adjustment.Hence, the results we provide are a validation exercise in the different context of non-visible parts of the spectrum without parameter finetuning.
To be consistent with previous evaluations on this dataset [11], we compare SiROC with PCVA and RCVA based on vegetation (VEG) and near-infrared (NIR) channels of Sentinel-2 as inputs.The score for DCVAMR is based on the full Sentinel-2 input images as the method was deliberately designed to incorporate all channels.
While SiROC achieves the top score in terms of specificity and precision with 90.81% and 74.23% respectively, it falls short of DCVAMR on sensitivity (88.70% vs 94.26%) and F1 score (80.85% vs 81.47%).DCVAMR seems to lean slightly more towards the change class whereas SiROC rather classifies a pixel as no change in unclear cases.SiROC is superior to PCVA and comparable to RCVA in performance for both, VEG and NIR channels as inputs.
The ablation scores underline that morphological profiles still help although the effects are smaller than in urban applications with an average difference in about 1-2 percentage points in all four criteria.Further, excluding the majority voting mechanism does not hurt performance but actually improves it slightly here.The vanilla HSR performs slightly worse but in the range of RCVA and better than PCVA on its own.Smaller benefits of including majority voting and morphological profiles could be linked to the fact that parameters for these operations were tuned in an urban RGB context.Although already quite effective, the accuracy of SiROC could likely be further improved with parameter finetuning.
Qualitative Results: Figure 6 presents pre and post RGB images (6a-6b), the ground truth 6c, and change predictions (6d-6h).The visual impression of change predictions confirms the quantitative results.Predictions are fairly accurate on this scene which suggests a comparably easy task relative to the more complex OSCD scenes.SiROC segments changing regions well and struggles with the varying field shapes only in rare instances.Similarly, the results of DCVAMR and RCVA are also fairly accurate with a slightly higher tendency to predict the change class.In comparison, the mask by PCVA produces some false positive regions.Overall, SiROC shows similar performance to highly effective methods also in the agriculture domain.

F. Results on Alpine Dataset
Quantitative Results: Results for the Alpine dataset can be found in Table V.Even though SiROC reaches the highest specificity, it does not quite pass the overall performance of DCVAMR from [11] on this dataset.Nevertheless, SiROC ranks highly also in sensitivity, precision and F1 score, particularly based on NIR inputs with total scores of 98.92%, 75.71% , 52.28% and 61.85%.RCVA with NIR inputs is comparable in performance but SiROC is the only method  that makes effective use of SWIR inputs compared to PCVA and RCVA.The ablation scores underline the effectiveness of morphological transformations with about a 20 pp. drop in F1 score compared to SiROC for both NIR and SWIR.Removing the ensembling leads to a notable drop in F1 scores, particularly with NIR inputs.
Qualitative Results: Figure 7 plots prediction masks for selected models for the Alpine dataset.In Panel 7a, the false color composite shows the annotated area of change affected (7b) by a fire in purple on the right.SiROC identifies this well although it is tempted to also classify a small number of false positives as change.While it is hard to control for seasonality in a bitemporal setting, SiROC (NIR) (7c) still excludes most other vegetation updates which are not the result of actual change here.The morphological profiles help on this scene to exclude spurious predictions (7d).Compared to SiROC (SWIR) (7e), SiROC (NIR) segments the changing area slightly better although the shape is identified more clearly by DCVAMR (7f).PCVA (NIR) (7g) seems to struggle slightly more with the shape of the burned area whereas the results of RCVA (NIR) (7h) look similar to the results of SiROC (NIR) which is in line with the quantitative scores of Table V.

IV. DISCUSSION
SiROC is an effective method for CD in medium-as well as high-resolution optical imagery which achieves competitive performance on four datasets.In the following, we elaborate on the intuition of SiROC's performance.When contrasting SiROC to image differencing methods, SiROC can be interpreted as an improvement over standard image differencing techniques because it does not assume the same changes in the acquisition conditions across time for the whole image.Rather, it allows for local changes in acquisition conditions.In standard CVA or RCVA, for example, an implicit assumption is that changes in the acquisition conditions across time affect each pixel similarly.SiROC releases this restriction and instead allows for local trends in regions of the image.If a pixel deviates from the local trend around it, it is likely to undergo a change in SiROC.In RCVA or CVA, one would compare this pixel against trends in the whole image and not against its surrounding only.This might be unrealistic in complex scenes where pixels values highly depend on local trends in the surroundings.This is for example the case when a new building casts a shadow on a previously illuminated pixel.Similarly, a cloudy pixel in t + 1 that was unobstructed in t might not necessarily be changing and is rather influenced by the local trend of a cloud rather than general image trends if large parts of the image are not obstructed by clouds.Hence, SiROC allows for a more granular analysis of deviations from trends in an image time-series because compared to previous methods it makes full use of multi-temporal information in close as well as distant neighbors.Although we compare our results to deep learning-based methods, our intention is rather to augment these models than replace them, especially with high-resolution images.SiROC provides an efficient and accurate way to obtain change labels that could also be infused into deep learning models.One application of SiROC could be in self-supervised learning where pseudo-labels are often obtained based on traditional image differencing techniques such as CVA [55].SiROC is not only superior in performance compared to image differencing.It also comes with a built-in, well-calibrated uncertainty of predictions.This could be especially beneficial in self-supervised settings since it automatically allows to discriminate pseudo-labels by confidence.For example, one could train only based on pseudo labels with high certainty and discard uncertain data points.Similarly, in some unsupervised methods such as MSDRL for VHR imagery, an initial pseudo-classification is separated by confidence where high confidence examples are used for training a classifier that subsequently obtains predictions for leftover uncertain pixels [19].In these methods, SiROC could also be used to obtain initial predictions and uncertainties to potentially improve not only the initial classification but maybe also the uncertainty categorization.The combination of deep learning-based methods and SiROC may hence open up new potential for CD methods.While we restrict our focus to CD with optical images here, the framework of SiROC may be extended for applications on other multitemporal CD problems in remote sensing as well.

V. CONCLUSION
We present SiROC, an efficient and accurate unsupervised method for CD in medium-and high-resolution optical images.SiROC is inspired by HSR which is used for exoplanet search in astronomy.It models a pixel of interest in t as a linear combination of its neighbors and applies this model to t + 1 to obtain a prediction for the pixel based on its neighbors.The difference of the prediction for t + 1 and the actual pixel value in t + 1 is interpreted as the change signal.If the prediction is far from the actual value, trends in the neighboring pixels divert from the difference in the pixel of interest over time which is seen as an indicator for change on the ground.
We refine and extend HSR in two major ways to apply it to optical satellite images as SiROC.First, we iterate over several, mutually exclusive neighborhoods and apply HSR with all of these neighborhoods as input to obtain a distribution of change predictions.We combine these predictions with majority voting which improves performance significantly and also returns a heatmap of votes per pixel which can be interpreted as a well-calibrated uncertainty.Second, we use morphological opening and closing at one spatial filter scale to transition from pixel-level to object-level predictions.
The results of SiROC are validated on four datasets.For urban change detection with medium-resolution images, we verify the effectiveness of our method on OSCD, which contains binary change annotations for 24 cities across the globe.SiROC sets a new state-of-the-art for unsupervised CD on OSCD which surpasses previous methods by 10 p.p. in terms of specificity, 2 p.p. in sensitivity, 11 p.p in precision, and 13 p.p. in F1 score.We further validate the performance of SiROC on high-resolution images with a dataset on the Beirut Harbor Explosion (BHED).Also on this dataset, SiROC surpasses the performance of competing methods and underlines its abilities to segment urban change accurately at several scales.Further, we provide two validation exercises on non-urban data with Sentinel-2 inputs.SiROC segments the effects of a fire in the Italian Alps accurately and in the range of competing methods.On the Agriculture dataset, SiROC falls short of DCVAMR in overall scores but still identifies the changing crop activity correctly.
While SiROC compares well against current deep learningbased unsupervised methods in CD, SiROC should rather be seen as a complement than a substitute to these methods.Since it provides an accurate way to predict change signals with a built-in, well-calibrated uncertainty it may be especially useful in conjunction with deep learning-based methods to generate pseudo-labels.Although we apply SiROC primarily to changes with multispectral data, the model may be applicable to other change detection problems as well which we plan to explore in future research.

ACKNOWLEDGMENT
We are grateful to the Center for Satellite Based Crisis Information (ZKI) of the German Aerospace Center for providing us with ground truth of the Beirut Explosion scene.

Fig. 1 .
Fig. 1.Half-Sibling Regression (HSR): Figure 1a: A set of neighbors is fit to a pixel of interest as a linear combination at time t).Figure 1b: At t + 1, the pixel values of the neighbors are used together with the coefficients obtained at time t to predict the pixel of interest in t + 1.Figure 1c: The predicted pixel values are compared with the actual pixel values at t + 1 to obtain a change signal.
Fig. 1.Half-Sibling Regression (HSR): Figure 1a: A set of neighbors is fit to a pixel of interest as a linear combination at time t).Figure 1b: At t + 1, the pixel values of the neighbors are used together with the coefficients obtained at time t to predict the pixel of interest in t + 1.Figure 1c: The predicted pixel values are compared with the actual pixel values at t + 1 to obtain a change signal.

Fig. 2 .
Fig. 2. Qualitative Comparison OSCD -Las Vegas.Figure 2 visualizes the number of change votes per pixel in SiROC (2a) and the corresponding binary predictions after (2b) and before morphological operations (2c).Competing models are in 2d-2g and the ground truth in 2h.SiROC predicts change regions and shapes of the ground truth well while competing methods struggle either with identifying the shapes (2d & 2e) or the areas of change (2f & 2g) for the Las Vegas Pair.

Fig. 4 .
Fig.4.Confidence-Performance Plots on four cities of the OSCD Dataset.On the x-axis, point in the image are sorted into buckets by SiROC confidence.For each of these buckets, the performance is estimated separately.Since performance is generally non-decreasing in confidence, the uncertainty measure is considered well calibrated.

Fig. 5 .
Fig.5.Qualitative Comparison Beirut Explosion.Figure5shows the Planetscope image pair (5a & 5b), change ground truth (5c) and model predictions (5d -5i) for the Beirut Explosion scene.SiROC with main parameters (5d) identifies the area of destroyed buildings around the epicenter correctly with few false positives although the shapes are lacking some granularity.Increasing the size of morphological operations improves accuracy but tends to fit one large blob with missing building shapes 5e.Excluding morphological operations increases false positives although the main changes in the center are still idenfitied well (5f).Competing methods struggle not only with the shape of change but also detect a number of false positives far away from the explosion (5g-5i).

Fig. 6 .Fig. 7 .
Fig.6.Qualitative Results Agriculture Dataset.Pre (6a) and post RGB (6b) image with changing agricultural fields.The ground truth (6c) shows high similarity to SiROC (6d) but also to DCVA (6f) and RCVA (6h) while PCVA (6g) has some false positive areas and spurious change predictions.Change segmentations without morphological profiles 6e for SiROC still works well which is line with the quantitative results of TableIV.
Table I reports specificity and sensitivity scores of SiROC and competing methods on the OSCD test set.Scores are averaged on the city level.SiROC with a voting threshold of v = 1/2 achieves a specificity of 88.31% with a sensitivity of 70.71% and 24.80 % precision as well as 36.72%F1-Score.This is a high score in all four categories by a significant margin.The difference to DCVAMR is about 6-13 percentage points (p.p.) depending on the category.DCVA achieves a sensitivity that is slightly below but close to SiROC but lacks behind in specificity, sensitivity and F1-score by more than 10 p.p. Compared to the best results of methods without deep learning-based feature extraction, SiROC gains about 12 p.p. in specificity, 7 p.p. in sensitivity 12 p.p. in precision and 15 p.p. in F1 on RCVA.