Effectiveness evaluation of DS-InSAR method fused PS points in surface deformation monitoring: a case study of Hongta District, Yuxi City, China

Abstract

InSAR technology provides a powerful tool for detecting large-scale surface deformation. In particular, the newly developed DS-InSAR method fused PS points has obvious advantages in monitoring bare land and vegetation covered areas, but at present, there is a lack of effectiveness evaluation of this method. Therefore, this paper takes Hongta District of Yuxi City as an example. The 29 scenes sentinel-1A data from January 2019 to December 2019 were processed and analyzed using the DS-InSAR method fused PS points. The research results show that the position and deformation trend of the inversion results of the three methods are highly consistent, the correlation between DS-InSAR method fused PS points and PS- InSAR and SBAS-InSAR methods is 0.9473 and 0.8583, respectively. But the spatial density of the measurement points (MPs) obtained by the DS-InSAR method of PS points is 11 times and 4 times that of the PS-InSAR method and SBAS-InSAR method respectively. It clearly shows that DS-InSAR method fused PS points has greater advantages than PS-InSAR and SBAS-InSAR methods in deformation mapping displacement, because the MPs density is higher. It is conducive to the detailed analysis of the spatio-temporal characteristics and deformation mechanism of deformation.

Keywords:

• DS-InSAR
• surface subsidence monitoring
• sentinel-1A
• feature decomposition
• phase optimization

1. Introduction

Surface subsidence is an environmental geological phenomenon that causes the gradual or sudden subsidence of the earth’s surface under the action of natural and man-made affiliation (Zheng et al. 2002; Liu et al. 2015). It is a slow and progressive process, once formed, it is difficult to recover (Yin et al. 2005). Surface subsidence is one of the most common geological hazards in many countries. For example, more than 40 states in the United States and nearly 20 provinces in China have experienced land subsidence to varying degrees (Dong et al. 2014). It was severely destructive, bringing serious environmental problems and large economic losses, such as damage to urban infrastructure, including bridges (Selvakumaran et al. 2018), buildings (Yang et al. 2016), railways (Karanam et al. 2021), highways (Zhu et al. 2021), and underground facilities (water supply pipelines, gas, power facilities, etc.) (Jiang et al. 2018; Hu et al. 2014; Garg et al. 2022).

Although the traditional monitoring methods have relatively reliable monitoring results at the monitoring points, their spatial resolution is low, it is difficult to obtain the deformation at non-monitoring points, and the monitoring cost is high, so it is difficult to carry out fine monitoring of wide-area surface deformation, which limits our comprehensive understanding of the deformation field and the spatio-temporal evolution characteristics of deformation in the entire study area (Luo et al. 2014). Synthetic Aperture Radar Interferometry (InSAR) is a new type of earth observation technology that has been developed rapidly in the past 30 years. It is widely used in various kinds of deformation monitoring owing to its high spatial resolution, all-day, all-weather and other characteristics (Zhu et al. 2017; Yuan et al. 2020; Zhang et al. 2022). Compared with conventional geodetic techniques such as global positioning system (GPS) and levelling, InSAR technology can obtain wide-area surface deformation time series with significantly improved spatial resolution (Bai et al. 2016).

In recent years, many advanced InSAR methods and some improved methods have emerged (Cao et al. 2016). The PS-InSAR technique is the first systematic time-series processing solution, whose main idea is to select one scene image from the image collection as the reference primary image, other images as secondary images and the reference primary image for differential interference processing (Ferretti et al. 2002a, 2002b). Select the point target with high signal noise ratio (SNR) and coherence, and achieve the extraction of deformation information by the optimal estimation is used to extract the deformation information and effectively overcome the decoherence and atmospheric effects of the images (Ferretti et al. 2011). With the widespread application of PS-InSAR technology, the drawback that its coherent targets are not effectively mined has become increasingly obvious, especially in non-urban areas where the PS point density is usually less than 10/km², making it difficult to obtain continuous surface deformation information (Samiei-Esfahany et al. 2016). Nevertheless, its technical ideas have been influencing the subsequent development of time-series processing techniques.

The SBAS-InSAR method, which was proposed almost simultaneously with PS-InSAR, fully considers the mining of distributed target information. In order to overcome the effect of decoherence, this technique usually selects only the interferometric pairs with relatively short baselines and performs singular value decomposition of the deconvolution results of a subset of interferometric pairs according to the redundant observations formed by the combination of image interferometric pairs to improve the accuracy of deformation monitoring (Berardino et al. 2002;  Yastika et al. 2019). Although this technique has been improved, it has not been able to escape from the dilemma of two-dimensional deconvolution, and the accuracy of deformation monitoring has not been improved fundamentally (Lauknes et al. 2011).

To overcome the shortcomings of the above InSAR methods. Ferretti et al. (2011) proposed the SqueeSAR technique, which determines whether a pixel is homogeneous by the similarity of each pixel intensity statistical distribution and performs spatial adaptive filtering based on the set of homogeneous pixels, and then uses the Phase-Linking technique to achieve the parity calculation of the winding phase (Sun et al. 2018). Since then, the DS-InSAR technique has been widely used, for example. Pawluszek-Filipiak et al. (2020) combined SBAS-InSAR and D-InSAR to study mining subsidence. Seyfallah Bouraoui (2013) merged PS and small baseline (SB) methods for oilfield, landslide, and seismic subsidence monitoring. Hu et al. (2014) analyzed the temporal and spatial characteristics of land subsidence in the Beijing area using SBAS-InSAR technology, revealing that the monitoring results of InSAR technology are consistent with the levelling results. Luo et al. (2014) adopted Multi-temporal InSAR (MT-InSAR) technology to monitor and analyze land subsidence in the Tianjin area, and verified it with GPS data. The results suggested that the monitoring results of the two monitoring methods were in good agreement. Lubna et al. (2018) employed ISBAS-InSAR technology to monitor the long-term surface movement of peatlands in northeastern Scotland, and the monitoring results were consistent with the trend of surface movement. Zhang et al. (2019) monitored and investigated the surface subsidence in Wuhan with SBAS-InSAR technology, and validated its results through levelling data. With the landslide in Danba County, Sichuan as an example, Zhang et al. (2018) and Dong et al. (2018) used the time-series InSAR technology to identify the hidden hazards of landslide disasters early and demonstrated the reliability of the results through GPS data. Zhao et al. (2019) conducted multi-temporal cataloging and long-term monitoring of potential loess landslides in Heifangtai, Gansu using InSAR technology, combined topographic maps and optical remote sensing images to further explore the instability modes of landslides, and obtained the accuracy of the landslide instability model through field surveys to verify the performance. Given the time-consuming issue of time-series analysis in acquiring images and data preprocessing in the study area, a new method of time-series analysis, LiCSBAS, was adopted to monitor surface deformation (Gao et al. 2020). The results unveiled that this method can considerably improve the efficiency of the solution. Yao et al. (2022) used a cross-comparison approach to study landslides along the Jinsha River in Gongjue County using the SBAS-InSAR and PS-InSAR techniques. Bui et al. (2021) monitoring recent land deformation detected in Hanoi, vietnam using InSAR technology and analyzing its relationship with groundwater level changes. Orhan (2021) investigate the causes of the land subsidence problem in Konya and its surroundings with the help of the SBAS-InSAR technique and auxiliary data, namely optic, Coordination of Information on the Environment (CORINE), and groundwater monitoring station data.

In summary, the DS-InSAR method fused PS point has been successfully applied to various surface deformation monitoring, but at present, only the DS-InSAR method fused PS points is used for surface deformation monitoring and analysis, lacking effective comprehensive evaluation and analysis. Therefore, this paper takes Hongta District of Yuxi City as the experimental object and uses the above the DS-InSAR method fused PS points to monitor its surface deformation, and The results are verified and analyzed with PS-InSAR and SBAS-InSAR results to highlight the advantages of DS-InSAR method fused PS point in surface deformation monitoring. The results of the study provide a scientific basis for the selection of the time-series InSAR technique.

2. Distributed scatterer InSAR method fused with PS

Distributed scatterers (DS) are implemented in the presence of random resolution scatterers with similar amplitudes, and none of them exhibit dominant behaviour. The DS candidate points (DSC) are obtained by comparing the similar features of the interference phase in the time dimension, and then homogenous filtering is performed on them. The phase is optimized in the complex domain to obtain the final DS points. This processing can effectively suppress the influence of noise and improve the signal-to-noise ratio (SNR), contributing to boosting the quality of the interferogram phase and obtaining deformation information with higher precision and enhanced spatial resolution. Compared with other time-series processing methods, this method can effectively improve the spatial density of coherent scatterer points. The key technologies of DS-InSAR involve (1) improving the estimation accuracy of the time series InSAR covariance matrix by performing spatial adaptive filtering on the pixel set of the homogeneous scatterer; (2) recovering the time series SAR image from the covariance matrix through the phase optimization algorithm to achieve the time-dimension phase adjustment.

Anderson-Darling (A-D) hypothesis test is an adaptive multi-view non-parametric test method improved based on the Kolmogorov-Smirnov (K-S) test [13]. Similar to the K-S test, the A-D test statistics are also dependent on the empirical cumulative distribution function (ECDF) of the amplitude between two pixels to judge their similarity. During the processing, a proper estimation window is set with each pixel of the image as the centre, and the consistency of the pixels in the window is judged by the statistical test method. It is assumed that N SAR images are given, one is selected as the primary image, and the remaining (N-1) images are used as auxiliary images. The amplitude empirical cumulative distribution function (ECDF) of any pixel p can be calculated by: (1) $$F_{A_p}\left(A\right)=\frac{1}{N}\sum _{i=1}^N{I}_{A_{p,i}\le A}$$(1) where $I_{A_{p,i}\le A}$ indicates a function. When $I_{A_{p,i}\le A},$ it is equal to 1, otherwise it is 0. For any two pixels $(p,q)$ in the statistical window, the AD test statistics can be expressed as: (2) $$T_{AD}=\frac{N}{2}\sum _A\frac{{\left(F_{A_p}\left(A\right)-F_{A_q}\left(A\right)\right)}^2}{F_{A_{p,q}}\left(A\right)\left(1-F_{A_{p,q}}\left(A\right)\right)}$$(2) where $F_{A_{p,q}}\left(A\right)$ denotes the empirical cumulative distribution function of two pixels p and q combined into one pixel. The A-D test statistic is more complicated than the K-S test statistic, and there is no fixed solution. Pettitt (1976) detailed the approximation of the A-D test and the selection of the threshold according to the significance level. Compared with the K-S method, the A-D test pays more attention to the tail of the distribution. Parizzi and Brcic (2011) demonstrated that the A-D test has better performance than other non-parametric tests for the homogeneous point selection in the synthetic aperture radar stack.

After the homogenous points are obtained by the A-D test, the eigendecomposition method is adopted to estimate the parameters of the covariance matrix of the homogenous points. The reflection signal received from the DS resolution unit is equal to the sum of the coherent information of the reflections from multiple elementary scatterers in the unit. Following the central limit theorem, its pixel values obey the complex circular Gaussian distribution. It is assumed that a set consists of N single-view complex SAR images. This set can be regarded as a complex random variable ${\mathrm{y}=[Z_1,Z_2,.,Z_N]}^T,$ and its probability density function (PDF) can be expressed as: (3) $$f\left(y\right)=\frac{1}{{\pi }^N\left|Q_y\right|}\exp \hspace{0.17em}\left(-y^H{Q}_y^{-1}y\right)$$(3) where $Q_y=E\left\{y{y}^H\right\}$ indicates the $N\times N$ complex covariance matrix, and H represents the Conjugate multiplication.

The complex coherence coefficient $r_{ij}$ between the two images $Z_i$ and $Z_j$ can be expressed by Chen and Meng (2020): (4) $$r_{ij}=\frac{E\left\{Z_i{Z_j}^H\right\}}{\sqrt{E\left\{\right|Z_i{|}^2\left\}E\right\{\left|Z_j{|}^2\right\}}}=\left|r_{ij}\right|\hspace{0.17em}\exp \hspace{0.17em}\left(j{\mathrm{\varphi }}_{ij}\right)$$(4) where ${\mathrm{\varphi }}_{ij}$ denotes the complex coherent phase value between the two scene images, and $\left|r_{ij}\right|$ signifies the absolute value of the coherence coefficient, which is a measure of the correlation between the two scene images. The noise component between the two scene images has relatively high similarity when the coherence is high. Hence, most of the noise components will cancel out during the interference phase. The complex correlation matrix between images is called the coherence matrix, defined as: (5) $$T=\left[\begin{array}{ccc}\begin{array}{cc}\begin{array}{c}1\\ r_{12}^H\end{array}& \begin{array}{c}{r}_{12}\\ 1\end{array}\end{array}& \cdots & \begin{array}{c}{r}_{1N}\\ r_{2N}\end{array}\\ \begin{array}{cc}⋮\mathrm{ }& ⋮\end{array}& \ddots & ⋮\\ \begin{array}{cc}{r}_{1N}^H& r_{2N}^H\end{array}& \cdots & 1\end{array}\right]$$(5)

The coherence matrix (5)(5) $$T=\left[\begin{array}{ccc}\begin{array}{cc}\begin{array}{c}1\\ r_{12}^H\end{array}& \begin{array}{c}{r}_{12}\\ 1\end{array}\end{array}& \cdots & \begin{array}{c}{r}_{1N}\\ r_{2N}\end{array}\\ \begin{array}{cc}⋮\mathrm{ }& ⋮\end{array}& \ddots & ⋮\\ \begin{array}{cc}{r}_{1N}^H& r_{2N}^H\end{array}& \cdots & 1\end{array}\right]$$(5) is a positive definite Hermitian matrix, which can be written as: (6) $$T=V\mathrm{\Lambda }{V}^H$$(6) where $\mathrm{\Lambda }=\mathrm{diag}({\mathrm{\lambda }}_1,{\mathrm{\lambda }}_2,{\mathrm{\lambda }}_3,\dots ,{\mathrm{\lambda }}_N)$ denotes a diagonal matrix composed of non-negative eigenvalues with diagonal elements, and $V=[v_1,v_2,v_3,\dots ,v_n]$ indicates its corresponding feature vector. The eigenvalues are sorted in descending order to obtain ${\lambda }_1\ge {\lambda }_2\ge {\lambda }_3,\dots ,\ge {\lambda }_N.$ T can be interpreted as the sum of $N$ scatterer coherence matrices, where each coherence matrix represents an independent scatterer. (7) $$\mathrm{T}=\sum _{i=1}^N{\mathrm{\lambda }}_i{T}_i=\sum _{i=1}^N{\lambda }_i{v}_i{v_i}^H$$(7)

When the following $(N-m)$ eigenvalues are small, the coherence matrix can be written approximately as: (8) $$\mathrm{T}=\sum _{i=1}^m{\mathrm{\lambda }}_i{T}_i=\sum _{i=1}^m{\lambda }_i{v}_i{v_i}^H$$(8)

Generally, the characteristic eigenvector corresponding to the maximum eigenvalue is its optimal phase value. Due to the uncertainty of the obtained eigenvectors (Zhu et al. 2017), there are multiple eigenvectors for the same eigenvalue. In other words, if the eigenvector corresponding to the eigenvalue $\lambda$ of matrix $A$ is $v,$ then: (9) $$Av=\lambda v$$(9)

For any non-zero complex constant $\mathrm{\alpha },$ there is: (10) $$A\left(\mathrm{\alpha }v\right)=\lambda \left(\mathrm{\alpha }v\right)$$(10)

Since $\mathrm{\alpha }v$ is also the eigenvector of $A,$ the eigenvector multiplied by a non-zero complex constant is an eigenvector. This will result in a constant offset between the actual phase and the estimated phase. Thus, the original interference phase value and estimated phase value of the primary image are set to zero in our study to calculate and eliminate the above offset.

After the final DS point is obtained, the 2D linear regression analysis of the unfolded interference phase is performed with PS to correct the error introduced by the external DEM. Subsequently, the residual atmospheric phase and nonlinear deformation are removed by the difference in time and space characteristics.

3. Research area and experimental data

3.1. Study area

The study area is located in the hinterland of central Yunnan, in the interval of 24°08′30″–24°32′18″ north latitude and 102°17′32″–102°41′37″ east longitude. It is connected to Jiangchuan District in the east and Tonghai County in the southeast. The southwest borders with Eshan Yi Autonomous County, and the north borders with the Jinning District of Kunming City. The transportation in the area is convenient. The 213 National Road and Kun-Yu Railway run through the north and south, forming a north-south transportation hub in Yunnan Province. The average annual total amount of water resources is 430 million cubic metres, of which groundwater accounts for 29%. The specific location of the study area is exhibited in Figure 1.

Figure 1. Study area location and DEM. The left picture shows the location of the study area and the right picture shows the DEM of the study area.

Figure 2. Stratum lithology types in the study area.

The study area belongs to the middle-north section of Yuxi north-south tectonic basin, which is a tectonic denudation low-mountain (hilly) topography, and is in the long-term bending area on the platform base, so the topography and geomorphology vary greatly in east-west direction, and the topography does not change much in north-south direction, but the stratigraphic lithology varies greatly, and the overall area is the intersection of Quaternary lake-phase alluvium, siltation and mountain margin, and there is bedrock outcropping only in local sections. Only the localized rock is revealed. According to relevant geological data, the formation of geodorrocysts in the study area mainly adopts artificial filling, sludge clay, powder clay, gravel soil, fine sandstone, muddy, and gray rock (Figure2), with poor stratosphere stability. The pleats in the study area are relatively developed, exhibit relatively fragile environmental geological conditions, and could induce corresponding geological disasters related to human activities.

3.2. Experimental data and processing

In this study, 29 scenic C-band descending sentienl-1A images in the study area are selected as experimental data, and their imaging time span is from January 2019 to December 2019. The polarization method is VV polarization, the wavelength is 5.63 cm, the revisit cycle is 12 days. The resolution is 5 m × 20 m, and each image consists of three sub-bands and multiple bursts. In order to remove the influence of terrain phase, the elevation data of the Shuttle Radar Topography Mission (SRTM-1) is selected as a reference, and the resolution is 30 m. The specific experimental data parameters are listed in Table 1.

Table 1. Experimental data parameters.

SAR sensor	Sentinel-1A
Orbit direction	Descending
Path	62
Band	C-band
Wavelength	5.6 (cm)
Repeat cycle	12 (days)
No. of images	29
Temporal coverage	Jan 2019-Dec 2019

The distributed scatterer InSAR method data processing flow of the PS point is illustrated in Figure 3. Since the azimuth registration accuracy of sentinel-1A satellite needs to reach 0.001 pixel during interference processing. On this basis, the precision track geometry enhanced spectrum set registration is performed (Wu et al. 2019). Based on the comprehensive consideration of the influence of temporal baseline, space baseline and doppler centroid frequency, the image data of April 22, 2019 is selected as the primary image (see the blue triangle in Figure 4), and other secondary images are registered to the primary image, and then the interferograms is generated by multiple primary images. When selecting PS points, the amplitude deviation threshold method is used to initially candidate PS points, and the amplitude deviation threshold is set to 0.4. The temporal coherence threshold is set to 0.35.

Figure 3. Distributed scatterer InSAR data processing flow of fusion PS. The red box is PS point selection process, the green box is DS point selection process.

Figure 4. Time-space baseline map. The blue triangle represents the primary image, the red rectangle represents the secondary image.

When the homogenous filter is processed, the filter window is set to 49 pixels, the orientation is 13 pixels, and the homogeneous pixel number threshold is minimized to 50 pixels. The coherent matrix and covariance matrix estimations are performed by interference graphics. The Feature decomposition (EVD) method is employed to decompose the coherent matrix, and phase optimization estimation is implemented. The original differential interferograms and the differential interference chart after phase optimization are provided in Figure 5. Figure 5(a and c) show the original differential interferograms for the shortest temporal baseline (12d) and the longest temporal baseline (252d), respectively. Figure 5(b) and (d) shows the differential interferograms after phase optimization for the shortest temporal baseline (12d) and the longest temporal baseline (252d), respectively. Figure 5 suggests that the phase optimization is significantly better compared to the original interferograms, especially when the temporal baseline is long. Besides, the non-urban area is more significant, as demonstrated outside the black rectangular frame in Figure 5, inside the urban area.

Figure 5. Differential interferogram before and after phase optimization. (a) and (c) are interferograms before phase optimization, and (b) and (d) are interferograms after phase optimization.

We choose bust6, burst7 and burst8 for data processing, and when using the DS-InSAR method of integrating PS points for data processing, it uses the joint PS and DS points for unified networking. Additionally, the corresponding parameter estimation algorithm is adopted to solve its phase observation value Wrap and error separation. Finally, the deformation rate and cumulative deformation in the study area are obtained.

4. Results analysis

4.1. Deformation rate analysis

Figure 6 shows the deformation rate obtained by PS-InSAR, SBAS-InSAR and DS-InSAR method integrating PS points. A positive value indicates that the monitoring scatterer is close to the line of sight (LOS) movement; a negative value reflects that the line of sight (LOS) is moved away. The red rectangular frame in the figure is where the Red Tower main city is located. The measurement points obtained by the PS-InSAR method are mainly distributed in urban areas, while in bare land and vegetation covered areas, the PS points are sparse, and he deformation position can only be roughly detected from the deformation results, which is difficult to reflect the specific deformation range and boundary (Figure 6(a)). Compared with the PS-InSAR method, the number of measurement points identified by SBAS-InSAR method has increased, but it can be seen from the figure that the range of deformation area is relatively vague, and the boundary of deformation area cannot be clearly determined (Figure 6(a)). The number of measurement points obtained by DS-InSAR method integrating PS points has increased significantly, especially in bare land and vegetation covered areas, and the deformation location and range boundary can be clearly determined according to deformation information. The DS-InSAR monitoring results suggest 4 deformation areas in the Hongta District of Yuxi City, which are distributed in the main city, located near Beicheng Street, Black Village – Yellow Mud Field (2 #), Daxia Street – Near Gaotang Street (3 #), Relien Street (4 #). Among them, the most severe deformation is 1 # region (Figure 6(d)), the Shenzhang District is the centre, north to Meiyuan Village, west to Lianchi Village, south to Liu Zhangying, and East to Honglong Road. The range is large, the deformation is more severe, and the maximum deformation rate is −71.58 mm/a. 2 #, 3 #, and 4 # deformation areas will not be significant, and the deformation rate is within -30mm/a.

Figure 6. Deformation rate map of InSAR monitoring results in various time series. (a) PS-InSAR results, (b) SBAS-InSAR results, (c) DS-InSAR results, and (d) 1# magnified rendering.

In addition, it is worth noting that there is an obvious deformation area in the north of Figure 6(b) (shown by the black rectangle), but the deformation area is not obvious in Figure 6(a) and (c). After analysis, it is considered that the deformation area is affected by mining, and the decoherence is serious in a long time interval. However, both PS-InSAR and DS-InSAR methods with PS points have interference pairs with long time intervals, so these two methods lose their monitoring performance. SBAS-InSAR adopts the combination strategy of short-time baseline, which can effectively obtain the deformation information of this area.

4.2. Key deformation regional analysis

As observed in Figure 6, the surface deformation rate of the study area is mainly distributed between −71.58 ∼ 31.4 mm/a. Near North City Street (1 #) is the largest in the monitoring area, and the deformation is the most remarkable. Therefore, the analysis is focussed on the subsidence area 1#, as shown inside the black dashed rectangle in Figure 7(a), and the local amplified subsidence rate maps of the 1# subsidence area from the monitoring results of PS-InSAR, SBAS-InSAR, and DS-InSAR technique with fused PS points are shown in Figure 7(b–d), respectively. Monitoring results of 1 # deformation area partially enlarged the deformation rate map. Besides, the DS-InSAR monitoring results can reflect the surface deformation details, which are more significant. The white dashed lines and white solid lines are presented in Figure 6(d). The figure suggests that the white dashed lines are more severe, and the deformation rate is less than −30 mm/a. Meanwhile, the deformation between white dashed lines and solid lines is relatively small, and the deformation rate is −30 mm/a ∼ −10 mm/a. There is almost no significant deformation of the white solid line outside.

Figure 7. Deformation rate map of 1# deformation area obtained by InSAR technology in each time series. (a) results of DS-InSAR, (b) results of PS-InSAR in Area 1#, (c) results of SBAS-InSAR in Area 1#, and (d) results of DS-InSAR in Area 1#.

Through statistical analysis of the number and spatial density of typical area measurement points (MPs) in Figure 7(b–d), as shown in Table 2. A total of 33,720 MPs were extracted by PS-InSAR technology, and their spatial density was 450/km². A total of 93,329 MPs were extracted by SBAS-InSAR, and their spatial density was 1244/km². The DS-InSAR technology with PS points extracted a total of 372,188 MPs with a spatial density of 4693/km². The statistical results show that the MPs density extracted by DS-InSAR method with PS points is about 11 times that of PS-InSAR method and 4 times that of SBAS-InSAR method.

Table 2. Number of MPs and their spatial density in the results of the three methods.

	Number of MPs	Spatial density (MPs/km^2)
	Descending	Descending
PS-InSAR	33,720	450
SBAS-InSAR	93,329	1244
DS-InSAR fused PS	372,188	4963

Figure 8 is a segmentation rate cross-sectional view drawn along the A-A 'and B-B' sections shown in Figure 7(d) for further illustrating the settling characteristics of the settling region. Figure 8(a) is a deformation rate along the cross-sectional view of A-A', and Figure 8(b) is a deformation rate along B-B'. It can be seen from Figure 8 that the variation trend of deformation rate obtained by the three methods is highly consistent. As the DS-InSAR method integrating PS points increases the spatial density of MPs, it has great advantages in drawing the details of deformation rate changes. There is an obvious jumping phenomenon at the boundary positions of deformed and non-deformed areas, as shown by the red vertical dotted line in Figure 8. Therefore, the deformation area boundary can be accurately drawn according to the details of deformation rate changes. The other two methods are difficult to reflect the details of surface deformation due to the small number of MPs obtained.

Figure 8. Profile of deformation rate obtained by three InSAR methods. (a) represents section A-A', (b) represents section B-B'.

4.3. Analysis of changes in time and space characteristics

In order to better demonstrate the spatio-temporal characteristics of the temporal deformation of the key subsidence area during January 2019–December 2019, the temporal cumulative deformation of the 1# subsidence area was obtained by taking the primary image (2019-04-22) as a reference and setting its cumulative deformation to 0. The cumulative deformation of the area is exhibited in Figure 9. The figure demonstrates that the cumulative variable and deformation range gradually increases with the increase in the temporal baseline, the cumulative deformation did not show a strictly linear variation with the temporal baseline; the maximum cumulative deformation reaches 57 mm. From 2019-04-22 to 2019-12-30, its cumulative deformation and temporal baseline did not present strict linear changes from 2019-04 to 2019-07. It is faster from 2019-04 to 2019-07, during which the maximum deformation amount is −59 mm. The study area exhibited a significant uplift phenomenon from October 19, 2019 to October 31, 2019 and a significant subsidence phenomenon until December 30, 2019, reaching the maximum cumulative subsidence value of −71 mm. During the monitoring period, the first is in the northern part of the deformation area, and the deformation centre gradually moves in the southwest direction.

Figure 9. Time series cumulative deformation map of 1# deformation area.

Thus, the length of the study area surface deformation and temporal baseline is not a strictly linear relationship, but exhibits seasonal changes. Figure 10 renders the timing of 3 points in the deformation area. The deformation in the study area is the most serious during April 2019 to July 2019 (see red rectangle in Figure 10). The maximum cumulative deformation amount is about −59 mm during this period. A rising phenomenon appeared in August and October 2019 (the black ellipse in Figure 10).

Figure 10. Time series cumulative deformation map of a single point in the deformation area.

4.4. Distributed Scatterer InSAR reliability analysis of fusion PS points

The ground deformation rate of multiple temporal InSAR technologies is correlated to illustrate the reliability of the distributed scatterer InSAR surface deformation monitoring results. The PS-InSAR is used as a reference value, and the nearest neighbour method is adopted. The DS-InSAR monitoring results of SBAS-InSAR and the fusion PS points are searched for in a certain spatial range (50 m). The Measurement PS Point is the measurement point acquired by multiple temporal InSAR technologies and correlates to the deformation rate of the same name (Figure 11). Among them, Figure 11(a) exhibits a PS-InSAR method to monitor the deformation rate coherence analysis results, with a total of 2375 coincident points and a correlation coefficient of 0.9302. Figure 11(b) shows the results of correlation analysis between the PS-InSAR method and the DS-InSAR method monitoring results of fused PS points, and a total of 5429 coincident points were identified with the correlation coefficient of 0.9473. Figure 11(c) shows the results of correlation analysis between the SBAS-InSAR method and the DS-InSAR method monitoring results of fused PS points, and a total of 2094 coincident points were identified with the correlation coefficient of 0.8583. The three InSAR monitoring methods are highly correlated, especially in areas with low deformation rate, while in areas with high deformation rate, the correlation is relatively low, as shown in the red ellipse in Figure 11. The monitoring results of PS-InSAR and SBAS-InSAR show that the deformation rate is obviously low, because the MPs density identified by PS-InSAR and SBAS-InSAR is low, and the distance between two adjacent points is too large, which leads to errors. Additionally, the coherent coefficient between PS-InSAR and the DS-InSAR of the fusion PS points is 0.9473, because the PS point is fused when the DS-InSAR processing of the PS point is fused. The method has a high coherence. The settling results obtained by the three methods have high consistency. However, the DS-InSAR monitoring method of fusion PS points is superior to the PS-InSAR and SBAS-InSAR methods in a large region.

Figure 11. The correlation map of deformation rate of the InSAR monitoring results in each time. (a) the deformation rate correlation between PS-InSAR and SBAS-InSAR, (b) the deformation rate correlation between PS-InSAR and DS-InSAR, and (c) the deformation rate correlation between SBAS-InSAR and DS-InSAR.

Figure 12 statistically illustrates the deformation rate difference histogram of the monitoring results obtained by three timing InSAR methods. As can be seen from the histogram of the difference in deformation rate, the difference in deformation rate between the coincident points of the monitoring results of the three time-series InSAR techniques is small, and the vast majority of them are distributed between -15mm/a and 15 mm/a. The number of coincident points with the absolute value of the difference in deformation rate between PS-InSAR and SBAS-InSAR, DS-InSAR of PS-InSAR and fused PS points, and SBAS-InSAR and fused PS The numbers of coincident points with absolute values of the difference in deformation rate between the DS-InSAR monitoring results of the points within 5 mm/a were 1088, 2710, and 780, accounting for 45.8%, 49.9%, and 37.2% of the number of coincident points screened, respectively. The sequential rate differences in each timing InSAR technology monitoring result suggest that the distributed scatterer InSAR monitoring results have high reliability.

Figure 12. The histogram of difference in the deformation rate of the points with the same mame among the InSAR monitoring results in each time series. (a) the histogram of deformation rate difference between PS-InSAR and SBAS-InSAR, (b) the histogram of deformation rate difference between PS-InSAR and DS-InSAR, and (c) the histogram of deformation rate difference between SBAS-InSAR and DS-InSAR.

From the perspective of data processing, the stability of the Delaunay triangulated irregular network can be significantly increased by combining PS and DS to construct Delaunay triangulated irregular network. Meanwhile, the reliability of phase unwrapping can be improved. Additionally, the quality of the interferogram is improved through the optimization of its phase. This is beneficial to obtain relatively reliable monitoring results of land subsidence.

4.5. Relationship between surface deformation and rainfall in research area

Normally, rainfall is one of the crucial factors leading to surface deformation. The five InSAR measuring points shown in Figure 13(a) in the study area are selected for the image of 2019-01-04 to analyze the formation of the research area and its relationship with the rainfall. The time series accumulated deformation amount is calculated. The accumulated deformation curves in Figure 13(b) and the same period of rainfall are plotted. The rainfall data are collected from GodDard Earth Sciences Data and Information Services Center (https://disc.gsfc.nasa.gov/). Figure 13 unveils that the amount of accumulated deformation at each point exhibits a certain seasonal change, and its change has a certain relationship with rainfall. Besides, the displacement deformation has a significantly accelerated sign before the rainy season of June 2019 and October 2019. The cumulative deformation during the rainy season of June to October reaches its peak, and its deformation speed slows down. After the rainy season, a short-term lifting phenomenon occurs after October 2019, and the subsequent movement significantly changes. The accelerated phenomenon in the rainy season and the time of the rainy season could be investigated based on the above analysis. The displacement speed is slower during the rainy season. This phenomenon could be explained as follows. Heavy rain causes a large amount of rain into the ground, and the soil water content increases as the rain is constantly accumulated. As a result, the soil porosity and the volume increase, and the surface variable speed slows down and even lifts. With the decrease in rainfall, the groundwater gradually loses, and the surface level movement will have a significant acceleration.

Figure 13. The relationship between cumulative deformation and rainfall in the study area. (a) the distribution map of selected InSAR measuring points, and (b) the relationship between cumulative deformation of InSAR measuring points and rainfall.

5. Discussion

In this paper, we apply PS-InSAR, SBAS-InSAR and DS-InSAR method combined with PS points to the same research area, and find that DS-InSAR method combined with PS points can effectively increase the number of measurement points and significantly improve the density of MPs. It also allows us to define and update the boundaries of the deformed and non-deformed regions more accurately.

Compared with PS-InSAR and SBAS-InSAR methods, the DS-InSAR method combining PS points can obtain more measurement points and better monitoring results. However, the pre-processing of this method is quite time-consuming and the computation time cost is about 4-5 times of other time-series InSAR methods, which is the main bottleneck of its application. For surface deformation monitoring over large areas, the time cost of DS-InSAR analysis may increase to a very high level. Therefore, how to utilize modern massively parallel computing techniques, such as graphics processing units (GPUs), to improve computational efficiency becomes a critical issue to be addressed. Another key issue worth investigating is how to optimize the DS-InSAR algorithm, especially the DS preprocessing process, in order to significantly reduce the storage space consumption.

6. Conclusion

In this paper, the sentinel-1A data is selected, and the distributed scatterer InSAR method of fusion PS is used to monitor and analyze the surface deformation of Yuxi City Hongta District in January 2019, and the reliability of the method in this paper was verified by the conventional PS-InSAR and SBAS-InSAR methods. The conclusions can be drawn as follows.

1. There are 4 main deformation areas in the Hongta District of Yuxi City, which is located in the deformation area near North City Street. The maximum deformation rate reaches −71.58 mm/a, as indicated in the cross-sectional analysis. The significant jump phenomenon reveals that the segmentation of the region may be induced by surface collapse or tomography. The deformation rates of the remaining three settling regions are relatively small, all within −30 mm/a.

2. At the time of time-sequential variable spatial feature analysis, the research area surface deformation showed seasonal changes, and the deformation in April and July 2019 was the most severe. The maximum deformation reached 59 mm. The uplift occurred twice in August 2019 and October 2019, and the slow subsidence occurred in other periods. Over time, the location of the subsidence centre moves from northeast to southwest.

3. The distributed scatterer InSAR method of fusion PS can be employed to obtain higher measurement point density than conventional PS-InSAR and SBAS-InSAR methods, especially compared to the PS-InSAR method. Moreover, its measurement point density is boosted. Thereby, the boundary between the deformation area and the non-deformation area can be delineated more precisely.

4. The monitoring results of the fused distributed scatterer InSAR method and conventional PS-InSAR and SBAS-InSAR methods are highly consistent on spatial position distribution, implying that their monitoring results are high. The correlation analysis demonstrates that in a region with a small deformation rate, the correlation is high, and the region-dependent correlation rate in the deformation rate is relatively low. Moreover, the deformation rate of the two methods of PS-InSAR and SBAS-InSAR are significantly smaller. This further verifies that the method can be utilized to more accurately acquire the deformation rate

Although the study area formation information has been successfully obtained, there are some shortcomings in this study. First, the conditions are limited, and the reliability and accuracy of this method are verified only by the monitoring results of the conventional timing InSAR method, instead of the measured data. Second, the results in this paper have been monitored and analyzed only in the study area, and its deformation mechanism is neglected. In order to better analyze the deformation mechanism of the study area, the corresponding deformation mechanism model will be established subsequently by combining the knowledge of other disciplines such as geohydrology and environment. To provide more accurate prediction and service to prevent potential geological hazards such as cave-in in the study area.

Acknowledgments

The sentinel-1A datasets were freely provided by the European Space Agency (ESA) and freely downloaded from the website https://vertex.daac.asf.alaska.edu. The NASA for providing the SRTM-1 DEM data and freely downloaded from the website http://rmw.recordist.com/index.html. The Precipitation Orbit Data (POD) was freely downloaded from the website https://qc.sentinel1.eo.esa.int/. The Rainfall data were freely downloaded from the website https://disc.gsfc.nasa.gov/.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

The data that support the findings of this study are available from the first author, Yongfa Li, upon reasonable request.

Additional information

Funding

This study was supported by National Natural Science Foundation of China(42161067 and 42004006).