基于单场降雨类型的降雨站点相似性评价方法

摘要

本发明公开了一种基于单场降雨类型的降雨站点相似性评价方法，属于水文水资源与数据挖掘技术交叉领域。首先提取降雨时间序列中的雨量项，经过数据预处理，选择和提取单场降雨特征，并将这些特征进行数据标准化处理；采用基于DBI指数的K-means自动聚类方法进行聚类，最后根据聚类结果计算降雨类型直方图进行相似性分析。本发明方法具有计算方法简单，但能包含降雨时间序列更多更细粒度信息，真正反映降雨站点相似性的优点。

说明书

技术领域

本发明涉及一种降雨站点相似性评价方法，尤其涉及一种基于单场降雨类型的降雨 站点相似性评价方法，，属于水文水资源与数据挖掘技术交叉领域。

背景技术

时间序列是一种在科学研究、商业应用中普遍存在的数据形式，如灾变预测，降雨 预测，金融数据和电力负载数据预测等。在水文领域中，对时间序列的研究由来已久， 通过研究水文时间序列的时空规律，对洪水暴雨等自然灾害进行预测，为流域水利建设 及社会建设提供决策依据。降雨是水资源的主要来源之一，水资源是社会经济发展的重 要支撑和保障条件，水利、农业以及生态系统都受其影响。因此，研究降雨时间序列对 社会经济可持续发展和农业水利事业布局具有重要的意义。

目前我国大部分地区雨量站分布稀疏，许多雨量站设站时间较晚，部分站点降雨资 料只有短短几年。为了更准确地认识水文现象和解决水资源问题，需要对缺失资料的站 点数据进行补充，也需要设立新的观测站点，去除旧的冗余站点。因此有必要对降雨站 点进行相似性研究，为水资源管理提供新的技术支撑。

以往对降雨时间序列的研究主要包括降雨预测、极值分布、周期分析和降水场的时 空分布等。例如，陈晓宏等将模式识别的聚类分析理论方法引入降雨空间特征分类研究， 进行降雨特征空间分布模式识别，预测无雨量站点的降雨模式。以年降雨量的变换数据 为基础，结合高程、气温和蒸发量作为特征参数，分析不同聚类数下的分类有效性，将 雨量站点类型分为六类，再经过模型验证，预测无记录的站点归类，为降雨资料的插补 延长提供了一种新的途径；张国建等利用最大年降水量与最小年降水量比值来分析降水 量的年际变化特征；苏布达等利用weibull分布分析了长江流域降水极值时间序列的分 布特征；王兆礼等采用Morlet小波函数对东江流域在不同时间尺度下的降雨周期和突 变点进行了研究；刘德地等将云模型引入到降雨量时空分布特性的研究中；欧春平等研 究了基于信息熵的流域水文要素时空变异；巴金福等根据黄河上游逐月降水资料分析了 汛期降雨量的变化趋势和空间分布特征及其变化。

降雨序列实质是一个不连续的时间序列，大多数此前的研究，都是基于日月或年降 雨量单纯累积的宏观统计分析，没有从更细粒度分析单场降雨的统计特征。而单场降雨 序列恰恰包含了更丰富的降雨信息，能反映降雨序列的本质。例如单独对两个站点一年 的降雨进行统计，同样降雨量总和为1000mm的降雨，可能是8场单场降雨为125mm， 也可能是20场单场降雨为50mm，甚至还可能是若干场不同降雨量，不同持续天数的 降雨。若单纯利用年降雨量的累计而不考虑单场降雨，则无法区分一年内两个站点降雨 类型的差异，更无法区分降雨持续的时间等信息；若考虑单场降雨的情况，则可以区分 出站点更多的降雨信息，除了单场降雨和，单场降雨日均值，大于或小于某个临界值的 降雨量，降雨天数等更多降雨信息。

发明内容

本发明所要解决的技术问题在于克服现有技术的不足，提供一种基于单场降雨类型 的降雨站点相似性评价方法，针对降雨时间序列不连续的特点，对单场降雨的各种特征 进行统计分析，通过研究基于单场降雨的降雨序列相似性，继而研究区域降雨的相似性。

本发明的基于单场降雨类型的降雨站点相似性评价方法，包括以下步骤：

步骤A、从降雨站点一段时期雨量历史数据中提取各单场降雨的以下统计特征量：单场 降雨雨量和、单场降雨天数、单场降雨日平均降雨量、单场降雨日降雨量最大值、单场 降雨日降雨量最小值、单场降雨日降雨量小于1.27mm的雨量之和、单场降雨日降雨量 小于1.27mm的降雨天数、单场降雨日降雨量大于50mm的雨量之和、单场降雨日降雨 量大于50mm的降雨天数；各单场降雨的以上统计特征量构成该单场降雨的特征向量；

步骤B、对所有单场降雨的特征向量进行均值为0，方差为1的标准化；

步骤C、对标准化后的单场降雨的特征向量进行聚类，从而将各单场降雨归入相应的类 别中；

步骤D、对聚类得到的各类别的单场降雨场次数目进行归一化处理，统计每个站点的降 雨类型直方图；

步骤E、根据下式计算任意两个降雨站点：M和N的降雨类型直方图的相似度Sim（M,N）， 相似度越大表示两个降雨站点的相似性越高：

$\mathrm{Sim}(M,N)=\underset{i=1}{\overset{K}{\Sigma }}\mathrm{MIN}(a_i/A,b_i/B)$

式中，K为步骤C得到的单场降雨的类别数，a_i和b_i分别为降雨站点M、降雨站点N 的第i类单场降雨的归一化后的场次；A和B分别为降雨站点M、降雨站点N的单场 降雨场次总数。

上述技术方案中，所述聚类可采用现有的各种聚类算法，例如，K-means（K-均值） 聚类、模糊聚类（fuzzy clustering）、谱聚类（spectral clustering）等；其中，K-means （K-均值）聚类算法由于其算法简单易懂，计算速度快，更适用于本发明的技术方案。 但传统的K-means聚类算法需要用户事先给定聚类数目K，因此具有一定的局限性， 为此，本发明在对标准化后的单场降雨的特征向量进行聚类时，优选采用基于DBI指 数（Davies-Bouldin Index，Davies-Bouldin指数）的K-means自动聚类算法，具体 如下：

步骤C1、设置聚类类别数K的初始值为MinK，MinK为预设的2到10之间的自然数；

步骤C2、判断K是否小于等于预设的阈值MaxK，MaxK为大于MinK的自然数；如 果是，则随机从样本数据中取K个聚类中心，转向步骤C2，否则转向步骤C7；

步骤C3、分别计算各个样本点到各聚类中心的距离，并且将该样本点归到与其距离最 近的类中；

步骤C4、将所有样本点归类结束之后，重新计算这K个类别的聚类中心；

步骤C5、比较新计算的K个聚类中心与步骤C2中选取的聚类中心是否相同，若不同， 则转向步骤C3；否则，转向步骤C6；

步骤C6、计算此时的DBI指数，并使K=K+1，转向步骤C2；

步骤C7、比较所有的DBI指数值，将DBI指数值最小的K所对应的聚类结果作为最 终的聚类结果输出。

降雨序列实质是一个不连续的时间序列，大多数此前的研究，都是基于日月或年降 雨量单纯累积的宏观统计分析，没有从更细粒度分析单场降雨的统计特征。而单场降雨 序列恰恰包含了更丰富的降雨信息，能反映降雨序列的本质。例如单独对两个站点一年 的降雨进行统计，同样降雨量总和为1000mm的降雨，可能是8场单场降雨为125mm， 持续3天的降雨，也可能是20场单场降雨为50mm，持续1~2天的降雨，甚至还可能 是若干场不同降雨量，不同持续天数的降雨。若单纯利用年降雨量的累计而不考虑单场 降雨，则无法区分一年内两个站点降雨类型的差异，更无法区分降雨持续时间等信息。 因此，相比现有技术，本发明方法具有计算方法简单，但能包含降雨时间序列更多更细 粒度信息，真正反映降雨站点相似性的优点。

附图说明

图1为本发明方法的流程示意图；

图2为N1降雨站点的单场降雨类型直方图；

图3为N2降雨站点的单场降雨类型直方图。

具体实施方式

下面结合附图对本发明的技术方案进行详细说明：

本发明的基于单场降雨类型的降雨站点相似性评价方法，如图1所示，按照以下步 骤：

步骤1、对降雨站点一段时间的降雨历史数据进行预处理，去除那些重复冗余的日降雨 信息；并按照降雨间隔小于两天为同一场降雨的原则，对单场降雨进行划分。

步骤2、降雨站点的单场降雨特征的选择和提取。

特征的选择是一个非常关键的步骤，它将直接影响聚类结果的合理性，进一步影响 相似性比较的有效性。经过相关文献的总结分析，本发明选取以下9个特征统计量构成 单场降雨的特征向量：单场降雨雨量和rsum、单场降雨天数n，单场降雨日平均降雨量 rave、单场降雨日降雨量最大值rmax、单场降雨日降雨量最小值rmin、单场降雨日降 雨量小于1.27mm的雨量之和r1.27sum、单场降雨日降雨量小于1.27mm的降雨天数 r1.27、单场降雨日降雨量大于50mm的雨量之和r50sum、单场降雨日降雨量大于50mm 的降雨天数r50。

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

步骤3、采用基于DBI指数的K-means自动聚类算法对单场降雨的特征向量进行聚类。

K-means聚类算法由于其算法简单易懂，计算速度快，通常被作为大样本聚类分析 的首选方案，是最普遍的聚类方法之一，广泛应用于数据挖掘研究中。传统的K-means 算法是将含有n个数据点（实体）的集合X={x₁,x₂,...,x_n}划分为K个类簇C_j的问题 （j=1,2,...,K）。首先随机选取K个数据点作为K个类簇的初始簇中心，集合中每个数 据点被划分到与其距离最近的簇中心所在的类簇中，形成初始的K个聚类分布。对已 分配的每一个类簇计算新的簇中心，然后继续进行数据分配的过程，这样迭代若干次后， 若簇中心不再发生变化，说明数据对象已经全部分配到自己所在的类簇中。但是传统的 K-means算法存在一定的局限性，如聚类数目K需要用户事先给定等。本发明优选采 用基于DBI指数的K-means自动聚类算法，该方法利用DBI指数对聚类有效性进行评 价，能自动确定最佳聚类个数，它解决了传统的K-means算法需要事先确定K值的缺 点。DBI指数是衡量聚类效果时常用的一个指数，由类内散布和类间散布的比值表示， 比值越小表示聚类效果越好，具体内容可参见文献[Davies D.L.,Bouldin D.W..A  cluster separation measure.1979.IEEE Trans.Pattern Anal.Machine Intell.1 (4).224-227]。具体的，本发明所采用的基于DBI指数的K-means自动聚类算法包括 以下步骤：

（1）准备好聚类所需的样本数据，以及聚类类别数K的范围，本发明设K初始 值为MinK，MinK一般可在2到10之间取值，阈值MaxK可在15到30 之间取值，具体与参与聚类的数据多少有关，参与聚类的数据越多，MaxK 也可以相应提高；

（2）判断K是否小于等于MaxK，如果是，则随机从样本数据中取K个聚类中心， 转向（3），否则转向（7）；

（3）分别计算各个样本点到各个初始聚类中心的距离，并且将该样本点归到与其 距离最近的类中；

（4）将所有样本点归类结束之后，重新计算K个类别的聚类中心；

（5）比较新计算的聚类中心与原来的聚类中心，若聚类中心不同，则转向（3）， 否则转向（6）；

（6）计算DBI指数的值，并使K=K+1，转向（2）；DBI指数的计算公式如下：

$\mathrm{DBI}=\frac{1}{K}\underset{i=1}{\overset{K}{\Sigma }}{\max }_{j\ne i}\left\{\frac{S_i+S_j}{d_{i,j}}\right\}---\left(1\right)$

其中，S_i表示第i个簇内样本与该簇中心之间的平均距离，S_j表示第j个簇内样本 与该簇中心之间的平均距离。d_i,j表示第i类样本与第j类样本之间的相异度（也即分 散程度）。K表示聚类的个数。

同时，S_i的计算公式如下：

$S_i=\frac{1}{\left|C_i\right|}\underset{x\in C_i}{\Sigma }\left|\right|x-{\nu }_i\left|\right|---\left(2\right)$

其中，C_i表示属于第i簇的样本个数,v_i表示第i个聚类的类别中心，x表示属于 第i簇的各个样本。

S_j的计算公式类似公式2。

d_i，j的计算公式如下：

d_i,j=||v_i-v_j‖                                                （3）

v_i表示第i个聚类的类别中心，v_j表示第j个聚类的类别中心。

式（1）~（3）中所用的距离都是指欧几里德距离。关于DBI指数更详细的内容 可参见文献[Davies D.L.,Bouldin D.W..A cluster separation measure.1979.IEEE  Trans.Pattern Anal.Machine Intell.1(4).224-227]。

（7）比较所有的DBI值，将DBI指数值最小的K所对应的聚类结果作为最终的 聚类结果输出。

为了保证每一个特征向量对聚类结果的贡献率不会因为数值的大小而有显著的差 异，本发明首先对所有特征向量进行均值为0，方差为1的标准化，然后再进行基于 DBI指数的K-means自动聚类。这样，各单场降雨根据其特征而被归入相应的类别。

在信息检索技术中，有一个词袋模型（Bag of words model）的概念，它是自然语 言处理和信息检索中的一种有效模型。在这种模型中，文本（段落或者文档）被看作是 无序的词汇集合，忽略语法甚至是单词的顺序。其概念来源于文本信息检索和分类任务 中的文档表示技术，此技术将文档表示为一个用词汇出现次数所组成的矢量。。为了便 于理解本发明，类似地，也可结合词袋模型的特点，将单场降雨对应的聚类结果类别号 看作一个单词，相同类别号的单场降雨落入相同的降雨类型袋子中，这种新的模型被命 名为降雨类型词袋模型（Bag of precipitation types）。由于可以用类似图像分类中直方图 的形式表示出来，即每一个类别号所对应的单场降雨场次为bin的高度值，每一个降雨 站点对应的各个类别单场降雨场次的组合便构成了一个直方图，即单场降雨类型直方 图。

步骤4、根据单场降雨类型评价任意两个降雨站点之间的相似性。

先将每一个降雨站点各类别对应的单场降雨场次进行归一化，然后利用直方图相交 的方法获取两个不同降雨站点之间的相似性。本发明中所用两站点直方图相交 (histogram intersection)算法如（式4）所示：

$\mathrm{Sim}(M,N)=\underset{i=1}{\overset{\mathrm{BIN}}{\Sigma }}\mathrm{MIN}(a_i/A,b_i/B)$（式4）

其中，M、N为两个不同的降雨站点；Sim（M,N）为降雨站点M、N之间的相似度， 其值越大表明相似性越高；BIN为单场降雨直方图中桶(bin)的个数（亦即聚类得到的单 场降雨类型数目K）；a_i和b_i分别为降雨站点M、降雨站点N的第i类单场降雨的归一 化后的场次；A和B分别为降雨站点M、降雨站点N的单场降雨场次总数。

两站点相似度越高，两站点相似度接近于1；两站点差异度越大，两站点相似度接 近于0。

为了验证本发明方法的效果，利用江苏省水文水资源勘测局提供的江苏省各地区雨 量站点2000年至2011年日降雨数据，进行了以下验证试验：

数据预处理阶段：首先挑选三块典型区域站点，分别以北部（N）、中部(M)，南部 (S)表示；然后将这三个区域中有降雨数据的站点数据提取出来，其对应的站点代号和 地名如下表1，其中北部站点6个，中部站点6个，南部站点8个，这些站点数据相对 较为完整。

表1.三块降雨区域所对应的站点

  北部站点名   站点代号   中部站点名   站点代号   南部站点名   站点代号   小塔山水库   N1   盱眙   M1   赤山新闸   S1   青口   N2   三河闸   M2   大溪水库   S2   石梁河水库   N3   龙王山水库   M3   高淳   S3   大兴镇   N4   桂五水库   M4   句容   S4   临洪   N5   化农水库   M5   溧阳   S5   小许庄   N6   红旗水库   M6   沙河水库   S6   天生桥闸   S7   王母观   S8

提取每个站点10年数据中的雨量项，去除那些重复冗余的日降雨信息。按照降雨 间隔小于两天为同一场降雨的原则，对单场降雨进行划分。再根据已经划分好的降雨场 次计算单场降雨和rsum，单场降雨天数n，单场降雨日平均降雨量rave，单场降雨日降 雨量最大值rmax和最小值rmin，单场日降雨量小于1.27mm雨量和r127sum与降雨天 数r127，单场日降雨量大于50mm的雨量和r50sum与降雨天数r50。

将这九个统计特征量作为单场降雨的特征向量，并且对三个区域的20个站点的所 有特征向量进行均值为0，方差为1的标准化，这样可以保证每一个特征值对聚类结果 的贡献率不会因为数值的大小而有显著的差异，再对这些单场降雨特征向量进行基于 DBI指数的K-means自动聚类，单场降雨类型被自动聚类成14类。

统计聚类结果，得到各站点的降雨类型直方图，例如北部N1、N2站点得到的单场 降雨类型直方图分别如图2、图3所示。

再利用降雨类型直方图相似性模型对站点的相似性进行分析，包括同一区域和不同 区域间站点的相似性分析，同一区域站点相似度如下表2、3、4所示。

表2北部区域内部站点相似度

  站点相似度   N1   N2   N3   N4   N5   N6   N1   1   0.903074   0.904808   0.883736   0.866606   0.787833   N2   0.903074   1   0.90049   0.892857   0.923903   0.826891   N3   0.904808   0.90049   1   0.914951   0.901389   0.808824   N4   0.883736   0.892857   0.914951   1   0.912932   0.79902   N5   0.866606   0.923903   0.901389   0.912932   1   0.807096   N6   0.787833   0.826891   0.808824   0.79902   0.807096   1

如表2为北部的6个站点间的相似度，表格对角线上的数值为“1”，这是由于对角 线处的相似度是同一个站点的比较，完全相同，所以相似度为1。

结合表2~4中数据分析如下：

在同一区域内，站点的相似度应该较高，N3与N1-4之间由于空间上较近，相似度 都高于0.9，符合模型的评价标准。N1与N2-5站点的相似度也都接近0.9，表明N1与 大部分区域内站点的相似度符合评价标准。由于N6在空间上偏离其他站点，所以与其 他站点的相似度相对较低，但N6与其他站点相似度仍然接近0.8。

表3中部区域内部站点相似度

  站点相似度   M1   M2   M3   M4   M5   M6   M1   1   0.88308   0.833835   0.834359   0.839372   0.719073   M2   0.88308   1   0.835478   0.81697   0.85122   0.754372   M3   0.833835   0.835478   1   0.84787   0.867531   0.720929   M4   0.834359   0.81697   0.84787   1   0.852214   0.770502   M5   0.839372   0.85122   0.867531   0.852214   1   0.836296   M6   0.719073   0.754372   0.720929   0.770502   0.836296   1

如表3为中部的6个站点间的相似度，M6站点由于地理上相对M1-4站点较远， 所以这部分相似度较低，略低于0.8，而与M6较近的M5站点，相似度仍然达到了0.8 以上。大部分站点相似性都在0.8以上，正确率几乎达到了100%。

表4南部区域内部站点相似度

  站点相似度   S1   S2   S3   S4   S5   S6   S7   S8   S1   1   0.920763   0.884571   0.890932   0.875543   0.924789   0.907639   0.897134   S2   0.920763   1   0.859095   0.865851   0.850829   0.904446   0.882869   0.879204   S3   0.884571   0.859095   1   0.890278   0.921974   0.890929   0.933477   0.826083   S4   0.890932   0.865851   0.890278   1   0.911865   0.892214   0.914057   0.850413   S5   0.875543   0.850829   0.921974   0.911865   1   0.901122   0.94101   0.84481   S6   0.924789   0.904446   0.890929   0.892214   0.901122   1   0.89213   0.896791   S7   0.907639   0.882869   0.933477   0.914057   0.94101   0.89213   1   0.849597   S8   0.897134   0.879204   0.826083   0.850413   0.84481   0.896791   0.849597   1

如表4中为南部区域站点间的相似度，数据显示南部站点的相似度普遍在0.85以 上，达到了高度的相似，相似度正确率也几乎达到了100%，这是由于南部站点靠近长 江，又有众多湖泊，区域降雨普遍较多，所以站点相似度较高。南部区域内部站点的高 度相似性验证了模型在区域内部相似性度量的可靠性。

不同区域站点相似性比较如下表5，6，7所示：

表5北部站点与中部站点相似性度

  站点相似度   M1   M2   M3   M4   M5   M6   N1   0.902992   0.837353   0.827357   0.82547   0.824003   0.709109   N2   0.882676   0.870861   0.886555   0.843642   0.847008   0.710314   N3   0.877564   0.825823   0.826942   0.864167   0.825509   0.711473   N4   0.862519   0.817446   0.807718   0.853725   0.832098   0.723207   N5   0.862576   0.824236   0.839599   0.816825   0.825518   0.678734   N6   0.790498   0.816485   0.843919   0.832157   0.900449   0.79452

如表5中为北部与中部站点间的相似度。北部和中部站点都位于长江以北，所以这 两个区域在降雨类型上呈现出一定的相似性。

表6北部站点与南部站点相似度

  站点相似度   S1   S2   S3   S4   S5   S6   S7   S8   N1   0.729325   0.742772   0.762788   0.78929   0.758561   0.74336   0.766783   0.664185   N2   0.755944   0.765987   0.79899   0.811591   0.797151   0.757487   0.797634   0.689482   N3   0.764918   0.766844   0.743071   0.807398   0.768543   0.763703   0.775881   0.711156   N4   0.747859   0.745411   0.747522   0.806823   0.764641   0.790696   0.743306   0.701138   N5   0.76567   0.753377   0.785245   0.837869   0.834148   0.782092   0.814349   0.731695   N6   0.789299   0.779099   0.778208   0.802734   0.787651   0.76668   0.796761   0.720472

如表6中为北部与南部站点相似度，由于南部和北部空间上相距相对较远，而且相 隔长江，降雨有较大的差异，站点的相似性普遍在0.8以下甚至只有0.6，相似度正确 性接近90%。

表7南部与中部站点相似度

  站点相似度   S1   S2   S3   S4   S5   S6   S7   S8   M1   0.751155   0.76593   0.803038   0.814678   0.800017   0.775839   0.806066   0.691108   M2   0.759015   0.758697   0.82703   0.797989   0.817501   0.780566   0.814873   0.693484   M3   0.741128   0.747001   0.770553   0.783148   0.779698   0.73611   0.789507   0.686717   M4   0.754152   0.760851   0.718512   0.787483   0.742252   0.748456   0.756946   0.71957   M5   0.783863   0.759731   0.808943   0.810009   0.791062   0.794841   0.797037   0.719664   M6   0.779466   0.731079   0.751476   0.73401   0.749433   0.76933   0.754903   0.698974

如表7中为南部与中部站点相似度，由于不同区域空间上的差异，超过85%以上 的站点相似性在0.8以下，验证了模型在不同区域内相似度的正确性。

从以上站点相似度比较的表格中，除个别站点外，同一区域内站点相似度普遍在 0.8甚至0.9以上，不同区域内站点的相似度普遍在0.8甚至接近0.6；地理空间相对较 远的北部与南部站点相似度小于地理空间相对较近的北部与中部站点相似度；而且中部 与南部站点相似性普遍大于中部与北部的站点相似度，这是由于以长江为分水岭，南北 区域的降雨有较大差异，这和实际情况比较吻合，从而验证了本发明方法的正确性。