蓄滞洪区洪水演进模拟及堤防溃决损失评估方法

蒋水华， 黄中发， 黄劲松， 张 颖， 江先河

(1.南昌大学 建筑工程学院， 江西 南昌 330031； 2.江西水利职业学院， 江西 南昌 330013)

1 研究背景

蓄滞洪区灾害模拟、风险分析、预警，以及洪水演进分析、评价、决策一直都是水文水资源领域和防灾减灾部门关注的重大课题[1]。蓄滞洪区在洪水来临期间降低下游洪灾损失的同时，也要考虑降低蓄滞洪区内进洪过程中的损失，而且随着时间的推移，蓄滞洪区内产业经济已得到发展，人口也有所增长，分洪时将造成较大损失，故有必要模拟蓄滞洪区洪水演进过程，估计对应的生命、经济和生态环境损失，从而指导防洪抢险。

洪水演进数值模拟是洪水风险分析及其损失评估的关键步骤，是进行洪水风险管理的中心环节[2]。童汉毅等[3]、Tucciarelli等[4]、李大鸣等[5]分别采用有限差分法、有限元法、有限体积法模拟蓄滞洪区洪水演进过程，取得了较好的数值结果，但是采用这些方法建立的模型推导和计算过程复杂，计算量较大。随着计算机和地理信息系统技术的发展，采用一维或二维水力学方法模拟洪水演进过程发展迅速，水文学方法得到较为广泛的应用[6]。魏凯等[6]、穆聪[7]等、王崇浩等[8]、郭凤清等[9]、朱世云等[10]、袁雄燕[11]等、Morales-Hernández等[12]和Bladé等[13]利用一维或二维水力学方法建立MIKE21 FM模型，分别模拟了蓄滞洪区、大坝溃坝、河湖和堤防溃决的洪水演进过程，结果表明二维水力学模拟方法可在短时间内直观地显示洪水淹没范围的变化，提供较详细的淹没区信息，可较好地根据数值模拟结果制定相关决策指导防洪抢险工作。黄琳煜等[14]基于MIKE FLOOD建立MIKE URBAN与MIAKE 21耦合模型，并构建了浦东川沙地区耦合模型来分析出现涝点和积水的原因，取得了较好的研究成果。王扬等[15]归纳总结了MIKE内部各个模型的基本原理、特点及其适用范围，并且指出了MIKE模型存在的局限性。

针对蓄滞洪区堤防的规模较大，而现有的防洪标准较低、堤身堤型不达标、区域内人口密度大、水保工程措施不到位，且区内产业严重失衡、种植业较多以及蓄洪排涝成本较大等因素引起的洪水演进过程中损失估算的难题。本文发展了堤防溃决洪水演进数值模拟方法，建立了基于蓄滞洪区洪水演进淹没数据的损失评估方法。依托鄱阳湖康山蓄滞洪区，采用MIKE21建立蓄滞洪区洪水演进数值模型，进而根据洪水演进过程模拟结果(淹没范围、水深和流速)估计堤防溃决造成的生命、经济和生态环境损失。

2 洪水演进数值模拟方法

丹麦MIKE21模型[8-10,16]是常用的洪水演进数值模拟软件，分为海岸水文学和海洋学、环境水文学、泥沙传输过程和波浪4个模块，其原理是采用二维水动力学数学模型模拟水流运动，直观得出水位、流量和流速等重要的水力指标随时间的变化关系。相比于其他模型，MIKE21模型的功能较强大，不仅可以进行软启动，还可以设置多种控制性水工结构，进行干、湿节点及单元设置等[11]，在我国乃至世界范围内得到了广泛使用，比如在我国南水北调工程和长江综合治理等工程中均得到了成功应用。

2.1 控制方程

MIKE21模型将河道水流看成不可压缩的牛顿液体，采用Navier-Stokes方程组来描述河道水流运动规律[17]。为了考虑河床底部摩阻作用和紊动影响，引进涡黏系数，建立平面二维水动力数值模型。模型控制方程包括平面二维流连续方程，计算表达式为:

(1)

式中：t为计算时间，s；x、y、z为右手Catesian坐标系；h为静止水深,m；u、v为流速在x、y方向上的分量，m/s；S为源汇项。x方向上的二维水流动量方程为:

(2)

2.2 边界条件

边界条件主要包含上、下边界条件、特殊边界条件以及动边界条件，具体如下：

(1)上边界条件。蓄滞洪区上边界条件一般选取边界处的流量或水位过程，根据研究目的不同，可选取历史最高水位或流量过程、设计水位或流量过程、实测时段水位或流量过程等作为对应的上边界条件。

(2)下边界条件。通常滞洪区在运行时为封闭区域，所以模型不设置下边界条件；若不是封闭区域，下边界条件一般为边界处的水位或流量过程。

(3)特殊边界条件。对起到排水作用的构筑物(涵管、涵洞等)和阻水作用的建筑物(道路、堤防等)进行概化，使得模拟更加符合客观实际。

(4)动边界条件。动边界是指计算区域中随水位变化的边界，为了保证模型计算的连续性，采用MIKE21模型中的干、湿边界来设置动边界条件[9]。一般来说，干湿水深分别采用系统默认值0.005 m与0.1 m可以满足计算精度，具体含义为当计算区域的水深小于0.005 m时，该边界条件则为干边界，不参与计算；当水深大于0.1 m时，该边界条件则为湿边界，重新参与计算。

2.3 参数选取

MIKE21模型参数包括数值参数和物理参数。数值参数一般取默认值，物理参数主要包括河床糙率系数、动边界计算参数和涡黏系数等。MIKE21水动力数学模型的主要参数为糙率系数，该参数是一个反映水流阻力的综合系数。研究区下垫面糙率系数可依据《洪水风险图编制导则》[18]给出的糙率变化范围取值或考虑研究区域河床现状进行率定取值。

2.4 溃口设置

蓄滞洪区分洪口门包括自然分洪、闸门分洪和爆破分洪这3种分洪方式，都须预留分洪口，根据国家防汛抗旱总指挥部《关于长江洪水调度方案的批复》(国讯[2011]22号)[19]，当湖口水位达到分洪水位时，采取相关工程措施将堤防预留的分洪口扒开，蓄滞洪区开始进洪。参照分洪口门的设置参数，进口形态接近流线型，故可采用流线型宽顶堰公式计算溃口流量[20]：

(3)

式中：Qb为溃口处出流流量，m3/s；Z为溃口处河道水位,m；m和σ分别为流量系数和淹没系数，根据文献[20-21]计算确定，没有考虑淹没系数与流量系数存在的不确定性；Zb为溃口顶部高程，m；B为溃口宽度，m。

3 堤防溃决损失评估方法

当发生超标准洪水时，为防洪抢险，将预留的堤防分洪口扒开，进而通过堤防分洪口的洪水会对下游蓄滞洪区造成一定的经济、生命、环境损失。生命损失是指堤防下游地区受溃堤洪水影响造成的死亡或受灾人口数量。经济损失是指溃堤洪水带来的直接经济损失和间接经济损失，可用货币量表示。生态环境损失指的是堤防溃决引起的资源破坏和环境污染对生态环境造成的影响。

关于生命损失，结合宋敬衖等[22]的研究成果，根据人口聚居特点结合行政区划将淹没区划分为若干子区域，使用风险人群积分算法得到实用的生命损失LOL计算公式如下：

(4)

式中：PARij为第i个区域中第j组风险人口数量；a为风险区域总数；c为风险人口分组数；IRij为第i个区域中第j组风险人群个体生命损失率，计算公式为:

(5)

式中：Pi为溃堤概率；fi为发生第i等级洪水个体死亡率，可根据洪水演进模拟的淹没范围、水深以及预警时间等因素查表得到；l为洪水等级数目[22]。此处fi为传递参量，需要根据淹没区受灾程度以及预警时间来确定。

经济损失主要包括提防的破坏、下游淹没区城镇的损害、工程收益的损失等。根据损失特征采用损失率方法估算直接经济损失，该方法适用于估算各类流动资产与固定资产的损失[22]，间接经济损失一般取直接经济损失的0.63倍[23-24]。

(6)

式中：S为按损失率计算得直接经济损失；Si为第i类财产损失；βijk、Wijk分别为第i类第j种财产在第k种淹没程度下的损失率和财产价值；e为财产分类总数；h为第i类财产分类总数；l为淹没水深等级。此处βijk、Wijk为传递参量，Wijk需要根据受灾程度选择相应研究范围，再到经济库行政面积中去查找当地经济数据确定；βijk需要根据淹没程度确定。

国际上对堤防溃决造成的环境损失研究较少，并且由于考虑的因素多且复杂，没有统一的风险标准，一般较难直接统计计算。为此，本文采取条件价值评估方法评估堤防溃决造成的环境损失。条件价值评估方法是评估环境非使用价值的主要方法[25]，利用该方法评价环境损失的本质就是通过问卷调查方式来调查风险区内人口愿意为改善该环境功能的支付意愿，或者是放弃该环境功能的补偿意愿，以此揭示环境对于被调查者的价值，即该区域环境资源的经济价值[25]，其中人均支付意愿E(元/a)和环境损失量F的计算公式为:

(7)

F=E·N·f

(8)

式中：Ai为每人愿意支付的金额；Di为个体选择该金额的概率；u为金额的具体数量；N为年限；f为人口数。此处E和F为传递参量，需要根据淹没历时来计算环境损失随时间的变化关系。

4 工程应用

4.1 工程概况

康山大堤兴建于1966年，属于四级堤防，全长36.25 km。设计洪水标准：保证水位22.50 m(吴淞高程，湖口站水位)，堤顶高24.55 m，内外边坡1∶3，防汛警界水位19.50 m。穿堤建筑物有梅溪咀闸(中型水闸，设计排涝流量为103 m3/s)、锣鼓山电排站(大(2)型泵站，设计排涝流量73.8 m3/s)、大湖口闸(中型水闸，设计排涝流量为160.5m3/s)、落脚湖电排站(中型泵站，设计排涝流量15.5 m3/s)、里溪闸(小型水闸，设计排涝流量为7.1 m3/s)[26]。

康山大堤内垸1985年被国务院批准列为国家蓄滞洪区，有效蓄洪量16.58×108m3，是江西省最大的蓄滞洪区，也是国家重点蓄滞洪区，担负着15×108m3分蓄洪任务。总集雨面积450.31 km2，蓄洪面积312.37 km2。区内有石口镇、三塘乡、大塘乡、瑞洪镇、康垦总场、康山乡等6个乡镇场和1个现代农业示范区，常住人口2万余户、10余万人，耕地面积117.21 km2，养殖水面面积139.01 km2。经统计，康山大堤保护面积343.4 km2，保护耕地面积94.2 km2，保护人口103 437人、22 560户，康山蓄滞洪区居民土地财产登记情况有：承包土地262.49 km2，包括农作物117.21 km2，水产类139.01 km2，经济林6.27km2；专业养殖家畜类11 919头，家禽类198 384只；住房19 432间，合计843 353 m2；农业机械7 244台，役畜类1 395头，居民主要耐用消费品17 250台[26]。

4.2 研究区地形概化及网格剖分

康山大堤溃堤洪水演进过程模拟的计算区域基础地理信息数据比例尺为1∶10000，行政区划数据比例尺为1∶50000。所收集的电子地图需按要求用国家大地2000坐标系(CGCS2000)，高程基准统一采用1985年国家高程基准。其中1∶10000的数字高程模型如图1所示。为了分析溃堤洪水演进过程，获得溃口流量变化过程和堤防保护区淹没信息，采用公式(1)和(2)建立康山蓄滞洪区水流运动平面二维非恒定数学模型。

模型计算范围为整个康山蓄滞洪区，以北面康山大堤的分洪口为进洪边界与退洪边界，以北面康山大堤及南面陆地为物理边界。由于计算区域地形复杂，为了合理布置网格，采用MIKE21模型中的FM非结构化网格(不规则的三角形网格)类型，利用有限体积法进行数值模拟。最大网格面积不超过0.01 km2，对居民区、建筑物密集区和种植区的计算网格适当进行加密处理。研究区域共剖分了58 095个网格节点和63 880个计算单元，网格中心点之间的最大间距为300 m，最小间距为26 m，蓄滞洪区的网格剖分如图2所示。所涉及的物理参数取值如表1和2所示，具体包括边界条件、初始条件、涡黏系数和经率定后的糙率等。需要说明的是，表1中采用1954年4月15日0时至1954年4月18日0时的洪水数据进行洪水过程模拟，这是因为1954年长江中下游发生了近100年最大洪水，洪水数据非常具有代表性。

表1 模型计算物理参数取值

表2 研究区域内不同下垫面糙率取值

图1 康山蓄滞洪区数字高程模型图(1∶10000)

图2 蓄滞洪区网格剖分图

4.3 洪水溃决过程及流量分析

康山蓄滞洪区在康山大堤预留有分洪口门。分洪口位于大堤桩号20+070～20+450 m之间，分洪口门宽380 m，在分洪口门两端布置深层水泥搅拌桩以防止口门无限制扩大。

根据国家防汛抗旱总指挥部《关于长江洪水调度方案的批复》(国讯[2011]22号)[19]，当鄱阳湖湖口水位达到20.59 m(吴淞高程22.50 m)时，对应的康山站水位20.68 m(吴淞高程22.55 m)，采用人工爆破的方式将预留溃口扒开，口门底宽为12.5 m，堤防溃决，蓄滞洪区开始进洪。口门爆破并冲刷扩大的过程按2 h计，最终口门底宽取300 m。参考《江西省鄱阳湖蓄滞洪区安全建设工程可行性研究报告》[26]关于规划口门底高程的选取，取口门底高程15.93 m(黄海高程)。溃决口门发展所需时间为口门初始形态发展到终止形态的时间，其中口门结构参数和口门发展时间需要结合当地或类似地区已发生过的实际溃决调查资料合理确定。不考虑区间降雨计算口门的发展过程，得到溃口流量变化过程曲线见图3。由图3可见，当防洪堤溃决后，溃口流量在2 h内很快达到最大值14 000 m3/s。24 h以后溃口流量开始下降，至48 h完成整个进洪过程，蓄滞洪区内水位达到20.68 m(黄海高程)，与外湖水位一致，溃口不再进流，流量为0。

图3 康山大堤溃口流量过程线

4.4 洪水演进过程模拟

为了便于分析，在蓄滞洪区的康山乡、瑞洪镇、大塘乡、三塘乡、石口镇古竹片区、康山垦殖场6个行政区的23个行政村(自然村)内布置了23个特征点，各特征点布置及溃口位置见图4。堤防溃决后研究区各典型时段的淹没深度见图5。

图4 康山蓄滞洪区特征点布置及溃口位置

通过模型计算得到堤防溃决后研究水域各典型时段的淹没水深，见图5。由图5可知，堤防溃决1h后，水流到达康山垦殖场和大湖管理局，12 h后，洪水可以推进到瑞洪镇片区的东源村，瑞洪镇片区的东三村由于地势较高，洪水需要在24 h左右才能达到淹没深度。相比之下，地势较高、远离龙口的大塘安全区在36 h后洪水才到达。由于康山垦殖场落脚湖分场和康山乡的临湖区域在口门附近，这两处的最大流速大于1.0 m/s，受洪水冲刷影响较大。瑞洪镇江家村附近部分地形高程小于20 m，当蓄滞洪区分洪36 h后，西面安全区开始受影响。由图5(c)可知，分洪结束时蓄滞洪区平均淹没水深达6.95 m，康山垦殖场及石口镇部分区域淹没水深达到8 m。

图5 堤防溃决后研究水域各典型时段的淹没深度

4.5 溃堤洪水淹没面积

康山大堤分洪口门溃决时，溃口宽380 m，口门底高程15.93 m。溃决前，康山蓄滞洪区内蓄洪底部水位为15.10 m，约69 km2(去除水域面积)被渍水淹没。整个溃堤洪水演进过程历时48 h，进洪总量为16.58×108m3，总的淹没面积为288 km2。

按照淹没严重程度对淹没区域进行分区，并统计各分区的淹没面积如表3所示。定义危险区淹没深度大于2.0 m，重灾区淹没深度在1.0～2.0 m之间，中灾区淹没深度在0.5～1.0 m之间，轻灾区淹没深度小于0.5 m，安全区洪水未到达区域，该区域是人员转移安置区。

由表3可知，轻灾区、中灾区、重灾区所占淹没面积的百分比较小，分别为0.66%、0.82%、1.77%。相比之下，危险区淹没面积为278.63 km2，占总淹没面积的96.75%，说明在这种工况下，康山蓄滞洪区受灾程度大且涉及范围广。将MIKE计算结果文件导入ArcGIS，计算得到洪水未到达区域的面积为31.3 km2，该区域亦为安全区，可作为人员转移安置区。

通过数值模拟可以得出不同淹没程度的面积分布情况，从而可为康山蓄滞洪区洪水风险分析、抗洪抢险以及人员转移安置提供较为准确的参考依据。

表3 溃堤洪水淹没面积统计表

4.6 特征点淹没水深、流速

以康山乡、瑞洪镇、大塘乡、三塘乡、石口镇古竹、康山垦殖场片区内23个特征点为代表，分析蓄滞洪区淹没水深及流速变化情况，如图6所示，从图6中可获得洪水到达特征点需要的时间、特征点的最大淹没水深、流速峰值和峰现时间等。

由图5和6可知，康山乡虽离溃口较近，但地势较高，直至堤防溃决后9 h 40 min洪水才达到府前村，最大淹没水深和流速峰值均出现在府前村，分别为2.65 m和0.0133 m/s，峰现时间为溃堤后30 h 20 min，大山、金山村未被洪水淹及。

在溃堤9 h 10 min后洪水达到瑞洪镇西岗村，瑞洪镇最大淹没水深出现在西岗村为4.55 m，最大流速出现在罗家村为0.0536 m/s，峰现时间为溃堤后24 h 30 min。其中该片区的把下村因地势较高，未被洪水淹及。

图6 不同片区特征点淹没水深及流速过程线

大塘乡在溃堤50 min后开始淹没，最大淹没水深和最大流速均出现在大湖村，分别为7.65 m和0.0731 m/s，峰现时间为2 h 40 min。其中该片区陈家塘村、南垅村因地势较高，未被洪水淹及。

石口镇古竹片区在溃堤50 min后开始淹没，最大淹没水深出现在古竹村为6.24 m；最大流速出现在西村为0.027 m/s，峰现时间为溃堤后28 h 50min。刘埠村和湖滨村由于地势较高，直至溃堤后48 h进洪完成，未被洪水淹及。

三塘乡因地处康山蓄滞洪区南部山区，离溃口较远，且地势较高，溃堤后23 h 40 min，洪水才达到三塘乡马山村，最大淹没水深3.67 m，最大流速为0.2635 m/s，峰现时间为溃堤后29 h。魏家、余家桥村未被洪水淹及，该未被洪水淹及的区域可作为避险转移的重点区域。

以上6个区域中康山垦殖场片区最先被淹没，其中插旗分场在溃堤后30 min开始淹没，里溪分场、甘泉洲分场在溃堤后40 min开始淹没。康山垦殖场片区最大淹没水深和最大流速均出现在里溪分场，最大淹没水深为8.32 m，最大流速为0.074 m/s，峰现时间为溃堤后2 h 30 min。

此外由图6可知，康山垦殖场片区最先被淹没，几乎都被淹没，且淹没水深平均超过2 m，故该片区人口物资需要在预警时间内迅速撤离至离该片区较近的石口镇古竹片区安全楼或大塘乡片区安全楼；古竹片区的安全楼可建立在刘埠村和湖滨村，这两个村未被淹及；大塘乡片区的安全楼可建立在陈家塘村和南垅村；瑞洪片区的把下村也未被淹没，可将该区域的西岗村和罗家村的人口和物资迁至把下村避险；康山乡片区的大山村和金山村从蓄洪开始到结束均未被淹没，可在此处建立安全楼以便康山乡片区人口物资避险；三塘片区最后被淹及，且仅有马山村被淹没，其余地区未遭到淹没，故该片区人口物资相对安全。

4.7 淹没区损失统计

基于上述获得的淹没范围、深度、流速、预警时间等信息采用公式(4)和(5)估计淹没区生命损失，步骤如下：(1)通过淹没信息划分灾区，统计计算各灾区的风险人口。(2)通过洪水历时、水深、淹没范围来查表确定死亡率，将死亡率与溃堤概率相乘得到损失率。(3)将步骤(1)和(2)计算结果相乘得到生命损失。获得的淹没区风险人口累积曲线和不同预警时间内的人口损失如图7(a)～7(c)所示。

图7 溃堤造成的研究区内生命、经济和环境损失图

经济损失估算步骤如下：(1)在经济数据库中查询康山蓄滞洪区的行政面积，求出模拟单元格中的受灾单元格相对于行政面积的比例；各类财产货币化后的价值与比例的乘积即为该网格单元内各类型的财产价值。(2)由洪水演进模拟得到的淹没信息(水深、流速和淹没历时)确定单元内各类财产的损失率。(3)将步骤(1)得到的财产价值与步骤(2)得到的损失率相乘，得到该单元格的财产损失，累积计算淹没区域财产损失即为该蓄滞洪区的经济损失。根据统计蓄滞洪区淹没范围内的直接经济组成成分，采用公式(6)计算得到区内直接经济损失和农林牧及家庭财产损失累积曲线如图7(d)～7(e)所示。淹没范围内所有单元的损失累加，得到直接经济损失为9.53×108元，间接损失取直接损失的63%，为6.00×108元，则康山大堤溃决总经济损失为15.53×108元。

环境损失通过问卷调查的方式评估，采用公式(7)和(8)计算得到湖区重点堤防整治工程居民的平均支付意愿为每人776.8～780元/a。目前康山蓄滞洪区共有居民29 895人，按人均776.8元/a，按20 a计算，堤防环境整治工程总价值评估为4.64×108元，即堤防溃决导致的社会环境损失为4.64×108元。进而按照洪水淹没进程得到环境损失累积曲线，如图7(f)所示。

基于蓄滞洪区洪水演进过程数值模拟得到的洪水淹没信息，采用第3节生命、经济和环境损失评估方法可得到风险人口、淹没-损失人口、经济损失和环境损失，据此可根据不同预警时间下的淹没损失制定相应的应急预案和具体的抗洪抢险措施。比如由图7(a)可知，该蓄滞洪区内的风险人口随着洪水历时的一个变化过程，可结合淹没范围信息制定准确详细的人口迁移方案，以减少淹没损失和降低抗洪成本。图7(b)～7(c)除了评估人口损失之外，还可根据不同预警时间编制抢险方案，制定预警级别和相应的响应措施。根据图7(e)可判断直接损失中损失比值最大的成分为农业经济损失，从而为提前采取相关止损措施提供参考依据。

5 结 论

本文发展了堤防溃决洪水演进数值模拟方法，依托鄱阳湖康山蓄滞洪区，采用MIKE21建立洪水演算数值模型模拟蓄滞洪区洪水演进过程，获得了最高水位条件下的淹没范围、洪水淹没深度及流速过程线，进而估算了堤防溃决造成的蓄滞洪区内的生命、经济和生态环境损失。主要结论如下：

(1)获得了康山蓄滞洪区历史最高水位条件下的溃堤淹没信息(包括淹没范围、水深和流速等)，并根据洪水淹没范围、水深划分了区内灾害等级，最后根据溃决洪水到达特征点所需的时间、特征点最大淹没水深、流速峰值和峰现时间确定了各片区安全楼的位置，以便制订人口和物资迁移方案，这些研究结果可以有效指导蓄滞洪区的防洪抢险和人口转移工作。

(2)建立了基于蓄滞洪区洪水演进淹没数据的损失评估方法，基于溃堤淹没数据得到了蓄滞洪区内各片区的三大损失累积曲线。并根据获得的风险人口、淹没区人口损失、经济损失和环境损失判断了风险严重程度，从而为堤防溃决风险评估、康山及类似蓄滞洪区防汛抢险决策方案的制定提供了重要的参考依据。