Abstract

The implementation of the South-to-North Water Diversion Project (SNWDP) has alleviated groundwater resource pressure in North China to some extent, resulting in a gradual deceleration of land subsidence and even rebound in some areas. To investigate the spatiotemporal evolution characteristics of land subsidence in the eastern plain of Beijing following the SNWDP, this study employs Ascending (ASC) and Descending (DES) InSAR data combined with a Strain Model (SM) to obtain a Three-Dimensional (3-D) deformation field from 2016 to 2018. Through analysis of the 3-D deformation characteristics and spatiotemporal evolution of land subsidence in this region from 2016 to 2018, the results reveal a shift in the distribution of subsiding areas after the South-to-North Water Diversion, with a marked decrease in subsidence rates in certain areas. The maximum subsidence rate in the Beijing area has decreased to 110 mm/yr, accompanied by horizontal deformation at a rate of 12 mm/yr. Additionally, by examining the spatial relationship between major active faults and subsidence deformation in this region, the study further elucidates the influence of fault activity on the spatial distribution of subsidence deformation.

Keywords

Land Subsidence, 3-D Deformation, Strain Model, InSAR

Share and Cite:

1. Introduction

Land subsidence is an environmental geological phenomenon characterized by a loss of ground elevation due to soil compression, leading to irreversible environmental and resource damage. Currently, more than 200 cities and regions worldwide are facing the challenge of land subsidence [1]-[6]. The threats posed by land subsidence to urban areas mainly fall into two categories: first, uneven subsidence damages infrastructure such as urban buildings and linear engineering works [7]; second, the loss of ground elevation exacerbates the risk of urban flooding [8] [9].

Land subsidence is influenced by a combination of natural and human factors. Natural factors include tectonic movements, natural soil consolidation, earthquakes, and climate change. Human factors involve the extraction of subsurface fluids (such as groundwater, geothermal fluids, and oil), mining of solid minerals (such as coal and metal ores), and the additional stress exerted on soil layers by surface engineering projects, which leads to soil layer compression and consolidation [10] [11]. Over-extraction of groundwater is considered the primary human-driven factor contributing to land subsidence in Beijing [12].

With the official opening of the South-to-North Water Diversion Project (SNWDP) and the implementation of groundwater extraction restrictions and bans, Beijing’s groundwater usage has significantly decreased, leading to a noticeable change in water use patterns and alleviate groundwater pressure. Groundwater levels in Beijing have stabilized and started to rise, subsidence rates have markedly declined, and the spatial distribution of subsidence has been gradually shrinking year by year, with even signs of uneven rebound in some areas. Therefore, investigating the spatiotemporal evolution characteristics of the 3-D deformation field of land subsidence under this new hydrological context is of great importance for understanding the occurrence, development, and evolution of land subsidence in Beijing in the next stage.

The main methods for monitoring land subsidence include leveling surveys, GNSS measurements, and Synthetic Aperture Radar Interferometry (InSAR) [13]-[15]. With the advancement of SAR satellite technology, InSAR has significant advantages in monitoring land subsidence over large areas and at high spatial resolutions. Additionally, InSAR offers certain advantages in measuring 3-D deformation [16]-[19]. At present, the approach to obtaining 3-D deformation using InSAR primarily relies on multi-orbit data, specifically combining ascending and descending orbits. This method constructs a system of observation equations by correlating data from the ascending (ASC) and descending (DES) orbits, or from multiple data sources. The 3-D deformation is then determined by applying the weighted least squares (WLS) method [20]. For example, combining ASC and DES orbits with GNSS and leveling data is currently the most direct measurement approach in engineering, yielding reliable deformation results. However, due to high costs, the density of station distribution often impacts the accuracy of measurements. Leveling surveys are also affected by limitations in manpower and financial resources. Another approach is to combine ASC-DES orbits and prior knowledge [21]-[27]. This approach uses specific assumptions or models tailored to the research context, minimizing the need for extra data when calculating 3-D deformation. However, these methods generally focus on rapid, severe deformation due to sudden geological events (such as earthquakes, landslides, mining subsidence, or iceberg movement), and pay less attention to gradual, long-term deformation. Therefore, this paper adopts a universal model, the Strain Model (SM), and combines it with InSAR ASC and DES orbits [28]-[32]. This approach acknowledges the interdependence among particles within a defined neighborhood radius on a continuous surface, with SM representing the strain relationship between particles to derive the 3-D deformation. This method partially addresses the issue of pixel-wise independence inherent in the WLS pixel-by-pixel or point-based solutions, providing broader applicability across various scenarios.

Based on the above, we first selected Sentinel-1A ASC data and Sentinel-1B DES data from 2016-2018 to obtain PS-InSAR results for each track. Then, we used the ASC and DES data in combination with the SM to jointly resolve the 3-D deformation field of Beijing after the South-to-North Water Diversion. The accuracy of the results was verified by comparing them with leveling and GNSS points. Finally, the relationship between 3-D deformation and faults is discussed by studying several major active faults in the region in order to investigate the influence of the distribution of faults on land subsidence.

2. Study Area and Dataset

The study area is located in the eastern plain region of Beijing, at the intersection of the Yanshan, North China Plain, and Taihang fault systems. Previous research has shown that extensive groundwater extraction led to a maximum subsidence rate of up to 14 cm/year from 2015 to 2017. The specific location of the study area is shown in Figure 1.

The SAR data utilized in this study comprises ASC and DES orbit data (tracks 142 and 47) acquired by Sentinel-1 A/B in interactive wide swath (IW) mode for the period 2016-2018. Specifically, 51 DES orbit datasets serve as the primary data, while 44 ASC datasets act as auxiliary data to validate the study’s findings, with both using vertical-vertical (VV) polarization. Detailed information on these datasets can be found in Table 1.

Table 1. Satellite information for the SAR data in this study.

Figure 1. Geographical location of the study area. (a) Sentinel-1 SAR image coverage and the distribution location of the validation data. (b) Approximate location of the study area. (c) The location of North China Plain and Beijing-Tianjin-Hebei Area in China.

This study also selected 15 leveling points and 13 GNSS points to verify vertical and horizontal deformation, respectively. By establishing buffer zones, PS points were compared with leveling and GNSS points to ensure the reliability of the deformation results.

3. Methodology

3.1. Construction of an ASC-DES Solution Model

In radar side-view geometry, InSAR observations capture a vector sum of the ground’s projections in three directions: east, north, and vertical. This allows the decomposition of InSAR’s line-of-sight (LOS) deformation into components along these directions. Representing these components as $d_e$ $d_n$ and $d_u$ for the east, north, and vertical directions respectively, the relationship between 3-D deformation and LOS deformation can be formulated as [33]

$d_{LOS}=d_u\cos \theta -\sin \theta \left[d_n\cos \left(\alpha -\frac{3}{2}\pi \right)+d_e\sin \left(\alpha -\frac{3}{2}\pi \right)\right]$ (1)

where $\theta$ is the satellite incidence angle, i.e., the angle between LOS direction and vertical direction; $\alpha$ is the satellite azimuth angle, i.e., the clockwise angle between the satellite heading and due north direction. To reduce the number of independent observations, the left-right view mode can be applied to decrease observation counts. This approach enables the construction of four equations from ASC-DES orbit data. By using the left-right view direction model, the azimuthal deformation coefficients are considered as positive and negative, respectively [29]

$\{\begin{array}{l}{d}_{LO{S}_{ASC}}=d_u\cos {\theta }_{ASC}-flag\cdot \sin {\theta }_{ASC}\left[d_n\cos \left({\alpha }_{ASC}-\frac{3}{2}\pi \right)+d_e\sin \left({\alpha }_{ASC}-\frac{3}{2}\pi \right)\right]\\ d_{LO{S}_{ASC}}=d_u\cos {\theta }_{ASC}-flag\cdot \sin {\theta }_{ASC}\left[d_n\cos \left({\alpha }_{ASC}-\frac{3}{2}\pi \right)+d_e\sin \left({\alpha }_{ASC}-\frac{3}{2}\pi \right)\right]\\ d_{LO{S}_{DES}}=d_u\cos {\theta }_{DES}-flag\cdot \sin {\theta }_{DES}\left[d_n\cos \left({\alpha }_{DES}-\frac{3}{2}\pi \right)+d_e\sin \left({\alpha }_{DES}-\frac{3}{2}\pi \right)\right]\\ d_{LO{S}_{DES}}=d_u\cos {\theta }_{DES}-flag\cdot \sin {\theta }_{DES}\left[d_n\cos \left({\alpha }_{DES}-\frac{3}{2}\pi \right)+d_e\sin \left({\alpha }_{DES}-\frac{3}{2}\pi \right)\right]\end{array}$ (2)

where ASC represents the ascending orbit observation, DES represents the descending orbit observation, and $flag$ is as follows

$flag=\{\begin{array}{l}1,\text{left-looking}\text{mode}\\ -1,\text{right-looking}\text{mode}\end{array}$ (3)

Converting the above equation into matrix form i.e.

$Ad=b$ (4)

where A is the coefficient matrix; d is the 3-D deformation ${\left[d_e{d}_n{d}_u\right]}^{\text{T}}$ ; b is the LOS value and its components. The coefficient matrix A and the variable b can be written as

$A=\left[\begin{array}{c}\begin{array}{ccc}-\sin {\theta }_{ASC}\sin \left(a_{ASC}-\frac{3}{2}\pi \right)& -\sin {\theta }_{ASC}\cos \left(a_{ASC}-\frac{3}{2}\pi \right)& \cos {\theta }_{ASC}\\ \sin {\theta }_{ASC}\sin \left(a_{ASC}-\frac{3}{2}\pi \right)& \sin {\theta }_{ASC}\cos \left(a_{ASC}-\frac{3}{2}\pi \right)& \cos {\theta }_{ASC}\end{array}\\ \begin{array}{ccc}-\sin {\theta }_{DES}\sin \left(a_{DES}-\frac{3}{2}\pi \right)& -\sin {\theta }_{DES}\cos \left(a_{DES}-\frac{3}{2}\pi \right)& \cos {\theta }_{DES}\\ \sin {\theta }_{DES}\sin \left(a_{DES}-\frac{3}{2}\pi \right)& \sin {\theta }_{DES}\cos \left(a_{DES}-\frac{3}{2}\pi \right)& \cos {\theta }_{DES}\end{array}\end{array}\right]$ (5)

$b=\left[\begin{array}{c}{d}_{LO{S}_{ASC}}\\ d_{LO{S}_{ASC}}\\ d_{LO{S}_{DES}}\\ d_{LO{S}_{DES}}\end{array}\right]$ (6)

3.2. Combining Strain Model with ASC-DES Solution Model

The Strain Model is a classical model in geophysics that reflects in detail the strain relationships between near-field points [34]. Suppose the 3-D coordinates and 3-D deformation of the central point $P^0$ are $x^0=\left[x^{e\left(0\right)}{x}^{n\left(0\right)}{x}^{u\left(0\right)}\right]$ and $d^0=\left[d^{e\left(0\right)}{d}^{n\left(0\right)}{d}^{u\left(0\right)}\right]$ , respectively. At the same time, a neighborhood is created with a radius of 1km from the point $P^0$ , and there are K points $P^i\left(i=1,2,3,\cdots ,k\right)$ in the neighborhood. Their 3-D coordinates and 3-D deformations are $x^i=\left[x^{e\left(i\right)}{x}^{n\left(i\right)}{x}^{u\left(i\right)}\right]$ and $d^i=\left[d^{e\left(i\right)}{d}^{n\left(i\right)}{d}^{u\left(i\right)}\right]$ respectively. $e,n,u$ represent the East-West (E-W), North-South (N-S) and Up-Down (U-D) directions respectively. Suppose ${\Delta }^i=x^i-x^0=\left[\Delta x^{e\left(i\right)}{x}^{n\left(i\right)}{x}^{u\left(i\right)}\right]$ is the difference in coordinates (relative coordinates) between the points $P^0$ and $P^i$ , then the deformation relations $d^0$ and $d^i$ between these two points can be expressed as

$d^i=H\cdot {\Delta }^i+d^0$ (7)

where $H=S+R$ is the displacement gradient matrix, and the symmetric and antisymmetric parts S and R can be written as

$S=\left[\begin{array}{ccc}{\xi }_{ee}& {\xi }_{en}& {\xi }_{eu}\\ {\xi }_{en}& {\xi }_{nn}& {\xi }_{nu}\\ {\xi }_{eu}& {\xi }_{nu}& {\xi }_{uu}\end{array}\right]$ (8)

$R=\left[\begin{array}{ccc}0& -{\omega }_{en}& {\omega }_{eu}\\ {\omega }_{en}& 0& -{\omega }_{nu}\\ -{\omega }_{eu}& {\omega }_{nu}& 0\end{array}\right]$ (9)

where $\xi$ and $\omega$ represent the parameters of the strain tensor and the rigid body rotation tensor, respectively. According to Equation (6) and (7), Equation (5) can be rewritten as [35]

$d^i=B_{sm}^i\cdot l$ (10)

where $l={\left[d^{e\left(0\right)}{d}^{n\left(0\right)}{d}^{u\left(0\right)}{\xi }_{ee}{\xi }_{en}{\xi }_{eu}{\xi }_{nn}{\xi }_{nu}{\xi }_{uu}{\omega }_{en}{\omega }_{eu}{\omega }_{nu}\right]}^{\text{T}}$ is the unknown parameter at the point $P^0$ ; $B_{sm}^i$ is the transformation matrix of $d^i$ and $l$ , which can be expressed as

$B_{sm}^i= {\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\begin{array}{ccc}\Delta x^{e\left(k\right)}& \Delta x^{n\left(k\right)}& \Delta x^{u\left(k\right)}\\ 0& \Delta x^{e\left(k\right)}& 0\\ 0& 0& \Delta x^{e\left(k\right)}\end{array}\begin{array}{ccc}0& 0& 0\\ \Delta x^{n\left(k\right)}& \Delta x^{u\left(k\right)}& 0\\ 0& \Delta x^{n\left(k\right)}& \Delta x^{u\left(k\right)}\end{array}\begin{array}{ccc}-\Delta x^{n\left(k\right)}& -\Delta x^{u\left(k\right)}& 0\\ \Delta x^{e\left(k\right)}& 0& -\Delta x^{u\left(k\right)}\\ 0& -\Delta x^{e\left(k\right)}& \Delta x^{n\left(k\right)}\end{array}\right]}_{3k\ast 12}$ (11)

In Section 2.1, we construct the multi-orbit solution model and obtain the coefficient matrix A and the LOS matrix b. Since there is only one orbital LOS deformation result, we replace the deformation values in matrix b with the LOS deformation values of the current orbitals owned, i.e., all of them are ascending or descending orbitals, and matrix A and matrix b are rewritten as follows

$A_{geo}^i={\left[\begin{array}{ccc}-\sin {\theta }_{ASC}\sin \left(a_{ASC}-\frac{3}{2}\pi \right)\cdot I^k& -\sin {\theta }_{ASC}\cos \left(a_{ASC}-\frac{3}{2}\pi \right)\cdot I^k& \cos {\theta }_{ASC}\cdot I^k\\ \sin {\theta }_{ASC}\sin \left(a_{ASC}-\frac{3}{2}\pi \right)\cdot I^k& \sin {\theta }_{ASC}\cos \left(a_{ASC}-\frac{3}{2}\pi \right)\cdot I^k& \cos {\theta }_{ASC}\cdot I^k\\ -\sin {\theta }_{DES}\sin \left(a_{DES}-\frac{3}{2}\pi \right)\cdot I^k& -\sin {\theta }_{DES}\cos \left(a_{DES}-\frac{3}{2}\pi \right)\cdot I^k& \cos {\theta }_{DES}\cdot I^k\\ \sin {\theta }_{DES}\sin \left(a_{DES}-\frac{3}{2}\pi \right)\cdot I^k& \sin {\theta }_{DES}\cos \left(a_{DES}-\frac{3}{2}\pi \right)\cdot I^k& \cos {\theta }_{DES}\cdot I^k\end{array}\right]}_{4k\ast 3k}$ (12)

$b^i= {\left[\begin{array}{c}{d}_{ASC\left(DES\right)}^k\\ d_{ASC\left(DES\right)}^k\\ d_{ASC\left(DES\right)}^k\\ d_{ASC\left(DES\right)}^k\end{array}\right]}_{4k\ast 1}$ (13)

where $I^k$ is the unit matrix of order k. $b^i$ represents the LOS deformation result for the current orbit. According to Equation (2), we can get

$A_{geo}^i\cdot d^i= b^i$ (14)

Combine Equation (8)-(12) to obtain

$b^i=A_{geo}^i\cdot B_{sm}^i\cdot l$ (15)

Let $B^i=A_{geo}^i\cdot B_{sm}^i$ to obtain

$b^i= B^i\cdot l$ (16)

Then according to the WLS method we can get $l$ (i.e., the 3-D deformation value of the current point $P^0$ ), $l$ can be written as

$l={\left(B^{\text{T}}WB\right)}^{-1}{B}^{\text{T}}Wb$ (17)

where W is a 4k * 4k weight matrix. Based on previous research results [36] [37], the expression is rewritten in this paper as

$W=\exp \left(\frac{s_i}{r}\right)\cdot I$ (18)

where $s_i$ is the Euclidean distance from point $P^i$ to point $P^0$ ; r is the radius of the neighborhood; and I is the 4k * 4k unit matrix.

Figure 2. Difference of different densities in a specific neighborhood. (a) The starting condition is satisfied and the number of points in the neighborhood is greater than 200. (b) The starting condition is satisfied and the number of points in the neighborhood is less than 200. (c) The starting conditions are not satisfied.

Due to the varying density of PS points obtained through PS-InSAR, the number of PS points within a specified neighborhood radius differs for each point, making it challenging to use a fixed principle to establish weights [38]. Therefore, this study adopted the most straightforward approach, using Euclidean distance for weighting. According to previous research, the SM model can yield stable parameter solutions when the number of pixels within a radius of 1 - 2 km exceeds 200 [35] [39]-[41]. However, PS point density may vary, and in some areas, the specified radius may not contain 200 PS points, resulting in unsolvable conditions, as shown in Figure 2.

In response to the above, we use a linear model to solve the problem that the PS points cannot satisfy the threshold. Assuming that the deformation occurs along a specific direction determined by the unit vector $\widehat{s}=\left(s_e,s_n,s_u\right)$ . Therefore, by specifying the additional equation that defines the parallelism between d and $\widehat{s}$ , i.e., $d\times \widehat{s}=0$ (with 0 being the null vector in the geometric space), Equations (5) and (6) can be re-expressed as

$\begin{array}{c}A=\left[\begin{array}{ccc}-\sin {\theta }_{ASC}\sin \left(a_{ASC}-\frac{3}{2}\pi \right)& -\sin {\theta }_{ASC}\cos \left(a_{ASC}-\frac{3}{2}\pi \right)& \cos {\theta }_{ASC}\\ \sin {\theta }_{ASC}\sin \left(a_{ASC}-\frac{3}{2}\pi \right)& \sin {\theta }_{ASC}\cos \left(a_{ASC}-\frac{3}{2}\pi \right)& \cos {\theta }_{ASC}\\ -\sin {\theta }_{DES}\sin \left(a_{DES}-\frac{3}{2}\pi \right)& -\sin {\theta }_{DES}\cos \left(a_{DES}-\frac{3}{2}\pi \right)& \cos {\theta }_{DES}\\ \sin {\theta }_{DES}\sin \left(a_{DES}-\frac{3}{2}\pi \right)& \sin {\theta }_{DES}\cos \left(a_{DES}-\frac{3}{2}\pi \right)& \cos {\theta }_{DES}\\ 0& s_u& -s_n\\ -s_u& 0& s_e\\ s_n& -s_e& 0\end{array}\right]\end{array}$ (19)

$b=\left[\begin{array}{c}{d}_{LO{S}_{ASC}}\\ d_{LO{S}_{ASC}}\\ d_{LO{S}_{DES}}\\ d_{LO{S}_{DES}}\\ 0\\ 0\\ 0\end{array}\right]$ (20)

Equivalently, assuming that $d_s$ is the (unknown) component of the displacement vector along the $\widehat{s}$ direction, the model constraints can also be expressed as

$d-d_s\widehat{s}=0$ (21)

The unknown displacement $d_s$ is

$a{d}_s=b$ (22)

$a=\left[\begin{array}{c}{\widehat{l}}_1\cdot \widehat{s}\\ {\widehat{l}}_2\cdot \widehat{s}\\ {\widehat{l}}_3\cdot \widehat{s}\\ {\widehat{l}}_4\cdot \widehat{s}\end{array}\right]=\left[\begin{array}{c}\cos {\beta }_1\\ \cos {\beta }_2\\ \cos {\beta }_3\\ \cos {\beta }_4\end{array}\right]$ (23)

where $\beta$ is the angle between the OS direction and $\widehat{s}$ . From this, it follows that

$d_s={\left(a^{\text{T}}a\right)}^{-1}{a}^{\text{T}}y=\frac{{{\displaystyle \sum }}_{m=1}^M{b}_m\cos {\beta }_m}{{{\displaystyle \sum }}_{m=1}^M{\cos }^2{\beta }_m}$ (24)

3.3. PS-InSAR Processing Flow

This study utilizes Sentinel-1 satellite data provided by the European Space Agency (ESA) and applies the Persistent Scatterer Interferometry (PS-InSAR) technique to monitor ground deformation in the study area. The specific workflow includes image preprocessing, interferogram generation, Persistent Scatterer (PS) point extraction, atmospheric phase removal, and deformation estimation. The data was processed and analyzed using the SARPROZ software, with the detailed workflow as follows: 1) Data Preparation and Preprocessing; 2) Image Co-Registration and Interferogram Generation; 3) Persistent Scatterer Point Extraction; 4) Atmospheric Phase Removal and Deformation Estimation; 5) Results Visualization and Validation.

4. Results and Analysis

4.1. InSAR Results

The results of the LOS deformation for the ASC and DES orbits are shown in Figure 3. The PS point density in the study area is approximately 50.75 points/km². Negative values indicate subsidence, while positive values indicate uplift. The InSAR results delineate the spatial distribution of surface deformation, revealing significant subsidence rates, with a maximum of 110 mm/yr in eastern Beijing.

Figure 3. LOS deformation result. (a) is the result of the DES orbit LOS deformation. (b) is the result of ASC orbit LOS deformation.

As shown in Figure 4, we calculated the deformation results in the E-W, N-S, and U-D directions and provided a display of the annual average rate of the 3-D deformation field. In the east-west direction, the colors range from red to blue, representing positive deformation (displacement of the surface to the east) and negative deformation (displacement of the surface to the west), respectively. In the north-south direction, the colors range from yellow to blue, representing positive deformation (surface displacement to the north) and negative deformation (surface displacement to the south). In the U-D direction, the colors also range from red to blue, corresponding to surface uplift and surface subsidence, respectively. The figure visualizes the significant spatial differences in the rate of surface deformation in each region.

Figure 4. 3-D deformation results. (a) E-W direction deformation. Positive values indicate east, and negative values indicate west direction. (b) N-S direction deformation. Positive values indicate south, and negative values indicate north direction. (c) U-D direction deformation. Positive values indicate uplift, and negative values indicate subsidence.

In order to verify the accuracy of the InSAR deformation results, GNSS measurement points and leveling measurement points are superimposed in the figure. These points provide the basic data support for verifying the vertical and horizontal deformation. The GNSS measurement points are mainly used to verify the horizontal deformation, while the level measurement points are used to verify the vertical deformation. This comparison of multi-source data can better ensure the reliability and accuracy of the InSAR observation results.

4.2. Accuracy Verification by Leveling and GNSS Points

To assess the accuracy of the InSAR three-dimensional deformation results, we selected GNSS and leveling points for validation. Leveling points are used to validate vertical deformation, while GNSS points are used to validate horizontal deformation, as leveling measurements are more sensitive to vertical changes, whereas GNSS is more sensitive to horizontal changes [42]. The validation results are shown in Figure 5, Figure 6. We used the cumulative deformation values from 2016 to 2018 to validate vertical deformation, resulting in a root mean square error (RMSE) of 12.7 mm, with an average annual error of 4.2 mm/year. For horizontal deformation, we used the average displacement rates from 2016 to 2018 for validation, obtaining an RMSE of 2.3 mm/year in the east-west direction and 3.3 mm/year in the north-south direction. In the case of north-south deformation, the R² value is relatively low, and the deformation fluctuations are not significant. This is mainly due to InSAR’s insensitivity to north-south deformation, which is a result of the physical imaging model of the SAR satellite [43]. In the validation of vertical deformation, some points had errors exceeding 10 mm, which may be attributed to the temporal mismatch between leveling points and persistent scatterer (PS) points. Additionally, the leveling and GNSS measurements themselves may also have inherent errors. InSAR has a high spatial resolution, while GNSS and leveling data are relatively sparse in sampling points. Due to the differences in spatial resolution between these two measurement methods, it is often necessary to establish buffer zones for data comparison, which may introduce errors. Furthermore, the acquisition times for InSAR data, GNSS, and leveling measurements are often inconsistent. For deformation such as land subsidence that change dynamically over time, differences in data collection times can lead to measurement errors. The inherent errors of different measurement methods are also critical factors. Leveling and GNSS measurements may be affected by operator errors, introducing additional measurement inaccuracies. Additionally, InSAR is susceptible to atmospheric delays, which can further affect measurement accuracy.

To better illustrate the advantages of the method results in this paper, we compared several other time-series 3-D deformation methods. Many of the existing 3-D deformation methods are based on DInSAR, making direct comparisons with MT-InSAR difficult. Consequently, we selected various mainstream time-series 3-D deformation estimation methods for comparison. These include utilizing Leveling data as constraints [18], combining GNSS observations [44], combining ASC and DES orbits using SM and Kalman filtering methods [45], using multi-source SAR satellites [46]. The comparison results are summarized in Table 2.

Table 2. RMSEs of different methods for different kinds of deformations (units: mm). The horizontal bars in the table indicate results that were not obtained for this experiment.

The analysis reveals that the KFInSAR method demonstrates the highest accuracy in U-D deformation, surpassing even the method constrained by Leveling data. This is because Kalman filter is well-suited to handle temporal discontinuities in the data. By predicting and updating the state of deformation iteratively, KFInSAR can mitigate the impact of missing or irregular data, maintaining the reliability of the results. In contrast, the combined Leveling data method need to take into account the spatial discontinuities of Leveling data, and careful consideration of the point distribution is required when combining with InSAR data. The combined GNSS method excels in N-S and E-W deformations due to the higher accuracy of GNSS horizontal deformations than vertical deformations. However, the problem of point space discontinuity in GNSS still exists in the process of combining with InSAR. The MSBAS method employs direct solutions derived from InSAR based on eight independent observations and obtained good results. But it involves a huge number of independent observations and high costs.

Figure 5. Vertical deformation verification results. R² is the correlation coefficient.

Figure 6. Horizontal deformation verification results. R² is the correlation coefficient.

5. Discussion

5.1. Characteristics of 3-D Surface Deformation Evolution in the Eastern Plains of Beijing

In the background of the South-to-North Water Diversion Project, groundwater resources in North China have been partially replenished. However, due to the long history of over-exploitation, the issue of land subsidence in this region remains severe. Figure 4 shows the distribution maps of ground deformation rates obtained using InSAR technology between 2016 and 2018, representing the line-of-sight (LOS) deformation and the decomposed three-dimensional deformation rates (E-W, N-S, and U-D directions), respectively.

In Figure 4, the deformation rate in the E-W direction is concentrated around the area between Beijing and Tianjin, with the region near Beijing showing the most significant deformation, forming a differential deformation field that trends from east to west. This indicates variations in lateral stress in the ground surface. The N-S direction exhibits relatively smaller deformation; however, notable differences are still observed in localized areas, particularly a slight trend of southward deformation near Beijing.

These results indicate that during the implementation of the South-to-North Water Diversion Project, while there has been some alleviation of groundwater pressure, Tianjin and its surrounding areas remain the most significantly subsiding regions due to historical issues related to groundwater extraction. The subsidence rate in this area exhibits a pronounced spatial heterogeneity, which may be closely related to factors such as stratigraphic conditions, soil structure, uneven recovery of groundwater levels, and human activities.

Since the introduction of water from the South-to-North Water Transfer Project to Beijing at the end of 2014, groundwater extraction and usage have gradually decreased. The proportion of annual water supply sourced from groundwater dropped from 52% in 2014 to 45% in 2016. Bringing in southern water and replacing private wells has helped mitigate the over-extraction of groundwater, thereby slowing the progression of land subsidence to some extent. However, the rate of subsidence in this region remains high [18]. The primary contributors to land subsidence in Beijing are the second and third confined aquifers. Horizontal deformation is primarily due to the tensile diffusion of vertical deformation from land subsidence, driven mainly by these two confined aquifers. Additionally, the spatial distribution of horizontal displacement may result from the combined effects of groundwater dynamics and other influencing factors.

5.2. 3-D Deformation and Major Active Faults

This section explores the relationship between three-dimensional deformation and faults by analyzing several major active faults in the Beijing area. Figure 7 shows six major active faults (F1 to F6) in the study area since the Holocene. To obtain a continuous three-dimensional deformation distribution, we applied Kriging interpolation to rasterize the data [47]. The results indicate a high correlation between the morphology of the subsidence area in eastern Beijing and the distribution of faults, suggesting that the fault zones significantly influence the spatial distribution characteristics of ground subsidence. Figure 8 displays the deformation rates along six profiles (A-A’, B-B’, ..., F-F’) in three directions.

In the E-W direction, F1 shows slight fluctuations, with a maximum displacement of approximately 1.5 mm/yr. The N-S direction remains relatively stable, while the U-D direction exhibits significant vertical subsidence, with a maximum subsidence rate of 12 mm/yr. For F2, a notable displacement gradient appears in the E-W direction, with a maximum displacement close to 3.5 mm/yr, indicating strong tensile stress in the east-west direction. The N-S direction shows slight fluctuations, while the U-D direction experiences significant subsidence, reaching a maximum of approximately 120 mm/yr. In F3, the E-W direction remains relatively stable, while there is a slight deformation trend in the N-S direction. The U-D direction shows significant subsidence, with a maximum value of about 60 mm/yr, indicating that subsidence in this area is mainly concentrated in the vertical direction. F4 displays distinct deformation in the E-W direction, with a maximum displacement of approximately 3.0 mm/yr, indicating prominent tensile stress in the east-west direction. The N-S direction shows minimal deformation, while the U-D direction exhibits significant subsidence, with a maximum rate of around 100 mm/yr. F5 shows noticeable fluctuations in the E-W direction, with a maximum displacement of approximately 4.0 mm/yr, suggesting strong horizontal stress in the east-west direction. The N-S direction shows slight fluctuations,

Figure 7. Distribution of the location of faults F1 - F6. The red line represents the fault zone. The yellow line represents the location of the profile.

Figure 8. Profiles of F1 - F6. The vertical axis represents the deformation rate and the horizontal axis represents the distance along profile. The green dashed line represents the location of the fault. The black dashed line represents the beginning and end position of the profile.

while the U-D direction experiences considerable subsidence, with a subsidence rate of about 80 mm/yr. F6 remains relatively stable in the E-W direction, while the N-S direction shows a noticeable displacement gradient. The U-D direction experiences pronounced subsidence, with a maximum subsidence rate of approximately 20 mm/yr.

According to Figure 7 and Figure 8, faults F2, F3, and F5 are located in the surrounding area of the eastern Beijing subsidence zone. Notably, F5 is clearly situated between two subsidence belts. In the directions of F2 and F3, linear boundaries of the subsidence belts can be observed, indicating a consistent spatial distribution of these faults. Nearby these faults, there is significant vertical deformation as well as horizontal tensile strain. In contrast, the vertical deformation and horizontal rates on either side of faults F1, F4, and F6 are not significant, but there is a notable relationship of positive and negative variations, demonstrating opposite movement trends on either side of the faults.

6. Conclusions

This study utilized the combined method of ascending-descending orbits and the Strain Model (SM) to obtain the time-series 3D deformation of the eastern Beijing area from 2016 to 2018. The results indicate that the root mean square error (RMSE) of the 3D deformation in the east-west, north-south, and vertical directions is 2.04 mm/yr, 3.3 mm/yr, and 4.02 mm/yr, respectively, confirming the reliability of the 3D deformation results presented in this study. Additionally, major active faults in the eastern Beijing Plain were monitored, and the following conclusions were drawn from the analysis:

1) Following the South-to-North Water Diversion Project, the ground subsidence rate in the eastern Beijing Plain has significantly slowed, with the maximum subsidence rate decreasing from 140 mm/yr to 110 mm/yr, accompanied by 12 mm/yr of horizontal deformation. In the upper Chaobai River region, where ecological water replenishment is substantial, there has been a marked rebound in ground elevation.

2) The distribution of faults has a notable impact on the spatial distribution of ground subsidence, particularly near active fault zones. Subsidence rates around fault zones are typically higher and exhibit a linear boundary along the fault direction, suggesting that fault zones may influence the compression and deformation behavior of the soil layers.

Acknowledgements

The authors acknowledge the open-source software SNAP for providing us with sentinel data incidence angle information. This work was supported by the National Natural Science Foundation of China (NO. 41930109/D010702, NO. 42371081/D0104, NO. 42371089/D0104, 42201081/D0104) and the Beijing Outstanding Young Scientist Program (BJJWZYJH01201910028032).

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References