Landslide Risk Assessments through Multicriteria Analysis

Abstract

Natural risks comprise a whole range of disasters and dangers, requiring comprehensive management through advanced assessment, forecasting, and warning systems. Our specific focus is on landslides in difficult terrains. The evaluation of landslide risks employs sophisticated multicriteria models, such as the weighted sum GIS approach, which integrates qualitative parameters. Despite the challenges posed by the rugged terrain in Northern Algeria, it is paradoxically home to a dense population attracted by valuable hydro-agricultural resources. The goal of our research is to study landslide risks in these areas, particularly in the Mila region, with the aim of constructing a mathematical model that integrates both hazard and vulnerability considerations. This complex process identifies threats and their determining factors, including geomorphology and socio-economic conditions. We developed two algorithms, the analytic hierarchy process (AHP) and the fuzzy analytic hierarchy process (FAHP), to prioritize criteria and sub-criteria by assigning weights to them, aiming to find the optimal solution. By integrating multi-source data, including satellite images and in situ measurements, into a GIS and applying the two algorithms, we successfully generated landslide susceptibility maps. The FAHP method demonstrated a higher capacity to manage uncertainty and specialist assessment errors. Finally, a comparison between the developed risk map and the observed risk inventory map revealed a strong correlation between the thematic datasets.

1. Introduction

A landslide is a mass wasting process that occurs on both natural and engineered slopes, involving the movement of rock, debris, or earth under the influence of gravity [1]. Factors such as human activities, climate, and geomorphology contribute to landslides, posing a global threat to human settlements and infrastructure. These events can devastate buildings, deform arable land, and cause significant natural changes. Nearly 9% of natural disasters worldwide are attributed to landslides [2]. With the global population growing and communities expanding, vulnerability to landslide disasters is increasing, increasing the risk of catastrophic losses. Slope instability, driven by various factors, poses a significant threat to people and their environments. To address landslides, we have referred to previous research published by leading authors in the field [3,4,5].

The Mila region, characterized by its dense population and multiple water dams, experienced the reactivation of the El Kherba landslide following the Mila earthquake of 7 August 2020. Recent papers [6,7] detail this event, describing a 26 million cubic meter landslide triggered by an exceptionally small earthquake. This recent earthquake has heightened concerns among administrators regarding the condition of the dams, prompting them to request a comprehensive qualitative and quantitative assessment from scientists. The University Center of Mila, in partnership with our department, has initiated a study on this phenomenon in the region, organizing two international conferences to address this issue. Our main focus is on evaluating the landslide risk in the Mila region, starting with an in-depth bibliographic study on the topic.

Researchers have proposed various approaches to understanding landslides, examining factors such as climate change [7], earthquake effects [8], dam construction [9], rapid urbanization, deforestation [10], and other anthropogenic factors. Examples from studies in the Haraz watershed and Brazil have shown the use of multicriteria analysis, such as the analytic hierarchy process (AHP), for mapping landslide susceptibility. Additionally, in Portugal, the AHP has been applied to study landslides in a coastal area from Figueira da Foz to Nazaré. These studies highlight the importance of multicriteria analysis and emphasize the use of remote sensing tools and specialized software for image processing and spatial data collection in assessing landslides globally [11,12,13,14].

Remote sensing and geographic information systems (GIS) are essential tools for quantifying land cover and land use (LCLU) changes on Earth [15,16,17,18] and for assessing landslide risk. Remote sensing applications use high-resolution satellite imagery to identify and monitor landslide-prone areas by mapping critical terrain features, such as slope and vegetation cover. Herein, our aim is to evaluate the cumulative factors contributing to this vulnerability, which are often difficult to change and can significantly affect exposure to natural hazards. Researchers are actively developing processes, methodologies, and techniques to effectively assess and manage risks associated with natural hazards. Mathematical techniques like fuzzy theory and multicriteria analysis, as explored in various publications [11,12,13,19,20,21], form the foundation of our study. We have combined the analytic hierarchy process (AHP) and fuzzy analytic hierarchy process (FAHP) approaches to use operational research tools, specifically multicriteria analysis. These methods have assessed the dangers connected with landslides in the Mila province. Our main objective has been to present a thorough strategy that makes use of both fuzzy and deterministic frameworks, resulting in a solid natural risk assessment model that has been rigorously tested.

The present research is divided into several sections as follows: Section 1 provides an overview of the research problem, covering various landslide types, causes, real-world instances, and methods for assessing landslide risk. Section 2 introduces essential tools, such as remote sensing (e.g., ArcGIS, ENVI), for spatial data collection, exploring operations research principles with a focus on multicriteria optimization. Section 3 delves into the AHP and FAHP methods, outlining matrices, implementation steps, and criteria for evaluating landslide hazard. Section 4 presents and discusses multicriteria analysis outcomes, highlighting high-risk landslide areas with preventive maps tailored to the Mila province and investigating preventive strategies for each high-risk zone. The paper concludes by summarizing the findings and providing future perspectives [14,15,16,17,18,22,23,24].

1.1. Study Area

The province of Mila is located in the northeast of Algeria, 50 km northwest of Constantine (Figure 1). It is located at the center of the Neogene basin of Mila, which is bordered to the north by the mountain ridges of M’cid Aicha and Sidi Driss, to the south by Djebel Ossmane and Djebel Grouz, to the east by the massifs of Djebel Chettaba, Akhal, and Kheneg, and to the west by Djebel Boucharef and Oukissène [25,26]. The Mila basin contains geological formations of clay and marl with detrital characteristics, thickening between 40 to 50 m, and interspersed with thin sandy layers varying from 10 to 20 cm in thickness. These clay-marl formations are vulnerable to soil shrinkage and swelling, leading to differential settlements that result in structural disorders, impacting the road network and buildings in the region [27].

1.2. Overview about Landslides

Landslides are the movement of soil on natural or artificial slopes. They take different forms depending on their shape and movement characteristics, including rotational, planar, and translational landslides. In this context, the term “mass movements of materials” includes soils, rocks, and fractured rocks [27]. The occurrence of landslides depends on factors such as slope, geology, vegetation, precipitation, and human activity. Different types of landslides are categorized according to their movement and geometry. The main types are rotational (with a curved shape), planar (with a flat surface), and translational (with a curved or wedge-shaped surface) [28].

Landslides can be caused by various factors, which can be divided into two main groups: natural causes, which are caused by geological, morphological, and climatic factors such as soil composition, topography, geological phenomena like earthquakes, and adverse climatic conditions; and anthropogenic causes, arising from human activities like mining, construction on unstable terrain, excavation at the base of slopes (e.g., Chaancun landslide in Dalian, China), and deforestation, as has been observed in the Nilgiris region in India, leading to increased landslide occurrences [11,12,20,21].

1.3. Landslides Worldwide

Landslides have caused destruction worldwide. The 1963 Vajont landslide in Italy, one of the deadliest in history, resulted in the loss of over 2000 lives [29]. Similarly, the 2014 Oso landslide in the United States led to 43 deaths, with heavy rainfall contributing to soil saturation and extensive devastation. In 2010, China was struck by a series of natural disasters, including severe flooding and landslides in Zhouqu region, which resulted in over 1200 fatalities and extensive damage to homes and infrastructure. El Salvador witnessed a catastrophic landslide triggered by an earthquake on the Pan-American Highway in 2001, resulting in around 450 deaths and 1200 missing people.

Algeria has experienced a considerable number of landslides, with the northern regions exhibiting a heightened susceptibility, largely attributed to their complex geomorphology and geological structures. Notable incidents include landslides in Bejaia and Tlemcen, which have caused significant damage and disruptions to urban regions and economic routes. The province of Mila is located in the north of Algeria, a region that is susceptible to seismic activity. A certain number of tremors occur there on a regular basis, although most of them are not felt by the population. On Friday, 7 August 2020, at 07:15, the Center for Research in Astronomy, Astrophysics, and Geophysics (CRAAG) recorded an earthquake with a magnitude of 4.9 degrees on the Richter scale in Hamala, within the province of Mila. This seismic event resulted in landslides in the Mila region.

Type of Approaches

The assessment of landslide hazard employs two principal methodologies: geomorphological field analysis and geographic information systems (GIS)-based parameter mapping. Geomorphological field analysis relies on expert on-site observations, offering a swift hazard evaluation across various scales, but necessitating extensive fieldwork. In contrast, GIS-based mapping integrates geological, geomorphological, hydrogeological, and other factors with specific weight assignments, expediting hazard assessment while introducing potential subjectivity concerns due to parameter weighting.

Quantitative methods are employed to analyze historical events in landslide-prone areas through bivariate and multivariate statistical analyses [30]. The former calculates density by merging parameters and landslide maps, a technique that is favored by geologists, but which is complex due to variable correlations. The latter involves the categorization of terrain into zones using multiple variables, employing map overlays and GIS. It may yield illogical results with extensive datasets.

The geotechnical approach assesses slope failure using deterministic and probabilistic methods. The former calculates a safety factor $\left(Fs\right)$ through slope stability analysis, marking instability if Fs < 1. The latter factors in variability for precise stability estimates. The modeling approach visually represents landslides through empirical, analytical, or numerical models, but may deviate from reality due to environmental homogeneity [31,32].

2. Preliminary Tools

This section presents the essential tools necessary for our study on landslide hazard assessment. The following section will discuss the utilization of remote sensing and specialized softwares (ENVI 5.3 and ArcGIS 10.2) employed for image processing (Landsat) and spatial data collection. Furthermore, we will introduce operations research concepts, with a particular emphasis on multicriteria optimization.

2.1. Image Processing

This investigation explores the field of remote sensing, which is a contactless method for capturing object characteristics using instruments like cameras, lasers, radars, or sonars, often deployed on remote platforms like satellites. The application of digital processing techniques, including multispectral classification, enables the generation of thematic maps portraying surface attributes. Remote sensing allows for the remote acquisition of environmental data using instruments such as aircraft, spacecraft, or satellites, providing valuable insights into our surroundings. Satellite images, derived from sensors aboard satellites, capture specific scenes with multiple bands, each linked to a distinct wavelength. These bands produce grayscale images with a pixel resolution of a few meters, and leveraging pixel spectral signatures ensures efficient data capture. Satellite images come in three types based on band count and spectral coverage: panchromatic, multispectral, and hyper-spectral. Multi-band sensors are crucial for collecting landscape data, especially in rural and peri-urban regions like agricultural lands, forests, mountains, and coastal areas, offering insights into cultivated fields and, at times, human habitation.

It is noteworthy that vegetation emits distinctive chlorophyll-related signals, particularly in the infrared and near-infrared spectra, emphasizing the importance of a broader spectral range for gathering relevant information. Furthermore, this study addresses image segmentation, a technique whereby images are divided into pixel groups based on similarity, creating connected components representing regions. In satellite image segmentation, methods use image statistics, either supervised with prior information or unsupervised with approaches like watershed transformation. Parametric techniques such as k-means, fuzzy k-means, and density-based clustering, while effective, require an estimated initial cluster count. Mean-shift clustering offers promise without such estimates. High-resolution satellite images, while informative, may overwhelm standard computers due to pixel count. Object-based image analysis overcomes this challenge by working with object groups rather than every pixel [11,12,13,14,15,16,17,18,33].

2.2. Multi-Criteria Decision-Making (MCDM)

Criteria in Decision-Making

A criterion, designated by function g over action set A with ordered values, reflects user preferences [34]. Criteria aid decision-makers by summarizing action evaluations along shared dimensions, essentially representing a “point of view”. Criteria are real-valued functions on potential actions, enabling comparisons as f_b ≥ f_a ⇒ bS_fa, which signifies that b is at least as good as a, under criterion f.

In decision support, criteria relate to the decision-maker’s value scale. Subjectivity in criterion selection arises from the decision-maker’s attention to specific aspects. Multi-criteria problems pose challenges due to complexity, strategic implications, conflicting criteria, uncertainty, and limited information. They lack a single objective solution, rendering the notion of an optimal solution ambiguous. Solutions depend on the decision-maker’s personality, context, problem formulation, and chosen decision support method, acknowledging subjectivity.

2.3. Multicriteria Analysis Methods

This section will examine two main approaches: the analytic hierarchy process (AHP) and the fuzzy analytic hierarchy process (FAHP). These methods have been widely applied in the landslide evaluation literature, and are therefore well suited to our study.

In multicriteria decision-making, the “partial aggregation” or “synthetic dominance” approach entails a certain degree of compromise with regard to result clarity, accommodating non-quantitative criteria and diverse units through pairwise comparisons. This approach is frequently used in a variety of decision-making contexts [35,36,37,38,39,40,41].

AHP structures and prioritizes criteria and sub-criteria using a comparison matrix, integrating stakeholder expertise for weight assessment. Matrix calculations determine criteria priorities, aiding quantitative evaluation for decision-making. Conversely, FAHP, an extension of AHP, considers uncertainty and subjectivity using fuzzy sets and membership functions, accommodating vague or imprecise information and generating flexible results [11]. Additionally, different multicriteria aggregation methods exist, including one that includes all performances in an aggregation function, assuming commensurability and transitivity of judgments. However, this may oversimplify nuances [36]. Another method is TOPSIS, aiming to select the best alternative by evaluating the distance to the ideal and negative ideal alternatives [42]. The SMART method uses an additive form for aggregating evaluations based on cardinal scales and articulated preferences.

2.3.1. The Analytic Hierarchy Process (AHP) Technique

The analytic hierarchy process (AHP), developed by Thomas Saaty in the 1970s, facilitates decision-making in complex scenarios with numerous criteria that are difficult to prioritize. AHP addresses these challenges by organizing problems into a hierarchical structure and using subjective stakeholder judgments in comparison matrices to determine the relative importance of criteria and alternatives. This method ultimately identifies the optimal solution based on predefined criteria. The key principles of AHP, as detailed by Saaty in his 1986 work, include decomposition, comparative judgments, and priority synthesis. AHP constructs hierarchies, establishes priorities, and performs logical coherence checks to effectively organize and evaluate the components of a problem, beginning with a clear hierarchical structure of criteria, sub-criteria, and alternatives.

Establishment of Priorities

The process starts by conducting pairwise comparisons at the same hierarchical level to determine the priorities of elements. These comparisons, based on logical and experiential evaluations, use a matrix framework to establish concrete judgments. Numerical values are assigned to the pairwise comparison matrix, representing the relative importance of one element over another for a specific property. The matrix is then normalized, and eigenvectors are calculated by averaging the column elements, yielding accurate normalized eigenvector values that represent criterion weights (see Algorithm 1).

Evaluating Consistency

AHP uses a consistency ratio to assess the coherence of judgments, with a ratio exceeding 0.1 suggesting the need for revisions. The consistency index (CI) quantifies this and is calculated using the eigenvalue method.

$$CI={\displaystyle \frac{{\lambda }_{max}-n}{\left(n-1\right)}}$$

(1)

where CI represents the consistency index, λ_max is the principal eigenvalue of the matrix, and n is the number of criteria or alternatives. The consistency ratio (CR) is computed as:

$$CR={\displaystyle \frac{CI}{RCI}}$$

(2)

CR signifies the consistency ratio, CI denotes the consistency index, and RCI stands for the random consistency indices. A technical description is given in the following subsection.

Algorithm 1: AHP Method

2.3.2. Fuzzy Analytical Hierarchy Process (FAHP)

The fuzzy analytical hierarchy process (FAHP) enhances the traditional AHP by incorporating fuzzy logic, thus refining the multi-criteria decision-making framework. Unlike the direct ranking system of AHP, FAHP utilizes linguistic variables and triangular fuzzy numbers for a more detailed and nuanced comparison process, as outlined in [40,41,42]. This approach, rooted in the pioneering works of van Laarhoven and Pedrycz [43], and further developed by Buckley [44] and Chang [45], places a strong emphasis on fuzzy priorities and advanced comparison methods. In this study, we adopt Buckley’s methodology for assigning relative weights to the criteria and alternatives (see Algorithm 2). The pairwise contribution matrix, utilizing fuzzy triangular numbers as described in Equation (3), encodes the preferences of decision-makers. Within this matrix, ${\tilde{d}}_k$ represents the preference of the k-th decision-maker for the i-th criterion over the j-th criterion, utilizing fuzzy triangular numbers for a more precise expression. The structure of this matrix is carefully detailed below:

$$\left(\begin{array}{ccc}{\overline{\mathrm{D}}}_{11}{\overline{\mathrm{d}}}_{12}& \cdots & {\overline{\mathrm{d}}}_{1n}\\ {\overline{\mathrm{d}}}_{121}{\overline{\mathrm{d}}}_{22}& \dots & {\overline{\mathrm{d}}}_{2n}\\ ⋮& \ddots & ⋮\\ {\overline{\mathrm{d}}}_{n1}{\overline{\mathrm{d}}}_{n2}& \cdots & {\overline{\mathrm{d}}}_{1n}\end{array}\right)$$

(3)

The steps to implement FAHP for calculating the criteria’s normalized weights are systematically outlined in the technical description below.

Algorithm 2: Algorithm for Finding Normalized Defuzzified Weights

3. Landslide Hazard Assessment Criteria and Methodology

This section outlines criteria for landslide hazard assessment, derived from remote sensing data and expert collaboration. In our study, we utilized satellite imagery from Landsat 8, which serves as a crucial data source for remote sensing analysis, distinct from GIS software tools such as ENVI and ArcGIS, which are used for processing and analyzing the imagery. We collected data on NDMI, slope aspect, land use, rainfall, NDVI, distance to hydrographic networks, network density, elevation, slope, and lithology. Analyzing these criteria enhances our understanding of landslide risk factors for informed risk management.

Heavy rain saturates soil, heightening water pressure, and triggering landslides. Studies have consistently linked rainfall to landslides [1]. Snowmelt also contributes, releasing water as temperatures rise, causing ground movements [46]. Altitude influences slope stability, correlating with elevation. Higher elevations mean lower temperatures and precipitation, impacting vegetation. Varying altitudinal zones exhibit landslide susceptibilities, influenced by shear strength and colluvium thickness. Proximity to water bodies shapes landslides, with drainage networks molding unstable regions. Dai and Lee’s research (2002) [47] underscores drainage’s profound effect [48,49,50]. Saturated materials impact slope stability, influenced by nearby drainage. River flows contribute to surface erosion and water level elevation. Groundwater influences flow patterns, vegetation, shear strength, and stability. During storm events, groundwater–river connections heighten landslide risk.

In essence, the proximity of water bodies and the interaction between surface and groundwater play a pivotal role in slope stability, impacting landslides. This understanding is crucial for effective risk assessment and proactive management. Slope orientation, or aspect, is vital in landslide susceptibility. Measured clockwise, aspect-related features influence moisture retention, vegetation, and soil strength. Certain orientations enhance precipitation concentration, affecting landslide initiation. The chosen methodology employs the analytic hierarchy process (AHP), integrating deterministic and fuzzy dimensions to handle uncertainty in landslide hazard assessment.

3.1. Hierarchical Decomposition

In this stage, we created a hierarchy of criteria, sub-criteria, and alternatives (Figure 2). Level 0 focuses on evaluating landslide hazard. The decision criteria included NDMI, aspect, land use, rainfall, NDVI, buffer, network density, altitude, slope, and lithology. Alternatives encompass very low, low, moderate, high, and very high.

3.2. Comparative Evaluation

At this stage, we employed AHP (deterministic) or FAHP (fuzzy) to compare elements within the same hierarchical level. The objective was to assess the relative importance of each element in relation to those at higher levels. We used appropriate scales, which were quantitative, qualitative, or expert-based depending on the study context, to evaluate the preference or significance between element pairs. To ensure consistency, the comparison matrices used adhered to fundamental principles. This step was crucial for constructing a well-organized hierarchy and making informed decisions, whether using the deterministic AHP or the fuzzy FAHP approach. Supplementary illustrations can be included if necessary.

3.2.1. Priority Determination

In this phase, we determined priorities by evaluating the relative significance of each element in the hierarchy, utilizing the assessments conducted in Step 2. Using the eigenvector method, priority values for each element were calculated, and the resulting normalized vector values provided a clear indication of the importance of criteria, where larger values denoted greater significance.

3.2.2. Priority Synthesis

Once priorities for all hierarchical criteria had been determined, these priorities were synthesized to calculate and evaluate landslide hazard using the following formula:

$$landslide\_hazard=\sum _{j=1}^K\sum _{i=1}^{n}PSC\left(i\right)\times PC\left(j\right)$$

(4)

k = number of criteria (in our case, 10 criteria).

n = number of sub-criteria for each criterion.

3.2.3. Landslide Hazard Assessment

Equation (4) quantifies landslide hazard (e) using the AHP-based method, initially utilized for determining criteria and sub-criteria weights. For each criterion j (1 ≤ j ≤ K), sub-criteria contribute with weighted values (${WSC}_i$ for sub-criterion i, ${WSC}_j$ for criterion j). These sub-criteria align with the earlier established hierarchy. Multiplying sub-criterion weights by corresponding criterion weights captures their relative significance in the landslide hazard assessment, highlighting how both criteria and sub-criteria influence the overall evaluation. Summing across relevant criteria $(j=1 \mathrm{to} j=K)$ yields the total value, indicating the overall landslide hazard for each zone.

Greater aggregation values signify heightened susceptibility to landslides. The equation quantitatively evaluates and compares different regions’ landslide susceptibility, considering criteria weights, sub-criteria, and predefined susceptibility categories. Using hazard values, outcomes are categorized as ‘very low’, ‘low’ ‘moderate’, ‘high’, or ‘very high’, offering a clear assessment of an areas’ potential landslide vulnerability.

4. Implementation and Results

Our application serves as a powerful decision-making tool for predicting potentially harmful events. Utilizing advanced multiple criteria analysis techniques, specifically AHP (analytic hierarchy process) and FAHP (fuzzy analytic hierarchy process), the application was crafted in Python to streamline decision-making, especially in uncertain scenarios. It boasts an intuitive user interface for defining criteria, result visualization to highlight factor importance, and a consistency-checking feature to identify irregularities in matrices. This quantitative methodology facilitates systematic and objective landslide hazard evaluations across diverse study areas. Grounded in AHP, our approach provides a robust framework for assessing landslide risks. By considering a broad range of criteria and addressing uncertainties, our application offers crucial insights for decision-making and preventive measures. These outcomes have the potential to significantly enhance landslide risk management, safeguarding communities and critical infrastructure from the devastating consequences of such events.

4.1. Results and Analysis of Generated Maps

This section is dedicated to showcasing and analyzing the outcomes derived from the resolution methodology that relies on the analytic hierarchy process (AHP) and the fuzzy approach (FAHP) for evaluating landslide hazard.

Firstly, maps for each of the ten criteria were obtained. These criteria encompassed various geomorphological, geological, climatic, and human aspects that play a crucial role in assessing landslide hazard. Remote sensing and digitization tools, along with specific methods, were utilized to generate the corresponding maps.

Subsequently, susceptibility maps were generated using both AHP and FAHP. These maps consolidated the weights assigned to each criterion and sub-criterion within the hierarchy, thereby reflecting their respective significance in the comprehensive evaluation of landslide hazard.

4.2. Data Preparation

In this section, we will outline the various steps taken to obtain the different maps required for our study.

Data Origins

The thematic layers were sourced from diverse origins, encompassing a Landsat 8 image taken on 30 March 2023, a digital elevation model at a 30 m resolution dated 4 April 2023 (Figure 3), and a geological map of Mila provided by CRAAG (The Algerian Center for Astronomy, Astrophysics, and Geophysics Research).

• Rainfall Analysis

We commenced the analysis by retrieving rainfall data from the provided source at power.larc.nasa.gov/data-access-viewer/ (accessed on 19 May 2024). Following that, we proceeded to import this rainfall data into ArcGIS, with a specific focus on the study area in question. Next, we employed interpolation techniques to estimate rainfall values for regions that lacked coverage from weather stations. Lastly, we utilized ArcGIS to create a thematic map that visually illustrated the distribution of rainfall within the study area, utilizing a graduated color scale to emphasize any variations (Figure 4).

• Altitude Analysis

We began by acquiring the digital elevation model (DEM) from usgs.gov, ensuring that it was in an appropriate format, e.g., TIFF. Next, we utilized ArcGIS’s ‘Clip’ tool to isolate the study area from the DEM. We verified that the DEM aligned with the coordinate system of the study area, and then transformed the DEM into the desired unit (e.g., meters) using ArcGIS’s ‘Times’ tool. Finally, we adjusted the map’s symbology to visually represent altitude variations, employing graduated colors or contour lines (Figure 5).

• Hydrological Network Distance Analysis (buffer)

We first secured reliable data for the study area. We imported these data into ArcGIS, incorporating it as a vector layer, and introduced the digital elevation model (DEM) as a raster layer. We utilized the ‘Euclidean Distance’ tool to compute distances from each DEM cell to the nearest point on the hydrological network, generating a new raster layer. Finally, we created a visual distance map in ArcGIS, using graduated colors to depict varying distances within the study area (see Figure 6).

• Slope aspect analysis

We began by downloading the DEM map in the appropriate format from usgs.gov. Then, we employed ArcGIS’s ‘Clip’ tool to focus on the study area, ensuring efficient data management. We verified the referencing to align the DEM with the study area’s coordinate system. Using the ‘Aspect’ tool in ArcGIS’s Spatial Analyst toolbox, we computed the slope aspect for each DEM cell, specifying angle measurements in degrees or cardinal directions. Finally, we processed the map and adjusted the symbology, utilizing either a graduated color palette or a windrose to visualize the slope directions within the study area (Figure 7).

• Normalized Difference Moisture Index (NDMI)

To analyze water presence, we downloaded the Landsat 8 image for the study area from the USGS website and performed radiometric and atmospheric corrections on the image. We then calculated water index values for each pixel in the image. An NDMI map was generated, assigning graduated colors to represent water presence within the study area (Figure 8).

• Lithology analysis

In the context of lithology analysis, the process commenced with the acquisition of a digital version of the geological map, primarily for digitization purposes. Following this, pertinent information was assigned to individual lithological units, and a map legend was established. Subsequently, the digitization phase was initiated, utilizing the appropriate digitization tool tailored to the specific feature type (e.g., point, line, polygon) to represent the diverse lithological units accurately. Ultimately, a map was generated by associating distinct colors or patterns with each lithological type present within the study area, facilitating the visualization of various rock or soil compositions across the region (Figure 9).

• NDVI Analysis

For NDVI analysis, we started by acquiring a Landsat 8 image for the study area from the USGS website. After applying radiometric and atmospheric corrections to enhance image quality [40], we created an NDVI map, representing variations in vegetation density across the study area using a graduated color scheme (Figure 10).

• Drainage Density Analysis

We obtained hydrological network data from trusted sources for the study area. After importing these data into ArcGIS as a vector layer, we converted the digital elevation model (DEM) into points using the ‘Raster to Point’ tool, creating a point layer that represented the entire DEM. Then, we employed the ‘Spatial Join’ tool to associate DEM points with adjacent hydrological network features, attributing each DEM point with the characteristics of nearby hydrological network elements. The calculation of drainage density was performed by dividing the count of DEM points linked to each hydrological network element by the area it covered within the study zone. Finally, a descriptive drainage density map was generated in ArcGIS, illustrating variations in drainage densities across the study area (Figure 11).

• Slope Analysis

We began by downloading the digital elevation model (DEM) map for the study area from the USGS website. Using the ‘Clip’ tool within ArcMap of ArcGIS, we then tailored the downloaded DEM to match the study area’s boundaries. To ensure proper alignment, we verified the referencing system to match the DEM with the correct coordinate system of the study area. Subsequently, we employed the ‘Slope’ tool found in the Spatial Analyst toolbox of ArcGIS to compute the slope for each cell within the DEM. Finally, we established a map presentation, fine-tuning the symbology settings to visualize varying slope values effectively through an appropriate color scale or class intervals. This graphical representation effectively delineated areas with both low and high slopes on the map (Figure 12).

• Land Use Analysis

In terms of land use analysis, the process began with obtaining the Landsat 8 image from the USGS website, as shown. Subsequently, rigorous efforts were made to correct radiometric and atmospheric anomalies within the image. Carefully chosen training samples were selected to represent diverse land use categories, including water, urban areas, bare soil, and vegetation, guided by their distinctive visual attributes as depicted in Figure 13. Spectral signatures specific to these land use types were then extracted. Employing supervised classification methods, namely maximum likelihood (ML) and support vector machine (SVM), the image pixels were classified effectively. In terms of performance, SVM achieved an impressive 99.1014% accuracy with a Kappa coefficient of 0.9840, while ML demonstrated a notable 97.5922% accuracy with a Kappa coefficient of 0.9580. Finally, a coherent land use map was generated based on SVM results, employing a variety of colors and symbols to represent distinct land use categories for straightforward interpretation (Figure 13).

4.3. Slide Inventory Map

Information about the location, date, and types of landslides or mass movements in a region can be obtained via a slide inventory map. The preparation methods for such maps vary based on project goals. This map plays a vital role in validating susceptibility map outcomes. It verifies whether areas identified as highly susceptible to landslides on susceptibility maps align with actual landslide occurrences, enhancing the credibility of susceptibility maps (Figure 14). Additionally, it aids in validating the reliability and effectiveness of landslide hazard assessment methods by comparing results with susceptibility maps. Furthermore, it can uncover areas with past landslides that were not considered high-risk on susceptibility maps, helping improve assessment models. Finally, the information it provides supports the planning of preventive measures and risk management policies, prioritizing actions in high-risk areas. The landslide inventory map for this study was obtained in collaboration with the CRAAG center (Center for Research in Astronomy, Astrophysics, and Geophysics).

4.4. Susceptibility Map Generation

To create the susceptibility maps, we followed a systematic process after developing the criteria maps. Initially, we divided each criterion map into representative subclasses in consultation with domain experts, enhancing precision in characterizing areas based on specific criteria-related traits. Subsequently, we assigned weights to these subclasses, computed during the AHP method application, reflecting the criteria’s relative importance in overall landslide susceptibility assessment. Utilizing ArcGIS’s Map Algebra tool, we aggregated criteria by sub-criteria, applying our established landslide susceptibility assessment (Equation (4)). This enabled us to create a comprehensive susceptibility map for the entire study area. Furthermore, we collaborated with specialists to categorize the susceptibility map into five classes, simplifying decision-making and interpretation. These classes ranged from “very low” to “very high” susceptibility. In the realm of land use analysis, the process began with obtaining the Landsat 8 image from the USGS website, as shown in Figure 3. Subsequently, rigorous efforts were made to correct radiometric and atmospheric anomalies within the image.

Carefully chosen training samples were selected to represent diverse land use categories, including water, urban areas, bare soil, and vegetation, guided by their distinctive visual attributes, as depicted in Figure 15 and Figure 16. Spectral signatures specific to these land use types were then extracted. Employing supervised classification methods, namely maximum likelihood (ML) and support vector machine (SVM), the image pixels were classified effectively. In terms of performance, SVM achieved an impressive 99.1014% accuracy with a Kappa coefficient of 0.9840, while ML demonstrated a notable 97.5922% accuracy with a Kappa coefficient of 0.9580. Finally, a coherent land use map was generated based on SVM results, employing a variety of colors and symbols to represent distinct land use categories for straightforward interpretation. Both the AHP and FAHP methods followed these steps for comparison. See the susceptibility maps below.

4.5. Results and Interpretation

• Susceptibility Class Analysis: Susceptibility maps indicate the vulnerability level of areas to landslides, typically categorized into five classes ranging from “very low” to “very high”. These classes help distinguish high-risk areas from those less exposed to landslides.
• Identification of Critical Zones: The most vulnerable areas to landslides are those classified as having “high” or “very high” susceptibility. In our case, northern communes of Mila, such as Tassadane Haddada, Minar Zarza, Tassala Lematai, Amira Arrese, Terrai Bainen, Chigara, Hamala, Grarem Gouga, Elaydi Barbes, and Derrahi Bousselah, were identified as critical areas to closely monitor.
• Analysis of Safe Zones: On the other hand, some communes like Ben Yahia Abderrahmane, Ain Mlouk, the north of Tajenanet, and Chalgoum Elaid exhibited low susceptibility to landslides. These areas are considered relatively safe and less exposed to landslide risks.
• Comparison between Methods: In this study, two methods, AHP and FAHP, were employed to generate landslide susceptibility maps. Overall, the results from both methods were satisfactory, as areas with past landslides were correctly identified as having high or very high susceptibility. However, an important observation was made: the FAHP method seemed to handle uncertainty and special assessment errors better. It offered a more realistic transition in values, evident in the gradual green color transitions on the map. In contrast, the AHP-produced map shows abrupt changes between susceptibility classes. All in all, we find that FAHP is more adept at handling uncertainties inherent in landslide risk assessment, primarily due to its ability to incorporate fuzzy logic into the decision-making process. This feature proves to be essential in environments with incomplete information, thereby enhancing the reliability of our susceptibility maps; this is clearly shown in the green part of Figure 15 and Figure 16.
• Influence Analysis: The overall susceptibility map considered all criteria, setting it apart from individual criterion maps. However, some criteria had a stronger influence than others, notably slope, elevation, vegetation, moisture, hydrographic network, and rainfall. These factors play a crucial role in landslide hazard assessment.
• Field Validation: Validation of susceptibility maps was carried out by comparing the results with the landslide inventory map. This verified the accuracy and reliability of the produced maps.
• Decision-Making: Susceptibility maps provide essential information to decision-makers for identifying the most landslide-vulnerable areas. They enable informed decisions in urban planning, territorial development, and risk management.
• Limitations and Uncertainties: Like any analytical method, susceptibility maps have associated limitations and uncertainties. In our case, the subjectivity of specialists in assigning weights to criteria and sub-criteria may have influenced results. However, efforts were made to minimize these uncertainties by using statistical data provided by CRAAG to determine weights.

The tables below present numerical results derived from the application of both the analytical hierarchy process (AHP) and fuzzy AHP methods.

4.5.1. AHP Results

In the analytic hierarchy process (AHP) analysis, each criterion was weighted to reflect its importance in the decision-making process. The normalized difference moisture index (NDMI) and the normalized difference vegetation index (NDVI) carried moderate weights of 0.098 and 0.085, respectively, suggesting that they are considered of moderate importance. Slope aspect, with a weight of 0.118, and buffer, with a weight of 0.126, are viewed as more crucial, but not the most critical factors. In contrast, rainfall and altitude were weighted significantly higher at 0.176 and 0.15, highlighting their substantial impact on decisions within the analyzed context, likely due to their relevance in environmental or geographical considerations. Conversely, land use and network density had lesser weights of 0.045 and 0.058, indicating that they are less influential. The slope and lithology were also considered, with weights suggesting they have a noticeable impact, but are not dominant in the decision matrix (

Table 1. Weights and coherence ratios (criterion vs. criterion) (AHP method).

Criterion	NDMI	Slope Aspect	Land Use	Rainfall	NDVI	Buffer	Network Density	Altitude	Slope	Lithology	CR
Weight	0.098	0.118	0.045	0.176	0.085	0.126	0.058	0.15	0.051	0.092	0.01

).

The coherence ratio (CR) of 0.01 demonstrates a high degree of consistency in the judgments made during the AHP process, underlining the reliability of the decision-making structure. This low CR indicates that the judgments across the criteria are consistent enough to support the decision framework, ensuring that significant factors like rainfall and altitude are given appropriate consideration. This methodical weighting and consistent assessment reflect a structured approach to prioritizing factors crucial for scenarios in environmental management or urban planning.

Table 2. Weights and coherence ratios of sub-criteria (AHP method).

S-Cri. (Land Use)	Weight	S-Cri. (Slope)	Weight
Bar Soil	0.042	0–5	0.143
Urban	0.347	5–10	0.571
Water	0.234	10–59	0.286
Vegetation	0.377	--	--
Coherence Ratio	0.04	Coherence Ratio	0.00

and

Table 3. Weights and coherence ratios of sub-criteria (AHP method).

S-Cri. (Rainfall)	Weight	S-Cri. (Buffer)	Weight
300–400	0.016	0–50	0.114
400–500	0.067	50–100	0.128
500–600	0.125	100–150	0.116
600–700	0.127	150–200	0.114
700–800	0.096	200–250	0.118
800–900	0.216	250–300	0.089
900–1000	0.071	300–350	0.178
1000–1200	0.051	350–400	0.048
1200–1400	0.177	400–450	0.04
1400–1600	0.018	450–500	0.033
1600–1800	0.018	500–550	0.011
More than 1800	0.018	550–600	0.011
Coherence Ratio	0.1	Coherence Ratio	0.07

of the AHP results present comprehensive weights and coherence ratios for a variety of environmental and geographical sub-criteria, reflecting their relative importance and consistency in evaluation. In land use analysis, urban areas were highlighted with the highest weight of 0.347 due to their significant role, closely followed by vegetation at 0.377, emphasizing a focus on green spaces or agricultural areas. Water and bar soil received lesser priority, with weights of 0.234 and 0.042, respectively, and a coherence ratio of 0.04 in this category ensured consistent judgments. Slope criteria were weighted progressively higher from 0.143 for gentle slopes to 0.286 for steeper slopes, addressing erosion and construction safety concerns with a perfect coherence ratio of 0.00. Rainfall and buffer evaluations also varied, with rainfall at 800–900 mm marked as most crucial at a weight of 0.216, while buffer distances showed decreasing importance with increasing distance, reflected in slightly elevated but acceptable coherence ratios. Sub-criteria such as NDMI and NDVI emphasize the importance of moderate levels of moisture and vegetation density. North and southwest slope aspects received higher prioritization due to their impact on sunlight exposure and wind patterns. Altitude gained more weight at higher elevations, indicating its crucial role in environmental assessments. The molassic series under lithology was notably weighted the highest, underlining its significant impact on land use due to soil properties. Network density was also emphasized to ensure optimal levels are maintained. Overall, the low coherence ratios across all categories highlight the reliability and structured approach of the AHP, crucial for effective decision-making in environmental management and urban planning (

Table 4. Weights and coherence ratios of sub-criteria (AHP method).

S-Cri. (NDMI)	Weight	S-Cri. (NDVI)	Weight
−0.6, −0.4	0.221	−0.44, −0.1	0.222
−0.4, −0.2	0.158	−0.1, 0.0	0.222
−0.2, 0.0	0.203	0.0, 0.2	0.15
0.0, 0.2	0.131	0.2, 0.4	0.142
0.2, 0.4	0.121	0.4, 0.6	0.205
0.4, 0.6	0.131	0.6, 0.8	0.04
0.6, 0.8	0.018	0.8, 0.84	0.018
0.8, 1.0	0.018	--	--
Coherence Ratio	0.02	Coherence Ratio	0.01
S-Cri. (Slope Aspect)	Weight	S-Cri. (Altitude)	Weight
Flat	0.077	65–200	0.043
N	0.072	200–400	0.103
N-E	0.077	400–600	0.047
E	0.077	600–800	0.018
S-E	0.085	800–1000	0.032
S	0.153	1000–1200	0.93
S-W	0.153	1200–1400	0.194
W	0.153	1400–1500	0.284
N-W	0.153	1500–1580	0.185
Coherence Ratio	0.00	Coherence Ratio	0.01
S-Cri. (Lithology)	Weight	S-Cri. (N. Density)	Weight
Alluvialscree land	0.136	0	0.014
Metamorphic Terrain	0.076	0–0.016	0.073
Micaschist Gneiss	0.332	0.016–0.435	0.111
Molassic Series	0.085	0.435–0.696	0.15
Yellow Sandy Marls Marl and Marly Limestone	0.146	0.696–0.991	0.168
Gypsum complex	0.142	0.991–1.351	0.168
Clays	0.015	1.351–2.029	0.09
Stony plateau alluvium	0.015	2.029–3.029	0.174
Flyschs clay microbreccia argillites	0.053	3.029–4.930	0.052
Coherence Ratio	0.08	Coherence Ratio	0.02

).

4.5.2. FAHP Results

Table 5. Weights and coherence ratios (criterion vs. criterion) (FAHP method).

Criterion	NDMI	Slope Aspect	Land Use	Rainfall	NDVI	Buffer	Network Density	Altitude	Slope	Lithology	CR
Weight	0.098	0.115	0.048	0.172	0.087	0.12	0.063	0.146	0.057	0.094	0.03

, utilizing the fuzzy analytic hierarchy process (FAHP), assigns weights and coherence ratios to key decision-making criteria. Rainfall and altitude were prioritized with weights of 0.172 and 0.146, respectively, due to their impact on water availability and climate. Slope aspect and buffer zones also held significant weights of 0.115 and 0.12, underscoring their environmental and developmental importance. Moderate weights were given to NDMI and NDVI for their roles in moisture and vegetation assessment, while land use, network density, slope, and lithology had lower weights, indicating their lesser but still pertinent effects. A coherence ratio of 0.03 confirmed the high consistency and accuracy of the evaluations, essential for effective decision-making (Table 5).

In

Table 6. Weights and coherence ratios of sub-criteria (FAHP method).

S-Cri. (Land Use)	Weight	S-Cri. (Slope)	Weight
Bare Soil	0.039	[0, 5]	0.15
Urban	0.34	[5, 10]	0.552
Water	0.248	[10, 595]	0.299
Vegetation	0.373	--	--
Consistency Ratio	0.06	Consistency Ratio	0.07
S-Cri. (Rainfall)	Weight	S-Cri. (Buffer)	Weight
300–400	0.016	0–50	0.118
400–500	0.055	50–100	0.12
500–600	0.137	100–150	0.131
600–700	0.141	150–200	0.118
700–800	0.101	200–250	0.121
800–900	0.214	250–300	0.09
900–1000	0.072	300–350	0.173
1000–1200	0.036	350–400	0.048
1200–1400	0.181	400–450	0.036
1400–1600	0.015	450–500	0.028
1600–1800	0.015	500–550	0.009
More than 1800	0.015	550–600	0.009
Consistency Ratio	0.09	Consistency Ratio	0.08

, weights and consistency ratios are assigned to various sub-criteria critical for environmental and urban planning decisions via the FAHP approach. In the land use category, urban areas were weighted most heavily at 0.34, reflecting their importance in urban development, followed by vegetation at 0.373, highlighting its significance in sustainability efforts. Water and bare soil were also considered with weights of 0.248 and 0.039, respectively, indicating their varied impacts on decision-making. The slope criterion differentiated importance across slope ranges, emphasizing the moderate slopes [15,16,17,18,19,20,21,22,23,24,25] with a weight of 0.552, due to their effects on erosion and safety. Rainfall weights varied from 0.016 to 0.214 across different intensity ranges, with the highest emphasis on the 800–900 mm range, crucial for flood risk management and agricultural planning. Buffer zones decreased in importance with distance, with 300–350 m highlighted at 0.173, reflecting its role in zoning for environmental protection. Consistency ratios for these criteria ranged from 0.06 to 0.09, underscoring a reasonable level of judgment consistency, which bolsters the reliability of the evaluations. This comprehensive analysis ensures stakeholders are well-informed of the prioritized factors essential for resource management and project planning in sensitive environmental and urban contexts.

All in all, susceptibility maps obtained through AHP and FAHP offer valuable insights for assessing landslide hazards. They help identify critical areas, compare different methods, inform decisions, and deliver a better understanding of factors influencing landslide susceptibility.

5. Conclusions and Remarks

In the present research, we have succeeded in presenting a quantitative geomorphological imaging study of the Mila region in Algeria. Using multicriteria analysis, a powerful operations research method, we carried out a thorough examination. In our investigation, the AHP and FAHP-based techniques proved to be helpful in determining the likelihood of landslides. Our results, derived from the susceptibility maps, have facilitated a more profound comprehension of the vulnerability of diverse geographical areas to landslides. These susceptibility maps were generated by considering multiple criteria, including slope, altitude, vegetation, moisture, and the hydrographic network. The mapping of these criteria was conducted with great precision using remote sensing and digitization tools, thereby ensuring an accurate and comprehensive representation of geomorphological, geological, climatic, and human factors.

We successfully delineated critical zones, which were identified as the most susceptible to landslides. This delineation empowers decision-makers to take appropriate measures for urban planning, land management, and risk mitigation. The comparison between AHP and FAHP revealed that both methods effectively identify high-risk areas, with FAHP offering superior management of uncertainty and accommodating subjective expert assessments. This leads to a more realistic evaluation of landslide susceptibility. It is important to acknowledge certain limitations and uncertainties associated with this approach, particularly the subjectivity in assigning criteria weights. However, we have mitigated these uncertainties by relying on statistical data for weight determination and by validating susceptibility maps through fieldwork.

In the future, further research should explore other regions, in particular those prone to earthquakes, and develop robust tools to predict such phenomena. A comparative study of various multicriteria optimization techniques like PROMETHEE, ELECTRE I, and ELECTRE II holds promise. A comparative analysis of different multicriteria analysis methods can significantly enhance our ability to predict landslide-prone areas.