湖泊入湖流量的预报方法

技术领域

本发明属于湖泊水环境的技术领域，具体涉及一种湖泊入湖流量的预报方法。

背景技术

湖泊入湖流量实测数据是湖泊运行管理的基础性数据。例如，在编制湖泊水文预报方案时，湖泊实测入湖流量作为已知数据，是率定、验证水文模型参数和评价预报方案效率、精度的基准，但湖泊实测入湖流量估计中存在的误差给水文预报工作带来了极大难度；在湖泊调度中，湖泊实测入湖流量是最基本的输入条件，湖泊洪水演算、湖泊控制运用意见的编制以及湖泊运行经济评价等均以湖泊入湖流量数据为基础，因此精准的入湖流量数据也是正确开展湖泊控制运用的基础。

目前湖泊入湖流量的主要方法有：1)采用湖泊水量平衡方程和湖泊水位-湖容曲线反推湖泊入湖流量；2)采用流域产流系数，通过计算水位与实测水位比较修正产流系数，进而反推河流入湖流量。对于有水文站控制的入湖河道，依据水文站控制断面实测的逐日流量、对应关键时长的降雨总量、上游集水域面积，利用线性回归法，建立控制水文站实测流量与集水域对应关键时长内总降雨线性关系式，求得流域单位面积降雨产流系数；对于无径流资料河道，可假设在时段内其单位面积降雨产流系数与流域内有水文站控制的入湖河道相等；3)基于水文模型的入湖流量预报方法；4)基于圣维南方程的水动力学方法；5)采用智能算法进行湖泊入湖流量预报，如人工神经网络等。

其中水动力学方法计算复杂，而且需要详细的湖泊水文地质资料；人工神经网络方法虽然具有较高的精度，但结构不具有唯一性，收敛性较慢；水量平衡方程、湖泊水位-湖容曲线反推方法存在的主要问题是湖泊入湖流量必须由湖泊水量平衡方程反推，易呈现较大幅度的振荡，甚至出现负值；流域产流系数法反推河流入湖流量只适用于入湖河道有水文站控制的湖泊。

基于水文模型的入湖流量预报方法，物理概念明确，有多种水文模型可供选择，适用面广，主要难点是率定和验证水文模型参数需要3至5年及以上可靠的入湖流量实测数据。建模所需主要数据包括流域降雨、蒸发以及用于模型率定和验证的3至5年及以上地表径流数据，其中降雨量、蒸发量来源于研究区内的雨量观测站和蒸发量观测站。对于入湖河道没有水文站控制的湖泊，入湖径流无实测资料，不能直接采用实测入湖流量数据来率定模型参数。在湖泊入流研究中，丰富的水位观测资料，以及水位与湖泊容量的关系，使得采用水量平衡法计算湖泊入流成为可能。由于难以准确获得包括农业用水、工业用水、生活用水、渗透等在内的详细的日常用水数据，基于湖泊水量平衡方程的入湖流量反推结果不合理，甚至可能是负值。因此，需要开发出一种新的预估湖泊入湖流量的预报方法以解决没有水文站控制的湖泊，入湖径流无实测资料，不能直接采用实测入湖流量数据来率定模型参数的问题。

发明内容

本发明的目的在于针对现有技术的不足之处，提供一种湖泊入湖流量的预报方法，较好地解决了入湖河道没有水文站控制的湖泊，入湖径流缺乏流量监测数据，无法率定和验证水文模型参数的难题。

为解决上述技术问题，本发明采用如下技术方案：

一种湖泊入湖流量的预报方法，包括如下步骤：

步骤1：收集流域水文数据、气象数据及湖泊控制运用资料；

步骤2：基于步骤1中获得的数据和湖泊控制运用资料，建立水量平衡方程和湖泊水位-湖容曲线，并根据水量平衡方程和湖泊水位-湖容曲线预估湖泊入湖流量；

步骤3：对步骤2预估得到的湖泊入湖流量进行数据处理以校正不合理数值；

步骤4：将数据处理后的入湖流量预估值设定为实测值，同时输入步骤1中得到的流域降雨和蒸发数据，采用降雨径流模型模拟流域内的降雨径流过程，以实测值为目标，率定和验证降雨径流模型参数，从而获得降雨径流模型的参数初值，并以率定后的参数为降雨径流模型的输入参数，模拟获得入湖流量模拟值；

步骤5：根据步骤4模拟的湖泊入湖流量值，结合步骤2中建立的水量平衡方程和湖泊水位-湖容曲线，预测由降雨径流导致湖泊水位的逐日变化，获得降雨导致的湖泊水位预报值日变化序列；

步骤6：对湖泊水位预报值日变化序列和湖泊水位实测值日变化序列变化曲线进行比较；

步骤7：优化降雨径流模型参数，使得湖泊水位预报日变化序列变化曲线逼近湖泊水位实测日变化序列，获得降雨径流模型参数的优化值；

步骤8：以步骤7中优化得到的参数作为降雨径流模型的输入参数，采用降雨径流模型进行长时间系列湖泊入湖流量预报。

进一步地，步骤2中，建立的水量平衡方程和湖泊水位-湖容曲线方程组成的方程组：

S(t+1)＝S(t)+[I(t)-D(t)-E(t)-DF

Z

式中，S(t+1)、S(t)分别t时段末、初湖泊蓄水量；I(t)为t时段入湖流量，D(t)为t时段出湖流量，E(t)为t时段湖面蒸发水量折算的流量，DF

进一步地，步骤3中对预估得到的湖泊入湖流量进行数据处理得方法为：当基于湖泊水量平衡方程的入湖流量反推结出现负值，将该负值当作零处理。

进一步地，步骤7具体包括以下子步骤：

7.1确定因子和因子水平

因子取NAM模型的多个主要参数，水平取步骤4确定的主要参数的初值和

7.2 Taguchi方法选直交表

选择符合条件的直交表，将多个主要参数分别代入直交表中，筛选出多种方案；

7.3降雨径流模拟

输入逐日降雨数据、逐日蒸发数据，并将数据处理后的入湖流量预估值作为实测值，通过降雨径流模拟，完成直交表中多种方案的计算；

7.4输出湖泊入湖流量预报值

分别输出多种方案对应的入湖流量降雨径流模拟值；

7.5基于水量平衡方程和湖泊水位-湖容曲线，预测由降雨、径流导致湖泊水位的逐日变化，获得降雨导致的湖泊水位预报值日变化序列；

7.6计算衡量湖泊水位预报值日变化序列和实测值日变化序列拟合程度的指标；

7.7信噪比分析

Taguchi优化中，根据确定性系数R

其中，SN为信噪比；yi为某目标试验值；n为试验次数；s

将确定性系数R

7.8确定NAM参数最优参数

在信噪比S/N分析中，参数的最佳水平也就是信噪最大的水平；换言之，每个因子选择最大信噪比相对应的那个水平，基于这样的原则，通过分析各因子的信噪比响应可以组成因子间的优化目标最佳组合方案。

进一步地，主要参数包括地表蓄水层最大含水量U

与现有技术相比，本发明的有益效果为：

(1)提供了一种以流域降雨、蒸发作为输入，湖泊入湖流量作为输出，湖泊水位作为中间变量的湖泊入湖流量计算方法，较好的解决了入湖河道没有水文站控制的湖泊，入湖径流缺乏流量监测数据，无法率定和验证水文模型参数的难题；

(2)时间分辨率显著提高，目前大部分湖泊入湖水量仅能估算到逐月水量，无逐日流量数据；

(3)使得流域各类水文信息价值得到充分挖掘，使得获得的逐日流量数据、水位数据保持了高度的相容性和协调性。

附图说明

图1为本发明实施例中研究区域地理位置及雨量站分布图；

图2为本发明实施例的流程图；

图3为本发明实施例中2007-2010年入湖流量预估值分布图；

图4为本发明实施例中2015-2016年入湖流量预估值分布图；

图5为本发明实施例中数据处理后2007-2010年入湖流量预估值分布图；

图6为本发明实施例中数据处理后2015-2016年入湖流量预估值分布图；

图7为本发明实施例中率定期(2007-2010年)实测径流量与模拟径流对比图；

图8为本发明实施例中验证期(2015-2016)实测径流量与模拟径流对比图；

图9为本发明实施例中2007年湖泊水位预报值日变化序列1和湖泊水位实测值日变化序列2变化曲线的比较图；

图10为本发明实施例中2008年湖泊水位预报值日变化序列1和湖泊水位实测值日变化序列2变化曲线的比较图；

图11为本发明实施例中2009年湖泊水位预报值日变化序列1和湖泊水位实测值日变化序列2变化曲线的比较图；

图12为本发明实施例中2010年湖泊水位预报值日变化序列1和湖泊水位实测值日变化序列2变化曲线的比较图；

图13为本发明实施例中2016年湖泊水位预报值日变化序列1和湖泊水位实测值日变化序列2变化曲线的比较图；

图14为本发明实施例中Taguchi优化方法流程图；

图15为本发明实施例中Taguchi法主效应图；

图16为本发明实施例中1965-2016年湖泊入湖流量预报图。

具体实施方式

下面将结合本发明实施例对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。

下面结合具体实施例对本发明作进一步说明，但不作为本发明的限定。

本实施例以湖北省四湖流域进行说明，湖北省四湖流域地处江汉平原腹地，南滨长江、北临汉江及东荆河，总面积11547.5km

表1洪涝灾害年长湖最高水位表 单位：m

本实施以长湖湖堤以上的集水区为研究区域。长湖湖区内共有4个雨量站和2个水文站(图1)。降雨、蒸发均为逐日观测资料，数据来自湖北省水文水资源局。使用资料如下：①长湖湖区内烟墩集、万城、拾桥和老合台4个雨量站和河溶水文站2007年-2010年和2015年-2016年逐日降雨资料；②习家口水文站2007年-2010年和2015年-2016年逐日水位资料；③刘岭闸、习家口闸、荆襄河闸、双店闸2007年-2010年和2015年-2016年逐日出湖流量资料；④习家口水文站2007年-2010年和2015年-2016年逐日蒸发量资料；⑤来源于湖北省水利厅2016年长湖内垸的分洪数据。

参见图2，本发明提供一种湖泊入湖流量的预报方法，包括以下步骤：

步骤1：收集流域水文数据、气象数据及其湖泊控制运用资料；

选择2007-2010年和20152016年烟墩集、河溶、万城、拾桥、老合台五个站的逐日雨量的算术平均值作为降雨量时间序列文件的来源；选择习家口站逐日蒸发量代表湖区蒸发量时间序列文件的来源；在本实施例中，湖泊控制运用资料为当《湖北省大型排涝泵站调度及主要湖泊控制运用意见(鄂政发〔2011〕74号)》。

步骤2：基于步骤1的数据，建立水量平衡方程和湖泊水位-湖容曲线，并根据水量平衡方程和湖泊水位-湖容曲线预估湖泊入湖流量，其中，建立的水量平衡方程和湖泊水位-湖容曲线方程组成的方程组：

S(t+1)＝S(t)+[I(t)-D(t)-E(t)-DF

Z

根据该步骤的方法，本实施例2007-2010年计算结果见图3，2015-2016年见图4。

步骤3：对步骤2预估得到的湖泊入湖流量进行数据处理以校正不合理数值；由于难以准确获得包括农业用水、工业用水、生活用水、渗透等在内的详细的日常用水数据，基于湖泊水量平衡方程的入湖流量反推结果有时候不合理，出现负值，处理方法是将负值当作零处理。数据处理后2007-2010年计算结果见图5，2015-2016年见图6。

步骤4：将数据处理后的入湖流量预估值设定为实测值，同时输入步骤1中得到的流域降雨和蒸发数据，采用降雨径流模型(NAM)模拟流域内的降雨径流过程，以实测值为目标，率定和验证降雨径流模型参数，从而获得降雨径流模型(NAM)的参数初值，并以率定后的参数为降雨径流模型的输入参数，模拟获得入湖流量模拟值；

根据掌握的水文和气象数据情况，以2007-01-01至2010-12-31作为模型率定期，2015-01-01至2016-12-31作为模型验证期。参数率定由计算机自动率定程序完成，率定的参数见表2。

表2 NAM模型考虑的主要参数和本研究初始取值

率定期实测径流量与模拟径流对比见图7，率定期模型的确定性系数R

验证期实测径流量与模拟径流对比见图8，率定期模型的确定性系数R

步骤5：根据步骤4模拟的湖泊入湖流量值，结合步骤2中建立的水量平衡方程(公式1)和湖泊水位-湖容曲线(公式2)，预测由降雨径流导致湖泊水位的逐日变化，获得降雨导致的湖泊水位预报值日变化序列1；

步骤6：以确定性系数R

2007-2010年湖泊水位预报值日变化序列1和湖泊水位实测值日变化序列2变化曲线的比较见图9-图112；其中，率定期模型的确定性系数R

2016年湖泊水位预报值日变化序列1和湖泊水位实测值日变化序列2变化曲线的比较见图13；其中，率定期模型的确定性系数R

步骤7：通过Taguchi优化方法，优化NAM模型参数，使得湖泊水位预报日变化序列1变化曲线逼近湖泊水位实测日变化序列2，获得NAM模型参数的优化值，见图14；

(1)确定因子和因子水平

因子取NAM模型的10个主要参数(U

表3 NAM参数及其因子水平

(2)Taguchi方法选直交表

选择符合条件的直交表，将多个主要参数分别代入直交表中，筛选出多种方案；

在本实施例中，符合条件的直交表为L

表4 NAM参数直交表

(3)降雨径流(NAM)模拟

输入逐日降雨数据、逐日蒸发数据，并将数据处理后入湖流量预估值作为实测值，通过降雨径流模拟，完成直交表中多个方案(本实施例是27个方案)的计算，分别输出多个方案对应的入湖流量模拟值；

(4)分别输出多个方案对应的入湖流量降雨径流(NAM)模拟值，在本实施例中输出27个方案的降雨径流模拟值；

(5)基于水量平衡原理和湖泊水位-湖容曲线，预测由降雨径流导致湖泊水位的逐日变化，获得降雨导致的湖泊水位日变化序列1的预报值。

(6)计算衡量湖泊水位预报值和实测值拟合程度的指标：R

本实施例，使用以下两个指标评价模型的精度：

1)确定性系数R

2)水量平衡指数％WBI，评估洪水期湖泊的水平衡，其值为0表示完全满足水平衡条件；

式中：Z

本实施例计算得到的R

表5输出确定性系数R

(7)信噪比(S/N)分析

Taguchi优化中，根据确定性系数R

The-Larger-the-better(目标值越大越优，LTB)：

其中，SN为信噪比；y

将确定性系数R

表6 R

(8)确定NAM参数最优参数

在信噪比(S/N)分析中，参数的最佳水平也就是信噪比最大的水平。观察图15，根据确定性系数R

U

CK2＝52.60；TOF＝0.068；TIF＝0.305；TG＝0.301；CKBF＝4565。

步骤8：以步骤7优化后的参数作为降雨径流模型的输入参数，采用降雨径流模型进行长时间系列湖泊入湖流量预报；

采用Taguchi方法优化的NAM参数，开展1965-2016年湖泊入湖流量预报工作，结果见图16。

以上仅为本发明较佳的实施例，并非因此限制本发明的实施方式及保护范围，对于本领域技术人员而言，应当能够意识到凡运用本发明说明书内容所作出的等同替换和显而易见的变化所得到的方案，均应当包含在本发明的保护范围内。