Three-Dimensional Surface Deformation Characteristics Based on Time Series InSAR and GPS Technologies in Beijing, China

Abstract

Excessive exploitation of the groundwater has resulted in obvious three-dimensional (3D) deformation features on the surface of the Beijing Plain. This paper, by combining Interferometric Synthetic Aperture Radar (InSAR) and Global Positioning System (GPS) technologies, has obtained time-series information of the 3D surface deformation in the Beijing Plain, analyzing its spatial distribution characteristics. On this basis, the relationship between different controlling factors with the 3D deformation of the surface has been analyzed as well. The following results are obtained: (1) From 2013 to 2018, the land subsidence, which generally showed the trend of slowing down, was mainly concentrated in the eastern, northern, and southern regions of Beijing Plain, with multiple subsidence centers. (2) Under the International Terrestrial Reference Frame 2005 (ITRF2005), the horizontal direction of all GPS points in the plain is basically the same, with the dominant movement direction being NE112.5°~NE113.8°. Under the Eurasian reference frame, the horizontal movement rate of GPS points significantly decreases. The movement rate and direction of each point are not characteristic of overall trend activity. (3) The distribution and extent of the 3D surface deformation in the Beijing Plain are controlled by the basement structure. Part of the subsided area corresponds to a Quaternary depression formed at the junction of active faults disrupting the area. Similarly, the distribution of horizontal deformation in the E-W and N-S directions of the plain is controlled by the regional basement structure comprising major faults bounding horizontal deformation. (4) Groundwater exploitation is the main cause of the 3D surface deformation in the Beijing Plain. The groundwater funnels of the second and third confined aquifer are in suitable agreement with the land subsidence. The horizontal movement in the Beijing Plain is either directed toward the center of the groundwater or the land subsidence funnel, and the deformation is directed from areas with higher to areas with lower groundwater levels.

1. Introduction

The surface deformation caused by many geophysical phenomena, to some extent, manifests the characteristics of 3D deformation, such as tectonic movements, active faults, earthquakes, volcanic eruptions, and land subsidence [1,2,3]. The accurate acquisition of the 3D deformation field information of these geophysical phenomena reveals the important role of geological data in the identification of their internal mechanism, movement process, and development trend [4,5,6]. Beijing is one of the few super-large cities in the world that use groundwater as the main source of water supply. The amount of groundwater extraction accounts for 50~70% of the total water supply of the city [7]. As a result, the groundwater resources have been over-exploited for a long time. From 1961 to 2018, the cumulative loss of groundwater resources was 9.6 billion cubic meters, which gave rise to serious surface deformation problems in the Quaternary strata of the plain [8,9]. In previous studies, most scholars concentrated on the vertical consolidation and compaction deformation of porous media caused by excessive groundwater extraction, which is land subsidence. By means of a large amount of observation data, they studied the acquisition of land subsidence monitoring information [10,11,12,13,14,15], the evolution law of land subsidence [16,17,18,19,20,21], as well as the causal mechanism of land subsidence [22,23,24,25,26,27,28,29], and up to now, they have obtained fruitful research results. However, there are few studies on the horizontal movement or deformation characteristics of the porous medium solid skeleton caused by the extraction of groundwater. The main reason is that the horizontal movement of porous media is not as intuitive and easy to monitor as the vertical consolidation and the compaction deformation [30,31]. A large amount of evidence shows that the horizontal movement and the deformation of the aquifer caused by pumping do exist, however [32,33,34,35].

Currently, Global Positioning System (GPS) and Interferometric Synthetic Aperture Radar (InSAR) methods are often used to monitor surface deformation in a large-scale area. The advantage of GPS technology lies in the high measurement accuracy and sensitivity of the ground displacement in the horizontal direction (E, N direction). With a longer observation period, the millimeter-level horizontal displacement on the ground can be detected by means of GPS technology [4]. However, the monitoring accuracy in the vertical direction is 2–3 times lower than that in the horizontal direction [6]. The advantage of InSAR technology lies in its large-scale measurement based on the image surface, which is characteristic of high spatial resolution, wide coverage, and spatial non-contact remote sensing, and that can effectively obtain the continuous surface deformation information, especially for large-scale surface deformation [36,37,38,39,40,41,42,43]. Time series InSAR analysis techniques, such as Permanent Scatters (PS) and Small Baseline Subset (SBAS) developed in recent years, have widely applied to surface deformation monitoring [43]. Due to the characteristics of side-view imaging of SAR sensors, InSAR is more sensitive to the vertical deformation of the surface, but it is difficult to detect horizontal deformations in the N-S direction. Because of this, the technology is restricted in the acquisition of 3D deformation information and the interpretation of corresponding geophysical mechanisms [6]. In previous studies, most scholars monitored the vertical deformation of the ground surface by using Envisat ASAR, Radarsat-2, Alos PALSAR, TerraSAR-X, and Sentinel-1 SAR data sets, and used InSAR technology [18,19,25,26,29], or monitored the crust horizontal movement characteristics by using GPS technology [44,45,46,47,48]. However, there are few studies on the combination of InSAR and GPS technology to obtain the 3D deformation characteristics of the surface caused by over-exploitation of groundwater. Therefore, so as to obtain complete 3D surface deformation information and further clarify the mechanism of surface deformation, it is very meaningful to combine the high-precision GPS horizontal displacement measurement with the high-density InSAR vertical deformation measurement [5,6,46].

This paper has obtained the time-series information of the 3D surface deformation field in the Beijing Plain area by using InSAR and GPS technology and analyzing its spatial distribution and evolution characteristics and laws. On this basis, the observation data of regional structure, groundwater flow field, and groundwater level are couple analyzed with the monitoring results of the surface 3D deformation field, with the relationship analyzed between different controlling factors with the 3D deformation of the surface. The research contents are mainly composed of the three following aspects. Firstly, by means of the Permanent Scatterers InSAR (PS-InSAR) and GPS technology, the vertical deformation in the plain area from 2013 to 2018 is obtained, and the evolution law is analyzed. Secondly, with the GPS observation data from 2013 to 2018 processed and calculated, the 3D coordinates of GPS points are obtained under the International Terrestrial Reference Frame 2005 (ITRF2005) and the Eurasian reference frame, respectively. Then the horizontal movement characteristics of GPS points under the two reference frames are analyzed. Thirdly, the different controlling factors are analyzed with the monitoring results of the 3D deformation, with the response relationship identified between the 3D deformation and multiple controlling factors. The results of this paper can be used to provide an important data basis and reference for further research on the causal mechanism of the 3D deformation and the construction of the water-soil 3D model in the Beijing Plain.

2. Description of the Study Area

Beijing is located on the northwestern edge of the North China Plain. The terrain is generally higher in its northwest part and lower in its southeastern part. It is surrounded by the Taihang Mountains in the west, the Yanshan Mountains in the north, and the inclined plain in the southeast. The geographical coordinates of Beijing are E115.7°~117.4°and N39.4°~41.6°, and its total area is 16410 km², the plain one of which is 6400 km² accounting for 39% of the total area [49]. The main structure of the Beijing Plain was formed during the Mesozoic Yanshan Movement and has been further transformed by the Himalayan Movement since the Cenozoic [50]. The depth of the underlying bedrock in the plain area varies greatly and is disrupted by several faults, including the NE(NNE) -SW(SSW) striking Huangzhuang-Gaoliying, Shunyi, Nanyuan-Tongxian, and Xia Dian-Mafang faults and the NW-SE direction Nankou-Sunhe fault (Figure 1). The Quaternary stratum is mainly formed by the alluvial-diluvial effect of five rivers, which are Yongding, Chaobai, Beiyun, Daqing, and Jiyun rivers [50]. On the top of the alluvial fan, the 20–40 m-thick Quaternary deposits comprise a single layer of sand or gravel or a thin layer of cohesive soil on the top of sand or gravel. In the middle and lower parts of alluvial fans and alluvial plains, the sediment thickness and the layers increase while the grain size decreases. The lithology becomes sand, gravel, and cohesive soil, alternately. Among them, the cohesive soil layer is the main stratum [24] (Figure 1).

Corresponding to the sedimentary characteristics of the Quaternary, the aquifer system in the Beijing Plain can be vertically divided into three main aquifer groups. The first aquifer system (unconfined aquifer and shallow confined aquifer) are the Holocene (Q4) and Upper Pleistocene (Q3) alluvial deposits, which reach bottom depths of about 25 m and 80~120 m, respectively. The second aquifer system (medium-deep confined aquifers) is a multi-layered structure of the Middle Pleistocene (Q2) strata, the lithology of which is mainly medium-coarse sand, partly with gravel, and the bottom depth of which is about 300 m. The third aquifer system (deep confined aquifers) is the Lower Pleistocene (Q1) stratum with a multi-layer structure, which is dominated by medium-coarse sand and gravel and is lain over the basement rocks [51] (Figure 2).

Corresponding to the division of Quaternary aquifer groups, the soil layers in the plain can be divided into three main compression layer groups. The first compression group (Q4 + Q3) is widely distributed in the Beijing Plain Area. Its lithology comprises Upper Pleistocene alluvial facies, alluvial silt, and cohesive soil, and its bottom depth is less than 100 m. The second compression group (Q2) is mainly distributed in the middle and lower parts of the Beijing Plain Area. Its lithology comprises Middle Pleistocene alluvial and alluvial lacustrine facies, silt, silty clay, and clay, and its bottom depth is less than 300 m. The third compression group (Q1) is mainly distributed in the center of several major sedimentary depressions. Its lithology comprises Lower Pleistocene alluvial and alluvial lacustrine facies, silty clay, and clay, and its roof depth is greater than 300 m [24,51] (Figure 2).

3. Data sets and Methods

3.1. Data Sets

3.1.1. SAR Data

54 RadarSAT-2 images that covered the Beijing Plain are selected as the data source for interference processing (Wide strip format). The acquisition parameters of SAR images are presented in

Table 1. Acquisition parameters of Radarsat-2 data sets.

SAR Sensor	Radarsat-2
Orbit direction	Descending
Orbit altitude	798 km
Band (wavelength)	C-band (5.6 cm)
Revisit cycle	24 days
Spatial resolution	30 m
Incidence angle (°)	27.8
Polarization	VV
Center location	40.20 116.40
Number of images	54
Temporal coverage	2013/01/16–2018/12/22

.

3.1.2. GPS and Leveling Data

The GPS observation data from 2013 to 2018 in the plain are selected for the network solution and the adjustment processing. The location of GPS points is shown in Figure 1. The GPS calculation each year contains 61 GPS and leveling integrated points (leveling points are constructed at the bottom of the GPS observation piers at the same time) and 14 Continuously Operating Reference Stations (CORS). Meanwhile, 40 primary leveling data during the same period are used to verify the accuracy of InSAR and GPS monitoring results. The geodetic elevation system is used for GPS survey, and the normal elevation system is used for leveling. There is an abnormal elevation between them. However, considering that the study area is small, when the vertical deformation (elevation difference) is compared, the elevation anomaly between them can be ignored. InSAR-derived result is the line of sight (LOS) deformation. Compared with the leveling results, the line of sight (LOS) deformation must be projected to the vertical direction.

3.1.3. Groundwater and Extensometer Data

Extensometer deformation data and corresponding groundwater-level data observed at Tianzhu subsidence monitoring station in Shunyi District were selected to analyze the relationship between layered soil deformation and groundwater level. The extensometers placed in the Tianzhu station add up to 10, which are buried in layers of different depths, respectively. The shallowest observation horizon of the extensometer is 2.4~35.5 m, and the deepest one is 219.0~238.1 m. The station boasts six groundwater-level monitoring wells altogether, which are arranged to monitor the water-level changes of aquifers at different depths. The shallowest observation horizon of monitoring wells is 27.5~31.0 m, and the deepest one is 245.0~308.1 m. At the same time, the water levels of different aquifer systems in 2017 were collected to analyze the relationship between 3D deformation.

3.2. Methods

The overall technical framework is shown in Figure 3. Firstly, by means of Permanent Scatterers Interferometry (PSI) technology, RadarSAT-2 images were processed to obtain vertical ground deformation information in the Beijing Plain. Meanwhile, the leveling data were used to validate the InSAR results. Secondly, with GAMIT/GLOBK software package, baseline vector calculation and network adjustment processing on GPS data were performed to obtain the horizontal deformation characteristics of GPS points under the ITRF2005 reference frame and the Eurasian reference frame, respectively. Then the control effect of the regional basement structure on the vertical and horizontal deformation in the plain area was studied by means of analyzing the basement structure and the 3D deformation. Finally, the impact of groundwater exploitation on the 3D surface deformation was assessed by analyzing the groundwater flow field, groundwater-level observation data, and the 3D surface deformation field.

3.2.1. Permanent Scatterers Interferometry (PSI) Technology

The basic idea of PS-InSAR is to select those pixels as the research object that maintain high coherence in the long-term sequence of SAR images, and the long-term sequence stability of them enables accurate phase information to be obtained. By decomposing the phase of each coherent target point, including the elevation, the orbit, and the phase change of the atmosphere, surface deformation can be obtained. This method can effectively overcome the influence of space and time decoherence and atmospheric delay of DInSAR and improve the monitoring accuracy of surface deformation, which can reach the millimeter level. PSI technology was first developed by Ferretti [10,37], and since then, it has largely been used to measure movements of Earth’s surface [52,53]. The main steps of the PSI processing chain are as follows: selection of a master image from the available series of SAR scenes, construction of a series of interferograms, selection of permanent scatterers (PS) of the radar signal according to amplitude dispersion index and coherence thresholds, and phase unwrapping. The differential interferometric phase Φ_PS of each PS in the interferogram can be expressed as the accumulation of five components [11]:

Φ_PS = Φ_def + Φ_top + Φ_atm + Φ_orb + Φ_n

(1)

In the formula, Φ_def is the deformation phase along the line of sight (LOS), Φ_top the topographic phase, Φ_atm the phase component due to atmospheric delay, Φ_orb the orbital error phase, and Φ_n the phase noise. The deformation phase Φ_def of PS can be separated from the other components. An external DEM with a resolution of 30 m is used to remove Φ_top from the interferometric phase. The phases due to atmosphere delay and noise are removed through filtering processes [54]. In this paper, GAMMA software is used to process 54-scene RadarSAT-2 images for time-series differential interference processing and PS point extraction. Interferometric image pairs with a baseline distance of more than 300 m are eliminated, and 30 m resolution SRTM data are used as the reference DEM. After removing the terrain phase, geocoding, and transforming the line of sight (LOS) to the vertical projection, the vertical deformation rate of the ground surface in the Beijing Plain area from 2013 to 2018 is obtained. The projection formula of the line of sight (LOS) to the vertical projection is as follows:

$$\mathrm{V}=\frac{{\mathrm{V}}_{\mathrm{LOS}}}{\cos \phi }$$

(2)

where V and ${\mathrm{V}}_{\mathrm{LOS}}$ represent the displacement rate along the vertical and LOS directions, respectively; $\phi$ is the incidence angle of the InSAR.

3.2.2. GPS Monitoring

The GPS network observation of deformation in Beijing has been carried out from August to October every year since 2008. The GPS observation level is class A. Each GPS point is observed continuously for 72 hours, and the CORS sampling interval is 30 seconds. Three bedrock reference points of BJFS, BJSH, and JIXN are selected as stable reference datums. The data calculation and adjustment processing are carried out from 3 bedrock reference points, 61 GPS and leveling integrated points, and 14 CORS stations by using the GAMIT/GLOBK software package. In this study, by using the GPS processing method mentioned above, the precise 3D coordinates of GPS points under the ITRF2005 framework in the plain area from 2013 to 2018 have been obtained. At the same time, the NNR-NUVEL1A model is used to subtract the Eurasian plate velocity from the horizontal velocity field under the ITRF frame in order to study the characteristics of the horizontal movement of GPS points much better. Finally, the horizontal velocity field of all GPS points in the Eurasian reference frame is obtained.

4. Results

4.1. Vertical Deformation

4.1.1. Regional Land Subsidence Characteristics

The time-series information of the vertical deformation of the Beijing Plain from 2013 to 2018 is obtained in this study by means of PSI and GPS technology, both of which have suitable consistency in obtaining the vertical deformation of the surface (Figure 4). Chaoyang and parts of Tongzhou, which are situated in the east of the Beijing Plain, are residential areas severely affected by land subsidence. Several major land subsidence centers are connected together, and the subsidence center velocity has been above 100 mm/y for many years. In 2014, the maximum land subsidence rate reached 143.20 mm/y. From 2015 to 2017, the maximum subsidence rate in the eastern part of the plain area exceeded 135.00 mm/y. The land subsidence in the eastern region, however, slowed down significantly in 2018, and the maximum subsidence rate declined to 118.60 mm/y. The distribution of land subsidence in the northern part of the plain was scattered, and the areas with severe subsidence were mainly distributed in the northern part of Haidian, southern Changping, and part of Shunyi. The rate of land subsidence in the northern plain began to decrease in 2016, along with the land subsidence area. This phenomenon became increasingly more obvious from 2017 to 2018. The land subsidence rate in the northern part of Haidian District was 60~90 mm/y from 2013 to 2015, and the subsidence rate fell to 30~50 mm/y in 2018. In the southern part of Changping, the land subsidence rate was 60~100 mm/y from 2013 to 2015, and the subsidence rate declined to 40~50 mm/y in 2018. In some areas of Shunyi, the land subsidence rate was 40~50 mm/y from 2013 to 2015 and decreased to 20~30 mm/y in 2018. The areas with severe land subsidence in the southern part of the plain were mainly located in Daxing District, especially the area bordering Hebei Province. The subsidence rate was always maintained at 50~70 mm/y from 2013 to 2018 (Figure 4).

From 2013 to 2018, the average subsidence rate, together with the area and volume of severe subsidence (subsidence rate greater than 50 mm/y), showed a decreasing trend in the Beijing Plain. The average land subsidence rate decreased from 21.7 mm/y in 2013 to 13.32 mm/y in 2018 and with a decrease of 8.38 mm/y (Figure 5a). The area of severe subsidence reduced from 572 km² in 2013 to 383 km² in 2018, with a decrease of 189 km² (Figure 5a). The total volume of the subsided area reduced from 13210 × 104 m³ in 2013 to 8120 × 10⁴ m³ in 2018, with a decrease of 5090 × 10⁴ m³ (Figure 5b).

4.1.2. Cumulative Land Subsidence Characteristics

The cumulative subsidence curves of six typical land subsidence centers (P1~P6) from 2013 to 2018 are obtained in this study by means of three methods, including PSI, GPS, and leveling. The locations are shown in Figure 4f. Now that the selected PS point and the GPS-leveling integrated point cannot completely overlap, the GPS-leveling integrated point is taken as the center, with the adjacent PS points selected within a radius of 100 m for analysis. As shown in the results, the cumulative land subsidence of each subsidence center is gradually increasing with time. From 2013 to 2018, the maximum accumulated subsidence at P1 reached 816.77 mm, which was the area with the largest accumulated land subsidence in Beijing. The maximum accumulated settlement at P2 was 792.10 mm, and the maximum accumulated land subsidence at P3~P6 749.00, 627.00, 456.00, and 437.00 mm, respectively. Through comparative analysis of PSI, GPS, and leveling measurement results, it is found that from 2013 to 2018, the deviation between the cumulative land subsidence obtained by PSI and the leveling measurement result was between 4 and 30 mm; the deviation between the cumulative land subsidence obtained by GPS and the leveling measurement result was between 5 and 50 mm, and the cumulative land subsidence obtained by these three methods had suitable consistency in the development trend (Figure 6).

4.1.3. Verification Accuracy of PSI and GPS Vertical Deformation

A total of 40 leveling points of 2017 are selected to validate InSAR and GPS vertical deformation results (

Table 2. Comparison of PSI, GPS vertical deformation, and leveling measurement results.

Benchmark Number	Annual Land Subsidence(mm)					Benchmark Number	Annual Land Subsidence(mm)
	PSI	GPS	Leveling	PSI- Leveling	GPS- Leveling		PSI	GPS	Leveling	PSI- Leveling	GPS- Leveling
BJ001	−42.5	−46.2	−38.4	−4.1	−7.8	BJ021	−25.6	−33.6	−18.3	−7.3	−15.3
BJ002	−38.0	−40.0	−32.0	−6.0	−8.0	BJ022	−18.3	−22.9	−15.2	−3.1	−7.7
BJ003	−39.2	−43.4	−42.9	3.7	−0.5	BJ023	2.20	−10.8	−1.1	3.3	−9.7
BJ004	−18.0	−30.0	−21.0	3.0	−9.0	BJ024	−45.0	−40.2	−52.2	7.2	12.0
BJ005	−19.8	−30.2	−25.2	5.4	−5.0	BJ025	2.10	−6.6	−1.6	3.7	−5.0
BJ006	−32.2	−25.1	−37.8	5.6	12.7	BJ026	−32.5	−40.3	−25.1	−7.4	−15.2
BJ007	−40.1	−54.2	−46.7	6.6	−7.5	BJ027	−5.4	−19.3	−9.7	4.3	−9.6
BJ008	−15.3	−19.1	−18.3	3.0	−0.8	BJ028	−24.0	−32.4	−20.0	−4.0	−12.4
BJ009	−10.2	−15.9	−6.0	−4.2	−9.9	BJ029	−29.0	−20.1	−31.4	2.4	11.3
BJ010	3.5	−1.7	−4.8	8.3	3.1	BJ030	−135.3	−124.2	−129.2	−6.1	5.0
BJ011	−134.0	−115.8	−126.9	−7.1	11.1	BJ031	−6.0	−4.0	−3.9	−2.1	−0.1
BJ012	−26.0	−37.3	−34.0	8.0	−3.3	BJ032	−16.9	−5.8	−11.5	−5.4	5.7
BJ013	−5.8	−4.9	−11.6	5.8	6.7	BJ033	2.1	−1.9	−2.9	5.0	1.0
BJ014	−19.8	−32.1	−16.3	−3.5	−15.8	BJ034	−36.2	−28.2	−33.9	−2.3	5.7
BJ015	−28.9	−35.2	−33.2	4.3	−2.0	BJ035	−33.0	−35.1	−33.1	0.1	−2.0
BJ016	−6.9	−19.2	−10.2	3.3	−9.0	BJ036	−30.0	−33.3	−32.9	2.9	−0.4
BJ017	−26.0	−35.7	−27.5	1.5	−8.2	BJ037	−12.5	−11.5	−15.1	2.6	3.6
BJ018	−6.2	−9.3	−3.5	−2.7	−5.8	BJ038	−119.4	−118.1	−118.4	−1.0	0.3
BJ019	−18.5	−12.0	−16.2	−2.3	4.2	BJ039	−57.6	−64.7	−60.5	2.9	−4.2
BJ020	−36.4	−23.6	−31.0	−5.4	7.4	BJ040	−8.3	−5.2	−9.5	1.2	4.3

). Through comparison, it is found that the InSAR and the leveling results show the same deformation characteristics. The difference range of them is between 2 and 8 mm, the mean square error is 4.33 mm, and there is a significant linear correlation, with an R² of 0.983 (Figure 7). Meanwhile, by comparing the GPS measurement with the leveling results, it is found that the difference between more than 80% of GPS points and the leveling points is less than 10 mm. Some points have a large mutual difference, the maximum of which is 15.80 mm, the mean square error is 7.56 mm, and the correlation coefficient (R²) is 0.946 (Figure 7). The verification results suggest that the accuracy of InSAR and GPS measurement is reliable, and they can reflect in detail the vertical deformation and spatial distribution characteristics of land subsidence in the Beijing Plain.

4.2. Horizontal Deformation

4.2.1. Horizontal Velocity of GPS under the ITRF2005 Reference Frame

Figure 8 shows the distribution of GPS horizontal displacement field in the Beijing Plain during 2013~2015, 2016, 2017, and 2018 under the ITRF2005 framework. The results show that the trend of all GPS points under the ITRF2005 framework was basically the same. The overall movement direction was SE, with the dominant movement direction being NE112.5°~NE113.8°, reflecting the continuous deformation feature under a unified continental dynamic environment. The movement speed of each GPS point in the N direction is −10.90~−19.73 mm/y, and the average value is −13.57 mm/y. The movement speed in the E direction was 27.12~36.19 mm/y, and the average value was 30.78 mm/y. The analysis results above are consistent with the overall crustal movement status of North China obtained by other scholars using GPS technology [45,47].

4.2.2. Accuracy Analysis of GPS Horizontal Velocity

Table 3. The horizontal velocity and accuracy of GPS in 2017 under the ITRF2005 reference frame.

Benchmark Number	Horizontal Velocity (mm/y)		Accuracy (mm/y)		Benchmark Number	Horizontal Velocity (mm/y)		Accuracy (mm/y)
	VE	VN	ΔE	ΔN		VE	VN	ΔE	ΔN
BJ001	28.82	−14.67	3.70	4.30	BJ024	30.36	−18.86	4.00	3.50
BJ002	32.62	−16.74	4.30	3.90	BJ025	29.10	−15.33	4.30	2.10
BJ003	28.29	−10.72	2.80	3.80	BJ026	31.51	−12.90	4.30	2.80
BJ004	28.69	−15.06	3.90	4.30	BJ027	34.60	−16.94	2.90	3.70
BJ005	27.59	−11.27	2.40	4.70	BJ028	27.69	−15.77	3.70	5.00
BJ006	27.29	−11.89	2.40	4.20	BJ029	32.73	−18.99	1.60	3.60
BJ007	29.48	−12.29	4.00	3.10	BJ030	30.76	−13.90	4.70	3.20
BJ008	28.40	−16.16	3.20	4.70	BJ031	34.41	−19.22	3.80	2.50
BJ009	30.92	−11.44	4.80	3.40	BJ032	32.53	−16.18	4.90	3.70
BJ010	31.08	−17.75	2.20	2.90	BJ033	31.57	−12.80	3.80	4.50
BJ011	31.76	−11.58	2.30	3.80	BJ034	34.99	−16.64	4.30	3.60
BJ012	33.22	−10.90	4.60	2.30	BJ035	33.87	−12.51	4.90	3.80
BJ013	32.65	−13.54	3.90	4.20	BJ036	31.53	−11.90	4.40	2.80
BJ014	28.19	−15.46	3.90	2.20	BJ037	31.12	−16.98	4.40	3.90
BJ015	28.14	−18.91	4.70	2.30	BJ038	33.85	−12.27	4.80	2.40
BJ016	32.39	−16.05	4.40	3.20	BJ039	28.82	−15.34	3.20	4.40
BJ017	29.20	−11.65	4.90	4.30	BJ040	32.83	−16.66	2.70	3.60
BJ018	28.36	−15.52	4.80	4.60	ZJWZ	28.58	−16.19	0.07	0.09
BJ019	28.94	−12.81	2.90	4.60	NLSH	27.80	−12.78	0.06	0.08
BJ020	32.36	−14.65	4.60	2.60	CHAO	28.92	−18.57	0.08	0.10
BJ021	32.63	−15.41	3.90	4.30	DSQI	34.12	−17.77	0.07	0.09
BJ022	31.40	−13.40	3.70	4.00	YUFA	27.32	−12.05	0.07	0.09
BJ023	29.59	−18.49	2.80	4.90	CHPN	30.10	−11.25	0.07	0.09

shows the horizontal velocity and internal accord accuracy of 40 GPS and level integration points and 6 CORS stations in 2017 under the ITRF2005 reference frame. According to the calculation results, it is shown that the error of the 40 GPS points in the velocity components of the E-W and S-N directions is basically between 2.00 and 5.00 mm/y; that of the 6 CORS stations in the velocity components of the E-W and S-N directions basically between 0.06 and 0.10 mm/y. The observation accuracy of the CORS stations is higher, which is mainly due to the better observation environment of the CORS stations and the longer time span of the observation data.

4.2.3. Horizontal Velocity of GPS under the Eurasian Reference Frame

The horizontal velocity under the ITRF2005 reference framework reflects the overall movement situation of the Beijing Plain under the global framework. Compared to global reference datums, the coordinate and velocity under this framework can be regarded as the absolute coordinate and velocity. Since the velocity under the framework contains the horizontal movement velocity of the entire plate and the horizontal strain within the region, the overall horizontal component is relatively large, which masks the relatively small deformation difference within the region. Therefore, the Eurasian plate velocity from the velocity under the ITRF frame needs to be deducted so as to obtain the internal horizontal velocity of all GPS points under the Eurasian reference frame. In the calculation process, three bedrock CORS stations, BJFS, BJSH, and JIXN, are selected as the stable reference datum, and the parameter values of the Euler vector in the Eurasian plate refer to the NNR-NUVEL1A model (longitude −112.27°, latitude 50.61°, angular velocity 0.23°/Ma).

As shown in Figure 9, the GPS points take the Eurasian plate as a reference, the horizontal velocity within the area is significantly reducing, and the inconsistency changes between the points are obvious; the movement rate and direction of each point show a certain disorder not manifesting the overall trend movement. From 2013 to 2015, 42% of GPS points moved horizontally in the SE direction, with a movement rate of 1.45~20.10 mm/y; 35% of GPS points in the NE direction, with a movement rate of 1.32~18.96 mm/y; 10% of GPS points in the SW direction, with a movement rate of 1.35~9.52 mm/y; and 13% of the GPS points in the NW direction, with a movement rate of 0.85~9.24 mm/y. The spatial distribution was roughly bounded by the northern line of the Daxing District, with a tensile deformation characteristic, and with a trend of NE to SW rotation in the northern plain area (Figure 9a). In 2016, 65% of GPS points moved horizontally in the SE direction, especially in the southern part of the plain area, with a movement rate of 1.35~18.25 mm/y. In Chaoyang, northern Tongzhou, and parts of Shunyi, the horizontal deformation was complex, some of which was indicative of tension and compression (Figure 9b). In 2017 and 2018, as shown in GPS observation results, the southern plain area mainly moved in the E and SE directions, and in Chaoyang, the northern part of Tongzhou, the southern Changping as well as part of Shunyi, the horizontal movement velocity was relatively large, without obvious directionality. The horizontal deformation, especially at the junction of several major active faults, was very complicated (Figure 9c,d). These areas often had severe land subsidence.

5. Discussion

5.1. The Control Effect of Base Structure on 3D Deformation

5.1.1. The Control Effect of Basement Structure on Vertical Deformation

By combining the land subsidence obtained by PSI with the main active faults, it is found that the distribution of land subsidence is structurally controlled in Beijing Plain. The subsided areas, part of which coincide with the Quaternary depression at the junction of several major active faults, are divided by several major active faults, which control the development of land subsidence (Figure 10a).

The NW-SE striking Nankou-Sunhe fault runs through the entire Beijing Plain, dividing the main subsidence area into two parts: the northeastern and the southwestern. The subsidence area in the northern part of Haidian and the southern part of Changping is mainly controlled by the northwestern section of the Nankou-Sunhe fault. The location of Shangzhuang’s subsided area in the north of Haidian is consistent with the Machikou Quaternary depression. The subsided area in the northern part of Changping is affected by both the NE-SW striking Huangzhuang-Gaoliying fault and the NW-SE striking Nankou-Sunhe fault. The distribution of subsided area in Shunyi is affected by both the NE-SW striking Shunyi fault and the northern section of NE-SW striking Huangzhuang-Gaoliying fault. The subsided area is mainly distributed in the Houshayu Quaternary Depression, which is bounded by Huangzhuang-Gaoliying and Shunyi fault. (Figure 10a).

The control effect of the fault on the land subsidence is mainly revealed by the offset of the fault and the different thickness of the Quaternary beds on both sides of the fault. According to the data from the Zhen 3 and Shun 3 drill holes of the upper and lower plates of the Huangzhuang-Gaoliying fault, respectively, the Zhen 3 hole in the northwestern part of the fault passes through the loose accumulation layer at 430 m to expose andesite, while the Shun 3 hole in the southeast plate only exposes the tuffaceous breccia at 645 m. Due to the influence of the early activities of the fault, the same sedimentary strata on both sides of the fault differed by more than 200 meters. (Figure 10b). According to the data from the Shun 4 and Shun 5 drill holes of the upper and lower plates of the Shunyi fault, respectively, the Shun 4 hole in the northwestern part of the fault passes through the loose accumulation layer at 474 m to expose dolomite, whereas no pre-Cenozoic bedrock has been found at a depth of 800 m in the Shun 5 hole on the southeast side of the fault. The same sedimentary strata on both sides of the fault differed by more than 200 meters. (Figure 10c). Meanwhile, the spatial distribution of land subsidence in Chaoyang and parts of Tongzhou are controlled by the southeast section of the Nankou-Sunhe fault.

It can be seen that the distribution and development trend of the main land subsidence areas in the Beijing Plain has an obvious corresponding relationship with the strike of active faults, and the structural control effect is very obvious. At present, the land subsidence in some areas has crossed the fault control and is directed parallel to the fault.

Four typical land subsidence profile lines E-E’, F-F’, G-G’, and H-H’ are drawn through the active faults in the paper, the positions of which are shown in Figure 10a. Section E-E’ passes through the northern section of the Huangzhuang-Gaoliying fault and the Shunyi fault. At the Huangzhuang-Gaoliying fault, the maximum differential land subsidence is 26.41 mm/y, and the maximum deformation gradient is 21.89 mm/km. At the Shunyi fault, the maximum differential land subsidence is 50.25 mm/y, and the maximum deformation gradient is 28.12 mm/km (Figure 11a). Section F-F’ passes through the Huangzhuang-Gaoliying fault, Nankou-Sunhe fault, Shunyi fault, and Nanyuan-Tongxian fault. At the Huangzhuang-Gaoliying fault, the maximum differential land subsidence is 28.23 mm/y, and the maximum deformation gradient is 16.50 mm/km. At the Nankou-Sunhe fault, the maximum differential land subsidence is 37.61 mm/y, and the maximum deformation gradient is 16.58 mm/km. At the Shunyi fault, the maximum differential land subsidence is 37.30 mm/y, and the maximum deformation gradient is 18.32 mm/km. At the Nanyuan-Tongxian fault, the maximum differential land subsidence is 94.44 mm/y, and the maximum deformation gradient reaches 46.24 mm/km (Figure 11b). Section G-G’ passes through the Nankou-Sunhe fault. The maximum differential land subsidence is 87.84 mm/y, and the maximum deformation gradient reaches 65.81 mm/km (Figure 11c). Section H-H’passes through the Nankou-Sunhe fault, the maximum differential land subsidence is 55.66 mm/y, and the maximum deformation gradient is 18.31 mm/km (Figure 11d). Through the analysis above, it can be known that in the areas where active faults pass through, the land subsidence profile shows obvious transitions or sudden changes, and the uneven settlement on both sides of the fault is very obvious. This is mainly due to the control effect of the basement structure.

5.1.2. The Control Effect of Basement Structure on Horizontal Deformation

The average horizontal velocity of GPS points was obtained by calculating the horizontal velocity of GPS points in the Eurasian reference frame from 2013 to 2018. The Kriging method was used to perform spatial interpolation of GPS points in the E-W direction and S-N direction to obtain continuous horizontal deformation of the surface (Figure 12). Figure 12a shows that the maximum deformation in the E direction was 14.5 mm/y, the maximum deformation in the W direction was 6.5 mm/y, and the deformation in the E direction was greater than the W direction. The distribution of the horizontal deformation in the E-W direction was controlled by the Nankou-Sunhe fault. The S-W side of the fault was dominated by E directional movement, and the N-E side was dominated by W directional movement. The surface on both sides of the fault showed a tendency to move toward each other, especially at the junction of Chaoyang and Shunyi. In the central part of Tongzhou, especially at the intersection of the Nankou-Sunhe fault and Xiadian-Mafang fault, both sides of the fault also showed the characteristics of moving toward each other. Some GPS points at the junction of Daxing and Hebei Province moved in W direction, and other areas mainly moved in E direction (Figure 12a).

Figure 12b shows that the maximum deformation in the N direction was 8.9 mm/y, the maximum deformation in the S direction was 11.7 mm/y, and the deformation in the S direction was greater than the N direction. In the northern section of the Huangzhuang-Gaoliying fault, the northern area was dominated by the S direction movement, and the southern area was dominated by the N direction movement. The surface on both sides of the fault showed a tendency to move toward each other. Near the Shunyi fault, the northern area was dominated by the N directional movement, and the southern area by the S direction movement. Near the Nanyuan-Tongxian Fault, the northern area was dominated by the S direction movement, while the southern area had a trend of the N direction movement (Figure 12b). This shows that the horizontal deformation of the ground surface also had obvious structural control characteristics. It should be noted that the GPS points were evenly distributed in the plain, and some monitoring points were far away from the fault. The horizontal deformation obtained by Kriging interpolation cannot represent the deformation characteristics of the fault itself.

5.2. Influence of Groundwater Exploitation on 3D Deformation

5.2.1. Influence of Groundwater Exploitation on Land Subsidence

According to previous studies, there is a clear correspondence between the development of land subsidence and the history of groundwater exploitation in the Beijing Plain. Excessive exploitation of groundwater is the main cause of ground subsidence in Beijing [8,9,28]. In this paper, the water level of different aquifers and the land subsidence velocity in 2017 are selected for analysis (Figure 13). The northern and western of the plain are located in the upper part of the alluvial fan, the stratum, which is mainly a single gravel layer structure and where the groundwater is mainly unconfined water, and the water table is generally at a relatively high level. Even if there is a groundwater-level funnel in these areas, it is not easy to cause land subsidence. For example, in the Miyun-Huairou-Shunyi groundwater-level funnel, the land subsidence is not obvious (Figure 13a). However, in the middle and lower parts of the plain, such as Chaoyang and Tongzhou, the stratum in these areas is an interactive structure of the cohesive soil layer and the sand layer, with the former being the main layer. When the groundwater-level funnel appears, the obvious land subsidence will result. Before the year 2000, the agricultural water exploitation layer was mainly concentrated in 50~100 m, and the industrial and domestic water exploitation layers in 100~180 m. The extraction amount of deep confined water from 180 to 300 m was very small, and the confined water with a roof depth of more than 300 m was not exploited. After 2000, due to the continuous decline of the shallow water table and the deterioration of water quality, the groundwater exploitation gradually developed to the deep layer. At the same time, groundwater exploitation has changed from unfined water or shallow confined water to medium-deep and deep confined water. Figure 13c,d illustrate that the groundwater-level funnel of the second confined aquifer (top and bottom buried depth of 100~180 m) and the third confined aquifer (top and bottom buried depth of 180~300 m) is the most consistent with the land subsidence funnel. This is consistent with the current main mining layer of groundwater in Beijing and is the main layer where land subsidence occurs.

Ten extensometers and six groundwater wells have been constructed in the Tianzhu monitoring station to respectively monitor soil deformation and water table changes in aquifers at different depths. The location of the Tianzhu station is shown in Figure 1, and the monitoring layer depth of extensometers and wells is shown in Figure 14a. It can be seen in Figure 14b–d that the confined water level of each layer showed a continuous downward trend from 2004 to 2018. During the whole year, the groundwater level was characteristic of seasonal fluctuation. With the decline of the groundwater level, the compression of each soil layer continued to increase. In contrast, with the rise of the groundwater level, the soil layer rebounded, or the amount of compression slowed down.

The extensometer F3-8 monitoring layer was 48.5~64.5 m, and the lithology was mainly fine sand and coarse sand. The average annual groundwater level dropped from −4.58 m in 2005 to −8.73 m in 2016, with a drop of 4.15 m. The extensometer F3-8 was compressed by 4.53 mm cumulatively. The groundwater level gradually began to increase in 2017, and the average annual water level rose to −7.16 m in 2018, with an increase of 1.57 m. The corresponding layered extensometer F3-8 rebounded by 0.31 mm. During a period of a year, the sand layer showed an almost synchronous compression or rebound with the decline or rise of the groundwater level. The soil layer was mainly characterized by elastic deformation (Figure 14b).

The extensometer F3-6 monitoring layer dominated by a silty clay layer was 82.3~102.0 m, containing a thin layer of fine sand with a thickness of 3.5 m. The average annual groundwater level dropped from −4.29 m in 2005 to −12.29 m in 2018, with a drop of 8.0 m. The extensometer F3-6 was compressed by 25.49 mm cumulatively. With the water level of the aquifer cyclically rising and dropping, the decline of the water level was more than its increase in each cycle, and the water lever had a trend of falling on the whole. Correspondingly, the soil layer was also subjected to loading and unloading, with the stress increasing each time. The effective stress on the soil layer increased continuously, and the soil layer was compressed continuously. It was only in 2008~2009, and in 2018 that the soil layer rebounded slightly. Moreover, the soil layer was mainly characterized by plastic deformation (Figure 14c).

The extensometer F3-4 monitoring layer, which was 117.0 m~148.5 m, was dominated by a fined sand layer containing a thin layer of silt with a thickness of 5.0 m. Before June 2018, in each cycle, the decline of the water level was more than the increase of it, with the overall trend downward. The maximum water-level drop was 15.00 m. This showed that the effective stress on the soil layer increased continuously, and the soil layer was compressed rapidly. After June 2018, the water level rose significantly, with the compression rate of the soil layer slowing down, but the compression still had a continuous increase, and the soil layer was mainly characterized by plastic deformation (Figure 14d).

5.2.2. Influence of Groundwater Exploitation on Horizontal Deformation

As shown in previous studies, the aquifer did not keep still during the process of groundwater extraction, the particle skeleton of which was dynamic and moving, but the movement velocity of which was lower than that of water. Whether the aquifer underwent vertical compaction or vertical deformation, the horizontal deformation of the aquifer always kept company with the pumping process [30,33,34,35]. With the monitoring results of 2017 taken as an example, the horizontal deformation under the Eurasian reference frame is analyzed with the land subsidence and the water level of different aquifers (Figure 15).

It can be seen in Figure 15a,b that in part of Haidian, Changping, and Shunyi in the northern plain, the horizontal deformation velocity of GPS points was 8.9~21.2 mm/y, and the direction of movement generally pointed to the center of land subsidence or groundwater funnel. It can be seen in Figure 15c,d that in Chaoyang and Tongzhou, where land subsidence developed severely, the distribution of land subsidence was more consistent with the second confined aquifer and (top and bottom buried depth of 100~180 m) and the third confined aquifer (top and bottom buried depth of 180~300 m); the horizontal movement velocity obtained by GPS points were 8.5~24.5 mm/y; the direction of movement generally pointed to the center of the land subsidence and the funnel of the groundwater, and the horizontal movement velocity gradually increased from the center of the funnel or the center of the groundwater to the outside. The main reason is that the edge of the groundwater funnel was the critical zone of stratum deformation, and there was a tensile strain zone in the direction pointing to the center of the groundwater funnel [31,33,34]. At the same time, in the western urban area located in the land subsidence funnel of Chaoyang and Tongzhou, Fengtai, and Daxing, the horizontal movement direction obtained by GPS points demonstrated the trend of moving from areas with higher groundwater levels to areas with lower groundwater level. In addition, there is also a situation in Figure 15 where the horizontal movement direction of some GPS points was inconsistent with the law mentioned above. The main reasons are as follows: Firstly, the surface deformation was caused by multiple wells pumping water at the same time, and small groundwater funnels may exist outside the groundwater funnel area in Figure 15. Secondly, the material composition of the aquifer was characteristic of heterogeneity and anisotropy, which resulted in uneven distribution of the material composition of the aquifer. Therefore, it can be further proved that the extraction of groundwater is also the main cause of the horizontal deformation of the Beijing Plain.

6. Conclusions

By combining InSAR and GPS technology, this study has obtained time-series information of the 3D surface deformation in the Beijing Plain, analyzing its spatial distribution characteristics and evolution laws. On this basis, the observation data of tectonic structure, groundwater flow field, and groundwater level have been couple analyzed with the monitoring results of the 3D deformation. Meanwhile, the relationship between different influencing factors with the 3D deformation of the surface has been analyzed as well. The results of this study can be used to provide an important data basis and reference for further research on the causal mechanism of the 3D deformation and the construction of the water-soil 3D model in the Beijing Plain. The main conclusions are as follows:

(1) Under the additional stress of the Quaternary system caused by pumping, the surface of the Beijing Plain shows significant 3D deformation characteristics. The deformation is mainly vertical, supplemented by horizontal. (2) From 2013 to 2018, the land subsidence was mainly concentrated in the eastern, northern, and southern regions of the Beijing Plain. There were multiple subsidence centers, and the land subsidence generally presented a decreasing trend. Chaoyang and parts of Tongzhou, which are situated in the east of the Beijing Plain, are residential areas severely affected by land subsidence. The subsided rate exceeded 100 mm/y for many years. The maximum settlement rate was 143.20 mm/y, the maximum accumulated settlement was 816.77 mm, and the uneven settlement was obvious.

(2) Under the ITRF2005 reference frame, the horizontal direction of all GPS points in the plain are basically the same, with the movement mainly in the SE direction and with the dominant movement direction NE112.5°~NE113.8°. The movement speed in the E direction is 27.12~36.19 mm/y and the average value 30.78 mm/y; the movement speed in the N direction is −10.90~−19.73 mm/y, and the average value −13.57 mm/y. Under the Eurasian reference frame, the horizontal movement rate of GPS points significantly decreases, with the inconsistency changes between the points obvious. The movement rate and direction of each point are not characteristic of overall trend activity. In particular, the areas with severe land subsidence at the junction of several major active faults are often the areas with high horizontal movement rates of GPS points.

(3) The distribution and the boundary of the 3D surface deformation in the Beijing Plain are obviously controlled by the basement structure. The land subsidence area is divided by several active faults, part of which coincides with the Quaternary depression at the junction of the active faults. In areas where active faults pass through, the land subsidence profile shows an obvious turn or abrupt change. The uneven land subsidence on both sides of the faults is very obvious. Similarly, the regional basement structure controls the distribution of horizontal deformation in the E-W direction and S-N direction of the plain. Several major faults become the boundary of horizontal deformation. There is an obvious difference in the direction of horizontal deformation on both sides of the faults.

(4) Groundwater exploitation is the main cause of the 3D surface deformation in the Beijing Plain. With the decline or rise of the groundwater level, the soil layer shows the characteristics of compression-rebound or slowing down of compression. The groundwater funnels of the second and third confined aquifer are in suitable agreement with the land subsidence, mainly contributing to it. The horizontal movement direction of GPS points generally moves toward the center of land subsidence or groundwater funnel or moves from areas with higher groundwater levels to areas with lower groundwater levels. This is mainly caused by the horizontal deformation component of the Quaternary aquifer system due to the extraction of groundwater.