Subsidence Monitoring of Fill Area in Yan’an New District Based on Sentinel-1A Time Series Imagery

Abstract

In recent years, many cities in the Chinese loess plateau (especially in Shanxi province) have encountered ground subsidence problems due to the construction of underground projects and the exploitation of underground resources. With the completion of the world’s largest geotechnical project, called “mountain excavation and city construction,” in a collapsible loess area, the Yan’an city also appeared to have uneven ground subsidence. To obtain the spatial distribution characteristics and the time-series evolution trend of the subsidence, we selected Yan’an New District (YAND) as the specific study area and presented an improved time-series InSAR (TS-InSAR) method for experimental research. Based on 89 Sentinel-1A images collected between December 2017 to December 2020, we conducted comprehensive research and analysis on the spatial and temporal evolution of surface subsidence in YAND. The monitoring results showed that the YAND is relatively stable in general, with deformation rates mainly in the range of −10 to 10 mm/yr. However, three significant subsidence funnels existed in the fill area, with a maximum subsidence rate of 100 mm/yr. From 2017 to 2020, the subsidence funnels enlarged, and their subsidence rates accelerated. Further analysis proved that the main factors induced the severe ground subsidence in the study area, including the compressibility and collapsibility of loess, rapid urban construction, geological environment change, traffic circulation load, and dynamic change of groundwater. The experimental results indicated that the improved TS-InSAR method is adaptive to monitoring uneven subsidence of deep loess area. Moreover, related data and information would provide reference to the large-scale ground deformation monitoring and in similar loess areas.

1. Introduction

Yan’an city of Shaanxi province is located between loess gullies, which greatly restricts the development of the local tourism. In recent years, to solve the problem of land resource scarcity, Yan’an New District (YAND) carried out China’s largest geotechnical project, called “mountain excavation and city construction,” in the area of collapsible loess plateau [1]. The large-scale excavation and filling works have a serious impact on the geological environment of loess, causing uneven ground settlement, slope instability, and other typical geological disasters [2]. Therefore, monitoring the distribution of ground subsidence in the new area and the influence of engineering construction on ground subsidence provides data support for long-term ground monitoring, planning and implementation of construction projects, and scientific prevention and control of ground subsidence in the area, which is of great significance for the sustainable development of urbanization.

Synthetic aperture radar interferometry based on two-view or multi-view synthetic aperture radar image data has been widely used in regional ground subsidence detection, early identification and monitoring of landslides and other geological hazards [3], and has gradually played an indispensable and important role [4] because of its features of all-day, all-weather, wide coverage and high accuracy of deformation monitoring [5]. The proposed small baseline subsets InSAR (SBAS-InSAR) approach [6,7,8,9] has greatly avoided the decoherence phenomenon caused by long temporal and spatial baselines and reduced the influence of topography on the differential and atmospheric effects. Such a method greatly improves the monitoring accuracy of the InSAR technique and thus provides an effective means to detect the spatial distribution and spatial and temporal evolution of regional ground subsidence. However, the traditional SBAS-InSAR approach has a disadvantage in adaptability in the natural surface and vegetation areas [10].

Nowadays, InSAR-based urban ground deformation monitoring is mainly focused on the regional ground deformation under the action of natural and artificial factors such as tectonic activities and over-abstraction of groundwater [11,12,13,14,15], while the research on ground deformation caused by large-scale geotechnical engineering and subsequent urban construction is relatively rare. The mountain excavation and city construction project of YAND mainly consists of excavation and filling as well as construction works, and the current research on its ground subsidence mechanism is mainly based on field level measurement data and borehole data [16], making it difficult to realize a comprehensive analysis of the deformation characteristics of the whole YAND.

In this paper, a total of 89 scenes of ascending-orbit Sentinel-1A image data were acquired from December 2017 to December 2020. Subsequently, the spatial distribution of surface deformation and its spatial and temporal evolution characteristics of YAND were extracted by improved TS-InSAR technology. Then, the main factors causing ground subsidence in the YAND were comprehensively analyzed, which provided ideas for solving the ground subsidence and slope instability problems caused by the construction of loess high-fill geotechnical projects. Related results may provide a reference for scientific site selection for urban expansion.

2. Methodology

The conventional Differential InSAR (DInSAR) [17] technique can detect large-scale deformation, but it suffers from severe atmospheric and out-of-coherence effects beyond a certain spatial and temporal baseline and cannot obtain information on the temporal evolution of deformation. The TS-InSAR [6,7,8,9,18,19] technique (i.e., SBAS-InSAR [6], PS-InSAR [18,19]) can make full use of SAR image data to monitor small-scale deformation without the influence of spatial and temporal decoherence and obtain the time-series information of ground deformation. However, the conventional SBAS-InSAR technique with only spatio-temporal baseline thresholds may also involve interferometric pairs with poor coherence and low quality in the solution, thus affecting the quality of the solution results. We optimize the number of observations involved in the TS-InSAR solution by calculating the long time-series interferometric coherence coefficients and extracting the interferometric pairs with higher coherence from them. The accuracy and reliability of the solution results are improved. In this paper, GAMMA software is used to preprocess the Sentinel-1A dataset, including the selection of the master image, image registration, combination of interferometric pairs and interferometric coherence calculation. The selection of highly coherent interference pairs is performed based on the Python language, and finally the deformation rates are calculated by inversion of the interferometric network using least squares (WLS) solutions in Mintpy software. The overall data processing flow is shown in Figure 1.

The main steps of the improved TS-InSAR data processing include interferometric coefficient setting, interferogram networking, network inversion, and coherent target (CT) selection.

DInSAR ground deformation detection is used to obtain spatial geometric information by extracting the phase difference between two radar waves, but the actual acquired phase signal contains real ground phase information and a large amount of noise signal. Therefore, the spatial coherence (SC) coefficient is an important parameter to evaluate the quality of the differential interferogram [20]. The expression model of the SC coefficient can be expressed by the following equation

$${\gamma }_{spt}=\frac{{\displaystyle {\sum }_{n=1}^N{\displaystyle {\sum }_{m=1}^M{u}_1(n,m)u_2{}^{*}(n,m)}}}{\sqrt{{\displaystyle {\sum }_{n=1}^N{\displaystyle {\sum }_{m=1}^M{\left|u_1(n,m)\right|}^2{\displaystyle {\sum }_{n=1}^N{\displaystyle {\sum }_{m=1}^M{\left|u_2(n,m)\right|}^2}}}}}}$$

(1)

The overall coherence of the interferogram pair will be effectively improved with a limited time baseline. A small baseline interferometric network with high coherence can then be formed. In the background of this principle, we use the SBAS-InSAR method to select CTs with high coherence. Finally, we apply the weighted least squares (WLS) estimation method to obtain the high-precision time-series deformation results of CTs. The combined approach of free networks adds redundant observations and thus improves the ability to be applied in low coherence regions. However, the interferogram network with a small baseline has a greater capability for fast decoherence in the time domain.

The traditional SBAS-InSAR technique forms differential interferogram pairs by selecting appropriate interferogram combinations of spatio-temporal baselines [21] and uses a mathematical model of phase observation to perform phase analysis of highly coherent target points. This combined method increases the temporal sampling rate and effectively overcomes the spatial uncorrelation problem caused by long spatial baselines [21] and extracts the atmospheric delay effects as well as topography and noise while removing the time-series deformation information of the ground surface and atmospheric delay effects as well as topography and noise [6]. The basic principle is as follows.

The acquired (N + 1) SAR images covering the same area are arranged in a certain time order (t₁,t₂, …,t_n), and all images are registered together to a master image, and M differential interferograms are generated by setting a suitable combination of temporal and spatial baseline thresholds, with M satisfying the following conditions.

$$\frac{N+1}{2}\le M\le \frac{N(N+1)}{2}$$

(2)

Suppose the interferogram j is generated by the differential interferential of two SAR images imaged at time $t_A$ and $t_B$ ($t_A$ < $t_B$), then the interference phase at pixel point (l, r) in the interferogram j can be expressed as [22]:

$$\delta {\phi }_j(l,r)=\phi (t_B,l,r)-\phi (t_A,l,r)$$

(3)

The above equation can also be expressed as:

$$\delta {\phi }_j(l,r)\approx {\phi }_{def,j}(l,r)+{\phi }_{topo,j}(l,r)+{\phi }_{atm,j}(l,r)+{\phi }_{noise,j}(l,r)$$

(4)

where $\phi (t_B,l,r)$ and $\phi (t_A,l,r)$ are the phase values at pixel (l,r) at moments $t_B$ and $t_A$, respectively, ${\phi }_{def,j}(l,r)$ is the deformation phase during the two images, ${\phi }_{topo,j}(l,r)$ is the topographic phase, ${\phi }_{atm,j}(l,r)$ is the phase delay error due to atmospheric effects, and ${\phi }_{noise,j}(l,r)$ is the phase information generated by random noise.

Each component in Equation (3) can be expressed as:

$${\phi }_{def,j}(l,r)=\frac{4\pi }{\lambda }\left[d(t_B,l,r)-d(t_A,l,r)\right]$$

(5)

$${\phi }_{topo,j}(l,r)=\frac{4\pi }{\lambda }\cdot \frac{B_{\perp }\Delta \mathrm{h}}{r\sin \theta }$$

(6)

$${\phi }_{atm,j}(l,r)={\phi }_{atm}(t_B,l,r)-{\phi }_{atm}(t_A,l,r)\infty$$

(7)

where $\lambda$ is the radar wavelength, $d(t_A,l,r)$ and $d(t_B,l,r)$ represent the LOS directional cumulative deformation information at times $t_A$ and $t_B$ with respect to the initial time, $B_{\perp }$ represents the vertical baseline value, $\Delta \mathrm{h}$ is the DEM elevation difference, r is the radar line-of-sight distance to the ground observer, and $\theta$ represents the radar incidence angle. After removing various interfering phases, the phase information of each interferogram can be represented by the following model:

$$\Delta \phi =A\phi +\Delta {\phi }_{\epsilon }$$

(8)

where $\Delta \phi ={\left[\Delta {\phi }^1,\dots ,\Delta {\phi }^M\right]}^T$ represents the interferometric phase of each interferogram; $\phi ={\left[{\phi }^2,\dots ,{\phi }^{N+1}\right]}^T$ is the temporal interferometric phase of the other images relative to the reference image (assuming the phase of the reference image is ${\phi }_1$), $\Delta {\phi }_{\epsilon }={\left[\Delta {\phi }_{\epsilon }{}^1,\dots ,\Delta {\phi }_{\epsilon }{}^M\right]}^T$ is the residual interferometric phase error, and A is an M × N matrix representing the combined interferometric patterns, consisting of 1, 0, and −1, where −1 represents the master image and 1 represents the slave image. The best estimate of the interferometric phase of the time series is calculated using the least square norm method [23].

$$\widehat{\phi }=\arg \min \Vert W^{1/2}(\Delta \phi -A\phi {)\Vert }_2={\left(A^{T}WA\right)}^{-1}{A}^{T}W\Delta \phi$$

(9)

where, $\widehat{\phi }$ is the best estimate of the time series interferometric phase, W is the Fisher In-formation Matrix (FIM) [24] diagonal weight matrix used in this paper. Equation (10) represents the model for calculating the Temporal Coherence (TC), which can be used to select CTs with high coherence. Subsequently, by estimating the WLS for the small baseline network, high accuracy time series deformation results for CTs in the study area can be obtained.

$${\gamma }_{\mathrm{tmp}}=\frac{1}{M}\left|{\displaystyle \sum _{i=1}^M\exp \left[\mathrm{k}\left(\Delta {\phi }_i-A{\widehat{\phi }}_i\right)\right]}\right|$$

(10)

where k is the imaginary unit.

3. Study Areas and Data Sources

3.1. Study Areas

Yan’an City in Shaanxi Province is located between 35°21′ and 37°31′N latitude and 107°41′ and 110°31′E longitude. It is a typical linear city located in the hilly and gully area of Loess Plateau, and its old city is densely distributed in the river valley in the shape of “Y”. Three river valleys spread out in a ray from the central point. The YAND is characterized by interlocking valleys and fragmented terrain. This means that the surface of the ground is full of rivers, and the erosion of rivers causes the surface of the ground to form thousands of ravines, and the complete surface of the ground becomes broken. The whole area is divided into the northern and southern valleys and the southeastern remnants of the plateau. The vegetation in the area is sparse and soil erosion is serious. The study area has a typical warm temperate continental monsoon climate. The climate is characterized by dry spring with little rain, hot and rainy summer, rapid cooling in autumn, and sparse rainfall in winter. The average annual water volume is 542.6 mm, and the dry period in recent years is often around 400 mm [25]. The average annual temperature is 9.4 °C.

A large number of old revolutionary sites are preserved in the main urban area of the narrow city of Yan’an. This has led to the problem that the preservation of historical and cultural relics and urban development are mutually constrained. To solve the shortage of land resources that restricts the economic and social development of Yan’an, the “mountain excavation and city construction” project was implemented to build a new urban area by “cutting mountains, filling ditches, creating land, and building cities” to organize the loess gully areas into urban construction lands. However, the speed and intensity of the construction of the YAND are far beyond the geological capacity, which has a significant impact on the geological environment of the loess area. Its engineering geology and hydrogeological conditions were significantly changed, which inevitably caused a large amount of post-engineering ground subsidence phenomenon. An overview of the study area is shown in Figure 2.

The construction of Yan’an New District North (Phase I) started on 17 April 2012, and the geotechnical works were basically completed in 2018. During the period, the rapid construction of building works started in 2014, and the new district carried out the relocation and resettlement work in 2017. With the accelerated urbanization, the type of land use in the new district is also changing rapidly.

3.2. SAR Datasets

In this paper, the SAR dataset acquired by the Sentinel-1A satellite, launched by the European Space Agency (ESA) on 3 April 2014, is selected for the spatial distribution of ground subsidence detection. Sentinel-1A is a C-band synthetic aperture radar imager with a revisit period of 12 d. It forms a dual-satellite system with the Sentinel-1B satellite with a revisit period of 6 d, which was launched on 26 April 2016, and enables all-weather, all-day imaging on a global scale.

A total of 89 scenes of Sentinel-1A radar images were selected in the time span of December 2017–December 2020, with the imaging mode of interferometric wide-area (IW) mode, polarization mode of VV + VH polarization, and spatial resolution of 5 m × 20 m (azimuth × distance). In addition, the digital elevation model of Shuttle Radar Topographic Mapping Mission (SRTM) with 30 m resolution was introduced to remove the topographic undulation error in differential interferometry. The specific data are shown in

Table 1. SAR datasets.

Sensor	Path No.	Number of Images	Date Range
Sentinel-1A	84	32	6 December 2017–25 December 2018
		28	25 December 2018–20 December 2019
		29	20 December 2019–26 December 2020

.

4. Results

4.1. Interferometric Coherence Analysis

As shown in Figure 3, the different colors of the connecting lines represent the average spatial coherence of the interferogram. The overall coherence of the interferogram is improved after refining the constraints of the spatio-temporal baseline. This provides assurance for the reliability and accuracy of the deformation results extracted in the collapsible loess area.

The traditional SBAS-InSAR technique has low interferogram quality in natural surface and vegetation areas. To optimize the interferometric combination for better interferometric phase quality, we calculated the average interferometric coherence coefficients of 2279 interferograms from 89 Sentinel-1A images between December 2017 and December 2020 to obtain Figure 4. As shown in Figure 4, the vertical axis represents the time span of the primary images, and the horizontal axis represents the time span of the secondary images. The temporal baseline threshold is set to 400 days and the perpendicular baseline is set to 200 m. The interferograms from May to November each year have high coherence, all above 0.7. The interferograms for other time periods generally have low coherence, as shown in Figure 4a. This causes an inaccurate result and accumulated errors propagation in the subsequent calculation procedures. Therefore, in terms of coherence, we set a time baseline of 49 days to obtain 330 interferograms with higher coherence coefficients, as shown in Figure 4b. However, this also leads to the loss of some of the high-coherence interferometric pairs. For making full use of high-coherence interferometric pairs, we added all high-quality interferograms into the network as beneficial redundant observations. As shown in the right subplot of Figure 4c, the interferograms with average spatial coherence higher than 0.6 are refined as the available interferograms. The number of interferograms are increased to 476, which means more redundant observations participate in the calculation of the network inversion. More valuably, the accumulated errors and the random residues can be mitigated to some extent.

4.2. Spatial-Temporal Characteristics of Surface Deformation in YAND Area

Figure 5 shows the subsidence rate map of YAND obtained using the TS-InSAR technique. The comparison between the improved TS-INSAR technique and PS-InSAR results shows that the subsidence areas detected by both methods remain basically the same, and the deformation rates are roughly distributed in the interval of −40~30 mm/yr. The accuracy and reliability of the improved TS-InSAR technique proposed in this paper are cross-validated. Figure 5a,b shows that from December 2017 to December 2020, most areas in the YAND remain relatively stable, with deformation rates mainly concentrated in the range of −10 to 15 mm/yr. However, there are three areas with significant ground deformation, namely the SHQ, YZD, and JSQ areas in the figure. All three subsidence areas are distributed in a strip, overlapping with the planned fill area of the YAND, and the maximum subsidence rate reaches −50 mm/yr. The SHQ area starts from North Xingzhi Road, passes through the whole section of South Xingzhi Road and West Shanghai Road to the whole section of East Shanghai Road, and extends to the section of Yan’an Wuyue Plaza and Zichang Road. Numerous construction works exist in this area. YZD is the area along Yanzhou Avenue. JSQ is the YAND south of the Northwest Comprehensive Survey and Design Institute and the Yan’an City Traffic Police Detachment. Land creation is currently underway. The time-series deformation spatial distribution map also shows that the cumulative surface deformation range of YAND decreases and decelerates from 2017 to 2020. This is due to the basic completion of the Pingshan land creation project in 2017 and the gradual completion of the subsidence return.

Figure 6a–c shows the surface deformation rate in the LOS direction for each year from December 2017 to December 2020. The results show that the surface subsidence in YAND is mainly distributed in two areas, SHQ and YZD, during the period from December 2017 to December 2018. The overall subsidence in SHQ area is banded and has further development with the road direction. The average subsidence rate in this area exceeded −50 mm/yr, and the maximum subsidence rate detected was −62 mm/yr. The YZD area was developing with an average rate of −80 mm/yr along Yanzhou Avenue. The maximum deformation rate reached −107 mm/yr and extended outward to the intersection of other highways. Among them, fewer deformation points were detected in the southern section of Yanzhou Avenue, and it was inferred from the historical optical images that the deformation caused by the frequent engineering disturbances and filling and excavation were beyond the detection range of InSAR technology and decoherence phenomenon occurred because the section was in the critical period of construction. The interferogram coherence was better and a new subsidence funnel was detected from 2018–2019, than 2017–2018. This was also in the JSQ region, where the subsidence rate was most significant, reaching an order of magnitude of −93 mm/yr. The 2019–2020 detections show that the original subsidence funnel maintained a subsidence rate of −20~15 mm/yr.

4.3. Analysis of Deformation Characteristics of Typical Regions

There are severe inhomogeneous deformations in the built-up area of the new district (SHQ), along Yanzhou Avenue (YZD), and in the eastern land-making area (JSQ). Local differential deformations can produce tensile stresses on the ground. These tensile stresses may lead to secondary disasters such as unstable building foundations, ground cracking or even collapse, and pipeline rupture [22]. In this paper, three typical uneven subsidence areas detected are analyzed based on time series InSAR results.

4.3.1. SHQ Region

Figure 7a shows that from 2017 to 2020, the average deformation rate in the SHQ area is −25 to 15 mm/yr. The deformation area underwent large-scale construction development. The most serious area of subsidence is located in the section of Yan’an Wuyue Plaza and Zichang Road Wuyue Plaza (point P1), where the deformation rate reaches −21.1 mm/yr. In addition, according to Figure 7b, point P1 shows a general sinking trend and a gradual increase in the subsidence rate. During the period, it kept a stable state from December 2017 to September 2018. Figure 7c–e shows that the Wuyue Plaza is in a state of planning and ground leveling, and large-scale building construction has not yet been carried out. From September 2018 to June 2019, the ground subsidence gradually increased and accumulated to 30 mm. This is because the Wuyue Plaza is under full-scale construction and the loess structure in the filling area is weak and the loess pores are large, which makes the contact area of loess particles small and the ability to resist external load becomes poor. Thus, the loess structure rapidly disintegrates, leading to collapse. Two obvious rebounds occurred in June, August, and September 2019, which were around 10 mm, probably due to the backfill of road reconstruction. From September 2019 to December 2020, the area continued to sink. At this time, the main body of the building and the road were completed, and the ground subsidence was still mainly related to the collapsibility of the loess. In addition, this deformation area overlaps with the filling area. Although the progress of the YAND city-building project slowed down in 2019 compared with the previous one, the rapid economic development in recent years led to the rapid rise of a large number of high-rise buildings. The weight of many buildings is loaded on top of the strata, leading to the growth of static load. The larger and denser the building, the more obvious the effect on an engineered ground subsidence [26]. Therefore, the main causes of severe local deformation are human factors such as land development and urban construction.

4.3.2. YZD Region

According to Figure 8a, the YZD area has experienced more serious subsidence, with an annual average deformation rate of −37~15 mm/yr. According to the cumulative subsidence at point P2 in Figure 8b, there is a general trend of gradual slowdown in the subsidence rate in the area. From December 2017 to March 2018, the subsidence fluctuates, which is related to the frequent dredging and filling works at the beginning. After March 2018, the overall trend decreased, and the accumulated subsidence reached 120 mm. Figure 8c–e shows the YZD area belongs to the hilly mountainous area in 2016, and in December 2017, the land creation was completed through the “mountain excavation and city construction” project. At the same time, the construction of Yanzhou Avenue was started through the government plan, and in March 2020, Yanzhou Avenue became the main traffic route for the construction of the outer part of the YAND. The area was almost covered with fill during the process of making a city on a flat hill. Among the fill works involved in the “mountain excavation and city construction” project, the maximum thickness of the fill exceeded 100 m. The distribution of the fill and its thickness obviously affected the distribution and size of ground subsidence. Previous studies monitored the ground subsidence of high fills in YAND in layers. The total subsidence of the fill body consists of both the subsidence of the fill body itself and the subsidence of the original foundation, and the subsidence of the fill body itself is the main part [27]. From

Table 2. Fill thickness and settlement of some measurement points in YAND [27].

NO.	Filling Depth/m	Thickness of Original Foundation/m	Total Subsidence/m	Subsidence of Original Foundation/m	Filling Subsidence/m	Subsidence Ratio of Original Foundation/%
1	31.3	25.1	0.672	0.409	0.263	60.8
2	47.7	16.5	1.182	0.227	0.955	19.2
3	50.2	7.8	1.260	0.354	0.906	28.1
4	45.5	15.9	1.272	0.471	0.801	37.0
5	21.4	39.5	0.249	0.010	0.149	40.1
Ave.	39.2	21.0	0.927	0.312	0.615	37.1

, the original foundation subsidence and the subsidence of the fill body itself account for 37% and 63% of the total subsidence, respectively, on average. Obviously, the average thickness of 39.2 m of the fill body’s own subsidence plays a dominant role in the total subsidence of the high fill during the construction period. Since the original foundation loess has certain structural properties, it produces less compression under the overlying load, while the fill body is a remodeled loess, and the bond strength between soil particles has been destroyed, which can still produce further compression under the long-term load. Therefore, the presence of large amounts of fill and high traffic density, dynamic and static loads associated with large trucks, and construction may be the main reasons for subsidence in this area [28].

4.3.3. JSQ Region

In addition to the area where the two points P1 and P2 are located, a subsidence funnel was also detected on the east side of the eastern traffic police detachment of the YAND (Figure 9a), and the range of the subsidence rate was mainly concentrated in the range of −50 to 15 mm/yr. Figure 9c–e shows that from March 2017 to October 2018, a large area of mountainous land was shifted and transformed into construction land. As a result, the land cover type changed, leading to the decoherence of the TS-InSAR detections. The change of land cover type due to the hill cutting is the main reason for the surface deformation. Figure 9b shows the trend of cumulative subsidence at point P3 in the area from December 2018 to December 2019, which shows that the cumulative subsidence at this point reached 100 mm in one year. The shear strength of the soil in the construction area changes when the dynamic load is greater than a certain vibration acceleration value. Under cyclic loading, the strain value of the soil increases with the increase of dynamic stress and increases with the increase of the number of load cycles. The sunken piles as well as the vehicle movements generate a fixed frequency of vibrations, and this regular dynamic load greatly affects the structure of large-thickness fills, roadbeds, and collapsible loess. This area is still in the construction phase and is more significantly affected than the living area where construction was completed. Therefore, the vibration from the impact of the sunken piles during construction combined with the dynamic load from the traffic cycles of large construction vehicles was the main factor causing ground subsidence.

5. Discussion

5.1. Influence of the Nature of the Fill Soil

Previous studies indicated that collapsibility is the typical characteristic of loess [29,30]. The typical hydrogeological profile before and after the project is shown in Figure 10, where GC indicates the location of the borehole [31]. The stratigraphy of the study area is mainly composed of Malan loess and Lishi loess and Jurassic sand mudstone. The fill is mainly composed of Malan loess and Lishi loess.

Large-scale land development in YAND and urban construction between 2017 and 2020 led to significant ground subsidence in the areas covered by these soft sediments. The change of land cover type destroys the original loess structure, which has a unique structure, and its soil particles will form high-strength cementation due to physicochemical action under long time contact. In cases where the original loess structure is destroyed, the cementation strength and connection strength between the particles will be significantly weakened, thus reducing the structural strength of the loess and causing hazards such as foundation instability. Although remodeled compacted loess does not have the unique large pore structure of original loess and does not have the primary conditions to produce collapsibility deformation, many experimental studies have found that compacted loess still has specific collapsibility [32]. There are even studies showing that under the same moisture content and dry density conditions, the repeated compaction of the soil layer during the filling process does not eliminate the collapsibility of the loess, and the collapsibility and compressibility of the compacted loess are higher than that of the original loess [33]. According to the literature [34], indoor experiments were conducted on in origin loess and fill loess samples in YAND. The soil samples at different locations in the YAND have different basic physical properties (see

Table 3. Differences in physical properties of original and fill loess [34].

Sampling Loess	Natural Density ρ/(g cm−3)	Moisture Content ω/%	Dry Density ρ/(g cm−3)	Radio Gs	Void Ratio e	Plastic Limit ωp/%	Liquid Limit ωL/%	Plasticity Index Ιp/%
Original loess	1.78	9.88	1.62	2.70	0.67	15.88	31.8	15.92
Filled loess 1	1.80	16.66	1.54	2.58	0.67	19.44	29.33	9.89
Filled loess 2	1.81	12.65	1.61	2.65	0.64	17.52	30.15	12.63
Filled loess 3	2.05	18.02	1.74	2.64	0.52	19.40	33.66	14.26

), which indicates that the physical properties of the loess were changed by remodeling during the filling project and that the filling project did not guarantee uniform physical properties of the fill. Therefore, there are differences in the mechanical properties of the fill at different locations, which affect the differential subsidence of the fill at different locations.

In addition, the physical properties and thickness of the fill soil are important influenced factors for the ground subsidence. With the increase of the fill soil, the pressure of the original loess below is more significant, and the ground subsidence has a relatively faster rate [33]. Pu [34] et al. conducted a series of experiments on original loess and fill loess samples in YAND, and the experimental results showed that the physical properties of the remodeled loess are not homogeneous. Compared with those of the original loess, the physical properties of the remodeled loess are different. There are differences in the mechanical properties of the fill at different locations, which induce the uneven ground subsidence of the fill at different locations [34]. The “mountain excavation and city construction” project in YAND destroyed the original natural state and local geological structure of loess. The ditches and valleys are filled with many loess areas, changing the original landscape hydrogeological conditions, leading to significant uneven ground subsidence.

5.2. Impact of Rapid Urban Construction on Surface Deformation

In the last four years, the YAND has accelerated its urbanization process. The earthworks in the YAND started in 2012, and with the completion of the filling works, the construction works also started in 2014, and the accelerated urbanization started in 2017, with construction loads acting on the soil, making the remodeled loess further compressed and thus accelerating the ground subsidence process [23]. In general, the subsidence area and the traffic road network show a strong correlation. The subsidence area is located in the northwest-southeast direction, consistent with the spatial distribution of the main traffic roads. Previous studies conducted buffer zone analysis of the road network in the subsidence area of the new district [34], and the experimental results also showed that the ground subsidence is mainly distributed within 100 m from the road, and its subsidence observation points account for 83.9% of all the observation points with deformation rate <−10 mm/a. Due to the relatively dense distribution of the road network in YAND, observation points with deformation rates in the range of −10 to −73 mm/a can be detected within 200 m from the road. The number of subsidence observation points keeps increasing as the distance from the road decreases, especially the observation points with higher subsidence rate (<−30 mm/yr). This indicates that the subsidence centered on the road is more severe. The long-term cyclic traffic load brought by the dense road network and its generated vibration load will accelerate the consolidation of its lower and peripheral fill, impacting the surface subsidence.

5.3. Groundwater Dynamic Change Impact

The dynamics of groundwater in the YAND are also one crucial influence factor of the ground subsidence. Previous studies showed that the industrial and urban living exploitation areas are mainly affected by anthropogenic influence. Excessive groundwater extraction over a long period causes the water level to drop [35] and then causes loess and foundation subsidence when the water level rises due to abundant water periods or underground pipeline blockage. Among the kinds of groundwater, there are several types (i.e., quaternary diving, bedrock pore water, bedrock fracture water) that significantly influence ground subsidence [36], especially to the loess hilly gully area. The groundwater types are quaternary pore diving and Jurassic bedrock fracture water. Quaternary pore diving is mainly distributed in the valley area, and bedrock fracture water is distributed throughout the area; the two are closely linked in the valley area and constitute a double-layer medium unified water-bearing body. After the filling and excavation project, the original hydrogeological conditions were destroyed. The compacted backfill loess is integrated with the original loess on both sides of the gully, which changes the original surface water runoff path and blocks the original groundwater drainage channel. Moreover, the permeability coefficient of natural loess is generally 0.02~0.30 m/d, and after compaction, the permeability coefficient is generally 0.001~0.110 m/d, and the permeability coefficient of individual points is very small, reaching 1.1 × 10⁻⁵ m/d [36]. The permeability coefficients of both the original loess and the compacted fill exhibit obvious discrete patterns, and the permeability of the compacted fill is 1 to 2 orders of magnitude lower than that of the original loess because the structure of the fill is destroyed during the remodeling and compaction process [36]. Although drainage measures such as drainage blind ditches were laid out in the fill project, the poor drainage conditions of the fill soil itself made the drainage blind ditches still drain less than the total amount of runoff before the construction of the project [36]. After the construction of the new area began, groundwater level monitoring points were deployed in Qiaoer Valley in the study area, and dynamic groundwater level monitoring work was conducted by hydrogeological boreholes between March 2014 and July 2015 [31]. Figure 11a shows the groundwater level change curve at a monitoring point in the middle reaches of the main ditch of Qiaoer Valley, and Figure 11b shows the groundwater level change curve at a monitoring point in the downstream of the main ditch of Qiaoer Valley. Figure 11a shows that the groundwater level in the middle reaches started to rise at 6 months after filling, rising about 20 cm and finally stabilizing at about 38.6 m. The monitoring point downstream of the main ditch decreases by 10 cm to 64.1 m in the first month of borehole monitoring, and then rises linearly, gradually stabilizing at 62.8 m after the tenth month, with a rise of 1.3 m. The loess is extremely water-sensitive, and the rising trend of the groundwater level has an infiltrative effect on the fill above. The large amount of water retained inside the fill caused the groundwater level in the fill area to rise continuously and cause further subsidence of the fill soil [37,38,39].

6. Conclusions

The geological conditions of the Loess Plateau are relatively fragile, especially in Shaanxi province. The uneven ground subsidence was caused by the huge scale of the mountain cutting project. In this paper, 89 Sentinel-1A radar images covering the study area were processed using the improved TS-InSAR technology. Information on the spatial distribution characteristics of ground subsidence before and after the construction of YAND was obtained, and the main factors of the image ground subsidence distribution were analyzed.

• Three significant subsidence areas were detected within the northern city-making area of YAND, and the deformation rates were concentrated in the range of −50~15 mm/yr, and the subsidence areas were in a northwest-southeast direction.
• There is uneven severe deformation in YAND, and the land subsidence along Wuyue Plaza, Yanzhou Avenue, and Eastern Development Zone is relatively significant. The maximum subsidence rates detected by improved TS-InSAR are −21.1 mm/yr, −37.5 mm/yr, and −50.2 mm/yr, respectively. The typical subsidence funnels appear in the construction area and fill area. From 2017 to 2020, the subsidence funnels enlarged and their subsidence rates accelerated.
• The backfilled compacted loess in YAND is compressible and collapsible, which is the basis for the occurrence of ground subsidence. The remodeling of loess and the change of its physical and mechanical properties in the filling project are the main intrinsic factors of ground subsidence in the filling area. The distribution of ground subsidence is determined by the filling and excavation works in the construction of the “mountain excavation and city construction” project. In addition, the thickness of the fill is the main controlling factor for the distribution and size of ground subsidence. In addition, human activities and changes in geological and environmental conditions also accelerate the development of ground subsidence.