基于EEMD-VMD-SSA-KELM模型的汛期日徑流預測

吳小濤，袁曉輝，袁艷斌，毛雅茜，肖加清

（1. 黃岡師范學院數(shù)學與統(tǒng)計學院，湖北 黃岡 438000； 2. 華中科技大學土木與水利工程學院，湖北 武漢 430074； 3. 華中科技大學數(shù)字流域科學與技術湖北省重點實驗室，湖北 武漢 430074； 4. 武漢理工大學資源與環(huán)境工程學院，湖北 武漢 430070）

0 引 言

徑流預測可以為水資源的開發(fā)和管理、防洪、水旱災害防治及社會可持續(xù)發(fā)展等提供決策支撐［1］。由于徑流受到下墊面、水文氣象和氣候等因素的影響，呈現(xiàn)出非線性、非平穩(wěn)等復雜特性，傳統(tǒng)的物理機制水文模型的預測精度往往不高。基于數(shù)據(jù)驅動的機器學習模型如神經(jīng)網(wǎng)絡模型［2，3］、最近鄰法［4］、支持向量機（Support Vector Machine，SVM）模型［5］等在充分挖掘歷史徑流數(shù)據(jù)內在規(guī)律的基礎上，忽略復雜的下墊面情況、氣象和氣候等因素，建立未來徑流與歷史徑流的映射關系，具有操作簡單、預測精度高等優(yōu)勢，是目前主要的研究方向。蘇輝東［6］等對同一組數(shù)據(jù)分別建立了分布式水文模型、BP神經(jīng)網(wǎng)絡模型和SVM 模型，通過模擬實驗發(fā)現(xiàn)，SVM 模型的預測精度最高。然而SVM 模型的預測精度依賴于核函數(shù)中的參數(shù)以及懲罰參數(shù)的選取，并且其計算復雜度隨樣本規(guī)模的增加而變大。Huang［7］提出了核極限學習機（Kernel Extreme Learning Machine，KELM）這種單隱層前饋神經(jīng)網(wǎng)絡，該方法具有結構簡單、運算速度快及泛化能力強等優(yōu)點。Liu等［8］證明了當樣本規(guī)模較大時，KELM 模型與SVM 模型的預測精度相差不大，但是KELM 模型的計算復雜度遠遠小于SVM模型。但是KELM 模型的預測精度依賴于核參數(shù)以及懲罰系數(shù)的選取，涂異等［9］提出用粒子群優(yōu)化（Particle Swarm Optimization，PSO）算法來優(yōu)化這兩個參數(shù)，但是PSO 的尋優(yōu)速度和優(yōu)化精度不高。麻雀搜索算法（Sparrow Search Algorithm，SSA）是薛建凱等［10］在2020年提出的一種新的智能優(yōu)化算法，與傳統(tǒng)的PSO、蟻群算法等相比，SSA具有簡單易實現(xiàn)，搜索精度高，收斂速度快，穩(wěn)定性好，魯棒性強等優(yōu)點［11］。本文采用SSA 優(yōu)化KELM 模型的核參數(shù)和懲罰系數(shù)。

利用數(shù)據(jù)分解算法如集合經(jīng)驗模態(tài)分解（Ensemble Empirical Mode Decomposition，EEMD）算法、變分模態(tài)分解（Variational Modal Decomposition，VMD）算法等先對徑流序列分解，可以有效提取隱含在徑流序列中的有用信息，降低徑流序列的非平穩(wěn)性。呂晗芳［5］、孫望良［12］等建立了基于VMD算法的徑流組合預測模型，李新華［13］等建立了基于小波包分解算法的徑流組合預測模型，實驗結果表明，組合了數(shù)據(jù)分解算法的模型較單一模型的預測精度高。文獻［14］設計了基于EEMD 算法和VMD算法各自特點的二階數(shù)據(jù)分解算法，稱為EEMD-VMD 算法。EEMD-VMD 算法先利用EEMD 算法分解徑流序列，再利用VMD算法進一步分解頻率最大的分量，降低該分量及徑流序列的不穩(wěn)定性。EEMD-VMD 算法較EEMD 算法、VMD 算法具有更優(yōu)的分解性能。

為了提高徑流預測精度，本文綜合考慮EEMD、VMD 處理非平穩(wěn)徑流序列的性能以及KELM 運算速度快、泛化能力強等優(yōu)點，設計了一種基于EEMD算法、VMD算法和SSA優(yōu)化KELM的組合徑流預測模型，簡稱EEMD-VMD-SSA-KELM 模型。即首先對原始徑流序列利用EEMD-VMD 算法分解；接著，對分解得到的每個分量分別建立SSA 優(yōu)化的KELM 模型進行預測；最后，重構分量的預測結果得到原始徑流序列的預測結果。本文將該模型應用于湖北宜昌寸灘水文站的汛期日徑流預測，并與傳統(tǒng)的BP 神經(jīng)網(wǎng)絡模型、最小二乘支持向量機（Least Squares Support Vector Machine，LSSVM）模型等進行比較，采用誤差指標進行評價分析，驗證該模型的適用性。

1 研究方法及模型建立

1.1 EEMD-VMD算法

EEMD算法依據(jù)序列的極值能有效分解出反映序列總體變化情況的趨勢分量、反映序列周期性等特征的細節(jié)分量和反映噪聲、干擾特征的隨機分量。趨勢分量和細節(jié)分量的波動頻率較隨機分量小很多，采用預測模型對其預測時，預測精度往往較高。但是隨機分量由于其波動頻率最大且規(guī)律性差，采用預測模型對其預測時，預測精度非常低，影響了序列的整體預測精度。由于VMD 算法依據(jù)序列的中心頻率能將序列分解為多個不同頻率的分量，故本文利用VMD算法對波動頻率最大的隨機分量再次分解，得到若干個頻率不同、較隨機分量更穩(wěn)定的分量。EEMD-VMD 算法能有效降低序列的不穩(wěn)定性，減小隨機分量對序列的整體預測精度造成的影響［14］。利用EEMD 算法分解序列得到的分量個數(shù)由序列的長度確定，而利用VMD算法分解序列得到的分量個數(shù)K須要先確定，本文采用文獻［15］中的最優(yōu)變分模態(tài)分解算法來確定K的值。

1.2 SSA優(yōu)化KELM 模型

SSA 是一種模擬麻雀種群在覓食、反捕食過程中的搜索行為的智能優(yōu)化算法。種群中適應值較高的麻雀被稱為發(fā)現(xiàn)者，負責搜索食物并提供食物所在的區(qū)域或方向；其余的麻雀被稱為加入者，通過發(fā)現(xiàn)者提供的信息覓食。適應值最差的麻雀由于得不到食物會飛出食物所在區(qū)域，另尋他處覓食。麻雀之間會互相競爭，因而其身份也會動態(tài)變化，但是麻雀總數(shù)和比例保持不變。另外，種群中的麻雀一旦發(fā)現(xiàn)捕食者會立刻發(fā)出警告信息，所有麻雀飛出覓食區(qū)域，發(fā)現(xiàn)捕食者的麻雀被稱為偵查者。SSA 模擬發(fā)現(xiàn)者的行為實現(xiàn)全局搜索，模擬加入者的行為實現(xiàn)局部探索，模擬適應值最差的麻雀以及麻雀的反捕食行為擴大搜索區(qū)域。發(fā)現(xiàn)者的位置更新策略［10］如式（1）所示。

式中：t為當前迭代數(shù)；Xt+1i為第i只麻雀在第t+1次迭代時的適應值，α∈(0，1)為隨機數(shù)；itermax為最大迭代次數(shù)；R2和ST分別表示預警值和安全值；Q為服從標準正態(tài)分布的隨機數(shù)；L為全一行向量。加入者的位置更新策略［10］如式（2）所示。

式中：Xtw為麻雀的最差位置；Xtp為發(fā)現(xiàn)者的最佳位置；A+為元素隨機取-1 或1 的行向量。偵查者的位置更新策略［10］如式（3）所示。

式中：Xtb為全局最佳位置；γ為步長控制參數(shù)；K∈(0，1)為隨機數(shù)；ε為保證分母不等于0 的常數(shù)；fi為當前麻雀的適應值；fb、fw分別為當前全局最優(yōu)適應值和最差適應值。

KELM 是一種改進的ELM 模型。對于含有n個樣本的數(shù)據(jù)集{(x1，t1)，(x2，t2)，…，(xn，tn)}，其中xi=[xi1，xi2，…，xid]∈Rd是輸入數(shù)據(jù)，i= 1，2，…，n，ti=[ti1，ti2，…，xik]∈Rk是輸出數(shù)據(jù)。ELM 模型假設隱含層有m個節(jié)點，其輸出函數(shù)為h(x) =[h1(x)，h2(x)，…，hm(x)]=H，隱含層與輸出層節(jié)點之間的連接權重為β=[β1，β2，…，βm]，則其輸出可表示為：

為了最小化訓練誤差ε和連接權重β，ELM 模型利用結構風險最小原理，構造二次規(guī)劃問題如下：

式中：C為懲罰系數(shù)；εi為第i個誤差變量，ε=[ε1，ε2，…，εn]。

進一步通過拉格朗日乘子αi，將式（5）轉化為如下無約束優(yōu)化問題：

其中，I是單位矩陣，y=[y1，y2，…，yn] 是期望輸出向量。

ELM中的h(x)為隨機映射，影響了ELM模型的穩(wěn)定性和預測精度，Huang提出采用核函數(shù)代替h(x)，將訓練樣本映射到高維空間，建立了基于核的ELM 模型即KELM 模型［7］。KELM 模型的核矩陣定義為Ω=HHT，其元素定義為Ω(i，j) =h(xi)h(xj)=K(xi，xj)，其中K(xi，xj) 是核函數(shù)。KELM 模型的輸出可表示為：

其中，核函數(shù)K(xi，xj)選擇高斯核函數(shù)，定義如下：

其中，σ為核參數(shù)，其取值直接影響到KELM 的泛化能力。另外，懲罰系數(shù)C會影響KELM 模型的預測精度，C的值偏小會產(chǎn)生較大的訓練誤差，偏大會產(chǎn)生過擬合。為了提高KELM 模型的泛化能力和預測精度，本文采用SSA 來優(yōu)化σ和C。選擇平均絕對百分比誤差作為SSA的適應值計算函數(shù)，即：

式中：Num為樣本總數(shù)；t*和y*分別為徑流的實測值和預測值。

SSA優(yōu)化σ和C的步驟為：

步驟1：初始化參數(shù)，包括最大迭代次數(shù)、種群規(guī)模，發(fā)現(xiàn)者數(shù)量、加入者數(shù)量和偵查者數(shù)量、預警值、麻雀的初始位置（即σ、C）及其優(yōu)化區(qū)間等。

步驟2：根據(jù)初始參數(shù)由式（4） ～（8）及式（10）計算得到麻雀個體的初始適應值，保存適應值最優(yōu)、最差的麻雀個體的位置Xb和Xw。

步驟3：由式（1）更新發(fā)現(xiàn)者位置，由式（2）更新加入者位置，由式（3）更新偵查者位置。

步驟4：計算更新位置后所有麻雀的適應值，更新Xb和Xw。

步驟5：判斷是否滿足迭代終止條件（最大迭代次數(shù)或誤差）。如果不滿足則返回到步驟3，滿足則尋優(yōu)結束，將適應值最優(yōu)的麻雀的位置σ和C作為輸出。

1.3 EEMD-VMD-SSA-KELM 預測模型

為了提高徑流的預測精度，本文設計了基于“分解-建模預測-重構”思想的徑流預測框架，如圖1 所示，具體預測步驟如下：

圖1 EEMD-VMD-SSA-KELM 預測模型框架圖Fig.1 Framework diagram of EEMD-VMD-SSA-KELM prediction model

步驟1：數(shù)據(jù)預處理。采用拉依達準則處理原始徑流序列中的異常值。

步驟2：EEMD分解。利用EEMD算法對預處理后的徑流序列進行分解，得到n個分量IMF1～ IMFn。

步驟3：VMD分解。利用VMD算法對頻率最大的分量IMF1進行分解，得到K個分量IMF11～ IMF1K以及一個殘差分量IMF1R。

步驟4：數(shù)據(jù)準備。將每個分量歸一化后劃分為訓練數(shù)據(jù)集和測試數(shù)據(jù)集。

步驟5：建模預測。基于每個分量的訓練數(shù)據(jù)集建立SSA優(yōu)化的KELM模型，并對測試數(shù)據(jù)集進行預測，將預測結果作反歸一化處理。

步驟6：重構。累加步驟5 中所有分量的預測結果，得到原始徑流序列的預測結果。

2 實例分析

2.1 數(shù)據(jù)來源

選取湖北宜昌寸灘水文站2006-2012年間的汛期日徑流序列進行實驗。寸灘水文站是國家基本水文站、中央報汛站，其防汛測報地位和功能十分突出。因為2006-2012 年間每年的1-5月以及10-12月的日徑流較小且波動不大，而6-9月處于汛期，其日徑流大且波動也較大，所以本文僅對6-9月的徑流進行預測。根據(jù)拉依達準則對2006-2012 年間每年的6-9 月的日徑流數(shù)據(jù)進行處理后得到854 個數(shù)據(jù)，汛期日徑流隨時間的變化情況如圖2所示。

圖2 汛期日徑流序列曲線圖Fig.2 Daily runoff sequence curve in flood season

2.2 數(shù)據(jù)分解

對上述徑流序列利用EEMD 算法進行分解，得到9 個分量IMF1～IMF9，如圖3所示。

圖3 EEMD 算法分解結果Fig.3 Decomposition results using EEMD algorithm

圖3中，IMF9是趨勢分量，波動頻率最小，IMF2～IMF8是細節(jié)分量，IMF1是隨機分量，波動頻率最大。利用VMD 算法對IMF1進一步分解，當模態(tài)個數(shù)K=9 時，殘差分量與IMF1的相關系數(shù)RK等于0.001，小于給定的閾值0.01，從而K的值取9。9 個分量和1 個殘差分量分別記為IMF11～IMF19、IMF1R，分解結果如圖4所示。

圖4 VMD算法分解結果Fig.4 Decomposition results using VMD algorithm

從圖4 可以看出，IMF11的波動頻率最小，IMF12～IMF19的波動頻率逐漸增加，但是這些分量的徑流值幾乎對稱分布在零的上下兩側，穩(wěn)定性較IMF1顯著增強。盡管IMF1R的波動頻率最大，但是其變化范圍非常小，其預測結果對序列的整體預測精度影響很小。

2.3 預測模型的輸入與輸出

對2012 年6-9 月每日的徑流進行預測，選擇2006-2011 年間每年6-9月的日徑流數(shù)據(jù)（共732個，編號為1～732）作為訓練數(shù)據(jù)，2012 年6-9 月的日徑流數(shù)據(jù)（共122 個，編號為733～854）作為測試數(shù)據(jù)。徑流序列或各個分量的預測模型的輸入維度根據(jù)徑流序列或分量的偏自相關系數(shù)（PACF）的延遲時間數(shù)來確定［16］。例如，徑流序列的PACF結果如圖5 所示，從圖5 可以看出，在所有這些延遲中，前3個延遲超過了95%置信區(qū)間對應的閾值，因此，對徑流序列直接進行預測的預測模型的輸入維度取3。

圖5 徑流序列的PACFFig.5 PACF of Runoff Series

各個分量的預測模型的輸入維度如表1所示。

表1 各個分量的預測模型的輸入維度Tab.1 Input dimension of prediction model of each component

2.4 預測結果分析

將本文提出的EEMD-VMD-SSA-KELM 模型與SSA-BP 模型、SSA-LSSVM 模型、SSA-KELM 模型、EEMD-SSA-KELM 模型和EEMD-VMD-PSO-KELM 模型作比較。前3 種模型中，采用SSA 分別優(yōu)化BP 神經(jīng)網(wǎng)絡模型的神經(jīng)元個數(shù)、LSSVM 模型的參數(shù)以及KELM 模型的參數(shù)。EEMD-SSA-KELM 模型和EEMD-VMD-PSO-KELM 模型與提出的模型類似，區(qū)別在于EEMD-SSA-KELM 模型只組合了EEMD 算法，提出的模型組合了EEMD-VMD 算法；EEMD-VMD-PSO-KELM 模型采用PSO優(yōu)化KELM 模型的參數(shù)，提出的模型采用SSA 優(yōu)化KELM 模型的參數(shù)。SSA-BP 模型、SSA-LSSVM 模型和SSA-KELM 模型的輸入維度均為3，EEMD-SSA-KELM 模型、EEMD-VMD-PSOKELM 模型以及EEMD-VMD-SSA-KELM 模型中的分量的預測模型的輸入維度如表1。6 種預測模型的預測曲線圖如圖6所示。

圖6 預測曲線圖Fig.6 Prediction curves

從圖6 可以看出，除了矩形區(qū)域對應的時間段（編號為777～792）外，6 種模型在其余時間的預測曲線整體上與實測值曲線變化趨勢類似。在所有時間點，EEMD-VMD-SSA-KELM模型和EEMD-VMD-PSO-KELM 模型的預測值與實測值無限接近，因此其預測精度最高，EEMD-SSA-KELM 模型的預測精度較高，SSA-BP 模型、SSA-LSSVM 模型和SSA-KELM 模型的預測精度最低。進一步，將矩形區(qū)域的預測曲線放大，如圖6中的小圖所示，可以看出，組合了數(shù)據(jù)分解算法的EEMD-SSAKELM 模型、EEMD-VMD-PSO-KELM 模型和EEMD-VMDSSA-KELM 模型的預測曲線較沒有組合數(shù)據(jù)分解算法的SSABP 模型、SSA-LSSVM 模型和SSA-KELM 模型的預測曲線更接近于實測值曲線。統(tǒng)計6 種預測模型的誤差指標MAPE、MRE、RMSE和R2如表2所示。

表2 預測誤差表Tab.2 Error values

由表2 可知，SSA-LSSVM 模型的4 種誤差指標值均劣于SSA-KELM 模型，說明KELM 模型的預測精度優(yōu)于LSSVM 模型。沒有組合數(shù)據(jù)分解算法的SSA-BP 模型、SSA-LSSVM 模型和SSA-KELM 模型的MAPE均高于6%，MRE均高于2%，RMSE均高于170 m3∕s，R2均低于97%，而組合了數(shù)據(jù)分解算法的3 種模型的MAPE均低于5%，MRE均低于1.5%，RMSE均低于110 m3∕s，R2均高于98%，說明組合了數(shù)據(jù)分解算法的預測模型的預測精度明顯優(yōu)于沒有組合數(shù)據(jù)分解算法的預測模型。進一步對比EEMD-SSA-KELM 模型和EEMD-VMD-PSO-KELM 模型、EEMD-VMD-SSA-KELM 模型的誤差指標值可以看出，EEMDSSA-KELM 模型的4 種誤差指標值均劣于后兩種模型，說明組合EEMD 算法的預測模型的預測精度劣于組合了EEMD-VMD算法的預測模型，即EEMD-VMD 算法較EEMD 算法具有更優(yōu)的分解性能。對比EEMD-VMD-PSO-KELM 模型和EEMDVMD-SSA-KELM 模型的4 種誤差指標值，除了MRE相同，均為0.59 % 之外，EEMD-VMD-SSA-KELM 模型的其余3 個誤差均優(yōu)于EEMD-VMD-PSO-KELM 模型，說明了SSA 較PSO 具有更優(yōu)的尋優(yōu)精度。6 種預測模型的4 種誤差指標的直方圖如圖7所示。

圖7 誤差指標直方圖Fig.7 Histogram of error index

由圖7（a）可知，對于MAPE指標，SSA-BP 模型的MAPE最高，等于6.76%，而EEMD-VMD-SSA-KELM 模型最低，只有1.78%，說明提出的模型的預測值的平均絕對百分比誤差最小。由圖7（b）可知，對于MRE指標，SSA-BP 模型、SSA-LSSVM 模型和SSA-KELM 模型的MRE均超過了2%，而EEMD-VMD-SSAKELM 模型的MRE最低，僅為0.59%，說明提出的模型的預測值和實測值之間的接近程度最高。由圖7（c）可知，對于RMSE指標，SSA-LSSVM 模型的RMSE最高，等于214.41 m3∕s，而EEMDVMD-SSA-KELM 模型的RMSE最低，僅為38.94 m3∕s，說明提出的模型的預測值與實測值之間的偏差最小。由圖7（d）可知，對于R2指標，EEMD-VMD-SSA-KELM 模型的R2達到了99.84 %，高于其余5 種模型，說明提出的模型的擬合程度最高。綜上所述，EEMD-VMD-SSA-KELM 模型的4個誤差指標MAPE、MRE、RMSE和R2均最優(yōu)，其預測精度遠優(yōu)于SSA-BP、SSA-LSSVM、SSA-KELM 和EEMD-SSA-KELM 等模型。EEMD-VMD-SSAKELM 模型能較準確的模擬復雜多頻信息的汛期日徑流變化趨勢，模型及方法可為水文預測及相關預測研究提供參考。

3 結 論

針對汛期日徑流序列非線性、非平穩(wěn)等特點導致徑流預測精度低的問題，本文提出了一種基于EEMD-VMD 算法和SSA優(yōu)化KELM 的組合徑流預測模型（EEMD-VMD-SSA-KELM），并將該模型應用于湖北宜昌寸灘水文站的汛期日徑流預測，得到如下結論。

（1）在本實驗中，KELM 模型的預測精度優(yōu)于LSSVM 模型，SSA 的優(yōu)化精度優(yōu)于PSO，EEMD-VMD 算法的分解性能優(yōu)于EEMD算法。

（2）組合了數(shù)據(jù)分解算法的EEMD-SSA-KELM 模型、EEMD-VMD-PSO-KELM 模型和EEMD-VMD-SSA-KELM 模型的預測精度明顯優(yōu)于沒有組合數(shù)據(jù)分解算法的SSA-BP 模型、SSA-LSSVM 模型和SSA-KELM 模型。但由于組合了數(shù)據(jù)分解算法的預測模型要對每個分量建立模型預測，因此其時間復雜度更高。

（3）EEMD-VMD-SSA-KELM 模型的預測精度遠優(yōu)于SSABP、SSA-LSSVM、SSA-KELM、EEMD-SSA-KELM 等模型，能較準確的模擬復雜多頻信息的汛期日徑流的變化趨勢。

本文的不足之處在于只對汛期日徑流序列進行了單步預測，在后續(xù)研究中，可以用提出的模型在單步預測的基礎上，進一步進行多步預測。另外，對于提出模型復雜度高的問題，可以考慮引入樣本熵，將復雜性相近的分量重組，減少分量個數(shù)，降低模型的時間復雜度。