Landslide Susceptibility Mapping of Landslides with Artificial Neural Networks: Multi-Approach Analysis of aBackpropagation Algorithm Applying the Neuralnet Package in Cuenca, Ecuador

Natural hazards generate disasters and huge losses in several aspects, with landslides being one of the natural risks that have caused great impacts worldwide. The aim of this research was to explore a method based on machine learning to evaluate susceptibility to rotational landslides in an area near Cuenca city, Ecuador, which has a high incidence of these phenomena, mainly due to its environmental conditions, and in which, however, such studies are scarce. The implemented method consisted of an artificial neural network multilayer perceptron (ANN MLP), generated with the neuralnet R package, with which, by means of different backpropagation algorithms (RPROP+, RPROP−, SLR, SAG, and Backprop), five landslide susceptibility maps (LSMs) were generated for the study area. A landslide inventory updated to 2019 and 10 conditioning factors, mainly topographical, geological, land cover, and hydrological, were considered. The results obtained, which were validated through the AUC-ROC value and statistical parameters of precision, recall, accuracy, and F-Score, showed a good degree of adjustment and an acceptable predictive capacity. The resulting maps showed that the area has mostly sectors of moderate, high, and very high susceptibility, whose landslide occurrence percentages vary between approximately 63% and 80%. In this research, different variants of the backpropagation algorithm were implemented to verify which one gave the best results. With the implementation of additional methodologies and correct zoning, future analyses could be developed, contributing to adequate territorial planning and better disaster risk management in the area.


Introduction
Natural hazards generate disasters and huge losses in several aspects [1,2]. These losses are linked to population growth, especially in urban areas, which have seen the increase of settlements and lifelines in hazardous areas and the propensity to suffer the effects of such hazards [3]. According to this, it is necessary to develop measures to mitigate the risk generated by these hazards, contributing to the development of Sustainable Development Goals (SDGs), such as SDG 11 "Sustainable cities and communities", which is related to the Sendai Framework, which promotes the reduction of disaster risk and losses caused by disasters [4]. Because hazards caused by natural phenomena occur spatially, their analysis must be spatial, which allows establishing methodologies such as the generation of susceptibility maps to certain risks [5].
Landslides are a hazard that cause major impacts worldwide [6,7]; they are considered destructive hazards due to their unpredictability and because their probability of occurrence generates better prediction results for the landslide susceptibility analysis in the study area, although it could also be applied in other areas. Recommended steps for a proper susceptibility analysis [21,57] were followed in the study zone, based on a landslide inventory and 10 conditioning factors. Finally, the performance of the implemented model was evaluated through a validation process and, with the use of GIS tools, corresponding LSMs were elaborated, which will be a useful element for adequate territorial planning and management.

Materials Study Area
The study area ( Figure 1) corresponds to the landslide inventory area (detailed below) and is located in a sector of Cuenca canton (Ecuador), located in Azuay province in the south of the country. It is a mountainous area of the Andes Mountains, with an elevation range that fluctuates between 2000 and 4500 m above sea level (m.a.s.l). This canton is one of the most important in the country; it has an area of approximately 3100 km 2 and a population of approximately 636,000 habitants [58]. In terms of topography, the canton has high areas and plains whose average altitude is 2853 m.a.s.l, specifically in the study area. Climatic conditions of the canton vary between dry periods (from June to November) and high levels of precipitation periods (from December to May) with an annual average of 940 mm [58]. Regarding interannual variability, the influence of climatic phenomena such as ENSO (El Niño Southern Oscillation) is apparently not significant on annual precipitation levels in the inter-Andean valleys of Ecuador [59]. In terms of geology, the area of Cuenca city corresponds to an intramountainous Andean basin. Materials consist of a volcanosedimentary series, with deposits of detrital materials, such as conglomerates, sands, silts, and clays, of Miocene-Pliocene age, interspersed with volcanic materials, such as andesitic and rhyolitic tuffs [40,58,60,61]. An important detail is that, in the study area, more than 75% of the landslides recorded in the inventory were rotational [55,58,59].
Cuenca city in Ecuador, evaluating different backpropagation algorithms (RPROP+, RPROP−, SAG, SLR, and Backprop) using the neuralnet package of the R language (free software), which allows the analysis of different approaches to an ANN multilayer perceptron (ANN-MLP). In this research, a different approach is applied in the sense that several backpropagation algorithms are analyzed using the package neuralnet. The studies reviewed have generally applied the backpropagation algorithm without analyzing its variants, which has been implemented in this research. With this, it is possible to verify which of these options generates better prediction results for the landslide susceptibility analysis in the study area, although it could also be applied in other areas. Recommended steps for a proper susceptibility analysis [21,57] were followed in the study zone, based on a landslide inventory and 10 conditioning factors. Finally, the performance of the implemented model was evaluated through a validation process and, with the use of GIS tools, corresponding LSMs were elaborated, which will be a useful element for adequate territorial planning and management.

Study Area
The study area ( Figure 1) corresponds to the landslide inventory area (detailed below) and is located in a sector of Cuenca canton (Ecuador), located in Azuay province in the south of the country. It is a mountainous area of the Andes Mountains, with an elevation range that fluctuates between 2000 and 4500 m above sea level (m.a.s.l). This canton is one of the most important in the country; it has an area of approximately 3100 km 2 and a population of approximately 636,000 habitants [58]. In terms of topography, the canton has high areas and plains whose average altitude is 2853 m.a.s.l, specifically in the study area. Climatic conditions of the canton vary between dry periods (from June to November) and high levels of precipitation periods (from December to May) with an annual average of 940 mm [58]. Regarding interannual variability, the influence of climatic phenomena such as ENSO (El Niño Southern Oscillation) is apparently not significant on annual precipitation levels in the inter-Andean valleys of Ecuador [59]. In terms of geology, the area of Cuenca city corresponds to an intramountainous Andean basin. Materials consist of a volcano-sedimentary series, with deposits of detrital materials, such as conglomerates, sands, silts, and clays, of Miocene-Pliocene age, interspersed with volcanic materials, such as andesitic and rhyolitic tuffs [40,58,60,61]. An important detail is that, in the study area, more than 75% of the landslides recorded in the inventory were rotational [55,58,59].

Methods
The methodology (Figure 2) was implemented through the following stages: (1) obtainment of the landslide inventory; (2) generation of the conditioning factors; (3) extraction of the training and test datasets; (4) implementation of the ANN-MLP algorithm and configuration of the hyperparameters; (5) results validation; and (6) generation of the landslide susceptibility maps.

Methods
The methodology (Figure 2) was implemented through the following stages: (1) obtainment of the landslide inventory; (2) generation of the conditioning factors; (3) extraction of the training and test datasets; (4) implementation of the ANN-MLP algorithm and configuration of the hyperparameters; (5) results validation; and (6) generation of the landslide susceptibility maps.

Obtaining the Landslide Inventory
The landslide inventory used in this research was elaborated in 2019 thanks to the joint work of the Institute for Sectional Regime Studies of Ecuador (IERSE) of the University of Azuay (Cuenca, Ecuador) with the University Federico II (Naples, Italy). The work was based on the landslide inventory conducted in 1994 by the PRECUPA project due to the "La Josefina" macro landslide. The inventory covered an area of 380 km 2 and was divided into 36 quadrants of 13.5 km 2 distributed equally between the northern and southern zones of the study area to facilitate the survey and verification of information. The inventory was based on fieldwork with the in situ presence of technical staff to obtain information. This task was carried out using the MARLI (Mobile Application for Regional Landslide Inventories) application [58]. Furthermore, it was necessary to perform photointerpretation tasks with satellite images (Planet) and orthophotos of the study area. Differential interferometry (DInSAR) was also applied with COSMO-Skymed (Constellation of Small Satellites for Mediterranean Observation) and Copernicus Sentinel-1 images to identify landslides in areas that were difficult access. It is important to mention that SAR data form a two-dimensional matrix containing phase and amplitude information related to backscattered electromagnetic radiation; therefore, SAR interferometry allows the extraction of the phase component [62]. The wide availability of SAR data, the improvement of processing techniques, and the reduction of satellite revisit times have allowed InSAR to be applied for the analysis of a wide range of ground motion phenomena [63]. In this sense, the specific InSAR application of this research (DInSAR), is detailed in the study of Miele et al. [61]; however, it is highlighted that the images used were COSMO-

Obtaining the Landslide Inventory
The landslide inventory used in this research was elaborated in 2019 thanks to the joint work of the Institute for Sectional Regime Studies of Ecuador (IERSE) of the University of Azuay (Cuenca, Ecuador) with the University Federico II (Naples, Italy). The work was based on the landslide inventory conducted in 1994 by the PRECUPA project due to the "La Josefina" macro landslide. The inventory covered an area of 380 km 2 and was divided into 36 quadrants of 13.5 km 2 distributed equally between the northern and southern zones of the study area to facilitate the survey and verification of information. The inventory was based on fieldwork with the in situ presence of technical staff to obtain information. This task was carried out using the MARLI (Mobile Application for Regional Landslide Inventories) application [58]. Furthermore, it was necessary to perform photointerpretation tasks with satellite images (Planet) and orthophotos of the study area. Differential interferometry (DInSAR) was also applied with COSMO-Skymed (Constellation of Small Satellites for Mediterranean Observation) and Copernicus Sentinel-1 images to identify landslides in areas that were difficult access. It is important to mention that SAR data form a two-dimensional matrix containing phase and amplitude information related to backscattered electromagnetic radiation; therefore, SAR interferometry allows the extraction of the phase component [62]. The wide availability of SAR data, the improvement of processing techniques, and the reduction of satellite revisit times have allowed InSAR to be applied for the analysis of a wide range of ground motion phenomena [63]. In this sense, the specific InSAR application of this research (DInSAR), is detailed in the study of Miele et al. [61]; however, it is highlighted that the images used were COSMO-Skymed and Sentinel 1A (upward mode) and B (downward mode), the latter being acquired between October 2016 and May 2019, and further processed using the coherent pixel technique-temporal phase coherence (CPT-TPC) approach. This process consisted of the selection of stable coherence scatterer (SCS) pixels, which were selected using the temporal phase coherence estimator (TPC). With this, the deformation of the selected pixels was evaluated and the results were geocoded. The use of InSAR has advantages since it allows measuring ground displacements with millimeter precision [61,64], ideal for detecting surface changes [62]. It allows easy monitoring of ground movements, considering the scale of analysis and intensity of deformation [63]. Furthermore, it has innovative applications that consist of combining them with machine learning algorithms, in which datasets are trained with InSAR data to improve susceptibility models [65] and have been used to assess the probability of occurrence of motion anomalies [63].
The inventory of landslides carried out with the MARLI application started from the need to generate a base inventory of landslides, since prior to this there was no inventory with the characteristics of precision, characterization, and completeness at acceptable levels. The only evidence of a landslide map is the one made in the PRECUPA project (1994-1998), which represents macro zones with a propensity for landslides [58].
The inventory update (2022) is at the project level to be delivered in November 2022, with the landslide risk management platform, which includes the landslide inventories and susceptibility maps generated, which will be published jointly between the University of Azuay (Cuenca-Ecuador), the Federico II University of Naples (Italy), the University of Jaen (Spain), and the National Risk Management Service (Ecuador). Based on the precept that every model is perfectible, with the update of the inventory, the field verification phase will allow the generation of more and new information to improve and validate these models, seeking the optimal factors of precision, completeness, and temporality, which are required for landslide risk management and governance [21,66].
In the study area, 710 landslides were recorded and classified according to their type of movement [11,12], being mainly rotational (529, approximately 75% of the recorded landslides), falls, rollovers, and flows, located mainly in the surroundings of the urban area of Cuenca city, with a landslide density per km 2 of 45.8 and an average landslide size of 2.2 ha. [58]. Since most recorded landslides were rotational, this research focused on analyzing them specifically, considering that each landslide typology usually occurs under different conditions [67,68] and that the factors affecting each typology may be different [69]. Figure 3 shows two examples of rotational landslides recorded in the study area.
pixel technique-temporal phase coherence (CPT-TPC) approach. This process consisted of the selection of stable coherence scatterer (SCS) pixels, which were selected using the temporal phase coherence estimator (TPC). With this, the deformation of the selected pixels was evaluated and the results were geocoded. The use of InSAR has advantages since it allows measuring ground displacements with millimeter precision [61,64], ideal for detecting surface changes [62]. It allows easy monitoring of ground movements, considering the scale of analysis and intensity of deformation [63]. Furthermore, it has innovative applications that consist of combining them with machine learning algorithms, in which datasets are trained with InSAR data to improve susceptibility models [65] and have been used to assess the probability of occurrence of motion anomalies [63].
The inventory of landslides carried out with the MARLI application started from the need to generate a base inventory of landslides, since prior to this there was no inventory with the characteristics of precision, characterization, and completeness at acceptable levels. The only evidence of a landslide map is the one made in the PRECUPA project (1994)(1995)(1996)(1997)(1998), which represents macro zones with a propensity for landslides [58].
The inventory update (2022) is at the project level to be delivered in November 2022, with the landslide risk management platform, which includes the landslide inventories and susceptibility maps generated, which will be published jointly between the University of Azuay (Cuenca-Ecuador), the Federico II University of Naples (Italy), the University of Jaen (Spain), and the National Risk Management Service (Ecuador). Based on the precept that every model is perfectible, with the update of the inventory, the field verification phase will allow the generation of more and new information to improve and validate these models, seeking the optimal factors of precision, completeness, and temporality, which are required for landslide risk management and governance [21,66].
In the study area, 710 landslides were recorded and classified according to their type of movement [11,12], being mainly rotational (529, approximately 75% of the recorded landslides), falls, rollovers, and flows, located mainly in the surroundings of the urban area of Cuenca city, with a landslide density per km 2 of 45.8 and an average landslide size of 2.2 ha. [58]. Since most recorded landslides were rotational, this research focused on analyzing them specifically, considering that each landslide typology usually occurs under different conditions [67,68] and that the factors affecting each typology may be different [69]. Figure 3 shows two examples of rotational landslides recorded in the study area.

Generation of Conditioning Factors
Adequately preparing the information that will allow the landslide susceptibility assessment is fundamental. Thus, information on landslides and conditioning factors of the

Generation of Conditioning Factors
Adequately preparing the information that will allow the landslide susceptibility assessment is fundamental. Thus, information on landslides and conditioning factors of the study area should be obtained, considering that the occurrence of future landslides based on certain conditions is likely to occur in areas where they have happened in the past [13] and that the selection of factors should be based on the characteristics of landslides and the environment where they occur [17]. The selection of the conditioning factors was based on the literature, analyzing systematic reviews [21,70] which, based on several publications related to LSMs, determined the most commonly used factors in this type of study. For this research, 10 conditioning factors were selected, divided into: (i) topographical factors (elevation, slope, aspect, curvature, stream power index (SPI), and topographic wetness index (TWI), derived from the digital elevation model (DEM)); (ii) geological factors (lithology, derived from the geological map); (iii) land cover factors (land cover, obtained from the land cover map; distance to road, obtained using a distance algorithm from the road network); and (iv) hydrological factors (distance to rivers, obtained from the hydrographic network of the study area) ( Table 1). The digital elevation model (DEM) used was obtained from the National Information System on Rural Land and Technological Infrastructure (SIGTIERRAS), with a spatial resolution of three meters and updated to 2010. Topographical factors were derived from this model using algorithms implemented in GIS software (QGIS). Furthermore, the land cover map was obtained from the same information system (SIGTIERRAS), with a scale of 1:25,000 and updated to 2018. Lithological information was obtained from the National Information System (SNI), with a scale of 1:100,000 and updated to 2005, while the information corresponding to road networks and hydrography was obtained from the Military Geographical Institute (IGM), with a scale of 1:25,000 and updated to 2013. A correlation analysis between factors was also performed in order to determine that the influence of these factors is high when implementing the model (this analysis is described in the Results section). Figure 4 shows the 10 conditioning factors used in this research and Table 2 shows additional information on them. It is important to mention that the ranges in this table were obtained based on the method of visualization of the rasters, with the features available for it in QGIS: Render type singleband pseudocolor, linear interpolation and equal interval mode. Table 3 shows some of the most important characteristics of the 10 conditioning factors used according to landslide susceptibility. interval mode. Table 3 shows some of the most important characteristics of the 10 conditioning factors used according to landslide susceptibility.  TWI, (f) curvature, (g) distance to roads, (h) distance to rivers, (i) lithology (lithology classes: 1sandy clays; 2-light laminated shales with gypsum; 3-silts, clays, sands, gravels, and blocks; 4heterogeneous mixture of fine materials and fine angular fragments without stratification; 5-heterogeneous mixture of fine materials and fine angular fragments (various sizes); 6-medium to coarse-grained tobaceous sandstones; 7-siltstones, shales, and fine-grained sandstones; 8-sandy silt, silty clay; 9-coarse andesitic conglomerates; 10-silts and clays; 11-red clays with sandstones and conglomerates alternation; 12-volcanic agglomerate; 13-sands, silts, clays, and conglomerates; 14-laminated mudstones, dark tobaceous sandstones; 15-variable proportion of silts, clays, sands, gravels, and blocks; 16-kaolinized tuffs and agglomerates; 17-massive siltstones), and (j) land cover.  TWI, (f) curvature, (g) distance to roads, (h) distance to rivers, (i) lithology (lithology classes: 1-sandy clays; 2-light laminated shales with gypsum; 3-silts, clays, sands, gravels, and blocks; 4-heterogeneous mixture of fine materials and fine angular fragments without stratification; 5heterogeneous mixture of fine materials and fine angular fragments (various sizes); 6-medium to coarse-grained tobaceous sandstones; 7-siltstones, shales, and fine-grained sandstones; 8-sandy silt, silty clay; 9-coarse andesitic conglomerates; 10-silts and clays; 11-red clays with sandstones and conglomerates alternation; 12-volcanic agglomerate; 13-sands, silts, clays, and conglomerates; 14-laminated mudstones, dark tobaceous sandstones; 15-variable proportion of silts, clays, sands, gravels, and blocks; 16-kaolinized tuffs and agglomerates; 17-massive siltstones), and (j) land cover.

Global Factor Conditioning Factor Data Range (Unit)
Topographical

Soil covers
Land cover

Global Factor Conditioning Factor Relevance
Topographical Elevation The probability of landslide occurrence is larger in areas where the elevation is higher [71].

Slope
Increasing slope decreases stability [55]. Its angle is considered as a controlling factor in landslide modeling [72].

Aspect
Indicates exposure to local climatic conditions [55] and atmospheric processes such as rainfall, wind, and solar radiation [73].

Curvature
Refers to a change in the slope gradient or aspect in a given direction [6]. Affects the control of water flow [33].

Stream power index (SPI)
Indicator of erosive processes caused by surface runoff. High SPI values indicate proximity to a stream. With low SPI values, a low susceptibility to landslide initiation is expected [55].

Topographic wetness index (TWI)
Indicator of saturated soil conditions during rainfall and sediment accumulation [55]. With a higher value of TWI, there is a greater tendency for the saturation of slope materials [73].

Soil covers
Land cover Each cover has characteristics and textures that influence landslide generation. It is also related to the degree of vegetation cover that influences the stability of slope materials [33].

Distance to roads
Related to the process of road construction, which, when developed in mountainous areas, causes impacts on slope stability [25].

Hydrological Distance to rivers
The closer this distance, the greater the probability of landslides [74]. It should be considered that a stream may be where landslides move to, which generates additional risk [75].

Extraction of Training and Test Datasets
Model performance should be evaluated using training and test datasets [35,76]. In this research, the dataset was divided into 70% training and 30% test, which is a common ratio [18] evidenced in several studies [28,29,32,40,47]. The main consideration for generating datasets in this research was to obtain points by random sampling in areas without landslides (labeled as 0) and areas with landslides (labeled as 1), based on rotational landslides and areas without landslides. Moreover, to generate the datasets, the numerical values of 10 conditioning factors were considered according to the location of random points in areas with and without landslides. The importance of these datasets lies in the possibility, with training data, to teach the model the classes to be predicted [24] and to verify its degree of fit [76], while with test data, model accuracy is verified [24], i.e., its predictive ability [76]. Susceptibility assessment should also consider the use of mapping units, which contain a set of specific terrain conditions [9] and can be grid cells, terrain units, topographic units, slope units, or unique-condition units. The grid cells method was used in this research.

Implementation of the ANN-MLP Algorithm
An artificial neural network (ANN) is a set of interconnected nodes useful for modeling problems with a complex relationship between analysis factors [35]. Its central processing unit is the neuron, which performs mathematical procedures to generate a result based on a set of input variables [77]. The application of an ANN in landslide susceptibility analysis is ideal because this phenomenon is dynamic and nonlinear [78]. The ANN architecture consists of a set of inputs (conditioning factors); a set of intermediate layers (hidden layers) that perform the processing; and an output layer [77] with the prediction result. Neural networks generally refer to supervised classification algorithms, which compare a given output with a predicted output, adapting the necessary parameters based on this comparison [79]. According to the literature, there are several neural network algorithms such as convolutional neural networks (CNNs), whose performance is adequate in visual image analysis [80], or recurrent neural networks (RNNs), useful when sequential information is available [81]; however, one of the most widely used is multilayer perceptron (MLP), which has been applied in several studies [31,33,35,82]. A perceptron is an individual neuron that allows classifying a set of inputs into one or two categories by means of a step function, which returns 1 if the weighted sum of inputs exceeds a threshold or otherwise returns 0 [77]. Another important feature is the activation function, which allows an ANN to work with nonlinear problems [77] and can be of linear, sigmoid, or logistic type (implemented in this research), hyperbolic tangent, or rectified linear unit (ReLU).
The MLP algorithm consists of a set of perceptrons organized in layers connected by synapses that are assigned a weight [79]. Connection weights between input layers, hidden layers, and the output layer are initialized and then updated using the backpropagation algorithm [35], whose basic operation consists of applying a chain rule to calculate the influence of each weight on the ANN with respect to an arbitrary error function E [83]: where w ij is the weight of neuron j to neuron i, S i is the output, and net i is the weighted sum of the inputs to neuron i. Once the partial derivative is calculated, the minimization of error function is achieved by simple gradient descent [83]: where is the learning rate, which relates to the time required for the network to reach convergence [83]. The use of hidden layers increases the model flexibility, though one hidden layer is usually sufficient to model continuous functions [84], considering that an excessive number of hidden neurons would cause overfitting problems; therefore, it is preferable to use the least number of hidden neurons for an adequate performance of the network [85]. Based on the above and considering that trial and error methods are used to determine the number of neurons in the hidden layer [86], the implementation of ANN in this research was based on an experimental procedure, elaborating several tests with different numbers of hidden neurons, which allowed obtaining a better performance of the network in terms of runtime, validation parameter values, and avoiding problems of non-convergence of the ANN, by configuring a hidden layer with three neurons. By increasing the number of hidden neurons or the number of neurons within this layer, problems such as an excessive runtime or lack of convergence were encountered.
ANN implementation in this research was performed using the neuralnet package of the R language. This package focuses on MLP [79] and allows the analysis of different implementation options of this algorithm. It is important to note that no research was found that obtained LSMs using different approaches to the MLP algorithm, since it is usually implemented without specifying the type of backpropagation algorithm used. Studies that applied neuralnet for classification problems (although not focused on landslide susceptibility analysis) with one hidden layer were also reviewed and obtained robust predictive performance [87,88].

Hyperparameter Settings
Hyperparameters are parameters that should be established before starting the ANN training process [89]. The main hyperparameter configured was the algorithm which allows the backpropagation algorithm to be chosen. It is possible to select between traditional backpropagation (Backprop); an algorithm with resilient backpropagation (RPROP), which simplifies the parameters needed for learning and can be RPROP+ (with weight backtracking) or RPROP− (without weight backtracking) [83]; or a modified globally convergent algorithm (GRPROP), based on RPROP− [90] and which is based on the modification of the learning rate SAG (smallest absolute gradient) or SLR (smallest learning rate) [90]. The remaining hyperparameters were selected according to the characteristics of the neuralnet package and were the same for each algorithm (Table 4). It is not possible to automate hyperparameters generation (tuning) when implementing an ANN as a classification problem, which, in this case, allows predicting zones with or without landslides, since implementing tuning is possible when using neuralnet only for regression problems [91]. Algorithms were implemented using RStudio Server and executed on a high-performance computer (HPC) from Corporación Ecuatoriana para el Desarrollo de la Investigación y la Academia (CEDIA).

Results Validation
The validation process is important for the model and derived susceptibility maps to have scientific meaning [26]. There are several methods for assessing both the fit of data and the predictive ability of the model. Commonly used methods to evaluate susceptibility models are the confusion matrix and area under the receiver operating characteristic (ROC) curve (AUC) [93,94]. The confusion matrix method is structured as true positive (TP), cases in which predicted and actual values are true; true negative (TN), cases in which predicted and actual values are false; false positive (FP), cases in which a value is predicted as true, but is false (Type 1 error); and false negative (FN), cases when a value is predicted as false but is true (Type 2 error) [77]. For susceptibility analysis, a threshold value is used as a basis for classifying the continuous susceptibility values into susceptible and non-susceptible. These values are compared with landslide distribution (presence or absence) and generate the matrix structure with TP, TN, FP, and FN values [95]. In this research, the threshold was defined was 0.5.
The AUC method allows the ROC curve to be quantitatively summarized [6] and thereby describes the ability of a model to correctly predict the occurrence or non-occurrence of an event [26], providing the model precision to predict landslide susceptibility [6] and indicating the effectiveness of the spatial forecast of the susceptibility model [94]. The AUC value usually varies between 0.5 and 1. A higher value indicates a better accuracy rate and hence reliability in the model [74,93,96]. A value less than 0.5 indicates that the prediction is a random guess [6]. Although the AUC value is a widely used method, implementing it as the only validation method to evaluate model performance is not recommended, since in some cases a high AUC value does not guarantee a high accuracy in spatial predictions [25]. For this reason, additional evaluation parameters were implemented in this research: sensitivity (recall), accuracy, positive predictive value (PPV or precision), and F-score, all obtained from the confusion matrix. Each parameter has particular characteristics, having in common that their values vary between 0 and 1, and was applied to the training and test data.
Sensitivity is the measure at which pixels corresponding to landslides are correctly classified as a landslide having occurred [35]. It is also called recall and its value is ideal as long as it is closer to 1 [77].
Accuracy is the ratio of landslide to non-landslide pixels that the model classifies correctly [35]. The closer its value is to 1, the better the performance the model presents [77].
Positive predictive value (PPV) refers to pixels that are correctly classified as landslides [35]. It is also known as precision and the closer its value is to 1, the more accurate the predictions will be [77].
Positive Predictive Value PPV (Precision) = TP FP + TP F-score (F1-score) is the harmonic mean between data classified as landslides and nonlandslides and provides a balanced view between them considering precision and recall values [99]; it is very useful when analyzing a scenario of uneven class distribution [8].

Obtaining Landslide Susceptibility Maps (LSMs)
At the end of ANN execution, the predicted values are obtained and it is possible to elaborate the susceptibility maps. For this purpose, it is necessary to establish the susceptibility levels based on these values, which were obtained using the quantile classification mode that defines the same number of elements in each class [100]. There were five susceptibility levels defined according to the classification method described above: very low; low; moderate (medium); high; and very high. LSMs were entirely elaborated with free GIS software (QGIS, GRASS, and SAGA), which was also applied in the implementation of the respective spatial analysis methods.

Correlation Analysis between Conditioning Factors
In this research, Spearman's correlation coefficient [55,101] was used to analyze the relationships between conditioning factors due to their non-normal distribution ( Figure 5).
The main advantage of this coefficient is that it is not affected by the data distribution [101]. Two factors are considered to be strongly correlated when the correlation coefficient exceeds 0.7 [102]. Since in the correlation matrix no factor exceeded the value of 0.7 (Table 5), all were considered for the model implementation, which has also been corroborated by analyzing the low levels of correlation between factors, which means that these factors have a high influence on the model [16].

Landslide Analysis According to Conditioning Factors
A rotational landslide occurrence analysis was performed according to the applied 10 conditioning factors. This analysis was performed using GIS tools in order to verify the frequency and density of landslides per square kilometer in relation to the conditioning factors ( Figure 6) and is similar to others performed previously [103,104], although in the present study all conditioning factors used were analyzed.
The classification of values of different classes for each conditioning factor was performed using the Pretty Breaks method, which generates sequences with "nice" intervals according to the data range [105] and is available in QGIS. Based on the range of values of each factor and according to the location of rotational landslides, the class intervals were regular, except for the initial and final classes. Moreover, according to the defined classes, a density analysis was performed to determine in which of them there was a higher occurrence of landslides per square kilometer.
Elevation: A total of 89.8% of rotational landslides in the study area occurred at a variable altitude between 2500 and 2900 m.a.s.l. The lowest number of events occurred between 2345 and 2400 m.a.s.l (0.2%). The highest landslide density value per km 2 was 2.49, located between 2800 and 2900 m.a.s.l.
Slope: Most of the rotational landslides (44.6%) occurred on slopes whose angle varied between 0 • and 10 • . Only 0.8% of landslides occurred on slopes between 40 • and 50 • and 0.2% occurred on slopes greater than 50 • . An important finding is that 54.4% of rotational landslides occurred on slopes between 10 • and 40 • . The highest landslide density value per km 2 was 1.62, located on slopes between 0 and 10 • .
Stream power index (SPI): A total of 49.5% of rotational landslides were located at SPI values less than 0 (specifically between −4.61 and 0). Only 1.5% of rotational landslides occurred at SPI values between 5 and 10. The highest landslide density value per km 2 was 1.77, located in negative SPI values (between −4.61 and 0).
Topographic wetness index (TWI): A total of 52% of rotational landslides were located at TWI values between 1 and 5. The smallest number of rotational landslides (0.6%) occurred at TWI values between 15 and 20. The highest landslide density value per km 2 was 1.92, located in TWI values between 15 and 20.
Curvature: Most rotational landslides (97.5%) occurred in a curvature range between −10 and 10. Specifically, 52.2% of landslides were recorded at curvatures between 0 and 10. The lowest number of events was recorded at curvature values above 10 (1.5%) and below −10 (1%). The highest landslide density value per km 2 was 2.67, located at curvature values between 10 and 16.8.
Distance to roads: A total of 95.7% of rotational landslides were located at distances less than 1500 m from a road network. Most of these (77.5%) were located at distances less than 500 m. A minimal number of rotational landslides (1.5%) occurred at distances between 2000 and 2500 m. The highest landslide density value per km 2 was 2.33, located at distances between 1000 and 1500 m.
Distance to rivers: A total of 88.7% of rotational landslides occurred at distances to a hydrographic stream less than 300 m. Most events (40.3%) occurred at distances between 0 and 100 m. The smallest number of rotational landslides (1.5%) occurred at distances to hydrographic streams greater than 500 m. The highest landslide density value per km 2 was 1.53, located at distances between 100 and 200 m.

Landslide Analysis According to Conditioning Factors
A rotational landslide occurrence analysis was performed according to the applied 10 conditioning factors. This analysis was performed using GIS tools in order to verify the frequency and density of landslides per square kilometer in relation to the conditioning factors ( Figure 6) and is similar to others performed previously [103,104], although in the present study all conditioning factors used were analyzed. The classification of values of different classes for each conditioning factor was performed using the Pretty Breaks method, which generates sequences with "nice" intervals according to the data range [105] and is available in QGIS. Based on the range of values of each factor and according to the location of rotational landslides, the class intervals were regular, except for the initial and final classes. Moreover, according to the defined classes, a density analysis was performed to determine in which of them there was a higher occurrence of landslides per square kilometer.
Elevation: A total of 89.8% of rotational landslides in the study area occurred at a variable altitude between 2500 and 2900 m.a.s.l. The lowest number of events occurred between 2345 and 2400 m.a.s.l (0.2%). The highest landslide density value per km 2 was 2.49, located between 2800 and 2900 m.a.s.l.
Slope: Most of the rotational landslides (44.6%) occurred on slopes whose angle varied between 0° and 10°. Only 0.8% of landslides occurred on slopes between 40° and 50° and 0.2% occurred on slopes greater than 50°. An important finding is that 54.4% of rotational landslides occurred on slopes between 10° and 40°. The highest landslide density value per km 2 was 1.62, located on slopes between 0 and 10°.
Stream power index (SPI): A total of 49.5% of rotational landslides were located at SPI values less than 0 (specifically between −4.61 and 0). Only 1.5% of rotational landslides Land cover: Most rotational landslides (32.1%) occurred in grasslands. The least number of rotational landslides occurred in herbaceous vegetation (2.5%), wasteland (0.9%), and anthropic infrastructures (0.6%). The highest landslide density value per km 2 was 1.46, corresponding to crops. Native forest and moorland covers did not register rotational landslides.
Lithology: Most rotational landslides (15.7%) were located in well-weathered thick and brecciated andesitic conglomerates, with sandstone intercalations, while the least number of events (5.2%) was located in silts, clays, sands, gravels, and blocks. The highest landslide density value per km 2 was 2.64, corresponding to clear laminated shales.

ANN Implementation and Performance
Once the high influence of all the conditioning factors in the model was verified, the neural network was implemented by configuring the necessary hyperparameters. Each of the applied algorithms presented different execution times (Table 6). The algorithm with the longest runtime was Backprop (560.7 s., approximately nine minutes) and it also presented the lowest performance according to the validation parameters explained below. The algorithm with the shortest runtime was RPROP− with 6.9 s., which was the one with the best results according to the ROC-AUC validation parameter with the test data.

Results Validation
The values obtained for AUC (Table 7) show that all the algorithms implemented in the model have a good fit with very similar values that presented a difference of 2.2%. The algorithm that presented the highest AUC value was SAG (0.889). As for test values (Figure 7), the results showed an acceptable performance and a difference of 5.4%, which in this research reflected a similar performance in the applied algorithms. The highest AUC value was obtained with the RPROP− algorithm (0.761); it is interesting to note that, without considering this algorithm, the difference between the remaining algorithms is minimal (0.7%). The values obtained for AUC (Table 7) show that all the algorithms implemented in the model have a good fit with very similar values that presented a difference of 2.2%. The algorithm that presented the highest AUC value was SAG (0.889). As for test values (Figure 7), the results showed an acceptable performance and a difference of 5.4%, which in this research reflected a similar performance in the applied algorithms. The highest AUC value was obtained with the RPROP− algorithm (0.761); it is interesting to note that, without considering this algorithm, the difference between the remaining algorithms is minimal (0.7%).  Statistical evaluation measures based on the confusion matrix with the training data ( Table 8) generally showed that the model had a good fit. When analyzing the PPV (precision) value, the highest value was obtained with SLR (0.996). For the sensitivity (recall) value, the highest value was obtained with Backprop (0.946). The RPROP− algorithm obtained the highest accuracy (0.85) and F-score (0.905) values. It is important to note that in PPV (precision) and sensitivity (recall) parameters the differences present a certain amplitude; in the first parameter there was a difference of 16.2%, while in the second there was a difference of 13.4%. The difference in accuracy was 3% and in F-score was 1.8%, which reflects a greater similarity when considering these parameters.  Statistical evaluation measures based on the confusion matrix with the training data ( Table 8) generally showed that the model had a good fit. When analyzing the PPV (precision) value, the highest value was obtained with SLR (0.996). For the sensitivity (recall) value, the highest value was obtained with Backprop (0.946). The RPROP− algorithm obtained the highest accuracy (0.85) and F-score (0.905) values. It is important to note that in PPV (precision) and sensitivity (recall) parameters the differences present a certain amplitude; in the first parameter there was a difference of 16.2%, while in the second there was a difference of 13.4%. The difference in accuracy was 3% and in F-score was 1.8%, which reflects a greater similarity when considering these parameters. Statistical evaluation measures with the test data (Table 9) showed that, based on PPV (precision), there was a difference of 18.8% between the values of each algorithm. The highest value was obtained with the SLR algorithm (0.964). For sensitivity (recall), the difference was smaller (5.7%) and the highest value was obtained with the Backprop algorithm (0.846). The difference in accuracy values was slightly smaller (5.4%); SLR obtained the highest value (0.775). For F-score, the difference was 5.9% and the highest value was obtained with the SLR algorithm (0.868).

Landslide Susceptibility Analysis
Due to the study area dimensions (X: 9427; Y: 6520), the total number of pixels in each map was 61,464,040. Since the study area contains several water bodies (mainly rivers), the corresponding pixels (19,448,146; 31.6%) whose value is constant were not assigned to any susceptibility level. With this consideration, the number of pixels of the study area (42,015,894) was obtained to determine the values corresponding to each susceptibility level, obtaining differences in the maps generated ( Figure 8). Based on the classification method used (quantiles), all algorithms obtained a similar percentage (around 20%) for each susceptibility level (Table 10), with the exception of SLR, which presented a high percentage at the moderate susceptibility level (44.79%) and a small percentage at the low level (3.29%). Furthermore, Backprop presented percentages of around 28% in very low and very high susceptibility levels. There was also a significant percentage of events in low and very low susceptibility areas, with SAG (36.1%), Backprop (35.3%), and RPROP− (34.4%). Figure 9 shows the percentage of rotational landslides assigned to each susceptibility level according to the implemented algorithms.

Discussion
The development of LSMs is a task that produces several uncertainties [13] and, in that context, the results obtained have shown that the level of fit and prediction of the implemented models is similar and acceptable according to the validation parameters used, with an average AUC value of 0.88 for the training dataset and 0.72 for the test dataset, which is similar to the results for the testing data obtained by Di Napoli et al. [16] in the specific implementation of an ANN (AUC = 0.74); Dou et al. [26] when applying an ANN evaluating fourteen conditioning factors (AUC = 0.73); and Pascale et al. [95] when implementing an ANN MLP with seven conditioning factors and one hidden layer (AUC = 0.75). However, the study by Mandal and Mondal [73], in which they implemented an ANN MLP, presented a superior performance to this research, considering its AUC value = 0.815. Regarding statistical evaluation based on a confusion matrix, the following values were also obtained on average: recall: 0.88 (training), 0.82 (test); accuracy: 0.84 (training), 0.75 (test); precision: 0.93 (training), 0.88 (test); and F-score: 0.90 (training), 0.84 (test). Bandara et al. [8], in their comparative study between different machine learning models for landslide susceptibility analysis, applied an ANN with precision = 0.82, recall = 0.61, and F-score = 0.70 for testing data, values which are lower than those obtained in this research. Similarity between the results described above and those obtained in this research reflects that, in the present study, the model implemented had a good fit and an acceptable predictive capacity, highlighting that, according to the highest value of AUC for the test data, the algorithm that seems to have the highest predictive accuracy when analyzing the susceptibility to rotational landslides in the study area is RPROP− (AUC = 0.761).
Considering other neural network models, Wang et al. [80], in their comparison of CNNs for LSMs in Yanshan

Discussion
The development of LSMs is a task that produces several uncertainties [13] and, in that context, the results obtained have shown that the level of fit and prediction of the implemented models is similar and acceptable according to the validation parameters used, with an average AUC value of 0.88 for the training dataset and 0.72 for the test dataset, which is similar to the results for the testing data obtained by Di Napoli et al. [16] in the specific implementation of an ANN (AUC = 0.74); Dou et al. [26] when applying an ANN evaluating fourteen conditioning factors (AUC = 0.73); and Pascale et al. [95] when implementing an ANN MLP with seven conditioning factors and one hidden layer (AUC = 0.75). However, the study by Mandal and Mondal [73], in which they implemented an ANN MLP, presented a superior performance to this research, considering its AUC value = 0.815. Regarding statistical evaluation based on a confusion matrix, the following values were also obtained on average: recall: 0.88 (training), 0.82 (test); accuracy: 0.84 (training), 0.75 (test); precision: 0.93 (training), 0.88 (test); and F-score: 0.90 (training), 0.84 (test). Bandara et al. [8], in their comparative study between different machine learning models for landslide susceptibility analysis, applied an ANN with precision = 0.82, recall = 0.61, and F-score = 0.70 for testing data, values which are lower than those obtained in this research. Similarity between the results described above and those obtained in this research reflects that, in the present study, the model implemented had a good fit and an acceptable predictive capacity, highlighting that, according to the highest value of AUC for the test data, the algorithm that seems to have the highest predictive accuracy when analyzing the susceptibility to rotational landslides in the study area is RPROP− (AUC = 0.761).
Considering other neural network models, Wang et al. [80], in their comparison of CNNs for LSMs in Yanshan (China) considering 16 conditioning factors, obtained different but similar AUC values for each applied algorithm, namely CNN-1D = 0.799; CNN-2D = 0.813; CNN-3D = 0.806; and LeNet-5 = 0.807. Although the performance obtained by CNNs is superior, it is not too far from that obtained in this research. On the other hand, Wang et al. [81], in their comparative study of RNNs for LSMs with 16 conditioning factors applied in Yongxin County (China), obtained better results than those obtained in this research according to the AUC value for each algorithm applied: RNN = 0.843; LSTM = 0.845; GRU = 0.860; and SRU = 0.838.
An important parameter that must be analyzed in the development of machine learning models is an algorithm's runtime. In this research, the implemented algorithms had notorious differences in the execution of the neural network, considering that the process was performed on an HPC. The RPROP+, RPROP−, and SAG algorithms had a very similar execution time (less than 10 s). RPROP− was the algorithm that had the fastest runtime with an approximate time of 7 s, while SLR and Backprop took longer compared to the aforementioned algorithms. Backprop took approximately nine minutes to run the ANN; however, its results in terms of AUC, accuracy, precision, and F-score values for the test data were the lowest performing. RPROP− is the algorithm that presented the best results in terms of AUC values for test data, and accuracy and F-score values for training data, in addition to presenting the shortest runtime. As they have a similar performance in terms of validation parameters, runtime may become an important parameter for choosing the algorithm with which the ANN could be implemented in future research.
Based on the test data and analyzing the different statistical validation parameters, it is interesting that the SLR algorithm stands out from most of them, even though its AUC (test) value was not the best (0.712). In this sense, the highest accuracy value was obtained by SLR (0.775), whose LSM presented a good percentage of zones with moderate susceptibility (44.79%). SLR also presented the highest precision (0.964) and F-score (0.868) values and showed a good coincidence between the points corresponding to rotational landslides and the zones of moderate and very high susceptibility, especially in the south of the study area. The above might suggest that the SLR algorithm is a good choice for generating LSMs. It is also interesting that Backprop presented the highest value for recall (0.846), but the lowest in the remaining parameters (accuracy = 0.721; precision = 0.776; F-score = 0.809). Moreover, the LSM generated with Backprop presented the highest percentage of zones with very high susceptibility (27.59%). These differences show the possible relationship between the validation parameters described above and the different percentages that each LSM presented with respect to landslide susceptibility levels.
In all LSMs, the highest susceptibility levels were observed to be "moderate" (SLR, SAG, and RPROP+) and "very low" (Backprop and RPROP−); however, all maps showed high and very high susceptibility levels in varying percentages (Figure 9). An interesting finding was that the performance according to the validation parameters can be classified as good in general, especially when analyzing the statistical parameters, and acceptable when analyzing the area under the ROC curve. When observing the LSM generated (Figure 8), it is evident that the part located southwest of the urban area of Cuenca shows moderate to very high susceptibility levels, with less coverage of low and very low susceptibility zones in all of them. Likewise, the great majority of the urban area presented a very low level of susceptibility on all maps, although there were some sectors of medium to very high susceptibility in the north of this area, since it has a higher elevation compared to the center and south of the city. Furthermore, it is interesting to note that SAG generated moderate coverages around the rivers that cross the study area and are notorious in the urban zone. In the southeast zone of the study area, all LSMs showed mostly moderate, high, and very high susceptibility levels, evidencing a good coincidence with the points corresponding to rotational landslides, especially in this zone and in the southwest zone.
The LSMs reflect that the study area has a higher proportion of sectors with moderate to very high susceptibility levels; however, it should be considered that this is a pioneering research, since no studies based on ANN-MLP or other machine learning algorithms for the generation of LSMs have been published in the study area, and although its findings will serve as a basis for the development of similar studies in the future, it is necessary to implement additional research to corroborate the results. In this context, it is suggested that future work should consider research to compare the performance between different models, such as classical (LR, logistic regression), frequently used in the development of LSMs (SVM or RF) and little used in this field (XGBoost); the implementation of methodologies such as principal component analysis (PCA) to determine the importance of conditioning factors; susceptibility analysis using higher values of spatial resolution, which has an important role in this regard [106,107]; the implementation of susceptibility analysis generating models with different conditioning factors or different numbers of them [26]; the use of other mapping unit methods, in order to verify the quality of the results obtained [17]; and analysis of different classification systems for LSM representation [108]. Sectors of zonification that present levels of susceptibility that imply risk should also be taken into account, in order to evaluate them exhaustively to determine the probability of occurrence of new events with a higher specificity level [109]. This can be addressed by specifically studying the sectors located southeast or southwest of the study area analyzed in this research or by focusing on the northern zone of the urban area of Cuenca city, which also presented areas with a variation in susceptibility level between moderate and very high.
Although the results obtained in this research, based on the validation parameters, indicate that the prediction levels are acceptable and serve as a starting point for susceptibility analysis in the study area, it is imperative to evaluate additional models due to the lack of similar research in this area. Furthermore, although the methodology implemented is solid, it is necessary to analyze other conditioning factors such as precipitation and to consider updating the landslide inventory to verify the variability of susceptibility levels. This opens up a research opportunity and therefore will generate valuable contributions, especially for studies comparing the performance of different models in order to determine the one that provides the best results, which would have a direct impact on the generation of useful elements that contribute to the formation of sustainable communities, since the generation of susceptibility maps should be oriented towards territorial interventions, providing reliable information about the probability that a phenomenon could negatively affect the population [109]. In this sense, the use of LSMs should be oriented towards their contribution in prevention decision making with the objective of mitigating the hazards inherent to landslides [26] and applying them in a practical way in land use planning or urban planning tasks [21]; however, their construction should be carried out carefully, based on the availability, quality, and accuracy of the data [32].

Conclusions
In this research, rotational landslide susceptibility maps of an area surrounding the city of Cuenca, Ecuador, were performed by applying different algorithms based on backpropagation (RPROP+, RPROP−, SLR, SAG, and Backprop) in an ANN-MLP implemented with the neuralnet R package and with the selection of 10 conditioning factors classified into topographical (elevation, slope, aspect, curvature, SPI, and TWI), geological (lithology), land cover (land cover and distance to roads), and hydrological (distance to rivers). This selection was based on a literature review; however, it is necessary to continue exploring different methods to determine the factors that contribute to the generation of more accurate prediction models. A landslide inventory in the area, updated to 2019, was also available, which is a fundamental element for the proposed evaluation. Conditioning factors were configured as input layers of the implemented ANN, which was executed based on training (70%) and test (30%) datasets to verify the fit and predictive capacity of the model, and with grid cells as the mapping unit. The results obtained were satisfactory since they showed a good fit and acceptable predictive capacity with average values of AUC = 0.72, precision = 0.88, accuracy = 0.75, and F-score = 0.84; however, analyzing the individual performance of the models based on the AUC value for the test data, the one that presented the best performance was RPROP− (0.761) and the one with the lowest performance was Backprop (0.707). However, for most statistical parameters, SLR performed well (accuracy = 0.775, precision = 0.964, and F-score = 0.868), which makes it an interesting option for the construction of LSMs. Since the models presented similar performances in terms of validation parameters, model runtime could also be considered for implementing ANN-MLP via neuralnet in future studies. In this sense, RPROP− was the algorithm that presented the fastest execution (6.9 s.), while Backprop took the longest to execute (approximately nine minutes). According to this, we can mention that all algorithms have presented acceptable results, except Backprop because of its lower performance in terms of validation parameters and its excessive runtime. The above evidence demonstrates the need for further research to corroborate these findings, mainly because this study is the first to apply a specific machine learning model for an evaluation of this type in the study area. It is also expected that the results obtained will be a good starting point to motivate further research of this type in the study area, especially in specific sectors, and considering additional methodological aspects such as implementing different machine learning models and comparing the results obtained; determining the importance of conditioning factors in the area; generating models with conditioning factors that present higher spatial resolution; or using different mapping units. The contribution of this type of research is of great importance for land use planning and management; however, it must be performed in the best possible way, considering above all the availability and quality of the data to be used. Funding: This research was funded by the University of Azuay (project No. 2020-0124) and by the project "Captura de Información Geográfica mediante sensores móviles redundantes de bajo coste. Aplicación a la gestión inteligente del territorio" (FEDER-UJA project No. 1265116).
Data Availability Statement: Not applicable.