Mechanism of Land Subsidence Mutation in Beijing Plain under the Background of Urban Expansion

: Under the background of over-exploitation of groundwater and urban expansion, the land subsidence in the Beijing Plain has dramatically increased recently, and has demonstrated obvious mutation characteristics. Firstly, this paper used the land-use transfer matrix (LUTM) to quantify the urban expansion of Beijing, from 1990 to 2015. Secondly, the gravity center migration model (GCM) and standard deviation ellipse (SDE) methods were employed in order to quantitatively reveal the response relationship between urban expansion and land subsidence in the study area. Finally, the research innovatively combines multi-disciplinary (remote sensing, geophysical prospecting, spatial analysis, and hydrogeology), to analyze the mechanism of land subsidence mutation in the Beijing Plain, at multiple scales. The results showed the following: 1. The development direction of the urban expansion and the subsidence bowl (subsidence rate > 50 mm/year) were highly consistent, with values of 116.8 ◦ and 113.3 ◦ , respectively. 2. At the regional scale, the overall spatial distribution of subsidence mutations is controlled by the geological conditions, and the subsidence mutation time was mainly in 2005 and 2015. The area where mutation occurred in 2005 was basically located in the subsidence bowls, and the correlation between the conﬁned water level and the subsidence rate was relatively high (r > 0.62). The area where the settlement mutation occurred in 2015, was mainly located outside the subsidence bowls, and the correlation between the conﬁned water level and the subsidence rate was relatively low (r < 0.71). 3. In the typical subsidence area, the subsidence mutation occurred mostly in the places where the stratigraphic density is reduced, due to human activities (such as groundwater exploitation). Human activities caused the reduction in stratigraphic density, at 20 m and 90 m vertical depth in urban and rural areas, respectively. 4. At the local scale, clusters of subsidence mutation were located in the fault buffer zone, with a lateral inﬂuence range of nearly 1 km in Tongzhou. The scattered settlement mutation is distributed as a spot pattern, and the affected area is relatively small, which basically includes high-rise buildings. scale is relatively large. The subsidence mutation of multiple circular or elliptic marks is the scattered type, without a continuous feature. The areas with relatively large subsidence mutation scales were located on both sides of the fault buffer zone, with the lateral inﬂuence range being approximately 1 km, and the subsidence mutation grid being distributed in strips on both sides of the fault zone. The scattered subsidence mutation is distributed in a spot pattern, and the scale is relatively small and basically includes high-rise buildings. It shows that high-rise buildings are one of the factors affecting the local uneven subsidence, but the lateral inﬂuence range is relatively small.


Introduction
Land subsidence refers to a geological phenomenon that the ground elevation decreases in a certain area, due to natural factors or human engineering activities [1]. At present, land subsidence has occurred in more than 200 countries and regions in the world, including central Mexico, Northern Italy, central Iran, the Bandung Basin in Indonesia, Las Vegas in the United States, and the Beijing-Tianjin-Hebei region of China [2][3][4][5][6][7][8][9]. It has Based on the InSAR, spatial analysis, and other technologies, previous researches focused on the spatial-temporal distribution characteristics of land subsidence. By combining the traditional measurement, the hydrogeology, and other data, they analyzed the land subsidence response relationship with groundwater level, dynamic/static loads, etc. [29,30]. However, in recent years, with the implementation of a series of water management measures, such as SNWDP, emergency pumping, artificial recharge, and self-provided well replacement, land subsidence in the Beijing Plain has changed rapidly [15]. Especially in Tongzhou, there are many new ground fissures, and the characteristics of subsidence mutation are obvious. Therefore, it is essential to study the mutation characteristics, and the formation mechanism, of the land subsidence in the Beijing Plain. By superimposing the acceleration of urbanization, and the efficient use of underground space and other factors, the differences and non-linearities in the evolution of Beijing′s land subsidence have become dramatically complicated [7]. In order to reveal the mechanism that are behind the complex factors, multi-disciplinary research from multiple perspectives is required.  applied the Mann-Kendall mutation test to acquire the information of the subsidence mutation [15]. However, the mechanism of the subsidence mutation was not analyzed with the measured data. It is essential to quantitatively analyze the causes of mutation in land subsidence, for scientific prevention and control of regional subsidence. Therefore, based on the mutation information of the subsidence in the Beijing Plain, from 2004 to 2015, the GCM and SDE were employed, to quantitatively reveal the response relationship between urban expansion and land subsidence in the study area, then innovatively combine multi-source data (new geophysical prospecting, aerial photogrammetry, and long-term dynamic measured data of confined water level); the mutation mecha- Based on the InSAR, spatial analysis, and other technologies, previous researches focused on the spatial-temporal distribution characteristics of land subsidence. By combining the traditional measurement, the hydrogeology, and other data, they analyzed the land subsidence response relationship with groundwater level, dynamic/static loads, etc. [29,30]. However, in recent years, with the implementation of a series of water management measures, such as SNWDP, emergency pumping, artificial recharge, and self-provided well replacement, land subsidence in the Beijing Plain has changed rapidly [15]. Especially in Tongzhou, there are many new ground fissures, and the characteristics of subsidence mutation are obvious. Therefore, it is essential to study the mutation characteristics, and the formation mechanism, of the land subsidence in the Beijing Plain. By superimposing the acceleration of urbanization, and the efficient use of underground space and other factors, the differences and non-linearities in the evolution of Beijing s land subsidence have become dramatically complicated [7]. In order to reveal the mechanism that are behind the complex factors, multi-disciplinary research from multiple perspectives is required.  applied the Mann-Kendall mutation test to acquire the information of the subsidence mutation [15]. However, the mechanism of the subsidence mutation was not analyzed with the measured data. It is essential to quantitatively analyze the causes of mutation in land subsidence, for scientific prevention and control of regional subsidence. Therefore, based on the mutation information of the subsidence in the Beijing Plain, from 2004 to 2015, the GCM and SDE were employed, to quantitatively reveal the response relationship between urban expansion and land subsidence in the study area, then innovatively combine multi-source data (new geophysical prospecting, aerial photogrammetry, and long-term dynamic measured data of confined water level); the mutation mechanism of subsidence in the Beijing Plain was analyzed on multi-scales. Firstly, on the regional scale, Remote Sens. 2021, 13, 3086 4 of 21 the geology and hydrogeology conditions were used to analyze the main controlling factors of subsidence mutation in the study site. Combined with the long-term measurements of confined water level, the relationship between the subsidence mutation and confined water level was acquired. Secondly, in the typical areas of land subsidence, seismic frequency resonance (SFR) was used to obtain the formation density information. Combined with the data of the groundwater pumping layer, etc., the information of formation density reduction was acquired at the location of subsidence mutation. Thirdly, on the local scale, combining the building height information obtained by photogrammetry technology, the response relationship between high-rise buildings and subsidence mutation was explored, on the background of a relatively single geological structure.

Study Area
Beijing is in the north of the North China Plain. In general, the part of the northwest is high, and the southeast is low [22]. The Beijing Plain is in the southeast of Beijing, covering an area of nearly 6236 km 2 , and accounting for 38% of the total area of Beijing. It has a continental monsoon climate, with distinctive seasons and an average annual temperature of 11.7 • C [10]. From 1978 to 2016, the average annual precipitation in Beijing was 540.2 mm/year, which was mainly concentrated in June to September, accounting for over 80% of the total annual precipitation [7].
The Beijing Plain is composed of alluvial and diluvial fan groups, which are formed by the combined action of five river systems (Juma River, Chaobai River, Yongding River, and Beiyun and Jiyun River) [10]. The quaternary loose sediments are distributed throughout the plain widely, and their thickness varies greatly in space. In general, from the piedmont to the alluvial plain, the thickness of the quaternary gradually increases, and the sediment particles gradually become finer [31]. The lithology is gradually transformed from a single sand and gravel layer to alternating sand, sand gravel, and clay soil. According to the groundwater utilization and hydrogeological conditions, there are four vertical-direction aquifer layers of groundwater in the Beijing Plain, including the phreatic aquifer (less than 50 m), the shallow confined water (50-100 m), the middle-deep confined aquifer (100-180 m), and the deep confined water (180-300 m) [7].
There are several NE (NNE) and NW faults in the Beijing Plain, and the main quaternary active faults include the Huangzhuang-Gaoliying fault, Nanyuan-Tongxian fault, Xiadian-Mafang fault, and Nankou-Sunhe fault ( Figure 1). The spatial distribution of land subsidence evolution in the Beijing Plain is controlled by active faults, which is closely related to the quaternary sediment thickness, lithology, and other factors. The serious subsidence is mainly located in the quaternary depression and the interlaced area of the alluvial fan. Figure 1b shows the Beijing subsidence bowl's spatial distribution.

Material
In order to investigate the mechanism of subsidence mutation in Beijing Plain, the data used in this article mainly include remote sensing data, hydrogeological data and new geophysical data. Among them, remote sensing data includes radar data, photogrammetry, optical remote sensing data and shuttle radar topography mission (SRTM) data. The hydrogeological data mainly includes the confined water level measurements at long-term observation hole and the information of groundwater bowl. The SFR data are 0-200 m underground density information of the Beijing typical land subsidence area. The flowchart of the processing is shown in Figure 2.

SAR Data
The used data in this article are 47 ascending ENVISAT ASAR datasets (from 2003 to 2010) and 48 descending RADARSAT-2 datasets (from 2010 to 2015) to derived regional deformation of the Beijing Plain from 2003 to 2015. The SBAS-InSAR and Quasi-PS InSAR were applied to derive time-series deformation based on the ENVISAT ASAR and RA-DARSAT-2 dataset, respectively. The specific processing parameters are shown in Guo et al. [15].

Land-Use Data
The land-use data comes from the Resource and Environmental Science Data Center of the Chinese Academy of Sciences (http://www.resdc.com), and the LUCC (land-use and land-cover change) classification system is adopted [32]. Figure 3 shows the land-use types of Beijing Plain from 1990 to 2015. Among them, Landsat TM/ETM remote sensing image data were used for remote sensing interpretation in 1990,1995,2000,2005 and 2010, and landsat-8 data were used for remote sensing interpretation in 2015. The land-use types

SAR Data
The used data in this article are 47 ascending ENVISAT ASAR datasets (from 2003 to 2010) and 48 descending RADARSAT-2 datasets (from 2010 to 2015) to derived regional deformation of the Beijing Plain from 2003 to 2015. The SBAS-InSAR and Quasi-PS In-SAR were applied to derive time-series deformation based on the ENVISAT ASAR and RADARSAT-2 dataset, respectively. The specific processing parameters are shown in Guo et al. [15].

Land-Use Data
The land-use data comes from the Resource and Environmental Science Data Center of the Chinese Academy of Sciences (http://www.resdc.com, accessed on 28 June 2021), and the LUCC (land-use and land-cover change) classification system is adopted [32].

Groundwater Level Data
In this study, the water level measurements at 6 confined water long-term observation wells and the distribution of groundwater bowl (2012) in the plain are collected from Beijing Institute of Hydrogeology and Engineering Geology. The distribution of monitoring wells is displayed in Figure 1, and the specific observation well information is in Table  1. In this paper, SFR is employed to detect the formation density information of 0-200 m in the typical land subsidence area of Beijing Plain. The specific profile position is shown in Figure 1. SFR uses the resonance relationship between the seismic wave and geological body to image the propagation function of the geological body, and then inverts the density of formation, which has a high vertical identification ability [7]. In 2020,

Groundwater Level Data
In this study, the water level measurements at 6 confined water long-term observation wells and the distribution of groundwater bowl (2012) in the plain are collected from Beijing Institute of Hydrogeology and Engineering Geology. The distribution of monitoring wells is displayed in Figure 1, and the specific observation well information is in Table 1. In this paper, SFR is employed to detect the formation density information of 0-200 m in the typical land subsidence area of Beijing Plain. The specific profile position is shown in Figure 1. SFR uses the resonance relationship between the seismic wave and geological body to image the propagation function of the geological body, and then inverts the density of formation, which has a high vertical identification ability [7]. In 2020, Guo et al. developed a geological model based on the boreholes and SFR data in TongZhou in Beijing, and founded that faults control the uneven pattern of the land subsidence.
This article attempts to employ SFR to detect the formation density information in typical subsidence areas and understand the subsidence mutation mechanism. The reasons are as follows: the causes of occurrence and evolution of land subsidence in Beijing Plain contain natural and human factors [30]. Natural factors include geological structure, soil secondary consolidation, etc. Human factors mainly include the excessive pumping of groundwater, the efficient utilization of underground space, and the effective stress of dynamic/static loads. The human factors can also be expressed by the alteration of density information of underground space, which can be interpreted by new geophysical technology. Moreover, SFR has a high depth of vertical identification, which basically covers the range of human intervention in underground space (groundwater extraction, building foundation, underground rail transit, etc.).
In this research, the subsidence gradient, spatial distribution of urban areas, rural areas and other information are comprehensively considered to determine the location of the transect where SFR technology is employed, to detect the 0-200 m density information of underground media along the profile of Tongzhou. The measuring points in this experiment are arranged by linear distribution. The test points are spaced 20 m apart, starting point is NO:1000, ending point is NO:16040, and the total length is m. A wide-frequency seismic geophone of 0.2-220 Hz was used, with a single acquisition time of 1 h and an actual recording length of 1800 s. The indoor data processing and imaging stages include multiple superposition of frequency domain data, noise elimination, burnishing treatment, ratio analysis and wave impedance change analysis. Finally, the seismic profile results were acquired, and the SFR results were validated with three hydrogeological boreholes.

Building Height Data
To further analyze the formation of mutation land subsidence in typical areas, this article collects the individual building height information from photogrammetric interpretation covering Tongzhou area. The data are from Beijing Information Resource Management Center, which is a point vector file, containing 35,215 pieces of building information. The property is the height information of the single building. The height accuracy error is within 0.3 m.

Methods
By employing the land-use transfer matrix (LUTM), the gravity center migration (GCM) and the standard deviation ellipse (SDE), the expansion process of Beijing's urbanization was quantitatively described, so as to provide a basis for the subsequent discussion on the response relationship between urban expansion and land subsidence.

Land-Use Transfer Matrix
LUTM is the application of Markov model in land use [33]. It can clearly reveal the changes in land-use types at the beginning and the end of the research, and reflect the direction and quantity of transfer. The specific transfer matrix formula is as follows: Among them, i and j represent n land-use types at the beginning and end of the research (i, j = 1, 2, . . . , n). Further, c ij represents the area of land class i (the beginning of the research) converted to land class j (the end of the research). Combined the spatial analysis, the LUTM of Beijing Eastern Plain was acquired from 1990 to 2015.

Gravity Center Migration Model
The GCM is a method to describe the spatial clustering degree of geographical phenomena, and can quantitatively express the spatial-temporal evolution characteristics of Remote Sens. 2021, 13, 3086 8 of 21 geographic phenomena [34]. Construction land is an essential indicator of human activities, and the more the frequent economic activities, the larger the construction land area is. In order to accurately grasp the trend of urban expansion and explore its response to land subsidence, this paper uses the GCM to calculate the barycenter track of construction land in the Beijing Eastern Plain (Beijing Chaoyang-Tongzhou, BCT) from 1990 to 2015 and the barycenter track of land subsidence bowl (the settlement rate greater than 50 mm/year) from 2005 to 2015.
The spatial barycenter coordinate of construction/land subsidence bowl is shown as follows: where X m and Y m are the space barycenter coordinates of the construction land/subsidence bowl in the m-th year, and C mi represents the area of the construction land/subsidence bowl in zone i. Further, X i and Y i represent the geometric barycenter coordinates of the zone i, and n represents the total number of zones in the study area.

Standard Deviation Ellipse Analysis
The standard deviation ellipse (SDE) was first proposed by Lefever (1926) to reveal the spatial distribution characteristics of geographical elements [35]. It can summarize the spatial distribution of elements, and identify the trend direction of geographical phenomena [36]. Among them, the long axis represents the gathering direction of the center of gravity, the short axis represents the discrete direction of the center of gravity, and the greater the oblateness of ellipse is, the stronger the directivity is.
The SDE of the gravity center of construction land (1990-2015) and the gravity center of subsidence bowl (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) were both calculated in Chaoyang-Tongzhou area of Beijing. The formula of the center coordinate of SDE is as follows: Among them, x i and y i are the coordinates of point i, and n represents the total number of points. The calculation formula of azimuth is as follows: The standard deviation of the x and y axes is as follows: The θ is the azimuth of the ellipse, expressed as the angle formed by clockwise rotation from north to the major axis of the ellipse.

Surface Displacements and Subsidence Mutation in Beijing Plain from 2004 to 2015
To validate the results, the vertical displacements were derived from the line-of-sight (LOS) displacements, which were measured from the PS-InSAR technique during the two time periods [37]. The benchmarks are compared with the average displacement of the PS points, within the radius of 150 m. From 2004 to 2010, the correlation coefficient between the level monitoring data and InSAR monitoring result was 0.99, and the standard deviation was 3 mm/year. From 2011 to 2013, the correlation was 0.94, and the standard deviation was 3.6 mm/year.
According to previous studies, this paper ignored the horizontal deformation of the Beijing Plain [37,38]. Also, with their method, the displacement in the LOS direction was converted into the vertical displacement [37]. As displayed in Figure 4, the subsidence in the Beijing Plain, from 2004 to 2015, was relatively serious, and the spatial difference of the subsidence was relatively obvious. The is the azimuth of the ellipse, expressed as the angle formed by clockwise rotation from north to the major axis of the ellipse.

Surface Displacements and Subsidence Mutation in Beijing Plain from 2004 to 2015
To validate the results, the vertical displacements were derived from the line-of-sight (LOS) displacements, which were measured from the PS-InSAR technique during the two time periods [37]. The benchmarks are compared with the average displacement of the PS points, within the radius of 150 m. From 2004 to 2010, the correlation coefficient between the level monitoring data and InSAR monitoring result was 0.99, and the standard deviation was 3 mm/year. From 2011 to 2013, the correlation was 0.94, and the standard deviation was 3.6 mm/year.
According to previous studies, this paper ignored the horizontal deformation of the Beijing Plain [37,38]. Also, with their method, the displacement in the LOS direction was converted into the vertical displacement [37]. As displayed in Figure 4, the subsidence in the Beijing Plain, from 2004 to 2015, was relatively serious, and the spatial difference of the subsidence was relatively obvious.    [15]. Among the single-year mutation grids, the number of grids that mutated in 2015 and 2005 was the largest, with 1344 and 915, respectively, showing certain clustering characteristics in space. Among them, the grids with the subsidence mutation occurred in 2015, and are primarily distributed in the middle and lower part of the alluvial fans; some of them are located on the outer edge of the ST subsidence bowl, and others are located in the non-subsidence area. The grids with a subsidence mutation in 2005 are mainly distributed in the Nankou diluvial fan (NKDF), and the intersection of the Chaobai River alluvial fan (CBRAF) and the Yongding River alluvial fan (YDRAF), and they are mostly located in the subsidence bowl. There were 152 grids with a subsidence mutation

Urban Expansion in BCT Region from 1990 to 2015
The land-use transfer matrix of the BCT area, from 1990 to 2015, is shown in Table 2. As shown in the table, 55.35% of the BCT area had the same land-use pattern from 1990 to 2015, while, for the rest, approximately 44.65% of that changed the pattern. Among them, the area that was occupied by construction land had increased, from 300.96 km 2 in 1990 to 867.73 km 2 in 2015, in particular, with an increase rate of 288.32%. The land-use type that decreased the most was cultivated land, which decreased from 1001.34 km 2 in 1990 to 470.02 km 2 in 2015, with a reduction rate of −53.06%. The type of land use that changed the most was the cultivated land, which converts into construction land, with an area of 541.74 km 2 . It can be observed that the BCT area expanded significantly from 1990 to 2015. In order to further describe the urban expansion intuitively, the percentage of BCT farmland and construction land in the total area in six periods (from 1990 to 2015), is calculated in this research, as shown in Figure 6. Among them, the proportion of cultivated land decreased from 74% in 1990 to 35% in 2015, and the proportion of construction land increased from 22% in 1990 to 64% in 2015.

Urban Expansion in BCT Region from 1990 to 2015
The land-use transfer matrix of the BCT area, from 1990 to 2015, is shown in Table 2. As shown in the table, 55.35% of the BCT area had the same land-use pattern from 1990 to 2015, while, for the rest, approximately 44.65% of that changed the pattern. Among them, the area that was occupied by construction land had increased, from 300.96 km 2 in 1990 to 867.73 km 2 in 2015, in particular, with an increase rate of 288.32%. The land-use type that decreased the most was cultivated land, which decreased from 1001.34 km 2 in 1990 to 470.02 km 2 in 2015, with a reduction rate of −53.06%. The type of land use that changed the most was the cultivated land, which converts into construction land, with an area of 541.74 km 2 . It can be observed that the BCT area expanded significantly from 1990 to 2015. In order to further describe the urban expansion intuitively, the percentage of BCT farmland and construction land in the total area in six periods (from 1990 to 2015), is calculated in this research, as shown in Figure 6. Among them, the proportion of cultivated land decreased from 74% in 1990 to 35% in 2015, and the proportion of construction land increased from 22% in 1990 to 64% in 2015. Remote Sens. 2021, 13, x FOR PEER REVIEW 11 of 21

The Result of Formation Density in Typical Subsidence Area via SFR
In order to focus on the analysis of the mechanism of the subsidence mutation in 2013, in the Tongzhou typical area, the density information of the underground formation (0-200 m) was acquired along the profile (Figure 7). Among them, the color code from red to blue represents the density information of the underground formation, and the white lines represent the hydrogeological boreholes. According to the well seismic profile, the hierarchical formation information that was obtained by SFR technology, is highly accurate. It should be noted that  As show in Figure 7, due to greater human intervention factors, the density of formation is relatively low at the crossroads compared with the density of the adjacent points. In rural areas, the decrease in formation density, caused by human factors (pumping groundwater, etc.), is relatively small, with a vertical influenced depth of approximately 20 m. In urban areas, the decrease in formation density, caused by human factors (pumping groundwater, underground space utilization, etc.), is dramatically serious, with the vertical influence depth reaching 90 m.

The Result of Formation Density in Typical Subsidence Area via SFR
In order to focus on the analysis of the mechanism of the subsidence mutation in 2013, in the Tongzhou typical area, the density information of the underground formation (0-200 m) was acquired along the profile (Figure 7). Among them, the color code from red to blue represents the density information of the underground formation, and the white lines represent the hydrogeological boreholes. According to the well seismic profile, the hierarchical formation information that was obtained by SFR technology, is highly accurate.

The Result of Formation Density in Typical Subsidence Area via SFR
In order to focus on the analysis of the mechanism of the subsidence mutation in 2013, in the Tongzhou typical area, the density information of the underground formation (0-200 m) was acquired along the profile (Figure 7). Among them, the color code from red to blue represents the density information of the underground formation, and the white lines represent the hydrogeological boreholes. According to the well seismic profile, the hierarchical formation information that was obtained by SFR technology, is highly accurate. It should be noted that  As show in Figure 7, due to greater human intervention factors, the density of formation is relatively low at the crossroads compared with the density of the adjacent points. In rural areas, the decrease in formation density, caused by human factors (pumping groundwater, etc.), is relatively small, with a vertical influenced depth of approximately 20 m. In urban areas, the decrease in formation density, caused by human factors (pumping groundwater, underground space utilization, etc.), is dramatically serious, with the vertical influence depth reaching 90 m. As show in Figure 7, due to greater human intervention factors, the density of formation is relatively low at the crossroads compared with the density of the adjacent points. In rural areas, the decrease in formation density, caused by human factors (pumping groundwater, etc.), is relatively small, with a vertical influenced depth of approximately 20 m. In urban areas, the decrease in formation density, caused by human factors (pumping groundwater, underground space utilization, etc.), is dramatically serious, with the vertical influence depth reaching 90 m.

Relationship between Urban Expansion and Land Subsidence
As shown in Figure 8, the points of the blue triangle were the gravity center migration of the construction land in the BCT area, from 1990 to 2015, and the points of the darkred circle were the gravity center migrations of the ground subsidence bowl (settlement rate > 50 mm/year), from 2004 to 2015. The blue line is the standard deviation ellipse of the construction land in the BCT area, from 1990 to 2015, and the red line is the standard deviation ellipse of the land subsidence bowl, from 2004 to 2015. Both have obvious directionality, and the standard deviation ellipticity of construction land is larger and more directional.

Relationship between Urban Expansion and Land Subsidence
As shown in Figure 8, the points of the blue triangle were the gravity center migration of the construction land in the BCT area, from 1990 to 2015, and the points of the dark-red circle were the gravity center migrations of the ground subsidence bowl (settlement rate > 50 mm/year), from 2004 to 2015. The blue line is the standard deviation ellipse of the construction land in the BCT area, from 1990 to 2015, and the red line is the standard deviation ellipse of the land subsidence bowl, from 2004 to 2015. Both have obvious directionality, and the standard deviation ellipticity of construction land is larger and more directional. The shift direction of the construction land gravity center is 116.8° from the north, and the subsidence bowl is 113.3°. It shows that the development direction of the land subsidence bowl is highly consistent with the urban expansion direction. Moreover, the urban gravity center of the BCT area was gradually expanding to the most severe subsidence region in the DB subsidence bowl Figure 8c. As can be observed from Table 3, the areas of the construction land and subsidence bowl both increased from 2005 to 2015, and the correlation coefficient was 0.81, indicating the high response relationship between land subsidence and urban expansion. Therefore, on the background of urban expansion, the subsidence mutation was analyzed in the Beijing Plain. The shift direction of the construction land gravity center is 116.8 • from the north, and the subsidence bowl is 113.3 • . It shows that the development direction of the land subsidence bowl is highly consistent with the urban expansion direction. Moreover, the urban gravity center of the BCT area was gradually expanding to the most severe subsidence region in the DB subsidence bowl Figure 8c. As can be observed from Table 3, the areas of the construction land and subsidence bowl both increased from 2005 to 2015, and the correlation coefficient was 0.81, indicating the high response relationship between land subsidence and urban expansion. Therefore, on the background of urban expansion, the subsidence mutation was analyzed in the Beijing Plain.

Analysis of Subsidence Mutation Mechanism at Regional Scale
In this section, the subsidence mutation mechanism of the Beijing Plain is analyzed at a regional scale. The multi-year subsidence mutations are mainly located in the middle-upper part of the alluvial fans, while the single-year mutations are mainly located in the middle-lower part of the alluvial fans Figure 5a. On both sides of the confined water bowl, the difference in subsidence mutations is relatively large Figure 5b. The spatial pattern of the land subsidence mutation in the Beijing Plain is mainly controlled by active faults. The grids with a single-year subsidence mutation are concentrated in 2005 and 2015, and are relatively concentrated in space.
At the same time, the spatial distribution of the bowls is relatively consistent between land subsidence and the second confined aquifer. Chen et al. (2016) indicated that the second confined aquifer (100-180 m) is the main contributor of ground subsidence in the Beijing Plain [27]. Hence, to further explore the response relationship between the subsidence mutation and groundwater level, this research collected six long-term observation wells (confined aquifers), of which the subsidence mutation years are 2005 and 2015. The specific groundwater observation wells information are shown in Table 4. The confined water observation well information, and its spatial distribution, are shown in Figure 1. This paper used basic geology, hydrogeology, and other data, to analyze the main control factors of the subsidence mutation in 2005 and 2015, at the regional scale. The confined aquifer bowl (2012) and the measured data of the confined groundwater level were applied, to analyze the response relationship between the subsidence mutation and confined groundwater level. Most of the grids with a subsidence mutation year in 2005 were located in the serious subsidence areas, which are mainly distributed in the Chaoyang-Tongzhou (CT, DB, CT, and CL) subsidence bowls and HS subsidence bowl. Among them, in the Chaoyang-Tongzhou area, the subsidence area with a mutation year in 2005 is located at the intersection of YDRAF and CBRAF. The quaternary stratum is relatively thick, the interbedding is dramatically increased, the particles became relatively smaller, and the compressibility became larger, which provides conditions for the development of land subsidence. In addition, affected by the exploitation of Huairou EGRR (operated in August 2003) and Pinggu EGRR (operated in August 2004), the land subsidence has a certain lag on the groundwater level, thus the land subsidence mutation in Chaoyang-Tongzhou happened in 2005.
The NKDF with a subsidence mutation year that occurred in 2005, belongs to the quaternary sag area, and the thickness of the quaternary system increased significantly, which provides the conditions for the occurrence of subsidence. Moreover, this area belongs to the same groundwater system as the middle-upper region of CBRAF, so the decrease in the confined groundwater level, caused by emergency pumping, was the main incentive for the subsidence mutation of this area in 2005.
According to the data from long-term confined groundwater observation stations Figure 9, in the three confined groundwater observation stations with a subsidence mutation in 2005, the confined groundwater levels all decreased, and the confined groundwater level of the S1, S2 and S3 stations decreased by 1.01 m, 5.51 m, and 1.41 m, respectively. The confined groundwater level, and the overall subsidence rate, were highly correlated, which were 0.62, 0.98, and 0.86, respectively.
Tongzhou area, the subsidence area with a mutation year in 2005 is located at the intersection of YDRAF and CBRAF. The quaternary stratum is relatively thick, the interbedding is dramatically increased, the particles became relatively smaller, and the compressibility became larger, which provides conditions for the development of land subsidence. In addition, affected by the exploitation of Huairou EGRR (operated in August 2003) and Pinggu EGRR (operated in August 2004), the land subsidence has a certain lag on the groundwater level, thus the land subsidence mutation in Chaoyang-Tongzhou happened in 2005.
The NKDF with a subsidence mutation year that occurred in 2005, belongs to the quaternary sag area, and the thickness of the quaternary system increased significantly, which provides the conditions for the occurrence of subsidence. Moreover, this area belongs to the same groundwater system as the middle-upper region of CBRAF, so the decrease in the confined groundwater level, caused by emergency pumping, was the main incentive for the subsidence mutation of this area in 2005.
According to the data from long-term confined groundwater observation stations Figure 9, in the three confined groundwater observation stations with a subsidence mutation in 2005, the confined groundwater levels all decreased, and the confined groundwater level of the S1, S2 and S3 stations decreased by 1.01 m, 5.51 m, and 1.41 m, respectively. The confined groundwater level, and the overall subsidence rate, were highly correlated, which were 0.62, 0.98, and 0.86, respectively.  Among them, station S1 is located in the middle-upper part of YDRAF, with good recharge conditions. The groundwater level generally showed good fluctuation from 2004 to 2015. However, from January 2004 to December 2005, the water level of the confined groundwater dropped sharply, from −1.16 m to −6.88 m, and the groundwater level decreased approximately linearly, far higher than the average groundwater drop level in 12 years (0.3 m). Thus, the above reasons are the main factors that caused the subsidence mutation in the grid of the S1 station in 2005.
Station S2 is located in the interlaced area of CBRAF and YDRAF, where the interbedded layers increased and formation particles became smaller. Among the six groundwater monitoring stations, the correlation coefficient between the confined groundwater level and subsidence rate is the highest, reaching 0.98. From 2005 to 2015, the fluctuation of the water level was relatively small, and showed a continuous downward trend. Among them, the confined groundwater level had decreased 5.51 m in 2005, with the largest decline in 12 years, which was the main inducement of the subsidence mutation in the S2 station in 2005.
Station S3 is located in the NKDF, which belongs to the quaternary sag area, and the quaternary thickness increases significantly. From 2005 to 2015, the fluctuation of the confined groundwater level was relatively strong, and generally showed a decreasing trend. The correlation coefficient between the confined groundwater level and the subsidence was relatively high, reaching 0.86. In 2005, the confined groundwater level of the S3 station decreased by 1.41 m. Compared with the S1 and S2 observation wells, the subsidence fluctuation in this area was slightly stronger.
To sum up, the geological condition was the main factor for the grid with the subsidence mutation in 2005, and the response relationship of the land subsidence rate and the confined groundwater level was relatively high. Due to emergency pumping, the rapid decline in the confined groundwater level was the main inducement of the subsidence mutation in the areas of the Beijing Plain in 2005.

The Land Subsidence Mutation in 2015
The subsidence area with the mutation in 2015 was mainly located in the middleupper part of the CBRAF, and the middle-lower part of the CBRAF and YDRAF. Among them, in the middle-lower part of the alluvial fan, with a subsidence mutation in 2015, the subsidence rate in this area was relatively lower from 2004 to 2015, and was located outside the subsidence bowl.
In the middle-upper part of the CBRAF, subsidence is relatively serious at the outer edge of the ST bowl. There are many water sources in this area, such as the 8th Waterworks, EGRR, etc. [1]. The target pumping layer is the quaternary loose rock hole water. Before the operation of SNDWP (on 12 December 2014), the long-term groundwater pumping in this region ensures the safety of the water supply in the Capital area, but induces relatively serious land subsidence [31]. After the operation of SNDWP, the daily pumping volume in the Chaobai River EGRR decreased from 30 × 10 4 /m 3 to 10 × 10 4 /m 3 , entering the stage of conservation. In addition, from August to November 2015, artificial recharge was carried out in the Niulanshan area, through natural river channels (the amount of recharge was 33.79 × 10 6 /m 3 ), and the groundwater level in some areas increased significantly. Therefore, groundwater reduction and artificial recharge were the main factors that caused subsidence mutation in the region in 2015.
According to long-term observations of confined groundwater Figure 10, in the three confined groundwater observation stations with a subsidence mutation in 2015, the confined groundwater level of the S4, S5, and S6 stations increased by 1.2 m, 1.97 m, and 0.59 m, respectively. Compared to the grids with a subsidence mutation in 2005, the correlation coefficient between the confined groundwater level and the total subsidence was slightly lower, with 0.54, 0.71, and 0.15, respectively.
The S4 and S5 stations are in the middle and upper part of the CBRAF, and both are located at the outer edge of the ST subsidence bowl. From 2005 to 2015, the confined groundwater level over both the stations fluctuated seasonally, and the overall trend was downward. The subsidence of the two observation stations also fluctuated seasonally, and generally showed a downward trend. According to Wang et al. (2014), if the land subsidence occurs at a location with good recharge conditions, after the groundwater level drops rapidly, the compression will slow down, once the groundwater level bounces back. Among them, the confined groundwater level of the S4 station increased by 1.2 m in 2015. In 2015, the water level of the S5 station increased by 1.97 m, which was the year with the highest water level height in 12 years, and the S5 station showed high correlation between the confined groundwater level and the subsidence rate, to some extent. Thus, the increase in the confined groundwater level is one of the main factors for the subsidence mutation in this area in 2015. Remote Sens. 2021, 13, x FOR PEER REVIEW 16 of 21 The S4 and S5 stations are in the middle and upper part of the CBRAF, and both are located at the outer edge of the ST subsidence bowl. From 2005 to 2015, the confined groundwater level over both the stations fluctuated seasonally, and the overall trend was downward. The subsidence of the two observation stations also fluctuated seasonally, and generally showed a downward trend. According to Wang et al. (2014), if the land subsidence occurs at a location with good recharge conditions, after the groundwater level drops rapidly, the compression will slow down, once the groundwater level bounces back. Among them, the confined groundwater level of the S4 station increased by 1.2 m in 2015. In 2015, the water level of the S5 station increased by 1.97 m, which was the year with the highest water level height in 12 years, and the S5 station showed high correlation between the confined groundwater level and the subsidence rate, to some extent. Thus, the increase in the confined groundwater level is one of the main factors for the subsidence mutation in this area in 2015.
The S6 station is in the middle-upper part of YDRAF. The accumulated subsidence in 12 years is only 15.57 mm. From 2004 to 2010, the confined groundwater level of the S6 station fluctuated slightly, and the groundwater level increased smoothly. The confined groundwater level continued to decrease from 2010 to 2012, and increased after 2012. The subsidence in this area is mainly elastic deformation, and the correlation between the confined water level and subsidence rate is the lowest.
In summary, most of the grids with a subsidence mutation in 2015 were located outside the edge of the ST subsidence bowl and the non-subsidence area. Among them, the relationship between the confined water level and the subsidence rate is relatively high in the middle-upper regions of the CBRAF. The increase in the confined groundwater level The S6 station is in the middle-upper part of YDRAF. The accumulated subsidence in 12 years is only 15.57 mm. From 2004 to 2010, the confined groundwater level of the S6 station fluctuated slightly, and the groundwater level increased smoothly. The confined groundwater level continued to decrease from 2010 to 2012, and increased after 2012. The subsidence in this area is mainly elastic deformation, and the correlation between the confined water level and subsidence rate is the lowest.
In summary, most of the grids with a subsidence mutation in 2015 were located outside the edge of the ST subsidence bowl and the non-subsidence area. Among them, the relationship between the confined water level and the subsidence rate is relatively high in the middle-upper regions of the CBRAF. The increase in the confined groundwater level is one of the main reasons for the subsidence mutation in 2015. In the non-subsidence area, the seasonal subsidence is relatively obvious, and the correlation between the confined groundwater level and land subsidence is minimal.

The Response Relationship between Formation Density and Subsidence Mutation in Tongzhou Typical Subsidence Area
The grids with subsidence mutations in 2005 and 2013 were spatially clustered, and most of the areas with a subsidence mutation in 2013 were in Tongzhou. The uneven subsidence is serious in urban, rural, and transitional areas in the Tongzhou District. In September 2012, the district became the administrative sub-center of the Capital, and the urban area grew speedily [7]. This research employed the SFR to detect the response relationship between the formation density and the subsidence mutation in this area (Figure 11a). In Figure 11b, the blue line is the position of the profile, the orange point is the boundary point of the subsidence mutation, and Figure 11c is the formation density information corresponding to the profile. It should be noted that from No:1000 to No:10800 there are rural areas, and the main groundwater pumping layer is approximately within 100 m. The No:10800 to No:16040 points are urban areas, and the main groundwater pumping layer is approximately within 300 m.
The grids with subsidence mutations in 2005 and 2013 were spatially clustered, and most of the areas with a subsidence mutation in 2013 were in Tongzhou. The uneven subsidence is serious in urban, rural, and transitional areas in the Tongzhou District. In September 2012, the district became the administrative sub-center of the Capital, and the urban area grew speedily [7]. This research employed the SFR to detect the response relationship between the formation density and the subsidence mutation in this area ( Figure  11a). In Figure 11b, the blue line is the position of the profile, the orange point is the boundary point of the subsidence mutation, and Figure 11c is the formation density information corresponding to the profile. It should be noted that from No:1000 to No:10800 there are rural areas, and the main groundwater pumping layer is approximately within 100 m. The No:10800 to No:16040 points are urban areas, and the main groundwater pumping layer is approximately within 300 m. Near Suzhuang Village, the subsidence evolution changed abruptly in 2013. From the geophysical profile Figure 11c, it can be observed that from the No:1000 to No:6800 points, there is a set of continuous formation 15-35 m underground, but the formation density decreases suddenly, from the No: 1100 to No:1520 points (block I area). The reason is that the area is the agricultural planting area Figure 11d, and excessive scattered pumps groundwater result in the decrease in the corresponding formation density.
The points from No:4960 to No:8080 are located near Liersi Village, with rural residents. Here, the subsidence mutation happened in 2013. From the geophysical profile Figure 11c, between 50-100 m underground, the density of formation decreased due to excessive groundwater exploitation (block II area), forming multiple groundwater bowls. Through field investigation, there are many orchards in this region, and the demand for groundwater water is relatively great. Overall, the reduction in the formation density, caused by the pumping of groundwater in rural areas, is approximately 20 m.
From the No:10800 to No:16040 points, the subsidence mutation happened in 2005. Among them, from the No:13520 to No:16040 points, there are urban Central Business District (CBD) areas, and the strata density has decreased significantly. Especially near Liyuan (block III area), from the No:14240 to No:15920 points, it formed multiple falling bowls, due to pumping groundwater. On the formation of 100 m underground (the main pumping layer), the vertical influence depth of formation density reduction, caused by pumping groundwater, can reach 90 m.
In summary, the subsidence mutation mainly occurred in the areas where the formation density was significantly reduced due to frequent human activities.

Explore the Response Relationship between High-Rise Buildings and Subsidence Mutation from Local Scale
In the area where the subsidence mutation occurred in Tongzhou in 2013, the overall settlement was low. The grids were mainly distributed in rural areas, with limited underground space utilization, and basically located outside the confined groundwater bowl. Therefore, this article collected the building height information that was interpreted by photogrammetry, and explored the response relationship between high-rise buildings and uneven subsidence on the background, at the local scale (the geological condition is relatively homogenous).
As shown in Figure 12, the 2013 Tongzhou subsidence mutation can be divided into the clustered type and scattered type. Among them, there are three blue rectangles that are marked as the clustered type, and subsidence mutation has a feature of being connected, and its scale is relatively large. The subsidence mutation of multiple circular or elliptic marks is the scattered type, without a continuous feature. The areas with relatively large subsidence mutation scales were located on both sides of the fault buffer zone, with the lateral influence range being approximately 1 km, and the subsidence mutation grid being distributed in strips on both sides of the fault zone. The scattered subsidence mutation is distributed in a spot pattern, and the scale is relatively small and basically includes high-rise buildings. It shows that high-rise buildings are one of the factors affecting the local uneven subsidence, but the lateral influence range is relatively small.

Conclusions
This article used the LUTM to quantify the urban expansion of Beijing, from 1990 to 2015. Then, GCM and SDE methods were employed to quantitatively reveal the response relationship between urban expansion and land subsidence in the Beijing Plain. Finally, based on the geological and hydrogeological conditions, combined multi-sources data, this article analyzes the mechanism of subsidence mutations in the Beijing Plain in 2005, 2013, and 2015, on multi-scales. The main conclusions are as follows: The direction of the urban expansion and the development direction of subsidence bowl in the Beijing Plain, showed obvious directionality and high consistency, with 116.8° and 113.3°, respectively. At the regional scale, the overall spatial distribution of subsidence mutations was controlled by the geological conditions. Most areas with a subsidence mutation in 2005 are located in the subsidence bowls of the Beijing Plain, and the primary controlling factors of the subsidence mutation are the geological background and excessive pumping of groundwater. In 2015, the subsidence mutations were mainly distributed in the middlelower parts of CBRAF and YDRAF. On the background of SNWDP, the increase in the confined groundwater level was one of the factors, which caused the subsidence mutation in 2015.
In the typical subsidence area of Tongzhou, subsidence mutation commonly occurred in the places where the density of formation decreased due to human activities. The vertical influence of formation density reduction, caused by human activities, such as

Conclusions
This article used the LUTM to quantify the urban expansion of Beijing, from 1990 to 2015. Then, GCM and SDE methods were employed to quantitatively reveal the response relationship between urban expansion and land subsidence in the Beijing Plain. Finally, based on the geological and hydrogeological conditions, combined multi-sources data, this article analyzes the mechanism of subsidence mutations in the Beijing Plain in 2005, 2013, and 2015, on multi-scales. The main conclusions are as follows: The direction of the urban expansion and the development direction of subsidence bowl in the Beijing Plain, showed obvious directionality and high consistency, with 116.8 • and 113.3 • , respectively. At the regional scale, the overall spatial distribution of subsidence mutations was controlled by the geological conditions. Most areas with a subsidence mutation in 2005 are located in the subsidence bowls of the Beijing Plain, and the primary controlling factors of the subsidence mutation are the geological background and excessive pumping of groundwater. In 2015, the subsidence mutations were mainly distributed in the middle-lower parts of CBRAF and YDRAF. On the background of SNWDP, the increase in the confined groundwater level was one of the factors, which caused the subsidence mutation in 2015.
In the typical subsidence area of Tongzhou, subsidence mutation commonly occurred in the places where the density of formation decreased due to human activities. The vertical influence of formation density reduction, caused by human activities, such as pumping groundwater, varies in depth between urban and rural areas. In the Tongzhou typical region, the influence depth is approximately 20 m in rural areas, and can reach approximately 90 m in urban areas.
At the local scale, in the case of a relatively homogeneous geological background, in 2013 the distribution of Tongzhou subsidence mutations could be classified into the following two types: clustered and scattered. The clustered type of subsidence mutation basically strips in the buffer zone of the fault zone, and the lateral influence range in the typical Tongzhou area is approximately 1 km. The scattered type of subsidence mutation is distributed as a spot pattern, the distance to the fault zone is relatively large, and the scale is relatively limited, basically including high-rise buildings.