A new sequential homogeneous pixel selection algorithm for distributed scatterer InSAR

ABSTRACT

Distributed scatterer interferometric synthetic aperture radar (DS InSAR) technology has been widely used in various fields. Homogeneous pixel selection is a crucial step in the use of DS InSAR, directly affecting the estimation precision and reliability of subsequent parameter calculations. The existing algorithms for selecting homogeneous pixels have inherent limitations, such as requiring many heterogeneous samples and strict requirements surrounding the required number of synthetic aperture radar (SAR) images. To address these problems, a new sequential selection algorithm for homogeneous pixels is proposed, based on the Baumgartner – Weiss–Schindler (BWS) test algorithm and dynamic interval estimation (DIE) theory. According to Monte Carlo simulation experiments, the average standard deviation (STD) of the mean of the rejection of the BWS-DIE algorithm under six sample conditions is 0.014. Compared with three existing algorithms, including the Kolmogorov‒Smirnov (KS), BWS and fast statistically homogeneous pixel selection (FaSHPS) algorithms, the BWS-DIE algorithm improves homogeneous pixel selection precision by 64.3%, 69.4% and 25.3%, respectively. In the real data experiment, 12 scenes of Advanced Land Observing Satellite-1 Phased Array type L-band Synthetic Aperture Radar (ALOS-1 PALSAR) data from February 2007 to March 2011 were used and the BWS-DIE multitemporal InSAR (MT InSAR) method based on the BWS-DIE algorithm was applied to surface subsidence monitoring in the western mining area of Xuzhou, Jiangsu Province, China. The experimental results show that, compared with the Stanford Method for Persistent Scatterers (StaMPS), the BWS-DIE MT InSAR method improves the ability to monitor the maximum subsidence by 12.3%, increases the point density by 5.7 times and decreases the root mean square error (RMSE) by 50%. In addition, new surface deformation patterns are found in the spatial-temporal evolution. The above experimental results show that the proposed BWS-DIE algorithm exhibits remarkable advantages in selection power and selection precision and is not limited by the number of SAR images. The proposed algorithm can further broaden the application scenarios for DS InSAR and provide high-quality and reliable monitoring data for subsequent scientific research.

KEYWORDS:

• Interferometric synthetic aperture radar
• distributed target
• homogeneous pixel selection
• deformation monitoring
• time series analysis

Introduction

Differential interferometric synthetic aperture radar (DInSAR) technology, based on synthetic aperture radar (SAR) phase information, can rapidly obtain surface deformation information across a large range with a high degree of precision. DInSAR has been widely applied in seismic deformation field monitoring (Massonnet et al. 1993, 1994; Lanari et al. 2010; Kim and Han 2023), volcanic eruption deformation monitoring (Massonnet, Briole, and Arnaud 1995; Jónsson et al. 1999; Albino et al. 2020), surface deformation monitoring in mining areas (Chen and Deng 2014; Chen et al. 2015, 2021, 2022; Wang et al. 2022), urban surface deformation monitoring (Chaussard et al. 2013; Zhang et al. 2016; Deng et al. 2017; Wu et al. 2019), landslide monitoring (Ferretti, Prati, and Rocca 2001; Zhang et al. 2018; Dai et al. 2022) and glacier movement monitoring (Goldstein et al. 1993; Li et al. 2018). However, due to factors such as temporal and spatial decorrelation, it is difficult to obtain high-precision measurement results from long-time series monitoring with traditional DInSAR technology. To overcome these limitations, some time series InSAR technologies have been subsequently proposed, such as persistent scatterer InSAR (PS InSAR) (Ferretti, Prati, and Rocca 2000, 2001) and small baseline subset InSAR (SBAS InSAR) (Berardino et al. 2002). Time series InSAR technologies have effectively improved the precision of surface deformation monitoring in long-term series monitoring. However, for non-urban areas, especially mining areas and areas that have been covered by vegetation and bare soil for an extended period, the point target density obtained by the existing method is too sparse. As a result, it is difficult to both provide sufficient deformation details for the deformation field and accurately retrieve the surface temporal-spatial evolution.

To improve monitoring point density, many scholars have proposed the distributed scatterer InSAR (DS InSAR) method (Ferretti et al. 2011; Jiang et al. 2014; Wang, Zhu, and Bamler 2012). Unlike PS point targets, distributed scatterer (DS) targets refer to point targets within a resolution cell with no dominant scatterer and high backscattering intensity (Lee and Pottier 2017). In areas with a large surface subsidence magnitude or high coverage of vegetation and bare soil, DS InSAR technology can substantially improve the density of monitoring point targets and is thus conducive to acquiring high-precision temporal-spatial surface evolution laws.

The recognition of homogeneous pixel sets is the basis for DS target extraction. Currently, there are two main types of recognition method for homogeneous pixel sets (Even and Schulz 2018): non-parametric testing methods and parametric testing methods. Non-parametric testing methods primarily include the two-sample Kolmogorov‒Smirnov (KS) test (Ferretti et al. 2011) and Baumgartner – Weiss–Schindler (BWS) test (Baumgartner, Weiß, and Schindler 1998). The advantages of the KS test algorithm are that the sampling distribution of the statistics is known and that the rejection domain can be calculated using an analytical formula. However, the KS test algorithm is insensitive to tail differences between samples (Stephens 1974), and the test power is relatively low for small samples (Marozzi 2009, 2013). The advantage of the BWS algorithm is that it focuses on the tail weights of the distribution functions and is therefore more sensitive to the differences between the sample’s tails (Jiang, Ding, and Li 2018), which allows it to obtain more power for small samples (Xu et al. 2022). Furthermore, the BWS algorithm can be easily applied and is at least as powerful as the commonly used non-parametric tests. However, the rejection domain of the testing statistics obtained by this method cannot be rendered in the form of a function but needs to be obtained applying a Monte Carlo simulation experiment; in addition, the algorithm’s calculation efficiency is low (Jiang, Ding, and Li 2018).

Another algorithm – fast statistically homogeneous pixel selection (FaSHPS) – is a well-known parametric testing method (Xu et al. 2022; Jiang et al. 2014), and its core idea is to assume that the amplitude of any pixel in the SAR image dataset follows a Rayleigh distribution (Jiang et al. 2016; Oliver and Quegan 2004). As the number of SAR images increases, the average pixel amplitude follows a Gaussian distribution in the time dimension (Jiang et al. 2014). Additionally, the algorithm uses confidence interval estimation theory to replace the traditional hypothesis testing method and determines whether the two pixels satisfy the same function distribution by using logical operations, thus achieving homogeneous pixel recognition. The merit of this algorithm is its high computational efficiency, and the experimental results are more reliable for small samples. However, in the case of many/large samples, if the average amplitude of the reference pixel is used as the true value to construct a confidence interval, the constructed interval may be biased, potentially leading to a reduction in the number of selected homogeneous pixel samples and a reduction in the reliability of homogeneous pixels.

To address the above problems, using the BWS algorithm and dynamic interval estimation theory, this study proposes the BWS dynamic interval estimation (BWS-DIE) method, which is a sequential selection algorithm for homogeneous pixels, i.e. combination of non-parametric testing methods and parametric testing methods. This method utilizes the advantages of the BWS algorithm (i.e. the proposed method can be easily applied and is at least as powerful as the commonly used non-parametric tests and more sensitive to the differences between the sample’s tails). Additionally, the proposed method is more sensitive to differences between the tails of the sample and is more powerful for small samples (Baumgartner, Weiß, and Schindler 1998). The proposed method also has the high computational efficiency of the interval estimation theory (Jiang et al. 2014, 2016). In this study, Monte Carlo simulation experiments were employed to test and evaluate the computational power and degree of precision of the BWS-DIE algorithm. In addition, based on the BWS-DIE algorithm, a new DS InSAR technology process was constructed and applied to surface deformation monitoring in the western mining area of Xuzhou, Jiangsu Province, China. The monitoring results were compared with the monitoring results obtained by employing the KS, BWS and FaSHPS algorithms and StaMPS, and the experimental results verified the feasibility and effectiveness of the proposed method.

Method

Principle of the BWS algorithm

The BWS algorithm is essentially a rank sum testing method. Its core idea is to use the variance of the square of the sample’s cumulative distribution function difference to weight the sample, with the weight depending on the true distribution of the sample. Because the sample’s true distribution is unknown, the statistics for testing are weighted by rank approximation (Baumgartner, Weiß, and Schindler 1998; Jiang, Ding, and Li 2018). Since the rejection domain of the statistics for the test cannot be calculated using function expressions, it is generally necessary to apply Monte Carlo simulation experiments to obtain the corresponding numerical solution.

Suppose that there is a dataset of N registered SAR images and that the corresponding SAR image size is s×s. Let x be the set of SAR image amplitudes arranged in the order of the SAR image acquisition time. Then, for any pixel P_j, the set of amplitudes is expressed as follows:

(1) $x\left(P_j\right)={\left[x_1\left(P_j\right),x_2\left(P_j\right),\dots ,x_i\left(P_j\right)\right]}^T,i=1,2,\dots \dots N;\text{breakj}=1,2,\dots \dots s2$(1)

where T is the transpose and x_i(P_j) is the amplitude of pixel P_j in the i-th SAR image.

To obtain the homogeneous pixel set for all pixels, suppose that each pixel has an estimation window of the same size, the estimation window’s center pixel is the reference pixel and the remaining pixels in the estimation window are the pixels to be estimated. For a fixed-size estimation window, assuming that the reference pixel is P₁ and that one of the pixels to be estimated is P₂, the corresponding time series amplitude samples of P₁ and P₂ are x(P₁) and x(P₂), respectively. According to the principle of the BWS algorithm (Jiang, Ding, and Li 2018), the null hypothesis is that P₁ and P₂ are homogeneous pixels, and the alternative assumption is that P₁ and P₂ are not homogeneous pixels. The BWS hypothesis testing method can be used to construct the test value B (Jiang, Ding, and Li 2018) as follows:

(2) $B_{P_1}=\frac{1}{2N^2}{\displaystyle {\sum }_{i=1}^N\frac{{\left(R_{x_i\left(P_1\right)}-2i\right)}^2}{\frac{i}{N+1}\cdot \left(1-\frac{i}{N+1}\right)}}{B}_{P_2}=\frac{1}{2N^2}{\displaystyle {\sum }_{i=1}^N\frac{{\left(R_{x_i\left(P_2\right)}-2i\right)}^2}{\frac{i}{N+1}\cdot \left(1-\frac{i}{N+1}\right)}}B=\frac{B_{P_1}+B_{P_2}}{2}$(2)

where $R_{x_i\left(P_1\right)}$ and $R_{x_i\left(P_2\right)}$ are the amplitude ranks of the i-th SAR image corresponding to pixels P₁ and P₂, respectively, after combining the time series samples x(P₁) and x(P₂); B is the statistics for the test.

If the statistics for test B fall within the rejection domain, the null hypothesis is rejected. At this time, P₁ and P₂ are not homogeneous pixels; otherwise, they are homogeneous pixels. The above steps are repeated for all pixels to be estimated in addition to the reference pixels in the estimation window; then, a set of homogeneous pixels for a specific reference pixel is obtained. The estimation window is then moved to traverse all of the reference pixels in the SAR image, and the above steps are repeated until the homogeneous pixel set ${\mathrm{\Omega }}_{BWSHP}$ for the SAR image dataset is obtained.

The above analysis reveals that the BWS test algorithm considers the tail weights of the distribution functions and can obtain more power for small samples. The algorithm does not require a subjective classification of the data and is easily implemented and suitable for the initial selection of homogeneous pixels.

Interval estimation theory

Interval estimation theory is a parameter estimation theory (Severini 1991) that was introduced into the DS point selection procedure by (Jiang et al. 2014). The core idea of interval estimation theory is to construct a corresponding confidence interval based on a certain confidence level by using the mean amplitude for a specified pixel. If the mean amplitude of a pixel to be estimated falls within the interval, the pixel is regarded as homogeneous (Jiang et al. 2014, 2016).

According to Equation (1)(1) $x\left(P_j\right)={\left[x_1\left(P_j\right),x_2\left(P_j\right),\dots ,x_i\left(P_j\right)\right]}^T,i=1,2,\dots \dots N;\text{breakj}=1,2,\dots \dots s2$(1) , the set of amplitudes for any reference pixel P_j is x(P_j), and the point estimate of the mean amplitude of pixel P_j is expressed as follows:

$\bar{x}\left(P_j\right)$=$(x_1\left(P_j\right)+x_2\left(P_j\right)+\dots +x_i\left(P_j\right))/N(3)$

In accordance with the central limit theorem, as the number of samples N increases, $\bar{x}\left(P_j\right)$ moves closer to the normal distribution. Therefore, the interval estimation of pixel P_j can be obtained as follows (Jiang et al. 2014):

(4) $P\left\{\bar{x}\left(P_j\right)-z_{1-\alpha /2}\cdot \sqrt{Var\left(P_j\right)/N}<\mu \left(P_j\right)<\bar{x}\left(P_j\right)+z_{1-\alpha /2}\cdot \sqrt{Var\left(P_j\right)/N}\right\}=1-\alpha$(4)

where $P\left\{\hspace{0.17em}\right\}\mathrm{i}\mathrm{s}\mathrm{t}\mathrm{h}\mathrm{e}$ probability; $\bar{x}\left(P_j\right)$and Var(P_j) are the mean and variance of the amplitudes of reference pixel P_j, respectively; α denotes the significance level; z_1-α/2 is the $1-\alpha /2$ quantile of the standard normal distribution; and μ(P_j) is the mean amplitude of the statistical population of reference pixel P_j. According to the research results in the literature (Lee and Pottier 2017), $\sqrt{Var\left(P_j\right)}/\bar{x}\left(P_j\right)$ can be approximated to 0.52, so Equation 9(9) ${\gamma }_{DSC}=\frac{1}{N\left(N-1\right)}{\sum }_{t=1}^N{\sum }_{r=t+1}^{N}ex{p}^{j\left({\mathrm{\varnothing }}_{rt}-({\mathrm{\varnothing }}_{or}-{\mathrm{\varnothing }}_{ot}\right))}$(9) can be rewritten as follows:

(5) $P_j\left\{\bar{x}\left(P_j\right)-z_{1-a/2}\cdot 0.52\cdot \bar{x}\left(P_j\right)/\sqrt{N}<\mu \left(P_j\right)\hspace{0.17em}<\bar{x}\left(P_j\right)\text{break}+z_{1-a/2}\cdot 0.52\cdot \bar{x}\left(P_j\right)/\sqrt{N}\right\}=1-\alpha$(5)

In accordance with Equation 9(9) ${\gamma }_{DSC}=\frac{1}{N\left(N-1\right)}{\sum }_{t=1}^N{\sum }_{r=t+1}^{N}ex{p}^{j\left({\mathrm{\varnothing }}_{rt}-({\mathrm{\varnothing }}_{or}-{\mathrm{\varnothing }}_{ot}\right))}$(9) , the estimated interval of the average amplitude for reference pixel P_j is $\left[\bar{x}\left(P_j\right)-z_{1-a/2}\cdot 0.52\cdot \bar{x}\left(P_j\right)/\sqrt{N}\bar{x}\left(P_j\right)+z_{1-a/2}\cdot 0.52\cdot \text{break}\bar{x}\left(P_j\right)/\sqrt{N}\right]$. If the mean amplitude value of the pixel to be estimated falls within the average amplitude estimation interval, the pixel is regarded as a homogeneous pixel of reference pixel P_j. Similarly, an estimation window is set for each reference pixel, and the window’s center pixel is the reference pixel. Then, the estimated interval is constructed using the average amplitude of the reference pixel to estimate the remaining pixels in the window. All of the pixels in the SAR image are then traversed to obtain the corresponding homogeneous pixel set ${\mathrm{\Omega }}_{IEH{P}_j}$.

The merit of this algorithm is its high computational efficiency, with the experimental results showing greater reliability for small samples. When the sample number is large, if the average amplitude of the reference pixel is considered the true value and is used to construct a confidence interval, that constructed interval may be biased. This bias leads to a reduction in the number of selected homogeneous pixel samples and a reduction in the reliability of homogeneous pixels. Therefore, for large sample numbers, the size of estimation windows should not be fixed and the corresponding confidence interval should be dynamically updated to obtain reliable homogeneous pixel selection results.

Sequential selection algorithm for homogeneous pixels combining the BWS test algorithm with dynamic interval estimation theory

The method proposed in this study is divided into two steps. The first step uses the BWS test algorithm to select local homogeneous pixels, and then the initial confidence interval is constructed according to those selected pixels. The second step is to adjust the size of the estimation window for homogeneous pixels, synchronously update the confidence interval and then use the final confidence interval to obtain the set of homogeneous pixels. This method utilizes the advantages inherent in the BWS algorithm (i.e. this method is easily applied and is at least as powerful as the commonly used non-parametric tests). Additionally, this method is more sensitive to the differences between the tails of the sample and more efficient when examining small samples. It also has the high computational efficiency of interval estimation theory.

Initial local homogeneous pixels selected based on the BWS algorithm

According to Section 2.1, for all reference pixels, there is an estimation window, as shown in the black box in (Figure 1). The window’s size is assumed to be D×D. Additionally, for any reference pixel, a BWS test window is set, as shown in the blue box in (Figure 1), and this window’s size is T×T. In the BWS test window, for any reference pixel $P_j$(j = 1, 2… T²), m pixels are selected (1 ≤m≤T²-1). Equation 9(9) ${\gamma }_{DSC}=\frac{1}{N\left(N-1\right)}{\sum }_{t=1}^N{\sum }_{r=t+1}^{N}ex{p}^{j\left({\mathrm{\varnothing }}_{rt}-({\mathrm{\varnothing }}_{or}-{\mathrm{\varnothing }}_{ot}\right))}$(9) is used to obtain the test value B_m. Then, a Monte Carlo simulation experiment is conducted to obtain the numerical solution of the rejection domain, and all the pixels to be discriminated are subject to the estimation of homogeneous pixels. Finally, all homogeneous pixel sets ${\mathrm{\Omega }}_{BWSH{P}_j}$ of reference pixel $P_j$ in the BWS test window are obtained. The blue pixels in (Figure 1) represent the homogeneous pixel set obtained after the BWS test.

Figure 1. Schematic diagram of BWS-DIE algorithm.

(6) $\left\{\begin{array}{c}{B}_{P_j}=\frac{1}{2N^2}{\sum }_{i=1}^N\frac{{\left(R_{x_i\left(P_j\right)}-2i\right)}^2}{\frac{i}{N+1}\cdot \left(1-\frac{i}{N+1}\right)}\\ B_{P_{\left(m\right)}}=\frac{1}{2N^2}{\sum }_{i=1}^N\frac{{\left(R_{x_i\left(P_{\left(m\right)}\right)}-2i\right)}^2}{\frac{i}{N+1}\cdot \left(1-\frac{i}{N+1}\right)}\\ B_m=\frac{B_{P_j}+B_{P_{\left(m\right)}}}{2}\\ m=1,2,\dots ,T^2-1\end{array}\right.$(6)

In Equation 9(9) ${\gamma }_{DSC}=\frac{1}{N\left(N-1\right)}{\sum }_{t=1}^N{\sum }_{r=t+1}^{N}ex{p}^{j\left({\mathrm{\varnothing }}_{rt}-({\mathrm{\varnothing }}_{or}-{\mathrm{\varnothing }}_{ot}\right))}$(9) , B_m denotes the statistics for the testing of m pixel(s) to be estimated in the BWS test window.

Constructing homogeneous point sets using local homogeneous points and dynamic confidence interval estimation

After obtaining the local homogeneous pixel set ${\mathrm{\Omega }}_{BWSH{P}_j}$of reference pixel $P_j$ in the BWS test window, the initial confidence interval is constructed using the mean of the pixel amplitude of all homogeneous pixels in the homogeneous pixel set ${\mathrm{\Omega }}_{BWSH{P}_j}$ in accordance with Equation 9(9) ${\gamma }_{DSC}=\frac{1}{N\left(N-1\right)}{\sum }_{t=1}^N{\sum }_{r=t+1}^{N}ex{p}^{j\left({\mathrm{\varnothing }}_{rt}-({\mathrm{\varnothing }}_{or}-{\mathrm{\varnothing }}_{ot}\right))}$(9) (Jiang et al. 2014), as follows:

(7) $\left\{\begin{array}{c}P\left\{E_1^j-z_{1-a/2}\cdot 0.52\cdot E_1^j/\sqrt{N}<\mu \left(P_j\right)<E_1^j+z_{1-a/2}\cdot 0.52\cdot E_1^j/\sqrt{N}\right\}=1-\alpha \\ E_1^j=\frac{1}{Nk}\sum _{k=1}^{\mathrm{S}({\mathrm{\Omega }}_{BSH{P}_j)}}\sum _{i=1}^N{x}_i\left(P_k\right)\end{array}\right.$(7)

where $P\left\{.\right\}$ is the probability; x_i(P_k) is the amplitude corresponding to pixel P_k in the homogeneous pixel set in the i-th SAR image; $\mathrm{S}({\mathrm{\Omega }}_{BWSH{P}_j)}$ is the number of homogeneous pixels in the homogeneous pixel set; $E_1^j$ is the mean amplitude of all pixels in the homogeneous pixel set ${\mathrm{\Omega }}_{BWSH{P}_j}$, where superscript j references pixel P_j and subscript 1 is the mean amplitude when the confidence interval is constructed for the first time; and μ(P_j) is the mean intensity of the statistical population. Additionally, the initial confidence interval obtainedis $\left[E_1^j-z_{1-a/2}\cdot 0.52\cdot E_1^j/\sqrt{N}{E}_1^j+z_{1-a/2}\cdot 0.52\cdot E_1^j/\sqrt{N}\right]$.

Next, the BWS test window is expanded outwards by one pixel to obtain the estimation window after the first expansion, as shown in the green box in (Figure 1). Then, the initial confidence interval obtained is employed to discriminate all pixels in the estimation window after the first expansion. The homogeneous pixel set that fell into the confidence interval is called ${\mathrm{\Omega }}_{BWS-DIE}^1$, and the pixels in the set were the union of the blue and green pixels in (Figure 1). Next, a new confidence interval $\left[E_2^j-z_{1-\frac{a}{2}}\cdot 0.52\cdot E_2^j/\sqrt{N}{E}_2^j+z_{1-\frac{a}{2}}\cdot 0.52\cdot E_2^j/\sqrt{N}\right]$ is constructed using all homogeneous pixels in set ${\mathrm{\Omega }}_{BWS-DIE}^1$, and the BWS test window is expanded outwards by two pixels as a whole to obtain a new interval estimation window, as shown in the yellow box in (Figure 1). The confidence interval obtained from the second time is used to discriminate all of the pixels in the window following the latest expansion, and a new homogeneous pixel set ${\mathrm{\Omega }}_{BWS-DIE}^2$ is obtained.

The above steps are repeated until the BWS test window is expanded to a given initial estimation window size (the default window size is 15 × 15), at which point the final homogeneous pixel set ${\mathrm{\Omega }}_{BWS-DIE}^{D-T}$ of reference pixel P_j can be obtained. These steps are repeated to traverse all of the pixels in the SAR image, and the homogeneous pixel set ${\mathrm{\Omega }}_{BWS-DIE}$of all pixels in the SAR image can be obtained. As the proposed algorithm combines the BWS algorithm and dynamic interval estimation theory, the method is later called BWS-DIE. A BWS-DIE algorithm flow chart is shown in (Figure 2).

Figure 2. Flow chart of BWS-DIE algorithm.

Multitemporal InSAR based on BWS-DIE algorithm

The homogeneous pixel set ${\mathrm{\Omega }}_{BWS-DIE}$, obtained using the BWS-DIE algorithm proposed in this study, belongs to the candidate DS point set. To determine the final DS point, the following conditions were initially used for the preliminary selection:

(8) ${\mathrm{\Omega }}_{DSC}={\mathrm{\Omega }}_{BWS-DIE}\left(P_j\right)>25,i=1,............s^2$(8)

where ${\mathrm{\Omega }}_{DSC}$ is the initial DS point selection result, and ${\mathrm{\Omega }}_{BWS-DIE}\left(P_j\right)$ is the number of homogeneous pixels of any reference pixel P_j in the homogeneous pixel set. Here, the threshold for the number of homogeneous pixels was set to 25.

The surface features in the distributed target pixel are homogeneous, and the echo signal is formed by the coherent superposition of the sub scatterers in the pixel. The phase stability of the distributed target pixel is poor since it is affected by time and space decorrelation. Therefore, when the homogeneous pixel recognition step is completed, an operation of phase optimization is needed. Commonly used phase optimization methods include the phase coherence weighting method based on a set of homogeneous pixels (Jiang, Ding, and Li 2018), the eigenvalue decomposition method based on the covariance matrix (Fornaro et al. 2014; Cao, Lee, and Jung 2015b), the least square method for phase estimation (Cao, Lee, and Jung 2015a) and maximum likelihood estimation (Monti Guarnieri and Tebaldini 2008). In this study, the eigenvalue decomposition method, based on the covariance matrix, was employed to optimize the original phase. After phase optimization of the distributed target was complete, the fitting degree γ_DSC (temporal coherence) was determined between the interferometric phases before and after phase optimization (Fornaro et al. 2014; Cao, Lee, and Jung 2015a), as follows:

(9) ${\gamma }_{DSC}=\frac{1}{N\left(N-1\right)}{\sum }_{t=1}^N{\sum }_{r=t+1}^{N}ex{p}^{j\left({\mathrm{\varnothing }}_{rt}-({\mathrm{\varnothing }}_{or}-{\mathrm{\varnothing }}_{ot}\right))}$(9)

where γ_DSC is the fitting degree (also called temporal coherence); φ_rt is the interferometric phase generated by image r and image t before phase optimization; φ_or is the interferometric phase generated by image o and image r following phase optimization; and φ_ot is the interferometric phase generated by image o and image t after phase optimization.

Here, the temporal coherence threshold was set to 0.75. A distributed target with good coherence and a high signal-to-noise ratio could also be selected. In this study, to improve point density, the PS point target was selected using the amplitude deviation index method (Ferretti, Prati, and Rocca 2001). In addition, the selected PS point and DS point targets were integrated for three-dimensional phase unwrapping (Hooper et al. 2004) to retrieve the deformation rate and time series deformation of all point targets. Since this method created multiple temporal InSAR using the BWS-DIE algorithm, it was named BWS-DIE MT (multiple temporal) InSAR, and its processing flow is shown in (Figure 3).

Figure 3. BWS-DIE MT InSAR processing flow.

Experiments

Simulation experiment

The Monte Carlo method uses random numbers extracted from a population as samples for experiments and employs the obtained statistical eigenvalues (e.g. mean, probability and distribution) as the numerical solutions to the question to be solved (Ripley 2009). In this study, Monte Carlo random experiments were conducted to test and compare the performances of the KS test algorithm, BWS test algorithm, FaSHPS algorithm and BWS-DIE algorithm. In the experiment, the size of the given grid was $15\times 15$, where the $15\times 8$ pixel was assigned to the homogeneous pixel $K_1$ and the pixels of the remaining $15\times 7$ were assigned to the heterogeneous pixel $K_2$. Additionally, a random number of Rayleigh distributions were added as noise to the sample populations of $K_1$ and $K_2$. In theory, the heterogeneous pixel percentage (HPP) can be expressed as (Jiang, Ding, and Li 2018) follows:

(10) $HPP=\left(\left(15\times 7+\alpha \times 15\times 8\right)/225\right)\times 100$(10)

The confidence level $\alpha$ in Equation 9(9) ${\gamma }_{DSC}=\frac{1}{N\left(N-1\right)}{\sum }_{t=1}^N{\sum }_{r=t+1}^{N}ex{p}^{j\left({\mathrm{\varnothing }}_{rt}-({\mathrm{\varnothing }}_{or}-{\mathrm{\varnothing }}_{ot}\right))}$(9) was always set to be 0.05, and the true power obtained was 49.3%. All of the pixels in the grid were tested one by one using the KS algorithm, BWS algorithm, FaSHPS algorithm and BWS-DIE algorithm, respectively. To avoid any influence from random errors, each algorithm was repeated 10,000 times during the experiment. As the ratio of K₁/K₂ increased, the power value of the alternative hypothesis at different contrast levels could be obtained (Jiang, Ding, and Li 2018; Jiang et al. 2017). To compare the performances of different homogeneous pixel selection algorithms for different sample numbers, the sample number M was increased from 10 to 60. The above experiment was then repeated and the final alternative hypothesis power values were obtained, as shown in (Figure 4) show the rejection values of the alternative hypotheses obtained using different homogeneous pixel selection algorithms for different sample numbers. In the figure, the upper and lower boundaries of the box plot are the upper quartile and the lower quartile (25% to 75%) of the 10,000 experimental results, respectively. The internal solid line is the median line, and the upper and lower lines parallel to the median line outside the box plot are outlier cutoff points. The selected cutoff values are [$Q_3$ +1.5IQR,$Q_3$-1.5IQR], where $Q_3$ and $Q_1$ are the upper and lower quartiles, respectively; IQR=$Q_3$-$Q_1$; the solid points inside the rectangular box are average values; and the internal dashed lines are the true values of the power.

Figure 4. Also show that, when the sample numbers differ, 4(a-f) as the contrast of K₁/K₂ increases, and when the difference between homogeneous pixels and heterogeneous pixels becomes increasingly obvious, the rejection means for all of the algorithms can converge to the vicinity of the real power. However, when the sample numbers differ, each algorithm performs differently. In (Figure 4(a-f)), when the contrast of K₁/K₂ is 1–1.4 (the power value does not reach the convergence state), the rejection mean of the KS algorithm is lower than the rejection means of the BWS, FaSHPS and BWS-DIE algorithms. This finding indicates that, when the difference between homogeneous pixels and heterogeneous pixels is not significant, the KS algorithm is at its lowest performance level. The BWS algorithm slightly outperforms the KS algorithm, but (Figure 4(a)) clearly indicates that its performance is relatively poor for small samples and that it is close to the real power value until the contrast of K₁/K₂ exceeds 2.2. This finding indicates that the BWS algorithm increases the type II error for small samples and is prone to contain a greater number of heterogeneous samples.

The rejection means of the FaSHPS and BWS-DIE algorithms are close and superior to those of the KS and BWS algorithms. Additionally, (Figures 4(a)) show that, when the rejection mean does not reach the convergence state, the rejection mean of the BWS-DIE algorithm is higher than that of the FaSHPS algorithm. This finding indicates that the results obtained by employing the BWS-DIE algorithm are more reliable. Compared with the rejection mean of the FaSHPS algorithm, the rejection mean of the BWS-DIE algorithm is closer to the true power value when the rejection mean reaches the convergence state; this finding further reveals that the results obtained by employing the BWS-DIE algorithm are more accurate. In addition, for further comparative analysis in this study, statistics on the standard deviation (STD) of the average rejection (convergence state) of the above four algorithms under six sample conditions and the runtime of the algorithms were collected. The statistical results are shown in (Table 1).

Table 1. Comparison of performance results for different homogeneous pixel algorithms. Note that all data were collected in the convergence state of power value (i.e. the contrast of K1/K2 is 3).

Algorithm	Sample Number	10	20	30	40	50	60
BWS-DIE	STD (pixel)	0.017	0.014	0.014	0.014	0.013	0.012
	Runtime (s)	9.6	14.7	19.7	28.4	32.2	37.3
KS	STD (pixel)	0.05	0.034	0.033	0.045	0.036	0.041
	Runtime (s)	464.4	407.1	503.5	506.8	549.8	601.6
BWS	STD (pixel)	0.043	0.047	0.045	0.047	0.045	0.047
	Runtime (s)	315.5	378.8	407.3	455.2	467.4	481.2
FaSHPS	STD (pixel)	0.024	0.019	0.018	0.018	0.017	0.017
	Runtime (s)	9.2	14.1	19.9	26.4	32.9	37.1

(Table 1) shows that the BWS-DIE algorithm has the highest precision when selecting homogeneous pixels for different samples. Compared to the KS, BWS and FaSHPS algorithms, the BWS-DIE algorithm increases average precision by 64.3%, 69.4% and 25.3%, respectively. Additionally, when the sample number is 60, the slowest runtime of the KS algorithm is 601.6 s, while at 481.2 s the BWS algorithm has the second-slowest runtime. The runtime of the FaSHPS algorithm is 37.1 s, which is essentially the same as the 37.3 s runtime of the BWS-DIE algorithm.

The above analysis indicates that, compared to the KS test algorithm, BWS test algorithm and FaSHPS algorithm, the BWS-DIE algorithm proposed in this study effectively reduces the homogeneous pixel selection error caused by the hypothesis test and greatly improves the reliability of homogeneous pixel recognition results. In addition, the proposed algorithm is insensitive to the number of samples, is highly efficient and can be applied to additional scenes.

Case study

Study area

To further verify the superiority and reliability of the proposed algorithm, this study used the coal mining subsidence area in western Xuzhou, Jiangsu Province, China as an experimental area, wherein the BWS-DIE MT InSAR method proposed in this study was employed to conduct deformation monitoring. Located in north-western Jiangsu Province, China, Xuzhou (Figure 5(a)) has a coal mining history exceeding 120 years, with a maximum annual output of approximately 25 million tons. The research area is located in the Tongshan and Quanshan Districts of Xuzhou. The terrain in this area is relatively flat, most of the ground vegetation is crops and the village represents the primary grouping of buildings. The research area includes three mines; namely, the Jiahe, Zhangxiaolou and Pangzhuang mines, covering the approximate area 117° E − 117°9“E, 34°18” − 34°22” (see Figure 5(b)). The mining information of the working panel is shown in (Table 2).

Figure 5. Location and mine distribution in research area, as modified from Chen et al. (2022) [remote sensing]. (a) Location map of research area. (b) Mine boundaries and distributions of some working panels. Red triangle represents Xuzhou – the research area. Red dots represent fourth-class leveling monitoring points for a total of ten points. Note that leveling data were collected using an OPTONS3 electronic level with accuracy of ±3 mm per kilometer round trip. (c), (d) and (e) are enlarged drawings of working panels of Zhangxiaolou, Jiahe and Pangzhuang mines, respectively..

Table 2. Working panel information.

Working panel no.	Mining period	Average mining thickness (m)	Average mining depth (m)
No. 1	2010/05–2010/11	1.89	1101
No. 2	2011/01–2011/10	1.89	1101
No. 3	2007/12–2008/07	1.74	1145
No. 4	2010/04–2011/02	1.90	1072
No. 5	2008/04–2008/10	1.80	1140
No. 6	2006/11–2007/06	1.86	1115
No. 7	2009/06–2010/01	2.20	1121
No. 8	2009/08–2011/05	1.80	1043
No. 9	2010/03–2010/07	1.75	1025
No. 10	2007/07–2008/03	1.8	1089
No. 11	2008/06–2009/07	2.2	1130
No. 12	2007/01–2007/08	1.9	1120
No. 13	2007/09–2007/12	1.86	1075
No. 14	2008/12–2009/03	2.85	480
No. 15	2007/8–2008/12	2.95	440
No. 16	2007/01–2007/07	3.41	448
No. 17	2007/08–2008/05	3	350
No. 18	2009/06–2010/05	3.2	475
No. 19	2010/03–2010/11	2.4	390
No. 20	2010/07–2011/08	2.1	370
No. 21	2009/06–2010/09	2.3	550
No. 22	2007/08–2008/12	2.5	585
No. 23	2009/04–2009/10	3.1	510
No. 24	2010/09–2011/04	3	430
No. 25	2010/05–2010/11	1.88	400
No. 26	2010/07–2011/11	2.8	418
No. 27	2008/03–2009/09	2.8	379
No. 28	2008/09–2009/05	2.8	427
No. 29	2007/10–2008/09	2.3	448
No. 30	2009/10–2010/09	2.5	547
No. 31	2009/10–2010/03	2.8	373
No. 32	2009/10–2010/03	2.42	391
No. 33	2009/10–2010/03	2.1	201
No. 34	2008/05–2009/04	2.1	326
No. 35	2007/08–2008/05	2.7	350
No. 36	2010/01–2010/05	2.2	340

The Jiahe mine is approximately 5.5 km long from east to west and 4.5 km wide from north to south, the Zhangxiaolou mine is approximately 5.0 km long from east to west and 3.4 km wide from north to south, and the Pangzhuang mine is approximately 6.1 km long from east to west and 3 km wide from north to south; resulting in a total area of approximately 60 km². In the previous literature, Fan et al. (2015) used conventional differential InSAR technology combined with the probability integration method to obtain surface subsidence information in the area during mining. However, they obtained only the surface deformation laws during part of the coal mining stage and the precision of the obtained monitoring results was limited due to the influence of spatial-temporal decorrelation. Additionally, Chen et al. (2022) employed SBAS InSAR technology to obtain the multidimensional spatial-temporal evolution law of the surface subsequent to the closure of the mines in the area but did not observe the spatial-temporal evolution law of surface deformation during mining. This is noteworthy, since the dynamic mining phase is often the stage at which surface deformation is most intense, and deformation at this time is more hazardous to surface buildings. Therefore, in this study, BWS-DIE MT InSAR technology was employed to retrieve the surface spatial-temporal evolution law during mining in an effort to provide a reference for disaster prevention and mitigation in mining areas.

Data

Twelve scenes of ALOS-1 PALSAR single-look complex (SLC) images covering the research area were selected. The resolutions of the range and azimuth directions were 4.68 m and 3.17 m, respectively, and the time span was from 20 February 2007 to 3 March 2011. The specific satellite parameters are shown in (Table 3) and the reference image acquisition date appears in bold red text. Notice that the incidence angle of radar is 38.72°with a fine beam dual polarization mode (HH and HV). The digital elevation model (DEM) used in the experiment was Shuttle Radar Topography Mission (SRTM) DEM 1 provided by NASA with a spatial resolution of 30 m × 30 m [NASA, 2015. NASA SRTM Version 3.0 Global 1 Arc Second Dataset. https://lpdaac.usgs.gov/documents/179/SRTM_User_Guide_V3.pdf]. The DEM was used to remove the topographic phase from the interferometric phase.

Table 3. SAR image parameters. The perpendicular and temporal baselines are calculated based on the reference image.

Date	Perpendicular baseline (m)	Temporal baseline (days)
2007/02/20	−1201	690
2008/01/08	−2632	368
2008/02/23	−3289	322
2009/01/10	0	0
2009/02/25	171	46
2009/10/13	1446	276
2010/01/13	1802	368
2010/02/28	2381	414
2010/05/31	2639	506
2010/10/16	3291	644
2011/01/16	3728	736
2011/03/03	4188	782

Experimental process and analysis of results

To mitigate the impact of spatial-temporal baseline decorrelation, the SAR image obtained on 10 January 2009 was set as the reference image in the experiment, and the spatial and temporal baseline thresholds were set to 300 d and 5000 m, respectively. A total of 24 interferometric pairs was obtained and the spatial-temporal baseline distribution of the interferometric pairs formed by the images is shown in (Figure 6). Then, all of the SAR images were registered to the reference image and the two-pass difference method was used to obtain the corresponding differential interferograms.

Figure 6. Spatial-temporal baseline distributions of interferometric pairs. The red dot represents the acquisition data of reference image.

Based on the co-registered time series SAR images, the corresponding intensity dataset was acquired. The initial estimation window size of the homogeneous pixel was set to $15\times 15$ and the BWS-DIE algorithm proposed in Section 2.3 was employed to obtain the set of homogeneous pixels. Thereafter, the eigenvalue decomposition method (Fornaro et al. 2014; Cao, Lee, and Jung 2015b), based on the covariance matrix, was applied to optimize the original phase of the homogenous pixel set, and the final DS point was determined according to Equations (8)(8) ${\mathrm{\Omega }}_{DSC}={\mathrm{\Omega }}_{BWS-DIE}\left(P_j\right)>25,i=1,............s^2$(8) and (9(9) ${\gamma }_{DSC}=\frac{1}{N\left(N-1\right)}{\sum }_{t=1}^N{\sum }_{r=t+1}^{N}ex{p}^{j\left({\mathrm{\varnothing }}_{rt}-({\mathrm{\varnothing }}_{or}-{\mathrm{\varnothing }}_{ot}\right))}$(9) ). Additionally, the amplitude dispersion index was used to select PS points in this area and 579,503 point targets (PS+DS) were obtained. Then, the annual average deformation rate of the research area was obtained using the BWS-DIE MT InSAR method proposed in Section 2.4, as shown in (Figure 7(a)). To compare the difference between the different methods, the StaMPS (Hooper et al. 2004; Hooper, Segall, and Zebker 2007) method was also applied to retrieve the annual average deformation rate of the research area, as shown in (Figure 7(b)).

Figure 7. Average annual subsidence velocity in line-of-sight (LOS) direction in research area. (a) Average annual subsidence velocity obtained using BWS-DIE MT InSAR method. (b) Average annual subsidence velocity obtained using StaMPS method. Note that in (a) and (b), white circles A, B, C, D, E, F, G and H indicate areas with drastic deformation; significant differences in results were obtained using both methods.

According to the results in (Figure 7(a)), three areas in the research area – namely, Areas A, C and E – exhibited high subsidence velocities. Among these areas, Area C in (Figure 7(a)) exhibited the highest subsidence velocity of 305 mm/yr, with Area A having the next highest at a maximum subsidence velocity of 280 mm/yr. Compared with the surface subsidence of Areas A and C, that for Area E was relatively gentle, with a maximum subsidence velocity of 195 mm/yr. As (Figure 5) and (Table 2) show, working panel Nos. 5–8, 6–30 and 4–15 were being mined directly below subsidence Areas A, C and E, respectively. As the DS points in these three areas were relatively dense, they reflect the spatial-temporal evolution patterns of the surface deformation above the working panel while coal was being mined. Compared with the results shown in (Figure 7(a)), only 85,915 point targets were obtained using the StaMPS method, as shown in (Figure 7(b)). The point density was reduced by 85.2% primarily because the farmland vegetation above the mining area was dense and the point selection approach of the StaMPS method is based on the principle of high coherence. As a result, many DS points were missed.

According to the results in (Figure 7(b)), from 20 February 2007 to 3 March 2011, the highest surface subsidence velocity was only 265 mm/yr. This occurred in Area C and represents a decrease of 13.1% compared with the monitoring results of BWS-DIE MT InSAR. Additionally, (Figure 7(b)) indicates that the number of monitoring points in subsidence Areas A, C and E was seriously insufficient, and only a few points in the locations yielded any drastic deformation. Therefore, the mining subsidence characteristics in the area could not be fully and accurately retrieved. Areas B, D, F, G and H, shown in (Figure 7), indicate areas where the results acquired by employing the BWS-DIE MT InSAR and StaMPS methods differed greatly, and (Figure 7(a)) reveals that Areas B, D, F, G and H all exhibited significant subsidence. However, only a few monitoring points were revealed in Areas B, D, F, G and H, which are shown in (Figure 7(a)), and no obvious subsidence occurred.

To quantify the pattern of surface deformation under mining conditions, time-dimensional integration was performed on the deformation rate to obtain the cumulative surface deformation at each mine, as shown in (Figures 8 and 9). According to (Figures 8(a-i)), and commencing on 20 February 2007, as underground mining activities progressed, the subsidence values in Areas A, C and E continuously increased and the subsidence range expanded. As of 3 March 2011, maximum surface subsidence reached −1303 mm in Area C of the Jiahe mine (as shown in Figure 8(i)); from 23 February 2008, subsidence began and continued to increase in Areas B and D. From 13 October 2009, subsidence began and continued to increase in Areas F and G, and as of 3 March 2011, subsidence values in Areas B, D, F and G have been increasing.

Figure 8. Time series cumulative subsidence obtained by BWS-DIE algorithm for research area in LOS direction. A, B and C indicate areas with significant subsidence.

Figure 9. Time series cumulative subsidence obtained by StaMPS method for research area in LOS direction. A, B and C indicate areas with significant subsidence.

According to (Figure 8), there were no working panels under Areas G and H but there were varying degrees of subsidence. Historical optical images show collapse lakes formed by coal mining collapses next to Areas G and H, and these lakes tend to expand. Therefore, it is speculated that the collapsed lake expansion may cause instability in nearby underground internal rock masses, resulting in surface subsidence in Areas G and H.

According to (Figures 9(a-i)), from 20 February 2007 to 3 March 2011, maximum surface subsidence was only −1160 mm, which is less than that obtained using the BWSDIE algorithm. In addition, (Figures 9(a-i)) show that point density was extremely sparse in Areas A, C and E, where surface deformation was the most drastic. In addition, Areas B, D, F and G showed no obvious deformation. Compared with the information presented in (Figures 8(a-i)), much important information is lost in (Figure 9), and it cannot fully and correctly retrieve the evolution patterns of ground deformation during the underground mining period.

Based on the above analysis, 579,503 points were obtained based on the BWS-DIE MT InSAR method. Compared with StaMPS, which obtained 85,915 point targets, this represents an increased point density of 5.7 times. The monitoring results also suggest that, compared to the StaMPS algorithm, the BWS-DIE algorithm improved the ability to monitor maximum subsidence by approximately 12.3%.

To verify the reliability of the monitoring results obtained using the BWS-DIE MT InSAR algorithm, the monitoring results acquired using the proposed method and the StaMPS algorithm at the same position as the leveling points in (Figure 5(b)) were extracted. To facilitate a comparison with the leveling monitoring results, the deformation results in the LOS direction were projected in the vertical direction, based on the SAR orbit geometric relationship (Chen et al. 2020). For time consistency, the leveling data were interpolated into the same acquisition period as the SAR image data (20 February 2007 to 16 October 2010) through spline interpolation.

A comparison of the results obtained using the BWS-DIE MT InSAR and StaMPS algorithms, and leveling measurement values, is shown in (Figure 10(a)). This figure indicates that the results obtained using the BWS-DIE MT InSAR algorithm are relatively close to the leveling value, while the results obtained using the StaMPS algorithm largely deviate from the leveling value at points 1 to 5. Additionally, using the leveling monitoring results as reference values, the RMSE and correlation coefficients of the results obtained using the BWS-DIE MT InSAR and StaMPS algorithms were calculated, and the specific results are shown in (Figures 10(a)). These figures further show the difference in precision between the two methods: the RMSE based on the BWS-DIE MT InSAR algorithm is 32 mm while that obtained using the StaMPS algorithm is 64 mm. Compared to the StaMPS algorithm, the BWS-DIE MT InSAR algorithm increases monitoring precision by 50%. Additionally, the correlation coefficient of the BWS-DIE MT InSAR algorithm is 0.93, while the correlation coefficient of the StaMPS algorithm is only 0.74; this finding further illustrates the reliability of the BWS-DIE MT InSAR algorithm.

Figure 10. Comparisons of different monitoring results. (a) Comparison of subsidence values obtained using BWS-DIE MT InSAR and StaMPS algorithms with fourth-order leveling results. (b) RMSE and correlation coefficients calculated from subsidence results and leveling data obtained using BWS-DIE MT InSAR algorithm. (c) RMSE and correlation coefficients calculated from subsidence results and leveling data obtained using StaMPS algorithm.

The main reason for the above results is that the point density of the results gained using the StaMPS algorithm was very sparse and many key points were missing. During verification, the deformation value of the corresponding point was primarily obtained by interpolation so that the obtained result deviated too much from the actual deformation; this result also demonstrates the importance of point density for accurately retrieving the evolution law of surface deformation.

Discussion

Results comparison and analysis of DS points selected using different homogeneous pixel algorithms

The above twelve scenes of SAR data were used as the test dataset, and representative homogeneous pixel selection algorithms (including the KS, BWS and FaSHPS algorithms) and the BWS-DIE algorithm proposed in this study were employed to conduct homogeneous pixel selection experiments on SAR datasets. For the KS, BWS and FaSHPS algorithms, the size of the homogeneous pixel estimation window was defined as $15\times 15$. For the BWS-DIE algorithm, the size of the initial estimation window was defined as $15\times 15$. The size of the BWS test window was defined as $7\times 7$, and the confidence level parameter $\alpha$was uniformly set to 0.05.

Generally, sparse vegetation and bare land contain many homogeneous pixels, while buildings and structures are mostly PS targets and contain fewer homogeneous pixels. Therefore, SAR intensity images can be employed to preliminarily identify areas containing many homogeneous pixels, thus assisting in verifying the accuracy of the homogeneous pixel estimation results. The SAR intensity image of the research area is shown in (Figure 11); most of the area covered by yellow circle I in (Figure 11) is a collection of dense, widely scattered buildings that are therefore mostly PS targets. Theoretically, the number of homogeneous pixels selected in this area is relatively small. In (Figure 11), the area covered by yellow circle K is mostly bare land, farmland and sparse vegetation, which are typical DS targets. The red dots J, L and M are located at the feature points of the buildings, roads and vegetation, respectively. (Figures 12(a-d)) show the results of selecting homogeneous pixels using the BWSDIE, KS, BWS and FaSHPS algorithms, respectively. The value of each pixel in the figure represents the number of homogeneous pixels selected for the reference pixels at the same position in the SAR image. (Figures 12(a-d)) show that, if the average number of homogeneous pixels from among all pixels is sorted, the number of homogeneous pixels selected using the KS algorithm is the largest, and the average number of homogeneous pixels for each reference pixel is 186. This value is followed by 166, 157 and 151, which are the average numbers obtained using the BWS, BWS-DIE and FaSHPS algorithms, respectively. In (Figure 12(a)), white circle I indicates a densely built area, and the selected results show few homogeneous pixels in this area. In contrast, the area overlain by white circle K is covered by vegetation and roads. The selected results show that the distribution difference of the homogeneous pixels between roads and vegetation areas is clear.

Figure 11. SAR intensity image of research area. Yellow circle K indicates typical area covered by bare land, farmland and sparse vegetation, and yellow circle I indicates area where dense buildings are located. J, L and M are feature points of buildings, roads and vegetation, respectively.

Figure 12. Comparison of results of selection numbers obtained using different homogeneous pixel algorithms. (a-d) are homogeneous pixel selection results obtained using BWS-DIE, KS, BWS and FaSHPS algorithms, respectively. Areas covered by white circles I and K are buildings and vegetation coverage areas, respectively. Color scale means number of homogeneous pixels for each pixel.

In (Figure 12(b)), the area covered by white circle K shows that roads and vegetation are not clearly distinguishable, and the selection results are partially incorrect. The area covered by white circle I should have few homogeneous pixels whereas the selection results contain many, which indicates that, although the KS algorithm selected the most homogeneous pixels, the results were unreliable and contained many heterogeneous pixels. Compared with (Figure 12(b)), the results of areas covered by white circles I and K in (Figure 12(c)) show that homogeneous pixels can be better discriminated and that some heterogeneous pixels can be eliminated. Compared to the areas covered by white circles I and K in (Figures 12(a)), the overall outline of homogeneous pixels is not sufficiently clear. This finding indicates that the results of the homogeneous pixel selection still contained heterogeneous pixels. In (Figure 12(a)), areas covered by white circles I and K have very similar results. The construction area covered by white circle I has very few homogeneous pixels, and the vegetation and roads in the area covered by white circle K are clearly visible.

To present a clearer illustration of the homogeneous pixel estimation results of the BWS-DIE, KS, BWS and FaSHPS algorithms, feature points J, L and M in (Figure 11) were selected as reference pixels in the experiment, and the homogeneous pixel estimation results for different reference pixels were obtained using the above algorithms, as shown in (Figure 13(a)). These figures show that the number of homogeneous points gained using the BWS-DIE algorithm in the building area is the lowest because the buildings are mostly permanent scatterers, and the number of homogeneous pixels should theoretically be close to zero.

Figure 13. Results comparing the number of homogeneous pixels according to different homogeneous pixel selection methods and different reference pixels. (a-d) show number of homogeneous pixels obtained using BWS-DIS, KS, BWS and FaSHPS algorithms, respectively, when reference pixel is J. Red dots indicate reference pixels, and the pixel is located in building area. Yellow dots indicate the results of selecting homogeneous pixels. (e-h) are number of homogeneous pixel samples obtained using BWS-DIS, KS, BWS and FaSHPS algorithms when reference pixel is L. Red dots indicate reference pixels and the pixel being located on a road; yellow dots indicate results of homogeneous pixel selection. (i-l) show number of homogeneous pixel samples obtained using BWS-DIS, KS, BWS and FaSHPS algorithms, respectively, when reference pixel is M. Red dots indicate reference pixels, and pixel is located in vegetation area. Yellow dots indicate results from homogeneous pixel selection.

The selection results obtained using the BWS-DIE algorithm matched the actual situation. However, the results of the KS algorithm (Figure 13(b)) contain many heterogeneous pixels and its findings are therefore not reliable. In (Figure 13(c)), the homogeneous pixels are relatively close to each other, and the selection exhibits some errors. (Figures 13) show that, when the estimation window contains roads and vegetation, the selection results in 15 (Figure 13(e)) can effectively distinguish between roads and vegetation coverage areas. The results of the KS test (Figure 13(f)) and BWS test (Figure 13(g)) show that the homogeneous pixel set contains many heterogeneous samples; compared to the BWSDIE algorithm (Figure 13(e)), the FaSHPS test result (Figure 13(h)) is not sufficiently smooth in terms of homogeneous pixels, indicating that some homogeneous pixels were missed. (Figures 13(i-l)) show that, in vegetation areas with more distributed targets, the performances of the BWS-DIE and FaSHPS algorithms are equivalent (Figures 13(i-l)) and the KS algorithm selects the fewest homogeneous pixels (Figure 13(j)).

In summary, it can be inferred that, although the KS and BWS algorithms selected many homogeneous pixels in some scenes (such as building areas), the homogeneous pixels selected by them contained many heterogeneous pixels and the results were unreliable. Additionally, in different scenes (including buildings, vegetation and roads), the homogeneous pixel results obtained using the BWS-DIE algorithm were most consistent with the actual scene.

The above experimental results were consistent with the results from the Monte Carlo simulation experiments in Section 3.1 (Figure 4(a)). That is, when small samples were tested, the mean rejection from the KS test was the lowest, while that from the BWS test was the next lowest. The STDs of rejection for these two test methods were much higher than those for the FaSHPS and BWS-DIE algorithms. Therefore, the test results from the KS and BWS algorithms contained many heterogeneous pixels and, as a result, the overall profile did not look sufficiently clear. In contrast, the MRE and STD of the BWS-DIE algorithm’s rejection mean were the lowest compared to those of the FaSHPS, BWS and KS algorithms. Therefore, the BWS-DIE algorithm could accurately identify homogeneous pixels and eliminate heterogeneous pixels, and (Figure 12(a)) looks clearer and smoother. The results of (Figure 12 and 13) also further illustrate that the BWS-DIE algorithm can effectively control type I and type II errors and contain all homogeneous pixels to the maximum extent. Compared to the KS, BWS and FaSHPS algorithms, the BWS-DIE algorithm obtained more robust results.

Comparative analysis of deformation monitoring results obtained using different homogeneous pixel selection algorithms

To investigate the influence of the aforementioned homogeneous pixel selection algorithms on the final deformation monitoring results, DS points were obtained using the homogeneous pixel algorithms mentioned in Section 4.1 (including the KS, BWS and FaSHPS algorithms) and then integrated into the StaMPS processing flow (Hooper et al. 2004; Hooper, Segall, and Zebker 2007). As a result, the annual average subsidence velocities of the research area were retrieved using different algorithms. As shown in (Figures 14(a-d)), white circle N indicates areas with the largest annual average subsidence velocities obtained by the algorithms for the research area along the LOS direction.

Figure 14. Annual average subsidence velocities obtained by different homogeneous pixel algorithms for research area along LOS direction. (a-d) show annual average subsidence velocities obtained by BWS-DIE, KS, BWS and FaSHPS algorithms for research area, respectively. White circle N indicates area with largest annual average velocity obtained by different homogeneous pixel algorithms.

(Figure 14) indicates that the research area’s overall deformation trend, obtained using different homogeneous pixel selection algorithms, fundamentally remains similar, but the differences in deformation velocity are distinct. For example, the maximum annual average subsidence velocity obtained using the BWS-DIE algorithm is 317 mm/yr (area N in Figure 14(a)), the maximum annual average subsidence velocity gained using the KS algorithm is 270 mm/yr (area N in Figure 14(b)), the maximum annual average subsidence velocity acquired using the BWS algorithm is 281 mm/yr (area N in Figure 14(c)) and the maximum annual average subsidence velocity obtained using the FaSHPS algorithm is 263 mm/yr (area N in Figure 14(d)). Additionally, (Figures 15(a-c)) were plotted to show the deformation velocities of the common points contained in the results obtained using the BWS-DIE and other algorithms. The blue dot indicates the smallest point density and the red dot indicates the largest. The correlation coefficients between the BWS-DIE and other algorithms and corresponding RMSEs were also calculated. (Figure 15(a)) reveals that the correlation between the results obtained using the BWS-DIE and KS algorithms is 0.91, and the RMSE is 14.4 mm. (Figure 15(b)) shows that the correlation between the results obtained using the BWS-DIE and BWS algorithms and the RMSE is superior to that obtained using the KS algorithm. Finally, (Figure 15(c)) shows that the correlation coefficient between the results of the BWS-DIE and FaSHPS algorithms is the highest, and the RMSE is 9.1 mm. The above results show that the common point results of the BWS-DIE algorithm were very similar to those of the KS, BWS and FaSHPS algorithms.

Figure 15. Average annual subsidence velocities and densities of common points of results of different algorithms. (a) LOS displacement velocity derived using KS algorithm versus BWS-DIE algorithm. (b) LOS displacement velocity derived using BWS algorithm versus BWS-DIE algorithm. (c) LOS displacement velocity derived using FaSHPS algorithm versus BWS-DIE algorithm. Results are colour-coded according to pixels’ smoothed density (Eilers and Goeman 2004), in which red represents the greatest dot density and blue represents the least dot density.

(Figure 16) shows histograms of the deformation velocities of all pixels in the results of the BWS-DIE, KS, BWS and FaSHPS algorithms. (Figure 16(a)) indicates that most of the results from the KS algorithm are between −250 mm/yr and 0 mm/yr, (Figure 16(b)) reveals that most of the results from the BWS algorithm are between −270 mm/yr and 0 mm/yr, and (Figure 16(c)) shows that most of the results from the FaSHPS algorithm are between 250 mm/yr and 0 mm/yr. In comparison, the results obtained using the BWS-DIE algorithm have a wider deformation velocity range, and the number of pixels of the BWSDIE algorithm was higher than that of other algorithms for most deformation velocity intervals. Notably, the BWS-DIE algorithm monitored deformation velocity between 310 mm/yr and −270 mm/yr, which results went unmonitored by the other algorithms. In short, with the BWS-DIE algorithm, more DS points were obtained, the deformation velocity range was wider, and new deformation patterns were found.

Figure 16. Histogram of deformation velocities derived by different methods. (a) BWSDIE algorithm versus KS algorithm. (b) BWS-DIE algorithm versus BWS algorithm. (c) BWS-DIE algorithm versus FaSHPS algorithm.

To retrieve the temporal-spatial evolution of the surface, time-dimension integration was performed on the deformation rate to obtain the cumulative surface deformation of each mine, as shown in (Figures 17(a-l)). For the sake of simplicity, the only deformation monitoring results provided in this study are from 20 February 2007 to January 10, 20 February 2009, 2007 to 28 February 2010 and 20 February 2007 to 3 March 2011. For the remaining deformation results, please see the supplementary file. (Figures 17(a-l)) show that the subsidence trends obtained using different algorithms are fundamentally similar. However, the subsidence results obtained by different homogeneous pixel algorithms are different in terms of maximum subsidence and local deformation patterns, as shown in the areas (Figure 17) covered by red circles N, O and P. 9. Figure 17(a)) shows the subsidence results from 20 February 2007 to 10 January 2009. During this stage, working panel Nos. 3, 5, 6, 10–17, 22, 27–29 and 34–35 began to be mined, and surface subsidence of varying degrees began to occur in the three mining areas. (Figure 17(c)) indicates that the maximum subsidence obtained by the BWS-DIE algorithm was −1382 mm and that it occurred (area N) between 20 February 2007 and 3 March 2011. This result is an increase of 19.3%, 20.6% and 17.3% when compared to the maximum subsidence values obtained by the KS (Figure 16(f)), BWS (Figure 16(l)) and FaSHPS (Figure 16(l)) algorithms, respectively.

Figure 17. Cumulative subsidence of research area for period February 20, 2007 to March 3, 2011, when different homogeneous pixel algorithms were employed. (a-c) show cumulative subsidence obtained by BWS-DIE algorithm; (d-f) show cumulative subsidence obtained by KS algorithm; (g-i) show cumulative subsidence obtained by BWS algorithm; and (j-l) show cumulative subsidence obtained by FaSHPS algorithm. Red circles O and P indicate areas with different local surface deformation characteristics obtained by algorithms from February 20, 2007 to January 10, 2009, and red circle N indicates area with largest surface subsidence value.

The areas covered by red circles O and P in (Figures 17(a,d,g,j)) show that the local surface deformation patterns retrieved by the algorithms are different from 20 February 2007 to 10 January 2009. The areas covered by these circles were also enlarged, as shown in (Figure 18 and 19). The results indicate that, with advancements in mining, the surface subsidence caused by working panel No. 15 continues to expand. However, the results obtained by the KS (Figures 18(d-i)) and FaSHPS (Figures 18(f-l)) algorithms show that the subsidence above working panel No. 15 is expanding no further, and thus some zero values appear in the results of the BWS algorithm (Figures 18(g-h)).

Figure 18. Time series cumulative subsidence results of areas covered by red dotted circle O from February 20, 2007 to January 10, 2009, when different algorithms were employed. Black polygon represents working panel No. 15..

Figure 19. Cumulative subsidence results of areas covered by red dotted circle P from February 20, 2007 to January 10, 2009, when different algorithms were employed.

(Figures 19(a-c)) indicate an ongoing subsidence trend in the north-eastern part of the collapsed lake from 20 February 2007 to 10 January 2009. (Figures 20(a-b)), while the results of other algorithms do not show subsidence trend in this area. (Figure 18 and 19) indicate that the BWS-DIE algorithm can provide additional details to facilitate the accurate retrieval of the temporal-spatial evolution of the surface, through which a greater number of deformation patterns can be discovered.

Figure 20. The Google images between March 2017 and October 2012. The irregular polygon formed by the blue curved line represents the area of the lake.

(Figure 21(a-h)) shows the time series comparison of the deformation monitoring results obtained by different algorithms and the leveling monitoring results, respectively. It is notice that all the monitoring results are interpolated using Weibull function (Yang, Kozak, and Smith 1978) to align in time with the leveling data. To quantitatively evaluate the monitoring precision obtained by the different algorithms, the RMSE and spatial correlation coefficients for each algorithm were calculated. In addition, statistics were collected on the number of DS points gained by each algorithm, the results of which are presented in (Figures 21(i-l)). (Figure 21(i)) shows that the BWS-DIE algorithm retrieved 546,727 DS points and an RMSE of 32 mm, and that it also gained the largest corresponding spatial correlation coefficient of 0.93. Compared to the performance of the other algorithms, the performance of the BWE-DIE algorithm on each index is optimal, indicating reliable performance. Although the KS algorithm selected a greater number of DS points (510,732), the results largely deviated from the leveling measurements (Figure 21(j)). Its RMSE is 102 mm, which is the largest, and the correlation coefficient is only 0.34, confirming that the results are unreliable. The number of DS points in (Figure 21(k)) are similar, but the RMSE and spatial correlation coefficients of the BWS algorithm (Figure 21(k)) are lower than those of the FaSHPS algorithm (Figure 20(l)). The above results once again demonstrate the superiority and reliability of the proposed algorithm in terms of DS points and deformation monitoring results. The runtime of different algorithms is shown in (Table 4) and it is notice that FaSHPS is most efficient, followed by BWS-DIE algorithm, which is inconsistent with (Table 1). This is because the BWS-DIE algorithm itself consists of two parts, namely the BWS algorithm and the dynamic interval estimation algorithm (similar to the FaSHPS algorithm), which has slightly more steps than the FaSHPS algorithm. In the Monte Carlo experiment, the sample size is relatively small (15 × 15), the difference in computational efficiency is not obvious, and even the BWS-DIE algorithm is slightly ahead of the FaSHPS algorithm. However, in the real experiment, the size of the SAR image (1000 × 600) is much larger than the simulated sample, so the BWS-DIE algorithm is much slower.

Figure 21. (a-h) is deformation time series comparison results; (i-l) is statistical comparison results between the different algorithms.

Table 4. Runtime of different algorithms.

Algorithm	BWS-DIE	KS	BWS	FaSHPS
Runtime (s)	145	532	383	27

Comparison with other method

Recently (Bao et al. 2021), proposed a fast and accurate distributed scatterer extraction algorithm which combines KS test and interval estimation algorithms (referred as KS-FaSHPS). The experiment results in (Bao et al. 2021) also show that KS- FaSHPS is able to fast identify the largest number of reliable DS compared with KS and FaSHPS. Here we try to compare the performance of the two methods, i.e. KS- FaSHPS and BWS-DIE. Therefore twelve scenes of ALOS PALSAR-1 data over Los Alamitos Army Airfield provided by (Jiang, Ding, and Li 2018) are employed to select homogeneous pixels using the BWS-DIE and KS- FaSHPS algorithms respectively, and the obtained results are shown in (Figures 22(a,b)). (Figures 22(a,b)) shows that visually the results obtained by the two algorithms are smooth and the outlines are clear, and the results of both algorithms are consistent in most areas (Figure 22(c)).

Figure 22. Comparison of homogeneous pixel selection results between BWS-DIE and KS- FaSHPS. It is noticed that red dots indicate reference pixels and yellow dots indicate results of homogeneous pixel selection; the reference pixels are located in road, bare land and building areas in Figures 22(d-g), (e-h) and (f-i), respectively.

To further compare the performance difference between the two methods, the reference pixels located in the road, bare land and building areas were selected for homogeneous pixel selection using BWS-DIE and KS- FaSHPS algorithms, and the obtained results are shown in (Figures 22(d-i)), respectively. (Figures 22(d,f,g,i)) indicate that the homogeneous pixel selection results of the BWS-DIE and KS- FaSHPS algorithms are basically consistent in road and building areas. However, the Type I error of the KS- FaSHPS algorithm in the bare land area is 14.2%, and the correct detection probability reaches 85.8%, compared with 10.6% and 89.6% of BWS-DIE. In addition, the KS- FaSHPS algorithm takes 112 s, whilst the BWS-DIE algorithm takes 78 s. The above comparison shows that the BWS-DIE algorithm overperforms the KS- FaSHPS algorithm in controlling type I errors and calculation efficiency.

Limitations and future work

The current experimental results show that twelve scenes of SAR images can obtain a large number of homogeneous pixels with reliable quality using the proposed method. It should also be noted that the influence of the number of SAR images on the results of homogeneous pixel selection still needs further research (Song et al. 2017). believes that methods such as KS test algorithm only consider temporal samples, which may be unreliable in the case of small data sets, as the samples of each pixel are too few for statistical goodness of fit test. So, they proposed an adaptive multilook approach which uses adaptive joint data vector comprising of temporal and spatial information for homogeneous pixels identification. Experimental results show that reliable results can be obtained with this method even if there are only six SAR images (Zhang et al. 2013). proposed a modified adaptive spatial filtering algorithm for accurate estimation of interferogram and coherence without resolution loss even in rural areas. The implemented method identifies the statistically homogenous pixels in a neighborhood based on the goodness-of-fit test. In their experiment, thirteen scenes of ALOS PALSAR image were used to select homogeneous pixels and monitor deformation in the Yellow River Delta region. Also in the case of small samples, the monitoring density reaches 92.68%. To overcome the insufficient use of complex signals in similarity testing (Wang et al. 2018), proposed an adaptive multilooking method based on complex patches (AMCP). The proposed algorithm takes full advantage of interferometric signals on the basis of the interferometric model for the similarity test and combines the amplitude information and the phase information. It is an efficient algorithm that can obtain reliable homogenous pixel selection results even for small data sets (six scenes SAR image in the experiment). Rather than using the amplitude as a proxy for phase stability, a method is developed by (Pepe, Mastro, and Jones 2020) to evaluate statistical homogeneity based on the interferometric phase directly. The developed approach relies on the application of the fundamentals of directional statistics theory.

The proposed method combines the advantages of parametric testing and non-parametric testing algorithms, but still belongs to the goodness-of-fit testing method in the time dimension. Moreover, compared with the above methods (Song et al. 2017; Pepe, Mastro, and Jones 2020), the proposed method only considers the amplitude information of SAR images, and ignores the interferometric phase signals which contain valuable and correlated information. In addition, the spatial correlation between pixels in the study area is not considered.

Next, we believe that a homogeneous pixel selection framework with regional characteristics should be established, rather than a universal framework which is difficult to achieve, and cannot achieve satisfactory results under different condition. For example, for the monitoring of mining subsidence, the correlation between pixels should be considered when selecting homogeneous pixels, as the law of mining subsidence in mining area is obvious, and generally the subsidence of basin center is the largest, and the edge is smaller. Moreover, there is also a specific functional relationship between subsidence and horizontal displacement (Chen et al. 2020, 2021). Therefore, the laws and characteristics should be considered when selecting homogeneous pixels in coal mining areas, and a comprehensive criterion for homogeneous pixel selection should be constructed in combination with the amplitude and phase information of SAR images to help obtain more density and accurate monitoring points.

In terms of the validation of homogeneous pixel selection results, only the number of homogeneous pixels, running time of the algorithm and deformation monitoring results are considered in the experiment. To more comprehensively validate the feasibility and accuracy of the selection results of homogeneous pixels, more approaches should be considered in the next work. For example, using the homogeneous pixel set for optimization of multitemporal interferometric phases, coherence estimation and intensity image denoising.

Conclusions

The selection of homogeneous pixel samples is a key step for DS InSAR technology, and its estimation precision directly affects the degrees of precision and reliability in subsequent parameter calculations. In this paper, a new sequential selection algorithm for homogeneous pixels is proposed. The algorithm combines the BWS test algorithm with dynamic interval estimation theory, and offers the following advantages: (1) compared with existing algorithms, the BWS-DIE algorithm can significantly improve the power and precision of homogeneous pixel selection; and (2) the algorithm is rarely affected by the number of SAR images and is applicable to scenarios with both large and small samples. In addition, a new DS InSAR method is constructed based on the BWS-DIE algorithm, which discovers new deformation patterns in the spatial-temporal evolution of the surface in the experimental area, by which the application scenarios for DS InSAR technology are further broadened.

To verify the reliability of the algorithm proposed in this study, simulation experiments and real case experiments were simultaneously conducted. The experimental results are as follows: (1) Compared to the KS, BWS and FaSHPS algorithms, the BWS-DIE algorithm increases homogeneous pixel selection precision by 64.3%, 69.4% and 25.3%, respectively. (2) Using the BWS-DIE MT InSAR method, the surface spatial-temporal evolution law of the western mining area in Xuzhou between 20 February 2007 and 3 March 2011 is retrieved. Additionally, new surface deformation patterns in the spatial-temporal evolution were found. (3) Compared with the StaMPS method, KS test algorithm, BWS test algorithm and FaSHPS algorithm, the BWS-DIE algorithm improves the ability to monitor maximum subsidence by 12.3%, 19.3%, 20.6% and 17.3%, respectively, and reduces RMSE by 50%, 68.6%, 56.8% and 45.8%, respectively.

In this study, only the statistical tests of the SAR image amplitude in the time dimension were considered. The properties of SAR images in the space dimension were not applied, and the use of SAR image polarization information was not considered. Combining the time-space property information with the polarization information to further improve the precision of homogeneous pixel estimation will be the priority of future work.

Data availability statement

The MATLAB codes for Monte Carlo simulation experiment, Phase optimization and BWS-DIE algorithm are available at https://github.com/Ming-subsidence/time-series-deformation.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by the [Natural Science Foundation of China] under Grant [number 41907240]; [Key Laboratory of Land Satellite Remote Sensing Application, Ministry of Natural Resources of the People’s Republic of China] under Grant [number KLSMNR-202203]; [the China Postdoctoral Science Foundation] under Grant [number 2019M663601];[A Project Funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions]; [Key Research and Development Program of Xuzhou] under Grant [number KC20180]; [the Shaanxi Province Science and Technology Innovation team] under Grant [number 2021TD-51]; [the Fundamental Research Funds for the Central Universities, CHD] under Grant [numbers 300102260301 and 300102261108]; [the European Space Agency through the ESA-MOST DRAGON-5 project] under Grant [number 59339]; and [International Cooperation and Exchanges National Natural Science Foundation of China] under Grant [number 41920104010].