(1.云南省水文水資源局 昭通分局，云南 昭通 657000; 2.長江科學(xué)院 農(nóng)業(yè)水利研究所，湖北 武漢 430010; 3.西藏自治區(qū)水土保持局，西藏 拉薩 850000)

基于DTW算法的云南昭通地區(qū)場次洪水相似性研究

任繼周1彭德才1喬偉2萍央3

(1.云南省水文水資源局昭通分局，云南昭通657000; 2.長江科學(xué)院農(nóng)業(yè)水利研究所，湖北武漢430010; 3.西藏自治區(qū)水土保持局，西藏拉薩850000)

利用隨機水文學(xué)思想分析場次洪水的歷史相似性對洪水預(yù)報、防洪調(diào)度等具有重要的現(xiàn)實意義。采用動態(tài)時間扭曲算法，對云南省昭通市綏江縣大汶溪流域1961～2014年洪峰流量進行了頻率分析，并對新華水文站場次洪水過程進行了相似性分析。分析得出，動態(tài)時間扭曲算法作為數(shù)據(jù)挖掘領(lǐng)域一種相似性分析的方法，能夠?qū)υ紨?shù)據(jù)進行充分挖掘，從而有效尋找出不同場次洪水中相似性最高的組合。研究成果可為場次洪水特征分析、洪水預(yù)報等提供參考。

場次洪水；相似性；動態(tài)時間扭曲；云南昭通

1 研究背景

洪水是我國江河的主要自然災(zāi)害，直接威脅著區(qū)域社會經(jīng)濟的發(fā)展和人民生命財產(chǎn)安全。目前，我國學(xué)者對洪水的分析主要集中在三方面：①洪水的頻率或重現(xiàn)期分析，如洪峰流量的頻率曲線擬合及洪峰、洪量、洪水歷時等多要素聯(lián)合概率分析等；②場次洪水過程的特征分析，如洪水過程形態(tài)、漲洪歷時等；③洪水過程與暴雨過程的關(guān)系，包括產(chǎn)流系數(shù)、徑流組成、洪水峰現(xiàn)時間滯后最大雨峰的時間等。從水文時間序列分析的角度來看，以往對場次洪水歷史相似性的研究較少。根據(jù)隨機水文學(xué)的思想，未來的洪水過程是歷史上某種洪水過程的重現(xiàn)，充分利用洪水信息，可對未來可能發(fā)生的洪水形勢進行預(yù)估，從而掌握洪水全貌[1]。因此，通過隨機水文學(xué)思想分析場次洪水的歷史相似性對洪水預(yù)報、防洪調(diào)度等具有重要的現(xiàn)實意義。

時間序列相似性研究中的動態(tài)時間扭曲算法(Dynamic Time Warping，簡稱DTW)最早起源于語音識別領(lǐng)域，隨后被廣泛應(yīng)用于各種領(lǐng)域進行數(shù)據(jù)挖掘，整個計算過程等同于動態(tài)尋優(yōu)[2]。該方法已被應(yīng)用于水文時間序列的相似性分析，并取得了較好的效果。因此，本文將該方法應(yīng)用于云南省昭通市綏江縣大汶溪流域，在對流域內(nèi)新華水文站1961～2014年洪峰流量進行頻率分析的基礎(chǔ)上，采用DTW方法對實測54場典型洪水過程進行場次洪水相似性分析，探尋歷史洪水的相似性，以期為場次洪水特征分析、洪水預(yù)報等提供參考。

2 資料與方法

2.1 數(shù)據(jù)資料

綏江縣位于昭通市北部，東與水富縣相鄰，南與鹽津縣、永善縣接壤，西、北面與四川省的屏山縣、雷波縣隔江相望；地理坐標介于東經(jīng)103°47′～104°16′、北緯28°21′～28°40′之間。縣境東西長約48.5 km、南北寬36 km，總面積 746.3 km2。綏江縣屬亞熱帶濕潤季風氣候，冬無嚴寒，夏無酷暑。特殊的氣候條件、地形地貌和人類活動的影響使綏江縣成為山洪地質(zhì)災(zāi)害的重災(zāi)區(qū)。大汶溪流域?qū)儆诮鹕辰野兜囊患壷Я鳎l(fā)源于綏江縣板栗鄉(xiāng)大堡頂，流經(jīng)板栗鄉(xiāng)后稱板栗河，主要支流有銅廠河、黃龍溪，全流域面積391 km2。新華水文站于1961年設(shè)立，位于流域下游，屬于出口控制站，控制流域面積324 km2，河長35.7 km，河道坡度29.0‰，流域平均高程1 214 m。

2.2 研究方法

DTW算法具體過程如下[3]：假定兩個時間序列為A=和B=，構(gòu)建一個k×s階的矩陣，該矩陣中存放的元素即為兩個時間序列對應(yīng)點間的距離，即矩陣中的第(i，j)個元素為點ai和點bj之間的距離d(xi，yj)，距離的計算通常采用d(xi,yi)=(xi-yi)2。上述k×s階矩陣中，元素的連續(xù)集構(gòu)成了扭曲路徑W，其中W的第m個元素被定義為wm=(i，j)m，因此可以得到一個路徑集：

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

M

(1)

扭曲路徑要求滿足連續(xù)性、邊界條件以及單調(diào)性條件限制。

(1)邊界條件：w1=(1，1)，wM= (k，s)。

(2)連續(xù)性：給定wm= (a，b)，wm-1=(a′，b′)，要求a-a′≤1和b-b′≤1，該連續(xù)性要求扭曲路徑必須是連續(xù)的，不能有跳過某點的情況。

(3)單調(diào)性：給定wm= (a，b)，wm-1=(a′，b′)，要求a-a′≥0和b-b′≥0。也就是說，在動態(tài)尋優(yōu)的過程中，不能再向相反的方向?qū)?yōu)。

從以上要求可知，對于兩個時間序列，有很多路徑都滿足以上條件，但最優(yōu)的路徑只有一條。為了尋找最佳路徑，需要定義從起點到終點的累積距離，并根據(jù)累積距離求得扭曲代價，累積距離最小則為最優(yōu)，最優(yōu)路徑滿足最小扭曲代價，如下式：

(2)

構(gòu)造累積距離矩陣r公式如下

(3)

初始條件設(shè)置為r(1,1)=d(p1,q1)。從兩個序列起始點(1，1)開始，根據(jù)公式(2)與公式(3)迭代計算累積距離，并最終得到最小累加值r(k,s)，該累加值即為時間序列P和Q的dDTW距離。

3 結(jié)果與討論

3.1 洪峰流量頻率分析

首先對新華水文站1961～2014年實測洪峰流量進行頻率分析。頻率分析采用國內(nèi)常用的P-Ⅲ型曲線進行。經(jīng)過適線，頻率分析結(jié)果如圖1所示。其中，新華水文站洪峰流量均值約為460 m3/s，變差系數(shù)Cv值為 0.51，Cs/Cv=3.0。從圖1可以看出，P-Ⅲ型曲線的擬合效果較好，且該水文站洪峰流量的頻率幾乎均大于5%，即新華水文站洪峰流量重現(xiàn)期基本小于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對應(yīng)的重現(xiàn)期約為2.5 a。

圖1 新華水文站洪峰流量P-Ⅲ型曲線

3.2 場次洪水相似性分析

在對場次洪水進行相似性分析的過程中，重點考慮場次洪水的形狀特征，即從洪水起漲到出現(xiàn)洪峰再到洪水退去的整個過程中。洪水流量的變化過程對研究當?shù)睾樗兓螒B(tài)、制定山洪預(yù)警方案有重要意義。在1961～2014這54 a中，每年的場次洪水不僅僅只有一場，有時連續(xù)兩場降雨間隔較短，導(dǎo)致由降雨引起的兩場洪水出現(xiàn)疊加的情況，不利于單次洪水過程的分析。因此，本文在場次洪水相似性分析的過程中，每年僅選取一場洪水進行研究。選取的原則為：每年選取一場由單次降雨引起的單場次洪水，若符合該類型的洪水超過一場，則在其中選取洪峰流量最大的那一場洪水作為研究對象。為描述洪水的變化過程，從洪水起漲算起，取洪水起漲后3，6，9，12，15 h及18 h對應(yīng)的洪峰流量共7個數(shù)字描述場次洪水過程，如表1中所列即為1961～1967年場次洪水待分析數(shù)據(jù)。

表1 代表年場次洪水流量過程 m3/s

準備好待分析場次洪水數(shù)據(jù)后，便可采用DTW算法對歷年場次洪水過程進行相似性分析。DTW算法的目的是得到不同年份場次洪水過程數(shù)據(jù)間的dDTW距離，并用一個方陣來存放所有子序列之間的dDTW距離，記為D_matrix：

D_matrix=

(4)

式中，f1，f2，…，fn為n個子序列，即子序列樣本數(shù)為n。矩陣的第(k，s)元素代表第k個子序列和第s個子序列之間的dDTW距離，即dDTW(fk，fs)。根據(jù)矩陣的定義可知，D_matrix是一個對稱矩陣，因為序列自身與自身完全相似，故矩陣對角線為0。根據(jù)DTW算法原理可知，在距離方陣中，數(shù)值越小則兩組數(shù)據(jù)越相似。計算得到的D_matrix矩陣如下所示。

(5)

從D_matrix的矩陣計算結(jié)果可以看出，不同年份之間的場次洪水相似性存在較大的差異。其中，計算得到最小的dDTW距離為75.8，對應(yīng)年份為1976年與1993年，說明這兩場洪水的相似度最高；計算得到最大的dDTW距離為2 184，對應(yīng)年份為1974年與2011年，說明這兩場洪水的相似度最低。

根據(jù)矩陣的計算結(jié)果，分別在矩陣中選擇4組dDTW距離最小的年份繪圖對比，所選擇的年份及其dDTW距離列表2如下，場次洪水過程對比如圖2所示。

表2 相似度最高的4組場次洪水對比

從圖2可以看出，1976年與1993年的場次洪水過程相似性最高，這兩場洪水過程線從起漲到洪峰，再到退水的整個過程非常相似，幾乎完全重合。其他3組場次洪水的過程線相似度也較高，不論是洪水過程線的形態(tài)還是洪峰流量，差別并不大。這說明動態(tài)時間扭曲算法能夠在眾多場次洪水數(shù)據(jù)中尋找出相似性最高的組合，其相似性搜索的效果顯著。同時，從圖2中的4幅圖中可以看出，新華水文站實測得到的場次洪水形態(tài)大致相似，屬于尖瘦峰型，即在洪水過程中，漲洪歷時較短且漲勢兇猛，退水過程相較漲水過程緩慢。此種類型洪水由于漲水快、洪量大，易造成山洪災(zāi)害，是山洪災(zāi)害中最需要預(yù)警的一種洪水類型。

圖2 相似性最高的4組場次洪水對比

根據(jù)D_matix計算結(jié)果，將矩陣各行(或各列)相加，即可得到每個年份場次洪水與其他年份場次洪水相似度的計量總和，如表3所示。從表3可以看出，ΣdDTW最小的年份是1966年，與其他年份dDTW距離之和為19 868.05；ΣdDTW最大的年份是1974年，該年份與其他年份dDTW距離之和為67 344.58。這說明，1966年的場次洪水過程與其他年份的相似度最高，而1974年的場次洪水過程與其他年份相似度最低。

表3 各年份dDTW距離之和

4 結(jié) 論

本文利用動態(tài)時間扭曲算法，在對云南省昭通市綏江縣大溪溝流域1961～2014年洪峰流量進行頻率分析的基礎(chǔ)上，對新華水文站場次洪水過程進行了相似性分析，主要結(jié)論如下：

(1)新華水文站洪峰流量均值約為460 m3/s，變差系數(shù)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。均值洪峰流量對應(yīng)的重現(xiàn)期約為 2.5 a。

(2)經(jīng)過動態(tài)時間扭曲算法的相似性搜索，1976年與1993年的場次洪水相似性最高，其dDTW距離為 75.8；1974年與2011年的場次洪水相似度最低，其dDTW距離為2 184。ΣdDTW最小的年份是1966年，說明1966年的場次洪水過程與其他年份的相似度最高；ΣdDTW最大的年份是1974年，說明1974年的場次洪水過程與其他年份相似度最低。

總之，動態(tài)時間扭曲算法作為數(shù)據(jù)挖掘領(lǐng)域一種相似性分析的方法，能夠?qū)υ紨?shù)據(jù)進行充分挖掘，從而有效尋找出不同場次洪水中相似性最高的組合。該方法可在水文領(lǐng)域繼續(xù)推廣，如應(yīng)用到暴雨過程與洪水過程及洪峰與洪量過程的相似性分析中。

[1] ?？? 流域場次暴雨洪水相似性分析的可拓模型構(gòu)建及應(yīng)用[D]. 南京：河海大學(xué), 2006.

[2] 劉衛(wèi)林, 董增川, 梁忠民,等. 暴雨洪水相似性分析及其應(yīng)用研究[J]. 中國農(nóng)村水利水電, 2007(2):132-135.

[3] 楊艷林, 葉楓, 呂鑫,等. 一種基于DTW聚類的水文時間序列相似性挖掘方法[J]. 計算機科學(xué), 2016, 43(2):245-249.

(編輯：朱曉紅)

2017-05-09

國家自然科學(xué)基金(51509010)；中央級公益性科研院所基本科研業(yè)務(wù)費資助項目(CKSF2016039/NS、CKSF2016028/NS)

任繼周，男，云南省水文水資源局昭通分局，高級工程師.

喬偉，女，長江科學(xué)院農(nóng)業(yè)水利研究所，工程師.

1006-0081(2017)08-0035-04

TV122

：A