基于多因子最近鄰抽樣回歸模型的徑流相似性預報

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

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

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

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

1 徑流相似性預報

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

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

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

(1)

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

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

(2)

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

(3)

(4)

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

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

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

(5)

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

1.2 滾動預報方式

為了延長預見期，實時接入滾動降雨預報信息，滾動預報未來徑流過程，提出3種滾動預報方式。

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

(6)

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

(7)

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

(8)

1.3 多模式自適應切換

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

2 應 用 實 例

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

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

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

2.2 未來1 d預報

2.2.1 預報方案設置

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

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

表1 未來1 d的預報方案Table 1 Forecast schemes with foreseen period of 1 day

2.2.2 預報結果分析

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

表2 預見期為1 d的預報精度評價Table 2 Forecast accuracy with foreseen period of 1 day

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

圖2 預見期為1 d的預報結果Fig.2 Forecast results with 1-day forecast period

2.3 未來1周預報

2.3.1 預報方案設置

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

表3 預見期為7 d的滾動預報方案Table 3 Forecast schemes with foreseen period of 7 days

2.3.2 預報結果分析

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

表4 預見期為7 d的預報精度評價Table 4 Forecast accuracy with foreseen period of 7 days

圖3 典型洪水過程的滾動預報結果Fig.3 Rolling forecast results of typical flood

2.4 多模式自適應預報

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

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

圖4 典型徑流過程的自適應預報結果Fig.4 Adaptive forecast results of typical runoff

表5 自適應模式的預報精度評價Table 5 Forecast accuracy of adaptive mode

3 結 語

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