Unsupervised Change Detection in Landsat Images with Atmospheric Artifacts : A Fuzzy Multiobjective Approach

A new unsupervised approach based on a hybrid wavelet transform and Fuzzy Clustering Method (FCM) with Multiobjective Particle Swarm Optimization (MO-PSO) is proposed to obtain a binary change mask in Landsat images acquired with different atmospheric conditions.The proposed method uses the following steps: (1) preprocessing, (2) classification of preprocessed image, and (3) binary masks fusion. Firstly, a photometric invariant technique is used to transform the Landsat images from RGB to HSV colour space. A hybrid wavelet transform based on Stationary (SWT) and Discrete Wavelet (DWT) Transforms is applied to the hue channel of two Landsat satellite images to create subbands. After that, mean shift clustering method is applied to the subband difference images, computed using the absolute-valued difference technique, to smooth the difference images. Then, the proposed method optimizes iteratively two different fuzzy based objective functions using MO-PSO to evaluate changed and unchanged regions of the smoothed difference images separately. Finally, a fusion approach based on connected component with union technique is proposed to fuse two binary masks to estimate the final solution. Experimental results show the robustness of the proposed method to existence of haze and thin clouds as well as Gaussian noise in Landsat images.


Introduction
Change detection is the process of identifying differences between multitemporal satellite images over time and it has been used in many important applications such as agriculture [1], geology [2], forestry [3], regional mapping and planning [4,5], oil pollution monitoring [6], and surveillance [7].In remote sensing context, many automatic change detection algorithms have been proposed to obtain changed and unchanged pixels between images.Amongst these methods, unsupervised approach has been widely and effectively used in change detection problem as it does not need any training data set and prior knowledge of images [8].Comprehensive surveys of unsupervised change detection methods can be found in [9][10][11].
Landsat imagery is one of the most frequently used remote sensing techniques to observe and monitor Earth's surface change due to its medium spatial resolution and frequency of data collections.In Landsat satellites, the images are acquired using passive optical and multispectral sensors, so that this remote sensing technique is very sensitive to any meteorology conditions [12,13].Due to this drawback, atmospheric artifacts can be depicted in the Landsat images which makes the change detection problem very challenging.Generally, in the literature, to find changed and unchanged pixels in medium resolution Landsat images, pixel-based change detection approach has been used and it mostly consists of two phases [14,15]: (1) Difference image estimation, where it is computed by comparing two images through a similarity metric method such as: image subtraction [14,16], Change Vector Analysis (CVA) [15,17], Correlation Coefficient (CC) [18], and Erreur Relative Globale Adimensionnelle de Synthese (ERGAS) [19].(2) Image analysis, which aims at obtaining binary change mask using clustering [14,20,21], Level set [17], population-based optimization [16,18,22], OTSU and Bayesian thresholding [19,23], and many others [6,[24][25][26].
Pixel-based unsupervised change detection methods have been applied into two different satellite image domains which are spatial domain and frequency domain.In Landsat change detection context, most of the existing methods are based on the spatial domain.For instance, Celik [14] proposes an unsupervised change detection technique for grayscale Landsat images.In this method, firstly, image subtraction is used to estimate similarity between two images.However, this technique is sensitive to noise and atmospheric condition changes so that it cannot be used on the raw images without any adjustment or modification.In [14], to improve the robustness of the change detection method against noise, Principle Component Analysis (PCA) is used to find the statistical properties of the difference image.Finally, -means clustering approach is applied to the PCA feature vectors to obtain change detection mask.Hao et al. [17] propose a change detection approach based on Expectation Maximization (EM) and level set methods.In this method, a difference image is initially generated using the CVA and is modeled using Gaussian Mixture Models (GMMs).The EM algorithm is employed to estimate the parameters of the statistical model.Then, an energy functional based on changed and unchanged areas in the difference image is defined and optimized using a level set method.However, in Landsat change detection problem, this method may provide high error detection rates as the CVA is highly dependent on a radiometric correction technique [29].Moreover, other unsupervised methods such as fuzzy -means [21], normalised cut clustering [20], and fuzzy active contour model [24] have also been applied to Landsat images to find binary change mask.However, these methods only find the significant changes between two input images and have lack of adaptability and flexibility to obtain weak changes.Therefore, optimizationbased change detection methods are proposed to find optimal binary change mask by obtaining both weak and strong changes on the difference image.For instance, Celik [22] uses Genetic Algorithm (GA) to cluster changed and unchanged pixels between two grayscale satellite images by minimizing a cost function based on the difference image computed using the image subtraction technique.Kusetogullari et al. [16] propose a parallel processing unsupervised change detection approach which simultaneously employed several Particle Swarm Optimizations (PSOs).This method is applied to the difference of three channels of RGB colour intensities to determine the optimal change detection mask.Most of the existing optimization-based change detection methods have two main drawbacks.(1) Optimization is performed on the raw difference images which increases computational cost to converge to the solution.(2) A cost function based on a weighted sum of the changed and unchanged objectives is used which may not provide optimal trade-off between changed and unchanged objective functions.In other words, these methods can easily get stuck in a local minimum due to those two issues.
Yavariabdi and Kusetogullari [30] proposed a multiobjective change detection method to obtain changed and unchanged pixels between two multitemporal multispectral Landsat images using multiobjective evolutionary algorithm (MOEA).The method minimizes two different objective functions which are not fuzzy using NSGA-II (Nondominated Sorting Genetic Algorithm II) and provides trade-off between the objective functions.Besides this, SSIM technique has been used to obtain the difference image which is optimized by creating and evolving binary change detection maps or possible solutions using NSGA-II.The results show the effectiveness of the NSGA-II over the existing change detection methods.However, the method has been performed in the spatial domain which is not robust to the noise on the satellite images.In this paper, a novel change detection method using Multiobjective Particle Swarm Optimization (MO-PSO) has been proposed which has been applied to the high and low frequency levels (subbands) of the two satellite images obtained using hybrid wavelet transform.Thus, the proposed approach is more robust to the noise.Another issue in [30] is that NSGA-II uses evolutionary techniques (crossover, mutation, selection, etc.) to find the most optimum result.Using a number of evolutionary techniques increases the computational time to achieve the best result but the proposed method MO-PSO is a computational based approach which achieves the results quicker and more efficiently.Furthermore, we develop fuzzy based objective functions in this paper.
Although unsupervised change detection methods are mostly performed in the spatial domain, it can also be carried out using the frequency representation of the Landsat images [23,25,31].For instance, Chen and Cao [25] propose an unsupervised change detection method based on the SWT and active contours algorithms which are robust to noise.Celik and Ma [23] propose a frequency-based change detection method for noisy grayscale Landsat images using the dual-tree complex wavelet transform and Bayesian thresholding approach.In this method, image subtraction is used to estimate similarity between corresponding highpass subbands.Next, the Bayesian thresholding method is used to find binary masks.Finally, a fusion based approach, which is the union of all the binary masks, is used to obtain the final result.However, this method suffers from two main problems: (1) neglecting the low-pass subband images which causes the removal of many details from images and (2) using simple fusion technique which can easily lead to high false alarm rate.
In most of the state-of-the-art methods, the effects of atmospheric artifacts on change detection problem are neglected as they assume that the Landsat images are radiometrically corrected.However, this is an important issue as variations of atmospheric conditions (i.e., existence of haze, fog, smoke, or thin cloud) in images can change the contrast and fade the colours in the image.Therefore, most of the existing change detection methods cannot be applied into the Landsat images which are captured with various atmospheric conditions as they can wrongly detect those artifacts as changed area(s).In this paper, we propose a new unsupervised change detection method which is robust to the varying atmospheric conditions in Landsat images.To achieve this, the proposed method uses a photometric invariant technique along with a hybrid wavelet transform and a density-based clustering (mean shift algorithm) to decrease the effect of atmospheric artifacts in the input images and then two cost functions based on image subtraction are generated and simultaneously minimized using a multiobjective population-based optimization.The method first uses a photometric invariant technique to convert RGB satellite images into the HSV colour space.Note that, only hue channel of Landsat images are used as they are invariant under both shadow and shading (i.e., illumination intensity changes) as well as highlights and specularities.Next, a hybrid wavelet method based on the SWT and the DWT is applied to the hue channel of the images, separately.Thus, for each satellite image, there will be eight subband images, four of which are obtained using the DWT and the rest using the SWT.Here, the LL (Low-Low subband) and HH (High-High subband) subbands of the DWT and the SWT of each Landsat image are only used for further processing.After that the mean shift clustering method is applied to the difference subband images obtained using the image subtraction technique to decrease the intensity variations in the images and preserve the high frequencies.This simply decreases the computational cost and improves the convergence rate of the population-based optimization algorithm which are the two main issues of the most conventional optimizationbased change detection methods (e.g., [16,18,22]).In the next step, for each smoothed difference subband image, two cost functions based on FCM are generated and optimized using MO-PSO.Here, two sets of change detection masks with different trade-off relationships between the costs are obtained and the median of each set is selected to obtain two binary change masks.Finally, in order to obtain the final binary change detection mask, a fusion method based on a connected component with union technique is used.The main contributions of the paper are summarized as follows: (i) Proposing a new frequency-based unsupervised change detection method for Landsat images which are captured with various atmospheric conditions.
(ii) Combining a photometric invariant technique with wavelet transforms to decrease the influence of atmospheric conditions on the change detection results.
(iii) Presenting two fitness cost functions based on FCM which are robust to noise, haze, and thin cloud(s) and being used in the MO-PSO.
(iv) Presenting a new procedure that focuses on decreasing the computational cost of the population-based optimization algorithms and improving their convergence rate.
The rest of the paper is organized as follows.Section 2 proposes our new change detection approach.Section 3 provides the results of our experiments and a discussion related to the quantitative and qualitative change detection results.Lastly, in Section 4, the paper is concluded.

Proposed Method
Generally, unsupervised change detection methods have been applied to the Landsat change detection problem with the following assumptions: (1) the data must be captured at anniversary dates, which greatly minimizes the effects of seasonal changes (i.e., leads to illumination variations) and (2) the data must be clear from clouds, haze, smoke, and extreme humidity as they introduce erroneous results to the change detection process.Therefore, in Landsat change detection problem, existence of such artifacts continues to be a major challenge.Note that, in this paper, we assume that the data sets are not captured at anniversary dates and may include thin cloud(s).One way to remove the effects of those artifacts is to use atmospheric correction and cloud removal algorithms [28,32].However, this strategy does not seem to remarkably benefit the change detection methods, and in some cases it even degraded accuracies.For instance, Figure 1 shows the influence of an atmospheric correction algorithm [28] on a change detection method.It is clear that this strategy partially removes the cloud from the image (Figure 1(a)), but it fades the colours in the image (Figure 1(b)).Figure 1(d) demonstrates that the change detection method [14] using Figures 1(b) and 1(c) incorrectly detects many unchanged pixels as changed pixels (red pixels).In order to solve this problem, we propose a new change detection method which is robust to variations of atmospheric conditions.

Preprocessing.
In the proposed method, three preprocessing approaches as shown in the first part of the block diagram in Figure 2 are applied to the coregistered Landsat images acquired from the same geographical area, but at two different time instances.In order to obtain the best change detection mask for the various atmospheric conditions, the following procedure has been used.Firstly, a photometric invariants technique is applied.Here, the strategy for using the photometric invariants is to convert RGB colour space into HSV colour space, which is obtained via conical transform.In this space, each colour represents as hue, saturation, and value.

Mathematical Problems in Engineering
Note that, most of the existing change detection methods use the grayscale or RGB colour images to find the change detection mask so that they cannot overcome the typical atmospheric artifacts occurring during Landsat acquisitions.Secondly, the hue channels of two images are transformed into the frequency domain using the SWT and the DWT, which makes the proposed method (1) robust to noise; (2) robust to global illumination variations; and (3) insensitive to existence of haze and thin clouds and (4) increase the accuracy of detection near boundaries as it is not based on spatially local image luminance differences.Finally, absolutevalued difference technique is used to find the differences between LL and HH subbands.

Colour Space Transformation.
Many photometric invariants such as normalised RGB, HSI, HSV, , and derivatives of the logarithmized colour channels have been proposed in literature [33].In general, the HSV colour transformation proves to be insensitive to shadow, shading, and specular edges [33] and has been used in different applications such as object segmentation [34], detection [35], tracking [36], and recognition in the images [37].More details about HSV colour transformation can be found in [38].In this colour space, the hue () resembles the pure colour, the saturation () shows the amount of gray in the colour, and the value () describes the brightness.In other words, the hue channel allows coping with shadow, shading, highlights, and specularities, the saturation channel is only invariant under shadow and shading, and the value channel cannot cope with any of these artifacts.Therefore, the hue channel is used in the proposed method.
Let   1 and   2 be two Landsat images with size of ℎ ×  pixels, where ℎ and  are height and width of images, respectively.They consist of three spectral bands (), i.e., , , and , which are used to obtain the hue channel () as follows: where  = {1, 2} denotes the corresponding Landsat image.

Hybrid Wavelet Transform.
A wavelet transform method decomposes an image into images or subbands which cover the frequency components or coefficients of the original image.Besides this, each obtained subband contains different details of information of the original image.For instance, the subband HH holds the detailed information (e.g., edges) and LL holds the most of the information of the image.In the change detection problem, unsupervised change detection methods can find the changed and unchanged pixels as well as undesired artifacts using the LL subband since it covers the most of the information of the image.In the corresponding problem, the challenging tasks are to detect the changed and unchanged pixels in the subbands and ignore both the noise and effect of the atmospheric variations (e.g., haze or thin cloud) on Landsat images.In order to find a better change detection result, a hybrid wavelet transform based on the DWT and the SWT is applied to the hue channel   of input images   .
Let the DWT and SWT decompose an image   into low and high frequency subband images, respectively.In this case, four subbands are obtained using the DWT ( LL, ,  HL, ,  LH, , and  HH, ) and four subbands are computed using the SWT ( LL, ,  HL, ,  LH, , and  HH, ).Note that, in the proposed method, only the LL and HH subbands ( LL, ,  HH, ,  LL, , and  HH, ) are used.Besides this, we empirically found out that the binary masks of LH and HL intersect with the binary mask of LL subband and they include many isolated pixels (noise) in nonintersection areas so that these subbands have not been used for further processing.
Unlike the SWT, the DWT uses downsampling to create the subband images; thus it is necessary to make the size of the selected subbands equal.To this end, the  LL, and  HH, subband images are interpolated using bicubic interpolation.Then, in order to strengthen the coefficients in the subband images which results in improving the image details, the mean of the subband coefficients (   ), where  denotes LL and HH, is used.Finally, the absolute-valued subtraction method is used to find two difference images.
2.2.Mean Shift Clustering.Generally, conventional unsupervised change detection methods use the natural difference image, which contains trivial details, to cluster pixels as changed and unchanged.The clustering techniques such as -means [14], EM-GMM [17], and fuzzy -means [21] used in change detection problem suffer from the following limitations, especially when the difference image contains trivial details: (1) the spatial structure and edge information of the difference image cannot be preserved and (2) pixels from disconnected regions of the difference image are grouped together if their feature space overlap.These limitations simply lead to low correct detection rates.To avoid these difficulties in the change detection problem, an edge preserving smoothing technique which aims to remove the trivial details while preserving major difference image

Combining sub-bands
Combining sub-bands structures must be used.To achieve this, mean shift algorithm is used as discontinuity preserving filtering on the difference subband images (see the second part of the block diagram in Figure 2).Mean shift is a nonparametric clustering algorithm, which is based on maxima of densities functions and has been applied to edge preserving smoothing, image segmentation, and tracking problems [39].In contrast to the other clustering techniques, mean shift clustering does not have embedded assumptions on the shape of the distribution and does not require the number of clusters.In this algorithm, a fixed size Parzen window is iteratively shifted to the closest stationary point along the gradient direction of the estimated density function.In the change detection process, each difference image pixel is changed by the mode of its intensity neighbours' values computed in the natural difference image.Consequently, this approach can significantly reduce the intensity variations while preserving salient features of overall difference image.For ℓ 1 × ℓ 2 sample pixels   (, ),  = 1, . . ., ℓ 1 , and  = 1, . . ., ℓ 2 , the mean shift vector is estimated as follows: where   is the center of the kernel for each subband and starts at an arbitrary value, (  ) is the densities estimation of sample pixels in the neighbourhood of   , ℓ 1 = 5 and ℓ 2 = 5 are the height and width of the kernel, respectively, ℏ > 0 is the bandwidth parameter, and () = exp(−(1/2)) denotes Gaussian kernel function.More specifically, the weighted function () adjusts the influence of each intensity value inside Parzen window according to its distance from   .Note that (3) is optimized using gradient ascent technique, so that at each iteration   is updated to (  ) until ‖ ℏ (  )‖ < , where  is a prespecified convergence threshold set to  = 10 −5 in our experiments.Note that the stationary points obtained via gradient ascent represent the modes of the density functions.All the pixel values associated with the same stationary point ( ℏ ) belong to the same cluster which leads to a filtered image denoted as .
In the proposed method, the mean shift approach has been applied to two different subband images and, thus, two smooth difference subband images are obtained.Unlike the other conventional change detection methods, the generation of binary change mask is not based on the pixels of the original difference image but on the results of mean shift algorithm.Combining mean shift algorithm with fuzzy- means clustering and Multiobjective Particle Swarm Optimization may bring several advantages.One advantage is that the mean shift algorithm significantly reduces the noise and the amount of smoothing for the abrupt changes in the edges as well as preservingthe salient features of the overall image.A second advantage is that since the number of basic image entities is much smaller than the original image, the computational complexity is considerably reduced.The next stage of the proposed approach is the generation of two clusters by clustering the feature vector space using the Multiobjective Particle Swarm Optimization algorithm.

Fuzzy Based Multiobjective Particle Swarm Optimization.
Change detection problem is an optimization problem and it has more than one objective function such as changed objective function and unchanged objective function.A multiobjective optimization problem can be formulated as follows: where  = ( 1 ,  2 ,  3 , . . .,   ) is the vector of decision variables and   () ( = 1, . . ., ) denotes the th objective to be minimized in the multiobjective optimization problem.
In most of the existing change detection algorithms, these objective functions are combined to be used as a single optimization function.Thus, the single-objective function can be used in Genetic Algorithm or Particle Swarm Optimization methods.This solution strategy may bring several issues: (1) improvement of one objective since it may cause deterioration for the other objective when the optimization method performs obtaining the global optimum result; (2) getting stuck with a local optimum result easily; (3) only one optimum result (e.g., local or global) being obtained at the end of the iteration.Therefore, using weighted sum of the objective functions is not an efficient approach to solve an optimization problem.Unlike the single-objective optimization solution approach, multioptimization solution approach provides nondominated solutions known as Pareto optimal solutions [27,40,41].Therefore, multiobjective optimization algorithms are used to (1) discover solutions as approximated to the optimal solutions as possible; (2) find solutions as uniform as possible in the obtained nondominated front; (3) determine solutions to cover the true Pareto Front (PF) as broad as possible.In order to achieve these three challenges, it is necessary to provide the best trade-off solutions or Pareto optimal solutions between the objectives.Pareto solution is defined as follows [27,41]: Multiobjective algorithms have been used and developed to resolve different optimization problems and they perform superiorly to resolve the multiobjective problems compared to the single optimization algorithms [27,40,41].In this paper, a new unsupervised Fuzzy Multiobjective Particle Swarm Optimization (Fuzzy MO-PSO) has been used to resolve the change detection problem.The pseudocode of Fuzzy MO-PSO is given in Algorithm 1.
First, the positions of particles,   , are initialized with binary values with the size of Landsat images  × W, where  is the height and  is the width of the image.In the algorithm,   is the position of th particle, and velocities of particles, V  , are set as 0s, where V  is the velocity of th particle.It is important to note that particles move or change positions to obtain the global optimum result based on the particle's own best results,    1, and   2, , and the overall best results of all particles,   1, and   2, , so far.The particle's own best results and the overall best results of all particles are updated in each iteration using the step numbers from (8) to (22) in Algorithm 1. Unlike the traditional PSO approach which is used in a single optimization problem [42], two different fitness values are obtained in the multiobjective problem.Therefore, there will be two best results for each particle,   1, and   2, , and two overall best results of all particles,   1, and   2, , with respect to fitness functions 1 and 2. The algorithm works by updating the positions of particles after applying ( 5) and ( 6), respectively.Then, the position values are used in (7) to create the binary change detection masks.Obtained updated particles are added in the new population .Then, both the current  and the new population  are joined; the resulting population, , is ordered according to Pareto solutions of particles by using ranking operator and a density estimator known as crowding distance [27].Thus, nondominated and dominated solutions are evaluated and the nondominated set of solutions are inserted into the population .These steps are repeated until the termination condition is fulfilled.
The position and velocity of th particle are estimated as follows: where  is an inertia coefficient;  1 and  2 are the weights of learning rate of a local optimum and learning rate of a global optimum, respectively;  1 and  2 are random parameters ranging from [0, 1].The position values are used in the following equation to create the binary change detection masks.
The purpose of the proposed method is to determine the set of Pareto optimal solutions.Each particle in each iteration contains two different fitness values.To solve the multiobjective change detection problem using the Fuzzy MO-PSO, it is necessary to estimate fitness values of each (1) Initialize the positions of particles with randomly generated binary values, and set velocities of particles to 0, the generation number  = 0, population size .(2)  ← initialized Binary Population Set (3)  ← 0 (4) while Stopping criteria is not reached do (5) For each particle  in  do (6) Compute the fuzziness fitness functions   1 and   2 using the binary valued positions.(7) end for (8) Compare the particle's current finesses   1 and   (23) Update the velocities of particles using Eq. ( 5) (24) Update the positions of particles using Eq. ( 6) (25) Insert the updated particles into the  (26) Create the binary change-detection mask using Eq. ( 7) (27)  ←  ∪  (28) Use the Ranking and crowding distance operator [27] Select Pareto solutions or non-dominated solutions (30) Insert selected Pareto solutions into the  (31)  =  + 1 (32) end while (33) Output1: Set of Pareto optimal solutions (34) Output2: Set of Binary change detection masks Algorithm 1: Pseudocode of the Fuzzy MO-PSO algorithm for solving the change detection problem.particle.To achieve this, two objective functions are used to estimate fitness values of the particles and these objective functions are formulated as follows: where  is the degree of fuzziness ( ≤ 1),    is the degree of membership of This iteration will stop when where  is a termination criterion,  is the iteration number, and ‖ ⋅ ‖ 2 denotes the ℓ 2 -norm distance.
The proposed Fuzzy MO-PSO algorithm provides a set of Pareto optimal solutions and change detection masks so that it is necessary to select the most optimal solution amongst all the obtained solutions.It is worth noting that both objectives (e.g., changed and unchanged) are equally important to find the most optimum solution.Let  be the number of Pareto optimal solutions and all of them are first ranked based on their priority level.For instance, the first ranked solution is the best result for the unchanged fitness function and the last ranked optimum result is the best result for the changed fitness function.Therefore, the medium optimal solution of overall Pareto solutions is selected as a trade-off solution in the proposed approach.In this case, the medium solution is at /2 if the  is even number or at ( + 1)/2 if the  is an odd number.The Fuzzy MO-PSO algorithm has been applied to two different subband images; then the proposed algorithm generates two different binary masks.

Fusion Method.
In the change detection literature, there are several approaches to fuse multiple binary images to obtain the final change detection mask.For instance, Celik and Ma [23] used union based fusion approach to combine different binary subband images.However, using this approach will increase the error detection rates.In [43], MRF based approach is used to combine two binary image maps but this technique is not suitable in our case as the binary masks should be equivalent to using this method.In this paper, we propose a new fusion approach based on connected component and union techniques to fuse the LL and the HH binary masks for final binary mask estimation.To achieve this, the connected component technique based on 4-neighbour connectivity is first used to remove isolated pixels from the HH subband binary mask.More specifically, each changed pixel in the binary mask is compared with its two vertical and two horizontal neighbouring pixels.If any of these neighbouring pixels is flagged as changed (foreground) then the pixel remains as changed; otherwise the pixel is flagged as unchanged pixel (background).Note that noises are generally considered as high frequency components and appear in HH subband binary mask.Therefore, the connected component strategy is only applied to the HH subband binary mask.Finally, the union operator is applied to fuse the obtained HH and the LL subband binary masks in order to generate the final result.

Results and Discussion
3.1.Data Sets.In order to quantitatively and qualitatively understand and analyze the effectiveness and robustness of the proposed method, various semisynthetic images and two real multitemporal medium resolution Landsat data sets are used.The first data set belongs to the surface of the Lake Milh in Iraq and the second data set corresponds to the surface of the Caspian sea in Kazakhstan.Both data sets are provided by United States Geological Survey (USGS) [44,45].It is worth noting that the ground truth change masks for both data sets are not available.Therefore, for the first data set, the ground truth is generated by visual inspection and manual labelling and for the second data set, a region of interest which includes different atmospheric conditions but not land cover changes is considered to generate a ground truth.Even though the ground truth in the second data set does not reveal the exact differences between two Landsat images, one can use it to analyze sensitivity of the change detection methods with existence of atmospheric artifacts.In addition, we evaluate the change detection methods with quantitative and qualitative measurements on two different sets of semisynthetic images, generated using Lake Milh (Figure 3(a)).The first set of semisynthetic images are obtained by artificially inserting haze and different shapes and levels of thin cloud(s) into the various regions in Figure 3(a) (see Figure 5).The experiments in these data sets can reveal the accuracy of the change detection methods in different atmospheric conditions.To make the problem more challenging, the second set of semisynthetic images are generated by inserting nine different Gaussian noise levels into Figure 5 In this data set, the region of interest which is shown with blue rectangle is considered for quantitative evaluations.Moreover, the ground truth change mask for the region of interest (blue rectangular area) only includes unchanged pixels as there is not any land cover changes in the region of interest (see Figure 3(f)).

Quantitative Error Measurements.
For quantitative comparison purposes, three standard measurements including false alarm rate ( FA ), missed detection rate ( MD ), and total error ( TE ) [14] are used.The false alarm rate can be defined as the number of true unchanged pixels that were mistakenly chosen as changed pixels.This error measurement can be estimated in percentage as  FA = FA/ 0 × 100, where  0 is the overall number of unchanged pixels in the ground truth mask.The missed detection rate can be described as the number of true changed pixels that were erroneously selected as unchanged pixels.This error measurement can be calculated in percentage as  MD = MD/ 1 ×100, where  1 is the overall number of changed pixels in the ground truth mask.The third measurement is based on the total error which can be defined as the total of false alarm and the missed detection rates.The total error in percentage can be manipulated as  TE = (FA + MD)/(0+1)×100.Note that these three quantitative error measurements are used for the first real Landsat and semisynthetic data sets whereas only  FA is used for the second real Landsat data set as there is no land cover changes in the region of interest (region of interest in Figures 3(d) and 3(e)).

Change Detection Methods.
The proposed change detection method is compared with one frequency domain and three spatial domain based change detection methods which are: DT-CWT-based [23], PCA- means-based [14], ERGASbased [19], and PSO-GA-based [18].The first method is the DT-CWT-based [23] change detection method; it is applied into the frequency domain which uses three levels of image decomposition for generating subband images and then, the EM algorithm is used to find the binary mask.Note that, it is not necessary to set any parameters for the EM algorithm and the same number of decomposition level preferred in [23] is used.The PCA- means-based method has two different parameters which are nonoverlapping blocks ℎ and the dimensions of the eigenvector space  where we use ℎ = 4 and  = 2 in our experiments.The ERGAS-based method does not need any parameter to set.In PSO-GA-based method, we use the parameters given in [18].In the proposed method, the parameters, which must be tuned, are as follows: (1) bandwidth parameter (ℏ) of the mean shift clustering method, (2) -mean fuzziness parameter , (3) population size , (4) acceleration coefficients  1 and  2 , and (5) inertia weight .In our experiments, the bandwidth parameter (ℏ) is empirically chosen as 1.8.The rest of the parameters belong to the Fuzzy MO-PSO algorithm which are affecting the generation of population and the convergence behaviour.
The fuzziness parameter  in ( 8) controls flexibility of the Fuzzy MO-PSO, which affects the performance of the algorithm.Here, to find the optimum , we evaluate the performance of the proposed method based on TE in Figures 3(a) and 3(b) when  varies from 2 to 20.Note that the smaller the TE is, the better the change detection performance is.The effect of fuzziness parameter  is shown in Figure 4.It is clear that, when  ranges from two to five, the change detection performance is superior.For the other  values, the performance of the approach decreases significantly, so that based on this sensitivity analysis the -mean fuzziness parameter  is set to be 2. Note that the Fuzzy MO-PSO is a population optimization strategy and utilizes the probability transition instead of the certainty rule.With the increasing of particles' populations size, the probability of finding the global optimal solution is increasing.However, increasing the size of the particles' populations causes an increase in the computational time.Therefore, in the proposed method, the population size () is set to be 40 in order to provide the best trade-off between quality of solutions and computational time.
The next important parameter is the inertia weight  in (5)  searching speed of the global optimum result will be very slow.Because of all these reasons, these parameters affect the performance of the Fuzzy MO-PSO to resolve the change detection problem.In the proposed approach, the Fuzzy MO-PSO is used with inertia weight  = 0.4 and acceleration coefficients  1 = 2 and  2 = 2.These parameters are selected empirically and the proposed method with these parameter values provides satisfactory results in terms of obtaining the best change detection mask and computational cost.

Tests on Semisynthetic Data Sets
3.4.1.Data Sets with Haze and Thin Cloud.Generally, multispectral satellite images (e.g., Landsat images) are captured via a passive optical sensor, by those means this remote sensing technique is highly sensitive to any change in meteorologic conditions.Due to this sensitivity, Landsat images may include clouds.Furthermore, the contrast of the Landsat images acquired on the same geographical area at different time instances can be different because of the existence of haze.Existence of clouds and haze in the Landsat images is particularly a problem for land cover change analysis as they can reduce the contrast of the observed area and fade the colours in the image, so that they can influence the change detection accuracy by mistakenly detecting cloud as changed area and/or failing to correctly detect changed/unchanged pixels.Here we generate four different semisynthetic data sets with different percentages of cloud cover and thickness to quantitatively and qualitatively measure the performance of the change detection algorithms with existence of atmospheric artifacts.To this end, Figure 3(a) is firstly chosen and then haze and various percentages of thin cloud(s) with different thickness are artificially superimposed onto the different area of the image (see Figure 5).The PCA- means-based method shows better performance, but it still mostly detects cloud as changed areas.The PSO-GA-based method shows less sensitivity to haze and thin cloud as it uses a cost function based on CC similarity metric and a hybrid PSO-GA optimization method, which are resulting in less sensitivity to intensity variations.In other words, this method can correctly detect thin cloud as unchanged area when the cloud does not cause large variation of intensities.For instance, the results in Figures 6 and 7(d) show that when the cloud(s) is outside the lake area, the method correctly detects them as unchanged pixels; however, when the clouds are exactly at the top of the lake, the PSO-GA method fails to detect the cloud pixels as unchanged as there is a huge intensity variation between two Landsat images.The ERGAS-based method shows a good performance over existence of haze, but it removes image features which are under the thin cloud areas (see Figures 6 and 7(c)).In contrast to all these methods, the proposed strategy provides promising results and the results show the robustness of the proposed method to the existence of haze and thin cloud on images and removes the influence of these artifacts effectively.
The quantitative results for all the semisynthetic data sets, which are shown in Figure 5, are tabulated in Table 1.The quantitative results show that the proposed method provides the highest accuracy as compared with the state-of-the-art methods.The FA, MA, and TE measurement results verify that the proposed method is less sensitive to existence      show that the accuracy of the proposed method slightly decreases as the thickness of the cloud increases.For instance, the last column in Table 1 shows that the proposed method obtains the lowest  FA ,  MA , and  TE with 2.83%, 4.03%, and 3.56%, respectively.The DT-CWT based method provides the poorest performance with 40.23%, 12.65%, and 22.91% for the  FA ,  MA , and  TE , respectively.

Data Sets with Gaussian Noise.
During the Landsat satellite image acquisition as well as the transmission and reception process, noise can be superimposed onto the informational signal.Generally, the noise on the images can be seen as random variations in signal intensity, and as the level of the noise increases, uncertainty in the Landsat images becomes greater.One of the most common noise models in Landsat satellite images is the Gaussian noise.Since the Landsat images can include both noise and atmospheric artifacts (e.g., thin cloud and haze), it is important to evaluate performance of the change detection methods with respect to all these artifacts.To achieve this, Figure 5(b), which includes haze and various thin clouds in different regions of the Lake Milh (Figure 3(a)), is selected and nine different levels of zero mean Gaussian noise are superimposed onto the image.To quantitatively demonstrate the noise levels in the image, the Peak Signal-to-Noise Ratio (PSNR), in decibel (dB), is used.The PSNR can be defined as follows: 2 ) , (10) where   1 is the Landsat image with all the artifacts and  1 is the Landsat image which can be seen in Figure 5(b).Note that as the PSNR value decreases, the noise level increases.
To test the performance of change detection methods, firstly, various noisy synthetic images using Figure 5(b) are generated and these images are considered as the first input image ( 1 ).Moreover, the second input image is shown in Figure 3(b).Note that to generate noisy synthetic images, Gaussian noise model with a mean of zero and nine different standard deviations  is used.In our experimental setup,  varies in the range of 0.001 and 0.9 with the step size of 0.0112.
Various experiments on semisynthetic data (Figure 5(b)) with various levels of Gaussian noise are conducted to assess the change detection methods.The results are tabulated in Table 2 and shown in Figure 8. Figure 8(a) illustrates the image which is used to generate noisy images (  1 ) (Figure 8(d)).The second input image ( 2 ) is shown in Figure 8(b).Figure 8(c) shows the ground truth change mask, by which quantitative evaluations of binary change detection masks can be made.It is worth noting that   1 is corrupted by increasing the levels of the Gaussian noise varying from 50 dB (the least noise level) to 10 dB (the highest noise level).Figure 8(e) shows that DT-CWT-based method is robust to the existence of noise and mostly can remove this artifact correctly and this is mostly due to the implication of wavelet transforms in the change detection problem.The qualitative results of the PCA- means change detection method (Figure 8(f)) show that this method is highly sensitive to the existence of noise and can wrongly detect noises as changed pixels, so that as the noise level increases (PSNR less than 35 dB), the accuracy rapidly decreases.In contrast with these two methods, the ERGAS-based method provides higher rate of correct detection.Note that even with high levels of noise interference, ranging from 30 to 10 dB, the ERGAS-based method correctly flags noise and atmospheric   artifacts as unchanged pixels (see Figure 8(g)).However, due to the use of OTSU's thresholding technique to generate final binary change mask, many small changes are discarded in this technique.The PSO-GA-based method also shows very good performance against Gaussian noise and can detect small changes in the images, but those with great risibilities (see Figure 8(h)).This is due to the fact that this method employs an optimization technique to generate the final change detection mask rather than OTSU's thresholding technique.However, this method is sensitive to existence of atmospheric artifacts as it uses CC similarity metric.It is clear from Table 2 and Figure 8(i) that the proposed technique consistently outperforms the other approaches against different noise levels and atmospheric variations.This is mostly because the  proposed method uses hue colour channel with wavelet transforms and mean shift algorithm which makes the proposed method robust to illumination changes, existence of clouds, and Gaussian noise.Furthermore, the proposed method uses a new and effective fuzzy multiobjective approach to iteratively optimize two fuzzy fitness functions to obtain the optimal change detection mask.

Tests on Real Data Sets.
In order to understand and analyze the performances of the change detection methods for the real data sets, Figures 3(a) and 3(b) and Figures 3(d) and 3(e) are used.Figure 9 illustrates the qualitative results of the change detection methods using the real data set with the atmospheric artifacts shown in Figures 3(d) and 3(e).According to the qualitative results, it is clear that the DT-CWT-based method provides the least accurate results and it is very sensitive to variation of haze and existence of thin cloud in the images.The PCA- means-based method shows better performance, but it still wrongly detects cloud areas mostly as changed regions.The PSO-GA-based method shows less sensitivity to haze and thin cloud as it uses CC similarity metric which is less sensitive to intensity variations.The ERGAS-based method shows good performance over existence of haze, but it detects the dense thin clouds as changed and removes image features under the light thin cloud areas.In contrast to all these methods, the proposed strategy provides promising results and the results show the robustness of the proposed method to the existence of haze and thin cloud on images and removes the influence of these artifacts effectively.
The quantitative results for the real data sets, which are shown in Figure 3, are tabulated in Table 3.For the first real data set (Figures 3(a) and 3(b)), the ERGAS method provides the poorest performance and it obtains the lowest  FA ,  MA , and  TE with 9.38%, 18.83%, and 16.48%, respectively.The PSO-GA change detection method gives good results as it uses optimization strategy to find the final change detection mask and it obtains 4.65%, 2.88%, and 4.04% for  FA ,  MA , and  TE , respectively.However, the method is not successful in finding the most optimal result as it uses weighted sum of fitness functions.On the other hand, the proposed method gives the best accuracy results compared to the state-of-theart change detection methods with 3.36%, 2.52%, and 3.26% for  FA ,  MA , and  TE , respectively.In the second real data set (Figures 3(d) and 3(e)), the complexity increases in the real data set as the haze and thin cloud are involved in the images.Note that the region of interest with the blue rectangular region as shown in Figure 3 is considered to evaluate the quantitative results and only  FA is estimated since there is no any change in the region of interest.For this data set, the DT-CWT based and PCA- means methods provide the least performances with 82.51% and 68.33%, respectively.The ERGAS-based method provides good result because it is more robust to the haze and thin cloud than the other methods but it is not enough to obtain the best change detection mask.On the other hand, the proposed method gives the most accurate result as the haze and thin cloud mostly have influence on performances of the other methods.The combination of preprocessing and postprocessing (Fuzzy MO-PSO) strategies make the proposed method more effective and robust than the state-of-the-art methods to resolve the change detection problem even if the complexity increases with different kinds of artifacts (e.g., haze, thin cloud, and noise).
3.6.Discussion.Generally, there are several reasons that the state-of-the-art methods fail to resolve the change detection problem between two images efficiently: (1) The methods use the gray/colour value constancy assumption, which makes them not robust to the atmospheric artifacts.(2) The PSO-GA, the PCA- means, and the ERGAS-based methods are applied into the spatial domain which decreases the accuracy of detection near boundaries.(3) The PSO-GA method is a single-objective optimization solution approach and uses weighted sum of the objective functions which limits the searching of optimum result.(4) The DT-CWT based method uses simple union fusion based approach and neglects the low subband coefficients which results in increasing the incorrect detection rate.On the other hand, the proposed method uses different preprocessing and postprocessing strategies than the other methods which improve the performance and provide the greatest accurate change detection results.In addition, by comparing the results, we can conclude that the proposed method provides the best results as using the hue channel of Landsat images with the hybrid wavelet transform make the proposed method strongly robust to existence of artifacts.Another reason that the proposed method provides the highest performance over the sophisticated scenarios is due to the use of mean shift algorithm as it smoothes the intensity variations which results in better convergence rate and decreasing computational cost of the Fuzzy MO-PSO.Furthermore, the proposed MO-PSO provides the optimum trade-off between the fuzzy objective functions for better discrimination.
In this paper, we have discussed the performance of the change detection methods if one of the Landsat images has haze and thin clouds.However, if both images have thin clouds and haze, then the change detection problem becomes even harder to resolve, especially if these atmospheric artifacts are in the mutual regions of the two different images (worst-case-scenario).If both Landsat images have thin clouds and haze in different regions, the proposed method will successfully resolve the corresponding problem since it uses several preprocessing techniques before applying the Fuzzy MO-PSO.For the worst-case-scenario, the proposed method may still provide the better result compared to the existing methods but it may not be successful in obtaining the most optimum result.

Conclusion
In this paper, a new unsupervised approach in Landsat images using hybrid wavelet transform and mean shift clustering along with Fuzzy MO-PSO has been proposed to find the optimal change detection mask.In the method, a hybrid wavelet transform based on the DWT and the SWT has been used to decompose the hue channel of the Landsat images.Thus, the influence of atmospheric variations (e.g., thin cloud and haze) and noise artifacts have been significantly reduced in Landsat change detection problem.The LL and HH subbands have been used to generate two difference images and they are smoothed using mean shift algorithm.Next, the Fuzzy MO-PSO method is used to iteratively optimize fuzzy based fitness functions to find the Pareto optimal solutions and change detection masks for two different difference subband images.The most ideal solution is selected amongst many Pareto solutions as a trade-off solution for each subband image.Finally, two binary masks are combined to estimate the final binary image using the fusion method based on the connected component and union technique.The proposed method has been compared with the state-of-theart change detection algorithms for real and semisynthetic data sets and the results verify robustness and effectiveness of the proposed method to the existence of haze and thin cloud as well as variations of Gaussian noise in Landsat images.

Figure 1 :
Figure 1: Influence of an atmospheric correction technique on change detection.(a) The first image with existence of haze and clouds.(b) The result of the atmospheric correction method [28].(c) The second Landsat image.(d) The result of PCA--means method [14].
Photometric invariant technique Photometric invariant technique

Figure 2 :
Figure 2: Block diagram of the proposed method.
ℏ (, ) in the cluster ,    is the center of the cluster  in   ℏ , and (  ℏ (, ),    ) is the distance measurement which estimates the square error between the intensity value of the   ℏ at the pixel (, ) and the cluster center    .Fuzzy partitioning is carried out through MO-PSO of the objective functions shown in (4), with the update of membership    (, ) and the cluster centers    by    = ∑  =1 ∑  =1    (, )    (, ) ∑  =1 ∑  =1    (, )     (, ) = ∑ 1 =0 (  (, ) −    ) (b) (for instance, see Figure 8(d)).The experiments on semisynthetic images imitate realistic Landsat satellite acquisition situations and clearly show the pros and cons of the change detection methods.The first Landsat data set [44] is with 1984 × 1451 pixels and shows the water surface of the Lake Milh.The images are acquired on July 16, 1995 (Figure 3(a)), and May 11, 2003 (Figure 3(b)), by the Landsat Multispectral scanner to evaluate the changes in the lake.The ground truth change detection mask for the whole area of the lake is shown in Figure 3(c).The second Landsat data set [45] is with 541×335 pixels and illustrates the north east area of the Caspian sea in 2009 (Figure 3(d)) and 2012 (Figure 3(e)) with existence of thin cloud.

Figure 5 :
Figure 5: Semisynthetic images that depict the  1 image (Figure 3(a)) with haze and various thin clouds.

Figures 6 and 7
Figures 6 and 7 illustrate the qualitative results of the change detection methods using Figures 5(a) and 3(b) and Figures 5(b) and 3(b), respectively.According to the qualitative results, it is clear that the DT-CWT-based method provides the least accurate results and it is very sensitive to variation of haze and existence of thin cloud in the images.The PCA- means-based method shows better performance, but it still mostly detects cloud as changed areas.The PSO-GA-based method shows less sensitivity to haze and thin cloud as it uses a cost function based on CC similarity metric and a hybrid PSO-GA optimization method, which are resulting in less sensitivity to intensity variations.In other words, this method can correctly detect thin cloud as unchanged area when the cloud does not cause large variation of intensities.For instance, the results in Figures6 and 7(d)show that when the cloud(s) is outside the lake area, the method correctly detects them as unchanged pixels; however, when the clouds are exactly at the top of the lake, the PSO-GA method fails to detect the cloud pixels as unchanged as there is a huge intensity variation between two Landsat images.The ERGAS-based method shows a good performance over existence of haze, but it removes image features which are under the thin cloud areas (see Figures6 and 7(c)).In contrast to all these methods, the proposed strategy provides promising results and the results show the robustness of the proposed method to the existence of haze and thin cloud on images and removes the influence of these artifacts effectively.The quantitative results for all the semisynthetic data sets, which are shown in Figure5, are tabulated in Table1.The quantitative results show that the proposed method provides the highest accuracy as compared with the state-of-the-art methods.The FA, MA, and TE measurement results verify that the proposed method is less sensitive to existence

Figure 8 :
Figure 8: Evaluation of the various methods on semisynthetic data set with various noise levels.(a) Synthetic data set with the existence of haze and thin clouds; (b) Lake Milh data set ( 2 ) shown in Figure 3(b); (c) ground truth change mask.In (d), the image, shown in (a), is corrupted with different levels of zero mean Gaussian noise (from second row to fourth row, PSNR = 40, 20, and 10 dB, respectively).(e)-(i) shows the final change detection masks generated by DT-CWT, PCA- means, ERGAS, GA-PSO, and proposed method, respectively.
we say that  dominates  (denoted by  ≺ ) if   ≤   for  = 1, . . ., .We say that two vectors  and  are nondominated if the value at position  in the vector  is greater than the value at position  of vector    ≥   for  = 1, . . ., .
Definition 3.Vector  is Pareto optimal solution if it is nondominated.

Table 1 :
Quantitative results for semisynthetic data sets with haze and thin cloud.

Table 2 :
Quantitative results for semisynthetic data sets with haze, thin cloud, and Gaussian noise.

Table 3 :
Quantitative results for the real data sets.