CN103106625A - 水库、闸泵群联合抑咸调度方法 - Google Patents

Abstract

本发明公开了一种水库、闸泵群联合抑咸调度方法，包括下述步骤：S1、建立基于EasyDHM的珠江流域分布式水文模型；S2、建立骨干水库群抑咸调度模型；S3、利用基于EasyDHM的珠江流域分布式水文模型计算出流域关键节点的流量过程，输入到水库群抑咸调度模型中，通过模型得到下游关键断面的流量控制过程；S4、利用闸泵群释淡补水压咸调度模型确定放水水闸、释淡补水压咸时机和放水压咸流量。本发明分析了珠江流域枯季径流特征，并对流域关键断面枯水期来水过程进行了预报；构建了流域骨干水库群优化调度模型，在优化调度方案的基础上实现了实时调度方案的耦合与嵌套和滚动修正，克服了确定性调度模型带来的不足，大大提高了调度期内水资源利用率。

Description

水库、闸泵群联合抑咸调度方法

技术领域

本发明涉及环保领域的水质保障技术，特别涉及河口咸潮抑制以及保障供水安全的水库、闸泵群联合抑咸调度方法。

背景技术

珠江下游地区是我国经济和社会高速发展的地区，在我国经济建设中具有举足轻重的战略地位。随着经济社会的发展，工农业生产和生活用水量将进一步增加，珠江河口地区的水体污染问题也日益严峻，饮水安全问题更加凸显，水资源已成为制约该地区社会经济发展的瓶颈。枯水期咸潮上溯严重加剧了水资源短缺的局面，凸现了珠江三角洲的饮水安全问题。尤其是，近来年咸潮活动越来越频繁，影响范围愈来愈大，且持续时间愈来愈长，咸潮上溯已成为港澳及珠三角地区饮水安全中亟待解决的突出问题。2006年至2012年枯水期咸潮再度肆虐，国家防汛抗旱总指挥部授权珠江防汛抗旱总指挥部先后实施珠江流域骨干水库调度和珠江枯季水量统一调度，以抵御咸潮。

实践证明，科学调配流域水资源是现阶段抑咸的有效途径。通过水库调度抑制咸潮上溯，并改善珠江三角洲地区水生态和水环境质量的非工程措施的实施，在珠江流域刚刚起步，在国内外其它地区也尚未见到比较系统的报道和科学研究工作。珠江流域的多次调水虽然大大缓解和基本抑制了近年来的咸潮危害，并对枯水季节的水质污染和生态危害起到了一定的抑制作用，但打破水库原有功能的调度也造成一定的损失。目前水资源调度模式与多数水库以防洪、灌溉、发电为主要功能的调度运行方式不相适应，不能满足咸潮上溯期饮用水的需要；此外，淡水资源分布与取水口布局不相适应，供水系统布局不合理，配套水源的建设滞后，导致淡水资源未能得到充分利用，解决枯水期珠江三角洲地区饮水安全问题已成为该地区经济发展和社会稳定的迫切要求。

开展以水库—闸泵群联合调度抑咸为核心的关键技术研究，系统研究抑咸的技术方法体系，科学调度，提高珠江河口河网区淡水资源利用效率，是一项十分重要和迫切的研究工作，对于保障珠江下游地区饮用水安全具有重要意义。

发明内容

本发明的目的在于克服现有技术的缺点与不足，提供一种以水库—闸泵群联合抑咸调度为核心的关键技术，通过流域和区域的抑咸水资源调配技术方法，建立咸潮抑制技术体系，通过科学调度，提高珠江河口河网区淡水资源利用效率。

本发明的目的通过下述技术方案实现：

本发明的水库、闸泵群联合抑咸调度方法，包含三个关键的模型，即基于EasyDHM的珠江流域分布式水文模型、骨干水库群抑咸调度模型和闸泵群释淡补水压咸调度模型，水库、闸泵群联合抑咸调度方法，包括下述步骤：

S1、建立基于EasyDHM的珠江流域分布式水文模型；

S2、建立骨干水库群抑咸调度模型，骨干水库群抑咸调度模型包括两个方面：一个是模型的表达形式，一个是模型的求解方法，模型的表达形式包括目标函数和约束条件；

S3、通过流域水文气象部门发布的实时水雨情信息，利用基于EasyDHM的珠江流域分布式水文模型计算出流域关键节点的流量过程，输入到水库群抑咸调度模型中，通过模型得到下游关键断面的流量控制过程；

S4、利用闸泵群释淡补水压咸调度模型确定放水水闸、释淡补水压咸时机和放水压咸流量。

优选的，步骤S1中，所述EasyDHM的珠江流域分布式水文模型实现步骤如下：

S21、基础数据准备：

选取DEM数字高程数据；设定最小水道集水面积阈值，生成河网水系栅格图层；首先把该土地利用类型数据转化为网格尺寸为1km×1km栅格数据格式，然后将土地类型转化为EasyDHM产流模型中相应的土地利用分类，最后根据EasyDHM产流模型的土地利用参数推求方法求各项产流参数；

S22、空间单元划分：

采用参数分区拓扑式率定方法，在基于DEM提取的数字流域基础上确定出每个水文站和水库的控制范围，进而划分出各个参数分区和水库分区；

S23、参数率定：

采用Nash效率系数来衡量模型模拟精度的好坏，其方程表达如下：

其中：R2是用于衡量模型模拟实际时间序列上径流量能力的纳什效率系数，纳什效率系数R2的变化范围是负无穷到1，若为1，则表示模拟的和实测的水文过程线达到完美拟合，R2小于0则表示平均实测径流量会是模拟径流量的一个较好预测。

优选的，所述的骨干水库群抑咸调度模型的调度目标包括抑咸调度效果、梯级电站发电经济效益、区域经济发展和生态环境，其目标函数描述为：

OBJ_Fun=Q_抑咸+E_发电+Q_航运+Q_生态

其中：OBJ_Fun表示目标函数，Q_抑咸表示最佳抑咸流量，E_发电表示发电量最大，Q_航运表示航道航运流量要求，Q_生态表示河道生态流量要求。

优选的，步骤S2中，所述目标函数包括发电效益目标、抑咸流量目标和生态环境目标，其具体如下：

发电效益目标：其表达式为

其中，A(m)为出力系数，QD(m,t)为第m个水库t时段的发电流量，H(m,t)为m个水库t时段的平均发电水头，△t为计算步长，M为电站个数，T为总时段数；

抑咸流量目标：在西江梧州站控制断面流量2100m³/s，石角控制断面流量250m³/s，思贤滘控制断面流量2500m³/s；

生态环境目标：西江梧州、北江石角以及西北江三角洲思贤滘三个控制断面的非汛期生态环境流量分别为1800m³/s、250m³/s、2200m³/s。

优选的，步骤S2中，所述模型约束条件具体如下：

（1）目标转换约束；

（2）水库水量平衡约束：

V(m,t+1)=V(m,t)+RW(m,t)-W(m,t)-LW(m,t)

W(m,t)=q(m,t)*△t

V(m,t)、RW(m,t)、W(m.,t)、LW(m,t)分别为第m个水库t时段库容、入库水量、出库水量和损失水量；q(m,t)为第m个水库t时段的出库流量；

（3）出库流量约束：

QDmin(m,t)≤QD(m,t)≤QDmax(m,t)

qmin(m,t)≤q(m,t)≤qmax(m,t)

QDmin(m,t)、QDmax(m,t)分别为第m个水库t时段最小、最大允许过机流量；qmin(m,t)、qmax(m,t)分别为第m个水库t时段最小、最大允许出库流量；

（4）出力约束：

Nmin(m,t)≤N(m,t)≤Nmax(m,t)

$\underset{i=1}{\overset{M}{\mathrm{\Σ}}}N(m,t)\≥\mathrm{NSUM}\min \left(t\right)$

N(m,t)、Nmin(m,t)、Nmax(m,t)分别为第m个水库t时段出力、允许最小和最大出力；NSUMmin(t)为梯级t时段允许最低总出力；

（5）水库库容约束：

Vmin(m,t)≤V(m,t)≤Vmax(m,t)

Vmin(m,t)、Vmax(m,t)分别为第m个水库t时段允许库容上下限；

（6）河道水量演进约束：

Q(i+1,t+1)=C₀Q(i,t+1)+C₁Q(i,t)+C₂Q(i+1,t)

∑C=1

Q(i,t)为第i个节点t时段的流量；

（7）变量非负约束。

优选的，步骤S2中，所述模型的求解方法具体为：

S51、根据区间来水及控制断面目标流量要求，考虑水量传播和水量损失因素，自下而上推求各库供水约束下限，并结合发电、生态约束拟定初始调度线；

S52、根据给定的天生桥一级计算期末消落水位，采取自上而下、逆时序的方法，推求各库时段最低、最高水位控制线；

S53、按初始调度线自上而下开始顺时序模拟梯级运行过程，计算时段末各库状态和梯级出力，若模拟结果经水位和出力辨识满足要求，则进入下一时段，否则，按一定规则加入反馈修正量，重新模拟时段运行过程，直到满足水位和出力辨识要求；

S54、重复步骤S53逐时段迭代模拟—反馈修正，直到计算期末，完成一轮迭代；

S55、进行目标辨识，若模拟期末天生桥一级水位与给定龙库期末水位之差满足要求，则结束；否则，形成修正量并反馈到输入端，从计算期初重新新一轮迭代，直至期末水位满足要求。

优选的，步骤S4中，所述放水水闸的确定是根据补水对象而定，具体是针对取水口释淡补水压下，分析调度期咸界位置、各取水口缺水紧迫程度以及压咸效果，综合确定；所述释淡补水压咸时机是根据补水对象不同时段淡水需求、咸潮上溯强度综合确定；所述放水压咸流量必须能够满足补水目标断面的最小抑咸流量需求，补水流量大小则通过闸门开度和外排泵站抽排流量来控制。

优选的，所述的闸泵群释淡补水压咸调度模型包括河口区二维潮流模块、网河区一维水动力数学模块以及闸泵群联合调控模块三个组成部分；

其中：

（1）河口区二维潮流模块采用贴体正交曲线坐标系下的二维潮流控制方程，形式如下：

连续方程

动量方程

$\frac{\∂\left(\mathrm{Hu}\right)}{\∂t}+\frac{1}{C_{\mathrm{\ζ}}{C}_{\mathrm{\η}}}\[\frac{\∂}{\∂\mathrm{\ζ}}\left(C_{\mathrm{\η}}\mathrm{Huu}\right)+\frac{\∂}{\∂\mathrm{\η}}\left(C_{\mathrm{\ζ}}\mathrm{Hvu}\right)+\mathrm{Hvu}\frac{\∂C_{\mathrm{\ζ}}}{\∂\mathrm{\η}}-H{v}^2\frac{\∂C_{\mathrm{\η}}}{\∂\mathrm{\ζ}}\]+\frac{\mathrm{gu}\sqrt{u^2+v^2}}{C^2}$

$+\frac{\mathrm{gH}}{C_{\mathrm{\ζ}}}\frac{\∂h}{\∂\mathrm{\ζ}}-\mathrm{fvH}=\frac{1}{C_{\mathrm{\ζ}}{C}_{\mathrm{\η}}}\[\frac{\∂}{\∂\mathrm{\ζ}}\left(C_{\mathrm{\η}}H{\mathrm{\σ}}_{\mathrm{\ζ\ζ}}\right)+\frac{\∂}{\∂\mathrm{\η}}\left(C_{\mathrm{\ζ}}H{\mathrm{\σ}}_{\mathrm{\ζ\η}}\right)+H{\mathrm{\σ}}_{\mathrm{\ζ\η}}\frac{\∂C_{\mathrm{\ζ}}}{\∂\mathrm{\η}}-H{\mathrm{\σ}}_{\mathrm{\η\η}}\frac{\∂C_{\mathrm{\η}}}{\∂\mathrm{\ζ}}\]$

$\frac{\∂\left(\mathrm{Hv}\right)}{\∂t}+\frac{1}{C_{\mathrm{\ζ}}{C}_{\mathrm{\η}}}\[\frac{\∂}{\∂\mathrm{\ζ}}\left(C_{\mathrm{\η}}\mathrm{Huv}\right)+\frac{\∂}{\∂\mathrm{\η}}\left(C_{\mathrm{\ζ}}\mathrm{Hvv}\right)+\mathrm{Huv}\frac{\∂C_{\mathrm{\η}}}{\∂\mathrm{\ζ}}-H{u}^2\frac{\∂C_{\mathrm{\ζ}}}{\∂\mathrm{\η}}]+\frac{\mathrm{gv}\sqrt{u^2+v^2}}{C^2}$

$+\frac{\mathrm{gH}}{C_{\mathrm{\η}}}\frac{\∂h}{\∂\mathrm{\η}}+\mathrm{fuH}=\frac{1}{C_{\mathrm{\ζ}}{C}_{\mathrm{\η}}}[\frac{\∂}{\∂\mathrm{\ζ}}\left(C_{\mathrm{\η}}H{\mathrm{\σ}}_{\mathrm{\ζ\η}}\right)+\frac{\∂}{\∂\mathrm{\η}}\left(C_{\mathrm{\ζ}}H{\mathrm{\σ}}_{\mathrm{\η\η}}\right)+H{\mathrm{\σ}}_{\mathrm{\ζ\η}}\frac{\∂C_{\mathrm{\η}}}{\∂\mathrm{\ζ}}-H{\mathrm{\σ}}_{\mathrm{\ζ\ζ}}\frac{\∂C_{\mathrm{\ζ}}}{\∂\mathrm{\η}}]$

式中：u、v分别为ζ、η方向的流速分量；h为水位；H为水深；g为重力加速度；f为柯氏力系数；系数C_ζ、C_η如下：

$C_{\mathrm{\ζ}}=\sqrt{x_{\mathrm{\ζ}}^2+y_{\mathrm{\ζ}}^2}$

$C_{\mathrm{\η}}=\sqrt{x_{\mathrm{\η}}^2{y}_{\mathrm{\η}}^2}$

σ_ζζ、σ_ηη、σ_ζη、σ_ηζ为应力项，其表达式如下：

${\mathrm{\σ}}_{\mathrm{\ζ\ζ}}=2v_t[\frac{1}{C_{\mathrm{\ζ}}}\frac{\∂u}{\∂\mathrm{\ζ}}+\frac{v}{C_{\mathrm{\ζ}}{C}_{\mathrm{\η}}}\frac{\∂C_{\mathrm{\ζ}}}{\∂\mathrm{\η}}]$

${\mathrm{\σ}}_{\mathrm{\η\η}}==2v_t[\frac{1}{C_{\mathrm{\η}}}\frac{\∂v}{\∂\mathrm{\η}}+\frac{u}{C_{\mathrm{\ζ}}{C}_{\mathrm{\η}}}\frac{\∂C_{\mathrm{\η}}}{\∂\mathrm{\ζ}}]$

${\mathrm{\σ}}_{\mathrm{\ζ\η}}={\mathrm{\σ}}_{\mathrm{\η\ζ}}=v_t[\frac{C_{\mathrm{\η}}}{C_{\mathrm{\ζ}}}\frac{\∂}{\∂\mathrm{\ζ}}\left(\frac{v}{C_{\mathrm{\η}}}\right)+\frac{C_{\mathrm{\ζ}}}{C_{\mathrm{\η}}}\frac{\∂}{\∂\mathrm{\η}}\left(\frac{u}{C_{\mathrm{\ζ}}}\right)]$

其中，v_t为紊动粘性系数，即：v_t=au_*H，式中：a为系数，u_*为摩阻流速，H为水深；

（2）网河区一维水动力数学模块采用一维圣维南方程组，方程如下：

连续方程 $B\frac{\∂Z}{\∂t}+\frac{\∂Q}{\∂x}=q$

动量方程 $\frac{\∂Q}{\∂t}+\frac{\∂}{\∂x}\left(\mathrm{\β}\frac{Q^2}{A}\right)+\mathrm{gA}(\frac{\∂Z}{\∂x}+S_f)+u_{l}q=0$

式中：Z为断面平均水位；Q、A、B分别为断面流量、过水面积、水面宽度；x、t为距离和时间；q为旁侧入流，负值表示流出；β为动量校正系数；g为重力加速度；S_f为摩阻坡降，采用曼宁公式计算，S_f=g/C²，C=h^1/6/n；u_l为单位流程上的侧向出流流速在主流方向的分量；

（3）闸泵群调控模块

闸泵群联合调度的总体目标是以外江典型水文过程条件，根据上游来水和河口咸潮运动情况，合理控制多汊河口联围内外闸泵群的启闭，适时将外江淡水资源引入联围，置换内河涌污水，尽可能蓄积最多的淡水，并于压咸期间释放淡水抑制咸潮上溯，保障取水口取水安全。

优选的，所述的闸泵群释淡补水压咸调度模型的目标函数为：

$\min C_p=(W_0\·C_0+\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\∫}_{t_0}^{t_1}(Q_{i,t}\·C_{i,t})\mathrm{dt})/(W_0+\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\∫}_{t_0}^{t_1}{Q}_{i,t}\mathrm{dt})$

式中：

C_p是联围内河涌污染物平均浓度，即换水结束时刻t₁时的联围内河涌污染物平均浓度；

W₀是换水起始时刻t₀时的联围内河涌蓄水量；

C₀是换水起始时刻t₀时的联围内河涌污染物平均浓度；

Q_i,t是t时刻第i座进/排水闸或泵的进/排水流量，以正负号表示流向，正值表示水从外江流入至联围内河涌，负值表示水从联围内河涌流出至外江，i=1，2，…，N，N为闸或泵总数；

C_i,t是t时刻第i座进/排水闸或泵的进/排水污染物浓度。

优选的，所述的闸泵群释淡补水压咸调度模型的约束条件为：

①换水历时约束

T_r≤T_rmax

②联围内河涌最高最低水位约束

${\underset{\‾}{Z}}_{k,t}\≤Z_{k,t}\≤{\stackrel{\‾}{Z}}_{k,t}$

③进水闸进水条件约束

$S_{i,t}\≤\stackrel{\‾}{S}$

Z_i,外,t>Z_i,内,t

④排水闸排水条件约束

Z_j,外,t<Z_j,内,t

⑤闸或泵运行状态持续时间约束

$T_s\&GreaterEqual;{\underset{\&OverBar;}{T}}_s$

⑥闸门启闭速度约束

${\underset{\&OverBar;}{T}}_d\&le;T_d\&le;{\stackrel{\&OverBar;}{T}}_d$

⑦闸或泵安全运行条件约束

$\left|{\mathrm{\&Delta;Z}}_{i,t}\right|\&le;\mathrm{\&Delta;}{\stackrel{\&OverBar;}{Z}}_i$

$\left|Q_{i,t}\right|\&le;{\stackrel{\&OverBar;}{Q}}_i$

式中：T_r为换水历时，T_r=t₁-t₀；T_rmax为允许的最长换水历时；Z_k,t、

分别为t时刻第k条内河涌（河段）水位、最高限制水位、最低限制水位，S_i,t、

分别为t时刻第i个进水闸或泵进水口的含氯度、取水允许的最高含氯度；Z_i,内,t、Z_i,外,t分别为t时刻第i个进水闸内、外水位；Z_j,内,t、Z_j,外,t分别为t时刻第j个排水闸内、外水位；T_s、

分别为闸或泵维持某一特定工作状态的时长和允许最短时长；T_d、

分别为闸门启闭单位高度的用时、允许的最小和最大用时；△Z_i,t、

分别为t时刻第i个闸内外水位差、闸门安全运行允许的最大水位差，Q_i,t、

分别为t时刻第i个闸或泵过水流量、设计最大过水流量。

本发明相对于现有技术具有如下的优点及效果：

1、本发明开创性地提出了水库、闸泵群联合调度抑咸技术体系理念，从流域和区域两个层面实现水资源调配，经多年的调度实践，证明有效。集成了流域分布式水文模型、枯季河道水流演进模型、流域骨干水库群优化调度模型、河口一二维潮流模型和联围闸泵群联合调度模型，构建了水库、闸泵群抑咸联合调度模型系统。

2、本发明建立了基于EasyDHM的珠江流域分布式水文模型，分析了珠江流域枯季径流特征，并对流域关键断面枯水期来水过程进行了预报；构建了流域骨干水库群优化调度模型，在优化调度方案的基础上实现了实时调度方案的耦合与嵌套和滚动修正，克服了确定性调度模型带来的不足，大大提高了调度期内水资源利用率。

3、本发明构建了量质同控的闸泵群联合抑咸调度模型，优化确定了联围内河涌引水冲污、开闸蓄淡、释淡抑咸等抑咸关键调度时机，实现了抑咸时长最大、工作量最省、可操作性最强、抑咸效果最好的最优闸泵群联合调控。

附图说明

图1是本发明实施例所在位置的示意图（中顺大围）；

图2是本发明实施例的骨干水库调度的节点概化图；

图3是本发明实施例的骨干水库调度模型自优化求解方法流程图；

图4是本发明实施例调度期间梧州断面预报流量过程；

图5是本发明实施例调度期间思贤滘断面预报流量过程；

图6是本实施例调度期间东河闸内外水位过程；

图7是本实施例调度期间西河闸内外水位过程；

图8是本实施例调度期间铺锦闸内外水位过程。

具体实施方式

下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。

实施例

本发明的水库、闸泵群联合抑咸调度方法包含三个关键的模型，即基于EasyDHM的珠江流域分布式水文模型、骨干水库群抑咸调度模型和闸泵群释淡补水压咸调度模型，是基于EasyDHM的珠江流域分布式水文模型、骨干水库群抑咸调度模型和闸泵群释淡补水压咸调度模型三个模型自上而下有机组成的整体，其具体包括如下步骤：

（1）建立基于EasyDHM的珠江流域分布式水文模型，EasyDHM产流模型是在Wetspa、SWAT、新安江等产流模型的基础上，进行集成创新而提出的一种产汇流模型，该模型在不同地区、不同水文地质条件下均能通用。模型实现步骤如下：

①基础数据准备：DEM数字高程数据来自于美国联邦地质调查局(USGS)的HYDRO1k（网址：http：//edcdaac.usgs.gov/gtopo30/hydro/）；设定最小水道集水面积阈值为150km²，生成河网水系栅格图层；土地利用源信息采用来自《中国资源环境遥感宏观调查与动态研究》课题的研究成果数据—全国分县土壤覆盖矢量数据。首先把该土地利用类型数据转化为网格尺寸为1km×1km栅格数据格式，然后将土地类型转化为EasyDHM产流模型中相应的土地利用分类，最后根据EasyDHM产流模型的土地利用参数推求方法求各项产流参数；土壤栅格数据主要采用的是全球FAO土壤数据；水文气象资料采用珠江流域范围内的170个气象站、14个水文站，气象资料范围从1954-2008年，共55年的降水、径流、水面蒸发日观测实测资料。

②空间单元划分：采用最小水道集水面积阈值为60km²时，珠江流域共划分为4649个子流域；在EasyDHM中为了体现模型参数在空间上的变异性，采用了参数分区拓扑式率定方法。在前面基于DEM提取的数字流域基础上确定出每个水文站和水库的控制范围，进而划分出各个参数分区和水库分区。

③参数率定：为了衡量模型模拟精度的好坏，引入Nash-Sutcliffe效率系数（即Nash效率系数），其方程表示如下：

其中：R2是用于衡量模型模拟实际时间序列上径流量能力的纳什效率系数。纳什效率系数R2的变化范围是负无穷到1。若为1，则表示模拟的和实测的水文过程线达到完美拟合。R2小于0则表示平均实测径流量会是模拟径流量的一个较好预测。

（2）建立骨干水库群抑咸调度模型，骨干水库群抑咸调度模型包括两个方面：一个是模型的表达形式，一个是模型的求解方法。模型的表达形式包括目标函数和约束条件，其中目标函数包括抑咸调度效果、梯级电站发电经济效益、区域经济发展和生态环境等。

平均发电水头，△t为计算步长，M为电站个数，T为总时段数。

抑咸流量目标：在西江梧州站控制断面流量2100m³/s，石角控制断面流量250m³/s，思贤滘控制断面流量2500m³/s。

生态环境目标：西江梧州、北江石角以及西北江三角洲思贤滘三个控制断面的非汛期生态环境流量分别为1800m³/s、250m³/s、2200m³/s。

模型约束条件：

①目标转换约束

②水库水量平衡约束

V(m,t+1)=V(m,t)+RW(m,t)-W(m,t)-LW(m,t)

W(m,t)=q(m,t)*△t

V(m,t)、RW(m,t)、W(m.,t)、LW(m,t)分别为第m个水库t时段库容、入库水量、出库水量和损失水量；q(m,t)为第m个水库t时段的出库流量。

③出库流量约束

QDmin(m,t)≤QD(m,t)≤QDmax(m,t)

qmin(m,t)≤q(m,t)≤qmax(m,t)

QDmin(m,t)、QDmax(m,t)分别为第m个水库t时段最小、最大允许过机流量；qmin(m,t)、qmax(m,t)分别为第m个水库t时段最小、最大允许出库流量。

④出力约束

Nmin(m,t)≤N(m,t)≤Nmax(m,t)

$\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}N(m,t)\&GreaterEqual;\mathrm{NSUM}\min \left(t\right)$

N(m,t)、Nmin(m,t)、Nmax(m,t)分别为第m个水库t时段出力、允许最小和最大出力；NSUMmin(t)为梯级t时段允许最低总出力。

⑤水库库容（水位约束）

Vmin(m,t)≤V(m,t)≤Vmax(m,t)

Vmin(m,t)、Vmax(m,t)分别为第m个水库t时段允许库容上下限。

⑥河道水量演进约束

Q(i+1,t+1)=C₀Q(i,t+1)+C₁Q(i,t)+C₂Q(i+1,t)

∑C=1

Q(i,t)为第i个节点t时段的流量。

⑦变量非负约束。

模型求解方法：模型求解采用自优化理论，其基本思路是：首先根据区间来水及控制断面目标流量要求，考虑水量传播和水量损失等因素，自下而上（逆流向）推求各库供水约束下限，并结合发电、生态等约束拟定初始调度线；然后再根据给定的天生桥一级计算期末消落水位，采取自上而下（顺流向）、逆时序（由计算期末到期初）的方法，推求各库时段最低、最高水位控制线；按初始调度线自上而下开始顺时序模拟梯级运行过程，计算时段末各库状态和梯级出力，若模拟结果经水位和出力辨识满足要求，则进入下一时段，否则，按一定规则加入反馈修正量，重新模拟时段运行过程，直到满足水位和出力辨识要求。如此逐时段迭代模拟—反馈修正，直到计算期末，完成一轮迭代。最后进行目标辨识，若模拟期末天生桥一级水位与给定龙库期末水位之差满足要求，则结束；否则，形成修正量并反馈到输入端，从计算期初重新新一轮迭代，直至期末水位满足要求。

本模型采用了三层辨识反馈结构。首先是各时段末的水库水位辨识，将水库时段末水位控制在最高与最低水位控制线之间，以保证不使供水破坏。若不满足辨识要求，该结构将返回一个修正量重新模拟系统运行。第二层是时段出力辨识模拟出力若达不到系统最小出力则返回一个修正量，如达到预定出力则继续下一时段模拟。第三层是目标辨识优化，根据天生桥一级水库预期期末水位和实际模拟水位，对调度期平均出力进行寻优。这三层辨识反馈只需给定允许误差，模型将自动迭代寻优，直到满足目标要求。

（3）通过流域水文气象部门发布的实时水雨情信息，利用基于EasyDHM的珠江流域分布式水文模型计算出流域关键节点的流量过程，输入到水库群抑咸调度模型中，通过模型得到下游关键断面的流量控制过程；

（4）闸泵群释淡补水压咸调度模型，联围释淡补水压咸方案主要内容是确定放水水闸、释淡补水压咸时机和放水压咸流量；放水水闸的确定需根据补水对象而定，具体方案实施中针对哪一取水口（取水河段）释淡补水压下，可分析调度期咸界位置、各取水口缺水紧迫程度以及压咸效果，综合确定。同样，释淡补水压咸时机需要根据补水对象不同时段淡水需求、咸潮上溯强度综合确定；放水压咸流量必须能够满足补水目标断面的最小抑咸流量需求，该流量的确定是闸泵群联合调度抑咸方案制定的关键点之一。在满足该最小流量需求的基础上，补水时间越长则抑咸效果越好，而补水流量大小则可以通过闸门开度和外排泵站抽排流量来控制。

所述的闸泵群释淡补水压咸调度模型，主要包括河口区二维潮流模块、网河区一维水动力数学模块以及闸泵群联合调控模块三个组成部分。

其中：

（1）河口区二维潮流模块采用贴体正交曲线坐标系下的二维潮流控制方程，形式如下：

连续方程

动量方程

$\frac{\&PartialD;\left(\mathrm{Hu}\right)}{\&PartialD;t}+\frac{1}{C_{\mathrm{\&zeta;}}{C}_{\mathrm{\&eta;}}}\&lsqb;\frac{\&PartialD;}{\&PartialD;\mathrm{\&zeta;}}\left(C_{\mathrm{\&eta;}}\mathrm{Huu}\right)+\frac{\&PartialD;}{\&PartialD;\mathrm{\&eta;}}\left(C_{\mathrm{\&zeta;}}\mathrm{Hvu}\right)+\mathrm{Hvu}\frac{\&PartialD;C_{\mathrm{\&zeta;}}}{\&PartialD;\mathrm{\&eta;}}-H{v}^2\frac{\&PartialD;C_{\mathrm{\&eta;}}}{\&PartialD;\mathrm{\&zeta;}}\&rsqb;+\frac{\mathrm{gu}\sqrt{u^2+v^2}}{C^2}$

$+\frac{\mathrm{gH}}{C_{\mathrm{\&zeta;}}}\frac{\&PartialD;h}{\&PartialD;\mathrm{\&zeta;}}-\mathrm{fvH}=\frac{1}{C_{\mathrm{\&zeta;}}{C}_{\mathrm{\&eta;}}}\&lsqb;\frac{\&PartialD;}{\&PartialD;\mathrm{\&zeta;}}\left(C_{\mathrm{\&eta;}}H{\mathrm{\&sigma;}}_{\mathrm{\&zeta;\&zeta;}}\right)+\frac{\&PartialD;}{\&PartialD;\mathrm{\&eta;}}\left(C_{\mathrm{\&zeta;}}H{\mathrm{\&sigma;}}_{\mathrm{\&zeta;\&eta;}}\right)+H{\mathrm{\&sigma;}}_{\mathrm{\&zeta;\&eta;}}\frac{\&PartialD;C_{\mathrm{\&zeta;}}}{\&PartialD;\mathrm{\&eta;}}-H{\mathrm{\&sigma;}}_{\mathrm{\&eta;\&eta;}}\frac{\&PartialD;C_{\mathrm{\&eta;}}}{\&PartialD;\mathrm{\&zeta;}}\&rsqb;$

$\frac{\&PartialD;\left(\mathrm{Hv}\right)}{\&PartialD;t}+\frac{1}{C_{\mathrm{\&zeta;}}{C}_{\mathrm{\&eta;}}}\&lsqb;\frac{\&PartialD;}{\&PartialD;\mathrm{\&zeta;}}\left(C_{\mathrm{\&eta;}}\mathrm{Huv}\right)+\frac{\&PartialD;}{\&PartialD;\mathrm{\&eta;}}\left(C_{\mathrm{\&zeta;}}\mathrm{Hvv}\right)+\mathrm{Huv}\frac{\&PartialD;C_{\mathrm{\&eta;}}}{\&PartialD;\mathrm{\&zeta;}}-H{u}^2\frac{\&PartialD;C_{\mathrm{\&zeta;}}}{\&PartialD;\mathrm{\&eta;}}]+\frac{\mathrm{gv}\sqrt{u^2+v^2}}{C^2}$

$+\frac{\mathrm{gH}}{C_{\mathrm{\&eta;}}}\frac{\&PartialD;h}{\&PartialD;\mathrm{\&eta;}}+\mathrm{fuH}=\frac{1}{C_{\mathrm{\&zeta;}}{C}_{\mathrm{\&eta;}}}[\frac{\&PartialD;}{\&PartialD;\mathrm{\&zeta;}}\left(C_{\mathrm{\&eta;}}H{\mathrm{\&sigma;}}_{\mathrm{\&zeta;\&eta;}}\right)+\frac{\&PartialD;}{\&PartialD;\mathrm{\&eta;}}\left(C_{\mathrm{\&zeta;}}H{\mathrm{\&sigma;}}_{\mathrm{\&eta;\&eta;}}\right)+H{\mathrm{\&sigma;}}_{\mathrm{\&zeta;\&eta;}}\frac{\&PartialD;C_{\mathrm{\&eta;}}}{\&PartialD;\mathrm{\&zeta;}}-H{\mathrm{\&sigma;}}_{\mathrm{\&zeta;\&zeta;}}\frac{\&PartialD;C_{\mathrm{\&zeta;}}}{\&PartialD;\mathrm{\&eta;}}]$

式中：u、v分别为ζ、η方向的流速分量；h为水位；H为水深；g为重力加速度；f为柯氏力系数；系数C_ζ、C_η如下：

$C_{\mathrm{\&zeta;}}=\sqrt{x_{\mathrm{\&zeta;}}^2+y_{\mathrm{\&zeta;}}^2}$

$C_{\mathrm{\&eta;}}=\sqrt{x_{\mathrm{\&eta;}}^2{y}_{\mathrm{\&eta;}}^2}$ σ_ζζ、σ_ηη、σ_ζη、σ_ηζ为应力项，其表达式如下：

${\mathrm{\&sigma;}}_{\mathrm{\&zeta;\&zeta;}}=2v_t[\frac{1}{C_{\mathrm{\&zeta;}}}\frac{\&PartialD;u}{\&PartialD;\mathrm{\&zeta;}}+\frac{v}{C_{\mathrm{\&zeta;}}{C}_{\mathrm{\&eta;}}}\frac{\&PartialD;C_{\mathrm{\&zeta;}}}{\&PartialD;\mathrm{\&eta;}}]$

${\mathrm{\&sigma;}}_{\mathrm{\&eta;\&eta;}}==2v_t[\frac{1}{C_{\mathrm{\&eta;}}}\frac{\&PartialD;v}{\&PartialD;\mathrm{\&eta;}}+\frac{u}{C_{\mathrm{\&zeta;}}{C}_{\mathrm{\&eta;}}}\frac{\&PartialD;C_{\mathrm{\&eta;}}}{\&PartialD;\mathrm{\&zeta;}}]$

${\mathrm{\&sigma;}}_{\mathrm{\&zeta;\&eta;}}={\mathrm{\&sigma;}}_{\mathrm{\&eta;\&zeta;}}=v_t[\frac{C_{\mathrm{\&eta;}}}{C_{\mathrm{\&zeta;}}}\frac{\&PartialD;}{\&PartialD;\mathrm{\&zeta;}}\left(\frac{v}{C_{\mathrm{\&eta;}}}\right)+\frac{C_{\mathrm{\&zeta;}}}{C_{\mathrm{\&eta;}}}\frac{\&PartialD;}{\&PartialD;\mathrm{\&eta;}}\left(\frac{u}{C_{\mathrm{\&zeta;}}}\right)]$

其中，v_t为紊动粘性系数，即：v_t=au_*H。式中：a为系数，u_*为摩阻流速，H为水深。

（2）网河区一维水动力数学模块采用一维圣维南方程组，方程如下：

连续方程 $B\frac{\&PartialD;Z}{\&PartialD;t}+\frac{\&PartialD;Q}{\&PartialD;x}=q$

动量方程 $\frac{\&PartialD;Q}{\&PartialD;t}+\frac{\&PartialD;}{\&PartialD;x}\left(\mathrm{\&beta;}\frac{Q^2}{A}\right)+\mathrm{gA}(\frac{\&PartialD;Z}{\&PartialD;x}+S_f)+u_{l}q=0$

式中：Z为断面平均水位；Q、A、B分别为断面流量、过水面积、水面宽度；x、t为距离和时间；q为旁侧入流，负值表示流出；β为动量校正系数；g为重力加速度；S_f为摩阻坡降，采用曼宁公式计算，S_f=g/C²，C=h^1/6/n；u_l为单位流程上的侧向出流流速在主流方向的分量。

（3）闸泵群调控模块

闸泵群联合调度的总体目标是以外江典型水文过程条件，根据上游来水和河口咸潮运动情况，合理控制多汊河口联围内外闸泵群的启闭，适时将外江淡水资源引入联围，置换内河涌污水，尽可能蓄积最多的淡水，并于压咸期间释放淡水抑制咸潮上溯，保障取水口取水安全。

目标函数：

$\min C_p=(W_0\&CenterDot;C_0+\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\&Integral;}_{t_0}^{t_1}(Q_{i,t}\&CenterDot;C_{i,t})\mathrm{dt})/(W_0+\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\&Integral;}_{t_0}^{t_1}{Q}_{i,t}\mathrm{dt})$

式中：

C_p——联围内河涌污染物平均浓度，即换水结束时刻t₁时的联围内河涌污染物平均浓度；

W₀——换水起始时刻t₀时的联围内河涌蓄水量；

C₀——换水起始时刻t₀时的联围内河涌污染物平均浓度；

Q_i,t——t时刻第i座进（排）水闸（或泵）的进（排）水流量，以正负号表示流向，正值表示水从外江流入至联围内河涌，负值表示水从联围内河涌流出至外江，i=1，2，…，N，N为闸（或泵）总数；

C_i,t——t时刻第i座进（排）水闸（或泵）的进（排）水污染物浓度。

整个闸泵群联合调度期间均不考虑联围内河涌天然降水量、人类用水取排量，以及污染物排放量。

约束条件：

①换水历时约束

T_r≤T_rmax

②联围内河涌最高最低水位约束

${\underset{\&OverBar;}{Z}}_{k,t}\&le;Z_{k,t}\&le;{\stackrel{\&OverBar;}{Z}}_{k,t}$

③进水闸进水条件约束

$S_{i,t}\&le;\stackrel{\&OverBar;}{S}$

Z_i,外,t>Z_i,内,t

④排水闸排水条件约束

Z_j,外,t<Z_j,内,t

⑤闸（或泵）运行状态持续时间约束

$T_s\&GreaterEqual;{\underset{\&OverBar;}{T}}_s$

⑥闸门启闭速度约束

${\underset{\&OverBar;}{T}}_d\&le;T_d\&le;{\stackrel{\&OverBar;}{T}}_d$

⑦闸（或泵）安全运行条件约束

$\left|{\mathrm{\&Delta;Z}}_{i,t}\right|\&le;\mathrm{\&Delta;}{\stackrel{\&OverBar;}{Z}}_i$

$\left|Q_{i,t}\right|\&le;{\stackrel{\&OverBar;}{Q}}_i$

式中：T_r为换水历时，T_r=t₁-t₀；T_rmax为允许的最长换水历时，越短越好；Z_k,t、

分别为t时刻第k条内河涌（河段）水位、最高限制水位、最低限制水位，

与堤岸高程有关，避免河水漫溢，

取决于景观、航运、生态等最低水位需求；S_i,t、分别为t时刻第i个进水闸（或泵）进水口的含氯度、取水允许的最高含氯度，一般地，水厂取水要求

不高于250mg/L；Z_i,内,t、Z_i,外,t分别为t时刻第i个进水闸内、外水（潮）位；Z_j,内,t、Z_j,外,t分别为t时刻第j个排水闸内、外水（潮）位；T_s、分别为闸（或泵）维持某一特定工作状态（如闸门全开或全关、泵站开启或关闭等）的时长和允许最短时长，

不宜过小，否则闸泵启闭频繁，调度操作难度大；T_d、

分别为闸门启闭单位高度的用时、允许的最小和最大用时，T_d与闸门设计有关，对调度模型水动力计算稳定性和闸门对水流方向的控制都有一定影响；△Z_i,t、

分别为t时刻第i个闸内外水位差、闸门安全运行允许的最大水位差，Q_i,t、

分别为t时刻第i个闸（或泵）过水流量、设计最大过水流量。

实施例基本情况：本发明选择中山市中顺大围作为示范工程实施地点，中顺大围内河网水系密布，区域具有一定的调蓄能力，而且联围内闸、泵众多，为水力调度调控提供了便利的条件；同时由于中顺大围濒临的磨刀门水道咸潮上溯严重，中山和珠海部分地区，尤其是澳门受咸潮影响较大，缓解这些地区的咸潮威胁具有极其重要经济和社会意义。凫洲水闸、东河水闸、西河水闸等众多水闸作为联围干堤上大型或重要水闸工程，无论是从调水过流能力，还是水闸所在的干流位置来看，对调配枯水期淡水资源都十分有利。开展以中顺大围为示范工程区的闸、泵群联合调度，对抑制咸潮保障供水，提高淡水资源利用率具有是十分重要的意义。

中顺大围位于珠江三角洲河网区南部，西濒西江干流磨刀门水道，东傍东海水道、马宁水道、小榄水道。地形上小下大，略呈三角形，总集水面积约709.36km²，是珠江三角洲五大重要堤围之一。因地跨中山、顺德两市，故名中顺大围，包括中山境内的古镇、小榄、东升、横栏、沙溪、大涌、坦背、板芙、港口、沙朗、张家边和石岐城区和顺德的均安。

中顺大围内主干河道有横贯联围中部的岐江河和与之相交的凫洲河、横琴海、中部排灌渠至狮滘河段以及东南部连接磨刀门水道和小榄水道的岐江河。围内有其他河涌140余条，总长约870km，除少数地处五桂山区的溪流是单向流外，其余绝大多数河流均受潮汐影响，是双向流。其他众多大小河涌、排水沟渠与主干河道相互交联，构成水系发达、结构复杂的联围内河网。凫洲河、横琴海、中部排灌渠、狮滘河河道上游接外江，下游和岐江河连接，全长约34km，河面宽度50～250m，水深约2～13m，是一条河床比较浅、河道比较窄的弱感潮河流。凫洲河口有凫洲河水闸，控制上游进入中顺大围的来水。岐江河横贯中山市中部，经城区向东出东河口水闸，汇入横门水道；向西南经渡头、板芙至西河口水闸，汇入磨刀门水道，全长约40km，河面宽80～200m，平均河宽150m，低潮时水深2～3m，可通航300～500吨位船舶，属感潮河段，是双向流动。

实施例技术方案：本发明实施技术通过2012年3月中顺大围闸泵群水质置换调度实际过程进行实施运用。利用外江涨落潮水位过程，根据模型计算确定的各闸门泵站启闭时间，形成内河涌有规则的可控流路，在内河涌污水置换基础上，进行了蓄水和补充调度，并同步进行原型观测，有效改善内河涌水质，提高下游取水口水质保证率。

技术方案的生成是以上游流域水文模型和骨干水库群联合调度计算预测的三角洲控制断面流量过程，以及三角洲河道潮水位变化过程为基础，确定中顺大围各闸外水位边界条件和下游咸界边界条件，通过前述闸泵群联合调度模型，确定调度方案，并分析预测调度实施效果。

1、基于EasyDHM的珠江流域水文预报。本次实施例通过2012年2月-3月的实际雨水情开展流域水文预报，利用基于EasyDHM的珠江流域分布式水文模型，对关键断面梧州、思贤滘进行短期来水预报。预报结果如图所示。

2、流域骨干水库群优化调度。通过骨干水库群优化调度模型，结合水文预报成果作为输入条件，计算得出各水库的调度方案，通过水库调度，得到下游控制断面的流量过程。

3、本次实施例的关键阶段是闸泵群联合调度的具体操作过程。通过控制各闸的开关，实际调度中顺大围内水量，以达到抑制咸潮的目标。图4-图8是本实施例调度期间流量过程和水位过程的示意图。

（1）调度方案

①调度目标

通过中顺大围主要水闸从外江引水冲污(2012年3月6日已完成)，并尽量多蓄淡水，于补水时刻通过西河水闸向外江补水压咸。

②参与调度水闸

中顺大围管理处可直接调度管理的西河闸、东河闸、铺锦闸，其他水闸实行日常调度。

③调度规则

根据历史多年过程和典型年过程，采用经验模型推求出2012年3月梧州、思贤滘预报流量，进行3天的滚动预报，采用日均实时流量数据进行模型校正，通过径流预报为下游闸泵群联合调度提供边界条件。

根据河口一二维联解潮流模型、潮汐表以及各闸外历史同期潮位过程综合分析确定的下游河口地区潮汐以及咸情预测过程。

根据上游骨干水库群联合调度计算预测的调度期三水马口流量，咸界位置处于东河水闸、西河水闸以下，同时根据河口一二维联解潮流模型、潮汐表以及各闸外历史同期潮位过程综合分析确定的闸外潮位预测过程，确定本次蓄淡、补水抑咸调度方案如下：

i)西河闸、东河闸、铺锦闸于3月17日、18日涨潮时段开闸蓄水，其他时段关闸；

ii)西河闸于3月18日落潮期开闸补水，东河闸、铺锦闸关闭或条件许可时开闸进水；

iii)中顺大围其他水闸自行调度运行。

实施例调度过程：2012年3月16日前，完成中顺大围闸泵群联合调度试验的径流、潮位、咸界预测，模型生成并确定最终调度方案；确定磨刀门水道和石岐河原型观测点，给定中顺大围外江三条测船经纬度定点坐标和中顺大围内测船巡测点；

2012年3月16日，调度组进驻中顺大围调度中心，将调度方案注入调度系统，并以调度指令的形式发送至各水闸管理处；测船准备就绪；

2012年3月17日，测量组携测量仪器设备到位，测船到位，进行设备调试和试测，于上午8时同时开展监测；各闸按调度指令开展蓄水调度；

2012年3月18日，5时许东河、西河、铺锦闸开闸进水；20时许西河闸开闸补水；至19日5时许调度完毕，转入中顺大围日常调度；

调度期间，车辆保障组及时将当天所取水样送回至珠科院中心试验室保存、分析；调度试验指导协调组于调度中心、水闸管理处、原型观测点以及中顺大围内各河段处巡查、指导，协调各闸门之间的调度以及调度与观测同步。

实施例调度成果分析：调度实施过程中，位于磨刀门水道西河水闸以下、南镇水厂取水口以上的三条测船每隔一小时进行一次表中底层流速和流向监测，同时于大涨大退时刻采集水样；石岐河内测船来回巡测，按既定采样点采集水样；中顺大围内河涌监测点、水闸内外监测点每五分钟采集一次水文数据，经通信网络系统自动发送至中心存储器，经调度中心系统处理为水位、流量、水量序列，并图形显示在控制中心。

图1～3分别为中顺大围闸泵群联合调度期间东河闸、西河闸、铺锦闸内外水位变化过程。

经过分析，通过东河、西河、铺锦水闸，于涨潮期间（特别是3月18日16时至20时大潮涨潮时段）开闸蓄水，中顺大围内水位显著增加（控制内河水位1.0m以下），其中东河闸内水位由3月17日大潮期的最高0.29m（3月17日21:45）增加至3月18日大潮期的最高0.7m（3月18日20:15），铺锦闸内水位同样由0.34m增加至0.74m，而西河水闸内水位于开闸释淡前最高增至0.55m。在3月18日大潮期间，东河闸利用高潮位从外江持续引水，最大引水流量达710.67m3/s，主要引水时段引水平均流量达342.76m³/s。西河闸开闸释淡补水抑咸，最大补水流量达827.32m³/s，主要排水时段补水平均流量达544.17m³/s，补水时段超过9h。在补水时段，西河水闸对磨刀门水道取水河段累计补水总量达1839.7万m³，平均每小时补水167.25万m³。同时，根据调度期间水样检测分析，石岐河水道沿程总氮、总磷等污染物指标随着调度蓄水增加呈不断下降趋势，内河涌水质改善较为明显。

可见，采用闸泵群联合调度的方法，利用联围内外河道有利的水文条件进行蓄水，并于抑咸期间开闸释淡补充外江淡水在技术上和实践中是可行的，不仅可以提高淡水资源利用率，而且改善了联围内河涌水质，示范工程下游水厂主要取水口的含氯度明显下降，取水保证率能提高到80%以上，有效得保障了供水安全，产生了较好的经济效益、社会效益及环境效益。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (10)

1.一种水库、闸泵群联合抑咸调度方法，其特征在于，基于EasyDHM的珠江流域分布式水文模型、骨干水库群抑咸调度模型和闸泵群释淡补水压咸调度模型，所述水库、闸泵群联合抑咸调度方法，包括下述步骤：

S1、建立基于EasyDHM的珠江流域分布式水文模型；

S2、建立骨干水库群抑咸调度模型，骨干水库群抑咸调度模型包括两个方面：一个是模型的表达形式，一个是模型的求解方法，模型的表达形式包括目标函数和约束条件；

S3、通过流域水文气象部门发布的实时水雨情信息，利用基于EasyDHM的珠江流域分布式水文模型计算出流域关键节点的流量过程，输入到水库群抑咸调度模型中，通过模型得到下游关键断面的流量控制过程；

S4、利用闸泵群释淡补水压咸调度模型确定放水水闸、释淡补水压咸时机和放水压咸流量。

2.根据权利要求1所述的水库闸泵群联合抑咸调度方法，其特征在于，步骤S1中，所述EasyDHM的珠江流域分布式水文模型实现步骤如下：

S21、基础数据准备：

选取DEM数字高程数据；设定最小水道集水面积阈值，生成河网水系栅格图层；首先把该土地利用类型数据转化为网格尺寸为1km×1km栅格数据格式，然后将土地类型转化为EasyDHM产流模型中相应的土地利用分类，最后根据EasyDHM产流模型的土地利用参数推求方法求各项产流参数；

S22、空间单元划分：

采用参数分区拓扑式率定方法，在基于DEM提取的数字流域基础上确定出每个水文站和水库的控制范围，进而划分出各个参数分区和水库分区；

S23、参数率定：

采用Nash效率系数来衡量模型模拟精度的好坏，其方程表达如下：

其中：R2是用于衡量模型模拟实际时间序列上径流量能力的纳什效率系数，纳什效率系数R2的变化范围是负无穷到1，若为1，则表示模拟的和实测的水文过程线达到完美拟合，R2小于0则表示平均实测径流量会是模拟径流量的一个较好预测。

3.根据权利要求1所述水库、闸泵群联合抑咸调度方法，其特征在于，所述的骨干水库群抑咸调度模型的调度目标包括抑咸调度效果、梯级电站发电经济效益、区域经济发展和生态环境，其目标函数描述为：

OBJ_Fun=Q_抑咸+E_发电+Q_航运+Q_生态

其中：OBJ_Fun表示目标函数，Q_抑咸表示最佳抑咸流量，E_发电表示发电量最大，Q_航运表示航道航运流量要求，Q_生态表示河道生态流量要求。

4.根据权利要1所述的水库闸泵群联合抑咸调度方法，其特征在于，步骤S2中，所述目标函数包括发电效益目标、抑咸流量目标和生态环境目标，其具体如下：

发电效益目标：其表达式为

其中，A(m)为出力系数，QD(m,t)为第m个水库t时段的发电流量，H(m,t)为m个水库t时段的平均发电水头，△t为计算步长，M为电站个数，T为总时段数；

抑咸流量目标：在西江梧州站控制断面流量2100m³/s，石角控制断面流量250m³/s，思贤滘控制断面流量2500m³/s；

生态环境目标：西江梧州、北江石角以及西北江三角洲思贤滘三个控制断面的非汛期生态环境流量分别为1800m³/s、250m³/s、2200m³/s。

5.根据权利要1所述的水库、闸泵群联合抑咸调度方法，其特征在于，步骤S2中，所述模型约束条件具体如下：

（1）目标转换约束；

（2）水库水量平衡约束：

V(m,t+1)=V(m,t)+RW(m,t)-W(m,t)-LW(m,t)

W(m,t)=q(m,t)*△t

V(m,t)、RW(m,t)、W(m.,t)、LW(m,t)分别为第m个水库t时段库容、入库水量、出库水量和损失水量；q(m,t)为第m个水库t时段的出库流量；

（3）出库流量约束：

QDmin(m,t)≤QD(m,t)≤QDmax(m,t)

qmin(m,t)≤q(m,t)≤qmax(m,t)

QDmin(m,t)、QDmax(m,t)分别为第m个水库t时段最小、最大允许过机流量；qmin(m,t)、qmax(m,t)分别为第m个水库t时段最小、最大允许出库流量；

（4）出力约束：

Nmin(m,t)≤N(m,t)≤Nmax(m,t)

$\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}N(m,t)\&GreaterEqual;\mathrm{NSUM}\min \left(t\right)$

N(m,t)、Nmin(m,t)、Nmax(m,t)分别为第m个水库t时段出力、允许最小和最大出力；NSUMmin(t)为梯级t时段允许最低总出力；

（5）水库库容约束：

Vmin(m,t)≤V(m,t)≤Vmax(m,t)

Vmin(m,t)、Vmax(m,t)分别为第m个水库t时段允许库容上下限；

（6）河道水量演进约束：

Q(i+1,t+1)=C₀Q(i,t+1)+C₁Q(i,t)+C₂Q(i+1,t)

∑C=1

Q(i,t)为第i个节点t时段的流量；

（7）变量非负约束。

6.根据权利要求1所述的水库、闸泵群联合抑咸调度方法，其特征在于，步骤S2中，所述模型的求解方法具体为：

S51、根据区间来水及控制断面目标流量要求，考虑水量传播和水量损失因素，自下而上推求各库供水约束下限，并结合发电、生态约束拟定初始调度线；

S52、根据给定的天生桥一级计算期末消落水位，采取自上而下、逆时序的方法，推求各库时段最低、最高水位控制线；

S53、按初始调度线自上而下开始顺时序模拟梯级运行过程，计算时段末各库状态和梯级出力，若模拟结果经水位和出力辨识满足要求，则进入下一时段，否则，按一定规则加入反馈修正量，重新模拟时段运行过程，直到满足水位和出力辨识要求；

S54、重复步骤S53逐时段迭代模拟—反馈修正，直到计算期末，完成一轮迭代；

S55、进行目标辨识，若模拟期末天生桥一级水位与给定龙库期末水位之差满足要求，则结束；否则，形成修正量并反馈到输入端，从计算期初重新新一轮迭代，直至期末水位满足要求。

7.根据权利要1所述的水库、闸泵群联合抑咸调度方法，其特征在于，步骤S4中，所述放水水闸的确定是根据补水对象而定，具体是针对取水口释淡补水压下，分析调度期咸界位置、各取水口缺水紧迫程度以及压咸效果，综合确定；所述释淡补水压咸时机是根据补水对象不同时段淡水需求、咸潮上溯强度综合确定；所述放水压咸流量必须能够满足补水目标断面的最小抑咸流量需求，补水流量大小则通过闸门开度和外排泵站抽排流量来控制。

8.根据权利要1所述的水库、闸泵群联合抑咸调度方法，其特征在于，所述的闸泵群释淡补水压咸调度模型包括河口区二维潮流模块、网河区一维水动力数学模块以及闸泵群联合调控模块三个组成部分；

其中：

（1）河口区二维潮流模块采用贴体正交曲线坐标系下的二维潮流控制方程，形式如下：

连续方程

动量方程

$\frac{\&PartialD;\left(\mathrm{Hu}\right)}{\&PartialD;t}+\frac{1}{C_{\mathrm{\&zeta;}}{C}_{\mathrm{\&eta;}}}\&lsqb;\frac{\&PartialD;}{\&PartialD;\mathrm{\&zeta;}}\left(C_{\mathrm{\&eta;}}\mathrm{Huu}\right)+\frac{\&PartialD;}{\&PartialD;\mathrm{\&eta;}}\left(C_{\mathrm{\&zeta;}}\mathrm{Hvu}\right)+\mathrm{Hvu}\frac{\&PartialD;C_{\mathrm{\&zeta;}}}{\&PartialD;\mathrm{\&eta;}}-H{v}^2\frac{\&PartialD;C_{\mathrm{\&eta;}}}{\&PartialD;\mathrm{\&zeta;}}\&rsqb;+\frac{\mathrm{gu}\sqrt{u^2+v^2}}{C^2}$

$+\frac{\mathrm{gH}}{C_{\mathrm{\&zeta;}}}\frac{\&PartialD;h}{\&PartialD;\mathrm{\&zeta;}}-\mathrm{fvH}=\frac{1}{C_{\mathrm{\&zeta;}}{C}_{\mathrm{\&eta;}}}\&lsqb;\frac{\&PartialD;}{\&PartialD;\mathrm{\&zeta;}}\left(C_{\mathrm{\&eta;}}H{\mathrm{\&sigma;}}_{\mathrm{\&zeta;\&zeta;}}\right)+\frac{\&PartialD;}{\&PartialD;\mathrm{\&eta;}}\left(C_{\mathrm{\&zeta;}}H{\mathrm{\&sigma;}}_{\mathrm{\&zeta;\&eta;}}\right)+H{\mathrm{\&sigma;}}_{\mathrm{\&zeta;\&eta;}}\frac{\&PartialD;C_{\mathrm{\&zeta;}}}{\&PartialD;\mathrm{\&eta;}}-H{\mathrm{\&sigma;}}_{\mathrm{\&eta;\&eta;}}\frac{\&PartialD;C_{\mathrm{\&eta;}}}{\&PartialD;\mathrm{\&zeta;}}\&rsqb;$

$\frac{\&PartialD;\left(\mathrm{Hv}\right)}{\&PartialD;t}+\frac{1}{C_{\mathrm{\&zeta;}}{C}_{\mathrm{\&eta;}}}\&lsqb;\frac{\&PartialD;}{\&PartialD;\mathrm{\&zeta;}}\left(C_{\mathrm{\&eta;}}\mathrm{Huv}\right)+\frac{\&PartialD;}{\&PartialD;\mathrm{\&eta;}}\left(C_{\mathrm{\&zeta;}}\mathrm{Hvv}\right)+\mathrm{Huv}\frac{\&PartialD;C_{\mathrm{\&eta;}}}{\&PartialD;\mathrm{\&zeta;}}-H{u}^2\frac{\&PartialD;C_{\mathrm{\&zeta;}}}{\&PartialD;\mathrm{\&eta;}}]+\frac{\mathrm{gv}\sqrt{u^2+v^2}}{C^2}$

$+\frac{\mathrm{gH}}{C_{\mathrm{\&eta;}}}\frac{\&PartialD;h}{\&PartialD;\mathrm{\&eta;}}+\mathrm{fuH}=\frac{1}{C_{\mathrm{\&zeta;}}{C}_{\mathrm{\&eta;}}}[\frac{\&PartialD;}{\&PartialD;\mathrm{\&zeta;}}\left(C_{\mathrm{\&eta;}}H{\mathrm{\&sigma;}}_{\mathrm{\&zeta;\&eta;}}\right)+\frac{\&PartialD;}{\&PartialD;\mathrm{\&eta;}}\left(C_{\mathrm{\&zeta;}}H{\mathrm{\&sigma;}}_{\mathrm{\&eta;\&eta;}}\right)+H{\mathrm{\&sigma;}}_{\mathrm{\&zeta;\&eta;}}\frac{\&PartialD;C_{\mathrm{\&eta;}}}{\&PartialD;\mathrm{\&zeta;}}-H{\mathrm{\&sigma;}}_{\mathrm{\&zeta;\&zeta;}}\frac{\&PartialD;C_{\mathrm{\&zeta;}}}{\&PartialD;\mathrm{\&eta;}}]$

式中：u、v分别为ζ、η方向的流速分量；h为水位；H为水深；g为重力加速度；f为柯氏力系数；系数C_ζ、C_η如下：

$C_{\mathrm{\&zeta;}}=\sqrt{x_{\mathrm{\&zeta;}}^2+y_{\mathrm{\&zeta;}}^2}$

$C_{\mathrm{\&eta;}}=\sqrt{x_{\mathrm{\&eta;}}^2{y}_{\mathrm{\&eta;}}^2}$

σ_ζζ、σ_ηη、σ_ζη、σ_ηζ为应力项，其表达式如下：

${\mathrm{\&sigma;}}_{\mathrm{\&zeta;\&zeta;}}=2v_t[\frac{1}{C_{\mathrm{\&zeta;}}}\frac{\&PartialD;u}{\&PartialD;\mathrm{\&zeta;}}+\frac{v}{C_{\mathrm{\&zeta;}}{C}_{\mathrm{\&eta;}}}\frac{\&PartialD;C_{\mathrm{\&zeta;}}}{\&PartialD;\mathrm{\&eta;}}]$

${\mathrm{\&sigma;}}_{\mathrm{\&eta;\&eta;}}==2v_t[\frac{1}{C_{\mathrm{\&eta;}}}\frac{\&PartialD;v}{\&PartialD;\mathrm{\&eta;}}+\frac{u}{C_{\mathrm{\&zeta;}}{C}_{\mathrm{\&eta;}}}\frac{\&PartialD;C_{\mathrm{\&eta;}}}{\&PartialD;\mathrm{\&zeta;}}]$

${\mathrm{\&sigma;}}_{\mathrm{\&zeta;\&eta;}}={\mathrm{\&sigma;}}_{\mathrm{\&eta;\&zeta;}}=v_t[\frac{C_{\mathrm{\&eta;}}}{C_{\mathrm{\&zeta;}}}\frac{\&PartialD;}{\&PartialD;\mathrm{\&zeta;}}\left(\frac{v}{C_{\mathrm{\&eta;}}}\right)+\frac{C_{\mathrm{\&zeta;}}}{C_{\mathrm{\&eta;}}}\frac{\&PartialD;}{\&PartialD;\mathrm{\&eta;}}\left(\frac{u}{C_{\mathrm{\&zeta;}}}\right)]$

其中，v_t为紊动粘性系数，即：v_t=au_*H，式中：a为系数，u_*为摩阻流速，H为水深；

（2）网河区一维水动力数学模块采用一维圣维南方程组，方程如下：

连续方程 $B\frac{\&PartialD;Z}{\&PartialD;t}+\frac{\&PartialD;Q}{\&PartialD;x}=q$

动量方程 $\frac{\&PartialD;Q}{\&PartialD;t}+\frac{\&PartialD;}{\&PartialD;x}\left(\mathrm{\&beta;}\frac{Q^2}{A}\right)+\mathrm{gA}(\frac{\&PartialD;Z}{\&PartialD;x}+S_f)+u_{l}q=0$

式中：Z为断面平均水位；Q、A、B分别为断面流量、过水面积、水面宽度；x、t为距离和时间；q为旁侧入流，负值表示流出；β为动量校正系数；g为重力加速度；S_f为摩阻坡降，采用曼宁公式计算，S_f=g/C²，C=h^1/6/n；u_l为单位流程上的侧向出流流速在主流方向的分量；

（3）闸泵群调控模块

闸泵群联合调度的总体目标是以外江典型水文过程条件，根据上游来水和河口咸潮运动情况，合理控制多汊河口联围内外闸泵群的启闭，适时将外江淡水资源引入联围，置换内河涌污水，尽可能蓄积最多的淡水，并于压咸期间释放淡水抑制咸潮上溯，保障取水口取水安全。

9.根据权利要1所述的水库、闸泵群联合抑咸调度方法，其特征在于，所述的闸泵群释淡补水压咸调度模型的目标函数为：

$\min C_p=(W_0\&CenterDot;C_0+\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\&Integral;}_{t_0}^{t_1}(Q_{i,t}\&CenterDot;C_{i,t})\mathrm{dt})/(W_0+\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\&Integral;}_{t_0}^{t_1}{Q}_{i,t}\mathrm{dt})$

式中：

C_p是联围内河涌污染物平均浓度，即换水结束时刻t₁时的联围内河涌污染物平均浓度；

W₀是换水起始时刻t₀时的联围内河涌蓄水量；

C₀是换水起始时刻t₀时的联围内河涌污染物平均浓度；

Q_i,t是t时刻第i座进/排水闸或泵的进/排水流量，以正负号表示流向，正值表示水从外江流入至联围内河涌，负值表示水从联围内河涌流出至外江，i=1，2，…，N，N为闸或泵总数；

C_i,t是t时刻第i座进/排水闸或泵的进/排水污染物浓度。

10.根据权利要1所述的水库、闸泵群联合抑咸调度方法，其特征在于，所述的闸泵群释淡补水压咸调度模型的约束条件为：

①换水历时约束

T_r≤T_rmax

②联围内河涌最高最低水位约束

${\underset{\&OverBar;}{Z}}_{k,t}\&le;Z_{k,t}\&le;{\stackrel{\&OverBar;}{Z}}_{k,t}$

③进水闸进水条件约束

$S_{i,t}\&le;\stackrel{\&OverBar;}{S}$

Z_i,外,t>Z_i,内,t

④排水闸排水条件约束

Z_j,外,t<Z_j,内,t

⑤闸或泵运行状态持续时间约束

$T_s\&GreaterEqual;{\underset{\&OverBar;}{T}}_s$

⑥闸门启闭速度约束

${\underset{\&OverBar;}{T}}_d\&le;T_d\&le;{\stackrel{\&OverBar;}{T}}_d$

⑦闸或泵安全运行条件约束

$\left|\mathrm{\&Delta;}{Z}_{i,t}\right|\&le;\mathrm{\&Delta;}{\stackrel{\&OverBar;}{Z}}_i$

$\left|Q_{i,t}\right|\&le;{\stackrel{\&OverBar;}{Q}}_i$

式中：T_r为换水历时，T_r=t₁-t₀；T_rmax为允许的最长换水历时；Z_k,t、

分别为t时刻第k条内河涌（河段）水位、最高限制水位、最低限制水位，S_i,t、

分别为t时刻第i个进水闸或泵进水口的含氯度、取水允许的最高含氯度；Z_i,内,t、Z_i,外,t分别为t时刻第i个进水闸内、外水位；Z_j,内,t、Z_j,外,t分别为t时刻第j个排水闸内、外水位；T_s、

分别为闸或泵维持某一特定工作状态的时长和允许最短时长；T_d、

分别为闸门启闭单位高度的用时、允许的最小和最大用时；△Z_i,t、

分别为t时刻第i个闸内外水位差、闸门安全运行允许的最大水位差，Q_i,t、

分别为t时刻第i个闸或泵过水流量、设计最大过水流量。

Images