An integrated neural network method for landslide susceptibility assessment based on time-series InSAR deformation dynamic features

ABSTRACT

We develop an integrated neural network landslide susceptibility assessment (LSA) method that integrates temporal dynamic features of interferometry synthetic aperture radar (InSAR) deformation data and the spatial features of landslide influencing factors. We construct a time-distributed convolutional neural network (TD-CNN) and bidirectional gated recurrent unit (Bi-GRU) to better understand the temporal dynamic features of InSAR cumulative deformation, and construct a multi-scale convolutional neural network (MSCNN) to determine the spatial features of landslide influencing factors, and construct a parallel unified deep learning network model to fuse these temporal and spatial features for LSA. Compared with the traditional MSCNN method, the accuracy of the proposed model is improved by 1.20%. The performance of the proposed model is preferable to MSCNN. The area under the receiver operating characteristic curve (AUC) of the testing set reaches 0.91. Our LSA results show that the proposed model clearly depicts areas with very high susceptibility landslides. Further, only 10.18% of the study area accurately covers 84.79% of historical landslide areas. Subjective consequences and objective indicators show that the proposed model that is integrated time-series InSAR deformation dynamic features can make full use of landslide characteristics and effectively improve the reliability of LSA.

KEYWORDS:

• Landslide susceptibility assessment
• InSAR; time distributed convolutional neural network; bidirectional gated recurrent unit; remote sensing

1. Introduction

Landslides damage or destroy houses, roads, bridges and other infrastructure, often resulting in casualties and property losses (Dai et al. 2020; He et al. 2023a). In addition, landslides destroy land ecosystem, causing environmental problems, such as soil erosion, water and soil loss and water pollution (Keefer and Larsen 2007). Early landslide positioning and prevention are the best measures available for reducing and avoiding disaster losses (Chawla et al. 2002; Mantovani et al. 2023). Landslide susceptibility assessment (LSA) explores the complex relationship between landslides and disaster-prone environments (He et al. 2021a). LSA determines the occurrence probability of unknown landslides in similar environments through comprehensive evaluation (Cantarino et al. 2019; Fan et al. 2021; Thomas et al. 2021). LSA can provide strong support to landslide disaster managers. However, the genetic mechanism of a landslide is complex, and some of the phenomenon’s characteristics have changed since the LSA method was developed (Novellino et al. 2021). After accurately prediction the position of a landslide, it is still a challenge to make comprehensive use of the landslide’s characteristics for reliable LSA modeling.

Existing LSA methods include probabilistic analysis based on landslide cataloguing (Wang et al. 2021), inference based on empirical knowledge (Niu 2021), a mathematical method based on physical mechanics (Chen, Miao, and Wu 2022) and a statistical method based on data (Huang et al. 2017). At present, the LSA method based on statistical analysis is the most commonly used because of its objectivity and data availability. Based on traditional linear statistical analysis, machine learning methods stand out by virtue of their ability to examine large amounts of data independently. Machine learning methods, such as random forest (RF) (Dou et al. 2019) and logistic regression (Zhang et al. 2019), have been extensively used in LSA. However, under the requirements of complex scenes or high precision, traditional machine learning algorithms cannot meet actual demand (He et al. 2021a; Zhao et al. 2022). Building on the neural networks present in machine learning, the neural network method effectively predicts complex nonlinear dynamic systems. It has been widely and successfully introduced into the field of LSA, including the convolutional neural network (CNN) (Wang, Fang, and Hong 2019; Gao et al. 2023a), recurrent neural network (RNN) and deep belief network (DBN) (Chen et al. 2020). The convolutional layer of CNN can extract multidimensional features from the input images and has good performance (Hakim et al. 2022). Gated recurrent unit (GRU) network of RNN variant has good performance in processing sequence data (Zhao et al. 2022). With the complexity of the environment, when faced with a limited sample, the ensemble learning model is also widely used in the LSA. For example, Wang et al. (2022) conducted the LSA based on the XGBoost ensemble learning model. Lv et al. (2022) combined CNN, DBN and ResNet models with the ensemble learning techniques of Stacking, Bagging and Boosting to generate the LSA.

The current LSA modeling process based on neural networking takes constant and landslide influencing factors, such as elevation, slope and rainfall as the model learning object (Budimir, Atkinson, and Lewis 2015). Interferometer synthetic aperture radar (InSAR) technology provides a new perspective for accurate and near-real-time landslide disaster research (Dai et al. 2016). InSAR was widely used to detect the slow movement of landslides and assess landslide susceptibility (Zhang et al. 2018; Kim et al. 2022). Meghanadh et al. (2022) input an InSAR one-dimensional line-of-sight (LOS) velocity map into the network model, thus improving the LSA. Yuan and Chen (2022) showed that the LSA was more precise when it considered InSAR deformation features. Gao et al. (2023b) systematically analyzed and proved that using InSAR deformation results improves the accuracy of LSA. Similarly, Yao et al. (2023) improved the LSA by adding InSAR deformation data into the LSA model. These studies show that InSAR data has demonstrated good value in the LSA. However, the existing studies about using InSAR deformation data in the LSA are mostly limited to static deformation velocity; the time-series InSAR dynamic features of landslide development have not been well applied or fully considered. From the perspective of feature learning, the existing methods offer inaccurate results due to the incomplete characterization of authentic landslide development. The issue of landslide variation limitations over time can be addressed by adding the time-series InSAR deformation feature to LSA models. Therefore, in this study, time-series InSAR deformation results are used as temporal dynamic features to express the dynamic development of landslides.

We propose an integrated neural network LSA method that integrates the temporal dynamic features of time series InSAR deformation information and the spatial features of landslide influencing factors for more accurate LSA. The main aims are: (1) to obtain the temporal dynamic features of time-series InSAR cumulative deformation. For this, a time-distributed CNN (TD-CNN) and bidirectional gated recurrent unit (Bi-GRU) are applied; (2) to obtain the spatial features of 10 landslide influencing factors. For this, a multi-scale CNN (MSCNN) is applied; and (3) to fuse temporal and spatial features. For this, a parallel unified deep learning network model that integrates TD-CNN, Bi-GRU and the MSCNN network model is constructed for LSA. We compare this approach with the traditional MSCNN model, using the subjective consequences and objective indicators to aid our evaluation.

2. Materials

2.1. Experimental scene

The main urban area of Lanzhou is the study area for the present research. The study area is a typical mountainous landform with high relief and a complex geographical environment (Figure 1). Lanzhou is located in the west of the Longxi Loess Plateau, a terrain high in the west and south and low in the northeast. The Yellow River flows from the southwest to the northeast, cutting through the mountains and forming a beaded river valley with alternating valleys and basins. The soil in Lanzhou is mainly weathered red sandstone and loess soils, which are softened and disintegrated by water; further, their structures are easily disturbed and destroyed during foundation pit excavation. With the implementation of China’s national western development strategy, the impacts of human engineering activities on the geological environment have intensified and potential geological hazards are increasing. A new survey shows that geological hazards in Lanzhou are characterized by frequent activities and many hidden hazards. Several geological disasters occur in the areas surrounding Lanzhou city, which greatly impacts on urban construction and the safety of people and property. As they are affected by factors such as tectonic movement, concentrated rainfall and human activities, landslide geological disasters are frequent and widespread. Landslides are mostly loess landslides (Wang et al., 2020). The hidden threat of a landslide in the study area not only constitutes a serious threat to life, property and important infrastructure; It also threatens the region’s ecological security and sustainable development (He et al. 2020). Therefore, it is necessary to investigate and manage these hidden danger areas early.

Figure 1. Experimental scene.

2.2. Data sources

2.2.1. Time series InSAR deformation data

Landslides are the product of displacement. InSAR data can reflect the time series characteristics of any landslide type. This study selects interference width (IW) mode Sentinel-1A data from an open-source data platform provided by the European space agency. In this study, the Environment for Visualizing Images (ENVI) SARscape platform is used to process the Sentinel-1A data from 2015 to 2020 in Lanzhou. We utilize the SBAS-InSAR technology to obtain time-series InSAR deformation information (Berardino et al. 2002; He et al. 2021b) (Figure 2). Due to the influence of satellite orbits and data cycles, there is data missing in the upper left corner of the study area; this is displayed as blank in Figure 2. However, this does not affect our experimental testing and analysis. When selecting positive and negative experimental samples, we avoid this region so that all the features included in the positive and negative samples are representative and complete and do not affect the model’s training and evaluation results. The parameters and methods for each step in the SBAS-InSAR experimental process are given in Table 1.

Figure 2. Time-series InSAR cumulative deformation results.

Table 1. Key parameters and methods of the SBAS-InSAR process in the study.

Number	Processing flow	Key parameters and methods
1	Data import and tailoring	Vector boundary clipping of study area
2	Connection graph generation	Spatial baseline threshold: 25% Time baseline threshold: 360 days; 3D unwinding
3	Interference processing	Goldstein filter; Minimum cost flow unwinding Azimuth Looks:1; Range Looks:4
4	Track refining and releveling	SRTM DEM (30m)
5	Deformation information inversion	Singular value decomposition method; Model method to eliminate atmospheric phase; High-pass filtering to remove nonlinear deformation phase; Low pass filtering to remove incoherent noise
6	Geocoding	Rate accuracy threshold: 2, elevation accuracy threshold: 5

To verify the reliability of InSAR deformation data, areas with large surface deformation are selected from the SBAS-InSAR results for field investigation and UAV shooting (Figure 2). We find that the areas with large surface deformation are mainly due to fractures, subsidence and collapse. Road cracks are very common, with an approximate width of 5 to 10 cm and a length of several meters. The UAV aerial film reveals that the gully is collapsed. Its facade reaches several meters, and the collapse border has toppled electric poles. It is surrounded by slides of different sizes. These observations indicate that the InSAR results of this study are reliable.

2.2.2. Landslide cataloging

Based on historical landslide data, InSAR deformation information, Google earth images and field surveys, we visually interpret 160 landslides in this study, a total landslide area is 10.328 km². A total of 832 landslide sites are obtained (Figure 1). The historical landslide data comes from historical records, which we obtained from relevant departments, over the period from 1990 to 2020. The historical landslide data mainly consists of landslide events and field investigations by professional personnel. The data is reliable and highly precise through visual interpretation. Landslide types in the study area include loess landslides and rainfall landslides; we select landslide influencing factors that mainly affect these two landslide types.

2.2.3. Landslide influencing factors

Based on existing studies and the actual situation of Lanzhou (He et al. 2021a; 2023a), the principle of systematicity, representativeness, hierarchy and operability, and the multi-collinearity analysis method between landslide influencing factors (Chen 2016), we select 10 landslide influencing factors. These cover topography, geology, hydrology, vegetation cover and human engineering activities across various categories. The 10 landslide influencing factors include altitude, slope, aspect, lithology, distance from faults, distance from rivers, land use, distance from roads, the Normalized difference vegetation index (NDVI) and rainfall. The construction results of the landslide influencing factors are presented in Figure 3.

Figure 3. Spatial distribution of landslide influencing factors.

Of the factors listed above, altitude and slope can provide accurate potential energy for predicting landslide occurrence (Pham et al., 2017a). Aspect is linked to rainfall and solar radiation (Chauhan et al. 2010). Lithology impacts slope structure and the shear strength of the soil (He et al. 2021a). Distance type is calculated by Euclidean distance; distance from faults affects the mechanical structure of landslides, while distance from rivers affects surface runoff on runoff erosion and soil erosion. Land use and distance from roads reflect human's intervention (Pu 2022). NDVI is an index used to measure the growth status of surface vegetation. When the NDVI value is small, it indicates poor vegetation growth and poor soil stability, which may equate to high landslide risk (Sun et al. 2013 ). A considerable amount of rainfall will increase the volume of surface water and groundwater, change the physical properties of soil, increase the permeability of soil and aggravate the occurrence of landslides. The landslide influencing factor data and resolution information used in this study are provided in Table 2.

Table 2. Landslide influencing factors data.

Data	Resolution/Scales	Content	Landslide factors	Data sources	Period	Website
DEM	30m	Digital elevation model	Altitude, slope, and aspect	United States Geological Survey	2000	http://earthexplorer.usgs.gov
Land use	10m	Different types of land use	Land use	Tsinghua University	2017	http://data.ess.tsinghua.edu.cn/
Rainfall	30m	Daily cumulative rainfall	Annual cumulative rainfall	Geographic remote sensing ecological network	2020	http://www.gisrs.cn/
NDVI	30m	Normalized difference vegetation index	NDVI	National Ecological data center resource sharing service platform	2015–2020	http://www.nesdc.org.cn
Lithology and faults	1:4,000,000	Different types of lithology, fault network calculated from faults	Lithology, distance from faults	United States geological service	2000	https://www.cgs.gov.cn/
Rivers and roads	1:1000,000	River and road network calculated from rivers and roads	Distance from rivers, distance from roads	Lanzhou Natural Resources Bureau	2020	http://zrzyj.lanzhou.gov.cn/

3. Methods

Figure 4 charts our proposed unified deep learning framework for LSA. First, 10 landslide influencing factors are chosen by multicollinearity analysis. Second, we construct landslide inventories, combining the 10 landslide influencing factors and time-series InSAR deformation data to construct training, validation and testing datasets. Third, we develop an integrated neural network method to obtain LSA results. Finally, we use subjective results and quantitative indicators to evaluate the performance of the proposed method and analyze the relationship between LSA and landslide factors.

Figure 4. Flowchart of the proposed integrated neural network framework for LSA.

3.1. Integrated neural network method for LSA

3.1.1. Data set construction

In total, 24 stages of time-series InSAR cumulative deformation information are taken every 96 days per quarter. To meet the input and learning requirements of deep learning models, data preprocessing is required first. We use dynamic normalization processing for time-series InSAR cumulative deformation information to retain the feature distribution and time series relationship of the original data (He et al. 2023b). InSAR cumulative deformation extreme values are taken as the basis for data normalization. The cumulative deformation value in each node is dynamically stretched into the relative space at each time to ensure the overall cumulative deformation feature distribution and to not change the deformation interval of a single node. Finally, time-series InSAR deformation data is stretched to range of 0 to 255.

We then normalize the 10 landslide influencing factors. This is owing to the different representation content of each landslide influencing factor and the various representation intervals present in the raster diagram. Landslide influencing factors are then equalized by histograms to facilitate data analysis and ensure learning efficiency. We then divide the 10 landslide influencing factors into two types: continuous data types (altitude, slope, distance from faults, distance from rivers, distance from roads, NDVI and rainfall) and discrete data types (aspect, lithology and land use). For continuous data, the distinctive approach is to linearly scale the pre-processed raster images of each continuous landslide factor to the range of 0 to 255 according to their extreme values, then store the single-band bitmap. For discrete data, according to the number of types in their attributes, we use different quantiles valued between 0 and 255 and use different values to represent the different types. We maintain a specified gap between the values, so that the original features can be retained, and the data can be distributed in the same range, which helps facilitate the subsequent normalization processes. This method preserves the feature relationship of each landslide factor image in a single band, ensures learning accuracy and accelerates data reading and processing. We also apply a unified linear normalization scheme before entering the model to facilitate the subsequent model training and evaluation.

Landslide sample construction adopts the sliding window strategy to extract the landslide neighborhood, taking the landslide occurrence center as the sample extraction basis (Zhao et al. 2022). According to the data resolution and the size of the study area, we use a 13 × 13 window size for slide cutting. This not only generates data redundancy but also represents sample characteristics. Further, the moving step is 1, which facilitates pixel-by-pixel prediction. Using this method, we obtain 832 landslide points. During dataset construction, the landslide samples are enhanced with data, and the limited landslide positive samples and corresponding factor images are flipped horizontally and vertically to extend the data set. After data enhancement, 2 667 landslide samples are obtained. We determine the location of the positive sample based on the historical landslide point and establish a buffer area on the basis of this point. We then select non-landslide samples outside the buffer area to avoid the interference of the two types of samples and ensure the accuracy of the selection of samples. We randomly select the negative samples outside the buffer area to avoid subjectivity. The number of negative samples is the same as that of the positive samples to train the model results without preference. A total of 3000 non-landslide points is randomly selected. The dataset is then split into training, validation and testing sets per the ratio 8:1:1. We use the training set to learn landslide characteristics and train the internal parameters of the model. We used the validation set to optimize the model hyperparameters. We use the testing set to check the final generalization errors and evaluate the robustness of the model.

3.1.2. GRU model

The recurrent neural network (RNN) model is widely used for temporal feature learning. The long short-term memory (LSTM) and GRU networks are optimized variants of RNN (Wang et al., 2020). LSTM includes an input gate, a forgetting gate and an output gate. To simplify LSTM, GRU structure combines the forgetting gate and input gate into an update gate, and it uses a single gate control unit (Zhang, Guo, and Xiao 2021). The GRU network comprises only an update gate and a reset gate (Figure 5).

Figure 5. GRU network structure.

The specific contents of the function are: r is the reset gate; z is the update gate; h(t-1) represents the output state at the previous time; h ̃(t) shows the hidden status information; and h(t) represents the output value.

The single-directional GRU model only receives the sequence feature information from front to back according to the input order and can only serve before the current time step for modeling. Some critical information from the time-series data may be ignored, thereby deteriorating the model’s expressiveness. We propose a bidirectional gated recurrent unit (Bi-GRU) to solve this problem. Bi-GRU uses forward and backward information from time series data for modeling, thus improving the expressiveness and robustness of the model. In the Bi-GRU network’s architecture, the hidden structure is composed of two layers of GRUs with forward and reverse inputs. Further, its output is given by the state of the two GRU units at the current moment (Figure 6). This structure not only retains the characteristics of GRU unit structure but also improves the limitations of a single input sequence ensuring full temporal feature learning. Therefore, this study uses Bi-GRU to learn time series InSAR deformation features.

Figure 6. Bi-GRU network structure.

3.1.3. MSCNN model

Landslide influencing factor data is composed of the various factors that can induce landslides. Feature learning is performed by extracting spatial neighborhood features. CNN can extract to spatial information (He et al. 2021a), but traditional CNN takes the remaining network layer as the output by default without considering the features of different levels of semantic information. This makes it difficult to consider the low-level semantic characteristics of shallow features. Should the spatial scale of feature maps change due to pooling, information will be lost. To combat this, the present study proposes using MSCNN to extract the spatial features from the landslide influencing factor dataset, construct multi-depth fusion feature maps and establish the feature learning process at different scales.

MSCNN uses multi-scale technology to add full connections after each pooling layer (Zhao et al. 2022). Figure 7 illustrates the MSCNN network structure. The ‘flattening’ layer is used to compress multidimensional features into one-dimensional vectors. The feature results extracted from three fully connected layers are then fused and activated through the ‘Concat’ layer to establish a learning process that considers shallow and deep features. The calculation is as follows: (1) $h_i=f(W^i{x}^i+b^i)$(1) (2) $h={\sum }_{i=1}^n{h}^i$(2)

Figure 7. MSCNN network structure.

Where $h^i$ represents the i th fully connected layer, $x^i$ represents the i th layer expansion feature, $W^i$ and $b^i$ represent weight matrix and bias, f(·) is the fully connected layer construction function. All fully connected layers are fused together to obtain h is a multi-scale feature fusion result.

3.1.4. Integrated neural network model for LSA based on Bi-GRU and MSCNN

Based on the advantages of Bi-GRU and MSCNN in learning temporal and spatial features, this study proposes an LSA method that uses the Bi-GRU and MSCNN models to extract temporal and spatial deep features from InSAR deformation information and the landslide influencing factor data. The restricted model construction ideas are shown in Figure 8. The network model learns the temporal features and spatial features about landslide occurrence in parallel and combines the two features to generate LSA.

Figure 8. The integrated neural network structure for LSA.

For deformation feature learning, InSAR deformation information images obtained with a fixed window size centered on the target point are sorted by time series. That is, the cumulative deformation of the point is displayed from the front to back over time. The sorted whole is then taken as the input to the deformation feature learning module. Here, the spatial information is aggregated into the time-dimension by time distribution CNN (TD-CNN), and the time series relationship is learned by Bi-GRU. The construction of the temporal dimension relies on the time distributed packaging layer under Keras (Luna-Alvarez et al. 2020), which stores temporally ordered deformation data in a common layer without affecting the CNN's convolution and pooling for spatial dimension feature aggregation. Using TD-CNN, the spatial deformation information at each time is aggregated to each time step in tensor form, creating a feature vector with both spatial features and temporal relationships. This vector is then passed to Bi-GRU to learn the overall deformation features. Bi-GRU accepts feature vectors with a step size of 1 and a total time dimension of 24, where each time step represents a fixed moment of deformation, and the vector represents the cumulative deformation relationship. The bidirectional learning process makes Bi-GRU learn not only changes of information from the forward time series, but also reverse learning to comprehensively summarize the entire change process.

To perform influencing factor feature learning, the window features of the same location are used as the input for this model. Here, we directly combine different influencing factors in the channel dimension through a ‘Concat’ operation to simultaneously learn multi-dimensional features. The spatial features of a landslide or the absence of one, and the interconnections between different influencing factors, are learned through multiple convolutional kernels. The feature scale is reduced through pooling layers. The more generalized and abstract features are gradually formed with the depth of the network. Multi-scale fusion is added to the simple CNN structure to fuse and splice the feature vectors formed by the feature maps under different depths. This ensures that the final obtained features have both specific neighborhood spatial relationships with the influencing factors and generalized high-dimensional features about the likelihood of a landslide to achieve accurate judgment.

Finally, the learning results of both types of features are output in the form of feature vectors. The two vectors are connected through the ‘Concat’. They are then integrated and mapped through two fully connected layers to achieve comprehensive feature utilization. The final output of the model is a single neuron, which represents the evaluation probability of landslide occurrence.

3.1.5. Model training

The proposed network model is implemented in Python under the TensorFlow framework, and the results are obtained on a host equipped with an Intel (R) Xeno (R) Silver 4214 processor and an NVIDIA Quadro P2200 graphics card. The deep learning experiment has high requirements for hardware; Table 3 shows their basic configuration. The configuration of important software is shown in Table 4. The MSCNN and the proposed model are trained for 64 epochs, taking 1 h, 34 min and 3 s, and 1 h, 40 min and 20 s, respectively.

Table 3. The basic system platform configuration.

Project	Operating system	CPU	Memory	Hard disk	GPU
Content	Microsoft Windows 10 Home version	Intel(R) Xeon(R) Silver 4214 CPU @2.20GHz	16G 2400MHz	TOSHIBA DT01ACA200 1.8TB	NVIDIA Quadro P2200

Table 4. The important software configuration.

Project	Graphic driver	CUDA version	Python	Keras	Tensorflow
Content	472.39	V11.4	3.6.13	V2.10.0	V2.3.0

The trial-and-error method is used in this study to adjust the model parameters several times before compiling the results. The exact parameter design of the proposed model is presented in Figure 8. Most of the network layer activation functions in the proposed model are expressed by a rectified linear unit (ReLU), which aims to avoid the insufficient response ability of simple linear relations to complex variables and uses nonlinear functions to approximate the output results and real data. The ReLU also combats the vanishing gradient problem and improves the efficiency of model training. We employ a Sigmoid function as the output activation function in the output layer, then quantify the probability of landslide susceptibility on a scale from 0 to 1.

We also add a BN layer after two convolutions of each feature graph to speed up training. Dropout operations are included in the hidden layers of the network to randomly drop a fraction of the neural units. This is to prevent the model from overfitting to the training data, thus enhancing the model’s generalization ability. For network optimization, we use an Adam optimizer to adaptively update the network weights to obtain more robust model results. We set the initial learning rate of the model to 0.001, then in the iteration process, the learning rate dynamically adjusts by observing the loss value of the validation set through the callback function. Specifically, when the performance is equal to every three iterations, the learning rate is attenuated by a multiple of 0.5 to better balance the convergence speed and quality.

For the output of LSA results, each value represents the probability of landslide occurrence at different locations. We treat this task, which involves mapping various features to a probability distribution range, as a regression problem. The choice of loss function significantly affects the model learning trend. For regression problems, the mean square error (MSE) is a commonly used loss function that provides a convenient way to measure the ‘mean error’ of the result from its true value. Therefore, we also use MSE as a loss function to measure the disparity between the predicted landslide probabilities and the actual existence of landslides. Smaller MSE values indicate more accurate predictions by the model.

3.2. Evaluation methods

3.2.1. Evaluation of Landslide influencing factors

3.2.1.1. Multicollinearity analysis

The multicollinearity of landslide factors significantly impacts the model, so it is necessary to conduct a multicollinearity analysis among the selected landslide factors. We use a variance inflation factor (VIF) and tolerance (T) to conduct collinearity tests on the 10 landslide influencing factors (He et al. 2023a). (3) $VIF={\displaystyle \frac{1}{1-A^2}={\displaystyle \frac{1}{T}}}$(3) Where, $A^2$ represents the complex correlation coefficient. VIF is greater than 10 and its reciprocal T is less than 0.1, which indicates that there is a multicollinearity problem.

3.2.1.2. Importance ranking of landslide factors

The degree to which each factor influences a landslide varies and thus requires further analysis. RF algorithm is commonly used in feature screening and other fields. RF can give a variable importance score (Archer and Kimes 2008). In this paper, we use the out-of-bag (OOB) error rate of RF to calculate the relative importance scores of different landslide factors.

3.2.1.3. Frequency ratio

Frequency ratio (FR) is expressed by calculating the probability of landslide occurrence according to different factor classification intervals. Here, the magnitude is equal to the ratio of the landslide area to the classification area ratio (He et al. 2022b). However, because historical landslides have a certain evolution progress, the landslide scope and actual area differ. This means it is not easy to accurately calculate the area of each landslide. Further, if there are many landslides, it will require significant work to calculate the landslide area. Therefore, we use relative frequency instead of the landslide area ratio, to better reveal the spatial distribution relationship between the landslide data and influencing factors. (4) $FR={\displaystyle \frac{N_{ij}}{N_r}/{\displaystyle \frac{A_{ij}}{A_r}}}$(4) Where $N_{ij}$ is the landslide area corresponding to the i landslide factor of the j th species, and $N_r$ is the whole landslide area. $A_{ij}$ is the region corresponding to the i th landslide factor of the j th species,$A_r$ is the whole study area.

3.2.2. Model accuracy analysis

The model must be able to accurately identify of landslide samples to ensure the reliability of the final identification results. To quantitatively evaluate the landslide identification results, we use common indexes from semantic segmentation problems. We use the confusion matrix, a basic and widely used method (He et al. 2021c), as illustrated in Figure 9.

Figure 9. Confusion matrix.

Here, TP indicates that the model correctly judged the landslide target; FN indicates that the landslide target was incorrectly judged as a non-landslide; FP indicates that the non-landslide target was incorrectly judged as a landslide; and TN indicates that the model correctly judged the non-landslide target. Using different data organizations for the four ratios, more indicators can be constructed to determine the model’s performance and accuracy. These indicators are described below.

Accuracy (ACC) describes the classification accuracy of a classifier. (5) $ACC={\displaystyle \frac{TP+TN}{TP+FP+FN+TN}}$(5) Precision measures the ratio between the correct target and the judged correct target; that is, the proportion of positive (true label is positive) samples among all the samples judged by the model as positive examples. (6) $\mathrm{Precision}={\displaystyle \frac{TP}{TP+FP}}$(6)

Recall measures the proportion of targets judged to be correct among the correct targets. (7) $\mathrm{Recall}={\displaystyle \frac{TP}{TP+FN}}$(7) F1-score is the weighted and harmonic average of the model’s accuracy and recall rates. Considering the importance of accuracy and recall, this F1-score is a good indicator of the model’s comprehensive performance, the greater the F1-score value, the better the model’s performance. (8) $F1\mathrm{\_}\mathrm{score}=2\times {\displaystyle \frac{\mathrm{Precision}\times \mathrm{Recall}}{\mathrm{Precision}\times \mathrm{Recall}}}$(8)

True negative rate (TNR) refers to specificity, describing the proportion of negative cases identified among all negative cases. (9) $\mathrm{Specificity}={\displaystyle \frac{TN}{FP+TN}}$(9) Finally, the Kappa coefficient checks consistency; specifically, whether the model’s judgment results are consistent with the actual classification results (He et al. 2021a). (10) $\mathrm{kappa}={\displaystyle \frac{N\sum _{i=1}^n{X}_{ii}-\sum _{i=1}^n(X_{i+}{X}_{+i})}{N^2-\sum _{i=1}^n(X_{i+}{X}_{+i})}}$(10)

Where, n is the total number of categories of the classified sample, that is, the sum of the columns of the confusion matrix; N represents the total number of samples; $X_{ii}$ is the number of samples in the $i^{th}$ row and $i^{th}$column, representing the samples that are correctly classified; $X_{i+}$ and $X_{+i}$ are the total number of samples in $i^{th}$ row and $i^{th}$column, respectively.

4. Results

4.1. Landslide influencing factor analysis

In this paper, the multicollinearity values of 10 landslide influencing factors are calculated using VIF and T in SPSS 26 software. The multicollinearity analysis results are presented in Figure 10. The VIF values of the 10 selected landslide influencing factors all exceed 1 and are less than 3; of these, the VIF value (2.189) for Distance from rivers is the highest. The lowest value is annual cumulative rainfall (1.036). The T value ranges from 0.4–1. The VIF and T value indicate that the 10 landslide influencing factors are independent of each other and have no collinearity problem among them. This indicates that the 10 landslide influencing factors examined in this study may be used in the model for learning and evaluation to ensure the effectiveness and accuracy of the data.

Figure 10. VIF and T values of landslide influencing factors.

The impact degree of each landslide factor varies, so we examine each further. The RF algorithm is widely used in feature screening and other fields. RF can provide a variable importance score while analyzing data. We use the OOB error rate of RF to calculate importance scores of different landslide influencing factors relative to the occurrence of landslides. These results are shown from high to low in Figure 11. Here, Distance from roads is valued as the most important, revealing the serious influence of human engineering activities on the occurrence of landslides in the Lanzhou urban area. The next most important are altitude and Distance from faults. The study area has considerable relief and a complex environment, making altitude and Distance from faults its the main landslide influencing factors. Rivers and rainfall similarly promote the occurrence of landslides by affecting soil structure. The study area is mainly loess, which is easily affected by rainfall. The other factors influence the occurrence of landslides in Lanzhou to varying degrees.

Figure 11. Ranking results of the 10 landslide influencing factors.

4.2. Model accuracy analysis

We create a comparative experiment. The objective of this comparison is to evaluate landslide susceptibility using the MSCNN network model and the proposed model. To control the variables, we verify that the data input forms and hyperparameters of MSCNN and the proposed model are consistent. Data from the validation and testing sets are compared and discussed.

Changes in MSE loss and ACC during the training process for the validation set reveal the learning effect and generalization ability of both models. Therefore, this study first uses the training process curve of the validation dataset for analysis (Figure 12), produced by Origin software. The training curves from the two models show that the gradual decline of MSE values and the gradual increase of ACC values reflect the learning abilities of the models. Here, the two models adopt the error back propagation method of the gradient descent algorithm as the training rule. The weight between neurons is adapted according to the error between the actual and predicted values of the training samples in the iterative process of gradient descent. The curve that eventually flattens shows that the two models both form a stable structure with decision-making ability.

Figure 12. Comparison of training process of validation set (a) MSCNN, (b) the proposed model.

In the training process, the MSCNN model exhibits a large oscillation in the initial iteration number. This suggests that the MSCNN model is sensitive to data anomalies in the beginning. It slightly overfits early features, but the two values return to their normal track later on, indicating that the MSCNN model has a robust learning ability. In contrast, both the MSE and ACC values of the proposed model stay with a normal change trend, indicating that the proposed model is not disrupted by accepted complex features. Further, the proposed model has a stable learning state for landslide factors and InSAR deformation data features and gradually increases. The proposed model becomes stable after a little more iteration, indicating that the proposed model parameters must be adjusted according to complex features and need a certain amount of learning time. The ACC value of the proposed model is 0.825, which is 1.2% higher than that of MSCNN (0.813). The MSE value is 0.065, which is 0.4% lower than that of MSCNN (0.069). The proposed model can thus better grasp landslide characteristics than the MSCNN model, and it produces better results.

The testing set data behave strangely in the model, showing the learning effect and generalization ability of the model more objectively. The testing set test results compare the predicted values and true values. Since the output result of the model is a probability distribution ranging from 0 to 1, the prediction results are divided into ‘landslide’ and ‘non-landslide’ with a threshold value of 0.5. Figure 13 shows the division results of the confusion matrix obtained by the two models.

Figure 13. Comparison of confusion matrix results of test set (a) MSCNN, (b) the proposed model.

In the confusion matrix, both MSCNN and the proposed model have a high true-positive rate (TPR) and true-negative rate (TNR) in the testing set, indicating the reliability of both models in the face of unfamiliar data classification. The proposed method has good classification results for both landslide and non-landslide data in the testing set, which further proves the learning ability of the proposed model (Figure 13(b)). Combined with the confusion matrix results of the two models, the classification effect of the true value of landslide is better than that of the true value of non-landslides, which proves that the proposed model’s learning ability of landslide features and performance. Per the confusion matrix results, there are some classification errors in both models, revealing areas of high landslide occurrence in the study area. However, the accurate evaluation of landslide susceptibility should be built on the accurate learning of landslide characteristics. The proposed method performs better with landslide samples and offers reasonable classification results in non-landslide areas. The corresponding indexes are calculated according to the resulting confusion matrix. Figure 14 shows the calculation results. The proposed method achieves better results in all aspects compared to the MSCNN method.

Figure 14. Comparison of evaluation indicators of both models.

In this paper, we generate the confusion matrix by taking the midpoint of the model prediction value (0.5) as a preset threshold, which is a typical classification evaluation method for landslide problems. However, the problem with LSA is not a simple binary classification problem. When the threshold changes, the TPR and FPR in the matrix may change, and the ratio under the variable threshold can be matched one-to-one. The receiver operating characteristic curve (ROC) can then be plotted. The ROC reflects continuous changes in data specificity and sensitivity and is commonly used in LSA studies (Pham et al., 2017b; Wang et al. 2015). The area under the curve (AUC) can directly reflect these results. Here, we calculated the ROC of the two models in the testing set. As shown in Figure 15, the AUC value of the proposed method reaches 0.9089, which is better than that of MSCNN model.

Figure 15. Comparison of AUC of both models.

We used the entire dataset (training, validation and testing sets) to calculate the MSE and RMSE values between the predicted values and true values (Figure 16). The MSCNN model and the proposed model obtained MSE values of 0.0487 and 0.0437, respectively. The proposed model thus offers the smallest difference between the predicted values and true values, further indicating that its superior predictive performance.

Figure 16. Errors of between the predicted value and the true value in training data (a) MSCNN, (b) the proposed model.

These results show that the proposed method exhibits excellent performance in all aspects compared to the traditional MSCNN method. The proposed model not only learns the complex temporal and spatial characteristics of landslides but also has a stable learning process and high robustness. Further, its performance in the testing set proves its generalization ability. We believe that the proposed model is suitable for evaluating of landslide susceptibility and can obtain more accurate and reliable results than presently available models.

4.3. LSA analysis

LSA is performed by re-assigning pixels using the trained model in ArcGIS 10.5. LSA results are divided into five levels by natural breakpoint (He et al. 2021a) using ArcGIS 10.5: very low, low, medium, high, and very high susceptibility. LSA results are presented in Figure 17. Figure 17(a) represents the LSA obtained by the MSCNN model. Figure 17(b) is the LSA obtained by the proposed model.

Figure 17. Results of LSA (a) MSCNN, (b) the proposed model.

As shown in Figure 17, the results of the two models differ but share some commonalities. The prediction locations of two methods are roughly the same for very high susceptibility areas, and most of the historical landslide points are recorded as having very high susceptibility, which proves the reliability of the results. More than half the areas predicted by the two methods are identified as very low susceptibility. Very high susceptibility is distributed along the mountain trend in the northern and southern mountains, which supports the actual events in the studied area. The proposed method has fewer very high susceptibility areas and exhibits more accurate location characterization.

We compare historical landslide distribution with our LSA results. We then evaluate our results objectively per the ratio of the historical landslide area to areas of each susceptibility level. The results are presented in Figure 18. In the LSA, previous studies show that most landslides are distributed in high and very high susceptibility areas, whereas low susceptibility areas have fewer landslides (Chauhan et al. 2010). In this study, the determined high and very high landslide susceptibility areas contain more than 80 percent of historical landslides, indicating that both models effectively evaluate landslide susceptibility in Lanzhou. Huang et al. (2020) show that high and very high susceptibility areas cover a small part of the area in the LSA. However, our results from the MSCNN model diverge from this, as the ratios of the high and very high susceptibility areas are higher and there exists a degree of landslide misjudgment. In contrast, the results from the proposed method are more reliable. The ratios of the high and very high susceptibility areas are lower, and the model correctly puts the largest historical landslides ratio in the very high susceptibility area. Here, only 10.18% of the study area accurately accounts for 84.79% of the historical landslides, indicating that the proposed model has strong feature learning ability. Therefore, compared with the MSCNN model, our proposed model exhibits better precision and offers LSA results with more important reference value.

Figure 18. Percentage of landslides and LSA of each category in two models (VL: very low, L: low, M: moderate, H: high, VH: very high).

To further understand the landslide development law in the study area, we use FR to analyze the distribution characteristics and environmental factors of the high susceptibility areas. FR represents the possible development degree of landslides in different intervals of different landslide factors, as well as the development relationship between the different types of a single landslide factor and landslides. According to the FR results (Figure 19), more than half of the areas in Lanzhou are below 2 000m above sea level, and most of the areas with high landslide susceptibility are distributed. Slope provides the necessary air conditions for landslide occurrence, and landslide occurrence gradually increases with the increase of slope. The maximum FR value (1.94) of slope in this study area is between 40° and 50°. The landslide susceptibility areas are allocated in all aspects without obvious regularity, and east aspect has the highest FR value (1.27). Among land use types, grasslands are the most susceptible to landslide occurrence, with forests exhibiting the lowest FR value (0.17). This indicates that high vegetation coverage has a definite inhibitory effect on landslide occurrence. Second is the impervious surface FR value (1.10), indicating that urban construction activities have a significant influence on landslide occurrence. Middle and Lower Pleistocene ice has the greatest influence on the occurrence of landslides, the lithology of which is released in the southern mountain areas of the Chengguan district. Distance from roads, faults and rivers all negatively correlate with the occurrence of landslides. The longer the distance, the smaller the FR value, and the smaller the contribution to landslides, which proves that landslide progress is indeed affected by human activities, active faults and soil erosion. More than 80 percent of the areas in Lanzhou have an average annual rainfall of less than 20 mm, and the areas with an annual rainfall of 220–240 mm have the highest FR value (2.00). The observed NDVI value further indicates that most landslides are distributed in areas with poor vegetation growth; the FR value (1.08) is highest here between −0.5 and 1.

Figure 19. Relationship between high susceptibility areas and environmental factors.

5. Discussion

5.1. Performance of proposed model frame

LSA predicts the probability of potential landslide occurrence, which is of profound significance to the prevention and control of landslide disasters (Reichenbach et al. 2018). Therefore, a reliable LSA model is of great importance (Zhu et al. 2020). The neural network method can accurately identify the implicit relationship between landslide influencing factors and landslide susceptibility (Huang et al. 2020), and it has achieved good results for LSA (He et al. 2021a; Zhao et al. 2022; Fang et al. 2020). However, due to the constantly changing nature of landslides, the time series dynamic features of the landslide development process have yet to be fully understood. This limits the accuracy of LSA and results in low LSA reliability.

Time series InSAR is commonly used to detect the slow movement of landslides. This data provides some time series deformation features (Dai et al. 2022; He et al., 2023b). Previous studies have shown that taking InSAR deformation features into account can improve the reliability of a landslide susceptibility prediction model (Yuan and Chen 2022; Yao et al. 2023; Gao et al., 2023b). However, studies to data use only a single image of an InSAR deformation rate as the feature input (Meghanadh et al. 2022), ignoring the multiple image features of the time series InSAR deformation process. In addition, the existing machine learning model introduces InSAR deformation information and landslide influencing factors together for training. Meanwhile, deep temporal features are difficult to mine and utilize so they have not been applied much. Further, our research proves that when temporal and spatial influencing factors are introduced to a model for training at the same time, the subsequent time series feature learning is insufficient, which produces in insufficient prediction accuracy (Figure 17(a)). Our aim, therefore, is to establish an integrated neural network combined model that learns InSAR deformation information temporal features and landslide influencing factor spatial features separately, then merges the temporal and spatial features to improve the model’s prediction accuracy of landslide susceptibility. Our results confirm this (Figures 16 and 17).

In our proposed model, we first extract InSAR deformation temporal features and landslide influencing factor spatial features, then fuse the high-dimensional features into deep learning classifiers to predict landslide susceptibility. This helps the model to focus on temporal and spatial features in different directions and different dimensions (Ji et al. 2020; He et al. 2021c; Wei et al. 2022; Chen et al. 2023). In our experiment, the temporal features of InSAR deformation are extracted by TD-CNN and multi-layer Bi-GRU, and the spatial features of the landslide influencing factors are extracted by MSCNN. As the network deepens layer by layer, the final learning vector feature combines with the deep layer, shallow layer and time series InSAR deformation feature. Learning is ‘Concat’ and merges in the channel dimension. The fusion result is then learned by the fully connected layer and output as the final LSA result. A visualization of this temporal and spatial feature extraction process is shown in Figure 20. Demonstrably, ‘landslide’ and ‘non-landslide’ areas can be effectively distinguished by integrating temporal and spatial features. The results of the LSA in this study are reliable and interpretable.

Figure 20. Feature extraction process and high-level feature vectors.

5.2. Comparative analysis with previous models

Machine learning has been widely used in the LSA, especially in support vector machine (SVM) (Pourghasemi and Kerle 2016) and multi-layer perceptron (MLP) (Yi et al. 2020) activities. To further explore the reliability of the deep learning model proposed in this paper, we compare our results with SVM and MLP results. In this experiment, temporal InSAR deformation information data and 10 landslide influencing factors are simultaneously input into SVM and MLP models for training, and LSA are obtained. Natural breakpoint method by ArcGIS 10.5 software is used to generate the LSA results of each level (Figure 21). Our results indicate a consistent trend between the LSA results obtained by the two classical machine learning models and the LSA results obtained by our proposed model (Figure 17). However, there are large differences in the susceptibility levels predicted. The results predicted by the SVM and MLP models account for a higher proportion of very high susceptibility, while the results obtained by our model account for a lower proportion of very high susceptibility. Previous studies show that the proportion of high susceptibility should be small in the LSA. Therefore, we can conclude that the LSA obtained by our proposed model is more reliable.

Figure 21. Results of LSA (a) SVM, (b) MLP.

The middle sections of the study area are urban, but the SVM and MLP models predict a moderate susceptibility level. This is inconsistent with the real value. The SVM and MLP models struggle to learn temporal features, because they input InSAR deformation data at the same time. This further proves that the simultaneous input of multiple time-series InSAR deformation images and landslide influencing factor data into existing models negatively affects accuracy. The framework proposed in this study is more reasonable, as the time-series InSAR deformation data and landslide influencing factor data features are learned separately. Only after this process are the space-time features fused. These results further demonstrate the superiority and reliability of the proposed model in this paper.

5.3. Limitations and future work

This study constructs a unified deep learning LSA model that integrates temporal and spatial features and InSAR deformation information. It benefits the LSA and provides data and technical support for disaster reduction and prevention in Lanzhou city. However, there are some limitations:

1. The InSAR technology used in this study is incoherent in some areas with dense vegetation or excessive deformation, leading to the absence of or inaccurate deformation value in some parts of the study area. This causes such features to be omitted from the analysis.

2. We use only the most precise surface deformation information as the temporal dynamic features of a landslide. In the actual development process of a landslide development, other factors not include in this study have temporal features as well; these include soil moisture, surface vegetation coverage and vegetation cover type. Subsequent experiments should make comprehensive use of multi-source dynamic features to achieve a fuller understanding of landslide characteristics and obtain more accurate LSA results.

5.4. Outlook

Based on the above analysis of high landslide susceptibility areas and environmental factors, landslides in Lanzhou are mainly distributed in areas with complex geological structures, elevated soil moisture content and frequent human engineering activities. From the perspective of landslide disaster prevention, this paper puts forward the following suggestions for avoiding landslide disaster losses:

1. Increase vegetation coverage: An appropriate increase in vegetation coverage effectively reduces soil erosion and surface erosion and improves soil stability and skid resistance. The government of Lanzhou should increase the level of vegetation coverage through afforestation, lawn construction and ecological engineering promotion; doing so will promote the improvement of the geological environment and ecological protection and reduce geological disaster occurrence.

2. Control of human engineering activities: Excessive human engineering activities destroy the integrity and stability of the land and exacerbate the occurrence of landslides. Due to the unique topography of Lanzhou, the relationship between man and land is complicated. While pursuing economic development, the government should strengthen its supervision of key human activities such as mining, construction and cutting mountains and land to control damage to the land.

3. Strengthen geological disaster monitoring: A sound geological disaster monitoring and early warning system should be established to actively identify geological disaster safety hazards in key districts and counties.

6. Conclusion

Based on a 24-scene time series InSAR cumulative deformation and 10 landslide influencing factors, we construct a unified deep learning LSA framework. The time-series InSAR cumulative deformation features are learned through TD-CNN and Bi-GRU, and landslide influencing factor spatial features are learned through MSCNN. Temporal and spatial multi-scale features are fused through channel fusion. An MSCNN model is also set up to compare and assess the validity of our proposed model. The ACC value of our proposed model increased by 1.2% and its MSE value decreased by 0.4% compared with that of the MSCNN model. This indicates that our proposed model better grasps landslide characteristics, and its learning process is stable and robust. On the testing set, the proposed model is larger than the MSCNN model in every index. Further, its AUC value reaches 0.9089, indicating that our proposed method has good generalization ability. Our proposed method accurately describes high susceptibility areas and accurately accounts for 84.79% of the historical landslides in 10.18% of the study area, demonstrating its high feature learning ability. According to these comprehensive results, the proposed model is an effective and feasible LSA method that can aid landslide representation learning and obtain accurate results. Our analysis of the LSA results reveals that the landslides in Lanzhou city mainly occur in areas with complex geological structures, high soil-water content and frequent human engineering activities.

Acknowledgments

We would like to express our great appreciation to the editors and three anonymous reviewers for constructive comments that helped improve the manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement (DAS)

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Additional information

Funding

This work was supported by the National Natural Science Foundation of China [grant number 42201459; 42261076], Science and Technology Project of Gansu Province [grant number 23JRRA881] and Science and Technology Research and Development Plan of China State Railway Group Co., Ltd. (grant number P2021G047).