基于经验小波消噪和神经网络融合的中长期径流预报方法与流程

本发明涉及水文学中的水文预报领域，特别涉及一种基于经验小波消噪和神经网络融合的中长期径流预报方法。

背景技术：

中长期径流预报(遇见期3天～1年)是近几十年来水文研究的热点问题之一，提高中长期径流预报精度对水资源调配、防洪、减灾等方面具有重要意义。现有中长期径流预报方法可以分为两类：物理模型和经验模型。物理模型一般通过模拟流域水循环的物理过程对未来径流进行预报。物理模型的构建需要获取较长时间段(遇见期)内的降水、蒸发和气温等气象数据，增加了采用物理模型进行中长期预报的难度。

中长期水文预报中涉及的各种经验模型主要包括时间序列模型和机器学习模型。然而，以往的研究表明，时间序列模型对非线性和非平稳性较强的径流时间序列预报能力有限。在各种机器学习方法中，人工神经网络(ann)包括反向传播神经网络(bpnn)、径向基函数(rbf)神经网络、广义回归神经网络(grnn)、elman神经网络和极限学习机(elm)等在中长期径流预报领域取得了丰硕的成果，但是单一ann模型的预报性能仍存在很大的改进空间。其中利用数据预处理技术对径流时间序列进行预处理是一种能有效提升中长期径流预报精度的方法。此外，另一种能够提升单一ann模性预报性能的方法是将多个ann模型的预报信息进行融合，综合考虑每个模型的优势，进而提高中长期径流预报精度。

技术实现要素：

发明目的：针对现有技术中存在的问题，本发明提供一种基于经验小波消噪和神经网络融合的中长期径流预报方法，能够有效提升中长期径流时间序列的预报精度。

技术方案：本发明提供了一种基于经验小波消噪和神经网络融合的中长期径流预报方法，其特征在于，包括以下步骤：步骤一：获取水文站点的历史实测径流数据，建立中长期径流时间序列，并将样本数据分为训练期和检验期；步骤二：采用ewt将径流时间序列分解为几个独立的经验模态，去除分解经验模态中频率最高的模态，并对剩余的经验模态进行线性求和，得到重构径流时间序列以消除原始径流时间序列的冗余噪声；步骤三：对经ewt处理后的径流时间序列进行相空间重构，构造相空间矩阵作为基预报模型的输入因子，建立ann基预报模型的输入、输出矩阵；步骤四：采用步骤三确定的训练期输入、输出矩阵分别对rbf、elm和elman三种神经网络模型进行训练；步骤五：将三种基预报模型的输出作为grnn模型的预报因子，再次建立输入、输出矩阵，对grnn模型进行训练，并将检验样本代入训练好的grnn模型，得到检验期的预测值；步骤六：采用四种评价指标包括均方根误差(rmse)、平均绝对误差(mae)、平均绝对百分误差(mape)和相关系数(r)评价对预报结果进行评价。

进一步地，在所述步骤二中，对于原始径流时间序列x(t)，ewt的分解过程如下：(1)根据快速傅立叶变换算法(fft)计算原始径流时间序列x(t)的傅立叶频谱f(ω)；(2)将傅立叶频谱f(ω)的频域[0,π]自适应的分割为k个带宽不等的频带，[0,ω1],[ω1,ω2],...,[ωk-1,π]，其中，小波滤波器的边界ωn,n＝1,2,...,n-1选取为两个连续局部极大值之间的中间频率；(3)根据meyer小波构造经验小波，确定经验小波函数ψk(ω)和经验尺度函数

式中，β(x)＝x⁴(35-84x+70x²-20x³)；(4)重构原始径流时间序列得到不同的经验模态。

优选地，在所述(4)中，原始径流时间序列通过下式进行重构：

式中，*为卷积运算，wx(0,t)为近似系数，wx(k,t)为细节系数，近似系数和细节系数分别通过下式给出：

式中，ψk(t)和ψk(t)分别表示经验小波函数和经验尺度函数，和分别表示ψk(τ-t)和的复共轭，f^-1(·)表示傅里叶逆变换，和分别表示ψ(ω)和的傅里叶变换；则经验模态uk(t)的定义如下：

优选地，rbf神经网络的隐含层使用高斯传递函数作为激活函数，高斯激活函数定义为：式中，x＝[x1,x2,...,xn]^t为n维的输入向量，ci＝[ci1,ci2,...,cin]^t为第i个隐含层神经元的中心，qi表示高斯函数的宽度，·为欧氏范数；则输出层第j个结点的响应为：式中，wij为第i个隐含层结点与第j个输出层结点之间的连接权重，h和m分别表示隐含层和输出层结点个数；极限学习机是一种单隐层前馈神经网络，给定k个训练样本{xk,yk},k＝1,2,...,k，其中x＝[x1,x2,...,xk]^t∈r^m为输入向量，y＝[y1,y2,...,yk]^t∈rⁿ为输出向量；则第k个训练样本的第j个隐含层神经元的输出表示为其中，j表示隐含神经元数目，g(·)表示激活函数，βj和aj分别表示输入层和隐含层之间想连接权重和偏置；elm模型第k个训练样本的输出可以表示为：式中，wj表示隐含层和输出层之间的权重；隐含层核映射矩阵表示如下：

elm网络的输出表示为：(11)；式中，w＝[w1,w2,...,wj]^t；elm模型的目标是找到最合适的w，使得网络输出值和实测值y之间的误差最小，通过求解以下优化问题得到隐含层与输出层之间的关联权重：

方程(12)的解根据moorepenrose(mp)广义逆求得：

式中，是矩阵d的广义逆矩阵；则elm模型第k个训练样本的输出表示为：elman神经网络是一种动态递归神经网络。除了输入层、隐含层和输出层外，elman神经网络还包括一个特殊的递归层，用于记忆隐含层神经元以前时刻的输出信息，然后将该信息用作隐藏层的输入。

进一步地，在所述步骤五中，假设已知随机向量x，y和联合概率密度函数p(x,y)，若x的观测值为x，则y关于x的回归表示为：式中，e[y|x]为给定输入x下输出y的预测值，密度函数p(x,y)未知，由parzen非参数估计进行估算：式中，x^j和y^j分别表示第j个样本的实测输入值和实测输出值，n是样本集的个数，d是x的维数，σ是平滑参数；根据式(15)和(16)，grnn模型的输出可以表示如下：

优选地，在所述步骤六中，rmse，mae，mape及r的计算公式如下：

式中，qp(i)为预测值；qo(i)为实测值；和分别表示为实测值和预测值的均值；n是样本集的个数。

有益效果：与现有技术相比，通过本发明建立的中长期径流预报方法，能够取得下列有益效果：

1)本发明采用经验小波变换对径流时间序列进行分解与重构，能够去除天然径流时间序列中的冗余噪声，降低径流时间序列的非线性波动对径流预报的干扰，进而提高预报模型的性能；

2)针对单一ann模型的预报精度难以达到预期水平，且每种模型各具优缺点，以三种ann预报模型的预报结果作为输入因子，采用grnn神经网络模型对月径流时间序列进行进一步预报，以综合不同ann预报模型的优势，实现各基模型预报信息的融合，从而获得精度更高的径流预报信息。

附图说明

图1为基于经验小波消噪和神经网络融合的中长期径流预报模型流程图；

图2为乌东德站径流时间序列的ewt分解以及原始与消噪径流时间序列的对比图；

图3为乌东德站检验期gne和ewt-gne模型的月径流预测结果对比图。

具体实施方式

下面结合附图对本发明进行详细的介绍。

实施方式1：

本发明以长江上游流域乌东德水文站月径流时间序列为实施例，进行实例仿真，以验证本发明的效果。图1为本发明提供的基于经验小波消噪和神经网络融合的中长期径流预报方法流程图，实施步骤如下：

步骤一：获取乌东德水文站1958年1月至2012年12月的历史实测月径流数据(660个样本数据点)，建立月径流时间序列，并将前525个月(1958年1月至2001年9月)的月径流实测值作为训练样本，后135个月(2001年10月至2012年12月)的月径流实测值作为检验样本。

步骤二：采用ewt将乌东德站月径流时间序列分解为几个独立的经验模态，去除分解经验模态中频率最高的模态，并对剩余的经验模态进行线性求和，得到重构径流时间序列。乌东德站径流时间序列的ewt分解以及原始与消噪径流时间序列的对比图如图2所示。

ewt的主要思想是通过对原始径流时间序列的fourier频谱进行自适应划分，并建立一组合适的小波滤波器筛选出原始径流时间序列的调幅调频(am-fm)成分。对于原始径流时间序列x(t)，ewt的分解过程如下：

(1)根据快速傅立叶变换算法(fft)计算原始径流时间序列x(t)的傅立叶频谱f(ω)。

(2)将傅立叶频谱f(ω)的频域[0,π]自适应的分割为k个带宽不等的频带，[0,ω1],[ω1,ω2],...,[ωk-1,π]，其中，小波滤波器的边界ωn,n＝1,2,...,n-1选取为两个连续局部极大值之间的中间频率。

(3)根据meyer小波构造经验小波，确定经验小波函数ψk(ω)和经验尺度函数

式中，β(x)＝x⁴(35-84x+70x²-20x³)。

(4)重构原始径流时间序列得到不同的模态。原始径流时间序列通过下式进行重构：

式中，*为卷积运算，wx(0,t)为近似系数，wx(k,t)为细节系数，近似系数和细节系数分别通过下式给出：

式中，ψk(t)和ψk(t)分别表示经验小波函数和经验尺度函数，和分别表示ψk(τ-t)和的复共轭，f^-1(·)表示傅里叶逆变换，和分别表示ψ(ω)和的傅里叶变换。

则经验模态uk(t)的定义如下：

步骤三：对经ewt处理后的消噪月径流时间序列进行相空间重构，构造相空间矩阵作为基预报模型的输入因子，建立ann基预报模型的输入、输出矩阵。

步骤四：采用步骤三确定的训练期输入、输出矩阵分别对rbf、elm和elman三种神经网络模型进行训练。

rbf神经网络是一种多层前馈神经网络(fnn)。rbf神经网络包括三层网络结构：由输入变量组成的输入层、对输入变量进行非线性变换的隐含层和生成网络线性响应的输出层。rbf神经网络的隐含层使用高斯传递函数作为激活函数。高斯激活函数定义为：

式中，x＝[x1,x2,...,xn]^t为n维的输入向量，ci＝[ci1,ci2,...,cin]^t为第i个隐含层神经元的中心，qi表示高斯函数的宽度，·为欧氏范数。则输出层第j个结点的响应为：

式中，wij为第i个隐含层结点与第j个输出层结点之间的连接权重，h和m分别表示隐含层和输出层结点个数。

极限学习机是一种单隐层前馈神经网络，给定k个训练样本{xk,yk},k＝1,2,...,k，其中x＝[x1,x2,...,xk]^t∈r^m为输入向量，y＝[y1,y2,...,yk]^t∈rⁿ为输出向量。则第k个训练样本的第j个隐含层神经元的输出可以表示为其中，j表示隐含神经元数目，g(·)表示激活函数，βj和aj分别表示输入层和隐含层之间想连接权重和偏置。elm模型第k个训练样本的输出可以表示为：

式中，wj表示隐含层和输出层之间的权重。

隐含层核映射矩阵可以表示如下：

elm网络的输出可以表示为：

式中，w＝[w1,w2,...,wj]^t

elm模型的目标是找到最合适的w，使得网络输出值和实测值y之间的误差最小。可以通过求解以下优化问题得到隐含层与输出层之间的关联权重：

方程(12)的解可以根据moorepenrose(mp)广义逆求得：

式中，是矩阵d的广义逆矩阵。

则elm模型第k个训练样本的输出可以表示为：

elman神经网络是一种动态递归神经网络。除了输入层、隐含层和输出层外，elman神经网络还包括一个特殊的递归层，用于记忆隐含层神经元以前时刻的输出信息，然后将该信息用作隐藏层的输入。

步骤五：将三种ann基预报模型的输出作为grnn模型的预报因子，再次建立输入、输出矩阵，对grnn模型进行训练，并将检验样本代入训练好的grnn模型，得到检验期的预测值。

grnn神经网络是rbf神经网络的一种特殊形式。假设已知随机向量x，y，和联合概率密度函数p(x,y)，若x的观测值为x，则y关于x的回归可以表示为：

式中，e[y|x]为给定输入x下输出y的预测值。密度函数p(x,y)通常是未知的，可以由parzen非参数估计进行估算：

式中，x^j和y^j分别表示第j个样本的实测输入值和实测输出值，n是样本集的个数，d是x的维数，σ是平滑参数。

根据式(15)和(16)，grnn模型的输出可以表示如下：

步骤六：采用四种评价指标包括均方根误差(rmse)、平均绝对误差(mae)、平均绝对百分误差(mape)和相关系数(r)评价对预报结果进行评价。

rmse，mae，mape及r的计算公式如下：

式中，qp(i)为预测值；qo(i)为实测值；和分别表示为实测值和预测值的均值；n是样本集的个数。

采用本发明所提的基于经验小波消噪和神经网络融合(ewt-gne)模型对乌东德水文站月径流时间序列进行预报，并与三种单一ann模型(rbf，elm和elman)以及基于grnn的神经网络融合预报模型(gne)进行对比，乌东德水文站月径流时间序列训练期和检验期的预报结果如表1所示，乌东德站检验期gne和ewt-gne模型的月径流预测结果对比图如图3所示。

表1乌东德站月径流预报结果误差统计

由表1可知，神经网络融合gne模型的rmse，mae和mape比三种单一ann预报模型小，且r比三种ann模型大，说明神经网络融合模型能够综合三种单一ann模型的优势，提高单一ann预报模型的性能。ewt-gne模型的预报效果优于gne模型，说明ewt能够有效地去除月径流时间序列中的冗余噪声，显著提高gne预报模型的精度。由图3可知，相比于gne模型，ewt-gne模型的预报流量与实测流量值更为贴近，进一步证明了ewt-gne模型的优越性。

上述实施方式只为说明本发明的技术构思及特点，其目的在于让熟悉此项技术的人能够了解本发明的内容并据以实施，并不能以此限制本发明的保护范围。凡根据本发明精神实质所做的等效变换或修饰，都应涵盖在本发明的保护范围之内。