Landslide Susceptibility Evaluation and Management Using Different Machine Learning Methods in The Gallicash River Watershed, Iran

This analysis aims to generate landslide susceptibility maps (LSMs) using various machine learning methods, namely random forest (RF), alternative decision tree (ADTree) and Fisher’s Linear Discriminant Function (FLDA). The results of the FLDA, RF and ADTree models were compared with regard to their applicability for creating an LSM of the Gallicash river watershed in the northern part of Iran close to the Caspian Sea. A landslide inventory map was created using GPS points obtained in a field analysis, high-resolution satellite images, topographic maps and historical records. A total of 249 landslide sites have been identified to date and were used in this study to model and validate the LSMs of the study region. Of the 249 landslide locations, 70% were used as training data and 30% for the validation of the resulting LSMs. Sixteen factors related to topographical, hydrological, soil type, geological and environmental conditions were used and a multi-collinearity test of the landslide conditioning factors (LCFs) was performed. Using the natural break method (NBM) in a geographic information system (GIS), the LSMs generated by the RF, FLDA, and ADTree models were categorized into five classes, namely very low, low, medium, high and very high landslide susceptibility (LS) zones. The very high susceptibility zones cover 15.37% (ADTree), 16.10% (FLDA) and 11.36% (RF) of the total catchment area. The results of the different models (FLDA, RF, and ADTree) were explained and compared using the area under receiver operating characteristics (AUROC) curve, seed cell area index (SCAI), efficiency and true skill statistic (TSS). The accuracy of models was calculated considering both the training and validation data. The results revealed that the AUROC success rates are 0.89 (ADTree), 0.92 (FLDA) and 0.97 (RF) and predication rates are 0.82 (ADTree), 0.79 (FLDA) and 0.98 (RF), which justifies the approach and indicates a reasonably good landslide prediction. The results of the SCAI, efficiency and TSS methods showed that all models have an excellent modeling capability. In a comparison of the models, the RF model outperforms the boosted regression tree (BRT) and ADTree models. The results of the landslide susceptibility modeling could be useful for land-use planning and decision-makers, for managing and controlling the current and future landslides, as well as for the protection of society and the ecosystem.

. Location of the study area.

Methodology
The current research demonstrates the application of different machine learning ensemble techniques in Landslide susceptibility mapping. Preparing the Landslide susceptibility map in the present study involved data collection, inventory map preparation, LCFs preparation, collinearity test, factors importance analysis using the RF and BRT models, LS mapping using RF, ADtree, and FLDA and validation of the models (Figure 2). Both primary and secondary data were used in this work. The LIM was prepared based on the field investigation in a different time (2018/05/12; 2018/08/10; and 2018/10/14) and Google Earth imagery. The topographic parameters of altitude, slope aspect convergence index (CI), stream power index (SPI) and topographical positioning index (TPI) were obtained from the 12.5 m × 12.5 m resolution ALOS PALSAR DEM that was downloaded from the Alaska Satellite Facility homepage. The enhanced thematic mapper plus (ETM+) satellite imagery of 15 m × 15 m resolution was obtained from USGS. The topographic map (1:50,000) and geological map (1:1,000,000) were obtained from the Cartographic and Geological Department, Iran. Rainfall data of different rain gauge stations was also collected from the Metrological Department of Iran. The LSMs were classified into five classes using the natural break method of Jenks. The Jenks method of optimization, also called the natural breaks classification method, is a method of data segmentation designed to determine the best value arrangement in different classes. This method of classification aims to decrease the average deviation from the mean class value while increasing the deviation from the mean of other classes. The method reduces the intra-class variance and maximizes the inter-class

Methodology
The current research demonstrates the application of different machine learning ensemble techniques in Landslide susceptibility mapping. Preparing the Landslide susceptibility map in the present study involved data collection, inventory map preparation, LCFs preparation, collinearity test, factors importance analysis using the RF and BRT models, LS mapping using RF, ADtree, and FLDA and validation of the models ( Figure 2). Both primary and secondary data were used in this work. The LIM was prepared based on the field investigation in a different time (2018/05/12; 2018/08/10; and 2018/10/14) and Google Earth imagery. The topographic parameters of altitude, slope aspect convergence index (CI), stream power index (SPI) and topographical positioning index (TPI) were obtained from the 12.5 m × 12.5 m resolution ALOS PALSAR DEM that was downloaded from the Alaska Satellite Facility homepage. The enhanced thematic mapper plus (ETM+) satellite imagery of 15 m × 15 m resolution was obtained from USGS. The topographic map (1:50,000) and geological map (1:1,000,000) were obtained from the Cartographic and Geological Department, Iran. Rainfall data of different rain gauge stations was also collected from the Metrological Department of Iran. The LSMs were classified into five classes using the natural break method of Jenks. The Jenks method of optimization, also called the natural breaks classification method, is a method of data segmentation designed to determine the best value arrangement in different classes. This method of classification aims to decrease the average deviation from the mean class value while increasing the deviation from the mean of other classes. The method reduces the intra-class variance and maximizes the inter-class variance [54,55]. The natural breaks classification method is an iterative process. That is, it repeats calculations using different breaks in the dataset to decide which set of the breaks have the smallest Remote Sens. 2020, 12, 475 5 of 29 variance in the class. There are three steps that need to be followed for this method: (1) calculate the sum of squared deviations from the class mean (SDCM); (2) calculate the sum of squared deviations from the array mean (SDAM); (3) after each SDCM is checked, decided to move one unit from a class with a larger SDCM to an adjacent class with a lower SDCM. New class deviations are then measured, and the process continues until a minimum value is reached by the sum of the deviations within the class [54]. The whole process was carried out using R studio, ArcGIS, ENVI and Weka.
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 30 variance [54,55]. The natural breaks classification method is an iterative process. That is, it repeats calculations using different breaks in the dataset to decide which set of the breaks have the smallest variance in the class. There are three steps that need to be followed for this method: (1) calculate the sum of squared deviations from the class mean (SDCM); (2) calculate the sum of squared deviations from the array mean (SDAM); (3) after each SDCM is checked, decided to move one unit from a class with a larger SDCM to an adjacent class with a lower SDCM. New class deviations are then measured, and the process continues until a minimum value is reached by the sum of the deviations within the class [54]. The whole process was carried out using R studio, ArcGIS, ENVI and Weka.

LIM
The LIM can be defined as the total set of landslide locations as it includes both past and present landslide points. The LIM is an essential part of the landslide susceptibility mapping. The relationship between the landslide inventory datasets and landslide conditioning factors (LCFs) is the most important for landslide susceptibility mapping [56]. The LIM contains both the old and new landslides that have occurred in the study area. A LIM can be prepared based on a combination of several data sources, including previous records, local field examinations, the perception survey with inhabitants and interpretation of satellite images [57]. Some researchers have used geomorphologic features visible in satellite imagery to detect landslides [58].
In this study, a LIM was prepared using historical records [59,60], aerial photos, field inquiry with GPS, Google Earth imagery and the PALSAR DEM. The PALSAR DEM is more accurate for detecting landslides in forested areas [61]. The resolution of the PALSAR DEM is 12.5 m × 12.5 m. Many landslides, which due to their remote locations aren't even detectable through field surveys, were extracted with acceptable resolution using Google Earth's optical images. A total of 249 landslides were randomly selected in the basin. Of the 249 landslides, 70% were randomly selected to be used for training, and 30% were used for validation. The research area's landslide identification results show that landslides in the study area are mainly translational earth flow slides, rotational slides and debris flows (Figure 3a). Of the 249 identified landslides 75 (30%) are creeps, 47 (19%) are rotational slides, 27 (11%) are earth flows, 87 (35%) are translational slides, and 13 (5%) are debris flows. The combined area covered by all the landslides in the study area is 16,213,651 m 2 . Of the identified

LIM
The LIM can be defined as the total set of landslide locations as it includes both past and present landslide points. The LIM is an essential part of the landslide susceptibility mapping. The relationship between the landslide inventory datasets and landslide conditioning factors (LCFs) is the most important for landslide susceptibility mapping [56]. The LIM contains both the old and new landslides that have occurred in the study area. A LIM can be prepared based on a combination of several data sources, including previous records, local field examinations, the perception survey with inhabitants and interpretation of satellite images [57]. Some researchers have used geomorphologic features visible in satellite imagery to detect landslides [58].
In this study, a LIM was prepared using historical records [59,60], aerial photos, field inquiry with GPS, Google Earth imagery and the PALSAR DEM. The PALSAR DEM is more accurate for detecting landslides in forested areas [61]. The resolution of the PALSAR DEM is 12.5 m × 12.5 m. Many landslides, which due to their remote locations aren't even detectable through field surveys, were extracted with acceptable resolution using Google Earth's optical images. A total of 249 landslides were randomly selected in the basin. Of the 249 landslides, 70% were randomly selected to be used for training, and 30% were used for validation. The research area's landslide identification results show that landslides in the study area are mainly translational earth flow slides, rotational slides and debris flows (Figure 3a). Of the 249 identified landslides 75 (30%) are creeps, 47 (19%)   Very little research is available on analyzing areas where the chance of landslides occurring is virtually zero, which is the definition of non-landslide areas [62]. Marchesini et al. [63] recommended a non-linear quantile model for identifying non-susceptible landslide areas based on a morphometric evaluation. However, non-landslide areas were delimited in the current study based on the availability of good inventory maps. It was necessary to generate non-landslide evidence for modeling as the landslide susceptibility analysis is based on a binary classification [64]. In the present research, we also used non-landslide data for generating the training and validation datasets by identifying non-landslide areas on Google Earth images where no building and slope failures were found. For both training and validation datasets, the number of non-landslide data equals the number of landslide data. Landslide data was assigned a value of "1," while non-landslide data was assigned a value of "0" for spatial landslide susceptibility modeling. Figure 1 shows the landslide locations of the training and validation datasets. The percentages of different landslide types along with some field photographs, are shown in Figure 3a-e. Very little research is available on analyzing areas where the chance of landslides occurring is virtually zero, which is the definition of non-landslide areas [62]. Marchesini et al. [63] recommended a non-linear quantile model for identifying non-susceptible landslide areas based on a morphometric evaluation. However, non-landslide areas were delimited in the current study based on the availability of good inventory maps. It was necessary to generate non-landslide evidence for modeling as the landslide susceptibility analysis is based on a binary classification [64]. In the present research, we also used non-landslide data for generating the training and validation datasets by identifying non-landslide areas on Google Earth images where no building and slope failures were found. For both training and validation datasets, the number of non-landslide data equals the number of landslide data. Landslide data was assigned a value of "1," while non-landslide data was assigned a value of "0" for spatial landslide susceptibility modeling. Figure 1 shows the landslide locations of the training and validation datasets. The percentages of different landslide types along with some field photographs, are shown in Figure 3a-e.

Preparing Landslide Conditioning Factors (LCFs)
A landslide susceptibility map cannot be accurate without understanding the effective LCFs. The topographical factors, i.e., elevation, slope, slope length (LS), plan curvature, convergence index (CI), stream power index (SPI), topographic position index (TPI), geomorphologic, i.e., distance to stream, topographic wetness index (TWI), geological, i.e., lithology, distance to fault, and environmental, i.e., land-use/land cover (LU/LC), soil type, normalized differential vegetation index (NDVI) were considered based on existing literature [19][20][21][22][23][24][25] and multi-collinearity analysis for the prediction of the landslide risk zones in the Gallicash river watershed. Determining the LCFs is the first step in mapping the LS. The Phased Array type L-band synthetic aperture radar (PALSAR) DEM is more reliable and accurate than the ASTER and SRTM DEM. So, in this study, we used the PALSAR DEM with 12.  Table 1. The geological map of the study area was prepared from the 1:100,000 geological map of Iran, (Figure 4m). The eight geological units in this watershed were identified (Table 1). Faults were extracted from the ETM+2002 satellite imagery using ENVI 4.7 software. The distance to the fault map was prepared using the Euclidian distance tool in ArcGIS (Figure 4i). The Sentinal-2 images with 30 m × 30 m resolution were used to construct the LULC ( Figure 4n) and NDVI maps (Figure 4o) using the supervised classification (maximum likelihood). For the assessment of the classification accuracy of the LULC map, 523 ground control points were used. The overall accuracy of the LULC map is 97.1% (Kappa coefficient = 0.971), which is acceptable. Different types of land-use, namely Afforest, Agriculture, Orchard, Dry farming, Water, Agriculture-Orchard, Rock and Urban area ( Figure 4n) have been identified in this watershed. The soil map was prepared from a projected soil map (scale 1: 50,000), which was obtained from the Agriculture and Natural Resources Center using the digitized process in the GIS environment ( Figure 4p). The study area is composed of different soil types, namely Rock Outcrops/Entisols, Inceptisols, Mollisols, and Alfisols ( Figure 4p and Table 1). The elevation, slope, LS, convergence index (CI), SPI, TPI, TWI were classified into five sub-layers with the aid of the natural break method (NBM) of Jenk's in GIS, as shown in Figure 4 and Table 1. The NDVI was categorized into the three sub-layers of <−0.201, −0.201 −0.36 and >0.36 using the NBM (Figure 4o and Table 1). The soil type, lithology and LU/LC are the categorical variables, with detailed description given in Figure 4 and Table 1.
Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 30 A landslide susceptibility map cannot be accurate without understanding the effective LCFs. The topographical factors, i.e., elevation, slope, slope length (LS), plan curvature, convergence index (CI), stream power index (SPI), topographic position index (TPI), geomorphologic, i.e., distance to stream, topographic wetness index (TWI), geological, i.e., lithology, distance to fault, and environmental, i.e., land-use/land cover (LU/LC), soil type, normalized differential vegetation index (NDVI) were considered based on existing literature [19][20][21][22][23][24][25] and multi-collinearity analysis for the prediction of the landslide risk zones in the Gallicash river watershed. Determining the LCFs is the first step in mapping the LS. The Phased Array type L-band synthetic aperture radar (PALSAR) DEM is more reliable and accurate than the ASTER and SRTM DEM. So, in this study, we used the PALSAR DEM with 12.  Table 1. The geological map of the study area was prepared from the 1:100,000 geological map of Iran, ( Figure  4m). The eight geological units in this watershed were identified (Table 1). Faults were extracted from the ETM+2002 satellite imagery using ENVI 4.7 software. The distance to the fault map was prepared using the Euclidian distance tool in ArcGIS (Figure 4i). The Sentinal-2 images with 30 m × 30 m resolution were used to construct the LULC ( Figure 4n) and NDVI maps (Figure 4o) using the supervised classification (maximum likelihood). For the assessment of the classification accuracy of the LULC map, 523 ground control points were used. The overall accuracy of the LULC map is 97.1% (Kappa coefficient = 0.971), which is acceptable. Different types of land-use, namely Afforest, Agriculture, Orchard, Dry farming, Water, Agriculture-Orchard, Rock and Urban area (Figure 4n) have been identified in this watershed. The soil map was prepared from a projected soil map (scale 1: 50,000), which was obtained from the Agriculture and Natural Resources Center using the digitized process in the GIS environment ( Figure 4p). The study area is composed of different soil types, namely Rock Outcrops/Entisols, Inceptisols, Mollisols, and Alfisols ( Figure 4p and Table 1). The elevation, slope, LS, convergence index (CI), SPI, TPI, TWI were classified into five sub-layers with the aid of the natural break method (NBM) of Jenk's in GIS, as shown in Figure 4 and Table 1. The NDVI was categorized into the three sub-layers of <−0.201, −0.201 −0.36 and >0.36 using the NBM (Figure 4o and Table 1). The soil type, lithology and LU/LC are the categorical variables, with detailed description given in Figure 4 and Table 1.
indicates the average angle between the aspect of adjacent cells and the direction to the central cell.

Testing Multi-Collinearity Problems
The multi-collinearity test may enhance the results of the models by helping in the selection of ideal factors for hazard mapping [77]. Multi-collinearity is the linear dependency that implies that there are two or more related variables in a dataset. The multi-collinearity test has been used for several purposes, such as landslide susceptibility mapping, soil and gully erosion susceptibility, groundwater potentiality mapping, etc. [1,[20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38]. The multi-collinearity was tested by applying the variance inflation factor (VIF) and tolerance [78]. The VIF is reciprocal to the tolerance values, calculated using the following Equations (1) and (2): where R 2 J is the coefficient of determination of a regression of explanator j on all the other explanators. The threshold value of tolerance and VIF is more than 0.1 and less than 10 [79]. Arabameri et al. [1] and Roy and Saha [33] used the multi-collinearity analysis to test the topographical, geomorphological, and environmental factors for LS mapping. In the current research, the multi-collinearity of the sixteen LCFs was tested using the SPSS software.

Applying Random Forest (RF)
RF is one of the main machine learning approaches [80]. The RF method generates multiple classification trees and aggregates them to compute a classification [81]. Hansen and Salamon [82] point out that ensemble classification trees are much more precise than any of their individual members. RF adds diversity among the classification trees by changing the data with alternate and modifying the explanatory factor sets arbitrarily over the different processes of tree induction. The two user-defined variables that are needed to grow a random forest are the number of trees (k) and the number of predictive factors used to divide the nodes (m). The numerical and categorical variables are used for conditioning factors. The out-of-bag (OOB) error is characterized as the percentage of the total number of out-of-bag items that are misclassified. The OOB error is a rational generalization error estimation. The OOB error is estimated at the time of construction of the RF model. Breiman [80] mentioned that the RF produces a limiting value of the error of generalization. The generalization errors often fall together as the number of trees grows. The k has to be set high enough to allow this convergence. The RF method calculates the value of a predictor variable by examining how much the OOB error decreases as OOB data is permuted for that variable while being constant for all other variables. The rise in the error of OOB corresponds to the value of the explanatory variable [83]. One of the major advantages of RF is the resistance to overtraining and the development of a large number of random forest trees where it doesn't really pose a risk of overfitting (e.g., each tree is a random experiment completed independently). There is no need to rescale, transform, or change the RF algorithm. For predictors, RF has resistance to outliers and manages the missing values automatically [76]. In this work, the RF model was computed using the "Randomforest" package in R 3.8, whereby, after a primary discursion, the number of trees was set to 1000 and the m samples at each node was chosen to be 3 to calculate the mutual contribution of attribute subsets while preserving gradual convergence during iterations. No training set is required to monitor the parameters [84,85]. The mean decrease in Gini error and accuracy were analyzed.

Applying Alternating Decision Tree (ADTree)
ADTree is one of the machine learning ensemble classifiers that incorporates the decision tree and thus improves the algorithm [49]. ADTree is related to direct decision tree generalizations in which each node is separated into two nodes, including a node of prediction and a node of splitter. For a mapping of the base ruler, the actual number consists of two conditions, namely a precondition (c 1 ) and a base condition (c 2 ). ADTree was detected the two real numbers a and b. The prediction is if c 1 ∩ c 2 or b if c 1 ∩ −c 2. − identified as negation (NOT). The two real numbers a and b are derived using Equations (3) and (4): where W(p) is the total weight of the p predicate training instances. The best c 1 and c 2 are calculated by minimizing the Z t (c 1 ,c 2 ) calculated using Equation (5): When R symbol is believed to be a set of basic rules, a new base rule will be defined as R t+1 = R t + r t .r t (x) showing the two predictive values of a and b at each predictive layer of T. X is a set of instances. An instance classification is the symbol of the sum of all predictive values in R t+1 : The algorithm begins by choosing the best constant forecasting for the whole dataset [49]. Typically, cross-validation is useful to build the selection [86].

Applying Fisher's Linear Discrimination Analysis (FLDA)
FLDA is another important machine learning approach. FLDA is a feature recognition method used in several fields [87]. Theoretically, m is used for the method as the recognize classes. X (i) j indicates the j-th class i training trial. X i corresponds to the mean of class i training trials, while X refers to complete training trials. The supposed learning trials M b and M w show a scatter matrix for class and with-class and can be determined using Equations (7) and (8) as well: where N i indicates the number of training trials in class i( m i=1 S i = S) and N is the total number of training trials carried out. FLDA aims at obtaining an optimal set of defining vectors for computing Φ d = (ϕ 1 , ϕ 2 , . . . . . . . . ϕ d ) transformation by optimizing the Fisher criterion, as indicated [88]: where T is the matrix's transpose. Because benchmark maximization allows the optimal discriminating vectors, a vector can be derived for one input instance, and the following vector can be used to categorize the subsequent instance [89].

Boosted Regression Tree (BRT) Model
BRT is the non-parametric machine learning method which establishes a relationship between dependent and independent variables [90]. It is a crucial means of measuring the importance of individual variables and to overcome grouping and forecasting problems [91]. The response to a parameter depends on the higher node values in the regression tree's hierarchical structure, and then the relationships between the factors are automatically built [92]. At first, the divisor variable is used to split the target variable in the BRT modeling. To select the best divisor variable, the data was classified into two groups [93]. The recursive aspect point segmentation splits each individual node into two nodes and ultimately separates these divided nodes. This differentiation ends when nodes become homogeneous [91]. The point of the regression analysis is to make a decision during postpone to the divisions. The tree algorithm splits all the data and can lead to overfitting if the divisions do not stop [94]. Trees are more complex and fragile, i.e., they have more errors in the new observations. BRT is calculated using the following Equation (10) [95]: where h(x; m) is described as a basic function of classification with parameter a and variables x, m is the stage of the model. BRT was computed using R 3.3.1 and the "Rattle" package [96].

Application of Frequency Ratio Model
The frequency ratio (FR) is a geospatial assessment tool and one of the important bivariate statistical methods [97]. The FR values are determined from the relationship between landslide locations and LCFs. The result of FR is useful and appropriate [98]. The FR can be described by Equation (11) [99]: where A is the number of landslide pixels in each class of LCFs, B is the total number of landslide pixels in the watershed, C is the total number of pixels of each sub-class of the LCFs, and D is the total number of pixels of the watershed.

Methods for Validating the Models
Validation is necessary for the justification of the LSMs created using different models. Without validation, the LSMs have no value or meaning [100]. The area under receiver operating characteristics (AUROC) curve was used to test the various models. The AUROC curve is a threshold-independent tool for the measurement of predictive performance [101][102][103][104][105][106]. The AUROC indicates the model's predictive accuracy. The AUROC is expressed as the percentage of prediction values [107]. Generally, the value of AUROC ranges from 0.5 to 1. The closer an AUROC value is to 1, the better the performance of the method [9,10]. The AUROC values were classified by Yeashiar and Topal [108] into five classes: low (0.5-0.6), moderate (0.6-0.7), good (0.7-0.8), very good (0.8-0.9) and excellent (0.9-1). Using all the learning and testing datasets, this approach was used to construct a curve, and subsequently, it was used to calculate the goodness of fit and model prediction accuracy [109]. The goodness of fit is superior in measuring the predictive performance of the models using the training and testing datasets [89]. In this work, the AUROC curve was constructed based on the training and testing landslide datasets. The AUROC curve plot consists of two variables, namely the true positive rate (TPR) and the false positive rate (FPR), called the sensitivity and specificity [110]. The TPR shows the percentage of landslide pixels accurately categorized as landslide occurrences, and the FPR indicates the percentage of non-landslide pixels wrongly identified as landslides [111]. The Formulas (12) and (13) were used to calculate the TPR and FPR: where TP (true positive) and TN (true negative) are accurately classified pixel numbers, while FP (false positive) and FN (false negative) are falsely classified pixel numbers. The AUROC reflects the reliability of simulations in correctly predicting landslide occurrence or non-occurrence [108]. The SCAI validation technique used in this study was developed by Süzen and Doyuran [112]. The SCAI is defined as the ratio between the percentage of pixels of the specific landslide susceptibility class and the percentage of the pixels of existing landslides in this specific landslide susceptibility zone. Two threshold-dependent measurement metrics, namely the efficiency (E) and true skill statistic (TSS), were also used to analyze the performance of the models [113]. The E and TSS can be calculated using Equations (14) and (15):

Considering Multi-Collinearity of Factors Contributing to Landslide Susceptibility
A multi-collinearity test was performed on the LCFs. Chen et al. [14] mentioned that the linear relationship between the parameters would minimize the model's predictability. The tolerance and variance inflation factors (VIF) were utilized to test the multi-collinearity. Tolerance and VIF values of <0.1 and above 10 indicate multi-collinearity problems. Arabameri et al. [1] and Roy and Saha [34] used Tolerance and VIF techniques to select suitable parameters for LSMs. In this study, different LCFs, namely elevation, slope, slope length, plan curvature, rainfall, CI, TPI, TWI, SPI, LULC, NDVI, lithology, soil type, distance from stream, distance from fault and distance to road were selected using the multi-collinearity test for the creation of LSMs using different models. The outcomes of the multi-collinearity analysis in this study ( Table 2) show that the effective factors of landslides have no multi-collinearity problem because all the LCFs have lower than threshold values of tolerance and VIF. The maximum tolerance value was identified for TPI (0.955), and the VIF value for TWI (4.55) was identified from the conditioning factors. Therefore, all the proposed landslide conditioning factors were considered to be suitable to produce the LSMs using RF, ADTree and FLDA models in Gallicash catchment area.

The Spatial Relationship Between Landslide Locations and Effective Factors by FR
The frequency ratio (FR) values were calculated based on the relationship between locations of landslides and the sixteen LCFs. The mathematical calculation of the FR value is presented in Table 3. An FR value over 1 indicates a high landslide susceptibility, whereas an FR value of 0 or close to 0 indicates that the sub-class of landslide conditioning factors does not affect landslide occurrence. The FR model was applied to the landslide conditioning factors using Equation (11).  The topographic factors are of great importance for the preparation of LSM [98]. The elevation has a strong relationship with landslides occurrences. The sub-classes <354 m of elevation dataset has the highest FR values (FR = 1.35), indicating the zone of highest landslide susceptibility. Precipitation is the primary cause of landslide occurrences of a region. The long duration of heavy precipitation directly affects landslides due to the loosening of soils, increase of surface runoff velocity and accelerated rate of sediment transmission. The FR results showed that the sub-layer of >851 mm, which is the highest rainfall in the study area, has the maximum FR value of 1.93, followed by the sub-classes of <425 mm (0.08), 425-561 mm (1.72) and 561-703 mm (0.45). In the case of the slope, the steeply sloping area is the most susceptible to landslide occurrences in comparison to gently sloping areas. In this study, The convergence index (CI) is the topographic factor that plays a role in the landslides occurrence of a region. The values of CI vary from +100 to −100. The sub-layers of <−32.54 in the present study have a maximum FR value of 1.39, preceded by −32.54 to −9.8 (FR = 1.15) and −9.8 to 9.1 (FR = 1.0), respectively. The regions that have negative values of CI are more prone to landslides than the areas that have positive CI values. The correlation between landslide locations and road distance indicates that the subset <500 m has the maximum FR value (2.24), which indicates an extreme landslide risk. Landslides frequently occur close to roads in this region. The distance sub-classes of 1000−1500 m, 1500-2000 m and >2000 m indicate a weaker relationship to landslide occurrence with greater distance to the roads. In terms of the distance to streams, close proximity to a stream has a strong control on landslide occurrence. The FR results show that the sub-layer of <110 m, which is the least distance from a stream, has a direct impact on the occurrence of landslides due to the high velocity and discharge of the streams. Comparatively, a greater distance to streams does not directly affect the occurrence of landslides. Eight LULC classes were identified in the study area, namely afforested area, agriculture land, orchard, dry farming, water, agriculture-orchard, rock outcrops. The majority of landslides occurred in the mountain rocks area (FR = 6.95) followed by agriculture-orchard (FR = 2.11), urban areas (FR = 1.87), the forest (FR = 0.47), agriculture (FR = 0.82) and dry farming (FR = 0.96) areas. A high vegetation density can reduce the landslide risk and also postpone the soil transmission and erosion, while a high concentration of landslides can be found in the low vegetation area, which is more susceptible to the occurrence of landslides. In this study, the sub-class of −0.201 to 0.36 has the maximum FR of 1.64, suggesting susceptibility to landslides. Therefore, the FR model depicts the reciprocal relationship of the landslide locations and landslide conditioning factors, taken as the research purpose to manage the landslide mitigation. The pixel-based study using FR is more accurate and appropriate for landslide susceptibility mapping.

Landslide Susceptibility Models
In this landslide susceptibility assessment, the dependent factors were expressed as binary variables, i.e., landslides and non-landslides. The numerical and categorical values of landslide conditioning factors were extracted for creating the LSMs based on the models using the program R. The landslide susceptibility mapping based on the random forest (RF) was constructed using the relative weight of the mean decrees in accuracy and the mean decrees in accuracy of the Gini index of LCFs (Table 7). The RF model results show that the very low, low, moderate, high and very high landslide susceptibility zones of the LSM cover 25.75%, 28.61%, 19.46%, 14.82% and 11.36% areas of the watershed, respectively ( Figure 5). Only 11.36% of the basin, mainly in the middle portion along the north-south line, has an even higher risk of landslides. Likewise, in the GIS setting, the LSM generated using the FLDA model was categorized into five sub-layers ( Figure 6). The results of the LSM based on the FLDA model show that 16.10% and 26.07% of the area have very high and high landslide susceptibility (LS), while 10.64%, 22.20%, and 24.99% of the area have very low, low, and moderate LS. The ADTree LS map was prepared by summarizing all prediction node values in the GIS environment (Figure 7). In this study, the ADTree model was classified into the five classes of very low, low, medium, high and very high according to the NBM of Jenk's. From total of study area, 31.27% has very low susceptibility to landslide occurance), 37.81% has low, 1.53% has moderate, 14.03 % has high and 15.37% has very high susceptibility to gully co-occurance.
conditioning factors were extracted for creating the LSMs based on the models using the program R. The landslide susceptibility mapping based on the random forest (RF) was constructed using the relative weight of the mean decrees in accuracy and the mean decrees in accuracy of the Gini index of LCFs (Table 7). The RF model results show that the very low, low, moderate, high and very high landslide susceptibility zones of the LSM cover 25.75%, 28.61%, 19.46%, 14.82% and 11.36% areas of the watershed, respectively ( Figure 5). Only 11.36% of the basin, mainly in the middle portion along the north-south line, has an even higher risk of landslides. Likewise, in the GIS setting, the LSM generated using the FLDA model was categorized into five sub-layers ( Figure 6). The results of the LSM based on the FLDA model show that 16.10% and 26.07% of the area have very high and high landslide susceptibility (LS), while 10.64%, 22.20%, and 24.99% of the area have very low, low, and moderate LS. The ADTree LS map was prepared by summarizing all prediction node values in the GIS environment (Figure 7). In this study, the ADTree model was classified into the five classes of very low, low, medium, high and very high according to the NBM of Jenk's. From total of study area, 31.27% has very low susceptibility to landslide occurance), 37.81% has low, 1.53% has moderate, 14.03 % has high and 15.37% has very high susceptibility to gully co-occurance.

Validation of Machine Learning Models
The success and predication curves ( Figure 8) were generated using the training and testing data sets by AUROC curves. The AUROC values of the success rate curves (SRC) are 0.896 (89%) for ADTree, 0.927 (92%) for FLDA and 0.972 (97%) for RF models, which means that the LSMs are very good to excellent. The AUC values of the prediction rate curve (PRC) are 0.822 (82%) for ADTree, 0.795 (79%) for FLDA and 0.985 (98%) for RF models, which is also very good to excellent for the resulting LSMs (Figure 8). Very similar results have been found for the success and prediction curves. Of the selected machine learning ensemble methods, the RF model produces a more accurate LSM.

Validation of Machine Learning Models
The success and predication curves ( Figure 8) were generated using the training and testing data sets by AUROC curves. The AUROC values of the success rate curves (SRC) are 0.896 (89%) for ADTree, 0.927 (92%) for FLDA and 0.972 (97%) for RF models, which means that the LSMs are very good to excellent. The AUC values of the prediction rate curve (PRC) are 0.822 (82%) for ADTree, 0.795 (79%) for FLDA and 0.985 (98%) for RF models, which is also very good to excellent for the resulting LSMs ( Figure 8). Very similar results have been found for the success and prediction curves. Of the selected machine learning ensemble methods, the RF model produces a more accurate LSM. The seed cell area index (SCAI) is a reliable validation technique. The validation results of the LSMs using the SCAI method are shown in Table 4. SCAI is the ratio of the percentage of each sublayer to the percentage of landslides that occur in the class [52]. If values of the SCAI, decreased from very low to very high classes of LS, the model is regarded as excellent [1,. The LSMs produced based on the ADTree, FLDA and RF models were categorized into five classes. These landslide susceptibility classes and the number of pixels of training and validation datasets were computed using the SCAI statistical method. The SCAI values of very high landslides susceptibility classes are 0.28, 0.34 and 0.14 for the FLDA, ADTree, and RF models, respectively ( Table 4).
The efficiency values of the FLDA, ADTree, and RF models using the testing dataset are 0.69, 0.77 and 0.69, respectively. The efficiency values of the FLDA, ADTree, and RF models using a training dataset are 0.74, 0.77 and 0.69, respectively. The TSS values of the FLDA, ADTree, and RF models using the testing dataset are 0.37, 0.55 and 0.39, respectively ( Table 5). The TSS values of the FLDA, ADTree, and RF models using the training dataset are 0.47, 0.53 and 0.38, respectively. The AUROC, SCAI, efficiency and TSS methods are a good representation of reality. The outcome of the validation also indicates a strong relationship between the current landslide locations (validation data set) and the model's expected landslide locations. Above all, as per the validation methods used here, the RF machine learning technique performs better than the other models for creating LSMs in the Gallicash watershed in Iran.   The seed cell area index (SCAI) is a reliable validation technique. The validation results of the LSMs using the SCAI method are shown in Table 4. SCAI is the ratio of the percentage of each sub-layer to the percentage of landslides that occur in the class [52]. If values of the SCAI, decreased from very low to very high classes of LS, the model is regarded as excellent [1,. The LSMs produced based on the ADTree, FLDA and RF models were categorized into five classes. These landslide susceptibility classes and the number of pixels of training and validation datasets were computed using the SCAI statistical method. The SCAI values of very high landslides susceptibility classes are 0.28, 0.34 and 0.14 for the FLDA, ADTree, and RF models, respectively (Table 4). The efficiency values of the FLDA, ADTree, and RF models using the testing dataset are 0.69, 0.77 and 0.69, respectively. The efficiency values of the FLDA, ADTree, and RF models using a training dataset are 0.74, 0.77 and 0.69, respectively. The TSS values of the FLDA, ADTree, and RF models using the testing dataset are 0.37, 0.55 and 0.39, respectively ( Table 5). The TSS values of the FLDA, ADTree, and RF models using the training dataset are 0.47, 0.53 and 0.38, respectively. The AUROC, SCAI, efficiency and TSS methods are a good representation of reality. The outcome of the validation also indicates a strong relationship between the current landslide locations (validation data set) and the model's expected landslide locations. Above all, as per the validation methods used here, the RF machine learning technique performs better than the other models for creating LSMs in the Gallicash watershed in Iran.

Model Performance and Comparison
The northern parts of Iran are frequently affected by landslides, which is hindering the economic growth of this region. Landslide susceptibility mapping may contribute to an effective solution for this region to reduce the financial losses and human deaths caused by the landslides. The landslide susceptibility mapping is an invaluable tool for deploying protection measures and managing sudden landslide occurrences in this region. To prepare an LSM, selecting suitable models and geo-environment factors is an important task to properly forecast areas that may eventually be affected by a landslide. In the forested area, it is difficult to investigate the topographic conditions and to map the landslide susceptibility. Due to the obstruction of a large canopy cover, ASTER and SRTM low-resolution images cannot provide sufficient information on topographical characteristics in forest areas [114]. The scarcity of the morphological information in the densely forested areas is creating difficulties in the environmental management in the study area. These problems can be solved using radar, which is an effective and reliable remote sensing technique that provides accurate information regarding natural resources. synthetic aperture radar (SAR) techniques leads to lightweight cost-effective imaging sensors of high resolution [115]. The ALOS PALSAR, ERS 1, 2, JERS Radar sat and Envisat have been used in this field. The ALOS PALSAR DEM is more suitable for collecting information from densely forested areas and also produces a more accurate DEM than ASTER and SRTM [108,114,116,117]. In this study, the PASLAR DEM with a resolution of 12.5 m × 12.5 m was used to extract the topographical and hydrological characteristics in the Gallicash watershed.
Recently, machine learning techniques have become more popular in the scientific community for environmental modeling as these methods are helpful in effectivley explaining the complex relationships between environmental predictors and responses. So far, many machine learning ensemble approaches, including RF, FLDA, BRT, AdaBoost, multiboost, bagging, CART, NBT, ADTree have been used for the evaluation of LSMs [1,[29][30][31]. different approaches and methods were developed and implemented around the world for the spatial assessment of environmental hazards.In this study, the RF, ADTree, and FLDA models are used for the assessment of the LSMs in the Gallicash watershed. The LSMs of the present study were prepared using the R program, remote sensing, and GIS. It is impossible to avoid the problem of overfitting because the sampling ratio of 70:30 was established for the training and test datasets without conducting a sampling ratio accuracy test, and despite a multi-collinearity check having been carried out, noise still exists in the landslide conditioning factors data. A key advantage of the machine learning algorithm used in this research is that it automates the task of searching multiple databases in order to collect crucial information. Such algorithms can compensate for particular theories that can be used together with automating the analysis of large amounts of data to support planning decisions.
The LSMs prepared using the RF, ADTree, and FLDA models were categorized into five separate landslide susceptibility classes, namely very low, low, moderate, high and very high, according to the four classification methods, namely quantile, natural breaks of the Jenk's classification, equal interval, and geometric interval. By comparing the results of each classification method and the distribution of landslides of the training and validation dataset on the high and very high classes of landslide susceptibility, it was noticed that the most reliable distribution was given by the natural break classification method. The main advantage of the natural break classification method is that there is no class bias, and the intra-class deviation is minimum and inter-class deviation is maximum. This is in line with the results of Arabameri et al. [1] since the natural break method is a reliable classifier in the mapping of landslide susceptibility. In the Gallicash watershed, the very high LS classes of the RF, ADTree, and FLDA models cover an area of about 11.36%, 15.37% and 16.10% (Table 4). To determine the goodness of fit and prediction accuracy of LSMs, the AUROC curve, SCAI, efficiency, and TSS was used in this study area. The AUROC curve produced the successive curve rate (SCR) and predication curve rate (PCR) using the training and validation datasets. The succession and the prediction values of the AUROC curve are shown in Figure 8. The AUROC, SCAI, efficiency and TSS methods have shown that all the models are suitable and efficient enough for producing the LSMs. Our results of these models are rational and justified by the works of Arabameri et al. [1], Dieu Tien Bui et al. [31], Lee et al. [40] and Chen et al. [14]. Above all, the study shows that the RF model is more suitable for the landslide susceptibility mapping than the ADTree and FLDA models. The severe landslide risk zones in this region have been found in the middle part along the north-south axis of the study area, likely due to the mountainous, elevated, slopping, high convergence index and topographical position of the landforms in this area.

Variable Contribution Analysis
In this work, a vital machine learning ensemble approach, namely the RF model, was used for landslide susceptibility mapping in the Gallcash river watershed. In the program R, the presence or absence of landslides, expressed as 0 and 1, were considered as the target variables, and the LCFs were considered as the independent variables. The point's values of the studies have been extracted LCFs regarding the target variables. The results of the RF model show that the out-of-bag error is 4.7 % ( Table 6). The RF algorithm calculates the importance of FCFs for the construction of LSM (  (Table 7). The TPI, distance to stream, elevation, and rainfall were shown to have a great effect on landslides, followed by plan curvature, lithology, SPI, CI, TWI, distance from road, NDVI, soil, distance from fault, slope and LULC. The relative importance of each of the LCFs was also calculated using the BRT machine learning technique. The results of the BRT model indicate that the TPI (0.5552) has the highest contribution to landslide susceptibility, followed by distance to stream (0.1148), elevation (0.0573), rainfall (0.0497), plan curvature (0.0419), lithology (0.0306), SPI (0.0270), CI (0.0225), TWI (0.0224), distance to road (0.0219), NDVI (0.0218), soil type (0.0165), distance to fault (0.0086), slope (0.0075), LULC (0.0017) and LS (0.0003). All the LCFs used in susceptibility modeling play a major role in making the area prone to landslides (Figure 9).

Conclusions
The presented study analyzes landslides and how they are impacted by tectonic activities, adverse topographical and hydrological conditions, land-use change, road construction and unnecessary activities of people in the Gallicash catchment, northern Iran. A total of 249 landslides were identified and analyzed to evaluate the LS conditions in the study area. LSMs based on an environmental hazards assessment is an important aspect of sustainable environmental management and planning. In this study, sixteen landslide conditioning factors and three reliable machine learning ensemble models, namely RF, ADTree and FLDA, were chosen to demonstrate the areas at risk of landslides. Landslide inventory datasets were considered to be the dependent variables and LCFs were presumed to be independent variables. The dependent variables are expressed as the binary values of 0 and 1, where 1 and 0 indicate the presence and absence of landslides. The high-quality data, such as the ALOS PALSAR DEM and Sentinel-2 image with a resolution of 12.5 m × 12.5 m were used for extracting the topographic, hydrological and geomorphological datasets to prepare LSMs. Of the different lithological segments in this study area, the H (Tre and TRJs lithological units) and B Figure 9. The relative importance of landslide conditioning factors using boosted regression tree model.

Conclusions
The presented study analyzes landslides and how they are impacted by tectonic activities, adverse topographical and hydrological conditions, land-use change, road construction and unnecessary activities of people in the Gallicash catchment, northern Iran. A total of 249 landslides were identified and analyzed to evaluate the LS conditions in the study area. LSMs based on an environmental hazards assessment is an important aspect of sustainable environmental management and planning. In this study, sixteen landslide conditioning factors and three reliable machine learning ensemble models, namely RF, ADTree and FLDA, were chosen to demonstrate the areas at risk of landslides. Landslide inventory datasets were considered to be the dependent variables and LCFs were presumed to be independent variables. The dependent variables are expressed as the binary values of 0 and 1, where 1 and 0 indicate the presence and absence of landslides. The high-quality data, such as the ALOS PALSAR DEM and Sentinel-2 image with a resolution of 12.5 m × 12.5 m were used for extracting the topographic, hydrological and geomorphological datasets to prepare LSMs. Of the different lithological segments in this study area, the H (Tre and TRJs lithological units) and B (Dp, and DCkh lithological units) geological sections have a strong influence on the occurrence of landslides. Similarly, alfisoil and inceptisoil soils contribute significantly to landslides occurring in this region. The machine learning software R Studio was used to run machine learning algorithms such as RF, ADTree, FLDA, etc. RS and GIS are the basic tools for performing the landslide susceptibility mapping. The integrations of R, RS and GIS have created a strong foundation for the construction of the LSMs using machine learning algorithms. The results of RF, ADTree and FLDA models show that 11.36% (RF), 15.37% (ADTree), and 16.10% (FLDA) of the area of the whole catchment belong to the very high landslide susceptibility class. Several common techniques, such as the ROC curve, SCAI, Efficiency and TSS, were used to verify the performance of the three chosen models. All models have strong LS simulation abilities. Finally, the RF model shows slightly better accuracy (AUC of SRC = 0.978, AUC of PRC = 0.985, the efficiency of VD = 0.77, efficiency of TD = 0.77, TSS of VD = 0.55 & TSS of TD = 0.53) than the other models in modeling the landslide susceptibility in this region as per the results of the validation techniques. The most critical zones of landslide susceptibility have been identified in the central part along the north-south axis in the Gallicash catchment. This could identify the relationships between the spatial distribution of landslide locations and landslide conditioning factors. The LSMs prepared based on the RF, ADTree, and FLDA may be utilized by decision-makers, planners, academic researchers, engineers, and government agencies to mitigate the loss of human life and property. The results of this study will also support sustainable development measures and management of the environment in the Gallicash catchment and landslides prone areas in general.