CN111080107B - 一种基于时间序列聚类的流域洪水响应相似性分析方法 - Google Patents

一种基于时间序列聚类的流域洪水响应相似性分析方法 Download PDF

Info

Publication number
CN111080107B
CN111080107B CN201911242133.9A CN201911242133A CN111080107B CN 111080107 B CN111080107 B CN 111080107B CN 201911242133 A CN201911242133 A CN 201911242133A CN 111080107 B CN111080107 B CN 111080107B
Authority
CN
China
Prior art keywords
flood
sequence
matrix
clustering
basin
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201911242133.9A
Other languages
English (en)
Other versions
CN111080107A (zh
Inventor
王帆
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China Institute of Water Resources and Hydropower Research
Original Assignee
China Institute of Water Resources and Hydropower Research
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China Institute of Water Resources and Hydropower Research filed Critical China Institute of Water Resources and Hydropower Research
Priority to CN201911242133.9A priority Critical patent/CN111080107B/zh
Publication of CN111080107A publication Critical patent/CN111080107A/zh
Application granted granted Critical
Publication of CN111080107B publication Critical patent/CN111080107B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/06—Resources, workflows, human or project management; Enterprise or organisation planning; Enterprise or organisation modelling
    • G06Q10/063—Operations research, analysis or management
    • G06Q10/0639—Performance analysis of employees; Performance analysis of enterprise or organisation operations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
    • G06F16/22—Indexing; Data structures therefor; Storage structures
    • G06F16/2228—Indexing structures
    • G06F16/2246—Trees, e.g. B+trees
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
    • G06F16/24—Querying
    • G06F16/245—Query processing
    • G06F16/2458—Special types of queries, e.g. statistical queries, fuzzy queries or distributed queries
    • G06F16/2465—Query processing support for facilitating data mining operations in structured databases
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
    • G06F16/24—Querying
    • G06F16/245—Query processing
    • G06F16/2458—Special types of queries, e.g. statistical queries, fuzzy queries or distributed queries
    • G06F16/2474—Sequence data queries, e.g. querying versioned data
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
    • G06F16/29—Geographical information databases
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q50/00—Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
    • G06Q50/10—Services
    • G06Q50/26—Government or public services
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A10/00—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE at coastal zones; at river basins
    • Y02A10/40—Controlling or monitoring, e.g. of flood or hurricane; Forecasting, e.g. risk assessment or mapping

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Business, Economics & Management (AREA)
  • Physics & Mathematics (AREA)
  • Databases & Information Systems (AREA)
  • General Physics & Mathematics (AREA)
  • Human Resources & Organizations (AREA)
  • General Engineering & Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Educational Administration (AREA)
  • Strategic Management (AREA)
  • Economics (AREA)
  • Software Systems (AREA)
  • Development Economics (AREA)
  • Tourism & Hospitality (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Marketing (AREA)
  • Fuzzy Systems (AREA)
  • Mathematical Physics (AREA)
  • Probability & Statistics with Applications (AREA)
  • General Business, Economics & Management (AREA)
  • Computational Linguistics (AREA)
  • Primary Health Care (AREA)
  • Health & Medical Sciences (AREA)
  • Remote Sensing (AREA)
  • General Health & Medical Sciences (AREA)
  • Game Theory and Decision Science (AREA)
  • Operations Research (AREA)
  • Quality & Reliability (AREA)
  • Alarm Systems (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种基于时间序列聚类的流域洪水响应相似性分析方法，包括以下步骤：1）数据的收集、处理与保存；2）洪水场次自动划分；3）生成洪水事件样本集合；4）基于洪水事件样本生成聚类树；5）基于聚类树提取各流域代表性洪水；6）基于代表性洪水进行流域洪水响应相似性分析。本发明中提出的方法能够充分利用流量时间序列数据所包含的信息，针对流域洪水响应的相似性进行判断，基于该方法分析的相似性能够有效提高洪水预报的准确性。

Description

一种基于时间序列聚类的流域洪水响应相似性分析方法
技术领域
本发明属于水利工程技术领域，尤其涉及防洪预报技术领域，具体为一种基于时间序列聚类的流域洪水响应相似性分析方法。
背景技术
目前，我国大江大河及其主要支流已经形成以堤防、水库和蓄滞洪区为主的防洪工程体系，防汛预警预报系统等非工程措施也逐步得到加强，基本能防御主要江河常遇洪水。然而对于5万多条中小河流，其分布广、数量多，自然地理、气候条件复杂多样，防洪能力总体落后，特别是近年来极端天气事件增多，中小流域暴雨事件增加，常常造成严重的洪涝灾害。中小河流洪水灾害造成的损失已成为我国洪涝灾害损失的主体。
洪水预报是非工程措施的重要组成部分，能够有效提高流域和区域的防灾减灾能力。但对于缺少径流资料的中小河流，传统的经验方法往往不再适用，目前广为采取的洪水预报方法是依据相似性和区域化的预报模型参数移植方法。这一方法非常依赖对水文相似性的准确判断，然而目前针对水文相似性判断的水文特征指标、水文模型参数、统计指标等方法存在片面化、不确定性强、针对性差等问题。例如水文特征指标法中经常采用的径流系数、流量历时曲线、季节指数、基流指数等，通常只能片面的描述流域的水文特征，且这种相似性的判断往往更适合于中长期径流过程而非针对短时段洪水响应过程。通过水文模型参数的相似性来定义和判断水文相似性，这种定义方式的潜在假定是率定后的模型参数相似反映了降雨径流转换行为的相似。但是这种相似受到模型结构的影响，可能会存在一种模型结构的定义下是相似的，在其他模型结构的定义下却不相似的情况，而且概念性模型的参数不确定性以及异参同效问题也对这种定义的可靠性带来隐患。通过统计指标定义相似，这一类方法通常关注洪峰等洪水行为的一部分特征，如以无量纲的洪水频率曲线为水文相似的衡量标准，也存在片面性的问题。
发明内容
本发明的目的在于克服以上技术缺陷，提出一种基于时间序列聚类的流域洪水响应相似性分析方法，区别于以往利用单一指数或者模型参数对水文相似进行判断。本发明的目的是通过以下技术方案实现的：
一种基于时间序列聚类的流域洪水响应相似性分析方法，包括以下步骤：
1)数据的收集、处理与保存：收集待分析流域出口水文站点的径流数据，获得流域出口流量时间序列，进行等时段处理获得等时段流域出口流量时间序列数据保存至数据库中，并使用流域编码进行标识；
2)洪水场次自动划分：根据流域编码依次读取数据库中各流域出口流量时间序列，以洪水的起涨和消落过程作为洪水事件的划分标准将其划分为独立的场次洪水；
3)生成洪水事件样本集合：基于步骤2)中提取的场次洪水生成洪水事件样本集合{S1，S2，...，Si，...，Sn}，其中包含n个子集，分别记录n个不同流域的场次洪水数据，子集Si为流域i的场次洪水样本集合，其元素为基于流域i的流量数据提取的且经过归一化处理的场次洪水，将归一化的场次洪水以时间序列的形式进行保存，同时建立子集及洪水事件索引；
4)基于洪水事件样本生成聚类树：以步骤3)集合中的子集为单元进行洪水事件层次聚类，生成n棵聚类树，每棵聚类树对应一个流域；
5)基于聚类树提取各流域代表性洪水：依次对n个流域的聚类树进行分析，根据样本总数量大小设置参数N，搜索每棵聚类树第1层至第N层的各个节点，计算各节点的聚类中心，作为该流域的代表性洪水，利用提取的所有代表性洪水生成训练集，并建立流域索引；
6)基于代表性洪水进行流域洪水响应相似性分析：利用步骤5)代表性洪水生成的训练集计算新的距离矩阵MatrixB，矩阵大小为(n×n)，n为流域个数，矩阵元素(i，j)为流域i与流域j的洪水响应相似性度量，矩阵元素(i，j)值d(i，j)的计算方法如下：
Figure BDA0002306538950000031
式中：N为步骤5)中设置的层数；di′为流域i与流域j在第i′层的代表性洪水的各种对应方式下DTW距离求和结果中最小的距离和值；
式(1)di＇具体的计算方法为：
di′＝min{sum(DDTW(Cim，Cjn))；m，n＝1，2，...，i′}，其中流域i与流域j在第i′层的代表性洪水各有i′个，共有Am n种对应方式，
Figure BDA0002306538950000041
其中n＝m＝i′，计算各种对应方式代表性洪水的DTW距离求和，di′为其中最小的DTW距离和值；Cim与Cjn为流域i与流域j在第i′层的代表性洪水，m、n为第i′层中各代表性洪水的索引；N为步骤5)中设置的层数。
完成MatrixB的计算后，依次搜索矩阵中非主对角线元素的最小值，最小值所在行列号对应的两个流域即为洪水响应最为相似的流域。
进一步的，步骤2)中采用算法从连续时间序列中自动提取洪水事件，具体方法为：
2-1基流分割，将径流序列中的基流部分和洪水部分划分开，通过数字滤波法得到流量过程的基流部分，如式(2)所示：
Figure BDA0002306538950000042
式中，bt为时刻t的基流，Qt为时刻t的径流，β为滤波系数；
通过滤波次数和滤波系数来控制滤波效果；获得基流序列后，通过式(3)获得洪水序列{q1，q2，q3，...，qt}，
qt＝Qt-bt (3)
2-2对洪水序列进行平滑处理，消除噪声项和异常点的影响，平滑滤波公式如式(4)：
Figure BDA0002306538950000051
2-3识别序列中的转折点：计算洪水序列的一阶差分序列，根据一阶差分序列的正负变换判断序列转折点的位置，并对极大值与极小值加以区分；对序列首尾的处理：若首尾值为极大值，则将其去掉；设定阈值Thmin，若首尾值小于序列均值除以阈值Thmin，则设定其为极小值；记录转折点序列及各值对应的极大值、极小值；
2-4识别洪水事件的开始、结束点：设定阈值Thslp，选定第一个极小值作为开始点并向后搜索另一极小值，依据转折点数组的一阶差分序列{d1，d2，d3，...，di，...，dt}进行判断，寻找满足式(5)要求的极小值点Mini作为结束点：
Mini-Min1＜Thslp·max(|d1|，|d2|，|d3|，...，|di|) (5)
2-5对步骤2-4中提取的洪水事件进行筛选与处理，具体步骤为：a.设置阈值Thpeak，对于一次独立的洪水事件，若峰值与序列起始点或结束点的差值小于阈值Thpeak，则认为本次洪水过程量级不足以纳入考虑范围；b.设定动态坡度阈值Thdy，对于一次独立的洪水事件，动态坡度为动态坡度阈值Thdy与极差的乘积，依据动态坡度删除洪水事件前部与后部的平坦部分；c.设定时间阈值ThΔT，对于一次独立的洪水事件，如持续时间小于阈值ThΔT，则认为本次洪水过程不足以纳入考虑范围。
进一步的，步骤3)中归一化的方法采用缩放法，如下式进行：
Figure BDA0002306538950000052
其中，xmax为时间序列中的最大值，xmin为时间序列中的最小值，xi与x′i分别为缩放前后的序列数值。
进一步的，步骤4)中聚类树生成过程中时间序列聚类分析时具体步骤为：
4-1.生成初始簇：将子集中的每一个元素作为一个初始簇；
4-2.计算基于一个子集的距离矩阵：矩阵大小为(m×m)，m为该子集中包含的洪水事件个数，矩阵的元素(i，j)为i簇与j簇的相似度，表示洪水事件i与洪水事件j的相似度，使用DTW距离作为相似性度量标准，距离越小则相似性越强；
4-3.基于步骤4-2中的距离矩阵对簇进行合并，找出距离最近的两个簇且进行合并，将聚类簇重新编号，并计算新簇与其他各簇的距离，更新距离矩阵；
4-4.重复步骤4-3直至所有的聚类簇合并为一个簇，由此生成一棵聚类树；
4-5.重复步骤4-2～4-4使基于样本集合中每个子集均生成一棵对应的聚类树。
进一步的，DTW距离的计算方法为：
对时间序列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)的累积距离计算公式为：
γ(i，j)＝d(xi，yj)+min{γ(i-1，j-1)，γ(i-1，j)，γ(i，j-1)} (7)
其中，i＝1，...，m；j＝1，...，n；为矩阵元素索引；给定初始条件γ(1，1)＝d(x1，y1)，以迭代计算得到累积距离矩阵，
Figure BDA0002306538950000071
即为时间序列X与Y的DTW距离。
进一步的，步骤5)中聚类中心的计算方法为：生成距离矩阵MatrixD，矩阵大小为(mi×mi)，mi为节点中包含的洪水事件个数，矩阵的元素(i，j)为洪水事件i与洪水事件j的DTW距离，首先计算距离矩阵MatrixD，然后计算其各行的和值，和值最小的行索引所对应的洪水事件即为聚类中心。
本发明的有益效果：
本发明提出一种基于时间序列聚类的流域洪水响应相似性分析方法，区别于以往利用单一指数或者模型参数对水文相似进行判断，本发明中提出的方法能够充分利用流量时间序列数据所包含的信息，针对流域洪水响应的相似性进行判断，基于该方法分析的相似性能够有效提高洪水预报的准确性。
下面结合附图及具体实施方式对本发明作进一步详细说明。
附图说明
图1为本发明方法整体流程图；
图2为洪水数据插值示意图；
图3为洪水场次提取示意图；
图4为时间序列的动态扭曲路径；
图5为聚类树示意图；
图6为聚类中心示意图；
图7为子集1的洪水样本情况；
图8为子集2的洪水样本情况；
图9为子集3的洪水样本情况；
图10为子集1聚类树；
图11为子集1代表性洪水；
图12为流域洪水响应相似性矩阵。
具体实施方式
实施例1
1)数据的收集、处理与保存
收集各待分析流域出口水文站点的洪水数据，根据机器学习对数据量的需求，洪水数据需要覆盖10年或10年以上。
将原始数据处理为等时段时间序列，若原始数据为非等时段数据，则需对数据进行插值处理，对于洪水数据建议采取线性内插，如图2所示，利用原始序列{Q1，Q2，Q3，...，Q7}插值获得等时段流量时间序列{Q′1，Q′2，Q′3，...，Q′12}。
将处理好的等时段流域出口流量时间序列数据保存至数据库，并使用流域编码进行标识。
2)场次洪水自动提取
根据流域编码依次读取数据库中各流域的出口流量时间序列，以洪水的起涨和消落过程作为洪水事件的划分标准将其划分为独立的场次洪水。由于机器学习对数据量的要求较大，依靠人工划分效率较低，因而采用算法从连续时间序列中自动提取洪水事件。具体方法为：
2-1.基流分割，将径流序列中的基流部分和洪水部分划分开，通过数字滤波法得到流量过程的基流部分，如下式所示。
Figure BDA0002306538950000091
其中bt为时刻t的基流，Qt为时刻t的径流，β为滤波系数，通过滤波次数和滤波系数来控制滤波效果。获得基流序列后，通过下式获得洪水序列{q1，q2，q3，...，qt}。
qt＝Qt-bt (3)
2-2.对洪水序列进行平滑处理，消除噪声项和异常点的影响，平滑滤波公式如下：
Figure BDA0002306538950000092
2-3.识别序列中的转折点，具体步骤为：计算洪水序列的一阶差分序列，根据一阶差分序列的正负变换判断序列转折点的位置，并对极大值与极小值加以区分。对序列首尾的处理：若首尾值为极大值，则将其去掉；设定阈值Thmin，若首尾值小于序列均值除以阈值Thmin，则设定其为极小值。记录转折点序列及各值对应的峰(极大值)、谷(极小值)标记。
2-4.识别洪水事件的开始、结束点，具体步骤为：设定阈值Thslp，选定第一个极小值作为开始点并向后搜索另一极小值，依据转折点数组的一阶差分序列进行判断{d1，d2，d3，...，di，...，dt}，寻找满足下式要求的极小值点Mini作为结束点：
Mini-Min1＜Thslp·max(|d1|，|d2|，|d3|，...，|di|) (5)
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，则认为是一次独立的洪水事件。
依据上述方法，得到n′个场次洪水序列{Qi1，Qi2，...，Qik，}及其时间标识序列{Ti1，Ti2，...，Tik，}，其中i＝1，...，n′，n′为洪水场次个数，k′为该场洪水对应的时段个数。
3)生成洪水事件样本集合
基于步骤2)中提取的场次洪水时间序列生成洪水事件样本集合{S1，S2，...，Si，...，Sn}，其中包含n个子集，分别记录n个不同流域的场次洪水数据。设子集Si为流域i的场次洪水样本集合，其元素为步骤2)中基于流域i的流量数据提取的且经过归一化处理n′场洪水，归一化的方法可以采用缩放法，如下式所示：
Figure BDA0002306538950000111
其中，xmax为时间序列中的最大值，xmin为时间序列中的最小值，xi与x′i分别为缩放前后的序列数值。
将归一化的场次洪水以时间序列的形式进行保存，同时建立子集及洪水事件索引。
4)基于洪水事件样本生成聚类树
以步骤3)集合中的子集为单元进行洪水事件层次聚类，生成n棵聚类树。基于索引遍历各子集，针对单个子集的洪水数据，时间序列聚类分析的具体步骤为：
4-1.生成初始簇，将子集中的每一个元素作为一个初始簇，对于一个具有m个元素的集合D＝{x1，x2，...，xm}，设定初始簇集合C＝{C1，C2，...，Cm}，其中Cj＝{xj}；
4-2.计算第一距离矩阵MatrixF，矩阵大小为(m×m)，m为子集中包含的洪水事件个数，矩阵的元素(i，j)为洪水事件i与洪水事件j的相似度，因而主对角线元素为0且为对称矩阵。使用DTW距离作为相似性度量标准，距离越小则相似性越强，DTW距离计算方法如下：
对时间序列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。
能够满足约束条件的路径有很多条，此处寻找扭曲代价最小的路径，即：
Figure BDA0002306538950000131
其中d(wk)为wk代表的两个对应元素间的距离。
根据动态规划思想，若点(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)的累积距离计算公式：
γ(i，j)＝d(xi，yj)+min{γ(i-1，j-1)，γ(i-1，j)，γ(i，j-1)} (7)
给定初始条件γ(1，1)＝d(x1，y1)，可以迭代计算得到累积距离矩阵。
Figure BDA0002306538950000132
即为时间序列X与Y的DTW距离，从点γ(m，n)出发反向搜索累积距离矩阵即可得到最佳匹配路径。
4-3.对簇进行合并，找出距离最近的两个簇Ci*和Cj*，合并Ci*和Cj*：Ci*＝Ci*∪Cj*，将聚类簇重新编号，删除当前距离矩阵的第j*行和第j*列，并计算新簇与其他各簇的距离，更新距离矩阵；
4-4.重复上一步骤直至所有聚类簇合并为一个簇，由此生成一棵聚类树，如图5所示。
4-5.重复步骤4-2～4-4使基于样本集合中每个子集均生成一棵对应的聚类树。
5)基于聚类树提取各流域代表性洪水
依次对n个流域的聚类树进行分析，提取聚类中心。对于单个聚类树，树的各节点均代表一类洪水事件，而每个节点的聚类中心即为此节点最具代表性的洪水事件，如图6所示。将根节点作为聚类树的第1层，则聚类树的第n层包含n个节点及n个聚类中心。将节点样本中与其他各元素的距离和值最小的元素作为节点的聚类中心，则对于节点i，其聚类中心的计算方法为：
生成距离矩阵MatrixD，矩阵大小为(mi×mi)，mi为节点中包含的洪水事件个数，矩阵的元素(i，j)为洪水事件i与洪水事件j的DTW距离。首先计算距离矩阵MatrixD，然后计算其各行的和值，和值最小的行索引所对应的洪水事件即为聚类中心。
根据样本总数量的大小设定参数N，计算第1层至第N层各节点的聚类中心，作为该流域的代表性洪水。利用提取的所有代表性洪水生成训练集，并建立流域索引。
6)基于代表性洪水进行流域洪水响应相似性分析
利用步骤5)中生成的代表性洪水训练集计算新的距离矩阵MatrixB，矩阵大小为(n×n)，n为流域个数，矩阵元素(i，j)为流域i与流域j的洪水响应相似性度量，因而矩阵为对称矩阵且主对角线元素为0。矩阵元素(i，j)值d(i，j)的计算方法如下：
Figure BDA0002306538950000151
式中：di′＝min{sum(DDTW(Cim，Cjn))；m，n＝1，2，...，i′}，其中流域i与流域j在第i′层的代表性洪水各有i′个，共有Am n种对应方式，
Figure BDA0002306538950000152
其中n＝m＝i′，计算各种对应方式代表性洪水的DTW距离求和，di′为其中最小的DTW距离和值；Cim与Cjn为流域i与流域j在第i′层的代表性洪水，m、n为第i′层中各代表性洪水的索引；
Figure BDA0002306538950000153
N为步骤5)中设置的层数；
完成MatrixB的计算后，依次搜索矩阵中非主对角线元素的最小值，最小值所在行列号对应的两个流域即为洪水响应最为相似的流域。
收集到我国黄河流域中游49个子流域出口水文站点的历史洪水数据，数据起止时间如下表所示，数据年限均在10年以上。
表1流域资料情况表
Figure BDA0002306538950000161
Figure BDA0002306538950000171
经过插值处理为等时段数据后，进行场次洪水的自动提取，各流域提取的洪水场次情况如下表所示：
表2流域场次洪水提取情况
Figure BDA0002306538950000172
Figure BDA0002306538950000181
对各场次洪水进行标准化处理，以样本子集1、2、3为例，三个子集的洪水样本情况如图7～9所示：
以各洪水样本子集为单元进行洪水事件层次聚类，共生成49棵聚类树，以第1子集的聚类树为例，如图10所示。
根据生成的聚类树，设定参数N＝2，提取各簇的聚类中心，即流域的代表性洪水，以第1个样本子集的代表性洪水为例，如图11所示：
根据代表性洪水生成大小为49×49的距离矩阵MatrixB，即为流域洪水响应相似性矩阵，如图12所示。
查找矩阵元素中的最小值所在行列号，得到流域3与流域5的洪水响应最为相似。

Claims (6)

1.一种基于时间序列聚类的流域洪水响应相似性分析方法，其特征在于：包括以下步骤：
1)数据的收集、处理与保存：收集待分析流域出口水文站点的径流数据，获得流域出口流量时间序列，进行等时段处理获得等时段流域出口流量时间序列数据保存至数据库中，并使用流域编码进行标识；
2)洪水场次自动划分：根据流域编码依次读取数据库中各流域出口流量时间序列，以洪水的起涨和消落过程作为洪水事件的划分标准将其划分为独立的场次洪水；
3)生成洪水事件样本集合：基于步骤2)中划分的场次洪水生成洪水事件样本集合，其中包含n个子集，分别记录n个不同流域的场次洪水数据，子集Si为流域i的场次洪水样本集合，其元素为基于流域i的流量数据提取的且经过归一化处理的场次洪水，将归一化的场次洪水以时间序列的形式进行保存，同时建立子集及洪水事件索引；
4)基于洪水事件样本生成聚类树：以步骤3)集合中的子集为单元进行洪水事件层次聚类，生成n棵聚类树，每棵聚类树对应一个流域；
5)基于聚类树提取各流域代表性洪水：依次对n个流域的聚类树进行分析，根据样本总数量大小设置参数N，搜索每棵聚类树第1层至第N层的各个节点，计算各节点的聚类中心，作为该流域的代表性洪水，利用提取的所有代表性洪水生成训练集，并建立流域索引；
6)基于代表性洪水进行流域洪水响应相似性分析：利用步骤5)代表性洪水生成的训练集计算新的距离矩阵MatrixB，矩阵大小为(n×n)，n为流域个数，矩阵元素(i，j)为流域i与流域j的洪水响应相似性度量，矩阵元素(i，j)值d(i，j)的计算方法如下：
Figure FDA0002557570170000021
式中：di′＝min{sum(DDTW(Cim，Cjn))；m，n＝1，2，...，i′}，其中流域i与流域j在第i′层的代表性洪水各有i′个，共有
Figure FDA0002557570170000022
种对应方式，计算各种对应方式代表性洪水的DTW距离求和，di′为其中最小的DTW距离和值；Cim与Cjn为流域i与流域j在第i′层的代表性洪水，m、n为第i′层中各代表性洪水的索引；N为步骤5)中设置的层数；
完成MatrixB的计算后，依次搜索矩阵中非主对角线元素的最小值，最小值所在行列号对应的两个流域即为洪水响应最为相似的流域。
2.根据权利要求1所述的基于时间序列聚类的流域洪水响应相似性分析方法，其特征在于：
步骤2)中采用算法从连续时间序列中自动提取洪水事件，具体方法为：
2-1基流分割，将径流序列中的基流部分和洪水部分划分开，通过数字滤波法得到流量过程的基流部分，如式(2)所示：
Figure FDA0002557570170000023
式中，bt为时刻t的基流，Qt为时刻t的径流，β为滤波系数；
通过滤波次数和滤波系数来控制滤波效果；获得基流序列后，通过式(3)获得洪水序列{q1，q2，q3，...，qt}，
qt＝Qt-bt (3)
2-2对洪水序列进行平滑处理，消除噪声项和异常点的影响，平滑滤波公式如式(4)：
Figure FDA0002557570170000031
2-3识别序列中的转折点：计算洪水序列的一阶差分序列，根据一阶差分序列的正负变换判断序列转折点的位置，并对极大值与极小值加以区分；对序列首尾的处理：若首尾值为极大值，则将其去掉；设定阈值Thmin，若首尾值小于序列均值除以阈值Thmin，则设定其为极小值；记录转折点序列及各值对应的极大值、极小值；
2-4识别洪水事件的开始、结束点：设定阈值Thslp，选定第一个极小值作为开始点并向后搜索另一极小值，依据转折点数组的一阶差分序列{d1，d2，d3，...，di，...，dt}进行判断，寻找满足式(5)要求的极小值点Mini作为结束点：
Mini-Min1＜Thslp·max(|d1|，|d2|，|d3|，...，|di|) (5)
2-5对步骤2-4中提取的洪水事件进行筛选与处理，具体步骤为：a.设置阈值Thpeak，对于一次独立的洪水事件，若峰值与序列起始点或结束点的差值小于阈值Thpeak，则认为本次洪水过程量级不足以纳入考虑范围；b.设定动态坡度阈值Thdy，对于一次独立的洪水事件，动态坡度为动态坡度阈值Thdy与极差的乘积，依据动态坡度删除洪水事件前部与后部的平坦部分；c.设定时间阈值ThΔT，对于一次独立的洪水事件，如持续时间小于阈值ThΔT，则认为本次洪水过程不足以纳入考虑范围。
3.根据权利要求1所述的基于时间序列聚类的流域洪水响应相似性分析方法，其特征在于：步骤3)中归一化的方法采用缩放法，如下式进行：
Figure FDA0002557570170000041
其中，xmax为时间序列中的最大值，xmin为时间序列中的最小值，xi与x′i分别为缩放前后的序列数值。
4.根据权利要求1所述的基于时间序列聚类的流域洪水响应相似性分析方法，其特征在于：步骤4)中聚类树生成过程中时间序列聚类分析时具体步骤为：
4-1.生成初始簇：将子集中的每一个元素作为一个初始簇；
4-2.计算基于一个子集的距离矩阵：矩阵大小为(m×m)，m为该子集中包含的洪水事件个数，矩阵的元素(i,j)为i簇与j簇的相似度，表示洪水事件i与洪水事件j的相似度，使用DTW距离作为相似性度量标准，距离越小则相似性越强；
4-3.基于步骤4-2中的距离矩阵对簇进行合并，找出距离最近的两个簇且进行合并，将聚类簇重新编号，并计算新簇与其他各簇的距离，更新距离矩阵；
4-4.重复步骤4-3直至所有的聚类簇合并为一个簇，由此生成一棵聚类树；
4-5.重复步骤4-2～4-4使基于样本集合中每个子集均生成一棵对应的聚类树。
5.根据权利要求4所述的基于时间序列聚类的流域洪水响应相似性分析方法，其特征在于：DTW距离的计算方法为：
对时间序列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)的累积距离计算公式为：
γ(i，j)＝d(xi，yj)+min{γ(i-1，j-1)，γ(i-1，j)，γ(i，j-1)} (7)
其中，i＝1,…,m；j＝1,…,n；为矩阵元素索引；给定初始条件γ(1，1)＝d(x1，y1)，迭代计算得到累积距离矩阵，
Figure FDA0002557570170000051
即为时间序列X与Y的DTW距离。
6.根据权利要求1所述的基于时间序列聚类的流域洪水响应相似性分析方法，其特征在于：步骤5)中聚类中心的计算方法为：生成距离矩阵MatrixD，矩阵大小为(mi×mi)，mi为节点中包含的洪水事件个数，矩阵的元素(i，j)为洪水事件i与洪水事件j的DTW距离，首先计算距离矩阵MatrixD，然后计算其各行的和值，和值最小的行索引所对应的洪水事件即为聚类中心。
CN201911242133.9A 2019-12-06 2019-12-06 一种基于时间序列聚类的流域洪水响应相似性分析方法 Active CN111080107B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911242133.9A CN111080107B (zh) 2019-12-06 2019-12-06 一种基于时间序列聚类的流域洪水响应相似性分析方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911242133.9A CN111080107B (zh) 2019-12-06 2019-12-06 一种基于时间序列聚类的流域洪水响应相似性分析方法

Publications (2)

Publication Number Publication Date
CN111080107A CN111080107A (zh) 2020-04-28
CN111080107B true CN111080107B (zh) 2020-09-15

Family

ID=70312978

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911242133.9A Active CN111080107B (zh) 2019-12-06 2019-12-06 一种基于时间序列聚类的流域洪水响应相似性分析方法

Country Status (1)

Country Link
CN (1) CN111080107B (zh)

Families Citing this family (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113344711A (zh) * 2021-04-12 2021-09-03 中国再保险(集团)股份有限公司 一种基于拉丁超立方抽样方法进行巨灾事件模拟的方法及计算机程序产品
CN113159408B (zh) * 2021-04-14 2023-11-21 交控科技股份有限公司 轨道交通站点客流预测方法及装置
CN113487070B (zh) * 2021-06-23 2023-10-10 中国长江三峡集团有限公司 洪水频率的分析方法、装置和计算机设备
CN114154417B (zh) * 2021-12-06 2022-06-07 中国水利水电科学研究院 基于深度学习框架的洪水预报模型及洪水预报方法
CN114240106B (zh) * 2021-12-06 2022-07-01 中国水利水电科学研究院 一种基于水文数据挖掘的流域洪水响应相似性分析方法
CN114580171B (zh) * 2022-03-03 2022-09-30 中国科学院地理科学与资源研究所 一种流域洪水类型辨识及其影响因子解析的方法
CN115063111B (zh) * 2022-06-24 2023-08-18 中国长江三峡集团有限公司 场次洪水识别方法、装置、电子设备及可读存储介质
CN117951540B (zh) * 2023-10-19 2024-10-18 中国水利水电科学研究院 一种基于特征立方体的流域相似性分析方法
CN117574778B (zh) * 2024-01-12 2024-03-29 河海大学 一种基于机器学习的相似场次洪水模式库构建方法

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108537247A (zh) * 2018-03-13 2018-09-14 河海大学 一种时空多元水文时间序列相似性度量方法

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10453026B2 (en) * 2015-03-06 2019-10-22 Walmart Apollo, Llc System and method for forecasting high-sellers using multivariate bayesian time series
CN104933268B (zh) * 2015-07-13 2018-09-04 国家电网公司 一种基于一维非恒定流数值模型的洪水分析方法
CN105069093B (zh) * 2015-08-05 2018-07-24 河海大学 一种基于嵌入式索引的水文时间序列相似性搜索方法
CN108846573B (zh) * 2018-06-12 2021-04-09 河海大学 基于时间序列核距离的流域水文相似性估计方法
CN109190261A (zh) * 2018-09-07 2019-01-11 中国水利水电科学研究院 一种一维河网概化与高性能一二维耦合的洪水分析方法
CN109325302A (zh) * 2018-10-08 2019-02-12 云南省水利水电勘测设计研究院 特小流域暴雨洪水分析模型构建方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108537247A (zh) * 2018-03-13 2018-09-14 河海大学 一种时空多元水文时间序列相似性度量方法

Also Published As

Publication number Publication date
CN111080107A (zh) 2020-04-28

Similar Documents

Publication Publication Date Title
CN111080107B (zh) 一种基于时间序列聚类的流域洪水响应相似性分析方法
CN111027763B (zh) 一种基于机器学习的流域洪水响应相似性分析方法
CN106650767B (zh) 基于聚类分析和实时校正的洪水预报方法
CN108304668B (zh) 一种结合水文过程数据和历史先验数据的洪水预测方法
CN111027764B (zh) 一种基于机器学习的适用于径流资料缺乏流域的洪水预报方法
CN104821082B (zh) 一种基于综合评价的短时交通流预测方法
CN117035201B (zh) 平原河网水工程集群多目标调度规则制定方法及系统
CN113705877B (zh) 基于深度学习模型的实时月径流预报方法
CN104090974B (zh) 展延水库后续来水的动态数据挖掘方法及系统
CN115495991A (zh) 一种基于时间卷积网络的降水区间预测方法
CN114897242B (zh) 一种自适应分割时段场次降雨径流的方法
CN111898831A (zh) 一种实时洪水概率预报实用化方法
CN110930282B (zh) 一种基于机器学习的局地降雨雨型分析方法
CN114154417A (zh) 基于深度学习框架的洪水预报模型及洪水预报方法
CN112561214A (zh) 一种自动识别场次洪水的方法及系统
CN114219252B (zh) 一种基于sce-ua算法的流域单位线分析方法
CN113779113B (zh) 基于雨洪时空过程相似性挖掘的洪水动态预估方法及系统
CN112036687A (zh) 一种梯级水库群防洪联合调度规则决策树获取方法
CN116911643B (zh) 城市内涝防治方案选择方法、系统、装置、存储介质
CN118013277A (zh) 具有时变权重的多模型组合径流预报方法
CN116167513A (zh) 基于单变量优化dmca模型的流域洪水响应时间计算方法
CN116502531A (zh) 一种基于多元线性回归模型的基流模拟方法
CN115688022B (zh) 一种基于最邻近算法的流域单位线实时优选方法
CN109284286B (zh) 一种从原始数据集中提取有效特征的方法
CN112215495A (zh) 一种基于长短时记忆神经网络的污染源贡献度计算方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant