Spatial Prediction of Landslide Susceptibility Using GIS-Based Data Mining Techniques of ANFIS with Whale Optimization Algorithm (WOA) and Grey Wolf Optimizer (GWO)

Abstract

The most dangerous landslide disasters always cause serious economic losses and human deaths. The contribution of this work is to present an integrated landslide modelling framework, in which an adaptive neuro-fuzzy inference system (ANFIS) is combined with the two optimization algorithms of whale optimization algorithm (WOA) and grey wolf optimizer (GWO) at Anyuan County, China. It means that WOA and GWO are used as two meta-heuristic algorithms to improve the prediction performance of the ANFIS-based methods. In addition, the step-wise weight assessment ratio analysis (SWARA) method is used to obtain the initial weight of each class of landslide influencing factors. To validate the effectiveness of the proposed framework, 315 landslide events in history were selected for our experiments and were randomly divided into the training and verification sets. To perform landslide susceptibility mapping, fifteen geological, hydrological, geomorphological, land cover, and other factors are considered for the modelling construction. The landslide susceptibility maps by SWARA, SWARA-ANFIS, SWARA-ANFIS-PSO, SWARA-ANFIS-WOA, and SWARA-ANFIS-GWO models are assessed using the measures of the receiver operating characteristic (ROC) curve and root-mean-square error (RMSE). The experiments demonstrated that the obtained results of modelling process from the SWARA to the SAWRA-ANFIS-GWO model were more accurate and that the proposed methods have satisfactory prediction ability. Specifically, prediction accuracy by area under the curve (AUC) of SWARA, SWARA-ANFIS, SWARA-ANFIS-PSO, SWARA-ANFIS-GWO, and SWARA-ANFIS-WOA models were 0.831, 0.831, 0.850, 0.856, and 0.869, respectively. Due to adaptability and usability, the proposed prediction methods can be applied to other areas for landslide management and mitigation as well as prevention throughout the world.

1. Introduction

Landslides are some of the most fatal natural disasters worldwide, which greatly influences human life and development of infrastructures in underdeveloped countries [1,2]. For instance, about 15,000 landslides happened in China between 2015 and 2017, which caused the economic loss of about one million US dollars. Therefore, it is urgent to conduct studies on landslide susceptibility mapping (LSM) and it is a primary and practical tool for mitigation and risk assessment [3].

In recent years, the assessment and management of landslides have resulted in significant reduction in losses in several countries [4,5], where correlation on developments of landslide and influencing factors is clarified. Persichillo et al. [6] found that the human activities have some influence on sediment connectivity of landslides and the intensity of their occurrence. Safran et al. [7] inferred that large landslides may locally inhibit incision by changing the river channels’ slope and width and riverbed characteristics. Zhao et al. [8] found that landslide movement is a non-uniform and non-rigid body motion which uses speckle and grayscale feature search methods. Gao et al. [9] found that a rotated ellipsoid trend surface and a random field of residuals were influenced by the spatial maximum rolling rainfall [10]. Rainfall can enlarge the weight of soil, alter the pore water pressure and decrease the strength of soil, and trigger other types of landslides [11,12,13].

Although researchers are making tremendous efforts, prediction of landslide with high accuracy is still challenging, especially in the regional scale [14,15]. The quality of landslide inventory map has improved significantly with the use of new remote sensing techniques. For example, Keyport et al. [16] analyzed pixel-level and object-level landslide detection methods using high resolution remote sensing images, where major landslides were easy to recognize with a few errors. However, quality of the landslide prediction is strongly dependent on the modelling approach [17], therefore, various data-driven approaches have been considered, including logistic regression [18], neural networks [19,20,21,22,23,24,25], support vector machine (SVM) [18,26,27,28,29,30,31], relevance vector machine [32], least squares SVM [33,34], decision trees [35,36,37,38,39], logistic model trees [40], random forest [38,41,42], gene expression programming [43]. However, no method or technique is the best for landslide modelling at the regional scale for all regions.

One of the difficulties in landslide modelling at the regional scale is to handle landslide influencing factors which are usually derived from various sources with different scales that may contain uncertainties and imprecisions [44]. This requires new modelling methods other than traditional data-driven approaches, which have abilities to deal with the above issues and enhance the model performance. Literature review shows that some approaches have been proposed, i.e., fuzzy k-means [45], fuzzy logic [46,47], neural fuzzy [48,49,50,51], and fuzzy k-nearest neighbor [44], and among them, the neural fuzzy, which is a hybrid of artificial neural network (ANN) with fuzzy logic, has proven its efficient for landslide modelling in term of prediction accuracy [52,53]. However, performance of the neural fuzzy model is strongly influenced by its weights for membership function and finding these optimal weights is still a challenging task [54]. Recent developments of machine learning (ML) have introduced new optimization algorithms, which could be used for optimizing weights for membership function of the neural fuzzy model. Furthermore, ML techniques have recently gained a good attention among environmental modeling research community as they are advantageous in efficiently capturing the complex relationship between the environmental predictors and the response, such as flood [55,56,57,58,59,60,61,62,63], earthquake [64,65], wildfire [66], sinkhole [67], droughtiness [68], gully erosion [69,70], groundwater [71,72,73,74] and land/ground subsidence [75], and landslide in this case [54,59,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102]. Nevertheless, investigation of new optimization algorithms and the neural fuzzy has not been carried out.

The aim of this study is, therefore, to patriotically fill this gap in literature through investigating potential application of the whale optimization algorithm (WOA) and the grey wolf optimizer (GWO) algorithm and the ANFIS algorithm for spatial probability of landslide occurrence in Anyang Country, China. The WOA and GWO algorithms have been proven better than other popular meta-heuristic algorithms, such as genetic algorithm (GA), ant colony optimization (ACO), and particle swarm optimization (PSO) for optimization in various real world problems [103,104]. In addition, the step-wise weight assessment ratio analysis (SWARA) method is used to evaluate the relationship between landslide and influencing factors. To our best knowledge, the hybrid framework of SWARA-ANFIS-WOA and SWARA-ANFIS-GWO has not been studied for landslide modelling. The prediction performance of the proposed methods was evaluated using several objective measures of root-mean-square error (RMSE), receiver operating characteristic (ROC) curve, and area under the ROC curve (AUC). Finally, the effectiveness of the proposed methods was assessed by comparing the landslide susceptibility maps with the historical landslide events.

2. Material and Methods

2.1. Introduction of the Study Area

Anyuan County, which has an area of about 2374.59 km², is located in the south of Jiangxi Province, China (Figure 1). Its longitude and latitude lie between 115°9′52′′ to 115°37′13′′ E and 24°52′18′′ to 25°36′52′′ N, respectively. The altitude of this area is from 132 m to 1180 m (msl). The climate of Anyuan County is of a subtropical monsoon region. The average temperature in the period of 1960–2016 and the humidity are 18.3 °C and 76%, respectively. Meanwhile, the annual average sunshine time and rainfall are about 72 and 108 days, respectively, and the average annual rainfall is 1562.1 mm.

There are two major rivers of Zhenjiang and Lianjiang in the Anyuan area and 37 reservoirs. The geological map of the study area demonstrates that more than 29 geologic groups and units can be observed, as shown in Table 1 and Figure 2, and the main lithological categories are carbonaceous shale and coal seam, feldspar quartz coarse sandstone, granite, residual tuff; crystalline chip, diorite. According to the fifth census data of China (http://www.ay.gov.cn/xxgk/tjxx), the population and the Gross Domestic Product (GDP) of Anyuan County are 398,614 people and about 8.7 billion US dollars per year. As for the population, 301,563 and 97,061 people are living in urban and rural areas, respectively.

2.2. Data Preparation

2.2.1. Landslide Inventory Map

The initial step of LSM is to produce an accurate and reliable inventory map, which is always obtained by field survey, image interpretation from remote sensing data and historical landslide records. In this work, the landslide inventory map of Anyuan County with 315 landslide locations provided by the Department of Land and Resources of the Jiangxi Province and the Jiangxi Meteorological Bureau. In this map, the smallest and largest landslides are 2.4 m² and 7,840 m², respectively. The distribution and impact of landslides with different sizes are greatly different. For instance, more than two-thirds of the total landslides are small-scale (<200 m²) landslides, which have influenced 2352 people and cause the economic loss of 0.1 million US dollars, followed by medium-scale (200–1000 m²) and large-scale (>1000 m²) landslides, which have influenced 479 and 55 people and caused the economic loss of 0.3 and 1 million US dollars, respectively. Meanwhile, heavy rainfall is another key factor for landslide occurrence.

2.2.2. Landslide Influencing Factors

In this study, fifteen landslide influencing factors, including altitude, slope, aspect, plan curvature, profile curvature, stream power index (SPI), sediment transport index (STI), topographic wetness index (TWI), lithology, distance to fault, distance to river, land use, normalized difference vegetation index (NDVI), soil and rainfall (Figure 3). A DEM of the Anyuan area was obtained from the ASTER GDEM Version 2 and we can produce the geomorphological factors of slope, altitude, aspect, curvature, SPI, STI and TWI by using ArcGIS and SAGA-GIS software. In Figure 3a,b, altitude was divided into seven categories of 132–282 m, 282–432 m, 432–582 m, 582–732 m, 732–882 m, 882–1032 m and >1032 m, and slope was classified to eight angles of <20.0°, 20.0–24.6°, 24.6–29.2°, 29.2–33.8°, 33.8–38.4°, 38.4–43.0°, 43.0–47.6° and >47.6°. As previous studies, aspect was classified to nine directions in Figure 3c, including east, south, west, north, northwest, northeast, southwest, southeast and flat. Both of the plan and profile curvatures were discriminated into three categories of <−0.001, (−0.001)–(0.001) and >0.001, as shown Figure 3d,e.

The factors of SPI, STI and TWI were divided into seven classes in Figure 3f–h, respectively. As for the geological factors, there are eight groups in the lithological map in Figure 3i, and distance to fault was initialized to six classes of 0–500 m, 500–1000 m, 1000–1500 m, 1500–2000 m, 2000–2500 m, and >2500 m in Figure 3j. As for the hydrological factor of distance to river, it was classified to seven intervals in Figure 3k, i.e., (1) 0–150, (2) 150–300, (3) 300–450, (4) 450–600, (5) 600–750, (6) 750–900, (7) >900. The land use and NDVI can be obtained using a Landsat 7 ETM+ satellite image acquired on 12 October 2001. The values of NDVI are in the range of [−1, 1] and were separated into six classes of <−0.05, (−0.05)–0.05, 0.05–0.15, 0.15–0.25, 0.25–0.35, 0.35–0.45, and >0.45 in Figure 3l. To obtain the land use information, the supervised classification algorithm of maximum likelihood was used with the overall accuracy of 91.2%, and the study area was classified to six objects of bare, forest, grass, residential, farmland, and water in Figure 3m. Another land cover factor of soil was provided by the Institute of Soil Science, Chinese Academy of Sciences (ISSCAS), China, into 5 different types of ATc (Cumulic anthrosols), ACu (Humic acrisols), ALh (Haplic alisols), ACh (Haplic acrisols), and RGc (Calcaric regosols) in Figure 3n. As for the factor of rainfall, the data from 23 rainfall stations for the period 1960–2015 were collected to produce the precipitation map in Figure 3o, and there are six levels at an interval of 100 mm in this map. Finally, all the factor maps were resampled with a spatial resolution of 25 m.

Figure 4 illustrates the proposed framework which mainly includes the four steps as follows. (1) The landslide inventory map of the study area is compiled and the influencing factors are selected; (2) The initial weights of each class of landslide condition factors evaluated by SWARA; (3) The landslide modelling is conducted on the study area using the proposed ANFIS-WOA, and ANFIS-GWO methods. (4) The experimental results are analyzed and evaluated.

2.3. The Initial Weights for Modelling Process

As a multiple-criteria decision analysis (MCDM) method, the step-wise weight assessment ratio analysis (SWARA) algorithm was firstly introduced for assessing the correlation between factors by computing a set of weights for each factor. There are two strategies for the SWARA algorithm. The first strategy is to analyse the different situations and each layer of influential factors (criteria) is prioritized based on the requirements and aims. The other strategy is experts’ knowledge, which is of great importance to the prioritization of each layer of the influencing factors [105]. For clarification, the SWARA algorithm is summarized in two steps:

Step 1: According to the relationship between the influencing factors, a decision-making model based on experts’ judgement is developed and then the criteria are prioritized and sorted in descending order.

Step 2: Compute weights for each factor.

First, each expert judges the prioritization of each criterion for each influencing factor. Then, the comparative importance of the average value (CIAV) is given by [106]:

(1)$IAV=\frac{{{\displaystyle \sum }}_i^N{D}_i}{N}$

where $N$ is the number of experts and $D_i$ denotes the offered ranks by the experts for each factor.

Then, the coefficient $A_i$ is defined as follows:

(2)$A_i=\{\begin{array}{ll}1& j=1\\ CIAV+1& j>1\end{array}$

where j indicate the number of the factors.

The re-calculated weight $D_j$ is computed as follows:

(3)$D_j=\frac{CIAV-1}{A_j}$

Finally, the relative weights of the criteria are represented as follows:

(4)$W_j=\frac{D_j}{{{\displaystyle \sum }}_{j=1}^m{D}_j}$

where $W_j$ denotes the relative weight of the jth criterion and m represents the total of the criteria.

2.4. Landslide Prediction Models

2.4.1. ANFIS

The ANFIS is derived from the ANN and fuzzy logic [107]. The ANN model does not obtain the output data from making decision. Although it has automatic learning capability and the fuzzy logic is inverse of ANN model [108], the ANFIS algorithm has the advantages of both the ANN and fuzzy logic models on producing input and output data in one framework [109], and it is known as a commonly used data-driven model for solving nonlinear issues [110] in the fields of data processing and fuzzy control [111].

This algorithm is based on the Takagi and Sugeno’s type by using two “If-Then” rules as follows [112]:

• Rule 1: If x₁ and x₂ are A₁ and B₁, respectively, then f₁ = p₁x₁ + q₁x₂ + r₁.

• Rule 2: If x₂ and x₂ are A₂ and B₂, respectively, then f₂ = p₂x₂ + q₂x₂ + r₂.

where A₁, A₂, B₁, and B₂, are the fuzzy sets, p_i, q_i and r_i are the parameters, ${\mathrm{x}}_1$ and ${\mathrm{x}}_2$ are input data and $f_1$ and $f_2$ are output data [113]. The ANFIS system includes six layers and the detailed descriptions on these steps can be referred in [112].

2.4.2. Whale Optimization Algorithm (WOA)

The WOA algorithm is a meta-heuristic optimization process that mimics humpback whales [103], whose brain has the same spindle cells in the cortex as humans’ [114]. WOA is inspired from the sole hunting method of colossal humpback whales [103] and is known as the bubble-net feeding method [115]. This method is based on the unique pattern for catching far more fishes at once. A group of whales come together and dive beneath the school of fishes by producing high pitch calls. At that moment, the fishes flee to the surface, where the whales release the distinctive bubbles along a circle of 9-shaped trail in an upward shrinking spiral around the fishes as an obstacle so that the fishes cannot swim, as shown in Figure 5. Finally, the whales ascend to the surface with their mouths open by the helix-shaped movement when the whale leader emits a hunting call [114]. The WOA can be formulated as following steps:

(1) Encircling prey

Encircling prey in the best position is the initial step in hunting prey. The humpback whales look for the suitable position of prey and renew their positions based on the optimal known solution. The solution can be represented by Equations (5) and (6):

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

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

where ${\overrightarrow{X}}^{*}$ and $\overrightarrow{X}$ are two position vectors and ${\overrightarrow{X}}^{*}$ indicates the optimal solution obtained so far, t is the current iteration. Furthermore, $\overrightarrow{A}$ and $\overrightarrow{C}$ specify the coefficient vectors that are given by:

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

(8)$\overrightarrow{C}=2\overrightarrow{r}$

where a is a linearly decreasing variable and r is a random vector in the range from 0 to 1.

(2) Bubble-net attacking method (exploitation phase).

This exploitation phase includes two processes given by:

(a) Shrinking encircling prey

In this mechanism, the values of $\overrightarrow{a}$ and $\overrightarrow{A}$ decrease in Equations (7) and (8). The vector $\overrightarrow{A}$ has a random value between [−a, a] and declines from 2 to 0. The new position can be reached using the original and current optimal agent positions.

(b) Spiral position updating

Let $\left(X,Y\right)$ and $\left(X^{*},Y^{*}\right)$ indicate the positions of the whale and prey and ${\overrightarrow{D}}^{\prime }$ the distance between them; the whales mimic the helix-shaped movement according to spiral equation as follows:

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

where ${\overrightarrow{D}}^{\prime }=|{\overrightarrow{X}}^{*}\left(t\right)-\overrightarrow{X}\left(t\right)|$ is constant, and $l$ is a random variable between [−1, 1], and b is constant. A mathematical model is required for the updating of the whales’ position in optimization process because they swim towards the prey in two mechanisms of shrinking circle and spiral-shaped path. The updating of the whales’ position is formulated as follows:

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

where p is a random variable ranging from 0 to 1.

(3) Search for prey (exploration phase)

This phase can perform a global search when the whales search a randomly chosen agent. This mechanism is used when $\overrightarrow{A}$ has a random value larger than 1 or less than −1. The search for prey can be modelled by the following equations:

(11)$\overrightarrow{D}=|\overrightarrow{C}.{\overrightarrow{X}}_{rand}-\overrightarrow{X}|$

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

where the position vector ${\overrightarrow{X}}_{rand}$ is chosen randomly from the whales between the current population. The flowchart of the WOA algorithm is illustrated in Figure 6.

2.4.3. Grey Wolf Optimizer (GWO)

The GWO algorithm is another meta-heuristic algorithm [104] and is derived from the hunting behavior of the grey wolves and the social hierarchy in nature [116]. The Grey wolves live in pack by a strict social dominant hierarchy and imitate the leadership hierarchy [117].

The grey wolves have various groups for activities including hunting the prey. The social dominant hierarchy of the grey wolves shown in Figure 7 has four forms: alpha (α), beta (β), delta (δ), and omega (ω). The α-wolves are the top of the social hierarchy as leaders for making decisions whereas other wolves follow them [118]. The β-wolves help the leaders and devise them. The δ-wolves situate the next level and obey the α- and β-wolves. Finally, the ω-wolves have to submit to all of them [117]. The GWO technique has an optimization process similar to other meta-heuristic algorithms through the collection of random candidate solutions [119], thus it is designed as a mathematically model for grey wolves. Specifically, α, β, δ and ω are the fittest, second optimal, third optimal and the rest of solutions, respectively [120].

The process of hunting behavior of GWO is shown in Figure 7 and this algorithm includes four phases described as follows:

(1) Encircling of prey

In the first phase, the grey wolves harass and encircle the prey during hunting. The parameter $D$ measures the distance between the grey wolf and the prey and is given by:

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

where t represents the current iteration, and ${\overrightarrow{X}}_P$ and $\overrightarrow{X}$ denote the position vectors of the prey and the grey wolves, respectively, and the coefficient vector is defined by $\overrightarrow{C}$ as follows:

(14)$\overrightarrow{C}=2.{\overrightarrow{r}}_1$

where $r_1$ is a random vector which is in the interval [0, 1]. The prey’s location can be represented as follows:

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

The value of the coefficient $\overrightarrow{A}$ is computed as follows:

(16)$\overrightarrow{A}=2a.r_2-a$

where a is a linearly decreasing variable and $r_2$ is a random vector between [0, 1] as $r_1$.

(2) Hunting

After the phase of encircling the prey, the hunting behavior is guided by α, β and δ, respectively, since they have compressive information about the prey’s position. This behavior is shown as follows:

(17)${\overrightarrow{D}}_{\alpha }=|{\overrightarrow{C}}_1.{\overrightarrow{X}}_{\alpha }-\overrightarrow{X}|,{\overrightarrow{D}}_{\beta }=|{\overrightarrow{C}}_2.{\overrightarrow{X}}_{\beta }-\overrightarrow{X}|,{\overrightarrow{D}}_{\delta }=|{\overrightarrow{C}}_3.{\overrightarrow{X}}_{\delta }-\overrightarrow{X}|$

(18)${\overrightarrow{X}}_1={\overrightarrow{X}}_{\alpha }-{\overrightarrow{A}}_1.\overrightarrow{D,}{\overrightarrow{X}}_2={\overrightarrow{X}}_{\beta }-{\overrightarrow{A}}_2.\overrightarrow{D,}{\overrightarrow{X}}_3={\overrightarrow{X}}_{\delta }-{\overrightarrow{A}}_3.\overrightarrow{D}$

where ${\overrightarrow{X}}_1$, ${\overrightarrow{X}}_2$ and ${\overrightarrow{X}}_3$ indicate the position vectors of α, β and δ, respectively. ${\overrightarrow{A}}_1,{\overrightarrow{A}}_2,{\overrightarrow{A}}_3$ and ${\overrightarrow{C}}_1$, ${\overrightarrow{C}}_2$, ${\overrightarrow{C}}_3$ are the coefficients that can be computed using (14) and (16). The position of a grey wolf in search space can be updated as follows:

(19)$\overrightarrow{X}\left(t+1\right)=\frac{{\overrightarrow{X}}_1+{\overrightarrow{X}}_2+{\overrightarrow{X}}_3}{3}$

The positions of other wolves update randomly according to the position of the prey.

(3) Attacking prey (exploitation)

The process of hunting finishes when the prey stops moving. Mathematically, $\overrightarrow{A}$ is linearly decreased from 2 to 0. The exploration trend happens when $|\overrightarrow{A}|<1$ and $|\overrightarrow{C}|<1$. At this moment, the wolves attack the prey.

(4) Searching of prey (exploration)

The Grey wolves track and chase the prey. The pursing prey is known as the exploration phase in GWO algorithm [121]. The parameters α, β, and δ have the duty of guidance roles in this process. If $|\overrightarrow{A}|>1$, it means the grey wolves diverge and scatter in different directions for searching of the prey. After finding it, they converge to attack [118]. The coefficient $\overrightarrow{C}$ provides a random weight for the prey while $|\overrightarrow{C}|>1$ and promotes the exploration phase. In addition, $\overrightarrow{C}$ models the natural obstacles in hunting for the grey wolves [122]. The flowchart map of the GWO algorithm is shown in Figure 8.

2.4.4. Particle Swarm Optimization (PSO)

PSO as a Meta-heuristic algorithm was introduced by Kennedy et al., 1995 [123,124]. The ability to optimize nonlinear problems and fast convergence as well as a few computations are the most remarkable features of this algorithm. These features make PSO significantly different from other evolutionary algorithms such as genetic algorithms. This algorithm stems from being used by the swarm intelligence of birds and fishes, which they apply it to discover the best means to find food. In this algorithm, each bird is implemented as a particle which in fact represents a solution to the problem. These particles are searched in n-dimensional space, where n is the number of problem parameters, to find the optimum answer for the problem. To this aim, particles are randomly scattered in the search space. Then in each iteration, based on Equations (20) and (21) each particle adjusts its location by finding the best location that it has ever been in and the best one adjacent to its neighbor.

Supposing $x_i^t=\left(x_{i1}^t,x_{i2}^t,\dots ,x_{in}^t\right)$ and $v_i^t=\left(v_{i1}^t,v_{i2}^t,\dots ,v_{in}^t\right)$ are respectively the location and rate of change in location of the i-th particle in t-th iteration, therefore, according to Equations (1) and (2), the location and rate of change in location of the i-th particle in t + 1-th iteration can be determined as follows:

(20)$v_i^{t+1}=\omega .v_i^t+c_1.r_1.\left(p_i^t-x_i^t\right)+c_2.r_2.\left(g_i^t-x_i^t\right)with-v_{max}\le v_i^{t+1}\le v_{max}$

(21)$x_i^{t+1}=x_i^t+v_i^{t+1}$

where, $x_i^t$ is the previous location of i-th particle, $p_i^t$ is the best location found by $g_i^t$ the i-th particle, ، is the best location found by other particles and $r_1,r_2$ are random numbers from 0 to 1. Moreover, the three parameters $\omega$, $c_1$ and C₂ are the cognitive coefficient, social coefficient, and inertia weight, respectively. It should be noted that there are diverse articles in order to set these parameters. In this paper, they are determined based on the following equations:

(22)$\omega =\frac{1}{2ln2}and{c}_1=c_2=0.5+ln2$

It is noteworthy that this algorithm continues until the best location found by each particle is the equal of the best position chosen by all the particles. In other words, all the particles are concentrated in one point of space, and in fact, the answer to the problem is optimized.

3. Results and Validation

3.1. Correlation between Landslides and Influencing Factors Based on SWARA

Table 2 lists the correlation between landslides and influencing factors based on SWARA. In this table, the SWARA weight for slope is highest of 0.19 for the class of >47.6°, indicating the highest possibility of landslide occurring. For aspect, the highest SWARA weight of 0.3 is obtained in the direction of south. In terms of altitude, the class 882–1032 m (0.33) was assigned with the highest SWARA weight of 0.33. The highest SWARA weights of 0.51 and 0.55 for plan and profile curvatures are achieved in the classes of >0.001 and <−0.001, respectively. For the three factors of SPI, STI and TWI, the highest SWARA weights of 0.32, 0.7, and 0.6 are reached in the classes of 800–1000, 30–40 (50–60), and 3–5, respectively. The lithology of the class C resulted in the best SWARA weight of 0.36. For the distance to fault and river, the highest SWARA weights 0.11 and 0.21 are computed for the distances of 0–500 and 600–750, respectively. The classes of forest and 0.35–0.45 for land use and NDVI were given the best SWARA weight of 0.53 and 0.4, respectively. The soil of ACu has the greatest SWARA weight of 0.54 and the SWARA weight for rainfall is highest of 0.21 in the classes of 1500–1600 mm. Also, one of the main advantages of this method is that considered the relative weights of each influencing factor compare to each other. Result shows that slope angle and Plan curvature have a highest and lowest impact on the landslide occurrences respectively (Table 3).

3.2. Application of the Integrated ANFIS Methods

To construct the training and verification sets, the same number of 315 non-landside grid cells were first randomly selected using ArcGIS. Then, 70% landside and non-landside grid cells were used for model building, and the remaining landside and non-landside grid cells were used for assessing the prediction performance [125,126,127,128]. It is better to noted that in this stage, only the weights of each factors by SWARA method is considered for integration. Once the training and verification sets were constructed, the SWARA, SWARA-ANFIS, SWARA-ANFIS-PSO, SWARA-ANFIS-WOA and SWARA-ANFIS-GWO methods were used to predict landslide susceptibilities for Anyuan County. In the proposed framework, the meta-heuristic algorithms were programmed using MATLAB software, and the parameters of membership functions were optimized [129,130,131,132]. The GWO, WOA, and PSO algorithms enhance the prediction accuracy of ANFIS. Firstly, a given weight is assigned to each class and each influencing factor, and then, an output is computed. Consequently, each weighted class is considered as input data for ANFIS to find the weight of each factor. In the next stage, GWO, WOA, and PSO based on ANFIS are implemented to optimize the obtained weight of ANFIS in the earlier stage. Figure 9 shows a mathematic example that how these algorithms run and lead to enhance the power prediction of ANFIS algorithm.

Figure 10, Figure 11 and Figure 12 show the accuracy of training and verification processes for the SWARA-ANFIS-GWO and SWARA-ANFIS-WOA, and SWARA-ANFIS-PSO methods, respectively. The RMSE was selected to assess the prediction performance of the three methods in both training and verification phases. The results demonstrate that the RMSE for SWARA-ANFIS-PSO, SWARA-ANFIS-WOA and SWARA-ANFIS-GWO methods were 0.099, 0.1, and 0.099 in the training phase and 0.120, 0.11, and 0.1 in the verification phase, respectively. Meanwhile, GWO can optimize the better parameters of membership functions than that of WOA in both the training and verification phase. To evaluate the efficiency of the two methods, the cost function values are calculated and plotted versus each iteration of convergence graph for the three proposed methods, as shown in Figure 13. The cost function values of ANFIS-GWO, SWARA-ANFIS-WOA and SWARA-ANFIS-PSO are stable after 17, 362, and 49 iterations, which validates that SWARA-ANFIS-GWO is more efficient than SWARA-ANFIS-WOA and SWARA-ANFIS-PSO.

3.3. Preparation of Landslide Susceptibility Mapping

To produce the final resultant maps, the landslide susceptibility index (LSI) for each pixel was used to indicate the landslide probability occurrence in a certain area. In this work, the LSI measure is estimated for the study area, and the resultant maps were reclassified to five levels using the natural break method [133]. Figure 14 shows five resultant maps by the SWARA, SWARA-ANFIS, SWARA-ANFIS-PSO, SWARA-ANFIS-WOA, and SWARA-ANFIS-GWO methods. Results show that according to the SWARA-ANFIS-WOA model, the moderate class has the largest area (20.80%), followed by low (20.32), very low (19.82%), high (19.68%), and very high (19.38) classes. For the SWARA-ANFIS-GWO model, the percentages are 37.66%, 32.39%, 12.22%, 11.52%, and 6.1% for the very low, low, moderate, high, and very high classes, respectively. For the SWARA model, the percentages are 19.60%, 20.22%, 19.95%, 19.95%, and 20.28% for the very low, low, moderate, high, and very high classes, respectively. For the SWARA-ANFIS model, the percentages are 0.24%, 32.23%, 20.60%, 17.02%, and 29.91% for the very low, low, moderate, high, and very high classes, respectively. For the SWARA-ANFIS-PSO model, the percentage are 15.72%, 22.55%, 22.86%, 18.50%, and 20.37% for the very low, low, moderate, high, and very high classes (Figure 15).

3.4. Validation of the Resultant Maps

To objectively evaluate the proposed methods, the ROC and AUC measures were used. In the ROC curve, the specificity and sensitivity are on the x- and y-axis, respectively [134,135]. Figure 16a,b demonstrate the ROC using the training and verification sets for the SWARA, SWARA-ANFIS, SWARA-ANFIS-PSO, SWARA-ANFIS-WOA, and SWARA-ANFIS-GWO methods, respectively. It can be observed that SWARA-ANFIS-WOA and SWARA-ANFIS-GWO can obtain the AUCs of 0.897 and 0.900 when using the training set and 0.856 and 0.869 when using the verification set. Furthermore, these observations are in accordance with the analysis with the RMSE results. Additionally, some methods and optimization algorithms such as SWARA, SWARA-ANFIS, and SWARA-ANFIS-PSO are used for comparison of the obtained results. Accordingly, results conclude that based on the training dataset the AUC for SWARA, SWARA-ANFIS, and SWARA-ANFIS -PSO are 0.836, 0.894, and 0.894, respectively. However, these results for validation dataset are 0.831, 0.831, and 0.850, respectively. The experimental results demonstrate that both the two proposed new hybrid evolutionary optimization algorithms can show the satisfactory performance and SWARA-ANFIS-GWO can achieve more accurate prediction accuracy than that of SWARA-ANFIS-WOA.

4. Discussion

These are no universal guidelines for choosing landslide influencing factors, the number and range of the classes of these factors [51]. Therefore, based on the geomorphological and geological characteristics of the study area, data availability and other landslide susceptibility studies, fifteen influencing factors were taken into account for modelling the process in this work. Meanwhile, the SWARA technique was used to assess the correlation between landslides and influencing factors. Accordingly, the SWARA weights revealed that most landslides occurred with a slope above 47.6°. In this case, the shear stress is dominated on the shear strength which results in landslide incidence using the gravity force. Also, the south-facing slope is more susceptible to landslide occurrence in comparison to the other directions. Since most of the rainfall in the study area occurred in this aspect, wet-and-dry cycles lead to decreasing the shear strength of soil and creating the landslide. This result is in concordance with the work in [51], which states that the rainfall is a crucial factor for landslide occurrence on south-facing slopes. When the altitude is lower than 432 m or higher than 1032 m, the frequency of landslide occurrence is very low. Although most of the landslides (148 landslides) occurred in the class of 432–582 m which covered an area of about 22.87%, the highest SWARA weight was obtained between 882 and 1032 m. However, at the middle altitude, 14 landslides occurred only have an area of 1.46% of the study area. The plan and profile curvature are less than -0.001, which represent that concave (negative) from of slope is the most susceptible to landslide occurrence due to having more concentration of water and erosion rather than convex (positive value) and flat (zero value). The SPI can indicate the erosion power of the stream, and the SWARA weights are larger as the SPI is increased. However, few landsides (only 3 landslides) occurred in the last class of this factor (800–1000). The STI classes of 30–40 and 50–60 obtained the highest SWARA weight, which showed the significant influence on the landslide occurrence. The SWARA weight of TWI increases when larger areas are involved in runoff generation. Since 200 landslide locations (63.49%) occurred in TWI between 3 and 5 (36.93% of the study area), this class of TWI was more effective in landslide occurrence. The most significant class of lithology was obtained for C unit (Kimberley rock) with the SWARA weight of 0.36. It was noted that only one landslide event occurred in this class of lithology, which covered the smallest section of the study area (less than 1%). The highest SWARA weight (0.21) was acquired for distance to fault less than 500 m, which indicated that about 19% of the landsides occurred due to fault activities in 13.63% of the study area. Results of correlation between distance to river and landside occurrence also depicted that the class of 600–750 m distance from the river was more effective on landslide incidence in the study area. The land use has indirect impact on landslide occurrence. Basically, forest class of land use was specified as the most susceptible class rather than the other classes due to having 243 landslides (77.14%) on 48.52% of the study area. The results also concluded that NDVI between 0.35 and 0.45, ACu type of soil, and rainfall between 1500 and 1600 mm were recognized as the most susceptible classes to landslide occurrence rather than other classes of influencing factors. In other words, these classes were more effective and had a significant role on landslide modelling process and preparing susceptibility maps.

Recently, hybrid machine learning model had been widely applied for landslide susceptibility modelling, for example, Pham et al. [136] used the novel hybrid model Bagging-based Naïve Bayes Trees (BAGNBT) at Mu Cang Chai district, located in northern Viet Nam, and the result showed that BAGNBT is a promising and better alternative method for landslide susceptibility modeling and mapping. Moayedi et al. [137] applied artificial neural network (ANN) optimized with particle swarm optimization (PSO) for the problem of landslide susceptibility mapping (LSM) prediction, and the result showed that PSO-ANN model showed higher reliability in estimating the LSM compared to the ANN. Tien Bui et al. [90] developed a new ensemble model which is a combination of a functional algorithm, stochastic gradient descent (SGD) and an AdaBoost (AB) Meta classifier namely ABSGD model to predict the landslides in the Sarkhoon watershed, located within the Zagros Mountains, Iran, and the result showed that the combined use of a functional algorithm and a Meta classifier prevents over-fitting, reduces noise, and enhances the power prediction of the individual SGD algorithm for the spatial prediction of landslides. To summarize, there is no algorithm that works perfectly for all optimization problems, and new algorithms must be applied and verified to determine the one that is most efficient [113]. It is also recommended to use some soft-computing benchmark algorithms and conventional methods to assess the validity of the proposed models. In order to manage and decrease the cost and destructive effects of landslide disasters, policy makers, governments, planners, and managers need to develop more reliable and precise landslide susceptibility maps.

5. Conclusions

Landslides are one of the most influential environmental challenges due to transiting a large amount of sediment to the closest stream, and loss of life and property as well. Therefore, landslide susceptible maps are the most practical tool to landslide mitigation. In this work, we introduce two novel hybrid machine leaning methods for LSM. Specifically, the ANFIS-based technique was improved by being combined with WOA, GWO, and PSO, and three proposed SWARA-ANFIS-WOA, SWARA-ANFIS-GWO and SWARA-ANFIS-PSO methods are introduced and the result is compared with SWARA and SWARA-ANFIS without any optimization process. For the proposed framework, the landslide inventory map of the study area is first compiled, and the influencing factors are selected. Then, the SWARA algorithm is used to obtain the initial weight of each class of landslide condition factors. Next, the landslide modelling is conducted on the study area using SWARA-ANFIS-WOA, SWARA-ANFIS-GWO, and SWARA-ANFIS-PSO. Finally, the landslide susceptibility maps are assessed using several objective measures. In general, the conclusion of this can be summarized as follows:

• (1). SWARA-ANFIS-GWO has the best result in terms of the convergence of objective function in comparison to SWARA-ANFIS-WOA.

• (2). SWARA-ANFIS-GWO outperforms SWARA-ANFIS-WOA and thus is a better tool for landslide susceptibility modelling.

• (3). The hybrid SWARA-ANFIS-GWO method has better goodness-of-fit and prediction accuracy and can produce more accurate reliable results. Therefore, it is recommended for LSM in other dangerous areas.

• (4). In the study area, the landslides are prone to occur under conditions of slop with >47.6°, south aspect, altitude with 882–1032 m, plan curvature with >0.001, profile curvature with <−0.001, SPI with 800–1000, TWI with 3–5, lithology of Lu jing Unit, distance to fault with 0–500 m, distance to river with 600–750 m, forest land-uses, NDVI with 0.35–0.45, soil of ACu and 1500–1600 mm rainfall.

Therefore, we recommend the ANFIS-GWO hybrid method as a promising technique to be used in other areas with more caution. It should be noted that the results of this study were obtained for the study area, which proposes that further studies should be considered in other areas with different geomorphological and geological characteristics. Since there are some uncertainties such as landslide inventory map, methods, raster resolution, and sample size, further attention should be paid to achieve reliable landslide susceptibility maps in regional scales.

Author Contributions

W.C., H.H., M.P., H.S., Y.W., A.S., S.P., A.A.A., K.K., S.P., F.R., S.L., A.J., D.T.B. and B.B.A. contributed equally to the work. W.C., H.H., Y.W. and S.L. collected field data and conducted the landslide susceptibility mapping and analysis. H.H., M.P., H.S., A.S., S.P., A.A.A., K.K., S.P., F.R. and W.C. wrote and revised the manuscript. S.L., A.J., D.T.B., B.B.A. and S.L. provided critical comments in planning this paper and edited the manuscript. All the authors discussed the results and edited the manuscript.

Funding

This research was supported by the National Natural Science Foundation of China (41807192, 61271408), Natural Science Basic Research Program of Shaanxi (Program No. 2019JLM-7, Program No. 2019JQ-094), project funded by Key Laboratory of Mine Geological Hazard Mechanism and Control, Shaanxi Provincial Natural Resources Department (Program No. KF2018-05), and Universiti Teknologi Malaysia (UTM) based on Research University Grant (Q.J130000.2527.17H84).

Conflicts of Interest

The authors declare no conflict of interest.

Figures and Tables

Figure 1. The study area and landslide locations.

Figure 2. The geological map of the Anyuan County in China.

Figure 3. Landslide influencing factor maps. (a) Altitude, (b) slope, (c) aspect, (d) plan curvature, (e) profile curvature, (f) SPI, (g) STI, (h) TWI, (i) lithology, (j) distance to fault, (k) distance to river, (l) land use, (m) NDVI, (n) soil and (o) rainfall.

Figure 3. Landslide influencing factor maps. (a) Altitude, (b) slope, (c) aspect, (d) plan curvature, (e) profile curvature, (f) SPI, (g) STI, (h) TWI, (i) lithology, (j) distance to fault, (k) distance to river, (l) land use, (m) NDVI, (n) soil and (o) rainfall.

Figure 4. Flowchart of the proposed framework.

Figure 5. Bubble-net feeding behavior of humpback whales [103].

Figure 6. Flowchart of the Whale Optimization Algorithm (WOA) algorithm.

Figure 7. The social hierarchy of grey wolves [104].

Figure 8. Flowchart of the grey wolf optimizer (GWO) algorithm.

Figure 9. An example scheme of modelling producer of evolutionary optimization algorithms.

Figure 10. The SWARA-ANFIS-WOA method. (a) Target and outputs for the training set; (b) target and outputs for the verification set; (c) MSE and RMSE for the training set; (d) frequency errors for the verification set; (e) MSE and RMSE for the verification set; (f) frequency errors for the verification set.

Figure 11. The SWARA-ANFIS-GWO method. (a) Target and outputs for the training set; (b) target and outputs for the verification set; (c) MSE and RMSE for the training set; (d) frequency errors for the verification set; (e) MSE and RMSE for the verification set; (f) frequency errors for the verification set.

Figure 12. The SWARA-ANFIS-PSO method. (a) Target and outputs for the training set; (b) target and outputs for the verification set; (c) MSE and RMSE for the training set; (d) frequency errors for the verification set; (e) MSE and RMSE for the verification set; (f) frequency errors for the verification set.

Figure 13. Convergence plot of the SWARA-ANFIS-GWO, SWARA-ANFIS-WOA and SWARA-ANFIS-PSO methods.

Figure 14. Landslide susceptibility maps. (a) SWARA-ANFIS-WOA, (b) SWARA-ANFIS-GWO, (c) SWARA, (d) SWARA-ANFIS, (e) SWARA-ANFIS-PSO.

Figure 14. Landslide susceptibility maps. (a) SWARA-ANFIS-WOA, (b) SWARA-ANFIS-GWO, (c) SWARA, (d) SWARA-ANFIS, (e) SWARA-ANFIS-PSO.

Figure 15. Percentages of different landslide susceptibility classes.

Figure 16. ROC curves for all methods using (a) the training dataset and (b) the verification dataset.

Table 1

Geological formations of the study area.

No.	Unit Name	Lithology
A	Hutian group	Dolomite, biogenic limestone, dolomitic limestone
		Two long (K-feldspar) granite
B	Zhong peng group	Red sandstone, siltstone, silty shale
	Zi shan group, Yun shan group	Quartz conglomerate and carbonaceous shale
C	Lu jing Unit	Kimberley rock
D	Lin shan group, Shui bei group	Carbonaceous shale and coal seam
	Luo ao group	Feldspar quartz coarse sandstone
	Ge xianshan group, Xishan pai unit	Granite
	Linsha group, Yuexing group	Granite
	Yu tian group, Ji long group	Residual tuff; crystalline chip
	Long zhou unit	Diorite
E	Shui jiang unit, Nankeng unit	Monzogranites
	Gan zhou group, Zhou tian group	Basalt, calcium mirabilite and halite
F	Yang jiaqiao group	Quartzite, siliceous rock
G	Da long group, Le ping group, Chang xing group	Coal seam, flint limestone, siliceous rock
	Tao xi group	Biotite plagioclase leptynite, sericite schist
H	Xia keng group	Granodiorite
I	Hu cheng unit	Tonalite
	Gui keng unit, Shang you unit	Granite
J	Le chang group	Slate, phyllite
	He zi unit	Peridotite, ophiolite, amphibolite
K	Bacun group, Niu jiaohe group	Silty slate, carbonaceous slate

Notes: A, B, C, D, E, F, G, H, I, J and K represent the lithology classes.

Table 2

Spatial relationship between landslides and influencing factors by the step-wise weight assessment ratio analysis (SWARA) algorithm.

Influencing Factors	Classes	No. of Pixels	No. of Landslide	SWARA Weight
Slope angle (degree)	<20.0	2656308	0	0.01
	20.0–24.6	460317	116	0.11
	24.6–29.2	312393	95	0.14
	29.2–33.8	188594	46	0.11
	33.8–38.4	101643	35	0.16
	38.4–43.0	46703	16	0.16
	43.0–47.6	18080	4	0.10
	>47.6	7394	3	0.19
Aspect	Flat	1200	0	0.00
	North	504881	12	0.03
	Northeast	433690	14	0.04
	East	479731	40	0.13
	Southeast	482847	70	0.23
	South	487601	97	0.30
	Southwest	431966	45	0.16
	West	480314	22	0.07
	Northwest	489202	15	0.04
Altitude (m)	132–282	739706	4	0.01
	282–432	1504606	75	0.06
	432–582	867013	148	0.22
	582–732	442338	39	0.11
	732–882	174022	35	0.26
	882–1032	55289	14	0.33
	>1032	8458	0	0.00
Plan curvature	<−0.001	1791191	146	0.43
	(−0.001)–(0.001)	171530	0	0.06
	>0.001	1828711	169	0.51
Profile curvature	<−0.001	1835414	184	0.55
	(−0.001)–(0.001)	33142	0	0.06
	>0.001	1922876	131	0.39
SPI	0–50	2957565	221	0.10
	50–100	372535	48	0.16
	100–200	224014	25	0.14
	200–400	129200	10	0.11
	400–800	66033	8	0.15
	800–1000	12254	3	0.32
	>1000	29831	0	0.01
STI	0–10	2366364	90	0.04
	10–20	797385	129	0.16
	20–30	300375	43	0.15
	30–40	129124	22	0.17
	40–50	67080	9	0.14
	50–60	39877	7	0.17
	>60	91227	15	0.16
TWI	<3	444	0	0.01
	3–5	1400515	200	0.61
	5–7	1708346	98	0.23
	7–9	514871	17	0.14
	>9	167256	0	0.01
Lithology	A	26836	0	0.01
	B	27889	4	0.12
	C	2422	1	0.36
	D	1159282	126	0.09
	E	710946	32	0.04
	F	76373	6	0.06
	G	508520	41	0.07
	H	6902	0	0.01
	I	142618	21	0.13
	J	719372	61	0.07
	K	410272	23	0.05
Distance to faults (m)	0–500	516832	60	0.21
	500–1000	466826	37	0.15
	1000–1500	436370	44	0.19
	1500–2000	381470	43	0.20
	2000–2500	329017	26	0.15
	>2500	1660917	105	0.11
Distance to rivers (m)	0–150	1001961	47	0.07
	150–300	799934	72	0.15
	300–450	703245	64	0.15
	450–600	543309	59	0.18
	600–750	388246	48	0.21
	750–900	220837	16	0.12
	>900	220837	16	0.11
Land use	Bare	656	0	0.00
	Forest	1839583	243	0.53
	Grass	748647	45	0.24
	Residential	189281	5	0.11
	Farmland	994815	22	0.09
	Water	18450	0	0.00
NDVI	<(−0.05)	110766	2	0.04
	(−0.05)–0.05	195838	7	0.07
	0.05–0.15	503739	18	0.07
	0.15–0.25	1185193	58	0.10
	0.25–0.35	1322962	141	0.20
	0.35–0.45	455283	88	0.40
	>0.45	17651	1	0.12
Soil	ATc	319882	2	0.02
	ACu	1035811	153	0.54
	ALh	117324	6	0.19
	ACh	2290342	154	0.25
	RGc	28073	0	0.00
Rainfall (mm)	<1100	85055	0	0.01
	1100–1200	427735	16	0.05
	1200–1300	1631694	74	0.06
	1300–1400	1040983	134	0.17
	1400–1500	433154	67	0.20
	1500–1600	78140	13	0.21
	1600–1700	38958	5	0.17
	>1700	55713	6	0.14

Table 3

Calculate weight of influencing factors by the SWARA algorithm.

Influencing Factors	CIAV	Aj	wj = (W(j−1))/Aj	Weight wj/Sigma wj
Slope angle (o)		1.000	1.000	0.105
Rainfall	0.05	1.050	0.952	0.100
Lithology	0.1	1.100	0.866	0.091
Soil	0.06	1.060	0.817	0.086
Land use	0.06	1.060	0.771	0.081
NDVI	0.07	1.070	0.720	0.076
Distance to Rivers (m)	0.07	1.070	0.673	0.071
Distance to faults (m)	0.03	1.030	0.653	0.069
TWI	0.05	1.050	0.622	0.066
SPI	0.08	1.080	0.576	0.061
Altitude	0.1	1.100	0.524	0.055
STI	0.1	1.100	0.476	0.050
Aspect	0.3	1.300	0.366	0.039
Profile curvature	0.3	1.300	0.282	0.030
Plan curvature	0.5	1.500	0.188	0.020