Monitoring land subsidence through improved CAESAR algorithm in time-series InSAR processing

Abstract

Time-series SAR interferometry (InSAR) combining permanent scatterer and distributed scatterer (DS), has been strongly developed in subsidence monitoring. It is known that the Component extrAction and sElection SAR (CAESAR) is a recently presented approach of selecting and filtering scattering mechanisms for DS. This article proposes an improved CAESAR algorithm in InSAR processing. Phase optimization is performed by eigenvalue decomposition and principal component analysis of the coherence matrix that constructed based on the identified homogeneous pixels. In addition, only the interferometric phases with low noise are used to calculate the goodness-of-fit value. The improved method has been tested for subsidence monitoring over a nonurban area located in Xiongxian, China using 25 Sentinel-1A images. The results show that the improved method can provide the high spatial density of deformation measurements with accuracy ensured. Noticeable land subsidence is revealed widely within the north of Xiongxian county, particularly in Daying, Mijiawun and Beishakou towns.

Keywords:

• Eigenvalue decomposition
• principal component analysis
• distributed scatterer
• permanent scatterer

1. Introduction

Land subsidence is a geological phenomenon in which the surface elevation is decreased under the effect of natural factors or human activities (Chen et al. 2019; Alam et al. 2022). It has the characteristics of slow formation and wide range of influence (Zhang et al. 2016). Land subsidence often causes damage to transportation, industrial and agricultural production, resources and the environment, which restricts the sustainable development of China’s society and economy (Liu 2018; Xu et al. 2021). The monitoring results of land subsidence are the basis for implementing prevention and control projects. Therefore, it is particularly important to carry out the research of time-series subsidence monitoring.

The traditional methods of subsidence monitoring mainly include leveling, total station and global navigation satellite systems (GNSS) measurement (Lu et al. 2019). Nevertheless, these techniques are based on point-by-point measurements, the acquired measurement points (MPs) of which are sparse and time-consuming to collect (Liu et al. 2020). As one of the earth observation technologies, the spaceborne interferometric synthetic aperture radar (InSAR) has the ability to acquire all-time data and provide large-scale surface deformation at a low cost (Robertson et al. 2020).

Permanent scatterer (PS) InSAR with single master image has a wide application in deformation monitoring (Ferretti et al. 2000; Wang and Zhu 2016; Costantini et al. 2016; Xu et al. 2021; Yazici and Gormus 2021). To realize a different methodology in paring acquisitions, Berardino et al. (2002) presented small baseline subset (SBAS) InSAR with multiple master images (Hooper 2008; Pepe et al. 2015; Liu et al. 2018). Compared with PS-InSAR, SBAS-InSAR can increase the temporal sampling rate by using all the acquisitions included in the different small baseline subsets and provide spatially dense deformation (Berardino et al. 2002; Du et al. 2018; Yu et al. 2021). Furthermore, Zhang et al. (2016) combined the advantages of PS and SBAS InSAR, and put forward Multiple-master Coherent Target Small-Baseline InSAR (MCTSB-InSAR), which was successfully applied to monitoring land subsidence over a large area (Zhang et al. 2019, 2021). MCTSB-InSAR uses multiple thresholds such as amplitude dispersion index and mean coherence coefficient to extract high coherent targets and local Delaunay triangulation to increase the redundancy of network connection (Zhang et al. 2016). However, it is difficult for these methods to obtain enough MPs in nonurban areas for detailed deformation information.

Previous studies have shown that DS generally correspond to natural targets, and nonurban areas are characterized by DS instead of PS (Ferretti et al. 2011). These targets, where many neighboring pixels share similar reflectivity values (Li 2014), are distributed scatterer (DS). Ferretti et al. (2011) defined DS for pixels larger than the thresholds of statistically homogeneous pixel (SHP) number and goodness-of-fit value, and used their optimized values to replace the phase values of the original SAR image in the SqueeSAR algorithm (Wang et al. 2016). They utilized the Kolmogorov–Smirnov (KS) test to identify SHP, and applied the phase triangulation method (PT) based on maximum likelihood estimation to optimize the phase. Deformation estimation of the SqueeSAR algorithm is the result of processing the selected DS jointly with the PS in the traditional PS InSAR module.

The key to DS research lies in DS selection and phase optimization. SHP identification is the basis for effective description of DS statistical distribution characteristics (Liu et al. 2019). (Parizzi and Brcic, 2011) compared nonparametric KS test, Anderson–Darling (AD) test and parametric Generalized Likelihood Ratio Test (GLRT). The results have showed that GLRT performs best among the three tests when the amplitude of SAR image satisfies the assumption of Rayleigh distribution. Nevertheless, these tests need the large amount of calculation. To overcome the difficulty of high computational cost, Jiang et al. (2015) put forward the Fast Statistically Homogeneous Pixel Selection (FaSHPS) algorithm that transformed the hypothesis testing problem from significant difference judgment to confidence interval estimation. This algorithm has been proved to be an effective tool for SHP identification (Jiang et al. 2016; Sun et al. 2018; Li et al. 2021). In addition, Cao et al. (2015) improved the PT algorithm and proposed equal-weighted and coherence-weighted PT algorithms for DS phase optimization. The two proposed algorithms increase the signal-to-noise ratio of phase estimation, but greatly reduce computational efficiency. To solve this problem, Fornaro et al. (2015) presented a new algorithm named Component extrAction and sElection SAR (CAESAR), which filtered the dominant scattering component of SAR signal (Cao et al. 2016; Liu et al. 2019).

It is known that the CAESAR algorithm uses covariance matrix eigen-decomposition-based algorithm to extract the power distribution of multi-baseline SAR data (Ansari et al. 2016; Verde et al. 2018). However, this algorithm has several shortcomings: (1) The covariance matrix is estimated based on the mean value of the pixels in the boxcar window. If the pixels in the window do not belong to the same feature, the estimation of covariance matrix is biased. (2) Working on covariance matrix is easily affected by the unbalances of backscattered power between SAR images. (3) In the process of DS selection, the goodness-of-fit value is calculated by the phases of all possible $N(N-1)/2$ interferometric pairs ($N$ is the number of SAR images), ignoring the adverse effect of low-quality interferometric phases.

In view of these shortcomings, this article makes the following improvements. Aiming at shortcomings (1) and (2), the coherence matrix rather than covariance matrix is used to reduce the effect of unbalanced backscattered power. The coherence matrix is constructed by homogeneous pixels identified by the FaSHPS algorithm. According to the shortcoming (3), we select the interferometric pairs with low noise and clear fringes to estimate the goodness-of-fit value. The focus of our research is to obtain as many reliable MPs as possible without loss of measurement accuracy in nonurban areas. We present a new method of coherence matrix eigen-decomposition, which optimizes the phase of multi-baseline SAR data on the basis of principal component analysis (PCA). The effectiveness of the improved method is proved by experiments on the Sentinel-1A datasets, which consists of 25 SAR images acquired by satellites over Xiongxian (China) from 22 January 2019 to 12 December 2019.

The rest of this article is organized as follows. Section 2 introduces the methods adopted in this work. Section 3 describes the study area and datasets. Section 4 shows the results of the experiments. Section 5 discusses the obtained results. Finally, Section 6 is the main conclusions of this study.

2. Study area and datasets

Xiongxian, belonging to Baoding City, Hebei Province, China, is situated in the northeast of Xiong’an New Area and the core of Beijing–Tianjin–Baoding region (Zhang et al. 2018). There are nine towns in Xiongxian county with a total area of about 524 km² (Liu and Xue 2019). Xiongxian county is the hometown of hot springs in China, rich in geothermal resources. The precipitation is about 600 mm and the elevation is within the range of 8 to14 m. The transportation is convenient with Tianjin–Baoding railway running through the whole county. The boundary of the study area is described by the blue rectangle in Figure 1, with geographic coordinates of 115°56′E-116°23′E, 38°52′N-39°12′N.

Figure 1. Location of Xiongxian county (outlined in black). The blue rectangle represents the outline of SAR images.

A total number of 25 Sentinel-1A images are acquired between 22 January 2019 and 12 December 2019 in ascending orbit with the polarization of transmit-vertical and receive-vertical (VV). The size of each pixel of the original Single-Look-Complex (SLC) images is 2.33 m in range and 13.96 m in azimuth, and the incident angle of radar are about 39.48°. The link to download Sentinel-1A images is provided (https://search.asf.alaska.edu/; https://dataspace.copernicus.eu/). In addition, Advanced Land Observing Satellite Digital Elevation Model (ALOS DEM) with a resolution of 30 m are applied for the generation of differential interferogram through the removal of the topographic phase contribution. Table 1 shows the basic parameters of SAR images.

Table 1. Details of the processed SAR data.

Items	Description
Satellite	Sentinel-1A
Imaging mode	TOPS
Orbit direction	Ascending
Polarization	VV
Wavelength (cm)	5.55
Azimuth resolution (m)	13.96
Range resolution (m)	2.33
Incidence angle (deg)	39.48
Number of images	25
Time spans	22 January 2019−12 December 2019

3. Methodology

3.1. CAESAR algorithm overview

CAESAR algorithm is applied to remove noise and optimize phase. The statistical characterization of DS can be described using the covariance matrix (Ferretti et al. 2011). The definition of covariance matrix $\mathbf{C}$ is (1) $$\mathbf{C}=E\left[\mathit{X}{\mathit{X}}^H\right]\approx \frac{1}{N_q}\sum _{X\in \mathit{\text{win}}}\mathit{X}{\mathit{X}}^H$$(1) where $E[\cdot ]$ is the expectation. H represents the conjugate transpose. $\mathit{X}={\left[x_1,x_2,\dots ,x_N\right]}^T$ is the original complex observation of $N$ SAR images. $\mathbf{C}$ is estimated from complex data by performing spatial averaging on $N_q$ pixels that belongs to the boxcar window $\mathit{\text{win}}.$

$\mathbf{C}$ is an $N\times N$ Hermitian matrix whose elements include interferometric phase. The eigenvalue decomposition (EVD) of $\mathbf{C}$ is a linear combination of outer products of the orthonormal eigenvectors that weighted by the corresponding eigenvalues in the following form (Fornaro et al. 2015; Verde et al. 2018). (2) $$\mathbf{C}=\sum _{i=1}^N{\lambda }_i{\mathit{u}}_i{\mathit{u}}_i^H$$(2) where the eigenvalue ${\lambda }_i$ are arranged in decreasing order. ${\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_N.$ $u_i$ is theeigenvector relative to the eigenvalue ${\lambda }_i.$ In line with PCA, the first principal component is the one associated with the first eigenvalue-eigenvector pair $({\lambda }_i,{\mathit{u}}_i).$ $\mathbf{C}$ is approximately a single scattering mechanism ${\mathbf{C}}_{\mathbf{1}}$ (where ${\mathbf{C}}_{\mathbf{1}}={\lambda }_1{\mathit{u}}_1{\mathit{u}}_1^H$). $\mathit{\theta }\text{=}{\left[{\theta }_1,{\theta }_2,\dots ,{\theta }_N\right]}^T$ is the optimal phase extracted from ${\mathit{u}}_1$ (Ansari et al. 2016).

The goodness-of-fit measure $\kappa ,$ has been used as a quality indicator of the CAESAR filtered phase. Indeed, $\kappa$ is generally considered as the temporal coherence of DS (Ferretti et al. 2011; Cao et al. 2015; Jiang and Guarnieri 2020). The value of $\kappa$ is between 0 and 1. It is an essential part of DS selection. (3) $$\kappa =\frac{2}{N(N-1)}\mathrm{\text{Re}}\left(\sum _{m=1}^N\sum _{n=m+1}^N{e}^{j{\varphi }_{m,n}}{e}^{-j({\theta }_m-{\theta }_n)}\right)$$(3) where $\mathrm{\text{Re}}(\cdot )$ denotes the real part of complex data. $j$ is the unit of imaginary part. ${\varphi }_{m,n}$ represents the interferometric phase between the original $m$th and $n$th SAR images. Most notably, $\kappa$ is calculated by $N(N-1)/2$ interferometric pairs before and after optimization. The larger the value of $\kappa$ is, the higher the coherence and the signal-to-noise ratio are.

3.2. Improved CAESAR in time-series InSAR processing

3.2.1. Improved CAESAR algorithm

In view of the shortcomings mentioned in the introduction, one modification aims to substitute the boxcar window with SHPs. FaSHPS algorithm is a widely used method to quickly identify SHP (Jiang et al. 2015).

Assuming that the amplitude $A(p)$ for pixel $p$ has a Rayleigh distribution with mean $u(p)$ and the variance $\mathit{\text{Var}}(p),$ $\sqrt{\mathit{\text{Var}}(p)}=(0.52\cdot u(p))/\sqrt{L}$ for $L-\mathit{\text{looks}}$ image. Based on the central limit theorem, the average amplitude $\overline{A}(p)$ of SAR images with size $N$ is supposed to follow a normal distribution with mean $u(p)$ and the variance $\mathit{\text{Var}}(p)/N.$ On the basis of the definition of standard normal distribution, the confidence interval of $\overline{A}(p)$ is (4) $$P\left\{u\left(p\right)-Z_{1-\frac{\alpha }{2}}\cdot \frac{0.52\cdot u\left(p\right)}{\sqrt{N\cdot L}}<\overline{A}(p)<u\left(p\right)+Z_{1-\frac{\alpha }{2}}\cdot \frac{0.52\cdot u\left(p\right)}{\sqrt{N\cdot L}}\right\}=1-\alpha$$(4) where $P\{\cdot \}$ is the probability, and $Z_{1-\frac{\alpha }{2}}$ is the $1-\alpha /2$ percentile of the standard normal probability density function (Sun et al. 2018; Li et al. 2021). It can be found that the pixels falling into the interval of above formula are homogeneous with the pixel $p$ when a significance level $\alpha$ is given. Moreover, pixels with SHP number greater than a certain value are regarded as the DS candidates.

Using the identified SHP, we improve formula (1). The covariance matrix ${\mathbf{C}}^{\mathbf{\prime }}$ can be written as (5) $${\mathbf{C}}^{\mathbf{\prime }}=E\left[\mathit{X}{\mathit{X}}^H\right]\approx \frac{1}{N_p}\sum _{X\in \Omega }\mathit{X}{\mathit{X}}^H$$(5) where $\Omega$ stands for the set of homogeneous pixels containing $N_p$ adjacent pixels with similar backscattering properties. As previously mentioned, researching on coherence matrix instead of covariance matrix can compensate for unbalanced backscattered power. The coherence matrix $\mathbf{T}$ can be given by (6) $$\mathbf{T}=E\left[\mathit{Y}{\mathit{Y}}^H\right]\approx \frac{1}{N_P}\sum _{Y\in \Omega }\mathit{Y}{\mathit{Y}}^H$$(6) where $\mathit{Y}={\left[y_1,y_2,\dots ,y_N\right]}^T$ is the normalized complex data and $E\left[{\left|y_i\right|}^2\right]=1$ (Cao et al. 2015; Wang et al. 2016). It can be known that the connection between $x_i$ and $y_i$ is $y_i=x_i/\sqrt{E\left[{\left|x_i\right|}^2\right]}.$

$\mathbf{T}$ is also an $N\times N$ Hermitian matrix with elements containing interferometric phase. The EVD of $\mathbf{T}$ can be obtained as (7) $$\mathbf{T}=\sum _{i=1}^N{v}_i{\mathit{\epsilon }}_i{\mathit{\epsilon }}_i^H$$(7) where the eigenvalue $v_i$ are arranged in decreasing order. $v_1\ge v_2\ge \cdots \ge v_N.$ ${\mathit{\epsilon }}_i$ is the eigenvector relative to the eigenvalue $v_i.$ As a result, the phase of ${\mathit{\epsilon }}_1$ is the optimal phase corresponding to the dominant scattering mechanism.

It can be seen from formula (3) that $\kappa$ is calculated from the phases of all possible combined interferometric pairs. To avoid the adverse effects of low-quality interferograms, not all the interferometric phases are involved in the following process of deformation estimation. Therefore, it is inaccurate to calculate the value of goodness-of-fit and evaluate the quality of phase optimization by Eq. (3)(3) $$\kappa =\frac{2}{N(N-1)}\mathrm{\text{Re}}\left(\sum _{m=1}^N\sum _{n=m+1}^N{e}^{j{\varphi }_{m,n}}{e}^{-j({\theta }_m-{\theta }_n)}\right)$$(3) . One improvement is to select those interferometric pairs with low noise and clear fringes to measure the goodness-of-fit.

Assuming the number of selected interferograms is $M$ $\left(\text{where}M<\frac{N(N-1)}{2}\right)$, the goodness-of-fit value ${\kappa }_{\mathit{\text{imp}}}$ can be expressed as (8) $${\kappa }_{\mathit{\text{imp}}}=\frac{1}{M}\mathrm{\text{Re}}\left(\sum _{s\in \Delta }\sum _{t\in \Delta }{e}^{j{\varphi }_{s,t}}{e}^{-j({\theta }_s-{\theta }_t)}\right)$$(8) where ${\varphi }_{s,t}$ indicates the interferometric phase between the original $s$th and $t$th SAR images. ${\theta }_s$ and ${\theta }_t$ are the optimal phases of the $s$th and $t$th SAR images, respectively. Most notably, $t=s+1.$ $\Delta$ represents the set of SAR images corresponding to $M$ interferometric pairs that we select.

DS selection is closely related to the goodness-of-fit value. ${\kappa }_{\mathit{\text{imp}}}$ is estimated from $M$ interferograms filtered by the improved CAESAR module, which theoretically avoids the effect of low-quality interferometric phase. Pixels with ${\kappa }_{\mathit{\text{imp}}}$ higher than a certain threshold are extracted from DS candidates, and the corresponding phase values of the original SAR images are substituted with their optimized values.

3.2.2. Time-series InSAR processing

The improved CAESAR algorithm is performed in the framework of MCTSB-InSAR to extend the analysis of coherent targets interferometry in nonurban areas. This method is defined as the improved CAESAR-InSAR in the following. Figure 2 shows the processing flowchart of the improved method for deformation estimation.

Figure 2. Flow chart of the improved method.

Time-series SAR images with size $N$ are used for fine registration. Interferometric pairs with free combination are generated according to the pre-set thresholds of temporal and perpendicular baseline, and then, $M$ interferograms with lower noise and clearer fringes are selected. After removing the flat-Earth and topographic phase through satellite orbit parameters and external DEM, $M$ differential interferograms and coherence images are generated.

Some coherent targets (PS and a part of the DS) can be extracted based on the threshold of amplitude dispersion index and mean coherence (Zhang et al. 2016). Furthermore, we increase the density of DS by using improved CAESAR algorithm, which is on the basis of coherence matrix eigen-decomposition. The selected DS and PS pixels are connected by the local Delaunay triangulation (Liu et al. 2008; Zhang et al. 2011), which has advantages in the redundancy of network connection.

Surface deformation can be divided into linear and nonlinear displacement (Mora et al. 2003; Zhang et al. 2016). We establish the observation equation among linear deformation rate, DEM error and residual phase of those two selected adjacent pixels. The least square algorithm is used to estimate the linear deformation rate and DEM error. Unlike SBAS, MCTSB-InSAR completes the deformation parameter calculation and phase unwrapping (PhU) at the same time by investigating the quadratic phase difference of two adjacent coherent targets in space (Zhang et al. 2018, 2019). GACOS atmospheric correction model is applied for the generation of atmospheric phase (Ding et al. 2008; Chen et al. 2017). The nonlinear deformation phase, noise phase and atmospheric phase are separated from the residual phase by temporal-spatial filtering algorithm (Li et al. 2021). In particular, the PhU consistency index $\gamma$ is evaluated to check the reliability of unwrapping (Verde et al. 2018). (9) $$\gamma \text{=}\frac{1}{M}\left|\sum _{i=1}^M{e}^{j{R}_i}\right|$$(9) where $R_i$ is the residual errors of the $i$th interferogram. The value of $\gamma$ varies from 0 to 1. The closer its value is to 1, the lower the residual errors produced during the PhU process. The coherent pixels with consistency index less than the certain threshold (set to 0.8 in this work) are masked.

The linear and nonlinear deformations are combined to form the line-of-sight (LOS) displacement, which can be converted into vertical subsidence by using the cosine value of the incident angle. The subsidence map in the SAR coordinate system is converted to that in the geographic coordinate system by geocoding.

4. Results

4.1. SHP identification and phase optimization results

For 25 Sentinel-1A images ($N=25$) in TOPS mode, two consecutive burst SLCs need to be concatenated into one larger burst SLC. It is necessary to deramp the SLCs for the azimuth spectrum ramp. The SLCs are co-registered through the enhanced spectral diversity (ESD) algorithm with an accuracy of a few thousands of a pixel in the azimuth direction. We set the number of looks in the azimuth and range directions to 2 and 8, respectively. After implementing multi-look of SLCs, 1135 × 1115 pixels are cropped for the following experiments.

Figure 3(a,b) is the average amplitude and SHP number of the study area, respectively. Homogeneous pixels are identified by FaSHPS algorithm when $\alpha =5\%.$ The value of each pixel in Figure 3(b) represents the number of SHP selected in a 15 × 15 window. SHP number is 0, which indicates that there are no homogeneous pixels in the window perhaps because of isolated PS or noise points. SHP number is 224, which demonstrates that all the pixels in the window are homogeneous. It can be seen from Figure 3 that SHP number is relevant to land cover type. The number of SHP in buildings and structures is relatively small, but different from that in farmland.

Figure 3. The map of average amplitude and SHP number in the study area, respectively.

It is well known that 25 SAR images can generate up to 300 interferometric pairs. Interferograms of multiple-master images are generated under the limit of temporal baseline less than 120 days and perpendicular baseline smaller than 150 meters, respectively. We select 100 interferograms ($M=100$) with less noise and clearer fringes afterwards. The distribution of interferometric pairs is shown in Figure 4.

Figure 4. Distribution of interferometric pairs. Black triangles stand for SAR acquisitions.

To compare the differences of CAESAR algorithm before and after improvement, we carry out two types of phase optimization according to the principle mentioned in Subsection 2.1 and Subsection 2.2.1. The goodness-of-fit values ($\kappa$ and ${\kappa }_{\mathit{\text{imp}}}$) are calculated to evaluate the quality of optimized phase, as shown in Figure 5. It can be found that the goodness-of-fit value of vegetation is smaller than that of buildings and structures. 88% of the pixels in Figure 5(b) show greater goodness-of-fit value than that in Figure 5(a). Most PS (such as pixels in buildings and structures) with goodness-of-fit value greater than 0.7 have less than 20 SHPs, so we regard pixels with SHP number more than 20 as DS candidates. The threshold of goodness-of-fit is set to 0.7 in this work. The proportion of pixels with ${\kappa }_{\mathit{\text{imp}}}\ge 0.7$ is higher than that of $\kappa \ge 0.7,$ which makes it possible to increase the number of DS.

Figure 5. The goodness-of-fit value of the study area obtained by (a) CAESAR algorithm and (b) improved CAESAR algorithm.

4.2. Deformation rate analysis

As stated in Subsection 2.2.2, the deformation estimation is performed through the improved CAESAR algorithm in the MCTSB-InSAR framework, which is defined as the improved CAESAR-InSAR. Similarly, the deformation estimation through the CAESAR module in MCTSB-InSAR processing chain, is named as CAESAR-InSAR. To testify the effectiveness of the improved method, we have compared the monitoring results of three methods (MCTSB-InSAR, CAESAR-InSAR and improved CAESAR-InSAR).

To evaluate improvements in the time-series processing, an important point is to assess the effectiveness of unwrapping. The PhU consistency index $\gamma$ of these three methods are evaluated, as shown in Figure 6. There are a larger number of occurrences when $\gamma$ is close to 1, especially in the improved CAESAR-InSAR processing. It indicates that the PhU improvement is more obvious by the improved method.

Figure 6. The percentage of the PhU consistency indexes.

The vertical deformation rate map in Figure 7 is obtained by three different methods. The location of land subsidence in Figure 7(a–c) are roughly the same. There is no land subsidence in the southwest of Xiongxian county. Noticeable land subsidence ($\le$ −10 mm year⁻¹) is situated in the northern part of Xiongxian county, particularly in Daying, Mijiawu and Beishakou towns.

Figure 7. The deformation rates derived from Sentinel-1A images with different methods. (a) MCTSB-InSAR result, (b) CAESAR-InSAR result and (c) improved CAESAR-InSAR result.

Twenty-nine thousand, seven hundred two points (PS and a part of the DS) are shown in Figure 7(a) with the maximum deformation rate of −116 mm year⁻¹. Most of the PS are concentrated in the area where buildings are gathered in Xiongxian county, and no MP is in the water area. There are 113,638 points in Figure 7(b) with the maximum rate value of −115 mm year⁻¹. Compared with Figure 7(a), the extra points in Figure 7(b) are the number of MP extracted by CAESAR algorithm. Additionally, 154,607 points are exhibited in Figure 7(c) with the maximum rate value of −116 mm year⁻¹. The extra points in Figure 7(c) are a proof of the increase in the number of MP selected by the improved CAESAR algorithm.

It is obvious that the number of MP in Figure 7(c) is more than that in Figure 7(a,b). Figure 7(c) shows more detail deformation than Figure 7(a,b) in the area of noticeable subsidence, for example Daying town. The density of MP in Figure 7(c) is 1.4 times of that in Figure 7(b), increasing greatly in Zhanggang town. In summary, the improved CAESAR-InSAR is effective in high spatial density of deformation measurements.

5. Discussion

5.1. Verification of monitoring results

During the monitoring period, we are unable to obtain the measured leveling data. To verity the monitoring results derived from different methods (MCTSB-InSAR, CAESAR-InSAR and improved CAESAR-InSAR), the correlation coefficient $r$ is calculated from the deformation rates of all the completely coincident points. The absolute value of $r$ (ranging from 0 to 1) measures the degree of close correlation between two different methods (Lu et al. 2019; Liu et al. 2020). Assuming that $\upsilon$ and $\mathit{\mu }$ represent the deformation rates estimated by two methods, their correlation coefficient $r$ can be expressed as: (9) $$r=\frac{\mathit{\text{Cov}}\left(\upsilon ,\mathit{\mu }\right)}{\sqrt{D\left(\upsilon \right)}\sqrt{D\left(\mathit{\mu }\right)}}$$(9) where $\mathit{\text{Cov}}(\cdot )$ and $D(\cdot )$ stand for the covariance and variance, respectively (He et al. 2021).

The value of $r$ is 0.988 in Figure 8(a), which demonstrates that the deformation rates in Figure 7(a) are in consistency with those in Figure 7(b). The correlation coefficient in Figure 8(b) is 0.986, which indicates that the deformation rates of Figure 7(c) is relatively close to that of Figure 7(b). That is, the monitoring results of improved CAESAR-InSAR method are as accurate as those of MCTSB-InSAR method and CAESAR-InSAR method. Overall, the improved method is an effective way to increase the density of deformation measurements with accuracy ensured.

Figure 8. The scatter diagram and regress line of completely coincident points’ deformation rates. (a) Correlation between Figure 7(a,b) and (b) correlation between Figure 7(b,c).

5.2. Damage caused by land subsidence

To obtain the detailed damage information caused by land subsidence, area statistics are performed according to the results of the improved CAESAR-InSAR. Referring to the existing literature (Hu et al. 2009; Wang et al. 2014), the grade of dangerous assessment index for land subsidence is different in different areas and actual situation. In this study, we consider the area with deformation rate less than −10 mm year⁻¹ as dangerous area, and then, classified these dangerous areas into four grades. A total of four groups are deformation rate values of $[-116,-100],$ $(-100,-50],$ $(-50,-30],$ $(-30,-10].$ It can be seen from Figure 7 that some public facilities have been suffered from land subsidence, including Tianjin–Baoding railway. The area of Xiongxian and the length along the railway line are statistically counted on the basis of the dangerous grading (see Table 2).

Table 2. Area statistics of Xiongxian and length statistics along the railway line.

Deformation rate (mm year^-1)	Grading name	Area (km^2)	Length (km)
(–30, −10]	Low dangerous	127.44	4.34
(–50, −30]	Medium dangerous	54.45	2.06
(–100, −50]	High dangerous	49.43	6.41
$\le$−100	Very high dangerous	1.33	0

The total subsidence area (less than −10 mm year⁻¹) affected by land subsidence is about 232.65 km². The railway runs through Zhugezhuang, Daying, Mijiawu and Shuangtang towns, with the maximum rate value of −94 mm year⁻¹. Except for Shuangtang town in the above four towns, the railway line suffered from surface deformation is 12.81 km long. Most of the severely affected railway line is located in Daying town. The damage will bring hidden dangers to people’s life and property safety. In light of this, related departments should give particular attention to this matter.

Based on the analysis of the existing data, land subsidence in Xiongxian county is mainly caused by the exploitation of geothermal and groundwater resources (Pang et al. 2017; Xia and Zhang 2017; Zhang et al. 2018). The geothermal distribution region of Xiongxian county that belongs to the southwest of Niutuozhen geothermal field, accounts for 60% of total area. The geothermal resources are shallow buried, with good water quality and great supply potential. Due to the increase of the exploitation of geothermal resources, the water level of geothermal field has decreased greatly. Most areas of Daying town have good conditions for geothermal development and a large number of geothermal wells, which has caused significant land subsidence in Daying town. As early as 2003, Xiongxian County began using geothermal energy through geothermal wells for industrial, living and tourism purposes. As of October 2018, Xiongxian County has collected a total of 52 geothermal wells. Land subsidence in Xiongxian County is closely related to geothermal exploitation (Li et al. 2022).

In addition, Xiongxian county, with more than 3000 plastic packaging enterprises, is the largest plastic packaging production base in northern China. Large-scale production of plastic packaging is likely to increase the groundwater exploitation and enhance the subsidence. Furthermore, water used in agriculture and aquaculture also increases the possibility of land subsidence. All the information provides the basis for taking corresponding prevention actions and is helpful to make land reclamation and ecological restoration plans.

5.3. Limitation and future work

In this work, we set the goodness-of-fit threshold in both CAESAR algorithm and improved CAESAR algorithm to 0.7. Most notably, the number of DS and the value of deformation are different with different goodness-of-fit thresholds. To be specific, the larger the value of goodness-of-fit is, the smaller the number of DS is, which cannot provide more detailed deformation. The smaller the value of goodness-of-fit is, the larger the number of DS is, which has a negative impact on the measurement accuracy. In view of this, the determination of goodness-of-fit threshold will be further analysed to help reasonably select DS pixels. Furthermore, we plan to use long time-series SAR images to enhance the real-time monitoring of land subsidence in the future.

6. Conclusions

A new method that based on improved CAESAR algorithm in time-series InSAR processing is presented to extend the analysis of coherent targets interferometry in nonurban areas. In the proposed strategy, the coherence matrix eigen-decomposition is used to optimize the phase of multi-baseline SAR data on the basis of PCA. The most contribution in our work is DS selection, which is a unique feature of the improved CAESAR-InSAR method. Specifically, only the phases of interferometric pairs with low noise and clear fringes are used to estimate the goodness-of-fit value, which is a quality indicator of phase optimization and an important parameter for DS selection.

The improved method is tested for land subsidence monitoring in nonurban areas of Xiongxian, China using 25 Sentinel-1A images. The correlation coefficient is calculated from the deformation rate of the completely coincident points obtained by three different methods (MCTSB-InSAR, CAESAR-InSAR and improved CAESAR-InSAR) to verify the monitoring results. Additionally, the area statistics of Xiongxian and the length along the Tianjin–Baoding railway line are statistically counted according to the dangerous grading, and the causes of land subsidence in Xiongxian are preliminary analysed.

Our experiments show that the number of MPs extracted by MCTSB-InSAR is much less than that obtained by CAESAR-InSAR and improved CAESAR-InSAR. The improved strategy can provide more deformation measurement details than other methods with measurement accuracy ensured, which proves that the improved CAESAR-InSAR method has great application potential in land subsidence.

There is noticeable land subsidence in the north of Xiongxian county, especially in Daying, Mijiawu and Beishakou towns. The main influencing factors of land subsidence are the exploitation of geothermal and groundwater resources. This research offers information to relevant departments for risk management purpose. Further analysis will be conducted on Xiongxian county to offer long-term displacements to relevant management departments.

Acknowledgement

We want to say our warmest thanks to reviewers for their constructive suggestions for this article.

Disclosure statement

No potential conflict of interest was reported by the authors.