水力学与水文学结合的骆马湖预报调度系统

胡文才

（沂沭泗水利管理局，江苏 徐州 221009）

0 引言

近年来由于水利工程的大量修建，改变了河流与湖泊特性，在做湖泊洪水预报调度中发现，单纯采用水文学的方法做洪水预报调度，其精度完全不能满足实际防洪调度需求，因此不断有专家学者将水力学和水文学的方法结合起来研制洪水预报调度系统。20 世纪 80 年代初意大利的 E.Todini 教授与我国合作进行了淮河王家坝至正阳关河道洪水演进与滞洪区洪水调度，采用水力学与水文学集合的方法；2003 年河海大学李致家教授等在南四湖流域分别采用一维和二维水力学方法进行应用研究，取得了一定的成果[1-4]。

本次在骆马湖洪水预报调度中区间产流采用经验法计算至入湖口，南四湖和沂河来水采用马斯京根法计算至骆马湖入湖口。沂河港上、中运河运河站和房亭河土山站的入流分别作为网格入流处理。采用水量平衡的方法进行骆马湖洪水调度，以确定各个闸的放水流量过程，作为水力学的下边界条件，然后根据二维水力学方法演算骆马湖内的洪水过程以确定各个网格点的流量和水位。整个骆马湖的初始水位场由骆马湖内 9 个水位站的实时水位进行资料插值获得。

1 基本原理

骆马湖二维水动力学模型，采用平面二维水流数学模型进行模拟计算，为反映骆马湖内外的河道及岸线边界，采用边界拟合网格，建立正交曲线坐标系下的二维水流运动数学模型，运用有限差分法离散水流运动基本方程组，通过交替隐式解法（ADI）法求解离散方程[5]，得到整个计算区域上的水位和流速。

1.1 二维水动力学原理

1.1.1 水流运动基本方程

假定压力沿水深按静水压力分布，在正交曲线坐标系 ξ-η下的二维不恒定流连续性方程为：

式中：u，v分别为ξ和η方向的流速分量；Z 为水位；H 为水深，H = Z + h，h 为基面以下水深；q 为包括雨量等旁侧流量；g 为重力加速度；C 为谢才系数，C =H1/6，n 为糙率；εξ、εη为紊动粘性系数；ρ、ρα分别为水和空气的密度；vf为水面上 10 m 处风速；CD为风剪切应力系数，CD= (1.1 + 0.0536vf)×10-3。

1.1.2 ADI 法

采用有限差分法离散方程，通过 ADI 法求解方程。

图1 为正交曲线坐标下的流速和水深点的交错布置，对以水位、水深点为中心的阴影计算单元，流速 u 布置在 Z 方向的边上（图中横箭头），流速v布置在 h 方向的边上（图中竖箭头），圆点处为水位、水深和糙率等的位置，可以避免将变量布置在同一位置而引起的计算波动。

图1 交错网格布置

根据 ADI 法，在 n ～（n + 1/2）时间步长内对式 (1)、(2) 联立求解，变量 v以 n 时刻值代入，其差分方程为：

式中：i，j 为空间变量；A，B，C，D 是对原来连续性和动量方程经过这个交替隐式解法差分后的系数，无具体物理意义。

数学计算流程如图 2 所示。

1.2 马斯京根法

马斯京根法是一个流行的集总式河道演算方法，它假设了 1 个可变化的流量～蓄量方程：

取水量平衡方程和槽蓄方程差分解，可得流量演算方程：

图2 数学模型计算流程

式中：W 为槽蓄量；k 为蓄量常数；x 为求槽蓄量 W时，河段楔蓄量所占的权重；C0、C1、C2为河道演算系数。

1.3 降雨径流相关法

区间产流计算采用经验方法——降雨径流相关法，即 P + Pa～R。

2 具体应用

2.1 流域概况及水文特征

骆马湖位于沂河末端，中运河东侧，是以防洪、蓄水为主，结合航运、发电、水产养殖等综合利用的多功能湖泊。骆马湖承接沂河干流、南四湖、邳苍地区 5.1 万 km2面积的来水，调蓄后主要由新沂河排入黄海。入湖主要河道有沂河、中运河，出湖主要河道有新沂河、中运河和六塘河。

骆马湖南北长 20 km，东西宽 16 km，周长 70 km。一般湖底高程 20.0 m，最低 19.0 m。正常蓄水位 23.0 m 时，湖面面积 375 km，容积 9.0 亿 m3；退守宿迁大控制后，当达到设计洪水位 25.0 m 时，包括骆马湖与中运河之间三角地带面积 18 km2，湖面面积为 450 km2，容积 15 亿 m3；校核洪水位 26.0 m 时，总库容 19.0 亿 m3[6]。

2.2 入湖水量计算

骆马湖入湖水量主要由以下 3 部分组成：1）沂河来水，沂河临沂站来水除去彭道口闸和江风口闸分泄流量，剩余水量全部进入骆马湖；2）运河来水，中运河来水包括南四湖下泄水量和运骆邳苍区间产水；3）房亭河来水，房亭河刘集闸站以上流域面积 933 km2的区间产水。各控制断面位置如图 3 所示。

图3 骆马湖入湖水量控制断面示意图

图中 1、2、3 分别为 3 个入湖控制断面，其中断面 1 流量过程由临沂站过程采用马斯京根法演算得出，断面 2、3 流量过程分别由 P + Pa～R 计算和马斯京根法演算结果叠加得出，马斯京根法河段演算系数如表1 所示。

表1 马斯京根法河段演算系数表

2.3 模型应用范围

根据地形和工程的情况，以骆马湖湖区范围建立模型，建立如图 4 所示的骆马湖二维水动力模型网格的概化地形图，共有 30×38 个，网格尺寸190～1200 m。骆马湖二维水动力模型湖区地形采用 2001年1:10000 地形图，高程采用废黄河零点基面。模型的边界条件为新沂河嶂山闸、六塘河洋河滩闸、中运河皂河闸、中运河运河镇、沂河入湖口的流量过程。湖区库容的计算范围在骆马湖湖区内。

2.4 初始条件

初始条件如下：

图4 骆马湖二维水动力模型概化地形

式中：u0、v0分别为初始流速在 ξ 和 η 上的分量，计算时取流速 u0= 0 和 v0= 0；Z0为初始水位，为一给定值，由骆马湖现有的 9 个水位站（皂河闸、洋河滩闸、嶂山闸、苗圩、窑湾、袁场、张宅、新店、晓店）预报开始时刻的实测资料插值获得整个骆马湖的初始水位场。模型经过一定时间的运行，初始条件的影响将会消除。

2.5 动边界处理

在河湖地区，滩地和沙洲随着水位的涨落淹没和露出，在数模计算中要进行动边界处理，根据水位的变化连续不断地修正边界，模拟滩地水边线的移动。本文采用富裕水深法，即在计算中判断每个网格的水深，当网格水深大于某一给定的小水深时，将单元开放，作为计算水域；反之，将单元关闭，置流速于零。

2.6 边界条件

模型开边界一般用水位或流量过程控制。在骆马湖二维水动力数学模型中，出、入湖边界条件均用流量控制，出、入湖方向为 3 进 4 出，数学模型边界条件如图 5 所示。模型的边界条件为新沂河嶂山闸、六塘河洋河滩闸、中运河皂河闸、中运河运河镇、沂河入湖口、老沂河分洪道入湖口的流量过程。闭边界采用不可入条件，即取边界的外法线方向流速 vn= 0。

2.7 计算调度

骆马湖湖区洪水演进预报采用 96 h的时段长度，与水文学模型相结合，根据实测降雨过程作洪水预报，通过水文学模型计算骆马湖流域的产汇流状况，得出骆马湖湖区的入流边界条件，分别是沂河、中运河、骆马湖区间等入流过程及骆马湖湖面产流过程。然后再根据不同的洪水调度方案，得出骆马湖湖区的出流边界条件，最后模拟预测 96 h 内骆马湖湖区水位、流量的变化状况。

计算时首先根据洋河滩闸的 t 时刻的水位 Zt，利用骆马湖水位-库容曲线，求出 t 时刻的骆马湖的库容 Vt；再利用骆马湖预报模型中产汇流及湖区降雨量得出的 △t 时间内的入湖量 △V，求出 (t + △t) 时刻骆马湖的库容 Vt+△t，Vt+△t=Vt+△V；

最后通过骆马湖水位-库容曲线，预估出 (t + △t)时刻骆马湖洋河滩闸的水位 Zt+△t。

图5 骆马湖数学模型边界条件

根据上述骆马湖洪水调度原则和方案，结合预估出的 (t + △t) 时刻骆马湖洋河滩闸的水位 Zt+△t，确定具体的骆马湖湖区洪水的调度规则，具体调度规则如表2 所示，再通过骆马湖二维水动力数学模型计算，预测当前骆马湖各种调度规则下，湖区的洪水水位及流量的变化过程。

表2 骆马湖湖区水动力洪水调度规则

2.8 实例应用

选择 2006、2008 和 2009 年，共 5 场洪水进行模拟预报调度，预报调度结果如表3 所示。由于篇幅问题示例图仅列 2009年2 场洪水预报调度过程图，具体如图 6、7 所示。

表3 实际应用成果统计表

图6 2009-07-17 洪水预报调度结果

图7 2009-08-17 洪水预报调度结果

模拟预报结果显示，水位预报误差最大为 0.13 m。预报最高水位出现时间，与实际最高水位出现时间相比，有的提前有的滞后，最大相差 1 d。出现这一结果主要有 3 个方面的原因：1）入湖过程滞后，模型计算港上站至骆马湖洪水传播时间为 6 h，但是实际传播时间为 12 h 左右；2）实际与计算调度过程不可能完全一致；3）预报入湖过程存在误差，入湖过程预报采用经验方法，其预报结果与实际入湖过程存在一定的误差。以 2009年7月17日洪水为例，预报入湖时间为 22日 5:00，但是实际入湖时间为 22日 11:00 左右，滞后 6 h；预报入湖洪峰流量为 5835 m3/s，但实际最大入湖洪峰流量为 4600 m3/s左右，预报入湖洪峰流量比实际入湖洪峰大1235 m3/s。由于实际入湖过程比预报过程滞后，导致计算结果最高水位出现时间比实际时间提前，预报最大入湖洪峰流量大于实际入湖洪峰流量，导致预报最高水位高于实际最高水位。

2.9 优缺点

该系统采用水力学与水文学结合的方法进行预报调度计算，湖内洪水演进采用水力学的方法，考虑了动库容的问题，解决了传统水文学方法计算时由于动库容带来的误差问题，使预报精度提高了。该系统水力学计算采用 Fortran 语言编程，计算速度提高了，目前做 1 次预报时间大约仅需要 2 min，与传统预报系统相比时间缩短了。

其他相似湖泊要应用该系统，需要考虑入湖流量和时间误差，以及实际出湖过程与预报过程不一致的情况。如果考虑上述几个因素带来的影响，进行预报结果修正，将会带 来更好的预报成果。

3 结语

本文将二维水流运动数学模型，有限差分法离散水流运动基本方程组，ADI 法求解离散方程的方法和水文学的方法结合应用于骆马湖洪水预报调度中；地形资料采用了 2001年1:10000 地形图，与实际地形基本吻合；采用 9 个水位控制站的实时水位，通过插值获取计算初始时刻骆马湖的水位场。在实际应用中发现二维水力学进行湖泊水力学演算，如果地形资料准确，湖泊内水位观测站点代表性较好，采用该方法能够较好地模拟各个测点的水位。在实时洪水预报和调度中，采用水文学和水力学结合的方法进行洪水预报调度，能取得令人满意的预报成果。

[1]李致家，董增川，梁忠民，等. 大流域洪水预报与洪水调度管理方法研究[J]. 河海大学学报，2004, 30 (1): 12-15.

[2]李致家. 通用一维河网不恒定流软件的研究[J]. 水利学报，1998 (8): 14-18.

[3]汪德灌. 计算水力学-理论与应用[M]. 南京：河海大学出版社，1989: 64-101.

[4]李致家，包红军，孔祥光，等. 水文学与水力学相结合的南四湖洪水预报模型[J]. 湖泊科学，2005, 17: 299-304.

[5]徐峰俊，刘俊勇. 伶仃洋海区二维不平衡非均匀输沙数学模型[J]. 水力学报，2003 (7): 16-23.

[6]郑大鹏. 沂沭泗防汛手册[M]. 徐州：中国矿业大学出版社，2003: 95-97.