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

: While surface deformations and their impact on buildings have been observed in the city of Lanzhou, it is di ﬃ cult to ﬁnd studies of surface deformation and the inﬂuential 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 signiﬁcant after June 2017and it continued to expand. The surface deformation of Lanzhou is a ﬀ ected 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 ﬁndings provide compelling data and theoretical support for disaster prevention and reduction in Lanzhou.


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. 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−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, largescale 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 0a−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 0c. 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.  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.

Study Area
Lanzhou is situated in the western part of China on a loess plateau. The urban area is about 1631.6 km 2 . 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.

Study Area
Lanzhou is situated in the western part of China on a loess plateau. The urban area is about 1631.6 km 2 . 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.

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. 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 75scene 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 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°C. 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.

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. 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].

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: Geosciences 2020, 10, 99 5 of 21 (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: where j ∈ (1, . . . , M); λ is the central wavelength of the signal; d(t B , x, r) and d(t B , x, r) are the cumulative shape of the radar line of sight for d(t 0 , x, r) = 0 at the time of t B and t A , respectively; ∆φ 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: the j th phase value of the interferogram can be expressed as: 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: Formula (5) is a matrix of M × 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.

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 si(i = 1，2，…，N). Interference pairs composed of the i and k SAR images can be recorded as , = * . Here, * stands for the complex conjugate operation. Assuming that the target p0 is the reference point, the interferometer phase∆∅ , , = ∠ , , will be dependent on the elevation, deformation, atmospheric disturbance and noise difference between the target p and the reference point p0. The phase difference caused by the elevation and line difference between the two targets can be expressed as: In the formulae, ∆ℎ , and ∆ , is the difference of elevation and deformation rate between target and . , and , represent the spatial baseline and time baseline perpendicular to the radar line of sight between SAR images i and k, respectively. λ is the radar signal wavelength, R is the oblique distance between the target and the satellite, and θ is the incident angle of SAR.
The target elevation and its subsidence rate were usually estimated by maximizing time coherence:

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 0 is the reference point, the interferometer phase ∆∅ i,k p, = ∠I i,k p,p0 will be dependent on the elevation, deformation, atmospheric disturbance and noise difference between the target p and the reference point p 0 . The phase difference caused by the elevation and line difference between the two targets can be expressed as: In the formulae, ∆h p,p0 and ∆v p,p0 is the difference of elevation and deformation rate between target p and p 0 . B i,k n and B i,k t represent the spatial baseline and time baseline perpendicular to the radar line of sight between SAR images i and k, respectively. λ is the radar signal wavelength, R is the oblique distance between the target and the satellite, and θ is the incident angle of SAR.
The target elevation and its subsidence rate were usually estimated by maximizing time coherence: where, M is the number of interferograms. ∆φ i,k p are the differential interferometer phase. ∆φ i,k H,p is the elevation correlation partial phase, and ∆φ i,k D,p is the terrain deformation correlation phase. The land deformation ∆φ i,k D,p is mainly composed of linear deformation phase and nonlinear deformation phase.
where, φ i,k non−linear is a nonlinear variable. In Formula (9), the maximization factor ξ = max ξ p was called time coherence and the time coherence of the target can be approximately expressed as 1 ∧ ξ p ≈ e −δ 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.
Geosciences 2020, 20, x FOR PEER REVIEW  7 of 22 where, M is the number of interferograms. ∆ , are the differential interferometer phase. ∆ , , is the elevation correlation partial phase, and ∆ , , is the terrain deformation correlation phase.
The land deformation ∆ , , is mainly composed of linear deformation phase and nonlinear deformation phase.
where, , is a nonlinear variable. In formula (9), the maximization factor ξ = max was called time coherence and the time coherence of the target can be approximately expressed as 1 ∧ ≈ ∅ ⁄ 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.

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: , (k= 1,2,…,n) (11)

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: where, UF k is the standard normal distribution, and the sequence of statistics is calculated by the order of time series x 1 , x 2 , ..., x n . s k is the cumulative number, E(s k ) is the mean value, and Var(S k ) is the variance. If you follow the time series x in reverse order x n , x n-1 , ... x 1 , 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.

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).
Geosciences 2020, 20, x FOR PEER REVIEW 8 of 22 Where, UFk is the standard normal distribution, and the sequence of statistics is calculated by the order of time series x1, x2 ..., xn. sk is the cumulative number, ( )is the mean value, and Var( ) is the variance.
If you follow the time series x in reverse order xn, xn-1, ... x1, the above process was repeated to calculate UBk. UFk and UBk curves are drawn, respectively. If the UFk and UBk curves intersect at a point between the borderline, then the moment at which they intersect is the abrupt point.

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 (0). To quantitatively analyze the differences between the PS-InSAR and SBAS-InSAR results, we counted the number of distribution and determined the cumulative percentage (0). 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 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). 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.

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. Geosciences 2020, 20, x FOR PEER REVIEW 10 of 22 Figure 6. The mean surface deformation rate map of Lanzhou.

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 (0) 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 (0).

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).

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.

Relationship between Land Surface Deformation and Geology
The geological data of the studied area showed five main types of geology (0): 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.
4. The Gansu Group, which consists of mudstone, sandy mudstone and glutenite mixed with marlstone.

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.

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.
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.
Geosciences 2020, 20, x FOR PEER REVIEW 14 of 22 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−37]. In summary, the geological characteristics of Lanzhou are important factors that impact surface deformation.

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).
Geosciences 2020, 20, x FOR PEER REVIEW 15 of 22 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 (0). Figure 11. Distribution of the deformation rate and geological 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.
Geosciences 2020, 20, x FOR PEER REVIEW  16 of 22 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.

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 0, 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 (0). Land cover change may be one

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. 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. 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)  Figure 13. Land cover changes in the study area. Figure 13. Land cover changes in the study area.

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.