基于n階多項(xiàng)式改進(jìn)退水曲線的枯水期徑流預(yù)測方法

董亞迪，廖勝利，燕志宇，張世欽

(1.大連理工大學(xué)水電與水信息研究所，遼寧大連116024；2.國網(wǎng)福建省電力有限公司，福建福州350003)

0 引 言

枯水期流域內(nèi)降水及來水量都較少，水資源供需矛盾突出，無論是跨流域水量分配，還是水電調(diào)峰調(diào)度都要求水庫特別是調(diào)節(jié)性能較差(季調(diào)節(jié)及以下)的水庫有比較準(zhǔn)確的日來水量預(yù)測；因此，開展枯水期徑流預(yù)測對充分發(fā)揮水庫效益及提高水電站發(fā)電調(diào)度水平具有重要意義[1-2]?？菟趶搅髦饕裳茨粼诹饔騼?nèi)的蓄水量消退以及枯水期降雨形成，影響因素非常復(fù)雜，主要包括以氣候和地質(zhì)為主的自然因素及工農(nóng)業(yè)用水等人為因素[3]，河流斷面徑流量呈逐漸消退趨勢但預(yù)測極為困難。

目前，枯水期徑流預(yù)測大多采用比較簡單的傳統(tǒng)方法，主要包括退水曲線法[4-5]、前后期徑流相關(guān)法[6]、回歸分析法[7]、人工神經(jīng)網(wǎng)絡(luò)法[8]和遺傳算法[9]等，尤以退水曲線法應(yīng)用最為普遍。趙銅鐵鋼等[10-11]采用隨機(jī)森林模型對長江上游屏山站和寸灘站枯水期徑流預(yù)報進(jìn)行了探討，研究結(jié)果顯示，隨機(jī)森林模型能有效選取水文—?dú)庀箢A(yù)報因子進(jìn)行徑流預(yù)報。近年來，隨著氣候條件的改變以及人為因素的影響，流域特性和枯水期徑流規(guī)律有了一定的變化[12]，因此為了提高預(yù)測精度，需要對傳統(tǒng)退水曲線法做出改進(jìn)。

本文針對退水曲線公式中消退系數(shù)隨時間逐漸遞增的關(guān)系，提出用n階多項(xiàng)式函數(shù)曲線擬合消退系數(shù)，構(gòu)建了改進(jìn)的退水曲線預(yù)測模型，通過比較各模型預(yù)測值平均偏差率選取最佳預(yù)測模型，并結(jié)合典型年法給出了枯水期徑流預(yù)測的具體步驟。最后將該模型應(yīng)用于瀾滄江流域功果橋水電站枯水期徑流預(yù)測，預(yù)測效果良好，模型可供其他電站或流域參考和應(yīng)用。

1 常規(guī)退水曲線

當(dāng)忽略枯水期少量的降雨時，枯水期徑流的來源主要是汛末滯留于流域內(nèi)的蓄水量，包括地下蓄水量和河網(wǎng)蓄水量。對于由地下水補(bǔ)給的河流，可以近似地認(rèn)為地下蓄水量與出流量之間為線性關(guān)系，因此其蓄量方程可以表示為

-Qg(t)=dWg(t)/dt

(1)

Wg(t)=kgQg(t)

(2)

整理后可得

dGg(t)/Gg(t)=-(1/kg)dt

(3)

求解式(3)可得

Qg(t)=Qg(0)e-t/kg

(4)

式中，Wg(t)為t時刻的地下蓄水量，m3；Qg(t)為t時刻的出流量，m3/s；Qg(0)為退水開始時刻t=0時的出流量，m3/s；ks為常數(shù)。

同理，由河網(wǎng)蓄水量補(bǔ)給的枯水期徑流，其蓄泄關(guān)系也近似呈線性，則出流量Qr(t)的消退規(guī)律是

Qr(t)=Qr(0)e-t/kr

(5)

式中，Qr(t)為t時刻的出流量，m3/s；Qr(0)為退水開始時刻t=0時的出流量，m3/s；kr為常數(shù)。

一般情況下，流域的退水過程是地下蓄水量和河網(wǎng)蓄水量共同補(bǔ)給的結(jié)果，可用一個總的退水公式表示。即

Q(t)=Q(0)e-t/K

(6)

式中，Q(0)、Q(t)分別為退水開始時刻t=0時的流量和t時刻的流量，m3/s；K為常數(shù)，其定義為

W(t)=KQ(t)

(7)

式中，W(t)為t時刻的蓄水量，m3。

式(7)表明，當(dāng)泄流流量恒定為Q(t)時，蓄水量W(t)泄完所需的時間正好為K。由于蓄水量分布在流域的各個位置上，距出口斷面的距離遠(yuǎn)近不同，匯集時間大小不等，因此K可解釋為流域水流平均匯集時間。

若簡化式(6)為

(8)

則

Cg=e-1/K

(9)

式中，Cg為常數(shù)。由于K為流域水流平均匯集時間，則Cg也可反映退水速率的快慢，稱為消退系數(shù)。式(8)即為常規(guī)的退水曲線公式。因此，掌握了Cg的變化規(guī)律就掌握了退水曲線的規(guī)律，分析退水曲線就是分析Cg的變化。

2 改進(jìn)退水曲線預(yù)測模型

2.1 常規(guī)方法存在的問題

傳統(tǒng)方法計算時，首先把退水曲線式(8)表達(dá)為遞推形式(Δt=1 d)

Q(t+1)=CgQ(t)

(10)

由式(10)知，消退系數(shù)Cg可直接由計算時段始末的兩個實(shí)測退水流量來確定。即

Cg=Q(t+1)/Q(t)

(11)

因此，將相鄰時段流量數(shù)據(jù)組成樣本系列，通過均值法或最小二乘法求出該流域的消退系數(shù)Cg。當(dāng)用最小二乘法來估計消退系數(shù)時，表達(dá)式為

(12)

式中，n為樣本系列中相鄰時段流量數(shù)據(jù)組的組數(shù)，即一組樣本系列中有n組觀測值；括號中i為觀測值序號，第二個數(shù)字為前后時段序號。即，Q(i，1)為第i組觀測值的前一時段流量，Q(i，2)為第i組觀測值的后一時段流量。

雖然最小二乘法考慮了一次退水過程每個時段的流量數(shù)據(jù)，但是當(dāng)把求得的Cg作為固定的常數(shù)代入式(8)進(jìn)行枯水期徑流預(yù)測時，會發(fā)現(xiàn)預(yù)測值偏差較大，結(jié)果往往達(dá)不到較好的效果。其原因是，此方法將Cg簡化成了一個恒定常數(shù)；然而，當(dāng)我們直接利用原流量系列Q(t)和Q(0)之間的函數(shù)關(guān)系反向推求Cg，不難發(fā)現(xiàn)Cg并非為一個恒定常數(shù)，而是隨時間t大致呈遞增趨勢。而且在最小二乘法中，Cg的推求是通過Q(t)和Q(t+1)的關(guān)系得到的；而在退水預(yù)測時，卻是運(yùn)用初始時刻的退水量即Q(0)和Q(t)的關(guān)系，這樣難免會造成Cg求解方法上的矛盾及較大的誤差。

2.2 基于n階多項(xiàng)式的改進(jìn)退水曲線

由于Cg并非是一個恒定常數(shù)，其隨時間t大致呈遞增的趨勢。根據(jù)其遞增的規(guī)律和線型，本文采用與傳統(tǒng)方法不同的假設(shè)，認(rèn)為Cg是隨時間變化的n階多項(xiàng)式函數(shù)。即用多項(xiàng)式函數(shù)擬合Cg。根據(jù)枯水期徑流資料，利用Q(t)和Q(0)之間的函數(shù)關(guān)系推求Cg系列值后，用n階多項(xiàng)式函數(shù)曲線進(jìn)行擬合，發(fā)現(xiàn)：當(dāng)階數(shù)大于等于4時，擬合曲線會出現(xiàn)振蕩，引起過擬合情況；隨著多項(xiàng)式階數(shù)增加，曲線擬合偏差會越來越大，徑流預(yù)測值曲線和實(shí)際流量曲線的偏差也會越來越大。因此，本文用多項(xiàng)式函數(shù)擬合Cg時，只考慮一階、二階和三階三種情況。即

(13)

(14)

(15)

將擬合結(jié)果代入式(8)可得改進(jìn)退水曲線公式(下式中的i式(12)已用，不能再表示另一個物理概念，需換符號j，以下相關(guān)均需改)

(16)

式中，j為階數(shù)，取1，2，3；Qi(t)為j階多項(xiàng)式對應(yīng)的t時刻預(yù)測流量值。

故本文采用改進(jìn)退水曲線式(16)進(jìn)行枯水期徑流預(yù)測。關(guān)于改進(jìn)退水曲線的最佳預(yù)測模型選擇標(biāo)準(zhǔn)，本文采用預(yù)測值偏差率式(17)來判斷。首先根據(jù)式(16)進(jìn)行枯水期徑流預(yù)測，然后根據(jù)式(17)計算每個退水曲線公式預(yù)測值對應(yīng)的平均偏差率，取平均偏差率最小值對應(yīng)的退水曲線公式作為此次預(yù)測的最佳預(yù)測模型。即

(17)

2.3 典型年選擇方法

對于某個特定流域，枯水期徑流消退規(guī)律在歷年中往往表現(xiàn)出較好的重復(fù)性?；趎階多項(xiàng)式的改進(jìn)退水曲線模型和常規(guī)退水曲線法一樣，仍然是一種利用歷年徑流變化規(guī)律來預(yù)測枯水期徑流的預(yù)測方法。因此，本文采用典型年法來獲取歷史徑流資料，作為枯水期徑流預(yù)測依據(jù)。

隨著人們生活方式的改變，人類活動的影響和氣候變化對流域枯水期徑流的影響不容忽視[11]。因此，在選擇典型年徑流資料時，這些影響因素需要考慮進(jìn)來。由于近20年內(nèi)流域特性和枯水期退水規(guī)律與預(yù)測年份相比往往更加相似；因此為了提高預(yù)測準(zhǔn)確度，本文將典型年選擇范圍縮短至近20年歷史徑流資料內(nèi)。

典型年流量資料選擇標(biāo)準(zhǔn)有兩個：第一，確定預(yù)測初始時刻退水流量，在歷史同期徑流系列中找出流量大小與預(yù)測年份較為相近的年份作為典型年候選年份(百分誤差不大于3%)。第二，計算預(yù)測初始時刻前10天平均流量，以及典型年候選年份中同期時段的平均流量，從而確定二者的百分誤差。然后將典型年候選年份中初始時刻流量百分誤差和前10天平均流量百分誤差相加，找出百分誤差和最小年份作為典型年。確定第二個選擇標(biāo)準(zhǔn)的原因是，流域預(yù)測初始時刻前期蓄水情況和徑流量對流域退水過程會有較大影響作用。

2.4 預(yù)測求解步驟

利用改進(jìn)退水曲線模型進(jìn)行枯水期徑流預(yù)測時，依據(jù)有兩個：預(yù)測年份初始時刻退水量Q(0)和利用典型年徑流資料擬合得到的Cg多項(xiàng)式函數(shù)。該公式并未考慮流域?qū)嶋H面臨的各種復(fù)雜氣候條件和人為因素。如果預(yù)測時段過長，突發(fā)性氣候和人類活動對退水過程的影響就會大大增加；因此預(yù)測誤差也會更大。長期的枯水期徑流預(yù)測準(zhǔn)確性和可信程度都偏低，使得水庫實(shí)際調(diào)度時必須考慮來水情況的各種變化。因此，為了提高預(yù)測準(zhǔn)確性，本文將一次預(yù)測周期確定為30 d。

本文利用改進(jìn)退水曲線模型進(jìn)行枯水期徑流預(yù)測的流程具體步驟如下：

(1)確定預(yù)測年份初始時間，根據(jù)典型年流量資料選擇標(biāo)準(zhǔn)確定典型年。

(2)根據(jù)典型年流量資料，確定消退系數(shù)Cg實(shí)際變化過程線，然后擬合得到Cg的一階、二階和三階多項(xiàng)式函數(shù)。

(3)將擬合到的Cg多項(xiàng)式表達(dá)式代入常規(guī)退水曲線式(8)，得到改進(jìn)退水曲線公式。

(4)將預(yù)測年份起始流量Q(0)代入改進(jìn)退水曲線式(16)，進(jìn)行逐日流量預(yù)測。

(5)計算每個改進(jìn)退水曲線公式預(yù)測值對應(yīng)的平均偏差率，取平均偏差率最小值對應(yīng)的退水曲線公式作為此次預(yù)測的最佳模型。最佳模型預(yù)測值即為最終預(yù)測成果。

3 實(shí)例研究

3.1 實(shí)例背景及資料

本文以瀾滄江流域?yàn)閷?shí)例進(jìn)行研究。瀾滄江在我國境內(nèi)長2 130 km，落差約5 000 m，流域面積17.4萬km2，水能資源蘊(yùn)藏量約3 656萬kW。全流域?qū)傥髂霞撅L(fēng)氣候，其顯著特點(diǎn)是干、濕兩季分明，一般5月至10月為濕季，11月至次年4月為枯水期，約85%以上降水量集中在濕季；而又以6月～8月為最集中，3個月降水量占全年降水量60%以上。枯水期降水量很小，枯水期徑流主要由流域前期蓄水補(bǔ)給，流量過程呈現(xiàn)較穩(wěn)定的消退趨勢。

功果橋水電站是瀾滄江干流水電基地中下游河段“兩庫八級”梯級開發(fā)方案的龍頭電站。本文以功果橋水電站壩址斷面1996年～2015年中的20年日平均流量資料作為歷史徑流數(shù)據(jù)庫選擇典型年，進(jìn)行枯水期徑流預(yù)測，并用2016年實(shí)際流量數(shù)據(jù)進(jìn)行成果檢驗(yàn)。

3.2 參數(shù)計算及選擇

以2016年11月11日為預(yù)測起始日期，起始流量為998.13 m3/s，前10天的平均流量為1 081.37 m3/s。根據(jù)典型年選擇標(biāo)準(zhǔn)在歷史數(shù)據(jù)中確定典型年為2013年。起始日流量為1 024.69 m3/s，誤差為2.5%，符合候選典型年選擇標(biāo)準(zhǔn)；前10天的平均流量為1 120.41 m3/s，誤差為3.4%，二者的誤差和為5.9%，在候選典型年系列中為最小值，符合典型年選擇標(biāo)準(zhǔn)。

根據(jù)典型年流量數(shù)據(jù)，對消退系數(shù)Cg進(jìn)行多項(xiàng)式函數(shù)擬合，可得Cg的一階、二階和三階多項(xiàng)式

(18)

(19)

0.004t+0.926 2

(20)

將擬合到的Cg表達(dá)式代入常規(guī)退水曲線式(8)，即可得改進(jìn)退水曲線公式。

3.3 計算結(jié)果及分析

圖1為典型年流量系列Cg曲線、多項(xiàng)式擬合曲線和最小二乘法Cg曲線。從圖1可以看出，典型年流量消退系數(shù)Cg隨時間t大致呈遞增的趨勢。當(dāng)采用傳統(tǒng)的最小二乘法處理時，將典型年相鄰時段流量數(shù)據(jù)組成樣本系列，根據(jù)最小二乘法式(12)可求出該流域消退系數(shù)Cg=0.971；此時，退水曲線公式中的Cg被視為一個恒定的常數(shù)，處理后的結(jié)果是實(shí)際的Cg和預(yù)測公式中的Cg有較大偏差。二階和三階多項(xiàng)式函數(shù)擬合效果要好于一階多項(xiàng)式函數(shù)。相比較三階多項(xiàng)式函數(shù)、二階多項(xiàng)式函數(shù)的擬合效果稍有優(yōu)勢。

圖1 典型年Cg曲線和擬合曲線

將2016年11月11日流量Q(0)=998.13 m3/s代入改進(jìn)退水曲線公式對11月12日至12月10日30天流量進(jìn)行預(yù)測；同時，將Q(0)=998.13 m3/s和通過最小二乘法求出的Cg=0.971代入常規(guī)退水曲線式(8)進(jìn)行預(yù)測，得到預(yù)測流量值(見圖2)。根據(jù)式(17)得不同預(yù)測方法下的平均偏差率：最小二乘法為10.85%；一階多項(xiàng)式為7.25%；二階多項(xiàng)式為2.84%；三階多項(xiàng)式為7.37%。

圖2 2016年實(shí)際流量曲線和預(yù)測流量曲線

從圖2的4條預(yù)測流量曲線可明顯看出，最小二乘法預(yù)測結(jié)果與實(shí)際流量的偏差最大，平均偏差率為10.85%；在前期大部分時間多項(xiàng)式函數(shù)預(yù)測結(jié)果和實(shí)際流量都非常接近，但尾部階段差異越來越大。與2016年實(shí)際流量數(shù)據(jù)行對比，二階多項(xiàng)式預(yù)測方法明顯優(yōu)于一階多項(xiàng)式和三階多項(xiàng)式預(yù)測方法。因此，本次預(yù)測的最佳模型為基于二階多項(xiàng)式的改進(jìn)退水曲線模型。

其中，11月29日和30日二階多項(xiàng)式預(yù)測流量與實(shí)際流量相比有些許的偏差。查看歷史降雨資料得知，原因是這段時間內(nèi)瀾滄江流域內(nèi)有一些局部降雨，從而引起實(shí)際流量略大于預(yù)測流量。整體來看，從11月中旬到12月中旬二階多項(xiàng)式的預(yù)測流量曲線與實(shí)際流量曲線基本吻合，平均偏差率僅為2.84%，預(yù)測效果良好。

4 結(jié) 語

本文分析了枯水期徑流預(yù)報的重要性，介紹了枯水期徑流的主要來源和影響因素。針對退水曲線公式中消退系數(shù)隨時間t逐漸遞增的關(guān)系，提出用n階多項(xiàng)式函數(shù)曲線擬合消退系數(shù)，從而構(gòu)建了基于n階多項(xiàng)式的改進(jìn)退水曲線預(yù)測模型，通過比較各模型預(yù)測值平均偏差率選取最佳預(yù)測模型，并結(jié)合典型年法給出了枯水期徑流測的具體步驟。

本文以瀾滄江流域功果橋水電站為實(shí)例，對2016年11月11日至12月10日進(jìn)行了30天的枯水期逐日徑流預(yù)測，通過比較各模型預(yù)測值平均偏差率，確定最佳模型為基于二階多項(xiàng)式的改進(jìn)退水曲線預(yù)測模型，其預(yù)測效果良好，可供其他電站或流域參考和應(yīng)用。