基于DTW算法的云南昭通地区场次洪水相似性研究

(1.云南省水文水资源局 昭通分局，云南 昭通 657000; 2.长江科学院 农业水利研究所，湖北 武汉 430010; 3.西藏自治区水土保持局，西藏 拉萨 850000)

基于DTW算法的云南昭通地区场次洪水相似性研究

任继周1彭德才1乔伟2萍央3

(1.云南省水文水资源局昭通分局，云南昭通657000; 2.长江科学院农业水利研究所，湖北武汉430010; 3.西藏自治区水土保持局，西藏拉萨850000)

利用随机水文学思想分析场次洪水的历史相似性对洪水预报、防洪调度等具有重要的现实意义。采用动态时间扭曲算法，对云南省昭通市绥江县大汶溪流域1961～2014年洪峰流量进行了频率分析，并对新华水文站场次洪水过程进行了相似性分析。分析得出，动态时间扭曲算法作为数据挖掘领域一种相似性分析的方法，能够对原始数据进行充分挖掘，从而有效寻找出不同场次洪水中相似性最高的组合。研究成果可为场次洪水特征分析、洪水预报等提供参考。

场次洪水；相似性；动态时间扭曲；云南昭通

1 研究背景

洪水是我国江河的主要自然灾害，直接威胁着区域社会经济的发展和人民生命财产安全。目前，我国学者对洪水的分析主要集中在三方面：①洪水的频率或重现期分析，如洪峰流量的频率曲线拟合及洪峰、洪量、洪水历时等多要素联合概率分析等；②场次洪水过程的特征分析，如洪水过程形态、涨洪历时等；③洪水过程与暴雨过程的关系，包括产流系数、径流组成、洪水峰现时间滞后最大雨峰的时间等。从水文时间序列分析的角度来看，以往对场次洪水历史相似性的研究较少。根据随机水文学的思想，未来的洪水过程是历史上某种洪水过程的重现，充分利用洪水信息，可对未来可能发生的洪水形势进行预估，从而掌握洪水全貌[1]。因此，通过随机水文学思想分析场次洪水的历史相似性对洪水预报、防洪调度等具有重要的现实意义。

时间序列相似性研究中的动态时间扭曲算法(Dynamic Time Warping，简称DTW)最早起源于语音识别领域，随后被广泛应用于各种领域进行数据挖掘，整个计算过程等同于动态寻优[2]。该方法已被应用于水文时间序列的相似性分析，并取得了较好的效果。因此，本文将该方法应用于云南省昭通市绥江县大汶溪流域，在对流域内新华水文站1961～2014年洪峰流量进行频率分析的基础上，采用DTW方法对实测54场典型洪水过程进行场次洪水相似性分析，探寻历史洪水的相似性，以期为场次洪水特征分析、洪水预报等提供参考。

2 资料与方法

2.1 数据资料

绥江县位于昭通市北部，东与水富县相邻，南与盐津县、永善县接壤，西、北面与四川省的屏山县、雷波县隔江相望；地理坐标介于东经103°47′～104°16′、北纬28°21′～28°40′之间。县境东西长约48.5 km、南北宽36 km，总面积 746.3 km2。绥江县属亚热带湿润季风气候，冬无严寒，夏无酷暑。特殊的气候条件、地形地貌和人类活动的影响使绥江县成为山洪地质灾害的重灾区。大汶溪流域属于金沙江右岸的一级支流，发源于绥江县板栗乡大堡顶，流经板栗乡后称板栗河，主要支流有铜厂河、黄龙溪，全流域面积391 km2。新华水文站于1961年设立，位于流域下游，属于出口控制站，控制流域面积324 km2，河长35.7 km，河道坡度29.0‰，流域平均高程1 214 m。

2.2 研究方法

DTW算法具体过程如下[3]：假定两个时间序列为A=和B=，构建一个k×s阶的矩阵，该矩阵中存放的元素即为两个时间序列对应点间的距离，即矩阵中的第(i，j)个元素为点ai和点bj之间的距离d(xi，yj)，距离的计算通常采用d(xi,yi)=(xi-yi)2。上述k×s阶矩阵中，元素的连续集构成了扭曲路径W，其中W的第m个元素被定义为wm=(i，j)m，因此可以得到一个路径集：

W={w1,w2,…,wm,…,wM|max(k,s)

M

(1)

扭曲路径要求满足连续性、边界条件以及单调性条件限制。

(1)边界条件：w1=(1，1)，wM= (k，s)。

(2)连续性：给定wm= (a，b)，wm-1=(a′，b′)，要求a-a′≤1和b-b′≤1，该连续性要求扭曲路径必须是连续的，不能有跳过某点的情况。

(3)单调性：给定wm= (a，b)，wm-1=(a′，b′)，要求a-a′≥0和b-b′≥0。也就是说，在动态寻优的过程中，不能再向相反的方向寻优。

从以上要求可知，对于两个时间序列，有很多路径都满足以上条件，但最优的路径只有一条。为了寻找最佳路径，需要定义从起点到终点的累积距离，并根据累积距离求得扭曲代价，累积距离最小则为最优，最优路径满足最小扭曲代价，如下式：

(2)

构造累积距离矩阵r公式如下

(3)

初始条件设置为r(1,1)=d(p1,q1)。从两个序列起始点(1，1)开始，根据公式(2)与公式(3)迭代计算累积距离，并最终得到最小累加值r(k,s)，该累加值即为时间序列P和Q的dDTW距离。

3 结果与讨论

3.1 洪峰流量频率分析

首先对新华水文站1961～2014年实测洪峰流量进行频率分析。频率分析采用国内常用的P-Ⅲ型曲线进行。经过适线，频率分析结果如图1所示。其中，新华水文站洪峰流量均值约为460 m3/s，变差系数Cv值为 0.51，Cs/Cv=3.0。从图1可以看出，P-Ⅲ型曲线的拟合效果较好，且该水文站洪峰流量的频率几乎均大于5%，即新华水文站洪峰流量重现期基本小于20 a一遇。100 a一遇的洪峰流量为1 251.86 m3/s，50 a一遇的洪峰流量为1 111.2 m3/s，20 a一遇的洪峰流量为921.7 m3/s，10 a一遇的洪峰流量为744.4 m3/s。该站点洪峰流量均值460 m3/s对应的重现期约为2.5 a。

图1 新华水文站洪峰流量P-Ⅲ型曲线

3.2 场次洪水相似性分析

在对场次洪水进行相似性分析的过程中，重点考虑场次洪水的形状特征，即从洪水起涨到出现洪峰再到洪水退去的整个过程中。洪水流量的变化过程对研究当地洪水变化形态、制定山洪预警方案有重要意义。在1961～2014这54 a中，每年的场次洪水不仅仅只有一场，有时连续两场降雨间隔较短，导致由降雨引起的两场洪水出现叠加的情况，不利于单次洪水过程的分析。因此，本文在场次洪水相似性分析的过程中，每年仅选取一场洪水进行研究。选取的原则为：每年选取一场由单次降雨引起的单场次洪水，若符合该类型的洪水超过一场，则在其中选取洪峰流量最大的那一场洪水作为研究对象。为描述洪水的变化过程，从洪水起涨算起，取洪水起涨后3，6，9，12，15 h及18 h对应的洪峰流量共7个数字描述场次洪水过程，如表1中所列即为1961～1967年场次洪水待分析数据。

表1 代表年场次洪水流量过程 m3/s

准备好待分析场次洪水数据后，便可采用DTW算法对历年场次洪水过程进行相似性分析。DTW算法的目的是得到不同年份场次洪水过程数据间的dDTW距离，并用一个方阵来存放所有子序列之间的dDTW距离，记为D_matrix：

D_matrix=

(4)

式中，f1，f2，…，fn为n个子序列，即子序列样本数为n。矩阵的第(k，s)元素代表第k个子序列和第s个子序列之间的dDTW距离，即dDTW(fk，fs)。根据矩阵的定义可知，D_matrix是一个对称矩阵，因为序列自身与自身完全相似，故矩阵对角线为0。根据DTW算法原理可知，在距离方阵中，数值越小则两组数据越相似。计算得到的D_matrix矩阵如下所示。

(5)

从D_matrix的矩阵计算结果可以看出，不同年份之间的场次洪水相似性存在较大的差异。其中，计算得到最小的dDTW距离为75.8，对应年份为1976年与1993年，说明这两场洪水的相似度最高；计算得到最大的dDTW距离为2 184，对应年份为1974年与2011年，说明这两场洪水的相似度最低。

根据矩阵的计算结果，分别在矩阵中选择4组dDTW距离最小的年份绘图对比，所选择的年份及其dDTW距离列表2如下，场次洪水过程对比如图2所示。

表2 相似度最高的4组场次洪水对比

从图2可以看出，1976年与1993年的场次洪水过程相似性最高，这两场洪水过程线从起涨到洪峰，再到退水的整个过程非常相似，几乎完全重合。其他3组场次洪水的过程线相似度也较高，不论是洪水过程线的形态还是洪峰流量，差别并不大。这说明动态时间扭曲算法能够在众多场次洪水数据中寻找出相似性最高的组合，其相似性搜索的效果显著。同时，从图2中的4幅图中可以看出，新华水文站实测得到的场次洪水形态大致相似，属于尖瘦峰型，即在洪水过程中，涨洪历时较短且涨势凶猛，退水过程相较涨水过程缓慢。此种类型洪水由于涨水快、洪量大，易造成山洪灾害，是山洪灾害中最需要预警的一种洪水类型。

图2 相似性最高的4组场次洪水对比

根据D_matix计算结果，将矩阵各行(或各列)相加，即可得到每个年份场次洪水与其他年份场次洪水相似度的计量总和，如表3所示。从表3可以看出，ΣdDTW最小的年份是1966年，与其他年份dDTW距离之和为19 868.05；ΣdDTW最大的年份是1974年，该年份与其他年份dDTW距离之和为67 344.58。这说明，1966年的场次洪水过程与其他年份的相似度最高，而1974年的场次洪水过程与其他年份相似度最低。

表3 各年份dDTW距离之和

4 结 论

本文利用动态时间扭曲算法，在对云南省昭通市绥江县大溪沟流域1961～2014年洪峰流量进行频率分析的基础上，对新华水文站场次洪水过程进行了相似性分析，主要结论如下：

(1)新华水文站洪峰流量均值约为460 m3/s，变差系数Cv值为 0.51，Cs/Cv=3.0。该水文站100 a一遇的洪峰流量为1 251.86 m3/s，50 a一遇的洪峰流量为1 111.2 m3/s，20 a一遇的洪峰流量为 921.7 m3/s，10 a一遇的洪峰流量为 744.4 m3/s。均值洪峰流量对应的重现期约为 2.5 a。

(2)经过动态时间扭曲算法的相似性搜索，1976年与1993年的场次洪水相似性最高，其dDTW距离为 75.8；1974年与2011年的场次洪水相似度最低，其dDTW距离为2 184。ΣdDTW最小的年份是1966年，说明1966年的场次洪水过程与其他年份的相似度最高；ΣdDTW最大的年份是1974年，说明1974年的场次洪水过程与其他年份相似度最低。

总之，动态时间扭曲算法作为数据挖掘领域一种相似性分析的方法，能够对原始数据进行充分挖掘，从而有效寻找出不同场次洪水中相似性最高的组合。该方法可在水文领域继续推广，如应用到暴雨过程与洪水过程及洪峰与洪量过程的相似性分析中。

[1] 牛俊. 流域场次暴雨洪水相似性分析的可拓模型构建及应用[D]. 南京：河海大学, 2006.

[2] 刘卫林, 董增川, 梁忠民,等. 暴雨洪水相似性分析及其应用研究[J]. 中国农村水利水电, 2007(2):132-135.

[3] 杨艳林, 叶枫, 吕鑫,等. 一种基于DTW聚类的水文时间序列相似性挖掘方法[J]. 计算机科学, 2016, 43(2):245-249.

(编辑：朱晓红)

2017-05-09

国家自然科学基金(51509010)；中央级公益性科研院所基本科研业务费资助项目(CKSF2016039/NS、CKSF2016028/NS)

任继周，男，云南省水文水资源局昭通分局，高级工程师.

乔伟，女，长江科学院农业水利研究所，工程师.

1006-0081(2017)08-0035-04

TV122

：A