陳江峰，龔書(shū)浩，李迎超

(1.河南理工大學(xué) 資源環(huán)境學(xué)院，河南 焦作 454003；2.河南省核工業(yè)放射性核素檢測(cè)中心, 河南 鄭州 450044； 3.黃河交通學(xué)院，河南 武陟 454950)

降雨量是水循環(huán)過(guò)程中最重要、最活躍的因素之一，降雨量的分布會(huì)直接影響到徑流等水文環(huán)境，同時(shí)降雨量也是描述一個(gè)地區(qū)氣候變化的關(guān)鍵指標(biāo)，一個(gè)地區(qū)降雨量的多少對(duì)當(dāng)?shù)貧夂颦h(huán)境、水資源分布、生態(tài)系統(tǒng)等都會(huì)產(chǎn)生重要影響。由于降雨等自然現(xiàn)象具有突發(fā)性，在發(fā)生時(shí)間上存在不規(guī)則性，降雨量的大小又受多種復(fù)雜因素的影響，因此很難對(duì)降雨量進(jìn)行精準(zhǔn)預(yù)測(cè)。近年來(lái)，國(guó)內(nèi)外研究人員發(fā)現(xiàn)，看似無(wú)序、不規(guī)則的降雨現(xiàn)象實(shí)際并不是完全隨機(jī)的，而是蘊(yùn)含著規(guī)律，比如自相似性，這是傳統(tǒng)理論難以解釋的，因此越來(lái)越多的學(xué)者開(kāi)始利用分形幾何學(xué)等理論來(lái)探討降雨現(xiàn)象的發(fā)生規(guī)律。

Hurst提出的重標(biāo)度極差法(R/S分析法)是分析時(shí)間序列的有力工具，是分形理論的重要方法之一[1]，其優(yōu)勢(shì)在于無(wú)論時(shí)間序列是正態(tài)分布還是非正態(tài)分布，分析結(jié)果的穩(wěn)定性均不受影響，因此迄今仍然是估計(jì)Hurst指數(shù)時(shí)最常用的方法[2]。R/S分析法以其能夠描述非線性時(shí)間序列復(fù)雜變化內(nèi)在規(guī)律的特點(diǎn)，在水文學(xué)和氣象學(xué)研究中有著廣泛應(yīng)用[3-4]。目前利用R/S分析法研究時(shí)間序列取得了很大進(jìn)展，但這類方法的研究與應(yīng)用對(duì)象均是時(shí)間序列的趨勢(shì)性分析，屬于定性預(yù)測(cè)。本研究在R/S趨勢(shì)性分析研究的基礎(chǔ)上，考慮時(shí)間序列長(zhǎng)度和相關(guān)性等因素，輔以計(jì)算機(jī)技術(shù)，實(shí)現(xiàn)了降雨量的定量預(yù)測(cè)。

1 理論基礎(chǔ)

1.1 β-Z-H三參數(shù)綜合趨勢(shì)分析

1.1.1 Mann-Kendall檢驗(yàn)

Mann-Kendall(M-K)檢驗(yàn)法是一種非參數(shù)統(tǒng)計(jì)檢驗(yàn)法，廣泛適用于非正態(tài)分布特征變量的趨勢(shì)分析。設(shè)x1,x2,…,xn為時(shí)間序列數(shù)據(jù)，n為時(shí)間序列數(shù)據(jù)的長(zhǎng)度，對(duì)于所有t、τ≤n，且t≠τ，xt和xτ的分布是不同的，M-K法定義統(tǒng)計(jì)量S為

(1)

其中

式中：xt、xτ分別為序號(hào)t、τ對(duì)應(yīng)的時(shí)間序列數(shù)據(jù)。

S的方差Var(S)計(jì)算公式為

Var(S)=[n(n-1)(2n+5)-

(2)

式中：i為時(shí)間序列數(shù)據(jù)出現(xiàn)的次數(shù)；ei為時(shí)間序列數(shù)據(jù)出現(xiàn)次數(shù)的個(gè)數(shù)。

(3)

在M-K檢驗(yàn)中，另一個(gè)非常有用的指數(shù)為肯德?tīng)栃甭师?，它表示時(shí)間序列趨勢(shì)大小，若β>0,則為上升趨勢(shì),若β<0,則為下降趨勢(shì)，若β=0，則說(shuō)明趨勢(shì)不明顯[6]。β值計(jì)算公式為

(4)

1.1.2 Hurst指數(shù)

R/S分析的基本思想是改變樣本序列的時(shí)間尺度，研究其在不同尺度范圍內(nèi)的統(tǒng)計(jì)規(guī)律，從而進(jìn)行大小時(shí)間尺度間的相互轉(zhuǎn)換[7-8]。該方法能從分形時(shí)間序列中區(qū)分出隨機(jī)和非隨機(jī)序列。通過(guò)Hurst指數(shù)的判定，可以判斷出時(shí)間序列的分形結(jié)構(gòu)和狀態(tài)持續(xù)性，為時(shí)序的復(fù)雜性演變提供一種有效的非線性科學(xué)預(yù)測(cè)方法[9]。

對(duì)于一個(gè)隨機(jī)過(guò)程樣本序列，如果其物理量的時(shí)間序列為x1，x2，…，xn，滿足

(5)

式中：τ為時(shí)間序列的序數(shù)；R(τ)為時(shí)間序列對(duì)應(yīng)的累積偏差的域；S(τ)為R(τ)對(duì)應(yīng)的標(biāo)準(zhǔn)偏差；R(τ)/S(τ)為重標(biāo)極差；H為Hurst指數(shù)，其值介于0到1之間[10],可以在雙對(duì)數(shù)坐標(biāo)系中根據(jù)[lnτ,ln(R/S)]的值用最小二乘法擬合求得。

若H=0.5，則表明現(xiàn)在和未來(lái)之間沒(méi)有相關(guān)性；若0

1.2 初始序列長(zhǎng)度

對(duì)一個(gè)具有狀態(tài)持續(xù)性的時(shí)間序列來(lái)說(shuō)，具有長(zhǎng)記憶效應(yīng)特征。從理論上講，今天所發(fā)生的一切將一直影響未來(lái)，即存在對(duì)初始條件的敏感性[12]，對(duì)時(shí)間序列做Vτ～lnτ分析，即

(6)

此時(shí)Vτ圖像對(duì)應(yīng)τ的循環(huán)變化長(zhǎng)度是時(shí)間序列持續(xù)性的平均長(zhǎng)度，是統(tǒng)計(jì)意義上的數(shù)值，而非時(shí)間序列的真正周期，也是預(yù)測(cè)初始序列長(zhǎng)度選取的可行值[13]。當(dāng)圖中有多個(gè)拐點(diǎn)時(shí)，遵循最大H值且擬合度較高的原則選擇突變點(diǎn)作為循環(huán)長(zhǎng)度的平均值[13]。

1.3 R/S定量預(yù)測(cè)方法的基本原理

H值反映了時(shí)間序列的趨勢(shì)性，若0

(1)對(duì)時(shí)間序列進(jìn)行R/S趨勢(shì)性分析，得到線性回歸方程ln(R/S)=a0lnτ+a1，將ln(τ+1)代入線性回歸方程，得到回歸方程下一點(diǎn)的最優(yōu)值ybest。

(2)將假設(shè)數(shù)據(jù)x(τ+1)作為下一點(diǎn)的新數(shù)據(jù)代入相關(guān)計(jì)算機(jī)程序中進(jìn)行運(yùn)算，得到新的預(yù)測(cè)點(diǎn)值ln(R/S)τ+1，判斷|ln(R/S)τ+1-ybest|<Δt是否成立(Δt可以根據(jù)需求賦值，以誤差最小、預(yù)測(cè)精度最高為準(zhǔn)則)。如成立，則此時(shí)的x(τ+1)即為定量預(yù)測(cè)點(diǎn)，程序停止運(yùn)行，輸出x(τ+1)，否則進(jìn)行下一步。

(3)若上步不滿足，則程序再次改變迭代數(shù)據(jù)，通過(guò)連續(xù)運(yùn)算輸出滿足要求的原始數(shù)據(jù)，即為誤差范圍內(nèi)距離ybest最近的原始數(shù)據(jù)，即為下一點(diǎn)的定量預(yù)測(cè)值x(τ+1)。

2 R/S定量預(yù)測(cè)實(shí)例

2.1 數(shù)據(jù)來(lái)源

阿壩藏族羌族自治州地處青藏高原東南緣、橫斷山脈北端與川西北高山峽谷的接合部，氣溫分布自東南向西北隨海拔升高而降低：海拔2 500 m以下的河谷地帶降水集中、蒸發(fā)快，為干旱、半干旱地帶；海拔2 500～4 100 m的坡谷地帶為寒溫帶，年平均氣溫1～5 ℃；海拔4 100 m以上為寒帶，終年積雪，長(zhǎng)冬無(wú)夏。本研究采用阿壩州小金水文觀測(cè)站實(shí)測(cè)年降雨量數(shù)據(jù)(1961—2000年)進(jìn)行分析和預(yù)測(cè)(圖1)。

圖1 阿壩州年降雨量變化曲線

2.2 年降雨量β-Z-H三參數(shù)綜合趨勢(shì)分析

對(duì)年降雨量數(shù)據(jù)進(jìn)行R/S趨勢(shì)性分析，得到H=0.663 1、R2=0.945 2，表明該時(shí)間序列具有正持續(xù)性，時(shí)間序列在未來(lái)的增減性與原來(lái)是相同的，但持續(xù)性較弱，如圖2所示。

將年降雨量數(shù)據(jù)帶入M-K檢驗(yàn)公式(1)～(4)中得到Z=0.338，β=1.532。β是評(píng)價(jià)時(shí)間序列歷史發(fā)展趨勢(shì)的指標(biāo)，Z可以反映其顯著程度，H能反映時(shí)間序列未來(lái)發(fā)展的持續(xù)性。β、Z、H都是對(duì)時(shí)間序列進(jìn)行分析的指標(biāo)，進(jìn)一步將三者綜合，則不同的數(shù)值對(duì)應(yīng)關(guān)系

圖2 年降雨量R/S趨勢(shì)預(yù)測(cè)

具有不同的指示意義。在α=0.05置信水平上，β-Z-H三參數(shù)綜合分析對(duì)降雨演化特征的指示作用如表1所示。

表1 β、Z、H綜合分析的指示作用

根據(jù)之前得出的Z=0.338，β=1.532，H=0.663 1，結(jié)合β、Z、H綜合分析指示作用表分析知，研究區(qū)年降雨量時(shí)間序列符合序號(hào)7的特征，即降雨量過(guò)去是增加的，未來(lái)將持續(xù)增加，但是趨勢(shì)不明顯，這與R/S趨勢(shì)分析得到的結(jié)果是相同的。

2.3 初始序列長(zhǎng)度

如果初始序列長(zhǎng)度過(guò)大，數(shù)據(jù)的波動(dòng)性將會(huì)增強(qiáng)，回歸方程的線性關(guān)系就會(huì)減弱，因此選取合適的初始序列長(zhǎng)度對(duì)預(yù)測(cè)結(jié)果的準(zhǔn)確性尤為重要。對(duì)上述時(shí)間序列做Vτ～lnτ分析(圖3)，時(shí)間序列循環(huán)變化長(zhǎng)度為12，因此選取12為時(shí)間序列預(yù)測(cè)的初始長(zhǎng)度。

圖3 降雨量時(shí)間序列長(zhǎng)記憶性分析

2.4 R/S定量預(yù)測(cè)

記1982—1993年降雨量數(shù)據(jù)為序列A，以后新的序列依次命名為序列B、C等。對(duì)該時(shí)間序列進(jìn)行R/S分析，得到回歸方程如圖4所示，H=0.548 17，R2=0.977，那么R/S趨勢(shì)預(yù)測(cè)點(diǎn)之間存在良好的線性關(guān)系，可以嘗試對(duì)該時(shí)間序列進(jìn)行定量預(yù)測(cè)，得到的預(yù)測(cè)值x(13) 為628 mm。

圖4 序列A的 R/S定量預(yù)測(cè)

將x(13)放入序列A中得到序列B，對(duì)序列B進(jìn)行同樣運(yùn)算，得到預(yù)測(cè)值x(14)=637 mm(圖5)。同理，得到序列C(圖6)，預(yù)測(cè)值x(15)=636 mm。

圖5 序列B的R/S定量預(yù)測(cè)

圖6 序列C的R/S定量預(yù)測(cè)

2.5 預(yù)測(cè)精度檢驗(yàn)

對(duì)得到的預(yù)測(cè)值與真實(shí)降雨量進(jìn)行精度檢驗(yàn)，結(jié)果見(jiàn)表2。

表2 初始序列長(zhǎng)度為12的定量預(yù)測(cè)值與真實(shí)值對(duì)比

將預(yù)測(cè)值與真實(shí)值對(duì)比可以看出，初次預(yù)測(cè)和第二次預(yù)測(cè)結(jié)果較為理想，特別是第一次預(yù)測(cè)精度為97.6%，但是第三次較前兩次誤差增大，說(shuō)明隨著預(yù)測(cè)時(shí)間序列長(zhǎng)度的增加，時(shí)間序列的波動(dòng)性增強(qiáng)，內(nèi)在規(guī)律性減弱?？傮w而言，模型預(yù)測(cè)的平均相對(duì)殘差為7.1%，平均預(yù)測(cè)精度為92.9%，預(yù)測(cè)效果顯著。

3 結(jié) 論

(1)輔以計(jì)算機(jī)技術(shù)，實(shí)現(xiàn)了R/S分析方法的定量預(yù)測(cè)，將非線性時(shí)間序列問(wèn)題進(jìn)行了線性求解，將該模型應(yīng)用于降雨量預(yù)測(cè)，平均精度高達(dá)92.9%。

(2)R/S定量預(yù)測(cè)與初始序列長(zhǎng)度有關(guān)，連續(xù)多次預(yù)測(cè)會(huì)使模型精度有所降低，因此對(duì)時(shí)間序列長(zhǎng)度的選擇進(jìn)行分析，顯得尤為重要。

(3)隨著人類活動(dòng)的日益增強(qiáng)，降雨量大小預(yù)測(cè)變得更為復(fù)雜，運(yùn)用R/S分析對(duì)降雨量進(jìn)行定量預(yù)測(cè)，可為非線性時(shí)間序列問(wèn)題的解決提供一種快速而有效的方法。