崔延華，王 理，吳文強

(1.廣西大藤峽水利樞紐開發(fā)有限責(zé)任公司，廣西 南寧 530000；2.水利部珠江水利委員會水文局，廣東 廣州 510611；3.中國水利水電科學(xué)研究院，北京 100038)

人類對水文過程的研究由來已久，隨著人類活動軌跡的擴展和科技水平的提高，不斷增加著新的內(nèi)容[1]。河川徑流量是水資源綜合開發(fā)利用、科學(xué)管理和優(yōu)化調(diào)度最重要的依據(jù)[2]，河川徑流的演變是水文學(xué)的重要研究內(nèi)容[3]。徑流量的科學(xué)預(yù)測對地區(qū)水資源合理分配和社會-經(jīng)濟-生態(tài)可持續(xù)發(fā)展具有重要作用[4]。徑流變化也具有明顯的灰色特征[5]。在中國，將灰色理論應(yīng)用于河流徑流量的預(yù)測開始于20世紀(jì)80年代，至今該理論的應(yīng)用已經(jīng)相當(dāng)成熟[6]。李建林等[7]基于灰色預(yù)測模型對黑河出山年徑流量進行了灰色預(yù)測。陳立華等[8]運用RVA法分析黔江中游水文情勢變化特征。目前,針對大藤峽出山徑流預(yù)報的研究較少，原因是大藤峽水利樞紐工程是國內(nèi)在建大型水利工程，興建時間晚，相關(guān)研究尚未涉足或正在開展，資料收集和數(shù)據(jù)分析需要時間過程。但是大藤峽水利樞紐工程作為黔江流域控制性水利工程，大藤峽出山徑流預(yù)報及預(yù)測精度準(zhǔn)確性都對大藤峽水利樞紐工程水情預(yù)報與水庫調(diào)度具有重要意義?；诖?，為了填補大藤峽出山徑流預(yù)報研究空白，本文利用R/S分形分析法能夠通過改變研究對象的樣本序列周期長度，有效解決樣本序列隨時間推移波動性大問題的優(yōu)點，將灰色預(yù)測模型與R/S分形分析法有效結(jié)合起來，預(yù)測大藤峽出山徑流量V，通過構(gòu)建R/S-灰色預(yù)測模型并對大藤峽出山徑流量V的預(yù)測結(jié)果及模型預(yù)測精度進行對比分析，同時對灰色預(yù)測模型與R/S分形分析法在大藤峽出山徑流演變規(guī)律及未來預(yù)測中的適用性進行分析討論，研究結(jié)果對大藤峽水利樞紐工程水情預(yù)報與水庫調(diào)度具有一定的參考價值。

1 研究方法

1.1 GM(1,1)模型

運用GM(1,1)模型可以對所研究時間序列對象的演變規(guī)律進行系統(tǒng)分析，并對其未來變化趨勢做出中長期預(yù)測?；诖?，本文建立灰色系統(tǒng)GM(1,1)模型[9]。

設(shè)原始時間序列為：

X(0)(k)={x(0)(1),x(0)(2),…,x(0)(n)}

(1)

對X(0)(k)作一次累加計算，得到：

X(1)(k)={x(1)(1),x(1)(2),…,x(1)(n)}

(2)

則GM(1,1)模型的微分方程為：

(3)

運用最小二乘法計算得出GM(1,1)模型參數(shù)a、u，即：

(a,u)T=(BTB)-1BTY

(4)

離散化所研究時間序列對象函數(shù)，進而對式(3)求解得：

x(1)(k+1)=[x(0)(1)-u/a]e-ak+u/a

(5)

經(jīng)過累減還原計算可得出其原始時間數(shù)據(jù)序列的估計值：

(6)

設(shè)大藤峽出山徑流量V實測序列見式(1)，其對應(yīng)的預(yù)測序列為：

(7)

則平均相對殘差值為：

(8)

精度為：

(9)

如果精度P≥0.8，說明模型精度滿足要求，反之要對模型進行優(yōu)化[10]。

1.2 R/S分析

R/S分析法也被稱為分形時間序列[11-12]。本文運用R/S分形分析法計算得出Hurst指數(shù)，進而判定大藤峽出山徑流量V和時間序列的平均循環(huán)周期T。分形時間系列滿足：

R(n)/S(n)=(an)H

(10)

式中R(n)/S(n)——重標(biāo)極差；n——增量區(qū)間長度；a——常數(shù)；H——Hurst指數(shù)[13]。

具體計算步驟如下。

步驟二計算Dm(n個期間)的累積離差Xk,m。

(11)

步驟三由Xk,m算出Dm的極差Rm，以及其對應(yīng)的標(biāo)準(zhǔn)差Sm。

步驟四Dm的重標(biāo)極差通過式(12)計算得出：

(R/S)m=Rm/Sm

(12)

步驟五對每個連續(xù)大藤峽出山徑流量V子時間序列，重復(fù)式(10)—(12)，則該連續(xù)大藤峽出山徑流量V子時間序列的重標(biāo)極差為：

(13)

步驟六要求n≥3，且n為3的倍數(shù)，時間序列{x(k)}符合式(10)要求，進而得到式(14)：

lg(R/S)n=Hlga+Hlgn

(14)

Hurst指數(shù)H的有效值范圍介于0～1[14]。為了定量表征大藤峽出山徑流量V時間序列的平均循環(huán)周期T，引入統(tǒng)計量Q(n)：

(15)

Q(n)隨Hurst指數(shù)H的變化而變化。Q(n)-lgn關(guān)系曲線圖上的拐點n就是大藤峽出山徑流量V時間序列的平均循環(huán)周期T。

1.3 R/S-GM(1,1)模型

長時間序列存在波動性問題，運用R/S-GM(1,1)模型處理波動性問題[15]，步驟如下。

步驟一截取大藤峽出山徑流(T-1)年數(shù)據(jù)作為初始值，標(biāo)記為Q={q1,q2,…,qT-1}。以時間序列Q為基礎(chǔ)建立GM(1,1)模型。如果GM(1,1)模型精度P>80%，則預(yù)測模型的殘差檢驗結(jié)果滿足要求；否則，要求對GM(1,1)模型進行優(yōu)化，使得P>80%；以下預(yù)測步驟都需要進行模型精度的檢驗。

步驟五按照步驟得出大藤峽出山年徑流T年預(yù)測值。

運用式(1)—(6)進行預(yù)測，運用式(7)—(9)進行驗證。

2 研究區(qū)及數(shù)據(jù)來源

2.1 研究區(qū)域概況

大藤峽水利樞紐位于珠江流域西江水系的黔江河段大藤峽谷出口的弩灘處，地理坐標(biāo)為東經(jīng)110°02′，北緯23°28′。壩址在北江橋上游6.6 km處。壩址以上控制流域面積為198 612 km2。大藤峽壩址以上流域近似呈東西向矩形，地勢呈西北高東南低，平均海拔高度約為1 450 m。壩址以上河流自西北向東南方向流動。黔江段河道狹窄彎曲，黔江流經(jīng)武宣盆地，自勒馬始，河流切割大瑤山形成著名的40多公里的大藤峽，峽谷兩岸山高坡陡，坡降為0.5‰～0.6‰，河床縱剖面起伏劇烈，個別地點江底高程低于珠江基面零點以下50余米。黔江出大藤峽谷后進入丘陵平原地帶。

2.2 數(shù)據(jù)來源

本研究的徑流資料采用桂平水文站1932—2021年逐年徑流量，并對上述資料進行了復(fù)核，結(jié)果表明各站資料基本合理，可以采用。大藤峽出山徑流1932—2021年逐年實測徑流量V序列見圖2。從圖2可知，大藤峽出山徑流1932—2021年逐年實測徑流量V序列的總體特征為：大藤峽出山徑流V在長期變化中呈現(xiàn)緩慢下降趨勢，大藤峽出山徑流V的變化幅度隨時間變化較為明顯。

3 結(jié)果分析

運用灰色模型，對大藤峽出山徑流1932—2021年共計90年的逐年實測序列進行預(yù)測分析。

3.1 徑流量的GM(1,1)灰色預(yù)測模型

首先，運用式(1)—(9)取大藤峽出山徑流1932—2015年的逐年實測數(shù)據(jù)構(gòu)建大藤峽出山徑流V的GM(1,1)預(yù)測模型，模型建立結(jié)果見表1。其中大藤峽出山徑流V為x1(0)。

表1 GM(1,1)模型的大藤峽出山徑流V預(yù)測結(jié)果 單位：億m3

累加大藤峽出山徑流量x1(0)得出累積大藤峽出山徑流量x1(1)，把累積大藤峽出山徑流量x1(1)代入式(3)求得a=-0.002 4，u=1 443.303，其中a為發(fā)展系數(shù)，u為灰作用量，進而將時間響應(yīng)函數(shù)離散化，利用初始條件對式(3)求解得：

(16)

經(jīng)累減還原可得大藤峽出山徑流量V數(shù)據(jù)序列的估計模型：

(17)

式中，k=1,2,…,n-1。

根據(jù)式(6)預(yù)測2016—2020年的逐年大藤峽出山徑流量V，其預(yù)測值見表1。由式(8)計算大藤峽出山徑流量GM(1,1)模型的平均相對殘差值，逐年殘差值見圖3，大藤峽出山徑流量GM(1,1)模型平均殘差值為0.156 2，則模型精度是84.38%；由表1可以看出，GM(1,1)模型的預(yù)測精度是86.28%，均大于80%，達到了灰色模型建立標(biāo)準(zhǔn)，預(yù)測效果良好。

3.2 徑流量的R/S模型分析

本文統(tǒng)計了大藤峽出山徑流1932—2021年逐年實測徑流量V序列，并對之進行變化趨勢分析。運用R/S模型方法分析大藤峽出山徑流V的平均循環(huán)周期T。

分別取n為3、5、6、9、10、15、18、30、45，依照式(10)—(14)對大藤峽出山徑流V序列進行分析計算，得到lgn與lg(R/S)的擬合關(guān)系，見圖4。

由圖4可以得出lg(R/S)與lgn之間的關(guān)系式為：lg(R/S)= 0.8049lgn-0.624,其中相關(guān)系數(shù)R2=0.988 9。因此，大藤峽出山徑流V序列的Hurst指數(shù)為0.804 9。H明顯大于0.5，這表明大藤峽出山徑流V序列為分形時間序列，其持續(xù)性和記憶性表現(xiàn)顯著。根據(jù)式(10)—(14)對大藤峽出山徑流V序列進行分析計算，得到lgn與lg(R/S)的擬合關(guān)系(圖4)。

由圖5可知：Q(n)-lgn關(guān)系曲線的趨勢性在箭頭所指處發(fā)生突然轉(zhuǎn)折，由拐點之前的上升趨勢突變?yōu)橄陆第厔?，而拐點處橫坐標(biāo)所對應(yīng)的時間長度n即為大藤峽出山徑流V序列的平均循環(huán)周期T。即大藤峽出山徑流V曲線圖上的拐點n就是大藤峽出山徑流量V時間序列的平均循環(huán)周期T。大藤峽出山徑流V序列在n=9時(箭頭所指處)出現(xiàn)一個突變點，這表明大藤峽出山徑流V序列的平均循環(huán)周期T=9。T反映了初始條件對大藤峽出山徑流V序列產(chǎn)生影響的平均循環(huán)周期長度，換言之，T表征了大藤峽出山徑流對初始條件的平均記憶長度，在9年后初始條件對大藤峽出山徑流V序列的影響作用將完全消失。

3.3 徑流量的R/S-GM(1,1)模型

R/S-GM(1,1)模型就是將R/S分形分析的結(jié)果應(yīng)用到GM(1,1)預(yù)測模型中，即在一個循環(huán)周期T內(nèi)進行大藤峽出山徑流量V的GM(1,1)模型預(yù)測，以減少由于時間序列的波動性引起的誤差。大藤峽出山徑流V的平均循環(huán)周期T為9a，因此，本文選取2008—2015年共計8a大藤峽出山徑流量V作為基礎(chǔ)數(shù)據(jù)，利用1.3節(jié)中的計算步驟，建立大藤峽出山徑流V的R/S-GM(1,1)預(yù)測模型，并對其進行精度驗證，預(yù)測結(jié)果見表2。

表2 R/S-GM(1,1)模型的大藤峽出山徑流V預(yù)測結(jié)果 單位：億m3

由表2可以看出，大藤峽出山徑流量R/S-GM(1,1)模型的模型精度為87.46%，該模型的預(yù)測精度為92.54%，R/S-GM(1,1)預(yù)測模型的精度大于90%，根據(jù)灰色模型精度檢驗等級，精度大于90%即為優(yōu)秀，因此，R/S-GM(1,1)預(yù)測模型不僅通過了殘差檢驗，達到了灰色模型建立標(biāo)準(zhǔn)，而且模型預(yù)測精度優(yōu)秀。

3.4 模型預(yù)測精度對比

由表3知，GM(1,1)模型和R/S-GM(1,1)模型的模型精度分別為84.38%和87.46%，預(yù)測精度分別為86.28%和92.54%，后者的模型精度和預(yù)測精度均高于前者。分析討論如下。

表3 不同模型的大藤峽出山徑流V預(yù)測精度對比 %

a)GM(1,1)模型以1932—2015年的實測大藤峽出山徑流V為基礎(chǔ)進行建模，由于所選取的大藤峽出山徑流原始數(shù)據(jù)序列時間跨度相對比較長，其帶來的波動性也相對比較大，因此影響GM(1,1)模型的模型精度和預(yù)測精度。

b)對比GM(1,1)模型，R/S-GM(1,1)模型第一步運用R/S分形分析得出大藤峽出山徑流序列的平均循環(huán)周期T等于9年，第2步選取連續(xù)9 a的實測大藤峽出山徑流V序列作為原始數(shù)據(jù)，在9 a的大藤峽出山徑流序列基礎(chǔ)上建立R/S-GM(1,1)模型，這樣做的好處是使得大藤峽出山徑流序列的時間跨度大幅減小，進而導(dǎo)致大藤峽出山徑流序列的波動性也隨之大幅降低，此外，R/S-GM(1,1)模型運用遞補法，在建立模型過程中用大藤峽出山徑流量的預(yù)測值代替實測值，目的是減緩大藤峽出山徑流序列的波動性，以上步驟提高了R/S-GM(1,1)模型的預(yù)測精度，使之高于GM(1,1)模型。該R/S-GM(1,1)模型的優(yōu)點是有效合理減小研究對象的序列長度，從而降低序列的波動性。綜上所述，R/S-GM(1,1)模型在GM(1,1)模型基礎(chǔ)上進一步提高了模型精度和預(yù)測精度。R/S-GM(1,1)模型適合用于預(yù)測未來中長期的大藤峽出山徑流量。

4 結(jié)論

大藤峽水利樞紐工程是紅水河綜合利用規(guī)劃中最末一個梯級，也是黔江干流乃至整個黔江流域唯一的水利樞紐，距離上游干流最近的橋鞏電站212 km，天然河道長，受上游電站人為調(diào)控影響較小，自然因素影響較大。因此，本文采用灰色理論模型方法，根據(jù)實測的大藤峽出山徑流量V序列建立灰色預(yù)測模型，結(jié)果表明，大藤峽出山徑流量V時間序列的分形特征表現(xiàn)顯著，對大藤峽出山徑流量V時間序列而言，R/S分形方法是一種有效的分析方法，R/S分形分析法中的2個重要參數(shù)是平均循環(huán)周期T和Hurst 指數(shù)H。Hurst指數(shù)H的作用是判定大藤峽出山徑流量V時間序列的狀態(tài)特性和分形結(jié)構(gòu)，平均循環(huán)周期T的作用是檢驗R/S分形分析法在使用過程中的穩(wěn)定性，判斷大藤峽出山徑流V序列的平均循環(huán)周期T是否存在，并據(jù)此估計平均循環(huán)周期T的長度。大藤峽出山徑流V時間序列的灰色特征表現(xiàn)顯著，在大藤峽出山徑流量的預(yù)測過程中，GM(1,1)模型的模型精度是84.38%，預(yù)測精度是86.28%；而在GM(1,1)模型的基礎(chǔ)上先進行R/S分形分析，再判定出大藤峽出山徑流序列的平均循環(huán)周期是9 a，基于此構(gòu)建的大藤峽出山徑流量R/S-GM(1,1)預(yù)測模型，其模型精度為87.46%，預(yù)測精度為92.54%，R/S-GM(1,1)模型相較于GM(1,1)模型的優(yōu)點是有效合理減小了研究對象的序列長度，從而降低了研究對象序列的波動性。R/S-GM(1,1)預(yù)測模型為大藤峽出山徑流量的中長期預(yù)測提供了一種新的方法。