月徑流預報建模方法對比分析
——以嘉陵江北碚站為例

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

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

0 引 言

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

因此，本文以嘉陵江流域北碚水文站為例，建立基于時間序列模型與成因驅動的徑流預報模型；同時應用時變權重方法，集合單一模型進行組合預報，通過對比單變量與多變量、單一與組合模型預報結果，分析引入多變量與組合模型是否有利于提高徑流預報容錯能力、穩定性與精度，研究成果可為徑流預報的建模策略提供參考。

1 研究方法

本文選擇時間序列法作為單變量預報方法，隨機森林、BP神經網絡與多元線性回歸作為多變量預報方法，并采用組合預報方法對單一模型的權重進行實時調整從而對月徑流預報結果進行優化。方法流程如圖1所示。

1.1 水文時間序列預報方法

時間序列分析方法最早由耶爾于1927年提出。水文時間序列預報就是依據水文過程在相鄰時間截口之間時序相依的統計特征，建立由序列前期值驅動的水文狀態量預報模型。對于平穩的時間序列，可以構建自回歸模型(Autoregressive Model，AR)、滑動平均模型(Moving Average Model，MA)與自回歸滑動平均模型(Auto-Regressive Moving Average Model，ARMA)[12]。自回歸滑動平均模型，記為ARMA(p，q)，具有以下結構：

(1)

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

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

1.2 成因驅動的多變量預報方法

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

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

BP神經網絡模型具有3層結構：輸入層、隱含層和輸出層。輸入層中每個神經元對應一項驅動因子，隱含層神經元個數通過“試錯法”調整，輸出層神經元個數為預報對象個數。將正向誤差反向傳播至網絡層，調整連接節點權重，不斷重復直至輸出結果誤差函數達到最小[17-18]。

多元線性回歸就是建立多個自變量與因變量之間的回歸方程，一般以誤差平方和為目標使用最小二乘法估計參數，數學表達式為

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

(2)

式中：αn為回歸系數；xn為驅動因子值；ε為模型誤差。

1.3 時變權重組合預測

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

(3)

(4)

1.4 預測性能評價指標

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

2 實例應用

2.1 研究區域與數據

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

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

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

2.2 單變量與多變量模型預報效果

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

表1 北碚站月徑流預報因子Tab.1 Monthly runoff forecast factors of Beibei Station

單變量與多變量預測模型的散點圖如圖5所示，4種模型的預測值都均勻分布在1∶1線的兩側，且在不同月份內同一模型都表現出相似的準確性與適應性。對于不同變量的預報模型，基于多變量的預報模型精度優于基于單變量的預報模型。相對而言，ARIMA預報徑流量的點據分布更為分散，說明其擬合程度相比多變量預報模型較差。基于兩種模式下的徑流預報的MAPE箱形圖如圖6所示。由圖可知，單變量ARIMA模型的MAPE值中位數約為40%，而3種多變量模型的MAPE值中位數均小于20%，說明多變量預測模型具有更好的預報精度。因此，從徑流成因角度盡可能地挖掘潛在的氣候驅動因子，增加預報模型的信息量，對提高徑流預報精度是有利的。

綜合對比3種多變量模型的預報性能，如表2所列。從預測精度上看，RF在訓練期的合格率最高、MRE與RMSE最低，在不同時段下合格率達到了82%以上、MRE在13%以下；BP在測試期的合格率最高、MRE與RMSE最低，在不同時段下合格率達到了76%以上、MRE在18%以下。總體上RF與BP模型表現均優于多元線性回歸模型。從數據離散程度上看，RF的R2最接近1，說明其擬合程度最優。從模型穩定性上看，BP的ΔRMSE最小，7月與8月測試期的RMSE小于訓練期，指示該模型具有較強的泛化能力。總之，在多變量預測模型中RF和BP神經網絡模型的預報性能優于傳統的多元線性回歸模型。

表2 不同模型預報精度Tab.2 Forecast accuracy of different models

2.3 組合模型預報效果

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

圖8比較了3種多變量模型和組合模型的預報性能。經組合預報后，訓練期合格率達到了89%以上、MRE在13%以下；測試期合格率達到了84%以上、MRE在15%以下。合格率在訓練期與測試期均得到了明顯提高，高于或接近于最優的單一模型；定量誤差在7月測試期時劣于BP神經網絡模型，在其他月份各時期都有明顯下降趨勢，低于或接近于單一最優模型。因此，時變權重組合預報方法能有效集合單一預測方法提供的信息，結合各模型的優點，進一步提高基于多變量的預測模型的精度。

3 結 論

本文以北碚站月徑流量預報為例，對比不同建模方法的預報性能，得到以下結論。

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

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

(3) 從單一與組合模型上看：時變權重能有效集合各單一模型的優勢，使組合后的預報精度在訓練期與測試期都高于或接近于最優的單一模型，在有效減小預報誤差的同時提高了模型的穩定性。

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