Characterizing the Topographic Changes and Land Subsidence Associated with the Mountain Excavation and City Construction on the Chinese Loess Plateau

: A mega project, Mountain Excavation and City Construction (MECC), was launched in the hilly and gully region of the Chinese Loess Plateau in 2012, in order to address the shortage of available land and create new ﬂat land for urban construction. However, large-scale land creation and urban expansion signiﬁcantly alters the local geological environment, leading to severe ground deformation. This study investigated the topographic changes, ground deformation, and their interactions due to the MECC project in the Yan’an New District (YND). First, new surface elevations were generated using ZiYuan-3 (ZY-3) stereo images acquired after the construction in order to map the local topographic changes and the ﬁll thickness associated with the MECC project. Then, the interferometric synthetic aperture radar (InSAR) time series and 32 Sentinel-1A images were used to assess the spatial patterns of the ground deformation in the YND during the postconstruction period (2017–2018). By combining the InSAR-derived results and topographic change features, the relationship between the ground deformation and large-scale land creation was further analyzed. The results indicated that the MECC project in the YND has created over 22 km 2 of ﬂat land, including 10.8 km 2 of ﬁlled area, with a maximum ﬁll thickness of ~110 m. Signiﬁcant uneven ground deformation was detected in the land-creation area, with a maximum subsidence rate of approximately 121 mm/year, which was consistent with the ﬁeld survey. The strong correlation between the observed subsidence patterns and the land creation project suggested that this recorded uneven subsidence was primarily related to the spatial distribution of the ﬁlling works, along with the changes in the thickness and geotechnical properties of the ﬁlled loess; moreover, rapid urbanization, such as road construction, can accelerate the subsidence process. These ﬁndings can guide improvements in urban planning and the mitigation of geohazards in regions experiencing large-scale land construction.


Introduction
The Loess Plateau, covering an area of approximately 440,000 km 2 in north-central China, is the most widely distributed, thickest loess area in the world [1]. With the rapid promotion of the 'western development strategy' and the 'Belt and Road' policy in China, several mega engineering projects have been underway on the Loess Plateau, such as the 'Mountain Excavation and City Construction' (MECC) [2], the 'Gully Control and Highland Protection' (GCHP) [3] and the 'Gully Land Consolidation' (GLC) [4]. The pace and intensity of these projects far exceed the geological bearing capacity and therefore [25,26], Hong Kong [12], Beijing [27,28], and Xi'an [29].
The main objectives of this study are (1) to characterize the spatiotemporal patterns of the topographic changes and land subsidence associated with the MECC project in the YND, and (2) to quantify the relationships between land subsidence and the MECC project in the YND. In this paper, we illustrate an integrated approach which uses space-borne optical photogrammetry, InSAR techniques and field surveys to investigate the topographic changes and ground deformation due to the large-scale MECC project and their interactions. First, space-borne optical photogrammetry was used to extract new surface elevations, and to map the local topographic and geomorphological changes caused by the MECC project. Moreover, the magnitude and spatial distribution of the ground deformation in the YND was derived from 32 Sentinel-1A images using the SBAS-InSAR method. Then, the relationships between the observed deformations and the large-scale MECC project were discussed in detail by combining the InSAR-derived results with the topographic change features and optical images. Finally, other possible natural and anthropogenic factors influencing ground settlement were analyzed and discussed. These results will facilitate a better understanding of the surface deformation induced by the MECC project, which can aid in urban planning and disaster prevention and mitigation in loess regions where large-scale land creation and rapid city construction projects occur.

Study Area
Yan'an city is located in the hilly-gully area of the south-central part of the Loess Plateau, China (Figure 2b). The old city is a typical linear city, with densely-distributed urban buildings along a Y-shaped river valley (Figure 2c). The narrow and long urban area of approximately 36 km 2 accommodates nearly 500,000 people, making the local construction land per capita far below the national standard for land utilization. Simultaneously, Yan'an city, as a revolutionary holy land of China, has over 140 historical relics [5]. However, in order to obtain more land for urban development, a large number of revolutionary sites have been occupied and destroyed, which intensified the competition between urban and landscape areas. With rapid urban development and population growth, the shortage of land resources for construction in the old city has further intensified, which has severely restricted the economic and social development of Yan'an city. In order to Yan'an, Shaanxi Province, is a city with a severe shortage of available land, with approximately 500,000 people living within a Y-shaped urban area of only 36 km 2 . In order to solve the problem of overcrowding in this historic city, the MECC project was started in April 2012 to create a new district of 78 km 2 , known as the Yan'an New District [6]. The MECC project in the YND is regarded as the largest geotechnical project in a collapsible loess region worldwide [2]. This unprecedented large-scale land creation and rapid urbanization has drawn great attention from scientists and society. Although the large-scale construction of the MECC project can create more available space for urban expansion, it could also significantly alter the original topography and create large loess fill areas, which may cause severe land subsidence and pose new challenges for urban development [5,6]. Several studies have investigated the mechanism and spatiotemporal patterns of land subsidence in the YND using geodetic methods, such as leveling [7], global positioning systems (GPSs), and interferometric synthetic aperture radar (InSAR) [5]. The results of these previous studies have highlighted the driving effects of human activities associated with large-scale land creation on the development of ground subsidence in the YND. For example, the in situ monitoring of filled loess embankments has indicated that the subsidence of the fill body itself accounts for the majority of the total subsidence, and the filling work related to land creation is the primary cause of land subsidence in the YND [7]. Wu et al. (2019) assessed the spatiotemporal patterns of the subsidence in the YND based on the small baseline subset InSAR (SBAS-InSAR) technique, and confirmed the ability of the SBAS-InSAR method to detect surface deformation in the land creation region; the results indicated that the land subsidence is mainly caused by earth filling and building loads [5]. However, the spatial characteristics and interrelations of the local topographic changes and land subsidence associated with land creation are not yet fully understood, largely due to the absence of high accuracy postconstruction topographic Remote Sens. 2021, 13,1556 3 of 20 data that could reflect the surface processes caused by land creation and provide detailed information on the fill thickness.
Remote sensing technology has proven to be an effective tool for investigating land use and topographic changes, and for monitoring land subsidence [8][9][10]. Space-borne optical photogrammetry provides the capability to extract new surface elevation, and is widely used to estimate local topographic changes due to glacial activity, geological hazards, or human activities [11]. In recent years, with the rapid development of the InSAR technique, it has become feasible and reliable to map surface displacement at very high spatiotemporal resolutions in land-creation areas. Compared with ground-based geodetic methods, such as spirit leveling and GPS, the InSAR technique has the unique capability of monitoring large-scale surface deformation with high accuracy and high spatial resolution under all weather conditions [12,13], and it has been widely applied to the detection and monitoring of various geohazards, including earthquakes [14], landslides [15][16][17], land subsidence [18], and other geologic processes [8,19]. In particular, advanced multitemporal InSAR (MT-InSAR) methods, such as persistent scatterers InSAR (PS-InSAR) [20,21] and small baseline subset InSAR (SBAS-InSAR) [22], have been rapidly developed to reduce the effects of decorrelation and uneven atmospheric delay, and to remove the topographic phase, which are intrinsic limitations of the traditional differential InSAR (D-InSAR) method. Accordingly, MT-InSAR methods have facilitated the monitoring of ground deformation in many countries and cities, such as Mexico City [23,24], Italy [25,26], Hong Kong [12], Beijing [27,28], and Xi'an [29].
The main objectives of this study are (1) to characterize the spatiotemporal patterns of the topographic changes and land subsidence associated with the MECC project in the YND, and (2) to quantify the relationships between land subsidence and the MECC project in the YND. In this paper, we illustrate an integrated approach which uses space-borne optical photogrammetry, InSAR techniques and field surveys to investigate the topographic changes and ground deformation due to the large-scale MECC project and their interactions. First, space-borne optical photogrammetry was used to extract new surface elevations, and to map the local topographic and geomorphological changes caused by the MECC project. Moreover, the magnitude and spatial distribution of the ground deformation in the YND was derived from 32 Sentinel-1A images using the SBAS-InSAR method. Then, the relationships between the observed deformations and the large-scale MECC project were discussed in detail by combining the InSAR-derived results with the topographic change features and optical images. Finally, other possible natural and anthropogenic factors influencing ground settlement were analyzed and discussed. These results will facilitate a better understanding of the surface deformation induced by the MECC project, which can aid in urban planning and disaster prevention and mitigation in loess regions where large-scale land creation and rapid city construction projects occur.

Study Area
Yan'an city is located in the hilly-gully area of the south-central part of the Loess Plateau, China (Figure 2b). The old city is a typical linear city, with densely-distributed urban buildings along a Y-shaped river valley (Figure 2c). The narrow and long urban area of approximately 36 km 2 accommodates nearly 500,000 people, making the local construction land per capita far below the national standard for land utilization. Simultaneously, Yan'an city, as a revolutionary holy land of China, has over 140 historical relics [5]. However, in order to obtain more land for urban development, a large number of revolutionary sites have been occupied and destroyed, which intensified the competition between urban and landscape areas. With rapid urban development and population growth, the shortage of land resources for construction in the old city has further intensified, which has severely restricted the economic and social development of Yan'an city. In order to realize the sustainable development of urbanization, and to protect the historical and cultural heritage of revolutionary sacred land, the MECC project was implemented to create flat land and build the YND by removing the tops of loess mountains to fill in the valleys [2] (Figure 1). is now characterized by relatively flat topography (Figure 2d), with an average elevation of 1100 m. The study area is mainly covered by Middle Pleistocene Lishi loess and Late Pleistocene Malan loess with strong collapsibility, which have been the primary sources of material for the MECC project. The climate is characteristic of a semi-humid and semiarid continental monsoon climate, with an average annual temperature of 10.3 °C [4] and an average annual precipitation of approximately 562 mm. Nearly 70% of the annual precipitation is concentrated in the period from June to September [5].

ZY-3 Satellite Images
ZiYuan-3 (ZY-3) is the first Chinese civilian multispectral satellite with high-resolution optical stereoscopic mapping capability; it was launched on 9 January 2012 [30]. This satellite can compose three-line-array (TLA) stereo images from three panchromatic cameras at different viewing angles (forward, nadir and backward). A scene of ZY-3 TLA data  Figure 2d). The original topography of the YND was a typical loess hilly and gully landform with many gullies. The terrain elevation ranged from 950 to 1260 m, and the valley depths ranged from 100 to 150 m. Following the completion of the land creation project, the YND is now characterized by relatively flat topography (Figure 2d), with an average elevation of 1100 m. The study area is mainly covered by Middle Pleistocene Lishi loess and Late Pleistocene Malan loess with strong collapsibility, which have been the primary sources of material for the MECC project. The climate is characteristic of a semi-humid and semiarid continental monsoon climate, with an average annual temperature of 10.3 • C [4] and an average annual precipitation of approximately 562 mm. Nearly 70% of the annual precipitation is concentrated in the period from June to September [5].

ZY-3 Satellite Images
ZiYuan-3 (ZY-3) is the first Chinese civilian multispectral satellite with high-resolution optical stereoscopic mapping capability; it was launched on 9 January 2012 [30]. This satellite can compose three-line-array (TLA) stereo images from three panchromatic cameras at different viewing angles (forward, nadir and backward). A scene of ZY-3 TLA data consists of images viewed from the nadir, backward, and forward angles with spatial resolutions of Remote Sens. 2021, 13, 1556 5 of 20 2.5 m, 3.5 m, and 3.5 m, respectively. The TLA stereo images of the ZY-3 satellite permit the generation of digital elevation models (DEMs) for the estimation of the terrain elevation and its changes [11]. The stereo-imagery pair composed of forward and backward view images are taken in almost identical lighting and weather conditions. Therefore, they are more applicable for DEM generation using automated DEM extraction techniques than the stereo pair taken by the crosstrack orbits [31]. Moreover, rational polynomial coefficients (RPCs) are provided together with the ZY-3 TLA data, providing information on the mathematical transformation between the geographic coordinates and the image coordinates [32]. In this study, three image data sets with cloudless and clear features on 17 November 2013, 12 May 2015, and 13 October 2017 were collected from the Satellite Surveying and Mapping Application Center (SASMAC) and clipped to extract the DEM, in order to estimate the topographic change and fill thickness in the study area. Detailed information on ZY-3 is presented in Table 1. The large-scale land creation projects and urban construction in the YND commenced in April 2012 and were still underway until August 2018, which induced great changes in the land surface. Therefore, in order to reduce the effect of the surface changes caused by land creation on the decorrelation of SAR images, and to acquire comprehensive and detailed surface deformation information on the entire YND, we collected 32 archived SAR images with a spatial resolution of 5 m × 20 m (range × azimuth) acquired from December 2018 to December 2019 by C-band Sentinel-1A satellites, with a sampling interval of 12 days. Moreover, precise orbit ephemerides (POEs) data from the European Space Agency were applied to the initial estimates of the interferometric baselines. In order to remove the topography-related phase from the interferometric phase and geocode the InSAR products, a shuttle radar topography mission (SRTM) with a spatial resolution of 30 m was applied as external DEM data in this study. Moreover, in order to explore the spatial association between the surface deformation and anthropogenic factors, multitemporal Google Earth images with a resolution of 0.48 m were utilized to aid in the interpretation and analysis of the results in this study. Detailed information about the SAR data used in this study is shown in Table 1, and their coverages are shown in Figure 2b.

Topographic Change Calculation
DEM extraction using stereo image pairs has developed into a relatively mature technology, which started from analytical stereo plotters in the 1980s [31]. With the development of digital image processing techniques, many commercial and open-source software programs have provided modules for DEM generation. In this study, the backward-forward image pairs provided by the ZY-3 TLA data were processed using the rational functional model (RFM) in TitanImage software in order to perform the extraction of the DEM (hereinafter ZY-3 DEM). The RFM is the major imaging geometric model currently employed for the block adjustment of most satellite images [33]. The RFM establishes a mathematical transformation between image coordinates and their corresponding object coordinates using RPCs [32]. However, DEM generation using RFM based solely on RPCs is biased, and needs ground control points (GCPs) to refine it, but we relied only on intensive and well-distributed tie points to enhance the solution because GCPs were unavailable [34]. Previous studies have documented that the DEM generated by ZY-3 stereo images can achieve high vertical accuracy with or without GCPs [30,31]. Hence, ZY-3 image data can be used to extract the DEM. In total, 61 evenly-distributed tie points were selected, and the ZY-3 DEM was generated with a spatial resolution of 5 m and the same projection information as the SRTM DEM.
Moreover, after we acquired the postconstruction ZY-3 DEM, the topographic changes in 2000-2013, 2000-2015 and 2000-2017 were estimated by calculating the difference between the SRTM and the ZY-3 DEMs ( Figure 3). The DEMs were matched and resampled to the same resolution prior to the calculation, as relative geolocation errors are common between DEMs [8,11]. The elevation change results in 2000-2017 were employed in order to discuss the uncertainty excluding the YND affected by the MECC project. The mean and standard deviation of the elevation difference in non-YND areas were used to characterize the relative elevation accuracy of the ZY-3 DEM [11].

Topographic Change Calculation
DEM extraction using stereo image pairs has developed into a relatively mature tec nology, which started from analytical stereo plotters in the 1980s [31]. With the develo ment of digital image processing techniques, many commercial and open-source softwa programs have provided modules for DEM generation. In this study, the backward-fo ward image pairs provided by the ZY-3 TLA data were processed using the rational fun tional model (RFM) in TitanImage software in order to perform the extraction of the DEM (hereinafter ZY-3 DEM). The RFM is the major imaging geometric model currently em ployed for the block adjustment of most satellite images [33]. The RFM establishes a mat ematical transformation between image coordinates and their corresponding object coo dinates using RPCs [32]. However, DEM generation using RFM based solely on RPCs biased, and needs ground control points (GCPs) to refine it, but we relied only on inte sive and well-distributed tie points to enhance the solution because GCPs were unavail ble [34]. Previous studies have documented that the DEM generated by ZY-3 stereo im ages can achieve high vertical accuracy with or without GCPs [30,31]. Hence, ZY-3 imag data can be used to extract the DEM. In total, 61 evenly-distributed tie points were s lected, and the ZY-3 DEM was generated with a spatial resolution of 5 m and the sam projection information as the SRTM DEM.
Moreover, after we acquired the postconstruction ZY-3 DEM, the topograph changes in 2000-2013, 2000-2015 and 2000-2017 were estimated by calculating the diffe ence between the SRTM and the ZY-3 DEMs ( Figure 3). The DEMs were matched an resampled to the same resolution prior to the calculation, as relative geolocation errors a common between DEMs [8,11]. The elevation change results in 2000-2017 were employe in order to discuss the uncertainty excluding the YND affected by the MECC project. Th mean and standard deviation of the elevation difference in non-YND areas were used characterize the relative elevation accuracy of the ZY-3 DEM [11].

InSAR-Derived Deformation Rates
The conventional differential InSAR (D-InSAR) technique uses two SAR images obtained at different times and covers a common region to retrieve surface displacement from the phase differences detected during consecutive SAR acquisitions [35]. However, the precision of D-InSAR is always affected by several factors, such as spatial and temporal decorrelation, atmospheric artifacts [36], and an incorrect DEM. In order to mitigate the effects of these issues, advanced multitemporal InSAR techniques-such as SBAS-InSAR and PS-InSAR [20,21]-have been developed to detect ground displacements.
The SBAS-InSAR technique is an advanced InSAR time series analysis method that relies on multimaster images to calculate ground deformation using interferometric pairs with shorter spatio-temporal baselines [22]. This technique effectively reduces the effects of spatial and/or temporal decorrelation and DEM errors on parameter estimation, which are common in the traditional D-InSAR processing technique [37]. SBAS-InSAR technology has been widely used for the monitoring of large-area surface deformation associated with anthropogenic activities, such as urbanization and groundwater extraction [5,9,38].
In order to acquire more reliable and higher-resolution time series deformation results in the land-creation area with a relatively high surface deformation rate, 32 scene Sentinel-1A images were processed using the SBAS-InSAR technique. First, all of the SAR images were accurately coregistered, and high-quality interferograms were generated based on small temporal and spatial baseline thresholds (less than 36 days and smaller than 100 m). External DEM data were utilized to simulate and remove the topographic phases. Then, multilooking and Goldstein filtering methods were applied in order to improve the signal-to-noise ratio and coherence of the interferograms. Phase unwrapping with the robust minimum cost flow (MCF) algorithm [39] was conducted on the coherent pixels for each wrapped interferogram.
The unwrapped interferometric phase at a pixel in the jth differential interferogram generated by two SAR images acquired at the start time t A and end time t B can be described as follows: where ϕ(t B ) and ϕ(t A ) are the phase signals at SAR acquisition dates t A and t B (t A < t B ), respectively; δϕ de f j represents the phase contributions from surface deformation along the line-of-sight (LOS) direction; δϕ topo j is the inaccurate DEM residual phase; δϕ atm j is the atmospheric phase caused by the difference in atmospheric delay between two observations; and δϕ noise j is the noise phase. In order to improve the accuracy and reliability of the deformation results, these nuisance terms should be effectively estimated and eliminated from differential interferograms. The residual topographic phase can be estimated by utilizing its linear relationship with perpendicular baselines and the least squares method [40]. The atmospheric phase was estimated and subtracted from the interferometric phase by using a low-pass spatial filter and a high-pass temporal filter. Finally, the deformation velocity and time series deformation results were retrieved using the singular value decomposition (SVD) method after removing the abovementioned residual phase components from the interferometric phase. The detailed flowchart of the time series InSAR processing method is depicted in Figure 3 Table 2 summarizes the fill thicknesses of the two main valley-filling areas. Generally, the spatial evolution of land creation and the corresponding topographic changes were clearly recorded, and three main stages could be roughly reflected. The first phase of land creation was completed before 17 November 2013, and was mainly concentrated in the Qiaoergou area, covering an area of approximately 10.6 km 2 . The thickness of the filling loess in the Qiaoergou area ranges from 20 m to 60 m, with an average of 36 m and a maximum of approximately 110 m, which is generally consistent with previous reports [7]. The second phase aimed mainly to create land in the small watershed in the western part of the YND, with an additional land area of 8.     The accuracy of the ZY-3 DEM was expressed as the mean and standard deviation of the elevation difference in the non-YND areas between the SRTM and ZY-3 DEMs from 2000 to 2017. Figure 5 delineates the difference map and the corresponding histogram of the elevation difference for the areas without land-creation work. The elevation differences are approximately normally distributed, and 65% of the rasters have absolute elevation differences less than 10 m. The mean value of the elevation difference is 3.7 m, while the standard deviation is 12.8 m. Compared with the excavation and filling of more than one hundred meters, this error has a limited impact on the estimation of the elevation changes in the YND. Basically, these statistical results show good agreement between the elevations of the two DEMs obtained from different sources. Some large local differences may be caused by various factors, such as the error contained in the preconstruction DEM, urban construction, gully land consolidation, and local land creation.

Spatial Pattern
In this study, SBAS-InSAR technology was employed to quantify the spatiotemporal pattern of the surface deformation in the YND from December 2017 to December 2018. According to previous studies, the vertical deformation in the YND is much larger than the horizontal deformation, and is the main component of surface deformation [7]. In addition, due to the flat terrain of the YND, the LOS deformation acquired from InSAR technology was simply interpreted as vertical deformations in this study [27]. Figure 6 shows the maps of the LOS deformation rates measured by SBAS from Sentinel-1A data stacks. The red color (negative values) in Figure 6 indicates that the coherent targets (CTs) are moving away from the satellite (i.e., subsidence), while the green color (positive values) indicates movement toward the satellite (i.e., uplift). Selected areas with concentrated deformation rates below −10 mm/year were analyzed in order to grasp the characteristics of the land subsidence in the YND. Figure 6 shows that the entire area of the old city remained relatively stable during the whole monitoring period, with deformation rates mainly in the range of −5 to 5 mm/year. However, there were approximately two regions located in the YND where severe subsidence was apparent, and the deformation rate was

Spatial Pattern
In this study, SBAS-InSAR technology was employed to quantify the spatiotemporal pattern of the surface deformation in the YND from December 2017 to December 2018. According to previous studies, the vertical deformation in the YND is much larger than the horizontal deformation, and is the main component of surface deformation [7]. In addition, due to the flat terrain of the YND, the LOS deformation acquired from InSAR technology was simply interpreted as vertical deformations in this study [27]. Figure 6 shows the maps of the LOS deformation rates measured by SBAS from Sentinel-1A data stacks. The red color (negative values) in Figure 6 indicates that the coherent targets (CTs) are moving away from the satellite (i.e., subsidence), while the green color (positive values) indicates movement toward the satellite (i.e., uplift). Selected areas with concentrated deformation rates below −10 mm/year were analyzed in order to grasp the characteristics of the land subsidence in the YND. Figure 6 shows that the entire area of the old city remained relatively stable during the whole monitoring period, with deformation rates mainly in the range of −5 to 5 mm/year. However, there were approximately two regions located in the YND where severe subsidence was apparent, and the deformation rate was primarily in the range of −10 to −60 mm/year. These two major subsidence zones with high deformation rates are located in the central Qiaoergou region and Gaojiagou region ( Figure 6). Field investigation photos Figure 6. Annual deformation rates in the Yan'an New District. The dashed line outlines the main deformation area. The black triangles correspond to the four typical sample areas used in the field surveys and the time series analysis (see Table  3 and Figure 7).

Temporal Evolution
The SBAS-InSAR method allowed for the temporal deformation evolution estimation for each of the CTs identified in the SAR imagery. Thus, we selected the four field survey points (P1 to P4) mentioned above on different parts of the YND, as well as a stable reference point P5 in the old city, as marked with triangles in Figure 6, in order to quantitatively analyze the temporal evolution of the ground deformation. The time series of the deformation from December 2017 to December 2018 measured by the Sentinel-1A data stacks is shown in Figure 7. A negative value indicates that the ground target was moving away from the satellite. The temporal evolution of the deformation at the four feature points (P1 to P4) shows similar near-linear trends with various velocities. P5-located in Figure 6. Annual deformation rates in the Yan'an New District. The dashed line outlines the main deformation area. The black triangles correspond to the four typical sample areas used in the field surveys and the time series analysis (see Table  3 and Figure 7).

Temporal Evolution
The SBAS-InSAR method allowed for the temporal deformation evolution estimation for each of the CTs identified in the SAR imagery. Thus, we selected the four field survey points (P1 to P4) mentioned above on different parts of the YND, as well as a stable reference point P5 in the old city, as marked with triangles in Figure 6, in order to quantitatively analyze the temporal evolution of the ground deformation. The time series of the deformation from December 2017 to December 2018 measured by the Sentinel-1A data stacks is shown in Figure 7. A negative value indicates that the ground target was moving away from the satellite. The temporal evolution of the deformation at the four feature points (P1 to P4) shows similar near-linear trends with various velocities. P5-located in Figure 6. Annual deformation rates in the Yan'an New District. The dashed line outlines the main deformation area. The black triangles correspond to the four typical sample areas used in the field surveys and the time series analysis (see Table  3 and Figure 7).

Temporal Evolution
The SBAS-InSAR method allowed for the temporal deformation evolution estimation for each of the CTs identified in the SAR imagery. Thus, we selected the four field survey points (P1 to P4) mentioned above on different parts of the YND, as well as a stable reference point P5 in the old city, as marked with triangles in Figure 6, in order to quantitatively analyze the temporal evolution of the ground deformation. The time series of the deformation from December 2017 to December 2018 measured by the Sentinel-1A data stacks is shown in Figure 7. A negative value indicates that the ground target was moving away from the satellite. The temporal evolution of the deformation at the four feature points (P1 to P4) shows similar near-linear trends with various velocities. P5-located in Figure 6. Annual deformation rates in the Yan'an New District. The dashed line outlines the main deformation area. The black triangles correspond to the four typical sample areas used in the field surveys and the time series analysis (see Table  3 and Figure 7).

Temporal Evolution
The SBAS-InSAR method allowed for the temporal deformation evolution estimation for each of the CTs identified in the SAR imagery. Thus, we selected the four field survey points (P1 to P4) mentioned above on different parts of the YND, as well as a stable reference point P5 in the old city, as marked with triangles in Figure 6, in order to quantitatively analyze the temporal evolution of the ground deformation. The time series of the deformation from December 2017 to December 2018 measured by the Sentinel-1A data stacks is shown in Figure 7. A negative value indicates that the ground target was moving away from the satellite. The temporal evolution of the deformation at the four feature points (P1 to P4) shows similar near-linear trends with various velocities. P5-located in  Table  3 and Figure 7).   Table 3 and Figure 7).
deformations were observed at P1 and P2, with magnitudes of tively. Notably, the subsidence at P1 to P4 still exhibited a lar deformation rate, which suggested that the YND will still unde sidence in the near future. Therefore, long-term, multiscale groun in order to further understand the detailed temporal evolution o

Relationships between the Topographic Changes and Land Subsid
As mentioned earlier, the large-scale land creation caused changes, forming large loess high-fill foundations with a maxim The results of the previous stratified settlement monitoring show the underlying original loess in the YND mainly occurred during while the postconstruction ground settlement mainly consisted [7]. Therefore, the distribution and thickness of the filling body h the distribution and magnitude of the land subsidence. The fil were spatially analyzed with InSAR-derived deformation inform 8. It is clear from Figure 8 that the InSAR-derived deformatio shows a high correlation between the detected subsidence and t loess. In areas with large difference values in fill thickness, sig differences can be identified. Moreover, relatively high subside were observed in the center of the Qiaoergou and Gaojiagou regi ized by a relatively greater thickness of filling loess. The Qiaoergou region (i.e., Region A1) is in the major urban area of the YND. In this region, the severe land subsidence is distributed in a zone from northwest to southeast along the filling areas, with deformation rates ranging from −10 to −45 mm/year. The deformation rates increased gradually from the northwest to southeast of the subsidence zone. The largest subsidence was detected at the southeastern end of this land subsidence zone, with a maximum deformation rate up to −73 mm/year, which is consistent with the results derived from a previous study, e.g., 70 mm/year in 2015-2019 [5]. The central area of the Gaojiagou region (i.e., Region A2) is the second largest land subsidence zone in the YND, with a deformation rate ranging from −10 to −75 mm/year and a maximum value of −121 mm/year. However, very sparse measurement points were detected in the southeastern end of the Gaojiagou area, which is most likely due to the decoherence caused by continuous and high-intensity engineering activities, such as land creation and rapid construction.
In addition to the land subsidence in the filling region, large-scale uplift was also observed in the mountain excavation region (i.e., Regions B1 and B2), mainly ranging from 10 to 20 mm/year, with a maximum deformation rate of 62 mm/year, which might be attributed to the released load of the excavated mountains [5]. Additionally, some slight to moderate deformation areas were also observed in the engineering slope in the western part of the YND (i.e., Regions C1 and C2), which was continuously deformed from December 2017 to December 2018. The main deformation rate ranged from −10 to −30 mm/year, with a maximum deformation rate of −38 mm/year.
Since high-precision ground-based levelling data were unavailable in this region, a series of detailed field investigations were conducted in April 2019 in order to verify the accuracy and reliability of the InSAR-derived deformation results. Table 3 shows the ground deformation characteristics obtained from the field survey based on the InSARderived deformation results. The locations of the field investigation photos are marked by black triangles in Figure 6. As seen in Table 3, a series of ground deformation occurrencessuch as road cracks and wall cracks-are distributed on the YND, which is consistent with the InSAR-derived results. For the Qiaoergou region, due to the gradual improvement of the infrastructure and major roads, the ground deformation characteristics in this region are more obvious than those in the Gaojiagou region, such as road pavement cracks, wall cracks, and surface drainage channel cracks (Table 3). However, in the Gaojiagou region, the completion time of the land creation project was relatively late, and there are almost no buildings; the area is mainly covered by woodland or wasteland. Therefore, although the ground deformation rate in this area is higher than that in the Qiaoergou region, the signs of deformation failure on the site are not obvious, and there is only obvious damage in the flood control ponds. Overall, good agreement between InSAR-derived results and field survey data was obtained, which suggests that the SBAS method can achieve reliable and accurate surface deformation measurements for land-creation areas. Additionally, the standard deviation of the InSAR results over a known nondeforming area (orange rectangle in Figure 6) was calculated as an indicator to further quantify the accuracy of our results. The standard deviation is approximately 4.3 mm/year, indicating that the InSAR measurements have reliable accuracy.

Temporal Evolution
The SBAS-InSAR method allowed for the temporal deformation evolution estimation for each of the CTs identified in the SAR imagery. Thus, we selected the four field survey points (P1 to P4) mentioned above on different parts of the YND, as well as a stable reference point P5 in the old city, as marked with triangles in Figure 6, in order to quantitatively analyze the temporal evolution of the ground deformation. The time series of the deformation from December 2017 to December 2018 measured by the Sentinel-1A data stacks is shown in Figure 7. A negative value indicates that the ground target was moving away from the satellite. The temporal evolution of the deformation at the four feature points (P1 to P4) shows similar near-linear trends with various velocities. P5-located in the old urban area-is basically stable, with a cumulative deformation within 2 mm during the whole monitoring period. The cumulative deformation of the four feature points (P1 to P4) exceeds 20 mm in only one year, which is enough to trigger and aggravate road pavement cracks and wall cracks ( Table 3). The maximum and minimum accumulative deformations were observed at P1 and P2, with magnitudes of −41 and −20 mm, respectively. Notably, the subsidence at P1 to P4 still exhibited a large subsidence trend and deformation rate, which suggested that the YND will still undergo significant land subsidence in the near future. Therefore, long-term, multiscale ground monitoring is required in order to further understand the detailed temporal evolution of the subsidence.

Relationships between the Topographic Changes and Land Subsidence
As mentioned earlier, the large-scale land creation caused significant topographic changes, forming large loess high-fill foundations with a maximum thickness of 110 m. The results of the previous stratified settlement monitoring showed that the settlement of the underlying original loess in the YND mainly occurred during the construction period, while the postconstruction ground settlement mainly consisted of fill body compression [7]. Therefore, the distribution and thickness of the filling body had a significant effect on the distribution and magnitude of the land subsidence. The fill thickness contour lines were spatially analyzed with InSAR-derived deformation information, as shown in Figure 8. It is clear from Figure 8 that the InSAR-derived deformation pattern over the YND shows a high correlation between the detected subsidence and the thickness of the filling loess. In areas with large difference values in fill thickness, significant land subsidence differences can be identified. Moreover, relatively high subsidence rates (>20 mm/year) were observed in the center of the Qiaoergou and Gaojiagou regions, which are characterized by a relatively greater thickness of filling loess. In order to further investigate the spatial correlation between the topographic changes and land subsidence, the elevation and deformation rates along profiles AA′ and BB' were extracted, as shown in Figure 9. Profile AA′ crosses the bottom of the Qiaoergou gully almost vertically. The subsidence along the profile increases with the increase in the fill thickness. The maximum deformation rate along the profile was approximately 45.6 mm/year, which is located in the middle of the profile with the largest thickness of filled loess. Profile BB′ extends from the middle to the top of the Qiaoergou gully in a northsouth direction. A similar phenomenon can be found in this profile, in which the maximum deformation is observed at the maximum filling loess thickness. In general, the subsidence rate is higher at positions with a greater thickness of filling loess, which indicates that the subsidence formed by the compression of the filling loess is relatively large under self-compression and external loading, because the filling loess is highly compressible. In order to quantitatively determine the correlation between the fill thickness and the land subsidence, the fill thickness at each CT in the Qiaoergou region was extracted and then statistically analyzed. Figure 10a presents the percentage of CTs with different deformation rates in different fill thicknesses. Almost all of the ground subsidence (94.7% of the total CTs with deformation rates below −10 mm/year) is concentrated in the area where the thickness of the fill is greater than 30 m. Moreover, as the thickness of the fill increases, the proportion of CTs with higher subsidence rates increases significantly; in particular, subsidence at rates exceeding 50 mm/year was only observed in the area where the thickness of the fill is greater than 70 m. In contrast, no subsidence or minor subsidence was Figure 9. Comparison of the subsidence with the thickness of filling loess along two profiles of which the positions are marked in Figure 8. A negative deformation means subsidence.
In order to further investigate the spatial correlation between the topographic changes and land subsidence, the elevation and deformation rates along profiles AA and BB were extracted, as shown in Figure 9. Profile AA crosses the bottom of the Qiaoergou gully almost vertically. The subsidence along the profile increases with the increase in the fill thickness. The maximum deformation rate along the profile was approximately 45.6 mm/year, which is located in the middle of the profile with the largest thickness of filled loess. Profile BB extends from the middle to the top of the Qiaoergou gully in a northsouth direction. A similar phenomenon can be found in this profile, in which the maximum deformation is observed at the maximum filling loess thickness. In general, the subsidence rate is higher at positions with a greater thickness of filling loess, which indicates that the subsidence formed by the compression of the filling loess is relatively large under self-compression and external loading, because the filling loess is highly compressible.
In order to quantitatively determine the correlation between the fill thickness and the land subsidence, the fill thickness at each CT in the Qiaoergou region was extracted and then statistically analyzed. Figure 10a presents the percentage of CTs with different deformation rates in different fill thicknesses. Almost all of the ground subsidence (94.7% of the total CTs with deformation rates below −10 mm/year) is concentrated in the area where the thickness of the fill is greater than 30 m. Moreover, as the thickness of the fill increases, the proportion of CTs with higher subsidence rates increases significantly; in particular, subsidence at rates exceeding 50 mm/year was only observed in the area where the thickness of the fill is greater than 70 m. In contrast, no subsidence or minor subsidence was observed in areas with low fill thickness, especially below 20 m. Figure 10b presents the average subsidence rate for different fill thicknesses. A positive correlation between the subsidence rate and the fill thickness was observed; that is, the subsidence rate increases significantly with increasing filling loess thicknesses. These phenomena suggest that the different magnitudes of land subsidence after land creation are largely controlled by the compaction process of the filling loess with different thicknesses.  Figure 11 shows the overlap between the ground deformation, the land-creation area, and the construction area. Clearly, there was an obvious high spatial correlation between the surface deformation and the land creation, indicating a high level of spatial control of the subsidence by the land creation. The density of the total CTs distributed in the nonfilled area is 1482/km 2 , while it decreased to 1191/km 2 in the filled area, which is mainly caused by the SAR image resolution and continuous engineering activities. However, the density of the CTs distributed in the filled area with deformation rates below −10 mm/year (indicating subsidence) reached 616/km 2 , while it decreased to 163/km 2 in the non-filled area.

Effect of Large-Scale Land Creation and Urban Construction
In addition, more than 81% of the total CTs with deformation rates below −10 mm/year (indicating subsidence) are concentrated in the infilling area, with a maximum deformation rate of −121 mm/year, which indicates that most of the land subsidence in the   Figure 11 shows the overlap between the ground deformation, the land-creation area, and the construction area. Clearly, there was an obvious high spatial correlation between the surface deformation and the land creation, indicating a high level of spatial control of the subsidence by the land creation. The density of the total CTs distributed in the non-filled area is 1482/km 2 , while it decreased to 1191/km 2 in the filled area, which is mainly caused by the SAR image resolution and continuous engineering activities. However, the density of the CTs distributed in the filled area with deformation rates below −10 mm/year (indicating subsidence) reached 616/km 2 , while it decreased to 163/km 2 in the non-filled area.
Remote Sens. 2021, 13, x FOR PEER REVIEW 17 of 21 Figure 11. Land creation and urban construction in relation to the subsidence in the YND. The rectangular box is enlarged and shown in Figure 12.

Relationship between the Land Subsidence and Road Network
Notably, the observed subsidence can also be affected by effects relevant to road construction. As of June 2018, the total length of the road network in the YND exceeded 80 km (Figure 11), with a distribution density of 3.64 km/km 2 , which had a significant impact on the development of the subsidence. Figure 12a shows that there is a strong spatial correlation between the observed subsidence and the dense road network. Severe land subsidence was observed along the road networks, which was consistent with road pavement cracks during the field survey (Table 3). Thus, subsidence is a major threat to vehicle safety.
For a further analysis of the spatial distribution relationship between the road network and the land subsidence, CTs at different distances from the center of the road were selected for the Qiaoergou filling area as an example, and then a statistical analysis of the subsidence rates corresponding to these CTs was executed. Figure 12b indicates the percentage of the number of largely-subsiding CTs compared to the total number of CTs within different distance ranges. In Figure 12b, the proportion of CTs with higher subsidence rates to the total number of CTs increases significantly as the distance from the center of the road decreases, indicating that the land subsidence around the center of the road is more significant. On the contrary, relatively slight land subsidence occurs only in the areas with larger distances from the center of the road, especially those larger than 150 m, suggesting that the ground deformation in this area is mainly caused by the natural consolidation process of the filled loess. This result indicates that the land subsidence in the YND is also somewhat affected by road construction. The long-term cyclic traffic load and its vibration load due to the dense road network will accelerate the consolidation process of the filling loess, thereby significantly exacerbating the progression of the land subsidence.  Figure 11. The background is a panchromatic ZY-3 image. (b) Statistics of the CTs at different distances from the road.

Natural Drivers of Subsidence
A complex geological environment is often an important factor that impacts land subsidence [41]. As mentioned earlier, almost all of the land subsidence occurred in the infilling area. The compressible and collapsible deformation of the filled loess may be the potential mechanism of large-scale land subsidence in the YND. Large-scale land creation, including the excavation of loess hills to fill valleys, will potentially modify the geotech- In addition, more than 81% of the total CTs with deformation rates below −10 mm/year (indicating subsidence) are concentrated in the infilling area, with a maximum deformation rate of −121 mm/year, which indicates that most of the land subsidence in the YND occurred in the filling area. On the other hand, there is a relatively weak correlation in the spatial distribution between the land subsidence and buildings. Generally, only 15.86% of the total CTs indicating subsidence are distributed in the building area, whereas 54.82% of these buildings are localized within the filled area. Compared with land creation, building construction has a significantly smaller impact on ground settlement, mainly because most of the building construction is located in relatively stable excavation areas. Therefore, it is not difficult to see that the spatial pattern of the land subsidence over the YND is bounded and controlled by large-scale land creation. Filled loess is formed by the remodeling and backfilling of original loess during filling work, which has a loose structure and weak mechanical properties. Therefore, in the land-creation area, especially the filling area, the loess-filled foundation is prone to compaction under its own weight and external load, which finally leads to significant subsidence.

Relationship between the Land Subsidence and Road Network
Notably, the observed subsidence can also be affected by effects relevant to road construction. As of June 2018, the total length of the road network in the YND exceeded 80 km (Figure 11), with a distribution density of 3.64 km/km 2 , which had a significant impact on the development of the subsidence. Figure 12a shows that there is a strong spatial correlation between the observed subsidence and the dense road network. Severe land subsidence was observed along the road networks, which was consistent with road pavement cracks during the field survey (Table 3). Thus, subsidence is a major threat to vehicle safety.
For a further analysis of the spatial distribution relationship between the road network and the land subsidence, CTs at different distances from the center of the road were selected for the Qiaoergou filling area as an example, and then a statistical analysis of the subsidence rates corresponding to these CTs was executed. Figure 12b indicates the percentage of the number of largely-subsiding CTs compared to the total number of CTs within different distance ranges. In Figure 12b, the proportion of CTs with higher subsidence rates to the total number of CTs increases significantly as the distance from the center of the road decreases, indicating that the land subsidence around the center of the road is more significant. On the contrary, relatively slight land subsidence occurs only in the areas with larger distances from the center of the road, especially those larger than 150 m, suggesting that the ground deformation in this area is mainly caused by the natural consolidation process of the filled loess. This result indicates that the land subsidence in the YND is also somewhat affected by road construction. The long-term cyclic traffic load and its vibration load due to the dense road network will accelerate the consolidation process of the filling loess, thereby significantly exacerbating the progression of the land subsidence.

Natural Drivers of Subsidence
A complex geological environment is often an important factor that impacts land subsidence [41]. As mentioned earlier, almost all of the land subsidence occurred in the infilling area. The compressible and collapsible deformation of the filled loess may be the potential mechanism of large-scale land subsidence in the YND. Large-scale land creation, including the excavation of loess hills to fill valleys, will potentially modify the geotechnical properties (such as the collapsibility and compressibility) of the original loess [38]. Therefore, in the infilling area, the infilling loess has the properties of weak cementation and a relatively loose structure [42], which makes it easy to experience substantial compression under the action of self-weight or an external load. In addition, many studies have shown that repeated loess compaction in layers during the filling work is not enough to completely remove the threat of the loess foundation's collapse. Conversely, the collapsibility of the compacted loess is even more serious than that of intact loess [38,43]. Therefore, under external loads (engineering vehicles and heavy buildings) and wetting (rainfall or irrigation), filling loess is more prone to compaction and collapsible deformation, ultimately leading to substantial land subsidence. In summary, considering the significant influence of the surface geological characteristics on the subsidence rate and magnitude, the spatial variability in filled loess in terms of thickness and lithological characteristics of soils largely determines the spatiotemporal pattern of the land subsidence in the YND [12,44].

Implications and Limitations
In this study, the InSAR-derived results showed that YND suffered a maximum subsidence deformation of more than 120 mm in only one year, which posed a great challenge to the urban construction of the YND. Severe land subsidence not only threatens the safety of buildings and infrastructures but also greatly restricts the application of the MECC project in the YND, and even affects the sustainable development of the city. Moreover, given the large-scale land creation, increasing urbanization of the YND and the related increase in external loads from high-rise buildings and traffic networks, the infilling loess is more prone to compaction and, eventually, accelerated subsidence [38]. Accordingly, the extensive research efforts on the postconstruction subsidence estimation and prediction should be continued and accelerated in the aspects of geological interpretation, geotechnical investigation, ground monitoring and modeling, and InSAR observation. Furthermore, InSAR analysis has been demonstrated to be a valuable tool for mapping land subsidence, which should be taken into account in future urbanization planning and technical solutions involved in avoiding infrastructure damage caused by subsidence [27]. Finally, the ongoing monitoring of land-creation areas is needed in order to update the information on the compression or consolidation of filling loess, and to identify early critical situations that may cause damage to the urban infrastructures in these areas [27]. Fortunately, our results indicate that the planning and construction in the YND are generally consistent with the rules for land creation construction in the YND, where construction work will be carried out first in the mountain excavation region and will not expand into the infilled areas until the ground stabilizes by natural settlement [45].
In this paper, the rates and patterns of postconstruction land deformation associated with MECC in the YND were investigated using advanced time series InSAR analysis combined with topographic change and optical remote sensing data; however, several limitations remain. These limitations consist in the lack of valid in situ measurements, such as levelling data and the latest monitoring well data, which would provide a powerful validation of the subsidence rates and patterns, and would be of assistance in urban planning and land subsidence mitigation.
In the near future, various in situ measurement data, such as ground and deep monitoring data, groundwater level changes, and geotechnical data, will be collected, and a further joint analysis of multidisciplinary data will be conducted in order to comprehensively understand the complex subsidence phenomena in the YND. This action will enable us to thoroughly understand the spatiotemporal evolution and mechanism of land subsidence in order to ensure the reliability of the long-term subsidence predictions and thereby assist decision making on land subsidence mitigation.

Conclusions
In this study, the local topographic changes and ground deformation caused by the large-scale MECC project in the YND, as well as their interactions, were investigated by jointly using space-borne optical photogrammetry, time series InSAR techniques, remote sensing images, and field surveys. First, new surface elevations were extracted from ZY-3 stereo images using space-borne optical photogrammetry, and then the topographic changes and fill thickness associated with the MECC project between 2000 and 2017 were derived from the ZY-3 DEM and SRTM DEM using the geodetic method. By October 2017, over 22 km 2 of flat land, including 10.8 km 2 of filled area, was created, with a maximum filling depth of approximately 110 m. Second, the ground deformation after the MECC project was retrieved by applying the SBAS-InSAR technique using the Sentinel-1 SAR dataset. Extensive land subsidence, with a maximum settlement rate of 121 mm/year, was detected in the land-creation area, which is consistent with the detailed field findings. Finally, the quantitative analysis of various driving forces of subsidence was conducted by combining the InSAR results with engineering activities, local geological data, and remote sensing data. Large-scale land creation determines the spatial and temporal patterns of land subsidence, while the thickness of the filled loess controls the magnitude of the land subsidence, and human activities, as well as changes in the geological environment, will also accelerate the process of subsidence. These results demonstrate that long-term, multiscale deformation monitoring is needed in order to update the information regarding the compression or consolidation of filling loess, thus ensuring the safe operation and sustainability of the MECC project.