An identification method of potential landslide zones using InSAR data and landslide susceptibility

Abstract

Landslides are destructive to property, infrastructure and people in potential landslide zones. Identifying potential landslides is an important step in landslide preparedness and will help develop sustainable landslide risk management. Interferometric synthetic aperture radar (InSAR) and landslide susceptibility assessment (LSA) have poor reliability in individually identifying potential landslide zones. This study proposes a threshold model for identifying potential landslide zones that fuses the InSAR two-dimensional (vertical direction and east-west direction) deformation rates and LSA results. The deformation rate threshold for this threshold model is |D_U| or |D_E|>10 mm/year (D_U and D_E are the vertical and the east-west deformation rates, respectively), with threshold levels of LSA set to high and very high susceptibility. The criterion of potential landslide zones is ((|D_U| or |D_E|>10 mm/year) ∩ (high or very high susceptibility of LSA)), and points with similar deformation and susceptibility are clustered by the K-means algorithm, and the potential landslide zones are obtained by elimination, smooth and speckle removal operations. The results showed that the InSAR two-dimensional deformation rates D_U and D_E were −32.71 − 12.72 mm/year and −14.88 − 24.81 mm/year, respectively, during 2015–2020 in Lanzhou city. The LSA showed that very low, low, medium, high, and very high susceptibility accounted for 55.36%, 10.54%, 21.37%, 9.63%, 3.1% of the total area, respectively. Using the proposed threshold model, 117 potential landslide zones were identified in Lanzhou. The overlap rate between potential landslide zones and the landslide inventory was 40.17%, indicating that about 40% of the potential landslide zones overlapped with the landslide inventory and that about 60% were new potential landslide zones in Lanzhou. The feasibility of the threshold model in identifying potential landslides was confirmed by field research and time-series InSAR analysis on typical areas (L1, L2, L3, and L4), which had large deformation variables and landslide features. The proposed method can quickly determine the spatial location of potential landslides, providing targeting data for landslide field investigations, technical support for rapid early landslide identification, and data support for landslide management and prevention in Lanzhou.

Keywords:

• Potential landslide zones
• SBAS-InSAR
• landslide susceptibility assessment
• threshold model

1. Introduction

A landslide is the phenomenon in which the material forming a slope (such as rock, sediment, artificial fill, or a combination of these) is moved downward and outward by the action of gravity and other influencing factors (Varnes 1978). Landslides can not only damage human life, property, and buildings but also cause secondary disasters such as weirs (Xu et al. 2008; Huang and Fan 2013; Dai et al. 2020). Between 1995 and 2014, a total of 3876 landslides killed 163658 people and injured 11689 people worldwide (Haque et al. 2019). Landslides are caused by a variety of factors, and there is a great deal of uncertainty about when and where they will occur, thus, it is important to accurately identify where landslides are likely to occur (Huang 2004; Novellino et al. 2021).

Lanzhou is located at the intersection of the Qinghai Tibet Plateau, Loess Plateau, and Inner Mongolia Plateau, with complex geological conditions, intense human engineering activities, and damage to the ecological environment, resulting in unstable slopes and increasing the likelihood of landslides. According to an official survey, the Linxia-Yongjing-Lanzhou section of the Yellow River was one of the main landslide distribution areas (The People’s Government of Gansu Province, 2017). In recent years, many scholars have conducted research on landslides in Lanzhou city (Lyu et al. 2018; Torizin et al. 2018; He et al. 2020; Shu et al. 2021). These studies were mainly based on landslide susceptibility evaluation, but they lacked related research on potential landslide identification. The landslides in Lanzhou moved as heavy rains fell, posing serious threats to the safety of local people’s lives, property, and vital infrastructure. Therefore, it is urgent to identify potential landslide zones in Lanzhou.

Small baseline subset interferometric synthetic aperture radar (SBAS-InSAR) (Berardino et al. 2002) overcomes the disadvantages of traditional deformation monitoring methods. The traditional methods are difficult and costly to deploy monitoring points and cannot achieve large-scale landslide monitoring (Berardino et al. 2002). SBAS-InSAR has been widely used to monitor ground surface deformation caused by earthquakes (Grandin et al. 2017), landslides (Dai et al. 2016, 2019), ground subsidence (Wang 2020; He et al. 2021a, 2021b), and volcanic movements (Lee and Lee 2018). SBAS-InSAR technology is suited for tracking small terrestrial deformations over large areas, which can be used to accurately identify the location and extent of early landslide deformation. Deformation traces are being used for the early identification of potential landslide zones (Tofani et al. 2013). Dai et al. (2022) used the SBAS-InSAR to obtain a wide-area potential landslide early identification the whole of Mao County, a mountainous area in western Sichuan (China), with a total of 41 potential landslides successfully detected. Zhu et al. (2021) used SBAS-InSAR to study the pre- and post-failure geomorphic changes and spatiotemporal evolution of loess landslide. Tian et al. (2022) studied the deformation and stability of loess landslides in the Yellow River Basin using SBAS-InSAR. SBAS-InSAR is used to analyse and conduct long-term monitoring of key landslides, which is an effective means of identifying deforming landslides, but it is not easy to directly identify deformed landslides.

Landslide susceptibility assessment (LSA) is a common method of landslide spatial prediction. Here, we used landslide inventory and landslide influencing factors to construct a model for predicting the spatial distribution of potential landslides . There are four main methods of LSA: physics-based models, viewpoint-driven models, statistical models, and machine learning models (Chang et al. 2019; Li et al. 2019; Pham et al. 2019; Thomas et al. 2021). Machine learning is a solution for addressing big-data spatial analytics in which the extent of the theoretical knowledge of a problem is incomplete and when statistical pre-assumptions are unreliable or unknown (Dou et al. 2020). It is an ideal technique for solving non-linear geo-environmental issues due to its robustness. Machine learning methods can quickly predict landslide susceptibility and have been widely used (He et al. 2020; Merghadi et al. 2020; Pham et al. 2020; Rodrigues et al. 2021). Multilayer perceptron (MLP) is widely used in LSA due to its parallel processing, good fault tolerance, and very strong adaptive capability. MLP can be used to the construct accurate landslide susceptibility maps (Zare et al. 2013; Yi et al. 2020; Tran et al. 2021). However, reliable spatial location identification of potential landslide zones is still a challenging task.

The combination of deformation monitoring and LSA can improve the quality of results (Zhao et al. 2019), and most studies of this kind tend to update the landslide susceptibility map by integrating InSAR data into the LSA process (Huang et al. 2020). We combined the two methods in an attempt to rapidly identify potential landslide zones. Based on the results of InSAR and MLP-based LSA, a threshold model is proposed to integrate landslide results from the two methods to achieve more accurate potential landslide zone extraction. The main objectives were to: (1)obtain the two-dimensional deformation rates of the two orbitals based on Sentinel1-A descending and ascending tracks images using SBAS-InSAR; (2)obtain the LSA based on landslide influence factors using the MLP model; and (3)construct the threshold model based on the two-dimensional deformation rates and the LSA levels, identifying the potential landslide zones of Lanzhou. Field validation and time-series InSAR vertical results were used to verify the feasibility of the new model. The method proposed in this study can provide technical support for rapid early potential landslide identification, and the research results can provide data support for landslide management and prevention in Lanzhou.

2. Methods

We initially monitored the ground surface deformation rates of Lanzhou based on Sentinel1-A descending and ascending tracks images using SBAS-InSAR. The LSA was then predicted using the MLP model. We combined the ground surface deformation rates with the LSA results to construct a threshold model to identify potential landslide zones in Lanzhou. The accuracy of the threshold model was verified based on the field investigation and time series InSAR results. The flow chart is shown in Figure 1.

Figure 1. Flow chart of this study.

2.1. SBAS-InSAR

The SBAS-InSAR technique (Berardino et al. 2002) uses the unique phase information of SAR data, which is combined with accurate SAR satellite imaging geometry and orbital ephemeris parameters to invert ground surface deformation. SBAS-InSAR is a representative of the time-series InSAR technique, which analyses SAR data acquired over a long time series in the same area. SBAS-InSAR obtains discrete targets that remain stable over a long period and that can monitor the ground surface deformation information at the millimetric level, offering significant advantages for monitoring slow surface deformation.

The basic principle of SBAS-InSAR is as follows:

In the period from $t_1$ to $t_{\mathrm{M}},$ $M$ scenes of SAR images of the same region are acquired; one of them is selected as the common master image, and then $\mathrm{N}$ scenes of interferograms are generated according to the principle of interferometric combination, which satisfies the following relationship. (1) $$\frac{M}{2}\le \mathrm{N}\le \frac{M(M-1)}{2}$$(1)

The $i (i = 1, 2, \dots , N)$ view interferogram is generated from the master image $A$ and the slave image $B.$ The interference phase produced at point $(x, r)$ can be expressed as: (2) $$\begin{array}{c}\Delta {\phi }_i\left(x,r\right)={\phi }_A\left(x,r\right)-{\phi }_B(x,r)\\ \approx \Delta {\phi }_{\mathit{def}}^i\left(x,r\right)+\Delta {\phi }_{\epsilon }^i\left(x,r\right)+\Delta {\phi }_{\alpha }^i\left(x,r\right)\\ \mathrm{ }+\Delta {\phi }_{\mathit{noi}}^i\left(x,r\right)\end{array}\mathrm{ }$$(2) where $t_A$、$t_B$ ($t_A$>$t_B$) are the SAR image acquisition times corresponding to the $i_{th}$ interferogram; $\Delta {\phi }_{\mathit{def}}^i\left(x,r\right)$ is the slope-up deformation corresponding to the moment from $t_B$ to $t_A;$ $\Delta {\phi }_{\epsilon }^i\left(x,r\right)$ is the topographic phase error; $\Delta {\phi }_{\alpha }^i\left(x,r\right)$ is the atmospheric phase error; and $\Delta {\phi }_{\mathit{noi}}^i\left(x,r\right)$ is the noise phase error.

Assuming that the rate of deformation between different interferograms is $v_{i,i-1},$ the cumulative deformation variables from $t_B$ to $t_A$ are expressed as: (3) $$\Delta {\phi }_{\mathit{def}}^i\left(x,r\right)=\frac{4\pi }{\lambda }\sum _{k=t_{B,i+1}}^{t_{A,i}}(t_k-t_{k-1})v_{k,k-1}$$(3)

The interferogram of the N SAR was obtained through phase unwrapping, deformation inversion, geocoding, and time-series deformation. The following diagram shows the flow of SBAS-InSAR (Figure 2).

Figure 2. SBAS-InSAR process framework.

The first step in this process is connection map generation, in which the super master image of the descending and ascending tracks is automatically selected as the image of May 28, 2019, with the time baseline set to 200 days. The Sentinel-1A image used in this study has at least one scene per month, and the overall pairing situation is good.

The second step is the interference process. Multi-looking with 4 × 1 in the range and the direction of the azimuth. SRTM DEM with 90 m resolution is used as auxiliary data to remove topographic phase. The interferogram is filtered by Goldstein method. The Minimum Cost Flow (MCF) is used at this step to phase unwrapping. After checking the quality of the interferograms, pairs with poor interference quality were removed.

Re-flattening is performed at the third step, which uses GCPs to estimate and remove the remaining phase constants and ramps from the unwrapped phase stack. We selected GCPs based on their deformation being close to zero, confirmed by the interpretation of the unwrapped phase; GCPs point are not in the phase jump region and are evenly distributed (Ghulam et al. 2015). Consequently, 38 GCPs were applied to refine orbital errors.

The fourth step, deformation inversion uses a linear model with good robustness to invert surface deformation. The matrix decomposition singular value decomposition (SVD) method is used to transform the interferogram residuals into the temporal SAR images to estimate the average displacement rate and residual topography. During the first inversion, the displacement rate and residual topography are first estimated. In the second inversion, based on the deformation rate obtained in the first inversion, high-pass filtering is carried out in the temporal dimension and low-pass filtering is carried out in the spatial dimension, so as to estimate and remove the atmospheric phase and obtain the final displacement result on the purer time series.

The SBAS-InSAR inversion results are geocoded in the fifth step, which removes the dummy value, and the deformation results are interpolated. SBAS-InSAR derived deformation is in the line-of-sight (LOS) direction. The LOS deformation directions cannot reflect the true deformation, so two-dimensional deformation (vertical and horizontal) can be solved by using the LOS deformation values of the ascending and descending orbits (Ferretti and Prati 2000). As Sentinel1-A is a near-polar satellite, it is not sensitive to monitor the north-south direction, so the vertical and east-west direction deformation can be obtained. The solution formula is as follows: (4) $$\{\begin{array}{c}{D}_{\mathit{LOSA}}=D_U·\cos {\theta }_A-D_E·\sin \left({\alpha }_A-3\pi /2\right)·\sin {\theta }_A\\ D_{\mathit{LOSD}}=D_U·\cos {\mathrm{\theta }}_D-D_E·\sin \left({\alpha }_D-3\pi /2\right)·\sin {\theta }_D\end{array}$$(4)

In the formula, θ is the incident angle of the satellite, and α is the azimuth angle of the flight direction of the satellite. Subscript A stands for ascending orbit, and subscript D stands for descending orbit.

2.2. LSA based on the MLP model

2.2.1. Construction of the landslide inventory and landslide influencing factors

In this study, we identified 152 loess landslides using Google Earth images combined with fieldwork, with the size of the landslides ranging from 0.001 km² to 0.52 km², reaching a total area of 7.18 km², and we constructed the landslide inventory (Figure 3). Based on the principles of systematicity, representativeness, hierarchy, and operability (Chen 2016), eight landslide influencing factors were selected (Table 1). Variance inflation factors and tolerances (VIF&T) were used to test the covariance of the influencing factors, ensuring that the factors were independent of each other. (Hair et al. 1998).

Figure 3. Landslide inventory in Lanzhou city.

Table 1. Data sources of eight landslide influencing factors.

Data name	Data source institutions	URL	Influencing factors	Data type	Resolution /Scale
SRTM DEM	United States Geological Survey (USGS)	earthexplorer.usgs.gov	Elevation	Raster	30m
			Slope
			Aspect
Land cover types	Tsinghua University	http://data.ess.tsinghua. edu.cn/	Land use type	Raster	30m
Lithology and fault	China Geological Survey (CGS)	https://www.cgs.gov.cn/	Lithology	Shapefile	1:4,000,000
			Distance from fault	Raster	30m
River and road	Department of Resources of Lanzhou	http://zrzyj.lanzhou.gov.cn/	Distance from river		30m
			Distance from road		30m

The VIF&T formula is as follows： (5) $$\mathit{VIF}=\frac{1}{1-A^2}=\frac{1}{T}$$(5) where $A^2$ is the variance between the independent variables, the variance inflation factor (VIF) is the ratio of the variance of the regression coefficients when the regression coefficient equation is not correlated with the independent variables, and tolerance (T) is the inverse of VIF. The larger the standard deviation, the higher the degree of multicollinearity. When VIF > 10 or T < 0.1, it indicates that there is a problem of collinearity in the selected data.

In this study, the covariance of the influencing factors was tested for covariance. VIF&T are obtained through computational analysis, as shown in Table 2.

Table 2. Collinearity evaluation table of evaluation factors.

Number	Landslide factors	Variance inflation factor (VIF)	Tolerances (T)
1	Slope	1.302	0.768
2	Elevation	1.978	0.506
3	Aspect	1.007	0.993
4	Lithology	1.391	0.719
5	Distance from fault	1.088	0.919
6	Distance from river	2.135	0.468
7	Distance from road	1.608	0.622
8	Land cover types	1.082	0.924

Table 2 shows that the VIF of the eight influencing factors selected in this study were between 1 and 3, and T is between 0.4 and 0.9. This shows that the VIF and T of the eight landslide influencing factors met the requirements of the critical values, there was no covariance problem and they had good independence (Zhou 2019).

2.2.2. MLP model

MLP is a feed-forward supervised learning network consisting mainly of the input layer, hidden layer, and output layer (Figure 4). The hidden layer consists of artificial neurons that form a full connection between each layer, and this neural network can approximate any continuous function in the real number range. The MLP network is a function of one or more predictors that minimis prediction error.

Figure 4. MLP network diagram.

MLP is expressed as (6) $$f\left(x\right)=G(b^{\left(2\right)}+W^{\left(2\right)}(\sigma (b^{\left(1\right)}+W^{\left(1\right)}x)\left)\right)$$(6) where $x$ denotes the input layer vector, $W$ is the weight, $b$ is the bias coefficient, $\sigma$ is the $\mathit{sigmoid}$ function, and $G$ is the $\mathit{softmax}$ function, The formula is: (7) $$\mathit{sigmoid}\left(x\right)=\frac{1}{1+e^{-a}}$$(7) (8) $$\mathrm{softmax}\left(x_j\right)=\frac{e^{x_j}}{\sum _{k=1}^K e^{x_k}}\mathrm{ }\text{for}\mathrm{ }\mathrm{j}=1,\dots ,\mathrm{K}$$(8)

2.3. Threshold model

The threshold model sets different thresholds for two-dimensional deformation rates and LSA levels. The zones that satisfy both the two-dimensional deformation rate threshold and the LSA threshold are potential landslide zones. The potential landslide zones are clustered and smoothed to eliminate errors. The detailed process of the threshold mode construction is shown in Figure 5.

Figure 5. The process of the threshold model construction.

The study area is divided into 30 × 30 m grids by the fishnet function in ArcGIS10.5. Each grid contains centrally located points, and the ground surface deformation rate data and LSA results are assigned to these points by extracting multi values to points function in ArcGIS10.5.

Surface deformation rates are one of the common methods for identifying landslides, and the deformation characteristics of landslides vary from region to region. According to the findings of other researchers in the same region, the surface deformation rate of −10 to 10 mm/year is stable, with a low probability of landslides (Li et al. 2021). We adopt this threshold in our study. The two-dimensional deformation rate threshold for the threshold model in this study is |D_U| or |D_E|>10mm/year (D_U and D_E are the vertical and the east-west deformation rates, respectively) (Figure 5). The LSA is set between 0 and 1, with 0 indicating a very low possibility of landslide, and 1 indicating a very high possibility of landslide. The natural breaks method maximises the similarity within each group and the difference between external groups, this method can separate different landslide possibilities (He et al. 2020; De Smith et al. 2015). This study used the natural breaks method to classify LSA into very low susceptibility, low susceptibility, medium susceptibility, high susceptibility, and very high susceptibility. The higher the LSA, the greater the probability of landslides (Pourghasemi et al. 2018). The landslide susceptibility level of the threshold model was set to very high and high susceptibility. Potential landslide zones in Lanzhou were identified by a threshold model, and points with similar deformation and susceptibility were clustered by a clustering method. There were no isolated points; each point had at least other points that were naturally adjacent to it, and the clustered points were converted to vectors by ArcGIS 10.5. The clustering method is the K-means algorithm (Jain 2010), as follows：

Cluster the set of points $=\left\{x_i\right\},i=1,\dots ,n$ into K classes, $C=\left\{c_k,k=1,\dots ,K\right\},$ and the K-means algorithm finds a partition such that the empirical mean of the clusters and the points in the clusters are minimised. ${\mu }_k$ is the mean of the clusters $C_k,$ and the root mean square error of the points $c_{k }$ is: (9) $$J\left(c_k\right)=\sum _{x_i\in c_k} {\Vert x_i-{\mu }_k\Vert }^2$$(9)

The objective of the K-means algorithm is to minimise the sum of squared errors in all K classes. (10) $$J\left(C\right)=\sum _{k=1}^K \sum _{x_i\in c_k} {\Vert x_i-{\mu }_k\Vert }^2$$(10)

To set a priori value of K-means, we combined field investigation, SBAS-InSAR and LSA results, revealing 200 was a priori value consistent with the actual situation. The points after clustering were split into several polygons, which have the disadvantage of discontinuous polygons below 1,000 m² and tiny polygons over. In this study, we used the elimination tool to merge the polygons below 1,000 m² with the neighbouring polygons with the longest shared boundary and to make the polygons continuous. The eliminated polygons were mainly rectangles and triangles, and the smoothing distance was set to 1,000 m for smoothing (Bodansky et al. 2001). Polygons smaller than 1,000 m² that could not be eliminated or smoothed were removed.

3. Study area

Lanzhou is located at the intersection of the Qinghai-Tibet Plateau, Loess Plateau, and Inner Mongolia Plateau. Lanzhou city is long from east to west and narrow from north to south (Figure 6). The city is flanked by the north and south hills, with the northern bank steeper and a single terrace, and the southern bank gently sloping and a mixed terrace. The geological structure of this area is complex, and the loess layer and gravel layer are widely distributed and mixed with a small part of the fine sand conglomerate, loam, and volcanic rock. The loess layer is mainly Malan loess with a loose structure and easy collapsibility. The city is located in an area with a temperate monsoon climate, with the precipitation concentrated in summer. The characteristics of concentrated precipitation and loess collapsibility easily predispose Lanzhou to landslides (Wang et al. 2019; Wang 2020).

Figure 6. The location of the study area.

Due to the constraints of the topography, urban construction had to expand to the slopes and loess plateau areas on both sides of the river valley. A series of strong human activities, such as slope cutting for engineering construction, land creation, irrigation in the surrounding mountains, irrigation of farmland, and mining of mineral resources, increased the frequency and possibility of landslides in Lanzhou (Li et al. 2014). A total of 227 landslides were recorded in the inventory during 2017–2018, with an area of approximately 1.4 × 10⁶ m², accounting for 0.35% of the whole study area (Shu et al. 2021). More than 10 landslide disasters have caused a large number of casualties and property losses (Su et al. 2019). Loess landslides are the main landslide types in Lanzhou. Loess landslides include bedrock contact landslides, palaeosol contact landslides, mixed landslides, and slides within four loess types. Loess is a typical unsaturated soil (Fredlund and Rahardjo 1993). It is very common to see that very steep slopes of loess can be stable for a long time, but they are vulnerable to failure when subjected to earthquake shock, rainfall, human factors, or irrigation (Tu et al. 2009). Both laboratory tests and field investigations have shown that loess deposits can accumulate substantial pre-failure strains as a consequence of gradually accumulating micro shears in the loess fabric (Zhang and Wang 2017; Shi et al. 2019; Zhang et al. 2020). Therefore, Lanzhou faces the huge threat of landslide disasters. It is necessary to identify potential landslide zones and contribute to disaster prevention and reduction in Lanzhou. The main urban area of Lanzhou, including Chengguan district, Anning district, Xigu district, and Qilihe district, was selected as the study area (Figure 6).

Landslides in Lanzhou are distributed in urban areas and in the wild, especially in the south of Lanzhou. SBAS-InSAR can obtain better coherence and more comprehensive deformation data than permanent scatterer interferometric synthetic aperture radar (PS-InSAR). Therefore, SBAS-InSAR technology was selected in this study to monitor the surface deformation of Lanzhou.

4. Results and analysis

4.1. Accuracy verification of SBAS-InSAR results

To verify the reliability of SBAS-InSAR results, the typical areas with large surface deformation is selected for field investigation (Figure 7). It is found that the areas with large surface deformation are mainly characterized by cracks, subsidence and collapse. Road cracks are common, ranging in width from 5 to 10 cm and in length from several meters. The crack in the wall was several centimeters wide. Through UAV aerial film, it is found that there is a gully obviously collapsed, the facade reaches several meters, the collapse boundary has a toppled power pole, and the surrounding landslide is accompanied by different sizes. It indicates that the SBAS-InSAR results of this study are credible.

Figure 7. Field investigation verification.

4.2. Ground surface deformation rates based on SBAS-InSAR

InSAR can monitor minor ground surface deformation, reflecting the movement characteristics of landslides. SBAS-InSAR technology was used to monitor ground surface deformation rates in the main urban area of Lanzhou based on Sentinel1-A descending and ascending tracks data from 2015 to 2020. SBAS-InSAR derived deformation is in the line-of-sight (LOS) direction. Two-dimensional vertical and east-west ground surface deformation rates can be solved by using the LOS deformation values of the ascending and descending orbits (Ferretti and Prati 2000). According to landslide research about Lanzhou (Li et al. 2014; Xue et al. 2016; Wu et al. 2021), the ground surface deformation rate of −10 - 10 mm/year is too small to cause landslides. Thus, regions of −10 – 10 mm/year are defined as stable, and regions of more than −10 – 10 mm/year are defined as unstable. Therefore, the threshold of two-dimensional vertical and east-west ground surface deformation rates in this study is set as ±10 mm. The vertical and east-west surface deformation rates based on descending and ascending tracks are shown in Figure 8.

Figure 8. Two-dimensional surface deformation rates in Lanzhou city, (a) vertical direction, (b) east-west direction.

The two-dimensional ground surface deformation rates of vertical and east-west directions were −32.71-12.72 mm/year and −14.88-24.81 mm/year, respectively, during 2015–2020. As shown in Figure 8, the SBAS-InSAR vertical results indicate that the unstable regions spread in Qingbaishi and Kyushu in Chengguan district, Shajingyi in Anning district, and Liuquan town in Xigu district, while the SBAS-InSAR east-west results indicate that the unstable regions are distributed in Qingbaishi in Chengguan district, Yingmenqian village in Qilihe district, and Liuquan town in Xigu district. The two-dimensional SBAS-InSAR results combined descending and ascending tracks compensated for the lack of single-track detection (Liu et al. 2018) and detected unstable regions with vertical and east-west directions deformation rates greater than 10 mm/year or deformation rates less than −10 mm/year. Unstable regions are important data sources for the threshold model of landslides.

4.3. LSA based on the MLP model

LSA is a prediction of the possible landslides’ spatial location based on landslide influencing factors. In this study, LSA was generated based on eight landslide influencing factors and landslide inventory using the MLP model in the main urban area of Lanzhou (Figure 9a). The accuracy of the prediction results was verified using the receiver operating characteristic (ROC) curve. ROC is the synthesis of a confusion matrix under different thresholds and is widely used in LSA (Kavzoglu et al. 2015). In the landslide inventory true positive (TP) was predicted to be a landslide, and false positive (FP) was predicted to be a non-landslide. ROC with the false positive rate (FPR) as abscissa and the true positive rate (TPR) as ordinate, reflects the continuous changes of data specificity and sensitivity. The area under the curve (AUC) of the ROC can directly reflect the results (Figure 9b). The AUC is greater than 0 and less than 1. An AUC closer to 1 indicates that the model is more likely to correctly classify landslide and non-landslides, and closer to 0 means that the model is less likely to correctly classify landslide and non-landslide (Chen et al. 2018). The AUC for this study was 0.839, with a standard deviation of 0.007, indicating good results for LSA using the MLP model.

Figure 9. LSA distribution and AUC curve.

LSA ranges between 0 and 1, where 0 indicates a very low possibility of landslide, and 1 indicates a very high possibility of landslide. This study used the natural breaks method (He et al. 2020; De Smith et al. 2015) to classify the LSA into very low susceptibility (0 − 0.34), low susceptibility (0.34 − 0.48), medium susceptibility (0.48 − 0.63), high susceptibility (0.63 − 0.82) and very high susceptibility (0.82 − 1.00), which accounted for 55.36%, 10.54%, 21.37%, 9.63%, 3.1% of the total area, respectively (Figure 9). High and very high susceptibility areas are located in the north and south hills of the main urban area of Lanzhou. From the perspective of landslide influencing factors, in the high and very high susceptibility areas of landslides, the stratigraphic situation is that of the upper Pleistocene wind deposits, the Xiliugou formation group, and the estuary group, consisting mainly of loess, mudstone, and siltstone. The land use types are grassland and woodland, with slopes ranging from 20° to 40° and close to rivers, roads, and faults. These regions have loose geological conditions, high slopes, and intensive human activities, which are more prone to landslides and other geological hazards. Therefore, the LSA experimental results based on the MLP model are reasonable.

According to the analysis results described above, monitoring the unstable areas based on SBAS-InSAR and the high and very high susceptibility areas of LSA based on the MLP model are consistent in spatial distribution. The unstable areas of SBAS-InSAR are delicate, and the high and very high susceptibility areas of LSA are rough. However, the two results did not identify potential landslide zones. To determine the possible landslide area in Lanzhou, this study constructs a threshold model to further extract the information from the unstable areas of SBAS-InSAR and high, very high susceptible areas of LSA to obtain potential landslide zones in Lanzhou.

4.4. Potential landslide zones based on the threshold model

We used the constructed threshold model to identify potential landslide zones in Lanzhou. The identified potential landslide zones are mainly blocky and triangular, therefore, the potential landslide zones are disposed of eliminating empty, smooth boundaries and removing speckles. A total of 117 potential landslide zones were identified in Lanzhou (Figure 10). From Figure 10, potential landslide zones were mainly distributed in Xujiashan National Forest Park in Chengguan district, Jinjiling and Bapanxia stations in Xigu district and Weiling town, Shifogou national forest park, Xiguyuan town, and Jingu town in Qilihe district in Lanzhou. The overlap rate of the potential landslide zones and the landslide inventory were calculated based on pixels, and we found that the overlap rate between the potential landslide zones and the landslide inventory was 40.17%, indicating that about 40% of the potential landslide zones overlapped with the landslide inventory, and about 60% were new potential landslides. These potential landslides were unstable and had high landslide sensitivity, suggesting that they had great landslide risks.

Figure 10. Distribution of the potential landslide zone in Lanzhou city (Landslide verification 1: L1and L2 belong landslide inventory and potential landslide zone; Landslide verification 2: L3 and L4 belong potential landslide zone).

The threshold model proposed in this study performs well in the Lanzhou area. The characteristics and inducing factors of landslides in different areas are regional, and the parameters of the deformation rate and LSA design of the threshold model will also be affected by regional factors. Thus, parameter learning and adjustment of the model will be needed in different regions.

To illuminate the relationship among the potential landslide zones, landslide inventory, and landslide influencing factors, we analysed the distribution among landslide influencing factors, as shown in Figure 11. The potential landslide zones and the landslide inventory were concentrated at an altitude of 1451–2100 m (Figure 11a), in the slope of 20 - 40° (Figure 11b), on the northern, eastern, and western aspects of the mountain range, especially on the western aspect, accounting for about 1/3 potential landslide zones of the total potential landslide zones, with obvious slope divergence (Figure 11c); Moreover, the zones were mainly distributed at 2000 m from the road (Figure 11d), 4000 m from the railway (Figure 11e), 1000 m from the river (Figure 11f) and 1000 m from the a fault (Figure 11g). The upper Pleistocene wind deposits and estuary group were also distributed (Figure 11h). To summarise the above features, the potential landslide zones and the landslide inventory distribution were characterised by areas of low elevation, sunny slopes, moderate slopes, soft lithology, and proximity to man-made facilities and rivers. It is worth noting that the potential landslide zones and the landslide inventory were quite different in the distribution of landslide influencing factors; for example, the landslide inventory had a relatively higher percentage between 1700 and 1900 m compared to potential landslide zones. The same differences can be observed for the east aspect and the distance from the railway and lithology. This may be explained by the fact that potential landslide zones are differentially assessed by the LSA and SBAS-InSAR results. InSAR monitors the movement characteristics of landslides, and LSA is a prediction of the possible landslides’ spatial location based on landslide influencing factors. The landslide inventory reports the historical landslides, which may be stationary and low-risk. Therefore, some differences are expected between the potential landslide zones and landslide inventory in the distribution of the factors influencing landslides.

Figure 11. Distribution of the potential landslide zone and the landslide inventory in landslide influencing factors (a is altitude, b is the slope, c is an aspect, d is the distance from the road, e is the distance from railway, f is the distance from the river, g is the distance from fault, h is lithology).

5. Discussion

The ground surface deformation rate obtained by InSAR technology can only reflect the surface displacement of landslides, and there is an error individually in identifying potential landslides from deformation data (Carlà et al. 2016). The results of the LSA can reflect the possibility of landslides from terrain, geological, hydrological, and human engineering factors, but it cannot accurately locate the spatial location of potential landslides (Huang et al. 2020). Accurately identifying landslides has become an important issue in landslide forecasting (Novellino et al. 2021). The research idea of this study integrated the advantages of the InSAR deformation rate reflecting landslide displacement and LSA reflecting landslide stability, combining their results to build a threshold model for identifying potential landslide zones.

The threshold model reported in this study identified potential landslide zones with the advantages of fast speed, good portability, and high identification accuracy. As long as InSAR deformation data and LSA results are available, the model can quickly determine the spatial location of potential landslides in the target area and provide target data for landslides field investigation. This study used a threshold model to identify 117 potential landslide zones in Lanzhou. Compared with existing studies about landslides in Lanzhou (Lyu et al. 2018; Torizin et al. 2018; Shu et al. 2021), this study combined InSAR deformation data and landslide influencing factors to realise individualised investigation of potential landslides, and convenient for landslide prevention and control. Shu et al. (2021) used the weighted frequency ratio model to calculate landslide susceptibility for each type of landslides. Lyu et al. (2018) evaluated the risks around Lanzhou by employing the analytic hierarchy process (AHP). Torizin et al. (2018) used the weights of evidence method to calculate statistical LSA in a dynamic environment. Wang et al. (2020) used the deep belief network to calculate landslide sensitivity of Sichuan Province, China, and realised a better prediction precision. Dou et al. (2020) explored the predictive performance power of different sampling techniques in landslide susceptibility mapping in the wake of increased usage of artificial intelligence. These methods are a very important inspiration to improve the accuracy of the landslide prediction model. However, these existing studies assessed only the level of landslide susceptibility and did not identify specific potential landslides. However, the threshold model constructed in this study, based on the InSAR results and LSA results, can rapidly identify potential landslide zones. The accuracy of the threshold model in identifying potential landslide zones was verified through field investigation and InSAR analysis.

To further validate the threshold model, the team selected four identified potential landslides for field validation in December 2020 and for InSAR time series analysis. L1 and L2 were typical areas in both potential landslide zones and the landslide inventory to ensure the accuracy of the landslides. L3 and L4 were typical areas in potential landslide zones but not in the landslide inventory, ensuring the accuracy of potential landslides in potential landslide zones.

5.1. Landslide verification 1 (typical areas L1 and L2)

L1 is located at the Liuquan tunnel in the Xigu district of Lanzhou, with a vertical deformation rate of −26.41 − 7.64 mm/year during 2015–2020, and an accumulated deformation of −39.68 − 8.64 mm. The accumulated deformation fluctuated widely, indicating disorderly movement of the slope surface displacements (Figure 12a and b). This area has large mountains and deep gullies, with mountain heights of more than 300 m and slopes of more than 20°. Through field research, we found that L1 had obvious landslide traces and small shallow landslides exist on the river edge under the erosion of the river (Tong and Schmidt 2016) as shown in Figure 12c, with dense landslides around the tunnel and obvious landslide backwalls (Figure 12d). This proves the correctness of potential landslide identification, and also verifies the rationality of InSAR monitoring results.

Figure 12. L1 field investigation and InSAR analysis (a is InSAR vertical deformation rate, b is InSAR time series deformation, c and d are landslide elements).

L2 is located on the hill behind the China petrol station in Taoshuping, Chengguan district of Lanzhou. On the front edge of the slope are large advertising production factories and a petrol station. L2 vertical deformation rate was −30.56-0.94 mm/year, and the accumulated deformation rate from 2015–2020 reached −150 mm; the deformation rate continued to increase (Figure 13a and b). The field study showed that landslides had already occurred on L2 (Figure 13c–e), with significant damage to the surrounding buildings and cracks in the walls (Figure 13c and f), indicating a history of landslides and the potential for landslides to occur.

Figure 13. L2 field investigation and InSAR analysis (a is InSAR vertical deformation rate, b is InSAR time series deformation, c-f are landslide elements).

5.2. Landslide verification 2 (typical areas L3 and L4)

L3 is located at the back of the Xinhang driving school in Chengguan district. The upper part of the landslide is a residential area, and the front edge of the slope is a floodway. The deformation rate of L3 was −31.12-0.13 mm/year during 2015–2020, and the accumulated deformation rate reached −150 mm (Figure 14a and b). The unstable slope was an artificial accumulation, mainly composed of slag, rubble, etc., with a length of 80 m and a width of 50 m. After field investigation (Figure 14c–e), we found the deformation characteristics of the landslide. There were several sinkholes in the upper and lower parts of the landslide and a considerable number of longitudinal cracks in the upper part of the landslide, with a width of 30 cm, which had the possibility of landslide (Highland and Bobrowsky 2008).

Figure 14. L3 field investigation and InSAR analysis (a is InSAR vertical deformation rate, b is InSAR time series deformation, c-e are landslide elements).

L4 is located in the first coal mine in Liushuwan, Qilihe district. The deformation data revealed that the deformation rate in this area was large, with a deformation rate of −31.68 − 8.98 mm/year and a cumulative deformation rate of −253.52 mm (Figure 15a–c). The upper part of the landslide body power poles tilted violently, with longitudinal cracks around the landslide. Many small landslides were found in the village, mostly distributed along roads or gullies (Figure 15d–f). Affected by slope movement, the wall cracks were serious and tens of centimetres wide, indicating obvious landslide features in L4 (Highland and Bobrowsky 2008; Guo et al. 2019).

Figure 15. L4 field investigation and InSAR analysis (a is InSAR vertical deformation rate, c is InSAR time series deformation, b, d-f are landslide elements).

These field investigations and InSAR analyses of L1, L2, L3, and L4 typical areas revealed that the potential landslide zones identified by the proposed threshold model had landslide characteristics. The suggests that the potential landslides identified by the threshold model in this study were correctly located and had a high degree of confidence. The proposed threshold model can quickly and accurately identify potential landslide zones and provide target areas for landslide investigation, thus improving the efficiency of landslide identification, reducing errors in manual landslide identification, and providing technical ideas and data support for disaster mitigation and prevention.

However, there are certain limitations to this study. InSAR technology is out of coherence in certain areas with dense vegetation or excessive deformation (Ciampalini et al. 2016). This may lead to no deformation values in the area, thus missing these areas. In future research, L-band SAR satellites or ground deformation monitoring equipment, such as GNSS, can be considered. Two situations may pose difficulty in assessing potential landslides. The first is the region with large deformation but low landslide sensitivity, such as the surface deformation caused by artificial activities (Wu et al. 2019). The second is the region with high landslide sensitivity but small deformation, which has the possibility of a landslide; however, it is difficult to judge whether the region will have a landslide based on the movement of the landslide. Future research should focus on fine-tuning the proposed model to better identify potential landslides in these two cases.

6. Conclusion

Potential landslide identification is an important step in landslide preparedness, and more accurate potential landslide identification will help to develop sustainable landslide risk management and reduce landslide risks. This study combined the SBAS-InSAR deformation rates of the Sentinel1-A descending and ascending tracks and the LSA based on the MLP model to construct a threshold model that could accurately identify potential landslide zones. The main conclusions are as follows:

1. The two-dimensional ground surface deformation rates of vertical and east-west directions were -32.71 – 12.72 mm/year and -14.88 – 24.81 mm/year, respectively, during 2015–2020. SBAS-InSAR results indicated that the unstable regions of Lanzhou were distributed in Qingbaishi and Kyushu in Chengguan district, Shajingyi in Anning district, Yingmenqian village in Qilihe district, and Liuquan town in Xigu district.

2. Very low, low, medium, high, and very high susceptibility regions of the LSA accounted for 55.36%, 10.54%, 21.37%, 9.63%, 3.1% of Lanzhou, respectively. High and very high susceptibility regions were located in the north and south hills of the main urban area of Lanzhou. These regions have loose geological conditions, high slopes, and intensive human activities, which are more prone to landslides and other geological hazards.

3. The threshold model was constructed based on the InSAR two-dimensional deformation rates and LSA results. Using the proposed threshold model, 117 potential landslide zones were identified in Lanzhou. The overlap rate between the potential landslide zones and the landslide inventory was 40.17%, indicating that about 40% of the potential landslide zones overlapped with the landslide inventory and that about 60% were new potential landslide zones.

4. The proposed threshold model was verified through field research and InSAR analysis. The feasibility of the threshold model in identifying potential landslides was confirmed by field research on typical areas L1, L2, L3, and L4, which were areas with large deformation variables and landslide features.

5. The distribution of the potential landslide zones and the landslide inventory among landslide influencing factors was statistically analysed in Lanzhou. The potential landslide zones and the landslide inventory were distributed at low altitudes, on sunny slopes, on moderate slopes, with soft lithology, and in close proximity to man-made facilities and rivers.

The proposed threshold model presented in this paper draws on the advantages of the InSAR deformation rate reflecting landslide displacement and landslide sensitivity evaluation reflecting landslide stability to rapidly determine the spatial location of potential landslides in the study area and to provide targeting data for landslide field investigation. The proposed threshold model is highly transposable and can be flexibly adapted to the individual circumstances of the target area to achieve accurate identification of potential landslides.

Acknowledgments

The authors would like to thank the United States Geological Survey (USGS) for providing access to the SRTM DEM, Tsinghua University for providing access to the land cover types, Geological Survey (CGS) for providing access to Lithology and fault, Department of Resources of Lanzhou for providing access to River and road. At the same time, we would like to thank the reviewers and editors for their valuable suggestions and their hard work.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

All data included in this study are available upon request by contact with the corresponding author.

Additional information

Funding

This work is funded by the Project Supported by the Open Fund of Key Laboratory of Urban Land Resources Monitoring and Simulation Ministry of Natural Resources (KF-2021-06-014), Tianyou Youth Talent Lift Program of Lanzhou Jiaotong University, Education Department of Gansu Province: Young Doctoral Fund Project (2022QB-058), State Key Laboratory of Geo-Information Engineering and Key Laboratory of Surveying and Mapping Science and Geospatial Information Technology of MNR, CASM (2022-03-03), Natural Science Foundation of Gansu Province (20JR10RA249) and Jiayuguan City 2021 Science and Technology Plan Projects (437069).