Application of Machine Learning to Debris Flow Susceptibility Mapping along the China–Pakistan Karakoram Highway

: The China–Pakistan Karakoram Highway is an important land route from China to South Asia and the Middle East via Pakistan. Due to the extremely hazardous geological environment around the highway, landslides, debris ﬂows, collapses, and subsidence are frequent. Among them, debris ﬂows are one of the most serious geological hazards on the Karakoram Highway, and they often cause interruptions to tra ﬃ c and casualties. Therefore, the development of debris ﬂow susceptibility mapping along the highway can potentially facilitate its safe operation. In this study, we used remote sensing, GIS, and machine learning techniques to map debris ﬂow susceptibility along the Karakoram Highway in areas where observation data are scarce and di ﬃ cult to obtain by ﬁeld survey. First, the distribution of 544 catchments which are prone to debris ﬂow were identiﬁed through visual interpretation of remote sensing images. The factors inﬂuencing debris ﬂow susceptibility were then analyzed, and a total of 17 parameters related to geomorphology, soil materials, and triggering conditions were selected. Model training was based on multiple common machine learning methods, including Ensemble Methods, Gaussian Processes, Generalized Linear models, Navies Bayes, Nearest Neighbors, Support Vector Machines, Trees, Discriminant Analysis, and eXtreme Gradient Boosting. Support Vector Classiﬁcation (SVC) was chosen as the ﬁnal model after evaluation; its accuracy (ACC) was 0.91, and the area under the ROC curve (AUC) was 0.96. Among the factors involved in SVC, the Melton Ratio ( MR ) was the most important, followed by drainage density ( DD ), Hypsometric Integral ( HI ), and average slope ( AS ), indicating that geomorphic conditions play an important role in predicting debris ﬂow susceptibility in the study area. SVC was used to map debris ﬂow susceptibility in the study area, and the results will potentially facilitate the safe operation of the highway.


Introduction
The China-Pakistan Karakoram Highway (KKH) is an important land route from China to South Asia and the Middle East via Pakistan [1,2]. Debris flows are one of the most serious geological hazards affecting the KKH [3]. A total of 150 debris flows occurred along the KKH from 2008 to 2011 and were investigated by Liao et al. [4]. They are rapid, surging flows of water-charged clastic sediments flowing in a steep channel [5]. For example, on 2 March 2015, a debris flow occurred in section 1607 of the KKH, on the border between Aktau county and Shufu county in Kashi Prefecture, suitable for the study area. Furthermore, the impact of different factors on the susceptibility of debris flow through a machine learning model needs to be analyzed, and the machine learning method needs lower capital and calculation costs compared with a physical model. Therefore, this paper chooses a machine learning method to map debris flow susceptibility along KKH.
The purpose of the present study was to develop a methodology for rapidly mapping debris flow susceptibility along the KHH, using remote sensing, GIS, and machine learning technology. First, the inventory of catchments which are prone to debris flows was obtained by visual interpretation of remote sensing images and field survey, and then, the factors influencing debris flows were collected and processed. Finally, debris flow susceptibility mapping was performed based on multiple machine learning methods.

Study Area
The study area is located in the northwest part of the Himalayan orogenic belt in northern Pakistan and northwest China, passing through the Karakoram, Himalayan, and Hindu Kush mountain ranges ( Figure 1). The KKH starts from Kashgar in Xinjiang, China, passing through Khunjerab (at the border between China and Pakistan), the Himalayas, west of the Karakoram Mountains, and the northern edge of the Pamir Plateau, reaching Thakot in northern Pakistan. The total length exceeds 1000 km. From Khunjerab to the northern city of Thakot in Pakistan, the decrease in altitude is up to 4500 m, and the landform characteristics are unique, mainly mountains, valleys, glaciers, and ice margins.
The KKH spans a wide area; the geological environment is complex. The sedimentary and metamorphic rocks exposed along the road mainly comprise Proterozoic, Paleozoic Carboniferous, Permian, Mesozoic Triassic, and Cretaceous strata. The lithologies are mainly mudstone, limestone, and slate, and there are large variations in lithology. Differential weathering and erosion lead to steep topography and widely distributed loose deposits [36], providing favorable terrain conditions and abundant materials for landslides.
The study region is located in the transitional zone between the monsoon zone of Central Asia and South Asia. In the northern part of the area, there is a high-elevation inland plateau climate, with low precipitation, strong solar radiation, low temperature, extensive glacier coverage, and pronounced freeze-thaw weathering. The southern foothills of the mountainous area are influenced by water vapor derived from the Indian Ocean, and the rainfall is greater, with the annual precipitation reaching 600-1000 mm [7]. The rainfall is concentrated in summer, and the high elevation mountainous terrain leads to pronounced small-scale climatic variability and extreme rainfall events. The climate along the highway has a pronounced vertical differentiation, with hot, dry valley, alpine, and glacial climates, with strong contrasts. The average annual rainfall in valley areas is less than 150 mm and the average annual rainfall of the high-altitude mountain areas, where the majority of the rainstorms occur, reaches 2000 mm.
This area is located in the active collision zone of the Indian and Asian plates [37]. The significant tectonic features are the main Karakoram thrust and Karakoram fault [38]. The main Karakoram thrust is the collision zone on the southern margin of the Eurasian plate and extends to Baltistan through the Hashupa, Shigar, and Shyok valleys, and is an active thrust fault with a large angle, along which many earthquakes occurred [39] (Figure 1).
As the geological environment around the highway is very hazardous, geological and geomorphological processes such as landslides, mudslides, collapses, and subsidence along the highway are frequent; they are major obstacles to the smooth operation of the highway, often blocking traffic and seriously affecting its normal operation [7]. Debris flows on both sides of the China-Pakistan highway are widely developed, and they can be divided into snowmelt debris flows and flood debris flows. The snowmelt debris flows are concentrated in the glacial areas of Kungai Mountain and Gongge Mountain, on both sides of the highway, with an altitude of more than 4000 m. The water source is mainly snowmelt and the loose material mainly consists of moraine and old debris flow deposits [6]. Flood type debris flows are widely distributed in the study area. Generally, the formation area is large, Figure 1. Location of the study area and distributions of faults, the KKH, earthquake sites, and the major settlements. (The Digital Elevation Model (DEM) data were provided by the International Scientific and Technical Data Mirror Site, Computer Network Information Center, Chinese Academy of Sciences, http://www.gscloud.cn; faults were provided by China Geological Survey; the earthquake sites were from earthquake database supplied by United States Geological Survey (USGS), and the sites with magnitude above 4 on the Richter scale since 2000 were selected, https://earthquake.usgs.gov/earthquakes).

Methods
The modeling process included several steps: parameter selection and preparation, data collection and processing, debris flow catchment preparation, model selection, model fitting, and model evaluation. A flow chart of the process is illustrated in Figure 2.

Methods
The modeling process included several steps: parameter selection and preparation, data collection and processing, debris flow catchment preparation, model selection, model fitting, and model evaluation. A flow chart of the process is illustrated in Figure 2.

Catchment Boundaries Division
Catchments were selected as mapping units because of their significant conceptual and operational advantages [40,41]. Based on GIS and Soil & Water Assessment Tool (SWAT) tools, a DEM of the study area with a 60-m resolution (the DEM data were provided by the International Scientific and Technical Data Mirror Site, Computer Network Information Center, Chinese Academy of Sciences, http://www.gscloud.cn) was used for hydrological analysis and watershed boundary division. Taking the primary tributaries of the water system along the China-Pakistan highway as the division standard of the catchment boundary, 653 catchments were divided after being sorted manually, with the drainage area ranging from 0.74 to 3990 km 2 ( Figure 3).

Inventory of Debris Flows
From 6 to 12 June 2019, our team carried out a survey on the disasters along the highway; a rough survey of the debris flow was carried out for safety reasons. However, field investigation can only identify the obvious debris flow fans, which may not be able to effectively identify the debris flow accumulation fan that has been damaged. Under the conditions of a large study area, few historical references, and difficult field conditions, remote sensing interpretation is a feasible and attractive approach [42]. Therefore, Google Earth images were used to avoid the time limit of field investigation. In order to minimize the problem of debris flow accumulation fans in some catchments being damaged and unrecognized, historical aerial images visible in Google Earth (the time is from 2005 to 2020) were interpreted; the ancient debris flow accumulation fans were not interpreted, and only the debris flow catchments that were still active in recent years were considered. According to Chevalier [31], criteria such as vegetation's change in the landscape, landslide scar(s), clear visibility of a torrent/gully/stream where roughness could be assessed, or the presence of potential deposition fans were considered.

Catchment Boundaries Division
Catchments were selected as mapping units because of their significant conceptual and operational advantages [40,41]. Based on GIS and Soil & Water Assessment Tool (SWAT) tools, a DEM of the study area with a 60-m resolution (the DEM data were provided by the International Scientific and Technical Data Mirror Site, Computer Network Information Center, Chinese Academy of Sciences, http://www.gscloud.cn) was used for hydrological analysis and watershed boundary division. Taking the primary tributaries of the water system along the China-Pakistan highway as the division standard of the catchment boundary, 653 catchments were divided after being sorted manually, with the drainage area ranging from 0.74 to 3990 km 2 ( Figure 3).

Inventory of Debris Flows
From 6 to 12 June 2019, our team carried out a survey on the disasters along the highway; a rough survey of the debris flow was carried out for safety reasons. However, field investigation can only identify the obvious debris flow fans, which may not be able to effectively identify the debris flow accumulation fan that has been damaged. Under the conditions of a large study area, few historical references, and difficult field conditions, remote sensing interpretation is a feasible and attractive approach [42]. Therefore, Google Earth images were used to avoid the time limit of field investigation. In order to minimize the problem of debris flow accumulation fans in some catchments being damaged and unrecognized, historical aerial images visible in Google Earth (the time is from 2005 to 2020) were interpreted; the ancient debris flow accumulation fans were not interpreted, and only the debris flow catchments that were still active in recent years were considered. According to Chevalier [31], criteria such as vegetation's change in the landscape, landslide scar(s), clear visibility of a torrent/gully/stream where roughness could be assessed, or the presence of potential deposition fans were considered. The identification results are shown in Figure 3. A total of 544 catchments prone to debris flow (DFs) were identified in the study area, and the remaining 109 catchments were not prone to debris flow (NDFs).
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 22 The identification results are shown in Figure 3. A total of 544 catchments prone to debris flow (DFs) were identified in the study area, and the remaining 109 catchments were not prone to debris flow (NDFs).

Selection of Factors Influencing Debris Flows
The selection and processing of prediction factors have an important influence on the accuracy of debris flow susceptibility mapping [34]. The factors influencing debris flow susceptibility were divided into three categories: geomorphological, materials, and triggering. The specific selection parameters and their bases are as follows.

Parameters Related to Geomorphological Conditions
The catchment area (CA) can reflect the sediment yield and concentration within the basin [17]. Channel length (CL) [17,43] reflects the potential length of debris flow movement and the mass of loose materials along its path. Catchment perimeter (CP) reflects the morphological characteristics of the catchment.
Average slope (AS) is the average value of the slope within the boundary of each catchment, which can reflect the overall slope of the catchment and slope stability. Catchment relief (CR) is the height difference between the top and the outlet of the catchment [44]. This can reflect the overall potential energy conditions of the catchment. The relief ratio (RR) reflects the overall steepness of the

Selection of Factors Influencing Debris Flows
The selection and processing of prediction factors have an important influence on the accuracy of debris flow susceptibility mapping [34]. The factors influencing debris flow susceptibility were divided into three categories: geomorphological, materials, and triggering. The specific selection parameters and their bases are as follows.

Parameters Related to Geomorphological Conditions
The catchment area (CA) can reflect the sediment yield and concentration within the basin [17]. Channel length (CL) [17,43] reflects the potential length of debris flow movement and the mass of loose materials along its path. Catchment perimeter (CP) reflects the morphological characteristics of the catchment.
Average slope (AS) is the average value of the slope within the boundary of each catchment, which can reflect the overall slope of the catchment and slope stability. Catchment relief (CR) is the height difference between the top and the outlet of the catchment [44]. This can reflect the overall potential energy conditions of the catchment. The relief ratio (RR) reflects the overall steepness of the channel, obtained by dividing CR by the CL [45]. AS, CR, and RR are important factors influencing debris flows, and they can provide sufficient potential energy conditions for the initiation and process of debris flow [44].
Cut density (CD) is the ratio of the channel length to CA. This reflects the geological and geomorphic conditions in the catchment, such as the weathering resistance and geomorphic development of rocks [46]. Drainage density (DD) is the ratio of RR to CP. Circularity ratio (CR2) is the ratio of CP to CL [47], reflecting the morphological characteristics of the catchment.
Hypsometric Integral (HI) reflects the slope distribution in the catchment [45]; it was calculated using the relief ratio method proposed by Pike and Wilson (1971) [48].
where H mean , H max , and H min are the average, maximum, and minimum elevations in the catchment, respectively. The Melton Ratio (MR) is the ratio of the CR to the square root of the CA. This reflects the dynamics and the susceptibility of the catchment to debris flows, and is widely used in the study of debris flows [16,43,[49][50][51].

Parameters Related to Soil Material and Geology Conditions
Catchment lithology determines the material background conditions for the formation of debris flows. As lithological data cannot be assigned to the catchment unit, in order to simplify the calculation, we used the solution proposed by Zhao et al. [35], which involves dividing the stratigraphic lithological data (vector maps published in December 2002, and supplied by China Geological Survey, 1:500,000) into five levels: very hard, hard, medium, soft, and very soft, according to the rock strength, and assigning values of 1-5, respectively ( Table 1). The vector data of strata were then transformed into grid data using a resolution of 30 × 30 m, so that each grid has an independent lithological strength value. Through the statistical analysis of each catchment, the average lithological strength value of each catchment is obtained as the formation lithology index (FLI) of that catchment (Figure 4a). The higher the FLI value, the softer the overall lithology of the catchment, and the greater the quantity of loose materials that can be provided. . It provides monthly precipitation data in mm/hr. The data product represents the average level of monthly precipitation. The monthly precipitation in 2017 is calculated and accumulated to obtain the annual precipitation in the study area (mm). In summary, a spatial database of catchments in the study area was established ( Table 2). Vegetation coverage is one of the most important parameters for evaluating DFs [28]. NDVI data cannot be assigned according to the catchment unit, and the average NDVI value of each catchment was used as the vegetation coverage index (VCI) of that catchment (Figure 4b). The higher the VCI, the greater the vegetation coverage. NDVI was derived from Landsat-8 images (June 2017) with a resolution of 30-m (Landsat-8 image courtesy of the U.S. Geological Survey).
Faults trigger earthquakes, producing discontinuities in rocks, and as a consequence, they provide rock debris that can be mobilized [24]. Therefore, the distance (km) from fault (DFF) is an important parameter. Based on fault distribution data (vector maps published in December 2002, and supplied by China Geological Survey), the average distance between each catchment and fault was calculated (Figure 4c). The lower the DFFs, the greater the influence of the fault and the higher the degree of fracture of the adjacent rock mass.
Chen et al. [52] found that debris flow occurrence was closely related to the impact of earthquakes and droughts, because they can increase the amount of loose materials. Wu et al. [53] evaluated the debris flow susceptibility of the Longxi River in the Wenchuan earthquake area. They summed the volumes of the source debris accumulations to reflect the impact of earthquakes on the debris flow susceptibility. Therefore, the average distance (km) from previous earthquake epicenters (DFE) is an important parameter. Based on the earthquake distribution data (the earthquake database was provided by USGS, and the sites with magnitude above 4 on the Richter scale since 2000 were selected, https://earthquake.usgs.gov/earthquakes), the average distance from the epicenter to each catchment was calculated (Figure 4d). The lower the DFE, the greater the earthquake influence, and the greater the slope instability.

Parameters Related to Triggering Conditions
Rainfall is the major factor triggering debris flow hazards [30]. In this study, the average annual rainfall of each catchment was calculated as the annual precipitation index (API) (Figure 4e). Annual precipitation was derived from TRMM_3B43 images with a resolution of 0. . It provides monthly precipitation data in mm/hr. The data product represents the average level of monthly precipitation. The monthly precipitation in 2017 is calculated and accumulated to obtain the annual precipitation in the study area (mm).
In summary, a spatial database of catchments in the study area was established ( Table 2). For some models, strongly correlated parameters contain some degree of redundancy, and substantial collinearity will lead to model instability. In this study, a parameter correlation matrix heat map ( Figure 5) was calculated based on the Seaborn Python visualization package. It was found that there were several pairs of variables which were highly correlated; for example, the following correlation coefficients were obtained: 0.92 for CA vs. CP; 0.88 for CL vs. CA; 0.99 for CL vs. CP; 0.92 for RR vs. DD; and 0.91 for RR and MR. Therefore, CL, CP, and RR were eliminated. Remote Sens. 2020, 12, x FOR PEER REVIEW 10 of 22

Resampling
For the sample data collected in this study, the ratio of DFs to NDFs is 5:1 ( Figure 6). This data imbalance will cause the learning machine to pay more attention to the classifications with DFs and hence, obtain an inaccurate model. In this study, SMOTE (the synthetic mineral oversampling technique) [54] was used to increase the number of NDF samples. This method randomly selects a nearest neighbor sample B from A (a sample in NDFs), and then, randomly selects a point C from the relationship between A and B, as a new minority sample. After resampling, the sample size of NDFs was increased to 544, and the ratio of NDFs to DFs was 1:1.

Resampling
For the sample data collected in this study, the ratio of DFs to NDFs is 5:1 ( Figure 6). This data imbalance will cause the learning machine to pay more attention to the classifications with DFs and hence, obtain an inaccurate model. In this study, SMOTE (the synthetic mineral oversampling technique) [54] was used to increase the number of NDF samples. This method randomly selects a nearest neighbor sample B from A (a sample in NDFs), and then, randomly selects a point C from the relationship between A and B, as a new minority sample. After resampling, the sample size of NDFs was increased to 544, and the ratio of NDFs to DFs was 1:1.

Resampling
For the sample data collected in this study, the ratio of DFs to NDFs is 5:1 ( Figure 6). This data imbalance will cause the learning machine to pay more attention to the classifications with DFs and hence, obtain an inaccurate model. In this study, SMOTE (the synthetic mineral oversampling technique) [54] was used to increase the number of NDF samples. This method randomly selects a nearest neighbor sample B from A (a sample in NDFs), and then, randomly selects a point C from the relationship between A and B, as a new minority sample. After resampling, the sample size of NDFs was increased to 544, and the ratio of NDFs to DFs was 1:1.

Data Standardization
Data standardization can speed up model convergence and improve model accuracy. Moreover, some machines are very sensitive to feature scales; therefore, in this study, we used a standard scaler algorithm (from Scikit-learn, https://scikit-learn.org) to standardize the features by removing the mean and scaling to unit variance. Scikit-learn is a Python library that provides a standard interface for implementing machine learning algorithms [55].

Generating a Cross-Validation Dataset
Using the cross-validation algorithm in Scikit-learn, 70% of the training data were randomly selected for model training, and the remaining 30% were used as test data to evaluate the model; this was repeated 10 times (Figure 7). In this way, the algorithm can use different training data subsets to build the model, and test data are used to evaluate model performance of the model and prevent the model from being over fitted.
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 22 Data standardization can speed up model convergence and improve model accuracy. Moreover, some machines are very sensitive to feature scales; therefore, in this study, we used a standard scaler algorithm (from Scikit-learn, https://scikit-learn.org) to standardize the features by removing the mean and scaling to unit variance. Scikit-learn is a Python library that provides a standard interface for implementing machine learning algorithms [55].

Generating a Cross-Validation Dataset
Using the cross-validation algorithm in Scikit-learn, 70% of the training data were randomly selected for model training, and the remaining 30% were used as test data to evaluate the model; this was repeated 10 times (Figure 7). In this way, the algorithm can use different training data subsets to build the model, and test data are used to evaluate model performance of the model and prevent the model from being over fitted.

Candidate Machine Selection
In order to find the most appropriate model to solve the research problem, we selected a variety of widely used machine learning methods. Typical algorithms were selected of each type and 22 model methods were finally selected.

Ensemble Methods
The concept of the Ensemble method is to combine several classifiers (or combine different parameters of an algorithm) to improve the effect of each single classifier. The classifiers can be divided into two types: average methods and boosting methods. AdaBoost, Gradient Tree Boosting (GDBT), Bagging, Random Forest, and Extra Trees were selected in this study.

Gaussian Processes (GP)
Gaussian Processes are a generic supervised learning method designed to solve regression and probabilistic classification problems.

Generalized Linear Models (GLM)
The Generalized Linear model is an extension of the linear model. The relationship between the mathematical expectation of the response variable and the prediction variable of the linear combination is established by the connection function. Logistic Regression (LR), Passive Aggressive, Ridge, Stochastic Gradient Descent (SGD), and Perceptron were selected in this study.

Naive Bayes (NB)
Naive Bayes classification is based on the concept of Bayes probability. Assuming that the attributes are independent of each other, the probability of each feature is obtained and the larger one is taken as the prediction result. Gaussian Naive Bayes and Bernoulli Naive Bayes were selected.

Candidate Machine Selection
In order to find the most appropriate model to solve the research problem, we selected a variety of widely used machine learning methods. Typical algorithms were selected of each type and 22 model methods were finally selected.

Ensemble Methods
The concept of the Ensemble method is to combine several classifiers (or combine different parameters of an algorithm) to improve the effect of each single classifier. The classifiers can be divided into two types: average methods and boosting methods. AdaBoost, Gradient Tree Boosting (GDBT), Bagging, Random Forest, and Extra Trees were selected in this study.

Gaussian Processes (GP)
Gaussian Processes are a generic supervised learning method designed to solve regression and probabilistic classification problems.

Generalized Linear Models (GLM)
The Generalized Linear model is an extension of the linear model. The relationship between the mathematical expectation of the response variable and the prediction variable of the linear combination is established by the connection function. Logistic Regression (LR), Passive Aggressive, Ridge, Stochastic Gradient Descent (SGD), and Perceptron were selected in this study.

Naive Bayes (NB)
Naive Bayes classification is based on the concept of Bayes probability. Assuming that the attributes are independent of each other, the probability of each feature is obtained and the larger one is taken as the prediction result. Gaussian Naive Bayes and Bernoulli Naive Bayes were selected.

Nearest Neighbors
The principle of the Nearest Neighbor method is to find a specified number of nearest sample points and then, use them to predict new points.

Support Vector Machines (SVM)
The basic concept of SVM is to solve the separation hyperplane, which can divide the training dataset correctly and provide the largest geometric interval. Support Vector Classification (SVC), Linear SVC, and Nu-SVC were selected.

Trees
The Tree classifier is a tree structure in which each internal node represents a judgment of an attribute, each branch represents an output of the judgment result, and finally, each leaf node represents a classification result. Decision Tree and Extra Tree were selected.

Discriminant Analysis
Discriminant Analysis is a method of multivariate statistical analysis which classifies the studied objects according to several observed indexes. Linear Discriminant and Quadratic Discriminant were selected.

eXtreme Gradient Boosting (XGBoost)
XGBoost is a boosting algorithm, and it is a type of lifting tree model. It implements the GBDT algorithm efficiently and makes many improvements to the algorithm, integrating numerous tree models to produce a strong classifier.

Model Fitting and Tuning
The training data in the cross-validation dataset were used to train the initial model. The models were then sorted according to the average accuracy score (ACC) of the test data in the cross-validation dataset ( Figure 8). The ACC (Formula 3) indicates the rate of correct assignment of all samples participating in the modeling (Table 3). It is evident from Figure 8 that the overall fitting effect of the integrated model is better than that of the other models. The highest score was achieved by Extra Trees, followed by Random Forest, Gaussian Process, Gradient Boosting, and XGBoost.
Among them (Table 3)  The top 12 models (ACC > 0.75) were selected for model optimization. Parameter grid and grid search cross-validation were used to fit the model, and the AUC (area under the receiver operating characteristic curve) score was used to search for the best super parameter. According to the optimal super parameters of each model given in Table 4, 10-fold cross-validation was again performed on the model training set, and the models were reordered according to the average accuracy score of the test data ( Figure 9). It can be seen that the overall accuracy of the models had been improved. After model optimization, the SVC model became the optimal model; its test data ACC is 0.91, and the average AUC of 10-fold cross-validation is 0.96 ( Figure 10). The AUC indicates the tradeoff between sensitivity and specificity [56]. From the perspective of the time needed to find the optimal parameters, the most time consuming model was Gradient Boosting (112,533.93 s), the least time consuming was Extra Tree (1.53 s), and the time consumption of SVC was 57.79 s. Therefore, SVC provided high efficiency combined with the highest accuracy.   The top 12 models (ACC > 0.75) were selected for model optimization. Parameter grid and grid search cross-validation were used to fit the model, and the AUC (area under the receiver operating characteristic curve) score was used to search for the best super parameter. According to the optimal super parameters of each model given in Table 4, 10-fold cross-validation was again performed on the model training set, and the models were reordered according to the average accuracy score of the test data ( Figure 9). It can be seen that the overall accuracy of the models had been improved. After model optimization, the SVC model became the optimal model; its test data ACC is 0.91, and the average AUC of 10-fold cross-validation is 0.96 ( Figure 10). The AUC indicates the tradeoff between sensitivity and specificity [56]. From the perspective of the time needed to find the optimal parameters, the most time consuming model was Gradient Boosting (112,533.93 s), the least time consuming was Extra Tree (1.53 s), and the time consumption of SVC was 57.79 s. Therefore, SVC provided high efficiency combined with the highest accuracy.

Results
The optimal model SVC (Super parameter: 'C' = 4, 'decision_function_shape' = 'ovo', 'gamma' = 0.5, 'probability' = True) was selected as the final model to map the landslide susceptibility of the study area. The SVC can only output classification results (i.e., 0/1), so the model cannot directly output the probability; however, the probability value is needed to map debris flow susceptibility. We used the method provided by Scikit-learn to set the probability to true when constructing the SVC function option, so that the model could output the classification probability (0-1). In the binary case, the probabilities are calibrated using Platt scaling (Platt "Probabilistic outputs for SVMs and comparisons to regularized likelihood methods"): logistic regression of the SVM's scores and fitting by an additional cross-validation of the training data [57].
The probability value of debris flow susceptibility was divided into five categories using the natural fracture method [58]: very low, low, moderate, high, and very high. The susceptibility mapping results ( Figure 11) show that debris flow susceptibility in large catchments (>100 km 2 ) is generally low, which may be because the fluvial processes take over the debris flow processes in these catchments. The areas of high and very high debris flow susceptibility are evenly distributed in the study area, and there is no obvious spatial distribution pattern. The proportion of debris flow susceptibility ranges from low to high, with values of 16.7%, 7.2%, 41.5%, 26.5%, and 8.1%, respectively. The areas with moderate susceptibility account for the highest proportion, while the areas of low susceptibility and very high susceptibility account for the lowest proportion. In order to help in understanding the magnitude of the problem as well as in early planning of the tasks, the lengths of the part of highway under low to high susceptibility were counted according to the fans interpreted in Google Earth, with lengths of 27 Figure 11. Susceptibility mapping results of the KKH area.

Discussion
The interpretability of the model can help determine the potential relationship between the influencing factors and debris flow susceptibility from the data, which is very important for understanding the factors influencing debris flow susceptibility in the study area. The SVC model is a "black box" process which cannot directly indicate the importance of the variables employed. Here, Figure 11. Susceptibility mapping results of the KKH area.

Discussion
The interpretability of the model can help determine the potential relationship between the influencing factors and debris flow susceptibility from the data, which is very important for understanding the factors influencing debris flow susceptibility in the study area. The SVC model is a "black box" process which cannot directly indicate the importance of the variables employed. Here, the importance of the variables was calculated using the Permutation Importance algorithm provided by the extension module of ELI5 (ELI5 is a Python library, which allows the visualization and debugging of various Machine Learning models using unified API. It has built-in support for several ML frameworks and provides a way to explain black-box models, https://eli5.readthedocs.io/en/latest/index.html). This determines the importance of parameters by measuring how to reduce the score when a feature is not available. This method is also known as "rank importance" or "mean decreased accuracy (MDA)" [59].
From the importance of the calculation results (Table 5), it can be seen that all of the parameters involved in the SVC have some degree of importance. Among them, MR is the most important. MR is a parameter proposed by Melton (1965) [60] to reflect the dynamic characteristics of the catchment and the debris flow susceptibility; it is widely used in the evaluation and classification of debris flows [43,50,51,61]. This is followed by DD, HI, and AS, which are related to the geomorphic conditions. This shows that the parameters related to these factors play an important role in the prediction of debris flow susceptibility in the study area. In the study of Xiong et al. [34], several geomorphological parameters (mean altitude, altitude difference, and groove gradient) were also the most important factors for the prediction of debris flow susceptibility. In the study of Zhang et al. [28], aspect was the most important factor, followed by rainfall, elevation, and slope curvature. In addition, there are three parameters related to material conditions: DFE, DFF, and VCI, followed by other factors. Among the factors related to soil material and geology conditions, the most important factor is the DFE, followed by DFF and VCI. Due to the lack of detailed data of snow cover and glacier melting, the annual rainfall was assumed to be the main trigger factor of debris flow susceptibility. From the parameter importance ranking (Table 5), we can see that the importance of API is not as high as we assumed. This may be due to the lack of data of glaciers, snow, and frozen soil. This is a defect of this paper that we hope to improve in future work, in order to obtain a better debris flow susceptibility map. The distributions of MR, DD, HI, and AS are shown in Figure 12. It can be seen that the proportion of DFs is higher when 0.15 < MR < 0.74, DD < 32, 0.45 < HI < 0.58, and AS < 15, AS > 33, which indicates that a debris flow is more likely to occur in these areas. Among these four parameters, MR has the strongest ability to distinguish DF and NDF. It can be seen from the figure that the differences in the peak distribution of DF and NDF in MR are obvious, but that this is not the case for the difference in the other three parameters. This indicates that MR is an important factor affecting debris flow susceptibility in the study area.

CA
0.0748 ± 0.0080 The distributions of MR, DD, HI, and AS are shown in Figure 12. It can be seen that the proportion of DFs is higher when 0.15 < MR < 0.74, DD < 32, 0.45 < HI < 0.58, and AS < 15, AS > 33, which indicates that a debris flow is more likely to occur in these areas. Among these four parameters, MR has the strongest ability to distinguish DF and NDF. It can be seen from the figure that the differences in the peak distribution of DF and NDF in MR are obvious, but that this is not the case for the difference in the other three parameters. This indicates that MR is an important factor affecting debris flow susceptibility in the study area.
Remote Sens. 2020, 12, x FOR PEER REVIEW 18 of 22 Figure 12. Distribution of MR, DD, HI, and AS. "1" indicates catchments prone to debris flow, and "0" indicates catchments not prone to debris flow.

Conclusions
We have applied a variety of machine learning methods to the mapping of debris flow susceptibility of the Karakoram Highway. Good prediction results were obtained. Comparison of super parameter optimization and model prediction accuracy indicated that the Support Vector Classification (SVC) model had the best performance, with an average test data accuracy of 0.91 and an average AUC of 10-fold cross-validation of 0.96. The search time for the optimal parameters of the model was 57.79 s, which shows that it combined high efficiency with the highest prediction accuracy. Hence, it provides a valuable tool for mapping debris flow susceptibility. MR was the most important predictor in the final model, followed by DD, HI, AS, DFE, DFF, VCI, and other factors. This shows that geomorphic conditions play an important role in the prediction of debris flow susceptibility. Among the factors related to soil material and geology conditions, the most important factor is the DFE, followed by DFF and VCI. The importance of API is not as high as we assumed, which may be due to the lack of data of glaciers, snow, and frozen soil.
The susceptibility mapping results show that debris flow susceptibility in large catchments is generally low, which may be because the fluvial processes take over the debris flow processes in these catchments. The areas of high and very high debris flow susceptibility are evenly distributed within Figure 12. Distribution of MR, DD, HI, and AS. "1" indicates catchments prone to debris flow, and "0" indicates catchments not prone to debris flow.

Conclusions
We have applied a variety of machine learning methods to the mapping of debris flow susceptibility of the Karakoram Highway. Good prediction results were obtained. Comparison of super parameter optimization and model prediction accuracy indicated that the Support Vector Classification (SVC) model had the best performance, with an average test data accuracy of 0.91 and an average AUC of 10-fold cross-validation of 0.96. The search time for the optimal parameters of the model was 57.79 s, which shows that it combined high efficiency with the highest prediction accuracy. Hence, it provides a valuable tool for mapping debris flow susceptibility. MR was the most important predictor in the final model, followed by DD, HI, AS, DFE, DFF, VCI, and other factors. This shows that geomorphic conditions play an important role in the prediction of debris flow susceptibility. Among the factors related to soil material and geology conditions, the most important factor is the DFE, followed by DFF and VCI. The importance of API is not as high as we assumed, which may be due to the lack of data of glaciers, snow, and frozen soil.
The susceptibility mapping results show that debris flow susceptibility in large catchments is generally low, which may be because the fluvial processes take over the debris flow processes in these catchments. The areas of high and very high debris flow susceptibility are evenly distributed within the study area, and there is no obvious spatial distribution pattern. The area of moderate debris flow susceptibility accounts for the highest proportion, while the areas of low and very high susceptibility account for the lowest proportion. Overall, we suggest that our approach can provide valuable inputs to policies, promoting the safe operation of highways.