Landslide Susceptibility Mapping for the Muchuan County (China): A Comparison Between Bivariate Statistical Models (WoE, EBF, and IoE) and Their Ensembles with Logistic Regression

: The main purpose of this study is to apply three bivariate statistical models, namely weight of evidence (WoE), evidence belief function (EBF) and index of entropy (IoE), and their ensembles with logistic regression (LR) for landslide susceptibility mapping in Muchuan County, China. First, a landslide inventory map contained 279 landslides was obtained through the ﬁeld investigation and interpretation of aerial photographs. Next, the landslides were randomly divided into two parts for training and validation with the ratio of 70 / 30. In addition, according to the regional geological environment characteristics, twelve landslide conditioning factors were selected, including altitude, plan curvature, proﬁle curvature, slope angle, distance to roads, distance to rivers, topographic wetness index (TWI), normalized di ﬀ erent vegetation index (NDVI), land use, soil, and lithology. Subsequently, the landslide susceptibility mapping was carried out by the above models. Eventually, the accuracy of this research was validated by the area under the receiver operating characteristic (ROC) curve and the results indicated that the landslide susceptibility map produced by EBF-LR model has the highest accuracy (0.826), followed by IoE-LR model (0.825), WoE-LR model (0.792), EBF model (0.791), IoE model (0.778), and WoE model (0.753). The results of this study can provide references of landslide prevention and land use planning for local government.


Introduction
As one of the most frequently-occurring geological disasters in the world, landslides have triggered a series of threats to human society including casualties, economic losses, infrastructure destruction, and geological environment problems [1][2][3]. Therefore, to reduce the losses, it is absolutely necessary to study the landslide susceptibility in a region [4,5]. According to the previous researches, landslide susceptibility can be roughly defined as the landslide occurrence probability in an area under the synergistic effect of a number of regional geological environmental factors [6,7]. Due to the large number and variability of landslide conditioning factors involved in the process, it is difficult to predict landslide-prone areas. What is known is that many geological and topographic conditions can influence the occurrence of landslides, such as formation lithology, faults, hydrology, attitude, slope angle, soil, and vegetation [8,9]. As landslide susceptibility mapping is the first and foremost step in landslide prevention, numerous researchers have been devoted to landslide susceptibility mapping in past years [10][11][12][13][14]. In general, the methods used in previous studies can be roughly divided into two types: qualitative and quantitative, for example, analytic hierarchy process (AHP) is the most commonly-used qualitative approach in landslide susceptibility mapping [14][15][16]. In recent

Landslide Conditioning Factors
According to the regional geological environment characteristics and previous studies [46][47][48], 12 landslide conditioning factors were selected, including altitude, plan curvature, profile curvature, slope angle, distance to roads, distance to rivers, topographic wetness index (TWI), normalized different vegetation index (NDVI), land use, soil, and lithology. A digital elevation model (DEM) with 30 × 30 m resolution provided by the Geospatial Data Cloud of Chinese Academy of Sciences (GSCloud) was introduced to generate a series of topographic factors, such as altitude, plan curvature, profile curvature, slope angle, slope aspect, and TWI [49]. The Landsite-8 image (LC81300402018099LGN00, April 9, 2018) with resolution of 30×30 m was obtained from the same place as DEM and was used to extract the NDVI map using ArcGIS software [50]. The land use map was extracted from the regional land use map at a scale of 1:100,000 provided by local Land and Resources Bureau. The distance to rivers and distance to roads were obtained from the topographic map with the scale of 1:50,000 that provided by the same place as land use map. The soil map was extracted from the regional soil map at a scale of 1:1,000,000 that provided by Institute of Soil Science, Chinese Academy of Sciences (ISSCAS) [51]. The lithology map was extracted from the regional lithology map with the scale of 1:200,000 that provided by National Geological Archives of China (NGAC) [52].. Moreover, satellite images from the Google Earth pro 7.1 were also applied to assist the research. All the data mentioned above were processed and used to generate the landslide conditioning factors in ArcGIS software.

Landslide Conditioning Factors
According to the regional geological environment characteristics and previous studies [46][47][48], 12 landslide conditioning factors were selected, including altitude, plan curvature, profile curvature, slope angle, distance to roads, distance to rivers, topographic wetness index (TWI), normalized different vegetation index (NDVI), land use, soil, and lithology. A digital elevation model (DEM) with 30 × 30 m resolution provided by the Geospatial Data Cloud of Chinese Academy of Sciences (GSCloud) was introduced to generate a series of topographic factors, such as altitude, plan curvature, profile curvature, slope angle, slope aspect, and TWI [49]. The Landsite-8 image (LC81300402018099LGN00, 9 April 2018) with resolution of 30 × 30 m was obtained from the same place as DEM and was used to extract the NDVI map using ArcGIS software [50]. The land use map was extracted from the regional land use map at a scale of 1:100,000 provided by local Land and Resources Bureau. The distance to rivers and distance to roads were obtained from the topographic map with the scale of 1:50,000 that provided by the same place as land use map. The soil map was extracted from the regional soil map at a scale of 1:1,000,000 that provided by Institute of Soil Science, Chinese Academy of Sciences (ISSCAS) [51]. The lithology map was extracted from the regional lithology map with the scale of 1:200,000 that provided by National Geological Archives of China (NGAC) [52]. Moreover, satellite images from the Google Earth pro 7.1 were also applied to assist the research. All the data mentioned above were processed and used to generate the landslide conditioning factors in ArcGIS software. As mentioned above, DEM was used to extract topographic related landslide conditioning factors, such as altitude, plan curvature, profile curvature, slope angle, slope aspect, and TWI. Altitude is one of the most frequently used factors in landslide susceptibility mapping [53][54][55]. The main reason is that the altitude can influence the topographic attributes which lead to spatial variability of different landscape processes. In addition, the altitude can also influence the vegetation distribution which has a close relationship with landslide. In this paper, the altitude was divided into eight categories, i.e. 290-500, 500-700, 700-900, 900-1100, 1100-1300, 1300-1500, 1500-1700, and 1700-1866 m (Figure 2a).
Plan curvature is considered as an important factor which controls the aggregation and separation of topography and has a direct influence on the velocity of water flow and thus erosion [41,56,57]. In this paper, plan curvature ranged from -26.54 to 26.21 and was reclassified into five categories, i.e. -26.54 to -3.17, -3.17 to -1.10, -1.10 to 0.56, 0.56 to 2.62, and 2.62 to 26.21 (Figure 2b).
Profile curvature is the curvature in vertical plane parallel of the slope section which also has an influence on the velocity of water flow and thus erosion [58][59][60]. In this study, the profile curvature ranged from -34. 72 Figure 2c).
Slope angle directly affects the stability of slope and it is considered as one of the most important factors in landslide susceptibility mapping [61][62][63]. In this study, the slope angle ranged from 0 to 77.14 • and was reclassified into eight categories, i.e. 0 • -10 Slope aspect affects the slope stability indirectly via hydrologic processes, such as the direction of rainfall and sunshine [64][65][66]. In addition, the slope aspect can also influence the weathering process and vegetation distribution [67,68]. In this study, the slope was classified into nine categories, i.e. flat(-1), north(337 Distance to roads is considered as an important factor in landslide susceptibility mapping as it is an anthropogenic factor which can change surface configuration, cause the slope to lose support and then induced the occurrence of landslide [69][70][71][72]. In this study, five different buffer zones were generated by the interval of 500 m in ArcGIS software, i.e. 0-500, 500-1000, 1000-1500, 1500-2000 and >2000 m (Figure 2f).
Distance to rivers affect the hydrologic processes of the slope. The closer to the river, the more the slope is eroded by the river. It can negatively affect the stability of slope via eroding the toe of slope [73][74][75][76]. In this study, five different buffer zones were classified by the interval of 200 m in ArcGIS software, i.e. 0-200, 200-400, 400-600, 600-800, and >800 m ( Figure 2g).
TWI is an important conditioning factor in landslide occurrence. It is a factor of soil moisture that has a profound influence on the most of landslides [77][78][79]. In this study, the TWI ranged from 0. 35 Figure 2h).
NDVI is a factor that indicates the vegetation growth status and vegetation coverage. The negative value shows the ground is cover by water, snow and so on. 0 means the ground is cover by rock or bare soil. A positive value indicates the ground is cover by vegetation, which increases with the increase of coverage. The higher the NDVI value, the lower the possibility of landslide occurrence [80][81][82]. In this study, the NDVI ranged from -0.13 to 0.58 and was reclassified into five categories, i.e. -0.13 to 0.13, 0.13 to 0.24, 0.24 to 0.30, 0.30 to 0.37, and 0.37 to 0.58 (Figure 2i).
The other important landslide-related factor in landslide susceptibility mapping is land use and it has been widely applied in previous researches [83][84][85]. In this study, the land use map was obtained through the interpretation of Landsat 8 image. Finally, the land use was reclassified into six categories, i.e. farm land, forest land, grass land, water, residential areas, and bare land (Figure 2j). Soil is the common component that constitutes of slope. Different soil has different physical and mechanical properties and it can affect the surface water infiltration and groundwater flow [86][87][88][89]. The study area has ten soil types (Figure 2k).
Lithology is also considered as an important landslide-related factor in landslide susceptibility mapping [90]. The variation of lithology may lead to the variation of strength and permeability of rock stratum. Usually, the landslide slides along a rock stratum with low strength and poor permeability [91][92][93][94]. The study area has ten lithology types (Figure 2l). Lithology is also considered as an important landslide-related factor in landslide susceptibility mapping [90]. The variation of lithology may lead to the variation of strength and permeability of rock stratum. Usually, the landslide slides along a rock stratum with low strength and poor permeability [91][92][93][94]. The study area has ten lithology types (Figure 2l).

Preparation of Training and Validation Datasets
Chung and Pham indicated that model validation process must be established on the basis of dividing the dataset [95,96]. Therefore, it was necessary to split the dataset into two parts in landslide susceptibility mapping. The first part was applied to build the models called training dataset and the rest was used to verify the model performances named validation dataset. However, there are no general rules for splitting the proportion of training and validation datasets [97]. According to the previous studies [98][99][100], the most commonly-used ratio (70/30) was introduced to select training and validation datasets, respectively. Finally, the training dataset was built by overlaying the 70% (195) of landslides onto 12 landslide conditioning factor layers. Oppositely, the rest of the landslides were used to construct the validation dataset. After the preparation of training and validation datasets, all the data analysis processes were carried out in ArcGIS software and the results were validated by Wilcoxon signed-rank test and ROC curve.

Weight of Evidence
Weight of evidence (WoE) is a spatial information integration model based on the Bayesian probability model, which can evaluate and predict objects with spatial significance and is one of the most important techniques in mineral potential assessment [39,101]. During the past years, the WoE model has been gradually applied in environmental evaluation, geological hazard and the survey of forest pest [102][103][104]. In this paper, the WoE model was introduced to landslide susceptibility research for its great suitability in analyzing the relationships between spatial distribution of landslides and conditioning factors.
The original description of mathematical formulation of WoE model was introduced by Bonham-Carter [105,106]. In WoE modeling, the weight of each conditioning factor class was calculated based on the occurrence of landslides within the area [107]. The weight in this paper can be separated into positive weight and negative weight and their calculation formulas are presented in Equation (1) and Equation (2), respectively.

Preparation of Training and Validation Datasets
Chung and Pham indicated that model validation process must be established on the basis of dividing the dataset [95,96]. Therefore, it was necessary to split the dataset into two parts in landslide susceptibility mapping. The first part was applied to build the models called training dataset and the rest was used to verify the model performances named validation dataset. However, there are no general rules for splitting the proportion of training and validation datasets [97]. According to the previous studies [98][99][100], the most commonly-used ratio (70/30) was introduced to select training and validation datasets, respectively. Finally, the training dataset was built by overlaying the 70% (195) of landslides onto 12 landslide conditioning factor layers. Oppositely, the rest of the landslides were used to construct the validation dataset. After the preparation of training and validation datasets, all the data analysis processes were carried out in ArcGIS software and the results were validated by Wilcoxon signed-rank test and ROC curve.

Weight of Evidence
Weight of evidence (WoE) is a spatial information integration model based on the Bayesian probability model, which can evaluate and predict objects with spatial significance and is one of the most important techniques in mineral potential assessment [39,101]. During the past years, the WoE model has been gradually applied in environmental evaluation, geological hazard and the survey of forest pest [102][103][104]. In this paper, the WoE model was introduced to landslide susceptibility research for its great suitability in analyzing the relationships between spatial distribution of landslides and conditioning factors.
The original description of mathematical formulation of WoE model was introduced by Bonham-Carter [105,106]. In WoE modeling, the weight of each conditioning factor class was calculated based on the occurrence of landslides within the area [107]. The weight in this paper can be separated into positive weight and negative weight and their calculation formulas are presented in Equation (1) and Equation (2), respectively.
where P is the probability. In represents the natural logarithm. B is the presence of potential landslide predictive factor, while B is the absence of potential landslide predictive factor. L is the presence of landslide. L is the absence of landslide. W + i means positive weight that shows a positive relationship between predictable variable and landslide and W − i indicates the correlation of them is negative. The difference between W + i and W − i is the weight contrast, which magnitude indicates all the spatial association between predictable variable and landslides. It's computational equation is presented in Equation (3).
where C i is the weight contrast of W + i and W − i .

Evidential Belief Function
The Dempster-Shafer theory (DST) of evidence proposed by Shafer in 1976 is regarded as a spatial integration model with mathematical representation. The framework of evidence belief function (EBF) model is developed from DST of evidence [108,109]. The EBF model usually applied as an effective approach to mineral potential mapping for its powerful ability of analyzing the incomplete data [110]. In addition, it also has the power to combine multiple sources of evidence. The model can evaluate how close the evidence proves rather than give the probabilities that the assumption is true [111]. EBF model is consisted of four basic functions, namely the degree of belief (Bel), the degree of disbelief (Dis), the degree of uncertainly (Unc) and the degree of plausibility (Pls). In recent years, EBF model has been widely and successfully applied in landslide susceptibility mapping [112][113][114]. In this paper, the basic functions of EBF model were determined by overlying the landslide inventory map on each landslide conditioning factor layer and the Bel was regard as the symbol of relationship between landslide and landslide conditioning factor. The equations of the above four functions can be expressed in Equations (4)- (8).
Bel + Unc + Dis = 1 (8) where Bel i is the degree of belief of ith conditioning factor. Dis i is the degree of disbelief of ith factor. Similarly, Unc i means the degree of uncertainty of ith factor. Pls shows the upper limits of probability. Oppositely, Dis shows the lower limits of probability.

Index of Entropy
The third model used in this paper is the index of entropy (IoE) model that is based on the principle of bivariate analysis [115,116]. It can calculate the weight of each input variable and the weight indicates the disorder of which parameter is the most relevant for the occurrence of landslide in a natural environment. The equations applied to calculate the weight of conditioning factors are expressed in Equations (9)- (14).
where a is the domain percentage. b is the landslide percentage. P ij is the probability density of class i of factor j. H jmax and H j are the entropy values of factor j. S j is the amount of classes factor j. I j is the information coefficient factor j. W j is the weight for the parameter as a whole.

Logistic Regression
Logistic regression (LR) model believes several physical parameters may affect the probability of landslide occurrence [117][118][119]. The LR model allows developing a multivariate regression correlation between a dependent variable and several independent variables [120,121]. The advantage of LR model is that through the addition of suitable link function to the usual linear regression model, the variables can be continuous, discontinuous, or both, and they do not to obey normal distribution [122,123]. In the present study, the main purpose of the LR model is to find the most suitable approach to obtain the relationship of the presence or absence of landslides with a series of independent variables. The relationship between the landslide occurrence and independent variables (landslide conditioning factors) can be described in Equation (15).
where P is the probability of landslide occurrence and its range is 0 to 1. Z is a linear sum of constants that obtained through the product of the independent variables and their respective coefficients. The interval of Z is -∞ to +∞ and the computational equation is shown in Equation (16).

Multicollinearity Analysis
Before the data analysis, the most important step is to verify whether the landslide conditioning factors are correlated with each other and the size of the correlation between them. If there exists a strong correlation between two or more factors, it will make it difficult to predict landslide occurrence [124]. Therefore, the multicollinearity analysis was introduced to illustrate the connections between the factors [125]. There are two commonly used statistical parameters in multicollinearity analysis, namely tolerance (TOL) and variance inflation factor (VIF), and they are a pair of reciprocals. According to previous studies [6,126], it can be considered that variables are mutually independent when the range of TOL value is 0.1 to 1. The multicollinearity analysis results of landslide conditioning factors under different models are calculated by SPSS software (Table 2) [127]. The results indicate all the factors are independent from each other. The calculation formulas of TOL and VIF are shown in Equation (16) and Equation (17), respectively. In addition, the results of correlation between landslides and conditioning factors calculated by WoE, EBF and IoE models are shown in Table 3.
where R i is the negative correlation coefficient of the ith independent variable that makes regression analysis on other independent variables.

Generating Landslide Susceptibility Maps
Generating landslide susceptibility maps is the final output step after the processes of independence test of conditioning factors, establishing, training, and validation models, and it can be roughly divided into two steps [128]. The first step is acquiring the landslide susceptibility indexes (LSIs) of all evaluation units. The second is regrouping the LSIs. For the first step, the LSI of each evaluation unit was calculated by probability distribution functions of evaluation models. For the second step, the LSIs were reclassified into several intervals by standard classification methods in ArcGIS software. Table 4 shows the coefficient of landslide conditioning factor in each ensemble model. When the coefficient is greater than 0, it indicates that the factor will promote the occurrence of landslide. On the contrary, a coefficient less than 0 means the factor will impede the occurrence of landslide. While 0 indicates the factor has no effect on the occurrence of landslide. It can be seen that only slope aspect impeded the occurrence of landslide and profile curvature performed best in the promotion factors, especially in IoE-LR model.
For the slope angle, landslides were most concentrate on the class 50 • -60 • which owned the Bel of 0.337, indicating it has the highest probability of landslide occurrence, followed by the class 20 • -30 • . For the land use, the Bel was highest for farm land (0.421) and lowest for grass land, water and bare land (they have the same Bel of 0). For the soil, the Bel was highest for the type of A (0.188). In the case of attitude, the Bel was 0.387 for the class 500-700 m. In the case of lithology, the Bel was 0.245 for the type B, indicating landslides were more likely to occur in this lithological condition. For the TWI, the Bel was the highest for class 0. 35-1.63. For the NDVI, the Bel was 0.286 for class 0.30-0.37. For the slope aspect, the Bel was the highest for northeast slopes (0.152) and lowest for flat slopes (0.0). In the case of distance to roads and rivers, the Bel decreases when distance to a road or river line increases. In the case of plan and profile curvature, the most of landslides occurred in class 0.56-2.62 and -1.39-0.88, respectively.
The process of LSI calculated by EBF model is illustrated as an example. As described above, Bel decides the LSI of EBF model, and the calculation equation is shown in Equation (19). Similarly, the equations for generating LSI of WoE and IoE models are shown in Equation (20) and Equation (21), respectively. In addition, the EBF-LR model is also illustrated as an example of ensemble models. Generally, the process of LSI produced by ensemble model can be roughly described as multiplying the coefficients of ensemble model by the weights obtained in bivariate statistical model. For example, the LSI of EBF-LR model was obtained by multiplying the coefficients of EBF-LR model and Bels generated by EBF model, as shown in Equation (22). Similarly, the equations for calculating LSI of WoE-LR and IoE-LR models are shown in Equation (23) and Equation (24), respectively. Finally, the natural break method was applied to reclassify the indexes into five classes, such as very low, low, moderate, high and very high, and the landslide susceptibility maps generated by EBF, WoE, IoE, EBF-LR, WoE-LR, and IoE-LR models are shown in Figure 3.

Model Validation and Comparison
Model validation is an important part of landslide susceptibility research which determines the accuracy of the study [26,129]. The landslide susceptibility maps will have no practical significances without model validation [130]. In the first place, the Wilcoxon signed-rank test was introduced to test the independence of all evaluation models (table 5). The 95% confidence interval of Wilcoxon signed-rank test is (-∝, -1.96) ∪ (1.96, +∝). It can be seen form the table 5 that all the evaluation models are independent from each other. Then, the receiver operating characteristic (ROC) curve was applied to assess the predictive ability of those models. The ROC curve is a widely-used statistical approach as it has the ability to evaluate the prediction accuracy of models quantitatively. The area under the ROC curve (AUC) shows the performances of landslide susceptibility evaluation models in this paper. Figure 4 represents the ROC curves of evaluation models and table 6 shows the parameters of ROC curves. It can be seen that EBF-LR model has the highest accuracy (0.

Model Validation and Comparison
Model validation is an important part of landslide susceptibility research which determines the accuracy of the study [26,129]. The landslide susceptibility maps will have no practical significances without model validation [130]. In the first place, the Wilcoxon signed-rank test was introduced to test the independence of all evaluation models ( Table 5). The 95% confidence interval of Wilcoxon signed-rank test is (-∝, -1.96) ∪ (1.96, +∝). It can be seen form the Table 5 that all the evaluation models are independent from each other. Then, the receiver operating characteristic (ROC) curve was applied to assess the predictive ability of those models. The ROC curve is a widely-used statistical approach as it has the ability to evaluate the prediction accuracy of models quantitatively. The area under the ROC curve (AUC) shows the performances of landslide susceptibility evaluation models in this paper. Figure 4 represents the ROC curves of evaluation models and Table 6 shows the parameters of ROC curves. It can be seen that EBF-LR model has the highest accuracy (0.826), followed by IoE-LR model

Model Validation and Comparison
Model validation is an important part of landslide susceptibility research which determines the accuracy of the study [26,129]. The landslide susceptibility maps will have no practical significances without model validation [130]. In the first place, the Wilcoxon signed-rank test was introduced to test the independence of all evaluation models (table 5). The 95% confidence interval of Wilcoxon signed-rank test is (-∝, -1.96) ∪ (1.96, +∝). It can be seen form the table 5 that all the evaluation models are independent from each other. Then, the receiver operating characteristic (ROC) curve was applied to assess the predictive ability of those models. The ROC curve is a widely-used statistical approach as it has the ability to evaluate the prediction accuracy of models quantitatively. The area under the ROC curve (AUC) shows the performances of landslide susceptibility evaluation models in this paper. Figure 4 represents the ROC curves of evaluation models and table 6 shows the parameters of ROC curves. It can be seen that EBF-LR model has the highest accuracy (0.

Discussions
Landslide is a kind of extremely dangerous geological disasters which has a significant influence on personal safety as well as property security and geological environment [131,132]. In this paper, based on LR model, three bivariate statistical models, namely EBF, WoE, and IoE models, were introduced to carry out the landslide susceptibility mapping in Muchuan County. The performances of those models were compared and validated by some statistical methods, such as ROC curve. Although those four models and their ensembles have been applied to landslide susceptibility mapping many times in recent years, they have the advantages of high accuracy and convenient operation. Althuwaynee et al. proposed a Geographic Information Systems (GIS)-based EBF model for data analysis and conducting the landslide susceptibility research in Kuala Lumpur city, Malaysia. They indicated that the prediction accuracy of EBF model is acceptable and landslide susceptibility map produced by it is trustworthy [38]. Wang et al. applied WoE model to landslide susceptibility mapping at Daguan County, Yunnan Province, China. The results believed that the method not only provided a reliable and easily interpreted result, but also simplified the modeling process and greatly improved the work efficiency [114]. Devkote et al. carried out the landslide susceptibility maps by IoE, LR, and certainty factor (CF) models at Mugling-Narayanghat road section, Nepal. After the validation and comparison processes, IoE model gained the highest AUC value, followed by LR and CF models [133]. As for the ensemble models, Chen et al. conducted a series of related studies and they reached the same conclusion as in this paper that the ensemble model obtained higher accuracy than their single models [58,130,134]. In the present study, a total of 279 landslides were identified and mapped in the landslide inventory map, of which 195 (70%) landslides were randomly selected to establish the landslide models and the rest of the 84 (30%) landslides were used to validate the veracity of models. As there are no general rules for selecting landslide conditioning factors, twelve factors such as altitude, plan curvature, profile curvature, slope angle, distance to roads, distance to rivers, TWI, NDVI, land use, soil, and lithology were prepared for the research according to the previous studies and local geological environment characteristics.
After the selection of landslide conditioning factors, all the factors were subjected to multicollinearity analysis so as to avoid the influences brought by linear factors. It can be seen from Table 2 that all the factors are independent from each other. W + and Wproduced by WoE model are a pair of weight parameters that named positive weight and negative weight. A positive weight means it may promote the occurrence of landslide, while negative weight means the opposite. The parameter of Bel generated by EBF model can be regard as the symbol that indicates the relationship between landslide and landslide conditioning factor. The parameter of W j generated by IoE model represents the importance of each factor. It can be seen from the Table 3 that slope angle has the highest weight (0.229) which means it is the most important factor, followed by land use (0.207), soil (0.194), altitude (0.190), lithology (0.141), TWI (0.08), distance to roads (0.062), slope aspect (0.051), NDVI (0.037), distance to rivers (0.022), plan curvature (0.009), and profile curvature (0.004). In addition, it is obvious to find that Bel and W + show the same results.
It can be seen from Figure 3 that the south and west sides of all the landslide susceptibility maps have the low and very low landslide susceptibility. From the parameter layers, it can be concluded that both the south and west sides have four similar characteristics, such as high altitude, far from the roads, high ratio of vegetation coverage and land use type is forest land. This may be the presence of cliffs with rock which is hard to weathering in high-altitude areas [133,135,136]. Besides, far from the roads means it is less affected by human activities [137,138]. In addition, high ratio of vegetation coverage and forest land indicates the vegetation growth status in this area is good. In other words, the better the growth status of vegetation, the lower the possibility of landslide occurrence [139][140][141].
The ROC curve, SE, and 95% CI were introduced to validate and compare the performances of evaluation models. According to the parameters of ROC curves, EBF-LR model performed best with the AUC value of 0.826, followed by IoE-LR, WoE-LR, EBF, IoE, and WoE models with the AUC values of 0.825, 0.792, 0.791, 0.778, and 0.753, respectively. Generally, all the models performed well in this study. Based on the LR model, combined with EBF, WoE and IoE models, the landslide susceptibility of Muchuan County was researched successfully. Last but not least, the models used in this study are worthwhile to apply in other areas where have the same geological environment characteristics as study area and the results obtained in this study can provide references for local government.

Conclusions
In this case research, WoE, EBF, IoE, and LR models were introduced to generate the landslide susceptibility maps for the Muchuan County, China. A total of 279 landslides were identified and mapped through the field investigation and interpretation of aerial photographs. The models were established by randomly selected 70% of the landslides and validated by the remaining landslides. The correlation between landslide occurrence and 12 landslide conditioning factors such as altitude, plan curvature, profile curvature, slope angle, distance to roads, distance to rivers, TWI, NDVI, land use, soil, and lithology was evaluated by the above models. Finally, six landslide susceptibility maps were classified into five classes, such as very low, low, moderate, high, and very high, by natural break method. The differences of spatial prediction ability of different models are compared by AUC values and the validation results indicate the EBF-LR model has the highest AUC value (0.826), which means EBF-LR model performed best in this research, followed by IoE-LR, WoE-LR, EBF, IoE and WoE models with the AUC values of 0.825, 0.792, 0.791, 0.778, and 0.753, respectively. Finally, those landslide susceptibility maps can provide references of landslide prevention and land use planning for local government.
Author Contributions: R.L. and N.W. contributed equally to the work. R.L. collected field data, conducted the landslide mapping and wrote the manuscript. N.W. provided critical comments in planning this paper and edited the manuscript. Both the authors discussed the results and edited the manuscript.