Comparative Analyses of Unsupervised PCA K-Means Change Detection Algorithm from the Viewpoint of Follow-Up Plan

In this study, principal component analysis and k-means clustering (PCAKM) methods for synthetic aperture radar (SAR) data are analyzed to reduce the sensitivity caused by changes in the parameters and input images of the algorithm, increase the accuracy, and make an improvement in the computation time, which are advantageous for scoring in the follow-up plan. Although there are many supervised methods described in the literature, unsupervised methods may be more appropriate in terms of computing time, data scarcity, and explainability in order to supply a trustworthy system. We consider the PCAKM algorithm, which is used as a benchmark method in many studies when making comparisons. Error metrics, computing times, and utility functions are calculated for 22 modified PCAKM regarding difference images and filtering methods. Various images with different characteristics affect the results of the configurations. However, it is evident that the PCAKM becomes less sensitive and more accurate for both the overall results and image results. Scoring by utilizing these results and other map information is a gap and innovation. Obtaining a change map in a fast, explainable, more robust and less sensitive way is one of the aims of our studies on scoring points in the follow-up plan.


Introduction
Change detection for temporal differential images is the implementation of an the algorithm/method to detect the changes that have occurred between two images obtained at different times from the same sensor, platform, and location. In other words, it is a process to divide the map into regions that are changed and unchanged.
Change detection algorithms are used in several areas such as video surveillance, remote sensing, medical diagnosis and treatment, civil infrastructure, underwater sensing, and driver assistance systems [1]. Different types of systems in remote sensing and aerial photography are used to detect changes between the scenes of the same location acquired at different times, which is also called remote sensing change detection [2]. Depending on the sensors, systems used, and the time-frequency of the images obtained, different tasks are assigned and executed. Such changes can trigger follow-up activities to determine the cause or type of change, such as triggering additional image requests [3,4], direct actions such as search and rescue missions [5], or influencing decisions made in the area, e.g., threat avoidance [6]. In all of these cases, quick, precise, and interpretable change detection is critical to deriving timely information and properly reacting in the subsequent follow-up. A synthetic aperture radar (SAR) sensor is extensively used in numerous areas [7] to obtain change maps since it is not affected by the weather, light, or flight altitude [8,9]. These images typically form the function for actions such as those listed above. Speckle noise [10][11][12], fuzzy edge of changed regions [12], and limited datasets [12] are the main challenges for the change detection of SAR images.
Future changes that are likely to occur can be predicted using spatial and temporal dynamics. Therefore, the follow-up information acquisition can be streamlined, by creating the foundation for the follow-up plan by scoring the predicted change detection map [13]. Follow-up activities require the consideration of the estimated change map's accuracy and computation time. Most follow-up planning may not allow for the large amounts of data that are necessary for supervised learning. Due to this, the follow-up planning benefits the adoption of unsupervised techniques that do not need training data and have quick computation times. In addition, classical methods generally provide more transparency, explainability, and interpretability than more complex ones do [14]. These properties support the system's trustworthiness [15][16][17], which is significant for any response to the detected change.
In this study, the change maps of different satellite images are calculated using the unsupervised change detection algorithm proposed by Celik [18] and called principal component analysis and k-means clustering (PCAKM). In situations where the follow-up plan needs to be made in short periods (such as disaster response, etc.) and training data are lacking, it is more appropriate to use unsupervised methods instead of supervised methods. In addition, it is not necessary for unsupervised change detection methods to specifically identify the kinds of changes in land use or cover that have occurred in the area of interest [19]. Depending on the change in the parameters used in PCAKM, the calculation times and performance results of the obtained change detection alter. Moreover, altering the inputs for Celik's algorithm affects results notably.
We produce several configurations as modified PCAKMs using different filters and DIs. The performance results obtained based on modified versions of Celik's algorithm are compared and examined to understand whether they are suitable for generating scores to form the foundation for planning follow-up detailed investigations or responses. As a result of these investigations, we seek to answer the following questions: • Is it possible to decrease sensitivity or increase consistency? • Is it possible to decrease computing time without decreasing accuracy?
These questions are critical to obtaining a modified method that is less affected by different PCAKM algorithm parameters and input image characteristics. The way to achieve this is to increase the average performance and reduce the variance of the results obtained. On the other hand, providing a decrease in computing time is important in real-time tasks. Change maps with a less sensitive method will then be input into the scoring stage for a follow-up plan. Change maps with lower error rates and variance will help the selection of the areas of interest (AOIs) generating the points of these AOIs and scoring. A gap and innovation in the follow-up plan is scoring using change map results and other map information. We aim to contribute to scoring points in the follow-up plan by focusing on obtaining the change map in a quick, explainable, more accurate, and less sensitive manner.
The paper is organized as follows. Section 2 provides a related literature review and the methods used in this paper. Section 3 explains the data, configurations, and performance metrics used in experiments, and shows the results. Comments and discussions on the results are included in Section 4. Finally, the paper is concluded in Section 4.2.

Methods
To find answers to the questions in Section 1, we make comparative analyses by performing unsupervised change detection and performance evaluation. Performance measurements were carried out using the unsupervised learning method for the change detection part. The method includes principal component analysis (PCA) and k-means clustering methods as described in [18].
The PCAKM method proposed by Celik [18] has been used as a benchmark comparison in many unsupervised and supervised SAR image change detection studies and continues to be used in the state-of-the-art research. Li et al. [20] compared PCAKM, Markov random field fuzzy c-means (MRFFCM), Gabor fuzzy c-means (GaborFCM), and Gabor two-layer classifier (GaborTLC). They determined that the Kappa coefficient (KC) difference between PCAKM and other methods for the benchmark data is at most 2.67%.
Gao et al. [13] proposed deep semi-non-negative matrix factorization (NMF) and a singular value decomposition network to compare with PCAKM, MRFFCM, GaborTLC, and deep neural networks with MRFFCM (D_MRFFCM). KC differences between PCAKM and methods that give better results for the three benchmark data are less than or equal to 6.77%, 7.57%, and 3.3%, respectively. In addition, PCAKM has a higher KC than MRFFCM and GaborTLC for two out of three data. Gamma deep belief network (gΓ-DBN) is proposed by Jia and Zhao [12] for comparison with PCAKM, convolutional-wavelet neural network (CWNN), deep belief network (DBN), and joint deep belief network (JDBN). According to their experimental results, PCAKM shows better performance than CWNN and DBN for one of the benchmark datasets in terms of KC. When all images are examined, improvements in KC for PCAKM are less than or equal to 7.98%. Wang et al. [21] presented a graph-based knowledge supplement network (GKSNet) to match against PCAKM, a neighborhood-based ratio and extreme learning machine (NR-ELM), Gabor PCA network (GaborPCANet), local restricted convolutional neural network (LR-CNN), transferred multilevel fusion network (MLFN), DBN, and deep cascade network (DCNet). In their study, the KC/F1-measure enhancements for PCAKM are less than or equal to 10.71%/9.42%, 2.75%/2.2%, 19.21%/22.57%, and 11.43%/9.75% for four different benchmark data, respectively. Even though the supervised methods provide these improvements in terms of accuracy, their run-time results show that there are reasonably high differences between PCAKM and other methods. The average run-times in seconds for PCAKM, NR-ELM, GaborPCANet, LR-CNN, MLFN, DBN, DCNet, and GKSNet are 2.3, 22.5, 442.8, 282.6, 187.6, 474.1, 509.6, and 144.92, respectively [21].
Considering its features such as being fast, not requiring learning data, and having a simple algorithm, we selected PCAKM as a benchmark method. It shows promising results for both unsupervised [13,18,20] and supervised methods [12,21].
Speckle is a type of grainy noise that occurs naturally in active radar, SAR, medical ultrasound, and optical coherence tomography images, and decreases their quality. Images of the same region taken at different times have different levels of speckle. Speckling creates difficulty in distinguishing opposite classes [22] since it increases the overlap of opposite-class pixels in the histogram of difference images. On the other hand, there is competitive interaction between altered regions and background regions due to a lack of past information, resulting in a fuzzy edge in the changed region that is difficult to discern. Another challenge is the lack of data, which is an issue for supervised learning.
Noise may develop as a result of the system's construction, illumination conditions, and image acquisition process. Numerous methods were proposed for speckle reduction or despeckling. Speckle reduction filters are classified as non-adaptive and adaptive filters [23]. Mean and median filtering methods are examples of non-adaptive techniques. On the other hand, Lee, Frost, Kuan, and G-MAP are adaptive filter examples. Qiu et al. [24] claimed that none of these filters consistently outperform others, in principle. Each filter has particular advantages and disadvantages according to the data. For this reason, choosing a more stable and consistent filter is important.
Moreover, speckling reduction techniques are categorized as the spatial domain, transform domain (or wavelet domain), non-local filtering, and total variational [23]. Specifically, anisotropic diffusion, bilateral filter (BF), fast non-local means filter (FNLMF), and guided filter (GF) are some other filters to reduce speckle noise [25]. Choi and Jeong [25] state that BF and the non-local mean filter (NLMF) have a low speckle noise reduction performance. In addition, non-linear methods such as BF and NLMF have poor computational time performance [25]. Partial differential Equation (PDE)-based algorithms including AD and adaptive window anisotropic diffusion also have a weak performance on speckle noise removal [25]. Some other conventional filtering methods such as discrete wavelet transform (DWT), Bayesian multiscale method in a non-homomorphic framework, and expectation maximization DWT perform poorly in terms of speckle noise removal, edge information preservation, and computing complexity [25]. We test the edge-protecting GF method, which has low computational complexity among the speckle noise elimination techniques considering the performance results for the SAR images in [25]. We also used the NLMF and BF methods to acquire the performance characteristics mentioned above.
In this study, we compared the PCAKM with its modified versions (different combinations of difference images and filters) in terms of accuracy and time performance. We consider whether there is a modified method with less sensitivity and higher accuracy for the change map to be used in any follow-up plan.

Original PCA K-Means Algorithm
The flow of the proposed original method is given in Figure 1 [18]. PCA and k-means methods [18] are utilized for the change detection part.  Firstly, the input image pairs I 1 and I 2 are converted into grayscale images. Then, the absolute difference image for the given image pair is calculated as Afterward, D 1 is divided into bs × bs non-overlapping blocks where bs is the length of one side of square blocks. After converting these blocks into row vectors, PCA is applied to these vector sets to obtain the orthonormal eigenvectors. In the next step, the feature vector space is created by projecting bs × bs overlapping blocks around each pixel onto the eigenvector space. Feature vector space is input for the k-means algorithm to get the change map. Using the k-means algorithm, the feature vector space is grouped into clusters. Then, each pixel is assigned to a cluster in a way that minimizes the distance between its feature vector and the cluster's mean vector. Briefly, we used two parameters bs and k as block width and cluster number, respectively. In Section 3, bs is between 2 and 8, whereas k is 2 and 3 for each image pair.

Other Difference Image Methods
The log-ratio difference image method, which is given in Equation (2), is utilized in many studies to reduce the multiplicative distortion effects of noise caused by speckle [10]. Moreover, Zhao et al. [26] produced the difference image via image regression as given in Equation (3) to avoid problems such as atmospheric condition changes, illumination variations, and sensor calibration [27]. The image regression method enhances the performance of the difference image, which is observed from direct subtracting. However, both the log-ratio and absolute log-ratio methods still do not perform well enough to eliminate speckle noise if the input becomes low-quality [10,27].
Definition 1 (Log-ratio difference image). Log-ratio image is the logarithmic transform of the image pair's division as where f = c * log((1 + p), 10), c ≈ 105 for all pixels p ∈ I 1 I 2 .
Definition 2 (Absolute log-ratio difference image). It is the absolute value of the log-ratio calculation as Zhang et al. [10] stated that the SAR images are contaminated by speckle noise, which has the multiplicative Goodman's model. The Nakagami distribution in Equation (4) is then used to represent the independently and identically distributed pixel amplitudes. The Nakagami distribution is where R s and I s are the reflectivity and pixel amplitudes in site s, respectively. Moreover, L is the equivalent number of looks, which is a parameter of multi-look SAR images, and represents the amount of averaging done to the SAR measurements both during the creation of the data and, on occasion, even after [28]. After several calculations to which Bayesian decision theory was applied, the difference image is given as in Equation (5), where it considers the knowledge that the speckles follow the Nakagami distributions.
Definition 3 (Nakagami log-ratio (NLR) difference image). It is a modified version of the log-ratio difference image given as Its absolute version can be written as Definition 4 (Modified NLR difference image 1). In this version of the NLR difference image, we use the squared values of each image given as ).
Its absolute value is Definition 5 (Modified NLR difference image 2). For this modified version of the NLR difference image, squares of each division are added to the NLR difference image itself as ).
The absolute value of it is Definition 6 (Improved ratio and log improved ratio difference image [29,30]). The improved ratio and its log transform version are given in Equations (11) and (12), respectively.

Non-Local Means Denoising
Basically, the color of a pixel is changed to an average of the colors of nearby pixels by non-local means denoising (NLMD) [31]. Since there is no justification for the closest pixels to a given pixel to be even close, it searches across a sizable chunk of the image for every pixel that resembles the pixel to be denoised. There are three parameters such as h, templateWindowsSize (tws), and searchWindowsSize (sws). The first one regulates the filter strength. If it is increased, then it removes the noise more precisely but removes the image details as well and vice versa. The tws parameter is the template patch's size in pixels, which is utilized to calculate weights. Lastly, sws is the window's size in pixels that is applied to estimate the weighted average for a specific pixel. We used OpenCV's recommended values for the last two parameters as 7 and 21, respectively. On the other hand, for h, we used 20 since SAR images contain a high degree of noise.

Bilateral Filter
In addition to using a (multiplicative) Gaussian filter component that is based on pixel intensity differences, the bilateral filter (BF) also employs a Gaussian filter in the space domain. Only pixels that are "spatial neighbors" are taken into account for filtering, owing to the Gaussian function of space. On the other hand, the Gaussian component used in the intensity domain makes sure that only the pixels with intensities close to the core pixel are taken into account when computing the blurred intensity value. BF is a method that preserves the edge information. We used 10 for the parameter denoisingWindowsize (dws), which is larger than the default value 3, which is similar to NLMD, and we consider that the SAR image has substantial noise.

Guided Filter
A guided filter is a smoothing light filter that preserves the edges. It filters out noise or texture while keeping sharp edges, just like a bilateral filter [32,33]. The GF is defined by the following Equations (13)-(15) as where (a k , b k ) are linear coefficients for a linear transform of the guidance image I at a pixel i with the input image p and supposed to be constant in a window w k (square window of a radius r) centered at the pixel k. Furthermore, µ k and σ 2 k are the mean and variance of I in w k , |w| is the numbers of pixels in w k , p k = 1 |w| ∑ i∈w k p i is the mean of p in w k , and is a regularization parameter penalizing large a k . Moreover, a i = 1 |w| ∑ k∈w i a k and b i = 1 |w| ∑ k∈w i b k are the overall average coefficients for windows that overlap with i where q i is the filtering output at a pixel i.

Truncated Singular Value Decomposition
Truncated singular value decomposition (TSVD) is a reduced rank approximation to any matrix A by selecting the first major singular values. We determine the subset of full components via the percentage of the total variance. Therefore, we utilize the var parameter that is a threshold for total variance. The reason for applying this method is to assess whether we can reduce the time performance without much loss of accuracy.

Data
Details of data used in experimental results are given in Table 1. For each image pair, there is a ground truth image for the change between the two images. The ground truth images are used to generate the confusion matrices and calculate the performance metrics mentioned above.
In Table 2, noise variance values based on the method in [34] are given. In Figure 2, all images with histograms are demonstrated. In Figure 3, Radon transforms between 0 and 180 degrees are illustrated. Radon transforms, which are also called sinograms, calculate image matrix projections over predetermined directions where lighter tones are more intense. It is apparent that there are different characteristics not just among each image pair but also between some image pairs across the set, as shown in Table 2, Figures 2 and 3.

Configurations
There are seven SAR image pairs and seven ground truths for change maps as data. We utilized 22 different configurations among which one is the original paper [18] as a benchmark and the others are modified versions of the original method. For each configuration, there are 98 (7 × 14) change detection results for each performance metric since we have two main parameters block size and number of clusters, which take values in the ranges of 2-8 and 2-3, respectively. We calculate each change detection result 1000 times and then obtain the minimum, maximum, and average calculation times. The accuracy results do not change for these 1000 experiments since all 22 configurations have deterministic skeletons.
All configurations are given in Table 3 with configuration numbers. Configurations containing more than one method are written according to the order of their implementation. The PCAKM algorithm is used after applying the written methods for any configuration. We select the radius of the square window (r) for GF as the block size (bs) parameter of the PCAKM algorithm. Explanations for other parameters in Table 3 are given in Section 2.

Performance Metrics
After calculating the change maps, performance metrics are estimated using the confusion matrix that is given in Table 4. Below are formulations for performance metrics using the true positive (tp), false positive (fp), false negative (fn), and true negative (tn) in the confusion matrix as where n = tp + tn + f p + f n. We use Kappa coefficient and f-measure as accuracy calculations. The range of the former is [−1, 1] and the latter has a range of [0, 1]. For both metrics, a higher value means better accuracy.
On the other hand, we estimate the utility functions by employing Kappa coefficient, f-measure, and average computing times. For each image pair, we have two utility values as where µ ij1 is the average of Kappa coefficient values, µ ij1 is the average of f-measure values, is the variance of f-measure values,t ij is the mean of average computing times, i is the image pair number, and j is the configuration number for i = 1, . . . , 7 and j = 1, . . . , 22. For each configuration, we have 14 results since we utilize the parameters block size and number of clusters, which take values in the ranges 2-8 and 2-3, respectively. Then, we use these 14 values for mean and variance calculations. On the other hand, we have 14 different average time calculations and each parameter pair result is calculated 1000 times. Then, we take the average of these 1000 calculation times and estimate the mean of 14 average calculation times.
In addition to the U1 and U2 utility values, we calculate the following utility values for overall images in a single configuration as where µ k1 and µ k2 are the average Kappa coefficient value and the average f-measure value of all 98 results (14 parameters combination for seven images), σ 2 k1 and σ 2 k2 are the average variances of seven different image variance results for each configuration,t k is the mean of seven images' average time computations (each image has 14 different averaged time results for 1000 experiments), and k is the configuration number for k = 1, . . . , 22. Since, as we mentioned in Section 3.1, each image pair has different characteristics according to noise variances, histograms, and Radon transforms, we take the average of seven different image variance results for each configuration.

Results
The best and the worst results and the mean and variance for error metrics (kc and fmeas) of 22 configurations, are given in Tables A1-A7 for each image pair, respectively. The bs and c demonstrate the block size and the number of clusters parameters. Furthermore, the "No" columns in the tables state configuration numbers. We calculate the U1 and U2 utility values in these tables by utilizing each configuration's average computing time. Image pairs, ground truth images, and the best change map result for all images are given in Figure A1 in Appendix B. The first two columns contain the image pairs. The third and fourth columns are the ground truths and best change map results, respectively.
Based on the image results given in Appendix A, Table 5 shows the order of configurations from largest to smallest utility values. On the other hand, Table 6 presents the overall mean and variance of error metrics for each configuration with average computing times regarding the configuration numbers. We calculated U3 and U4, where the mean and variance values are the average values of seven different image results. Bold values show the highest mean and lowest variance estimations for error metrics.  Table 7 demonstrates the ranking of U3 and U4 values in terms of configuration numbers. If we look at Tables A1-A7, it is evident that the absolute versions of difference images do not increase the mean accuracy values for any configurations and images. There is no systematic decrease or increase for variance and average time performances as well.
We check the configuration pairs 2&3, 9&11, and 9&12 to observe the effects of TSVD. Applying TSVD decreases the average computing times for the Ottawa, Yellow River 2, Yellow River 4, San Francisco, and Bern datasets. It increases the mean values for Ottawa and the Yellow River 4, decreases the variance for Ottawa and Yellow River 2, and increases variance for Yellow River 3. Otherwise, increases and decreases show variability.
Furthermore, if we consider configuration pairs 2&4, 9&10, 15&16, and 18&19, utilizing GF increases the average computing times for all images. However, there is no regular increase and decrease path for mean and variance since difference images and images affect the procedure. In addition, utilizing TSVD and GF with D4 does not expose ordered changes for mean, variance, and time if we regard each image pair. Furthermore, using both NLMD and NLMD with GF increases the average computing times for D2 in a noticeable way due to NLMD. Nevertheless, changes in mean and variance across all images are not in harmony, which is the case also for BF for D2.
Applying different methods generally illustrates different mean, variance, and time change effects as difference images and images are altered. Similarly, the ranking of U1 and U2 values changes depending on the image since each image has different characteristics, as mentioned in Section 3.1. Therefore, at this point, it would be more accurate to focus on the utility values obtained by considering the overall results. Because even if we use the same sensor, the dataset to be obtained may have different characteristics depending on the environment and other factors. Table 6 presents the overall accuracy results and utility values. It shows the overall mean and variance of error metrics for each configuration considering all the images' results. Note that employing TSVD increases the average means and decreases the average variances and computing times for D2 and D4 where averages are calculated from seven different image pair outcomes. On the other hand, GF-added configurations produce higher mean values and computing times except configuration 5, but they generate variances changing in different directions. Additionally, although NLMD shows an improvement in mean values, it gives worse results in variance and time performances for D2. On the other hand, BF increases variance values even if it displays an improvement in time and mean outcomes for D2.

Overall Results
Considering that we work with inputs with different characteristics, the U3 and U4 values of the overall results are considered. If there is no concern about time performance, U3 is calculated considering the high accuracy and low variance that is desired for consistency on different images and parameters. U4 is obtained when the time performance is also taken into account. Table 7 clarifies the ranking of U3 and U4 in terms of configuration number. It is apparent that using one of TSVD, NLMD, GF, or some combination thereof raises the U3 value. Additionally, a higher U3 value is obtained when the absolute value is taken, but D8 (D9 is its absolute version) is an exception. Furthermore, D6 produces a higher U3 value than D8, D4, D2, and D1 do. Difference images with more information seem to work better. As an exception, D8 (a combination of D4 and D6) give a lower U3 value than D6 but still has a close value to D6. Nevertheless, D8 generates a better U3 value than D4, D2, and D1. When we using average time calculations, of course, the rankings change.
Since there are configurations with close utility values, more than one method can be selected and applied to scoring points for the follow-up plan. For example, utility values can be normalized between 0 and 1, and those above the value obtained by subtracting a certain percentage from the highest value can be selected. As such, configurations with high accuracy and low variance are selected for a set containing data with different properties. If time performance is also important, it is also counted. After determining the U3 and U4 values that fall within a certain percentage or above the threshold value, the configurations that are common to both sets can be selected.
According to all outcomes, we can answer the questions in Section 1. We find that it is possible to decrease the sensitivity (i.e., increase consistency). On the other hand, accuracy improves while the computation time is reduced for some, but not all configurations. However, no configuration works faster than the original method (config. 1) in terms of average calculation time. Despite this, in Table 6, there are configurations with high accuracy and average time calculation values that are close to the original method's result.

Conclusions
In this study, we compared the original PCAKM and its modified versions. All the configurations we use are deterministic, so the results are robust. In addition, none of them need the large training datasets required in supervised methods. Unsupervised methods, which work much faster than supervised methods, also stand out in terms of explainability. Today, issues such as explainability, interpretability, and transparency contribute to a trustworthy system [14], which is important for all stakeholders. Trustworthiness is an important concept to ensure that no undesirable consequences of AI systems occur during deployment.
Since PCAKM has more than one parameter combination and the analyses have different image types, it seems reasonable to look at the error metrics from an overall examination. As such, we have more consistent (i.e., less sensitive) information about the mean, variance, and time calculation performances of the error metrics. Since the difference between the error values to be obtained for all parameter combinations will be less, it will be more beneficial to use the combination of all of them. It is apparent that difference image and noise reduction makes a significant difference in the obtained results in terms of accuracy.
In the future, we plan to use the obtained results for point scoring in the follow-up activity, which may affect the road map of different agents [37,38]. A more consistent unsupervised method will help assign specific scores of interest to points on the map in a fast and efficient manner. For example, using different layer information such as the vegetation index, scoring can be done on change map information depending on the followup activity. In Figure A2, the figure on the left is an image taken from Google Maps and the figure on the right is the vegetation index [39] map (VIM) we calculate for this image. When a change map is produced for the image taken from Google Maps, the information to be obtained by overlapping the change map and VIM can play an important role in the scoring method according to the follow-up plan. As per the follow-up plan, VIM can belong to the first of the selected times for the change map, or it can belong to the second. In other words, it is determined by the follow-up action taken. Examples could be to investigate specific parts of the road infrastructure after an earthquake or flood.
Other maps similar to VIM can be used as labels to be illustrated on GIS. Maps that can be used for various situations (weather-related maps, information maps from the user or the planner, etc.) are merged with the change map as different layers to determine scores for the follow-up plan. Our next step will be to develop follow-up plan types and important layer maps that will contribute to the planning and scoring methods for each follow-up plan. At this point, it is worth noting the fact that SAR images are not affected by factors such as time zone and weather. Therefore, they offer an advantage in matters such as disaster response in terms of seeing the big picture. In addition, employing the proposed unsupervised method provides a robust, fast, explainable, less sensitive, and more accurate solution. These features will bridge the gap between scoring in AOIs for the follow-up plan that needs to produce a quick and estimated change map. We aim to develop an innovative scoring method for the follow-up action by merging change map results and other relevant map information as significant layers.
In addition, future work will aim to enhance the proposed change detection method with other unsupervised and supervised methods for different sensor types such as optical and thermal. In this way, the purpose of this is to obtain different layer maps by classifying the image [40]. These different layers will be used in scoring for the follow-up planning.
Author Contributions: The concept was developed through rigorous discussion among all the authors. D.K.K. took care of the coding parts for all computational aspects after a thorough discussion with P.N. All the authors were equally involved in the manuscript preparation. All authors read and agreed to the published version of the manuscript.
Funding: This research received no external funding.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available upon request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:          Figure A1. First two columns display SAR image pairs, third column shows ground truth change map, and fourth column illustrates predicted best change map results.
Google Maps image and its vegetation index map is illustrated in Figure A2.