A landslide susceptibility assessment method based on auto-encoder improved deep belief network

2.1 Study area

Zhouqu County is located south Gansu Province, the middle reaches of the Bailong River basin, and includes the Qinghai–Tibet Plateau, Loess Plateau, and Sichuan Basin [27], and the eastern edge of the Min and Di Mountain system of the West Qinling Mountains and the Qinghai–Tibet Plateau. The county is one of the most serious areas of geological hazards in China [28]. Zhouqu County is adjacent to Wudu County in the east, DiBu County in the west, Jiuzhaigou County in the south, and Tangchang County in the north. Zhouqu County is about 99.4 km long from east to west and 88.8 km wide from north to south, with a total area of 3,010 km². The geographical location is shown in Figure 1. The topography of Zhouqu County is high in the northwest and low in the southeast, with a difference in elevation reaching 3,305 m. Rivers running through the county include the Bailong, Gongba, and Boyu, distributed in a dendritic pattern with tributaries such as Nanyu Ditch and Quwa Ditch. Fault activity is widely distributed, providing favorable conditions for landslides. Heavy rainfall is more frequent in summer, causing sudden groundwater level rise, washing away farmland, causing soil erosion, and triggering geological disasters, such as landslides and mudslides. The vigorous construction of roads in recent years and the large amount of earth and rock excavations on the slopes have significantly increased the probability of various geological hazards in Zhouqu County.

Figure 1

Location of the study area.

2.2 Historical landslide dataset

Landslide identification is the basis of landslide susceptibility assessment. The dense vegetation on the surface of Zhouqu County makes it challenging to identify historical landslides using optical image texture features. The small baseline subset interferometric synthetic aperture radar (SBAS-InSAR) technique is widely used for landslide identification [29,30,31 32,]. Therefore, this study is based on European Space Agency Copernicus Sentinel-1A satellite SAR images and uses SBAS-InSAR technology combined with Google Earth imagery to perform historical landslide identification. In this article, we obtained surface deformation information from Zhouqu County based on SBAS-InSAR technology. The density of SBAS-InSAR coherent points in the region was calculated to be 347/km², which proved the reliability of the deformation results obtained using SBAS-InSAR technology. Combined SBAS-InSAR deformation results, topography, and Google Earth images for historical landslide identification. The specific rules are as follows:

1. Deformation points with a deformation rate greater than 10 mm/year or a deformation rate less than −10 mm/year.

2. Deformation points with slope values greater than 5°.

3. Deformation points that are spatially distributed in aggregation.

The points satisfying the above criteria were loaded onto the Google image. The first to be determined was whether the aggregated points had a circle chair, ellipse, rectangular shape, and irregular polygon shape holistically from the aggregated points. If so, we made a local judgement, observing whether there were features such as a landslide body, a landslide back wall, and a landslide boundary.

Using the method above, 84,128 deformation points were screened, 188 landslides were identified, with areas ranging from 0.001 to 4.194 km², and a total area of 40.32 km², as shown in Figure 2. Most of these landslides were distributed along river valleys, and most were crossed by faults. As a result, gully erosion and tectonic activities could significantly influence landslides in the area; a small number of landslides were distributed on roads, demonstrating that human activities also contribute to landslides in the area. Comparing the identified landslides with the catalogued landslides in the 2003 geological hazard map of Zhouqu County, 44 landslides with existing data were identified in this study, accounting for 80% of the total number of landslides in the 2003 geological hazard map, indicating that the landslide data identified in this study are credible.

Figure 2

Location distribution of landslides in Zhouqu County: (a) Zhongpai landslide, (b) Suoertou landslide, and (c) discharge slope landslide.

In addition, there are three landslides that have been in news because of their huge harm in study area, named Zhongpai landslide, Suoertou landslide and Discharge slope landslide, respectively. Taking into account time, money, and the pandemic, we select the Google image of the three landslides as a field view (Figure 2(a–c)).

2.3 Landslide factor dataset

2.3.2 Landslide factor screening

Numerous factors affecting landslide instability, among which there are inevitably large co-linear trends, which will increase the instability of subsequent landslide susceptibility assessment models and reduce their prediction accuracy [34]. Therefore, the factors influencing landslides with large co-linear trends need to be removed. The stepwise regression method allows the independent variables to be entered sequentially into a specific method according to the degree of influence of the independent variables on the dependent variable. A test was performed for each independent variable entered so that the final filtered independent variables had no significant linear correlation and had the maximum degree of influence on the dependent variable [35]. In this study, based on the order of importance of each landslide influence factor by geographic detectors [36], each factor was input into logistic regression, SVM, random forest, and MLP models in turn. The area under the curve (AUC) value obtained for each factor was counted, and the landslide factor corresponding to the AUC value after it started to decrease was removed.

The ranking of the importance of the 17 landslide impact factors in Zhouqu County based on geographic detectors is shown in Figure 3. The order is distance to road (q = 0.3447) > elevation (q = 0.29) > distance to river (q = 0.2739) > land use type (q = 0.2628) > MNDWI (q = 0.1701) > NDVI (q = 0.1289) > lithology (q = 0.1110); NDBI (q = 0.0843) > Distance to fault (q = 0.0812) > Aspect (q = 0.0513) > Slope (q = 0.0111) > SDS (q = 0.0109) > TWI (q = 0.0099) > RDLS (q = 0.0085) > Profile curvature (q = 0.0066) > curvature (q = 0.0046) > planar curvature (q = 0.0004).

Figure 3

Ranking of importance of landslide influence factors.

The statistical plots of AUC values for each method of logistic regression, SVM, random forest, and MLP are shown in Figure 4. Both the logistic regression and SVM models had the highest AUC values of 0.889 and 0.92, respectively, after adding the 14th landslide influence factor (Figure 4a and b). The random forest model had the highest AUC value of 0.951 after inputting the 10th landslide influence factor (Figure 4c). The MLP reached the highest AUC value of 0.924 after inputting the 9th landslide influence factor (Figure 4d). Combining the results of the four models, nine landslide impact factors were screened out: distance to road, elevation, distance to river, land use classification, MNDWI, NDVI, lithology, NDBI, and distance to fault for subsequent landslide susceptibility assessment in Zhouqu County. The spatial distribution of the nine landslide factors is shown in Figure 5.

Figure 4

AUC values of each model in landslide factor screening.

Figure 5

Spatial distribution of landslide impact factors after screening. In (g), APC, CSR, IPR, MSR, and SSR represent acid plutonic rock, sedimentary carbonate rock, intermediate plutonic rock, mixed sedimentary rock, and siliciclastic sedimentary rock, respectively.

3.1 AE

AE is a neural network with the same input and learning objectives, and its structure is divided into two parts: encoder and decoder. The structure of AE is shown in Figure 6: given an unlabeled input dataset { x n } n = 1 N , where x n ∈ R m × 1 , h n represents the hidden encoder vector computed from x n , x ∼ n represents the decoder vector of the output layer, then the encoding process of AE can be expressed as:

where f is the encoding function, w ₁ is the weight matrix of the encoder, and b ₁ is the bias vector of the encoder.

Figure 6

AE model structure diagram.

The decoding process can be expressed as follows:

where g is the decoding function, w ₂ is the weight matrix of the decoder, and b ₂ is the bias vector of the decoder.

The parameter set of the AEs is optimized to minimize reconstruction error:

In the formula, L represents the loss function, and the calculation formula is as follows:

3.2 Construction of an improved DBN network model based on AE

The basic component of the DBN network is the RBM model. However, the training of the RBM relies on the edge probability and conditional probability of the visual layer, which makes the RBM training process difficult [25]. In the AE training process, compared with the loss function of RBM, which is probability-based maximum likelihood, its partial derivative is easier to calculate directly [26]. Therefore, this research will use the AE model to replace the RBM model in the DBN network, complete the model improvement by modifying layers and adjusting parameters in AE model, and add label data and classifier to the hidden layer of AE at the top of the model to complete the construction of the landslide susceptibility assessment model based on improved DBN. The structure of the model is shown in Figure 7. The structure and hyperparameters of the model are shown in Table 1.

Figure 7

Structure diagram of the improved DBN landslide susceptibility assessment model.

In this study, a four-layer DBN is constructed, including three AE and an output layer to train the existing data. The entire training process was divided into two steps:

1. Pretraining AE, parameters such as the number of neurons, learning rate, and batch update numbers, must be continuously adjusted according to research needs during the training process. The number of three-layer neurons in AE was 128, 64, and 32, respectively. The optimizer was Adam, and the learning rate was 0.0001, 0.00025, and 0.0005, respectively. The epochs were set to 32, 48, and 64, and the batch size were 256, 256, and 128, respectively. The loss function adopted Logcosh, the evaluation index adopted R ², and the nonlinear activation function was not used in the three-layer AE.

2. To train the entire network, nine input layer and output layer nodes were selected, representing nine landslide factors. After each layer of AE training, the output layer and its parameters were discarded, and then connect the following layers. The number of network nodes in the last layer was one, representing the predicted value of unit landslide susceptibility. The activation function used Sigmoid, the optimizer used Adam, the loss function used balanced cross entropy (BCE), the evaluation index was AUC, the learning rate was 0.001, the epoch was 512, the batch size was 128, 80% of the valid data was used as the training set, and 20% of the valid data was used as the test set.

The model training process proposed in this study can be described as follows:

1. The landslide influence factors are preprocessed and are used as the initial input layer of the model.

2. The number of neurons in the hidden layer of an AE model is determined, and the gradient descent method is used to solve the weights and bias values from the input layer to the hidden layer and the hidden layer to the output layer. After training, the output layer and its corresponding weights and biases are removed, the encoding parameters of the input layer and the hidden layer are retained, and the hidden layer is used as the input layer of the next AE model.

3. Step (2) is repeated until the topmost AE model is trained.

4. the hidden layer of the topmost AE model and the label data is used as the input layer of the classifier, train the classifier, and fine-tune the parameters between the input layer and the hidden layer of each previous AE model so that the obtained classification results are optimal for completing the training of the entire model.

3.3 Evaluation index

1. Logcosh loss function

   Logcosh is a loss function for regression tasks that is smoother than mean squared error [37]. Logcosh is the logarithm of the hyperbolic cosine of the prediction error. The Logcosh equation is as follows:

   (9) Logcosh ( y , y p ) = ∑ i n log ( cosh ( y i p − y i ) ) .

   In the formula, y _i represents the landslide sample, y i p prepresents the predicted value, n represents the amount of data in the test data set.

2. Coefficient of determination

   The coefficient of determination R ² is used to evaluate the fit of the regression model coefficients after linear regression of the model [38]. The calculation method is as follows:

   (10) y ̅ = ∑ t = 1 N y t N ,

   (11) R 2 = 1 − ∑ t = 1 N ( y t − y ˆ t ) 2 ∑ t = 1 N ( y t − y ̅ t ) 2 .

   In the formula, y _t is the landslide sample, y ˆ t is the predicted value, and N is the amount of data in the test data set. The numerator represents the sum of the squared differences between the true value and the predicted value, that is, the error produced by using the model prediction, the denominator represents the sum of the squared differences between the true value and the mean value, that is, the error produced by using y = y ̅ . The accuracy of the model is judged according to the R ² value, and its value range is [0, 1]. The R ² value is close to 0, indicating that the model-fitting effect is poor, and on the contrary, the model error is small. Generally speaking, the larger the R ² value, the better the model fitting effect.

3. BCE loss function

   The BCE loss function is a binary loss function that can effectively fit the distribution [39]. The loss function is defined as follows:

   (12) BCE ( x ) i = − ∑ i = 0 n ( y i log y ˆ i + ( 1 − y i ) log ( 1 − y ˆ i ) ) / 2 ,

   where y ˆ i is the probability that the model predicts that the sample is a positive example, y _i is the sample label, if the sample is a positive example, the value is 1; otherwise, the value is 0.

4. AUC curve

To verify and improve the predictive ability of the DBN model, the AUC of the receiver operating characteristic curve (ROC) is used as the evaluation index for the predictive ability of the binary model [21]. AUC also considers the classifier’s ability to classify positive and negative samples and can still make a reasonable evaluation of the ability of the classifier in the case of sample imbalance and poor sampling. According to the confusion matrix of the binary classification problem, we can divide the prediction results of the model into four categories, as follows:

• The predicted and actual results are both positive; this is called true positive.

• The predicted results are positive, and the actual results are negative; this is called a false positive.

• The predicted and actual results are both negative; this is called the true negative.

• The predicted results are negative, but the actual results are positive; this is called a false negative.

The horizontal axis of the ROC curve corresponds to the false positive rate (FPrate), and the calculation formula is [13], which refers to the proportion of samples with positive prediction results among all samples whose true category is negative. The vertical axis corresponds to the true positive rate (TPrate), and the calculation formula is [14], which refers to the proportion of samples with positive prediction results among all samples whose true category is positive.

when FPrate and TPrate are equal, the model cannot distinguish between positive and negative samples. Currently, the ROC curve is a straight line, and the AUC = 0.5. Therefore, in the model evaluation, the AUC of the classification results of the binary model should be greater than 0.5, and the closer to 1, the better.

4.1 Model capability evaluation

In this study, the data set was divided into training set and test set according to 8:2. Figure 8(a–c) show the change curves of the improved DBN network model in the training phase of layer-by-layer feature extraction for the evaluation metric Logcosh loss function value and the evaluation metric R ² value. Figure 8(d) shows the change curves of the improved model in the training phase for the BCE loss function value and the evaluation metric binary classification accuracy.

Figure 8

Evaluation index curves of the training set and verification set.

As shown in Figure 8(a–c), in the layer-by-layer training of feature extraction, the loss value of AE in each layer can reach near 0 at the 2nd epoch. The R ² was almost the highest at the 10th epoch, and the R ² on the training set was 0.98, 0.92, and 0.88, respectively. The R ² on the test set reached 0.97 and above, indicating that the proposed model has a feature encoding extraction phase. Figure 8(d) shows the curves of the classification ability of the model in the training set. In Figure 8(d), we can see that on the 256th epoch, the loss value of the training set and the loss value function BCE of the validation set reached 0.26 and 0.24, respectively. The accuracy of the training set and the validation set reached 0.90 and 0.92, respectively, indicating that the improved model performed better and that the results of the landslide susceptibility evaluation were reliable.

4.2 Comparative analysis of landslide susceptibility evaluation results

To further analyze the reliability of the improved DBN method, the results of landslide susceptibility evaluation achieved by this study method were compared with those achieved by the hierarchical analysis composite index, SVM, and MLP methods. In this study, we used the natural breakpoint method of the ArcGIS platform to classify landslide susceptibility evaluation results acquired by the above four methods into non-susceptibility, low susceptibility, moderate susceptibility, and high susceptibility areas [21]. The visualization results are shown in Figure 9.

Figure 9

Landslide susceptibility evaluation map of Zhouqu County: (a) SIM, (b) SVM, (c) MLP, and (d) improved DBN.

ArcGIS software was used to calculate the area proportions of the four types and the proportions of the four types of landslide points (Figure 10). In addition, the AUC of the ROC was calculated (Figure 11).

• The high-susceptibility areas obtained from the improved DBN method accounted for 9.2% of the total study area; landslide points accounted for 76.6% of the total hazard points; the less susceptible areas were 78.4% of the total area; and landslide points accounted for 1.1% of the total landslide points, with an obvious district delineation and an AUC value of 0.90.

• The areas of high susceptibility obtained based on the SIM method were 22.1% of the total area; landslide points accounted for 45.7% of the total landslide points; the less susceptible area accounted for 14.0% of the total area; and landslide points accounted for 2.7% of the total landslide points. The distinguishing effect was second, and its AUC value was 0.66.

• The high-susceptibility area based on the SVM method accounted for 14.3% of the total area; landslide points accounted for 74.5% of the total hazard points, the less susceptible area accounted for 70.1% of the total area, and landslide points accounted for 5.3% of the total landslide points, with an obvious distinction, and its AUC value was 0.85.

• The high-susceptibility area obtained based on the MLP method accounted for 13.2% of the total area; landslide points accounted for 70.8% of the total hazard points; the less susceptible area accounted for 64.9% of the total area; landslide points accounted for 3.7% of the total landslide points, with an obvious distinction, and its AUC value was 0.86.

Figure 10

Ratios of landslide susceptibility assessment to historical landslides in Zhouqu County.

Figure 11

AUC curves of the four methods.

Among the landslide susceptibility evaluation maps of the four methods, the improved DBN, SVM, and MLP methods had the highest percentage of non-susceptibility areas and the lowest percentage of moderate susceptibility areas. Meanwhile, the landslide susceptibility evaluation results from the improved DBN model showed the lowest percentage of high-susceptibility areas and the highest percentage of non-susceptibility areas. However, the percentage of various SIM predictions was relatively balanced, and the prediction accuracy was generally low. Generally, the category of historical landslides is basically in the high susceptibility area (Figure 10), and previous studies have shown that the high susceptibility area should cover a small area in the landslide susceptibility evaluation [40]. It was pointed out that, in landslide susceptibility evaluation, the frequency of landslides in high-susceptibility areas should be the highest, while the frequency of landslides in low-susceptibility areas should be the lowest [41], which indicates that the improved DBN model susceptibility evaluation in this study has the best results and reliable method. The northeastern part of Zhouqu County has the highest probability of sliding, and areas assessed using all four methods reflected very high susceptibility. However, the SVM and MLP models have a wide range of prediction results for the high susceptibility area, which indicates that the SVM and MLP models cannot avoid the overfitting problem [21]. However, the DBN model had the best prediction results, avoided the model over-fitting problem, and had a higher prediction accuracy. Historical landslides are mainly concentrated in high susceptibility areas, while no historical landslides exist in non-susceptibility areas, indicating that the proposed model correctly predicts most landslides in highly susceptible areas, indicating that the improved DBN method based on the AE model can accurately predict the number of landslides in areas with high susceptibility. In addition, the improved DBN model was compared with the hierarchical analysis composite index, MLP, and SVM. Figure 11 shows the calculated ROC curves, and it can be seen that the proposed model had the highest AUC value of 0.97, while the AUC values of SIM, MLP, and SVM were lower than those of the proposed method. Of the four methods, the improved new method has the best prediction effect on landslide susceptibility areas.