Artificial intelligence approaches for spatial prediction of landslides in mountainous regions of western India

The prediction of landslide is a complex task but preparing the landslide susceptibility map through artificial intelligence approaches can reduce life loss and damages resulting from landslides. The purpose of this study is to evaluate and compare the landslide susceptibility mapping (LSM) using six machine learning models, including random forest (RF), deep boost (DB), stochastic gradient boosting (SGB), rotation forest (RoF), boosted regression tree (BRT), and logit boost (LB) in the mountainous regions of western India. The landslide inventory map consisting of 184 landslide locations has been divided into two groups for training (70% data set) and validation (30% data set) purposes. Fourteen landslide triggering factors including slope, topographical roughness index, road density, topographical wetness index, elevation, slope length, drainage density, stream power index, geomorphology, rainfall, soil, lithology, lineament density, and normalized difference vegetation index have been considered using the boruta approach for the LSM. The results reveal that the RF model has the highest precision in terms of area under curve (0.88; 0.89), kappa (0.62; 0.50), accuracy (0.81; 0.77), and specificity (0.86; 0.86) both in the study region and secondary region, respectively. Hence, it can be concluded that the RF is an effective and promising technique as compared to DB, SGB, RoF, BRT, and LB for landslide susceptibility assessment in the research area as well as in regions having similar geo-environmental configuration.


Introduction
Globally, the frequency of natural disasters has been dramatically increased over the last few decades (UNISDR 2015). Each year, a loss of billions of dollars and thousands of lives due to natural disasters, such as earthquakes, tsunamis, cyclones, floods, and landslides are recorded (Yesilnacar and Topal 2005;Yilmaz 2009). Almost there is no region in this earth which is free from the impact of such devastating natural hazards (Suzen and Doyuran 2003). The occurrence of landslides constitutes about 17% of the natural disasters worldwide as reported by the Center for Research on the Epidemiology of Disasters (Pourghasemi et al. 2012;Nohani et al. 2019). Yilmaz (2009) and Pourghasemi et al. (2012) believe that landslide events will tend to rise in the next decades owing to deforestation, land-use change, and changing climatic condition. Many scientists have assessed landslide hazards and delineated the high-risk zones in past few decades. In general, landslide mapping techniques can be divided into two major methods, i.e., direct and indirect. In case of direct mapping, geomorphologists determine the degree of susceptibility based on the experience and knowledge of terrain characteristics, hillslopes etc. (Van et al. 2003). However, indirect mapping utilizes different statistical or high-performance deterministic machine learning or hybrid models to predict the landslide probable areas on the basis of previous landslide distribution. In the preceding two decades, multitudes of LSM studies have been undertaken with increasing popularity of Geographic Information System (Aleotti and Chowdury 1999;Guzzetti et al. 2006).
The pioneering works of the landslide mapping were done in 1970s (Carrara et al. 1977;Brabb et al. 1978). Since then, various approaches and methods have been put 720 Page 2 of 20 forwarded and employed for LSM, including heuristic terrain and susceptibility zoning, geomorphological mapping, statistical modelling, physically based numerical modelling etc. (Montgomery and Dietrich 1994;Hansen et al. 1995;Aleotti and Chowdhury 1999;Guzzeti et al. 1999;Lee et al. 2001; Van et al. 2008;Reichenbach et al. 2018). However, the above approaches are no longer capable, because they were restricted to large amounts of field data in a smaller area. From the beginning of the current century, statistical modelling with high-resolution GIS data has been in use in different geographical regions for LSM, including frequency ratio (Yilmaz 2009;Nohani et al. 2019), logistic regression (Pourghasemi et al. 2013;Zhao and Chen 2020), weight of evidence (Lee and Choi 2004;Regmi et al. 2010). The statistical methods are less accurate compared to the machine learning methods, because it can overcome the statistical limitations and predict non-linear structures (Tehrany et al. 2013;Costache and Bui 2019;Prasad et al. 2020). Hence, in recent years, several machine learning models have been developed and applied successfully for LSM, including random forest (Kim et al. 2018;Arabameri et al. 2019), rotation forest (Hong et al. 2018;Pham et al. 2020a), support vector machine (Chen et al. 2017;Phong et al. 2019), boosted regression tree (Youssef et al. 2016;Park and Kim 2019), naive bayes (Bui et al. 2012;Pham et al. 2020a), adaboost (Hong et al. 2018;Wu et al. 2020), and artificial neural network Yilmaz 2009). Such methods have also been employed in various spatial prediction mapping, viz. flood susceptibility Pham et al. 2020b;Shahabi et al. 2020), land subsidence Oh et al. 2019), gully erosion (Amiri et al. 2019;Pourghasemi et al. 2020;Chen et al. 2021b), piping erosion (Hosseinalizadeh et al. 2019;Chen et al. 2021aChen et al. , 2021bChen et al. , 2021c, land-use change (Aburas et al. 2019;Jamali 2020), and groundwater potentiality (Rahmati et al. 2016;Naghibi et al. 2019). Besides statistical and machine learning models, novel models with high accuracy have been also used in landslide mapping Nguyen et al. 2017;Abedini et al. 2019;Chen and Li 2020;Pham et al. 2020a).
In the aforementioned studies, the performance of the applied models has been assessed with respect to particular study regions. It is imperative to show the model applicability in another region for general findings (Prasad et al. 2020). To the best of our knowledge, in landslide susceptibility mapping, the models have been validated by applying them in more than one geographic region. Thus, the present study aims to map the landslide-prone zones in the western India through different approaches of artificial intelligence and to validate the models in another region. In spite of being a high landslide prone area, the present area of investigation has not received much attention from the researchers. The current study is the first attempt to prepare the landslide susceptibility map hereof using machine learning approaches. The LSM is foremost mitigative step to reduce the havoc fatalities and infrastructural damages as a dire consequence of frequent severe landslides along the precarious Ghats escarpment triggered by torrential rainfall during the Indian Summer Monsoon. In addition to the comparison of different machine learning approaches, the scientific delineation of landslide-prone regions will be helpful for proper planning and management of the study area.

Study area
The study area lying approximately between 17°50 to 18°19 North latitudes and 73°18 to 73°40 East longitudes, comprises an approximate area of 1180 km 2 (Fig. 1). The study region represents a hilly topography ranging from 4 to 1404 m above the mean sea level. The average slope of the area is about 15°, though some areas are close to the vertical slope. The study area shows a general tilting facing towards the Arabian Sea.
Geologically, the area spans from upper Cretaceous to lower Eocene, which comes under the Sahyadri group of formation. From the geomorphological perspective, more than fifty percent of the area belongs to structural origin highly and moderately dissected plateau. The climatic condition of the study region is subtropical monsoon type. The minimum and maximum average temperatures of the area are 25 °C and 30 °C, respectively. The Arabian branch of south-west monsoon brings huge amount of cloud that leads to intense orographic rainfall along the hilly slopes. The area receives about 3000 mm annual mean rainfall during the monsoon seasons (June-September) (Department of Agriculture, Maharashtra State). The principal river Savitri, descends from the Mahabaleshwar peak of Western Ghats and flows for about 110 km before it empties into the Arabian sea.

Landslide inventory map
Inventory data is crucial for any natural hazard mapping especially to understand the landslide distribution and effective variables (Kadavi et al. 2018;Panahi et al. 2020). In the present study, the landslide locations were collected from the Public Works Department (PWD) office, Google Earth historical images and intensive field survey (Fig. 2). A total of 184 landslide locations were taken into account for this study. The data set was divided into two groups 70% (129 locations) as training data for model fitting and the rest 30% (55 locations) as testing data for model predictive capability assessment (Fig. 1). Simultaneously, an equal number of non-landslide localities were recognized from the topographical map, field survey, and satellite images for the purpose of model establishment. Likewise, the inventory map of another region was generated (Fig. 1).

Landslide triggering factors
The important step for the LSM is to prepare a spatial database of landslide conditioning factors (Reichaenbach et al. 2018;Nohani et al. 2019). The predisposing variables of the landslide were selected after extensive literature survey and field investigation (Hong et al. 2018;Wang et al. 2019). These factors are: slope, elevation, slope length (LS), topographical wetness index (TRI), stream power index (SPI), topographical wetness index (TWI), drainage density (DD), road density (RD), lineament density (LD), normalized difference vegetation index (NDVI), rainfall, lithology, geomorphology, and soil ( Fig. 3a-n). The terrain variables (elevation, slope, LS, SPI, TWI, and TRI) were extracted from the SRTM digital elevation model (DEM) with 30 m spatial resolution. Each conditioning variable was transformed into a spatial database of 30 × 30 m size and the grid of the study region was generated by 1337 columns and 1740 rows (1,312,960 pixels; 1180 km 2 ).
The slope is one of the most important triggering variables, since it directly influences the landslide occurrences (Yalchin 2008). The slope map of the study region was categorized into five classes: 0-7.19, 7.19-15.25, 15.25-24.18, 24.18-35.98, and 35.98-73.40 degrees. The altitude of the area is classified into five classes. LS represents slope steepness and slope length of the topography, which affects the soil loss potential (Naghibi et al. 2019). The LS values are extracted by the given mathematical formula (Moore and Burch 1986): where A and θ indicate the flow accumulation and slope values in degree, respectively.
The ratio between the land surface and flat surface in a specific area is defined as TRI (Stambaugh and Guyette 2008). TRI map is derived from the following formula: where maximum and minimum altitude are written as max and min, respectively.
(1) LS = (A × cellsize∕22.13) 0.4 × (sin ∕0.0896) 1.3 (2) TRI = Abs(max 2 − min 2 ) SPI and TWI are vital hydrological factors for LSM. The erosive power and discharge in a specific area are expressed by SPI Chapi et al. 2017). TWI is defined as the ratio between the catchment area and slope (Gallant 2000;Chapi et al. 2017). Both variables were calculated from the catchment area (A) and slope (β) by the given formulae : The stability of slope is controlled by the saturation level of material, which is indirectly related to the drainage density (Yalchin 2008). The construction of roads beside high slope area reduces burden on both the heel of slope and topography that causes landslide (Yalchin 2008). Lineament is the linear or curvilinear feature on the earth surface (Das et al. 2018). It represents weak zones, where the slope instability leads to the landslide. Lineament map of the study area was prepared from the shaded relief maps. The density maps of drainage, road, and lineament were prepared in Arc GIS 10 environment using spatial analyst tools.
To prepare the NDVI thematic layer, band 4 (red band) and band 5 (infrared) of Landsat 8 OLI (https:// earth explo rer. usgs. gov/) image (date: 29.10.2018) were used. The NDVI value was extracted as follows: where NIR and R denote the infrared and red bands of the electromagnetic spectra, respectively. The lithology map obtained from the Geological Society of India (GSI) was reclassified into four broad categories: Karla, Indrayani, Diveghat, and Purandgarh formation on 1:250,000 scale. The geomorphology map of the study area was reclassified into five major groups, namely, structural origin moderately dissected plateau (SoMDp), structural origin highly dissected plateau (SoHDp), structural origin low dissected plateau (SoLDp), denudational origin moderately dissected lower plateau (DoMDLp), denudational origin pediment pediplain complex (DoPPc) on 1:50,000 scale based on Landsat 8 image, NRSC base map, topographical maps and field survey. The thematic map of soil was modified using the base map of the National Bureau of Soil Survey and Land Use planning (NBSSLUP) on 1: 500,000 scale. The soil map was categorized into three broad classes: Typic Ustorthents, Lithic Ustorthents, and Udic Rhodustalfs. Rainfall has been the most influential factor in landslide occurrences especially in the Indian sub-continent due to arrival of the southwest monsoon. The rainfall data of the last 5 years (2013-2018) was collected from the rain gauge stations of the department of agriculture, Maharashtra state. Based on this data, rainfall map was classified into four classes using the inverse distance weighted method. Among the classification methods, natural break method was extensively applied for classification of landslide variables (Phong et al. 2019;Chen and Chen 2021). In the present study, natural break method was used for the classification of the causative factors of landslide in ArcGIS platform.

Methodology
The methodology flow-chart of the present study is presented in Fig. 4.

Boruta and multicollinearity feature selection methods
Variable selection is the essential function in natural hazard mapping. The factor determination method improves the model performance by increasing the data quality and reducing noise, overfitting, dimensionality of feature space and preventing redundancy problems ). There are two popular features choosing methods for modeling including wrapper and filter accessible. In the present study, the boruta algorithm was used for feature selection based on wrapper method. The classification in boruta approach is implemented by voting of multiple unbiased weak classifier decision trees. The scale of importance of the variable is obtained as the loss of classification accuracy due to the random alternate of factor values between aims. The algorithm works as follows: (i) At the first step, it adds randomness to the dataframe by shuffling each feature which are called shadow feature (ii) Then, it trains the RF classifier to measure the Z score values and checks whether the attribute has a higher Z score than the maximum Z score among shadow attributes (MZSA). (iii) Finally, the features having value less than MZSA are rejected for being irrelevant and vice versa (Kursa and Rudnicki 2010;Amiri et al. 2019;Prasad et al. 2021). Shadow minimum, mean, and maximum values represent the Z score of a shadow attribute. On the other hand, tentative variables lying close to the maximum shadow (maximum Z score), are confirmed or rejected based on the mean shadow (mean Z score).
In addition, multicollinearity analysis was performed to verify the results of boruta method. Multicollinearity is the popular statistical method to find out whether any positive correlation between the conditioning factors exists or not. Researchers like Chen et al. (2018a), Pourghasemi et al. (2020) showed that the influencing variables containing the variance inflation factors (VIF) value more than 5 and tolerance (TOL) value less than 0.2 were subjected to multicollinearity, while in the studies of Nguyen et al. (2017), Wang et al. (2019), the VIF value > 10 or TOL value < 0.1 of the controlling factors were rejected.

Landslide events and variable inter-relationship by Weight of evidence (WoE) method
The WoE is a quantitative approach used to forecast the occurrence of hazard events (Armas, 2012). Based on Bayes rule, the approach was primarily propounded for mineral potential mapping (Bonham-carter 1994;Chen et al. 2018a).
With the march of time, WoE approach was popularized in other spatial prediction mapping (Rahmati et al. 2015;Aghdam et al. 2017;Chen et al. 2018a;Arabameri et al. 2018). In this study, WoE method was used to determine the correlation between the landslide events and effective parameters. The primary purpose of this model is to measure the positive (W + i ) and negative weights (W − i ): where P and ln stand for the probability and natural logarithm, respectively. B and B indicate the presence and absence of desired sub-class of landslide predisposing factors, respectively, whereas A and Ā are the presence and absence of landslide events, respectively. (2001), is a nonparametric method for classification and regression analysis, The RF algorithm creates many decision trees to enhance the performance of the model (Breiman 2001; Rahmati et al. 2016). Each decision tree is originated by bootstrap samples and leaves around one-third of samples for validation purpose using the out of bag (OOB) error. The advantages of the model are: i) large data set handling, ii) no overfitting, iii) OOB error estimation, iv) low bias and low variance, v) no requirement of prior data transformation and rescaling, and vi) higher prediction performance. To tune the model, the algorithm requires two important criteria: number of trees (n tree ) and number of features (m try ). Deep boost is a new ensemble learning algorithm formulated by Cortes, Mohri, and Syed (2014). The method is updated and some new solutions that are left about the theory underlying AdaBoost. It consists of deep decision trees or members of complex families. The method is a capacity-conscious criterion to achieve high accuracy without overfitting the data. The advantage of the algorithm is to minimize the corresponding learning bound and the mixture weight assigned to each sub-family (Cortes et al. 2014;Chen et al. 2021a, b, c). The algorithm can extend by examining non-differentiable convex surrogate losses. As far as we are concerned, the DB model has been used for the first time in natural hazard mapping.

Random forest introduced by Breiman
Stochastic gradient boosting is an ensemble machine learning technique that merges both the advantages of bagging and boosting methods and is used for classification and regression purpose (Friedman 2001(Friedman , 2002Zhou et al. 2016). SGB is the modification of adaboost and adabtive bagging method. SGB has many advantages, such as i) it uses random sub-sample (selected without replacement) of the data set at each stage of boosting process instead of the whole data set; ii) at each step, small trees are generated rather than large entire classification trees; iii) since small trees are developed, SGB method is highly resistant to overfitting; and iv) the steepest gradient algorithm of the model is resistant to outliers (Friedman 2002;Lawrence et al. 2004;Moisen et al. 2006).
Rotation forest is an ensemble learning method devised by Rodriguez et al. (2006) to enhance the predictive performance of individual classifiers. In this method, every single classifier is merged with one another and each classifier may be trained using different i) training data sets; ii) subsets of features; iii) parameters of the classifier; or iv) classifier models (Polikar et al. 2010;Nguyen et al. 2017). The model is a combination of the random subspace and bagging method with principal component analysis (PCA) for the ensemble classifier. By applying PCA technique, the training data set was randomly classified into K subsets (Rodriguez et al. 2006;Nguyen et al. 2017;Naghibi et al. 2019).
Boosted regression tree utilizes both machine learning and statistical techniques for classification and regression problems (Youssef et al. 2016;Naghibi et al. 2016). BRT algorithm uses boosting technique to combine a large number of simple tree for maximizing the predictive capacity (Elith et al. 2008;Naghibi et al. 2019). The advantages of this model are: i) it replaces the missing data using surrogates, ii) it can deal with any type of predictive variable (numeric, categorical, binary), and iii) the predictive value get affected by monotone transformation, and different scales of measurement. To run the BRT model, three parameters are required to optimize it, namely, interaction depth, shrinkage or learning rate, and number of trees. Interaction depth represents the number of nodes and the shrinkage indicates the contribution of trees.
Logit Boost developed by Friedman and Tibshirani (2000), uses boosting method to reduce the bias and variance (Oh et al. 2019). The LB method is a minor variation of popular boosting technique (Adaboost). In this method, the additive logistic regression function is used to classify the data set (Tehrany et al. 2018). The main advantages of the LB technique are: i) it can handle noisy data; ii) it has the ability to reduce the training errors; and iii) it is used for the multiclass tasks (Oh et al. 2019). The vector values of the model define the input variables and two output groups (landslide and non-landslide events).

Validation of the models
Validation is an imperative stage in machine learning modelling for the evaluation of the study, without which the models have not any scientific significance (Rahmati et al. 2016;Naghibi et al. 2016;Prasad et al., 2020). In the validation data set, dependent variables are not considered for verifying the model efficiency to predict the location of landslides. In the present study, the models were validated based on a set of quantitative parameters, such as accuracy, kappa index, sensitivity, specificity, and area under the receiver operating characteristics curve (AUROC). These quantitative criteria were assessed using the possible consequences (PC), namely, true positive TP), false positive (FP), true negative (TN), and false negative (FN). The detailed information of these quantitative parameters is available in Chapi et al. (2017).
The ROC curve is a scientific method to examine the overall performance of the models (Rahmati et al. 2016;Chapi et al. 2017;Bui et al. 2019). The AUROC is a quantitative index for model capability to predict correctly the occurrence or non-occurrence of specific events. Area under curve (AUC) value lying in between 0 to 1, was classified as poor (0.5-0.6), average (0.6-0.7), good (0.7-0.8), very good (0.8-0.9), and excellent (0.9-1) (Chen et al. 2018b). The higher AUC value, more will be the model perfection (Chapi et al. 2017;Naghibi et al. 2017;Das 2019Das , 2020. AUROC value is calculated as follows: where the P and N refer to the total numbers of landslide occurrences and non-landslide occurrences in the study area, respectively.

Feature selection and important landslide triggering variables
The selection of significant conditioning factors in machine learning techniques is of vital importance to improve model accuracy by eliminating impertinent factors Bui et al. 2019). In the current research, the boruta technique was used for the relevant variable selection. From the nineteen landslide predisposing factors, boruta evaluation elected fourteen suitable variables (slope, TRI, RD, TWI, elevation, LS, DD, SPI, geomorphology, rainfall, soil, lithology, LD, and NDVI), whereas five factors including standard curvature, plan curvature, profile curvature, land-use, and aspect were removed (Fig. 5a). These rejected variables were also excluded in the study of Chen et al. 2018b. The boruta approach is capable of determining the mean, median, minimum, and maximum importance of the factors displayed in Table 1. In congruence with the result of boruta algorithm, the result of multicollinearity analysis presented in Table 2, showed that the standard, plan and profile curvature were disqualified, since there VIF value exceeded 10 and TOL fell below 0.1.
After choosing the significant factors that are crucial for the modelling, the variable importance function of RF method was used to measure the feature importance for LSM in the study area. The measured importance values of the landslide conditioning variables were 100 ( (soil), and 0 (lithology). These results of the variable importance analysis revealed that the slope, TRI, elevation, and TWI had more influence on landslides, followed by LS, DD, RD, rainfall, SPI, NDVI, LD, geomorphology, soil, and lithology (Fig. 5b). From the literature, it was noticed that the contribution of the conditioning factors in a particular region varies with the use of different measuring methods (Chen et al. 2018b). In the case of present study, variable importance results of both boruta and RF methods were almost same that confirmed the predictive ability of the algorithms.

Spatial interaction of landslides and conditioning factors using WoE model
The spatial interconnection between landslide events and affecting variables applying WoE method is presented in Table 3. In this table, C/SC value is the final weight used for the above relationship. The C/SC value of the present study was lying between -9.60 and 9.11. The higher positive value shows strong correlation between landslide events and conditioning factors, whereas negative weight indicates poor relation. Regarding slope, increasing slope angle is directly related to the landslide events in a certain range Chen and Chen 2021). This fact was reflected from the current research, where most of the landslide happened in the higher slope angle classes viz. 24°-34° and 35°-73° with 8.87 and 7.69 C/SC values, respectively, whereas lower slope angle held negative weights. The elevation is another significant factor of landslide events as shown in the study of Wang et al. (2019) and Arabameri et al. (2019). The analysis of WoE for landslide events and elevation expressed that the moderate elevation (430-713 m) recorded the maximum C/ SC value (9.11). On the other side, lower (4-154 m) and higher (920-1404 m) elevation classes had negative weights. The increasing slope length is favorable for water infiltration into sub-soil that enhances the possibility of slope failure ). The LS value was highest (7.61) in the class of 0.42-2.10, indicating strong relation with landslide occurrences. The weight of the C/SC went upward with increasing topographic roughness. The C/SC values   . The longer period of monsoon rainfall in the study area had direct impact on landslide occurrence because of more runoff caused soil loosening and sediment flow. The maximal rainfall class (619-656 mm) in the area under study indicates that around fifty percent of landslides took place in this group. Regarding lithology, the highest weights were obtained from the bdd (3.32) and bdp (2.47) groups. Since, the rock of these groups belongs to 350-450 m elevation (landslide occurrence zone), it can be inferred that they are more weathered. For geomorphology, structural origin highly dissected plateau has strongly correlated with the landslide events of C/SC value 6.98. This sub class also obtained maximum score by AHP method in the research of Bera et al. 2019.
In the case of soil, Udic Rhodustalfs soil was only positive weight (6.21) among the soil groups.

Model construction
In machine learning models, researchers applied 10 crossvalidation method to enhance the precision of classification and to reduce the noise . Therefore, the present work adopted this method with five repetitions using 70% of the landslide and non-landslide points. These points along with their extracted fourteen conditioning factor values comprised the training data set, whereas the remaining points were used for validation purpose. The exact pixel values of the derived parameters influenced the validation result due to the higher range of the values. To solve this problem, values were transformed into 0 to 1 range with the help of min-max normalization. After the model establishment with training data set, testing data set was utilized to verify the performance and predictive capability of each model. This complete procedure was done using the R and Arc GIS software.

Performance and comparison of the models
In this study, AUC and selected statistical measurements were used for evaluation of the models. These criteria considered both training and testing results for model comparison and application. The training data set defines the model fitness, whereas the testing data set represents the predictive capability of the model Prasad et al. 2021). In the case of training data set, among the six models, the RF performed excellently with the highest value of 1 quantified by accuracy, Kappa, sensitivity, specificity, and AUC in both the primary and secondary study area (Table 4). Landis and Koch (1977) determined that the kappa value > 0.8 means the perfect match of measured value with observed value. The value of 1 for sensitivity and specificity represent that all landslide and non-landslide points were correctly classified. The AUC values of the remaining models for primary and secondary regions were 0.89, 0.99 (DB); 0.87, 1 (SBM); 0.95, 0.97 (RoF); 0.82, 0.89 (BRT); and 0.86, 1 (LB), respectively. The highest AUC value of RF method in both the region manifested the better goodness of fit compared to the rest. The established models were corroborated with the testing data set to check the predictive ability of the models. The results of validation data set revealed that RF method held the highest as measured by accuracy (0.81), kappa (0.62), specificity (0.86), and AUC (0.88) compared to other models (Fig. 6a, Table 5). It is also observed that the RF model appeared to be ideal in the secondary region based on accuracy (0.77), kappa (0.50), specificity (0.86), and AUC (0.89) (Fig. 6b, Table 5

Landslide susceptibility map construction
The preparation of the landslide susceptibility maps is essential to know the possibility of future landslide occurrences in an area. After the training processes, each pixel of the study area was computed for the prediction of landslide. These values ranging from 0 to 1 indicate low to high chances of landslide occurrences. In general, the classification of the derived probability values includes automatic and user-defined methods (Chen et al. 2021a). However, there are several standard techniques in GIS to reclassify the probability maps, namely, natural breaks, min-max normalizations, standard deviations, equal intervals, geometrical, and quantiles. A specific method is selected based on the research objective and distribution of data (Chapi et al. 2017). It was learned from the previous studies that many researchers had employed the natural break method for generating the susceptibility maps (Hong et al. 2018;Park and Kim 2019;Panahi et al. 2020;Chen et al. 2021a). This classification method has advantage of no class bias as well as minimum intra-class and maximum inter-class deviations (Arabameri et al. 2020). In the present work, the distribution of probability values of the models abruptly changed. For this reason, the landslide susceptibility maps were classified using the natural break method in ArcGIS environment into four groups, including low, moderate, high, and very high (Fig. 7). The low and moderate prone zones of the entire study area account for 22.74%, 40.06 (RF); 30.08%, 40.94% (DB); 38.99%, 30.37% (SBM); 54.59%, 28.75% (RoF); 57.13%, 15.93% (BRT); and 9.02%, 28.8 (LB), respectively. The high susceptible classes of RF, DB, SBM, RoF, BRT, and LB models constitute 26.06%, 20.13%, 19.37%, 13.02%, 13.51%, and 44.3% of the study region, respectively. The share of 10.58% (RF), 8.83% (DB), 11.25% (SBM), 3.61% (RoF), 13.42% (BRT), and 17.85% (LB) belong to very high vulnerable zones (Fig. 8). The lowest minimal variations of high and very high susceptibility classes indicate the high model precision (Naghibi et al. 2017;Prasad et al. 2020). In the present study, minuscule difference between RF, DB, and SBM, models with regard to high and very high landslide prone zones ensured the models high accuracy (Fig. 8). The low susceptibility classes of the RoF and BRT models occupied much larger area than the moderate susceptibility classes in compared to other models (except LB model). This is mainly on account of the fact that the degree of prediction varies model to model for their different algorithm. Regarding LB model, susceptible classes did not correspond with the classes of other models because of less precision in prediction. The study pointed out that the north-eastern, eastern, and south-eastern parts of the research area are more vulnerable to landslides owing to the steepness, roughness, and wetness of these mountainous regions.

Variable contribution analysis
Although many studies have been accomplished worldwide on LSM, still there is no specific set of variables (Van et al. 2006;Chen and Li 2020). The influencing factors vary from one study region to another due to heterogeneity of the earth. In the current study, a total of nineteen factors were chosen consulting recent studies carried out in the different regions with more or less same elevation as of the study area under consideration and subsequently fourteen relevant factors were selected based on boruta algorithm for LSM (Hong et al. 2018;Kadavi et al. 2018;Nohani et al. 2019;Arabameri et al. 2019;Wang et al. 2019). It is worth mentioning that the result of boruta approach almost matched with  TP  40  42  39  39  34  35  10  12  11  11  11  9  TN  50  47  44  44  49  46  24  21  22  21  21  22  FP  8  11  14  14  9  12  4  7  6  7  7  6  FN  13  11  14  14  19  18  6  4  5  5  the result of multicollinearity which reflects the robustness of boruta algorithm in feature selection. It can be seen that the topographical factors, such as slope, TRI, and elevation exerted more influence on landslide occurrences compared to the categorical factors (lithology, soil, and geomorphology) in the study area. Researcher inferred that topographical variables are the most responsible for landslide incidents (Nakileza and Nedala 2020; Arabameri et al. 2020). This is because of the topographical factors mainly slope angle which plays a key role in slope stability and controls the shear and normal strength on the plane of weakness. The study area is basically more or less homogeneous bedrock region (Sahyadri group) with a thin veneer of top soil. Therefore, the factors lithology and soil did not have much influence on landslide occurrence. Most of the slides in this area that occurred within accessible locations are rock falls or debris flow. On the other side, WoE method expressed the influence variability of sub classes of each landslide governing factor. Based on the result, the highest C/SC values noted from the slope and elevation sub classes revealed that the landslide occurrences mostly took place into moderate elevation (430-713 m) with higher slope angle (24-73°) probably due to more weathering at this height of the study area.

Previous works and future work proposal
Mathematically, landslide susceptibility is the likelihood of the occurrence of slope failures in location, while the geo-environmental condition is fixed (Guzzeti et al. 2005). LSM is generally done to predict the potential geographic locations of future landslides (Guzzeti et al. 1999;. In recent years, machine learning approaches have become more popular and outstripped other conventional methods in landslide susceptibility and hazard mapping Arabameri et al. 2020). For the mapping of potential landslide locations, a large number of machine learning methods have been developed, artificial neural network (Lee et al. 2001), support vector machine (Dou et al. 2015), random forest (Miner et al. 2010), neuro-fuzzy (Pradhan 2013), naïve bayes (Bui et al. 2012), quadratic discriminant (Rossi et al. 2010), adaptive neuro-fuzzy inference systemsatin bowerbird optimizer (Chen et al. 2021a) just to name a few. However, there is no universal model for spatial mapping of landslide but the applied model can be appropriate for similar geo-environmental settings Prasad et al. 2020). Therefore, the present study aims to select suitable method for LSM in the mountainous regions of western India. In this research, six machine learning approaches, including RF, DB, SBM, RoF, BRT, and LB were compared in the study area as well as in the secondary region to choose the significant model. Among these models, RF model showed the highest precision, consistent performance, and minimal difference between training and validation results, followed by DB, SBM, RoF, BRT, and LB for LSM in both the regions. The robustness of the RF method is also manifested in other spatial prediction mappings with very high precision (Catani et al. 2013;Naghibi et al. 2019;Costache and Bui 2019;Amiri et al. 2019;Prasad et al. 2020). A study by Chen et al. (2018b) in Chongren County, China, located at about an elevation (1-1230 m) as of the current study area, RF model proved to be the most suitable with respect to the performance of other applied models in landslide prediction. Arabameri et al. (2018) found that the RF model was the best fit in the Gallicash River Watershed, Iran. Bera et al. (2019) worked on LSM in mountainous region of western India. In this study, accuracy of fuzzy gamma operator (81%) and analytical hierarchy process (AHP) (69%) models were lesser than the present research. In the above research works, the researchers employed several methods including novel model only to one respective study area, which restricted the applicability of the models. To widen applicability of the models, it is important to apply the selected models in alike geographic environment (Prasad et al. 2020). In the present study, the outcomes of the models were almost same for both the primary and secondary regions. Therefore, the comparability of the models remained consistent with the either regions under investigation. From the results, it can be stated that the RF model is appropriate for LSM in the western India and also in regions with similar geo-environmental conditions. Likely, in future, more researches are required to conduct in different geographical areas around the world for the assessment of applicability of machine learning models. Furthermore, it can be suggested that the integration of SAR and optical data sets may improve the precision of the models and help to identify the landslide points in a larger area.

Advantages and disadvantages of the models
The machine learning method is applied in natural hazard and resource mapping because of its unique advantageous capability of automatically performing the task of searching several data sets aimed at obtaining necessary information . It is also noted that the predictive accuracy of the model, to some extent, depends on the quality of data and the establishment of the model from training data set (Nohani et al. 2019). In present study, an attempt was made to find out the probable reasons behind the better performance of RF method and comparatively poor achievement of LB model. RF method suitably fitted resulting in reduction of the error for validation data set in the primary and secondary regions. The reason for better performance of the RF method may be due to composition of multiple decision trees and interconnection capability of conditioning variables and non-linearity (Catani et al. 2013;Prasad et al. 2020). Besides, RF method escapes the overfitting of data set which is common problem in artificial intelligence approaches. In the case of LB method, the performance is relatively less in either region. The underperformance of LB model may be because of i) steady property indicated by logistic loss; ii) any change of output value does not affect the loss value; iii) oversimplification in model establishment can reduce the validation accuracy; and iv) problem in node spit gain and value fitting (Sun et al. 2014).

Conclusions
The Landslide susceptibility mapping through high performing artificial intelligence models is an indispensable step for downscaling future loss of properties and human lives. Therefore, in this study, six machine learning approaches (RF, DB, SGB, RoF, BRT, and LB) were employed for effective LSM. To confirm the model applicability, the selected methods were tested in another region as well. The major findings of the study are as follows: i) Among the nineteen landslide causing factors, fourteen (slope, elevation, LS, SPI, TRI, TWI, RD, LD, DD, NDVI, rainfall, lithology, soil, and geomorphology) were selected and remaining five (aspect, land-use, plan curvature, profile curvature and standard curvature) were discarded with the help of boruta approach for LSM.
ii) From the attribute evaluation of WoE method, it was found that the landslide events mostly occurred in the moderately elevated area within higher slope angle, where the topographical roughness is maximum. As per the variable importance of RF method, topographical factors were most important compared to the other factors for the LSM.
iii) Out of six methods, RF model is the most desirable for its superlative performance in model fitness and validation in both the primary and secondary regions. Besides, novel application of the second most suitable model DB in LSM has been successful which widens its scope in natural hazard mapping. iv) According to landslide susceptibility classes of RF method, about 11% and 26% of the study area come under the very high and high landslide vulnerable zones, respectively.
The outcomes of this study can be useful to the planners for planning and management of future landslides events. To meet the purpose, RF method can be opted for preparing LSM in alike geographic environment.