一种基于多模型组合的中长期径流集合预报方法

摘要

本发明公开了一种基于多模型组合的中长期径流集合预报方法，涉及水文预测技术领域；该方法采用多种机器学习算法构建中长期径流预报模型并将其用作弱学习器，以此为基础提出了一种基于多模型组合的集成模型构建方法。同时，通过参数扰动搜寻等效预报构建预报集合，进行集合预报。相较于目前常用确定性预报方法，针对其中存在的部分缺陷进行了改进，提升了中长期预报的精度和泛化能力。同时，通过概率预报定量化描述了预报的不确定性，提升了预报的准确性和参考价值。

权利要求

1.一种基于多模型组合的中长期径流集合预报方法，其特征在于，包括以下步骤：S1，气候系统指数‑径流相关性分析：选取多项气候系统指数历史数据与待预报流域历史径流数据，采用相关性分析法确定径流量与气候系统指数之间的相关性系数；选取相关性系数绝对值最大的前20项气候系统指数作为初选因子，形成初选因子矩阵；

S2，得到初选因子矩阵后，采用主成分分析进一步对步骤S1中得到的初选因子矩阵进行降维处理，提取其中影响径流过程的关键影响因子，形成关键影响因子矩阵；

S3，对步骤S2得到的关键影响因子矩阵进行归一化处理；将归一化处理后的关键影响因子矩阵与对应年份的预报流域历史径流数据结合，构建影响因子‑径流数据集；

S4，训练弱学习器：基于SVR、KNN和BP神经网络算法，构建基于智能算法的样本预报模型，并使用步骤S3中构建的影响因子‑径流数据集对弱学习器进行训练；

S5，构建特征预报‑径流数据集：基于步骤S4中训练的弱学习器构建集成模型，对三个弱学习器的预报数据进行归纳整理，形成特征预报‑径流数据集，用于所述集成模型的训练和验证；

S6，构建集成模型：采用多元回归算法，率定回归权重并建立回归模型，以此构建集成模型；

S7，构建预报集合：基于集合预报对预报数量的要求，需在步骤S6得到的集成模型的基础上，利用参数扰动算法构建等效参数集，通过等效参数集扩展预报数量，构建预报集合；

步骤S6具体包括：

为每个弱学习器分别设定一个权重，采用PSO算法搜索各弱学习器的回归权重和回归偏置的最优值，以集成模型在训练集上的均方根相对误差RMSRE和整体波动性EZ的线性组合作为PSO算法的适应度评估函数，其中均方根相对误差 整体波动性 其中，σi表示第i个预报样本的相对误差，yi为模型的预测值；

寻得最优回归权重后，通过权重线性组合的形式构建回归模型，其回归公式为其中αi为第i个弱学习器的回归权重，ri为第i个弱学习器的预报结果，β为模型的回归偏置，n为弱学习器的个数。

2.根据权利要求1所述的基于多模型组合的中长期径流集合预报方法，其特征在于，还包括以下步骤：

S8，计算集合预报输出：集合预报通过整理、提炼预报集合上的预报信息，整合所有预报计算集合预报输出；通过统计所有集合成员给出的预报结果，以预报集合上的集合均值作为数值预报输出，以预报集合上的置信区间作为概率预报输出；

S9，预报结果展示：以真实径流数据为基准，定量评价集合预报在测试集上的数值预报结果，并输出数值预报结果；以图表的形式展示预报集合在测试集上的置信区间。

3.根据权利要求1所述的基于多模型组合的中长期径流集合预报方法，其特征在于，步骤S1中计算径流量y和各气候系统指数x的相关性系数的公式为：

4.根据权利要求1所述的基于多模型组合的中长期径流集合预报方法，其特征在于，步骤S3中对关键因子采用公式(2)进行归一化处理：式中， 为关键影响因子的归一化值，xij为关键影响因子的原始值，μj为第j列的列均值， σj为第j列的方差， n为矩阵的行数。

5.根据权利要求1所述的基于多模型组合的中长期径流集合预报方法，其特征在于，步骤S4具体包括：

S41，训练PSO‑SVR模型：采用粒子群优化算法确定SVR模型中误差惩罚参数C、核参数σ和不敏感损失系数ε三个参数的最优值；

S42，构建KNN模型：采用手动试算方法确定KNN方法中K的取值，并使用训练集作为样本库，构建KNN模型；

S43，训练神经网络模型：采用循环迭代穷举法确定BP神经模型中的网络结构参数，设定模型的隐藏层数不大于4，隐藏层的总节点数不大于使用kolmogov定理计算得到的最大隐层节点总数；依据最大隐藏层数设置三个循环遍历所有网络结构，每个循环依据网络层数和网络隐层节点总数设定循环结构，迭代中按照设定的网络结构训练神经网络，以平均相对误差评价模型在训练集上的拟合效果，选择平均相对误差最优的网络结构参数作为神经网络模型的网络结构参数。

6.根据权利要求1所述的基于多模型组合的中长期径流集合预报方法，其特征在于，步骤S5具体包括：统计每个弱学习器的预报结果，将其按照年份由远及近的顺序排成一列，而后整合三个弱学习器的预报结果，构建特征预报矩阵，矩阵中每行包括三个特征，且同一行的特征均为同一时段不同模型的预报结果；

选定特征预报矩阵中S年的预报数据构建特征预报‑径流数据集，将数据集中前N年的真实径流量与对应时间的特征预报矩阵组合作为训练集，数据集中后M年的真实径流量与对应时间的特征预报矩阵组合作为测试集，其中M，N的限制条件为S＝N+M，N＞2*M，S、N、M均为正整数。

7.根据权利要求1所述的基于多模型组合的中长期径流集合预报方法，其特征在于，步骤S7具体包括：

S71，构建等效参数集：为避免预报集合中预报同质化，在使用PSO对回归权重进行寻优的过程中，改变适应度评估函数中RMSRE和EZ的权重，以此寻找集成模型的等效参数；

S72，构建预报集合：在步骤S71的基础上，以MAPE值和EZ作为评价指标，分别对使用等效参数构建的集成模型进行定量评价；由定量评价结果，依照设定的流程选择预报集合的成员。

8.根据权利要求7所述的基于多模型组合的中长期径流集合预报方法，其特征在于，步骤S72具体包括：

S721，以预报结果的MAPE从小到大的顺序对集成方案进行排序，若两个方案的MAPE相等则按照Ez从小到大的方式进行排序；

S722，由步骤S721中的排序结果，取首尾集成方案的MAPE值，计算两者MAPE的差值；

S723，以两者MAPE的差值作为模型间准确度偏差，若偏差值大于10％，则去除当前排位最后的集成模型，而后返回步骤S722重新计算首尾集成模型MAPE的差值；若偏差值小于

10％，则认为模型间精确度偏差满足条件；

S724，统计当前集成方案的数量，采用当前的集成方案集合作为预报集合。

说明书

一种基于多模型组合的中长期径流集合预报方法

技术领域

[0001] 本发明涉及水文预报技术领域，尤其涉及一种基于多模型组合的中长期径流集合预报方法。

背景技术

[0002] 受气候、气象、下垫面、人类活动等诸多因素的综合影响，流域中长期径流预报过程具有一定的时空不确定性。传统的中长期径流预报方法如物理成因法、数理统计方法、回
归分析方法和现代新兴的预报方法如模糊分析法、灰色系统方法、神经网络等方法常以气
象因子作为备选因子，从中为特定流域选取适当的预报因子，以此学习预报因子与流域径
流之间关系。传统的径流预报方法受预报模式和数据的影响，往往难以解决单一确定性预
报因预报模式的高度非线性和径流过程的混沌性造成的预报技巧严重受限的问题，其预报
在可靠性和准确定方面仍有一定的不足。

发明内容

[0003] 本发明的目的在于提供一种基于多模型组合的中长期径流集合预报方法，从而解决现有技术中存在的前述问题。

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

[0005] 一种基于多模型组合的中长期径流集合预报方法，包括以下步骤：

[0006] S1，气候系统指数‑径流相关性分析：选取多项气候系统指数历史数据与待预报流域历史径流数据，采用相关性分析法计算径流量与气候系统指数之间的相关性系数；选取
相关性系数绝对值最大的前20项作为初选因子，形成初选因子矩阵；

[0007] S2，提取关键影响因子：进一步对步骤S1中得到的初选因子矩阵进行降维处理，提取其中影响径流过程的关键影响因子，生成关键影响因子矩阵；

[0008] S3，构建影响因子‑径流数据集：对步骤S2得到的关键影响因子矩阵进行归一化处理，将归一化后的关键影响因子矩阵与选定的待预报流域历史径流数据结合，构建影响因
子‑径流数据集；

[0009] S4，训练样本预报模型：采用SVR、KNN和BP神经网络等机器学习算法，构建基于智能算法的弱学习器，并使用步骤S3中构建的影响因子‑径流数据集训练弱学习器；

[0010] S5，构建特征预报‑径流数据集：对三个样本预报模型的预报进行归纳整理，构建特征预报‑径流数据集，用于所述集成模型的训练和验证；

[0011] S6，构建集成模型：基于多元回归方法，设定回归权重和偏置建立回归模型集成多个弱学习器，以此构建集成模型；

[0012] S7，构建基于等效参数集的预报集合：基于集合预报对预报数量的要求，在步骤S6得到的集成模型的基础上，采用等效参数构建等效集成方案，通过等效参数集扩展预报数
量，以此构建预报集合。

[0013] 优选地，还包括以下步骤：

[0014] S8，计算集合预报输出：集合预报通过整理、提炼预报集合上的的预报信息，整合所有预报计算集合预报输出；具体的，以预报集合上的集合均值作为数值预报输出；以预报
集合上预报的95％置信区间作为概率预报输出。

[0015] S9，预报结果展示：以真实径流数据为基准，定量评价集合预报在测试集上的数值预报结果，并输出数值预报结果；以图表的形式展示预报集合在测试集上的置信区间。

[0016] 优选地，步骤S1中计算径流量y和各气候系统指数x的相关性系数的公式为：

[0017] 优选地，步骤S3中对关键影响因子采用公式(2)进行归一化处理：

[0018]

[0019] 式中， 为归一化的关键影响因子，xij为该元素的原始值，μj为第j列的列均值，σj为第j列的方差， n为矩阵行数。

[0020] 优选地，步骤S4具体包括：

[0021] S41，训练PSO‑SVR模型：采用粒子群优化算法搜索SVR模型中误差惩罚参数C、核参数σ和不敏感损失系数ε三个参数的最优值；

[0022] S42，构建KNN模型：采用手动试算方法确定KNN方法中K的取值，并使用训练集作为样本库，构建KNN模型；其中K取值一般为不大于5的整数。

[0023] S43，训练神经网络模型：采用循环迭代穷举的方法确定BP神经模型中的网络结构参数，包括网络层数和各隐藏层节点数，具体过程为：

[0024] 采用循环迭代穷举法确定BP神经模型中的网络结构参数，设定模型的隐藏层数不大于4，隐藏层的总节点数不大于使用kolmogov定理计算得到的最大隐层节点总数；依据最
大隐藏层数设置三个循环遍历所有网络结构，每个循环依据网络层数和网络隐层节点总数
设定循环结构，迭代中按照设定的网络结构训练神经网络，以平均相对误差评价模型在测
试集上的拟合效果，选择平均相对误差最优的网络结构参数作为神经网络模型的网络结构
参数。

[0025] 更优选地，步骤S41中具体包括：使用C、σ、ε的值作为粒子群中粒子的位置的坐标值，根据设定的粒子群规模随机初始化粒子群，计算每个粒子在当前位置的适应度，通过比
较适应度确定粒子的个体最佳位置和粒子群的全局最佳位置，并更新每个粒子的速度和位
置，迭代至满足结束条件或者满足最优输出结果，得到位置最优的里的位置坐标值即为SVR
模型参数C、σ、ε的最佳值。

[0026] 优选地，步骤S5具体包括：统计每个弱学习器的预报结果，将其按照年份由远及近的顺序排成一列，而后整合三个弱学习器的预报结果，构建特征预报矩阵，矩阵中每行包括
三个特征，且同一行的特征均为同一时段不同模型的预报结果；

[0027] 选定特征预报矩阵中S年的预报数据构建特征预报‑径流数据集，将数据集中前N年的真实径流量与对应时间的特征预报矩阵组合作为训练集，数据集中后M年的真实径流
量与对应时间的特征预报矩阵组合作为测试集，其中M，N 的限制条件为S＝N+M，N>2*M，S、
N、M均为正整数。

[0028] 优选地，步骤S6具体包括：

[0029] 分别为每个弱学习器设定一个权重，采用PSO算法对各弱学习器的权重和回归偏置进行寻优，以集成模型在训练集上的均方根相对误差RMSRE和整体波动性EZ的线性组合作
为PSO算法的适应度评估函数，其中均方根相对误差 整体波
动性

[0030] 寻得最优回归权重后，通过模型、权重线性组合的形式构建集成模型，其回归公式为 其中αi为第i个弱学习器的回归权重，ri为第i个弱学习器的预报结
果，β为模型的回归偏置，n为弱学习器的个数。

[0031] 优选地，步骤S7具体包括：

[0032] S71，构建等效参数集：为避免集合中预报同质化，在使用PSO对回归权重进行寻优的过程中，改变适应度评估函数中RMSRE和EZ的权重，以此寻找集成模型的等效参数；

[0033] S72，构建预报集合：在步骤S71的基础上，以MAPE值和EZ作为评价指标，对基于等效参数构建的集成模型进行定量评价；由定量评价结果，依照设定的流程选择集合的预报
成员。

[0034] 更优选地，步骤S72具体包括：

[0035] S721，以MAPE从小到大的顺序对集成方案进行排序，若两个方案的MAPE 相等则按照Ez从小到大的方式进行排序；

[0036] S722，由步骤S721中的排序结果，取首尾集成方案的MAPE值，计算两者 MAPE的差值；

[0037] S723，以两者MAPE的差值作为模型间准确度偏差，若偏差值大于10％，则去除当前排位最后的集成模型，而后返回步骤S722重新计算首尾集成模型 MAPE的差值；若偏差值小
于10％，则认为模型间精确度偏差满足条件；

[0038] S724，统计当前集成模型的数量，使用当前的集成方案集合作为预报集合。

[0039] 本发明的有益效果是：

[0040] 本发明公开了一种基于多模型组合的中长期径流集合预报方法，该方法采用多种机器学习算法构建中长期径流预报模型并将其用作弱学习器，以此为基础提出了一种基于
多模型组合的集成模型构建方法。同时，通过参数扰动搜寻等效预报构建预报集合，进行集
合预报。相较于目前常用确定性预报方法，针对其中存在的部分缺陷进行了改进，提升了中
长期预报的精度和泛化能力。同时，通过概率预报定量化描述了预报的不确定性，提升了预
报的准确性和参考价值。

附图说明

[0041] 图1是本发明提供的一种基于多模型组合的中长期径流集合预报方法的技术路线；

[0042] 图2是是实施例1中基于PSO‑SVR的样本预报模型构建流程图；

[0043] 图3是具体实施方式中丹江口水库8月月径流集合预报模型在全部数据集上的概率预报结果图。

具体实施方式

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

[0045] 实施例

[0046] 本实施例提供一种基于多模型组合的中长期径流集合预报方法，包括以下步骤：

[0047] S1，气候系统指数‑径流相关性分析：选取多项气候系统指数历史数据与待预报流域历史径流数据，采用相关性分析法计算径流量与气候系统指数之间的相关性系数；选取
相关性系数绝对值最大的前20项作为初选因子，形成初选因子矩阵；

[0048] S2，提取关键影响因子：进一步对步骤S1中得到的初选因子矩阵进行降维处理，提取其中影响径流过程的关键影响因子，生成关键影响因子矩阵；

[0049] S3，构建影响因子‑径流数据集：对步骤S2得到的关键影响因子矩阵进行归一化处理，将归一化后的关键影响因子矩阵与选定的待预报流域历史径流数据结合，构建影响因
子‑径流数据集；

[0050] S4，训练样本预报模型：采用SVR、KNN和BP神经网络等机器学习算法，构建基于智能算法的弱学习器，并使用步骤S3中构建的影响因子‑径流数据集训练弱学习器；

[0051] S5，构建特征预报‑径流数据集：对三个样本预报模型的预报进行归纳整理，构建特征预报‑径流数据集，用于所述集成模型的训练和验证；

[0052] S6，构建集成模型：基于多元回归方法，设定回归权重和偏置建立回归模型集成多个弱学习器，以此构建集成模型；

[0053] S7，构建基于等效参数集的预报集合：基于集合预报对预报数量的要求，在步骤S6得到的集成模型的基础上，采用等效参数构建等效集成方案，通过等效参数集扩展预报数
量，以此构建预报集合；

[0054] S8，计算集合预报输出：集合预报通过整理、提炼预报集合上的的预报信息，整合所有预报计算集合预报输出；具体的，以预报集合上的集合均值作为数值预报输出；以预报
集合上预报的95％置信区间作为概率预报输出。

[0055] S9，预报结果展示：以真实径流数据为基准，定量评价集合预报在测试集上的数值预报结果，并输出数值预报结果；以图表的形式展示预报集合在测试集上的置信区间。

[0056] 本实施例中的步骤S1中计算径流量y和各气候系统指数x的相关性系数的公式为：

[0057] 本实施例中步骤S3中对关键影响因子采用公式(2)进行归一化处理：

[0058]

[0059] (2)式中， 为归一化的关键影响因子，xij为该元素的原始值，μj为第j 列的列均值， σj为第j列的方差， n为矩阵行数。

[0060] 本实施例中，步骤S4具体包括：

[0061] S41，训练PSO‑SVR模型：采用粒子群优化算法搜索SVR模型中误差惩罚参数C、核参数σ和不敏感损失系数ε三个参数的最优值；

[0062] 具体为：使用C、σ、ε的值作为粒子群中粒子的位置的坐标值，根据设定的粒子群规模随机初始化粒子群，计算每个粒子在当前位置的适应度，通过比较适应度确定粒子的个
体最佳位置和粒子群的全局最佳位置，并更新每个粒子的速度和位置，迭代至满足结束条
件或者满足最优输出结果，得到位置最优的里的位置坐标值即为SVR模型参数C、σ、ε的最佳
值。

[0063] S42，构建KNN模型：采用手动试算方法确定KNN方法中K的取值，并使用训练集作为样本库，构建KNN模型；其中K取值一般为不大于5的整数。

[0064] S43，训练神经网络模型：采用循环迭代穷举的方法确定BP神经模型中的网络结构参数，包括网络层数和各隐藏层节点数，具体过程为：

[0065] 采用循环迭代穷举法确定BP神经模型中的网络结构参数，设定模型的隐藏层数不大于4，隐藏层的总节点数不大于使用kolmogov定理计算得到的最大隐层节点总数；依据最
大隐藏层数设置三个循环遍历所有网络结构，每个循环依据网络层数和网络隐层节点总数
设定循环结构，迭代中按照设定的网络结构训练神经网络，以平均相对误差评价模型在测
试集上的拟合效果，选择平均相对误差最优的网络结构参数作为神经网络模型的网络结构
参数。

[0066] 本实施例中步骤S5具体包括：统计每个弱学习器的预报结果，将其按照年份由远及近的顺序排成一列，而后整合三个弱学习器的预报结果，构建特征预报矩阵，矩阵中每行
包括三个特征，且同一行的特征均为同一时段不同模型的预报结果；

[0067] 选定特征预报矩阵中S年的预报数据构建特征预报‑径流数据集，将数据集中前N年的真实径流量与对应时间的特征预报矩阵组合作为训练集，数据集中后M年的真实径流
量与对应时间的特征预报矩阵组合作为测试集，其中M，N 的限制条件为S＝N+M，N>2*M，S、
N、M均为正整数。

[0068] 本实施例中步骤S6中具体包括：分别为每个弱学习器设定一个权重，采用 PSO算法对各弱学习器的权重和回归偏置进行寻优，以集成模型在训练集上的均方根相对误差
RMSRE和整体波动性EZ的线性组合作为PSO算法的适应度评估函数，其中均方根相对误差
整体波动性

[0069] 寻得最优回归权重后，通过模型、权重线性组合的形式构建集成模型，其回归公式为 其中αi为第i个弱学习器的回归权重，ri为第i个弱学习器的预报结
果，β为模型的回归偏置，n为弱学习器的个数。

[0070] 本实施例中步骤S7具体包括：

[0071] S71，构建等效参数集：为避免集合中预报同质化，在使用PSO对回归权重进行寻优的过程中，改变适应度评估函数中RMSRE和EZ的权重，以此寻找集成模型的等效参数；

[0072] S72，构建预报集合：在步骤S71的基础上，以MAPE值和EZ作为评价指标，对基于等效参数构建的集成模型进行定量评价；由定量评价结果，依照设定的流程选择集合的预报
成员。

[0073] 步骤S72具体包括：

[0074] S721，以MAPE从小到大的顺序对集成方案进行排序，若两个方案的MAPE 相等则按照EZ从小到大的方式进行排序；

[0075] S722，由步骤S721中的排序结果，取首尾集成方案的MAPE值，计算两者 MAPE的差值；

[0076] S723，以两者MAPE的差值作为模型间准确度偏差，若偏差值大于10％，则去除当前排位最后的集成模型，而后返回步骤S722重新计算首尾集成模型 MAPE的差值；若偏差值小
于10％，则认为模型间精确度偏差满足条件；

[0077] S724，统计当前集成模型的数量，使用当前的集成方案集合作为预报集合。

[0078] 具体实施例

[0079] 本实施例提出了一种基于多模型组合的中长期径流集合预报方法，如图1所示，首先分析130项气候系统指数与历史径流的相关性关系并从中提取初选因子矩阵，而后通过
主成分分析法对初选因子矩阵进行降维，得到关键影响因子矩阵。再对其进行零均值归一
化处理，然后结合历史径流数据构建样本模型数据集。将样本模型数据集分为两类，一类为
训练集，一类为测试集，并使用训练集分别训练基于SVR、KNN和BP神经网络的样本模型，三
个样本模型的训练过程如下：

[0080] 训练SVR模型的流程如图2所示，采用粒子群优化算法确定SVR模型中 C、σ、ε三个参数的值，其大致步骤如下：设定粒子群规模，并为每个粒子随机设定初始位置，而后计算
粒子的个体最佳位置和粒子群的全局最佳适应度，更新每个粒子的速度和位置，不断迭代
直至达到预期的预测精度或者达到迭代次数上限，得到位置最好的粒子的坐标值即为SVR
模型参数C、σ、ε的最优解；将得到的C、σ、ε三个参数代入SVR模型中，使用训练集训练模型。

[0081] 训练KNN模型时采用手动试算法确定KNN模型中的关键参数K值，其大致步骤如下：依据K的取值范围[2,5]，设定初始K值。而后基于设定的K值，以训练集为样本库构建KNN模
型，然后使用预设的评价指标对模型进行评价。迭代所有K的取值，选择评价指标最高的KNN
模型作为样本模型，该样本模型的K值作为模型的最佳参数。

[0082] 训练BP神经网络模型时迭代穷举法确定BP神经网络模型中的隐藏层数目 Nhidden和隐层节点数 其大致步骤如下：根据kolmogov定理计算网络最大的隐层节点总数，
而后使用嵌套迭代划定范围内的所有网络结构。依据可能的网络层数设置循环数，对
进行迭代，由确定网络结构参数和预设的超参数训练神经网络模型，然后使用预设的评价
指标对模型进行评价。迭代所有的网络结构，选择评价指标最高的神经网络模型作为样本
模型，该样本模型的网络结构参数作为模型的最佳参数。

[0083] 样本模型训练完成后，计算样本预报模型在测试集上的预测精度，分析误差。由三个样本模型给出的预报构建集成模型数据集，将集成模型数据集分为两类，一类为训练集，
一类为测试集，并使用训练集训练基于多元回归的集成模型。训练集成模型时，采用粒子群
优化算法对样本模型的权重进行寻优，采用寻优得到的最优权重构建集成模型。在粒子群
优化算法寻优的过程中，通过改变适应度评估函数中评价指标的权重，设定不同的适应度
评估函数，寻找集成模型的等效参数。而后通过等效参数训练集成方案，并由集成方案构建
预报集合。最后，综合预报集合中所有方案给出的预报，计算数值预报结果和概率预报结
果，对未来径流进行预报。

[0084] 具体实施例

[0085] 本实施例中，对汉江流域丹江口水库8月入库径流进行预报，按照本发明提供的方法进行实施，该方法包括以下步骤：

[0086] 1)气候系统指数‑径流相关性分析：选取多项气候系统指数历史数据与待预报流域历史径流数据，采用相关性分析法确定径流量与气候系统指数的相关性系数，利用公式
计算径流量y和各气候系统指数x
的相关性系数，选取相关性系数绝对值最大的前20项作为初选因子。

[0087] 本具体实施方式中从国家气候中心官方网站获取历年百项气候系统指数集 (即88项大气环流指数、26项海温指数和16项其他指数)作为历史多项气候系统指数数据，将该
历史多项气候指数数据与丹江口水库多年历史8月径流数据进行相关性分析，对相关性分
析结果按相关性系数绝对值从大到小进行排序。月预报选取前半年的气象因子进行相关行
分析，每月130项因子，即选取780项气候系统指数因子参与相关计算，由于数量过多不全部
展示，取排序结果的前20项气候系统指数因子(初选因子)进行展示，如下表所示：

[0088]

[0089]

[0090] 2)提取预报因子：本实施例中使用的主成分分析法对步骤1)得到的初选因子矩阵进行降维，设定选取因子的总贡献度设为95％，从初选因子矩阵中提取关键影响因子矩阵
如下表所示：

[0091]

[0092]

[0093] 3)构建影响因子‑径流数据集：采用零均值标准化方法对步骤2)得到的关键影响因子矩阵进行归一化处理，利用公式对关键影响因子矩阵中的元素进行处理：
式中， 为归一化的关键影响因子，xij为该元素的原始值，μj为当前列的列均值，
σj为当前列的方差，  n为矩阵的行数；将归一
化处理后得到的关键影响因子矩阵与选定的待预报流域历史径流数据结合，构建影响因
子‑径流数据集。本实例采用丹江口水库近38年 (1981年～2018年)的入库径流数据作为样
本径流，结合关键影响因子，构建影响因子‑径流数据集并对数据集进行划分，即选用前22
年(1981年～2002年)的历史径流数据与对应的标准化关键影响因子构建训练集，选用后16
年(2003年～2018年) 的历史径流数据与对应的标准化关键影响因子构建测试集。经归一
化处理后得到的训练集如下表所示：

[0094]

[0095]

[0096] 经归一化处理后得到的测试集如下表所示：

[0097]

[0098] 4)训练样本预报模型：基于SVR、KNN和BP神经网络算法，构建基于智能算法的样本模型，并利用步骤3)中构建的影响因子‑径流数据集对弱学习器进行训练。

[0099] 4‑1)采用图2中的流程训练PSO‑SVR模型：本实施例中用于优化SVR模型参数C、σ、ε的粒子群优化算法具体步骤包括：

[0100] 4‑1‑1)确定C、σ、ε三个参数的取值范围,即粒子位置和速度的取值范围。本实施例中这三个值的取值范围均采用现有方法中的相关数据，其中粒子位置的最大值设定为
(100,200,100)，最小值设定为(0.01,0.01,0.01)；粒子速度的调整幅度最大值设定为(10,
1,10),最小值设定为(‑10,‑1,‑10)。

[0101] 4‑1‑2)初始化粒子群，即设定粒子群的规模、迭代次数、位置和速度等：本实施例中粒子群规模设定为300，粒子群迭代次数设定为1000，粒子的位置和速度在步骤4‑1‑1)中
设定的取值范围中取随机值作为粒子的初始位置和初始速度。

[0102] 4‑1‑3)确定适应度评估函数，并使用适应度函数计算每个粒子的适应度值：本实施例选取确定性系数作为适应度评估函数，使用如下所示的公式作为评价函数计算各粒子
的适应度值： 式中，DC为确定性系数， yc(i)为预测值，yo(i)
为实测值， 为多年实测值的平均值，n为训练集样本年数。

[0103] 4‑1‑4)确定每个粒子的个体最佳位置：将每个粒子当前位置的适应度与其历史最佳位置localbest(即局部最优解)的适应度值作比较，确定当前的最佳位置 localbest。

[0104] 4‑1‑5)确定整个粒子群的全局最佳位置。将每个粒子当前最佳位置的适应度与整个群体的当前最佳位置globalbest(即全局最优解)的适应度值作比较，确定当前的最佳位
置globalbest。

[0105] 4‑1‑6)更新粒子的速度和位置：本实施例中根据公式：vi+1＝wvi+ c1r1(localbest‑xi)+c2r2(globalbest‑xi)更新粒子速度，根据公式： xi+1＝xi+vi+1更新粒子位置；式中，i代表
迭代次数，xi表示第i次迭代时粒子的位置，vi表示第i次迭代时粒子的速度，r1,r2表示为
(0,1)之间的两个随机数，c1c2表示的是增速因子，其值一般取为2，w为动态权重因子，其取
值范围为[0.4,0.9], 其根据式 动态更新权重w，其中wini
和wend分别为动态权重因子的起始值和终止值，其值分别为0.9和0.4,G为PSO算法的迭代次
数，g 为当前的迭代次数。

[0106] 4‑1‑7)判断算法是否满足结束条件，不满足则转到4‑1‑3)步骤；满足则输出最优结果，此时的全局优解即是三个参数的最优解；

[0107] 经过计算，得到SVR的C、σ、ε三个参数的全局最优解为：100、95.81、0.01；将参数带入SVR模型中进行训练，得到预报模型，并利用预报模型在测试集上进行径流预报，得到的
预报结果如下表所示：

[0108]

[0109] 4‑2)构建KNN模型：本实施例构建KNN样本预报模型的具体步骤包括：

[0110] 4‑2‑1)基于K值构建样本模型：由于K取正整数且取值范围一般为[2,5],本实施例中通过人工设定的方式遍历所有K值，并基于选定的K值，结合训练集构建 KNN模型。

[0111] 4‑2‑2)使用评价指标对模型进行评价：本实施例采用以下公式对模型进行评价：式中yi为模型的预报值， 为实测的径流值。

[0112] 4‑2‑3)选择最优K值：经过4‑2‑3)中对所有KNN模型的评估后，选取MAPE 值最小的模型作为最佳模型，采用该模型的K值作为最优K值。

[0113] 4‑2‑4)基于最佳K值构建KNN模型。基于选择的最佳K值，结合训练集构建 KNN模型。

[0114] 经过评估后，得到最优的K值为5，使用训练集构建KNN模型样本库，并利用KNN模型在测试集上进行径流预报，得到的预报结果如下表所示：

[0115]

[0116] 4‑3)训练神经网络模型：采用循环迭代对比的方法确定BP神经模型中的网络结构参数(网络层数和各层节点数)，设定模型的隐藏层层数不大于4，隐藏层的总节点数不大于
使用kolmogov定理计算最大隐层节点数；基于网络层数设定循环个数，各层节点数作为循
环对象，迭代所有网络结构，以训练结果的MAPE 作为评价指标，选择MAPE最大的网络结构
参数作为神经网络模型的网络结构参数。本实施例构建BP神经网络样本预报模型的具体步
骤包括：

[0117] 4‑3‑1)使用kolmogov定理计算最大隐层节点总数。根据kolmogov定理使用以下公式计算最大隐层节点总数： 式中，M为网络节点总数，Mi为第i
层中的节点数目，N为网络层数，n为输入节点数,且当N≥2时，单层神经元数目满足以下公
式：M≤2n+1。本实施例中模型的输入节点数为5，故由kolmogov定理计算得到网络中最大的
隐层节点数为11。

[0118] 4‑3‑2)基于循环结构迭代所有网络结构。本实施例中设计的最大隐藏层数为 3，故设置3个循环分别对隐层数为1、2、3的网络结构进行迭代。在循环中由设定网络层数设定
循环的嵌套层数，即单隐层网络设定一层循环，2隐层网络设定双层嵌套循环，以此类推。在
循环的最里层，依据循环设定的网络结构构建神经网络结构框架，使用影响因子‑径流数据
集的训练集训练结构测试模型。

[0119] 4‑3‑3)使用评价指标评价所有结构测试模型：本实施例采用以下公式对模型进行评价： 式中yi为模型的预报值， 为实测的径流值。

[0120] 4‑3‑4)对比分析确定网络结构：经过4‑3‑3)中对所有神经网络模型的评估后，选取MAPE值最小的模型作为最佳模型，采用该模型的网络结构作为模型的预设结构训练样本
预报模型。

[0121] 经过迭代评估后，得到神经网络的网络结构参数为：9‑38‑1，隐藏层层间使用Dropout函数随机消除50％节点，迭代600次，采用训练集训练神经网络模型，并利用得到模
型在测试集上进行径流预报，得到的预报结果如下表所示：

[0122]

[0123] 5)构建特征预报‑径流数据集：对三个样本预报模型的预报进行归纳整理，构建特征预报‑径流数据集，用于集成模型的训练和验证。统计每个弱学习器给出的预报结果，将
其按照年份由远及近的顺序排成一列，而后整合三个弱学习器的预报结果构建特征预报矩
阵，本实施例中的特征预报矩阵中的每行有三个元素，且所有元素均为同一时段的预报结
果。本实例采用丹江口水库近38年(1981 年～2018年)的8月历史入库径流数据作为样本径
流，结合样本模型的预报结果，构建样本模型数据集并对数据集进行划分，即选用前22年
(1981年～2002年)的历史真实径流数据与对应时间的特征预报矩阵组合构建训练集，选用
后16年(2003 年～2018年)的历史真实径流数据与对应时间的特征预报矩阵组合构建测试
集。

[0124] 经整理后得到的训练集如下表所示：

[0125]

[0126] 经整理后得到的测试集如下表所示：

[0127]

[0128] 6)构建集成模型：本实施例中用于集成模型权重的粒子群优化算法具体步骤包括：

[0129] 6‑1)确定回归权重的取值范围,即粒子位置和速度的取值范围。本实施例中这三个值的取值范围均采用现有方法中的相关数据，其中粒子位置的最大值设定为(1,1,1)，最
小值设定为(0,0,0)；粒子速度的调整幅度最大值设定为(0.1,0.1,0.1), 最小值设定为(‑
0.1,‑0.1,‑0.1)。

[0130] 6‑2)初始化粒子群，即设定粒子群的规模、迭代次数、位置和速度等。本实施例中粒子群规模设定为300，粒子群迭代次数设定为1000，粒子的位置和速度在步骤6‑1)中设定
的取值范围中取随机值作为粒子的初始位置和初始速度。

[0131] 6‑3)确定适应度评估函数，并使用适应度函数计算每个粒子的适应度值。本实施例选取确定性系数作为适应度评估函数，使用如下所示的公式作为适应度评价函数计算各
粒子的适应度值： 式中， fitness为适应度系数，
为RMSRE的权重， 为EZ的权重，且 RMSRE均方根相对误差，
Ez为整体波动性，

[0132] 6‑4)确定每个粒子的个体最佳位置。将每个粒子当前位置的适应度与其历史最佳位置localbest(即局部最优解)的适应度值作比较，确定当前的最佳位置 localbest。

[0133] 6‑5)确定整个粒子群的全局最佳位置。将每个粒子当前最佳位置的适应度与整个群体的当前最佳位置globalbest(即全局最优解)的适应度值作比较，确定当前的最佳位置
globalbest。

[0134] 6‑6)更新粒子的速度和位置。本实施例根据式vi+1＝wvi+c1r1(localbest‑ xi)+c2r2(globalbest‑xi)更新粒子速度，根据式xi+1＝xi+vi+1更新粒子位置；式中，i代表迭代次数，xi
表示第i次迭代时粒子的位置，vi表示第i次迭代时粒子的速度，r1,r2表示为(0,1)之间的两
个随机数，c1c2表示的是增速因子，其值一般取为2，w为动态权重因子，其取值范围为[0.4,
0.9],其根据式 动态更新权重w，其中wini和wend分别为动态
权重因子的起始值和终止值，其值分别为0.9和0.4，G为PSO算法的迭代次数，g为当前的迭
代次数。

[0135] 6‑7)判断算法是否满足结束条件，不满足则转到6‑3)步骤；满足则输出最优结果，此时的全局优解即是三个权重的最优解；本实施例中结束条件设定为算法的迭代次数达到
上限1000次。

[0136] 7)构建基于等效参数集的预报集合：由于集合预报对预报数量的要求，需在步骤6)得到的集成模型的基础上，以等效参数构建集成方案，通过等效参数集扩展预报数量，构
建预报集合。

[0137] 7‑1)构建等效参数集：为避免集合中预报同质化，在使用粒子群优化算法对回归权重进行寻优的过程中，改变适应度评估函数中RMSRE和EZ的权重，以此寻找集成模型的等效
(t)
参数。本实施例中定义以下公式计算RMSRE的权重：ω ＝(ωini‑ωend)(Gk‑g)/Gk+ωend；式
中，ωini＝0.8为初始权重，ωend＝0.2为终止权重，Gk为设定的目标等效参数的数目，g为当
寻优的参数方案的标号。使用以下公式计算RZ的权重： 式中，
为RMSRE的权重， 为EZ的权重。通过不断调整调整寻优目标函数，改变PSO寻得的各模型回
归权重和偏置，构建等效参数集。

[0138] 7‑2)构建预报集合：在步骤7‑1)的基础上，以MAPE和EZ作为评价指标，对基于等效参数构建的集成模型进行定量评价。由定量评价结果，由设定的流程选择集合的预报成员。
预报成员的选择流程如下：

[0139] 7‑2‑1)以MAPE从大到小的顺序对集成方案进行排序，若两个方案的MAPE 相等则按照Ez从小到大的方式进行排序。

[0140] 7‑2‑2)由排序结果，取首尾集成方案的MAPE值，计算两者MAPE的差值。

[0141] 7‑2‑3)以两者MAPE的差值作为模型间准确度偏差，若偏差值大于50％，则去除当前排位最后的集成模型，而后返回步骤7‑2‑2)重新计算首尾集成模型 MAPE的差值；若偏差
值小于10％，则认为模型间精确度偏差满足条件；

[0142] 7‑2‑4)统计当前集成模型的数量，使用当前的集成方案集合作为预报集合；

[0143] 由步骤7计算统计得到集合预报成员(等效参数组)如下表所示：

[0144]

[0145]

[0146] 8)计算集合预报产品：

[0147] 8‑1)数值预报产品的计算。本实施例中采用集合均值作为数值预报产品，通过以下公式进行计算： 式中，Vpred为集合预报给出的数值预报结果，Vi为预
报集合中第i个预报，N为预报集合中预报成员的数目。

[0148] 8‑2)概率预报产品的计算。本实施例中采用95％置信区间作为概率预报产品，具体的步骤如下：

[0149] 8‑2‑1)计算样本均值。由步骤8‑1)计算方法可得到预报集合上的样本均值。

[0150] 8‑2‑2)计算标准误差。通过以下公式计算集合上预报的标准误差： 式中，ste为标准误差，std为标准差， N为预报集合中的预报成员数
目。

[0151] 8‑2‑3)计算上、下界：通过如下公式计算置信区间的上界：upLimit＝Z*ste+ mean；式中，Z为置信区间的设定值，ste为标准误差，mean为均值。通过如下公式计算置信区
间的下界：downLimit＝‑Z*ste+mean；式中，Z为置信区间的设定值，ste为标准误差，mean为
均值。

[0152] 9)预报结果展示：以测试集中真实径流数据对比集合预报给出的数值预报产品，进行定量评价，并输入数值预报结果；以图表的形式展示集合预报给出的概率预报产品。

[0153] 将训练得到的预报成员集合带入集合预报模型，并使用测试集进行测试，在全部数据集上的数值预报结果折线图如图3所示，在测试集上得到的数值预报结果如下表所示：

[0154]

[0155]

[0156] 将训练得到的预报成员集合带入集合预报模型，并使用测试集进行测试，在全部数据集上的概率预报结果折线图如图3所示。

[0157] 通过采用本发明公开的上述技术方案，得到了如下有益的效果：

[0158] 本发明公开了一种基于多模型组合的中长期径流集合预报方法，该方法采用多种智能方法构建中长期径流预报模型并将其用作弱学习器，以此为基础提出了一种基于多模
型组合的集成模型构建方法。同时，通过参数扰动搜寻等效预报构建预报集合，进行集合预
报。相较于目前常用确定性预报方法，针对其中存在的部分缺陷进行了改进，提升了中长期
预报的精度和泛化能力。同时，通过概率预报定量化描述了预报的不确定性，提升了预报的
准确性和参考价值。

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