CN103729550B - 基于传播时间聚类分析的多模型集成洪水预报方法 - Google Patents

Abstract

本发明公开了一种基于传播时间聚类分析的多模型集成洪水预报方法，属于水文预报技术领域。首先采用派生的动态时间弯曲匹配方法进行洪水过程相似性分析，估计上下游各站点的流量传播时间，并通过对流量传播时间进行聚类分析将样本分解为若干簇，然后分别对子流量序列建立SVM回归模型模拟洪水形成过程，最后再将这些子模型合并成一个综合模型。将该方法的综合预测结果与常规条件下的单一模型和基于流量聚类的模型预测结果相比较，结果显示该模型综合表现更佳。

Description

基于传播时间聚类分析的多模型集成洪水预报方法

技术领域

本发明涉及一种洪水预报方法，尤其涉及一种基于传播时间聚类分析的多模型集成洪水预报方法，属于水文预报技术领域。

背景技术

在水文时间序列预测分析中，径流预报是一类经典的水文问题，对水资源调配管理和防洪减灾调度决策有着重要的意义。然而由于上游站点流量以及气候与下垫面等因素的综合影响，流量往往表现出复杂的非线性和非平稳性特征。

目前已经有很多进行水文时间序列预测的方法，其中最为简单的是线性模型预测方法，如AR、ARMA等。但应用线性模型预测非线性时间序列很难取得良好的效果。非线性模型由于其本身非线性特性适用于非线性时间序列建模，如人工神经网络、二次回归模型、混沌模型、支持向量机回归模型等。由于单个水文预报模型往往强化了水文预报的某些方面而忽视了另一些方面，文献[张驰，周惠成，李伟.基于数据分析技术的水文组合预报应用研究，大连理工大学学报，2007，47(2)：246-251]、文献[Kunhui Lin,Qiang Lin,ChangleZhou,et al.Time Series Prediction Based on LinearRegression and SVR.ThirdInternational Conference on Natural Computation,ICNC 2007:688-691]采用多模型结合预报的方法，从多角度模拟流域状况，得到更符合实际情况、更加理想的预测效果。在实际应用中，不同的洪水预报模型有不同的适用范围，甚至在同一流域，不同时间情况下，最适合采用的模型也应是不同的。目前，已有一些研究采用不同的方法对流量或水位过程进行分解聚类，提高神经网络预报精度。尹雄锐等(尹雄锐，张翔，夏军.基于聚类分析的人工神经网络洪水预报模型研究，四川大学学报(工程科学版)，2007，39(3)：34-40)应用模糊C均值和自组织映射网络对洪水流量聚类成不同的类别，然后采用相同的输入分别建模、预测。胡铁松等(胡铁松，丁晶.径流长期分级预报的Kohonen网络方法.水电站设计，1997(6)：13(2))提出了径流长期分级预报的一种新的模式识别方法，通过Kohonen自组织神经网络对历史样本的学习，识别出了蕴含在样本中径流级别与其因子之间的规律性，并证明了该方法的有效性。王玲和黄国如(王玲，黄国如.基于径流分类的日径流量预测神经网络模型[J].灌溉排水，2002，21(4)：45-48)通过对将时间分为枯水期、2个湿润期和一个丰水期，进而达到对流量分类的目的，建立的综合神经网络模型较单一的神经网络模型有更好的精度。

不同的流量，水流速度不一样，导致洪水传播时间有所变化，且不同的时间段雨量大小变化较大，最后导致降雨径流预报输入输出关系的侧重点有所不同，如果只用单一的模型，很难将各部分的数据同时拟合得很好。

发明内容

本发明所要解决的技术问题在于克服现有技术不足，提供一种基于传播时间聚类分析的多模型集成洪水预报方法，能有效的提高整体预测能力，特别是对洪水期流量较大的情况下的预测精度有显著的提高。

本发明的基于传播时间聚类分析的多模型集成洪水预报方法，首先根据历史流量/水位数据建立包含多个子模型的洪水预报综合模型，然后利用所建立的洪水预报综合模型进行洪水预报；具体包括以下步骤：

步骤1、对于目标站点的历史流量/水位数据中的每个流量/水位数据，利用各上游站点与目标站点之间的流量传播时间构造对应于该流量/水位数据的传播时间向量，目标站点的历史流量/水位数据中的所有流量/水位数据所对应的传播时间向量构成传播时间向量集合；

步骤2、对所述传播时间向量集合进行聚类，得到k个簇，并对聚类得到的每个簇分别进行以下处理：

对于每一个上游站点，统计该簇中所包含的该上游站点与目标站点之间的流量传播时间的分布情况，并将出现频次最低的部分流量传播时间剔除，剩余的流量传播时间作为该簇中该上游站点的可用流量传播时间，最终得到该簇中各上游站点的可用流量传播时间；

步骤3、构建k个训练样本集，k个训练样本集与步骤2得到的k个簇一一对应，k为聚类得到的类别数；其中任意一个训练样本集中的训练样本按照以下方法得到：

选取一组目标站点的历史流量/水位数据分别作为该训练样本集中各训练样本的输出；对每一个训练样本的输出，根据其所属训练样本集所对应的簇中各上游站点的可用流量传播时间，从各上游站点的历史流量/水位数据中确定相应的流量/水位数据，并结合相应的雨量输入信息及目标站点的其它输入历史流量/水位数据，构成该训练样本的输入；步骤4、利用所构建的k个训练样本集各自对预测模型进行训练，得到k个预测子模型，这k个预测子模型共同构成洪水预报综合模型；

步骤5、根据测试样本按照以下方法从洪水预报综合模型中选择相应的预测子模型：对测试样本中目标站点的流量/水位数据，利用各上游站点与目标站点之间的流量传播时间构造对应于该流量/水位数据的传播时间向量，并从步骤2中聚类得到的各个簇中选出簇中心与该传播时间向量的距离最小的簇，距离最小的簇所对应的预测子模型即为所选择的预测子模型；

步骤6、以所述测试样本作为输入，利用所选择的预测子模型进行目标站点的洪水预报。

优选地，所述各上游站点与目标站点之间的流量传播时间利用派生动态时间弯曲算法获取。进一步地，在利用派生动态时间弯曲算法获取各上游站点与目标站点之间的流量传播时间时，针对上游站点流量/水位的特征点，若该特征点是极小点，则将该特征点所匹配的目标站点若干流量/水位数据点中出现极小值的那一个作为上游站点流量/水位的特征点的唯一匹配；若是极大点，则将该特征点所匹配的目标站点若干流量/水位数据点中出现极大值的那一个作为上游站点流量/水位的特征点的唯一匹配；对于非极值点出现多个匹配情况，则采用与该点之前最近的极值点传播时间最相近的传播时间。

优选地，所述聚类使用基于DBI指标的K-means聚类方法。

优选地，步骤2中所述将出现频次最低的部分流量传播时间剔除，具体是指将每个上游站点所对应的流量传播时间中出现频次最低的10％的流量传播时间剔除。

本发明采用基于传播时间来划分模型，通过聚类将整个序列分解为若干个小组，对于不同传播时间分布的流量建立不同的子模型，并根据样本传播时间分布选择相应的子模型进行洪水预报。相比现有技术，本发明方法简单有效，不仅可以准确估算出各模型的洪峰传播时间，还有利于帮助准确确定各模型输入，最终的预测精度比较理想。

附图说明

图1为本发明方法的流程示意图，其中左侧为训练过程流程，右侧为测试过程流程；

图2为淮河潢川站与王家坝站2006年一段洪水过程的DDTW匹配结果；

图3为DBI指标随聚类类别数目变化的波动情况；

图4(a)～图4(c)依次为息县、潢川、班台三个站点的传播时间直方图。

具体实施方式

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

本发明针对现有的单一时间序列预测模型的不足，以及洪水预报时经常出现的高流量峰值预测不准以及预测“延时”的现象，提出一种基于传播时间聚类分析的多模型集成洪水预报方法，和常规的流量分类预测思想不同，本发明采用基于传播时间来划分模型，通过聚类将整个序列分解为若干个小组，对于不同传播时间分布的流量建立不同的子模型。

本发明基于传播时间聚类分析的多模型集成洪水预报方法与现有技术相同，都包括训练阶段和测试阶段，其流程如图1所示，具体包括以下步骤：

步骤1、对于目标站点的历史流量/水位数据中的每个流量/水位数据，利用各上游站点与目标站点之间的流量传播时间构造对应于该流量/水位数据的传播时间向量，目标站点的历史流量/水位数据中的所有流量/水位数据所对应的传播时间向量构成传播时间向量集合。

洪水预报，不仅需要准确预报洪峰的值，更要能准确掌握洪水到达的时间。所谓流量传播时间是指上游站点的某时刻流量传到下游站点所需的时间。流量传播时间的确定是进行流量预报、洪水演进过程分析、水量调度等的基础，也是流量预报的一项重要内容。以往的研究大多是根据经验由水文业务人员提供的平均传播时间直接确定输入信息来预测洪峰流量的值，但是实际上随着自然和社会条件的变化，洪水传播时间也会逐渐发生变化，特别是不同量级洪水的传播时间差异更大。

近年来，有一些学者对流量传播时间的算法进行研究，已有的算法包括：相关水位法、灰色关联分析法、分段积分法、相关分析法等。而在实际洪水传播过程中，影响因素较多，如流速、河床粗糙率等。用水文学上的传统理论方法去计算洪水传播时间非常复杂，而且效果也不好。对洪水传播时间的研究其实可以归纳为时间序列研究中的相似性搜索问题。时间序列相似性搜索常用的一个经典算法是动态时间弯曲匹配(Dynamic TimeWarping,DTW)，它是基于动态规划的一种模式匹配方法，解决了欧式距离在进行相似性度量时时间序列长度必须一致的问题，同时具备时间轴和幅度上的伸缩能力。

如果两个序列只在时间轴上局部不一致，DTW可以取得理想的匹配效果。当两个序列同时在Y轴上不一致时，如果是序列全局的不一致，比如不同的均值，不同的幅度等，可以通过偏移量转化，幅度标准化等一些预处理解决。但当两个序列在Y轴上出现局部的不一致时，DTW的匹配就出现了问题。为了使序列的特征点(如峰、谷等)很好的匹配，文献(李士进，张晓花，万定生等.基于DTW的测站水位影响关系估计，江南大学学报(自然科学版)，2007，6(6)：678-682)引入派生动态时间弯曲(Derivative Dynamic Time Warping,DDTW)。本发明优选采用派生动态时间弯曲算法(DDTW)来确定各上游站点与目标站点之间的流量传播时间。在介绍DDTW之前，我们先给出DTW的匹配原理。

设有两个时间序列Q和C，长度分别为n和m。

Q＝(q₁,...,q_i,...,q_n)

C＝(c₁,...,c_j,...,c_m) (1)

为利用DTW将两个时间序列对准，首先构造一个n行m列矩阵M，矩阵中的元素(i,j)为两时间序列数据中对准点q_i和c_j之间的距离d(q_i,c_j)，计算公式如下：

d(q_i,c_j)＝(q_i-c_j)² (2)

弯曲路径W是矩阵中邻近元素的集合。

W＝ω₁,ω₂,...,ω_k,...ω_K max(m,n)≤K≤m+n-1 (3)

它是序列Q与C之间的一个映射。DTW距离取弯曲路径总长度的最小值，即

$d_{DTW}(Q,C)=min\[\sqrt{\underset{k=1}{\overset{K}{\Σ}}{\mathrm{\ω}}_k}\]---\left(4\right)$

最佳路径可以由时间起始点(1，1)到终点(m，n)之间的局部最优解通过递归获得，公式如下：

$\begin{array}{c}\mathrm{\γ}\left(i,j\right)=d\left(q_i,c_j\right)+\\ min\{\mathrm{\γ}(i-1,j-1),\mathrm{\γ}(i-1,j),\mathrm{\γ}(i,j-1)\}\end{array}---\left(5\right)$

式中γ(i,j)为累加距离，由当前对准点的距离和相邻点的累加DTW距离计算得到。

弯曲路径W满足如下条件：

(1)边界条件：ω₁＝(1,1)，ω_K＝(m,n)。弯曲路径从左下角出发终止于右上角。

(2)连续条件：假如要计算点(i,j)，它的前驱点必须是(i-1,j-1)，(i-1,j)或(i,j-1)。

(3)单调性：ω_k＝(a,b)，那么ω_k-1＝(a',b')必须满足a-a'≥0，b-b'≥0。

(4)斜率限制：弯曲路径不能太陡，防止很短的序列匹配很长的序列。

传统的DTW方法中，n行m列的矩阵中元素(i,j)表示两个点q_i和c_j之间的距离d(q_i,c_j)(如公式2)。而在DDTW中，d(q_i,c_j)是q_i和c_j的派生距离，由q_i派生出D_x[q]如下：

$D_x\[q\]=\frac{(q_i-q_{i-1})+\left(\right(q_{i+1}-q_{i-1})/2)}{2}---\left(6\right)$

它是两个斜率的平均，分别为当前数据与其左邻居所在直线的斜率，当前数据左右邻居所在直线的斜率。经验上，采用这个派生距离进行DTW匹配比任何只考虑两个数据点的方法具有更高的鲁棒性。

下面给出基于DDTW计算上、下游站点之间每个样本流量传播时间的方法。

采用DDTW方法对这两测站的水位或流量过程线进行相似性匹配。具体做法如下：分别计算当前洪水过程该样本之前一段子序列两站水位或流量的DDTW距离(本具体实施方式中取15个时刻点)，并记录两站水位或流量线之间的匹配关系，即上游的某一时刻水位或流量经DDTW匹配后，与下游哪一时刻或哪几个时刻相匹配。在计算DDTW相似距离时，需要结合水文知识，假设不存在回游的情况，则上游的某一时刻只能和下游该时刻以后的某一时刻相匹配，因此对DDTW作了改进限制其弯曲路径，即上游的点只能和下游该时刻之后的点相匹配。

由于上下游的水位或流量影响关系在极值点处比较明显，即上游的某时刻出现了水位或流量的极小点或是极大点，而一般的水位或流量影响关系不大，因此，需要提取上游水位或流量过程线的特征点，方法如下：采用一个滑动窗口，窗口大小为R(R是奇数)。当窗口的中心滑动到每一点时，若它是局部极小点，则窗口中左端的(R-1)/2个点都大于该点，并且这些点的值从左至右递减，窗口中右端的(R-1)/2个点也都大于该点，并且这些点的值从左至右递增。局部极大点的获取方法类似。

实际的水文关系中，上下游的影响关系是唯一对应的，但由于DDTW算法本身的原因，上游的某一时刻可能对应于下游的多个时刻，或者是下游的某一时刻对应于上游的多个时刻。为得到唯一的匹配关系，本发明处理的方法是：针对上游站点流量/水位的特征点，若该特征点是极小点，则将该特征点所匹配的目标站点若干流量/水位数据点中出现极小值的那一个作为上游站点流量/水位的特征点的唯一匹配；若是极大点，则将该特征点所匹配的目标站点若干流量/水位数据点中出现极大值的那一个作为上游站点流量/水位的特征点的唯一匹配；对于非极值点出现多个匹配情况，则采用与该点之前最近的极值点传播时间最相近的传播时间。这样，就得到了上下游的唯一匹配关系。

最后，对这些匹配关系进行统计。对于上下游水位或流量线的每组匹配关系，分别计算匹配时间偏差(即上游站点过程线某时刻样本与下游站点过程线匹配点的时刻差)，该时间偏差即为上下游测站之间在该样本点的流量传播时间。图2给出了淮河潢川站与王家坝站2006年一段洪水过程DDTW匹配结果。

假设目标站点的流量是由m个上游站点的流量汇流而成，则目标站点的某一个流量或水位数据样本对应着m个流量传播时间，这m个流量传播时间即可构成一个对应于该流量/水位数据的m维的传播时间向量；目标站点的一系列流量/水位数据所对应的传播时间向量构成传播时间向量集合。

步骤2、对所述传播时间向量集合进行聚类，得到k个簇，并对聚类得到的每个簇分别进行以下处理：

对于每一个上游站点，统计该簇中所包含的该上游站点与目标站点之间的流量传播时间的分布情况，并将出现频次最低的部分流量传播时间剔除，剩余的流量传播时间作为该簇中该上游站点的可用流量传播时间，最终得到该簇中各上游站点的可用流量传播时间。

由于不同流量时期，水文情况有所差异，传播时间也不同，因此本发明将各类情况区分开来，分别为其建立预测(预报)模型。模型的划分可以归结为一个聚类问题。聚类是一种将给定的输入集合进行划分，使得其生成的子类具有相同的属性的方法。然而，聚类的初始中心和聚类数目一直是聚类过程中的一个重要环节，无论采用什么算法进行聚类，最终的聚类结构都需要进行有效性验证。

一般地，衡量聚类有效性需要通过两个方面来衡量，一是要求类内部尽可能地紧凑，同一数据簇内数据应有较高的紧密度；二是要求类与类之间的距离尽可能地远，不同簇间的分离度应尽可能地大。目前，计算聚类有效性的指标有很多，常用的主要有Silhouetteindex、Dunns based index、Davies-Bouldin index、基于类间相关系数的聚类有效性指标和引进惩罚因子的V_xie指标。Davies-Bouldin Index(DBI)指标与其他方法相比，对异常值和边界点不敏感，同时适用于聚类个数超过两类的时候同样能指导正确的分类。所以本发明选用了DBI指标来评价聚类的结果。DBI的计算公式如下所示：

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

其中，S_i表示第i类样本之间的相似度，S_j表示第j类样本之间的相似度；d_i,j表示第i类样本与第j类样本之间的相异度；K表示聚类的个数。

同时，S_i的计算公式为：

d_i,j的计算公式为：

d_i,j＝||v_i-v_j|| (9)

其中，C_i表示属于第i类的样本个数；v_i表示第i个聚类的类别中心；x表示属于第i类的各个样本。当类内相似性小，类大相异性大的时候，DBI值越小，表示聚类效果越好，因此，好的聚类应该是类内距离与类间距离的比值尽可能小，我们通过计算不同聚类数k时的DBI值，将最小DBI值对应的聚类数k作为最佳聚类结果。

本发明所采用的基于DBI指标的K-means聚类方法步骤如下：

(1)准备好聚类所需的样本数据，这里设K初始值为5，Max K＝22；

(2)判断是否K<＝Max K，如果是，则随机从训练样本中取K个聚类中心，转向(3)，否则转向(7)；

(3)分别计算各个样本点到各个聚类中心的距离，并且将该样本点归到与其距离最近的簇中去；

(4)将所有样本点归类完之后，重新计算K个簇的聚类中心；

(5)比较新计算的聚类中心与原来的聚类中心之间距离，如果该距离大于某个阈值，则用新的聚类中心替代原来的聚类中心，转向(3)，否则转向(6)；

(6)计算并保存DBI的值，然后K＝K+1，转向(2)；

(7)比较所有的DBI值，将最小的DBI值所对应的K作为最后输出；

基于DBI的K-means聚类算法自动的确定了最佳聚类个数，解决了普通K均值算法需要事先确定K值的缺点。由于传播时间本身就具有很大的不确定性，把聚类分析运用到洪水预报的洪水分类中能使分类更切合实际，所得到的结果比传统的分类方法也更合理，也为最终的的多模型预报中的子模型建立奠定了基础。

通过上述基于DBI的K-means聚类算法将传播时间向量集合分解为k种类型，即得到每一个训练样本所对应的类别(簇)。如果参与聚类的是m维上游站点的传播时间，那生成的k个子类同样是关于m个站点的，并且，各类的聚类中心就代表着各类中各站点的主要传播时间。但研究发现如果只取一个主要传播时间作为每个上游站点的输入，则由于聚类得到的每个簇中各样本传播时间不尽完全相同从而会降低后续训练模型的预测准确度。因此本发明进一步提出通过统计每个簇中每个站点传播时间的局部分布(例如采用直方图分析方法)，并将出现频次最低的部分流量传播时间剔除，剩余的流量传播时间作为该簇中该上游站点的可用流量传播时间，最终得到该簇中各上游站点的可用流量传播时间。本发明选取包含90％以上的传播时间作为该簇对应上游站点的输入，即将每个上游站点所对应的流量传播时间中出现频次最低的10％的流量传播时间剔除。

步骤3、构建k个训练样本集，k个训练样本集与步骤2得到的k个簇一一对应，k为聚类得到的类别数；其中任意一个训练样本集中的训练样本按照以下方法得到：

选取一组目标站点的历史流量/水位数据分别作为该训练样本集中各训练样本的输出；对每一个训练样本的输出，根据其所属训练样本集所对应的簇中各上游站点的可用流量传播时间，从各上游站点的历史流量/水位数据中确定相应的流量/水位数据，并结合相应的雨量输入信息及目标站点的其它输入历史流量/水位数据，构成该训练样本的输入。

通过步骤2的处理，即可得到分别用于训练k个子模型的k个训练样本集，每个训练样本的输出以及各上游站点的流量/水位输入数据均已确定，此时还需要确定训练样本中的降雨输入信息以及目标站点的预见期前的历史流量/水位数据。由于降雨信息并未参加聚类，无法通过直方图分布的方法确定，本发明通过降雨序列与目标站点流量/水位序列的相关系数分析方法确定雨量输入信息，即通过降雨自身不同预见期的序列以及累加雨量信息与目标站点流量/水位序列相关系数取最大值时的两个降雨量和累加降雨量作为降雨的输入(两个序列的相关系数计算方法可参见文献：朱跃龙,李士进,范青松,万定生.基于小波神经网络的水文时间序列预测[J].山东大学学报(工学版).

2011,41(04):119-124)。由此即最终确定了各预测子模型所需的训练样本集。

步骤4、利用所构建的k个训练样本集各自对预测模型进行训练，得到k个预测子模型，这k个预测子模型共同构成洪水预报综合模型。

本发明可采用现有的各种预测模型，例如人工神经网络模型(钟登华,王仁超,皮钧.水文预报时间序列神经网络模型[J].水利学报.1995(02):69-75)或支持向量机回归模型(董辉,傅鹤林,冷伍明.支持向量机的时间序列回归与预测[J].系统仿真学报.2006,18(07):1785-1788)。以支持向量机(support vector machine,SVM)回归模型为例，用支持向量回归模型进行流量预测时，输入是上游各测站的流量和降雨信息，目标输出则是下游站点的当前流量。目标站点的流量可以写成各上游站点前期流量及降雨信息的函数，目标站点预测可以用如下表达式表示：

Q＝f(Q-t_i,P_i-τ_i,S) (10)

这里，Q代表的是目标站点输出，f为待训练的未知非线性函数，P_i(i＝1,2,...,n)为上游第i个站点的当前流量序列，t_i(i＝1,2,...,n)为目标站点的预见期，τ_i(i＝1,2,...,n)为第i个站点流量汇至下游目标站点的传播时间，Q-t_i为目标站点前t_i时流量，P_i-τ_i表示第i站点前τ_i小时流量，S代表雨量信息。

步骤5、根据测试样本按照以下方法从洪水预报综合模型中选择相应的预测子模型：对测试样本中目标站点的流量/水位数据，利用各上游站点与目标站点之间的流量传播时间构造对应于该流量/水位数据的传播时间向量，并从步骤2中聚类得到的各个簇中选出簇中心与该传播时间向量的距离最小的簇，距离最小的簇所对应的预测子模型即为所选择的预测子模型。

待得到各个训练子模型后，就要将新输入的测试样本数据“分配”到对应的模型中去。确定测试样本所对应子模型的基本思想是训练的子模型必须测试传播时间相同类型的测试样本。也就是说，测试样本与训练样本存在传播时间的“相似性”。本发明中对象之间的相异程度是基于欧式距离来计算的，距离越近就越相似。因此，对于新输入的测试样本，根据该样本所对应的传播时间向量，分别计算其与各子模型的传播时间聚类中心的距离，选择距离最近的子模型作为该样本的子模型。

步骤6、以所述测试样本作为输入，利用所选择的预测子模型进行目标站点的洪水预报。

采用步骤5所确定的对应的子模型对该新输入样本进行未来时刻流量预测。若新输入样本同时与多个子模型距离相近，还可采用加权组合的方法来进行预测。

下面以一个实例来验证本发明效果。

王家坝站作为淮河入皖第一个水文站，控制流域面积30630km²，下属3个雨量站、3个地下水位站、2个水质监测断面，承担王家坝淮干、官沙湖钐岗分洪道、地理城等7个断面的流量监测任务，并且每年的七八月份都会成为防洪防汛工作的重点。通过王家坝闸，可以有效削减淮河洪峰，减轻淮河中游压力。因此，及时掌握王家坝的洪水过程对防洪决策具有重大意义。

本实例以王家坝站汛期日流量作为预测对象，选取与之相关的上游三个观测站：班台、息县、潢川(以下分别简称为：xx、hc、bt)，此三站分别为上游三大支流上的站点。数据时间跨度为1998年-2008年，具体为1998年到2008年每年的6月1日到9月30日，每天采样4个数据，间隔为6小时。共80场次大中小型洪水，选取其中1998-2006年的洪水时期数据作为训练集，2007-2008的洪水时期数据作为测试集。

结果评价指标：

对于时间序列预测，常用的评价标准有均方误差(MSE)、确定性系数(DC)、相关系数(correlation coefficient)。

1.均方误差(MSE)

$MSE=\frac{1}{n}\underset{i=1}{\overset{n}{\&Sigma;}}{(p_i-q_i)}^2---\left(11\right)$

上式中，p_i、q_i分别代表预测值和实际值。

2.确定性系数(DC)

洪水预报过程与实测过程之间的吻合程度可用确定性系数作为指标，其计算公式如下：

$DC=1-\frac{\underset{i}{\overset{n}{\&Sigma;}}{\&lsqb;y_c\left(i\right)-y_0\left(i\right)\&rsqb;}^2}{\underset{i=1}{\overset{n}{\&Sigma;}}{\&lsqb;y_0\left(i\right)-\stackrel{\&OverBar;}{y_0}\&rsqb;}^2}---\left(12\right)$

其中，y₀(i)-实测值；y_c(i)-预测值；-实测值的均值；n-资料序列长度；

3.相关系数(correlation coefficient)

$r=\frac{\underset{i=1}{\overset{n}{\&Sigma;}}(p_i-\stackrel{\&OverBar;}{p})(q_i-\stackrel{\&OverBar;}{q})}{\sqrt{\underset{i=1}{\overset{n}{\&Sigma;}}{(p_i-\stackrel{\&OverBar;}{p})}^2}\sqrt{\underset{i=1}{\overset{n}{\&Sigma;}}{(q_i-\stackrel{\&OverBar;}{q})}^2}}---\left(13\right)$

上式中，p_i、q_i分别代表预测值和实际值，分别是实际值和预测值的统计平均值。

确定多模型的输入

首先采用DDTW估计上下游站点之间的流量传播时间，然后进行聚类。聚类方法采用基于DBI的K-means方法来划分模型，本实例将最小类设置为5，最大类设置为22，DBI指标值随类别数变化的波动情况如图3所示，以间隔为单位的时间聚类中心结果如表1所示。

从图3可以看到，最小DBI对应的类别数目为7，因此本实例中传播时间的最佳聚类数目可以确定为7。比较各类的聚类中心如表1所示：

表1基于DBI的K均值聚类中心

分析表1的聚类中心结果，各类的传播时间中心都代表了各类别主要的传播时间，且有显著性区别，从而为后面的子模型确定不同的输入信息。

得到聚类结果后，对每个模型各个站点的传播时间作直方图分布分析。息县、潢川、班台三个站点的传播时间(传播间隔)直方图分别如图4(a)～图4(c)所示。息县、潢川、班台各站点在范围[18,48]范围内传播时间所占的比例如表2所示：

表2类1中息县、潢川、班台[18,48]范围内传播时间的比例

	xx	Hc	bt
				18h	0	0	0
24h	0	1.55％	0
				30h	0	7.77％	0
36h	0	12.44％	29.53％
				42h	20.21％	15.03％	16.58％
48h	79.79％	63.21％	53.89％

表2中，比例为0的代表该类的传播时间中没有出现过此种传播时间，选取涵盖该类中各站点90％信息的传播时间作为输入信息。因此，对于类1，根据直方图分布确定的上游三个站点的输入信息如表3所示，其中xx-42h表示息县站42小时前的流量值，hc-36h表示潢川站36小时前的流量值，bt-36h表示班台站36小时前的流量值，其他依此类推。据此原理，同样处理其他各类信息。

表3类1中息县、潢川、班台的输入信息

站点	输入信息
		息县xx	xx-42h、xx-48h
潢川hc	hc-36h、hc-42h、hc-48h
		班台bt	bt-36h、bt-42h、bt-48h

由于雨量信息是个复杂因素，这里采用相关性分析法来确定每个子模型中雨量信息的输入。分别计算每类中王家坝目标流量与前2-3天累加雨量信息(pa++2-3d，即(pa-24h)+...+(pa-72h))、前3-4天累加雨量信息(pa++3-4d，即(pa-72h)+...+(pa-96h))、pa-24h、pa-30h、pa-36h、pa-42h、pa-48h、pa-54h、pa-60h、pa-66h、pa-72h的相关系数。本实例选取相关系数最大的两个时刻作为雨量的输入信息。

按照上面介绍的确定上游各站点输入信息及每类雨量输入信息的方式确定最终的各子模型输入信息，因此本实例所确定的最终输入输出信息如表4所示,其中wjb-0h表示王家坝站当前时刻的流量，wjb-24h表示王家坝站24小时前的流量，pa++2-3d为前2-3天累加雨量信息，即(pa-24h)+...+(pa-72h))，pa-24h为24小时前降雨量，其他依此类推。

表4基于传播时间分析的各模型输入输出信息

对每类训练集进行SVM回归建模，并用建得的模型对对应的测试点进行模拟预测，其预测结果如表5所示。

表5本发明方法

从表5可以看出，即使是训练数据较少的子模型，其测试集的预测效果也能保持与实际值相吻合，如模型1的训练数据和对应的测试数目都是所有子模型中最少的，但其表现的预测效果却是所有子模型中最好的，预测精度高达0.9967，可见，在本发明的方法中，模型的预测效果并没有受到样本个数的影响。

为了说明本发明方法相对于现有技术的优势，比较了现有基于流量聚类的多模型预测方法(参见文献[尹雄锐，张翔，夏军.基于聚类分析的人工神经网络洪水预报模型研究，四川大学学报(工程科学版)，2007，39(3)：34-40]、文献[周晓阳，张勇传.洪水的分类预测及优化调度.水科学进展，1997(6)：8(2)]与本发明方法的区别。为了便于比较，同样采用7-模型的流量模型对比。根据已有资料统计，一般情况下，息县、潢川、班台流量、区间面降雨汇至下游王家坝站的传播时间分别为30小时、30小时、24小时、60小时。根据上述条件，使用wjb-24h、xx-30h、hc-30h、bt-24h的流量数据参与聚类。聚类完成后，再确定每个子模型所需要的输入信息。保持上游各站点的输入信息不变，通过分别计算每类中王家坝目标流量与前2-3天累加雨量信息(pa++2-3d)、前3-4天累加雨量信息(pa++3-4d)的相关系数选取相关性最高的作为雨量的输入信息。最终确定的基于流量聚类的各子模型输入输出信息如表6所示。

表6基于流量聚类的各子模型的输入输出信息

同时，建立常规条件下的基于所有流量数据的单一模型。具体的输入信息为：wjb-24h、xx-30h、hc-30h、bt-24h、pa-60h；输出信息为：wjb-0h。

根据以上信息分别建立基于流量聚类的多模型与常规条件下的单一模型，并与本发明方法相比较，即基于所有流量数据的单模型、基于汛期数据的流量聚类的多模型预测以及基于传播时间聚类分析的多模型集成洪水预报方法(本发明方法)进行对比，各方法预测结果及曲线对比如表7所示：

表7单模型与多模型整体预测结果对比

对比发现：单模型的方法在某些峰值点上会异常的高出，而基于流量聚类的多模型方法在这些峰值点的预测误差有了明显的降低，但是对特大洪峰的预测误差仍然较大，没有得到很好的改善；基于传播时间的多模型方法又要优于传统的基于流量聚类的多模型预测，提高了特大洪峰点的峰值预测精度。综合各项指标，多模型综合预测结果要优于所有流量数据的单模型预测结果，并且基于传播时间聚类的多模型综合预测方法在高流量级别上较其他两个模型有了明显改善。

实验证明本发明方法简单有效，不仅可以准确估算出各模型的洪峰传播时间，还有利于帮助准确确定各模型输入。最终的预测精度比较理想。

Claims (8)

1.基于传播时间聚类分析的多模型集成洪水预报方法，首先根据历史流量/水位数据建立包含多个子模型的洪水预报综合模型，然后利用所建立的洪水预报综合模型进行洪水预报；其特征在于，具体包括以下步骤：

步骤1、对于目标站点的历史流量/水位数据中的每个流量/水位数据，利用各上游站点与目标站点之间的流量传播时间构造对应于该流量/水位数据的传播时间向量，目标站点的历史流量/水位数据中的所有流量/水位数据所对应的传播时间向量构成传播时间向量集合；

步骤2、对所述传播时间向量集合进行聚类，得到k个簇，并对聚类得到的每个簇分别进行以下处理：

对于每一个上游站点，统计该簇中所包含的该上游站点与目标站点之间的流量传播时间的分布情况，并将出现频次最低的部分流量传播时间剔除，剩余的流量传播时间作为该簇中该上游站点的可用流量传播时间，最终得到该簇中各上游站点的可用流量传播时间；

步骤3、构建k个训练样本集，k个训练样本集与步骤2得到的k个簇一一对应，k为聚类得到的类别数；其中任意一个训练样本集中的训练样本按照以下方法得到：

选取一组目标站点的历史流量/水位数据分别作为该训练样本集中各训练样本的输出；对每一个训练样本的输出，根据其所属训练样本集所对应的簇中各上游站点的可用流量传播时间，从各上游站点的历史流量/水位数据中确定相应的流量/水位数据，并结合相应的雨量输入信息及目标站点的预见期前的历史流量/水位数据，构成该训练样本的输入；

步骤4、利用所构建的k个训练样本集各自对预测模型进行训练，得到k个预测子模型，这k个预测子模型共同构成洪水预报综合模型；

步骤5、根据测试样本按照以下方法从洪水预报综合模型中选择相应的预测子模型：对测试样本中目标站点的流量/水位数据，利用各上游站点与目标站点之间的流量传播时间构造对应于该流量/水位数据的传播时间向量，并从步骤2中聚类得到的各个簇中选出簇中心与该传播时间向量的距离最小的簇，距离最小的簇所对应的预测子模型即为所选择的预测子模型；

步骤6、以所述测试样本作为输入，利用所选择的预测子模型进行目标站点的洪水预报。

2.如权利要求1所述基于传播时间聚类分析的多模型集成洪水预报方法，其特征在于，所述各上游站点与目标站点之间的流量传播时间利用派生动态时间弯曲算法获取。

3.如权利要求2所述基于传播时间聚类分析的多模型集成洪水预报方法，其特征在于，在利用派生动态时间弯曲算法获取各上游站点与目标站点之间的流量传播时间时，针对上游站点流量/水位的特征点，若该特征点是极小点，则将该特征点所匹配的目标站点若干流量/水位数据点中出现极小值的那一个作为上游站点流量/水位的特征点的唯一匹配；若是极大点，则将该特征点所匹配的目标站点若干流量/水位数据点中出现极大值的那一个作为上游站点流量/水位的特征点的唯一匹配；对于非极值点出现多个匹配情况，则采用与该点之前最近的极值点传播时间最相近的传播时间。

4.如权利要求1所述基于传播时间聚类分析的多模型集成洪水预报方法，其特征在于，所述聚类使用基于DBI指标的K-means聚类方法。

5.如权利要求1所述基于传播时间聚类分析的多模型集成洪水预报方法，其特征在于，所述预测模型为支持向量机回归模型。

6.如权利要求1所述基于传播时间聚类分析的多模型集成洪水预报方法，其特征在于，步骤2中所述将出现频次最低的部分流量传播时间剔除，具体是指将每个上游站点所对应的流量传播时间中出现频次最低的10％的流量传播时间剔除。

7.如权利要求1所述基于传播时间聚类分析的多模型集成洪水预报方法，其特征在于，所述簇中心与传播时间向量的距离为欧氏距离。

8.如权利要求1所述基于传播时间聚类分析的多模型集成洪水预报方法，其特征在于，所述训练样本中的雨量输入信息利用相关系数分析方法确定。