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

: 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 signiﬁcantly 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 conﬁned 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.


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].     Figure 1).

Figure 2.
Hydrogeologic cross-section of Beijing alluvial-proluvial fan along the north to south A-A1 (the location is shown in Figure 1). 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.

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.

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.

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 Remote Sens. 2021, 13, 3964 6 of 25 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.
Remote Sens. 2021, 13, x FOR PEER REVIEW 6 of 26 eling 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.

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]: 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

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]: 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 Remote Sens. 2021, 13, 3964 7 of 25 area from 2013 to 2018 is obtained. The projection formula of the line of sight (LOS) to the vertical projection is as follows: where V and V LOS represent the displacement rate along the vertical and LOS directions, respectively; ϕ is the incidence angle of the InSAR.

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.

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 2 in 2013 to 383 km 2 in 2018, with a decrease of 189 km 2 ( Figure 5a). The total volume of the subsided area reduced from 13210 × 104 m 3

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 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 GPSleveling 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).
Remote Sens. 2021, 13, x FOR PEER REVIEW 10 of 26 cumulative land subsidence obtained by these three methods had suitable consistency in the development trend ( Figure 6).

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). 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 2 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 2 ) 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.

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). 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 2 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 2 ) 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.      Table 3 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.   Table 3 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.  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).

Accuracy Analysis of GPS Horizontal Velocity
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. sion 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.

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).  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). by more than 200 meters. (Figure 10c). Meanwhile, the spatial distribution of land sub-sidence 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. 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. 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.

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). 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.   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.

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  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 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)

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).

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 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.

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 •~N E113.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.