Improved Difference Images for Change Detection Classifiers in SAR Imagery Using Deep Learning

Figure 2 illustrates the overall architecture of the proposed method. It replaces the conventional method that is illustrated in Figure 1 image differencing step. The idea of the proposed method is to improve the SAR image comparability by considering the acquisition conditions of the SAR images. The proposed method utilizes a mapping transformation function that creates artificial SAR images in the requested acquisition conditions. The mapping transformation function $\mathcal{F}$ is a neural network model that is trained to predict the SAR image $I_{t}$ at the time $t$. The neural network output $\hat{I}_{t}$ is the artificial SAR image that is created in the acquisition conditions of $I_{t}$, therefore it should be more comparable to the $I_{t}$ than previous SAR images from the location that might have been captured in different acquisition conditions.

The model input consists of three distinct features, which are: The previous SAR images from the location; the acquisition conditions of the SAR images (including at time $t$); and the digital elevation map from the location. The objective of the neural network model is to learn to replicate the SAR image at the time $t$. The only information from the time $t$ in the model input are the image acquisition conditions of the $I_{t}$. This means that for the model to be able to replicate the $I_{t}$, it needs to learn to map the information contained in the previous SAR images and the digital elevation map to the image acquisition conditions of the $I_{t}$. With an ideal model that could perfectly replicate the $I_{t}$, the $\hat{I}_{t}$ and $I_{t}$ would be identical if nothing has changed between the image acquisition of the $I_{t-1}$ and $I_{t}$, however the $\hat{I}_{t}$ would be missing the change if something had changed after the previous image acquisition since the information of the change is not included in the model input data. In practice the SAR images include random noise that is impossible to replicate accurately, and the acquisition conditions are not accurate enough for perfect replication of the $I_{t}$, therefore the $\hat{I}_{t}$ only approximates the $I_{t}$.

The objective function of the model training to produce artificial SAR images is defined as follows:

The intuitive description of the $\hat{I}_{t}$ is that the neural network-based mapping transformation function produces a prediction how the $I_{t}$ should look like based on previous information about the location and the actual imaging conditions of the $I_{t}$. The produced image $\hat{I}_{t}$ can be used with the actual image $I_{t}$ to create the difference image $\hat{I}_{DI}$ by using a simple algebraic operation like subtraction, ratio, or log ratio. Generating the difference image is the standard method of conducting change detection, especially when using unsupervised methods [17].

Conventional methods of producing the difference image often use only one of the previously captured images with the most recent image to generate the image e.g. $I_{DI}=g(I_{t},I_{t-y})$ [26]. This method has the previously discussed drawbacks of noise and imaging conditions affecting the final difference image quality. By using the proposed mapping transformation function, the predicted image $\hat{I}_{t}$ is used in the place of the previously captured image to generate the difference image e.g. $\hat{I}_{DI}=g(I_{t},\hat{I}_{t})$. Recall that $\hat{I}_{t}$ is a representation of $I_{t}$ based on geospatial information from the time $t-1$ and earlier, therefore it is missing all the changes that have happened after that time. The predicted image $\hat{I}_{t}$ does not contain noise and the mapping transformation function can correct the acquisition condition mismatch between the images, therefore the proposed method should produce better quality difference images when comparing it to the conventional method.

SAR imaging is sensitive to the soil moisture content of the imaged area [19]. Change in the soil moisture level changes the dielectric constant of the soil, and that way changes the SAR backscatter intensity. Often the soil moisture content changes should be ignored by the change detection system. Otherwise, the system would notify changes after every rainy day. This is one of the advantages of the proposed method. By adding weather to the model input acquisition condition parameters, the mapping transformation function can learn to construct the $\hat{I}_{t}$ in the actual weather conditions of $I_{t}$ and should correctly model the changes in the soil moisture changing the backscatter intensity. Therefore, the false positive changes, that are potentially caused by weather condition changes, should be reduced.

In addition of weather, the acquisition condition parameters also include the imaging angle and identify the satellite that captured the image. A location is imaged by one of the sentinel satellites with an interval ranging from a few days to about a week. The satellite does not capture the image from the same angle every time. The satellite can be in ascending or descending orbit during the image acquisition and the incidence angle can vary between the overpasses. The ascending or descending orbit changes the look direction of the satellite, and that way has a considerable affect to the resulting image. The Sentinel-1 satellites are right-looking. When the satellite is descending from North to South it is imaging to the direction of West, and for ascending passes it is imaging to the direction of East [27]. Various 3D features, like forest edges, lake banks and hills are sensitive to the look direction, therefore the imaging angle is an important parameter when computing the difference image. When using an image differencing method where only one previous image is used for difference image computation, the imaging angle of the most recent image can restrict what previous images can be used to produce the difference image. Seasonal changes, like foliage growth or change in snow cover, means that the most optimal image for the differencing would be the most recent previous image, however different imaging angles can limit the usage of the most recent images. This problem is not present with the proposed method. The model input includes $n$ previous images and their imaging angle information. The model output image $\hat{I}_{t}$ is produced using the actual acquisition conditions of $I_{t}$. The model can use all the information from all $n$ input images, despite the input including images from different look directions, and the produced image $\hat{I}_{t}$ represents an image that is acquired from the same angle as $I_{t}$.

Figure 3 illustrates the architecture of the neural network-based mapping transformation function. The architecture is based on the well-known U-Net neural network architecture [28]. The previous $n$ SAR images, and the digital elevation map (DEM) are stacked to construct the input. The previous images and the DEM are all from the same location. The images are projected to the same resolution and the pixels across the different images are aligned to match the same geographical position. The U-Net architecture is constructed from encoder and decoder units. The encoder takes the input and compresses the input image stack to the latent space by using a set of downsampler blocks that half the input resolution using convolution layers with stride $2\times 2$. The encoder stacks enough downsampler blocks so that the input image stack is compressed to $1\times 1$ resolution in image height and width dimensions. The image acquisition conditions vector, that contains the information of the acquisitions conditions for the $n$ input images and the target image, is concatenated to the latent vector as described at the end of section II-C. The resulting vector is then fed to the decoder that decodes the vector back to the dimensions of a normal SAR image outputting the $\hat{I}_{t}$. The decoder is constructed from upsample blocks that double the width and height dimensions using transposed convolution layers with stride $2\times 2$. The decoder has same amount of upsampler blocks as the encoder has downsampler blocks. The number of filters, that are used in the upsampler and downsampler blocks, can be configured for every block individually, except for the final upsample block that has the same number of filters as the SAR image has bands. The encoder and decoder layers are connected with skip connections that help the model in producing the output by not forcing the model to pack all the information to the latent vector. Instead, the information can flow from the input to the output by skipping most of the layers in the architecture. This is a standard method in U-Net style architectures.

A dataset is needed for the training of the neural network-based mapping transformation function. As discussed previously, the mapping transformation function input is composed from the previously taken SAR images; the acquisition conditions of the previous and the most recent SAR image; and the digital elevation map from the location. The objective of the model is to learn to predict the most recent SAR image based on the input, therefore the most recent SAR image is the target in the training dataset. This means that the training dataset does not require any labelled data making the learning process of the proposed method unsupervised and economical to implement. The dataset can be generated directly from available data sources without the need of human labelling for the data. The dataset is available at the Fairdata.fi service [29].

The SAR images for the dataset were acquired from the ESA Copernicus Open Access Hub [30]. High resolution Interferometric Wide Swath (IW) Ground Range Detected (GRD) products were used in this study [31]. The images were captured between March 2020 and August 2021 from the area illustrated in the Figure 4. All images from the time frame that included the area were downloaded from the Copernicus Open Access Hub. The images were preprocessed using the Sentinel-1 Toolbox from the Sentinel Application Platform (SNAP) [32], by applying the data preprocessing workflow described by Filipponi in [33]. The optional noise filtering step was applied to the dataset using the Refined Lee filter from the SNAP toolkit. The more accurate AUX_POEORB precise orbit files were used in the Apply Orbit File step. The AUX_POEORB files are available 20 days after the image acquisition [34], and since the processing was done in spring 2022, the more accurate orbit files were available for all images. The proposed workflow in [33] uses the SRTM Digital Elevation Database in the Range Doppler Terrain correction step, however the database does not cover the area from where the dataset was created, therefore the Copernicus 30m Global DEM was used that does cover the area. The SNAP toolkit can automatically download the required DEM files during preprocessing and the Terrain Correction step supports multiple different DEM sources, including the Copernicus 20m Global DEM, thus the change was trivial to implement. The preprocessed images were saved as GeoTIFF files and uploaded to PostgreSQL ²²2https://www.postgresql.org/ database that was using the PostGIS ³³3https://postgis.net/ extension. Using a relational database as the storage backend simplified the dataset generation process since all the data was available in one place and queryable with SQL.

Although the Copernicus 30m Global DEM was used in the SAR image terrain correction preprocessing step, the product was not used for the mapping transformation function input. Instead, we used more accurate DEM from National Land Survey of Finland (NLS). NLS provides the DEM in multiple different resolutions of which the most accurate 2m grid DEM was used [37]. The data is open access and distributed under Attribution 4.0 International (CC BY 4.0) license ⁴⁴4https://creativecommons.org/licenses/by/4.0/. The DEM was downloaded in GeoTIFF format and uploaded to the same PostgreSQL database with the SAR images.

As discussed before, the image acquisition condition data included information about the weather when the images were captured. This data was acquired from Finnish Meteorological Institute (FMI) that provides daily weather observations that are interpolated to $1\times 1$ km grid [38]. The interpolation method is described by Aalto et al. in [39]. The data is distributed in NetCDF format and uploaded once a month. Daily mean temperature, daily precipitation sum, and snow depth data was downloaded from the time range. The daily observations were extracted from the NetCDF files, converted to daily GeoTIFF rasters, and uploaded to the same PostgreSQL database with the SAR images and DEM.

The final data samples were created by sampling random locations from the area and random dates from the time range. For training dataset, the time range was limited to the time before 20th of June in 2021, and for the test dataset the time was limited after the date. An assumption was made that the samples do not have any changes between the acquisitions $I_{t-1}$ and $I_{t}$. This assumption is not likely true for all of the samples, however the total dataset size is created to be large enough so that the samples that have changes between the two last acquisition should be marginally small fraction of the total dataset and neural networks can adapt to noisy data [40]. The image size was set to $512\times 512$ pixels, and number of previous images was set to $4$. The spatial resolution of a high resolution IW GRD product is $20\times 22$ meters, and the images are distributed with $10\times 10$ meter pixel spacing [41]. The geographical dimensions of the images were set to $3\times 3$ km making the pixel size $3000/512\approx 5.9$ meters. This is higher resolution than the original 10 meter pixel size, therefore the information loss is minimized during processing. For each random location and date, the target SAR image $I_{t}$ was the next SAR image from the location that was available after the date. The input SAR images $I_{t-4},I_{t-3},I_{t-2},I_{t-1}$ were the SAR images from the four previous acquisitions from the location that were captured before the $I_{t}$. The SAR images and the DEM was queried from the PostgreSQL database and the rasters were projected to the same projection window with the same $512\times 512$ resolution and $3\times 3$ km spatial dimensions using GDAL library [42]. The gdal.Translate function was used for the projection with nearest neighbor resampling algorithm. After the projection, all pixels were geographically aligned across all images and the images could be stacked to construct the input image stack. The Sentinel-1 satellites use Interferometric Wide swath mode with dual polarization over the land areas thus one SAR image has two bands [12]. Both bands are used in all input images and the target image. That makes the input image stack to have $1+4\cdot 2=9$ channels (DEM has one channel and every SAR image has two bands/channels), and the model output image has two bands.

The acquisition conditions were composed from the following features:

1. 1.

   Mean temperature of the acquisition date

2. 2.

   Snow depth in the acquisition date

3. 3.

   Satellite orbit direction during the acquisition (Ascending/Descending)

4. 4.

   Incidence angle

5. 5.

   Satellite id (Sentinel-1A or Sentinel-1B)

6. 6.

   Precipitation amount in the acquisition date and three previous dates

In addition of imaging conditions, such as weather and imaging angle, the satellite that captured the image is also added to the acquisition condition vector. The satellite id is encoded to $1$ if the image is captured by the S1A satellite and $0$ if the image is captured by the S1B satellite. The imaging instrumentation is not necessarily identical in both satellites and the model might learn to use this information to create more accurate images. All other features were scalar values from the acquisition date except for precipitation that is a vector with values for four different days. Since the moisture content of the soil has known effect to the signal, and moisture can linger long times in the soil, it was decided to include the precipitation amounts from multiple days to the acquisition conditions. Taking the precipitation amounts from the previous $4$ days was a somewhat arbitrary decision with a reasoning that the neural network can learn to ignore the precipitation amounts from previous days if they have no use. The features were flattened to the final vector with dimensionality of $|D|=9$.

The final generated dataset had had around $230,000$ training samples, and around $9,000$ test samples.

The performance of the proposed method was measured using experimentation. The main contribution of this paper is to offer a new strategy for computing the difference image. Existing methods generally use a strategy where the difference image is computed using $I_{DI}=g(I_{t-y},I_{t})$, where the $g$ is the differencing function, $I_{t-y}$ is one of the previous images from the location captured at some previous date, and $I_{t}$ is the most recent image from the location. The proposed method uses the neural network output $\hat{I}_{t}$ in place of the $I_{t-y}$ to compute the difference image $\hat{I}_{DI}=g(\hat{I}_{t},I_{t})$. The mapping transformation function factors in the imaging conditions of $I_{t}$ when generating the $\hat{I}_{t}$, therefore the $\hat{I}_{DI}$ should be higher quality when compared to $I_{DI}$. The difference image is generally further used in the change detection system to detect the changes by applying a classifier to the difference image. The classifier outputs a change map indicating the pixels that contain the detected changes. By using identical classifier to classify the difference images generated by the two different methods and comparing the classifying accuracy of the resulting change maps, the quality of the two difference images can be measured.

Different neural network parameters were experimented with, and the best results were achieved with the parameters shown in the Table I. Mean squared error was used as the loss function, and AdamW [51] was used as the optimizer. The neural network architecture was implemented using TensorFlow deep learning framework [52]. The training was conducted on one NVIDIA V100 GPU with batch size of $200$, and training time of around 30 hours.

Figure 6c demonstrates the model performance for one of the test samples. Figure 6a shows the real SAR image that the model tries to predict. Figure 6b illustrates the difference between the real SAR image and the model output with a heat map where lighter color indicates a greater error. The predicted image is very close to the real SAR image except for lack of noise that is purely random and impossible for the model to predict. Likewise, the lower right corner of the image has an area that has greater error in the prediction. The error is located in a lake, therefore the error can be a result of waves that are likewise impossible to predict.

The proposed method depends on that the mapping transformation function adapts the predicted image $\hat{I}_{t}$ based on the imaging conditions of $I_{t}$. To verify that the model genuinely uses the image acquisition conditions to produce the $\hat{I}_{t}$, the model was experimented to produce outputs with manually modified imaging condition vector $D_{t}$. Figure 6d and Figure 6f image pair illustrates model outputs where the $D_{t}$ is modified to have opposite orbit directions. Figure 6e illustrates the difference between the images. The lake banks and the upper left corner of the image, where there is a small hill, have large differences between the two generated images. All locations, where there are greater differences between the images, are 3D features. The Sentinel-1 satellites have different look directions on ascending and descending orbit directions. Therefore, the scattering of the radar signal is different and the difference is most noticeable on 3D features. Since the differences are so clearly located on the 3D features in the image the model is clearly factored in the orbit direction when generating the output. This verifies that the imaging conditions are used by the model to produce the $\hat{I}_{t}$ in the imaging conditions of $I_{t}$.

The same experiment was conducted by modifying the precipitation amounts in Figure  6g and Figure 6i. The difference between the generated images is shown in the Figure 6h. This time the difference between the generated images is focused on swamp, meadow, and agricultural land areas in the image. The forest areas have only small differences between the images. In forest areas, the radar signal is scattered back by the forest foliage where the moisture does not affect the scattering properties as much as the open areas. In open areas, the radar signal hits the ground where the soil moisture content is altered more by the rain, thus changing the backscatter intensity. This experiment suggests that the model uses the precipitation information correctly when generating the output image.

The conventional method of computing the difference image is to use one of the previous SAR images that is captured at some preceding date with the most recent image to produce the difference image $I_{DI}=g(I_{t-y},I_{t})$. There are multiple different strategies when selecting the previous image. The simplest strategy is to select the previous image that is preceding the image that was captured most recently. This strategy has the advantage that the least amount of time has elapsed between the images, therefore the number of natural changes, like foliage growth or soil moisture changes, are minimized. However, the problem is that the previous image has very likely different incidence angle and it might have been captured from different orbit direction (ascending/descending). To make sure that we compare the proposed method to the best conventional method, three different previous image selection strategies were compared to identify the best strategy. The threshold classifier was used to compare the quality of the difference images that were produced using the different strategies. The strategies have different trade offs between the elapsed time and imaging angle:

1. Method 1:

   Closest incidence angle and the same orbit direction.

2. Method 2:

   Most recent previous image with the same orbit direction.

3. Method 3:

   Most recent previous image preceding the target image ($I_{t-1}$).

Figure 7 illustrates the comparison of the three different methods using ROC curve plots. Table II shows the results in a list format by displaying the AUC metrics. The strategy where the previous image is captured from the same orbit direction and has the closest incidence angle with the $I_{t}$ is the best $I_{t-y}$ selection strategy. From this on forward, the Method $1$ is always used when referring to the conventional method of computing difference image.