TS-InSAR analysis for monitoring ground deformation in Lanzhou New District, the loess Plateau of China, from 2017 to 2019

Journal Pre-proofs TS-InSAR Analysis for Monitoring Ground Deformation in Lanzhou New District, the Loess Plateau of China, from 2017–2019 Yi He, Youdong Chen, Wenhui Wang, Haowen Yan PII: DOI: Reference:

S0273-1177(20)30786-9 https://doi.org/10.1016/j.asr.2020.11.004 JASR 15054

To appear in:

Advances in Space Research

Received Date: Revised Date: Accepted Date:

6 June 2020 11 October 2020 2 November 2020

Please cite this article as: He, Y., Chen, Y., Wang, W., Yan, H., TS-InSAR Analysis for Monitoring Ground Deformation in Lanzhou New District, the Loess Plateau of China, from 2017–2019, Advances in Space Research (2020), doi: https://doi.org/10.1016/j.asr.2020.11.004

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2020 COSPAR. Published by Elsevier Ltd. All rights reserved.

TS-InSAR Analysis for Monitoring Ground Deformation in Lanzhou New District, the Loess Plateau of China, from 2017–2019 Yi He1,2,3*, Youdong Chen1,2,3, Wenhui Wang1,2,3, Haowen Yan1,2,3* 1

Faculty of Geomatics, Lanzhou Jiaotong University, Lanzhou, China National-Local Joint Engineering Research Center of Technologies and Applications for National Geographic State Monitoring, Lanzhou, China 3 Gansu Provincial Engineering Laboratory for National Geographic State Monitoring, Lanzhou, China 2

* Correspondence: Yi He [email protected] Abstract Large-scale land creation projects involving the cutting of mountains to infill gullies for construction have been carried out in Lanzhou New District (LZND). However, there is an urgent need for comprehensive and detailed research on the spatiotemporal evolution of ground deformation in LZND. Based on Sentinel-1A SAR data, combined with the urban geological background, the ground deformation in LZND from 2017 to 2019 was analysed. Two independent, multi-temporal techniques, persistent scatterers interferometry (PS-InSAR) and the small baseline subset (SBAS-InSAR), were used to calculate the deformation time series, and the results were cross-verified. The time seriesmonitoring results of the SBAS and PS calculations exhibited strong consistency in LZND and verified the high reliability of the experimental results. The results showed the whole surface of the LZND from March 2017 to October 2019 maintained stability, and the deformation rate was primarily in the range of -10 to 10 mm/year. However, ground deformation in the Xicha area was evident. The maximum annual deformation rates monitored by SBAS-InSAR and PS-InSAR were -52.48 mm/year and -56.35 mm/year, respectively. The most typical deformation areas include the built-up area and the land creation area. The surface subsidence area is concentrated in the filling area. The ground deformation range of LZND kept expanding and accelerating from 2017-2019. Land creation, urban construction, geology and precipitation were the primary factors contributing to local severe ground deformation. The results of this study provide reference for the regional urban planning in LZND. Keywords: Sentinel-1A image; land creation; road network; field investigation; loess 1. Introduction As a result of population growth and the acceleration of urbanization, the functional areas of traditional cities gradually spread out, and new districts are generated (Ortega et al.,1999, Li et al., 2020). As an important way to control ‘urban disease’, new areas or new towns can reduce urban pressure and optimize the regional spatial pattern. Typical examples of this phenomenon are New York New Town, Paris New Town, Tokyo New Town and Xiongan New District (Motagh et al., 2007). In China, the national new area was a newly developed metropolitan area established in the early 1990s as a comprehensive functional area that undertakes the strategic task of national development, reform and opening up. However, with the construction of new districts and the impact 1

of human activities, different degrees of surface deformation have been detected in new districts in different regions – for example, the Xiongan New District in the Hebei province (Chaussard et al., 2013), the Pudong New District in Shanghai (Erkens et al.,2015) and the Nansha New District in the Guangdong province (Dong et al., 2018). For a long time, cities on the Loess Plateau have suffered a shortage of flat land for urban construction (Wang et al., 2018). Therefore, earthworks, such as digging mountains and filling valleys, have become an important means of obtaining flat land for urban expansion (Liu and Li, 2014). However, these processes can lead to unpredictable environmental damages and geo-hazards (Li et al.,2014, Juang et al., 2019). For example, significant uneven surface deformation exists in Yan’an New District, with the average rate being from −70 to 30 mm/year (Wu et al., 2019). Therefore, deformation caused by large-scale flat land creation needs to be monitored constantly. Currently, levelling, global navigation satellite systems (GNSS) and interferometric synthetic aperture radar (InSAR) measurements are widely used in urban surface deformation monitoring (Abidin et al., 2016, Chen et al., 2018). Levelling and GNSS can realize single-point, high-precision monitoring of local ground subsidence, but the spatial resolution is low and the cost is high (Abidin et al., 2016). InSAR technology (such as persistent scatterers interferometry (PS-InSAR) and the small baseline subset (SBASInSAR)) can obtain the surface deformation information of a large space with millimetre accuracy (Galloway and Burbey, 2011). PS-InSAR can accurately characterize slowmoving and abundant ground feature information, and SBAS-InSAR can more effectively detect surface deformation information in unstable or mountainous areas (Chaussard et al.,2014, Wasowski and Bovenga , 2014). In regions with unstable or limited radar targets, such as those with rapid urbanization and large-scale land development, the comprehensive use of PS-InSAR and SBAS-InSAR technologies can effectively describe the complex deformation process (Chen et al., 2018). In addition, the proliferation of synthetic aperture radar (SAR) satellites of varying temporal and spatial resolution offers great opportunities for researchers to combine their observations with long-term geological hazard monitoring (Armas et al., 2017). The European space agency (ESA) has provided the Sentinel-1 satellite, an earthobservation satellite, which can acquire C-band, all-weather, continuous, high-resolution SAR images at all hours of the day; this technology is of great help for the rapid monitoring of land-movement phenomena (European Space Agency, 2012). In addition, regarding the Sentinel-1 width of interference (IW) product, the stripe width is about 250 kilometers. Made using fewer images to monitor the quality of a large area, it can simultaneously maintain phase coherence and relatively high spatial resolution, and can be in 12 days or 6 days in a short period of time to visit, to monitor and quantify surface elevation changes in different environments (European Space Agency, 2012, Di Traglia ,2018). The Sentinel-1 IW product has a short revisit time and a large band width, which has been widely used to provide a very high level of support for surface deformation monitoring applications (Wu et al., 2019, Chen et al., 2018, He et al., 2020). Lanzhou New District (LZND) is a significant connecting point between the Silk Road Economic Belt and the Eurasian ‘Continental Bridge’. With the acceleration of urbanization and industrialization in LZND, the water consumption of various industries will increase greatly. The project ‘Lead into the Qin Engineering’, which transfers a large amount of water from Datong River into the Qinwangchuan Basin, is the main source of the water supply in LZND. The groundwater is a vital emergency reserve source which 2

is very important to the safety of the water supply and the sustainable development and utilization of water resources. LZND is approved by the State Council in August 2012, large-scale land creation projects involving cutting the mountains to infill the gullies to produce flat land for construction have been carried out in LZND. The large-scale earthworks and complex site conditions might induce significant surface deformation in LZND (Chen et al., 2018). There is an urgent need for comprehensive and detailed research on the spatiotemporal evolution of surface deformation after the land creation in LZND. In this manuscript, we applied SBAS-InSAR and PS-InSAR techniques to obtain time series land deformation information about LZND based on 40 Sentinel-1A (C-band) images from 2017 to 2019. The spatial and temporal characteristics were analysed. The effects of land creation, urban construction, geology and precipitation on ground deformation were discussed. The results of this study provide reference for the investigation of regional land subsidence. 2. Materials and Methods 2.1 Study Area LZND is located at 103°29′22″–103°49′56″E, 36°17′15″–36°43′29″N, which is a new national district under the jurisdiction of the Gansu province (Figure 1). LZND is in the Qinwangchuan Basin, situated in the middle of the symbiotic belt of Lanzhou (Gansu province), Xining (Qinghai province) and Yinchuan (Ningxia province). The study area is in the northwest of the Longxi Loess Plateau. It is the confluence of the Qinghai-Tibet Plateau, the Mongolian Plateau and the Loess Plateau, and it is a hilly landform typical of the Loess Plateau. LZND is developed by the geomorphology of the Pingchuan, Liangmao gully and river valley. The mountains around the basin are mainly covered by late-Pleistocene loose aeolian loess, which is the main material used to fill valleys and form flat land for construction. LZND is the fifth national new area in China and the first national new district in northwest China approved by the State Council in August 2012. Its development planning period is 2011-2030, the period from 2013 to 2017 is the main construction period. However, the temporal and spatial characteristics of surface deformation after construction are still unknown. The study area has a typical, semi-arid, continental climate – it is cold and dry in winter and windy and rainy in spring. The annual average temperature and precipitation are 6.9 °C and 300–350 mm, respectively. With the acceleration of urbanization and industrialization in LZND, the change of land use types is very distinct.

3

Figure1. Location and topographic characteristics of the study area 2.2 Data In study area and study period (2017-2019), there are a total of 56 scene Sentinel1A images (C-band, wavelength 5.6 cm) of the Terrain Observation with Progressive Scans (TOPS) mode of the same descending orbit provided by the European Space Agency. But the data from March 2017 to March 2018 are only 13 images, it is basically one image per month. The existing image data is relatively rough, while the image data is relatively dense in the later period (March 2018 to October 2019). To study the uniformity of the image data during the entire time period, 1-2 scene images are guaranteed in per month. Finally, the data of 40 scene descending Sentinal-1A orbit are used in this paper. The incident angle is 39.5°, the spatial resolution is 5×20 m (azimuthal × range direction), the observation mode is interference width (IW), and the polarization mode is single polarization (VV). The data collection period was from 20 March 2017 to 24 October 2019, as shown in Table 1. The 30 m resolution DEM (Digital Elevation Model) from Shuttle Radar Topography Mission (SRTM) data published by USGS was used to remove the topographic phase and encode it. In order to achieve build-up data, two Gaofen-1 remote sensing images from 2017 and 2019 can be downloaded from China Centre for Resources Satellite Data and Application (http://218.247.138.119:7777/DSSPlatform/index.html). Table 1. Sentinel-1A data for Lanzhou New District NO. 1 2 3

Acquisition Date 20 March 2017 25 April 2017 19 May 2017

NO. 21 22 23

4

Acquisition Date 23 September 2018 5 October 2018 29 October 2018

4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

24 June 2017 23 August 2017 16 September 2017 22 October 2017 21 December 2017 2 January 2018 7 February 2018 19 February 2018 8 April 2018 26 May 2018 7 June 2018 19 June 2018 1 July 2018 25 July 2018 6 August 2018 30 August 2018 11 September 2018

24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40

10 November 2018 22 November 2018 4 December 2018 28 December 2018 9 January 2019 21 January 2019 2 February 2019 26 February 2019 22 March 2019 3 April 2019 27 April 2019 9 May 2019 21 May 2019 14 June 2019 20 July 2019 18 September 2019 24 October 2019

2.3 Methods 2.3.1 SBAS-InSAR The SBAS-InSAR technique is a kind of time series InSAR (TS-InSAR) that uses multiple master images to complete interferometry (Berardino et al., 2002). The principle is an analysis of the phase-coherent target (CT) in order to acquire the timing deformation. In a synthetic aperture radar image, the value for each pixel is the coherent sum of contributions from all scatterers within the associated ground resolution element (Hooper 2008). However, changes in the motion and angle of the scatterers cause decorrelation (Zebker and Villasenor, 1992, Zhang et al., 2020). Decorrelation is reduced by spectral filtering in range and discarding the non-overlapping doppler frequencies in azimuth (Hooper, 2008). Furthermore, the synthetic aperture radar image also includes multilooking interferograms and pixels whose phases, when thus filtered, decorrelate little over short time intervals. These pixels, which are referred to as slowly-decorrelating filtered phase (SDFP) pixels, are the targets of the SBAS-InSAR method. There are four main steps in SBAS-InSAR processing using SARscape 5.5. The detailed process is shown in Figure 2. Connection graph generation

Super master image selection

Interference process

SBAS inversion

Geocoding

Registration

Remove orbital error

Dummy removal

Coherent graph generation

Remove residual atmospheric phase

Relax interpolation

De-flattening filtering

SVD

Geocoding

Phase unwinding

Estimate deformation rate

Time-series deformation results

Baseline estimation

Figure2. SBAS-InSAR basic process framework

5

The first step is the generation of the connection graph, including automatic super master image selection and baseline estimation. The 19 June 2018 image is automatically selected as the super master image, and the time baseline is 200 days (Figure 3). The second step is the interference process. A total 140 differential interferograms were generated. Multi-looking with 4×1 in the range and in the direction of the azimuth and the Goldstein filtering method (Zebker and Villasenor, 1992) were used for improving the signal-to-noise ratio (SNR) of the interferograms. The progress phase unwrapping of the wrapped phase signal is performed with a three-dimensional (3D) phase unwrapping approach (Hooper and Zebker, 2007). Some interferogram pairs with low coherence and unwrapped phase errors were discarded (Chen et al., 2013). The third step is the SBAS-InSAR inversion process. The remaining phase content and phase ramps were calculated by selecting and refining the stable ground control points (GCPs). We selected GCPs on the basis of their deformation being close to zero, confirmed by the interpretation of the unwrapped phase (Ghulam et al., 2015). As a consequence, 40 GCPs were applied to refine the orbital error and finish the deformation inversion steps. Then, the initial deformation were estimated by the linear model, and the residual terrain was removed. Finally, the Singular Value Decomposition (SVD) method was used to retrieve the time series of land surface deformation by combining multiple small baseline data. After the average deformation velocity was estimated, the influence of the atmosphere could be removed by high-pass filtering in the time domain and lowpass filtering in the spatial domain. In order to improve the monitoring accuracy, iterative deformation estimation and atmospheric effect removal were carried out. The fourth step is the geocoding process. The first and second SBAS inversion results were geocoded. The dummy value was removed and the deformation results were interpolated. InSAR-derived deformation are in the line-of-sight (LOS) direction, because the surface after land creation in the study area is relatively flat, and the vertical deformations dominate ground deformation, therefore, the deformation rate of LOS dLOS converts to the vertical direction dv (He et al., 2020). According to the following incidence angle θ: dLOS

dv = cosθ

(1)

Figure3. A network of interferogram pairs obtained from Sentinel-1A using SBASInSAR 6

2.3.2 PS-InSAR The PS-InSAR theory and technique was proposed by Ferretti et al (2001). The basic principle is extracted after the differential interference processing of the PS point in phase and then deducted from the original differential interference figure linear deformation and elevation error component to get the residual phase diagram. Atmospheric delay and phase error of orbit were removed, the deformation linear and nonlinear deformation superposition of deformation time series is obtained. The biggest distinctions between this technique and SBAS-InSAR technique are that the PS technique targets resolution cells dominated by a single scatterer, while the former targets cells with many scatterers, none of which dominate. There are four main steps in PS-InSAR processing using SARscape 5.5. The specific process is shown in Figure 4. Connection graph generation

Super master image selection

Interference process

PS inversion

Geocoding

Registration

Remove atmospheric phase

Deformation rate estimation

Coherent graph generation

Detection of PS

Precision estimation

PS point selection

Residual decomposition

Geocoding

De-flattening

Deformation rate

Time-series deformation results

Baseline estimation

Figure4. PS-InSAR basic process framework The 19 June 2018 image is selected as the super main image, and the time baseline is 300 days (Figure 5). The first inversion uses a linear model to derive the residual height and the deformation velocity. A sub-region is 25 km2. Each sub-region was processed separately by selecting a reference point based on the amplitude dispersion index. The results of all sub-regions are mosaiced to output a complete result. The second inversion uses the products of the first linear model from the previous step to estimate the atmospheric phase components. The atmospheric filter uses a temporal high-pass filter (365 days) and a spatial low-pass filter (1200 m) to remove the atmospheric phase components. Finally, the deformation rate and time series deformation are obtained. The final inversion results are geocoded. Ultimately, the deformation rate of the line-of-sight (LOS) dLOS converts to the vertical direction dv.

Figure 5. A network of interferogram pairs obtained from Sentinel-1A using PS-InSAR

7

2.3.3 The theoretical precision evaluation In the process of SBAS-InSAR and PS-InSAR data processing, the principal errors occur in the domains of data orbit accuracy, coherence, atmospheric phase error and phase unwrapping. Therefore, the theoretical accuracy of the data processing process need be evaluated. The interference measurement process is the core of the SBAS-InSAR and PSInSAR technologies. The theoretical accuracy of the experimental results was evaluated by the phase variance (Ramon Hanssen, 2001). The higher the phase variance value is, the lower the experiment precision was. The formula used for the precision calculation is: σ2φ =

1 ― γ2 2γ2

(2)

[rad2]

Where σ2φ is the phase variance and γ is the interferometric coherence coefficient. The interferometric coherence coefficient (γ) is calculated as follows: N

γ=

M

|∑i = 1∑j = 1μ1(i,j)μ2(i,j)| N

M

2 N

M

∑i = 1∑j = 1|μ1(i,j)| ∑i = 1∑j = 1|μ2(i,j)|

2

(3)

M and N are the size of the data block to calculate the coherence, i and j are the row and column numbers in the data block, μ 1 and μ 2 are the complex values at the image coordinates (i, j) in the main pair image data block,  denotes conjugate multiplication. 3. Result 3.1 Comparison of ground deformation results between SBAS-InSAR and PS-InSAR The accuracy comparison of ground deformation results between SBAS-InSAR and PS-InSAR includes the theoretical accuracy of the data processing process, the comparison of the results and the field investigation. Field validation also is discussed in a subsequent analysis. 3.1.1 The theoretical accuracy of the data processing process The external DEM data was used for the registration of the master image and slave image. Then, combination images from each pair were used for the differential interferometry, and the interference graph and correlation coefficient graph were obtained. The external DEM data were removed for the terrain phase, and the Goldstein filtering method was used to eliminate noise. In the phase progress, the unwrapping signal was performed with a 3D phase unwrapping approach. Finally, the errors in data processing were evaluated according to the variables of interference coherence (γ) and wavelength (λ). The interference measurement process is the core of the SBAS-InSAR and PS-InSAR technologies. Surface variations can easily cause loss of coherence that corrupts interferometric signals (Zebker and Villasenor, 1992; Dai et al.,2018). A coherence analysis was performed to check whether the conventional time series InSAR (PS-InSAR, SBAS-InSAR) can be performed with the Sentinel-1 images described above in this study. Figure 6 shows the temporal coherence map of some interferometric pairs based on SBAS-InSAR and Figure 7 shows the coherence coefficient based on PSInSAR, which formed by the 40 Sentinel-1 images and their mean coherence coefficients. 8

The interferograms with coherence higher than 0.5 are regarded as high-quality ones to easily perform correct phase unwrapping (Dai et al.,2018). The overall coherence coefficients of all the interferograms is very higher than 0.5, therefore, the interferograms in the study area have the better coherence. The interferograms are enough to form a robust spatial-temporal network for calculation and successful subsidence measurement. The phase variance σ2φ was calculated for the theoretical accuracies. The SBASInSAR and PS-InSAR achieved favorable coverage in the study area, with the phase variance of 0-0.31 and 0-0.03, respectively. The results showed that the theoretical accuracies of the two technologies were in good agreement, and the calculated phase variance σ2φ was small for both, indicating that the experimental data processing parameters were designed reasonably and the theoretical accuracy was high.

Figure 6. Coherence map of some interferometric pairs based on SBAS-InSAR

Figure 7. Coherence coefficient based on PS-InSAR 9

3.1.2 Comparison of ground deformation results LOS directional variables were obtained by the two sequential processing methods of SBAS-InSAR and PS-InSAR (Figure 8 a, b). Due to the lack of data of the measured level of LZND in the same period, this paper compared the monitoring results of SBASInSAR and PS-InSAR, and the field validation was performed. From the perspective of the spatial distribution (Figure 8 a, b), the results of SBAS-InSAR and PS-InSAR were consistent in most areas. In order to better compare the results of the two technologies, three typical regions (regions A, B and C) were selected (Figure 8). The ArcGIS software was used to extract the deformation values of SBAS-InSAR and PS-InSAR corresponding to the two InSAR values of the same location, and a pearson correlation analysis was carried out (Figure 9). The correlation coefficients of regions A, B and C were 0.79, 0.72, 0.82, respectively, which indicated that the two time-series monitoring results had a strong consistency, and also verified the high reliability of the experimental results. The research team conducted a field investigation on the typical deformation areas in Lanzhou New District and found that the InSAR measurement results were consistent with the field investigation, indicating that the experimental results of this study were credible (Figure 8). However, there were some differences in the spatial distribution between SBASInSAR and PS-InSAR, which were mainly due to the two different implementation methods, especially the different selection processes of the connection graphs, which may have led to the deviation of the deformation rate values. The SBAS-InSAR retains discretetargets, and the backscatter response of radar echoes is formed by the backscattered signals of all targets entering the resolution unit, the SBAS-InSAR method obtains the deformation result of a surface, which is continuously in spatial. In contrast, the PS-InSAR method only considers the response of each individual scatterer in the time interval, which has a higher spatial resolution and can better describe the details of surface deformation. As mentioned above, the difference in spatial distribution is due to the difference of two methods. PS-InSAR is more capable of describing the details of surface deformation than SBAS-InSAR. The surface deformation details of PS-InSAR are more obvious in space: the deformation map of PS-InSAR is spatially discontinuous (the surface deformation has jump in spatial distribution) and the point density of PS-InSAR in mountainous areas is less than the point density of SBAS-InSAR. As the coherence of the SBAS-InSAR technology increases, so does the space coverage. This is especially true in areas where rapid sedimentation occurs – the PS-InSAR method is unable to ‘follow’ this trend, as it loses potential candidates for persistent scattering which resulted SBAS-InSAR is better than PS-InSAR in low coherence area, the most typical area is Xicha, the deformation location is roughly equivalent, but the coherence of SBAS-InSAR is better than PS-InSAR and decoherence area is less. This may be because the study area belongs to a mountainous area, where there are fewer PS scattering body points, resulting in a larger error than that of SBAS. SBAS-InSAR can more effectively detect surface deformation information in unstable or mountainous areas (Chaussard et al., 2014, Wasowski and Wasowski, 2014). In the study area, the deformation rate of SBAS-InSAR inversion was -52.48–31.13 mm/year, and the inversion PS-InSAR diagram showed that the deformation rate was 56.35–50.68 mm/year. The deformation rate difference is caused by different GCP points selected in PS-InSAR and SBAS-InSAR. In the course of the experiment, it was found that the GCP points of PS-InSAR and SBAS-InSAR were not universal, so we chose different GCP points based on the same criteria in the two experiments. The results obtained by the SBAS-InSAR method are more stable than those obtained by the PSInSAR method. SBAS-InSAR smooths the deformation during the phase unwrapping 10

process and the jump of surface deformation is smaller than that of a single point, which caused the results of SBAS-InSAR method are more stable than those obtained by the PS-InSAR method (Figure 9).

Figure 8. Ground deformation results between SBAS-InSAR and PS-InSAR from 2017–2019. (a) SBAS-InSAR, (b) PS-InSAR

11

Figure 9. Cross-comparison analysis of SBAS-InSAR and PS-InSAR for the A, B, and C in the same location as shown in Figure 7. (a) and (b) for A, (c) and (d) for B, (e) and (f) for C 3.2 Spatial-temporal characteristics of ground deformation in LZND Figure 8 (a) and (b) show that the deformation trends in LZND obtained by the methods of SBAS-InSAR and PS-InSAR from 2017 to 2019 were basically similar. The whole surface of the LZND from March 2017 to October 2019 maintained stability, and the deformation rate was primarily in the range of -10 to 10 mm/year. However, ground deformation in the Xicha area was evident. The maximum annual settlement rates monitored by SBAS-InSAR and PS-InSAR were 52.48 mm/year and 60.52mm/year, respectively. There was a land uplifting area in the residential area on the southeast of Zhongchuan. The spatial distribution of two monitoring results has a high consistency. The most typical areas include the built-up area (P1, P2) and the land-creation area (P3, P4). The time accumulation sequence diagram also shows that the spatial consistency of the results of SBAS-InSAR and PS-InSAR was relatively high. From 2017 to 2019, the cumulative surface deformation range in Xicha area of LZND kept expanding and accelerating (Figure 10). The major land creation projects have been basically completed, and the surface deformation experienced a two-year uneven subsidence, with an increasing rate from 2017 to 2019.

12

Figure 10. The time series deformation based SBAS-InSAR (a, b, c, d) and PSInSAR (e, f, g, h). 3.3 Analysis of deformation characteristics of typical regions There was serious non-uniform deformation in Qinchuan, Zhongchuan and Xicha areas of LZND, and the local differential deformation will produce tensile stress on the ground, which may cause secondary disasters, such as exposed building foundations, cracking walls and pipeline ruptures (Wu et al., 2019). The four typical regions of uneven subsidence were selected (P1, P2, P3 and P4 in Figure 8) and analysed based on the results of SBAS-InSAR. 3.3.1 P1 region For the P1 region, the average deformation was -13.67–11.15 mm/year from 2017 to 2019 (Figure 11 a). The most severe deformation was close to the mountains, where the main earth filling operation and the building was carried out (Figure 11 b, c). Moreover, Figure 11d shows that the cumulative deformation increased annually. It is found that the sinking reached to 10 mm from June 2017 to October 2017, which might be attributed to the packed loess characterized by weak structure, loose arrangement of soil grains and poor resistance to external load in filling area. The uplift was found from October 2017 to April 2018, which might be attributed to the landfill in building reconstruction. There was a sinking from April 2018 to July 2018, the reason may be that this stage had abundant precipitation, the soil in the loess area was loose and precipitation can easily cause subsidence. In the period from July 2018 to May 2019, the curve was stable. A significant sinking of 25 mm was observed from May 2019 to October 2019, which might be attributed to the building construction. Furthermore, many residential and commercial buildings were built along the deformation region. Thus, land creation and urban construction may be the primary human factors contributing to local severe deformation.

13

Figure 11. (a) The deformation region in LZND superimposed on the Google Earth image; (b, c) the change from 2017 to 2019 superimposed on the Google Earth image; (d) the cumulative deformation of the P1 region. 3.3.2 P2 region

Figure 12. (a) The deformation region in LZND superimposed on the Google Earth image; (b, c) the ground sinks and the walls crack; (d) the evolution in the P2 region; (e, f, g) the surface change based on the Google Earth image. 14

As can be seen from Figure 12, the P2 region presented subsidence (Figure 12 a), with an annual rate of change of -13.34–12.44 mm/year from 2017 to 2019. The field survey results also confirmed that there was a serious subsidence disaster in this region: the ground was sinking and the walls were cracked (Figure 12 b, c). Figure 12 d showed that since 20 March 2017, the settlement had been relatively fast. From October 2017 to February 2019, the settlement fluctuated, and then showed an upward trend, followed by a downward trend in August 2019. As can be seen from the image (Figure 12 e, f, g), the P2 region belonged to the cultivated land area in 2016. On 20 March 2017, P2 had been developed into the industrial base of Lanzhou railway port, which was continuously expanded from 2017 to 2019. Therefore, the high traffic density, vibration associated with large trucks, static load and construction may be the main factors of deformation in this region (He et al., 2020, Zhou et al., 2017, Dixon et al., 2006). 3.3.3 P3 region

Figure 13. (a) The deformation region in LZND superimposed on the Google Earth image; (b, c) the ground sinks and land creation; (d) the evolution in the P3 region; (e, f, g) the surface change based on the Google Earth image. It can be seen from Figure 13 that the settlement of the P3 region was obvious (Figure 13 a), and the subsidence rate from 2017 to 2019 was -15.14–23.75mm/year. Field investigation results also confirmed the existence of serious subsidence hazards in the area (Figure 13 b, c). As can be seen in Figure 13 (d), the cumulative subsidence showed an increasing trend. As can be seen from Figures 13 (e), (f) and (g), a large area of mountain was cut from 2013 to 2017, and construction land increased from 2017 to 2019. Land cover types had changed significantly from 2013 to 2019. The change of land cover 15

was one of the causes of surface deformation. A large number of construction projects, including the development of underground space and the excavation of building foundation pits, led to the underground extraction of liquid support and solid support and the destruction of the pressure balance between rock and soil, leading to ground deformation (Qin et al., 2017). 3.3.4 P4 region From 2013 to 2019, this region experienced a rapid period of construction. The subsidence was serious, and the maximum subsidence rate was 52.48 mm/year (Figure 14 a). The field investigation also confirmed that the ground beneath them in the area fell and that there were cracks in the ground (Figure 14 b). At the same time, some areas were under construction (Figure 14 c). The cumulative subsidence showed the increasing trend (Figure 14 d). As can be seen from the Figures 14 (e), (f) and (g), a large area of the mountain was cut to infill gullies, and construction land increased from 2017 to 2019. Generally speaking, when the loess is in the dry state, the compression deformation under the external load is an important part of the loess deformation. On the other hand, when the loess is wet, the loess soil layer is prone to collapsible deformation under the action of gravity, and the deformation degree gradually decreases over a long period of time (Wu et al., 2019). In addition, a large number of buildings added a large external load to the surface, so the process of ground deformation was accelerated (Chen et al., 2018).

Figure 14. (a) The deformation region in LZND superimposed on the Google Earth image; (b, c) the ground sinks and the land creation; (d) the evolution in the P4 region; (e, f, g) the surface change based on the Google Earth image. 16

4. Discussion 4.1 Influence of geology on ground deformation The geological data of the studied area included the Upper Pleistocene Aeolian Group and the Holocene Aeolian Group as the two main types of geology (Figure 15). As can be seen from Figure 15, the spatial distribution characteristics of surface deformation in LZND are controlled by lithology. The Upper Pleistocene Aeolian Group mainly consists of loess, and the Holocene Aeolian Group mainly consists of sub-sand soil, a sub-clay sand layer and a gritty layer lens body. It is well known that loess has collapsibility (Wu et al., 2019, Yuan and Wang, 2009). During the period of 2017–2019, large-scale land development and rapid urbanization in Xicha town led to obvious surface deformation in the areas covered by these soft sediments. The land use destroyed the original loess structure and caused serious instability in the loess foundation (Jiang et al., 2014). In addition, the large number of buildings increased the external pressure on the surface, accelerating the compressible deformation of the loess soil and building foundations (Wu et al., 2019, Cui, 2018). The subsidence distribution was also influenced by the presence and thickness of compressible sediments (Chen et al., 2017). Quaternary loose sediments are the main compressible sediments in LZND. The thickness of compressible sediments controls the spatial distribution of land subsidence. In general, the greater the thickness of soft soil, the higher the subsidence rate (Juang et al., 2019). There is a high correlation between surface deformation and earthwork volume (Wu et al., 2019). This phenomenon is well explained by the large-scale land creation projects caused by mountain cutting and valley filling. A large number of earthworks have destroyed the natural state and local geological structure of the original loess, produced a large number of loess filling areas, and changed the original landscape and hydrological conditions (Wu et al., 2019).

Figure 15. The relationship between geology and ground deformation 17

4.2 Influence of land creation and build-up on ground deformation Figure 16 shows the relationships land creation and ground deformation. It can be clearly seen from the Figure 16 a that the surface of the excavation area appears to be uplifted. As mountain excavation removes the upper part of the mountain, and the remaining, lower part of the mountain releases the rebound stress in the process of removal. This phenomenon gradually leads to the uplift of the ground surface in the hilly excavation area until the equilibrium state is reached. In addition, Figure 16 b shows that the surface subsidence area is concentrated in the filling area. As gully filling changes the natural state of the original loess. The surface of the original loess layer in the gully region is covered with the filled loess layer (Li et al, 2014, Wu et al., 2019). The packed loess is characterized by weak structure, loose arrangement of soil grains and poor resistance to external load, resulting in poor engineering performance. The natural development of a layer of loess is a time-consuming process, including compression, foldable deformation, integration and compaction (Guo, 1988). Therefore, the change of surface structure and subsidence are not completely stable. In addition, the dynamic compaction method can only weaken the collapsibility and improve the strength of the loess filling to a certain extent (Li et al., 2014). It cannot completely change the physical characteristics of collapsible loess, resulting in incomplete compaction of the loess foundation, which in turn makes land subsidence more likely to occur. Moreover, due to the collapsibility of loess, the structure of the soil will be destroyed and the strength will be reduced rapidly when it meets water. The compaction degree of the loess in the filled area is far less than that of the original loess in the natural state, which leads to the loess in the filled area is more likely to settle. In addition, there are serious stability problems in unnatural connections.

Figure 16. The relationships between land creation and ground deformation. (a) Excavation area, (b) Filled area

18

In the past six years, LZND has experienced rapid urbanization. Earthwork began in 2013, and the construction of major roads, buildings and infrastructure in the new area was basically completed in 2017. According to the survey on the change of floor area in LZND from 2013 to 2019, there is a high correlation between the subsidence phenomenon and the construction (Wu et al., 2019). From 2013 to 2017, the study area experienced a rapid construction period with severe subsidence. Figure 17 shows the relationship between build-up and ground deformation. It indicates that some uplifts and subsidence existed in the build-up area. This phenomenon might due to the fact that many buildings were rebuilt and the buildings added a lot of load over the ground. Therefore, build-up has a strong influence on ground deformation.

Figure 17. The relationship between build-up and ground deformation 4.3 Influence of precipitation on surface deformation Precipitation is one of the inducing factors for widespread land deformation in the loess area (Zhang, 2016). The soil in the loess area is loose and precipitation can easily cause collapse. Surface water infiltration into the loess foundation damaged its internal structure, causing the collapse of subgrade and building foundations, which in turn led to the deformation of the foundation surface (Wu et al., 2019). The filling areas in LZND consists of loess; due to the low density and content of filling loess, the structural strength of soil particles is poor (Wang, 2018). As can be seen in Figure 18, the cumulative precipitation has a great relationship with cumulative surface deformation. The greater the cumulative precipitation, the greater the cumulative deformation, pearson correction is 0.84. Rain water reduces the cohesive strength and matric suction of loess clay, resulting in increased softening and collapsibility of the loess. Therefore, precipitation is also a major factor affecting surface deformation in LZND.

19

Figure18. The relationship between cumulative precipitation and cumulative surface deformation. 5. Conclusions Based on Sentinel-1A SAR data, combined with the urban geological background, ground deformation in LZND from 2017 to 2019 was analysed. Two independent multitemporal techniques, SBAS-InSAR and PS-InSAR, were used to calculate the deformation time series, and the results were cross-verified. The conclusions are as follows: (1)The two time series-monitoring results of SBAS-InSAR and PS-InSAR had a strong consistency in LZND and verified the high reliability of the experimental results. (2)There is serious, non-uniform deformation in LZND. Ground deformation in the Xicha area was evident. The maximum annual subsidence rates monitored by SBASInSAR and PS-InSAR were 52.48 mm/year and 50.68 mm/year, respectively. The most typical areas include the built-up area and the land creation area. From 2017 to 2019, the ground deformation range of LZND kept expanding and accelerating. (3)Land creation, build-up, geology and precipitation may be the primary human factors and natural factors contributing to local severe deformation in LZND. (4)Although there are limitations of the level points in this study, the ground deformation characteristics of LZND were evaluated by combining SBAS, PS technology and field research, which provided support for the effective planning of land use and the mitigation of relevant risks. Setting up level point validation in the future will help validate and improve these initial results. 6. Conflict of Interest The authors declare no conflict of interest. 7.Acknowledgments The authors would like to thank the European space agency and the China Centre for Resources Satellite Data for providing data. 8.Funding This work is funded by the National Key R&D Program of China (2017YFB0504201), The Youth fund of LZJTU, No.2017002; LZJTU EP, No201806, Tianyou Youth talent lift program of Lanzhou Jiaotong University, the National Scientific Foundation of China (Grants Nos. 41371435, 41761088, 41761014, 41801395, 41761082 and 41861059) and the National Key Research and Development Program of China (2016YFC0803106), "InSAR monitoring of land subsidence in main urban area of Lanzhou city" by department of education of Gansu Province(2019A-043); China 20

postdoctoral science foundation (2019M660092XB), Lanzhou Jiaotong UniversityTianjin University Innovation Project Fund Project (2020055), Gansu Provincial Department of Transportation Project (2020-11). 9.References Abidin, H.Z., Andreas, H., Gumilar, I., Yuwono, B.D., Murdohardono, D., Supriyadi, S., 2016. On Integration of Geodetic Observation Results for Assessment of Land Subsidence Hazard Risk in Urban Areas of Indonesia. In IAG 150 Years, Springer: Cham, Switzerland. 143, 435–442. Armas, I., Mendes, D.A., Popa, R.G., Gheorghe, M., Popovici, D., 2017. Long-term ground deformation patterns of Bucharest using multi-temporal InSAR and multivariate dynamic analyses: A possible transpressional system. Scientific Reports. 7, 43762. Berardino, P., Fornaro, G., Lanari, R., Sansosti, E., 2002. A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms. IEEE Transactions on geoscience and remote sensing. 40(11): 2375-2383. Chaussard, E., Wdowinski, S., Cabral-Cano, E., Amelung, F., 2014. Land subsidence in central Mexico detected by ALOS InSAR time-series. Remote Sensing of Environment. 140(1):94-106. Chaussard, E.Amelung, F., Abidin, H., Hong, S., 2013. Sinking cities in Indonesia: ALOS PALSAR detects rapid subsidence due to groundwater and gas extraction. Remote Sensing of Environment. 128, 150-161. Chen, B., Gong, H., Li, X., Lei, K., Zhu, L., Gao, M., Zhou, C., 2017. Characterization and causes of land subsidence in Beijing, China. International Journal of Remote Sensing. 38(3): 808–826. Chen, F., Lin, H., Zhou, W., Hong, T., Wang, G., 2013. Surface deformation detected by ALOS PALSAR small baseline SAR interferometry over permafrost environment of Beiluhe section, Tibet Plateau, China. Remote Sensing of Environment. 138, 10-18. Chen, G, Zhang, Y., Zeng, R., Yang, Z., Chen, X., Zhao, F., Meng, X., 2018. Detection of land subsidence associated with land creation and rapid urbanization in the chinese loess plateau using time series insar: A case study of Lanzhou new district. Remote Sensing. 10(2): 270. Cui, Z., D., 2018. Physical Model Test of Layered Soil Subsidence Considering Dual Effects of Building Load and Groundwater Withdrawal. In Land Subsidence Induced by the Engineering-Environmental Effect. Springer: Singapore, pp. 169–191. ISBN 978981-10-8040-1. Dai, K, R., Liu, G. X., Li, Z. H., et al., 2018. Monitoring Highway Stability in Permafrost Regions with X-band Temporary Scatterers Stacking InSAR. Sensors. 18(1876):1-17.

21

Di, Traglia, F., Nolesini, T., Ciampalini, A., Solari, L., Frodella, W., Bellotti, F., Fumagalli, A., De Rosa, G., Casagli, N., 2018. Tracking morphological changes and slope instability using spaceborne and ground-based SAR data. Geomorphology. 300, 95–112. Dixon, T.H., Amelung, F., Ferretti, A., Novali, F., Rocca, F., Dokka, R., Sella, G., Kim, S-W., Wdowinski, S., Whitman, D., 2006. Space geodesy: Subsidence and flooding in New Orleans. Nature. 441(7093): 587–588. Dong, S., Samsonov, S., Yin, H., Huang, L., 2018. Two-Dimensional Ground Deformation Monitoring in Shanghai Based on SBAS and MSBAS InSAR Methods. Journal of Earth Science. 29(4): 960-968. Erkens, G., Bucx, T., Dam, R., Lange, D., Lambert, J., 2015. Sinking coastal cities. Proceedings of the International Association of Hydrological Sciences. 372, 189-198. European Space Agency (ESA)., 2005. Mission Requirements Document for the European Radar Observatory Sentinel-1, ES-RS-ESA-SY-0007. European Space Agency: Paris, France. 1, 1–31. Attema, E., Bertoni, R., Bibby, D., Carbone, A., Cosimo, G.D., Geudtner, D., Giulicchi, L., Lokas, S., Navas-traver, I. and Ostergaard, A., 2012. Sentinel-1: ESA's Radar Observatory Mission for GMES Operational Services. In ESA Communications report ESA SP-1322/1. European Space Agency. Ferretti, A., Prati, C., Rocca, F., 2001. Permanent scatterers in SAR interferometry. IEEE Transactions on geoscience and remote sensing. 39, 8-20. Galloway, D.L., Burbey, T.J., 2011. Regional land subsidence accompanying groundwater extraction. Hydrogeol. 19(8), 1459–1486 Ghulam A, Grzovic M, Maimaitijiang M, Sawut, M., 2015. InSAR monitoring of land subsidence for sustainable urban planning. CRC Press: Boca Raton, FL, USA,. Guorui, G., 1988. Formation and development of the structure of collapsing loess in China. Engineering Geology. 25(2-4): 235–245. Hanssen, R.F., 2001. Radar interferometry: data interpretation and error analysis (Vol. 2). Springer Science & Business Media. ISBN: 0-7923-6945-9. He, Y., Wang, W., Yan, H., Zhang, L., Chen, Y., Yang, S., 2020. Characteristics of Surface Deformation in Lanzhou with Sentinel-1A TOPS. Geosciences. 10(3): 99. Hooper, A., Zebker, H, A., 2007. Phase unwrapping in three dimensions with application to InSAR time series. JOSA A. 24(9): 2737-2747. Hooper, A., 2008. A multi‐temporal InSAR method incorporating both persistent scatterer and small baseline approaches. Geophysical Research Letters. 35(16):1-5. Jiang, M., Zhang, F., Hu, H., Cui, Y., Peng, J., 2014. Structural characterization of natural loess and remolded loess under triaxial tests. Engineering Geology. 181, 249–260.

22

Juang, C.H., Dijkstra, T., Wasowski, J., Meng, X., 2019. Loess geohazards research in China: Advances and challenges for mega engineering projects. Engineering Geology. 251, 1–10. Li, P., Qian, H., Wu, J., 2014. Environment: Accelerate research on land creation. Nature. 510(7503): 29–31. Li, Y., Sun, Q., Ji, X., Xu, L., Lu, C., Zhao, Y., 2020. Defining the Boundaries of Urban Built-up Area Based on Taxi Trajectories: a Case Study of Beijing. Journal of Geovisualization & Spatial Analysis. 4(1):8. Liu, Y., Li, Y., 2014. Environment: China’s land creation project stands firm. Nature. 511(7510): 410. Minderhoud, P.S.J., Coumou, L., Erban, L.E., Middelkoop, H., Stouthamer, E., Addink, E.A., 2018.The relation between land use and subsidence in the Vietnamese Mekong delta. The Science of Total Environment. 634, 715–726. Motagh, M., Djamour, Y., Walter, T, R., Wetzel, H., Zschau, J., Arabi, S., 2007. Land subsidence in Mashhad Valley, northeast Iran: results from InSAR, levelling and GPS. Geophysical Journal International. 168(2): 518-526. Ortega‐Guerrero, A., Rudolph, D, L., Cherry, J, A., 1999. Analysis of long‐term land subsidence near Mexico City: Field investigations and predictive modeling. Water Resources Research. 35(11): 3327-3341. Pei, X., Zhang, F., Wu, W., Liang, S., 2015. Physicochemical and index properties of loess stabilized with lime and fly ash piles. Applied Clay Science. 114, 77–84. Qin, X, Q., Yang, M, S., Liao, M, S., Wang, H, M., Yang, T, L., 2017. Exploring Temporal-Spatial Characteristics of Shanghai Road Networks Settlement with Multitemporal PSInSAR Tecnique. Geomat. Geomatics and Information Science of Wuhan University. 42, 170–177. Wang, F., 2018. Experimental Study on Deformation Characteristics of Structural Loess under High Stress in Yan’an Area. Master’s Thesis, Xi’an University of Science and Technology, Xi’an, China. Wang, S., Fu, B., Chen, H., Liu, Y., 2018. Regional development boundary of China’s Loess Plateau: Water limit and land shortage. Land Use Policy.74, 130–136. Wasowski, J., Bovenga, F., 2014. Investigating landslides and unstable slopes with satellite multi temporal interferometry, current issues and future perspectives. Engineering Geology. 174(8):103-138. Wu, Q., Jia, C., Chen, S., Li, H., 2019. SBAS-InSAR Based Deformation Detection of Urban Land, Created from Mega-Scale Mountain Excavating and Valley Filling in the Loess Plateau: The Case Study of Yan’an City. Remote Sensing. 11(14), 1673. Yuan, Z, X., Wang, L, M., 2009. Collapsibility and seismic settlement of loess. Engineering Geology. 105(1-2): 119–123. 23

Zebker, H, A., Villasenor, J., 1992. Decorrelation in interferometric radar echoes. IEEE Transactions on geoscience and remote sensing. 30(5): 950-959. Zhang, J., Sang, L., Li, X., Wang, H., Li, Y., 2020. Design and Implementation of Raw Data Compression System for Subsurface Detection SAR Based on FPGA. Journal of Geovisualization & Spatial Analysis. 4(1):2. Zhang, Y., 2016. Study on the Formation Mechanism of Loess Collapse Induced by Rainfall in Yan’an Area. Master’s Thesis, Xi’an University of Science and Technology, Xi’an, China. Zhou, Y, Y., Chen, M., Gong, H, L., Li, X, J., Yu, J., Zhu, X, X., 2017. The Subsidence Monitoring of Beijing-Tianjin High-speed Railway Based on PS-In SAR. Journal of Geo-Information Science. 19(10), 1393–1403.

24

The authors declare no conflict of interest.

25