耦合梯级水库应急调度与溃口洪水过程演进分析

Flood Routing Process Analysis Coupling Cascade Reservoir Emergency Operation and Breach Flow

摘要: 现有研究没有综合考虑流域梯级水库群溃坝应急调度与超标洪水河道演进模拟。本文将防洪调度和溃坝洪水演进计算结合起来，同时基于一维圣维南方程组，结合有限体积法、追赶法等方法，提出了一种适用于梯级水库群洪水演进的通用快速计算模型。该模型以调度结果作为洪水演进计算的边界条件，并根据实测资料对模型参数进行率定，实现了对洪水在河道内传播过程的精确模拟。以汉江流域断面实测数据与泉河流域的两座梯级水库为例，验证了所建立模型的合理性与准确性。结果表明，本文建立的快速计算模型能够实现对溃坝洪水河道演进全过程的精确模拟，并计算溃坝洪水影响范围与到达时间，有利于流域防洪减灾保障体系的构建，具有实际推广和应用价值。

Abstract: The existing researches have not synthetically considered dam-break emergency operation and exceeding flood routing simulation of cascade reservoirs. This paper coupled of the flood control operation and dam-break flow process calculation. Based on the one-dimensional Saint-Venant equations, Finite Volume and Chase-after methods, a general flood routing rapid simulation model for cascade reservoirs was proposed. The operation results were used as the boundary conditions for flood routing evolution calculation, and the parameters were calibrated based on the measured data to achieve an accurate simulation of the flood routing process. Taking the measured data of Hanjiang basin and two cascade reservoirs in Quan River basin as examples, the rationality and accuracy of the developed model were verified. The results indicated that the established rapid calculation model realizes the accurate simulation of the whole process of dam break flood channel routing and calculates the impact range and arrival time of dam-break floods, which is conducive to the construction of the basin flood control and disaster reduction security system, and has practical promotion and application value.

1. 引言

我国的水资源存在着时空分布极不均匀的现象 [1] ，极端气象事件呈现频率趋多、强度增大、范围变广等特征。由此诱发的洪涝灾害往往对流域内居民的生命财产安全与地区经济发展造成极大威胁。梯级水库群在流域防洪减灾中起着十分重要的作用 [2] ，但同时，由于大坝溃决往往会对下游人民群众安全和生态环境等方面造成巨大破坏，流域梯级水库群的溃坝应急调度与洪水演进分析需要得到重视 [3] [4] 。

现有的溃坝研究大多集中于溃坝风险分析 [5] 、避险方案研究 [6] 等，溃坝洪水演进模拟多数采用MIKE、HEC-GeoRAS、WCA-2D等平台或软件进行计算。夏军强等人 [7] 采用有限体积法求解所建立的二维水动力学模型，模拟了溃坝洪水的流动过程；吴天蛟等人 [8] 以MIKE11为基础，结合坝前水位与出库流量两种边界条件，针对三峡库区的洪水演进进行研究；贺娟等人 [9] 利用Google地球和HEC-RAS等计算工具，对长河坝电站水库大坝及下游区域进行溃坝洪水演进模拟，最后通过HEC-GeoRAS分析研究区域的洪水淹没范围及流速分布；王欣等人 [10] 基于MIKE FLOOD将MIKE11模型和MIKE21模型进行动态耦合，对溃坝洪水在深圳市龙华新区民治水库及下游片区的演进过程进行仿真模拟。陈祖煜等人 [11] 建立了能够快速计算土石坝溃决洪水的双曲线冲刷侵蚀模型。在溃口洪水流量过程的研究方面，国外学者Ritter [12] 、Stoker [13] 、Hunt [14] 与国内学者林秉南 [15] 、谢任之 [16] 等人先后在数学理论方面进行了研究；同时，美国、奥地利等国家 [17] ，长江水利委员会、黄河水利委员会水利科学院等 [18] 机构先后对溃坝问题进行了建模与试验，并总结了溃口洪水最大流量的计算公式。

当前的溃坝洪水演进研究大多利用各种已开发的软件，进行溃决水库下游库区的洪水演进模拟及灾害预测，缺乏对水库群防洪调度、溃口流量过程计算及洪水河道演进模拟的联合研究。基于此，本文耦合水库群防洪调度子模型与溃口洪水计算子模型，同时基于一维圣维南方程组，结合有限体积法、追赶法等方法，提出了一种适用于梯级水库群洪水演进的通用快速计算模型。以泉河流域红岩一级、红岩二级及大峡水库为研究对象，验证模型的合理性。结果表明，研究建立的耦合模型能实现流域梯级水库群溃坝应急调度、溃口洪水计算及溃坝洪水河道演进全过程的精确模拟。

2. 模型原理

2.1. 河道洪水演进模型

本文选择圣维南方程组作为基本控制方程建立了一维溃坝洪水演进模型，并选择直观且理论成熟有限体积法进行模型求解。基本控制方程包括质量守恒方程与动量守恒方程，具体公式如下：

$\{\begin{array}{l}B\frac{\partial h}{\partial t}+\frac{\partial Q}{\partial x}=q\\ \frac{\partial Q}{\partial t}+\frac{\partial }{\partial x}\left(\frac{Q^2}{A}\right)+gA\frac{\partial h}{\partial x}+gA\frac{n^{2}Q\left|Q\right|}{A{R}^{4/3}}=0\end{array}$ (1)

式中，B为明槽宽；h为水位；Q为流量；t为时间；x为流程；q为旁侧入流；A为过水断面面积；g为重力加速度；n为糙率；R为水力半径。

使用有限体积法离散后的连续性方程表达式如下：

$\{\begin{array}{l}-Q_{i-1}^{j+1}+G_i{h}_i^{j+1}+Q_i^{j+1}=G_i{h}_i^j+Q_{si}\\ -F_i{h}_i^{j+1}+D_i{Q}_i^{j+1}+F_i{h}_{i+1}^{j+1}=Q_i^{j^{\prime }}+Q_{si}\end{array}$ (2)

式中， $D_i=1+\frac{2^{7/3}{n}^i{n}^{i+1}g\Delta t{Q}_i^j}{\left(A_i^j+A_{i+1}^j\right){\left(R_i^j+R_{i+1}^j\right)}^{4/3}}$ ； $F_i=g\frac{\Delta t}{\Delta {x^{\prime }}_i}\cdot \frac{A_i^j+A_{i+1}^j}{2}$ ； $G_i=g\frac{\Delta x_i{B}_i^j}{\Delta t}$ ；i为河道断面编号；j为时刻数。

本研究采用有限体积法离散连续性方程，使用隐式方法离散动量方程的水位梯度项、使用欧拉–拉格朗日方法求解动量方程的对流项，并运用追赶法进行方程组求解。以河段上边界流量过程与河段下边界水位过程作为边界条件进行计算。同时，根据河段中间测站实测水位流量过程进行模型参数率定，洪水演进过程计算的流程图如图1所示。

Figure 1. Flowchart of the flood routing simulation model

图1. 洪水演进模型计算流程图

根据河道洪水演进、水库群防洪调度、溃口流量计算三个子模型，建立水库防洪调度–溃口洪水计算–溃坝洪水河道演进耦合模型，具体计算过程为：

1) 将梯级水库防洪调度和溃口洪水计算进行耦合，通过防洪调度得到不同调度方案下各个水库的坝前水位过程；

2) 将调度结果作为洪水演进模型的边界条件进行输入。运行河道洪水演进模拟，计算不同边界条件下河道内洪水的传播情况。

耦合模型流程图如图2所示。

Figure 2. Flowchart of the coupling model

图2. 耦合模型流程图

2.2. 防洪应急调度模型

以水库在各调度时段出库流量与入库流量之差的绝对值最小为水库群防洪应急调度模型的目标函数，具体公式如下所示：

$F_{i,j}=\min \left|Q_{ij}^{out}-Q_{ij}^{in}\right|$ (3)

式中， $F_{i,j}$ 为第i个水库在j时段的目标函数； $Q_{ij}^{in}$ 与 $Q_{ij}^{out}$ 分别为第i个水库在j时段的入库与出库流量。

考虑水库自身运行条件限制及其承担的下游防洪任务等，本文所建立的防洪调度模型需要满足水量平衡约束、上下游水力联系约束、出库流量约束、水位约束、出力约束及参数非负约束等。其中，水量平衡约束与上下游水力联系约束公式如下所示：

1) 水量平衡约束：

$V_{i,j+1}=V_{ij}+\left(Q_{ij}^{in}-Q_{ij}^{out}\right)\Delta t$ (4)

式中， $V_{i,j+1}$ 与 $V_{i,j}$ 分别为第i个水库在j、 $j+1$ 时段的库容； $\Delta t$ 为调度时段时长。

2) 上下游水力联系约束：

$Q_{ij}^{in}={\displaystyle \underset{k=1}{\overset{U_i}{\sum }}{Q}_{kj}^{out}}+B_{ij}$ (5)

式中， $Q_{ij}^{in}$ 为第i个水库在j时段的入库流量； $Q_{kj}^{out}$ 为第k个水库在j时段的出库流量； $B_{ij}$ 为第i个水库在j时段的上游区间来水流量； $U_i$ 为第i个水库的直接上游水库集合。

2.3. 溃口流量计算模型

选取考虑大坝溃决形式为在坝长和坝高方向同时发生瞬间局部溃坝的计算公式、美国水道试验站的修订公式及黄委水科院通过溃坝实验得到的计算公式 [19] ，以上述三式求平均值得到溃口最大流量 $Q^{b\max }$ ，并以四次抛物线概化溃口流量过程。具体计算过程如下：

$Q_1^{b\max }\text{=}\frac{8}{27}\sqrt{g}{\left(\frac{B}{b}\right)}^{1/4}\left(\frac{H_0-h^{\prime }}{H_0-0.827}\right)b\sqrt{H_0}\left(H_0-h^{\prime }\right)$ (6)

$Q_2^{b\max }=\frac{8}{27}\sqrt{g}{\left(\frac{B{H}_0}{bh}\right)}^{0.28}b{h}^{1.5}$ (7)

$Q_3^{b\max }\text{=}\frac{8}{27}\sqrt{g}{\left(\frac{B}{b}\right)}^{0.4}{\left(\frac{11H_0-10h}{H_0}\right)}^{0.3}b{h}^{1.5}$ (8)

$Q^{b\max }=\frac{Q_1^{b\max }+Q_2^{b\max }+Q_3^{b\max }}{3}$ (9)

式中， $b=k{\left(W^{1/2}{B}^{1/2}{H}^{1/2}\right)}^{1/2}$ ，为溃口宽度； $Q^{b\max }$ 为溃口最大流量；H₀为大坝溃决前上游水深； $h^{\prime }$ 为溃决后剩余大坝高度；H为溃坝前坝前水深；h为残坝坝顶至溃坝前水面的距离，也称为决口顶上水深，即 $h=H_0-h^{\prime }$ ；W为溃坝前溃口以上的库区蓄水量；B为坝顶长度；k为土质系数；g为重力加速度。具体溃口流量过程如图3所示。

Figure 3. Flow process line of dam site

图3. 坝址流量过程线

Figure 4. Distribution map of dam embankment in Quanhe basin

图4. 泉河流域坝堤分布图

3. 流域基本情况

泉河流域位于湖北省竹溪县境内中心腹地南部地区，河流全长82.2 km，流域面积894.6 km²。干流河势陡峭，落差集中，水利资源理论蕴藏量9.61万kw。汛期在每年的4月至10月，汛期径流量约占全年径流量的85.1%。流域洪水主要由汛期暴雨形成，洪水过程以单峰为主，时间与暴雨发生时间基本一致。流域内开发有红岩一级、红岩二级、大峡、龙滩和白沙河水电站等五座梯级水电站，其分布图如图4所示。其中，红岩一级水电枢纽为流域梯级开发的第一级电站，水库为中型水库，坝址以上流域面积340.4 km²，占全流域38.1%；红岩二级电站为第二级电站，为引水式电站；大峡水利枢纽为泉河流域梯级开发的第三级电站，控制流域面积482.7 km²，占全流域54%。

根据已有断面地形资料，考虑到红岩一级和大峡电站有一定的蓄水能力及闸门调节能力，红岩二级作为引水式电站，无防洪能力，研究以泉河流域红岩一级–大峡水库河段为研究对象，以此实现在流域突发超标洪水情境下，梯级水库群防洪调度–溃口计算–洪水演进全场景模拟。研究可为流域突发溃坝事故时的快速应急响应提供参考。

4. 结果分析

4.1. 河道演进子模型结果分析

汉江为长江第一大支流，流经区域广阔、水量丰沛且变化较大，流域内水文测站分布密集，水文资料丰富且各测站数据变化较为明显，具有可研性与代表性，且汉江流域洪水河道演进模型的构建与本文主要研究区域——泉河流域一致。因此选取汉江流域对本文提出的河道演进子模型进行准确性验证，将沙洋水文站至汉口水文站部分河段划分为249个断面，其中仙桃站位于96号断面，汉川站位于140号断面。以沙洋站流量过程作为上边界条件，汉口水位过程为下边界条件，对照仙桃站与汉川站进行结果验证。沙洋站~仙桃站，仙桃站~汉川站，汉川站~汉口站各河段各站点区间糙率分别为0.045，0.047，0.043。仙桃站与汉川站计算结果如图5所示。

Figure 5. Water level process of the Hankou and Xiantao station

图5. 仙桃站与汉川站水位过程对比

表1给出汉川站与仙桃站水位计算与实测结果，仙桃站的水位计算均值为22.05 m，与实测均值的差为0.05 m，峰值差为0.04 m；汉川站的水位计算均值为16.84 m，与实测均值的差为−0.04 m，峰值差为0.09 m。说明所建立的河道演进子模型，能够准确地模拟河段各断面的水位过程。

Table 1. Water level calculation results of Xiantao and Hanchuan control sections

表1. 仙桃和汉川站控制断面水位计算结果

Figure 6. Water level process of the section 1, 9 and 19

图6. 控制断面1、9、19的水位过程

Figure 7. Discharge process of the section 1, 9 and 19

图7. 控制断面1、9、19的流量过程

4.2. 耦合模型结果分析

以泉河流域红岩一级至大峡水库河段为研究对象，将该河段划分为19个控制断面，其中，断面1位于红岩一级水库下游，断面19位于大峡水库上游。以流域放大后的千年一遇72 h洪水过程为模型输入，首先，将防洪调度子模型和溃口洪水计算子模型进行耦合，计算得到各个水库的下泄流量过程和坝前水位过程；在此基础上，以红岩一级水库出库流量过程为上边界条件，以大峡水库的坝前水位过程为下边界条件进行洪水演进计算，得到各断面下的水位与流量过程。其中，由于河段缺乏实测糙率资料，根据经验假定河段糙率为0.045。断面1、9、19的水位与流量过程如图6、图7所示。

由耦合模型得到的断面1与19的水位过程可知，在模拟过程中红岩一级水库未发生溃坝，而大峡水库在第23小时左右发生溃决，坝前水位迅速降低，并在一段时间后保持较低水位。在洪水演进计算结果中，大峡控制断面的流量过程也在同一时刻迅速上升，峰现到达时间延后两个小时，且峰值流量大于红岩一级的出库流量峰值。结果表明，本文建立的耦合模型能实现对各控制断面水位、流量过程的快速、准确模拟。

5. 结论

本文建立了水库防洪调度–溃口洪水计算–溃坝洪水河道演进耦合模型，以泉河流域两座有调节能力的水库间河段为例，以流域超标洪水为模型输入，模拟了该河段19个控制断面在72 h内的洪水演进过程，主要得到了以下结论：

1) 以汉江流域249个断面为例，将洪水演进模型结果与仙桃站、汉川站实测数据进行对照。由结果可知，两个站点的水位过程趋势基本吻合；且峰值差分别为0.04 m与0.09 m，均值差分别为0.05 m与−0.04 m。结果表明，所建立的河道演进模型能够较为准确地模拟洪水传播过程。

2) 当泉河流域遭遇放大后的千年一遇洪水过程时，由耦合模型的模拟结果可知，红岩一级水库未发生溃坝事故，而大峡水库在23 h左右发生了溃决，各断面的水位与流量过程合理、准确。

3) 本文基于一维圣维南方程组，结合有限体积法、追赶法等方法，提出了一种适用于梯级水库群洪水演进的快速计算模型；并结合防洪调度、溃口流量计算与河道洪水演进三个子模型，考虑梯级水库调度特点及洪水演进模型边界条件，实现了对各控制断面水位流量过程的高效模拟，研究可估算溃坝洪水的影响范围和到达时间，为溃坝洪水应急处置提供依据。

基金项目

十四五国家重点研发计划(2021YFC3200400)；西藏自治区科技计划项目(XZ202301YD0044C)。