基于統(tǒng)計(jì)模型的西江枯季中長(zhǎng)期徑流預(yù)報(bào)研究

藍(lán)羽棲，張 尹，農(nóng)振昌，韋永江

(1.水利部珠江水利委員會(huì)水文局，廣東 廣州 510611；2.龍灘水電開(kāi)發(fā)有限公司龍灘水力發(fā)電廠，廣西 河池 547300)

流域水安全保障是維持經(jīng)濟(jì)社會(huì)高質(zhì)量發(fā)展的重要一環(huán)，在新形勢(shì)新背景下，保障流域水安全對(duì)經(jīng)濟(jì)社會(huì)的發(fā)展有著至關(guān)重要的作用[1]。然而，隨著經(jīng)濟(jì)社會(huì)的快速發(fā)展，受咸潮、水污染和過(guò)度開(kāi)發(fā)等問(wèn)題的威脅，枯水期水資源供需矛盾愈發(fā)突出[2]。西江作為珠江流域來(lái)水的主要來(lái)源，其中長(zhǎng)期來(lái)水預(yù)報(bào)對(duì)于珠江下游特別是珠江三角洲枯水期意義非凡。2021年汛期，西江來(lái)水為1946年以來(lái)同期第四枯，主汛期(6—8月)“當(dāng)汛不汛”，來(lái)水均持續(xù)偏少。且流域沒(méi)有出現(xiàn)編號(hào)洪水，西江梧州站最高水位為近3年最低，受降雨偏少影響，2021年汛期，西江流域骨干水庫(kù)來(lái)水較多年同期偏少3～4成。面對(duì)前期不容樂(lè)觀的來(lái)水形勢(shì)，此時(shí)，可靠的中長(zhǎng)期徑流預(yù)報(bào)成為了實(shí)現(xiàn)西江枯水期水庫(kù)群聯(lián)合優(yōu)化調(diào)度、提高水資源利用效率的關(guān)鍵前提，對(duì)抗旱保供水工作具有深刻的科學(xué)價(jià)值。為提高西江水量調(diào)度水平，探求流域內(nèi)主要站點(diǎn)和重要水庫(kù)中長(zhǎng)期徑流預(yù)報(bào)方法具有重要意義。

針對(duì)中長(zhǎng)期徑流預(yù)報(bào)，目前通行的方法是在分析水文要素自身的演變規(guī)律或挖掘與徑流相關(guān)的前期水文氣象資料的基礎(chǔ)上，構(gòu)建徑流時(shí)間序列模型或前期水文氣象要素與預(yù)報(bào)月徑流的映射關(guān)系，從而提供科學(xué)的、預(yù)見(jiàn)期較長(zhǎng)的徑流預(yù)測(cè)結(jié)果。近年來(lái)，針對(duì)統(tǒng)計(jì)模型在水文領(lǐng)域中的應(yīng)用已有了諸多研究，隨著計(jì)算機(jī)的應(yīng)用與發(fā)展，機(jī)器學(xué)習(xí)也逐步被引入到中長(zhǎng)期預(yù)報(bào)中。目前常用的統(tǒng)計(jì)學(xué)方法有隨機(jī)森林法[3-4]、神經(jīng)網(wǎng)絡(luò)算法[5]、季節(jié)自回歸等[6]，如胡義明等[7]基于AdaBoost模型、隨機(jī)森林模型和支持向量機(jī)模型在淮河流域進(jìn)行了月徑流預(yù)報(bào)，比較分析了各模型的預(yù)報(bào)效果和特點(diǎn)，發(fā)現(xiàn)AdaBoost模型整體上優(yōu)于其他2個(gè)機(jī)器學(xué)習(xí)模型；酈于杰等[8]將支持向量回歸機(jī)應(yīng)用于漢江皇莊站的中長(zhǎng)期徑流預(yù)報(bào)，并對(duì)預(yù)報(bào)結(jié)果進(jìn)行了不確定性分析，實(shí)現(xiàn)了高精度的定值預(yù)報(bào)并以置信區(qū)間的方式量化了預(yù)報(bào)的可靠度；常新雨等[9]采用灰色關(guān)聯(lián)分析法篩選預(yù)報(bào)因子，構(gòu)建了基于數(shù)據(jù)驅(qū)動(dòng)的深度神經(jīng)網(wǎng)絡(luò)、Elman神經(jīng)網(wǎng)絡(luò)和支持向量機(jī)3種模型并運(yùn)用于黃龍灘水庫(kù)旬月徑流預(yù)報(bào)，發(fā)現(xiàn)支持向量機(jī)在洪峰預(yù)報(bào)上誤差更小。上述研究大多采用的是機(jī)器學(xué)習(xí)方法，也均取得了較好的預(yù)報(bào)效果。然而，盡管機(jī)器學(xué)習(xí)方法學(xué)習(xí)能力強(qiáng)、容錯(cuò)性能高，但對(duì)數(shù)據(jù)容量具有較高要求，相比之下統(tǒng)計(jì)相關(guān)分析類(lèi)和時(shí)間序列類(lèi)方法原理簡(jiǎn)單、應(yīng)用方便、耗時(shí)短，因此后者仍是實(shí)際中長(zhǎng)期徑流預(yù)報(bào)應(yīng)用的常用方法。鑒于此，本文基于均生函數(shù)、周期分析、多元逐步回歸3種統(tǒng)計(jì)模型對(duì)西江流域枯季月徑流和整個(gè)枯水期徑流開(kāi)展應(yīng)用研究。

1 研究區(qū)域及資料

西江是珠江流域的主干流，發(fā)源于云南省曲靖市烏蒙山余脈的馬雄山東麓，自西向東流經(jīng)云南、貴州、廣西、廣東4省(自治區(qū))，至廣東佛山市三水區(qū)思賢滘，全長(zhǎng)2 075 km，平均坡降0.58‰，流域面積353 120 km2，占珠江流域總面積的77.8%。干流從上而下由南盤(pán)江、紅水河、黔江、潯江及西江5個(gè)河段組成。本文以龍灘水庫(kù)和梧州水文站為代表進(jìn)行枯季徑流中長(zhǎng)期預(yù)報(bào)方法的應(yīng)用研究。

龍灘水庫(kù)位于西江流域干流紅水河河段(圖1)，是廣西境內(nèi)最大的水電工程，也是西江防洪工程體系的控制性工程。龍灘水庫(kù)以發(fā)電為主，兼具防洪、航運(yùn)等綜合效益。壩址以上流域面積98 500 km2，占紅水河流域面積的71.2%。本文選取龍灘水庫(kù)1958—2020年天然徑流枯水期6個(gè)月(10月、11月、12月、次年1月、次年2月、次年3月)的數(shù)據(jù)，分別對(duì)各月、整個(gè)枯水期建立統(tǒng)計(jì)模型進(jìn)行徑流模擬，其中1958—2004年采用龍灘實(shí)測(cè)入庫(kù)流量，2005—2020年采用經(jīng)上游天生橋一級(jí)水庫(kù)、光照水電站還原后的龍灘天然入庫(kù)流量。

梧州水文站位于西江干流與支流桂江匯合口以下約3 km處(圖1)，東經(jīng)111°20′、北緯23°28′，是西江流域的出口控制站，集水面積327 046 km2。本文選取梧州站1963—2020年天然徑流枯水期6個(gè)月的數(shù)據(jù)，其中1963—2005年采用梧州實(shí)測(cè)流量，2006—2020年采用經(jīng)上游天生橋一級(jí)水庫(kù)、光照水電站、龍灘水庫(kù)、百色水庫(kù)4座水庫(kù)還原后的梧州天然流量。

圖1 研究區(qū)域

收集了國(guó)家氣候中心1956—2020年的130項(xiàng)氣候監(jiān)測(cè)指數(shù)資料(由于需要利用前期氣候指數(shù)資料，故起始年份向前推2年)、龍灘水庫(kù)和梧州以上流域的前期降雨流量資料(4—9月)作為多元逐步回歸的初選因子集。

2 研究方法

2.1 均生函數(shù)

均生函數(shù)預(yù)報(bào)方法的基本思想是假定事物過(guò)去存在的趨勢(shì)會(huì)延伸到未來(lái)，在分析時(shí)間序列變量的基礎(chǔ)上構(gòu)建預(yù)報(bào)模型，通過(guò)向外延伸時(shí)間趨勢(shì)進(jìn)而確定未來(lái)預(yù)報(bào)結(jié)果[10-11]。

設(shè)任意長(zhǎng)度的時(shí)間序列X(t) (t=1，2，…，N)，構(gòu)造其均生函數(shù)如下：

(1)

由此得到L(L=lmax=[N/2])階均生矩陣，通過(guò)對(duì)均生函數(shù)作周期外延進(jìn)一步構(gòu)造出外延矩陣，最后從中篩選出B個(gè)相關(guān)性強(qiáng)的序列，并基于此構(gòu)建多元回歸模型進(jìn)行模擬預(yù)測(cè)。

2.2 周期分析

周期分析法的基本思想是將變量隨時(shí)間的變化當(dāng)做是由多個(gè)不同周期的周期波疊加而成的[12-13]。在給定置信度條件下，采用F檢驗(yàn)的方法識(shí)別周期波，并對(duì)所識(shí)別的各周期波進(jìn)行外延，線(xiàn)性疊加后可用于預(yù)測(cè)。設(shè)水文要素序列x(t)，周期分析法基本計(jì)算式如下：

(2)

式中Pi(t)——第i個(gè)周期波序列；l——識(shí)別到的周期波個(gè)數(shù)；ε(t)——誤差項(xiàng)。

2.3 多元逐步回歸

采用氣象因子作為預(yù)報(bào)因子進(jìn)行多元回歸計(jì)算的方式較為常見(jiàn)[14-15]。通過(guò)對(duì)在國(guó)家氣候中心網(wǎng)站下載的88項(xiàng)大氣環(huán)流指數(shù)、26項(xiàng)海溫指數(shù)和16項(xiàng)其他氣象指數(shù)，共計(jì)130項(xiàng)天氣學(xué)因子與對(duì)應(yīng)的水文要素進(jìn)行分析綜合，同時(shí)考慮前期流域降雨、江河來(lái)水等因素與預(yù)報(bào)對(duì)象具有較大相關(guān)性，本文將其與上述130項(xiàng)天氣學(xué)因子共同作為預(yù)報(bào)因子，用逐步回歸分析方法與預(yù)報(bào)對(duì)象建立方程，并從中篩選出物理意義明確、統(tǒng)計(jì)貢獻(xiàn)顯著的預(yù)報(bào)因子，對(duì)預(yù)測(cè)對(duì)象進(jìn)行定量預(yù)測(cè)。多元逐步回歸基本計(jì)算式如下：

y=b0+b1·x1+b2·x2+b3·x3+…+bn·xn

(3)

式中y——預(yù)報(bào)對(duì)象；bn——回歸系數(shù)；xn——預(yù)報(bào)因子；n——因子個(gè)數(shù)。

2.4 精度評(píng)價(jià)指標(biāo)

2.4.1變幅誤差(M)

采用GB/T 22482—2008《水文情報(bào)預(yù)報(bào)規(guī)范》[16]中針對(duì)中長(zhǎng)期定量預(yù)報(bào)的精度評(píng)定方法，即以多年同期實(shí)測(cè)變幅的20%作為許可誤差，當(dāng)預(yù)報(bào)值與實(shí)測(cè)值間的絕對(duì)誤差小于許可誤差時(shí)判定為合格(M<20%)，否則不合格。M計(jì)算方法如下：

(4)

2.4.2平均絕對(duì)百分比誤差(MAPE)

該指標(biāo)反映了預(yù)報(bào)值相對(duì)于實(shí)測(cè)值的偏離程度，MAPE越接近于0，則模擬效果越好。MAPE計(jì)算方法如下：

(5)

3 預(yù)報(bào)結(jié)果分析

以龍灘水庫(kù)和梧州水文站為例，采用上述均生函數(shù)、周期分析和多元逐步回歸3種統(tǒng)計(jì)方法分別構(gòu)建徑流預(yù)報(bào)模型，預(yù)報(bào)時(shí)段為當(dāng)年10月至次年3月各月(6個(gè)時(shí)段)和整個(gè)枯水期(1個(gè)時(shí)段)，即2個(gè)斷面各有3×7個(gè)預(yù)報(bào)模型?？紤]到率定期與驗(yàn)證期的長(zhǎng)度比一般符合3∶1左右的原則，對(duì)于龍灘水庫(kù)，以1958—2004年作為模型率定期，2005—2020年作為模型驗(yàn)證期；對(duì)于梧州站，以1963—2005年作為率定期，2006—2020年作為驗(yàn)證期。采用變幅誤差和平均絕對(duì)百分比誤差2個(gè)指標(biāo)對(duì)各模型預(yù)報(bào)結(jié)果進(jìn)行精度評(píng)定，龍灘水庫(kù)、梧州站各月和整個(gè)枯水期的徑流預(yù)報(bào)精度情況見(jiàn)表1、2。

表1 龍灘入庫(kù)各月、枯水期徑流預(yù)報(bào)精度統(tǒng)計(jì) %

表2 梧州站各月、枯水期徑流預(yù)報(bào)精度統(tǒng)計(jì) %

將2個(gè)斷面統(tǒng)一考慮，由表1、2可知，在率定期內(nèi)，3種模型10月至次年3月各月和整個(gè)枯水期(10月至次年3月)的平均合格率均超過(guò)了75%，其中均生函數(shù)模型合格率在91.5%～100%，平均合格率為97.4%，多元逐步回歸模型在80.9%～100%，平均合格率為92.8%，周期分析較其他2個(gè)模型略低，2個(gè)斷面14個(gè)模型的合格率在55.0%～97.5%，平均合格率為78.8%；3種模型對(duì)應(yīng)平均MAPE值均在30%以?xún)?nèi)，其中均生函數(shù)、多元逐步回歸2個(gè)模型平均MAPE值分別為15.3%和19.2%，均在20%以?xún)?nèi)，而周期分析模型相對(duì)較高，平均MAPE值為29.1%。綜合來(lái)看，3個(gè)模型在率定期均具有良好的預(yù)報(bào)效果，其中均生函數(shù)模型和多元逐步回歸模型明顯優(yōu)于周期分析模型。

相對(duì)于率定期而言，驗(yàn)證期的預(yù)報(bào)精度略有降低，其中均生函數(shù)模型驗(yàn)證期平均合格率為89.6%，周期分析模型為74.8%，分別較率定期下降了7.8%和4.0%，而多元逐步回歸模型不降反升，平均合格率為94.5%，較率定期上升了1.7%；3種模型對(duì)應(yīng)平均MAPE值較率定期總體存在不同程度的增漲，相似地，其中均生函數(shù)模型和周期分析模型均增漲10%，而多元回歸模型僅增漲5%，說(shuō)明多元回歸模型具有較強(qiáng)的泛化能力，模型穩(wěn)定性更強(qiáng)，這可能是因?yàn)橄啾惹皟煞N模型僅分析自身演變規(guī)律，多元回歸模型使用前期水文氣象因子作為依據(jù)，獲得了更多的預(yù)報(bào)信息。精確到各月時(shí)發(fā)現(xiàn)，龍灘水庫(kù)次年二月的MAPE值均較高，但合格率卻不低，這是因?yàn)?010年2月龍灘天然來(lái)水是1958年以來(lái)最枯，較多年同期偏少92%，所以該年MAPE值異常偏高，導(dǎo)致平均MAPE值也更高。

為了更直觀地查看和對(duì)比3種模型的預(yù)報(bào)效果，圖2、3分別給出了整個(gè)時(shí)間序列內(nèi)上游龍灘水庫(kù)和下游梧州站預(yù)報(bào)值和實(shí)測(cè)值的散點(diǎn)關(guān)系。

a)均生函數(shù)預(yù)報(bào)

a)均生函數(shù)預(yù)報(bào)

由圖2、3可以明顯看到，在率定期與驗(yàn)證期內(nèi)，均生函數(shù)模型和多元逐步回歸模型的散點(diǎn)均聚攏分布在45°線(xiàn)(圖中紅線(xiàn))附近，且線(xiàn)兩側(cè)散點(diǎn)個(gè)數(shù)比重相當(dāng)，說(shuō)明2個(gè)模型對(duì)枯水期各月流量均具有較好的預(yù)報(bào)效果。相比之下周期分析模型45°線(xiàn)兩側(cè)雖也分布均勻，但散點(diǎn)較離散，且隨著量級(jí)增加，預(yù)報(bào)的偏差越大，上半部分散點(diǎn)大多分布于紅線(xiàn)以下，即在流量量級(jí)較大時(shí)周期分析模型存在預(yù)報(bào)偏低現(xiàn)象。由此可知，均生函數(shù)模型和多元回歸模型的預(yù)報(bào)精度整體較周期分析模型更高，在高流量值的預(yù)報(bào)上誤差更小。另外，通過(guò)對(duì)比龍灘入庫(kù)和梧州站可以看到，前者各模型散點(diǎn)分布更聚集，究其原因，可能與梧州站以上眾多水庫(kù)的調(diào)度調(diào)節(jié)有關(guān)，隨著水庫(kù)的增建和人類(lèi)活動(dòng)的影響增加，經(jīng)4座水庫(kù)還原后的梧州站還原流量并不是真正的天然流量，序列資料一致性遭到一定程度的破壞，導(dǎo)致模型精度降低。

4 結(jié)論

可靠的中長(zhǎng)期徑流預(yù)報(bào)是流域開(kāi)展科學(xué)的枯水期水量調(diào)度工作的重要依據(jù)?；诰瘮?shù)、周期分析、多元逐步回歸3種統(tǒng)計(jì)模型，本文對(duì)西江上游龍灘水庫(kù)、下游梧州水文站進(jìn)行了自10月至次年3月各月以及整個(gè)枯水期的徑流預(yù)報(bào)，得出的結(jié)論如下。

a)3種模型均呈現(xiàn)出較好的預(yù)報(bào)效果，率定期和驗(yàn)證期平均合格率均在75%以上，MAPE值基本在30%以?xún)?nèi)，其預(yù)報(bào)結(jié)果可為西江流域水量調(diào)度方案的編制提供技術(shù)參考。其中均生函數(shù)和多元逐步回歸的預(yù)報(bào)精度相近，且明顯高于周期分析，特別是在對(duì)極值的預(yù)報(bào)能力方面。

b)與均生函數(shù)和周期分析相比，多元逐步回歸模型穩(wěn)定性更強(qiáng)。前兩者驗(yàn)證期平均合格率較率定期均有所下降，后者不降反升，且仍保持在90%以上，MAPE值的增漲也最緩慢。

c)在預(yù)報(bào)天然流量時(shí)，序列資料的一致性程度會(huì)影響預(yù)報(bào)精度。就模型預(yù)報(bào)與實(shí)測(cè)的整體擬合情況而言，上游龍灘水庫(kù)較下游梧州站擬合得更好，這可能與梧州水文站作為西江控制站受到眾多水庫(kù)共同調(diào)控的影響有關(guān)。

d)由于計(jì)算思路不同，統(tǒng)計(jì)模型各有優(yōu)劣，目前還不存在一個(gè)具有絕對(duì)優(yōu)勢(shì)的通用模型，因此可以考慮通過(guò)融合多種模型的方式，達(dá)到優(yōu)化預(yù)報(bào)效果的目的。