陳 雪 怡，陳 元 芳，王 文 鵬，邱 鵬

(1.河海大學(xué) 水文水資源學(xué)院，江蘇 南京 210098； 2.重慶市水文監(jiān)測總站，重慶 401147)

0 引 言

準(zhǔn)確的中長期水文預(yù)報(bào)對流域水旱災(zāi)害防御、水資源綜合利用、河湖水環(huán)境保護(hù)等水利業(yè)務(wù)工作具有重要意義[1]。根據(jù)預(yù)報(bào)變量的不同，徑流預(yù)報(bào)可以劃分為兩類[2]：一類是基于單變量的水文預(yù)報(bào)，即發(fā)掘水文要素自身演變規(guī)律，常見的方法有時(shí)間序列法[3-4]、周期疊加法等；另一類是基于多變量的水文預(yù)報(bào)，即發(fā)掘大氣環(huán)流、海陸狀況等成因變量與水文要素之間的聯(lián)系，常見的方法有隨機(jī)森林[5]、人工神經(jīng)網(wǎng)絡(luò)、支持向量機(jī)[6]等。近年來，這兩類方法在中長期徑流預(yù)報(bào)中都有廣泛的應(yīng)用[7]。Fan等[8]建立自回歸求和滑動(dòng)平均模型預(yù)測了塔里木河干流徑流量，并分析了溫度與降水對其的影響。Wang等[9]將經(jīng)驗(yàn)?zāi)B(tài)分解法與時(shí)間序列法相結(jié)合，對黃河上游進(jìn)行長期徑流預(yù)測。林劍藝等[10]將支持向量機(jī)應(yīng)用于中長期預(yù)報(bào)中，在建模過程中引入徑向基核函數(shù)簡化非線性問題的求解，并應(yīng)用SCE-UA算法辨識(shí)支持向量機(jī)的參數(shù)。趙文秀等[11]采用隨機(jī)森林模型篩選預(yù)報(bào)因子，將篩選后的因子作為RBF神經(jīng)網(wǎng)絡(luò)的輸入層進(jìn)行長期徑流預(yù)報(bào)。綜上，國內(nèi)外學(xué)者對兩類預(yù)報(bào)方法均進(jìn)行了大量探索，但何種方法更優(yōu)并沒有統(tǒng)一的結(jié)論。實(shí)際上，在具體流域的預(yù)報(bào)工作中，仍有必要綜合考慮預(yù)報(bào)要素的水文特性、預(yù)報(bào)因子的可獲得性、預(yù)報(bào)方法的穩(wěn)健性等因素，通過定量對比分析選取合適的方法建立預(yù)報(bào)模型，實(shí)現(xiàn)精度可靠的中長期水文預(yù)報(bào)。

因此，本文以嘉陵江流域北碚水文站為例，建立基于時(shí)間序列模型與成因驅(qū)動(dòng)的徑流預(yù)報(bào)模型；同時(shí)應(yīng)用時(shí)變權(quán)重方法，集合單一模型進(jìn)行組合預(yù)報(bào)，通過對比單變量與多變量、單一與組合模型預(yù)報(bào)結(jié)果，分析引入多變量與組合模型是否有利于提高徑流預(yù)報(bào)容錯(cuò)能力、穩(wěn)定性與精度，研究成果可為徑流預(yù)報(bào)的建模策略提供參考。

1 研究方法

本文選擇時(shí)間序列法作為單變量預(yù)報(bào)方法，隨機(jī)森林、BP神經(jīng)網(wǎng)絡(luò)與多元線性回歸作為多變量預(yù)報(bào)方法，并采用組合預(yù)報(bào)方法對單一模型的權(quán)重進(jìn)行實(shí)時(shí)調(diào)整從而對月徑流預(yù)報(bào)結(jié)果進(jìn)行優(yōu)化。方法流程如圖1所示。

1.1 水文時(shí)間序列預(yù)報(bào)方法

時(shí)間序列分析方法最早由耶爾于1927年提出。水文時(shí)間序列預(yù)報(bào)就是依據(jù)水文過程在相鄰時(shí)間截口之間時(shí)序相依的統(tǒng)計(jì)特征，建立由序列前期值驅(qū)動(dòng)的水文狀態(tài)量預(yù)報(bào)模型。對于平穩(wěn)的時(shí)間序列，可以構(gòu)建自回歸模型(Autoregressive Model，AR)、滑動(dòng)平均模型(Moving Average Model，MA)與自回歸滑動(dòng)平均模型(Auto-Regressive Moving Average Model，ARMA)[12]。自回歸滑動(dòng)平均模型，記為ARMA(p，q)，具有以下結(jié)構(gòu)：

(1)

式中：xt為平穩(wěn)時(shí)間序列；εt為白噪聲序列；φp，θq為模型參數(shù)，常用尤爾-沃爾克法估計(jì)；對于中心化序列xt，φ0=0。

對于非平穩(wěn)的水文序列，通過差分運(yùn)算將原序列轉(zhuǎn)化為平穩(wěn)序列，再通過白噪聲檢驗(yàn)后與ARMA模型組合，得到自回歸滑動(dòng)平均求和模型(Autoregressive Integrated Moving Average Model，ARIMA)，記為ARIMA(p，d，q)。

1.2 成因驅(qū)動(dòng)的多變量預(yù)報(bào)方法

基于成因驅(qū)動(dòng)的多變量預(yù)報(bào)模型是從徑流成因的角度構(gòu)建能夠反映驅(qū)動(dòng)因子與預(yù)報(bào)要素之間關(guān)系的多元相關(guān)預(yù)報(bào)模型，應(yīng)用此類方法時(shí)驅(qū)動(dòng)因子的篩選與預(yù)報(bào)模型的選擇至關(guān)重要。首先根據(jù)徑流的物理成因選擇88項(xiàng)大氣環(huán)流指數(shù)、26項(xiàng)海溫指數(shù)、16項(xiàng)其他氣象指數(shù)和月降水量作為初選影響因子，計(jì)算前一水文年各因子與徑流量的相關(guān)系數(shù)，剔除無物理背景、不顯著相關(guān)的影響因子并進(jìn)行排序。對排序前50的影響因子構(gòu)建隨機(jī)森林模型，計(jì)算未被袋裝法選中的樣本(OOB數(shù)據(jù))打亂前后均方誤差的增加量(IncMSE)來定量評價(jià)因子的重要性[13]，取IncMSE排位靠前的因子作為預(yù)報(bào)模型輸入變量。本文對比分析的模型為隨機(jī)森林、BP神經(jīng)網(wǎng)絡(luò)與多元線性回歸。

隨機(jī)森林(Random Forest，RF)最早由Breiman于21世紀(jì)初提出[14]。隨機(jī)森林由多棵決策樹構(gòu)成，其基本思想為：從原始樣本中隨機(jī)多次抽取樣本，對每個(gè)重抽樣本都建立決策樹模型，組合多棵決策樹對應(yīng)變量的估計(jì)結(jié)果，當(dāng)應(yīng)變量的綜合估計(jì)結(jié)果在樣本訓(xùn)練集和驗(yàn)證集表現(xiàn)俱佳時(shí)，就構(gòu)成了隨機(jī)森林模型[15-16]。

BP神經(jīng)網(wǎng)絡(luò)模型具有3層結(jié)構(gòu)：輸入層、隱含層和輸出層。輸入層中每個(gè)神經(jīng)元對應(yīng)一項(xiàng)驅(qū)動(dòng)因子，隱含層神經(jīng)元個(gè)數(shù)通過“試錯(cuò)法”調(diào)整，輸出層神經(jīng)元個(gè)數(shù)為預(yù)報(bào)對象個(gè)數(shù)。將正向誤差反向傳播至網(wǎng)絡(luò)層，調(diào)整連接節(jié)點(diǎn)權(quán)重，不斷重復(fù)直至輸出結(jié)果誤差函數(shù)達(dá)到最小[17-18]。

多元線性回歸就是建立多個(gè)自變量與因變量之間的回歸方程，一般以誤差平方和為目標(biāo)使用最小二乘法估計(jì)參數(shù)，數(shù)學(xué)表達(dá)式為

y=α0+α1x1+α2x2+…+αnxn+ε

(2)

式中：αn為回歸系數(shù)；xn為驅(qū)動(dòng)因子值；ε為模型誤差。

1.3 時(shí)變權(quán)重組合預(yù)測

權(quán)重組合預(yù)測將“好而不同”的單一模型預(yù)報(bào)結(jié)果加權(quán)求和，以獲得更準(zhǔn)確的預(yù)報(bào)效果。根據(jù)權(quán)重是否隨時(shí)間改變，可以分為定常權(quán)重與時(shí)變權(quán)重兩類。定常權(quán)重是指保持各單一模型的權(quán)重不變進(jìn)行組合，該類方法的研究較成熟，但在某些時(shí)段組合預(yù)測精度較單一模型更差。時(shí)變權(quán)重則隨著不同模型性能的變化而進(jìn)行權(quán)重調(diào)整，通過動(dòng)態(tài)預(yù)報(bào)有效規(guī)避定常權(quán)重組合預(yù)測常介于最優(yōu)與最劣模型間的問題[19]。本文選擇時(shí)變權(quán)重方法中的時(shí)變誤差平方和(TSSE)方法，根據(jù)預(yù)報(bào)前序時(shí)段的誤差來確定局部預(yù)報(bào)時(shí)段的動(dòng)態(tài)權(quán)重，公式如下：

(3)

(4)

1.4 預(yù)測性能評價(jià)指標(biāo)

采用合格率(QR，預(yù)報(bào)值與實(shí)際值相對誤差在±20%之間為合格)、平均絕對百分比誤差(MAPE)、平均相對誤差(MRE)、均方根誤差(RMSE)和決定系數(shù)(R2)5項(xiàng)指標(biāo)作為評價(jià)模型預(yù)報(bào)精度的指標(biāo)。采用均方根誤差變化量(ΔRMSE)作為評價(jià)模型穩(wěn)定性的指標(biāo)。具體計(jì)算公式見文獻(xiàn)[19]。

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

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

本文預(yù)報(bào)對象為嘉陵江北碚站汛期月徑流量，其月徑流量數(shù)據(jù)來源于《中華人民共和國水文年鑒》。計(jì)算武都、略陽、廣元、萬源、閬中、巴中、達(dá)縣、遂寧、高坪區(qū)和沙坪壩10個(gè)氣象站(見圖2)的算術(shù)平均降水量作為代表降水量，數(shù)據(jù)來源于中國氣象數(shù)據(jù)網(wǎng)。收集全球逐月百項(xiàng)氣象指數(shù)，數(shù)據(jù)來源于國家氣候中心。本文以1979～2006年為訓(xùn)練期，以2007～2019年為測試期。

由于嘉陵江中下游地區(qū)的大規(guī)模水土流失綜合治理、水庫水電站等大型水利設(shè)施的修建等人類活動(dòng)對下墊面條件的影響，嘉陵江的實(shí)測徑流序列發(fā)生了突變，不能代表天然徑流狀況。因此考慮方差修正的預(yù)置白方法和Mann-Kendall檢驗(yàn)對北碚站訓(xùn)練期降水序列與徑流序列進(jìn)行分析與一致性修正[20-22]，發(fā)現(xiàn)北碚站汛期(6～8月)降水量沒有明顯突變，而徑流量在1986年發(fā)生了突變，如圖3所示。

繪制訓(xùn)練期徑流量-降水量雙累計(jì)曲線圖，根據(jù)突變前雙累計(jì)曲線回歸方程對突變后徑流進(jìn)行修正，如圖4所示。按照修正前后總徑流量的比例對各月徑流量同倍比縮放，生成還原后的天然月徑流量序列。對修正后的月徑流量序列再進(jìn)行突變檢驗(yàn)，未發(fā)現(xiàn)新的突變點(diǎn)，滿足一致性條件。

2.2 單變量與多變量模型預(yù)報(bào)效果

以經(jīng)過還原得到的北碚站6～8月天然徑流量作為目標(biāo)變量，構(gòu)建預(yù)報(bào)模型。對于單變量預(yù)報(bào)模型，參數(shù)p、q根據(jù)d階差分序列的自相關(guān)圖和偏相關(guān)圖初選，再結(jié)合模型AIC準(zhǔn)則進(jìn)行確定，分別為ARIMA(2，2，0)、ARIMA(6，2，0)和ARIMA(1，1，0)。對于多變量預(yù)報(bào)模型，對初選因子按照IncMSE指標(biāo)進(jìn)行排序，選擇指標(biāo)值明顯大于其他因子的因子集作為預(yù)報(bào)因子，6月和8月選定前8項(xiàng)因子，7月選定前7項(xiàng)因子，作為模型的輸入變量，其具體名稱如表1所列。經(jīng)過反復(fù)測試，增加剩余因子不足以明顯改善解釋能力，反而增加模型復(fù)雜度，降低模型的泛化能力。構(gòu)建隨機(jī)森林模型時(shí)，對隨機(jī)森林模型的兩個(gè)待選參數(shù)進(jìn)行優(yōu)化，即子決策樹數(shù)量Ntree與回歸樹節(jié)點(diǎn)劃分?jǐn)?shù)Mtry。挑選使OOB數(shù)據(jù)泛化誤差最小時(shí)的Mtry值作為最優(yōu)參數(shù)，確定Mtry后代入使模型誤差基本穩(wěn)定的Ntree值作為最優(yōu)參數(shù)。本文中各月隨機(jī)森林模型Mtry分別取4，4和3，Ntree均取1 000。構(gòu)建BP神經(jīng)網(wǎng)絡(luò)模型時(shí)，隱含層神經(jīng)元個(gè)數(shù)用試錯(cuò)法比選，使得訓(xùn)練誤差和測試誤差最小，網(wǎng)絡(luò)結(jié)構(gòu)分別為8-4-1，7-5-1和8-5-1。

表1 北碚站月徑流預(yù)報(bào)因子Tab.1 Monthly runoff forecast factors of Beibei Station

單變量與多變量預(yù)測模型的散點(diǎn)圖如圖5所示，4種模型的預(yù)測值都均勻分布在1∶1線的兩側(cè)，且在不同月份內(nèi)同一模型都表現(xiàn)出相似的準(zhǔn)確性與適應(yīng)性。對于不同變量的預(yù)報(bào)模型，基于多變量的預(yù)報(bào)模型精度優(yōu)于基于單變量的預(yù)報(bào)模型。相對而言，ARIMA預(yù)報(bào)徑流量的點(diǎn)據(jù)分布更為分散，說明其擬合程度相比多變量預(yù)報(bào)模型較差?；趦煞N模式下的徑流預(yù)報(bào)的MAPE箱形圖如圖6所示。由圖可知，單變量ARIMA模型的MAPE值中位數(shù)約為40%，而3種多變量模型的MAPE值中位數(shù)均小于20%，說明多變量預(yù)測模型具有更好的預(yù)報(bào)精度。因此，從徑流成因角度盡可能地挖掘潛在的氣候驅(qū)動(dòng)因子，增加預(yù)報(bào)模型的信息量，對提高徑流預(yù)報(bào)精度是有利的。

綜合對比3種多變量模型的預(yù)報(bào)性能，如表2所列。從預(yù)測精度上看，RF在訓(xùn)練期的合格率最高、MRE與RMSE最低，在不同時(shí)段下合格率達(dá)到了82%以上、MRE在13%以下；BP在測試期的合格率最高、MRE與RMSE最低，在不同時(shí)段下合格率達(dá)到了76%以上、MRE在18%以下?？傮w上RF與BP模型表現(xiàn)均優(yōu)于多元線性回歸模型。從數(shù)據(jù)離散程度上看，RF的R2最接近1，說明其擬合程度最優(yōu)。從模型穩(wěn)定性上看，BP的ΔRMSE最小，7月與8月測試期的RMSE小于訓(xùn)練期，指示該模型具有較強(qiáng)的泛化能力?？傊?，在多變量預(yù)測模型中RF和BP神經(jīng)網(wǎng)絡(luò)模型的預(yù)報(bào)性能優(yōu)于傳統(tǒng)的多元線性回歸模型。

表2 不同模型預(yù)報(bào)精度Tab.2 Forecast accuracy of different models

2.3 組合模型預(yù)報(bào)效果

目前能夠公開獲取且滿足時(shí)效性要求的徑流成因要素還比較有限，相關(guān)預(yù)報(bào)模型的建模依據(jù)仍以小容量樣本為主，模型泛化容易出現(xiàn)過擬合情況。根據(jù)統(tǒng)計(jì)學(xué)原理，集成單一模型的預(yù)報(bào)結(jié)果，可以減少模型泛化的方差，有效降低預(yù)報(bào)不確定性。因此，選擇對序列趨勢模擬較準(zhǔn)確的3種多變量模型進(jìn)行組合預(yù)報(bào)，根據(jù)前文權(quán)重計(jì)算方法得到權(quán)重取值如圖7所示。總體而言，RF模型的預(yù)報(bào)誤差較小，組合權(quán)重占優(yōu)，但隨著預(yù)報(bào)時(shí)刻推移，各模型的組合權(quán)重會(huì)沿程發(fā)生變化。多元線性回歸模型的權(quán)重呈現(xiàn)出平穩(wěn)的特征，而兩種機(jī)器學(xué)習(xí)模型的權(quán)重在各月都表現(xiàn)出非平穩(wěn)特征，且在訓(xùn)練期和測試期都有明顯變化趨勢。RF模型的權(quán)重呈現(xiàn)先增大后減小的趨勢，說明其預(yù)報(bào)誤差在訓(xùn)練期呈現(xiàn)減小趨勢，到了測試期開始呈現(xiàn)增大趨勢，而BP模型則與之相反，這也驗(yàn)證了RF模型泛化能力差而BP模型穩(wěn)定性強(qiáng)的特點(diǎn)。

圖8比較了3種多變量模型和組合模型的預(yù)報(bào)性能。經(jīng)組合預(yù)報(bào)后，訓(xùn)練期合格率達(dá)到了89%以上、MRE在13%以下；測試期合格率達(dá)到了84%以上、MRE在15%以下。合格率在訓(xùn)練期與測試期均得到了明顯提高，高于或接近于最優(yōu)的單一模型；定量誤差在7月測試期時(shí)劣于BP神經(jīng)網(wǎng)絡(luò)模型，在其他月份各時(shí)期都有明顯下降趨勢，低于或接近于單一最優(yōu)模型。因此，時(shí)變權(quán)重組合預(yù)報(bào)方法能有效集合單一預(yù)測方法提供的信息，結(jié)合各模型的優(yōu)點(diǎn)，進(jìn)一步提高基于多變量的預(yù)測模型的精度。

3 結(jié) 論

本文以北碚站月徑流量預(yù)報(bào)為例，對比不同建模方法的預(yù)報(bào)性能，得到以下結(jié)論。

(1) 從單變量與多變量模型選擇上看：多變量預(yù)報(bào)法優(yōu)于單變量預(yù)報(bào)法。多變量預(yù)報(bào)模型的MAPE中位數(shù)均小于20%，說明從徑流成因角度挖掘潛在預(yù)報(bào)因子、增加樣本信息量，是提高月徑流預(yù)報(bào)精度的有效途徑之一。

(2) 從不同多變量模型表現(xiàn)上看：RF與BP模型在定量誤差與模型穩(wěn)定性方面均表現(xiàn)較好，說明針對成因機(jī)制復(fù)雜，樣本空間非線性可分的預(yù)報(bào)對象，非線性映射模型常優(yōu)于線性回歸模型。

(3) 從單一與組合模型上看：時(shí)變權(quán)重能有效集合各單一模型的優(yōu)勢，使組合后的預(yù)報(bào)精度在訓(xùn)練期與測試期都高于或接近于最優(yōu)的單一模型，在有效減小預(yù)報(bào)誤差的同時(shí)提高了模型的穩(wěn)定性。

本文探討了天然徑流量預(yù)報(bào)的建模方法，但在實(shí)際預(yù)報(bào)業(yè)務(wù)中，仍需進(jìn)一步從水保措施治理面積、大中型水利工程建設(shè)數(shù)量、用水量、水庫蓄水量等因素，定量分析人類活動(dòng)對月徑流量的影響及其與時(shí)間的關(guān)系，將天然徑流量的預(yù)報(bào)結(jié)果與人類活動(dòng)影響導(dǎo)致的變化量相結(jié)合，得到月徑流過程的實(shí)際變化，以提高預(yù)報(bào)模型的實(shí)用性。