Integrated Machine Learning Approaches for Landslide Susceptibility Mapping Along the Pakistan–China Karakoram Highway

Abstract

The effectiveness of data-driven landslide susceptibility mapping relies on data integrity and advanced geospatial analysis; however, selecting the most suitable method and identifying key regional factors remains a challenging task. To address this, this study assessed the performance of six machine learning models, including Convolutional Neural Networks (CNNs), Random Forest (RF), Categorical Boosting (CatBoost), their CNN-based hybrid models (CNN+RF and CNN+CatBoost), and a Stacking Ensemble (SE) combining CNN, RF, and CatBoost in mapping landslide susceptibility along the Karakoram Highway in northern Pakistan. Twelve geospatial factors were examined, categorized into Topography/Geomorphology, Land Cover/Vegetation, Geology, Hydrology, and Anthropogenic Influence. A detailed landslide inventory of 272 occurrences was compiled to train the models. The proposed stacking ensemble and hybrid models improve landslide susceptibility modeling, with the stacking ensemble achieving an AUC of 0.91. Hybrid modeling enhances accuracy, with CNN–RF boosting RF’s AUC from 0.85 to 0.89 and CNN–CatBoost increasing CatBoost’s AUC from 0.87 to 0.90. Chi-square (χ²) values (9.8–21.2) and p-values (<0.005) confirm statistical significance across models. This study identifies approximately 20.70% of the area as from high to very high risk, with the SE model excelling in detecting high-risk zones. Key factors influencing landslide susceptibility showed slight variations across the models, while multicollinearity among variables remained minimal. The proposed modeling approach reduces uncertainties, enhances prediction accuracy, and supports decision-makers in implementing effective landslide mitigation strategies.

1. Introduction

Landslides are among the most frequent and destructive natural hazards, particularly in mountainous regions, where they pose significant risks to human life, ecosystems, and infrastructure and even lead to the abandonment of towns [1,2]. These events involve the downward movement of debris, soil, and rocks due to gravity, classified by material type (e.g., mud, rock, soil, debris) and movement mode (toppling, flowing, sliding) [3]. The causes of landslides are multifactorial, resulting from natural triggers such as tectonic activity, heavy rainfall, snowmelt, and earthquakes, as well as human-induced factors like deforestation, urbanization, and infrastructure development [4,5]. Over the past few decades, anthropogenic activities, particularly in developing countries, have exacerbated the frequency and severity of landslides. Rapid population growth and unregulated development in mountainous regions have increased vulnerabilities, resulting in significant human casualties and economic losses [6,7,8]. The growing threat requires the urgent need for effective risk management strategies to identify hazardous areas and mitigate landslide impacts [9,10].

The Karakoram Highway (KKH), a critical transport route connecting Pakistan and China within the China–Pakistan Economic Corridor (CPEC), is highly vulnerable to landslides due to its steep topography, seismic activity, and extreme climatic conditions. Constructed between 1974 and 1978 and operational since 1979, the KKH traverses rugged, high-altitude terrain characterized by loose debris and frequent geological hazards such as landslides, rockfalls, and avalanches [11,12]. Past events, such as the 2010 Attabad landslide that blocked the KKH and created a massive dam, illustrate the catastrophic consequences of these hazards [13,14]. Landslide risk assessment along the Karakoram Highway (KKH) is particularly challenging due to harsh environmental conditions, limited data availability, and technical complexities.

Recent advances in Artificial Intelligence (AI), particularly in Machine Learning (ML), have significantly improved landslide susceptibility modeling, offering more accurate tools to identify landslide risk zones and support effective disaster mitigation [15,16,17]. For instance, Liang et al. [18] compared traditional models like Information Value (IV) and Logistic Regression (LR) with advanced methods such as CatBoost and Convolutional Neural Networks (CNNs), demonstrating the superiority of ML approaches in Tibet. Seydi et al. [19] evaluated ML models like Support Vector Machines (SVMs), Random Forest (RF), and Cascade Forest Models (CFMs) for flood-prone zones, showcasing their versatility. Ensemble models, such as stacking techniques, further enhance accuracy; Yousefi et al. [20] used optimized SVMs, RF, and XGBoost for high-accuracy landslide maps in Kermanshah, Iran. Qin et al. [21] integrated SVMs, BPNN, and RF optimized with the Sparrow Search Algorithm (SSA) to improve susceptibility predictions. In deep learning, LeCun, Bengio, and Hinton [22] advanced image recognition in remote sensing with DNNs, while Hochreiter and Schmidhuber [23] developed LSTM-based RNNs for analyzing environmental data, such as Himalayan weather patterns. Park et al. [24] studied landslide susceptibility mapping in Korea’s Jumunjin area using decision tree models, and Rumelhart, Hinton, and Williams [25] demonstrated Multi-Layer Perceptrons (MLPs) effectiveness in processing geospatial data.

To address the limitations of traditional landslide susceptibility assessments, which often rely on conventional statistical models that may not fully capture the complex interactions between geo-environmental factors [26], recent research has focused on AI-driven hybrid modeling to improve landslide susceptibility mapping (LSM). AI-powered machine learning classifiers such as support vector machines (SVMs) and random forests (RFs) have demonstrated significant effectiveness in LSM tasks [27,28], offering robust capabilities in feature selection, dimensionality reduction, and handling complex datasets with high accuracy and interpretability [29]. Convolutional neural networks (CNNs) offer an alternative, automating feature extraction to enhance LSM model performance [30]. In hybrid LSM models, CNNs are used for feature extraction before classification by ML algorithms, combining the strengths of both methods and boosting model robustness [31,32].

This study investigates and compares machine learning techniques for assessing landslide susceptibility along the Karakoram Highway. Its application to LSM provides advancements by integrating AI-driven methodologies, including Convolutional Neural Networks (CNNs), Random Forest (RF), and Categorical Boosting (CatBoost)—a combination not previously explored. A range of machine learning models—including traditional, hybrid, and ensemble methods—are used to analyze terrain vulnerabilities. Twelve region-specific conditioning factors, such as slope, aspect, lithology, curvature, and land use, are incorporated into the analysis. Advanced feature selection techniques, including Geodetector and Information Gain Ratio (IGR), along with Variance Inflation Factor and correlation analyses, are employed to identify and refine the most influential factors. The findings provide valuable insights into landslide susceptibility, contributing to efforts aimed at addressing landslide hazards along the highway.

2. Materials and Methods

2.1. Study Area and Geological Setting

The research area, strategically located along the Karakoram Highway (KKH) in northern Pakistan, spans a 426 km stretch with a 10 km radius buffer zone, forming a critical part of the China–Pakistan Economic Corridor (CPEC), as shown in Figure 1 [33]. Despite its strategic importance, this region is highly vulnerable to hydro-climatological and geological hazards, particularly landslides, due to extreme weather conditions, active tectonics, and high erosion rates. The area experiences harsh winters, moderate summers, annual rainfall of 120–130 mm, and significant temperature fluctuations ranging from −21 °C to 25 °C [34]. Geologically, the terrain comprises diverse formations spanning the Paleozoic, Mesozoic, and Cenozoic eras, shaped by significant tectonic and magmatic processes. The Chilas Complex, near Chilas, consists of Cretaceous mafic-ultramafic rocks (K) within the Kohistan Arc [35,36], while the Kohistan Batholiths contain Tertiary igneous rocks (Ti), as depicted in Figure 2, including granodiorite, granite, and diorite from Eocene–Oligocene magmatism. Metasedimentary rocks, such as Lower Paleozoic rocks (Pzl) and Jurassic metamorphic–sedimentary rocks (Jms), reflect deposition and Himalayan orogenic metamorphism, while the Kamila Amphibolite Complex marks Jurassic arc magmatism [37,38,39,40]. Quaternary deposits (Qs), formed through glacial and fluvial processes, contribute to loosening soils and steepening slopes, heightening landslide susceptibility [41]. Key settlements within this dynamic landscape, including Chilas, Gilgit, and surrounding areas, require continuous monitoring to mitigate geological risks.

2.2. Landslide Inventory Mapping

Landslide inventory mapping is a critical step in landslide susceptibility analysis, as shown in Figure 3, with the accuracy of the inventory map significantly influencing model performance [42]. Detailed inventory maps, including information on landslide type, extent, and location, provide essential data for predicting future landslides [43]. In this study, Landslide sites were identified using high-resolution Sentinel-2 satellite imagery (10 m), provided by the European Space Agency (ESA) and Google Earth images, cross-referenced with literature references [38,44,45,46,47], NASA geohazard records, and field investigations to ensure data reliability and accuracy (Figure 4). The inventory comprises 272 polygons representing landslides in the study area, categorized into translational and rotational types. Most landslides were very small (74.65%, <1 × 10⁵ m²), followed by small (14.70%, 1 × 10⁵–3 × 10⁵ m²), medium (8.88%, 3 × 10⁵–5 × 10⁵ m²), large (1.55%, 5 × 10⁵–1 × 10⁶ m²), and very large (0.22%, >1 × 10⁶ m²). The majority of landslides had volumes between 10⁴ m³ and 10⁶ m³. Using ESRI ArcGIS, an equal number of non-landslide polygons were randomly selected based on the methodology of Yu et al. [48], assuming these areas were free from landslides. However, this approach may risk contamination, as some landslide-prone areas may be mistakenly included in the negative samples. Future studies will investigate more robust sampling methods, such as those proposed by Hu et al. [49] and Khabiri et al. [50], to improve the accuracy of negative sample selection. Furthermore, for model development, 70% of the data were used for training, while 30% were reserved for validation.

2.3. Conditioning Factors

Landslide occurrences are influenced by various factors that directly or indirectly contribute to their triggering mechanisms [51]. The criteria for evaluating landslide susceptibility mapping depend on the unique characteristics of the study area, landslide types, primary causes, available data, and chosen assessment methods [52]. Although no universally accepted guidelines exist for selecting conditioning factors [53], this study identifies twelve factors based on a comprehensive review of relevant literature and region-specific data, as shown in Figure 5. These factors are grouped into topographic (elevation, slope, Topographic Wetness Index, curvature, aspect), environmental (NDVI, rainfall, land use), hydrological (proximity to streams), geological (lithology, proximity to faults), and anthropogenic (proximity to roads) categories.

Elevation, slope, and aspect significantly influence landslide susceptibility, particularly in mountainous regions where they affect rainfall dynamics and terrain stability [54]. Lithology and proximity to faults also play crucial roles, with less permeable or fractured rock increasing landslide risks [55]. Additionally, the Normalized Difference Vegetation Index (NDVI) is vital for assessing vegetation health, which impacts water infiltration and slope stability, making it a key factor in landslide studies [56]. The Topographic Wetness Index (TWI) identifies regions susceptible to water accumulation and soil saturation [57]. It is calculated using the following equation [58]:

$$TWI=\mathrm{l}\mathrm{n}\left({\displaystyle \frac{A_s}{\mathrm{t}\mathrm{a}\mathrm{n}\beta }}\right)$$

(1)

where $A_s$ represents the specific catchment area, and $\mathrm{t}\mathrm{a}\mathrm{n}\beta$ denotes the slope angle.

Slope curvature affects water dynamics, with convex slopes facilitating runoff and concave slopes retaining water. Increased rainfall heightens landslide risk by enhancing soil moisture [59]. Changes in land use impact soil stability, while proximity to streams and roads exacerbates landslide susceptibility due to erosion, flooding, and structural disturbances from traffic and construction [60].

The Alaska Satellite Facility (ASF), located at the University of Alaska Fairbanks, Fairbanks, AK, USA, offers the ALOS-PALSAR Digital Elevation Model (DEM) with a 12.5 m resolution, accessible at https://search.asf.alaska.edu/ (accessed on 7 July 2024), which is essential for analyzing landslide susceptibility by deriving factors such as slope, aspect, curvature, the Topographic Wetness Index (TWI), and elevation, etc., as shown in

Table 1. Lists the conditioning factors, sources, and spatial resolutions used for landslide susceptibility assessment.

Sl. No.	Factors	Source	Spatial Detail
1	Slope	DEM	12.5 m
2	Aspect	DEM	12.5 m
3	Curvature	DEM	12.5 m
4	Landcover	USGS Global Land Cover (https://www.usgs.gov/centers/eros/science) (accessed on 7 July 2024)	12.5 m (resampled)
5	Lithology	National Library of Australia (http://nla.gov.au/nla.obj-233602228) (accessed on 7 July 2024)	12.5 m (resampled)
6	Proximity to fault	National Library of Australia (http://nla.gov.au/nla.obj-233602228) (accessed on 7 July 2024)	12.5 m (resampled)
7	Proximity to road	Natural Earth Data (https://www.naturalearthdata.com/downloads/10m-cultural-vectors/) (accessed on 7 July 2024)	12.5 m (resampled)
8	Topographic Wetness Index (TWI)	DEM	12.5 m
9	Normalized Difference Vegetation Index (NDVI)	(https://apps.sentinel-hub.com/eo-browser/) (accessed on 7 July 2024)	12.5 m (resampled)
10	Elevation	DEM	12.5 m
11	Rainfall	Climate Hazards Group InfraRed Precipitation with Station data	12.5 m (resampled)
12	Proximity to stream	Natural Earth Data (https://www.naturalearthdata.com/downloads/10m-cultural-vectors/) (accessed on 7 July 2024)	12.5 m (resampled)

. Additional parameters were incorporated using resampled data from various sources, including the USGS for global land cover, the National Library of Australia for lithology and proximity to faults, Natural Earth Data for proximity to roads and proximity to streams, the Climate Hazards Group for rainfall data, and Sentinel-2 imagery for NDVI analysis. Thematic maps for these conditioning factors were created using ArcGIS software version 10.8. To ensure consistency, all datasets were resampled to maintain the same 12.5 m resolution as the DEM. The data were standardized and normalized to streamline the analysis, with each factor categorized into different levels of influence. All maps were aligned to the WGS 1984 datum and projected in the UTM zone 42 coordinate system.

2.4. Machine Learning Models

2.4.1. Convolutional Neural Network

Convolutional Neural Networks (CNNs), a powerful Artificial Intelligence (AI) tool, are highly effective for spatial data processing, making them suitable for landslide susceptibility mapping. This study utilized TensorFlow, an open-source machine learning framework developed and maintained by Google LLC (Mountain View, CA, USA). The convolutional neural network (CNN) architecture employed consisted of three 3 × 3 convolutional layers with 32, 64, and 128 filters, each followed by 2 × 2 max-pooling layers. Two dense layers (128 and 64 neurons) employed ReLU activation, while the output layer used a sigmoid function for binary classification. AI-powered data augmentation techniques, including random rotations and scaling, enhanced model robustness and mitigated overfitting, enabling accurate landslide prediction, as shown in the following equation [61]:

$$C_k={\sum }_{i=1}^{M}f\left(w_k\ast x_i+b_k\right),k=\mathrm{1,2},\dots ,K$$

(2)

where the output $C_k$ of the convolutional operation for the $k$-th kernel sums all $M$ input vectors, where $f$ is the activation function, $w_k$ is the weight matrix, $x_i$ is the $i$-th input vector, $b_k$ is the bias term, and $\ast$ denotes convolution.

2.4.2. Random Forest

Random Forest (RF), an AI-driven ensemble learning method, is effective for landslide susceptibility mapping by combining predictions from multiple decision trees to improve accuracy and robustness. Random Forest (RF), implemented using Scikit-learn, an open-source library developed by the Scikit-learn team and maintained by INRIA, France, was applied to spatial and geological datasets with a 70%–30% train–test split to ensure balanced evaluation. An AI-based random sampling of features and data subsets during training reduced overfitting and enhanced generalization. Final predictions were averaged across all trees, providing stable and reliable results, as shown in the following equation [62]:

$$G(x)={\displaystyle \frac{1}{N}}{\sum }_{i=1}N{g}_i(x)$$

(3)

In this context, $G(x)$ represents the final aggregated prediction for the input $x$, $N$ is the total number of models in the ensemble, and $g_i(x)$ is the prediction from the $i$-th model.

2.4.3. Categorical Boosting

Categorical Boosting (CatBoost), an advanced Artificial Intelligence (AI) based gradient boosting algorithm developed by Yandex (headquartered in Moscow, Russia), efficiently handles categorical features through ordered boosting, optimizing predictive performance. CatBoost was implemented using its dedicated library and applied to spatial and categorical datasets for landslide susceptibility mapping. AI-optimized hyperparameter tuning, including adjustments to the number of trees, learning rate, and tree depth, along with L2 regularization and bagging techniques, was employed to reduce overfitting. The dataset was split into training (70%) and testing (30%) subsets to ensure balanced evaluation and robust, reliable predictions, as shown in Equation (4) [63].

$$F(x)=F_0(x)+\sum _{m=1}M{\alpha }_m{T}_m(x)$$

(4)

where $F(x)$ is the final prediction for input $x$, starting from $F_0(x)$, the initial prediction. The final result sums the weighted predictions from $M$ trees: ${\sum }_{m=1}^M{\alpha }_m{T}_m\left(x\right).$

2.4.4. CNN-Based Hybrid Modeling

This hybrid modeling approach for Landslide Susceptibility Mapping (LSM) utilizes Artificial Intelligence (AI) by combining Convolutional Neural Networks (CNNs) for advanced feature extraction with traditional machine learning techniques for robust classification. As shown in Figure 6, the process includes two stages: feature extraction and classification [64,65]. In the feature extraction stage (Figure 6a), spatial and geological data are preprocessed for analysis and input into a Convolutional Neural Network (CNN) using architectures like VGG16 or ResNet. Using AI-driven transfer learning, CNN extracts high-level spatial features, transforming raw data into compact, informative feature vectors. In the classification stage (Figure 6b), the refined feature vectors are passed into machine learning classifiers such as Random Forest (RF), Gradient Boosting, and CatBoost. These models perform binary classification to distinguish between landslide-prone and non-prone areas, with output probabilities indicating the susceptibility of each region. The final susceptibility predictions are integrated into a Geographic Information System (GIS) to generate accurate, high-resolution landslide susceptibility maps.

2.4.5. Stacking Ensemble

Stacking Ensemble is an advanced Artificial Intelligence (AI) based method that combines multiple base models to enhance predictive performance. Unlike bagging and boosting, it employs a two-tier framework where a meta-learner analyzes base model predictions to generate the final output, enabling hierarchical learning. The process involves two key stages, as follows.

(1)

Training Base Models

In this study, the base models selected were Convolutional Neural Networks (CNNs) for automatic feature extraction, Random Forest (RF) for handling high-dimensional data, and Categorical Boosting (CatBoost) for processing categorical variables. The stacking architecture was constructed by splitting the dataset into two disjoint subsets: 70% for training the base models and 30% for generating out-of-fold (OOF) predictions through K-fold cross-validation. K-fold cross-validation was employed for each base model, dividing the dataset into $k$ subsets. Each base model was trained on $k-1$ folds, with the remaining fold used for validation, ensuring that every data point was utilized for both training and validation. This process enabled independent predictions on unseen data, reducing overfitting and mitigating data leakage. The OOF predictions, derived from the validation subsets, served as input features for training the meta-learner in the second-level stacking phase. These predictions from all base models were combined to form the second-level dataset, where each data point comprises the base model predictions on the validation subset, which were then used to train the meta-learner.

(2)

Training Meta-Learner

The second-level dataset was constructed by collecting AI-generated OOF predictions from all base models for the validation subset. Alongside these predictions, the true labels $\hat{y}$ from the validation set were included. This dataset was represented as a matrix $Z$, where each element $h_i(x_j)$ denotes the prediction of the $i$-th base model for the $j$-th data point. The matrix $Z$ is structured as follows:

The second-level dataset was constructed by collecting the out-of-fold (OOF) predictions from all base models for the validation subset. Alongside these predictions, the true labels $\hat{y}$ from the validation set were included. This dataset was represented as a matrix $Z$, where each element $h_i(x_j)$ denotes the prediction of the $i$-th base model for the $j$-th data point. The matrix $Z$ is structured as follows:

$$Z=\left[\begin{array}{cccc}{h}_1\left(x_1\right)& h_1\left(x_2\right)& \cdots & h_1\left(x_n\right)\\ h_2\left(x_1\right)& h_2\left(x_2\right)& \cdots & h_2\left(x_n\right)\\ ⋮& ⋮& \ddots & ⋮\\ h_M\left(x_1\right)& h_M\left(x_2\right)& \cdots & h_M\left(x_n\right)\end{array}\right]$$

(5)

where $n$ represents the number of data points in the validation subset, and $M$ is the number of base models. Each row corresponds to predictions for a single data point across all base models, while each column represents predictions of a specific base model for all data points.

For the second-level model, Logistic Regression (LR) was chosen as the meta-learner due to its simplicity, effectiveness in combining model outputs, and recommendations from previous studies [66,67]. The second-level dataset $Z$, along with the true labels $\hat{y}$, was used to train the meta-learner. To ensure robust evaluation, the dataset was split into 70% for training and 30% for validation. Hyperparameter tuning was not performed, as LR typically performs well with its default settings. After training, the meta-learner combines the predictions from the base models to make final predictions for new, unseen data points. For a given input $x$, the final prediction $\hat{y}$ is computed using the following equation [68].

$$\hat{y}=LR(h^1(x),h^2(x),\dots ,h^M(x))$$

(6)

where the final prediction $\hat{y}$ by the stacking ensemble model is made by the meta-learner $LR$ (Logistic Regression) which combines the base models’ predictions $(h^1(x),h^2(x),\dots ,h^M(x))$ for input $x$.

2.5. Multicollinearity and Feature Selection

To improve the accuracy of landslide susceptibility modeling, we addressed multicollinearity among factors such as slope, elevation, and rainfall by analyzing each input variable independently. The Variance Inflation Factor (VIF) and Tolerance (TOL) were employed to assess multicollinearity among the 12 input variables. VIF measures how much the variance of regression coefficients increases due to correlations between predictors, while its inverse, TOL, shows the level of independence among predictors. Typically, VIF values > 10 and TOL values < 0.1 signal significant multicollinearity issues. The formulas for VIF and TOL are as follows [69]:

$$VI{F}_j={\displaystyle \frac{1}{1-R_j^2}}={\displaystyle \frac{1}{TOL}}(j=\mathrm{1,2},3,\dots ,n)$$

(7)

where $R_j^2$ is the determination coefficient for the $j\mathrm{th}$ variable.

Feature selection in machine learning enhances performance, accelerates computation, improves interpretability, and simplifies the handling of high-dimensional data. In this study, the Information Gain Ratio (IGR) and GeoDetector methods were applied to identify the most relevant features. IGR, an improvement on the Information Gain method, evaluates feature importance by accounting for both information gain and dataset entropy. Higher IGR values indicate more significant factors for predicting landslide susceptibility, calculated as follows [70]:

$$\mathrm{I}\mathrm{G}\mathrm{R}\left(x\right)={\displaystyle \frac{F\left(N\right)}{F\left(N\right)·\mathrm{I}\mathrm{G}\left(Y,x\right)}}={\displaystyle \frac{F\left(N\right)-F\left(N|x\right)}{F\left(N\right)}}$$

(8)

where $F\left(N\right)$ is the empirical entropy of the dataset $N$, $F\left(N|x\right)$ is the empirical conditional entropy of factor $x$ given the dataset $N$, and $\mathrm{I}\mathrm{G}\left(Y,x\right)=F\left(N\right)-F\left(N|x\right)$ is the Information Gain of factor $x$.

GeoDetector (GD) is a method to detect spatial differentiation and identify driving factors, such as landslides. It measures how well a factor $x$ explains the spatial pattern of landslides $y$ using the $q$ value. A higher $q$ value indicates a more important factor and can be represented by the following equation [71]:

$$q=1-{\displaystyle \frac{N{\sigma }^2}{{\sum }_{x=1}^L{N}_x{\sigma }_x^2}}$$

(9)

where $L$ is the number of subcategories, $N_x$ and $N$ are the unit counts in the $x$-th category and overall, respectively, and ${\sigma }_x^2$ and ${\sigma }^2$ are the variances of $y$ in the $x$-th category and overall, respectively; values range from 0 to 1, with higher values indicating stronger explanatory power.

The Weight of Evidence (WoE) approach, based on Bayesian statistical theory, examines the interaction between spatial characteristics and key factors influencing landslides. Using a logarithmic model, it assesses both positive ($\gamma$) and negative ($\delta$) associations and combines these to determine their cumulative effect ($\theta$) on landslide probability as illustrated below [72].

$$\gamma =\ln \left[P\left(A|B\right)\cdot P\left(A|B^{\prime }\right)\right]$$

(10)

$$\delta =\ln \left[P\left(A^{\prime }|B\right)\cdot P\left(A|B^{\prime }\right)\right]$$

(11)

$$\theta =\gamma +\delta$$

(12)

where $P$ is the event probability, $A$ represents influential factors, $A^{\prime }$ represents non-factors, $B$ denotes landslide-prone areas, and $B^{\prime }$ denotes areas less prone to landslides.

2.6. Hyperparameter Optimization

The models used for landslide susceptibility mapping were optimized using the grid search method, ensuring robust training and evaluation by selecting the most effective hyperparameter combinations outlined in

Table 2. Hyperparameter optimization values for machine learning models.

Model	Hyperparameter Values
CNN	Learning rate = 0.01, epochs = 100, batch size = 500
RF	Number of trees = 800, mtry = 5, max depth = 10, min samples leaf = 5
CatBoost	Number of trees = 500, learning rate = 0.01, depth = 6, min samples leaf = 5

. For the Convolutional Neural Networks (CNNs), the number of epochs was set to 100, indicating that the entire dataset was passed through the network 100 times during training to ensure sufficient iterations for convergence [73], while the learning rate was explored in the range of from 0.0001 to 0.01 with a step size of 0.001, and the batch size was tested at 32, 64, and 128 to balance training time and model accuracy. For the Random Forest (RF) model, 800 trees were used, with the mtry (maximum number of features) explored in the range of from 2 to 10 with a step size of 1, the maximum depth tested from 5 to 20, and the minimum samples per leaf varied from 1 to 10, with a step size of 1, to optimize model performance. The CatBoost model was configured with 500 trees, a learning rate explored in the range of from 0.001 to 0.1 with a step size of 0.005, a depth tested from 4 to 10 with a step size of 1, and a minimum of five samples per leaf, to optimize model performance [74,75]. The number of trees in both models reflects their complexity and capacity to train effectively on the dataset, with hyperparameters chosen to balance model complexity and prevent overfitting, as shown in Figure 7.

2.7. Model Validation

The evaluation of landslide susceptibility models utilized three key metrics. The Receiver Operating Characteristic (ROC) curve, measured by the Area Under the Curve (AUC), evaluates classification effectiveness, with values close to 1 indicating high accuracy and values near or below 0.5 suggesting poor performance [76]. The F-Score, the harmonic mean of precision and recall, ranges from 0 to 1, with higher values reflecting better performance [77]. The True Skill Statistic (TSS), ranging from −1 to +1, compares model predictions against random and ideal results, where +1 denotes perfect accuracy and 0 or less indicates no better than random guessing [78]. These metrics were applied to a test dataset comprising 30% of the total data for model validation.

Additionally, the McNemar test was employed to assess statistically significant differences in binary outcomes between methods [79]. A chi-square value > 3.841 and a p-value < 0.05 indicate a significant difference, highlighting variation in performance between the methods [80].

3. Results

3.1. Evaluation of Influencing Variables

In this study, multicollinearity among the conditioning factors was analyzed, with results presented in

Table 3. Multicollinearity analysis for variables contributing to landslides.

Variable	VIF	TOL
Rainfall	3.98	0.25
Lithology	3.63	0.28
Elevation	3.02	0.33
NDVI	2.11	0.48
Slope	1.59	0.63
TWI	1.54	0.65
Proximity to Road	1.52	0.67
Proximity to Stream	1.51	0.68
Landcover	1.41	0.71
Proximity to fault	1.08	0.92
Curvature	1.04	0.96
Aspect	1.03	0.98

. All factors exhibited low multicollinearity, as indicated by VIF values below 4, confirming their independent contributions to the models. Aspect showed the lowest VIF (1.03) and highest TOL (0.98), while Rainfall exhibited the highest VIF (3.98) and lowest TOL (0.25). These results demonstrate that none of the twelve factors exhibited significant collinearity, validating their suitability as input variables for the models.

As illustrated in Figure 8, the Information Gain Ratio and GeoDetector methods were used to identify key factors influencing landslide susceptibility. Based on these analyses, Proximity to Roads, Slope, and Elevation are the most influential factors for landslide susceptibility, with Proximity to Roads showing the highest impact (q = 0.211, IGR = 0.069), followed by Slope (q = 0.077, IGR = 0.046) and Elevation (q = 0.138, IGR = 0.028). Other factors, like Rainfall, Land Use, Proximity to Streams, and Topographic Wetness Index, had moderate impacts with comparatively lower values, such as Land Use (q = 0.066, IGR = 0.016) and Rainfall (q = 0.037, IGR = 0.024). The least influential factors were Aspect (q = 0.009, IGR = 0.004) and Curvature (q = 0.007, IGR = 0.003), indicating minimal effects on susceptibility. The analysis shows that, despite varying levels of importance, each factor enhances the accuracy of landslide prediction models, supporting their inclusion in further analysis.

3.2. Landslide Susceptibility Assessment

In mapping landslide susceptibility, various predictive models are used to predict landslide occurrence probabilities. The study area is classified into five levels of landslide susceptibility, from very low to very high, using the quantile method. This approach uses selected landslide-contributing factors as input variables, with the model output classifying each grid unit as either landslide (1) or non-landslide (0), closely reflecting the region’s geographical characteristics.

Table 4. Classification of susceptibility and corresponding landslide percentages determined by machine learning models.

RF
Susceptibility Level	Pixel Count	Area Ratio (%)	Landslide Ratio (%)
Very Low	3,487,332	41.67	2.8
Low	2,314,819	27.66	14.95
Moderate	1,194,213	14.27	15.89
High	945,971	11.3	35.51
Very High	426,181	5.09	30.84
CatBoost
Susceptibility level	Pixel count	Area ratio (%)	Landslide ratio (%)
Very Low	1,300,925	15.55	1.87
Low	4,984,946	59.57	21.5
Moderate	1,318,324	15.75	24.3
High	331,589	3.96	16.82
Very High	432,732	5.17	35.51
CNN
Susceptibility level	Pixel count	Area ratio (%)	Landslide ratio (%)
Very Low	1,633,416	19.52	0.93
Low	3,021,970	36.11	9.35
Moderate	1,979,029	23.65	18.69
High	1,274,154	15.23	30.84
Very High	459,947	5.5	40.19
CNN–RF
Susceptibility level	Pixel count	Area ratio (%)	Landslide ratio (%)
Very Low	1,712,661	20.47	2.47
Low	2,500,588	29.88	9.88
Moderate	2,107,671	25.19	17.28
High	1,325,556	15.84	43.21
Very High	722,040	8.63	59.26
CNN–CatBoost
Susceptibility level	Pixel count	Area ratio (%)	Landslide ratio (%)
Very Low	2,130,078	25.45	1.87
Low	2,565,896	30.66	6.54
Moderate	1,878,321	22.45	20.56
High	1,124,562	13.44	29.91
Very High	669,659	8	41.12
SE
Susceptibility level	Pixel count	Area ratio (%)	Landslide ratio (%)
Very Low	1,300,925	15.55	1.87
Low	2,198,087	26.27	4.67
Moderate	2,189,121	26.16	14.95
High	1,485,091	17.75	22.43
Very High	1,195,292	14.28	56.07

shows that the Very Low susceptibility category consistently covers a significant portion of the study area across all models, ranging from 15.55% to 41.67% of the region, with landslide occurrence rates between 0.93% and 2.8%. In contrast, the Low susceptibility category shows significant spatial variation, occupying from 26.27% to 59.57% of the area, with landslide frequencies from 4.67% to 21.5%. The Moderate susceptibility zones cover from 14.27% to 26.16% of the area, experiencing landslide frequencies of from 14.95% to 24.3%. High susceptibility zones, with some variation, span from 3.96% to 17.75% of the area and exhibit landslide frequencies between 16.82% and 43.21%. Lastly, the Very High susceptibility category, though covering the smallest area (from 5.09% to 14.28%), consistently shows the highest landslide frequencies across models, ranging from 30.84% to 59.26%.

Additionally, the analysis, as shown in Figure 9, highlights distinct trends in risk categorization across models, stemming from inherent differences in their algorithms, feature sensitivities, and decision-making processes. The Random Forest model generally classifies larger portions of the area as very low risk, reflecting its tendency to generalize more broadly, while CatBoost designates smaller areas within this category due to its ability to handle categorical variables effectively and provide more granular classifications. For low-risk zones, CatBoost covers the most extensive area, whereas CNN–RF identifies a relatively smaller portion. In the moderate-risk category, variability among models is notable, with CNN–RF assigning the largest area and Random Forest covering a smaller extent. Hybrid models like CNN–RF and CNN excel at identifying spatial patterns, which may lead to broader high-risk designations in areas with strong feature correlations. For high-risk zones, CNN–RF and the Convolutional Neural Network designate substantial areas, indicating higher landslide risks, while CatBoost identifies a smaller region. In the very high-risk category, CNN–RF allocates the largest area, corresponding to the highest observed landslide frequency, while CatBoost and CNN identify smaller high-risk regions. The Random Forest and SE models show smaller extents but still indicate elevated landslide frequencies in these zones, further emphasizing the variability and unique strengths of the different modeling approaches. The geographic distribution of landslide susceptibility indicates that Passu consistently ranks as a high-risk area across all models, whereas Chilas and Babusar are classified as comparatively lower-risk zones. High to very high-risk levels are concentrated in the northern regions, particularly near highways and rivers, likely due to rugged terrain and close proximity to these features. In contrast, the central region shows moderate to high-risk levels, indicative of relatively stable geological conditions. The southernmost areas are generally categorized as low to moderate risk, attributed to increased geological stability and greater distance from steep slopes.

3.3. Validation and Model Comparison

The Area Under the Receiver Operating Characteristic Curve (AUC-ROC), True Skill Statistic (TSS), and F-score values derived from the test data demonstrate a satisfactory level of predictive performance for the landslide susceptibility models, as shown in

Table 5. Performance indicators for machine learning models in landslide risk analysis.

Model	AUC Score	F-Score	TSS Metric
CNN	0.84	0.58	0.63
RF	0.85	0.63	0.69
CatBoost	0.87	0.71	0.74
CNN–RF	0.89	0.73	0.80
CNN–CatBoost	0.90	0.83	0.85
Stacking Ensemble	0.91	0.87	0.89

. Notably, the Stacking Ensemble model achieved superior predictive accuracy compared to other models, with an AUC of 0.91, an F-score of 0.87, and a TSS of 0.89. Among the hybrid models, CNN–CatBoost outperforms CNN–RF, showing increases of 0.03 in AUC, 0.12 in F-score, and 0.11 in TSS. Additionally, the hybrid models consistently showed higher accuracy than the standalone CNN across various evaluation metrics. As illustrated in Figure 10, hybrid models clearly outperform standalone classifiers by effectively combining their strengths, resulting in improved predictive performance. For example, the CNN–RF hybrid improves RF’s AUC from 0.85 to 0.89, F-score from 0.63 to 0.73, and TSS from 0.69 to 0.80. Similarly, the CNN–CatBoost hybrid enhances CatBoost’s AUC from 0.87 to 0.90, F-score from 0.71 to 0.83, and TSS from 0.74 to 0.85, demonstrating significant performance improvements. The chi-square (${\chi }^2$) and p-values comparing the models are presented in

Table 6. Outcomes of the McNemar analysis.

Pairwise Comparison	${\mathit{\chi }}^2$	p-Values	Statistical Significance
CNN vs. CNN–CatBoost	7.8	0.02	Yes
CNN vs. CNN–RF	12.5	0.0005	Yes
CNN vs. CatBoost	10.1	0.001	Yes
CNN vs. RF	4.2	0.04	Yes
CNN vs. SE	15.3	0.0001	Yes
CNN–CatBoost vs. SE	19.2	0.00002	Yes
CNN–RF vs. CNN–CatBoost	21.1	0.00001	Yes
CNN–RF vs. SE	17.6	0.00005	Yes
CatBoost vs. CNN–CatBoost	8.4	0.007	Yes
CatBoost vs. CNN–RF	9	0.005	Yes
CatBoost vs. SE	13.3	0.0002	Yes
RF vs. CNN–CatBoost	14.4	0.0001	Yes
RF vs. CNN–RF	5.5	0.04	Yes
RF vs. CatBoost	6.1	0.03	Yes
RF vs. SE	11.7	0.002	Yes

. All pairwise model comparisons reveal statistically significant differences, indicating that the observed performance variations are not due to chance but instead reflect meaningful distinctions in the models’ underlying structures and working processes. The ${\chi }^2$ values for each comparison range from 4.2 to 21.1, consistently surpassing the critical threshold of 3.841. Similarly, p-values range from 0.04 to 0.00001, all below the conventional significance level of 0.05, thereby confirming the statistical significance of the observed differences. These results demonstrate that the models utilize distinct approaches to pattern recognition and decision-making, leading to significant variations in performance.

3.4. Analysis of Key Factors in Predictive Models

The spatial patterns of landslide susceptibility are significantly influenced by the characteristics of conditioning variables, making their importance a crucial aspect of assessment within machine learning models. In this study, we used the information gain method to identify the primary contributing factors for each model. For the CNN model, the top contributing factors were proximity to roads (26%), rainfall (22%), and slope (15%). In contrast, the RF model ranked slope (25%), elevation (22%), and proximity to roads (19%) as the most influential factors. For the CatBoost model, proximity to faults (22%), slope (20%), and proximity to roads (18%) were identified as the highest-priority factors. These findings align with previous studies. For example, Kasahara et al. [81] found that, in the Güneysu district of Rize, tea farming areas on steep slopes (from 30° to 40°) have a from 3.5 to 9.1 times higher likelihood of landslides compared to forested regions. Additionally, Ye et al. [82] used the CatBoost algorithm to evaluate the importance of conditioning factors in Fujian Province, China, where distance to faults was identified as the most influential factor in landslide occurrence. As shown in Figure 11, the most significant factors for the CNN, CatBoost, and RF models are proximity to roads, elevation, slope, and distance to faults. Conversely, curvature, land cover, and aspect were identified as the least important factors across all models. This is consistent with findings by Youssef and Pourghasemi [83], who found minimal impact of profile curvature and land use/land cover on landslide susceptibility in the Abha basin, Saudi Arabia. Similarly, Kavzoglu and Teke [84], in their research in the Macka district of Trabzon Province, Turkey, identified elevation and slope as the most significant factors, while plan curvature and NDVI were the least significant.

3.5. Factors Impacting Landslide Susceptibility

The influence of various sub-classes of conditioning variables on landslide occurrence is illustrated in Figure 12. Elevation significantly influences landslide occurrences, with 17.8% of landslides occurring within the 2201–3450 m range, with a Weight of Evidence (WoE) of 0.88, while the 4701–5950 m range shows lower susceptibility with only 10.2% of landslides. Slope steepness also affects landslide risk, with slopes between 32.71° and 48.88° accounting for 13% of landslides in the study area. Proximity to geological fault zones, roads, and streams within 500 m further increases susceptibility, showing WoE values of 0.5, 0.53, and 0.4, and landslide rates of 17.8%, 19.7%, and 15%, respectively. Land cover, aspect, TWI, and curvature further contribute to landslide susceptibility. Snow-covered regions, especially those prone to freeze–thaw cycles, show a higher landslide risk, with a WoE value of 1.6 and accounting for 19.4% of landslides, while areas with waterbodies exhibit minimal risk. Aspect is another contributing factor, with south-facing slopes (289°–360°) and west-facing slopes (0°–72°) being more prone to landslides, showing WoE values of 0.3 and 0.2, and landslide rates of 10.5% and 10.25%, respectively. Moreover, high TWI values (21.4–26.3) correspond to a 16.3% landslide rate. Curvature shows an influence, with steeper areas ranging from 23.4° to 35.2° accounting for 12.1% of landslides, while flatter regions show lower susceptibility. Rainfall, NDVI, and lithology influence landslide susceptibility. Areas receiving 427.8–480.4 mm of rainfall annually show a higher landslide rate of 15.5%, whereas regions with lower rainfall (322.4–375.0 mm/year) have a reduced rate of 7%. Sparse vegetation (NDVI from −0.29 to −0.15) correlates with a higher landslide rate of 12%, while dense vegetation (NDVI from 0.65 to 0.85) is associated with a lower rate of 6.4%. Lithology also plays a role, with limestone (Ti) formations showing a 12% landslide rate, whereas more stable formations, such as Precambrian rocks (pC), exhibit lower susceptibility.

4. Discussion

4.1. Assessing Performance of Models

This study offers significant insights into understanding and analyzing landslide susceptibility. The analysis demonstrated that the Artificial Intelligence (AI) driven stacking ensemble and hybrid model were the most effective, outperforming other approaches such as random forest, CatBoost, and CNN. The superior performance of the AI-powered stacking ensemble lies in its ability to integrate predictions from diverse base models through a meta-model, capturing complex patterns and interactions often lacking in analyses performed using individual models or simpler ensembles. This method reduces errors and biases while enhancing robustness and generalization by addressing model-specific weaknesses and capturing the fundamental patterns and relationships within the data. These findings are supported by the work of Lee and Lee [85], which showed that stacking ensembles combining CatBoost, XGBoost, and Random Forest significantly outperformed individual models in predicting landslides using 30 variables across climate, topography, and environment. Studies consistently show that AI-driven ensemble-based approaches improve the reliability of high-resolution landslide susceptibility mapping (LSM) [44], emphasizing their growing importance in the field.

Moreover, CNNs have demonstrated significant effectiveness in landslide susceptibility mapping, with studies such as that of Jiang et al. [86] showing their superiority over traditional methods in Hongxi River Basin, Sichuan Province, China. CatBoost excels in handling categorical data and modeling non-linear relationships, while RF provides robustness through ensemble decision-making. Hybrid AI models like CNN–RF and CNN–CatBoost effectively combine CNN’s spatial feature extraction with the enhanced decision boundaries and overfitting control of RF and CatBoost, resulting in superior predictive accuracy and reliability. In Yongxin County, China, Fang et al. [55] reported performance gains from hybrid models, including a 0.015 AUC increase with CNN–SVM, a 0.027 improvement with CNN–RF, and a 0.063 AUC boost with CNN–LR, which also achieved an 8.72% accuracy increase. Similarly, Aslam et al. [56] observed significant enhancements in Pakistan’s Muzaffarabad and Mansehra districts, where CNN–SVM achieved a 0.01 AUC gain over SVM, CNN–RF improved by 0.027, and CNN–LR achieved the highest gains with a 0.054 AUC boost and an 8.72% accuracy increase, demonstrating the effectiveness of hybrid approaches in diverse regions.

This study demonstrated the effectiveness of an AI-enabled stacking ensemble framework combining CNN, RF, and CatBoost, however, we recognize the potential for further advancements through the incorporation of additional machine learning or deep learning architectures. Emerging AI models, such as Long Short-Term Memory (LSTM) networks and Transformers, could provide unique benefits, particularly for datasets with temporal or sequential characteristics. Although this was beyond the scope of the current study, integrating such models represents a promising direction for future research to evaluate their impact on landslide susceptibility mapping. While focused on the Karakoram Highway, our research demonstrates the generalizability of its approach to landslide susceptibility mapping by offering a transferable framework and comparative insights applicable to other regions. Moreover, the enhanced performance of hybrid models (CNN–CatBoost, CNN–RF) and the Stacking Ensemble highlights their adaptability, as these models can be recalibrated to local conditions, making them a robust tool for diverse geographical contexts.

4.2. Analysis of Factors Contributing to Landslides

Landslide susceptibility in the study area arises from a complex interplay of natural and human-induced factors that contribute to geological instability. Key factors influencing AI-driven susceptibility models include proximity to roads, rainfall, slope, elevation, and proximity to faults. Although some factors may individually have limited impact, their effects are significantly magnified when combined with other geo-environmental or anthropogenic influences. Li et al. [87] highlight that human activities such as road construction, agricultural reclamation, and deforestation exacerbate landslide risk by altering geomorphology and destabilizing slopes. Elevation is particularly critical, as mid-elevation areas face heightened vulnerability due to steep slopes, fractured lithology, and increased human activity, all of which undermine stability. Similarly, Shahi et al. [88] report that a combination of steep topography, lithological diversity, and human interventions significantly amplifies landslide risks, as observed in Nepal’s Jajarkot region. Slope orientation, or aspect, further influences susceptibility, with south- and west-facing slopes being more prone to landslides due to intensified weathering from solar radiation. Bachri et al. [89] demonstrate this relationship in East Java, Indonesia, where slope, geology, and rainfall inform both risk assessments and land-use planning strategies. Proximity to fault zones, roads, and streams also heightens risk by destabilizing slopes through tectonic activity, erosion, and human-induced vibrations, particularly in areas with loose soils, as noted by Wu et al. [90]. Land cover and vegetation play a pivotal role, with sparse vegetation increasing vulnerability due to a lack of root systems to stabilize soil, while dense vegetation provides natural reinforcement. Balocchi et al. [91] confirm this, emphasizing the stabilizing effect of dense vegetation in reducing landslide risk. Lithology is another important factor, as rock types such as limestone are more susceptible to weathering and erosion, while more stable formations, like Precambrian rocks, exhibit greater resistance. Marino et al. [92] show the heightened vulnerability of limestone compared to more resilient rock types. Land surface curvature also affects susceptibility, with concave slopes being generally more vulnerable due to water accumulation, although specific local conditions can elevate risk on convex slopes as well. Zhou et al. [93] note that, while concave slopes typically present higher susceptibility, convex slopes may also become unstable in certain contexts. Despite the recognized importance of soil characteristics, their omission remains a common limitation in landslide studies, including our analysis. Soil properties, such as texture, composition, drainage capacity, and cohesion, are crucial for determining slope stability and the likelihood of landslides, and incorporating detailed soil data in future research will enhance understanding of landslide dynamics. However, this gap is not unique to our research, as Yavuz Ozalp et al. [75], Hussain et al. [94], and Ali et al. [95] also excluded soil type as a variable in their landslide susceptibility mapping studies, despite its potential significance.

4.3. Landslide Prevention and Policy Implications

Landslide Susceptibility Mapping (LSM) is essential for identifying landslide-prone areas and predicting incidents that threaten both environmental stability and socioeconomic development [96]. By identifying zones prone to landslides, this study enables stakeholders—including policymakers, government agencies, land-use planners, environmental organizations, transportation and utility companies, agricultural and forestry managers, and the general public—to make informed decisions that can mitigate landslide impacts. AI-enhanced LSM supports targeted policies and proactive measures to protect vulnerable areas, thereby enhancing environmental and economic resilience. Our analysis of six models demonstrates their effectiveness in identifying vulnerable zones, particularly in northern regions near highways and rivers, where rugged terrain and proximity to these features heighten susceptibility. High- to very high-risk areas require careful planning for agricultural, residential, and infrastructure projects. Community-based disaster risk reduction with real-time monitoring is crucial for timely warnings and effective evacuations. As Que et al. [97] emphasize, public participation in disaster mitigation improves with awareness, trust in authorities, and resource availability, particularly in proactive communities. Coordinated efforts should prioritize data sharing, aligned disaster responses, and investments in training and early warning systems. High-elevation regions, where over 90% of the land consists of middle to high mountain formations, require strategic planning to avoid development in unstable areas, alongside proactive monitoring of surrounding mountains to mitigate landslide risks in nearby valleys. Conversely, for settlements on plains intersected by fault lines, mitigation efforts should focus on robust drainage systems and avoid residential development in areas prone to active geological movement, rather than prioritizing elevation. This research employs an AI-enabled stacking ensemble model utilizing a unique combination of Convolutional Neural Networks (CNNs), Random Forest (RF), and Categorical Boosting (CatBoost), along with the development of a CNN–CatBoost hybrid model, representing a previously unexplored advancement in landslide susceptibility mapping. The AI-driven integration of these models enhances predictive accuracy and robustness, addressing the challenge of model selection in regions with complex geological and environmental factors. This work contributes to natural hazard prediction, offering valuable insights for improving disaster risk management and resilience in landslide-prone areas.

5. Conclusions

Machine learning techniques have become integral to landslide susceptibility mapping (LSM), demonstrating strong predictive capabilities across diverse regions, although a consensus on the optimal methods and influencing factors for specific areas remains lacking. This study evaluates the predictive performance of six machine learning models along the landslide-prone Karakoram Highway in northern Pakistan, using a landslide inventory of 272 occurrences identified through high-resolution Sentinel-2 imagery, Google Earth analysis, and field validation. The majority of landslides (74.65%) were very small (<1 × 10⁵ m²), with very large ones (>1 × 10⁶ m²) being rare (0.22%). Using 12 conditioning factors, such as slope, land use, and lithology, the dataset was divided into 70% for model training and 30% for validation. This research employs a stacking ensemble model that uniquely integrates CNN, RF, and CatBoost alongside a CNN–CatBoost hybrid, representing a previously unexplored advancement for landslide susceptibility assessments. Single-classifier models (CNN, RF, and CatBoost), hybrid models (CNN–RF and CNN–CatBoost), and the stacking ensemble were evaluated, with all models effectively generating landslide susceptibility maps, while the stacking ensemble and hybrid models demonstrated superior performance in F1-scores, TSS, and AUC metrics, providing a robust framework for LSM in the region. Moreover, minimal multicollinearity was observed among the 12 conditioning factors, with rainfall showing the strongest interaction (VIF = 3.98). Key factors influencing landslide susceptibility, including proximity to roads, rainfall, slope, elevation, and proximity to faults, were analyzed using GeoDetector and Information Gain Ratio (IGR) to evaluate variable independence, spatial explanatory power, and their contribution to landslide occurrences. This study highlights critical areas for targeted mitigation, providing a robust framework for integrating susceptibility maps into land use planning and landslide risk reduction, thereby supporting sustainable development in mountainous regions. Despite limitations in dataset scale and scope, future research should incorporate additional factors, higher-resolution data, and diverse geographical contexts to enhance LSM, emphasizing interdisciplinary approaches to address landslide risks effectively.