CN106168991A - 一种基于水动力数值模拟的感潮河网潮位预报方法 - Google Patents

Abstract

本发明提供了一种感潮河网内潮位预报的方法，综合运用了计算流体力学(CFD)方法和水位预报方法，解决了感潮河网内水位受多汊道分流及外海潮汐影响而难以准确预报的问题，提高感潮河网潮位预报的适用范围和精度，属于水动力模拟和水文预报技术领域。主要包括以下步骤，1)建立感潮河网区及上游河道一维水动力数学模型；2)利用调和分析预测河口控制站逐时潮位过程并利用预报前一日实测数据进行预报值校正；3)然后采用加权移动平均方法推求上游控制站点流量过程；4)利用河网区及上游干流区一维水动力数学模型求解感潮河网区各断面潮位过程，实现河网区任意断面潮位预报。

Description

一种基于水动力数值模拟的感潮河网潮位预报方法

技术领域

本发明涉及水动力模拟和水文预报技术领域，特别涉及一种基于水动力数值模拟的感潮河网潮位预报方法。

背景技术

随着河口及三角洲经济建设的不断发展，对水文预报工作提出了更高要求，如何快速准确地进行潮位预报，成为当前急需解决的问题。防汛决策、水利枢纽设计参考、闸泵群启闭方案优化等对感潮网河区的潮位预报需求日益增加。感潮河网的水位变化过程受上游径流、下游潮汐、河汊分流、区间入流及降水等多种因素共同影响，其水文特性既有内陆河流的特性，表现出洪、中、枯水期的季节变化，又有海洋的某些特征，表现出大、中、小潮的变化和潮汐引起水位涨落与流速流向的变化，此外，还受陆地气候的影响，上游径流不仅年内有明显的变化，而且在年际间变化也很不相同；受河道地形的影响也较大，情况十分复杂，感潮河网内潮位预报是长期以来河口及三角洲治理开发中的一个难点问题。

传统的潮位预报多为高水位预报，一般采用相应水位法，即由某站水位、传输距离推求另一站水位，其预报有效期受传输时间左右，无法满足长期、中期预报需求；另有分离过程法将预报站水位分离成上、下游水位相关的两部分分别进行预报，其不足之处与相关法相同，主要是预报时间较为受限；神经网络法预报水位则主要存在可能由于局部极值难收敛、振幅不大时较精准、高水位预报误差大的问题。为满足三角洲对潮位预报需求，发展一种基于数值模拟的能实现感潮河网内任意水域潮位快速预报的方法是十分必要的。

发明内容

本发明的目的在于克服河网地区现有潮位预报方法误差较大、适用范围小的不足，提供一种基于水动力数学模型的感潮河网任意水域潮位快速预报及实时校正方法，提高感潮河网潮位预报的适用范围和精度。

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

一种基于水动力数值模拟的感潮河网潮位预报方法，包括以下步骤：

步骤S1、根据预报位置选择合适的研究范围，确定河口下边界及上游控制站点与边界；

步骤S2、对研究范围内河道进行断面剖分和地形插值；

步骤S3、建立上游控制站至下游河口边界之间的网河区的一维非恒定流水动力数学模型；

步骤S4、采用实测水文数据对一维非恒定流水动力数学模型进行率定和验证；

步骤S5、根据河口控制站点调和参数进行潮位过程预报计算：

步骤S6、根据实测数据对河口潮位边界进行校正；

步骤S7、采用加权移动平均方法推求上游控制站点流量过程；

步骤S8、根据实时实测数据对上游流量边界进行校正；

步骤S9、输入预报上游边界流量和下游河口边界潮位过程，进行网河区的一维非恒定流水动力数学模型计算，求得网河区内各断面预潮位过程。

进一步地，所述步骤S3、建立上游控制站至下游河口边界之间的网河区的一维非恒定流水动力数学模型具体过程如下：

步骤S31、采用圣维南方程组作为一维潮流数学模型，方程如下：

连续方程：

动量方程：

将相关支流汇入或流出点作为网河区汊点，汊口点水流按以下连接条件控制：

流量连接条件：

水位连接条件：Z_i,j＝Z_m,n＝……＝Z_l,k；

式中，Z为断面平均水位；Q，A，R分别为断面流量、过水面积、水力 半径；B为水面宽度；X,T为距离和时间；G为重力加速度；谢才系数C＝R^1/6/n，n为糙率；

步骤S32、采用四点加权Preissmann固定网格隐式差分格式将方程离散，采用一维河网三级联解算法求解方程，具体求解方法为：进行网格划分，以S代表流量Q和水位Z，则S在Δx河段、Δt时段内的加权平均量及相应偏导数可分别表示为：

式中：θ为加权系数，θ一般取0.5～1.0。

进一步地，所述步骤S5、根据河口控制站点调和参数进行潮位过程预报计算具体过程如下：

步骤S51、有分潮参数的直接根据预报时间计算潮位过程；

步骤S52、无分潮参数，但有长于一年的实测历史数据的，调和分析各分潮参数，然后重复步骤S51；

步骤S53、分潮参数和历史数据均无则参考临近站点简化处理。

进一步地，所述步骤S51、有分潮参数的直接根据预报时间计算潮位过程的具体计算公式如下：

式中：为t时刻的潮位预报值，m；A₀为从某基准面算起的平均海平面高度，m；f_i为第i个分潮的交点因子；H_i为第i个分潮的振幅，m；σ_i为分潮角速度；(v0+u)_i为分潮的天文初相角；g_i为分潮迟角；i为分潮编号；N为分潮总数；H_i、g_i亦称为调和常数。

进一步地，所述步骤S6、根据实测数据对河口潮位边界进行修正具体过程如下：

步骤S61、在[t0-26,t0]范围内搜索判别出高高潮位Z_H ^PRE、低低潮位Z_L ^PRE，计算潮幅Z_AM ^PRE，Z_AM ^PRE＝Z_H ^PRE+Z_L ^PRE；

步骤S62、同样在[t0-26,t0]范围内搜索出实测潮位过程中的高高潮位Z_H ^MESU、低低潮位Z_L ^MESU，计算潮幅Z_AM ^MESU，Z_AM ^MESU＝Z_H ^MESU+Z_L ^MESU；

步骤S63、根据步骤S61及步骤S62结果对预报潮位中轴线进行调整，使得Z_AM ^PRE＝Z_AM ^MESU；

步骤S64、计算中轴线偏移量DIF_Z_ZERO，DIF_Z_ZERO＝(Z_AM ^PRE-Z_AM ^MESU)/2.0；

步骤S65、对于预报潮位过程值Z^PRE逐一进行中轴线偏移校正，

Z^PRE＝Z^PRE-DIF_Z_ZERO；

步骤S66、计算潮幅偏差DIF_Z_AM，DIF_Z_AM＝Z_AM ^PRE/Z_AM ^MESU；

进一步地，所述步骤S7中加权移动平均方法的公式如下：

$Q_{pree}^t=Q_A^t*2+Q_{pree}^{t-1}$

式中：为t时刻预报流量，m³/s；为第t时刻的流量趋势率的修正量，m³/s；为第t-1时刻的预报流量，m³/s。

进一步地，所述流量趋势率的修正量的求解步骤如下：

步骤S71、对于每一段预报期(t0～tn)内的首、末趋势率： 式中：分别为预报期内首时刻t0、末时刻tn对应流量趋势率；Q_T0-1和Q_T0-2分别为预报前1、前2小时的实测流量；

步骤S72、按线性插值法计算每一段预报期内首末时刻以外其它各时刻趋势率，

式中：为预报期内首末时刻以外其它各时刻趋势率；T为预报期总时长。

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

本发明公开了一种基于水动力数学模型的感潮河网任意水域潮位快速预报及实时校正方法，综合运用了计算流体力学(CFD)方法和水位预报方法，解决了感潮河网内水位受多汊道分流及外海潮汐影响而难以准确预报的问题，提高感潮河网潮位预报的适用范围和精度。

附图说明

图1是本发明公开的的一种基于水动力数值模拟的感潮河网潮位预报方法 的流程步骤图；

图2是一维河网水动力模型断面剖分示意图；

图3(a)是本发明在中顺大围调度系统中一维河网水动力模型研究范围图；

图3(b)是本发明在中顺大围调度系统中中顺大围区域的研究范围图；

图4(a)是98.6洪水期三沙口站水位过程验证成果图；

图4(b)是98.6洪水期勒竹站水位过程验证成果图；

图4(c)是05.1枯水期小榄站水位过程验证成果图；

图4(d)是05.1枯水期三沙口站水位过程验证成果图；

图4(e)是05.1枯水期竹银站水位过程验证成果图；

图4(f)是05.1枯水期亭角站水位过程验证成果图；

图5为河口控制点实测潮位、预报潮位、预报修正潮位对比图；

图6为上游梧州站实测流量、预报流量、预报修正流量对比图；

图7(a)东河水闸实测潮位与预报潮位过程对比图；

图7(b)西河水闸实测潮位与预报潮位过程对比图。

具体实施方式

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

实施例

请参见图1，图1是本实施例中公开的一种基于水动力数值模拟的感潮河网潮位预报方法的流程步骤图。图1所示的一种基于水动力数值模拟的感潮河网潮位预报方法，具体包括以下步骤：

步骤S1、根据预报位置选择合适的研究范围，确定河口下边界及上游控制站点与边界；

步骤S2、对研究范围内河道进行断面剖分和地形插值；

步骤S3、建立上游控制站至下游河口边界之间的网河区的一维非恒定流水动力数学模型；具体步骤如下：

步骤S31、一维潮流数学模型采用圣维南方程组，方程如下：

连续方程：

动量方程：

相关支流汇入或流出点作为网河区汊点，汊口点水流按以下连接条件控制：

流量连接条件：

水位连接条件：Z_i,j＝Z_m,n＝……＝Z_l,k；

式中，Z为断面平均水位；Q，A，R分别为断面流量、过水面积、水力半径；B为水面宽度；X,T为距离和时间；G为重力加速度；谢才系数C＝R^1/6/n，n为糙率。

步骤S32、方程离散采用四点加权Preissmann固定网格隐式差分格式，采用目前应用广泛的一维河网三级联解算法，网格布置图2所示。

具体求解方法为：对于图2所示的网格划分，以S代表流量Q和水位Z，则S在Δx河段、Δt时段内的加权平均量及相应偏导数可分别表示为：

式中：θ为加权系数，θ一般取0.5～1.0。

步骤S4、采用实测水文数据对一维非恒定流水动力数学模型进行率定和验证；

步骤S5、根据河口控制站点调和参数进行潮位过程预报计算，具体过程如下：

步骤S51、有分潮参数的直接根据预报时间计算潮位过程；

步骤S52、无分潮参数，但有长于一年的实测历史数据的，调和分析各分潮参数，然后重复步骤S51；

步骤S53、分潮参数和历史数据均无则参考临近站点简化处理。

式中：为t时刻的潮位预报值，m；A₀为从某基准面算起的平均海平 面高度，m；f_i为第i个分潮的交点因子；H_i为第i个分潮的振幅，m；σ_i为分潮角速度；(v0+u)_i为分潮的天文初相角；g_i为分潮迟角；i为分潮编号；N为分潮总数；H_i、g_i亦称为调和常数。

步骤S6、根据实测数据对河口潮位边界进行修正；

步骤S7、采用加权移动平均方法推求上游控制站点流量过程，并同时进行修正；

式中：为t时刻预报流量，m³/s；为第t时刻的流量趋势率(修正量)，m³/s；为第t-1时刻的预报流量，m³/s。

步骤S8、输入预报上游边界流量和下游河口边界潮位过程，进行网河区的一维非恒定流水动力数学模型计算，求得网河区内各断面预潮位过程。

其中，步骤S6所述潮位修正步骤如下：

步骤S61、在[t0-26,t0](预报起点前一个潮期，约26h)范围内搜索判别出高高潮位Z_H ^PRE、低低潮位Z_L ^PRE，计算潮幅Z_AM ^PRE；

步骤S62、同样在[t0-26,t0]范围内搜索出实测潮位过程中的高高潮位Z_H ^MESU、低低潮位Z_L ^MESU，计算潮幅Z_AM ^MESU；

步骤S63、根据步骤S61及步骤S62结果对预报潮位中轴线进行调整，使得(Z_H ^PRE+Z_L ^PRE)＝(Z_H ^MESU+Z_L ^MESU)；

步骤S64、计算中轴线偏移量DIF_Z_ZERO，

DIF_Z_ZERO＝((Z_H ^PRE+Z_L ^PRE)-(Z_H ^MESU+Z_L ^MESU))/2.0；

步骤S65、对于预报潮位过程值Z^PRE逐一进行中轴线偏移校正，

Z^PRE＝Z^PRE-DIF_Z_ZERO；

步骤S66、计算潮幅偏差DIF_Z_AM，DIF_Z_AM＝Z_AM ^PRE/Z_AM ^MESU；

其中，步骤S7所述流量修正中的求解步骤如下：

步骤S71、对于每一段预报期(t0～tn)内的首、末趋势率： 式中：分别为预报期内首时刻t0、末时刻tn对应流量趋势率；Q_T0-1和Q_T0-2分别为预报前1、前2小时的实测流量。

步骤S72、按线性插值法计算每一段预报期内首末时刻以外其它各时刻趋势率，

式中：为预报期内首末时刻以外其它各时刻趋势率；T为预报期总时长。

实施例二

以本发明在珠江三角洲河网中顺大围调度系统外江潮位预报实例及附图对本发明作进一步详细的描述。

(1)根据中顺大围的位置特性，确定潮位预报模型研究范围——包了西、北江三角洲、东江三角洲及广州出海水道、潭江水道；上边界取以下控制水文(位)站：西江梧州、北江石角、老鸦岗(流溪河)、麒麟咀(增江)、东江博罗、潭江石咀站；下边界取至八大口门控制站点：大虎、南沙、冯马庙、横门和磨刀门灯笼山、黄金、西炮台、官冲，见附图3(a)和图3(b)；

(2)对河道进行断面剖分和地形插值，共设4140个断面，采用2005年实测地形；

(3)建立包括珠江三角洲网河区的一维非恒定流水动力数学模型；

(4)采用1998年6月、1999年7月、2005年1月实测水文资料对珠江三角洲网河区一维非恒定流水动力数学模型进行率定和验证，部分成果见图4(a)～图4(f)；

(5)根据河口控制站点调和参数进行潮位过程预报计算，8个下边界控制站点无分潮参数资料，其分潮参数均由1～2年实测历史数据调和分析而得；

(6)根据实测数据对河口潮位边界进行校正，以灯笼山为例，实测值、预测值、预测校正值如图5所示，所有控制站点潮位预报误差统计见表1、表2所示；

(7)采用加权移动平均方法推求上游控制站梧州站流量过程，流量预报采用15小时滚动预报；

(8)根据实测数据对上游梧州站流量边界进行预报及校正，实测流量、预报流量及经校正后的预报流量过程如图6所示。

表1控制站点潮位预报误差统计(2014.05.01.0:00～5.31.23:00一个月)

表2控制站点2～4个月预报期潮位预报误差统计

(9)输入预报的上游梧州站流量过程和河口八大口门潮位过程，进行一维非恒定流水动力数学模型计算，求得中顺大围内石岐河与外江相接的两个主要控制闸——东河水闸、西河水闸预报潮位过程，如图7(a)、图7(b)所示。

综上所述，本发明公开的一种基于水动力数值模拟的感潮河网潮位预报方法，首先，建立感潮河网区及上游河道一维水动力数学模型；其次，利用调和分析预测河口控制站潮位过程并利用预报前一日实测数据进行预报值校正；然后采用加权移动平均方法推求上游控制站点流量过程；最后利用河网区及上游干流区一维水动力数学模型求解感潮河网区各断面潮位过程，实现河网区任意断面潮位预报。

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

Claims (7)

1.一种基于水动力数值模拟的感潮河网潮位预报方法，其特征在于，包括以下步骤：

步骤S1、根据预报位置选择合适的研究范围，确定河口下边界及上游控制站点与边界；

步骤S2、对研究范围内河道进行断面剖分和地形插值；

步骤S3、建立上游控制站至下游河口边界之间的网河区的一维非恒定流水动力数学模型；

步骤S4、采用实测水文数据对一维非恒定流水动力数学模型进行率定和验证；

步骤S5、根据河口控制站点调和参数进行潮位过程预报计算：

步骤S6、根据实测数据对河口潮位边界进行校正；

步骤S7、采用加权移动平均方法推求上游控制站点流量过程；

步骤S8、根据实时实测数据对上游流量边界进行校正；

步骤S9、输入预报上游边界流量和下游河口边界潮位过程，进行网河区的一维非恒定流水动力数学模型计算，求得网河区内各断面预潮位过程。

2.根据权利要求1所述的一种基于水动力数值模拟的感潮河网潮位预报方法，其特征在于，

所述步骤S3、建立上游控制站至下游河口边界之间的网河区的一维非恒定流水动力数学模型具体过程如下：

步骤S31、采用圣维南方程组作为一维潮流数学模型，方程如下：

连续方程：

动量方程：

将相关支流汇入或流出点作为网河区汊点，汊口点水流按以下连接条件控制：

流量连接条件：

水位连接条件：Z_i,j＝Z_m,n＝……＝Z_l,k；

式中，Z为断面平均水位；Q，A，R分别为断面流量、过水面积、水力半径；B为水面宽度；X,T为距离和时间；G为重力加速度；谢才系数C＝R^1/6/n，n为糙率；

步骤S32、采用四点加权Preissmann固定网格隐式差分格式将方程离散，采用一维河网三级联解算法求解方程，具体求解方法为：进行网格划分，以S代表流量Q和水位Z，则S在Δx河段、Δt时段内的加权平均量及相应偏导数可分别表示为：

$\left\{\begin{array}{c}\frac{\∂S}{\∂t}=\frac{S_{i+1}^{n+1}+S_i^{n+1}-S_{i+1}^n-S_i^n}{2\mathrm{\Δ}t}\\ \frac{\∂S}{\∂x}=\mathrm{\θ}\frac{S_{i+1}^{n+1}+S_i^{n+1}}{{\mathrm{\Δx}}_i}+(1-\mathrm{\θ})\frac{S_{i+1}^n-S_i^n}{{\mathrm{\Δx}}_i}\\ S=\frac{1}{4}(S_{i+1}^{n+1}+S_i^{n+1}+S_{i+1}^n+S_i^n)\end{array}\right.$

式中：θ为加权系数，θ取值区间为0.5～1.0。

3.根据权利要求1所述的一种基于水动力数值模拟的感潮河网潮位预报方法，其特征在于，

所述步骤S5、根据河口控制站点调和参数进行潮位过程预报计算具体过程如下：

步骤S51、有分潮参数的直接根据预报时间计算潮位过程；

步骤S52、无分潮参数，但有长于一年的实测历史数据的，调和分析各分潮参数，然后重复步骤S51；

步骤S53、分潮参数和历史数据均无则参考临近站点简化处理。

4.根据权利要求3所述的一种基于水动力数值模拟的感潮河网潮位预报方法，其特征在于，

所述步骤S51、有分潮参数的直接根据预报时间计算潮位过程的具体计算公式如下：

$z_{pree}^t=A_0+\underset{i=1}{\overset{N}{\Σ}}{f}_i{H}_{i}cos\[{\mathrm{\σ}}_i+{(v0+u)}_i-g_i\]$

式中：为t时刻的潮位预报值，m；A₀为从某基准面算起的平均海平面高度，m；f_i为第i个分潮的交点因子；H_i为第i个分潮的振幅，m；σ_i为分潮角速度；(v0+u)_i为分潮的天文初相角；g_i为分潮迟角；i为分潮编号；N为分潮总数；H_i、g_i亦称为调和常数。

5.根据权利要求1所述的一种基于水动力数值模拟的感潮河网潮位预报方法，其特征在于，

所述步骤S6、根据实测数据对河口潮位边界进行修正具体过程如下：

步骤S61、在[t0-26,t0]范围内搜索判别出高高潮位Z_H ^PRE、低低潮位Z_L ^PRE，计算潮幅Z_AM ^PRE，Z_AM ^PRE＝(Z_H ^PRE+Z_L ^PRE)；

步骤S62、同样在[t0-26,t0]范围内搜索出实测潮位过程中的高高潮位Z_H ^MESU、低低潮位Z_L ^MESU，计算潮幅Z_AM ^MESU，Z_AM ^MESU＝(Z_H ^MESU+Z_L ^MESU)；

步骤S63、根据步骤S61及步骤S62结果对预报潮位中轴线进行调整，使得Z_AM ^PRE＝Z_AM ^MESU；

步骤S64、计算中轴线偏移量DIF_Z_ZERO，DIF_Z_ZERO＝(Z_AM ^PRE-Z_AM ^MESU)/2.0；

步骤S65、对于预报潮位过程值Z^PRE逐一进行中轴线偏移校正，

Z^PRE＝Z^PRE-DIF_Z_ZERO；

步骤S66、计算潮幅偏差DIF_Z_AM，DIF_Z_AM＝Z_AM ^PRE/Z_AM ^MESU。

6.根据权利要求1所述的一种基于水动力数值模拟的感潮河网潮位预报方法，其特征在于，

所述步骤S7中加权移动平均方法的公式如下：

$Q_{pree}^t=Q_A^t*2+Q_{pree}^{t-1}$

式中：为t时刻预报流量，m³/s；为第t时刻的流量趋势率的修正量，m³/s；为第t-1时刻的预报流量，m³/s。

7.根据权利要求6所述的一种基于水动力数值模拟的感潮河网潮位预报方法，其特征在于，

所述流量趋势率的修正量的求解步骤如下：

步骤S71、对于每一段预报期(t0～tn)内的首、末趋势率： 式中：分别为预报期内首时刻t0、末时刻tn对应流量趋势率；Q_T0-1和Q_T0-2分别为预报前1、前2小时的实测流量；

步骤S72、按线性插值法计算每一段预报期内首末时刻以外其它各时刻趋势率，

式中：为预报期内首末时刻以外其它各时刻趋势率；T为预报期总时长。