Polarimetric synthetic aperture radar change detection for specific land cover types

Abstract

This paper presents a supervised polarimetric synthetic aperture radar (PolSAR) change detection method applied to specific land cover types. For each pixel of a PolSAR image, its target scattering vector can be modeled as having a complex multivariate normal distribution. Based on this assumption, the joint distribution of two corresponding vectors in a pair of PolSAR images is derived. Then, a generalized likelihood ratio test statistic for the equality of two likelihood functions of such joint distribution is considered and a maximum likelihood distance measure for specific land cover types is presented. Subsequently, the Kittler and Illingworth minimum error threshold segmentation method is applied to extract the specific changed areas. Experiments on two repeat-pass Radarsat-2 fully polarimetric images of Suzhou, China, demonstrate that the proposed change detection method gives a good performance in determining the specific changed areas in PolSAR images, especially the areas that have changed to water.

Keywords:

• change detection
• generalized likelihood ratio test
• KI threshold segmentation
• microwave remote sensing

1. Introduction

Recently, many methods have been proposed for change detection in single-channel synthetic aperture radar (SAR) image, and several features can be utilized, such as area, shape, average intensity, texture, and underlying land cover. Compared with single-channel SAR data, polarimetric synthetic aperture radar (PolSAR) data contain both phase and amplitude information from radar returns transmitted in two different polarizations, signifying that more scattering information can be utilized when PolSAR data are available for change detection (Lee and Pottier 2009).

SAR change detection involves four major steps: pre-processing filtering and co-registration, difference map extraction, threshold segmentation, and image fusion. Among these steps, difference map extraction is one of the most important sectors. Recently, some excellent studies have been performed in the field of difference map extraction. Many polarimetric features have been applied, such as polarization state conformation, optimal polarization theory, and fully polarimetric clutter models (Conradsen et al. 2003; Qong 2004; Brisco et al. 2013).

Taking advantage of the discrimination of different scattering mechanisms, terrain and land-use classification has become the most important application of PolSAR (Lee and Pottier 2009). However, in terms of difference map extraction, we still cannot distinguish the land type before and after a change from the difference map generated by the PolSAR change detection method. Furthermore, in certain applications, we do not need all of the changed areas, but only those from specific class w_i to specific class w_j. For example, only the areas that changed into urban areas are useful for research on city expansion, and in this case, other changed areas will be regarded as a false alarm. Therefore, the method of applying the class information of PolSAR to change detection deserves more specific studies. The objective of this paper is to propose a distance measure to detect the areas changing from specific class w_i to specific class w_j and then apply this distance to PolSAR change detection.

2. Proposed distance measure

In this section, we first derive the joint distribution of two random complex vectors. To extract the changed areas from specific class w_i to specific class w_i, a generalized likelihood ratio test (GLRT) is built and a maximum likelihood (ML) distance measure is derived.

2.1. The joint distribution of two random complex Vectors

We say that a p-dimensional random complex vector X follows a complex multivariate normal distribution with mean 0 and covariance matrix Σ, that is,

(1)

if the probability density function (PDF) is (Goodman 1963):

(2)

where denotes the matrix determinant, etr (·) = exp (tr(·)) represents the exponential trace operator, and ^H indicates the complex conjugation ^* and matrix transpose^T. Then the Hermitian random p × p matrix is complex Wishart distributed with n degrees of freedom and is denoted as W ∈ CW (n, p, Σ). The ML estimator of Σ is W/n (Conradsen et al. 2003).

Let us partition

(3)

into q- and (p − q)-component subvectors, respectively. And W and Σ are partitioned as:

(4)

Thus, X⁽¹⁾ ∈ CN (q, Σ₁₁) and X⁽²⁾ ∈ CN (p−q, Σ₂₂).

Note that the independence of X⁽¹⁾ and X⁽²⁾ is not determined, and their joint distribution cannot be calculated directly by multiplying their respective PDF. Here, we make a nonsingular linear transformation to the subvectors

(5)

where the components Y⁽¹⁾ are uncorrelated with the components Y⁽²⁾ = X⁽²⁾. The matrix B should satisfy the equation:

(6)

where indicates the mathematical expectation. Thus, and the vector Y = MX with

(7)

is a nonsingular transform of X, and therefore, Y has a complex normal distribution with mean 0 and covariance matrix

(8)

Let , we have Y⁽¹⁾ ∈ CN (q, Σ_11.2) and (Anderson 2003).

Since Y⁽¹⁾ and Y⁽²⁾ are independently normally distributed, their joint density is:

(9)

The density of X⁽¹⁾ and X⁽²⁾ can then be obtained from this expression by substituting for y⁽¹⁾ and x⁽²⁾ for y⁽²⁾ (the Jacobian of this transformation being 1).

(10)

Under the assumption that the secondary SAR data are independent and identically distributed, the likelihood function of this joint distribution is:

(11)

Note that |Σ| = |Σ_11.2||Σ₂₂| and

(12)

The likelihood function can be recast as a function of only one parameter Σ.

(13)

2.2. Proposed distance measure

The well-known PolSAR scattering vector, based on the complex Pauli spin matrix, is given by

(14)

For two repeat-pass fully PolSAR images, a six-element complex scattering target vector k₆ can be formed by stacking the target vectors with

(15)

where k₁ and k₂ are the target vectors belonging to the two images. Then, the Polarimetric SAR Interferometry Coherency (PolInSAR) matrix is expressed by

(16)

The matrices T₁₁ and T₂₂ are the conventional polarimetric coherency matrices, describing polarimetric properties for each individual image, while T₁₂ is a non-Hermitian complex matrix, which contains polarimetric and interferometric correlation information between k₁ and k₂. It is commonly assumed that the elements of the scattering vector are circular complex multivariate Gaussian with zero mean. Consequently, the unnormalized PolInSAR 2q × 2q matrix W₆ = nT₆ follows a complex Wishart distribution, denoted as W₆ ∈ CW (n, 2q, Σ₆), where n is the sample size used in the average and q is the dimension of k₁ or k₂ (Ferro-Famil, Pottier, and Lee 2001).

Because k₁ ∈ CN (q, Σ₁₁) and k₂ ∈ CN (q, Σ₂₂), the likelihood function of the joint distribution of k₁ and k₂ can be obtained from equation (13).

(17)

Defining the coherency matrix of w_i as T_I, the coherency matrix of w_j as T_J, and the unnormalized PolInSAR matrix of w_i and w_j as W_s = nT_s ∈ (n_s, 2q, Σ_s), we now consider the null hypothesis H₀ : Σ₆ = Σ_s, which states that pixel 1 and the corresponding pixel 2 of these two repeat-pass PolSAR images belong to the areas changing from specific class w_i to specific class w_j, against the alternative hypothesis H₁ : Σ₆ ≠ Σ_s. The GLRT for testing H₀ is (Anderson 2003):

(18)

Substituting the ML estimators of Σ into equation (18), we have

(19)

Maximizing the GLRT is equivalent to minimizing the following proposed distance:

(20)

The null hypothesis H₀ is rejected if the measure d₁ is too large. However, under normal circumstances, it is expected that the extracted changed areas will be represented by a large digital number in the difference map. Therefore, for better display, a deformation

(21)

can be used, where i represents the pixel label for all of the image. Note that d_P has the same value range and opposite monotonicity with d₁, and consequently, the areas changing from specific class w_i to class w_j will have a higher d_P than the other changed areas and those areas that have not changed.

2.3. The reduced distance under uncorrelated assumption

Due to the calculation of matrix T₁₂, the proposed distance d_P in equation (21) can only be obtained from two single-look PolSAR images. However, for data compression and preliminary speckle reduction, the PolSAR data are frequently multi-look processed by spatial averaging at the expense of loss of spatial resolution. In this case, the uncorrelated assumption of the target vectors, k₁ and k₂, should be applied to estimate the value of d_P.

When we assume that k₁ and k₂ are mutually independent, their joint distribution can be obtained directly by multiplying f (k₁) and f (k₂). Therefore, the likelihood function in equation (13) becomes a function of the parameters Σ₁₁ and Σ₂₂:

(22)

And the null hypothesis for testing the equality of such likelihood functions is . According to the derivation method of Section II-B, the distance d₁ reduces to:

(23)

Then the deformation distance in equation (21) becomes

(24)

Note that d₂ can be obtained just from multi-look complex PolSAR data instead of single-look data, which increases the operability of the algorithm. However, compared with the distance d₁, d₂ employs only the polarimetric information, discarding the interferometric information contained within T₁₂. Therefore, when the land cover types of the two images are very close, which means that T₁₂ cannot be ignored, the distance d_R will cause inevitable errors.

3. Change detection method

Based on the proposed distance, a supervised PolSAR detection method for specific land cover type changes is presented. First, two repeat-pass fully PolSAR images should be co-registered in advance. Then, the training areas are selected for class w_i and class w_j from the two images respectively, and their coherency matrices are calculated. By applying the proposed measure d_P, the change from class w_i to class w_j is obtained.

When the difference map is extracted, a binary map is constructed to obtain an image of the change region. The most common solution to this problem is to select the threshold in the histogram of the difference image automatically. The Kittler and Illingworth (KI) minimum error threshold segmentation method is widely used in image segmentation based on Bayesian decision theory (Kittler and Illingworth 1986). Moser and Serpico (2006) propose a generalized KI algorithm, which gives a good performance in the field of change detection. Here, we apply this improved KI method to extract the changed areas from specific class w_i to specific class w_j.

After KI threshold segmentation, we could observe that the segmented binary image had many small unconnected changed areas, and these areas are usually not considered for change detection. In addition, there are some black holes in the changed areas, which increase the difficulty of calculating the changed areas and extracting boundary and reduce the final accuracy of change detection. Therefore, a morphologic operation, connected components extraction, is applied to remove the small unconnected changed areas and small black holes (Gonzalez and Woods 2007). A binary change mask is achieved by morphologic operation and we label the changes in the original images.

4. Experimental results and discussion

To illustrate the change detection capability of the derived GLRT distance measure, two repeat-pass Radarsat-2 fully PolSAR images taken at Shihu Lake area of Suzhou, China, acquired separately on 9 April 2009 and 15 June 2010, were used for our experiment. The spatial resolution is 10 × 8 m, and the incidence angle ranges from 38.37° to 39.85°. There are urban areas, lawns, forest, lake, mountain, and bare land in this area, and thus, it is very suitable for verifying the detectability of specific land cover types. A 7 × 7 neighborhood sliding window size is selected for estimating the polarimetric coherency matrix.

The Pauli decomposition maps are presented in Figure 1. According to field investigations, as well as two Google Earth images of this area, which have a similar imaging time with the experimental data (acquired separately on 15 March 2009 and 19 June 2010), we find three major feature change classes between the two images, labeled from T1 to T3 in Figure 1b. Class T1 indicates the change from bare lands to urban areas, class T2 denotes the change from urban architecture to bare lands, and class T3 indicates the change from lawns or bare lands to water. These three changed areas are delineated as the reference images for accuracy validation.

Figure 1. Experimental data-set. (a) Pauli decomposition map, 9 April 2009. (b) Pauli decomposition map, 15 June 2010. Class T1 indicates the change from bare lands to urban areas, class T2 denotes the change from urban architecture to bare lands, and class T3 indicates the change from lawns to water.

Figure 2 shows the difference maps using the proposed distance d_P and its reduced form d_R, respectively. It can be observed that all of the changed areas marked in Figure 1b have been detected successfully, and the difference maps generated by d_P have less noise and higher contrast than those of d_R. For further analysis, the comparison of these two measures for three change types is shown in Figure 3. When d₂ is equal to d₁, the curve lies along the diagonal line. The average of d₁ and d₂ is shown in Figure 3d. In general, the test statistic for the full PolInSAR matrix in equation (19) is only slightly better than that for the uncorrelated case. We may conclude that the correlation information between these two images is small.

Figure 2. The top images are the difference maps using d_p for three change types: (a) class T1, (b) class T2, and (c) class T3. The bottom images are the difference maps using d_R for three change types: (d) class T1, (e) class T2, and (f) class T3. Class T1 indicates the change from bare lands to urban areas, class T2 denotes the change from urban architecture to bare lands, and class T3 indicates the change from lawns to water.

Figure 3. The comparison of d₁ and d₂ for three change types: (a) class T1, (b) class T2, and (c) class T3. On the abscissa are the d₁ values, and the ordinate contains the d₂ values. (d) The average of d₁ and d₂ for all change types.

The independence of class w_i and class w_j can be verified by:

(25)

If these two classes are orthogonal, V is almost 1. The V values for the three change types are 0.959, 0.908, and 0.997. Therefore, class w_i and class w_j can be regarded as mutually independent, and d_R is an approximation of d_P.

After the extraction of the difference map, we use KI threshold segmentation and connected components extraction to divide the difference map into two classes: changed areas and no-change areas. Figure 4 presents the corresponding color-coded image, with red indicating T1, green denoting T2, and blue indicating T3. The two most important indices for change detection, namely detection rate and false alarm rate of changed areas, are shown in Table 1. Note that all change classes have been detected with a high detection rate and a low false alarm rate and that the detection capability of class T3 gives a better performance than the two other change classes. A possible reason for this is that the sample selection for water areas is easier than the other land cover types, especially for mixed lawn and urban areas in class T2. Accordingly, this method has great potential for application in the monitoring of floods.

Figure 4. Color code for each binary mask.

Table 1. Accuracy of change detection.

Change class	Detection rate (%)	False alarm rate (%)	Kappa coefficient
T1	88.03	6.55	0.905
T2	87.11	6.27	0.901
T3	91.10	3.28	0.936
Total	91.20	8.04	0.911

5. Conclusions

In this paper, we have developed a supervised PolSAR change detection method applied to specific land cover types. To detect the areas changing from class w_i to class w_j between two repeat-pass PolSAR images, a distance measure d_P is derived based on a generalized likelihood ratio test, and the specific changed areas are noted to have a higher d_P than the other changed areas and no-change areas. When we assume that the target vectors belonging to the two images are mutually independent, a reduced distance d_R can be obtained. d_R has fewer constraints than d_P and can be estimated only from multi-look complex PolSAR data. For PolSAR change detection, class w_i and class w_j have a slight correlation, and thus, d_R can be promoted as a substitute of d_P when there are no single-look PolSAR data available.

The effectiveness of the proposed change detection method is illustrated by two Radarsat-2 fully polarimetric SAR images taken at Shihu Lake area of Suzhou, China. The results show that the change areas of all change classes can be detected successfully with a high detection rate and a low false alarm rate. Moreover, the application of this method in the monitoring of floods is worthy of further investigation.

Additional information

Funding

This work was supported in part by the National High-Tech Research and Development Program of China under grant number [2011AA120404] and in part by the National Natural Science Foundation of China under grant numbers [4133176] and [41371352].