張 勇

（甘肅省臨澤縣平川水利管理所，甘肅臨澤 734200）

蒸散發(fā)量是水文循環(huán)的重要環(huán)節(jié)，準(zhǔn)確的預(yù)算流域蒸散發(fā)量能深入了解流域內(nèi)陸—?dú)庵g水熱交換和水分循環(huán)過程，對于流域水資源規(guī)劃管理和利用極為重要。通過了解流域參考作物蒸散發(fā)量的非線性變化特征及其特征的變化規(guī)律，對預(yù)測黃河流域參考作物蒸散發(fā)量未來幾年可能的非線性變化、指導(dǎo)黃河流域綠色農(nóng)業(yè)參考作物種植區(qū)域生態(tài)結(jié)構(gòu)的變化進(jìn)行科學(xué)良性調(diào)整、合理開發(fā)和利用可調(diào)配的種植地區(qū)水土資源、促進(jìn)種植區(qū)域綠色農(nóng)業(yè)和種植地區(qū)經(jīng)濟(jì)與生態(tài)環(huán)境的良性平衡發(fā)展均有重要意義[1]。

1 研究背景

近年來，在我國的水文科學(xué)中已經(jīng)開始了廣泛應(yīng)用混沌理論。在我國，丁晶最先深入研究混沌理論，他論述了將混沌運(yùn)動特性分析的理論和方法，并把它廣泛應(yīng)用于暴雨和洪水的時(shí)間序列，證明了在洪水過程中的混沌現(xiàn)象[2]。本次研究以黃河流域?yàn)檠芯繉ο螅孟嚓P(guān)氣象數(shù)據(jù)，通過Penman-Monteith公式計(jì)算出參考作物蒸散發(fā)量，進(jìn)而分析得出黃河流域參考作物蒸散發(fā)量的時(shí)間和空間變化特征，最后運(yùn)用量化方法描述出不同時(shí)空黃河流域參考作物蒸散發(fā)量的混沌特性。

2 研究區(qū)域概況

黃河流域總面積752 443 km2，涉及9個(gè)?。▍^(qū)），黃河流域自然和生態(tài)條件復(fù)雜多樣，東西走向的跨越構(gòu)成我國地形的三大階梯，由于區(qū)域內(nèi)地形以及氣候條件的復(fù)雜多樣性，造成了水資源的時(shí)空分布不均勻，洪澇災(zāi)害頻發(fā)，水資源的供需不平衡，人和自然的矛盾日益突出，嚴(yán)重影響了流域內(nèi)經(jīng)濟(jì)的可持續(xù)發(fā)展。

3 研究的方法

3.1 Penman-Monteith公式法

采用FAO1998年推薦的Penman-Monteith公式法計(jì)算逐日參考作物蒸散發(fā)量[3]：

式中，ET0參考作物蒸散發(fā)量為逐日參考作物蒸散發(fā)量（mm/d）；Rn為冠層表面凈輻射[MJ/（m2·d）]；G為土壤熱通量[MJ/（m2·d）]，一般情況下為0；T為平均溫度（℃）；Δ為飽和水汽壓溫度曲線的斜率（KPa/℃）；μ2為高度2.0 m處風(fēng)速（m/s）；es為飽和水汽壓（KPa）；ea為實(shí)際水汽壓（KPa）；γ為濕度常數(shù)（KPa/℃）。

3.2 相空間重構(gòu)理論

一般時(shí)間序列主要是在時(shí)間域中進(jìn)行模型研究，混沌時(shí)間序列必須在相空間中模型研究，因此說相空間重構(gòu)是分析時(shí)間序列的關(guān)鍵環(huán)節(jié)[4]。

自相關(guān)函數(shù)法：x(t)的時(shí)間序列x1，x2，…，xn，自相關(guān)函數(shù)表達(dá)式為：

G-P算法（飽和關(guān)聯(lián)維數(shù)法）：本文選擇飽和關(guān)聯(lián)維數(shù)法來確定相空間的嵌入維數(shù)，其原理如下：對應(yīng)于嵌入維數(shù)M的相空間，關(guān)聯(lián)積分C(r)定義為重構(gòu)相空間點(diǎn)對{Yi，Yj}中，距離||Yi-Yj||小于給定正數(shù)r的數(shù)目在所有相點(diǎn)中占的比例：

其中，M=N-(m-1)τ，為總的相點(diǎn)數(shù)；f為Heaviside函數(shù)，表達(dá)式如下：

r的選取滿足下式：

對式（5）進(jìn)行變換:

在實(shí)際應(yīng)用中，經(jīng)常給定一組M值，選取適當(dāng)?shù)膔，繪出lnC(r)-lnr曲線，求曲線當(dāng)中直線段的斜率，即為關(guān)聯(lián)維數(shù)D。通過觀察關(guān)聯(lián)維數(shù)D隨嵌入維數(shù)M變化的過程圖，找到圖中變化趨于飽和的點(diǎn)，此時(shí)對應(yīng)的M為最佳嵌入維數(shù)，對應(yīng)的D為飽和關(guān)聯(lián)維數(shù)。有無飽和關(guān)聯(lián)維數(shù)的存在是判斷系統(tǒng)是否具有混沌特性的一個(gè)指標(biāo)，也是定量描述混沌特性的重要指標(biāo)。在同一區(qū)域，同一個(gè)時(shí)間段內(nèi)，時(shí)間序列的飽和關(guān)聯(lián)維數(shù)越大，混沌特性越明顯。

3.3 最大Lya punov指數(shù)

采取小數(shù)據(jù)量法計(jì)算最大Lyapunov指數(shù)，首先對時(shí)間序列進(jìn)行快速傅里葉變換（FFT），求出平均周期P。重構(gòu)相空間{Yj，j=1，2，…，M}，在重構(gòu)相空間的基礎(chǔ)上，尋找給定軌道上每個(gè)點(diǎn)的最近鄰近點(diǎn)，即

對于每一個(gè)離散步長i，求出所有j的lndj(i)的平均y(i)，即

在y（i）與i的關(guān)系圖中用最小二乘法擬合，直線的斜率即為最大Lyapunov指數(shù)[5]。如果系統(tǒng)具有混沌特性，則Lyapunov指數(shù)大于零。

4 實(shí)例應(yīng)用及結(jié)果分析

4.1 參考作物蒸散發(fā)量的時(shí)間序列變化特征

利用相關(guān)氣象數(shù)據(jù)通過Penman-Monteith公式計(jì)算出參考作物蒸散發(fā)量，以此作為混沌識別的基礎(chǔ)序列。以瑪曲站為代表站繪制時(shí)序圖，從圖1可以看出，參考作物蒸散發(fā)量的時(shí)間序列的過程圖表現(xiàn)出強(qiáng)烈的季節(jié)周期性和波動性，但是僅從圖像并不能判斷出此序列是否為混沌序列，需要借助于其他方法來判斷參考作物蒸散發(fā)量的時(shí)間序列中是否存在混沌現(xiàn)象[6]。

圖1 瑪曲站參考作物蒸散發(fā)量的時(shí)間序列變化過程

4.2 參考作物蒸散發(fā)量混沌特性識別

由于混沌理論在水文領(lǐng)域的發(fā)展歷程較短，沒有任何一位學(xué)者明確提出某一種混沌特性識別的方法，因此使用單一的方法進(jìn)行時(shí)間序列的識別可靠性不高[7]。本文對序列的識別主要采用的是定量分析中的飽和關(guān)聯(lián)維數(shù)法和最大Lyapunov指數(shù)法進(jìn)行對比識別。

4.2.1 飽和關(guān)聯(lián)維數(shù)法混沌特性識別 由公式（2），計(jì)算出黃河流域瑪曲站點(diǎn)1970—2009年的參考作物蒸散發(fā)量年月均值的自相關(guān)函數(shù)值，并且作出自相關(guān)函數(shù)值隨延遲時(shí)間變化的圖像，選擇自相關(guān)函數(shù)值第一次經(jīng)過零點(diǎn)的時(shí)間作為延遲時(shí)間。如圖2所示。

圖2 瑪曲站參考作物蒸散發(fā)量時(shí)間序列的自相關(guān)函數(shù)圖

從圖3中可以看出，自相關(guān)函數(shù)值第1次經(jīng)過零點(diǎn)的延遲時(shí)間為3，即瑪曲站的延遲時(shí)間τ=3。

圖3 瑪曲站lnC(r)-lnr關(guān)系圖

由自相關(guān)函數(shù)的計(jì)算結(jié)果延遲時(shí)間τ=3，設(shè)維數(shù)M=2，3，4，…，16，在MATLAB軟件中運(yùn)行G-P算法的相關(guān)程序，輸出lnC（r）-lnr的曲線圖，如圖3所示，從上向下的曲線分別是維數(shù)M=2，3，4，…，16，從lnC（r）-lnr關(guān)系圖中確定出直線段的斜率，即為關(guān)聯(lián)維數(shù)D，并作出關(guān)聯(lián)維數(shù)D隨嵌入維數(shù)M的變化的過程圖，如圖4所示。

圖4 瑪曲站m-D關(guān)系圖

由圖4可知，當(dāng)M≥6時(shí)達(dá)到飽和，增長趨勢變?yōu)槠骄?，所以嵌入維數(shù)M=6時(shí)對應(yīng)的飽和關(guān)聯(lián)維數(shù)D2=2.15；這種關(guān)聯(lián)維數(shù)隨著嵌入維數(shù)的增大而達(dá)到飽和的這個(gè)狀態(tài)可以說明參考作物蒸散發(fā)的時(shí)間序列存在確定性動力系統(tǒng)，判斷出此時(shí)間序列為混沌序列，反之，為隨機(jī)序列[8]。

4.2.2 最大Lyapunov指數(shù)法混沌特性識別 由于本文數(shù)據(jù)量很小，所以選擇了易于操作且計(jì)算精度相對較高的小數(shù)據(jù)量法。首先利用快速傅里葉變換估算出時(shí)間序列的平均周期P，取最高譜線所對應(yīng)的下標(biāo)求出平均周期。運(yùn)用MATLAB中的相關(guān)程序計(jì)算最大Lyapunov指數(shù)，計(jì)算結(jié)果見圖5。最后利用最小二乘法擬合得到的直線的斜率就是參考作物蒸散發(fā)量時(shí)間序列的最大Lyapunov指數(shù)λ?，斍镜淖畲驦yapunov指數(shù)λ為0.007。

圖5 瑪曲站時(shí)間序列y（i）和i的關(guān)系圖

4.3 參考作物蒸散發(fā)量混沌特性分析

黃河流域瑪曲站點(diǎn)的混沌特性的分析方法采取以下兩種方法，一是飽和關(guān)聯(lián)維數(shù)法，二是最大Lyapunov指數(shù)法，計(jì)算成果如表1所示。

分析表1可知，黃河流域瑪曲站均存在飽和關(guān)聯(lián)維數(shù)且都為分?jǐn)?shù)維，并且有最大Lyapunov指數(shù)λ＞0，充分說明黃河流域參考作物蒸散發(fā)量時(shí)間序列存在混沌現(xiàn)象。

5 結(jié)論

本文選取了黃河流域內(nèi)的瑪曲氣象站點(diǎn)1970—2009年的逐日氣象資料，采用Penman-Monteith公式計(jì)算出參考作物蒸散發(fā)量的時(shí)間序列。對參考作物蒸散發(fā)量時(shí)間序列的混沌識別，主要在相空間重構(gòu)的基礎(chǔ)上，通過飽和關(guān)聯(lián)維數(shù)法和最大Lyapunov指數(shù)法兩種量化方法，得出了時(shí)間序列的關(guān)聯(lián)維數(shù)隨著嵌入維數(shù)的增加而達(dá)到飽和狀態(tài)且為分?jǐn)?shù)維，并且均有最大Lyapunov指數(shù)λ＞0，計(jì)算結(jié)果均表明了黃河流域參考作物蒸散發(fā)量時(shí)間序列存在混沌特征。這不僅可以了解黃河流域參考作物蒸散發(fā)量的非線性變化規(guī)律，而且為預(yù)測流域內(nèi)參考作物蒸散發(fā)量的未來變化趨勢提供依據(jù)。