CHANGE DETECTION IN MULTI-TEMPORAL IMAGES USING MULTISTAGE CLUSTERING FOR DISASTER RECOVERY PLANNING

Change detection analysis on multi-temporal images using various methods have been improved by many researchers in the field of spatial data analysis and image processing. Change detection analysis has many benefit for real-world applications such as medical image analysis, valuable material detector, satellite image analysis, disaster recovery planning, and many others. Indonesia is one of the most countries that encounter natural disaster. The most memorable disaster happened on December 26, 2004. Change detection is one of the important parts of management planning for natural disaster recovery. This article presents the fast and accurate result of change detection on multi-temporal images using multistage clustering. There is three main step for change detection in this article, the first step is to find the image difference of two multi-temporal images between the time before the disaster and after a disaster using operation log ratio between those images. The second step is clustering the difference image using Fuzzy C means divided into three classes. Change, unchanged, and intermediate change region. Afterword the last step is to cluster the change map from fuzzy C means clustering using k means clustering, divided into two classes. Change and unchanged region. Both clustering’s based on Euclidian distance. The accuracy achieved using this method is about 98.80%.


Introduction
Change Detection is a method in the field of spatial data analysis used for analyzing the area change of a location at a different time [1]. Based on the process method, change detection is divided into three categories, namely supervised, semi-supervised, and unsupervised change detection [2]. Supervised change detection is one of the classification-based change-detection which separates the changed and

Muhamad Soleh et.al, Change Detection In Multi-Temporal Images Using Multistage Clustering For
Disaster Recovery Planning111 unchanged areas that require training data. Determination of change by using clustering techniques which known cluster the number is called semi-supervised detection, while unsupervised detection is one of the change detection techniques that use clustering method without knowing the number of its clusters. Based on the input dataset resolution, change detection is categorized into two, both are objectbased change detection and pixel-based change detection [2]. Object-based detection is a technique to analyze the change of an area based on the change of objects found in the area being observed. Each object located in that area is categorized by segmentation of a certain pixel range value, whereas pixel-based change detection uses datasets with small to medium resolution. The dataset is used to determine the area of change and unchanged. The area is classified based on the maximum pixel value and the minimum pixel value on the change map.
Some researchers with the topic of interest on remote sensing have done some research related to change detection, using either supervised [3], semisupervised [1] or unsupervised methods [2]. The dataset used for change detection research is the dataset of remote sensing image. Each remote sensing uses either active sensor such as Synthetic Radar Apparatus (SAR) or passive sensor or example like Optical Landsat Satellite. The most pixel-based change detection consists of two processes. The first step is to get the image difference from multi-temporal images. The second step is the filtering area that undergoes changes by using binary masking. The most common methods used in the first phase to find image different those are image subtraction [1], image log ratio [3], change vector analysis [4], correlation coefficient [5], similarity structure index measurement (SSIM) [2] or other methods. As for the binary masking process to get change and unchanged areas, most researchers use optimization [2,[5][6][7], thresholding [8], clustering [1,[9][10], and classification methods [3,[11][12].
To find the best method to use, most researcher decide based on both types of the dataset and the purpose of the study. Based on the result of some comparing research on change detection method, log ratio is better used for the remote sensing dataset which used an active sensor, while for the passive sensor better-using image subtraction method [3]. In 2013, Ghosh et al proposed a new method of change-detection research using the Change Vector Analysis (CVA) to see the inter-band vector changes in multispectral images. CVA combined with Gibbs Markov Random Fields by sampling the CVA results to determine the difference image based on the previous pixel relationship. The binary masking process proposed by Ghosh is an optimization method using expectation maximization by applying Hopfield-type neural networks [12].
In 2016, Feng Gao et al proposed a method of change detection using deep learning based. Feng Gao et al used three different datasets, these are San Francisco, Ottawa, and Yellow River Datasets. The difference image gets by using log ratio. The result then transformed into orthogonally using Gabor wavelet as feature extraction, feature vector produced by Gabor wavelet used by Fuzzy C Means Clustering to produce change map which has three clusters consist of change, unchanged, and intermediate change. The binary masking in this research uses PCA-Net as feature extraction. Then classified into change and unchanged area using Support Vector Machine (SVM) [3].
In 2017 A.Yavariabdi publishes an article on change detection topic using the Structure Similarity Index Measurement (SSIM) and Multi-objective Evolutionary Algorithm (MOEA). The SSIM used to get the difference image which can handle luminance, contrast, and the structure of the image input, while the MOEA Non-dominated Sorting GA (NSGA II) is used for binary masking for determining the changed region then enhanced by using Markov Random Field Fusion [2].
One of the useful implementation in the real world by doing this research related to change detection is for disaster recovery planning management [13]. There are several factors to consider in recovery after natural disasters happened, such as the number of personnel to be prepared, the supply of food needed, the expected time required, the availability of access road to the location and so on. Indonesia is one country among other countries that often happened the natural disasters.
Based on the data reported by the national disaster management agency, Republic of Indonesia, during the period of January to June 2017, Indonesia has an event 1368 natural disaster. The number of housing damage is 18,983 units. The death toll and disappeared reached 227 people and the number of injured and displaced around 1,710,973 people. Total of Disasters, Victims, and the Impact until June 2017 can be seen in table 1 and table 2 [14]. The purpose of this study is to find the effective and efficient method for change detection. Hopefully, it is can be used for disaster recovery planning.

Methods
The dataset used in this experiment consists of the Ottawa Dataset and the San Francisco Dataset. The Ottawa dataset is obtained from the RADARSAT SAR sensor, where two images were taken before and after the flood [3]. San Francisco Dataset obtained from the ERS-2 SAR sensor, the dataset shows the state of San Francisco city to see changes in urban flow that occurred in the city in the ten- using Absolute Log Ratio, optimize difference image using Fuzzy C Means Clustering, and Binary Masking using K-Means clustering.
Change detection is one of the methods that can use remote sensing to detect the change in the landscape over time. Its method can separate an interest region of changing area over time. There are several methods that have been implemented to find the change region on multi-temporal images. Absolute log ration is one of the most useful basic methods to separate changing region from two inputs with different time at the same location. Absolute log ratio is an extended method of image ratio. With the log-ratio operator, the multiplicative of speckle noise can be transformed in an additive noise component. Furthermore, the range of variation of the ratio image will be compressed and thereby enhances the lowintensity pixels [16].
The logarithmic scale is characterized by enhancing the low-intensity pixels while weakening the pixels in the areas of high intensity; therefore, the distribution of two classes (changed and unchanged) could be made more symmetrical. However, the information of changed regions that are obtained by the log-ratio image may not be able to reflect the real changed trends in the maximum extent because of the weakening in the areas of high-intensity pixels. Therefore, to obtain the change map from multitemporal images, we need to optimize by combining with the other method. A clustering algorithm is one of the methods that can reflect the real changed and unchanged region. Specifically, this method incorporates the information about spatial context to the corresponding objective function for reducing the effect of speckle noise [17].
In this research we use multiple clustering. This combine the Fuzzy C Means and K Means Clustering. FCM used for pre-classification. The main purpose of pre-classification is to find pixels that have high probabilities of being correctly classified into changed or unchanged classes. In fuzzy clustering, each point has a probability of belonging to each cluster, rather than completely belonging to just one cluster as it is the case in the traditional k-means. K-means clustering used for binary masking for producing change map. Kmeans classify the intermediate region area from FCM into change and unchanged area.

Difference image using Log Ratio
In the first stage, two multi-temporal satellite images 1 and 2 in the same location as input images. Difference image ( ) is obtained by using the Log Ratio.
Difference image is used to see the difference between two image inputs. The result of difference image then processed to extract the feature using multiplication operation between and Gabor kernel. Gabor wavelets concern of frequencies and orientations, where the filters used in the wavelet Gabor are composed of two components, real and imaginary components that represent orthogonal directors.
Where is Gabor real component and is Gabor imaginer component with is wavelength waves, is spatial aspect ratio, is standard deviation, and is wave phase. Where ′ and ′ denoted as the equation 4 and 5.
By taking the magnitude of the real component and imaginer Gabor filter it will generate a feature vector. = [ 1 , . . . , ]. Where M and N represent the length of the width of . Feature Vector obtained from Gabor then used as input for the clustering process using Fuzzy C Means Clustering.

Optimize Difference image using Fuzzy C Means Clustering
At this stage, the feature vector is processed three processes. The first process is to divide the feature vector into two clusters representing potential areas as change region ( ) and unchanged ( ), but the number of clusters used as many as 5 clusters. The second process as the same as in the first process, but the number of clusters used is 5 clusters ( 1 , 2 , 3 , 4 , 5 ). While the third process is to divide the five clusters into 3 clusters that represent the change area ( ), intermediate change ( ) and unchanged ( ). Based on the value of a certain objective function. In the Fuzzy C Means algorithm the clustering process is done by minimizing the objective function value during the iteration.

Where
_ is change map fuzzy with c is cluster, x is feature vector [ 1 , . . . , ], m is fuzzifier, m ϵ R with m >=1. Where (Membership Value) denoted as the equation 7.
CM_F is a fuzzy change map consisting of three class labels, namely change, intermediate, and unchanged. CM_F is then transformed orthogonally using Principe Component Analysis (PCA) to obtain a new feature vector. The process of PCA transformation in CM_F also consists of three processes including 1). Make h x h non-overlapping block from difference image; 2) Creating an eigenvector space using PCA on the h x h nonoverlapping block already created; 3) Create a feature vector of CM_F by projecting h x h nonoverlapping block into Eigen vector space. The process of obtaining an eigenvector on PCA uses the Singular Value Decomposition (SVD) principle of the feature vector input, to produce a feature vector output resulting from a different image projection that has been restored into its original dimension.

= .
() Where T is Feature Vector Output, X is Feature Vector Input, and W Eigenvector. Using SVD method, X value can be written as the equation 9.
Where U is the orthogonal unit vectors and is the Rectangular diagonal matrix of positive numbers. So the value of T can be rewritten as the equation 10.
The result of the vector T feature is then used for the clustering process. Clustering process aims to determine the change and unchanged region and then produce a change map (CM).

Binary Masking using K-Means clustering
Basically, the k-means clustering principle is the same as fuzzy c means clustering. Both methods use distance and centroid as cluster boundary and majority vote as new raster determination.
is Change Maps, k is numbers of cluster, s is a set of clusters from input data and is a centroid.

Results and Analysis
The experimental results of this study consist of two parts. Performance measurement test and Implementation of the proposed method.

Performance test of the proposed method
The first section we will discuss the performance test of the proposed method using baseline dataset obtained from previous research [3]. The results of the proposed methods are compared with preexisting methods such as PCA-k-Means [1], MRF-FCM [9], Gabor-TLC [10], and D_MRF-FCM [11]. The parameters used for measuring the performance of proposed algorithm are Miss Alarm (MA), False Alarm (FA), Overall Error (OE), and Accuracy (ACC). MA is the total number of pixels of change region that is not successfully detected as a change region, FA is the total number of pixels in the unchanged region detected as pixels in the change region, OE is the sum of FA and MA, while the ACC is the ratio between OE and the total pixels then minus 1 and multiplied by 100%.
Based on the results of the performance obtained in Table 3, it is known that the proposed method is one of the promising methods for performing the change detection. So, it`s can be processed to be applied to various applications that can solve some real-life problems. One of them is disaster recovery planning. Compared to the previous research, the proposed method has a good quality in Accuracy, False Alarm, and Overall Error. We also have provided other calculation to measure the performance of the proposed method. False alarm rate (PFA), Missed detection rate (PMD), percentage correct classification (PCC), and Kappa coefficient (KC) can be used for the indicators to measure the effectiveness of different methods [3].
In the reference ground-truth map, the actual numbers of pixels belonging to the unchanged and changed classes are denoted by NU and NC, respectively. Then, PFA is FA/NC × 100% and PMD is MD/NC × 100%, in which MD refers the number of changed pixels are detected as unchanged. The PCC is (NU + NC − FA −MD) / (NU + NC) × 100%. The KC is a statistic which measures inter-   rater agreement for qualitative (categorical) items. It is generally thought to be a more robust measure than simple percent agreement calculation, as κ considers the possibility of the agreement occurring by chance [10]. Table 3 shows the result of the proposed method compared to the others. For the San Francisco dataset, the result shows 0.15, 13.15, 98.80, and 87.45 for PFA, PMD, PCC, and KC respectively. Otherwise, for the Ottawa dataset, the result shows 1.64, 5.93, 97.67, and 93.66 for PFA, PMD, PCC, and KC respectively. Table 4 shows PCA-Net proposed by Feng Gao et al performance better results than the proposed method. But in computing time, PCA-Net is not very effective. To perform a change detection analysis process that requires real-time results, especially in disaster recovery planning based on the experiments that have been run the time required for running PCA-Net is around 2274.69 to 3828.17 seconds. While on the proposed method, the running time required only about 9.42 to 28.51 seconds. The ratio of computational time between PCA-Net and the proposed method is very much different. By using multistage clustering, the result of change detection is near the performance of PCA-Net, even the proposed method is very efficient in terms of its algorithm complexity. PCA-Net composed of three steps: the first two-stage are PCA filter convolutions and the last stage is hashing and histogram generation. PCA filter convolutions composed of several layers. PCA-Net consist of training and testing phase. The sample of change and unchanged images were generated from the images. Those needs more time to classify the change and unchanged area.
In Figure 3, it shows the visualization for each method used in this research. There are two regions which indicated changed and unchanged areas from multi-temporal images. Several methods are still not optimal for performing changes map in the different image. There are still too many speckle noise in some areas. Based on the conventional method of change detection, the proposed method shows the more effective result to find the changing area. Although the result still has a speckle noise but compare with the other method, the speckle noise from the proposed method is less than the other conventional method for Ottawa and San Francisco dataset. The best method for reducing the speckle noise is based on deep learning technique. The visualization of the image also proves that the values contained in table 3 and 4 correspond to the visualization results shown in Figure 2.

Implementation of Proposed Method on Aceh Tsunami Dataset
This section we will discuss the implementation of the proposed method using new datasets. The research consists of preprocessing datasets, dividing the regions into three clusters, binary masking in the intermediate change area into two clusters, thus testing the image quality based on structure, luminance, and contrast. Preprocessing on the dataset is done by cropping in the region of interest. This cropping process is an important part of this research. The number of pixels on both data sets should be the same so the results can be compared. The characteristics of the dataset are made differently, this is intended to differentiate the changing and unchanged regions.
The dataset used is from remote sensing images located in Lhoknga, Aceh Province taken on January 10, 2003 (before the tsunami) and 29 December 2004 (3 days after the tsunami), another dataset is in Gleebruk, Aceh Province taken on 12 April 2004 (before the tsunami) and on 2 January 2005 (some after the tsunami). The input images used in this research based on the result of preprocessing data with the size of 705 x 705 pixels, so the total pixel used is 497,025 pixels. The characteristics of Lhoknga area used are residential areas on the seafront, while the Gleebruk area is a middle residential area.
The next step is to calculate the value of similarity index measurement (SSIM) in each input dataset. The SSIM used in this study is the method proposed by Z. Wang for image quality assessment [15]. The SSIM value in the range of 0 -1. The closer the two-image resemblance, the SSIM value is close to 1, whereas if the two comparable images have a slight resemblance, the value of SSIM is shown to be close to 0. SIM value is obtained by using equation 12.
( , ) = (2 + 1 )(2 + 2 ) ( 2 2 + 1 )(2 2 2 + 2 ) Where is the average of x, is the average of y, 2 is the variance of x, 2 is the variance of y, is the covariance of x and y, 1 is ( 1 ) 2 ; 1 2 is ( 2 ) 2 are two variables to stabilize the division with weak denominator with is the dynamic range of the pixel-values and 1 is 0.01; 2 is 0.03 by default.
The value of SSIM then used as a baseline in analyzing the result of change detection by using the proposed method. The relationship between the SSIM value and research of change detection is certain closely related. The closer the resemblance of the two image inputs, the region undergoing the change after the tsunami disaster will be less, whereas if the value of SSIM approaches 0, then it indicates a larger area change in the image. Figure 3 shows the visualization of the proposed method using Lhoknga and Gleebuk area before and after Tsunami happened in Indonesia, Especially in Aceh Province. For each dataset, consist of two images. Figure 3a and 3b show the images of Lhoknga and Gleebuk dataset before and after Tsunami. Figure 3c and  For examples such as finding safety area for mitigation process, selected area for emergency reconstruction building, or effective food distribution and first aid personnel on the victim etc. Table 5 shows the comparison between Lhoknga and Gleebruk datasets using the proposed method. In the dataset located in Lhoknga, the total of all pixels, there has been a change is about 160.658 pixels. While on the Gleebruk location the total pixels that change reached 357.332 pixels. From both area, we suggest to the disaster recovery team to concern more in Gleebruk area which has more damage than Lhoknga, although we consider there are other factors may be influence.
Based on the structural similarity index measurement (SSIM) on both image inputs, it is known that in the Gleebruk area, the similarity values show the result is very small, even to the point of 0.10. This is because the image before and after the tsunami has increasingly changed. Whereas in the Lhoknga dataset, the result of SSIM value is still in the range of 0.38. It shows there is still similarity between the image before and after the tsunami. SSIM uses a comparison of structure, luminance, and contrast on two images.

Conclusion
The Multistage clustering method proposed in this article has succeeded in detecting a region's changes in multi-temporal images with satisfactory performance results. The accuracy achieved using this method is about 98.80%. It was better than another result and close to deep learning-based method, with more efficient time complexity. The proposed method has also been implemented in Tsunami dataset that ever happened in Indonesia precisely located in Lhoknga and Gleebruk, Aceh Province. The future work of this research will implement the deep learning technique such as the convolutional neural network to find the change and unchanged region form multi-temporal images.