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

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 Appl. Sci. 2019, 9, 3755; doi:10.3390/app9183755 www.mdpi.com/journal/applsci Appl. Sci. 2019, 9, 3755 2 of 32 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.


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- [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,. 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.

Introduction of the Study Area
Anyuan County, which has an area of about 2374.59 km 2 , 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. 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.   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.

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 2 and 7,840 m 2 , 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 2 ) landslides, which have influenced 2352 people and cause the economic loss of 0.1 million US dollars, followed by medium-scale (200-1000 m 2 ) and large-scale (>1000 m 2 ) 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.

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     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 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.

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]: where is the number of experts and denotes the offered ranks by the experts for each factor. Then, the coefficient is defined as follows: where j indicate the number of the factors. The re-calculated weight is computed as follows:

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]: 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: where j indicate the number of the factors. The re-calculated weight D j is computed as follows: Finally, the relative weights of the criteria are represented as follows: where W j denotes the relative weight of the jth criterion and m represents the total of the criteria.

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 1 and x 2 are A 1 and B 1 , respectively, then f 1 = p 1 x 1 + q 1 x 2 + r 1 . Rule 2: If x 2 and x 2 are A 2 and B 2 , respectively, then f 2 = p 2 x 2 + q 2 x 2 + r 2 .
where A 1 , A 2 , B 1 , and B 2 , are the fuzzy sets, p i , q i and r i are the parameters, x 1 and 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].

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): where → X * and → X are two position vectors and → X * indicates the optimal solution obtained so far, t is the current iteration. Furthermore, → A and → C specify the coefficient vectors that are given by: 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 ⃗ and ⃗ decrease in Equations (7) and (8). The vector ⃗ 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 ( , ) and ( * , * ) indicate the positions of the whale and prey and ⃗ the distance between them; the whales mimic the helix-shaped movement according to spiral equation as follows: where ⃗ = ⃗ * ( ) − ⃗ ( ) is constant, and 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: 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 ⃗ has a random value larger than 1 or less than −1. The search for prey can where 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: 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 → A has a random value larger than 1 or less than −1. The search for prey can be modelled by the following equations: where the position vector → X rand is chosen randomly from the whales between the current population. The flowchart of the WOA algorithm is illustrated in Figure 6.

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

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].

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: where t represents the current iteration, and → X P and → X denote the position vectors of the prey and the grey wolves, respectively, and the coefficient vector is defined by → C as follows: where r 1 is a random vector which is in the interval [0, 1]. The prey's location can be represented as follows: The value of the coefficient → A is computed as follows: 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: where → X 1 , → X 2 and → X 3 indicate the position vectors of α, β and δ, respectively.
→ 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: The positions of other wolves update randomly according to the position of the prey.  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 | → 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 → C provides a random weight for the prey while | → C| > 1 and promotes the exploration phase. In addition, → C models the natural obstacles in hunting for the grey wolves [122]. The flowchart map of the GWO algorithm is shown in Figure 8.

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.
. . , v t in 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: where, x t i is the previous location of i-th particle, p t i is the best location found by g t i 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 ω, c 1 and C 2 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: 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. 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).

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.
Appl. Sci. 2019, 9, x 18 of 34 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.

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

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).

Validation of the Resultant Maps
To objectively evaluate the proposed methods, the ROC and AUC measures were used. ROC curve, the specificity and sensitivity are on the x-and y-axis, respectively [134,135].

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   Percentage of susceptibility class

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

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

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.

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.