13

(1.北京師范大學(xué) 水科學(xué)研究院，北京 100875； 2.北京師范大學(xué) 城市水循環(huán)與海綿城市技術(shù)北京市重點(diǎn)實(shí)驗(yàn)室，北京 100875； 3.北京師范大學(xué) 地下水污染控制與修復(fù)教育部工程研究中心，北京 100875)

1 研究背景

混沌是具有確定性的動(dòng)力學(xué)系統(tǒng)受到非線性變量的影響而產(chǎn)生的一種沒有規(guī)則的、貌似隨機(jī)的現(xiàn)象?；煦缪芯渴加?0世紀(jì)70年代左右，美國(guó)氣象學(xué)家Lorenz[1]發(fā)表了《確定性非周期流》，發(fā)現(xiàn)在確定性的非線性方程中存在混沌運(yùn)動(dòng)，隨后又在不同的動(dòng)力學(xué)方程中也研究出了混沌現(xiàn)象。水文學(xué)領(lǐng)域的混沌研究起步于20世紀(jì)80年代后期，Hense[2]在研究月降雨量序列的關(guān)聯(lián)維數(shù)時(shí)，同樣發(fā)現(xiàn)了混沌現(xiàn)象存在于水文系統(tǒng)中。

30年來，許多專家學(xué)者利用混沌理論對(duì)水文學(xué)領(lǐng)域的各個(gè)方面進(jìn)行了初步的研究，取得了豐碩的成果。迄今為止，水文學(xué)領(lǐng)域的混沌研究仍處于探索階段，還有許多的問題有待于研究。Jayawardena 和Lai[3]于1994年對(duì)香港3個(gè)站點(diǎn)的日降雨序列和2個(gè)站點(diǎn)的日徑流數(shù)據(jù)進(jìn)行了研究分析，為香港的日降雨序列和日徑流序列存在混沌性提供了比較可信的證據(jù)，雖然研究中非線性預(yù)測(cè)方法的精度不高，但對(duì)比傳統(tǒng)的自回歸移動(dòng)平均法具有獨(dú)特的優(yōu)勢(shì)。Sivakumar等[4-7]研究了新加坡6個(gè)觀測(cè)站點(diǎn)不同序列長(zhǎng)度的日降雨數(shù)據(jù),證明了降水過程混沌性的存在，并使用代理數(shù)據(jù)法判定了日降雨時(shí)間序列是非線性的，還對(duì)以往研究中預(yù)測(cè)精度較低的原因進(jìn)行了分析。Waelbroeck等[8]研究了墨西哥特拉克斯卡拉地區(qū)的降雨資料，發(fā)現(xiàn)利用累計(jì)10 d的降雨資料進(jìn)行暴雨預(yù)測(cè)更有效。

對(duì)于降雨過程，特別是暴雨過程，學(xué)者普遍認(rèn)為存在混沌性。王紅瑞等[9]認(rèn)為水文系統(tǒng)中混沌的提出是正確的，但對(duì)于受噪聲干擾的時(shí)間序列而言，現(xiàn)有方法所得到的結(jié)論僅可表明其可能存在混沌特征。王文等[10]通過對(duì)水文時(shí)間序列的最大可預(yù)報(bào)時(shí)間長(zhǎng)度、混沌特征參數(shù)估計(jì)和數(shù)據(jù)量大小等問題的討論，指出目前水文系統(tǒng)中混沌的研究尚缺少謹(jǐn)慎合理的結(jié)論和解釋。一系列基于混沌特性識(shí)別基礎(chǔ)的預(yù)測(cè)方法的出現(xiàn)表明，混沌理論的應(yīng)用十分有必要。Domenico等[11]以及Wang等[12]認(rèn)為徑流序列中存在隨機(jī)成分的構(gòu)成以及對(duì)初始條件的敏感，即混沌特性對(duì)徑流的動(dòng)態(tài)變化起決定作用，在這種情況下徑流預(yù)測(cè)應(yīng)該更側(cè)重簡(jiǎn)化模型并考慮混沌理論的應(yīng)用。Zounemat-Kermani[13]對(duì)日徑流序列和月徑流序列的研究表明其都有混沌性存在的可能，在中期徑流預(yù)報(bào)中非線性局部逼近法要優(yōu)于ARIMA法，在短期預(yù)報(bào)中2種方法都是可以接受的。

安康地處漢江上游，是南水北調(diào)中線工程水源地重要的水量來源，其降水量的變化在一定程度上影響著漢江徑流量的變化。相關(guān)研究表明，漢江上游年降水量呈不顯著減少[14]。尤其是進(jìn)入21世紀(jì)以來，隨著氣候變化和人類活動(dòng)的雙重影響，漢江上游典型地區(qū)安康極端降水等問題日益突出[15]?；诖?，本文利用混沌理論對(duì)漢江上游安康降水量進(jìn)行特征分析，研究其降水系統(tǒng)動(dòng)態(tài)變化作用下的徑流預(yù)報(bào)，研究結(jié)果可為南水北調(diào)中線工程水量變化研究提供及時(shí)的數(shù)據(jù)和理論支撐。

2 研究區(qū)概況

安康市地處陜西省東南部，境內(nèi)河流屬長(zhǎng)江流域漢江水系，處在漢江上游。河網(wǎng)密度為1.43 km/km2，年平均徑流量為10.65億m3。過境客水145.72億m3，水資源總量共計(jì)252.29億m3。漢江因其水量大、水質(zhì)好而成為國(guó)家南水北調(diào)中線工程水源涵養(yǎng)地[16]。其中安康水文站是漢江上游控制站，位于陜西省安康市漢濱區(qū)西南城郊(109°12′E，32°25′48″N)，觀測(cè)場(chǎng)海拔290.8 m。

3 降水特征分析

根據(jù)1929—2014年安康站的降水?dāng)?shù)據(jù)可知，其平均年降水量為788.76 mm，1938年降水量最高，達(dá)到1 392 mm。1929—2014年中年降水量<500 mm的有5 a，均出現(xiàn)在1950年以前。經(jīng)過分析相關(guān)數(shù)據(jù)可知，1950年以前數(shù)據(jù)缺失較多，準(zhǔn)確性較差，不能較好地表現(xiàn)安康站降水特征的變化，故單獨(dú)統(tǒng)計(jì)分析1950—2014年安康站的降水?dāng)?shù)據(jù)，見圖1，其平均年降水量為838.27 mm。其中，2010年降水量最高，為1 232 mm，1999年降水量最低，為526 mm。1950—2014年中共計(jì)13 a全年降水量>1 000 mm，共計(jì)3 a全年降水量<600 mm，分別是1957年、1959年、1999年。

圖1 1950—2014年安康站歷年逐月降水量

通過圖2可知，安康站降水主要集中在每年6—9月份，這4個(gè)月的降水總量占全年的近60%，7月份最多。其中2010年7月份單月降水量達(dá)到455 mm，平均每月146 mm，最低為1959年7月份降水，僅21 mm。

圖2 1950—2014年安康站各月份平均降水量

若將降水過程考慮為隨機(jī)過程，則降水時(shí)間序列{xn}可以寫為

記：“四個(gè)全面“中，全面依法治國(guó)是很重要的一點(diǎn)；但現(xiàn)在，有法不依、執(zhí)法不嚴(yán)的情況時(shí)有發(fā)生。您認(rèn)為全面依法治國(guó)還存在哪些困境？該如何破局？

xn=fn+zn。

(1)

式中：噪聲{zn}是一個(gè)離散的隨機(jī)過程且期望E(zn)=0；fn是趨勢(shì)函數(shù)f(t)在tn時(shí)刻的趨勢(shì)值。假定噪聲{zn}與趨勢(shì)fn無關(guān)，趨勢(shì)函數(shù)f(t)具有變化較慢的單調(diào)性。對(duì)于趨勢(shì)函數(shù)，并非考慮時(shí)間間隔而進(jìn)行提取，而是結(jié)合時(shí)間序列值域變化的空間進(jìn)行分段的線性近似所得[17]。

對(duì)1950—2014年降水量數(shù)據(jù)，利用文獻(xiàn)[17]所提供的去趨勢(shì)算法，得到如圖3虛線所示的趨勢(shì)結(jié)果，可以看出安康站1950—2014年歷年總降水量以及趨勢(shì)變化情況整體上較為波動(dòng)，可以分為4個(gè)大周期：1950—1955年呈現(xiàn)上升趨勢(shì)，1955—1960年呈現(xiàn)下降趨勢(shì)；1961—1965年呈現(xiàn)上升趨勢(shì)，1966—1972年呈現(xiàn)下降趨勢(shì)；1973—1984年呈現(xiàn)上升趨勢(shì)，1985—1999年呈現(xiàn)下降趨勢(shì)；2000—2010年呈現(xiàn)上升趨勢(shì)，最近幾年略有回落。

圖3 安康站1950—2014年降水量以及趨勢(shì)變化

4 混沌特性識(shí)別與可預(yù)測(cè)時(shí)長(zhǎng)分析

降水過程等水文現(xiàn)象不能用微分方程組來描述，無法直接得到其相空間，故研究降水時(shí)間序列的混沌性首先需要重構(gòu)相空間。混沌特性分析是在確定非線性系統(tǒng){x(t)}的延遲時(shí)間τ、嵌入維度m后重構(gòu)相空間X(t)，對(duì)X(t)進(jìn)行混沌特征識(shí)別，并判斷系統(tǒng)是否具有混沌特征，進(jìn)而對(duì)其可預(yù)報(bào)性和可預(yù)報(bào)時(shí)間尺度進(jìn)行分析。

4.1 相空間重構(gòu)

相空間重構(gòu)是水文線性混沌性識(shí)別和非線性預(yù)測(cè)的基礎(chǔ)，Grassberger和Procaccia[18]提出了用延遲坐標(biāo)法對(duì)混沌時(shí)間序列{x(t)}(t=1,2,…，N)(N指時(shí)間序列的長(zhǎng)度)進(jìn)行相空間重構(gòu)，相空間中的點(diǎn)可以表示為

X(t)=[x(t),x(t-τ),…,x(t-(m-1)τ)]T。

(2)

式中：m為嵌入維數(shù)；τ為延遲時(shí)間，本研究中單位為月。

在重構(gòu)相空間時(shí)，參數(shù)嵌入維數(shù)m和延遲時(shí)間τ的選取很關(guān)鍵。目前選取τ的方法有多種，包括自相關(guān)函數(shù)法、廣義相關(guān)函數(shù)法等。嵌入維數(shù)m的選擇方法也有多種，常用的有飽和關(guān)聯(lián)維數(shù)法、偽最近鄰點(diǎn)法等。

有些研究人員認(rèn)為m和τ是有關(guān)聯(lián)的[19]，對(duì)于某個(gè)實(shí)際的時(shí)間序列，其嵌入窗寬τw=(m-1)τ相對(duì)固定，近似于相軌線的平均周期。由此人們提出了一些計(jì)算m和τ的聯(lián)合算法，如時(shí)間窗口法、C-C(Corelation integral)法等。

4.1.1 通過自相關(guān)函數(shù)法確定延遲時(shí)間τ

時(shí)間序列{x(t)}(t=1,2,…，N)的自相關(guān)系數(shù)rτ為

(3)

式中μ為時(shí)間序列的均值。以自相關(guān)函數(shù)首次通過0點(diǎn)時(shí)對(duì)應(yīng)的τ值作為延遲時(shí)間τ。

4.1.2 通過G-P算法求得關(guān)聯(lián)維數(shù)確定最佳嵌入維數(shù)

基于已經(jīng)完成的重構(gòu)相空間X(t)(t=1,2,…,T)，若給定一個(gè)鄰域半徑r，則可以定義關(guān)聯(lián)積分Cm(r)為

(4)

式中：i、j是不同的t值；T是重構(gòu)相空間序列的長(zhǎng)度；總的相點(diǎn)數(shù)M=t-(m-1)τ；H為階躍函數(shù)，當(dāng)x>0時(shí)，H(x) = 1，當(dāng)x≤0時(shí)，H(x) = 0；‖·‖為范數(shù)。

通常選取嵌入維數(shù)m≥2Dm+1，其中Dm是動(dòng)力系統(tǒng)的關(guān)聯(lián)維數(shù)[20]。飽和關(guān)聯(lián)維數(shù)法，即G-P算法基于已經(jīng)完成的重構(gòu)相空間X(t)，通過求得關(guān)聯(lián)積分來確定飽和關(guān)聯(lián)維數(shù)。

4.1.3 利用C-C法聯(lián)合求延遲時(shí)間τ和嵌入維數(shù)m

4.2 混沌特性識(shí)別

Wolf法是基于相軌跡、相平面、相體積等演化來估計(jì)Lyapunov指數(shù)的，對(duì)于重構(gòu)相空間X(t)，取初始點(diǎn)X(t0)，并設(shè)其與最鄰近點(diǎn)X0(t0)的距離為L(zhǎng)0。直到t1時(shí)刻，兩點(diǎn)間距離超過某一規(guī)定值ε>0，即L0=|X(t1)-X0(t1)|>ε，則保留X(t1)，并在X(t1)鄰近另外找一點(diǎn)X1(t1)，使L1=|X(t1)-X1(t1)|<ε，并且與X1(t1)的夾角要盡可能小。重復(fù)上述過程，直至X(t)到達(dá)時(shí)間序列的終點(diǎn)N，這時(shí)追蹤演化過程總的迭代次數(shù)為M，則最大Lyapunov指數(shù)λi為

(5)

4.3 可預(yù)測(cè)時(shí)間尺度

目前較常使用最大Lyapunov 指數(shù)λi的倒數(shù)來定義系統(tǒng)的最大可能預(yù)報(bào)時(shí)間尺度Tf，即

(6)

圖4 rτ-τ關(guān)系

5 運(yùn)算結(jié)果及討論

圖5 C-C方法計(jì)算的延遲時(shí)間和時(shí)間窗

設(shè)定一系列m，取值從1～21變化，計(jì)算所對(duì)應(yīng)的Dm。對(duì)于不同的嵌入維數(shù)m取值，作一簇lnCm(r)-lnr曲線。根據(jù)G-P關(guān)聯(lián)積分法，除去斜率為0和∞的直線區(qū)段，確定其無標(biāo)度區(qū)的最佳擬合直線的斜率為關(guān)聯(lián)維數(shù)Dm。假若隨著嵌入維數(shù)m的增加，關(guān)聯(lián)維數(shù)Dm趨于飽和，此刻的Dm稱為飽和關(guān)聯(lián)維數(shù)，作為奇怪吸引子維數(shù)的一種估計(jì)。對(duì)安康站1950—2014年的降水時(shí)間序列數(shù)據(jù)進(jìn)行處理后所得到的關(guān)聯(lián)維數(shù)Dm-嵌入維數(shù)m的曲線如圖6，發(fā)現(xiàn)關(guān)聯(lián)維數(shù)Dm并非趨于飽和，而是一直在增加。所以G-P法并不能得到安康站降水時(shí)間序列數(shù)據(jù)存在混沌特性的結(jié)果，但也并非說明其不存在隨機(jī)時(shí)間序列的可能。

圖6 關(guān)聯(lián)維數(shù)Dm隨嵌入維數(shù)m的變化情況

通過C-C關(guān)聯(lián)積分法求得延遲時(shí)間τ和嵌入維數(shù)m，對(duì)降水時(shí)間序列重構(gòu)空間，得到最大Lyapunov指數(shù)λi=0.14，說明安康站降水是存在混沌特性的，則對(duì)于1950—2014年安康站降水系統(tǒng)的最大可能預(yù)報(bào)時(shí)間尺度Tf≈7，說明基于現(xiàn)有65 a 780個(gè)月的降水?dāng)?shù)據(jù)，最大可能預(yù)報(bào)的時(shí)間長(zhǎng)度是7個(gè)月。由于降雨與徑流之間具有直接的因果關(guān)系，則基于一定序列長(zhǎng)度的安康站漢江降雨數(shù)據(jù)，可以估計(jì)安康站漢江中期徑流預(yù)報(bào)的最大可能時(shí)間尺度也應(yīng)為7個(gè)月，中長(zhǎng)期徑流預(yù)報(bào)的數(shù)據(jù)時(shí)間尺度的選擇更應(yīng)以月為單位[22]。但基于怎樣尺度的降水-徑流數(shù)據(jù)以及混沌導(dǎo)致的降水與徑流的延遲效應(yīng)值得進(jìn)一步探討。

6 結(jié) 論

(1)基于安康站1950—2014年65 a來降水量趨勢(shì)變化考察其降水特征，可以大致分為4個(gè)上升下降的大周期，分別為1950—1960年、1961—1972年、1973—1999年、2000—2014年。

(2)利用C-C法來確定對(duì)780個(gè)月的降水時(shí)間序列相空間重構(gòu)所需的延遲時(shí)間以及嵌入維數(shù)，并同時(shí)對(duì)比自相關(guān)函數(shù)所求的延遲時(shí)間，2種方法的結(jié)果略有差異。G-P關(guān)聯(lián)維法所得到的結(jié)果表明，關(guān)聯(lián)維數(shù)Dm并沒有隨著嵌入維數(shù)的增加趨于飽和，但也并非否決了降水時(shí)間序列存在混沌的可能性。最大Lyapunov 指數(shù)λi=0.14，說明安康站降水是存在混沌特性的，對(duì)于1950—2014年安康站降水系統(tǒng)的最大可能預(yù)報(bào)時(shí)間尺度Tf≈7，即根據(jù)65 a共780個(gè)月的降水?dāng)?shù)據(jù)最大可能預(yù)測(cè)到接下來7個(gè)月的降水量。

(3)多個(gè)方法所求的結(jié)果不一致，表明了混沌特性存在的可能性需要進(jìn)一步分析。在此基礎(chǔ)上將多種方法的預(yù)測(cè)值進(jìn)行相互比較，不同方法使用時(shí)所考慮的噪聲對(duì)混沌特性識(shí)別的影響可能導(dǎo)致了所求得結(jié)果的差異性。