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

: 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, ﬂood disasters in China have increasingly a ﬀ ected the sustainable development of metropolitan areas. Risk assessment is an essential step in ﬂood management and disaster mitigation, which provide a quantitative measure of ﬂood risk. However, the di ﬃ culty of ﬂood 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 ﬂood 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 ﬂood 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 ﬂuvial 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 ﬁeld investigations. The veriﬁed risk assessment results were spatially consistent with the historical ﬂood disaster records and loss positions, indicating the accuracy and reliability of the risk assessment map using the FSE approach. the IPCC (Intergovernmental Panel on Climate Change) TAR (Third Assessment Report) and AR5 (Fifth Assessment Report) methods, FSE has signiﬁcant advantages in handling uncertainty, complexity, and the non-linear relationship between indices and risk grades. This study provides a novel quantitative method for ﬂood risk assessment in metropolitan areas and practical implications for urban ﬂood management.


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.

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.

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 2 ): 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, a large code value indicates a strong infiltration capacity.

•
Gross domestic product (GDP, yuan/km 2 ): 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 2 ): This indicator reflects the spatial distribution of human exposure in the research area in 2015.

•
Grain output (GO, kg/hm 2 ): 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 2 ): 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. 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).

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.

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.

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.

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.

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.

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, uij is the membership degree of the indicator i at the grade j. It is not difficult to find that the actual data xi of an indicator that has two membership degrees (ui2 and ui3). However, it can be derived that xi 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.

Comprehensive Evaluation
Once the membership functions are determined, a membership function matrix (Rk) is generated according to Equations (A1)-(A6) that combines the actual values of all indicators. The resulting vector Bk is the product of the weight vector W and the fuzzy relation matrix Rk, and its mathematical derivation is as follows: where W is the weight vector of indicators (its calculation process is described in detail below). The element of Bk 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.

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

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.

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

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 2 ) 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 2 ) 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 2 ) 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 2 , 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 Figures 6 and 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.

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 (Tables 6 and 7). 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 (Table 8). 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.

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 where 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).
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).
Therefore, the calculation formula of two IPCC methods are as follows: 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). 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 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.

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.

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.

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.

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.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Membership Functions of Flood Risk Assessment in the BTH Area
Here, x i is the original data of evaluation units in each indicator.

Appendix B. Detailed Calculation Steps of Combined Weight for Risk Assessment
(1) The judgment matrix X = (x ij ) m×n is constructed with m evaluation indicators and n flood risk assessment levels.
x m1 x m2 · · · x mn (2) Considering different indicators have different orders of magnitude and dimensions, it is necessary to standardize indicators so that they are under the same evaluation system. Normalized processing for positive and negative indicators use the formula (16)-(17), respectively.
where, max(x ij ) and min(x ij ) represent the maximum and minimum values of an indicator, respectively.
(3) The entropy value e j of each indicator can be calculated by the following formula: where, p ij = y ij / m i=1 y ij and p ij lnp ij is defined to zero when p ij is equal to zero.
(4) We use the entropy value to calculate the entropy weight of the evaluation indicator: where, j = 1, 2, . . . , n and satisfy n j=1 z j = 1.
(5) Then, the analytic hierarchy process moderate to calculate the indicator weights and the consistency check is performed using the following formula. Detailed calculation steps of this method can refer to the research of [16,20].
where, CR represents the consistency ratio and CI indicates the consistency index. RI is the random consistency index that is computed by Saaty from a sample of 500 matrixes randomly generated [25]. (6) Finally, combining the weight results of the AHP and EW methods, formula (16) is used to calculate the combined weight of the risk indicators.
where, h j and z j represent the subjective and objective weights calculated by using AHP and EW methods, respectively.