基于相對關聯距離熵的水文變異點識別*

杜 建，陳曉宏，陳志和

(中山大學 水資源與環(huán)境研究中心∥廣東省高校水循環(huán)和水安全重點實驗室，廣東 廣州 510275)

受人類活動和氣候變化的影響，水文特征的時間變異是客觀存在的，認識一個區(qū)域水文要素變化規(guī)律，修正區(qū)域水文時間序列長序列變化的非一致性，從而可為洪水預報、水資源開發(fā)以及水利工程設計施工和管理提供科學依據。目前在水文時間序列中用于確定可能變異點的方法，如有序聚類法、Mann-Kendall檢驗法、秩和檢驗法、R/S分析法［1］、貝葉斯方法［2］、啟發(fā)式分割算法［3］和信息測度模型［4］等，這些方法大多數是從概率統(tǒng)計方面來研究，卻很少有從系統(tǒng)內部的特征變化來識別水文變異點。

由于降雨、徑流等水文過程受眾多因素的影響而造成強烈的變異性，表現出偽隨機特征，致使傳統(tǒng)確定性數學模型對這些水文過程的模擬遇到了很大的困難，而混沌理論的出現為研究這種高度復雜的系統(tǒng)提供了新的思路［5］，所以從20世紀80年代開始，混沌理論和方法在降雨、徑流、洪澇、地下水、水質以及降雨－徑流分析及預測等水文水資源領域的應用研究快速發(fā)展［6－7］。

水文變異點在水文系統(tǒng)狀態(tài)發(fā)生改變時表現出來，如何提取系統(tǒng)內部狀態(tài)變化的信息成為識別水文變異點的關鍵。由時間序列重構相空間的過程可知:重構相空間中點的軌跡是原時間序列動力性質的一種反映。如果系統(tǒng)的狀態(tài)發(fā)生了變化，軌跡的形狀及大小必然發(fā)生相應的變化。如能找到反應軌跡的形狀及大小變化的參數，則可從系統(tǒng)內部特征的變化找到水文變異點。文獻 ［8］曾應用相對關聯距離熵來描述統(tǒng)重構后相空間中軌跡形狀及大小的變化，將其應用于內燃機氣門漏氣故障診斷研究中，取得了較好的效果。

本文嘗試應用混沌理論對水文變異點進行診斷，在重構水文時間序列的基礎上，根據反映系統(tǒng)狀態(tài)變化的相對關聯距離熵的變化來對水文序列時間變異點進行識別，以期對水文變異分析做進一步探討。

1 水文時間序列的相空間重構

混沌學從理論走向應用很大程度上取決于Packard 和 Takens提出的相空間重構理論［9－10］。相空間重構把系統(tǒng)中每一點的實質性信息都抽象出來，揭示了其它方式看不見的運動模式，就象紅外線照片可以揭示非感官所及的圖形和細節(jié)一樣，這對于研究復雜的水文系統(tǒng)是一個很有力的工具［11］。相空間重構能保持水文時間時序所對應的原動力系統(tǒng)內在結構的幾何不變性，系統(tǒng)中任一分量的演化都是由與之相互作用著的其它分量所決定的，因此這些相關分量的信息就隱含在任一分量的發(fā)展過程之中［12］。時間延遲技術是時間序列相空間重構最常用的方法，即為了重構一個等價的狀態(tài)空間，只需考察一個分量，并將它在某些固定的延遲時間間隔點上的測量作為新維處理，它們確定了某個多維狀態(tài)空間中的一點，重復這一過程并測量相對于不同時間的各延遲量，就可以產生出許多這樣的點，它可以將系統(tǒng)的許多性質保存下來，即用系統(tǒng)的一個觀察量可以重構出原動力系統(tǒng)模型。

水文系統(tǒng)可以看作是具有n個變量的動力系統(tǒng)［13］，用下述方程表述:

系統(tǒng)的時間演變的狀態(tài)空間可由坐標x(t)加上其(n－1)階導數所構成的n維相空間來表示，即

Ruelle提出用離散的時間序列x(t)和它的 (m－1)時滯位移構建一個新的m維相空間 (即嵌入相空間)，以代替這種連續(xù)變量x(t)和它的導數所反映的系統(tǒng)狀態(tài)空間，即

其中，τ為滯時，也稱延遲時間。對于某一水文站可觀測的徑流時間序列x1，x2，…，xn，選定τ和m，則其重構嵌入相空間可表示為:

其中，N=n-(m-1)τ

從水文徑流時間序列重構相空間的過程可知:由時間序列重構的相空間中點的軌跡是原時間序列動力性質的一種反映。軌跡的形狀及大小均反映了原系統(tǒng)的動力性質，如果系統(tǒng)的狀態(tài)發(fā)生了變化，那么其外在表現必然發(fā)生相應的變化。

2 徑流序列的混沌性識別

Lyapunov指數是一個譜系 λi(i=1，2，…，m－1)，若譜系中只有最大者為正值，則系統(tǒng)為一維混沌;若譜系中具有兩個或兩個以上正的Lyapunov指數，則系統(tǒng)為多維混沌或超混沌，否則系統(tǒng)為非混沌。因此，一般用最大Lyapunov指數來進行混沌分析。最大Lyapunov指數是判斷混沌存在的一個重要依據，對于時間序列 x1，x2，…，xn，其求取方法如下［13－14］:先按上節(jié)所述方法進行相空間重構，然后根據式 (5)求得歐氏空間意義上Yi的最近鄰點及兩點間距離，最后根據式 (6)求最大Lyapunov指數λ。

其中N=n-(m-1)τ。

根據以上方法，求得東江流域博羅站月徑流系列的最大Lyapunov指數為0.245，從而說明了該水文時間序列存在混沌特性。

3 基于相對關聯距離熵的水文變異識別方法

變異點是水文序列發(fā)生急劇變化的一種表現形式，當水文系統(tǒng)從一種狀態(tài)轉化為另一種狀態(tài)時表現出來，水文系統(tǒng)狀態(tài)的改變將導致水文時間序列的隨機性和混沌強弱程度發(fā)生變化，進而使得重構相空間中相軌跡點的分布情況發(fā)生改變，從而相空間點分布的疏密程度也會發(fā)生改變。如何從重構的軌跡中提取反應水文狀態(tài)變化的信息，將成為水文變異診斷成功與否的一個關鍵環(huán)節(jié)，而相對關聯距離熵正是體現了這一變化的特征參數。

相對關聯距離熵的計算方法如下［8］:

1)有水文時間序列 {xi}(i=1，2，…，n)，以嵌入維數m和延遲時間間隔τ重構相空間，按照上述重構相空間的方法得到重構相空間中的點{Yi}(i=1，2… N，N=n-(m － 1)τ)。

2)求相空間中系統(tǒng)相軌跡點對之間的距離

3)求相空間中系統(tǒng)相軌跡點對之間的距離的最大值

4)對相空間中系統(tǒng)相軌跡點對之間的距離dij進行歸一化處理

5)定義絕對關聯距離熵

由CH的表達式可知，定義的關聯距離熵反映了系統(tǒng)相軌跡上點在相空間中分布的疏密程度，是系統(tǒng)在重構相空間中軌道性質的一種反映。對于變點，是觀測值的分布按某種規(guī)律變化，而在變點處發(fā)生突變，改換成另一種規(guī)律，在水文時間序列重構相空間的表現是相軌跡點的疏密程度發(fā)生了改變，而CH是一個能很好反映這種變化的參數。但是CH反映的是一種絕對大小，與相點個數有關，由此可知這與數據的個數n以及相空間選取的嵌入維數m和延遲時間間隔τ有關，而m和τ的選擇是人為控制的，不具有客觀性，不利于水文要素的變異點識別，在此基礎上提出相對關聯距離熵的概念。由于CH反映了動力系統(tǒng)在重構相空間中軌跡上點分布的疏密程度，點分布越均勻則CH的絕對值越大，而點分布越集中則CH的絕對值越小。因此當軌跡上點均勻分布時CH值應為最大，令其為CHmax，在此基礎上，對CH進行歸一化處理，這樣得到的熵稱為相對關聯距離熵。

經處理后得到的Hd是反映系統(tǒng)性質的參量，與N無關，Hd是一個無量綱的量，有利于水文變異的診斷，系統(tǒng)性質的改變將導致重構相空間中軌跡點分布疏密程度的不同，從而CH不同，相應的Hd不同，Hd是反映系統(tǒng)性質的參量。

因相空間重構對數據長度有一定要求，假設變異的點在j年，其中j={8，9，…，n}，取j時刻前的時間序列xj，在重構徑流相空間的基礎上，按照以上方法求出時間序列 xj的 Hd(xj)，計算 Hd(xj)的增值比d Hd(xj)，即

找出dHd(xj)最大的那一點，則視這一點為變異的最大點，相應便可診斷出變異的年份。

4 實例研究

如果流域受到自然及人類大規(guī)模活動影響，從而流域的自然條件、下墊面條件會發(fā)生改變，這些會使該流域水文時間序列的特征數相對關聯距離熵發(fā)生變化。因此，可以通過分析流域水文時間序列的相對關聯距離熵的變化，找出該時間序列的變異點。本文以廣東省東江流域博羅站1956－2005年的逐月徑流量時間序列為例來驗證本方法的有效性。

假設變異的點為 j={8，9，10，…，n}，取 j時刻前的時間序列xj，在重構徑流相空間的基礎上，按照以上方法求出時間序列xj的Hd(xj)，其中m=4，τ =8，計算Hd(xj)的增值比

d Hd(xj)，找出dHd(xj)最大的那一點，則視這一點為變異的最大點，相應便可診斷出變異的年份，分析的計算結果如表1所示。

表1 博羅站徑流序列相對關聯熵與相對關聯距離熵增值比Table 1 Relative correlation distance entropy and increment ratio of Hd(xj)of runoff series at Boluo Station

從表1可見，相對關聯距離熵Hd在1.09－1.14之間，變化的幅度較小，這主要相對關聯距離熵反映的是相空間主對角線方向點分布的疏密程度，使得重構相空間中點的疏密程度改變不是特別明顯，但是從表1可以看出，1963－1972年的Hd在1.09～1.11之間交替變化，1973－1979年的Hd在1.12左右，1980－1983年的Hd在1.13左右，1984－2005年的Hd在1.14左右，由此可見1963－1972年的變化較其他年份變化劇烈，在1963－1972年徑流序列發(fā)生變異的可能性最大，這主要因為在60－70年代東江流域開始進行大規(guī)模的水利建設，在此期間新豐江水庫、東深供水等水利工程相繼建成，這些人類活動的影響可能導致徑流序列的狀態(tài)發(fā)生變化。從圖1年增值比d Hd(xj)的變化趨勢看，1963－1972變化比較劇烈，1973－1985年變化次之，1986－2002年變化緩慢，從2003年開始又有變化較大的趨勢，這與東江流域的主要水利工程的建設的過程是一致的，東江流域的三大水庫中新豐江水庫 (1962年建成)、楓樹壩水庫 (1973年建成)、白盆珠水庫 (1985年建成)的相繼發(fā)揮作用，使得博羅站的徑流發(fā)生變化。從圖1中的年增值比dHd(xj)的變化可以明顯的看出博羅站徑流量在1972、1982年發(fā)生較大變化1972年明顯變化最大，所以1972年可以診斷為東江流域博羅站徑流量的變異點，這與文獻 ［15］東江流域水文特征在1972年發(fā)生變異的結論相一致。

圖1 博羅站徑流序列Hd(xj)的增值比的變化Fig.1 Change of the increment ratio of Hd(xj)of runoff series

為了近一步探討徑流序列變異的原因，分析博羅站年徑流量與博羅站面雨量的相關關系 (圖2)，在變化前 (1956－1972年)，降雨徑流的相關系數為0.969，變化后 (1973－2000年)降雨徑流的相關系數為0.937，而新豐江水庫與博羅站年徑流的相關關系 (如圖3所示)在變化前 (1961－1972年)的相關系數為0.895，變化后 (1973－2005年)的相關系數為0.933，由此可看出，博羅站徑流量在受降雨影響的同時，越來越受新豐江水庫的調節(jié)作用的影響。從圖1還可以看出除1972年變化較大外，1982年也變化較大。為了探討1982年變化較大的原因，分析1973－1982年博羅站與楓樹壩水庫年徑流量的相關系數為0.789(圖4)，而在1982變化后的相關系數為0.856，而此時 (1973－1982年)博羅站年徑流與面雨量的相關性0.966，而在1982變化后博羅站年徑流與面雨量的相關性為0.951。

同樣表明，博羅站的徑流受到楓樹壩水庫的調節(jié)作用。1972年變化比1982年變化劇烈是因為新豐江是東江流域最大的水庫，總庫容為138.96億m3，具有多年調節(jié)能力，而楓樹壩水庫的總庫容為19.32億m3，相比與新豐江水庫，楓樹壩水庫的調節(jié)性能小的多。綜上所述，東江流域降雨徑流變化過程關系密切，降雨對徑流的影響有隨著年代推移而減弱的趨勢。這說明，變化環(huán)境下，徑流主要受降雨影響的同時，愈來愈受到水利工程的影響，在徑流演變過程中的作用越來越不容忽視。

為了與傳統(tǒng)的水文變異方法比較，對東江博羅站的徑流序列做Mann-Kendall突變點檢驗，如圖5所示。

圖5 博羅年徑流量Mann-kendall統(tǒng)計量曲線Fig.5 M-K trends of runoff series at Boluo Station

從曲線的交點來看，有多處交點，位置在1963、1966、1972、2002年，但是無法判斷哪些交點是比較顯著的變點。Mann-Kendall方法對于趨勢變化明顯的單變點情形既簡單而且直觀有效，但若對多變點情形誤用可能得出錯誤結論，必需在先驗的了解數據規(guī)律后，對多變點情形分段應用該法。而基于相對關聯距離熵的水文變異點的識別方法可以清晰的判斷出比較顯著的變點。

關聯距離熵反映了重構相空間中吸引子軌道的分布及形狀，體現了系統(tǒng)的本質特征，從而具有較好的對系統(tǒng)性質表達能力。根據相對關聯距離熵的變化來反應系統(tǒng)平穩(wěn)性和隨機性的變化，水文變異點是在從一種狀態(tài)轉化為另一種狀態(tài)時表現出來，在相空間中的體現是相空間點軌跡點的疏密程度發(fā)生變化，而相對關聯距離熵正是體現了這種變化?；陉P聯距離熵的水文時序變異點的識別從系統(tǒng)的動力特征的角度，探索從系統(tǒng)內部的特征的變化來辨識水文系統(tǒng)的狀態(tài)的跳躍，從而識別出水文變異點。

基于相對關聯距離熵的水文變異點識別方法與其他水文變異點識別方法相比，能清晰的識別多個變點的情況，計算簡單，從系統(tǒng)的動力特征的角度，識別了水文變異點。從關聯距離熵的定義可知，關聯距離熵反映的是系統(tǒng)在相空間中沿主對角線方向點分布的疏密程度，但無法反映系統(tǒng)在相空間中其它方向的分布情況，使得在應用中關聯距離熵的變化不是特別明顯，如何從重構的吸引子中提取更全面、能反映系統(tǒng)的參數需要進一步探討。

5 結論

本文以東江博羅站50年月徑流量序列為例子，應用混沌理論中相對關聯距離熵對博羅站月徑流序列進行水文變異識別。首先將月徑流量構成多維相空間，展示系統(tǒng)的動力特征，并基于此計算序列的相對關聯距離熵，分析相對關聯距離熵的變化來辨識系統(tǒng)的變化，從而識別出水文變異點。本方法從系統(tǒng)的動力特征的角度，探索從系統(tǒng)內部的特征的變化來識別水文變異點，作為水文變異點識別的新方法是可行的。水文系統(tǒng)狀態(tài)的改變，在相空間點的疏密程度也會發(fā)生改變，但是變化程度在相空間中沿主對角線方向反映不是特別明顯，如何從重構的吸引子中提取更全面、能反映系統(tǒng)的參數尚需要進一步探討研究。

水文時間序列既有確定性的部分，又有隨機性的部分，探索水文的變化規(guī)律，把握水文變化特征，更可靠的的預測未來，選用相關研究方法是非常重要的。本文將混沌理論中相對關聯距離熵引入到水文變異點的識別是一次有益的嘗試，尚需要進一步研究和完善。

［1］王孝禮，胡寶清，夏軍.水文時序趨勢與變異點的R/S分析法［J］.武漢大學學報:工學版，2002，35(2):10 －12.

［2］熊立華，周芬，肖義等.水文時間序列變點分析的貝葉斯方法［J］.水電能源科學，2003，21(4):39－41.

［3］陳廣才，謝平.基于啟發(fā)式分割算法的水文變異分析研究［J］.中山大學學報:自然科學版.2008，47(5):122－125.

［4］肖宜，穆宏強，申明亮.差異信息理論在水文時間序列變異點診斷中的應用［J］.中國農村水利水電，2002，11:28－30.

［5］丁濤，周惠成，黃健輝.混沌水文時間序列區(qū)間預測研究［J］.水利學報，2004(12):15 － 20.

［6］王紅瑞，宋宇，劉昌明等.混沌理論及在水科學中的應用與存在的問題［J］.水科學進展，2004(5):400－407.

［7］李榮峰，沈冰，張金凱.基于相空間重構的水文自記憶預測模型［J］.水利學報，2006(5):583 －588.

［8］夏勇，趙紅.基于關聯距離熵診斷方法研究.振動與沖擊［J］，2003，22(2):76 －77.

［9］ROSSLER O E.An Equation for Continuous Chaos［J］.Phys Lett，1976，57:13 －15.

［10］PACKARD N H，CRUTCHF IELD J P，FARMER J D，et al.Geometry from a time series［J］.Phys Rev Lett，1980，45(9):712 － 716.

［11］趙永龍，丁晶，鄧育仁.混沌分析在水文預測中的應用和展望［J］.水科學進展，1998(6):181－186.

［12］呂金虎，陸君安，陳士華.混沌時間序列分析及其應用［M］.武漢:武漢大學出版社.2002.

［13］石教智，陳曉宏.林汝顏.東江流域降水時問序列的混沌特征分析［J］.中山大學學報:自然科學版，2006，45(4):111 －115.

［14］LAI Y C，LERNER D.Efective scaling regime for corn.puting the correlating dimension for chaotic time series［J］.Physica D，1998(11):10 －18.

［15］劉德地.變化環(huán)境下水資源優(yōu)化配置［D］.廣州:中山大學，2008.