Analysis of surface deformation and driving forces in Lanzhou

3.1 Basic theories of SBAS-InSAR

SBAS-InSAR was proposed by Berardino et al. [21] and Lanari et al. [22]. SBAS-InSAR is a time-series analysis method that combines data to obtain short space baseline differential interferogram datasets. These differential interferograms can overcome spatial decorrelation phenomena. Using singular value decomposition (SVD) to solve the deformation rate, isolated SAR data sets separated by large spatial baselines can be connected to improve the time sampling rate of the observed data. The high-density temporal and spatial information of SBAS-InSAR can effectively eliminate the atmospheric effect phase, making the measurement results more accurate [23].

The basic principle is as follows:

M-scene SAR images of the same region are obtained in the time period from t₁ to t_M, one of which is selected as the common main image, and then n-scene interferogram is generated according to the interference combination principle, which satisfies the following relationship:

The ith-scene (i = 1, 2,…,N) interferogram generated from the main image called A and the minor image named B, and the interference phase generated at point (x, r) can be expressed as follows:

where tA and tB (tA > tB) are the acquisition time of SAR image corresponding to the ith interferogram; Δφdefi(x,r) is the deformation on the slope range corresponding to tBtA⋅Δφεi(x,r) is the terrain phase error; Δφαi(x,r) is the atmospheric phase error; and Δφnoii(x,r) is the noise phase error.

Assuming that the deformation rate between different interferometric graphs is vi,i−1, the cumulative shape variables of tB to tA can be expressed as follows:

Three-dimensional phase unwinding of the interferograms of N-scene SAR images can be used to calculate the deformation rates of different SAR image acquisition times.

This paper used 32 Sentinel-1A SLC images covering the study area from March 2015 to January 2017. The experimental platform of this article is ENVI5.3. There are six steps of SBAS-InSAR in ENVI. The first step is the connection graph generation. This step defines the combination of pairs (interferograms) that will be processed by the SBAS. Given N acquisitions, the maximum theoretical available connections are (N*(N − 1))/2. The super master will be automatically chosen among the input acquisitions. Image 2016/02/13 was automatically selected as the super master image, with a maximum time baseline of 200 days, the range looks of 4, azimuth looks of 1. The super master is the reference image of the whole process, and all the processed slant range pairs will be co-registered on this reference geometry. The second step is interferometry, which is to generate a stack of unwrapped interferograms. All of the interferograms are finally co-registered on the super master geometry and ready for the refinement and re-flattening tool and the SBAS inversion kernels. To increase the SNR of the interferograms and provide a more reliable coherence estimation, the multilooking is 4:1. The unwrapping method for the SBAS is the Delaunay MCF, this method works well for the connection of groups of high coherence pixels to other isolated high coherence groups. The third step is refinement and re-flattening. This step is executed to estimate and remove the remaining phase constant and phase ramps from the unwrapped phase stack. The fourth step is an inversion to the first step. This step implements the SBAS inversion kernel that retrieves the first estimate of the displacement rate and the residual topography. Moreover, a second unwrapping is done within this stage on the input interferograms to refine and improve the input stack because of the next step. We chose the most robust inversion model: linear model. Coherence thresholds is an important criterion for evaluating the quality of interference [23,50,51]. The coherence threshold in this step is 0.75. This step will get the estimated deformation rate. The fifth step is inversion second step. After the retrieval of the displacement time-series first estimation, a custom atmospheric filtering is performed on these preliminary results to recover the final and cleaned displacement time series. The purpose of the atmospheric filter is to smooth the displacement temporal signature respecting some physical properties of the atmosphere. This filter is implemented with a low-pass spatial filter, combined with a high-temporal pass filter. The sixth step is geocoding, geocoding converts results to the geographic coordinate system. 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 θ.

3.2 Basic theories of geo-detector

The geographic detector model (geo-detector) is a statistical method proposed by Wang [24,52], which can detect spatial variability and reveal driving forces. The core idea of this method is: if a factor has an important influence on the appearance of a phenomenon in space, then the factor should have a similar spatial distribution as the phenomenon. This method can not only detect the influence of a single factor but also judge the strength, direction, and linearity of the interaction across multiple factors. The geo-detector includes four detectors: differentiation and factor detection, interaction detection, risk area detection, and ecological detection.

Differentiation and factor detection can detect to what extent a factor X explains the spatial differentiation of attribute Y through the following expression:

where h = 1, L is the strata of variable Y or factor X (that is classification or partitioning); N_h and N are the number of units in layer h and the whole region, respectively. The variables σh2 and σ2 are variances of the Y values of the h layer and the entire area, respectively. The range of q is [0, 1]. The larger the value, the more obvious the spatial differentiation of Y is. If the stratification is generated by the independent variable X, the larger q value indicates the stronger explanatory power of the independent variable X on the attribute Y, and vice versa. The q value means that X explains 100 × q% of Y.

Interaction detection (that is, to identify the interaction between different risk factors Xs) combines evaluation factors X1 and X2. It is increased or reduced when the explanatory power of the dependent variable Y. The evaluation method first calculates the q value of Y, caused by two kinds of factors X₁ and X₂, respectively: q(X₁) and q(X₂). It then calculates their interaction q-value: q(X₁ ∩ X₂) and compares q(X₁), q(X₂), and q(X₁ ∩ X₂).

The relationship between the two factors can be divided into the following categories (Table 2):

Risk zone detection uses the t-statistic to determine whether there is a significant difference in the mean value of attributes between the two subregions.

where Y¯h indicates the properties within the subdomain h (averaged), nh is the number of samples in subregion h, and Var is the variance. The t statistic approximately obeys deformation’s t distribution, and the calculation method of the degree of freedom is as follows:

Null hypothesis H0: Y¯h=1 = Y¯h=2, If H0 is rejected at confidence level α. Two child attributes show that there are significant differences between regions.

Ecological detection is used to compare whether there is a significant difference in the influence factor X₁ and X₂ on the spatial distribution of attribute Y, which is measured by the statistic F.

where NX1 and NX2 represent the sample sizes of factors X₁ and X₂, respectively, and SSWX1 and SSWX2 represent the sum of intra-layer variances of the layers formed by X₁ and X₂, respectively. L1 and L2 represent the number of layers of variables in X₁ and X₂, respectively. Null hypothesis H0: SSWX1 = SSWX2. If H0. is rejected at the α significance level, indicating a significant difference in the effect of factors X₁ and X₂ on the spatial distribution of the attribute Y. Geo-detector was used to interpret the interpretation of single-factor and multifactor effects.

The uses of geographical detectors are as follows:

1. Data collection and arrangement: these data include dependent variable Y and independent variable data X. The independent variable is type quantity, the independent variable is discretized by the Jenks Natural Breaks (Jenks).

2. The sample (Y, X) was read into the geographic detector software, and then the software was run. The results mainly consisted of four parts: differentiation and factor detection, interaction detection, risk area detection, and ecological detection. Differentiation and factor detection, interaction detection, and ecological detection were analyzed in this paper.

3.3 Data processing

Aerial survey elevation was used to verify the accuracy of SBAS-InSAR. To simplify data processing, this paper used the same shapefile to cut elevation data from the aerial survey and SBAS-InSAR, to calculate the average of aerial survey elevation and SBAS-InSAR elevation data, and to compute the root-mean-square error (RMSE). The Department of Resources of Lanzhou offered eight aerial survey elevation sites, as shown in Figure 2.

In this paper, to simplify the work of analysis, the study area is divided into 34 grids by finishnet in ArcGIS 10.6, as seen in Figure 4a. The density of the road network and the built-up area was calculated for each net and analyzed in the next section; 30 random points in each grid were generated according to the divided study area. Using the generated random points to extract the attributes of deformation rate, temperature, precipitation, road network density, land cover types, and built-up area, but some random points can’t extract attributes for no attribute, there are 788 random points remaining after removing invalid points (Figure 4a). The extracted attributes were used to detect the spatial differentiation of surface deformation in the main urban area of Lanzhou. The Independent variable is a numerical quantity, they need to be discretized. We used the Jenks Natural Breaks (Jenks) method to reclassify road network density, built-up area, land cover types, precipitation, and temperature. Road network density and built-up area were reclassified into 20 categories due to their large differences in values. The land cover types of data were processed in seven categories. The precipitation and temperature were reclassified into six and eight categories, respectively. The X data referred to density of road network, built-up area, land cover types, precipitation, and temperature. The Y data referred to the deformation rate. X and Y data were imported into the geo-detector for calculation.

Figure 4

Results of data processing: (a) grids and random points, (b) road network, (c) density of road network, (d) built-up area, (e) built-up area in grids, (f) land cover types, (g) precipitation, (h) temperature.

Road maps were obtained from the Department of Resources of Lanzhou, mainly including urban highways, highways, state roads, pedestrian paths, nine grade roads, provincial roads, railways, county roads, and township roads, to a total of 510.978 km. The roads were merged into a new layer (Figure 4b). The density of the road network was calculated by road length. The density of the road network was calculated by the grid and the road length, the road length in each grid was calculated, and then the road length was divided by the grid area, finally, the road network density was obtained (Figure 4c).

Data for the built-up area (Figure 4d) were extracted from Landsat 8 OLI. First, experiments to monitor change in the built-up area Lanzhou from 2015 to 2017 were completed using Landsat 8 OLI images during the same period. After pre-processing the images, they were classified through the Random Forest method, and the ground objects were divided into a built-up area, roads, green spaces, water bodies, and others. Second, the classification results were input in the change monitoring process, monitoring results of the changes in the main urban area of Lanzhou from 2015 to 2017 were obtained, where the converted into the built-up area were extracted as the built-up area in this paper. To facilitate the calculation, the area of the construction area in the grid was calculated as the area of the construction area of each grid (Figure 4e).

This paper used the 2017 global 10 m resolution land cover types map (Figure 4f) released by Tsinghua University and the deformation rate map to analyze the relationship between land cover types and surface deformation in the main urban area of Lanzhou. We used the study area vector to cut the cover types map into seven types: cropland, forest, grass, shrub, water, impervious, and bare land. Cropland and shrub accounted for a relatively small scale, so cropland, forest, grass, and shrub were merged into vegetation for the convenience of research. Bodies of water in SBAS-InSAR lose coherence in the deformation rate graph, so no research was conducted on them in this paper. Therefore, the relevant land cover types in this study were as follows: vegetation, impervious, and bare land. The deformation rate of the three land types was obtained by using the three types of land to cut the deformation rate map separately.

Precipitation (Figure 4g) and temperature (Figure 4h) were applied to verify the impact of meteorological factors on the surface deformation of the main urban area of Lanzhou. Precipitation and temperature grid data were obtained from the Center for Climatic Research, Department of Geography, University of Delaware, Newark. The spatial resolution of the grid data is 0.25 degrees. The precipitation and temperature raster maps of the study area were obtained by the interpolation of grid data through inverse distance weighting (IDW). Then, the average monthly precipitation and temperature of the study area were compared with surface deformation to calculate a correlation.

4.1 Precision evaluation

Since level data could not be obtained, the paper used local surveying and mapping department elevation and field data to verify the SBAS-InSAR results. The mapping accuracy of local surveying and mapping department elevation is 1:1,00,000. There were two ways to evaluate the accuracy of the result: horizontal accuracy and elevation accuracy, which were evaluated separately. In practical applications, only the elevation accuracy needs to be evaluated. The RMSE calculation is simple and easy to understand, and it can describe the dispersion degree of terrain parameters and true values from the whole [53]. Therefore, the RMSE measures of the two groups of elevation were compared and analyzed in this paper. As Figure 5 demonstrates, the results showed that the difference between the two groups of data was very small (between −2 and 2), and the RMSE was 1.17, indicating that the results of this study have high reliability.

Figure 5

The elevation evaluation.

The research team went to the field to investigate the deformation of Lanzhou (Figure 6). According to the results, Country Garden and Jiuzhou are more severely deformed, and the Lanzhou west station, which is a high-speed railway station in Lanzhou, was not as serious. The researchers found the deformation of these locations to be consistent with the InSAR results, and the deformation of Country Garden and Jiuzhou was identifiable to the human eyes. The types of deformation mainly were cracks, subsidence, and collapses. In particular, the road cracks are very common, with a width of about 5–10 cm and a length of several meters. Wall crack width is several centimeters, land subsidence tens of centimeters (Figure 6).

Figure 6

Field evaluation.

4.2 Deformation results

Based on SBAS-InSAR technology, the time-series deformation map and deformation rate map of the study area from March 2015 to January 2017 were obtained.

Figures 5 and 6 show, respectively, the deformation rate and the time-series deformation of the main urban area of Lanzhou from March 2015 to January 2017. The maximum deformation rate was −26.50 mm year⁻¹, and the maximum rate of increase was 9.80 mm year⁻¹. The accumulative deformation was −60.14 mm. From the perspective of its spatial distribution (Figure 7), the main urban area of Lanzhou was stable, but a few regions were unstable. The deformation of the Chengguan district was mostly concentrated in the area around the Nanhuan road, Dongfanghong square, Jiuzhou, and Country Garden. The deformation of Qilihe district was mainly in the Dachaiping and Yujiaping areas. The deformation of the Anning district was mainly in the Lanzhou North Freight Yard, and the deformation of the Xigu district was mainly in the Liuquan Town.

Figure 7

Deformation rate of the study area.

Figure 8 shows the rate of deformation over time. It can be concluded that the first deformation began on 2015/03/14 in Lanzhou Country Garden, Jiuzhou, North Freight Yard, Yanjiaping, and Dachaiping. As time progressed, deformation in these areas gradually grew, and the range of deformation gradually expanded. By 2016/07/30, the uplifting trend of Xigu district had intensified, and some areas with large deformation in the central city (Dongfanghong Square) had begun to stand out. By 2016/10/16, the uplifting trend of Xigu district had slowed down. The time-series deformation peaked by 2017/01/20. The partial deformation of the central city was further aggravated.

Figure 8

Time-series deformation.

In this paper, the raster deformation rate maps with a coherence of 0.70 in SBAS-InSAR results were converted into vector points, covering a total of 415,893 vector points in the study area (see Table 3), and the vector points with a deformation rate of −5 to 5 mm year⁻¹ accounted for the vast majority of the deformation rate (96.90%). Only a small number of points were between −5.00 to −26.50 mm year⁻¹ and 5.00–10.00 mm year⁻¹, indicating the main urban area of Lanzhou was stable from March 2015 to January 2017, but there were also some regions with large deformation, which deserve further study.

5.1 The analysis of differentiation and factor detection

Table 4 describes the driving coefficient of each driving force and its explanatory force. The driving coefficient q is the highest in the built-up area and the lowest in land cover types. The p-value represents a significant test. The smaller the P, the higher the accuracy of the data. Therefore, the built-up area and the density of the road network are the main driving forces for surface deformation in Lanzhou from 2015 to 2017. The deformation rate for the built-up area is interpreted as 40.10%, while the interpretation for the density of road network interpretation is 39.65%. The deformation rate for temperature and precipitation is interpreted as 12.90% and 15.30%, respectively, but the influence of these factors on the deformation rate of Lanzhou cannot be ignored. The actual influence of temperature and precipitation higher than the experimental value, since the resolution of the meteorological data, is insufficient leading to differences in the spatial distribution of the meteorological data. This paper further analyzed the spatial distribution and cause of each driving force.

5.1.1 Density of road network and surface deformation

It is essential to analyze the relationship between surface deformation and density road network for road surface deformations that have a significant effect on the speed profile of vehicles and traffic flow conditions [3,4]. This paper studied the density road network to explain the reason for surface deformation. The density of the road network is between 1.54 and 11 km km² (Figure 9a). The density of the road network has 19 areas between 5 and 11 km. They are concentrated in the Chengguan district and Qilihe district. There were 19 time-series deformation of grids greater than 20 mm (absolute value). The density of the road network is greater than 5, and the time series form variables are greater than 20 mm at area intersections (Figure 9c), amount to a total of 12 areas (5, 9, 10, 12, 13, 15, 16, 17, 18, 20, 27, 29). The density of the road network is likely to the major cause of Nanhuan Road, Dongfanghong Square, and Dachaiping’s deformation, as shown in Figure 9a. Soil deformation and stratum movement are caused by the loading and unloading on the ground, which may affect the surface structure, the relatively concentrated ground load is an important factor of road deformation [10,54,55].

Figure 9

Density of road network and time-series deformation: (a) density of road network; (b) time series deformation; (c) the area where the density of road network >5 km/km² and time series d > 20 mm; (d) the road network in Lanzhou.

5.1.2 Built-up area and surface deformation

Urbanization is the focus of many Chinese scholars [56,57]. Population density and built-up areas are often important indicators of urban sprawl. After reviewing the statistical yearbooks in Lanzhou from 2015 to 2017 [58,59], we found that economic output increased from 20.093 to 252.354 billion yuan from 2015 to 2017, urban population density increased from 3,514 people km² to 3,576 people km²; the urbanization process was fast. To analyze the urbanization effect on surface deformation in Lanzhou, this paper analyzes the spatial relationship between the built-up area and the surface deformation. The relationship reflects the relationship between urbanization and surface deformation.

As shown in Figure 10d, the built-up area of Lanzhou is 19.38 km². Using Fishnet, the study divides these areas into 34 grids: the area of the built-up area of each grid is shown in Figure 10a, the cumulative deformation is shown in Figure 10b, with time-series deformations larger than 20 mm and the built-up area larger than 0.8 km² selected. As shown in Figure 10c, the surface deformation of the 8, 10, 18, 26, and 27 regions may be related to changes in the built-up area [60]. Figure 10c/8 and Figure 10c/10 show the Lanzhou north freight yard and Nanhuan road, respectively, which also have a large built-up area and serious deformation. The built-up area is also an indispensable cause of surface deformation in the region. The construction of a large number of projects, including the development of underground spaces 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 of the rock and soil, leading to surface deformation [10].

Figure 10

Built-up area and time-series deformation: (a) built-up area from change detection; (b) time series deformation in Lanzhou; (c) the area where the built-up area (2015–2017) >0.80 km² and time-series deformation >20 mm; (d) the area of the built-up area in Lanzhou.

5.1.3 Land cover types and surface deformation

Figure 11 shows that impervious land (Figure 11a) accounts for the largest section of the study area, followed by vegetation (Figure 11b), the urban area is relatively evenly distributed, and other areas are symmetrically distributed. In general, vegetation coverage in Lanzhou is low. Finally, the bare land is distributed mainly in Jiuzhou and Country Garden (Figure 11c). After observing the optical image, the reason for this distribution is continuous land-creation projects.

Figure 11

The land cover types with deformation rate: (a) deformation rate in impervious; (b) deformation rate in vegetation; (c) surface deformation in bar land; (d) the deformation rate in Lanzhou.

The area of impervious land (Figure 11a) accounts for the largest sectionof Lanzhou, and the deformation rate is between −5 and 5 mm year⁻¹ (Figure 11a). Surface deformation varies in the North Freight Yard, Dachaiping, Yujiaping, Jiuzhou, and Country Garden Large, between −5 and −15 mm year⁻¹, and the deformation rate of Jiuzhou and Country Garden varies between −15 and −26.50 mm year⁻¹. The vegetation area is small and the distribution is relatively uniform, with a deformation rate mainly concentrated between −5 and 5 mm year⁻¹. Surface deformation in vegetation is relatively small and relatively stable (Figure 11b). Bare land is mainly distributed in North Freight Yard, Jiuzhou, and Country Garden (Figure 11c), and the surface deformation of these areas is also serious. The land cover types are more obvious in bare land. The deep reason is the continue of land-creation. According to relevant scholars, human settlements [61] and industrial areas are inextricably linked to surface deformation. Because groundwater exploitation in human activities is also serious, leading to a decline in groundwater level. Groundwater decline and surface deformation are closely related [29,62].

5.1.4 Meteorological factors and surface deformation

Precipitation and deformation values are also related (Figure 12a). In the winter between 2015 and 2016, there was less precipitation in the study area and more deformation, whereas in the summer (July, August, and September), precipitation was great and deformation as small, especially in 2016, when precipitation was more pronounced. Heavy precipitation supplements groundwater, thus reducing surface deformation [35].

Figure 12

Time-series deformation and average monthly precipitation and temperature. (a) Average time-series deformation and the monthly mean precipitation. (b) Average time-series deformation and the monthly mean temperature.

The temperature rises from February to August and falls from September to January (Figure 12b). The highest temperature is about 20°C in the summer and the lowest temperature is about −10°C in the winter. As the temperature rises, the surface sinks, and with the decrease of temperature, the surface shows an upward trend. This phenomenon is mainly due to the seasonal temperature decrease in the Lanzhou area: freezing soil causes the volume to expand. As the temperature rises, the frozen soil gradually melts and the volume shrinks, leading to surface deformation [8,38].

Furthermore, to quantitatively study the relationship between the time-series deformation and meteorological factors, the correlation between precipitation and temperature and time-series deformation is analyzed, through a linear equation and correlation coefficient (R), as shown in Figure 13. Time-series deformation has a clear negative correlation with precipitation and temperature. The correlation coefficient (R) of the precipitation and time-series deformation is −0.61 (Figure 13a), and the correlation coefficient (R) between temperature and time-series deformation is −0.583 (Figure 13b). The correlation between precipitation and time-series deformation is stronger than temperature, indicating that precipitation has a greater impact on surface deformation.

Figure 13

The correlation between precipitation, temperature, and deformation. (a) Correlations between average time series deformation and monthly mean precipitation. (b) Correlations between average time series deformation and monthly mean temperature.

5.2 The analysis of ecological detection

Figure 14 depicts the result of ecological detection, which means the difference in the combined effects of the driving forces on surface deformation: the effects of temperature, precipitation, density of road network, land cover types, and built-up area show varying levels of influence on surface deformation are significantly different. There is no major difference in the influence of the built-up area and density of road network on surface deformation, while there is a significant difference in the influence of built-up area and temperature, precipitation, and land cover types on surface deformation, and there is also a significant difference in the influence of road and temperature, precipitation, and land cover types. In addition, there is no significant difference in the influence of temperature, precipitation, and land cover types, and no significant difference in the influence of precipitation and land cover types. From the above, it can be concluded that the impact of the built-up area and density of the road network on surface deformation is consistent. Temperature and precipitation are both meteorological factors, and their effects on surface deformation are consistent. It can be found that the influences of land cover types and temperature and precipitation on surface deformation are also consistent (Figure 14).

Figure 14

The difference in the combined effects of the driving forces. N: there is no significant difference in the influence of two driving forces; Y: there is a significant difference in the influence of two driving forces.

5.3 The analysis of interaction detection

Table 5 describes the results of interaction detection. Interaction detection identified the interaction between different risk factors Xs. From Table 5, it can be concluded that the spatial distribution and differentiation of surface deformation in the main urban area of Lanzhou is not only affected by a single driving factor but a result of the interaction of multiple driving factors, whose combined effect is more significant than any single driving factor. The interactive explanatory power between the built-up area and land cover types is 0.55, which demonstrates, the greatest impact on surface deformation, followed by the density of road network ∩ land cover types and the density of road network ∩ temperature, which are 0.54 and 0.52, respectively. Compared with the explanatory power of the single driving factor, all driving factors have an enhanced effect on the spatial distribution and differentiation characteristics of surface deformation after mutual interaction, and the effect is not independent. The interaction between the driving factors is a complex superposition relationship, rather than a simple operational relationship. It is worth noting that the interaction patterns between land cover types and density of road network, temperature, precipitation, and built-up area are a nonlinear enhancement, because their explanatory power is quite different from that of other driving factors [63].