基于多模型融合的中長期徑流集成預測方法

朱非林，陳嘉乙，張 咪，徐向榮，鐘平安

(河海大學水文水資源學院，江蘇 南京 210024)

0 引 言

中長期徑流預報在水資源開發調度與管理、水利工程運行與維護以及水旱災害防御等領域中扮演著重要的技術角色。近年來，隨著極端洪澇干旱事件的頻發，它也成為水文領域研究應用中備受關注的難點和熱點問題。然而，受全球范圍內氣候和大氣環流等復雜要素的綜合影響，中長期徑流預報缺乏可靠的氣象預報信息，導致預報結果的不確定性較大。因此，如何提高預報精度已成為當前研究的重點[1]。

近年來，中長期水文預報研究取得了諸多進展。根據研究方法的不同，預報模型可分為考慮徑流時間序列演變特征的序列演變方法和考慮預報因子與徑流數量關系的因子關系方法[2]。通過機器學習方法建立預報模型成為當前主流的研究途徑，在單一模型預測方面，岳兆新等[3]基于改進深度信念網絡模型對雅礱江流域進行中長期徑流預測，實現更好的預測速度和精度；楊文發等[4]使用流域多尺度水文預報模型，探討了流域水文預報中涉及的不確定性因素以及水文氣象耦合方面的不匹配性；邵健偉等[5]提出了基于SpringBoot框架的中長期水文預報系統，實現了中長期水文預報業務化運行。在多模型集成預測方面，李宏亮[6]采用加權平均法組合3種不同預報模型實現水文集合預報；周研來等[7]采用Gamma Test數據驅動模型與長短期記憶神經網絡結合3種機器學習算法的綜合模型，以解決在變化環境下降雨-洪水過程統計的非線性、隨機性和時變性問題；Tan等[8]采用EEMD-ANN混合方法建立適合長江流域汛期月徑流的預測模型；Liu等[9]將隱馬爾可夫模型與高斯混合回歸相結合，證明其可以處理多模態和異方差數據。然而，單一預報模型雖然可以通過特征提取、參數優化等方式提高預測精度并縮短訓練時間，但其模型結構相對固定，參數選擇范圍有限，考慮因子無法適應不同流域。現有集成預測方法大多研究多個確定性模型的加權組合，但由于不同模型數據結構、參數形式和變量分布具有復雜性，簡單的隨機加權容易導致過擬合現象。

鑒于此，本文提出了一種基于多模型融合的水庫中長期徑流集成預測方法。采用Stacking融合算法，選取ARMA、BP、LSTM、RF和SVR等5個異質預測模型進行融合，并通過超參數尋優方法對這些模型進行參數優化。以黃河上游龍羊峽水庫作為研究對象，驗證了該方法在提高中長期徑流預報精度方面的有效性。

1 研究方法

圖1為基于多模型融合的中長期徑流集成預測框架流程圖。首先，采用相關性分析確定與徑流顯著相關的因子，并利用主成分分析降低數據維度，完成預報因子的篩選。然后，使用不同累積貢獻率的主成分組對應的基學習器進行訓練，并采用超參數優化尋優，得到各基學習器的最優主成分結果。最后，通過K折交叉驗證處理最優主成分結果，訓練元學習器預測模型，得到多模型融合的預測結果。

圖1 基于多模型融合的中長期徑流集成預測框架

1.1 預報因子的篩選

1.1.1 初始預報因子集

中長期徑流變化的復雜性使其影響因子的類型較多。本文從影響中長期徑流變化的機制出發，選擇前期徑流變量、大氣環流因子、水文要素因子3類因子作為研究數據。具體包括：1959年～2016年龍羊峽水庫逐月入庫徑流量；1959年～2016年中國氣象局國家氣候中心88項大氣環流指數逐月數據；1959年～2016年青海省西寧氣象站逐月降水數據。

1.1.2 相關性分析

相關性分析是一種常用的統計分析方法，用于研究兩個或多個變量之間的相關性或關聯性，甄別對徑流變化貢獻大、作用強的影響因子。本文采用Spearman相關系數來衡量徑流量與影響因子之間的相關程度，其公式為

(1)

1.1.3 主成分分析

主成分分析(PCA)是一種常用的多元統計分析方法，該方法通過對一組變量的線性組合來尋找變量間的主要模式和結構，將多個變量轉化為少數幾個主成分，在降低維度的同時使保留的信息量最大化。

假設原始樣本為X={x1，x2，…，xm}，xi∈Rn，樣本數m個，每個樣本有n個特征維度。對樣本進行標準化變換，得標準化矩陣Z。

(2)

對標準化矩陣Z求相關系數矩陣，rij為兩變量之間的相關系數，相關系數矩陣R為

R=[rij]n

(3)

解樣本相關系數矩陣R的特征方程|R-λIn|=0，λ為特征根，In為特征向量，按不同方差貢獻率a確定各因子的信息利用率，即

(4)

(5)

1.2 單一預測模型

本文選擇在水文預報領域經過大量實踐驗證的5種具有典型性、代表性、異質性的單一模型，分別是屬于監督學習類的支持向量回歸模型(Support vector regression，SVR)、BP神經網絡模型(BP neural network，BP)；屬于集成學習類的隨機森林模型(Random forest，RF)；屬于時間序列分析類的自回歸滑動平均模型(Autoregressive moving average model，ARMA)；屬于深度學習類的長短期記憶神經網絡模型(Long short-term memory，LSTM)。

1.2.1 自回歸滑動平均模型

ARMA模型是一種基于時間序列分析的預測模型，結合了自回歸(AR)和移動平均(MA)模型，可被用于描述時間序列中的自相關和移動平均關系。

AR模型yt=c+∑(φiyt-i)+εt

(6)

MA模型yt=μ+εt+∑(θtεt-i)

(7)

ARMA模型

yt=c+∑(φiyt-i)+εt+∑(θtεt-i)

(8)

式中，yt為時間序列的觀測值；c為常數項；φi和θt分別為AR和MA模型的系數；εt為誤差項；μ為時間序列的均值。

1.2.2 BP神經網絡模型

BP預測模型是一種基于人工神經網絡的機器學習算法，用于解決回歸和分類問題。它通過反向傳播算法來訓練神經網絡，不斷調整權重和偏置使預測值與實際值之間的誤差最小化。模型結構如圖2所示。

圖2 BP神經網絡模型結構

1.2.3 長短期記憶神經網絡模型

長短期記憶網絡LSTM是循環神經網絡(RNN)的一種改進算法。LSTM的細胞單元包括輸入門、遺忘門、輸出門，其中輸入門用于控制信息輸入，遺忘門用于控制是否遺忘歷史序列信息，輸出門用于控制信息輸出。

遺忘門中輸入xt與狀態記憶單元St-1、中間輸出ht-1共同決定狀態記憶單元遺忘部分。輸入門中的xt分別經過sigmoid和tanh函數變化后共同決定狀態記憶單元中保留向量。中間輸出ht由更新后的St與輸出ot共同決定，計算公式為

ft=σ(Wfxxt+Wfhht-1+bf)

(9)

it=σ(Wixxt+Wihht-1+bi)

(10)

gt=σ(Wgxxt+Wghht-1+bg)

(11)

ot=σ(Woxxt+Wohht-1+bo)

(12)

St=gt⊙it+St-1⊙ft

(13)

ht=φ(St)⊙ot

(14)

式中，ft、it、gt、ot、ht和St分別為遺忘門、輸入門、輸入節點、輸出門、中間輸出和細胞單元的狀態；Wfx、Wfh、Wix、Wih、Wgx、Wgh、Wox、Woh分別為相應門與輸入xt和中間輸出ht-1相乘的矩陣權重；bf、bi、bg、bo分別為相應門的偏置項；⊙為向量中元素按位相乘；σ為sigmoid函數變化；φ為tanh函數變化。

1.2.4 隨機森林模型

隨機森林是一種基于Bagging思想和特征子空間思想的集成學習模型，以多個由Bagging技術訓練得到的決策樹作為基學習器，組合單個決策樹的輸出結果，投票決定最終分類結果。隨機森林模型具有較強的非線性模擬能力，不易出現過擬合現象，對數據集特征的挖掘的魯棒性強[10]。

通過K輪訓練，以Bagging抽樣技術從原始訓練集{(xi，yi)，xi∈X，yi∈Y，i=1，2，…，N}中產生K個決策樹{h(X，θk)，k=1，2，…，K}，其中N為樣本容量，θk為由Bagging思想和特征子空間思想產生的隨機變量序列。用它們構成集成學習模型，以投票法產生最終分類決策，即

(15)

式中，H(x)為組合分類模型；h(X，θi)為單個決策樹模型；Y為目標變量；I為示性函數。

1.2.5 支持向量回歸模型

支持向量回歸(SVR)是一種用于預測連續數值型目標變量的機器學習模型。它基于支持向量機(SVM)的原理，在高維特征空間中找到一個超平面，最小化預測值與實際值之間的誤差。SVR目標函數(誤差最小化)為

(16)

式中，ω為超平面的權重向量；x為輸入特征向量；y為實際值；f(x)為預測值；ε為控制模型容忍誤差的參數；C為正則化參數。

1.3 超參數優化

超參數的選擇與估計是機器學習模型在實際應用中的關鍵問題，模型泛化性能的優劣依賴于對超參數的合理選擇[11]。超參數優化是一種通過尋找最優超參數值來優化模型性能的方法，目標是在超參數空間中找到使目標函數最小化或最大化的超參數組合。

常見的超參數優化方法有網格搜索、隨機搜索、貝葉斯優化等。其中，隨機搜索要求樣本點集足夠大，貝葉斯優化算法容易陷入局部最優值。而網格搜索在采用較大的搜索范圍及較小的補償時，有更高概率獲得全局最優值[12]。本文主要采用GridSearchCV方法進行超參數尋優。該方法是一種基于網格搜索的交叉驗證方法，將超參數空間劃分為多個網格，然后對于每個網格對應的超參數組合進行模型訓練和評估，選擇最優結果[13]。

1.4 多模型集成預測

1.4.1K折交叉驗證

K折交叉驗證方法是一種用于提高分類效果的統計分析方法。將樣本量大小為N的數據集劃分為大小近乎相同但互不交叉的K份數據集，并以不同方式組合為訓練集與測試集。“交叉”意為某種組合方式中訓練集的樣本，在另一組合方式中可能成為測試集，重復使用數據。其中訓練集用于訓練模型，測試集用于評估模型性能。該方法可使兩層使用的訓練數據不同，在一定程度上防止過擬合。

1.4.2 Stacking融合算法

Stacking融合算法被廣泛應用于多種學習器的組合。該算法分為2層，先訓練第一層中的多個基學習器，產生的結果作為輸入訓練第二層中的元學習器，得到最終輸出。基學習器的選擇對算法的預測精度有重要影響，學習能力較強的基學習器可提升算法的預測精度；同時，基學習器間的異質性可彌補不同基學習器的短板，進一步提升模型的泛化能力。元學習器的訓練集由基學習器的輸出產生，直接使用基學習器來產生元學習器時，容易造成過擬合。選擇較簡單的元學習器可有效降低過擬合風險。本文采用多元線性回歸(MLR)模型對5種單一模型預測結果進行擬合，MLR模型是多元回歸分析中最基礎，最簡單的模型，適合作為Stacking融合算法中的第二層學習器。算法原理如圖3所示。

圖3 Stacking融合算法原理示意

本文采用的Stacking融合算法中，第一層為K折交叉驗證方法下的5種單一模型預測結果，將輸入數據訓練集劃分為5折對每個單一模型進行交叉驗證，驗證后組合的訓練集作為單一模型的輸入項，得到5組預測結果。第二層將5組預測結果輸入多元線性回歸模型，分析每組輸入項的權重分配，得到最終預測結果并與單一模型結果進行比較。算法結構如圖4所示。

圖4 雙層Stacking融合算法結構示意

1.5 模型評價指標

本文選取4種常用評價指標來衡量模型預測精度，分別為平均絕對誤差(Mean absolute error，MAE)、均方誤差(Mean square error，MSE)、平均絕對百分比誤差(Mean absolute percentage error，MAPE)和確定性系數(R-square，R2)。各項指標計算公式如下

(17)

(18)

(19)

(20)

式中，MAE、MSE、MAPE指標值越小，R2指標越接近1，說明模型擬合效果好，預測結果誤差小；反之，模型擬合效果差，預測結果誤差大。

2 實例分析

2.1 研究區域概況

龍羊峽位于青海省共和縣與貴南縣交界處的黃河上游干流河段，龍羊峽水電站是黃河上游第一座大型梯級電站，水庫控制流域面積13.1萬km2，調節庫容193.6億m3，總庫容247億m3。庫區海拔高程2 600～3 000 m，多年平均流量650 m3/s。

2.2 因子相關性分析

為分析各項大氣環流指數與龍羊峽水庫月入庫徑流在時間尺度上的變化，確定預報因子的前滯期，將龍羊峽水庫月入庫徑流與前1～4個月的88項大氣環流指數進行相關性分析，結果如圖5所示。龍羊峽水庫月入庫徑流與前1個月的大氣環流指數相關性最大，前滯期為2個月、3個月時，相關性明顯隨滯時增加而減小，前滯期為4個月時部分環流指標相關性由正相關變為負相關，但大多數相關系數絕對值仍明顯小于前滯期為1個月時對應的相關系數。故選取大氣環流指數影響因子前滯期為1個月。

圖5 不同前滯期龍羊峽月平均徑流與大氣環流指標相關性比較

利用相關性分析確定水文要素因子的最佳前滯期，選取前滯期分別為1～4個月，對不同前滯期龍羊峽水庫月入庫徑流量與不同前滯期西寧站月降水量作相關性分析，如表1所示。由表1可知，龍羊峽月徑流量與西寧站降水量在前滯期為1個月時相關性最大，隨著前滯期增長，相關系數逐步減小，前滯期4個月時由正相關轉為負相關，但相關系數絕對值仍小于前滯期1個月時的相關系數。故選取水文影響因子的前滯期為1個月。將各項影響因子的月數據前滯1個月，與龍羊峽水庫當月入庫徑流量做相關性分析，選擇相關系數大于0.6的影響因子作為篩選結果。

表1 不同前滯期龍羊峽月平均徑流與西寧站降水量相關性比較

2.3 因子綜合篩選

對使用相關性分析篩選出的36項影響因子進行主成分分析，結果如表2所示。以不同的累計貢獻率為界限來確定輸入項：以80%的累計貢獻率為下限，提取主成分4個(CA1～CA4)，累計貢獻率達81.85%；以85%的累計貢獻率為下限，提取主成分6個(CA1～CA6)，累計貢獻率達87.10% ；以90%的累計貢獻率為下限，提取主成分8個(CA1～CA8)，累計貢獻率達90.46%；以95%的累計貢獻率為下限，提取主成分13個(CA1～CA13)，累計貢獻率達95.79%。根據主成分的荷載計算分別得到4組主成分序列(CA4、CA6、CA8、CA13)作為中長期徑流預報模型的輸入因子。

表2 因子綜合篩選結果 %

2.4 單一模型預測結果

以默認參數建立5個單一模型，劃分輸入因子前495項為訓練集，后200項為測試集，使用GridsearchCV方法進行超參數尋優，得到最優參數下的不同主成分組合的預測結果，對其精度指標進行比較，以最優主成分組合得到的結果作為單一模型的最優預測結果。每個模型的預測結果對比如圖6所示，精度指標如表3所示。

表3 單一模型預測精度指標(K折交叉驗證前)

圖6 單一模型預測值與實測值對比

超參數優化前后預測結果對比(以MAPE指標為例)如表4所示。

表4 超參數優化前后單一模型預測結果精度對比 %

由表3和表4可知，不同的單一模型對龍羊峽徑流預報效果具有一定的差異性；超參數尋優方法對模型精度有一定提升，但是僅局限于模型參數的優化，模型本身存在的缺點導致的精度誤差無法避免，預測結果有待進一步優化。

2.5 多模型集成預測

第一層Stacking將5個單一模型分別采用5折交叉驗證進行預測，輸入因子使用最優主成分組合，模型參數采用超參數優化后的最優參數，預測結果對比如圖7所示，精度指標如表5所示。

表5 單一模型預測精度指標(K折交叉驗證后)

圖7 K折交叉驗證后單一模型預測結果與實際值的對比

對比表4和表5可知，K折交叉驗證后MAE、MSE、MAPE指標值減小，R2指標提高，說明K折交叉驗證后一定程度上降低了過擬合，有效提高了預測精度。

第二層Stacking以單一模型的預測結果作為MLR模型的輸入因子，通過超參數尋優得到單一模型的最優權重，得到最終預測結果對比，如圖8所示，精度指標與單一模型K折交叉驗證后對比如表6所示。

圖8 Stacking融合算法預測結果與實際值的對比

由表6可知，基于多模型融合的集成預測方法的MAE、MSE、MAPE指標值最小，R2指標最接近1，說明模型擬合效果好，預測結果誤差小，預測效果相較于單模型有較大的精度提升。

3 結 論

本文提出了一種基于多模型融合的集成預測方法，通過結合5個異質模型，建立了龍羊峽水庫中長期徑流預測模型。同時，采用超參數尋優方法對各模型參數進行了優化，將得到的預測結果與單一模型的預測結果進行比較，得出以下結論：

(1)采用超參數優化方法確定模型的最優參數，可以有效提高融合后的集成模型的預測精度。

(2)本文使用了Stacking融合算法，并采用了K折交叉驗證策略，有效避免了模型訓練中出現的過擬合現象，進一步提升了預測精度。

(3)基于Stacking融合算法，將ARMA、BP、LSTM、RF和SVR等5個異質模型進行了融合。與單一模型相比，集成預測方法取得了更高的精度，R2值從0.71提升至0.82。