Wildland Fire Susceptibility Mapping Using Support Vector Regression and Adaptive Neuro-Fuzzy Inference System-Based Whale Optimization Algorithm and Simulated Annealing

Abstract

Fires are one of the most destructive forces in natural ecosystems. This study aims to develop and compare four hybrid models using two well-known machine learning models, support vector regression (SVR) and the adaptive neuro-fuzzy inference system (ANFIS), as well as two meta-heuristic models, the whale optimization algorithm (WOA) and simulated annealing (SA) to map wildland fires in Jerash Province, Jordan. For modeling, 109 fire locations were used along with 14 relevant factors, including elevation, slope, aspect, land use, normalized difference vegetation index (NDVI), rainfall, temperature, wind speed, solar radiation, soil texture, topographic wetness index (TWI), distance to drainage, and population density, as the variables affecting the fire occurrence. The area under the receiver operating characteristic (AUROC) was used to evaluate the accuracy of the models. The findings indicated that SVR-based hybrid models yielded a higher AUROC value (0.965 and 0.949) than the ANFIS-based hybrid models (0.904 and 0.894, respectively). Wildland fire susceptibility maps can play a major role in shaping firefighting tactics.

1. Introduction

In the past decades, wildfires and forest fires have been one of the major and most pervasive hazards in destroying natural ecosystems. Every year, millions of hectares of rangelands and forests are destroyed by fire worldwide [1]. Lightning, global warming and climate change, deforestation, land management decisions, insufficient precipitation, hot winds, litter accumulation, and friction between dry litter are among the factors that can cause natural fires in rangelands and forests [2,3]. Major climatic factors contributing to increased forest fires include the spread of drought, rising global temperatures, and increased incidence of foehn winds [4,5]. Every year, approximately 4 million hectares of forests worldwide are damaged by fire [6]. Damage caused by fire can be minimized by identifying areas with high fire susceptibility, implementing fire prevention, and taking fire safety measures [7]. To identify fire-susceptible zones, it is necessary to determine the factors affecting the occurrence of fire, such as fuels, topographic and climatic conditions, and human factors [8]. Consequently, the relationship between these factors and the occurrence of fire must be determined. For this purpose, it is necessary to record areas where fires occurred in the past and then match them with the layers of the factors affecting fire susceptibility to determine the relationship between them [9].

Statistical modeling is performed using parametric or nonparametric methods [10,11]. There are several default parameters and limitations in modeling the relationships between variables in parametric models. Considering a default distribution such as a normal distribution for the response variables, linearity of the proposed relationship, equality of error variance, and independence of the variables are among the limitations of parametric models [12,13]. If the actual data do not satisfy the default conditions, it will not be possible to use these models in practice, or their application will be accompanied by considerable error [12]. In addition, none of the parametric models can model complex and nonlinear relationships [14]. The limited amount of input data and the sensitivity of most parametric models to missing data or outliers are among the other limitations [15]. Consequently, parametric modeling techniques such as linear combination and logistic regression used in previous studies lack the necessary precision and efficiency because of the complexity of the relationship between forest fire locations and factors influencing them [8,16]. Therefore, precise and new nonparametric models for predicting and spatial mapping of forest fire susceptibility, including artificial neural networks or data mining methods, as well as machine learning such as support vector machines (SVMs), can better model the complex relationships between factors influencing fire occurrence. To date, many models have been developed in this category and have been used in spatial and non-spatial fields [17,18,19,20,21,22,23].

The adaptive neuro-fuzzy inference system (ANFIS) and SVM are two of the most powerful and successful models used in modeling various phenomena, such as forest fires and flood susceptibility [24,25,26,27,28]. The present study used ANFIS and SVM combined with the two meta-heuristic methods, i.e., the whale optimization algorithm (WOA) and simulated annealing (SA), to prepare forest fire susceptibility maps for the Jerash Governorate (Jordan). A total of four hybrid models (SVM-WOA, SVM-SA, ANFIS-WOA, and ANFIS-SA) were obtained, compared, and evaluated. To the best of our knowledge, no such study with these models has been conducted to date for forest fire susceptibility mapping. Such studies are important in a country like Jordan, where the percentage of cultivated land to the total area of the Kingdom is only about 3%, and the crops affected by fires, especially in the summer and fall seasons of every year, are considered a major source of income for farmers in the region affected by the fires.

2. Study Area and Data

2.1. Study Area

Jerash Governorate is located in the northwestern part of the Hashemite Kingdom of Jordan and is the smallest Jordanian governorate among all 12 governorates with an area of 420 km². It is approximately 45 km north of Amman, the capital. Jerash city (32°16′ N, 35°53′ E) is the center of the Jerash Governorate (Figure 1). The total population of the governorate is 262,100 people, with a population density of 639.6 people/km².

Jerash is one of the richest regions in terms of biodiversity, especially in the field of agriculture in Jordan, where many crops, especially olives, are cultivated. It also cultivates fruits, grains, and field crops, and has pastures and natural forests. The climate of Jerash is a Mediterranean climate, which is cold to moderate in winter and hot in summer, but it is locally considered as one of the most temperate climates in Jordan, so that the temperature in summer may reach 40 °C, while snow may cover many of the heights of Jerash in some winters. Jerash is considered a mountainous region in Jordan, with altitudes ranging from 300 to 1247 m above sea level.

This rich biodiversity in the study area leads to multiple human activities, and this is a major cause of fires. Hikers, tourists, shepherds, and local people who move in the area intensively, and some of them deal with dry weeds and crops during the harvest season in a manner of indifference and lack of adequate awareness of the basics of environmental conservation, which leads to the occurrence of fires.

The majority of these fires occur in the summer and autumn seasons, and the number of fires increases significantly in the years in which the winter season is later than September to October or November in some years, as dry weeds and crop residues remain after the end of the harvest season without any care from the field owners.

2.2. Data Description

To determine the degree of susceptibility to wildfire in the area being studied, data were collected from past wildfire outbreaks and factors that act as predictors. The selection of factors affecting forest fire occurrence is based on a review of the literature, as well as the characteristics of the study area and data availability. The General Directorate of Civil Defense of Jerash Province reported 109 incidents related to fire, and 109 additional randomly sampled points were created outside the fire zones. This resulted in a broad range of points that were separated by 50 m or more, so that no neighboring points were chosen in the same pixel. The random generation and distribution of sampling points created a model of fire and non-fire incidents (Figure 2). Once the sampling points were agreed, they were divided into training and testing datasets to facilitate model training and validation [29]. Wildfire susceptibility mapping was used to calculate the success rate values using the 70% training dataset (76 fire locations) and the 30% validation dataset comprising 33 fire locations [30].

When assessing and mapping forest fires, various factors provide quantitative tools that determine different elements, such as elevation, aspect, slope, land use, rainfall, wind speed, temperature, the normalized difference vegetation index (NDVI), soil texture, radiation, population density, distance to roads, topographic wetness index (TWI), and distance to drainage systems.

2.2.1. Elevation

The study area elevation ranges between −227 and 1225 m above sea level (Figure 2a). The topographic features were derived from a 30-m open-source digital elevation model (DEM). This DEM was downloaded from the International Agriculture Research Consortium for Spatial Information Shuttle Radar Topographic Mission (SRTM) data [31] (http://earthexplorer.usgs.gov, accessed on 25 December 2018). The DEM data were processed using the ArcGIS 10.8 tool. Elevation is among the most critical parameters in most spatial modeling and analyses [32]. Fire behavior is affected by a number of topographic variables. Erten et al. [33] and Ireland and Petropoulos [34] stated that one of the key factors that determines the seriousness and size of fires is elevation. Ramos-Neto and Pivello [35] added that, at high altitudes, fires spread quickly because flames are fanned by the wind, weeds tend to dry quickly in such places, and fire fighters encounter problems when trying to access these areas. Elevation explains the shifts in temperature that follow the adiabatic lapse rate at which temperature falls with altitude, particularly where the land’s surface properties and material features are multifaceted.

2.2.2. Slope

The intensity of a fire is also affected by slope, as steep slopes result in greater preheating of fuels, and an ever-increasing, faster rate of spread as the fire moves up-slope. In addition, areas with higher slopes can only be tackled using complicated and challenging methods of fire control [33,36]. The slope position determines the likelihood of a fire-moving uphill or retreating downhill, and backing fires are more commonly found on lower slope positions, spread more slowly, and have shorter flame lengths. In addition, slope position and elevation are associated with different types of vegetation, as upper slopes contain thickened woody plants with small, hard leaves, whereas the lower slopes contain significant amounts of medic vegetation. Figure 2b shows the slope map obtained from the DEM.

2.2.3. Aspect

Aspect plays a key role in the amount of solar radiation in an area, as well as moisture availability, both of which directly influence the behavior of fires, and have an indirect effect because they determine the type of vegetation found in a location, as well as its density. The aspect map was derived from the DEM of the area and the results were subsequently classified and divided into nine categories: “north, northeast, flat, west, northwest, east, southeast, south, and southwest”.

2.2.4. Wind Speed

Landform trends and wind speed also play a major role in shaping firefighting tactics and the occurrence of fires. Wind speed is also a very important factor in the creation and spread of fires.

2.2.5. TWI

The topographic elevation, aspect, and slope all determine the amount of solar radiation reaching a particular location. Conversely, TWI is a factor related to the structure of the landscape and offers an explanation of the topographic states and settings that play a part in shaping the pattern of the spatial soil and creating its texture [36]. As the catchment area expands and the slope gradient falls, the values of TWI increase in parallel. The increase in the TWI values concurrently raises the moisture levels in the soil and bolsters the greenness of the vegetation, thereby causing a drop in the intensity and range of the forest fire.

2.2.6. Distance to Drainage

The vegetation close to the drainage networks also holds more moisture, and this serves to make the plants very green, and in the process acts as a brake on its spread. The distance from the drainage networks was calculated using Euclidean distance and generated in ArcGIS 10.8.

2.2.7. Temperature

Eight rain gauges and meteorological stations, with records spanning a quarter of a century, supplied the climatic variables, namely average yearly rainfall, average air temperature, average solar radiation, and average wind speed, using the “inverse distance weighted (IDW)” interpolation method. In this way, raster data layers of the particular indices in the GIS environment were produced (Figure 2). Temperature impacts atmospheric air mass, which dictates the relative humidity, air mass, and moisture content of the soil. As the temperature rises, so does the tendency for fire rage and this is particularly true if the increase in temperature is linked to a fast, dry wind. Conversely, the relationship between radiation and temperature is positive, and radiation increases the likelihood of fire accidents, particularly in regions where the land is covered by dry grass and field crops.

2.2.8. Radiation

Radiation differs from place to place, and these variations are determined by a range of factors, such as texture, soil type, and degree of humidity.

2.2.9. Rainfall

Clearly, precipitation results in higher levels of moisture in the soil and encourages the growth of green vegetation, both of which limit the intensity and spread of forest fires [37].

2.2.10. Population Density and Distance to Roads

The majority of wildfires are caused by human actions, including agricultural pursuits, hunting, or any acts that allow or depend on naked flames touching inflammable woodland biomass [38]. Human presence and activity in forest areas will have a major influence on the likelihood of wildfires breaking out, and we therefore recommend combining two features: population density and the distance to a road to tackle anthropic parameter factors. While the wildfire area will be at risk from nearby populations and their habitats, the existence of local roads will help to flag and supply paths that can be used for fire suppression responses.

2.2.11. Land Use Land Cover

Combustibility factors refer to the ability of vegetation to spread fire and burn up fuels, and its characteristics therefore reflect the degree to which it is consumed by heat masses [39]. The land use/cover of the area plays a part in promoting fire to various degrees, depending on the make-up of the cover and the amount of human contact. For example, the risk of fire rises when agricultural land is situated close to forests or orchards whose undergrowth and dry trunks can be used as fuel, thus promoting fire spreading. The Landsat 8 Operational Land Imager (OLI) imagery data for the study area on 25 December 2018, were downloaded from the United States Geological Survey (USGS) website [40] (http://glovis.usgs.gov/, accessed on 25 December 2018). The image was categorized into nine LULC classes: bare soil, water bodies, vegetation, rural, urban, soil rocks, rocks, valley, and sand (Figure 2g), using the maximum likelihood classifier algorithm [41]. An error matrix was created for accuracy assessments by comparing the results of remote sensing analysis with reference sample points [41]. A random number generator was used to yield the random x and y coordinates of the total number of 196 reference sites within the study area. All locations were then either visited in the field or assessed using the Google Earth map service system. The results of the accuracy assessment showed that the overall accuracy was 85.5%, and the Kappa coefficient was 0.79.

2.2.12. NDVI

NDVI measures the surface cover and density of vegetation, which is an essential piece of information because it highlights the amount of potential fuel available for feeding and spreading fire.

2.2.13. Soil Texture

Knowing the soil texture type is useful for uncovering the kind of vegetation that can be planted and plays a major role in drawing up the land use/cover map, which is why a soil texture map of the study area was compiled. The soil map of the study area was produced by the Ministry of Agriculture (MoA) at a scale of 1:250,000. The map was digitized and processed using ArcGIS 10.8. The soil texture of the area was divided into four categories: loam, silt loam, clay loam, and silty clay loam.

3. Methodology

This section describes the methods used to build the proposed hybrid models as below:

3.1. ANFIS

ANFIS was first used by Jang in 1993 [42]. This model has been widely employed in various scientific fields. ANFIS has the advantages of both fuzzy inference systems and neural networks. In other words, ANFIS can be considered a powerful and efficient model for predicting and modeling because it allows the modelers to use adaptive training algorithms to train fuzzy system parameters, thus enabling them to exploit the advantages of neural networks and fuzzy systems [42].

The Takagi–Sugeno fuzzy system, which has a linear relationship with its output and parameters that can be estimated by combining the least squares with backward propagation of errors based on gradient descent, is the most common fuzzy inference system capable of being used in an adaptive network [43]. In these systems, the structure of a model is selected first with defined factors that are proportional to the inputs, degree of membership, and rules [44]. A part of the dataset was then selected and entered into the training phase. In this phase, the values of the parameters and factors in the model must be close to their actual values by minimizing the amount of error. Data not used in the training phase were used to validate the model.

Figure 3 presents the architecture of the ANFIS model as the structure of a progressive network with five layers with two inputs, one rule, and one output. In the first layer (input), the degree of belonging of each input to various fuzzy intervals is determined by the user. Multiplying the input values of each node with each other yields the weight of the second layer. The relative weights of the rules were calculated for the third layer.

The fourth layer is the layer of rules in which the relative weights of the rules (w_i) are computed. The last layer is the network output, which aims to minimize the difference between the output obtained from the network and the actual output.

3.2. SVR

This model was extracted from Vapnik’s statistical learning theory for working with regression problems and estimating and predicting data [46]. When using the SVR model to solve nonlinear problems, a kernel function must first be defined to map the input data into a higher-dimensional space [47]. In this study, the radial basis function (RBF) kernel function was used owing to its high performance in previous research [48,49]. Equation (1) shows this function.

$$k\left(x,y\right)=\exp \left(-\frac{1}{{\gamma }^2}\parallel x-y\parallel 2\right)$$

(1)

where x, y, and γ are the input vectors of factors, the response or dependent variable, and the RBF kernel parameter, respectively. The problem of SVR model optimization is solved using Equation (2) [46]:

$$\begin{gathered} \frac{1}{2}{w}^{T}w+C\sum _{i=1}^n({\xi }_i+{\xi }_i^{*}) \\ \mathrm{Subject} \mathrm{to}:\left\{\begin{array}{c}{y}_i-\left(w^T{x}_i+b\right)\le \epsilon +{\xi }_i\\ \left(w^T{x}_i+b\right)-y_i\le \epsilon +{\xi }_i^{*}\\ {\xi }_i,{\xi }_i^{*}\ge 0\end{array}\right. \end{gathered}$$

(2)

In Equation (2), $\epsilon$ is the deviation from the hyper-plane, ${\xi }_i,{\xi }_i^{*}$ are noise levels of the classes, C is a constant, and w is the vector orthogonal to the hyper-plane. By solving the above equation by employing the Lagrangian method and using the Karush–Kuhn–Tucker (KKT) conditions, the following equation is obtained:

$$f\left(x\right)={\displaystyle \sum }_{i=1}^n({\alpha }_i-{\alpha }_i^{*})K\left(x_i\right)+b$$

(3)

In this equation, ${\alpha }_i \mathrm{and} {\alpha }_i^{*}$ are Lagrange multipliers, and b is derived from Equations (4) and (5):

$$b=y_i-w^T{x}_i-\epsilon \mathrm{for} 0<{\alpha }_i<C$$

(4)

$$b=y_i-w^T{x}_i+\epsilon \mathrm{for} 0<{\alpha }_i^{*}<C$$

(5)

The quality of the SVR model depends on selecting a suitable kernel function and the proper setting of its parameters. In the presented equations and the RBF kernel function, the C and γ parameters must be tuned. The C parameter, which is used to prevent overfitting, is called the regulation parameter. The γ parameter controls the degree of nonlinearity of the RBF kernel function.

3.3. WOA

Optimizing a system means minimizing or maximizing a function that is a measure of system performance [50,51,52]. Optimization methods have been used in various fields [53,54,55,56]. One of the most important functions of these methods is to determine the optimal values of the parameters in various mathematical problems. The WOA is a novel metaheuristic algorithm for solving optimization problems that mimics the hunting behavior of humpback whales [57]. A species of whales, the humpback whale, utilizes a specific hunting mechanism known as the bubble-net (Figure 4). In this strategy, the whales dive to the 12-m depth and blow bubbles in a shrinking circle around the prey, trapping them. The whales then swim upward and reduce the bubble-net diameter, ultimately hunting a group of fish. Inspired by the spiral bubble net, the authors in [57] presented the WOA meta-heuristic algorithm.

Similar to other population-based algorithms, this algorithm includes a set of solutions: the exploitation phase, where the whales surround the prey once spotted. The algorithm considers the position of the best solution as far as the prey location. Therefore, whales move to the best location using Equation (6) [57]. D in Equation (6) is also calculated using Equation (7).

$$\overrightarrow{X}\left(t+1\right)=\overrightarrow{X^{*}}\left(t\right)-\overrightarrow{A}.\overrightarrow{D}$$

(6)

$$\overrightarrow{D}=\left|\overrightarrow{C}.\overrightarrow{X^{*}}\left(t\right)-\overrightarrow{X}\left(t\right)\right|$$

(7)

$$\overrightarrow{A}=2\overrightarrow{a}.\overrightarrow{r}-\overrightarrow{a}$$

(8)

$$\overrightarrow{C}=2.\overrightarrow{r}$$

(9)

The term $\overrightarrow{X}\left(t+1\right)$ denotes the whale position in the next iteration, $\overrightarrow{X^{*}}$ is the position of the best solution, and $\overrightarrow{A}$ and $\overrightarrow{C}$ are the coefficient vectors. The symbol || returns the absolute value, and the dot product is denoted by a dot. In Equation (8), $\overrightarrow{r}$ is a random vector within the range [0,1], and $\overrightarrow{a}$ in Equation (8) decreased from two to zero in each iteration. In Equation (6), which updates the position of whales in each iteration, the $\overrightarrow{A}$ and $\overrightarrow{C}$ vectors can be adjusted to allow whales to move to different positions around the best solution. The hunting strategy of humpback whales is composed of two phases. They reduce the size of the circle around the prey and swim upward in a spiral path. To model the shrinking circle, the parameter $\overrightarrow{a}$ is considered in Equation (8), which reduces from two to zero during the different iterations of the algorithm. Consequently, parameter $\overrightarrow{A}$, which assumes a value within the range [−a, a], is reduced.

To model the spiral motion, the distance between the whale and prey (the best solution so far) was first measured, and then the spiral motion of whales was simulated using the spiral equation:

$$\overrightarrow{X}\left(t+1\right)=\overrightarrow{D^{\prime }}.e^{bl}.\cos \left(2\pi l\right)+\overrightarrow{X^{*}}\left(t\right)$$

(10)

where $\overrightarrow{D^{\prime }}=\left|\overrightarrow{X^{*}}\left(t\right)-\overrightarrow{X}\left(t\right)\right|$ is the distance between the whale and pray, b is a constant that determines the shape of the logarithmic spiral, $l$ is a random number in the range [−1,1], and is the dot product. The probability $p$ is considered in the algorithm because the shrinking circles and traveling of the spiral path are simultaneously carried out by the whales. The search agent chooses to perform a shrinking circle or spiral motion with a 50% probability. Hence, the position of the search agents is updated using Equation (11) [57]:

$$\overrightarrow{X}\left(t+1\right)=\left\{\begin{array}{cc}\overrightarrow{X^{*}}\left(t\right)-\overrightarrow{A}.\overrightarrow{D}\hfill & if p<0.5\hfill \\ \overrightarrow{D^{\prime }}.e^{bl}.\cos \left(2\pi l\right)+\overrightarrow{X^{*}}\left(t\right)\hfill & if p\ge 0.5\hfill \end{array}\right.$$

(11)

Exploration phase: For the algorithm to be efficient and converge to a global optimum, both the exploitation and exploration phases should be utilized. In this algorithm and for the exploration phase, instead of selecting the best solution (search agent), the search agents randomly select another search agent and move toward it. Vector $\overrightarrow{A}$ was used for this purpose. In case $\left|\overrightarrow{A}\right|<1$, the search agent moves toward the best solution, and if $\left|\overrightarrow{A}\right|\ge 1$, it moves toward a random solution (Equations (12) and (13); [57]):

$$\overrightarrow{X}\left(t+1\right)=\overrightarrow{X_{rand}}-\overrightarrow{A}.\overrightarrow{D}$$

(12)

$$\overrightarrow{D}=\left|\overrightarrow{C}.\overrightarrow{X_{rand}}-\overrightarrow{X}\right|$$

(13)

The WOA can be generally described as follows: an initial population of search agents is randomly initialized, and the best solution is determined by assessing the solutions. The search agents then attempt to improve the solutions based on Equation (14), which expresses the overall structure of the model.

$$\overrightarrow{X}\left(t+1\right)=\left\{\begin{array}{c}shirink circle\left\{\begin{array}{cc}move toward best solution \left(eq.1\right)& if \left|\overrightarrow{A}\right|<1\hfill \\ move toward random solution \left(eq.7\right)& if \left|\overrightarrow{A}\right|\ge 1\hfill \end{array}\right. if p<0.5\\ move along spiral shape path\left(eq.5\right) if p\ge 0.5\end{array}\right.$$

(14)

The search agents continue this process until the termination criterion is met, at which point the global optimum is probabilistically determined by the algorithm.

3.4. SA

SA is one of the most popular meta-heuristic algorithms used for optimization problems [58]. This algorithm simulates the process of melting and cooling metals to alter their physical properties (referred to as annealing in the field of metallurgy). At high temperatures, the metal atoms move freely, but once the temperature drops, the movement of atoms is restrained and they are held together in a regular crystalline structure with a minimum energy level, providing the metal with a high physical strength [58,59]. These modifications prepare the metal for high workloads.

In optimization problems, the objective is to find the solution minimizing/maximizing the target function. For instance, in the minimization case, the SA algorithm regards the minimization of the objective function as minimizing the energy level in the annealing process. In the controlled cooling process, if a better solution exists in the vicinity of the current solution in terms of minimizing the objective function, it is selected by the algorithm; otherwise, the algorithm selects the worst solution by evaluating a set of criteria. This approach allows the algorithm to escape local minima [58,59]. The probability of selecting the worst solution is calculated as follows [59]:

$$p=\exp \left(\frac{-\Delta E}{K.T}\right) pϵ\left[0,1\right]$$

(15)

$$\Delta E=E_{next state}-E_{current state}$$

(16)

where T is the temperature, $\Delta E$ is the change in energy level, and K is a constant. At the beginning of the search, when the temperature is high, the algorithm is more likely to accept worse solutions to explore different parts of the search space [60]. As the temperature drops, the search space is gradually reduced, and the tendency of the algorithm gravitates more toward accepting better solutions and disregarding the worse ones, allowing it to focus on the better regions of the search space [60]. As the process continues, the algorithm ultimately converges to a good solution. Hence, the probability of selecting worse solutions helps the algorithm to search for different parts of the search space without being trapped in local optima, ultimately converging to the global optimum [58,59].

This algorithm can solve different optimization problems that are regarded as unsolvable by traditional methods. As its most important advantage, this algorithm escapes local minima and tends to find the global optimum. However, the optimization duration can be long for large-scale problems [61,62].

3.5. Hybrid Models

ANFIS and SVM are among the most widely used and powerful models in machine learning and data mining, and both have been used in various spatial and non-spatial problems. However, both algorithms need to tune their parameters to obtain accurate and powerful models. In SVR, the model’s performance is highly dependent on the C and gamma parameters, and it needs to be fine-tuned. In the ANFIS algorithm, two sets of parameters, used in the adaptive layers of the model (the first and fourth layers in Figure 2), require fine-tuning. In this study, we used the SA and WOA metaheuristic algorithms to tune these parameters and obtain accurate models.

3.6. Frequency Ratio (FR)

FR determines the quantitative relationship between the occurrence of fires and each class of variables [63]. The ratio of fire occurrence in each class is an influential factor in relation to the total number of fires. The FR was calculated using Equation (17):

$$FR=\frac{\left(\raisebox{1ex}{$A$}\!\left/ \!\raisebox{-1ex}{$B$}\right.\right)}{\raisebox{1ex}{$C$}\!\left/ \!\raisebox{-1ex}{$D$}\right.}$$

(17)

where A is the number of fires occurring in each class, B is the total number of fires, C is the number of pixels in each class, and D is the total number of cells in the study area.

3.7. Feature Selection Process

The selection of effective features in modeling is one of the most important steps [64]. Before modeling, feature selection was performed to eliminate irrelevant factors based on the two models of ANFIS and SVR and meta-heuristic algorithms of SA and WOA. For feature selection using the SA algorithm, we defined an initial solution as a binary vector $F=〈f_1, f_2,\dots , f_{14}〉$, $f_i\in \left\{0, 1\right\}$, where zero means that the variable is not present in the model, while one means that the variable is used in modeling. The initial solution was randomly generated and assumed to be the optimal solution. The neighboring solutions are created using the initial solution by applying a slight change in the vector, for example, by changing one or two bits. In the next step, the solutions are evaluated using a fitness function. In this study, the SVM and ANFIS models were created using the variables corresponding to each solution. Then, the RMSE metric was used as the cost function to evaluate these models, and the algorithm tries to select the features that minimize this metric. In this process, the three major parameters of the SA algorithm, namely, the initial temperature, temperature reduction rate (TRR), and termination condition, also need to have appropriate values. Various studies have shown that the initial temperature should be sufficiently high to allow for sufficient transactions to occur. In this study, we set this parameter to 100, while the TRR was set to 0.99. In addition, the termination condition was selected to be after 200 iterations. For feature selection using WOA, as in other meta-heuristic algorithms, we need to define which form the solutions would take. In this population-based algorithm, each solution is modeled as a one-dimensional vector with $N$ binary variables, where zero means that the feature is not present in the model, and one means that the feature is used in modeling. As in the SA algorithm, first, an initial set of solutions (whale population) is generated. Then, the RMSE metric was used to evaluate the fitness values of the SVM and ANFIS models corresponding to each solution. Once the solution corresponding to the best fitness value denoted as $X^{*}$ is considered, the algorithm enters a loop. The whales’ positions, that is, the solutions, are updated at each iteration. To do this, at each iteration, the values of the A, a, C, l, and p parameters (introduced in the equations explained in Section 3.3) are updated. In this algorithm, if A > 1, the exploration phase is executed by searching the neighborhood of a random solution. However, if A < 1, the exploitation phase is executed, and the neighborhood of the current best solution is searched. This process continues until the stop criterion is satisfied. All models were executed using the MATLAB software. In addition, owing to the small size of the training dataset, the implementation of the models was not time consuming.

3.8. Relative Operating Characteristics (ROC)

The ROC [63] was used to evaluate the performance of the models. The ROC curve is a graphical representation of the trade-off between the negative and positive error rates for each possible sample. This metric is a curve whose vertical and horizontal components are calculated using Equations (18) and (19) [65]:

$$X=1-\left[\frac{True Negative \left(TN\right)}{True Negative \left(TN\right)+False Positive \left(FP\right)}\right]$$

(18)

$$Y=\left[\frac{True Positive \left(TP\right)}{True Positive \left(TP\right)+False Negative \left(FN\right)}\right]$$

(19)

The TP, TN, FP, and FN values were obtained from the confusion matrix (

Table 1. Confusion matrix.

		Model
Reality		Positive	Negative
	Positive	True Positive (TP)	False Negative (FN)
	Negative	False Positive (FP)	True Negative (TN)

) with the definition of the threshold between zero and one.

In addition, the precision–recall curve is another evaluation metric for machine learning models. This metric is suitable for small and skewed datasets [66]. Precision and recall are calculated using Equations (20) and (21).

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

(20)

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

(21)

4. Results

Table 2. Spatial Relationship Between each class of the factors and forest fire occurrence by the FR.

Factors	Classes	No. of Pixels	No. of Fires	FR Weights
Elevation (m)	<318	40,604	4	0.57
	318–533	100,277	13	0.75
	533–724	113,580	32	1.63
	724–917	98,876	10	0.59
	>917	94,337	18	1.11
Aspect	Flat	3155	0	0
	North	45,764	4	0.51
	Northeast	40,550	16	2.29
	East	57,810	12	1.21
	Southeast	68,725	5	0.42
	South	86,407	12	0.8
	Southeast	55,591	9	0.94
	West	52,336	14	1.56
	Northwest	37,333	5	0.78
Slope	0–5	81,146	11	0.79
	5–15	267,384	44	0.96
	15–30	92,246	22	1.39
	>30	6895	0	0
Rainfall	300–350	113,527	8	0.41
	350–400	154,167	30	1.13
	250–300	721	0	0
	400–450	145,311	28	1.12
	450–500	33,945	11	1.88
NDVI	−0.46–−0.12	9921	4	2.34
	−0.12–−0.05	30,689	16	3.03
	−0.05–0.00	87,641	22	1.5
	0.00–0.04	164,081	23	0.81
	0.04–0.28	155,339	12	0.45
TWI	<9	121,081	27	1.3
	9–13	217,055	32	0.86
	>13	109,535	18	0.96
Distance to rivers	0–150	133,092	22	0.96
	150–300	107,127	20	1.09
	300–450	87,220	9	0.6
	450–600	60,126	12	1.16
	>600	60,106	14	1.35
Temperature (C)	15.80–17.38	102,494	19	1.08
	17.38–18.48	101,010	10	0.58
	18.48–19.48	122,338	34	1.62
	19.48–20.48	85,621	13	0.88
	20.48–22.50	36,208	1	0.16
Distance to roads	0–150	336,475	67	1.16
	150–300	77,584	7	0.52
	300–450	23,018	3	0.76
	450–600	7010	0	0
	>600	3584	0	0
Wind speed (m/s)	<2	44,027	14	1.83
	2–4	297,820	45	0.87
	4–6	99,681	17	0.98
	>6	2231	1	2.58
LULC	Water	1612	0	0
	Bare soil	159,679	48	1.75
	Vegetation	256,798	25	0.57
	Urban	2020	2	5.76
	Rural	20,120	1	0.29
	Soil rocks	3953	0	0
	Rocks	1049	0	0
	Valley	1296	0	0
	Sand	1144	1	5.08
Radiation (watt/m2)	5.43–5.64	18,272	3	0.95
	5.64–5.70	81,265	17	1.22
	5.70–5.75	180,178	43	1.39
	5.75–5.82	167,956	14	0.48
Soil texture	Loam	4698	1	1.24
	Silty Clay	481	0	0
	Silty Loam	261,932	44	0.98
	Clay Loam	180,560	32	1.03
Population density	0.29–2.36	227,182	31	0.79
(person/km2)	2.36–5.46	124,104	26	1.22
	5.46–9.90	413,38	8	1.13
	9.90–15.28	39,838	11	1.61
	15.28–26.65	15,209	1	0.38

present the results of FR method. In the elevation factor, the class 533–724 m and areas with height below 318 m have the highest and lowest weight, respectively. In the aspect factor, flat class has the least weight, whereas the areas with northeast direction have the greatest impact on forest fire occurrence. Likewise, the values of FR weights for all classes of the factors can be seen in Table 2. After building the training and validation datasets and preparing the factors affecting fire occurrence, the modeling was performed using four hybrid models. Feature selection was performed prior to the modeling procedure to prevent correlations among the factors (

Table 3. Optimal features selected by SVR-SA, SVR-WOA, ANFIS-SA, and ANFIS-WOA models.

Explanatory Variables	SVR-SA	SVR-WOA	ANFIS-SA	ANFIS-WOA
Elevation	☑	☑	☑	☑
Slope aspect	☑	-	-	-
Slope	☑	☑	☑	☑
Land use	☑	-	☑	☑
TWI	☑	☑	☑	☑
NDVI	☑	☑	☑	☑
Distance to drainage	☑	☑	☑	☑
Distance to roads	-	-	-	-
Soil texture	-	-	☑	-
Wind speed	-	☑	-	☑
Solar radiation	☑	☑	☑	☑
Rainfall	☑	☑	☑	☑
Temperature	☑	☑	☑	☑
Population density	☑	☑	☑	☑

). In this study, for feature selection, we followed the approach of [26,67]. The major factors selected by the SVR-WOA model were solar radiation, population density, TWI, distance to drainage, temperature, rainfall, NDVI, land use, slope, and elevation. However, the factors selected by SVR-SA were radiation, population density, TWI, distance to drainage, temperature, rainfall, NDVI, land use, slope, aspect, and elevation. Furthermore, the factors selected by ANFIS-WOA consisted of radiation, population density, TWI, distance to drainage, temperature, rainfall, NDVI, land use, slope, wind speed, and elevation, whereas the ANFIS-GA selected radiation, population density, TWI, distance to drainage, temperature, rainfall, NDVI, land use, slope, soil, and elevation (Table 3).

After the selection of the main factors, they were entered into the ANFIS and SVR models. As discussed earlier, SA and WOA were employed to determine the premise and consequent parameters in the ANFIS model, as well as the C and γ parameters in SVR. Figure 5 and Figure 6 show the real and modeled values of the training and validation data in addition to the RMSE and MSE for the ANFIS-WOA and ANFIS-SA models. Accordingly, the ANFIS-WOA model had a smaller RMSE (0.317) in the training phase than the ANFIS-SA model (0.324). However, in the validation phase, the RMSE of the ANFIS-SA model (0.371) was slightly lower than that of the ANFIS-WOA (0.376).

Figure 7 presents the best fitness (best RMSE) and mean fitness (mean RMSE) for each solution in each generation for the SVR-WOA and SVR-SA models. In the SVR-WOA, 20 generations and a total number of 500 iterations were set. Similarly, 500 iterations are selected for SVR-SA. Figure 7 shows that SVR-SA outperformed SVR-WOA in terms of the best fitness (best RMSE).

The ROC is one of the key metrics used to evaluate the model accuracy in various fields. Figure 8 shows the AUROC values for the four hybrid models in the training and test datasets. Accordingly, the SVR-WOA and SVR-SA models obtained the highest AUROC in the validation phase. In addition, Figure 9 shows the precision-recall curves for the models in the training and test datasets. Accordingly, SVR-based models have higher average precision (AP) than ANFIS-based models.

The AUROC value of SVR-SA in the training phase (0.927) and in the validation phase (0.966) clearly shows that unlike the other three models, this model is less accurate in the training phase than in the validation phase. This can be due to the small size of the training and test samples, as well as the sampling method [27]. After modeling, all cells in the study area with features obtained by feature selection phase were prepared by combined analysis in a GIS environment and were entered into the four models, and the probability of forest fire occurrence for the entire area was calculated. After obtaining the output for all cells, join and look up analyses were used in the GIS environment to prepare forest fire susceptibility maps. Figure 10 shows the wildfire susceptibility maps for the four models. These maps are classified into five classes: very low, low, moderate, high, and very high.

Table 4. Overlap of high and very high susceptibility classes with each land use class.

Land Use Class	Pixel No.
Water	0 (0.00%)
Bare soil	67,614 (55.23%)
Vegetation	50,032 (40.87%)
Urban	387 (0.32%)
Rural	3117 (2.55%)
Soil	552 (0.45%)
Rocks	67 (0.05%)
Valley	143 (0.12%)
Sand	183 (0.15%)

shows the areas occupied by the five classes of the four models. As can be seen in Table 4, the two land use classes of bare soil and vegetation have the highest number of high- and very high-risk cells, and only 0.32% of these two classes occurred in the urban land use class. On the other hand, there were 5785 urban cells in the study area, of which only 387 cells (6.68%) were in high and very high classes. Bare lands are almost covered by dry weeds in the summer and autumn seasons, which make them a target of fires due to human activities in the area, which is consistent with the results of many similar studies.

5. Discussion

Considered important sources of information, wildfire susceptibility maps can help managers and decision-makers to prevent fires or mitigate potential risks. This study proposed four hybrid models, namely SVR-WOA, SVR-SA, ANFIS-WOA, and ANFIS-SA, for wildfire susceptibility in the Governance of Jerash in Jordan. In this study, we concluded that the SVR performed better. Similar results have been reported in other studies [68,69]. Nevertheless, such conclusions may vary from one area to another; therefore, researchers should compare SVR with other models, especially ANFIS, in terms of performance in order to obtain more accurate modeling outputs.

The results obtained in the present study include factors influencing wildfire occurrence in the Jerash Governorate in Jordan, as well as the wildfire susceptibility maps. A study of the important factors identified by feature selection algorithms revealed that temperature, distance to drainage, TWI, population density, solar radiation, elevation, slope, land use, NDVI, and rainfall were the 10 most important factors among the 14 factors. These 10 factors were considered for wildfire susceptibility mapping. Generally, the effect of elevation on fire occurrence is related to a set of other factors that are human factors, such as the development of residential areas, population density, and intensity of land use changes or factors related to the environment and climate, such as rainfall, temperature, moisture, and evapotranspiration [70]. In general, the probability of fire occurrence is much higher at lower to medium elevations in a region compared to the other elevation classes because of their relatively higher temperature, lower moisture, and easier human access. In the study region, all elevation classes (−227 to 1225 m), including forests and agricultural lands, along with the mentioned factors, provide the conditions for fire occurrence by making available a sufficient amount of combustible materials. Many studies have proved the effect of steep slope areas on the occurrence of fire [71]. Nevertheless, the present research indicates that slopes of over 30% exhibit lower susceptibility to fires because they are covered with rocks and are devoid of combustible materials. Moreover, Oliveira et al. [72] showed that areas with high temperatures and low rainfall increased fire risk if combustible materials were available. In the present study, areas with high temperatures and low rainfall had the highest risk of fire. Land use is also considered as one of the most important factors influencing fire occurrence. In fact, the presence of sufficient combustible materials was the main reason for fire occurrence in rangelands, forests, and agricultural lands. In addition to natural fires in these areas, human factors play a substantial role in fire occurrence by causing intentional fires to remove crop residues, repelling harmful and destructive animals in agriculture, and accelerating the regeneration of herbaceous species in the rangelands. It is noteworthy that population density is an important human factor in fire occurrence in the present study. Distance to drainage was also considered an important factor among the 14 studied factors. Some studies have attributed this factor to tourist hubs. However, studies in the region suggested that this factor was not an important factor in the occurrence of fire in the present study. Areas distant from the drainage network faced a high risk of fire, and the solar radiation factor was closely related to forest fuel moisture. In the present study, areas with high solar radiation intensity had a low susceptibility to fires. TWI and NDVI factors influenced fire occurrence which is also corroborated by [73].

The combination of machine learning models with meta-heuristic algorithms has been repeatedly used by researchers. However, the number of studies comparing ANFIS-based hybrid models with SVM-based hybrid models is small. As the SVM algorithm has obtained better results in many studies [68,74], the development of hybrid models with this algorithm and its comparison with other algorithms such as ANFIS can help modelers and planners in choosing the optimal model. Therefore, the present study showed that the combined SVR-WOA and SVR-SA algorithms are more efficient in forest fire susceptibility mapping.

Overfitting occurs when the learning performance is very good, but the performance on the validation dataset is poor. In addition, under-fitting occurs when learning performance not good and it is also not good on the validation data set. Therefore, as can be seen from Figure 8, these two problems did not occur for all models. The AUROC in the training and validation phases for all models had close values. For the SVR-SA model, the AUROC in the validation phase was greater than that in the training phase, but there was not much difference between them. [27] showed that there are three reasons for this case: (1) feature selection process, (2) small size of training and test data sets, and (3) the sampling technique used for training and validation datasets. For more details, please refer to [27].

6. Conclusions

This study employed four hybrid models, namely SVR-WOA, SVR-SA, ANFIS-WOA, and ANFIS-SA, to analyze forest fire susceptibility in the Jerash Province of Jordan by providing forest fire probability maps. Fourteen parameters were selected as factors affecting the occurrence of fire in the study area. Among the considered variables, solar radiation, population density, TWI, distance to drainage, temperature, rainfall, NDVI, land use, slope, and elevation were identified by all hybrid models as variables affecting the occurrence of fire. The results indicated that there is a small difference in the AUC-ROC between the SVR-WOA, SVR-SA, ANFIS-WOA, and ANFIS-SA models. However, according to the research findings, the SVR-based hybrid algorithms outperformed the ANFIS-based hybrid algorithms. FFSMs can be used to help managers control and prevent fires. In high-risk areas, measures such as building water tanks, prohibiting construction, and building emergency aid stations can be adopted to improve the situation.