SAR change detection based on intensity and texture changes

ISPRS Journal of Photogrammetry and Remote Sensing 93 (2014) 123–135

Contents lists available at ScienceDirect

ISPRS Journal of Photogrammetry and Remote Sensing journal homepage: www.elsevier.com/locate/isprsjprs

SAR change detection based on intensity and texture changes Maoguo Gong ⇑, Yu Li, Licheng Jiao, Meng Jia, Linzhi Su Key Laboratory of Intelligent Perception and Image Understanding of Ministry of Education, Xidian University, Xi’an 710071, PR China

a r t i c l e

i n f o

Article history: Received 8 December 2013 Received in revised form 1 April 2014 Accepted 10 April 2014 Available online 9 May 2014 Keywords: Change detection Multivariate generalized Gaussian model Robust principal component analysis Graph cuts Synthetic aperture radar

a b s t r a c t In this paper, a novel change detection approach is proposed for multitemporal synthetic aperture radar (SAR) images. The approach is based on two difference images, which are constructed through intensity and texture information, respectively. In the extraction of the texture differences, robust principal component analysis technique is used to separate irrelevant and noisy elements from Gabor responses. Then graph cuts are improved by a novel energy function based on multivariate generalized Gaussian model for more accurately fitting. The effectiveness of the proposed method is proved by the experiment results obtained on several real SAR images data sets. Ó 2014 International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS) Published by Elsevier B.V. All rights reserved.

1. Introduction The task of change detection is to analyze two registered remote sensing images acquired in the same geographical area at two different times (Bazi et al., 2005; Blaschke, 2010). With the development of earth observation technology, there is much multitemporal synthetic aperture radar (SAR) data available for monitoring applications (Hussain et al., 2013; Oliver and Quegan, 2004). However, the presence of the speckle noise in SAR images (Schubert et al., 2013; Villasensor et al., 1993) makes the related studies difficult. Generally, there are two pivotal steps in the process of unsupervised change detection (Bruzzone and Prieto, 2002). One is to generate a difference image; the other is to use effective methods for analyzing the difference image and identifying the pixels as changed or not. Both steps can affect the final result. For the first step, a difference image can be generated by applying the subtraction operator or the ratio operator to the intensity of the two multitemporal images pixel by pixel (Celik, 2010). The ratio operator is widely used because it is more adaptable to the statistics of SAR images and robust to radiometric errors (Rignot and van Zyl, 1993). Additionally, the ratio image is expressed in a logarithmic scale. But the effectiveness of such a pixel-based operator to suppress speckle noise is limited, because it ignores the spatial information. When changes occur between two multitemporal SAR images on their corresponding land cover, only the changes in ⇑ Corresponding author. Tel.: +86 029 8820 2661. E-mail address: [email protected] (M. Gong).

gray cannot lead to an accurate result. Rignot proposed the mean ratio detector by considering the changes in first-order statistics (Rignot and van Zyl, 1993). Although it has made use of local information, the method may fail to detect changes. The reason for this is that when changes occur at the texture level, the mean value may still be preserved. For SAR images, texture, which is a representation of the spatial relationship of gray levels in an image, provides important characteristics for surface and object identification (He et al., 2010). It can reveal the intrinsic spatial variability of SAR images, which is a valuable feature in discriminating the changes between land-cover types. As mentioned in (Ulaby et al., 1986), different surfaces can be characterized by distinct textural features, for example, rocks, sea-ice, sea-water, vegetation as well as urban areas. If water could be seen as textureless (Ulaby et al., 1986), the vegetation category belongs to medium texture class, and the urban category represents a high texture class. Thus, when one type of land cover changes to another, the texture difference will be reflected. Therefore we can compare the changes between the textures of the two images directly. For the above mentioned reasons, in this paper, textural feature is selected to generate the texture difference image, which can compensate for the deficiency of intensity difference image. For the second mentioned step, its aim is to identify the changes between the two co-registered SAR images. In the following paragraphs, some of the most popular methods proposed in the past decades will be reviewed. In 2000, Bruzzone proposed an automatic threshold technique based on the Bayes theory (Bruzzone and Prieto, 2000). It obtains

http://dx.doi.org/10.1016/j.isprsjprs.2014.04.010 0924-2716/Ó 2014 International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS) Published by Elsevier B.V. All rights reserved.

124

M. Gong et al. / ISPRS Journal of Photogrammetry and Remote Sensing 93 (2014) 123–135

Image I1

Image I2

IDI

TDI

MGG-based graph cuts

Change-Detection Map Fig. 1. General scheme of the proposed approach.

the optimal threshold using Expectation–Maximization (EM) iterative algorithm under two assumptions. One is that the pixels of difference image are spatially independent. The other is that the conditional density functions of the two classes (changed and unchanged) can be modeled by Gaussian distributions. Considering the spatial–contextual information of each pixel, Bruzzone used a MRF approach to improve the accuracy (Bruzzone and Prieto, 2000). In 2005, Bazi et al. produced the change-detection map based on Kittler–Illingworth (K&I) threshold selection criterion under the generalized Gaussian (GG) assumption (Bazi et al.,

2005). Although they achieve good results for lots of images, when changed and unchanged classes are strongly overlapped or their statistical distribution cannot be modeled accurately, thresholding methods provide poor results. So the main disadvantage of these automatic methods is the difficulty in modeling changes and unchanges in the difference image correctly. Change detection can also be seen as a segmentation issue to discriminate the pixels in the difference image between changed and unchanged classes. So we can refer to segmentation approaches and adapt them to our studies. For example, Bazi et al. used level set approach in change detection (Bazi et al., 2010); Celik et al. proposed a technique for change detection in optical satellite images using principal component analysis and the k-means clustering algorithm (it is also attractive for SAR images). In 2012, our group used an enhanced fuzzy clustering in change detection and gained good performance (Gong et al., 2012). The advantage of these methods is that they do not need to establish statistical models. However, these algorithms do have some limitations. For example, each algorithm has its own bias for the optimization of different criteria. In addition, they are sensitive to their initializations which may induce a tendency to get stuck to local optima. Besides the above mentioned methods, recent graph cuts attract more and more attention in the field of computer vision. It is an energy-based segmentation technique, which aims to minimize an energy function by seeking the minimum cut in a weighted graph. By utilizing both boundary and regional information, it can guarantee continuity and lead to smoothness of the contours. Furthermore, it can find the global minimum. In this paper, we

Fig. 2. Multitemporal images relating to Ottawa. (a) Image acquired in July 1997 during the summer flooding. (b) Image acquired in August 1997 after the summer flooding.

Table 1 Change detection results of the Ottawa data set obtained by different methods. The best value for each metric is shown in bold. Method

Difference images

False-alarm rate (%)

Miss-alarm rate (%)

PCC

Kappa

K&I

IDI TDI ITDIs

3.05 0.63 1.11

0.95 1.44 0.98

0.9600 0.9793 0.9791

0.8575 0.9208 0.9219

FCM

IDI TDI ITDIs

0.20 0.28 0.14

2.02 1.56 1.97

0.9778 0.9815 0.9789

0.9125 0.9283 0.9167

G_GC

IDI TDI ITDIs

1.77 0.42 1.75

0.61 1.83 0.56

0.9762 0.9776 0.9769

0.9138 0.9126 0.9157

MGG_GC

ITDIs

1.21

0.59

0.9820

0.9335

M. Gong et al. / ISPRS Journal of Photogrammetry and Remote Sensing 93 (2014) 123–135

make it adapt to our change detection studies and prove its effectiveness by experiments. Graph cuts were first proposed by Boykov and Jolly (Boykov et al., 2001) in 2001. Its energy function consists of two components, data constraint and smoothness constraint, respectively. Data constraint is designed as a penalty for assigning the corresponding label to a pixel, and it is computed by using the negative log-likelihood mostly. To get the value of this term, it is often assumed that the data follow the Gaussian or Multivariate

Gaussian model for different purposes (Pelizzari and BioucasDias, 2007a, 2007b). However, in practice these models assumption cannot fit with the real data distribution, especially for SAR images (Bazi et al., 2005). The multivariate generalized Gaussian model (MGG) has the ability to approximate a gigantic class of statistical distributions, and a more exact probability density function can be obtained by adjusting the shape parameter. Thus, in this paper, MGG is proposed to model the two classes (changed and unchanged) for more accurately fitting.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

125

(k)

Fig. 3. Final change-detection maps obtained by different approaches relating to the Ottawa data set. (a) K&I based on IDI. (b) K&I based on TDI. (c) K&I based on ITDI. (d) FCM based on IDI. (e) FCM based on TDI. (f) FCM based on ITDIs. (g) G_GC based on IDI. (h) G_GC based on TDI. (i) G_GC based on ITDIs. (j) MGG_GC based on ITDIs. (k) Ground truth.

126

M. Gong et al. / ISPRS Journal of Photogrammetry and Remote Sensing 93 (2014) 123–135

(a)

(b)

Fig. 4. Multitemporal images relating to the city of Bern. (a) Image acquired in April 1999 before the flooding. (b) Image acquired in May 1999 after the flooding.

Table 2 Change detection results of the Bern data set obtained by different methods. The best value for each metric is shown in bold. Method

Difference images

False-alarm rate (%)

Miss-alarm rate (%)

PCC

Kappa

K&I

IDI TDI ITDIs

1.78 0.02 0.44

0.10 0.43 0.11

0.9812 0.9956 0.9945

0.5482 0.7906 0.8064

FCM

IDI TDI ITDIs

0.47 0.14 0.22

0.33 0.21 0.15

0.9920 0.9965 0.9963

0.7000 0.8591 0.8556

G_GC

IDI TDI ITDIs

0.66 0.18 0.22

0.06 0.18 0.15

0.9929 0.9964 0.9963

0.7700 0.8561 0.8562

MGG_GC

ITDIs

0.16

0.17

0.9967

0.8705

This work contributes a new method of constructing IntensityTexture Difference Images (ITDIs), including Intensity Difference Image (IDI) and Texture Difference Image (TDI). ITDIs are two difference images, which can exploit both gray and texture changes. To get the change-detection map based on ITDIs, the graph cuts technique is adopted. And in formulating its energy function, MGG is used first as a more precise statistical model. This paper is organized as follows: Section 2 describes the new method for change detection with ITDIs in details. Section 3 presents experimental results, where different data sets are used to demonstrate the effectiveness. Finally, the discussion and conclusion are presented in Section 4. 2. Methodology Consider two co-registered intensity SAR images acquired over the same geographical area at two different times. They are expressed as I1 and I2 of size H  W. X = {xc, xu} is the set of classes corresponding to changed and unchanged. In order to reduce speckle noise, many well-known adaptive filters can be used before generating the difference images, e.g., the Lee, the Kuan, the Frost, and the Gamma-MAP filter (Bovolo and Bruzzone, 2005). In this paper, the Lee filter is applied to the SAR images before generating the difference images. Then, as shown in Fig. 1, the proposed method has two main steps: (1) Generation of ITDIs; (2) Automatic analysis of ITDIs. 2.1. Generation of ITDIs In this work, two difference images are used. One is IDI, usually generated by applying the log-ratio operator on SAR images I1 and

I2. The log-ratio operator is widely used in literatures. IDI can be obtained by the following function

    I1  IDI ¼ Log : I2 

ð1Þ

The other is TDI, which presents texture difference and can make use of the information about the neighborhood of each pixel to compensate the deficiency of IDI. The following section will focus on how to construct TDI in details.

2.1.1. Generating TDI by Gabor filters Texture provides important characteristics for surface and object identification (Min et al., 2010; Torres-Torriti and Jouan, 2001; Li and He, 2006). When changes occur in the land covers of the two multitemporal SAR images, the textures also change. Therefore, texture can be utilized to improve change detection reliability. There are many different texture analysis techniques, such as Co-occurrence matrix (GLCM) (Kandaswamy et al., 2005b), Fourier power spectrum (Florindo and Bruno, 2012), Markov random field (Deng and Clausi, 2004, 2005), and Gabor filters (Jain and Farrokhnia, 1991). Among them, Gabor filter banks have strong discriminating power in classification problems (Kandaswamy et al., 2005a; Manjunath and Ma, 1996). Some researchers have made a comparison of Gabor and other texture features (Manian and Vasquez, 1997; Torres-Torriti and Jouan, 2001). The results indicate that Gabor filters can achieve comparable or better performance at very low computation cost. The Gabor representation of an image can be obtained by convolving the image with a two-dimensional Gabor kernel. Let G(x,

M. Gong et al. / ISPRS Journal of Photogrammetry and Remote Sensing 93 (2014) 123–135

y) be Gabor kernel (Kruizinga and Petkov, 1999), and the convolutions of IDI with Gabor kernel are given by

GDIðx; yÞ ¼ IDIðx; yÞ  Gðx; yÞ:

ð2Þ

With Nk scales and Nh orientations, we will get N sets of Gabor responses of IDI(x, y), where N = NkNh. When each Gabor GDIi is stacked to one column, a Gabor difference image matrix is achieved by

M GDI ¼ ðGDI1 ; GDI2 ;    ; GDIi ;    ; GDIN Þ:

ð3Þ

127

2.1.2. Separating the noises in TDI Using MGDI to represent texture changes is extremely redundant. For example, if 5 scales and 12 orientations of Gabor filters are selected, each pixel in the image is now transformed to a 60 dimensional vector. Moreover, in many cases some irrelevant or noisy elements are involved in Gabor responses (Kim and Hong, 2009), which can lead to a poor change detection result. We can benefit from the advantage of large scale Gabor responses to suppress speckle noise, because large scale makes the similar areas affected by speckle noise show similar texture. However, artifact

Fig. 5. Final change-detection maps obtained by different approaches relating to the Bern data set. (a) K&I based on IDI. (b) K&I based on TDI. (c) K&I based on ITDI. (d) FCM based on IDI. (e) FCM based on TDI. (f) FCM based on ITDIs. (g) G_GC based on IDI. (h) G_GC based on TDI. (i) G_GC based on ITDIs. (j) MGG_GC based on ITDIs. (k) Ground truth.

128

M. Gong et al. / ISPRS Journal of Photogrammetry and Remote Sensing 93 (2014) 123–135

Fig. 6. Multitemporal images relating to the coastline of the Yellow River Estuary. (a) Image acquired in June 2008. (b) Image acquired in June 2009.

Fig. 7. Multitemporal images relating to the farmland area of the Yellow River Estuary. (a) Image acquired in June 2008 before the planting. (b) Image acquired in June 2009 after the planting.

Fig. 8. Multitemporal images relating to the river of the Yellow River Estuary. (a) Image acquired in June 2008. (b) Image acquired in June 2009.

errors may occur at the same time. So selecting the optimal features is important but also difficult. Usually, it is based on a trial and error method or time-consuming selection process (Jain and Karu, 1996; Ma and Manjunath, 2000).

For the motivation of dimension reduction and error features removal, we can refer to the low-rank subspace theory (Candès et al., 2011; Candès and Recht, 2009; Keshavan et al., 2010). In 2009, Wright et al. proposed Robust Principal Component Analysis

129

M. Gong et al. / ISPRS Journal of Photogrammetry and Remote Sensing 93 (2014) 123–135 Table 3 Change detection results of the coastline area of the Yellow River Estuary obtained by different methods. The best value for each metric is shown in bold. Method

Difference images

False-alarm rate (%)

Miss-alarm rate (%)

PCC

Kappa

K&I

IDI TDI ITDIs

2.41 0.02 0.19

0.01 0.34 0.21

0.9758 0.9964 0.9960

0.6263 0.9065 0.9027

FCM

IDI TDI ITDIs

1.56 0.65 0.59

0.02 0.00 0.01

0.9842 0.9935 0.9940

0.7207 0.8654 0.8727

G_GC

IDI TDI ITDIs

0.41 0.34 0.32

0.02 0.10 0.01

0.9957 0.9956 0.9967

0.9052 0.9004 0.9262

MGG_GC

ITDIs

0.12

0.12

0.9977

0.9442

Fig. 9. Final change-detection maps obtained by different approaches relating to the coastline of the Yellow River Estuary. (a) K&I based on IDI. (b) K&I based on TDI. (c) K&I based on ITDI. (d) FCM based on IDI. (e) FCM based on TDI. (f) FCM based on ITDIs. (g) G_GC based on IDI. (h) G_GC based on TDI. (i) G_GC based on ITDIs. (j) MGG_GC based on ITDIs. (k) Ground truth.

130

M. Gong et al. / ISPRS Journal of Photogrammetry and Remote Sensing 93 (2014) 123–135

(RPCA) (Wright et al., 2009) aiming to recover a low-rank matrix L from highly corrupted measurements M = L + S, where S represents sparse errors with arbitrarily large magnitudes. P Let kLk :¼ i ri ðLÞ, which denotes the nuclear norm of matrix L. P It is a convex relaxation of rank function. And kSk1 = ij|Sij| denotes the l1-norm of S. The RPCA aims at decomposing M into the low-rank matrix L and the sparse S by minimizing

kLk þ kkSk1

s:t: M ¼ L þ S;

ð4Þ

where parameter k is used to balance the effects of the two parts. Gabor difference images (GDIi, i = 1, 2  N) are elements in RN (N is the number of filters). These Gabor difference images have similar structure. So we assume that they come from the same subspace. RPCA can decompose MGDI into two components. One is the low-rank part of MGDI, which presents the good features recovered from GDIs, named as LGDI. The other is the sparse part of MGDI denoting the irrelative and error features that we want to get rid of, named as SGDI. Although LGDI can capture the spatial characteristics of textured regions and contains more effective information, the size of LGDI is as same as that of the MGDI. To make it applicable to the modelbased segmentation method and avoid the problems caused by high dimension, we map LGDI to one dimension. Finally, the stack column formation of TDI is obtained by

TDISC ¼

N X 2 jLGDI ð:; iÞj

Each pixel p e P in ITDIs has two features, which represent intensity difference and texture difference, respectively. The goal of change detection is to assign p to xc (changed class) or xu (unchanged class). Graph cuts are used here. Since the composite feature vectors consist of the intensity feature and the texture feature, the unsupervised graph cuts should be based on multivariate finite mixture model. The process of graph cuts is to minimize an energy function (Boykov and FunkaLea, 2006; Boykov and Kolmogorov, 2004; Kolmogorov and Zabin, 2004). The energy function consists of two terms, and it is defined by

X X Dp ðlp Þ þ V p;q ðlp ; lq Þ  dlp –lq p2P

ð6Þ

p;q2N



1 if lp –lq 0 if lp ¼ lq

ð7Þ

:

Cðn=2Þbl

eF ;

ð9Þ

ib l 1h ðzp  ml ÞT R1 ; l ðzp  ml Þ 2

ð10Þ

Prðzp jhl Þ ¼

F¼

2.2. Automatic analysis ITDIs

ð8Þ

where zp is a two dimension feature vector (zIp ; zTp ) of pixel p, which are the values of IDI and TDI, respectively. Pr (zp|hl) is the probability of pixel p belonging to l, and hl is a set of parameters of the distribution model associating with label l. These composite feature vectors are assumed to be generated from a multivariate finite mixture density model. Among the possible models, GG distribution is a particularly attractive candidate (Bazi et al., 2005). However, a direct multivariate extension of univariate GG distribution does not exist. The MGG used here is a particular case of the Kotz-type distribution defined as (Verdoolaege and Scheunders, 2011)

ð5Þ

:

LGDI(:, i) is the ith column of LGDI. After reshaping TDISC into the original size of the input images, we get TDI.

where dlp –lq ¼

Dp ðlp Þ /  ln Prðzp jhl Þ;

!1=2

i¼1

EðlÞ ¼

Dp(lp) is data constraint term and Vp,q(lp, lq) is smoothness constraint term. Dp(lp) derives from the observed data and measures the cost of assigning the label l (changed or unchanged) to the pixel p. To obtain more accurate object boundaries, Vp,q(lp, lq) specifies the penalty of the discontinuity between p and q of similar attributes in N , where N  P  P is the neighborhood of pixel p. The value of the individual Dp(lp), which is treated as the penalty of assigning pixel p to label l, is computed using the negative log-likelihood

n 2bl

2

p Cðn=2bl ÞjRl j1=2 n 2

where n is the dimensionality of the probability space, and here we set n = 2. The terms ml, Rl, and bl are the mean vector, the covariance matrix, and the shape parameters, respectively. jRl j is determinant of Rl, and CðÞ is Gamma function which is defined as R1 CðxÞ ¼ 0 et tx1 dt. The shape parameter bl determines the decay rate of the density. Note bl = 1 yields the multivariate Gaussian density function and bl = 1/2 yields the multivariate Laplacian density function. The smaller the value of the shape parameter bl is, the more peaked the distribution will become. To estimate parameter bl (Verdoolaege and Scheunders, 2011), we can initialize bl = 1, and use

ml ¼ min ml

Rl ¼

Ll  bl X ðzp  ml ÞT ðzp  ml Þ ;

ð11Þ

i¼1

Ll  bl bl X ðzp  ml ÞT R1 l ðzp  ml Þ Ll i¼1

!1=bl ð12Þ

to get the optimal ml and Rl. Ll is the number of samples in this case. With

Table 4 Change detection results of the farmland area of the Yellow River Estuary obtained by different methods. The best value for each metric is shown in bold. Method

Difference images

Miss-alarm rate (%)

PCC

Kappa

K&I

IDI TDI ITDIs

False-alarm rate (%) 3.54 0.00 1.90

0.44 2.93 1.10

0.9602 0.9707 0.9701

0.7130 0.6569 0.7471

FCM

IDI TDI ITDIs

14.25 2.24 2.38

0.53 0.16 0.64

0.8522 0.9760 0.9698

0.3641 0.8148 0.7615

G_GC

IDI TDI ITDIs

2.13 0.37 1.63

0.61 2.04 0.62

0.9726 0.9759 0.9776

0.7800 0.7511 0.8133

MGG_GC

ITDIs

0.78

1.05

0.9817

0.8320

131

M. Gong et al. / ISPRS Journal of Photogrammetry and Remote Sensing 93 (2014) 123–135

Ll Ll Wð1=bl Þ þ bl b2l   Ll ib l X 1h T 1  ðzp  ml ÞT R1 ðz  m Þ lnððz  m Þ R ðz  m ÞÞ p p p l l l l l 2 i¼1 ¼ 0; ð13Þ we can get the new bl. W denotes the digamma function. These Eqs. (11) - (12) and (13) can be solved recursively (Do and Vetterli, 2002). The term V p;q ðlp ; lq Þ should be interpreted as a penalty for a discontinuity between pixel p and q, and it is defined as

V p;q ðlp ; lq Þ / exp 

kzp  zq k22 2r 2

! 

1 ; dðp; qÞ

ð14Þ

where k  k2 is the distance between two feature vectors zp and zq. Parameter r is the variance of the quantities kzp  zq k2 for all pairs of p and q in N , and it acts as a term to normalize the vectors’ distance. dðp; qÞ is the distance between pixel p and q on the grid of image in the given neighborhood system. When V p;q ðlp ; lq Þ is small, p and q are likely assigned to different labels. On the contrary, when it is big, p and q are bound to be assigned to the same label. Change detection can be seen as the process of finding the global minimum of the energy E(l), which can be solved by graph cuts

Fig. 10. Final change-detection maps obtained by different approaches relating to the farmland area of the Yellow River Estuary. (a) K&I based on IDI. (b) K&I based on TDI. (c) K&I based on ITDI. (d) FCM based on IDI. (e) FCM based on TDI. (f) FCM based on ITDIs. (g) G_GC based on IDI. (h) G_GC based on TDI. (i) G_GC based on ITDIs. (j) MGG_GC based on ITDIs. (k) Ground truth.

132

M. Gong et al. / ISPRS Journal of Photogrammetry and Remote Sensing 93 (2014) 123–135

Table 5 Change detection results of the river area of the Yellow River Estuary obtained by different methods. The best value for each metric is shown in bold. Method

Difference images

False-alarm rate (%)

Miss-alarm rate (%)

PCC

Kappa

K&I

IDI TDI ITDIs

0.68 0.38 0.42

1.20 1.00 1.02

0.9812 0.9862 0.9856

0.6805 0.7618 0.7523

FCM

IDI TDI ITDIs

1.65 1.71 1.38

0.34 0.20 0.27

0.9801 0.9809 0.9835

0.7375 0.7552 0.7775

G_GC

IDI TDI ITDIs

3.18 1.37 1.09

0.16 0.22 0.29

0.9666 0.9840 0.9862

0.6365 0.7858 0.8055

MGG_GC

ITDIs

0.98

0.31

0.9872

0.8165

algorithm (Boykov and Funka-Lea, 2006). Let g ¼ hm; ei be a weighted graph, where m are a set of nodes and e are a set of directed edges that connect them. There are also two additional terminal nodes S (source) and T (sink) that represent ‘‘changed’’ and ‘‘unchanged’’ labels. Normally, there are two types of edges in the graph: N-links and T-links. Edges between pixels are called N-links and they represent a neighborhood system in the image. T-links connect pixels with terminal nodes S and T. Both N-link and T-link are assigned some non-negative weight. The cost of a N-link is normally derived from the term V p;q ðlp ; lq Þ. The cost of T-link is usually derived from the data term Dp ðlp Þ. The algorithm introduced in (Boykov and Kolmogorov, 2004) can be used to find the optimized cut.

3. Experimental results To demonstrate the effectiveness of our proposed method, we will display the numerical results on three different data sets. All the experiments aim at verifying the following points: (1) the proposed ITDIs can lead to better performance compared with IDI and TDI; (2) change detection with graph cuts based on MGG has a higher accuracy compared with traditional ones. With regards to the first point, three change detection methods (K&I, FCM and graph cuts) are applied to IDI, TDI and ITDIs in order to verify the effectiveness of ITDIs. The details of K&I based on univariate model and multivariate model can be referred to (Bazi et al., 2005; Delong and Junbin, 2009; Hao and Zhu, 2005). In the formulation of graph cuts’ energy function, Gaussian distribution and multivariate Gaussian distribution (denoted as G_GC) are respectively used. The second point is to compare the MGG-based graph cuts (MGG_GC) with traditional change detection methods under the same situation, where they all use ITDIs. To get TDI, several parameters should be chosen. In the step of getting N GDIs, the bank of Gabor filters are set with five scales (2, 4, 6, 8 10) and twelve orientations (0°, 15°, 30°, 45°, 60°, 75°, 90°, 105°, 120°, 135°, 150°, 165°). Phase offset u = 0, spatial aspect ratio c = 0.5, and bandwidth b = 2 are used. The interpretation of these parameters can be referred to (Jain and Farrokhnia, 1991). In the step of decomposing M GDI composed by these 60 Gabor responses with RPCA, Alternating Direction Methods (AMD) is used to get the results. The sparsity ratio is set as pffiffiffiffiffi 0.15, rank as 1, and regularization parameter as 1= m, as recommended by (Wright et al., 2009), where m is the number of rows of matrix MGDI . The quantitative analysis of change detection results are based on the following values: false-alarm rate, miss-alarm rate, PCC, and Kappa (Congalton, 1991; Moser and Serpico, 2006), of which Kappa is a measure of accuracy or agreement based on the difference between the error matrix and chance agreement.

3.1. Results on the Ottawa data set The first experiment is tested on two multitemporal SAR images of 290  350 pixel-sized each, which are acquired by a RADARSAT SAR sensor and provided by Defence Research and Development Canada Ottawa. These images were registered by the automatic registration algorithm from A.U.G. Signals Ltd. that is available through the distributed computing at . There are two regions (water and land) in these images. The images are shown in Fig. 2. As shown in Table 1, the change detection results of different methods based on IDI, TDI and ITDIs are compared separately. For K&I method, TDI makes a lower value in false-alarm rate. This phenomenon can be seen when Fig. 3(a) and (b) are compared. TDI has more ability to resist speckle noise. However, if we focus on the region marked by red1 rectangle in Fig. 3(b), TDI will induce some details missing. Luckily, ITDIs find a tradeoff between noise resistance and detail preservation. Fig. 3(c) can prove this, although PCC of ITDIs is little worse than that of TDI. For FCM, IDI and TDI resist noise approximately. However, comparing the red rectangle marked in Fig. 3(d) and (e) with the ground truth, it can be found that IDI leads to higher miss-alarm and TDI leads to higher falsealarm. ITDIs also do not give a progress in Kappa and PCC. The change detection accuracy of G_GC is not better than those of FCM since Gaussian distribution sometimes cannot fit with the real data distribution. With our extension, MGG_GC gives an outstanding result compared with other methods. Its PCC is the highest and Kappa reaches to 0.9335. A visual comparison showed in Fig. 3 confirms that MGG_GC based on ITDIs can effectively reflect the real changes of the land cover, and the change-detection map we get is much closer to the ground-truth image. 3.2. Results on the Bern data set The experiment is tested on two multitemporal SAR images of 301  301 pixel-sized each acquired by the European Remote Sensing 2 satellite SAR sensor over an area near the city of Bern, Switzerland, in April and May 1999, respectively. Between the two dates, the river Aare flooded parts of the cities of Thun and Bern and the airport of Bern entirely. The two images are shown in Fig. 4. The changed areas only occupy a small percentage of the whole image compared with the unchanged ones. In Table 2, the performance is illustrated. Firstly, it can be seen that TDI-based methods (K&I, FCM, or G_GC) have a lower falsealarm rate. Moreover, most methods using TDI work as well as ITDIs. This means that TDI gives primary contribution for this data set. Secondly, with the addition of IDI, ITDIs give lower value in 1 For interpretation of color in Fig. 3, the reader is referred to the web version of this article.

M. Gong et al. / ISPRS Journal of Photogrammetry and Remote Sensing 93 (2014) 123–135

miss-alarm rate and have more ability to discover weak changes and maintain fine details. Thirdly, from the last row in Table 2, it can be seen that ITDIs-based MGG_GC gives the best results of all. Its false-alarm rate equals 0.16%, miss-alarm rate equals 0.17%, PCC equals 0.9967, and Kappa equals 0.8705. The change-detection maps achieved by these methods are illustrated in Fig. 5. From Fig. 5(b), (e), and (h), it can be seen that TDI-based methods have very powerful ability to resist noise. Fig. 5(g) shows that traditional univariate graph cuts may cause more distortion. In particular, the combination of ITDIs and

133

MGG_GC showed in Fig. 5(j) gives the best change-detection map with a tradeoff between detail and noise. 3.3. Results on the Yellow River Estuary data set The data set used in the experiment consists of two SAR images acquired by Radarsat-2 at the region of Yellow River Estuary in China in June 2008 and June 2009, respectively. The original size of these two SAR images acquired by Radarsat-2 is 7666  7692. The changed areas are scattered and small, and cannot be

Fig. 11. Final change-detection maps obtained by different approaches relating to the river of the Yellow River Estuary. (a) K&I based on IDI. (b) K&I based on TDI. (c) K&I based on ITDI. (d) FCM based on IDI. (e) FCM based on TDI. (f) FCM based on ITDIs. (g) G_GC based on IDI. (h) G_GC based on TDI. (i) G_GC based on ITDIs. (j) MGG_GC based on ITDIs. (k) Ground truth.

134

M. Gong et al. / ISPRS Journal of Photogrammetry and Remote Sensing 93 (2014) 123–135

discovered obviously in the whole image. Here three different changed areas are used in the experiments. Especially, the image acquired in 2008 is four-look, and the one acquired in 2009 is single-look. So the influence of speckle noise on the later image is much greater than the former one. Besides, the resolution of this data set is about 8 m, which is much higher than Ottawa and Bern data sets. The three typical areas can represent the changes of different geographic sites of Yellow River Estuary. Fig. 6 shows the multitemporal images of coastline area, where the changes occurring on the surface of the sea along the coast only occupy a small part. Fig. 7(a) and (b) shows a block of farmland which is landlocked, and the changed areas are relatively large and regular. Fig. 8 shows a section of the river, and it is an inland water area. The changes between the two images lie in the borderline of the river, and the shapes of the changed area are relatively irregular.

3.3.1. Results on the coastline of the Yellow River Estuary Table 3 shows the performance of different methods on the coastline of Yellow River Estuary. Firstly, these three methods can be compared under the same conditions. From the second, fifth, and eighth rows, IDI-based graph cuts are prior to K&I and FCM, especially in false-alarm rate. This phenomenon can be observed in Fig. 9(a), (d) and (g). Graph cuts approach has more ability to resist noise, which can contribute to the smoothness constraint term in its energy function. For TDI, K&I and graph cuts give better performance. K&I causes a higher value in miss-alarm and graph cuts cause a higher value in false-alarm. The visualized comparison can be seen in Fig. 9(b) and (h). As for ITDIs, from the fourth, seventh and tenth rows, it can be seen that graph cuts are superior in the PCC and Kappa. Secondly, focus is placed on K&I, FCM and graph cuts separately to see how much progress is made when using ITDIs compared with IDI and TDI. From the second row to the fourth, it can be seen that TDI-based and ITDIs-based K&I both make progress in PCC and Kappa. Combined with TDI, K&I has more ability to resist noise. For FCM, IDI has a higher false-alarm, which is expressed by noises in Fig. 9(d). TDI-based FCM has the lowest miss-alarm. However, Fig. 9(e) shows the changed region detected by this means is much larger than that of the ground-truth data. And ITDIs-based FCM find a trade between these two points. As for graph cuts, although IDI gives a relatively good result, from the eighth row to the tenth it can be observed that graph cuts also make a further progress when using ITDIs. This phenomenon can be observed when we compare the differences of Fig. 9(g–i).

Fig. 11 (continued)

Finally, from the last row of Table 3, MGG_GC gives the best result of all. The change-detection map got by it has smoother edges and is much closer to the ground-truth image, as showed in Fig. 9(j). 3.3.2. Results on the farmland of the Yellow River Estuary The performance of the ten different methods on the farmland area of Yellow River Estuary is listed in Table 4. In general, TDI has a good characteristic to resist speckle noise. However, its performance is not always better than IDI. For example, for the method of K&I and G_GC, the value of Kappa is worse compared with IDI when we use TDI alone. From Fig. 10(b) and (h) it can be noted that some changed regions are missed. The value of Kappa of TDI-based FCM is high attributing to its low false alarm and low miss alarm. However, Fig. 10(e) shows that the quality of the figure is not very good, because some detected regions have merged together. When ITDIs is used, most methods find a tradeoff between detail and noise, and cause better performance. In addition, ITDIs-based MGG_GC gives the highest accuracy, which can be observed obviously in Fig. 10(j). 3.3.3. Results on the river of the Yellow River Estuary Table 5 shows the performance of the different methods on the river of Yellow River Estuary. For this area, TDI makes a progress in PCC and Kappa compared with IDI. Except K&I, ITDIs-based FCM and graph cuts have further growth. MGG_GC gives the best result of all. Fig. 11 gives visualized results. From Fig. 11(a–c), we can find TDI-based and ITDIs-based K&I detect more changed pixels. When focusing on the region marked by red rectangle in Fig. 11(d–f), it can be found that although TDI-based FCM has more ability of resisting noise, it cannot grasp the shape of the changed region exactly. ITDIs-based FCM performs better. Fig. 11(g) and (h) are produced by IDI-based and TDI-based graph cuts, respectively. In the former figure, some distortion can be observed. As for the latter, some details miss. Above all, Fig. 11(i) and (j) are much closer to the ground-truth image. 4. Conclusions In this paper, a novel unsupervised change-detection approach specifically oriented to the analysis of multitemporal SAR images has been presented. This approach is based on independent construction of IDI and TDI. The automatic analysis of these two difference images is carried out by the extension of graph cuts. Three novel methodological contributions characterize this work compared to traditional change-detection techniques: (1) it generates two difference images exploiting both intensity and texture information; (2) RPCA is used to find the low rank subspace and separate the irrelative and noisy information in Gabor features; (3) graph cuts are used for change detection and have been extended. With regards to the first contribution, although there are several methods to get difference image, most researches are based on a single difference image. The model proposed by this paper is attractive, because it gives a new idea to make use of the merits of IDI and TDI. The second contribution is the first use RPCA to decompose Gabor responses. Low rank property of these Gabor responses is used to separate some irrelevant and noisy elements. Inspired by this step, more outstanding features can be obtained in the following research. As mentioned in the third contribution, in the formulation of graph-cuts’ energy function MGG is used for more accurately fitting. The experimental results show the effectiveness of the proposed approach. As expected, the combination of IDI and TDI can effectively remove noise and give better performance. In addition,

M. Gong et al. / ISPRS Journal of Photogrammetry and Remote Sensing 93 (2014) 123–135

it has been proved that MGG-based graph cuts can lead to higher change detection accuracy. However, because it is difficult to gather reliable ground truth data for diverse land cover changes, the data sets used in the experiments are all connected with the changes caused by the extension of water surfaces. In our future work, we will try to demonstrate the proposed method on diverse land cover changes. Acknowledgements The authors wish to thank the editors and anonymous reviewers for their valuable comments and helpful suggestions which greatly improved the paper’s quality. This work was supported by the National Natural Science Foundation of China (Grant No. 61273317), the National Top Youth Talents Program of China, the Specialized Research Fund for the Doctoral Program of Higher Education (Grant No. 20130203110011) and the Fundamental Research Fund for the Central Universities (Grant No. K5051202053). References Bazi, Y., Bruzzone, L., Melgani, F., 2005. An unsupervised approach based on the generalized Gaussian model to automatic change detection in multitemporal SAR images. IEEE Trans. Geosci. Remote Sens. 43, 874–887. Bazi, Y., Melgani, F., Al-Sharari, H.D., 2010. Unsupervised change detection in multispectral remotely sensed imagery with level set methods. IEEE Trans. Geosci. Remote Sens. 48, 3178–3187. Blaschke, T., 2010. Object based image analysis for remote sensing. ISPRS J. Photogramm. Remote Sens. 65, 2–16. Bovolo, F., Bruzzone, L., 2005. A detail-preserving scale-driven approach to change detection in multitemporal SAR images. IEEE Trans. Geosci. Remote Sens. 43, 2963–2972. Boykov, Y., Funka-Lea, G., 2006. Graph cuts and efficient ND image segmentation. Int. J. Comput. Vision 70, 109–131. Boykov, Y., Kolmogorov, V., 2004. An experimental comparison of min-cut/maxflow algorithms for energy minimization in vision. IEEE Trans. Pattern Anal. Mach. Intell. 26, 1124–1137. Boykov, Y., Veksler, O., Zabih, R., 2001. Fast approximate energy minimization via graph cuts. IEEE Trans. Pattern Anal. Mach. Intell. 23, 1222–1239. Bruzzone, L., Prieto, D.F., 2000. Automatic analysis of the difference image for unsupervised change detection. IEEE Trans. Geosci. Remote Sens. 38, 1171– 1182. Bruzzone, L., Prieto, D.F., 2002. An adaptive semiparametric and context-based approach to unsupervised change detection in multitemporal remote-sensing images. IEEE Trans. Image Process. 11, 452–466. Candès, E.J., Recht, B., 2009. Exact matrix completion via convex optimization. Found. Comput. Math. 9, 717–772. Candès, E.J., Li, X., Ma, Y., Wright, J., 2011. Robust principal component analysis? J. ACM (JACM) 58, 11. Celik, T., 2010. Change detection in satellite images using a genetic algorithm approach. IEEE Geosci. Remote Sens. Lett. 7, 386–390. Congalton, R.G., 1991. A review of assessing the accuracy of classifications of remotely sensed data. Remote Sens. Environ. 37, 35–46. Delong, Z., Junbin, Z., 2009. Minimum error thresholding based on two dimensional histogram. In: Computer Science and Information Engineering, 2009 WRI World Congress on. IEEE, pp. 169–175. Deng, H., Clausi, D.A., 2004. Gaussian MRF rotation-invariant features for image classification. IEEE Trans. Pattern Anal. Mach. Intell. 26, 951–955. Deng, H., Clausi, D.A., 2005. Unsupervised segmentation of synthetic aperture radar sea ice imagery using a novel Markov random field model. IEEE Trans. Geosci. Remote Sens. 43, 528–538. Do, M.N., Vetterli, M., 2002. Wavelet-based texture retrieval using generalized Gaussian density and Kullback–Leibler distance. IEEE Trans. Image Process. 11, 146–158. Florindo, J.B., Bruno, O.M., 2012. Fractal descriptors based on Fourier spectrum applied to texture analysis. Phys. A: Stat. Mech. Appl. 391 (20), 4909–4922.

135

Gong, M., Zhou, Z., Ma, J., 2012. Change detection in synthetic aperture radar images based on image fusion and fuzzy clustering. IEEE Trans. Image Process. 21, 2141–2151. Hao, Y., Zhu, F., 2005. Fast algorithm for two-dimensional Otsu adaptive threshold algorithm. J. Image Graph. 4, 014. He, C., Zhao, Y., Wei, A., 2010. Land-use/land-cover change detection by using the extended change-vector analysis. In: Information Science and Engineering (ICISE), 2010 2nd International Conference on. IEEE, pp. 3809–3812. Hussain, M., Chen, D., Cheng, A., Wei, H., Stanley, D., 2013. Change detection from remotely sensed images: From pixel-based to object-based approaches. ISPRS J. Photogramm. Remote Sens. 80, 91–106. Jain, A.K., Farrokhnia, F., 1991. Unsupervised texture segmentation using Gabor filters. Pattern Recogn. 24, 1167–1186. Jain, A.K., Karu, K., 1996. Learning texture discrimination masks. IEEE Trans. Pattern Anal. Mach. Intell. 18, 195–205. Kandaswamy, U., Adjeroh, D.A., Lee, M.-C., 2005a. Efficient texture analysis of SAR imagery. IEEE Trans. Geosci. Remote Sens. 43, 2075–2083. Kandaswamy, U., Adjeroh, D.A., Lee, M., 2005b. Efficient texture analysis of SAR imagery. IEEE Trans. Geosci. Remote Sens. 43, 2075–2083. Keshavan, R.H., Montanari, A., Oh, S., 2010. Matrix completion from noisy entries. J. Mach. Learn. Res. 99, 2057–2078. Kim, J.-S., Hong, K.-S., 2009. Color–texture segmentation using unsupervised graph cuts. Pattern Recogn. 42, 735–750. Kolmogorov, V., Zabin, R., 2004. What energy functions can be minimized via graph cuts? IEEE Trans. Pattern Anal. Mach. Intell. 26, 147–159. Kruizinga, P., Petkov, N., 1999. Nonlinear operator for oriented texture. IEEE Trans. Image Process. 8, 1395–1407. Li, Y., He, M., 2006. Texture-based segmentation of high resolution SAR images using contourlet transform and mean shift. In: Information Acquisition, 2006 IEEE International Conference on. IEEE, pp. 201–206. Ma, W.-Y., Manjunath, B., 2000. EdgeFlow: a technique for boundary detection and image segmentation. IEEE Trans. Image Process. 9, 1375–1388. Manian, V., Vasquez, R., 1997. A framework for SAR image classification: comparison of co-occurrence and a Gabor based method. In: Geoscience and Remote Sensing, 1997. IGARSS’97. Remote Sensing-A Scientific Vision for Sustainable Development, 1997 IEEE International. IEEE, pp. 335–337. Manjunath, B.S., Ma, W.-Y., 1996. Texture features for browsing and retrieval of image data. IEEE Trans. Pattern Anal. Mach. Intell. 18, 837–842. Min, W., Shu-dao, Z., Heng, B., Ning, M., Song, Y., 2010. SAR water image segmentation based on GLCM and wavelet textures. In: Wireless Communications Networking and Mobile Computing (WiCOM), 2010 6th International Conference on. IEEE, pp. 1–4. Moser, G., Serpico, S.B., 2006. Generalized minimum-error thresholding for unsupervised change detection from SAR amplitude imagery. IEEE Trans. Geosci. Remote Sens. 44, 2972–2982. Oliver, C., Quegan, S., 2004. Understanding Synthetic Aperture Radar Images. SciTech Publishing. Pelizzari, S., Bioucas-Dias, J., 2007a. Oil spill segmentation of SAR images via graph cuts. In: Geoscience and Remote Sensing Symposium, 2007. IGARSS 2007. IEEE International. IEEE, pp. 1318–1321. Pelizzari, S., Bioucas-Dias, J.M., 2007b. Bayesian oil spill segmentation of SAR images via graph cuts. Pattern Recognition and Image Analysis. Springer, pp. 637–644. Rignot, E.J., van Zyl, J.J., 1993. Change detection techniques for ERS-1 SAR data. IEEE Trans. Geosci. Remote Sens. 31, 896–906. Schubert, A., Faes, A., Kääb, A., Meier, E., 2013. Glacier surface velocity estimation using repeat TerraSAR-X images: wavelet-vs. correlation-based image matching. ISPRS J. Photogramm. Remote Sens. 82, 49–62. Torres-Torriti, M., Jouan, A., 2001. Gabor vs. GMRF features for SAR imagery classification, image processing, 2001. In: Proceedings. 2001 International Conference on. IEEE, pp. 1043–1046. Ulaby, F.T., Kouyate, F., Brisco, B., Williams, T.L., 1986. Textural information in SAR images. IEEE Trans. Geosci. Remote Sens., 235–245. Verdoolaege, G., Scheunders, P., 2011. Geodesics on the manifold of multivariate generalized Gaussian distributions with an application to multicomponent texture discrimination. Int. J. Comput. Vision 95, 265–286. Villasensor, J., Fatland, D.R., Hinzman, L.D., 1993. Change detection on Alaska’s North Slope using repeat-pass ERS-1 SAR images. IEEE Trans. Geosci. Remote Sens. 31, 227–236. Wright, J., Ganesh, A., Rao, S., Ma, Y., 2009. Robust principal component analysis?. In: Proceedings of the Conference on Neural Information Processing Systems (NIPS’09).