CN111553513A - 一种基于二次分解与回声状态网络的中长期径流预测方法 - Google Patents

Abstract

本发明公开了一种基于二次分解与回声状态网络（ESN）的中长期径流预测方法，属于径流预测技术领域。该方法包括：步骤一：获取径流序列x(t)，根据径流序列的数据情况，将其分为训练样本和测试样本；步骤二：利用自适应噪声完备经验模态分解（CEEMDAN）将径流序列分解为若干固有模态函数（IMF）和一个趋势项（Res）；步骤三：使用变分模态分解（VMD）方法对频率最高的IMF分量进行二次分解，得到若干变分模态（Mode）；步骤四：将两次分解出的子序列分别输入ESN进行预测，并将各预测结果重构得到最终预测值。本发明解决了对一次分解中高频分量预测误差较大的问题，进一步提高ESN的预测精度。

Description

一种基于二次分解与回声状态网络的中长期径流预测方法

技术领域

本发明涉及一种基于二次分解与回声状态网络(ESN)的中长期径流预测方法，属于径流预测技术领域。

背景技术

目前径流预测方法使用较为广泛的是主要以人工神经网络为代表的人工智能方法，但传统的递归神经网络具有记忆渐消的缺点。因此，具有良好短期记忆功能、简单训练算法以及容易确定网络结构的ESN获得越来越多研究人员的关注。然而，径流序列的复杂性影响了ESN模型的参数率定及学习效率，单一模型不能充分挖掘并学习径流序列中所包含的全部信息，预测精度提升空间有限，从而出现了“分解-预测-重构”模式的组合预测方法，通过将径流序列分解成若干子序列的方式，简化原始径流序列。但这种模式仍有一定缺陷，分解出的最高频分量包含了原始径流序列的主要特征，具有波动强烈且振动幅值大的特点，序列仍较为复杂，导致最高频分量的预测精度不足，进而对整个径流序列预测效果产生较大的影响。

因此，寻求新的组合模式来减小最高频分量预测误差，以期提高整体预测精度，已成为水文预报领域的热点问题。

发明内容

本发明旨在提供一种基于二次分解与回声状态网络(ESN)的中长期径流预测方法，解决了一次分解中最高频分量仍波动强烈的问题，进一步提高了性能较好的模型预测精度，对中长期径流预测具有很好的参考价值。

本发明提供了一种基于二次分解与回声状态网络(ESN)的中长期径流预测方法，将CEEMDAN-VMD与ESN模型组合使用，首先通过CEEMDAN方法，将径流序列分解出一系列包含不同尺度局部信息的子序列，以解决径流的复杂性影响模型参数确定及预测效率；再采用VMD对子序列中频率最高的分量进行二次频域特性分解，降低最高频分量的波动程度，在此基础上，采用预测性能较好的ESN模型，从而综合提高预测精度，对中长期径流预测具有借鉴意义。

上述内容使用CEEMDAN可有效地将复杂的径流序列分解为包含不同尺度时域信息的子序列，可避免EMD分解出现的模态混叠问题，与集合经验模态分解(EEMD)相比，提高了分解效率，同时解决互补集合经验模态分解(CEEMD)添加不同噪声信号产生不同数量模态分解的问题。该方法实现了分解信号的近似完美重构，解决了不同加噪信号产生不同数量模态的缺陷；然而其频率最高的IMF分量包含主要的原始序列频率信息，波动较为强烈，序列仍较为复杂，导致预测精度较低，借鉴分解方法对原始径流序列的处理思路，采用VMD方法对最高频分量进行频域特性分解，进一步提取最高频分量的不同频域信息，降低序列复杂程度，以期进一步提高ESN模型预测精度；而使用的ESN模型，其作为传统循环神经网络(RNN)的改进方法，具有独特的动态储备池结构，通过调整储备池中神经元状态的数值特征，使得系统具有良好的短期记忆功能及非线性逼近能力，且简化了网络训练过程，提高计算效率。

上述基于二次分解与回声状态网络的中长期径流预测方法，具体包括以下步骤：

步骤一：获取径流序列x(t)，根据径流序列的数据情况，将其分为训练样本和测试样本；

步骤二：利用CEEMDAN将径流序列分解为若干IMF分量和一个Res；

步骤三：使用VMD方法对频率最高的IMF分量进行二次分解，得到若干Mode；

步骤四：将两次分解出的子序列分别输入ESN进行预测；

步骤五：将各子序列的预测结果重构，得到最终径流序列预测值。

下面对各步骤具体说明如下：

步骤二所述利用CEEMDAN将径流序列分解为若干IMF分量和一个Res，具体包括以下步骤：

(1)对原始径流序列x(t)添加服从标准正态分布的白噪声序列：

x_m(t)＝x(t)+ε₀σ_m(t)m＝1,2,…,M

式中，x_m(t)为加入白噪声后的序列；ε₀为适用于x_m(t)的信噪比；M为集总次数，即实现M次经验模态分解EMD；σ_m(t)为服从标准正态分布的白噪声；

(2)对时间序列x_m(t)进行M次EMD分解，并计算总体平均值，将其定义为原始信号的第一个IMF分量，即：

式中，c_m为EMD分解得到的第m个IMF分量；c₁为第一个IMF分量；

(3)计算一阶残差r₁(t)：

r₁(t)＝x(t)-c₁(t)

(4)对一阶残差添加白噪声后再次进行EMD分解，求取平均值得到第二个IMF分量c₂(t)：

式中，ε₁为适用于r₁(t)的信噪比；E₁(·)为EMD算法产生第1个IMF分量的算子；

(5)对后面的每个阶段，k＝2,3,…,K得到的第k个余量信号r_k(t)与第(3)步相同，再重复第(4)步计算第k+1个IMF分量，即：

r_k(t)＝r_k-1(t)-c_k(t)

式中，ε_k为适用于r_k(t)的信噪比；E_k(·)为EMD算法产生第k+1个IMF分量的算子；

(6)直至残差不能再分解，从而求得最终残差R(t)：

(7)则原始序列x(t)的表达式为：

式中，K代表在给定精度下原始序列分解出的IMF分量总个数。

步骤三所述使用VMD方法对频率最高的IMF分量进行二次分解，得到若干Mode，具体包括以下步骤：

(1)采用中心频率观察法确定VMD分解出的变分模态Mode的个数r，令r＝1,2,…,R，采用VMD对IMF中的最高频分量进行二次分解，若当r＝R+1时，VMD分解出的Mode中心频率接近，则r＝R为最终选取的Mode个数；

(2)对于IMF分量中的最高频分量采用Hilbert变换，得到原始Mode的解析信号，并得到单边频谱：

式中，δ(t)为单位脉冲函数；j表示虚数单位；u_r(t)为第r个Mode；*表示卷积；

(3)将各个Mode的频谱调解到基带，即对

乘以

来调整整个u_r(t)的估算中心频率，将各个Mode调解到基频带；

(4)利用约束变分模型求解每个Mode分量的带宽：

同时满足：

式中，ω_k为各个Mode的角频率；c₁(t)为被二次分解的最高频分量；

表示对t求偏导；

(5)通过引入二次惩罚因子α和拉格朗日乘法算子λ(t)求解约束变分模型的最优解，即可得到最优Mode；其中二次惩罚因子α能保证信号的重构精度，拉格朗日乘法算子λ(t)能加强约束。

步骤四所述将两次分解出的子序列分别输入ESN进行预测，具体包括以下步骤：将步骤二、步骤三分解出的各子序列，分为与原始径流序列长度相同的训练样本及测试样本，再将各子序列训练样本输入ESN进行模型训练，寻找最优参数，最后采用建立好的模型预测各子序列的测试样本。

进一步地，所述的ESN建立过程为：

(1)训练阶段

①初始化ESN网络：

对输入连接权值Wⁱⁿ、反馈连接权值W^back和内部状态向量W进行初始化，设定储备池的内部各个参数，构建动态储备池，并更新内部状态向量：

x(n)＝f(Wⁱⁿu(n)+Wx(n-1)+W^backy(n))n＝2,3,…,N

式中，u(n)＝u₁(n),u₂(n),…,u_P(n)为ESN的输入向量；x(n)＝x₁(n),x₂(n),…,x_N(n)为ESN的状态向量；y(n)＝y₁(n),y₂(n),…,y_L(n)为ESN的输出向量；P为输入向量个数；N为储备池规模；L为输出向量个数；f为储备池处理单元的内部神经元激活函数，常用tanh函数；

②收集内部状态向量B与相应输出信号T，以此求得输出权重W^out：

W^out＝B^-1T

③根据求出的输出权重W^out建立高效的ESN预测模型；

(2)测试阶段根据训练阶段得到网络输出权重W^out求得最终预测结果：

y(n)＝f^out(W^outu(n),x(n),y(n-1))n＝2,3,…,N

式中，y(n)＝y₁(n),y₂(n),…,y_L(n)为输出向量；L为输入向量个数；f^out为输出处理单元的内部神经元的激活函数，常用tanh函数。

本发明将CEEMDAN-VMD与ESN模型组合使用，带来的有益效果：

(1)使用CEEMDAN可有效地将复杂的径流序列分解为包含不同尺度时域信息的子序列，可避免EMD分解出现的模态混叠问题，与集合经验模态分解(EEMD)相比，提高了分解效率，同时解决互补集合经验模态分解(CEEMD)添加不同噪声信号产生不同数量模态分解的问题。

(2)采用VMD方法对最高频分量进行频域特性分解，进一步提取最高频分量的不同频域信息，降低序列复杂程度，以期进一步提高ESN模型预测精度；

(3)使用的ESN模型，其作为传统循环神经网络(RNN)的改进方法，具有独特的动态储备池结构，通过调整储备池中神经元状态的数值特征，使得系统具有良好的短期记忆功能及非线性逼近能力，且简化了网络训练过程，提高计算效率。

附图说明

图1为本发明的流程图；

图2为ESN的结构图。

具体实施方式

下面结合附图以及具体实施方式进一步说明本发明。

本发明提供了一种基于二次分解与回声状态网络的中长期径流预测方法，该方法包括：

步骤一：获取径流序列x(t),根据径流序列的数据情况，将其分为训练样本和测试样本；

步骤二：利用CEEMDAN将径流序列分解为若干IMF分量和一个Res；

步骤三：使用VMD方法对频率最高的IMF分量进行二次分解，得到若干Mode；

步骤四：将两次分解出的子序列分别输入ESN进行预测；

步骤五：将各子序列的预测结果重构，得到最终径流序列预测值。

图1为本发明的流程图，根据图1所示，本发明基于二次分解与回声状态网络的中长期径流预测方法流程介绍如下：

该方法主要包括一次分解、二次分解、预测及重构三个部分：(一)利用CEEMDAN将径流序列分解成一系列包含不同尺度时域信息的IMF分量及一个Res，通过快速傅里叶变换(FFT)对各IMF分量进行频谱分析，检验CEEMDAN对径流序列信息的提取效果；(二)将IMF分量中的最高频分量通过VMD进一步分解，得到若干Mode，再次通过FFT对各Mode分量进行频谱分析，检验VMD对最高频分量的信息提取及简化效果；(三)将两次分解得到的子序列通过ESN预测及重构，得到最终径流序列预测值。

所述的第(一)部分中，CEEMDAN是以EMD方法为基础，是对EEMD方法的重要发展，该方法解决了EMD的模态混叠，EEMD不能精确重建，且计算效率较低，以及CEEMD分解过程中添加不同噪声信号产生不同数量模态分解等问题。CEEMDAN分解过程通过以下步骤实现：

(1)对原始径流序列x(t)添加服从标准正态分布的白噪声序列：

x_m(t)＝x(t)+ε₀σ_m(t)m＝1,2,…,M

式中，x_m(t)为加入白噪声后的序列；ε₀为适用于x_m(t)的信噪比；M为集总次数，即实现M次EMD分解；σ_m(t)为服从标准正态分布的白噪声。

(2)对时间序列x_m(t)进行M次EMD分解，并计算总体平均值，将其定义为原始信号的第一个IMF分量，即：

式中，c_m为EMD分解得到的第m个IMF分量；c₁为第一个IMF分量。

(3)计算一阶残差r₁(t)：

r₁(t)＝x(t)-c₁(t)

(4)对一阶残差添加白噪声后再次进行EMD分解，求取平均值得到第二个IMF分量c₂(t)：

式中，ε₁为适用于r₁(t)的信噪比；E₁(·)为EMD算法产生第1个IMF分量的算子。

(5)对后面的每个阶段，k＝2,3,…,K得到的第k个余量信号r_k(t)与第(3)步相同，再重复第(4)步计算第k+1个IMF分量，即：

r_k(t)＝r_k-1(t)-c_k(t)

式中，ε_k为适用于r_k(t)的信噪比；E_k(·)为EMD算法产生第k+1个IMF分量的算子。

(6)直至残差不能再分解，从而求得最终残差R(t)：

(7)则原始序列x(t)的表达式为：

式中，K代表在给定精度下原始序列分解出的IMF分量总个数。将分解后的IMF分量进行FFT检验CEEMDAN提取不同尺度信息以简化原始径流序列的效果。

所述的第(二)部分中，将IMF分量中的最高频分量通过VMD进一步分解，通过以下步骤实现：

(1)采用中心频率观察法确定VMD分解出的变分模态Mode的个数r，令r＝1,2,…,R，采用VMD对IMF中的最高频分量进行二次分解，若当r＝R+1时，VMD分解出的Mode中心频率接近，则r＝R为最终选取的Mode个数；

(2)对于IMF分量中的最高频分量采用Hilbert变换，得到原始Mode的解析信号，并得到单边频谱：

式中，δ(t)为单位脉冲函数；j表示虚数单位；u_r(t)为第r个Mode；*表示卷积。

(3)将各个Mode的频谱调解到基带，即对

乘以

来调整整个u_r(t)的估算中心频率，将各个Mode调解到基频带；

(4)利用约束变分模型求解每个Mode分量的带宽：

同时满足：

式中，ω_k为各个Mode的角频率；c₁(t)为被二次分解的最高频分量；

表示对t求偏导。

(5)通过引入二次惩罚因子α和拉格朗日乘法算子λ(t)求解约束变分模型的最优解，即可得到最优Mode。其中二次惩罚因子α可保证信号的重构精度，拉格朗日乘法算子λ(t)可加强约束。

将分解后的Mode再次进行FFT变换，检验VMD进一步提取最高频分量中不同频率信息的效果。

所述的第(三)部分中，将各IMF分量、Res及Mode通过ESN预测及重构，得到最终径流序列预测值，具体过程为：

(1)对于每个分量计算不同时间延迟下的偏自相关系数(PACF)，选择PACF值在95％置信水平之外的历史径流序列作为预测输入变量。

(2)训练阶段

基本的ESN拓扑结构如图2所示。图2中左侧为P个输入节点；中间是储备池网络，由N个内部节点以及稀疏的节点连接权值矩阵W构成；右侧是L个输出节点。图中虚线箭头表示不同情况下还可能存在的连接，即进行误差调节时的信息反馈。在ESN训练阶段中，只有输出权重会改变。具体训练过程如下：

1)初始化ESN网络。对输入连接权值Wⁱⁿ、反馈连接权值W^back和内部状态向量W进行初始化，设定储备池的内部各个参数，构建动态储备池，并更新内部状态向量：

x(n)＝f(Wⁱⁿu(n)+Wx(n)+W^backy(n))n＝2,3,…,N

式中，u(n)＝u₁(n),u₂(n),…,u_P(n)为ESN的输入向量；x(n)＝x₁(n),x₂(n),…,x_N(n)为ESN的状态向量；y(n)＝y₁(n),y₂(n),…,y_L(n)为ESN的输出向量；P为输入向量个数；N为储备池规模；L为输出向量个数；f为储备池处理单元的内部神经元激活函数，常用tanh函数。

2)收集内部状态向量B与相应输出信号T，以此求得输出权重W^out：

W^out＝B^-1T

3)根据求出的输出权重W^out建立高效的ESN预测模型。

(3)测试阶段根据训练阶段得到网络输出权重W^out求得最终预测结果：

y(n)＝f^out(W^outu(n),x(n),y(n-1))n＝2,3,…,N

式中，y(n)＝y₁(n),y₂(n),…,y_L(n)为输出向量；L为输入向量个数；f^out为输出处理单元的内部神经元的激活函数，常用tanh函数。

(4)将各分量及趋势项的预测结果重构，即为原径流序列测试样本最终预测值。

Claims (5)

1.一种基于二次分解与回声状态网络的中长期径流预测方法，其特征在于：包括以下步骤：

步骤一：获取径流序列x(t)，根据径流序列的数据情况，将其分为训练样本和测试样本；

步骤二：利用CEEMDAN将径流序列分解为若干IMF分量和一个Res；

步骤三：使用VMD方法对频率最高的IMF分量进行二次分解，得到若干Mode；

步骤四：将两次分解出的子序列分别输入ESN进行预测；

步骤五：将各子序列的预测结果重构，得到最终径流序列预测值。

2.根据权利要求1所述的基于二次分解与回声状态网络的中长期径流预测方法，其特征在于：步骤二所述利用CEEMDAN将径流序列分解为若干IMF分量和一个Res，具体包括以下步骤：

(1)对原始径流序列x(t)添加服从标准正态分布的白噪声序列：

x_m(t)＝x(t)+ε₀σ_m(t) m＝1,2,…,M

式中，x_m(t)为加入白噪声后的序列；ε₀为适用于x_m(t)的信噪比；M为集总次数，即实现M次经验模态分解EMD；σ_m(t)为服从标准正态分布的白噪声；

(2)对时间序列x_m(t)进行M次EMD分解，并计算总体平均值，将其定义为原始信号的第一个IMF分量，即：

式中，c_m为EMD分解得到的第m个IMF分量；c₁为第一个IMF分量；

(3)计算一阶残差r₁(t)：

r₁(t)＝x(t)-c₁(t)

(4)对一阶残差添加白噪声后再次进行EMD分解，求取平均值得到第二个IMF分量c₂(t)：

式中，ε₁为适用于r₁(t)的信噪比；E₁(·)为EMD算法产生第1个IMF分量的算子；

(5)对后面的每个阶段，k＝2,3,…,K得到的第k个余量信号r_k(t)与第(3)步相同，再重复第(4)步计算第k+1个IMF分量，即：

r_k(t)＝r_k-1(t)-c_k(t)

式中，ε_k为适用于r_k(t)的信噪比；E_k(·)为EMD算法产生第k+1个IMF分量的算子；

(6)直至残差不能再分解，从而求得最终残差R(t)：

(7)则原始序列x(t)的表达式为：

式中，K代表在给定精度下原始序列分解出的IMF分量总个数。

3.根据权利要求1所述的基于二次分解与回声状态网络的中长期径流预测方法，其特征在于：步骤三所述使用VMD方法对频率最高的IMF分量进行二次分解，得到若干Mode，具体包括以下步骤：

(1)采用中心频率观察法确定VMD分解出的变分模态Mode的个数r，令r＝1,2,…,R，采用VMD对IMF中的最高频分量进行二次分解，若当r＝R+1时，VMD分解出的Mode中心频率接近，则r＝R为最终选取的Mode个数；

(2)对于IMF分量中的最高频分量采用Hilbert变换，得到原始Mode的解析信号，并得到单边频谱：

式中，δ(t)为单位脉冲函数；j表示虚数单位；u_r(t)为第r个Mode；*表示卷积；

(3)将各个Mode的频谱调解到基带，即对

乘以

来调整整个u_r(t)的估算中心频率，将各个Mode调解到基频带；

(4)利用约束变分模型求解每个Mode分量的带宽：

同时满足：

式中，ω_k为各个Mode的角频率；c₁(t)为被二次分解的最高频分量；

表示对t求偏导；

(5)通过引入二次惩罚因子α和拉格朗日乘法算子λ(t)求解约束变分模型的最优解，即可得到最优Mode；其中二次惩罚因子α能保证信号的重构精度，拉格朗日乘法算子λ(t)能加强约束。

4.根据权利要求1所述的基于二次分解与回声状态网络的中长期径流预测方法，其特征在于：步骤四所述将两次分解出的子序列分别输入ESN进行预测，具体包括以下步骤：将步骤二、步骤三分解出的各子序列，分为与原始径流序列长度相同的训练样本及测试样本，再将各子序列训练样本输入ESN进行模型训练，寻找最优参数，最后采用建立好的模型预测各子序列的测试样本。

5.根据权利要求4所述的基于二次分解与回声状态网络的中长期径流预测方法，其特征在于：所述的ESN建立过程为：

(1)训练阶段

①初始化ESN网络：

对输入连接权值Wⁱⁿ、反馈连接权值W^back和内部状态向量W进行初始化，设定储备池的内部各个参数，构建动态储备池，并更新内部状态向量：

x(n)＝f(Wⁱⁿu(n)+Wx(n-1)+W^backy(n)) n＝2,3,…,N

式中，u(n)＝u₁(n),u₂(n),…,u_P(n)为ESN的输入向量；x(n)＝x₁(n),x₂(n),…,x_N(n)为ESN的状态向量；y(n)＝y₁(n),y₂(n),…,y_L(n)为ESN的输出向量；P为输入向量个数；N为储备池规模；L为输出向量个数；f为储备池处理单元的内部神经元激活函数，常用tanh函数；

②收集内部状态向量B与相应输出信号T，以此求得输出权重W^out：

W^out＝B^-1T

③根据求出的输出权重W^out建立高效的ESN预测模型；

(2)测试阶段根据训练阶段得到网络输出权重W^out求得最终预测结果：

y(n)＝f^out(W^outu(n),x(n),y(n-1)) n＝2,3,…,N

式中，y(n)＝y₁(n),y₂(n),…,y_L(n)为输出向量；L为输入向量个数；f^out为输出处理单元的内部神经元的激活函数，常用tanh函数。

Images