基于粒子群优化的K-means径流丰枯年份预报方法

摘要

本发明公开了一种基于粒子群优化的K‑means径流丰枯年份预报方法，涉及水文预报技术领域。该方法，首先选定多年径流量数据作为历史样本，并分为训练集和验证集，通过粒子群优化算法对训练集中的样本进行寻优，得到多组预报因子组合，将各组预报因子组合分别用于K‑means聚类算法中进行训练集样本的聚类，对得到的聚类结果分别计算准确率，选取准确率最高的结果以及对应的预报因子组合，经过检测集检验合格后，利用该预报因子组合来进行预报年份径流丰枯年的预报，极大的提高了径流丰枯年的预报精度、合理性以及实用性，可以作为中长期径流丰枯年份预报的一种有效的方法。

说明书

技术领域

本发明涉及水文预报技术领域，尤其涉及一种基于粒子群优化的K-means径流丰枯年份预报方法。

背景技术

在水文学中，径流是指降雨及冰雪融水或者在浇地的时候在重力作用下沿地表或地下流动的水流。受天气系统和下垫面因素等的综合影响,河川径流呈现出复杂的年内和年际变化,根据径流的这种变化，可以分为丰水年和枯水年，由于径流变化的波动具有强烈的随机性，所以，对径流时序的丰枯状态变化进行预报,是认识水文动力系统演化特征的重要内容,对水资源的开发和利用等具有重要意义。

目前,径流丰枯划分方法主要有灰色分类评价法、距平百分率及模糊聚类法等。其中，聚类是一个将数据集中在某些方面相似的数据成员进行分类组织的过程。而K-means聚类是最著名的划分聚类算法，是一种基于样本间相似性度量的间接聚类方法，属于非监督学习方法，由于简洁和效率使得他成为所有聚类算法中最广泛使用的。此算法以k为参数，把n个对象分为k个簇，以使簇内具有较高的相似度，而且簇间的相似度较低。相似度的计算根据一个簇中对象的平均值(被看作簇的重心)来进行。此算法首先随机选择k个对象，每个对象代表一个聚类的质心。对于其余的每一个对象，根据该对象与各聚类质心之间的距离，把它分配到与之最相似的聚类中。然后，计算每个聚类的新质心。重复上述过程，直到准则函数收敛。

K-means的基本思想是通过预先设定的K值及每个类别的初始质心对相似的数据点进行划分，并通过划分后的均值迭代优化获得最优的聚类结果。基本流程图如图1所示，该方法主要包括以下步骤：

1)在数据集范围内，任意选择K个对象作为初始聚类中心；

2)遍历样本数据集合，对于每个样本，计算并比较其到K个聚类中心的距离，并按照距离最近原则将其分配到与其最相似的簇中；

3)计算所有簇中样本的均值作为新的聚类中心；

4)判断聚类中心是否发生变化，若是则跳到步骤2继续进行迭代，否则算法执行结束，返回聚类结果；

实验结果表明，上述聚类方法容易陷入局部最优，从而无法得到最优的分类结果，进而影响对径流丰枯年份预报的准确度。

发明内容

本发明的目的在于提供一种基于粒子群优化的K-means径流丰枯年份预报方法，从而解决现有技术中存在的前述问题。

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

一种基于粒子群优化的K-means径流丰枯年份预报方法，包括以下步骤：

S1，选取河川S年的径流量数据作为历史样本，将历史样本中前N年的径流量数据及对应年份的环流指数数据作为训练集，后M年的径流量数据及对应年份的环流指数数据作为检验集，其中，S＝N+M，N＞M，S、N、M均为正整数，并统计历史样本的实际丰枯年份；

S2，基于S1得到的训练集，采用粒子群优化算法进行多次迭代，得到K-means聚类算法的多个预报因子组合，所述预报因子组合包括5-6个粒子；

S3，对于每个预报因子组合，均利用该预报因子组合，根据K-means聚类算法，确定训练集的丰枯聚类结果，其中，K为2，包括丰水聚类和枯水聚类；

S4，对于S3中得到的每个丰枯聚类结果，计算其准确率；

S5，选取准确率最高的丰枯聚类结果以及对应的预报因子组合；

S6，根据S5确定的预报因子组合，根据K-means聚类算法，确定检验集的丰枯聚类结果；

S7，依据实际丰枯年份，对检验集的丰枯聚类结果进行误差分析，若误差不在设定范围内，则调整粒子群优化算法中的执行次数，并返回至S2；若误差在设定范围内，则将预报年份作为样本加入到检验集中，重复S6，得到预报年份的丰枯结果。

优选地，S2包括以下步骤：

S201，确定粒子群优化算法的如下初始条件：将S1得到的训练集样本的所有气候指数做为初始种群规模，将每个气候指数作为一个粒子，设定列数的取值范围；

S202，选定适应度函数，计算初始化种群中各粒子适应度函数值；

S203，在粒子群优化算法过程中对初始化种群各粒子的适应度值依次进行评估；

S204，执行次数均达到预定次数后，则满足终止条件，评估结束；对评估后剩下的粒子进行适应度值的比较，选出适应度值最大的粒子；

S205，确定每个粒子的个体最佳位置；

S206，确定整个粒子群的全局最佳位置；

S207，按照如下公式更新粒子的速度和位置：

v_i+1＝wv_i+c₁r₁(p_best-x_i)+c₂r₂(g_best-x_i)

x_i+1＝x_i+v_i+1

式中：i表示迭代代数，x_i表示第i次迭代时粒子所在位置，v_i表示第i次迭代时粒子的速度，r₁,r₂为(0,1)之间的两个随机数，c₁,c₂表示的是增速因子，它们的取值均大于0，一般都取为2，w为权重因子，取值范围是(0,1)。

优选地，S3包括以下步骤：

S301，遍历训练集，对于每个样本，根据预报因子组合计算并比较其到K个聚类中心的距离，并按照距离最近原则将其分配到与其最相似的簇中；

S302，计算所有簇中样本的均值作为新的聚类中心；

S303，判断聚类中心是否发生变化，若是则跳到S301继续进行迭代，否则算法执行结束，返回聚类结果。

优选地，S4中，按照如下公式计算其准确率：

其中，fitness为目标函数，n_f为丰水聚类年份中实际丰水年份的个数，n_k为枯水聚类年份中实际枯水年份的个数，N_f丰水聚类中的年份个数，N_k枯水聚类中的年份个数，为实际丰水年份个数占丰水聚类年份个数的比例，为实际枯水年份个数占枯水聚类年份个数的比例。

本发明的有益效果是：本发明实施例提供的基于粒子群优化的K-means径流丰枯年份预报方法，首先选定多年径流量数据作为历史样本，并分为训练集和验证集，通过粒子群优化算法对训练集中的样本进行寻优，得到多组预报因子组合，将各组预报因子组合分别用于K-means聚类算法中进行训练集样本的聚类，对得到的聚类结果分别计算准确率，选取准确率最高的结果以及对应的预报因子组合，经过检测集检验合格后，利用该预报因子组合来进行预报年份径流丰枯年的预报，极大的提高了径流丰枯年的预报精度、合理性以及实用性，可以作为中长期径流丰枯年份预报的一种有效的方法。

附图说明

图1是K-means聚类算法流程示意图；

图2是本发明提供的预报方法流程示意图；

图3是丹江口水库1981-2016年历年每月130项气候指数集部分截图，其中，第一列为丹江口水库1981-2016年历年年平均径流值，气候指数为12个月数据，共计130*12＝1560列。

具体实施方式

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

径流预报属于水文预报范畴，是应用水文学的一个重要组成部分，它是建立在掌握客观水文规律的基础上，预见未来径流变化的一门应用科学技术，是水资源调度、水利防汛和抗旱科学实施的前提。径流预报按预见期可分短期径流预报和中长期径流预报，划分的标准一般以流域汇流时间为界，凡预报的预见期小于或等于流域汇流时间的称为短期预报，预报的预见期大于流域汇流时间的称为中长期预报。

受气候、气象、下垫面、人类活动等诸多因素的综合影响，流域中长期径流预报过程具有一定的时空不确定性。从趋势性、周期性、随机性、区域性等方面分析径流的时空演变规律是水文学的重要研究内容之一，也是中长期径流预报的基础。从诸多水文-气象因子中为特定流域选取合适的预报因子集合，并构建预报因子集合与流域径流之间的关系，是中长期径流预报的难点。

本发明实施例中，通过粒子群优化算法得到最优的预报因子组合，并将其用于K-means聚类算法中，得到最优的聚类结果。提高了径流丰枯年的预报精度、合理性以及实用性，所以，本发明实施例提供的方法可以作为中长期径流丰枯年份预报的一种有效的方法。

如图2所示，本发明实施例提供了一种基于粒子群优化的K-means径流丰枯年份预报方法，包括以下步骤：

S1，选取河川S年的径流量数据作为历史样本，将历史样本中前N年的径流量数据及对应年份的环流指数数据作为训练集，后M年的径流量数据及对应年份的环流指数数据作为检验集，其中，S＝N+M，N＞M，S、N、M均为正整数，并统计历史样本的实际丰枯年份；

S2，基于S1得到的训练集，采用粒子群优化算法进行多次迭代，得到K-means聚类算法的多个预报因子组合，所述预报因子组合包括5-6个粒子；

S3，对于每个预报因子组合，均利用该预报因子组合，根据K-means聚类算法，确定训练集的丰枯聚类结果，其中，K为2，包括丰水聚类和枯水聚类；

S4，对于S3中得到的每个丰枯聚类结果，计算其准确率；

S5，选取准确率最高的丰枯聚类结果以及对应的预报因子组合；

S6，根据S5确定的预报因子组合，根据K-means聚类算法，确定检验集的丰枯聚类结果；

S7，依据实际丰枯年份，对检验集的丰枯聚类结果进行误差分析，若误差不在设定范围内，则调整粒子群优化算法中的执行次数，并返回至S2；若误差在设定范围内，则将预报年份作为样本加入到检验集中，重复S6，得到预报年份的丰枯结果。

其中，S2可以包括以下步骤：

S201，确定粒子群优化算法的如下初始条件：将S1得到的训练集样本的所有气候指数做为初始种群规模，将每个气候指数作为一个粒子，设定列数的取值范围；

S202，选定适应度函数，计算初始化种群中各粒子适应度函数值；

S203，在粒子群优化算法过程中对初始化种群各粒子的适应度值依次进行评估；

S204，执行次数均达到预定次数后，则满足终止条件，评估结束；对评估后剩下的粒子进行适应度值的比较，选出适应度值最大的粒子；

S205，确定每个粒子的个体最佳位置；

S206，确定整个粒子群的全局最佳位置；

S207，按照如下公式更新粒子的速度和位置：

v_i+1＝wv_i+c₁r₁(p_best-x_i)+c₂r₂(g_best-x_i)

x_i+1＝x_i+v_i+1

式中：i表示迭代代数，x_i表示第i次迭代时粒子所在位置，v_i表示第i次迭代时粒子的速度，r₁,r₂为(0,1)之间的两个随机数，c₁,c₂表示的是增速因子，它们的取值均大于0，一般都取为2，w为权重因子，取值范围是(0,1)。

其中，粒子群优化(Particle Swarm Optimization,PSO)算法是Kennedy和Eberhart受人工生命研究结果的启发、通过模拟鸟群觅食过程中的迁徙和群聚行为而提出的一种基于群体智能的全局随机搜索算法。PSO初始化为一群随机粒子(随机解)，然后通过迭代找到最优解。在每一次迭代中，粒子通过跟踪两个极值来更新自己；第一个就是粒子本身所找到的最优解，这个解称为个体极值；另一个极值是整个种群目前找到的最优解，这个极值是全局极值。另外也可以不用整个种群而只是用其中一部分作为粒子的邻居，那么在所有邻居中的极值就是局部极值。

S3可以包括以下步骤：

S301，遍历训练集，对于每个样本，根据预报因子组合计算并比较其到K个聚类中心的距离，并按照距离最近原则将其分配到与其最相似的簇中；

S302，计算所有簇中样本的均值作为新的聚类中心；

S303，判断聚类中心是否发生变化，若是则跳到S301继续进行迭代，否则算法执行结束，返回聚类结果。

S4中，可以按照如下公式计算其准确率：

其中，fitness为目标函数，n_f为丰水聚类年份中实际丰水年份的个数，n_k为枯水聚类年份中实际枯水年份的个数，N_f丰水聚类中的年份个数，N_k枯水聚类中的年份个数，为实际丰水年份个数占丰水聚类年份个数的比例，为实际枯水年份个数占枯水聚类年份个数的比例。

具体实施例

本发明实施例提供了一种基于粒子群优化的K-means丹江口水库径流丰枯年份预报方法，采用如下步骤进行实施：

1)选定历史资料：将历史的多项环流指数数据和丹江口水库历史径流数据作为K-means算法的输入样本；本实施例从国家气候中心官方网站下载36年(1981-2016年)每月130项气候指数集(如图3所示)，并全部作为K-means算法的输入样本；

2)选定丹江口水库S年的径流量数据作为历史样本，将历史样本中前N年的径流量数据及对应年份的环流指数作为训练集，后M年的径流量数据及对应年份的环流指数作为检验集，S＝N+M，N＞M，S、N、M均为正整数；本实施例采用丹江口水库近36年(1981-2016年)年径流量数据(如表1所示)作为历史样本，对该样本进行分类，即将历史样本前24年的年径流量数据作为训练集，将历史样本后12年的年径流量数据作为检验集；依据《水文情报预报规范》划分丰枯水年，采用距平百分率(36年(1981-2016年)的年平均径流值为1087.5m³/s，大于该值的为丰年，小于该值的为枯年)统计丰枯年份及丰枯年份的数量，结果如表1所示，从表中可以看出，丰年包括15年，枯年包括21年。

表1 36年(1981-2016年)年径流量数据(m³/s)及丰枯年份

3)利用步骤2)得到的训练集对K-means算法进行训练，采用粒子群优化算法确定K-means算法中的输入样本中的预报因子，使得分类结果达到最优：将训练集的每月的每个气候指数指数作为一个粒子，分别计算每个粒子的适应度值，对粒子适应度值进行评估以后得到满足条件的适应度值最大的粒子；具体步骤包括：

3-1)确定粒子群优化算法的初始条件：步骤2)得到的训练集的所有气候指数代表初始种群规模，将每个气候指数作为一个粒子，设定预报因子的取值范围；

3-2)选定适应度函数，计算得到的初始化种群中各粒子适应度函数值；

3-3)在粒子群优化算法过程中对初始化种群各粒子的适应度值依次进行评估；

3-4)执行次数均达到预定次数后，则满足终止条件，评估结束；对评估后剩下的粒子进行适应度值的比较，选出适应度值最大的粒子；

3-5)确定每个粒子的个体最佳位置；

3-6)确定整个粒子群的全局最佳位置；

3-7)更新粒子的速度和位置：

v_i+1＝wv_i+c₁r₁(p_best-x_i)+c₂r₂(g_best-x_i)

x_i+1＝x_i+v_i+1

式中：i表示迭代代数，x_i表示第i次迭代时粒子所在位置，v_i表示第i次迭代时粒子的速度，r₁,r₂为(0,1)之间的两个随机数，c₁,c₂表示的是增速因子，它们的取值均大于0，一般都取为2，w为权重因子，取值范围是(0,1)；

4)确定K个聚类中心，计算各粒子到聚类中心的距离；具体步骤如下：

4-1)遍历样本数据集合，对于每个样本，计算并比较其到K个聚类中心的距离，并按照距离最近原则将其分配到与其最相似的簇中；本实施例选定K值为2，即对应径流丰枯类别，统计每类中丰枯年份的个数，并计算各粒子到聚类中心的欧几里得距离，在进行距离计算时，利用得到的预报因子组合进行计算，比如预报因子组合如果包括5个维度，则计算5个维度的距离，并进行比较；

4-2)计算所有簇中样本的均值作为新的聚类中心；

4-3)判断聚类中心是否发生变化，若是则跳到步骤4-1)继续进行迭代，否则算法执行结束，返回聚类结果；

5)经过对丹江口水库历年每月气候指数集进行粒子群优化，以丰枯年份对应在K-means对应类别的准确率为目标函数计算最优的预报因子组合，同时得出丰枯年的预测结果，目标函数为：

其中fitness为目标函数，n_f为丰水聚类年份中实际丰水年份的个数，n_k为枯水聚类年份中实际枯水年份的个数，N_f丰水聚类中的年份个数，N_k枯水聚类中的年份个数，为实际丰水年份个数占丰水聚类年份个数的比例，为实际枯水年份个数占枯水聚类年份个数的比例。

经过计算，确定最优的预报因子组合为预报年份前1年9月NINO W区海表温度距平指数、预报年份前1年5月北半球极涡中心纬向位置指数、预报年份前1年7月西太平洋暖池面积指数、预报年份前1年11月冷空气次数、预报年份前1年12月南海副高脊线位置指数。根据最优预报因子得到的聚类结果的预测准确率为85.3％。

6)将步骤2)中选定的检验集与步骤5)中得到的初步预测结果进行比较，分析误差，若误差不在设定范围内，则调整粒子群优化算法中的执行次数，重新返回步骤3)；若误差在设定范围内，则输出最终预测结果。

7)根据确定的最优的预报因子组合，对预报年份的径流丰枯结果进行预测。

本实施例提供的方法克服了目前常用方法中线性分析选取预报因子方法的不恰当性，以及K-means算法容易陷入局部最优的缺点，所提出的方法编写简单、参数设置少、全局寻优能力强，可以有效避免计算量大、耗时长、需要的训练样本个数多等问题。本发明旨在提高预报精确度和预报效率，是一种高效的中长期径流丰枯年份预报方法，可以应用在实际的中长期径流预报中。

通过采用本发明公开的上述技术方案，得到了如下有益的效果：本发明实施例提供的基于粒子群优化的K-means径流丰枯年份预报方法，首先选定多年径流量数据作为历史样本，并分为训练集和验证集，通过粒子群优化算法对训练集中的样本进行寻优，得到多组预报因子组合，将各组预报因子组合分别用于K-means聚类算法中进行训练集样本的聚类，对得到的聚类结果分别计算准确率，选取准确率最高的结果以及对应的预报因子组合，经过检测集检验合格后，利用该预报因子组合来进行预报年份径流丰枯年的预报，极大的提高了径流丰枯年的预报精度、合理性以及实用性，可以作为中长期径流丰枯年份预报的一种有效的方法。

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