0 引言

随着城市化进程不断加快，气候条件及下垫面条件不断变化，造成城市地区洪水峰高量大，洪涝灾害更为频繁，影响范围也愈加广泛^[1]。我国目前正面临着从“控制洪水”向“洪水管理”的转变^[2]，采取工程措施和非工程措施相结合的方式防范洪涝灾害已经成为我国防洪排涝事业发展的必然趋势。作为一项典型的非工程措施，洪水风险分析以信息化、自动化为基础，将地理信息和防汛信息综合考虑，在预测洪水演进路线、避险转移、土地利用、保护人民生命财产安全方面具有重要的技术支撑作用。

近年来，防汛信息化建设尤其受到国内外学者关注。尤其是洪水风险图系统开发与应用、水动力学模型及洪水风险分析等方面^[3-5]，国内外很多学者开展了相关研究。例如：Horritt等^[6]利用HEC-RAS和Telemac-2D模型分别进行了洪水淹没模拟计算，对比了两者的精度，结果表明，对于该河段，HEC-RAS和Telemac-2D模型都可以根据流量或淹没区数据进行校准，并对淹没区进行良好的预测；Liang等^[7]采用动态链接库技术(DLL)将一、二维数值模型相结合，对浅水流进行了模拟，并且对二维模型进行了改进，改进后的模型能够较好地反映地形、洪水、干旱过程；Patro等^[8]采用一维-二维耦合水动力模型Mike 11对印度马哈纳迪河流域三角洲地区的洪水淹没程度和淹没深度进行了模拟，通过验证，模型精度良好；Dhruvesh等^[9]对印度苏拉特市进行了基于HEC-RAS的一维/二维耦合水动力模拟，通过实例研究证明了一维/二维耦合水动力模型在洪水淹没测绘中的适用性良好；张大伟等^[10]建立了基于天然河道溃堤模拟的水动力学模型，通过典型算例验证了其合理性，并将其应用于松花江哈尔滨段；陈文龙等^[11]建立了基于侧向联解的一维-二维耦合水动力学数值模型，通过大量算例验证得出，模型能够有效地模拟出任意形状溃口的溃堤洪水演进过程；苑希民等^[12]运用干湿水深理论对一、二维耦合水动力学数值模型进行优化，并将其应用于黄河青铜峡河西灌区，较为真实地再现了洪水在计算区域内的演进过程与淹没范围。由于计算机计算速度及基础资料缺失的限制，早期的研究模型较为简单，耦合形式松散。随着计算机技术的不断发展，现有模型已经能够接近还原真实水文循环过程^[13]。

本文构建了能够反映真实地形的浙江省台州市灵江下游洪水风险分析模型，通过水量交互对其进行耦合，考虑线状地物的阻水作用，通过干湿水深和糙率分区对其进行优化，并结合社会经济资料和地形数据对其进行洪水影响分析，以期对灵江流域防洪避险、分洪利用及洪水调度等提供参考依据。

1 水动力学模型方法及其耦合原理 1.1 一维水动力学模型

灵江流域位于浙江省台州市。灵江属于独流入海河流，江水注入东海。流域地形复杂、河网密布，水系间情况复杂、相关性大。一维水动力模型可以较为准确地对流域内河网水动力过程进行快速模拟，因此本研究采用一维水动力模型。一维水流运动由连续性方程与动量方程组成的Saint-Venant方程组进行描述：

式中：A为河道过水断面面积，m²；t为过水时间，s；Q为河道流量，m³/s；x为沿流向的水平坐标；q_i为外侧入流流量，m³/s；α为断面动量修正系数；g为重力加速度，取9.81m/s²；y为断面水位高度，m；J_f为摩阻坡降；u为外侧入流在河道方向上的流速，m/s。

在河道交汇处各河段依据水量平衡关系进行连接:

式中：Q_b^c+1为第c+l时段流入节点b的外加流量，m³/s；L(b)为与河道交汇节点b相连接的河段总数；B为节点总数；Q_{b, j}^c+1为第c+l时段河段j流入节点b的流量，m³/s；ΔV为河道交汇处蓄水量变化量，m³。

河道恒定流的模拟可釆用曼宁公式：

式中：n为河道糙率；R为河道水力半径，m；i为底坡坡度。

本文对上述方程离散采用Preissman四点隐式差分格式。该格式在运算过程中具有计算精确、过程稳定、适应性强等优点，有利于提高计算复杂地形情况下的模型稳定性^[14]。

1.2 二维水动力学模型

河道发生溃决后，水流漫溢进入防洪保护区中，此时水流呈现明显的二维特性，水流流态更为复杂^[15]。此时，一维模型不再适用于模拟漫溢后进入防洪保护区的水流。二维水动力学模型能够比一维模型提供更加丰富的水力要素信息，能更直观地显示洪水演进过程^[16]。在缺乏动态洪水数据时，可以通过一维-二维耦合模型对一维模型进行率定^[17]。在现代计算机技术与计算方法的快速发展下，二维水动力模型是常用的分析技术方法之一，越来越多的水利工程采用二维水动力模型对任意组合、不同规模的洪水进行洪水淹没模拟。

二维水流运动同样也通过连续性方程与动量方程进行描述：

式中：h为断面水深，m；E，F分别为x，y方向上的单宽流量，m/s；S₀为源项；m；J_fx，J_fy分别为x，y方向摩阻坡降；i_0x ，i_0y分别为x，y方向底坡坡度。

科氏力和紊动项的影响在本文中影响较小，故未进行分析。为使水量守恒条件在计算域中得到满足，本文运用守恒格式对水流连续方程进行离散，时间、空间上均采用向中间差分的蛙跳格式，并且使用交错网格对物理量进行计算^[18]。

溃堤水流具有水面梯度大、水面存在间断等特点，为了能够准确模拟溃堤水流的复杂流态，本文采用有限体积法进行非结构化网格离散求解。本文采用近似Roe格式的Riemann解对界面通量进行求解，其对激波有较好的捕捉能力，应用十分广泛；采用特征分解法对底坡源项进行数值离散，以保证计算格式的和谐性；采用半隐式格式对摩阻源项进行数值离散，以保证计算格式的稳定性^[10]。

1.3 一维-二维模型耦合

本文一维河网模型与二维浅水模型的耦合通过溃口处的水力连接条件来实现。在溃口位置用堰流公式计算溃口处的流量与水位，通过水量平衡关系来实现一维-二维模型的耦合与衔接。模型耦合的关键在于河道与保护区的水力信息交互。一维河网模型计算得到断面的水力学参数，二维水动力学模型计算得到网格平均水力学参数。溃口处的网格边界与河道断面相衔接。一维河网模型为二维防洪保护区模型提供流量值Q作为边界条件，在二维计算区域水流入口处设置流量值Q，二维模型网格中水位互不相等，将二维模型计算获得的网格平均水位值Z₂反馈给一维模型, 以此来进行下一时段的计算，通过这种方式来实现一维-二维模型的耦合^[19]。具体耦合原理见图 1。

2 模型构建与验证 2.1 研究区域概况

本文研究区域灵江下游位于浙江省东南部，地势自西北向东南倾斜，毗邻东海，历史上洪灾频发，对两岸居民生命财产安全危害巨大。研究区域上游受始丰溪、永安溪来水影响，义城港支流汇入研究区域中游，同时下游受东海高潮位顶托。夏季常受热带风暴和台风影响，突降暴雨，形成洪灾。研究区域位置和地形如图 2所示。

2.2 模型构建 2.2.1 一维模型构建

根据流域河道地形、水文气象等资料，构建灵江下游一维河网模型，进行洪水演进过程计算。本文对研究区范围内主要河道进行建模，涉及河道包括始丰溪、永安溪、义城港、灵江，其余较小支流由于缺少河道断面资料等原因，采用集中入流方式汇入，即将支流流量加到主河道对应节点上。一维河道概化断面间距在200~2 000 m之间，共184个断面，概化河道总长约176 km。研究区一维河网模型上边界采用永安溪罗渡站、始丰溪百步站两个水文站作为流量边界，下游采用灵江入海口附近海门潮位站作为水位边界。研究区一维河网概化模型如图 3所示。

2.2.2 二维模型构建

使用地理信息系统软件勾勒出研究区轮廓，初步计算研究面积，划分计算网格。通过大量级洪水多次试算，划定可能淹没的最大范围作为研究区域，进行网格剖分。选取计算网格时，选取易于获得高质量网格、易于适应复杂地形的非结构化网格。本文考虑了线性地物(如高速公路、铁路、省道、堤防等)的阻水和导水作用，将堤防、研究范围作为外部约束条件，铁路、高速路等作为内部约束条件，并将其周围网格进行适当加密，以充分考虑其对水流运动和洪水演进过程的影响。在人口密集、经济发达的城市中心地区，为了能够准确分析洪水淹没范围和淹没水深，确定洪水损失，网格也进行适当加密，网格边长为100 m。最终确定计算方案的网格计算单元数为10 515，结点数为9 239，平均边长为123 m，网格平均面积为0.004 km²，是比较合理的网格划分结果。根据数字高程模型和高精度的遥感影像给网格赋高程值。进行高程插值后的研究区二维计算网格如图 4所示。

2.3 模型参数选取与率定

永安溪、始丰溪均为山溪性河流，河底多为砂石、卵石，根据水力学相关参数取值规范^[20]，根据高分辨率遥感图结合河道断面实际地形情况，河道糙率选取0.030~0.065，义城港及灵江均为平原河网，河道糙率选取0.028~0.035。二维网格利用高清遥感影像分辨下垫面土地利用类型进行赋值。

选取2015年“苏迪罗”洪水进行模型率定。进行模型参数率定时，模型输入“苏迪罗”洪水发生时各边界水位流量条件，计算出西门站水位过程与实测值。模型计算起始时刻为20150809T0:00，计算结束时刻为20150811T20:00。模型计算值与实测值对比如图 5所示，可以看出模型计算参数能够较为准确地反映洪水演进过程。

2.4 模型验证

模型验证选取2007年“韦帕”洪水。在以2007年“韦帕”洪水进行率定和验证时，模型输入洪水期间罗渡站流量、百步站的流量和海门站的实测潮位，将计算出的各典型断面水位与实测值做比较，来推断模型参数的合理性。通过水文年鉴搜集柏枝岙站、沙段站及各海门实测资料，由柏枝岙站和沙段站流量资料根据比拟法推求出罗渡、百步站的流量。本研究中将柏芝岙站和沙段站作为参证站，根据式(8)推求出罗渡站和百步站流量：

式中：Q_s为设计站流量，m³/s；Q_c为参证站流量，m³/s；F_s为设计站面积，km²；F_c为参证站面积，km²；m为流域面积指数。

将罗渡、百步站的流量作为一维-二维耦合模型上边界条件，海门站潮位资料作为该模型下边界条件，用各典型断面实测洪水作为验证。模型计算结果与实测水位对比关系如表 1所示。

由表 1可见，永安溪(柏枝岙)、始丰溪(沙段)和灵江沿程(其余6个站点)特征点验证计算水位与实测水位相差幅度均在0.10 m以内，说明模型模拟情况较好，边界条件、参数选取以及断面概化均十分合理，可以作为规划方案的计算和分析。

3 耦合模型计算结果分析 3.1 合理性分析

各频率设计洪水条件下计算得到的典型断面水位见表 2。选取模型计算水位与流域相关规划或工程设计报告进行了对比分析，计算结果见表 2。由表 2看出，在频率分别为20%，10%，5%，2%，1%下：永安溪柏枝岙站断面水位计算结果与永安溪治理可研报告计算结果偏差分别为0.46%、-0.24 %、-0.69 %、-1.10 %、-1.53 %，均不超过±5.00%，偏差范围不大，断面计算结果合理；始丰溪沙段站断面水位计算结果与始丰溪综合治理可研报告计算结果偏差分别为-0.12 %、-1.56 %、0.16 %、0.73 %、-0.38 %，偏差范围不大，该断面计算结果合理；灵江西门站断面计算结果与椒江综合规划报告计算结果偏差分别为-3.34 %、-0.58 %、-2.52 %、-3.09 %、-2.41 %，偏差在合理范围内，该断面计算结果合理。

根据各频率典型断面水位模型计算结果与工程报告计算结果对比，偏差均在合理范围内(均不超过±5.00%)，可认为该模型反映各量级、各种组合洪水淹没情况是可信的，模型计算结果合理可靠，可应用于研究区洪水风险分析。

3.2 洪水演进过程分析

由于研究区水系沿岸堤防情况复杂，故本文选取100年一遇超标准设计洪水为例进行详细分析，溃口洪水方案选取永安溪10年一遇括苍镇马横溪段右岸溃口洪水进行分析。同时，分析5种不同重现期设计洪水的演进趋势及淹没水深情况。其中，不同重现期设计洪水只考虑漫溢，不考虑溃堤。其余洪水工况不再赘述。

3.2.1 灵江100年一遇设计洪水演进分析

灵江100年一遇设计洪水演进过程如图 6所示。

灵江遭遇100年一遇洪水时，对应西门站水位为11.60 m。如图 7所示：杜岐、花街、上百岩附近最先被淹没；然后，在洪水漫延6 h左右，就会影响后泾村、横山前村、上金、管岙、涌泉、汛桥、下洋峙村、石村、杨梅村、石鼓、梅浦村、沿岸；在24 h以后，会影响西洋、渡头范、两水村、小峧头、峰山头、官园、小溪、青潭头、车门桥、松山、下缸窑、白水洋、张家渡、应家山头、双楼、更楼、留贤、白毛、长潭、河头等。最大淹没水深达到8.43 m，最大淹没历时将会达到41.01 h，淹没面积92.16 km²，影响人口4.77万人。

3.2.2 永安溪10年一遇括苍镇马横溪段右岸溃口洪水分析

根据险工险段及历史溃口调查情况显示，括苍镇马横溪段右岸历史上曾经发生过河堤溃决，结合地形条件及工程建设情况，选择该处作为本次模型溃口进行分析。溃口具体位置见图 8。

永安溪如果遭遇10年一遇洪水，永安溪括苍镇马横溪段右岸河堤出现溃口。溃口宽度为10 m，溃堤底高程为12 m，溃堤水位为12 m，溃决方式采用瞬时全溃，溃口最大流量为58.84 m³/s，淹没面积2.50 km² ，影响人口0.07万人。

本次溃口洪水演进过程如图 9所示。

由图 10可知：在第30小时溃口处保护区水位开始上涨，此时溃口开始溃决; 到第35小时溃口处保护区水位达到14.5 m，流量达到58 m³/s; 之后随着河道水位降低，溃口处流量也逐渐减小; 在溃口溃决发生17 h后，即t = 47 h时，河道水位与保护区水位基本持平。保护区洪水回流至河道中。

3.2.3 各重现期洪水演进趋势对比

本文通过对不同频率水文资料的收集、计算，得到了重现期分别为5、10、20、50、100 a时的水文边界条件，输入模型中进行计算，得到了不同重现期下洪水演进范围，如图 11所示。

由图 11可以看出：遭受洪水侵袭时，相比灵江沿岸，始丰溪、永安溪沿岸淹没范围较小，这是由于始丰溪、永安溪两条河流均为山溪性河流，两岸地势较高，发生小量级洪水时，两岸基本不受影响；灵江下游城区防洪体系较为完善，堤防标准达到50年一遇，但模型计算显示，当城区遭遇100年一遇洪水时，城区也会受到影响。

4 结论及建议

1) 本文基于一维-二维耦合水动力学模型对灵江中下游典型城市进行洪水风险综合分析。通过构建研究区水系一维河网模型与防洪保护区二维模型及其耦合模型，对研究区5、10、20、50、100年一遇设计洪水和括苍镇马横溪段溃口洪水的演进过程进行模拟计算，利用历史洪水方案对其进行率定验证，保证了模型的计算效率和准确性。采用非结构化网格以提高模型对复杂地形的适应性。通过对道路、堤防等线状地物的概化，充分考虑了构筑物的导水及阻水作用，建立了能够反映真实地形的水动力学模型，提高了模型的合理性。

2) 为了能够快速准确地得到实时洪水风险要素淹没结果，可采用人工神经网络方法进行数据交换和一维-二维模型之间的连接，对模型进行进一步优化，以提高模型计算速度。

3) 高精度的基础地形资料和水文资料是洪水风险分析的基础，因此应当及时收集下垫面及防洪工程建设情况资料，定期对洪水风险分析计算结果进行复核，并加强防汛信息化建设。

4) 将防洪工程措施与非工程措施相结合，防汛部门应当洪水风险分析暴露出的现状防洪问题采取工程措施进行解决。永安溪、始丰溪均为山溪性河流，坡陡流急，河道比降大，两岸缺乏防护设施，河道冲刷严重。建议加快建设主要水系堤防达标工程，并对已建堤防进行堤基防渗及堤脚加固处理。