Landslide Susceptibility Evaluation Using Hybrid Integration of Evidential Belief Function and Machine Learning Techniques

In this study, Random SubSpace-based classification and regression tree (RSCART) was introduced for landslide susceptibility modeling, and CART model and logistic regression (LR) model were used as benchmark models. 263 landslide locations in the study area were randomly divided into two parts (70/30) for training and validation of models. 14 landslide influencing factors were selected, such as slope angle, elevation, aspect, sediment transport index (STI), topographical wetness index (TWI), stream power index (SPI), profile curvature, plan curvature, distance to rivers, distance to road, soil, normalized difference vegetation index (NDVI), land use, and lithology. Finally, the hybrid RSCART model and two benchmark models were applied for landslide susceptibility modeling and the receiver operating characteristic curve method is used to evaluate the performance of the model. The susceptibility is quantitatively compared based on each pixel to reveal the system spatial pattern between susceptibility maps. At the same time, area under ROC curve (AUC) and landslide density analysis were used to estimate the prediction ability of landslide susceptibility map. The results showed that the RSCART model is the optimal model with the highest AUC values of 0.852 and 0.827, followed by LR and CART models. The results also illustrate that the hybrid model generally improves the prediction ability of a single landslide susceptibility model.


Introduction
Landslides are large-scale movements of rocks, mud, and gravel from the top to the bottom of the mountain [1]. According to statistics, five percent of the world's natural disasters were landslides from 1994 to 2013 [2]. Landslides pose serious risks to human life, environment, resources, and property. Therefore, for the sake of decreasing the hazards caused by landslides, landslide susceptibility evaluation is becoming a common topic.
Landslide susceptibility refers to the probability of landslide occurring in an area on account of local geological environmental factors [3,4]. Geographic information system (GIS) has been diffusely applied for evaluating landslide susceptibility in recent decades, and many methods have been proposed. These methods mainly include deterministic methods, traditional statistical methods, and machine learning technologies [4,5].
In deterministic models, the stability of the slope is in the form of a calculated safety factor, which is suitable for small watersheds [6,7]. Statistical models, such as evidential belief function (EBF) [8,9], Preparation landslides inventory includes the location of past and recent landslides [56]. It also o understand the type and timing of landslides [57]. This study constructed a la y including 263 landslide locations by consulting aerial photos and collecting h e events, of which most of the landslides were slides (201), the others included 62 fa ng to an analysis in the GIS environment, the largest landslide was more than 1 × 10 7 llest landslide was nearly 120 m 3 . Rainfall and human engineering activities, such a l construction, road engineering construction and water conservancy engineering, gering factors of landslides. To establish and verify the landslide model, the landslid ly divided into two parts: (1) 70% were used to construct the training dataset; g 30% were used to generate validation dataset. ed on literature review and data availability [42,56,59,60], 14 landslide influencing d in this study. Slope angle, elevation, aspect, sediment transport index (STI), topogr index (TWI), stream power index (SPI), profile curvature, plan curvature, distance to to roads, soil, NDVI, land use, and lithology were considered. Slope angle, elevation and slope aspect are often used in landslide susceptibility mapping [61][62][63][64][65]. In this study, slope angles were divided into six classes:  Table 1). The elevation of the study area was reclassified into seven categories with an interval of 100 m, such as <1000 m, 1000-1100 m, 1100-1200 m, 1200-1300 m, 1300-1400 m, 1400-1500 m, and >1500 m ( Figure 3b, Table 1). Slope aspect was divided into nine types as follows: flat, north, northeast, east, southeast, south, southwest, west, and northwest ( Figure 3c, Table 1).

Data Preparation
The landslides inventory includes the location of past and recent landslides [56]. It also allows people to understand the type and timing of landslides [57]. This study constructed a landslide inventory including 263 landslide locations by consulting aerial photos and collecting historical landslide events, of which most of the landslides were slides (201), the others included 62 falls [58]. According to an analysis in the GIS environment, the largest landslide was more than 1 × 10 7 m 3 , and the smallest landslide was nearly 120 m 3 . Rainfall and human engineering activities, such as urban and rural construction, road engineering construction and water conservancy engineering, are the main triggering factors of landslides. To establish and verify the landslide model, the landslides were randomly divided into two parts: (1) 70% were used to construct the training dataset; (2) the remaining 30% were used to generate validation dataset.
Based on literature review and data availability [42,56,59,60], 14 landslide influencing factors were used in this study. Slope angle, elevation, aspect, sediment transport index (STI), topographical wetness index (TWI), stream power index (SPI), profile curvature, plan curvature, distance to rivers, distance to roads, soil, NDVI, land use, and lithology were considered.
where α is the specific watershed area, and β is the slope. In this study, STI values were classified into five classes, such as <10, 10-20, 20-30, 30-40, and >40 ( Figure 3d, Table 1). TWI affects slope material by affecting soil moisture and groundwater flow. It can be represented by Equation (2) as follows [68]: where α is the area drained per unit contour length at a point, and β is the slope. In the study area, the TWI was also divided into five categories: 1-2, 2-3, 3-4, 4-5, and >5 ( Figure 3e, Table 1). The SPI is an important hydrologic factor, which shows the process of potential flow erosion [68]. It can be represented by Equation (3) as follows [68]: where α is the specific watershed area, and β is the slope. In the study area, the SPI was also divided into five categories: 1-2, 2-3, 3-4, 4-5 and >5 ( Figure 3f, Table 1). Profile curvature and plane curvature are important topographic factors. Profile curvature is defined as the curvature at any point is perpendicular to the slope [69]. Profile curvature is divided into five types: (−7.  Table 1). Plan curvature is described as the curvature of the contour formed by the intersection with the plane [70].  Table 1).
The distance to rivers and distance to roads are important factors in landslide susceptibility modeling [71,72]. The river causes the bedrock to produce enough topographic change through erosion. This process makes the slope more prone to block failure [73]. Distance to rivers was divided into five classes with an interval of 200 m (Figure 3i, Table 1). The natural conditions of the slope were destroyed during the construction of transportation facilities, and oil and gas development, and other human engineering activities [74]. In the present study, distance to roads was grouped into five buffing zones using an interval of 100 m (Figure 3j, Table 1).
With the increase of soil depth, the runoff velocity decreases. The unstable nature of shallow is more likely to lead to landslides [75]. There are four types of soil, such as: cultivated loessal soils, alluvial soils, red clay soils and water ( Figure 3k, Table 1).
NDVI can quantitatively estimate vegetation growth and biomass by measuring surface reflectance [76]. The NDVI value is calculated based on the following formula: where IR is the infrared band, and R is the red band. The NDVI value was reclassified into five categories: −0.15-0.11, 0.01-0.04, 0.04-0.07, 0.07-0.09, 0.09-0.31 ( Figure 3l, Table 1). Land use is affected by environmental changes and plays an important role in slope stability. In this paper, land use is divided into six categories, including farmland, forestland, grasslands, water bodies, residential areas and others ( Figure 3m, Table 1).
Lithological units greatly influence the landslide occurrence [77]. The lithologic units were divided into five groups ( Figure 3n, Table 1): the first group is quaternary (loess, silt), the second group is tertiary (mudstone, conglomerate), the third group is cretaceous (arkose), the fourth group is Jurassic (shale, sandstone, mudstone, conglomerate), and the fifth group is Triassic (mudstone, sandstone, conglomerate).
Finally, all the 14 landslide influencing factors were converted to the same resolution of 30 m × 30 m.

Evidential Belief Function (EBF)
The EBF is an evidence algorithm which is a mathematical method based on bivariate statistics [78]. The EBF is a flexible model, because it not only accepts uncertainty but also absorbs the belief of many sources. The EBF model is made up of four mathematical functions: the degree of belief (Bel), the degree of disbelief (Dis), the degree of uncertainty (Unc), and the degree of plausibility (Pls). 0 to 1 is the range of these functions [79]. This study uses Bel values to represent the correlation between landslides and factors, and Bel values can be represented by the following expression: where Bel i is the extent of belief of ith influencing factor.

Classification and Regression Tree (CART)
The CART is an effective method because it has attested to be a technique for dealing with difficult classification problems [80]. CART is a decision tree proposed by Breiman [81]. The data is processed in a recursive form by the CART model [82]. In this process, the value of internal node features is "yes" and "no". Categorical and continuous dependent variables (regression) are predicted by classification and regression trees.
If the dependent variable is a category scale, CART will generate a classification tree, and if the dependent variable is a continuous data, CART will generate a regression tree [83]. The classification tree using CART can be constructed in four main steps: (1) building tree, (2) stopping tree building, (3) pruning tree, and (4) selecting optimal tree [84]. The information that needs to be processed is represented by CART in an intuitive way, and the metrics between the predictive variables do not affect the model results [53].
CART has many advantages as a classifier that it is not only a method suitable for numeric data types, but also generates an invariant for the transformation of independent variables. In addition, CART does not have to pick a variable first [85].

Random Subspace (RS)
The RS was proposed by Ho in 1998, which is defined as a classical integration algorithm [86]. The RS can establish the training subset, which is randomly selected from the original training set [87][88][89][90][91][92]. Also, the RS can modify the training data set in the feature space. RS is also a very effective integration algorithm because it can solve the data sets with redundant features and over-fitting problems. The RS enables the training data to maintain the highest accuracy. However, it cannot be ignored is that the growth of training data increases the complicacy of accuracy [93].
The main goal of RS is to collect the feature set of the high-dimensional feature space into the low-dimensional subspace, and then construct a classifier to classify the class based on the subspace. The final result is obtained by the majority voting rule [94]. More specifically, let the n-dimensional vector X in be the n features (landslide influencing factors) of the training sample X i . In RS, randomly select r < n features from the n-dimensional data set of the original space X to obtain the random subspace of r dimension. The modified training data set X e contains the r-dimensional training object X e i . Finally, the classifier is constructed in the RS X e and the majority voting rule is adopted, as follows [91]: where y ∈ (0, 1) is a decision of the classifier, δ ij is the Kronecker symbol, C a (x) are the generated classifiers (a = 1,2, . . . , A) [91].

Logistic Regression (LR)
Logistic regression can analyze a series of problems whose results are affected by one or more factors. The factors influencing the results are called independent variables, which can be continuous, discrete, or a combination of the two types. Logistic regression coefficients can figure out the influence of independent variables on landslide occurrence [95,96]. LR is defined as the follows: where A is the linear combination, β 0 is the intercept, β 1 , β 2 , . . . β n are the coefficients of logistic regression, and Y 1 , Y 2 , . . . Y n are the independent variables [97]. Using Equation (3), the landslide probability P is expressed as:

Correlation Analysis of Influencing Factors
In the present study, Bel value was used to represent the correlation between landslide and various landslide influencing factors. The results (Table 1) showed that Bel value in the south is the highest (0.181), followed by Bel value in the southwest (0.165), and Bel value in the southeast (0.151). At an altitude of more than 1500 m (0.305) greater regional influence landslide occurred right. The results indicate that slopes between 40-50 have a greater impact on the occurrence of landslides (Bel = 0.309). In terms of plan curvature, Bel value is the highest for the class of 1.44-7.56 (0.279).
In terms of profile curvature, the class of 0.58-1.97 has the highest value (0.288), so it has a greater impact on the landslide occurrence. TWI value between 2 and 3 has the highest Bel value, which is 0.388. Areas within 200 m to rivers are more prone to landslides (Bel = 0.590) because rivers provide wetter soil for landslides. Similarly, areas within 100 m to the road networks have a greater impact on landslides occurrence (Bel = 0.313). This is because the soil near the road becomes more loosely structured during human activity. For the NDVI, the area between 0.07 and 0.09 has a higher impact on landslide occurrence. Residential areas have higher impact on landslide events (Bel = 0.385). In terms of lithology, Bel value of the fourth group is 0.31, indicating that Jurassic rocks were more sensitive to landslide occurrences. The area covered by red clay has a greater impact on landslide occurrence (Bel = 0.614). The highest Bel values of SPI and STI are 0.234 and 0.219, respectively, and the corresponding ranges are 20-30 and 30-40, respectively. The regions within these ranges have a greater impact on the landslide occurrence.
The multicollinearity test of landslide influencing factors is very important in landslide susceptibility mapping. The two most widely used indices in multicollinearity analysis are tolerance and variance inflation factors (VIF). If the tolerance value is <0.1 or the VIF value is >10, it indicates that there is serious multicollinearity among landslide influencing factors [98][99][100]. Tolerance and VIF values of 14 landslide influencing factors show that the distance to the river has the smallest tolerance value of 0.715 (>0.1), and the largest VIF value is 1.399 (<10) ( Table 2). Therefore, there is no multicollinearity among the 14 landslide influencing factors. At the same time, this study also evaluated the importance of 14 landslide influencing factors using the 10-fold cross-validation correlation attribute evaluation method (CAE) in Weka software [101]. The evaluation results are sorted in descending order according to the average merit ( Table 3

Application of Hybrid and Benchmark Model
RS integration can not only use random subspace to construct and aggregate the basic classifier, but also make the basic classifier easier to train than in the smaller subspaces. Therefore, the performance of the base classifier in the RS may be better than the original feature space, and the feature-to-instance ratio is significantly improved. The RSCART model is constructed by using Weka software [102] through optimization and classification steps. In the optimization step, RS is trained to obtain the sub-data set, which is randomly divided by the original data set. The training data set is divided by the iterative method. In the classification step, combined with the CART algorithm, the optimal training data set obtained in the previous step is used for spatial prediction of landslide [28]. In the process of RSCART model construction, the optimal number of iterations was 28, the optimal number of execution slots was 1. Finally, the landslide susceptibility map is developed by using the RSCART model and reclassified into five classes using the natural break method [103,104]. The very high class occupies a very small area of 9.27%. The moderate, very low, low, and high areas are accounting for 30 According to Equation (8), logistic regression coefficients of all landslide influencing factors are positive (Table 4), which indicate that all influencing factors are positively correlated with the landslide occurrence. Finally, the landslide susceptibility map was divided into five categories according to the natural break method. It can be seen that the area percentage with low susceptibility is the largest (29.42%), followed by moderate, very low and high (22.13%, 22%, and 14.38%, respectively). The area proportion of very high is the smallest (12.06%) (Figures 4c and 5) optimal training data set obtained in the previous step is used for spatial prediction of landslide [28]. In the process of RSCART model construction, the optimal number of iterations was 28, the optimal number of execution slots was 1. Finally, the landslide susceptibility map is developed by using the RSCART model and reclassified into five classes using the natural break method [103,104]. The very high class occupies a very small area of 9.27%. The moderate, very low, low, and high areas are accounting for 30.56%, 13.31%, 27.22% and 19.63%, respectively (Figures 4a and 5).

Validation and Comparison of Models
The landslide susceptibility map should have the ability to verify with existing landslide data and predict future landslides [105]. Therefore, in this study, the ROC curve and the AUC are used to assess the prediction capability of models [106][107][108]. The best models tend to have the highest AUC among the models studied [4,109]. ROC curves and AUC values of the training dataset of the three models are shown in Figure 6. The RSCART model has the highest AUC value (0.852), followed by LR model (0.797) and CART model (0.793). The ROC curves and AUC values of the validation dataset are shown in Figure 7. The prediction accuracy of RSCART model (AUC = 0.827) is higher than that of the LR model (AUC = 0.758) and CART model (AUC = 0.749). Therefore, the comparison of AUC values indicates that the RSCART model is the best of the three models. In addition, the landslide susceptibility map can be verified by calculating the landslide density. LD is defined as the ratio of the percentage of landslide points in each susceptibility classification to the percentage in each susceptibility classification [110]. The landslide density calculation results of the three models are shown in Table 5

Comparison of Landslide Susceptibility Maps
This study also quantitatively compares the susceptibility values on each pixel to reveal the systematic spatial pattern of the differences between susceptibility maps, following the methodology proposed by [111]. The RSCART model was selected as the benchmark because it has a higher AUC value than the other two models. The susceptibility map of the baseline model is used in a GIS system to pair with the remaining models and subtract their values to define their differences ( Figure 8). The

Comparison of Landslide Susceptibility Maps
This study also quantitatively compares the susceptibility values on each pixel to reveal the systematic spatial pattern of the differences between susceptibility maps, following the methodology proposed by [111]. The RSCART model was selected as the benchmark because it has a higher AUC value than the other two models. The susceptibility map of the baseline model is used in a GIS system to pair with the remaining models and subtract their values to define their differences ( Figure 8). The values of the comparison map are divided into three levels: "underestimation", "approximation" and "overestimation". The values of both comparison maps are broken at −0.2 and 0.2. The percentage of each grade in the total area is shown in Table 6. At the same time, to explore the key factors influencing the susceptibility difference, overestimation and underestimation statistics were performed for each category of each impact factor (Tables 7 and 8). For each class, we calculate "A" as a percentage of the total area of each class. "B" is the ratio of the underestimated (overestimated) pixels found in this class to the total underestimated (overestimated) pixels, "B-A", as the difference ratio between the two maps, can be used to identify key class of underestimation (overestimation) anomaly clustering. According to the "B-A" value defined for each class, the class with the highest degree of imbalance was identified (Table 8). To be able to clearly illustrate the relationship between the most imbalanced classification and the underestimated or overestimated area, the visual inspection is required. Underestimations of "RSCART-LR" driven by slope angle (40 • -50 • ) (Figure 9a), overestimations of "RSCART-LR" driven by slope angle (<10 • ) (Figure 10a), Underestimations of "RSCART-CART" driven by distance to rivers (0-200 m) (Figure 9d) and overestimations of "RSCART-LR" driven by slope angle (<10 • ) ( Figure 10d). As shown in Figure 9a, almost all the underestimated pixels are in the classification with slope of 40-50, and the percentage is 98.55%. In Figure 10a, 99.35% of the overestimation pixels are in the slope of 0 • -10 • . In Figure 9d, all underestimation pixels are clustered in a distance to rivers of 0-200 m. In Figure 10d, 76.87% of the overestimation pixels are in the classification with a slope less than 10 • . degree of imbalance was identified (Table 8). To be able to clearly illustrate the relationship between the most imbalanced classification and the underestimated or overestimated area, the visual inspection is required. Underestimations of "RSCART-LR" driven by slope angle (40°-50°) ( Figure  9a), overestimations of "RSCART-LR" driven by slope angle (<10°) (Figure 10a), Underestimations of "RSCART-CART" driven by distance to rivers (0-200 m) (Figure 9d) and overestimations of "RSCART-LR" driven by slope angle (<10°) (Figure 10d). As shown in Figure 9a, almost all the underestimated pixels are in the classification with slope of 40-50, and the percentage is 98.55%. In Figure 10a, 99.35% of the overestimation pixels are in the slope of 0°-10°. In Figure 9d, all underestimation pixels are clustered in a distance to rivers of 0-200 m. In Figure 10d, 76.87% of the overestimation pixels are in the classification with a slope less than 10°.    Table 7. Statistics on underestimation pixels and overestimation pixels of RSCART-LR.    Table 7. Statistics on underestimation pixels and overestimation pixels of RSCART-LR.

Discussion
The selection of landslide influencing factors will affect the quality of landslide susceptibility analysis [112,113]. In this study, 14 landslide influencing factors were selected. The EBF model is used to analyze the correlation between the subclasses of landslide influencing factors and landslide, and to reclassify the landslide influence factors. Then, through the multicollinearity analysis, it was found

Discussion
The selection of landslide influencing factors will affect the quality of landslide susceptibility analysis [112,113]. In this study, 14 landslide influencing factors were selected. The EBF model is used to analyze the correlation between the subclasses of landslide influencing factors and landslide, and to reclassify the landslide influence factors. Then, through the multicollinearity analysis, it was found that there was no multicollinearity among all landslide influencing factors (Table 2), and all factors contributed to the model. Finally, the results of the CAE model's importance analysis and Bel values were used to analyze the various landslide influencing factors. According to the results of the importance analysis, the distance to rivers (AM = 0.378) is the most important landslide influencing factor in the study area. The results of the Bel value are similar to those found in previous studies, with a higher probability of landslides near the river area [114]. Different geomorphic parts of the river cause different external forces and stress distributions on the slope, so the failure modes of the slope are also different. The distance to roads (AM = 0.173, Bel = 0.313) has similar results. The closer we get to the road, the higher the probability of a landslide [115]. This is easy to understand because road construction can destabilize the slope by breaking the support of the slope foot [116]. As for the slope (AM = 0.213), the possibility of landslides will increase with the increase of the slope, because the gravity load and stress of the material forming the slope will increase [117]. It can be seen from the Bel value that landslides are most likely to occur in areas with a slope of 40-50. Different types of lithology (AM = 0.181) will result in different slope strengths [118]. According to the Bel value, it can be seen that shale, sandstone, mudstone, and conglomerate have a certain effect on the occurrence of landslides. Elevation (AM = 0.172) has always been a key factor in landslide susceptibility mapping [119,120]. The weathering and shear strength of rocks at different elevations are also different [119,120]. According to the Bel value, it can be concluded that the possibility of landslides of 1500 m-1574 m in the study area is higher. Therefore, it can be judged that landslides are more likely to occur in higher research areas. TWI (AM = 0.171) and SPI (AM = 0.154) are usually indispensable constraints in landslide susceptibility modeling [110,121,122]. Aspect (AM = 0.143, Bel = 0.181) in the south are more likely to cause landslides. This may be due to the more intense solar radiation when the slope is facing south, which may result in different weathering degree. Previous studies have shown that the susceptibility of different types of soil (AM = 0.143) to landslides is also different [101]. According to the Bel value, red clay is the key factor leading to landslides in the study area. The shear strength of red clay soils during rainwater infiltration is reduced, making slopes more prone to landslides [123]. Profile curvature (AM = 0.138) can affect the stress distribution of the slope, and the variations in curvature may be useful to identify depletion and accumulation zones [124]. NDVI (AM = 0.103) is an important influencing factor for landslides. Many studies have shown that plants play an active role in the occurrence of landslides because their root systems can increase soil strength and reduce water infiltration [125][126][127]. In addition, for those factors whose AM value is less than 1, many scholars have studied their relationship with the occurrence of landslides, and these factors should not be ignored when mapping landslide susceptibility [128][129][130].
From the results we can observe that the number of overestimation and underestimation is limited (the highest proportion is 0.057). However, it has a great influence on the value of AUC. This is because overestimation and underestimation are not randomly distributed, but there are some spatial patterns [111]. From Figure 2 we can see that the overestimations of the two comparisons have some similarities in spatial distribution. From Table 9 we can find that in the two comparisons, the overestimation has almost the same class of imbalance factors. Slopes less than 10 are classified as the most unbalanced class. From Figure 10, we can find that the classification with a slope of less than 10 has a high degree of overlap with the spatial distribution of the river. Areas with a slope of less than 10 are all closer to the river. For the "RSCART-LR" comparison map, Alluvial soils is also an unbalanced class. In Figure 10c, we can find that the spatial distribution of alluvial soil also overlaps with the river. This is because alluvial soil develops. The classification with the distance to the road less than 100 is also the overestimated imbalanced class of the two comparison maps (Figure 10b,f), and it has a high spatial overlap with the overestimated pixels. The class with a slope of 40-50 is the most unbalanced class in the "RSCART-LR" comparison, which contains 98.55% of the underestimated pixels (Figure 9a). At the same time, "RSCART-LR" was also significantly affected by the distance to the river less than 200 ("B-A" = 71.3%) (Figure 9b). In the comparison of "RSCART-CART", the classification with the river distance less than 200 included all the underestimated pixels ( Figure 9d) and was also driven by the slope classification of 40-50 ("B-A" = 41.42%). From Figure 9, we can see that the overlap of the spatial distribution of the underestimated pixels with a distance from the river less than 200 is extremely high. The remaining imbalance classes do not significantly dominate the underestimated spatial distribution. In summary, the integrated model RSCART can better use the landform information related to the river and the information closer to the road than the other two models. However, it is easy to ignore the influence of the river itself and the high slope on the susceptibility to landslides. Table 9.
Most imbalanced classes driving the spatial distribution of underestimations and overestimation.

Comparison Maps Imbalanced Classes
Underestimation RSCART-LR slope, 40-50; STI, 10-20; TWI,1.11-2; SPI, 30-40; distance to rivers, 0-200; distance to roads, >400 RSCART-CART slope, 40-50; elevation, 1000-1100; profile, 0.58-1.97; distance to rivers, 0-200; distance to roads, >400 Overestimation RSCART-LR slope, <10; STI, 0-10; TWI, 2-3; SPI, 0-10; distance to roads, 0-100; soil, group2; RSCART-CART slope, <10; STI, 0-10; SPI, 0-10; distance to roads, 0-100 Some studies have proved that LR model is an excellent model for landslide susceptibility research. Polykretis and Chalkias shows that the performance of LR model is better than that of weight of evidence and artificial neural networks in landslide susceptibility mapping in drainage basin of Selinous River [131]. Oh et al. conclude that the LR model has the highest prediction and training accuracy compared with EBF and support vector machine (SVM) in the region surrounding Yongin, South Korea [132]. According to the results of Luc Yen district by Pham et al., the performance of CART model is better than that of LR model [28]. The results also show that the performance of CART model is lower than that of LR model. In addition, some studies have shown that RS is an excellent ensemble algorithm. Hong et al. indicate that the integrated model RSSVM has the optimal performance [55]. In this study, training (goodness of fit) and validation (prediction accuracy) data sets were used to compare the proposed integration model with the benchmark model. The goodness of fit results of the landslide model show that RSCART (AUC = 0.852) is superior to CART (AUC = 0.793) and LR (AUC = 0.797). The validation results show that RSCART model has better prediction accuracy than CART (AUC = 0.749) and LR (AUC = 0.758) models. The RSCART model has the advantages of RS integration and CART classifier. When integrated with the RS model, the performance of the CART model, whose performance is lower than the benchmark LR model, is improved significantly. This is because RS integration constructs the benchmark model by using the random subspace, and the performance of the base classifier in the random subspace is optimized. Therefore, the RSCART model based on machine learning hybrid method is more effective than the single basic model. In addition, the landslide susceptibility figure can be verified by calculating the landslide density (Table 5). In general, higher LD values are associated with higher landslide susceptibility levels. As can be seen from the result of LD, the results of the three models are consistent. The LD value of the very high type is the highest, followed by high, moderate, low, and very low. Therefore, the three base maps of landslide susceptibility are reliable, and the optimal model RSCART also has the highest LD value. Therefore, the model used in this study to generate the landslide susceptibility map has reference significance for the study area. In addition, the selection of factors and the construction of models are also of reference values for similar studies.

Conclusions
In this study, the combined CART with RS and two benchmark models (LR and CART) were used to draw three landslide susceptibility maps for Zichang County. The correlation between influencing factors and landslide was evaluated using multicollinearity analysis and EBF method. The AUC value was used to test the performance of the three models. The validating results show that the RSCART model has the highest performance, followed by the LR model and CART model. Finally, this paper also uses a method to quantitatively compare the susceptibility values of each pixel to reveal the systematic spatial pattern of the differences between susceptibility maps. In conclusion, the landslide susceptibility maps compiled in this study are useful for land use and decision-making in landslide-prone areas. In addition, this study also proves the superiority of hybrid model in landslide susceptibility modeling.
Author Contributions: Y.L. and W.C. contributed equally to the work. Y.L. and W.C. collected field data and conducted the landslide susceptibility mapping and analysis. Y.L. and W.C. wrote and revised the manuscript. All the authors discussed the results and edited the manuscript. All authors have read and agreed to the published version of the manuscript.