技术领域
[0001] 本
发明属于水利工程技术领域，尤其涉及防洪预报技术领域，具体为一种基于时间序列聚类的流域洪水响应相似性分析方法。
背景技术
[0002] 目前，我国大江大河及其主要支流已经形成以堤防、水库和蓄滞洪区为主的防洪工程体系，防汛预警预报系统等非工程措施也逐步得到加强，基本能防御主要江河常遇洪水。然而对于5万多条中小河流，其分布广、数量多，自然地理、
气候条件复杂多样，防洪能
力总体落后，特别是近年来极端天气事件增多，中小流域暴雨事件增加，常常造成严重的洪涝灾害。中小河流洪水灾害造成的损失已成为我国洪涝灾害损失的主体。
[0003] 洪水预报是非工程措施的重要组成部分，能够有效提高流域和区域的防灾减灾能力。但对于缺少径流资料的中小河流，传统的经验方法往往不再适用，目前广为采取的洪水预报方法是依据相似性和区域化的预报模型参数移植方法。这一方法非常依赖对水文相似性的准确判断，然而目前针对水文相似性判断的水文特征指标、水文模型参数、统计指标等方法存在片面化、不确定性强、针对性差等问题。例如水文特征指标法中经常采用的径流系数、流量历时曲线、季节指数、基流指数等，通常只能片面的描述流域的水文特征，且这种相似性的判断往往更适合于中长期径流过程而非针对短时段洪水响应过程。通过水文模型参数的相似性来定义和判断水文相似性，这种定义方式的潜在假定是率定后的模型参数相似反映了降雨径流转换行为的相似。但是这种相似受到模型结构的影响，可能会存在一种模型结构的定义下是相似的，在其他模型结构的定义下却不相似的情况，而且概念性模型的参数不确定性以及异参同效问题也对这种定义的可靠性带来隐患。通过统计指标定义相似，这一类方法通常关注洪峰等洪水行为的一部分特征，如以无量纲的洪水
频率曲线为水文相似的衡量标准，也存在片面性的问题。
发明内容
[0004] 本发明的目的在于克服以上技术
缺陷，提出一种基于时间序列聚类的流域洪水响应相似性分析方法，区别于以往利用单一指数或者模型参数对水文相似进行判断。本发明的目的是通过以下技术方案实现的：
[0005] 一种基于时间序列聚类的流域洪水响应相似性分析方法，包括以下步骤：
[0006] 1)数据的收集、处理与保存：收集待分析流域出口水文
站点的径流数据，获得流域出口流量时间序列，进行等时段处理获得等时段流域出口流量时间序列数据保存至
数据库中，并使用流域编码进行标识；
[0007] 2)洪水场次自动划分：根据流域编码依次读取数据库中各流域出口流量时间序列，以洪水的起涨和消落过程作为洪水事件的划分标准将其划分为独立的场次洪水；
[0008] 3)生成洪水事件样本集合：基于步骤2)中提取的场次洪水生成洪水事件样本集合{S1，S2，...，Si，...，Sn}，其中包含n个子集，分别记录n个不同流域的场次洪水数据，子集Si为流域i的场次洪水样本集合，其元素为基于流域i的流量数据提取的且经过归一化处理的场次洪水，将归一化的场次洪水以时间序列的形式进行保存，同时建立子集及洪水事件索引；
[0009] 4)基于洪水事件样本生成聚类树：以步骤3)集合中的子集为单元进行洪水事件层次聚类，生成n棵聚类树，每棵聚类树对应一个流域；
[0010] 5)基于聚类树提取各流域代表性洪水：依次对n个流域的聚类树进行分析，根据样本总数量大小设置参数N，搜索每棵聚类树第1层至第N层的各个
节点，计算各节点的聚类中心，作为该流域的代表性洪水，利用提取的所有代表性洪水生成训练集，并建立流域索引；
[0011] 6)基于代表性洪水进行流域洪水响应相似性分析：利用步骤5)代表性洪水生成的训练集计算新的距离矩阵MatrixB，矩阵大小为(n×n)，n为流域个数，矩阵元素(i，j)为流域i与流域j的洪水响应相似性度量，矩阵元素(i，j)值d(i，j)的计算方法如下：
[0012]
[0013] 式中：N为步骤5)中设置的层数；di′为流域i与流域j在第i′层的代表性洪水的各种对应方式下DTW距离求和结果中最小的距离和值；
[0014] 式(1)di＇具体的计算方法为：
[0015] di′＝min{sum(DDTW(Cim，Cjn))；m，n＝1，2，...，i′}，其中流域i与流域j在第i′层的代表性洪水各有i′个，共有Amn种对应方式， 其中n＝m＝i′，计算各种对应方式代表性洪水的DTW距离求和，di′为其中最小的DTW距离和值；
Cim与Cjn为流域i与流域j在第i′层的代表性洪水，m、n为第i′层中各代表性洪水的索引；N为步骤5)中设置的层数。
[0016] 完成MatrixB的计算后，依次搜索矩阵中非主对
角线元素的最小值，最小值所在行列号对应的两个流域即为洪水响应最为相似的流域。
[0017] 进一步的，步骤2)中采用
算法从连续时间序列中自动提取洪水事件，具体方法为：
[0018] 2-1基流分割，将径流序列中的基流部分和洪水部分划分开，通过数字滤波法得到流量过程的基流部分，如式(2)所示：
[0019]
[0020] 式中，bt为时刻t的基流，Qt为时刻t的径流，β为滤波系数；
[0021] 通过滤波次数和滤波系数来控制滤波效果；获得基流序列后，通过式(3)获得洪水序列{q1，q2，q3，...，qt}，
[0022] qt＝Qt-bt (3)
[0023] 2-2对洪水序列进行平滑处理，消除噪声项和异常点的影响，平滑滤波公式如式(4)：
[0024]
[0025] 2-3识别序列中的转折点：计算洪水序列的一阶差分序列，根据一阶差分序列的正负变换判断序列转折点的
位置，并对极大值与极小值加以区分；对序列首尾的处理：若首尾值为极大值，则将其去掉；设定
阈值Thmin，若首尾值小于序列均值除以阈值Thmin，则设定其为极小值；记录转折点序列及各值对应的极大值、极小值；
[0026] 2-4识别洪水事件的开始、结束点：设定阈值Thslp，
选定第一个极小值作为开始点并向后搜索另一极小值，依据转折点数组的一阶差分序列{d1，d2，d3，...，di，...，dt}进行判断，寻找满足式(5)要求的极小值点Mini作为结束点：
[0027] Mini-Min1＜Thslp·max(|d1|，|d2|，|d3|，...，|di|) (5)
[0028] 2-5对步骤2-4中提取的洪水事件进行筛选与处理，具体步骤为：a.设置阈值Thpeak，对于一次独立的洪水事件，若峰值与序列起始点或结束点的差值小于阈值Thpeak，则认为本次洪水过程量级不足以纳入考虑范围；b.设定动态坡度阈值Thdy，对于一次独立的洪水事件，动态坡度为动态坡度阈值Thdy与极差的乘积，依据动态坡度删除洪水事件前部与后部的平坦部分；c.设定时间阈值ThΔT，对于一次独立的洪水事件，如持续时间小于阈值ThΔT，则认为本次洪水过程不足以纳入考虑范围。
[0029] 进一步的，步骤3)中归一化的方法采用缩放法，如下式进行：
[0030]
[0031] 其中，xmax为时间序列中的最大值，xmin为时间序列中的最小值，xi与x′i分别为缩放前后的序列数值。
[0032] 进一步的，步骤4)中聚类树生成过程中时间序列
聚类分析时具体步骤为：
[0033] 4-1.生成初始簇：将子集中的每一个元素作为一个初始簇；
[0034] 4-2.计算基于一个子集的距离矩阵：矩阵大小为(m×m)，m为该子集中包含的洪水事件个数，矩阵的元素(i，j)为i簇与j簇的相似度，表示洪水事件i与洪水事件j的相似度，使用DTW距离作为相似性度量标准，距离越小则相似性越强；
[0035] 4-3.基于步骤4-2中的距离矩阵对簇进行合并，找出距离最近的两个簇且进行合并，将聚类簇重新编号，并计算新簇与其他各簇的距离，更新距离矩阵；
[0036] 4-4.重复步骤4-3直至所有的聚类簇合并为一个簇，由此生成一棵聚类树；
[0037] 4-5.重复步骤4-2～4-4使基于样本集合中每个子集均生成一棵对应的聚类树。
[0038] 进一步的，DTW距离的计算方法为：
[0039] 对时间序列X＝{x1，x2，...，xi，...，xm}和Y＝{y1，y2，...，yi，...，yn}，通过扭曲路径W来表示时间序列X与Y间的映射关系，W＝{w1，w2，...，wk，...，wK}，max(n，m)≤k≤n+m-1，其中：m、n为时间序列X和时间序列Y的长度，K为扭曲距离的长度；W的第k个元素记为wk＝(i，j)，表示时间序列X的第i个元素与时间序列Y的第j个元素的对应关系；构建一个m×n阶矩阵，矩阵元素(i，j)为两个时间序列点xi和点yj之间的距离d(xi，yj)＝(xi-yj)2，定义点(i，j)的累积距离计算公式为：
[0040] γ(i，j)＝d(xi，yj)+min{γ(i-1，j-1)，γ(i-1，j)，γ(i，j-1)} (7)[0041] 其中，i＝1，...，m；j＝1，...，n；为矩阵元素索引；给定初始条件γ(1，1)＝d(x1，y1)，以
迭代计算得到累积距离矩阵， 即为时间序列X与Y的DTW距离。
[0042] 进一步的，步骤5)中聚类中心的计算方法为：生成距离矩阵MatrixD，矩阵大小为(mi×mi)，mi为节点中包含的洪水事件个数，矩阵的元素(i，j)为洪水事件i与洪水事件j的DTW距离，首先计算距离矩阵MatrixD，然后计算其各行的和值，和值最小的行索引所对应的洪水事件即为聚类中心。
[0043] 本发明的有益效果：
[0044] 本发明提出一种基于时间序列聚类的流域洪水响应相似性分析方法，区别于以往利用单一指数或者模型参数对水文相似进行判断，本发明中提出的方法能够充分利用流量时间序列数据所包含的信息，针
对流域洪水响应的相似性进行判断，基于该方法分析的相似性能够有效提高洪水预报的准确性。
[0045] 下面结合
附图及具体实施方式对本发明作进一步详细说明。
附图说明
[0047] 图2为洪水数据插值示意图；
[0048] 图3为洪水场次提取示意图；
[0049] 图4为时间序列的动态扭曲路径；
[0050] 图5为聚类树示意图；
[0051] 图6为聚类中心示意图；
[0052] 图7为子集1的洪水样本情况；
[0053] 图8为子集2的洪水样本情况；
[0054] 图9为子集3的洪水样本情况；
[0055] 图10为子集1聚类树；
[0056] 图11为子集1代表性洪水；
[0057] 图12为流域洪水响应相似性矩阵。
具体实施方式
[0059] 1)数据的收集、处理与保存
[0060] 收集各待分析流域出口水文站点的洪水数据，根据
机器学习对数据量的需求，洪水数据需要
覆盖10年或10年以上。
[0061] 将原始
数据处理为等时段时间序列，若原始数据为非等时段数据，则需对数据进行插值处理，对于洪水数据建议采取线性内插，如图2所示，利用原始序列{Q1，Q2，Q3，...，Q7}插值获得等时段流量时间序列{Q′1，Q′2，Q′3，...，Q′12}。
[0062] 将处理好的等时段流域出口流量时间序列数据保存至数据库，并使用流域编码进行标识。
[0063] 2)场次洪水自动提取
[0064] 根据流域编码依次读取数据库中各流域的出口流量时间序列，以洪水的起涨和消落过程作为洪水事件的划分标准将其划分为独立的场次洪水。由于机器学习对数据量的要求较大，依靠人工划分效率较低，因而采用算法从连续时间序列中自动提取洪水事件。具体方法为：
[0065] 2-1.基流分割，将径流序列中的基流部分和洪水部分划分开，通过数字滤波法得到流量过程的基流部分，如下式所示。
[0066]
[0067] 其中bt为时刻t的基流，Qt为时刻t的径流，β为滤波系数，通过滤波次数和滤波系数来控制滤波效果。获得基流序列后，通过下式获得洪水序列{q1，q2，q3，...，qt}。
[0068] qt＝Qt-bt (3)
[0069] 2-2.对洪水序列进行平滑处理，消除噪声项和异常点的影响，平滑滤波公式如下：
[0070]
[0071] 2-3.识别序列中的转折点，具体步骤为：计算洪水序列的一阶差分序列，根据一阶差分序列的正负变换判断序列转折点的位置，并对极大值与极小值加以区分。对序列首尾的处理：若首尾值为极大值，则将其去掉；设定阈值Thmin，若首尾值小于序列均值除以阈值Thmin，则设定其为极小值。记录转折点序列及各值对应的峰(极大值)、谷(极小值)标记。
[0072] 2-4.识别洪水事件的开始、结束点，具体步骤为：设定阈值Thslp，选定第一个极小值作为开始点并向后搜索另一极小值，依据转折点数组的一阶差分序列进行判断{d1，d2，d3，...，di，...，dt}，寻找满足下式要求的极小值点Mini作为结束点：
[0073] Mini-Min1＜Thslp·max(|d1|，|d2|，|d3|，...，|di|) (5)
[0074] 2-5.对步骤2-4中提取的洪水事件进行筛选与处理，具体步骤为：a.设置阈值Thpeak，对于一次独立的洪水事件，若峰值与序列起始点或结束点的差值小于阈值Thpeak，则认为本次洪水过程量级不足以纳入考虑范围；b.设定动态坡度阈值Thdy，对于一次独立的洪水事件，动态坡度为阈值Thdy与极差的乘积，依据动态坡度删除洪水事件前部与后部的平坦部分；c.设定时间阈值ThΔT，对于一次独立的洪水事件，如持续时间小于阈值ThΔT，则认为本次洪水过程不足以纳入考虑范围。如图3所示，横纵坐标分别代表时间与流量，Q2-Q1大于阈值Thslp·max(|d1|，|d2|，|d3|，...，|di|)，则不作为结束点，而Q3-Q1小于阈值Thslp·max(|d1|，|d2|，|d3|，...，|di|)且T3-T1大于时间阈值ThΔT，则认为是一次独立的洪水事件。
[0075] 依据上述方法，得到n′个场次洪水序列{Qi1，Qi2，...，Qik，}及其时间标识序列{Ti1，Ti2，...，Tik，}，其中i＝1，...，n′，n′为洪水场次个数，k′为该场洪水对应的时段个数。
[0076] 3)生成洪水事件样本集合
[0077] 基于步骤2)中提取的场次洪水时间序列生成洪水事件样本集合{S1，S2，...，Si，...，Sn}，其中包含n个子集，分别记录n个不同流域的场次洪水数据。设子集Si为流域i的场次洪水样本集合，其元素为步骤2)中基于流域i的流量数据提取的且经过归一化处理n′场洪水，归一化的方法可以采用缩放法，如下式所示：
[0078]
[0079] 其中，xmax为时间序列中的最大值，xmin为时间序列中的最小值，xi与x′i分别为缩放前后的序列数值。
[0080] 将归一化的场次洪水以时间序列的形式进行保存，同时建立子集及洪水事件索引。
[0081] 4)基于洪水事件样本生成聚类树
[0082] 以步骤3)集合中的子集为单元进行洪水事件层次聚类，生成n棵聚类树。基于索引遍历各子集，针对单个子集的洪水数据，时间序列聚类分析的具体步骤为：
[0083] 4-1.生成初始簇，将子集中的每一个元素作为一个初始簇，对于一个具有m个元素的集合D＝{x1，x2，...，xm}，设定初始簇集合C＝{C1，C2，...，Cm}，其中Cj＝{xj}；
[0084] 4-2.计算第一距离矩阵MatrixF，矩阵大小为(m×m)，m为子集中包含的洪水事件个数，矩阵的元素(i，j)为洪水事件i与洪水事件j的相似度，因而主对角线元素为0且为对称矩阵。使用DTW距离作为相似性度量标准，距离越小则相似性越强，DTW距离计算方法如下：
[0085] 对时间序列X＝{x1，x2，...，xi，...，xm}和Y＝{y1，y2，...，yi，...，yn}，通过扭曲路径W来表示时间序列X与Y间的映射关系，如图4所示，W＝{w1，w2，...，wk，...，wK}，max(n，m)≤K≤n+m-1，W的第k个元素记为wk＝(i，j)，表示时间序列X的第i个元素与时间序列Y的第j个元素的对应关系。扭曲路径的选取有三个约束条件：扭曲路径始于矩阵的起始元素，结束于对角元素，即w1＝(1，1)，wK＝(m，n)；扭曲路径每一步都是连续的，即对于wk＝(a，b)，wk-1＝(a′，b′)，要求a-a′≤1且b-b′≤1；扭曲路径在时间轴上是单调的，即对于wk＝(a，b)，wk-1＝(a′，b′)，要求a-a′≥0且b-b′≥0。
[0086] 能够满足约束条件的路径有很多条，此处寻找扭曲代价最小的路径，即：
[0087]
[0088] 其中d(wk)为wk代表的两个对应元素间的距离。
[0089] 根据动态规划思想，若点(i，j)在最佳路径上，那么从点(1，1)到点(i，j)的子路径也是局部最优解，即从点(1，1)到点(m，n)的最佳路径可以由起始点(1，1)到终点(m，n)之间的局部最优解递归搜索获得，因而可以方便地找到这个最佳路径。具体步骤为：首先构建一个m×n阶矩阵，矩阵元素(i，j)为两个时间序列点xi和点yj之间的距离d(xi，yj)＝(xi-yj)2。定义点(i，j)的累积距离计算公式：
[0090] γ(i，j)＝d(xi，yj)+min{γ(i-1，j-1)，γ(i-1，j)，γ(i，j-1)} (7)[0091] 给定初始条件γ(1，1)＝d(x1，y1)，可以迭代计算得到累积距离矩阵。即为时间序列X与Y的DTW距离，从点γ(m，n)出发反向搜索累积距离矩阵即可得到最佳匹配路径。
[0092] 4-3.对簇进行合并，找出距离最近的两个簇Ci*和Cj*，合并Ci*和Cj*：Ci*＝Ci*∪Cj*，将聚类簇重新编号，删除当前距离矩阵的第j*行和第j*列，并计算新簇与其他各簇的距离，更新距离矩阵；
[0093] 4-4.重复上一步骤直至所有聚类簇合并为一个簇，由此生成一棵聚类树，如图5所示。
[0094] 4-5.重复步骤4-2～4-4使基于样本集合中每个子集均生成一棵对应的聚类树。
[0095] 5)基于聚类树提取各流域代表性洪水
[0096] 依次对n个流域的聚类树进行分析，提取聚类中心。对于单个聚类树，树的各节点均代表一类洪水事件，而每个节点的聚类中心即为此节点最具代表性的洪水事件，如图6所示。将根节点作为聚类树的第1层，则聚类树的第n层包含n个节点及n个聚类中心。将节点样本中与其他各元素的距离和值最小的元素作为节点的聚类中心，则对于节点i，其聚类中心的计算方法为：
[0097] 生成距离矩阵MatrixD，矩阵大小为(mi×mi)，mi为节点中包含的洪水事件个数，矩阵的元素(i，j)为洪水事件i与洪水事件j的DTW距离。首先计算距离矩阵MatrixD，然后计算其各行的和值，和值最小的行索引所对应的洪水事件即为聚类中心。
[0098] 根据样本总数量的大小设定参数N，计算第1层至第N层各节点的聚类中心，作为该流域的代表性洪水。利用提取的所有代表性洪水生成训练集，并建立流域索引。
[0099] 6)基于代表性洪水进行流域洪水响应相似性分析
[0100] 利用步骤5)中生成的代表性洪水训练集计算新的距离矩阵MatrixB，矩阵大小为(n×n)，n为流域个数，矩阵元素(i，j)为流域i与流域j的洪水响应相似性度量，因而矩阵为对称矩阵且主对角线元素为0。矩阵元素(i，j)值d(i，j)的计算方法如下：
[0101]
[0102] 式中：di′＝min{sum(DDTW(Cim，Cjn))；m，n＝1，2，...，i′}，其中流域i与流域j在第i′层的代表性洪水各有i′个，共有Amn种对应方式，其中n＝m＝i′，计算各种对应方式代表性洪水的DTW距离求和，di′为其中最小的DTW距离和值；Cim与Cjn为流域i与流域j在第i′层的代表性洪水，m、n为第i′层中各代表性洪水的索引；
N为步骤5)中设置的层数；
[0103] 完成MatrixB的计算后，依次搜索矩阵中非主对角线元素的最小值，最小值所在行列号对应的两个流域即为洪水响应最为相似的流域。
[0104] 收集到我国黄河流域中游49个子流域出口水文站点的历史洪水数据，数据起止时间如下表所示，数据年限均在10年以上。
[0105] 表1流域资料情况表
[0106]
[0107]
[0108] 经过插值处理为等时段数据后，进行场次洪水的自动提取，各流域提取的洪水场次情况如下表所示：
[0109] 表2流域场次洪水提取情况
[0110]
[0111]
[0112] 对各场次洪水进行标准化处理，以样本子集1、2、3为例，三个子集的洪水样本情况如图7～9所示：
[0113] 以各洪水样本子集为单元进行洪水事件层次聚类，共生成49棵聚类树，以第1子集的聚类树为例，如图10所示。
[0114] 根据生成的聚类树，设定参数N＝2，提取各簇的聚类中心，即流域的代表性洪水，以第1个样本子集的代表性洪水为例，如图11所示：
[0115] 根据代表性洪水生成大小为49×49的距离矩阵MatrixB，即为流域洪水响应相似性矩阵，如图12所示。
[0116] 查找矩阵元素中的最小值所在行列号，得到流域3与流域5的洪水响应最为相似。
