A Novel Change Detection Method Based on Statistical Distribution Characteristics Using Multi-Temporal PolSAR Data

Unsupervised change detection approaches, which are relatively straightforward and easy to implement and interpret, and which require no human intervention, are widely used in change detection. Polarimetric synthetic aperture radar (PolSAR), which has an all-weather response capability with increased polarimetric information, is a key tool for change detection. However, for PolSAR data, inadequate evaluation of the difference image (DI) map makes the threshold-based algorithms incompatible with the true distribution model, which causes the change detection results to be ineffective and inaccurate. In this paper, to solve these problems, we focus on the generation of the DI map and the selection of the optimal threshold. An omnibus test statistic is used to generate the DI map from multi-temporal PolSAR images, and an improved Kittler and Illingworth algorithm based on either Weibull or gamma distribution is used to obtain the optimal threshold for generating the change detection map. Multi-temporal PolSAR data obtained by the Radarsat-2 sensor over Wuhan in China are used to verify the efficiency of the proposed method. The experimental results using our approach obtained the best performance in East Lake and Yanxi Lake regions with false alarm rates of 1.59% and 1.80%, total errors of 2.73% and 4.33%, overall accuracy of 97.27% and 95.67%, and Kappa coefficients of 0.6486 and 0.6275, respectively. Our results demonstrated that the proposed method is more suitable than the other compared methods for multi-temporal PolSAR data, and it can obtain both effective and accurate results.


Introduction
Change detection is an important remote sensing technology that is used to identify the changes of the Earth's surface through multi-temporal images of the same geographical area observed at different times [1]. On the one hand, due to the impact of environmental factors and social development, changes occur all the time. On the other hand, as a result of the development of satellite systems, a huge number of remote sensing images can now be acquired to detect these changes. Owing to the explosive increase behaviors of a change detection DI map. Therefore, in order to obtain a precise threshold, Weibull or gamma distribution is used to modify the K&I algorithm.
Above all, we propose a novel change detection method based on the omnibus likelihood ratio test statistic and improved K&I algorithm using the Weibull or gamma distribution in this paper. The appropriate DI map is generated by an omnibus likelihood ratio test statistic from the multi-temporal covariance (or coherency) matrix. Moreover, to estimate the distribution properties of the DI map, statistical histogram curve-fitting is utilized to determine the distribution model. Furthermore, Weibull or gamma distribution is used to improve the K&I algorithm, as well as obtain the optimal threshold. Finally, an accurate change detection map is generated using the DI map and the optimal threshold. The rest of this paper is organized as follows: the principle and detailed procedures of the proposed method are described in Section 2. Section 3 introduces the case study. Section 4 provides the experimental results and analyses. Potentials and limitations of the proposed method are briefly discussed in Section 5. Finally, our conclusions are drawn in Section 6.

The Model of PolSAR Data
Fully polarimetric information is used in the proposed method, which can be expressed by Equation (1) when the target reciprocity condition is satisfied.
where h and v represent horizontal and vertical polarization, respectively. The scattering vector S hv represents the vertical transmitting and horizontal receiving polarizations. T denotes the transpose operation.
After multi-look processing, as an example of a covariance matrix, the backscattered information matrix of the ground target follows a Wishart distribution and can be expressed as follows: After image preprocessing, the co-registered and equal-sized PolSAR images from the same geographical region can be represented by X t1 (i, j), X t2 (i, j), . . . , X tn (i, j) at time t1, t2, . . . , tn, respectively. Because of the different data obtained over a long time interval, we can assume that the temporal PolSAR measurements (X t1 , X t2 . . . X tn ) are independent and can be expressed as follows: . . .
where p represents the matrix dimension of X t1 , X t2 , . . . , X tn , and m 1 , m 2 , . . . , m n represents the number of looks of X t1 , X t2 , . . . , X tn , respectively. The probability density function (PDF) of the X ti is expressed by Equation (4). Σ Xti denotes the dispersion matrix of the i-th PolSAR image and is calculated by Equation (5).
Comparing X t1 , X t2 , . . . , X tn , when these dispersion matrices satisfy Σ Xt1 = Σ Xt2 = . . . = Σ Xtn , this means that the ground targets in these regions at different dates show the characteristic of being unchanged, and they can be defined by the null hypothesis (H 0 ). In contrast, when these dispersion matrices satisfy Σ Xti Σ Xtj and i, j denote the arbitrary time, this means that at least one change happened in these regions, and it can be defined by the alternative hypothesis (H 1 ).
The joint densities of omnibus test statistics based on MLE can be expressed by f (Σ Xt1 , Σ Xt2 , . . . , Σ Xti , ε), where the parameter ε is the set of the probability function. Furthermore, the omnibus test statistic (Q) using the likelihood ratio is expressed by Equation (6).
where L( ) and f ( ) are the likelihood function and frequency function, respectively, and Ω = H 0 ∪H 1 . When combined with Equations (3) and (5),Σ can be expressed as follows: If we assume that Σ Xt1 = Σ Xt2 = . . . =Σ Xtn =Σ, and Equation (4) is input into Equation (6), then Q can be expressed as follows: When combined with Equation (7) and assuming m 1 = m 2 = . . . = m n = m, Equation (8) can be simplified as follows: After a logarithm operation, the DI map (d) with the omnibus test statistic can be denoted by Equation (10). After generating the DI map from the multi-temporal PolSAR images, an automatic thresholding method is introduced. To obtain good change detection results, an appropriate thresholding algorithm is needed to separate the DI into the changed and unchanged classes. We recommend the K&I algorithm, which obtains the optimal threshold value by minimizing the classification error [38]. As the statistical distribution function is considered in this algorithm, it is best to establish which distribution model is most compatible with the DI map. Typically, Gaussian and generalized Gaussian distribution models are used in remote sensing imagery. Due to the impact of speckle noise, the negative exponential model is widely used in SAR imagery [47]. Furthermore, the Weibull distribution can be compatible with PolSAR data over a wide resolution range, and the PDF involves a scale parameter θ and a shape parameter γ, as follows: The gamma distribution is also a two-parameter distribution model, which can degenerate to an exponential or chi-squared distribution in some circumstances. Moreover, this distribution model is fit for PolSAR imagery with a medium resolution. The PDF of the gamma distribution involves a rate parameter θ, a shape parameter γ, and a gamma function Г(γ), as follows: A visual way to determine if the distribution model is suitable or not is distribution function fitting. After we make the distribution function fit in some homogeneous areas via MLE operation, we can observe the goodness of fit by comparing it with the statistical histogram. The homogeneous areas can be lakes or vegetation areas. Moreover, as one of important types of surface coverage, urban area was also selected to analyze its distribution property.

(b) The K&I Algorithm Based on Weibull or Gamma Distribution
In order to obtain a precise threshold, both Weibull and gamma distributions are separately used to modify the K&I algorithm. As an extension of Bayesian theory, the traditional K&I thresholding method can be expressed as shown in Equation (13), where L, T, and h(d l ) denote the numbers of possible gray levels, the threshold, and the histogram of the DI map, respectively. J( ) and c( ) denote the criterion function and cost function, respectively. Under the condition of d l and the threshold T, P(ω i |d l ,T) (i = u,c) denotes the posterior probability of unchanged (u) and changed (c) classes, respectively. The traditional K&I method assumes that the class-conditional distribution follows a Gaussian, generalized Gaussian, Rayleigh, or exponential distribution, but this assumption cannot accurately reflect the distribution of the DI map. After we analyzed the statistical property of the DI, we found that it was more compatible with the two-parameter distribution models, i.e., Weibull or gamma distributions, which is later confirmed in the experimental part. In order to improve the threshold selection, an improved model based on Weibull or gamma distribution is used to describe the statistical behaviors of the changed and unchanged classes in the DI. We, therefore, propose a new and improved version of the K&I algorithm. If we assume that the distribution model is a Weibull distribution, by combining Equations (11) and (13), the new criterion function can be shown as follows: The symbols θ u , γ u , θ c and γ c are the scale and shape parameters of P u (d l | u, T) and P c (d l | c, T), respectively, which are the Weibull distribution functions estimated with the pixels in the DI map by MLE. After expanding the above equation, the criterion is changed to When the gamma distribution is considered, we can generate another criterion function by combining Equations (12) and (13). Following the same operation as above, the K&I algorithm based on a gamma distribution is as follows: In the application of this algorithm, we firstly calculate the statistical histogram of the DI. We then compute criterion J(T) at each histogram level. Finally, we choose the histogram level where criterion J(T) is the minimum as the optimal threshold, i.e., T * = argmin{J(T): T = 0, 1, . . . , L − 1}. In change detection, the area is marked as a changed one if its pixel value is higher than the optimal threshold in the DI.

The Workflow of the Proposed Method
The entire steps in the proposed approach can be demonstrated as follows: • Step (1): The preprocessing that involves co-registration and removal of speckle noise for the test data is firstly conducted. In this study, the precision of the co-registration was less than one pixel, and the refined Lee filter [47] was used to suppress the speckle noise.

•
Step (2): The omnibus test statistic method is used to generate the DI map from PolSAR images at different dates.

•
Step (3): Statistical histogram curve-fitting is utilized to determine the function distribution model in homogeneous and unchanged regions. Weibull and gamma distributions are separately applied to improve the K&I algorithm.

•
Step (4): The pixels in position (i, j) at different dates should be determined as unchanged or changed. If d ij < T, the index of this pixel is equal to zero, which denotes that the pixels in this position are similar or unchanged. Otherwise, the corresponding pixels are deemed as different or changed, and the index is equal to one. After completing the indexes of all the pixels in the multi-temporal images, the change detection map is produced. Figure 1 shows the detailed algorithm flowchart.

Evaluation Criterion
For the quantitative evaluation, the false alarm (FA) rate, the total errors (TE), the overall accuracy (OA), and the Kappa coefficient are used to verify the performance of the proposed method. When the ground truth is available, these indicators can be expressed as follows: where true positives (TP) denote the number of changed points which are detected correctly, true negative (TN) denotes the number of unchanged points which are detected correctly, false positives (FP) denote the number of unchanged points which are detected incorrectly, and false negatives (FN) denote the number of changed points which are detected incorrectly. Moreover, the numbers of changed points and unchanged points are denoted by N c and N u , respectively.

Study Area and Background
As a central city in China, the city of Wuhan is located in the eastern part of Hubei province and lies in the middle reaches of the Yangtze River ( Figure 2). It belongs to the north subtropical monsoon climate zone, and it has abundant rainfall. There are many rivers and lakes in Wuhan, which account for one-quarter of the total city area. Due to urbanization and population growth, Wuhan is one of the most rapidly developing cities in China. As a result, how to detect changes accurately is important for both regional economic development planning and dynamic monitoring. Because of the construction of the East Lake Greenway and the 50-year return period rainfall that occurred between 2015 and 2016, changes in urban facilities and inundation detection are the emphasis of this research. Compared with optical sensors, SAR sensors are a better solution to detect these changes in the scenario of constant rainfall. Two PolSAR images of Wuhan from this period were, thus, used to monitor the changes in urban facilities and those caused by inundation.

Datasets and Image Preprocessing
Because this flooding event happened suddenly and receded rapidly, the timing of SAR data acquisition is extremely critical. Thanks to the C-band Radarsat-2 sensor (with a repeat cycle of 24 days), we acquired near-real-time PolSAR data on 6 July 2016. Furthermore, the archived pre-event PolSAR data were used to detect changes caused by urbanization and inundation in Wuhan. Table 1 shows the parameters of the experimental datasets. In order to ensure the precision of the change detection results, the preprocessing of the experimental datasets included radiometric calibration, image co-registration, and speckle filtering, all of which were completed using two free open-source software packages (NEST and PolSARpro). After the preprocessing, the pixel values at the different dates were directly related to the radar backscattering coefficients and had the same original resolution. Moreover, for the co-registered images, speckle noise was removed by refined Lee filtering based on a 7 × 7 window.

Results and Analyses
The Pauli-RGB images (|S hh -S vv | for red (R), |S hv | for green (G), and |S hh + S vv | for blue (B)) of Wuhan in 2015 and 2016 are shown in Figure 3. In order to verify the efficiency of the proposed method, the regions labeled by the two red boxes (region 1 is East Lake and region 2 is Yanxi Lake) in Figure 3a were selected to provide a detailed assessment. These regions cover urban areas, water bodies, forest, and other urban facilities. To establish which distribution model is most compatible with the data, three typical surface features were chosen to analyze the statistical distribution of the DI map: water, a vegetated area, and an urban area. The exact regions are shown in Figure 3b. The statistical histograms of the chosen regions in the DI map are shown in Figure 4. In addition, distribution models based on Rayleigh, exponential, Gaussian, Weibull, and gamma distributions are denoted by the different color curves in Figure 4, and these curves were fitted by the MLE method with pixels in these regions. In Figure 4, it is apparent that the Weibull distribution (red curve) and gamma distribution (green curve) are more suitable than the other function distributions in the water region (Figure 4a   The first research region was East Lake, for which the image size was 800 × 500 pixels. The land cover of this region is mainly greenway, urban areas, water bodies, and grassland. Due to the greenway re-construction that took place in this period, we focus on a detailed assessment of this area. The Pauli-RGB images at the different dates and the reference data are shown in Figure 5. Furthermore, the reference map was produced by field surveys conducted and visual interpretation from Google Earth. The two DI maps generated by the log-ratio method using single-channel SAR information and the omnibus test statistic method using fully polarimetric information are shown in Figure 6. Both of these methods have high values in changed regions. However, the DI map based on the log-ratio method also has a high response in the unchanged areas, which increases the difficulty of distinguishing changed and unchanged regions. Unlike the log-ratio method using single-channel SAR data, the omnibus test statistic method employs more polarimetric information to obtain the DI map, which results in a low response in the unchanged areas. Therefore, the contrast between unchanged and changed regions when using the proposed method is more obvious than for the log-ratio method. It is obvious from Figure 6a, b that the proposed method can obtain a DI map that is more stable and precise. The results of the comparison experiment based on the different methods are shown in Figure 7. Because the contrast of the DI map obtained using the single-channel SAR data is not obvious, the result obtained using the log-ratio and K&I method has many wrong detections in the unchanged regions (e.g., water body). Comparing Figure 7a, e, the result based on the omnibus test statistic and K&I method shows a better performance than the result of the log-ratio and K&I method. This proves that using more polarization information from different channels of PolSAR can obtain a DI map that is more accurate and has greater contrast. In order to show the effects of the different threshold selection methods, the methods based on significance levels of 5% and 1%, the constant false alarm rate (CFAR) algorithm, and the traditional K&I algorithm were used in the contrast experiment. The results of the change detection (Figure 7b, c, d, e) demonstrate that using the traditional K&I method results in fewer false alarms than the other methods. Although the traditional K&I method shows a better performance, the Gaussian assumption only fits a symmetric distribution, and the DI map obtained using the fully polarimetric information and omnibus test statistic method does not fit the symmetric distribution assumption. The change detection maps based on the different distribution fitting methods for K&I are shown in Figure 7e-h. Comparing Figure 7e, f, due to one more parameter being used, the K&I method based on a generalized Gaussian distribution obtains fewer false detections in the East Lake area, and it obtains more detail information in the greenway area than the K&I method based on a Gaussian distribution. This proves that adding one more parameter in the Gaussian model can improve the accuracy and efficiency of the change detection results. However, the generalized Gaussian assumption, as an extension of the Gaussian distribution, is still designed for a symmetric distribution. Due to this drawback, both the Gaussian and the generalized Gaussian distributions are not a good fit for the DIs from PolSAR images. This results in the K&I method being unable to accurately extract the optimal threshold. Comparing Figure 7f a, g, the result for the gamma distribution shown in Figure 7h shows the best performance, with fewer false alarms and wrong detections than the other function distributions. Moreover, the K&I method based on a Weibull distribution also performs better than the Gaussian and generalized Gaussian distributions. This proves that, to generate a change detection result with more accuracy and efficiency, choosing a suitable function distribution is more important than adding more parameters. The quantitative evaluation of these different methods is given in Table 2, where the proposed method shows the best performance in all four indicators (FA (%), TE (%), OA (%), Kappa). This confirms that our method is superior to the other methods.  Figure 8 shows the effect of the gray level on the K&I algorithm based on different function distributions. Because the generalized Gaussian distribution is similar to the Gaussian distribution and uses only one more parameter, the generalized Gaussian distribution, the Weibull distribution, and the gamma distribution were used to verify the performance of the proposed method. The different function distributions are shown in Figure 8 by the green curve (generalized Gaussian distribution), blue curve (Weibull distribution), and red curve (gamma distribution). In these figures, it can be seen that the generalized Gaussian distribution obtains the best performance when the gray level is equal to 2500, but the precision is lower than our method based on Weibull or gamma distribution. Furthermore, Weibull distribution is easily affected by the different gray levels, and it obtains the best performance when the gray level is equal to 2500. All the results indicate that the proposed method with gamma distribution in different gray levels performs better than the K&I method with other function distributions. The proposed method is also more stable.

Analysis of Inundation and Urban Facilities Changes for Yanxi Lake from 2015 to 2016
The second research region was Yanxi Lake, for which the image size was 600 × 600 pixels. The land cover of this region is mainly urban areas, water bodies, and grassland. Due to the impact of continuous heavy rainfall and construction, the main changes in this area were in the water areas and urban areas. The Pauli-RGB images for the different dates and the ground reference data are shown in Figure 9. The two DI maps generated by the log-ratio method using single-channel SAR information and the omnibus test statistic method using fully polarimetric information are shown in Figure 10. Both of these methods have high values in changed regions. However, the DI map based on the log-ratio method also shows a high response in the unchanged areas, which increases the difficulty of distinguishing changed and unchanged regions. Unlike the log-ratio method using single-channel SAR data, the omnibus test statistic method uses more polarimetric information to obtain the DI map, which results in a low response in the unchanged areas. Therefore, the contrast between the unchanged and changed regions is more obvious than when using the log-ratio method. Moreover, the visual performance shown in Figure 10a,b proves that the proposed method can obtain a DI map that is more stable and precise. The results of the comparison experiment based on the different methods are shown in Figure 11. Comparing Figure 11a, e, the results based on the omnibus test statistic and K&I method show a better performance than the results of the log-ratio and K&I method. This again proves that using more polarization information from different PolSAR channels can obtain a DI map that is more accurate and has greater contrast. The results of the change detection (Figure 11b-h) show that using the K&I method with different distribution results in fewer false alarms than for the other methods. Moreover, the results based on the proposed method obtains fewer false detections in the Yanxi Lake area, and it obtains more detailed information in the urban area. This again proves that, to produce a change detection result with more accuracy and efficiency, choosing a suitable function distribution is more important than adding more parameters. The quantitative evaluation for the different methods is provided in Table 3, where the proposed method based on gamma distribution shows the best performance in all four indicators (FA (%), TE (%), OA (%), Kappa). Moreover, Weibull distribution also shows a better performance than the Gaussian and generalized Gaussian distributions. This confirms that the proposed method is effective and shows a significant improvement in detecting the changes caused by urbanization and inundation over the other methods.  Figure 12 shows the effect of the gray level on the K&I algorithm based on the different function distributions. In the figures, it can be seen that the generalized Gaussian distribution obtains the best performance when the gray level is equal to 3000, but the precision is lower than the other function distributions. The Weibull distribution is stable when the gray levels are above 1500. These results indicate that the proposed method with gamma distribution performs better than the K&I method with other function distributions for the case of different gray levels. Again, the proposed method is also more stable than the other methods.

Discussion
In this section, the potentials and limitations of the proposed method are discussed. The proposed method using multi-temporal PolSAR data has better performance than other methods. In order to verify the usability of our approach, we made further experiments based on multi-temporal single and dual polarization data in Yanxi lake. Single and dual polarization information extracted from the complex target vector in Equation (1) are applied in the proposed method. The change detection results in Yanxi Lake are shown in Figure 13. The results illustrate that the proposed method can also be applied in single or dual polarization data. Moreover, the results of dual polarization data based on Weibull (Figure 13b3, b4) or gamma distribution (Figure 13c3,c4) have good performance. However, compared with the results of dual polarization, the results of single-channel SAR data have more omission detection in Figure 13b1(2), c1 (2). This proves that using more polarization channel information can obtain more accurate change detection results. Comparing Figure 13b1(2), c1(2), the results with single polarization data based on generalized Gaussian distribution have a better performance in Figure 13a1(2). This proves that the symmetric distribution assumption is the most compatible with DI map from the single polarization. Furthermore, quantitative evaluation of results from different polarizations is provided in Figure 14 by the green curve (generalized Gaussian distribution), blue curve (Weibull distribution), and red curve (gamma distribution). The proposed method based on gamma distribution using quad polarization data shows the best performance in TE, OA, and Kappa than others. Moreover, the accuracy (red curve) improves as the polarization increases in Figure 14b-d. This confirms that using more channel information produces more stable and precise results. In addition, the proposed method based on Weibull distribution (blue curve) with dual polarization data has better performance (TE and OA) than single and quad polarization. Although the Kappa is lower than that from quad polarization, the difference is small. These demonstrate that the proposed method based on Weibull distribution is suitable in dual polarization. It proves that our approach has great potential to be applied in Sentinel-1 data with dual polarization datasets.
Although the proposed method has a good performance in change detection and can be applied in both single and dual polarization data, it still has some limitations. Firstly, the proposed method can detect the water changes in two dimensions, but it cannot detect the water-level changes. To solve this problem, external data are applied to improve the change detection accuracy. Secondly, the proposed method can obtain the unchanged and changed information, but it cannot obtain the types of land-cover changes. To solve this problem, classification techniques can be combined with our method to get the types of land-cover changes. Thirdly, multi-temporal time-series observations are ideal for detecting changes accurately. Because of limited data, our research only uses the bi-temporal PolSAR data to detect the changes in urban facility and inundation. The application of the proposed method based on time-series observations from different sensors (e.g., Sentinel-1) will be explored in further work. Figure 13. Change detection results using the HH polarization based on omnibus test statistic and K&I based on (a1) a generalized Gaussian distribution, (b1) a Weibull distribution, and (c1) a gamma distribution; the VV polarization based on omnibus test statistic and K&I based on (a2) a generalized Gaussian distribution, (b2) a Weibull distribution, and (c2) a gamma distribution; the (HH, HV) polarization based on omnibus test statistic and K&I based on (a3) a generalized Gaussian distribution, (b3) a Weibull distribution, and (c3) a gamma distribution; the (VV, VH) polarization based on omnibus test statistic and K&I based on (a4) a generalized Gaussian distribution, (b4) a Weibull distribution, and (c4) a gamma distribution.

Conclusions
In this paper, we showed that the combination of an omnibus test statistic using a covariance matrix and the K&I method based on Weibull or gamma distribution can improve the performance of change detection for multi-temporal PolSAR images. The proposed method for generating DI maps is appropriate for PolSAR images, and it can obtain an accurate intermediate result. Due to the solid mathematical fundamentals of Bayesian theory, the traditional K&I method shows better performance than the other methods. This illustrates that the threshold selection method based on K&I is stable and effective in change detection. After analyzing three typical land covers, i.e., water, vegetation, and urban areas, in the difference image, we find that the variable in the difference image more likely follows the Weibull or gamma distribution. Thus, we proposed the K&I change detection algorithm based on Weibull or gamma distribution. Moreover, the proposed method can be applied to single/dual/quad polarization data, while the different polarization data with different distributions can generate stable and precise change detection results. Based on the proposed method, the changes caused by urban development and inundation between 2015 and 2016 in Wuhan, China, were detected accurately. Over the East Lake and Yanxi Lake regions where the reference map was available, our proposed method obtained the best performance with FA of 1.59% and 1.80%, TE of 2.73% and 4.33%, OA of 97.27% and 95.67%, and Kappa of 0.6486 and 0.6275, respectively. These confirmed that the proposed method is superior to the other methods, and it is effective in detecting changes caused by urbanization and inundation.
Author Contributions: J.Z. identified the difficulties addressed in this research, proposed the new approach, and wrote and revised the paper. Y.C. completed the partial equation derivation. J.Y. provided the multi-temporal Radarsat-2 data. Y.N. helped to preprocess the datasets. Z.L. revised the paper and gave some useful advice. P.L. gave some key advice. All authors have read and agreed to the published version of the manuscript.