A Hybrid GIS Multi-Criteria Decision-Making Method for Flood Susceptibility Mapping at Shangyou, China

Abstract

Floods are considered one of the most disastrous hazards all over the world and cause serious casualties and property damage. Therefore, the assessment and regionalization of flood disasters are becoming increasingly important and urgent. To predict the probability of a flood, an essential step is to map flood susceptibility. The main objective of this work is to investigate the use a novel hybrid technique by integrating multi-criteria decision analysis and geographic information system to evaluate flood susceptibility mapping (FSM), which is constructed by ensemble of decision making trial and evaluation laboratory (DEMATEL), analytic network process, weighted linear combinations (WLC) and interval rough numbers (IRN) techniques in the case study at Shangyou County, China. Specifically, we improve the DEMATEL method by applying IRN to determine connections in the network structure based on criteria and to accept imprecisions during collective decision making. The application of IRN can eliminate the necessity of additional information to define uncertain number intervals. Therefore, the quality of the existing data during collective decision making and experts’ perceptions that are expressed through an aggregation matrix can be retained. In this work, eleven conditioning factors associated with flooding were considered and historical flood locations were randomly divided into the training (70% of the total) and validation (30%) sets. The flood susceptibility map validates a satisfactory consistency between the flood-susceptible areas and the spatial distribution of the previous flood events. The accuracy of the map was evaluated by using objective measures of receiver operating characteristic (ROC) curve and area under the curve (AUC). The AUC values of the proposed method coupling with the WLC fuzzy technique for aggregation and flood susceptibility index are 0.988 and 0.964, respectively, which proves that the WLC fuzzy method is more effective for FSM in the study area. The proposed method can be helpful in predicting accurate flood occurrence locations with similar geographic environments and can be effectively used for flood management and prevention.

1. Introduction

In recent decades, extreme weather events and meteorological disasters have frequently occurred in the context of global warming. Climate change is considered a critical factor in the growth and development of global warming [1,2,3]. As the most serious meteorological disasters, floods frequently occur around the world and cause casualties and property losses [4,5,6]. According to statistics, the current losses from floods among all natural disasters comprise 40%. The serious influences of floods on natural ecosystems and human activities have become an important factor to restrict sustainable development of societies and economies [4,7]. Therefore, the assessment and regionalization of flood disaster risks are becoming increasingly important and urgent. Although it is a very difficult task to prevent floods, we can predict and compensate for the disaster. To predict the probability of a flood, an essential step is to map flood susceptibility.

Flooding is a complicated phenomenon and there are many human factors and natural factors in the occurrence and development of floods. Among them, climate change plays an important role for extreme flood occurrence [8,9]. For instance, climate change may influence land use and result in flood hazards [10]. In recent years, the impact of climate change on flood risk and occurrence has been studied at a global scale [11,12]. Flood susceptibility mapping (FSM) can identify and predict future flood hazards based on statistical or deterministic methods. Mapping the areas that are susceptible to historic disaster locations is crucial for flood prevention and management [13]. Furthermore, it can be the basis for decision makers (DMs) to reduce the damage caused by disasters. To produce susceptibility maps, prediction methods combine different conditioning factors and weights the importance of these factors using subjective decision-making rules based on the experience of experts.

China has experienced severe floods since ancient times. In recent decades, the frequency of floods and disaster losses have greatly increased [12]. Jiangxi province is located in the south of China and characterized with a subtropical, warm, and humid monsoon climate. Specifically, this area is influenced by the continental monsoon climate, and nearly twenty climate-warming, extreme-weather and disaster-climate events occurred, especially floods, which have caused great losses to the national economy [14]. According to the Ministry of Water Resources of the People’s Republic of China, flood disasters have become very serious in recent years. For instance, the number of floods reached 31 during 2003–2008, and the total number of deaths was 114. Therefore, flood disaster risks assessment and regionalization must be conducted to avoid and reduce economic losses in the study area. Many scholars around the world have studied flood disaster risks and presented various methods to deal with the problem [7,15,16,17]. These commonly used methods can be mainly divided into three categories: (1) hydrological methods of WetSpa [18], SWAT [19] and HYDROTEL [18] and hydrodynamic approaches based on the shallow water equations initialized by rainfall [9,13,20,21,22,23]. However, the time cost of parameter setting and model construction in these methods is very high. (2) Statistical and data-driven methods. In the past decade, FSM has been conducted using remote sensing data and geographic information system (GIS). Examples of the GIS-based methods in flood contributions, including weights of evidence [24], logistic regression (LR) [25], analytic hierarchy process (AHP) [26], frequency ratios (FR) [27,28]. These methods are often based on statistical assumptions and the selection of flood conditioning factors is a challenging problem in FSM. (3) Machine learning methods which can cope with complex nonlinear problems without statistical assumptions [22], such as neuro-fuzzy logic [23], artificial neural networks (ANNs) [29,30], decision trees (DTs) [31], support vector machines (SVMs) [19,24,32]. The poor projection from nonlinear structures due to different data ranges in the datasets is the main shortcoming of such methods. Therefore, the optimal method to examine flood susceptibility is a complicated process and this problem is still debated. As mentioned previously, GIS has been widely used to perform flood modeling and risk mapping [31,33]. At present, the combination of multi-criteria decision analysis (MCDA) and GIS is a fundamental strategy to deal with the problem of risk reduction and flood management [1,12,26,34,35,36,37,38,39]. The integration of GIS and MCDA provides decision makers the single option for geographical region and different weight coefficients for different options. It means that this integration can create remarkable capabilities which complement each other. On the one hand, we can use the GIS environment to manipulate, store, manage, analyze and visualize spatial data. On the other hand, MCDA provides a rich mixture of procedures, techniques and algorithms to design and structure the decision problems and evaluates the alternative decisions [40]. Using MCDA-GIS for flood spatial prediction application can improve the understanding of uncertainties surroundings.

Recently, the integration of remote sensing (RS) and GIS have been used to map flooded areas [39,41,42,43,44,45]. It is suitable to perform flood detection using remote sensing images. For prediction, the detection is used to produce new ground truth data that are in agreement with the recent remote sensing data. Furthermore, these RS data can be used for image classification to obtain land use/cover information using well-known classifiers and the extraction of normalized difference vegetation index (NDVI), etc. Meanwhile, the risk analysis of floods is unimaginative without the support of GIS due to its powerful geostatistical tools. Following the recent contributions previously mentioned, we combine the two techniques for FSM. The main objective of this work is to investigate a novel hybrid MCDA-GIS technique to evaluate flood hazards, which is constructed by ensemble of decision making trial and evaluation laboratory (DEMATEL), analytic network process (ANP) and interval rough numbers (IRN) techniques. Specifically, we improve the DEMATEL method by applying IRN to determine connections in the network structure based on criteria and to accept imprecisions during collective decision making. The application of IRN can eliminate the necessity of additional information to define uncertain number intervals. Therefore, the quality of the existing data during collective decision making and experts’ perceptions that are expressed through an aggregation matrix can be retained. The flood susceptibility map is the first step in the development of flood risk management and is expected to be used by local governments for effective management planning purposes. In this paper, the FSM application of the IRN-DEMATEL-ANP method is presented in the case study at Shangyou County, China. To validate the effectiveness of the proposed method, the resultant susceptibility map was compared with historical flood locations. Moreover, the commonly used objective measures of receiver operating characteristic (ROC) curve and area under the curve (AUC) [46,47] were used for evaluation.

The remainder of the paper is structured as follows. Section 2 briefly describes the study area of Shangyou County, China and some available data. Section 3 introduces the proposed method and the phases of this method. Section 4 reports the evaluation of flood-prone areas at Shangyou County. The last section presents some concluding remarks and the future work.

2. Study Area and Data

2.1. Description of Study Area

The Shangyou area is located in Ganzhou in the southern Jiangxi Province. This area is located in the northern hilly district of Nanling between longitudes 114°00′ E and 114°40′ E and latitudes 25°42′ N and 26°01′ N. The area and population of the Shangyou district are almost 1543.0 km² and 3.22 million people, respectively. According to the Shangyou Meteorological Bureau, this district belongs to a subtropical monsoon climate. The average annual sunshine time is 1708.3 h, which is an intermediate level in Jiangxi Province, and the average annual temperature is 18.6 °C. January and July are the coldest and hottest months in each year, with an average of −2.7 °C and 38.0 °C, respectively. During 1959–2016, the average annual rainfall was 1483.4 mm and it was between 933.7 and 2147.6 mm from 1959 to 2014. The rainfall in the study area can greatly vary in spring and summer and the annual rainy season is usually from April to August. Specifically, the average annual rainfall from June to August is approximately 1055.8 mm and the largest daily rainfall can be over 105 mm during this time, while the monthly rainfall from September to January is only 69.3 mm, which indicates that the rainfall in these months is very small.

2.2. Available Data

To produce the flood inventory map is an important step in FSM. The optimal method to create flood inventory maps is still debated [30]. In this study, field surveys and remote sensing images were applied to create a flood inventory map, as shown in Figure 1, where “Nonflood” represents locations without historical flood occurrences. The triangles and circles in Figure 1 denote historical flood locations (HFLs) and locations without historical flood, respectively. Figure 2 illustrates some flood field photos of the study area: https://www.baidu.com. The selection of flood conditioning factors is another key issue that has been studied by many researchers. In this study, the factors of altitude, slope, plan curvature, topographic wetness index (TWI), sediment transport index (STI), NDVI, distance from rivers, rainfall, land cover/use, lithology, and soil were used to analyze the flood susceptibility. Detailed information regarding these flood conditioning factors can be found in

Table 1. Available data used in the susceptibility assessment, including the data sources and associated factor classes for FSM in the study area.

Sub-Classification		Source of Data	GIS Data Type		Scale or Resolution
Spatial database	Data layers		Spatial database	Derived map	Spatial database
Flood inventory	Flood inventory	Jiangxi Meteorological Bureau and Department of Civil Affairs of Jiangxi province	Point and polygon	-	-
Topographic map	Slope	ASTER GDEM Version 2	GRID	Slope gradient (in degrees)	30 m
	Elevation		GRID	Elevation	30 m
	TWI		GRID	Topographic wetness index	30 m
	SPI		GRID	Stream power index	30 m
	STI		GRID	Sediment transport index	30 m
River	Drainage network		ARC/INFO Line coverage	Drainage network	30 m
Soil	Soil	Institute of Soil Science, Chinese Academy of Sciences	Polygon	Soil	1:1,000,000
Geology Map	Lithology types	China Geology Organization	ARC/INFO coverage	Lithology	1:200,0000
Land-use type	Land use	Landsat 7 ETM + images	ARC/INFO GRID	Land use	30 m
Normalized difference vegetation index	NDVI		ARC/INFO GRID	NDVI	30 m
Rainfall	Rainfall	Jiangxi Meteorological Bureau	GRID	Precipitation map (mm)	1:50,000

.

3. Methodology

3.1. Methodological Background

Spatial MCDA consists several procedures, including the use of both geographic data, the preferences of the decision maker (DM), the manipulation of these data and preferences according to specified decision rules [40,48]. This method exploits the capabilities of GIS in the management of spatial data and the flexibility of MCDA to combine real spatial information (e.g., slope and land use) with value-based information (e.g., expert opinions, standards, and surveys) [49]. The main advantage of integrating GIS and MCDA is that their specific capabilities can complement each other. The methodological hierarchy of this paper is based on a MCDA-GIS structure. The proposed methodology to define flood hazard zones at Shangyou County, China includes the main steps from a methodological perspective as shown in Figure 3.

In this study, we used ArcGIS 10.3 software to construct and collect all flood conditioning factors in the GIS database of ESRI, which includes a completely new approach to the formulation of spatial databases. Additionally, all the GIS analyses, including production and classification of conditioning factor maps, conditioning factor standardization with fuzzy logic, conditioning factor map reclassification, weighted linear combinations (WLC), aggregation methods and creation of the final flood susceptibility maps, were conducted by using the ArcGIS 10.3 software. All these functions that were used in these GIS analyses are included in the 3D and spatial-analyst extensions. Meanwhile, FSI values and HFLs were produced from validation analysis and comparisons, and a database in Microsoft Excel with Visual Basic support and the XLSTAT Add-in were used for the DEMETAL-ANP method.

When solving real problems, the available conditioning factors have not the same contribution to flood occurrence, thus the DMs should define importance of each conditioning factor by using appropriate weight coefficients (weights). The WLC method requires normalization of the weights. After determining weight coefficients of the conditioning factors, the hybrid IRN-DEMATEL-ANP method is used to calculate the normalized weight conditioning factors, and then apply the WLC method to obtain final flood susceptibility maps.

3.2. IRN-DEMATEL-ANP Method

IRN-DEMATEL-ANP method is ensemble by IRN, DEMATEL and ANP methods, which is used to calculate the normalized weight conditioning factors. We introduce this method followed by the order of integration.

3.2.1. DEMATEL Method

The DEMATEL method was originally developed by the Science and Human Affairs Program of the Battelle Memorial Institute of Geneva to study complex and intertwined problematic groups [50]. It is a comprehensive method in both the design and analysis of structural methods that are characterized by causal relationships between complex factors [51]. The result of this method is a sum of the direct and indirect effects of all the factors that are transferred to and received by other factors. Meanwhile, the DEMATEL method is based on graph theory and used to identify the dependent factors and degree of dependence between them. Specifically, it enables the visual planning and solving problems so that all the relevant factors can be classified into causal and consequential factors to better understand their interrelations. Moreover, it can improve our understanding of the complex structure of the analysed problem and define relationships both between factors and between the structure level and strength of factors [52].

In this paper, we improve the DEMATEL method by applying IRN [53] to accept imprecisions during collective decision making. The application of IRN can eliminate the necessity of additional information to define uncertain number intervals [54,55,56,57,58]. Therefore, the quality of the existing data during collective decision making and experts’ perceptions that are expressed through an aggregation matrix can be retained. Previous studies proved that the original DEMATEL method can be modified to comply with their problems [59,60,61,62]. The improved IRN-DEMATEL method in this paper was motivated by the work of Gigović et al. [38]. We will introduce IRN-DEMATEL in Section 3.2.3. The basic idea of IRN is necessary to be illustrated first.

3.2.2. IRN Method

The interval rough numbers technique is used to deal with uncertainty. IRN consist of both lower and upper limit interval rough sequences, which give the range of this uncertainty. The lower limit and upper limit rough sequence contain detail object classes, and object classes are calculated by criterion value. This is a construction process from top to bottom.

The theory of IRN is introduced in this subsection and the calculation details are provided in Section 4.2. Before definition, the criterion number k and the expert number m are known. For each pair-wise criterion, each expert gives an interval $\left[z_{ij}^{eL},z_{ij}^{e\prime U}\right]$, where $z_{ij}$ indicates the degree of the ith criterion C_i affecting the jth criterion C_j, the superscripts eL and e’U represent the expert lower and upper interval limits, respectively. And the expert serial number e is from 1 to m. The object classes are defined as:

$$\underset{\_}{Lim}\left(z_{ij}^{eL}\right)=\frac{1}{M_L}{\displaystyle \sum z_{ij}^L}, \mathrm{subject} \mathrm{to} z_{ij}^L\le z_{ij}^{eL}$$

(1)

$$\overline{Lim}\left(z_{ij}^{eL}\right)=\frac{1}{M_U}{\displaystyle \sum z_{ij}^U}, \mathrm{subject} \mathrm{to} z_{ij}^U\ge z_{ij}^{eL}$$

(2)

where M_L and M_U are the number meeting the constraints, L and U are the index of expert meeting $L,U\in \left[1,m\right]$. It is worth noting that object classes count the numbers larger or less than the $z_{ij}^{eL}$ value given by the current expert e. For the $z_{ij}^{e\prime U}$ given by the current expert e, the definition form of $z_{ij}^{e\prime U}$ is similar to $z_{ij}^{eL}$, we just define $z_{ij}^{eL}$ as demonstration in this paper.

The rough sequences define as:

$$RN\left(z_{ij}^{eL}\right)=\left[\underset{\_}{Lim}\left(z_{ij}^{eL}\right),\overline{Lim}\left(z_{ij}^{eL}\right)\right]$$

(3)

$$RN\left(z_{ij}^{e\prime U}\right)=\left[\underset{\_}{Lim}\left(z_{ij}^{e\prime U}\right),\overline{Lim}\left(z_{ij}^{e\prime U}\right)\right]$$

(4)

IRN consists of this two rough sequences, and the IRN can be defined as

$$IRN\left(z_{ij}^e\right)=\left[RN\left(z_{ij}^{eL}\right),RN\left(z_{ij}^{e\prime U}\right)\right]$$

(5)

In other word, IRN contains a total of four values in two intervals.

3.2.3. IRN-DEMATEL Method

IRN-DEMATEL can provide the network relationship map among criteria, and the CER diagram is constructed by mapping all coordinate sets to visualize the complex interrelationship. For the subsequent analysis, the matrix definitions need to be stated. The direct-relation fuzzy matrix Z consists of the comparison of the pare-wise criteria value $z_{ij}$ provided by each expert, and the definition has completed before. The means of triangular fuzzy number ${\tilde{z}}_{ij}^e$ is the stack vector of min, average and max values that form the average matrix ${\tilde{Z}}^e$.

$${\tilde{z}}_{ij}^e=({\tilde{z}}_{ij,e}^{(l)},{\tilde{z}}_{ij,e}^{(m)},{\tilde{z}}_{ij,e}^{(r)})$$

(6)

where

$${\tilde{z}}_{ij,e}^{(l)}=\underset{M}{\min }\left\{z_{ij,e}^{(l)}\right\},M=\{1,2,\mathrm{....},e,\dots ,m\}$$

(7)

$${\tilde{z}}_{ij,e}^{(m)}=\frac{1}{m}{\displaystyle \sum _{t=1}^m{\tilde{z}}_{ij,t}^{(m)}}$$

(8)

$${\tilde{z}}_{ij,e}^{(r)}=\underset{M}{\max }\left\{z_{ij,e}^{(r)}\right\},M=\{1,2,\mathrm{....},e,\dots ,m\}$$

(9)

The normalized initial direct-relation matrix D is calculated by summing and normalize average matrix ${\tilde{Z}}^e$, the element $d_{ij}$ in the D defines as:

$$d_{ij}=\frac{{\tilde{z}}_{ij}^e}{Q}=(\frac{{\tilde{z}}_{ij,e}^{(l)}}{q^{(r)}},\frac{{\tilde{z}}_{ij,e}^{(m)}}{q^{(m)}},\frac{{\tilde{z}}_{ij,e}^{(r)}}{q^{(l)}})$$

(10)

$$Q=\max \left({\displaystyle {\sum }_{j=1}^n{\tilde{z}}_{ij}^e}\right)=\left(q^{(l)},q^{(m)},q^{(r)}\right)$$

(11)

Q is the sum of ${\tilde{Z}}^e$ and D can be separated into submatrix $\left(D_1,D_2,D_3\right)$. The total-relationship matrix is defined as:

$$T=\underset{w\to \propto }{\lim }\left(D+D^2+\dots +D^w\right)=D{\left(I-D\right)}^{-1}$$

(12)

where I is the identity matrix satisfied ${\lim }_{w\to \propto }(D+D_s+\dots D_s^2+D_s^k)={(I-D_s)}^{-1}$.

Each IRN consists of two rough sequences, and every rough sequence includes an upper and lower approximation, so the normalized matrix of average perception $D={\left[IRN(d_{ij})\right]}_{n\times n}$ can be divided into four sub-matrices, i.e., $D=\left(\left[D^L,D^U\right],\left[D^{\prime L},D^{\prime U}\right]\right)$, where $D^L={\left[\underset{\_}{Lim}(d_{ij})\right]}_{n\times n}$, $D^U={\left[\overline{Lim}(d_{ij})\right]}_{n\hspace{0.17em}\times n}$, $D^{\prime L}={\left[\underset{\_}{Lim}(d_{ij}^{\prime })\right]}_{n\times n}$ and $D^{\prime U}={\left[\overline{Lim}(d_{ij}^{\prime })\right]}_{n\hspace{0.17em}\times n}$. Moreover, $\underset{m\to \infty }{\lim }{\left(D^L\right)}^m=O\hspace{0.17em}$, $\underset{m\to \infty }{\lim }{\left(D^U\right)}^m=O\hspace{0.17em}$, $\underset{m\to \infty }{\lim }{\left(D^{\prime L}\right)}^m=O\hspace{0.17em}$ and $\underset{m\to \infty }{\lim }{\left(D^{\prime U}\right)}^m=O\hspace{0.17em}$, where O denotes a zero matrix. For a more effective explanation, we calculate in these four sub-parts:

$$\begin{array}{l}\underset{m\to \infty }{\lim }\hspace{0.17em}\hspace{0.17em}\left(I+D^L+D^{2L}+\dots +D^{mL}\right)={\left(I-D^L\right)}^{-1}\hspace{0.17em}\\ \underset{m\to \infty }{\lim }\hspace{0.17em}\hspace{0.17em}\left(I+D^U+D^{2U}+\dots +D^{mU}\right)={\left(I-D^U\right)}^{-1}\hspace{0.17em}\\ \underset{m\to \infty }{\lim }\hspace{0.17em}\hspace{0.17em}\left(I+D^{\prime L}+D^{2\prime L}+\dots +D^{m\prime L}\right)={\left(I-D^{\prime L}\right)}^{-1}\hspace{0.17em}\\ and\\ \underset{m\to \infty }{\lim }\hspace{0.17em}\hspace{0.17em}\left(I+D^{\prime U}+D^{2\prime U}+\dots +D^{m\prime U}\right)={\left(I-D^{\prime U}\right)}^{-1}\hspace{0.17em}\end{array}\}$$

(13)

Therefore, the matrix of the total influences T can be obtained by calculating the following elements:

$$\begin{array}{l}{T}^L=\underset{m\to \infty }{\lim }\hspace{0.17em}\hspace{0.17em}\left(I+D^L+D^{2L}+\dots +D^{mL}\right)={\left(I-D^L\right)}^{-1}={\left[\underset{\_}{Lim}(t_{ij}^L)\right]}_{n\times n}\hspace{0.17em}\\ T^U=\underset{m\to \infty }{\lim }\hspace{0.17em}\hspace{0.17em}\left(I+D^U+D^{2U}+\dots +D^{mU}\right)={\left(I-D^U\right)}^{-1}={\left[\underset{\_}{Lim}(t_{ij}^U)\right]}_{n\times n}\hspace{0.17em}\\ T^{\prime L}=\underset{m\to \infty }{\lim }\hspace{0.17em}\hspace{0.17em}\left(I+D^{\text{'}L}+D^{2\text{'}L}+\dots +D^{m\text{'}L}\right)={\left(I-D^{\text{'}L}\right)}^{-1}={\left[\underset{\_}{Lim}(t_{ij}^{\text{'}L})\right]}_{n\times n}\hspace{0.17em}\\ and\\ T^{\prime U}=\underset{m\to \infty }{\lim }\hspace{0.17em}\hspace{0.17em}\left(I+D^{\text{'}U}+D^{2\text{'}U}+\dots +D^{m\text{'}U}\right)={\left(I-D^{\text{'}U}\right)}^{-1}\hspace{0.17em}={\left[\underset{\_}{Lim}(t_{ij}^{\text{'}U})\right]}_{n\times n}\end{array}\}$$

(14)

where $D^L={\left[\underset{\_}{Lim}(d_{ij})\right]}_{n\times n}$, $D^U={\left[\overline{Lim}(d_{ij})\right]}_{n\hspace{0.17em}\times n}$, $D^{\prime L}={\left[\underset{\_}{Lim}(d_{ij}^{\prime })\right]}_{n\times n}$ and $D^{\prime U}={\left[\overline{Lim}(d_{ij}^{\prime })\right]}_{n\hspace{0.17em}\times n}$.

The sub-matrices $T^L$, $T^U$, $T^{\prime L}$ and $T^{\prime U}$ represent the rough-interval matrix of the total influences $T=\left(\left[T^L,T^U\right],\left[T^{\prime L},T^{\prime U}\right]\right)$. A total-relationship matrix can be defined based on Equations (13) and (14):

$$T=\left[\begin{array}{cccc}IRN(t_{11})& IRN(t_{12})& \cdots & IRN(t_{1n})\\ IRN(t_{21})& IRN(t_{22})& \cdots & IRN(t_{2n})\\ ⋮& ⋮& \ddots & ⋮\\ IRN(t_{n1})& IRN(t_{n2})& \cdots & IRN(t_{nn})\end{array}\right]$$

(15)

where $IRN(t_{ij})=\left[RN(t_{ij}^L),RN(t_{ij}^U)\right]$ is a IRN that is used to express the indirect effects of factor i on factor j. Then, the matrix T reflects the inter-dependence of each pair of factors.

The total-relationship matrix T, which are denoted as vectors R and C, respectively, and have rank n × 1:

$$IRN(R_i)={\left[{\displaystyle \sum _{j=1}^{n}IRN(t_{ij})}\right]}_{n\times 1}={\left[\left(\left[{\displaystyle {\sum }_{j=1}^n{t}_{ij}^L},{\displaystyle {\sum }_{j=1}^n{t}_{ij}^U}\right],\left[{\displaystyle {\sum }_{j=1}^n{t}_{ij}^{\prime L}},{\displaystyle {\sum }_{j=1}^n{t}_{ij}^{\prime U}}\right]\right)\right]}_{n\times 1}$$

(16)

$$IRN(C_i)={\left[{\displaystyle \sum _{i=1}^{n}IRN(t_{ij})}\right]}_{1\times n}={\left[\left(\left[{\displaystyle {\sum }_{i=1}^n{t}_{ij}^L},{\displaystyle {\sum }_{i=1}^n{t}_{ij}^U}\right],\left[{\displaystyle {\sum }_{i=1}^n{t}_{ij}^{\prime L}},{\displaystyle {\sum }_{i=1}^n{t}_{ij}^{\prime U}}\right]\right)\right]}_{1\times n}$$

(17)

The value R_i denotes the sum of the i-th row of matrix T and shows the total direct and indirect effects that criterion I delivers to other factors. Similarly, the value C_i is the sum of the j-th column of matrix T and represents the total direct and indirect effects that factor j receives from other factors. In cases where i = j, the equation (R_i + C_i) indicates the effect of the factors and the equation (R_i − C_i) indicates the intensity of the factors compared to others [35].

The threshold value α is calculated as the mean of all elements in T, Equation (18):

$$\alpha =\frac{{\displaystyle {\sum }_{i=1}^n{\displaystyle {\sum }_{j=1}^n\left[IRN(t_{ij})\right]}}}{N}$$

(18)

where N denotes the number of matrix elements (15).

And the CER diagram is constructed by mapping all coordinate sets of (R_i + C_j, R_i − C_j) to visualize the complex interrelationship. The procedure for this method is described in the following sections. The IRN-DEMATEL algorithm is shown as follows:

Algorithm 1: IRN-DEMATEL

Input: The expert pairwise comparison matrices Z Output: CER diagram

Step 1: Analysis of factors by experts. Step 2: Calculation of the average matrix ${\tilde{Z}}^e$. Step 3: A normalized initial direct-relationship matrix $D={[IRN\left(d_{ij}\right)]}_{n\times n}$ can be obtained based on the average matrix Z. Step 4: The total-relationship matrix $T={[IRN\left(t_{ij}\right)]}_{n\times n}$. Step 5: Calculate the sums of the rows Ri and columns Cj of the total-relationship matrix T [35]. Step 6: Set a threshold value (α) and construct a CER diagram.

3.2.4. IRN-DEMATEL-ANP Method

The ANP method considers the dependence and feedback of the criteria, which avoids the hierarchical constraints. In fact, ANP is a generalization of the AHP method. The calculation of the relative weights of criteria with traditional ANP mandates that the interdependence levels of factors are treated as reciprocal values. On the contrary, the interdependence levels of factors do not have reciprocal values when using the DEMATEL method, which is closer to real circumstances [63,64]. The following section deals with a novel approach, namely, the integration of the IRN-DEMATEL method into the ANP method, i.e., the IRN-DEMATEL-ANP method, which is shown as follows:

Algorithm 2: IRN-DEMATEL-ANP

Input: Total-relation matrix T Output: Weighted super matrix Wα

Step 1: Develop an unweighted super matrix. Step 2: Create a normalized total-influence matrix for criteria Tcα., which Tcα is the normalized matrix of T Step 3: Calculating the elements of the unweighted super matrix W, which is the transpose of Tcα. Step 4: Develop a weighted normalized super matrix Wα, which is TDα multiplication by W. The weighted normalized super matrix Wα can be calculated by the normalized influence matrix T with respect to the perspectives. TDα is the result of Tcα divided by dimension. Step 5: Find the limit of the weighted super matrix Wα, which multiply the weighted super matrix by itself multiple times results in the limit super matrix. The weight of each evaluation criterion is solved.

The IRN-DEMATEL construct the total-relation matrix T and we need ANP to confirm the final weight. Because the T is calculated in Algorithm 1. The Computation process of T is omitted and T is as the Input in Algorithm 2. A network model for the ANP method should be defined prior to the development of the unweighted super matrix. This network model can be defined based on the total-relationship matrix and CER diagram.

The unweighted super matrix is created when each level with the total degree of influence from the total-relationship matrix T is normalized by IR’DEMATEL. To normalize this matrix, we must determine the sum of the column elements of the matrix.

$$\begin{array}{ll}& \begin{array}{l}\begin{array}{cccc}{D}_1& D_2& & D_n\\ c_{11}{c}_{12}\dots c_{1m1}& c_{21}{c}_{22}\dots c_{2m2}& \dots & c_{n1}{c}_{n2}\dots c_{nmn}\end{array}\end{array}\\ \begin{array}{cc}& \\ T_c=\begin{array}{cc}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}{D}_1\\ \\ \end{array}\\ \\ \end{array}\\ D_2\\ \end{array}\\ \\ \\ D_n\end{array}& \begin{array}{c}\begin{array}{l}{c}_{11}\\ c_{12}\\ ⋮\end{array}\\ \begin{array}{l}{c}_{1m1}\\ c_{21}\\ c_{22}\end{array}\\ \begin{array}{l} ⋮\\ c_{2m2}\end{array}\\ \begin{array}{l}{c}_{n1}\\ c_{n2}\\ ⋮\\ c_{nmn}\end{array}\end{array}\end{array}& \end{array}& \left[\begin{array}{cccccccccccccccc}{T}_c^{11}& & & & & T_c^{12}& & & & & \dots & & & & & T_c^{1n}\\ & & & & & & & & & & & & & & & \\ & & & & & & & \\ & & & & & & & \\ T_c^{21}& & & & & T_c^{22}& & & & & \dots & & & & & T_c^{2n}\\ & & & & & & & & & & & & & & & \\ & & & & & & & & & & & & & & & \\ ⋮& & & & & ⋮& & & & & \ddots & & & & & ⋮\\ & & & & & & & & & & & & & & & \\ & & & & & & & & & & & & & & & \\ T_c^{n1}& & & & & T_c^{n2}& & & & & \dots & & & & & T_c^{nn}\end{array}\right]\end{array}$$

(19)

where the matrix $T_c^{11}$ contains factors from group D₁ and influences with respect to the factors from group D₁, the matrix $T_c^{12}$ (20) contains factors from the group (criteria) D₂ and influences with respect to the factors from group D₂, etc.

$$T_c^{12}=\left[\begin{array}{ccccc}IRN(t_{c^{12}}^{12})& \dots & IRN(t_{c^{1j}}^{12})& \dots & IRN(t_{c^{1m1}}^{12})\\ ⋮& & ⋮& \dots & ⋮\\ IRN(t_{c^{i1}}^{12})& \dots & IRN(t_{c^{ij}}^{12})& \dots & IRN(t_{c^{1m1}}^{12})\\ ⋮& & ⋮& & ⋮\\ IRN(t_{c^{m12}}^{12})& \dots & IRN(t_{c^{m1j}}^{12})& \dots & IRN(t_{c^{m1m1}}^{12})\end{array}\right]$$

(20)

Normalization of total-influence matrix for criteria T_c^α is conducted once T_c is developed. The criteria total-influence matrix T_c yields T_c^α after normalization. The normalized matrix T_c^α is shown below Equation (21):

$$\begin{array}{ll}& \begin{array}{l}\begin{array}{cccc}{D}_1& D_2& & D_n\\ c_{11}{c}_{12}\dots c_{1m1}& c_{21}{c}_{22}\dots c_{2m2}& \dots & c_{n1}{c}_{n2}\dots c_{nmn}\end{array}\end{array}\\ \begin{array}{cc}& \\ T_c^{\alpha }=\begin{array}{cc}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}{D}_1\\ \\ \end{array}\\ \\ \end{array}\\ D_2\\ \end{array}\\ \\ \\ D_n\end{array}& \begin{array}{c}\begin{array}{l}{c}_{11}\\ c_{12}\\ ⋮\end{array}\\ \begin{array}{l}{c}_{1m1}\\ c_{21}\\ c_{22}\end{array}\\ \begin{array}{l} ⋮\\ c_{2m2}\end{array}\\ \begin{array}{l}{c}_{n1}\\ c_{n2}\\ ⋮\\ c_{nmn}\end{array}\end{array}\end{array}& \end{array}& \left[\begin{array}{cccccccccccccccc}{T}_c^{\alpha 11}& & & & & T_c^{\alpha 212}& & & & & \dots & & & & & T_c^{\alpha 1n}\\ & & & & & & & & & & & & & & & \\ & & & & & & & \\ & & & & & & & \\ T_c^{\alpha 21}& & & & & T_c^{\alpha 22}& & & & & \dots & & & & & T_c^{\alpha 2n}\\ & & & & & & & & & & & & & & & \\ & & & & & & & & & & & & & & & \\ ⋮& & & & & ⋮& & & & & \ddots & & & & & ⋮\\ & & & & & & & & & & & & & & & \\ & & & & & & & & & & & & & & & \\ T_c^{\alpha n1}& & & & & T_c^{\alpha n2}& & & & & \dots & & & & & T_c^{\alpha nn}\end{array}\right]\end{array}$$

(21)

The total-influence matrix T_c describes the interdependence among the dimensions and criteria, so we can transpose the normalized total-influence matrix T_c^α by the dimensions based on the basic concept of ANP, resulting in the unweighted super matrix W = [T_c^α]’, as shown in Equation (22):

$$\begin{array}{ll}& \begin{array}{l}\begin{array}{cccc}{D}_1& D_2& & D_n\\ c_{11}{c}_{12}\dots c_{1m1}& c_{21}{c}_{22}\dots c_{2m2}& \dots & c_{n1}{c}_{n2}\dots c_{nmn}\end{array}\end{array}\\ \begin{array}{cc}& \\ W=(T_c^{\alpha })\prime =\begin{array}{cc}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}{D}_1\\ \\ \end{array}\\ \\ \end{array}\\ D_2\\ \end{array}\\ \\ \\ D_n\end{array}& \begin{array}{c}\begin{array}{l}{c}_{11}\\ c_{12}\\ ⋮\end{array}\\ \begin{array}{l}{c}_{1m1}\\ c_{21}\\ c_{22}\end{array}\\ \begin{array}{l} ⋮\\ c_{2m2}\end{array}\\ \begin{array}{l}{c}_{n1}\\ c_{n2}\\ ⋮\\ c_{nmn}\end{array}\end{array}\end{array}& \end{array}& \left[\begin{array}{cccccccccccccccc}{W}^{11}& & & & & W^{12}& & & & & \dots & & & & & W^{1n}\\ & & & & & & & & & & & & & & & \\ & & & & & & & \\ & & & & & & & \\ W^{21}& & & & & W^{22}& & & & & \dots & & & & & W^{2n}\\ & & & & & & & & & & & & & & & \\ & & & & & & & & & & & & & & & \\ ⋮& & & & & ⋮& & & & & \ddots & & & & & ⋮\\ & & & & & & & & & & & & & & & \\ & & & & & & & & & & & & & & & \\ W^{n1}& & & & & W^{n1}& & & & & \dots & & & & & W^{nn}\end{array}\right]\end{array}$$

(22)

where the matrix W¹¹ denotes the values of the factor influences from group D₁ in relation to the other factors from group D₁.

The elements of the weighted normalized super matrix W^α can be obtained by multiplying the elements from the unweighted super matrix W and the corresponding elements from the normalized total-influence matrix $T_D^{\alpha }$. The elements of the normalized total-influence matrix $T_D^{\alpha }$ can be obtained by normalizing the total-influence matrix T_D, as shown in Equation (23):

$$T_D^{\alpha }=\left[\begin{array}{ccccc}IRN(t_D^{\alpha 11})& \dots & IRN(t_D^{\alpha 1j})& \dots & IRN(t_D^{\alpha 1n})\\ ⋮& & ⋮& \dots & ⋮\\ IRN(t_D^{\alpha i1})& \dots & IRN(t_D^{\alpha ij})& \dots & IRN(t_D^{\alpha in})\\ ⋮& & ⋮& & ⋮\\ IRN(t_D^{\alpha n1})& \dots & IRN(t_D^{\alpha nj})& \dots & IRN(t_D^{\alpha nn})\end{array}\right]$$

(23)

where $IRN(t_D^{\alpha j})=IRN(t_D^j)/IRN(d_i)$ and the value of $IRN(d_i)$ can be obtained from $IRN(d_i)={\displaystyle {\sum }_{j=1}^{n}IRN(t_D^{ij})}$.

Once the elements of the matrix $T_D^{\alpha }$ are obtained, the elements of the new weighted super matrix W^α can be calculated. The elements of the matrix W^α can be obtained by multiplying the normalized total-influence matrix with dimensions $T_D^{\alpha }$ and unweighted super matrix W. It’s worth noting that the weighted super matrix W^α need multiply itself several time in Step 5. We multiply the rough weighted super matrix by itself multiple times to obtain the rough-limit super matrix, and then the weight of each evaluation criteria can be obtained. The rough-interval weighted super matrix can be raised to the limiting powers until the super matrix has converged and become a long-term stable super matrix to obtain the global priority vectors, or IRD’ANP influence weights, such as $\underset{k\to \infty }{\lim }=W^k$, where W denotes the limit super matrix and k denotes any number. Furthermore, we aggregate the criteria after determining their weight coefficients [52].

4. Results

The procedures of the conditioning factor selection and zoning of flood hazards at Shangyou County include six experts with experience in the fields of risk management, hydrology, spatial planning and environmental protection. These experts’ interviews were used to collect data that were further processed and their opinions were aggregated.

4.1. Conditioning Factor Selection

The selection of flood conditioning factors is very important and complex, and many scientists have different viewpoints regarding this topic. According to previous literature, topography, climate, human activity and soil are the key factors for occurrence and development of floods [24,25]. As mentioned in Section 2.2, eleven conditioning factors were selected to analyse the flood susceptibility. Figure 4C1–C11 shows all the factors that were uniformly transformed into a grid spatial size of 30 × 30 m and a grid of the Kelantan area was constructed.

4.2. MCDA-GIS Evaluation

This phase involves the standardization, experts’ work, weighting, summary analysis, and aggregation and validation of all the conditioning factors to be considered in the decision making. Since all the data were collected in various ways with different formats, the first step of MCDA is to normalize these data that can be used for comparison [38]. For this aim, the technique of the fuzzy concept was used according to the literature and experience of the experts. The fuzzy logic concept is flexible and suitable for modelling data in which there is no exact boundary of the set, determined by 0 or 1 [65]. In such cases, the affiliation of objects to a set is defined based on the degree of belonging to one of the functions of sigmoidal, J-shaped, linear or user-defined. The membership functions that are used often depend on the characteristics of the input data and the decisions and experience of experts. In this study, we used the discrete classification in which experts directly determined the values of the elements of fuzzy sets, as listed in

Table 2. Fuzzy standardization of the conditioning factors.

Conditioning Factors	Fuzzy Membership Function	Control Points/Value Points	Final Utility
Altitude (C1)	Linearly monotonically decreasing	c = 200 m; d = 800 m	0–200 m: equal to 1; 200–800 m: between 0 and 1; more than 800 m: equal to 0
Slope (C2)	Linearly monotonically decreasing	c = 5°; d = 25°	0°–5°: equal to 1; 5°–25°: between 0 and 1; more than 25°: equal to 0
Curvature (C3)	Linearly monotonically decreasing	c = -10; d = 10	0–−10: equal to 1; −10–10: between 0 and 1; more than 10: equal to 0
TWI (C4)	Linearly monotonically increasing	a = 4; b = 12	0–4: equal to 0; 4–12: between 0 and 1; more than 12: equal to 1
STI (C5)	Linearly monotonically decreasing	c = 1; d = 50	0–1: equal to 1; 1–50: between 0 and 1; more than 50: equal to 0
NDVI (C6)	Linearly monotonically decreasing	c = 0; d = 50	−1–0: equal to 1; 0–0, 6: between 0 and 1; more than 0, 6: equal to 0
Distance from river (C7)	Linearly monotonically decreasing	c = 100 m; d = 1000 m	0–100 m: equal to 1, 100–1000 m: between 0 and 1; more than 1000 m: equal to 0
Rainfall (C8)	Linearly monotonically increasing	a = 1000 mm; b = 2000 mm	0–1000 mm: equal to 0; 1000–2000 mm: between 0 and 1; more than 2000 mm: equal to 1
Land cover use (C9)	Discrete categorical data	Water = 1; Residential = 0.9; Bare soil = 0.7; Grass = 0.5; Farmland = 0.3; Forest = 0.1
Lithology (C10)	Discrete categorical data	D = 0.9; A = 0.8; I = 0.7; C = 0.6; G = 0.5; H = 0.4; B = 0.3; E = 0.2; F = 0.1
Soil (C11)	Discrete categorical data	WR = 0.9; ATc = 0.8; RGc = 0.7; CMo = 0.6; LVh = 0.5; ACh = 0.3; ACu = 0.2; Alh = 0.1

. It should be noted that elements of the conditioning factors of land cover/use, lithology and soil have categorical values. For the other conditioning factors, which represented values of gradual change between locations, the elements of the set were standardized by using the fuzzy concept based on the linear membership function. A scale ranging from zero to one byte was used for fuzzification, where zero indicates the lowest hazard and one indicates the most dangerous element of the set value in relation to the likelihood of flooding occurrence.

After standardization, the experts should define the significance measures of the conditioning factors by using the suitable coefficient weights (weights) or the conditioning factors weights. In this work, the IRN-DEMATEL method was used by the experts to analyse the factors and six experts were considered in this research. The following scale was used during evaluation: 1—very low influence, 2—low influence, 3—moderate influence, 4—high influence, and 5—very high influence. All the experts participated in the evaluation of the clusters and conditioning factors. Once the evaluation by the experts was completed, six pairwise comparison matrices were obtained, as shown in

Table 3. Expert comparison matrix in the pairwise conditioning factors.

Expert 1
	C1	C2	C3	C4	C5	C6	C7	C8	C9	C10	C11
C1	(0:0)	(3:5)	(3:5)	(2:5)	(2:4)	(3:4)	(3:5)	(2:5)	(2:4)	(1:4)	(4:4)
C2	(3:3)	(0:0)	(4:5)	(1:4)	(1:5)	(2:5)	(3:4)	(3:4)	(3:4)	(2:4)	(4:5)
C3	(5:5)	(3:3)	(0:0)	(4:4)	(4:5)	(3:5)	(2:4)	(3:4)	(3:3)	(2:5)	(4:4)
C4	(3:3)	(5:5)	(5:5)	(0:0)	(3:5)	(2:4)	(3:4)	(3:4)	(4:4)	(3:4)	(3:5)
C5	(3:5)	(3:3)	(3:5)	(3:5)	(0:0)	(2:5)	(3:4)	(3:4)	(4:5)	(1:4)	(3:5)
C6	(4:4)	(4:4)	(3:5)	(2:5)	(2:5)	(0:0)	(2:4)	(2:4)	(5:5)	(2:5)	(4:5)
C7	(4:4)	(4:4)	(4:5)	(2:5)	(2:4)	(2:3)	(0:0)	(4:4)	(2:4)	(2:4)	(3:4)
C8	(4:4)	(4:5)	(5:5)	(1:5)	(1:5)	(4:4)	(2:5)	(0:0)	(2:4)	(1:4)	(3:4)
C9	(3:3)	(2:2)	(4:4)	(1:4)	(1:5)	(3:5)	(4:4)	(2:4)	(0:0)	(2:4)	(3:3)
C10	(4:5)	(2:2)	(2:4)	(5:5)	(2:5)	(2:4)	(4:4)	(2:5)	(2:4)	(0:0)	(3:4)
C11	(1:2)	(1:3)	(2:4)	(2:4)	(1:4)	(3:4)	(3:4)	(3:5)	(2:4)	(1:4)	(0:0)
Expert 6
	C1	C2	C3	C4	C5	C6	C7	C8	C9	C10	C11
C1	(0:0)	(2:4)	(2:3)	(1:3)	(1:5)	(2:5)	(2:3)	(1:3)	(1:5)	(2:5)	(3:4)
C2	(3:4)	(0:0)	(4:5)	(1:5)	(1:3)	(2:3)	(4:4)	(2:3)	(2:5)	(2:5)	(3:4)
C3	(3:5)	(3:5)	(0:0)	(3:5)	(3:3)	(3:4)	(2:3)	(2:3)	(2:5)	(3:4)	(4:5)
C4	(3:4)	(4:5)	(4:5)	(0:0)	(2:4)	(2:3)	(3:3)	(4:4)	(4:5)	(2:3)	(2:4)
C5	(4:5)	(3:5)	(3:4)	(2:4)	(0:0)	(1:3)	(4:5)	(4:5)	(3:4)	(2:5)	(3:4)
C6	(3:5)	(3:4)	(3:4)	(1:4)	(1:3)	(0:0)	(2:3)	(1:5)	(4:5)	(3:3)	(4:4)
C7	(3:5)	(3:5)	(4:4)	(2:4)	(2:3)	(2:4)	(0:0)	(3:5)	(1:5)	(3:5)	(3:3)
C8	(3:4)	(4:4)	(4:5)	(1:4)	(1:3)	(4:5)	(1:3)	(0:0)	(1:5)	(2:5)	(4:5)
C9	(2:4)	(1:3)	(4:5)	(1:5)	(1:3)	(2:4)	(3:5)	(1:3)	(0:0)	(3:3)	(2:5)
C10	(4:5)	(1:3)	(1:5)	(4:5)	(1:4)	(1:3)	(4:5)	(1:3)	(1:3)	(0:0)	(4:4)
C11	(1:3)	(1:5)	(1:3)	(1:5)	(1:3)	(4:4)	(2:3)	(4:5)	(1:5)	(2:5)	(0:0)

.

Table 3 lists expert comparison matrix in the pairwise conditioning factors. It can be seen that the values of i and j are different, thus the experts expressed uncertainty when defining the influences of the conditioning factors during evaluation. In accord with the implementation of the IRN-DEMATEL method, the initial comparison matrices with pairwise conditioning factors are transformed into interval rough. The IRN consists of two rough sequences. In the following, we show the formation of rough individual sequences for a single position in conditioning factor matrices. The determination of the elements of the interval rough comparison matrix elements z¹, z², …, z⁶ representing expert in is demonstrated by calculating elements at the position C2–C3. Two rough sequences at the position C2–C3 that constitute IRN are obtained for each matrix $z^{\mathrm{e}}$. Two classes of objects $z_{ij}^L$ and $z_{ij}^U$ are defined for the position C2–C3 from the comparison matrices. Each class includes six elements:

$$z_{2,3}^L=\text{{}4;5;4;4;3;4\text{}}$$

$$z_{2,3}^U=\text{{}5;5;5;4;4;5\text{}}$$

Applying Equations (1) and (2), we can produce rough sequences for every object class. For the first class, we obtain

Expert 1 rough sequence for C2-C3	$\underset{\_}{Lim}\left(z_{2,3}^{1L}\right)=\underset{\_}{Lim}\left(4\right)=\frac{1}{5}\left(4+4+4+3+4\right)=3.8$
	$\overline{Lim}\left(z_{2,3}^{1L}\right)=\overline{Lim}\left(4\right)=\frac{1}{5}\left(4+5+4+4+4\right)=4.2$
Expert 2 rough sequence for C2-C3	$\underset{\_}{Lim}\left(z_{2,3}^{2L}\right)=\underset{\_}{Lim}\left(5\right)=\frac{1}{6}\left(4+5+4+4+3+4\right)=4.0$
	$\overline{Lim}\left(z_{2,3}^{2L}\right)=\overline{Lim}\left(5\right)=\frac{1}{1}\left(5\right)=5.0$
…	…
	…
Expert 6 rough sequence for C2-C3	$\underset{\_}{Lim}\left(z_{2,3}^{6L}\right)=\underset{\_}{Lim}\left(4\right)=\frac{1}{5}\left(4+4+4+3+4\right)=3.8$
	$\overline{Lim}\left(z_{2,3}^{6L}\right)=\overline{Lim}\left(4\right)=\frac{1}{5}\left(4+5+4+4+4\right)=4.2$

For second object class,

Expert 1 rough sequence for C2-C3	$\underset{\_}{Lim}\left(z_{2,3}^{1U}\right)=\underset{\_}{Lim}\left(5\right)=\frac{1}{6}\left(5+5+5+4+4+5\right)=4.67$
	$\overline{Lim}\left(z_{2,3}^{1U}\right)=\overline{Lim}\left(5\right)=\frac{1}{4}\left(5+5+5+5\right)=5.0$
Expert 2 rough sequence for C2-C3	$\underset{\_}{Lim}\left(z_{2,3}^{2U}\right)=\underset{\_}{Lim}\left(5\right)=\frac{1}{6}\left(5+5+5+4+4+5\right)=4.67$
	$\overline{Lim}\left(z_{2,3}^{2U}\right)=\overline{Lim}\left(5\right)=\frac{1}{4}\left(5+5+5+5\right)=5.0$
…	…
	…
Expert 6 rough sequence for C2-C3	$\underset{\_}{Lim}\left(z_{2,3}^{6U}\right)=\underset{\_}{Lim}\left(5\right)=\frac{1}{6}\left(5+5+5+4+4+5\right)=4.67$
	$\overline{Lim}\left(z_{2,3}^{6U}\right)=\overline{Lim}\left(5\right)=\frac{1}{4}\left(5+5+5+5\right)=5.0$

Thus, rough sequences that constitute IRN can be obtained:

Expert 1 IRN for C2-C3	$RN(z_{2,3}^{1L})=\left[3.8,4.2\right]$,$RN(z_{2,3}^{1U})=\left[4.67,5.0\right]$
	$\Rightarrow IRN(z_{2,3}^1)=\left(\left[3.8,4.2\right],\left[4.67,5.0\right]\right)$
Expert 2 IRN for C2-C3	$RN(z_{2,3}^{2L})=\left[4.0,5.0\right]$,$RN(z_{2,3}^{2U})=\left[4.67,5.0\right]$
	$\Rightarrow IRN(z_{2,3}^2)=\left(\left[4.0,5.0\right],\left[4.67,5.0\right]\right)$
…	…
	…
Expert 6 IRN for C2-C3	$RN(z_{2,3}^{6L})=\left[3.8,4.2\right]$,$RN(z_{2,3}^{6U})=\left[4.67,5.0\right]$
	$\Rightarrow IRN(z_{2,3}^6)=\left(\left[3.8,4.2\right],\left[4.67,5.0\right]\right)$

The IRN of other comparison matrices are obtained by applying the same method as shown in

Table 4. Interval rough expert comparison matrix in the pairwise conditioning factors.

Expert 1
	C1	C2	C3	C4	C5	...	C11
C1	[(0,0),(0,0)]	[(2.5,3),(3,67,4,33)]	[(2.5,3),(3,3.75)]	[(1.5,2),(3,3.75)]	[(1.5,2),(4.25,5)]	...	[(3,4),(4,4.25)]
C2	[(3,4),(3.75,4)]	[(0,0),(0,0)]	[(3.8,4.2),(4.47,5)]	[(1,1,25),(3,75,5)]	[(1,25,2),(3,4)]		[(3,5,4),(3,67,4,33)]
C3	[(3,5),(4.25,5)]	[(3,3.5),(3,5)]	[(0,0),(0,0)]	[(3.25,4),(4.25,5)]	[(3,3.25),(3,4.25)]		[(4,4),(4.75,5)]
C4	[(3,3.5),(3,4)]	[(4,5),(4.75,5)]	[(4.25,5),(4.75,5)]	[(0,0),(0,0)]	[(2,2.25),(3.5,4.67)]		[(2.67,3.33),(4,4.5)]
C5	[(4,4.5),(3,5)]	[(3,3.25),(3.5,5)]	[(3,3.25),(3.67,4.33)]	[(2.25,3),(4,4.5)]	[(0,0),(0,0)]		[(3,3.25),(3.67,4.33)]
C6	[(3.5,5),(3.75,4)]	[(4,4.25),(3,4)]	[(3,3.25),(4,4.5)]	[(1.25,2),(3.67,4.33)]	[(1,1.25),(3.4,25)]		[(4,4.25),(4,4.5)]
C7	[(3.25,4),(4.25,5)]	[(2,75,4),(4.25,5)]	[(3.33,4.33),(4,4.25)]	[(1,2),(4,4.5)]	[(2,2),(3,4)]		[(3,3.5),(3,3.5)]
C8	[(3,4),(3.5,4)]	[(3.67,4.33),(4,4.5)]	[(3.5,5),(4.25,5)]	[(1,1),(3.33,4.33)]	[(1,1),(2.5,4)]		[(3,4),(4.5,5)]
C9	[(2.75,4),(2.75,4)]	[(1.25,2),(2,3)]	[(4,4.5),(4.5,5)]	[(1,1.5),(3.5,5)]	[(1.5,2),(2.5,4)]		[(2.5,3),(3.25,5)]
C10	[(3.67,4.33),(4.5,5)]	[(1.5,2),(2.25,3)]	[(1.25,2),(3.5,5)]	[(4.25,5),(4.5,5)]	[(1,1.25),(2.67,4.5)]		[(3,3.25),(3.75,4)]
C11	[(1,1.5),(2,3)]	[(1,1.25),(2.75,5)]	[(1.5,2),(2.33,3.5)]	[(1.25,2),(3.5,5)]	[(1,1),(2.33,3.33)]		[(0,0),(0,0)]
Expert 6
	C1	C2	C3	C4	C5	...	C11
C1	[(0,0),(0,0)]	[(2,2.5),(4,5)]	[(2,2.5),(3.75,5)]	[(1,1.5),(3.75,5)]	[(1,1.5),(3.5,4.67)]	...	[(2.67,3.33),(4,4.25)]
C2	[(2.67,3.33),(3,3.75)]	[(0,0),(0,0)]	[(3.8,4.2),(4.67,5)]	[(1,1.25),(3.33,4.33)]	[(1,1.25),(4,5)]		[(3,3.5),(4,5)]
C3	[(2.33,3.67),(4,4.25)]	[(3,3.5),(2.33,3.67)]	[(0,0),(0,0)]	[(3,3.25),(3.5,4.67)]	[(3,3.25),(4.25,5)]		[(4,4),(4,4.75)]
C4	[(2.67,3.33),(3,3.5)]	[(4,4.75),(3.67,4.33)]	[(3.5,4.67),(4,4.75)]	[(0,0),(0,0)]	[(2,2.25),(4.25,5)]		[(2,3),(4.5,5)]
C5	[(3.67,4.33),(4.5,5)]	[(2.5,4),(3,3.25)]	[(3,3.25),(4,5)]	[(2,2.25),(4.5,5)]	[(0,0),(0,0)]		[(3,3.25),(4,5)]
C6	[(2.5,4),(3,3.75)]	[(2.67,3.33),(4,4.25)]	[(3,3.25),(4.5,5)]	[(1,1.25),(4,5)]	[(1,1.25),(4.25,5)]		[(4,4.25),(4.5,5)]
C7	[(2.5,3.67),(4,4.25)]	[(2.33,3.33), (4,4.25)]	[(3.75,5),(4,4.25)]	[(2,2),(4.5,5)]	[(2,2),(3.67,4.33)]		[(3,3.5),(3.5,4)]
C8	[(2.67,3.33),(3,3.5)]	[(4,4.5),(4,5)]	[(3.4,5),(4,4.25)]	[(1,1),(3.75,5)]	[(1,1),(3.5,5)]		[(3.67,4.33),(4,4.5)]
C9	[(2.33,3.33),(2,2.75)]	[(1,1.25),(1.67,2.33)]	[(4,4.5),(4,4.5)]	[(1,1.5),(3,4.5)]	[(1.5,2),(3.5,5)]		[(2,2.5),(2.67,3.67)]
C10	[(4,5),(4,4.5)]	[(1,1.5),(1.5,2.67)]	[(1,1.25),(3,4.5)]	[(4,4.25),(4,4.5)]	[(1,1.25),(3.25,5)]		[(3.25,4),(3.75,4)]
C11	[(1,1.5),(1.67,2.33)]	[(1,1.25),(2,4)]	[(1,1.5),(2.75,4)]	[(1,1.25),(3,4.5)]	[(1,1),(2.75,4)]		[(0,0),(0,0)]

as follows.

The interval rough matrices that refer to the responses are aggregated on the next level of the IRN-DEMATEL method. The mean IRN values are obtained based on the conditioning factor response matrices from Table 3 according to Equations (6)–(9). Therefore, the interval rough average matrix of can be obtained as listed in

Table 5. Interval rough average matrix of the conditioning factors.

	C1	C2	C3	C4	...	C11
C1	[(0,0),(0,0)]	[(2.25,2.75),(3.83,4.67)]	[(2.25,2.75),(3.38,4.38)]	[(1.25,1.75),(3.38,4.38)]	...	[(2.83,3.67),(4,4.25)]
C2	[(2.83,3.67),(3.5,3.75)]	[(0,0),(0,0)]	[(3,3.5),(3.54,4.67)]	[(1,1.25),(3.54,4.67)]		[(3.25,3.75),(3.83,4.67)]
C3	[(2.67,4.33),(4.13,4.63)]	[(2.67,4.33),(3,3.5)]	[(0,0),(0,0)]	[(3.13,3.63),(3.88,4.83)]		[(4,4),(4.38,4.88)]
C4	[(2.83,3.67),(3,3.5)]	[(4.38,4.88),(3.83,4.67)]	[(4.38,4.88),(3.88,4.83)]	[(0,0),(0,0)]		[(2.33,3.17),(4.25,4.75)]
C5	[(3.83,4.67),(3.75,4.75)]	[(3,3.25),(3,4.5)]	[(3,3.25),(3.83,4.67)]	[(2.13,2.63),(4.25,4.75)]		[(3,3.25),(3.83,4.67)]
C6	[(3.4,5),(3.38,3.88)]	[(2.83,3.67),(4,4.25)]	[(3,3.25),(4.25,4.75)]	[(1.13,1.63),(3.83,4.67)]		[(4,4.25),(4.25,4.75)]
C7	[(2.88,3.83),(4.13,4.63)]	[(2.54,3.67),(4.13,4.63)]	[(3.54,4.67),(4,4.25)]	[(1,2),(4.25,4.75)]		[(3,3.5),(3.25,3.75)]
C8	[(2.83,3.67),(3.25,3.75)]	[(3.83,4.67),(4,4.5)]	[(3.25,4.75),(4.13,4.63)]	[(1,1),(3.54,4.67)]		[(3.33,4.17),(4.25,4.75)]
C9	[(2.38,2.88),(2.54,3.67)]	[(1.13,1.63),(1.83,2.67)]	[(4.25,4.75),(4,4.5)]	[(1,1.5),(3.25,4.75)]		[(2.25,2.75),(2.96,4.33)]
C10	[(4.25,4.75), (3.83,4,67)]	[(1.25,1.75),(1.88,2.83)]	[(1.13,1.63),(3.25,4.75)]	[(4.13,4.63),(4.25,4.75)]		[(3.13,3.63),(3.75,4)]
C11	[(1,1.5),(1.83,2.67)]	[(1,1.25),(2.38,4.5)]	[(1.25,1.75),(2.54,3.75)]	[(1.13,1.63),(3.25,4.75)]		[(0,0),(0,0)]

as follows.

Once the average matrix of the conditioning factors is obtained, the determination of the initial direct-relationship matrix can be conducted. In this step, the initial direct-relationship matrix of the conditioning factors is transformed into the total-relationship matrix of the conditioning factors, as shown in

Table 6. The total relationship matrix of the conditioning factors.

	C1	C2	C3	C4	...	C11
C1	[(0.22,0.35),(0.56,0.78)]	[(0.29,0.36),(0.62,0.94)]	[(0.29,0.36),(0.70,1.00)]	[(0.17,0.24),(0.72,1.02)]		[(0.34,0.46),(0.76,1.01)]
C2	[(0.34,0.49),(0.65,0.86)]	[(0.27,0.34),(0.55,0.87)]	[(0.36,0.43),(0.72,1.03)]	[(0.19,0.26),(0.74,1.05)]		[(0.41,0.53),(0.78,1.04)]
C3	[(0.39,0.57),(0.68,0.88)]	[(0.41,0.47),(0.62,0.95)]	[(0.33,0.40),(0.66,0.95)]	[(0.28,0.35),(0.77,1.05)]		[(0.47,0.58),(0.80,1.04)]
C4	[(0.43,0.60),(0.66,0.85)]	[(0.48,0.55),(0.65,0.94)]	[(0.49,0.56),(0.75,1.01)]	[(0.22,0.29),(0.69,0.95)]		[(0.55,0.67),(0.83,1.03)]
C5	[(0.41,0.58),(0.71,0.91)]	[(0.40,0.48),(0.67,0.99)]	[(0.41,0.48),(0.79,1.07)]	[(0.25,0.33),(0.82,1.09)]		[(0.47,0.61),(0.86,1.08)]
C6	[(0.37,0.55),(0.67,0.87)]	[(0.40,0.47),(0.64,0.94)]	[(0.39,0.45),(0.76,1.03)]	[(0.21,0.29),(0.77,1.05)]		[(0.43,0.57),(0.82,1.04)]
C7	[(0.37,0.54),(0.67,0.85)]	[(0.41,0.49),(0.61,0.90)]	[(0.42,0.48),(0.73,0.99)]	[(0.21,0.30),(0.76,1.01)]		[(0.44,0.57),(0.79,0.99)]
C8	[(0.35,0.51),(0.62,0.85)]	[(0.39,0.47),(0.60,0.94)]	[(0.40,0.47),(0.69,1.02)]	[(0.20,0.27),(0.71,1.04)]		[(0.43,0.56),(0.75,1.02)]
C9	[(0.32,0.49),(0.56,0.79)]	[(0.29,0.38),(0.52,0.85)]	[(0.37,0.45),(0.65,0.96)]	[(0.19,0.27),(0.65,0.97)]		[(0.39,0.54),(0.68,0.97)]
C10	[(0.36,0.51),(0.62,0.83)]	[(0.30,0.38),(0.55,0.86)]	[(0.30,0.38),(0.67,0.97)]	[(0.27,0.34),(0.71,0.98)]		[(0.40,0.52),(0.72,0.97)]
C11	[(0.35,0.52),(0.61,0.79)]	[(0.37,0.44),(0.63,0.90)]	[(0.35,0.42),(0.72,0.96)]	[(0.22,0.30),(0.73,0.98)]		[(0.37,0.49),(0.70,0.90)]

.

The values of total direct and indirect effects that the jth conditioning factor has received from other conditioning factors and transferred to others can be obtained by summing the elements of the total-relationship matrix of the clusters/conditioning factors at each row and column. These values and the threshold value (α) of the total-relationship matrix are used to define CER diagram as shown in Figure 5, which is created to visualize the complicated causal relationship of the conditioning factors into a visible structural model.

The elements of matrix T in Table 5 with values above the threshold value α are identified and mapped on the diagram, as shown in Figure 5, where the x–axis and denotes IRN C_i + R_i and the y–axis denotes IRN C_i − R_i. These values are used to demonstrate the relationship between two factors. The arrows that denote a cause-and-effect membership are oriented from factors with values below α to elements with values above α.

The weight coefficients of the conditioning factors were input into the IRN-ANP method based on the CER diagram. First, the elements of the interval-rough unweighted and weighted super matrices were calculated based on the total-relationship matrix of the conditioning factors. The total-influence matrix T was included in the ANP method and the unweighted and weighted super matrices were obtained. The influential weights of the stable matrix were defined once the unweighted and weighted super matrices were computed. The weights of the conditioning factors were obtained based on the values in the weighted super matrix, as shown in

Table 7. Weight coefficients of the conditioning factors.Conditioning factor.

	Weight Coefficient	Rank	Crisp Weight Coefficient
Altitude (C1)	[(0.034,0.331),(0.040,0.337)]	4	0.1002
Slope (C2)	[(0.079,0.253),(0.071,0.272)]	6	0.0910
Curvature (C3)	[(0.034,0.303),(0.031,0.357)]	5	0.0972
TWI (C4)	[(0.015,0.237),(0.026,0.320)]	7	0.0790
STI (C5)	[(0.020,0.142),(0.053,0.186)]	9	0.0538
NDVI (C6)	[(0.023,0.224),(0.023,0.257)]	8	0.0708
Distance from rivers (C7)	[(0.022,0.368),(0.011,0.497)]	3	0.1184
Rainfall (C8)	[(0.091,0.378),(0.089,0.380)]	2	0.1266
Land cover/use (C9)	[(0.093,0.501),(0.091,0.656)]	1	0.1776
Lithology (C10)	[(0.012,0.107),(0.013,0.139)]	11	0.0360
Soil (C11)	[(0.023,0.153),(0.023,0.170)]	10	0.0496

.

To validate the effectiveness of the proposed method, two popular methods that were proposed recently were selected for comparison, including the crisp DEMTEL-ANP [29] and fuzzy DEMTEL-ANP methods [35]. A symmetric form of triangular fuzzy numbers was used to calculate the weight coefficients by applying the DEMTEL-ANP method. The comparison results are shown in Figure 6. All the three methods generated sequences of weight coefficients that were characterized by similar rank of C9 > C8 > C7 > C1 > C3 > C2 > C4 > C6 > C5 > C11 > C10 with different values.

Figure 6 shows each interval number with two colour shades (dark and light): the darker shade denotes the upper and lower ranges, while the lighter shade means the intersection of two rough sequences of the IRN. The crisp DEMTEL-ANP method computes weight coefficients by applying crisp numbers. Thus, any uncertainty and vagueness from the group decision making process can be neglected. Meanwhile, uncertainties from group decision making in the fuzzy DEMTEL-ANP and IRN-DEMTEL-ANP methods are represented by various dimensions of fuzzy and IRN of weight coefficients. These interval values are the result of various mechanisms that are employed to treat uncertainty and subjectivity. However, the fuzzy DEMTEL-ANP method handles uncertainty by means of fuzzy sets with previously defined boundaries, which cannot be extended or narrowed, the boundaries in IRN are flexible and can be adjusted to accommodate uncertainties in the data. Moreover, the previously defined boundaries in the fuzzy DEMTEL-ANP method increase the subjectivity in the group-decision-making process because the boundaries are defined based on subjective assessments. This phenomenon can significantly affect the degree of uncertainty, which is represented in the size of an interval, unlike in the IRN-DEMTEL-ANP method. Therefore, the proposed method can efficiently measure uncertainties during conditioning factor evaluation and reflect the perception of a decision maker.

4.3. Aggregation of Weighted Linear Combinations

WLC is the most commonly used method that is compensatory, meaning that a low result in one conditioning factor can be compensated by high results in another, which is desirable for this specific decision problem. Thus, WLC was selected as an aggregation method and it multiplied each fuzzy standardized conditioning factor map (i.e., each raster cell of 30 × 30 m) with the weights of the conditioning factors, obtaining different variations from the AHP method, and then summarizing the results. The mathematical expression for calculating the convenience index in the WLC method is as follows:

$$S={\displaystyle \sum w_i{x}_i}$$

(24)

where S is the suitability index, w_i is the normalized value of the factor weight and x_i is the criterion score of factor i. The WLC method was used to create a map of the aggregation of the conditioning factors on the final map, which is displayed in the same range of fuzzy values from 0 to 1, based on the employed conditioning factors and determination of their weights according to the scenario in the previous section. Finally, the FSI scales were calculated using the defuzzification algorithm of the standard deviation method using the Reclass Spatial Analyst Tool in ArcGIS 10.3 software. Each cell is classified into five categories and has been assigned a new value from 1 to 5, representing the FSI scales. The final results of the flood hazard assessment are provided in Figure 7.

Table 8. Areas of the classes from the final flood susceptibility map.

FSI		Area		Number of Cells (30 × 30 m)
		(km2)	%
FSI 1	Very low	118.4	7.7	131,544
FSI 2	Low	434.8	28.2	483,080
FSI 3	Moderate	588.4	38.1	653,747
FSI 4	High	329.8	21.4	366,514
FSI 5	Very high	71.6	4.6	79,534

shows the five classes of the final flood susceptibility map with the covered area of the study territory. According to the results in Table 8, the area with the largest FSI scale (FSI 5) was 71.6 km² or 4.6% of the study area. Additionally, 329.8 km² in the study area exhibited a high FSI scale (FSI 4). These areas were mainly located in the eastern and southern regions, which have lower altitudes, lowland relief and areas near river flows.

4.4. Model Validation

The validation of flood susceptibility maps is an essential step during modelling. The prediction performance of the WLC method and the FSI scales were evaluated by using a non-dependent threshold approach of ROC [28]. The AUC is a synthetic index that is calculated for ROC curves and has been generally used to evaluate the accuracy of flood susceptibility maps [25,27]. The AUC is the probability that a positive event is classified as positive by the test given all the possible values of the test [66,67,68,69,70]. In addition, the model of this paper is established by expert knowledge, all historical flood samples are considered.

First, historical flood locations were used as the validation set to evaluate the accuracy of the prediction methods. The validation samples that represented flooded and non-flooded areas were integrated into a single vector. We filled the attributes of these validation samples with the values of the standardized conditioning factors and the values of the final maps from the WLC fuzzy method and the FSI scales. Histograms of the final map values showed a high correlation between the locations of the validation samples and the flood and non-flood areas. For the WLC fuzzy method, the maximum and minimum values for the flood validation samples were 0.874 and 0.582, respectively, and the standard deviation was 0.046, as shown in Figure 8a. In Figure 8b, the maximum and minimum values for the non-flood validation samples were 0.597 and 0.182, respectively, and the standard deviation was 0.106. In addition, the FSI scales in the final map ranged from 3.114 to 5 and 1 to 3.508 for the flood and non-flood validation samples as shown in Figure 8c,d, respectively. The ROC curve corresponds to the graphical representation of the specificity and sensitivity for the various possible threshold values. The ROC curves and specificity/sensitivity diagrams for the WLC fuzzy model and FSIs are shown in Figure 9.

The final results of the eleven flood susceptibility maps (maps generated through IRN-DEMANTEL-ANP method) were validated by comparing them with the historical flood locations, using the success-rate and prediction-rate methods. The success-rate results were obtained based on a comparison of the flood grid cells in the training dataset with the eleven flood susceptibility maps. The success rate measures how the flood analysis results fit the training dataset. This method divides the area of a flood susceptibility map into five classes, ranging from the highest to the lowest FSI values (Figure 7). The success-rate method uses the training dataset, therefore, it might not be a suitable method for assessing the prediction capacity of the flood susceptibility models.

The prediction rate can explain how well the flood susceptibility models and flood conditioning factors predict flood occurrence. In this study, the prediction-rate results were obtained by comparing the flood grid cells in the validation dataset with the eleven flood susceptibility maps.

The area under the prediction-rate ROC curve was calculated. Figure 9 shows the prediction-rate results of the eleven flood susceptibility maps obtained from IRN-DEMATEL-ANP models.

According to the experimental results, all the flood susceptibility maps had the most acceptable and representable appearance with AUC above 0.95. The WLC fuzzy technique exhibited the best cross-validated performance, followed by the FSI based technique. For instance, the WLC fuzzy technique achieved a very high validation accuracy with an AUC value of 0.988, as shown in Figure 9a,b. Meanwhile, both the visual assessment and quantitative validation through the ROC curve indicated that the FSI based technique was an excellent prediction method with an AUC value of 0.964, as shown in Figure 9c,d.

5. Discussions

The main objective of this work is to investigate the use a novel hybrid MCDA-GIS technique to evaluate FSM, which is constructed by ensemble of DEMATEL, ANP and IRN technique in the case study at Shangyou County, China. In this paper, we improve the DEMATEL method by applying IRN to determine connections in the network structure based on criteria and to accept imprecisions during collective decision making. The application of IRN can eliminate the necessity of additional information to define uncertain number intervals. Therefore, the quality of the existing data during collective decision making and experts’ perceptions that are expressed through an aggregation matrix can be retained.

The considered conditioning factors have not the same contribution to flood occurrence, thus the decision makers should define importance of each conditioning factor by using appropriate weight coefficients. The WLC method requires normalization of the weights. After determining weight coefficients of the conditioning factors, the hybrid IRN-DEMATEL-ANP method is used to calculate the normalized weight conditioning factors. Pair-wise comparisons of the criteria from the IRN-DEMATEL method serve as the input parameters of the ANP method. The weight coefficients of the criteria are obtained as the output values from the ANP method. Finally, WLC method is subsequently used in the GIS environment to obtain the final flood susceptibility map.

According to the importance determination of different used factors for FSM, results of the current study showed that the most important conditioning factor are land cover/use, followed by rainfall, distance from river, altitude and curvature.

Land cover/use is a very important factor in recognizing sensitive regions prone to flooding. Strong correlation of the Land cover/use and flooding is not doubtable as each land use has its own impact on increasing or decreasing the speed of water and therefore generating the flood. Vegetated areas are most often represent protection zones that make the land less prone to flood occurrence [31]. Conversely, urban areas increase surface runoff due to the possession of large impermeable surfaces.

Also the effect of the rainfall was appraised as a main conditioning factor in flood occurrence. The higher rainfall amount usually increases a chance of flood prone area [71]. In this work, the available 1:50,000 rainfall data was considered in the experiments. The main reasons can be summarized as follows: First of all, each factor was converted in the forms of spatially-defined layers of maps with grid size of 30 × 30 m from the ASTER GDEM Version 2 and remote sensing data, along with geological, meteorological and soil data that were provided by different government departments for further analysis. In fact, the rainfall data with the resolution of 1:50,000 was suitable for this study based on our experiments. Second, it is not appropriate for FSM under the assumption that the higher the resolution, the better the result. Since we obtain a pixel-based susceptibility map, using high resolution rainfall data may be useless for practical uses of the map. Moreover, it is not statistically reasonable in most cases if the resolution of other input data cannot accord with the rainfall one. Finally, there is no meteorological data in high resolution and only the 1:50,000 rainfall data are publicly available. However, even if we obtain different factor layers with a very high resolution, the available rainfall data can be interpolated to match the other data using the ArcGIS tool to support FSM.

In addition, one of the conditioning factors, which have high significant impact on flooding, is distance from the river. River levels will increase due to the heavier precipitation during flooding, causing an overflow of water into areas closest to the river bank [24]. The role of curvature in flood occurrence in the study area was also quite remarkable. In the case of flood mostly occurred in the area that has flat curvature. This was confirmed by many different authors [12,24,29,31,34], who also confirmed that rainfall, land cover/use, curvature, elevation, and distance from river had higher impacts on flooding.

On the other hand, in this study, STI, soil type and lithology had lowest importance for flood occurrence. This is partly in agreement with findings reported in Chapi et al. [72] and Shafizadeh-Moghadam et al. [73], even though distace from river was reported as the main effective factor [72,73].

In order to assess the performance of the used IRN-DEMATEL-ANP method, the defined flood susceptibility zones are compared with historical locations of flood events. For evaluation of the method, the area under the curve (AUC) method was used and prediction rate curves were calculated. The capability of IRN-DEMATEL-ANP method was evaluated using a non-dependent threshold approach: the receiver operating characteristic (ROC) curve. Based on overall estimates, the proposed approaches have shown the most acceptable results for mapping the flood occurrence in the study area.

Besides, the proposed modification of the DEMATEL and ANP using RN makes it possible to take into account doubts that occur during the expert evaluation of criteria thus bridging the existing gap in the methodology in the treatment of uncertainty based on RN. In such a way, the quality of the existing data in the collective decision making process can be retained, as well as the experts’ perception, which is expressed through the aggregation matrix. In the surveyed literature, no works on the problem of FSM, taking into account the interdependencies between the criteria in the framework of the rough number-based MCDM approach, have been found. The suggested method was validated through a case study at Shangyou County.

Taking into account the causal relationship between the proposed criteria for FSM, it may be stated that the DEMATEL methodology can support decision making by deriving a visual graph showing the degree of their influence on the final result [74]. In more realistic conditions, the analytical network process (ANP) is also capable of handling interdependencies between the criteria, but it assumes an equal cluster weight to obtain the weighted supermatrix [75,76]. To overcome this shortcoming, DEMATEL-ANP is applied to find the influential weights of the criteria based on the network relationship map [77]. The ANP technique is used to determine the weights of the performance criteria based on the total relation matrix formed by the DEMATEL, thus avoiding pairwise comparisons of the criteria required for ANP analysis [53]. Calculating the relative weights of criteria using traditional ANP means that the levels of interdependence of the factors are treated as reciprocal values. In contrast, in using the DEMATEL method, the levels of interdependence of factors do not have reciprocal values, which is closer to real circumstances [42,43,78,79,80,81]. Because all of the above, implementation DEMATEL method in ANP model gives more objective insight into weight coefficient values.

The results obtained by the proposed method can be considered promising when the reliability and stability of the ranking results are checked. The proposed model can identify the mutual influence of the criteria, thereby aiding purchasing managers to better understand the performance-related issues and to devise the appropriate improvement strategies. Furthermore, it allows for applying the developed framework for case studies from the perspective of particular areas and comparing the results by examining possible differences and their causes.

6. Conclusions

This work validates the application of a novel MCDA-GIS framework for FSM in the case of Shangyou County, China. To accommodate uncertainties in the MCDA process, the IRN technique is used to handle this problem. The application of IRN in the MADA process is presented through the hybrid IRN-DEMTEL-ANP method. The final susceptibility map was obtained using the WLC method. Some conclusions can be drawn as follows. First, the proposed method can conduct FSM by combining both qualitative and quantitative methods. It can make full use of experts’ knowledge, which is different from statistical learning methods. Second, the validation process was performed based on the comparison of the historical flood locations to the different flood susceptible zones on the final map. The higher AUC values demonstrate the effectiveness of the proposed method. Specifically, according to the resultant FSI scales, 71.6 km² (4.6%), 329.8 km² (21.4%) and 588.4 km² (38.1%) of the study area were labelled with ‘very high’, ‘high’ and ‘moderate’, respectively, which proves that Shangyou County is mostly located in the flood-prone areas. These high susceptible zones locate at the eastern and southern areas which are near river flows with lower altitudes and lowland reliefs. In the future, we will take some useful ideas on the relationship between planning and flood assessment using the IRN-DEMTEL-ANP method, following the contributions of [82,83]. Moreover, the future improvement of the proposed method may be based on classical fuzzy sets and the intervals of fuzzy numbers when determining the parameters of the conditioning factors. Besides, the implementation of individual phases of the Best-Worst method can be involved in the phases of the DEMATEL or ANP methods.