Optimal slope units partitioning in landslide susceptibility mapping

ABSTRACT In landslide susceptibility modeling, the selection of the mapping units is a very relevant topic both in terms of geomorphological adequacy and suitability of the models and final maps. In this paper, a test to integrate pixels and slope units is presented. MARS (Multivariate Adaptive Regression Splines) modeling was applied to assess landslide susceptibility based on a 12 predictors and a 1608 cases database. A pixel-based model was prepared and the scores zoned into 10 different types of slope units, obtained by differently combining two half-basin (HB) and four landform classification (LCL) coverages. The predictive performance of the 10 models were then compared to select the best performing one, whose prediction image was finally modified to consider also the propagation stage. The results attest integrating HB with LCL as more performing than using simple HB classification, with a very limited loss in predictive performance with respect to the pixel-based model.


Introduction
Mapping units (MUs), which are defined as the portion of terrain where the geo-environmental conditions differ from the adjacent units across distinct boundaries (Carrara et al., 1995;Guzzetti et al., 2006Guzzetti et al., , 1999;;Hansen, 1984;van Westen et al., 1997), are essential both in predicting performance and output results, playing a very important role to determine the design of landslide susceptibility maps.The most used MUs are grid-cells units and slope-units (Reichenbach et al., 2018).Usually, grid-cells are directly derived through a Digital Elevation Model (DEM) and the resolution of the predictor variables is assumed as corresponding to that of the DEM pixels.Grid-cell partitioning is very easy and fast and, generally, a very good performance is obtained for modeling, especially when landslide initiation points are to be detected, as in the case of flow landslides (Cama et al., 2017;Lombardo et al., 2015;Rotigliano et al., 2011;Van Den Eeckhaut et al., 2009).However, gridcells units can result as too local to diagnostically express unstable conditions especially when larger slide failures are to be predicted.Moreover, pixelbased final maps frequently result hard to read and not friendly for land use planners.In fact, these maps are often mosaics of cells (generally, with metric resolution), each having a specific susceptibility, without any constraint on the spatial coherence or connectivity between the adjacent pixels in a single slope.Thereby, the territorial planning becomes very complex and complicated thus causing even an unwillingness from administrators to exploit landslide susceptibility maps as an effective tool for risk analysis and management.
Instead, the slope-units (SLU), which are defined as the slop sectors limited by drainage and water divide lines, are preferred since it is assumed the complete landslide kinematic (initiation, propagation and accumulation) occurs inside.Nevertheless, the SLU partitioning criteria may adversely affect the predictive efficacy of the susceptibility models, resulting in too large units unsuitable for picking the initiation zone and discriminating the landslides specific dynamics.Furthermore, a non-optimal SLU construction would made the proxy variable characterization into the slope-units as too smoothed (Amato et al., 2019;Ba et al., 2018;Rotigliano et al., 2011Rotigliano et al., , 2012) ) so that applying zonal statistics to predictors could be very misleading for the geo-environmental slope conditions featuring: e.g.averaging the steepness could smooth gradient jumps along the slope, which play a very important role in landslide triggering.
Thus, grid-cells units seem to be the better MU for modeling whilst a dimensionally adequate slope-unit appear the best way to obtain the final landslide susceptibility maps.Rarely, the attention of researchers has been focused to convert optimal results achieved by pixel-based modeling into successful and easy-toapply slope-based maps (Domènech et al., 2020).
In this research, several approaches to obtain a suitable landslide susceptibility map were tested.
In a test area (the Imera Settentrionale river basin, in the north sector of Sicily, Italy), a systematic rotational/translational slides inventory was obtained by using high resolution LIDAR images and a set of 12 layers expressing selected physical-environmental predictors prepared.By using MARS (Multivariate Adaptive Regression Splines), a grid-cells statistical susceptibility model was obtained, which was then integrated with 10 types of slope-units, produced through different partitioning criteria, involving half-basin extraction and landform classification subdivision.Finally, the 10 SLU-based susceptibility layers were derived by applying zonal statistic tools to dissolve the pixel-based scores.In this way, optimal slope-units partitioning scheme were explored, by comparing predictive performances and suitability of the final maps.

Study area
The Imera Settentrionale river basin extends for 343 km 2 in the northern sector of Sicily (Italy, Figure 1(a); Figure (a) of Main Map).The outcropping rocks mainly consists of carbonate, siliceous-carbonate and siliciclastic successions, which form tectonic units having imbricate geometric structures, as a consequence of the compressive phase that built up the Sicilian chain since the Oligocene (Morticelli et al., 2015).During the Late Miocene the Sicilian chain was associated with the development of wide to quite narrow syntectonic basins partly located above the growing thrust sheets (Gugliotta et al., 2013), filled with terrigenous deposits.
The area is affected by several instability phenomena (Agnesi et al., 2000;Di Maggio et al., 2017) related to the climatic conditions but also to the recurrent seismic events that occurred in the Madonie mountains area (Agnesi et al., 2005).The neotectonic activity, is also responsible for the big deep-seated gravitational slope deformation (Agnesi et al., 1997;Di Maggio et al., 2014) in the southern part of the basin.The largest part of the basin, where flysch sediments, multi-cycles sequences of alluvial clastic sediments and marine pelites outcrop, is densely marked by water erosion landforms (rills, gullies, pipes), somewhere giving rise to typical Italian 'calanchi' badlands systems (Buccolini et al., 2012;Cappadonia et al., 2011Cappadonia et al., , 2016;;Pulice et al., 2012), and landslides, which are typically of flow and slide type.
The mean annual rainfall (750 mm) and temperature values of about 16°C allow us to classify the climate of Imera Settentrionale river basin as Mediterranean type, with rainfall concentrated mainly in few of the winter semester days, while summer period is characterized by almost drought conditions.

Landslide inventory and landslide conditioning factors
Landslide recognition was carried out using high resolution (0.25 m) LIDAR (LIght Detection And Ranging) images taken in 2012 by ARTA (Assessorato Regionale al Territorio e all'Ambiente).The inventory includes 1608 rotational/translational slides that, on the whole, affect an area extended for about 26 km 2 (∼8% of the study area).The obtained landslides inventory (Figure 1 (b)) was checked with a randomly hot spots field survey which were carried out in 2019; for each of the mapped landslide polygon the highest point along the crown (Landslide Identification Point -LIP) was extracted, which several researches (e.g.Cama et al., 2015;Lombardo et al., 2014Lombardo et al., , 2016;;Rotigliano et al., 2011) highlighted as an effective diagnostic site for landslide susceptibility evaluation.
The local environmental conditions which directly influence the instability of the slopes are mainly related to morphology (slope angle, orientation slope, curvature, elevation, roughness etc.) and others predisposing characteristics (e.g.lithology, soil use, tectonic condition, landform classification).Any kind of controlling factor can be considered as a potential predictor and can be introduced into the susceptibility analysis.In this research, a set of 10 topographic predictors was derived by processing with GIS hydro-morphological tools a 8 m cell digital elevation model (DEM) obtained from a LIDAR survey in 2008 by ARTA: elevation (ELE), landform classification (LCL), steepness (STP), aspect (expressed as NORTHerness and EASTerness), plan (PLN) and profile (PRF) curvatures, topographic wetness index (TWI), terrain ruggedness index (TRI) and stream power index (SPI).
A bedrock lithology map (LITO) was prepared by grouping the different outcropping lithologies (Figure 2) in light of their expected mechanical behavior.The boundaries between the different units were verified and adapted on the basis of remote and field surveys.The Corine Land Cover 2018 (USE) was also used to classify land use of the Imera Settentrionale river basin.

MARS modeling
The Multivariate Adaptive Regression Splines (MARS; Friedman, 1991) stochastic method was applied to model landslide susceptibility.This method, which already proved to be suitable for landslide modeling (e.g.Conoscenti et al., 2015Conoscenti et al., , 2016;;Rotigliano et al., 2018Rotigliano et al., , 2019;;Vargas-Cuervo et al., 2019), allows us to optimize the linear regression obtained between the dependent and the independent variables.The MARS regression function can be written as: where y is the dependent variable (the outcome) predicted by the function f(x), α is the model intercept and β i the coefficients of the h i basis functions, N is the number of basis functions.
In this research, the MARS models were implemented through the 'earth' package (Milborrow, 2019) of Rstudio software.

Pixel-based model building and validation
Each 8 m pixel was classified as stable or unstable (positive or negative cases) depending on whether or not it hosts at least one LIP and the local values of all the proxy predictors assigned.Thereby, a R dataframe was obtained and exploited to extract a suite of onehundred balanced (positive/negative) datasets, each including the whole set of positive pixels and an equal number of randomly extracted (without replacement) subsets of negatives.Each dataset was then submitted to a balanced random partition furnishing a 75% subset for calibration and a 25% unknown subset for blind validation tests.
The model performance was analyzed in terms of a cut-off independent metric (by AUC: the Area Under the receiver operating Curve) as well as, based on an optimal cut-off (Youden, 1950), in terms of TPR and TNR (true positive and negatives rates, respectively) and related accuracy.By exploiting the availability of the one-hundred replicates, each accuracy metric was assessed also in terms of precision.
Each model was tested in predicting the validation test, first, and the whole study area, then.In fact, in evaluating the suitability of the final maps it is of great importance to verify to what extent the forced false positive generation (all the unstable status pixels are inside the calibration/validation subsets) is balanced by the performance in terms of true negative prediction.

Slope Units partitioning and scoring
By using two main criteria, 10 types of slope-unit partitioning schemes were adopted.
Firstly, two classic hydro-morphological units were obtained by means of contributing area-based procedures: in this way, by exploiting 2000 and 5000 contributing area thresholds, a 2000 half basin (2000_HB_SLU) and a 5000 half basin (5000_HB_SLU) SLU-layers were obtained.
Then, the TPI based Landform Classification (LCL; Guisan et al., 1999) SAGA tool was exploited to obtain four different LCL maps from the DEM, by changing the inner/outer radius as 100/1000, 100/2000, 500/1000 and 500/2000 meters.LCL divides and classifies the study area in geomorphological-classes considering the relative vertical position of each pixel.
By intersecting the 2000_HB_SLU and 5000_HB_SLU shapes with the four LCL maps, eight new types of slope-units were obtained (Table 1), in which each half basin unit was split into several areas having a different geomorphological classification (Figure 3).
For each of the 10 types of obtained slope-unit layers, MEAN (average) and MSTD (mean plus one standard deviation) of the pixel-based scores were zoned for producing 10 derived prediction images.A positive (unstable) status was then assigned to the slope-units if containing at least one LIP.In this way, each slope-unit was so characterized by a predicted (derived by MEAN and MSTD scoring) and an observed status, so that Receiver Operating Curve (ROC)-plots and Confusion Matrixes, through a new Youden Index cut-off, were analyzed.

Results
Figure 4(a), which shows the ROC-plot of the onehundred pixel-based models, attests for excellent (Hosmer & Lemeshow, 2000) AUC mean values (0.89), associated to very low standard deviations (0.01), and a 0.48 Youden Index cut-off.A similar high performance in terms of prediction accuracy (Figure 4(b)) was also obtained from the confusion matrix both for the balanced pixel-based model and the entire area (0.79-0.81).Coherently, both specificity and sensitivity do not change from the balanced datasets (0.83) to the entire area (0.85).On the other hand, while similar (about 0.80) Positive and Negative Predicted Values (PPV and NPV) were achieved by the balanced test, very different results were obtained for the whole area (0.001 and 0.999, respectively).
Figure 5(a) shows a graphical resume of the most important selected predictors, by using the nsubsets criterion (Conoscenti et al., 2016;Rotigliano et al., 2019).
The categorical variables (LCL, LITO, USE) were disassembled into specific classes to analyze their different impact.In this way, 23 out of 45 variables have contributed to discriminate the positive to negative cases but only seven can be considered very important.TRI, SPI, PRF, ELE, SLO and PLN are the most important DEM-derived variables.Between the categorical variables, 'lithoid units' and 'calcareous and clayey marls' lithology classes, 'non-irrigated arable land' and 'sclerophyllous vegetation' land use classes and all the LCL classes played a very important role to discriminate the stable/unstable pixels.
The binary correlation between the continuous variables (Figure 5(b)) highlights a very high negative correlation between TWI and SLO, TWI and TRI whilst a positive one between SLO and TRI.
In Table 2, the MEAN-and the MSTD-derived LCL_SLU AUCs for the 10 types of slope-units are showed.The LCL_SLU AUCs perform always slightly better by using the MSTD-derived score (>0.84) rather than the MEAN-derived score (>0.82).A similar behavior is observed also for the HB_SLU models but down-shifted toward less performing results (between 0.74 and 0.78).
The Youden Index cut-offs are lower than the pixelbased models ones, with a more marked decreasing when passing either from HB_SLU to LCL_SLU or from the MSTD-to the MEAN-derived score.
By using the MSTD-derived score, very high sensitivity values arise for all the LCL_SLU models (>0.90) together with a lower specificity (about 0.7) resulting in a less performing accuracy (about 0.7).On contrary, the 2000_HB_SLU and the 5000_HB_SLU performances, in spite of a similar (0.81 and 0.78) sensitivity, are affected by a marked decreasing in specificity (0.65 and 0.66) resulting in a just satisfactory accuracy (about 0.7).
The principal behaviors of the HB_SLU and LCL_SLU models hold also using the MEAN-derived score but all the performances result lower than MSTD-derived score: the accuracy fall down (<0.7) due to a marked specificity decreasing (<0.7, with an extreme value equal to 0.56 for 5000_HB_SLU) while the sensitivity attains good or optimal values (between 0.80 and 0.95).
In order to explore the general coherence between pixel-and SLU-based scores for the 10 different partitioning adopted criteria, the dispersion of the pixel scores zoned into the SLUs was analyzed (Figure 6).It is worth to note that a less variability arises for HB_SLUs split by LCLs with a minimum internal radius (100 m), whilst the maximum variability arises by using the HB_SLUs (especially with the 5000_HB_SLU).
In Figure 7 a comparison between the source pixelbased and the final LCL_SLU susceptibility maps is presented for a selected representative sector.In order to objectively reclassify the susceptibility maps, three cut-off values were directly obtained by analysing the ROC-plot.A first main cut-off (m cut-off) was first  identified on the ROC plot, as the probability value resulting in the maximum difference between false positive (FP) and true positive (TP) rate (Youden, 1950).The same procedure was then applied to the two semi-plots which were obtained by splitting the ROC plot using the m cut-off.In this way, the two secondary cut-off values (l and h) were identified.

Discussion and conclusions
The pixel-based model shows excellent performances both in balanced and in the entire validation area.In fact, identical excellent AUC values are obtained for the two tests and the other performance metrics reveal only a slight variation.This proves that the false positive generation in the whole area is widely balanced by the true negative successful prediction.Moreover, the one-hundred replicates demonstrate both high accuracy and precision.In light of the above, the false positive cases potentially are to be expected as future positives (in other words, future landslide initiations) rather than type I errors.
It is worth to note that the variables importance analysis shows that the LCL classes play a very important role in discriminating the landslide initiation areas.This suggests a further slope-units subdivision through the LCL classes as useful for the performant and correct identification of landslide initiation zones and, as a consequence, for developing a zonal environmental subdivision suitable for optimizing the landslide susceptibility map output.
In fact, the results confirmed that the LCL_SLU models produced homogeneous and more easy to interpret prediction images, with a predictive performance higher than HB_SLU, suffering only from a limited performance lowering with respect to the pixelbased model, almost totally dependent on type I errors (false positives); at the same time, the LCL_SLU models resulted in the highest skill in finding/predicting the unstable pixels (sensitivity).In particular, the best LCL_SLU model (HB_SLU: 5000 + LCL: 100-2000) resulted in very high performances (AUC = 0.84; sensitivity = 0.95; specificity = 0.71) and much more readable maps.Moreover, the best LCL_SLU model is   characterized by a very low pixel-score variability inside each slope-unit, attesting also for a coherent and stable susceptibility zonation.
Once the best performing susceptibility map was selected, the final step was to connect the propagation stage to the landslide initiation obtained from its prediction images.To this aim, the same morphodynamic meaning LCL_SLU was exploited after having dissolved the LCL_SLU smaller than 20,000 m 2 (the 3rd quantile of the landslide frequency distribution).A degree of fit plot was then prepared to validate the new integrated susceptibility map in predicting the source landslide polygons, which actually represent the real landslide propagation/arrest areas for the calibration inventory.This plot was then compared to that obtained from a pixel-based model, which differently was calibrated using the same landslide polygons (Figure 8).The comparison between the two results attested a very more performing behavior of the LCL_SLU integrated model.

Software
Google Earth Pro was employed to detect and map the landslides.GRASS GIS was utilized to partition the study area in SLU by using the r.watershed processing.Saga GIS 7.0.0 was employed to produce the LCL maps and to obtain the zonal statics analyzed.Quantum GIS 3.8.2Zanzibar was used to clip the SLUs with the LCL maps and to print all maps displayed in this research.The modeling and the cut-off independent and dependent matrices were carried out by using the Rstudio software.

Disclosure statement
No potential conflict of interest was reported by the authors.

Figure 4 .
Figure 4. (a) Roc plot of pixel-based model (the one-hundred replicates are plotted in orange while the averaged ROC is in red); (b) Summary of the averaged confusion matrices for the pixel-based model validations (balanced and entire area schemes).

Figure 5 .
Figure 5. (a) Variables importance for the pixel-based model: the predictors with an overall higher than 25 are plotted in cyan while those with an overall lower than 25 are in light-cyan; (b) Binary correlation between the continuous variables: the radius of circle isproportional to the intensity of correlation (between 0 and 1) while the color is representative to the direction (yellow for positive correlation, green for negative one); the X symbol corresponds to a non-significant correlation (p-value <.01).

Figure 6 .
Figure 6.Boxplot of the standard deviation pixel scores zoned into the SLUs.
All the obtained results show that by integrating the classical pixel-based modeling (Figure (b) of Main Map) on the LCL_SLU mapping units (Figure (c) of Main Map), it is possible to optimize the results of the landslide susceptibility evaluation.The same LCL_SLUs allow also to involve the propagation stage into the predictive images, obtaining easy to interpret high performing integrated susceptibility maps (Figure (d) of Main Map).

Table 1 .
Characteristics of the LCL_SLU partitioning scheme adopted.

Table 2 .
Summary of the validation metrics for the 10 susceptibility models.