Next Article in Journal
MAT: GIS-Based Morphometry Assessment Tools for Concave Landforms
Previous Article in Journal
Observation of Cloud Base Height and Precipitation Characteristics at a Polar Site Ny-Ålesund, Svalbard Using Ground-Based Remote Sensing and Model Reanalysis
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Three-Dimensional Surface Displacement of the Eastern Beijing Plain, China, Using Ascending and Descending Sentinel-1A/B Images and Leveling Data

1
Key Laboratory of the Ministry of Education Land Subsidence Mechanism and Prevention, Capital Normal University, Beijing 100048, China
2
College of Resources Environment and Tourism, Capital Normal University, Beijing 100048, China
3
Observation and Research Station of Groundwater and Land Subsidence in Beijing-Tianjin-Hebei Plain, Beijing 100048, China
4
Beijing Laboratory of Water Resources Security, Capital Normal University, Beijing 100048, China
5
Beijing Institute of Hydrogeology and Engineering Geology, Beijing 100195, China
*
Author to whom correspondence should be addressed.
Remote Sens. 2021, 13(14), 2809; https://doi.org/10.3390/rs13142809
Submission received: 3 June 2021 / Revised: 9 July 2021 / Accepted: 13 July 2021 / Published: 17 July 2021
(This article belongs to the Section Environmental Remote Sensing)

Abstract

:
Surface displacement is an common environmental geological phenomenon in the Beijing Plain. Research on surface displacement in the Beijing Plain has mainly focused on vertical surface displacement, whereas the horizontal displacement has scarcely been studied. To investigate the 3-D surface displacement in the Beijing Plain, we construct a leveling-constrained multidirectional PS-InSAR 3-D surface displacement estimation method to obtain the 3-D surface displacement information. The results show that the surface displacement in the study area during 2016–2018 was mainly vertical displacement with two main northern and southern subsidence centers; the vertical displacement ranged from −150 mm/year (down) to 5 mm/year (up), and the east–west horizontal displacement ranged from 20 mm/year (east) to 22 mm/year (west). Validation results show that the 3-D surface displacement estimation results agree well with leveling data and GPS data, indicating the reliability of the 3-D surface displacement datasets. The 3-D surface displacement results show that horizontal displacement is obvious in the areas with a large vertical displacement in the eastern Beijing Plain. Additionally, the horizontal displacement is directed toward the center of vertical displacement. The compressive strain is observed close to the centers of vertical displacement, whereas tensile strain occurs far from the centers of vertical displacement. The main cause of the 3-D surface displacement in the study area is the long-term groundwater overexploitation, especially deep groundwater exploitation. The spatial and temporal extents of displacement do not exactly match the locations of the groundwater sinks in different aquifers; instead, geological structures and stratigraphic/lithological conditions may have a combined effect. Moreover, the spatial and temporal distributions of surface displacement are closely related to ground fissure activity, and both influence each other.

1. Introduction

Surface displacement, which refers to the environmental geological phenomenon in which the surface of the Earth’s crust is displaced [1], can be triggered by both natural factors, such as volcanic activity, glacial drift, seismicity, and landslides, and anthropogenic factors, such as mining activities and groundwater extraction. Accordingly, surface displacement constitutes a key issue that must be considered in the planning and construction of cities and large linear projects (e.g., high-speed railroads and large-span bridges; therefore, it is important to effectively and accurately monitor surface displacement. To date, surface displacement caused by the extraction of underground fluids has caused disasters globally, especially in some developed cities and regions, such as Mexico City [2], the Kanto region of Japan [3], Indonesia [4], the Beijing–Tianjin–Hebei region [5,6,7,8,9,10,11,12,13,14,15,16], and the Yangtze River Delta region [17,18,19] in China, where urban safety is under serious threat.
Land subsidence in Beijing was first discovered in 1935 and continues to develop rapidly [20,21]. In 2016, the maximum cumulative vertical displacement in the Beijing Plain reached 1864 mm, and the maximum annual displacement rate exceeded 110 mm/year. The land subsidence in Beijing is distributed mainly throughout the eastern suburbs of Bailizhuang-Dajiuting and the northeastern suburbs of Laiguangying in Chaoyang District and Tongzhou District, Changping Shahe-Baxianzhuang in Changping District, Daxing Yufa-Lixian in Daxing District, and Shunyi Pinggezhuang in Shunyi District. The continuous development of land subsidence has led to a large number of inaccurate leveling points and subsidence, which has aggravated urban flooding and affected the layout and planning of urban construction.
Conventional surface displacement monitoring usually relies on leveling surveys and the Global Positioning System (GPS). While GPS technology can continuously monitor surface displacement, the measurement density is limited by the distribution of expensive ground equipment. Furthermore, leveling measurements are limited by human, material and financial resources: the general deployment of a leveling survey consists of relatively few points, sparse routes, long monitoring cycles, and very low spatial and temporal resolutions. Consequently, it remains difficult to meet the rapid and large-scale monitoring needs of the modern disaster prevention and mitigation of surface displacement. As a new type of spaceborne geodesic measurement tool developed in the past three decades, interferometric synthetic aperture radar (InSAR) has numerous advantages, including the ability to acquire all-day, all-weather measurements over a large area (tens to hundreds of kilometers) with high accuracy (centimeter to millimeter scale) and high spatial resolution (tens to several meters); hence, InSAR technology has been increasingly recognized by experts and scholars and is widely used to monitor earthquakes [22], volcanoes [23], landslides [24], glacial drift [25], tectonic movements [26], groundwater extraction [27], mining [28], etc.
However, most existing studies have employed single-satellite synthetic aperture radar (SAR) image data to extract one-dimensional (1-D) displacement information in the line-of-sight (LOS) direction by differential interferometric synthetic aperture radar (D-InSAR) or persistent scatterer interferometric synthetic aperture radar (PS-InSAR). Nevertheless, only 1-D displacement information can be retrieved in the radar LOS direction by PS-InSAR; that is, three-dimensional (3-D) displacement information is not available. In recent years, through the efforts of some scholars at home and abroad, some methods have been developed to convert 1-D LOS displacement measurements into 3-D displacement fields, and real 3-D displacement information sufficient to represent the surface has been successfully obtained. This information is widely used to investigate surface displacement, volcanic activity, earthquakes, mining [29] and glacial drift. These methods can be broadly classified into one of the following two categories.
The first category is the InSAR 3-D surface displacement estimation technique, which mainly uses weighted least squares to retrieve the 3-D surface displacement field by combining LOS or azimuth measurements from at least two different orbits. Multidirectional LOS vector results can be obtained by the D-InSAR and multidirectional InSAR techniques for solving this 3-D problem [30,31] and can be combined with azimuth measurements obtained by offset tracking or multiple-aperture InSAR (MAI) techniques [32,33,34,35,36,37,38]. Similarly, a priori knowledge can reduce the requirement for SAR datasets [39,40]. The second type of method is to obtain the 3-D surface displacement field mainly by combining InSAR measurements with external data, especially GPS data [41,42,43,44,45,46]. This approach can obtain high-precision displacement results in all three dimensions, but the accuracy of the results depends on the surface displacement characteristics within the study area and the number and spatial distribution of GPS observation points.
Based on the accurate acquisition of 3-D surface displacement, many scholars worldwide have conducted research on the multidisciplinary retrieval of 3-D displacement information in the fields of hydrology, geology and others to investigate the causes and evolutionary mechanisms of 3-D displacement and ground fissures. For example, Qu et al. [47], Liu et al. [48] and Yang et al. [49] studied the relationship between ground subsidence and ground fissures in the Fenwei Basin and confirmed the dominance of vertical displacement in the surface displacement induced by groundwater overexploitation; they further verified that areas with large east–west horizontal displacement corresponded to areas of active ground fissures with large vertical displacement. In a study of aquifer dynamics in the Salt Lake Valley, Hu et al. [50] found that east–west horizontal displacement and vertical displacement occurred simultaneously and were both seasonal in nature, which was also reported in the Santa Clara Valley by Chaussard et al. [51].
Although many studies on the eastern Beijing Plain have shown that the surface displacement therein is mainly vertically oriented and slow with small amounts of horizontal displacement, the long-term cumulative effect of this displacement is not negligible, so it is necessary to accurately obtain the regional 3-D displacement field. Our research group has been engaged in the study of regional ground subsidence for a long time; as a consequence, some levelling data have been accumulated with a high point monitoring accuracy. Moreover, among the various available 3-D surface displacement estimation methods, multidirectional InSAR offers the best range and accuracy in estimating the 3-D surface displacement in the study area. Therefore, this paper constructs a leveling-constrained multidirectional PS-InSAR 3-D surface displacement estimation method to obtain the 3-D surface displacement field of eastern Beijing, verifies that this method can accurately reflect the evolutionary characteristics of this 3-D displacement field through leveling and GPS data, and ultimately clarifies the evolutionary and influencing mechanisms of the 3-D displacement field in combination with regional hydrogeological data. This method will provide technical support and a scientific basis for preventing and controlling regional surface displacement and its secondary hazards.
The remainder of this paper is organized as follows. Section 2 presents the background of the study area. In Section 3, the datasets and the method of processing the data are described. In Section 4, we employ the proposed leveling-constrained multidirectional PS-InSAR 3-D surface displacement estimation method to obtain the 3-D surface displacement field and verify the reliability of the InSAR-obtained results with ground level and GPS measurements. Furthermore, the spatial distribution characteristics of the surface displacement field are explored by combining spatial analysis with the estimation results of the 3-D surface displacement. In Section 5, we discuss the correlations between the 3-D surface displacement and the groundwater level and geological formations. Section 6 summarizes our conclusions.

2. Study Area and Datasets

2.1. Study Area

The study area (Figure 1) is located in the eastern part of the Beijing Plain within the warm-temperate zone with a semi-humid and semi-arid continental monsoon climate and an average annual temperature of 11–12 °C. Precipitation occurs in the area with an extremely uneven temporal distribution, with summer precipitation accounting for approximately 70% of the annual precipitation.
The Beijing Plain is located at the junction of the Yanshan, North China Plain and Taihang Mountain fault systems and, therefore, features a complex geological environment and active neotectonic activity [52]. Three main faults run through the study area: the Shunyi–Liangxiang fault, the Nanyuan–Tongxian fault and the Nankou–Sunhe fault. These faults, most of which are normal and play a key role in controlling the geological development of rift basins and depressions in this area [52], strike mainly northeast (NE)–southwest (SW) and northwest (NW)–southeast (SE). The main characteristics of these faults are shown in Table 1.
The groundwater system in the study area comprises contributions from the Yongding River, Wenyu River and Chaobai River. Quaternary sediments are widely distributed throughout the plain area, and the Quaternary aquifers in the study area are divided into four main aquifer groups according to the groundwater replenishment conditions, depth, drainage conditions, groundwater extraction levels and genesis. The lithology of each aquifer is listed in Table 2.

2.2. Dataset

The SAR dataset used in this paper consists of lift-orbit images acquired by Sentinel-1A/B, two Earth observation satellites in the Copernicus program (formerly known as the Global Monitoring for Environment and Security (GMES) program) of the European Space Agency that carry C-band synthetic aperture radars with a revisit period of 6 days. This study utilizes 44 ascending and 66 descending images (Tracks 142 and 47) acquired under vertical-vertical (VV) polarization in interferometric wide swath (IW) mode between 2016 and 2018. The spatial resolution of Sentinel-1 in IW Terrain Observation with Progressive Scans (TOPS) mode is 5 m × 20 m with a width of 250 km (https://sentinel.esa.int/, accessed on 16 July 2020). The detailed parameters of the SAR data are shown in Table 3.

3. Principle of Retrieving 3-D Displacement Data from Multidirectional PS-InSAR Measurements and Leveling Data

Direct observations of the surface from space do not represent the real surface displacement but are instead the sum of the projections of the surface displacement in the due east, due north and vertical directions in the radar LOS. Therefore, the displacement results obtained from InSAR measurements can be decomposed into the components of the above three directions according to the spatial relationship of the observation geometry. The LOS-oriented displacement observed by InSAR can be expressed as [55]:
d L O S = d u · cos θ i + d n · sin θ i sin α i 3 π / 2 d e · sin θ i cos α i 3 π / 2 + C
where d L O S , d n , d e , and d u denote the LOS direction and the displacement in the due north, due east, and vertical directions, respectively; θ is the central incidence angle; α is the angle (measured clockwise) between the satellite heading and the due north direction; C, representing the systematic error of the LOS direction result, is a constant caused by the displacement of the selected reference point not being zero.
There is a large difference in the vertical displacement rates between the two sides of a ground fissure found in an area with active faulting, and the accompanying horizontal displacement is not negligible. Assuming that the rate of change of vertical displacement is positively correlated with the horizontal displacement rate and considering that the slope can be used for each image to calculate the maximum rate of change of the value in the direction from that image to its neighboring image, the slope of the vertical displacement rate is used to characterize its rate of change and thus to assess the magnitude of horizontal displacement.
The slope depends on the rate of change (increment) of the surface in the horizontal ( d z / d x ) and vertical ( d z / d y ) directions starting from the central image element. The basic algorithm used to calculate the slope takes the values of the central image element and its eight adjacent images to determine the horizontal and vertical increments. These adjacent elements are determined using the letters a through i (as shown in Figure 2), where e indicates the element for which the slope is currently being calculated.
Then, the rates of change of the image element e in the x and y directions can be calculated by the following equations [56]:
d z d x = c + 2 f + i a + 2 d + g / 8 * c e l l s i z e
d z d y = g + 2 h + i a + 2 b + c / 8 * c e l l s i z e
The slope of the central image element e can be calculated by substituting the rate of change in the x and y directions:
s l o p e d e g r e e s = tan 1 d z d x 2 + d z d y 2 × 180 π
According to Equation (4), the slope of the whole area can be calculated by combining the vertical displacement rate of the study area obtained by PS-InSAR. The approximate vertical displacement rate of the study area can be calculated by the following equation while ignoring the horizontal (north–south and east–west) components [57]:
d u d L O S ÷ cos θ i
The kriging interpolation method is chosen to interpolate the point-like vertical displacement rates onto a regular grid. According to the calculated slopes of the vertical displacement rate, the leveling measurements from the area where horizontal displacement does not easily occur are selected to remove the systematic error C .
The ground leveling measurements screened in the above way are chosen to remove the systematic error of the LOS results. Since the horizontal displacement at the screened level can be neglected, Equation (1) can be simplified as:
d L O S = d u · cos θ i + C
According to Equation (6), the leveling results are projected upward in the LOS direction, the average error of the InSAR results in each direction is calculated, and the average error is considered to be the systematic error caused by the selection of the reference point. Then, the systematic error is removed for the InSAR results from the ascending and descending tracks.
From Equation (1), if we want to solve for the 3-D displacement field by using the InSAR results, even if the systematic error C is removed, we still need at least three different InSAR LOS directions [58]:
d u d e d n = Γ 1 Γ 2 Γ 3 d L O S , 1 d L O S , 2 d L O S , 3 =   a 1     b 1     c 1   a 2     b 2     c 2 a 3     b 3     c 3 1 d L O S , 1 d L O S , 2 d L O S , 3
where:
a i = cos θ i   ,   i = 1 , 2 , 3
b i = sin θ i sin ( α i 3 π / 2 )   ,   i = 1 , 2 , 3
c i = sin θ i cos ( α i 3 π / 2 )   ,   i = 1 , 2 , 3
The previously obtained regular grid results, the central incidence angle, and the azimuth angle in each direction are input into Equation (8) grid by grid for the calculation, and all the grids are traversed to obtain the vertical and east–west horizontal displacement components in the study area with high precision. The main processing steps are shown in Figure 3.

4. Results and Analysis

4.1. Obtaining the LOS Displacements by Implementing PS-InSAR

To measure the time-series surface displacement in the eastern Beijing Plain, the lift-track SAR data were processed separately by implementing PS-InSAR using SARProZ processing software (refer to https://www.sarproz.com/, accessed on 24 September 2018). The main processing steps are as follows: master image selection, SAR data alignment, digital elevation model (DEM) simulation, interferogram generation, PS point selection using an amplitude stability index (ASI) threshold, atmospheric phase (APS) estimation and removal, and temporal surface displacement estimation. Among them, the DEM used for removing the terrain phase is 90 m resolution Shuttle Radar Topography Mission (SRTM) (http://dds.cr.usgs.gov/srtm/, accessed on 28 September 2018) data. In addition, only PS points with an ASI greater than 0.75 were used to ensure high coherence and stability.
A total of 1,066,192 PS points were identified from 110 SAR images by the PS-InSAR technique. The density of PS points in the study area is approximately 450.77 points/km2. The left and right plots of Figure 4 present the LOS-oriented displacement rates of the ascending and descending tracks, respectively; negative values indicate subsidence, while positive values indicate uplift. The trends of the two LOS-oriented displacement results are in good agreement; the differences in the values may be due to different reference points, random errors, and different central incidence angles.
The InSAR results also reveal the heterogeneity of the spatial distribution of surface displacement in this region. As shown in Figure 4b, the areas with greater displacement are located mainly in southwestern Shunyi District, eastern Chaoyang District, and northern Tongzhou District. In addition, the maximum LOS displacement rate of −120 mm/year is observed in the Chaoyang subsidence area, and the tendency of the boundaries of several subsidence areas, such as the Chaoyang and Tongzhou subsidence areas, to join together is noted.
The PS-InSAR LOS results of the descending orbit were substituted into Equation (5), following which we obtained the vertical displacement rate field with a resolution of 50 m × 50 m by using kriging interpolation. In addition, Equations (2)–(4) were applied to calculate the slope of the vertical displacement rate throughout the study area (Figure 4). As shown in Figure 5, the slope is low in most of the study area. Comparing the spatial distribution of the slope with the LOS-oriented displacement map in Figure 4 reveals that the slope is relatively high at the margins of the subsidence areas and that the magnitudes of the slope are positively correlated with the subsidence rates of the corresponding subsidence centers. The frequency distribution of the subsidence slopes in the study area was statistically analyzed. Based on a priori knowledge, the horizontal displacement in the study area is not negligible where the slope is in the top 20% (i.e., where the slope exceeds 2.45 × 10−3). Since the subsidence slope is positively correlated with the horizontal displacement rate, the subsidence slope values at the leveling points were extracted, and the leveling points located at slope values lower than 2.45 × 10−3 were selected for the following accuracy evaluation and systematic error removal.
To evaluate the accuracy of the InSAR results, the InSAR-obtained displacement results were compared with 11 ground leveling measurements. The mean value of the PS point measurements near each benchmark (distance of less than 100 m) was calculated as the corresponding InSAR measurement result assuming no significant change in the ground elevation within this radius. The displacement information acquired by PS-InSAR is along the SAR observation direction (LOS direction). Thus, through previous studies and a priori knowledge, the SAR results in the study area with a small displacement slope were assumed to be projected in the vertical direction with respect to the corresponding incidence angle, and the leveling points were screened using a map of the LOS direction displacement slope. The leveling measurements were then projected onto the two LOS directions to verify the InSAR results. The comparison results (Figure 6a,b) show that the InSAR and leveling measurement results are basically consistent, but some systematic errors are observed.
Further statistical calculations (Table 4) show that the absolute errors of the ascending and descending results are within 10 mm, the average error is within 3 mm, and the root-mean-square-error (RMSE) is within 6 mm, indicating that the results are reliable. The average error was used to remove the systematic error from the LOS results of the ascending and descending orbits to carry out the next step of estimating the 3-D surface displacement field.

4.2. Retrieving the Vertical and East–West Displacement Components by the Proposed Approach

SAR satellites usually fly in polar orbits that are nearly parallel to the north–south direction, so they are less sensitive to the horizontal displacement in the north–south direction than to that in the vertical and east–west directions. Moreover, the horizontal displacement in the north–south direction is not the main displacement component in the study area; thus, if the north–south displacement is substituted into the calculation, it will reduce the impacts of the results in the other two directions. Consequently, the displacement in the north–south direction is ignored in the proposed method to improve the calculation accuracy in the other two directions. It follows that only the displacement results in the two LOS directions are needed for the proposed method. Then, Equation (7) can be simplified as follows:
d u d e = a 1   b 1 a 2   b 2 1 d L O S , 1 d L O S , 2
With the proposed leveling-constrained multidirectional PS-InSAR 3-D surface displacement estimation method, we decomposed the leveling-constrained LOS displacement time series and velocity fields into the vertical and east–west components by combining the PS-InSAR results from the Sentinel-1 ascending and descending orbits with the observed imaging geometry. By obtaining the imaging geometry (i.e., satellite heading and central incidence angle) of the lift-orbit SAR dataset, we were able to calculate the vertical and east–west horizontal displacement fields. Figure 7a,b show the corresponding directional components. In the plot of the vertical displacement rate (Figure 7a), negative values indicate subsidence, while positive values indicate uplift. The results show that the areas with extensive linear features are mainly located in southwestern Shunyi District, eastern Chaoyang District, and northern Tongzhou District, which is in good agreement with the spatial distribution of LOS displacement, indicating that the main displacement direction in this area is vertical. The maximum vertical displacement rate of −150 mm/year occurs in the Chaoyang subsidence area.
In the east–west horizontal displacement rate map (Figure 7b), negative values indicate the horizontal movement of the ground to the west, while positive values indicate horizontal movement to the east. The displacement rate in most areas is within ±5 mm/year. However, in some areas, there is non-negligible east–west horizontal displacement (e.g., areas A,B in Figure 7). The maximum westward displacement rate in the study area is 22 mm/year, located in the eastern part of area B, and the maximum horizontal eastward displacement rate is 20 mm/year, located in the southeastern part of area B.
The accuracies of the vertical and east–west horizontal components of the 3-D displacement estimation results were evaluated using ground leveling measurement data and GPS data, respectively. A comparison between the vertical displacement field and the interpolated lift-track LOS results (Figure 8) shows generally good agreement between the InSAR and measured results. The 11th leveling point was removed in the previous section due to the large slope at this location; the accuracy of the resolved vertical displacement at this point is slightly higher than that of the PS-InSAR results. We extracted the east–west horizontal displacement rate at this point, which is −6 mm/year, which may be one of the reasons for the difference in the accuracy of the 3-D displacement results.
Further statistical analysis (Table 5) reveals that the maximum absolute error of the estimated vertical displacement in the 3-D surface displacement field is 9.46 mm, the average error is 1.24 mm, and the RMSE is 4.29 mm. Likewise, the maximum absolute error of the east–west horizontal displacement is 4.65 mm, the average error is 2.88 mm, and the RMSE is 3.40 mm. Moreover, the accuracy in the vertical direction is higher than that in the east–west direction. The improvement of the vertical displacement accuracy may be due to the removal of systematic errors. In addition, there is generally good agreement between the east–west horizontal displacement data and the GPS results, but some systematic errors remain. We believe that these systematic errors may be due to the measurement range: the InSAR results are averaged over a 50 × 50 m range, while the GPS results are point measurements.

5. Discussion

In this section, we obtained deformation information every six months and made several profiles through ground fractures and faults. The spatial and temporal evolution of surface deformation in the eastern plain of Beijing is investigated. We also discuss the relationship between surface deformation and groundwater level, fractures and fractures.

5.1. Time Series Evolution of Surface Displacement in the Eastern Beijing Plain

Figure 9 and Figure 10 present the evolution of the vertical cumulative displacement, the east–west horizontal cumulative displacement spanning the period from January 2016 to September 2018. To observe the development of regional displacement, the annual average subsidence rate contours are plotted at 50 mm/year intervals. During the research period, the vertical deformation and the east–west horizontal deformation were gradually developing. The maximum value of the eastward horizontal displacement appeared on the west side of the area where the vertical displacement was less than −100 mm, while the maximum value of the westward horizontal displacement appeared on the east side of the area (also where the vertical displacement was less than −100 mm). In other words, the horizontal displacement was spatially corresponded to the extreme values of the corresponding vertical displacement areas, and the displacement direction was oriented toward the center of those corresponding vertical displacement areas. Comparing the vertical and horizontal displacement time series results reveals that the cumulative horizontal displacement increased with increasing vertical displacement in the corresponding displacement zone, i.e., they were positively correlated.
To further investigate the spatial characteristics of the east–west horizontal displacement, profile a–a′ was established through the Shunyi–Liangxiang fault in area A (Figure 7), and profile b–b′ was constructed through the subsidence center in area B. Figure 11 shows magnified views of the surface displacement in these two areas (marked with black rectangles in Figure 7).
The vertical displacement rates and east–west horizontal displacement rates on profiles a-a′ and b-b′ were obtained, and the slope of the vertical displacement rate and the absolute value of the east–west horizontal displacement rate were calculated. The results are shown in Figure 11. Figure 11e demonstrates that the east–west horizontal displacement rate was within ±4 mm/year in most areas of profile a–a′, but at f1, the east–west horizontal displacement rate reached −8 mm/year, and the vertical displacement rate at f1 abruptly changed. From profile b–b′, we found that the trends of the vertical and east–west horizontal displacement along the profile direction are similar, and the maximum values all appear near I. Additionally, comparing the displacement results at I and II suggests that even if the vertical displacement rates vary widely among the different areas, when the slope of the vertical displacement rate changes to a similar extent, the east–west horizontal displacement exhibits a similar response. However, in the area between III and IV, the east–west horizontal displacement rate does not change significantly even though the slope of the vertical displacement rate increases significantly. According to Figure 11, profile b–b′ is oriented approximately east–west, and no significant vertical displacement occurs in the area along profile b-b′ between III and IV, but Figure 11c displays a region with large vertical displacement in the southern part of the area between III and IV, and in the north–south. Therefore, it is presumed that the slope change in this section corresponds mainly to the horizontal displacement in the north–south direction; that is, the horizontal displacement in the area between III and IV is mainly oriented north–south.
According to the right-hand panels of Figure 11e, the displacement at I and II is eastward, while that at III and IV is westward. The area between II and III exhibits the largest vertical displacement rates on profile b–b′, so the east–west horizontal displacement direction tends to point toward the corresponding displacement center of the vertical displacement zone. Moreover, by analyzing the displacement information along both profiles a–a′ and b–b′, the area closest to the vertical displacement center displays extrusion strain, while tensile strain characterizes the areas far from the vertical displacement center. Liu et al. confirmed this finding by conducting indoor simulated pumping experiments [59]. In the horizontal displacement zone, the east–west horizontal displacement rate increases with increasing distance from the center of the corresponding vertical displacement area, and the horizontal displacement rate reaches its maximum value in the area with the largest vertical displacement slope, but the horizontal displacement rate gradually decreases with distance from this area.

5.2. Correlation between the Surface Displacement and Groundwater Level

As shown in Figure 7a, Jinzhan (area B) suffered the most severe cumulative subsidence in the Beijing Plain, and in recent years, the surface displacement rate in the Jinzhan area has continued to exceed 100 mm/year. The extraction and utilization of groundwater has been gradually reduced since the introduction of water from the south into Beijing at the end of 2014; hence, the proportion of the groundwater supply to the annual water supply decreased from 52% in 2014 to 45% in 2016. The introduction of southern water into Beijing and the replacement of private wells have helped to alleviate the excessive extraction of groundwater and have slowed the continued development of ground subsidence to a certain extent. However, the rate of ground subsidence in this region remains high. The spatial relationship between the groundwater level and vertical displacement in 2018 is shown in Figure 12. In the four subplots (a, b, c, and d), the contours indicate the water levels of the unconfined aquifer, the first confined aquifer, the second confined aquifer, and the third confined aquifer, respectively. Figure 12 demonstrates that the groundwater levels in the second and third confined aquifers are relatively low compared to those of the unconfined aquifer and the first confined aquifer; moreover, the displacement zones are distributed in the vicinity of the groundwater sink, and thus, the groundwater being exploited in this region originates mainly from the second and third confined aquifers [6].
Currently, the land use types in the Jinzhan area are mainly urban land and urban greenery. Although groundwater is no longer being overexploited following the replacement of private wells and the relocation and vacating of factories [60], environmental applications and residents living in this area continue to rely on southern water and tap water. Nevertheless, the groundwater in the area, especially the deep groundwater, remains at low levels because it is more difficult for deep sources of groundwater to recharge and recover. The overexploitation of groundwater has led to a decrease in groundwater level, a decrease in pore water pressure and an increase in the effective stresses in the aquifers. As deep groundwater cannot be recharged quickly, the aquifer skeleton becomes compressed with continued extraction, which triggers surface displacement.
Figure 13 shows the spatial relationships between the groundwater level and horizontal displacement in the east–west direction. In the four subplots (a, b, c, and d), the contours represent the water levels of the unconfined aquifer, the first confined aquifer, the second confined aquifer, and the third confined aquifer, respectively. Most of the horizontal displacement is distributed at the edge of the groundwater sink, and the shape of the horizontal displacement zone is most similar to the outlines of the second and third confined aquifer contours. However, comparing Figure 1 and Figure 12 indicates that the spatial pattern of the displacement zone matches the distribution of the compressible layer thickness well. Although the north-northeast-trending Shunyi–Liangxiang fault and northwest-trending Nankou–Sunhe fault are located near the area, the spatial distribution of the displacement zone does not follow the obvious controls of these faults, which indicates that the surface displacement in the area is less likely to be influenced by fault activity.
The Jinzhan area has a thick compressible layer, which, together with the overextraction of groundwater, has led to the formation of an elliptical vertical displacement field. In addition, we found a symmetrical east–west horizontal displacement zone relative to the center of subsidence, which may be due to the strong dynamic water pressure that formed in the seepage direction during the migration of groundwater from the periphery to the subsidence center under continuous groundwater extraction; moreover, the kinetic energy generated by the dynamic water pressure on solid particles produces an obvious viscous drag effect on the aquifer skeleton [59]. The accumulation of this viscous drag effect throughout the aquifer creates an area of concentrated tensile strain within the overlying soil layer. In addition, a drop in the water table makes the water-bearing formations deform vertically, resulting in the 3-D displacement of the overlying clay layer.
To further discuss the relationship between the groundwater level and surface displacement in each aquifer, profile e–e′ was established through the displacement center, and its location is shown in Figure 12. Figure 14 plots the vertical displacement rates and east–west horizontal displacement rates, respectively, versus the water level in each layer. The morphology of the groundwater sink in each aquifer does not exactly match the displacement area. Hence, we tentatively suggest that the groundwater level is responsible for the displacement in this area, but the spatial spread of displacement may be the result of the joint action of groundwater level variations and other factors.

5.3. Correlations between the 3-D Surface Displacement and Fissures and Active Faults

To discuss the relationships between the active faults and 3-D displacement of the surface, we superimposed the surface displacement with the surface traces of the main active faults and ground fissures that pass through the study area (Figure 7). Obvious differential displacement occurs on both sides of some of the faults and ground fissures, and the shapes of these displacement zones are highly consistent with the trends of the faults, which indicates that the spatial distribution of surface displacement is controlled by the faults and ground fissures to some extent in the eastern part of the Beijing Plain. To better understand the differences in displacement between the two sides of the faults and ground fissures, we plotted profiles across three vertical faults (the left panel of Figure 11e), as shown in Figure 15. In the vicinity of these ground fissures and faults, we detected obvious changes in the slope of vertical displacement, and tensile strain was discovered in the east–west direction. On some parts of the faults and ground fissures, we noted obvious differences in displacement on both sides, and the shapes of these displacement areas are highly consistent with the trends of these linear structures, which indicates that the spatial distribution of surface displacement is controlled by faults and ground fissures to some extent in the eastern Beijing Plain.
As shown in Figure 11, there are obvious differences between the ground subsidence rates and east–west horizontal displacement rates in the areas northwest and southeast of the airport; specifically, the ground subsidence rate to the northwest is larger than that to the southeast, and both sides are deformed horizontally to the west; however, extrusion strain occurs on the northwest side, whereas tensile strain occurs on the southeast side. The boundary between these two sides is essentially along the Shunyi–Liangxiang fault, indicating a correlation between the surface displacement and the fissures and active faults in this area. Such differential subsidence has continued for many years, and obvious cumulative effects can already be seen, as they have triggered the appearance of ground fissures [57].
In Section 5.2, the relationships between the groundwater level and surface displacement are discussed, and it is concluded that the overexploitation of groundwater is the main driving force for the surface displacement in the eastern Beijing Plain. This inhomogeneous surface displacement is closely related to the occurrence and development of ground fissures in the eastern Beijing Plain. Taking f1 (the ground fissure near the airport in Shunyi District) in area A as an example, the long-term overexploitation of the confined aquifer approximately 70–200 m below the surface in Shunyi District has led to a decrease in the groundwater level, a decrease in pore water pressure, and an increase in the effective stress in the aquifer. As deep groundwater cannot be recharged quickly, the aquifer skeleton becomes compressed, triggering surface displacement. The Nanfashin displacement center to the north and the Houshayu displacement center to the west also experienced obvious surface displacement and were significantly affected by uneven local displacement. The period of ground fissure activity coincides with the period in which ground subsidence rapidly developed. A borehole investigation showed that the vertical slip rate along the Shunyi–Liangxiang fault is approximately 0.03 mm/year, which is much smaller than the displacement rate detected by InSAR, indicating that the activity of F1 cannot fully explain the recent ground fissure activity. Generally, ground fissures formed by uneven displacement triggered by groundwater pumping will show distribution characteristics related to the height of the groundwater level; accordingly, subsidence tends to be more severe corresponding to lower groundwater levels. However, there is no significant difference in the groundwater table on either side of f1 (Figure 12). The early mid Quaternary Shunyi–Liangxiang fault is a normal fault, and thus, thick, fine-grained, loose and highly compressible sediments formed on the hanging wall of this fault; furthermore, the Quaternary sediment thickness varies greatly on both sides of the fault. In this case, the fault controls the spatial distribution of differential displacement, and the differential displacement between the foot and hanging walls leads to the formation of ground fissures.
The left panels of Figure 11e show that the slopes of vertical displacement and the rates of east–west horizontal displacement increase closer to the Shunyi ground fissure in the vicinity of the airport, and the slopes and rates reach their peak at the main fissure. The peak horizontal displacement is −10 mm/year along profile a-a′, and the maximum differences between the two sides of the fissure are 28 mm/year in the vertical displacement rate and 6/year in the east–west horizontal displacement rate, with the horizontal displacement oriented toward the center of the vertical displacement area. Similar characteristics are also found for the Songzhuang ground fissure. According to the left panels of Figure 15 the vertical displacement slope and east–west horizontal displacement rate are larger closer to the ground fissure and reach their peak at the main fissure. The peak horizontal displacement along profile c–c′ is −7 mm/year, the maximum difference in the vertical displacement rate between the two sides of the ground fissure is 27 mm/year, and the maximum difference in the east–west horizontal displacement rate is 13 mm/year.
As shown in Figure 7, the spatial distribution of the Songzhuang fissure is consistent with that of the Nanyuan–Tongxian fault. The overall trend is 50°~30° east of north, and the fault surface tends to dip toward the northwest with an inclination angle of approximately 50°. The northwest side is the hanging wall, and the southeast side is the foot wall, which is uplifted, making this fault a normal fault. Under the influence of this fault, the sediment thicknesses and sedimentation rates of different Quaternary periods on both sides of the fault are significantly different, with the Quaternary sediment thickness on the foot wall being 553 m and that on the hanging wall being 328 m [61]. Although this fault is currently active and the ground fissure orientation is very consistent with the fault strike, comparing the InSAR-monitored displacement information with the activity rate of the fault reveals that the fault activity cannot fully explain the displacement pertaining to the ground fissure. The causes of the Shunyi fissure near the airport and the Songzhuang fissure are similar in that they were both created by the large difference in sediment thickness between the two sides of the fault exacerbated by the activity of the normal fault and by the basic reliance on the extraction of deep groundwater for industrial and agricultural purposes, resulting in the occurrence of differential displacement and the generation of fissures over many years. However, the specific characteristics regarding the origins of different ground fissures vary, as they may have been caused by differences in stratigraphic conditions, the locations of groundwater extraction wells, and/or the extraction time and extraction volume. The two sides of a ground fissure usually exhibit different surface displacement, which can change or weaken the local displacement trend. On the one hand, a large local subsidence slope disrupts the surface tension of the soil and promotes relative movement between the two sides of the ground fissure. On the other hand, the formation of ground fissures interrupts the flow of groundwater, restricts and weakens the horizontal spreading of depressions, and intensifies local subsidence. Thus, surface displacement and ground fissures promote each other.

6. Conclusions

In this paper, we take the eastern Beijing Plain as the study area and study the evolutionary characteristics and influencing factors of the regional 3-D displacement field. First, we select 110 ascending and descending track Sentinel-1A/B images and apply the PS-InSAR technique to obtain 3-D surface displacement information from 2016 to 2018. After obtaining single-track SAR-influenced displacement information using the PS-InSAR technique, we then construct a leveling-constrained multidirectional InSAR 3-D surface displacement estimation method to obtain the vertical and east–west horizontal surface displacement information in the study area from 2016 to 2018. The leveling data and GPS data are also used to verify the accuracy of the InSAR results. The results show that the maximum values of the vertical displacement, eastward horizontal displacement and westward horizontal displacement in the study period are −150 mm/year, 20 mm/year and 22 mm/year, respectively. The InSAR results are highly consistent with the leveling measurements and GPS measurements, indicating that our results are reliable. Moreover, comparing the results of the vertical displacement field before and after the estimation suggests that the accuracy of the vertical displacement field after the leveling-constrained multidirectional PS-InSAR 3-D surface displacement estimation method is higher.
A spatial overlap analysis and profile analysis were performed to analyze the evolutionary characteristics of the 3-D surface displacement time series in the study area. The results show that the surface displacement within the study area is mainly vertical displacement during the study period, but there is obvious horizontal displacement oriented toward the center of vertical displacement in the areas with large vertical displacement gradients. Hence, if the displacement is large or ground faults are active, the horizontal component of the surface displacement cannot be ignored. Extrusion strain occurs closer to the centers of vertical displacement, while tensile strain occurs far from the centers of vertical displacement. In the horizontal displacement zones, the east–west horizontal displacement rate is greater farther from the corresponding center of vertical displacement, and when the vertical displacement slope reaches its peak, the horizontal displacement rate also reaches its maximum value, after which the horizontal displacement rate gradually decreases with increasing distance.
On the basis of the above findings, the correlations between the surface displacement and the spatial distributions of the groundwater level, ground fissures and faults are discussed by combining the InSAR 3-D displacement results with regional hydrogeological data. We discover that the morphology of the groundwater sink in each aquifer is in good (but not perfect) agreement with the displacement area. The spatial elongation of surface displacement along active faults is highly consistent with the fissure direction, suggesting that faults control the displacement direction by influencing the thickness and composition of sediments on both sides of the lineament. The analysis results further demonstrate that although the fluctuation in the groundwater level is the main factor influencing the surface displacement, the spatial distribution of displacement is the result of the joint action of groundwater level variations and differences in the compressible layer thickness, fault distribution and other factors. Moreover, analyses of the regional vertical and east–west displacement characteristics and the ground fissure spreading patterns reveal that differential surface displacement easily leads to the generation of ground fissures, and the presence of ground fissures can induce local surface displacement, intensify local subsidence, and limit the horizontal spreading of groundwater sinks. Therefore, the spatial and temporal variations in surface displacement are related to the activity of ground fissures, and they promote each other.

Author Contributions

S.Z. designed the experiments, implemented the algorithm, analyzed the results, and wrote the paper. B.C. made important suggestions on writing the paper. H.G., M.S. and C.Z. provided crucial guidance. K.L. provided the data. All authors have read and agreed to the published version of the manuscript.

Funding

This work was founded by Beijing Outstanding Young Scientist Program (No. BJJWZYJH01201910028032), Beijing Youth Top Talent Project, National Natural Science Foundation of China (No. 41930109/D010702, 41771455/D010702), National “Double-Class” Construction of University Projects, Beijing Postdoctoral Research Foundation (No. 2018M641407).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Liu, G.; Zhang, R.; Li, T.; Yu, B.; Li, T.; Jia, H.; Nie, Y. Extracting 3D ground deformation velocity field by multi-platform persistent scatterer SAR interferometry. Chin. J. Geophys. 2012, 55, 2598–2610. [Google Scholar]
  2. Osmanoğlu, B.; Dixon, T.H.; Wdowinski, S.; Cabral-Cano, E.; Jiang, Y. Mexico City subsidence observed with persistent scatterer InSAR. Int. J. Appl. Earth Obs. 2011, 13, 1–12. [Google Scholar] [CrossRef]
  3. Ishitsuka, K.; Tsuji, T.; Matsuoka, T.; Nishijima, J.; Fujimitsu, Y. Heterogeneous surface displacement pattern at the Hatchobaru geothermal field inferred from SAR interferometry time-series. Int. J. Appl. Earth Obs. 2016, 44, 95–103. [Google Scholar] [CrossRef]
  4. Gumilar, I.; Abidin, H.Z.; Hutasoit, L.M.; Hakim, D.M.; Sidiq, T.P.; Andreas, H. Land Subsidence in Bandung Basin and its Possible Caused Factors. Procedia Earth Planet. Sci. 2015, 12, 47–62. [Google Scholar] [CrossRef] [Green Version]
  5. Zhang, Y.; Gong, H.; Gu, Z.; Wang, R.; Li, X.; Zhao, W. Characterization of land subsidence induced by groundwater withdrawals in the plain of Beijing city, China. Hydrogeol. J. 2014, 22, 397–409. [Google Scholar] [CrossRef]
  6. Chen, B.; Gong, H.; Chen, Y.; Li, X.; Zhou, C.; Lei, K.; Zhu, L.; Duan, L.; Zhao, X. Land subsidence and its relation with groundwater aquifers in Beijing Plain of China. Sci. Total Environ. 2020, 735, 139111. [Google Scholar] [CrossRef] [PubMed]
  7. Chen, M.; Tomás, R.; Li, Z.; Motagh, M.; Li, T.; Hu, L.; Gong, H.; Li, X.; Yu, J.; Gong, X. Imaging Land Subsidence Induced by Groundwater Extraction in Beijing (China) Using Satellite Radar Interferometry. Remote Sens. 2016, 8, 468. [Google Scholar] [CrossRef] [Green Version]
  8. Gao, M.; Gong, H.; Chen, B.; Zhou, C.; Chen, W.; Liang, Y.; Shi, M.; Si, Y. InSAR time-series investigation of long-term ground displacement at Beijing Capital International Airport, China. Tectonophysics 2016, 691, 271–281. [Google Scholar] [CrossRef]
  9. Guo, J.; Zhou, L.; Yao, C.; Hu, J. Surface Subsidence Analysis by Multi-Temporal InSAR and GRACE: A Case Study in Beijing. Sensors 2016, 16, 1495. [Google Scholar] [CrossRef] [Green Version]
  10. Zhu, L.; Gong, H.; Li, X.; Wang, R.; Chen, B.; Dai, Z.; Teatini, P. Land subsidence due to groundwater withdrawal in the northern Beijing plain, China. Eng. Geol. 2015, 193, 243–255. [Google Scholar] [CrossRef]
  11. Ng, A.H.; Ng, A.H.; Ge, L.; Ge, L.; Li, X.; Li, X.; Zhang, K.; Zhang, K. Monitoring ground deformation in Beijing, China with persistent scatterer SAR interferometry. J. Geodesy. 2012, 86, 375–392. [Google Scholar] [CrossRef]
  12. Zeng, D.; Ke, Y.; Gong, H.; Chen, B.; Zhu, L. Preliminary research on land subsidence prediction method in Beijing. In Proceedings of the 2016 IEEE International Geoscience and Remote Sensing Symposium (IGARSS), Beijing, China, 10–15 July 2016; pp. 5951–5954. [Google Scholar] [CrossRef]
  13. Chen, B.; Gong, H.; Li, X.; Lei, K.; Zhang, Y.; Li, J.; Gu, Z.; Dang, Y. Spatial-temporal characteristics of land subsidence corresponding to dynamic groundwater funnel in Beijing Municipality, China. Chin. Geogr. Sci. 2011, 21, 753–764. [Google Scholar] [CrossRef]
  14. Zhou, C.; Gong, H.; Chen, B.; Guo, L.; Gao, M.; Chen, W.; Liang, Y.; Si, Y.; Wang, J.; Zhang, X. Spatiotemporal characteristics of land subsidence in Beijing from small baseline subset interferometric synthetic aperture radar and standard deviational ellipse. In Proceedings of the 2016 4th International Workshop on Earth Observation and Remote Sensing Applications (EORSA), Guangzhou, China, 4–6 July 2016. [Google Scholar] [CrossRef]
  15. Luo, Q.; Perissin, D.; Zhang, Y.; Jia, Y. L- and X-Band Multi-Temporal InSAR Analysis of Tianjin Subsidence. Remote Sens. 2014, 6, 7933–7951. [Google Scholar] [CrossRef] [Green Version]
  16. Yang, Y. Effectiveness of InSAR monitoring of land subsidence in Beijing. Shanghai Land Resour. 2013, 34, 21–24. [Google Scholar]
  17. Yin, J.; Yu, D.; Wilby, R. Modelling the impact of land subsidence on urban pluvial flooding: A case study of downtown Shanghai, China. Sci. Total Environ. 2016, 544, 744–753. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  18. Chen, J.; Wu, J.; Zhang, L.; Zou, J.; Liu, G.; Zhang, R.; Yu, B. Deformation Trend Extraction Based on Multi-Temporal InSAR in Shanghai. Remote Sens. 2013, 5, 1774–1786. [Google Scholar] [CrossRef] [Green Version]
  19. Luo, Y.; Ye, S.; Wu, J.; Wang, H.; Jiao, X. A modified inverse procedure for calibrating parameters in a land subsidence model and its field application in Shanghai, China. Hydrogeol. J. 2016, 24, 711–725. [Google Scholar] [CrossRef]
  20. Jia, S.; Wang, H.; Zhao, S.; Luo, Y. A tentative study of the mechanism of land subsidence in Beijing. Urban Geol. 2007, 2, 20–26. [Google Scholar]
  21. Luo, Y.; Jia, S.; Yang, Y.; Wang, R.; Lei, K. Prevention and control of land subsidence in Beijing. Shanghai Land Resour. 2014, 35, 110–113, 141. [Google Scholar]
  22. Rossi, M.; Peltzer, G.; Adragna, F.; Carmona, C.; Massonnet, D.; Feigl, K.; Rabaute, T. The displacement field of the Landers earthquake mapped by radar interferometry. Nature 1993, 364, 138–142. [Google Scholar] [CrossRef]
  23. Lu, Z.; JR, C.W.; Power, J.; Dzurisin, D.; Masterlark, T. Interferometric synthetic aperture radar studies of Alaska volcanoes. In Proceedings of the IEEE International Geoscience and Remote Sensing Symposium, Toronto, ON, Canada, 24–28 June 2002; Volume 1, pp. 191–194. [Google Scholar] [CrossRef] [Green Version]
  24. Hilley, G.E.; Bürgmann, R.; Ferretti, A.; Novali, F.; Rocca, F. Dynamics of slow-moving landslides from permanent scatterer analysis. Science 2004, 304, 1952–1955. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  25. Goldstein, R.M.; Engelhardt, H.; Kamb, B.; Frolich, R.M. Satellite radar interferometry for monitoring ice sheet motion: Application to an antarctic ice stream. Science 1993, 262, 1525–1530. [Google Scholar] [CrossRef] [PubMed]
  26. Hudnut, K.W.; Peltzer, G.; Bawden, G.W.; Stein, R.S.; Thatcher, W. Tectonic contraction across Los Angeles after removal of groundwater pumping effects. Nature 2001, 412, 812–815. [Google Scholar] [CrossRef]
  27. Bell, J.W.; Amelung, F.; Ferretti, A.; Bianchi, M.; Novali, F. Permanent scatterer InSAR reveals seasonal and long-term aquifer-system response to groundwater pumping and artificial recharge. Water Resour. Res. 2008, 44. [Google Scholar] [CrossRef] [Green Version]
  28. Carnec, C.; Delacourt, C. Three years of mining subsidence monitored by SAR interferometry, near Gardanne, France. J. Appl. Geophys. 2000, 43, 43–54. [Google Scholar] [CrossRef]
  29. Fan, H.; Wang, L.; Wen, B.; Du, S. A New Model for three-dimensional Deformation Extraction with Single-track InSAR Based on Mining Subsidence Characteristics. Int. J. Appl. Earth Obs. 2021, 94, 102223. [Google Scholar] [CrossRef]
  30. Wright, T.J.; Parsons, B.E.; Lu, Z. Toward mapping surface deformation in three dimensions using InSAR. Geophys. Res. Lett. 2004, 31. [Google Scholar] [CrossRef] [Green Version]
  31. Gray, L. Using multiple RADARSAT InSAR pairs to estimate a full three-dimensional solution for glacial ice movement. Geophys. Res. Lett. 2011, 38, L05502. [Google Scholar] [CrossRef]
  32. Chang, Z.; Wang, Y.; Qian, S.; Zhu, J.; Wang, W.; Liu, X.; Yu, J. An approach for retrieving complete three-dimensional ground displacement components from two parallel-track InSAR measurements. J. Geodesy 2020, 94, 1–2. [Google Scholar] [CrossRef]
  33. Fialko, Y.; Simons, M.; Agnew, D. The complete (3-D) surface displacement field in the epicentral area of the 1999 M(w)7.1 Hector Mine earthquake, California, from space geodetic observations. Geophys. Res. Lett. 2001, 28, 3063–3066. [Google Scholar] [CrossRef] [Green Version]
  34. Raucoules, D.; de Michele, M.; Malet, J.P.; Ulrich, P. Time-variable 3D ground displacements from high-resolution synthetic aperture radar (SAR). application to La Valette landslide (South French Alps). Remote Sens. Environ. 2013, 139, 198–204. [Google Scholar] [CrossRef] [Green Version]
  35. Jung, H.S.; Lu, Z.; Won, J.S.; Poland, M.P.; Miklius, A. Mapping Three-Dimensional Surface Deformation by Combining Multiple-Aperture Interferometry and Conventional Interferometry: Application to the June 2007 Eruption of Kilauea Volcano, Hawaii. IEEE Geosci. Remote Sens. Lett. 2011, 8, 34–38. [Google Scholar] [CrossRef]
  36. Gourmelen, N.; Kim, S.W.; Shepherd, A.; Park, J.W.; Sundal, A.V.; Björnsson, H.; Pálsson, F. Ice velocity determined using conventional and multiple-aperture InSAR. Earth Planet. Sci. Lett. 2011, 307, 156–160. [Google Scholar] [CrossRef] [Green Version]
  37. He, L.; Wu, L.; Liu, S.; Wang, Z.; Su, C.; Liu, S. Mapping Two-Dimensional Deformation Field Time-Series of Large Slope by Coupling DInSAR-SBAS with MAI-SBAS. Remote Sens. 2015, 7, 12440–12458. [Google Scholar] [CrossRef] [Green Version]
  38. Baek, W.K.; Jung, H.S. Precise Three-Dimensional Deformation Retrieval in Large and Complex Deformation Areas via Integration of Offset-Based Unwrapping and Improved Multiple-Aperture SAR Interferometry: Application to the 2016 Kumamoto Earthquake. Engineering 2020, 6, 927–935. [Google Scholar] [CrossRef]
  39. Joughin, I.R.; Kwok, R.; Fahnestock, M.A. Interferometric estimation of three-dimensional ice-flow using ascending and descending passes. IEEE Trans. Geosci. Remote Sens. 1998, 36, 25–37. [Google Scholar] [CrossRef] [Green Version]
  40. Mohr, J.J.; Madsen, S.N.; Reeh, N. Three-dimensional glacial flow and surface elevation measured with radar interferometry. Nature 1998, 391, 273–276. [Google Scholar] [CrossRef]
  41. Gudmundsson, S.; Sigmundsson, F.; Carstensen, J.M. Three-dimensional surface motion maps estimated from combined interferometric synthetic aperture radar and GPS data. J. Geophys. Res. Solid Earth 2002, 107, 11–13. [Google Scholar] [CrossRef]
  42. Samsonov, S.; Tiampo, K. Analytical Optimization of a DInSAR and GPS Dataset for Derivation of Three-Dimensional Surface Motion. IEEE Geosci. Remote Sens. Lett. 2006, 3, 107–111. [Google Scholar] [CrossRef]
  43. Guglielmino, F.; Nunnari, G.; Puglisi, G.; Spata, A. Simultaneous and Integrated Strain Tensor Estimation From Geodetic and Satellite Deformation Measurements to Obtain Three-Dimensional Displacement Maps. IEEE Trans. Geosci. Remote Sens. 2011, 49, 1815–1826. [Google Scholar] [CrossRef]
  44. Vollrath, A.; Zucca, F.; Bekaert, D.; Bonforte, A.; Guglielmino, F.; Hooper, A.J.; Stramondo, S. Decomposing DInSAR Time-Series into 3-D in Combination with GPS in the Case of Low Strain Rates: An Application to the Hyblean Plateau, Sicily, Italy. Remote Sens. 2017, 9, 33. [Google Scholar] [CrossRef] [Green Version]
  45. Catalao, J.; Nico, G.; Hanssen, R.; Catita, C. Merging GPS and Atmospherically Corrected InSAR Data to Map 3-D Terrain Displacement Velocity. IEEE Trans. Geosci. Remote Sens. 2011, 49, 2354–2360. [Google Scholar] [CrossRef] [Green Version]
  46. Liu, X.; Hu, J.; Sun, Q.; Li, Z.; Zhu, J. Deriving 3-D Time-Series Ground Deformations Induced by Underground Fluid Flows with InSAR: Case Study of Sebei Gas Fields, China. Remote Sens. 2017, 9, 1129. [Google Scholar] [CrossRef] [Green Version]
  47. Qu, F.; Zhang, Q.; Lu, Z.; Zhao, C.; Yang, C.; Zhang, J. Land subsidence and ground fissures in Xi’an, China 2005–2012 revealed by multi-band InSAR time-series analysis. Remote Sens. Environ. 2014, 155, 366–376. [Google Scholar] [CrossRef]
  48. Liu, Y.; Zhao, C.; Zhang, Q.; Yang, C. Complex surface deformation monitoring and mechanism inversion over Qingxu-Jiaocheng, China with multi-sensor SAR images. J. Geodyn. 2018, 114, 41–52. [Google Scholar] [CrossRef]
  49. Yang, C.; Lu, Z.; Zhang, Q.; Zhao, C.; Peng, J.; Ji, L. Deformation at longyao ground fissure and its surroundings, north China plain, revealed by ALOS PALSAR PS-InSAR. Int. J. Appl. Earth Obs. 2018, 67, 1–9. [Google Scholar] [CrossRef]
  50. Hu, X.; Oommen, T.; Lu, Z.; Wang, T.; Kim, J. Consolidation settlement of Salt Lake County tailings impoundment revealed by time-series InSAR observations from multiple radar satellites. Remote Sens. Environ. 2017, 202, 199–209. [Google Scholar] [CrossRef]
  51. Chaussard, E.; Milillo, P.; Bürgmann, R.; Perissin, D.; Fielding, E.J.; Baker, B. Remote Sensing of Ground Deformation for Monitoring Groundwater Management Practices: Application to the Santa Clara Valley During the 2012–2015 California Drought. J. Geophys. Res. Solid Earth 2017, 122, 8566–8582. [Google Scholar] [CrossRef]
  52. Hu, L.; Dai, K.; Xing, C.; Li, Z.; Lu, Y. Land subsidence in Beijing and its relationship with geological faults revealed by Sentinel-1 InSAR observations. Int. J. Appl. Earth Obs. 2019, 82. [Google Scholar] [CrossRef]
  53. Xu, X.; Wu, W.; Zhang, X.; Ma, S.; Ma, W.; Gu, M.; Jiang, W. Crust Newly Tectonic Deformation and Earthquake in Captial Zone; Science Press: Beijing, China, 2002. [Google Scholar]
  54. BBG (Beijing Bureau of Geology and Mineral Exploration and Development). Groundwater in Beijing; Land Press of China: Beijing, China, 2008. [Google Scholar]
  55. Zhang, R. Modeling and Deformation Estimating with Multi-Platform Persistent Scatterer Radar Interferometry Based on Multi-Level Networking. Ph.D. Thesis, Southwest Jiaotong University, Chengdu, China, 2012. [Google Scholar]
  56. Burrough, P.A.; McDonell, R.A. Principles of Geographical Information Systems; Oxford University Press: New York, NY, USA, 1998; p. 190. [Google Scholar]
  57. Gao, M.; Gong, H.; Li, X.; Chen, B.; Zhou, C.; Shi, M.; Guo, L.; Chen, Z.; Ni, Z.; Duan, G. Land Subsidence and Ground Fissures in Beijing Capital International Airport (BCIA): Evidence from Quasi-PS InSAR Analysis. Remote Sens. 2019, 11, 1466. [Google Scholar] [CrossRef] [Green Version]
  58. Hu, J. Theory and Method of Estimating Three-Dimensional Displacement with InSAR Based on the Modern Surveying Adjustment. Ph.D. Thesis, Central South University, Changsha, China, 2013. [Google Scholar]
  59. Liu, H. Horizontal Aquifer Movement Induced from Ground Water Pumping and Its Relation with Ground Fissure. Master’s Thesis, Chang’an University, Xi’an, China, 31 May 2007. [Google Scholar]
  60. Ye, C.; Tian, F.; Luo, Y.; Wang, X.; Tian, M.; Cui, W.; Wang, L.; Lei, K. Control zoning and prevention measures on land subsidence in Beijing. J. Nanjing Univ. 2019, 55, 440–448. [Google Scholar]
  61. Zhao, L.; Li, Y.; Cui, W.; Luo, Y.; Zhang, Y.; Tian, F.; Lei, K.; Qiao, L.; Sha, T.; Tian, M.; et al. Disaster characteristics and influence factors for ground fissure at Songzhuang village in Beijing. J. Eng. Geol. 2018, 26, 1600–1610. [Google Scholar]
Figure 1. Geographical location of the study area. (a) Sentinel-1 SAR image coverage. (b) Sketched map of the study area. F1: Shunyi–Liangxiang fault. F2: Nanyuan–Tongxian fault. F3:Nankou–Sunhe fault.
Figure 1. Geographical location of the study area. (a) Sentinel-1 SAR image coverage. (b) Sketched map of the study area. F1: Shunyi–Liangxiang fault. F2: Nanyuan–Tongxian fault. F3:Nankou–Sunhe fault.
Remotesensing 13 02809 g001
Figure 2. Diagram of the numbering of adjacent image elements.
Figure 2. Diagram of the numbering of adjacent image elements.
Remotesensing 13 02809 g002
Figure 3. Framework of the 3-D surface displacement estimation method. (LOS: Line of sight.).
Figure 3. Framework of the 3-D surface displacement estimation method. (LOS: Line of sight.).
Remotesensing 13 02809 g003
Figure 4. Annual LOS displacement maps from Sentinel-1. (a) PS-InSAR result of the ascending dataset; (b) PS-InSAR result of the descending dataset.
Figure 4. Annual LOS displacement maps from Sentinel-1. (a) PS-InSAR result of the ascending dataset; (b) PS-InSAR result of the descending dataset.
Remotesensing 13 02809 g004
Figure 5. Slope of the vertical displacement rate. 1–11 in the image is the order of benchmark.
Figure 5. Slope of the vertical displacement rate. 1–11 in the image is the order of benchmark.
Remotesensing 13 02809 g005
Figure 6. Comparison of the InSAR-derived surface displacement rates with those derived from GPS and leveling measurements. (a,b) Comparison of the ascending and descending InSAR LOS displacements, respectively, with the leveling measurements.
Figure 6. Comparison of the InSAR-derived surface displacement rates with those derived from GPS and leveling measurements. (a,b) Comparison of the ascending and descending InSAR LOS displacements, respectively, with the leveling measurements.
Remotesensing 13 02809 g006
Figure 7. Surface displacement maps from 2016 to 2018 by using ascending and descending SAR data: (a) vertical component of the displacement field; (b) east–west component of the displacement field. (f1: Shunyi fissure; f2: Songzhuang fissure. A: Beijing Capital International Airport (BCIA); B: Jinzhan. F1: Shunyi–Liangxiang fault. F2: Nanyuan–Tongxian fault. F3:Nankou–Sunhe fault. a–a’, b–b’, c–c’, d–d’ are the displacement profiles within depressions are marked by black dashed lines).
Figure 7. Surface displacement maps from 2016 to 2018 by using ascending and descending SAR data: (a) vertical component of the displacement field; (b) east–west component of the displacement field. (f1: Shunyi fissure; f2: Songzhuang fissure. A: Beijing Capital International Airport (BCIA); B: Jinzhan. F1: Shunyi–Liangxiang fault. F2: Nanyuan–Tongxian fault. F3:Nankou–Sunhe fault. a–a’, b–b’, c–c’, d–d’ are the displacement profiles within depressions are marked by black dashed lines).
Remotesensing 13 02809 g007
Figure 8. Comparison of the InSAR-derived surface displacement rates with those from GPS and leveling measurements. (a,b) Comparison of the ascending and descending InSAR LOS displacements, respectively, with leveling measurements from 2018. (c) Comparison of the vertical displacement rates with leveling measurements from 2018. (d) Comparison of the east–west displacement rates with GPS measurements from 2017.
Figure 8. Comparison of the InSAR-derived surface displacement rates with those from GPS and leveling measurements. (a,b) Comparison of the ascending and descending InSAR LOS displacements, respectively, with leveling measurements from 2018. (c) Comparison of the vertical displacement rates with leveling measurements from 2018. (d) Comparison of the east–west displacement rates with GPS measurements from 2017.
Remotesensing 13 02809 g008
Figure 9. Time series maps showing the vertical cumulative displacement. (a: Extremum of west horizontal displacement; b: Extremum of east horizontal displacement).
Figure 9. Time series maps showing the vertical cumulative displacement. (a: Extremum of west horizontal displacement; b: Extremum of east horizontal displacement).
Remotesensing 13 02809 g009
Figure 10. Time series maps of the east–west cumulative displacement. (a: Extremum of west horizontal displacement; b: Extremum of east horizontal displacement).
Figure 10. Time series maps of the east–west cumulative displacement. (a: Extremum of west horizontal displacement; b: Extremum of east horizontal displacement).
Remotesensing 13 02809 g010
Figure 11. Magnified views of the highlighted areas. (a,b) Beijing Capital International Airport (BCIA); (c,d) Jinzhan. (e) shows vertical and east–west horizontal displacement rates (2016–2018) and displacement rate gradients along the two profiles whose positions are marked in Figure 7. I, II, III, and IV are representative points on profile b-b′. The violet dashed lines indicate the locations of ground fissures and faults. F1: Shunyi–Liangxiang fault. F3:Nankou–Sunhe fault. a–a′, b–b′ are the displacement profiles within depressions are marked by black dashed lines The unit of the slope is 10−3°.
Figure 11. Magnified views of the highlighted areas. (a,b) Beijing Capital International Airport (BCIA); (c,d) Jinzhan. (e) shows vertical and east–west horizontal displacement rates (2016–2018) and displacement rate gradients along the two profiles whose positions are marked in Figure 7. I, II, III, and IV are representative points on profile b-b′. The violet dashed lines indicate the locations of ground fissures and faults. F1: Shunyi–Liangxiang fault. F3:Nankou–Sunhe fault. a–a′, b–b′ are the displacement profiles within depressions are marked by black dashed lines The unit of the slope is 10−3°.
Remotesensing 13 02809 g011
Figure 12. Spatial relationships between the groundwater level and vertical displacement. The contours in (ad) represent the water levels of the unconfined aquifer, the first confined aquifer, the second confined aquifer and the third confined aquifer, respectively. e–e’ is displacement profile within depressions are marked by black dashed lines.
Figure 12. Spatial relationships between the groundwater level and vertical displacement. The contours in (ad) represent the water levels of the unconfined aquifer, the first confined aquifer, the second confined aquifer and the third confined aquifer, respectively. e–e’ is displacement profile within depressions are marked by black dashed lines.
Remotesensing 13 02809 g012
Figure 13. Spatial relationships between the groundwater level and east–west displacement. The contours in (ad) represent the water levels of the unconfined aquifer, the first confined aquifer, the second confined aquifer and the third confined aquifer, respectively. e–e’ is displacement profile within depressions are marked by black dashed lines.
Figure 13. Spatial relationships between the groundwater level and east–west displacement. The contours in (ad) represent the water levels of the unconfined aquifer, the first confined aquifer, the second confined aquifer and the third confined aquifer, respectively. e–e’ is displacement profile within depressions are marked by black dashed lines.
Remotesensing 13 02809 g013
Figure 14. Relationships between the surface displacement and groundwater levels in (from top to bottom) the unconfined aquifer, the first confined aquifer, the second confined aquifer, and the third confined aquifer. Profile e–e′ is plotted in Figure 8.
Figure 14. Relationships between the surface displacement and groundwater levels in (from top to bottom) the unconfined aquifer, the first confined aquifer, the second confined aquifer, and the third confined aquifer. Profile e–e′ is plotted in Figure 8.
Remotesensing 13 02809 g014
Figure 15. Vertical and east–west horizontal displacement rates (2016–2018) and displacement rate gradients along two profiles whose positions are marked in Figure 3. I, II, III, and IV are representative points on the profile. The brown dashed lines indicate the locations of ground fissures and faults. The unit of the slope is 10−3°.
Figure 15. Vertical and east–west horizontal displacement rates (2016–2018) and displacement rate gradients along two profiles whose positions are marked in Figure 3. I, II, III, and IV are representative points on the profile. The brown dashed lines indicate the locations of ground fissures and faults. The unit of the slope is 10−3°.
Remotesensing 13 02809 g015
Table 1. Characteristics of the main faults in the eastern Beijing Plain based on Xu et al. [53]. Q2, Q3 and Q4 indicate the middle Pleistocene, late Pleistocene and Holocene, respectively. R: right-lateral, L: left-lateral. See Figure 1 for a map of the fault locations.
Table 1. Characteristics of the main faults in the eastern Beijing Plain based on Xu et al. [53]. Q2, Q3 and Q4 indicate the middle Pleistocene, late Pleistocene and Holocene, respectively. R: right-lateral, L: left-lateral. See Figure 1 for a map of the fault locations.
Fault NameShunyi–Liangxiang
Fault
Tongxian–Nanyuan
Fault
Nankou–Sunhe
Fault
IDF1F2F3
Fault strikeNNENNENW
Fault propertyNormal fault (R)Normal faultNormal fault (L)
Active timeQ3–4Q3Q4
Mean slip rate
(mm/year)
0.150.750.3
Table 2. Division of aquifers in the Beijing Plain [54].
Table 2. Division of aquifers in the Beijing Plain [54].
AquiferMajor LithologyDepth Range (m)
The unconfined aquiferSilt, silty sand, and sandy clay0–50
The first confined aquiferMultiple types of gravel, sand and clay soil50–100
The second confined aquiferMultiple types of gravel, sand and clay soil100–180
The third confined aquiferMainly sand180–300
Table 3. Satellite information for the data used in this study.
Table 3. Satellite information for the data used in this study.
SatelliteSentinel-1ASentinel-1A/B
BandC
Orbit directionAscendingDescending
Heading (°)−13.22−166.59
Incidence angle (°)43.8034.07
Track14247
PolarizationVertical-vertical
Image modeInterferometric wide swath
Number of images4466
Date range14 January 2016–12 September 20189 January 2016–11 September 2018
Table 4. Statistics for validating the PS-InSAR results.
Table 4. Statistics for validating the PS-InSAR results.
DataValidation DataAbs. Max
(mm)
Abs. Min
(mm)
Avg.
(mm)
RMSE
(mm)
Ascending. LOSBenchmark8.110.272.065.50
Descending. LOSBenchmark9.330.93−2.385.30
Table 5. Statistics for validating the leveling-constrained multidirectional PS-InSAR 3-D surface displacement estimation results.
Table 5. Statistics for validating the leveling-constrained multidirectional PS-InSAR 3-D surface displacement estimation results.
DataValidation DataAbs. Max
(mm)
Abs. Min
(mm)
Avg.
(mm)
RMSE
(mm)
Asc. LOSBenchmark14.110.272.766.54
Des. LOSBenchmark9.070.93−2.694.66
VerticalBenchmark9.460.321.244.29
East–west horizontalGPS4.650.212.883.40
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Zhang, S.; Chen, B.; Gong, H.; Lei, K.; Shi, M.; Zhou, C. Three-Dimensional Surface Displacement of the Eastern Beijing Plain, China, Using Ascending and Descending Sentinel-1A/B Images and Leveling Data. Remote Sens. 2021, 13, 2809. https://doi.org/10.3390/rs13142809

AMA Style

Zhang S, Chen B, Gong H, Lei K, Shi M, Zhou C. Three-Dimensional Surface Displacement of the Eastern Beijing Plain, China, Using Ascending and Descending Sentinel-1A/B Images and Leveling Data. Remote Sensing. 2021; 13(14):2809. https://doi.org/10.3390/rs13142809

Chicago/Turabian Style

Zhang, Shunkang, Beibei Chen, Huili Gong, Kunchao Lei, Min Shi, and Chaofan Zhou. 2021. "Three-Dimensional Surface Displacement of the Eastern Beijing Plain, China, Using Ascending and Descending Sentinel-1A/B Images and Leveling Data" Remote Sensing 13, no. 14: 2809. https://doi.org/10.3390/rs13142809

APA Style

Zhang, S., Chen, B., Gong, H., Lei, K., Shi, M., & Zhou, C. (2021). Three-Dimensional Surface Displacement of the Eastern Beijing Plain, China, Using Ascending and Descending Sentinel-1A/B Images and Leveling Data. Remote Sensing, 13(14), 2809. https://doi.org/10.3390/rs13142809

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop