CN103886187B - 一种基于数据同化的河道水沙实时预测方法 - Google Patents

Abstract

本发明涉及一种基于数据同化的河道水沙实时预测方法，属于水利工程技术领域。首先收集待预测河道的地形数据、上下游的边界条件数据以及河道的断面数据，构建一维非恒定流和非平衡输沙模型，模型进行求解；接收实时观测数据时，对实时观测数据进行水沙模型的同化，并将同化值作为计算的初始场；根据初始场信息和边界条件数据计算未来的水位、流量以及含沙量的变化。本发明方法采用基于同化实时水位、流速、含沙量等数据，实现水位流速预测、洪水预报等功能，具有针对性强，功能齐全，方便实用等特点，可应用在大江大河的河道洪水实时预报中，为实际防汛指挥工作提供决策支持。

Description

一种基于数据同化的河道水沙实时预测方法

技术领域

本发明涉及一种基于数据同化的河道水沙实时预测方法，属于水利工程技术领域。

背景技术

在我国常见的十多种自然灾害中，洪灾发生频繁，其影响范围之广，造成损失之大，均居前列。洪水对人民生命财产、国民经济建设构成严重威胁，影响社会、经济的稳定和发展。因此，若能对河道水位、流量、含沙量等重要的行洪参数进行正确预测，可以有效的调整水库调度策略以及居民的疏散方案，最大程度的减少洪灾所造成的损失。

传统的预测系统主要采用历史数据回归的方法，如神经网络、支持向量机等对河道的水沙现状进行预测。这些方法忽略了河道水沙演进过程中的动力学特性和守恒规律，因此此类系统只适合于瞬时预测，无法辅助决策，达不到防洪减灾的目的。近年来，利用水沙数学模型对河道状态的演变进行数值计算成为了重要的预测手段，并在航道规划、整治和淤积预报等方面起到了关键的作用。然而，利用该模型进行实时水位、流量、含沙量等方面的预报时出现了误差大、精度低的问题，其主要原因是：(1)传统水沙方程存在的误差；(2)未知量不封闭或者条件参数不封闭；(3)模型参数不能适应边界条件的变化。

将原型观测资料纳入到水沙数学模型方程中，用于提高模型预测精度成为系统研究和开发的新方向。但是，直接将实时观测值代入水沙数学模型会导致方程计算失稳和整个计算区域的不和谐。数据同化技术综合考虑了模型结构误差、进出口边界条件误差和观测值误差，融合实时观测资料校正模型的预测值，提高模型预测精度，同时将目前的规划模型发展为实时校正模型，从而为水库调度提供决策支持、防洪减灾提供保障依据。

发明内容

本发明的目的是提出一种基于数据同化的河道水沙实时预测方法，采用基于同化实时的水位、流速、含沙量等数据，预测河道实时状态的系统，用于支持水库调度决策，并为保障防洪减灾提供依据。

本发明提出的基于数据同化的河道水沙实时预测方法，包括以下步骤：

(1)采集河道待预测段的信息，该信息包括：河道待预测段的河道底高程E和河宽B、河道待预测段上游边界的流量Q_up和含沙量S_up、河道待预测段下游边界的水位Z_down，各个断面的初始流量场Q_初始、初始水位场Z_初始以及初始含沙量场S_初始，将上述信息按采样时间的先后顺序分别进行排列，得到一张河道待预测段信息表；

(2)根据上述信息，求解以下一维非恒定流非平衡输沙模型：

$\frac{\∂A}{\∂t}+\frac{\∂Q}{\∂x}=q_l$

$\frac{\∂Q}{\∂t}+\frac{\∂}{\∂x}\left(\frac{Q^2}{A}\right)+gA\frac{\∂Z}{\∂x}+g\frac{Q\left|Q\right|}{C^{2}AR}=0$

$\frac{\∂\left(AS\right)}{\∂t}+\frac{\∂\left(QS\right)}{\∂x}+\mathrm{\α}\mathrm{\ω}B(S-S_{*})=0$

得到河道待预测段各断面的水位计算值Z_预报、流量计算值Q_预报和含沙量计算值S_预报，其中，A为河道待预测段各断面面积，g为重力加速度，R为河道待预测段各断面的水力半径，q_l为河道待预测段的旁侧流量，C为谢才系数，由曼宁公式计算求得n根据河道实际情况进行经验取值，α为恢复饱和系数，取值范围为0.01-10；ω为泥沙沉速，B为河道待预测段各断面的宽度，S_*为河道待预测段各断面的平均挟沙力，t为时间，x为河道待预测段中相邻两个预测断面之间的距离；

(3)实时采集待测河道各断面的信息，包括水位Z_测量、流量Q_测量和含沙量S_测量，根据实时采集待测河道各断面的信息，采用集合卡尔曼滤波方法，建立以下数据同化水沙模型：

Q_同化＝Q_预报+K(Q_测-HQ_预报)

Z_同化＝Z_预报+K(Z_测-HZ_预报)

S_同化＝S_预报+K(S_测-HS_预报)

其中，H为观测算子，一般情况下取1，K为卡尔曼增益，K＝P/(P+PS)，P为数据同化水沙模型计算误差，采用蒙特卡罗方法计算得出n次模拟的P_i计算得出，PS为观测值误差，根据《河流流量测验规范》，流量和含沙量的误差为测量值的5％，水位误差为3cm；

求解上述数据同化水沙模型，得到河道待预测段的上游边界的同化流量Q_同化、同化水位Z_同化及同化含沙量S_同化；

(4)根据设定的河道水沙预测周期，将上述步骤(3)得到的Q_同化、Z_同化、S_同化作为河道待预测段各断面的初始流量场Q_初始、初始水位场Z_初始以及初始含沙量场S_初始，集合步骤(1)的待测河道上游边界的流量Q_up、含沙量S_up、下游边界的水位Z_down，利用步骤(2)的一维非恒定流非平衡输沙模型，进行河道水沙预测计算，得到河道待预测段未来时刻各断面的未来水位Z_预报、未来流量Q_预报和未来含沙量S_预报。

本发明提出的基于数据同化的河道水沙实时预测方法，其优点是：

1、本发明的基于数据同化的河道水沙实时预测方法，针对实际工程需要，结合已有的实时水沙监测手段，构建水沙实时数据接收系统，实时获取水位、流量和含沙量等数据，将实测的水位、流量和含沙量等数据考虑进模型中，使得原有的水沙数学模型的应用范围从工程的设计和规划领域拓展到河道水情沙情的实时预报领域。

2、本发明采用先进的数据同化方法，同时考虑实测数据的测量误差和水沙模型的计算误差，对实测数据融入水沙模型的过程进行优化，使得模型预报的初始场取得最优值，从而有效地提高了水沙模型实时预报河道水情沙情的预报精度。

3、本发明方法系统地提出了实时水沙数据接收、水沙数学模型演算、数据同化初始场、水沙状态量预报等全部模块，完善了水沙模型实时预报系统框架，可实现水位流速预测、洪水预报等功能，具有针对性强，功能齐全，方便实用等特点，可应用在大江大河的河道洪水实时预报中，为实际防汛指挥工作提供决策支持。

附图说明

图1是本发明方法中涉及的河道待预测段的示意图，图1中，1是河道待测段的上游边界，2是河道待测段中的预测断面，3是河道待测段的下游边界。

图2是本发明方法中计算差分法涉及的四点偏心隐格式示意图。

具体实施方式

本发明提出的基于数据同化的河道水沙实时预测方法，包括以下步骤：

(1)采集河道待预测段的信息，该信息包括：河道待预测段的河道底高程E和河宽B、河道待测段的上游边界的流量Q_up和含沙量S_up、河道待测段的下游边界的水位Z_down，各个断面的初始流量场Q_初始、初始水位场Z_初始以及初始含沙量场S_初始，将上述信息按采样时间的先后顺序分别进行排列，得到一张河道待预测段信息表；图1所示是河道待预测段的示意图，其中，1是河道待测段的上游边界，2是河道待测段中的预测断面，3是河道待测段的下游边界。

(2)根据上述信息，求解以下一维非恒定流非平衡输沙模型：

$\frac{\∂A}{\∂t}+\frac{\∂Q}{\∂x}=q_l---\left(a\right)$

$\frac{\∂Q}{\∂t}+\frac{\∂}{\∂x}\left(\frac{Q^2}{A}\right)+gA\frac{\∂Z}{\∂x}+g\frac{Q\left|Q\right|}{C^{2}AR}=0---\left(b\right)$

$\frac{\∂\left(AS\right)}{\∂t}+\frac{\∂\left(QS\right)}{\∂x}+\mathrm{\α}\mathrm{\ω}B(S-S_{*})=0---\left(c\right)$

得到河道待预测段各断面的水位计算值Z_预报、流量计算值Q_预报和含沙量计算值S_预报，

其中，A为河道待预测段各断面面积(m²)，g为重力加速度(m/s²)，R为河道待预测段各断面的水力半径，q_l为河道待预测段的旁侧流量(m³)，C为谢才系数，由曼宁公式计算求得n根据河道实际情况进行经验取值，α为恢复饱和系数，取值范围为0.01-10；ω为泥沙沉速(m/s)，B为河道待预测段各断面的宽度(m)，S_*为河道待预测段各断面的平均挟沙力(kg/m³)，t为时间(s)，x为河道待预测段中相邻两个预测断面之间的距离(m)；

上述一维非恒定流非平衡输沙模型的求解方法可以有多种，本发明的一个实施例中，采用四点偏心隐格式法求解上述一维非恒定流非平衡输沙模型：

四点偏心隐式差分格式的特点是围绕矩形网格中的一点M来取偏导数并进行差商逼近。其网格的距离步长Δx可以是不等距的，而时间步长一般是等间距的，点M距离已知时层为θΔt，距未知时层为(1-θ)Δt，网格形式如图2所示，M处于距离步长中间。设每一矩形网格内函数f呈直线变化，则M点的函数值f可以由四个顶点的函数值表示如下：

根据四点偏心隐格式，上述一维非恒定流非平衡输沙模型中的方程(a)和(b)可差分为：

a_1jΔQ_j+b_1jΔZ_j+c_1jΔQ_j+1+d_1jΔZ_j+1＝E_1j

a_2jΔQ_j+b_2jΔZ_j+c_2jΔQ_j+1+d_2jΔZ_j+1＝E_2j

式中，

$a_{1j}=-\frac{4\mathrm{\θ}\mathrm{\Δ}t}{\mathrm{\Δ}x(B_j^n+B_{j+1}^n)}$

$b_{1j}=1-\frac{4\mathrm{\θ}\mathrm{\Δ}t(Q_{j+1}^n-Q_j^n)}{\mathrm{\Δ}x{(B_{j+1}^n+B_j^n)}^2}\·\frac{{\mathrm{dB}}_j^n}{{\mathrm{dz}}_j^n}$

$c_{1j}=\frac{4\mathrm{\θ}\mathrm{\Δ}t}{\mathrm{\Δ}x(B_j^n+B_{j+1}^n)}$

$d_{1j}=1-\frac{4\mathrm{\θ}\mathrm{\Δ}t(Q_{j+1}^n-Q_j^n)}{\mathrm{\Δ}x{(B_{j+1}^n+B_j^n)}^2}\·\frac{{\mathrm{dB}}_{j+1}^n}{{\mathrm{dz}}_{j+1}^n}$

$E_{1j}=-\frac{4\mathrm{\Δ}t}{\mathrm{\Δ}x}\frac{(Q_{j+1}^n-Q_j^n)}{(B_j^n+B_{j+1}^n)}$

$a_{2j}=1-\frac{4\mathrm{\θ}\mathrm{\Δ}t}{\mathrm{\Δ}x}\left(\frac{Q_j^n}{A_j^n}\right)+2g\mathrm{\θ}\mathrm{\Δ}t\frac{A_j^n\left|Q_j^n\right|}{{\left(K_j^n\right)}^2}$

$b_{2j}=\frac{\mathrm{\θ}\mathrm{\Δ}t}{\mathrm{\Δ}x}\[\frac{2{\left(Q_j^n\right)}^2{B}_j^n}{{\left(A_j^n\right)}^2}-g(A_{j+1}^n+A_j^n)+g(z_{j+1}^n-z_j^n)\·B_j^n\]+g\mathrm{\θ}\mathrm{\Δ}t\·\frac{Q_j^n\left|Q_j^n\right|}{{\left(K_j^n\right)}^2}\[B_j^n-\frac{2A_j^n}{K_j^n}\·\frac{{\mathrm{dK}}_j^n}{{\mathrm{dz}}_j^n}\]$

$c_{2j}=1+\frac{4\mathrm{\θ}\mathrm{\Δ}t}{\mathrm{\Δ}x}\left(\frac{Q_{j+1}^n}{A_{j+1}^n}\right)+2g\mathrm{\θ}\mathrm{\Δ}t\frac{A_{j+1}^n\left|Q_{j+1}^n\right|}{{\left(K_{j+1}^n\right)}^2}$

$d_{2j}=\frac{\mathrm{\θ}\mathrm{\Δ}t}{\mathrm{\Δ}x}\[-\frac{2{\left(Q_{j+1}^n\right)}^2{B}_{j+1}^n}{{\left(A_{j+1}^n\right)}^2}+g(A_{j+1}^n+A_j^n)+g(z_{j+1}^n-z_j^n)B_{j+1}^n\]+g\mathrm{\θ}\mathrm{\Δ}t\·\frac{Q_{j+1}^n\left|Q_{j+1}^n\right|}{{\left(K_{j+1}^n\right)}^2}\[B_{j+1}^n-\frac{2A_{j+1}^n}{K_{j+1}^n}\·\frac{{\mathrm{dK}}_{j+1}^n}{{\mathrm{dz}}_{j+1}^n}\]$

$E_{2j}=\frac{\mathrm{\Δ}t}{\mathrm{\Δ}x}\[-\frac{2{\left(Q_{j+1}^n\right)}^2}{\left(A_{j+1}^n\right)}+\frac{2{\left(Q_j^n\right)}^2}{A_j^n}-g(A_{j+1}^n+A_j^n)(z_{j+1}^n-z_j^n)\]-g\mathrm{\Δ}t\[\frac{A_{j+1}^n{Q}_{j+1}^n\left|Q_{j+1}^n\right|}{{\left(K_{j+1}^n\right)}^2}+\frac{A_j^n{Q}_j^n\left|Q_j^n\right|}{{\left(K_j^n\right)}^2}\]$

配以进口边界条件和下游边界条件方能构成求解封闭方程组，

a₁₁ΔQ₁+b₁₁Δz₁+c₁₁ΔQ₂+d₁₁Δz₂＝E₁₁ ①

a₁₂ΔQ₂+b₁₂Δz₂+c₁₂ΔQ₃+d₁₂Δz₃＝E₁₂ ②

a₂₁ΔQ₁+b₂₁Δz₁+c₂₁ΔQ₂+d₂₁Δz₂＝E₂₁ ③

a₂₂ΔQ₂+b₂₂Δz₂+c₂₂ΔQ₃+d₂₂Δz₃＝E₁₂ ④

ΔQ₁＝f₁(t) ⑤

Δz₃＝f₂(t) ⑥——边界条件

利用以上的6个方程就可以求解如图1所示的三个断面河段的六个未知数。

写成矩阵形式：

$\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ a_{11}& b_{11}& c_{11}& d_{11}& 0& 0\\ a_{21}& b_{21}& c_{21}& d_{21}& 0& 0\\ 0& 0& a_{12}& b_{12}& c_{12}& d_{12}\\ 0& 0& a_{22}& b_{22}& c_{22}& d_{22}\\ 0& 0& 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}\mathrm{\Δ}{Q}_1\\ {\mathrm{\Δz}}_1\\ {\mathrm{\ΔQ}}_2\\ {\mathrm{\Δz}}_2\\ {\mathrm{\ΔQ}}_3\\ {\mathrm{\Δz}}_3\end{array}\right]=\left[\begin{array}{c}{f}_1\left(t\right)\\ E_{11}\\ E_{21}\\ E_{12}\\ E_{22}\\ f_2\left(t\right)\end{array}\right]$

求解线性方程组即可求出计算时段内所有断面的流量Q_预报和水位Z_预报。

上述一维非恒定流非平衡输沙模型中的悬沙输移方程(c)差分格式如下：

$S_i^{j+1}=\frac{\mathrm{\Δ}t{\left({\mathrm{\αB\ωS}}_{*}\right)}_i^{j+1}+{\left(AS\right)}_i^j+\frac{\mathrm{\Δ}t}{{\mathrm{\Δx}}_{i-1}}\left|{\left(QS\right)}_{i-1}^{j+1}\right|}{A_i^{j+1}+\mathrm{\Δ}t{\left(\mathrm{\α}B\mathrm{\ω}\right)}_i^{j+1}+\frac{\mathrm{\Δ}t}{{\mathrm{\Δx}}_{i-1}}\left|Q_i^{j+1}\right|}$

将方程(a)(b)计算出的各断面流量Q_预报和水位Z_预报带入上述悬沙输移方程(c)的差分方程后，即可算出各断面的S_预报。

(3)实时采集待测河道各断面的信息，包括水位Z_测量、流量Q_测量和含沙量S_测量，根据实时采集待测河道各断面的信息，采用集合卡尔曼滤波方法，建立以下数据同化水沙模型：

Q_同化＝Q_预报+K(Q_测-HQ_预报)

Z_同化＝Z_预报+K(Z_测-HZ_预报)

S_同化＝S_预报+K(S_测-HS_预报)

其中，H为观测算子，一般情况下取1，K为卡尔曼增益，K＝P/(P+PS)，P为数据同化水沙模型计算误差，采用蒙特卡罗方法计算得出n次模拟的P_i计算得出，PS为观测值误差，根据《河流流量测验规范》，流量和含沙量的误差为测量值的5％，水位误差为3cm；

求解上述数据同化水沙模型，得到河道待预测段的上游边界的同化流量Q_同化、同化水位Z_同化及同化含沙量S_同化；

(4)根据设定的河道水沙预测周期，将上述步骤(3)得到的Q_同化、Z_同化、S_同化作为河道待预测段各断面的初始流量场Q_初始、初始水位场Z_初始以及初始含沙量场S_初始，集合步骤(1)的待测河道上游边界的流量Q_up、含沙量S_up、下游边界的水位Z_down，利用步骤(2)的一维非恒定流非平衡输沙模型，进行河道水沙预测计算，得到河道待预测段未来时刻各断面的未来水位Z_预报、未来流量Q_预报和未来含沙量S_预报。

Claims (1)

1.一种基于数据同化的河道水沙实时预测方法，其特征在于该方法包括以下步骤：

(1)采集河道待预测段的信息，该信息包括：河道待预测段的河道底高程E和河宽B、河道待预测段上游边界的流量Q_up和含沙量S_up、河道待预测段下游边界的水位Z_down，各个断面的初始流量场Q_初始、初始水位场Z_初始以及初始含沙量场S_初始，将上述信息按采样时间的先后顺序分别进行排列，得到一张河道待预测段信息表；

(2)根据上述信息，求解以下一维非恒定流非平衡输沙模型：

$\frac{\∂A}{\∂t}+\frac{\∂Q}{\∂x}=q_l$

$\frac{\∂Q}{\∂t}+\frac{\∂}{\∂x}\left(\frac{Q^2}{A}\right)+gA\frac{\∂Z}{\∂x}+g\frac{Q\left|Q\right|}{C^{2}AR}=0$

$\frac{\∂\left(AS\right)}{\∂t}+\frac{\∂\left(QS\right)}{\∂x}+\mathrm{\α}\mathrm{\ω}B(S-S_{*})=0$

得到河道待预测段各断面的水位计算值Z_预报、流量计算值Q_预报和含沙量计算值S_预报，

其中，A为河道待预测段各断面面积，g为重力加速度，R为河道待预测段各断面的水力半径，q_l为河道待预测段的旁侧流量，C为谢才系数，由曼宁公式计算求得n根据河道实际情况进行经验取值，α为恢复饱和系数，取值范围为0.01-10；ω为泥沙沉速，B为河道待预测段各断面的宽度，S_*为河道待预测段各断面的平均挟沙力，t为时间，x为河道待预测段中相邻两个预测断面之间的距离；

(3)实时采集待测河道各断面的信息，包括水位Z_测量、流量Q_测量和含沙量S_测量，根据实时采集待测河道各断面的信息，采用集合卡尔曼滤波方法，建立以下数据同化水沙模型：

Q_同化＝Q_预报+K(Q_测-HQ_预报)

Z_同化＝Z_预报+K(Z_测-HZ_预报)

S_同化＝S_预报+K(S_测-HS_预报)

其中，H为观测算子，H的取值为1，K为卡尔曼增益，K＝P/(P+PS)，P为数据同化水沙模型计算误差，采用蒙特卡罗方法计算得出n次模拟的P_i计算得出，PS为观测值误差，根据《河流流量测验规范》，流量和含沙量的误差为测量值的5％，水位误差为3cm；

求解上述数据同化水沙模型，得到河道待预测段的上游边界的同化流量Q_同化、同化水位Z_同化及同化含沙量S_同化；

(4)根据设定的河道水沙预测周期，将上述步骤(3)得到的Q_同化、Z_同化、S_同化作为河道待预测段各断面的初始流量场Q_初始、初始水位场Z_初始以及初始含沙量场S_初始，结合步骤(1)的待测河道上游边界的流量Q_up、含沙量S_up、下游边界的水位Z_down，利用步骤(2)的一维非恒定流非平衡输沙模型，进行河道水沙预测计算，得到河道待预测段未来时刻各断面的未来水位Z_预报、未来流量Q_预报和未来含沙量S_预报。