Next Article in Journal
Improving the Completion of Weather Radar Missing Data with Deep Learning
Previous Article in Journal
The All-Solid-State Narrowband Lidar Developed by Optical Parametric Oscillator/Amplifier (OPO/OPA) Technology for Simultaneous Detection of the Ca and Ca+ Layers
Previous Article in Special Issue
Coexistence of a Marginal Mountain Community with Large-Scale and Low Kinematic Landslide: The Intensive Monitoring Approach
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Evaluation of Susceptibility by Mass Movements through Stochastic and Statistical Methods for a Region of Bucaramanga, Colombia

by
Joaquín Andrés Valencia Ortiz
1,*,
Antonio Miguel Martínez-Graña
1 and
Lenny Mejía Méndez
2
1
Department of Geology, Faculty of Sciences, University of Salamanca, Plaza de los Caídos s/n, 37008 Salamanca, Spain
2
Geophysical Institute (GPI), Karlsruhe Institute of Technology (KIT), Hertzstr. 16, D-76187 Karlsruhe, Germany
*
Author to whom correspondence should be addressed.
Remote Sens. 2023, 15(18), 4567; https://doi.org/10.3390/rs15184567
Submission received: 17 June 2023 / Revised: 8 September 2023 / Accepted: 14 September 2023 / Published: 16 September 2023

Abstract

:
Mass movements are one of the hydrometeorological phenomena with the most negative impacts on the study area, and their evaluation through the calculation of susceptibility provides a tool of vital importance within territorial planning and disaster risk management on natural and anthropic environments. Their evaluation took algorithms designed within stochastic and statistical methods, such as the artificial neural network, the bivariate statistical method, and the logistic regression method, which integrate inherent variables (geoenvironmental characterization) against events or dependent variables. This correlation simulates regions with a probability of occurrence of mass movement under training or weight assignment. Its construction for this study took, as a basis, 50% of the events (test) and 50% of the events (validation) randomly and with equivalent area distribution against the inherent variables. As a result, it was observed that the bivariate method presented a good performance in spatial prediction. This model presents values of AUC = 82.2% (test) and AUC = 76.9% (validation), grouping a total of 591 events of the 856 events in the high category (69%). In turn, from a second evaluation carried out by this method to each hydrographic basin, a condition was established in the area (50 km2) for coherent results at a level of analysis 1:25,000, based on the idea that the variables do not present changes greater than 20% in their attributes, added to a knowledge of the area evaluated.

Graphical Abstract

1. Introduction

The analysis of mass movements (MM), among other hydrometeorological phenomena in recent decades, has been undergoing a great boom due to its impact on anthropogenic elements especially. Although these natural phenomena are not recent and are closely linked to Earth dynamics, in recent decades, their studies have increased due to the frequency that has been occurring and its consequences. A study conducted by the World Health Organization between 1998 and 2017 determined that mass movements affected 4.8 million people, as well as a total of 18,000 human losses [1,2]. Relative to other natural hazards, mass movements account for 4.9% of all natural disasters and 1.3% of all deaths from natural disasters between 1990 and 2015 in the world [3]. These conditions are largely presented by unplanned territorial development and incorrect land use. In turn, this problem is due to an implicit complexity in the understanding and representation of the geoenvironmental and sociocultural conditions that determine the degree of stability of the terrain in the occurrence of a mass movement (Colombian Geological Service—SGC, [4]).
Within the context of the study area, mass movements have generated great economic losses in infrastructure since they have partially or totally affected roads, bridges, and homes, as well as caused human losses. According to the database of the Disaster Information Management System—DesInventar, [5], these damages and losses have occurred especially between the years of 1939, 1972, 1979, 2010, 2011, 2013, 2015, and 2020, with rainfall being the most important reported triggering factor for the generation of mass movements, especially in the years of 2010, 2011, and 2020. In turn, according to a report made by the Institute of Geoscientific Research and Information (Ingeominas) in 2002, the region where the study area is located ranks fourth nationally with the largest number of records of mass movements [6], information which is complemented by a study carried out by the SGC at a scale of 1:100,000, where it is observed that the study area presents a high hazard to mass movements [7].
To understand this problem a little more thoroughly, it is necessary to study the conditions where mass movements occur that are considered as a displacement of a mass of rock, debris, or land downhill [8]. These conditions or inherent elements of the terrain are part of a geoenvironmental characterization that determines the influence or not of unstable areas. Within this geoenvironmental characterization, variables such as geology, geomorphology, land cover, and land uses are considered, which are relevant factors in the generation of mass movements [9,10,11]. For the correlation of these inherent variables with the areas that have historically generated a degree of instability (mass movements) mathematical methods are used where algorithms or simulations are used to estimate this degree of instability or susceptibility to mass movements, where its result is a non-temporal spatial vector element that considers the classification, volume (or area), and spatial distribution of mass movements [12,13,14]. These methods are grouped according to their scale, availability, quality, data accuracy, susceptibility resolution, and the results required in heuristic, stochastic, statistical, and deterministic methods [15,16]. According to the above, for a geoenvironmental characterization and its subsequent susceptibility model per MM at an intermediate level of detail (scale 1:25,000 or 1:10,000), stochastic and statistical methods are taken into consideration for the creation of the model [17]. The stochastic and statistical methods are proposed for the present study due to their performance in the spatial prediction of MMs and have been frequently studied in the construction of viable predictive models.
For the study area, mass movements in recent years have generated great effects on roads, bridges, and towns, especially in the winter season of 2010–2011 and 2019–2020. The occurrence of these events has manifested itself largely towards the northeast of the city of Bucaramanga and towards the southeast of the city of Piedecuesta. Where these regions have generated movements such as landslides and flows that have significantly affected the social and economic aspects of the region, even compromising the physical integrity of its inhabitants [18]. Part of this problem is due to a process of transformation of the land cover and land use in regions of urban expansion that are not suitable for the development of infrastructure projects of various types, mainly commercial and residential. Therefore, the evaluation of the territory from the present methods that integrate inputs available for free acquisition, such as the digital elevation model (DEM), base topography, and satellite images, which, in turn, serve as a basis in the construction of the variables for the generation of the susceptibility model, contributing to the knowledge to reduce and mitigate the impact of possible natural disasters in the near and/or distant future.
With this approach, an evaluation of the conditions of instability categorized by the calculation of susceptibility by mass movements would improve the spatial understanding of surface dynamics, giving clarity as a prediction of possible unstable environments. This is the main idea in the realization of the present study, which is based on the identification of these susceptible regions as a point of help to the control and monitoring of conservation practices and the good use of land cover to avoid an increase in unstable regions, added to an interest within urban planning. At the same time, it seeks to provide new data on how to build susceptibility models based on the correlation between variables within a regional framework at a scale of 1:25,000, this being an integral advance within regional and national policies in the knowledge of disaster risk management that are contrasted with the objectives set within the Sendai Framework for Disaster Risk Reduction, 2015–2030 [19].

2. Regional Setting

The study area is located north of the capital of Colombia (Bogotá), between the cities of Bucaramanga, Floridablanca, Piedecuesta, and Tona (Santander). This region has a topography between 682 m and 3920 m altitude, with average altitudinal of 2301 m (relative relief of 3238 m) and an ascending surface of 124 m per kilometer, with a very rugged topography. The study area has an area of 404 km2 and is subdivided into the following 4 hydrographic basins: Tona River (196.5 km2), Rio Frio (45.7 km2), Rio del Hato (31.1 km2), and Río de Oro (77.7 km2), added to regions without definition of hydrographic basin (52.8 km2) (Figure 1).
The region has a slope between 0° and 80°, with an average of 25.71° and dominance between the range of 23° and 32°, with surfaces that are mostly oriented towards the west–northwest (247.5°–337.5°). The drainage pattern is defined mainly subparallel, parallel to subdendritic (Figure 2). The region has a varied climate, from warm arid to humid cold, passing through humid temperate and humid high paramo (Institute of Hydrology, Meteorology and Environmental Studies—IDEAM, [20]). Where precipitation values ranging from 66 to 140 mm with temperatures between 23 °C and 29 °C are estimated for the city of Bucaramanga, which is at an altitude of 950 m, and precipitation values between 14 and 100 mm with temperature between 8 °C and 14 °C for the city of Tona, which is at an altitude of 1950 m [20].
In the geological context, the region presents old rocks of the Precambrian supereon that make up the basement, which is composed of schists, gneiss, and migmatites (PC). Later, in the pre-Devonian period, there were the following two units: PD1 and PD2. (Figure 3). The PD1 unit is in unconformity with the PC unit. This unit is made up of phyllites, schists, and quartzites in smaller quantities of metavolcanic-type rocks [21]. The PD2 unit presents quartz-monzonitic gneiss, granodioritic and tonalitic gneiss, diorite. For the Permian period (P), sedimentary rocks are formed in the conglomerate type with limestone pebbles. In the Triassic period, there were igneous rocks of hornblendic biotitic tonalite composition (TR1) and sedimentary rocks, such as sandstone and shale (TR2) [21].
In the Triassic—Jurassic period (JTR), igneous rocks were created with biotitic quartz-monzonitic composition, quartz-monzonite, and granite. For the Jurassic (JR), granitic rocks were formed, biotitic type alaskite and sedimentary rocks of conglomeratic sandstone, conglomerate, and siltstone type (J) developed [21]. For the Cretaceous period, three units are described (K1, K2, and K3). Unit K1 consists of sedimentary rocks with a composition of quartz sandstone, siltstone, and conglomerate sandstone. The K2 unit consists of fine-grained limestone, shale, and sandstone, and the K3 unit presents thin-bedded black shales [21]. Finally, for the Quaternary period, glacial, hydrogravitational, and alluvial deposits developed. The glacial deposits of tillitas (Qg), the hydrogravitational deposits of colluvial, hillside, and collapse type (Qc), and, finally, the deposits of alluvial type (Qal) [21]. In the structural context, the study area presents a control generated mainly by the stroke of the Bucaramanga Fault, which is a regional type of structure with a combined component in the course (sinestral) and in the dip (reverse) with an NNW–SSE direction [22,23]. There is also the Tona Fault, which is an inferred structure with a S–SSW direction that extends through the metamorphic rocks of the pre-Devonian, and the faults, Picacho and Sevilla that are inferred structures, parallel to each other and with a S–SSW direction that extend over the rocks of the pre-Devonian, Jurassic, and Cretaceous [21]. The Picacho Fault ends against the Tona Fault, while the Seville Fault extends to the Bucaramanga Fault. From a geomorphological point of view, the region is characterized by an abrupt relief, with straight and steep slopes, rounded and sharp blades, deep and rectilinear drains, a good development of the drainage network, and deep valleys in the form of “V”, typical of a relief with a strong structural control [24]. The region has very large areas where the interaction of processes, such as sedimentation in low areas and active erosion in mountain fronts, occurs [24].

3. Methodology and Materials

For the construction of the susceptibility map by mass movements for each proposed method, the respective cartographs will be elaborated at a working scale of 1:25,000 with a minimum cartographic area of 10,000 m2. To this end, a methodological scheme was proposed that easily describes the acquisition of a result that provides answers effectively and efficiently within a spatial distribution that is prone to generate unstable areas (Figure 4). Since understanding the dynamic behavior of these areas from the geoenvironmental characterization (variables), provides the opportunity to have a vital resource within the planning of the territory and risk management plans minimizing the impacts generated by the activation of areas with medium to high susceptibility. The inputs for the construction of the respective variables such as geology, morphology, morphometry, and morphodynamics for the calculation of susceptibility will be acquired from data from free databases at the national level. The base cartography will be obtained from the Agustín Codazzi Geographic Institute (IGAC) at a scale of 1:25,000 [25], in turn, the digital elevation model (DEM) will be implemented with a resolution of 12.5 m [26], which will be obtained from the portal ASF—ALOS PALSAR. The geological map will be acquired from the Colombian Geological Service—SGC, [27]. The land cover and layers of the forest–non-forest change relationship (temporary, evaluated years 1990 and 2019) will be obtained from IDEAM [28,29]. These inputs provide within each model a good resolution for the proposed interpretation scale. The construction of the layers and the calculations will be developed under the platform of the ArcGIS V10.8 program. This variable is constructed from a photointerpretation process that integrates satellite images (Google Earth®, Google, Mountain View, CA, USA), DEM, and base cartography.

3.1. Construction of Inherent and Dependent Variables

Following the scheme of Figure 4, for the geology attribute and its superficial geological unit (UGS) variable, is taken the description made by Hermelín, cited in SGC, [4], where the geological materials that emerge on the surface are a set (homogeneous) that come from the same origin and retain in general the same physical characteristics and geomechanical behavior of the parent rock, up to a few tens of meters below the surface of the ground. The UGS is classified according to its origin in rocks and soils. The rock is described by its degree of resistance according to its incidence to structural elements that control the mechanics of the rock massif and are classified as low, medium, high, and very high resistance [30,31]. Soils are defined as loose or soft deposits of natural origin formed on the earth’s surface by physical, chemical, and biological processes that act to produce a material rich in organic matter with characteristic horizons (layers) at shallow depths [30]. These soils are classified into the following: soils derived from rocks in situ by weathering process and transported soils that have been generated by deposits with different genesis [32].
For the geomorphological attribute, three variables will be constructed (morphogenetic, morphometric, and morphodynamic). The morphogenetic variable is elaborated from a photointerpretation process that integrates satellite images (Google Earth®), DEM, and base cartography at a resolution level of geomorphological subunit (SGMF) [33,34,35]. Its creation is based on establishing the causes and processes that modeled the relief from the endogenous and exogenous factors of the terrestrial dynamics, added to the degradation and gradation of the terrain product of the processes of weathering, erosion, and transport over time [36,37]. The correlation of these elements in the SGMF determines the morphogenetic environments that, in turn, for the present study are focused on estimating the degree of contribution in the generation of mass movements. For the morphometric variable, a processing will be carried out, which is based on a calculation that integrates the topographic base and the DEM, from this, different classes of variables will be obtained that will henceforth be described as parameters, which enter the calculation of the susceptibility model [38]. The parameters of the morphometric variable are selected for their degree of influence on slope instability processes [39,40]. From this, for the study area only the slope and rugosity will be considered due to their influence on the generation of mass movements, in addition to the fact that these parameters provide important information about topographic contrasts.
The slope parameter evaluates the relationship between the degree of inclination of the slopes and the type of rock or soil [41]. This relationship is given by the angle of friction, permeability, and cohesion between particles [42,43]. In turn, within the analysis of susceptibility by mass movements, the slope can be a determining variable to identify potential breaks or areas of instability since it is a geometric factor directly related to the tangential and normal shear stresses of the soil [44]. The rugosity parameter determines the variation of a surface described by the vectors calculated in a three-dimensional space (spherical distribution variable), where the slope and orientation of the terrain are taken as a basis [45]. This variable defines well the forms of the terrain, such as the limits of the slopes and hillside both in the valleys and in the crests; in turn, the areas of high rugosity are more prone to generate mass movements [45].
For the morphodynamic variable that is also described as the dependent variable corresponds to the inventory of mass movements present in the study area. This variable is a function of the geoenvironmental characteristics described above and is one of the most important elements in the calculation of susceptibility by mass movements since it gives an idea about the failure mechanisms, causative factors, frequency of occurrence, volume, and damage caused [46]. The construction of this variable is based on a photointerpretation process carried out on images arranged on the Google Earth® platform and that the study area has an optimal resolution to determine the characteristics of type, geometry, distribution, and material of the mass movements [47,48]. The definition and characterization of mass movements will be carried out as initially described by Varnes [49], and which is extended by PMA-GCA [50] and Hungr et al. [51]; in turn, the aspects described by Skempton & Hutchinson [52] for landslide and flow type mass movements, which can be classified as superficial and deep, are taken into account.
For the land cover attribute, the variables of current land cover and forest change—non-forest (temporary) will be generated. The current coverage variable is taken from IDEAM [29] for the year 2018 and its construction is based on the CORINE land cover methodology [53], where inputs that present a higher temporal and spatial frequency, generated with medium- and high-resolution images (Sentinel 1–2 and Planet Scope) are taken into account. This methodology describes the type of vegetation present at different levels depending on the scale of interpretation of the terrain, being level 3 and lower, which will be taken as a reference for this study [4]. The current land cover is an important variable in the study of susceptibility to mass movements since, depending on the type of cover, it offers better soil protection and contributes to the dissipation of rainfall energy [54]. In turn, it helps the terrain in the protection of erosion processes that can subsequently affect the generation of mass movements, especially on those surfaces that present an unfavorable slope in the stabilization of the hillside [55,56].
The temporal variable will be constructed from the layers of forest–not forest, which is an input generated by the IDEAM, and from which the years 1990 and 2019 will be taken [28]. This input was built through the aspects discussed in the Warsaw framework for reducing emissions due to deforestation and forest degradation (REDD+) (United Nations—UN, [57]), in addition to the aspects contained in the study carried out by Csiszat, et al., [58] on the changes generated in the forest surfaces through satellite monitoring, which is based on the medium- and high-resolution images Sentinel 1–2 and Planet Scope. The change in the cover evaluated in the present study by the relationship forest–not forest, can be given as a similarity of the relationship of change in cover of secondary forests to pastures and/or cropland, being this type of relationship a factor of increase in areas prone to mass movements [59,60]. These changes generate effects on the stability of the slopes due to a variation in the hydrological and mechanical dynamics of the environment. Within the mechanical factors are described the reinforcement of the soil due to the adhesion of the roots, the erosion control, and the resistance of the soil soil-roots system [61,62]. The hydrological effects are related to the occurrence of mass movements due to the low interception of rainfall by land cover, the loss of soil moisture by evapotranspiration, and the effects on hydraulic conductivity [55,61].

3.2. Methods for Calculating Susceptibility by Mass Movements

From the geoenvironmental characterization of the inherent and conditioning elements, the calculation of susceptibility by mass movements for landslide and flow events will be generated. This calculation is made by means of two methods (stochastic and statistical), which are in accordance with the scale of study and the degree of detail required according by Soeters [16], Guzzetti [63], Cascini [17] and Corominas et al. [46]. The stochastic method is based on variables taken as random data, which are related to each other with probabilistic functions [64,65]. Within the existing stochastic methods for the calculation of susceptibility, the artificial neural network (ANN) algorithm is taken due to its versatility in the generation of this type of models [66,67,68]. Under the ANN algorithm, the backpropagation (BP) function is taken, which presents a better performance in the solution of this type of problems [69,70,71].
For statistical methods, its basis is to establish combinations between the conditioning elements and the mass movements generated in the past in a statistical way, in the same way, it creates a quantitative predictive model on the areas that do not present mass movements, but that have similar conditions [16]. Within the statistical methods, there is bivariate and multivariate statistical analysis. Bivariate statistical analysis is based on the Bayesian theorem, which indicates the probability that event A occurred because event B occurred [72,73]. With this relationship, a combination of the densities of occurrence of these mass movements and the maps of each of the conditioning factors associated with the movements is sought, based on a spatial relationship, distribution of the areas affected by mass movements and distribution of the conditioning factors that cause the movements [4,74,75].
Multivariate statistical analysis studies and examines the set of statistical techniques that analyze the simultaneous effect of multiple variables. Multivariate statistical methods are used to analyze the aggregate behavior of more than one random variable [76]. Within the multivariate statistical methods, there are two that have been widely used, discriminant analysis (AD) and logistic regression (LR) [77]. For the present study, the logistic regression algorithm is taken due to its versatility in the simulation of input data that can be quantitative (continuous or discrete), qualitative, or a combination of both [78,79]. This algorithm has shown good results in the generation of susceptibility maps by mass movements [78,80,81,82]. In addition, this algorithm establishes relationships to find the model that best fits the relationship between the presence or absence of mass movements and the variables that condition it, thus evaluating the susceptibility of the area in terms of probability of occurrence to a hillside movement [83]. For the calculation, the studies and guidelines proposed by Gómez et al. [82], Aguayo [84] and Aguayo & Lora [85].

3.3. Validation of the Susceptibility Model

The validation of the models will be carried out from the construction of the receiver operating characteristic (ROC) curve, which is the graphical representation of the sensitivity to specificity for a binary classifier system as the discrimination threshold is varied, and from this its area under the curve (AUC) is estimated as a statistical measure of the success and prediction rates of each model [86,87]. The results obtained by each of the methods (stochastic and statistical) will be expressed in three categories (high, medium, and low) according to their degree of incidence present by the correlation of the variables and the unstable points in a non-temporal spatial condition. For the correlation of the data, a partition is made in percentiles of 100 categories, each category of the 1% is crossed with each other to estimate its degree of relationship in percentage [86]. From this correlation, two results are obtained in percentage, one that measures the goodness of fit of the model and the other its degree of prediction. The latter is measured by the crossing between the frequency of values of the calculation against the frequency of values of 50% of MMs that were selected for validation and that are not part of the initial calculation of the model. The purpose of this analysis is to observe if each model predicts coherently about the areas that did not enter the calculation (MM) and their degree of confidence, which must present an acceptable percentage of adjustment greater than 70% in both cases (test and validation) [4]. With the behavior obtained from the statistical correlation performed, a selection of the percentiles is made to establish the three categories of susceptibility following the criteria raised by van Westen, [74] where the percentile that groups > 75% of mass movements for high susceptibility, 25% of MMs for medium susceptibility, and <2% of MMs for low susceptibility is determined, this grouping being carried out in each method proposed.

4. Results

4.1. Description of Geoenvironmental Elements

For this stage of processing the susceptibility model, a correct correlation and correspondence in scale and resolution must be obtained. According to Casini, [17] the evaluation of susceptibility within a scheme of urban planning and zoning of regions may lead to a problem with anthropogenic factors that must be subject to large-scale cartographic elements. As well, it is relevant to establish a parameter of a minimum cartographic area that ensures a geometric construction relevant to the work scale, this is in order to ensure that the cartographic content of each variable preserves the criteria with which they were generated [88]. To this end, a minimum cartographic area was defined for the 1:25,000 scale of 10,000 m2 that is applied to each attribute generated for the variables of geology, morphogenesis, morphodynamics, and land cover. For the UGS variable, a cartographic adjustment was made through a process of photointerpretation of the input obtained from the UGS. As a result, for the study area, a total of 35 UGS were characterized with the parameters described above. These units were classified into rocks with very high strength presenting a total area of 103.5 km2 (25.7%), high strength with 151.8 km2 (37.7%), medium strength with 107.1 km2 (26.6%), low resistance with 9.6 km2 (2.4%), residual soil with area of 15.4 km2 (3.8%), and transported soil with an area of 15.4 km2 (3.8%). For the SGMF variable, a total of 38 SGMF was obtained through photointerpretation. These units were classified as follows: denudational ambience with a total area of 332.1 km2 (82.3%), structural, 60.9 km2 (15.1%), glacial, 6.8 km2 (1.7%), fluvial, 3.4 km2 (0.8%), and anthropogenic, 0.6 km2 (0.1%).
For the slope, the ranges in which mass movements are more affected were classified [89]. From this, a statistical analysis was generated where the standard deviation of the data is taken as a reference to obtain the classification of the ranges of the slopes. The present statistical classification divided the slope into 8 classes, each one, with a range of 10° (Table 1). Although it is observed that the slope for the study area has a wide range (min.: 0° and max.: 74.8°), when evaluated with mass movements, the largest amount of area in which the unstable zones occur is between 30° and 40° (39.05% of the MM area). These areas are described as being steep to very steep slopes where there are entrenched processes of the denudational type, with a high possibility of the development of erosive processes where chaotic granular deposits of little thickness can be generated [90].
For the rugosity in the study area, it is described that surfaces with a high rugosity present the greatest amount of mass movements (40% of the MM area), followed by surfaces with a medium rugosity (38% of the MM area), a condition consistent with what was described above. For the morphodynamic variable (MM inventory) of the study area, a total of 856 MMs were mapped, with a total area of 1.72 km2 representing 0.4% of the total area of the region (max. area: 49,903.4 m2; min. area: 33.5 m2 affected by an event), with a frequency of 2 movements per km2 (Figure 5). These movements are described as, traslational debris landslide with 379 events (1.09 km2, 63.3%), traslational earth landslide with 205 events (0.17 km2, 10.1%), debris flow with 113 events (0.24 km2, 14.1%), and earth flow with 159 events (0.21 km2, 12.5%). It is also described that 364 events (0.39 km2, 22.6%) are of surface type and 492 events (1.33 km2, 77.4%) are of deep type (>2 m thick). Performing an analysis of the basins, it is observed that the Tona River basin presents the largest number of events (323 MM, 0.51 km2), followed by the Rio Frio basin with 150 MMs (0.26 km2), Río de Oro basin with 95 MMs (0.15 km2), and Río del Hato basin with 78 MMs (0.07 km2); furthermore, the zones that do not have a representation of a hydrographic basin (W/B) present a total of 210 MMs (0.73 km2), being the zones that present the greatest extension of mass movements over the entire study area. With the present inventory, a selection of MMs is made for the calculation of susceptibility (test) and its subsequent validation. This condition is established in order not to generate trends or biases within the results obtained and are applied for each method. To this end, 50% MMs are selected for test and 50% for validation following the following premises:
  • Selection of MM randomly;
  • Randomly selected MMs should have equivalent areas in both test and validation;
  • The MMs must be spatially distributed;
  • Within 50% of each selection, you must have 25% of each type of MM described (4 types of MM).
Figure 5. Spatial distribution of mass movements presents in the study area.
Figure 5. Spatial distribution of mass movements presents in the study area.
Remotesensing 15 04567 g005
On the other hand, for these MMs, there is a correlation between them with respect to the detonating source since, from the geometric construction and its characterization, a clear pattern associated with the rainfalls is distinguished, where 693 MMs (81%) of these detonating sources could be established. The other 163 MMs (19%) are reported in databases, but do not have a specific description or the condition of the triggering factor associated with the event. From this condition and added to the characteristic of the speeds of displacement of the different types of movements described by Cruden and Varnes [91], the mass movements such as landslides and flows are grouped for the calculation of the susceptibility of each model.
For the variable of land cover that was acquired from IDEAM for the year 2018, it is taken as a direct variable for the calculation of susceptibility due to the way it was built. The following three aspects were considered for its use: 1. The coverage described for the year 2018 has not undergone considerable changes due to anthropogenic factors in the last decade. 2. The input is sufficiently detailed for the susceptibility calculation. 3. The extension of each of the generated coverages does not simplify the susceptibility function for mass movements. With these aspects, we have a layer with 26 covers where the coverage in high secondary vegetation presents the largest extension affected by mass movements (0.38 km2, 22.06%), followed by the mosaic coverage of crops, pastures, and natural spaces with an extension 0.27 km2 (15.54%) and the high-density forest of the mainland with an extension of 0.25 km2 (14.6%).
The temporal variable describes the behavior of change generated between the years 1990 and 2019 (29 years) for areas where tree cover predominates with a minimum canopy density of 30% and a minimum canopy height in situ of 5 m, with a minimum area of 0.01 km2 (forest layer). The areas with coverage other than natural forest are called non-forest areas (non-forest layer) [92]. With the forest–non-forest layers for each year, a cross is made between them, and by means of a relationship in a combination matrix, the classes of persistence, gain, and request of the cover are determined. This relationship allows us to identify the influence of the time variable as a function of mass movements for a period of 29 years. From this relationship it is observed that the greatest amount of mass movements are found in areas that present a persistence of the type of cover (forest and/or non-forest) with an area of 1.34 km2 (78%), followed by areas with a loss of cover (forest to non-forest) with an area of 0.35 km2 (20.3%) and, finally, areas that present a gain in the type of cover (not forest to forest) with an area of 0.03 km2 (1.7%). From this analysis, it would be expected that the areas with a loss of cover would have the greatest amount of mass movements, but this condition generated over areas of persistence may be associated with intrinsic factors of the terrain other than the type of coverage. These factors can be related to the inclination of the land, the type of rock present, the supersaturation of the ground by infiltration, or the structural control that can exert a fragile behavior, giving rise to a degradation of the rock massif by the geomechanical aspects of the rocks. Another important aspect to consider is the variation of forest–non-forest in a period of 29 years, where there is a loss of approximately 60.4 km2 (2.08 km2/year). With this rate of change and without variations, an extrapolation is made based on the year 1990 to which it is projected with a loss of 50%. This scenario would be generated by the year 2052 and with a loss of 100% by the year 2112. Based on this analysis, the conditions of instability of the terrain can be extrapolated to the future since the regions described as lost present 20.3% of the area in mass movements, a condition that can increase in a critical way.

4.2. Calculation of Susceptibility

From the geoenvironmental characterization carried out, the calculation of susceptibility by mass movements is generated considering the conditions described in Section 3 and following the methodological scheme proposed (Figure 4). The first susceptibility model was performed by means of the stochastic method under the algorithm of ANN in function BP. For the analysis of the ANN model, a total of 20 architectures were built, where, for each of these architectures, the parameters of the modeling percentage, test percentage, validation percentage, number of neurons, and the type of training algorithm were modified. This is in order to estimate by way of training and its performance in the degree of prediction is measured by the result of the quadratic regression, performance curve, training status, and confusion matrix, thus defining the best architecture in the susceptibility model. As part of the processing to obtain a better performance, the variables were normalized to guarantee a correct correlation between them and the points where mass movements are presented [71]. Of the 20 architectures, the 6-50-1 architecture (6 variables, 50 neurons, and 1 output) was taken, which has a configuration of 75% data for training, 15% data for validation, and 10% data for testing, in a Levenberg–Marquardt learning algorithm. An important aspect to mention is that, in all simulations, in their training, a standard arrangement was left in the number of iterations (1000), a goal equal to zero (0), and a maximum of failures of 25, in turn, the values taken for training, validation, and testing were distributed internally by the program in 50% of the data of the variables and MMs. These percentages are part of the training form of the algorithm, thus generating a second selection of data, which guarantees a greater degree of confidence when reviewing the goodness of fit of the model [71]. This selection of percentages can be generated manually (expert criteria) or by the predefined configuration of the program. For the present study, it was decided to leave the values predefined by the program since the data entered in the simulation only represent 50%, and making a decrease in the data in the training does not guarantee an ideal model. The 6-50-1 architecture was taken due to its performance in each of the results, which present quadratic regression values equal to 0.94265 and an overall accuracy of 94.2% within the simulation, added to its performance in training and performance, conditions that in the other architectures describe considerable variations (Figure 6).
The result of the simulation of the 6-50-1 architecture under the ANN algorithm is processed by GIS for conversion of a numerical factor (ASCII) to a vector characterization and later for classification by categories. This classification is performed by means of the statistical construction by percentiles described above. From this result, it is obtained that the high category of susceptibility presents an area of 139.41 km2 (34.52%) of the total study area, and on this area, 501 MMs were grouped, which generate an affected area of 1 km2 (58.57%). The medium category of susceptibility has an area of 179.06 km2 (44.34%), groups a total of 309 MMs with an area of 0.63 km2 (36.14%). Finally, the low category has an area of 85.32 km2 (21.14%) and groups a total of 46 MMs, and an affected area of 0.09 km2 (5.29%) (Figure 7). The final susceptibility model represented by the three categories was validated by the ROC curve with a 0.695, describing an optimal goodness of fit.
For the present model, a close relationship is observed between the variables that entered the simulation and the regions that were categorized into the three susceptibility classes. Although, in a specific way, the variables of the slope and rugosity show a better correlation with the areas that present instability and variables such as coverage and temporal present a moderate to high agreement. On the other hand, the ROC (0.695) is below the standards (0.7 or 70%), but the vector characterization and the number of MMs (50% test) that are presented in the high category give an acceptable goodness of fit, a condition that will be evaluated later with the validation of all models with the other 50% of MMs. With the present result, a description was made regarding hydrographic basins in order to determine the percentage of affectation in each basin. To this end, it is observed that, by limiting the model to the basin, regions such as the Rio Frio basin describe a better relationship between the category and the number of MMs (evaluation in high category), on the contrary, the Tona River basin presents a low prediction in the category of high susceptibility and its relationship with the MMs is moderate (Table 2). The reason for these changes in the form of prediction will be described in more detail at the end of the three methods of susceptibility assessment to have a more complete picture.
For the susceptibility model generated by means of the statistical method, the calculation of the bivariate statistical analysis is performed first. Due to the form of construction of the present model, it is not possible to establish a relationship in which a variable presents more or less influence on the final result, as described for the ANN, since the bivariate statistical analysis works according to the areas, a weight is established in each attribute within the variable. But you can make a description of the attribute that can give a greater or lesser weight in the variable. As a result, for the present model, it is observed that the high category of susceptibility has an area of 105.65 km2 (26.25%) of the study area, and on this were grouped 591 MMs with an affected area of 1.19 km2 (69.21%) of the total area of MMs. The medium category of susceptibility presents an area of 164.65 km2 (40.91%), groups a total of 225 MMs, with an area of 0.45 km2 (26.23%). Finally, the low category of susceptibility has an area of 132.16 km2 (32.83%), groups a total of 40 MMs with an affected area of 0.08 km2 (4.55%) (Figure 8).
Within the evaluation of the weights of evidence, it is observed that, the slope variable in the range of 50°–70° contributes to the occurrence of MMs, and the slope between 0° and 20° does not contribute to the occurrence of MMs. An important aspect to highlight in this evaluation is the range of the slope between 50° and 70° since, as mentioned, the present method works on the relationship of the areas and the MMs. Observing this range of the slope, we have a region that presents very little extension, but with a large number of MMs. This is mentioned since in the literature there is a description of a critical range that generates instability on the slopes, which is between 30° and 40° [93,94], which is consistent with what is described in Table 1, but differs from the weight of evidence generated by the bivariate method.
The other contributing variables are (UGS) low-resistance shale rock (PD1) and quartz-monzonite (JTR), (SGMF) very steep erosive slope and pressure loin, (land cover) rocky outcrops and rounds of water bodies of urban areas, (temporary) loss of coverage, and (morphometric) high rugosity. The variables that do not contribute are (UGS) the rock of moderate sandstone resistance (K1) and the glacial transported soil (Qg), (SGMF) the structural slope and the glacio-fluvial cones, (land cover) open shrub and dense grassland with shrubs, (temporary) coverage gain, and (morphometric) low rugosity. From the previous description, it is observed that all the attributes are consistent with the conditions of instability of the surfaces, making an exception of the slope. The final susceptibility model presents a ROC equal to 0.822 that describes a very good goodness of fit. Performing an analysis of the basins, it is observed that the W/B regions describe the best relationship of the category against the number of MMs (evaluation in high category), followed by the Rio Frio basin, on the contrary, the Tona River basin presents a low prediction in the category of high susceptibility and its relationship with the MMs is low (Table 3).
Finally, the susceptibility model is made from the multivariate method created under the LR algorithm. For this model, a total of 16 simulations were generated by modifying the conditions of the dependent and covariate variables, the amount of data entered into the model, and the shape of the data (normalized or in native format). For each simulation performed, the behavior of Cox and Snell’s R2, Nagelkerke’s R2, the classification percentage, and the results of the classification table as a function of Wald and Sigmoidal were observed [81,82]. Due to the method of classification of the dependent and covariate variables, which does not comply with the assumption of normality in the distribution of the values, a nonparametric test was performed to observe the degree of correlation that extends between the data. For this analysis, the data were evaluated by means of the Spearman correlation to measure the relationship between the two groups from the arrangement generated for each of the 16 simulations (Table 4).
With each arrangement and its simulation, the one that obtained the best performance in the degree of prediction was evaluated and taken, being the arrangement that presents 5 unnormalized variables, 3 of ordinal type (UGS, land cover, and temporary coverage), and 2 of scalar type (slope and rugosity), with values of Cox and Snell R2 equal to 0.240 and the Nagelkerke R2 equal to 0.321, a classification percentage of 69.4%, Wald values greater than zero and sigmoidal without the presence of values greater than 0.5, and data obtained after the last correlation generated. This simulation, in its first Spearman correlation test, showed that the SGMF variable did not present a good correlation since it obtained a result of 0.03 below the standard (0.05), so it was not taken for the final simulation. From the result of the present model, it is obtained that the high category of susceptibility presents an area of 122.3 km2 (30.37%) of the total study area, and on this area, 257 MMs were grouped with an affected area of 0.50 km2 (29%) of the total area of MMs. The medium category of susceptibility presents an area of 146 km2 (36.28%), groups a total of 325 MMs with an area of 0.69 km2 (40.32%) of the total area of MMs. Finally, the low category has an area of 134.2 km2 (33.34%), groups a total of 274 MMs with an affected area of 0.53 km2 (30.68%) of the total area of MMs (Figure 9). The final susceptibility model has a ROC equal to 0.530, which describes a good goodness of fit. With the present susceptibility model, an analysis is carried out by hydrographic basins, from which it is observed that regions such as the Tona and Río de Oro rivers basins describe the best relationship of the category against the number of MMs (evaluation in high category), on the contrary, the Río Frio basin presents a low prediction in the category of high susceptibility and its relationship with the MMs is low (Table 5).
With the results generated in each of the models proposed for the calculation of susceptibility by mass movements over the study area, considerable contrasts are observed in the prediction of spatial form. Based on the fact that the calculation is given by the interaction of the variables or inherent conditions of the area with the points where ground movements have been generated, in order to extrapolate these scenarios to regions where no events have occurred but that may eventually be unstable in the future [71,94]. These models describe different ways of relating these two conditions, as mentioned above, but their result at the level of spatial analysis determines those sectors that are susceptible to changes, but at a level of procedural control, such as vector construction, they differ from each other. The final models of each method were subjected to their own tests given for each algorithm and the program executed, thus ensuring the selection of the simulation that best represents the relationship between the conditioning factors and the dependent variable; however, each of the final models was given three tests to contract the effectiveness in predicting the data. These three tests are part of the validation in which 50% of the mass movements that were not entered in the initial calculation (creation of susceptibility model for each proposed method) are taken as a basis. These movements were selected according to the specifications described above. The purpose of the present validation was to observe if the model really predicts on the regions that were not entered in the execution of the algorithm. The first validation test was based on the ROC curve and its AUC. The second test was based on comparing the agreement observed between each set of data obtained between the test and the validation of each model, for this comparison the statistical measure of Cohen’s kappa coefficient was taken. The third test consisted of the degree of classification generated from a confusion matrix, and from this, the values of accuracy, precision, recall, and harmonic mean (F1-score) were obtained. For the construction of the confusion matrix, 100% MM (856 MMs) was taken, which was contrasted against the classification of susceptibility represented in three categories for each model. For this process, MM data and susceptibility categories are taken by areas represented in pixels (Figure 10).
For the susceptibility model created under the ANN algorithm, AUC = 69.5% (test) and AUC = 61.7% (validation), Kappa index (0.592—moderate), accuracy = 63.86%, precision = 0.67, recall = 0.89, and F1-score = 0.76 are described. In the model created under the bivariate algorithm, values of AUC = 82.2% (test) and AUC = 76.9% (validation), Kappa index (0.718—good), accuracy = 73.76%, precision = 0.77, recall = 0.93, and F1-score = 0.84 are described, and for the model generated by the RL algorithm values of AUC = 53.1% (test) and AUC = 51.1% (validation), Kappa index (0.951—very good), accuracy = 59.68%, precision = 0.47, recall = 0.78, and F1-score = 0.58 are described. With the construction of the AUC curves for each method, the Kappa index and the analysis from the confusion matrix and its results, it is observed that the calculation of susceptibility through a bivariate statistical analysis presents a good performance in the prediction both at the level of the calculation by each algorithm, and in the unstable regions that did not enter the model. An important aspect to highlight is that all simulations were performed with 50% of the MMs, which shows that the results present here have a high predictive power independent of the method executed. On the other hand, it is clear to mention that, although the ANN generates an estimate more consistent with this type of prediction models, the correlation of data by means of the bivariate method with weights of evidence simulated the points not entered into the model (test) in a better way. This condition may be associated with the cartographic construction of the variables, which, although these are according to the scale of study, the quantity and differentiation can generate non-optimal trends for the ANN and LR methods.
Once the result of each model from the vector construction was seen, the calculation performed by the bivariate method preserves the spatial relationship of both the mapped events and the areas that have a low susceptibility to mass movements since these regions have surfaces that have a slope <10° and UGS of the hard rocks type little weathered. In general, from the considerations described above, a prediction and quality of the model can be established regardless of the method used, with only 50% of the MM inventoried, but with the clarification that the variables that are really affecting the study area and their correct categorization must be adequately defined.

5. Discussion

The scenarios simulated from each method proposed to establish the degrees of susceptibility of a territory to the possibility of collapse of a surface show these conditions of spatial distribution from the relationships between the inherent elements of the terrain and their correlation with the instability events that have been triggered by some phenomenon of internal or external type of the earth [9,10]. This is seen from a conception of temporal space or hazard type since they are factors of great importance at the time of establishing a territorial zoning that helps to delimit regions where natural and anthropic environments are exposed and vulnerable. From this perspective, and only taken in account the first step for the construction and analysis of a territorial zoning, the susceptibility to mass movements as a spatial element describes those environments where land use planning at an adequate scale will determine those potentially susceptible sites so that people as decision makers and rulers of the day can have tools and basic inputs at the time of territorial planning [95].
Within this context, a correct geoenvironmental characterization made from each defined variable and its corresponding description by attributes that are integrated and modeled by stochastic and statistical methods for this particular study is an important basis for these structures of the territorial organization since it is described in a quantitative and qualitative way within the cartography of an environment, regions with the possibility of instability. However, the construction of a susceptibility model depends to a large extent on expert judgment and its scale of interpretation since the mathematical methods adopted work based on the data entered and their corresponding correlation of the information. However, a model can be viable and have the same degree of perdition if the geometry of the study area is modified, where the same number of variables generated is correlated to the established work scale. For this case, an exercise was carried out where the variables previously constructed were taken but limited to the geometry of the basin and modeled by means of the bivariate method. Subsequently, with this idea, a model of susceptibility by mass movements was generated to observe the prediction behavior applying the conditions described for the previous models but based on the area of the basin (4 basins) of the study area, applying the bivariate method. This method was selected due to its good development in the calculation of susceptibility in the previous chapter, added to its versatility and correspondence between the independent or inherent variables and the dependent. The procedure to perform this test is based on the considerations raised in Section 3 and described in Section 4 for the construction of the inherent variables and the dependent variable that were limited to the geometry of the river basin. For the dependent variable, 50% of MMs was selected for the test and 50% MMs for validation in each basin. This calculation was made to observe the behavior of the susceptibility model created regionally and later contrasted with the model generated by the river basin. With the above, the calculation of susceptibility and its subsequent correlation with the results obtained in the study area was carried out (Figure 11).
From the data obtained from this evaluation, the following observations are described:
  • The results obtained from the susceptibility calculation carried out in each basin have a prediction greater than 70% (AUC curve) in the test phase. This demonstrates a good development in the spatial prediction of regions with different susceptibility ranges. As it is also appreciated that the AUC curves of the validation present a percentage higher than 60%, indicating that the areas with the events that were not included in the calculation of the test were classified correctly;
  • Within the calculation of susceptibility and taking into consideration the results obtained from the AUC curves, a condition can be established in the construction of the inherent variables, which should not necessarily be restricted to the geometry of the hydrographic basin. Since, generally, susceptibility models are limited to this geometry. From the analysis carried out in the present study, the regions of the model can be expanded by duly safeguarding their limits in scale and level of interpretation. This can be established if you have your own knowledge of the region evaluated and experience in the field of study (expert criteria). As well, if variables such as geology, geomorphology, and land cover, do not present a change of more than 20% in their attributes over the areas that enter the calculation of the susceptibility model;
  • From the present study, it is also indicated that the hydrographic basins that have a surface area less than 50 km2 do not generate an ideal model of susceptibility. In these cases, it is convenient to make a model with a greater extension. But with the dimension that must be considered with those described in the previous point and the calculation must be performed with the algorithm of bivariate statistical analysis.
One of the most relevant conditions between the two approaches to the calculation of susceptibility is the variation that occurs between the result of the categories between the total area (regional model) and the evaluation by river basins. This variation is manifested in the form of prediction of each of the categories that presents a high variability in its distribution by extension (area), which is described in Table 6. On the one hand, areas with a low susceptibility category are being underestimated, as shown in the Tona River basin, since, in the calculation generated for the regional model, it presents a greater extension in area. This condition is also evidenced in the average susceptibility of the Rio Frio, Rio del Hato, and Rio de Oro rivers basins. On the other hand, it overestimates the low susceptibility category in the Rio Frio and Rio del Hato rivers basins, which greatly increases the areas of this category in these basins. In turn, the relationship of the spatial distribution of the susceptibility categories in the model by basins tends to generate a concentration and homogenization of the categories in certain areas, as evidenced southeast of the Rio Frio basin, on the high susceptibility category (Figure 11). To understand this aspect, it must be established that the construction of the variables for the model by basin part of the cartography was carried out for the regional model and was simply limited to the geometry of the basin. Therefore, although there is a variation in the size of the area for the evaluation of the weights of evidence between the two models (regional and by basin), the variability of the categories should not be so considerable. This aspect is also evident in the central region of the Río de Oro basin, where it homogenizes the present area to a single (low) category. As mentioned above, these aspects where variability, overestimation and/or underestimation, concentration, and homogenization are generated in the evaluation of the categorized regions of susceptibility must have an important control from the point of view of the expert criterion since it can induce incorrect interpretations about the model. In addition, a model can be very good as long as it is on the standards, but its interpretation and implementation must be subject to regional conditions and experts who have knowledge of the area evaluated and its geoenvironmental elements.
On the other hand, regarding the evaluation carried out of the susceptibility that has been well-described as a spatial element of the inherent factors of the terrain, it can be extended to a spatiotemporal or hazard level, where the triggering factors are included in temporal terms. In this context, the study by Valencia Ortiz et al., [96] which was focused on this region and evaluated the conditions of rainfall and seismic activity as triggers of the MMs. For this study, a characterization of the rainfall was carried out, obtaining a threshold of 53 mm for 24 h rainfall and 158 mm for accumulated rainfall of 15 previous days. These thresholds have a high probability of being generated for the months of May and October in a return period of 25 years. For seismic activity, the present study mentions that events of magnitude ≥ 4 Mw may be generating unstable areas that are subsequently detonated by rain in a combined action. From this, the descriptions made of susceptibility models that are in a spatial function can be characterized spatiotemporally, where the regions valued with a medium or high susceptibility would have a high probability of being detonated with the aforementioned characteristics. Based on this, with the information described in the study by Valencia Ortiz et al., [96], added to the evaluation carried out in this study on the susceptibility by mass movements, there would be a more robust tool to support policies and objectives within the risk management and urban planning of the study area since the combination of these two aspects would indicate the areas susceptible to collapse and their triggering factor in terms of temporal probability; in turn, it can be an input for the design of emergency plans and early warnings, added to an environmental diagnosis on the conditions of the surfaces.

6. Conclusions

With the methods proposed to determine the categories (low, medium, and high) of susceptibility by mass movements, it was observed that the bivariate statistical analysis algorithm presents the best performance in the spatial prediction of instability events. This algorithm reached a degree of measurement prediction of the ROC curve and its AUC = 82.2% (test) and AUC = 76.9% (validation), grouping a total of 591 MM of the 856 MM in the high category (69%). In turn, it is described that the basins of the Rio Frio and Rio del Hato present the highest percentages of susceptibility in the high category evaluated with the bivariate method. This result was obtained from the correlation of the variables of the surface geological unit (UGS—geology), geomorphological subunit (SGMF—morphogenesis), slope, rugosity, land cover, and temporal cover (1990–2019) against the inventory of mass movements. These results were achieved with a model that takes 50% of the MMs to generate the test and 50% of the MMs for validation, and its construction was carried out by taking the MMs randomly and with a total distribution in an equivalent area. On the other hand, from a second susceptibility model generated by the bivariate statistical calculation on each river basin, it was obtained that this presents a high degree of variability with respect to the regional model, being the regional model the one that best describes the prediction of the categories of susceptibility by MM. The calculation made by river basins for the present study tends to generate variability, overestimation and/or underestimation, concentration, and homogenization of susceptibility categories in certain areas. This behavior is due to the shape of the variable and its construction of the attributes since if the interpretation scale for both models is considered (1:25,000), some attributes have very large surfaces that generalize the regions and skew the values or weights assigned. In order not to make a generalization of the results, which leads to an incorrect interpretation of the generated model, from the present study, it is established that at a level of analysis 1:25,000, hydrographic basins that present an extension of less than 50 km2 should be evaluated by increasing their study area. This extension in the area is feasible if it is taken into consideration that the variables within the study area do not present changes greater than 20% in the number of attributes evaluated, added to a knowledge of the region. The latter is conditioned to a model generated by the bivariate statistical analysis method since the evaluation with other methods is subject to the results and their degree of prediction.

Author Contributions

Conceptualization, data collection, data analysis, writing draft and final manuscript, generation of geology and land cover maps, calculation of susceptibility by neural networks and bivariate, J.A.V.O. Concept, writing draft manuscript, A.M.M.-G. Generation of the geomorphology map, calculation of the susceptibility by logistic regression and creation, revision, and adjustments of the document, L.M.M. All authors have read and agreed to the published version of the manuscript.

Funding

Grant 131874B-I00 funded by MCIN/AEI/10.13039/501100011033.

Data Availability Statement

Different programs were used to generate each of the susceptibility models. For the multivariate method, the IBM SPSS Statistics 26 program was used. For the artificial neural network method, the MatLab R2022a program was used. For the bivariate model, the Microsoft 365 program (Excel template) was used. Finally, the integration of the models and graphical outputs were generated using the ArcMap V10.8 program. All these programs under academic license of the University of Salamanca.

Acknowledgments

This research was assisted by the GEAPAGE research group (Environmental Geomorphology and Geological Heritage) of the University of Salamanca. The first author thanks the Foundation for the Future of Colombia (COLFUTURO) for the scholarship for doctoral stay.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. World Health Organization. Landslides. 2019. Available online: https://www.who.int/health-topics/landslides#tab=tab_1 (accessed on 25 June 2022).
  2. EM-DAT. Centre for Research on the Epidemiology of Disasters. Available online: https://www.emdat.be/ (accessed on 25 June 2022).
  3. Froude, M.; Petly, D. Global Fatal Landslide Ocurrencen 2004 to 2016. Nat. Hazards Earth Syst. Sci. Discuss. 2018, 18, 2161–2181. [Google Scholar] [CrossRef]
  4. SGC. Guía Metodológica para la Zonificación de Amenaza por Movimientos en Masa Escala 1:25.000; Servicio Geológico Colombiano: Bogota, Colombia, 2017. [Google Scholar]
  5. DesInventar.org. Disaster Information Management System. Available online: https://www.desinventar.org/ (accessed on 10 April 2022).
  6. Ingeominas. Catálogo nacional de movimientos en masa. In Subdirección de Amenazas Geoambiéntales; Ingeominas: Bogotá, Colombia, 2002. [Google Scholar]
  7. SGC. Las amenazas por movimientos en masa de Colombia, una visión a escala 1:100.000. In Grupo de Evaluación de Amenaza por Movimientos en Masa; Servicio Geológico Colombiano: Bogotá, Colombia, 2017. [Google Scholar] [CrossRef]
  8. Cruden, D.M. A simple definition of a landslide. Bull. Int. Assoc. Eng. Geol. 1991, 43, 27–29. [Google Scholar] [CrossRef]
  9. Popescu, M. Landslide Causal Factors and Landslide Remedial Options. In Proceedings of the 3rd International Conference on Landslides, Slope Stability and Safety of Infra-Structures, Singapore, 11–12 July 2002; pp. 61–81. [Google Scholar]
  10. van Westen, C.J.; Castellanos, E.; Kuriakose, S.L. Spatial data for landslide susceptibility, hazard, and vulnerability assessment: An overview. Eng. Geol. 2008, 102, 112–131. [Google Scholar] [CrossRef]
  11. Liang, L.; Tao, C.; Jie, D.; Antonio, P. A hybrid ensemble-based deep-learning framework for landslide susceptibility mapping. Int. J. Appl. Earth Obs. Geoinf. 2022, 108, 102713. [Google Scholar] [CrossRef]
  12. Fell, R.; Corominas, J.; Bonnard, C.; Cascini, L.; Leroi, E.; Savage, W. Guidelines for landslide susceptibility, hazard and risk zoning for land-use planning, on behalf of the JTC-1 Joint Technical Committee on Landslides and Engineered Slopes. Eng. Geol. 2008, 102, 85–98. [Google Scholar] [CrossRef]
  13. Martinez-Graña, A.M.; Goy, J.L.; Zazo, C. Geomorphological applications for susceptibility mapping of landslides in natural parks. Environ. Eng. Manag. J. 2016, 2, 15. [Google Scholar] [CrossRef]
  14. Rui, L.; YueKai, D.; Deliang, S.; Haijia, W.; Qingyu, G.; Shuxian, S.; Mingyong, L. Insights into spatial differential characteristics of landslide susceptibility from sub-region to whole-region cased by northeast Chongqing, China. Geomat. Nat. Hazards Risk 2023, 14, 25. [Google Scholar] [CrossRef]
  15. Soeters, R.; van Westen, C.J. Slope instability recognition, analysis and zonation. In Landslide Types and Processes; Special Report; Turner, A.K., Schuster, R.L., Eds.; National Research Council 247, National Academy of Sciences: Washington, DC, USA, 1996; Volume 247, pp. 129–177. [Google Scholar]
  16. Huanga, F.; Cao, Z.; Guo, J.; Jiang, S.; Li, S.; Guo, Z. Comparisons of heuristic, general statistical and machine learning models for landslide susceptibility prediction and mapping. Catena 2020, 191, 104580. [Google Scholar] [CrossRef]
  17. Cascini, L. Applicability of landslide susceptibility and hazard zoning at different scales. Eng. Geol. 2008, 102, 164–177. [Google Scholar] [CrossRef]
  18. Ingeominas. Zonificación Sismogeotécnica Indicativa del Área Metropolitana de Bucaramanga; Ingeominas, Subdirección de Amenazas Geoambientales: Bogotá, Colombia, 2001. [Google Scholar]
  19. UNDRR. Sendai framework for Disaster Risk Reduction 2015–2030. In United Nations (ONU); United Nations Office for Disaster Risk Reduction: Geneva, Switzerland, 2015. [Google Scholar]
  20. IDEAM. Geoportal—IDEAM. Consulta y Descarga de Datos Hidrometeorológicos. Available online: http://dhime.ideam.gov.co/atencionciudadano/ (accessed on 5 April 2022).
  21. Ward, D.E.; Goldsmith, R.; Jimeno, A.; Cruz, J.; Restrepo, H.; Gómez, E. Mapa geológico de Colombia, cuadrángulo h-12 Bucaramanga. Planchas 109 Rionegro—120 Bucaramanga. Cuadrángulo H-13 Pamplona. In Planchas 110 Pamplona—121Cerrito; Ingeominas: Bogotá, Colombia, 1973. [Google Scholar]
  22. Velandia, F.; Bermúdez, M. The transpressive southern termination of the Bucaramanga fault (Colombia): Insights from geological mapping, stress tensors, and fractal analysis. J. Struct. Geol. 2018, 115, 190–207. [Google Scholar] [CrossRef]
  23. Siravo, G.; Fellin, M.G.; Faccena, C.; Maden, C. Transpression and the build-up of the Cordillera: The example of the Bucaramanga fault (Eastern Cordillera, Colombia). J. Geol. Soc. 2020, 177, 14–30. [Google Scholar] [CrossRef]
  24. Osorio Naranjo, J.A.; Hernández Moreno, C.; Torres Jaimes, E.M.; Botero Santa, P.A.; Diederix, H. Modelo Geodinámico del Macizo de Santander; Ingeominas: Bogotá, Colombia, 2008. [Google Scholar]
  25. IGAC. Agustín Codazzi Geographical Institute. Open Data. Available online: https://geoportal.igac.gov.co/contenido/datos-abiertos-cartografia-y-geografia (accessed on 9 April 2022).
  26. ASF, A.P. Alaska Satellite Facility. Dataset: ASF DAAC 2015, ALOS PALSAR_Radiometric_Terrain_Corrected_Low_Res. Available online: https://asf.alaska.edu/data-sets/derived-data-sets/alos-palsar-rtc/alos-palsar-radiometric-terrain-correction/ (accessed on 9 April 2022).
  27. Servicio Geológico Colombiano (SGC). Estado de la Cartografía Geológica. Available online: http://srvags.sgc.gov.co/Flexviewer/Estado_Cartografia_Geologica/ (accessed on 9 April 2022).
  28. Instituto de Hidrología, Meteorología y Estudios Ambientales (IDEAM). Open Data Geographic Information. Find File: Land Cover. Available online: http://www.ideam.gov.co/capas-geo (accessed on 9 April 2022).
  29. Instituto de Hidrología, Meteorología y Estudios Ambientales (IDEAM). Open Data Geographic Information. Find File: Forest, Not Forest. Available online: http://www.ideam.gov.co/web/guest/galeria-de-mapas (accessed on 9 April 2022).
  30. Selby, M.J. Hillslope Materials and Processes, 2nd ed.; Oxford University Press: Oxford, UK, 1993. [Google Scholar]
  31. Kanki, M.A. Critical issues in soft rocks. J. Rock Mech. Geotech. Eng. 2014, 6, 186–195. [Google Scholar] [CrossRef]
  32. Hermelín, M. Suelos, Rocas y Formaciones Superficiales. DYNA 1985, 106, 25–29. [Google Scholar]
  33. Verstappen, H. Geomorfología Aplicada al Estudio de los Riesgos Naturales; Apuntes de Clase ITC: Holanda, The Netherlands, 1987. [Google Scholar]
  34. Carvajal, J.H. Primeras Aproximaciones a la Estandarización de la Geomorfología en Colombia; Instituto Colombiano de Geología y Minería: Bogotá, Colombia, 2008. [Google Scholar]
  35. Carvajal, J.H. Propuesta de Estandarización de la Cartografía Geomorfológica en Colombia; Servicio Geológico Colombiano: Bogotá, Colombia, 2012. [Google Scholar]
  36. Carvajal, J.H. Caracterización de la metodología geomorfológica adaptada por ingeominas. In Documento Interno Sometido a Discusión y Modificaciones; Instituto Colombiano de Geología y Minería: Bogotá, Colombia, 2002. [Google Scholar]
  37. Zinck, J.A. Geopedología: Elementos de Geomorfología para Estudios de Suelos y de Riesgos Naturales. In Enschede: International Institute for Geo-Information Science and Earth Observation (ITC); ITC Special Lecture Notes Series; University of Twente: Enschede, The Netherlands, 2012. [Google Scholar]
  38. Pike, R.J. Geomorphometry: Progress, pratice, and prospect. Z. Geomorphol. Suppl. 1995, 101, 221–238. [Google Scholar]
  39. Gallant, J.C.; Hutchinson, M.F. Digital terrain analysis. In Guidelines for Surveying Soil and Land Resources; McKenzie, N.J., Grundy, M.J., Ringrose-Voase, R.W.A.J., Eds.; Australian Soil and Land Survey Handbook Series; CSIRO: Melbourne, Australia, 2008; Volume 2. [Google Scholar]
  40. Olaya, V. Chapter 6 Basic Land-Surface Parameters. Dev. Soil Sci. 2009, 33, 141–169. [Google Scholar] [CrossRef]
  41. Wilson, J.; Gallant, J. Terrain Analysis: Principles and Applications; John Wiley & Sons, Inc.: Wiley, NY, USA, 2000. [Google Scholar]
  42. Saha, A.K.; Gupta, R.P.; Sarkar, I.; Arora, M.K.; Csaplovics, E. An approach for GIS-based statistical landslide susceptibility zonation—With a case study in the Himalayas. Landslides 2005, 2, 61–69. [Google Scholar] [CrossRef]
  43. Yalcin, A.; Reis, S.; Aydinoglu, A.; Yomralioglu, T. A GIS-based comparative study of frequency ratio, analytical hierarchy process, bivariate statistics and logistics regression methods for landslide susceptibility mapping in Trabzon, NE Turkey. Catena 2011, 85, 274–287. [Google Scholar] [CrossRef]
  44. Santacana, N. Análisis de la susceptibilidad del terreno a la formación de deslizamientos superficiales y grandes deslizamientos mediante el uso de sistemas de información geográfica. In Aplicación a la Cuenca Alta del río Llobregat; Tesi Doctoral, UPC; Departament d’Enginyeria del Terreny, Cartogràfica i Geofísica: Barcelona, Spain, 2001. [Google Scholar]
  45. Felicisimo, A.M. Modelos Digitales del Terreno; Pentalfa: Oviedo, Spain, 1994; Available online: http://www.etsimo.uniovi.es/~feli (accessed on 16 November 2016).
  46. Corominas, J.; Van Westen, C.; Frarrini, P.; Cascini, L.; Malet, J.P.; Fotopoulou, S.; Smith, J.T. Recommendations for the quantitative analysis of landslide risk. Bull. Eng. Geol. Environ. 2014, 73, 209–263. [Google Scholar] [CrossRef]
  47. Sato, H.P.; Harp, E.L. Interpretation of earthquake-induced landslides triggered by the 12 May 2008, M7. 9 Wenchuan earthquake in the Beichuan area, Sichuan Province, China using satellite imagery and Google Earth. Landslides 2009, 6, 153–159. [Google Scholar] [CrossRef]
  48. Guzzetti, F.; Mondini, A.C.; Cardinali, M.; Fiorucci, F.; Santangelo, M.; Chang, K. Landslide inventory maps: New tools for an old problem. Earth-Sci. Rev. 2012, 112, 42–66. [Google Scholar] [CrossRef]
  49. Varnes, D.J. Slope movement types and processes. In Landslides: Analysis and Control; Schuster, R.L., Krizek, R.J., Eds.; Transportation and Road Research Board, National Academy of Science: Washington, DC, USA, 1978; pp. 11–33. [Google Scholar]
  50. PMA-GCA. Proyecto Multinacional Andino (PMA): Geociencias para las Comunidades Andinas (GCA). In Movimientos en Masa en la Región Andina: Una Guía para la Evaluación de Amenazas; Publicación Geológica Multinacion, Servicio Nacional de Geología y Minería: 2007; 2007. [Google Scholar]
  51. Hungr, O.; Leroueil, S.; Picarelli, L. The Varnes classification of landslide types, an update. Landslides 2014, 11, 167–194. [Google Scholar] [CrossRef]
  52. Skemton, A.W.; Hutchinson, J.N. Stability of natural slopes and embankment foundations. In Proceedings of the Seventh lntemational Conference on Soil Mechanics and Foundation Engineering 4; Sociedad Mexicana de Mecánica de de Suelos: México City, Mexico, 1969; pp. 291–340. Available online: https://trid.trb.org/view/125702 (accessed on 7 April 2022).
  53. Bossard, M.; Feranec, J.; Otahel, J. CORINE Land Cover Technical Guide: Addendum 200; European Environment Agency: Copenhagen, Denmark, 2000; Volume 40. [Google Scholar]
  54. Charman, P.V.; Murphy, B.W. Soils: Their Properties and Management, 2nd ed.; Melbourne y Oxford; Oxford University Press: Oxford, UK, 2000. [Google Scholar]
  55. Suárez, N.J. Deslizamientos y Estabilidad de Taludes en Zonas Tropicales; Universidad Industrial de Santander: Bucaramanga, Colombia, 1998. [Google Scholar]
  56. Mugagga, F.; Kakembo, V.; Buyinza, M. Land use changes on the slopes of Mount Elgon and the implications for the occurrence of landslides. Catena 2012, 90, 39–46. [Google Scholar] [CrossRef]
  57. UN. United Nations—Climate Change. Available online: https://unfccc.int/topics/land-use/resources/redd-documents (accessed on 13 April 2022).
  58. Csiszar, I.A.; Justice, C.O.; Goldammer, J.G.; Lynham, T.; Groot, W.D.; Prins, E.M.; Stephens, G. The GOFC-GOLD fire mapping and monitoring theme: Assessment and strategic plans. In Remote Sensing and Modeling Applications to Wildland Fires; Springer: Berlin/Heidelberg, Germany, 2013; pp. 341–372. [Google Scholar] [CrossRef]
  59. Feranec, J.; Jaffrain, G.; Soukup, T.; Hazeu, G. Determining changes and flows in European landscapes 1990–2000 using CORINE land cover data. Appl. Geogr. 2010, 30, 19–35. [Google Scholar] [CrossRef]
  60. Reichenbach, P.; Busca, C.; Mondini, A.C.; Rossi, M. The Influence of Land Use Change on Landslide Susceptibility Zonation: The Briga Catchment Test Site (Messina, Italy). Environ. Manag. 2014, 54, 1372–1384. [Google Scholar] [CrossRef]
  61. Van Beek, L.H.; Van Asch, T.J. Regional assessment of the effects of land-use change and landslide hazard by means of physically based modeling. Nat. Hazards 2004, 30, 289–304. [Google Scholar] [CrossRef]
  62. Pisano, L.; Zumpano, V.; Malek, Ž.; Rosskopf, C.M.; Parise, M. Variations in the susceptibility to landslides, as a consequence of land cover changes: A look to the past, and another towards the future. Sci. Total Environ. 2017, 601, 1147–1159. [Google Scholar] [CrossRef]
  63. Guzzetti, F.; Carrara, A.; Cardinali, M.; Reichenbach, P. Landslide hazard evaluation: A review of current techniques and their application in a multi-scale study, Central Italy. Geomorphology 1999, 31, 181–216. [Google Scholar] [CrossRef]
  64. Gardiner, C.W. Handbook of Stochastic Methods; Springer: Berlin, Germany, 1985; Volume 3. [Google Scholar]
  65. Barrera, D. Modelos Deterministicos y Probabilisticos; Universidad José María Vargas: Caracas, Venezuela, 2016. [Google Scholar]
  66. Simpson, P.K. Artificial Neural System—Foundation, Paradigm, Application and Implementation; Pergamon Press: New York, NY, USA, 1990. [Google Scholar]
  67. Anderson, D.; McNeill, G. Artificial Neural Networks Technology; Kaman Sciences Corporation: Bloomfield, CT, USA, 1992. [Google Scholar]
  68. De Brio, B.; Sanz, A. Redes Neuronales y Srosos; Alfaomega: México, México, 2007. [Google Scholar]
  69. Pradhan, B.; Lee, S.; Buchoithner, M.F. A GIS-based back-propagation neural network model and its cross-application and validation for landslide susceptibility analyses. Comput. Environ. Urban Syst. 2010, 34, 216–235. [Google Scholar] [CrossRef]
  70. Lee, S.; Park, I.; Choi, J.K. Spatial prediction of ground subsidence susceptibility using an artificial neural network. Environ. Manag. 2012, 49, 347–358. [Google Scholar] [CrossRef]
  71. Valencia Ortiz, J.A.; Martínez-Graña, A.M. A neural network model applied to landslide susceptibility analysis (Capitanejo, Colombia). Geomat. Nat. Hazards Risk 2018, 9, 1106–1128. [Google Scholar] [CrossRef]
  72. Bonham-Carter, G.F.; Agterberg, F.P.; Wright, D.F. Integration of geological datasets for gold exploration in Nova Scotia. Photogramm. Eng. Remote Sens. 1988, 54, 1585–1592. [Google Scholar] [CrossRef]
  73. Bonham, G. Geographic Information Systems for Geoscientists: Modelling with GIS (No. 13); Pergamon, Ed.; Elsevier: Amsterdam, The Netherlands, 1994; Volume 13. [Google Scholar]
  74. van Westen, C. Guidelines for the Generation of 1:50.000 Scale Landslide Inventory, Susceptibility Maps, and Qualitative Risk Maps, Illustrated with Case Studies of the Provinces Thanh Hoa and Nghe An; University of Twente: Enschede, The Netherlands, 2013. [Google Scholar]
  75. Chen, W.; Chai, H.; Sun, X.; Wang, Q.; Ding, X.; Hong, H. A GIS-based comparative study of frequency ratio, statistical index and weights-of-evidence models in landslide susceptibility mapping. Arab. J. Geosci. 2016, 9, 204. [Google Scholar] [CrossRef]
  76. Kendall, M.G. Multivariate Analysis; Griffin: London, UK, 1975. [Google Scholar]
  77. Carrara, A.; Cardinali, M.; Detti, R.; Guzzetti, F.; Pasqui, V.; Reichenbach, P. GIS techniques and statistical models in evaluating landslide hazard. Earth Surf. Proc. Land 1991, 16, 427–445. [Google Scholar] [CrossRef]
  78. Bernknopf, R.; Cambell, R.; Brookshire, D.; Shapiro, C. A probabilistic approach to landslide hazard mapping in Cincinnati, Ohio, with applications for economic evaluation. Bull Int. Assoc. Eng. Geol. 1988, 25, 39–56. [Google Scholar] [CrossRef]
  79. Fiuza, P.D.; Rodríguez, P.J. La regresión logística: Una herramienta versátil. Nefrologia 2000, 20, 495–500. [Google Scholar]
  80. Baeza, C.; Corominas, J. Assessment of shallow landslide susceptibility by means of multivariate statistical techniques. Earth Surface Processes and Landforms. J. Br. Geomorphol. Res. Group 2001, 26, 1251–1263. [Google Scholar] [CrossRef]
  81. Kavzoglu, T.; Kutlug, S.E.; Colkesen, I. An assessment of multivariate and bivariate approaches in landslide susceptibility mapping: A case study of Duzkoy district. Nat. Hazards 2015, 76, 471–496. [Google Scholar] [CrossRef]
  82. Gómez, N.F.; Mendez, L.M.; Valencia Ortiz, J.A. Zonificación de la susceptibilidad por movimientos en masa, escala 1:10.000 para el municipio de Suratá, Santander, aplicando métodos estadísticos multivariables (regresión logística). In Trabajo de Grado; Universidad Industrial de Santander: Bucaramanga, Colombia, 2017. [Google Scholar]
  83. Lee, S.; Pradhan, B. Landslide hazard mapping at Selangor, Malaysia using frecuency ratio and logistic regression models. Landslides 2007, 4, 33–41. [Google Scholar] [CrossRef]
  84. Aguayo, M. Como Hacer una Regresión Logística con SPSS “Paso a Paso” (I). Available online: http://metodos-avanzados.sociales.uba.ar/wp-content/uploads/sites/216/2014/03/Regres_log_AGUAYO-otros.pdf (accessed on 16 April 2022).
  85. Aguayo, M.; Lora, E. Como Hacer una Regresión Logística Binaria “Paso a Paso” (II): Análisis Multivariante. Available online: https://documentop.com/como-hacer-una-regresion-logistica-binaria-fabis_5a0eb71e1723dd59aa821c5b.html (accessed on 16 April 2022).
  86. Dahal, R.K.; Hasegawa, S.; Nonomura, A.; Yamanaka, M.; Dhakal, S.; Paudyal, P. Predictive modelling of rainfall-induced landslide hazard in the Lesser Himalaya of Nepal based on weights-of-evidence. Geomorphology 2008, 102, 496–510. [Google Scholar] [CrossRef]
  87. Vakhshoori, V.; Zare, M. Landslide susceptibility mapping by comparing weight of evidence, fuzzy logic, and frequency ratio methods. Geomat. Nat. Hazards Risk 2016, 7, 1731–1752. [Google Scholar] [CrossRef]
  88. Salitchev, K.A. Cartografía; Pueblo y Educación: La Habana, Cuba, 1979. [Google Scholar]
  89. Verstappen, H.T.; van Zuidam, R.A. The ITC System of Geomorphological Survey; ITC Publication Number 10; ITC: Enschede, The Netherlands, 1992. [Google Scholar]
  90. van Zuidam, R.A. Aerial Photointerpretation in Terrain Analysis and Geomorphological Mapping; Smiths Publishers: Tha Hague, The Netherlands, 1986. [Google Scholar]
  91. Cruden, D.M.; Varnes, D.J. Landslide Types and Processes. Special Report. In Lanslides: Investigation and Mitigation; Turner, A.K., Schuster, R.L., Eds.; National Research Council 247; National Academy of Sciences: Washington, DC, USA, 1996; pp. 36–75. [Google Scholar]
  92. IDEAM. Monitoreo de la Superficie Cubierta por Bosque Natural. Available online: http://www.ideam.gov.co/web/ecosistemas/superficie-cubierta-por-bosque-natural?inheritRedirect=true (accessed on 13 April 2022).
  93. Li, Y.; Mo, P. A unified landslide classification system for loess slopes: A critical review. Geomorphology 2019, 340, 67–83. [Google Scholar] [CrossRef]
  94. Valencia Ortiz, J.A.; Martínez-Graña, A.M. Morphometric Evaluation and Its Incidence in the Mass Movements Present in the Chicamocha Canyon, Colombia. Sustainability 2023, 15, 1140. [Google Scholar] [CrossRef]
  95. Bathrellos, G.D.; Gaki-Papanastassiou, K.; Skilodimou, H.D.; Papanastassiou, D.; Chousianitis, K.G. Potential suitability for urban planning and industry development using natural hazard maps and geological–geomorphological parameters. Environ. Earth Sci. 2012, 66, 537–548. [Google Scholar] [CrossRef]
  96. Valencia Ortiz, J.A.; Martínez-Graña, A.M. Calculation of precipitation and seismicity thresholds as triggers for mass movements in the region of Bucaramanga, Colombia. Ecol. Indic. 2023, 152, 110355. [Google Scholar] [CrossRef]
Figure 1. Location map of the study area.
Figure 1. Location map of the study area.
Remotesensing 15 04567 g001
Figure 2. (a) Map of the slope. (b) Orientation map.
Figure 2. (a) Map of the slope. (b) Orientation map.
Remotesensing 15 04567 g002
Figure 3. Geological map of the study area. Modified from the geological map at scale 1:100,000, Colombian Geological Service, sheet 109 [21].
Figure 3. Geological map of the study area. Modified from the geological map at scale 1:100,000, Colombian Geological Service, sheet 109 [21].
Remotesensing 15 04567 g003
Figure 4. Methodological scheme for the calculation of susceptibility by mass movements. 100 K: scale 1:100,000, SGC: Colombian geological service, DEM: digital elevation model, UGS: surface geological unit, 25 K: scale 1:25,000, MM: mass movement, LR: logistic regression.
Figure 4. Methodological scheme for the calculation of susceptibility by mass movements. 100 K: scale 1:100,000, SGC: Colombian geological service, DEM: digital elevation model, UGS: surface geological unit, 25 K: scale 1:25,000, MM: mass movement, LR: logistic regression.
Remotesensing 15 04567 g004
Figure 6. Simulation results for the 6-50-1 architecture. (a) Quadratic regression, (b) training status, (c) performance curve, and (d) confounding matrix.
Figure 6. Simulation results for the 6-50-1 architecture. (a) Quadratic regression, (b) training status, (c) performance curve, and (d) confounding matrix.
Remotesensing 15 04567 g006
Figure 7. Susceptibility map by mass movements generated by the artificial neural network algorithm with a 6-50-1 architecture and classified into three categories (high, medium, and low).
Figure 7. Susceptibility map by mass movements generated by the artificial neural network algorithm with a 6-50-1 architecture and classified into three categories (high, medium, and low).
Remotesensing 15 04567 g007
Figure 8. Susceptibility map by mass movements generated by bivariate statistical analysis and classified into three categories (high, medium, and low).
Figure 8. Susceptibility map by mass movements generated by bivariate statistical analysis and classified into three categories (high, medium, and low).
Remotesensing 15 04567 g008
Figure 9. Susceptibility map for mass movements generated by the logistic regression algorithm (LR) and classified into three categories (high, medium, and low).
Figure 9. Susceptibility map for mass movements generated by the logistic regression algorithm (LR) and classified into three categories (high, medium, and low).
Remotesensing 15 04567 g009
Figure 10. ROC curve and its evaluation of the AUC for each proposed method and confusion matrix generated on the final result of each susceptibility model in its three categories.
Figure 10. ROC curve and its evaluation of the AUC for each proposed method and confusion matrix generated on the final result of each susceptibility model in its three categories.
Remotesensing 15 04567 g010
Figure 11. Susceptibility models generated by the bivariate statistical analysis method for each basin.
Figure 11. Susceptibility models generated by the bivariate statistical analysis method for each basin.
Remotesensing 15 04567 g011
Table 1. Classification of slopes and areas without and affected by mass movements.
Table 1. Classification of slopes and areas without and affected by mass movements.
Ranking Classification
CountArea km2MinMaxRangeMeanStd
10,9821.71574.897°74.897°31.42936°9.7631
Slope
RangeWithout MMWith MM
CountArea km2%CountArea km2%
0–10°167,86126.236.522010.031.83
10–20°595,62893.0723.1210720.179.76
20–30°937,852146.5436.4034570.5431.48
30–40°662,061103.4525.7042890.6739.05
40–50°183,42328.667.1216320.2614.86
50–60°25,3793.970.992910.052.65
60–70°38650.600.15390.010.36
>70°2360.040.0110.000.01
Table 2. Correlation of the susceptibility model by the ANN method against the MM inventory for each basin. Evaluation performed by number of pixels.
Table 2. Correlation of the susceptibility model by the ANN method against the MM inventory for each basin. Evaluation performed by number of pixels.
BasinSusceptibility by MMQuantity of MM per Category
Low%Medium%High%Low%Medium%High%
Tona383,50030.5547,71143.5326,80626.02176.7130340.1172653.2
Río Frio 90703.1114,28439.1169,13857.8342.069742.093056.0
Rio del Hato19,81210.0108,57754.670,58135.54210.218645.018544.8
Rio de Oro77,55215.6210,29442.3209,12742.1374.235940.648955.3
W/B56,14616.6165,16548.9116,59034.52515.2142529.8310564.9
W/B: Region that does not have a basin expression.
Table 3. Correlation of the susceptibility model by the Bivariate method against the MM inventory for each basin. Evaluation performed by number of pixels.
Table 3. Correlation of the susceptibility model by the Bivariate method against the MM inventory for each basin. Evaluation performed by number of pixels.
BasinSusceptibility by MMQuantity of MM per Category
Low%Medium%High%Low%Medium%High%
Tona637,55850.8441,65535.2174,78513.942313.0113635.0168451.9
Río Frio 30631.0172,97659.1116,43339.800.065639.5100560.5
Rio del Hato18570.9119,90160.377,17538.841.018945.822053.3
Rio de Oro196,39139.6190,92938.5108,12221.8738.235640.245651.5
W/B69572.1128,35038.3199,66459.600.054411.4423588.6
W/B: Region that does not have a basin expression.
Table 4. Spearman correlation of the variables taken for the mass movement susceptibility model with the logistic regression (LR) method.
Table 4. Spearman correlation of the variables taken for the mass movement susceptibility model with the logistic regression (LR) method.
CorrelationMMLANDTEMPENRUGUGS
MMSpearman correlation1.000
Sig. (bilateral)
N11,994
LANDSpearman correlation−0.391 **1.000
Sig. (bilateral)0.000
N11,99411,994
TEMSpearman correlation−0.289 **−0.046 **1.000
Sig. (bilateral)0.0000.000
N11,99411,99411,994
PENSpearman correlation0.056 **−0.198 **−0.316 **1.000
Sig. (bilateral)0.0000.0000.000
N11,99411,99411,99411,994
RUGSpearman correlation0.107 **−0.057 **−0.176 **0.141 **1.000
Sig. (bilateral)0.0000.0000.0000.000
N11,99411,99411,99411,99411,994
UGSSpearman correlation−0.108 **0.252 **0.041 **−0.383 **−0.150 **1.000
Sig. (bilateral)0.0000.0000.0000.0000.000
N11,99411,99411,99411,99411,99411,994
**. The correlation is significant at the 0.01 level (bilateral).
Table 5. Correlation of the susceptibility model by the LR method against the MM inventory for each basin. Evaluation performed by number of pixels.
Table 5. Correlation of the susceptibility model by the LR method against the MM inventory for each basin. Evaluation performed by number of pixels.
BasinSusceptibility by MMQuantity of MM per Category
Low%Medium%High%Low%Medium%High%
Tona359,62928.7479,23038.2415,13733.187927.1124038.2112334.6
Río Frio 85,43729.2135,11646.271,91824.659335.764438.842325.5
Rio del Hato51,10125.786,77443.661,05630.710826.313131.917241.8
Rio de Oro76,19315.4246,07849.7173,17035.018020.440946.329433.3
W/B70,55921.1141,08242.1123,32836.898920.7261254.7117724.6
W/B: Region that does not have a basin expression.
Table 6. Variation of the susceptibility calculation by MM correlated with the results of the total area (regional) and by basins.
Table 6. Variation of the susceptibility calculation by MM correlated with the results of the total area (regional) and by basins.
Study AreaBasin
SusceptibilityArea (km2)%Area (km2)%Variation (%)
Low99.6250.8484.9243.34−14.75
Tona RiverMedium69.0135.2274.2137.877.53
High27.3113.9436.8118.7934.77
Total195.94 195.94
Low0.481.053.838.37699.51
Río Frio RiverMedium27.0359.1421.7647.61−19.50
High18.1939.8120.1244.0210.57
Total45.70 45.70
Low0.290.937.4023.792448.73
Río del Hato RiverMedium18.7360.2712.9041.49−31.16
High12.0638.7910.7934.72−10.51
Total31.08 31.08
Low30.6939.6434.5744.6612.67
Río de Oro RiverMedium29.8338.5424.8132.05−16.85
High16.8921.8218.0323.296.74
Total77.41 77.41
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Valencia Ortiz, J.A.; Martínez-Graña, A.M.; Méndez, L.M. Evaluation of Susceptibility by Mass Movements through Stochastic and Statistical Methods for a Region of Bucaramanga, Colombia. Remote Sens. 2023, 15, 4567. https://doi.org/10.3390/rs15184567

AMA Style

Valencia Ortiz JA, Martínez-Graña AM, Méndez LM. Evaluation of Susceptibility by Mass Movements through Stochastic and Statistical Methods for a Region of Bucaramanga, Colombia. Remote Sensing. 2023; 15(18):4567. https://doi.org/10.3390/rs15184567

Chicago/Turabian Style

Valencia Ortiz, Joaquín Andrés, Antonio Miguel Martínez-Graña, and Lenny Mejía Méndez. 2023. "Evaluation of Susceptibility by Mass Movements through Stochastic and Statistical Methods for a Region of Bucaramanga, Colombia" Remote Sensing 15, no. 18: 4567. https://doi.org/10.3390/rs15184567

APA Style

Valencia Ortiz, J. A., Martínez-Graña, A. M., & Méndez, L. M. (2023). Evaluation of Susceptibility by Mass Movements through Stochastic and Statistical Methods for a Region of Bucaramanga, Colombia. Remote Sensing, 15(18), 4567. https://doi.org/10.3390/rs15184567

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop