三峽水庫中長期徑流預(yù)測及不確定性分析研究

黃華平，酈于杰，王 棟，靳高陽

（1.中水珠江規(guī)劃勘測設(shè)計(jì)有限公司，廣州 510610；2.浙江省水利水電勘測設(shè)計(jì)院，杭州 310002；3.長江水利委員會水文局，武漢 430010）

0 引 言

中長期徑流預(yù)報(bào)是指基于水文現(xiàn)象的變化規(guī)律，以已知信息作為輸入條件，采用成因分析及數(shù)學(xué)建模等手段對預(yù)見期超過3 天的徑流過程進(jìn)行定量或定性的預(yù)測［1］。作為一種非工程措施，中長期徑流預(yù)報(bào)在水旱災(zāi)害防治、水庫調(diào)度運(yùn)用、水資源優(yōu)化管理及生態(tài)環(huán)境保護(hù)等方面均發(fā)揮了顯著作用［2］。

目前，關(guān)于中長期徑流預(yù)報(bào)已有較多研究，依據(jù)原理的不同，可劃分為以下兩類—基于水文循環(huán)過程驅(qū)動和基于相關(guān)因子數(shù)據(jù)驅(qū)動。前者一般是指將氣象要素預(yù)測數(shù)據(jù)作為輸入項(xiàng)，驅(qū)動水文模型來對徑流過程進(jìn)行預(yù)測。如劉甜等［3］以CFS降水?dāng)?shù)據(jù)為輸入項(xiàng)，驅(qū)動SWAT 模型對丹江口水庫入庫流量進(jìn)行了相關(guān)預(yù)測研究；Singh 等［4］采用ESP 方法獲取未來氣象數(shù)據(jù)，將其導(dǎo)入Topnet 模型，預(yù)測了新西蘭南島四個典型流域的徑流過程。這一類方法具有較為堅(jiān)實(shí)的理論基礎(chǔ)，但預(yù)測過程中需要搜集大量數(shù)據(jù)來驅(qū)動模型運(yùn)轉(zhuǎn)，且降水預(yù)測的不確定性在水文模擬過程中將進(jìn)一步放大，導(dǎo)致預(yù)測結(jié)果往往差強(qiáng)人意，故而該類方法相關(guān)應(yīng)用較少。而基于相關(guān)因子數(shù)據(jù)驅(qū)動是指在分析前期徑流或氣候要素與未來徑流過程相關(guān)性基礎(chǔ)上，建立統(tǒng)計(jì)模型來實(shí)現(xiàn)徑流中長期預(yù)測。該類研究在早期一般采用灰度理論、小波分析、方差分析等數(shù)理統(tǒng)計(jì)方法［5-7］。隨著人工智能與數(shù)據(jù)挖掘技術(shù)的快速發(fā)展，大量機(jī)器學(xué)習(xí)算法被應(yīng)用于中長期徑流預(yù)測研究中，包括神經(jīng)網(wǎng)絡(luò)模型、隨機(jī)森林算法、支持向量機(jī)及深度信念網(wǎng)絡(luò)等［8-11］。但上述研究多數(shù)僅提供了徑流的確定性預(yù)測結(jié)果，缺乏預(yù)測不確定性方面的考慮。

本文將LGB算法與HUP模型相耦合，應(yīng)用于三峽水庫逐月流量中長期預(yù)測中。研究首先采用相關(guān)系數(shù)檢驗(yàn)及逐步回歸法篩選出合理的預(yù)測因子；其次，以預(yù)測因子為輸入項(xiàng)，驅(qū)動LGB 算法對逐月徑流過程進(jìn)行預(yù)測；最終，采用HUP 模型對預(yù)測結(jié)果的不確定性進(jìn)行相應(yīng)分析。

1 計(jì)算方法

1.1 輕量梯度提升樹

輕量梯度提升樹（Light Gradient Boosting Machine，LGB）是微軟于2017年提出的一種以梯度提升決策樹（Gradient Boosting Decision Tree，GBDT）算法為基礎(chǔ)的輕量化串行增強(qiáng)算法［12］。這類算法的基本思想是將訓(xùn)練過程階梯化，每一輪訓(xùn)練以擬合前一輪計(jì)算結(jié)果的殘差值為目標(biāo)，直到殘差值足夠小或迭代次數(shù)達(dá)到最大值為止。具體模型結(jié)構(gòu)如式（1）所示，每輪僅訓(xùn)練一棵決策樹來擬合前期決策樹之和的殘差，然后將新的決策樹加入前期決策樹中構(gòu)成新模型。

與GBDT 算法相比，LGB 算法優(yōu)勢主要反映在以下4 個方面：

（1）采用直方圖算法對輸入數(shù)據(jù)進(jìn)行離散化，減小噪聲影響的同時(shí)降低運(yùn)行所需內(nèi)存，模型訓(xùn)練速度顯著提升；

（2）采用梯度單側(cè)采樣法，保留大梯度樣本和部分小梯度樣本，減小訓(xùn)練數(shù)據(jù)且不改變樣本分布，提高了模型計(jì)算效率；

（3）通過將互斥特征停留在不同直方圖中構(gòu)造特征束，對特征束進(jìn)行排序，從而達(dá)到合并容忍范圍內(nèi)互斥特征的目的，大幅降低遍歷的特征數(shù)量；

（4）采用了一種更為高效的按葉生長策略，在分裂前會遍歷所有葉節(jié)點(diǎn)，選擇具有最大增益的節(jié)點(diǎn)分裂，可降低分裂誤差，得到更高的模擬精度。

1.2 水文不確定性處理器

Krzysztofwicz R［13］于1999年提出了貝葉斯預(yù)報(bào)系統(tǒng)（Bayesian Forecasting System，BFS），用于解決水文預(yù)報(bào)中的不確定性問題。該系統(tǒng)包含降水不確定性處理器（PUP）和水文不確定性處理器（HUP），分別用于處理降水輸入產(chǎn)生的不確定性和除降水外其他所有不確定性。其中，HUP模塊基本原理如下所述：

記H0為已知實(shí)測流量過程，Hn與Sn分別為不同預(yù)見期n對應(yīng)的實(shí)際流量和預(yù)報(bào)流量。依據(jù)貝葉斯原理及實(shí)測流量過程，可推斷出預(yù)報(bào)值Sn=sn條件下，實(shí)際流量Hn=hn的后驗(yàn)密度函數(shù)，如式（2）。

式中：h0為已知實(shí)測流量過程；g(hn|h0)為hn的先驗(yàn)分布，與h0相關(guān)；fn(sn|hn，h0)為sn關(guān)于hn的似然函數(shù)，反映了水文模型的預(yù)報(bào)能力。

假設(shè)每個目標(biāo)點(diǎn)Tl(l=1,2,3,…,m×n)的面積均為Δm×Δn,如果目標(biāo)點(diǎn)被覆蓋,則目標(biāo)點(diǎn)的聯(lián)合感知概率為1,覆蓋面積為Δm×Δn,否則為0,所以目標(biāo)點(diǎn)Tl的覆蓋面積可表示為Il×Δm×Δn;同樣,整個區(qū)域T的總面積為AS=(m×n)(Δm×Δn)。節(jié)點(diǎn)部署后節(jié)點(diǎn)所覆蓋的面積占部署區(qū)域總面積的比值稱為節(jié)點(diǎn)覆蓋率ψ[16],計(jì)算如下:

依據(jù)實(shí)測流量的后驗(yàn)密度函數(shù)φn(hn|sn，h0)，結(jié)合觀測數(shù)據(jù)及任意時(shí)刻的預(yù)報(bào)流量，可以估計(jì)相應(yīng)時(shí)刻實(shí)際值的條件分布函數(shù)，從而對預(yù)報(bào)結(jié)果的不確定性分析進(jìn)行相應(yīng)分析。

1.3 精度評價(jià)指標(biāo)

3種預(yù)測精度指標(biāo)（納什效率系數(shù)、相關(guān)系數(shù)及平均絕對誤差）及三種不確定性區(qū)間指標(biāo)（覆蓋率、相對帶寬及平均偏移程度）被用于評價(jià)模型確定性預(yù)測結(jié)果及不確定性置信區(qū)間，計(jì)算公式下如下：

（1）納什效率系數(shù)（NSE）。

（2）相關(guān)系數(shù)（CC）。

（3）平均絕對誤差（MAPE）。

（4）覆蓋率（CR）。

（5）相對帶寬（RB）。

（6）平均偏移程度（RD）。

2 實(shí) 例

本文搜集了三峽水庫1965-2016年逐月流量數(shù)據(jù)，考慮到不同月份徑流特性與影響成因存在差異，研究對12個月份分別建立LGB 預(yù)測模型，并將GBDT 算法作為對比參照來驗(yàn)證前者精度，在確定性預(yù)報(bào)基礎(chǔ)上，采用HUP 模型對預(yù)測結(jié)果不確定性進(jìn)行了相關(guān)分析。

2.1 預(yù)報(bào)因子篩選

將國家氣象局提供的130 項(xiàng)氣候系統(tǒng)指數(shù)（88 項(xiàng)大氣環(huán)流指數(shù)、26項(xiàng)海溫指數(shù)和16項(xiàng)其他指數(shù)（http：//cmdp.ncc-cma.net/Monitoring/cn_index_130.php）與前期徑流過程作為預(yù)測因子，考慮其與當(dāng)月徑流過程間的遙相關(guān)性，將徑流發(fā)生前12個月氣候系統(tǒng)指數(shù)與前一個月徑流量納入預(yù)測因子初選范疇，采用相關(guān)系數(shù)顯著性檢驗(yàn)及逐步回歸法對上述因子進(jìn)行篩選，具體步驟如下：

（1）搜集各月徑流過程對應(yīng)的所有相關(guān)變量，計(jì)算各相關(guān)變量與當(dāng)月徑流過程間的相關(guān)系數(shù)；

（2）依據(jù)相關(guān)系數(shù)顯著性檢驗(yàn)表，設(shè)置置信度為0.05，篩選出相關(guān)系數(shù)顯著的變量作為備選因子；

（3）采用逐步回歸法對備選因子進(jìn)行分析，依據(jù)不同因子的方差貢獻(xiàn)率，篩選出10 個獨(dú)立性強(qiáng)，相關(guān)系數(shù)高的備選因子作為最終預(yù)測因子。

其中，1月及8月選定的預(yù)測因子如表1所示。

表1 典型月份對應(yīng)預(yù)測因子統(tǒng)計(jì)表Tab.1 Selected predictors for two typical months（January and August）

2.2 模型率定與驗(yàn)證

以1965-2001年為率定期，2002-2016年為驗(yàn)證期，將逐月預(yù)測因子導(dǎo)入模型中，對三峽水庫逐月徑流過程進(jìn)行模擬預(yù)測。模擬過程中，為防止出現(xiàn)過擬合問題，研究采用留一交叉驗(yàn)證法對模型參數(shù)進(jìn)行充分率定，并將實(shí)測系列與模擬系列間MAPE值最小作為率定目標(biāo)，預(yù)測結(jié)果如圖1、2及表2所示。

圖1 三峽水庫1965-2016年逐月流量過程對比Fig.1 The observed and simulated monthly runoff series from 1965 to 2016 for the Three Gorges Reservoir

表2 模型率定期及驗(yàn)證期精度評價(jià)Tab.2 Performance Indices for the calibration and validation periods

由表2 可知，率定期內(nèi)，GBDT 和LGB 模型關(guān)于三峽水庫逐月流量的模擬結(jié)果精度較為接近，其中MAPE值分別為16.1%和15.7%，CC值均為0.91，NSE值分別為0.84 和0.86；而驗(yàn)證期內(nèi)，除CC值較為接近外，其他兩種指標(biāo)差異顯著增大，其中MAPE值分別為25.8%和23.0%，NSE值分別為0.71 和0.75。整體來看，兩種模型均具有較好的預(yù)測精度，與GBDT 模型相比，LGB 模型在率定期雖然與前者精度差異不大，但在驗(yàn)證期精度明顯更高。

圖1 提供了三峽水庫實(shí)測流量與模擬流量過程對比，其結(jié)果表明：率定期內(nèi)，豐水年的模擬值偏小，尤其對于1997 和1998年兩場特大洪水而言，汛期模擬值與實(shí)測值間絕對偏差超過30%，而對于平水年和枯水年的模擬結(jié)果略偏大；驗(yàn)證期內(nèi)，兩種模型模擬結(jié)果整體大于實(shí)測系列，但偏差程度不大。圖2提供兩種模型模擬結(jié)果與實(shí)測系列間散點(diǎn)分布圖，依據(jù)其特征可以得出同樣的結(jié)論，率定期內(nèi)，兩種模型位于下半部分的散點(diǎn)均勻地分布在45°線兩側(cè)，而上半部分散點(diǎn)則略低于45°線，驗(yàn)證期下半部分散點(diǎn)分布特征與率定期類似，而上半部分散點(diǎn)則略高于45°線。對比兩種模型結(jié)果可發(fā)現(xiàn)，圖1 與圖2 均表明：與GBDT模型相比，LGB模型對于三峽水庫逐月流量過程的模擬性能更佳，特別是對汛期而言，其模擬結(jié)果與實(shí)測值的偏差程度更小，對應(yīng)散點(diǎn)也更為接近45°線。

圖2 三峽水庫實(shí)測與模擬月流量散點(diǎn)圖Fig.2 Scatter plots of observed and predicted monthly runoff for the Three Gorges Reservoir

2.3 不確定性分析

HUP 模型處理過程中，需先采用正態(tài)分位數(shù)變換將實(shí)測值與預(yù)測值分別轉(zhuǎn)換至正態(tài)空間中，然后依據(jù)貝葉斯理論與線性似然假設(shè)來推求實(shí)際流量的后驗(yàn)概率密度函數(shù)，從而對預(yù)報(bào)結(jié)果進(jìn)行不確定性分析。與上述預(yù)測過程類似，研究將1965-2001年作為率定期，2002-2016年作為驗(yàn)證期，來對HUP 模型后處理精度進(jìn)行相應(yīng)驗(yàn)證。

依據(jù)HUP模型確定的實(shí)際流量后驗(yàn)概率密度分布，可以提取逐月徑流過程對應(yīng)的50%分位數(shù)流量值（Q50值）。表3和圖3 分別提供了Q50 值的精度指標(biāo)及其與實(shí)測過程的對比。兩者結(jié)果表明：與LGB 模型預(yù)報(bào)結(jié)果相比，Q50 值與實(shí)測系列間擬合程度更佳，各精度指標(biāo)在率定期及驗(yàn)證期也顯著優(yōu)于LGB 模型結(jié)果。為進(jìn)一步分析兩者精度差異的年內(nèi)變化特征，研究還計(jì)算了驗(yàn)證期內(nèi)實(shí)測值與模擬值間的相對誤差，如圖4 所示。由圖4，不難看出，兩者在枯水期（1-4月及12月）相對誤差普遍偏小，大體上處于-35%～35%之間，而在汛期，模擬值較實(shí)測值明顯偏大，在部分年份的8-10月內(nèi)，相對誤差甚至達(dá)到了90%，而對比圖4（a）與圖4（b），不難看出與LGB 模型預(yù)測結(jié)果相比，Q50 值在大部分月份的預(yù)測精度呈現(xiàn)一定程度提高，特別是對汛期而言，改善效果更為顯著。

圖3 三峽水庫1965-2016年逐月流量觀測系列與Q50系列對比圖Fig.3 The observed runoff and Q50 values series for the Three Gorges Reservoir

圖4 驗(yàn)證期逐月流量模擬結(jié)果相對誤差圖Fig.4 The relative error of predicted results for the validation period

表3 HUP模型率定期及驗(yàn)證期精度評價(jià)（Q50預(yù)報(bào)值）Tab.3 Performance Indices for the calibration and validation periods（Q50 predictions）

除Q50 預(yù)測值外，研究還提取了預(yù)測結(jié)果的90%不確定性區(qū)間，其中驗(yàn)證期不確定性區(qū)間如圖5 所示。為定量評估預(yù)測結(jié)果的不確定性，研究選取了覆蓋率（CR），相對帶寬（RB）和平均偏移（RD）幅度分別對率定期及驗(yàn)證期的不確定性區(qū)間進(jìn)行相應(yīng)分析，對應(yīng)指標(biāo)結(jié)果如表4 所示。由表4 及圖5 結(jié)果，不難看出提取的90%不確定性區(qū)間覆蓋了率定期與驗(yàn)證期絕大多數(shù)的實(shí)測點(diǎn)據(jù)，兩者覆蓋率分別達(dá)到0.86 和0.83，均較為接近0.9，說明提取的置信區(qū)間是較為可靠的。而相對帶寬和平均偏移幅度結(jié)果表明，與率定期相比，驗(yàn)證期的相對帶寬和平均偏移幅度更大，說明后者的預(yù)測精度要差于前者，且對應(yīng)不確定性也更大。

圖5 三峽水庫月流量不確定性分析結(jié)果Fig.5 Uncertainty analysis of the predicted monthly runoff sequences of the Three Gorges Reservoir

表4 不確定性區(qū)間精度評價(jià)Tab.4 Performance Indices for uncertainty confidence intervals

3 結(jié) 論

（1）預(yù)測精度指標(biāo)結(jié)果表明，GBDT 模型與LGB 模型對于三峽水庫逐月流量過程的預(yù)測結(jié)果在率定期和驗(yàn)證期內(nèi)均取得較好效果，但后者的精度要顯著優(yōu)于前者。

（2）采用HUP 模型對LGB 模型預(yù)測結(jié)果進(jìn)行后處理，提取的Q50 值系列對應(yīng)的MAPE、CC及NSE值均優(yōu)于LGB 模型預(yù)測結(jié)果，說明HUP模型具有一定程度的校正能力。

（3）90%不確定性區(qū)間在率定期及驗(yàn)證期的覆蓋率分別為0.86和0.83，較為接近0.9，且相對帶寬和平均偏移幅度較小，說明基于HUP模型提取的90%不確定性區(qū)間是較為合理可靠的。□