Characteristics of Surface Deformation in Lanzhou with Sentinel-1A TOPS

Abstract

While surface deformations and their impact on buildings have been observed in the city of Lanzhou, it is difficult to find studies of surface deformation and the influential factors in the recent decades. This study was designed to detect the spatial position of these surface deformations and to understand the mechanism behind them. Sentinel-1A TOPS model image data acquired in descending orbits between March 2015 and May 2019 were processed by using Small Baseline Subsets Interferometry (SBAS)-Interferometric Synthetic Aperture Radar (InSAR) technology, and then compared with geology, tectonical aspects of the study area and land cover types in Lanzhou. The results reveal that the land surface deformation is uneven in this city, and seven surface deformation regions were detected in the following areas: the north freight yard, Jiuzhou, Country garden, Donggang, Yanjiaping, Zhongxinping and Liuquan town. The land surface deformation rate in Lanzhou ranges from −82.13 mm/year to 19.31 mm/year. Time-series land surface deformation analysis showed that deformation increased over time in major deformation regions. Surface deformation expansion was significant after June 2017and it continued to expand. The surface deformation of Lanzhou is affected by natural factors (geology and geological faults) and human activities (land cover types/land cover changes). Local geological conditions control the location of the surface deformation process. These findings provide compelling data and theoretical support for disaster prevention and reduction in Lanzhou.

1. Introduction

Land subsidence is a slow and long-term geological hazard [1,2]. In many cities across the globe, the ground is subsiding at an annual average rate of several to tens of centimeters [3,4,5]. In most cities, land subsidence is due to human activities [6]. Urban land subsidence has caused the deformation of road surfaces and subgrades, even damaging buildings and urban infrastructure, resulting in significant casualties and economic losses [7]. Therefore, to maintain the sustainability of urban areas, it is very important to understand the deformation pattern and evolution of municipal ground [8].

To monitor the surface deformation accurately, academic researchers and industry personnel have studied and developed technology that can be utilized to monitor surface deformation. Currently, the technologies used for monitoring surface deformation mainly include leveling, the Global Positioning System (GPS) and Interferometric Synthetic Aperture Radar (InSAR). Among these, the precision of leveling and GPS technology for ground monitoring meets the requirements needed to monitor urban surface deformation, but the monitoring process is time-consuming, which is not conducive to the rapid emergency response required to address urban disasters [1,9]. GPS observes the shape variable in the three-dimensional direction, while InSAR technology can obtain the deformation on the oblique path. Moreover, with the increase in the temporal resolution of microwave data, InSAR technology can quickly obtain a large range of surface deformation information and measure the surface displacement of large areas (regional and local scale) with millimeter accuracy over a long period of time [10], and it is widely used in urban ground subsidence monitoring [11]. The development of InSAR technology has led to the development of different techniques, such as Differential Interferometry (D-InSAR), Permanent Scatterer Interferometry (PS-InSAR) and Small Baseline Subsets Interferometry (SBAS-InSAR) [12,13,14]. Among these, SBAS-InSAR is most effective in detecting surface deformation information for unstable or mountainous areas [4,15]. Moreover, SBAS-InSAR technology is an improved InSAR time-series analysis tool that effectively solves baseline time discontinuous problems [16]. SBAS-InSAR technology overcame limitations of the conventional InSAR, mostly due to decorrelation and compensated phase erroneous contributions owing to atmospheric phase delay, inaccurate topographic model and uncertain satellite orbits.

Lanzhou is a city in China that is placed in the region that has the thickest loess deposit in the world. However, loess geology is unstable, and various natural environmental factors are very sensitive to its impact. Under the guidance of Chinese policies, such as reform and opening up, Western development and “One Belt and One Road”, land use in Lanzhou has significantly changed [17], and the level of urban construction and urbanization has significantly improved [18]. Still, large-scale mountain-cutting for land use has led to severe surface deformation. Some of the observed surface deformations and their impact on buildings are shown in Figure 1a–d. Ground collapses and sinks are clearly seen in Figure 1a,b. Evidence of deformation due to creating new land areas by cutting through the mountains is shown in Figure 1c. Figure 1d shows the presence of a significant number of cracks and openings in the wall of a building. In a city such as Lanzhou, which is prone to natural disasters, it is necessary to monitor surface deformation changes.

This study was designed to detect the spatial position of surface deformation in Lanzhou and understand its mechanism. Toward that end, Sentinel-1A TOPS model image data covering the region were chosen. Sentinel-1A data that were acquired in descending orbits between March 2015 and May 2019 were processed using SBAS-InSAR. The velocity fields obtained with SBAS-InSAR were compared with the geology, tectonic characteristics of the study area and the types of land cover found in Lanzhou.

2. Data and Methods

2.1. Study Area

Lanzhou is situated in the western part of China on a loess plateau. The urban area is about 1631.6 km². The following municipal areas of Lanzhou were investigated in this study: Chengguan District, Qilihe District, Xigu District, Anning District and the land-cutting area of Gaolan County. The location of the study area is presented in Figure 2.

The terrain of the area is high in the west and low in the east. The area has a temperate continental climate; the temperature difference is large about 20 ℃. The annual amount of precipitation is low about 300 mm. The Yellow River runs through the main urban area, and the city is situated in the plain of the eroded and accumulated Yellow River valley; thus, it is a typical valley city situated in a mountainous region. Due to its special topography, landform, geology, climate conditions and aggressive human activities, geological disasters occur frequently in Lanzhou, such as landslides, building collapse, ground collapse, etc. Because Lanzhou is located in the loess region, concentrated precipitation may lead to rapid destruction of soil-soaked ground, resulting in significant deformation.

2.2. Data Sources

A total of 75-scene descending orbits obtained from Sentinel-1A TOPS SAR images (C-band) of Lanzhou acquired from 14 March 2015 to 4 May 2019 were selected to estimate the average vertical surface subsidence velocity and subsidence time-series (https://search.asf.alaska.edu/#/). The Sentinel 1A data are shown in

Table 1. Sentinel-1A data used in the study area.

Acquisition Date	Acquisition Date	Acquisition Date	Flight Direction
14 March 2015	2 November 2016	28 December 2017	Descending
7 April 2015	9 November 2016	9 January 2018	Descending
1 May 2015	26 November 2016	21 January 2018	Descending
25 May 2015	3 December 2016	2 February 2018	Descending
18 June 2015	27 December 2016	26 February 2018	Descending
12 July 2015	20 January 2017	10 March 2018	Descending
5 August 2015	6 February 2017	9 May 2018	Descending
29 August 2015	13 February 2017	21 May 2018	Descending
22 September 2015	15 March 2017	2 June 2018	Descending
16 October 2015	20 April 2017	14 June 2018	Descending
9 November 2015	2 May 2017	26 June 2018	Descending
3 December 2015	26 May 2017	8 July 2018	Descending
20 January 2016	7 June 2017	13 August 2018	Descending
13 February 2016	19 June 2017	6 September 2018	Descending
8 March 2016	13 July 2017	30 September 2018	Descending
1 April 2016	25 July 2017	12 October 2018	Descending
25 April 2016	6 August 2017	24 October 2018	Descending
19 May 2016	18 August 2017	5 November 2018	Descending
12 June 2016	30 August 2017	4 January 2019	Descending
30 July 2016	11 September 2017	28 January 2019	Descending
23 August 2016	23 September 2017	5 March 2019	Descending
16 September 2016	5 October 2017	29 March 2019	Descending
21 September 2016	10 November 2017	10 April 2019	Descending
27 September 2016	22 November 2017	22 April 2019	Descending
16 October 2016	4 December 2017	4 May 2019	Descending

. Shuttle Radar Topography Mission (SRTM) with 90 m resolution provided by the National Aeronautics and Space Administration (NASA) was used to correct the topography-related phase and geocode interferograms. SBAS-InSAR was used to process the 75-scene descending Sentinel 1A images (Table 1) [12]. PS-InSAR was used to process 32-scene descending Sentinel 1A images covering the study area from March 2015 to January 2017 to perform a cross-validation of the differences between the SBAS-InSAR and PS-InSAR results [13,15]. In our experiment, these were both incorporated into the ENVI 5.5 platform through the SARscape 5.5 module. Most of the time baselines of the data were within one month, and all the spatial baselines of the image pairs complied with the requirements.

Geological data were purchased from the Geo-cloud of China Geological Survey (http://geocloud.cgs.gov.cn/#/portal/home), which includes the location of geological faults and the types of geology found the region. The land cover types were downloaded from Tsinghua University’s global 10 m resolution land cover maps [19].

2.3. Methods and Processing

2.3.1. SBAS-InSAR Algorithm

The SBAS-InSAR was proposed by Berardino [12] and Lanari et al. [20]. SBAS-InSAR time-series analysis method combines the data, and it can be used to obtain the short space baseline differential interferogram dataset. These differential interferograms can overcome spatial decorrelation phenomena. Using Singular Value Decomposition (SVD) to solve the deformation rate, isolated SAR datasets separated by large spatial baselines can be connected to improve the time sampling rate of the observed data. Its high-density temporal and spatial information can effectively eliminate the atmospheric effect phase, resulting in more accurate measurements [16]. The detailed principles are as follows:

(1) N+1 SAR images covering the study area are arranged by followed the time series $t_0$ to $t_N$, then one of them is select as the super master image, and other slave SAR images are co-registered to this super master image. M multi-visual differential interferograms are generated based on N+1 SAR images. M satisfies the following inequality:

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

(1)

(2) The $j_{th}$ differential interferogram is generated from the slave SAR image obtained at the time of $t_A$ and the master SAR image obtained at the time of $t_B$ (t_B > t_A). The interference phase in a pixel of azimuth and range coordinates (x,r) can be expressed as:

$$\begin{array}{c}\delta {\varphi }_j\left(x,r\right)={\varphi }_B\left(x,r\right)-{\varphi }_A\left(x,r\right)\\ \approx \frac{4\pi }{\lambda }\left[d\left(t_B,x,r\right)-d\left(t_B,x,r\right)\right]+\Delta {\varphi }_{topo}^j\left(x,r\right)+\Delta {\varphi }_{atm}^j\left(t_B,t_A,x,r\right)\hfill \\ +\Delta {\varphi }_{noise}^j\left(x,r\right)\hfill \end{array}$$

(2)

where j $\in$ (1, …, M); $\lambda$ is the central wavelength of the signal; $d\left(t_B,x,r\right)$ and $d\left(t_B,x,r\right)$ are the cumulative shape of the radar line of sight for $d\left(t_0,x,r\right)=0$ at the time of $t_B$ and $t_A$, respectively; $\Delta {\varphi }_{topo}^j\left(x,r\right)$ is residual topographic phase in differential interferograms.$\Delta {\varphi }_{atm}^j\left(t_B,t_A,x,r\right)$ is atmospheric delayed phase; $\Delta {\varphi }_{noise}^j\left(x,r\right)$ is coherent noise such as system thermal noise.

(3) To obtain the sedimentation sequence with physical significance, the phase in the Formula (2) is expressed as the product of the average phase velocity and time between two obtained times:

$$v_j=\frac{{\varphi }_j-{\varphi }_{j-1}}{t_j-t_{j-1}}$$

(3)

the $j_{th}$ phase value of the interferogram can be expressed as:

$${\displaystyle \sum }_{k=t_{A,j}+1}^{t_{B,j}}\left(t_k-t_{k-1}\right)v_k=\delta {\varphi }_j$$

(4)

That is, the velocity of each period needs to be integrated over the time interval between the master and slave images. The form of matrix is:

$$Bv=\delta \varphi$$

(5)

Formula (5) is a matrix of M $\times$ N. Because the differential interferograms with small baseline sets adopt the multi-principal image strategy, matrix B is prone to rank defect. Using the SVD method, the generalized inverse matrix of matrix B can be achieved, and then the minimum norm solution of the velocity vector can be obtained. Finally, shape variables of each time period can be obtained by integrating the velocities in each time period.

This study used the 75-scene Sentinel-1A SLC images covering the study area from March 2015 to May 2019. The 22 September 2015 image was automatically selected as the super master image, and other slave SAR images were co-registered and resembled to the master image; its maximum time baseline is 200 days, its range lookup is 4 and its azimuth looks is 1. Meanwhile, interferometer pairs with low coherence were removed. Finally, the connection graph was obtained (Figure 3). The interference process mainly includes interferogram generation, de-flattening, filtering, coherent graph generation and phase unwinding. SRTM DEM (Digital Elevation Model) was utilized to remove the terrain phase and registration. The Goldstein filtering method was utilized to remove noise. The polynomial refinement method was used to flatten the image. The Minimum Cost Flow (MCF) method was used only for phase unwrapping of the wrapped figure. Linear model with better robustness was used to estimate the deformation rate and to remove the atmospheric phase to obtain the final deformation rate. Ultimately, the deformation rate of the Line-of-Sight (LOS) converts to the vertical direction.

2.3.2. PS-InSAR Algorithm

The D-InSAR algorithm establishes an interference phase model based on multiple SAR image sequences covering the same area by selecting PS points, and settlement results are achieved by interpreting the model [21]. PS-InSAR was proposed by Ferretti [10]; The PS-InSAR method utilizes a time-series of radar images to identify highly coherent points, the so-called persistent scatterers (PS) [10]. It uses phase stable points to extract deformation information. Because of its coherence, it is rarely affected by time and space baselines; thus, it is possible to prevent the limitation of time and space incoherence of traditional D-InSAR technology [22]. Basically, by using the same area of the SAR images, PS-InSAR uses space borne radar to detect the stability of the point target. Because PS-InSAR technology relies on the point target to analyze the change in the time-series, the point target in the time-series is nearly unaffected by the speckle noise in the SAR image. Moreover, after a long period of time, very high coherence and scattering properties remain stable. In the present study, the temporal and spatial changes of these stable points were analyzed, and the phase composition of each point, including the elevation error, deformation phase and errors caused by atmospheric disturbance, was decomposed. Finally, through iterative regression analysis of the interference phase of these stable points, the time-series variables were calculated, and the accumulation of ground settlement was obtained [23].

PS-InSAR technology is mainly concerned with the discrete spatial distribution of PS points. The i-th image in the time series data sets containing N single look complex (SLC) SAR images was marked as s_i (i = 1, 2, …, N). Interference pairs composed of the i and k SAR images can be recorded as $I^{i,k}=s_i{s}_k^{*}$. Here, * stands for the complex conjugate operation. Assuming that the target p₀ is the reference point, the interferometer phase $\Delta {\varnothing }_{p,}^{i,k}=\angle I_{p,p0}^{i,k}$ will be dependent on the elevation, deformation, atmospheric disturbance and noise difference between the target p and the reference point p₀. The phase difference caused by the elevation and line difference between the two targets can be expressed as:

$$\Delta {\varphi }_{H,p,p_0}^{i,k}=\frac{4\pi }{\lambda }\frac{1}{R\sin \theta }\Delta h_{p,p_0}{B}_n^{i,k}$$

(6)

$$\Delta {\varnothing }_{D,p,p_0}^{i,k}=\frac{4\pi }{\lambda }\Delta v_{p,p_0}{B}_t^{i,k}$$

(7)

In the formulae, $\Delta h_{p,p0}$ and $\Delta v_{p,p0}$ is the difference of elevation and deformation rate between target $p$ and $p_0$. $B_n^{i,k}$ and $B_t^{i,k}$ represent the spatial baseline and time baseline perpendicular to the radar line of sight between SAR images i and k, respectively. $\lambda$ is the radar signal wavelength, R is the oblique distance between the target and the satellite, and $\mathsf{\theta }$ is the incident angle of SAR.

The target elevation and its subsidence rate were usually estimated by maximizing time coherence:

$$\Delta 1_h^{\wedge }{}_p,\Delta 1_v^{\wedge }{}_p=\arg \left\{\max \left(\left|{\xi }_p\right|\right)\right\}$$

(8)

$${\xi }_p=\frac{1}{M}{{\displaystyle \sum }}_{i,k}{e}^{j\left(\Delta {\varphi }_p^{i,k}-\Delta {\varphi }_{H,p}^{i,k}-\Delta {\varphi }_{D,p}^{i,k}\right)}$$

(9)

where, M is the number of interferograms. $\Delta {\varphi }_p^{i,k}$ are the differential interferometer phase. $\Delta {\varphi }_{H,p}^{i,k}$ is the elevation correlation partial phase, and $\Delta {\varphi }_{D,p}^{i,k}$ is the terrain deformation correlation phase.

The land deformation $\Delta {\varphi }_{D,p}^{i,k}$ is mainly composed of linear deformation phase and nonlinear deformation phase.

$${\varphi }_{D,p}^{i,k}={\varphi }_{linear}^{i,k}+{\varphi }_{non-linear}^{i,k}=\frac{4\pi }{\lambda }\left(\Delta v_p{B}_t^{i,k}+Def{o}_{non-linear}\right)$$

(10)

where, ${\varphi }_{non-linear}^{i,k}$ is a nonlinear variable. In Formula (9), the maximization factor $\xi =\max \left(\left|{\xi }_p\right|\right)$ was called time coherence and the time coherence of the target can be approximately expressed as $1_{\xi }^{\wedge }{}_p\approx e^{-{\delta }_{\varnothing }^2/2}$ a function of the deviation of the phase residual. Finally, the variance of elevation and deformation rate is calculated by PS-InSAR technology.

This study also used 32-scene Sentinel-1A SLC images covering the investigated area from March 2015 to January 2017. Firstly, a super master image was automatically selected by the program, the 20 January 2016 image was selected as the super master image, and other slave SAR images were co-registered and resembled to the master image, the coordinate system was converted to the master image coordinate system. Secondly, the connection graph was generated (Figure 4) and the interference was carried out. The interferometer process mainly includes interferogram generation, PS points selection, de-flattening, filtering, phase unwrapping, geometric ground control points (GCPs) selection, atmospheric phase removal, and PS points estimation and geocoding. DEM was applied to flatten the terrain phase. If there was an error in the geographical position of the DEM, GCPs were needed to correct the SAR images and DEM. The coherence threshold was 0.85 to improve the coherence in interferometer procedure. Then, the number of GCPS selected in no deformation areas was 40 to remove the phase of atmosphere and noise. Meanwhile, the interferometer pairs with low coherence were removed. The Goldstein filtering method was used to remove noise. The polynomial refinement method was used to flatten the image and eliminate the orbital phase. The three dimensions (3D) method was used for phase unwrapping of the wrapped figure. Finally, we could achieve the time series of PS points. Ultimately, the deformation rate of the line-of-sight (LOS) converts to the vertical direction.

2.3.3. Mann–Kendall Test

The Mann–Kendall (M–K) test is a non-parametric statistical method that was proposed and developed by Mann and Kendall [24]. Its advantages are that it does not require the samples to follow a specific distribution, and it is not impacted by a small number of outliers. It is more suitable for type variables and order variables. The advantage of this method is that the calculation is simple, and the start time of the mutation can be specified. This mutation monitoring method is widely used [25]. The method is as follows:

The time series is random and independent, the defined statistics are expressed as:

$${UF}_k=\frac{s_k-E\left(s_k\right)}{\sqrt{Var\left(S_k\right)}},(k=1,2,\dots ,n)$$

(11)

Where, UF_k is the standard normal distribution, and the sequence of statistics is calculated by the order of time series x₁, x₂, …, x_n. s_k is the cumulative number, $E\left(s_k\right)$ is the mean value, and $Var\left(S_k\right)$ is the variance.

If you follow the time series x in reverse order x_n, x_n-1, … x₁, the above process was repeated to calculate UB_k. UF_k and UB_k curves are drawn, respectively. If the UF_k and UB_k curves intersect at a point between the borderline, then the moment at which they intersect is the abrupt point.

3. Results and Analysis

3.1. Validation of the Results

To quantitatively determine the accuracy of the land surface deformation, a cross-validation of the differences between the PS-InSAR and SBAS-InSAR results were used. The two sets of SAR data had overlapped part of the time series from March 2015 to January 2017, thus the differences during the overlapped timing sequence were used to test these two result sets. The Raster Calculator was utilized to calculate the differences in the PS-InSAR and SBAS-InSAR results (Figure 5).

To quantitatively analyze the differences between the PS-InSAR and SBAS-InSAR results, we counted the number of distribution and determined the cumulative percentage (

Table 2. Statistical analysis results.

Interval Value	Counts	Cumulative Percentage
−75	1	0.0049%
−65	3	0.0148%
−55	0	0.0000%
−45	5	0.0246%
−35	15	0.0739%
−25	67	0.3302%
−15	627	3.0902%
−5	8182	40.3253%
5	9843	48.5116%
15	1403	6.9147%
25	114	0.5619%
35	15	0.0739%
45	11	0.0542%
55	3	0.0148%
65	0	0.0000%
75	1	0.0049%

). The verification results were mainly distributed from −5 mm to 5 mm, which accounted for a cumulative percentage of 88.83%. The average value of the differences between the PS-InSAR and SBAS-InSAR results covering the study area was 0.14 mm, the maximum and minimum values were 72.51 mm and −75.87 mm, respectively. It was found that most of the validation results in the study area were less than 10 mm, indicating that the two sets of data were well tested, and ground deformation rate calculations were reliable. However, the recorded differences between SBAS-InSAR and PS-InSAR results may be explained by errors from the two algorithms. The PS density is usually low in vegetated, forested and low-reflectivity areas [26], while the SBAS performs better in nonurban vegetated areas, and also in areas with high deformation rates [4]. Furthermore, the research team went to the field to investigate surface deformation in Lanzhou (Figure 5). Country Garden and Jiuzhou were more severely deformed, and the types of deformation mainly were cracks, subsidence and ground collapses. Surface deformation locations were consistent with the InSAR results.

3.2. Spatial Distribution of the Surface Deformation

InSAR is more sensitive to vertical than horizontal movements. We convert LOS (d_LOS) into vertical displacement (d_v) for every time-series using the Sentinel-1A incidence angle (θ = 39.58°): d_v = d_LOS/cosθ. In the rest of this paper we use vertical deformation rates converted from the observed LOS rates. The mean deformation rate map was generated using the SBAS-InSAR and Sentinel 1-A images (Figure 6). The land surface deformation rate of Lanzhou ranges from −82.13 mm/year to 19.31 mm/year. Surface deformation rates of most areas (87.23%) range from −10 mm/year to 10 mm/year, which indicates that the land surface is stable in this city. However, some serious surface deformation locations were detected in the following regions: the north freight yard (P1), Jiuzhou (P2), Country Garden (P3), Donggang (P4), Yanjiaping (P5), Zhongxinping (P6) and Liuquan (P7) (Figure 6). Among these seven regions, P3 was identified as having the largest deformation, and its surface deformation range was also the largest. A more detailed surface deformation analysis and explanation for regions P1−P7 will be presented in Section 4.

3.3. Time-Series Land Surface Deformation Map

The 10 images have been selected among 75 images at intervals of about six months to show the time-series land surface deformation and reduce the number of images presented in Figure 7. At the beginning of March 2015, there were no deformations in Lanzhou; however, in December 2015 surface deformation increased over time in the major deformation regions, including Country Garden (P3), Jiuzhou (P2), Liuquan (P7) and Yanjiaping (P5). In December 2016, surface deformation was aggravated in Country Garden (P3), Jiuzhou (P2), the north freight yard (P1), Liuquan (P7), Yanjiaping (P5), Zhongxinping (P6) and Donggang (P4). In December 2017, the deformation and deformation range of the Xigu District increased rapidly. Six months later, this trend was holding back and it was gradually eliminated. In May 2019, the ultimate deformation was formed and distributed in the regions under study. More details on the distribution of the surface deformation for Regions P1−P7 are shown in Figure 8.

3.4. Mutation Analysis of Surface Deformation for Regions P1−P7

Statistics of UF and UB were calculated, and the UF and UB curves were plotted. M-K mutation test results (Figure 9) show that mutations of deformation occurred in P1 and P2 in July 2017 and September 2017, respectively, while the mutations of deformation occurred in P3 in September 2018. They occurred in P4 in May 2017, and in P5, P6 and P7 in June 2017. The significance level for all these regions was P > 0.05, indicating that surface deformation expansion was significant after June 2017, and it continued to expand. Surface deformation M–K test results are in good agreement with the time-series analysis results for Regions P1−P7 (Figure 8).

4. Discussion

The reasons for land surface deformation are complex due to the complex geological environment in Lanzhou and human activities (land cover/land cover changes) in that city. We explored how geology, geological faults, land cover types and land cover changes caused the observed land surface deformation in the region.

4.1. Relationship between Land Surface Deformation and Geology

The geological data of the studied area showed five main types of geology (Figure 10):

1. The Upper Pleistocene Aeolian Group, which mainly consists of loess.

2. The Xiliugou Group, which mainly consists of red massive loose sandstone, grey-white fine conglomerate and glutenite.

3. The Hekou Group, which mainly consists of sandstone, conglomerate, glutenite, argillaceous siltstone and mudstone.

4. The Gansu Group, which consists of mudstone, sandy mudstone and glutenite mixed with marlstone.

5. The Holocene Aeolian Group, which mainly consists of sub-sand soil, a sub-clay sand layer and a gritty layer lens body.

Bozzano et al. [27] indicated that local geological conditions control the magnitude of the deformation process. In this section, we discuss the relationship between land surface deformation and geology in regions P1−P7.

In Region P1, the geology is part of the Holocene Aeolian Group and the Gansu group, which consist of materials, such as sub-sand soil, a sub-clay sand layer, a gritty layer lens body, mudstone, sandy mudstone and glutenite mixed with marlstone. The gravel layer is dominated by medium gravel, with moderate weathering and good roundness [28]. The gravel layer has high shear strength, fast compaction under external load, little deformation and high bearing capacity [29]. The shear strength of the mudstone and sandstone is linked to their lithology and water content. Rock masses with large water content are weak, and they show plastic failure and have low shear strength [30]. The hydrological properties of the materials are mainly swelling and disintegration after encountering water. The repeated changes of the wet and dry cycle will destroy the mudstone [31]. Therefore, Region P1 had serious surface deformation. In Region P2, the geology is mainly part of the Xiliugou Group and the Upper Pleistocene Aeolian Group, which consist of materials, such as red massive loose sandstone, grey-white fine conglomerate, glutenite and loess. Due to the larger natural moisture content, high compressibility, low strength, and poor structure of soft clay, structures built on soft clay subgrade are more prone to displacement and instability [32]. The shear strength of sandstone is linked to its lithology and water content. Rock masses with considerable water content are weak, and they show plastic failure and have low shear strength [30]. It is well-known that loess has collapsibility [33,34]; therefore, surface deformation was seen in Region P2. In Region P3, the geology is part of the Hekou Group and the Upper Pleistocene Aeolian Group, which consist of materials such as sandstone, conglomerate, glutenite, argillaceous siltstone, mudstone and loess. The shear strength of sandstone is related to its lithology and water content. The hydrological properties of sandstone are mainly swelling and disintegration after encountering water, and deformation tends to occur when the water content is brought up. Moreover, the loess is prone to collapse, so these factors will inevitably lead to regional deformation. In P4, P5, P6 and P7, the geology is part of the Upper Pleistocene Aeolian Group, which consists of loess that is prone to collapsing. The main hazard associated with collapsible loess is the damage of the large pore structure of the soil after it comes into contact with water; this results in significant collapsible deformation, causing a building’s foundation to sink [35,36,37]. In summary, the geological characteristics of Lanzhou are important factors that impact surface deformation.

4.2. Surface Deformation in Relation to Geological Faults

To verify the relationship between land surface deformation and geological faults, we drew seven profiles that crossed the fault lines. These profiles are placed in Regions P1−P7 (Figure 11).

We analyzed the surface deformation rate along seven profile lines, denoted as (1) A-A’, (2) B-B’, (3) C-C’,(4) D-D’, (5) E-E’, (6) F-F’ and (7) G-G’ (Figure 12). These lines cross the geological faults.

There are five processes for the A-A’ profile (Figure 12a):

1. The surface deformation rate is approximately −3 mm/year.

2. The surface deformation rate showed a fluctuating trend.

3. The surface deformation rate increased along the profile to about −6 mm/year.

4. The surface deformation rate decreased to about −2 mm/year.

5. The surface deformation rate rapidly increased to the maximum value of −20 mm/year.

Near the end of the profile line, the surface deformation rate showed a fluctuating upward trend before decreasing to be almost stable.

There are four processes for the B-B’ profile (Figure 12b):

1. The surface deformation rate showed a fluctuating trend.

2. The surface deformation rate increased along the profile to about −14 mm/year.

3. The surface deformation rate reduced to about 2 mm/year.

4. The surface deformation rate is increased to the maximum value of −22 mm/year.

There are three processes for the C-C’ profile (Figure 12c):

1. The surface deformation rate showed a fluctuating trend.

2. The surface deformation rate is increased to the maximum value of −30 mm/year.

3. The surface deformation rate decreased to about 5 mm/year.

There were three processes for the D-D’ profile (Figure 12d):

1. The surface deformation rate showed a fluctuating trend.

2. The surface deformation rate increased to the maximum value of −20 mm/year.

3. The surface deformation rate decreased to about −4 mm/year.

For the E-E’ profile (Figure 12e), the surface deformation rate showed a fluctuating increase tendency. For the F-F’ profile, the surface deformation rate showed an increase to about −8 mm/year and decrease to 1 mm/year in the cycle. For the G-G’ profile (Figure 12f), the surface deformation rate showed a fluctuating increase trend. It is obvious that the surface deformation rate of the A-A’, C-C’ and D-D’ profiles is the highest in the geological faults, while the surface deformation rate of the B-B’, E-E’, F-F’ and G-G’ profiles showed the maximum in faults.

The geological activity in the fault layer is actively involved and it can easily slide [38]. The location of faults has a clear effect on the spatial distribution of land subsidence, probably due to the abrupt discontinuity in horizontal hydraulic conductivity related to the structure of the basement and the vertical differential compaction of heterogeneous geological materials [39]. Hence, geological faults make a strong impact on land surface deformation. This also demonstrates that the location of the surface subsidence is closely related to its geological fault in the city of Lanzhou.

4.3. Surface Deformation in Relation to Land Cover/land Cover Changes

The overexploitation of ground water was the major cause of surface deformation [2]. However, due to the difficulty in accessing water extraction bookkeeping records from local water management agencies, we are not in a position to obtain the groundwater level changes. Instead, we utilized land cover types as a proxy for the general amount of groundwater extracted, the percentage of groundwater use in the area is: agricultural (77%), urban (13%) and industry (10%) [4]. We used the world land cover types map released by Tsinghua University, and we combined these with optical images to distinguish land cover types. The seven regions investigated in this study have different degrees of deformation and diverse land cover types.

In P1, P2, P3, P5 and P6, the land cover type is urban, which includes residences, infrastructure and a network of underground. As shown in Figure 13, that the type of land cover changed from 2015 to 2019, due to land creation (P2,P3) and construction (P5,P6). We concluded that the time of land cover changed in accordance with the time of mutation of deformation (Figure 9). Land cover change may be one reason for land surface deformation and the reason for mutation [40]. The construction of a large number of projects, including the development of underground space and the excavation of building foundation pits, resulted in the extraction of underground liquid supports, the excavation of solid supports and the destruction of the stress balance between the rock and soil, resulting in the formation of land surface deformation [41]. In P4, the cover is industry, as saw in the area known as the Xinglong Industrial Park that is home to many logistics companies and warehouses. Static load may be the main factor for deformation [42,43]. In P7, the land cover consists of a small number of houses and large stretches of farmland, which belongs to agricultural areas. In this region, new buildings were erected in 2015–2019, and the time of land cover changed in accordance with the time of the mutation (Figure 9, Figure 13). The annual amount of local rainfall is insufficient to support the region’s agricultural production, so underground water needs to be mined, which irrigated large stretches of farmland that require a large amount of water. Rain-fed and agricultural land-use types showed higher subsidence rates, due to rain-fed and irrigation had a high impact of the seasonal water table lowering [44]. To sum up, construction and irrigation may be the reasons for deformation [23,45] in region P7.

The observed surface deformation rates discussed above being summarized in

Table 3. Summary of the observations in the seven studied areas: locations, maximum deformation rate, correlation between velocity gradients and locations of pre-existing geological faults, relationship to land cover types (agricultural, industrial and urban).

Locations	Max Deformation (mm)	Correlation with Main Geology	Correlation with Faults	Correlation with Land Cover Types
P1	−115.07	Holocene Aeolian Group	Fault-limited	Urban
P2	−183.11	Xiliugou Group	Fault-limited	Urban
P3	−211.92	Hekou Group and Upper Pleistocene Aeolian Group	Fault-limited	Urban
P4	−132.07	Upper Pleistocene Aeolian Group	Fault-limited	Industry
P5	−132.54	Upper Pleistocene Aeolian Group	Fault-limited	Urban
P6	−109.04	Upper Pleistocene Aeolian Group	Fault-limited	Urban
P7	−142.20	Upper Pleistocene Aeolian Group	Fault-limited	Agriculture

. The land cover changes are illustrated in Figure 13. The time of the land cover changes are in accordance with the time of the mutation, and notarized land cover changes for land creation and construction may be the cause for the mutation. Regions P1−P7 are limited by geological faults, and there are geological differences on both sides of the fault lines, indicating that land surface deformation is pertinent to geological faults. Therefore, the surface deformation of Lanzhou is related to geology, geological faults and land cover/land cover changes.

5. Conclusions

This study was designed to detect the spatial position of surface deformation in Lanzhou and understand its mechanism using a total of 75-scene descending orbits obtained from Sentinel-1A TOPS SAR (C-band) model image data. The data were acquired in descending orbits between March 2015 and May 2019 and processed using SBAS-InSAR. Cross-validation among the PS-InSAR, SBAS-InSAR results and field investigation was performed. The relationship between surface deformation and geology, geological faults and land cover types in the city of Lanzhou was analysed in detail. The main conclusions are as follows: The results obtained by SBAS-InSAR, PS-InSAR and investigation in field work showed good agreement. The results were cross-validated to ensure the accuracy. The land surface deformation rate of Lanzhou ranges from −82.13 mm/year to 19.31 mm/year, and most of the areas are stable. However, some serious surface deformation locations were detected in the following regions: the north freight yard, Jiuzhou, Country Garden, Donggang, Yanjiaping, Zhongxinping and Liuquan. Time-series land surface deformation analysis showed that deformation increased over time in the major deformation regions. Surface deformation expansion was significant after June 2017 and the surface deformation continued to expand. Surface deformation in Lanzhou is affected by geology, geological faults and land cover/land use changes, geological structure affects the surface deformation, in that the surface deformation changes obviously across the faults. And spatial distribution of deformation is highly related to land creation and construction, industrial areas and agricultural areas. Future studies will need to go to quantitative methods to further analyze the mechanisms of surface deformation in Lanzhou city, such as the relationship between surface deformation and groundwater, the periodic characteristics of signals using more data.