Assessment of Landslide Susceptibility Using Integrated Ensemble Fractal Dimension with Kernel Logistic Regression Model

The main aim of this study was to compare and evaluate the performance of fractal dimension as input data in the landslide susceptibility mapping of the Baota District, Yan’an City, China. First, a total of 632 points, including 316 landslide points and 316 non-landslide points, were located in the landslide inventory map. All points were divided into two parts according to the ratio of 70%:30%, with 70% (442) of the points used as the training dataset to train the models, and the remaining, namely the validation dataset, applied for validation. Second, 13 predisposing factors, including slope aspect, slope angle, altitude, lithology, mean annual precipitation (MAP), distance to rivers, distance to faults, distance to roads, normalized differential vegetation index (NDVI), topographic wetness index (TWI), plan curvature, profile curvature, and terrain roughness index (TRI), were selected. Then, the original numerical data, box-counting dimension, and correlation dimension corresponding to each predisposing factor were calculated to generate the input data and build three classification models, namely the kernel logistic regression model (KLR), kernel logistic regression based on box-counting dimension model (KLRbox-counting), and the kernel logistic regression based on correlation dimension model (KLRcorrelation). Next, the statistical indexes and the receiver operating characteristic (ROC) curve were employed to evaluate the models’ performance. Finally, the KLRcorrelation model had the highest area under the curve (AUC) values of 0.8984 and 0.9224, obtained by the training and validation datasets, respectively, indicating that the fractal dimension can be used as the input data for landslide susceptibility mapping with a better effect.


Introduction
Landslides are regarded as one of the most destructive and frequently occurring natural disasters in the world. Globally, landslides cause about 1200 deaths and 3.5 billion dollars of loss each year [1]. China is a high-incidence region for landslides. Every year, it is reported that around 8935 landslides occur in China and about 350 people lose their lives due to landslides. Due to the diversity of the geological environment, the vagaries of climate, and the uneven distribution of the population, the spatial distribution of landslide risk in China is not uniform, which increases the obstructions in landslide control [2].

Methodology
To build the landslide susceptibility model and obtain the LSM, there were four main steps in the present research: (1) Data preparation, including landslide inventory and a description of landslide predisposing factors; (2) landslide predisposing factor analysis, based on a series of indexes and methods; (3) landslide modeling using the KLR model, the KLR box-counting model, and the KLR correlation model; and (4) the models' performance evaluation.

Landslide Inventory
The landslide inventory map, which reflects the relationship between predisposing factors and landslide distribution, is considered as the most crucial and essential phase before landslide susceptibility modeling [46]. Generally, it can obtain an inventory of the landslide location, category, occurrence date, size, volume, and active state [47]. In this study, the landslide inventory map was produced using existing literature and reports, field survey data, and the results from the interpretation of aerial photographs ( Figure 1). There were 316 landslides including four debris flows, 295 rainfall-induced slides, and 17 falls in the landslide inventory map [48], and the largest plane proportion of landslides was approximately 11.4 × 10 4 m 2 , the minimum area was about 295 m 2 , and the average proportion was 61 m 2 . The centroid method was applied to convert these landslide pattern spots into points to represent landslide locations. For subsequent landslide susceptibility modeling, the same number of non-landslides locations were randomly generated on the landslide inventory map. Then, a total of 632 points were divided into two parts according to the ratio of 70%:30%, with 70% (442) of the points used as the training dataset to train the models, and the remaining, namely the validation dataset, were applied for validation.

Methodology
To build the landslide susceptibility model and obtain the LSM, there were four main steps in the present research: (1) Data preparation, including landslide inventory and a description of landslide predisposing factors; (2) landslide predisposing factor analysis, based on a series of indexes and methods; (3) landslide modeling using the KLR model, the KLRbox-counting model, and the KLRcorrelation model; and (4) the models' performance evaluation.

Landslide Inventory
The landslide inventory map, which reflects the relationship between predisposing factors and landslide distribution, is considered as the most crucial and essential phase before landslide susceptibility modeling [46]. Generally, it can obtain an inventory of the landslide location, category, occurrence date, size, volume, and active state [47]. In this study, the landslide inventory map was produced using existing literature and reports, field survey data, and the results from the interpretation of aerial photographs ( Figure 1). There were 316 landslides including four debris flows, 295 rainfall-induced slides, and 17 falls in the landslide inventory map [48], and the largest plane proportion of landslides was approximately 11.4 × 10 4 m 2 , the minimum area was about 295 m 2 , and the average proportion was 61 m 2 . The centroid method was applied to convert these landslide pattern spots into points to represent landslide locations. For subsequent landslide susceptibility modeling, the same number of non-landslides locations were randomly generated on the landslide inventory map. Then, a total of 632 points were divided into two parts according to the ratio of 70%:30%, with 70% (442) of the points used as the training dataset to train the models, and the remaining, namely the validation dataset, were applied for validation.

Landslide Predisposing Factors
The reasons behind the causes of landslide occurrence are complicated; so far, there have been no consistent comments with regard to the determination of landslide predisposing factors. In this case study, thirteen types of landslide predisposing factors, including slope aspect, slope angle,

Landslide Predisposing Factors
The reasons behind the causes of landslide occurrence are complicated; so far, there have been no consistent comments with regard to the determination of landslide predisposing factors. In this case study, thirteen types of landslide predisposing factors, including slope aspect, slope angle, altitude, lithology, mean annual precipitation (MAP), distance to rivers, distance to faults, distance to roads, normalized differential vegetation index (NDVI), topographic wetness index (TWI), plan curvature, profile curvature, and terrain roughness index (TRI), were employed according to observations in the wild and previous studies on the study area [49]. In addition, a 30 m-resolution digital elevation model (DEM) was used to extract the slope aspect, slope angle, altitude, TWI, TRI, and the plan curvature and profile curvature layers using ArcGIS tools. The lithology and MAP layers were produced based on 1:100,000 geological map and meteorological data collected from the local government. The GF-2 remote sensing image and the 1:50,000 topographical map were applied to construct the distance to rivers, distance to faults, distance to roads, and NDVI layers.
Slope aspect is a significant factor for slope stability and landslide distribution [50]. Different slope aspects receive different light radiation, which influences the water content of the soil. In this study, slope aspect was classified into nine directions using the natural break method as follows: flat, north, northeast, east, southeast, south, southwest, west, and northwest, respectively ( Figure 2a).
In general, the probability of landslide occurrence increases with the increase of slope angle, which may influence the slope shear stress, and is still considered as one of the essential landslide predisposing factors by many scholars [51]. In this study, slope angle was classified into five sections using the natural break method as follows: 0-10. Altitude is also an important predisposing factor for landslide occurrence [52]. The change of altitude affects the magnitude of slope stress and affects the potential energy of the landslide. Using the natural break method, the altitude value in the study area was classified into five ranges as follows: 848-1037.6823; 1037.6823-1128.4000; 1128.4000-1210.8706; 1210.8706-1298.8392; and 1298.8392-1549 m, respectively ( Figure 2c).
Lithology is considered as the material basis of landslide development and occurrence. The weathering resistance and strength of rock and soil are determined by the types of lithology. On the other hand, the type and feature of landslides differ depending on the combination of rock mass with different properties, hardness, and structure [53]. According to geologic ages and lithofacies ( Table 1), all of the geological units were reclassified into eight categories ( Figure 2d).
Tectonic movement is not only one of the important factors in evaluating the regional geological stability, but is also a pivotal factor in landslide occurrence [54]. For this study, the value of distance to faults was employed to quantify the impact of faults on landslide occurrence and was reclassified into five ranges as follows: 0-2000; 2000-4000; 4000-6000; 6000-8000; and >8000 m, respectively ( Figure 2e).
River erosion plays a key role in the development of landslides. Many scholars believe that the effect of erosion on landslide stability is mainly reflected in the weakening of resistance of the landslide front and the increase of the free surface [55]. Therefore, the value of distance to rivers was employed to quantify the impact of river erosion on landslide development and was reclassified into five ranges according to the field observations and local conditions as follows: 0-200; 200-400; 400-600; 600-800; and >800 m, respectively ( Figure 2f).
Human activity is a primary factor that triggers landslides, as road construction is mainly the performance of human activities. The excavation of the slope and the earthwork accumulation during the construction process changes the local geological environment, which will directly or indirectly trigger a landslide [56]. In this study, the value of distance to roads was used as one of the condition factors and reclassified into five ranges: 0-200, 200-400; 400-600; 600-800; and >800 m, respectively ( Figure 2g).
Rainfall is considered to be an important factor in landslides because the study area is covered by a large area of loess and the structure will become loose after the loess is immersed in water [57]. In this study, the value of MAP was employed to represent the influence of rainfall on landslides. The MAP was divided into six sections according to the intervals of 20 mm/yr as follows: <520; 520-540; 540-560; 560-580; 580-600; and >600 mm/yr, respectively ( Figure 2h).
Vegetation plays a positive role in the stability of landslides and can improve the shear strength of the soil, while increasing the stability of the slope [58]. According to the observations of extensive field investigation, the more vegetation there is, the lower the number of landslides. In light of this, the value of the NDVI, which reflects the degree of vegetation coverage, was reclassified into four ranges based on the natural break method as follows: −0.9315-0.0776; 0.0776-0.4087; 0.4087-0.5742; and 0.5742-2.8915, respectively (Figure 2i).
The slope stability can be influenced by the shape of the slope, which can be evaluated by its profile curvature and plan curvature [59]. The profile curvature is defined as the curvature of a contour line generated by the intersection of the vertical plane with the surface, whereas the plan curvature is defined as that with the horizontal plane [60]. In this study, the profile curvature was classified into five ranges using the natural break method: -15.1897 to -1.5337; -1.5337 to -0.4607; -0.4607-0.5146; 0.5146-1.8802; and 1.8802-9.6837, respectively ( Figure 2j). Then, a similar method was applied to divide the plan curvature into five ranges: -9.7777 to -1.8107; -1.8107 to -0.5629; -0.5629-0.3009; 0.3009-1.2608; and 1.2608-14.6991, respectively ( Figure 2k).
TWI is commonly used to reflect the water condition in soil [61]. The value of TWI was calculated through the DEM using Equation (1) where A denotes for the specific catchment's region, and B is the value of slope angle in the study area. TRI was applied to reflect the fluctuation in the surface and the extent of erosion [62]. In the present research, TRI was calculated using Equation (2) and was classified into five ranges: -4508 to -1874; -1874 to -176; -176-57; 57-2398; and 2398-10,418, based on the natural break method (Figure 2m).
Entropy 2019, 21 5 ranges based on the natural break method as follows: −0.9315-0.0776; 0.0776-0.4087; 0.4087-0.5742; and 0.5742-2.8915, respectively ( Figure 2i). The slope stability can be influenced by the shape of the slope, which can be evaluated by its profile curvature and plan curvature [59]. The profile curvature is defined as the curvature of a contour line generated by the intersection of the vertical plane with the surface, whereas the plan curvature is defined as that with the horizontal plane [60]. In this study, the profile curvature was classified into five ranges using the natural break method: -15.1897 to -1.5337; -1.5337 to -0.4607; -0.4607-0.5146; 0.5146-1.8802; and 1.8802-9.6837, respectively ( Figure 2j). Then, a similar method was applied to divide the plan curvature into five ranges: -9.7777 to -1.8107; -1.8107 to -0.5629; -0.5629-0.3009; 0.3009-1.2608; and 1.2608-14.6991, respectively ( Figure 2k).
TWI is commonly used to reflect the water condition in soil [61]. The value of TWI was calculated through the DEM using Equation (1)

TWI ln tan
where A denotes for the specific catchment's region, and B is the value of slope angle in the study area. TRI was applied to reflect the fluctuation in the surface and the extent of erosion [62]. In the present research, TRI was calculated using Equation (2) and was classified into five ranges: -4508 to -1874; -1874 to -176; -176-57; 57-2398; and 2398-10,418, based on the natural break method (Figure 2m).

Frequency Ratio
The input data required for the classification model used in this study were of the numerical type; however, the slope aspect and lithology are nominal variables, so it was necessary to use frequency ratio (FR) data to assign values for these three predisposing factors. The frequency ratio is defined as the ratio of the area where landslides have occurred to the total study region and is also the ratio of the landslide occurrence probabilities to the non-landslide occurrence for a given attribute [63]. The FR data can be calculated according to the following formula: where X and Y are the number of landslides in a domain for each class and the number of pixels in a domain for each class, respectively. X and Y stand for the number of total landslides and pixels in the study area, respectively.
In the current research, the slope aspect and lithology factors assigned by the FR values and the remaining 11 predisposing factors, with the original numerical data, were defined as dataset 1 , which was used to run the KLR model.

Box-Counting Dimension
The spatial distribution of landslides is commonly considered to be not uniform, but is instead clustered at different scales. The fractal dimension originating from Mandelbrot's fractal theory is the value that quantitatively measures the degree of spatial clustering of the landslides. There are many techniques to calculate the fractal dimension, such as the slit island method, box-counting method, and the semi-variance method [64]. The first technique applied in the current research was the box-counting method, and it was employed to calculate the box-counting dimension, which could be used as the input data for landslide susceptibility modeling.
The box-counting method is applicable to both point datasets and can also be used for the calculation of fractal dimension in the two-dimensional and three-dimensional space. The principle of this method is to use a square segmentation plane with side length ε to calculate the number of grids containing landslide points N(ε), then change the value of ε to re-divide the plane and calculate the number of grids corresponding to the distribution of landslide points to obtain the sequence of landslide point pair (ε, N(ε)). In the case where the value ε is reasonable, if the aforementioned sequences satisfy or approximately satisfy Equation (4), the box-counting dimension (D 1 ) is considered to exist.
Through the python circumstance, the values of the box-counting dimension for each predisposing factor were measured and are shown in Table 4. In addition, 13 predisposing factors assigned by the box-counting dimension values were named as dataset 2 , which was used to run the KLR box-counting model.

Correlation Dimension
The second fractal dimension used as input data for landslide susceptibility modeling was the correlation dimension. The correlation dimension reveals the spatial fractal characteristics and regional differences of landslides from the perspective of the distance between the landslide points, and also reflects the degree of fragmentation of the geomorphological types in the study area. The calculation principle of the correlation dimension is to assume that the number of landslide points is N, then set Entropy 2019, 21, 218 9 of 23 a critical value r, determine the landslide point pair where the distance is less than r, and calculate its proportion in all landslide point pairs (N 2 ), as shown in the following formula: If r is set too large, then all points are less than r and C(r) = 1. Therefore, the value of r is gradually increased and the corresponding C(r) is calculated to obtain a set of sequences. If the above sequences satisfy or approximately satisfy Equation (7), the correlation dimension (D 2 ) is considered to exist.
Similarly, the values of box-counting dimension for each predisposing factor were measured based on the python circumstance (Table 4). A total of 13 predisposing factors, assigned by the correlation dimension values, were named as dataset 3 , which was used to run the KLR correlation model.

Multicollinearity Diagnosis
The premise of establishing a regression model is that each explanatory variable is independent of each other. If there is a strong linear correlation between the explanatory variables, it is considered that there is a multicollinearity problem among predisposing factors. The multicollinearity problem may lead to instability in the calculation of regression parameters, which will cause a major error in the results [65]. For these reasons, it is necessary to detect the potential multicollinearity problem between factors. In this study, two indicators obtained from the linear regression analysis, namely variance inflation factors (VIF) and tolerance (TOL), were employed to detect the potential multicollinearity problem. The VIF > 4 or TOL < 0.25 indicates a multicollinearity problem [66].

Selection of Predisposing Factors
In the process of landslide susceptibility modeling, not all predisposing factors have a positive influence on the accuracy of the classification modeling. In order to obtain a more accurate and reliable classification result, all of the predisposing factors needed to be filtered by estimating their contribution to the classification model [67]. For this reason, by calculating the information gain ratio (IG) of each predisposing factor to complete the filter process in this study, and the factors whose values of information gain ratio that are equal to or approximately equal to 0 must be excluded before landslide susceptibility modeling. The information gain ratio can be calculated using the following formulas: where D is the training dataset; Entropy(D) denotes the entropy of the training dataset; and y stands for the number of species in D. p k represents the proportion of category k in D. Then, the training dataset was divided into D v (v = 1, 2, 3, . . . , m) using s, which represents one of the predisposing factors, and we calculated the Gain(D, s) using Equation (9).
Entropy 2019, 21, 218 10 of 23 The information gain ratio for predisposing factor s is computed as: where IV(s) can be obtained by Equation (11).

Description of the KLR Model
The classification model selected in the current research to construct the landslide susceptibility modeling was a kernel logistic regression model (KLR). KLR is considered as a kernel version of logistic regression [68]. The main principle of the KLR model is to use a kernel function to perform logistic regression operations in high-dimensional feature space on data that are difficult to divide in the current dimensional space [69]. In this study, we took the landslide predisposing factors as input vector x and used a kernel function ϕ to complete the non-linear transformation of x. Accordingly, the non-linear form of logistic regression can be carried out as follows: where w and b are preferred by minimizing a cost function to represent the optimal parameters of the model, and p is the probability of landslide occurrence. The logit form of Equation (12) can be written as: The aforementioned kernel function is defined as the inner product between the images of vectors in the feature space. K There are several kernel functions that have been suggested such as the polynomial kernel, the linear kernel, the radial basis function (RBF), and the sigmoid kernel [70]. In the present research, the kernel function used for modeling was the RBF kernel, which can be written as follows: The kernel sensitivity is controlled by the turning parameter δ [71].

Statistical Index
In this study, the cut-off values were used in the final landslide susceptibility mapping to reclassify the landslide susceptibility index (LSI) into one of the response levels; however, the phenomenon of misclassification always exists in the LSM [72]. In order to evaluate the performance of classification models, six statistical indexes including the positive predictive rate (PPR), negative predictive rate (NPR), sensitivity, specificity, accuracy (ACC), and kappa index were employed as the assessment criteria, and these statistical indexes have frequently been used in many studies [39,73,74]. The PPR, NPR, sensitivity, specificity, and ACC can be calculated based on four basic indexes: the true positive (TP), true negative (TN), false positive (FP), and false negative (FN), as follows: where TP and TN denote the number of pixels which are correctly classified and FN and FP represent the number of pixels which are incorrectly classified. The kappa index can express the reliability of the classification model, and its calculation process is as follows: where n represents the total pixels of the training datasets [75].

The Receiver Operating Characteristic (ROC) Curve
Model comparison is considered as a significant step in landslide susceptibility modeling. In this study, the ROC curve, which is considered to be the most popular and widely used method of comparison models in landslide susceptibility modeling, was applied for assessing the classification model [76]. The x-axis and y-axis of the ROC curve are 1-specificity and sensitivity, respectively. The model comparison was undertaken by measuring the value of the area under the ROC curve (AUC), and the calculation formula of AUC is as follows: where P and N denote for the total number of landslides and non-landslides in the study area, respectively.

Multicollinearity Diagnosis
In order to detect the potential multicollinearity problems between landslide predisposing factors, the VIF and TOL of dataset 1 , dataset 2 , and dataset 3 were obtained through linear regression modeling [77]. For dataset 1 , it was observed from Table 2 that the maximum VIF value (1.7055) and the minimum TOL value (0.5863) belonged to the distance to rivers. For dataset 2 , the maximum VIF value (1.2358) and the minimum TOL value (0.8092) belonged to the distance to faults. For dataset 3 , the slope angle had the maximum VIF value and the minimum TOL value, which were 1.2546 and 0.7971, respectively. As a result, the VIF and TOL values of 13 predisposing factors were not within the range of VIF > 4 or TOL < 0.25, indicating that there were no potential multicollinearity problems in dataset 1 , dataset 2 , and dataset 3 .

Predisposing Factors Optimization
In this study, the contribution of predisposing factors for the classification model was quantified by calculating the average merit (AM) as the average IG values using the 10-fold cross-validation. As shown in Table 3, it was obvious that 13 predisposing factors in dataset 2 and dataset 3 had a positive contribution to build the classification model (AM > 0). In contrast, the AM values of the TWI, profile curvature, and TRI in dataset 1 were equal to 0, which means that these three predisposing factors in dataset 1 had no predictive ability in landslide susceptibility modeling. For this reason, the TWI, profile curvature, and TRI were abandoned from dataset 1 . The FR values of slope aspect and lithology factors and the classification of all predisposing factors are shown in Table 4. The FR value reveals the density of the landslide distribution, and the higher the FR value, the greater the density of the landslide distribution. In the case of slope aspect, the maximum value of FR (1.9024) appeared in the southeast, followed by the south (1.7262), and the east (1.1692), while the minimum FR value was north (0.5845). For lithology, category D had the highest FR value (19.5595), followed by category F with the FR value of 5.0326.  Dataset 1 was used as the input data to run the KLR model. The LSI values ranged from 0.0001 to 0.9999. Then, ArcGIS software was applied to visualize the LSI, which should be divided into different ranges to generate the LSM [78]. There are different types of classification schemes such as natural break, quantile, interval, standard deviation, and geometrical interval in ArcGIS software. In the current research, according to the geometrical interval method, the LSI of KLR model was divided into five categories: very low (0.0015-0.2404); low (0.2405-0.3931); moderate (0.3932-0.5615); high (0.5616-0.7494); and very high (0.7495-0.9674). The final LSM of the KLR model is shown in Figure 3a.

Integration of the KLR Model and Fractal Dimension
The acquired box-counting dimensions of dataset 2 and the correlation dimensions of dataset 3 are listed in Table 4. For slope angle, the highest box-counting dimension (0.4924) appeared in the section of 18.6711-25.7839 • , and the maximum correlation dimension (0.6981) also appeared in this section. In terms of the slope aspect, the class of west had the highest box-counting dimension Similarly, dataset 3 was also employed as the input data to run the KLR correlation model. The LSI values of KLR correlation model were in the range of 0.0001-0.9999. Then, the LSM of the KLR correlation model was produced by dividing the LSI values into five categories using the geometrical interval method (Figure 3c). The final threshold segmentations of LSI were as follows: very low (0.0866-0.3878); low (0.3879-0.5159); moderate (0.5160-0.5704); high (0.5705-0.6986), and very high (0.6987-0.9998).
Entropy 2019, 21 15 research, according to the geometrical interval method, the LSI of KLR model was divided into five categories:

Model Performance
In order to evaluate the performance of the classification models, six statistical indexes including PPR, NPR, sensitivity, specificity, ACC, and kappa index were calculated using the training datasets from dataset 1 , dataset 2 , and dataset 3 . As shown in Table 5, the KLR correlation model yielded the highest PPR, NPR, and ACC of 87.84%, 80.09%, and 83.97%, respectively. For sensitivity, the KLR correlation model showed the best performance for the classification of landslides (sensitivity = 81.59%), followed by the KLR box-counting model (sensitivity = 78.30%), and the KLR model (sensitivity = 66.03%). In terms of specificity, the KLR box-counting model showed the best performance for the classification of non-landslides (specificity = 86.76%), followed by the KLR correlation model (specificity = 87.44%), and the KLR model (sensitivity = 81.67%). Moreover, according to the criteria of the kappa index given from [79]: poor (<0); slight (0-0.2); fair (0.2-0.4); moderate (0.4-0.6); substantial (0.6-0.8); and perfect (0.8-1.0), the KLR box-counting model (kappa index = 0.7657) and the KLR correlation model (kappa index = 0.7828) expressed a substantial reliability. Unfortunately, the KLR model (kappa index = 0.5966) only showed a moderate reliability.

Model Validation
In this study, the results of model validation using the validation datasets from dataset 1 , dataset 2 , and dataset 3 are shown in Table 6. The maximum PPR (86.67%), NPR (90.59%), and ACC (88.42%) appeared in the KLR correlation model. For sensitivity, the KLR correlation model expressed the best performance for the classification of landslide (sensitivity = 91.92%), followed by the KLR box-counting model (sensitivity = 83.67%), and the KLR model (sensitivity = 70.41%). For specificity, the KLR box-counting model showed the best performance for the classification of non-landslide (specificity = 85.87%), followed by the KLR correlation model (specificity = 84.62%), and the KLR model (specificity = 79.35%). Furthermore, the kappa indexes of the KLR box-counting model, KLR correlation model, and KLR model were 0.8400, 0.8785, and 0.7336, respectively, indicating a substantial reliability between the reality and models.

Model Comparison
In this study, the model comparison was completed using the AUC value from the ROC curve. Figure 4a shows the final ROC curves and AUC values produced by the training datasets. The KLR correlation model expressed the maximum AUC value of 0.8984, followed by the KLR box-counting model with the AUC value of 0.8828, and the KLR model with the AUC value of 0.8352.
Additionally, the ROC curves and AUC values produced by the validation datasets are shown in Figure 4b. The KLR correlation model showed the maximum AUC value of 0.9224, followed by the KLR box-counting model with the AUC value of 0.9203, and the KLR model with the AUC value of 0.8605.

Discussion
The calculated box-counting dimensions and correlation dimensions in this study are listed in Table 4. The value range of the box-counting dimensions was between 0.9261 and 4.6410, while the correlation dimensions ranged from 1.4166 to 6.1590. Although the dimensions of the two fractal methods were different, it can be observed from Figure 5 that the overall trend of variation in the fractal was roughly the same. This indicates that the spatial distribution features of the landslide measured by the two fractal methods were relatively stable and the results more reliable. On the other hand, using the fractal dimension to optimize the predisposing factors may become a new approach that needs to be explored in future research.

Discussion
The calculated box-counting dimensions and correlation dimensions in this study are listed in Table 4. The value range of the box-counting dimensions was between 0.9261 and 4.6410, while the correlation dimensions ranged from 1.4166 to 6.1590. Although the dimensions of the two fractal methods were different, it can be observed from Figure 5 that the overall trend of variation in the fractal was roughly the same. This indicates that the spatial distribution features of the landslide measured by the two fractal methods were relatively stable and the results more reliable. On the other hand, using the fractal dimension to optimize the predisposing factors may become a new approach that needs to be explored in future research. methods were different, it can be observed from Figure 5 that the overall trend of variation in the fractal was roughly the same. This indicates that the spatial distribution features of the landslide measured by the two fractal methods were relatively stable and the results more reliable. On the other hand, using the fractal dimension to optimize the predisposing factors may become a new approach that needs to be explored in future research. Before building the classification models, the potential multicollinearity problems of dataset1, dataset2, and dataset3 were detected. All predisposing factors in these three datasets were independent of each other; however, the difference between dataset1, dataset2, and dataset3 can also be seen from Table 2. In terms of dataset1, the TOL values of altitude, distance to roads, distance to rivers, and NDVI were less than 0.7, which seems to indicate that these four factors had a tendency to have multicollinearity problems [80]. Moreover, if these four factors are excluded, it may affect the diversification of the input data. In contrast, the TOL values of all factors in dataset2 and dataset3 were greater than 0.7, which means that each factor had strong independence as the input data. In addition, from the results of the factor optimization shown in Table 3, three factors including TWI, profile Before building the classification models, the potential multicollinearity problems of dataset 1 , dataset 2 , and dataset 3 were detected. All predisposing factors in these three datasets were independent of each other; however, the difference between dataset 1 , dataset 2 , and dataset 3 can also be seen from Table 2. In terms of dataset 1 , the TOL values of altitude, distance to roads, distance to rivers, and NDVI were less than 0.7, which seems to indicate that these four factors had a tendency to have multicollinearity problems [80]. Moreover, if these four factors are excluded, it may affect the diversification of the input data. In contrast, the TOL values of all factors in dataset 2 and dataset 3 were greater than 0.7, which means that each factor had strong independence as the input data. In addition, from the results of the factor optimization shown in Table 3, three factors including TWI, profile curvature, and TRI in dataset 1 were excluded, but all predisposing factors in dataset 2 and dataset 3 were retained. In summary, dataset 2 and dataset 3 , which were constructed by the fractal dimension, can maintain a multiplicity of predisposing factors, while dataset 1 cannot. The basic classification model used in this study was the KLR model, which is considered as one of the state-of-the art advanced machine learning algorithms [81,82]. Meanwhile the KLR model has been used in landslide susceptibility mapping with high accuracy. However, an exploration of improving the KLR model has seldom been carried out. We used the fractal dimension as the input data of the KLR model for the first time, and the grid search method was applied to ensure that the parameters in the RBF kernel function were optimal at the same time. For model evaluation and comparison, the KLR correlation model constructed by dataset 3 performed the best, and its AUC values generated by the training dataset and validation dataset were the highest in the three models. Furthermore, the AUC values generated by the KLR model were significantly smaller than the other two models, which may be caused by the excessive difference in the dimension of the original data.

Conclusions
With the increasing threat of landslides to human beings, the prediction of landslide occurrence is particularly important. Landslide susceptibility mapping is considered as one of the preliminary steps to predict landslide occurrence, the main aim of which is to divide a specified region into multiple classes that range from stable to unstable ones. In this study, to obtain the landslide susceptibility map (LSM), thirteen predisposing factors (i.e., slope aspect, slope angle, altitude, lithology, mean annual precipitation (MAP), distance to rivers, distance to faults, distance to roads, normalized differential vegetation index (NDVI), topographic wetness index (TWI), plan curvature, profile curvature, and terrain roughness index (TRI)) were selected. Then, the KLR model and two hybrid models, namely the KLR box-counting model and the KLR correlation model generated with box-counting dimension and correlation dimension as input data, were used to perform landslide susceptibility mapping in the Baota District, Yan'an City, China.
From the final results, the classification results of all classification models were relatively reliable. For statistical evaluation methods, the performances of the two hybrid models were better than the KLR model. For the result of model comparison, the KLR correlation model had the highest values for landslide susceptibility mapping.
As the final conclusion, the results in the present study proved that using the fractal dimension as input data to build the hybrid model is feasible for landslide susceptibility mapping in the study area, and could provide a reference for local landslide prevention and decision making.