張玉國
( 黑龍江省水文局,哈爾濱150010)
水文系統(tǒng)是一個復雜系統(tǒng),又是一個動態(tài)的非線性復合系統(tǒng)[1]。在自然環(huán)境和人類活動的強烈干擾下,河川徑流呈現(xiàn)出高度的復雜性和非線性。對徑流序列進行復雜性分析,可以反映區(qū)域水文系統(tǒng)影響因子的數(shù)量多少,可以充分挖掘其內(nèi)在信息與聯(lián)系,是認識水文系統(tǒng)動態(tài)演變特征的重要途徑。目前,復雜性測度方法主要有分形理論[2]、符號動力學[3]、混沌理論[4]、C1和C2復雜度[5]、漲落復雜度[6]等。但或多或少存在著所需數(shù)據(jù)量大、抗噪聲能力弱等缺點。有些水文時間序列的長度可能較短,所以十分需要用較短數(shù)據(jù)就能有效測度序列復雜性特征的算法。2000年,Richman J S 等[7?提出了樣本熵,它具有計算不依賴數(shù)據(jù)長度、較高一致性、對缺失數(shù)據(jù)不敏感等優(yōu)點[8],是一種較好的時間序列復雜性測度方法,已經(jīng)廣泛地應用于各個領域。因此,本文采用樣本熵方法對阿城站多年月徑流序列的復雜性特征進行測度,從而揭示其動態(tài)變化規(guī)律,以期為流域水資源的合理配置提供科學依據(jù)。
樣本熵是一種有別于近似熵的不計入自身匹配的統(tǒng)計量,是對近似熵的改進。樣本熵表示非線性動力學系統(tǒng)產(chǎn)生新信息的概率,主要用來定量地刻畫系統(tǒng)的規(guī)則度及復雜度。樣本熵值越大,序列自我相似性越低,產(chǎn)生新信息的概率越高,序列越復雜[9]。可以用SampEn( m,r,N) 來表示樣本熵,其中,N 為長度,r 為預先選定的相似容限,m 為預先選定的模式維數(shù)。樣本熵的具體算法[10-11]如下:
設原始數(shù)據(jù)為x(1) ,x(2) ,…,x( N) ,共N 點。
1) 按序號連續(xù)順序組成一組m 維矢量,從Xm(1) 到Xm( N - m +1) ,其中:

2) 定義X( i) 與X( j) 之間的距離d[Xm( i) ,Xm( j) ]為兩者對應元素中差值最大的一個,即:

求其對所有i 的平均值,即:

4) 將維數(shù)增加1,變?yōu)閙+1 維矢量,重復式(1)~式(3) 的步驟后,得:

式中:Bm( r) 和Bm+1( r) 分別為m 點和m +1 點兩序列相似的概率,則序列理論上的樣本熵為:

當N 為有限時,得出序列樣本熵估計值,即:

參數(shù)m、r 的選擇是樣本熵估計的關鍵,但目前為止沒有最佳標準,根據(jù)以往研究[8~9,12],通常取m=2,r = 0.2SD,SD 為原始序列的標準差。利用Matlab R2009a軟件進行編程計算。
本文以阿城站為例,對其1961—2010年逐月實測徑流序列進行復雜性分析。阿城站位于黑龍江省阿城市阿城鎮(zhèn),是松花江支流阿什河下游控制站,集水面積2 313 km2。阿什河屬于半山區(qū)性河流,上游為山區(qū),平均坡降1/500 左右,中游坡降為1/1 000 ~1/1 500,下游河道進入平原,坡降為1/1 500 ~1/1 800。此流域地處松花江中游,地形東南高、西北低,上游山巒起伏,生長著大片的次生林;中下游為丘陵平源,土壤肥沃,植被茂盛。本區(qū)屬中寒溫帶大陸性季風氣候區(qū),春季回暖快,多風少雨,易發(fā)生春旱;夏季溫熱,暴雨集中,易發(fā)生洪水;秋季涼爽降溫快,霜凍寒流來的早;冬季寒冷干燥。本站多年平均徑流量為4.537 ×108m3;多年平均年降水量512.0 mm,主要集中在汛期6—9月份,占年降水量的82.0%。
為描述阿城站月徑流序列復雜性的動態(tài)變化,以其1961—2010年的逐月實測徑流時間序列為樣本,以120個月的數(shù)據(jù)長度為滑動窗口,沿徑流時間序列H( t) ( t=1,2,…,N) 移動,滑動步長為1,直至序列結(jié)束。對每個窗口內(nèi)的數(shù)據(jù)計算其樣本熵值,根據(jù)計算值繪制徑流序列樣本熵值的動態(tài)變化曲線,見圖1。

3) 給定閾值r,對每一個1≤i≤N -m 值,統(tǒng)計d[Xm( i) ,Xm( j) ]<r 的數(shù)目( 模板匹配數(shù)) 及此數(shù)目與距離總數(shù)N-m-1 的比值,記作Bim( r) :

圖1 阿城站月徑流序列樣本熵動態(tài)變化曲線
由圖1 可知,在1961—2010年,阿城站月徑流序列SampEn 值具有較為顯著的周期性變化特征,開始呈現(xiàn)逐漸減小的趨勢,然后呈現(xiàn)明顯增加的趨勢,造成這種發(fā)展趨勢的可能原因有降水、下墊面條件、蒸發(fā)、氣溫、土壤蓄水狀況等水文地質(zhì)條件及人類活動如修建排水工程及河道護坡、修建水庫及水閘等。同時可知,阿城站月徑流序列復雜度具有明顯的隨時間變化而變化的特征,SampEn 值由最初的平穩(wěn)波動狀態(tài)轉(zhuǎn)變是逐漸持續(xù)減小的總體發(fā)展態(tài)勢,這說明在流域自然環(huán)境相對穩(wěn)定的狀態(tài)下,人類生產(chǎn)活動對水文系統(tǒng)的影響程度較大。
由圖2 可知,阿城站月徑流序列與其相對應的SampEn 值之間存在著某種相關關系。具體來說,在以120個月為滑動窗口的條件下,阿城站月徑流量滑動平均值與其相對應的SampEn 值之間存在著負相關關系,而且在兩者波峰處存在著非常明顯的對應關系。若月徑流量滑動平均值增加,則與其相對應的SampEn 值就會減小,說明徑流在向有序方向發(fā)展,徑流序列復雜性正在減弱;反之,若月徑流量滑動平均值減小,則與其相對應的SampEn 值就會增大,說明徑流在向無序方向發(fā)展,徑流序列復雜性正在加強。

圖2 阿城站月徑流量與其樣本熵關系曲線
本文運用樣本熵方法分析了阿城站月徑流序列的復雜性特征,初步得出以下結(jié)論:
1) 樣本熵的動態(tài)分析結(jié)果表明,阿城站月徑流序列的復雜性具有明顯的時間差異,SampEn 值的波動較大且具有顯著的周期性特征,這反映了不同時期的水利工程建設、護坡建設等人類生產(chǎn)活動對水文系統(tǒng)動力學結(jié)構(gòu)產(chǎn)生了重大影響,從而導致了月徑流序列復雜性的明顯波動。
2) 阿城站月徑流量滑動平均值與其相對應的SampEn 值之間存在著負相關關系,而且兩者波峰處的對應關系非常明顯,這種變化特征及對應關系對于當?shù)貜搅餍蛄械母呔阮A測研究具有重大意義。
3) 樣本熵的計算不依賴數(shù)據(jù)長度,只需較少的數(shù)據(jù)就能反映時間序列的非線性特征,而且對缺失數(shù)據(jù)不敏感,算法方便易行,計算結(jié)果具有較高一致性,為徑流序列及其它水文時間序列的復雜性研究提供了一種新的方法。
[1]馮國章,宋松柏,李佩成. 水文系統(tǒng)復雜性的統(tǒng)計測度[J]. 水利學報,1998(11) :76 -81.
[2]姜林,艾波,漆濤. 分形理論在軟件復雜度中的應用[J].計算機應用,2010,30(10) :2730 -2734.
[3]劉函林,黃華,劉壙彬. 腦電信號分析的實用符號動力學方法研究[J]. 生物醫(yī)學工程學雜志,2010,27(2):407-410.
[4]曹蕾,于國榮. 長江日流量時間序列混沌特性研究——相空間嵌入維數(shù)和嵌入滯時的聯(lián)合確定[J]. 水利科技與經(jīng)濟,2011,17(7) :9 -11.
[5]童勤業(yè),孔軍,徐京華. 腦電復雜性分析的新方法[J]. 中國生物醫(yī)學工程學報,1998,17(3) :222 -226.
[6]李素蘭. 腦電時間序列分析的非線性動力學方法[D].杭州:浙江大學,1998.
[7]Richman J S,Moorman J R. Physiological time-series analysis using approximate entropy and sample entropy[J]. Am J Physiol Heart Circ Physiol,2000,278:H2039 -H2049.
[8]白冬梅,邱天爽,李小兵. 樣本熵及在腦電癲癇檢測中的應用[J]. 生物醫(yī)學工程學雜志,2007,24(1) :200 -205.
[9]彭濤,陳曉宏,莊承彬. 基于樣本熵的東江月徑流序列復雜性分析[J]. 生態(tài)環(huán)境學報,2009,18(4) :1379 -1382.
[10]葛家怡,周鵬,趙欣. 睡眠腦電時間序列的非線性樣本熵研究[J]. 電子器件,2008,31(3) :972 -975.
[11]Alcaraz R,Rieta J J.A review on sample entropy applications for the non -invasive analysis of atrial fibrillation electrocardiograms[J].Biomedical Signal Processing and Control,2010,5(1) :1 -14.
[12]蘇文勝,王奉濤,朱泓等. 基于小波包樣本熵的滾動軸承故障特征提取[J]. 振動、測試與診斷,2011,31( 2) :162 -166,263.