譚喬鳳，陳 然，朱 陽，胡立春，聞 昕

(1.河海大學(xué)水利水電學(xué)院，江蘇 南京 210098； 2.國電大渡河流域水電開發(fā)有限公司，四川 成都 610041)

徑流預(yù)報(bào)是防洪減災(zāi)、水資源保障、電力生產(chǎn)等流域管理和調(diào)度決策的關(guān)鍵依據(jù)，也是全球變化下水文水資源領(lǐng)域研究的前沿?zé)狳c(diǎn)[1-2]。目前廣泛運(yùn)用的徑流預(yù)報(bào)方法主要分為過程驅(qū)動和數(shù)據(jù)驅(qū)動兩類[3]。前者從徑流的物理成因出發(fā)，通過一系列含參數(shù)的數(shù)學(xué)物理方程描述產(chǎn)匯流過程，且每個參數(shù)都有明確的物理意義，但是當(dāng)流域自然地理環(huán)境或產(chǎn)匯流條件復(fù)雜時，模型參數(shù)率定困難[4-5]；后者是以建立輸入輸出數(shù)據(jù)最優(yōu)數(shù)學(xué)關(guān)系為目標(biāo)的黑箱子方法[6-7]，無需考慮徑流形成的物理機(jī)制，預(yù)報(bào)精度依賴于可靠且海量的樣本數(shù)據(jù)[8]。

近年來，由于人工智能和大數(shù)據(jù)技術(shù)的快速發(fā)展，流域內(nèi)現(xiàn)代化的水文測站網(wǎng)絡(luò)的建立與完善，以及歷史水文氣象數(shù)據(jù)的積累，為突破徑流預(yù)報(bào)瓶頸開辟了新思路和新途徑。借助大數(shù)據(jù)技術(shù)，可對成千上萬條歷史數(shù)據(jù)實(shí)現(xiàn)全面、多層次的分析，隱藏于數(shù)字背后的水文規(guī)律可以“挖掘”出來。實(shí)際上，對于某一特定流域，在一定的地理環(huán)境條件下，制約降雨的主導(dǎo)天氣系統(tǒng)和降雨類型會反復(fù)出現(xiàn)，在相似天氣系統(tǒng)條件下所產(chǎn)生的降雨(或暴雨)過程及其徑流(或洪水)過程也將是相似的。當(dāng)具有較長時間、較多場次的歷史降雨、徑流過程時，采用數(shù)據(jù)挖掘技術(shù)，找出相似降雨-徑流過程，“參考過去預(yù)測未來”，即根據(jù)歷史相似降雨、徑流提供的信息預(yù)測未來徑流，這對延長徑流預(yù)報(bào)預(yù)見期，提高流域水庫群精細(xì)化管理及精益化調(diào)度水平具有重要意義。

1 徑流相似性預(yù)報(bào)

1.1 多因子最近鄰抽樣回歸模型

研究將基于實(shí)時降雨、歷史降雨和徑流的相似性特征，利用多因子最近鄰抽樣回歸模型(NNBR)在歷史降雨產(chǎn)流過程中搜索相似過程，預(yù)測后期徑流過程。NNBR無需對研究對象的相依形式和概率分布作假定，能避免確定函數(shù)關(guān)系所帶來的不確定性問題，是一類數(shù)據(jù)驅(qū)動的非參數(shù)模型[9]。其基本原理是：客觀世界的發(fā)生發(fā)展存在一定聯(lián)系，未來的運(yùn)動軌跡與歷史具有相似性，即未來發(fā)展模式可以從已知眾多模式中去尋求。對于確定控制斷面的徑流預(yù)報(bào)，首先識別出與其相關(guān)性較高的前滯降雨、徑流因子，建立降雨、徑流特征指標(biāo)矢量：

Xj=(X1,j,X2,j,…,Xp,j，Xp+1,j,Xp+2,j,…,Xp+q,j)

(1)

式中：Xj——預(yù)報(bào)因子；p——降雨特征指標(biāo)個數(shù)；q——徑流特征指標(biāo)個數(shù),共有m個特征指標(biāo)，m=p+q。

考慮到降雨和徑流在數(shù)值量級上存在差異，引入降雨影響權(quán)重α，提出降雨、徑流綜合歐氏距離的計(jì)算公式：

(2)

根據(jù)上述思想，用式(3)表示多因子NNBR的基本形式：

(3)

(4)

模型涉及的參數(shù)有：最近鄰數(shù)、特征矢量維數(shù)、抽樣權(quán)重、降雨影響權(quán)重。

a. 最近鄰數(shù)K：模型根據(jù)已有的特征矢量，尋找歷史中最近鄰特征矢量的個數(shù)，也是模型的抽樣個數(shù)。K的選取主要考慮兩個因素：(a)降雨、徑流的歷史數(shù)據(jù)長度；(b)預(yù)報(bào)的徑流量級在歷史中出現(xiàn)的頻率。首先根據(jù)歷史數(shù)據(jù)長度預(yù)選一個最近鄰數(shù)，當(dāng)模型建立之后，測試不同K值下的預(yù)報(bào)效果，選取預(yù)報(bào)效果最好的值作為最終的K值。

b. 特征矢量維數(shù)：徑流特征矢量維數(shù)由徑流時間序列在不同滯時下的自相關(guān)性確定；降雨特征矢量維數(shù)需要結(jié)合降雨、徑流物理成因，分析降雨與預(yù)報(bào)斷面徑流響應(yīng)時間來確定，這個時間不僅與降雨量級有關(guān)，還與降雨中心的地理位置到預(yù)報(bào)斷面的距離有關(guān)。

(5)

d. 降雨影響權(quán)重：研究結(jié)合遺傳算法等現(xiàn)代啟發(fā)式算法，以預(yù)報(bào)效果最優(yōu)為目標(biāo)，優(yōu)化降雨影響權(quán)重α。其中，預(yù)報(bào)效果評價指標(biāo)主要包括納什系數(shù)(NS)、均方根誤差(RMSE)、平均絕對誤差(MAE)和平均相對誤差(MARE)[10]。

1.2 滾動預(yù)報(bào)方式

為了延長預(yù)見期，實(shí)時接入滾動降雨預(yù)報(bào)信息，滾動預(yù)報(bào)未來徑流過程，提出3種滾動預(yù)報(bào)方式。

a. 方式1：不考慮降雨預(yù)報(bào)和徑流預(yù)報(bào)信息，直接根據(jù)已經(jīng)發(fā)生的實(shí)測降雨和徑流信息尋找降雨、徑流相似過程，以預(yù)報(bào)t～t+f時段的徑流過程為例。

(6)

b. 方式2：考慮降雨預(yù)報(bào)信息滾動更新降雨相似性，同時不斷加入預(yù)報(bào)的徑流作為輸入滾動更新徑流相似性以預(yù)報(bào)下一時段的徑流。以預(yù)報(bào)t～t+f時段的徑流過程為例：

(7)

c. 方式3：考慮降雨預(yù)報(bào)信息滾動更新降雨相似性，但是不將上一時段預(yù)報(bào)的徑流作為輸入，預(yù)測下一時段徑流。

(8)

1.3 多模式自適應(yīng)切換

針對流域在漲退水等不同階段的產(chǎn)匯流特性，建立可根據(jù)實(shí)時水雨情自適應(yīng)切換的降雨、徑流輸入模式，以進(jìn)一步提高徑流預(yù)報(bào)精度。例如，根據(jù)考慮前期降雨滯時長度的不同，分為短降雨滯時模式和長降雨滯時模式。短降雨滯時模式只考慮短期內(nèi)的實(shí)測降雨，例如只考慮前1 d的降雨；長降雨滯時模式考慮更長時期內(nèi)的實(shí)測降雨，例如考慮前3 d之內(nèi)的降雨。短降雨滯時模式以退水期預(yù)報(bào)為主，長降雨滯時模式以漲水期預(yù)報(bào)為主，模型將根據(jù)前期徑流、前期降雨、預(yù)報(bào)降雨等因素自適應(yīng)切換預(yù)報(bào)模式。另外，歷史庫中不同量級徑流出現(xiàn)頻率不同，對于大量級洪水過程，歷史上出現(xiàn)的相似降雨和徑流過程較少，模型將對徑流進(jìn)行實(shí)時追蹤，當(dāng)徑流數(shù)值到達(dá)某一個量級時，自動減小最近鄰數(shù)以保證抽樣樣本具有相似性。

2 應(yīng) 用 實(shí) 例

2.1 研究區(qū)域與數(shù)據(jù)

大渡河流域是國家規(guī)劃的十三大水電基地之一，提高徑流預(yù)報(bào)精度，延長預(yù)見期，對流域水庫群精細(xì)化管理及精益化調(diào)度具有重要意義。本文以大渡河干流上游丹巴斷面為研究對象。斷面以上共有21個雨量遙測站點(diǎn)，3個流量站(綽斯甲、足木足、馬爾康)，站點(diǎn)分布情況見圖1。為降低模型的輸入維度，將預(yù)報(bào)斷面以上的流域降雨站點(diǎn)分為5個子區(qū)域，并建立各分區(qū)點(diǎn)雨量到面雨量的空間映射關(guān)系，在體現(xiàn)流域降雨空間分布格局的前提下降低降雨數(shù)據(jù)輸入維度。區(qū)域1包含色達(dá)、壤塘、一林場3個雨量站；區(qū)域2包含班瑪、燈塔2個雨量站；區(qū)域3包含俄熱、太陽河2個雨量站；區(qū)域4包含日部、草登、足木足、馬爾康、刷經(jīng)寺5個雨量站；區(qū)域5包含丹東、邊爾、布科、東谷、沙沖、牦牛、撫邊河、小金、達(dá)維9個雨量站。研究收集了各站點(diǎn)2008—2019年的逐日歷史實(shí)測降雨、徑流資料，并以2008—2016年的數(shù)據(jù)構(gòu)建歷史降雨、徑流樣本，對2017—2019年的徑流進(jìn)行相似性預(yù)測，以評價模型未來1 d和未來7 d逐日的預(yù)測效果。

圖1 大渡河流域丹巴以上水文站點(diǎn)分布Fig.1 Distribution of hydrological stations above Danba in Dadu River Basin

2.2 未來1 d預(yù)報(bào)

2.2.1 預(yù)報(bào)方案設(shè)置

預(yù)見期為1 d的預(yù)報(bào)因子分為3類：上游雨量(點(diǎn)或面雨量)、上游徑流、丹巴本站徑流。根據(jù)相關(guān)性分析，上游前1 d雨量與丹巴徑流相關(guān)性最高，因此可以初步取降雨滯時為1 d。丹巴徑流自相關(guān)系數(shù)在前5 d分別為：0.96、0.87、0.79、0.69、0.61。前4 d滯時下的相關(guān)性系數(shù)明顯下降，因此丹巴本站前期徑流滯時取3 d。另外，研究測試了預(yù)報(bào)結(jié)果為絕對值和相對值兩種輸出模式。絕對值模式根據(jù)相似性樣本下一時段徑流的絕對值預(yù)報(bào)實(shí)時徑流；相對值模式根據(jù)相似性樣本下一時段徑流相對于前一時段徑流的變化量預(yù)報(bào)實(shí)時徑流。

根據(jù)不同的輸入因子和輸出模式，研究設(shè)置了多種預(yù)報(bào)方案，如表1所示。其中，方案1～4采用絕對值模式，方案1與方案2是為了對比點(diǎn)雨量和面雨量輸入對預(yù)報(bào)效果的影響。方案3與方案2是為了對比上游支流站點(diǎn)徑流對預(yù)報(bào)效果的影響。方案4的預(yù)報(bào)因子只考慮丹巴本站前期徑流，可測試只考慮徑流自身時間序列的預(yù)報(bào)效果。方案5～10采用相對值模式，6種不同滯時組合下預(yù)報(bào)方案形成對比。值得注意的是，方案5與方案2的預(yù)報(bào)因子輸入完全相同，但是方案2采用絕對值模式，而方案5采用相對值模式，因此這兩種方案可直接對比2種輸出模式的預(yù)報(bào)效果。

表1 未來1 d的預(yù)報(bào)方案Table 1 Forecast schemes with foreseen period of 1 day

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

方案1～10的預(yù)報(bào)精度評價結(jié)果見表2。從表2可知，采用面雨量輸入的方案2各項(xiàng)指標(biāo)均優(yōu)于采用點(diǎn)雨量輸入的方案1。采用相對值輸出的方案5～10的預(yù)報(bào)效果整體優(yōu)于采用絕對值輸出的方案1～4，因此預(yù)報(bào)模型輸入取面雨量、輸出取相對值更適合本案例。方案2與方案3的精度對比可知，是否考慮前1 d上游3站的徑流對丹巴預(yù)報(bào)效果影響不大。與方案3相比，方案4的預(yù)報(bào)精度明顯降低，說明前期降雨對丹巴徑流影響顯著。方案5～10共6種預(yù)報(bào)方案的預(yù)報(bào)效果整體接近，在NS等于0.96的方案7～9中，方案7的RMSE雖然比方案8、方案9略高，但是MARE最低，因此確定方案7是預(yù)見期1 d的最優(yōu)方案。

表2 預(yù)見期為1 d的預(yù)報(bào)精度評價Table 2 Forecast accuracy with foreseen period of 1 day

從2017—2019年3 a的預(yù)報(bào)結(jié)果來看，當(dāng)預(yù)見期為1 d時，模型的整體預(yù)報(bào)精度較高，但是在個別洪峰處的預(yù)報(bào)效果較差。以2018汛期(5—10月)為例(圖2)，在年最大洪峰處的預(yù)報(bào)值明顯低于實(shí)際值，可能的原因有兩個：一是NNBR的預(yù)報(bào)結(jié)果是歷史相似樣本的加權(quán)平均值，容易出現(xiàn)極大值預(yù)報(bào)偏小和極小值預(yù)報(bào)偏大的情況；二是歷史相似的大量級洪水場次較少，導(dǎo)致模型抽樣得到的樣本相似性較弱。

圖2 預(yù)見期為1 d的預(yù)報(bào)結(jié)果Fig.2 Forecast results with 1-day forecast period

2.3 未來1周預(yù)報(bào)

2.3.1 預(yù)報(bào)方案設(shè)置

為了延長徑流預(yù)報(bào)預(yù)見期，實(shí)時接入7 d降雨預(yù)報(bào)信息，預(yù)報(bào)流域未來1周的徑流過程。由于缺失歷史降雨預(yù)報(bào)資料，用實(shí)測降雨代替降雨預(yù)報(bào)信息。在預(yù)見期為1 d的預(yù)報(bào)方案基礎(chǔ)上，提出3個滾動預(yù)報(bào)方案。3種方案的預(yù)報(bào)因子均是上游5個區(qū)域前1 d面降雨和丹巴本站前3 d、前2 d、前1 d徑流，并且均采用相對值輸出模式。3種方案分別對應(yīng)章節(jié)1.2所介紹的3種滾動預(yù)報(bào)方式，具體方案如表3所示。

表3 預(yù)見期為7 d的滾動預(yù)報(bào)方案Table 3 Forecast schemes with foreseen period of 7 days

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

滾動預(yù)報(bào)方案1～3的預(yù)報(bào)精度見表4。滾動方案2在7 d內(nèi)預(yù)報(bào)效果最好，滾動方案1與滾動方案3的預(yù)報(bào)效果很接近。且隨著預(yù)見期的增加，預(yù)報(bào)精度顯著下降。圖3顯示了滾動方案2對2019年汛期典型洪水的預(yù)報(bào)過程圖，黑色線為實(shí)測徑流過程，每一個日期圖例表示以該日期為預(yù)報(bào)起點(diǎn)的未來7 d的徑流預(yù)報(bào)過程，預(yù)報(bào)結(jié)果逐日滾動更新。典型徑流過程的預(yù)報(bào)結(jié)果與實(shí)測值較為接近，整體預(yù)報(bào)精度較高。

表4 預(yù)見期為7 d的預(yù)報(bào)精度評價Table 4 Forecast accuracy with foreseen period of 7 days

圖3 典型洪水過程的滾動預(yù)報(bào)結(jié)果Fig.3 Rolling forecast results of typical flood

2.4 多模式自適應(yīng)預(yù)報(bào)

在進(jìn)行未來1周徑流預(yù)測時發(fā)現(xiàn)當(dāng)前1 d的降雨P(guān)t-1延長至前3 d降雨P(guān)t-3,Pt-2,Pt-1時，能夠更好地預(yù)測大渡河丹巴斷面的汛期漲水過程。通過分析Pt-1和Pt-3,Pt-2,Pt-1兩種輸入的滾動預(yù)測效果，研究確定了模型自適應(yīng)判斷條件，當(dāng)滿足判斷條件時，模型將從1 d模式切換至3 d模式；否則，默認(rèn)采用1 d模式進(jìn)行滾動預(yù)報(bào)。同時，研究發(fā)現(xiàn)，當(dāng)前期徑流量級超過3 500 m3/s時，應(yīng)適當(dāng)減少最近鄰數(shù)，推薦取4。

采用多模式自適應(yīng)切換對2017—2019年的滾動預(yù)報(bào)效果見表5?？傮w精度相比于原滾動方案2有顯著提升。預(yù)見期為3 d的NS大于0.9，MARE小于10%；預(yù)見期為7 d的NS大于0.8，差MARE小于15%。與原滾動方案2相比，在徑流發(fā)展的不同階段自適應(yīng)預(yù)報(bào)效果均有提升，其中漲水階段的預(yù)報(bào)效果提升最為明顯(圖4)。自適應(yīng)預(yù)報(bào)帶來的精度提升主要源于最近鄰樣本相似性增強(qiáng)和累計(jì)誤差的減小。通過判斷當(dāng)前降雨、徑流狀態(tài)，選擇不同的預(yù)報(bào)模式尋找相似性樣本，可以有效提高預(yù)報(bào)精度，減小累計(jì)誤差。

圖4 典型徑流過程的自適應(yīng)預(yù)報(bào)結(jié)果Fig.4 Adaptive forecast results of typical runoff

表5 自適應(yīng)模式的預(yù)報(bào)精度評價Table 5 Forecast accuracy of adaptive mode

3 結(jié) 語

基于多因子最近鄰抽樣回歸模型開展徑流相似性預(yù)報(bào)，預(yù)報(bào)效果整體較好。同時，針對流域在漲退水等不同階段的產(chǎn)匯流特性，建立可根據(jù)實(shí)時水雨情自適應(yīng)切換的降雨、徑流輸入模式，進(jìn)一步提高了徑流預(yù)報(bào)精度。在大渡河丹巴斷面的應(yīng)用效果表明預(yù)報(bào)效果達(dá)到預(yù)期，3 d預(yù)見期的NS大于0.9，MARE小于10%；7 d預(yù)見期的NS大于0.8，MARE小于15%。未來可以對洪峰處的徑流過程做進(jìn)一步深入研究，優(yōu)化不同應(yīng)用環(huán)境下的模型參數(shù)，優(yōu)選特定徑流過程的相似性樣本，全面提升徑流預(yù)報(bào)精度。