CN115456422A - 一种基于计算水动力学的灌区配水计划动态预演修正方法 - Google Patents

Abstract

本发明公开了一种基于计算水动力学的灌区配水计划动态预演修正方法，包括构建基本地理空间矢量数据图层；对构建的基本地理空间矢量数据图层进行网格剖分；构建用于渠系调控模拟计算的双向扩展的渠系数据结构；根据灌溉渠系中的汊点位置、水工建筑物位置，将渠系划分为单渠段的集合；收集灌溉渠系配水计划，对地理信息进行面矢量的空间建模；采用基于标量耗散有限体积法的一维水动力学模型构建电子渠系；利用电子渠系分别对灌溉渠系整体和独立渠段的运行过程进行动态预演修正；采用基于条件触发的模型预测控制算法对动态预演修正结果进行调控优化。本发明更全面的考虑了渠系运行过程中的水流运动物理规律，具有低成本、高解释性的特点。

Description

一种基于计算水动力学的灌区配水计划动态预演修正方法

技术领域

本发明涉及灌溉渠道配水技术领域，具体涉及一种基于计算水动力学的灌区配水计划动态预演修正方法。

背景技术

灌区是一种强人类活动影响区，水资源的时空分布及运动与自然流域显著不同。国内灌区运行管理经验性强、不够科学，现有的基于最优化理论的配水模型对渠系物理约束考虑不足，无法充分发挥水工建筑物的协同调控潜力，导致灌溉渠道配水计划与实施过程割裂开来，供需分离，操作性失水、弃水现象突出，存在输水过程调控不到位，基层管水员劳动强度大，末级量水设备失效等问题。

目前，现有技术主要有两类技术路线：一类是采用运筹学理论构建配水优化编组模型，对水资源时空分配进行优化调度，该方法在可供水量或调度目标发生较大变化时具有较高的参考价值，能够适时调整调度计划的主要发力方向，但仍存在过于概化渠道系统的动力学特征的问题，无法定量表征渠道系统的非线性和迟滞性，不能为调度人员提供精确地决策依据；另一类技术路线是采用更具有可解释性的数值仿真技术对渠道系统调度方案开展测试与评价。然而，由于数值仿真技术在水利方面多应用于洪水演进模拟、工程设计或渠道控制算法性能测试，相关算法与模型并没有针对具有大量闸门、分水口及放水涵的灌溉渠系整体计算进行深入的研发与调教，导致相关研究始终没有摆脱原有应用场景的束缚，无法直接应用于灌溉渠系的实时运行调度。

发明内容

为了克服现有技术的不足之处，针对灌区中存在的渠系复杂、控制性非控制性环境交错分布、各类内边界条件数量众多等特点，本发明提出一种基于计算水动力学的灌溉配水计划动态预演修正方法。能够在配水计划制定阶段完成相对充分的调控工作，克服供需分离问题，形成更符合物理规律、更易于实施、更符合需求的配水计划，挖掘现有工程的协同调控潜力，有效减少渠系损失浪费，提高渠系的运行管理水平，降低基层劳动强度。

为了达到上述发明目的，本发明采用的技术方案为：

一种基于计算水动力学的灌区配水计划动态预演修正方法，包括以下步骤：

S1、构建基本地理空间矢量数据图层，并划分用水单元；

S2、对构建的基本地理空间矢量数据图层进行网格剖分，收集渠系的设计数据和设计运行条件数据，并将其输入网格模型；

S3、构建用于渠系调控模拟计算的双向扩展的渠系数据结构；

S4、根据灌溉渠系中的汊点位置、水工建筑物位置，将渠系划分为单渠段的集合，利用渠系数据结构将渠系组织为渠段集合；

S5、收集灌溉渠系配水计划，对地理信息进行面矢量的空间建模，叠加于步骤S2得到的网格模型，并将用水单元的配水信息输入渠系数据结构；

S6、采用基于标量耗散有限体积法的一维水动力学模型构建电子渠系；

S7、利用电子渠系分别对灌溉渠系整体和独立渠段的运行过程进行动态预演修正；

S8、采用基于条件触发的模型预测控制算法对动态预演修正结果进行调控优化。

可选地，步骤S3具体包括以下分步骤：

S31、采用线段表示每个双向扩展的渠段，将每个分水节点作为渠段端点，并且在相邻渠段端点之间采用弧线箭头表示两个端点之间的关系以及信息交换，构建数据结构；

S32、获取各个渠道名称与属性的对应关系，根据数据结构分类存储各渠道及渠段的信息；

S33、采用汊点处局部节点拓扑关系处理方法，确定各渠道之间的关系，构建双向扩展的渠系数据结构。

可选地，步骤S6具体包括以下分步骤：

S61、利用地表水位和流量为因变量的矩阵形式的一维圣维南方程组描述明渠水流运动过程；

S62、沿河道纵向进行空间离散，并采用标量耗散有限体积法对一维圣维南方程组进行时空离散，得到空间节点上关于地表水位和流量的线性代数方程组。

可选地，所述空间节点上关于地表水位和流量的线性代数方程组具体为：

其中，r_i为时空步长比值，

为第n+1个时间离散节点中空间离散节点i-1在第p次收敛时的系数，

为第n+1个时间离散节点中空间离散节点i-1在第p+1次收敛时的断面面积，

为第n+1个时间离散节点中空间离散节点i在第p次收敛时的系数，

为第n+1个时间离散节点中空间离散节点i在第p次收敛时的单位长度入渗率，ω为实虚步长比值，

Δτ为任意节点处于湿边界状态的持续时间，

为第n+1个时间离散节点中空间离散节点i在第p次收敛时的系数，

为第n+1个时间离散节点中空间离散节点i在第p次收敛时的系数，

为第n+1个时间离散节点中空间离散节点i在第p次收敛时的流量，

为第n个时间离散节点中空间离散节点i在第p次收敛时的单位长度入渗率，

为第n+1个时间离散节点中空间离散节点i在第p+1次收敛时的流量，g为重力加速度，

为第n+1个时间离散节点中空间离散节点i在第p次收敛时的水力半径项，

为第n个时间离散节点中空间离散节点i在第p次收敛时的流量，Δt为时间步长，

为第n+1个时间离散节点中空间离散节点i在第p次收敛时的水位，Δx_i为第i空间离散节点处的空间步长，

为第n+1个时间离散节点中空间离散节点i在第p次收敛时的系数，

为第n个时间离散节点中空间离散节点i在第p次收敛时的流量；

边界条件为：

h＝HIq＝0(t＝0)

q_n＝Q_kn(t_kn<t<t_(k+1)n)

q_η＝Q_kη (t_kη<t<t_(k+1)η)

其中，h为目标节点在目标时刻的水深，HI为最小水深约束，q为目标节点在目标时刻的水深；(t_kn,t_(+1)n)为第n个闸门的第k个配水时段；q_n为第n个闸门的流量值；Q_kn为第n个闸门在时段k内的调控流量值；h_η是自由出流条件下出流边界节点η的地表水深值；q_η为出流边界节点η处的流量值；Q_kη为出流边界节点η处在时段k内的排水流量值；(t_kη,t_(k+1)η)为出流边界节点η的第k个排水时段。

可选地，步骤S7具体包括以下分步骤：

S71、将电子渠系作为控制对象，采用设定采用周期在电子渠系运行过程中持续检测渠道系统的运行状态数据是否满足设定需求；

S72、在渠道运行过程中，渠道水位偏离预定目标时，采用前馈控制方法调节各区段上游闸门过流量；

S73、获取初始入流条件后，开始渠段首部流量的反演计算过程，在每次修正时仅考虑渠段首部流量调控动作开始的初态与调控扰动传递到分水建筑物的末态，利用水动力模型演进渠段内部的动力学过程。

可选地，所述采用前馈控制方法调节各区段上游闸门过流量对应的系统模型具体为：

x(k+1)＝A(k)x(k)+B_u(k)u(k)+B_d(k)d(k)

y(k)＝C(k)x(k)+D(k)u(k)

e(k)＝(y_r(k)-y(k))

其中，x(k+1)为k+1时刻状态向量；A(k)为状态矩阵；x(k)为k时刻状态向量；B_u(k)为控制向量；u(k)为输入向量；B_d(k)为扰动矩阵；d(k)为扰动向量；y(k)为输出向量；C(k)和D(k)分别为观测矩阵和输入输出矩阵；y_r(k)为期望输出向量；e(k)为k时刻系统输出向量y(k)和期望输出向量y_r(k)的偏差；H_aim和Q_aim为需求确定的状态变量；H^k为k时刻的所有节点的水深组成的向量，Q^k为k时刻的所有节点的流量组成的向量。

可选地，所述反演计算过程中每次正向演进计算根据目标分水建筑物处的模拟值与需水值偏差来确定下一次正向演进时渠段首部的调整幅度，表示为

s.t.

x(k+1)＝A(k)x(k)+B_u(k)u(k)+B_d(k)d(k),

u(t)-q_up,k<δsq_up,k

其中，J为描述目标点处流量是否满足目标要求的性能指标；q_n,k为第k个分水建筑物处的目标流量；q_sn,k(t_f)为第k个分水建筑物处在时刻t_f时的模拟流量；x₀为初始条件，即某次正向演进中渠首的恒定流量；(τ_nk)^T为水力响应时间；Δq_nr,k(t_f)为第k个分水建筑物目标流量q_n,k与模拟流量q_sn,k(t_f)之差；τ_foremost为优先级最高的分水建筑物需要渠段上游控制建筑物提前调度的时间；q_up,k为渠段首部现状流量；(0,1)^T为控制向量，表明对流量进行操作；Ψ(x(t_f),t_f)为末态目标集；error为容忍误差；u(t)为输入向量，即控制律，其形式与状态向量相匹配，为(ΔH,ΔQ),这里仅对渠段首部流量进行控制，因此取控制向量为(0,1)^T；p_k为第k个分水建筑物的权重系数；X为沿渠道方向的空间坐标；δs为渠段上游控制建筑物每次动作的最大幅度。

可选地，步骤S8具体包括以下分步骤：

S81、构建调控优化的触发条件；

S82、在每一次模型预测推演前，根据调控目的规划调控信息传输路径；

S83、以需水点模拟流量与需求流量的偏差和各渠段中控制点末态水位与设计水位的偏差作为优化目标，以控制建筑物过流量变化幅度为控制变量，根据需求侧变化对各个控制建筑物进行调控。

可选地，所述根据调控目的规划调控信息传输路径具体为：

首先从需求点开始，逐级向上游检索，同一渠道内沿程已确定需求的分水口作为已知分水参与计算，而不再进行其下级渠道的反演计算；然后将点所属渠道级别相同且相交的渠段需求也作为已知分水参与计算，从需水点至渠道系统首部，逐渠段、逐渠道不断向上游递推，最终确定此次调控优化的传输路径。

可选地，所述根据需求侧变化对各个控制建筑物进行调控表示为：

s.t.

其中，J(k)为调控优化过程中第k时间步的性能指标；N_c为控制点数量；e(k+j)为状态变量变幅向量；

为(k+j)时间步时控制点i处水位的变化幅度；k为调控优化过程中的第k时间步；j为预测时域对应的时间步数，k+j即为预测时域内的末态；Q为n×n维状态加权矩阵，n＝N_c；P为m×m维状态加权矩阵；N_n为用水单元数量；

为(k+j)时间步时用水单元o处流量的变化幅度；

为(k+j)时间步时控制点i处水位；

为(k+j)时间步时用水单元o处流量的变化幅度；δH_min,i和δH_max,i分别为

控制点处水位波动范围的上下限；δQ_min,o和δQ_max,o分别为

用水单元处流量调节量的上下限；H_min,i和H_max,i分别为

控制点处渠道系统运行水位上下限；Q_min,o和Q_max,o分别为

用水单元处闸门过流量能力上下限。

本发明具有以下有益效果：

(1)本发明将渠道系统运行控制仿真技术应用于灌溉渠系调度决策支持，弥补了现有调度决策支持模型对闸控约束下变断面、变流态渠系输配水过程物理机制描述不足的问题。

(2)本发明为了将灌区尺度渠系水动力学仿真模型与模型预测控制算法集成于一体，提出了特征识别算法与定量表征方法，构建了灌溉渠系输配水预演及自适应调控模型，对缝补现有渠道运行控制仿真研究中“仿真”与“控制”之间的割裂关系做出了有益探索。

(3)本发明可以基于高精度水动力学模型进行配水计划预演及动态修正，在配水计划制定阶段完成相对充分的调控工作，更全面的考虑了渠系运行过程中的水流运动物理规律，具有低成本、高解释性的特点。

附图说明

图1为本发明实施例中一种基于计算水动力学的灌区配水计划动态预演修正方法的流程示意图；

图2为本发明实施例中数据结构示意图；

图3为本发明实施例中渠道系统空间分布示意图；

图4为本发明实施例中双向扩展的渠系数据结构示意图；

图5为本发明实施例中灌溉渠系输配水预演及自适应调控模型结构示意图；

图6为本发明实施例中灌溉渠系输配水预演及自适应调控模型计算流程示意图；

图7为本发明实施例中调控信息传输路径示意图；

图8为本发明实施例中渠首供水模拟与实测流量过程及总需水流量过程示意图。

具体实施方式

下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

本发明的技术路线主要包括以下两部分内容：

(1)建立渠系树状结构特征识别与定量表征方法。以灌溉渠系及其运行调度为主要研究对象，分析渠道系统空间拓扑关系，并结合网格剖分技术细节，拟设计一种专用于灌区尺度树状渠系运行控制仿真的特征识别与定量表征方法，将渠道基本特征整理为方便仿真计算的形式，为编制水动力仿真模型与控制模型耦合程序提供基础要素支撑。

(2)构建灌溉渠系输配水预演及自适应调控模型。深入分析渠道系统运行控制规律，拟提出适用于灌溉渠系运行控制仿真的事件触发式模型预测控制算法，并结合专为灌区尺度树状渠系运行控制仿真设计的特征识别和定量表征方法，构建灌溉渠系输配水预演及自适应调控模型。

如图1所示，本发明实施例提供了一种基于计算水动力学的灌区配水计划动态预演修正方法，包括以下步骤S1至S8：

S1、构建基本地理空间矢量数据图层，并划分用水单元；

在本发明的一个可选实施例中，本发明根据灌区地理空间信息，通过地理信息系统软件进行空间建模，得到基本地理空间矢量数据图层，根据实际情况划分用水单元。

S2、对构建的基本地理空间矢量数据图层进行网格剖分，收集渠系的设计数据和设计运行条件数据，并将其输入网格模型；

在本发明的一个可选实施例中，本发明利用预处理软件对基本的地理空间矢量数据图层进行网格划分，收集渠系的设计数据和设计运行条件数据，其中渠系的设计运行条件数包括渠道长度、纵坡、糙率、横纵断面参数、设计水位、设计流量、安全运行范围等，并将上述信息赋予相关网格节点。

S3、构建用于渠系调控模拟计算的双向扩展的渠系数据结构；

在本发明的一个可选实施例中，本发明的步骤S3具体包括以下分步骤：

S31、采用线段表示每个双向扩展的渠段，将每个分水节点作为渠段端点，并且在相邻渠段端点之间采用弧线箭头表示两个端点之间的关系以及信息交换，构建数据结构；

S32、获取各个渠道名称与属性的对应关系，根据数据结构分类存储各渠道及渠段的信息；

S33、采用汊点处局部节点拓扑关系处理方法，确定各渠道之间的关系，构建双向扩展的渠系数据结构。

具体而言，本发明设计一种专用于灌区尺度树状渠系运行控制仿真的特征识别与定量表征方法，将渠道基本特征整理为方便仿真计算的形式，即构建专用于渠系调控模拟计算的双向扩展的渠系数据结构(Double Extension Pool List-DEPL)，具体步骤如下：

S31、确定调控模拟计算对数据结构的需求；

调控模拟计算的需求包括：

①能够将具有空间属性的渠段组织为能够利用空间矢量数据快速定位的渠段集合。

②能够根据节点拓扑关系索引计算域内的渠道及建筑物参数。

③尽量避免遍历，能够与数值计算过程快速交互，且便于并行化处理。

S32、根据调控模拟计算需求，结合代数拓扑学，构建专用的数据结构，以承担复杂约束下的渠系调控模拟计算的数据支撑工作。

如图2所示，为数据结构示意图。其中p_u、p_c、p_dw、p_b分别代表上游方向渠段、当前渠段、下游方向渠段、分支渠段；相邻节点间的弧线箭头代表其之间的关系以及信息交换(箭头所指节点属于箭头另一侧节点。

每个双向渠段保存如下对象：

①上游方向(Upward)渠段下端端点；

②下游方向(Downward)渠段上端端点；

③分支(Branch)渠段上端端点；

④目标渠段(Current)内节点。

该数据结构的主要思想是：①每个分水节点都为渠段端点；②每个渠段都是双向扩展的，其端点之间都存在对应关系。

为了进一步介绍数据结构，考虑某一简单的渠道系统，干渠部分的蓝红相间部分表示渠道存在变断面情况(即Q1,Q2,Q3分别代表具有不同横断面参数的渠段)，其空间分布如图3所示。

上述渠道系统中各渠道名称与属性的对应情况如表1所示。

表1某渠道系统属性表

S33、当我们获得表1的数据时，无法得知各个渠道之间的隶属关系，仅能初步区分其层级关系，此时需要结合数据结构来分类存储各渠道及渠段的信息，数据结构的树状表示如图4所示，采用类似于前述的汊点处局部节点拓扑关系处理方法，确定各渠道之间的关系。完成该步骤后，无论哪一个节点在计算过程中触发请求，DEPL都能够快速定位，再借助合适的方法规划调控信息传输路径即可对数字空间中的渠道系统进行自适应调控。

S4、根据灌溉渠系中的汊点位置、水工建筑物位置，将渠系划分为单渠段的集合，利用渠系数据结构将渠系组织为渠段集合；

S5、收集灌溉渠系配水计划，对地理信息进行面矢量的空间建模，叠加于步骤S2得到的网格模型，并将用水单元的配水信息输入渠系数据结构；

在本发明的一个可选实施例中，本发明收集灌溉渠系配水计划(包括用水单元的位置、面积、需水量、需水时间、停水时间、设计引水流量等等)，利用地理信息系统软件进行面矢量的空间建模，叠加于2)中的网格模型，并将用水单元的配水信息输入渠系数据结构(DEPL)。

S6、采用基于标量耗散有限体积法的一维水动力学模型构建电子渠系；

在本发明的一个可选实施例中，本发明的步骤S6具体包括以下分步骤：

S61、利用地表水位和流量为因变量的矩阵形式的一维圣维南方程组描述明渠水流运动过程；

S62、沿河道纵向进行空间离散，并采用标量耗散有限体积法对一维圣维南方程组进行时空离散，得到空间节点上关于地表水位和流量的线性代数方程组。

具体而言，本发明从灌溉渠系配水的需求出发，选用基于标量耗散有限体积法的一维水动力学模型作为电子渠系，具体步骤如下：

S61、基于断面流速平均假设，圣维南方程组(Saint-Venant Equations-SVE)能够合理地描述渠道水流运动过程，在实际工程中的应用也非常广泛，本文取水位ξ和流量Q为因变量的矩阵形式的一维圣维南方程组描述明渠水流运动过程，控制方程如下：

其中，U为因变量向量；F为对流通量；S_ξ为水面坡度项；S_f为渠/沟底摩阻项；S_in为入渗项。

上式中向量的解析表达式如下：

其中，A为断面面积；Q为流量；u为流速；g为重力加速度；ξ＝h+z_b为水位；h为水深；x为空间坐标；n为曼宁糙率系数；R为水力半径；I为单位长度入渗率(I＝i_cP)；i_c是地表入渗率。

S62、沿河道纵向进行空间离散，并采用标量耗散有限体积法对一维圣维南方程组进行时空离散，得到空间节点上关于地表水位ξ和流量Q的线性代数方程组，表示为：

其中，r_i为时空步长比值，

为第n+1个时间离散节点中空间离散节点i-1在第p次收敛时的系数，

为第n+1个时间离散节点中空间离散节点i-1在第p+1次收敛时的断面面积，

为第n+1个时间离散节点中空间离散节点i在第p次收敛时的系数，

为第n+1个时间离散节点中空间离散节点i在第p次收敛时的单位长度入渗率，ω为实虚步长比值，

Δτ为任意节点处于湿边界状态的持续时间，

为第n+1个时间离散节点中空间离散节点i在第p次收敛时的系数，

为第n+1个时间离散节点中空间离散节点i在第p次收敛时的系数，

为第n+1个时间离散节点中空间离散节点i在第p次收敛时的流量，

为第n个时间离散节点中空间离散节点i在第p次收敛时的单位长度入渗率，

为第n+1个时间离散节点中空间离散节点i在第p+1次收敛时的流量，g为重力加速度，

为第n+1个时间离散节点中空间离散节点i在第p次收敛时的水力半径项，

为第n个时间离散节点中空间离散节点i在第p次收敛时的流量，Δt为时间步长，

为第n+1个时间离散节点中空间离散节点i在第p次收敛时的水位，Δx_i为第i空间离散节点处的空间步长，

为第n+1个时间离散节点中空间离散节点i在第p次收敛时的系数，

为第n个时间离散节点中空间离散节点i在第p次收敛时的流量；

其中系数α_i、β_i、γ_i、η_i、η_mass,i、η_momentum,i的表达式如下：

边界条件设置如下：

初始条件：

h＝HIq＝0(t＝0)

入流条件及内边界条件：

q_n＝Q_kn(t_kn<t<t_(k+1)n)

出流条件：

q_η＝Q_kη (t_kη<t<t_(k+1)η)

其中，h为目标节点在目标时刻的水深，HI为最小水深约束，q为目标节点在目标时刻的水深；(t_kn,t_(k+1)n)为第n个闸门的第k个配水时段；q_n为第n个闸门的流量值；Q_kn为第n个闸门在时段k内的调控流量值；h_η是自由出流条件下出流边界节点η的地表水深值；q_η为出流边界节点η处的流量值；Q_kη为出流边界节点η处在时段k内的排水流量值；(t_kη,t_(k+1)η)为出流边界节点η的第k个排水时段。

S7、利用电子渠系分别对灌溉渠系整体和独立渠段的运行过程进行动态预演修正；

在本发明的一个可选实施例中，本发明的步骤S7具体包括以下分步骤：

S71、将电子渠系作为控制对象，采用设定采用周期在电子渠系运行过程中持续检测渠道系统的运行状态数据是否满足设定需求；

S72、在渠道运行过程中，渠道水位偏离预定目标时，采用前馈控制方法调节各区段上游闸门过流量；

S73、获取初始入流条件后，开始渠段首部流量的反演计算过程，在每次修正时仅考虑渠段首部流量调控动作开始的初态与调控扰动传递到分水建筑物的末态，利用水动力模型演进渠段内部的动力学过程。

具体而言，本发明为实施渠系运行过程的动态预演修正，利用两个独立的水动力学模型分别承担灌溉渠系整体计算和独立渠段的迭代反演计算，均采用圣维南方程描述。

S71、对于单渠段系统，取t^k＝kT为采样周期(T由时间步长和需求确定)，在电子渠系运行过程中持续监测其是否满足相关需求。

S72、渠道运行过程中，渠道水位偏离预定目标时，需要对相关渠段的上游端节制闸施加适宜的控制动作，使渠道系统在保证渠道安全稳定的前提下满足下游用水需求。

本发明采用前馈控制方法调节各区段上游闸门过流量，相应的系统模型表达式如下：

x(k+1)＝A(k)x(k)+B_u(k)u(k)+B_d(k)d(k)

y(k)＝C(k)x(k)+D(k)u(k)

e(k)＝(y_r(k)-y(k))

其中，x(k+1)为k+1时刻状态向量；A(k)为状态矩阵；x(k)为k时刻状态向量；B_u(k)为控制向量；u(k)为输入向量；B_d(k)为扰动矩阵；d(k)为扰动向量；y(k)为输出向量；C(k)和D(k)分别为观测矩阵和输入输出矩阵；y_r(k)为期望输出向量；e(k)为k时刻系统输出向量y(k)和期望输出向量y_r(k)的偏差；H_aim和Q_aim为需求确定的状态变量；H^k为k时刻的所有节点的水深组成的向量，Q^k为k时刻的所有节点的流量组成的向量。

S73、获取初始入流条件后，即可开始渠段首部流量的反演计算过程。为了提高计算效率和稳定性，反演计算在每次修正时仅考虑渠段首部流量调控动作开始的初态与调控扰动传递到分水建筑物的末态(假设初态与末态的渠道系统皆处于稳定状态，判断渠道系统达到稳态的条件需要根据具体目标渠段确定)，利用水动力模型演进渠段内部的动力学过程，以考虑渠道系统的迟滞性、非线性特点。反演过程中涉及到多次系统模型的正向演进计算，每次正向演进计算根据目标分水建筑物处的模拟值与需水值偏差来确定下一次正向演进时渠段首部的调整幅度，表示为

s.t.

x(k+1)＝A(k)x(k)+B_u(k)u(k)+B_d(k)d(k),

u(t)-q_up,k<δsq_up,k

其中，J为描述目标点处流量是否满足目标要求的性能指标；q_n,k为第k个分水建筑物处的目标流量；q_sn,k(t_f)为第k个分水建筑物处在时刻t_f时的模拟流量；x₀为初始条件，即某次正向演进中渠首的恒定流量；(τ_nk)^T为水力响应时间；Δq_nr,k(t_f)为第k个分水建筑物目标流量q_n,k与模拟流量q_sn,k(t_f)之差；τ_foremost为优先级最高的分水建筑物需要渠段上游控制建筑物提前调度的时间；q_up,k为渠段首部现状流量；(0,1)^T为控制向量，表明对流量进行操作；Ψ(x(t_f),t_f)为末态目标集；error为容忍误差；u(t)为输入向量，即控制律，其形式与状态向量相匹配，为(ΔH,ΔQ),这里仅对渠段首部流量进行控制，因此取控制向量为(0,1)^T；p_k为第k个分水建筑物的权重系数；X为沿渠道方向的空间坐标；δs为渠段上游控制建筑物每次动作的最大幅度。

S8、采用基于条件触发的模型预测控制算法对动态预演修正结果进行调控优化。

在本发明的一个可选实施例中，本发明的S8具体包括以下分步骤：

S81、构建调控优化的触发条件；

S82、在每一次模型预测推演前，根据调控目的规划调控信息传输路径；

S83、以需水点模拟流量与需求流量的偏差和各渠段中控制点末态水位与设计水位的偏差作为优化目标，以控制建筑物过流量变化幅度为控制变量，根据需求侧变化对各个控制建筑物进行调控。

具体而言，渠道系统运行控制仿真与实际调度工作中的逻辑相匹配，包括基本的水量平衡原理、调度单元的划分、渠道运行和控制概念、基本的运行要求、供求变化等等。为了对渠道系统开展清晰的调控工作，常将灌溉系统分解为多个单元(子系统)，分别进行平衡性分析以及自动控制。本研究始终以整体性为落脚点，采用模型预测控制算法(MPC)，并利用仿真模型持续地对灌溉渠系整体的运行过程进行模拟，在数字空间实现对渠道系统运行全过程的预警、预演、预案，模型架构如图5所示。

本发明旨在提出一种基于高效、高精度灌区尺度渠系水动力学仿真模型的自适应调控模型，以满足灌区调度人员在实际工作中提出的定时、定量的决策支持需求。以前述提出的树状渠系处理方法和调控方法为基础，结合MPC对其进行组织与拓展，构建灌溉渠系输配水预演及自适应调控模型(Self-adaption Flow Routing Model–SAFRM)，该模型能够根据仿真模型的实时模拟结果计算并实施调控动作，并施加于仿真模型，其计算流程如图6所示.

在本发明中，针对灌溉渠系调度决策问题对MPC算法进行了修改，将周期性滚动优化过程改为有条件的启动滚动优化，并在洪水条件下考虑启用泄水闸。构建的调控优化的触发条件包括：

条件一：在纯仿真环境下，所有取水以及外生变量皆为已知，因此可以提前对各个扰动计算其需要的调度时间T_pre。在计算前可设置配水初期持续时间T_ini，当T_pre<T_ini，模型忽略骤升、骤降的用水需求，持续保证大流量灌溉；当T_pre>T_ini时，根据渠道系统仿真结果实时进行调控计算，得到最优控制律u^*(t)，使渠道系统从已知初态x₀(t)过渡到控制目标要求的末态x(t_f)，并使目标函数minJ达到预期范围。

条件二：实时监测各单元的质量守恒情况，若单元内水量不守恒(如输入降雨、泄洪等情况)，开始周期性执行步骤(1)和(2)，并允许启动泄水闸。

本发明构造的事件触发式算法是一种改编的模型预测控制算法(CentralizedModel Predictive Control-MPC)。算法的主要目的是尽量减少无谓的模型预测推演，提高对数字空间渠道系统运行过程的调控效率，提高整体计算速度。为了实现这个目的，我们需要在每一次模型预测推演前根据调控目的，规划好主要的调控信息传输路径，也即根据需求确定主要的控制对象。一般来说，其主要内容为：(1)从需求点开始，逐级向上游检索，同一渠道内沿程已确定需求的分水口作为已知分水参与计算，而不再进行其下级渠道的反演计算(锁定)；(2)点所属渠道级别相同且相交的渠段需求也作为已知分水参与计算。从需水点至渠道系统首部，在该原则指导下逐渠段、逐渠道不断向上游递推，最终可确定此次调控信息的传输路径。为了实现算法内容(1)定功能，结合优化算法基本理论中的“剪枝”思想，提出了“锁”机制，将已经过反演计算且不属于此次调控计算重点的渠道简化为出流边界。

为了进一步解释事件触发式算法，参照如图7所示的渠道系统，说明完整的调控计算过程。假设此时“斗渠3-Q8”提出了需求，其斗口即为“需求点”。此时，“需求点”下游的所有渠道皆被锁定，相应的下游需水被简化为出流边界，汇于渠段端点5。另一个方向，从“需求点”开始向上游递推同时进行反演计算，当检索到“支渠-Q5”与“干渠-Q3”交汇处(渠段端点4)时，依据“锁”机制，将“干渠-Q3”下游部分简化为出流边界。完成上述步骤，即可开始“干渠-Q2”反演计算，其下游端点4的初始条件由简化的“干渠-Q3”出流边界和渠段端点4和5之间的反演结果共同确定。同理，渠段端点2和3处也用类似的方法处理，最终可以反演得到渠段端点1的入流边界数据。

为了提高配水计划合理性并充分满足灌区各用水单元需求，本发明将渠道系统运行条件作为优化问题的约束条件，如建筑物过流能力以及渠道运行安全范围。为了在模型中实现对灌溉渠系调度工作的决策支持，建立了相应的目标函数，持续在每一次调控计算中构建相应于当前状态的优化问题。

本发明选择流量作为控制目标，目标函数的主要评价对象是：(1)需水点模拟流量与需求流量的偏差；(2)各渠段中控制点末态水位与设计水位的偏差；以控制建筑物过流量变化幅度为控制变量，根据需求侧变化对各个控制建筑物进行调控；约束条件分类如下：(1)控制变量范围约束：水工建筑物的过流能力范围；(2)状态变量范围约束：渠段全局安全运行范围。

将调控优化过程中第k时间步的优化问题表示为：

s.t.

其中，J(k)为调控优化过程中第k时间步的性能指标；N_c为控制点数量；e(k+j)为状态变量变幅向量；

为(k+j)时间步时控制点i处水位的变化幅度；k为调控优化过程中的第k时间步；j为预测时域对应的时间步数，k+j即为预测时域内的末态；Q为n×n维状态加权矩阵，n＝N_c；P为m×m维状态加权矩阵；N_n为用水单元数量；

为(k+j)时间步时用水单元o处流量的变化幅度；

为(k+j)时间步时控制点i处水位；

为(k+j)时间步时用水单元o处流量的变化幅度；δH_min,i和δH_max,i分别为

控制点处水位波动范围的上下限；δQ_min,o和δQ_max,o分别为

用水单元处流量调节量的上下限；H_min,i和H_max,i分别为

控制点处渠道系统运行水位上下限；Q_min,o和Q_max,o分别为

用水单元处闸门过流量能力上下限。

下面以具体案例对本发明的方案进行说明。

以山西大禹渡灌区南干渠作为研究区域，该渠道系统包含3个渠段，52个支(斗)分水口，研究对象时间段内25个分水口有灌水任务，相关的配水数据如表1所示。

具体步骤包括：

利用步骤S1构建基本地理空间矢量数据图层，根据实际情况划分用水单元，这里以分水口为单位进行用水单元划分；

利用步骤S2对矢量数据图层进行网格剖分，收集设计数据以及运行条件数据，并将其输入网格模型；

步骤S3构建的灌溉渠系调控模拟计算的专用数据结构(DEPL)具有泛化性，可通用于各种树状渠系，因此不必重新构建，直接进行步骤S4和S5，利用(DEPL)将渠系组织为渠段集合；

步骤S6构建的电子渠系组织方法具有泛化性，可通用，在网格模型之上即可进行研究区域的水动力学实时计算；

执行步骤S7和S8，将电子渠系作为控制对象，在求解圣维南方程的时间循环中持续获取渠系运行数据，当无法满足渠系运行需求或即将满足某用水单元的用水需求时，逐级对渠段进行迭代反演计算，并将计算结果作为调控动作施加于电子渠系。

基于上述步骤，对山西大禹渡灌区南干渠某年春浇部分时段输配水过程进行了仿真，大禹渡灌区的调度管理方式主要是根据用水户代表提出的用水需求制定配水计划，当需求出现频繁的变化时，根据管理经验对南干渠首部流量的短时间突变进行过滤，在保证用水需求得到满足的同时尽量减少渠首闸门调控次数。模型计算得到的南干首部流量过程与实际情况较为契合，平均相对误差为8.5％，且在配水初期(TIME∈(0,20000))提供了持续时间较长、流量值较大的调度方案，以减少渠道系统运行过程中的渗漏损失；在配水中前期(TIME∈(300,850))过滤了多次浮动的用水需求，持续保证配水中前期的大流量灌溉状态。由此可知模型能够根据用水需求调控水源供水过程，在用水需求较为旺盛的前期尽量给出了较高的配水流量以减少渗漏损失。

为了进一步评价模型计算结果，下面对渠首调度过程中各控制建筑物及控制点的水位变化进行分析。由图8可知，此次春浇的供水流量较小，因此渠道始终保持低水位运行，模型在整个仿真过程中不断地根据需求流量过程调节控制建筑物过流量，各渠段的控制点及节制闸前水位也随着调控过程不断变化。大部分时间中，各控制点及节制闸前水位始终处于最低水位和最高水位之间。

本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

本发明中应用了具体实施例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (10)

1.一种基于计算水动力学的灌区配水计划动态预演修正方法，其特征在于，包括以下步骤：

S1、构建基本地理空间矢量数据图层，并划分用水单元；

S2、对构建的基本地理空间矢量数据图层进行网格剖分，收集渠系的设计数据和设计运行条件数据，并将其输入网格模型；

S3、构建用于渠系调控模拟计算的双向扩展的渠系数据结构；

S4、根据灌溉渠系中的汊点位置、水工建筑物位置，将渠系划分为单渠段的集合，利用渠系数据结构将渠系组织为渠段集合；

S5、收集灌溉渠系配水计划，对地理信息进行面矢量的空间建模，叠加于步骤S2得到的网格模型，并将用水单元的配水信息输入渠系数据结构；

S6、采用基于标量耗散有限体积法的一维水动力学模型构建电子渠系；

S7、利用电子渠系分别对灌溉渠系整体和独立渠段的运行过程进行动态预演修正；

S8、采用基于条件触发的模型预测控制算法对动态预演修正结果进行调控优化。

2.根据权利要求1所述的一种基于计算水动力学的灌区配水计划动态预演修正方法，其特征在于，步骤S3具体包括以下分步骤：

S31、采用线段表示每个双向扩展的渠段，将每个分水节点作为渠段端点，并且在相邻渠段端点之间采用弧线箭头表示两个端点之间的关系以及信息交换，构建数据结构；

S32、获取各个渠道名称与属性的对应关系，根据数据结构分类存储各渠道及渠段的信息；

S33、采用汊点处局部节点拓扑关系处理方法，确定各渠道之间的关系，构建双向扩展的渠系数据结构。

3.根据权利要求1所述的一种基于计算水动力学的灌区配水计划动态预演修正方法，其特征在于，步骤S6具体包括以下分步骤：

S61、利用地表水位和流量为因变量的矩阵形式的一维圣维南方程组描述明渠水流运动过程；

S62、沿河道纵向进行空间离散，并采用标量耗散有限体积法对一维圣维南方程组进行时空离散，得到空间节点上关于地表水位和流量的线性代数方程组。

4.根据权利要求1所述的一种基于计算水动力学的灌区配水计划动态预演修正方法，其特征在于，所述空间节点上关于地表水位和流量的线性代数方程组具体为：

其中，r_i为时空步长比值，

为第n+1个时间离散节点中空间离散节点i-1在第p次收敛时的系数，

为第n+1个时间离散节点中空间离散节点i-1在第p+1次收敛时的断面面积，

为第n+1个时间离散节点中空间离散节点i在第p次收敛时的系数，

为第n+1个时间离散节点中空间离散节点i在第p次收敛时的单位长度入渗率，ω为实虚步长比值，

Δτ为任意节点处于湿边界状态的持续时间，

为第n+1个时间离散节点中空间离散节点i在第p次收敛时的系数，

为第n+1个时间离散节点中空间离散节点i在第p次收敛时的系数，

为第n+1个时间离散节点中空间离散节点i在第p次收敛时的流量，

为第n个时间离散节点中空间离散节点i在第p次收敛时的单位长度入渗率，

为第n+1个时间离散节点中空间离散节点i在第p+1次收敛时的流量，g为重力加速度，

为第n+1个时间离散节点中空间离散节点i在第p次收敛时的水力半径项，

为第n个时间离散节点中空间离散节点i在第p次收敛时的流量，Δt为时间步长，

为第n+1个时间离散节点中空间离散节点i在第p次收敛时的水位，Δx_i为第i空间离散节点处的空间步长，

为第n+1个时间离散节点中空间离散节点i在第p次收敛时的系数，

为第n个时间离散节点中空间离散节点i在第p次收敛时的流量；

边界条件为：

h＝HIq＝0(t＝0)

q_n＝Q_kn(t_kn<t<t_(k+1)n)

q_η＝Q_kη(t_kη<t<t_(k+1)η)

其中，h为目标节点在目标时刻的水深，HI为最小水深约束，q为目标节点在目标时刻的水深；(t_kn，t_(k+1)n)为第n个闸门的第k个配水时段；q_n为第n个闸门的流量值；Q_kn为第n个闸门在时段k内的调控流量值；h_η是自由出流条件下出流边界节点η的地表水深值；q_η为出流边界节点η处的流量值；Q_kη为出流边界节点η处在时段k内的排水流量值；(t_kη，t_(k+1)η)为出流边界节点η的第k个排水时段。

5.根据权利要求1所述的一种基于计算水动力学的灌区配水计划动态预演修正方法，其特征在于，步骤S7具体包括以下分步骤：

S71、将电子渠系作为控制对象，采用设定采用周期在电子渠系运行过程中持续检测渠道系统的运行状态数据是否满足设定需求；

S72、在渠道运行过程中，渠道水位偏离预定目标时，采用前馈控制方法调节各区段上游闸门过流量；

S73、获取初始入流条件后，开始渠段首部流量的反演计算过程，在每次修正时仅考虑渠段首部流量调控动作开始的初态与调控扰动传递到分水建筑物的末态，利用水动力模型演进渠段内部的动力学过程。

6.根据权利要求5所述的一种基于计算水动力学的灌区配水计划动态预演修正方法，其特征在于，所述采用前馈控制方法调节各区段上游闸门过流量对应的系统模型具体为：

x(k+1)＝A(k)x(k)+B_u(k)u(k)+B_d(k)d(k)

y(k)＝C(k)x(k)+D(k)u(k)

e(k)＝(y_r(k)-y(k))

其中，x(k+1)为k+1时刻状态向量；A(k)为状态矩阵；x(k)为k时刻状态向量；B_u(k)为控制向量；u(k)为输入向量；B_d(k)为扰动矩阵；d(k)为扰动向量；y(k)为输出向量；C(k)和D(k)分别为观测矩阵和输入输出矩阵；y_r(k)为期望输出向量；e(k)为k时刻系统输出向量y(k)和期望输出向量y_r(k)的偏差；H_aim和Q_aim为需求确定的状态变量；H^k为k时刻的所有节点的水深组成的向量，Q^k为k时刻的所有节点的流量组成的向量。

7.根据权利要求6所述的一种基于计算水动力学的灌区配水计划动态预演修正方法，其特征在于，所述反演计算过程中每次正向演进计算根据目标分水建筑物处的模拟值与需水值偏差来确定下一次正向演进时渠段首部的调整幅度，表示为

s.t.

x(k+1)＝A(k)x(k)+B_u(k)u(k)+B_d(k)d(k)，

u(t)-q_up，k＜δsq_up，k

其中，J为描述目标点处流量是否满足目标要求的性能指标；q_n，k为第k个分水建筑物处的目标流量；q_sn，k(t_f)为第k个分水建筑物处在时刻t_f时的模拟流量；x₀为初始条件，即某次正向演进中渠首的恒定流量；(τ_nk)^T为水力响应时间；Δq_nr，k(t_f)为第k个分水建筑物目标流量q_n，k与模拟流量q_sn，k(t_f)之差；τ_formost为优先级最高的分水建筑物需要渠段上游控制建筑物提前调度的时间；q_up，k为渠段首部现状流量；(0，1)^T为控制向量，表明对流量进行操作；Ψ(x(t_f)，t_f)为末态目标集；error为容忍误差；u(t)为输入向量，即控制律，其形式与状态向量相匹配，为(ΔH，ΔQ)，这里仅对渠段首部流量进行控制，因此取控制向量为(0，1)^T；p_k为第k个分水建筑物的权重系数；X为沿渠道方向的空间坐标；δs为渠段上游控制建筑物每次动作的最大幅度。

8.根据权利要求1所述的一种基于计算水动力学的灌区配水计划动态预演修正方法，其特征在于，步骤S8具体包括以下分步骤：

S81、构建调控优化的触发条件；

S82、在每一次模型预测推演前，根据调控目的规划调控信息传输路径；

S83、以需水点模拟流量与需求流量的偏差和各渠段中控制点末态水位与设计水位的偏差作为优化目标，以控制建筑物过流量变化幅度为控制变量，根据需求侧变化对各个控制建筑物进行调控。

9.根据权利要求8所述的一种基于计算水动力学的灌区配水计划动态预演修正方法，其特征在于，所述根据调控目的规划调控信息传输路径具体为：

首先从需求点开始，逐级向上游检索，同一渠道内沿程已确定需求的分水口作为已知分水参与计算，而不再进行其下级渠道的反演计算；然后将点所属渠道级别相同且相交的渠段需求也作为已知分水参与计算，从需水点至渠道系统首部，逐渠段、逐渠道不断向上游递推，最终确定此次调控优化的传输路径。

10.根据权利要求1所述的一种基于计算水动力学的灌区配水计划动态预演修正方法，其特征在于，所述根据需求侧变化对各个控制建筑物进行调控表示为：

e(k+j)＝x(k+j)-x(k+j-1)

s.t.

其中，J(k)为调控优化过程中第k时间步的性能指标；N_c为控制点数量；e(k+j)为状态变量变幅向量；

为(k+j)时间步时控制点i处水位的变化幅度；k为调控优化过程中的第k时间步；j为预测时域对应的时间步数，k+j即为预测时域内的末态；Q为n×n维状态加权矩阵，n＝N_c；P为m×m维状态加权矩阵；N_n为用水单元数量；

为(k+j)时间步时用水单元0处流量的变化幅度；

为(k+j)时间步时控制点i处水位；

为(k+j)时间步时用水单元o处流量的变化幅度；δH_min，i和δH_max，i分别为

控制点处水位波动范围的上下限；δQ_min，o和δQ_max，o分别为

用水单元处流量调节量的上下限；H_min，i和H_max，i分别为

控制点处渠道系统运行水位上下限；Q_min，o和Q_max，o分别为

用水单元处闸门过流量能力上下限。

Images