一种集合经验模态分解及遥相关的中长期径流预报方法与流程

本发明涉及径流预报技术领域，尤其涉及一种集合经验模态分解及遥相关的中长期径流预报方法。

背景技术：

中长期径流预报是水资源配置方案、水库水量调度计划最重要的数据基础和制定依据，是实现水资源高效利用、防洪减灾的关键技术，径流预报的精度同时又给库区来水、调度风险带来不确定性。径流序列受大气变化、下垫面情况、人类活动等多因素影响，径流变化多样，影响因素复杂，呈非稳态，具有一定的周期性，随机性和趋势性，故准确预报径流并掌握其物理机理具有较大的难度。

目前，常用的中长期径流预报方法主要有：基于数据统计分析的自回归模型预报和气象因子遥相关分析预报两大类方法。其中预报模型通常使用人工神经网络、支持向量机、线性回归等智能优化算法对径流进行预测。现有的两类方法建立的预报模型都将径流序列假定为稳态序列，从而导致预报精度不高。另外，自回归预报模型忽略了径流过程中的水文物理机理，使用智能算法对径流进行自相关分析，会导致预报过拟合，出现失真现象。

技术实现要素：

本发明的目的在于提供一种集合经验模态分解及遥相关的中长期径流预报方法，从而解决现有技术中存在的前述问题。

为了实现上述目的，本发明采用的技术方案如下：

一种集合经验模态分解及遥相关的中长期径流预报方法，包括如下步骤：

s1，提取径流序列及响应的预报因子数据，从径流序列q中剔除趋势项trend生成新的无趋势变化的径流序列q'，将滞后时段t的遥相关数据d组成一系列的气象因子集合x(d,t)；

s2，采用集合经验模态分解的方法将径流序列q'分解为满足hilbert变换要求的n层本征模态函数imf(1),imf(2),…,imf(n)和一层残余项res，组成一组自变量集合y＝[imf(n),res]；

s3，确定自变量集合y相关的预报因子f(x,n+1)，包括如下步骤：

s301，对不同的气象因子x(d,t)与每层本征模态函数和残余项res进行相关性分析，得到相关因子；

s302，根据得到的相关因子，判断x(d,t)和y＝[imf(n),res]中每组序列的相关性，通过t检验，初步筛选每组自变量的相关性较大的气象因子；

s303，根据逐步回归方法，挑选显著影响因子，最终确定每组自变量的预报因子f(x,n+1)，完成对n阶自变量y的预报因子识别工作；

s4，利用筛选出来的预报因子f(x,n+1)和已知径流分解后的imf(1),imf(2),…,imf(n)、res构建基于集合经验模态分解与遥相关数据的径流预报模型；

s5，根据所述径流预报模型对下一时段的中长期径流进行预报。

优选地，s4包括如下步骤：

s401，构建径流预报模型的训练、验证样本集：根据样本数据长度，确定模型训练期和检验期长度，并将输入因变量f(x,n+1)为已识别的气象因子及确定的滞后期，自变量y为n层本征模态函数imf(1),imf(2),…,imf(n)、残余项res分为训练样本集和检验期样本集两段；

s402，基于已获得的训练样本集，建立并训练人工神经网络模型，进而预报检验期的自变量y；

s403，对人工神经网络模型训练期和检验期对自变量y的预报效果进行评价，若评价效果符合设定的阈值，则径流预报模型建立完成，该模型用于对未来时段径流进行预报，否则，重新调整人工神经网络的模型参数，重新建模直至模型完成建立。

优选地，s5具体为：

使用验证好的模型对下一时段的自变量y进行预报，同时对下一时段的趋势项trendf进行预测，其预测径流的公式为：

其中，qf为预测的下一时段径流量，imff(j)为第j层本征模态函数分量的预测值，resf为预测的残余量。

优选地，s403中，所述对人工神经网络模型训练期和检验期对自变量y的预报效果进行评价，，评价的指标包括纳什效率系数、相对误差和合格率；

所述纳什效率系数nash按照如下公式进行计算：

其中，qo为实测径流量，qf为预测径流量，qo为实测径流量均值，nash越接近1，预报越精准；

所述相对误差mare按照如下公式进行计算：

其中，n为qo的序列长度，其余参数含义同上式。mare越接近0，说明实测与预报值越接近，预报效果越精准，常认为mare＜20％时，效果较好；

所述合格率qr按照如下公式进行计算：

其中，n为合格预报次数，m为预报总次数，当qr＞80％时，认为预报效果较好。

优选地，s1中，所述滞后时段t的遥相关数据d，包括130项遥相关数据，其中，88项环流指数，26项海温数据，16项其他因子数据；滞后时段t主要考虑滞后1～12个月，共同组成一组130×12的气象因子矩阵x(d,t)。

优选地，s2中，所述n层本征模态函数项imf(j)和一层残余项res的关系式为：

其中，q'为连续多年的除趋势项的径流序列，imf(j)为第j层本征模态函数分量，n为本征模态函数的分量数。

优选地，s301中，所述相关性分析，通过皮尔逊相关系数判断，具体采用如下公式：

其中，r(x,y)为相关因子，cov(x,y)为x和y的协方差，var[x]为x的方差var[y]为y的方差，|r(x,y)|越接近1，x和y相关性越大。

本发明的有益效果是：本发明实施例提供的集合经验模态分解及遥相关的中长期径流预报方法，与现有技术相比，充分考虑了径流序列的非稳态性，使用集合经验模态分解的方法将径流转化为多组稳态序列，实现了水文序列稳态化，为常规的预报方法提供了最基础的数据条件。同时，针对常规自适应中长期径流预报在物理机理方面的不足，提出采用遥相关数据分析的方法，分析多组稳态序列的气象影响因子，能充分剥离不同频率序列的物理因子，挖掘径流序列隐藏的物理机制，建立的预报模型提高了中长期径流预报的准确性。

附图说明

图1为本发明提供的集合经验模态分解及遥相关的中长期径流预报方法流程示意图；

图2为本发明采用的集合经验模态分解法流程示意图；

图3为基于集合经验模态分解方法对径流序列进行分解后的第一层序列示意图；

图4为基于集合经验模态分解与遥相关分析方法预测中长期径流的效果示意图。

具体实施方式

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

针对中长期径流序列的非稳态的特点，为了克服自相关预报缺乏物理机理的不足，预报方法的基础数据假设的缺陷，本发明提供一种基于集合经验模态分解及遥相关分析的中长期径流预报方法，有效地将水文序列稳态化，并且通过遥相关数据分析，识别了影响径流序列的气象因子，挖掘了其物理意义，提高了中长期径流预报的准确性。

如图1所示，本发明实施例提供了一种集合经验模态分解及遥相关的中长期径流预报方法，包括如下步骤：

s1，提取径流序列及响应的预报因子数据，从径流序列q中剔除趋势项trend生成新的无趋势变化的径流序列q'，将滞后时段t的遥相关数据d组成一系列的气象因子集合x(d,t)；

s2，采用集合经验模态分解的方法将径流序列q'分解为满足hilbert变换要求的n层本征模态函数imf(1),imf(2),…,imf(n)和一层残余项res，组成一组自变量集合y＝[imf(n),res]；

s3，确定自变量集合y相关的预报因子f(x,n+1)，包括如下步骤：

s301，对不同的气象因子x(d,t)与每层本征模态函数和残余项res进行相关性分析，得到相关因子；

s302，根据得到的相关因子，判断x(d,t)和y＝[imf(n),res]中每组序列的相关性，通过t检验，初步筛选每组自变量的相关性较大的气象因子；

s303，根据逐步回归方法，挑选显著影响因子，最终确定每组自变量的预报因子f(x,n+1)，完成对n阶自变量y的预报因子识别工作；

s4，利用筛选出来的预报因子f(x,n+1)和已知径流分解后的imf(1),imf(2),…,imf(n)、res构建基于集合经验模态分解与遥相关数据的径流预报模型；

s5，根据所述径流预报模型对下一时段的中长期径流进行预报。

其中，s4可以包括如下步骤：

s401，构建径流预报模型的训练、验证样本集：根据样本数据长度，确定模型训练期和检验期长度，并将输入因变量f(x,n+1)为已识别的气象因子及确定的滞后期，自变量y为n层本征模态函数imf(1),imf(2),…,imf(n)、残余项res分为训练样本集和检验期样本集两段；

s402，基于已获得的训练样本集，建立并训练人工神经网络模型，进而预报检验期的自变量y；

s403，对人工神经网络模型训练期和检验期对自变量y的预报效果进行评价，若评价效果符合设定的阈值，则径流预报模型建立完成，该模型用于对未来时段径流进行预报，否则，重新调整人工神经网络的模型参数，重新建模直至模型完成建立。

s5具体可以为：

使用验证好的模型对下一时段的自变量y进行预报，同时对下一时段的趋势项trendf进行预测，其预测径流的公式为：

其中，qf为预测的下一时段径流量，imff(j)为第j层本征模态函数分量的预测值，resf为预测的残余量。

s403中，所述对人工神经网络模型训练期和检验期对自变量y的预报效果进行评价，，评价的指标包括纳什效率系数、相对误差和合格率；

所述纳什效率系数nash按照如下公式进行计算：

其中，qo为实测径流量，qf为预测径流量，qo为实测径流量均值，nash越接近1，预报越精准；

所述相对误差mare按照如下公式进行计算：

其中，n为qo的序列长度，其余参数含义同上式，mare越接近0，说明实测与预报值越接近，预报效果越精准，常认为mare＜20％时，效果较好；

所述合格率qr按照如下公式进行计算：

其中，n为合格预报次数，m为预报总次数，当qr＞80％时，认为预报效果较好。

s1中，所述滞后时段t的遥相关数据d，可以包括130项遥相关数据，其中，88项环流指数，26项海温数据，16项其他因子数据；滞后时段t主要考虑滞后1～12个月，共同组成一组130×12的气象因子矩阵x(d,t)。

s2中，所述n层本征模态函数项imf(j)和一层残余项res的关系式可以为：

其中，q'为连续多年的除趋势项的径流序列，imf(j)为第j层本征模态函数分量，n为本征模态函数的分量数。

s301中，所述相关性分析，通过皮尔逊相关系数判断，具体可以采用如下公式：

其中，r(x,y)为相关因子，cov(x,y)为x和y的协方差，var[x]为x的方差var[y]为y的方差，|r(x,y)|越接近1，x和y相关性越大。

具体实施例

本实施例中，选取雅砻江锦屏一级水库的入库月径流作为实施例，对本发明内容的效果进行验证，具体按照如下步骤进行实施：

步骤1，提取径流序列及响应的预报因子数据。

在本发明的实施例，气息数据为同期130项(88项环流指数，26项海温数据，16项其他因子数据)气候因子所组成的矩阵。本发明的案例中，为了更好地分析径流数据的同期变化，选用雅砻江锦屏一级1954年1月-2010年12月的入库月流量作为径流序列[q1,q2,...,q12]。

为了减少干扰，还原径流序列[q1,q2,...,q12]的稳态序列，先将其剔除趋势项trend生成新的无趋势变化的径流序列[q′1,q′2,...,q′12]＝[q1-trend1,q2-trend2,...,q12-trend12]作为自变量。

本发明考虑了气象因子影响的滞后性，需要将滞后期加入因变量范围，故在本发明的实施例中，滞后期选择1～12个月，与130个气象因子共同组成一组130×12的气象因子矩阵x(d,t)，如下式所示：

其中，t为自变量的预报时刻，x(d,t)为t时刻的备选预报气象因子矩阵，di,t-j为第i∈(1,130)个气象因子滞后期为j∈(1，12)的备选预报气象因子。

根据本发明的实施例，由于考虑气象因子的滞后期1～12个月，1954年1月至12月的气象因子的滞后期，仅能推求1955年1月的径流，故1954年全年的径流数据不参与预报模型建模。

步骤2，径流序列[q′1,q'2,...,q′12]的稳态化。

如图2所示，确定集合经验模态分解方法的基本参数：噪声方差(nstd＝0.2)、噪声组数(ne＝100)、迭代次数(maxiter＝500)，通过集合经验模态分解的方法将径流序列[q′1,q'2,...,q′12]分别分解为满足hilbert变换要求的n阶稳态序列的本征模态函数和1组残余项构成的1～12月的12组因变量集合y1＝[imf1(1),imf1(2),...,imf1(n),res1]，y2＝[imf2(1),imf2(2),...,imf2(n),res2]，...，y12＝[imf12(1),imf12(2),...,imf12(n),res12]，每组序列均为离散函数，为不同频率的线性或非线性的序列。

每一组关于第i月的imfi(1),imfi(2),...,imfi(n),resi的因变量集合yi，具体关系式如下：

其中，yi为第i月的因变量序列，imfi(j)为第i月第j层的本征模态序列，resi为第i月残余序列。

根据本发明的实施例，分解1955年～2010年的[q′1,q'2,…,q′12]自变量，由于篇幅所限，仅展现1955年～2010年7月q'7的分解结果的第一层imf7(1)，如图3所示。

步骤3，确定自变量y＝[y1,y2,...,y12]相关的预报因子fi,j(x,n+1)。

(3-1)对气象因子x(d,t)与12个月每个月的imfi(1),imfi(2),...,imfi(n),resi的n+1个数据序列进行相关分析，按照下式计算相关系数；

(3-2)设定α＝0.025的显著水平，对序列和气象因子的相关系数进行t检验，计算如下式，通过假设检验的因子被认为显著相关。将通过t检验的因子根据相关系数降序排序，选取前15个为初步筛选的备选预报因子；

其中，n为资料样本数，r为皮尔逊相关系数相关系数。

(3-3)根据逐步回归方法，挑选显著影响因子，确定每个月的imfi(1),imfi(2),...,imfi(n),resi的最终预报因子fi,j(x,n+1)，完成对自变量y＝[y1,y2,...,y12]的预报因子识别工作，其最终的选取的预报因子。

步骤4，分析12个月每个月的imfi(1),imfi(2),...,imfi(n),resi与之相对应的预报因子fi,j(x,n+1)的相关性，分析影响自变量集合y＝[y1,y2,...,y12]的的气象因子，挖掘径流序列包含的每个月imfi(1),imfi(2),...,imfi(n),resi所隐藏的物理意义，仅分析影响的因子，故不考虑滞后性。对全年分解出的6层函数数列imf(1),imf(2),...,imf(5),res与，以及不同月份全部本征模态序列和残余序列分别与气象因子的相关系数进行统计分析，发掘不同月份，不同分层的数据序列的气象机理。

步骤5，构建基于集合经验模态分解与遥相关数据的径流预报模型。利用筛选出来的预报因子fi,j(x,n+1)和imfi(1),imfi(2),...,imfi(n),resi预报下一时段的径流,包括以下子步骤：

(5-1)采用10折交叉检验，构建模型的训练、验证样本集，根据数据长度，确定模型训练期和检验期长度，并将输入因变量fi,j(x,n+1)为已识别的气象因子及确定的滞后期，将自变量imfi(1),imfi(2),…,imfi(n),resi分为训练样本集和检验期样本集两段；

(5-2)基于已获得的训练样本集，建立并训练人工神经网络模型，进而预报检验期的自变量y＝[y1,y2,...,y12]；

(5-3)使用纳什效率系数、相对误差、合格率来评估人工神经网络模型训练期和检验期对自变量y＝[y1,y2,...,y12]的预报效果。本发明实施例中，模型的预报结果如图4所示，评价结果如下表，可以看出，预报效果较好。

(5-4)使用验证好的模型对未来自变量y＝[y1,y2,...,y12]进行预报，同时对每个月下一时段的趋势项trendi,f进行预测，其预测的径流如下式：

其中qi,f为第i月预测的下一时段径流量，imfi,f(j)为第i月第j层本征模态函数分量的预测值，resi,f为预测的残余量。

通过采用本发明公开的上述技术方案，得到了如下有益的效果：本发明实施例提供的集合经验模态分解及遥相关的中长期径流预报方法，与现有技术相比，充分考虑了径流序列的非稳态性，使用集合经验模态分解的方法将径流转化为多组稳态序列，实现了水文序列稳态化，为常规的预报方法提供了最基础的数据条件。同时，针对常规自适应中长期径流预报在物理机理方面的不足，提出采用遥相关数据分析的方法，分析多组稳态序列的气象影响因子，能充分剥离不同频率序列的物理因子，挖掘径流序列隐藏的物理机制，建立的预报模型提高了中长期径流预报的准确性。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视本发明的保护范围。