Coupling Different Machine Learning and Meta-Heuristic Optimization Techniques to Generate the Snow Avalanche Susceptibility Map in the French Alps

Abstract

The focus of this study is to introduce a hybrid predictive framework encompassing different meta-heuristic optimization and machine learning techniques to identify the regions susceptible to snow avalanches. To accomplish this aim, the present research sought to acquire the best-performed model among nine different hybrid scenarios encompassing three different meta-heuristics, namely particle swarm optimization (PSO), gravitational search algorithm (GSA), and Cuckoo Search (CS), and three different ML approaches, i.e., support vector classification (SVC), stochastic gradient boosting (SGB), and k-nearest neighbors (KNN), pertaining to different predictive families. According to diligent analysis performed with regard to the blinded testing set, the PSO-SGB illustrated the most satisfactory predictive performance with an accuracy of 0.815, while the precision and recall were found to be 0.824 and 0.821, respectively. The F1-score of the predictions was found to be 0.821, and the area under the receiver operating curve (AUC) was obtained to be 0.9. Despite attaining similar predictive success via the CS-SGB model, the time-efficiency analysis underscored the PSO-SGB, as the corresponding process consumed considerably less computational time compared to its counterpart. The SHapley Additive exPlanations (SHAP) implementation further informed that slope, elevation, and wind speed are the most contributing attributes to detecting snow avalanche susceptibility in the French Alps.

1. Introduction

Avalanches are natural disasters caused by interactions between the atmosphere, hydrosphere, and biosphere, typically occurring in high- and mid-altitude alpine regions [1]. They can severely impact people, ecosystems, infrastructure, and landscapes in the areas below these mountainous environments. Regarding their evolving mechanism, snow masses suddenly start to descend on steep slopes with the driving force of gravity [2]. Avalanches result in substantial economic losses due to damage to infrastructure, disruption of transportation routes, and the cost of rescue and recovery efforts, but most importantly, numerous fatalities and injuries. Additionally, they can have detrimental effects on the environment and public health, with potential consequences such as water contamination and an increased risk of disease outbreaks in affected areas. According to data from the European Avalanche Warning Services, there have been 146 fatalities recorded in France from 2018 to the present, with an additional 635 fatalities across Europe [3]. This highlights the critical need for feasible, effective, reasonable, and scientifically grounded investigations to address these challenges. Evidence suggests divergent approaches in evaluating not only the reasons behind their occurrences but also their detrimental impacts on widespread considerations. In this context, creating snow avalanche susceptibility maps and understanding the causes of avalanches have become key focuses for researchers to develop effective mitigation strategies.

The physical mechanism of avalanches involves a complex interplay of factors, including snowpack conditions, terrain characteristics, and weather patterns. Typically, avalanches are triggered by the release of a weak layer within the snowpack, which fails under the weight of new snow or external forces such as skiers, snowmobilers, or explosive triggers used in avalanche control measures [4]. They typically bring out devastating forces, such as dragging down trees and boulders and leading to the blown-down of immeasured snow masses to ground level. Therefore, understanding the triggers of snow avalanches is crucial for assessing risks and developing effective safety measures to prevent potential disasters in mountainous areas. Since the presence of many factors impacting/triggering the occurrence of snow avalanches, deterministic approaches may result in hard-to-understand and complex solutions. On the other hand, advanced computational approaches (such as machine learning techniques or data-driven modeling strategies) ensuring the diligent investigation of a diverse range of conditioning factors and their large volume of datasets are regarded as promising alternatives to holistically interpret the sophisticated nature of snow avalanches. Accordingly, several researchers have recently acknowledged the significant contribution of these accelerated techniques in such climatologic events, e.g., prediction of temperature [5], precipitation [6], humidity [7], solar radiation [8], etc. Likewise, regarding the nature-induced disasters that are particularly related to meteorological conditions, [9,10,11] utilized the corresponding techniques for modeling flash floods, pluvial floods, thunder, and avalanches, respectively.

In recent years, various advanced modeling techniques have been proposed in the relevant literature. In this regard, machine learning algorithms have been shown by researchers to be valuable tools for creating models that replicate real-world challenges and mimic natural conditions. For instance, Yariyan et al. [12] aimed to optimize the machine learning algorithms to produce snow avalanche susceptibility maps in the Zarrinehroud and Darvan watersheds in Iran. They combined four machine learning models, namely radial basis function, multilayer perceptron, fuzzy art-map, and self-organizing map (SOM), with three statistical algorithms, i.e., frequency ratio, statistical index, and weights of evidence, alongside the k-means clustering. According to their analysis, the combination of k-means clustering and SOM (k-means&SOM) outperformed its counterparts, offering a viable solution to detect the regions potentially experiencing snow avalanches. The authors also found the most influential criteria as slope, TWI (topographic wetness index), land use, WEI (wind exposition index), and distance from the stream, respectively, whereas the LS (length slope) and VRM (vector ruggedness measure) provided limited contributions to the predictions. Additionally, Bian et al. [13] created an ensemble machine learning model to discover the snow avalanche susceptibility in the central regions of Shaluli Mountain, Sichuan Province, China. They employed four integrated models (i.e., EBF-LR, EBF-MLP, CF-LP, CF-MLP) obtained by crossing two statistical models, such as evidence confidence function (EBF) and certainty coefficient (CF), and two machine learning models, namely logistic regression (LR) and multilayer perceptron (MLP). In conclusion, CF-MLP yielded the most accurate results, followed by the frequency ratio (FR), which was the benchmarking approach. Similar to the findings of Yariyan et al. [12], they found TWI as one of the most determinant factors, while they also overestimated the contributions of elevation, NDVI, and aspect to the overall prediction scheme. From different aspects encompassing divergent machine learning approaches pertaining to different families, Wen et al. [2] compared the performance of support vector classification (SVC), k-nearest neighbors (KNN), classification and regression tree (CART), and multilayer perceptron (MLP) in determining the avalanche-prone areas in the snowy mountains of the Qinghai–Tibet Plateau. The results showed that all the trained models had good prediction capabilities, but the SVM provided the most robust and accurate estimates based on various performance evaluation metrics. Concentrating on the effectiveness of tree-based machine learning algorithms, Iban and Bilgilioglu [14] explored the avalanche susceptibility levels in the Province of Sondrio, Italy. The authors integrated the Shapley Additive exPlanations (SHAP) technique into the developed machine-learning predictive schemes to extract the contribution of each avalanche conditioning factor on the estimations performed. In this way, the explainability of the established black-box framework has been augmented. As a result, among divergent tree-based algorithms, the extreme gradient boosting (XGBoost) method outperformed its counterparts with the highest predictive accuracies, and it was followed by the gradient boosting (GB), light gradient boosting machine (LGBM), natural gradient boosting (NGBoost), random forest (RF), and adaptive boosting (Adaboost). Regarding the importance level of the attributes included, they found that elevation played the most crucial role, while maximum temperature, slope, and wind speed were the 2nd, 3rd, and 4th significant criteria, respectively. On the other hand, the best-performed model underestimated the role of not only TPI and TRI but also plan and profile curvatures, proximity to streams and roads, and solar radiation.

Recent studies and reports from relevant institutions have highlighted that snow avalanches are significant natural hazards in Alpine areas [14,15,16], emphasizing the need for comprehensive prediction frameworks to investigate such incidents. Hence, the present study employed three different machine learning techniques from three predictive families, i.e., stochastic gradient boosting (SGB) from tree-based models, KNN from distance-based approaches, and SVC from the kernel-based algorithms, in conjunction with three meta-heuristic optimization strategies, namely particle swarm algorithm (PSO), gravitational search algorithm (GSA), and cuckoo search algorithm (CS), in order to obtain the most representative snow avalanche susceptibility map in Alpine regions within the French borders. The major contributions of this study indicating where its novelty comes from are as follows:

• Upon examining the literature, this research is the initial attempt regarding the utilization of the SGB technique, and additionally, it is the first attempt in terms of the integration of such a comprehensive hyperparameter tuning strategy in snow avalanche susceptibility mapping.
• Along with implementing several pre-processing steps for spinning up to the predictive analysis, the present study aimed to discover the role of one of the nascent explainable artificial intelligence techniques, namely the SHAP, in post-processing attempts. To the best of the authors’ knowledge, such a well-rounded and systematic approach has not been acknowledged in the pertinent literature with regard to the designation of susceptible regions to snow avalanches.
• Susceptibility mapping with such a holistic model consisting of a total of 17 variables (including topographic, geological, meteorological, and land use factors) has not yet been implemented in snow avalanche analyses. In this sense, hybrid predictive modeling for sensitive ecosystems like the French Alps is important for creating real-time decision support systems by considering geographic and climatic factors in light of current and constantly varying data.

Hence, it is believed that the current research not only contributed to the relevant literature by disclosing the most effective hybrid prediction strategy but also to the practical implications by addressing the attributes primarily contributing to the occurrences of snow avalanches.

2. Research Framework

The overarching objective of the present research is to delineate regions that are susceptible to snow avalanches through advanced soft computational approaches. In this vein, as a case study, the French Alps, which have been subjected to serious avalanche incidents and their consequences, were chosen. Hence, three different algorithms from three distinct machine learning families were employed. The SVC was chosen from the kernel-based approaches, while the SGB and KNN were selected among tree-based and distance-based algorithms, respectively. This research aimed to evaluate the performance of three optimization techniques (PSO, GSA, and CS) combined with three machine learning algorithms due to the complexity of machine learning algorithms and their hyperparameters. In this regard, a total of 9 scenarios are accomplished, allowing us to perform holistic comparisons in identifying the snow avalanche susceptibility in the French Alps. Additionally, several performance evaluation metrics were considered to assess the predictive ability of the models established. Such that the accuracy metric was selected as an objective function during the calibrations, and some others, including precision, recall, and F1-score, were obtained for assessing the overall generalization ability of the classification tasks. Furthermore, not only widely acknowledged ROC plots but also confusion matrices were generated to visually compare the established hybrid predictive frameworks.

This study also examined the role of explainable artificial intelligence techniques in determining the reasons behind the occurrence of avalanche incidents. To accomplish this, a game-theoretical SHAP algorithm was integrated into the best-performing machine learning predictive strategy. Finally, the snow avalanche susceptibility map with regard to the French Alps comprising a total of 14 massifs next to Italy and Switzerland was delineated, and susceptibility levels across the focalized study domain were categorized into five classes (i.e., very low, low, moderate, high, and very high) for providing significant insights to different purposes, including granular risk assessments, enhanced decision-making schemes, risk prioritization, etc., to be conducted in the following attempts. The corresponding steps implemented within the current research are graphically represented in Figure 1.

3. Study Area

This study covers the southern reaches of the French Alps, a significant segment within the broader Alpine region situated in southeastern France, Europe. While the expansive Alpine Region encompasses seven countries—Austria, France, Germany, Italy, Liechtenstein, Slovenia, and Switzerland (A European Union Strategy for the Alpine Region, n.d.)—the adjacent countries to our specific study area are Switzerland and Italy. The study area is located within the North parallels of 43°29′44.39″ N to 46°46′5.63″ N and the East meridians of 4°13′23.41″ E to 8°53′45.26″ E. The focusing study domain covers 21,387 km², with maximum and minimum elevations of 4803 m and 152 m, respectively (Figure 2).

The French Alps are known for their varied topography, shaped over thousands of years by tectonic activity and glacial erosion, featuring high peaks, rugged mountain ranges, and large plateaus [17]. The region also experiences significant climatic differences, with mild summers and cold, icy winters. The environment at high altitudes creates distinct microclimates that impact the distribution of flora and patterns of precipitation. According to the region’s physical characteristics, the climate in the Alpine area varies greatly [16]. Climates at lower elevations are typically drier and warmer than those at higher altitudes in the vicinity. Most precipitation above 1500 m occurs as snow during the winter months. At heights of approximately 1500 m above sea level, snow cover normally lasts from mid-November to the end of April [18]. Average January temperatures in the valley bottoms are between −5 and 4 °C, but higher altitudes—especially at lower elevations—can see temperatures as high as 8 °C [19]. On the other hand, average July temperatures range from 15 to 24 °C. Temperature inversions are responsible for the extended stagnant air and fog that valleys frequently experience, especially in the fall and winter [20].

According to Habersack and Piégay [21], the French Alps serve as a vital that produces a multitude of rivers and streams, such as the Rhône, Isère, and Durance. The hydrological network of the area is influenced by glacial meltwater, which comes from glaciers and permanent snowfields. This water shapes riverine landscapes and supports a variety of aquatic species. Therefore, understanding how avalanches work is crucial for reducing hazards, preserving biodiversity in alpine environments, and protecting cultural heritage sites vulnerable to avalanche risks.

4. Materials

4.1. Avalanche Inventory Mapping

Using machine learning for snow avalanche forecasting requires a thorough understanding of past avalanche events worldwide, making the creation of avalanche inventory maps essential. However, due to the short-winded nature of avalanches, continuous observation proves challenging. The majority of available avalanche records stem from human-reported incidents, serving as fundamental datasets for predicting future events and constructing avalanche susceptibility maps. Recently, many researchers have studied avalanche dynamics in different areas [1,2,14], highlighting its key role in avalanche risk assessment and management.

This study examines the dataset comprising reported avalanche records sourced from the database of the Data Avalanche Association (www.data-avalanche.org/explore (accessed on 1 February 2024)). The Data Avalanche Association is dedicated to cataloging avalanche incidents worldwide. Based on the avalanche event database maintained by the Data Avalanche Association, a total of 3610 point-based snow avalanche occurrences were documented within the French Alpine Region spanning the years 2000 to 2023.

The Data Avalanche Association website utilizes the JSON data format, facilitating the retrieval of avalanche datasets with a well-designed Python script. This dataset encompasses comprehensive details of recorded avalanches, encompassing coordinates, departure and arrival altitudes, incident timestamps, avalanche movement orientation, and additional pertinent information for each individual avalanche incident. The historical dataset comprises point-based avalanche records diligently sourced from the Data Avalanche database, aided by satellite imagery from Google Earth. Utilizing the widely recognized geographical information system ArcGIS version 10.3 [22], a snow avalanche inventory map is generated from the retrieved dataset. Consequently, an equivalent number of non-avalanche points are randomly distributed across the study area using the Create Random Points function within the ArcGIS 10.3 software as a rule of thumb in susceptibility assessments of different disastrous phenomena [23,24].

4.2. Avalanche Triggering Factors

The complex and unpredictable nature of snow avalanches makes them difficult to study using deterministic approaches. Instead, forecasting and identifying key factors, followed by training machine learning algorithms to evaluate these factors, offers a more effective way to explore avalanches. Several criteria influence the mechanisms affecting avalanche formation, e.g., rainfall, snowfall, elevation, temperature, wind speed, solar radiation, snow depth, snow density, slope, aspect, land use and cover, solar radiation, proximity to faults, and various topographic factors such as the topographic position index (TPI), topographic wetness index (TWI), and topographic roughness index (TRI). Although some factors are beyond the scope of this study due to limited access to relevant datasets, this research focuses on incorporating 17 key avalanche conditioning factors into explainable machine learning models. Thus,

Table 1. List of the snow avalanche triggering factors.

Criterion	Definition	Direct/Indirect Impacts on Snow Avalanches	Data Source
Elevation	The height above sea level of a certain location.	Elevation influences snowpack characteristics, avalanche initiation, and runout.	United States Geological Survey (USGS) [26]
Slope	The steepness of the terrain is typically expressed as an angle or percentage.	Slope impacts the initiation of the snow avalanches.	Retrieved from the Digital Elevation Model.
Aspect	The compass direction that a slope face.	Aspect affects snow accumulation and melting rates, influencing snowpack stability and avalanche occurrence.	Retrieved from the Digital Elevation Model.
Profile Curvature	The curvature of the terrain profile along a slope.	Terrain profile curvature affects snow deposition, wind redistribution, and snowpack stability.	Retrieved from the Digital Elevation Model.
Plan Curvature	The curvature of the terrain is perpendicular to the slope direction.	Plan curvature influences snow distribution and wind loading patterns.	Retrieved from the Digital Elevation Model.
Land use/Land cover (LULC)	The classification and mapping of surface cover types in a geographic area.	LULC affects snow accumulation, stability, and avalanche occurrence.	Coordination of Information on the Environment (CORINE) [27]
Topographic position index (TPI)	A measure of a location’s relative position within a landscape.	TPI influences snow distribution and avalanche behavior as it reflects terrain morphology.	Retrieved from the Digital Elevation Model.
Topographic wetness index (TWI)	A measure of topographic moisture conditions, which is calculated from the ratio of upslope contributing area to the tangential slope.	TWI influences soil moisture, vegetation distribution, and snowmelt rates, which in turn affect snow stability and avalanche potential.	Retrieved from the Digital Elevation Model.
Topographic ruggedness index (TRI)	A measure of terrain roughness or variability.	TRI influences snow distribution, wind transport, and avalanche behavior as it represents the complexity and variability of terrain morphology.	Retrieved from the Digital Elevation Model.
Lithology	The study of the physical and chemical properties of rocks and soil.	Lithology affects snowpack stability and avalanche release through its influence on terrain roughness, slope stability, and snowpack composition.	Hengl [28]
Rainfall	Precipitation is in the form of liquid water falling from the atmosphere.	Rainfall can destabilize the snowpack by increasing water infiltration and percolation, weakening snow layers, and promoting avalanche release.	Thiemig et al. [25]
Wind Speed	The velocity of air movement is typically measured at a certain height above the ground.	Wind speed influences snow transport patterns, deposition, and loading, affecting avalanche release and propagation.	Thiemig et al. [25]
Minimum Temperature	The lowest temperature is recorded within a specific period.	Minimum temperatures affect snow metamorphism, stability, and avalanche conditions.	Thiemig et al. [25]
Maximum Temperature	The highest temperature is recorded within a specific period.	Maximum temperatures influence snowmelt rates, snowpack settlement, and avalanche conditions.	Thiemig et al. [25]
Solar Radiation	The energy received from the sun is typically measured as solar irradiance or insolation.	Solar radiation drives snowmelt, consequently influencing snowpack stability and avalanche conditions.	Thiemig et al. [25]
Snow Depth	The vertical thickness of the snowpack was measured from the ground surface.	Snow Depth is a fundamental indicator of an avalanche potential.	Tao et al. [29]
Distance to Faults	The proximity of a location to geological faults, fractures, or seismic zones.	Faults and geological structures influence terrain stability and snowpack characteristics, potentially serving as release zones for avalanches.	Basili et al. [30]

is presented not only to explain the distinct characteristics of the influential factors incorporated into but also their definitions and sources. Additionally, existing research on snow avalanche susceptibility mapping shows that datasets related to hydrological factors typically focus on the winter months, particularly January. Hence, this study encompasses the hydro-meteorological data (i.e., rainfall, wind speed, minimum and maximum temperatures, solar radiation, and the snow depth) from the three winter months (namely December, January, and February) as the heavy snowfall or rainfall patterns may span over the corresponding temporal periods. The resolution of the remote sensing databases is of critical importance in the accuracy of machine learning models. This study utilized a high-resolution gridded meteorological dataset provided by Thiemig et al. [25]. This dataset, namely European Meteorological Observations with a spatial resolution of 5 km (EMO-5), is based on historical and real-time observations obtained by integrating data from 18,964 ground weather stations across Europe. Additionally,

Table 2. Descriptive statistics of the conditioning factors.

Attribute	Unit	Mode	Min	Mean	Max	Std. Dev.	Skewness	Kurtosis
Elevation	m	1880	174	1945.05	4389	690.10	−0.26	−0.36
Slope	-	16.98	0.00	27.17	72.52	11.52	−0.10	−0.28
Aspect	-	45.00	−1.00	184.02	359.82	105.94	−0.06	−1.28
Profile Curvature	-	−0.07	−5.06	0.00	4.35	0.57	−0.01	6.25
Plan Curvature	-	0.00	−3.52	0.01	4.12	0.45	0.30	5.17
LULC *	-	-	1	-	41	-	-	-
TPI	-	−10.12	−235.26	15.26	319.79	78.61	0.37	0.34
TWI	-	6.44	2.85	6.08	15.25	1.60	1.52	3.37
TRI	-	0.40	0.00	0.47	0.81	0.10	−0.48	0.71
Lithology *	-	-	2	-	96	-	-	-
Rainfall	mm	0	0	87.37	161.42	27.75	0.29	−0.05
Wind Speed	m/s	0	0	2.20	4.60	0.46	−0.17	3.26
Minimum Temperature	°C	0	−20.87	−8.06	4.37	3.86	0.31	−0.47
Maximum Temperature	°C	0	−14.62	−0.28	10.84	3.59	0.11	−0.18
Solar Radiation	Joule/m2	0	0	6,102,559.36	7,438,510	817,567.33	−2.67	19.41
Snow Depth	cm	43.23	2.68	38.05	80.03	22.17	0.27	−0.75
Distance to Faults	km × 10−3	0	0	0.23	0.94	0.15	0.86	1.05
Target	-	-	0.00	-	1.00	-	-	-

Note: * refers to the categorical attributes.

displays the descriptive statistics of the employed snow avalanche conditioning factors in a way that provides insights regarding the nature of the 17 independent variables. This study further ensured their respective maps in Figure 3.

5. Methods

5.1. Pre-Processing

5.1.1. Data Encoding

Although machine learning applications typically rely on numerical data, real-world datasets often contain categorical variables. In such cases, data pre-processing becomes necessary before applying predictive algorithms. However, these pre-processing steps should enhance the meaningfulness of the overall process. Various approaches are used to convert categorical data into numerical form, with the method chosen depending on which best preserves the original characteristics of the data. In pertinent literature, two of the most commonly used techniques are label encoding and one-hot encoding. Label encoding directly converts categorical variables into numerical values, where each category is represented by a unique number. However, in cases where the hierarchical relationship does not exist, other techniques can yield superior performance compared to label encoding. In such scenarios, one-hot encoding is often more effective, as it better reflects the structure of the original data. One-hot encoding works by transforming categorical data into binary vectors, assigning a separate binary column to each category. While this increases the dataset’s dimensionality by adding new columns, it is particularly advantageous when dealing with categorical variables that have nominal relationships. One-hot encoding, therefore, offers an efficient solution for machine learning algorithms. In this study, the one-hot encoding technique was applied to transform the categorical variables of LULC and lithology, which were among the input features.

5.1.2. Data Scaling

Dealing with variables varying within significantly different ranges and addressing the challenges due to the presence of extreme values in datasets play an important role in enhancing prediction accuracy. To overcome the corresponding challenge, the relevant literature recognizes the utilization of scaling techniques, e.g., min–max scaling, standardization, robust scaling, and logarithmic transformations. Min–max scaling is performed by considering the minimum and maximum values in attributes interested, essentially ensuring that all data take values between 0 and 1. On the other hand, in standardization, the z-score is calculated for the respective attribute. Robust scaling, unlike other techniques, takes the median value and interquartile range of the dataset into account to realize the scaling mechanism. This provides significant advantages, especially when dealing with outlier instances. Additionally, the relationship between data points is maintained in this regard, and the underlying structure of the data remains intact.

5.1.3. Data Splitting Strategy

The selection of the train/test splitting rationale is of significance in machine learning implementations. In this sense, the pertinent literature contains diverse mechanisms, including 50/50, 60/40, 70/30, 80/20, and 90/10 for training/testing sets. However, dividing the entire test into two different groups, 70% and 30% for training and testing, respectively, can be regarded as the most widely endorsed strategy. Hence, this study allocated 70% of the training set, while the remaining 30% is ensured as the blinded testing set. Additionally, within the present research, a rigorous 5-fold cross-validation is employed to delineate the optimum hyperparameter sets during the optimization attempts. It is worth noting that the corresponding technique not only facilitates robust predictive efforts but also enhances the generalization capacity of the classification analysis conducted, mitigating the risks associated with overfitting and underfitting phenomena.

5.2. Processing

5.2.1. Optimization Algorithms

Optimization algorithms are computational methods designed to find the best possible solution to a problem from a set of solutions. These algorithms iteratively explore the solutions to identify the optimal solution, which generally maximizes or minimizes an objective function while satisfying a predefined set of search spaces. Given their use in various fields for efficient decision-making and problem-solving, three meta-heuristic optimization techniques are employed in this study.

Particle Swarm Optimization (PSO) is an artificial intelligence technique that ensures the extraction of approximate solutions that pose significant challenges to solving minimization and maximization problems. It was formulated by Kennedy and Eberhart [31] with the inspiration and observation of a food-seeking flock of birds. PSO mainly focuses on the position and velocity information of randomly distributed individuals in the swarm. Individuals in the swarm are commonly called particles; meanwhile, the swarm itself is called the population [32]. Through each iteration, the position and velocity of each particle change with respect to the modifications made in the input parameters. Afterward, the pre-update and post-update position and velocity of a particle are compared to find the best-performing particle in the iterations performed. The solution lies between the overall information of the swarm and the particles themselves. Thus, the optimization process terminates once the position of the particles does not alter anymore or the computational limits are met [33,34].

This research also used the Gravitational Search Algorithm (GSA) in determining the hyperparameters of the employed machine learning algorithms. The GSA was first introduced by Rashedi et al. [35] as a stable metaheuristic optimization algorithm. Mechanism underlying the GSA based on the law of gravity and interactions between masses. In this algorithm, searcher agents are a set of masses that interact with each other based on Newtonian gravity and the laws of motion [36]. According to its working mechanism, agents having a heavier mass attract more particles than less heavier ones. Subsequently, the position of the agent with the highest number of particles leads to the optimal solution in the search space [37].

Lastly, the Cuckoo’s Search (CS) algorithm, proposed by Yang and Deb [38], was used in conjunction with the predictive models. Although the GSA brought in inspiration coming from the breeding habits of cuckoo birds, Rajabioun [39] formulated the CS algorithm as an optimization technique. In nature, cuckoo birds leave their eggs in another bird species’ nests, and the host bird accepts the cuckoo egg. When the egg matures, it immigrates to another habitat for breeding. The newly mature egg follows the same steps that its parents have already taken. In the optimization context, a set of cuckoo birds is defined with a certain number of eggs according to the maximum and minimum limits defined in the algorithm [40]. Additionally, the maximum distance a cuckoo bird can take is defined in the algorithm as well. Similar to the breeding nature of cuckoo birds, the algorithm iterates each generation of cuckoo bird groups and tracks breeding locations. By determining the distance between the center point and cuckoo bird locations and updating the center point, the optimal solution is yielded as the position of the center point [41]. These iterations go on until the location of the center point stabilizes or computational limitations exceed [42]. The hyperparameters of the utilized optimization algorithms are provided in

Table 3. The utilized parameters of optimization algorithms.

Algorithm	Parameter	Abbreviation	Value	Reference
PSO	Population Size *	$P$	500	-
	Number of Population *	$NP$	250	-
	Cognitive component	$C_1$	2.8	Karaguzel et al. [43]
	Social component	$C_2$	1.45	Ibrahim et al. [44]
	Inertial weight	$w$	0.3	Singh et al. [45]
	Minimal velocity	$v_{Min}$	0.1	Anter and Hassenian [46]
	Maximal velocity	$v_{Max}$	0.9	Anter and Hassenian [46]
GSA	Population Size *	$P$	500	-
	Number of Population *	$NP$	250	-
	Gravitational Constant	$G_0$	50	Amin [47]
	Number of Masses	$D$	20	Koc et al. [48]
CS	Population Size *	$P$	500	-
	Number of Population *	$NP$	250	-
	Fraction	$P_a$	0.25	Shehab et al. [49]
	Step size	$a$	1	Zhang et al. [50]

Note: * refers to the user-defined parameters.

.

5.2.2. Machine Learning Algorithms

Machine learning algorithms are computational techniques used to enable computers to teach themselves from a given dataset and improve their performance on a given task over time. With the help of these algorithms, computers can identify patterns, make predictions, or make decisions without being explicitly programmed for each specific task. In this study, three different machine learning algorithms (i.e., SVC, SGB, and the KNN) are employed, and their brief descriptions are provided in the following paragraphs.

Support vector classification (SVC), proposed by Vapnik [51], is a universal structural learning process based on statistical learning theory. Its reliance on structural risk minimization distinguishes it from many algorithms (such as neural networks). SVC applies the principle of structural risk minimization, which has been shown to achieve better performance in terms of experimental risk, to predict a single and optimal separating plane in the hidden feature space using quadratic programming [52]. By mapping input vectors into multidimensional feature space, SVC transforms many problems with complex structures into simple forms, such as the application of linear discriminant functions. The SVC accurately classifies data by assigning them to their respective zones. This classification process relies on the use of a straight line known as a hyperplane. The optimal hyperplane, referred to as the margin, is strategically positioned to maximize the distance between two distinct classes [53]. When determining the hyperplane that linearly separates the instance space, only the boundary values significantly influence its placement. Changes in the remaining data points do not affect the positioning of the hyperplane. Even with very small datasets, SVM tends to face fewer challenges with regard to the overfitting phenomena. In the SVC application, kernel functions play a crucial role in transforming input data into a higher-dimensional space, enabling the model to find a hyperplane that effectively separates data points of different classes. Kernel functions achieve this by implicitly mapping the original feature space into a higher-dimensional space where the data points become linearly separable [54]. The most commonly used Kernel functions are Linear Kernel, Polynomial Kernel, Sigmoid Kernel, Radial Based Kernel (RBF), and Pearson VII Kernel (PUK) Function. Choosing the most appropriate function is highly effective in the success of classification. The selection of the kernel function depends on the data characteristic, and it can also be selected separately by trial and error as well as using different optimization algorithms. This study incorporated meta-heuristics into the machine learning algorithms to tune these hyperparameters. Studies show that SVC can handle various kernel functions, works well with small datasets [55], and finds the best hyperplane to separate classes, improving generalization on unseen data [56]. However, for large datasets, the algorithm can be computationally intensive, requiring significant memory and storage resources.

Friedman [57] introduced stochastic gradient boosting (SGB) by adding bootstrap aggregation to gradient boosting [58], making it more resistant to overfitting and improving accuracy and robustness [59,60]. In a typical SGB implementation, weak learners (or weak predictors) performing slightly better than the random chance are used to generate the ensemble models by embracing a step-wise approach [61]. In other words, simple trees are sequentially placed, utilizing the gradient of the loss function from the preceding tree to accentuate focus on poorly modeled observations. At each iteration, a random subsample of the training dataset (without replacement) is utilized as input, in which this randomness incorporated into the training process enhances model robustness and reduces overfitting. The SGB is considered more computationally efficient than traditional gradient boosting, especially with large datasets, as it processes only a subset of the data at each iteration [62,63]. Rather than developing individual complex trees, relatively small trees are amalgamated by averaging their weighted predictions, promoting both model interpretability and performance. Gradient boosting stands out for its ability to improve generalizability by optimizing various differentiable loss functions. However, introducing randomness can increase variance in predictions, affecting the stability of the analysis [64]. The SGB algorithm also requires significant attention in model construction, particularly when taking its sensitivity to its hyperparameters into account.

The k-nearest neighbor, proposed by Fix and Hodges Jr. [65], is one of the commonly applied classification techniques utilized in supervised data mining algorithms. The class assignments for observations are identified by assessing the proximity of their nearest neighbors based on a predetermined/specified k-value. In the KNN, the classification decision for a certain data point is influenced by the classes of its k nearest neighbors in the attribute space. Being a non-parametric method, it does not assume any underlying probability distributions for the data in order to perform the predictions. Hence, it only relies on the distances between data points to conduct the estimations. The algorithm aims to accurately assign individuals or objects to predefined classes or groups by leveraging their inherent properties [66]. Additionally, it offers classification for new observations. When classifying a new observation, the algorithm identifies its closest k neighbors—those with the highest similarity—within the dataset used for learning. This dataset, which forms the basis for model creation, is referred to as the training part. The k-nearest neighbor method offers numerous advantages, including its ability to provide interpretable and impactful results, its capability to handle missing observations in continuous variables, its flexibility in evaluating missing observations in categorical variables, and its competitiveness in terms of its computational efficiency [67]. This method accommodates response variables that can be categorical, continuous, or a blend of both while requiring minimal assumptions due to its non-parametric nature. However, the algorithm also presents certain drawbacks, such as necessitating determining the number of nearest neighbors to consider and being sensitive to the choice of distance measure. Additionally, there is uncertainty regarding the selection of the appropriate distance measure, which can influence the algorithm’s performance. It is also argued that the algorithm needs qualified datasets in order to yield satisfactory outcomes, especially where the input space is extended and the number of instances is increased [68]. The search space of the utilized machine learning algorithms’ hyperparameters and the total combinations scanned are illustrated in

Table 4. Parameter ranges in ML methods.

ML Method	Parameters	Parameter Ranges	Step	Count	Total Combination
SVC	Kernel Function	Polynomial, Radial Basis, Sigmoid	-	3	600
	Gamma	2−15–23	-	10
	C	0–200	10	20
SGB	Number of trees	0–500	50	10	500
	Learning rate	0.0025–0.015	0.0025	5
	Maximum Depth	0–10	1	10
KNN	Metric	Euclidean, Manhattan, Chebyshev, Minkowski	-	4	200
	Number of neighbors	0–50	1	50

.

5.2.3. Performance Evaluation

In classification analysis with machine learning models, there are several performance indicators that allow us to assess the predictive success of the utilized algorithms. It is especially worth mentioning that the vast majority of the corresponding metrics are grounded on four statistical measures: (i) true positive (TP) depicting the number of correctly classified instances in the positive label, (ii) true negative (TN) indicating the number of correctly classified instances in the negative label, (iii) false positive (FP) reflecting the number of incorrectly classified instances in the negative label, and (iv) false negative (FN) representing the number of incorrectly classified instances in the positive label (

Table 5. Possible cases in binary classification.

		Predicted Cases
		Yes	No
Observed Cases	Yes	True Positive (TP)	False Negative (FN)
	No	False Positive (FP)	True Negative (TN)

). Based on the aforementioned four cases, this research computed four different performance metrics, i.e., accuracy, precision, recall, and F1-score.

Accuracy, which is the widely accepted indicator, can be calculated through the total number of correctly classified items divided by the total number of items. Hence, Equation (1) is used for the corresponding calculation:

$$\mathrm{Accuracy}={\displaystyle \frac{\mathrm{TP}+\mathrm{TN}}{\mathrm{TP}+\mathrm{FN}+\mathrm{FP}+\mathrm{TN}}}$$

(1)

On the one hand, precision is obtained through the number of correctly classified positive items divided by the total number of positively classified items, as presented in Equation (2). On the other hand, recall is the ratio of correctly classified positive items to the total number of observed positive items, as indicated in Equation (3). It can, therefore, be inferred that the first focalizes on the false positive instances, whereas the latter concentrated on the false negative instances.

$$\mathrm{Precision}={\displaystyle \frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FP}}}$$

(2)

$$\mathrm{Recall}={\displaystyle \frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}}$$

(3)

In addition, as being one of the most used performance assessment criteria in the pertinent literature, the F1-score combines the information received from both precision and recall by computing the harmonic mean of them.

$$\mathrm{F}1-\mathrm{score}={\displaystyle \frac{2\times \mathrm{Precision}\times \mathrm{Recall}}{\mathrm{Precision}+\mathrm{Recall}}}$$

(4)

The present research further incorporated the area under receiver operating curve (AUROC) into the model evaluations in order to judge the predictive accuracy holistically. It is based on the calculation of the area under the ROC curve encompassing false positive rate (i.e., 1-specificity) and true positive rate (i.e., sensitivity) on the x- and y-axis, respectively. The lower and upper boundaries of the AUROC are 0.5 and 1, respectively, in which approaching 0.5 shows low predictive success, while 1 highlights the perfect match. Equation (4) expresses how AUROC is computed:

$$AUC=1-{\displaystyle \frac{1}{m^{+}{m}^{-}}}{\displaystyle \sum }_{x^{+}ϵ M^{+}} {\displaystyle \sum }_{x^{-}ϵ M^{-}}\left(\left(f\left(x^{+}\right)<f\left(x^{-}\right)\right)+{\displaystyle \frac{1}{2}}\left(f\left(x^{+}\right)=f\left(x^{-}\right)\right)\right)$$

(5)

where $m^{+}$ and $m^{-}$ are positive and negative instances, $M^{+}$ and $M^{-}$ are set of all positive and negative instances, and $x^{+}$ and $x^{-}$ are positive and negative classes, respectively. In the formula, $f\left(x\right)$ denotes the outcome of prediction model sample $x$.

Along with the well-known accuracy, precision, recall, F1-score, and AUC, this research further extended the models’ evaluation using two additional performance measures, namely the Matthews Correlation Coefficient (MCC) and Cohen’s Kappa index. The MCC is generally employed in ML implementations to evaluate the performance of binary classification attempts. It takes the true and false positives and true and false negatives into account, making it a balanced metric that remains reliable and comprehensive. Hence, the MCC is formulated as follows:

$$\mathrm{MCC}={\displaystyle \frac{\left(\mathrm{TP}\times \mathrm{TN}\right)-\left(\mathrm{FP}\times \mathrm{FN}\right)}{\sqrt{\left(\mathrm{TP}+\mathrm{FP}\right)\left(\mathrm{TP}+\mathrm{FN}\right)\left(\mathrm{TN}+\mathrm{FP}\right)\left(\mathrm{TN}+\mathrm{FN}\right)}}}$$

(6)

Representing a correlation coefficient between −1 and +1, despite arbitrary limits, 0–0.19 is regarded as very weak, 0.2–0.39 as weak, 0.40–0.59 as moderate, 0.6–0.79, as strong and 0.8–1 as very strong correlation as a general assessment. Additionally, Cohen’s Kappa index is regarded as a statistical index used to assess the degree of agreement between two categorical datasets. As a robust index providing inter-classifier reliability, it takes agreement occurring by random chance into account, making it a useful tool for assessing whether the predictions of the utilized model are better than a random guess. Equation (7) illustrates the simple computation of Cohen’s Kappa index.

$$\mathsf{\kappa }={\displaystyle \frac{{\mathrm{p}}_0-{\mathrm{p}}_{\mathrm{e}}}{1-{\mathrm{p}}_{\mathrm{e}}}}$$

(7)

where ${\mathrm{P}}_0$ depicts the observed proportion of agreement between the observed and predicted classes, while ${\mathrm{P}}_{\mathrm{e}}$ is the expected proportion of agreement by chance. In general, Kappa values under 0.2 represent poor agreement, while 0.2–0.39, 0.40–0.59, and 0.6–0.79 denote fair, moderate, and substantial agreement, respectively. Values equal to and greater than 0.8 are the indicators of a strong agreement between two classes.

5.3. Post-Processing

Improving the explainability and interpretability of predictive outcomes has long been a key research focus, particularly due to the black-box nature of many machine learning techniques [69,70]. Several methodologies have been proposed, with SHapley Additive exPlanations (SHAP) becoming one of the most widely adopted due to its numerous advantages. SHAP provides a detailed understanding of feature importance by attributing contributions to each feature for individual predictions. It is model-agnostic, meaning it can be used to interpret any machine learning model regardless of its underlying architecture. Additionally, it handles complex models with high-dimensional feature spaces and non-linear relationships, including interactions between features [71,72,73]. In addition, SHAP ensures the generation of intuitive and visually appealing explanations via several graphical representations (e.g., summary plots, individual instance plots, dependence plots, etc.).

The SHAP, pioneered by Shapley [74], is a model-independent technique that is formulated based on the game theory approach. It assumes that the outcome of the estimation model can be explained by an additive attribution rationale highlighting the linear incorporation of the independent variables into the predictive scheme [75]. Let an original machine learning model $f\left(\mathrm{x}\right)$ comprise input variables $\mathrm{x}=\left({\mathrm{x}}_1,{\mathrm{x}}_2,\dots ,{\mathrm{x}}_i\right)$, in which the number of predictors is denoted as $i$.

$$f\left(\mathrm{x}\right)=g\left({\mathrm{x}}^{\prime }\right)={\mathsf{\varphi }}_0+{\displaystyle \sum }_{i=1}^M{\mathsf{\varphi }}_i{x}_i^{\prime }$$

(8)

In the equation, ${x_i}^{\prime }$ and $g\left({\mathrm{x}}^{\prime }\right)$ refer to the simplified predictor and the explanation model, respectively. Also, $M$ denotes the number of instances considered in predictions and ${\mathsf{\varphi }}_0$ is the constant value when all inputs are missing. It is worth noting that ${\mathsf{\varphi }}_0$ accounts for different values with regard to the predictions made in subsequent steps, where they can either increase or decrease the predicted value of $g()$. According to Lundberg and Lee [76], who proposed the SHAP algorithm, the equation can only provide a single solution with local accuracy , ensuring that the output of the function is the sum of the feature attributions, missingness ensuring that no importance is assigned to missing features, and the consistency that changes a larger impact feature will not decrease the attribution assigned to the corresponding feature are satisfied. Hence, the Shapley values ${\mathsf{\varphi }}_i$ can be computed through:

$${\mathsf{\varphi }}_{\mathrm{i}}\left(\mathrm{f},\mathrm{x}\right)={\displaystyle \sum }_{{\mathrm{z}}^{\prime }\subseteq {\mathrm{x}}^{\prime }}{\displaystyle \frac{\left|{\mathrm{z}}^{\prime }\right|!\left(\mathrm{M}-\left|{\mathrm{z}}^{\prime }\right|-1\right)!}{\mathrm{M}!}}\left[{\mathrm{f}}_{\mathrm{x}}\left({\mathrm{z}}^{\prime }\right)-{\mathrm{f}}_{\mathrm{x}}\left({\mathrm{z}}^{\prime }\setminus \mathrm{i}\right)\right]$$

(9)

in which $\left|{\mathrm{z}}^{\prime }\right|$ depicts the number of non-zero entries in setting ${\mathrm{z}}^{\prime }$ satisfying the abovementioned conditions, where ${\mathrm{z}}^{\prime }\subseteq {\mathrm{x}}^{\prime }$.

As being a model-agnostic approach, Shapley values can be attained by integrating it into different machine learning families, and accordingly, kernel SHAP, deep SHAP, and tree SHAP can be achieved with regard to employed methodologies. The present research adopted the SHAP algorithm concerning the best-performed model (i.e., the SGB, a tree-based machine learning algorithm) as found and demonstrated in the results, and therefore extracted the attributes’ contributions using the tree-SHAP. Using the coalition game theory elucidating the multifaced challenging interactions [77,78], the Shapley values ensuring stable, transparent, and consistent feature rankings were attained. The outcomes were visualized through the SHAP summary plots, where not only the importance levels of the model predictors but also their individual impacts on the outcomes attained with respect to the alterations in actual feature values are provided.

6. Results and Discussion

This section covers the training results of the hyperparameter tuning processes (Section 6.1), validation of the obtained findings with respect to the testing set and the predictive performances (Section 6.2), and the attributes’ contributions to the performed susceptibility analysis using the explainable machine learning (Section 6.3).

6.1. Training and Validation Results

In the present research, a holistic predictive framework was established with the integration of three meta-heuristics (i.e., PSO, GSA, and CS) for conducting hyperparameter optimization processes and three machine learning algorithms (i.e., SVC, SGB, and KNN) to perform the binary classification analysis. To accomplish this, a total of nine different scenarios were employed. The details for the training and validation processes (i.e., 5-fold cross-validation) are provided in

Table 6. Comparison of optimization results.

Scenario	Best Candidate (G; C)	Duration (s)	Mean Training Accuracy	Mean Testing Accuracy	Optimum Hyperparameters
PSO-SVC	G: 6 C:54	1834	0.8416	0.7956	Kernel: RBF Gamma: 0.0078125 C: 150
PSO-SGB	G: 186 C:63	2499	0.8933	0.8161	Number of trees: 250 Learning rate: 0.0125 Maximum Depth: 8
PSO-KNN	G: 4 C:64	603	0.8276	0.7953	Metric: Manhattan Number of neighbors: 11
GSA-SVC	G: 3 C:74	1936	0.8528	0.7982	Kernel: RBF Gamma: 0.003125 C: 10
GSA-SGB	G: 4 C:57	2533	0.8859	0.8149	Number of trees: 250 Learning rate: 0.0075 Maximum Depth: 10
GSA-KNN	G: 5 C:106	763	0.8276	0.7953	Metric: Manhattan Number of neighbors: 11
CS-SVC	G: 11 C:308	8685	0.8528	0.7982	Kernel: RBF Gamma: 0.003125 C: 10
CS-SGB	G: 11 C:325	8561	0.8933	0.8161	Number of trees: 250 Learning rate: 0.0125 Maximum Depth: 8
CS-KNN	G: 4 C:88	1281	0.8276	0.7953	Metric: Manhattan Number of neighbors: 11

Notes: Calculations were performed using Google Colaboratory (Browser-based Python code execution platform). G: Generation, C: Candidate.

. This table contains the best generation and candidate composition for each hybrid model, the duration required for training and validating each model, and mean training accuracy and standard deviation for training based on the adopted scenarios. In addition, optimum hyperparameters are also provided in the table for all predictive frameworks.

Table 6 is quite revealing that the best training performance was achieved by both PSO-SGB and CS-SGB models with the same accuracy of 0.8933. According to the table, these models were followed by the GSA-SGB. It shows that the SGB model outperformed its counterparts, i.e., SVC and KNN, in all meta-heuristic integrations. Concerning the comparison of the SVC and KNN techniques, the SVC yielded superior performance regardless of the meta-heuristics that they were incorporated into. Such that the accuracies were found as 0.8416 and 0.8276 for PSO-SVC and PSO-KNN, respectively, while the GSA-SVC and GSA-KNN produced accuracies of 0.8528 and 0.8276, respectively. Likewise, in conjunction with the CS algorithm, the predictions made with the SVC resulted in an accuracy of 0.8528, which is greater than the one obtained through the KNN (0.8276) approach. Also, as for the KNN technique, the three meta-heuristics ensured the same hyperparameter configurations by taking the distance metric as “Manhattan” and the number of neighbors as “11” into account.

In machine learning, accuracy is not the only focus; the practical applicability of the models is also important. This depends on the computational efficiency of the predictive frameworks. Therefore, this research evaluated the feasibility of the adopted hybrid algorithms based on their computational efficiency. In this sense, the duration of each analysis is recorded, and the evaluations regarding the computational expense are totally fair as all the analyses were conducted via Google Colaboratory, which is a browser-based Python code execution platform. From this facet, one can conclude that the KNN is the most computationally efficient technique with lesser time requirements, whereas the accuracies attained through the KNN were not as satisfactory as the predictive performances achieved through the other two ML models. What also stands out in Table 6 is the similar concern between the most determinant two models, i.e., PSO-SGB and CS-SGB. Here, two meta-heuristics assigned the optimal model configuration for the SGB algorithm with exactly the same hyperparameter values. Thus, the number of threes is identified as 250, while the learning rate and maximum depth were found as 0.0125 and 8, respectively. Overall, by also concerning the trade-offs between computational efficiency and predictive accuracy, this research pinpoints the superiority of the hybrid PSO-SGB model (Table 6).

In addition, this research presents the convergence graphs of the employed hybrid models for the validation set in Figure 4, Figure 5 and Figure 6 with regard to the PSO, GSA, and CS, respectively. The y-axis in these figures shows the change in the objective function (i.e., accuracy as one of the most widely used performance indicators in binary classification). Furthermore, in the figures positioned in the left column, the x-axis displays the number of iterations, while those in the right column represent the candidate solutions. Similar to the training results, the PSO-SGB is the best-performed model (Figure 4), as it produced the highest accuracy (0.816), followed by the PSO-SVC (0.796) and PSO-KNN (0.795). The figure further demonstrates that, although the PSO-KNN model converged at the 4th generation and 64th candidate, the PSO-SGB model yielded the most accurate results by converging the optimal solution at the 186th generation and 63rd candidate. What is striking about Figure 5 is that the SGB model tuned by the GSA (i.e., GSA-SGB) not only converged faster with 3rd generation at 57th candidate but also produced the most accurate outcomes with 0.815 accuracy compared to its counterparts, i.e., GSA-SVC (0.798) and GS-KNN (0.795). Likewise, from Figure 6, one can conclude that the SGB model hybridized with the CS demonstrates the highest predictive power within the validation set with an accuracy of 0.816, followed by the CS-SVC (0.798) and CS-KNN (0.795). However, it could also be noted that the SGB model optimized through the CS algorithm converged to the optimum solution at the highest number of generation and candidate solutions at 186th and 63rd, respectively.

6.2. Avalanche Susceptibility Mapping with Respect to the Testing Results

Along with the training and validation results with respect to accuracy scores, the present research further evaluated the performance of the established models using different performance indicators, such as precision, recall, F1-score, and AUROC. Hence, the outputs of this holistic evaluation scheme are summarized in

Table 7. Summary of the model performances.

	Performance Measures
	Training					Testing
Scenario	Precision	Recall	F1-Score	MCC	Kappa	Precision	Recall	F1-Score	MCC	Kappa
PSO-SVC	0.8560	0.8508	0.8503	0.6803	0.6753	0.8115	0.8089	0.8085	0.6126	0.6089
PSO-SGB	0.8908	0.8877	0.8875	0.7906	0.7876	0.8237	0.8214	0.8211	0.6541	0.6805
PSO-KNN	0.8372	0.8308	0.8299	0.6725	0.6653	0.8026	0.7936	0.7922	0.5956	0.5896
GSA-SVC	0.8560	0.8508	0.8503	0.6803	0.6753	0.8115	0.8089	0.8085	0.6126	0.6089
GSA-SGB	0.8882	0.8851	0.8849	0.7852	0.7557	0.8228	0.8205	0.8202	0.6473	0.6391
GSA-KNN	0.8372	0.8308	0.8299	0.6725	0.6653	0.8026	0.7936	0.7922	0.5956	0.5896
CS-SVC	0.8560	0.8508	0.8503	0.6803	0.6753	0.8115	0.8089	0.8085	0.6126	0.6089
CS-SGB	0.8908	0.8877	0.8875	0.7906	0.7876	0.8237	0.8214	0.8211	0.6541	0.6805
CS-KNN	0.8372	0.8308	0.8299	0.6725	0.6653	0.8026	0.7936	0.7922	0.5956	0.5896

Note: Avalanche Support = 5041.

. These results are acquired concerning the best candidate solutions. The table shows that the SGB integrated with the PSO was the most powerful method in detecting avalanches in the French Alpes. Such that this hybrid model produced the highest precision (0.8908 for training and 0.8237 for testing), recall (0.8877 for training and 0.8214 for testing), and F1-score (0.8875 for training and 0.8211 for testing). Likewise, the highest MCC and Cohen’s Kappa values were observed using the PSO-SGB model with 0.7906 and 0.7876 for training, respectively. For the testing set, this model yielded the corresponding measures as 0.6541 and 0.6805, respectively. These results indicated that the best-performed model demonstrated a strong correlation and substantial agreement between the observed and predicted instances. Similar outcomes can also be obtained from the confusion matrices illustrated in Figure 7. For instance, the PSO-SGB model (Figure 7b) predicted 86.23% of the regions that experienced avalanche incidents correctly, while this hybrid framework also captured 78.08% of the non-avalanche points. This model was followed by the GSA-SGB, leading to acquiring satisfactory prediction performance with a precision of 0.8115, as well as 0.8089 and 0.8085 in terms of recall and F1-score, respectively (Table 7). The selection of the kernel is of critical importance in configuring the SVC models. Therefore, among the candidate solutions, accuracy values with regard to the other two kernel functions are derived for the testing set. Such that for the SVC models configured using the PSO, we concluded that the highest testing set accuracy values for the polynomial and sigmoid functions were attained as 0.781 and 0.742, respectively. For the models configured via the GSA, they computed as 0.780 and 0.742, respectively, and for those optimized through the CS, the accuracy values pertaining to the corresponding kernels were found as 0.780 and 0.742 as well.

On the other hand, the predictions made by using the KNN technique produced the least accurate results in detecting the avalanche and non-avalanche points, as can be seen from Figure 7c,f,i. Although these integrated algorithms predicted the avalanche points slightly better than the other hybrid models (87.81%), they underestimated the non-avalanche regions with significant differences. Such that they were only capable of detecting 70.99% of the regions not susceptible to avalanche incidents. The corresponding performance made the KNN-based models least effective according to the overall estimations. Despite exactly the same results achieved via the PSO-SGB and CS-SGB, the computational efficiency of the PSO-SGB with lesser processing time (Table 6) showed its effectiveness among all the established hybrid frameworks.

As indicated in Figure 8, the outcomes were also assessed using the AUROC metric, which is one of the most reliable indicators in binary classification analysis. Associatively, among the PSO-based models (Figure 8a), the PSO-SGB outperformed its counterparts with an AUROC of 0.900, while the 2nd and 3rd models were the PSO-SVC and PSO-KNN with 0.876 and 0.871 AUROC values, respectively. In addition, the GSA-based analysis provided similar outputs, such that the predictions made by means of the SGB technique yielded the highest accuracy with an AUROC of 0.899 (Figure 8b). As also mentioned earlier, CS-SGB ensured the highest predictive success among not only CS-based models but also all hybrid approaches. Upon examining the prediction results based on the testing set, CS-SGB was followed by the CS-KNN and CS-SVC models offering 0.871 and 0.849, respectively.

The respective results significantly corroborate with the existing body of knowledge. Such that as a result of implementing different methodologies, Bian et al. [13] concluded that the machine learning models ensure high predictive capacity compared to other techniques, such as statistical techniques like frequency ratio. Likewise, Wen et al. [2] showed the extensive capability of machine learning models to identify snow avalanche-prone areas, highlighting the superiority of the SVC over the KNN algorithm. On the other hand, the comparison of ensemble machine learning models has long received significant attention from divergent disciplines. In specific to spatial modeling of snow avalanches, Akay [79] found the AdaBoost M1 ensemble of random trees, which is a tree-based hybrid machine learning algorithm, as the most powerful predictive model in comparison to several machine learning techniques (such as J48, random forest, bagging, random subspace, etc.).

Furthermore, the major objective of the current research is to generate an avalanche susceptibility map for the French Alpes. To do this, the best-performing model was utilized, and the susceptibility levels were determined in this regard. Hence, the final map was generated through the PSO-SGB hybrid framework, as demonstrated in Figure 9. The figure categorizes the susceptibility levels of the region into five: very low, low, moderate, high, and very high. Based on this figure, one can conclude that slightly more than a quarter of the region has very low susceptibility to avalanche phenomena, and these regions are mainly distributed along the west and south parts of the focalized region. Also, nearly 20% of the entire region is under low avalanche susceptibility, in accordance with the results obtained through the predictive analysis. The figure clearly demonstrates that the susceptibility of the region to avalanche events tends to increase from west to east parts of the French Alpes. Such that 16% of the entire region concentrated on the mid-part posed moderate susceptibility to avalanches. When it comes to the mid-east and east parts of the region, susceptibility levels can mostly be identified as high and very high. For instance, nearly 15% of the region is found to be under a high susceptibility level, whereas regions covering more than 20% of the study area are characterized as having very high susceptibility to avalanche events. In addition, the east part of the French Alpes is clearly found to have very high avalanche susceptibility, which is further distributed to some northern and eastern parts.

6.3. Model Interpretability

The present study integrated the SHAP algorithm into the SGB optimized with the PSO model, which outperformed its counterparts, considering its predictive ability and computational efficiency. As a result, the interpretability of the corresponding hybrid approach is carried out via the SHAP summary plot, as depicted in Figure 10. It is important to mention that the x-axis depicts the Shapley values, ensuring the assessment of the importance levels of the attributes given in the y-axis. The predictors are sorted top to bottom based on their mean absolute Shapley values on the left-hand side of the y-axis, and the right-hand side demonstrates that the red dots represent higher actual values of these predictors and towards the blue color; their values tend to decrease. In this regard, the slope is the most determinant factor in identifying the regions that are susceptible to snow avalanches, while it was followed by elevation, wind speed, maximum and minimum temperatures, and solar radiation, respectively. According to the figure, the higher the slope, the higher the snow avalanche susceptibility of the region. This outcome supports the relevant findings, as most of the literature highlights the significance of slope incline and slope stability in determining snow avalanches [80,81], and the studies showed that snow avalanche incidents occur on slopes between 30° and 45° [82,83]. For instance, Yariyan et al. [84] found that the slope is the most determinant criterion among the 20 decision layers in delineating the snow avalanche susceptibility maps. Likewise, an increase in elevation leads to an increase in the susceptibility of the region to snow avalanches. The elevation, which was found as the 2nd most important feature in this study, has also been shown among the most determinant factors by the scholars, such that Bian et al. [13] and Iban and Bilgilioglu [14] identified it as the most important factors concerning different topographical and meteorological attributes in different regions, such as the Shaluli Mountains (China) and Province of Sondrio (Italy), respectively.

This study yielded contradictory evidence against the relevant literature regarding the relationship between wind speed and snow avalanche susceptibility. Although the high wind speeds can contribute to the number of avalanches by redistributing snow, creating wind-loaded slopes, and forming dense slab layers, there is a nuanced relationship that needs to be diligently investigated based on site-specific conditions. In this study, an inverse proportional relationship between wind speed and snow avalanches was found in the French Alp, which can be attributed to the interrelationship among wind speed and other determinant factors (e.g., elevation, maximum and minimum temperatures, and snow depth). For instance, windy conditions can be observed in the west, northern, and southern east parts of the focused region (Figure 3l), and in these regions, elevations are relatively lower (Figure 3l), temperatures are relatively higher (Figure 3m,n), and the snow depth is lower (Figure 3p) compared to other parts. The stability of the snowpack is a key factor in avalanche risk, and even in windy conditions, a well-bonded, stable snowpack can result in a low chance of avalanches [85]. The contributions of meteorological inputs on determining the snow avalanche-prone regions since the occurrences of such events can further vary according to the seasonality. In this sense, Yang et al. [86] indicated the changes in the dynamic processes of snow avalanche hazards with regard to meteorological factors, such as temperature, especially concerning the spring season.

Furthermore, regarding the temperature aspect, the increase in both maximum and minimum temperatures results in decreased snow avalanche susceptibility, unsurprisingly (Figure 10). These outcomes are consistent with those provided in the pertinent literature. Such that Iban and Bilgilioglu [14] implemented the SHAP algorithm to delineate the contributions of different factors on the snow avalanche occurrences and found that the hotter a pixel is, the less likely an avalanche may occur. In accordance with the findings of this study, the researchers also found maximum and minimum temperatures among the most influential attributes, ranking 2nd and 5th, respectively. What is also striking about the SHAP outcomes is that the existence of coniferous forests (denoted by LULC_24 according to the CORINE dataset) offered some degree of protection against snow avalanches in the French Alps. This is because coniferous forests, with their dense canopy and vertical structure, can intercept and retain snowfall while providing anchor points for the snowpack, helping to stabilize it and reduce the risk of slab avalanches [87]. From another aspect concerning different LULC characteristics in the focalized region, the presence of natural grasslands (i.e., LULC_26) and bare rocks (i.e., LULC_31) increases the regions’ susceptibility to snow avalanches. Such that natural grasslands tend to provide less protection level against snow avalanches compared to dense forests due to their minimal vegetation, leading to less retaining capability to snow packs [88]. Consequently, by offering less resistance to avalanche flows, such regions may allow snow masses to travel further downslope before dissipating. Likewise, despite a multidimensional relationship between the presence of bare rocks and the occurrences of snow avalanches, bare rocks typically create localized areas where snow can accumulate, forming cohesive slabs that are prone to release as avalanches [86]. Additionally, concave terrain features, such as gullies or couloirs formed by bare rocks, can channelize avalanche movement and amplify its destructive potential. Hence, the corresponding outcomes with regard to the alterations in LULC characteristics of the regions interested in underscore the importance of diligent investigations to assess the probability of snow avalanche occurrences.

7. Concluding Remarks, Limitations, and Implications

Snow avalanches, as one of the most disruptive natural disasters, lead to serious outcomes in living life and a sustainable environment with various facets. These challenges from such damaging incidents make detailed examination crucial. With past records of snow avalanches now more accessible, recent literature has shifted towards using data-driven techniques. As a result, researchers have adopted various methods, including machine learning, which relies on the availability of relevant avalanche data. Therefore, the present research underscored the use of advanced methodologies for exploring the snow avalanche susceptibility levels in the French Alps. To realize this objective, not only the performance comparison of different ML algorithms was conducted, but also the integral role of the meta-heuristic optimization attempts in tuning the hyperparameters of the employed ML techniques was investigated. To identify the most effective predictive framework for generating snow avalanche susceptibility maps in the French Alps, this study conducted a comprehensive comparison of nine different hybrid models with three ML models (i.e., SVC, SGB, and KNN) and three meta-heuristics (i.e., PSO, GSA, and CS). Among the ML techniques adopted, the results highlighted the predictive capacity of the SGB-based models, outperforming both SVC- and KNN-based hybrid structures. Regarding the comparison made within meta-heuristics, the ML models tuned by the PSO and CS algorithms yielded similar outcomes, whereas those optimized through the GSA provided lower accuracy. The predictive ability is not the only concern of this research; rather, the computational efficiency of the adopted techniques quantified with the training times is also considered. According to the overall evaluations, the SGB optimized via the PSO (i.e., PSO-SGB) provided the most satisfactory classification performance with regard to the divergent indicators, including accuracy, precision, recall, F1-score, and AUROC.

This study also addressed a longstanding challenge in data-driven methods: their black-box nature. To overcome this, the game theory-based SHAP algorithm was used to interpret the underlying structure of the best-performing model. Accordingly, the model-agnostic SHAP technique was integrated into the PSO-SGB, and the explainability of the predictions was augmented with a SHAP summary plot, ensuring a well-rounded graphical representation of the interactions between predictors and snow avalanche points. On the one hand, the results showed the contribution of certain geomorphological and meteorological conditioning factors in assessing snow avalanche susceptibility levels. In this regard, slope was found to be the most influential criterion in detecting avalanche-prone regions, followed by elevation, wind speed, maximum and minimum temperatures, and solar radiation. On the other hand, the model underestimated the contribution of most of the lithology and land use/land cover classes (except the presence of coniferous forests lowering the susceptibility of regions to snow avalanches).

Although the present attempt focused on a comprehensive comparison of various predictive frameworks in assessing the snow avalanche susceptibility in the French Alps, it still has some points that need to be urged upon in future studies. For instance, to elucidate the global applicability of the proposed model, its effectiveness can be examined with respect to different regions across the globe. In addition, the utilization of ML algorithms can pose one of the limitations of this research, paving the way to investigation of the performance of deep learning techniques in follow-up attempts, especially for large-scale investigations requiring extensive datasets. Likewise, using a similar approach, different aspects of the model interpretations can be questioned by generating SHAP dependence, force, and decision plots. Additionally, not only different techniques in dealing with the black-box nature of the ML algorithms can be employed, but also inferences with regard to the extraction of the triggering factors of snow avalanche events can be compared to acquire far-reaching outcomes. Overall, the findings reported within this study provided important insights to decision-makers in the region of interest regarding the detection of potential snow avalanche incidents and assisted policy-makers in fostering the mitigation measures to deal with their adverse consequences.