Landslide Risk Assessments through Multicriteria Analysis

1. Introduction

Landslides, stemming from factors like human activities, climate, and geomorphology, pose a global threat to human settlements and infrastructure. Despite their frequency, landslides are inadequately studied compared to other disasters, possibly due to mitigating effects from concurrent events. The expanding global population increases vulnerability to landslide disasters, especially as communities and regions develop, heightening the risk of catastrophic losses. The instability of slopes, influenced by human activities or not, poses inherent dangers to people and their surroundings. Our focus is on evaluating the cumulative decisions contributing to this vulnerability, often challenging to alter, and with the potential to significantly increase or decrease 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, explored in various publications (see [1,2,3,4,5,6]), form the framework for our study.

Researchers worldwide, as exemplified by studies in the Haraz watershed and Brazil, have utilized multicriteria analysis, such as the Analytic Hierarchy Process (AHP), to map landslide susceptibility. In Portugal, the AHP was employed to study landslides in a coastal area from Figueira da Foz to Nazaré. These studies underscore the importance of multicriteria analysis and highlight the use of remote sensing tools and specialized software for image processing and spatial data collection in assessing landslides globally. This article focuses on the Wilaya of Mila in Algeria, exploring its geological features and past landslide history. The occurrence of a magnitude 4.9 earthquake on the 7th of August 2020 in the Mila region underscores the importance of seismic and geotechnical factors.

To address landslides, we draw on previous research by Pourghasemi, de Brito, and Ramos, focusing on Iran, Brazil, and Portugal, respectively (see [7,8,9]). Operational research techniques, particularly multicriteria analysis, will be utilized, incorporating both the Analytic Hierarchy Process (AHP) and Fuzzy Analytic Hierarchy Process (FAHP) methodologies. These approaches will evaluate landslide hazards and associated risks in the Wilaya of Mila. Our primary goal is to propose a comprehensive approach using both deterministic and fuzzy frameworks, leading to a robust natural risk assessment model validated through rigorous testing. The paper is structured as follows:

Section 1 provides an overview of the research problem, encompassing various landslide types, causes, real-world instances, and methods for assessing landslide risk. The subsequent section introduces essential tools, such as remote sensing (e.g., ArcGIS, ENVI, Landsat), 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 Mila province and investigating preventive strategies for each high-risk zone. The paper concludes by summarizing findings and providing future perspectives [10,11,12,13,14,15,16,17].

1.1. Generalities on Landslides

Landslides involve the movement of soil on natural or artificial slopes. They come in various types based on shape and movement characteristics, including rotational, planar, and translational landslides. The term “mass movements of materials” in this context encompasses soils, rocks, and fractured rocks (defined in [18]). Landslide occurrence depends on factors like slope, geology, vegetation, precipitation, and human activity.

Different types of landslides are categorized by 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) (see [19]).

Landslides can result from various factors categorized into two main groups: natural causes, driven 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 observed in the Nilgiris region in India, leading to increased landslide occurrences, [2,3,5].

1.2. Landslides Worldwide

Landslides have inflicted severe damage globally. The 1963 Vajont landslide in Italy, one of the deadliest in history, claimed over 2,000 lives [20]. Similarly, the 2014 Oso landslide in the United States led to 43 fatalities, with heavy rainfall contributing to soil saturation and extensive devastation.

China faced foods and landslides in Zhouqu in 2010, claiming more than 1,200 lives and causing 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 also seen numerous landslides, with the northern regions particularly susceptible due to their complex geomorphology and geological structures. Notable incidents include landslides in Bejaia and Tlemcen, causing significant damage and disruptions to urban areas and economic routes. The wilaya of Mila is located in the north of Algeria, a region that is prone to seismic activity. Permanently, a certain number of tremors occur there, although most of them are not felt by the population. On Friday, August 7, 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 wilaya of Mila. This seismic event resulted in landslides in the Mila region.

1.2.1. Type of Approaches

Landslide hazard assessment employs two primary methodologies: geomorphologicalfield 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 analyze historical events in landslide-prone areas through bivariate and multivariate statistical analyses (see [21]). The former calculates density by merging parameters and landslide maps, favored by geologists but complex due to variable correlations. The latter categorizes 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 (see [22,23]).

2. Preliminary Tools

This section introduces the essential tools necessary for our study on landslide hazard assessment. We will discuss remote sensing and the specialized software utilized for image processing and spatial data collection (Landsat, ENVI, and ArcGIS). Furthermore, we will introduce operations research concepts, with a particular emphasis on multicriteria optimization.

2.1. Image Processing

In this investigation, we explore Remote Sensing, a contactless method for capturing object characteristics using instruments like cameras, lasers, radars, or sonars, often deployed on remote platforms like satellites. Employing 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 agriculture, forests, mountains, and coastal areas, offering insights into cultivated fields and, at times, human habitation. Notably, 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. Additionally, the study addresses Image Segmentation, a technique dividing images 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 suchas 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 ([24]).

2.2. Multi-Criteria Decision-Making (MCDM)

2.2.1. Criteria in Decision-Making

A criterion, denoted by function g over action set A with ordered values, reflects user preferences (see [25]). 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.2.2. Complexity in Multicriteria Decision-Making

Multicriteria decision-making is inherently complex, characterized by strategic challenges arising from conflicting criteria, uncertainty, and data limitations. Navigating mathematical ambiguity is a central challenge, as no single optimal solution typically prevails across all criteria simultaneously. Unlike conventional optimization problems, multicriteria contexts demand a nuanced approach shaped by individual attributes, decision-making context, problem formulation, and chosen support method. Objective solutions in multicriteria ranking problems are rare, and these problems vary based on decision consequences, alternative modeling, and decision context ([26]). Multicriteria analysis methods aim to assist decision-makers in identifying the most suitable solution, considering subjective elements such as personality, circumstances, problem framing, and the chosen method, with emphasis on the approach used and the desired outcome [27]. Roy’s framework categorizes multicriteria decision-making into four problematic categories: Choice, Sorting, Ranking, and Description Problematics, each serving distinct decision-making purposes.

2.3. Multicriteria Analysis Methods

In this section, we focus on two main approaches: the Analytic Hierarchy Process (AHP) and the Fuzzy Analytic Hierarchy Process (FAHP). Widely used in landslide evaluation literature, these methods suit our study.

In multicriteria decision-making, the “partial aggregation” or “synthetic dominance” approach sacrifices result clarity, accommodating non-quantitative criteria and diverseunits through pairwise comparisons. Commonly used methods like ELECTRE and PROMETHEE adapt to various decision contexts [28,29,30,31].

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 (see [4]). 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 (see [27]). Another method is TOPSIS, aiming to select the best alternative by evaluating the distance to the ideal and negative ideal alternatives (Yoon & Hwang, 1981). 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), pioneered by Thomas Saaty in the 1970s, assists in decision-making within intricate scenarios featuring numerous hard-to-prioritize criteria. AHP simplifies these challenges by structuring problems hierarchically, employing subjective judgments from stakeholders in comparison matrices to ascertain the relative importance of criteria and alternatives, ultimately yielding the optimal solution based on predefined criteria. Its fundamental principles, outlined by Saaty in his 1986 work, encompass decomposition, comparative judgments, and priority synthesis. AHP systematically constructs hierarchies, determines priorities, and ensures logical coherence checks, effectively organizing and evaluating problem components. The process begins by establishing a precise hierarchical structure encompassing criteria, sub-criteria, and alternatives.

Establishment of Priorities

The initial step in determining element priorities within a decision problem involves conducting pairwise comparisons at the same hierarchical level. These comparisons are executed through logical and experiential evaluations, utilizing a matrix framework. This process establishes concrete comparative judgments, gathering relevant data to assess relationships between elements and criteria. Numerical values are assigned to populate the pairwise comparison matrix, representing the relative importance of one element compared to another concerning a specific property.

To ensure the relevance of the comparisons made in the matrix, we normalize it. After normalizing the pairwise comparison matrix, we calculate eigenvectors by averaging the column elements. This approach yields accurate values for normalized eigenvectors, representing criterion weights with improved precision.

2.3.1. Evaluating Consistency

AHP method employs a consistency ratio to assess the coherence of evaluations. A ratio exceeding 0.1 suggests that evaluations might be somewhat arbitrary and may requirerevisions (see [32]). The Consistency Index (CI) is used to quantify this, calculated using the eigenvalue method

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

(1)

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

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

(2)

CR signifies the consistency ratio, CI denotes the consistency Index, and RCI stands for the Random Consistency Indices.

2.3.2. Fuzzy Analytical Hierarchy Process (F-AHP)

The Fuzzy Analytical Hierarchy Process (F-AHP) enhances the traditional AHP by incorporating fuzzy logic, thus refining the multi-criteria decision-making framework. Unlike the direct ranking system of AHP, F-AHP utilizes linguistic variables and triangular fuzzy numbers for a more detailed and nuanced comparison process, as outlined in [33]. This approach, rooted in the pioneering works of van Laarhoven and Pedrycz [34], and further developed by Buckley [35] and Chang [36], 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. The pairwise contribution matrix, utilizing fuzzy triangular numbers as described in Equation (3), encodes the preferences of decision-makers. Within this matrix, d˜krepresents 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 F-AHP for calculating the criteria’s normalized weights are systematically outlined below.

3. Landslide Hazard Assessment Criteria and Methodology

This section outlines criteria for landslide hazard assessment, derived from remote sensing and expert collaboration. Employing satellite analysis and tools like ArcGIS, ENVI, and Landsat, we collect 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 consistently link rainfall to landslides ([1]). Snowmelt also contributes, releasing water as temperatures rise, causing ground movements (see [37]).

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) ([38]) underscores drainage’s profound effect ([39]). 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, water body proximity and surface-groundwater interactions 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 create a hierarchy of criteria, sub-criteria, and alternatives (see Figure 1). Level 0 focuses on evaluating landslide hazard. Decision criteria include 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 employ AHP (deterministic) or FAHP (fuzzy) to compare elements within the same hierarchical level. The objective is to assess the relative importance of each element in relation to those at higher levels. We use appropriate scales, which can be quantitative, qualitative, or expert-based, depending on the study context, to evaluate the preference or significance between element pairs. To ensure consistency, comparison matrices must adhere to fundamental principles. This step is crucial forconstructing 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.1.1. Priority Determination

In this phase, we determine 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 are calculated, and the resulting normalized vector values provide a clear indication of the importance of criteria, where larger values denote greater significance.

3.2.2. Priority Synthesis

• Once priorities for all hierarchical criteria have been determined, these priorities aresynthesized 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, initiallyutilized 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 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 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 is 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 part 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 are obtained. These criteria encompass 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, are utilized to generate the corresponding maps.

Subsequently, susceptibility maps are generated using both AHP and FAHP. These maps consolidate 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 to obtain the different maps required for our study:

4.2.1. Data Origins

The thematic layers have been sourced from diverse origins, encompassing a Landsat 8 image taken on March 30, 2023 (depicted in Figure 2), a Digital Elevation Model at a 30m resolution dated April 4, 2023 (illustrated in Figure 2), and a Geological map of Mila provided by CRAAG (The Algerian Center for Astronomy, Astrophysics, and Geophysics Research).

• For rainfall Analysis

We commence the analysis by retrieving rainfall data from the provided source at power.larc.nasa.gov/data-access-viewer/. Following that, we proceed to import this rainfall data into ArcGIS, with a specific focus on the study area in question. Next, we employ interpolation techniques to estimate rainfall values for regions that lack coverage from weather stations. Lastly, we utilize ArcGIS to create a thematic map that visually illustrates the distribution of rainfall within the study area, utilizing a graduated color scale to emphasize any variations. See the obtained map in Figure 3. For altitude analysis Start by acquiring the Digital Elevation Model (DEM) from usgs.gov, ensuring it’s in the appropriate format, like TIFF. Next, utilize ArcGIS’s ‘Clip’ tool to isolate the study area from the DEM. Verify that the DEM aligns with the coordinate system of the study area. Transform the DEM into the desired unit (e.g., meters) using ArcGIS’s ‘Times’ tool. Finally, adjust the map’s symbology to visually represent altitude variations, employing graduated colors or contour lines (see Figure 3).

In hydrological network distance analysis (buffer), we first secure reliable data for the study area. We import this data into ArcGIS, incorporating it as a vector layer, and introduce the Digital Elevation Model (DEM) as a raster layer. We utilize 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 create a visual distance map in ArcGIS, using graduated colors to depict varying distances within the study area (Figure 3).

For slope aspect analysis, begin by downloading the DEM map in the appropriate format from usgs.gov. Then, employ ArcGIS’s ‘Clip’ tool to focus on the study area, ensuring efficient data management. Verify the referencing to align the DEM with the study area’s coordinate system. Use the ‘Aspect’ tool in ArcGIS’s Spatial Analyst toolbox to compute slope aspect for each DEM cell, specifying angle measurements in degrees or cardinal directions. Finally, process the map and adjust the symbology, utilizing either a graduated color palette or a windrose, to visualize the slope directions within the study area (see Figure 3).

• For NDMI

To analyze water presence, we download the Landsat 8 image from the USGS website for the study area (Figure 2) and we perform radiometric and atmospheric corrections on the Landsat 8 image. We calculate then water index values for each pixelhtp in the image. An NDMI map is generated, assigning graduated colors to represent water presence within the study area (see Figure 3).

In the context of lithology analysis, the process commences with the acquisition of a digital version of the geological map, primarily for digitization purposes. Following this, pertinent information is assigned to individual lithological units, and a map legend is established. Subsequently, the digitization phase is 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 is 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 (refer to Figure 3).

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

For drainage density analysis, we obtain hydrological network data from trusted sources for the study area. After importing this data into ArcGIS as a vector layer, we convert the Digital Elevation Model (DEM) into points using the ‘Raster to Point’ tool, creating a point layer that represents the entire DEM. Then, we employ 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 is performed by dividing the count of DEM points linked to each hydrological network element by the area it covers within the study zone. Finally, a descriptive drainage density map is generated in ArcGIS, illustrating variations in drainage densities across the study area (see Figure 4).

• For slope analysis

We begin by downloading the Digital Elevation Model (DEM) map for the study area from the USGS website. Using the ‘Clip’ tool within ArcMap of ArcGIS, wethen tailor the downloaded DEM to match the study area’s boundaries.

To ensure proper alignment, we verify the referencing system to match the DEM with the correct coordinate system of the study area. Subsequently, we employ the ‘Slope’ tool found in the Spatial Analyst toolbox of ArcGIS to compute the slope for each cell within the DEM. Finally, we establish 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 delineates areas with both low and high slopes on the map (refer to Figure 5).

For Land Use Analysis

In the realm of land use analysis, the process began with obtaining the Landsat 8 image from the USGS website, as shown in Figure 2.

Subsequently, rigorous efforts weremade 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 7. 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 (see Figure 5).

4.3. Slide Inventory Map

Information about the location, date, and types of landslides or mass movements in a region. Its preparation methods 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. 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 is obtained in collaboration with the CRAAG center (Center for Research in Astronomy, Astrophysics, and Geophysics). See Figure 7 for more details.

4.4. Susceptibility Map Generation

To create susceptibility maps, we follow a systematic process after developing the criteria maps. Initially, we divide each criterion map into representative subclasses in consultation with domain experts, enhancing precision in characterizing areas based on specific criteria-related traits. Subsequently, we assign 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 aggregate criteria by sub-criteria, applying our established landslide susceptibility assessment equation 4. This enables us to create a comprehensive susceptibility map for the entire study area. Furthermore, we collaborate with specialists to categorize the susceptibility map into 5 classes, simplifying decision-making and interpretation. These classes range 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 2. Subsequently, rigorous efforts were made to correct radiometric and atmospheric anomalies within the image.

Figure 6. Slide Inventory Map.

Figure 6. Slide Inventory Map.

Carefully chosen training samples were selected to represent diverse land use categories, including water, urban areas, bare soil, and vegetation, guided by their distinctive visualattributes as depicted in Figure 7. 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 (see Figure 5). Both the AHP and FAHP methods follow these steps for comparison. See the susceptibility maps below.

4.5. Results and Interpretation

• Susceptibility Class Analysis: Susceptibility maps indicate the vulnerabilitylevel of areas to landslides, typically categorized into five classes ranging from “verylow” to “very high.” These classes help distinguish high-risk areas from those lessexposed to landslides.
• Identification of Critical Zones: The most vulnerable areas to landslides arethose classified as having “high” or “very high” susceptibility. In our case, northerncommunes of Mila, such as TassadaneHaddada, Minar Zarza, TassalaLematai, Amira Arrese, TerraiBainen, Chigara, Hamala, GraremGouga, Elaydi Barbes, andDerrahiBousselah, are identified as critical areas to closely monitor.
• Analysis of Safe Zones:On the other hand, some communes like Ben YahiaAbderrahmane, Ain Mlouk, the north of Tajenanet, and ChalgoumElaid exhibit lowsusceptibility to landslides. These areas are considered relatively safe and less exposed tolandslide risks.
• Comparison between Methods: In this study, two methods, AHP and FAHP, were employed to generate landslide susceptibility maps. Overall, results from bothmethods are satisfactory, as areas with past landslides are correctly identified ashaving high or very high susceptibility. However, an important observation wasmade: the FAHP method seems to handle uncertainty and specialist assessmenterrors better. It offers a more realistic transition in values, evident in the gradualcolor transitions on the map. In contrast, the AHP-produced map shows abruptchanges between susceptibility classes.
• Influence Analysis: The overall susceptibility map considers all criteria, setting it apart from individual criterion maps. However, some criteria have astronger influence than others, notably slope, elevation, vegetation, moisture, hydrographic network, and rainfall. These factors play a crucial role in landslidehazard assessment.
• Field Validation: Validation of susceptibility maps is carried out by comparing theresults with the landslide inventory map. This verifies the accuracy and reliabilityof the produced maps.
• Decision-Making: Susceptibility maps provide essential information to decisionmakers 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 can influence results. However, efforts have been 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.

Figure 7. Susceptibility Map Generation.

Figure 7. Susceptibility Map Generation.

4.5.1. AHP Results

In the Analytic Hierarchy Process (AHP) analysis, each criterion is weighted to reflect its importance in the decision-making process. The Normalized Difference Moisture Index (NDMI) and the Normalized Difference Vegetation Index (NDVI) carry moderate weights of 0.098 and 0.085, respectively, suggesting 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 are 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 have lesser weights of 0.045 and 0.058, indicating they are less influential. The Slope and Lithology are also considered, with weights suggesting they have a noticeable impact but are not dominant in the decision matrix.

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 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

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

Table 2 and Table 3 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 are 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 receive lesser priority, with weights of 0.234 and 0.042 respectively, and a coherence ratio of 0.04 in this category ensures consistent judgments. Slope criteria are 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 vary, with rainfall at 800-900 mm marked as most crucial at a weight of 0.216, while buffer distances show 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 South-West slope aspects receive higher prioritization due to their impact on sunlight exposure and wind patterns. Altitude gains more weight at higher elevations, indicating its crucial role in environmental assessments. The molassic series under lithology is notably weighted the highest, underlining its significant impact on land use due to soil properties. Network density is 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.

4.5.2. FAHP Results

Table 5, utilizing the Fuzzy Analytic Hierarchy Process (FAHP), assigns weights and coherence ratios to key decision-making criteria. Rainfall and Altitude are 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 hold significant weights of 0.115 and 0.12, underscoring their environmental and developmental importance. Moderate weights are given to NDMI and NDVI for their roles in moisture and vegetation assessment, while Land Use, Network Density, Slope, and Lithology have lower weights, indicating their lesser but still pertinent effects. A coherence ratio of 0.03 confirms the high consistency and accuracy of the evaluations, essential for effective decision-making.

In Table 6,weights and consistency ratios are assigned to various sub-criteria critical for environmental and urban planning decisions via FAHP approach. In the Land Use category, Urban areas are 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 are also considered with weights of 0.248 and 0.039, respectively, indicating their varied impacts on decision-making. The Slope criterion differentiates importance across slope ranges, emphasizing the moderate slopes [5,10] with a weight of 0.552 due to their effects on erosion and safety. Rainfall weights vary 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 decrease in importance with distance, with 300-350 m highlighted at 0.173, reflecting its role in zoning or environmental protection. Consistency ratios for these criteria range 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.

In conclusion, susceptibility maps obtained through AHP and FAHP offer valuable insights for assessing landslide hazard. They help identify critical areas, compare different methods, make informed decisions, and gain a better understanding of factors influencing landslide susceptibility.

5. Conclusions and Remarks

In this paper, we present a quantitative geomorphological imaging study of the Mila region in Algeria. We conducted a detailed analysis using a robust operations research tool, Multicriteria analysis. The methodology based on AHP and FAHP has demonstrated its effectiveness in assessing landslide hazards in our study. The results derived from susceptibility maps have provided us with a deeper understanding of various areas’ vulnerability to landslides. These susceptibility maps were generated by considering multiple criteria, including slope, altitude, vegetation, moisture, and the hydrographic network. We meticulously mapped these criteria using remote sensing and digitization tools, ensuring accurate and comprehensive representation of geomorphological, geological, climatic, and human factors.

We successfully delineated critical zones, identifying them 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 has 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.

Looking ahead, future research should explore other regions, particularly 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.