Imaging Land Subsidence Induced by Groundwater Extraction in Beijing (China) Using Satellite Radar Interferometry

Beijing is one of the most water-stressed cities in the world. Due to over-exploitation of groundwater, the Beijing region has been suffering from land subsidence since 1935. In this study, the Small Baseline InSAR technique has been employed to process Envisat ASAR images acquired between 2003 and 2010 and TerraSAR-X stripmap images collected from 2010 to 2011 to investigate land subsidence in the Beijing region. The maximum subsidence is seen in the eastern part of Beijing with a rate greater than 100 mm/year. Comparisons between InSAR and GPS derived subsidence rates show an RMS difference of 2.94 mm/year with a mean of 2.41  ̆ 1.84 mm/year. In addition, a high correlation was observed between InSAR subsidence rate maps derived from two different datasets (i.e., Envisat and TerraSAR-X). These demonstrate once again that InSAR is a powerful tool for monitoring land subsidence. InSAR derived subsidence rate maps have allowed for a comprehensive spatio-temporal analysis to identify the main triggering factors of land subsidence. Some interesting relationships in terms of land subsidence were found with groundwater level, active faults, accumulated soft soil thickness and different aquifer types. Furthermore, a relationship with the distances to pumping wells was also recognized in this work.


Introduction
Located at the northern end of the North China Plain, Beijing is ranked as the 5th most water-stressed city in the world [1], and groundwater is the main water source for industrial, agricultural and household activities.With its rapid urban growth, there has been increasing water demand in Beijing.Previous studies [2,3] reveal that the Beijing region has been suffering from land subsidence due to over-exploitation of groundwater since 1935, and more seriously, the rate and extent of land subsidence shows an increasing trend.Land subsidence is a severe geohazard threating the safety of the public and urban infrastructure.Hence, continuous monitoring of land subsidence is critical for detecting potential hazards and designing compensation strategies.
Interferometric Synthetic Aperture Radar (InSAR) is a powerful tool widely used to map land subsidence over wide regions with high spatial-temporal resolutions.InSAR time series techniques such as Persistent Scatterer Interferometry (PSI) [4,5] and the Small Baseline Subset (SBAS) [6][7][8] can minimise the limitations of traditional InSAR (e.g., spatial and temporal decorrelation and atmospheric effects) and have been demonstrated to be able to map ground displacements with precision comparable to traditional geodetic techniques such as leveling [9][10][11][12].In this work, the Small Baseline InSAR technique has been employed to investigate land subsidence in the Beijing region and its relationship with different triggering and conditioning factors.
Gong et al. [13] explored the application of PSI method to monitor land subsidence in Beijing with C-band Envisat ASAR data acquired from 2003 to 2006.Their results showed that the spatial pattern of the subsidence is influenced by Quaternary faults, and land subsidence often develops in areas with clay layer thicker than 50 m.Li et al. [14] employed SBAS InSAR to investigate city subsidence in Beijing with two adjacent tracks of C-band Envisat ASAR images, reporting a maximum subsidence rate reaching 100 mm/year in Chaoyang district.Chen et al. [15] studied the spatial and temporal evolution of the groundwater depression funnels in the Beijing plain and the 2D and 3D spatial evolution of land subsidence from a model based on the combination of piezometric, GPS, and radar data.Liang et al. [16] monitored the accumulated crustal deformation and its characteristics in Beijing and its surrounding areas by traditional D-InSAR method with L-band SAR data.They reported that the maximum deformation rate along Line of Sight (LOS) was about 125 mm/year between 2007 and 2010.Ng et al. [17] applied PS InSAR to 41 Envisat ASAR images (acquired from 2003 to 2009) and 24 ALOS PALSAR images (acquired from 2007 to 2009) to investigate the land deformation rate maps over Beijing showing that the vertical displacement rates were in the range of ´115 to 6 mm/year.Zhu et al. [18] analyzed the spatial relationship between land subsidence and three factors (groundwater level, compressible soil thickness and building areas) by means of remote sensing and Geographical Information System (GIS) tools and used the Back Propagation Neural Network (BPN) model combined with Genetic Algorithm (GA) to simulate regional distribution of land subsidence.Hu et al. [19] used 52 Envisat ASAR images acquired from 2003 to 2010 to derive land deformation information in Beijing, and found a maximum subsidence velocity of up to 110 mm/year.Zhang et al. [20] utilized stress-strain analysis and oedometer tests to characterize the hydraulic and mechanical properties of the five hydrogeologic units at different depths in the Beijing plain.The results revealed that the second (64.5-82.3m) and third (102-117 m) aquitards contributed 39% of the total compression deformation, and the second (82.3-102 m) and third (117-148 m) confined aquifers exhibited elasto-plastic mechanical behaviour.Zhu et al. [21] used SAR images, optical images and hydrogeological data to study the land subsidence in the northern Beijing plain.PS InSAR was applied to C-band Envisat ASAR images acquired from 2003 to 2010 to estimate land subsidence on this area.The results indicated that the largest subsidence rates reached 52 mm/year and the silty clay layers contribute to the larger land subsidence.
Compared with medium resolution SAR imagery such as C-band Envisat ASAR and L-band ALOS PALSAR data, recent advanced imaging radar sensors such as the TerraSAR-X SAR system (X-band radar wavelength and a shorter revisiting interval) can provide higher spatial resolution and geometric accuracy imagery, which makes it possible to monitor land subsidence at a finer scale [22,23].Chen et al. [24] utilized C-band Envisat ASAR data and X-band COSMO-SkyMED images to detect land subsidence rates in Shanghai and validated the results by spirit leveling measurements.The results indicated the location of a significant subsidence funnel and the root mean square error of differences of the point targets between two band SAR images was 3 mm/year.Luo et al. [25] extracted land deformation information from L-band ALOS PALSAR data between 2007 and 2010 and X-band TerraSAR-X data from 2009 to 2010.The obtained results indicated that the subsidence rate was in the range of ´90 mm/year to ´10 mm/year from TerraSAR-X data and ´190 to ´10 mm/year from PALSAR data, which was in good agreement with leveling and GPS data.
Although existing investigations have tried different ways of monitoring land subsidence in the Beijing region [13][14][15][16][17][18][19][20][21], the consistency of the land subsidence monitoring results in this densely populated area in the world acquired by multi-wavelength InSAR time series analysis and validation by GPS data is rarely covered.Furthermore, the comprehensive and quantitative spatio-temporal analysis to identify the main conditioning and triggering factors of land subsidence still need to be systematically researched further.This paper aims to use TerraSAR-X (X-band) stripmap images collected from 2010 to 2011 and Envisat ASAR (C-band) images acquired between 2003 and 2010 to provide a new insight into the spatial and temporal distribution characteristics and main conditioning and triggering factors of land subsidence in the Beijing plain.The joint cross analysis of radar-derived deformation data and existing geo-information will help to improve the knowledge of the mechanisms that govern land subsidence to be used for the development of suitable groundwater management policies in this area.The paper is organized as follows: Section 2 describes the geographical and geological setting of the study area.SAR images and InSAR processing used in this work is described in Section 3. Section 4 shows the main InSAR-derived products, validation of the results and the spatial-temporal characteristics of land subsidence over the study area.Section 5 discusses and analyses the effect of different conditioning and triggering factors on land subsidence.The main conclusions are summarized in Section 6.

Description of the Beijing Basin
Beijing is the political, cultural and economic centre of China.As one of the most populated cities in the world, its population has reached 20 million [26].As one of the most water-scarce cities in the world, the major water source for the Beijing municipality is groundwater, accounting for about two thirds of water use [27].
Beijing is located in northern China (Figure 1a), at 39 ˝28 1 -41 ˝05 1 north latitude and 115 ˝25 1 -117 ˝30 1 east longitude, with an area of 16,807 square kilometers, of which 6390 square kilometers are the plain.The northern and north-eastern parts of Beijing are dominated by (Figure 1b) the Jundu Mountains, while the western part of Beijing is framed by the Taihang Mountains.The southeastern part of Beijing is alluvial-pluvial plains crossed by five rivers (Figure 1b): Chaobai, Yongding, Wenyu, Ju and Juma.The terrain of Beijing is high in the west-north area but low in south-east part, sloping from the front of the mountain to the south-east.The elevations range from around 60-80 m.a.s.l. in the front of mountains to around 20-60 m.a.s.l. in the plains.
The study area is characterized by a monsoon-influenced semi-arid and semi-humid continental climate [20].The average annual temperature is around 11.7 ˝C.Annual precipitation, mainly concentrated in summer months (from June to September), is around 570 millimetres and constitutes the main source of water resources in Beijing [27].
The thickness of quaternary deposits, whose lithologies change from single gravel to multilayer structures of clay with sand interbeds and whose grain-size turns from coarse to fine, varies from a few meters over the mountain-front area to hundreds of meters in the plain zone [20,21].Groundwater levels change from deep to shallow and the groundwater varies greatly from phreatic to multilayer confined water throughout the study area, increasing from tens of meters in the mountain-front zone to several hundred meters in the central or southeast plain area [28].
The sediment, forming the principal aquifer system units, predominantly consists of gravel, fine-to-coarse sand, silt and clays (Figure 1c). Figure 1b,c illustrate the gradual transition from fine-to-coarse sand and gravel in the proximal fan to fine sand, silt or clay in the mid and distal-fan areas.The Quaternary aquifers of Beijing plain can be divided into three aquifer groups [2,20,29]: (a) The first aquifer group (late Pleistocene-Holocene) is widespread over the Beijing plain.The depths of the bottom of this group are less than 100 m.In this aquifer group, the phreatic aquifer with single gravel structure located in the top area of the alluvial-diluvial fan and over the middle and lower parts of the alluvial-proluvial fan, distributes phreatic aquifer and the shallow confined water with a multilayer structure; (b) The second aquifer group (middle Pleistocene) is mainly located over the middle and lower segments of the alluvial-diluvial fan, which is the middle and deep confined aquifers with depths of the roof from 80 to 100 m, and depths of the bottom plate around 300 m.This aquifer group presents a multilayer structure, which is mainly made up of medium-coarse sand and small amounts of gravel; (c) The third aquifer group (early Pleistocene) consists on deep confined aquifers with the depth of the top over 300 m.This aquifer group, mainly distributed over the north-east part and the south-east part of the centre of the depression land in Beijing plain, is mainly composed of medium-coarse sand and gravel, with clay aquiclude thicker than 30 m covering the top.
The unconfined aquifer and the first confined aquifer are intensely exploited for irrigation providing about 68% of the total groundwater pumpage of the aquifer system.However, the deeper aquifer units are not significantly exploited [20].The Quaternary aquifers of Beijing plain can be divided into three aquifer groups [2,20,29]: (a) The first aquifer group (late Pleistocene-Holocene) is widespread over the Beijing plain.The depths of the bottom of this group are less than 100 m.In this aquifer group, the phreatic aquifer with single gravel structure located in the top area of the alluvial-diluvial fan and over the middle and lower parts of the alluvial-proluvial fan, distributes phreatic aquifer and the shallow confined water with a multilayer structure; (b) The second aquifer group (middle Pleistocene) is mainly located over the middle and lower segments of the alluvial-diluvial fan, which is the middle and deep confined aquifers with depths of the roof from 80 to 100 m, and depths of the bottom plate around 300 m.This aquifer group presents a multilayer structure, which is mainly made up of medium-coarse sand and small amounts of gravel; (c) The third aquifer group (early Pleistocene) consists on deep confined aquifers with the depth of the top over 300 m.This aquifer group, mainly distributed over the north-east part and the south-east part of the centre of the depression land in Beijing plain, is mainly composed of medium-coarse sand and gravel, with clay aquiclude thicker than 30 m covering the top.
The unconfined aquifer and the first confined aquifer are intensely exploited for irrigation providing about 68% of the total groundwater pumpage of the aquifer system.However, the deeper aquifer units are not significantly exploited [20].Corresponding to the classification of the aquifers, three compressible layers responsible of land subsidence can be recognized in Beijing [2,29].The first compressible layer is widely distributed over the Beijing plain.The depth of the bottom of this layer, which is composed of normally consolidated silt, silty clay and clay with plastic or hard plastic states, is less than 100 m.The second compressible layer is widely distributed over the middle and lower parts of the alluvial-proluvial fan in Beijing.The main lithologics of this layer are silt, silty clay and clay.The depth of the bottom over the south-west part of Beijing is less than 150 m, while the depth of the bottom over the east and north parts of Beijing can reach a depth of 280 m.The compressible over-consolidated layers with plastic or hard plastic clays account for 60~80 percent of the overall thickness.The third compressible layer is distributed mainly over the centre of the depression of the Beijing plain.The depth of the roof of this layer is larger than 270 m.This layer is mainly composed of over-consolidated clays with a solid state [28].

InSAR Processing
The SAR dataset used in this study consists of 41 images from Envisat-ASAR acquired from 2003 to 2010 and 14 images from TerraSAR-X (TSX) acquired from 2010 to 2011.The Envisat ASAR images were collected in Stripmap mode from descending track with VV polarisation, and the TSX images were obtained in Stripmap mode from ascending track with HH polarisation (Figure 1).The interferograms were generated from single-look complex (SLC) images using DORIS (Delft object oriented radar interferometric software) [30].Among all the total amount of possible interferograms formed by pairs of images, only 91 interferograms from Envisat ASAR images and 36 interferograms from TerraSAR-X images were selected for further processing.The selection of interferograms was restricted to those interferograms with spatial and temporal baselines smaller than 1070 m and 1500 days for Envisat ASAR images, and temporal baselines smaller than 200 m and 200 days for TerraSAR-X images.The external Shuttle Radar Topography Mission (SRTM) DEM with 90 m resolution was used to remove the topographic component of the interferometric phase and geocode the interferograms.
In this study, land subsidence information was obtained by using a small baseline subset of interferograms with the StaMPS (Stanford Method for Persistent Scatterer/Multi-Temporal InSAR) technique [31].The processing includes three main steps: Firstly, a network of multi-master differential interferograms was created by taking into account small spatial and temporal baselines and reduced Doppler centroid frequency differences to minimize the spatial and temporal decorrelation and topographic errors for the group of small baseline interferograms.Figure 2 shows the network of the small baseline interferograms used for the SBAS analysis of the land subsidence in Beijing.
Corresponding to the classification of the aquifers, three compressible layers responsible of land subsidence can be recognized in Beijing [2,29].The first compressible layer is widely distributed over the Beijing plain.The depth of the bottom of this layer, which is composed of normally consolidated silt, silty clay and clay with plastic or hard plastic states, is less than 100 m.The second compressible layer is widely distributed over the middle and lower parts of the alluvial-proluvial fan in Beijing.The main lithologics of this layer are silt, silty clay and clay.The depth of the bottom over the south-west part of Beijing is less than 150 m, while the depth of the bottom over the east and north parts of Beijing can reach a depth of 280 m.The compressible over-consolidated layers with plastic or hard plastic clays account for 60~80 percent of the overall thickness.The third compressible layer is distributed mainly over the centre of the depression of the Beijing plain.The depth of the roof of this layer is larger than 270 m.This layer is mainly composed of over-consolidated clays with a solid state [28].

InSAR Processing
The SAR dataset used in this study consists of 41 images from Envisat-ASAR acquired from 2003 to 2010 and 14 images from TerraSAR-X (TSX) acquired from 2010 to 2011.The Envisat ASAR images were collected in Stripmap mode from descending track with VV polarisation, and the TSX images were obtained in Stripmap mode from ascending track with HH polarisation (Figure 1).The interferograms were generated from single-look complex (SLC) images using DORIS (Delft object oriented radar interferometric software) [30].Among all the total amount of possible interferograms formed by pairs of images, only 91 interferograms from Envisat ASAR images and 36 interferograms from TerraSAR-X images were selected for further processing.The selection of interferograms was restricted to those interferograms with spatial and temporal baselines smaller than 1070 m and 1500 days for Envisat ASAR images, and temporal baselines smaller than 200 m and 200 days for TerraSAR-X images.The external Shuttle Radar Topography Mission (SRTM) DEM with 90 m resolution was used to remove the topographic component of the interferometric phase and geocode the interferograms.
In this study, land subsidence information was obtained by using a small baseline subset of interferograms with the StaMPS (Stanford Method for Persistent Scatterer/Multi-Temporal InSAR) technique [31].The processing includes three main steps: Firstly, a network of multi-master differential interferograms was created by taking into account small spatial and temporal baselines and reduced Doppler centroid frequency differences to minimize the spatial and temporal decorrelation and topographic errors for the group of small baseline interferograms.Figure 2 shows the network of the small baseline interferograms used for the SBAS analysis of the land subsidence in Beijing.Secondly, after filtering the azimuth and range spectra to reduce decorrelation effects due to geometry and non-overlapping Doppler frequencies, the slowly decorrelating filtered phase (SDFP) pixels whose phase shows slow decorrelation over short time intervals were selected.For computational efficiency, the SDFP candidates were initially selected by setting a threshold for the amplitude difference dispersion, which is the standard deviation of the amplitude difference between the master and slave divided by the mean amplitude.Afterwards, an iteratively conducted phase analysis yielded the phase stability estimation of each candidate to generate the final set of SDFP pixels [31].The amplitude dispersion index threshold was set at 0.6 to generate the largest set of candidate pixels.
Thirdly, after the iteratively robust process of SDFP selection and correcting for spatially uncorrelated look angle errors, a three-dimensional phase unwrapping was applied on the final sets of SDFP pixels by using a statistical cost flow algorithm [32,33].The time series of displacement for each SDFP pixel was derived through a least-squares inversion method [34].This technique estimates displacement time series without prior consideration of a temporal model for the deformation and thus allows the derivation of temporally varying deformation processes.From these values, long-wavelength atmospheric effects and orbit error are estimated from the SDFP pixels using a spatial-temporal filter method.

InSAR-Derived Land Subsidence Maps
Figure 3 shows the result of deformation pattern along line of sight (LOS) for the multi-temporal SBAS processing of the available Envisat ASAR dataset.In this area, a high concentration of SDFP pixels (i.e., an average density 400 SDFP pixels per square kilometre) is recognized mainly due to the dense urbanization of the scene.The Beijing plain area exhibits high displacement rates with an average velocity of up to 110 mm/year in the LOS direction from 2003 to 2010.
Secondly, after filtering the azimuth and range spectra to reduce decorrelation effects due to geometry and non-overlapping Doppler frequencies, the slowly decorrelating filtered phase (SDFP) pixels whose phase shows slow decorrelation over short time intervals were selected.For computational efficiency, the SDFP candidates were initially selected by setting a threshold for the amplitude difference dispersion, which is the standard deviation of the amplitude difference between the master and slave divided by the mean amplitude.Afterwards, an iteratively conducted phase analysis yielded the phase stability estimation of each candidate to generate the final set of SDFP pixels [31].The amplitude dispersion index threshold was set at 0.6 to generate the largest set of candidate pixels.
Thirdly, after the iteratively robust process of SDFP selection and correcting for spatially uncorrelated look angle errors, a three-dimensional phase unwrapping was applied on the final sets of SDFP pixels by using a statistical cost flow algorithm [32,33].The time series of displacement for each SDFP pixel was derived through a least-squares inversion method [34].This technique estimates displacement time series without prior consideration of a temporal model for the deformation and thus allows the derivation of temporally varying deformation processes.From these values, long-wavelength atmospheric effects and orbit error are estimated from the SDFP pixels using a spatial-temporal filter method.

InSAR-Derived Land Subsidence Maps
Figure 3 shows the result of deformation pattern along line of sight (LOS) for the multi-temporal SBAS processing of the available Envisat ASAR dataset.In this area, a high concentration of SDFP pixels (i.e., an average density 400 SDFP pixels per square kilometre) is recognized mainly due to the dense urbanization of the scene.The Beijing plain area exhibits high displacement rates with an average velocity of up to 110 mm/year in the LOS direction from 2003 to 2010.  Figure 4 illustrates the spatio-temporal evolution of the accumulated LOS surface displacements at approximately one-year intervals derived from Envisat ASAR images.It is clear that land subsidence gradually increases with time, mainly in those areas located at the northeast part of the city in which subsidence rates are higher (Figure 3). Figure 4 illustrates the spatio-temporal evolution of the accumulated LOS surface displacements at approximately one-year intervals derived from Envisat ASAR images.It is clear that land subsidence gradually increases with time, mainly in those areas located at the northeast part of the city in which subsidence rates are higher (Figure 3).

InSAR Data Validation
To evaluate the quality of the InSAR derived land subsidence rate map from Envisat ASAR, a comparison with in situ instrumental data acquired by means of GPS measurement was conducted.Six GPS points that approximately cover the studied period (i.e., from 2007 to 2010) were utilized as checkpoints.The GPS data were projected along the LOS according to [35]. Figure 5 shows the scatter plot of GPS and InSAR solutions over the six GPS points.Note that InSAR-derived displacement rates considered for the comparison of both techniques are here calculated while averaging displacement rate values from all SDFP pixels contained in a 50 m vicinity for the regional comparison of deformation patterns.It can be seen that the maximum discrepancy (MaxD) is 5.34 mm/year, and the minimum one (mD) is 0.14 mm/year.The mean absolute discrepancies (MD) at six checkpoints is 2.41 mm/year, the standard deviation of the discrepancies is ˘1.84 mm/year and the root mean square (RMS) is 2.94 mm/year.The data correlation between GPS data and InSAR measurement is over 0.98.This fact indicates that the overall SBAS results agree well with GPS, suggesting the reliability of InSAR-derived land subsidence rates.

InSAR Data Validation
To evaluate the quality of the InSAR derived land subsidence rate map from Envisat ASAR, a comparison with in situ instrumental data acquired by means of GPS measurement was conducted.Six GPS points that approximately cover the studied period (i.e., from 2007 to 2010) were utilized as checkpoints.The GPS data were projected along the LOS according to [35]. Figure 5 shows the scatter plot of GPS and InSAR solutions over the six GPS points.Note that InSAR-derived displacement rates considered for the comparison of both techniques are here calculated while averaging displacement rate values from all SDFP pixels contained in a 50 m vicinity for the regional comparison of deformation patterns.It can be seen that the maximum discrepancy (MaxD) is 5.34 mm/year, and the minimum one (mD) is 0.14 mm/year.The mean absolute discrepancies (MD) at six checkpoints is 2.41 mm/year, the standard deviation of the discrepancies is ±1.84 mm/year and the root mean square (RMS) is 2.94 mm/year.The data correlation between GPS data and InSAR measurement is over 0.98.This fact indicates that the overall SBAS results agree well with GPS, suggesting the reliability of InSAR-derived land subsidence rates.Ten thousand samples of InSAR-derived subsidence rates between multi-wavelength SAR data (i.e., C-and X-bands) were also compared for further validation.Figure 6 shows the correlation and differences between the two data sets used in this work.A high correlation coefficient of 0.92 was observed with an RMS difference of 7.48 mm/year, further indicating the reliability of InSAR retrievals of land subsidence rate maps.
The overall good agreement at least suggests (1) the study region exhibited similar linear motions during these two time intervals; and (2) it is sensible to use InSAR derived subsidence rates (rather than displacement time series) for further analysis (e.g., to identify controlling factors and etc.).Ten thousand samples of InSAR-derived subsidence rates between multi-wavelength SAR data (i.e., C-and X-bands) were also compared for further validation.Figure 6 shows the correlation and differences between the two data sets used in this work.A high correlation coefficient of 0.92 was observed with an RMS difference of 7.48 mm/year, further indicating the reliability of InSAR retrievals of land subsidence rate maps.
The overall good agreement at least suggests (1) the study region exhibited similar linear motions during these two time intervals; and (2) it is sensible to use InSAR derived subsidence rates (rather than displacement time series) for further analysis (e.g., to identify controlling factors and etc.).

Spatio-Temporal Analysis of Land Subsidence
There are more than 2.5 million high coherent points extracted in the study region with a total area of around 8000 km 2 .The reference point, located in the centre of Beijing, shows a relative stable state without significant subsidence appearance, as shown in Figure 7.The settlement amounts are relatively minor over the downtown, west, south-east and south-west parts of Beijing city with a mean subsidence rate of around 10 mm/year.The fillings covering the West part of Beijing mainly consists of coarse-grain and medium-grain sediments (Figure 1) from proximal fans with low compressibility [20], resulting in small settlements.Nevertheless, the sediments over the east, north and south parts of Beijing mainly consist of high compressibility fine-grain sediments [20] with thick quaternary deposits (Figure 1).The aquifers units of these areas have been intensely over-exploited, generating severe land subsidence.
The main subsidence bowls are distributed over the Chaoyang, Changping, Shunyi and Tongzhou Districts.It can be observed in Figures 7 and 8 that the subsidence bowls located in the east, north-east and north parts of Beijing have been gradually connected, indicating the serious situation in these areas.Consequently, this bowl, which appeared for the first time in 1975 [15], has gradually increased in depth and extension until the present.

Spatio-Temporal Analysis of Land Subsidence
There are more than 2.5 million high coherent points extracted in the study region with a total area of around 8000 km 2 .The reference point, located in the centre of Beijing, shows a relative stable state without significant subsidence appearance, as shown in Figure 7.The settlement amounts are relatively minor over the downtown, west, south-east and south-west parts of Beijing city with a mean subsidence rate of around 10 mm/year.The fillings covering the West part of Beijing mainly consists of coarse-grain and medium-grain sediments (Figure 1) from proximal fans with low compressibility [20], resulting in small settlements.Nevertheless, the sediments over the east, north and south parts of Beijing mainly consist of high compressibility fine-grain sediments [20] with thick quaternary deposits (Figure 1).The aquifers units of these areas have been intensely over-exploited, generating severe land subsidence.
The main subsidence bowls are distributed over the Chaoyang, Changping, Shunyi and Tongzhou Districts.It can be observed in Figure 7 and Figure 8 that the subsidence bowls located in the east, north-east and north parts of Beijing have been gradually connected, indicating the serious situation in these areas.Consequently, this bowl, which appeared for the first time in 1975 [15], has gradually increased in depth and extension until the present.The Xianninghou-Shuangqiao subsidence area, located in the eastern suburbs of Chaoyang District, is in the maximum settlement centre from the East bowl.The maximum subsidence rate in this area for the period 2003-2010 is 110 mm/year (Figure 7).The maximum subsidence from the Northeast bowl is located in Jinzhan-Dongba area with a subsidence rate over 80 mm/year (Figure 7).The Laiguangying subsidence area over Chaoyang District is one of the areas affected by severe land The Xianninghou-Shuangqiao subsidence area, located in the eastern suburbs of Chaoyang District, is in the maximum settlement centre from the East bowl.The maximum subsidence rate in this area for the period 2003-2010 is 110 mm/year (Figure 7).The maximum subsidence from the Northeast bowl is located in Jinzhan-Dongba area with a subsidence rate over 80 mm/year (Figure 7).
The Laiguangying subsidence area over Chaoyang District is one of the areas affected by severe land subsidence with displacement rates over 80 mm/year.As mentioned earlier, the subsidence bowls over the east, northeast and north parts of Beijing have gradually merged, indicating the serious subsidence situation in these districts.In the north part of the Haidian-Changping subsidence area, the maximum displacement rate reaches 60 mm/year.The Dongbalizhuang-Dajiaoting subsidence area (including Sanjianfang, Guanzhuang) over the eastern suburbs of Chaoyang District measured settlements exceeding 50 mm/year (with local values up to 80 mm/year) from 2003 to 2010.The subsidence rate in the Shunyi subsidence area reaches 40 mm/year.
Figure 7 and Figure 8 clearly indicate that the bowls of subsidence in Beijing present a growth trend.The mean subsidence rate around the centres of subsidence in Chaoyang district, located in the eastern part of the city, reaches 80 mm/year, and a maximum subsidence rate higher than 100 mm/year.
Figure 8 illustrates the temporal evolution of the accumulated displacements with an approximately one-year interval, showing that the magnitudes of the annual accumulated displacement vary from one year to another.However, the subsidence rate clearly decreases in 2008.For example, comparing the average accumulated subsidence for a group of SDFP pixels located within the main deformation area (39.898 ˝N, 116.572 ˝E) in different years (December 2007-September 2008 with annual displacement of ´99.86 mm and September 2008-October 2009 with annual displacement of ´119.77mm), a decrease of the annual accumulated subsidence of 16.6% is observed.

Discussion
In this section, different InSAR products including mean displacement rates and displacement time-series are used for the analysis of InSAR-derived surface displacement in the Beijing plain and compared with different triggering and conditioning factors.

Relationships between Land Subsidence and Groundwater Level
The over-exploitation of groundwater in the study region has caused a regional drop of piezometric levels, decreasing pore pressure and increasing effective stresses according to Terzaghi's [36] theory.Consequently, when groundwater is pumped from the aquifer system, both the aquifers and the aquitards that constitute the aquifer system undergo deformation, but to different degrees, as most of the permanent subsidence occurs due to the irreversible compression or consolidation of aquitards [37].Thus, the deformation of an aquifer-system can be inelastic (the aquifer-system undergoes permanent rearrangement resulting in irreversible compaction and only a part of the total deformation is recovered when water level rises) or elastic (when water level rises, the skeleton of the aquifer-system expands, recovering from the previous occurred deformation) depending on whether or not the applied stresses are beyond its previous maximum level (preconsolidation stress threshold) [37].
When land subsidence data are used in conjunction with good well logs and water-level data from adjacent observation wells, the deformation history can provide the basis for stress-strain analysis [38] and inverse modeling that defines the average compressibility and vertical hydraulic conductivity of the aquitards [39].In the study region, hundreds of observation wells are used to monitor the evolution of groundwater levels.Figure 9 shows the relationship between piezometric level and land subsidence for some piezometers from the study region.It is worth noting that the displacements measured by InSAR on the ground surface are the addition of the deformations occurred in the different aquitard and aquifer layers underlying the ground surface and, thus, the stress-strain curves shown in Figure 9 correspond to the whole aquifer system.The analysis of the different stress-strain relationship clearly indicate that the aquifer system presents different behaviour (Figure 9).
If we analyse the general trend of water level and InSAR time series from P1 (see location in Figure 3) we can realise that from June 2003 to February 2007 the water level fell from 36.1 to 26.87 m.a.s.l.(i.e., ´9.23 m) and the accumulated subsidence was up to 38.2 mm.However, after this period, the water level recovered up to 33.5 m.a.s.l. in May 2008 (i.e., it increased +6.6 m, near the original water level) and the ground surface also rose 43.0 mm (the previous deformation was recovered).So, according to the general trend of both time series (water level and displacements), we can see that they both follow a similar trend, i.e., they mainly exhibit an elastic behaviour as when the support provided by the fluid is removed or partially removed (i.e., when water level falls), the support provided by the pore pressure is gradually transferred to the soil skeleton, causing subsidence.Conversely, when the pore pressure recovers (i.e., increase due to water level recovery), the support provided by the skeleton is gradually transferred to the pore pressure, which in turn leads to soil expansion (i.e., uplift) [37].
Remote Sens. 2016, 8, 468 13 of 21 behaviour as when the support provided by the fluid is removed or partially removed (i.e., when water level falls), the support provided by the pore pressure is gradually transferred to the soil skeleton, causing subsidence.Conversely, when the pore pressure recovers (i.e., increase due to water level recovery), the support provided by the skeleton is gradually transferred to the pore pressure, which in turn leads to soil expansion (i.e., uplift) [37].On the other hand, wells P10, P7, P15 exhibit a predominantly inelastic behaviour as the aquifer system presents a continuous subsidence associated with the gradual water level decrease (Figure 9).As the fine sediment layers of the study area are normally consolidated [21] (i.e., the present effective overburden pressure is the maximum pressure that the soil was subjected to in the past), then the geomechanical behaviour of the aquifer-system is characterized by the inelastic storage coefficient (Skv).Skv is a dimensionless parameter which is proportional to the deformability of the aquifer system (i.e., the largest the Skv value, the higher the settlement induced by a water level change).This parameter has been calculated as the inverse slope of the stress-displacement data trend line, according to the graphical methodology proposed by Riley [38], for the observational wells P4 to P18 shown in Figure 9, providing values varying from 4.2 ˆ10 ´3 in P4 to 5.6 ˆ10 ´2 in P18 and confirming that the deformability of the aquifer apparently increases from the north towards the south of the studied area.

Relationships between Land Subsidence and Geological Structures
The distribution and development trend of the land subsidence over the study region are controlled by geological structures, shown in Figure 10a.Tectonic movements cause tectonic units to decline entirely and slowly under the regional stress field and, complementarily, the relative displacement of the hanging walls and the footwalls induces differences in the land subsidence patterns on both sides of the fault [40].There are several active high angle normal faults in the Beijing plain.Huangzhuang-Gaoliying, Shunyi and Nanyuan-Tongzhou faults present a SW-NE direction.The Huangzhuang-Gaoliying fault dips to the SE with an angle of 55 ˝~75 ˝.The Shunyi fault with a SE dipping angle is a boundary fault controlling the Houshayu Quaternary depression, and the Nanyuan-Tongzhou fault dips towards the NW.The Nankou-Sunhe fault is the unique active NW-SE oriented fault in the Beijing basin.The Northwestern segment of the Nankou-Sunhe fault (from Nankou to Baishan) dipping SW controls the Machikou Quaternary depression with a maximum thickness of Quaternary sediments about 600 m, and the Southeastern segment of the Nankou-Sunhe fault (from Baishan to Tongzhou) dipping towards NE controls the depocenter of Dongba, in the southern part of the Shunyi Quaternary depression [27,[41][42][43][44][45][46].In the Beijing basin, Dongbalizhuang-Dajiaoting subsidence area over the eastern subsidence zone is collectively controlled by the Shuyi fault and the Nanyuan-Tongzhou fault.The Laiguangying subsidence area over northeast subsidence district is mainly controlled by the Shuyi fault and the northwest subsidence zone is mostly led by the Nankou-Sunhe and Huangzhuang-Gaoliying faults, for which the main deformation direction is basically parallel to the fracture direction.
The examples of four profiles of the average subsidence rate for Nankou-Sunhe fault can be found in Figure 10b.Profiles AA 1 and BB 1 parallel to the Nankou-Sunhe fault and Profiles CC 1 and DD 1 perpendicular to the Nankou-Sunhe fault with a length of about 10 km are chosen to analyze the spatial characteristics of land subsidence along this fault zone.It can be found that the ground surface displacement along Profile AA 1 located at the north of the Nankou-Sunhe fault shows a relative stable displacement rate of around 25 mm/year.However, the profile of average subsidence rate for BB 1 located in the south of the fault exhibits a significant increasing trend with the maximum subsidence rate exceeding 70 mm/year.The profiles of the average displacement rate for CC 1 and DD 1 perpendicular to the Nankou-Sunhe fault also show a high ground surface displacement gradient along the profile with maximum values up to 77 mm/year, indicating the markedly differential settlement existing on both sides of the fault.

Relationships between Land Subsidence and Compressible Soil Thickness
As is well known, when piezometric levels decrease, the thicker the accumulated compressible soil is, the greater the land subsidence rates are [47,48].Land subsidence over the Beijing plain mainly develops in the sand, gravel and clayey layer of Quaternary system.Among the biggest contributors to land subsidence is clayey layer.Researches have revealed that given the same underground water exploitation situation, the clayey units (aquitards made of fine soils) are two or three times more deformable than those of sand and gravel (coarse soils) [2]. Figure 11a shows the thickness distribution of the clayey layer existing in the first 100 m depth from surface over the Beijing plain area which is superimposed on the land subsidence rate map. Figure 11b shows the statistics of InSAR subsidence rates for the considered soft soil ranges.As can be seen in Figure 11a, the subsidence bowls, which present displacement rates ranging from −110 to −60 mm/year (i.e.,

Relationships between Land Subsidence and Compressible Soil Thickness
As is well known, when piezometric levels decrease, the thicker the accumulated compressible soil is, the greater the land subsidence rates are [47,48].Land subsidence over the Beijing plain mainly develops in the sand, gravel and clayey layer of Quaternary system.Among the biggest contributors to land subsidence is clayey layer.Researches have revealed that given the same underground water exploitation situation, the clayey units (aquitards made of fine soils) are two or three times more deformable than those of sand and gravel (coarse soils) [2]. Figure 11a shows the thickness distribution of the clayey layer existing in the first 100 m depth from surface over the Beijing plain area which is superimposed on the subsidence map. Figure 11b shows the statistics of InSAR subsidence rates for the considered soft soil ranges.As can be seen in Figure 11a, the subsidence bowls, which present displacement rates ranging from ´110 to ´60 mm/year (i.e., about 123329 SDFP pixels within this area from ASAR data), are mostly situated in the zone with a clayey layer thickness of around 50 to 70 m (i.e., about 98190 SDFP pixels within this zone from ASAR data), accounting for 79.6% of the whole amount of subsidence.

Relationships between Land Subsidence and the Aquifer Structure
The structure of the aquifer system in the Beijing plain consists of [20]: (a) coarse-grain sediments with single sandy gravel; (b) median-grain sediments with multiple layers of pebbles and gravels; and (c) fine-grain sediments with multiple sand layers and gravel layers.Figure 12 shows the statistics of InSAR subsidence data for the different types of aquifers systems.It can be seen that the subsidence bowls with displacement values varying from −110 to −60 mm/year are located over predominantly fine-grain sediments, accounting for 8% of the total subsidence rate data.

Relationships between Land Subsidence and the Aquifer Structure
The structure of the aquifer system in the Beijing plain consists of [20]: (a) coarse-grain sediments with single sandy gravel; (b) median-grain sediments with multiple layers of pebbles and gravels; and (c) fine-grain sediments with multiple sand layers and gravel layers.Figure 12 shows the statistics of InSAR data the different types of aquifers systems.It can be seen that the subsidence bowls with displacement values varying from ´110 to ´60 mm/year are located over predominantly fine-grain sediments, accounting for 8% of the total subsidence rate data.

Relationships between Land Subsidence and Main Groundwater Extraction Areas
Groundwater extraction is the main water supply source in the study region [21].Groundwater is extracted from different depths for agricultural, industrial and domestic uses.In this section, we evaluate the effect of groundwater extraction from the Huairou emergency wells drilled in 2003 located along the banks of Huai, Sha and Yanxi Rivers and the well field located along the Chaobai River which were drilled in the 70s (Figure 13).
To evaluate the influence of the pumping wells on land subsidence, different buffers at specified distances were created as shown in Figure 13a.The subsidence rates varied from −35 to −5 mm/year in the buffer areas.The statistics of InSAR subsidence rates against the distance to the wells are plotted in Figure 13b.It shows that the proportion of maximum subsidence rates from −35 to −25 mm/year decrease as the distance to the normal pumping wells increases.A similar relationship is found for the emergency pumping wells (Figure 13b).It indicates that, as expected, land subsidence rates decrease with the distance to pumping wells according to the shape of the cone of depression that develops around active pumping wells.

Relationships between Land Subsidence and Main Groundwater Extraction Areas
Groundwater extraction is the main water supply source in the study region [21].Groundwater is extracted from different depths for agricultural, industrial and domestic uses.In this section, we evaluate the effect of groundwater extraction from the Huairou emergency wells drilled in 2003 located along the banks of Huai, Sha and Yanxi Rivers and the well field located along the Chaobai River which were drilled in the 70s (Figure 13).
To evaluate the influence of the pumping wells on land subsidence, different buffers at specified distances were created as shown in Figure 13a.The subsidence rates varied from ´35 to ´5 mm/year in the buffer The of InSAR subsidence rates against the distance to the wells are plotted in Figure 13b.It shows that the proportion of maximum subsidence rates from ´35 to ´25 mm/year decrease as the distance to the normal pumping wells increases.A similar relationship is found for the emergency pumping wells (Figure 13b).It indicates that, as expected, land subsidence rates decrease with the distance to pumping wells according to the shape of the cone of depression that develops around active pumping wells.

Conclusions
In this study, land subsidence due to over-extraction of groundwater in the Beijing region was investigated with SBAS InSAR using 41 Envisat ASAR and 14 TSX images.The results reveal that the Beijing region has experienced significant ground subsidence from 2003 to 2010 with a maximum accumulative displacement of 790 mm.The SBAS results have been validated by GPS measurements with a mean difference of 2.41 ± 1.84 mm/year and a RMS of 2.94 mm/year, demonstrating that SBAS InSAR can effectively monitor and detect complicated settlements.A high correlation coefficient of 0.92 and a standard deviation of 7.48 mm/year between both TerraSAR-X and Envisat ASAR InSAR-derived subsidence rates was observed, indicating the reliability of InSAR derived land subsidence rate maps.
The spatiotemporal analysis of land subsidence indicates an increasing trend in the rate and extent of land subsidence.The join spatial and/or temporal analysis of InSAR data and conditioning

Conclusions
In this study, land subsidence due to over-extraction of groundwater in the Beijing region was investigated with SBAS InSAR using 41 Envisat ASAR and 14 TSX images.The results reveal that the Beijing region has experienced significant ground subsidence from 2003 to 2010 with a maximum accumulative displacement of 790 mm.The SBAS results have been validated by GPS measurements with a mean difference of 2.41 ˘1.84 mm/year and a RMS of 2.94 mm/year, demonstrating that SBAS InSAR can effectively monitor and detect complicated settlements.A high correlation coefficient of 0.92 and standard of 7.48 mm/year between both TerraSAR-X and Envisat ASAR InSAR-derived subsidence rates was observed, indicating the reliability of InSAR derived land subsidence rate maps.
The spatiotemporal analysis of land subsidence indicates an increasing trend in the rate and extent of land subsidence.The join spatial and/or temporal analysis of InSAR data and conditioning and triggering factors shows that land subsidence is correlated with groundwater levels, active faults, different soft soil ranges and aquifer types.Furthermore, a relationship between land subsidence and the distance to the pumping wells has been found.
The analysis of different relationships shows that the aquifer-system exhibits different behaviour, such as quite elastic behaviour and predominantly inelastic behaviour.Additionally, the analysis of the variation of the inelastic storage coefficient (Skv) of the aquifer system computed from the stress-strain curves shows that the deformability of the aquifer system increases from the north towards the south.
The distribution and development trends of the land subsidence over the study area are obviously controlled by its geological structures.Among the biggest contributors to land subsidence are the clayey layers.The subsidence bowls located in the east, north-east and north are mostly placed over a zone with a clayey layer underneath, with a thickness around 50 to 70 m, and those subsidence bowls are located over predominantly fine-grain sediments.The InSAR data statistics and the distance to the pumping wells show that land subsidence rates are higher near the pumping wells, as expected according to the shape of a well's cone of depression.
To summarize, accurate InSAR-derived data from two different sensors have allowed us to perform a comprehensive spatial and temporal study of the land subsidence in Beijing basin from 2003 to 2011.Complementarily, these data and existing geo-information have been used to perform a study of the role of the main conditioning and triggering factors controlling land subsidence mechanisms in this area, providing useful information for improved development of future models for the prediction of land subsidence affecting this area.

Figure 1 .
Figure 1.(a) Location of the study area; (b) geology of the Beijing plain (China); (c) hydrogeological cross-section of the studied aquifer (adapted from [20]).

Figure 1 .
Figure 1.(a) Location of the study area; (b) geology of the Beijing plain (China); (c) hydrogeological cross-section of the studied aquifer (adapted from [20]).

Figure 2 .
Figure 2. Networks of the small baseline interferograms obtained from (a) Envisat ASAR and (b) TerraSAR-X images used for SBAS InSAR analysis for land subsidence in Beijing.

Figure 2 .
Figure 2. Networks of the small baseline interferograms obtained from (a) Envisat ASAR and (b) TerraSAR-X images used for SBAS InSAR analysis for land subsidence in Beijing.

Figure 4 .
Figure 4. Spatio-temporal evolution of InSAR-derived accumulated LOS displacement in Beijing from Envisat-ASAR images.Note that for simplicity, the displacement patterns of only 8 of the 41 acquisition dates are shown (i.e., one per year).

Figure 4 .
Figure 4. Spatio-temporal evolution of InSAR-derived accumulated LOS displacement in Beijing from Envisat-ASAR images.Note that for simplicity, the displacement patterns of only 8 of the 41 acquisition dates are shown (i.e., one per year).

Figure 5 .
Figure 5.Comparison of GPS and Envisat InSAR results for points G1-G6.See GPS stations location in Figure 3.

Figure 5 .
Figure 5.Comparison of GPS and Envisat InSAR results for points G1-G6.See GPS stations location in Figure 3.

Figure 7 .
Figure 7. (a) Subsidence bowls located in the NE of Beijing; (b) extension of the subsidence bowls.Mapped InSAR results along LOS have been derived from Envisat ASAR images acquired from 2003 to 2010.

Figure 7 .
Figure 7. (a) Subsidence bowls located in the NE of Beijing; (b) extension of the subsidence bowls.Mapped InSAR results along LOS have been derived from Envisat ASAR images acquired from 2003 to 2010.
Remote Sens. 2016, 8, 468 11 of 21 subsidence with displacement rates over 80 mm/year.As mentioned earlier, the subsidence bowls over the east, northeast and north parts of Beijing have gradually merged, indicating the serious subsidence situation in these districts.

Figure 8 .
Figure 8. LOS displacements difference in Beijing obtained from Envisat ASAR images between the acquisition dates selected in Figure 4.

Figure 8 .
Figure 8. LOS displacements difference in Beijing obtained from Envisat ASAR images between the acquisition dates selected in Figure 4.

Figure 9 .
Figure 9.Time series of piezometric level and land subsidence measured by means of InSAR (left); and strain-stress curves, for different areas of Beijing (right).The location of the observation wells is shown in Figure 3.Note that the wells are ordered from north (top) to south (bottom).A1 to A4 refers to different aquifers units described in Section 2.

Figure 9 .
Figure 9.Time series of piezometric level and land subsidence measured by means of InSAR (left); and strain-stress curves, for different areas of Beijing (right).The location of the observation wells is shown in Figure 3.Note that the wells are ordered from north (top) to south (bottom).A1 to A4 refers to different aquifers units described in Section 2.
Remote Sens. 2016, 8, 468 16 of 21 about 123329 SDFP pixels within this area from ASAR data), are mostly situated in the zone with a clayey layer thickness of around 50 to 70 m (i.e., about 98190 SDFP pixels within this zone from ASAR data), accounting for 79.6% of the whole amount of subsidence.

Figure 11 .
Figure 11.(a) Distribution of clayey layers on the study area [28]; (b) Statistics of InSAR-derived land subsidence rate data for the different soft soil ranges.

Figure 11 .
Figure 11.(a) Distribution of clayey layers on the study area [28]; (b) Statistics of InSAR-derived land subsidence rate data for the different soft soil ranges.

Figure 12 .
Figure 12.(a) Distribution and (b) statistics of InSAR derived data for the different types of aquifers.

Figure 12 .
Figure 12.(a) Distribution and (b) statistics of InSAR derived data for the different types of aquifers.

Figure 13 .
Figure 13.(a) Location of the Huairou emergency wells and the well field located along the Chaobai River overlapped to the InSAR subsidence rate map; (b) Statistics of InSAR-derived data against the distance to the Chaobai River (left) and Huairou (right) pumping wells.

Figure 13 .
Figure 13.(a) Location of the Huairou emergency wells and the well field located along the Chaobai River overlapped to the InSAR subsidence rate map; (b) Statistics of InSAR-derived data against the distance to the Chaobai River (left) and Huairou (right) pumping wells.