崔延華，王 理，吳文強(qiáng)

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

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

1 研究方法

1.1 GM(1,1)模型

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

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

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

(1)

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

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

(2)

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

(3)

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

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

(4)

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

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

(5)

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

(6)

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

(7)

則平均相對(duì)殘差值為：

(8)

精度為：

(9)

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

1.2 R/S分析

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

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

(10)

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

具體計(jì)算步驟如下。

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

(11)

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

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

(R/S)m=Rm/Sm

(12)

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

(13)

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

lg(R/S)n=Hlga+Hlgn

(14)

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

(15)

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

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

長(zhǎng)時(shí)間序列存在波動(dòng)性問題，運(yùn)用R/S-GM(1,1)模型處理波動(dòng)性問題[15]，步驟如下。

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

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

運(yùn)用式(1)—(6)進(jìn)行預(yù)測(cè)，運(yùn)用式(7)—(9)進(jìn)行驗(yàn)證。

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

2.1 研究區(qū)域概況

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

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

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

3 結(jié)果分析

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

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

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

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

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

(16)

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

(17)

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

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

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

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

分別取n為3、5、6、9、10、15、18、30、45，依照式(10)—(14)對(duì)大藤峽出山徑流V序列進(jìn)行分析計(jì)算，得到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序列為分形時(shí)間序列，其持續(xù)性和記憶性表現(xiàn)顯著。根據(jù)式(10)—(14)對(duì)大藤峽出山徑流V序列進(jìn)行分析計(jì)算，得到lgn與lg(R/S)的擬合關(guān)系(圖4)。

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

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

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

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

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

3.4 模型預(yù)測(cè)精度對(duì)比

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

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

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

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

4 結(jié)論

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