药湖联圩防洪保护区属于赣江水系，有石岗、松湖等多个重要乡镇以及昌樟公路等重要交通设施，土地肥沃，水土资源以及矿产资源十分丰富，是当地重要的商品粮生产基地。同时区内水文条件错综复杂，降水时空分布不均，夏多冬少；流域洪水主要由暴雨形成，洪水季节与暴雨季节相一致。加上药湖联圩水系繁多，圩堤设计标准偏低，保护区地形平坦，地势低洼，且同时受赣江和锦江洪水的双重影响，汛期时圩内险象环生，防洪任务异常繁重，不仅严重制约了当地社会经济的发展，同时对保护区内8万居民的生命财产安全也构成极大的威胁。而对于这样重要的区域，其洪水风险评估仍是一片空白，因此进行药湖联圩洪水风险分析，对于保护当地居民生命财产安全、促进保护区经济社会的可持续发展，具有十分重要的现实意义。

近年来，国内外学者对洪水数值模拟技术开展了广泛研究，并取得了较大进展。苑希民等^[1-3]基于二维浅水方程，建立了模拟河道和灌区洪水演进的二维水动力模型。在此基础上，李云等^[4]将一维、二维洪水嵌套模型成功应用于淮河临淮岗段。姜晓明等^[5]运用基于黎曼近似解的溃堤洪水一、二维耦合数学模型对松花江干流胖头泡溃堤洪水进行了模拟计算。张大伟等^[6]建立了能够对天然河道溃堤水流进行模拟的一维、二维耦合数学模型，并已成功应用于松花江哈尔滨段松北分洪区溃堤洪水的模拟。李大鸣等^[7]采用有限体积法建立了河道、滞洪区联合计算的一、二维衔接洪水演进数学模型及洪水仿真模拟系统，并对大清河滞洪区五洼联合调度方案进行了实时模拟。王炜^[8]将Delft 3D-FLOW应用于洪泛区的二维洪水演进数值模拟并与一维模型相衔接，对洪水造成的损失做了初步评估。可以发现，由于一、二维耦合水动力模型具有计算速度快、精度高的优点，已被广泛应用于溃堤洪水模拟中^[9]。

本文将建立一、二维耦合洪水演进模型，一维河道和二维洪泛区水流运动分别采用四点隐式Preissmann格式以及中心格式的有限体积法进行离散求解，模拟极端洪水过程获取相关洪水风险要素并进行洪水风险分析与评估。

1 模型构建 1.1 河道一维模型

一维河道(河网)的洪水运动用Saint-Venant方程组描述：

$\partial Q/\partial x + \partial A/\partial t = q$

(1)

$\frac{{\partial Q}}{{\partial t}} + \frac{\partial }{{\partial x}}\left( {\alpha \frac{{{Q^2}}}{A}} \right) + gA\frac{{\partial Z}}{{\partial x}} + gA\frac{{\left| Q \right|Q}}{{{K^2}}} = q{v_x}$

(2)

式中：Q为断面流量；x为河长；t为时间；A为过水面积；Z为断面水位；α为动量修正系数；K为流量模数；q为旁侧入流；v_x为入流沿水流方向的速度分量。其上、下游边界的控制条件一般采用水位过程控制、流量过程控制、流量与水位关系控制等形式。

对Saint-Venant方程组(1)和(2)采用四点隐式的Preissmann格式进行离散求解：

$\partial f/\partial t = \left( {f_{i + 1}^{n + 1} + f_i^{n + 1} - f_{i + 1}^n - f_i^n} \right)/\left( {2\Delta t} \right)$

(3)

$\partial f/\partial x = \theta \left( {f_{i + 1}^{n + 1} - f_i^{n + 1}} \right)/\Delta x + \left( {1 - \theta } \right)\left( {f_{i + 1}^n - f_i^n} \right)\Delta x$

(4)

$f = \left( {f_{i + 1}^{n + 1} + f_i^{n + 1} + f_{i + 1}^n + f_i^n} \right)/4$

(5)

式中：θ为权重系数，可取0.5＜θ＜1.0。

1.2 洪泛区二维模型

对于宽浅河道和洪泛区内大范围的溃堤洪水运动，可以采用浅水方程来描述，控制方程的守恒形式如下:

$\frac{{\partial \mathit{\boldsymbol{U}}}}{{\partial t}} + \nabla \cdot \mathit{\boldsymbol{F = }}\frac{{\partial \mathit{\boldsymbol{U}}}}{{\partial t}} + \frac{{\partial \mathit{\boldsymbol{E}}}}{{\partial x}} + \frac{{\partial \mathit{\boldsymbol{G}}}}{{\partial y}} = \mathit{\boldsymbol{S}}$

(6)

$\begin{array}{l} \mathit{\boldsymbol{U = }}{\left[ {h\quad hu\quad hv} \right]^{ - 1}};\mathit{\boldsymbol{F = }}Ei + Gj;\mathit{\boldsymbol{E = }}{\left[ {hu\quad h{u^2} + \frac{{g{h^2}}}{2}\quad huv} \right]^{ - 1}};\\ \mathit{\boldsymbol{G = }}{\left[ {hu\quad huv\quad h{u^2} + \frac{{g{h^2}}}{2}} \right]^{ - 1}};\\ \mathit{\boldsymbol{S = }}\left[ {0\quad gh({S_{ox}} - {S_{fx}})\quad gh({S_{oy}} - {S_{fy}})} \right];{S_{ox}} = - \frac{{\partial z}}{{\partial x}};\\ {S_{fx}} = - \frac{{\partial z}}{{\partial y}};{S_{fx}} = \frac{{{n^2}u\sqrt {{u^2} + {v^2}} }}{{{h^{\frac{4}{3}}}}};{S_{fy}} = \frac{{{n^2}v\sqrt {{u^2} + {v^2}} }}{{{h^{\frac{4}{3}}}}}。\end{array}$

式中：h为水深；t为时间变量；u，v分别为x，y方向的流速；S_ox，S_oy分别为x，y方向的坡底源项；S_fx，S_fy分别为x，y方向的摩阻项；g为重力加速度；z为底高程；n为曼宁系数。

本文采用了中心格式的有限体积法对控制方程式(6)进行离散，控制体取为任意的三角形单元，应用格林公式可得如下形式：

$\int_{{A_i}} {\frac{{\partial \mathit{\boldsymbol{U}}}}{{\partial x}}{\rm{d}}A + \int_{{\mathit{\Gamma }_\mathit{i}}} {\left( {\mathit{\boldsymbol{F}} \cdot \mathit{\boldsymbol{n}}} \right)} } {\rm{d}}\mathit{\Gamma = }\int_{{A_i}} {\mathit{\boldsymbol{S}}{\rm{d}}A} $

(7)

式中：A_i为i单元的计算面积；Γ_i为i控制体的边界条件；n为边界的外法线方向单位向量；F为界面通量。针对溃堤水流的流态复杂，水面梯度大等特点，利用Roe格式的近似Riemann解对界面通量进行数值求解^[10]。

1.3 模型耦合

将一维河网模型与二维洪水模型进行耦合，构建最终的洪水分析模型。选定侧堰流公式来实现溃口上下游水流信息的交互，具体为：

$Q = {C_{\rm{d}}}b\sqrt {2g} {H^{\frac{3}{2}}}$

(8)

式中：Q为侧堰流量；H为侧堰首端河渠断面水深；b为侧堰堰宽；C_d为流量系数。当建立溃堤水流模型时，可根据湍流方程推求流量系数的近似方程：

${C_{\rm{d}}} = \sqrt {\left( {1 - m} \right)} {d^{0.067}}/\left( {n\sqrt {{D_{\rm{s}}}} } \right)$

(9)

式中：m为模块限制系数；d为平均水深；D_s为水流源汇距离；n为曼宁系数。在溃口处二维计算单元一般通过多个网格点与一维计算单元连接。由于一维模型计算结果中的水力参数是物理量的断面平均值，二维模型计算出的是变量各网格中心处的节点值，因此在溃口连接处需要对一维、二维模型的交换数据进行转化和衔接。一维模型为二维模型提供流量值Q作为二维模型的边界条件，将Q值分布到二维计算单元的各节点上，考虑到在连接处二维计算网格的水位值并不相等，因此取各个计算网格的平均水位值Z返回给一维模型，以进行下一时段的计算。

2 药湖联圩区洪水风险分析 2.1 区域概况与地形概化

药湖联圩防洪保护区位于锦江下游南岸，丰城市与南昌市新建区交界处，面积为166 km²，地势总体西北高、南东低，内有锦江、药湖、浠湖、佳湖、上塘导托、圳头导托等多条水系，其中锦江主河长305 km，河面宽为200~500 m，流域面积383 km²，计算区域水系(不含导托)分布见图 1。药湖联圩受东亚季风影响，形成了亚热带季风气候，区域内热量丰富、雨水充沛，光照充足。区域雨洪同期，洪水最早出现在3—4月，但洪量值相对较小，5—6月进入主汛期，约占全年降水量的50%。

药湖联圩防洪保护区外洪模型主要考虑锦江溃堤洪水对保护区的影响，对锦江河道构建一维模型，洪泛区构建二维模型。计算区域采用非结构化三角形网格进行剖分，以计算域外边界、区域内堤防、阻水建筑物、较大河渠、主要公路、铁路作为依据，使网格的大小随地形地势和阻水建筑物的分布灵活确定，而且尽可能地将影响水流的阻水建筑物作为网格边界，充分反映计算域的特征。但是，保护区内部的道路、堤防需概化为线状阻水建筑物，内部河道则概化其河道堤防、下挖河底地形，并对网格适当加密，在二维地形中充分反映其特征。计算区域模型共生成计算网格约16万个，网格边长不超过300 m，最大三角形面积为20 000 m²。司家闸溃口附近部分网格概化如图 2所示。

2.2 计算条件与糙率确定

根据历史溃口及险工险段分析，考虑可能发生的最危险情况，通过实地查勘并结合地方防汛部门专家的建议，选择锦江堤段的鸡鸣洲(溃口1)、松湖堤段中段大王庙(溃口2)、药湖堤段司家闸(溃口3)为计算溃口，如图 1所示。堤防设计水位25.37 m，溃口宽度分别为170，110，190 m，溃口堤顶高程分别为30.0，28.5，28.4 m。锦江上游边界条件为高安站洪水流量过程，下游边界为设计条件下市汊相应洪水水位过程，利用锦江高安站1953—2013年61年实测洪峰流量，采用年最大值法选择实测洪峰流量系列，以此为基础进行洪峰流量频率分析，选择2010年洪水为典型洪水，按照洪峰控制的同倍比放大法得到高安站20年一遇设计洪水过程线见图 3。分析赣江樟树至外洲河段河道设计水面线，可得到市汊20年一遇设计水位为22.68 m。本文选择计算时段为2010-05-21T08:00/2010-05-28T08:00，计算方案为锦江20年一遇洪水，初始时刻计算区域淹没水深为0，计算步长30 s。

根据计算区域地形、地貌的实际情况，结合以往规划设计资料和经验值分析确定保护区糙率。计算区域地表采用不同下垫面不同的糙率分区，尽可能反映洪水演进的真实情况，网格糙率取值见表 1。

2.3 洪水演进计算结果

图 4为3个溃口处分洪流量变化过程。比较3条曲线可以看出，溃口1的洪峰流量最大，但溃口入流总量最小，溃口2入流总量最大，溃口3次之，这与区域地形分布相符。通过提取模型中水量平衡分析的内容，对药湖联圩防洪保护区内3个不同溃口的水量平衡进行了分析，具体见表 2。总体而言，模型对于水量平衡计算的相对误差基本控制在0.1%以内，误差相对较小，说明该模型模拟溃堤洪水的精度较高。

以溃口3为例，在锦江20年一遇洪水过程条件下，设定险工段司家闸发生瞬时溃堤，计算时长为7 d。从开始模拟时刻2010年5月21日8:00到23日3:40之间没有发生溃决，锦江水位正常升高，保护区水位(溃口附近的计算网格的水位)等于地面高程，溃口流量为零；在23日3:40，溃口3发生溃决，洪水从锦江涌入保护区，溃口流量由零瞬时达到最大，锦江水位降低，保护区内水位增高；溃决开始到24日18:00这一时段内，锦江水位高于保护区内水位，洪水从锦江流入保护区，但水位差逐渐减小，溃口流量逐渐减小；在24日18:00，保护区内水位反过来高于锦江水位，溃口流量出现负值，洪水从保护区回退入锦江，这是因为与保护区相比，锦江河道槽蓄量小，行洪能力大，退水时水位下降得更快；在25日17:20，保护区内水位降到溃口底高程齐平时，此时，洪水无法通过溃口从保护区退入锦江，溃口流量降为零。

2.4 洪水风险分析

洪水风险是指发生由洪水造成不同程度损失的可能性^[11-12]。根据3个溃口洪水分析计算结果及区域社会经济情况，并结合保护区地物分布，统计各方案的淹没面积、淹没耕地面积、受影响人口总数、最大淹没水位、最大淹没水深等洪水风险数据，计算结果见表 3。

综合分析药湖联圩防洪保护区锦江来洪洪水风险状况，可以发现区内洪水淹没范围大小和影响人口数量与溃口位置有关，从高到低依次为大王庙溃口洪水风险、司家闸溃口洪水风险、鸡鸣洲溃口洪水风险。具体原因分析为：

(1) 保护区内松湖镇地势平坦低洼，依靠铁臂圩与锦江松湖堤段保护，但铁臂圩防御标准仅10年一遇，在大王庙发生溃口时，洪水将会漫过铁臂圩，淹没范围最大，影响人口最多。

(2) 司家闸溃口位于保护区下游，由于该区域经常遭遇洪涝灾害，人口、工矿企业非常少，其发生溃口时洪水不会影响到松湖镇地区，淹没范围较小，影响人口相对较少。

(3) 鸡鸣洲发生溃堤时，南边山地和东边圳头导托排洪渠对洪水有非常大阻挡作用，淹没区域比较集中，淹没范围有限，影响人口最小。

另外，从防洪保护区内的淹没分布来看，基本遵从地形高低的原则，由于进洪周期比较长，进洪量一般都比较大，保护区内的淹没水深基本都大于3 m。

以溃口3为例，通过对计算结果进行统计分析，得到淹没过程中最大水深分布、淹没历时分布以及洪水前锋到达时间分布情况。将淹没水深图中满足水深≥1.0 m条件的区域作为危险区，选择具有容纳该区域内所有需避险人口的坚固楼房或公共设施(学校、体育馆、办公楼等)，可维持过水期间避险人员基本生活卫生条件以及常规的医疗需求的区域作为安置区。在分析药湖联圩防洪保护区地形特性基础上，对保护区进行了实地调研，发现保护区整体地形特征为保护区西侧、南侧为自然高地，且地形梯度很小。结合危险区的分析，可以将危险区周边的自然高地作为避险转移安置场所，并根据最短路径原则规划避险转移路径。最终将计算成果以洪水风险图的形式表达，部分区域成果如图 5所示。研究成果对防洪决策部门制定救灾抢险方案具有重大参考价值。

3 结语

本文对药湖联圩防洪保护区建立了溃堤洪水一、二维耦合水动力模型，利用糙率分区、河道地形还原等技术优化了模型，提高了模型精度。以锦江20年一遇洪水为基础，对计算区域进行了模拟，结果表明，此模型的计算能满足水量平衡，具有较高精度。通过模型计算得到结果，结合区域社会经济情况做了洪水风险分析，可以为防汛指挥部门安排抗洪救灾、避险转移及选择安置区等提供重要的信息。

与传统模型相比，该模型具有很好的多功能性和适用性，可应用于溃堤及其任意组合方案的洪水演进模拟并具有很好的稳定性，能满足洪水风险分析基本要求。但溃口过程是瞬时溃决，没有考虑溃口展宽，模型可以实现入流边界固定时溃口处内外水流的交换，但不能做到实时动态交换。因此，建议在地形数据能够满足计算精度要求的情况下，通过堤防决口处高程与宽度渐变关系反映溃口渐变溃决，模拟溃口处水流的复杂变化状态，同时模拟洪水在河道与洪泛区的交换过程及其在洪泛区内的演进情况。