The Application of the Hybrid GIS Spatial Multi-Criteria Decision Analysis Best – Worst Methodology for Landslide Susceptibility Mapping

The main goal of this article is to produce a landslide susceptibility map by using the hybrid Geographical Information System (GIS) spatial multi-criteria decision analysis best–worst methodology (MCDA-BWM) in the western part of the Republic of Serbia. Initially, a landslide inventory map was prepared using the National Landslide Database, aerial photographs, and also by carrying out field surveys. A total of 1082 landslide locations were detected. This methodology considers the fifteen conditioning factors that are relevant to landslide susceptibility mapping: the elevation, slope, aspect, distance to the road network, distance to the river, distance to faults, lithology, the Normalized Difference Vegetation Index (NDVI), the Topographic Wetness Index (TWI), the Stream Power Index (SPI), the Sediment Transport Index (STI), annual rainfall, the distance to urban areas, and the land use/cover. The expert evaluation takes into account the nature and severity of the observed criteria, and it was tested by using two scenarios: the different aggregation methods of the BWM. The prediction performances of the generated maps were checked by the receiver operating characteristics (ROCs). The validation results confirmed that the areas under the ROC curve for the weighted linear combination (WLC) and the ordered weighted averaging (OWA) aggregation methods of the MCDA-BWM have a very high accuracy. The results of the landslide susceptibility assessment obtained by applying the proposed best–worst method were the first step in the development of landslide risk management and they are expected to be used by local governments for effective management planning purposes.


Introduction
Landslides are the most dangerous natural geologic processes that cause different types of damage and also affect people, industries, and the environment, especially in times of dramatic climate change effects on the one hand, and urban sprawl and land consumption on the other.Throughout the world, hundreds of thousands of houses and buildings have been destroyed and many people have been injured by and died due to landslides.Because of this, landslide susceptibility mapping provides us with one of the crucial pieces of information in spatial planning for landslide susceptibility areas and is very important for the economic, cultural, environmental, and social sustainability of human beings.As a consequence, the study of landslide susceptibility mapping is rapidly becoming the focus of major scientific research, engineering study, and practices throughout the world [1].Many researchers have pursued work with the intention of predicting and preventing landslide hazards by using a wide variety of methods.A landslide susceptibility map represents the areas with the potential for landslides in the future by combining some of the critical factors that contributed to the occurrence of past landslides [2].
The spatial prediction of landslides using the production of landslide susceptibility maps is the first important step for landslide disaster and hazard management.The spatial probability of a landslide vulnerability can be expressed as the probability of the spatial occurrence of slope failures with a set of geo-environmental conditions.However, due to the complex nature of landslides, producing a reliable spatial prediction of landslide susceptibility is not easy [3].
Landslides are caused by various factors, such as topography, geology, soil characteristics, forest conditions, and climate variables [4].By applying various Geographical Information System (GIS)-based techniques and spatial data, different methods can be used to ultimately determine the level of landslide hazard and susceptibility [5].A GIS model can be used to combine a set of maps or factors using a function in order to produce the final map [6].The function can take numerous structures, including discriminate analysis, conditional analysis, linear regression, and multiple regression, etc. [7,8].There are several different approaches to the production of the landslide susceptibility and hazard map of a region (Table 1).
The main difference between the present study and the approaches described in the aforementioned references in Table 1 is that a GIS-based best-worst multi-criteria decision analysis (GIS-MCDA-BW) by different parameters was applied, its results being compared for landslide susceptibility in the western part of the Republic of Serbia.Specifically, we used a novel BW method to determine the connections in the network structure based on the criteria and accepted imprecisions during collective decision-making.Therefore, the quality of the existing data during collective decision-making and the experts' perceptions expressed through an aggregation matrix can be retained.The final susceptibility map was obtained by using the weighted linear combination (WLC) and the ordered weighted averaging (OWA) aggregation methods.

Deterministic methods
Deterministic methods are mainly based on the geotechnical and groundwater properties of the rock and soil of unstable areas.In this case, specific mathematical models are used to find the factor of the safety of unstable slopes, and slope stability models are used to determine landslide hazard [9].

Safety factor method
The safety factor method uses slope-displacement-simulated models, which are based on identifying the most dangerous sliding surface in order to calculate the factor for analyzing the slope stability [10][11][12].

Probabilistic methods
The probabilistic approach considers whether future environmental conditions will meet the requirements for a landslide identified in previous landslides.Thus, the probabilistic analysis considers the statistical relationships between the historical landslide locations and the landslide conditioning factors [13][14][15][16].

Analytic hierarchy process (AHP)
The AHP mainly depends on the knowledge of experts, who assign a priority to each parameter and establish sub-criteria from pairwise comparisons.The process is based on the three principles: decomposition, comparative judgment, and the synthesis of data [16][17][18].

Weighted linear combination (WLC) method
The WLC method starts with a comparison of the data-layers corresponding to the landslide controlling parameters and the landslide inventory map, and involves the computation of the landslide density so as to assign primary-level weights for each class of a particular parameter.The final steps of this method are a combination of all the weighted layers into a single map, and the classification of the scores of this map into landslide susceptibility [19,20].

Spatial multicriteria decision analysis (MCDA)
The MCDA can be defined as a decision aid and a mathematical tool allowing the comparison of different alternatives or scenarios according to many, often conflicting, criteria in order to guide the decision-maker towards a judicious choice [16,21].

Index of entropy (IoE) method
The entropy method has been widely used to determine the weighted index of natural hazards and carry out integrated environmental assessments of natural processes.The entropy of a landslide represents the factors influencing its development, and its value can be used to calculate the objective weights of the index system [22,23].

Bivariate analysis
In bivariate analysis, each individual factor is combined with a landslide distribution map, and the weight values based on the landslide densities are calculated for each parameter class [9,[24][25][26].

Multivariate analysis
In multivariate analysis, many causative factors are sampled, and for each of the sampling units, the presence or absence of landslides is also determined.This analysis allows the estimation of the relative weights of each contributing factor by means of statistical procedures.There is a trend towards using multivariate statistical analysis, such as discriminant analysis, factor analysis, logistic regression analysis, and conditional analysis [2,27].

Statistical index (SI) method
The SI method, bivariate statistical analysis, is considered as the simplest and quantitatively suitable method in landslide susceptibility mapping.In this method, the weighting value for each categorical unit is defined as the natural logarithm of the landslide density in the class divided by the landslide density in the whole studied area [28,29].

Weights of evidence (WoE) method
Weights of evidence (WoE), based on Bayesian Bayes' theorem and the assessment of the relationship between the spatial distribution of the areas affected by landslides and the spatial distribution of the conditioning factors causing landslides, is one of bivariate models [30][31][32].

Fuzzy logic method
The idea of using a fuzzy approach in landslide susceptibility mapping is to consider the pixels on any causal factor layer as susceptible to landslides.Pixel values can be numeric and range from 0 (i.e., "not susceptible") to 1 (i.e., "susceptible") [33][34][35].

Neuro-fuzzy method
The neuro-fuzzy method is the neural network that is functionally equivalent to the fuzzy inference model.It can be trained to develop "if-then" fuzzy rules and determine the membership functions for the input and the output variables of the system.One of the neuro-fuzzy inference systems is the adaptive neuro-fuzzy inference system (ANFIS) [41][42][43].

Support vector machines method (SVM)
The SVM method is a training algorithm based on the non-linear transformations that use the classification based on the principle of structural risk minimization, which has performed well in the test phase.The SVM model performs this process according to the three main concepts: the margin, the support vector, and the kernel [39,44,45].

Evidential belief function method
The EBFs are the compound of Bel (the degree of belief), Dis (the degree of disbelief), Unc (the degree of uncertainty), and Pls (the degree of plausibility).
The main parts of the theory are represented by Bel = a lower probability and Pls = an upper probability [45,46].

Decision tree method
The decision tree (DT) approach is a recently developed probabilistic approach based on the multivariate methods that are mainly used for classification schemes [47][48][49].

Naïve Bayes (NB) method
The NB classifier is a classification system based on Bayes' theorem assuming that all attributes are fully independent given the output class, called the conditional independence assumption [45,50].

Frequency ratio method
The frequency ratio is the ratio of the area where landslides have occurred and the total study area, also being the ratio of the probabilities of the landslide occurrence to non-occurrence for a given attribute [26,51,52].

Random forest method
Random forest is an ensemble-learning technique.It generates many classification trees aggregated so as to compute a classification.The Random forest algorithm has resistance to outliers in predictors and automatically handles the missing values.The random forest technique estimates the importance of a predictive variable [53].

Study Area
The study area covers approximately 5067 km 2 between the latitudes 43 • 46 26" to 44 • 56 10" N, and the longitudes 19 • 06 46" to 20 • 14 56" E, and is located in the western part of the Republic of Serbia (Figure 1).The study area consists of Serbian Mačva and Kolubara Districts, and the Tara National Park.The elevation of the study area ranges from 0 to 1586 m above mean sea level (MSL).

Study Area
The study area covers approximately 5067 km 2 between the latitudes 43°46′26″ to 44°56′10″ N, and the longitudes 19°06′46″ to 20°14′56″ E, and is located in the western part of the Republic of Serbia (Figure 1).The study area consists of Serbian Mačva and Kolubara Districts, and the Tara National Park.The elevation of the study area ranges from 0 to 1586 m above mean sea level (MSL).The slope angles of the area range from 0° to as much as 89°.The total annual rainfall in the different parts of the study region ranges from 623 to 1055 mm.Based on the records obtained from the Serbian Meteorological Service, the maximum rainfall occurs between March and June.The study area is heterogeneous in terms of the terrain complexity, which consists of plain and mountainous areas.The slope angles of the area range from 0 • to as much as 89 • .The total annual rainfall in the different parts of the study region ranges from 623 to 1055 mm.Based on the records obtained from the Serbian Meteorological Service, the maximum rainfall occurs between March and June.The study area is heterogeneous in terms of the terrain complexity, which consists of plain and mountainous areas.

Data
The detection of the accurate location of a landslide and the creation of the historical database record is a very important primary process for performing landslide susceptibility analysis.However, the detection of the landslide location is often a challenging and time-consuming process [54].The National Landslide Database was created by the Ministry of Mining and Energy of the Republic of Serbia within the BEWARE (Beyond Landslide Awareness) project, which represents a standardized post-event landslide database closely involving the local community of 27 municipalities, preparing them to cope with landslide hazard events in the future.The BEWARE (GIS) web portal is a platform for reporting interactive landslide events and unifying landslide data records.
The landslide inventory map is an important step in landslide susceptibility mapping.Therefore, field survey, high-resolution images, and aerial photo interpretation methods were applied in this study to generate the landslide inventory map.The landslide conditioning factor is yet another key topic, which has been researched many scientists.Hence, the elevation, the slope, the aspect, the topographic wetness index (TWI), the sediment transport index (STI), the stream power index (SPI), the normalized difference vegetation index (NDVI), the distance to the river, the distance to roads, the distance to the urban areas, rainfall, the land cover/use, lithology, and the soil type were used in this study in order to analyze landslide susceptibility.Detailed information regarding the landslide conditioning factors' sources are given in Table 2.

Methodological Background
The methodological hierarchy in this paper is based on the GIS multi-criteria decision analysis (GIS-MCDA) structure.This approach uses the capabilities of GIS in the management of geospatial data and the flexibility of the MCDA to combine factual pieces of information (e.g., land use, slope, aspect, TWI, etc.) with value-based information (e.g., expert opinion, standards, surveys, etc.) [55].The main advantage of integrating GIS and the MCDA can be seen in their specific capabilities supplementing each other [56].From a methodological point of view, the defined landslide susceptibility areas in the western part of the Republic of Serbia include the following main steps (Figure 2).

Methodological Background
The methodological hierarchy in this paper is based on the GIS multi-criteria decision analysis (GIS-MCDA) structure.This approach uses the capabilities of GIS in the management of geospatial data and the flexibility of the MCDA to combine factual pieces of information (e.g., land use, slope, aspect, TWI, etc.) with value-based information (e.g., expert opinion, standards, surveys, etc.) [55].The main advantage of integrating GIS and the MCDA can be seen in their specific capabilities supplementing each other [56].From a methodological point of view, the defined landslide susceptibility areas in the western part of the Republic of Serbia include the following main steps (Figure 2).The primary advantages of the BW method for the landslide susceptibility mapping [57][58][59] suggested by the authors are as follows: (1) compared with the analytic hierarchy process (AHP) method, which is most commonly used in the literature to determine weight coefficients [60], it requires significantly fewer pairwise comparisons (the AHP method requires n(n−1)/2 comparisons, the BW 2n−3 comparisons); ( 2) the values of the weight coefficients obtained by the BW method are more reliable because comparison in the BW method is carried out with a higher consistency ratio compared with the AHP method; (3) while for the majority of MCDM models (e.g., the AHP) the consistency ratio is a test of whether the criteria comparison is consistent or not, in the BW method the consistency ratio is used to determine the level of confidence since the BW outputs are always consistent; (4) the BW model only uses integers for the criteria pairwise comparison, as opposed to the other MCDM methods (e.g., the AHP), which also require the use of fractional numbers [61].
Since the method is of a very recent date, the literature has, so far, only contained the traditional (crisp) BW method [57][58][59]62] and the modification of the BW method performed by using fuzzy numbers [63].The following section presents the algorithm for the BW method that includes the following steps: The primary advantages of the BW method for the landslide susceptibility mapping [57][58][59] suggested by the authors are as follows: (1) compared with the analytic hierarchy process (AHP) method, which is most commonly used in the literature to determine weight coefficients [60], it requires significantly fewer pairwise comparisons (the AHP method requires n(n−1)/2 comparisons, the BW 2n−3 comparisons); ( 2) the values of the weight coefficients obtained by the BW method are more reliable because comparison in the BW method is carried out with a higher consistency ratio compared with the AHP method; (3) while for the majority of MCDM models (e.g., the AHP) the consistency ratio is a test of whether the criteria comparison is consistent or not, in the BW method the consistency ratio is used to determine the level of confidence since the BW outputs are always consistent; (4) the BW model only uses integers for the criteria pairwise comparison, as opposed to the other MCDM methods (e.g., the AHP), which also require the use of fractional numbers [61].
Since the method is of a very recent date, the literature has, so far, only contained the traditional (crisp) BW method [57][58][59]62] and the modification of the BW method performed by using fuzzy numbers [63].The following section presents the algorithm for the BW method that includes the following steps: Step 1. Determining a set of evaluation criteria.This starts from the assumption that the process of decision-making involves m experts.In this step, the experts consider a set of evaluation criteria and select the final set of the criteria C = {c 1 , c 2 , . . .c n }, where n represents the total number of the criteria.
Step 2. Determining the most significant (the best) and the least significant (the worst) criteria.The experts decide on the best and the worst criteria from the set of criteria C = {c 1 , c 2 , . . .c n }.If the experts decide on two or more criteria as the best, or the worst, the best and the worst criteria are arbitrarily selected.
Step 3. Determining the preferences of the most significant (the most influential) criteria (B) from within the set C over the remaining criteria from the defined set.Under the assumption that there are m experts and n criteria under consideration, each expert should determine the degree of the influence of the best criterion B on the criteria j (j = 1, 2, . . ., n).This is how a comparison between the best criterion and the other criteria is obtained.The preference of the criterion B compared to the j-th criterion defined by the e-th expert is denoted by a e Bj (j = 1, 2, . . ., n; 1 ≤ e ≤ m).The value of each pair a e Bj takes a value between 1 and 9 [58].
As a result, the Best-to-Others (BO) vector is obtained: where a e Bj represents the influence (preference) of the best criterion B on (over) criterion j, whereby a e BB = 1.This is how the fuzzy BO matrices A 1 B , A 2 B , . . ., A m B are obtained for each expert.
Step 4. Determining the preferences of the criteria from the set C over the worst criterion (W) from the defined set.Each expert should determine the degree of the influence of the criterion j (j = 1, 2, . . ., n) in relation to the criterion W. The preference of the criterion j in relation to the criterion W defined by the e-th expert is denoted as a e jW (j = 1, 2, . . ., n; 1 ≤ e ≤ m).The value of each pair a e jW takes a value between 1 and 9.As a result, the Others-to-Worst (OW) vector is obtained: where a e jW represents the influence (preference) of the criterion j in relation to the criterion W, whereby a e WW = 1.This is how the fuzzy OW matrices A 1 W , A 2 W , . . ., A m W are obtained for each expert., the matrices of the aggregated sequences of experts are formed as follows: where a e Bj = a 1 Bj , a 2 Bj , . . ., a m Bn represent the sequences of a e Bj , by means of which the relative significance of the criterion B is described in relation to the criterion j.
Finally, by applying Equation (3), the average BO matrix is obtained: where e represents the e-th expert (e = 1, 2, . . ., m), a e Bj represents the sequence that describes the relative significance of the criterion B in relation to the criterion j.Thus, the averaged BO matrix of the average responses A B is thus obtained: Step 6. Determining the OW matrix of the experts' average responses.Based on the OW matrices of the experts' responses A e W = a e jW 1×n , as with the BO matrices, the matrices of the experts aggregated sequences are formed: where a e jW = a 1 jW , a 2 jW , . . ., a m nW represents the sequences of a e jW , by means of which the relative significance of the criterion j is described in relation to the criterion W.
Finally, by applying Equation ( 6), the average OW matrix is obtained: where e represents the e-th expert (e = 1, 2, . . ., m), a e jW represents the sequence describing the relative significance of the criterion j in relation to the criterion W. Thus, the averaged OW matrix of the average responses A W is obtained: Step 7. The calculation of the optimal values of the weight coefficients of the criteria [w 1 , w 2 , . . . ,w n ] from the set C. The goal is to determine the optimal value of the evaluation criteria, which should satisfy the condition that the difference in the maximum absolute values ( 8) for each value of j is minimized.In order to meet these conditions, the solution that satisfies the maximum differences according to the absolute value w B w j − a Bj and w j w W − a jW should be minimized for all the values of j.The previously defined limits will be presented in the following min-max model: where w B represents the weight of the best criterion, w W the weight of the worst criterion, and w j the weight of the j-th criterion.Model (10) is equivalent to the following model: where w j represents the optimum values of the weight coefficients, w B and w W represent the weight coefficients of the best and the worst criteria, respectively, whereas a jW and a Bj , respectively represent the values from the average OW and BO vectors (see Equations ( 4) and ( 7)).
By solving Model (11), the optimal values of the weight coefficients of the evaluation criteria [w 1 , w 2 , . . ., w n ] and ξ * are obtained.

The Consistency Ratio of the Best-Worst Methodology (BWM)
The consistency ratio is a very important indicator, by means of which the consistency of the pairwise comparison of the criteria in the BO and OW matrices is checked.Definition 1.The comparison of criteria is consistent when the condition a Bj × a jW = a BW is fulfilled for all of the criteria j, where a Bj , a jW and a BW , respectively, represent the preference of the best criterion over the criterion j, the preference of the criterion j over the worst criterion, and the preference of the best criterion over the worst criterion.
However, when comparing criteria, some pairs of the criteria j may not prove to be completely consistent.Therefore, the next section defines the consistency ratio (CR), which provides us with information on the consistency of the comparison between the BO and the OW vectors.In order to show how the CR is determined, we start from the calculation of the minimum consistency when comparing criteria, which is explained in the following section.
As previously indicated, the pairwise comparison of the criteria is carried out based on the predefined scale in which the highest value is 9, or any other maximum from the scale defined by the decision-maker.The consistency of the comparison decreases when a Bj × a jW is less or greater than a BW , i.e., when a Bj × a jW = a BW .It is clear that the greatest inequality occurs when a Bj and a jW have the maximum values equaling a BW , which continues to affect the value of ξ.Based on these relationships, the following conclusion can be drawn: As the largest inequality occurs when a Bj and a jW have their maximum values, the value ξ needs to be subtracted from a Bj and a jW , and a BW needs to be added.Thus, Equation ( 13) is obtained: Given the fact that the minimum consistency a Bj = a jW = a BW applies, Equation ( 14) is presented as follows: By solving Equation (14) for the different values of a BW , the maximum possible values of ξ can be determined, which is the CI for the BWM.Since the values of a BW are obtained on the basis of the experts' aggregated decisions, and given the fact that these change the interval, the values of ξ are impossible to predefine.The values of ξ depend on uncertainties in decisions since such uncertainties change the average value.As explained in the algorithm for the BWM average values of a Bj , a jW and a BW change depending on uncertainties in evaluating the criteria.
Based on CI, the consistency ratio (CR) is obtained: The CR takes values from interval [0, 1], where values closer to zero show a high consistency, whereas the values of the CR closer to one show a low consistency.

Results
This section is divided into the subsections.It provides a concise and precise description of the aim of this study.It is aimed at proposing a reliable GIS-MCDA BW methodology for the landslide susceptibility mapping, which could serve as a useful tool for preventing and reducing the landslide hazard for spatial planners to create spatial policies and systems for landslide management.
The various thematic data layers representing the landslide conditioning factors, such as the slope, the aspect, the elevation, lithology, the land use/cover, the distance to faults, the distance to rivers, the distance to roads, the topographic wetness index (TWI), the stream power index (SPI), the sediment transport index (STI), rainfall, the distance to urban areas, the soil type, and the normalized difference vegetation index (NDVI), were prepared.These factors fall under the three categories of the conditioning factors that make the area susceptible to movement without actually initiating a landslide; thus, these factors are considered to be responsible for the occurrence of landslides in the regions for which pertinent data can be collected from the available resources and from the field.The triggering factor, such as rainfall, sets the movement off by shifting the slope from a marginally stable to an actively unstable area [1].Furthermore, the attributes of the ground in terms of landslide susceptibility were considered in the present study.

Landslide Inventory Map
The mapping of the actual landslides in the study area is essential for describing the relationship between the landslide distribution and the conditioning factors.In order to produce a detailed and reliable landslide inventory map, extensive field surveys and observations were performed in the study area.A total of 1082 landslides were identified and mapped in the study area by evaluating the aerial photos supported by the field survey.Also, we used the National Landslide Database of the Republic of Serbia within the BEWARE project, which represents a standardized post-event landslide database.The landslide standard spatial database contains data such as the exact location of the landslide, the width, length, depth, and volume of the landslide, the average slope data of the terrain on which the landslide is located, the general information on the landslide process, the type of the occurrence, the trend of the movement, the type of the triggered material, water content, the movement speed, the activity, the mode of the movement, as well as the terrain, vulnerability, and the landslide photos.

Conditioning Factors
The landslide hazard evaluation criteria selection is an important step in the analysis.It is essential to identify the landslide conditioning factors in order to create a reliable landslide susceptibility map.Based on the experts' opinions and longer observations from the field, this study adopted the fifteen criteria that are an important cause of the landslide susceptibility assessment.
The selected criteria, together with a brief description of the same, are given in Table 3, and they are shown in Figures 3-5.
The lithology classes are defined by the different types of the geological formation and they are accounted for in Table 4.

Category
Factor Description Topography (Figure 3) The elevation is a significant landslide conditioning factor because it controls several geologic and geomorphologic processes [64].An elevation map is prepared from the 20 × 20 m digital elevation model (1: 25,000 scale with 10-m contour intervals) and grouped into 6 classes.

Topo2
The slope is widely used in landslide susceptibility studies since it is directly connected with the movement of landslide materials [49].Specifically, shear stresses on the slope material increases with the slope gradient and landslides are generally expected to occur on the steepest slopes.

Topo3
The aspect affects the slope material in an indirect relationship because the aspect determines the exposure of a landscape to rainfall and solar radiation, and therefore, to the propensity of vegetation to grow, which in turn affects the soil stability.

Topo4
The topographic wetness index (TWI) describes the effect of topography on the location and size of the saturated areas of the runoff generation.It is defined as [65]: TWI = ln (A S /tan β), where A S is the catchment area and β is the slope angle measured in degrees.

Topo5
The stream power index (SPI) is the measure of the erosive power of flowing water based on the assumption that discharge is proportional to the specific catchment area.The stream power index was calculated based on the formula given by Moore [66].SPI = A S × tanβ, where A S is the area of the specific catchment and β is the local slope gradient measured in degrees.

Topo6
The sediment transport index (STI) describes the tendency of the flow and can be used to depict the location of a potential erosion.It is calculated by using the following formula: STI = (A S /22.13) 0.6 × (sinβ/0.0896) 1.3 , where A S is the area of the specific catchment and β is the local slope gradient measured in degrees.

Env1
The soil type reflects the textures and compositions of the soil materials affecting the landslide occurrence [67].The soil map was constructed from the Basic Engineering National Soil Map at the scale 1:000,000, and was classified into fine-silt, course-loamy, fine-loamy, mixed-loamy, skeletal-loamy.

Env2
The drainage system of any area plays an important role in the slope stability particularly with respect to toe cutting and the bank erosion.The distance to the river was created by using a topographical map and was calculated based on the Euclidean distance method in ArcGIS 10.4 and the obtained distances were classified into (<500), (500-1000), (1000-2000), (2000-3000), and (>3000) m classes.

Env3
Lithology.The underlying geology is part of the most significant factors for landslide modeling [68].Different geology formations have different compositions and structures which contribute to the strength of the material.The stronger rocks give more resistance to the driving forces as compared to the weaker rocks and, hence, are less prone to landslides.The lithology structure of the study area includes 18 classes.

Env4
Distance to fault.The distance from the faults is calculated at 100 m intervals by using the geological map.Faults are the tectonic breaks that usually decrease the rock strength.These dislocations may be responsible for triggering a large number of landslides.

Env5
The normalized difference vegetation index (NDVI).The NDVI map was produced from the Lands at 8 OLI imagery showing the surface vegetation coverage and density in an image.

Soc1
Landslides may occur on the road and on the side of the slopes affected by roads.The distance to roads was created by using a topographical map and calculated based on the Euclidean distance method in ArcGIS 10.4, and the obtained distances were classified into the (<500), (500-1000), (1000-2000), (2000-3000), and (>3000) m classes.

Soc3
Land use/cover is considered to be a factor in environmental protection.The data on the land use/cover were taken on the basis of the Corine Land Cover 2006 (CLC2006) database, collected within the framework of the European Commission's CORINE (Coordination of Information on the Environment) program.The land use also plays a significant role in the stability of the slope.The land covered with a forest regulates the continuous water flow and water infiltrates regularly, whereas the cultivated land affects the slope stability due to the saturation of the covered soil.The description of the soil types based on the codes is shown in Table 5.The description of the land cover/land use based on the codes is shown in Table 6.

GIS MCDA-BW Methodology
This phase involves the standardization, expert work, weighting, summary analysis, the aggregation of all of the criteria to be considered in the decision-making process, and their validation as well.
Given the fact that the data were collected in different ways and that they have different formats, the first step implies that all the datasets have to be standardized and expressed in the units that can be compared.Based on the literature and the experts' experiences, the fuzzy concept applied in this study was used to standardize the data criteria.The fuzzy logic concept is flexible and suitable for the data modeling in which there are no exact boundaries of the elements belonging to the set, determined as either zero or one [69].In such cases, the elements belonging to the set are defined on the basis of the degree of affiliation to a function (sigmoidal, J-shaped, linear or user-defined).Which membership functions will be used depends on the nature of the data and the experts' opinions.
In this case, with the criteria set whose elements have categorical values (the land cover use), the discrete classification in which the experts directly determined the values of the elements of the fuzzy sets was applied.With respect to the other criteria, which are the values of the gradual change from one location to another, the elements of the set were standardized by using the fuzzy concept based on the linear membership function.The scale ranging from 0 to 1 byte was used for the purpose of fuzzification, where zero stands for the least susceptibility, and one represents the most dangerous element of the set value in relation to the likelihood of the occurrence of the landslide (Table 7).After standardization, it is necessary for decision-makers to define the significant factors of particular criteria by using the appropriate coefficient weights (weights) or the criteria weights.The BWM was used for the analysis of the factors carried out by the experts.This research included six experts.Having defined the evaluation clusters/criteria within the framework of the clusters and each group of the criteria, the experts also determined the best (B) and the worst (W) clusters/criteria.On this basis, the experts determined the BO and the OW vectors, in which the preferences of the B and the W over the clusters/criteria were considered for the remaining clusters/criteria from the defined set.The evaluation of the clusters/criteria was carried out by using the scale [1,9]: 1-a very low influence; 2-a low influence; . . .; 8-a high influence; 9-a very high influence.The BO and the OW vectors are presented in Table 8.C32 (Distance to roads) 2; 2; 3; 2; 3; 2; 3 C31 (Land use/cover) 2; 3; 3; 3; 3; 4; 4 C33 (Distance to urban areas) 2; 3; 3; 3; 3; 3; 2 C32 (Distance to roads) 2; 3; 2; 3; 3; 2; 3 Using Equations ( 4) and ( 7) we obtain the average BO and OW vectors (Table 9).On the basis of the BO and the OW vectors for each group of the clusters/criteria (Table 8), the optimal values of the weight coefficients of the clusters/criteria were calculated.Table 9 displays the four BO or OW vectors, which means that the four models (12) were formed for calculating the optimal values of the weight coefficients of the clusters/criteria: Based on the above models, the optimal values of the weight coefficients were obtained as demonstrated in Table 10.Table 10 presents the global and the local values of the weight coefficients of the criteria.The global weights of the criteria were obtained by multiplying the weight coefficients of the clusters by the weight coefficients of the sub-criteria.The global weight criteria continue to be used for the evaluation of the alternatives in the multi-criteria model [70].
By solving Model (11), the values of ξ * are obtained, which are as follows: ξ * Clusters = 0.75792, ξ * C Topo = 0.02467, ξ * C Envir = 0.09899 and ξ * C Soc = 0.66228.The values of ξ * are used to determine the consistency ratio, Equation (16).Since the value of a BW was obtained on the basis of the experts' aggregated decisions, the values of the consistency index ξ are impossible to predefine.By applying Equation (15), the values of the consistency index (ξ) were defined, as shown in Table 11.As seen in Table 11, the CR values are satisfactory.

Aggregation by Applying a Weighted Linear Combination
The weighted linear combination (WLC) method is used in the process of the criteria map aggregation.In addition, it is compensatory, meaning that the low scores in one criterion can be compensated for by the high scores in another, which is desired for this particular decision problem [71].For these reasons, the WLC was selected as the aggregation method.The weighted linear combination (WLC) method multiplies each fuzzy standardized criteria map (i.e., each raster cell 20 × 20 m) by the criteria weights, thus obtaining different variations from the GIS-MCDA-BW method, and then sums the results.The mathematical expression for calculating the suitability index in the WLC is given as follows ( 17): where S is the suitability index, w i is the normalized value of the factor weight, and x i is the criterion score of the factor i.
Based on the adopted criteria and the determination of their weights, and according to the three scenarios described in the previous part, the WLC is used to execute the aggregation map of the criteria in the final landslide susceptibility map, which is presented in the same fuzzy value range from 0 to 1. Finally, the landslide susceptibility index is calculated by using the defuzzification algorithm by the standard deviation method from the reclassify Spatial Analyst Tools ArcGIS 10.4 software release [72].Based on this, each cell is classified into five categories and receives a new value from 1 to 5, representing the Landslide Susceptibility Index (LSI).The results of the landslide susceptibility assessment are given in the maps displayed in Figure 6.Value 1 is the area with the least probability of a landslide occurrence, whereas Value 5 represents the areas with the highest probability of a landslide occurrence.

Aggregation by Applying Ordered Weighted Averaging
Ordered weighted averaging (OWA) is a class of multicriteria operators, which was given quantifier-guided aggregation [73].In this method, the first set of weights controls the relative contribution of a specific criterion, whereas the second set of weights controls the order of the aggregation of the weighted criteria [74].
The OWA provides a tool for generating a wide range of decision strategies in a decision strategy analysis by applying a set of order weights to the criteria ranked in order on a location-by-location (object-by-object) basis.The order weights process is central to the OWA combination procedures.The number of order weights is equal to the number of the criteria and must sum to one.The position of the set of order weights can be identified in the multi-criteria decision analysis based on the tradeoff and risk concepts [75,76].
For a given set of n criterion (attributes) maps, an OWA operator can be defined as the following function [77]:

Aggregation by Applying Ordered Weighted Averaging
Ordered weighted averaging (OWA) is a class of multicriteria operators, which was given quantifier-guided aggregation [73].In this method, the first set of weights controls the relative contribution of a specific criterion, whereas the second set of weights controls the order of the aggregation of the weighted criteria [74].
The OWA provides a tool for generating a wide range of decision strategies in a decision strategy analysis by applying a set of order weights to the criteria ranked in order on a location-by-location (object-by-object) basis.The order weights process is central to the OWA combination procedures.The number of order weights is equal to the number of the criteria and must sum to one.The position of the set of order weights can be identified in the multi-criteria decision analysis based on the trade-off and risk concepts [75,76].
For a given set of n criterion (attributes) maps, an OWA operator can be defined as the following function [77]: where:

Validation
The validation of landslide susceptibility maps (LSMs) is an essential step in the modeling process.The capability of the WLC and the OWA methods was evaluated by applying a nondependent threshold approach: the receiver operating characteristic (ROC) curve.The area under the

Validation
The validation of landslide susceptibility maps (LSMs) is an essential step in the modeling process.The capability of the WLC and the OWA methods was evaluated by applying a non-dependent threshold approach: the receiver operating characteristic (ROC) curve.The area under the curve (AUC) is the synthetic index calculated for the ROC curves, and the same has generally been used in several types of research studies to evaluate the accuracy of the landslide susceptibility map [78].The AUC is the probability that a positive event is classified as positive by doing the test given all the possible values of the test.The ROC curves are generated by XLSTAT (statistical software for Microsoft Excel) and represent the evolution of the proportion of true positive cases (also called sensitivity) as a function of the proportion of false positive cases (corresponding to 1 minus the specificity).
Furthermore, the model of this article is based on expert knowledge and all of the historical landslide occurrence samples are taken into account.In this study, the historical landslide sites were randomly divided into the training sets (70% of the total) and the validation sets (30%).
The final results of the fifteen landslide susceptibility maps (the maps generated by the WLC and the OWA aggregation method of the best-worst model) were validated by being compared with the existing landslides locations, using the success rate and the prediction rate methods.The success-rate results were obtained based on the comparison of the landslide grid cells in the training dataset with the fifteen landslide susceptibility maps.The success rate measures how the results of the landslide susceptibility analysis fit the training dataset.This method divides the area of the landslide susceptibility map into five classes, ranging from the highest to the lowest LSI values (Figure 7).The success rate method uses the training dataset, so it may not be an appropriate method to evaluate the predictability of the landslide susceptibility models.
The prediction rate can explain how well the landslide susceptibility models and the landslide conditioning factors predict a landslide occurrence.In this study, the prediction rate results were obtained by comparing the landslide grid cells in the validation dataset with the fifteen landslide susceptibility maps.
The area under the prediction rate ROC curve (AUC) was calculated.Figure 8 shows the prediction rate results of the fifteen landslide susceptibility maps obtained from the WLC and the OWA aggregation methods for the GIS-MCDA-BW methodology.
ISPRS Int.J. Geo-Inf.2019, 8 FOR PEER REVIEW 8 The final results of the fifteen landslide susceptibility maps (the maps generated by the WLC and the OWA aggregation method of the best-worst model) were validated by being compared with the existing landslides locations, using the success rate and the prediction rate methods.The successrate results were obtained based on the comparison of the landslide grid cells in the training dataset with the fifteen landslide susceptibility maps.The success rate measures how the results of the landslide susceptibility analysis fit the training dataset.This method divides the area of the landslide susceptibility map into five classes, ranging from the highest to the lowest LSI values (Figure 7).The success rate method uses the training dataset, so it may not be an appropriate method to evaluate the predictability of the landslide susceptibility models.
The prediction rate can explain how well the landslide susceptibility models and the landslide conditioning factors predict a landslide occurrence.In this study, the prediction rate results were obtained by comparing the landslide grid cells in the validation dataset with the fifteen landslide susceptibility maps.
The area under the prediction rate ROC curve (AUC) was calculated.Figure 8 shows the prediction rate results of the fifteen landslide susceptibility maps obtained from the WLC and the OWA aggregation methods for the GIS-MCDA-BW methodology.According to validation results, all of the final landslide susceptibility maps were considered to have the most acceptable and representable appearances (AUC > 0.9).The OWA aggregation technique exposed the overall best cross-validated performance, followed by the weighted linear combination method.With its value AUC = 0.941, the OWA aggregation technique is confirmed to have a very high accuracy of validation (Figure 8b).Also, both the visual assessment and the quantitative validation-through the ROC curve-agreed with the WLC aggregation technique as an excellent performing model approach with an AUC value 0.905 (Figure 8a).According to validation results, all of the final landslide susceptibility maps were considered to have the most acceptable and representable appearances (AUC > 0.9).The OWA aggregation technique exposed the overall best cross-validated performance, followed by the weighted linear combination method.With its value AUC = 0.941, the OWA aggregation technique is confirmed to have a very high accuracy of validation (Figure 8b).Also, both the visual assessment and the quantitative validation-through the ROC curve-agreed with the WLC aggregation technique as an excellent performing model approach with an AUC value 0.905 (Figure 8a).

Discussion
The main goal of this work was to investigate the use a novel MCDA-GIS BW technique for the purpose of evaluating landslide susceptibility mapping, which is constructed by the WLC and the OWA aggregation methods for the GIS-MCDA-BW methodology in the case study of Western Serbia.
The considered conditioning factors did not have the same contribution to a landslide occurrence, for which reason decision-makers should define the importance of each conditioning factor by using appropriate weight coefficients.The WLC and the OWA methods require weight normalization.After determining the weight coefficients of the conditioning factors, the hybrid GIS MCDA best-worst method was used to calculate the normalized weight conditioning factors.Finally, the WLC and the OWA methods were subsequently used in the GIS environment in order to obtain the final landslide susceptibility map.
According to the determination of the importance of the differently used factors for landslide susceptibility mapping, the results of the current study have shown that the most important conditioning factor is the slope, which is only followed by annual rainfall, the elevation, and the TWI [49,[79][80][81].
As the slope increases, the strain on the soil or another unconsolidated material generally increases as well.Due to the generally lower shear stress associated with low gradients, mild slopes are expected to have a low frequency of landslides.However, steep natural slopes resulting from the excavation of rocks may not be susceptible to severe landslides.For the most part, landslides occurred in the eastern and south-eastern facing of the slope, which is in accordance with the research reported by Pourghasemi et al. [31].
In fact, landslide susceptibility can be considered as a comprehensive assessment of the basic conditions of a disaster and a comprehensive measurement of the existing landslide characteristics, without considering the dynamic predisposing factors, such as extreme climatic conditions and human engineering activities, and without involving the problem of when a landslide will occur.Therefore, it is different from a timely warning or sequence analyses, without taking the time, frequency, and intensity of local landslides into consideration.The additional objective of this paper was to find spatial trends in the relationship between landslide susceptibility mapping and annual rainfall in the study area.However, each landslide event has a different precipitation and duration.Therefore, the uniformity of extreme rainfall and its duration cannot be taken as a conditioning factor.Finally, extreme rainfall and its duration will be used as the key factor in landslide susceptibility mapping in a future study.
In addition to the foregoing, the TWI is one of the conditioning factors having a highly significant impact on landsliding.The TWI describes the effect of topography on the location and size of the saturated areas of the runoff generation.The TWI shows the diversity and complexity of the topographical landslide surface, highlighting the areas of preferential landslide drainage and relatively dry areas within its borders.
On the other hand, the distance to rivers and the distance to urban areas referred to in this study had the lowest importance for a landslide occurrence.This is partly in agreement with the findings reported in Hong et al. [82] and Akgun et al. [14], even though the distance from faults was reported as the main effective factor [48,[83][84][85][86][87].
According to the determination of the importance of the used methodology, compared with the AHP method as the most commonly used heuristic and index-based approach in the literature referred to for the purpose of determining the weight coefficients, the best-worst method requires significantly fewer pairwise comparisons.In addition, the values of the weight coefficients obtained by the BW method are more reliable because comparison in the BW method is carried out with a higher consistency ratio compared with the AHP method [88].Furthermore, the consistency ratio for the majority of the MCDM models (e.g., the AHP) is the test of whether the comparison of criteria is consistent or not; in the BW method, the consistency ratio is used to determine the level of confidence since the outputs from the BW are always consistent.Ultimately, the BW model only uses integers for the pairwise comparison of criteria as opposed to the other MCDM methods (e.g., the AHP), which also require the use of fractional numbers.

Conclusions
A landslide is considered to be one of the most disastrous hazards throughout the world.Every year, it results in the death of many people and large economic loss.In order to reduce the landslide-resulting loss, many scholars have joined forces in conducting research studies on the cause and prediction of landslides.The landslide susceptibility map is the first step in the development of landslide risk management and local governments are expected to use it for effective management planning purposes.Therefore, the GIS-MCDA-BW technique was used in this study in the section dealing with the determination of the weights, namely with the significance of each criterion, and the WLC and the OWA techniques were used to sum the weights and finally to identify the landslide susceptibility areas.The connection between the tools was possible in the ArcGIS software environment by using the ESRI firm's ArcGIS 10.4 software.
For the landslide susceptibility mapping, field survey, high-resolution images, and aerial photo interpretation methods were used to make the landslide inventory map.A total of 15 conditioning factors were prepared, such as elevation, slope, aspect, TWI, STI, SPI, NDVI, distance to faults, distance to roads, distance to the river, distance to urban areas, rainfall, land cover/use, lithology, and soil type.The final landslide susceptibility maps were confirmed to have a very high accuracy of the validation represented by using the AUC value.
It is possible to draw the following conclusions.First, the proposed method can use both qualitative and quantitative methods to perform landslide susceptibility maps.Expert knowledge can fully be used, which differs from statistical learning.Second, the validation process was performed based on the comparison of historical landslide locations with the different landslide-susceptible zones on the final map.The higher AUC values demonstrate the effectiveness of the proposed method.Specifically, according to the resultant LSI scales, 776.3 km 2 (15%), 1126 km 2 (21.8%), and 1120 km 2 (21.7%) of the study area were labelled with "the highest", "high", and "moderate" susceptibility, respectively, which proves that the tested area is mostly located in the landslide-prone areas.These highly susceptible zones are located on the most rugged slopes, which are the area with the greatest rainfall and with a higher elevation and highland reliefs.
The main recommendations for the further application of this model are those pointing out its accessible mathematical apparatus and the stability (consistency) of its solution, as well as the possibility of combining it with other methods, especially in the section related to the determining of the weight criteria.A future improvement of the GIS MCDA model may be introduced by applying a multiplicative full consistency method [84,85] or the intervals of fuzzy numbers in order to determine the parameters of the conditioning factors.In addition to the fuzzy technique, grey theory is yet another suitable tool to be applied to the conditions of uncertainty and ambiguity.

Figure 1 .
Figure 1.The location of the study area (Mačva and Kolubara Districts, and the Tara National Park).

Figure 1 .
Figure 1.The location of the study area (Mačva and Kolubara Districts, and the Tara National Park).

Figure 2 .
Figure 2. The flowchart of the applied methodology

Figure 2 .
Figure 2. The flowchart of the applied methodology.

Step 5 .
Determining the BO matrix for the experts' average answers.Based on the BO matrices of the experts' answers A e B = a e Bj 1×n

Figure 4 .
Figure 4.The environmental factors related to landslides: (a) the soil type, (b) the distance to the river, (c) lithology, (d) the distance to faults, (e) the NDVI, (f) rainfall.

Figure 5 .
Figure 5.The social factors related to landslides: (a) the distance to roads, (b) the distance to urban areas, (c) the land use/cover.Figure 5.The social factors related to landslides: (a) the distance to roads, (b) the distance to urban areas, (c) the land use/cover.

Figure 5 .
Figure 5.The social factors related to landslides: (a) the distance to roads, (b) the distance to urban areas, (c) the land use/cover.Figure 5.The social factors related to landslides: (a) the distance to roads, (b) the distance to urban areas, (c) the land use/cover.

6 Figure 6 .
Figure 6.The final landslide susceptibility map created by applying the WLC aggregation method.

Figure 6 .
Figure 6.The final landslide susceptibility map created by applying the WLC aggregation method.
ν 2 , . . ., ν n ]-the set of the order weights Ai = [a i1 , a i2 . . ., a in ]-the set of the standardized criterion value z ij = [z i1 , z i2 . . ., z in ]-the sequence obtained by reordering the criterion values a i1 , a i2 . . ., a in .By using the reclassify tool in Spatial Analyst Tools ArcGIS 10.4 software, each cell of the final map is classified into five categories and receives a new value from one to five, representing the LSI.The results of the landslide hazard assessment are given in the maps shown in Figure 7. Value 1 is the area with the least probability of a landslide occurrence, whereas Value 5 represents the areas with the highest probability of a landslide occurrence.ISPRS Int.J. Geo-Inf.2019, 8 FOR PEER REVIEW 7 By using the reclassify tool in Spatial Analyst Tools ArcGIS 10.4 software, each cell of the final map is classified into five categories and receives a new value from one to five, representing the LSI.The results of the landslide hazard assessment are given in the maps shown in Figure 7. Value 1 is the area with the least probability of a landslide occurrence, whereas Value 5 represents the areas with the highest probability of a landslide occurrence.

Figure 7 .
Figure 7.The final landslide susceptibility map created by using the OWA aggregation method.

Figure 7 .
Figure 7.The final landslide susceptibility map created by using the OWA aggregation method.

Table 1 .
Summarized several landslide susceptibility and hazard map approaches.

Table 2 .
The data used in the susceptibility assessment, the data sources, and the associated factor classes for the landslide susceptibility mapping in the study area.

Table 3 .
The conditioning factors' description.

Table 4 .
The types of the geological formation of the study area.

Table 5 .
The description for the soil type.

Table 6 .
The description for the land cover/land use.

Table 7 .
The fuzzification of the conditioning factors.

Table 8 .
The Best-to-Others (BO) and the Others-to-Worst (OW) vectors obtained by the experts' evaluation.

Table 9 .
The average BO and OW vectors.

Table 10 .
The optimal values (weights) of the criteria.

Table 11 .
The consistency index and the consistency ratio of the best-worst methodology (BWM).