Multi-Hazard Exposure Mapping Using Machine Learning Techniques: A Case Study from Iran

Abstract

Mountainous areas are highly prone to a variety of nature-triggered disasters, which often cause disabling harm, death, destruction, and damage. In this work, an attempt was made to develop an accurate multi-hazard exposure map for a mountainous area (Asara watershed, Iran), based on state-of-the art machine learning techniques. Hazard modeling for avalanches, rockfalls, and floods was performed using three state-of-the-art models—support vector machine (SVM), boosted regression tree (BRT), and generalized additive model (GAM). Topo-hydrological and geo-environmental factors were used as predictors in the models. A flood dataset (n = 133 flood events) was applied, which had been prepared using Sentinel-1-based processing and ground-based information. In addition, snow avalanche (n = 58) and rockfall (n = 101) data sets were used. The data set of each hazard type was randomly divided to two groups: Training (70%) and validation (30%). Model performance was evaluated by the true skill score (TSS) and the area under receiver operating characteristic curve (AUC) criteria. Using an exposure map, the multi-hazard map was converted into a multi-hazard exposure map. According to both validation methods, the SVM model showed the highest accuracy for avalanches (AUC = 92.4%, TSS = 0.72) and rockfalls (AUC = 93.7%, TSS = 0.81), while BRT demonstrated the best performance for flood hazards (AUC = 94.2%, TSS = 0.80). Overall, multi-hazard exposure modeling revealed that valleys and areas close to the Chalous Road, one of the most important roads in Iran, were associated with high and very high levels of risk. The proposed multi-hazard exposure framework can be helpful in supporting decision making on mountain social-ecological systems facing multiple hazards.

1. Introduction

Data and reports on Nat. Hazards show an increasing frequency of such disasters worldwide [1]. In addition, the number of people affected and the economic losses resulting from natural disasters have increased worldwide except in some developed countries, which have experienced a decreasing trend since 1900 [2,3]. Globally, communities in general are facing a wide range of natural hazards. These hazards cannot be eliminated totally, but their effects can be minimized by applying new methods and providing data support for the decision-making process [4,5].

Mountainous areas are characterized as social-ecological systems and, because of their geotechnical and hydrological characteristics, are associated with elements of danger. For example, mountainous regions are relatively active hydrologically and geophysically and there is considerable topography-driven variation in vegetation, moisture, and energy [6,7]. Mountainous environments commonly experience a wide range of natural disasters, such as avalanches [8,9,10], landslides [11,12,13,14], floods [5,15,16,17,18], debris flows [19,20,21], soil erosion [22,23,24], and rockfalls [25,26,27]. Therefore, planners need hazard susceptibility maps to cope with natural disasters in mountainous areas, particularly flash floods, avalanches, and rockfalls [8,25,28]. These three most common hazards can damage transportation systems, property, and danger human life in risk areas [28]. However, even in areas with a wide range of natural risks, previous studies have mainly focused on one type of hazard [8,29,30,31,32]. As a result, integrated multi-hazard analysis has been under-emphasized in natural disaster management and planning. More importantly, most previous studies have focused solely on forecasting and controlling Nat. Hazards [33], neglecting the human effects and the exposure component [4]. In fact, the exposure of mountain social–ecological systems to geohazards has increased because of a combination of the dynamic bio-geophysical environment and the challenges of land use planning [34,35]. Separate spatial studies of different risk types in mountainous regions could in fact result in greater underestimation of risks in total. A multi-risk assessment technique could thus be very useful for considering the interactive reactions of different hazards on surrounding areas [36,37].

Multi-hazard exposure mapping is an important key tool in enhancing resilience in the context of Nat. Hazards in mountain areas and can play an important role in ensuring the long-term sustainability of local ecosystems and communities [38,39,40]. Resilient social–ecological systems have the ability to adapt and utilize relevant knowledge in order to mitigate natural disasters by self-organizing and developing safe infrastructures based on nature-based solution [41]. Fortunately, there is growing recognition of the importance of multi-hazard exposure analysis, and hazard-based attitudes and policies are slowly changing to risk-based management in different parts of the world [1,42,43,44]. In an early study, Bell and Glader [4] developed a general method for analyzing multi-hazards in mountainous areas for multiple processes, focusing on avalanches, rockfalls, and debris flow hazards. In a later study, Gruber and Mergili [36] performed a regional-scale analysis of multi-hazards and risk indicators in high mountainous regions of Tajikistan. Schmidt et al. [45] introduced a generic framework for multi-risk modeling developed in New Zealand and developed a prototype software that is capable of “plugging in” various natural hazards. Recently, Emmer [24] provided a review on natural disasters worldwide and found that research on hydro-meteorological disasters prevails over the research on geological and geomorphic disasters. However, their framework does not provide any guidance on producing a single natural hazard map.

Avalanches, rockfalls, and floods frequently occur in mountainous regions of Iran and cause severe damage to buildings, people, and natural environments [11,12,46,47]. This study investigates multi-hazard exposure using different machine learning approaches due to the following reasons: (1) Machine learning (ML) is a subfield of artificial intelligence where models can learn and improve themselves based on historical events; (2) ML models can easily identify trends and patterns in a large volume of data and involve continuous improvement during operation, which lets them make better decisions [48]; (3) they are capable of handling data that are multi-dimensional and multi-variety [30]; and (4) they can directly extract knowledge of natural disaster processes based on previous disaster occurrences and geo-environmental factors without human intervention, thus, they do not need experts’ experiences and judgements to determine the importance of predictive variables [49]. To evaluate this framework in an area with multiple hazards, the Asara watershed (Iran) was selected because it has become considerably disaster-prone in the past decade and has experienced a disproportionately high number of catastrophic natural disasters. In work to develop a global methodology for multi-hazard exposure mapping in mountainous regions based on historical disaster events, the novelty of this study lies in constructing an integrated framework based on different state-of-the-art machine learning algorithms and using several morphometric, geo-environmental, and topo-hydrological factors for a multi-hazard exposure analysis.

Specific objectives of the study were to: (1) Produce an exposure map based on proximity to residential areas, main roads, rural roads, and power transmission lines; (2) apply three advanced machine learning models including BRT (Boosted Regression Trees), GAM (Generalized Additive Model), and SVM (Support Vector Machine) to produce a multi-hazard map (avalanche, flood, and rockfall hazards); (3) evaluate the performance of each model using threshold-dependent and threshold-independent methods; and (4) produce an integrated multi-hazard exposure map. The intention was for the practical framework of multi-hazard exposure mapping developed in the study to be transferable to other mountainous areas around the world.

2. Study Area

The study area, the Asara watershed, is in Alborz province, which is located in the Alborz mountains of Iran (Figure 1). The watershed occupies an area of 1094.9 km², with a range of elevation from 1332 m.a.s.l at the watershed outlet to 4323 m.a.s.l at the headwaters. Mean annual precipitation in the study area is 265 mm, of which 40% falls as snow. Chalous Road, which is one of the most important transfer lines to the north of Iran, passes through the Asara watershed with 79.54 km. Traffic on the road is generally very heavy, especially on the weekends and on state holidays.

In a land use perspective, the watershed is almost entirely covered by rangeland (95%), but it has 46 villages with around 5000 residents and many restaurants (near the main road). However, only 0.33% of the watershed is residential area. In this mountainous watershed, there are markedly high frequencies and magnitudes of natural hazards, but the level of exposure and risk are also increasing due to socio-economic changes such as the increase of population [50,51]. Some rivers and watercourses in this semiarid region generate flashy stormflows that rapidly peak following brief and intense storms which cause flashy floods in the area. The geological structure and geomorphological characteristics of the study area is considerably complicated which often prone to rockfall and snow avalanche events occurring [52].

3. Methodology

In this section, we provide a step-by-step account of our machine learning approach to creating a multi-hazard exposure map for the study watershed, using an exposure map and a multi-hazard map. A flowchart of the methodology is presented in Figure 2.

The first step of the study is to record the location of disaster occurrences (snow-avalanches, rockfalls, and floods) over previous years (2010–2018). Next, the hazard susceptibility map of each disaster was created and then they were integrated using a weighted linear combination (WLC) technique (i.e., the same weights were considered for all hazard maps) to produce a multi-hazard map. The accuracy of each hazard map was evaluated using statistical criteria including the area under the receiver operating characteristic curve (AUC) and true skill score (TSS). To prepare a multi-hazard exposure map, an exposure map should be created.

3.1. Multi-Hazard Modelling

3.1.1. Predictive Factors

The first step of the work was hazard modeling for avalanches, rockfalls, and floods in the Asara watershed, beginning with selection of predictive factors. However, there are no universal guidelines on choosing influential factors for each natural disaster. Therefore, based on the relevant literature [5,47,53,54,55,56,57,58,59], different topo-hydrological, geo-environmental, and morphometric factors were selected for testing (Figure 3).

A digital elevation model (DEM) was produced for the study area from the airborne LIDAR (Light Detection and Ranging) system that is effective and reliable means of terrain data collection in large areas with cloudy weather conditions. The generated triangular irregular network (TIN) was converted to an ArcGIS grid of 1 m pixel resolution using a TOPOGRID algorithm. In order to acquire homogeneity with other predictive maps, and to increase the efficiency in terms of manipulation and storage, the generated high-resolution DEM was coarsened to 10 m resolution. Scaling from high- to coarse-resolution images is a common approach in spatial analysis and modeling, especially in large-scale regions, and often results in decreasing computational load and time [10,33,44,45]. Therefore, a spatial resolution of 10 m was used to calculate topography related predictive factors including topographic position index (TPI), terrain ruggedness index (TRI), topographic wetness index (TWI), length-slope (LS), relative slope position (RSP), vector ruggedness measure (VRM), aspect, slope degree, distance from stream, and profile curvature. TPI factor can be computed as Equation (1):

$$\mathrm{TPI}=\frac{E_{pixel}}{E_{surrounding}}$$

(1)

where E_pixel and E_surrounding are the altitude of the pixel and the mean altitude of the neighbor pixels, respectively. TWI is a topo-hydrological factor and reflects the wetness potential of each pixel. It can be calculated as below (Equation (2)):

$$\mathrm{TWI}=ln\left(\frac{\alpha }{tan\beta }\right)$$

(2)

where a is an area which drains to a point and β depicts the slope degree at the pixel. TRI also shows the topographic nature of a watershed and can be calculated as Equation (3):

$$\mathrm{TRI}=\sqrt{\left|x\right|\left(ma{x}^2-mi{n}^2\right)}$$

(3)

where x shows the elevation of each neighbor cell to cell (0,0) (m). In addition, min and max reflect the smallest and largest elevation value among nine neighbor pixels, respectively. RSP is also represented as a ratio, ranging from 0 (valley floor) to 100 (ridge top). It can be defined as Equation (4):

$$\mathrm{RSP}=\left(\frac{z\left(s\right)-z{\left(s\right)}_v}{z{\left(s\right)}_r-z{\left(s\right)}_v}\right)\times 100+0.5$$

(4)

where z(s), z(s)_r, z(s)_v are elevation, ridge elevation, and valley elevation, respectively. It is worth noting that TPI, TWI, TRI, and RSP layers were produced in SAGA-GIS software.

In addition, a geological map with 1:100,000 scale was obtained from the Iranian Department of Geology. The details of each geological class described are summarized in

Table 1. Geology of the study area.

Code	Formation	Description	Area (%)
COm	Mila	Dolomite platy and flaggy limestone containing trilobite; sandstone and shale	2.12
Eabtd	-	Andesitic and basaltic volcanic tuff deictic	5.34
Ek	Ziarat	Reef-type limestone and gypsiferous marl	0.03
Ek	Karaj	Well bedded green tuff and tuffaceous shale	53.20
Ek.a	Asara	Calcareous shale with subordinate tuff	6.41
Ekgy	-	Gypsum	5.04
Eksh	Shale	Greenish-black shale, partly tuffaceous with intercalations of tuff	7.96
Jl	Lar	Light grey, thin—bedded to massive limestone	2.03
Mm,s,l	-	Marl, calcareous sandstone, sandy limestone and minor conglomerate	0.72
P	-	Undifferentiated Permian rocks	0.87
pC-Cs	Soltanieh	Thick dolomite and limestone unit, portly cherty with thick shale intercalations	1.06
pCk	Kahar	Dull green grey salty shales with subordinate intercalation of quarzitic sandstone	3.13
Pr	Ruteh	Dark grey medium—bedded to massive limestone	2.81
Qft2	-	Low level pediment fan and valley terrace deposits	1.24
TRe	Elikah	thick bedded to pinkish shaly limestone with worm tracks and well to thick	0.75
TRJs	Shemshak	Dark grey shale and sandstone	7.30

. In addition, Landsat-8 Operational Land Imager (OLI) images were used to create land use map in 2018. Atmospheric corrections were carried out using the FLAASH module in ENVI software. Furthermore, the Gram-Schmidt Pan Sharpening module was applied for fusion of panchromatic and multispectral satellite images. In accordance with the resolution of other DEM-extracted layers, the produced land use map was resampled with a spatial resolution of 10 m in ENVI software. Next, a supervised classification method with the maximum likelihood algorithm was employed. In order to evaluate the accuracy of the produced land use map, not only field investigations with a global positioning system (GPS) receiver were conducted, but also points from Google Earth were also used.

According to confusion matrices, the overall accuracy was 92.40%, which indicates a high classification accuracy.

According to the literature, different predictive factors were selected for each disaster type, as indicated in

Table 2. The spatial database of predictive factors.

Factor	Natural Hazard Type			Data Type and Scale
	Snow-Avalanche	Rockfall	Flood
Elevation	✓	✓	✓	Grid (10m)
Topographic position index (TPI)	✓	✓	✓	Grid (10m)
Terrain Ruggedness Index (TRI)	✓	✓	✓	Grid (10m)
Topographic wetness index (TWI)	-	-	✓	Grid (10m)
Length-slope (LS)	✓	✓	-	Grid (10m)
Relative slope position (RSP)	✓	✓	-	Grid (10m)
Vector Ruggedness Measure (VRM)	✓	✓	-	Grid (10m)
Wind exposition index (WEI)	✓	-	-	Grid (10m)
Aspect	✓	✓	-	Grid (10m)
Slope degree	✓	✓	✓	Grid (10m)
Distance from stream	-	-	✓	Grid (10m)
Profile curvature	✓	✓	✓	Grid (10m)
Lithology	✓	✓	✓	Polygon (1:100,000)
Land use	✓	✓	✓	Polygon (1:100,000)

[26,27,53,54,55,56,57,58,59]. As can be seen in Table 2, some factors are relevant for all hazards.

3.1.2. Multi-Hazard Inventory

Based on field surveys in different seasons and available documents at the Road Organization (RO) of Iran and the Regional Water Company (RWC), the locations of 292 disaster events involving snow avalanches (58), rockfalls (101), and floods (133) have been recorded (2010–2018). Unfortunately, there is no information on disasters before 2010. The RWC have been using Sentinel-1-based processing and ground-based information for flood detection and monitoring. This method has different advantages, as have been described in the literature [9,12]. It operates at specific wavelengths not limited by a lack of illumination or cloud cover. In addition, it can acquire data over a site under all weather conditions during night or daytime. Moreover, it is suitable for gathering information in flooded areas because of its high revisit frequency (i.e., high temporal resolution).

For each disaster type, the collected observed data were randomly divided into two groups, one for model training (70%) and one for model validation (30%), as suggested elsewhere [30,60]. The locations of each hazard type and some field photographs of the damages of hazards are shown in Figure 1. As can be seen, field photographs of the study watershed clearly indicate the sites of avalanches, rockfalls, and floods and the consequences and damage.

3.1.3. Application of Models

Generalized Additive Model (GAM)

GAM is a non-parametric regression method proposed by Hastie and Tibshirani [61], which is a combination of the Generalized Linear Model and an additive model. Basically, GAM is formed by two functions, a link function and a smoothed function, in which the first one is used to build the relation of the mean of the dependent variable with the smoothed function of the independent variables. Thus, incorporating the two mentioned functions enables the GAM not only to work well with various distribution types of data, but it also has ability to deal with highly non-linear problems.

Consider a dataset D = $\left({\mathit{x}}_i,y_i\right)$, ${\mathit{x}}_i\in {\mathrm{R}}^{\mathrm{n}}$ and $y_i$ is the output class i = 1, …, cl, the general form of the GAM is depicted in the following equation:

$$G\left(y\right)= \alpha +f_1\left(x_1\right)+\dots +f_m\left(x_m\right)$$

(5)

where g(.) is the link function, f(.) is the smooth function; $x_i$ is the independent variable; y is the expected value; and $\alpha$ is the intercept.

The process of training the GAM model can be separated into three main steps. The first one is to carry out the data analysis to identify the distribution curve for the dependent variable and the link function, whereas the second one is to determine the parameters of the GAM model, and the last one is to optimize the model.

Boosted Regression Tree (BRT)

BRT has been proposed by Friedman [62]. It is relatively new and powerful machine learning algorithm for regression and classification, which combines a regression tree and a boosting algorithm. The motivation of the BRT is that the performance of the regression tree as a single base learner is not satisfactory, but the performance can be boosted with a series of base learners [63].

Using the dataset D, the boosting algorithm improves the regression tree model, $F\left(x\right)$ by adding an estimator, $h\left(x\right)$ to derive a new BRT model, $F_{new}\left(x\right)$ as in Equation (6). This is an iteration process and the number of iterations (M) plays an important role in the performance of the final BRT model

$$F_{new}\left(x\right)= F\left(x\right)+\gamma h\left(x\right)$$

(6)

where $\gamma \in \left(0,1\right)$ is the learning rate which is used to control the problem of over-fitting.

The loss function is established as below:

$$L=\frac{1}{2} {\left[y-F\left(x\right)\right]}^2$$

(7)

With each iteration, a new tree add to the original model must ensure the reduction of the loss function. The training phase of the BRT will be completed when the predefined number of iterations is reached. Other detailed descriptions of the BRT can be found in Elith et al. [64] and Chung [65].

Support Vector Machine (SVM)

SVM is one of the most powerful supervised machine learning algorithms for both classification and regression, which was developed in the framework of statistical learning theory [66]. The working principal of the SVM model is that, first, it transforms the original input data to a new feature space, and then, in the new space, a decision function is established through an optimal hyper plane. The decision function (Equation (8)) will be used to classify new data:

$$g\left(x\right) = sign({\displaystyle \sum }_{i=1}^n{y}_i{\alpha }_{i}K\left({\mathit{x}}_i,{\mathit{x}}_j\right)+b)$$

(8)

where ${\alpha }_i$ is the Lagrange multiplier; $K\left({\mathit{x}}_i,{\mathit{x}}_j\right)$ is the kernel function. In this research radial basis kernel function (RBF, Equation (9)) was used:

$$K\left({\mathit{x}}_i,{\mathit{x}}_j\right)=exp\left(-\gamma \parallel {\mathit{x}}_i-{\mathit{x}}_j\parallel {}^2\right)$$

(9)

The detail explanations of the SVM are given by [67]. All three models were performed in R software. Each disaster (i.e., avalanche, rockfall, flood) was modeled separately using each of these models.

3.1.4. Accuracy Assessment

The goodness-of-fit and predictive performance of each model were assessed using the training group and validation group of data, respectively. In order to quantitatively evaluate the accuracy of the models in the training and validation steps, threshold-independent and -dependent methods were applied.

Among threshold-dependent methods, the true skill score (TSS) was selected because it is a useful evaluation metric for assessing the accuracy of the model [68]. TSS is an alternative for the kappa, while still keeping all its advantages. Similar to the ROC curve, TSS is calculated based on Sensitivity and Specificity as follows (Equation (10)):

$$TSS=\mathrm{Sensitivity}+ \mathrm{Specificity}-1$$

(10)

Both Sensitivity and Specificity factors can be computed as below:

$$\mathrm{Sensitivity}=\frac{TP}{TP+FN}$$

(11)

$$\mathrm{Specificity}=\frac{TN}{FP+TN}$$

(12)

where FN (false negative) and FP (false positive) are the numbers of pixels erroneously classified, whereas TN (true negative) and TP (true positive) are the number of pixels that are correctly classified [46]. All FN, FP, TN, and TP are the fundamental elements of the contingency table (also termed the confusion matrix).

The receiver operating characteristic (ROC) curve is a plot of the Sensitivity versus “1—Specificity” of a diagnostic test. The area under the ROC curve (AUC) is the most commonly used method for scrutinizing the performance of the model and it has been known as a threshold-independent evaluation metric [69,70]. AUC values vary from 0 to 1, where a value of 0 reflects a completely inaccurate test and a value of 1 implies a perfectly accurate test [49]. According to the common classification of the model accuracy, an AUC under 0.7 is unacceptable, 0.7–0.8 is considered good, 0.8–0.9 is considered very good, and more than 0.9 is excellent [14].

3.2. Exposure-Related Factors

Main road network, rural road network, residential areas, and power transmission lines in the study watershed were defined as vulnerable components. The Euclidian distance tool was used to prepare maps of exposure-related factors in the ArcGIS environment. Based on the analytical hierarchy process (AHP) method [71], the weight of each factor was determined. In this step, questionnaires were prepared and distributed to ten experts (hydrologists, geologists, etc.) within and outside Iran. In order to assess the consistency of the comparisons, the consistency ratio (CR) should be calculated. A CR value less than 0.1 indicates that comparisons are consistent and acceptable. In addition, a rate (R) was assigned to each class of exposure-related factor, according to the order of their importance and exposure, using an expert knowledge-based method [71]. Following Adiat et al. [1], R values of 1, 2, 3, 4, and 5 were assigned to classes indicating, respectively, very low, low, medium, high, and very high exposure. These rates were then normalized for use in the spatial analysis of exposure.

3.3. Multi-Hazard Exposure Mapping

In order to generate the final multi-hazard exposure map, avalanche hazard (AH), rockfall hazard (RH), and flood hazard (FH) maps with the highest accuracy from the three machine learning models were selected. Next, a multi-hazard map was produced by combining the most accurate AH, RH, and FH maps (multi-hazard = AH + RH + FH). Finally, the exposure and multi-hazard maps were fed into a risk equation (Risk = Hazard × Exposure) [71,72] and a multi-hazard exposure map for the Asara watershed was generated.

4. Results

4.1. Hazard Maps

Hazard maps for the Asara watershed based on the BRT, GAM, and SVM models were created for avalanches (Figure 4), rockfalls (Figure 5), and floods (Figure 6). According to Figure 4, the snow avalanche susceptibility maps of the study area generated by the three models had the same overall spatial pattern, with high avalanche susceptibility in the central part of the study area. However, the spatial detail of each model also differs: High susceptible class in the GAM result covered more areas than other two models. It is clear that the high susceptibility is mostly in avalanche slopes in the northern aspect, which is in agreement with the results of Kumar et al. [43] who conducted avalanche susceptibility mapping in Nubra valley region, Indian Himalaya. In fact, the shady slopes become instable in the winter, while the sun facing aspects remain more stable.

As can be seen in Figure 5, BRT and SVM models resulted in relatively similar pattern, although the predicted susceptibility degree in the GAM model is a bit different. The northern part of the study area has the highest snow avalanche susceptibility, which also has the highest slope degree. In the flood susceptibility mapping, all the three models showed the relatively same prediction (Figure 6): Flood susceptible areas are mostly in areas with high drainage density. In general, high drainage density sites can be found in different parts of the study area.

Table 3. Goodness-of-fit and predictive performance of models for three hazard types.

Hazard Type	Model	Goodness-of-Fit (Training Step)		Predictive Performance (Validation Step)
		AUC (%)	TSS	AUC (%)	TSS
Snow-avalanche	BRT	91.2	0.73	89.6	0.71
	GAM	84.4	0.68	83.1	0.61
	SVM	94.1	0.82	92.4	0.72
Rock-fall	BRT	92.6	0.75	88.5	0.67
	GAM	92.2	0.72	90.3	0.69
	SVM	95.3	0.84	93.7	0.81
Flood	BRT	96.5	0.83	94.2	0.8
	GAM	89.7	0.72	86.9	0.67
	SVM	91.9	0.75	89.3	0.69

summarized the results of the goodness-of-fit for each disaster type and also each model. In snow avalanche modeling, SVM had the highest accuracy (AUC = 94.1%, TSS = 0.82), followed by the BRT (AUC = 91.2%, TSS = 0.73). However, the least accuracy belonged to the GAM (AUC = 84.4%, TSS = 0.68). According the accuracy of the rockfall models in the training step, SVM outperformed other models based on the both AUC (95.3%) and TSS (0.84) criteria. Considering the AUC criterion, the variation in the goodness-of-fit of BRT (92.6%) and GAM (92.2%) models was not high. TSS values for BRT and GAM indicated that BRT (TSS = 0.75) performed better than GAM (TSS = 0.72). In flood hazard modeling, BRT produced a better prediction (AUC = 96.5%, TSS = 0.83), followed by SVM (AUC = 91.9%, TSS = 0.75), and GAM (AUC = 89.7%, TSS = 0.72). However, the goodness-of-fit of the model illustrates how well the model fit the training data set. The prediction capacity and generalization abilities of the model cannot be judged using the goodness-of-fit of the model because it is determined by the training data set that were already used to build the model [13].

The validation results showed that the SVM model had the highest performance for predicting snow avalanche (AUC = 92.4%, TSS = 0.72), followed by BRT (AUC = 89.6%, TSS = 0.71) (Table 3). However, GAM was the least accurate of the three with an AUC value of 83.1% and a TSS value of 0.61. For rockfall mapping, SVM was the best model based on the AUC (93.7%) and TSS (0.81); GAM (AUC = 90.3%, TSS = 0.69) and BRT (AUC = 88.5%, TSS = 0.67) were second and third, respectively. In addition, the results of validation in flood hazard modeling demonstrated that BRT (AUC = 94.2%, TSS = 0.8) outperformed SVM (AUC = 89.3%, TSS = 0.69) and GAM (AUC = 86.9%, TSS = 0.67).

Considering the best model for each hazard type (i.e., SVM for snow avalanche and rockfall, and BRT for flood), a weighted integration method was used to conduct the multi-hazard susceptibility analysis which the AUC values in the validation step are the assigned weight. The multi-hazard susceptibility map of the study area is shown in Figure 7.

4.2. Exposure Map

Figure 8 represents exposure factors and

Table 4. Normalized rates (NR) for classes of the vulnerability factors.

Vulnerability Factors	Class	Rate (R)	Normalized Rate (NR)
Distance from residential area (m)	0–50	5	5/15 = 0.333
	50–100	4	4/15 = 0.266
	100–250	3	3/15 = 0.2
	250–500	2	2/15 = 0.133
	>500	1	1/15 = 0.066
	Total	15	-
Distance from main road (m)	0–50	5	0.416
	50–100	4	0.333
	100–250	2	0.166
	>250	1	0.083
Distance from rural road (m)	0–50	5	0.555
	50–100	3	0.333
	>100	1	0.111
Distance from power transmission lines (m)	0–50	5	0.625
	50–100	2	0.25
	>100	1	0.125

illustrates the classes of each exposure factor and their normalized rate (NR) values. Results of the AHP based exposure analysis showed that the distance from the residential area (DFRA) is the most important factor, whereas the distance from power transmission lines (DFPTL) had the least importance among exposure-related factors (

Table 5. A matrix of pair-wise comparisons of four criteria for determining weights in AHP method.

	DFRA	DFMR	DFRR	DFPTL	Weight
DFRA	1	3	5	7	0.546
DFMR		1	5	6	0.304
DFRR			1	3	0.099
DFPTL				1	0.055

DFRA: distance from residential area; DEMR: distance from main road; DFRR: Distance from rural road; DFPTL: Distance from power transmission lines.

). It is worth noting that the consistency ratio (CR) was 0.07 in this study. Therefore, it shows that there is a reasonable level of consistency in the pairwise comparison because it is less than 0.1. Based on the maps of exposure-related factors (main roads, rural roads, residential areas, and power transmission lines), the final exposure map was produced (Figure 9). It indicates that residential areas and roads are the most important elements in terms of the exposure analysis.

4.3. Multi-Hazard Exposure Map

Figure 10 represents the multi-hazard exposure map in the Asara watershed. Based on the multi-hazard exposure classification (equal interval method), the watershed was classified into five classes: ery low risk, low risk, medium risk, high risk, and very high risk. There is significant spatial variation in multi-hazard exposure. High risk areas cover 1.88% and very high cover 0.08% of the study area, whereas 52.12%, 36.27%, and 9.65% of study watershed were classified as very-low, low, and medium risk areas, respectively.

5. Discussion

The results of this study indicated that the predictive ability of models in Nat. Hazards is highly affected by the type of the hazard. In this study, we found that SVM was the most accurate model for analyzing avalanche and rockfall hazards in the Asara watershed, while BRT had the highest accuracy for spatially predicting flood hazard. Therefore, a machine learning model that has a high predictive performance for a natural disaster (e.g., snow avalanche) may showed a lower performance for another natural disaster (e.g., flood). The SVM-based classifiers automatically determine classification areas by statistical analysis of observed training datasets. In addition, the SVM model does not offer direct access to classification borders and thus it is not easy to check the classification, so users should carefully select positive and negative events. The results obtained in the present study confirm previous findings on the ability of SVM to predict avalanche hazards [73,74] and rockfall hazards [25]. The results of this study confirm previous findings that BRT is one of the most accurate models for mapping flood-prone areas [57]. Tree-based models automatically fit non-linear relationships to flood area input factors [75]. Therefore, since Nat. Hazards present a high level of complexity, especially in a complex system such as a mountainous watershed, multi-hazard analysis should be undertaken using different advanced machine learning techniques [17,76]. ML models are critical in helping landscape managers to understand the complex factors controlling Nat. Hazards and more efficiently plan the mitigation measures [71,76].

Production of an exposure map based on AHP was one of the important steps in our multi-hazard exposure mapping approach. Analytical hierarchy process offers a flexible, step-by-step, transparent way to combine factors that play an important role for exposure in different dimensions [77,78].

The present study identified the possible high-risk areas in one of the most important mountainous watersheds (in terms of travel and traffic) in Iran. In general, the results should not be viewed as definite risk maps, but should be used as basis for planning risk management and implementing ecological engineering measures to the high-risk areas. The multi-hazard exposure map (Figure 10) showed that areas associated with a high and very high levels of risk represented a very small proportion of the Asara watershed. However, the areas with high and very high risk levels were mostly located in the valleys which are also the concentration points of infrastructure and population.

These are the areas where the preventive actions should be prioritized. In case of floods, it is necessary to protect lower lands that are exposed to flood hazards by constructing impounding dams upstream of the larger villages [79]. The dams will not only prevent flooding but also provide water during the dry season. For villages that are located headwaters and are less exposed to flood risk, a construction of flood protection walls should reduce the risk considerably. To reduce the rock-fall risk, structural protection measures in the areas with high risk, existing infrastructure, and residential areas need to be implemented. These measures can include the removal of unstable rock elements, slope reshaping or stabilization of unstable rock masses by the use of vegetation or steel wire mesh. New buildings should be carefully planned in all high risk multi-hazard areas. They should be sited and designed according to the character of the local landscape and very high risk areas should be avoided [80]. Moreover, all aforementioned engineering measures should be accompanied with raising the people’s awareness in the high risk areas.

6. Conclusions

A multi-hazard exposure study was conducted based on three natural disasters (avalanches, rockfalls, and floods) in a mountainous region in Iran. The SVM, BRT, and GAM models were used to predict disaster prone areas for a very important Iranian watershed with a high level of vehicular traffic. Results indicated that the SVM model had the highest accuracy for avalanches and rockfalls hazards, while BRT demonstrated the best performance for spatially predicting flood hazard. Therefore, the hazard models applied differed in their accuracy of prediction for the different hazards, which should be borne in mind when selecting a model for risk mapping. Furthermore, results revealed that hazard type can affect the ranking of the models. The approach presented and the results obtained could be very useful for multi-risk mapping in mountainous areas around the world. The multi-hazard exposure analysis helps to spatially prioritize the areas of implementation of most cost-efficient engineering and nature-based solution practices in high risk areas. Therefore, it helps the communities to adapt to natural disasters. Further research should be done to investigate the importance of predictive factors in multi-hazard modeling.