InSAR Time-Series Analysis of Land Subsidence under Different Land Use Types in the Eastern Beijing Plain , China

In the Beijing plain, the long-term groundwater overexploitation, exploitation, and the utilization of superficial urban space have led to land subsidence. In this study, the spatial–temporal analysis of land subsidence in Beijing was assessed by using the small baseline subset (SBAS) interferometric synthetic aperture radar (InSAR) technique based on 47 TerraSAR-X SAR images from 2010 to 2015. Distinct variations of the land subsidence were found in the study regions. The maximum annual land subsidence rate was 146 mm/year from 2011 to 2015. The comparison between the SBAS InSAR results and the ground leveling measurements showed that the InSAR land subsidence results achieved a precision of 2 mm. In 2013, the maximum displacement reached 132 and 138 mm/year in the Laiguangying and DongbalizhuangDajiaoting area. Our analysis showed that the serious land subsidence mainly occurred in the following land use types: water area and wetland, paddy field, upland soils, vegetable land, and peasant-inhabited land. Our results could provide a useful reference for groundwater exploitation and urban planning.


Introduction
Land subsidence has become a global geological issue which threatens the human living environment.Compared with the traditional leveling and global positioning system (GPS) measurement methods, the interferometric synthetic aperture radar (InSAR) technique can obtain a wide spatial range of surface deformation information with millimeter-scale precision [1,2].Persistent scatterer interferometry (PSI) and small baseline subset (SBAS)-advanced InSAR time series techniques-can minimize the limitations of traditional InSAR, such as the atmospheric errors, the orbit errors, and the terrain errors [3][4][5][6][7][8][9].These two methods (PSI and SBAS) have been successfully applied to study the land subsidence based on the time-series analysis on point targets [10][11][12].
Amelung et al. [13] used InSAR technology to acquire the land subsidence information in Las Vegas; they found that the subsidence rates declined when the ground-water levels rose, and that the spatial destitution of the subsidence is controlled by both faults and clay thickness.Raucoules et al. [14] assessed European Remote Sensing 1 and 2 (ERS-1and ERS-2) images to investigate the land subsidence in the city of Prato (Italy); they suggested that the differential interferometric synthetic aperture radar (DInSAR) technique could play a significant role in the land use planning and natural risk assessment in forthcoming years.Strozzi et al. [15][16][17] acquired the land subsidence in Mexico City by using the DInSAR technique.The DInSAR technique was applicable for operational monitoring of land subsidence in their study.Gong et al. [18] used 15 Envisat ASAR images from 2003 to 2006 to derive the land subsidence information in the Beijing plain based on PSI technology.Zhang et al. [19] combined GPS and InSAR techniques to monitor and analyze land subsidence and ground fissures of Xi'an, Shanxi.Osmanoglu et al. [20] and Bock et al. [21] acquired deformation of Mexico City (2004City ( -2006) ) and Venice (2003)(2004)(2005)(2006)(2007) based on PSI and GPS technique.Chen et al. [22] quantitatively analyzed the relationship between the load construction changes and the land subsidence, and the study results showed a positive correlation between the load density and the homogeneity of subsidence-especially in the areas which have a high sedimentation rate.Regional land subsidence is a frequent environmental problem which needs to be examined in long-term systematic research.For example, Zhang et al. [23] detected the hydraulic and mechanical properties of hydrogeologic units at different depths by using stress-strain analyses and oedometer tests based on the data from leveling surveys, borehole extensometers, and multilayer monitoring of groundwater levels.These results can help to improve the assessment and the prediction of land subsidence.SBAS InSAR has been applied in Xi'an, Shaanxi to study the relationships between land subsidence and fissures [24].They detected four main land subsidence zones in the Xi'an Hi-tech Zone, Xi'an Qujiang New Zone, Chanba River Economic Zone, and Hujiamiao, with an average land subsidence rate of about 50 mm/year during 2005 to 2012.Good correlations of the ground water-level decline and the land subsidence were found in their studies.Meanwhile, Zhu et al. [25] applied PSI on the ENVISAT ASAR images from 2003 to 2010 to investigate the land subsidence in the northern Beijing plain.They reported that the maximum deformation rate reached 52 mm/year, with a cumulative maximum sinking 342 mm in Houshayu, located in the southwestern part of the Beijing plain.Luo et al. [26] analyzed the historical and current data of the land subsidence in the Beijing plain to reveal that the land subsidence center moved from the west to the east and also moved from the north to the south with the increased amount of land subsidence.
Beijing, in Northern China, is the capital city of China, and the world's third most populous city.Almost two-thirds of the urban water supply comes from groundwater [27].In recent years, the water consumption has sharply increased due to the rapid expansion of the urban population [28].Groundwater extraction was measured as 2.6 × 10 9 m 3 /year, with an overexploitation of approximately 1 × 10 6 m 3 /year [21].This long-term overexploitation of groundwater caused a substantial decline in head water, and consequently induced land subsidence.By the end of 2010, the land subsidence area reached 4.2 × 10 3 km 2 , and 66% of the Beijing plain has been affected by land subsidence (>50 mm), with a maximum sinking of 1233 mm [25].
Previous studies have generally focused on time series InSAR technology to obtain deformation information combined with groundwater numerical simulation to reveal the mechanism of land subsidence.However, the combination of the InSAR technology and optical remote sensing technology to investigate the characteristics of the land subsidence in the different land use types has been rarely studied.In this manuscript, we applied the Stanford method for persistent scatterer (StaMPS) SBAS technique [10] to obtain long-term time-series land displacement information in the eastern Beijing plain based on 47 TerraSAR-X (X-band) images from 2010 to 2015.The geographical and geological characteristics of the Beijing plain have been assessed to examine the forcing factors of the land subsidence.Our analysis was based on high-accuracy deformation results, which were evaluated by the ground truth leveling data.

Study Area
Beijing is the capital city of China, and is the political, cultural, educational, and international center with over 20 million residences.In Beijing, the major urban water consumption comes from groundwater.Beijing is bounded by the Taihang Mountains to the west and the Jundu Mountains to the northeast, which is located in the longitudes 115 • 25 E-117 • 30 E and latitudes 39 • 26 N-41 • 03 N. The city is adjacent to Hebei Province and Tianjin City in the south.Beijing has 16 districts and 2 counties.The total area of Beijing is about 1.641 × 10 4 km 2 with an approximately 6.39 × 10 3 km 2 plain area.Beijing was built on the alluvial fan of the YongDing River, which originates from Shanxi Province.
The study area (Figure 1) is a typical warm temperate semi-arid continental monsoon climate where the annual mean temperature is about 10-12 • C. The precipitation in Beijing is unevenly distributed throughout the year; the rainfall is concentrated during June-September, accounting for 60-85% of the whole year's rainfall.

Study Area
Beijing is the capital city of China, and is the political, cultural, educational, and international center with over 20 million residences.In Beijing, the major urban water consumption comes from groundwater.Beijing is bounded by the Taihang Mountains to the west and the Jundu Mountains to the northeast, which is located in the longitudes 115°25′E-117°30′E and latitudes 39°26′N-41°03′N.The city is adjacent to Hebei Province and Tianjin City in the south.Beijing has 16 districts and 2 counties.The total area of Beijing is about 1.641 × 10 4 km 2 with an approximately 6.39 × 10 3 km 2 plain area.Beijing was built on the alluvial fan of the YongDing River, which originates from Shanxi Province.
The study area (Figure 1) is a typical warm temperate semi-arid continental monsoon climate where the annual mean temperature is about 10-12 °C.The precipitation in Beijing is unevenly distributed throughout the year; the rainfall is concentrated during June-September, accounting for 60-85% of the whole year's rainfall.The quaternary aquifers of the Beijing plain have three aquifer groups [29,30].The first aquifer group is the Holocene and late Pleistocene strata, the depths of which are less than 100 m to the bottom.This aquifer has two divided aquifers.One is mainly distributed in the lower part of the alluvial fan with a depth of about 25 m to the bottom.The lithologies include fine sand, silt, silty sand, and sandy clay.Another one has a depth of 80-100 m to the bottom.The bedrock forms the bottom boundary in some areas.The second aquifer group is the middle Pleistocene strata with depths to the roof of 80 to 100 m and depths to the bottom of around 300 m.Both the lithologies and the middle Pleistocene strata are multilayer structure and include multiple types of gravel, sand, and clay soil.The groundwater here is also used as agricultural, domestic, and industrial water.The third aquifer group is the early Pleistocene strata, and the depth to the top is more than 300 m.The groundwater exploitation is less intensive in this layer because its compressibility is low.
From the 1970s to the 2000s, land subsidence in the Beijing plain area developed significantly.By the end of 2012, more than 2/3 of the area in the Beijing plain was affected by land subsidence.A total of five land subsidence bowls were formed, which were named as the DongbalizhuangDajiaoting, Laiguangying, ChangpingShahe-Ba Xianzhuang, DaxingYufa-Lixian, and Shunyi-Ping Gezhuang.In this paper, we selected the DongbalizhuangDajiaoting and Laiguangying land subsidence areas The quaternary aquifers of the Beijing plain have three aquifer groups [29,30].The first aquifer group is the Holocene and late Pleistocene strata, the depths of which are less than 100 m to the bottom.This aquifer has two divided aquifers.One is mainly distributed in the lower part of the alluvial fan with a depth of about 25 m to the bottom.The lithologies include fine sand, silt, silty sand, and sandy clay.Another one has a depth of 80-100 m to the bottom.The bedrock forms the bottom boundary in some areas.The second aquifer group is the middle Pleistocene strata with depths to the roof of 80 to 100 m and depths to the bottom of around 300 m.Both the lithologies and the middle Pleistocene strata are multilayer structure and include multiple types of gravel, sand, and clay soil.The groundwater here is also used as agricultural, domestic, and industrial water.The third aquifer group is the early Pleistocene strata, and the depth to the top is more than 300 m.The groundwater exploitation is less intensive in this layer because its compressibility is low.
From the 1970s to the 2000s, land subsidence in the Beijing plain area developed significantly.By the end of 2012, more than 2/3 of the area in the Beijing plain was affected by land subsidence.A total of five land subsidence bowls were formed, which were named as the DongbalizhuangDajiaoting, Laiguangying, ChangpingShahe-Ba Xianzhuang, DaxingYufa-Lixian, and Shunyi-Ping Gezhuang.In this paper, we selected the DongbalizhuangDajiaoting and Laiguangying land subsidence areas located in eastern of Beijing as the study area.The area is the future regional planning and construction center of Beijing, including the Central Business District (CBD), CBD eastward expansion, and the Tongzhou New City planning area.Moreover, a number of rail transit lines have been built or are under construction in the area.

Data Source
The TerraSAR-X sensor was launched on 21 June 2007.It has a revisiting time of 11 days and operates in the X-band of HH, VV, HV, and VH polarization modes with a 31 mm wavelength [31].In this paper, we used 48 SAR images from ascending track, acquired by the TerraSAR-X satellite between April 2010 and November 2015.The tracking number is 157 and the incidence angle is 33.23 • with the look direction of NE.The spatial and temporal baselines ranged from 596 m to 10 m and from 33 days to 924 days, respectively.Doris software (Delft object-oriented radar interferometric software) was applied to perform the interferometric process, and included the following steps: master image selection, topographic phase removal, and interferogram generation.In the study, we employed the SBAS (Stanford method for persistent scatterers, StaMPS) technology to acquire land subsidence information.The topographic phase was removed by using the 90 m resolution Shuttle Radar Topography Mission (SRTM) DEM data from the United States Geological Survey (USGS).The spatial adaptive filter and temporal filter were applied to separate the atmospheric phase, noise phase, and nonlinear deformation phase.Then, two-dimensional spatial high-pass filtering was used to remove the orbital linear ramp and other long-wavelength noise from the differential interferograms.A total of 69 interferogram pairs were built with spatial and temporal baselines smaller than 300 m and 300 days, respectively (Figure 2).
Remote Sens. 2017, 9, 380 4 of 16 located in eastern of Beijing as the study area.The area is the future regional planning and construction center of Beijing, including the Central Business District (CBD), CBD eastward expansion, and the Tongzhou New City planning area.Moreover, a number of rail transit lines have been built or are under construction in the area.

Data Source
The TerraSAR-X sensor was launched on 21 June 2007.It has a revisiting time of 11 days and operates in the X-band of HH, VV, HV, and VH polarization modes with a 31 mm wavelength [31].In this paper, we used 48 SAR images from ascending track, acquired by the TerraSAR-X satellite between April 2010 and November 2015.The tracking number is 157 and the incidence angle is 33.23° with the look direction of NE.The spatial and temporal baselines ranged from 596 m to 10 m and from 33 days to 924 days, respectively.Doris software (Delft object-oriented radar interferometric software) was applied to perform the interferometric process, and included the following steps: master image selection, topographic phase removal, and interferogram generation.In the study, we employed the SBAS (Stanford method for persistent scatterers, StaMPS) technology to acquire land subsidence information.The topographic phase was removed by using the 90 m resolution Shuttle Radar Topography Mission (SRTM) DEM data from the United States Geological Survey (USGS).The spatial adaptive filter and temporal filter were applied to separate the atmospheric phase, noise phase, and nonlinear deformation phase.Then, two-dimensional spatial high-pass filtering was used to remove the orbital linear ramp and other long-wavelength noise from the differential interferograms.A total of 69 interferogram pairs were built with spatial and temporal baselines smaller than 300 m and 300 days, respectively (Figure 2).

Small Baseline Interferometry
Small baseline interferometry (SBAS) was originally proposed by Kampes et al. [8].The method aims to minimize the separation both in the time period and the Doppler frequency range of acquisition pairs to increase the correlation of the interferograms' formation.Hopper et al. [24] developed the single-look slowly-decorrelating filtered phase (SDFP) for the short time interval application.Thus, the isolated SDFP pixels surrounded by associated pixels can be identified for processing [7,9].
In our studies, we first created the interferograms in our study sites.Interferogram pairs with perpendicular, temporal, and Doppler separations less than the threshold values were selected for processing.All the available data were selected for processing to ensure a network of image pairs without isolated clusters (Figure 2).In this paper, the maximum perpendicular absolute baseline

Small Baseline Interferometry
Small baseline interferometry (SBAS) was originally proposed by Kampes et al. [8].The method aims to minimize the separation both in the time period and the Doppler frequency range of acquisition pairs to increase the correlation of the interferograms' formation.Hopper et al. [24] developed the single-look slowly-decorrelating filtered phase (SDFP) for the short time interval application.Thus, the isolated SDFP pixels surrounded by associated pixels can be identified for processing [7,9].
In our studies, we first created the interferograms in our study sites.Interferogram pairs with perpendicular, temporal, and Doppler separations less than the threshold values were selected for processing.All the available data were selected for processing to ensure a network of image pairs without isolated clusters (Figure 2).In this paper, the maximum perpendicular absolute baseline difference (300 m) and maximum absolute temporal difference (300 days) were applied according to the data availability and the environmental settings in the study site.A total of 69 interferograms were generated to examine the forcing factors of land subsidence in the Beijing plain.Subsequently, we selected all the SDFP pixels which showed little loss of correlation in the short time period.These processes greatly benefited the SBAS operation at the highest spatial resolution, and enabled the identification of isolated SDFP pixels which were surrounded by the decorrelated pixels.The temporal coherence threshold for the pixels characteristic was 0.7, which was applied in most InSAR applications.Then, we applied the least-squares method in the phase unwrapping.The spatial filters were also assessed to reduce the phase noise.In addition, the spatial-correlated look angle (SCLA) error phase and atmosphere and orbit error (AOE) phase were simultaneously estimated after phase unwrapping.Finally, we derived the time series displacement information with high spatial resolution in Beijing plain.

SBAS-Derived Land Subsidence Map
A deformations map in the Beijing plain was generated to illustrate the spatial distribution of the ground surface subsidence.Figure 3 shows the result of the annual average deformation rate from 2011 to 2015.The results indicate that the spatial distribution of the land subsidence varies widely.The major land subsidence was concentrated in the eastern Beijing regions.In the period of 2011-2015, we detected a large amount of SDFP pixels (122,272) in the land subsidence spatial-temporal dynamics analysis.The maximum land subsidence rate reached 146 mm/year, demonstrating that significant land subsidence was captured.We selected the area located in the east of the Beijing plain-where the land subsidence is more serious (the region of the red line in Figure 3)-for the analysis of the land subsidence.
Remote Sens. 2017, 9, 380 5 of 16 difference (300 m) and maximum absolute temporal difference (300 days) were applied according to the data availability and the environmental settings in the study site.A total of 69 interferograms were generated to examine the forcing factors of land subsidence in the Beijing plain.Subsequently, we selected all the SDFP pixels which showed little loss of correlation in the short time period.These processes greatly benefited the SBAS operation at the highest spatial resolution, and enabled the identification of isolated SDFP pixels which were surrounded by the decorrelated pixels.The temporal coherence threshold for the pixels characteristic was 0.7, which was applied in most InSAR applications.Then, we applied the least-squares method in the phase unwrapping.The spatial filters were also assessed to reduce the phase noise.In addition, the spatial-correlated look angle (SCLA) error phase and atmosphere and orbit error (AOE) phase were simultaneously estimated after phase unwrapping.Finally, we derived the time series displacement information with high spatial resolution in Beijing plain.

SBAS-Derived Land Subsidence Map
A deformations map in the Beijing plain was generated to illustrate the spatial distribution of the ground surface subsidence.Figure 3 shows the result of the annual average deformation rate from 2011 to 2015.The results indicate that the spatial distribution of the land subsidence varies widely.The major land subsidence was concentrated in the eastern Beijing regions.In the period of 2011-2015, we detected a large amount of SDFP pixels (122,272) in the land subsidence spatial-temporal dynamics analysis.The maximum land subsidence rate reached 146 mm/year, demonstrating that significant land subsidence was captured.We selected the area located in the east of the Beijing plain-where the land subsidence is more serious (the region of the red line in Figure 3)-for the analysis of the land subsidence.

SBAS Accuracy Assessment
We compared the satellite-derived land subsidence results with the ground truth leveling data to evaluate the accuracy of SBAS-derived deformation from TerraSAR-X data [32].A total of six leveling benchmarks measuring precise annual deformation data from 2011 to 2013 were assessed.The result was showed in Table 1.The location of these benchmarks and reference points are evenly distributed in our study regions (Figure 1).Since the deformation information is along the SAR's line-of-sight (LOS) direction, we projected the subsidence results along with a vertical direction with respect to the corresponding incidence angles [27].SDFP points surrounding within 100 m of the benchmarks were collected.Then, the average displacement values around the benchmarks from 2011 to 2013 were calculated for the evaluation.The accuracy evaluation of our results is presented in Figure 4.The SBAS-derived results matched with the ground leveling measurements with high accuracy.The error ranged from 2 mm to 12 mm.The correlation coefficient (95% confidence level) was 0.997, and the R-squared was 0.994.The mean absolute deviation was 6 mm/year, and the standard deviation was ±3 mm/year.Evaluation of the results suggested that the land subsidence in the Beijing plain was successfully derived with high accuracy.This information provides reliable sources for the analysis of the impact factors on the land subsidence in our study.

SBAS Accuracy Assessment
We compared the satellite-derived land subsidence results with the ground truth leveling data to evaluate the accuracy of SBAS-derived deformation from TerraSAR-X data [32].A total of six leveling benchmarks measuring precise annual deformation data from 2011 to 2013 were assessed.The result was showed in Table 1.The location of these benchmarks and reference points are evenly distributed in our study regions (Figure 1).Since the deformation information is along the SAR's line-of-sight (LOS) direction, we projected the subsidence results along with a vertical direction with respect to the corresponding incidence angles [27].SDFP points surrounding within 100 m of the benchmarks were collected.Then, the average displacement values around the benchmarks from 2011 to 2013 were calculated for the evaluation.The accuracy evaluation of our results is presented in Figure 4.The SBAS-derived results matched with the ground leveling measurements with high accuracy.The error ranged from 2 mm to 12 mm.The correlation coefficient (95% confidence level) was 0.997, and the R-squared was 0.994.The mean absolute deviation was 6 mm/year, and the standard deviation was ±3 mm/year.Evaluation of the results suggested that the land subsidence in the Beijing plain was successfully derived with high accuracy.This information provides reliable sources for the analysis of the impact factors on the land subsidence in our study.Comparisons of SBAS-derived land subsidence rates and leveling measurement rates from 2011 to 2013.The location of benchmarks is marked as "BJ021", "BJ031", "BJ054", "BJ070", "BJ076" and "BJ135" in Figure 1.Comparisons of SBAS-derived land subsidence rates and leveling measurement rates from 2011 to 2013.The location of benchmarks is marked as "BJ021", "BJ031", "BJ054", "BJ070", "BJ076" and "BJ135" in Figure 1.

Time Series Land Subsidence in Eastern Beijing
We focused on the two severe land subsidence areas of the Laiguangying and DongbaluzhuangDajiaoting in Eastern Beijing for further analysis.Figure 5 shows the land deformation from 2011 to 2015 by in these two regions.The land subsidence was spatially concentrated more in the east regions than in the west regions.In the Laiguangying area, the maximum land subsidence rates were 123, 115, 132, 129, and 119 mm/year from 2011 to 2015.In the DongbalizhuangDajiaoting area, the maximum rates were 132, 116, 138, 124, and 115 mm/year from 2011 to 2015.The various land use types in the study area had impacts on the uneven development of the land subsidence.Meanwhile, there were similar temporal patterns in the two land subsidence areas (Figure 6).Indeed, the land subsidence rate in the Laiguangying and DongbalizhuangDajiaoting regions decreased from 2011 to 2012 and from 2013 to 2015.However, the land subsidence rate increased from 2012 to 2013.The research area is located in a major area which suffered groundwater exploitation.However, after the implementation of the South-North Water Transfer Project, the reduction in the groundwater exploitation may have slowed down the land subsidence rate from 2013 to 2015.The water transfer project was first implemented on 27 December 2002.From 2008, Beijing city received water from the South-North Water Transfer Project.The implementation of the South-North Water Transfer Project greatly released the water consumption pressures in the Beijing city.The project also significantly contributed to the land subsidence management in the Beijing plain.

Time Series Land Subsidence in Eastern Beijing
We focused on the two severe land subsidence areas of the Laiguangying and DongbaluzhuangDajiaoting in Eastern Beijing for further analysis.Figure 5 shows the land deformation from 2011 to 2015 by in these two regions.The land subsidence was spatially concentrated more in the east regions than in the west regions.In the Laiguangying area, the maximum land subsidence rates were 123, 115, 132, 129, and 119 mm/year from 2011 to 2015.In the DongbalizhuangDajiaoting area, the maximum rates were 132, 116, 138, 124, and 115 mm/year from 2011 to 2015.The various land use types in the study area had impacts on the uneven development of the land subsidence.Meanwhile, there were similar temporal patterns in the two land subsidence areas (Figure 6).Indeed, the land subsidence rate in the Laiguangying and DongbalizhuangDajiaoting regions decreased from 2011 to 2012 and from 2013 to 2015.However, the land subsidence rate increased from 2012 to 2013.The research area is located in a major area which suffered groundwater exploitation.However, after the implementation of the South-North Water Transfer Project, the reduction in the groundwater exploitation may have slowed down the land subsidence rate from 2013 to 2015.The water transfer project was first implemented on 27 December 2002.From 2008, Beijing city received water from the South-North Water Transfer Project.The implementation of the South-North Water Transfer Project greatly released the water consumption pressures in the Beijing city.The project also significantly contributed to the land subsidence management in the Beijing plain.

Land Subsidence Characteristics in Different Land Use Types
Correspondence analysis was applied to detect the spatial characteristics of land subsidence in the different land use areas.Correspondence analysis is a multivariate statistical approach for exploring the cross-tabulated categorical data, such as social meaning variables.The technique generates a contingency table composed of categorical variables, and it represents the table in a two-dimensional graph for representing and interpreting the relationships between categorical variables [25].We applied correspondence analysis to investigate the relationship between the land use types and land subsidence based on SPSS 22 [33].The land use classification data (2013) were provided by the Beijing Municipal Environmental Protection Bureau.The land use classification data are presented in Table 2.The land subsidence rate classification is presented in Table 3 (the minus signs represent subsidence).We first calculated the number of SDFP pixels in the different land subsidence classes in the different land use types.Then, we calculated the ratio of SDFP pixels of each deformation region in different land use types.This ratio was used as the frequency in SPSS for correspondence analysis.The frequency here represents the probability of the occurrence of land subsidence rate in order to judge the land subsidence characteristics in different land use type areas.The frequencies are shown in Table 4.

Land Subsidence Characteristics in Different Land Use Types
Correspondence analysis was applied to detect the spatial characteristics of land subsidence in the different land use areas.Correspondence analysis is a multivariate statistical approach for exploring the cross-tabulated categorical data, such as social meaning variables.The technique generates a contingency table composed of categorical variables, and it represents the table in a two-dimensional graph for representing and interpreting the relationships between categorical variables [25].We applied correspondence analysis to investigate the relationship between the land use types and land subsidence based on SPSS 22 [33].The land use classification data (2013) were provided by the Beijing Municipal Environmental Protection Bureau.The land use classification data are presented in Table 2.The land subsidence rate classification is presented in Table 3 (the minus signs represent subsidence).We first calculated the number of SDFP pixels in the different land subsidence classes in the different land use types.Then, we calculated the ratio of SDFP pixels of each deformation region in different land use types.This ratio was used as the frequency in SPSS for correspondence analysis.The frequency here represents the probability of the occurrence of land subsidence rate in order to judge the land subsidence characteristics in different land use type areas.The frequencies are shown in Table 4.The results are shown in Figure 7.The steps to understanding Figure 7 are as followings: the X-axis and Y-axis are illustrated beginning at the zero point in Figure 7a-g.Then, we built the vectors and extended lines from the central point to each type of land subsidence class.Figure 7a-g represents the first land subsidence class to the seventh land subsidence class.The perpendicular lines from each land use type to the vectors and extended lines were also built.The distance of the vertical point to the positive vector represented the likelihood of the land subsidence in certain land use type (closer means higher chances).
We found that the correlations were the highest between the first and fourth land subsidence rates and the fifth land use type (Figure 7a-g).Meanwhile, the second and sixth land subsidence rate categories were most correlated with the second land use type.The third land subsidence rate was related to the seventh land use type.The fifth land subsidence rate often occurs in the first land use type, and the seventh land subsidence shows the strongest correlation with the sixth land use type.These sequences of the relationship between land use types and land subsidence rates are shown in Table 5.The first column of land use types has the strongest correspondence with land subsidence, which was water area and wetland, vegetable land, paddy field, upland soils, and peasant-inhabited land.We have ignored the urban construction land use type because it corresponds with a relatively small land subsidence rate.
Table 5.The sequence of the relationship between land use types and land subsidence rate.From left to right, the correlation between land subsidence and land use types decreases gradually.WA (water area and wetland); L (large industrial zones airport and other special land); U (urban construction land); PE (peasant-inhabited land); WS (wasteland); PA (paddy field and upland soils).

Relationships between Land Subsidence and Faults
The spatial distribution of faults affects the distribution and development of the land subsidence in the Beijing plain.There are three active faults across the study area, which are the Liangxiang-Qianmen-Shunyi, Nanyuan-Tongxian, and Nankou-Sunhe faults (Figure 8).The direction of the Liangxiang-Qianmen-Shunyi and Nanyuan-Tongxian faults is SW to NE.The direction of the Nankou-Sunhe fault is NW to SE.We found that both sides of the Liangxiang-Qianmen-Shunyi fault are in the Laiguangying land subsidence area, and the Nanyuan-Tongxian fault is in the DongbalizhuangDajiaoting land subsidence area, with the clear distribution of the land subsidence.The two profiles of the Liangxiang-Qianmen-Shunyiand Nankou-Sunhe (profile A-A1) and Nanyuan-Tongxian (profile B-B1) faults were conducted to observe how the deformation varies on both sides of the faults.The Figure 9 illustrated that some peak values and valley values appeared on both sides of the faults.On the east side of the Liangxiang-Qianmen-Shunyi fault (F1), the land subsidence rate increased significantly with the maximum subsidence rate over 90 mm/year.However, the land subsidence rate on both sides of the Nankou-Sunhe (F2) and Nanyuan-Tongxian (F3) faults shows distinct different subsidence with the maximum values of 98 and 106 mm/year.We can conclude from the results that the geological structures of the Beijing plain have remarkable control of the spatial distribution of land subsidence.

Relationships between Land Subsidence and Faults
The spatial distribution of faults affects the distribution and development of the land subsidence in the Beijing plain.There are three active faults across the study area, which are the Liangxiang-Qianmen-Shunyi, Nanyuan-Tongxian, and Nankou-Sunhe faults (Figure 8).The direction of the Liangxiang-Qianmen-Shunyi and Nanyuan-Tongxian faults is SW to NE.The direction of the Nankou-Sunhe fault is NW to SE.We found that both sides of the Liangxiang-Qianmen-Shunyi fault are in the Laiguangying land subsidence area, and the Nanyuan-Tongxian fault is in the DongbalizhuangDajiaoting land subsidence area, with the clear distribution of the land subsidence.The two profiles of the Liangxiang-Qianmen-Shunyiand Nankou-Sunhe (profile A-A1) and Nanyuan-Tongxian (profile B-B1) faults were conducted to observe how the deformation varies on both sides of the faults.The Figure 9 illustrated that some peak values and valley values appeared on both sides of the faults.On the east side of the Liangxiang-Qianmen-Shunyi fault (F1), the land subsidence rate increased significantly with the maximum subsidence rate over 90 mm/year.However, the land subsidence rate on both sides of the Nankou-Sunhe (F2) and Nanyuan-Tongxian (F3) faults shows distinct different subsidence with the maximum values of 98 and 106 mm/year.We can conclude from the results that the geological structures of the Beijing plain have remarkable control of the spatial distribution of land subsidence.

Relationships between Land Subsidence and Ground Water Level in Different Land Use Types
Our study indicated that the land subsidence mainly occurs in four land use types: water area and wetland, vegetable land, paddy field and upland soils, and peasant-inhabited land.The spatial distribution of land subsidence versus the land use types is shown in Figure 8.We analyzed the area of each land use type and the characteristics of land subsidence.We found that the largest area of land use type associated with land subsidence was the water area and wetland, with the area of 46.69 km 2 and 522 SDFP pixels.The land use type of vegetable has the smallest area, with 3.41 km 2 and 275 SDFP pixels.From 2011 to 2015, the maximum land subsidence rate in these four land use types was all above 120 mm/year (121, 125, 126, and 127 mm/year) (Table 6).Similarly, the mean land subsidence rate in these four land use types from 2011 to 2015 all exceeded 40 mm/year.Previous studies have proven that the overexploitation of the ground water was the major cause of land subsidence in the Beijing plain.Because detailed information about ground water pumping in recent years is not publically accessible, we obtain the groundwater level changes from 2011 and 2013.We generated thirty random SDFP pixels from four land use types for the investigation.Figure 10 shows a good relationship between the groundwater level and land subsidence.The decline of the groundwater level generally led to the increase of the ground subsidence.Furthermore, we investigated the time series of groundwater level and land subsidence in Figure 10 (right).Additionally, groundwater level and land subsidence have a strong correlation in the time series.There is a relatively stronger correlation between the groundwater level and the land subsidence in the vegetable land (Figure 10

Relationships between Land Subsidence and Ground Water Level in Different Land Use Types
Our study indicated that the land subsidence mainly occurs in four land use types: water area and wetland, vegetable land, paddy field and upland soils, and peasant-inhabited land.The spatial distribution of land subsidence versus the land use types is shown in Figure 8.We analyzed the area of each land use type and the characteristics of land subsidence.We found that the largest area of land use type associated with land subsidence was the water area and wetland, with the area of 46.69 km 2 and 522 SDFP pixels.The land use type of vegetable has the smallest area, with 3.41 km 2 and 275 SDFP pixels.From 2011 to 2015, the maximum land subsidence rate in these four land use types was all above 120 mm/year (121, 125, 126, and 127 mm/year) (Table 6).Similarly, the mean land subsidence rate in these four land use types from 2011 to 2015 all exceeded 40 mm/year.Previous studies have proven that the overexploitation of the ground water was the major cause of land subsidence in the Beijing plain.Because detailed information about ground water pumping in recent years is not publically accessible, we obtain the groundwater level changes from 2011 and 2013.We generated thirty random SDFP pixels from four land use types for the investigation.Figure 10 shows a good relationship between the groundwater level and land subsidence.The decline of the groundwater level generally led to the increase of the ground subsidence.Furthermore, we investigated the time series of groundwater level and land subsidence in Figure 10 (right).Additionally, groundwater level and land subsidence have a strong correlation in the time series.There is a relatively stronger correlation between the groundwater level and the land subsidence in the vegetable land (Figure 10  (a2-d2) is the time series relationship between the groundwater level and land subsidence.

Relationships between Land Subsidence and Compressible Deposits
The uneven distribution of the land subsidence was also strongly affected by the presence and thickness of compressible layers [28,34].Thus, we correlated the areas of different thickness of the compressible deposits (50 to 210 m) with the average land subsidence from 2011 to 2015 (Figure 11).As seen from Figure 11 (left), the maximum land subsidence rate occurred in the compressible thickness of around 150 to 200 m.We found that the average subsidence rate decreased from the compressible thickness of 50 m to 260 m.Consequently, we proposed that the land subsidence is likely to occur in the region where compressible deposits exceed 100 m.The reason would be that

Relationships between Land Subsidence and Compressible Deposits
The uneven distribution of the land subsidence was also strongly affected by the presence and thickness of compressible layers [28,34].Thus, we correlated the areas of different thickness of the compressible deposits (50 to 210 m) with the average land subsidence from 2011 to 2015 (Figure 11).As seen from Figure 11 (left), the maximum land subsidence rate occurred in the compressible thickness of around 150 to 200 m.We found that the average subsidence rate decreased from the compressible thickness of 50 m to 260 m.Consequently, we proposed that the land subsidence is likely to occur in the region where compressible deposits exceed 100 m.The reason would be that the stratum of the compressible thickness between 100 to 260 is the middle Pleistocene, and the corresponding aquifer is the main mining layer of groundwater in the Beijing plain [35].
Remote Sens. 2017, 9, 380 14 of 16 the stratum of the compressible thickness between 100 to 260 is the middle Pleistocene, and the corresponding aquifer is the main mining layer of groundwater in the Beijing plain [35].

Conclusions
Using a regional SBAS time series survey with TerraSAR-X InSAR imagery for 2010-2015, we detected the line-of-sight (LOS) ground surface deformation in the Beijing plain.Over the period of 2010-2015, we detected 122,272 SDFP pixels.The maximum land subsidence rate reached 146 mm/year in the study regions.The comparisons of the InSAR results with ground truth leveling measurements showed that the best precision of our results was 2 mm.Meanwhile, the correlation coefficient (the confidence level at 95%) was 0.997 and the R-squared was 0.994.We selected the Laiguangying and DongbalizhuangDajiaoting areas in the Eastern Beijing plain to study the land subsidence.The land subsidence varied greatly in the study area from 2011 to 2015.The maximum land subsidence rate reached 127 mm/year.The minimum land subsidence rate was 16 mm/year.
We analyzed the spatial characteristics of land subsidence in different land use types.A more serious ground subsidence rate ranging from −127 to 20 mm/year was more likely to occur in the following order: water area and wetland, paddy field, upland soils, vegetable land and peasant-inhabited.The different land use regions generally related to the different ground water consumption levels.
We have discussed the relationships between the multiple forcing factors (faults, groundwater level, and compressible thickness) on the land subsidence rates.The geological structure affects the land subsidence, in that the ground subsidence changes obviously across the faults.The land subsidence increased with the decrease of the groundwater levels.In addition, the land subsidence is also affected by compressible thickness.The subsidence is more likely to occur in the areas where the compressible thickness is above 100 m.
In this paper, we applied the X-band (TerraSAR-X) SAR data to investigate the land subsidence caused by human activities.In future research, we will use the X-band SAR data to distinguish natural and anthropogenic subsidence [36].Furthermore, we will use the multi-band SAR data to analyze the spatial distribution characteristics of land subsidence [37].

Figure 1 .
Figure 1.The geographic location of the study area.The red star is the location of reference points for interferometric synthetic aperture radar (InSAR) and leveling measurements.The red pushpins (BJ021, BJ031, BJ054, BJ070, BJ076, BJ135) in (a) and (b) indicate the locations of the leveling benchmarks.In (a), the black box represents the TerraSAR-X data spatial coverage.

Figure 1 .
Figure 1.The geographic location of the study area.The red star is the location of reference points for interferometric synthetic aperture radar (InSAR) and leveling measurements.The red pushpins (BJ021, BJ031, BJ054, BJ070, BJ076, BJ135) in (a) and (b) indicate the locations of the leveling benchmarks.In (a), the black box represents the TerraSAR-X data spatial coverage.

Figure 2 .
Figure 2. A network of interferogram pairs obtained from TerraSAR-X used in small baseline subset (SBAS) InSAR.

Figure 2 .
Figure 2. A network of interferogram pairs obtained from TerraSAR-X used in small baseline subset (SBAS) InSAR.

Figure 3 .
Figure 3. Mean land subsidence rate based on the TerraSAR-X data from 2011 to 2015.The red lines are the study area boundary.L and D are the Laiguangying and DongbalizhuangDajiaoting land subsidence areas, respectively.

Figure 3 .
Figure 3. Mean land subsidence rate based on the TerraSAR-X data from 2011 to 2015.The red lines are the study area boundary.L and D are the Laiguangying and DongbalizhuangDajiaoting land subsidence areas, respectively.

Figure 5 .
Figure 5. Displacement information between 2011 and 2015 in the Laiguangying and DongbalizhuangDajiaoting land subsidence areas measured by the SBAS technique using TerraSAR-X data.(a-e) is annual deformation from 2011 to 2015, respectively.

Figure 5 .
Figure 5. Displacement information between 2011 and 2015 in the Laiguangying and DongbalizhuangDajiaoting land subsidence areas measured by the SBAS technique using TerraSAR-X data.(a-e) is annual deformation from 2011 to 2015, respectively.

Figure 6 .
Figure 6.Time series of the land displacement in the study area.L indicates the Laiguangying land subsidence area, and D indicates the DongbalizhuangDajiaoting land subsidence area.

Figure 6 .
Figure 6.Time series of the land displacement in the study area.L indicates the Laiguangying land subsidence area, and D indicates the DongbalizhuangDajiaoting land subsidence area.

Figure 7 .
Figure 7. Correspondence analysis chart of land subsidence rate versus the land use types.(a-g) The figures of the corrections between nine land use types and seven land subsidence rate classifications.The solid black arrows indicate vectors and the arrows of the direction of the positive vector.The black dotted lines are the extension of the vectors.The red solid arrows indicate perpendicular lines from each land use type to the vectors.The blue solid arrows are the closest to the vectors.

Figure 7 .
Figure 7. Correspondence analysis chart of land subsidence rate versus the land use types.(a-g) The figures of the corrections between nine land use types and seven land subsidence rate classifications.The solid black arrows indicate vectors and the arrows of the direction of the positive vector.The black dotted lines are the extension of the vectors.The red solid arrows indicate perpendicular lines from each land use type to the vectors.The blue solid arrows are the closest to the vectors.

Figure 8 .
Figure 8. Land use types are indicated with different contour colors: the blue is water area and wetland, the dark green is vegetable land, the light green is paddy field and upland soils area, and the purple is peasant-inhabited land.The five gray patterns divide the compressible deposits thickness into four classes with a thickness ranging from 50 to 100 m, 100 to 150 m, 150 to 200 m, and 200 to 260 m.A-A1 and B-B1 are two profiles perpendicular to F1, F2, and F3.

Figure 8 .
Figure 8. Land use types are indicated with different contour colors: the blue is water area and wetland, the dark green is vegetable land, the light green is paddy field and upland soils area, and the purple is peasant-inhabited land.The five gray patterns divide the compressible deposits thickness into four classes with a thickness ranging from 50 to 100 m, 100 to 150 m, 150 to 200 m, and 200 to 260 m.A-A1 and B-B1 are two profiles perpendicular to F1, F2, and F3.

Figure 9 .
Figure 9. Land subsidence rate with profiles A-A1 and B-B1.The red line is the location of the faults.
(b2), R = 0.988) and in the peasant-inhabited (Figure 10(c2), R = 0.748) land than in the paddy fields and in the upland soils (Figure 10(a2), R = 0.491), water area and wetland areas (Figure 10(d2), R = 0.462).In Beijing plain, the large proportion of ground water consumption belongs to the agricultural uses.These four land use types have generally large water consumption.The different land use types could indicate different levels of

Figure 9 .
Figure 9. Land subsidence rate with profiles A-A1 and B-B1.The red line is the location of the faults.
(b2), R = 0.988) and in the peasant-inhabited (Figure 10(c2), R = 0.748) land than in the paddy fields and in the upland soils (Figure 10(a2), R = 0.491), water area and wetland areas (Figure 10(d2), R = 0.462).In Beijing plain, the large proportion of ground water consumption belongs to the agricultural uses.These four land use types have generally large water consumption.The different land use types could indicate different levels of ground water consumption.Overall, changes in the groundwater level are still the main cause of land subsidence.Remote Sens. 2017, 9, 380 13 of 16ground water consumption.Overall, changes in the groundwater level are still the main cause of land subsidence.

Figure 10 .
Figure 10.(a1-d1) show the relationship between ground water level and the land subsidence from 2011 to 2013 in paddy field and upland soils, vegetable land, peasant-inhabited land areas, and water area and wetland.Negative values indicate the decrease of groundwater level and land subsidence.(a2-d2) is the time series relationship between the groundwater level and land subsidence.

Figure 10 .
Figure 10.(a1-d1) show the relationship between ground water level and the land subsidence from 2011 to 2013 in paddy field and upland soils, vegetable land, peasant-inhabited land areas, and water area and wetland.Negative values indicate the decrease of groundwater level and land subsidence.(a2-d2) is the time series relationship between the groundwater level and land subsidence.

Figure 11 .
Figure 11.Distribution of land subsidence rate in the different compressible thickness layers.

Table 1 .
Comparisons of annual average land subsidence rate (mm/year) between SBAS-derived results and ground leveling data at the 6 benchmarks.SDFP: single-look slowly-decorrelating filtered phase.

Table 1 .
Comparisons of annual average land subsidence rate (mm/year) between SBAS-derived results and ground leveling data at the 6 benchmarks.SDFP: single-look slowly-decorrelating filtered phase.

Table 2 .
Land use classification and statistics for the different land use types.

Table 3 .
Land subsidence rate classification.Minus signs represent subsidence.

Table 2 .
Land use classification and statistics for the different land use types.

Table 3 .
Land subsidence rate classification.Minus signs represent subsidence.

Table 4 .
The fractions of SDFP pixels of each deformation in each land use type.

Table 6 .
Characteristics of land subsidence in different land use types.

Table 6 .
Characteristics of land subsidence in different land use types.