CN104732092B - 一种基于聚类的水文降雨一致区分析方法 - Google Patents

Abstract

本发明为一种基于聚类的水文降雨一致区分析方法，公开了一种基于雨量站降雨特征的降雨区域相似分析方法，首先分割日降雨量时间序列，提取单场降雨时间序列；然后选择和提取单场降雨特征向量，并对单场降雨特征进行标准化，并对单场降雨特征向量进行PCA变换，基于PCA变换特征值的贡献率保留单场降雨的主要特征，实现单场降雨特征的降维，然后基于改进K‑means聚类方法对单场降雨进行聚类，得到单场降雨类型；最后构造雨量站的降雨特征多元时间序列，基于S_PCA度量降雨特征多元时间序列的相似程度，从而建立雨量站相似模型。

Description

一种基于聚类的水文降雨一致区分析方法

技术领域

本发明涉及一种能够进行降雨区域相似的方法，具体是一种基于雨量站降雨特征的降雨区域相似分析方法，属于水文水资源和数据挖掘交叉的领域。

背景技术

降雨是水资源的主要构成之一，水资源是社会和经济发展的基础。为了进行流域水资源科学管理和调度，需要全面了解流域水资源现状，特别是降雨信息。目前我国大部分地区雨量站分布稀疏，虽然有3400多个国家级水文站点，但布局不甚合理，大多部署在经济发达地区。随着水利信息化的不断推进，新增设了许多雨量站，但由于许多雨量站设站时间较晚，部分站点降雨资料只有短短几年。为了更准确地获取基础水文数据，科学管理水资源，需要解决以下问题：第一、缺失资料站点的数据补齐。处于同一降雨一致区域内的站点之间的降雨数据可以互相补充，如A站点有20年的降雨数据，B站点有30年的降雨数据，若它们在近20年内的降雨类型相似，则B站可以补齐A站前 10年的数据；第二、站点的合理规划。去除旧的冗余站点，新增必须的测站。地理位置相近的两个站点若降雨相似度极高，则可以考虑去除其中的一个冗余站点(若两站点相似度很高但地理位置较远，一般不考虑去除冗余站点的问题)；第三，对相似类型的降雨站点进行分类，作为一个重要的参考依据，以便于对降雨站点进行统一有效的管理。

目前对降雨站点进行相似性几乎没有相关研究，也没有其他方式为水资源管理提供这方面的技术支撑。主要是在于如何对降雨数据进行科学地标准化并精确聚类，提供具有指导性的建议存在较大困难。

发明内容

发明目的：本发明提供一种基于聚类的水文降雨一致区分析方法，为补齐雨量站缺失数据，以及水文站点的合理布设提供科学依据。

技术方案：为了实现上述发明目的，本发明的一种基于聚类的水文降雨一致区分析方法包括如下步骤：

(1)对雨量站监测的汛期的日雨量时间序列，根据单场降雨分割规则进行分割，

获得单场降雨时间序列；

(2)统计所述单场降雨时间序列的特征量，获得单场降雨n维特征向量；

(3)对所述单场降雨n维特征向量进行主成份变换，保留特征值贡献率超过阈值

的p个维度，并向p维子空间重新投影，获得单场降雨p维特征向量；

(4)基于改进K-means方法对单场降雨p维特征向量进行聚类；

(5)基于雨量站所有汛期的降雨类型特征向量建立基于监测时间的降雨特征序列；

(6)采用相似度量建立相似性模型。

同一场降雨有一定的持续时间，不同的降雨场次之间有一个多日无雨的时间段(无雨期)，因此日降雨时间序列的数据是不连续的。为了分析各雨量站监测到的单场降雨的特征，需要将雨量站监测到的单场降雨量时间序列提取出来。根据实际的降雨观察，同一场降雨过程中，有可能某一天的降雨量有可能变小，甚至雨量可以忽略，但是从降雨的本质上来看，这一场雨实际上并没有结束，甚至存在一种情况，有连续两天由于雨量值实在太小，以至于难以观测到，但是这样的无雨间隔不会太长。所述单场降雨分割规则是将雨量监测站无降雨记录的间隔天数大于两天的相邻降雨记录为不同的单场降雨。

本发明对降雨特征进行了精心的挑选。它将直接影响单场降雨聚类结果的合理性，进一步影响相似性比较的有效性。所述单场降雨n维特征向量为：

{p，d，a，dmax，dmin，p127，dp127，p50，dp50}，

其中，p表示单场降雨总雨量，d表示单场降雨持续天数，a表示日平均降雨量，dmax表示日降雨量最大值，dmin表示日降雨量最小值，dp127表示日降雨量小于1.27mm的雨量，dp127表示以及降雨天数，p50表示日降雨量大于50mm的总雨量，dp50表示降雨天数。

单场降雨总雨量p和与降雨天数d可以反映单场降雨总量和持续时间；日平均降雨量a能够反映单场降雨的平均强弱情况；日降雨量最大值pmax和日降雨量最小值pmin 能够反映单场降雨的极值情况；我国气象部门规定，暴雨为24h雨量超过50mm的降雨，所以选择50mm/d作为日降雨极大值的统计量；根据Munger指数，日降雨量小于1.27mm 可以判定干旱天气，所以选择1.27mm/d作为降雨的极小值的统计量。

由于每一个特征值在数值上大小不同，单位不同，造成了每一个特征分量对聚类结果的贡献不同，所以需要对特征向量进行标准化处理，以达到聚类算法对输入数据的要求。本发明对所述单场降雨n维特征向量进行z-score标准化处理。z-score标准化基于原始数据的均值(mean)和标准差(standard deviation)进行数据的标准化，标准化后数据的均值为0，标准差为1。

获取单场降雨n维特征向量后，这n个维度之间会存在相互的关联，因此存在冗余信息，同时较多的维度对聚类结果也会产生不利的影响。为了解决这一问题，所述步骤 (3)中利用PCA变换对单场降雨n维特征向量进行降维。从而保留主要信息，利用提取后的低维特征向量进行聚类，这样可以提高聚类的准确性。

假设所有雨量站的所有单场降雨特征向量构成单场降雨数据集X_m×n，X为行向量，即每一行为描述一个场次的单场降雨的特征向量，m为单场降雨特征个数，n个所有的单场降雨场次数。对X_m×n进行PCA变换，实际上就是对X_m×n的相关矩阵或者协方差矩阵进行奇异值分解。由于数据已经采用标准化，因此采用相关矩阵或者协方差矩阵得到的PCA结果相同。

本发明进行PCA变换的步骤为：

(31)计算X_m×n的相关矩阵M；

(32)进行PCA变换，得到特征向量矩阵C_m×m和特征值矩阵L_1×m，C_m×m的每列表示一个特征向量，是一个单位向量。保留特征值贡献率超过阈值，取前p个特征向量，这 p个特征向量互相正交，构成原始向量空间的一个子空间C_m×1…p；

(33)将单场降雨n维特征向量X_n×m投影到p维子空间，获得单场降雨p维度特征向量，投影的矩阵变换为X’＝X_n×m·C_m×1…p。X’中的元素在步骤(4)中作为被聚类的对象。

其中，阈值可选择80％或95％。

每个雨量站都监测到许多的单场降雨，若直接使用这些单场降雨进行相似性比较，处理过程繁琐，效率低。本发明的思路为：将单场降雨进行聚类后，每场次的单场降雨采用所属聚类标号标识，用于后继处理。所述步骤(4)改进K-means聚类方法按照顺序步骤是：

(41)设定聚类数k，以及最大的聚类数目k_min、最小的聚类数目k_max，并设定k 的初始值为k_min；

(42)若k大于k_max，则转到步骤(44)，否则转(43)；

(43)采用传统K-means聚类算法；

(44)计算类内平均距离和d_in与类间距离和d_out的比值定义r，

r＝d_in/d_out (1)，

其中，

d_i表示第i个类的类内平均距离：

N_i为类i中成员的个数，c_i为类i的中心，x_ij类表示类i中的第j个对象；

(45)k增加1；

(46)比较k在范围[k_min，k_max]内聚类时的r值，将r最小的聚类结果以及k作为最后输出，获得优化聚类数k。

传统的K-means算法是将含有n个数据点(实体)的集合X＝{x₁，x₂，...，x_n}划分为k个类簇C_j的问题(j＝1，2，...，k)。首先随机选取k个数据点作为k个类簇的初始簇中心，集合中每个数据点被划分到与其距离最近的簇中心所在的类簇中，形成初始的k个聚类分布。对已分配的每一个类簇计算新的簇中心，然后继续进行数据分配，这样迭代若干次后，若簇中心不再发生变化，说明数据对象已经全部分配到各自所在的类簇中。传统的K-means聚类方法需要首先指定聚类个数，这个对数据特征不了解的情况下，不容易实现。

本发明对K-means聚类算法进行改进，首先指定聚类数目范围内，然后算法递增聚类数目并聚类，计算聚类性能指标，最后以聚类性能指标最优的聚类作为最终聚类结果。聚类结束后，每个单场降雨采用其所属于的1…k之间的某个类标号进行标识，被聚到同一类中的单场降雨标识相同。

然后构造雨量站降雨特征多元时间序列。在对单场降雨序列进行聚类后，得到所有单场降雨的类型(即所属于类的标号)，雨量站监测的所有单场降雨都可以归到指定的类型。将雨量站监测到的降雨按年进行描述，采用单场降雨年特征向量描述雨量站年降雨特征。

最后基于雨量站单场降雨特征的多元时间序列进行相似性建模。所述步骤(6)基于监测时间的降雨特征序列建立S_PCA模型，计算两个雨量站之间的S_PCA距离：

其中，L、M分别为A、B两个雨量站基于监测时间的降雨特征序列，L^T、M^T分别为其转置矩阵，用于矩阵相乘后求矩阵的迹，d为维数，cos²θ_ij为A雨量站的第i个特征向量与B雨量站的第j个特征向量夹角的余弦平方值；S_PCA越大，表示两个雨量站的降雨特征越相似。

有益效果：本发明降雨区域相似分析方法在对单场降雨聚类前，对单场降雨的特征向量进行降维，利用改进的聚类方法自动确定最优聚类数目，避免人为指定带来的误差，提高聚类准确性，同时在雨量站相似建模时，采用雨量站监测到的逐年的降雨特征的多元时间序列描述雨量站的降雨特征，能够有效的体现雨量站监测的降雨的时间变化特性，隐含了本区域降雨的时间上的特征变化过程。在水文信息处理过程中，一般用其所测的雨量代表其附近的面状区域的降雨。本发明实际上是通过雨量站相似来得到降雨区域相似，从而提供降雨区域相似分析，为缺失资料站点的数据补齐、合理设置布设水文站点等提供科学依据。

附图说明

图1是本发明的流程示意图；

图2是本发明提取日降雨两时间序列的流程图；

图3是本发明改进K-means聚类流程图；

图4是某雨量站降雨特征多元时间序列。

具体实施方式

下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

如图1所示，一种基于聚类的水文降雨一致区分析方法，具体包括如下步骤：

步骤110：获得时段降雨时间序列

《实时雨水情数据库表结构与标识符SL323-2011》中日降雨量和时段降雨量混合存储在时段降雨量表中，因此需要从时段降雨提取出逐日降雨量信息。

步骤120：提取日降雨信息

时段降时间序列中包括两种逐日降雨信息，一是从前一天的06：00到当天的06：00，二是从前一天的08：00到当天的08：00。06：00是为了防汛会商的需求提前报的日雨量，一般常规采用08：00的日雨量作为前一天的日降雨。这里某个雨量站提取日降雨信息的规则如下：针对某日的所有的时段雨量，若存在08：00时刻的日降雨量，则采用该值作为前一天的日降雨量；若不存在08：00的日降雨量，但存在06：00的日降雨量，则采用 06：00作为前一天的日降雨量；若不存在08：00和06：00日雨量信息，则设定前一天的日雨量为0。

如图2提取日降雨时间序列流程图所示所示：

步骤121：从时段降雨时间序列中读入一个时段降雨值；

步骤122：判断读入的是否为08：00的时间值，若是则保留该值，转步骤121，若不是转步骤123；

步骤123：判断读入的是否为06：00的时间值，若是则保留该值，转步骤121，若不是转步骤121；

步骤124：判断时段降雨时间序列中，是否还有数据，若有，则转步骤121，否则，结束日降雨时间序列的提取，得到日降雨量时间序列。

时段降雨数据预处理后形成雨量站的日降雨量时间序列。

步骤130：划分单场降雨时间序列

同一场降雨有一定的持续时间，不同的降雨场次之间有一个多日无雨的时间段(无雨期)，因此日降雨时间序列的数据是不连续的。为了分析各雨量站监测到的单场降雨的特征，需要将雨量站监测到的单场降雨量时间序列提取出来。根据实际的降雨观察，同一场降雨过程中，有可能某一天的降雨量有可能变小，甚至雨量可以忽略，但是从降雨的本质上来看，这一场雨实际上并没有结束，甚至存在一种情况，有连续两天由于雨量值实在太小，以至于难以观测到，但是这样的无雨间隔不会太长。单场降雨分割规则是将雨量监测站无降雨记录的间隔天数大于两天的相邻降雨记录为不同的单场降雨。

举例来说，如果7月9日不降雨，但7月8日和7月10日两条降雨记录仍然视为同一场降雨序列中的信息。而如果7月9日、7月10日都没有降雨记录，7月8日和7 月11日的两条降雨记录，则应该分属于两场降雨，基于该的原则对降雨记录进行单场降雨进行划分。

步骤140：单场降雨特征提取

特征选择非常关键，它将直接影响单场降雨聚类结果的合理性，进一步影响相似性比较的有效性。本事实例提取单场降雨9维特征向量为：

{p，d，a，dmax，dmin，p127，dp127，p50，dp50}，

其中，p表示单场降雨总雨量，d表示单场降雨持续天数，a表示日平均降雨量，dmax表示日降雨量最大值，dmin表示日降雨量最小值，dp127表示日降雨量小于1.27mm的雨量，dp127表示以及降雨天数，p50表示日降雨量大于50mm的总雨量，dp50表示降雨天数。

单场降雨总雨量p和与降雨天数d可以反映单场降雨总量和持续时间；日平均降雨量a能够反映单场降雨的平均强弱情况；日降雨量最大值pmax和日降雨量最小值pmin 能够反映单场降雨的极值情况；我国气象部门规定，暴雨为24h雨量超过50mm的降雨，所以选择50mm/d作为日降雨极大值的统计量；根据Munger指数，日降雨量小于1.27mm 可以判定干旱天气，所以选择1.27mm/d作为降雨的极小值的统计量。

步骤150：单场降雨特征标准化

由于每一个特征值在数值上大小不同，单位不同，造成了每一个特征分量对聚类结果的贡献不同，所以需要对特征向量进行标准化处理，以达到聚类算法对输入数据的要求。本发明对单场降雨n维特征向量进行z-score标准化处理。z-score标准化基于原始数据的均值(mean)和标准差(standard deviation)进行数据的标准化，标准化后数据的均值为0，标准差为1。

假设属性A的均值为m，标准差为v，则属性A中原始值x使用z-score标准化到 x′的公式为x’＝(x-m)/v。

本发明中的数据标准化，是对所有雨量站所有的单场降雨特征向量中的特征量进行标准化，如对p的标准化，是所有雨量站所有的单场降雨中的p构成的数据序列的标准化。

步骤160：单场降雨特征PCA变换降维

获取单场降雨9维特征向量后，这9个维度之间会存在相互的关联，因此存在冗余信息，同时较多的维度对聚类结果也会产生不利的影响。为了解决这一问题，采用主成分变换法，利用PCA变换对单场降雨9维特征向量进行降维。从而保留主要信息，利用提取后的低维特征向量进行聚类，这样可以提高聚类的准确性。

假设所有雨量站的所有单场降雨特征向量构成单场降雨数据集X_m×n，X为行向量，即每一行为描述一个场次的单场降雨的特征向量，m为单场降雨特征个数，n个所有的单场降雨场次数。对X_m×n进行PCA变换，实际上就是对X_m×n的相关矩阵或者协方差矩阵进行奇异值分解。由于数据已经采用标准化，因此采用相关矩阵或者协方差矩阵得到的PCA结果相同。

本发明进行PCA变换的步骤为：

(161)计算X_m×n的相关矩阵M；

(162)进行PCA变换，得到特征向量矩阵C_m×m和特征值矩阵L_1×m，C_m×m的每列表示一个特征向量，是一个单位向量。保留特征值贡献率达到95％的阈值，取前p个特征向量，这p个特征向量互相正交，构成原始向量空间的一个子空间C_m×1…p；

(163)将单场降雨n维特征向量X_n×m投影到p维子空间，获得单场降雨p维度特征向量，投影的矩阵变换为X’＝X_n×m·C_m×1…p。X’中的元素在步骤(170)中作为被聚类的对象。

PCA变换的过程可利用如Matlab等软件实现，该步骤重点是在PCA分析后，将原始特征向量投影到新的p维子空间中，得到单场降雨的p维特征向量，后继的分析采用单场降雨的p维特征向量。

步骤170：基于改进K-means聚类划分单场降雨类型

每个雨量站都监测到许多的单场降雨，若直接使用这些单场降雨进行相似性比较，处理过程繁琐，效率低。如图3所示，针对单场降雨9维特征向量，现将单场降雨进行聚类后，每场次的单场降雨采用所属聚类标号标识，用于后继处理：

(171)设定聚类数k，以及最大的聚类数目k_min、最小的聚类数目k_max，并设定k 的初始值为k_min；

(172)若k大于k_max，则转到步骤(176)，否则转(173)；

(173)采用传统K-means聚类算法；

(174)计算类内平均距离和d_in与类间距离和d_out的比值定义r，

r＝d_in/d_out (1)，

其中，

d_i表示第i个类的类内平均距离：

N_i为类i中成员的个数，c_i为类i的中心，x_ij类表示类i中的第j个对象；

(175)k增加1，回到步骤(172)；

(176)比较k在范围[k_min，k_max]内聚类时的r值，将r最小的聚类结果以及k作为最后输出，获得优化聚类数k。

聚类结束后，每个单场降雨采用其所属于的1…k之间的某个类标号进行标识，被聚到同一类中的单场降雨标识相同。

步骤180：构造雨量站降雨特征多元时间序列

构造雨量站降雨特征多元时间序列。在对单场降雨序列进行聚类后，得到所有单场降雨的类型(即所属于类的标号)，雨量站监测的所有单场降雨都可以归到指定的类型。将雨量站监测到的降雨按年进行描述，采用单场降雨年特征向量描述雨量站年降雨特征。

若对所有的单场降雨聚类后，形成的最优聚类数目为k，某雨量站Y年监测到的单场降雨中，属于类1的单场降雨有m₁场，属于类2的单场降雨有m₂场，属于类i的单场降雨有m_i场，则雨量站Y年监测到的降雨的特征向量M为：

M＝{m₁，m₂，...，m_k}，

进一步将其表达为基于时间序列的降雨特征向量S，行为监测年份，维度数为聚类数目。

S＝{s_ti|监测数据的开始年＜＝t＜＝监测数据结束年，1＜＝i＜＝k}，

图4为某个雨量站的降雨特征向量。

步骤190：雨量站降雨特征相似模型

最后基于雨量站单场降雨特征的多元时间序列进行相似性建模。基于监测时间的降雨特征序列建立S_PCA模型，具体流程是：

(191)对雨量站降雨特征多元时间序列进行PCA变换，获取雨量站PCA变换的特征向量；

(192)计算A，B两个雨量站之间的S_PCA距离：

其中，L，M分别为A，B两个雨量站的监测时间降雨特征序列，d为维数。从几何上来看，θ_ij其实是A雨量站的第i个特征向量与B雨量站的第j个特征向量夹角的余弦值。S_PCA越大，表示两个雨量站的降雨特征越相似。

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

Claims (4)

1.一种基于聚类的水文降雨一致区分析方法，其特征在于包括如下步骤：

(1)对雨量站监测的汛期的日雨量时间序列，根据单场降雨分割规则进行分割，获得单场降雨时间序列；所述单场降雨分割规则是将雨量站无降雨记录的间隔天数大于两天的相邻降雨记录为不同的单场降雨；

(2)统计所述单场降雨时间序列的特征量，获得单场降雨n维特征向量；

(3)对所述单场降雨n维特征向量进行主成份变换，保留特征值贡献率超过阈值的p个维度，并向p维子空间重新投影，获得单场降雨p维特征向量；

(4)基于改进K-means方法对单场降雨p维特征向量进行聚类；

(5)基于雨量站所有汛期的降雨类型特征向量建立基于监测时间的降雨特征序列；

(6)采用相似度量建立相似性模型，基于监测时间的降雨特征序列建立S_PCA模型，

计算两个雨量站之间的S_PCA距离：

其中，L、M分别为A、B两个雨量站基于监测时间的降雨特征序列，L^T、M^T分别为其转置矩阵，用于矩阵相乘后求矩阵的迹，d为维数，cos²θ_ij为A雨量站的第i个特征向量与B雨量站的第j个特征向量夹角的余弦平方值；

所述步骤(4)改进K-means聚类方法的步骤是：

(41)设定聚类数k，以及最大的聚类数目k_max、最小的聚类数目k_min，并设定k的初始值为k_min；

(42)若k大于k_max，则转到步骤(46)，否则转(43)；

(43)采用传统K-means聚类算法；

(44)计算类内平均距离和d_in与类间距离和d_out的比值定义r，

r＝d_in/d_out (1)，

其中，

d_i表示第i个类的类内平均距离：

N_i为类i中成员的个数，c_i为类i的中心，x_ij表示类i中的第j个对象；

(45)k增加1，回到步骤(42)；

(46)比较k在范围[k_min，k_max]内聚类时的r值，将r最小的聚类结果以及k作为最后输出，获得优化聚类数k。

2.根据权利要求1所述的方法，其特征在于，所述单场降雨n维特征向量为：

{p，d，a，dmax，dmin，p127，dp127，p50，dp50}，

其中，p表示单场降雨总雨量，d表示单场降雨持续天数，a表示日平均降雨量，dmax表示日降雨量最大值，dmin表示日降雨量最小值，p127表示日降雨量小于1.27mm的雨量，dp127表示日降雨量小于1.27mm的降雨天数，p50表示日降雨量大于50mm的总雨量，dp50表示日降雨量大于50mm的降雨天数。

3.根据权利要求1或2所述的方法，其特征在于：对所述单场降雨n维特征向量进行z-score标准化处理。

4.根据权利要求1所述的方法，其特征在于：所述步骤(3)中利用PCA变换对单场降雨n维特征向量进行降维。