Level set incorporated with an improved MRF model for unsupervised change detection for satellite images

ABSTRACT This study proposes the use of a level set incorporated with an improved Markov random field (MRF) model in unsupervised change detection for satellite images. MRF provides a means of modelling the spatial contextual information in the level set, and an edge indicator function is introduced into the MRF model to control the contribution of local information in the boundary areas to change detection. On the basis of the improved MRF model, local label relationships and edge information are considered in the level set energy functional to conduct a novel local term and attract the contours into desired objects. By merging the novel energy term, the proposed approach not only reduces noise but also obtains accurate outlines of the changed regions. Experimental results obtained with Landsat 7 Enhanced Thematic Mapper Plus and SPOT 5 data sets confirm the superiority of the proposed model when compared with state-of-the-art change detection methods.


Introduction
Change detection is the process of identifying changes in multitemporal satellite images that are acquired over the same geographical area at different points in time (Hussain, Chen, Cheng, Wei, & Stanley, 2013;Lu, Mausel, Brondízio, & Moran, 2004). Over the past few years, various change detection methods for different applications have been proposed. These methods can be divided into two kinds, namely, supervised methods and unsupervised methods (Bruzzone & Prieto, 2000). Unsupervised methods perform change detection by conducting a direct comparison of two multitemporal images without incorporating any additional information (Celik, 2009). The present study focuses on unsupervised change detection for multitemporal satellite images. Unsupervised change detection techniques are mainly developed on the basis of the so-called "difference images", which are generated using image differencing, image rationing, principal component analysis, or change vector analysis (CVA) (Lu et al., 2004). A few automatic analysis methods for identifying changes in difference images have been applied to unsupervised change detection; these methods include fuzzy c-means (Ghosh, Mishra, & Ghosh, 2011), Markov random field (MRF) (Wang et al., 2013), particle swarm optimization (Kusetogullari, Yavariabdi, & Celik, 2015) and support vector machine (Bovolo, Bruzzone, & Marconcini, 2008).
In recent years, the active contour model (ACM) based on level set methods has gained popularity in the field of change detection for satellite images because it can achieve robust image segmentation and naturally handle topological changes (Bazi, Melgani, & Al-Sharari, 2010;Li, Gong, & Liu, 2015). The level set method was first introduced into the ACM by Osher and Sethian (1988) to represent a contour as the zero level set of a function, and the level set function evolved according to a partial differential equation. Chan and Vese (2001) improved this model as the Chan-Vese (CV) model based on the Mumford-Shah energy functional (Mumford & Shah, 1989). The overall premise of the CV model is to avoid using gradient information and consequently reduce the complexity of optimization problems. On the basis of the CV model, Bazi et al. (2010) proposed multiresolution level set (MLS) methods for detecting changes using multispectral satellite images. In another study, wavelet transformation was adopted to obtain a multiresolution representation of the difference image and in turn reduce the effect of noise contamination; then, the ACM was applied to the transformed images for change detection (Celik & Ma, 2011). Ardila, Bijker, Tolpekin, and Stein (2012) introduced ACMs with a localized data fitting energy term for identifying tree-crown image regions in high resolution multitemporal images. A fuzzy ACM and genetic algorithm were used in a previous study to analyse the difference image and reduce speckle noise for change detection in synthetic aperture radar images (Shi, Wu, Paul, Jiao, & Gong, 2014). Hao, Shi, Zhang, and Li (2014) developed an expectation-maximization (EM)-based level set method for change detection, which did not require initial contours. Li, Shi, Myint, Lu, and Wang (2016) proposed a semi-automated change detection approach based on a fast level set algorithm for landslide inventory mapping. Normally, level set methods based on the CV model assume that the label of each spatial location in an image is independent of others, thereby ignoring neighbouring label relationships among image pixels during the evolutions of the level set curves. Thus, change detection results may contain much salt-and-pepper noise or lose detailed changes such as accurate outlines. Within this context, the present work proposes a level set algorithm incorporating an improved MRF model for unsupervised change detection from satellite images. In order to generate accurate change detection results, a novel local energy term considering neighbouring label relationships and edge information is introduced into the energy functional on the basis of the improved MRF model.
The flowchart of our algorithm is shown in Figure 1. First, the CVA method is conducted on the multitemporal satellite images to generate a difference image. Second, image downsampling is performed to generate the multi-scale representation of the difference image. Next, the initial curves of the level set are located and the level set function is initialized as a signed distance function. The level curves can then evolve with the MRF-based local energy constraints in the multi-scale framework. Finally, the change map can be obtained by the improved level set algorithm.

Proposed change detection method
Consider two coregistered multispectral satellite images,

Change detection by level set
Change detection using level set methods aims to find an optimal contour C, which splits the difference image I : Ω ! < into non-overlapping regions associated with changed classesw 1 and unchanged classes w 2 . In the Bayesian framework, the optimal segmentation is given by the maximum a posteriori estimate, which is described as follows: (1) where L ¼ ðl 1 ; l 2 ; . . . ; l n Þ is a label set of the difference image and n is the pixel number. PðIÞ is a prior probability distribution of the class labels of the difference image and PðIjLÞ is the probability density function of the pixel values in the difference image.
In the variational framework, maximizing the posterior probability in Equation (1) is equivalent to minimizing the following energy functional (Brox & Weickert, 2006): where ϕ is the level set function, such that ϕðx; yÞ > 0 if ðx; yÞ 2 Ω 1 and ϕðx; yÞ < 0 if ðx; yÞ 2 Ω 2 ; H denotes the Heaviside step function, i.e. HðzÞ ¼ 1 if z ! 0 and HðzÞ ¼ 0 if otherwise. Pðw i Þ is a prior probability of class w i (i ¼ 1; 2) and pðIjw i Þ represents the probability density function, given class w i (i ¼ 1; 2). The parameter μ is a constant that controls the trade-off between the goodness of fit and the length of the contour. The minimization of the energy function is achieved by solving the associated Euler-Lagrange equation. The evolving equation of level set curves is given by where δðϕÞ denotes the Dirac function. The CV model assumes that the image is a piecewise constant function (Chan & Vese, 2001). If the prior probabilities Pðw 1 Þ and Pðw 2 Þ of both regions are considered equal, the associated evolving equation can be obtained as follows: where c 1 and c 2 are the intensity averages of image I inside and outside the contour C, respectively.

Level set incorporating MRF
In this study, the Markov prior information was incorporated into CV model to reduce spurious noisy regions. MRF is an effective approach to exploit the spatial contextual information within a Bayesian framework (Bruzzone & Prieto, 2000;Szeliski et al., 2008). Using the conditional distribution of pixel labels based on the Markov property, the prior probability distribution of class labels can be expressed as where N x denotes the local neighborhood of the pixel x (x‚N x ), Uðl i Þ is the Gibbs energy function and Z is a normalizing factor. In the MRF model, if the boundary pixels are not properly controlled, the result will reveal over-smooth boundaries and lose significant details (Tso & Olsen, 2005). Thus, in this study, the edge information is used to control the contribution of the spatial information in the boundary areas to change detection. In addition, the Potts model is utilized to model each of the conditional distributions of pixel labels in MRF. Then, Uðx i Þ is given as follows: where β is a constant that controls the influence of the spatial contextual information in the MRF model. g Represents an edge indicator function, and gðxÞ ¼ 1 if the pixel x is recognized as an edge; otherwise, it is equal to 0. The function V is given by Equation (7): By merging the prior probability distributions of the two class labels expressed as Equation (5), the evolving formulation of Equation (3) can be written as follows: Then, the final evolving formulation can be expressed as follows: where the term À I À c 1 j 2 þ I À c 2 j 2 is the global data fitting term derived from the CV model and the term ωðÀUðw 1 Þ þ Uðw 2 ÞÞ refers to the proposed MRF-based local penalty term. ω is set to c 1 À c 2 j j, which is calculated adaptively according to the intensity averages of image I inside and outside the contour during the evolution of level curves. When the value of ω is very small, which means the intensity averages of the changed and unchanged class are close to each other, the global term will dominate to attract the contours to the objects rapidly. As the value of ω increases, the local term will tend to control the neighbouring label relationships and generate an accurate change detection result.

Implementation
Multi-scale analysis was performed to reduce the search space and avoid being trapped into a local minimum (Bazi et al., 2010). The multi-scale images were generated by downsampling the difference image. The proposed model was then applied to the images with a coarse-to-fine evolving strategy. The main steps of the proposed level set algorithm are summarized as follows: (1) Decompose the difference image into K scales.
(2) Let k ¼ K and initialize the level set function ϕ.

Experimental results and discussion
To assess the performance of the proposed model for change detection, two satellite image data sets were tested in the experiments. Three standard error measures were used to evaluate the results (Yetgin, 2012): (1) missed alarm rate (P m ): the number of changed pixels in the changed detection map that were incorrectly classified as unchanged pixels over the total number of changed pixels in the ground reference map; (2) false alarm rate (P f ): the number of unchanged pixels in the changed detection map that were incorrectly classified as changed pixels over the total number of unchanged pixels in the ground reference map; and (3) total error rate (P t ): the total number of detection errors caused by either missed or false detections over the total number of pixels.
In the experiments, the initial curves of the level set were small circles that evenly covered the entire difference image. The optimal Canny edge detector (Canny, 1986) was used in the proposed model to extract the edge information. We compare the proposed model with the MLS method, CV model, undecimated discrete wavelet transform and active contour (UDWTAC) approach (Celik & Ma, 2011), MRF model (Bruzzone & Prieto, 2000) and EM-based threshold method to prove the effectiveness of the proposed model.

Experiment on data set 1
The first data set consisted of two satellite images acquired by the Landsat 7 Enhanced Thematic Mapper Plus (ETM+) in August 2001 and August 2002 over the Liaoning Province of China. Both images were coregistered, and the histogram matching method was implemented for the relative radiometric correction. Figure 2(a,b) shows band 4 of the 2001 and 2002 images, respectively. The ground reference map was created manually on the basis of the analysis of multitemporal input images. This map is shown in Figure 2(c). Figure 3 shows the change detection results generated by the MLS, CV, UDWTAC, MRF, EM and the proposed model. In the proposed model, we set the parameter of the level set μ ¼ 0:1, the time step Δt = 0.1 and β = 3.2. In addition, the β value of the MRF method is set to 0.8, and the CV model, MLS and UDWTAC are implemented with the same parameters of μ = 0.4 and Δt = 0.1. The visual comparison between the different changed results confirms that the proposed model yields the closest change map to the ground reference. As marked by the rectangles in Figure 3, the noise is reduced, and accurate outlines of the changed regions are detected in the proposed model. Figure 3(a-c) shows that the results of the other level set methods contain significant noise and lost certain detailed change information. Although the MRF and EM yield homogenous regions more effectively than the other level set methods, the former generates over-smooth results on the boundary areas, and the latter produces much "salt-and-pepper" noise, as shown in Figure 3(d,e), respectively. Figure 4 presents the variations in the missed detection rates, false alarm rates and total error rates when different β values are used in the proposed model. The change maps are generated using β values ranging from 0 to 3.6 with a step of 0.2. As the β value increases, the miss detection rates continuously decrease, while the false alarm rates first slightly drop and then increase. As a result, the total error rates decline and reach the minimum at the β value of 3.2. The change detection result with β ¼ 0 contains more error pixels than the result with the β value larger than 0, which reveals that the MRF-based local energy can improve the accuracy of change detection.
To quantify the effectiveness of the proposed model, we compare the computations of P m , P f and P t based on the change maps and the reference map, as depicted in Table 1. The values of P m , P f and P t obtained by the proposed model are 21.23%, 0.58% and 4.62%, respectively. Thus, the proposed model can obtain more accurate change detection results than the other methods in this experiment.

Experiment on data set 2
The second data set represented two 2.5 m spatial resolution images acquired by the SPOT 5 satellite over the Wuqing district in Tianjin, China, in 2008 (t1, shown in Figure 5(a)) and 2009 (t2, shown in Figure 5(b)). The t1 image was registered to the t2 image, and the histogram matching method was used in the relative radiometric correction. Figure 5(c) displays the ground reference map that was collected manually by comparing the multitemporal images.
In the proposed model, we set the parameter of the level set μ ¼ 0:1, the time step Δt = 0.1 and β = 2. The MLS, CV and UDWTAC are also implemented to generate change detection maps with the same parameters μ= 0.4 and Δt = 0.1. Furthermore, the β value of the MRF method is set to 0.8 for change detection. All change maps are displayed in Figure 6. From a qualitative point of view, the change map provided by the proposed model appears to be closer to the ground reference, compared with change detection maps. The rectangles in Figure 6 show that MLS, CV and UDWTAC yield several false alarms and generate certain noise in large change blocks. By contrast, the proposed model reduces the false alarms and retains detailed changes in the homogenous regions, as shown in Figure 6(f).
The variations in the three indices (i.e. missed detection rate, false alarm rate and total error rate) with different β values in the proposed model are displayed in Figure 7. Different change maps produced by β values ranging from 0 to 3.6 with a step of 0.2 are used to analyse the effects of different β values on change detection results. The false alarm rates decrease with an increase in the β value, while the miss detection rates first remain nearly constant and then increase. The total error rates decline and reach the minimum when the change map is generated by the β value of 2. As shown in Figure 7(c), the change detection method generates the largest number of error pixels when the β value equals 0, which proves the effectiveness of the local energy term. Table 2 illustrates the quantitative error measures obtained by all the change detection methods used in this research. For the proposed model, the P m , P f and P t values calculated by comparing the generated change map with the ground reference map are equal to 6.97%, 8.24% and 8.17%, respectively. The proposed model generates the most accurate change detection result among all the methods in the experiment.

Discussion
Since the proposed level set model incorporates MRF-based local information to fully describe the mutual relationships of pixel labels in the spatial neighbourhood system, the noise in the change detection map can be reduced. In addition, an edge indicator function is used to avoid the excessive use of local information in the detailed areas and the proposed model can consequently retain accurate outlines of the changed regions. The experimental results confirm that the accuracy of the change detection result is adequately improved by the proposed model.
The penalty coefficient β in MRF controls the influence of the local term in the level set evolution functional. For the change detection result with the elimination of the local penalty term in the level set evolution, namely, β ¼ 0, there exist more error     pixels than the result with a larger β value. This indicates that the MRF-based local energy term can guide the level curve evolution correctly and improve the accuracy of change detection.

Conclusion
This study develops a novel change detection approach for satellite images using a level set method incorporated with an improved MRF model. Both local information considering neighbouring label relationships and edge information are introduced into the level set energy functional to conduct a local penalty term. Due to the introduction of the novel MRF-based term, the proposed model can reduce noise and while preserving the outlines in boundary areas. The experiments conducted on the multitemporal Landsat 7 ETM+ and SPOT 5 data sets confirm the superiority of the proposed approach in generating accurate change detection results when quantitatively and qualitatively compared with state-of-the-art change detection methods. Apart from spatial contextual information, other information, such as texture information and morphology characteristics, can be used in this method to improve the accuracy of change detection in future work.