Surface Subsidence Analysis by Multi-Temporal InSAR and GRACE: A Case Study in Beijing

The aim of this study was to investigate the relationship between surface subsidence and groundwater changes. To investigate this relationship, we first analyzed surface subsidence. This paper presents the results of a case study of surface subsidence in Beijing from 1 August 2007 to 29 September 2010. The Multi-temporal Interferometric Synthetic Aperture Radar (multi-temporal InSAR) technique, which can simultaneously detect point-like stable reflectors (PSs) and distributed scatterers (DSs), was used to retrieve the subsidence magnitude and distribution in Beijing using 18 ENVISAT ASAR images. The multi-temporal InSAR-derived subsidence was verified by leveling at an accuracy better than 5 mm/year. Based on the verified multi-temporal InSAR results, a prominent uneven subsidence was identified in Beijing. Specifically, most of the subsidence velocities in the downtown area were within 10 mm/year, and the largest subsidence was detected in Tongzhou, with velocities exceeding 140 mm/year. Furthermore, Gravity Recovery and Climate Experiment (GRACE) data were used to derive the groundwater change series and trend. By comparison with the multi-temporal InSAR-derived subsidence results, the long-term decreasing trend between groundwater changes and surface subsidence showed a relatively high consistency, and a significant impact of groundwater changes on the surface subsidence was identified. Additionally, the spatial distribution of the subsidence funnel was partially consistent with that of groundwater depression, i.e., the former possessed a wider range than the latter. Finally, the relationship between surface subsidence and groundwater changes was determined.


Introduction
Surface subsidence is the main regional environmental geological disaster in the plain areas of China that may negatively affect the safety of infrastructure and human life [1]. In China, more than 70,000 km 2 in seventeen provinces have been affected to varying degrees by surface subsidence [2,3]. Beijing, as an international metropolis located in the northwest edge of the North China Plain, is one of the most serious areas in China suffering from surface subsidence [4]. Surface subsidence in Beijing was first reported in the 1950s, and it has developed rapidly during recent decades. By the end of 2010, the area of surface subsidence, with a cumulative subsidence of more than 50 mm in Beijing, exceeded 4200 km 2 , and the maximum cumulative subsidence in the significant subsidence area reached 1233 mm [5]. Meanwhile, the annual average subsidence rates have exceeded 100 mm/year in some areas with serious subsidence, such as in Chaoyang and Tongzhou Districts. More than in June 2015. The plain area in Beijing is mainly the piedmont alluvial-proluvial plain formed by the joint action of the Yongding River, Chaobai River, Wenyu River, Ju River and other rivers. Multiple subsidence areas (e.g., Chaoyang and Tongzhou subsidence areas) are situated in the plain area. The range of the black box illustrated in Figure 1a is the study area, covering most of the Beijing Plain; the central latitude and longitude are 39°50′ N and 116°31′ E, respectively. Figure 1b shows that the amplitude image of ENVISAT ASAR covers the study area with a total area of approximately 6956.77 km 2 (approximately 82.73 km in the range and 84.09 km in the azimuth). Since the 1960s, population growth, urban construction, and industrial and agricultural development in Beijing have caused increasing exploitation of the groundwater, resulting in a gradually expanding surface subsidence range. The exploitation of groundwater in Beijing was 5.  [40]. During 1999 to 2009, the groundwater exploitation remained larger than the groundwater recharge in Beijing. The annual average over-exploitation reached 9.55 × 10 8 m 3 /year between 1999 and 2007 [41], causing a continuous and rapid decline in the groundwater level and accelerating the surface subsidence in Beijing.

Data
In this paper, we selected eighteen level 0 descending ENVISAT ASAR images acquired from 1 August 2007 to 29 September 2010, over Beijing to estimate the line-of-sight (LOS) average surface subsidence velocity and subsidence time series. The incidence angle of the ASAR image over the study area ranges from 17.9° in the near range to 22.8° in the far range. All of the ASAR images are in VV polarization mode with a resolution of 7.804 m in slant range and 4.750 m in azimuth. The three arc-second Shuttle Radar Topography Mission (SRTM) DEM provided by the National Aeronautics and Space Administration (NASA) was used to remove the topographic phases. Meanwhile, the DORIS orbit data released by the European Space Agency (ESA) was adopted to improve the orbit data precision of the ASAR images. To validate the results, twenty-six leveling points provided by the Beijing Institute of Surveying and Mapping were used to validate the multi-temporal InSARderived subsidence results, the leveling data were acquired from 2007 to 2010.
As mentioned in Section 2.1, large groundwater loss has occurred in Beijing, and the groundwater changes are mainly concentrated on the plain area of Beijing. Although the study area  [40]. During 1999 to 2009, the groundwater exploitation remained larger than the groundwater recharge in Beijing. The annual average over-exploitation reached 9.55 × 10 8 m 3 /year between 1999 and 2007 [41], causing a continuous and rapid decline in the groundwater level and accelerating the surface subsidence in Beijing.

Data
In this paper, we selected eighteen level 0 descending ENVISAT ASAR images acquired from 1 August 2007 to 29 September 2010, over Beijing to estimate the line-of-sight (LOS) average surface subsidence velocity and subsidence time series. The incidence angle of the ASAR image over the study area ranges from 17.9 • in the near range to 22 orbit data released by the European Space Agency (ESA) was adopted to improve the orbit data precision of the ASAR images. To validate the results, twenty-six leveling points provided by the Beijing Institute of Surveying and Mapping were used to validate the multi-temporal InSAR-derived subsidence results, the leveling data were acquired from 2007 to 2010.
As mentioned in Section 2.1, large groundwater loss has occurred in Beijing, and the groundwater changes are mainly concentrated on the plain area of Beijing. Although the study area is much smaller than the spatial resolution of GRACE global harmonic solutions, the true amplitudes of the large groundwater changes in the study area can be retrieved by GRACE solutions with the help of hydrologic modeling. In this study, we adopted the monthly gravity fields from the Release (RL) 05 solution provided by the Center for Space Research (CSR) at the University of Texas at Austin from August 2007 to September 2010 to estimate the regional TWS changes by averaging the values in grids over the Beijing area. TWS changes at 6 grid points derived from GRACE in the region from 115 • -118 • E and 39 • -41 • N with an area of~56,000 km 2 were selected. Soil moisture (SM) changes from the NOAH version of the Global Land Data Assimilation System (GLDAS) [42] at a spatial resolution of 1 • × 1 • were used to obtain regional groundwater storage variations by subtracting them from the GRACE-based TWS changes.

Multi-Temporal InSAR Technique
The multi-temporal InSAR technique combines the approach of selecting high-coherence points in the PS InSAR technique and the superiorities of SBAS InSAR that retrieve high-accuracy deformation information based on multiple-master interferograms and perform better in detecting distributed scatterers. This approach first uses the PS InSAR and SBAS InSAR techniques to select the high-coherence points from images and then retrieves the deformation information via joint processing of the high-coherence point data sets. Therefore the spatial sampling of deformation information is improved. In this study, the subsidence velocity map and time series in the Beijing area were obtained by processing the 18 ENVISAT ASAR images with the multi-temporal InSAR technique. The main steps are as follows ( Figure 2): (1) The image acquired on 14 October 2009, was selected as the super master image for the interferometric combinations, and all the slave images were co-registered and resampled to the super master image. (2) Seventeen interferograms were generated for PS analysis ( Figure 3a); meanwhile, 71 pairs of small baseline differential interferograms were processed using SBAS InSAR with a perpendicular baseline shorter than 800 m and a temporal baseline less than 400 days ( Figure 3b). Supposing interferogram j is generated from a master image and a slave image acquired at times t B and t A (t B > t A ), after removing the flat earth effect and topographic phase, the interferometric phase of a pixel located at coordinates (x,r) in interferogram j can be expressed as follows [22,43]: where φ(t B ,x,r) and φ(t A ,x,r) are the phase values of the SAR images at times t B and t A , respectively. φ def,j (x,r) represents the deformation phase between times t A and t B in the LOS direction. φ atm,j (x,r) refers to the atmospheric phase error. φ orb,j (x,r) is the residual phase due to orbit inaccuracies.
φ θ,j (x,r) is the residual phase due to look angle error, and φ noise,j (x,r) is the random noise phase. (3) PSs were identified from the interferograms generated by PS InSAR using the algorithm proposed by Hooper et al. [43]. We retrieved DSs from small baseline interferograms in the same way that PS pixels were retrieved, i.e., following the algorithm of Hooper et al. [43]. The spatially-uncorrelated look angle error term, which includes contributions from both spatially-uncorrelated height errors and deviations of the pixel's phase center from its physical center, was estimated and removed after selecting the PSs and DSs. After calculating the equivalent small baseline interferometric phase for the PSs by recombining the single-master interferometric phase, the small baseline interferometric phases from both selected PSs and DSs were combined. If a pixel occurred in both data sets, a weighted mean value for the phase was calculated by summing the complex signal from both data sets [28]. (4) The phase was unwrapped by a minimum-cost flow algorithm. Then, the spatially-correlated look angle error and other spatially-correlated noise were estimated and subtracted from the differential phase. Next, the deformation rates of the study area were obtained by a least-squares algorithm because no isolated interferogram clusters in the analysis. After obtaining the deformation rates, according to the time span between ASAR images, the corresponding deformation time series was derived. obtained by a least-squares algorithm because no isolated interferogram clusters in the analysis. After obtaining the deformation rates, according to the time span between ASAR images, the corresponding deformation time series was derived.

Calculation of TWS and Groundwater Storage Changes from GRACE
Based on the method of Wahr et al. [44], TWS changes can be recovered from the GRACE temporal gravity field, and the calculation model of TWS changes can be expressed by the equivalent water height (EWH) as follows [44][45][46]: where a (63,781,363.3 m) and ave (5517 kg/m 3 ) are the average radius and density of the Earth, respectively; water is the density of water (1000 kg/m 3 ); is the fully normalized associated Legendre function; kl represents the degree-l Love numbers; ΔClm and ΔSlm are the variations in the spherical harmonic coefficients of the gravity field with respect to the average gravity field for the period ranging from August 2007 to September 2010; and θ and λ are the co-latitude and longitude, respectively.
Monthly gravity fields from RL05 (at degree-l and order-m 60) provided by the CSR for the obtained by a least-squares algorithm because no isolated interferogram clusters in the analysis. After obtaining the deformation rates, according to the time span between ASAR images, the corresponding deformation time series was derived.

Calculation of TWS and Groundwater Storage Changes from GRACE
Based on the method of Wahr et al. [44], TWS changes can be recovered from the GRACE temporal gravity field, and the calculation model of TWS changes can be expressed by the equivalent water height (EWH) as follows [44][45][46]: where a (63,781,363.3 m) and ave (5517 kg/m 3 ) are the average radius and density of the Earth, respectively; water is the density of water (1000 kg/m 3 ); is the fully normalized associated Legendre function; kl represents the degree-l Love numbers; ΔClm and ΔSlm are the variations in the spherical harmonic coefficients of the gravity field with respect to the average gravity field for the period ranging from August 2007 to September 2010; and θ and λ are the co-latitude and longitude, respectively.

Calculation of TWS and Groundwater Storage Changes from GRACE
Based on the method of Wahr et al. [44], TWS changes can be recovered from the GRACE temporal gravity field, and the calculation model of TWS changes can be expressed by the equivalent water height (EWH) as follows [44][45][46]: where a (63,781,363.3 m) and ρ ave (5517 kg/m 3 ) are the average radius and density of the Earth, respectively; ρ water is the density of water (1000 kg/m 3 ); P lm is the fully normalized associated Legendre function; k l represents the degree-l Love numbers; ∆C lm and ∆S lm are the variations in the spherical harmonic coefficients of the gravity field with respect to the average gravity field for the period ranging from August 2007 to September 2010; and θ and λ are the co-latitude and longitude, respectively. Monthly gravity fields from RL05 (at degree-l and order-m 60) provided by the CSR for the period of August 2007 to September 2010 were used to calculate the TWS changes. The data processing, which followed the methods of Luo et al. [47], included the following steps: (1) In GRACE observed harmonics coefficients, the degree-2 zonal harmonics were not well estimated because the GRACE orbit geometry is less sensitive to the coefficient of the gravity field. Therefore, degree-2 zonal C20 time series were replaced by analyzed Satellite Laser Ranging (SLR) data [48]; (2) Due to the correlated errors among certain spherical harmonics coefficients and other errors found in GRACE spherical harmonic coefficients, a hybrid filtering scheme that combined de-correlation filter P3M6 (at spherical harmonics orders 6 and above, a degree 3 polynomial is fitted via least squares and is removed from even and odd coefficients pairs) [49] and a 300 km Fan filter [50] were applied to reduce noise in the GRACE data. The computed model can be expressed as follows [51,52]: where W l and W m are the smoothing kernel functions associated with order and degree, respectively. That is: ; r is the filtering radius; and ∆ C lm and ∆ S lm are the variations in the spherical harmonic coefficients of the gravity field in Equation (2) after decorrelation filtering.
The GLDAS-based SM was initially processed in the same way as the GRACE data (i.e., model grids were transferred to spherical harmonic coefficients, truncated at degree and order 60, and a de-correlation filter and a 300 km Fan filter were applied) for consistency with the GRACE resolution. Then, groundwater storage changes were computed by subtracting SM changes from GRACE TWS changes [53].

Results
The LOS average subsidence velocity map in the study area derived by the multi-temporal InSAR technique is shown in Figure 4a, where the negative values indicate that the surface is moving away from the satellite (i.e., subsidence in the LOS direction) and the positive values indicate surface uplift. The standard deviations of the average subsidence velocities in the study area are shown in Figure 4b. The standard deviations were obtained by computing the deviations of the linear fitting of the velocities. If a PS/DS point showed a strong non-linear motion, it resulted in a large residual with respect to the linear model, i.e., in a high standard deviation. Based on the collected leveling data, a stable point called R (Figure 4a) was selected as a reference point in the study area; all of the monitoring points' average subsidence velocities in the deformation velocity field are related to the reference point.
As shown in Figure 4a, in the study area, uneven subsidence is prominent. The five significant subsidence areas are located in Changping, Shunyi, Chaoyang, Tongzhou and Daxing Districts. Additionally, a pronounced subsidence area located in Langfang City can be identified, which is adjacent to Beijing and has a maximum velocity exceeding 75 mm/year. The main subsidence areas in Beijing are distributed in the Chaobai River, Wenyu River and Yongding River alluvial-proluvial fan plains. The Changping, Shunyi, Chaoyang and Tongzhou subsidence areas are mainly distributed along the Wenyu River, Chaobai River, Qing River, and Tonghui River basins, and these subsidence areas gradually connect into a continuous area and constitute the "North" subsidence area. The "North" subsidence area gradually expands toward the northern part of Shunyi, the eastern part of Chaoyang, and the eastern and southern parts of Tongzhou.  In addition, a new subsidence area has formed east of the study area, i.e., the Yanjiao Town subsidence area, and some of the subsidence velocities in this area have reached 50 mm/year. The Daxing subsidence area is mainly distributed in the Yongding River coast and constitutes the "South" subsidence area. The subsidence range of the "South" subsidence area is expanding markedly in a northward direction.
The surface subsidence in downtown Beijing (e.g., Haidian, Fengtai and Shijingshan Districts and the central urban area) is small and relatively stable, and most of the subsidence velocities range from −10 mm/year to 10 mm/year. Compared with the Changping, Chaoyang and Tongzhou subsidence areas, subsidence velocities in the Shunyi and Daxing subsidence areas are lower. The most serious subsidence area in Daxing District is Lixian Town, where the maximum subsidence velocity exceeds 70 mm/year. The maximum subsidence velocity in the Shunyi subsidence area reaches 60 mm/year. Multiple subsidence funnels have formed in Changping District, among which the Shahe and Shangzhuang Town subsidence centers are more serious, with an annual subsidence velocity exceeding 80 mm/year. The subsidence in Laiguangyin, Sunhe and Jinzhan Townships in Chaoyang District are substantial; the maximum subsidence velocity in the Sunhe subsidence center exceeds 80 mm/year, and the most serious subsidence area is the Jinzhan subsidence area in which the maximum subsidence velocity exceeds 130 mm/year. Guanzhuan, Sanjianfang and Dougezhuang Townships are the most severe subsidence areas in Tongzhou District. Subway Batong Line, Subway Line 6 and Huitong River run through these areas along an east-west direction, and these areas include several highways and railways. In addition, the groundwater exploitation in these areas is serious. Therefore, it is expected that the surface subsidence in Tongzhou District is not only affected by the surface dynamic load and the utilization of underground space but also by the groundwater over-exploitation and geological structure. Multiple subsidence centers are located in the area, which includes the metro, highways, railway and rivers, and the maximum annual average subsidence velocities range between 140 mm/year and 150 mm/year. Figure 4b shows that the standard deviations of the average velocities from all monitoring points derived by the multi-temporal InSAR technique range from 0.2 mm/year to 8.5 mm/year; the monitoring points with larger standard deviations are mainly distributed in the serious subsidence areas (e.g., Tongzhou, Chaoyang and Changping Districts); the standard deviations of the average velocities of most monitoring points distributed in the downtown area of Beijing are less than 2 mm/year. The statistical analysis of the standard deviations of all monitoring points ( Figure 5) suggests that the standard deviations of the average velocities of 96.81% monitoring points are less than 5 mm/year and 88.46% are less than 3 mm/year; the standard deviation of the subsidence results In addition, a new subsidence area has formed east of the study area, i.e., the Yanjiao Town subsidence area, and some of the subsidence velocities in this area have reached 50 mm/year. The Daxing subsidence area is mainly distributed in the Yongding River coast and constitutes the "South" subsidence area. The subsidence range of the "South" subsidence area is expanding markedly in a northward direction.
The surface subsidence in downtown Beijing (e.g., Haidian, Fengtai and Shijingshan Districts and the central urban area) is small and relatively stable, and most of the subsidence velocities range from −10 mm/year to 10 mm/year. Compared with the Changping, Chaoyang and Tongzhou subsidence areas, subsidence velocities in the Shunyi and Daxing subsidence areas are lower. The most serious subsidence area in Daxing District is Lixian Town, where the maximum subsidence velocity exceeds 70 mm/year. The maximum subsidence velocity in the Shunyi subsidence area reaches 60 mm/year. Multiple subsidence funnels have formed in Changping District, among which the Shahe and Shangzhuang Town subsidence centers are more serious, with an annual subsidence velocity exceeding 80 mm/year. The subsidence in Laiguangyin, Sunhe and Jinzhan Townships in Chaoyang District are substantial; the maximum subsidence velocity in the Sunhe subsidence center exceeds 80 mm/year, and the most serious subsidence area is the Jinzhan subsidence area in which the maximum subsidence velocity exceeds 130 mm/year. Guanzhuan, Sanjianfang and Dougezhuang Townships are the most severe subsidence areas in Tongzhou District. Subway Batong Line, Subway Line 6 and Huitong River run through these areas along an east-west direction, and these areas include several highways and railways. In addition, the groundwater exploitation in these areas is serious. Therefore, it is expected that the surface subsidence in Tongzhou District is not only affected by the surface dynamic load and the utilization of underground space but also by the groundwater over-exploitation and geological structure. Multiple subsidence centers are located in the area, which includes the metro, highways, railway and rivers, and the maximum annual average subsidence velocities range between 140 mm/year and 150 mm/year. Figure 4b shows that the standard deviations of the average velocities from all monitoring points derived by the multi-temporal InSAR technique range from 0.2 mm/year to 8.5 mm/year; the monitoring points with larger standard deviations are mainly distributed in the serious subsidence areas (e.g., Tongzhou, Chaoyang and Changping Districts); the standard deviations of the average velocities of most monitoring points distributed in the downtown area of Beijing are less than 2 mm/year. The statistical analysis of the standard deviations of all monitoring points ( Figure 5) suggests that the standard deviations of the average velocities of 96.81% monitoring points are less

Validation with Leveling
To verify the reliability of the multi-temporal InSAR-derived results, we compared and analyzed the multi-temporal InSAR-derived results with the leveling data provided by the Beijing Institute of Surveying and Mapping; the location distribution of all leveling points are illustrated in Figure 6. Several steps were performed for the validation: (1) the overlapping period data of the leveling and ENVISAT ASAR image data were selected; (2) the vertical deformation of all leveling points were projected into the LOS direction; and (3) the initial subsidence on 1 August 2007, was assumed to be zero. Because there are few opportunities for a leveling point and the corresponding PS point to be located in the same place, two different methods were adopted to extract the subsidence of the PS point corresponding to the nearest leveling point. In the first method (Kriging method), the spherical variogram model (detailed information shown in Appendix A) was computed based on the multitemporal InSAR-derived PS point subsidence. Based on the spherical variogram model, Kriging was employed to estimate the subsidence of the corresponding PS point located in the same location with respect to the leveling point. For the second method (average method), the average subsidence rates of the PS points that lie within 200 m of the leveling point were used as the subsidence of the corresponding PS point. We adopted the Kriging and average methods to extract the average subsidence velocity of a PS point located at the same location as a leveling point; the results are shown in Figure 7. Figure 7a shows that the deformation changes of the three curves are approximately the same, indicating that the average subsidence velocities extracted by the two methods agree well with the results obtained from the leveling data. Moreover, the differences between multi-temporal InSAR and leveling are within ±10 mm/year, as shown in Figure 7b. We quantitatively compared the average subsidence velocities of the 26 points (i.e., BM1-BM26) derived using the multi-temporal InSAR technique and leveling (results listed in Table 1). For the Kriging method, the mean error is 1.44 mm/year, the rootmean-square error (RMSE) is 4.97 mm/year, the maximum difference is 8.33 mm/year and the minimum difference is 1.01 mm/year, while the mean error, RMSE, MAX and MIN achieved from the average method are 1.90 mm/year, 4.78 mm/year, 8.60 mm/year, and 0.18 mm/year, respectively. The accuracies of the two methods are similar. Therefore, the validation results indicate that the average subsidence velocities obtained by the multi-temporal InSAR technique are highly consistent with those derived by leveling, and the multi-temporal InSAR technique can satisfy the subsidence monitoring requirements in the study area.

Validation with Leveling
To verify the reliability of the multi-temporal InSAR-derived results, we compared and analyzed the multi-temporal InSAR-derived results with the leveling data provided by the Beijing Institute of Surveying and Mapping; the location distribution of all leveling points are illustrated in Figure 6. Several steps were performed for the validation: (1) the overlapping period data of the leveling and ENVISAT ASAR image data were selected; (2) the vertical deformation of all leveling points were projected into the LOS direction; and (3) the initial subsidence on 1 August 2007, was assumed to be zero. Because there are few opportunities for a leveling point and the corresponding PS point to be located in the same place, two different methods were adopted to extract the subsidence of the PS point corresponding to the nearest leveling point. In the first method (Kriging method), the spherical variogram model (detailed information shown in Appendix A) was computed based on the multi-temporal InSAR-derived PS point subsidence. Based on the spherical variogram model, Kriging was employed to estimate the subsidence of the corresponding PS point located in the same location with respect to the leveling point. For the second method (average method), the average subsidence rates of the PS points that lie within 200 m of the leveling point were used as the subsidence of the corresponding PS point. We adopted the Kriging and average methods to extract the average subsidence velocity of a PS point located at the same location as a leveling point; the results are shown in Figure 7. Figure 7a shows that the deformation changes of the three curves are approximately the same, indicating that the average subsidence velocities extracted by the two methods agree well with the results obtained from the leveling data. Moreover, the differences between multi-temporal InSAR and leveling are within ±10 mm/year, as shown in Figure 7b. We quantitatively compared the average subsidence velocities of the 26 points (i.e., BM1-BM26) derived using the multi-temporal InSAR technique and leveling (results listed in Table 1). For the Kriging method, the mean error is 1.44 mm/year, the root-mean-square error (RMSE) is 4.97 mm/year, the maximum difference is 8.33 mm/year and the minimum difference is 1.01 mm/year, while the mean error, RMSE, MAX and MIN achieved from the average method are 1.90 mm/year, 4.78 mm/year, 8.60 mm/year, and 0.18 mm/year, respectively. The accuracies of the two methods are similar. Therefore, the validation results indicate that the average subsidence velocities obtained by the multi-temporal InSAR technique are highly consistent with those derived by leveling, and the multi-temporal InSAR technique can satisfy the subsidence monitoring requirements in the study area.   To analyze the temporal evolution of surface subsidence in the study area, we selected eight leveling points (i.e., BM3, BM4, BM8, BM10, BM12, BM14, BM18 and BM25 illustrated in Figure 6) and the corresponding PS points to compare and analyze their subsidence time series. The results are shown in Figure 8.
As shown in Figure 8,   To analyze the temporal evolution of surface subsidence in the study area, we selected eight leveling points (i.e., BM3, BM4, BM8, BM10, BM12, BM14, BM18 and BM25 illustrated in Figure 6) and the corresponding PS points to compare and analyze their subsidence time series. The results are shown in Figure 8.
As shown in Figure 8,  The leveling data were acquired only once during the summers of 2007, 2008, 2009 and 2010. To analyze the temporal evolution of surface subsidence in the study area, we selected eight leveling points (i.e., BM3, BM4, BM8, BM10, BM12, BM14, BM18 and BM25 illustrated in Figure 6) and the corresponding PS points to compare and analyze their subsidence time series. The results are shown in Figure 8. Figure 8, the multi-temporal InSAR-derived results show nonlinear subsidence. The selected points show different amounts of subsidence depending on their location. Leveling point BM3 and the corresponding PS point are located in the downtown area of Beijing (Figure 6), where the accumulated subsidence values of the two different measuring methods are less than 30 mm (Figure 8a). Meanwhile, the subsidence is nearly steady after April 2009, which is related to the stringent restrictions on groundwater exploitation in the urban area and the small compressibility of most surface sediments (e.g., coarse particle sand, gravel, and cobble) in the western part of Beijing. Among the selected points, BM10 and BM18 are located in the most serious subsidence area, the Tongzhou District ( Figure 6). The subsidence of these pairs continues to increase, and the accumulated subsidence has exceeded 200 mm at both locations. At BM10 and BM18, the results derived by multi-temporal InSAR and leveling coincide well, as shown in Figure 8d,g. For other points, the multi-temporal InSAR-derived results are also consistent with the leveling-derived results. The comparative analysis indicates that the subsidence values acquired by leveling agree well with the multi-temporal InSAR-derived results.  (Figure 8a). Meanwhile, the subsidence is nearly steady after April 2009, which is related to the stringent restrictions on groundwater exploitation in the urban area and the small compressibility of most surface sediments (e.g., coarse particle sand, gravel, and cobble) in the western part of Beijing. Among the selected points, BM10 and BM18 are located in the most serious subsidence area, the Tongzhou District ( Figure 6). The subsidence of these pairs continues to increase, and the accumulated subsidence has exceeded 200 mm at both locations. At BM10 and BM18, the results derived by multi-temporal InSAR and leveling coincide well, as shown in Figure 8d,g. For other points, the multi-temporal InSAR-derived results are also consistent with the leveling-derived results. The comparative analysis indicates that the subsidence values acquired by leveling agree well with the multi-temporal InSAR-derived results.  Figure 6 between the subsidence time series of multi-temporal InSAR (blue squares) and leveling (purple triangles) data.

Comparison between Subsidence Changes and Groundwater Changes
To analyze the correlation between surface subsidence and groundwater changes in Beijing, the monthly gravity fields from RL05 (at degree-l and order-m 60) between August 2007 and September 2010 were first used to estimate the TWS changes in the study area. The long-term trend in TWS changes estimated from GRACE is shown in Figure 9. The mean long-term trend of TWS changes in the study area is −8.0 mm/year (in terms of EWH). The time series of the regional TWS changes in the study area were achieved by averaging the values in the study area. This result is illustrated in Figure  10a (blue curve) and shows pronounced seasonal variations.
The TWS changes estimated from the GRACE data include primarily groundwater changes, surface water changes, soil moisture (SM) changes, and ice and snow water changes. Beijing is located in a warm temperate and semi-humid continental monsoon climate zone. Thus, the snow water component is typically small, and its contribution to the long-term TWS trend in the Beijing area is negligible. Nearly all of the major rivers in Beijing are dammed for municipal and industrial use; no reservoirs are located in the study area. The surface water reservoir storage changes in the study area have little influence on long-term TWS changes during the study period [54]. Because our aim is to derive the long-term trend in groundwater changes, the impact of the surface water component on TWS changes is beyond our consideration.  Figure 6 between the subsidence time series of multi-temporal InSAR (blue squares) and leveling (purple triangles) data.

Comparison between Subsidence Changes and Groundwater Changes
To analyze the correlation between surface subsidence and groundwater changes in Beijing, the monthly gravity fields from RL05 (at degree-l and order-m 60) between August 2007 and September 2010 were first used to estimate the TWS changes in the study area. The long-term trend in TWS changes estimated from GRACE is shown in Figure 9. The mean long-term trend of TWS changes in the study area is −8.0 mm/year (in terms of EWH). The time series of the regional TWS changes in the study area were achieved by averaging the values in the study area. This result is illustrated in Figure 10a (blue curve) and shows pronounced seasonal variations.
The TWS changes estimated from the GRACE data include primarily groundwater changes, surface water changes, soil moisture (SM) changes, and ice and snow water changes. Beijing is located in a warm temperate and semi-humid continental monsoon climate zone. Thus, the snow water component is typically small, and its contribution to the long-term TWS trend in the Beijing area is negligible. Nearly all of the major rivers in Beijing are dammed for municipal and industrial use; no reservoirs are located in the study area. The surface water reservoir storage changes in the study area have little influence on long-term TWS changes during the study period [54]. Because our aim is to derive the long-term trend in groundwater changes, the impact of the surface water component on TWS changes is beyond our consideration.
The Yan Mountains and Xi Mountain are located in the northern and western parts of Beijing (Figure 9), respectively, and the area of these mountains accounts for approximately 60 percent of the area of Beijing. In Figure 4a, it was shown that the subsidence velocities are significantly larger than the surface uplift. Additionally, because the mountains have limited capacity to store groundwater and groundwater exploitation in Beijing is mainly concentrated in the plain area, the TWS change estimated from the GRACE data in the selected area results largely from groundwater depletion in the plain region of Beijing. Although the area for the multi-temporal InSAR data is much smaller than the footprint of GRACE data, the prior conclusion means that it is reasonable to compare the GRACE-based groundwater trend with the multi-temporal InSAR-derived surface subsidence trend. The Yan Mountains and Xi Mountain are located in the northern and western parts of Beijing (Figure 9), respectively, and the area of these mountains accounts for approximately 60 percent of the area of Beijing. In Figure 4a, it was shown that the subsidence velocities are significantly larger than the surface uplift. Additionally, because the mountains have limited capacity to store groundwater and groundwater exploitation in Beijing is mainly concentrated in the plain area, the TWS change estimated from the GRACE data in the selected area results largely from groundwater depletion in the plain region of Beijing. Although the area for the multi-temporal InSAR data is much smaller than the footprint of GRACE data, the prior conclusion means that it is reasonable to compare the GRACEbased groundwater trend with the multi-temporal InSAR-derived surface subsidence trend.  Figure  10a (purple curve). After removing the SM changes, the TWS changes in the Beijing area were found to be dominated by groundwater, which agrees with the intensive groundwater consumption reported in previous studies [55,56]. Therefore, groundwater changes were estimated as the difference between the GRACE-based TWS changes and the GLDAS-based SM changes, which decreased gradually during the study period (black curve in Figure 10b). The red line in Figure 10b, which was derived by the linear fitting method, reflects the long-term trend in groundwater changes in the study area. We took an average of all multi-temporal InSAR pixels to obtain the average subsidence time series for the study area, with the result (blue curve) shown in Figure 10c. The average subsidence presents a nonlinear declining tendency with seasonal variability. We adopted the linear fitting method to deduct the seasonal variability in the average subsidence and obtained the relatively long-term trend of surface subsidence in the study area, as shown in Figure 10c (red dotted line). Figure 10d shows that both the groundwater changes and subsidence changes in the study area exhibit a continuously decreasing tendency from August 2007 to September 2010. In addition, the GRACE-based groundwater decreases at a rate of 5.3 mm/year (in terms of EWH), while the multitemporal InSAR-derived surface subsidence trend is 3.9 mm/year (in terms of subsidence). The longterm decreasing trend in both the groundwater changes and average subsidence show a relatively high consistency. A comparative analysis shows that the groundwater decline and surface subsidence in the Beijing area are consistent. The groundwater over-exploitation seriously affects the stability of the surface in Beijing.  Figure 10a (purple curve). After removing the SM changes, the TWS changes in the Beijing area were found to be dominated by groundwater, which agrees with the intensive groundwater consumption reported in previous studies [55,56]. Therefore, groundwater changes were estimated as the difference between the GRACE-based TWS changes and the GLDAS-based SM changes, which decreased gradually during the study period (black curve in Figure 10b). The red line in Figure 10b, which was derived by the linear fitting method, reflects the long-term trend in groundwater changes in the study area. We took an average of all multi-temporal InSAR pixels to obtain the average subsidence time series for the study area, with the result (blue curve) shown in Figure 10c. The average subsidence presents a nonlinear declining tendency with seasonal variability. We adopted the linear fitting method to deduct the seasonal variability in the average subsidence and obtained the relatively long-term trend of surface subsidence in the study area, as shown in Figure 10c (red dotted line). Figure 10d shows that both the groundwater changes and subsidence changes in the study area exhibit a continuously decreasing tendency from August 2007 to September 2010. In addition, the GRACE-based groundwater decreases at a rate of 5.3 mm/year (in terms of EWH), while the multi-temporal InSAR-derived surface subsidence trend is 3.9 mm/year (in terms of subsidence). The long-term decreasing trend in both the groundwater changes and average subsidence show a relatively high consistency. A comparative analysis shows that the groundwater decline and surface subsidence in the Beijing area are consistent. The groundwater over-exploitation seriously affects the stability of the surface in Beijing.

The Correlation between Surface Subsidence and Groundwater Depression
Groundwater is an important factor for maintaining soil stress balance. Groundwater overexploitation may destroy the original soil stress balance and lead to surface subsidence. To further analyze the correlation between the surface subsidence and groundwater exploitation in the Beijing area, we compared and analyzed the spatial distribution of the average subsidence velocity from August 2007 to September 2010 and the groundwater depression cone of the main mining layer in 2010 in the Beijing area; the result is shown in Figure 11. The groundwater depression is mainly distributed in the eastern and northeastern parts of downtown Beijing, and the more serious subsidence areas (i.e., the Shunyi, Chaoyang and Tongzhou subsidence areas) are primarily located in the groundwater depression cone. The subsidence funnel and the groundwater depression are consistent to some extent, although the distribution range of the subsidence funnel is larger than that of the groundwater depression. The surface subsidence in Beijing is mainly from the deformation of the compressible deposits. Under the same stress condition, areas with higher compressible deposits are more prone to surface subsidence. The compressible sediments in the western part of Beijing (e.g., downtown) are mainly coarse-particle sand, gravel, and cobble with small compressibility, while the aquifers in the eastern, northern and southern regions (e.g., Changping, Chaoyang, Tongzhou and Daxing) are mainly composed of silt clays and fine sand characterized by thick compressible sediments [57]. The uneven surface subsidence is strongly influenced by the thickness variability of the compressible sediments, and the serious surface subsidence is generally distributed over thick compressible sediments. Thus, although the formation and development of surface subsidence in the Beijing area are mainly caused by groundwater over-exploitation, the formation and development of surface subsidence are also affected by the thickness of the compressible deposits, the hydrogeological condition, stratum lithology, and shallow surface space utilization.

The Correlation between Surface Subsidence and Groundwater Depression
Groundwater is an important factor for maintaining soil stress balance. Groundwater over-exploitation may destroy the original soil stress balance and lead to surface subsidence. To further analyze the correlation between the surface subsidence and groundwater exploitation in the Beijing area, we compared and analyzed the spatial distribution of the average subsidence velocity from August 2007 to September 2010 and the groundwater depression cone of the main mining layer in 2010 in the Beijing area; the result is shown in Figure 11. The groundwater depression is mainly distributed in the eastern and northeastern parts of downtown Beijing, and the more serious subsidence areas (i.e., the Shunyi, Chaoyang and Tongzhou subsidence areas) are primarily located in the groundwater depression cone. The subsidence funnel and the groundwater depression are consistent to some extent, although the distribution range of the subsidence funnel is larger than that of the groundwater depression. The surface subsidence in Beijing is mainly from the deformation of the compressible deposits. Under the same stress condition, areas with higher compressible deposits are more prone to surface subsidence. The compressible sediments in the western part of Beijing (e.g., downtown) are mainly coarse-particle sand, gravel, and cobble with small compressibility, while the aquifers in the eastern, northern and southern regions (e.g., Changping, Chaoyang, Tongzhou and Daxing) are mainly composed of silt clays and fine sand characterized by thick compressible sediments [57]. The uneven surface subsidence is strongly influenced by the thickness variability of the compressible sediments, and the serious surface subsidence is generally distributed over thick compressible sediments. Thus, although the formation and development of surface subsidence in the Beijing area are mainly caused by groundwater over-exploitation, the formation and development of surface subsidence are also affected by the thickness of the compressible deposits, the hydrogeological condition, stratum lithology, and shallow surface space utilization. Figure 11. Distribution of the main groundwater depression (black curve) [58] and the LOS average subsidence velocity in the study area.

Conclusions
In this study, we adopted the multi-temporal InSAR technique to derive subsidence information for Beijing from 1 August 2007 to 29 September 2010. The surface subsidence and the relationship between the surface subsidence and groundwater changes in Beijing were analyzed by combining ENVISAT ASAR data, leveling data, GRACE data, etc. The main conclusions are as follows: (  Figure 11. Distribution of the main groundwater depression (black curve) [58] and the LOS average subsidence velocity in the study area.

Conclusions
In this study, we adopted the multi-temporal InSAR technique to derive subsidence information for Beijing from 1 August 2007 to 29 September 2010. The surface subsidence and the relationship between the surface subsidence and groundwater changes in Beijing were analyzed by combining ENVISAT ASAR data, leveling data, GRACE data, etc. The main conclusions are as follows: (1) The surface subsidence in Beijing is notably uneven. Five significant subsidence areas that exist in Beijing are Changping, Shunyi, Chaoyang, Tongzhou and Daxing, among which Tongzhou is the most serious subsidence area, with a maximum velocity exceeding 140 mm/year. The subsidence in the downtown area of Beijing is relatively stable, where most subsidence velocities are less than 10 mm/year. Furthermore, the surface subsidence presents nonlinear deformation. (2) A statistical analysis of the standard deviations of the average velocities indicated that 96.81% of the monitoring point standard deviations are less than 5 mm/year, and the standard deviation of the multi-temporal InSAR-derived subsidence in the study area is 1.99 mm/year. In addition, the multi-temporal InSAR-derived results were validated with leveling data. The mean error and RMSE are 1.44 mm/year and 4.97 mm/year, respectively, demonstrating that the multi-temporal InSAR technique is effective for surface subsidence monitoring in Beijing. (3) The groundwater changes derived from the GRACE data in Beijing show a decreasing tendency and profound obvious seasonal variability. Based on the average subsidence of the study area, the long-term decreasing trends in groundwater and average subsidence are consistent. In addition, the spatial distribution of the subsidence funnel only partially overlaps the groundwater depression region. The formation and development of the surface subsidence in Beijing are seriously affected by groundwater over-exploitation.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
In this study, the Kriging method was implemented through ArcGIS 10.2 software and its geostatistical analyst extension to estimate the subsidence of a PS point corresponding to the nearest leveling point. Based on the multi-temporal InSAR-derived PS point subsidence, the spherical variogram model γ (h) was obtained as follow: where h is the distance in meters. In the spherical variogram model, the nugget, sill and major range are 0, 1.077 and 402.338, respectively. Figure A1 shows the spherical variogram model. providing leveling data. The authors also thank TUDelft for DORIS, JPL/Caltech for ROI PAC, and Andy Hooper for StaMPS/MTI.
Author Contributions: All authors contributed to the manuscript and discussed the results. Jiming Guo and Lv Zhou developed the idea that led to this paper. Lv Zhou performed all data processing and analyses, and contributed to the manuscript of the paper. Jiming Guo provided critical comments and contributed to the final revision of the paper. Chaolong Yao performed the GRACE data processing. Jiyuan Hu revised the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
In this study, the Kriging method was implemented through ArcGIS 10.2 software and its geostatistical analyst extension to estimate the subsidence of a PS point corresponding to the nearest leveling point. Based on the multi-temporal InSAR-derived PS point subsidence, the spherical variogram model   where h is the distance in meters. In the spherical variogram model, the nugget, sill and major range are 0, 1.077 and 402.338, respectively. Figure A1 shows the spherical variogram model. Figure A1. The spherical variogram model.

Appendix B
The acronyms used in the text are listed in Table B1.  Figure A1. The spherical variogram model.

Appendix B
The acronyms used in the text are listed in Table B1. Table B1. List of acronyms used in the text.