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

Abstract

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 field 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, profile curvature, slope angle, distance to roads, distance to rivers, topographic wetness index (TWI), normalized different 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.

1. 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 years, machine learning method has been gradually applied in landslide susceptibility mapping researches, such as artificial neural network (ANN) [17,18,19], support vector machine (SVM) [20,21,22], logistic model tree (LMT) [23,24], rotation forest (RF) [25,26], classification and regression tree (CART) [27,28], adaptive neuro-fuzzy inference systems (ANFIS) [29,30], and genetic algorithm (GA) [31,32]. Furthermore, statistical approach is another widely-used model which can also be divided into two types: bivariate and multivariate. In statistical approaches, the weight of each class of every factor was calculated by overlying the landslide inventory map and landslide conditioning factors map [33,34]. The frequently-used statistical approaches are frequency ratio (FR) [35,36,37], logistic regression (LR), evidential belief function (EBF) [38,39], weights of evidence (WoE) [40,41], certainty factor (CF) [42,43], and information value (IV) [44,45].

In this paper, three bivariate statistical methods, namely weight of evidence (WoE), evidence belief function (EBF), index of entropy (IoE), coupled with logistic regression (LR) were introduced to carry out the landslide susceptibility mapping in Muchuan County, which lies in the southeast of Sichuan Province, China. All the approaches were first applied in Muchuan County and the main purpose of this article is to obtain the landslide susceptibility maps through the above models. The main difference between this study and literatures referenced above is to achieve accurate landslide susceptibility maps for Muchuan County through the combination and comparison of those models. At last, the results were validated by the area under the receiver operating characteristic (ROC) curve.

2. Geological and Geomorphological Setting

The study area (Muchuan County) is located in Leshan city, Sichuan Province, China. It lies between longitudes 103°21′ to 104°37′E, and latitudes 28°49′ to 29°31′N (Figure 1). It covers an area of 1408 km². The attitude of the study area ranges from 290 to 1866 m above the sea level. The annual rainfall is 1332 mm and it mainly concentrate from June to September. Geologically, the study area located in the Muchuan–Mabian arcuate fold belt of southwest of Yangtze para-platform where folds are developed more than faults. The study area is covered by various types of lithological units such as mudstone, sandstone, siltstone, and coal measure strata (

Table 1. Formations.

Group	Lithology	Geologic Ages
A	Brick red massive allochemical rock, sandstone sandwiched mudstone and siltstone Brick red thin-thick layer silty fine-grained arkose, mudstone Brick red massive allochemical rock, sandstone sandwiched mudstone and siltstone	Cretaceous
B	Grayish-purple arkose, siltstone sandwiched mud shale and coquina Bright red mudstone, sandwiched with the same color sandstone and siltstone Grayish-purple arkose, siltstone sandwiched mud shale and coquina	Jurassic
C	Yellow-gray feldspar-quartz sandstone interbedded with purple-red mudstone	Jurassic
D	Magenta mudstone, quartz sandstone, and siltstone sandwiched biosparite and marl	Jurassic
E	Gray sandstone, siltstone, and mudstone Grayish-yellow debris-feldspar, siltstone, mudstone, and coal	Triassic
F	Yellow-gray medium-thick dolomite sandwiched limestone, gypsum salt, and salt-soluble breccia	Triassic
G	Limestone, dolomite, and shale	Triassic
H	Yellow-green siltstone sandwiched with mudstone and coal	Permian
I	The upper part is limestone and dolomite and the lower part is shale sandwiched siltstone	Permian
J	Grayish-green dense amygdaloidal basalt sandwiched picrite, tuff sand mudstone, shed coal and siliceous rock	Permian

).

3. Materials and Methods

3.1. 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 to 29.66 and was reclassified into five categories, i.e. –34.72 to –4.17, –4.17 to –1.39, –1.39 to 0.88, 0.88 to 3.66, and 3.66 to 29.66 (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°, 10°–20°, 20°–30°, 30°–40°, 40°–50°, 50°–60°, 60°–70°, and 70°–77.14° (Figure 2d).

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.5°–360°, 0–22.5°), northeast(22.5°–67.5°), east(67.4°–112.5°), southeast(112.5°–157.5°), south(157.5°–202.5°), southwest(202.5°–247.5°), west(247.5°–292.5°), and northwest(292.5°–337.5°) (Figure 2e).

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 to 14.58 and was reclassified into five categories, i.e. 0.35 to 1.63, 1.63 to 2.41, 2.41 to 3.36, 3.36 to 4.81, and 4.81 to 14.58 (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).

3.2. 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.

3.3. 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.

$$W_i^{+}=\mathrm{In}\frac{p\left\{B\right|L\}}{P\left\{B\right|\overline{L}\}}$$

(1)

$$W_i^{-}=\mathrm{In}\frac{p\left\{\overline{B}\right|L\}}{P\left\{\overline{B}\right|\overline{L}\}}$$

(2)

where P is the probability. In represents the natural logarithm. B is the presence of potential landslide predictive factor, while $\overline{B}$ is the absence of potential landslide predictive factor. L is the presence of landslide. $\overline{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).

$$C_i=W_i^{+}-W_i^{-}$$

(3)

where C_i is the weight contrast of $W_i^{+}$ and $W_i^{-}$.

3.4. 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=\frac{Be{l}_1+Be{l}_2+\dots +Be{l}_n}{1-{\displaystyle \sum _{i=2}^{n}Be{l}_{i-1}Di{s}_i-Di{s}_{i-1}Be{l}_i}}$$

(4)

$$Dis=\frac{Di{s}_1+Di{s}_2+\dots +Di{s}_n}{1-{\displaystyle \sum _{i=2}^{n}Be{l}_{i-1}Di{s}_i-Di{s}_{i-1}Be{l}_i}}$$

(5)

$$Unc=\frac{{\displaystyle \sum _{i=2}^n(Un{c}_{i-1}Un{c}_i+Be{l}_{i-1}Un{c}_i+Be{l}_{i}Un{c}_{i-1}+Di{s}_{i-1}Un{c}_i+Di{s}_{i}Un{c}_{i-1})}}{1-{\displaystyle \sum _{i=2}^{n}Be{l}_{i-1}Di{s}_i-Di{s}_{i-1}Be{l}_i}}$$

(6)

$$Pls=Bel+Unc$$

(7)

$$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.

3.5. 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).

$$P_{ij}=\frac{b}{a}$$

(9)

$$(P_{ij})=\frac{P_{ij}}{{\displaystyle \sum _{j=1}^{S_j}{P}_{ij}}}$$

(10)

$$H_j=-{\displaystyle \sum _{i=1}^{S_j}(P_{ij})}{\log }_2(P_{ij})j=1,2,3,\dots ,n$$

(11)

$$H_{j\max }={\log }_2{S}_j$$

(12)

$$I_j=\frac{H_{j\max }-H_j}{H_{j\max }},I=(0,1)j=1,2,3,\dots ,n$$

(13)

$$W_j=I_j\times P_{ij}$$

(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.

3.6. 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).

$$P=\frac{1}{1+\exp (-Z)}$$

(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).

$$Z=\alpha +{\beta }_1{x}_1+{\beta }_2{x}_2+{\beta }_3{x}_3+\dots +{\beta }_n{x}_n$$

(16)

where α is a constant, β_i(i = 1,2,3,…,n) is the coefficient, and x_i(i = 1,2,3,…,n) is the independent variable.

4. Results and Analysis

4.1. 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. Analysis. WOE: weight of evidence; EBF: evidence belief function; IoE: index of entropy; TOL: tolerance; VIF: variance inflation factor.

Number	Factors	WOE		EBF		IoE
		TOL	VIF	TOL	VIF	TOL	VIF
1	Slope aspect	0.921	1.086	0.922	1.085	0.924	1.083
2	Altitude	0.835	1.197	0.659	1.517	0.659	1.517
3	Land use	0.861	1.162	0.830	1.206	0.828	1.208
4	Lithology	0.691	1.447	0.577	1.733	0.577	1.732
5	NDVI	0.980	1.020	0.953	1.050	0.954	1.048
6	Plan curvature	0.898	1.114	0.894	1.118	0.905	1.105
7	Profile curvature	0.912	1.097	0.896	1.116	0.929	1.076
8	Distance to rivers	0.981	1.020	0.969	1.032	0.973	1.028
9	Distance to roads	0.764	1.309	0.706	1.416	0.707	1.414
10	Slope angle	0.869	1.151	0.802	1.247	0.818	1.222
11	Soil	0.798	1.253	0.717	1.394	0.723	1.384
12	TWI	0.916	1.091	0.860	1.163	0.863	1.158

) [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. Between landslides and conditioning factors using WoE, EBF, and IoE models. Bel: belief; Dis: the degree of disbelief: Unc: the degree of uncertainly; Pls: the degree of plausibility.

Factors	Class	No. of Landslide	No. of Pixels in Domain	W+	W-	C	Bel	Dis	Unc	Pls	Wj
Altitude (m)	290–500	42	727325	0.035	−0.009	0.044	0.270	0.125	0.605	0.875	0.190
	500–700	110	1330115	0.394	−0.352	0.746	0.387	0.089	0.524	0.911
	700–900	38	640772	0.061	−0.014	0.076	0.277	0.125	0.598	0.875
	900–1,100	5	354758	−1.375	0.081	−1.456	0.066	0.137	0.797	0.863
	1,100–1,300	0	246634	0.000	0.073	0.000	0.000	0.000	1.000	1.000
	1,300–1,500	0	144601	0.000	0.042	0.000	0.000	0.000	1.000	1.000
	1,500–1,700	0	50150	0.000	0.014	0.000	0.000	0.000	1.000	1.000
	1,700–1,866	0	2344	0.000	0.001	0.000	0.000	0.000	1.000	1.000
Plan curvature	−26.54–−3.17	6	125417	−0.153	0.005	−0.159	0.165	0.199	0.636	0.801	0.009
	−3.17–−1.10	39	592578	0.166	−0.037	0.203	0.228	0.191	0.582	0.809
	−1.10–0.56	75	1660025	−0.211	0.158	−0.369	0.156	0.232	0.612	0.768
	0.56–2.62	62	909888	0.200	−0.081	0.282	0.236	0.182	0.582	0.818
	2.62–26.21	13	208791	0.110	−0.007	0.118	0.215	0.196	0.588	0.804
Profile curvature	−34.72–−4.17	9	170680	−0.056	0.003	−0.059	0.199	0.201	0.600	0.799	0.004
	−4.17–−1.39	43	777064	−0.008	0.002	−0.010	0.209	0.201	0.590	0.799
	−1.39–0.88	76	1299477	0.048	−0.029	0.077	0.221	0.195	0.585	0.805
	0.88–3.66	57	1012294	0.010	−0.004	0.014	0.212	0.200	0.588	0.800
	3.66–29.66	10	237184	−0.280	0.018	−0.297	0.159	0.204	0.637	0.796
Slope angle (°)	0–10	0	738575	0.000	0.237	0.000	0.000	0.000	1.000	1.000	0.229
	10-20	84	1233686	0.200	−0.128	0.328	0.168	0.110	0.722	0.890
	20-30	63	836797	0.300	−0.117	0.417	0.185	0.111	0.703	0.889
	30-40	29	433617	0.182	−0.029	0.210	0.165	0.122	0.714	0.878
	40-50	11	186715	0.055	−0.003	0.058	0.145	0.125	0.730	0.875
	50-60	8	58411	0.899	−0.025	0.924	0.337	0.122	0.541	0.878
	60-70	0	8690	0.000	0.002	0.000	0.000	0.000	1.000	1.000
	70–77.14	0	208	0.000	0.000	0.000	0.000	0.000	1.000	1.000
Slope aspect	Flat	0	1192	0.000	0.000	0.000	0.000	0.000	1.000	1.000	0.051
	North	28	508458	−0.013	0.002	−0.015	0.124	0.111	0.764	0.889
	Northeast	33	489376	0.190	−0.035	0.225	0.152	0.107	0.741	0.893
	East	26	507513	−0.085	0.014	−0.099	0.115	0.113	0.772	0.887
	Southeast	20	402282	−0.115	0.014	−0.129	0.112	0.113	0.775	0.887
	South	24	393289	0.090	−0.012	0.102	0.138	0.110	0.753	0.890
	Southwest	16	351360	−0.203	0.020	−0.223	0.103	0.113	0.784	0.887
	West	20	414052	−0.144	0.018	−0.162	0.109	0.113	0.778	0.887
	Northwest	28	429177	0.157	−0.024	0.181	0.147	0.108	0.744	0.892
Distance to roads (m)	0–500	69	827671	0.402	−0.167	0.569	0.271	0.167	0.562	0.833	0.062
	500–1,000	48	611007	0.343	−0.091	0.433	0.255	0.180	0.564	0.820
	1,000–1,500	36	493759	0.268	−0.052	0.320	0.237	0.188	0.575	0.812
	1,500–2,000	20	366762	−0.022	0.003	−0.025	0.177	0.198	0.625	0.802
	> 2,000	22	1197500	−1.110	0.300	−1.410	0.060	0.267	0.674	0.733
Distance to rivers (m)	0–200	57	1032619	−0.010	0.004	−0.014	0.208	0.201	0.591	0.799	0.022
	200–400	69	844083	0.382	−0.160	0.543	0.308	0.171	0.521	0.829
	400–600	29	700716	−0.298	0.063	−0.361	0.156	0.213	0.630	0.787
	600–800	22	496233	−0.229	0.033	−0.263	0.167	0.207	0.626	0.793
	> 800	18	423048	−0.271	0.032	−0.303	0.160	0.207	0.633	0.793
TWI	0.35–1.63	104	1318729	0.347	−0.289	0.635	0.401	0.152	0.447	0.848	0.080
	1.63–2.41	66	1221212	−0.031	0.016	−0.048	0.275	0.206	0.519	0.794
	2.41–3.36	20	653109	−0.599	0.099	−0.698	0.156	0.224	0.621	0.776
	3.36–4.81	4	243618	−1.223	0.051	−1.274	0.083	0.213	0.703	0.787
	4.81–14.58	1	60031	−1.208	0.012	−1.221	0.085	0.205	0.710	0.795
NDVI	-0.13–0.13	2	110968	−1.130	0.022	−1.151	0.077	0.206	0.717	0.794	0.037
	0.13–0.24	15	350282	−0.264	0.026	−0.290	0.184	0.206	0.610	0.794
	0.24–0.30	41	864378	−0.162	0.048	−0.210	0.204	0.211	0.585	0.789
	0.30–0.37	87	1306577	0.177	−0.123	0.300	0.286	0.178	0.536	0.822
	0.37–0.58	50	864494	0.036	−0.012	0.049	0.249	0.199	0.553	0.801
Land use	Farm land	140	1951094	0.252	−0.449	0.701	0.421	0.108	0.471	0.892	0.207
	Forest land	54	1459960	−0.411	0.216	−0.627	0.217	0.210	0.572	0.790
	Grass land	0	42187	0.000	0.012	0.000	0.000	0.000	1.000	1.000
	Water	0	25016	0.000	0.007	0.000	0.000	0.000	1.000	1.000
	Residential areas	1	16239	0.099	0.000	0.100	0.362	0.169	0.469	0.831
	Bare land	0	2203	0.000	0.001	0.000	0.000	0.170	0.830	0.830
Soil	Type 1	31	440582	0.232	−0.038	0.271	0.188	0.096	0.716	0.904	0.194
	Type 2	12	208103	0.033	−0.002	0.036	0.154	0.099	0.747	0.901
	Type 3	4	96481	−0.296	0.007	−0.304	0.111	0.100	0.789	0.900
	Type 4	9	147720	0.088	−0.004	0.093	0.163	0.099	0.738	0.901
	Type 5	1	41304	−0.834	0.007	−0.841	0.065	0.100	0.835	0.900
	Type 6	0	3783	0.000	0.001	0.000	0.000	0.000	1.000	1.000
	Type 7	0	27855	0.000	0.008	0.000	0.000	0.000	1.000	1.000
	Type 8	54	1441599	−0.398	0.207	−0.605	0.100	0.123	0.777	0.877
	Type 9	0	68835	0.000	0.020	0.000	0.000	0.000	1.000	1.000
	Type 10	84	1020437	0.389	−0.218	0.608	0.220	0.080	0.700	0.920
Lithology	A	26	396261	0.163	−0.023	0.185	0.159	0.098	0.743	0.902	0.141
	B	43	424981	0.596	−0.120	0.715	0.245	0.089	0.666	0.911
	C	59	765703	0.323	−0.113	0.437	0.186	0.089	0.724	0.911
	D	17	252004	0.190	−0.016	0.207	0.163	0.099	0.738	0.901
	E	38	686719	−0.008	0.002	−0.010	0.134	0.100	0.766	0.900
	F	1	267436	−2.702	0.074	−2.777	0.009	0.108	0.883	0.892
	G	10	467170	−0.958	0.091	−1.048	0.052	0.110	0.839	0.890
	H	0	97599	0.000	0.028	0.000	0.000	0.000	1.000	1.000
	I	1	46018	−0.942	0.008	−0.951	0.053	0.101	0.846	0.899
	J	0	92808	0.000	0.027	0.000	0.000	0.000	1.000	1.000

.

$$\mathrm{TOL}=\frac{1}{\mathrm{VIF}}$$

(17)

$$\mathrm{VIF}=\frac{1}{1-R_i^2}$$

(18)

where R_i is the negative correlation coefficient of the ith independent variable that makes regression analysis on other independent variables.

4.2. 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. The coefficients of evaluation models.

Landslide Conditioning Factors	WoE-LR	EBF-LR	IoE-LR
Slope aspect	0.780	6.477	15.569
Altitude	1.016	8.123	11.110
Land use	0.271	1.560	2.973
Lithology	0.484	3.951	3.567
NDVI	0.720	2.710	17.231
Plan curvature	0.628	4.169	82.950
Profile curvature	2.878	22.826	1239.456
Distance to rivers	1.045	5.298	52.206
Distance to roads	0.467	2.634	7.665
Slope angle	3.995	20.993	12.357
Soil	0.204	2.520	2.100
TWI	0.774	3.204	11.377
Slope aspect	−1.856	−17.300	−17.304

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.

$$\begin{array}{cc}\hfill LS{I}_{EBF}=& (Altitud{e}_{Bel})+(Plancurvatur{e}_{Bel})+(Profilecurvatur{e}_{Bel})+\hfill \\ & (Slopeangl{e}_{Bel})+(Slopeaspec{t}_{Bel})+(Distancetoroad{s}_{Bel})+\hfill \\ & (Distancetoriver{s}_{Bel})+(TW{I}_{Bel})+(NDV{I}_{Bel})+(Landus{e}_{Bel})+\hfill \\ & (Soi{l}_{Bel})+(litholog{y}_{Bel})\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill LS{I}_{WoE}=& (Altitud{e}_C)+(Plancurvatur{e}_C)+(Profilecurvatur{e}_C)+\hfill \\ & (Slopeangl{e}_C)+(Slopeaspec{t}_C)+(Distancetoroad{s}_C)+\hfill \\ & (Distancetoriver{s}_C)+(TW{I}_C)+(NDV{I}_C)+(Landus{e}_C)+\hfill \\ & (Soi{l}_C)+(litholog{y}_C)\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill LS{I}_{IoE}=& (Altitude\ast 0.190)+(Plancurvature\ast 0.009)+\hfill \\ & (Profilecurvature\ast 0.004)+(Slopeangle\ast 0.229)+\hfill \\ & (Slopeaspect\ast 0.051)+(Distancetoroads\ast 0.062)+\hfill \\ & (Distancetorivers\ast 0.022)+(TWI\ast 0.080)+(NDVI\ast 0.037)+\hfill \\ & (Landuse\ast 0.207)+(Soil\ast 0.194)+(lithology\ast 0.141)\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill LS{I}_{EBF-LR}=& (Altitud{e}_{Bel}\ast 8.123)+(Plancurvatur{e}_{Bel}\ast 4.169)+\hfill \\ & (Profilecurvatur{e}_{Bel}\ast 22.826)+(Slopeangl{e}_{Bel}\ast 20.993)+\hfill \\ & (Slopeaspec{t}_{Bel}\ast (-17.300))+(Distancetoroad{s}_{Bel}\ast 2.634)+\hfill \\ & (Distancetoriver{s}_{Bel}\ast 5.298)+(TW{I}_{Bel}\ast 3.204)+(NDV{I}_{Bel}\ast 2.710)+\hfill \\ & (Landus{e}_{Bel}\ast 1.560)+(Soi{l}_{Bel}\ast 2.520)+(litholog{y}_{Bel}\ast 3.951)\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill LS{I}_{WoE-LR}=& (Altitud{e}_C\ast 1.016)+(Plancurvatur{e}_C\ast 0.628)+\hfill \\ & (Profilecurvatur{e}_C\ast 2.878)+(Slopeangl{e}_C\ast 3.995)+\hfill \\ & (Slopeaspec{t}_C\ast (-1.856))+(Distancetoroad{s}_C\ast 0.467)+\hfill \\ & (Distancetoriver{s}_C\ast 1.045)+(TW{I}_C\ast 0.774)+(NDV{I}_C\ast 0.720)+\hfill \\ & (Landus{e}_C\ast 0.271)+(Soi{l}_C\ast 0.204)+(litholog{y}_C\ast 0.484)\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill LS{I}_{IoE-LR}=& (Altitude\ast 2.111)+(Plancurvature\ast 0.747)+\hfill \\ & (Profilecurvature\ast 4.958)+(Slopeangle\ast 2.830)+\hfill \\ & (Slopeaspect\ast (-0.883))+(Distancetoroads\ast 0.475)+\hfill \\ & (Distancetorivers\ast 1.149)+(TWI\ast 0.910)+(NDVI\ast 0.638)+\hfill \\ & (Landuse\ast 0.615)+(Soil\ast 0.407)+(lithology\ast 0.502)\hfill \end{array}$$

(24)

4.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. Wilcoxon signed-rank test for evaluation models.

Model	M1	M2	M3	M4	M5	M6	M7	M8
Z	−12.651	−12.968	−14.119	−15.345	−13.752	−15.018	−15.077	−16.006
Asymp. Sig. (2-tailed)	0	0	0	0	0	0	0	0
Model	M9	M10	M11	M12	M13	M14	M15
Z	−14.900	−10.849	−12.203	−10.313	−2.577	−4.809	−3.680
Asymp. Sig. (2-tailed)	0	0	0	0	0	0	0

M1: WoE–EBF, M2: IoE–EBF, M3: EBF_LR–EBF, M4: WoE_LR–EBF, M5: IoE_LR–EBF, M6: IoE–WoE, M7: EBF_LR–WoE, M8: WoE_LR–WoE, M9: IoE_LR–WoE, M10: EBF_LR–IoE, M11: WoE_LR–IoE, M12: IoE_LR–IoE, M13: WoE_LR–EBF_LR, M14: IoE_LR–EBF_LR, M15: IoE_LR–WoE_LR.

). 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. ROC curves. AUC: Area under curve; SE: standard error; CI: confidence interval.

Variables	AUC	SE	95% CI
EBF model	0.791	0.0226	0.747 to 0.830
WoE model	0.753	0.0243	0.707 to 0.795
IoE model	0.778	0.0233	0.734 to 0.819
EBF-LR model	0.826	0.0207	0.784 to 0.862
WoE-LR model	0.792	0.0226	0.748 to 0.831
IoE-LR model	0.825	0.0208	0.784 to 0.861

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 (0.825), WoE-LR model (0.792), EBF model (0.791), IoE model (0.778), and WoE model (0.753). The other statistical methods show the similar results, such as EBF-LR model has the smallest standard error (SE) and 95% confidence interval (CI).

5. 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 W^- produced 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.

6. 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.