Leveraging autoencoders in change vector analysis of optical satellite images

Various applications in remote sensing demand automatic detection of changes in optical satellite images of the same scene acquired over time. This paper investigates how to leverage autoencoders in change vector analysis, in order to better delineate possible changes in a couple of co-registered, optical satellite images. Let us consider both a primary image and a secondary image acquired over time in the same scene. First an autoencoder artificial neural network is trained on the primary image. Then the reconstruction of both images is restored via the trained autoencoder so that the spectral angle distance can be computed pixelwise on the reconstructed data vectors. Finally, a threshold algorithm is used to automatically separate the foreground changed pixels from the unchanged background. The assessment of the proposed method is performed in three couples of benchmark hyperspectral images using different criteria, such as overall accuracy, missed alarms and false alarms. In addition, the method supplies promising results in the analysis of a couple of multispectral images of the burned area in the Majella National Park (Italy).


Introduction
The Earth's surface is constantly changing due to anthropogenic and natural causes like the progression of desert areas, deforestation, glacier movements, fires or earthquakes (Alberti et al., 2003). Monitoring these changes over time may provide valuable information on the transformation of the Earth's environment paving the way for a better decision policy on the minimisation of the risk of disasters (Michel et al., 2012). In particular, the rapid development of multispectral (MS) and hyperspectral (HS) technology has recently unleashed the potential of change detection (CD) methods in a wide range of remote sensing applications Also in this study, we consider autoencoders for data restoring. In particular, we use the output of the decoder function of the primary image-specific autoencoder to restore both the primary image and the secondary image. Note that this is not directed to operate data denoising only. In principle, the autoencoder trained on the primary samples can contribute to recovering denoised samples of the pixels in the primary image, as well as denoised samples of the unchanged pixels in the secondary image, but it should see changed samples of the secondary image as anomalies, and so reconstruct them badly. So, the idea is to exploit autoencoders to disclose patterns that better delineate the pixels of the sensed scene where a change has occurred over time. We take advantage of these patterns by completing the CVA on the restored data (rather than on the original, acquired data). In particular, we compute the Spectral Angle Mapper (SAM) distance pixel-by-pixel between the restored spectral data vectors of the primary and secondary image, respectively. This distance quantifies the spectral change range at each pixel of the scene. The Otsu's algorithm (Otsu, 1972) is, finally, adopted to separate the foreground regions, where a change occurred, from the unchanged background.
We evaluate the proposed method performing the CD analysis of various, benchmark, bi-temporal, co-registered HS images collected in various urban and rural scenarios. As the change information is available on these datasets, the empirical study can verify the accuracy of the proposed CD method. In addition, we perform the evaluation of the viability of the proposed method in delineating the burnt area of bi-temporal, co-registered MS images acquired with Sentinel-2 in the area of Majella National Park (Italy).
The paper is organised as follows. The related works are presented in Section 2. The basic concepts are introduced in Section 3, while the proposed CD method is illustrated in Section 4 and the implementation is described in Section 5. The findings in the evaluation of the proposed strategy with benchmark HS data are discussed in Section 6, while the achievements in the analysis of the MS data of the burnt area in the Majella National Park (Italy) are illustrated in Section 7. Finally, Section 8 draws conclusions and proposes future developments.

Related work
Since obtaining a large number of labelled samples for supervised training is usually timeconsuming and labour-intensive, remote sensing research devotes significant research effort for the formulation of CD methods in unsupervised machine learning.
In the unsupervised machine learning paradigm (Hussain et al., 2013;Bruzzone & Prieto, 2000), changes are commonly detected by resorting to the CVA strategy that computes a measure of similarity (or distance) between co-located pixels of a couple of images and uses a threshold-based approach to identify a distance threshold to separate changed pixels from the unchanged background. Various similarity (or distance) measures have been investigated for CVA methods (e.g., Appice et al., 2020;Falini et al., 2020;Seydi & Hasanlou, 2017;Yang & Mueller, 2007). The threshold to detect the changes is estimated by resorting to the spectral data (i.e. in a data-driven manner) (Lu et al., 2010;Najafi et al., 2017;Penglin et al., 2012) by leveraging probabilistic information extracted from the distribution of the (distance or similarity) measure among the pixels. A well-known approach commonly used for the threshold determination is Otsu's method (Otsu, 1972;Sahoo et al., 1988). In (López-Fandiño et al., 2019), Otsu's algorithms is evaluated in combination with SAM and Watershed algorithm. Alternatively, clustering algorithms are adopted Appice et al., 2020), in order to separate distances (or similarities) of changed pixels from the unchanged background.
Algebra-based methods, similarity-based methods, as well as distance-based methods belong to the threshold-based family of CD approaches. In particular, algebra-based CD methods use mathematical operations (such as image differencing or image ratio) on images taken at different times to generate a change matrix output (Ilsever & Unsalan, 2012). Similarity-based CD methods resort to the computation of a similarity measure (e.g., correlation measure) between a pair of spectral vectors (Choi et al., 2010). Distance-based CD methods are founded on a spectral distance measure (e.g. SAM, Z-score Information Divergence) computed between the spectral vectors on corresponding pixels (Choi et al., 2010). A recent study  adopts a spectral-spatial distance for CVA. In addition, it introduces an iterative upgrade of the traditional distance-based approach by accounting for a representation of the possible change iteratively learned through classification. Due to the lack of ground change information, classification is supervised with pseudo-labels yielded on the spectral-spatial distance information via clustering.
A different unsupervised perspective performs the change detection on image combination or transformation. In Deng et al. (2008), the Principal Component Analysis (PCA) is used to extract the difference between two images suppressing correlated information and highlighting variance in multi-temporal data. The change is identified in the second component, while the first component is assumed to be the sum of the common information. In Gao et al. (2016), PCA is used as a convolutional filter to determine the representative neighbourhood features from each pixel and generate change matrices with less noise spots. Gabor wavelets and fuzzy c-means are utilized to select interested bi-temporal pixels that have a high probability of being changed or unchanged. Then, new image patches centred at interested pixels are generated and a PCANet model-a deep learning network with its convolution filter banks chosen from PCA filters-is trained using these patches. Finally, pixels of bi-temporal images are classified by the trained PCANet model.
In , an autoencoder architecture is trained on the pixelwise difference computed between the spectral data of two images. By assuming that the spectral differences compressed at the bottom encoder layer preserve the hidden information to separate changed pixels from the unchanged background, encoded differences are coupled to distance information and processed through clustering to separate changed pixels from the unchanged background. The data compression ability of the encoder layer of an autoencoder is also investigated in Kalinicheva et al. (2018). Each image is encoded to an equal-sized compressed feature representation and the output of the subtraction operation between the encoded images is analysed for the change detection. In Kalinicheva et al. (2019), a convolution autoencoder is trained on the patches of a time series of images. The reconstruction error of each patch is analysed to discriminate changed pixels from the background.
Supervised CD methods (Khanday, 2016) are based on the availability of ground change information (often acquired by human intervention) and use a classification framework, in which the ground truth is used to learn a classifier. The spectral, spatial or proper combination of this information is used to build a measure able to detect the change and aid in the classifier decision. ANN (Clifton, 2003;Helmy & El-Taweel, 2010) and Expectation Maximization (EM) (Ming et al., 2014) algorithms fall in this category. Although both ANN and EM are based on different concepts (basically, the former is based on nonlinear regression, the latter on a maximum likelihood with unknown parameters), they provide a binary decision like the classifiers. In Seydi and Hasanlou (2017), a supervised method is illustrated. It acquires a sample of change labels on a scene, in order to determine the optimal threshold for determining the remaining labels with a distance-based change detection approach.
In , changes are identified by using a trained classifier to directly classify data from multiple periods (i.e., multi-date classification or direct classification) and comparing multiple classification maps (i.e. post-classification comparison). In Planinšič and Gleich (2018), a logistic regression layer is trained to perform supervised fine-tuning and classification on the autoencoder denoised representation of image time series feature extracted within tunable Q discrete wavelet transform. Finally, the transfer learning-based structure has been recently investigated to alleviate the lack of training samples and optimize the training process in a semi-supervised scenario. Transfer learning uses training in one domain to enable better results in another domain and, specifically, the lower to midlevel features learned in the original domain can be transferred as useful features in the new domain performing a fine-tuning according to a few labelled samples (Kerner et al., 2019;Larabi et al., 2019) 3 Preliminary concepts A MS/HS sensor records reflected light in tens (MS)/hundreds (HS) of narrow frequencies covering the visible, near-infrared and shortwave infrared bands of a wavelength λ (also called spectrum). The spectrum is an M-dimensional feature vector (spectral feature vector), so that λ is spanned on M numeric spectral features (bands) λ 1 , λ 2 , . . . , and λ M .
Let X and Y be two co-registered MS/HS images-digital images of an observed Earth's scene, which are produced at different time points using the same MS/HS sensor mounted on aircraft or satellites. Note that if the searched changes concern the vegetation, images should be acquired in the same season to avoid the focus on changes related to different phenological stages of the vegetation. X is denoted as the primary image, while Y is denoted as the secondary image. Every image (see Fig. 1) is an hyper-cube of size U × V × M, which represents a collection of spectral vectors measured on an M-dimensional spectrum λ over a grid of U ×V pixels. Every pixel (u, v) is a region of around a few square meters of the Earth's surface, which is a function of the sensor spatial resolution. X(u, v) and Y(u, v) are one-dimensional real-valued spectrum sections of hyper-cubes X and Y, respectively, indexed by spatial coordinates u and v within the sensor resolution of the camera.
Every pixel of a scene for which bi-temporal MS/HS data are acquired can, in principle, be labelled according to an unknown binary target function, whose range is a finite set of two distinct labels, i.e. "changed" and "unchanged". According to this function, a change matrix C can be associated with the bi-temporal image couple X and Y. In particular, C is a two-dimensional set of U × V change values with every value C(u, v) representing the change label of the pixel indexed by the spatial coordinates u and v. A CD method takes as input images X and Y to learn C.
In this paper, an unsupervised CVA method based on autoencoder is proposed. An autoencoder is a deep learning ANN trained to attempt to copy its input to its output (Goodfellow et al., 2016). It can be viewed as being composed of two functions: an encoder f -mapping the input vector x to a hidden representation h via a deterministic mapping h = f (x), parameterized by θ f -and a decoder g-mapping back the resulting hidden representation h to a reconstructed vector in the input space x = g(h), parameterized by θ g . The functions g and f correspond to two different ANNs combined in a single one, whose parameters {θ f , θ g } are simultaneously learned by minimizing a loss function

The proposed methodology
In this section we describe ORCHESTRA-an unsupervised CVA method enhanced with autoencoder information. The method takes as input the bi-temporal images, X (primary image) and Y (secondary image), and learns a change matrix C. Figure 2 shows the block diagram of ORCHESTRA.
Initially, we train the autoencoder architecture g · f on the pixel spectral vectors acquired with the primary image X. Since the activation produced by the top layer in the decoder network g corresponds to a reconstructed feature vector in the same M-dimensional spectral input space of the autoencoder, we consider this output feature vector as new learned features of the spectrum λ. The CVA of X and Y is then completed in this new feature space. According to these premises, g · f is used to restore the pixel spectral vectors of both X and Y and build the image reconstructions X and Y so that, for each pixel (u, v), Some inherent remarks can be formulated on the reconstructions X and Y . As the autoencoder g · f is trained on the pixel spectral vectors of X, we expect that X well reconstructs X as it mainly performs a denoising transformation of X. We also expect that g · f well reconstructs the spectral vectors of Y associated with the unchanged pixels, while it poorly reconstructs the spectral vectors of Y associated with the changed pixels. As a consequence, the spectral vector reconstructions of unchanged pixels in Y should be more similar to the corresponding spectral vector reconstructions restored in X than reconstructions associated with changed pixels. This conjecture (that is experimentally verified in Section 6) inspires the idea of computing the distance between the proposed autoencoder transformation of the original images, in order to better disentangle the differences between changed and unchanged pixels.
Then we compute pixelwise the distance between X and Y by resorting to the algorithm SAM that is commonly used in CVA methods (e.g. Appice et al., , 2020Lopez-Fandino et al., 2018). As pointed out in Seydi and Hasanlou (2017), the computation of SAM is independent of the number of spectral bands and insensitive to sunlight. Let us consider pixel (u, v), SAM(u, v) measures the angle between the bi-temporal reconstructed spectral vectors associated with (u, v) in both X and Y . This angle is computed as follows: Subsequently, we perform the Otsu's algorithm to automatically determine the upper threshold θ otsu of SAM distances for separating pixels of the study scene into background ("unchanged" pixels with low SAM range) and foreground ("changed" pixels with high SAM range). In particular, we assign pixels (u, v) with SAM(u, v) higher than θ otsu to the label "changed", while we assign the remaining pixels to the label "unchanged". The Otsu's algorithm is an adaptive, non-parametric and unsupervised threshold algorithm introduced in Otsu (1972). It is commonly used in image binarization problems to turn a single intensity threshold that separates pixels into two classes. The threshold is determined by minimising the intra-class intensity variance defined as a weighted sum of variances of the two classes. 1 In this paper, we assume that the SAM distances, computed pixelwise in the study scene, are represented in an histogram with L equal-width bins (levels) denoted as [1, . . . , L]. Let η i be the number of pixels at level i, so that L i=1 η i corresponds to the total number of pixels in the scene, i.e. L i=1 η i = UV . Based upon these premises, the probability of each level i is computed as p i = η i UV . The Otsu's algorithm identifies the optimal threshold level θ otsu , in order to divide the pixels of the processed scene into the background class C 1 , spanned over the SAM levels [1, 2, . . . , θ otsu ], and the foreground class C 2 , spanned over the SAM levels [θ otsu + 1, . . . , L], respectively. The optimal θ otsu is searched for minimizing the intra-class variance that is defined as a weighted sum of variances of the two classes: where σ 2 1 (θ) ad σ 2 2 (θ) are the variance computed on the two classes separated by θ. The weights w 1 (θ) and w 2 (θ) are the probabilities of the two classes, which are computed as follows: Further considerations concern the fact that the direct application of the Otsu's algorithm for change labelling will neglect the spatial arrangement of pixels. It may occasionally yield spurious assignments of pixels to classes. To avoid this issue, we may apply the principle of local auto-correlation congruence of objects (Appice et al., 2016(Appice et al., , 2017Du et al., 2012;, according to which detected clusters, comprising changed objects, generally expand across contiguous areas . Based on this principle, we may decide to change the assignment of pixels that strongly disagree with surrounding assignments. This mainly corresponds to performing a spatial-aware correction of the change assignment defined with Otsu's threshold. This correction assigns each pixel to the label that originally groups the majority of its neighbouring pixels (see Fig. 3). Formally, where c (u, v) and u(u, v) count how many pixels, falling in neighbourhood (u, v), are labelled as "changed" and "unchanged", respectively, with the Otsu's threshold. The neighbourhood (u, v) is a set of pixels surrounding (u, v) in the study scene. As in Guccione et al., 2015), we consider a Fig. 3 Majella Park bi-temporal scene: the separation of the scene into black changed pixels and white unchanged background, which is computed via the Otsu's algorithm (Fig. 3a); the spatial correction of the change labels assigned through clustering (Fig. 3b) square-shaped neighbourhood. Let R be a positive, integer-valued radius, the square-shaped neighbourhood (u, v) of pixel u, v is defined as follows:

Implementation details
ORCHESTRA is implemented in Python 3.8. A pre-processing step is performed to scale spectral data in the range [0, 1] and process spectral bands with values in comparable ranges. The autoencoder is developed in Keras 2.4.3 2 with TensorFlow 3 as the back-end. The set-up of learning rate and batch size is decided by resorting to the tree-structured Parzen estimator algorithm, as implemented in the Hyperopt library (Bergstra et al., 2013). This hyper-parameter optimization is done by using 20% of the entire training collection as a validation set. Therefore, we automatically choose the configuration of learning rate and batch size, which achieves the best validation loss in training the autoencoder. The values of learning rate and batch size explored with the tree-structured Parzen estimator, are defined as follows: learning rate varies in the range [0.00001, 0.01] and batch size ranges among 32, 64, 128, 256 and 512. The autoencoder architecture comprises 5 fully-connected (FC) layers of 128 × 64 × 32 × 64 × 128 neurons when trained with HS data and 3 fully-connected (FC) layers of 8×4×8 neurons when trained with MS data. Both architectures comprise a dropout layer to prevent overfitting. The mean squared error (mse) is used as the loss function. The classical rectified linear unit (ReLu) (Glorot et al., 2011) is selected as the activation function for each hidden layer, while for the last layer the Linear activation function is used. The number of epochs is set equal to 150, retaining the best models achieving the lowest loss on the validation set.
For the autoencoder architectures, both the number of layers and the number of neurons per layer are selected by taking into account the size of the spectral feature vector of each imagery dataset. In particular, the HS images are spanned on a spectral feature vector with either 224 or 242 spectral bands (see Table 1), while the MS images are spanned on a spectral feature vector with 13 spectral bands. As the MS data are simpler than the HS data, the autoencoder architecture adopted to process the MS images is simpler than the architecture to process the HS data images. On the other hand, we also account for the principle that a high number of layers may cause an increase of the computational effort that may not be rewarded with a gain in accuracy (Uzair & Jamil, 2020). With regard to the number of neurons, we follow the guidelines reported in Vanhoucke et al. (2011) and select the number of neurons in each hidden layer as a power of two, in order to improve the speed in the computation of the neural network. In fact, most of the computation time spent training an ANN is devoted to performing matrix multiplication. This is computed as a SIMD (single instruction, multiple data) operation in CPUs by using a batch size that is a power of 2. Finally, the threshold-based step is performed using the implementation of Otsu's algorithm from skimage.filters.threshold otsu, 4 with the number of levels L = 256.

Experimental evaluation and discussion
In this study we consider three co-registered, bi-temporal HS datasets (see Section 6.1) acquired in both rural and urban environments. For these datasets, the ground truth change information is available to validate the accuracy of ORCHESTRA. In particular, the accuracy performance is evaluated with the Overall Accuracy (OA), the number of Missed Alarms (MA -changed pixels assigned to the unchanged background) and the number of False Alarms (FA -unchanged pixels labelled as changed). These metrics are commonly considered in remote sensing for the evaluation of change detection methods. In addition, we measure the residual error of autoencoders (mean squared error on restored HS data) on both the primary image and the secondary image to explore the ability of the autoencoder in HS data reconstruction. The results, achieved on each dataset, are discussed in Section 6.2.

HS data
We consider three public available datasets 5 -Hermiston, Santa Barbara and Bay Area. Each data set comprises a couple of co-registered HS images of a scene, as well as ground truth information of the change occurred in the sensed scene. A brief description of the datasets is reported in Table 1.
In Hermiston, the study areas cover an irrigated agricultural field. This area provides a benchmark agricultural scene, which has been frequently used in the evaluation of the accuracy of HS CD methods (e.g. Appice et al., , 2020Lopez-Fandino et al., 2017. The land-cover types are soil, irrigated fields, rivers, buildings and types of cultivated land and grassland. In this dataset, the bi-temporal HS images were acquired with the HYPER-ION sensor. This is a space-borne system carried on the EO-1 satellite, which includes 242 spectral bands, covering wavelengths between 400 nm and 2.5 μm. The spectral range is divided into two intervals: the VNIR range (that includes 70 bands with wavelengths ranging from 356 to 1058 nm) and the SWIR range (that consists of 172 bands with wavelengths between 852 and 2577 nm). The spectral and spatial resolution of this sensor is about 10 nm and 30 m, respectively, over a 7.5-km strip. The Hermiston scene was monitored in the years 2004 and 2007 with the sensor over Hermiston City, Umatilla County, Oregon, USA. Each HS image of the dataset consists of 390 × 200 pixels acquired across 242 spectral bands.
In both Santa Barbara and Bay Area, the study areas cover an urban suburb in California. Both datasets have been already used in the evaluation of the HS CD methods illustrated in . In both the datasets, images were acquired by using the AVIRIS sensor. This is an optical sensor that delivers calibrated images of the upwelling spectral radiance in 224 contiguous spectral bands with wavelengths from 400 to 2500 nm. The spectral and spatial resolution of this sensor are about 10 nm and 4 m, respectively. The Santa Barbara scene was monitored in the years 2013 and 2014 with the sensor over the Santa Barbara region (California). It consists of 984 × 740 pixels and includes 224 spectral bands. The Bay Area scene was monitored in the years 2013 and 2015 with the sensor surrounding the city of Patterson (California). It consists of 600 × 500 pixels and includes 224 spectral bands.

Results
We start evaluating how the autoencoder trained on the primary image discloses knowledge that may contribute to separate changed pixels from unchanged pixels. To this aim, we explore how the autoencoder g ·f trained on the image of the couple assigned to the primary role can accurately reconstruct the primary image, while badly reconstructing the changed pixels of the secondary image. We evaluate two configurations defined by assigning the role of the primary image to (1) the oldest image and (2) the newest image of the couple, respectively. Table 2 reports the mean squared error (mse) computed comparing pixelwise each image to its reconstruction restored through the trained autoencoder. In both configurations, the autoencoder trained on the primary image reconstructs worse the secondary image, getting a poor restore of spectral vectors of changed pixels. This can be seen in Fig. 4a, b and c that depict the maps of the squared errors computed pixelwise on the reconstructions of the images acquired in the Hermiston dataset. The reconstructions are done with the autoencoder configuration (1) that is trained considering the oldest image acquired in 2004 as the primary image. These maps highlight that the changed area is already delineated from the poorly reconstructed pixels in the secondary image acquired on 2007. This supports our  hypothesis that the autoencoder transformation can disclose a representation of the spectral data, which contributes to better disentangle the change.
We proceed with measuring how the autoencoder can actually improve the accuracy of the CVA strategy. Table 3 reports the accuracy metrics of both ORCHESTRA and its baseline (CVA), that is defined by implementing the basic CVA with SAM and Otsu's algorithm on the original data (i.e. without the autoencoder architecture). Results show that both configurations of ORCHESTRA-(1) and (2)-outperform CVA. Interestingly, the highest accuracy (OA) is always achieved with the configuration of ORCHESTRA that maximizes the ratio of the mse computed on the reconstruction of the secondary image on the mse computed on the reconstruction of the primary image ( mse (Y,Y ) mse(X,X ) ) reported in Table 2). This defines a promising criterion to automatically select the best configuration of ORCHESTRA in an unsupervised manner. Final considerations concern the spatial correction that is beneficial except for Hermiston. (1) denotes the configuration of ORCHESTRA with the oldest image as primary image.
(2) denotes the configuration of ORCHESTRA with the newest image as primary image. (*) marks the configuration of ORCHESTRA that maximises mse(Y,Y ) mse(X,X ) in Table 2. Results are reported without spatial correction, as well as with spatial correction (R = 3). The highest OA is in bold Finally, we analyse the accuracy of few CVA methods that have been defined in the recent literature and evaluated on the same datasets. Table 4 reports the OA results. The compared methods also use SAM and spatial information for final label assignment. In addition, Lopez-Fandino et al. (2018) and López-Fandiño et al. (2019) 6 introduce the watershed analysis,  resorts to autoencoder for dimensional reduction, while  uses an iterative combination of clustering and classification. ORCHESTRA performs closely to competitors on Hermiston. It outperforms and López-Fandiño et al., 2019 on both Santa Barbara and Bay Area. On the other hand, the iterative procedure defined in Appice et al. (2020) may be considered for a future upgrade of ORCHESTRA.
Upon the completion of this comparative analysis, we perform the Friedman-Nemenyi statistical test (Demšar, 2006) on Hermiston, Santa Barbara and Bay Area. This test ranks the compared CVA methods for each dataset separately, so the best performing method is given rank of 1, the second best rank 2 and so on (Demšar, 2006). Figure 5 ranks the CVA methods according to the result of the Friedman-Nemenyi statistical test done on OA. The results of the test confirm that ORCHESTRA enables the construction of the change matrix that achieves the highest OA by having  as runner-up.

Majella national park analysis
Wildfires generate significant and complex environmental changes such as physical and chemical variations of soils, structural changes of vegetation, changes in ecological processes and ecosystem services (Meng & Zhao, 2017). Satellite MS data are traditionally exploited for monitoring burnt areas and wildfire effects. In this paper, we analyse the ability of ORCHESTRA in detecting environmental changes (e.g. physical and chemical variations of soils, structural changes of vegetation, changes in ecological processes and ecosystem services) caused by wildfires in MS images. In particular, we process two co-registered Sentinel2 L1C images acquired on both August 16, 2017 (Fig. 6a) and September 15, 2017 (Fig. 6b), in the area of the Morrone Mountain (within the Majella National Park, Italy). This area was burnt in a wildfire started on August 19, 2017 and lasted 25 days, which burnt more than 2,000 ha of an inaccessible area covered by coniferous forest and gorse. The processed MS images are composed of 1494 × 1338 pixels, with pixel resolution equal to 10m/pixels and MS resolution equal to 13 spectral bands (Aiello et al., 2019).

Fig. 5
Comparison based on the Friedman-Nemenyi test of OA computed on the change matrices built using ORCHESTRA and the related methods (Lopez-Fandino et al., 2018;López-Fandiño et al., 2019;Appice et al., , 2020. The authors of (Lopez-Fandino et al., 2018) and (López-Fandiño et al., 2019) test the same CVA approach We perform a preliminary analysis calculating the Normalized Burn Ratio (NBR) index on both the pre-fire and post-fire MS images. This index is commonly used to highlight burnt areas. Formally, where the reflectance in the mid-infrared band (SW I R), that is sensitive to the water content of both soil and vegetation, increases after a fire. On the other hand, the near-infrared band (NI R) declines in reflectance after a fire due to the decrease of the phytomass chlorophyllcontent. So, following the the conclusions drawn in Key and Benson (2006), we are able to assess the fire severity in a study area by measuring the difference between the NBR index calculated on both the pre-fire and a post-fire satellite images: In fact, this difference is correlated with the magnitude of changes caused by fires on the vegetation (Key & Benson, 2006). By assuming that the unburnt areas have similar spectral behaviour in two satellite images acquired before and after a fire, dNBR measures values around zero in unburnt areas, while it measures positive values in burnt areas. Figure 7 delineates the fire borders (red line) detected in the study area with the dNBR analysis conducted as described in Key and Benson (2006).
Although dNBR is one of the well-performing indexes in the detection of burnt areas over large fire zones with open forests and woodlands (Tran et al., 2018), it suffers from Burnt area detected by dNBR (red line) and ORCHESTRA (yellow zone), respectively. New burnt areas detected (blu circles), the false positives areas (orange circle) and the burnt areas detected independently of clouds (green circle) by ORCHESTRA a few limits. It is influenced by the fact that unburnt areas do not remain static over time, but they naturally undergo changes, passing from more or less dry/humid conditions over time. The parameters of the dNBR analysis need to be reviewed in each scenario based on several factors, e.g. seasonality of images, closeness of image acquisition to the fire event.
In addition, the dNBR computation is sensitive to variations in the soil brightness (Epting et al., 2005a), the type of vegetation (Epting et al., 2005b) and the density of the vegetation (Lentile et al., 2009). Finally, both clouds and their shadows can worsen the scenario when the dNBR analysis is done on large areas.
In this study, we explore which limits of dNBR analysis may be overcome by performing the CVA strategy with ORCHESTRA. To this aim, we consider the configuration of ORCHESTRA that handles the pre-fire image as X and the post-fire image as Y. This configuration allows us to train the autoencoder that maximizes the ratio of the mse computed on the reconstructed images. We apply the correction with R = 10. In particular, we focus the attention on: (1) the correctness of the detected fire borders; (2) the ability to detect any new burnt area correctly detected, as well as the presence of false-alarm areas; (3) the robustness of the performance to possible clouds. To reduce computational effort, we use the Corine Land Cover 2018 classification and then we analyse only pixels belonging to "Forest and semi-natural areas". Figure 7 highlights the advantages achieved with the CVA completed with ORCHESTRA. Blu circles underline that ORCHESTRA is able to detect newly burnt areas that are undetected with the dNBR. The green circle is a zoom-in to show the capability of ORCHESTRA to avoid changes that are due to the presence of clouds. Finally, we note that only one polygon (orange circle) is detected as a false alarm. We can conclude that, also for this particular dataset, ORCHESTRA shows good potential to reach more effective identification of burnt areas.

Conclusion
This paper describes a CVA method for analyzing a couple of optical satellite images (i.e., primary image and secondary image) acquired over time on the same scene, in order to separate pixels where a change occurs from the unchanged background in the scene. In particular, the proposed method takes advantage of autoencoders to identify spectral patterns that may aid in better disentangling changed pixels from unchanged ones. First an autoencoder is trained on the primary image and used to restore both the primary and the secondary image. Then the SAM distance is computed pixel-by-pixel between the restored images as a measure of the spectral change. Finally, the Otsu's algorithm is used on the computed distances to isolate the changed pixels, which are the pixels that measure the highest distance.
The novelty of the proposed CVA method is the specific use of an autoeconder architecture to transform the spectral data to compare, in order to enhance the spectral changes resulting in processed data. This is different from the common use of the autoencoders for data dimensionality reduction. Specifically, we base on the considerations that the autoencoder trained on the primary image should restore both the pixels of the primary image and the unchanged pixels of the secondary image accurately. Instead, it should see changed pixels of the secondary image as anomalies and reconstruct them badly. Therefore, computing a distance between restored spectral data measured at the same pixel aids in better delineating possible changes in the scene.
The experiments are performed by processing three couples of satellite HS images, collected either in a benchmark agricultural scene or in an urban scene. These experiments prove that the autoencoder component of the methodology contributes to the gain in detection accuracy. These experiments also reveal that the proposed method is able to provide competitive accuracy, compared to recent state-of-the-art CVA methods (comprising recent methods with autoencoders). In fact, with the encouraging performance of the proposed method, precise land-use and land-cover (or cropping pattern) changes may be identified. In addition, the method supplies promising results in the analysis of a couple of satellite MS images of a burnt area in the Majella National Park (Italy).
Some directions for further work are still to be explored. For example, appropriate classification algorithms may be studied to discriminate among different change types. The performances of various distance measures may be considered for the CVA. In addition, we plan to study the performance of the autoencoder-enhanced distance measures within a deep metric learning framework (e.g. Siamese network or Triplet network). Finally, we intend to investigate different autoencoder architectures, e.g. convolutional autoencoders, in the spectral data reconstruction.