**Figure 2.** The process of evaluating spatial agreement among various machine learning technique-based landslide susceptibility maps and optimizing the landslide prediction map. LSM=landslide susceptibility maps; KNN = K-Nearest Neighbor; MLP = Multi-Layer Perceptron; RF = Random Forest; SVM = Support Vector Machine.

The area receives the annual mean precipitation of 4288 mm [31]. The heavy rainfall triggers both flash floods and landslides in this area [11,12]. The majority of the historical landslides in Bangladesh occurred in this region [11]. For instance, a major landslide triggered by heavy rain in June 2017 killed at least 156 people in the south eastern hilly region of Bangladesh where the study area is located [32]. Unplanned urbanization, rapid growth of population, hill cutting, and deforestation are associated with the recent increase in landslide hazards [11,31]. Notably, Rohingya refugee camps, especially the Kutupalong camp of Ukhia sub-district (Figure 3) is located in areas that are highly susceptible to landslides. The Kutupalong camp is considered as the most densely populated refugee settlement area in the world, where around 75,000 people live per km2 [31,33]. Any catastrophic landslide will cause significant damage to human lives and assets. Hence, an accurate assessment of landslide susceptibility is paramount for developing a plan for landslide risk management.

#### *2.2. Landslide Inventory Mapping*

Landslide inventory mapping is one of the essential steps for landslide prediction and susceptibility mapping. This study utilized the landslide inventory map developed by Ahmed, Rahman, Sammonds, Islam, and Uddin [11] (Figure 3). They developed the latest landslide inventory map of the Cox's Bazar district by retrieving the historical landslide information from newspapers and various organizations and later verified those with global positioning system (GPS) and reconnaissance surveys. This study also used information about landslide movement type, its distribution and style, rate of flow, damage, the volume of displacement, material, and the reason for movement by preparing a landslide investigation form collected from Ahmed, Rahman, Sammonds, Islam, and Uddin [11]. A total of 1262 sample locations were used, where the number of landslide and non-landslide locations was 670 and 592, respectively. To develop the models, it is necessary to obtain non-landslide cells (where landslides did not occur). From the existing literature, Huang et al. [34] identified three methods for obtaining non-landslide grid cells: (i) the seed cell procedure; (ii) randomly selecting non-landslide locations from the landslide free areas; and (iii) non-landslide locations selected in areas with a slope lower than 2◦ . This study followed the second approach to select random locations within the study area,

where landslides did not occur. These cells provided the models with the necessary data during the training stage [35,36]. The sample locations were split into two classes: (i) 60% locations (52% landslide and 48% non-landslide locations) were used to train the machine learning-based landslide prediction models, and (ii) 40% testing locations (54% landslide and 46% non-landslide locations) were employed to evaluate the performance of the machine learning models (Figure 2).

**Figure 3.** Landslide inventory map of this study. Data source: [11].

#### *2.3. Landslide Conditioning Factor*

The performance of LSMs depends on the choice of landslide conditioning factors. Numerous studies on LSM have been conducted based on machine learning techniques [1,16,18,23,26,37,38], with various combinations of landslide conditioning factors being used. However, the selection of factors should be (i) based on their degree of affinity with landslide locations, (ii) measurable, (iii) non-redundant, and (iv) based on the knowledge of geomorphological characteristics of the area under study [2]. Based on the knowledge obtained from the literature, as well as, expert knowledge on the study area, a total of 12 variables were selected in this current study (Table 1). Areas with an elevation of less than 5 m, as well as, waterbodies and sandy sea beach areas (waterbody and restricted in Figure 4) were excluded from the LSMs [39].

**Figure 4.** *Cont*.

**Figure 4.** Landslide conditioning factors. SPI = Stream Power Index. NDVI = Normalized Difference Vegetation Index.

Topographical and hydrological parameters including aspect, elevation, slope, curvature, and Stream Power Index (SPI) are important factors that limit the density and spatial extent of landslides [2,37,38,40]. Raster maps of aspect, elevation, curvature, slope, and SPI were derived at 30-m spatial resolution from the Advanced Land Observing Satellite (ALOS) Digital Elevation Model (DEM) [28] (Figure 4a–e). Elevation influences landslides primarily by affecting different biophysical parameters and anthropogenic activities. Only a limited number of studies, conducted on a specific basin, found that landslides occur at certain elevations [41]. Elevation can determine the spatial variability of landslides because it is affected by geological tectonics [37]. It can also influence the occurrence of landslides by impacting other causative factors such as slope, curvature, and SPI [42]. Aspect, indicating the direction of slope [43], indirectly influences the distribution of landslide locations by affecting the general physiographic trend of the area and/or the main precipitation direction [37,40]. Slope angle is considered as one of the most influential factors for the occurrence of landslides, as it affects the concentration of moisture and the level of pore pressure, as well as, controls regional hydraulic continuity [2,40]. All of these processes influence slope instability [8]. Curvature is also considered as a landslide influencing factor that directly controls the velocity of water flow, delimiting erosion [8,40]. SPI also determines the erosion potential of the surface [43] and is considered as an essential predictor of landslides [37,38]. Areas with high SPI values indicate a higher erosion potential, while negative values suggest no predicted erosion [44,45]. In this study, a layer of SPI was derived using the following equations in GIS:

$$SPI = A\_s \times \tan \beta \tag{1}$$

where *As* and β indicates the specific catchment area (m2 /m) and slope gradient, respectively [43].

It is widely conceived that various geological factors significantly influence the occurrence of landslides, as these factors often lead to a difference in strength and permeability of rocks and soils [2]. This study considered the three geological factors of surface geology, soil type, and soil texture (Figure 4i–k). Digital geologic and geophysical data of Bangladesh were collected from the U.S. Geological Survey [46]. The surface geology map of the study area includes a total of 11 classes: water (H2O), Bhuban formation (Miocene, Tb), Dupi Tila formations undivided (QTdd), valley alluvium and colluvium (ava), Girujan clay (Pleistocene and Neogene, QTg), Tipam Sandstone (Neogene, Tt), Boka Bil formation (Neogene, Tbb), beach and dune sand (csd), marsh clay and peat (ppc), Dupi Tila formation (Pleistocene and Pliocene, QTdt), and Dihing formation (Pleistocene and Pliocene, QTdi) (Figure 4i). Primary-level parameters such as soil type and soil texture are essential predictors of landslides. These parameters determine the amount of moisture content indicating the degree of stability of the soil [25,37,47,48]. Soil type and soil texture data were collected from the Bangladesh Agricultural Research Council [49].

Other anthropogenic, environmental, and locational factors considered in this study include distance to stream, land cover, normalized difference vegetation index (NDVI), and distance to road (Figure 4f–h,l). The land cover and NDVI maps of the year 2020 were prepared using Landsat satellite images based on the Google Earth Engine Platform. The land cover map was developed applying a supervised classification technique with the Random Forest algorithm. In the case of southern Bangladesh, a recent study demonstrated that this method has a higher classification accuracy compared to other land cover classification techniques [50]. The land cover map contains five classes: bare land, built-up area, crop land, vegetation, and waterbody (Figure 4g). Proximity to roads explains the locations of landslides, as the artificial and natural slopes adjacent to a road are sensitive to this hazard [51]. Road-cuts, excavation, and additional load can induce anthropogenic instability of the soil, promoting landslides [2,5]. A layer of distance to road network was developed using the Euclidian distance algorithm. Likewise, the location of areas with respect to natural drainage channels can also demonstrate the locations of landslides [11], as streams may change the stability of an area by eroding the slopes [5,51]. In this study, distance to stream networks was derived from the ALOS DEM. Again, by applying the Euclidian distance algorithm, a map of distance to stream was generated.

**Table 1.**Landslide conditioning factors used in this study.

#### *2.4. Multi-Collinearity Analysis of Landslide Conditioning Factors*

The selected landslide causative factors could be subject to multi-collinearity; hence, it is necessary to estimate the correlation of independent variables before modelling landslide susceptibility [8]. To eliminate the factors susceptible to multi-collinearity, this study determined variance inflation factors (VIF) [53] of 12 selected landslide conditioning factors using R [54]. VIF is a well-known method to determine the multi-collinearity of landslide conditioning factors [8,55]. A VIF value of a variable exceeding 5 indicates potential serious multicollinearity [53,55]. In this study, the selected landslide conditioning factors yielded VIF values < 2.8, indicating the absence of potential multi-collinearity (Table 1).

#### *2.5. Landslide Susceptibility Modelling*

#### 2.5.1. Pre-Processing

Using the binary locations (landslide and non-landslide), values of the selected 12 conditioning factors were extracted in a geographic information system (GIS) environment. As evident in Table 1, eight were continuous variables, while the remaining four variables had discrete characteristics. In order to represent discrete (categorical) variables semantically, they must be considered as a composite feature (where the number of generated binary features and the number of categories are equal). These discrete variables were encoded using a one-hot encoding scheme [31], implying that multiple binary features were generated to represent a single discrete feature. The number of one-hot encoded features depends on the number of variable classes. For instance, there are 11 categories in the geology variable. If a landslide location was found in a geology class, a value 1 was encoded to the class, while the other 10 classes were encoded as 0. This data pre-processing method was applied for all other discrete variables. For each variable, mean and standard deviation were calculated. The mean of each variable was then subtracted from the corresponding value in a variable and divided by the standard deviation. This reduces training time since optimization routines have a smaller parameter space to traverse.

#### 2.5.2. Hyperparameter Optimization

Hyperparameter optimization can improve the accuracy of machine learning algorithm-based models. The process aims to select the optimal hyperparameter values according to the evaluation index [56]. Three approaches are frequently used for optimizing hyperparameters: grid search, random search, and Bayesian optimization [57]. This current study applied the grid search technique along with 5-fold cross-validation on the training set to perform hyperparameter optimization. Hyperparameters that provide the best performance were chosen for final training and testing samples of respective machine learning models. For instance, the optimal number of neighbors of five in the KNN (Table 2) indicates that values of landslide conditioning factors corresponding to a landslide location were compared against the values of landslide predictors of five other sample locations, to obtain the most reliable prediction. Table 2 summarizes the hyperparameters, their search range and optimal values of the four models.

**Table 2.**Hyperparameters, search range, and optimal values of the machine learning-based landslide susceptibility models.

#### 2.5.3. Machine Learning Models

(1) K-Nearest Neighbor (KNN)

The KNN algorithm classifies an instance (landslide or non-landslide) that is mostly represented within its (*k*) neighbors. The parameter *k* is often a small positive integer [58]. The proximity between the samples is measured using a distance metric. The distance metric indicates how similar or different are the profiles of conditioning factors for any given two samples. Data points with similar conditioning factors will have a small feature distance between them. Though the model is simple in terms of hyperparameters, it becomes computationally expensive as the number of samples becomes large. The landslide susceptibility associated with a certain set of values of conditioning factors is determined by calculating its distance to each training data point (in high-dimensional feature space). The *k* nearest data points are used to determine the landslide susceptibility. The dominant susceptibility class within those *k* nearest neighbors (i.e., the class with the highest number of members in the *k* members) becomes the class membership of the new data point [14,59].

#### (2) Multi-Layer Perceptron (MLP)

Multi-Layer Perceptron (MLP) is a type of neural network with one or more hidden layers. Due to the presence of a hidden layer, the internal representations (as higher-order intermediate features) can be learned. Each layer consists of one or more neurons; the outputs (activations *a*) can be represented by Equation 2 [31]. The output of the *i*th neuron in the *j*th layer was obtained by calculating the sum of activations from the previous layer (*i* − 1) weighted by the parameters of layer *i* and then passing into an activation function *f*. Considering that there are several types of activation functions (sigmoid, hyperbolic tangent, rectified linear unit), the choice of activation functions is discussed in Section 2.5.3 (Hyperparameter Optimization). In this study, since a total of 23 features were derived by one-hot encoding during the pre-processing step, the first input layer of the MLP had 23 neurons. The resultant map was represented in terms of the probability of landslide occurrence.

$$a\_i^j = f(\sum\_{k=0}^n \omega\_k^j a\_k^{j-1}) \tag{2}$$

where *f* is the activation function, ω *j k* is the weight of *k*th neuron in layer *j*, *a j*−1 *k* is the activation of neuron k in layer *j* − 1 (the previous layer), *j* is the layer index, *i* is the neuron index, and *n* is the number of neurons in layer *j*.

#### (3) Random Forest (RF)

Random forest (RF) is considered as a powerful ensemble-learning method that can be applied for classification, regression, and unsupervised learning [18]. This method has been widely applied in landslide susceptibility mapping [18,56,60]. Ensemble models generally train several weak learners and then take their aggregated outputs to obtain more reliable predictions. The RF algorithm builds weak learners in the form of decision trees. It estimates the mean of outputs of the individual weak learners, as shown in Equation (3). Each weak learner (*b*) corresponds to a function *fb* (*x*). The RF uses bootstrap aggregating where the weak learners train parallelly [31].

$$\mathcal{F}(\mathbf{x}) = \frac{1}{B} \sum\_{b=1}^{B} f\_b(\mathbf{x}) \tag{3}$$

where *F*ˆ(*x*) is the ensembled prediction from weak learners, *B* is the total number of weak learners, *b* is the weak learner index, and *fb* (*x*) is the function for *b*th weak learner.

#### (4) Support Vector Machine (SVM)

Support Vector Machine (SVM) is also a widely used machine learning algorithm in landslide susceptibility mapping [8,16,23,26]. This supervised learning method separates the classes with a decision surface that maximizes the margin of class boundaries [8]. The training locations, closest to the optimal hyperplane, are called support vectors [23]. Suppose each sample location has M number of features, the objective of the SVM algorithm is to find a hyperplane in M dimensional feature space that separates the samples of different classes. A hyperplane is R(M−1) dimensional in RM. A hyperplane in R 2 is a line, a hyperplane in R3 is a plane, and so on. This hyperplane functions as a decision boundary, which determines the label of a sample (i.e. landslide or non-landslide). The margin around the hyper lane indicates that value exceeding 1 denotes a positive sample (landslide), and a value equal to −1 denotes a negative sample (non-landslide).  , *Xn*) is the vector of landslide affecting factor and *Yj* (*Y*1, *Y*2) is the vector of landslide (1) or non-landslide (0) event, the optimal hyperplane can be found by solving Equation (4) [26].

$$f(\mathbf{x}) = \text{sign}\left[\sum\_{i=1}^{n} \alpha\_i Y\_j k(X\_i X\_i) + k\right] \tag{4}$$

where *k* is the offset from the origin of the hyperplane, *n* is the total number of factors that affects landslide, α*i* is the positive real constant, and *k*(*X*, *Xi* ) is the Kernel function. To classify the binary events (landslide or non-landslide), the condition to solve Equation (4) was assumed as below:

$$\mathcal{Y}\_j[\boldsymbol{\omega}^T \boldsymbol{\varphi}(\mathbf{x}\_i) + \mathbf{c}] \ge 1 \iff \begin{cases} \boldsymbol{\omega}^T \boldsymbol{\varphi}(\mathbf{x}\_i) + \mathbf{c} \ge 1 \text{ if landside events occur} \left(\boldsymbol{Y}\_j = 1\right) \\\ \boldsymbol{\omega}^T \boldsymbol{\varphi}(\mathbf{x}\_i) + \mathbf{c} \le 0 \text{ if landside events not occur} \left(\boldsymbol{Y}\_j = 0\right) \end{cases} \tag{5}$$

where *w* is the weight vector and ϕ(*xi* ) is the total number of factors that affects landslide.

#### 2.5.4. Performance Evaluation Methods

The performance of landslide susceptibility models was evaluated using a well-known method called receiver operating characteristic (ROC) curve and subsequent area under the curve (AUC) [1,11,23,31,43]. The ROC curves were developed using the 40% sample testing data. The ROC curve indicates the performance of a binary classifier system, representing sensitivity as a function of the false positive rate (1-specificity).

The sensitivity of a model is the ratio of the number of true positives to the sum of the number of true positives and false negatives. The specificity is the ratio of the number of true negatives to the sum of the number of true negatives and false positives. The ROC curve can be developed by plotting sensitivity in the y-axis against the cumulative distribution function of the false positive rate in the x-axis. The estimated AUC value can be categorized as poor (0.5–0.6), average (0.6–0.7), good (0.7–0.8), very good (0.8–0.9), and excellent (0.9–1) [1,43,60]. Besides, various statistical indices such as overall accuracy, precision, recall, and F1-score were estimated by developing a confusion matrix [1,11,43].

#### *2.6. Evaluation of Spatial Agreement and Optimizing Prediction Map*

To evaluate the inter-model agreeability, a pixel-wise agreement between two machine learning algorithms was estimated. Therefore, Pearson's correlation coefficient was estimated for a total of six possible combinations of machine learning model-based landslide susceptibility maps. Here, Pearson's correlation coefficient indicates the covariance of landslide predictions, obtained by using two algorithms, divided by the product of their standard deviations. The correlation coefficient can range from +1 to −1, where values zero as indicating no agreement and ±0.29 as low degree, ±0.30–±0.49 as moderate degree, ±0.50 to <±1 as high degree, and ±1 as perfect agreement [60].

Following the evaluation of spatial agreement, an optimized landslide prediction map was developed combining susceptibility maps generated by applying the four machine learning algorithms. The combined map was developed by following a methodology proposed by Rossi, Guzzetti, Reichenbach, Mondini, and Peruccacci [22], where they established a logistic regression (LR) model. The LR model included binary landslide and non-landslide locations as the dependent variable and the results of the four landslide susceptibility models as the independent variables. The obtained regression coefficients were incorporated in Equation (6) [43] in GIS to derive the probability (*P*) of landslides in the study area.

$$P = \frac{1}{1 + e^{-z}} \tag{6}$$

where *z* is the linear combination of independent variables which was estimated using the following equation:

$$z = \theta\_0 + \theta\_1 \mathbf{x}\_1 + \theta\_2 \mathbf{x}\_2 + \dots + \theta\_n \mathbf{x}\_n \tag{7}$$

where θ0 is the intercept of the model, θ*i* (*i* = 1, 2, . . . , *n*) indicates the regression coefficient of independent variables, and *xi* (*i* = 1, 2, . . . , *n*) represents the *n* number of independent variables. Validation of the resultant combined model was performed by developing the ROC curve by using the 40% testing data.

#### **3. Results**

#### *3.1. Landslide Susceptibility Modelling*

#### 3.1.1. Landslide Prediction

Figure 5 shows landslide susceptibility maps of Cox's Bazar district developed by applying the four machine learning algorithms. The generated landslide probability maps were classified into five categories each by applying the Jenks natural breaks classification method in GIS: (i) very low (0–0.1), (ii) low (0.11–0.3), (iii) medium (0.31–0.5), (iv) high (0.51–0.85), and (v) very high (0.86–1). As evident in Figure 6, the proportion of landslide susceptible area varied from one model to another. Among all methods, the SVM resulted in the highest proportion of area (38.7%) susceptible to the landslide of 'high' and 'very high' severity, while the Random Forest (RF) algorithm yielded a relatively lower proportion (23.1%) of landslide susceptible area. Likewise, the ratio of the population exposed to 'high' and 'very high' landslide susceptible zones varied for different algorithms. For all the four methods, the percentage of landslide exposed population ranged between 34% to 48% (Figure 6).

**Figure 5.** *Cont*.

**Figure 5.** Landslide susceptibility maps obtained by four machine learning algorithms: (**a**) K-Nearest Neighbor (KNN), (**b**) Multi-Layer Perceptron (MLP), (**c**) Random Forest (RF), and (**d**) Support Vector Machine (SVM).

**Figure 6.** Landslide susceptible area obtained using five models: K-Nearest Neighbor (KNN), Multi-Layer Perceptron (MLP), Random Forest (RF), Support Vector Machine (SVM), and combined model. Blue dots indicate the proportion of people exposed to 'high' and 'very high' susceptible zones.

#### 3.1.2. Evaluation of Models' Performance

To evaluate the performance of various landslide susceptibility models, a performance matrix was derived using the test samples (40% of the total data) (Table 3). The performance evaluation indices indicated a very high prediction accuracy of all the models. In terms of overall accuracy, the RF classifier resulted in the highest accuracy (96.63%), followed by the MLP (95.45%), SVM (94.06%), and KNN (90.69%). However, the overall accuracy is a universal metric, hence, it does not indicate which specific classes were being inaccurately classified. To obtain further insights into the agreement

between the observed and modelled locations (landslide and non-landslide)—precision, F1-score, and recall were estimated (Table 3). The RF classifier achieved the best accuracy with respect to all performance indicators. The MLP followed closely and consistently in terms of all indicators. In relation to the estimated AUC values, the RF classifier yielded the highest accuracy (0.962), followed by the MLP (0.960), SVM (0.935), and KNN (0.927) (Figure 7). The relatively greater values of performance indicators of RF and MLP can be attributed to their ability to learn complex relationships between geospatial characteristics of an area and the occurrence of landslides [17,18,60].

**Model Overall Accuracy Precision F1-score Recall Non-Landslide Landslide Non-Landslide Landslide Non-Landslide Landslide** KNN 0.9069 0.9227 0.9227 0.9015 0.9015 0.8811 0.8811 MLP 0.9545 0.9547 0.9547 0.9528 0.9528 0.9508 0.9508 RF 0.9663 0.9633 0.9633 0.9652 0.9652 0.9672 0.9672 SVM 0.9406 0.9385 0.9385 0.9385 0.9385 0.9385 0.9385

**Table 3.** Performance evaluation indicators of the machine learning based landslide susceptibility models.

**Figure 7.** Receiver operating characteristic (ROC) curves of the five models: K-Nearest Neighbor (KNN), Multi-Layer Perceptron (MLP), Random Forest (RF), Support Vector Machine (SVM), and combined model.

#### *3.2. Spatial Agreement of Various Methods*

This study developed a correlation matrix by comparing pixel-wise landslide probabilities between various methods (Figure 8) to evaluate the extent of agreement of one landslide susceptibility model over another. Although the values of AUC were very similar for various methods (Figure 7), a substantial difference in the agreement was observed in LSMs obtained using the different techniques. Overall, the correlation coefficient ranges from 0.69 to 0.85 (Figure 8). The combinations of SVM-RF resulted in the highest degree of the agreement, while the KNN-SVM yielded the lowest agreement.

**Figure 8.** Correlogram to show the agreement between the five landslide susceptibility models: K-Nearest Neighbor (KNN), Multi-Layer Perceptron (MLP), Random Forest (RF), Support Vector Machine (SVM), and combined model.

#### *3.3. Aggregated Landslide Susceptibility Mapping*

− − Since spatial heterogeneity in landslide prediction exists between the different machine learning-based approaches, an aggregated susceptibility map combining the outputs of all algorithms would minimize the uncertainty of individual methods. In this study, a regression-based approach was adopted. A multivariate logistic regression (LR) was established incorporating the binary landslide inventory data as the dependent variable and the results of the four landslide susceptibility models as independent variables. The outcome of the LR model is summarized in Table 4. Among the four models, the MLP, RF, and SVM were statistically significant (*p*-value < 0.05). The coefficient of determinants (R2 ) of 0.80 indicates a very good model performance. In relation to the estimated regression coefficients, the RF model had the highest degree of agreement with the landslide inventory, followed by the MLP, SVM, and KNN. The pattern of influence of various models in predicting landslides corresponds to their level of accuracy in terms of their respective AUC values (Figure 7).

−

The estimated LR coefficients of the four models and intercept were incorporated in Equation (6) to derive the combined LSM. Again, the resultant aggregated map was categorized into five classes applying the Jenks Natural Break algorithm (Figure 9). The combined susceptibility map yielded the highest AUC value (0.965) compared to the single susceptibility forecasts (Figure 7). About 26.8% of the total study area was within the 'high' and 'very high' landslide susceptible zones, where approximately 21.7% of total population inhabit (Figure 6). In respect to spatial agreement, the combined LSM resulted in greater spatial agreement with the all four models, with the correlation coefficient ranging between 0.85 and 0.92 (Figure 8).

**Figure 9.** (**a**) Combined landslide susceptibility map of Cox's Bazar district, (**b**) ratio of landslide susceptible zones (high and very high) in various sub-districts (Upazila), and (**c**) landslide susceptibility in the Rohingya refugee camps of Ukhia sub-district.

The extent of landslide susceptible areas varies in different sub-districts (*Upazila*) of Cox's Bazar. Teknaf Upazila is the most susceptible, where more than 8% of the total study area was susceptible to landslides of 'high and 'very high' severity (Figure 9b). A substantial proportion of area (7% of the study area) in Ukhia sub-district was also susceptible. The Rohingya refugee camps in this area were located within high and very high landslide susceptible zones. Various recent studies also found that changes in the geomorphological, hydrological, and anthropogenic environments due to the Rohingya influx caused their settlement areas vulnerable to landslides [11,31].

#### **4. Discussion**

The spatial disagreement in prediction among various techniques creates challenges in selecting the most suitable susceptibility map for managing landslide hazards [9,21,22]. The current study seeks to address this challenge by estimating the extent of spatial agreement, as well as, proposing a method to combine landslide susceptibility maps and thus incorporating the valid results of various models. The study focused on the Cox's Bazar district of Bangladesh, which is well known as being vulnerable to landslide disasters [11,12,31]. First, LSMs were developed by applying four machine learning algorithms—K-Nearest Neighbor (KNN), Multi-Layer Perceptron (MLP), Random Forest (RF), and Support Vector Machine (SVM)—featuring hyperparameter optimization. In comparison to the existing studies on LSM of Cox's Bazar district [11,12,31], this current study employed up-to-date data of landslide conditioning factors. In addition, it restricted low-lying areas (waterbodies and elevation < 5 m) in susceptibility mapping, otherwise, the resultant maps would have been prone to overestimation of landslide susceptible zones, as was the case of some recent studies. While evaluating the models' performance, all of them yielded very high prediction accuracy, with the AUC values ranging between 0.927 to 0.962. The results of various recent studies have also ascertained that different machine learning-based models yielded high accuracy in predicting landslides [1,15,17,26,27,61].

This study hypothesized that different susceptibility models can result in different LSMs, despite incorporating similar landslide inventory data. The assessment of spatial agreement between various models revealed spatial heterogeneity in landslide predictions, with the estimated pixel-wise correlation coefficients of landslide probability between various models ranging from 0.69 to 0.85. The spatial distribution of landslide susceptibility obtained in this study was also different than that of a recent study conducted in Cox's Bazar district of Bangladesh [11]. This highlights the uncertainty in landslide predictions of various models, despite their considerable prediction accuracy in terms of the AUC values. Most of the existing studies on machine learning-based LSM had a major focus on identifying the most suitable method for predicting this natural phenomenon [1,8,18,23,60], while little attention has been given in analyzing uncertainties resulting from the spatial disagreement in landslide prediction [9]. The current study is the first case study-based contribution to investigate this major gap in the existing literature.

This study further developed a combined LSM by integrating the results of the four machine learning-based landslide predictions, adopting a method proposed by Rossi, Guzzetti, Reichenbach, Mondini, and Peruccacci [22]. The result indicates an improvement in landslide prediction accuracy. Existing studies, which applied multiple machine learning algorithms to map landslide susceptibility, mainly evaluated different methods based on quantitative measures [8,16,18]. Whilst quantitative measures of model fit are useful, they are not conclusive in determining the efficacy and reliability of susceptibility assessment [22]. A combined landslide susceptibility map that this study developed would help to minimize the uncertainties of individual methods.

#### **5. Conclusions**

Predicting a complex natural phenomenon such as a landslide is a challenging task and subject to considerable uncertainties. An accurate prediction of landslides is the prerequisite for managing this hazard. In this study, the spatial association in landslide prediction between various machine learning-based models was analyzed to quantify the spatial agreement of predicted landslide susceptibility. By addressing uncertainties in various models, this study also developed a landslide susceptibility map combining the outcomes from various models. The results indicate an improvement in landslide prediction compared to the individual models.

Despite achieving an improved result in landslide prediction, this study has some limitations that could be addressed in future results. Landslide inventory data used in this study was developed based on various secondary sources and validated through fieldwork [11]. Scarcity of data, including detailed landslide inventory on the study area, made it difficult to model a landslide more accurately. Accuracy of the LSM results depended on input parameters used, particularly the DEM. The ALOS

DEM of 30-m resolution used in this study had a low root-mean-square error (1.78 m) in vertical accuracy and was considered to be the most accurate freely available DEM [43,62]. However, for future research, high-resolution DEM could be employed to improve the existing landslide susceptibility modelling frameworks.

This study is an attempt to integrate results of multiple machine learning-based landslide susceptibility models to minimize uncertainties and improve landslide predictions. The modelling framework used in this study could be transferred to other landslide-susceptible regions. Landslide susceptibility maps can enable urban planners in identifying suitable areas for urban development [63]. The combined landslide susceptibility map of Cox's Bazar district could be useful to policymakers and practitioners in sequencing and prioritizing interventions in managing landslides. The proposed model is an advancement in the existing landslide susceptibility models that intends to predict landslides more accurately. The results of this model could be utilized in improving the existing landslide early warning system [11], to strengthen landslide disaster risk mitigation strategies to support for the resilient future of inhabitants of the study area.

**Author Contributions:** Conceptualization, M.S.G.A. and B.A.; methodology, M.S.G.A., M.S.R., and N.A.; software, M.S.G.A., M.S.R., N.A., and R.M.R.; validation, M.S.G.A., M.S.R., and N.A.; formal analysis, M.S.G.A., M.S.R., and N.A.; investigation, M.S.G.A., B.A., R.M.R., and M.F.R.; resources, M.S.G.A., M.F.R., and B.A.; data curation, M.S.G.A. and B.A.; writing—original draft preparation, M.S.G.A.; writing—review and editing, M.S.G.A. and B.A.; visualization, M.S.G.A. and B.A.; supervision, B.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

# *Article* **Evaluating the E**ff**ects of Digital Elevation Models in Landslide Susceptibility Mapping in Rangamati District, Bangladesh**

#### **Yasin Wahid Rabby 1, \* , Asif Ishtiaque 2 and Md. Shahinoor Rahman 3**

Received: 16 July 2020; Accepted: 20 August 2020; Published: 22 August 2020

**Abstract:** Digital elevation models (DEMs) are the most obvious data sources in landslide susceptibility assessment. Many landslide casual factors are often generated from DEMs. Most studies on landslide susceptibility assessments rely on freely available DEMs. However, very little is known about the performance of different DEMs with varying spatial resolutions on the accurate assessment of landslide susceptibility. This study compared the performance of four different DEMs including 30 m Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model (GDEM), 30–90 m Shuttle Radar Topographic Mission (SRTM), 12.5 m Advanced Land Observation Satellite (ALOS) Phased Array Type L band Synthetic Aperture Radar (PALSAR), and 25 m Survey of Bangladesh (SOB) DEM in landslide susceptibility assessment in the Rangamati district in Bangladesh. This study used three different landslide susceptibility assessment techniques: modified frequency ratio (bivariate model), logistic regression (multivariate model), and random forest (machine-learning model). This study explored two scenarios of landslide susceptibility assessment: using only DEM-derived causal factors and using both DEM-derived factors as well as other common factors. The success and prediction rate curves indicate that the SRTM DEM provides the highest accuracies for the bivariate model in both scenarios. Results also reveal that the ALOS PALSAR DEM shows the best performance in landslide susceptibility mapping using the logistics regression and the random forest models. A relatively finer resolution DEM, the SOB DEM, shows the lowest accuracies compared to other DEMs for all models and scenarios. It can also be noted that the performance of all DEMs except the SOB DEM is close (72%–84%) considering the success and prediction accuracies. Therefore, anyone of the three global DEMs: ASTER, SRTM, and ALOS PALSAR can be used for landslide susceptibility mapping in the study area.

**Keywords:** landslide susceptibility; Bangladesh; digital elevation model; random forest; modified frequency ratio; logistic regression

#### **1. Introduction**

Local terrain conditions, including terrain relief, hydrology, geology, and land use, are crucial for the assessment of landslide susceptibility [1]. These local features are often termed as "causal factors" [2]. Identifying these causal factors is considered as the steppingstone of landslide susceptibility assessment. Landslide susceptibility represents the likelihood of landslide occurrence in an area. It assumes that future landslides may occur in previous landslide locations where the causal factors already created conducive environments for triggering landslides [2–4]. Several natural and anthropogenic

factors can trigger landslides, thus called "triggering factors," such as volcanic activity, groundwater excavation, prolonged rainfall, rapid snow melting, hill cutting, deforestation, land-use change, etc. [5–8]. Moreover, landslide susceptibility assessment relies on the characteristics of landslide inventory—a detailed register of distribution and characteristics of past landslides [9,10]. Therefore, the success of landslide susceptibility assessment largely depends on the selection of causal factors and the quality of landslide inventories.

Landslide susceptibility can be assessed qualitatively or quantitatively [11]. Qualitative approaches of landslide susceptibility assessment are based on experts' judgments on causal factors. However, the mathematical relationships between landslide locations and casual factors are utilized in quantitative approaches [6]. Frequently, mixed methods and semi-quantitative methods are adopted to process experts' opinions in qualitative assessments. The widely-used methods in this domain are the analytical hierarchy process (AHP) [12], fuzzy logic [13], and GIS-based AHP [14]. In bivariate techniques, each causal factor is divided into a set of classes, and landslide locations are compared with each class. Thus, the bivariate relationship is established between landslide occurrence and one factor-class at a time [4]. The commonly used bivariate methods are frequency ratio, the weight of evidence, fuzzy logic, evidential belief function, and statistical index [3,15]. In contrast, the multivariate statistical methods determine the relationship between landslide occurrence and multiple causal factors. Examples of multivariate methods are logistic regression, adaptive regression spline, general additive models, and simple decision trees [16].

One of the limitations of bivariate and multivariate models is that they are constrained by normality and collinearity assumptions. Compared to these models, machine learning-based models are relatively less limited by these assumptions [17] and, therefore, can consider the nonlinear nature of landslides [18]. Some argue that machine learning-based models such as the random forest, gradient boosting, and support vector machines often outperform both bivariate and multivariate statistical models [19,20]. While the selection of methods is essential, landslide susceptibility assessment also depends on the types and quality of causal factors, mapping unit, and the scale of investigation [21]. Many causal factors are often derived from analyzing satellite imageries and topographic models, including land cover, elevation, slope, aspect, and hill cut [7]. Thus, these derived causal factors are often impacted by the spatial resolution of sources, geometric error, and instrument or sensor type. In landslide susceptibility mapping, digital elevation models (DEMs) often replace topographic maps to derive the most important causal factors (e.g., slope, topography, aspect). Moreover, many developing countries may not have topographic maps. Thus, landslide studies from these countries usually rely on the free of charge DEMs derived from remote sensors.

DEM is the digital representation of the earth's surface. It is widely used in various research areas in which topography plays an important role, such as hydrological modeling, geomorphological analysis, and feature extraction, landslide susceptibility and hazard assessment, erosion susceptibility, and glacier monitoring [22]. DEMs are often generated using data obtained from different remote sensors, including optical imaging sensors, light detection, and ranging (LiDAR), and synthetic aperture radar (SAR) [23]. The qualities of DEM-derived factors often depend on the spatial resolution of DEMs. Therefore, the choice of DEM is important for the assessment of landslide susceptibility [1]. To this day, a few attempts have been made to compare the performance of DEMs with different spatial resolutions in landslide susceptibility assessment. For instance, Dietrich et al. [24] compared different DEMs and found a similarity in performance irrespective of spatial resolution in identifying moderate landslide class (see also [25]). They argued that the resolution of DEM might not be very important to represent the slope failures. Similarly, Tian et al. [26] contended that finer resolution does not essentially lead to higher accuracy in landslide susceptibility assessment (see also [1,27]).

These past studies further indicated that the performance of DEMs is context-dependent meaning that the performance of a DEM in a region may not be assumed to be similar in another region [28–30]. They also argued that DEMs with fine spatial resolution may not necessarily have better performance over coarse resolution DEMs. Therefore, it is an utmost need to have comparative assessments of DEMs

in various contexts. We found that even though some parts of Bangladesh are vulnerable to landslides, no study investigated the relative performance of different DEMs in landslide susceptibility assessments. Against this backdrop, this study contextualized landslide susceptibility in Bangladesh and compared the performance of different DEMs and modeling techniques. The study area is selected from Bangladesh because the hilly southeastern parts of the country encounter landslides almost every year that often claim tens of lives [31–34]. Because of unavailability of LiDAR, the majority of the landslide susceptibility-related studies in Bangladesh used 30m Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model (GDEM), 30–90 m Shuttle Radar Topographic Mission (SRTM), and 12.5 m Advanced Land Observation Satellite (ALOS) Phased Array type L-band Synthetic Aperture Radar (PALSAR) to obtain various topographic factors [7,31,32,35]. Survey of Bangladesh (SOB) has developed a DEM (25 m) for the whole of Bangladesh and no study has used this DEM in landslide susceptibility assessment. Therefore, the landslide susceptibility maps in Bangladesh are influenced by the usage of different DEMs and which DEM provides the most accurate susceptibility maps are largely unexplored [36]. As many causal factors are derived from DEM datasets, the selection of appropriate DEM is crucial for landslide susceptibility assessment. This study aims to compare the relative performance of four DEMs: ASTER GDEM (30 m), ALOS PALSAR (12.5 m), SRTM (30 m), and Survey of Bangladesh (SOB) DEM (25 m). This study further aims to compare three quantitative landslide susceptibility assessment techniques: modified frequency ratio (bivariate method), logistic regression (multivariate method), and random forest (machine learning method) [4,8,29,37].

#### **2. Study Area**

Rangamati, a hilly southeastern district of Bangladesh is selected as the study area because of the regular landslide occurrences in this area. More than 100 people died, and 12,000 families suffered losses due to landslides in this district [38]. Most landslides in Rangamati take place in three upazilas (sub-districts): Rangamati Sadar, Kaptai, and Kawkhali [35]. As such, the study area is narrowed down to these three Upazilas (Figure 1). The combined geographical area of these three sub-districts is 1145 km2 and more than 40% of it is forested [39]. The geology of this area comprises Dhihing, Dupi tila, Girujan clay, Bhuban, Bokabil, and Tipam sandstone (Figure A3f of Appendix C)). The bedrock and soil structure of the areas are not stable, which makes the hills highly prone to landslides [40]. Climatologically, this area falls under a tropical monsoon climate, and the annual average temperature varies from a maximum of 36.5 degrees to a minimum of 12.5 degrees Celsius, and annual rainfall is 2673 mm with mean humidity level 71.6% [39].

**Figure 1.** Study area and landslide and non-landslide (pseudo-absence points) locations used ‐ ‐ in modeling.

#### **3. Synopsis on Data Utilization**

#### *3.1. Landslide Inventory*

Landslide inventory records different information, including the exact location, size, type, time of occurrence, causalities, trigger, and causes [2]. Rabby and Li [33] prepared and published [34] landslide inventory of the Chittagong hilly area, Bangladesh. They used participatory field mapping to prepare this inventory [41]. In our study, we selected 168 landslides from their inventory as it has been advised to use more than one method in landslide inventory preparation [42], we analyzed available Google Earth images on the Google Earth platform to map more landslides in the study area. We used the method proposed by Rabby and Li [33] for Google Earth mapping and mapped 93 landslides that occurred from January 2001 to January 2019. Therefore, in our study, we used a total of 261 (168 + 93) landslide locations. The mean size of the landslide was 274.2 m2 . The smallest and the largest dimensions of the landslides were about 14.6 m2 and 3422.4 m2 , respectively.

#### *3.2. Landslide Causal Factors*

In our study, we used a total of 15 causal factors (Table 1) to produce landslide susceptibility maps (see Appendices A–C). Of these factors, seven factors were derived from DEM: elevation, slope, plan curvature, profile curvature, topographic wetness index (TWI), stream power index (SPI), and aspect. As we are comparing the performance of different DEMs, we derived each of these factors from four different DEMs: ASTER, SRTM, ALOS PALSAR, and SOB. The rest of the eight factors namely land use/land cover, land use/land cover change, geology, distance from the road networks, distance from the fault lines, distance from the drainage networks, rainfall, and normalized difference vegetation index (NDVI) were collected from different datasets. Maps of different causal factors had different resolutions, but for the convenience of comparison, we kept the 30 m resolution as the standard for landslide susceptibility maps. In the following sub-section, we provide a brief overview of the causal factors that we used in this study. We classify these factors into several classes primarily using Jenks Natural Break method in ArcGIS 10.7, unless otherwise mentioned.

#### 3.2.1. Elevation

A change in elevation can bring changes in geomorphology, vegetation, and rate of erosion in an area and thus alters the landslide susceptibility [8]. We derived elevation from ASTER, SRTM, ALOS PALSAR, and SOB DEMs (Figure A1a–d of Appendix A) and divided them into five classes (see Table A1 of Appendix D).

#### 3.2.2. Slope

The slope is one of the most critical factors of landslides. Generally, with the increase of slope, shear stress increases, and therefore landslide susceptibility increases [43,44]. Like elevation, we derived slopes from four different DEMs using the slope tool in ArcGIS 10.7 (Figure A1e–h of Appendix A) and divided them into five classes using the Jenks natural break method (see Table A1 of Appendix D). We found different maximum slope values for different DEMs-ASTER (51.89º), SRTM (61.24º), ALOS PALSAR (65.36º), and SOB (46.4º). As these values were different for the same study area, the five classes of slopes (see Table A1 of Appendix D) were different.

#### 3.2.3. Aspect

The direction of precipitation, sunlight, and wind depend on aspect, and therefore it has effects on the growth of vegetation, rate of erosion, and thickness of soil [45]. From the DEMs, four aspect maps (Figure A1i–l of Appendix A) were prepared and were divided into ten classes (see Table A1 of Appendix D).

#### 3.2.4. Plan Curvature and Profile Curvature

Curvature is the rate of change of slope over time in an area. We used the four different DEMs to produce four plan and profile curvature maps (Plan: Figure A2a–d; Profile: A2e–h of Appendix B). Profile and plan curvatures were divided into three classes: concave, convex, and flat. Among these classes, concave slopes are more prone to landslides because water cannot disperse equally on these slopes [5].

#### 3.2.5. Topographic Wetness Index and Stream Power Index

Topographic wetness index (TWI) increases with the decrease of the slope; therefore, it is inversely related to landslide susceptibility [4,43]. Stream power index (SPI) represents the erosion power of streams. SPI is directly related to slopes; in a steeper slope, SPI will be higher, representing more erosion power while in a flat alluvial plain, SPI is low [45]. TWI and SPI maps were derived from four DEMs using Equation (1) and Equation (2) (TWI: Figure A2i–l of Appendix B; SPI: Figure A2m–p of Appendix B)

$$\text{TWI} = \text{Ln}\left(\frac{\text{A}}{\tan \alpha}\right) \tag{1}$$

$$\text{SPI} = \text{A} \times \tan{\alpha} \tag{2}$$

where, A = Area of a specific catchment and α = Slope gradient of the specific area. We divided TWI and SPI into five classes (Table A1 of Appendix D)

#### *3.3. Rainfall*

The intensity and duration of rainfall controls the initiation of landslides [44]. We used the mean annual rainfall of five weather stations of Bangladesh Meteorological Department (BMD) to prepare the rainfall map using the Kriging interpolation method (Figure A3a of Appendix C). We later divided it into five classes (Table A1 of Appendix D).

#### *3.4. Distance-Based Causal Factors*

Distance from the road networks, drainage networks, and fault lines were the three distance-based causal factors in this study. We used the Euclidean distance tool in ArcGIS 10.7 to derive the distance of landslides from the targeted features: road, drainage, and fault lines (Figure A3b–d of Appendix C) and divided the distances into five classes (Table A1 of Appendix D). Distance from the road networks is one of the most critical factors. The undercutting of slopes during road construction and the vibration created by vehicles damage the slope stability [43]. Drainage network indicates the zone of erosion in an area, and erosion is indirectly linked with the landslide susceptibility [4]. Fault lines indicate the geomorphological discontinuity in an area. Near the fault lines, the shear strength of rock is minimum. Therefore, areas near to the fault lines are prone to landslides [44].

#### *3.5. Normalized Di*ff*erence Vegetation Index (NDVI)*

NDVI indicates the growth of vegetation and biomass of an area [46]. Generally, the probability of the occurrence of the landslide on the naturally vegetated surface is lower than the bare lands [8]. We used Landsat 8 level 2 imagery of 11/10/2017 to prepare the NDVI map (Figure A3e) and divide it into five classes (Table A1 of Appendix D).

#### *3.6. Geology*

The strength of rock and soil permeability depends on the geology of an area. Therefore, geology has an impact on landslide susceptibility [43]. In this study, we used the geological map (Figure A3f) (1:100,000) of the Bangladesh Geological Survey (BGS). There are eight types of geologic formation found here: Dihing and Dupi tila formation; Bokabil formation; Bhuban formation; Tipam Sandstone; Valley alluvium and colluvium; Dihing formation, Girujan clay, and waterbodies.

#### *3.7. Land Use land Cover*

Land use/land cover and land use/land cover change are the two most crucial landslide causal factors in the study area [35,47]. Ahmed [31] and Rabby and Li [33] found that the rate of change of land use/land cover in our study area is high compared to other adjacent areas. In this study, we used two Landsat imageries to analyze the land use/land cover change (Landsat 5: date of Acquisition: 24/12/1998; Landsat 8: Acquisition Date: 29/11/2018). We used supervised maximum likelihood classifier to classify the 1998 and 2018 images into four land use classes: bare land, vegetation, built-up, and water bodies (Figure A3g of Appendix C). We later employed post-classification change detection techniques to analyze the land use/land cover changes between 1998 and 2018 (Figure A3h of Appendix C).

#### **4. Methodology**

To compare the effects of four different DEMs: ASTER, SRTM, ALOS PALSAR, and SOB on landslide susceptibility maps, we used a bivariate method: modified frequency ratio (MFR), a multivariate method: logistic regression (LR) and a machine learning method: random forest (RF). We assessed two scenarios: (a) considering only DEM-based seven causal factors in the models, and (b) considering all 15 causal factors, including the DEM-based factors, in the models. As we used three methods: MFR, LR, and RF on four different DEMs under two scenarios, therefore, the outcome would be twenty-four landslide susceptibility maps. For legibility, we used different acronyms in later sections (see Table 2).

**Table 2.** Acronyms used for different models.

#### *4.1. Training and Validation Dataset*

We divided the 261 landslide locations randomly into training (75%) and validation (25%) datasets. Bivariate MFR is a one-class classification method where non-landslide locations or absence of landslides are not required [48]. On the other hand, for multivariate LR and machine learning-based RF, the selection of non-landslide locations (pseudo-absence points) is essential [4]. Any place that does not have landslide can be considered as non-landslide. We randomly selected 261 non-landslide locations (Figure 1) (pseudo absence points) from the study area [49]. We split these non-landslide locations into training (196) and validation (65) data sets. In total, we had 392 (196: landslide locations; 196: non-landslide locations) data points for training and 130 (65: landslide locations; 65: non-landslide locations) data points for testing the LR and RF models.

#### *4.2. Modified Frequency Ratio (MFR)*

MFR is an improved version of the widely used frequency ratio (FR) method [35,50]. Lee and Talib [51] proposed FR, which assesses the spatial relationship between the landslide locations and the landslide causal factors [52]. In the FR method, each of the causal factors must be divided into subclasses or categories; for example, in this study, we split slope into five categories using the Jenks natural break method (Table A1 of Appendix D). We calculated the FR values using Equation (3), and later these FR values of each of the subclasses of causal factors were used in the MFR model.

$$\text{FR}\_{\text{iĭ}} = \frac{\text{N}\_{\text{iġ}}/\text{N}}{\text{M}\_{\text{i}}/\text{M}} \tag{3}$$

where FRij = Frequency Ratio of jth Subclass of Factor i

Nij = Total area of the landslide pixels within the jth subclass of factor i

N = Total area of landslide pixels in the study area

Mi j = Total area of the pixels in the jth subclass of factor i

M = Total area of the study area

FR > 1 means association of landslides with that subclass. In other words, there is a probability of occurrences of landslides in that subclass. FR < 1 means no association [53].

For calculating the MFR, then we normalized the FRs using Equation (4).

$$\text{Rf}\_{\text{ij}} = \frac{\text{FR}\_{\text{ij}}}{\sum \text{FR}\_{\text{i}}} \tag{4}$$

where Rfij = Relative frequency of jth subclass of factor i

FRij = Frequency ratio of the jth subclass of factor i

P FRi = Sum of the frequency ratios of factor i

Later, we calculated the prediction rate (PR) using Equation (5). In the FR model, the overall contribution of causal factors to the occurrence of landslides is not measured. Only the subclass wise contribution is measured [6,50]. In MFR, we can measure the overall contribution because PR indicates overall association of a causal factor. The lowest value of PR is 1 and the higher the PR value the stronger is the association of causal factor with the landslides [6].

$$\text{PR}\_{\text{i}} = \frac{(\text{MaxRf}\_{\text{i}} - \text{MinRf}\_{\text{i}})}{(\text{MaxRf}\_{\text{i}} - \text{MinRf}\_{\text{i}})\text{min}} \tag{5}$$

where PRi = Prediction rate of factor i

MaxRfi = Maximum relative frequency of factor i

MinRfi = Minimum relative frequency of factor i

(MaxRfi − MinRfi ) min = Lowest difference between maximum and minimum relative frequency of all the factors

To calculate the landslide susceptibility ondex (LSI) and to produce the landslide susceptibility maps Equation (6) was used

$$\text{LSI} = \sum\_{i=1}^{n} \text{Rf}\_{\text{ij}} \times \text{PR}\_{\text{i}} \tag{6}$$

where LSI = Landslide susceptibility index

Rfij = Relative frequency of jth subclass of factor i

PRi = Prediction rate of factor i

For landslide susceptibility mapping, we used the reclassify tool in ArcGIS 10.7 to reclassify the categories of causal factors according to the Rf values. Later, we multiplied each of the reclassified raster layers with the prediction rates and summed up to produce the final landslide susceptibility maps.

#### *4.3. Logistic Regression (LR)*

Logistic regression (LR) is one of the most widely-used multivariate statistical methods in landslide susceptibility mapping [4,6,54–56]. An LR model predicts the presence of landslides using the binary landslide data (presence and absence of landslides or landslides and non-landslides) and their relationship with the landslide causal factors [56,57]. Here, landslide and non-landslide locations are dependent variables, and causal factors are independent variables. These independent variables can be numerical or categorical [4]. Equation (7) is used in LR model

$$\text{Logit}\left(\mathbf{Y}\right) = \ $1\mathbf{X}1 + \$ 2\mathbf{X}2 + \ $3\mathbf{X}3 + \mathbf{....} + \mathbf{....} + \mathbf{....} + \mathbf{....} + \dots + \mathbf{<\color{red}{251}{$ \mathbf{X} \ $}} + \$ 3\mathbf{X}i + \mathbf{e} \tag{7}$$

where Y = The presence of landslides

Xi = ith Causal factor ßi = Regression coefficient of the ith causal factor e = Error We used Equation (8) to determine the probability.

$$\mathbf{P} = \frac{\exp \mathbf{Y}}{1 + \exp \mathbf{Y}} \tag{8}$$

We used the R software environment for the forward stepwise LR method. We multiplied the raster layers of the statistically significant causal factors with the coefficients and summed up using Equation (7) in the R software environment. Finally, we used Equation (8) to produce the landslide susceptibility maps.

#### *4.4. Random Forest Classification (RF)*

The random forest method was developed by Breiman [58] and is an ensemble learning method [19]. Lately, the use of the RF method for landslide susceptibility mapping has increased due to its high performance in predicting landslide locations [37,59].

This method uses bootstrapping techniques to generate a bunch of classification trees based on subsets of observations [27]. There is high variance among the individual trees, and therefore classification based on a single tree is unstable and prone to overfitting [37]. Random forest is improved over commonly used tree-based methods, such as a decision tree or bagged tree because it decorrelates the trees. RF uses ensembles of trees and lets each tree define the class membership, and finally, the respective class is assigned based on the highest votes [27,37]. Since the bootstrapping method is used, a set of data is not used in the model training stage and this set of data is known as out-of-bag (OOB) [27]. These OOB data are used to calculate the mean decrease of accuracy and Gini coefficient [37]. The accuracy and Gini coefficient are used in variable selection and ranking [19]. They also provide the statistical weights or variable importance of each of the predictors used in the model [27].

There are several advantages of using RF methods, such as rescaling and transformation of data are not essential; missing data and outliers can be ignored [60]. Moreover, it can deal with both numerical and categorical data, and the use of a dummy variable is not required [37]. In this study, we developed the Random Forest model in the R software environment using the "randomForest" package [61].

#### *4.5. Multicollinearity Diagnostics*

In the LR model, multicollinearity can bring inaccuracies in variance and unsuitability in estimates. On the other hand, in the RF model, it can affect the variable importance [62]. Therefore, multicollinearity diagnostics: variance inflation factor (VIF) and tolerance were used before using the causal factors in LR and RF models. Since VIF values were <10 and tolerance were <0.3, causal factors were independent and were used in these two models [63]

#### *4.6. Model Validation and Comparison*

Success and prediction rate curves were used to test the performance of the susceptibility models. The training data set was used for the success rate curve, and the validation data set was used for the prediction rate curve [6]. "The area under the curve" (AUC) of success rate indicates how well the model fits the training data while the AUC of prediction rate suggests how well the model will predict future landslides [15,64]. AUC value ranges from 0–1 or 0%–100% and it can be grouped into the following categories: 0.50–0.60 (fail); 0.60–0.70 (poor); 0.70–0.80 (fair); 0.80–0.90 (good), and 0.90–1.00 (excellent) [53]. We also used two non-parametric tests: Friedman and Wilcoxon Signed Rank test to assess whether there are any significant differences in performances between the susceptibility [65–68]. Friedman's test is used for multiple comparisons. This test determines whether there is any significant difference in performance in multiple models [29], while the Wilcoxon Signed Rank test is used for pairwise comparison of susceptibility models and, therefore, can indicate which models are significantly different [68].

However, statistical performance assessments such as success and prediction rate curves cannot show the level of agreement among the models. Therefore, we used convergent validation through the coverage based cross-comparison [69,70]. The MFR gives landslide susceptibility index while LR and RF provide the probability of landslides. We reclassified them into five susceptibility zones: very low, low moderate, high, and very high using the Jenks natural break method. Later, we used the raster calculator in the ArcGIS platform to subtract the reclassified model from one another. The outcome can be any integer, but only "zero" will indicate the areas which were classified into the same susceptibility zones by two compared models. We calculated the percentage of area under this "zero" class and this percentage indicates the spatial convergence or agreement between two compared models [69].

#### **5. Results**

#### *5.1. Susceptibility Assessment Using Modified Frequency Ratio (MFR)*

We found that the prediction rates (PRs) of the slope, aspect, and elevation derived from three DEMs (ASTER, SRTM, and ALOS PALSAR) are similar. The PRs for these three factors are around 2.25, 1.0, and 3.0, respectively (see Table A1 of Appendix D). Compared to these three DEMs, the SOB DEM showed different PRs. Our analysis revealed that for slope, aspect, and elevation, the SOB DEM-based PRs are 1.79, 2.37, and 2.41, respectively. We further found that with the increase of slope, the probability of landslide increases. The relatively safer zones are in areas below 8º slope where frequency ratio (FR) < 1. For ASTER, SRTM, and ALOS PALSAR DEMs, we found that the slope class 14–23º had the highest probability of landslides. However, for SOB DEM, the highest probability of landslide is in the 8–14º slope class. In the case of TWI, we found that ALOS PALSAR DEM has the highest PR (4.30) followed by ASTER DEM (PR = 3.19) and SOB DEM (PR = 1.61). For all the four DEMs, the probability of landslides was higher in areas where TWI is less than 6. The SPIs derived from four DEMs have lower PRs compared to other causal factors and the class-wise weight (FR values)

showed the same sort of pattern. We further found that for plan curvature, ALOS PALSAR has the highest PR (4.30), while for profile curvature, SRTM has the highest PR (4.68).

In general, we found no specific pattern of FR and PR values for these seven factors. We observed that causal factors derived from either ALOS PALSAR or SRTM have higher PR values, and SOB has the lowest PR among the four DEMs. The causal factors derived from different DEMs did not have a significant impact on the FR and PR values of MFR. This finding is similar to the findings of Chang et al. [29]. Most of the topographic factors are the first derivative of the DEMs other than the TWI and SPI. Thus, for TWI, DEMs have more substantial impacts than other factors. It is because there can be a small difference among the DEMs, and when the second derivative is used, these become pronounced [29].

The class-wise FR values of the eight factors that are not derived from the DEMs are generally similar; however, the PR values are different (see Table A1 of Appendix D). The causal factors derived from the DEMs played a crucial role in determining the PR values for these eight factors. For ASTER (0.141), SRTM (0.134), and ALOS PALSAR (0.139), the lowest difference was between the maximum and minimum relative frequency of aspect. While for SOB (0.160), it was SPI. Since these values are slightly different, it affected the PR values.

#### Landslide Susceptibility Maps (MFR)

We produced Landslide Susceptibility Indices (LSIs) of the four MFR models using Equation (6) for DEM-based causal factors. The LSI of MFR\_ASTER\_DEM ranged from 994.6 to 8388.1. The LSIs for MFR\_SRTM\_DEM LSIs for MFR\_SRTM\_DEM; MFR\_ALOS\_DEM and MFR\_SOB\_DEM ranged from 591.4 to 9458.8 and 527.8 to 10056.5 and 638.8 to 6130.0, respectively. LSI does not have a unit. It is the product of relative frequency and prediction rate Equation (6), and both of these do not have units. The greater the LSI, the greater is the landslide susceptibility, and the smaller the LSI value, the lower is the susceptibility [35,44]. The ranges of LSIs indicate that MFR\_SRTM\_DEM and MFR\_ALOS\_DEM models had a comparatively broader range than the rest of the two models. This happened because of the variable FR and PR values. For ASTER DEM, the highest PR value was for plan curvature, while for other DEMs it was for profile curvature (see Table A1 of Appendix D). For all DEM-based factors, SOB had the lowest PRs among the four models and therefore it affected the LSIs. Later, we used the same Equation (6) to produce four MFR models based on 15 causal factors. The LSI of MFR\_ASTER DEM ranged from 1613.00 to 20370.10. The LSIs for MFR\_SRTM, MFR\_ALOS, and SOB DEMs ranged from 1314.40 to 22300.34 and 1234.95 to 22180.24 and 1995.7 to 17316.9, respectively. The LSIs of MFR\_SRTM and MFR\_ALOS DEMs have a comparatively broader range than the rest of the two models. The highest PR value for ASTER was 6.25 and for SOB, it was 5.64 (Table A1 of Appendix D) for distance to the road network, while for SRTM and ALOS PALSAR, it was 6.61 and 6.38, respectively. For other causal factors (Table A1 of Appendix D), SOB had the lowest PRs among the four, and therefore LSI was the lowest. As mentioned earlier, the FR values varied for seven topographic factors derived from four different DEMs, and these factors had impacts on the PR of eight common factors. Ultimately these variations defined the LSI of the susceptibility maps.

Since different models had different LSIs, Rescale by Function tool in ArcGIS was used to normalize the LSIs into a 0.0–1.0 scale. Later, we used the Jenks natural break method to classify the normalized LSIs into five susceptibility zones: very low, low, moderate, high, and very high. Generally, the spatial appearances of the landslide susceptibility maps have similarities with the map of causal factors that have a higher contribution to landslides. In this study, this contribution is shown by PR and FR values (see Appendices A–C). We found that the spatial appearance of seven causal factors derived from SOB was different from the spatial appearance of seven causal factors derived from ASTER, SRTM, and ALOS PALSAR. ASTER, SRTM and ALOS PALSAR based susceptibility maps (Figure 2a–c) show a comparatively lesser percentage of area as very low or low susceptibility zones than the SOB-based landslide susceptibility map (Figure 2d). However, the SRTM based map shows comparatively more areas as high and very high susceptibility zones that the other maps. We found when all factors were

considered, areas near to the road network are highly susceptible to landslides (Figure 3a–d), because the PR values of the distance from the road networks were the highest.

‐ **Figure 2.** Landslide susceptibility maps produced using Modified Frequency Ratio (MFR), Logistic Regression (LR) and Random Forest (RF) models (seven DEM-based factors).

**Figure 3.** Landslide susceptibility maps produced using MFR, LR, and RF models (all 15 factors).

#### *5.2. Susceptibility Assessment using Logistic Regression (LR)*

‐ For the DEM-based factors, the LR model detected two to five statistically significant factors (see Table A2 of Appendix E). Elevation and slope were the two common statistically significant causal factors for ASTER, SRTM, and ALOS PALSAR based models. Since the ALOS PALSAR based model, the highest number of causal factors was chosen, DEM had the highest impact on the landslide susceptibility map. Odds ratio (Table A2 of Appendix E) shows that slope was the most important factor of landslides for ASTER and SRTM based models, while aspect came out as the most crucial factor for ALOS PALSAR and SOB based models.

When all (15) causal factors were used, we found a total of four to eight statistically significant causal factors (see Table A2 of Appendix E). Slope, elevation, SPI, and aspect were the significant DEM-based causal factors for LR\_SRTM and LR\_ALOS based models. For LR\_SOB, aspect was the statistically significant DEM-based causal factor. When 15 causal factors were used, the model assessed the interaction of DEM-based factors with the common eight factors. Therefore, when only the DEM-based causal factors were used, some factors came out statistically significant (e.g., SPI for

LR\_ASTER\_DEM). When all the factors were used in the LR\_ASTER based model, SPI came out insignificant. For LR\_SRTM and LR\_ALOS, three DEM-based causal factors were detected as statistically significant; therefore, for these two landslide susceptibility maps, DEM would have more impact than the LR\_ASTER and LR\_SOB susceptibility maps.

#### Landslide Susceptibility Maps (LR)

We used the Jenks natural break method to classify the probability of landslides into five zones: very low, low, moderate, high, and very high. The spatial appearances of the LR\_SOB\_DEM model (Figure 2h) have a different appearance than the other three maps. LR\_SRTM\_DEM (Figure 2f) and LR\_ALOS\_DEM (Figure 2g) have an almost identical spatial appearance. While in LR\_ASTER\_DEM models, the slope had a higher coefficient (ß = 0.26) (Table A2 of Appendix E) than the SRTM (ß = 0.14) and ALOS PALSAR ((ß = 0.05). Therefore, in the LR\_ASTER\_DEM map, the areas with steeper slopes, mainly in the mid-north and mid-south of the study area, were classified as high to very highly susceptible. In LR\_SRTM\_DEM and LR\_ALOS\_DEM maps, these areas were classified either as moderate or high susceptibility zones. In LR\_SOB\_DEM, only elevation and aspect were two significant factors. Therefore, the susceptibility map took the shape of the map of these two factors (see Appendix A).

When all (15) causal factors were used, the spatial appearance of LR\_SOB (Figure 3h) was different from the other three maps (Figure 3e–g). LR\_SOB map was influenced by the distance from the fault lines, distance from the road networks, and land use/land cover. Although aspect was a significant factor, the coefficient value of aspect was similar to other significant factors, and distance from the fault lines (ß = 1.07) (Table A2 of Appendix E) had a higher coefficient value than aspect. Therefore, most of the study area was classified as low or very low susceptibility zones. In the LR\_ASTER model (Figure 3e), the slope had the highest coefficient (ß = 0.31) (Table A2 of Appendix E), and therefore, areas with steeper slopes were classified as high or very high susceptibility zones. But in LR\_SRTM and LR\_ALOS models, the slope had comparatively lower coefficient values than ASTER. As a result, some areas were classified as moderate susceptibility zones in these two maps (Figure 3f–g). In LR\_SRTM and LR\_ALOS susceptibility maps, common factors such as distance from the road networks and fault lines, land use/land cover, and land use/land cover change did not have higher coefficient values than the DEM-based causal factors. That is why, unlike LR\_SOB, the spatial appearance of the susceptibility maps did not follow the appearance of the maps of common factors.

#### *5.3. Susceptibility Assessment using Random Forest (RF)*

RF\_ASTER\_DEM, RF\_SRTM\_DEM, and RF\_ALOS\_DEM models detected (Figure 4a) slope and RF\_SOB\_DEM model identified the aspect as the most critical factor. When all 15 factors were used in the models, RF\_ASTER and RF\_SRTM (Figure 4b) detected slope and the rest of the two models detected distance from the road network as the most important causal factor. For RF\_ASTER, RF\_SRTM, and RF\_ALOS, DEM-based factors such as elevation, TWI, and aspect had higher importance in the models than the common factors. But in the RF\_SOB model DEM-based factors had less importance than the common factors. There is no similarity among the models in detecting the importance of DEM-based causal factors. For example, In RF\_SOB, the slope was ranked as one of the least important factors, but for other models, it was ranked as the most important factor (Figure 4b).

‐

‐ ‐

‐

 ( **Figure 4.** Variable importance plots of random forest model: (**a**) Digital Elevation Model (DEM)based causal factor (7) used in the models. (**b**) All (15) factors used in the models. Pl=plan curvature; Pr= profile curvature; LC= land use/land cover; LLC= land use/land cover change; DF= distance from the fault lines; DD= distance from the drainage networks; DR= distance from the road networks.

Landslide Susceptibility Maps (RF)

Like MFR and LR, we used the same method to classify the probability of landslides into five susceptibility zones. The spatial appearance of the RF\_SOB\_DEM susceptibility map (Figure 2l) was different from the susceptibility maps of the other three models (Figure 2i–k). For, RF\_ASTER\_DEM (Figure 2i) areas in the mid-north to mid-south were classified as high or very high susceptibility zones. While in RF\_SOB\_DEM and RF\_ALOS\_DEM the same areas were classified as moderate susceptibility zones. In RF\_ASTER\_DEM, slope was the most critical factor. Similarly, in RF\_SRTM\_DEM and RF\_ALOS\_DEM slope was the most crucial factor, but in these two models, the contribution of the slope (Figure 4a) to the model is lesser than the RF\_ASTER\_DEM model. In RF\_ASTER\_DEM, the difference of variable importance between slope and other factors was comparatively higher than the other models, the effect of slope on the susceptibility map was visible.

RF\_ASTER; RF\_SRTM, RF\_ALOS; and RF\_SOB models (Figure 3i–l), spatial appearances were different from each other. Since in the RF\_SOB model, distance from the road networks was the most crucial factor, areas near to roads were classified as high or very high susceptibility zones. Distance from the road network was not ranked as the most critical factor in the other three models.

#### *5.4. Validation and Comparison of Landslide Susceptibility Maps*

#### 5.4.1. Success and Prediction Rate Curves

#### DEM-Based Causal Factors

When only DEM-based seven causal factors were used for MFR models, among all DEMs the MFR\_SRTM\_DEM model gave the superior performance for both success (AUC = 80.73%) and prediction (AUC = 77.37%) rate curves (Figure 5a,b). The MFR\_SOB\_DEM model appears to perform the weakest in assessing success and prediction. The AUCs of success rate curves (Figure 5a) showed that MFR\_SRTM\_DEM falls under the good category while MFR\_ASTER\_DEM and MFR\_ALOS\_DEM fall under the fair category. But MFR\_SOB\_DEM falls under the poor category. AUCs of prediction rate curves (Figure 5b) show that all models other than the MFR\_SOB\_DEM gave fair performances.

‐ **Figure 5.** Success and prediction rate curves (seven DEM-based factors): (**a**) MFR (success); (**b**) MFR (prediction); (**c**) LR (success); (**d**) LR (prediction); (**e**) RF (success); (**f**) RF (prediction).

For LR, LR\_ALOS\_DEM outperformed the other three models (Figure 5c,d). LR\_SOB\_DEM presented the weakest performance among the four models and thus fell under the fail category. The AUCs of success and prediction rates (Figure 5c,d) show that the other three models are under the fair category. For RF models, we got similar results as the LR model. RF\_ALOS\_DEM outperformed other models, and RF\_SOB\_DEM was the least accurate model. RF\_SRTM\_ALOS\_DEM and RF\_SRTM\_DEM gave an almost similar performance.

We found that when all 15 causal factors are used, different models showed variable performances. For the MFR model, the use of 15 causal factors decreases the predictive performance on an average by 5% of landslide susceptibility maps based on ASTER, SRTM, and ALOS PALSAR (see Figure 6a,b).

models, DEM-based causal factors can give better prediction performance and use of non-DEM-based factors can reduce the accuracy. Inclusion of more DEM-based causal factors and more landslide locations may increase the accuracy of the models in the study area. ‐

**Figure 6.** Success and prediction rate curves (15 factors): (**a**) MFR (success); (**b**) MFR (prediction); (**c**) LR (success); (**d**) LR (prediction); (**e**) RF (success); (**f**) RF (prediction).

For the LR model, the use of 15 factors increased the accuracy of the model. For three global DEMs success rates increased by around 5.0% but for SOB it increased by 21.7%. On the other hand, prediction rates showed the same trend as for three global DEMs the increase of performance was around 3.5% but for SOB it was 17.0%. It proves that the use of common factors increased the accuracy substantially for SOB DEM.

Like the LR model, for RF models, the use of 15 causal factors increased the accuracy of the model. For ASTER and SRTM the increase of the success rate was around 7%. For ALOS PALSAR the success rate increased by 12.4%. Here again, SOB had the highest increase (22.6%) in success rate. The increase in prediction rate was not as high as the success rate. For three global DEMs prediction rates increased by around 2%–3% and for SOB the prediction rate increased by 10.1%. Machine learning algorithms such as the random forest learn the behavior or the training data. Therefore, the increase of success rate due to inclusion of new variables was high. Since machine learning algorithms learn the behavior of the training data it fails to predict the validation or unknown data [60]. Therefore, in our study, the increase of prediction rate is around 50% lower than the increase of prediction rate for RF models.

#### 5.4.2. Spatial Comparison of Landslide Susceptibility Maps

Spatial convergence indicates how much area is classified into same susceptibility zones. When seven DEM-based factors were used in the MFR model, MFR\_SOB\_DEM had 30% of spatial convergence while the remaining DEMs had around 40% (Table 3). As we discussed before, the landslide

susceptibility maps of MFR\_SOB\_DEM has a different spatial appearance (Figure 2d) than the rest of the three susceptibility maps and these results (Table 2) support the previous discussion of our study. For the LR models, Table 3 shows a similar trend. Spatial convergences of the LR\_ASTER\_DEM, LR\_SRTM\_DEM, and LR\_ALOS\_DEM were around 44%, while the LR\_SOB\_DEM showed approximately 19% of spatial convergence. For the RF models, Table 3 shows the similarities with the findings of MFR and LR.

**Table 3.** Spatial comparison and convergence analysis of landslide susceptibility maps.

When all factors were considered for modeling, spatial convergence between the DEMs (Table 3) increased around 40% for MFR models. While for the LR and RF models, the spatial convergence was approximately 25% and 12%, respectively. In the MFR model, all causal factors were used, while in the LR model, significant causal factors were used and in the RF model, 2–3 causal factors, for example, profile and plan curvatures had comparatively low or no variable importance in the model.

The results of Friedman tests (Table 4) show that in both the scenarios (a. seven DEM-based causal factors used, and b. 15 causal factors used) for MFR and RF models *P* < 0.05. It means at least one of the landslide susceptibilities models had significantly different performance than the rest of the models. While for LR models, when seven DEM-based causal factors were used, at least one of the models was statistically different in performance than the rest of the models. When all factors were used in LR models, there was no statistically significant difference in performance between the models. Wilcoxon signed-rank test conducted the pairwise comparison. The results (Table 4) show that in scenario two, the performances of landslide susceptibility maps produced using SRTM and ALOS PALSAR did not have a statistically significant difference. Other than that, all the performances of the MFR based landslide susceptibility maps were statistically (α = 0.008 after Bonferroni correction) different from each other. When seven causal factors were used for LR models, the performance of SOB based models was significantly different (Table 5) from all other models. But when all factors are used, these differences become insignificant. It indicates that eight common factors overshadow the effect of SOB based causal factors. In the MFR model, it did not happen since it did not consider the interaction of causal factors.

**Table 4.** Result of the Friedman Test for landslide susceptibility maps.

In the case of RF models, RF\_ALOS\_DEM (Table 5) was statistically different from the other models. But when all causal factors were used the difference of performance became insignificant for RF\_SRTM model. RF used a more complex algorithm than the LR and MFR models. Therefore, the Wilcoxon Signed-Rank test gave different results for RF models than MFR and LR models.

**Table 5.** Comparison of landslide susceptibility maps based on three different DEMs using Wilcoxon signed-rank test.

\* = Significant after Bonferroni correction (*P* < 0.008).

#### **6. Discussion**

This paper evaluates the suitability of three available global DEMs: ASTER, SRTM, and ALOS PALSAR and a local DEM: SOB for landslide susceptibility mapping in Rangamati district, Bangladesh. Causal factors derived from ASTER and ALOS PALSAR DEM have been used in landslide susceptibility mapping in different parts of the Chittagong hilly areas, Bangladesh [7,9,31–35]. Since the study areas of these studies were different, we could not compare them to find out which DEM gives the best accuracy in the prediction of landslide susceptibility [36]. Our study showed that three global DEMs outperformed the local DEM in all three landslide modeling scenarios.

In the first scenario, only DEM-based causal factors were used in modeling and for MFR models, MFR\_SRTM\_DEM outperformed the other three models. The difference of AUCs of both the success and prediction rate curves between \_SRTM\_DEM and ALOS\_DEM was around 3%, indicating a similarity in predictions. For the processing of ALOS PALSAR DEM, SRTM GL1 data is used for radiometric correction [71]. Therefore, the quality of ALOS PALSAR DEM depends on the quality of SRTM DEM. SOB DEM-based MFR did not show a good performance for the study area and the prediction performance can be improved when a more representative landslide inventory with more landslide locations is used. We found that the ASTER DEM-based MFR model showed a weaker performance than the other two open-source global DEMs. Our result is consistent with other studies that utilized ASTER DEM [29,72,73]. It may happen because ASTER DEM contains many artifacts such as the presence of peaks in the flat terrain, and it ultimately affects the landslide susceptibility map [29]. The poor performance of the SOB DEM can be attributed to the interpolation methods that were used to extrapolate elevations from the available spot heights in the hilly parts of Bangladesh [74]. It affected the accuracy and quality of DEM in the hilly parts of Bangladesh. On the other hand, for LR and RF, ALOS PALSAR based models outperformed the rest of the models. Here, again the difference between LR\_ALOS\_DEM and LR\_SRTM\_DEM was low for success and prediction. In all cases, SOB based models gave the worst performance and causal factors derived from SOB DEM cannot explain the landslide susceptibility of the study area. For example: In LR\_SOB\_DEM, the LR model used two significant causal factors: elevation and aspect (Table A2 of Appendix E). The low coefficient (ß) values of these two factors indicate that these two causal factors cannot adequately explain the landslide susceptibility of the study area.

In the second, scenario, for MFR models, the use of 15 causal factors increased the prediction accuracy for MFR\_SOB. But for the other three models, it reduced accuracy. It indicates that causal factors derived from three global DEMs were capable enough to explain the landslide susceptibility of the study area. MFR is a bivariate model, and it does not consider the interaction of the causal factors. Moreover, unlike LR and RF models, it does not require non-landslide (pseudo absence point) in modeling [44]. When all causal factors were used, PRs of some of the common causal factors were comparatively higher than the PRs of the DEM-based causal factors. For example, PRs of distance from the road networks were around 6.00 (Table A1 of Appendix D), while PRs of the DEM-based causal factors ranged from 1.00–4.43. Therefore, the higher PRs of the common causal factors overshadowed the DEM-based factors [6,35,44]. On the other hand, the quality of the SOB DEM was poor and was not capable of explaining the landslide susceptibility of the study area. That is why, in the MFR\_SOB model, the prediction accuracy increased, but the increase was very low. For LR and RF models, in second scenarios, the prediction accuracy increased by 2%–3% for global DEM-based models, while for SOB based models, it exceeded 10%. Both LR and RF models consider the interaction of causal factors, and therefore, in these two models, the effects of common factors on prediction performance were revealed better than the MFR models.

The findings of this study have similarities with other research where the suitability of different global and local DEMs was evaluated for landslide susceptibility mapping [26–28]. In most of the studies, ASTER DEM-based landslide susceptibility maps were outperformed by the other global DEM-based landslide susceptibility maps. On the other hand, local DEM-based landslide susceptibility maps have better prediction accuracy than the global DEMs [1,29]. In these studies, local DEMs

were mainly light detection and ranging (LiDAR)-based DEMs, which generally have very high spatial resolution compared to global DEMs [29]. SOB DEM was prepared under the project titles "Improvement of Digital Mapping System of Survey of Bangladesh", where the main aim was to prepare a 1:5000 scale DEM of major cities in Bangladesh. This SOB project prepared DEM for the hilly areas of Bangladesh, including our study area using the local spot heights. Different interpolation methods were used to prepare 25m DEM from these local spot heights [73]. Therefore, in hilly areas of Bangladesh, the quality of SOB DEM is not good enough, and the application of this DEM in geomorphological studies such as landslide susceptibility mapping can give questionable results similar to what we got in our study. Our study also pointed out the importance of the preparation of very high-resolution DEMs using LiDAR for the hilly areas of Bangladesh. It will help in various geomorphological studies, including landslide susceptibility mapping. As we did not find a substantial difference among global DEMs, any global DEM can be utilized for landslide susceptibility mapping in Bangladesh in the absence of very high-resolution DEMs.

For an in-depth study, we utilized non-parametric tests. Chang et al. [30] used non-parametric tests to detect a significant difference in performance in landslide susceptibility maps prepared using different DEM-derived causal factors. Their study did not find any significant difference in performance for two machine learning methods: RF and support vector machines. In our research, we did not see any significant difference in performance for RF and LR based models. But for the MFR model, we found a significant difference in performance. It indicates that the effect of DEM-based causal factors on the performance of bivariate models is more than on the multivariate and machine learning models. Thus, we suggest using multivariate and machine learning methods and any one of the global DEMs for landslide susceptibility mapping in Bangladesh.

#### **7. Conclusions**

This paper assesses the effects of the DEM-derived causal factors on the landslide susceptibility maps produced using the bivariate (e.g., MFR), multivariate (e.g., LR), and machine learning (e.g., RF) models. In this study, we tested two scenarios: a. susceptibility assessment with only seven DEM-based causal factors; b. inclusion of other 8 causal factors along with DEM-derived factors. The success and prediction rate curves showed SRTM DEM outperformed under an MFR model in both scenarios. Our analysis revealed that ALOS PALSAR DEM provided the best prediction accuracy while using both LR and RF models in landslide susceptibility mapping. For all models and scenarios, the SOB DEM does not perform well compared to other DEMs.

The prediction accuracies of landslide susceptibility mapping using only DEM-derived factors is low compared to the utilization of all casual factors. Although SOB DEM has a poor performance in susceptibility assessment with only a DEM-derived factor, the accuracy is significantly improved when other non-DEM-based factors are added. Therefore, we argue that causal factors derived from the SOB DEM have limited influence in landslide susceptibility assessment for the study area. Besides, for the LR and RF models, the use of all the causal factors increased the prediction performance. Spatial convergence analysis showed that three global DEM-based models have similar accuracies and performed far better than SOB DEM-based models. Therefore, we recommend that SOB DEM should not be used for landslide susceptibility assessment in Bangladesh. Although ASTER DEM-based models showed the weakest performance among three global DEMs, the difference of performance was negligible. Therefore, we recommend any one of these three global DEMs can be utilized for landslide susceptibility assessment in Bangladesh.

Our study also highlights that, for the MFR model, DEM had the highest impact on the accuracy of landslide susceptibility assessment as Wilcoxon rank tests showed that the performance of susceptibility maps was significantly different. For the LR and RF models, the effect of DEM was less significant. We contend that while using the bivariate model, we must be careful about the quality of DEM. Two scenarios in this study helped to understand the impact of DEMs in landslide susceptibility assessment. Moreover, we used multiple metrics to evaluate the accuracies of susceptibility assessment

including AUC, spatial convergence analysis, and non-parametric test. Although the use of one performance measure is a common practice, the utilization of various measures helped this research to understand the impacts of DEM and makes the conclusion robust.

The main limitation of our study is that we cannot ascertain if adding more DEM-derived factors such as terrain roughness could improve the result. The impacts of DEM on landslide susceptibility maps are not universal and may vary from place to place. Therefore, we cannot conclude that a specific DEM can be better than others. Advanced machine learning and deep learning methods can be used to check whether these algorithms can reduce the differences in prediction performance of landslide susceptibility maps prepared using different DEM-derived causal factors.

**Author Contributions:** Conceptualization, Y.W.R; methodology, Y.W.R.; M.S.R.; analysis, Y.W.R.; writing, Y.W.R. A.I.; M.S.R.; review and editing, A.I. and M.S.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding

**Acknowledgments:** The authors thank Bayes Ahmed (Lecturer, Institute for Risk and Disaster Reduction, University College London) for his assistance in obtaining some datasets.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

**Figure A1.** Landslide Causal Factors: (**a**) Elevation (ASTER); (**b**) Elevation (SRTM); (**c**) Elevation (ALOS PALSAR); (**d**) Elevation (SOB); (**e**) Slope (ASTER); (**f**) Slope (SRTM); (**g**) Slope (ALOS PALSAR); (**h**) Slope (SOB); (**i**) Aspect (ASTER); (**j**) Aspect (SRTM); (**k**) Aspect (ALOS PALSAR); (**l**) Aspect (SOB).

#### **Appendix B**

**Figure A2.** Landslide Causal Factors: (**a**) Plan Curvature (ASTER); (**b**) Plan Curvature (SRTM); (**c**) Plan Curvature (ALOS PALSAR); (**d**) Plane Curvature (SOB); (**e**) Profile Curvature (ASTER); (**f**) Profile Curvature (SRTM); (**g**) Profile Curvature (ALOS PALSAR); (**h**) Profile Curvature (SOB); (**i**) TWI (ASTER); (**j**) TWI (SRTM); (**k**) TWI (ALOS PALSAR); (**l**) TWI (SOB); (**m**) SPI (ASTER); (**n**) SPI (SRTM); (**o**) SPI (ALOS PALSAR); (**p**) SPI (SOB).

#### **Appendix C**

**Figure A3.** Landslide Causal Factors: (**a**) Rainfall; (**b**) Distance from the Road Networks; (**c**) Distance from the Drainage Networks; (**d**) Distance from the Fault Lines; (**e**) NDVI; (**f**) Geology; (**g**) Land Use / Land Cover; (**h**) Land Use/ Land Cover Change.

#### **Appendix D**

**Table A1.** Spatial relationship between Causal Factors and Landslides.ASTER DEM; B= SRTM DEM; C= ALOS PALSAR DEM D= SOB DEM.

**Table A1.** *Cont.*

**Table A1.** *Cont.*

**Table A1.** *Cont.*

**Table A1.** *Cont.*

**Table A1.** *Cont.*

#### **Appendix E**

1 = Land use/land cover change; 2 = distance from the drainage networks; 3 = distance from the road networks;

4 = land use/land cover; 5 = distance from the fault lines.

#### **References**

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

# **Landslide Characterization Applying Sentinel-1 Images and InSAR Technique: The Muyubao Landslide in the Three Gorges Reservoir Area, China**

**Chao Zhou 1,2 , Ying Cao 3 , Kunlong Yin 3, \*, Yang Wang 3 , Xuguo Shi 1 , Filippo Catani 4 and Bayes Ahmed 5**

Received: 21 August 2020; Accepted: 13 October 2020; Published: 16 October 2020

**Abstract:** Landslides are a common natural hazard that causes casualties and unprecedented economic losses every year, especially in vulnerable developing countries. Considering the high cost of in-situ monitoring equipment and the sparse coverage of monitoring points, the Sentinel-1 images and Interferometric Synthetic Aperture Radar (InSAR) technique were used to conduct landslide monitoring and analysis. The Muyubao landslide in the Three Gorges Reservoir area in China was taken as a case study. A total of 37 images from March 2016 to September 2017 were collected, and the displacement time series were extracted using the Stanford Method for Persistent Scatterer (StaMPS) small baselines subset method. The comparison to global positioning system monitoring results indicated that the InSAR processing of the Muyubao landslide was accurate and reliable. Combined with the field investigation, the deformation evolution and its response to triggering factors were analyzed. During this monitoring period, the creeping process of the Muyubao landslide showed obvious spatiotemporal deformation differences. The changes in the reservoir water level were the trigger of the Muyubao landslide, and its deformation mainly occurred during the fluctuation period and high-water level period of the reservoir.

**Keywords:** landslide deformation; InSAR; reservoir water level; Sentinel-1; Three Gorges Reservoir area (China)

#### **1. Introduction**

Landslides are one of the most common types of natural hazards that cause serious economic losses, casualties, and damages to buildings, critical infrastructures and industrial settlements [1]. The impact of landslides is particularly high in rural mountainous areas, where land management is difficult to achieve, and risk mitigation is scarce due to the large distances, difficult logistics and harsh climatic conditions. As an example, on 2 May 2014, a large-scale landslide occurred in the Badakhshan Province in Northeastern Afghanistan due to continuous heavy rainfall. It resulted in nearly 2700 deaths and more than 300 houses were buried [2]. Again, a giant landslide occurred in

the Himachal Pradesh region of Northern India on 13 August 2017 that killed 47 people and nearly a 300-m highway was completely destroyed.

Monitoring and early warning systems are effective methods to reduce the risk of landslides [3], but the application of them on each single, potentially unstable slope is often impossible in mountain areas, due to the high spatial frequency of affected slopes. In such cases, remote sensing may help to reduce the cost and time for the application of mitigation measures, because it can provide a preliminary low-cost assessment of the severity of the slope instability and allow for the prioritization of critical cases. The use of remote-sensing data is one of the hotspots of landslide research at present. Among all the physical manifestations of a mass movement, the surface deformation is the most intuitive and comprehensive to measure and use for hazard assessment. It is a critical indicator for developing landslide early warning systems. Therefore, the implementation of deformation monitoring is of great significance for landslide prevention and mitigation [4,5]. The majority of fatal landslide events occur in less developed countries, such as China, India, Nepal, Bangladesh, etc. [1]. Traditional deformation monitoring equipment, such as an inclinometer, global positioning system (GPS), etc., are well-suitable, but their high costs limit their applications in underdeveloped areas.

Synthetic Aperture Radar Interferometry (InSAR) is an effective surface deformation monitoring method with a wide range and high accuracy. Commercial radar imaging can be expensive and mainly used for landslide hazard prevention in developed areas, including landslide identification [6,7], monitoring and early warning [8–10] and risk assessment [11]. With the free distribution of Sentinel-1 satellite images, InSAR has become a low-cost landslide monitoring technique. Small Baselines Subset (SBAS) InSAR is a time series analysis method of InSAR and can achieve effective monitoring in rural and urban areas. It has been gradually applied in landslide-prone areas, which provides new technical means for landslide monitoring and early warning [12–16].

The Three Gorges Reservoir area (TGRA) has been impounded since 2003. The periodical scheduling of the reservoir water level from 145 m to 175 m has been carried out since 2009. Affected by heavy rainfall and reservoir water level scheduling, a large number of new landslides and reactivation of old landslides have occurred [17]. Until now, more than 500 landslides have undergone varying degrees of deformation. However, due to the high cost of dedicated landslide monitoring, the latter was carried out only in 254 landslides in the TGRA [18], which cannot meet the needs of disaster prevention and mitigation.

The Muyubao landslide located in the TGRA was taken as a case study in this paper. Based on the Sentinel-1 radar images, the SBAS InSAR analysis method was applied to extract the information of landslide movement. The characterization of its temporal and spatial deformation was conducted as well. Moreover, the evolution mechanism of the Muyubao landslide was analyzed with the combined consideration of its engineering geological conditions and triggering factors, including rainfall and fluctuation of reservoir water level. This study aims to explore the feasibility and effectiveness of Sentinel-1 images and InSAR technique in landslide monitoring and analysis in the TGRA as a case example from an underdeveloped landslide-prone region.

#### **2. Case Study: Muyubao Landslide**

#### *2.1. Geological Conditions*

The Muyubao landslide occurred in Shazhenxi Town, Zigui County of Hubei Province on the bank of the Yangtze River and 56 km away from the Three Gorges Dam. The Muyubao landslide is delimited laterally by N-S-oriented trenches, while the upper portion is straight and smooth. It is of a chair-like shape, with smooth topography in the middle and lower parts and steeper terrain in the upper part (Figure 1). The elevation of the Muyubao landslide ranges from 120 m.a.s.l to 425 m.a.s.l, with an average slope angle of about 20◦ . It has an average thickness of 50 m and an estimated volume of 9 Hm3 . The entire landslide involves an area of 180 × 104 m2 , with a maximum longitudinal dimension of 1200 m and an average width of 1500 m (Figure 2).

**Figure 1.** (**a**) The location of the Three Gorges Reservoir area (TGRA), (**b**) location of the Muyubao landslide, and (**c**) geomorphological delimitation of the Muyubao landslide.

**Figure 2.** Topographical map of the Muyubao landslide.

The Muyubao landslide developed within a flysch formation dipping in the same direction and roughly parallel to the slope. The stratigraphy is made by the typical succession of quartz sandstone and sand-mudstone layers of the Jurassic Xiangxi formation (Figure 3). The sliding body is mainly composed of two parts. The surface layer is made up by loose deposits composed of alluvial sub-clay, gravel layers with clay, colluvial and eluvial sub-clay and collapsed block stones. The lower part of the sliding body is composed of highly disrupted layered Quartz sandstone, which is relatively intact in the western upper portion and weathered and fractured in the upper and lower portions.  Schematic geological cross-section I-I′ of the Muyubao landslide. Please see the location of I-I′ in Figure 2.

 The upper portion of the sliding body has a linear detachment surface in the profile map, with an inclination of about 25◦ and a thickness of 60~90 m. The front part is an uplift platform formed by shear sliding, with the thickness of 80~120 m, and the inclination angle of the rock layer is about 27◦ (Figure 3). The thickness values here refer to the central part of the landslide, as illustrated in the profile of Figure 3, which may decrease towards the boundaries. The landslide is still active, as it has experienced several reactivation events in the recent past [19,20]. This makes it a priority case for monitoring, since, in case the movement would develop towards sudden failure, it will endanger the lives and property of 140 households (500 people) in the landslide area and threaten the safety of the road and shipping of the Yangtze River.

#### *2.2. Field Investigation of Landslide Deformation*

‐ After the reservoir impounding to 175 m in the TGRA in September 2009, the deformation of the Muyubao landslide began to accelerate. Macro-deformations appeared in the landslide mass. In 2011, the road surface across both the boundaries of the landslide was damaged, and the roadbed on the west boundary was broken in widths of 20 cm (Figure 4a). The cracks in the wall of the village buildings on the eastern side of the landslide continued to be stretched (Figure 4b). In the investigation of 2014, a local collapse was found on the west side of the landslide (Figure 4c). The road across the boundaries was severely damaged, showing tensile cracks (Figure 4d). There was also a collapse under the road in the middle of the west side of Muyubao landslide, with a length of 50~70 m, a width of 10~15 m and a thickness of 5~10 m (Figure 4e). A tensile trench was formed by a partial collapse in the east side of the landslide, which is about 60 m long and 15 m wide (Figure 4f).

**Figure 4.** Deformations on the Muyubao landslide: (**a**) road damage on the west boundary, (**b**) deformation of buildings, (**c**) a local collapse on the west side, (**d**) road damage on the east side, (**e**) a collapse in the middle of the west side and (**f**) a tensile trench on the east side.

#### **3. Methods**

#### *3.1. Small Baselines Subset InSAR Analysis*

‐ ‐ With the development of the InSAR technique and the increasing of synthetic aperture radar (SAR) images, Ferretti et al. [21] proposed the method called Persistent Scatterer InSAR (PSInSAR) based on the differential InSAR technique. The main idea of the PSInSAR is to use multiple SAR images covering the same area to analyze the stability of the amplitude and phase and, then, to identify the pixels that are less-affected by the spatiotemporal decorrelation. To obtain accurate deformation information, the components of the phase need to be jointly analyzed and modeled to remove errors. In the PSInSAR method, a unique master image is selected from all images to generate interferograms. The strategy will result in the presence of long-space baseline interference pairs. It produces a low

density of the target reflectors in case of nonurban areas. Hence, in rural areas with few artificial buildings, the permanent scattering pixels may be very sparse [22].

The small baseline subset method was initially proposed to overcome the problem of decorrelation by making full use of the interferograms with both small temporal baselines and short perpendicular baselines [23,24]. In rural areas, plenty of pixels have no dominated scattering characteristics in the SAR image. The pixels showing a slow decorrelating filtered phase (SDFP) are widely distributed in natural terrain and maintained good coherence during a short time interval. The Stanford Method for Persistent Scatterer (StaMPS) was proposed in 2004 [25]. In the StaMPS SABS method, these kinds of pixels are identified and analyzed to monitor surface displacement. The SDFP pixels are selected through their phase characteristics. In order to reduce the calculation burden, an initial subset of pixels containing almost all SDFP pixels are firstly selected through amplitude analysis. The amplitude dispersion value is calculated as the indicator of phase stability. The wrapped phase is composed of spatially correlated phase and spatially uncorrelated look angle error. The spatially correlated phase, including ground deformation, elevation error and orbit error, is estimated by the bandpass filtering in the frequency domain. The spatially uncorrelated look angle error mainly consists of the spatially uncorrelated elevation error. It is estimated through its correlation with the perpendicular baseline. The residual, obtained by removing these two terms from the wrapped phase, gives an estimation of the decorrelation noise γ*x* (Equation (1)), which indicates the stability of the pixel phase:

$$\gamma\_{\mathbf{x}} = \frac{1}{N} \left| \sum\_{1}^{N} \exp \left\{ \sqrt{-1} (\psi\_{\mathbf{x},i} - \widetilde{\psi}\_{\mathbf{x},i} - \Delta \widehat{\phi}\_{\boldsymbol{\theta},\mathbf{x},i}^{\boldsymbol{u}}) \right\} \right| \tag{1}$$

where *N* is the number of interferograms, ψ*x*,*i* is the wrapped phase, φˆ *u* θ,*x*,*i* is the spatially uncorrelated look angle error and ψe*x*,*i* is the spatially correlated term. The final SDFP pixels are selected through the threshold analysis of γ*x*. Then, these SDFP pixels can be decomposed with the three-dimensional phase unwrapping method [26]. The unwrapped phase on a given pixel can be expressed as:

$$\psi\_{\mathbf{x},i} = \mathcal{W} \{ \phi\_{\mathbf{D},\mathbf{x},i} + \phi\_{\mathbf{A},\mathbf{x},i} + \phi\_{\mathbf{O},\mathbf{x},i} + \phi\_{\mathbf{T},\mathbf{x},i} + \phi\_{\mathbf{N},\mathbf{x},i} \} \tag{2}$$

where ψ*x*,*i* is the unwrapped phase, and φ*D*,*x*,*i* , φ*T*,*x*,*i* , φ*A*,*x*,*i* , φ*O*,*x*,*i* and φ*N*,*x*,*i* are the phase components due, respectively, to ground deformation, topographic error, atmospheric disturbance, inaccurate orbit information and noise. In the StaMPS SBAS analysis, a theoretical framework for three-dimensional phase unwrapping was proposed. These different phase components were estimated through iterative filtering with consideration of their characteristics in spatial-temporal domains. The deformation phase can be obtained by removing the other components from the unwrapped phase. Afterwards, a time series of deformation can be obtained from the phase using singular value decomposition.

#### *3.2. Time Series InSAR Processing of Muyubao Landslide*

Sentinel-1 is a two-satellite constellation with the prime objective of land. It is the first satellite developed by the European Commission and the European Space Agency (ESA) for the Copernicus Programme. This satellite was launched in April 2014. After half a year of trial operation, it began to be gradually used after October 2014. One goal of the mission is to provide C-Band SAR data continuity following the retirement of ERS-2 and the end of the Envisat mission. In this study, 37 images acquired by Sentinel-1a, in the period from 10 March 2016 to 13 September 2017 (Figure 5), were collected from the ESA [27]. We collected the precise orbit parameters of the Sentinel-1 in the ESA as well [28]. The angle of the orbit with the north is 12.16◦ , and the average incidence angle of three sub-swaths are 34.01◦ , 39.39◦ and 40.02◦ , respectively.

‐ **Figure 5.** The time distribution of sentinel-1 images. ‐

‐ ‐ ‐ Interferograms are generated applying the two-track differential method [29]. The SRTM (Shuttle Radar TopographyMission) DEM (digital elevation model) covering the study area is utilized to generate single-look differential interferograms. Temporal and perpendicular thresholds are, respectively, set to 90 days and 1000 m for generating interferograms. Considering that these thresholds will generate many interferograms and result in a heavy computational burden, each image was mostly combined to three subsequent images for interferometric processing. Data processing was carried out in COMET-LiCSAR [30]. As a result, 97 interferograms were produced as the input for the SBAS analysis (Figure 6). The amplitude difference dispersion was applied to selected coherent targets, and it was set as 0.6. Subsequently, the time series displacement was obtained by applying the StaMPS SBAS method. Moreover, in order to reduce the effects of residual atmospheric artefacts, the obtained time series displacement within the Muyubao landslide has been referenced to a stable neighbor point (31◦1 ′52.62"N, 110◦29′57.48"E) [15,31]. ‐ ‐ ′ ″ ′ ″

**Figure 6.** Baseline network of the interferograms.

‐ ‐ ‐ ‐ The direction of displacement acquired by the SBAS method is the line-of-sight (LOS) of radar satellite. Before applying the measurement to the landslide analysis, the LOS displacements should be projected along real slope components to get actual vectors. Hilley et al. [32,33] proposed a widely used projection method; the LOS displacement can be projected into the direction of the steepest descent slope. In this case study, a field investigation showed that the slip direction was 16◦ . Hence, the LOS displacement was projected onto the downward direction of along the sliding direction.

#### **4. Results**

#### *4.1. InSAR Results of Muyubao Landslide*

As stated in Section 3.2, the displacement time series of the Muyubao landslide was extracted by applying the time series InSAR technique and projected onto the direction of the downward sliding direction. Figure 7 shows the average deformation speed of the Muyubao landslide. The strong deformation mainly occurred at the eastern upper portion of the landslide, with deformation rates up to 100 mm/year. In the direction perpendicular to the sliding direction, the deformation of the eastern side was slightly stronger than the western side of the Muyubao landslide.

**Figure 7.** The spatial distribution of the Interferometric Synthetic Aperture Radar (InSAR) monitoring points, with deformation speed and the spatial distribution of the feature points. Note: Global Positioning System is abbreviated to GPS, and Persistent Scatterer InSAR is abbreviated to PSI.

‐ ‐ ‐ ‐ ‐ ‐ ‐ ‐ Some points in different parts of the Muyubao landslide were selected to further analyze the spatial deformation characteristics. The locations and cumulative displacement time series curves are shown in Figures 7 and 8, respectively. It can be seen that the cumulative displacements of PSI-04 and PSI-05 show the active part of the landslide (Figure 7). The cumulative displacement of these two points from March 2016 to September 2017 exceeded 250 mm (Figure 8d,e). The deformations of the monitoring points in the middle part were smaller, with the cumulative displacement about 150 mm (Figure 8f,h). The deformation rates were the smallest (about 50 mm/y) at the western front part of the landslide (Figures 7 and 8i–k). Compared to the linear monitoring curves at the upper portion of the landslide, the monitoring curves near the lower portion showed slight step-like characteristics (PSI-01, PSI-02, PSI-06, PSI-09 and PSI-10), which were directly affected by periodic rainfall and reservoir scheduling [34,35]. In summary, the deformation of the Muyubao landslide seems to show obvious

differences in space. The deformation on the eastern upper portion of the landslide was the strongest, followed by the middle part, and the western part showed the smallest deformation.

**Figure 8.** Time series displacement of InSAR and triggering factors.

#### *4.2. The Deformation Analysis of Muyubao Landslide*

The Three Gorges Reservoir oscillates regularly between 145 m and 175 m, with annual cycles (Figure 9). It can be seen from Figure 9 that, from January to February 2017, the TGRA remained at a normal water storage level of 175 m. Under the uplift pressure of the reservoir water, the Muyubao landslide was in a state of continuous deformation, the southeast edge of which showed a larger displacement rate (Figures 10f and 11f). From March to April 2017, with the slow decline of the reservoir level, the landslide gradually stabilized as a whole, and the displacement rates of the monitoring points were less than 10 mm/month (Figures 10g and 11g). From May to June 2017, the TGRA was experiencing its rainy season, while the reservoir level began to decline rapidly. The displacement rates of all monitoring points of the landslide increased gradually, some of which reached 30 mm/month (Figure 11h). Due to the hysteresis characteristics of rainfall infiltration, the displacement rates of the upper potion were greater than the middle part, which grew gradually with the increase of rainfall infiltration capacity and infiltration depth (Figure 10h). The landslide showed slow deformation when

the reservoir water kept at the level of 145 m from July to August 2017. The infiltration lag effect of the previous rainfall, combined with the large continuous rainfall in the reservoir area, triggered the strong deformation period of the year 2017. Some of the monitoring displacement rates reached 40 mm/month (Figure 11c,i).

The tensile cracks at the upper portion of the landslide provides a passage for the infiltration of continuous heavy rainfall. Therefore, the monitoring points of the eastern upper portion showed the largest cumulative displacements and displacement rates, which gradually decreased towards the front part (Figure 10i).

**Figure 9.** The operation curve of the Three Gorges Reservoir area (2010–2018).

**Figure 10.** Deformation velocity of the Muyubao landslide (time-yymmdd), the subfigures of ( ‐ **a**–**i**) are refer to the different time range of InSAR monitoring.

‐ **Figure 11.** Deformation velocity value distribution of the Muyubao landslide (time-yymmdd), the subfigures of (**a**–**i**) are refer to the different time range of InSAR monitoring.

#### **5. Discussion**

#### *5.1. The Reliability Analysis of InSAR Monitoring in the Muyubao Landslide*

Interferometry techniques have a huge potential to monitor landslides, as demonstrated by worldwide examples [12,15,32]. The monitoring accuracy of InSAR is affected by various factors, such as the atmosphere and terrain. Therefore, it is necessary to evaluate the reliability of the time series displacement before it is applied for landslide analysis. In the InSAR monitoring of the Muyubao landslide, the reliability analysis of the results can be carried out from two aspects.

‐ ‐ ‐ ‐ (a) Comparison of InSAR monitoring results at different locations. The monitoring point PSI-12 outside the Muyubao landslide was selected for comparative analysis. According to the InSAR monitoring results (Figure 12), the monitoring location of PSI-12 did not deform during the monitoring period, which is consistent with the results of the field investigation. While the monitoring points (PSI-01~PSI-11) within the landslide were deforming during the same period (Figure 8). By comparing the InSAR results of the monitoring points within and outside the landslide, it demonstrates that the InSAR monitoring in this study can identify the deformation and nondeformation zones.

‐ **Figure 12.** The InSAR monitoring results of PSI-12. ‐

‐ ‐ ‐ ‐ ‐ ‐ (b) Comparison of monitoring results with different techniques. GPS is a reliable monitoring technique of ground deformation, which has been widely used in various fields [36,37]. In order to verify the accuracy of InSAR monitoring, the comparison was carried out between the monitoring results of GPS and PSI points at similar locations of the Muyubao landslide (the locations are shown in Figure 7). Considering that the automatic GPS monitoring equipment was installed in June 2016, the monitoring results from June 2016 to September 2017 were used for comparison. The results of GPS-01-PSI-08 and GPS-02-PSI-07 are shown in Figures 13 and 14, respectively. The cumulative displacement of InSAR and GPS at similar locations are highly consistent in magnitude. It indicates that the InSAR monitoring of this study is reliable and can be utilized for landslide analysis.  Results comparison of GPS-01 and PSI-08.

 Results comparison of GPS-02 and PSI-07. ‐ ‐

#### *5.2. The Formation Mechanism of the Muyubao Landslide*

Based on the analysis of geological conditions and monitoring data (Figure 8), the deformation of the Muyubao landslide showed a large spatial difference. Large deformation mainly occurred in the middle and upper portions, while the deformation of the uplift at the front part was smaller (Figures 7 and 10). The profile with the largest deformation is shown in Figure 15. It can be seen that the deformation of the entire section is relatively small in the period of July-September 2016. It became larger in the period of January–March 2017. At the front edge of the profile near the Yangtze River, the Muyubao landslide deformed slightly in both periods. A field investigation was conducted for further analysis of the landslide deformation mechanism.

By analyzing the results of the field investigation, it can be seen that the deformation and failure mode is controlled by the lithology condition and the structure of the slope. Under the long-term gravitational process, the slope was creeping down along the weak mudstone layer. The front stress accumulated gradually, which resulted in the bending and uplifting in the middle and front stratum of the slope. Once the potential slip surface of the bending part connects, the accumulated energy may suddenly release, which may cause collapse and high-speed sliding. In conclusion, the evolution process can be divided into four stages: slipping stage, bending and uplifting stage, local shear stage and completely connecting stage of the shearing surface (Figure 16a) [19].

The bending and uplifting in the front stratum limited the deformation of the front part, which occurred when the slope was creeping down along the lines of weakness. Due to the topography and geomorphology conditions, the uplift in the front increases the sliding resistance effect of the landslide. When rainfall occurs, the gullies distributed on the surface of the landslide provides favorable channels for the rapid discharge of surface water, thus reducing the threat of water seepage to the landslide stability. Meanwhile, in the front part, the sliding surface is of a gentle slope, and the rock mass is bent and dips inside, which forms a sliding resistance for the landslide as well. In short, under natural conditions, the geological characteristics are beneficial to the overall stability of the landslide. However, due to the combined influence of factors in the TGRA, the Muyubao landslide has been undergoing a continuous creep deformation.

 Spatiotemporal characteristics of the deformation of the Muyubao landslide: (**a**) the velocity map of Muyubao landslide, and (**b**) the velocity of geological cross-section I-I′ .

**Figure 16.** (**a**) Conceptual model of the landslide evolution process (revised from Deng et al. [19]), and (**b**) a photo of the front part of the Muyubao landslide.

‐

‐

‐

‐

#### *5.3. The Relationship between Landslide and Influencing Factors*