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

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

0 引 言

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

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

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

1 研究方法

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

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

1.1 預(yù)報因子的篩選

1.1.1 初始預(yù)報因子集

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

1.1.2 相關(guān)性分析

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

(1)

1.1.3 主成分分析

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

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

(2)

對標準化矩陣Z求相關(guān)系數(shù)矩陣，rij為兩變量之間的相關(guān)系數(shù)，相關(guān)系數(shù)矩陣R為

R=[rij]n

(3)

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

(4)

(5)

1.2 單一預(yù)測模型

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

1.2.1 自回歸滑動平均模型

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

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為常數(shù)項；φi和θt分別為AR和MA模型的系數(shù)；εt為誤差項；μ為時間序列的均值。

1.2.2 BP神經(jīng)網(wǎng)絡(luò)模型

BP預(yù)測模型是一種基于人工神經(jīng)網(wǎng)絡(luò)的機器學(xué)習(xí)算法，用于解決回歸和分類問題。它通過反向傳播算法來訓(xùn)練神經(jīng)網(wǎng)絡(luò)，不斷調(diào)整權(quán)重和偏置使預(yù)測值與實際值之間的誤差最小化。模型結(jié)構(gòu)如圖2所示。

圖2 BP神經(jīng)網(wǎng)絡(luò)模型結(jié)構(gòu)

1.2.3 長短期記憶神經(jīng)網(wǎng)絡(luò)模型

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

遺忘門中輸入xt與狀態(tài)記憶單元St-1、中間輸出ht-1共同決定狀態(tài)記憶單元遺忘部分。輸入門中的xt分別經(jīng)過sigmoid和tanh函數(shù)變化后共同決定狀態(tài)記憶單元中保留向量。中間輸出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分別為遺忘門、輸入門、輸入節(jié)點、輸出門、中間輸出和細胞單元的狀態(tài)；Wfx、Wfh、Wix、Wih、Wgx、Wgh、Wox、Woh分別為相應(yīng)門與輸入xt和中間輸出ht-1相乘的矩陣權(quán)重；bf、bi、bg、bo分別為相應(yīng)門的偏置項；⊙為向量中元素按位相乘；σ為sigmoid函數(shù)變化；φ為tanh函數(shù)變化。

1.2.4 隨機森林模型

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

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

(15)

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

1.2.5 支持向量回歸模型

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

(16)

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

1.3 超參數(shù)優(yōu)化

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

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

1.4 多模型集成預(yù)測

1.4.1K折交叉驗證

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

1.4.2 Stacking融合算法

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

圖3 Stacking融合算法原理示意

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

圖4 雙層Stacking融合算法結(jié)構(gòu)示意

1.5 模型評價指標

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

(17)

(18)

(19)

(20)

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

2 實例分析

2.1 研究區(qū)域概況

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

2.2 因子相關(guān)性分析

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

圖5 不同前滯期龍羊峽月平均徑流與大氣環(huán)流指標相關(guān)性比較

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

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

2.3 因子綜合篩選

對使用相關(guān)性分析篩選出的36項影響因子進行主成分分析，結(jié)果如表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%。根據(jù)主成分的荷載計算分別得到4組主成分序列(CA4、CA6、CA8、CA13)作為中長期徑流預(yù)報模型的輸入因子。

表2 因子綜合篩選結(jié)果 %

2.4 單一模型預(yù)測結(jié)果

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

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

圖6 單一模型預(yù)測值與實測值對比

超參數(shù)優(yōu)化前后預(yù)測結(jié)果對比(以MAPE指標為例)如表4所示。

表4 超參數(shù)優(yōu)化前后單一模型預(yù)測結(jié)果精度對比 %

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

2.5 多模型集成預(yù)測

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

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

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

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

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

圖8 Stacking融合算法預(yù)測結(jié)果與實際值的對比

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

3 結(jié) 論

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

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

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

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