Spatiotemporal Evolution of Land Subsidence in the Beijing Plain 2003–2015 Using Persistent Scatterer Interferometry (PSI) with Multi-Source SAR Data

Abstract

Land subsidence is one of the most important geological hazards in Beijing, China, and its scope and magnitude have been growing rapidly over the past few decades, mainly due to long-term groundwater withdrawal. Interferometric Synthetic Aperture Radar (InSAR) has been used to monitor the deformation in Beijing, but there is a lack of analysis of the long-term spatiotemporal evolution of land subsidence. This study focused on detecting and characterizing spatiotemporal changes in subsidence in the Beijing Plain by using Persistent Scatterer Interferometry (PSI) and geographic spatial analysis. Land subsidence during 2003–2015 was monitored by using ENVISAT ASAR (2003–2010), RADARSAT-2 (2011–2015) and TerraSAR-X (2010–2015) images, with results that are consistent with independent leveling measurements. The radar-based deformation velocity ranged from −136.9 to +15.2 mm/year during 2003–2010, and −149.4 to +8.9 mm/year during 2011–2015 relative to the reference point. The main subsidence areas include Chaoyang, Tongzhou, Shunyi and Changping districts, where seven subsidence bowls were observed between 2003 and 2015. Equal Fan Analysis Method (EFAM) shows that the maximum extensive direction was eastward, with a growing speed of 11.30 km²/year. Areas of differential subsidence were mostly located at the boundaries of the seven subsidence bowls, as indicated by the subsidence rate slope. Notably, the area of greatest subsidence was generally consistent with the patterns of groundwater decline in the Beijing Plain.

1. Introduction

Widespread surface deformation caused by over-exploitation of groundwater is a worldwide geological and environmental problem [1]. The potential consequences of urban land subsidence mainly include damage to the utility infrastructures, buildings, railroads, highways and bridges [2,3,4], as well as flooding in coastal cities [5,6]. Obtaining accurate information on the spatiotemporal evolution of subsidence is important to reduce the impact of the associated hazards.

Over the past decades, surface deformation monitoring has been rapidly improved thanks to conventional methods, such as precise leveling and GPS (Global Positioning System) [7,8,9], which are, however, time consuming, point-based and lacking in fine details. Since the beginning of the new millennium, space-borne observation techniques based on Interferometric Synthetic Aperture Radar (InSAR) have been applied to measure and map subsidence. The capability of InSAR to measure ground movement over large areas with an accuracy on the order of millimeters has significantly improved knowledge of the phenomenon. As a traditional InSAR method, Differential Interferometric SAR (D-InSAR) has been effective in arid regions and sparse vegetation areas [10,11]. However, the method has challenges for monitoring subtle surface deformation due to temporal and spatial decorrelation and atmospheric disturbances [12,13].

To overcome the limitations of conventional InSAR methods, Persistent Scatterers InSAR (PSI) techniques have been presented in recent years. Ferretti et al. (2001) developed the concept of Permanent Scatterer Interferometry, and Werner et al. (2003) improved Interferometric Point Target Analysis (IPTA), which can reduce spatial and temporal decorrelation, reduce the error that occurs as a result of atmospheric delay, and enhance temporal and spatial resolution [14,15,16]. Hooper et al. (2004) proposed the Stanford Method for Persistent Scatterers (StaMPS) to identify Persistent Scatterers (PS) pixels and analyzed the phase terms, which established a PS recognition algorithm by amplitude dispersion characteristics and spatially correlated features of interferometric phase [17,18]. PSI techniques have been applied in many regions to measure time series surface deformation by resolving general problems in interferometry, e.g., the difficulties in identifying long-term stable targets, increasing the available quantity of image pixels and improving the temporal resolution of the interference [19,20,21].

Land subsidence has been documented in many major cities in the Chinese mainland, including Shanghai [22,23], Tianjin [24,25] and Xi’an [20]. Beijing, the capital of China, is one of the largest cities in China with a population of 21,520,000 people and an area of 16,410 km². It suffers from a serious water shortage. As reported by the China Geological Survey, groundwater has accounted for more than two-thirds of Beijing’s water supply in recent decades, and groundwater withdrawals have induced a rapid increase in land subsidence [26].

Land subsidence in Beijing has been studied using a variety of monitoring techniques. Gong et al. [27] employed the PS-InSAR method to monitor land subsidence in Beijing with ENVISAT ASAR data from 2003 to 2006. Their results showed that differential subsidence is affected by the quaternary sediment thickness. Zhang et al [28] integrated a leveling survey, borehole extensometers and multilayer monitoring of groundwater to characterize the land subsidence of the Beijing Plain. The results showed that the second (64.5–82.3 m) and third (102–117 m) aquitards contributed 39% of the total subsidence. Chen et al. [29] employed ENVISAT ASAR and TerraSAR-X images to investigate land deformation in Beijing region, and they found that the maximum subsidence occurs in the eastern part of Beijing with a rate more than 100 mm/year. Li et al. [30] analyzed land subsidence in Beijing with two adjacent tracks of C-band ENVISAT ASAR data from 2003 to 2006 using the InSAR method. Their results showed that the maximum subsidence rate is 110 mm/year. Zhu et al. [31] collected land subsidence from PSI, leveling survey, hydrogeological and Landsat TM images to detect the spatial and temporal features of land subsidence, and the results indicated that subsidence is unevenly distributed, but continuously increased from 2003 to 2010. Hu et al. [32] employed the techniques of multi-temporal InSAR with C-band ENVISAT ASAR images to produce a linear deformation rate map. They found there were three large subsidence funnels in the Beijing Plain area that are separately located in the Chaoyang and Shunyi districts. Zhang et al. [33] used multiple SAR stacks to characterize subsidence in the Beijing–Tianjin–Hebei region from 1992 to 2014, but their study lacks spatiotemporal analysis of the subsidence characteristic in Beijing. Ng et al. [7] used PSI techniques with ENVISAT ASAR and ALOS PALSAR images from 2003 to 2009 to investigate the ground deformation in Beijing, and the results indicate that the vertical deformation rates were in the range of −115 mm/year to 6 mm/year.

In summary, although previous investigations have employed different methods of monitoring land subsidence in Beijing [7,27,28,29,30,31,32,33,34], few studies have investigated the spatiotemporal evolution of land subsidence in the Beijing Plain using long time series data. This work aims to improve knowledge of the land subsidence in the Beijing Plain and provide a quantitative, long-term analysis of land displacement at regional and local scales. In our study, spatiotemporal evolution of the ground displacements is quantified by the PSI technique on multi-source SAR data (ENVISAR ASAR images from 2003 to 2010, RADARSAT-2 images from 2011 to 2015 and TerraSAR-X images from 2010 to 2015). Time series deformation is analyzed together with the geological structures using Geographic Information System (GIS) spatial analysis. The expansion of the subsidence region in different directions is quantified by the Equal Fan Analysis Method (EFAM). Last, the potential causes of the ground movements, as well as potential future advances in monitoring land subsidence in the Beijing Plain, are discussed.

2. Study Area and Data Sets

2.1. Description of the Beijing Plain

Beijing (115°25′–117°35′E and 39°28′–41°05′N) is located in the extreme northern part of the North China Plain, and is bordered by the Taihang Mountains to the west and Jundu Mountains to north and northeast. The study area includes Beijing and the surrounding plain, with a total area of approximately 6300 km². The Beijing Plain is formed by the alluvial deposit of five rivers: Yongding, Chaobai, Wenyu, Jijun and Juma rivers, with a slight decrease in elevation from northwest to southeast. The elevations of the Beijing Plain ranges from about 20–60 m [29]. Climatically, Beijing has a monsoon-influenced semi-arid and semi-humid continental climate, characterized by high temperatures in the summer, low temperatures in the winter and abundant precipitation. The average annual precipitation is approximately 570 mm, concentrated in the summer season from June to September [35]. The thickness of quaternary deposits, with lithologies that include single gravel to multilayer structures of clay with sand interbeds and a grain-size that ranges from coarse to fine sand, varies from a few meters over the mountain-front area to hundreds of meters in the plain zone [28,31]. The basement rocks are Cambrian–Ordovician in age [28]. The North China Plain has been previously thought to be dominated in the Cenozoic by eastwest extension along north-northeast trending normal faults, creating a half-graben and host system [36]. However, recent work using GPS suggests interseismic deformation across a broad zone, about 1100 km wide, of left-lateral shear, with compensating right-lateral strike-slips along the NNE trending faults [37].

There are several active high angle normal faults in the Beijing Plain (Figure 1). The Nankou–Sunhe fault (Fa-1) is the single NW-striking active fault, which was precisely determined by shallow seismic method and drilling data [38]. The Nankou–Sunhe fault has shown both stick-slip and creep-slip characteristics, with a vertical movement rate of 0.31–0.67 mm/year in the Holocene. Huangzhuang–Gaoliying (Fb-1), Liangxiang–Shunyi (Fb-2), Nanyuan–Tongxian (Fb-3), Lixian–Niubaotun (Fb-4) faults follow a NE–SW direction. The Huangzhuang–Gaoliying fault dips to the SE with an angle of 55°–75°. The northern part of this fault is the largest, hidden and active fault in the Beijing Plain. The fault moved through the Holocene at an average rate of 0.035 mm/year [39]. Ground fissures with a throw of 15–30 cm have been observed in Beiqijia County, north of Beijing, and are assumed to be surface expressions of the Huangzhuang–Gaoliying fault [38]. The Liangxiang–Shunyi and Nanyuan–Tongxian faults dip to NW with an angle of 60°~80° and 50°~75°, respectively [40,41]. The northern part of the later fault is active fault which the rate of neotectonic activity is around 0.75 mm/year [42]. Overall, the Beijing Plain is generally divided into elongated areas by these active faults.

2.2. Datasets Used

Ground movements were measured with multi-source spaceborne SAR datasets: 48 C-band ENVISAT-ASAR, 40 C-band RADARSAT-2 and 51 X-band TerraSAR-X images, all collected in stripmap mode. The ENVISAT-ASAR images were acquired between 10 December 2003 and 7 August 2010. The RADARSAT-2 images were collected from 22 November 2010 to 14 December 2015. The TerraSAR-X images collected from 13 April 2010 to 11 December 2015. ENVISAT-ASAR and RADARSAT-2 data were used to monitor the land deformation information during 2003 to 2015, and the TerraSAR-X data were used to cross-validate the results. The coverage of the SAR images is shown in Figure 1. Additional more parameters of the SAR datasets are summarized in

Table 1. The parameters of the SAR data used.

Parameter	ENVISAT ASAR	RADARSAT-2	TerraSAR-X
Band	C	C	X
Wavelength (cm)	5.6	5.5	3.1
Polarization	VV	VV	HH
Orbit directions	Descending	Descending	Ascending
Track No.	2218	60115	8
Incidence angle (°)	22.9	27.6	33.1
Heading (°)	−164	−168.8	−10
Spatial resolution (m)	30	30	3
No. of images	48	40	51
Date range	June 2003–August 2010	June 2011–December 2015	April 2010–December 2015

.

The DORIS precise orbit state vectors of the ENVISAT-ASAR images were provided by the European Space Research Institute Help Desk of the United Space in Europe (ESA), and used to calculate the initial estimation of interferometric baselines. Shuttle Radar Topography Mission (SRTM) DEM with a spatial resolution of 90 m is used to remove the topographic phase and geocode interferograms. The groundwater level contours was provided by Beijing Water Authority, which is used for comparing the relationship between groundwater and the subsidence.

3. PSI Methods and Data Processing

3.1. Selection of the Master Image and Co-Registration

The PSI technique employs a unique public master image for obtaining the interferometric phase. The choice of the best public master image from the multi-temporal images is important. The first procedure is to order the images by date of acquisition. Then, one image is selected as the unique master image, while the other images are used as the slave images. To guarantee the resulting quality of the interferograms, we select the optimal master image by maximizing the Joint Correlation Index (JCI) [43] of all the images, based on:

$${\rho }^m=\frac{1}{N}{\displaystyle \sum }_{k=1}^{N}c\left(B_{\perp }^{k,m},B_c\right)c\left(T^{k,m},T_c\right)c\left(f_{DC}^{k,m},f_c\right),$$

(1)

where the function c is defined as

$$c(x,a)=\{\begin{array}{ll}1-\frac{\left|x\right|}{a}& (x<a)\\ 0& (x\ge a)\end{array}.$$

(2)

In Equation (1), ${\rho }^m$ denotes the JCI value when $m$ is used as the master image, $B_{\perp }^{k,m}$, $T^{k,m}$ and $f_{DC}^{k,m}$ mean the Spatial Baseline (SB), the Temporal Baseline (TB) and the Doppler Centroid Difference (DCD) respectively between the image $k$ and image $m$. In Equation (2), $a$ denotes the evaluation index of SB, TB or DCD. We set the maximum SB, TB and DCD of all the interferograms as their respective critical values. We assumed every image can be the master image and calculated all of the JCI values by Equation (1).

The image corresponding to the maximum JCI value was chosen as the optimal master image. By using this method the images acquired on 27 June 2007 (ENVISAT ASAR), 2 August 2013 (RADARSAT-2) and 15 December 2013 (TerraSAR-X) were chosen as the master images.

The accurate co-registration of SAR imagery is a key prerequisite for any change detection; all the SAR images have to be co-registered into the same space with sub-pixel accuracy. With precise orbits and reference DEM data, the slave SAR images were co-registered on sampling grids of the selected master image by maximizing correlation of amplitude data between SAR acquisitions.

3.2. Interferogram Formation and PS Pixel Identification

In our study, the interferograms were generated from single-look complex (SLC) by using Delft Object-oriented Radar Interferometric Software (DORIS). The temporal and perpendicular baselines of the interferograms are shown in Figure 2. Then we used the 90 m resolution SRTM DEM to remove the topographic phase contribution from the interferometric phase and geocode the interferograms.

PSI techniques identify Permanent Scatterers (PS) based on the amplitude dispersion index of radar signals. Obtaining sufficiently large amount of PS pixels is critical to the success of this method. We used the software StaMPS/MTI (Stanford Method for Persistent Scatterers/Multi-Temporal InSAR) to identify PS pixels and extract deformation information [18]. StaMPS/MTI methods identify the PS candidates using the amplitude dispersion index (ADI):

$$D_A=\frac{{\sigma }_A}{{\mu }_A}.$$

(3)

In Equation (3), $D_A$ is the ADI, ${\sigma }_A$ and ${\mu }_A$ are the standard deviation (SD) and the mean of a series of amplitude values. In general, a threshold value of $D_A$ is adopted for potential PS candidates, and the number of PS candidates increases with the threshold value. Generally an ADI value of 0.4 or higher is chosen, as a higher threshold is better. In this study, we adopted $D_A$ with a value of 0.5 to generate the largest set of PS candidate pixels and guarantee the quality of them [44]. Then PS pixels were selected from the PS candidate pixels based on the noise levels.

The phase ${\varphi }_{x,i}$ of the $x^{th}$ pixel in the $i^{th}$ topographically corrected interferogram can be expressed as the sum of the phase changes due to surface deformation in the satellite line-of-sight (LOS) direction ${\varphi }_{def,x,i}$, the atmosphere delay ${\varphi }_{atm,x,i}$, the residual phase due to orbit inaccuracies ${\varphi }_{orb,x,i}$, look angle error ${\varphi }_{\epsilon ,x,i}$, and the noise term due to variability in scattering, thermal noise, co-registration errors uncertainty in the position of the phase center in azimuth ${\varphi }_{n,x,i}$ [18].

$${\varphi }_{x,i}={\varphi }_{def,x,i}+{\varphi }_{atm,x,i}+{\varphi }_{orb,x,i}+{\varphi }_{\epsilon ,x,i}+{\varphi }_{n,x,i}.$$

(4)

Using the error characteristics of the individual terms in Equation (4), a measure of the variation of the residual phase for a pixel $\gamma$ as is defined:

$${\gamma }_x=\frac{1}{N}\left|{\displaystyle \sum _{i=1}^N\exp \left\{\sqrt{-1}\left({\varphi }_{x,i}-{\tilde{\varphi }}_{x,i}-{\stackrel{\wedge }{\varphi }}_{\epsilon ,x,i}\right)\right\}}\right|.$$

(5)

In Equation (5), $N$ is the number of the interferograms, ${\tilde{\varphi }}_{x,i}$ is a wrapped estimate of the spatially correlated parts of each of the terms on the right side of Equation (5) and ${\stackrel{\wedge }{\varphi }}_{\epsilon ,x,i}$ denotes the spatially uncorrelated part of ${\varphi }_{\epsilon ,x,i}$. Through spatial–temporal filtering, long-wavelength atmospheric effects and orbit error are removed.

3.3. Phase Unwrapping and Time-Series Deformation Retrieval

Spatial phase unwrapping was applied to recover unambiguous phase values. In order to unwrap correctly, the spatially uncorrelated contribution needs to be subtracted before unwrapping. We use the statistical-cost approach to unwrapping unambiguous phase values [45]. The spatial and temporal band pass filtering was applied to remove the spatially correlated contributions. Finally, the displacement time series for each PS pixel was obtained by least-squares inversion. Time series deformation for each PS pixel was calibrated to a presumably stable reference point.

4. Results and Spatial Analysis

4.1. PSI-Derived Results

The three sets of SAR images were processed by PSI methods and the deformation of the Beijing Plain over the three time periods was calculated. The numbers of PS pixels obtained over the study area were 329,675 for the ENVISAT ASAR images, 773,905 for the TerraSAR-X images and 377,213 for the RADARSAT-2 images. The mean deformation velocity mosaic of the three time spans is shown in Figure 3, respectively. In order to generate a long period information of cumulative deformation, the deformation rate in LOS has been projected into the up–down direction as subsidence rate by Equation (6) [33]:

$${\upsilon }_{u-d}(x,y)={\upsilon }_{LOS}(x,y)/\cos {\theta }_{(x,y)}.$$

(6)

In Equation (6), ${\upsilon }_{u-d}(x,y)$ represents the up–down subsidence rate, ${\upsilon }_{LOS}(x,y)$ is SAR LOS deformation rate, ${\theta }_{(x,y)}$ is SAR viewing angle at position (x, y). Of course, this is only under the ideal condition that the ground just have the deformation of uplift or subsidence. For this study, this assumption is reasonable, as the main deformations in the Beijing Plain are ground subsidence and uplift [46].

The average deformation velocity ranges from −136.9 to +15.2 mm/year from 2003 to 2010 detected by ENVISAT ASAR (Figure 3a), −149.4 to +8.9 mm/year during 2011 to 2015 by RADARSAT-2 (Figure 3b), and −154.7 to 21.4 mm/year during 2010 to 2015 by TerraSAR (Figure 3c). All of the deformation results are relative to the same reference point, which is located at the Tiananmen Square (116.394°E, 39.917°N). It is clear that the main subsidence area is distributed in the east and north parts of the Beijing Plain.

As the locations of the PS pixels generally cannot be exactly the same between ENVISAT ASAR and RADARSAT-2 results, we use Kriging to interpolate all of the discrete PS points into raster images [44]. As the ENVISAT ASAR and RADARSAT-2 SAR images are limited, there is a small time gap (about several months) between the two data sets. We ignored the deformation during the gap in calculating the cumulative deformation. For this study, this assumption is reasonable, because the gap is very short compared to the entire range of analyzed years (2003–2015). Through raster processing, we calculated cumulative deformation from 2003 to 2015 (Figure 12). It is clear that land subsidence occurred across most of the Beijing Plain, increasing in extent and depth in every year during this period, from 2003 to 2015. In the year of 2015, the maximum cumulative subsidence in the Beijing Plain was already more than 1500 mm. From the maximum cumulative subsidence of each year, we can see that there was a gradual acceleration of the deformation rates in the analyzed years.

4.2. The Validation of the PSI Results

There are many methods to validate the deformation results detected from PSI, such as GPS, leveling [47], field photos of ground/building cracks and the estimated heights [48,49]. Considering the limited available data, we chose the leveling measurements and cross-validation of different SAR results to evaluate the quality of the results from PSI.

4.2.1. Comparison between PSI and Leveling Measurements

We compared the time series of deformation from PSI and leveling measurements at Pinggezhuang station (PGZ) and Beiqijiazhen (BQJZ). We chose the nearest PS point of each leveling station for the comparison. Unfortunately, at PGZ there are only leveling measurements from 2008 to 2010 and at BQJZ from 2011 to 2013 (Figure 4). Although a single, long time series data set to cover the whole time span studied is not available, the measured deformation from PSI and Leveling techniques are in general agreement. The difference between the two methods are almost smaller than ±5 mm at most points, and the standard deviation of the difference of ENVISAT ASAR and leveling, RADARSAT-2 and leveling are 3.99 and 6.45, respectively, which indicates that the PSI results are reliable and can meets the accuracy requirements for ground deformation monitoring.

4.2.2. Cross-Validation between Different SAR Results

The PSI outcomes of ENVISAT ASAR and RADARSAT-2 were also validated by statistically comparing the results of TerraSAR-X in the area the two data sets intersect. Because the position of PS from different SAR data sets are generally not located at the same point, we firstly select ASAR and RADARSAT-2 PS points by using a regular grid, where the size of each numbered grid cell is about 500 × 500 m. As the resolution of TerraSAR-X is three meters, there are many TerraSAR-X PS points near the PS points of ASAR or RADARSAT-2. In this study, only the nearest TerraSAR-X PS target of each PS target from ASAR or RADARSAT-2 is chosen for the comparison. After this initial selection, there are 727 PS pairs of ASAR and TerraSAR-X PS points (Figure 5a), as well as 875 PS pairs of RADARSAT-2 and TerraSAR-X PS points (Figure 5c) marked in the grid. The resulting statistical analysis shows there is high correlation between the PSI outcomes from ASAR and TerraSAR-X PS points (Figure 5b), as well as RADARSAT-2 and TerraSAR-X PS points (Figure 5d) during their common time period.

The R-squared of the regression line between ENVISAT ASAR and TerraSAR-X, and the regression line between RADARSAT2 and TerraSAR-X, is 0.86 and 0.97 respectively. The overall good agreement of the results from the different SAR data sets indicates that it is appropriate to use PSI derived subsidence rates from ENVISAT ASAR and RADARSAT-2 for further analysis.

5. Discussion

5.1. Spatiotemporal Evolution of Subsidence Bowls

The mean deformation rate map (Figure 6) was produced from the time-series data of PS points during the study period (2003–2015). It is clear that the main land subsidence region is located in the East and north parts of the Beijing Plain (Haidian and Tongzhou district). The deformation rate ranges from −138.8 to +13.9 mm/year. This serious land subsidence with a rapid rate and a wide range of values has been observed in previous studies. Our deformation results from PSI are similar with the land subsidence results detected by IPTA (Interferometric Point Target Analysis) method [50], SBAS (small baseline subset) method [51,52], and Multi-Temporal InSAR [26]. There is a similar maximum subsidence rate about 140 mm/year during 2003 to 2010 in the Beijing Plain.

As there is a clear distinction across the −40 mm/year deformation line in the map shown in Figure 6, we used the −40 mm/year deformation line as a boundary of the subsidence regions. Combining with the −60 mm/year deformation line, we found there are mainly seven subsidence bowls in the Beijing Plain: Chaoyang–Tongzhou (CT), Tongzhou (TZ), Chaoyang–Jinzhan (CJ), Chaoyang–Laiguangying (CL), ShunYi (SY), ChangPing (CP) and Haidian (HD) Subsidence Bowl. All of these subsidence bowls are around the downtown of Beijing to the north and east sides of it, except for the SY subsidence bowl which is located in the Shunyi district. CT, TZ, CJ and CL subsidence bowls are approximately circular in shape, whereas SY, CP and HD subsidence bowls are close to elliptical. SY subsidence bowl is elongated in the north and south directions, CP towards southeast and northwest directions, almost parallel to the Fb-2 fault, and HD is elongated in a northeast to southwest direction.

It is evident that some subsidence bowls are separated by the major faults in the plain: CT, TZ and CJ are divided by the Fa-1 and Fb-3 faults; Fb-2 crosses in the middle of CJ and CL; the boundaries of CL are Fa-1, Fb-1 and Fb-2; and both CL and SY are located between Fb-1 and Fb-2, but are separated by Fa-1, while there are also some exceptions, such as Fb-3 fault passes through the CT subsidence bowl. That means there may be some other geological structure reasons lead to this phenomenon. On the whole, these spatial characteristics confirm that the influence of faults on land subsidence in the Beijing Plain cannot be ignored.

In order to analyze the development of the subsidence bowls, points (A–G, Figure 6) located almost in the central part of each of the main seven subsidence bowls, were selected for time series analysis of the deformation (Figure 7).

The largest subsidence (more than 100 mm/year) occurred in the central part of the CT subsidence bowl, near point A. Similarly, the land subsidence velocities of TZ (point B) and CJ (point C) subsidence bowls were more than 80 mm/year. A notable feature of Figure 7 is that the year of 2007 was a clear turning point, before which the cumulative subsidence of all points (A–G) was mostly stable following an initial period of subsidence in early 2004. In contrast, the cumulative subsidence increased rapidly after the year of 2007. The other important turning point is the year 2012, when the subsidence speed of D and E slowed down relative to the other subsidence bowls, resulting in lower overall subsidence compared to the other bowls by the end of 2015, and the cumulative subsidence of the two bowls was smaller than the other bowls in the end of 2015.

We examined deformation profiles through time in the northwest–southeast (AA’) and southwest–northeast (BB’) directions. The cumulative deformation of the profiles are shown in Figure 8 from 2003 to 2015. Profile AA’ which is almost along Nankou–Sunhe fault across subsidence bowls CP, CL, CJ and TZ. It is apparent that Fb-1 is located between CP and CL, Fb-2 is located between CL and CJ, and Fb-3 is located between CJ and TZ, which means that the subsidence area is divided into regions along the northwest–southeast direction. From the profile BB’ we can see that there was almost no subsidence in the left part (i.e., southwest), but there was large subsidence (about 500 mm from 2003 to 2015) in the right part (i.e., northeast). We can conclude that there are generally two parts (stable and subsidence areas) along the southwest–northeast direction.

5.2. The Area of Different Mean Deformation Rates Region

From the GIS spatial analysis, we discovered that the most serious part of land subsidence in the Beijing Plain is located in the east part of Chaoyang district, a 26.51 km² region with a subsidence rate of more than 100 mm/year. The area of subsidence rate with more than 40 mm/year is 536.06 km² large (Figure 9).

The Equal Fan Analysis Method (EFAM) was used to define quantitatively the expansion of the subsidence region in different directions. This method has previously been used to study the evolution of city development [53,54]. EFAM is an effective way to explore the changing features in different directions of an object. In our study, we use this method to expose the spatial evolution of land subsidence in the Beijing Plain as a whole. First, we made Tiananmen Square as the center point, and drew a circle with a 50 km radius, covering most of the subsidence area in the Beijing Plain. Then we divided the circle into sixteen equal fans by the EFAM method. Each sector has its own direction, such as N, NNE, NE, NEE, E, and so on. The area of subsidence in each fan represents the expansion of the subsidence in the corresponding direction.

The sixteen radial subsidence regions are shown in the Figure 10a and listed in

Table 2. The Expanding Area in Different Directions.

Direction	Expanding Area (km2)	Rate (km2/Year)	Direction	Expanding Area (km2)	Rate (km2/Year)
N	24.46	1.88	S	0	0
NNE	56.32	4.33	SSW	0	0
NE	63.26	4.87	SW	0	0
NEE	71.81	5.52	SWW	0	0
E	146.87	11.30	W	0	0
SEE	95.66	7.36	NWW	0	0
SE	27.53	2.12	NW	20.11	1.55
SSE	0	0	NNW	26.93	2.07

, Figure 10b summarizes that the area expansion in different directions using a radar map. The maximum extension is in the E direction, at about 146.83 km², and the expansion rate is nearly 11.30 km²/year. The land subsidence area (314.34 km²) in the E, NEE and SEE directions occupy approximate 66% of the total subsidence area. The average of the expansion in the northeast direction is relatively higher than it on the west and south directions. This means that Shunyi and Tongzhou districts are particularly affected by land subsidence and Chaoyang district is the most affected region of land subsidence. In contrast, there is almost no subsidence in the south and west part of the Beijing Plain. The results of subsidence area expansion observed by EFAM are not only similar to the previous findings of Chen et al. (2016), the subsidence bowls in the east, northeast and north parts of Beijing have gradually merged [29], but also show the average speed of the expansion in the northeast direction during 2003 to 2015.

5.3. Differential Subsidence Zone

Differential subsidence, which implies a large difference in deformation within a short distance, may be particularly hazardous for buildings and infrastructure [55,56]. It is therefore important to locate differential subsidence areas. In order to recognize the area of differential land subsidence, we adopted the concept of slope in geographic spatial analysis. The Subsidence Rate Slope (SRS) was defined as:

$$I_s=\frac{\Delta V}{l}.$$

(7)

In Equation (7), $\Delta V$ represents the change of land subsidence rate, and $l$ stands for the distance.

We assume that the region with a larger difference rate than 7 mm/year over 100 m in the horizontal direction is a differential subsidence area. Thus, we normalize the value of SRS to the range 0 to 1, and divide the differential subsidence region into different levels for a better distinction by using the subsidence slope map (Figure 11). Not surprisingly, the differential subsidence areas are mainly distributed on the boundary of the subsidence bowls. These regions expose larger structures to differential subsidence rates, and this may lead to structural failure. In addition, when we added the subway lines on the differential subsidence map, we found the east parts of subway line Six and line One plus Batong are mostly located in the differential subsidence area. The Airport Express line is also goes across the differential subsidence region.

5.4. Subsidence and Groundwater

Long-term over-exploitation of groundwater could cause a regional drop of piezometric levels, decreasing pore pressure and increasing effective stresses [56]. Previous research showed that the predominant cause of subsidence in the Beijing Plain is intense groundwater extraction, rather than the other driving mechanisms [31].

Combining ground water level contours from 2003 to 2015 and the cumulative deformation extracted by PSI during the period, the correlation between subsidence response patterns and phreatic groundwater flow field was analyzed comprehensively (Figure 12). A comparison of 2015 groundwater levels with 2004 levels shows water level declines of 15–30 m throughout the north plain, with a 50 m maximum decline occurring in the central part of the Beijing Plain.

In the year of 2004, the groundwater elevations ranged from 10 to 30 m and there was only a small groundwater bowl on the northwest part of the plain due to long term over exploitation of groundwater. However, the groundwater contours in 2015 indicate a more complex pattern than those of 2003 with many groundwater bowls in the middle of the Beijing Plain and a minimum value of the groundwater of −20 m. The cumulative deformation also had great change over the same period with a maximum subsidence of more than 1500 mm. The most serious and extensive subsidence in the study area occurred in the area with the greatest groundwater depression, particularly in the east and north parts of the Beijing Plain. Furthermore, it also shows that the depression of the groundwater and the response of subsidence in the Beijing Plain is generally consistent. The Chaoyang–Jinzhan, Chaoyang–Laiguangying, ShunYi and ChangPing subsidence bowls were generally located in the same area of water level decrease. However, the spatial distribution of the subsidence zones are not fully aligned with the groundwater depression area. There are also some differences between them, the Chaoyang–Tongzhou and TongZhou subsidence bowls are not located in the zone of declining groundwater.

It is worth noting that the displacements on the ground surface are the total deformation that occurred in the different underground aquitards and aquifers. For more detailed analysis about the relationship between subsidence and groundwater, we compared the time series analysis of deformation measured by PSI and groundwater level change (Figure 13).

Four representative monitoring wells (BLG, DJ, MJQ and MT) were selected from study area to analyze the temporal response of land subsidence through a comparison with the time series deformation of the nearest PS point. BLQ and DJ wells are located in the Chaoyang–Tonghzou and Chaoyang–Laiguangying subsidence bowls respectively, while MJQ and MT well are located in the south and southeast part of study area, which are a little far from the subsidence area (marked in Figure 6). The phreatic water levels of BLQ and MJQ have stabilized rephase and water level of MT recovered by as much three meters during the study period. The water levels of confined aquifers continued to decline with the maximum 18 m decrease in the BLQ, and the corresponding subsidence was about 300 mm for the period 2005–2010. The water levels of the phreatic and confined aquifers at DJ location were generally have the same trend, dropping from 2005 to 2008, and recovering slightly after 2008. The comparison shows a close correlation between water-level fluctuations and seasonal deformation response of the location of MJQ. Seasonal displacement was also found at the MT well. Furthermore, the deformation show an obvious delayed response (about three months) to the seasonal head change which can be attributed to the delayed drainage and fluid-pressure equilibration of the low-permeability aquitard or interbed. The displacements at BLQ and DJ wells continued at relative stable rates through most of the period, even though the water level recovered with a relative large magnitude (DJ). The BLQ and DJ wells show a mainly inelastic behavior as the deformation presents a continuous compaction with the groundwater level decline. The inelastic aquifer-system deformation may be primarily caused by two factors: residual compaction and secondary compression of the aquitards.

6. Conclusions

In this study, land subsidence due to over-extraction of groundwater in the Beijing Plain was investigated with the PSI method by using multi-source SAR (ENVISAT ASAR, RADARSAT-2 and TerraSAR-X) images. The PSI results were compared with conventional leveling measurements and cross-validation, which showed that the PSI monitoring results replicate both the seasonal and long-term trends of the conventional benchmark results with an accuracy of about five millimeters. The displacement map reveal reveals that the Beijing Plain has experienced severe land subsidence from June 2003 to December 2015, with a maximum cumulative subsidence of about 1508 mm.

From the average deformation velocity fields during 2003 to 2015, seven major subsidence bowls were identified: Chaoyang–Tongzhou (CT), Tongzhou (TZ), Chaoyang–Jinzhan (CJ), Chaoyang–Laiguangying (CL), ShunYi (SY), ChangPing (CP) and Haidian (HD) Subsidence Bowls, in the east and north parts of the Beijing Plain.

The spatiotemporal analysis of land subsidence indicates an increasing trend in the rate and extent of land subsidence. A profile across the northeast-trending Quaternary faults (Figure 6) reveals that the distribution of subsidence bowls is clearly controlled by the faults. The Equal Fan Analysis Method was used to detect the expansion in different directions of the subsidence region, showing that the maximum expansion direction is towards the east, with a speed of 11.30 km²/year.

The comparisons of land subsidence and head changes of multi-aquifer system between 2003 and 2015 shows that the main subsidence region is distributed in the same area with the area of groundwater depression. The comparison between the time-series deformation and water-level fluctuations shows the different response patterns of the aquifer system. Seasonal displacement was also found at the MT well. However, the deformation show an obvious delayed response to the seasonal head change which can be attributed to the delayed drainage and fluid-pressure equilibration of the low-permeability aquitard or interbed.

To summarize, mapping the land subsidence by PSI techniques with multi-source SAR data has allowed us to characterize the spatial and temporal pattern of the land subsidence in the Beijing Plain during the study period. The accuracy of deformation could be further improved by increasing the SAR data, and by combining data from different SAR viewing geometries and multiple satellites. The greater resolution of the seasonal signals may be detected and can be used to estimate the storage change of aquifer systems. Time-series displacements of PSI may also be used to identify the depth-dependent mechanisms of land subsidence in conjunction with conventional hydrogeophysical methods.