Land Subsidence Monitoring and Building Risk Assessment Using InSAR and Machine Learning in a Loess Plateau City—A Case Study of Lanzhou, China

: As a representative city located in the Loess Plateau region of China, Lanzhou is affected by various environmental and engineering factors, such as precipitation, earthquake subsidence, and building construction, which all lead to frequent geological disasters. Obtaining information on land subsidence over a long time series helps us grasp the patterns of change in various types of ground hazard. In this paper, we present the results of using Interferometric Synthetic Aperture Radar (InSAR) to monitor land subsidence in the main urban area of Lanzhou from 26 October 2014 to 12 December 2021. The main inﬂuential factors leading to subsidence were analyzed and combined via machine learning simulation to assess the land subsidence risk grade distribution of a building unit. The results show that the annual average deformation rate in Lanzhou ranged from − 18.74 to 12.78 mm/yr. Linear subsidence dominated most subsidence areas in Lanzhou during the monitoring period. The subsidence areas were mainly distributed along the Yellow River, the railway, and villages and towns on the edges of urban areas. The main areas where subsidence occurred were the eastern part of Chengguan District, the railway line in Anning District


Materials and Methods
The four main urban areas of Lanzhou, Gansu Province, China, were used as the study area. The vegetation cover in these areas is low, and C-band data could be effectively used to meet our research needs. We selected 150 scenes covering the study area from 26 October 2014 to 12 December 2021, Sentinel-1A Descending Single-Look Complex (SLC) as the dataset (25 December 2022, https://search.asf.alaska.edu); the imaging mode was Interferometric Wide Swath (IW), and the polarization was VV polarization. In addition, topographic phase removal, using 30m resolution SRTM DEM data published by NASA, was used in the experiment; the used SAR image parameters are listed in Table 2. The precision orbit data were removed using POD precision orbiting ephemeris data published by the European Space Agency (ESA). Digitized 1:500,000 Chinese geological maps, Chinese road and river network vector data (10 March 2023, https://kmap.ckcest.cn/), and (10 March 2023, http://data.cma.cn/) annual precipitation raster data for Gansu Province provided by the National Meteorological Science Data Centre were collected during clas-

Materials and Methods
The four main urban areas of Lanzhou, Gansu Province, China, were used as the study area. The vegetation cover in these areas is low, and C-band data could be effectively used to meet our research needs. We selected 150 scenes covering the study area from 26 October 2014 to 12 December 2021, Sentinel-1A Descending Single-Look Complex (SLC) as the dataset (25 December 2022, https://search.asf.alaska.edu); the imaging mode was Interferometric Wide Swath (IW), and the polarization was VV polarization. In addition, topographic phase removal, using 30m resolution SRTM DEM data published by NASA, was used in the experiment; the used SAR image parameters are listed in Table 2. The precision orbit data were removed using POD precision orbiting ephemeris data published by the European Space Agency (ESA). Digitized 1:500,000 Chinese geological maps, Chinese road and river network vector data (10 March 2023, https://kmap.ckcest.cn/), and (10 March 2023, http://data.cma.cn/) annual precipitation raster data for Gansu Province provided by the National Meteorological Science Data Centre were collected during classification of the spatial risk level of ground subsidence in Lanzhou. Relevant data on the groundwater level in Lanzhou were obtained via digitization from the article "Impact from subway project on groundwater environment in Lanzhou fault basin" [66]. This study used the trial version of ENVI 5.6 SARscape for time series deformation analysis.

SBAS-InSAR Technology
In this study, data were processed using SBAS-InSAR, a method proposed by Berardino et al. to fill in the gaps of the PS-InSAR method. While the PS-InSAR method is suitable for non-urban or densely vegetated areas, the SBAS-InSAR method is more suitable for complex terrain. The SBAS-InSAR method is applied for distributed target analysis and is based on the principle of the phase analysis of coherent targets to obtain temporal distortion. The principle of processing flow can be described as follows: N+1 SLC images were sorted in temporal order (t 0 , t 1 , . . . , t N ), the SAR images were aligned and divided into sets according to the spatio-temporal baseline, and differential interference is performed on the images between the sets to obtain M interferometric relative, with M [36] satisfying: In this study, we obtained 596 differential interferometric pairs ( Figure 2). The yellow dot is the super primary image, the green dot is the secondary image, and the red dot is the image outside the constraints. However, since the absence of some time periods in the pre-2016 images previously resulted in low relative coherence, we removed these low-coherence and poorly spread interferometric pairs, yielding final time series results with a significantly larger time series interval before 2016 than after 2016.
The phase composition equation ϕ int,x for the image element x is: 6 of 32 where ϕ de f ,x is the deformation phase, ϕ topo,x is the terrain phase, ϕ f lat,x is the flat earth phase, ϕ noise,x is the system thermal noise, ϕ atm,x is the atmospheric delay phase, ϕ orb,x is the phase due to orbital error, h is the elevation error, B ⊥ is the vertical baseline length, R is the line-of-sight to slope distance, λ is the radar wavelength, and θ is the angle of incidence. The POD data were used to reduce systematic errors due to orbital errors, and to subtract the terrain phase modeled by the DEM to produce a de-flattened differential interferogram, so that the differential interferometric phase ϕ di f ,x of the image element x is: where ϕ di f ,x is the differential interference phase, ϕ ε,x is the topographic residual phase, and ε is the DEM error. The deformation phase ϕ de f ,x consists of the linear deformation phase ϕ linear,x and the non-linear deformation phase ϕ nonlinear,x : where v is the average deformation rate, and T is the monitoring period.
The differential interference phase model ϕ di f ,x can be rewritten as: where ϕ res,x is the residual differential phase, consisting of the system thermal noise, nonlinear deformation, and atmospheric delay phase: The Goldstein filter was used for phase filtering, and Delaunay 3D and MCF were used for phase unwrapping [67]. ϕ di f ,x can be calculated based on manually selected high-coherence control points for the orbit and phase offset to remove orbit residuals and phase ramps; v was solved using the least squares parity method. ϕ res,x was solved using the singular value decomposition method, and then, the singlepoint deformation phase equation was determined, the elevation error solved, and the deformation rate obtained [68]. Finally, the atmospheric phase and nonlinear deformation were retrieved via the residual phase, the deformation time series in the monitoring period was solved, and the vertical surface deformation information of Lanzhou was obtained.

PS-InSAR Technology
PS-InSAR uses statistical methods to select the permanent scatterer (PS) computational points capable of maintaining high coherence and strong scattering over long periods, and to model the time series of the image element phases of all PS points to obtain deformation displacement values and deformation rates. In contrast to the continuous surface results of SBAS-InSAR, this method gives non-continuous points.
Using the acquisition of N+1 SAR images of the same area at the monitoring time, a primary image was selected based on maximizing the relative correlation of interferences; the primary image was aligned with the secondary image, and N interferometric pairs The Goldstein filter was used for phase filtering, and Delaunay 3D and MCF were used for phase unwrapping [67].
, can be calculated based on manually selected high-coherence control points for the orbit and phase offset to remove orbit residuals and phase ramps; v was solved using the least squares parity method.
, was solved using the singular value decomposition method, and then, the single-point deformation phase equation was determined, the elevation error solved, and the deformation rate obtained [68]. Finally, the atmospheric phase and nonlinear deformation were retrieved via the residual phase, the deformation time series in the monitoring period was solved, and the vertical surface deformation information of Lanzhou was obtained.

PS-InSAR Technology
PS-InSAR uses statistical methods to select the permanent scatterer (PS) computational points capable of maintaining high coherence and strong scattering over long periods, and to model the time series of the image element phases of all PS points to obtain

Random Forest Model
In this study, we used a random forest model for spatial hazard classification. The random forest model is a supervised integrated machine learning model developed based on the weak assessor CART decision tree. This decision tree is based on the core idea of constructing multiple mutually independent assessors, and then, averaging or applying the majority voting principles of their predictions to determine the outcome of the integrated assessors [69,70].
Evaluation factors must be screened and quantified to construct a quantitative land subsidence hazard evaluation map. We analyzed the land subsidence influential factors according to the land subsidence monitoring result in Lanzhou and selected eight indicators as evaluation factors: average settlement rate, cumulative settlement, lithology, building height, distance from railway, distance from rivers, distance from faults, and average annual precipitation, in combination with expert opinions [71,72]. The natural break method was used to divide the original continuous evaluation factors into subclasses [73] and construct input variables for the land subsidence hazard evaluation model; afterward, the FR algorithm was used to analyze and calculate the influence of each subclass of condition factors on land subsidence [74,75]. As a bivariate statistical method, the method reflects the correlation between land subsidence and each predicted layer, as shown in Equation (17) [76]. Table 3 lists the values calculated in this study. The interference phase of a particular SAR image pixel in PS-InSAR is: where ϕ de f is the deformation phase, ϕ topo is the residual topographic phase, ϕ f lat is the phase of the reference ellipsoid, ϕ atm is the atmospheric delay phase, and ϕ noise is the noise phase. The terrain phase was simulated using SRTM DEM to remove the flat earth phase and generate a differential interferogram to calculate the amplitude departure index.
where ∆h error represents the DEM error, and ϕ de f can be subdivided into linear and nonlinear deformation phases: where ϕ linear is the linear deformation phase, ϕ non−linear is the non-linear phase, ∆v is the linear deformation rate, and T is the time interval from the main image. A threshold value of 0.75 was set based on the amplitude departure index to obtain 452619 PS candidate points. The phase obtained at the PS point according to the Delaunay triangulation network is expressed as: where ϕ di f is the value of the unwinding phase. Assuming that t |ϕ res | < π, the joint calculation of multiple differential interferograms yields the incremental elevation error and linear deformation rate; the expression is as follows: where γ values are in the range [0, 1] and represent an indicator for judging PS points, N − 1 is the number of interferograms, and ϕ is the phase of PS points on a given interferogram. The mean displacement rates of the PS candidates were estimated with DEM correction factors, the atmospheric phase was removed, and the final mean displacement rates and cumulative surface deformation variables were calculated and geocoded.

Random Forest Model
In this study, we used a random forest model for spatial hazard classification. The random forest model is a supervised integrated machine learning model developed based on the weak assessor CART decision tree. This decision tree is based on the core idea of constructing multiple mutually independent assessors, and then, averaging or applying the majority voting principles of their predictions to determine the outcome of the integrated assessors [69,70].
Evaluation factors must be screened and quantified to construct a quantitative land subsidence hazard evaluation map. We analyzed the land subsidence influential factors according to the land subsidence monitoring result in Lanzhou and selected eight indicators as evaluation factors: average settlement rate, cumulative settlement, lithology, building height, distance from railway, distance from rivers, distance from faults, and average annual precipitation, in combination with expert opinions [71,72]. The natural break method was used to divide the original continuous evaluation factors into subclasses [73] and construct input variables for the land subsidence hazard evaluation model; afterward, the FR algorithm was used to analyze and calculate the influence of each subclass of condition factors on land subsidence [74,75]. As a bivariate statistical method, the FR method reflects the correlation between land subsidence and each predicted layer, as shown in Equation (17) [76]. Table 3 lists the FR values calculated in this study.
where %Subsidence (+) refers to the ratio of the settlement area in each sub-category of evaluation factors, and %domain (+) refers to the ratio of each sub-category of evaluation factors to the total area.  Figure 4 shows the annual average vertical deformation velocity obtained using SBAS-InSAR for the main urban area of Lanzhou from October 2014 to December 2021. Positive values (blue) indicate that the target moves along the ground in the vertical bulge direction. Conversely, negative values (red) indicate that the target moves vertically downward along the ground, representing land subsidence. During the study period, extensive land subsidence occurred in the main urban area of Lanzhou, with an average annual deformation rate ranging from −18.74 to 12.78 mm/yr. According to the technical standard of the China Geological Survey (DD2014-11) (DD2014-11, 2014), it belongs to a lower degree of land subsidence, and therefore, the surface subsidence in the main urban area of Lanzhou is relatively stable overall.  The land subsidence pattern in Lanzhou is uneven; there is a northwest-southeast directional subsidence zone, with no large subsidence funnels but several small spread- The land subsidence pattern in Lanzhou is uneven; there is a northwest-southeast directional subsidence zone, with no large subsidence funnels but several small spread-out subsidence funnels. Subsidence mostly occurs along the Yellow River, railway lines, and villages and towns at the edge of urban areas. The deformation situation can be summarized as one uplifted area and four subsidence areas, namely, subsidence area I along the railway in Anning District, subsidence area II in the south of Xigu District, subsidence area III in the urban area of Qilihe District, subsidence area IV in the eastern part of Chengguan District, and the uplifted area in the western part of Chengguan District. The area share of deformation by the administrative district is shown in Table 4. Among the above four subsidence areas, the area share of the subsidence areas is ranked as follows: III Qilihe District > I Anning District > IV Chengguan District > II Xigu District. The largest and most severe subsidence area occurs in III Qilihe District. The most severe subsidence mainly occurs at the edges of the four urban areas, mainly in the following locations: Baidaooping Village (Figure 4 Figure 3), Liuquan Town (Box K in Figure 3), and along the railway, and its range has been increasing since then. By 14th of May 2018, the lateral expansion of subsidence almost ceased, but the related deformation continued in the subsidence areas. At the same time, an uplift occurred on the western side of IV Chengguan District and the northern part of the II Xigu urban area. We selected cross-sectional points located in the four administrative regions of severe deformation in Figure 6 to study their time series deformation sequences. Figure 7a-d show the time series deformation plots for these cross-sectional points, and Figure 7e-h show the point line plots of the time series change for some subsidence points. Overall, the time series results for the nonsubsidence areas of Lanzhou show smooth seasonal fluctuations, while most time series deformation characteristics of subsidence areas show uniform subsidence. However, there are still several subsidence sites with slightly different time series subsidence patterns. For example, during the monitoring period, the subsidence rate at points in I Anning District slowed down after May 2017, while the subsidence rate at some points in IV Chengguan District increased after March 2017, which might be related to the subsidence mechanism of each point.

Discussion of the Causes of Land Subsidence in Different Regions
Details of the area of significant subsidence in Lanzhou are shown in Figure 4. Factors contributing to subsidence include groundwater extraction, river erosion, lithology, fault distribution, underground pressure, and pressure from high-rise buildings. The significant subsidence in areas S, H, J, K, and G is mainly due to the lithology of loose Quaternary deposits-loess. Due to the large void structure of the loess, the loess undergoes collapse deformation upon exposure to wet loads because the volume of hydraulic consolidation is reduced. This process leads to excessive deformation and instability of the foundation [77,78]. Fortunately, the area is situated at a low building height, mostly between 1 and 18 m. Thus, to prevent excessive loading of the building on the loess, buildings taller than 20 m are unsuitable for the site, as this could aggravate the ground settlement and cause disasters such as foundation collapse. Areas A, B, C, E, T, V, L, and S have railways and subways passing through them, all of which have undergone major or minor deformation. Railways are prone to roadbed deformation due to differences between the power of the track superstructure and the intensive dynamic loads of railway transport or adverse effects, such as climatic and engineering geological factors [79]. The deformation in this area can be caused by the railway construction, especially by the large subsidy of Lanzhou West Railway Station and Lanzhou North Railway Station. Station settlement is more likely to occur due to the dense population and the number of tracks and trains arriving and departing daily. Lanzhou West Station is much larger than Lanzhou North Station, and the subsidence in this area is also more significant, with a maximum deformation of −6.18 mm/yr; so, extra attention has to be paid to monitoring the condition of roadbed deformation in this section and predicting its development in time. In addition, the western part of Chengguan District, where A, B, C, and S are located, has become the second largest area of subsidence in Lanzhou due to the influence of railways, faults, rivers, and buildings, with a maximum settlement rate of −11.68 mm/yr. This area needs special attention from the railway department regarding repair and maintenance, and is not suitable for building buildings that are too tall. D, E, F, W, and I are all located in areas with great building height and density, and all of these areas have experienced significant subsidence or formed subsidence funnels. Commonly, surface construction loads cause compressive land subsidence [54]. The construction of dense building clusters and underground cross-river tunnels, and the fact that the soil layer is mostly soft, can cause an engineering disturbance phenomenon in the soil, and the stress balance of the foundation soil layer can be easily disrupted, leading to vertical deformation. These areas are mostly located in commercial zones or upcoming residential developments in the central city, with a tendency to build high-rise buildings. Therefore, the government should control the construction of floors in these areas; otherwise, they are prone to collapse due to unstable foundations. Areas M, N, and F are located on one side of the fault. Subsidence in this area is characterized by subsidence on one side of the fault and an uplift on the other side, especially in the eastern part of area F where most of the urban area is uplifted; it is assumed that the subsidence in this area is mainly caused by geological and tectonic movement.

Accuracy Evaluation of Land Subsidence Monitoring
This study also verifies the experimental accuracy of Land Subsidence Monitoring by comparing the vertical directional variable results obtained using the SBAS-InSAR and PS-InSAR processing methods, indicating that the deformation point results obtained using these two methods are consistent in most areas. To better compare the results of the two methods, we used all the points of the same name obtained for the sedimentation rate, and generated a histogram, Figure 8. At the same time, considering the coarse differences in the reference data, and based on the trend of the subsidence field reflected in the InSAR monitoring results, we removed some points with significant differences, and selected 7729 subsidence rate points that were evenly distributed throughout the study area for 'point-by-point' comparison of the observations; finally, we constructed a correlation which is shown in Figure 9.

Discussion of the Causes of Land Subsidence in Different Regions
Details of the area of significant subsidence in Lanzhou are shown in Figure 4. Factors contributing to subsidence include groundwater extraction, river erosion, lithology, fault distribution, underground pressure, and pressure from high-rise buildings. The significant subsidence in areas S, H, J, K, and G is mainly due to the lithology of loose Quaternary deposits-loess. Due to the large void structure of the loess, the loess undergoes collapse deformation upon exposure to wet loads because the volume of hydraulic consolidation  A Gaussian fit was used in the histogram in Figure 8. We found that the mutual differences between the two datasets conform to a normal distribution with a mean value of −1.64 mm/yr. However, PS-InSAR generally gives slightly lower values, probably due to the different spatial resolutions of the two. The R-Squ values represent the goodness of fit, and the Adj. R-Squ values represent the adjusted goodness of fit, measuring the degree of merit of the Gaussian regression. The closer the value to 1, the better the fit result. The final fitted equation of the normal function obtained is:  A Gaussian fit was used in the histogram in Figure 8. We found that the mutual differences between the two datasets conform to a normal distribution with a mean value of −1.64 mm/yr. However, PS-InSAR generally gives slightly lower values, probably due to the different spatial resolutions of the two. The R-Squ values represent the goodness of fit, and the Adj. R-Squ values represent the adjusted goodness of fit, measuring the degree of merit of the Gaussian regression. The closer the value to 1, the better the fit result. The final fitted equation of the normal function obtained is: A Gaussian fit was used in the histogram in Figure 8. We found that the mutual differences between the two datasets conform to a normal distribution with a mean value of −1.64 mm/yr. However, PS-InSAR generally gives slightly lower values, probably due to the different spatial resolutions of the two. The R-Squ values represent the goodness of fit, and the Adj. R-Squ values represent the adjusted goodness of fit, measuring the degree of merit of the Gaussian regression. The closer the value to 1, the better the fit result. The final fitted equation of the normal function obtained is: The parameters in Equation (18) are y 0 = 128 ± 37, x c = −1.64 ± 0.04, w = 2.16 ± 0.08, A = 8217 ± 314. Figure 9 shows the linear correlation regression analysis of the two datasets. The linear regression equation y = a + bx is obtained, where a = 0.63 ± 0.02 and b = 1.04 ± 0.01, with a correlation coefficient ρ = 0.9343. The closer the correlation coefficient to 1, the higher the correlation between PS-InSAR and SBAS-InSAR, and the correlation between them in this study reached 93.43%. If the effect of "b" is not considered, the magnitude of "a" indicates the overall deviation between the two datasets. A significant linear correlation exists between the monitoring results obtained from PS-InSAR and SBAS-InSAR measurements.
Regarding statistical significance, the data obtained using the SBAS-InSAR and PS-InSAR processing methods are highly consistent, with a mean error of 1.51 mm/yr. The validation results prove the reliability of the sedimentation rate results obtained in this study.

The Response of Subsidence to Groundwater Level Changes
Water level data from the three groundwater level monitoring points were taken to analyze the relationship between the groundwater level and land deformation. Overall, the time series of surface displacement at the three monitoring points shows a nearly linear distribution, proving that subsidence in the study area was mainly inelastic subsidence, which agrees with the findings of Mathu et al. [80]. The three monitoring points are all located in areas with insignificant subsidence rates and with the water table buried at depths of 9.72-10.01 m at Q62, 43.02-44.31 m at Q96, and 66.65-71.77 m at Q20. The relationship between the water table and land subsidence is plotted in Figure 10. The land subsidence at the three monitoring sites is more related to the decline in the groundwater level, showing that during the groundwater drawdown phase before September 2016 and September 2017, the land also experienced varying degrees of accelerated subsidence.
The parameters in Equation (18) are 0 = 128 ± 37, = −1.64 ± 0.04, = 2.16 ± 0.08, = 8217 ± 314. Figure 9 shows the linear correlation regression analysis of the two datasets. The linear regression equation y = a + bx is obtained, where a = 0.63 ± 0.02 and b = 1.04 ± 0.01, with a correlation coefficient = 0.9343. The closer the correlation coefficient to 1, the higher the correlation between PS-InSAR and SBAS-InSAR, and the correlation between them in this study reached 93.43%. If the effect of "b" is not considered, the magnitude of "a" indicates the overall deviation between the two datasets. A significant linear correlation exists between the monitoring results obtained from PS-InSAR and SBAS-InSAR measurements.
Regarding statistical significance, the data obtained using the SBAS-InSAR and PS-InSAR processing methods are highly consistent, with a mean error of 1.51 mm/yr. The validation results prove the reliability of the sedimentation rate results obtained in this study.

The Response of Subsidence to Groundwater Level Changes
Water level data from the three groundwater level monitoring points were taken to analyze the relationship between the groundwater level and land deformation. Overall, the time series of surface displacement at the three monitoring points shows a nearly linear distribution, proving that subsidence in the study area was mainly inelastic subsidence, which agrees with the findings of Mathu et al. [80]. The three monitoring points are all located in areas with insignificant subsidence rates and with the water table buried at depths of 9.72-10.01 m at Q62, 43.02-44.31 m at Q96, and 66.65-71.77 m at Q20. The relationship between the water table and land subsidence is plotted in Figure 10. The land subsidence at the three monitoring sites is more related to the decline in the groundwater level, showing that during the groundwater drawdown phase before September 2016 and September 2017, the land also experienced varying degrees of accelerated subsidence.  The land subsidence at Q62 is more significant and highly affected by the groundwater table. The groundwater level change curve at Q20 is almost synchronous with the time series curve of surface deformation at this point, while the overall surface deformation at Q96 is small, and its groundwater level is also less variable.
It is worth noting that the magnitude of land subsidence before the completion of the groundwater table decline phase (green straight lines in Figure 10b-d) is arranged in the order Q62 > Q20 > Q96. The difference in the magnitude of subsidence is perhaps related to the groundwater table burial depth and surface building loads. Generally, the shallower the water table depth, the more susceptible the land to water table deformation [81]. This is demonstrated by the fact that the ground settlement rate is also more significant when the water table drops rapidly at Q62, and it is much greater than at Q20 and Q96. However, the time series deformation diagrams for Q96 and Q20 do not conform to this pattern. It can be noted that the difference in the water table burial depth between points Q96 and Q20 is not significant compared to Q62. However, there is a massive difference in the height of the surface buildings at Q96 and Q20, i.e., the height of the surface buildings at Q20 is 60 m. In comparison, the surface buildings at Q96 are only 11 m high, so the settlement at Q20 is greater than at Q96 during groundwater level fall. In addition, the subsidence is severe in the southern part of Qilihe District, located in the Lanzhou fault basin. The Lanzhou Downfaulte basin is an essential water source for the urban area to supply tap water with groundwater due to the abundance of phreatic water containing loose rock types, and it is also a key area for groundwater extraction. However, the groundwater richness of this fault basin gradually decreases from north to south, so the subsidence in the southern part of the region may be due to the weaker groundwater richness and groundwater over-exploitation. The severe decline in groundwater levels during the pre-2018 phase was due to the large amount of groundwater extracted in Lanzhou because of the construction of Metro Line 1; consequently, the rate of surface subsidence in some areas was significantly greater before 2018 than after 2018. Therefore, the relevant authorities should strictly control the groundwater level in the future to prevent severe surface subsidence in Lanzhou, which could damage buildings and underground pipelines.

The Response of Subsidence to Geological Faults and Lithology
The distribution of geological formations and lithology in Lanzhou is shown in Figure 11. Several linear boundaries can be observed in the subsidence rate map of Lanzhou, implying a potential link between subsidence and geological faults. The Leitan River Fault and the Jinchengguan Fault in the study area move in opposite directions on either side of the fault. Figure 12 shows the subsidence rates for eight profiles near the Leitan River Fault and the Jinchengguan Fault. It is evident from Figure 12b,d that subsidence is relatively severe in the western part of the Jinchengguan Fault and the Leitan River Fault, with subsidence rates varying from about −8 to −2 mm/year. In the eastern part of the Jinchengguan Fault and the Leitan River Fault, ground uplift rates vary from 1 to 2 mm/year, suggesting an influence of fault traces on subsidence. Particularly in the eastern part of the Leitan River Fault, we can see a clear blue zone of abnormal uplift, which we can attribute to geological faults.
These findings lead us to several conclusions: (1) The subsidence of the Jinchengguan Fault and the Leitan River Fault has changed dramatically. (2) Spatial subsidence distribution is location-dependent and varies along the same fault. For example, the subsidence on the eastern and western sides of the northern part of the Leitan River Fault is similar. This may be due to the discontinuity in force conductivity associated with the fault structure and different vertical compactions of the non-homogeneous thickness of the compressible soils deposited on either side of the fault [82,83]; (3) The Jinchengguan Fault and the Leitan River Fault may currently be active. These faults should be considered in future urban planning, and technical solutions should be applied to avoid infrastructure damage due to subsidence discrepancies. Table 5. Surface geological features, such as lithology, are essential factors influencing subsidence [84], and apparent lithological control of the spatial subsidence distribution in Lanzhou was detected. As shown in Figure 11, the significant small-scale subsidence funnels in Lanzhou are mainly located in the loess-covered area ( Table 5, Number 1), with a maximum deformation rate of −18.74 mm/yr. There is also a small amount of significant subsidence in sandstone, conglomerate, sand, gravel, muddy siltstone, mudstone, and clay areas ( Table 5, Number 2), with a maximum deformation rate of −15.28 mm/yr. Medium-amplitude large-area subsidence occurs mainly in areas with clayey sandy soil, sandy clay interbedded with sand, and gravel layer lenses ( Table 5, Number 3), with a maximum deformation rate of −7.63 mm/yr. Wet sink loess is the most widespread deposit in Lanzhou, accounting for 82.92% of the total area, mainly occurring in the central part of Qilihe District, the northern and southern parts of Chengguan District, and the eastern and northwestern parts of Xigu District, where the land cover mainly originates from towns and agricultural areas. About 26.94% of the loess has subsided, mostly in the form of small funnels; moreover, all buildings areas of loess have subsided, mainly in Qilihe District and Xigu District. This happens because the soil structure is rapidly destroyed when wet sinking loess is wetted by water under pressure, producing additional subsidence and a rapid strength reduction, so the loess is less stable and has a weaker bearing capacity as a building foundation, becoming prone to subsidence [78]. The main urban area of Lanzhou is mainly composed of moderately compressible soils ( Table 5, Number 3), distributed over 2.09% of the total area. This type of soil subsidence mainly comprises large areas of medium amplitude and occurs over an area of approximately 43.81%. In addition, areas with clay, sandstone, mudstone, or conglomerate geology experience moderate or significant subsidence, with rates ranging from −7.63 to −15.28 mm/yr in areas with urban and agricultural land cover.  These findings lead us to several conclusions: (1) The subsidence of the Jinchengguan Fault and the Leitan River Fault has changed dramatically. (2) Spatial subsidence distribution is location-dependent and varies along the same fault. For example, the subsidence on the eastern and western sides of the northern part of the Leitan River Fault is similar. This may be due to the discontinuity in force conductivity associated with the fault structure and different vertical compactions of the non-homogeneous thickness of the compressible soils deposited on either side of the fault [82,83]; (3) The Jinchengguan Fault and the Leitan River Fault may currently be active. These faults should be considered in future urban planning, and technical solutions should be applied to avoid infrastructure damage due to subsidence discrepancies.

The statistical relationship between lithology and subsidence is shown in
The statistical relationship between lithology and subsidence is shown in Table 5. Surface geological features, such as lithology, are essential factors influencing subsidence [84], and apparent lithological control of the spatial subsidence distribution in Lanzhou was detected. As shown in Figure 11, the significant small-scale subsidence funnels in Lanzhou are mainly located in the loess-covered area ( Table 5, Number 1), with a maximum deformation rate of −18.74 mm/yr. There is also a small amount of significant subsidence in sandstone, conglomerate, sand, gravel, muddy siltstone, mudstone, and clay areas ( Table 5, Number 2), with a maximum deformation rate of −15.28 mm/yr. Mediumamplitude large-area subsidence occurs mainly in areas with clayey sandy soil, sandy clay interbedded with sand, and gravel layer lenses (Table 5, Number 3), with a maximum deformation rate of −7.63 mm/yr. Wet sink loess is the most widespread deposit in Lanzhou, accounting for 82.92% of the total area, mainly occurring in the central part of Qilihe District, the northern and southern parts of Chengguan District, and the eastern and northwestern parts of Xigu District, where the land cover mainly originates from towns and agricultural areas. About 26.94% of the loess has subsided, mostly in the form of small funnels; moreover, all buildings areas of loess have subsided, mainly in Qilihe District and Xigu District. This happens because the soil structure is rapidly destroyed when wet sinking loess is wetted by water under pressure, producing additional subsidence and a rapid

The Response of Subsidence to Rail Transit
A 1000 m buffer zone was established for the railway transit in Lanzhou, and the deformation information within the buffer zone was extracted, as shown in Figure 13. Linear boundaries along the railway line are visible on the Lanzhou subsidence map, particularly in Figure 13 (A,B,D,E), suggesting a possible link between subsidence and railway construction. The combination of static load originating from the weight of trains and dynamic load, due to the dynamic response during traveling, exerts continuous pressure on the railway [85]. Table 6 shows the level of railway subsidence in Lanzhou, ranking the severity of railway subsidence according to location, length, and the percentage of subsidence. The railway exhibits a severe subsidence effect in the case of the following scenarios: high subsidence value, passing through a densely populated area, or when subsidence accounts for a high proportion of the railway. During the monitoring period, all railways in Lanzhou experienced a certain degree of settlement at specific sections. The railways with more severe subsidence are the Baolan Railway (A), the Longhai Railway (B), the Lanyu Railway (B), and the Lanzhou North Ring Road (C), while the rest of the railways have fewer subsidence sections and less subsidence. The Longhai and Lanyu Railways in B show the most severe subsidence. It is subsided in the eastern part of Chengguan District, with maximum subsidence rates of −11.68 and −10.67 mm/yr. The Lanzhou North Ring Road (Figure 13 (C)), passing through Chengguan, Anning, and Xigu Districts, is also severely subsided. This route passes through the extremely severe subsidence of Baidao Ping, Shangping Village, and Lanzhou North Railway Station, with a maximum subsidence rate of −6.18 mm/yr; the remaining sections do not experience severe subsidence. Lanzhou Metro exhibits an uneven deformation pattern, with alternating subsidence and deformation. The subsidence amplitude is more significant than the uplift amplitude, with subsidence accounting for 49.71% of the deformation. The most severe area settles at a rate of −4.13 mm/yr. Figure 13 shows a positive correlation between the distance to the railway and subsidence. If the subsidence level exceeds the maximum allowed for the railway, the elastic deformation of trains at high speed combined with the deformation of the roadbed can be directly reflected in the track, potentially severely affecting the safety of railway operation [86]. Therefore, the railway construction department should consider and monitor quantitative subsidence, especially on the Longhai Railway, Lanyu Railway, and Baolan Railway in the eastern part of Chengguan District, and the Lanzhou North Ring Road in Anning District, to prevent track deformation from threatening the safety of high-speed trains.

The Response of Subsidence to Building Loads
Information on building heights and structures [87] in Lanzhou is based on digitizing Gaofen satellite images and field measurements. A detailed table of building heights and structures in Lanzhou is shown in Table 7. Buildings in Lanzhou are predominantly low-rise and multi-rise, which account for 74.19% of all buildings; high buildings are dominantly mid-rise and high-rise, and account for 25.12%, while super high-rise buildings account for the lowest percentage of 0.69%. The high-rise buildings in Lanzhou are mainly located in Chengguan District, the eastern part of Qilihe Urban District, the southern part of Xigu District, and the eastern part of Anning District. The building structures in Lanzhou are shown in Figure 14; brick-and-mortar houses and frame houses account for 98.15% of the total number of buildings, and the corresponding building structure distribution includes 0.18% old houses, 33.33% frame houses, 64.82% brick-and-mortar houses, 1.66% brickand-timber houses, and 0.00% civil engineering houses. The FR value is calculated using Equation (17) in Section 3.2. The FR reflects the correlation between the subsidence and each building structure, with a higher FR representing a higher probability of subsidence occurring under a particular building structure, i.e., old housing (FR = 1.43), frame houses (FR = 0.89), brick-and-mortar houses (FR = 1.07), brick-and-timber houses (FR = 0.52), and civil engineering houses (FR = 1.89). Thus, civil engineering houses and old houses are most susceptible to subsidence, although brick-and-mortar and frame houses also exhibit a certain level of subsidence susceptibility. This can be explained by the poor construction of old and civil engineering houses, which are prone to structural instability upon exposure to external forces. While brick-and-mortar houses have the highest density, frame houses correspond more to high-rise buildings. As shown in Figure 15, most large, high-rise, and super high-rise buildings have undergone subsidence, since the amplitude of the load imposed on the ground by the building is another factor contributing to the consolidation and high displacement rates [88]. Previous studies showed that taller buildings are more susceptible to subsidence, but a high factor of influence was also found in groups of buildings between 9 and 17 m in height, which is perhaps related to building density (this category has the highest building density). Underground space development and the excavation of building pits are severe in areas with high building density. The extraction of underground liquid support and the excavation of solid support may induce a stress disbalance in geotechnical soil, readily leading to surface deformation [89]. In summary, building structure, height, and density all influence the subsidence phenomenon. When a settlement exceeds the bearing limit of the building, wall cracking, house water seepage, house tilt collapse, and other disasters that threaten human lives and property safety may readily appear [90].

Future Land Subsidence Risk Assessment
The weights of factors calculated using the random forest model are 0.498 for the land subsidence rate, 0.213 for cumulative land subsidence, 0.054 for lithology, 0.062 for the distance to faults, 0.035 for the distance to rivers, 0.026 for the distance to railways, 0.036 for building height, and 0.076 for precipitation. The higher the weight, the more influence the factor has in terms of land subsidence. The ArcGIS spatial analysis method was used to superimpose the factors and obtain a comprehensive subsidence risk evaluation score. The evaluation score was divided into five categories of very low, low, medium, high, and very high sensitivity, according to the natural breakpoint method, to obtain the hazard zoning map shown in Figure 16. Here, we obtain static hazard zoning, assuming that the city will not expand in the future, the building density does not change, and groundwater abstraction does not increase. However, the building density and groundwater abstraction in Lanzhou are expected to increase with population, which could exacerbate the risk in the area. According to the risk map, the total areas of very high, high, medium, low, and very low risk are 22.02, 54.46, 97.54, 114.95, and 53.00 km 2, respectively. The possible very high-risk areas for subsidence are distributed along the railway line in I Anning District, II Southern Xigu, III Qilihe urban area, and IV Eastern Chengguan District, which is almost consistent with our subsidence monitoring results in Lanzhou from 2014 to 2021. This may be related to the greater weighting of the two influencing factors-the land subsidence rate and cumulative land subsidence. In addition, new settlements will be generated by lateral expansion based on the current settlement in the future.

Future Land Subsidence Risk Assessment
The weights of factors calculated using the random forest model are 0.498 for the land subsidence rate, 0.213 for cumulative land subsidence, 0.054 for lithology, 0.062 for the distance to faults, 0.035 for the distance to rivers, 0.026 for the distance to railways, 0.036 for building height, and 0.076 for precipitation. The higher the weight, the more influence the factor has in terms of land subsidence. The ArcGIS spatial analysis method was used to superimpose the factors and obtain a comprehensive subsidence risk evaluation score. The evaluation score was divided into five categories of very low, low, medium, high, and very high sensitivity, according to the natural breakpoint method, to obtain the hazard zoning map shown in Figure 16. Here, we obtain static hazard zoning, assuming that the city will not expand in the future, the building density does not change, and groundwater abstraction does not increase. However, the building density and groundwater abstraction  Figure 17. The areas at very high risk of subsidence in each administrative district are ranked as follows: IV Chengguan District > II Xigu District > III Qilihe District > I Anning District. The areas at very high risk and high risk for subsidence in Chengguan District are larger than in other districts, making Chengguan District the place at greatest risk of future subsidence. As the population of Chengguan District is as high as 1.484 million, almost 2-7 times that of the other districts, this finding is of great significance for the state's subsidence management in this area. Specific very high subsidence risk locations are shown in Table 8. Our subsidence risk maps and risk tables for the investigated districts are an important resource for developing the most effective low-risk management strategies and future urban development plans [51]. Policymakers can use this when assessing subsidence hazard zones in urban development to identify very high-risk and high-risk hazard-prone areas within each district. consistent with our subsidence monitoring results in Lanzhou from 2014 to 2021. This may be related to the greater weighting of the two influencing factors-the land subsidence rate and cumulative land subsidence. In addition, new settlements will be generated by lateral expansion based on the current settlement in the future. The statistic data for risk level in each administrative district are shown in Figure 17. The areas at very high risk of subsidence in each administrative district are ranked as follows: Ⅳ Chengguan District > Ⅱ Xigu District > Ⅲ Qilihe District > Ⅰ Anning District. The areas at very high risk and high risk for subsidence in Chengguan District are larger than that of the other districts, this finding is of great significance for the state's subsidence management in this area. Specific very high subsidence risk locations are shown in Table  8. Our subsidence risk maps and risk tables for the investigated districts are an important resource for developing the most effective low-risk management strategies and future urban development plans [51]. Policymakers can use this when assessing subsidence hazard zones in urban development to identify very high-risk and high-risk hazard-prone areas within each district. We counted a total of 51,163 houses in the very high-risk area. According to the administrative division, 6383 houses in Anning District, 20,068 houses in Chengguan District, 15,553 houses in Qilihe District, and 9171 houses in Xigu District are at very high risk of subsidence, accounting for 12. 48, 39.22, 30.40, and 17.93% of the number of houses at very high risk of subsidence, respectively. According to the classification of building structures, 1 civil house, 213 brick-and-timber houses, 43,180 brick-and-mortar houses, 7325 frame houses, and 435 old houses are located in very high-risk areas. About 44.47% of brick-and-timber houses and 51.36% of old houses are located in very high-risk areas. These two types of house are structurally unstable, and land subsidence can aggravate the impact on them, highlighting the need for them to be demolished and rebuilt or strengthened. According to the building height classification, 3242 high-rise buildings and 730 super high-rise buildings are in very high-risk areas. About 14.09% of high-rise buildings and 52.78% of super high-rise buildings are at very high risk of subsidence. Lanzhou is not suitable for further super-tall buildings, and the government needs to focus on assessing whether these super-tall buildings are cracking or leaning.   43,180 brick-and-mortar houses, 7325 frame houses, and 435 old houses are located in very high-risk areas. About 44.47% of brick-and-timber houses and 51.36% of old houses are located in very high-risk areas. These two types of house are structurally unstable, and land subsidence can aggravate the impact on them, highlighting the need for them to be demolished and rebuilt or strengthened. According to the building height classification, 3242 high-rise buildings and 730 super high-rise buildings are in very high-risk areas. About 14.09% of high-rise buildings and 52.78% of super high-rise buildings are at very high risk of subsidence. Lanzhou is not suitable for further super-tall buildings, and the government needs to focus on assessing whether these super-tall buildings are cracking or leaning.
Upon comparing the Lanzhou subsidence risk map ( Figure 16) with the Lanzhou subsidence monitoring map (Figure 4), we found differences between the two, and we added six new high-and very high-risk zones from n1 to n6. In addition, there are several risk assessment areas where errors exist. The predicted subsidence area in part H was significantly smaller than the actual monitoring results because no monitoring points were identified at this location, resulting in no subsidence rate or cumulative settlement factor values included in the calculation, ultimately resulting in lower predicted values; in reality, it this area should have an extremely high value. The western part of Chengguan District (n6) was assessed as a high-risk subsidence area, contrary to the current deformation monitoring results. In Section 4.3.2, we discussed that the uplift in the western part of Chengguan District near n6 might be due to the fault's movement. We cannot predict the trend in the fault movement from these data, but it is worth noting that the area in the transition between the uplifted area in the western part of Chengguan District and the subsidence area in the east may lead to damage to and the collapse of buildings due to differential subsidence [91,92]. That is, the accumulation of uneven settlement over a distance is a current hazard in the Chengguan District of Lanzhou, which may endanger the structural integrity of buildings.

Conclusions
In this study, we used the SBAS-InSAR method to analyze and interpret a dataset of 150 Sentinel-1 SAR scenes in Lanzhou from 26 October 2014 to 12 December 2021 to describe the current patterns, rates, impact factors, and long-term predictions of land subsidence in Lanzhou. The main conclusions are as follows: (1) The average annual deformation rate in Lanzhou ranges from −18.74 to 12.78 mm/yr.
The land subsidence pattern is uneven, and it occurs mainly in the form of small subsidence funnels. The land subsidence dynamics in Lanzhou involve horizontal expansion followed by vertical development. During the monitoring period, Lanzhou was dominated by linear subsidence. The overall deformation situation in Lanzhou can be summarized as one uplift area and four subsidence areas, namely, the subsidence area along the railway in Anning District, the southern Xigu subsidence area, the Qilihe urban subsidence area, the eastern Chengguan District subsidence area, and the western Chengguan District uplift area. Most of the significant subsidence occurs along the Yellow River and railway and in villages and towns at the edge of the urban area. (2) Factors contributing to subsidence in Lanzhou include groundwater extraction, river erosion, stratigraphic lithology, fault distribution, underground pressure, and pressure from high-rise buildings. The fall in the groundwater table inevitably leads to land subsidence, but the extent of the subsidence is mainly influenced by the water table burial depth and surface building loads. A correlation between land subsidence and geological faults in Lanzhou can be established. In the presence of active faults, the land movement pattern is characterized by an uplift on one side of the fault and subsidence on the other side; this phenomenon occurs both in the Leitanhe Fault and the Jinchengguan Fault. The spatial subsidence distribution in Lanzhou is under apparent lithological control; significant subsidence mainly occurs in areas covered by loess, and the remaining large subsidence areas are almost in more compressible soils. The land subsidence in Lanzhou is strongly linked to railway lines, with subsidence mostly distributed along the railway. Moreover, Lanzhou is prone to subsidence funnels near high-rise buildings, and a high building density also causes land subsidence. (3) The existing subsidence in Lanzhou will continue to occur, and may even trigger the cracking and collapse of buildings. The area at very high risk of future subsidence in Lanzhou is 22.02 km 2 , while the area at high risk of subsidence is 54.47 km 2 . A total of 51,163 houses are at very high risk of subsidence, and they are mainly located in Chengguan District and Qilihe District, accounting for 39.22% and 30.40% of the buildings at very high risk of subsidence, respectively. In total, 44.47% of brick-andtimber houses, 51.36% of old houses, and 52.78% of super-tall buildings are at very high risk of subsidence. In addition, new subsidence zones are expected to emerge due to the active use of railways, high-rise buildings, ground lithology, faults, rivers, precipitation, and other factors. This study maps the subsidence risk of each building and provides a list of areas with very high subsidence risk, which should be the focus of government efforts to prevent and control possible accidents in these buildings. (4) In future urban planning and disaster prevention, the findings of this systematic study-i.e., that the loess-covered areas of Lanzhou are not suitable for mediumrise, high-rise, and super high-rise buildings, while the central city of Lanzhou is not suitable for too many super high-rise buildings-should be considered. Since the ground is highly susceptible to subsidence where the water table is shallow, the changes in groundwater levels in the area should be controlled. At the same time, attention should be paid to buildings in the transition area between the western uplift area and the eastern subsidence area in Chengguan District to prevent their cracking and collapse due to differential settlement. In addition, the authorities should pay attention to the list of buildings at high risk of settlement and the areas at very high risk of settlement that were determined and drawn up in this paper, and focus on preventing possible accidents related to these buildings.
Author Contributions: All authors contributed to the manuscript and discussed the results. Y.X. drafted the manuscript and was responsible for the data processing, and the analysis and interpretation of the results. Z.W. proposed the ideas for the thesis, designed the structure and contributed to the final revision of the thesis. H.Z. revised the manuscript and contributed to the final revision of the thesis. J.L. provided the code for the random forest model and directed the risk assessment section. Z.J. collected the literature. All authors have read and agreed to the published version of the manuscript.
We would like to express our gratitude to the editors and reviewers for their helpful comments on earlier versions of this paper.

Conflicts of Interest:
The authors declare no conflict of interest.