SBAS-InSAR Based Deformation Detection of Urban Land, Created from Mega-Scale Mountain Excavating and Valley Filling in the Loess Plateau: The Case Study of Yan’an City

: Yan’an new district (YND) is one of the largest civil engineering projects for land creation in Loess Plateau, of which the amount of earthwork exceeds 600 million m 3 , to create 78.5 km 2 of ﬂat land. Such mega-scale engineering activities and complex geological characteristics have induced wide land deformation in the region. Small baseline subset synthetic aperture radar interferometry (SBAS-InSAR) method and 55 Sentinel-1A (S-1A) images were utilized in the present work to investigate the urban surface deformation in the Yan’an urban area and Yan’an new airport (YNA) from 2015 to 2019. The results were validated by the ground leveling measurements in the YNA. It is found that signiﬁcant uneven surface deformation existed in both YND and YNA areas with maximum accumulative subsidence of 300 and 217 mm, respectively. Moreover, the average subsidence rate of the YND and YNA areas ranged from − 70 to 30 mm / year and − 50 to 25 mm / year, respectively. The present work shows that the land deformation su ﬀ ered two periods (from 2015 to 2017 and from 2017 to 2019) and expanded from urban center to surrounding resettlement area, which are highly relevant with urban earthwork process. It is found that more than 60% of land subsidence occurs at ﬁlled areas, while more than 65% of surface uplifting occurs at excavation areas. The present work shows that the subsidence originates from the earth ﬁlling and the load of urban buildings, while the release of stress is the major factor for the land uplift. Moreover, it is found that the collapsibility of loess and concentrated precipitation deteriorates the degree of local land subsidence. The deformation discovered by this paper shows that the city may su ﬀ er a long period of subsidence, and huge challenges may exist in the period of urban maintaining buildings and infrastructure facilities.


Introduction
Urban land subsidence is a widespread phenomenon associated with the rapid urbanization around the world [1][2][3]. Different reasons, including reclamation in coastal areas [4][5][6], industrial development and urban construction in delta areas [7][8][9] and the overexploitation of groundwater in plains areas [10,11] may cause urban ground deformation. Land subsidence is a complex geotechnical disaster with specific characteristics such as slow formation, long duration and irreversibility [12,13]. It should be indicated that land subsidence in large areas might cause settling funnels, seawater backflow, surface collapse and crustal displacement [14]. Moreover, land subsidence in urban areas could cause foundation collapse, underground pipeline destruction and road damage [15]. Therefore, the anticipation, detection and mitigation of urban land subsidence is of great significance [16].

Study Area
Yan'an city (36 • 10 33"N to 37 • 2 5"N, 109 • 14 10"E to 110 • 30 43"E) is located at the middle region of the Yellow River in the central and southern areas of the Loess Plateau, China. Yan'an has a continental monsoon climate. The annual average precipitation in Yan'an is about 560 mm, where more than 70% of the total annual precipitation is concentrated in June to September [42]. Moreover, most of rainstorms occur in August, which is one of the main factors of the soil erosion in the study area [43]. For example, the extreme rainfall in July 2013 caused serious loess collapse, loess landslide, land subsidence and other geo-hazards in Yan'an [44]. The main landscapes in Yan'an are hills and valleys with an average elevation of 1200 m. Late Pleistocene aeolian Malan loess and Holocene alluvial loess with strong collapsibility are widely distributed, which are the primary deposits causing engineering construction complexity [45]. Figure 1 shows the old Yan'an urban area located at the intersection of three mountains (Phoenix Hill, Qingliang Hill and Pagoda Hill) and two rivers (the Yanhe and Nanchuan rivers). The old city has a typical Y-shaped structure with a total length of over 20 km, but is just hundreds of meters wide, which causes available land to be an extremely limited resource among a huge population. Furthermore, the new district (YND) is located at the north of Qingliang Hill, covering an area about 5 km from south to north, and about 6 km from west to east, which has been established by the local government for anticipated urban developments.
Remote Sens. 2019, 11, x FOR PEER REVIEW 3 of 20 Yan'an city (36°10′33″N to 37°2′5″N, 109°14′10″E to 110°30′43″E) is located at the middle region of the Yellow River in the central and southern areas of the Loess Plateau, China. Yan'an has a continental monsoon climate. The annual average precipitation in Yan'an is about 560 mm, where more than 70% of the total annual precipitation is concentrated in June to September [42]. Moreover, most of rainstorms occur in August, which is one of the main factors of the soil erosion in the study area [43]. For example, the extreme rainfall in July 2013 caused serious loess collapse, loess landslide, land subsidence and other geo-hazards in Yan'an [44]. The main landscapes in Yan'an are hills and valleys with an average elevation of 1200 m. Late Pleistocene aeolian Malan loess and Holocene alluvial loess with strong collapsibility are widely distributed, which are the primary deposits causing engineering construction complexity [45]. Figure 1 shows the old Yan'an urban area located at the intersection of three mountains (Phoenix Hill, Qingliang Hill and Pagoda Hill) and two rivers (the Yanhe and Nanchuan rivers). The old city has a typical Y-shaped structure with a total length of over 20 km, but is just hundreds of meters wide, which causes available land to be an extremely limited resource among a huge population. Furthermore, the new district (YND) is located at the north of Qingliang Hill, covering an area about 5 km from south to north, and about 6 km from west to east, which has been established by the local government for anticipated urban developments. The construction of YND is the largest geotechnical project in a collapsible loess area around the world, and the volume of this earthwork is twice of that of the Three Gorges Dam. The earthwork was commenced April 2012, and was partly completed by the end of June 2015, with remaining earthwork postponed due to unpredictable environmental problems. Based on the official reports, the real amount of earthwork was about 600 million m 3 , of which about 330 million and 270 million m 3 were related to excavating and filling works, respectively. The maximum depth of filling exceeded 112 m, and the real created flat land for the development and construction was approximately 22.7 km 2 , including 10.4 km 2 of excavating and 12.3 km 2 of filling ( Figure 2) [41].
Yan'an new airport is located in Liulin Town, Baota District, 13 km southwest of Yan'an city center. It was built between August 2013 and November 2018, with a track 3 km long and 45 m wide. It should be indicated that, similar to YND, the land for the new airport was also acquired by mountain excavating and valley filling [46]. The construction of YND is the largest geotechnical project in a collapsible loess area around the world, and the volume of this earthwork is twice of that of the Three Gorges Dam. The earthwork was commenced April 2012, and was partly completed by the end of June 2015, with remaining earthwork postponed due to unpredictable environmental problems. Based on the official reports, the real amount of earthwork was about 600 million m 3 , of which about 330 million and 270 million m 3 were related to excavating and filling works, respectively. The maximum depth of filling exceeded 112 m, and the real created flat land for the development and construction was approximately 22.7 km 2 , including 10.4 km 2 of excavating and 12.3 km 2 of filling ( Figure 2) [41].
Yan'an new airport is located in Liulin Town, Baota District, 13 km southwest of Yan'an city center. It was built between August 2013 and November 2018, with a track 3 km long and 45 m wide. It should be indicated that, similar to YND, the land for the new airport was also acquired by mountain excavating and valley filling [46].

Sentinel-1A Data
In order to obtain a comprehensive and detailed surface deformation of the research region, 55 Sentinel-1A images are utilized from the period of February 2015 to January 2019, with a sampling interval of approximately 30 days. The mean and maximum absolute baseline are 51 and 121 m, respectively. Table 1 shows the specific information of the data set. Moreover, precise orbital ephemeris data and global terrain data (shuttle radar topography mission (SRTM)) are applied into the orbit refinement and terrain phase removal, respectively.

Basic Theory of SBAS-InSAR
The SBAS-InSAR technique is an analysis method for multi-master-image InSAR time series that utilizes interferometric pairs of small time-space baselines to generate the surface deformation. This technique overcomes the low coherence of some interferograms induced by single super-master images and reduces the requirement of the quantity of SAR images with respect to PS-InSAR [29]. SBAS-InSAR is widely applied in land subsidence detection [26,27,31], and is based on the following equations [25]: Equation (1) represents the quantity range of M differential interferograms generated by N + 1 SAR images at the same area with ordered time (t 0 , t 1 · · · , t N ). Moreover, Equation (2) represents the composition of the interferometric phase of interferogram j (generated by two images at t B and t A ) in the pixel (x, r), where x and r are the azimuth and the range coordinate, respectively. Phase ∆φ disp is induced by the change of distance between target and radar along the line of sight (LOS). Furthermore, phases ∆φ topo , ∆φ orb , ∆φ atm and ∆φ noise are generated by the terrain, satellite orbit error, atmosphere effect (especially the effect of tropospheric delay) and other noise, respectively. SBAS-InSAR removes the residual components from interferometric phase to obtain deformation phase ∆φ disp . Then, a system of M equations with N unknowns can be organized as the following matrix: where A donates an M × N coefficient matrix, φ = [φ(t 1 ), φ(t 2 ), · · · , φ(t N )] T is the N × 1 vector of unknown deformation phase values with measurement points, and δφ = [δφ(t 1 ), δφ(t 2 ), · · · , δφ(t N )] T is the N × 1 vector of unwrapped phase values. Deformation velocity can be obtained by rewriting Equation (3) into the following equation: where B is an M × N coefficient matrix and v donates the mean phase velocity as follows: Finally, the least square (LS) or singular value decomposition (SVD) algorithms can be applied to Equations (4) and (5) to calculate the deformation velocity [31,47]. The cumulative deformation can be generated by performing the calculation on aforementioned equations over a specific time span.

Data Processing
For SBAS-InSAR process, primary errors come from data quality (orbit accuracy, coherence quality, etc.), atmosphere phase error and phase unwrapping. The data processing procedure can be divided into four parts as follow: • The first part is the basis of the SBAS-InSAR. By selecting a definite temporal-spatial baseline threshold (time baseline 120 days, space baseline 121 m), the data sets were divided into different subsets and 201 interferometric pairs were generated ( Figure 3).

•
The second part is the interferometric flow step. Firstly, the external DEM data was utilized for auxiliary co-registration with accuracy of the sub-pixel level [48]. Then, the differential interferometry process was conducted for each pair to obtain interferograms and a coherence coefficient map. In addition, we removed the topographical phase by external DEM data, eliminated the noise phase using a Goldstein filter [49] and unwrapped the phase using the minimum cost flow (MCF) method [50]. Finally, Equation (6) was employed to calculate the theoretical precision with the interferometric coherence γ and the wavelength λ.
The third part is the key step in generating surface deformation. In this research, ground control points (GCPs) were exploited to eliminate the change of phase caused by satellite orbit errors and remove the residual topographic phase. Once the quality inspection was approved, the linear model and residual phase were used to estimate the surface deformation. The effect of atmosphere delay was alleviated by a spatial-temporal filter (time filter window set to 365 days, space filter window set to 1800 m) to obtain the accumulated land deformation and average deformation rate along the LOS direction [51].

•
The forth part is the process of generating land subsidence (vertical deformation). It is assumed that surface deformation mainly originated by vertical components while ignoring the influence of horizontal components. Vertical deformation S was calculated via combination with incident angle θ i and displacement d in the LOS direction, based on the relationship shown as Equation (7) [52]. Finally, the surface displacement data were geocoded to obtain the result in the WGS84 coordinate system.

Results
The SBAS-InSAR achieved favorable coverage in the study area, with theoretical precision of 0-1.8 mm/year, 0-2.5 mm/year and 0-4 mm/year in the old urban, YND and YNA areas, respectively. Figure 4a indicates that the whole ground surface of the old urban area from February 2015 to January 2019 maintained stability, and the deformation rate was primarily in the range of −6 to 7 mm/year. Moreover, Figure 4b shows severe surface deformation within the YND region, varying from −70 to 30 mm/year. On the other hand, uneven surface deformation with a range of −50 to 25 mm/year was found in the south and north end of the runway in the YNA area. It should be indicated that the positive and negative values represent land uplift and subsidence, respectively. In addition, a precise ground leveling survey in YNA was utilized to validate obtained results from the SBAS-InSAR scheme, where 7.4 mm/year of the root mean square error (RMSE) was calculated. Figure 4 illustrates the results of the surface deformation in YND. The major land sinking was concentrated in filling areas along valleys with average and maximum subsidence rates of 20 and 70 mm/year, respectively. There was a northwest to southeast subsidence zone along the area between Shanghai Road and Xingzhi Road, which was spread between Yan'an Grand Theater and Yan'an People Park. The most significant subsidence occurred in the southeast end of the subsidence zone with cumulative land subsidence reaching 300 mm. The majority of land uplifting was performed at the mountain-excavated area, where the average and maximum uplift rates were 15 and 30 mm/year, respectively. The most significant uplifting was in the residential area on the northwest of YND and the surroundings of Yan'an Learning Academic, with cumulative uplift of 100 mm. Mega-scale earthwork of serious collapsible loess might be the main reason contributing to this phenomenon.

Results
The SBAS-InSAR achieved favorable coverage in the study area, with theoretical precision of 0-1.8 mm/year, 0-2.5 mm/year and 0-4 mm/year in the old urban, YND and YNA areas, respectively. Figure 4a indicates that the whole ground surface of the old urban area from February 2015 to January 2019 maintained stability, and the deformation rate was primarily in the range of −6 to 7 mm/year. Moreover, Figure 4b shows severe surface deformation within the YND region, varying from −70 to 30 mm/year. On the other hand, uneven surface deformation with a range of −50 to 25 mm/year was found in the south and north end of the runway in the YNA area. It should be indicated that the positive and negative values represent land uplift and subsidence, respectively. In addition, a precise ground leveling survey in YNA was utilized to validate obtained results from the SBAS-InSAR scheme, where 7.4 mm/year of the root mean square error (RMSE) was calculated. Figure 4 illustrates the results of the surface deformation in YND. The major land sinking was concentrated in filling areas along valleys with average and maximum subsidence rates of 20 and 70 mm/year, respectively. There was a northwest to southeast subsidence zone along the area between Shanghai Road and Xingzhi Road, which was spread between Yan'an Grand Theater and Yan'an People Park. The most significant subsidence occurred in the southeast end of the subsidence zone with cumulative land subsidence reaching 300 mm. The majority of land uplifting was performed at the mountain-excavated area, where the average and maximum uplift rates were 15 and 30 mm/year, respectively. The most significant uplifting was in the residential area on the northwest of YND and the surroundings of Yan'an Learning Academic, with cumulative uplift of 100 mm. Mega-scale earthwork of serious collapsible loess might be the main reason contributing to this phenomenon.  Figures 5 shows that a significant subsidence zone developed from the northwest towards southeast between the Xingzhi Road and Shanghai Road. Figure 5a shows that the most severe subsidence is located at narrow strips between the Xingzhi Road and the Shanghai Road, close to the Jinming Street (region A). It is observed that the maximum and average cumulative subsidence are 306 and 180 mm, respectively. Moreover, Figure 5b indicates that the velocity of subsidence decreases annually. The average rate of subsidence was approximately 100 mm/year from 2015 to 2017, while it was approximately 50 mm/year from 2017 to 2019. Furthermore, Figure 5a shows another severe subsidence region along the subsidence zone at the other end of the narrow strips, close to Zichang Road (region B). It is found that the maximum and average cumulative subsidence are 195 and 150 mm, respectively. Moreover, both maximum and average accumulations increase linearly with the maximum rate of 50 mm/year.

Spatial Distribution of the Land Subsidence
The subsidence zone is located in the central area of YND, where the main earth filling operation was carried out. The filling height is dozens of meters, with a maximum height of nearly 100 m. Figure  6 shows that there a negative correlation exists between the subsidence and elevation. Furthermore, it is observed that many residential and commercial buildings are built along the subsidence zone. Thus, human activities, including land creation and urban construction, might be the primary human factors contributing to local severe subsidence.  Figure 5 shows that a significant subsidence zone developed from the northwest towards southeast between the Xingzhi Road and Shanghai Road. Figure 5a shows that the most severe subsidence is located at narrow strips between the Xingzhi Road and the Shanghai Road, close to the Jinming Street (region A). It is observed that the maximum and average cumulative subsidence are 306 and 180 mm, respectively. Moreover, Figure 5b indicates that the velocity of subsidence decreases annually. The average rate of subsidence was approximately 100 mm/year from 2015 to 2017, while it was approximately 50 mm/year from 2017 to 2019. Furthermore, Figure 5a shows another severe subsidence region along the subsidence zone at the other end of the narrow strips, close to Zichang Road (region B). It is found that the maximum and average cumulative subsidence are 195 and 150 mm, respectively. Moreover, both maximum and average accumulations increase linearly with the maximum rate of 50 mm/year.
The subsidence zone is located in the central area of YND, where the main earth filling operation was carried out. The filling height is dozens of meters, with a maximum height of nearly 100 m. Figure 6 shows that there a negative correlation exists between the subsidence and elevation. Furthermore, it is observed that many residential and commercial buildings are built along the subsidence zone. Thus, human activities, including land creation and urban construction, might be the primary human factors contributing to local severe subsidence.   Figure 7a shows that there is an apparent uplifting trend in the north of the Yan'an Learning Academic (region C), where the large-scale mountain excavation project was launched. Figure 7b indicates that the maximum and average cumulative uplifts are 85 and 58 mm, respectively. Moreover, it is found that the uplift reached to 62 mm from February 2015 to October 2016, and a significant lifting of 30 mm is observed from May 2017 to August 2017, which might be attributed to the building construction. It is illustrated that the ground surface shows a steady trend during the interval of these two periods and after August 2017. The most significant land uplift is detected at the   Figure 7a shows that there is an apparent uplifting trend in the north of the Yan'an Learning Academic (region C), where the large-scale mountain excavation project was launched. Figure 7b indicates that the maximum and average cumulative uplifts are 85 and 58 mm, respectively. Moreover, it is found that the uplift reached to 62 mm from February 2015 to October 2016, and a significant lifting of 30 mm is observed from May 2017 to August 2017, which might be attributed to the building construction. It is illustrated that the ground surface shows a steady trend during the interval of these two periods and after August 2017. The most significant land uplift is detected at the    Figure 8 shows the evolution of cumulative surface deformation in the YND area from February 2015 to January 2019. It should be noted that two significantly different periods were observed. From 2015 to 2017, there were rapid land subsidence and uplift with maximum rates of 100 and 50 mm/year, respectively. The deformation region gradually expanded from the urban center to the surrounding areas, and the major land subsidence was concentrated in the subsidence zone ( Figure  5a), while the major uplift was in the area surrounding the Yan'an Learning Academic (Figure 7a). The deformation might be attributable to the compressible and collapsible deformation, during which the air among the soil would be gradually discharged as the soil densifies. From 2017 to 2019, several new resettled residential quarters around YND experienced obvious land subsidence, and the urban center suffered, with maximum rate of 50 mm/year. The major earthwork was finished in the middle of 2015, and the deformation experienced a two-year uneven subsidence, with a decreasing rate since the end of 2017 and an uplift that was nearly zero during the same time, meaning that the load release had completed.  Figure 8 shows the evolution of cumulative surface deformation in the YND area from February 2015 to January 2019. It should be noted that two significantly different periods were observed. From 2015 to 2017, there were rapid land subsidence and uplift with maximum rates of 100 and 50 mm/year, respectively. The deformation region gradually expanded from the urban center to the surrounding areas, and the major land subsidence was concentrated in the subsidence zone (Figure 5a), while the major uplift was in the area surrounding the Yan'an Learning Academic (Figure 7a). The deformation might be attributable to the compressible and collapsible deformation, during which the air among the soil would be gradually discharged as the soil densifies. From 2017 to 2019, several new resettled residential quarters around YND experienced obvious land subsidence, and the urban center suffered, with maximum rate of 50 mm/year. The major earthwork was finished in the middle of 2015, and the deformation experienced a two-year uneven subsidence, with a decreasing rate since the end of 2017 and an uplift that was nearly zero during the same time, meaning that the load release had completed.    Figure 9a shows that the deformation velocity map derived from SBAS-InSAR method covers most area of the runway and the surrounding slopes in YNA. The obtained results show that there was uneven surface deformation between the north-south direction of slopes and runway from February 2015 to January 2019. Figure 9b indicates that the most significant subsidence occurred in the northeast of runway, where large amounts of filling loess was applied to fill valleys. On the other hand, Figure 9c indicates that the maximum and average cumulative subsidence amounts during the study period were 217 and 144 mm, respectively. Moreover, a linear trend of 50 mm/year is observed. It is found that the subsidence at the both ends of the runway is 30 mm/year. The west slope slightly uplifts with a maximum rate of 20 mm/year, where the mountains were excavated. The local uneven surface deformation may be attributable to the instability of loess foundation in the filling area and the released load of excavated mountains.  Figure 9a shows that the deformation velocity map derived from SBAS-InSAR method covers most area of the runway and the surrounding slopes in YNA. The obtained results show that there was uneven surface deformation between the north-south direction of slopes and runway from February 2015 to January 2019. Figure 9b indicates that the most significant subsidence occurred in the northeast of runway, where large amounts of filling loess was applied to fill valleys. On the other hand, Figure 9c indicates that the maximum and average cumulative subsidence amounts during the study period were 217 and 144 mm, respectively. Moreover, a linear trend of 50 mm/year is observed. It is found that the subsidence at the both ends of the runway is 30 mm/year. The west slope slightly uplifts with a maximum rate of 20 mm/year, where the mountains were excavated. The local uneven surface deformation may be attributable to the instability of loess foundation in the filling area and the released load of excavated mountains.

Validation of SBAS-InSAR Results
Twenty-one highly precise leveling points in YNA ( Figure 10) were used to verify the validity of SBAS-InSAR results. This research compared the results of SBAS and leveling data at these leveling points. It should be noted that leveling results were derived from differences between the last and the first measurement, while SBAS results were extracted from the pixels, including the location of the leveling point. Figure 11 shows a comparison of the two sets from August 2015 to June 2016. Results show that the two data sets have reasonable consistency. Moreover, the RMSE and correlation coefficient value of the two sets are 5.9 mm (7.4 mm/year) and 0.96, respectively. The largest and mean differences between the two sets are 12.4 mm and 3.8 mm, respectively. Although leveling data of subsidence points in YND are absent, the data source and deformation process in YND and YNA are consistent, which indicates that the SBAS-InSAR monitoring in the YND also yielded relatively reliable results.

Validation of SBAS-InSAR Results
Twenty-one highly precise leveling points in YNA ( Figure 10) were used to verify the validity of SBAS-InSAR results. This research compared the results of SBAS and leveling data at these leveling points. It should be noted that leveling results were derived from differences between the last and the first measurement, while SBAS results were extracted from the pixels, including the location of the leveling point. Figure 11 shows a comparison of the two sets from August 2015 to June 2016. Results show that the two data sets have reasonable consistency. Moreover, the RMSE and correlation coefficient value of the two sets are 5.9 mm (7.4 mm/year) and 0.96, respectively. The largest and mean differences between the two sets are 12.4 mm and 3.8 mm, respectively. Although leveling data of subsidence points in YND are absent, the data source and deformation process in YND and YNA are consistent, which indicates that the SBAS-InSAR monitoring in the YND also yielded relatively reliable results.

Human Activities Contributing to the Land Subsidence
In the aforementioned sections, it was found that land creation and urban construction in Yan'an induced a wide area of surface deformation. Land creation has destroyed the original loess structure of the local area, which has caused severe instabilities in filling loess foundation. In addition, large numbers of buildings have increased the external pressure on surface and accelerated the compressible deformation of the loess layer and building foundations. Figure 12 shows the overlap between the surface deformation area and earthwork region. It is indicated that there is a high correlation between the surface deformation and earthwork. More than 60% of the land subsidence occurred in the filling region, and more than 65% of the land uplift occurred in the mountain excavation region. The large-scale land creation project caused by excavating mountains and filling valleys explains this phenomenon appropriately. On the one hand, over 22 km 2 of flat land was created to relieve the serious pressure caused by rapid population growth [41]. On the other hand, the large amount of earthwork destroyed the natural state and local geological structure of the original loess, producing many filled loess regions and changed the original landscapes and hydrological conditions. Figure 12d shows that the mountain excavation removed the upper part of the mountain and that the remaining part underneath the mountain released a rebound stress due to the process of mass removal. This phenomenon gradually led to surface uplift in the hilly excavation region until a balanced state was reached. In addition, valley filling changed the natural state of the local original loess. In the gully area, the surface of the original loess layer was covered with infilling loess stratum with maximum thickness up to 112 m [37]. The filling loess has the characteristics of weak structure, a relatively loose arrangement of soil particles and poor capacity for bearing external load, resulting in poor engineering features [53]. While the natural development of a loess layer is a time consuming process, including compression, collapsible deformation, consolidation and compaction (the latest loess layer has experienced a long development period of 5000 years) [45], the filling loess stratum in the YND area only developed for three years (the first phase of earthwork in YND commenced in April 2012 and was completed in June 2015). Thus, the structural change and collapsible subsidence of the surface are not completely stable. In addition, the dynamic compaction method used in the consolidation of the foundation in filling area could only weaken the collapsibility and enhance the intensity of filling loess to a certain extent [54]. It cannot completely change the physical characteristics of the collapsible loess, which led the loess foundation to not be fully tamped and more prone to ground subsidence. Moreover, there are serious stability problems in the unnatural connection of loess in filling and excavation areas, which can cause various geological disasters such as loess landslides, loess collapse and land subsidence [38].

Human Activities Contributing to the Land Subsidence
In the aforementioned sections, it was found that land creation and urban construction in Yan'an induced a wide area of surface deformation. Land creation has destroyed the original loess structure of the local area, which has caused severe instabilities in filling loess foundation. In addition, large numbers of buildings have increased the external pressure on surface and accelerated the compressible deformation of the loess layer and building foundations. Figure 12 shows the overlap between the surface deformation area and earthwork region. It is indicated that there is a high correlation between the surface deformation and earthwork. More than 60% of the land subsidence occurred in the filling region, and more than 65% of the land uplift occurred in the mountain excavation region. The large-scale land creation project caused by excavating mountains and filling valleys explains this phenomenon appropriately. On the one hand, over 22 km 2 of flat land was created to relieve the serious pressure caused by rapid population growth [41]. On the other hand, the large amount of earthwork destroyed the natural state and local geological structure of the original loess, producing many filled loess regions and changed the original landscapes and hydrological conditions. Figure 12d shows that the mountain excavation removed the upper part of the mountain and that the remaining part underneath the mountain released a rebound stress due to the process of mass removal. This phenomenon gradually led to surface uplift in the hilly excavation region until a balanced state was reached. In addition, valley filling changed the natural state of the local original loess. In the gully area, the surface of the original loess layer was covered with infilling loess stratum with maximum thickness up to 112 m [37]. The filling loess has the characteristics of weak structure, a relatively loose arrangement of soil particles and poor capacity for bearing external load, resulting in poor engineering features [53]. While the natural development of a loess layer is a time consuming process, including compression, collapsible deformation, consolidation and compaction (the latest loess layer has experienced a long development period of 5000 years) [45], the filling loess stratum in the YND area only developed for three years (the first phase of earthwork in YND commenced in April 2012 and was completed in June 2015). Thus, the structural change and collapsible subsidence of the surface are not completely stable. In addition, the dynamic compaction method used in the consolidation of the foundation in filling area could only weaken the collapsibility and enhance the intensity of filling loess to a certain extent [54]. It cannot completely change the physical characteristics of the collapsible loess, which led the loess foundation to not be fully tamped and more prone to ground subsidence. Moreover, there are serious stability problems in the unnatural connection of loess in filling and excavation areas, which can cause various geological disasters such as loess landslides, loess collapse and land subsidence [38].

Rapid Urban Construction
In the past five years, Yan'an has experienced a rapid urbanization progress. Figure 13 indicates that large numbers of buildings and roads are built to expand the urban space. Since the commencement of earthwork in April 2012, the construction of major roads, buildings and infrastructure in the new district has been completed, concluding in 2018. An investigation of the changing building area in YND from 2013 to 2019 indicates that there is a high correlation between the subsidence phenomenon and building construction, with a correlation coefficient of 0.88 ( Figure  14). Before 2015, few buildings, including Yan'an Municipal Government and Yan'an New District Management Committee were built, and surface deformation was only slightly observed. Between 2015 and 2016, the area experienced a rapid construction period affecting more than 1,300,000 m 2 of the building area from the urban center to the surrounding area. At the same time, severe subsidence with maximum accumulation of 100 mm was detected. From 2017 to 2019, the velocity of urban

Rapid Urban Construction
In the past five years, Yan'an has experienced a rapid urbanization progress. Figure 13 indicates that large numbers of buildings and roads are built to expand the urban space. Since the commencement of earthwork in April 2012, the construction of major roads, buildings and infrastructure in the new district has been completed, concluding in 2018. An investigation of the changing building area in YND from 2013 to 2019 indicates that there is a high correlation between the subsidence phenomenon and building construction, with a correlation coefficient of 0.88 ( Figure 14). Before 2015, few buildings, including Yan'an Municipal Government and Yan'an New District Management Committee were built, and surface deformation was only slightly observed. Between 2015 and 2016, the area experienced a rapid construction period affecting more than 1,300,000 m 2 of the building area from the urban center to the surrounding area. At the same time, severe subsidence with maximum accumulation of 100 mm was detected. From 2017 to 2019, the velocity of urban construction slowed down, and the maximum subsidence rate reduced to 50 mm/year. Generally, when loess is dry, compression deformation under an external load includes a large part of loess deformation. On the other hand, when loess is wet, the loess layer is prone to collapsible deformation under self-weight, and the degree of deformation gradually declines over a relatively long period [54]. Furthermore, large amounts of buildings add great external load to the surface and accelerate the deformation process significantly [13]. In the period of rapid urban construction, the land subsidence is also more visible and comprehensive in YND.
Remote Sens. 2019, 11, x FOR PEER REVIEW 16 of 20 construction slowed down, and the maximum subsidence rate reduced to 50 mm/year. Generally, when loess is dry, compression deformation under an external load includes a large part of loess deformation. On the other hand, when loess is wet, the loess layer is prone to collapsible deformation under self-weight, and the degree of deformation gradually declines over a relatively long period [54]. Furthermore, large amounts of buildings add great external load to the surface and accelerate the deformation process significantly [13]. In the period of rapid urban construction, the land subsidence is also more visible and comprehensive in YND.   construction slowed down, and the maximum subsidence rate reduced to 50 mm/year. Generally, when loess is dry, compression deformation under an external load includes a large part of loess deformation. On the other hand, when loess is wet, the loess layer is prone to collapsible deformation under self-weight, and the degree of deformation gradually declines over a relatively long period [54]. Furthermore, large amounts of buildings add great external load to the surface and accelerate the deformation process significantly [13]. In the period of rapid urban construction, the land subsidence is also more visible and comprehensive in YND.

Natural Factors Contributing to the Land Subsidence
A related study has shown that the loess layer can result in a collapsible deformation when encountering external water, which might cause building instability and ground subsidence [54]. From the field research on humidification, consolidation and compression of loess in Yan'an, the strong collapsibility of loess is the basis of local widespread land deformation. The reasons are described as the following: (1) The filling loess in Yan'an consists mainly of fine soil-while the clay content is about 15%, the silt content is about 76%, and the content of fine sand is about 9%. Moreover, the wet and dry densities are between 1.39 to 1.60 g/cm 3 and 1.24 to 1.38 g/cm 3 , respectively. The low density and content of filing loess has led to a poor structural intensity of soil particles [55].
(2) The loess is loose and porous, and the porosity is as high as 40-50%. In addition to the intergranular pores, a variety of unique micropores, including joints, insect pores, radioactive pores, plant root pores and so on have also developed, leading to a strong water absorption ability in the local loess [54].
(3) The compactness level of the foundation in Yan'an new district is relatively low. It should be indicated that the compactness of the 1-12 m foundation is 64-72%, and the compactness of the 13-33 m foundation is 72-77% [56].
In addition, most precipitation in Yan'an is in the form of heavy rain, which is one of the inducing factors for widespread land deformation in the local area [42]. Rainwater reduces the cohesive intensity and matrix suction of loess clay, causing the loess to soften and an increase in collapsibility. Extreme precipitation (such as a rainstorm) produces a large amount of surface water in a short period of time that exceeds the water storage capacity of the soil. Surface water infiltrating into the loess foundation destroys the internal structure of the loess foundation, which induces the collapse of roadbed and building foundations and leads to the deformation of the surface.

Conclusions
In this study, the SBAS-InSAR method combined with Sentinel-1A data was utilized to calculate the urban land subsidence during the rapid construction period of Yan'an city from 2015 to 2019. The spatial distribution features and evolution of surface deformation in YND were analyzed. Moreover, the reasons for surface deformation and its correlation with urban construction were discussed. The conclusions can written as follows: • The SBAS-InSAR method and Sentinel-1A data achieved favorable results in the application of urban surface deformation monitoring in the Yan'an area, demonstrating the applicability of multi-temporal InSAR technology in the field of high-precision and large-scale land subsidence monitoring.

•
Uneven surface deformation was observed in YND and YNA areas from 2015 to 2019. The average deformation rate and maximum cumulative subsidence in YND were between −70 and 30 mm/year and 300 mm, respectively. Moreover, the area experienced a rapid subsidence period from 2015 to 2017 and a slow subsidence period from 2017 to 2019, with maximum subsidence rates of 100 and 50 mm/year, respectively. The average rate of surface deformation in YNA is between −50 and 25 mm/year, and the maximum cumulative subsidence is 217 mm. The strong collapsibility of loess and concentrated precipitation in local areas are the main natural factors contributing to surface deformation.

•
There was a high correlation between the deformation distribution and the earthwork. In other words, the land creation project carried out by mountain excavation and filling valleys and rapid urban construction are the primary human factors causing surface deformation.