A Method to Automatically Detect Changes in Multitemporal Spectral Indices: Application to Natural Disaster Damage Assessment

: This paper presents a new method, based on clustering and thresholding, to automatically perform binary change detection in multitemporal spectral indices. The method is denoted as Bu ﬀ er-From-Cluster Approach (BFCA). To estimate the distributions of changed and unchanged pixels, as needed for the purpose of a reliable thresholding of a spectral index, a clustering algorithm is preliminarily applied to identify image objects possibly corresponding to areas where signiﬁcant changes occurred. Then, a bu ﬀ er zone is created around the selected cluster to identify unchanged areas surrounding changed ones. The cluster and the bu ﬀ er zone are jointly analyzed to estimate the distributions of changed and unchanged pixels and to verify that they can be distinguished from each other. Finally, the results of thresholding and clustering are combined to generate the binary change map. The BFCA has been conceived to map the extent of the areas a ﬀ ected by a natural disaster like wildﬁre. To validate the proposed method, burned area maps produced by applying the BFCA to spectral indices derived from Sentinel-2 data have been compared to maps produced by the Copernicus Emergency Management Service. For testing the multi-hazard detection capability, the same kind of exercise has been carried out for a ﬂooding test case too. The positive results of the comparison have conﬁrmed the e ﬀ ectiveness of the proposed method.


Introduction
The frequency and the destructivity of natural hazards is continuously increasing. They leave in their wake a trail of injury, death, loss of livestock, property damage, and economic loss. This implies the need for improving the capability to manage these events. The availability of a synoptic map of the affected area in a short time strongly supports a rapid response to natural disasters [1]. Satellite earth observation (EO) data represent the most useful source of information to quickly map the extension of the area affected by a natural hazard [2].
Operational services like the rapid mapping component of the Copernicus Emergency Management Service (CEMS) currently provide emergency managers with maps of natural disasters at different spatial resolutions depending on the kind of EO data used to produce the maps. CEMS is an on-demand service, being triggered on request by authorized users (e.g., National Civil Protection Agencies), which is often at least partially based on the visual interpretation of a skilled operator. However, the present availability of Sentinel-1/2/3 data allows for routinely producing maps of areas affected by disasters (without the need of waiting for an activation). An unsupervised, fully automatic algorithm able to rapidly map the disaster's extent is needed to accomplish this routine production, especially if many EO data must be sequentially processed (e.g., when working at large spatial scales).
"Buffer-From-Cluster Approach" (BFCA). The BFCA has been conceived to work with one, or more than one, multitemporal spectral index (difference between the corresponding pixels in two spectral index images). The use of more than one spectral index generally adds complementary information for the change detection task (e.g., [28]), thus likely improving the accuracy of the result of the BFCA.
The BFCA combines clustering, thresholding, and region growing. With respect to the approach proposed in [27], that has been specifically designed for BA mapping, the BFCA can be applied to map the extent of other kinds of disasters like floods. In addition, the BFCA does not need any information furnished by ancillary products and the thresholds discriminating changed and unchanged pixels are determined using one of the automatic thresholding algorithms available in the literature [14]. With respect to the processor designed in [24] for AVHRR data, which uses the thresholding method developed in [14] and the RGA, the BFCA estimates the distribution functions representing changed and unchanged pixels to reliably apply automatic thresholding. In the processing chain designed in [22], which combines clustering, thresholding and region growing, a SBA similar to that proposed in [4,15] is used to determine the threshold. Conversely, the BFCA directly analyzes portions of an image of a multitemporal spectral index (MTSI) that likely contains comparable amounts of changed and unchanged pixels.
Burned area mapping has been chosen as main application of the BFCA and, to assess the aforementioned possibility to apply it even to other kinds of natural disasters, flood mapping has been tested as well. In both cases, Sentinel-2 (S2) data have been used to derive the MTSIs. This paper describes the various phases of the BFCA and presents the results of its application to S2-derived MTSIs useful to map BAs and floods. Validation has been accomplished by comparing the BFCA-derived maps to maps produced by the CEMS in the period 2018-2020, which have been therefore considered as the benchmark.

The Buffer-From-Cluster Approach
The flowchart of the new BFCA is shown in Figure 1. Each operation will be described in detail hereafter. Remote Sens. 2020, 12 The first step of the BFCA consists of the application of the Iterative Self-Organizing Data Analysis Technique Algorithm (ISODATA) to the input data. The ISODATA is an unsupervised clustering algorithm that separates groups of objects in a scene. A typical scene consists of regular and/or irregular regions arranged in a patchwork manner, each containing one class of surface cover

Clustering and Creation of the Buffer Zone
The first step of the BFCA consists of the application of the Iterative Self-Organizing Data Analysis Technique Algorithm (ISODATA) to the input data. The ISODATA is an unsupervised clustering algorithm that separates groups of objects in a scene. A typical scene consists of regular and/or irregular regions arranged in a patchwork manner, each containing one class of surface cover type (e.g., burned or flooded areas, croplands, forests, manmade structures, etc.). These homogeneous regions are the objects in the scene. Objects within each cluster should be as close to each other as possible and as far from other objects in other clusters as possible. The ISODATA begins with arbitrary cluster means and then clusters the pixels according to the minimum spectral distance technique. It is an iterative procedure, in which each iteration recalculates means and reclassifies pixels with respect to the new means. The process ends when a maximum number of iterations has been performed, or a maximum percentage of unchanged pixels between two iterations has been reached. The ISODATA does not keep a fixed number of clusters. In particular, clusters having a very large number of members are split, while if two cluster means are very close in terms of spectral distance (e.g., they represent an unnecessary or an injudicious division of the data), they are merged [29]. Clusters that contain so few points as to be meaningless (e.g., that they would not give acceptable statistics estimates if used in training a maximum likelihood classifier) are discarded and their members are reassigned to other clusters.
Since the ISODATA splits and merges clusters, it just requires specifying a range for the number of clusters. The optimum number of clusters is not known, so it is generally chosen to be conservatively high [29]. Hence, although here the goal is to distinguish changed pixels from unchanged ones, a maximum number of 10 clusters has been chosen as done in [22]. In Section 4.3, the impact of this choice on the BFCA-derived maps will be analyzed. A cluster can be considered as a sort of unknown class of surface cover type, and it is the task of the user to label the classes after the application of the clustering [29].
Let us assume that we are searching for positive values of the considered MTSI, that is, we are searching for significant increases of the spectral index. In this case, the objects belonging to the cluster presenting the highest median of the pixel values of the considered MTSI are selected and labelled as the class of changed pixels. Median value is used instead of the mean value being less sensitive to outliers. The reliability of the selection of the pixels belonging to the class of changed pixels can be improved by adding some constraints, e.g., by removing pixels having negative values of the MTSI (if present). If more than one MTSI is used as input data, the ISODATA clustering is sequentially applied to each index and the corresponding selected clusters are then intersected (i.e., only pixels included in overlapping clusters are maintained, thus basically applying a logical AND). The set of pixels derived in this way will be denoted hereafter as clustering-derived changed area.
The second step of the BFCA aims at searching for a reliable threshold separating changed and unchanged pixels. A buffer zone of buffering distance 50 pixels is initially created around the clustering-derived changed area. The buffer zone likely contains unchanged pixels. Hence, it is expected that, by merging the buffer zone with the clustering-derived changed area, the histogram of the pixel values of the considered MTSI will tend to be bimodal. If the clustering-derived changed area contains less than 30% of the pixels of the population formed by the union of the buffer zone and the clustering-derived changed area itself, the buffering distance is halved (25 pixels). The opposite (doubling) is done if the buffer zone contains less than 30% of the pixels of the total population. The processes of halving (or doubling) the buffering distance continues until at least 30% of the total population is included in each class. Note that, in [19,30], it is recommended that the smaller population should represent at least 10% of the other one. Hence, 30% is a conservative value chosen to be sure that the populations of changed and unchanged pixels are balanced. In Section 4.3, a discussion on the choice of 50 pixels as initial buffering distance will be provided.

Bimodality Check
Once the merging between the buffer zone and the clustering-derived changed area is carried out, the BFCA checks the bimodality of the histogram of the pixel values by evaluating the bimodality coefficient [31,32] and the Ashman's D coefficient [33].
The former, used also in [22,34], is defined as: where P is the population size, k is the kurtosis and γ is the skewness of the distribution. Values of BC greater than 5/9 (∼0.555) indicate a bimodal distribution.
The Ashman's D coefficient, used in [19], quantifies how well two Gaussian distributions are separated, thus assuming a Gaussian mixture model (GMM) for the distribution of changed and unchanged pixels. It is defined as [30]: where µ 1 and µ 2 are the mean values of the distributions and σ 1 and σ 2 are the corresponding standard deviations. A distribution has two peaks if D > 2. In order to carry out a reliable application of the D coefficient, a non-linear least squares fit to a Gaussian function is computed for the histograms. D is then determined based on the µ i and σ i (i = 1:2) parameters of the Gaussian functions fitting the histograms [19]. For a MTSI, the bimodality check is passed if BC > 5/9 and D > 2. Otherwise, the buffering distance is modified by halving or doubling it, as done in the previous step, depending on which population (clustering-derived changed area, or buffer zone) has the major number of samples. Minimum and maximum buffering distances of 3 and 150 pixels, respectively, have been assumed. When the buffering distance reaches the aforementioned limit values, the bimodality check stops.
If more than one MTSI is considered to perform the change detection through BFCA, it can be assumed that, overall, the bimodality check is passed if it is passed for the majority of the MTSIs. Otherwise, it can be deduced that changed pixels have not been found.

Thresholding and Region Growing
If, for a given MTSI, the bimodality check has been passed, the Otsu's automatic thresholding method [14] is applied to the corresponding histogram. If more than one MTSI is considered and the bimodality check is passed only for the majority of the MTSIs, empirical thresholds found in the literature can be used for the remaining indices (see Section 2.4).
Once a threshold (Th) has been determined, it is applied within the framework of a RGA, as done in [22][23][24]26,34], in order to account for the spatial context that is not considered if changed pixels are simply identified based on a threshold value, without any tolerance. To derive the seed region and the tolerance from the analysis of the distributions of changed and unchanged pixels, Th and µ 2 − 2σ 2 are used, where µ 2 and σ 2 are the mean value and the standard deviation of the Gaussian distribution fitting the histogram of the clustering-derived changed area. The seed region is represented by the pixels whose value is greater than the maximum between Th and µ 2 − 2σ 2 . The tolerance is the minimum between Th and µ 2 − 2σ 2 . Note that, as done in Section 2.1.1, here it is assumed that increases of the considered index are searched for. Otherwise, the seed region is represented by the pixels whose value is less than the minimum between Th and µ 2 + 2σ 2 , while the tolerance is the maximum between Th and µ 2 + 2σ 2 .
The application of thresholding and RGA generates the thresholding-derived changed area. If more than one MTSI is used, the corresponding thresholding-derived changed areas are intersected.

Final Classification
The final step of the BFCA consists of combining the clustering-derived and thresholding-derived changed areas. The combination is based on the following scheme: A.
Pixels belonging to both clustering-derived and thresholding-derived changed areas are classified as changed. B.
Pixels belonging only to the clustering-derived changed area are classified as changed if the corresponding cluster includes seed pixels. C.
Pixels belonging only to the thresholding-derived changed area are classified as changed if they are located in a buffer zone of buffering distance 50 pixels created around the changed area determined according to points A and B.

Data
The Sentinel-2 multispectral instrument measures the reflectance (ρ) in 13 spectral bands with spatial resolutions ranging from 10 to 60 m. For this study, Level-2A (L2A) S2 products, freely available from the Copernicus Open Access Hub have been used. L2A surface reflectance products are corrected for the atmospheric effects and include a Scene Classification (SCL) map, which is useful to mask clouds, snow, and water bodies. Each Level-2A product is available as a 100 × 100 km 2 tile in cartographic geometry (UTM/WGS84 projection). Only ρ data at 10 and 20 m resolution have been processed to generate the considered MTSIs. The 20-m grid has been chosen as reference, so that, before applying the BFCA, the 10 m S2 bands have been resampled to 20 m (5490 × 5490 pixels for each L2 tile). A land cover map derived from the most recent version (2018) of the CORINE land cover [35] has also been used to mask targets where BAs mapping is unreliable, as urban areas.
CEMS-derived maps of BA and flood extent have been used for validation purposes. They have been generated in the framework of the CEMS activations listed in Table 1 (flood mapping) and Table 2 (BA mapping). Note that since the use of multispectral data for flood mapping is often hampered by the presence of clouds, only one test case concerns flood mapping. Nonetheless, a case of a big flood mapped by CEMS by using also S2 data, namely the flood that hit Northern Queensland (Australia) in February 2019, has been found. The area affected by the flood was so wide that four S2 tiles have been processed (and the corresponding maps have been mosaicked) to generate a map to be compared to the CEMS-derived one. They are listed in Table 3. Most of the reference BA maps have been derived by CEMS by processing pre-fire and post fire S2 data. In these cases, the same images used by CEMS have been processed to test the BFCA. For the EMSR371, EMSR374, EMSR377 and EMSR443 case studies, CEMS produced the BA maps by analyzing Pleiades and SPOT-6/7 images as post-fire data. For the EMSR371 and EMSR377 case studies, S2 post-fire data acquired on the same days of the post-fire images used by CEMS have been selected to apply the BFCA. For the EMSR374 test case, the S2 image was acquired one day later than the SPOT one; for the EMSR443 test case, the S2 image was acquired two days later than the SPOT one. The S2 tiles used to generate the BFCA-derived BA maps are listed in Table 4.  Table 3. S2-L2A tiles used in this study for flood mapping.

CEMS Activation Code
Sentinel-2 Data

CEMS Activation Code
Sentinel-2 Data

Selection of Fire Spectral Indices for the Buffer-From-Cluster Approach
The three spectral indices chosen to test the BFCA for the BA mapping application have been selected by accounting for the results of the analysis of the separability parameter [36], similarly to what was done in [37,38]. The separability between two thematic classes (burned/changed pixels and unburned/unchanged pixels in this case) is defined as: where µ 1 and µ 2 are the mean values of the samples belonging to the classes and σ 1 and σ 2 are the corresponding standard deviations. By comparing Equations (2) and (3), it can be noted that the definition of S is quite similar to that of the Ashman's coefficient. The separability between two classes is generally considered good if S > 1 and the higher S, the higher the separability [36][37][38]. S has been evaluated for the differences between the pre-fire and post-fire values of the following indices: Normalized Difference Vegetation Index (NDVI) [39] and Normalized Burn Ratio (NBR) [40] used in [22,41], and Normalized Burned Ratio 2 (NBR2) [42] and Mid-Infrared Burned Index (MIRBI) [43], used in [23,38]. These indices are defined as: where NIR indicates the near infrared band, RED indicates the red band, and SWIR_L and SWIR_S are the short wave infrared long reflectance and short wave infrared short reflectance, respectively. The differences between the pre-fire and post-fire values of the aforementioned indices will be indicated hereafter as dNDVI, dNBR2, dNBR, and dMIRBI.
Besides dNBR, even its relativized forms proposed in the literature, namely the relativized dNBR (RdNBR) [44] and the relativized burn ratio (RBR) [45] have been considered in the separability analysis. These variables are defined as: To perform the separability analysis, the S2 data listed in Table 2 have been used. For each case study, the BA class has been derived by gathering the pixels labelled as burned by CEMS. The class of unburned pixels has been generated by considering both a buffer zone around the BAs and by randomly collecting pixels in areas labelled as unburned by CEMS. When building up the class of unburned pixels, attention has been paid to include in each class at least 30% of the total sample pixels and to exclude clouds, snow and water bodies. The result of the analysis is shown in Figure 2. Looking at Figure 2, it can be noted that the MTSIs that present the highest S between the classes of burned and unburned pixels are dNBR2 and dNBR. These MTSIs provide different information about a target, because NBR2 accounts for the SWIR bands while NBR considers SWIR_L and NIR bands. Conversely, RdNBR and RBR basically provide the same information as dNBR. Even dMIRBI presents a quite high value of S. Although MIRBI accounts for the SWIR bands like NBR2, it has been demonstrated in [38] that the distributions of burned and unburned classes are different for dNBR2 and dMIRBI. Hence, dNBR2, dNBR and dMIRBI have been selected as the MTSIs used to generate the maps of BA through the BFCA.

Application of the Buffer-From-Cluster Approach to Burned Area Mapping
The operations described in this section have been carried out for each test case. They have been implemented using the IDL programing language, and the ENVI routines that can be launched by means of specific IDL instructions. Firstly, the ISODATA clustering has been applied to the maps of dNBR2, dNBR and dMIRBI (excluding masked pixels). For burned pixels, the post-fire values of NBR and NBR2 tend to decrease with respect to the pre-fire ones [37,38,44,46], while the opposite occurs for MIRBI [23,38]. Hence, the clusters presenting the highest median value of dNBR and dNBR2 and the lowest median value of dMIRBI have been selected and then intersected, generating the clustering-derived changed (burned) area. Moreover, pixels presenting a negative value of dNBR or dNBR2, or a positive value of dMIRBI have been excluded from the clustering-derived changed area. To further increase the reliability of the clustering-derived changed area, even pixels whose NBR2postfire values were greater than the tile mean value of NBR2post-fire or whose MIRBIpost-fire values were less than the tile mean value of MIRBIpost-fire have been removed as suggested in [38].
Successively, for each MTSI, the bimodality check has been performed. It has been marked as passed if it has been passed for at least two out of three MTSIs. It has been verified that this happened for all the test cases listed in Table 2 and that, for most of the test cases, bimodal histograms were obtained for all the MTSIs. For bimodal histograms, the Otsu's method has been applied to determine the threshold Th, as described in Section 2.1.3. For the remaining non-bimodal histograms, the empirical thresholds proposed by in [38] for dNBR2 (0.05) and dMIRBI (−0.25) and in [46] for dNBR (260) have been used.
After having applied the RGA, the thresholding-derived BA has been combined with the clustering-derived one as described in Section 2.1.4 to produce, for each test case, a map of BAs. All Looking at Figure 2, it can be noted that the MTSIs that present the highest S between the classes of burned and unburned pixels are dNBR2 and dNBR. These MTSIs provide different information about a target, because NBR2 accounts for the SWIR bands while NBR considers SWIR_L and NIR bands. Conversely, RdNBR and RBR basically provide the same information as dNBR. Even dMIRBI presents a quite high value of S. Although MIRBI accounts for the SWIR bands like NBR2, it has been demonstrated in [38] that the distributions of burned and unburned classes are different for dNBR2 and dMIRBI. Hence, dNBR2, dNBR and dMIRBI have been selected as the MTSIs used to generate the maps of BA through the BFCA.

Application of the Buffer-From-Cluster Approach to Burned Area Mapping
The operations described in this section have been carried out for each test case. They have been implemented using the IDL programing language, and the ENVI routines that can be launched by means of specific IDL instructions. Firstly, the ISODATA clustering has been applied to the maps of dNBR2, dNBR and dMIRBI (excluding masked pixels). For burned pixels, the post-fire values of NBR and NBR2 tend to decrease with respect to the pre-fire ones [37,38,44,46], while the opposite occurs for MIRBI [23,38]. Hence, the clusters presenting the highest median value of dNBR and dNBR2 and the lowest median value of dMIRBI have been selected and then intersected, generating the clustering-derived changed (burned) area. Moreover, pixels presenting a negative value of dNBR or dNBR2, or a positive value of dMIRBI have been excluded from the clustering-derived changed area. To further increase the reliability of the clustering-derived changed area, even pixels whose NBR2 post-fire values were greater than the tile mean value of NBR2 post-fire or whose MIRBI post-fire values were less than the tile mean value of MIRBI post-fire have been removed as suggested in [38].
Successively, for each MTSI, the bimodality check has been performed. It has been marked as passed if it has been passed for at least two out of three MTSIs. It has been verified that this happened for all the test cases listed in Table 2 and that, for most of the test cases, bimodal histograms were obtained for all the MTSIs. For bimodal histograms, the Otsu's method has been applied to determine the threshold Th, as described in Section 2.1.3. For the remaining non-bimodal histograms, the empirical thresholds proposed by in [38] for dNBR2 (0.05) and dMIRBI (−0.25) and in [46] for dNBR (260) have been used.
After having applied the RGA, the thresholding-derived BA has been combined with the clustering-derived one as described in Section 2.1.4 to produce, for each test case, a map of BAs. All the burned pixels grouped in patches of ground smaller than 10,000 m 2 (25 pixels) have not been labelled as burned, thus assuming 1 ha as minimum mapping unit for BA mapping from S2.

Water Index and Application of the Buffer-From-Cluster Approach to Flood Mapping
Two water indices are commonly used to map floodwater using multispectral data. They are the Normalized Difference Water Index (NDWI) [47] and the Modified Normalized Difference Water Index (MNDWI) [48]. Positive values of NDWI and MNDWI correspond to the presence of water. In this case, we have decided to test the BFCA using only one MTSI. The MNDWI has been selected and increases of MNDWI are searched for by BFCA to map flooded areas. Hence, in this case dMNDWI is computed as the difference between post-flood and pre-flood values of the MNDWI. The latter is defined as: After having applied the ISODATA clustering, the clusters presenting the highest median value of dMNDWI have been selected to represent the clustering-derived changed (flooded) area. Pixels presenting a negative value of dMNDWI have been excluded from the clustering-derived changed area. The flood maps have been generated by combining the clustering-derived changed area with the thresholding-derived one.

Validation
To quantitatively represent the results of the comparison between BFCA-derived and CEMSderived maps, a confusion matrix has been computed for each case study, assuming the CEMS-derived maps as ground truth. A confusion matrix summarizes the performance of a classification algorithm. From the confusion matrices, four accuracy metrics have been derived: overall accuracy (OA), kappa coefficient (κ), errors of commission (ε comm ), and omission (ε om ) for the class of changed areas. It must be considered that, if the extension of the changed area is small, the class of changed pixels is underrepresented compared to the other class (unchanged pixels). This implies that very high values of OA are often obtained. To tackle this problem, limited geographic areas have been considered to derive the confusion matrices. The boundaries of the geographic areas have been determined using the "area of interest" shapefile included in the CEMS data (see [22] for more details). Figure 3 shows a full example of the results of the BFCA for the Calci test site (EMSR316, see Table 2). It includes the intermediate results of the various steps of the BFCA and the final BA map. Looking at Figure 3, it can be noted that, for each MTSI, the cluster presenting the highest median of the pixel values (red color in the left panels) has been labelled as the class of changed pixels. The upper central panel shows that with a buffer zone of buffering distance 50 pixels, the clustering-derived changed area contains less than 30% of the pixels of the population formed by the union of the buffer zone and the clustering-derived changed area itself. The same occurs with a buffering distance of 25 pixels. Conversely, with a buffer zone of buffering distance 12 pixels (blue color in the upper central panel of Figure 3) bimodal histograms have been obtained, as confirmed by Figure 4, where the histograms of the pixel values of the population generated by merging the buffer zone with the clustering-derived changed area are presented for dNBR2, dNBR, and dMIRBI. The Gaussian fitting functions are also shown. It can be seen that, in this case, the GMM is suitable to describe the distribution of the pixel       Looking at the upper right panel of Figure 3, it can be deduced that a very noisy map would have been obtained by using only the ISODATA clustering. Conversely, an underestimation of the BA would have resulted from the simple application of the thresholding and region growing techniques. This can be inferred by looking also at the lower central panel of Figure 3, which includes the boundary of the CEMS-derived BA too. Figure 5 shows, for the EMSR342 test case, the histograms of the pixel values of dMNDWI for the clustering-derived changed area and the buffer zone (S2 tile T54KVE). The latter histogram has been generated with a final value of the buffering distance equal to 25 pixels. A threshold value equal to 0.59 has been obtained. Even in this case the GMM properly describes the distribution of the pixel values. Figure 6 presents the result of the mosaicking of the flood maps produced by applying the BFCA to the four tiles listed in Table 3.

Results
Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 21 Figure 5 shows, for the EMSR342 test case, the histograms of the pixel values of dMNDWI for the clustering-derived changed area and the buffer zone (S2 tile T54KVE). The latter histogram has been generated with a final value of the buffering distance equal to 25 pixels. A threshold value equal to 0.59 has been obtained. Even in this case the GMM properly describes the distribution of the pixel values. Figure 6 presents the result of the mosaicking of the flood maps produced by applying the BFCA to the four tiles listed in Table 3.    Figure 5 shows, for the EMSR342 test case, the histograms of the pixel values of dMNDWI for the clustering-derived changed area and the buffer zone (S2 tile T54KVE). The latter histogram has been generated with a final value of the buffering distance equal to 25 pixels. A threshold value equal to 0.59 has been obtained. Even in this case the GMM properly describes the distribution of the pixel values. Figure 6 presents the result of the mosaicking of the flood maps produced by applying the BFCA to the four tiles listed in Table 3.   The results of the comparison between BFCA-derived and CEMS-derived maps, expressed in terms of the accuracy metrics introduced in Section 2.6 are synthetically presented in Figure 7 (OA and κ) and Figure 8 (ε comm and ε om ) for the BA mapping application.
The results of the comparison between BFCA-derived and CEMS-derived maps, expressed in terms of the accuracy metrics introduced in Section 2.6 are synthetically presented in Figure 7 (OA and κ) and Figure 8 (εcomm and εom) for the BA mapping application.  The OA always exceeds 91% with an average value (over the 13 test cases considered here) of 97.5%, while κ is in the range [0.80-0.97] with an average value of 0.88. With respect to the errors, εcomm ranges between 1.1% and 29.8% with an average value of 10.3%, while εom is in the range [0.2%-26.3%] with an average value equal to 9.5%. We have also calculated the errors for all the pixels in the 13 tiles together obtaining in this case εcomm = 4.3% and εom = 11.3%.

Discussion
This paper presents a new method to automatically perform binary change detection in MTSIs conceived to map the extent of the area affected by a natural disaster like wildfire. Several approaches dealing with unsupervised binary change detection in MTSIs are available in the literature (see the Introduction for some references). Most of them are generally quite complex using complicated data modeling and parameter estimation. The proposed method represents a rapid approach that aims at joining computational efficiency and good accuracy. The results of the comparison between BFCA-derived and CEMS-derived maps, expressed in terms of the accuracy metrics introduced in Section 2.6 are synthetically presented in Figure 7 (OA and κ) and Figure 8 (εcomm and εom) for the BA mapping application.  The OA always exceeds 91% with an average value (over the 13 test cases considered here) of 97.5%, while κ is in the range [0.80-0.97] with an average value of 0.88. With respect to the errors, εcomm ranges between 1.1% and 29.8% with an average value of 10.3%, while εom is in the range [0.2%-26.3%] with an average value equal to 9.5%. We have also calculated the errors for all the pixels in the 13 tiles together obtaining in this case εcomm = 4.3% and εom = 11.3%.

Discussion
This paper presents a new method to automatically perform binary change detection in MTSIs conceived to map the extent of the area affected by a natural disaster like wildfire. Several approaches dealing with unsupervised binary change detection in MTSIs are available in the literature (see the Introduction for some references). Most of them are generally quite complex using complicated data modeling and parameter estimation. The proposed method represents a rapid approach that aims at joining computational efficiency and good accuracy. The OA always exceeds 91% with an average value (over the 13 test cases considered here) of 97.5%, while κ is in the range [0.80-0.97] with an average value of 0.88. With respect to the errors, ε comm ranges between 1.1% and 29.8% with an average value of 10.3%, while ε om is in the range [0.2%-26.3%] with an average value equal to 9.5%. We have also calculated the errors for all the pixels in the 13 tiles together obtaining in this case ε comm = 4.3% and ε om = 11.3%.
As for the flood mapping application, the following values of the values of the accuracy metrics have been obtained: OA = 91.1%, κ = 0.80, ε comm = 12.1%, ε om = 0.5%.

Discussion
This paper presents a new method to automatically perform binary change detection in MTSIs conceived to map the extent of the area affected by a natural disaster like wildfire. Several approaches dealing with unsupervised binary change detection in MTSIs are available in the literature (see the Introduction for some references). Most of them are generally quite complex using complicated data modeling and parameter estimation. The proposed method represents a rapid approach that aims at joining computational efficiency and good accuracy.

Computational Efficiency and Accuracy
As for the computational efficiency, although, before carrying out the automatic thresholding, the BFCA uses ISODATA clustering, the latter is not significantly time consuming even if it is applied sequentially to three MTSIs, as done for the BA mapping application, mainly because just an iteration is used to cluster the pixels. Even the fitting of the histograms of changed and unchanged pixels with a Gaussian function is not time consuming (a maximum of 20 iteration has been fixed). It takes about 10 min to process a S2 image tile (5490 × 5490 pixels) using an i7-7800X 3.5GHz processor with 128 GB RAM.
As for the accuracy, looking at Figures 7 and 8, it can be deduced that a fairly good agreement between CEMS-derived and BFCA-derived BA maps has been obtained for the case studies considered here. In the literature, using S2 data and fully automatic BA mapping algorithms, omission errors equal to around 40% [49] and 26.5% [38] and commission errors of 19.3% [38] and 25% have been obtained [49]. Hence, average values around 10% for both errors can be considered as an indication that the BFCA works well in the identification of BAs from S2 data.
The case studies listed in the first five rows of Table 1 (EMSR306-368) were also analyzed in [22], where a SBA was used in combination with clustering. A significant reduction of ε om has been achieved through the BFCA, although this implied a slight increase of ε comm . A significant improvement of the results has been obtained in particular for the EMSR331 case study. The agreement with the CEMS-derived BA maps has increased also for the EMSR345 and EMSR 368 test cases, while for the EMSR306 and EMSR316 test cases, the decrease of ε om has been compensated by an increase of ε comm . It must be pointed out that the improvement of the performances of the BA mapping is due not only to the use of the BFCA, but also to a more accurate choice of the MTSIs, performed through the analysis of the separability between classes. To support this hypothesis, the algorithm proposed in [22] has been applied to the MTSIs used for BA mapping. For some test cases, e.g., EMSR316, values of the accuracy metrics very similar to those obtained using the BFCA have been achieved. For other test cases, a worsening of the performances with respect to those achieved through the BFCA has been observed. For instance, for the EMSR435 test case, ε om has increased from 23.1% (BFCA) to 24.2% and for the EMSR428 test case it has increased from 26.3% (BFCA) to 27.3%. The great importance of selecting suitable variables to detect changes in any change detection problem is pointed out in [3]. The case studies listed in Table 2 include fires having very different extents from about 1 km 2 (EMSR401 test case) to about 170 km 2 . The use of a medium-high resolution instrument like S2 contributed to the quite low values of ε om (e.g., 10.3% for the EMSR401 test case) achieved when fire burned few km 2 of surface. It can be expected that the BFCA might underestimate small-size fires using MTSIs derived from coarser resolution data. In [27,38,50], missed detection of fires having a size less than 1 km 2 has been ascribed to the use of data having a coarse spatial resolution (250-500 m).
The highest value of ε comm (29.8%) has been obtained for the EMSR377 case study. However, in this case the CEMS-derived maps have been produced using the data from a different satellite (SPOT-6/7). This may have an impact on the agreement between the BA maps especially if a small area (3.4 km 2 for EMSR377, see Table 2) is affected by fire, because in this case the fire likely had a quite low intensity and an ambiguous spectral signature.
The highest value of ε om (26.3%) has been obtained for the EMSR428 case study. Figure 9 shows an RGB color composite of the pre-fire (red) and post-fire (green and blue) dNBR2 data used to map the Tasarte fire (EMSR428). It can be noted that some areas classified as burned by CEMS appear in cyan tone; this indicates that NBR2 increased rather than decrease, as expected in a BA, thus explaining the large value of ε om . Note that even a Pleaides post-fire image acquired on the same day as the S2 one has been analyzed by CEMS to produce the BA map. This could at least partially explain the discrepancy, considering that Pleiades do not have SWIR bands. For what concerns flood mapping, it has been chosen to test the performances of the BFCA using just one MTSI, obtaining again a good agreement with reference data. The high OA (91.1%) is significant for the EMSR342 case study, because in this case the class of changed/flooded pixels is not underrepresented with respect to the other one. Using S2-derived dMNDWI data at a resolution of 20 m, commission errors ranging from 2% to 7.6% (lower than the εcomm value obtained through the BFCA for the EMSR342 case study) and omission errors in the range [9-42.1%] (considerably higher than the εom value obtained through the BFCA for the EMSR342 case study) have been reported in [51], but analyzing a very different study area (the Venice coastland) and only very small portions of the images (2 × 2 km 2 ). Even in this case, it is expected that combining more MTSIs could lead to the achievement of better accuracies, as discussed in [52], where quite high omission and commission errors (generally higher than those obtained through the BFCA for the EMSR342 case study) were reported by applying crisp thresholds to single MTSIs.
It must be pointed out that, while the effect of a wildfire can be assessed even several days after the event (e.g., reforestation generally takes a long time), the temporal scale of flood events can be very fast. For instance, the possibility to map flash floods vanishes few hours after the occurrence of the event. Even if flood mapping using multispectral instruments like S2 can provide good performances, microwave instruments, capable to operate in almost all-weather conditions and during both day-time and night-time, are generally more suitable to map floods also because of the sensitivity of the microwave radiation to water [53]. Hence, flood mapping is more commonly carried out using synthetic aperture radar (SAR) (e.g., [54,55]), which joins the aforementioned characteristics of microwave sensors to a high spatial resolution (that ranges from hundreds of meters to less than 1 m).

Combination of Different Image Processing Techniques
The combination of different image processing techniques is not new for what concerns BA mapping (e.g., [56]), flood mapping (e.g., [25,34]), and change detection in general (e.g., [7], where clustering was used in combination with principal components analysis). Here, clustering is basically used to estimate two distribution functions, one representing the class of changed pixels and the other one representing the background. Nonetheless, it has also a role in determining the final classification because objects belonging to the clustering-derived changed area are finally classified as changed if they include seeds of the RGA (see Section 2.1.4). Although, for BA mapping, the clustering-derived For what concerns flood mapping, it has been chosen to test the performances of the BFCA using just one MTSI, obtaining again a good agreement with reference data. The high OA (91.1%) is significant for the EMSR342 case study, because in this case the class of changed/flooded pixels is not underrepresented with respect to the other one. Using S2-derived dMNDWI data at a resolution of 20 m, commission errors ranging from 2% to 7.6% (lower than the ε comm value obtained through the BFCA for the EMSR342 case study) and omission errors in the range [9-42.1%] (considerably higher than the ε om value obtained through the BFCA for the EMSR342 case study) have been reported in [51], but analyzing a very different study area (the Venice coastland) and only very small portions of the images (2 × 2 km 2 ). Even in this case, it is expected that combining more MTSIs could lead to the achievement of better accuracies, as discussed in [52], where quite high omission and commission errors (generally higher than those obtained through the BFCA for the EMSR342 case study) were reported by applying crisp thresholds to single MTSIs.
It must be pointed out that, while the effect of a wildfire can be assessed even several days after the event (e.g., reforestation generally takes a long time), the temporal scale of flood events can be very fast. For instance, the possibility to map flash floods vanishes few hours after the occurrence of the event. Even if flood mapping using multispectral instruments like S2 can provide good performances, microwave instruments, capable to operate in almost all-weather conditions and during both day-time and night-time, are generally more suitable to map floods also because of the sensitivity of the microwave radiation to water [53]. Hence, flood mapping is more commonly carried out using synthetic aperture radar (SAR) (e.g., [54,55]), which joins the aforementioned characteristics of microwave sensors to a high spatial resolution (that ranges from hundreds of meters to less than 1 m).

Combination of Different Image Processing Techniques
The combination of different image processing techniques is not new for what concerns BA mapping (e.g., [56]), flood mapping (e.g., [25,34]), and change detection in general (e.g., [7], where clustering was used in combination with principal components analysis). Here, clustering is basically used to estimate two distribution functions, one representing the class of changed pixels and the other one representing the background. Nonetheless, it has also a role in determining the final classification because objects belonging to the clustering-derived changed area are finally classified as changed if they include seeds of the RGA (see Section 2.1.4). Although, for BA mapping, the clustering-derived changed area is generated by the sequential application of the ISODATA to different MTSIs, the intra-class variability in each MTSI can be high, so that ISODATA may not produce a very accurate classification. Looking at the upper right panel of Figure 3, it can be seen that ISODATA is generally able to identify the BA, but many false alarms may be raised. On the other hand, the application of the thresholding, even if performed in the framework of an RGA, is very sensitive to uncertainties in the determination of the seeds and the tolerances of the RGA. Since the BFCA fits the distributions of changed and unchanged pixels by assuming a GMM, an imperfect fitting can imply imperfect estimations of the thresholds and the tolerances.
To further underline the importance of combining different techniques, BA maps have also been produced by using a thresholding algorithm with thresholds derived from the SBA and applied within the framework of an RGA. The individual maps derived from the application of SBA, thresholding and RGA to each MTSI have been combined through a logical AND (i.e., intersected) and through a logical OR (i.e., merged). For the first case (map intersection) an increase of the omission error has been generally observed. For instance, for the EMSR316 test case, ε om has increased from 3.1% (BFCA) to 25.4% and for the EMSR428 test case it has exceeded 50%. As expected, for map merging, the opposite occurs. For instance, for the EMSR316 test case, ε comm has increased from 15.2% (BFCA) to 24.3%.
The values of the accuracy metrics shown in Figures 7 and 8 suggest that the combination of clustering and thresholding succeeds in compensating the errors related to the application of each technique alone. Figure 10 shows, for one MTSI (dNBR2), the ISODATA outputs obtained by selecting different values for the maximum number of classes: 5, 10, and 15. The EMSR316 test case has been used to carry out the exercise (hence the central panel of Figure 10 is the same as the upper left panel of Figure 3). The pixels belonging to the class having the highest median of the dNBR2 pixel values (red color in Figure 10) are almost the same in the three maps. For the BFCA, it is fundamental that the majority of changed pixels are identified by the ISODATA because it is basically used to estimate the distribution of this population. Hence, if among the classes of unchanged pixels discriminated by the ISODATA, which can for instance represent different types of land cover, oversegmentation or undersegmentation occur, this does not represent a problem for the accuracy of the BFCA. changed area is generated by the sequential application of the ISODATA to different MTSIs, the intraclass variability in each MTSI can be high, so that ISODATA may not produce a very accurate classification. Looking at the upper right panel of Figure 3, it can be seen that ISODATA is generally able to identify the BA, but many false alarms may be raised. On the other hand, the application of the thresholding, even if performed in the framework of an RGA, is very sensitive to uncertainties in the determination of the seeds and the tolerances of the RGA. Since the BFCA fits the distributions of changed and unchanged pixels by assuming a GMM, an imperfect fitting can imply imperfect estimations of the thresholds and the tolerances.

Maximum Number of Clusters and Initail Buffering Distance
To further underline the importance of combining different techniques, BA maps have also been produced by using a thresholding algorithm with thresholds derived from the SBA and applied within the framework of an RGA. The individual maps derived from the application of SBA, thresholding and RGA to each MTSI have been combined through a logical AND (i.e., intersected) and through a logical OR (i.e., merged). For the first case (map intersection) an increase of the omission error has been generally observed. For instance, for the EMSR316 test case, εom has increased from 3.1% (BFCA) to 25.4% and for the EMSR428 test case it has exceeded 50%. As expected, for map merging, the opposite occurs. For instance, for the EMSR316 test case, εcomm has increased from 15.2% (BFCA) to 24.3%.
The values of the accuracy metrics shown in Figures 7 and 8 suggest that the combination of clustering and thresholding succeeds in compensating the errors related to the application of each technique alone. Figure 10 shows, for one MTSI (dNBR2), the ISODATA outputs obtained by selecting different values for the maximum number of classes: 5, 10, and 15. The EMSR316 test case has been used to carry out the exercise (hence the central panel of Figure 10 is the same as the upper left panel of Figure   3). The pixels belonging to the class having the highest median of the dNBR2 pixel values (red color in Figure 10) are almost the same in the three maps. For the BFCA, it is fundamental that the majority of changed pixels are identified by the ISODATA because it is basically used to estimate the distribution of this population. Hence, if among the classes of unchanged pixels discriminated by the ISODATA, which can for instance represent different types of land cover, oversegmentation or undersegmentation occur, this does not represent a problem for the accuracy of the BFCA. Besides the selection of the maximum number of classes, another arbitrary choice done to implement the BFCA is the initial value of the buffering distance (50 pixels). Looking at the upper central panel of Figure 3, it can be noted that the initial buffer zone (grey color) covers a large part of Besides the selection of the maximum number of classes, another arbitrary choice done to implement the BFCA is the initial value of the buffering distance (50 pixels). Looking at the upper central panel of Figure 3, it can be noted that the initial buffer zone (grey color) covers a large part of the spatial subset around the area affected by the fire shown in the figure. This is due to the quite large number of small objects belonging to the class of changed pixels. In principle, for each object belonging to the class of changed pixels, a different dimension of the buffering distance, proportional to the size of the object, should be used. However, this would imply a significant decrease of the computational efficiency, especially if many objects belong to the class of interest. In the BFCA the buffer is mainly used to derive a bimodal histogram for the population formed by the selected changed and unchanged pixels. The adaptive modification of the buffering distance is an efficient way to mitigate this problem. Its increase/decrease ensures that two distinguishable modes are present in the histograms as confirmed by Figure 4 (and Figure 5 for flood mapping). Hence, it is not worthwhile to decrease the computational efficiency of the BFCA if this goal is achieved without generating different buffer sizes for objects of different sizes.

Maximum Number of Clusters and Initail Buffering Distance
In [27], unburned pixels around burned ones are searched for within a strip whose length has been empirically set between 10 and 20 km in order to determine a suitable threshold to identify BA. This strip is placed around a buffer having a distance of 1875 m. In [27], no adaptive modification of the strip and the buffer size is performed. It is worth noting that, for a pixel size of 20m, a buffering distance of 50 pixels corresponds to 1 km. Using 250-m resolution data, like the MODIS highest resolution ones from which BA maps are produced in [27], a range between 10 and 20 km corresponds to a buffering distance of 40-80 pixels. Hence, using the BFCA for coarse spatial resolution data (250-500 m), an initial buffering distance greater than that chosen for S2 (e.g., 80 pixels) can be used.

Conclusions and Perspectives
The new Buffer-From-Cluster Approach (BFCA) to automatically perform binary change detection in multitemporal spectral indices has been presented. It has been conceived to map the extent of areas affected by a natural hazard like wildfire. The BFCA combines different image processing techniques such as clustering and thresholding, which are applied within the framework of a region growing algorithm. Nonetheless, it is quite simple in terms of computation and effective in identifying changes caused by the occurrence of a disaster. It is therefore suitable for real-time applications, especially if a large amount of multispectral data has to be processed, for instance for routinely producing maps of areas affected by disasters.
The effectivity of the BFCA in the identification of changes caused by the occurrence of a disaster has been verified by comparing BFCA-derived burned area maps, generated using Sentinel-2 data, with maps produced by the CEMS. Values around 10% for omission and commission errors have been obtained considering the average over all the case studies. They never exceeded 30%, while the overall accuracy always exceeded 91%.
The values of the accuracy metrics seem to indicate that the BFCA provides results better than those found in the literature for automatic burned areas mapping using Sentinel-2 data and fully automatic algorithms. However, more tests are necessary to confirm this indication. For this purpose, it is planned to make available the BA mapping algorithm derived from the BFCA in the WASDI platform (www.wasdi.net), where registered users will be able to test it on their own Sentinel-2 data. In addition, starting from summer 2020, this algorithm will update a previous version implemented to provide, for the Italian territory, a daily service of burned areas mapping based on Sentinel-2 data [22].
Although burned areas mapping has been chosen as the main application of the BFCA, the latter has been conceived to be applied also to other disasters like floods. A test for this kind of application gave encouraging results (omission error of 0.5%, commission error of 12%). Even the algorithm for flood mapping using Sentinel-2 data based on the BFCA will be made available through the WASDI platform. However, it must be pointed out that the inability of multispectral instruments to provide data when cloud cover is present or during night-time strongly limits their effectivity for flood mapping. Therefore, future work will evaluate the possibility to adapt the BFCA to flood mapping from SAR, which is the most suitable sensor for this application.