CN107633141A - 一种一维水动力模型泵站特性曲线系数的辨识方法 - Google Patents

Abstract

本发明公开了一种一维水动力模型泵站特性曲线系数的辨识方法，涉及水动力模型参数识别技术领域。方法包括：基于圣维南方程，建立目标明渠的一维明渠水动力模型，并率定待辨识泵站前、后连接渠道的糙率；待辨识泵站边界以连续方程和泵站特性曲线为相容方程，与S1中的圣维南方程耦合，形成包含待辨识泵站边界的一维渠道‑泵站水动力模型；构建以待辨识泵站的三个泵站特性曲线系数为状态向量的状态空间方程包括状态方程：Xⁿ⁺¹＝ΦⁿXⁿ+wⁿ⁺¹和量测方程：Yⁿ⁺¹＝HXⁿ⁺¹+vⁿ⁺¹；选取泵站机组角度或转速不变时段的观测数据，调用集合卡尔曼滤波器，同化曲线系数。本发明方法自动辨识任意工况下的曲线系数，有效提高河网水动力模型的数值模拟的精度和适应性。

Description

一种一维水动力模型泵站特性曲线系数的辨识方法

技术领域

本发明涉及水动力模型参数识别技术领域，尤其涉及一种一维水动力模型泵站特性曲线系数的辨识方法。

背景技术

泵站是常见的渠系水工建筑物之一，其在灌溉排涝、跨流域调水、生产生活供水等方面发挥着重要作用。泵站特性曲线直接反映了泵站各相关水力要素的联系，是一维水动力模型中常用的泵站边界相容条件。因此，精确辨识泵站特性曲线系数，对提高一维明渠水动力模型的精度具有极其重要的价值。

在包含泵站的河渠如平原河网、梯级泵站调水工程等一维水动力模拟中，水泵特性曲线系数常由模型水泵试验得到。但是，由于模型水泵和原型水泵之间存在比尺效应，导致水泵特性曲线系数与其真值存在差异。又由于泵站特性曲线随水泵机组开关、运行工况的变化而变化，而导致水泵特性曲线系数具有时变性。目前，现有技术中并没有能够实时辨识泵站特性曲线系数的方法。因此，需要一种研究泵站特性曲线系数的自动辨识方法，以提高含泵站河网水动力模型的精度和适应性。

发明内容

本发明的目的在于提供一种一维水动力模型泵站特性曲线系数的辨识方法，从而解决现有技术中存在的前述问题。

为了实现上述目的，本发明所述一维水动力模型泵站特性曲线系数的辨识方法，所述方法包括：

S1，基于圣维南方程，建立目标明渠的一维明渠水动力模型，并以此模型为基础，率定目标明渠中待辨识泵站的泵站前连接渠道的糙率和泵站后连接渠道的糙率；

S2，待辨识泵站边界以连续方程和泵站特性曲线为相容方程，与S1中的圣维南方程耦合，形成包含待辨识泵站边界的一维渠道-泵站水动力模型，求解所述一维渠道-泵站水动力模型得到所述目标明渠的水动力过程；

所述连续方程为公式(1)：

Q_j＝Q_j+1 (1)

其中，j表示待辨识泵站进水侧的断面编号，j+1表示待辨识泵站出水侧的断面编号；

所述泵站特性曲线为公式(2)：

其中，Q_j为泵站进水侧流量，Q_j+1为泵站出水侧流量，Z_j为进水侧水位，Z_j+1为出水侧水位，A，B，C均为泵站特性曲线系数；

S3，所述目标明渠的水动力过程的基础上，构建以待辨识泵站的三个泵站特性曲线系数为状态向量的状态空间方程，所述状态空间方程包括状态方程和量测方程；

所述状态方程为公式(3)：

Xⁿ⁺¹＝ΦⁿXⁿ+wⁿ⁺¹ (3)；

所述量测方程为公式(4)：

Yⁿ⁺¹＝HXⁿ⁺¹+vⁿ⁺¹ (4)；

公式(3)中，Xⁿ⁺¹为n+1时刻的状态向量，Xⁿ⁺¹＝[Aⁿ⁺¹Bⁿ⁺¹Cⁿ⁺¹]^T；Φⁿ为n时刻到n+1时刻的状态转移矩阵；Xⁿ为n时刻的状态向量，Xⁿ＝[AⁿBⁿCⁿ]^T，T表示行向量转置列向量的转置符号；wⁿ⁺¹为n+1时刻的扰动矩阵；

公式(4)中，Yⁿ⁺¹为n+1时刻的观测向量，即n+1时刻待辨识泵站的水力要素观测数据；H为一维渠道-泵站水动力模型，记为观测算子；vⁿ⁺¹为n+1时刻的观测噪声向量；

S4，在泵站机组角度或转速不变时段，获取目标明渠中待辨识泵站现场的水力要素观测数据并在该水力要素观测数据的基础上，采用集合卡尔曼滤波器同化状态空间方程，得到同化后的泵站特性曲线系数。

优选地，所述Φⁿ为单位矩阵I。

优选地，量测方程的观测矩阵H为非线性观测算子。

本发明所述方法以一维水动力模型中的泵站特性曲线系数为研究对象，构建以所述系数为状态向量的状态空间方程，采用集合卡尔曼滤波技术，获取同化系数，为提高河网水动力数值模拟的精度提供了科学依据和技术支撑。

本发明的有益效果是：

(1)现有技术一般采用模型试验获取泵站特性曲线系数，缺点为：试验周期长。而本发明所述方法利用当时水力要素的观测数据，调用同化模型，实现了在一维明渠水动力模拟过程中实时调整泵站特性调整曲线系数。

(2)现有技术采用模型试验方法获得的泵站特性曲线系数，存在比尺效应，直接应用于一维渠道-泵站水动力模型中，模拟精度低。而本发明采用观测数据反演得到泵站特性曲线系数，显著提高水动力模拟结果的精度。

(3)现有技术一般只能得到特定工况的泵站特性曲线系数，而本发明可以得到任意工况下的泵站特性曲线系数。

附图说明

图1为本发明方法流程图；

图2为渠道-泵站系统示意图；

图3为集合卡尔曼滤波同化泵站特性曲线系数流程图；

图4为埝头泵站进水侧水位模拟结果；

图5为埝头泵站进水侧水位模拟结果；

图6为埝头泵站流量模拟结果。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施方式仅仅用以解释本发明，并不用于限定本发明。

关于本申请所述方法更详细的解释说明：

图1为本发明方法流程图，本发明所述一维水动力模型泵站特性曲线系数的辨识方法，具体技术方案如下：

S1，基于圣维南方程(Saint Venant方程)，建立目标明渠的一维明渠水动力模型，并率定目标明渠中待辨识曲线系数泵站的泵站前连接渠道的糙率和泵站后连接渠道的糙率；一维明渠水动力模型用于渠道模拟，根据渠道模拟结果判断率定糙率是否合理。

不考虑旁侧入流情况下，明渠一维非恒定流的计算通常采用Saint Venant方程组描述，该方程包括连续方程和动量方程：

其中，x和t分别为空间坐标和时间坐标；A为渠道处的过流面积；Q为渠道流量；Z为水位；α为动量校正修正；g为重力加速度；C为谢才系数；R为水力半径。

采用Pressimann四点隐式差分格式离散Saint Venant方程，并用追赶法求解。该方法属于已有技术。

步骤S1中，率定目标明渠中待辨识曲线系数泵站的泵站前连接渠道的糙率和泵站后连接渠道的糙率时，以渠池上游边界水位的模拟值与实测值之差最小为目标。如图2为渠道-泵站系统示意图，假设渠道顺直，无其他局部建筑物。选取上游侧渠池和下游侧渠池，分别进行渠道糙率率定。率定的方法可以是人工试错法、优化算法等。各相邻渠道糙率率定完成后，可以认为渠道-泵站系统水动力模拟误差基本来自泵站边界的参数。

S2，以连续方程和泵站特性曲线为相容条件，与Saint Venant方程耦合求解。

在不考虑泵站内部流量损失的情况下，泵站入口流量和泵站出口流量基本相等，满足连续性条件。泵站内部多水泵机组运行工况变化时，水泵特性曲线会随之变化，导致泵站综合特性曲线也会随之变化。泵站特性曲线表征的是泵站的水位差和流量关系，属于泵站的固有属性。因此，以连续方程(式2)和泵站特性曲线(式3)作为相容条件，与SaintVenant方程耦合，求解含泵站河网的水动力过程。

Q_j＝Q_j+1 (2)

式中，Q_j为泵站进水侧流量，Q_j+1为泵站出水侧流量，Z_j为进水侧水位，Z_j+1为出水侧水位，A，B，C均为泵站特性曲线系数。若泵站中仅开一台机组，则系数与水泵机组系数相同；若泵站开多台机组，则曲线系数为多水泵的组合转角系数。运行工况变化时，组合系数随之发生变化。

步骤3，选取泵站进水侧和出水侧的水力要素观测数据作为观测量，构建以泵站特性曲线系数A，B，C为状态向量的状态空间方程，即状态方程和量测方程。

泵站特性曲线系数对水动力模拟结果的精度有重要的影响，特别是泵站作为内边界时，模拟误差会向上下游传递，影响更为显著。

水力要素观测数据为相同观测断面的水位或流量或不同观测断面的水位或流量。同时以泵站特性曲线系数A，B，C为状态向量，以X表示，构建状态空间方程。

状态方程见式(4)：

Xⁿ⁺¹＝ΦⁿXⁿ+wⁿ⁺¹ (4)

式中，Xⁿ⁺¹为n+1时刻的状态向量，Xⁿ⁺¹＝[Aⁿ⁺¹Bⁿ⁺¹Cⁿ⁺¹]^T；Φⁿ为n时刻到n+1时刻的状态转移矩阵；Xⁿ为n时刻的状态向量，Xⁿ＝[AⁿBⁿCⁿ]^T，T表示行向量转置列向量的转置符号；w^{n +1}为n+1时刻的扰动矩阵。

量测方程由一维渠道-泵站水动力模型转换得到，见式(5)：

Yⁿ⁺¹＝HXⁿ⁺¹+vⁿ⁺¹ (5)

式中，Yⁿ⁺¹为n+1时刻的观测向量，即n+1时刻泵站j的水力要素观测数据；H为一维渠道-泵站水动力模型,记为观测算子，H用于表示同一时刻的状态向量与观测向量的映射关系，在本申请所述方法中，n+1时刻的状态向量Xⁿ⁺¹为一维渠道-泵站水动力模型的输入参数；计算HXⁿ⁺¹得到目标明渠的水动力过程；vⁿ⁺¹为n+1时刻的观测噪声向量。

式(4)和式(5)中，假设wⁿ⁺¹和vⁿ⁺¹相互独立，并且分别服从均值为0，协方差矩阵为Qⁿ⁺¹和Rⁿ⁺¹的多元正态分布。

S4，选取泵站机组角度或转速不变时段观测数据，调用集合卡尔曼滤波器，同化曲线系数。

泵站特性曲线随水泵机组运行工况(如转角、转速等)的变化而改变，因此需要选取工况稳定时段的观测数据，逐次同化各工况下的曲线系数。

在一维水动力模型泵站曲线系数同化中，状态向量与观测向量的映射关系由一维明渠-泵站水动力模型H表示，H为非线性观测算子。考虑观测算子的非线性特点，采用集合卡尔曼滤波(EnKF)方法同化曲线系数。

假设待辨识的参数个数为M(M为3的倍数),集合个数为N。若状态向量在第n时刻的基本值为已知值，获得其在第n+1时刻扰动分析值的步骤如下：

(1)计算状态变量在第n+1时刻的扰动预报值：

当n＝0时，即初始时刻的由泵站模型试验系数加扰动生成。的总体作为基本集合。的总体称为预报集合。

(2)计算第n+1时刻的扰动新息(观测值减去预报值)：

(3)计算状态变量在第n+1时刻的扰动分析值：

其中，的整体称为分析集合。Pⁿ⁺¹H^T和HPⁿ⁺¹H^T分别定义为：

式(9)和式(10)中，X^n+1，f为扰动预报值的平均值，即：

(4)计算状态变量的最终同化值：

将n+1时刻的分析集合传递给该时刻的基本集合按时间顺序逐次执行滤波过程(1)-(3)，当n+1＝T(T为观测时段长度)时，取

X^T,a即为最终稳定的泵站同化曲线系数。

实施例

本实施例以密云水库调蓄工程埝头泵站为例，按照本申请描述方法，识别泵站曲线系数并根据泵站工况进行实时更新。

密云水库调蓄工程是南水北调中线北京段的配套工程，该工程将中线沿线省市的剩余水量加压输送至密云水库，能够增加密云水库蓄水量，回补密怀顺水源地，改善当地的生态环境。其中，团城湖到怀柔段，共包含6座泵站，以明渠相连接。泵站是系统的主要调控建筑物，准确识别泵站特性曲线系数，提高一维水动力数值模拟的精度，可以避免水资源的浪费，确保整个系统的安全运行。

研究时段为2016年4月26日5时-2016年4月27日19时。前柳林-埝头段：上游边界为前柳林泵站出水侧流量过程，下游边界为埝头泵站进水侧水位过程，率定目标为前柳林泵站出水侧水位过程模拟值与实测值之差最小，采用人工试错法进行率定，该段糙率取为0.015。埝头-兴寿段：上游边界为埝头泵站出水侧流量过程，下游边界为兴寿泵站进水侧水位过程，率定目标为埝头泵站出水侧水位过程模拟值与实测值之差最小，采用人工试错法进行率定，该段糙率取为0.031。

前柳林出水侧-兴寿进水侧组成的渠道-泵站系统中，仅含埝头泵站，在研究时段埝头泵站特性曲线初始系数为A＝-0.0337935，B＝-0.0419652，C＝5.42052。状态向量X即为埝头泵站曲线系数A，B，C组成的列向量。观测向量Y选择埝头泵站对应时刻进水侧水位、出水侧水位和流量的观测值。

埝头泵站在研究时段，机组工况稳定，没有进行转角或转速的调节。通过数据同化，得到该时段的同化曲线系数为：A＝-0.001699，B＝-0.031272，C＝2.169805。

分别采用初始系数和同化系数对渠道-泵站系统进行一维水动力数值模拟，模拟结果图4-6。结果表明，通过本发明方法进行泵站特性曲线参数辨识，能够显著提高一维水动力模型的精度。

通过采用本发明公开的上述技术方案，得到了如下有益的效果：

(1)现有技术一般采用模型试验获取泵站特性曲线系数，缺点为：试验周期长。而本发明所述方法利用当时水力要素的观测数据，调用同化模型，实现了在一维明渠水动力模拟过程中实时调整泵站特性调整曲线系数。

(2)现有技术采用模型试验方法获得的泵站特性曲线系数，存在比尺效应，直接应用于一维明渠-泵站水动力模型中，模拟精度低。而本发明采用观测数据反演得到泵站特性曲线系数，显著提高水动力模拟结果的精度。

(3)现有技术一般只能得到特定工况的泵站特性曲线系数，而本发明可以得到任意工况下的泵站特性曲线系数。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视本发明的保护范围。

Claims (3)

1.一种一维水动力模型泵站特性曲线系数的辨识方法，其特征在于，所述方法包括：

S1，基于圣维南方程，建立目标明渠的一维明渠水动力模型，并以此模型为基础，率定目标明渠中待辨识泵站的泵站前连接渠道的糙率和泵站后连接渠道的糙率；

S2，待辨识泵站边界以连续方程和泵站特性曲线为相容方程，与S1中的圣维南方程耦合，形成包含待辨识泵站边界的一维渠道-泵站水动力模型，求解所述一维渠道-泵站水动力模型得到所述目标明渠的水动力过程；

所述连续方程为公式(1)：

Q_j＝Q_j+1 (1)

其中，j表示待辨识泵站进水侧的断面编号，j+1表示待辨识泵站出水侧的断面编号；

所述泵站特性曲线为公式(2)：

<mrow> <msub> <mi>Z</mi> <mrow> <mi>j</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>-</mo> <msub> <mi>Z</mi> <mi>j</mi> </msub> <mo>=</mo> <msubsup> <mi>AQ</mi> <mrow> <mi>j</mi> <mo>+</mo> <mn>1</mn> </mrow> <mn>2</mn> </msubsup> <mo>+</mo> <msub> <mi>BQ</mi> <mrow> <mi>j</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>+</mo> <mi>C</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> <mo>;</mo> </mrow>

其中，Q_j为泵站进水侧流量，Q_j+1为泵站出水侧流量，Z_j为进水侧水位，Z_j+1为出水侧水位，A，B，C均为泵站特性曲线系数；

S3，所述目标明渠的水动力过程的基础上，构建以待辨识泵站的三个泵站特性曲线系数为状态向量的状态空间方程，所述状态空间方程包括状态方程和量测方程；

所述状态方程为公式(3)：

Xⁿ⁺¹＝ΦⁿXⁿ+wⁿ⁺¹ (3)；

所述量测方程为公式(4)：

Yⁿ⁺¹＝HXⁿ⁺¹+vⁿ⁺¹ (4)；

公式(3)中，Xⁿ⁺¹为n+1时刻的状态向量，Xⁿ⁺¹＝[Aⁿ⁺¹Bⁿ⁺¹Cⁿ⁺¹]^T；Φⁿ为n时刻到n+1时刻的状态转移矩阵；Xⁿ为n时刻的状态向量，Xⁿ＝[AⁿBⁿCⁿ]^T，T表示行向量转置列向量的转置符号；w^{n +1}为n+1时刻的扰动矩阵；

公式(4)中，Yⁿ⁺¹为n+1时刻的观测向量，即n+1时刻待辨识泵站的水力要素观测数据；H为一维渠道-泵站水动力模型，记为观测算子；vⁿ⁺¹为n+1时刻的观测噪声向量；

S4，在泵站机组角度或转速不变时段，获取目标明渠中待辨识泵站现场的水力要素观测数据并在该水力要素观测数据的基础上，采用集合卡尔曼滤波器同化状态空间方程，得到同化后的泵站特性曲线系数。

2.根据权利要求1所述水动力模型泵站特性曲线系数的辨识方法，其特征在于，所述Φⁿ为单位矩阵I。

3.根据权利要求1所述水动力模型泵站特性曲线系数的辨识方法，其特征在于，量测方程的观测矩阵H为非线性观测算子。