黄河水库下游河段洪水影响规律数值模拟研究

孙继鑫，王 剑，杨少雄，杨波，陈光照，侯精明

(1.西安航天天绘数据技术有限公司,西安 710100;2.航天恒星科技有限公司，北京 100094；3.西安理工大学省部共建西北旱区生态水利国家重点实验室,西安 710048)

0 前 言

洪水是破坏性最强的自然灾害之一，严重危害人民生命财产安全和道路、桥梁等基础设施[1-3]。统计数据表明洪涝灾害已成为危及中国三分之二地区的重大灾害，造成的年均经济损失约占全国经济总产值的3.15%[4]。在防洪实践工作中，仅利用工程措施并不能完全抵御洪水威胁，与非工程措施相结合是防洪减灾工作的有效方法之一[5-6]。随着气候变化，近年来中小流域强降雨事件时有发生，形成较大洪水，尤其是水库区域，更可能诱发溃坝，造成严重的洪涝灾害。作为重要的非工程措施之一，提前制定避洪方案能使受灾区域更为科学地部署防汛抢险，组织群众应急避洪转移，提高避洪和应对洪水灾害的综合能力[7-8]。

近年来，国内外学者在洪水数值模拟方法和模型应用方面开展了深入探索，目前洪水分析常采用方法有水文学法、历史水灾法和水力学法[9-10]。水文学法是通过洪水频率分析，在获取不同洪水频率水位、流量和河段洪量的基础上，依据区域地形地貌分析流域产汇流特征，采用水文学模型计算得到不同频率洪水的淹没范围和淹没水深分布图。刘登嵩等[11]基于新安江等多种水文模型对浙江金华江和龙泉溪流域进行洪水预报与演进过程模拟；郝晓博等[12]基于NAM水文模型对福州市降雨洪水演进过程进行模拟与分析；李瑞等[13]基于河北雨洪模型对黄河上游小流域1994—2010年的10场洪水过程进行模拟,并分别以峰值加权均方根误差、加权误差及效率系数为目标函数，研究了不同目标函数对模型模拟结果的影响。历史水灾法是依据历史洪水系列计算相应的洪水频率和洪量等，在此基础上提取典型场次洪水的洪水淹没范围和淹没水深等，通常用于因洪致涝等水灾成因比较复杂的区域的洪灾风险评估，具有客观、简便和实用等特点。由于历史防洪形势与当前防洪形势通常存在一定差异，所得结果通常需要通过水文学和水力学模拟分析进行完善，难以独立地完成洪灾风险评估[14]。在洪水分析工作中，详细的历史洪水数据常常出现缺失乃至无法获取的情况,而水动力模拟方法由于其需要的参数、基础资料较少，且能够模拟水力要素在复杂地形上的变化过程，生成准确的洪水演进变化计算结果，近年来受到越来越多研究者的关注[15]。目前已开发了诸多商业软件及自主开发模型，这些模型在国内外进行了许多实例应用[16-20]，如张兆安等[21]采用耦合溃口演变过程的水动力数值模型,模拟分析了不同淤积程度下淤地坝溃坝洪水演进过程；朱世云等[22]采用MIKE21 FM模型进行水动力数值模拟对岳阳市洪水演进不同时段的淹没情况、洪水流态和进洪量进行了模拟；程坤等[23]采用MIKE21模型对藏木水电站大坝溃坝洪水过程进行了模拟，计算得到下游河道关键洪水要素，为评估河道下游两岸的受灾程度做基础；付成威等[24]基于实时动态耦合的水力学一、二维模型模拟了谷堆圩蓄滞洪区的溃堤洪水演进情况，得到了一二维耦合模型相较于传统模型精度及计算效率更高的结论；戴荣等[25]将HEC-HMS模型应用于新疆天山北坡地区各流域的洪水计算，证明了该模型适用于新疆天山北坡地区洪水计算；贺娟等[26]通过HEC-GeoRAS分析长河坝电站水库大坝下游溃坝洪水演进过程洪水淹没范围及流速分布。以上研究多集中于对河道洪水演进过程的模拟，缺少对沿程水力要素变化规律的量化分析及对避洪决策的研究。

近年来，吉林市黄河水库进行了规划整治，并建立了相对完善的防洪体系。但是依然存在着部分工程老化损坏严重，水库下游饮马河河道行洪能力不足、部分河道冲刷严重、防洪标准降低等问题。本文以水力学方法为理论基础，搭建黄河水库一、二维耦合洪水数值模型，根据当地潜在的洪涝风险因素情况设置多套洪水情景模拟方案，并结合社会经济数据分析灾害影响并确立避洪减灾方案，旨在增强当地人民群众及决策者的避灾减灾能力，为洪灾损失评估和避险转移策略制定提供技术支撑。

1 研究方法

洪水风险分析模型由水动力数值模型、大坝溃决分析模型、灾害影响评估分析模型、避洪转移决策模型构成。一、二维水动力数值模型用于洪水演进过程计算，确定洪水淹没要素(淹没范围、淹没水深和淹没历时等)，灾害影响评估分析模型及避洪转移决策模型用于避灾方案制定。

1.1 一、二维耦合水动力数值模型

本文所用水动力学模型为HydroMPM-FloodRisk模型，该模型可实现明渠、排水河道、水工构筑物及二维坡面流的模拟，主要应用于同时出现一、二维流动特征的数值模拟，实现一维河道计算和平面二维流动计算的无缝连接。河道非恒定流的水动力学模拟基于圣维南方程，方程守恒形式如下：

(1)

(2)

公式(1)～(2)中：q为侧向入流，m3/s；Q为流量，m3/s；s为距离坐标，m；V为断面平均流速，m/s；h为水深，m；A为过水断面面积，m2；B为河道断面宽度，m；i为渠底坡降；R为水力半径。

二维水动力计算模型基本原理如下：

(3)

(4)

公式(3)～(4)中：H为水深，m；Z为水位，m；q为连续方程中的源汇项；M与N分别为x和y方向的垂向平均单宽流量，m2/s；u和v分别为垂向平均流速在x与y方向的分量，m/s；n为曼宁系数；g为重力加速度，m/s2。

一维模型采用Abbott六点隐式格式离散控制方程组，该离散格式在每一个网格节点按顺序交替计算水位和流量。因此能够在较大库伦数的条件下保证计算稳定，实现大时间步长计算以节省计算时间。二维模型采用有限体积法离散求解模型控制方程，保证了水量和动量在计算域内守恒。并使用三角形非结构网格单元离散计算域，更有利于拟合复杂边界线，如堤防、桥、闸等；应用干湿网格判别法，能够高效高精度地处理潮滩动边界。本文在计算时，以设计流量过程作为输入条件，设定库区及下游无降雨过程，故在建模时并未耦合产汇流模型，在计算时河道洪水采用一维水力学方法进行洪水分析计算；淹没区洪水演进过程采用二维水动力学方法，模拟已设定计算方案对应的洪水演进过程，同时建立一、二维洪水耦合计算模型实现河道和淹没区的洪水演进过程计算。

1.2 避洪转移决策模型

按照《避洪转移图编制技术要求(试行)》[27]规范要求，选取各方案最大量级洪水作为避洪转移分析对象，分析洪水计算结果，提取满溢分流形成的最大淹没水深及淹没范围，绘制包络图，以形成包络图的淹没水深、淹没范围、洪水流速作为洪水风险图层，绘制基本洪水风险图。在基本风险图基础上，增加危险区及避洪转移范围图层、转移单元图层、安置场所图层、避洪转移路线图层，绘制形成避洪转移图。

2 模型构建及参数选择

2.1 研究区概况

黄河水库位于第二松花江流域饮马河水系，距烟筒山承德村李塘房屯0.5 km，流域面积784 km2，河道长度62.5 km，河道平均比降1.43‰，黄河水库是磐石市境内规模最大且集灌溉、防洪、发电、旅游等功能于一体，各项设施基本配套的中型水库。库区下游区间包括黄榆乡、山河街道、烟筒山镇、齐家镇和金家满族乡等，人口、经济较为集中。

2.2 模型构建

黄河水库下游洪水风险分析区至金家乡芦家村长岭水文站处，河段一、二维耦合计算模型范围包括黄河水库～长岭水文站河道两岸洪水影响区域，河段控制点位于长岭水文站，采用最大量级设计洪水试算得到最大淹没范围，在此淹没边界基础上适当外拓并参照地形等高线划定二维计算区，区域面积142.83 km2。分别对该段河道和两岸区域建立一维水动力模型和二维水动力模型，同时实现两者的动态耦合。

2.2.1一维河道水动力模型构建

根据饮马河地形地貌特征，一维模型构建的河段(黄河水库坝下至长岭水文站)全长53.6 km，共设置21个关键断面，断面间距变化范围为100～1 500 m。河道断面和河道断面布置分别如图1、2所示。

一维非恒定流模型的边界条件包括上及下边界条件。上边界条件选用流量过程，下边界条件选用水位过程或水位流量关系曲线，本次计算一维模型边界如表1所示。参考《吉林省磐石市饮马河重点段治理工程初步设计报告》[28]内容论述，曼宁系数设定为0.033。综合黄河水库下游实际情况及现场查勘结果，考虑模型稳定及运算效率等多种因素，设定黄河水库下游(黄河水库～金家乡芦家村)河道一维水动力模型计算初始水深为0.2 m。

表1 一维河道洪水计算边界

2.2.2二维水动力模型构建

模型采用非结构网格技术对模拟范围进行网格划分，最大网格面积不超过0.003 km2，重要地区、地形变化较大地区的计算网格适当加密。网格总数54 980个，计算总面积218 km2。在本次洪水分析过程中，二维水流模型的边界条件分为两类。第一类包括与一维模型耦合处的边界，具体包括与黄河水库下游河道的侧向连接处边界，此类边界均为动水位边界，由模型自动耦合计算；第二类是二维模型模拟区域外边界，由于在确定建模范围时已考虑了河道洪水边界，模型计算范围内区域与区域外不存在水量交换，因此确定为固边界。考虑到编制区域内地物种类的不同且分布零散，为保证二维模型计算精度，曼宁系数按编制区域土地利用情况进行分区，并参照其他类似区域糙率设定，无实测资料的地区可根据水力学手册确定，或参考采用相似条件地区的糙率。本地区缺少溃堤洪水演进实测资料，故按照《洪水风险图编制技术细则(试行)》及《水力计算手册(第二版)》[29-30]所列的合理范围内取值，如表2所示。

表2 洪水风险区域曼宁系数取值

2.2.3一、二维水动力耦合模型构建

采用侧向连接的耦合方式，实现河道一维水动力模型与二维水动力模型的耦合，实时耦合计算河道洪水漫溢淹没风险。侧向连接方式为通过河道断面标注堤顶等效为堰，堰顶高程及堰宽以该处断面左右堤顶高程及宽度为准；沿一维河道边界线在二维区域设定耦合线，耦合线由具体坐标确定，在运行期，耦合线映射到网格单元的边上，确定耦合的网格单元。模型通过比较二维网格与相应里程处一维河道内断面水深，利用堰流公式计算通过侧向连接的流量，实现一维河道堤顶过流水流漫溢进入二维平面区域和一二维水流的动态交互。耦合模型计算时间步长设定为较小值60 s。

3 模型验证

本文所用水动力学模型为HydroMPM-FloodRisk模型，该模型2014年被国家防汛抗旱总指挥部办公室批准列入重点地区洪水风险图编制项目软件名录，模型精度较高，计算结果的准确性得到了广泛的认可和验证[31-33]。

为保证一二维耦合水动力数值模型在本次模拟中结果的可靠性，从水量平衡、流场分布、淹没水深、淹没面积进行对比分析。对50年一遇洪水模拟计算方案进行了水量平衡统计分析，表3列出了洪水模拟计算方案模拟得到的保护区初始水量、进洪量、出水量以及最终蓄水量。水量平衡误差为-0.48×106m3，均低于106数量级，满足水量平衡要求。

表3 水量平衡对比 /(×106 m3)

对比图3所示的DEM整体高程与淹没区域流场分布模拟结果，可以明显看出流场分布均匀一致，流速较大的区域集中在坡度变化剧烈，高低落差较大的地形区，洪水演进的趋势遵循地形高程从高到低的原则，洪水态势较准确，洪水模拟流速合理，淹没区域与DEM地形所示的地势低洼区相符。模拟结果表明模型计算较为合理，且精度较高，可应用于洪水风险计算。

4 黄河水库下游河段洪水研究规律

黄河水库下游(黄河水库～金家乡芦家村)河道洪水来源主要有上游来洪、黄河水库至金家乡芦家村区间洪水、溃坝洪水。黄河水库下游现有堤防1处，堤防位于黄榆乡平埠子村-茶壶嘴村饮马河河段，防洪标准为20年一遇，堤防长度为10.5 km，河道安全泄洪量为820 m3/s。由于堤坊工程缺乏维护，局部堤段存在隐患，致使防洪标准降低。为了量化黄河水库下游水力要素与来流量的规律，分别对10、20、50年一遇上游来洪进行模拟分析，并对下游关键断面处最高及最低水位进行统计整理，并基于回归分析法法对其进行拟合。

4.1 设计洪水方案

采用SL 44-2006《水利水电工程设计洪水计算规范》[34]中所列的Pearson-Ⅲ曲线进行频率计算，利用矩法初步估算水文连序系列和不连序系列的均值、Cv、Cs等统计参数值，频率分析成果见表4。

表4 长岭站设计洪水成果

根据洪水量级分析，推求出长岭水文站的10、20、 50年一遇的设计洪水过程。依据SL 44-2006《水利水电工程设计洪水计算规范》，设计洪水采用放大典型洪水过程线的方法，并选择能够反映洪水特性、对工程防洪运用较不利的大洪水作为典型。长岭水文站的2017年7月19日至7月24日的实际洪水过程峰高量大，能够反映黄河水库下游大洪水的洪水特性。因此，选择此次洪水作为典型洪水，流量放大后的设计洪水过程线见图4。

4.2 洪水过程水力要素模拟计算

4.2.1不同洪水条件下最大淹没水深模拟结果分析

基于HydroMPM-FloodRisk水力学模型对黄河水库下游不同重现期来洪条件下洪水演进过程进行模拟，在模拟时，对研究区域内高于0.5 m的道路等线状地物进行适当概化。当洪水未达到阻水建筑物顶高程时，线状地物起阻水作用，区域不过水；当洪水超过线状地物顶部高程时，水流以漫溢的形式通过。模型计算时间步长为60 s。根据模拟结果，得到10、20、50年一遇来洪条件下黄河水库下游最大淹没水深如图5所示，图中以色带的颜色深浅表征洪水水深程度。

由图5可知，随着上游来洪重现期的增加，在研究区域内，淹没面积和最大淹没水深均增加。在各重现期条件下，均产生不同程度的漫溢，但随着重现起的增大，漫溢的面积和水深逐渐增大，为了量化不同重现期条件下研究区域水位的变化规律，需提取重要断面的水力要素，并对其进行分析。

4.2.2不同洪水条件下沿程断面水力要素结果分析

模拟黄河水库下游段10、20、50年一遇洪水演进过程，根据河道上断面入流流量过程，结合各断面水力参数，采用非恒定流水动力计算方法推算各断面水位及流量过程，统计计算结果，提取河道计算断面最高、最低水位。如表5所示。绘制河道最高水面线，如图6所示。

表5 不同重现期洪水下游断面水位

基于黄河水库下游各监测断面在10、20、50年一遇洪水工况下的水位数据，拟合以里程为自变量，断面水位为因变量的线性方程，如图7所示。可以通过线性方程获知在10、20、30年一遇工况下不同里程的河道断面水位变化情况，从而实现通过少量断面数据获知河道沿程水位变化情况，极大地缩短了研究人员的数据处理和分析决策时间。

4.3 避洪转移决策分析

基于研究区50年一遇洪水演进模拟结果，分析淹没水深风险信息，提取最大淹没水深和最短到达时间数据，绘制最大水深最短时间范围包络图，确定淹没范围即为危险区，淹没区内居民地将被定为转移单元。当黄河水库下游计算区遭遇50年一遇洪水时，黄河水库至金家乡芦家村发生河道洪水漫溢，不同淹没水深等级全区资产影响统计见表6。

表6 50年一遇洪水损失统计结果

采用GIS技术开展空间分析，根据转移单元人数、安置区容纳能力、距离远近、行政区界等因素，确定转移单元和安置区的对应关系。在尽可能避免交通(特别是交汇路口)拥堵、避免跨县级以上行政区安置，且各安置人口密度相当的情况，按照就近原则，结合道路分布情况，确定转移方向及路径，建立转移单元与安置区的对应关系，避洪转移方案如图8所示。

5 结 论

针对中小流域洪水分析资料匮乏，本文采用了一、二维耦合水动力模型精细化模拟黄河水库下游洪水演进过程。结合社会经济基础数据，进行了灾情分析和损失评估并制定避洪转移方案。文中综合考虑了黄河水库洪源、防洪工程以及河道防洪标准，以不同重现期设计洪水方案为例，对水库下游建立一、二维耦合水动力模型模拟河流洪水漫溢分流演进过程，并提取了关键断面水力要素值进行量化分析，且基于标准洪水风险划分及避洪转移要求，绘制了50年一遇工况条件下研究区域避洪转移图，形成结论如下：

(1) 中小流域普遍存在径流观测资料缺失或不足的情况，而本次采用一、二维耦合水动力模型能够良好地避免受到制约。通过对模拟结果分析，精度满足要求，具有一定的可靠性，可为其他缺乏径流观测资料中小流域的洪水风险分析及避洪决策制定提供参考。

(2) 通过对不同洪水重现期条件的洪水演进过程模拟，随着洪水重现期的增大，淹没水深及淹没面积逐渐增大，且在各重现期洪水条件下，河道两岸沿程多处发生河道漫溢现象。

(3) 通过对不同洪水重现期条件下断面水力要素的提取与分析，发现河道沿程水位呈现线性变化规律，且随着洪水重现期的增大，沿程最低与最高水位均逐渐增加。如在水库下游4 000 m处，当洪水重现期为10年一遇时，最低水位为252.84 m，最高水位为255.43 m，当洪水重现期增大至50年一遇时，最低水位增高至253.34 m，最高水位升高至256.89 m。

(4) 根据洪水模拟计算结果，进行洪灾影响分析评估，确定了防洪保护区在遭受上述洪水时需转移单元及人口数量、规划安置场所、制定转移路线，避洪方案的合理性分析。

本文提出的水库洪水风险分析及避洪决策方法能够较好地预测洪水演进过程，可据此制定正确的避洪转移方案。未来研究中将完善该洪水风险分析方法，开展实时洪水风险计算分析，结合水库调度及风险管理开发基于动态洪水风险图的防汛决策系统，同时，还可具有灾情评估分析、避洪转移分析等功能，所有结果在系统平台上实时表现，应用灵活，效果直观。