Flood Risk Assessment Based on Fuzzy Synthetic Evaluation Method in the Beijing-Tianjin-Hebei Metropolitan Area, China

Abstract

Flooding is one of the most devastating natural events and leads to enormous and recurring loss of life, properties, and resources around the globe. With climate change and accelerating urbanization, flood disasters in China have increasingly affected the sustainable development of metropolitan areas. Risk assessment is an essential step in flood management and disaster mitigation, which provide a quantitative measure of flood risk. However, the difficulty of flood risk zoning is dealing with the uncertainty of the evaluation process and the complicated non-linear relationship between indicators and risk levels. To address this issue, a fuzzy synthetic evaluation (FSE) method based on combined weight (CW) was utilized in this paper to generate flood risk maps at a grid-scale (1 × 1 km). For the case study in the Beijing-Tianjin-Hebei metropolitan area (BTH) in China, fourteen indicators were selected to construct the flood risk assessment model based on the FSE approach integrated with ArcGIS. The research demonstrates that moderate, high, and very high risk zones are distributed in the southeast fluvial plain of the BTH area, accounting for 31.36% of the total land area. Meanwhile, low and very-low risk zones occupy 68.64% of the total land area, and are primarily located in the high plateau and mountain regions in the northwest. We analyzed the risk level of each county and proposed risk mitigation measures based on field investigations. The verified risk assessment results were spatially consistent with the historical flood disaster records and loss positions, indicating the accuracy and reliability of the risk assessment map using the FSE approach. Compared with the IPCC (Intergovernmental Panel on Climate Change) TAR (Third Assessment Report) and AR5 (Fifth Assessment Report) methods, FSE has significant advantages in handling uncertainty, complexity, and the non-linear relationship between indices and risk grades. This study provides a novel quantitative method for flood risk assessment in metropolitan areas and practical implications for urban flood management.

1. Introduction

Flooding is generally considered to be one of the most common and destructive natural disasters, and has been responsible for major losses to human and social economies [1,2]. A preliminary survey by the International Disaster Database (EM-DAT) shows that from 1900 to 2016, floods have killed approximately 7 million people and caused more than $700 billion in damage worldwide [3]. In China, metropolitan areas have become the strategic cores and future development orientations of new urbanization, and they play an essential role in globalization and national economic development. However, acceleration of urbanization and intense human activities has led to dramatically increased flood risks in the past few decades [4,5,6]. Worse still, most of the flooded areas in China have been primarily located in metropolitan areas and urban agglomerations with large populations and advanced economies [7]. Previous studies have focused on a global [8], national [9], watershed [10], single-city [11] and community scales [12], but few studies have been conducted at the scale of metropolitan area [13]. Floods have become a major threat to the sustainable development of metropolitan regions in China [14]. Within this context, identifying high risk areas at the regional level and adopting substantive optimum flood management strategies is of considerable significance to metropolitan areas.

Flood risk assessment is an essential means of flood management and risk mitigation that demonstrates significant practical values in floodplain management, disaster warning and evacuation, flood insurance and decision-making [15,16]. In recent decades, numerous methods have been proposed and applied to flood risk assessment. Traditionally, these approaches can be classified into four primary types: GIS (Geographic Information System) and remote sensing methods (GRSM) [17], scenario simulation analysis (SSA) [18], the historical disaster mathematical statistics method (HDMS) [19], and multi-criteria analysis (MCA) [20,21]. GRMS requires high-resolution remote sensing images, and thus has not been widely used due to higher costs and difficulties in data acquisition [22]. SSA uses hydrological/hydraulic models for flood routing simulation, although the impacts of hazard-formative environments and hazard-affected bodies are generally not considered in risk evaluation results (except for high-resolution topographic and meteorological data). HDMS uses a substantial amount of historical disaster data, such as affected population, collapsed buildings and other socio-economic information, but does not accurately reflect the spatiotemporal variability of flood risk. Additionally, due to discrepancies in risk perception among experts, different risk indicator division criteria will generate several risk maps, which causes the MCA method to have high uncertainty.

Currently, the main challenge facing flood risk assessment is the fuzziness and uncertainty of risk estimation and the complicated relationships between risk indicators and risk levels [1,23]. The aforementioned traditional risk assessment methods cannot handle the uncertainty and ambiguity efficiently of the evaluation process. To solve this problem, we proposed a fuzzy synthetic evaluation (FSE) method based on combined weight (CW). FSE is a method that uses fuzzy sets and fuzzy logic theory to make an overall evaluation of complicated things affected by multiple factors [23,24]. Therefore, FSE can efficiently tackle fuzziness or uncertainty in risk evaluation [23]. FSE has been widely used in various research fields, but it is rarely utilized in flood risk assessment research. In this study, we use the FSE method to perform flood risk assessments in the Beijing-Tianjin-Hebei metropolitan area (BTH) to construct an accurate and reliable risk map.

Furthermore, since different indicators exhibit different contributions to flood risk, the assignment of indicator weights is crucial to the risk assessment results [15,25]. The two current primary weighting methods are subjective weighting and objective weighting. The analytic hierarchy process (AHP) and entropy weight (EW) method are typical representatives of subjective and objective weighting methods, respectively [26,27]. Nevertheless, AHP is significantly affected by the perception of experts, and cannot adequately reflect the multiple data provided by risk indicators [15]. On the contrary, EW excessively depends on the objective data of risk indicators without considering the relative importance of each indicator to the risk outcome [15]. Therefore, to avoid the limitations of the single weight method, a combined weight (CW) method with both advantages was proposed and applied in this research.

The main objectives of this study are to develop a flood risk assessment model for metropolitan areas based on the FSE approach, to verify the feasibility of the results using historical flood data and casualty data, to put forward a combined weight (CW) method based on subjective (AHP) and objective (EW) weighting methods to improve the scientificity of weight distribution, and finally to generate risk maps of the BTH area from a grid (1 × 1 km) and propose corresponding mitigation measures combined with field investigations at the regional level.

2. Study Area and Data

2.1. Study Area

The BTH area is located in northern China and comprises the two central cities and eight prefecture cities of Hebei Province, including Beijing, Tianjin, Chengde, Zhangjiakou, Qinhuangdao, Tangshan, Baoding, Langfang, Cangzhou, and Shijiazhuang (Figure 1). This region consists of three geomorphic units, including the Bashang Plateau, the Yanshan Mountains and Taihang Mountains, and the vast alluvial plain (part of the North China Plain). The BTH area has the typical characteristics of a warm temperate semi-humid and semi-arid continental monsoon climate [28]. The mean annual precipitation is about 500 mm, of which about 80% occurs from June to August [28]. This region is located in the Haihe River and Luanhe River basins, with a dense river network and artificial channels.

The BTH area is the political, economic, and cultural center of China, and has experienced a rapid urbanization process since the 1980s [29,30]. The population of the area is 91.2 million (2017), accounting for 7% of China’s total population. The regional GDP is $127.19 billion, consisting of 9% of the national total GDP (2017). Urban expansion, intense human activity, and rainstorm weather have commonly contributed to severe flooding in the area. On 21 July 2012, a 60-year return period of extreme rainstorms hit the BTH area, with a maximum rainfall of 541 mm in Hebei Town, Fangshan District, Beijing. This extreme rainstorms resulted in wide waterlogging in urban areas and torrential flooding in mountain regions, with 115 human casualties and a direct economic loss of $4.88 billion [31,32]. To cope with the flood and reduce potential loss, a systematic flood risk assessment of the BTH area is of great practical significance.

2.2. Indicator Selection and Data Analysis

Regional risk indicator variables differ due to differences in physical geography and socio-economic characteristics [29]. An indicator that shows a high degree of impact on flood risk in a specific area may not rank similarly in other regions [33]. According to the conceptual model of risk assessment and the actual conditions in the BTH area [34,35,36,37], 14 risk indicators were determined. Detailed descriptions and data sources for each indicator are as follows:

• Digital elevation model data (DEM, m): DEM can reflect the absolute elevation and topographic fluctuations of the study area. In general, flooding is more likely to occur at lower altitudes than at higher elevations. The original resolution of DEM is 90 × 90 m.
• Slope data (SL, degree): The slope reflects the topographic change and runoff velocity. In general, steep slopes of the mountain areas promote runoff generation, whereas low-lying level regions are prone to severe flooding and waterlogging due to poor drainage. SL was extracted from the DEM data using GIS technique.
• Rainfall intensity data (RI, normalized value): By analyzing the daily precipitation data of 126 precipitation stations in the study area from 1981 to 2010, the maximum three-day precipitation was calculated to characterize rainfall intensity [38]. Then, spatial interpolation was performed using Kriging interpolation and the RI was generated.
• Rainstorm frequency (RF, times/year): Rainstorm refers to a heavy rainfall event with precipitation of more than 50 mm within 24 hours in China. The rainfall frequency of each rain gauge station was calculated based on the daily precipitation data from 1981 to 2010.
• Drainage density (DD, km/3 km²): Based on the linear density analysis technology of ArcGIS, the drainage density was analyzed using the river vector data in the study area [16]. The search radius was 3 km and a data layer of 1 × 1 km was generated.
• Normalized difference vegetation index (NDVI): This indicator can effectively reflect the distribution of vegetation coverage on the spatial scale of the study area. The vegetation not only regulates the surface runoff but also prevents soil erosion [39].
• Soil texture (ST): Different soil textures have various infiltration capacities and a specific impact on surface runoff. In this study, we used the soil infiltration capacity according to the code values in the Harmonized World Soil Database (HWSD). The soil texture classification criteria in this database were divided into 13 categories by the US Department of Agriculture (USDA) [1]. As shown in

  Table 1. Soil texture (ST) and its infiltration capacity (IC).

  ST	Sand	Loamy Sand	Sandy Loam	Sandy Clay Loam	Loam	Sand Clay	Silt Loam
  IC	13	12	11	10	9	8	7
  ST	Silt	Clay loam	Silty clay loam	Clay	Silty clay	Clay (heavy)
  IC	6	5	4	3	2	1

  , a large code value indicates a strong infiltration capacity.
• Gross domestic product (GDP, yuan/km²): GDP indicates the prosperity degree of the national economy of a country or region, which is an essential indicator of flood vulnerability. In this paper, we collected the GDP data from 2015 with a resolution of 1 × 1 km.
• Population density (PD, people/km²): This indicator reflects the spatial distribution of human exposure in the research area in 2015.
• Grain output (GO, kg/hm²): Food is a special commodity and an essential strategic reserve related to people’s livelihoods. This data mainly considers the grain yield of five crops (wheat, corn, rice, soybean, and potato) in the study area.
• Road network density (RND, km/km²): The strong hydrodynamic conditions during the flood will inevitably destroy roads. Furthermore, secondary geological disasters caused by floods, such as landslides and debris flows in mountainous areas, pose severe threats to road and rail network systems.
• Per capita disposable income (PDI, 10,000 yuan/people): PDI can reflect the overall socio-economic level of a region. The higher the PDI, the stronger the disaster resistance and the faster the post-disaster recovery and reconstruction.
• The ratio of children and the elderly (RCE, %): Children (0 to 14 years old) and older people (65 years and older) are more vulnerable to flooding due to age and some physiological reasons.
• Average schooling years (ASY, year): ASY is an important indicator to measure the degree of popularization of higher education in a country or region (over six years old). Generally, ASY is positively related to people’s awareness of disaster prevention and reduction.

Out of the 14 indexes, there are seven hazard indicators (DEM, SL, RI, RF, DD, NDVI, and ST) and seven vulnerability indicators (GDP, PD, GO, RND, PDI, REC, and ASY), respectively. The 14 indicators were converted into a grid data layer of 1 × 1 km using the GIS technique (Figure 2), with a total of 180,816 grids. The data sources and original resolutions of the 14 indices are shown in

Table 2. Data sources and original resolutions of 14 risk indicators.

Indices	Time	Original Resolution	Data Sources
DEM	/	90 × 90 m	Computer Network Information Center of the Chinese Academy of Sciences (http://www.gscloud.cn/)
SL	/	90 × 90 m	Same as above
RI	1981–2010	Vector data	National Meteorological Information Center (http://data.cma.cn/)
RF	1981–2010	Vector data	Same as above
DD	2017	Vector data	National Geomatics Center of China (http://www.ngcc.cn/ngcc/)
RND	2017	Vector data	Same as above
NDVI	2018	1 × 1 km	Resource and Environment Data Cloud Platform of the Chinese Academy of Sciences (http://www.resdc.cn/)
GDP	2015	1 × 1 km	Same as above
GO	/	1 × 1 km	Same as above
ST	/	1 × 1 km	Food and Agriculture Organization of the United Nations (http://www.fao.org/home/en/)
PD	2015	1 × 1 km	Peking University Geographic Data Platform (https://geodata.pku.edu.cn/)
PDI	2017	County unit	Beijing Regional Statistical Yearbook (2018), Tianjin Statistical Yearbook (2018), Hebei Economic Yearbook (2018)
RCE	2017	County unit	Same as above
ASY	2000	County unit	Fifth population census data of the National Bureau of Statistics of the People’s Republic of China (2000)

Notes: digital elevation model data (DEM); normalized difference vegetation index (NDVI); slope data (SL); rainfall intensity data (RI); rainstorm frequency (RF); drainage density (DD,); soil texture (ST); gross domestic product (GDP); population density (PD); grain output (GO); road network density (RND); per capita disposable income (PDI); ratio of children and the elderly (RCE); average schooling years (ASY).

.

3. Methodology

3.1. Overall Framework of the Evaluation Procedure

We adopted a fuzzy synthetic evaluation (FSE) to eliminate the possible fuzziness or uncertainty that may exist in flood risk assessment in the BTH area. The structure and process of this paper consist of four parts: data collection and preprocessing, indicator selection and analysis, fuzzy synthetic evaluation, and verification analysis. Figure 3 shows the detailed procedure of flood risk assessment in the BTH area based on the FSE method.

3.2. Fuzzy Synthetic Evaluation (FSE)

Fuzzy Synthetic Evaluation (FSE) is a kind of method that uses fuzzy mathematic theory to make an overall evaluation of complicated things that are affected by multiple factors [23]. The method has been applied to flood risk assessment to eliminate the fuzziness or uncertainty that may exist in the assessment process [40]. Furthermore, FSE converts qualitative evaluation into quantitative evaluation by utilizing fuzzy set theory and fuzzy logic, thereby allowing it to obtain precise evaluation results [23,24,25,26,27,28,29,30,31,32,33,34,35]. Three significant analytical processes are indispensable to implementing the FSE approach: fuzzification, fuzzy inference, and defuzzification, respectively. Fuzzification is the process of converting the determined values of the input quantity of a fuzzy controller into corresponding fuzzy linguistic variables. Fuzzy inference process combines membership functions (MFs) and control rules to output a fuzzy relation matrix. Defuzzification is the process that maps fuzzy sets to corresponding crisp sets based on different fuzzy operators in order to obtain final results. In this paper, the above three processes correspond to fuzzy classification, membership functions, and comprehensive evaluation, respectively.

3.2.1. Fuzzy Classification

During the application of FSE in risk assessment, flood risk was divided into five grades according to actual requirements, namely, very low, low, moderate, high, and very high, respectively. Discontinuous points (d_i) were then obtained, using the statistics of standard deviation or mean value from the original values of each indicator [23,24].

The interval values of several indicators, such as SL and DEM, were determined based on the actual situation of the research area. Discontinuities of the SL refer to the five-grades division method proposed by the Agricultural Division Committee of China in the Technical Regulation of Land Use Status Survey of 1984. Considering that the slopes of the coastal plain and the lacustrine plain in the east of the research area are almost zero, a 0.2-degree discontinuity was inserted to identify these highly hazardous regions. Thus, discontinuities of the SL were classified as 0.2°, 2°, 6°, 10°, and 20°. Moreover, interval values of DEM were ultimately divided into 5, 50, 500, 1000, and 1500 mm, which were based on an overall analysis of the topography and geomorphology of the study area and combined with geomorphic classification standards of China.

The critical values of GO and ST were generated by the method of natural breaks (Jenks), which is based on ArcGIS. Jenks is a clustering method that sets boundary values at a position where the difference is relatively evident according to the inherent natural group features of the data. It aims to reduce the variance in each category and maximize the variance between different classes. The discontinuous points (d_i) of all indicators in the BTH area are shown in

Table 3. Fuzzy classification of risk assessment indicators in the BTH area.

Flood Disaster Risk	Indicators	Fuzzy Classification (Discontinuities)					Indicator Attributes
		d1	d2	d3	d4	d5
Hazard	Digital elevation model data (DEM, m)	1500.00	1000.00	500.00	50.00	5.00	negative
	Slope data (SL, degree)	20.00	10.00	6.00	2.00	0.20	negative
	Rainfall intensity data (RI, Normalized value)	0.21	0.37	0.53	0.69	0.85	positive
	Rainstorm frequency (RF, times/year)	0.50	1.50	2.50	3.50	4.50	positive
	Drainage density (DD, km/3 km2)	0.00	0.12	0.24	0.36	0.48	positive
	Normalized difference vegetation index (NDVI)	0.90	0.81	0.72	0.63	0.54	negative
	Soil texture (ST)	13.00	12.00	11.00	9.00	7.00	negative
Vulnerability	Gross domestic product (GDP, yuan/km2)	50.00	200.00	1000.00	3000.00	5000.00	positive
	Population density (PD, people/km2)	1.00	200.00	400.00	800.00	1600.00	positive
	Road network density (RND, km/km2)	0.00	0.55	1.10	1.65	2.20	positive
	Grain output (GO, kg/hm2)	967.94	2581.17	4463.27	6668.02	8872.77	positive
	Per capita disposable income (PDI, 10000 yuan/people)	2.64	2.29	2.01	1.76	1.55	negative
	Ratio of children and the elderly (RCE, %)	23.00	25.00	26.00	27.00	28.00	negative
	Average schooling years (ASY, year)	9.40	8.80	8.20	7.60	7.20	negative

The indicator attributes refer to the relationships between the indicators and the flood risk degrees. Positive indicators indicate a higher flood risk degree, and negative indicators indicate the opposite.

.

3.2.2. Membership Function and Comprehensive Evaluation

A fuzzy set quantizes the fuzziness through membership functions (MFs), eliminating the uncertainty and obtaining a fuzzy evaluation matrix [41,42]. Therefore, determining suitable membership functions is critical to risk assessment results. There are many types of membership functions, including triangular waveforms, trapezoidal waveforms, Gaussian waveforms, bell-shaped waveforms, sigmoidal waveforms, S-curve waveforms, and so on. According to previous studies, if a system has a significant dynamic variation in a short period, it is more appropriate to utilize a triangular or trapezoidal waveform [23]. Therefore, considering that floods are short-duration events, it is ultimately necessary to select the trapezoidal and triangular waveforms in order to define the piecewise functions.

As shown in Figure 4, u_ij is the membership degree of the indicator i at the grade j. It is not difficult to find that the actual data x_i of an indicator that has two membership degrees (u_i2 and u_i3). However, it can be derived that x_i has only one membership degree according to the principle of maximum membership degree in the fuzzy operator. Thus, FSE converts fuzziness and uncertainty into accuracy and certainty through membership functions based on fuzzy set and fuzzy logic theory. As shown in Appendix A, Equations (A1)–(A3) and (A4)–(A6) are used to calculate the fuzzy membership of positive and negative indicators, respectively.

3.2.3. Comprehensive Evaluation

Once the membership functions are determined, a membership function matrix (R_k) is generated according to Equations (A1)–(A6) that combines the actual values of all indicators. The resulting vector B_k is the product of the weight vector W and the fuzzy relation matrix R_k, and its mathematical derivation is as follows:

$$B_k=W\times R_k=\begin{array}{c}[\end{array}{w}_1,w_2,\dots ,w_i]\times \left[\begin{array}{cccc}{u}_{11}& u_{12}& \dots & u_{15}\\ u_{21}& u_{22}& \dots & u_{25}\\ \dots & \dots & \ddots & \dots \\ u_{i1}& u_{i2}& \dots & u_{i5}\end{array}\right]=\begin{array}{c}[\end{array}{b}_1,b_2,b_3,b_4,b_5]$$

(1)

where W is the weight vector of indicators (its calculation process is described in detail below). The element of B_k is the comprehensive membership degree, which indicates the extent to which each evaluation unit (or grid) belongs to the five risk assessment levels. In this paper, risk assessment grade is ultimately determined by the maximum membership. Therefore, the comprehensive flood risk level (CFRL) of each evaluation unit can be obtained by inputting Equation (2) in the GIS environment.

$$\begin{array}{l}CFRL=Con(\max (b_1,b_2,b_3,b_4,b_5)=b_1,1,\\ Con(\max (b_1,b_2,b_3,b_4,b_5)=b_2,2,\\ Con(\max (b_1,b_2,b_3,b_4,b_5)=b_3,3,\\ Con(\max (b_1,b_2,b_3,b_4,b_5)=b_4,4,5))))\end{array}$$

(2)

3.3. Calculation of the Weights for Indices Based on the CW Method

The commonly used weight distribution methods mainly include subjective weight (entropy weight, EW) and objective weight (analytic hierarchy process, AHP). Although the EW method can genuinely reflect the useful information offered by each evaluation index, it relies excessively on actual data. On the contrary, the AHP method is proven to determine weights comprehensively by considering the subjective attributes of the data, but it may not accurately reflect the useful information provided by the data [15]. In order to address the contradiction between objective and subjective weight and make full use of the advantages of both, a combined weight (CW) method combining EW and AHP was proposed and applied in our flood risk assessment research. The theoretical knowledge of the entropy is beyond the scope of this paper, but a comprehensive and detailed introduction can be found in a previous work of Shannon [43]. The main calculation steps of the CW are summarized in Appendix B.

4. Results and Discussions

4.1. Weight Calculation Results

Using the weight calculation steps listed in Section 3.3., we used the CW method to calculate the weights of 14 risk indicators. Initially we constructed judgment matrixes for each criterion layer by using the pairwise comparison method, which we based on importance rankings of the risk indicators by 10 experts in the field of risk assessment (

Table 4. Judgment matrixes in the analytic hierarchy process (AHP) method.

Indicator			Hazard				Vulnerability
Hazard			1				4/3
Vulnerability			3/4				1
λmax = 2; CI = 0; RI = 0; CR = 0
Indicator	RI	RF		DEM	DD	SL		NDVI	ST
RI	1	1		2	2	3		4	5
RF	1	1		2	2	3		4	5
DEM	1/2	1/2		1	2	3		3	4
DD	1/2	1/2		1/2	1	2		3	4
SL	1/3	1/3		1/3	1/2	1		2	3
NDVI	1/4	1/4		1/3	1/3	1/2		1	3
ST	1/5	1/5		1/4	1/4	1/3		1/3	1
λmax = 7.215; CI = 0.036; RI = 1.320; CR = 0.027
Indicator	GDP	PD		GO	RND	PDI		RCE	ASY
GDP	1	1		3	4	5		7	8
PD	1	1		3	4	5		7	8
GO	1/3	1/3		1	3	4		5	7
RND	1/4	1/4		1/3	1	3		5	5
PDI	1/5	1/5		1/4	1/3	1		4	5
RCE	1/7	1/7		1/5	1/5	1/4		1	3
ASY	1/8	1/8		1/7	1/5	1/5		1/3	1
λmax = 7.574; CI = 0.096; RI = 1.320; CR = 0.073

). The largest eigenvalue (λ_max) of the judgment matrix and its corresponding maximal eigenvector could then be calculated by the weighted sums approach. The judgment matrix was considered to be satisfactory only when the consistency ratio (CR) was less than 0.1 [44]. After verification, all CRs in this paper were less than 0.1 (Table 4), which fully demonstrates that the weight calculation results based on AHP were reasonable and reliable. The calculation result of the AHP weight is shown in Table 4. Subsequently, we used the Model Builder in ArcGIS to construct the calculation procedure of the CW method, and the weight of the AHP and EW of the 14 risk indicators were calculated, respectively. Based on this, the combined weights were finally generated (

Table 5. Weight outputs based on different methods.

Indicator	AHP (Analytic Hierarchy Process)	EW (Entropy Weight)	CW (Combined Weight)
DEM	0.0999	0.0062	0.0058
SL	0.0485	0.0024	0.0011
RI	0.1450	0.0224	0.0303
RF	0.1450	0.0192	0.0260
DD	0.0777	0.0224	0.0162
NDVI	0.0343	0.0496	0.0159
ST	0.0206	0.0504	0.0097
GDP	0.1313	0.3133	0.3843
PD	0.1313	0.3414	0.4188
RND	0.0438	0.0651	0.0266
GO	0.0712	0.0954	0.0635
PDI	0.0279	0.0018	0.0005
RCE	0.0142	0.0078	0.0010
ASY	0.0094	0.0028	0.0002

).

Table 5 indicates the inconsistencies between AHP and EW weights. DEM and SL were considered essential factors influencing flood risk in AHP weightings and thereby are assigned larger weighting values. However, the weights of both were relatively small in the EW method; this is mainly because the southeast part of the study area is a vast alluvial plain (part of the North China Plain), while the west and north are mountainous areas. In terms of spatial heterogeneity, compared with other risk indicators (such as GDP, RI, and RF), DEM’s and SL’s variability are relatively moderate. According to the basic principle of the EW method, since DEM and SL provide less useful information on risk assessment, the weights of both are inevitably relatively low. On the contrary, GDP and PD have larger weights in both methods, which is because the social economy and population are the primary hazard-bearing bodies, while on the other hand the spatial distribution of GDP and PD in the study area is hugely uneven. Therefore, it is scientifically reasonable for DEM and SL to be assigned as larger weights.

Furthermore, this paper proposes a combined weighting method based on the AHP and EW methods that not only considers the contributions and relative significance of different risk indicators to the final risk level, but also makes full use of the useful information offered by the risk indicators. The combination weighting method is put forward to take advantage of subjective and objective weights and to make the weighting results more reasonable.

4.2. Distribution Characteristics of Flood Risk in the BTH Area

Based on the calculation approach mentioned above, the spatial distribution map of flood hazard, flood vulnerability, and flood risk in the BTH area were generated at a grid scale (1 × 1 km). As shown in Figure 5c, the comprehensive flood risk in the research area is roughly centered on the Beijing area, demonstrating a “petal-like” distribution pattern. Very high risk zones are mainly located in the central, eastern and southwestern regions of the study area, including most of the counties or districts in Beijing, Tianjin, and Tangshan and accounting for approximately 22.23% (24,871 km²) of the total area. Furthermore, very high risk areas also lie in river valleys and the piedmont alluvial plain areas of the Taihang and Yanshan mountainous regions, most of which are scattered in dots and bands. High risk zones alternate with the very high risk zones distributed predominantly in the eastern and southern parts of the study area covering approximately 7.25% (13,104 km²) of the total area. Specifically, these areas include some districts and counties in Shijiazhuang, Tangshan, and Qinhuangdao. Moreover, the research results clearly show that a total of 29.48% of the regions are under very high or high flood risk, and that these are regions with dense populations and developed social economies.

The moderate risk zones mainly lie in the transition zone between mountains and plains, accounting for about 9.50% (17,179 km²) of the total area. The central and western regions of Shijiazhuang, the central and eastern regions of Baoding, and the southern areas of Qinhuangdao are all moderate risk zones. In addition, the very low and low risk zones are spread over vast areas of the northern and western parts of the study area, including Zhangjiakou, Chengde, western parts of Beijing, Baoding, Shijiazhuang, and most of Cangzhou. These areas are mainly located in the Taihang and Yanshan mountainous areas, which have sparse populations and relatively backward economies. The surface areas of the very low and low risk zones are 111,886 and 13,776 km², occupying 61.02% and 7.62% of the total area, respectively. Moreover, other risk levels are distributed sporadically across a wide range of very low and low risk areas.

According to disaster risk theory, the flood risk in an area is the result of the interaction of physical geographical elements and socio-economic factors [45]. Therefore, spatial distribution maps of flood risk that consider hazard or vulnerability indices alone are neither integrated nor objective. In order to compare the differences between these indices, we generated a flood hazard map (Figure 5a) and a flood vulnerability map (Figure 5b), respectively, based on the FSE approach. Although the spatial distribution of hazards and the comprehensive flood risk map had similar trends, the differences were noticeable in certain regions that are significantly affected by the socio-economic indicators. There are areas with low hazard intensity that may also have high flood risk. For example, most areas in Baoding and Shijiazhuang in the southwest of the study area are the high and very high risk zones, although their flood hazards are moderate or even low. This is due to their dense populations, transportation, property, and crops, and small floods would cause severe casualties and flood losses.

In contrast, hazards that strike in regions with low exposure and vulnerability will not become disasters, or at least there will be no catastrophic flood disasters. The primary reason is that these areas are sparsely populated and economically backward (such as Qinglong Manchu Autonomous County in the north of Qinhuangdao), and although the probability of being inundated is relatively high, there are few casualties and property losses. Therefore, comprehensive flood risk zoning should consider not only hazardous factors such as precipitation, topography, rivers, and so on, but also vulnerability factors such as population, property, crops, and so on. Such a risk zoning map could be more representative of the actual flood risk of the study area.

To provide a more specific spatial distribution of flood risk from administrative units and to identify the susceptibility of each county or district towards risk, the flood risk of each county or district was classified into five categories based on the GIS technique. As shown in Figure 6 and Figure 7, we can better identify higher risk counties or districts using this method, and we can determine which regions should pay significant attention and urgently improve their capabilities of disaster prevention and mitigation. Risk information at the county level provides scientific guidance for the government and non-governmental organizations to prioritize the allocation of flood control resources. Moreover, it offers accurate risk backgrounds for urban planning, evacuation drills, and insurance companies.

According to the estimation results, the risk level of the southeast region is higher than that of the northwest region in general. High and very high risk zones are mainly distributed in densely populated and economically developed plain areas, while very low and low risk zones are mainly located in sparsely populated and economically backward mountainous regions. Moderate risk zones primarily lie in the transition zone between mountains and plains due to their moderate populations and relatively developed economies.

4.3. Results Validation

In order to verify the overall flood risk assessment results, 126 historical flood records and the submerged areas of two historical floods in the Haihe River Basin were applied, most of which were provided by water conservancy departments. A few typical flood data were obtained through field investigations. During field surveys, high-precision GPS (PromarkTM 500) was used to record point data of historical floods to ensure accuracy. The historical floods were converted into vector format and the inundation areas were obtained by using the “Extract by Mask” tool in ArcGIS (Figure 8). Subsequently, the accuracy of the risk assessment results was verified through a comparative analysis with historical flood data.

The black diamond symbols in Figure 8a show the locations where 126 historical floods occurred, including both measured and field investigations. The blue triangle symbol in Figure 8a represents the locations of casualties in historical floods in data released by the governments of Beijing and Hebei Province from a death toll exceeding 110. Through overlay analysis, we precisely calculated the risk level of the grid cells where the historical flood coordinate points were located. Based on the above approach, we found that 88.1% of the historical flood records and 95.5% of the locations of casualties were distributed in the moderate to very high risk zones, and that less than 4% of the historical flood records and only 4.5% of casualties were located in the very low and low risk areas (

Table 6. Results verification using historical flood records.

Risk Levels	Very High	High	Moderate	Low	Very Low
Historical floods (%)	34.13	15.87	38.10	8.73	3.17

and

Table 7. Death location of a historical flood and its corresponding risk levels.

Risk Levels	Very High	High	Moderate	Low	Very Low
Dead points (%)	77.50	11.70	6.30	2.70	1.80

). This consistency suggests that the flood risk assessment results using the FSE method were convincing.

Due to the difficulty in obtaining historical flood data, the historical flood hydrological records collected in this paper mainly detail the Taihang and Yanshan mountains. The reliability of the risk assessment results in mountainous areas was proved by superposition analysis. Subsequently, we used the submerged area of two historical floods in the Haihe River Basin to further verify the risk assessment results. Figure 8b,c display the submerged area (blue lines) of the historic floods of 1939 and 1963, respectively. The flood in 1939 resulted in most of the southeastern part of the study area being submerged, covering approximately 31,713 km² and accounting for 17.54% of the total area. The submerged area of the 1963 flood was mainly located in the south-central part of the study area, with an area of 24,342 km² occupying 13.46% of the total area. The results of layer overlay analysis are shown in Table 6 and Table 7. For the 1939 flood, 66% of the inundation area was consistent with moderate to very high risk zones. The results of the superimposed analysis of the 1963 flood are generally in accordance with the 1939 flood, and 65.80% of the regions overlap with moderate to very high risk zones (

Table 8. Results validation using the inundation areas of the 1939 and 1963 floods.

Historical Floods	Inundated Area (km2)	Percentage of Total Area (%)	Overlapping Ratio with Different Risk Levels (%)
			Very High	High	Moderate	Low	Very Low
1939	31713	17.54	34.90	15.50	15.60	6.40	27.70
1963	24342	13.46	22.30	17.20	26.30	7.30	27.00

). Compared with the 1939 flood, high and very high risk zones in the 1963 flood were 12.60% more than in the 1939 flood, and moderate risk areas were 10.70% less than the 1939 flood (Table 8). Additionally, the low and very low risk areas of the two historical floods were both about 27% (Table 8). The principal reason is that these areas (such as the coastal area of the research area and parts of Cangzhou) are sparsely populated and economically backward, and although the probability of being inundated is relatively high, few casualties and property losses would occur here. Spatially, the inundation area data support the accuracy of the flood risk assessment result.

Consistent matching of the historical flood records with the overall flood risk results, as demonstrated through the above analysis, powerfully illustrates the reliability and accuracy of the FSE method proposed in this paper.

4.4. Comparison with other Methods

There are multiple conceptual models and methodologies related to flood risk assessment, although risk assessment is a global concept [16]. In order to further analyze and demonstrate the advantages of the proposed method in this paper, the IPCC TAR (Third Assessment Report) and AR5 (Fifth Assessment Report) methods were compared to the FSE method.

In the IPCC TAR approach, risk is assessed using a multiplicative formula of hazard and its probability of occurrence along with its vulnerability [46,47]. Furthermore, vulnerability is evaluated from socio-economic indicators and in a linear form of exposure, sensitivity, and adaptive capacity. The formula according to TAR [46,47] is described as

Risk = Hazard × Vulnerability

(3)

where

Vulnerability = Exposure + Sensitivity − Adaptive Capacity

(4)

In the IPCC AR5 approach, risk is defined as the results of the interaction of hazards with vulnerability and exposure of human and natural systems [48], and its conceptual model expression is Equation (5).

Risk = Hazard × Exposure × Vulnerability

(5)

IPCC AR5 gives a detailed definition of related notions of the risk involved in the above formula [48]. According to its description, vulnerability encompasses a variety of concepts and elements including sensitivity or susceptibility to harm, and lack of capacity to cope and adapt [48]. The vulnerability is described as Equation (6).

Vulnerability = Sensitivity − Adaptive Capacity

(6)

Therefore, the calculation formula of two IPCC methods are as follows:

$$Ris{k}_{TAR}=[{\displaystyle \sum (H\ast W_H)}]\times [{\displaystyle \sum (E\ast W_E)}+{\displaystyle \sum (S\ast W_S)}-{\displaystyle \sum (AC\ast W_{AC})]}$$

(7)

$$Ris{k}_{AR5}=[{\displaystyle \sum (H\ast W_H)}]\times [{\displaystyle \sum (E\ast W_E)]}\times [{\displaystyle \sum (S\ast W_S)}-{\displaystyle \sum (AC\ast W_{AC})]}$$

(8)

where Risk_TAR and Risk_AR5 are the flood risk of the TAR and AR5 approaches, respectively; the values of the H, E, S, and AC represent the indicators of hazard, vulnerability, sensitivity, and adaptive capacity, respectively; W_H, W_E, W_S, and W_AC indicate the weights of H, E, S, and AC, respectively.

To implement the two approaches, we reclassified flood risk indicators based on the judgments of experts and the definitions of terminology for conceptual models proposed by IPCC [46,47,48] (

Table 9. Distribution of indicators among different domains in IPCC methods.

Flood Risk	Scenario 1	Scenario 2
	Index Distribution in TAR and AR5	Index Distribution in AR5
Hazard	RI, RF	RI, RF
Exposure	GDP, PD, GO, RND	GDP, PD, RND
Sensitivity	DEM, SL, DD, NDVI, ST	GO, DEM, SL, DD, NDVI, ST
Adaptive Capacity	PDI, RCE, ASY	PDI, RCE, ASY

). Flood risk is determined by indicators in four domains (hazard, exposure and adaptive capacity) and corresponding weights according to the IPCC calculation formula. However, the IPCC AR5 method does not restrict any fixed boundary for assigning a specific indicator to a particular domain. In this study, two scenarios were created by applying the AR5 in respect to the opinions of four experts. In Scenario 2 of the AR5 (Table 9), the indicator grain output (GO, the red font in Table 9) is moved from the exposure domain to the sensitivity domain, according to the perception of two other experts.

The ‘Raster Calculator’ under the ‘Spatial Analyst’ module in ArcGIS was utilized to normalize the risk indicator layer by using Equations (A8) and (A9). The same weighting method as the aforementioned CW was applied to calculate the indicator weights of the two scenarios in the IPCC methods, respectively. Ultimately, flood risk zoning maps of the study area were constructed based on two IPCC approaches (Figure 9).

As demonstrated in Figure 9, the overall patterns of the two IPCC methods were similar to the FSE approach, and they all exhibited risk being relatively higher in the southeast of the study area than in the northwest. However, there were noticeable differences in certain regions. Compared to the FSE method, the two IPCC methods had significantly lower estimates of very high and high risk areas in the study area. Figure 5c identifies the urban regions of Beijing, Tianjin, Tangshan, and Shijiazhuang as very high risk zones, but these areas were classified as moderate risk or low risk areas in AR5 and TAR methods (Figure 9a,b). In the mountain valleys in the northwest, such as the Zhangjiakou urban area, concentration of population and property coupled with a frequent occurrence of the flash flooding will inevitably lead to severe flooding losses. Thus, Figure 5c demonstrates that the FSE approach can better identify these high to very high risk areas in river valley areas. However, the two IPCC methods showed apparent deficiencies in the identification of high risk zones in mountainous areas. Moreover, the calculation results of the two IPCC methods were relative index values rather than the final flood risk zoning map. Therefore, different classification standards will generate multiple ‘apparent real’ risk maps, which will also increase the uncertainty of evaluation results. We used the classification methods of standard deviation (Figure 9a) and defined interval (Figure 9c) to obtain two flood risk zoning maps for the IPCC AR5 approach. Compared with Figure 9a, Figure 9c shows a significant difference in the spatial distribution of risk levels.

Based on the above analysis, there are considerable uncertainties in the IPCC methodology. According to the risk definitions of the two IPCC approaches, flood risk entirely depended on the magnitude of hazard, exposure, sensitivity, and adaptive capacity. Hence, the distribution of risk indicators in the four domains (i.e., hazard, vulnerability, sensitivity, and adaptive capacity) inevitably affected the final risk assessment results. Since IPCC does not limit any fixed boundaries for assigning specific indicators to one particular domain, different researchers may switch indicators from one domain to another, which results in uncertainties in the evaluation procedure. Considering the opinions of experts, two different scenarios were created based on the AR5 method. As shown in Figure 9a,d, the risk maps generated were significantly different after moving GO from the exposure domain to the sensitivity domain. In this way, we can use the opinions of different experts to generate a large number of ‘apparently true’ risk maps theoretically. Nevertheless, from the perspective of experts, all of these maps were ‘real true’ risk maps in the study area, since experts think that their perceptions are correct. It is challenging to identify which risk map is the ‘real truth’ risk map among the number of risk assessment results. Furthermore, the risk evaluation results calculated by the two IPCC methods could not directly demonstrate the risk level results. We need to classify index values to get the final risk map. As a result, multiple ‘apparently real’ risk maps will also be constructed due to the different index classification methods. All these factors commonly contribute to the uncertainty of the two IPCC approaches.

In contrast, when the FSE approach is applied to assess risk, it can minimize various uncertainties by quantifying risks that arise from human perception [23]. The FSE method is based on the fuzzy set and fuzzy logic theories, which can transform vagueness into certainty in the risk assessment procedure. Flood risk results depend on the grade interval of each indicator and the waveform of the membership functions in the FSE approach. Furthermore, the grade intervals were divided according to the inherent properties of each indicator’s values and the actual situation of the study area. The waveform pattern of the membership functions was selected based on the intrinsic natural properties of floods. Therefore, in addition to considering the subjective weighting method in the weight calculation process, there is little scope for a human to affect the risk output. Moreover, compared with the two IPCC methods, only one risk map can be generated based on the FSE approach combined with a specific set of risk indicators, and this can be safely considered as the ‘real truth’ risk map of the study area. Subsequently, it should be noted that when we used the FSE method to conduct flood risk assessment research in this paper, we adopted the risk conceptual framework proposed by the United Nations Department of Humanitarian Affairs (UNDHA) and divided the indicators into two categories: hazard and vulnerability. The primary purpose of this method is to choose a more comprehensive and reasonable risk indicator, and further understand the impact of natural and socio-economic factors on flood risk, without influencing the final risk zoning result.

4.5. Flood Risk Mitigation

One crucial purpose of flood risk assessment is the identification of areas with potential flood risk, so that sustainable flood management strategies can be undertaken to minimize the damage caused by floods. Based on the risk zoning map constructed by the FSE approach combined with flood field investigations, we propose adopting a series of flood management and mitigation measures in high risk regions at the level of regional risk prevention. These activities include both engineering and non-engineering.

4.5.1. Engineering Measures

Engineering measures mainly refer to all kinds of buildings and facilities that serve to prevent, temporarily store, divert, and drain floodwater [49]. During our field survey, it was necessary to carry out engineering measures in this region.

(1) Dredging the channel in highly aggrading areas severely sedimentary reach: As a river flows out of the mountains and enters the plains, the hydrodynamic environment is weakened due to the widening of the river channel and the slowing down of the slope. Sediment aggradation raises the river bed and increases risk of flooding. Our field investigations revealed that many rivers in Shijiazhuang and Baoding are seriously silted up, such as the Hutuo River, Dasha River, and the Baiyangdian river system. Thus, it is essential to restore the natural condition in these ‘hotspots’ of aggradation through regular strategic dredging.

(2) Upgrading flood protection structures and remove dangerous hydraulic engineering facilities: Flood protection embankments (artificial levees) have been built along the riverine due to flat terrain and poor drainage. However, Baoding, Langfang, and the southwest mountainous areas of Beijing have comparatively low fortification criteria, and can only withstand floods with a recurrence period of 10–20 years (or even lower). It has been found that there are many dangerous sluices sick-dangerous sluices and reservoirs that have been in disrepair for a long time in the study area, such as the Yanghe reservoir (large, Zhangjiakou) and Xiyanghe Reservoir (medium, Qinhuangdao), and this poses a serious threat to downstream towns. Therefore, upgrading and removing dangerous water dam projects is necessary to improve flood prevention capacity.

(3) Delineate flood detention areas and excavate the diversion channels: Due to the superimposed effect of the flood peaks of the tributaries, the convergence area of the fan-shaped river network generally forms floodplains, such as the Baiyangdian area in Xiong’an New District. We recommend establishing flood storage and detention areas in these confluence regions. Furthermore, the field survey found that there are many tributaries in the upstream of Daqing River, Yongding River, and Ziya River, but there are few flood discharge channels downstream. Considering the issue of flood discharge, it is necessary to excavate the diversion channels in Xiong’an New District and Tianjin.

(4) Use of paleochannels as conduits for water flow during flood events: The paleochannels are old courses left after river migration. Field investigations found that there were a large number of paleochannels in Shijiazhuang, Baoding, Cangzhou, Tianjin, and Beijing. Most of the riverbeds are 1~3 m below the flat land, and these can be used as natural conduits for water flow during flood events. Hence, it is strongly recommended to select several sites as pilot projects in high risk zones to reactivate the vital role of paleochannels in flood control.

4.5.2. Non-Engineering Measures

All means to reduce flood disaster losses in floodplains through laws, policies, administrative management, economics, and technologies other than flood-control engineering are collectively referred to as non-engineering measures. In this study, we propose the following suggestions combined with the actual situation.

(1) Flood control regulations: Field investigations found that the river channels in Shijiazhuang, Baoding, and Langfang were chaotic, and that villages and farmland occupied the floodplains that severely affected the flood discharge capacities of river channels. We recommend the formulation and implementation of flood zoning management policy.

(2) Land use planning policy: The land use in the research area should be planned reasonably along with a floodplain zoning management policy. The regions that are classified as high risk zones, especially very high risk buffer zones of the river channel, should be designated for low-occupancy land-use planning. In addition, residential and commercial land should be planned to be used in areas with high altitudes, with urban planning relatively far away from river channels.

(3) Constructing a flood risk transferring mode-flood insurance: Flood insurance is an important means of flood risk management in floodplains. Essentially, flood insurance is a social or collective economic compensation method for property damage caused by flood disasters. We suggest that the government establish a flood insurance system for people living in areas with high risk levels.

(4) Capacity-building for flood resilience: The construction of flood resistance capacity involves many aspects. Real-time and accurate hydrological monitoring, forecasting, and early warning are all essential for mitigating losses in all river systems. Furthermore, evacuation drills should be conducted regularly to improve awareness of flood control and disaster reduction by community citizens. All these measures can reduce flood risk and economic losses.

5. Conclusions

Rapid urbanization and climate change have led to dramatically increased flood risks in metropolitan areas over the past few decades. Flood risk assessment as an essential means to identify regional flood risk levels is crucial for risk management in the future.

In this research, an FSE method based on combined weight was proposed and applied to risk evaluation in the BTH area. The spatial distribution of flood risk levels in the BTH area was identified from a grid scale (1 × 1 km). Research demonstrated that moderate, high and very high risk zones were mainly distributed in the east, central and southwest regions (31.36%), respectively, while low and very low risk zones were located in the plateau and mountainous regions in the northwest areas (68.64%), respectively. The area proportion of each risk level was counted at a county scale, which is of great significance for reducing flood damage and resource allocation in the BTH area.

To avoid the limitations of a single weighting approach, the combined weight (CW) approach was proposed integrating AHP with EW methods and provide a new feasible method for scientifically assigning risk indicator weights. Additionally, corresponding mitigation measures for high and very high risk zones were put forward at the regional level and combined with field investigations. This is of great significance to improving urban disaster prevention capabilities and reducing casualties and economic losses.

The FSE approach had distinct advantages compared with two IPCC methods. The FSE method is based on fuzzy set and fuzzy logic theory and can eliminate or reduce fuzziness and uncertainty in risk assessment procedures. Risk output depends on the grade interval of each indicator and the waveform of the membership function, which is entirely determined by the inherent properties of the value of each indicator and the characteristics of the flood event. Only one risk zoning map can be generated that can be safely considered to be the ‘real truth’ risk map of any study area when using the FSE approach and combining it with a specific set of risk indicators. In future study, the selection and accuracy of data indicators need to be refined (e.g., using a resolution of 100 × 100 m), and related algorithms (e.g., deep learning) should be used to construct a risk evaluation model to further improve the accuracy and objectivity of the risk evaluation results.