Landslide susceptibility zonation using statistical and machine learning approaches in Northern Lecco, Italy

This study deals with landslide susceptibility mapping in the northern part of Lecco Province, Lombardy Region, Italy. In so doing, a valid landslide inventory map and thirteen predisposing factors (including elevation, slope aspect, slope degree, plan curvature, profile curvature, distance to waterway, distance to road, distance to fault, soil type, land use, lithology, stream power index, and topographic wetness index) form the spatial database within geographic information system. The used predictive models comprise a bivariate statistical approach called frequency ratio (FR) and two machine learning tools, namely multilayer perceptron neural network (MLPNN) and adaptive neuro-fuzzy inference system (ANFIS). These models first use landslide and non-landslide records for comprehending the relationship between the landslide occurrence and predisposing factors. Then, landslide susceptibility values are predicted for the whole area. The accuracy of the produced susceptibility maps is measured using area under the curve (AUC) index, according to which, the MLPNN (AUC = 0.916) presented the most accurate map, followed by the ANFIS (AUC = 0.889) and FR (AUC = 0.888). Visual interpretation of the susceptibility maps, FR-based correlation analysis, as well as the importance assessment of predisposing factors, all indicated the significant contribution of the road networks to the crucial susceptibility of landslide. Lastly, an explicit predictive formula is extracted from the implemented MLPNN model for a convenient approximation of landslide susceptibility value.


Introduction
Landslide is a ubiquitous environmental disaster that is responsible for nearly 17% of life losses and exclusively accounts for about 5% of the natural catastrophes worldwide (Kjekstad and Highland 2009;Pham et al. 2019). Different definitions and classification systems have been proposed for landslides (Li and Mo 2019). Cruden (1991), for instance, 1 3 explained landslide as mass movements of rock, earth or debris down a slope. It may not be feasible to stop or control the landslide. Nevertheless, utilizing state-of-the-art technologies like remote sensing (Scaioni et al. 2014) and decision support systems for recognizing landslide-prone areas is an effective way toward mitigating the corresponding damages (Thai Pham et al. 2019).
Landslide is a dynamic phenomenon (Pandey et al. 2021) whose susceptibility assessment requires taking the effect of several parameters into consideration (Hua et al. 2021). Due to this reason, experts have tested various methodologies such as statistical analysis (Reichenbach et al. 2018), decision-making (Yoshimatsu and Abe 2006), and soft computing (Huang and Zhao 2018) for susceptibility (and hazard) assessment of landslide (Catani et al. 2005). Tonini et al. (2020), for instance, presented a robust statistical method for spatiotemporal analysis of inventory maps toward a better visualization and quantitative assessment of recent Italian landslides. Yalcin (2008) used three models including analytical hierarchy process (AHP), statistical index (SI), and weighting factor for landslide susceptibility mapping at Ardesen Region, Turkey. Regmi et al. (2014) examined the efficiency of frequency ratio (FR), SI, and weights-of-evidence models for mapping the susceptibility of landslide in central regions of Nepal Himalaya. It was shown that the FR with a success rate and predictive accuracy of 76.8% and 75.4%, respectively, performed better than two other models. The applicability of decision-making models has been investigated in various studies (Maqsoom et al. 2021;Mirdda et al. 2020). Pourghasemi et al. (2012) produced the landslide susceptibility map of the Safarood Basin, located in the Northern part of Iran, using index of entropy and conditional probability techniques. Both models presented higher than 82% accuracy which indicates the reliability of the maps.
Machine learning models have greatly assisted engineers in coping with many nonlinear problems, including the evaluation of environmental phenomena like snow avalanches , forest fire , flood (Costache et al. 2020), groundwater potential (Naghibi et al. 2016), gully erosion (Gayen et al. 2019), etc. Recent decades have witnessed the increasing popularity of various machine learning models. Logistic regression (LR) (Ayalew and Yamagishi 2005), boosting algorithms (Sahin 2020), support vector machine (Pourghasemi et al. 2013), tree-based approaches , and neighbor-based methods (El-Magd et al. 2021) are some of the well-known machine learning models that have been effectively used for landslide-related analysis, particularly susceptibility assessments.
Artificial neural network (ANN) and adaptive neuro-fuzzy inference system (ANFIS) are other notions of soft computing that have powerfully served for exploring nonlinear engineering parameters. These models use sophisticated algorithms to dig and learn the pattern of the intended parameters.  and Lee et al. (2003) tested and professed the feasibility of ANNs for landslide susceptibility mapping. Similar evaluations have been reported by Oh and Pradhan (2011) and Vahidnia et al. (2010) for the ANFIS. Jacinth Jennifer and Saravanan (2021) manifested a successful application of ANNs for susceptibility estimation of landslide in Idukki District of India. It was observed that the ANN models with a more complex structure achieved above 90% accuracy of prediction. Further applications of the MLPNN model for predicting environmental threats can be found in (Avand et al. 2021;Mohajane et al. 2021;Yi et al. 2020).
Hybridizing machine learning models with metaheuristic optimization algorithms has led to developing sophisticated techniques for landslide-related simulations (Moayedi et al. 2019;Nguyen et al. 2021). Panahi et al. (2020) optimized two conventional models, namely ANFIS and support vector regression using two metaheuristic techniques of bee algorithm (Bee) and the gray wolf optimizer (GWO). A comparison between the conventional and improved models revealed the higher accuracy of the hybrid tools. Mehrabi et al. (2020) trained an ANFIS model using four metaheuristic approaches for landslide susceptibility modeling in Qazvin County of Iran. They showed that genetic algorithm (GA) is a suitable algorithm for optimizing the parameters of membership functions in ANFIS. In research by Hong et al. (2020), the GA showed nice performance also for supervising multilayer perceptron neural network (MLPNN).
Comparative studies have always released valuable findings toward a convenient selection of landslide predictive models (Nhu et al. 2020;Zhu et al. 2018). In many works, authors have declared the better efficiency of soft computing models compared to traditional statistical tools. Park et al. (2013) compared the performance of ANN with FR, AHP, and LR for landslide susceptibility mapping in Inje Region, Korea. According to the respective values of area under the curve (AUC) 0.789, 0.794, 0.794, and 0.806 obtained for the AHP, FR, LR, and ANN, the superior accuracy of the ANN was deduced. A similar effort and conclusion were reported by Yilmaz (2009) for a case study from Kat landslides in Tokat City of Turkey. Sadighi et al. (2020) showed the ANFIS model outperforms ANN for landslide susceptibility modeling at Tajan Watershed, Northern Iran. The AUC values were 0.902 and 0.866. However, a hybrid ANFIS (coupled with imperialist competitive algorithm) with an AUC equal to 0.966 was found to be superior over both regular ANFIS and ANN. Lucchese et al. (2021) employed and compared the ANN with ANFIS for the same purpose in Rolante River Basin, Brazil. Based on the calculated AUCs, 0.8886 for ANFIS and 0.9409 for ANN, the later model could achieve a considerably larger accuracy.
The present research aims to produce applicable landslide susceptibility maps for the northern part of Lecco Province, Italy. This country, owing to its relief and lithological and structural features, is characterized by particularly high landslide risk (Trigila and Iadanza 2008). Having a look at the existing literature, while some studies have dealt with landslide susceptibility prediction in different parts of the Lombardy Region (Fabbri and Patera 2021;Sterlacchini et al. 2011), Lecco demands to receive proper analysis for mitigating the risk of this catastrophe.
Three different methodologies consisting of one traditional statistical method, namely the FR, and two artificial intelligence models, namely MLPNN and ANFIS, are employed in this research to predict landslide susceptibility value (LSV) all over the area. The above studies have professed the high efficiency of these three models (i.e., FR, MLPNN, and ANFIS) in the field of landslide susceptibility assessment. In the following, the case study and used database are introduced and spatial interactions are investigated. Then, the models are executed to produce and interpret the susceptibility maps. It is followed by accuracy assessment along with comparative validation of the results, and finally, the study ends up with introducing a neural-based LSV predictive formula.

Study area
The study area is the northern part of the Lecco Province located in the Lombardy Region, Italy. Figure 1 shows the exact location of this area. It lies within the longitude 09° 15′ to 09° 32′ E and latitude 45° 50′ to 46° 09′ N adjacent to the Como Lake. The area is roughly 473 km 2 and referring to the 2007 census, the province of Lecco has nearly 340,000 inhabitants (Parente et al. 2013). This region is under a warm and temperate climate with an annual rainfall around 1360 mm. According to online meteorological data (www. en. clima te-data. org; www. world weath eronl ine. com), the average monthly temperature ranges from 0° C in January to 19.1 °C in August. In the selected part, the altitude varies from 197.3 to 2608.6 m. The range of slope starts from 0 and peaks at around 87° close to the summits. According to the soil map, Cambisol is the dominant soil type covering around 45% of the area. Geologically, the area includes 43 units, out of which, the largest coverage is reported for Limestones and Dolomite. Considering the utilization of land, most of the area is classified as forest and seminatural.

Landslide inventory map
Inventory maps are essential prerequisites for spatial analysis of environmental threats like landslides (Can et al. 2019;Varnes 1984). They illustrate the distribution of the past events within the area of interest and contain also information regarding the date, outlines, and characteristics of occurrence Singh and Kumar 2018). The inventory map used in this study merges the information of two published inventory maps. As shown in Fig. 1, a total of 92 landslides are identified in the geographic information system (GIS). Out of these events, 82 landslides (18 areal events + 64 single events) are taken from dataset 1, the inventory map prepared by Calvello and Pecoraro (2018), and 10 landslides are taken from dataset 2, the database provided by the Inventario dei Fenomeni Franosi in Italia (IFFI https:// www. geopo rtale. regio ne. lomba rdia. it/ en/ metad ati?p_ p_ id= detai lShee tMeta data_ WAR_ gptme tadat aport let&p_ p_ lifec ycle= 0&p_ p_ state= norma l&p_ p_ mode= Fig. 1 The study area and distribution of the landslides view&_ detai lShee tMeta data_ WAR_ gptme tadat aport let_ uuid=% 7B1D4 AAE9F-EB7B-4E3E-AB8A-EE29A 22115 93% 7D).
In dataset 1, single events refer to the records of only one landslide, while areal events refer to the records of multiple landslides that were triggered by the same cause in the same geographic area. Extreme weather conditions have triggered most of these multiple landslides in one or more days. The news typically marks the affected area and focuses on the most severe landslide event(s). In this sense, rare reports can be found to completely represent an areal landslide (Calvello and Pecoraro 2018).
The landslides taken from dataset 1 occurred between early 2010 and late 2019, while the events taken from dataset 2 took place between 2000 and 2011. Considering the severity, the landslides are mostly classified as C2: severe events with injured persons and/ or evacuations, and C3: minor events which did not cause any physical harm to people (Calvello and Pecoraro 2018).
To create the spatial database of this study, a total of 92 non-landslide points are randomly generated within the areas without landslides. To do so, the pixels in the very close vicinity of the landslide points were excluded (using appropriate buffers), and afterward, the remaining area which was surely devoid of historical landslides was subjected to creating 92 random points. A few random distributions were tried to attain a reasonable nonlandslide map. Producing the same number of non-landslide points (relative to landslide points) is a well-accepted approach for creating the database of landslide, as well as similar phenomena like flood and wildfire (Avand et al. 2021;Avand and Moradi 2021;Moayedi et al. 2020;Yousefi et al. 2020).

Predisposing factors
The susceptibility of landslide is a function of several parameters that can affect the occurrence of this phenomenon. Therefore, proper selection of these parameters is of great importance (Dou et al. 2020;Li et al. 2021). Moreover, the reliability of data is another key parameter as it contributes to the quality of data (Mandal et al. 2021). In the present study, thirteen landslide predisposing factors, namely elevation, slope aspect, slope degree, plan curvature, profile curvature, distance to waterway, distance to road, distance to fault, soil type, land use, lithology, stream power index (SPI), and topographic wetness index (TWI), are taken into the equation to predict the LSV. For preparing these layers, the essential required layers were digital elevation model (DEM), the shapefile of linear phenomena (i.e., waterways, roads, and faults), as well as the shapefiles of soil type, land use, and geology map. All mentioned layers were downloaded from the website of Territorial Information of Region Lombardy (Geoportal of the Lombardy Region: (www. geopo rtale. regio ne. lomba rdia. it/ en/ home)) whose data is under either an IODL 2.0 https:// www. dati. gov. it/ iodl/2. 0/ or CC-BY-NC-SA 3.0 https:// creat iveco mmons. org/ licen ses/ by-nc-sa/3. 0/ it/ deed. it licenses. The GIS layers were then clipped for the area of interest for data processing and subsequent analysis. According to the source metadata, the DEM layer has been provided with 5 × 5 m spatial resolution using various resources such as local topographical data, 1 × 1 m resolution Lidar surveys, and the former edition of the 20 × 20 m regional DEM (metadata https:// www. geopo rtale. regio ne. lomba rdia. it/ en/ metad ati?p_ p_ id= detai lShee tMeta data_ WAR_ gptme tadat aport let&p_ p_ lifec ycle= 0&p_ p_ state= norma l&p_ p_ mode= 1 3 view&_ detai lShee tMeta data_ WAR_ gptme tadat aport let_ uuid=% 7BFC0 6681A-2403-481F-B6FE-5F952 DD48B AF% 7D).
Apart from the elevation layer that is represented by the DEM, the layers of slope aspect, slope degree, plan curvature, and profile curvature were directly produced from the DEM. Figure 2a shows the elevation map. The altitude values range from 197.3 to 2608.6 m which were classified into six groups including < 200, (200-700), (700-1200), (1200-1700), (1700-2200), and > 2200 m. The slope aspect, which illustrates the direction of slope face, has local influences on inducing instabilities of slopes (Chawla et al. 2019). For instance, two aspects can be in different situations regarding the effect of wind (e.g., temperature), rainfall (e.g., humidity and loosening of the soil), and sunlight (e.g., vegetation and melting snow). Figure 2b shows the slope aspect. Based on the GIS classification, this layer has the following groups: Flat, North, North-East, East, South-East, South, South-West, West, and North-West. Famously, the slope is one of the most relevant factors in landslide susceptibility assessment which demonstrates the rate of change in altitude (Mathew et al. 2009;Mokarram and Zarei 2018). The produced slope map is presented in Fig. 2c. Based on this map, the gentlest and steepest terrains are represented by the slopes of 0° and 86.7°, respectively. This layer was classified into < 15, (15-25), (25-35), (35-45), and > 45° (Tangestani 2009). Plan curvature is an indicator of flow acceleration and erosion/deposition rate and profile curvature can impact the variation of flow velocity down the slope (Kalantar et al. 2018). Understanding the combination of plan and profile curvature is substantial in analyzing the probability of landslide because it delineates the behavior of flow (Moore et al. 1991). Distance to linear features (i.e., waterways, road, and fault) has been among the most crucial predisposing factors for landslide susceptibility modeling that have been regarded in many previous efforts (Ozdemir and Altural 2013;Razavi-Termeh et al. 2021;. In general, the proximity to watercourses leads to increasing the erosion in the slope toes and also increasing the saturation of materials (Zhao et al. 2019). Higher groundwater and pore water pressure are other related consequences that may cause a higher likelihood of landslide near waterways. Road construction is a complicated task that changes the topography of the area. In many cases, this operation is associated with physical deformations in the slopes, and a larger probability of slope failure. Generally, distance to the faults is representative of tectonic activities, and the degree of weathering and fracturing in rock (Ali et al. 2021). Considering the drop of shear strength in fractured textures, this factor is inversely proportional to the probability of landslide. Figures 2f-h depicts the maps of distance to waterways, distance to road, and distance to fault, respectively. The distribution of previous events indicates that the majority of landslides have occurred along the named linear phenomena. The distance classes of (0-50), (50-100), (100-150), and > 150 m were applied for the waterways and roads, while the map of distance to fault was grouped into (0-200), (200-400), (400-700), (700-1000), and > 1000 m (Tangestani 2009).
Soil type is another important parameter that contributes to the occurrence of landslides through characteristics like permeability and cohesiveness (Avtar et al. 2011;Mandal et al. 2018). Seven soil categories are detected in the soil type map which is presented in Fig. 2i. These categories are originally Cambisols, Regosols, Fluvisols, Umbrisols, Cambisols podzolici, Leptosols, Luvisols that are represented by A, B, …, G in the corresponding map, respectively. The land use map was cropped from the CORINE Land Cover (CLC) layer. In the CLC legend, the utilizations are categorized into five major classes: (a) artificial surfaces, (b) agricultural areas, (c) forest and seminatural areas, (d) water lands, and (e) water bodies. Each major class is further divided into two levels of classification. Consequently, a three-digit number indicates the land use of each subclass (e.g., 312 stands for the land use category that Fig. 2 The map of predisposing factor: a elevation, b slope aspect, c slope degree, d plan curvature, e profile curvature, f distance to waterway, g distance to road, h distance to fault, i soil type, j land use, k lithology, l SPI, and m TWI 1 3 is the 2nd subclass belonging to the 1st class of the 3rd major category). More details regarding this layer can be found in (www. geopo rtale. regio ne. lomba rdia. it/ en; www. land. coper nicus. eu). Figure 2j displays the land use layer for the intended area. This map comprises 19 classes which are detailed in Table 1.
Lithology is a crucial parameter in landslide-related assessments as it can affect the formation and evolution of landslides, as well as the type and scale of this phenomenon Xiao et al. 2019;Yalcin et al. 2011). Based on the lithology map (scale: 1:250,000) shown in Fig. 2, a total of 43 geological units can be found in this area. As explained in Table 2, these layers are named A, B,…, Z, AA, AB, …, AQ.
The SPI and TWI are two important hydrological parameters that are frequently used in landslide-related analysis (Nsengiyumva et al. 2019;. The SPI denotes the erosive potential of the streams (Kumar and Anbalagan 2016) and TWI is an indicator of soil moisture contents that contribute to the occurrence of landslide (Pokharel et al. 2021). Given β as the steepness of terrain and SCA as specific catchment area, Eqs. 1 and 2 express the formulation of these predisposing factors (Ghasemain et al. 2020):

Correlation assessment
For quantifying the relationship between the landslides and predisposing factors, the FR model is employed. Since the FR is one of the models used in this study, it is well explained in the next section. In this method, based on statistical analysis, one value is calculated for each subclass which determines its correlation with the landslide occurrence. Thus, the larger the FR is, the bigger the contribution of the subclass is (Huang et al. 2021). Also, the FR values equal to 1 reflect an average correlation, while those below 1 and above 1 indicate a lower and higher correlation, respectively (Akgün and Bulut 2007).
All layers were prepared in the ArcMap with a cell size of 5 × 5 m. It resulted in layers with 18,908,558 pixels. The FR values were calculated by intersecting each layer with the whole landslide points. Table 3 gives the results of this process. In the elevation map, the FR values of the first two classes are meaningfully greater than others. In the aspect layer, the FR of South, South-West, and West is larger than 1. Considering the slope map, half of the landslides are located in slopes < 15° which cover approximately 16% of the area. It resulted in an FR value of 3.05. As for the plan curvature and profile curvature layers, the biggest FRs (i.e., 4.15 and 9.43, respectively) are calculated for the fifth and first subclasses, respectively. Concerning the linear features, while large FRs are obtained for pixels far from the waterways and faults, the pixels close to roads have significantly larger FRs. The biggest FR in this regard is 6.28 for pixels which are at maximum 50 m off the roads. The subclasses labeled as C from the soil type layer (FR = 3.62), O from the land use layer (FR = 9.88), and AG from the lithology layer (FR = 5.69) are characterized by the largest correlation with the landslides. The smallest SPI subclass and the largest TWI subclass acquired the highest FR which were 1.01 and 2.16, respectively.  Here, one may argue that some of these correlation results are not in consistent agreement with theoretical expectations. For example, while steep terrains are generally more prone to landslide, the FR values obtained for the low-slope classes were comparably larger. The same can be raised for places far from waterways which were a little more correlated with landslide occurrence compared to nearby areas. In response, it is well known that the occurrence of landslide is a result of interrelated effects. In other words, the association of several predisposing factors determines the susceptibility of a region. Therefore, in the above example, gentle slopes might have experienced more landslides under the effect of other parameters (e.g., lying on fractured rocks or change in the land use). This is why an importance analysis is performed in the following sections to elucidate the contribution of the predisposing factors.

Methodology
The scheme taken for fulfilling the objective of this study is depicted in Fig. 3. After proper provision of landslide predisposing factors and inventory map, the FR model delineates the spatial correlation between the landslide and predisposing factors. The database is then divided into the training and testing subsets. Utilizing the training set, the FR, MLPNN, and ANFIS models are executed to calculate the LSVs over the study area in order to produce the susceptibility maps. The accuracy of the produced maps is evaluated with the help of testing points. Following this, a comparison points out the most accurate model. Finally, a neural-based equation is derived from the MLPNN model to be used for conveniently approximating the LSV.
The mechanism of the employed models (i.e., FR, MLPNN, and ANFIS) is explained in the following:

Frequency ratio
The FR model is a simple bivariate statistical approach that enables the user to acquire a quantitative representation from the spatial relationship between the landslide and predisposing factors (Termeh et al. 2018). It is a broadly used tool for probabilistic assessments of natural hazards in which multi-classified maps are involved (Bonham-Carter 1994).
In the FR method, each subclass is distinguished by a weight, for calculating which, two questions should be regarded. Assuming the landslide susceptibility problem, the questions are: (a) what percentage of the landslides are included in this subclass? and (b) what percentage of the whole area does this subclass cover? Eq. (3) can be written as follows:    in which l represent the number of landslides included in the subclass of interest, L is the number of all landslides, a symbolizes the area of the subclass of interest, and A is the whole area. The LSV of each pixel is finally calculated by summing the FR values obtained for all predisposing factors.

MLPNN
An MLPNN is a specific powerful type of ANNs that is distinguished by its layered structure (Hornik 1991). Generally speaking, ANNs are deemed as simulated biological neural networks which are capable of exploring complex engineering problems (Hornik et al. 1989). Provided a numerical dataset, the network uses part of the data for pattern recognition (i.e., training) and the rest is exposed to the obtained knowledge as testing data. Since the model has not met with the testing data before, the testing performance can represent the generalization power of the model. The components of an MLPNN are neurons (A.K.A nodes) that are connected by synapses. As Fig. 4 depicts the MLPNN used in this study, the connection is handled by weights (black arrows). These weights, as well as some bias terms (blue arrows), are tunable through the training process. Utilizing specific training schemes (e.g., backpropagation technique associated with the Levenberg-Marquardt (LM) algorithm (Moré 1978)), the neurons perform calculations in the form of Eq. 4 to establish a nonlinear dependency between the LSV and predisposing factors. where W and b signify the involved weight and bias, respectively. Also, PF is the predisposing factor and g() represents the activation function of the neuron.

ANFIS
An ANFIS is a hybrid tool composed of the intelligent computational strategy of the ANN and fuzzy logic. This model was designed by Jang (1993). Taking the advantage of if-then rules with respect to human experience, fuzzy logic attempts to map nonlinear complexities into scalar formats. The calculations in a fuzzy-based model draw on three major processes, namely fuzzification, fuzzy inference, and defuzzification. Based on this idea, crisp values are transformed into a linguistic fuzzy variable for feeding a FIS. The FIS then uses implication operations to apply so-called elements "fuzzy rules" to fuzzy variables. Lastly, the outcome of this process is converted into a crisp format again (i.e., defuzzification procedure) (Alajmi and Almeshal 2020). Equivalent to the weights and biases in an ANN, the parameters of the fuzzy membership functions are tunable items in the ANFIS.

Accuracy indicators
Plotting the receiving operating characteristic (ROC) curve, along with computing the area beneath it (i.e., the AUC index) is a recognized accuracy evaluation approach in such studies Zabihi et al. 2018). The ROC diagram draws specificity on the x-axis versus sensitivity on the y-axis. Sensitivity denotes the proportion of correctly classified landslide pixels, while specificity is expressed as the proportion of correctly classified nonlandslide pixels. Having TP, TN, FP, and FN as true positive, true negative, false positive, and false negative, respectively, specificity and sensitivity are calculated as follows (Chen et al. 2017;Hong et al. 2018;Pham et al. 2020): Fig. 4 The structure of the used MLPNN Note that the TP and TN stand for the numbers of correctly classified pixels, whereas FP and FN signify the numbers of erroneously classified pixels.
Moreover, two error indicators, namely mean square error (MSE) and mean absolute error (MAE) are used to calculate the error of prediction. Let LSV i expected . and LSV i modelled . represent the real (i.e., 0 and 1) and predicted LSVs, respectively, the MSE and MAE are defined as follows: in which the number of landslide points is represented by Z (which equals 128 for the training dataset and 56 for the testing dataset).

Results
To fulfill the purpose of the study, statistical and intelligent models are used to analyze the relationship between the landslide and its predisposing factors, and subsequently, predict the susceptibility for the area of interest. The results of the research are presented in this section. First, an importance analysis is carried out to evaluate the contribution of each predisposing factor. Next, the implementation of the models is explained, and after producing the landslide susceptibility maps, the accuracy of the models is assessed and compared.

Importance analysis
Importance assessment is a significant task in landslide-related simulations as it investigates the role of predisposing factors on the occurrence of this phenomenon. In this sense, statistical approaches are promisingly used for optimizing the problem through selecting the most contributing factors (Ciurleo et al. 2016). Applying some machine learning tools is another effective idea for importance assessment. Random forest (RF) (Breiman 1996) is one of these tools that can provide predictor importance estimation. In this work, a permutation-based importance approach is applied through training an RF (i.e., a bagged ensemble composed of 200 regression trees) in the MATLAB environment. This model recognizes independent factors as influential when the performance of the model drops as a consequence of permuting the values of the factors. In other words, the performance of the model shows little (or no) sensitivity to permutations applied to non-influential factors. The overall judgment is based on importance values (IVs) that are obtained by averaging the importance of predictors that originate from the same factor . Figure 5 shows the IVs obtained for predisposing factors in the form of column charts. It is clearly seen that distance to road (IV = 2.62) plays the most important role in this dataset. After that, elevation (IV = 0.56), slope degree (IV = 0.41), and slope aspect (IV = 0.38) have the largest contributions.

Model implementation
For all three models, the training data are used for acquiring the knowledge regarding the landslide pattern and the susceptibility maps are accordingly produced. More clearly, based on the existing events, the model conducts specific calculations to explore the landslide pattern and applies this knowledge to the whole area for producing the landslide susceptibility map.

FR
For producing the landslide susceptibility map using FR, the training events were intersected with the classified map of the predisposing factors to calculate the l ∕ L . and a ∕ A . ratios (see Eq. 3). After calculating the FR value for each subclass, the susceptibility map was obtained as the sum of all weight layers. It is worth noting that the FR values used for producing this map were different from those presented in Table 3, bacause the values in that table are calculated when all landslides are concerned.

Artificial intelligence models
Implementing intelligent models like the MLPNN and ANFIS entails providing appropriate numerical data to feed their networks. The values of thirteen predisposing factors were extracted to a total of 184 points (having 92 landslides and 92 non-landslides). Moreover, each point received a target value which was either 0 or 1 depending on being a nonlandslide or landslide point, respectively. After separating the 56 testing points, a training dataset composed of 128 samples was provided for the MLPNN and ANFIS. They explore

MLPNN
Due to the number of inputs (i.e., 13) and the unique target parameter, it is immediate that the proposed MLPNN should have 13 input neurons and one output neuron. In contrast, determining the number of hidden neurons is a challenging task in ANN models. Trial and effort, coupled with the user's experience, is a well-tried approach for properly determining this parameter. In this study, 50 MLPNNs with the architecture 13 × x × 1 (x = 1, 2, 3,…, 50) were considered where each one was executed 5 times for assessing the repeatability of results. This process revealed that three hidden neurons give the best results. Hence, an MLPNNs with the architecture 13 × 3 × 1 was chosen among a total of 250 tested networks. This model was trained by the LM algorithm which is among the most powerful techniques for this objective. In the training process, based on the influence of predisposing factors on the occurrence of landslides, the weights and biases of the MLPNN were adjusted. Almost in all cases, the divergence of the error (with 6 times tolerance) stopped the training process. Figure 6 shows the training results. In this figure, the predicted LSVs are compared to the expected ones (i.e., 0 and 1). For each sample, an error value gives the direct difference between these two values, and the frequency of these errors is analyzed in the form of a histogram chart. The calculated values of MSE (0.1115) and MAE (0.2623)

ANFIS
Similar to the MLPNN, there are some parameters that should be determined for achieving a reliable implementation of the ANFIS. The number of clusters in a fuzzy inference system can significantly affect the quality of learning. This parameter was tested to be 4, 5, 6, 7, and 8. (The values outside this range did not yield reasonable training.) Each model was repeated 10 times to check the repeatability of results.
Finally, out of 50 tested models, the ANFIS with six clusters was chosen due to the lowest training error provided by this architecture. Concerning other parameters of the ANFIS, the number of repetitions (i.e., epochs) was set to be 1000. A few greater values were tested as well, but the network experienced negligible changes. Between the hybrid and backpropagation optimization methods, the latter method was preferred based on a trialand-error effort. Figure 7 exhibits the training performance of the ANFIS. As is seen, the network has acquired a good understanding of the behavior of landslides. The responses of the network corresponding to the landslide data tend to 1, and likewise, those corresponding to the non-landslide data tend to 0. Based on the MAE and MSE values equaling 0.2278 and 0.0932, respectively, the training error is in an acceptable range. These errors are also smaller than the MLPNN model.

Susceptibility maps and interpretation
Three landslide susceptibility maps are prepared. As explained, the calculations of the FR model were carried out in the ArcGIS. But the procedure for producing the susceptibility maps of the MLPNN and ANFIS was different. For all pixels, the values corresponding to thirteen predisposing factors were converted to ASCII format. Then, they were given to the trained networks of MLPNN and ANFIS as new environmental conditions. The models predicted an LSV for each pixel and the results were imported back into ArcGIS to produce the susceptibility maps. The next step was classification of the maps to propose susceptibility levels within the studied area. To this end, the Natural Break classification technique was applied and yielded five categories characterizing very low, low, moderate, high, and very high susceptibility. It is worth mentioning that the Natural Break (A.K.A Jenks optimization method) is a well-tried classifier for the maps of natural hazards and other phenomena (Pourtaghi et al. 2015;Tehrany et al. 2019). For a certain number of classes, it aims to detect breaks that maximize between-class differences and minimizes within-class variance (Ahmed et al. 2021). Figure 8a-c presents the resulted maps. It can be seen that these maps are of good reliability, due to compatibility with the occurred landslides. The location of many previous events has been characterized by high or very high susceptibility. An appreciable point deduced from these maps is the high susceptibility of the northern part of the case study (around the dense road network) where no landslide has been reported. In contrast, some vast central and eastern regions are represented by green colors exhibiting low susceptibility.
To have a better perception of the susceptible areas, Fig. 8d shows very high susceptible places jointly identified by the FR, MLPNN, and ANFIS models. According to this map, as well as the 3D view shown in Fig. 8e, significant parts of the coastline of Lake Como (Lecco side) are deemed as crucially susceptible. For instance, a view of the municipality of Bellano is shown in Fig. 8f using Google Earth images.
Moreover, Fig. 9a illustrates very high susceptible areas predicted by the MLPNN model. In other components of this figure, Google Earth views are shown to magnify some specific parts. Figure 9b shows potential steep slopes above populated regions, Fig. 9c demonstrates open-pit mines exposed to very high susceptibility, Fig. 9d shows local and main roads at central parts which are under very high susceptibility, and Fig. 9e shows susceptible slopes above Moregge District and road tunnels there.
Radar diagrams in Fig. 10 report the percentage of the area that is covered by susceptibility classes. According to this figure, the majority of the area has been categorized as low and very low susceptible. The FR, MLP, and ANFIS have determined 4.52%, 12.6%, and 8.12% of the area as very high susceptible, respectively. Table 4 demonstrates the percentage of the training and testing landslides found in each susceptibility class. Around 48%, 66%, and 55% of the training points, and around 46%, 68%, and 57% of the testing point (respectively, by the FR, MLPNN, and ANFIS) are estimated to be in areas under very high susceptibility.
Back to Fig. 8, a susceptible line is detected in the low-height central area. In accordance with the map of predisposing factors, these areas are mostly overlaid with the road networks. The high contribution of the road networks was also inferred from the importance assessment (see Fig. 5) where the greatest IV was obtained for the factor distance to road. One reason for this high contribution can be anthropologically-induced instabilities, 1 3 i.e., mass movements initiated by road cutting operations (Nefeslioglu et al. 2008). It is known that building roads on the slopes lead to losing toe support, and consequently, increasing the strain behind the mass that causes the advent of cracks (Zhao et al. 2019). It may necessitate applying proper mitigation measures along the constructed roads, and likewise, performing specific investigations prior to upcoming road constructions.
A notable discrepancy with expectation is the high susceptibility of gentle terrains. Crossing Fig. 8d with the classified slope map revealed that the slope corresponding to 84.74% of the susceptible areas is below 15°. This observation can also be supported by  Table 3) where the biggest FR was observed for this group (< 15°) of slope layer. Hence, the areas with gentle slopes may receive equal, and even higher attention, compared to steep areas.
The same procedure was repeated for the layers of land use, soil type, and lithology. Areas distinguished with land use code 112 (i.e., discontinuous urban fabric with largest FR = 8.56) have 81.33% overlay with Fig. 8d. Likewise, Cambisols and Fluvisols are soil types that contain around 47 and 31% of the crucially susceptible areas, respectively. Concerning the geological units, approximately 46% of Fig. 8d falls into the K-labeled unit (i.e., Conoids) having an FR value of 3.73.

Validation and comparison
It was explained that this study employs the AUC as a common accuracy indicator for all three models. Figure 11 depicts the ROC diagram obtained for the (a) learning and (b) prediction phases of the FR, MLPNN, and ANFIS. Following the training results in Sect. 4.2, the training AUCs of the FR, MLPNN, and ANFIS were 0.928, 0.927, and 0.953, respectively. It can be inferred that the ANFIS comprehended the landslide behavior most accurately.
The landslide points that were selected as the testing data were used to examine the prediction accuracy of the maps. Having a glance at Fig. 11b, all three models have produced susceptibility maps with reliable accuracy. Based on the calculated AUC values, the line corresponding to the MLPNN has the largest area underneath. With 91.6% accuracy, the MLPNN was the most accurate model, followed by the FR and ANFIS with 88.8% and 88.9% accuracy, respectively. The prediction results are also compared in terms of sensitivity and specificity. The sensitivity and specificity of a model indicate true positive and false negative rates in assessing the classification accuracy, which in this case, respectively, correspond to classifying the landslide and non-landslide points . With 92.86% sensitivity, the performance of the ANFIS was superior to FR (sensitivity = 85.71%) and MLPNN   The calculated specificity value for the FR and ANFIS was 85.71%. Referring to the above comparison of the classification accuracies, the MLPNN could predict the susceptibility of landslide more successfully than the FR and ANFIS. Here, two intelligent models (i.e., the MLPNN and ANFIS) are also compared in terms of the MSE and MAE which reflect the error of prediction. The MSE and MAE of the ANFIS were 0.1425 and 0.2939, while these values were 0.1139 and 0.2720 for the MLPNN. Both criteria demonstrate that the neural model was stronger for predicting the landslide pattern in unseen environmental conditions. Table 5 summarizes the obtained accuracy indicators (MSE, MAE, and AUC associated with sensitivity and specificity).

An explicit LSV formula
In Sect. 3.2, the prediction mechanism of the MLPNN was mathematically described. Therefore, after training the model, its neural configuration exposes an explicit predictive formula that can be used for predicting the LSV. The MLPNN was structured as 13 × 3 × 1, denoting 13 input neurons, 3 hidden neurons, and 1 output neuron (see Fig. 4). In consequence, there are (13 × 3 =) 39 weights that connect the input and hidden neurons and (3 × 1 =) 3 weights that connect the hidden and output neurons. Also, each of the hidden and output neurons owns a bias term for creating their equations. Altogether, the formula is composed of 46 parameters. Equation 9 expresses the LSV formula developed in the unique output neuron.
where LO 1 , LO 2 , and LO 3 stand for the local outputs (i.e., the outputs of three hidden neurons). These values can be calculated by Eqs. 10, 11, and 12 wherein distance to waterway, distance to road, and distance to fault are abbreviated as DTW, DTR, and DTF, respectively.
(9) where Tansig() is the activation function that is employed by the hidden neurons for producing the local outputs (g() in Eq. 4). This function is expressed by Eq. 13.
This formula, owing to the desired accuracy of the MLPNN, can give a dependable approximation of the LSV in compatible areas. To apply this formula, the same input configuration (i.e., exactly these thirteen predisposing factors) should be used, and also, the classification of data should be in agreement with this study. For instance, class number 4 should be considered for Land use class D (i.e., mixed forest).
It should be regarded that, since the calculations are not computer-supervised here, Eqs. 10, 11, and 12 must be fed by normalized input data (e.g., within [− 1, 1]). Accordingly, the released LSVs can be different from the designated susceptibility extent which is 0 for not susceptible and 1 for quite susceptible. However, since susceptibility mapping requires relative LSVs, the outputs of this formula do not necessarily need to return into the unnormalized range. In other words, the map can be directly produced using the calculated values (similar to the LSVs obtained for the FR model).

Problem and solution
Susceptibility prediction is an effective solution for dealing with environmental threats (Sun et al. 2020;Wei et al. 2021). In the case of landslides, prediction-oriented efforts are of great value to engineers and decision-makers toward providing appropriate mitigation strategies and land use planning . For instance, authorities may seek to identify safe locations for the development of new habitats or facilities like roads, transmission, energy stations, etc. (Moradi et al. 2019). Up to now, a wide variety of modeling tools have been proposed to model the susceptibility of landslide all over the world. Concerning the methodologies that require prior landslide events as primary information, statisticalbased model and artificial intelligence are among the most popular ones.
This study employed the FR and two intelligent models (i.e., MLPNN and ANFIS) for landslide susceptibility prediction in the province of Lecco (Northern half). Using highresolution data led to producing detailed landslide susceptibility maps. Moreover, accuracy assessments revealed that all three maps can yield a reliable approximation of the susceptible areas. The AUC values corresponding to the prediction phase were 0.888, 0.916, and 0.889 that profess high robustness of all three models against the complexity and nonlinearity of the given problem.

Data and inventory map
The used inventory map was formed by single landslide events and some landslides referred to as areal events. As elucidated earlier (Sect. 2.2), an areal landslide may comprise more than one event that, due to the same triggering cause, have occurred in the same geographical area. Hence, the attributed coordinates are indicative and, by focusing on major landslides, they give an approximation of the affected region. Given the number and importance of these events, this portion of data was taken into the analysis as released by FraneItalia catalog. Furthermore, it is worth remembering that, in such analysis, the landslide predisposing factors are treated as classified data. More clearly, all pixels grouped in one class play an equal role in the occurrence of landslide (e.g., all pixels with the slopes between 15 and 25° fell into one class and received the same weight in terms of this particular predisposing factor). On the other hand, for a given areal landslide, multiple events that have the same triggering factor within identical geographical areas are very likely to fall into the same class as their major representative point. Thus, they would be identically treated. This reasoning, besides explanations in Sect. 2.2, may answer the question "How reasonable is in this work to represent areal landslides by a single point?" According to Calvello and Pecoraro (2018), data collection was associated with limitations like uncertainties in the number of areal events. For this reason, the only estimated areal event (i.e., the event characterized based on the operators' inference) in the study area was eliminated to ensure a completely certified dataset is exposed to the models. From a computational point of view, it helps to keep the analysis (e.g., weight assignments in the MLPNN model) safe from undesirable uncertainties. However, improving post-event investigations by related organizations can be an effective remedy to shortcomings like this.

A comparative evaluation
From a comparative viewpoint, spatial assessments using bivariate statistical approaches (e.g., FR) are associated with disadvantages like being time-consuming (Yilmaz 2010) and disregarding the interaction between variables (in determining the weight) (Ozdemir 2011). This is while both of these weaknesses are nicely overcome when artificial intelligence techniques are applied. They are capable of automatically tuning the variables in a very efficient way. Nonetheless, a significant difficulty in using intelligent models lies in the necessity of converting the GIS data into pure numerical formats like ASCII. Choosing appropriate hyperparameters is another issue that should be carefully taken care of (Youssef and Pourghasemi 2021). Determining the number of hidden neurons in MLPNN (or the number of clusters in ANFIS) is an example of this challenge that, in the present research, was solved by taking the advantage of the trial-and-error method.
Both MLPNN and ANFIS were trained using the backpropagation scheme in which the algorithm propagates backward for rectifying the parameters with respect to the calculated error. The backpropagation is a capable algorithm for landslide-related evaluations which has been successfully tried by authors like  and Wu et al. (2013). Although the ANFIS achieved a more accurate training (MSEs of 0.0932 vs. 0.1115 and AUCs of 0.953 vs. 0.927), the testing performance of the MLPNN was superior (MSEs of 0.1139 vs. 0.1425 and AUCs of 0.916 vs. 0.889). It can reflect the higher flexibility of the MLPNN when it is applied to unseen environmental conditions. Furthermore, the training process of the MLPNN took less than one second (0.0612 s), while the time elapsed for the ANFIS was 5.9285 s under an Intel Core i7 64-bit operating system with 16 gigs of RAM. It points out the MLPNN to be much faster than ANFIS in solving the given landslide susceptibility problem. Altogether, the MLPNN was introduced as the most efficient landslide predictive model, based on which, an explicit LSV predictive formula was extracted.

Comparison with literature
Based on the existing inventory maps, landslide is considered a frequent major natural disaster in Italy (Trigila et al. 2010). The produced susceptibility maps were in satisfying agreement with historical landslide events of the Lecco Province. Nonetheless, the outcomes are compared to the results of earlier studies conducted in different scales. From a regional point of view, Lombardy has experienced many landslides, due to which, many researchers have regarded landslide susceptibility modeling in different parts of this region (e.g., Valtellina Valley (Van Den Eeckhaut et al. 2012)). Yordanov and Brovelli (2020) produced the landslide susceptibility map of the Val Tartano District (Southern Sondrio) using the SI, LR, and RF approaches. Val Tartano is in the right-hand-side proximity of the area examined in the present work. They also investigated the effect of partitioning ratio for forming the training and testing databases. A 70/30 ratio was suggested as the most suitable proportion which is the one used in the present study. In a risk assessment effort conducted by Lari et al. (2009) over the Lombardy Region, significant parts of Lecco, particularly coastlines on both sides of the Lake Como and some central parts, fell into high, very high, and extremely high hydrogeological risk (i.e., the risk of landslides along with floods and snow avalanches) (see Fig. 6a of the cited paper). These hazardous parts have many places in common with landslide susceptible areas recognized in this study.
The findings of this research also compromise with similar studies carried out at national/ international levels. It was here discussed that, noticing the significant contribution of road networks and gentle terrains in the occurrence of landslides, the study area requires appropriate attention for risk mitigation. In the attention level map suggested by Trigila and Iadanza (2008), the northern part of Lombardy (including the present studied area) demands high and very high attention with respect to the risk of landslide and utilization of land (see Fig. 18 of the cited paper).
The results of Europe-wide landslide susceptibility assessment fulfilled by Van Den Eeckhaut et al. (2012) showed that the vicinity of the Alp Mountains is considerably more susceptible than the rest. Magnification on the obtained map (see Fig. 6d of the cited paper) illustrates that this vicinity includes the north of Italy, and above all, the province of Lecco.
Furthermore, Europe-wide and national-wide landslide susceptibility mappings carried out by Günther et al. (2008) and Trigila et al. (2013), respectively, are in partial accordance with this study. Günther et al. (2008) presented a probabilistic susceptibility map in which the probability values between 0.8 and 1.00 are most prevalent in the north of Lombardy (see Fig. 2 of the cited paper). Similar evaluations can be inferred from the landslide indices calculated by Trigila et al. (2013) (see Fig. 7 of the cited paper). Also, other global-scale analyses can validate this inference as well (Hong et al. 2007;Stanley and Kirschbaum 2017).
Once again, compared to the case of this study, the maps regarded for the above validation were drawn for much wider extents, i.e., national, continental, and global scales. They mostly provide information for decision-making in such levels (e.g., allocating budgets). Hereupon, they cannot be expected to reflect the same details as the local maps do. For instance, the exact distribution of susceptible/not susceptible properties and habitats needs local-scale investigations.

Future efforts
Concerning the results of this study, future projects are recommended to focus on a profound landslide risk analysis through assessing the hazard, exposure, and vulnerability over the study area. Proposing mitigation policies (e.g., monitoring (Hojat et al. 2019) and early warning systems (Pecoraro et al. 2019)) that may reduce the risk of landslide over the populated areas and valuable assets may be another viable subject.
Moreover, several ideas can be regarded to improve prediction efficiency. The betterment can happen in terms of accuracy and problem complexity. Employing optimization techniques for both training of the intelligent models and optimizing the number of predisposing factors would be worth trying. Comparative studies are also of high interest for recognizing more effective intelligent models.

Conclusions
Italy, and more particularly Lecco Province, is a landslide-prone area. In this research, reliable susceptibility maps were produced using three robust methodologies (the MLPNN and ANFIS as machine learning models and FR as a statistical model) and high-resolution spatial data. The validation of the results reflected good accuracy for all implemented models. However, the MLPNN emerged as a more generalizable predictive tool. The results illustrated a reasonable susceptibility zonation of landslide. Significant areas were jointly recognized as highly susceptible by all three models. There are valuable assets and populated areas falling into this level of susceptibility. Considering the contribution of predisposing factors, it was discussed that slope degree and distance 1 3 to road factors are expected to be highly regarded for taking mitigation measures. All in all, the findings of this paper can be interesting for mitigating the risk of landslide through land use planning and decision-making within the studied area.