1. 引言
《减轻灾害风险全球评估报告》指出 [1],到2050年,估计全球人口的40%将生活在经历过严重水患的江河流域,尤其是在非洲和亚洲。对地形复杂且自然灾害频发的中国而言,洪灾更是面临的最主要的自然灾害 [2]。特别是新中国成立以后,随着泛洪区人口、社会财富等不断的聚集,洪灾给我国造成了大量的人员伤亡和财产损失。据统计,1950~2011年期间,洪涝灾害导致我国近28万人死亡,房屋倒塌1.2亿间,直接经济损失约6万亿元,农作物成灾面积3.35亿公顷 [3]。然而,作为一种非工程措施,洪灾风险评价可以为减轻各种洪水灾害损失提供有力支撑。
近年来,国内众多学者在洪灾风险评价方面做了大量的研究 [4] [5]。根据这些研究所用数据的不同,其评价方法主要可分为如下四种 [2] [6]:基于历史灾情数理统计法、基于指标体系的加权指数法、基于情景模拟分析法和基于遥感数据分析法。虽然这些方法在识别洪灾风险中起到了一定的作用,但是,也都有着各自的应用局限性。例如:基于历史灾情数理统计法要求大量的历史灾情数据、不能反映洪灾风险的空间差异性;基于指标体系的加权指数法通常不考虑洪水的淹没特征,另外,指标权重的计算方法有很多种,比如常用的有层次分析法、熵权法等,这也一定程度上造成了评价指标权重的不确定性;基于情景模拟分析法不考虑承灾体和孕灾环境等。此外,对于上述方法来说,当前也没有统一的洪水风险等级区化标准。为了能够综合指标体系法和情景模拟的优点,改善指标权重,Hongshi Xu等人 [7] 在情景模拟的基础上,构建了风险评价指标体系,并应用乘法合成法组合了层次分析法权重和熵权法权重,提出了改进熵权法,计算了各评价指标的权重,最后以此评价了海口地区的城市洪水风险。但是,李刚等人 [8] 的研究表明乘法合成法在一定情况下会出现“倍增效应”。而由级差最大化合成法组合的主、客观权重既不会出现乘法合成法那样的“倍增效应”,又能对评价对象进行有效的区分。因此,本文拟构建综合考虑情景模拟结果的基于指标体系的洪灾风险评价方法,并通过级差最大化合成法求各评价指标的组合权重。
本文在现有研究的基础上,旨在从致灾因子、孕灾环境和承灾体三方面来构建综合考虑情景模拟结果的基于指标体系的溃堤洪水风险评价方法,并研究评价指标权重的计算方法,探索K-Means聚类算法在洪水风险等级区划中的应用,最后以此作为方法框架对浑太胡同进行溃堤洪水风险评价研究。结果表明,该方法方能够对洪灾风险进行有效的评价。
2. 研究区与数据
浑太胡同防洪保护区位于辽宁省,“浑太”指的是浑河和太子河,如图1所示,自小北河起到下游的三岔河,浑河和太子河大体呈平行状向西南流去,在两河的堤防间形成了一块封闭的类似胡同状的区域(以下简称:浑太胡同),该区域的面积约为300 km2。

Figure 1. Sketch map of Hun-Tai flood prevention area and rivers
图1. 浑太胡同防洪保护区及河流示意图
浑太胡同涉及到的乡镇级行政单位包括辽阳县境内的小北河镇、唐马寨镇及海城市境内的高坨镇、温香镇。根据2012年的统计数据 [9],区域内人口、经济等数据如表1所示:

Table 1. Social and economic statistics of Hun-Tai flood prevention area
表1. 浑太胡同防洪保护区社会经济数据
模型输入的数据主要有:DEM数据、糙率数据和入流边界条件。研究区的地形数据来源于松辽委,如图1所示,DEM的分辨率为30 m,本文对研究区域进行了网格剖分,共划分了大小平均约为10,000 m2的网格31,273个(下文所指的评价对象);参考对研究区洪水模拟的已有成果及研究区的实际情况,将居民地、水田、旱田、树丛和空地5种不同土地利用类型的糙率值分别设置为0.070、0.050、0.060、0.065和0.035;对于入流边界条件,根据研究区堤坝和洪水特点,本次风险评价选取迟拓子与马家堡子组合溃口作为典型溃口,溃口位置如图1所示,并将组合溃口100年一遇的流量过程作为模型的入流边界条件。各溃口的流量过程如图2所示:

Figure 2. The 100-year flow process of the combination rupture
图2. 组合溃口100年一遇流量过程
3. 研究方法
本文基于指标体系和情景模拟的风险评价方法,首先,利用HEC-RAS软件构建研究区溃堤洪水演进模型来获取致灾因子指标,通过灾损评估得到承灾体洪灾经济脆弱性指标,并利用ArcGIS处理研究区地形数据得到孕灾环境指标;其次,由级差最大化合成法对层次分析法权重和熵权法权重进行优化组合,得到各评价指标的组合权重,并计算各评价对象的综合得分;最后,通过K-Means聚类算法对各评价对象的综合得分进行聚类分析,可得到研究区的洪水风险图。洪水灾害风险评价的思路如图3所示:

Figure 3. The steps of flood risk assessment
图3. 洪水灾害风险评价
3.1. 洪水二维演进数值模拟
蓄滞洪区的洪水要素分析方法可分为水文学法和水动力学法 [10],其中,水文学法主要有时段单位线法、线性或非线性水库法等,而水动力学法可通过求解浅水方程来进行洪水演进模拟,浅水方程计算公式如下 [11]:
连续方程:
(1)
动量方程:
(2)
式中,H为水面高程(m);t为时间(s);h为水深(m);V为水流速度(m/s);q为旁侧入流量(m2/s);g为重力加速度(m/s2);
为运动粘度(m2/s);
为摩擦系数;f为科里奥利系数;k为单位矢量。
HEC-RAS作为一款水动力学软件,其2D模块便是基于浅水方程原理开发的。该软件的主要功能包含以下四个方面:① 恒定流水面计算;② 1D和2D非恒定流模拟;③ 可变边界的泥沙输移计算;④ 水温和水质传输模拟。利用以上功能,HEC-RAS既可以模拟桥梁、溢洪道、堰、闸门、滞洪池等水工建筑物对水流的影响,也可以对洪水进行二维演进数值模拟。
3.2. 指标构建与权重计算
1) 指标体系
根据研究区的特点,参考国内外关于洪水灾害风险评价的研究成果,本文在选择评估指标时,综合考虑了洪水对自然、经济、社会等方面的影响,从致灾因子、孕灾环境和承灾体脆弱性三个方面出发,确定了用以评价研究区溃堤洪水风险的指标体系。
指标体系可以划分为目标层、准则层和指标层,其中,目标层是以浑太胡同溃堤洪水风险评价为总目标,准则层包含致灾因子、孕灾环境和承灾体。致灾因子指标包括:最大淹没水深(MD)、洪水到达时间(AT)、洪水淹没时间(DT);孕灾环境指标包括:地表高程(DEM)、地形坡度(SL);承灾体指标:洪灾GDP损失(LOP)。根据指标与洪水灾害风险的相互关系,可以将指标划分为正向指标或负向指标,风险随正向指标的增大而增加,随负向指标的增大而减小。指标体系如表2所示:

Table 2. The index system of risk assessment
表2. 风险评价指标体系
2) 指标权重
指标权重的计算方法可以分为主观赋权法、客观赋权法和主客观组合赋权法 [12] [13]。李刚等人 [8] 的研究表明,综合考虑主观权重和客观权重,计算评价指标的主客观组合权重,更能对评价对象进行合理的评估或对备选方案进行有效的区分。因此,本文将采用级差最大化合成法来对由层次分析法计算所得的主观权重和由熵权法 [14] 计算所得的客观权重进行优化组合,以得到各评价指标的组合权重。
级差最大化组合赋权模型以各评价对象的区分度最大为目标函数、以各指标权重取值区间与指标权重和等于1为约束条件来对各指标的权重进行优选,其中各评价对象的区分度可以用其综合得分的方差s2来度量。记k个评价对象的综合得分为
,k个评价对象n个评价指标的标准化矩阵为X,n个评价指标的权重为
,
为各评价对象得分的平均值,则有:
(3)
第i个指标的权重
的取值范围
由层次分析法计算所得的权重和熵权法计算所得的权重共同确定,进而,组合权重优化模型的数学表达式如(4)所示:
(4)
3.3. 基于K-Means的风险等级聚类分析
K均值聚类算法 [15] (k-means clustering algorithm,简称K-Means)是由James MacQueen于1967年提出的一种迭代求解的聚类分析算法,该方法被广泛应用于数据分析、信号处理和机器学习等领域。在利用K-Means进行聚类分析时,关键的一步是确定聚类数目。为此,本文先采用手肘法确定聚类数目,然后再通过K-Means聚类算法对研究区溃堤洪水风险进行聚类分析。应用手肘法确定聚类数目的核心指标是误差平方和(SSE, sum of the squared errors),其计算公式如下:
(5)
式中,k为聚类数目,
为第i个样本簇,p为
的样本点,
为
的质心。SSE随k值变化的规律是:k取值越大时,样本数据便被划分的越细,每个样本簇的聚合程度就越高,进而误差平方和SSE变小。根据k与SSE的这种关系便可确定K-Means的聚类数目。因k与SSE的关系曲线形似手肘,故而该方法被称为手肘法。
4. 评价结果与分析
4.1. 模拟结果合理性分析
2016年,中国水利水电科学研究院等单位对浑太地区的洪水风险进行了研究。本文根据水科院等单位的研究成果报告 [9] 来对HEC-RAS模型的参数进行率定。本文利用HEC-RAS对如下四个溃口进行了洪水演进模拟:迟拓子、马家堡子、刀把子和西高溃口,当各溃口均输入与报告相同的百年一遇流量过程时,HEC-RAS模型模拟所得的最大淹没范围分别为:284.43 km2、312.56 km2、159.15 km2和226.31 km2,与报告中已知淹没范围的相对离差分别为:0.06%、1.90%、−2.06%和−0.51%。由此可见,模拟结果与已知结果基本吻合,所以认为HEC-RAS模型的参数设置合理,模拟结果合理。
4.2. 风险评价指标量化
通过HEC-RAS对研究区的溃堤洪水演进进行模拟,可直接获得致灾因子指标MD、AT和DT。基于研究区洪水淹没计算结果,结合不同财产的灾前价值,根据分类财产损失率与水深的关系 [9],可估算出当迟拓子与马家堡子组合溃口发生100年一遇的流量过程时研究区的直接经济损失情况。间接经济损失 [16] 是指洪水次生灾害和衍生灾害所造成的损失,主要包括因采取各种防洪措施而增加的费用、因骨干交通路线的中断造成相关生产加工企业原材料中断而停产或绕道取材增加的成本、因农产品减产给相关加工企业造成的损失等。由此可见,间接经济损失涉及面广,计算较为繁杂,因此本文将按直接经济损失的20%来估算间接经济损失 [17]。根据计算所得的各种经济损失,便可得到研究区内各类承灾体的洪灾经济脆弱性指标LOP。研究区当迟拓子与马家堡子组合溃口发生100年一遇的溃堤流量过程时,各指标的空间分布情况如图4所示:


Figure 4. Spatial distributions of risk indices when 100-year flow process occurred at combination rupture
图4. 迟拓子与马家堡子组合溃口发生100年一遇流量过程时研究区各指标空间分布
4.3. 风险评价指标权重
为了计算层次分析法权重,本文通过参考其它相关研究成果和专家的意见,对比了各评价指标的相对重要程度,形成了用以计算指标层次分析法权重的判断矩阵,判断矩阵如表3所示。

Table 3. The judgment matrix for AHP
表3. 层次分析法判断矩阵
由各种权重计算方法求得的评价指标权重结果见表4。从表中的计算结果可以看出,由层次分析法和熵权法计算所得的LOP指标权重均排序第一,熵权法计算得到的SL指标权重远高于由层次分析法计算所得的权重;级差最大化合成法综合考虑了主、客观权重两方面的因素,所得权重中,LOP、MD指标权重分别排序第一、第二。

Table 4. The weight results of risk evaluation indices
表4. 风险评价指标权重计算结果
4.4. 洪水风险评价结果
图5为各评价对象的综合得分经手肘法分析后的结果。从图中可以看出,对于级差最大化合成法来说,当k大于4时,随着k值的增加,SSE的下降幅度越来越小。为了方便区划研究区洪水风险等级,可将K-Means聚类数目k定为5,也即研究区溃堤洪水的风险等级可划分为5级:I最高、II较高、III中等、IV较低、V最低。将k和存有各评价对象综合得分的数据文件输入到K-Means计算程序中进行聚类分析,由分析结果可知,当迟拓子与马家堡子组合溃口发生100年一遇的溃口流量时,各级洪水风险的取值范围分别为:I (0.5012, 1)、II (0.3812, 0.5012)、III (0.3369, 0.3812)、IV (0.2591, 0.3369)、V (0, 0.2591)。当迟拓子与马家堡子组合溃口发生100年一遇溃口流量过程时,图6中的(a)为对应求得的研究区洪水风险图。

Figure 5. The relation curve between SSE and k
图5. SSE与k的关系曲线
(a) 级差最大化法 (b) 层次分析法 (c) 熵权法
Figure 6. The flood risk map of the study area obtained by different weight calculation methods
图6. 由不同权重计算法求得的研究区洪水风险图
4.5. 评价结果分析
1) 评价结果合理性
图6中(b)和(c)分别为由层次分析法和熵权法求得的研究区洪水风险图。对比图6中的(a)和(b),可以看出由级差最大化组合权重求得的风险图和由层次分析法所求得的风险图整体趋势基本相同,但(b)图在研究区上游的中等、较高风险区域面积小于(a)图的;对比图6中的(a)和(c),可以看出(a)和(c)的上半部分基本相同,但下半部分差别较大;而图6中的(b)和(c)整体差别较大。由此可见,利用级差最大化合成法组合层次分析法权重和熵权法权重,求各评价指标的组合权重,能够改善单一权重法的不足,使得评价结果能够兼顾主、客观两方面的因素。
2) 洪水风险分析
从图6中的(a)可以看出,研究区溃堤洪水风险整体呈现由上游到下游增大的趋势,这与溃堤洪水水深由上游到下游呈增大的趋势一致;研究区洪水风险等级最高的区域多为居民居住地,较高风险区域多为下游受淹没的水田和旱田,而风险最低的区域多处于研究区上游,这主要与研究区上游地势较高、水深较浅有一定的关系。表5总结了研究区各级洪水风险区域的面积和占比情况,从表可以看出,中等风险区域和较高风险区域的面积和占到了研究区总面积的89.64%。

Table 5. Area and the proportion of different risk levels in the study area (km2)
表5. 研究区各风险等级区域统计(面积单位:km2)
3) 评价指标分析
图7为不同风险级别下各指标的平均取值,图8为图6中(a)所示区域各指标的取值情况,各区域均由程序随机选取得到,其中,1和2为最低风险区域,3和4为较低风险区域,5和6为中等风险区域,7和8为较高风险区域,9和10为最高风险区域。从图7可以看出,风险最高区域的洪灾经济损失最高,从前边分析知,这些区域多是人们生产和生活的聚集地,在研究区内呈现点状分布。从图7和图8可以看出,研究区整体的洪水风险等级、不同区域的洪水风险等级不依某一单一指标取值的增减而呈现增大或减小的趋势,体现了不同区域的风险等级是各评价指标综合作用的结果,比如,虽然区域4、5的DEM取值相同,但二者的风险级别不同。

Figure 7. Mean index values under different risk levels
图7. 不同风险等级下各指标的平均取值

Figure 8. Values of each indicator for representative regions
图8. 代表区域各指标取值
4) 权重组合方法对比
通过乘法合成法 [8] 对表4的层次分析法权重和熵权法权重进行组合,得到指标MD、AT、DT、DEM、SL和LOP的组合权重分别为:0.2440、0.0299、0.1121、0.0562、0.0177、0.5401,所得指标组合权重由于“倍增效应”的影响,使得LOP指标的组合权重远高于层次分析法权重或熵权法权重,因此弱化了其它指标在评价体系中的作用,权重结果的合理性也难以解释。从图5可以看出,对于同样的研究对象,当k小于3时,由乘法合成法所得权重求得的SSE值高于由级差最大化合成法所得权重求得的SSE值。
5. 结论
本文利用HEC-RAS模型对研究区当迟拓子与马家堡子组合溃口发生100年一遇溃堤流量时进行了洪水演进模拟,构建了综合考虑情景模拟结果的基于指标体系的溃堤洪水风险评价方法,利用级差最大化合成法,求得了用以评价研究区溃堤洪水风险的各评价指标的组合权重,并利用K-Means聚类算法对研究区洪水风险进行了聚类分析。得出以下结论:
1) 当迟拓子与马家堡子组合溃口发生100年一遇的溃口流量过程时,研究区的洪水风险整体呈现由上游到下游增大的趋势,风险最高的地方多为居民聚居地,中等风险和较高风险区域占到了研究区总面积的89.64%。
2) 利用乘法合成法组合主、客观权重时,可能会出现倍增效应,使得某一指标的组合权重远高于由层次分析法或熵权法所得的该指标的权重,因此弱化了其他指标在评价体系中的作用。
3) 由级差最大化合成法得到的组合权重能够对评价对象进行有效区分,且各区域的风险高低不依某一单一指标的增减而降低或升高,而是各指标综合作用的结果。
基金项目
该研究是在中央高校基本科研业务费(DUT18RC(3)072)的资助下完成的。