Characterization of Susceptible Landslide Zones by an Accumulated Index

Abstract

In order to characterize the landslide susceptibility in the central zone of Guerrero State in Mexico, a spatial model has been designed and implemented, which automatically generates cartography. Conditioning factors as geomorphological, geological, and anthropic variables were considered, and as a detonating factor, the effect of the accumulated rain. The use of an inventory map of landslides that occurred in the past (IL) was also necessary, which was produced by an unsupervised detection method. Before the design of the model, an analysis of the contribution of each factor, related to the landslide inventory map, was performed by the Jackknife test. The designed model consists of a susceptibility index (SI) calculated pixel by pixel by the accumulation of the individual contribution of each factor, and the final index allows the susceptibility cartography to slide in the study area. The evaluation of the obtained map was performed by applying an analysis of the frequency ratio (FR) graphic, and an analysis of the receiver operating characteristic (ROC) curve was developed. Studies like this can help different safeguarding institutions, locating the areas where there is a greater vulnerability according to the considered factors, and integrating disaster attention management or prevention plans.

1. Introduction

The landslides are phenomena of great importance due to their spatial dimension effects on socio-ecological systems. Typically, the landslides are present in areas where hydrometeorological phenomena also occur [1, 2, 3]. In recent years, an increase in landslide risk has been recorded due to the increase of population settlements in vulnerable areas such as coastal and mountain regions.

This problematic and global situation has aroused the interest of governments and academics, and significant efforts have been made to characterize and identify the areas that have the potential to suffer landslides. Under a territorial approach, studies and susceptibility models have been developed, where factors related to instability processes and their spatial distribution are analyzed and linked to the landslide inventories.

The remote sensing techniques and the management of geospatial information through geographic information systems (GIS) currently play a leading role in numerous studies in which statistical or semi-statistical methods are applied to the modeling and generation of landslide susceptibility cartography [4].

In the last four decades, different probabilistic methods have been developed to predict, model, assess, and produce cartography about the risk related to landslides. Van Westen et al. [5] propose three scale levels in risk maps:

1. The qualitative-heuristic approach, with empirical recognition and small-scale maps (1: 100,000–1: 250,000)

2. The statistical approach, to determine the causal factors in the quantitative susceptibility mapping (scale 1: 25,000–1: 50,000)

3. The deterministic approach, for detailed large-scale studies (1: 2000–1: 10,000)

Guzzetti et al. [6] describe that the most important proposed methods can be grouped mainly in the following categories:

1. Geomorphological risks mapping

2. Analysis of landslide inventories

3. Heuristic or indexing methods

4. Functional and statistical models

5. Geotechnical models or based on the physical structure

The literature about the evaluation of susceptibility to landslide events is vast; Carrara [7], Guzzetti [8], Chung et al. [9], Baeza and Corominas [10], Hervás et al. [11], Ayalew and Yamagishi [12], Lee [13], Mancini et al. [14], Kayastha et al. [15], and Dou et al. [16], among others, use the statistical methods of multivariate analysis (discriminant analysis, processes of analytical hierarchies or logistic regression). According to Feizizadeh and Blaschke [17], the use of these methods has increased in landslide susceptibility mapping; however, the results are sensitive to the quality of the dependent variable, for instance, the inventory of previous landslide events [4, 18, 19].

Other authors have analyzed damaged scenarios with a geostatistical approach applying kernel density maps to evaluate geological, geomorphological, and risk data [20]. Lazzari and Danese proposed a multi-temporal kernel density estimation approach for local landslide susceptibility evaluation and forecasting, based on spatial statistics techniques and in particular on kernel density estimation, considering the interaction between landslides that are located close to each other, as a second-order effect in landslide spatial distribution [21].

Due to the complexity of landslides, many authors have focused on nondeterministic methods to assess the susceptibility. Lee et al. [22], Gómez and Kavzoglu [23], Park et al. [24], and Conforti et al. [25] have applied data mining methods using artificial neural networks, while Ercanoglu and Gokceoglu [26] and Pistocchi et al. [27] have applied the fuzzy logic methods.

Ward et al. [28], Terlien et al. [29], Fall et al. [30], and Jia et al. [31] have adopted a mechanical approach, using deterministic and/or numerical methods, which have been successful for the analysis of the stability of a slope [32].

In this study, we apply a heuristic method, which considers the relationship between the location of landslide events previously occurred and the geomorphological and geological environment [33, 34, 35, 36]. In the heuristic methods, the experts’ criteria are essential, due and based on their expertise, and they determine the factors that influence the instability [37]; here each factor is weighted according to its importance/influence on the landslide triggering. A final map is obtained by superimposing the thematic layers related to each factor. Examples of the application of this type of indirect methods are described below.

Irigaray [38, 39] proposes the critical slope method, with which the susceptibility of lithology in slope intervals is estimated. The critical slope per lithological unit is obtained from the relation of the increase of the percentage occupied by the surface of rupture of the class of immediate lower slope and the increase of slope of both classes.

Van Westen [40, 41] developed a procedure that calculates the landslide density, assigning a weight to the variables (factors), through the difference with the total density and finally adding the weighted densities. An important feature is that this method contemplates different types of landslides and activity stages, considering the slid zones.

Sarkar et al. [42] performed a similar procedure, but without assigning weights, they calculate the landslide density for each class and factor, obtaining a susceptibility index through the weighted accumulation of the individual indices corresponding to each factor included and thus obtaining the final susceptibility map.

Barredo et al. [43] apply two knowledge-based methods. The direct method performs a very detailed geomorphological cartography using specifically coded polygons that were evaluated one by one by experts to assess the type and risk degree. The indirect method is based on an indexing technique, where the conditioning factors are combined using multicriteria evaluation techniques.

Hervás and Barredo [44] and Hervás et al. [11] propose an indirect method of indexing to evaluate the landslide susceptibility, through the pondered linear sum of the weights of the considered factors and classes, for the assignation of the final weight corresponding to each conditioning factor. In their study, they applied the analytical hierarchy method, creating a matrix with the relative value judgments between pairs of conditioning factors of instability. The resulting susceptibility index was presented in a five-category thematic map.

Ayalew et al. [34] applied the linear weight combination method to determine classes of control parameters (factors): lithology, slope, terrain orientation, elevation, profile, and flat curvature, considering the order of importance in the process of landslides and weighing the impact of one parameter against another.

Wati et al. [35] applied a heuristic approach with the weighted scoring method, where the judgment of an expert determines the weighting and scoring, which represents the effect of each factor in the process of a landslide. The higher the weight and the score, the higher will be the influence of a particular factor to trigger a landslide. The weight designation to the corresponding factors was aided by the Integrated Land and Water Information System (ILWIS) software, and they used a rank of values, 1 to 4. The results were presented in a thematic susceptibility map organized in five classes whose ranges were determined by the Natural Breaks method [45] based on the accumulated weight score.

The work presented here is only one phase of a global susceptibility analysis project in the Guerrero State, which includes landslides, floods, earthquake damage, and environmental effects due to industrial activity, among other phenomena that occur in the study area.

Under this background, the present work consists of the analysis, design, and implementation of a spatial model to characterize the landslide susceptible zones in the central area of Guerrero State in Mexico.

2. Study area

The study area covers a space of 3300 km², in the central zone of Guerrero State in Mexico, and consists of a mountainous region, with elevations ranging from 280 m to 3540 m above mean sea level and slopes over 40 (Figure 1).

Figure 1.

Study area, the south-central region of the state of Guerrero, Mexico. Source: Adaptation of ASTER image from Dec 13 2013.

According to the Mexican Institute of Statistics and Geography (INEGI), for the year 2013, temperature variations were recorded from 14.3°C (December) to 28.3°C (May) [46].

The records of precipitation of the National Meteorological Service (SMN-CONAGUA) for the year 2013 indicate that the ranges vary from 800 mm as minimum average to 2100 mm as maximum recorded from June to September [47]. The area was 74.8% covered by forest (coniferous, mesophilic, and mixed), 14.1% by deciduous forest, 7.8% by agricultural land, 3.2% by induced vegetation, and 0.1% by human settlements and urban areas [48]. According to the 2010 Population and Housing Census, there are 187 cities and towns, in which there are 15,230 homes inhabited by 59,098 people [49].

Geologically, the area is located physiographically in the Sierra Madre del Sur [50] and has a variety of metamorphic rock composition consisting of schists and gneisses of biotite and quartzite, outcrops of deposited limestone, metavolcanic rocks with sedimentary influence, siltstones, sandstones, conglomerates, and carbonate rocks. Rhyolitic rocks are also found as a result of Oligocene-Miocene volcanism. The youngest rocks correspond to alluvial deposits present in the margins and riverbeds [51, 52]. Figure 2 shows the geological map of the study area, integrated from two geological charts and ten geological-mining charts, referenced to 2012.

Figure 2.

Geological map of the study area [53].

The study area is interesting because of its geographic location, topographic, and geological conditions and also due to the presence of extraordinary hydrometeorological phenomena in recent years, which have triggered massive landslides that have severely affected the population and infrastructures of their communities. The events of September 2013 were particularly notable during tropical depression No. 13 in the Pacific Ocean and subsequent simultaneous hurricanes: Manuel in the Pacific and Ingrid in the Gulf of Mexico caused significant floods and landslides on the coast of the Guerrero State [54] (Figure 3).

Figure 3.

A massive landslide occurred on La Pintada community. September 2013, in the study area. Source: Personal adaptation from https://www.jornada.com.mx/2013/09/24/ciencias/a03a1cie.

In the study area, there is an important concentration of inhabitants. From the 187 localities distributed in the area, 182 are considered as rural (with less than 2500 inhabitants), and their environmental conditions, the materials of their homes, and characteristics are similar to La Pintada.

3. Dataset and methodology

An overview of the process to characterize, through an accumulated index, the areas susceptible to landslides is shown in Figure 4.

Figure 4.

Methodology overview.

As can be seen, the proposal is made up of five stages:

• Dataset. Consists of the acquisition of primary data and the generation of maps derived from it.

• Pre-processing. In this stage, additional processes are applied to the derived information, such as the categorization of maps with continuous data, the generation of kernel density maps, or the combination of maps with related information.

• Jackknife test. In this stage, the factors considered are analyzed to determine the contribution of each of them to the global susceptibility model.

• Accumulated index of landslide susceptibility. Here, the proposed model is implemented to build the landslide susceptibility map.

• Accuracy assessment. At this stage, the map generated is evaluated.

Each of the stages is explained in detail throughout the work.

3.1 Landslide inventory

An image dated on August 12, 2014, from Google Earth online platform was used, just the closest one available after the meteorological events of September 2013 previously described. Photointerpretation techniques were applied to identify and digitize a sample of 617 polygons of which 602 (17,385 pixels) represent real landslide areas and 15 (17,632 pixels) represent non-landslide zones to integrate an inventory validated through field visits in the study area (Figure 5A).

Figure 5.

Landslide inventory. (A) Landslide polygons digitized in Google earth validated. (B) landslide inventory map resulting from unsupervised detection methods [53].

Due to the conditions of access and insecurity prevailing in the study area, not all the digitized polygons were field verified. The field verification was in 321 polygons (53.3%) in secure zones near to the largest cities and close to the most important roads. The field verification indicated that 320 were translational landslides, and only 1 was a debris avalanche.

The size of the slid polygons in the integration of the landslide inventory was 450 m², which is the area that corresponds to two pixels, due to the spatial resolution of the ASTER images used, which is 15 m.

A landslide inventory map (Figure 5B) was used to assess the resulting susceptibility landslide model. This inventory was automatically produced through unsupervised detection methods applied to ASTER-derived products [53].

The described landslide inventory produced by unsupervised detection methods consists of 216,174 pixels that represent the landslide zones caused by the extraordinary rains of the simultaneous hurricanes Ingrid and Manuel in September 2013. These pixels were divided into three equal parts, of which 2/3 (144,116 pixels) were used to define the degree of contribution of the factor terrain orientation and flat curvature to the susceptibility model.

The other part (72,058 pixels) was complemented by a similar number of randomly selected pixels of the inventory map that were identified as non-landslides and reserved for the stages of evaluation of the influence of factors on the occurrence of landslides and to assess the accuracy of the landslide susceptibility map obtained in this work.

3.2 Conditioning and detonating factors

In this work, the quantity and type of the analyzed factors were the results of a selection process supported on the expertise or researchers based on specific characteristics of each factor and the conditions of the study area. Thus, it is important to emphasize that finally, the explanatory factors are considered, depending on the characteristics of the study area and the availability of corresponding information. The factors considered, as well as the processes applied for their integration, are described as follows.

Based on the data obtained from the topographic maps, courtesy of National Institute of Statistics and Geography of Mexico (INEGI), a 15 m grid cell digital elevation model (DEM) was generated; then, maps of the angular slope, terrain orientation, drainage network, and flat curvature were produced from DEM to be used as factors in the susceptibility modeling.

From the topographic maps, the road infrastructure data was obtained. The lithology map was obtained from the geological-mining charts, while the map of the structural elements was obtained from the geological chart courtesy of INEGI and the Mexican geological service (SGM).

From the soil chart of Guerrero State (INEGI), the textural class map of the study area was obtained. Also, a map of accumulated precipitation was generated from the data bank of the daily precipitation from January to September 2013, obtained from the meteorological stations located in or near to the study area (Figure 6). With the data of accumulated precipitation, textural class, and terrain slopes, a new combined product was generated that was called the potential rain effect map.

Figure 6.

Meteorological stations considered in the study. Source: Personal adaptation from ESRI and SMN-CONAGUA.

According to the kind of data of each factor, in some cases, additional processing was performed to be used in the modeling of susceptibility. For the terrain slopes (S), in concordance with Alcantara-Ayala [55], the slope values smaller than 5° were considered as flat zones, and then without effect on the susceptibility index, while for the upper slope values, a limit value through the statistical parameters of all slope data was determined, calculated by the mean plus twice the standard deviation (μ + 2σ). The resulting value was 45.3°, and from this value onwards, it was considered that the effect on the susceptibility to slip is the same, and thus all the upper-value slopes were equalized considering them with the same weight (Figure 7A) [56].

Figure 7.

Explanatory factors of landslide susceptibility in the central zone of Guerrero state. (A) Terrain slope. (B) Potential erosion index. (C) Drainage network density. (D) Road density. (E) Structural element density. (F) Terrain orientation. (G) Flat curvature. (H) Potential rain effect.

About lithology, a susceptibility degree was defined for each lithological unit. The map was categorized according to Aramburu Maqua and Escribano Bombín [57], who refers to the degree of cohesion, potential erosion, and mechanical behavior of each type of rock.

In this sense, the potential erosion index (PEI) was used, since it refers to erosion by surface runoff due to hydrological factors, which are considered in this study as the primary trigger of landslides. Based on the analysis of the lithological content of the rocks present in the area, five classes were defined that indicate the susceptibility degree, related to the potential erosion of the rock type. According to the categorization based on the PEI [57] and the identified rock types, four susceptibility classes were determined: (1) very low, (2) low, (3) medium, and (5) very high (Figure 7B). It is important to note that the classification corresponds to the cited works consulted, and according to the analysis performed, class (4) high does not exist in the study area.

Regarding the drainage network maps, road infrastructure, and geological structures, it was considered to use density maps, since if a point in the space next to an “element” (runoff, road, or geological structure) has a certain degree of susceptibility, then this should be greater if the same point is close to two or more elements; therefore, a density map would represent the most susceptible areas to slide according to its position with respect to one or more of the elements considered.

From each map (drainage network, road infrastructure, and geological structures), the density maps were generated by the kernel density estimator (KDE), which calculates a magnitude of influence of the “elements” per unit area according to the neighborhood of the pixels [20, 58]. The KDE can be applied on specific entities or lines. KDE implies the use of parameters, among which the search radius (bandwidth) is one of the most important and represents the value of the radius within which the density will be estimated.

The default search radius (SR) is computed according to the spatial distribution of input dataset to be analyzed using a spatial variant of Silverman’s rule of thumb that is robust to spatial outliers (i.e., points that are far away from the rest of the points), through the equation

$$\mathit{SR}=0.9\ast \mathit{min}\left(\mathit{SD}\sqrt{\frac{1}{ln\left(2\right)}\ast \mathit{Dm}}\right)\ast n^{-2}$$ E1

Conceptually, an adjustment of a uniform curved surface is made on each “element.” The value of KDE is higher on the line and decreases as the pixels are farther from the feature, even reaching the zero value. To calculate the density of each pixel, the values of all the superposed kernel surfaces are added [59] (Figure 7C–E).

The terrain orientation (TO) was classified according to Rawat et al. [60], excluding the value of −1 since it represents flat areas and does not contribute any degree of susceptibility to landslides in the model. The rest of the orientation values were organized in eight classes, to which a cross-tabulation was applied with 2/3 of the slides taken from the ground truth (landslide inventory map) to determine the frequency of each class of orientation of the terrain in which the pixels are slid.

From the crossing of information of the class of terrain orientation and the pixels of the landslide inventory map, the map of orientations was reclassified according to the frequency of the slid pixels in each class, being the order of frequency and the degree of susceptibility assigned (Table 1, Figure 7F).

N.	Class	Azimuth terrain orientation (°)	Frequency (%)	Reclassification
1	North	337.5–22.5	2.48	1
2	Northeast	22.5–67.5	5.42	3
3	East	67.5–112.5	13.95	5
4	Southeast	112.5–157.5	19.74	7
5	South	157.5–202.5	25.36	8
6	Southwest	202.5–247.5	17.68	6
7	West	247.5–292.5	11.31	4
8	Northwest	292.5–337.5	3.60	2

Table 1.

Reclassification of terrain orientation according to slid pixel frequency.

Source: [56].

The flat curvature (FC) is an important topographic feature commonly used in landslide susceptibility studies [61]. The slopes can be subdivided into regions of flat curvature, lateral concave, lateral convex, and linear; however, statistical analysis of the flat curvature and landslide datasets indicates that slopes with linear flat curvature have the highest probability of landslides, and this probability decreases as the slopes become more concave or convex [62].

The flat curvature map contains values from −25,798 to 30,944, which were adjusted according to the empirical information of the landslide inventory map, using a similar methodology to that applied in the terrain orientation; first, the curvature values were analyzed, and then the data were cross analyzed with the landslide inventory information, using two-third of the slid pixels considered as ground truth. Thus, the sampled slid pixels show curvature values from −12,427 to 5486; to these range values, a classification was applied in three ranges by the Natural Breaks method [45]; and for each class, the frequency of slid pixels was obtained (Table 2).

Class	Minimum	Maximum	Pixels %
1	−12.427	−1.328	10.8
2	−1.328	0.358	72.1
3	0.358	5.486	17.0

Table 2.

Flat curvature and relative frequency from the sampled slide pixels.

Source: [56].

Most of the sampled pixels (72.1%) are included in the range (−1.328 to 0.358) that corresponds to class 2; therefore this class was assigned with the highest degree of susceptibility, and according to the frequency values, to class 1 the medium susceptibility was assigned, and to class 3 the low susceptibility was assigned. For curvature values outside of the analyzed ranges [values lower than −12,427 (convex areas) and higher than 5486 (concave zones)], a null susceptibility was assigned for effects of the development of the susceptibility model; this is because they were not linked to slid pixels, according to the empirical analysis (Figure 7G).

The potential rain effect (PRE) map (Figure 7H) did not require any empirical analysis.

3.2.1 Evaluation of the influence of factors on the occurrence of landslides: the jackknife test

Once the information of the considered factors in the study was integrated, the Jackknife test was applied, which allows identifying the importance degree of each factor by executing the models in a separate way, with different sets of factors and excluding the factor analyzed. The Jackknife test provides alternative estimators of landslide patterns as a function of the factors, allowing so the identification of the significance of each of them, by their contribution degree and by their importance regarding the global results of the model [63, 64, 65, 66].

In this study, the Jackknife test graph was produced as a useful tool to assess the individual contribution of each factor used and also to assess the contribution of each range or category of each factor to the susceptibility landslide.

3.3 Accumulated index of landslide susceptibility

The applied model consists of the integration of an accumulated global index of landslide susceptibility, calculated from the individual contributions of each factor in the following described stages.

3.3.1 Factors standardization

Each factor reports continuous or discrete values. Therefore, a standardization process was applied; in this way, individual contribution values were assigned to the susceptibility model, calculated between 0 (null contribution) and 1 (maximum contribution); and thus, evaluate the corresponding contribution in the global accumulated susceptibility index. So, all the factors will have the same weight in the model.

The standardization process was made according to the data of each factor. Thus, for instance, for the PEI factor which is represented in a thematic map of five categories of contribution to susceptibility ((1) very low, (2) low, (3) medium, (4) high, and (5) very high), the standardized value was determined by finding the relationship of each pixel in connection with the value corresponding to the highest category, resulting in the PEI factor with values between 0 and 1.

For the standardization of some other factors, such as the flat curvature or the terrain orientation, the landslide inventory map was used to empirically obtain range values of curvature or orientation in slid areas (ground truth). A schematic example of the terrain orientation standardization process is shown in (Figure 8).

Figure 8.

Terrain orientation standardization process.

The standardization of the factors registered in maps of continuous values was determined by finding the relationship between the value of each pixel and the highest pixel value recorded.

All the models and processes described in this paper were developed using the Dinamica EGO software [67, 68, 69].

3.3.2 Landslide susceptibility map through the susceptibility index (SI)

The final value of the SI is integrated by the accumulation of the individual contribution of each factor considered in the study [70], calculated by the equation:

$$\mathit{SI}=S_S+{\mathit{\text{TO}}}_S+{\mathit{DD}}_S+{\mathit{FC}}_S+{\mathit{PEI}}_S+{\mathit{SED}}_S+{\mathit{RD}}_S+{\mathit{PRE}}_S$$ E2

where

$\mathit{SI}=$ Cumulated landslide susceptibility index.

$S_S=$ Standardized terrain slope.

$\mathit{\text{TO}}=$ Standardized terrain orientation.

${\mathit{DD}}_s=$ Standardized drainage network density.

${\mathit{FC}}_S=$ Standardized flat curvature.

${\mathit{PEI}}_S=$ Standardized potential erosion index.

${\mathit{SED}}_S=$ Standardized geological structural elements density.

${\mathit{RD}}_S=$ Standardized roads density.

${\mathit{PRE}}_S=$ Standardized potential rain effect.

Under the premise that each factor provides a degree of susceptibility with standardized values ranging from 0 (null) to 1 (high), the final map obtained represents the pixel by pixel accumulated index of landslide susceptibility of the analyzed area, according to the considered factors. It is important to note that from the obtained results of the Jackknife test, the accumulated SI model can be modified by excluding some factors due to its degree of contribution.

3.4 Accurate assessment of the accumulated index of landslide susceptibility

The accuracy in the susceptibility mapping refers to the degree of agreement between the prediction of areas susceptible to landslides obtained by the applied methods and the location of actual landslides, according to the information considered as reference data or ground truth in the landslides inventory.

The implementation of the accumulated SI model produces a susceptibility map with continuous values, in which, according to their ranks, applying natural breaks were categorized into five classes: (1) very low, (2) low, (3) medium, (4) high, and (5) very high. On the other hand, there is a landslide inventory map organized into two categories: (1) non-landslides and (2) landslide. So, to determine the accumulated SI model accuracy, two methods were used:

1. The radio frequency graph (RFG). The frequency ratio was graphed, taking the slid pixels from the inventory sample reserved for the validation stage and superimposing them on the landslide susceptibility map resulting from the applied SI method. We identify and quantify the coincidence frequency of slide pixels in both: the map produced by the model and landslide inventory map, for each class considered in the model (very low, low, medium, high, and very high) [71]. In an ideal landslide susceptibility map, the frequency value should be increased from a low susceptibility to a very high susceptibility [72, 73].

2. The analysis of the receiver operating characteristic (ROC) curve. It is a standard method to evaluate the accuracy of a diagnostic test [74]. The ROC curve is a comparative graphical representation of ordered pairs between the false-negative and false-positive rates, resulting in a diagnosis for classes of values. Conventionally, the graph shows the rate of the pixels diagnosed as false-positive (false-positive rate (FPR)) on the X-axis (Eq. (3)) and the rate of the pixels diagnosed as true positive (true-positive rate (TPR)) on the Y-axis (Eq. (4)):

$$X=\mathit{FRP}=1-\left[\frac{\mathit{TN}}{\mathit{TN}+\mathit{FP}}\right]$$ E3

$$Y=\mathit{TRP}=1-\left[\frac{\mathit{TP}}{\mathit{TP}+\mathit{FN}}\right]$$ E4

where

$\mathit{TN}$ represents the true-negative pixels.

$\mathit{FP}$ represents false-positive pixels.

$\mathit{TP}$ represents the true-positive pixels.

$\mathit{FN}$ represents false-negative pixels.

The area under the ROC curve (AUC) characterizes the quality of the predicted system, describing its ability to anticipate the occurrence or nonoccurrence of predefined events. The value of the AUC varies from 0.5 to 1.0. If a model does not predict the occurrence of landslides better than randomly, then the AUC would be equal to a value of 0.5, while a ROC curve with an AUC value of 1 represents the perfect prediction. The quantitative-qualitative relationship between the AUC and the prediction accuracy can be classified as follows: values from 0.9 to 1, excellent; from 0.8 to 0.9, very good; from 0.7 to 0.8, good; from 0.6 to 0.7, medium; and finally from 0.5 to 0.6, poor [75].

In the application of accuracy assessment, one-third of the total of pixels identified as a landslide in the inventory map (72,058 pixels) were used, which were omitted from the previous stages of training and preparation of the factors. This information was supplemented with an equal number of pixels randomly selected from the landslide inventory map that were identified as non-landslides.

4. Results and discussion

4.1 Analysis of the jackknife test

As we explained before, the Jackknife test allows identifying the degree of importance of each considered factor, executing the model in a separate way with different sets of factors and excluding each factor in turn. Figure 9 shows the results of the Jackknife test for the factors used in the present study.

Figure 9.

Jackknife test result graph.

The graph shows in blue the individual contribution of each factor in the implementation of the model; in cyan color, the gain of the model is shown, excluding the factor analyzed; and finally, the bar in red represents the gain of the model including all the factors.

According to the graph, we can see that the factor with the highest individual contribution corresponds to the potential erosion index, followed by terrain orientation and terrain slope factors. Also, we can see that the factor with the lowest contribution to the model corresponds to the flat curvature and can also be verified that the gain of the model is not affected if the flat curvature is omitted in the execution of the model.

One of the most important contributions of the Jackknife test is the possibility of identifying the degree of contribution and importance of each of the analyzed factors. Thus, we can determine which of them are essential in the reproduction of the pattern represented by the random sample of the landslides observed in the ground truth and so determine its degree of importance concerning the results reported by the global model. Thus, it is also possible to identify those factors that do not substantively contribute individual gains to the global model and observe the gain that is achieved when they are excluded. It is also possible to consider eliminating them in the model implementation, knowing previously the result that would be obtained when doing so. So, it is theoretically possible that the accuracy of the results achieved by the global model remains the same, and it is even likely that there will be improvements in the overall result by excluding them from the model in a new performance.

According to the results of the Jackknife test (Figure 9), it may be possible to propose the elimination of the flat curvature in the application of the model, since it is the factor that records the lowest individual contribution to the global result and at the same time is the factor that if omitted in the execution would not affect the overall result of the model.

4.2 Accumulated SI model

The implementation of the cumulative global model resulted in the landslide susceptibility map shown in Figure 10.

Figure 10.

Map of the cumulative susceptibility index model.

The resulting map is represented in continuous values, which according to their ranks was categorized into five classes; thus, the areas identified as (1) very low and (2) low susceptibility can be seen in blue tones; in shades of yellow are the areas with category (3) medium and in orange and red the areas evaluated by the model with categories (4) high and (5) very high susceptibility.

A concentration of areas with the highest susceptibility levels can be seen in the north and center of the study area, while the zones with low and very low susceptibilities are concentrated towards the south.

Analyzing the behavior of the factors in the zones identified with high susceptibility, we can note that in these areas the slopes have a terrain orientation to the south, southeast, and southwest (Figure 7F), which favor the occurrence of landslides in the study area. On the other hand, in those same areas identified as having high susceptibility, there are slopes greater than 30° (Figure 7A), which further favors, in addition to the orientation, the presence of landslides.

It should be noted that the SI method considers the logical, empirical, and even intuitive effect designated to each factor, that is, higher terrain slopes correspond to greater susceptibility; greater proximity to faults also corresponds to greater susceptibility. In other words, the weight is not assigned to the factors according to the presence of landslides but mainly is assigned as a function of the characteristics and factor attributes, although, as previously explained, in some specific cases the contribution of susceptibility index was determined empirically according to the presence of landslides.

Figure 11 shows a zoom to the northeast zone of the study area, which is characterized by having recorded a significant presence of landslides occurred in September 2013.

Figure 11.

Massive landslide zone occurred in 2013.

In the figure, an acceptable degree of accuracy is observed in the susceptibility map generated by the SI model, since in general, the occurrence of landslides coincides with the areas categorized with high or very high susceptibility.

It can be observed that most of the landslides occurred to the north of the stream are over an area identified by the model in the very high susceptibility class, and the few landslides that did not meet this condition occurred in another area identified as high susceptibility. On the other hand, also it can be seen that most of the landslides located to the south of the stream occurred over areas identified as high susceptibility, but it is also possible to see landslides recorded in the ground truth that occurred in areas identified in the medium susceptibility category.

4.3 Accuracy assessment of the accumulated SI model

4.3.1 Analysis of the radio frequency graph (RFG)

Following the methodology for the accuracy assessment, the sampled pixels categorized as landslides from the map inventory reserved for validation were superimposed on the resulting susceptibility map from the implemented model, and then the frequency of the occurrence of landslides was obtained on each map category, building the graph of the frequency radius (Figure 12).

Figure 12.

Frequency ratio graph.

It can be observed that the pattern of the resulting graph shows that the frequency of the slid pixels of the ground truth is lower in the very low susceptibility category, and this value gradually increases over the following categories to finally reach the highest frequency value in the very high category; in concordance with Pradhan and Lee [72], Pourghasemi et al. [73] and Pradhan and Lee [76] in their work conclude that in an ideal map of landslides susceptibility, the value of the slid frequency should gradually increase from a low susceptibility zone to very high one.

If we group the extreme categories, the very low and low classes in one hand, and the high and very high classes on the other, we confirm that the susceptibility map shows an accumulated frequency of 15.4% of pixels slid in the very low and low susceptibility categories, and on the other end (high and very high susceptibility categories), we can see that the susceptibility map shows an accumulated frequency of 68.1% of pixels slid. This fact confirms the positive trend of the results obtained in congruence with the previously cited works.

4.3.2 Analysis of the receiver operating characteristic (ROC) curve

From the cross-validation between the categories included in the SI model, and the random sample of the ground truth taken from the landslide inventory map, including for this case both landslides and non-landslides, the graph corresponding to the ROC curve was constructed, and the area under the curve (AUC) was calculated (Figure 13).

Figure 13.

ROC curve and AUC from the accumulated SI model.

As previously mentioned, the AUC describes the ability of a model to predict the occurrence or nonoccurrence of predefined events. According to Chen et al. [77], the higher AUC represents the higher performance of the model analyzed.

The AUC value reported by the model is 0.79, indicating that its predictive capacity is good, under the classification values consulted in the work of Yesilnacar [75], Devkota et al. [78], and Kim et al. [79].

According to the results, the SI accumulated model, in the way that has been implemented in the present work, has proven to be a valid method to characterize the areas susceptible to the occurrence of landslides in the study area. However, one of the main problems consists in assigning the correct weight value to each of the explanatory factors considered. According to Barredo et al. [43] and Clerici et al. [80], the heuristic methods use selective criteria that require expert knowledge to be applied properly, which implies a substantial degree of subjectivity insofar as each factor is assigned with a certain degree of a priori importance.

It is important to highlight the design of the model to be implemented and the explanatory factors to be considered depending on the specific characteristics of the study area and the availability of the corresponding data. In other words, the results shown here may be different if the described methodology and the same amount and type of explanatory factors are applied in a study area with characteristics different from those considered in the present study. However, the methodological proposal can be taken as a starting point and as a guide for the specific design of a model applied to a different study area and adjusted according to the available information and the topographic, geological, hydrological, or environmental conditions prevailing in that area.

5. Conclusions

According to the Jackknife test, lithology is the factor that records the most significant individual contribution to the global susceptibility, followed by the orientation and slope of the terrain. On the other hand, the factor with the least contribution to the model is the flat curvature.

According to the results obtained by the SI method, the heuristic approach proved to be valid and straightforward, in agreement with similar works [43], to assess the susceptibility of landslides on a medium scale such as the central zone of the State of Guerrero.

Cartographic information about landslide susceptibility is necessary for areas prone to landslides since it can be used as a support to the establishment of early warning systems for localities and residents in risk areas. It can also be a support to the institutions and organizations responsible for safeguarding in the management of events of this type, locating the areas with a higher vulnerability according to the factors considered and integrating care or disaster prevention plans considering areas with higher priority than others. It is also important to mention that this work can be complemented with the definition of evacuation routes, identification of care centers, and identification and capacity of shelters, among other possible actions, all this putting as a fundamental focus, the knowledge of the territory through the application of geographic information technologies.