Land Subsidence and Ground Fissures in Beijing Capital International Airport (BCIA): Evidence from Quasi-PS InSAR Analysis

: Land subsidence is a global environmental geological hazard caused by natural or human activities. The high spatial resolution and continuous time coverage of interferometric synthetic aperture radar (InSAR) time series analysis techniques provide data and a basis for the development of methods for the investigation and evolution mechanism study of regional land subsidence. Beijing, the capital city of China, has su ﬀ ered from land subsidence for decades since it was ﬁrst recorded in the 1950s. It was reported that uneven ground subsidence and fractures have seriously a ﬀ ected the normal operation of the Beijing Capital International Airport (BCIA) in recent years before the overhaul of the middle runway in April 2017. In this study, InSAR time series analysis was successfully used to detect the uneven local subsidence and ground ﬁssure activity that has been gradually increasing in BCIA since 2010. A multi-temporal InSAR (MT-InSAR) technique was performed on 63 TerraSAR-X / TanDem-X (TSX / TDX) images acquired between 2010 and 2017, then deformation rate maps and time series for the airport area were generated. Comparisons of deformation rate and displacement time series from InSAR and ground-leveling were carried out in order to evaluate the accuracy of the InSAR-derived measurements. After an integrated analysis of the distribution characteristics of land subsidence, previous research results, and geological data was carried out, we found and located an active ground ﬁssure. Then main cause of the ground ﬁssures was preliminarily discussed. Finally, it can be conducted that InSAR technology can be used to identify and monitor geological processes, such as land subsidence and ground ﬁssure activities, and can provide a scientiﬁc approach to better explore and study the cause and formation mechanism of regional subsidence and ground ﬁssures.


Introduction
Land subsidence occurs in over 150 countries all over the world, particularly in densely populated metropolitan areas, such as big cities in China [1][2][3], Japan [4], Italy [5], Mexico [6,7], and Indonesia [8,9], with extremely serious consequences. Beijing, one of the largest cities in the world, has suffered from land subsidence since the 1950s [10]. According to public reports from the China Geological Survey (CGS), the area with a sinking rate of >50 mm/year in Beijing plain was over 770 km 2 in 2017, with a maximum sinking rate of 157 mm/year, and the maximum cumulative subsidence in Beijing Plain reached 1864 mm in 2017. Differential land subsidence can result in ground fissures and can lead to damage to buildings, dams, roads, high-speed railways, airports [11] and other urban facilities (e.g., pipelines, flyovers, subway lines, and underground structures).
Beijing Capital International Airport (BCIA) is located in Shunyi District, northeast of downtown Beijing, and the airport's platform has been experiencing land subsidence since 2003 [12]. In particular, uneven land subsidence occurred in the northern area of the BCIA in recent years [13,14], which has seriously affected the safe operation of the airport. The airport handled more than 1700 flights and 270,000 passengers per day in 2017. Geological hazard such as land subsidence and ground fissure in an airport can lead to damage of ground facilities such as terminals, runways, and electronic equipment and can cause personal security and economic losses [11]. Therefore, monitoring the temporal and spatial evolution of long-term land subsidence is critical to implementing disaster warnings and remedial measures for such critical civil infrastructure to minimize potential damage. Considering the potential damage that ground fissures can cause to the environment, numerous techniques (e.g., geodetic survey and geophysical exploration technology) have been used to monitor and study the formation mechanism of ground fissures [15][16][17]. Geodetic methods, such as leveling and global positioning systems (GPS), tend to focus on deformation monitoring of ground fissures, while geophysical methods are commonly used to detect the location and depth of underground fractures.
Over the past few decades, interferometric synthetic aperture radar (InSAR) has been developed as an effective geodetic surveying technique and used for investigating geophysical phenomena, including seismic deformation [18,19], volcano dynamics [20,21], earthquakes [22,23], landslides [24,25], and other geologic processes. The InSAR technique can also monitor the ground deformation caused by human activities (e.g., mining [26,27], delta sinking [28], and land subsidence [29][30][31][32]). Ground fissure is commonly a secondary environmental geological problem caused by land subsidence [33] and the InSAR technique can provide a new tool for investigating ground fissure activity [34,35]. Compared with the conventional ground-based measurement methods on sparse points, such as differential GPS and leveling, InSAR can determine a planar observation on a large scale with high precision and low cost [30]. Since the revisit period of the space-borne radar system is currently less than two weeks (e.g., for TerraSAR-X/TanDEM-X and Sentinel-1A/1B, the orbit recurrent period can be shortened to 11 and 6 days, respectively), this technique can also provide sufficient data for long-term observations. However, the conventional InSAR technique is fatally limited by temporal and spatial decorrelation and atmospheric delay caused by different satellite positions, long acquisition intervals, and atmospheric fluctuations [36], respectively. In this case, the deformation monitoring accuracy provided by D-InSAR is only from decimeters to centimeters and there are still great limitations in the application of urban infrastructure monitoring. Hence, persistent scatterer interferometry (PS-InSAR) [37,38] has been developed to overcome the incoherence limitation over time intervals and large areas of the conventional InSAR technique. This advanced technique represents a specific class of multi-temporal InSAR (MT-InSAR) techniques, which allows us to restore the process of ground deformation at millimeter accuracy level by identifying the target points with stable scattering characteristics (PS points) from multiple synthetic aperture radar (SAR) images of the same area during the SAR acquisitions [39,40]. The PS-InSAR technique has obvious advantages in urban land subsidence monitoring, since buildings and artificial facilities provide high-density PS points and greatly improve the signal to noise ratio (SNR) of the interferograms. The purpose of this paper is to provide a better identification and understanding of local subsidence and ground fracture activity in the airport area by using MT-InSAR technique. In this paper, long-term land subsidence at BCIA, Beijng, was obtained by performing an advanced multi-temporal InSAR approach. The dataset consists of 63 TSX/TDX scenes from April 2010 to December 2017. This study has two main purposes. Firstly, an advanced PS-InSAR approach-Quasi-PS InSAR (QPS-InSAR) [41]-was applied to investigate the ground deformation of the BCIA, then extra leveling data was used to evaluate the accuracy of the InSAR-derived measurements. Finally, the characteristic of the newly discovered ground fissures was described and an integrated analysis of hydrogeological data and InSAR-derived measurements was carried out to discuss the main cause of the ground fissures.

Study Site
Beijing, located in the north of China mainland, covers a region of between 39.4 • N 115.7 • E and 41.6 • N 117.4 • E, with a total area of over 16,000 km 2 and a population of 21.542 million (2018). The west, north, and northeast of the city are surrounded by mountains and the southeast is a plain (named Beijing Plain) toward the Bohai sea. The Beijing Plain, with elevation of between 20 and 40 m, is a typical piedmont alluvial diluvial plain formed from sediments carried by the Yongding, Chaobai, and Wenyu rivers. Beijing has a typical warm temperate semi-humid continental monsoon climate, with four distinct seasons; the annual average temperature is approximately 10-12 • C and the mean annual precipitation is about 600 mm (from 1949 to 2018).
Beijing currently has three civil airports, namely the Beijing Capital International Airport, Nanyuan Airport, and the recently completed Daxing International Airport. The BCIA is located in the northeast of Beijing, 25 km from downtown ( Figure 1a). The airport was built in 1958 and has been in operation for more than 60 years. As shown in Figure 1b, which covers an area of 14.8 km 2 , the airport has three terminals (namely T1, T2, and T3, respectively), three runways (01/19, 18L/36R, and 18R/36L; detailed information is shown in Table 1), and 314 aircraft stands. In 2018, the airport handled a total of 614,022 air traffic movements, 100,983,290 passengers, and over 2,074,000 tons of cargo, representing year-on-year decreases of 2.8%, 5.4%, and 2.2%, respectively (statistics were provided by the civil aviation administration of China).
The continuous over-exploitation of groundwater in Beijing Plain has resulted in two main subsidence zones-the northern zone and southern zone. The northern subsidence zone caused by groundwater withdrawal includes several sinking centers (Figure 1a), such as Laiguangying and Balizhuang-dajiaoting in Chaoyang funnel, Nanfaxin in Shunyi funnel, Liyuan-taihu in Tongzhou funnel, and Shahe-baxianzhuang in Changping funnel. The BCIA is located nearby the Nanfaxin sinking center. The airport has a history of land subsidence problems and has been studied in detail [13]. been in operation for more than 60 years. As shown in Figure 1b, which covers an area of 14.8 km 2 , the airport has three terminals (namely T1, T2, and T3, respectively), three runways (01/19, 18L/36R, and 18R/36L; detailed information is shown in Table 1), and 314 aircraft stands. In 2018, the airport handled a total of 614,022 air traffic movements, 100,983,290 passengers, and over 2,074,000 tons of cargo, representing year-on-year decreases of 2.8%, 5.4%, and 2.2%, respectively (statistics were provided by the civil aviation administration of China).   Figure 2 is the distribution map of the quaternary sediments (QS) and major faults, superimposed with compressible thickness (CT) contours and recorded epicenters around the study site. The quaternary sediments with good compressibility near the study area are about 150-250 m thick, which was proved to provide geological conditions for the generation and development of local land subsidence [42]. The Shunyi-liangxiang fault is one of the major faults crossing the study area, with a total length of 35 km and a trend of~30-40 NE. According to the geological survey report of Shunyi county, the Shunyi-liangxiang fault runs through the whole quaternary period stratum and the staggered distance decreases successively from the bottom to the top. The fractured surface of the Shunyi-liangxiang fault extends to the ground and there is an obvious offset at a depth of 5 m, with a maximum vertical dislocation of 14 cm. Also, the early activity of Shunyi-liangxiang faults resulted in an obvious difference in thickness between the upper side (NW) and the lower side (SE), about 300 m, and uneven soil structure on both sides of the fault. , and indicate the quaternary alluvium sediments, the quaternary alluvial-diluvial sediments, the quaternary eolian sediment, the quaternary diluvialalluvium sediments, the quaternary marine sediments, the quaternary lacustrine sediment, and the quaternary diluvial sediment, respectively.

SAR Data Set and Elevation Data
In this study, 63 TerraSAR-X/TanDem-X (TSX/TDX) SAR images in stripmap mode was collected between April 2010 and December 2017 from ascending orbit, with HH polarization. The spatial resolution of the TSX/TDX stripmap mode image is better than 3 m on the ground, which is more suitable for identifying permanent scatterers (PS) in urban areas. The external digital elevation model (DEM) used for MT-InSAR processing to remove the topographic phase and flatten effects was obtained from the TerraSAR-X add-on for Digital Elevation Measurements. The open downloaded TanDEM-X 90m DEM is a product derived from the global DEM with a 0.4 arcsec posting, with an absolute vertical accuracy of 10 m and an relative vertical accuracy of 2 m, which is better than that of SRTM DEM (more information is available at: https://geoservice.dlr.de/web/dataguide/tdm90/). As shown in Figure 3, the TSX/TDX dataset has perpendicular baselines ranging from −178.8 to 297.9 m and temporal baselines ranging from 44 to 1573 days, with respect to the master image dated 2014/08/14. indicate the quaternary alluvium sediments, the quaternary alluvial-diluvial sediments, the quaternary eolian sediment, the quaternary diluvial-alluvium sediments, the quaternary marine sediments, the quaternary lacustrine sediment, and the quaternary diluvial sediment, respectively.

SAR Data Set and Elevation Data
In this study, 63 TerraSAR-X/TanDem-X (TSX/TDX) SAR images in stripmap mode was collected between April 2010 and December 2017 from ascending orbit, with HH polarization. The spatial resolution of the TSX/TDX stripmap mode image is better than 3 m on the ground, which is more suitable for identifying permanent scatterers (PS) in urban areas. The external digital elevation model (DEM) used for MT-InSAR processing to remove the topographic phase and flatten effects was obtained from the TerraSAR-X add-on for Digital Elevation Measurements. The open downloaded TanDEM-X 90m DEM is a product derived from the global DEM with a 0.4 arcsec posting, with an absolute vertical accuracy of 10 m and an relative vertical accuracy of 2 m, which is better than that of SRTM DEM (more information is available at: https://geoservice.dlr.de/web/dataguide/tdm90/). As shown in Figure 3, the TSX/TDX dataset has perpendicular baselines ranging from −178.8 to 297.9 m and temporal baselines ranging from 44 to 1573 days, with respect to the master image dated 2014/08/14.

MT-InSAR Processing
The Quasi-PS technique is an advanced PS-InSAR approach, which is able to extract information of partially coherent targets as the weight of phase estimation [43]. Since the specific information carried by interferograms can be different from point to point, the coherent interferograms subset of each point must be considered [44]. Hence, in QPS-InSAR processing, only the appropriate target-dependent subset of interferograms is used to estimate the target height and displacement. In this study, the quasi-PS processing was carried out using the SARPROZ software (refer to https://www.sarproz.com/). Within the standard single master PS-InSAR processing in SARPROZ, the master scene for coregistration was selected by synthetically taking into account space-time baselines, Doppler centroid frequency different, and the weather information (such as precipitation, temperature, atmospheric pressure, humidity, and cloud cover) at the time of the SAR image acquisition, to minimize the decorrelations of the single master interferograms. After co-registration, the external DEM was imported and a series of interferograms were generated by InSAR processing. Following this, persistent scatterer candidates (PSCs) were chosen by applying a threshold on the amplitude stability index (ASI) for the first time. An atmospheric phase screen (APS) was then estimated and removed the base on a given reference point. Finally, relative height, absolute height, displacement, and deformation rate were collectively estimated from the filtered interferogram subset and TS (time series) deformation was reconstructed simultaneously.
Afterward, only PS points with high temporal coherence (better than 0.75 in this study) were selected for further analysis, to ensure a high coherence and stability within the SAR-image acquisition time-span. Since previous studies pointed out that Beijing Plain presented a low horizontal movement of 1.57-1.93 mm/year [45,46] and vertical displacement is a prerequisite for the comparison with ground-leveling measurements, the displacement from line-of-sight (LOS) was converted to vertical by using Equation (1): where dv is the vertical displacement (mm), dLOS is the line-of-sight (LOS) displacement (mm), and θ is the incidence angle (35.28° for TSX/TDX stripmap mode in this study).

Measurement Calibration
The InSAR technique measures the movement relative to the reference point, so it is crucial to select a stable point as the reference point to estimate the displacement and deformation rate. During QPS processing, points meeting the following requirements can be selected as a reference point: The peak of residual height histogram should be at value 0, which indicates the reference point is on the

MT-InSAR Processing
The Quasi-PS technique is an advanced PS-InSAR approach, which is able to extract information of partially coherent targets as the weight of phase estimation [43]. Since the specific information carried by interferograms can be different from point to point, the coherent interferograms subset of each point must be considered [44]. Hence, in QPS-InSAR processing, only the appropriate target-dependent subset of interferograms is used to estimate the target height and displacement. In this study, the quasi-PS processing was carried out using the SARPROZ software (refer to https://www.sarproz.com/). Within the standard single master PS-InSAR processing in SARPROZ, the master scene for co-registration was selected by synthetically taking into account space-time baselines, Doppler centroid frequency different, and the weather information (such as precipitation, temperature, atmospheric pressure, humidity, and cloud cover) at the time of the SAR image acquisition, to minimize the decorrelations of the single master interferograms. After co-registration, the external DEM was imported and a series of interferograms were generated by InSAR processing. Following this, persistent scatterer candidates (PSCs) were chosen by applying a threshold on the amplitude stability index (ASI) for the first time. An atmospheric phase screen (APS) was then estimated and removed the base on a given reference point. Finally, relative height, absolute height, displacement, and deformation rate were collectively estimated from the filtered interferogram subset and TS (time series) deformation was reconstructed simultaneously.
Afterward, only PS points with high temporal coherence (better than 0.75 in this study) were selected for further analysis, to ensure a high coherence and stability within the SAR-image acquisition time-span. Since previous studies pointed out that Beijing Plain presented a low horizontal movement of 1.57-1.93 mm/year [45,46] and vertical displacement is a prerequisite for the comparison with ground-leveling measurements, the displacement from line-of-sight (LOS) was converted to vertical by using Equation (1): where d v is the vertical displacement (mm), d LOS is the line-of-sight (LOS) displacement (mm), and θ is the incidence angle (35.28 • for TSX/TDX stripmap mode in this study).

Measurement Calibration
The InSAR technique measures the movement relative to the reference point, so it is crucial to select a stable point as the reference point to estimate the displacement and deformation rate. During QPS processing, points meeting the following requirements can be selected as a reference point: The peak of residual height histogram should be at value 0, which indicates the reference point is on the ground and the peak of estimated velocity histogram should also be at value 0, which indicates the reference point is relatively stable.
The reference point of the QPS-InSAR analysis in this study was set near Xuxinzhuang in Tongzhou District (Figure 1a), corresponding to a PS point with a high ASI value (0.98). Since the absolute deformation of the reference point was unknown, the InSAR measurement was calibrated by utilizing the observation records provided by the subsidence monitoring station located in the subsidence center of Pinggezhuang. Unfortunately, the location of the station did not correspond to any PS points, hence the average measurement of the radar targets located within a distance of 100 m from the monitoring station was used instead. Then, the displacement time series of each PS points were updated accordingly using Equation (2). The calibration factor v diff has been quantified as −4.2 mm/year in this study.
where d(t i ) is the corrected displacement at the time, t i d insar is the InSAR-derived displacement at the time, t i v diff is the calibration factor calculated from station observations, and ∆t = t i − t 0 is the time-span (in years) from the first (April 2010) to the i-th scene.

InSAR-derived Measurements
By applying the QPS technique, a total of 389,515 PS points were identified from more than 180 interferograms (as shown in Figure 3), among which about 9000 PS points were located within the BCIA area, with a density of about 640 points/km 2 . Figure 4 illustrates all the PS points detected by the QPS-InSAR approach, colored with the corresponding vertical displacement (cm). The displacements were relative to a reference point in Figure 1a and calibrated by Equation (2). In this study, negative values in all the figures indicate downward movement of the ground, while positive values indicate uplift. Figure 4 shows the overall displacements of the whole coverage from the SAR data frame spanning 2010 to 2017, derived from TSX/TDX dataset by using the QPS technique. The results indicated that significant land subsidence occurred in the eastern Beijing Plain. The unevenness of the spatial distribution was also revealed by the InSAR-derived displacement map. The serious subsidence regions were mainly located in the southwest of Shunyi District, the east of Chaoyang District, and the north of Tongzhou District. In particular, the maximum displacement of −117.4 cm and maximum sinking rate of 152.9 mm/year were both observed in Chaoyang funnel. In addition, it was noticed that multiple subsidence regions, such as the Chaoyang and Tongzhou funnels, joined together with large spatial coverage [47].  In order to show the deformation pattern more clearly, we zoomed in and recolored partial details of the airport area. As shown in Figure 5, most of the PS points were identified on angular targets (e.g., buildings), but many of them were in areas without angular reflection structures, such as runways and platforms in the north of the airport. Predictably, various degrees of land subsidence were observed in the airport area and the figures showed that the land deformed at different rates during April 2010 to December 2017, between −69.6 mm/year (sinking) and +8.1 mm/year (uplift). We drew a boundary, A-A' (marked as the red dashed line in Figure 4), roughly to divide the region into two subregions with distinct subsidence characteristics, according to the distribution of the deformation rate. The northwest showed a clear trend of subsidence (e.g., 18R/36L and 18L/36R runway, terminal T1 and T2, and the Northern Airport Expressway), while the southwest showed little ground motion. However, the deformation rates generated from QPS-InSAR measurements also pointed out that, even at very short distances, there was relative spatial variability on the terminal T2 buildings and north 18L/36R runway. Details will be presented in the next section. In order to show the deformation pattern more clearly, we zoomed in and recolored partial details of the airport area. As shown in Figure 5, most of the PS points were identified on angular targets (e.g., buildings), but many of them were in areas without angular reflection structures, such as runways and platforms in the north of the airport. Predictably, various degrees of land subsidence were observed in the airport area and the figures showed that the land deformed at different rates during April 2010 to December 2017, between −69.6 mm/year (sinking) and +8.1 mm/year (uplift). We drew a boundary, A-A' (marked as the red dashed line in Figure 4), roughly to divide the region into two subregions with distinct subsidence characteristics, according to the distribution of the deformation rate. The northwest showed a clear trend of subsidence (e.g., 18R/36L and 18L/36R runway, terminal T1 and T2, and the Northern Airport Expressway), while the southwest showed little ground motion. However, the deformation rates generated from QPS-InSAR measurements also pointed out that, even at very short distances, there was relative spatial variability on the terminal T2 buildings and north 18L/36R runway. Details will be presented in the next section.

Validation With Ground Leveling Measurements
The InSAR-derived observations in this study were compared with the ground-leveling measurements, including deformation rate and displacement time series, to verify the reliability of the results from QPS-InSAR processing. The averages of the measurements from pixels lying nearby (with a distance of less than 100 m) each benchmark were calculated as the corresponding InSAR measurements by assuming that the ground did not change significantly within this distance. As shown in Figure 6, a cross-comparison was performed between InSAR and leveling measurements by linear regression analysis. The deformation rates from leveling and InSAR measurements show good agreement, with a correlation coefficient of 0.96 and a standard error of 1.33 mm/year.
To further validate our results, the InSAR-derived displacement time series were compared with ground-leveling measurements at the selected two benchmarks (marked by black triangles in Figure  3). Figure 7 compares the time series plots of measurements from the two leveling sites and InSAR. It can be observed that the displacement time series of the two methods are in good agreement. The root mean square (RMS) errors for InSAR measurements at benchmark #1 and #2 were 6.5 mm and 3.8 mm, respectively. These comparisons suggest that the MT-InSAR technique is effective and reliable in extracting ground deformation information. In addition, since the observing frequency of leveling (in years) is much lower than that of the spaceborne SAR system (almost monthly), it is not possible to describe the deforming process in detail over the whole period.

Validation With Ground Leveling Measurements
The InSAR-derived observations in this study were compared with the ground-leveling measurements, including deformation rate and displacement time series, to verify the reliability of the results from QPS-InSAR processing. The averages of the measurements from pixels lying nearby (with a distance of less than 100 m) each benchmark were calculated as the corresponding InSAR measurements by assuming that the ground did not change significantly within this distance. As shown in Figure 6, a cross-comparison was performed between InSAR and leveling measurements by linear regression analysis. The deformation rates from leveling and InSAR measurements show good agreement, with a correlation coefficient of 0.96 and a standard error of 1.33 mm/year.
To further validate our results, the InSAR-derived displacement time series were compared with ground-leveling measurements at the selected two benchmarks (marked by black triangles in Figure 3). Figure 7 compares the time series plots of measurements from the two leveling sites and InSAR. It can be observed that the displacement time series of the two methods are in good agreement. The root mean square (RMS) errors for InSAR measurements at benchmark #1 and #2 were 6.5 mm and 3.8 mm, respectively. These comparisons suggest that the MT-InSAR technique is effective and reliable in extracting ground deformation information. In addition, since the observing frequency of leveling (in years) is much lower than that of the spaceborne SAR system (almost monthly), it is not possible to describe the deforming process in detail over the whole period.

Ground Fissures Identified By InSAR
As described above, the land subsidence in the capital airport area was obviously zoned and the deformation rates in the northwest and southeast regions are quite different. The uneven land subsidence in the airport area has lasted for several years and there will definitely be accumulation effect, e.g., damage to runways, buildings, and electric equipment. In order to ascertain whether the uneven ground deformation has affected the normal activities of the airport, we collected some relevant materials, including news articles, media reports, papers, and internal reports that cannot be

Ground Fissures Identified By InSAR
As described above, the land subsidence in the capital airport area was obviously zoned and the deformation rates in the northwest and southeast regions are quite different. The uneven land subsidence in the airport area has lasted for several years and there will definitely be accumulation effect, e.g., damage to runways, buildings, and electric equipment. In order to ascertain whether the uneven ground deformation has affected the normal activities of the airport, we collected some relevant materials, including news articles, media reports, papers, and internal reports that cannot be

Ground Fissures Identified By InSAR
As described above, the land subsidence in the capital airport area was obviously zoned and the deformation rates in the northwest and southeast regions are quite different. The uneven land subsidence in the airport area has lasted for several years and there will definitely be accumulation Remote Sens. 2019, 11, 1466 11 of 17 effect, e.g., damage to runways, buildings, and electric equipment. In order to ascertain whether the uneven ground deformation has affected the normal activities of the airport, we collected some relevant materials, including news articles, media reports, papers, and internal reports that cannot be made public. As reported, construction of the 18L/36R runway started in 1954 and was completed in 1958 and it has been used for more than 20 years since 1996 when asphalt concrete was added. In recent years, the differential land subsidence caused relative ground fluctuations within a short distance, which affects the safe operation of the aircraft when passing through. The runway was repaired partially in 2013 and 2015, but was completely overhauled in 2017. The 28-day complete overhaul in April 2017, involving an area of 200,000 m 2 , was partly aimed at repairing runway cracks and fluctuations caused by uneven subsidence. In the meantime, as field investigations show, the terminal T2 that was put into use at the end of 1999 also suffered from land subsidence and local ground fissures.
In order to locate the ground fissures from the InSAR-derived measurements, we calculated the gradient of land subsidence rate in the study area. Thus, the general location of the fracture was determined by profile analysis along the direction in which the InSAR-derived deformation rate changed the most. The analysis indicated that the ground fissures were distributed roughly along the~35-55 • N direction. The primary direction of the ground fissures was identified as 43.1 • N and is marked as the red dashed line in Figure 8.
Remote Sens. 2019, 11, x FOR PEER REVIEW 11 of 17 made public. As reported, construction of the 18L/36R runway started in 1954 and was completed in 1958 and it has been used for more than 20 years since 1996 when asphalt concrete was added. In recent years, the differential land subsidence caused relative ground fluctuations within a short distance, which affects the safe operation of the aircraft when passing through. The runway was repaired partially in 2013 and 2015, but was completely overhauled in 2017. The 28-day complete overhaul in April 2017, involving an area of 200,000 m 2 , was partly aimed at repairing runway cracks and fluctuations caused by uneven subsidence. In the meantime, as field investigations show, the terminal T2 that was put into use at the end of 1999 also suffered from land subsidence and local ground fissures. In order to locate the ground fissures from the InSAR-derived measurements, we calculated the gradient of land subsidence rate in the study area. Thus, the general location of the fracture was determined by profile analysis along the direction in which the InSAR-derived deformation rate changed the most. The analysis indicated that the ground fissures were distributed roughly along the ~35-55°N direction. The primary direction of the ground fissures was identified as 43.1°N and is marked as the red dashed line in Figure 8.
In this case, we took the segment (about 4 km) where the ground fissure was located as the centerline (zero-scribed line in Figure 9) and extended the buffer zone of 6 km on both sides. Then, an area of about 12 × 4 km 2 was obtained (marked as the black box in Figure 8). We plotted the cumulative vertical displacement of the points within the black box onto a profile sketch and colored each point according to the deformation rate.  As shown in Figure 9, the minimum value of the vertical displacement from the points on the SE (lower) side of the ground fissure is −141.3 mm, with an average settlement of 250.9 mm and a sinking rate of 33.2 mm/year; the maximum value of the vertical displacement from the points on the NW (upper) side of the ground fissure is −680.7 mm, with an average settlement of 473.5 mm and a sinking rate of 61.7 mm/year. The maximum relative displacement on both sides of the ground fissure is over 530 mm, while the difference in average settlement and sinking rate on both sides reach 222.6 mm and 28.5 mm/year, respectively. The distribution characteristics of ground subsidence can more or less provide clues to the activity of the new ground fissures, especially after 2010, since runway 18L/36R has dropped more than 200 mm in the north relative to the sinking surface from 2003 to 2013, as shown in our previous study [13]. The differences between the sinking rate of points on both sides of the ground fissure have In this case, we took the segment (about 4 km) where the ground fissure was located as the centerline (zero-scribed line in Figure 9) and extended the buffer zone of 6 km on both sides. Then, an area of about 12 × 4 km 2 was obtained (marked as the black box in Figure 8). We plotted the cumulative vertical displacement of the points within the black box onto a profile sketch and colored each point according to the deformation rate.
As shown in Figure 9, the minimum value of the vertical displacement from the points on the SE (lower) side of the ground fissure is −141.3 mm, with an average settlement of 250.9 mm and a sinking rate of 33.2 mm/year; the maximum value of the vertical displacement from the points on the NW (upper) side of the ground fissure is −680.7 mm, with an average settlement of 473.5 mm and a sinking rate of 61.7 mm/year. The maximum relative displacement on both sides of the ground fissure is over 530 mm, while the difference in average settlement and sinking rate on both sides reach 222.6 mm and 28.5 mm/year, respectively.
The distribution characteristics of ground subsidence can more or less provide clues to the activity of the new ground fissures, especially after 2010, since runway 18L/36R has dropped more than 200 mm in the north relative to the sinking surface from 2003 to 2013, as shown in our previous study [13]. The differences between the sinking rate of points on both sides of the ground fissure have further increased. As illustrated in Figure 10, the sinking rate on the roof from the northern and southern waiting halls of terminal T2 were 37.6 mm/year (P1) and 7.1 mm/year (P2), respectively. Even more than that, the InSAR-derived measurements showed the relative deformation rate of runway 18L/36R exceeded 47 mm/year, with a relative displacement of 365.4 mm during April 2010 to December 2017. The time series displacements and deformation rates of the selected P1, P2, P3, and P4 are listed on the right panel in Figure 10. Meanwhile, the activity of the newly identified ground fissure in the capital airport area is increasing and it requires constant attention.
As shown in Figure 9, the minimum value of the vertical displacement from the points on the SE (lower) side of the ground fissure is −141.3 mm, with an average settlement of 250.9 mm and a sinking rate of 33.2 mm/year; the maximum value of the vertical displacement from the points on the NW (upper) side of the ground fissure is −680.7 mm, with an average settlement of 473.5 mm and a sinking rate of 61.7 mm/year. The maximum relative displacement on both sides of the ground fissure is over 530 mm, while the difference in average settlement and sinking rate on both sides reach 222.6 mm and 28.5 mm/year, respectively. The distribution characteristics of ground subsidence can more or less provide clues to the activity of the new ground fissures, especially after 2010, since runway 18L/36R has dropped more than 200 mm in the north relative to the sinking surface from 2003 to 2013, as shown in our previous study [13]. The differences between the sinking rate of points on both sides of the ground fissure have

Discussion
As one of the busiest airports in the world, it is difficult to carry out underground detection or other related field work for the BCIA area. In this case, the cause of ground fissures can only be preliminarily discussed through limited hydrogeological data. In order to preliminarily analyze the genesis of ground fissures, we collected limited geological data and plotted them on the interpolated deformation rate map ( Figure 11). As can be seen, the ground fissure is located in the southern section of the Shunyi-liangxiang faults and the strike of ground fissure is basically the similar as the active fault.
Remote Sens. 2019, 11, x FOR PEER REVIEW 13 of 17 further increased. As illustrated in Figure 10, the sinking rate on the roof from the northern and southern waiting halls of terminal T2 were 37.6 mm/year (P1) and 7.1 mm/year (P2), respectively. Even more than that, the InSAR-derived measurements showed the relative deformation rate of runway 18L/36R exceeded 47 mm/year, with a relative displacement of 365.4 mm during April 2010 to December 2017. The time series displacements and deformation rates of the selected P1, P2, P3, and P4 are listed on the right panel in Figure 10. Meanwhile, the activity of the newly identified ground fissure in the capital airport area is increasing and it requires constant attention.

Discussion
As one of the busiest airports in the world, it is difficult to carry out underground detection or other related field work for the BCIA area. In this case, the cause of ground fissures can only be preliminarily discussed through limited hydrogeological data. In order to preliminarily analyze the genesis of ground fissures, we collected limited geological data and plotted them on the interpolated deformation rate map ( Figure 11). As can be seen, the ground fissure is located in the southern section of the Shunyi-liangxiang faults and the strike of ground fissure is basically the similar as the active fault. Figure 11. The distribution of the quaternary sedimentary types, groundwater level (GWL), and active faults, with interpolated deformation rate (mm/year) in the study area. The black outlines show airport features, including terminals and airport boundaries. A-A' indicates the location of the ground fissures identified by InSAR measurements, while the B-B' indicates the same from Ma and Du, 2017 [48].
The borehole investigation showed that the vertical slip rate of the fracture was about 0.03 mm/year [49], far less than the deformation rate detected by InSAR, which indicates that the activity of Shunyi-liangxiang faults cannot fully explain recent ground fissure activity. The borehole investigation showed that the vertical slip rate of the fracture was about 0.03 mm/year [49], far less than the deformation rate detected by InSAR, which indicates that the activity of Shunyi-liangxiang faults cannot fully explain recent ground fissure activity.
Generally, the ground fissures formed by uneven deformation caused by groundwater pumping will present the distribution characteristics related to the height of groundwater level (GWL); where GWLs are low, subsidence tends to be more severe. However, the distribution of GWL and land subsidence in the airport area did not show this correlation (as shown in Figure 11). It is preliminarily proved that the variation of GWL is not the main reason for the occurrence and development of the ground fissures.
In Shunyi County, the long-term overexploitation of the confined aquifer of~70-200 m underground resulted in the decline of the confined GWL, the decrease of pore water pressure and the increase of effective stress in the aquifer. Since the groundwater cannot be replenished in time, the skeleton of the aquifer was compressed and realized as land subsidence of the ground surface. The activity of ground fissures is related to the ground subsidence centers in Shunyi County. As mentioned in the second part of this study, the airport is located at the southern edge of the newly formed Nanfaxin sinking center and is adjacent to the Houshayu sinking center in the west, which is obviously affected by local uneven subsidence. The ground fissure is located at the southwest edge of the Houshayu-nanfaxin subsidence area and the relative average of subsidence on its both sides reaches 22.26 cm, which is in line with the field investigation (a maximum displacement of 30 cm was obtained) implemented by CGS in late 2017; the period of ground fissure activity (identified by InSAR) coincides with the rapid development of land subsidence (since 1999). The normal faulting of the Shunyi-liangxiang fault in the early and middle quaternary period resulted in the development of thick, fine-grained, loose deposits with high compressibility in the lower side of the fault. In this case, the differential subsidence of the northern and southern sides of the fault resulted in the formation of the faulted zone in the airport area.
The overhaul of the 18L/36R runway lasted 28 days in April 2017, during which the number of daily flights decreased by more than 300 and the number of passengers sent was cut by about 15 percent. Furthermore, the airport had to change from a three-runway operation to a two-runway operation, which directly led to a decrease in runway capacity during rush hour, an extension of night flight time, and an increase in flight delay rate. It can be concluded that uneven ground subsidence and fractures seriously affected the normal operation of the airport.
The quaternary sediments at the location of the new ground fissure are about 400 m thick and the fractures penetrate from the bedrock to the ground surface; engineering measures barely eliminate the existence and development of ground fissures [48]. Uneven land subsidence and tectonic activities of Shunyi-liangxiang faults continuously sustain and enhance the ground fissures, which may induce geological hazards. Therefore, in order to ensure the safe operation of the airport, InSAR technology can be used for long-term and near-real-time monitoring of the airport area in cooperation with ground inspection. In this way, early warnings and engineering measures can be taken in time when abnormal conditions are found.

Conclusions
In this study, we performed an InSAR time series analysis for investigating the long-term uneven land subsidence and intensifying ground fissures since 2010 in Beijing Capital International Airport, China. An advanced InSAR technique, referred to as QPS-InSAR, was performed on TerraSAR-X/TanDem-X data set, then spatial and time series measurements were obtained during the period of investigation from 2010 to 2017. The vertical displacements have been estimated to have a maximum value of −1174.1 mm, with a sinking rate of 152.9 mm/year. Comparison of the deformation rates from the InSAR and ground-leveling measurements shows the correlation coefficient is 0.96 with a standard deviation of 1.33 mm/year. After an integrated analysis of the distribution characteristics of land subsidence in Beijing Capital International Airport, previous research results, and geological data, we found and located an active ground fissure. Then, the main cause of the ground fissures was preliminarily discussed. We believe that the uneven land subsidence relating to the difference in sediment thickness on both sides of the Shunyi-liangxiang fault is the main factor of the formation of the new ground fissures, rather than the normal faulting itself and the variation of GWL. Due to the special environment of the airport where the ground fissures are located, we suggest the management of the airport combine spaceborne InSAR technology and ground surveys to monitor and forecast ground fissures, in order to provide technical guarantee for the safe operation of the airport. It can be concluded that InSAR technology can be used to identify and monitor geological processes such as land subsidence and ground fissure activities and can provide a scientific approach to better explore and study the cause and formation mechanisms of regional subsidence and ground fissures. Furthermore, it is worth remarking that this refreshes our understanding of the complex cause mechanisms of local land subsidence in the Beijing Plain area.