基于變分模態分解和引力搜索算法的徑流預測模型研究

李曉春， 包苑村， 羅軍剛， 左崗崗

(1.陜西省江河水庫工作中心， 陜西 西安 710018；2.西安理工大學 西北旱區生態水利國家重點實驗室，陜西 西安 710048)

月徑流時間序列受到多種因素的影響和制約，當前流行的數據驅動模型能很好地捕捉徑流時間序列的非線性關系，但由于月徑流屬于小樣本，使得數據驅動模型容易陷入過擬合以及局部最小值等情形[1]。支持向量回歸機(SVR)采用結構風險最小化準則設計學習機器，折衷考慮經驗風險和置信范圍，在小樣本的情況下具有更好的泛化能力[2]。與此同時，徑流時間序列的高度非線性和不穩定性，決定了需要將SVR與其它模型方法相結合，才能更好地提升預測精度。黃巧玲等[3]將耦合離散小波變換(DWT)與支持向量回歸機(SVR)結合，建立了月徑流預測的小波支持向量機耦合模型(WSVR)，并應用于涇河張家山站的月徑流預測中。結果表明，WSVR模型可以有效提高單一支持向量機模型的預測精度。周有榮等[4]利用同熱傳遞搜索(SHTS)算法，對混合核支持向量機(SVM)關鍵參數和混合權重系數進行優化，結果表明，SHTS算法尋優精度高于TLBO、GWO優化算法。李祥蓉[5]利用靜電放電算法(ESDA)優化混合核SVM關鍵參數和混合權重系數，研究結果表明，混合核ESDA-SVM模型在預測精度、泛化能力等方面均優于對比模型，具有較好的實際應用價值。王遷等[6]將粒子群優化算法(PSO)引入到SVR模型中，建立了PSO-SVR模型，實現了對SVR的RBF核函數的三個參數的自動優選。實驗表明，PSO-SVR模型較ANN模型穩定性更強、可信度更高，且具有更好的應用價值。梁浩等[7]在優選多元線性回歸(MLR)、人工神經網絡(ANN)和支持向量機(SVM)單一預報模型的基礎上，分別基于經驗模態分解(EMD)、集合經驗模態分解(EEMD)和小波分解(WD)構建了多種混合模型。結果表明，混合預測模型的預測精度均高于單一模型。

上述研究成果成功將多種優化算法與模型組合，在實例應用中取得了不錯的效果。但大多數研究方法在實驗開始時就將整個徑流時間序列分解，實際帶入了徑流時間序列的未來信息[8]。同時，許多優化算法易陷入局部最優[9]。因此，本文提出基于變分模態分解(VMD)、引力搜索(GSA)與支持向量回歸機(SVR)的組合模型。VMD分別對訓練集數據和測試集數據進行分解，GSA對SVR的參數進行全局尋優。將本模型應用于渭河流域臨潼站與咸陽站的月徑流預測中。實例研究表明，該模型在有效提高預測精度的同時，也更加符合實際預測過程。

1 研究方法

1.1 Mann-Kendall(M-K)檢驗

在月徑流的預測中，通常訓練集與測試集的比例為7∶3左右。由于月徑流屬于小樣本，按比例分配后，測試集的數據量過少，缺乏實際研究意義，因此引入Mann-Kendall(M-K)檢驗。

M-K檢驗是一種檢驗水文時間序列的趨勢以及突變點的非參數統計檢驗[10]。給定月徑流時間序列變量x={x1,x2,…,xn}，則統計量Sk為：

(1)

其中：

(2)

M-K檢驗假設x有獨立的觀測值，如果沒有趨勢存在，這些觀測值分布相同。在此假設下，Sk的均值和方差分別為：

E(Sk)=0

(3)

(4)

x的檢驗統計量定義為：

(5)

在顯著性水平1-α(Z1-α)下，通過比較UFn與標準正態變量來檢測x的趨勢，α(0<α<0.5)是M-K檢驗錯誤拒絕原假設H0的容忍概率。當|UFn|≥Z1-α時，拒絕H0；否則，接受H0。在拒絕H0的條件下，如果|UFn|大于(小于)零，x隨時間呈單調增加(減少)的趨勢。

1.2 變分模態分解(VMD)

鑒于月徑流時間序列的高度不穩定性，引入電力系統的信號分解模式[11]，將月徑流序列分解為多個分量，以便模型更好地學習月徑流的變化規律。

VMD是一種新的信號分解模式，對原始信號采用非遞歸和變分模式分解的方法[12]，將輸入信號分解成若干個子模式，即周期性的分量(IMF)和一個殘差(R)。VMD可以手動設置分解的模態數K，對噪聲具有較好的魯棒性，能顯著降低計算的復雜度。因此，VMD可以對訓練集和測試集分別進行分解，以避免模型在訓練過程中混入測試集的信息，同時也能保證訓練集與測試集的特征維度一致。

VMD的核心是一個受限的變分問題，將非平穩信號f分解為K個有限帶寬的模態分量，為了保證每個模態分量的估計帶寬之和最小，須使所有模態之和與原始信號相等。Dragomiretskiy與Zosso在2014年提出了該受限變分問題[13]：

(6)

式中：uk是模態函數的集合；ωk是第k個模態的中心頻率；K是模態數；δ(t)是狄拉克分布；?是卷積運算；uk(t)是各子序列的模態函數； e-jωkt是復平面上模態函數中心頻率的向量描述。

1.3支持向量回歸機(SVR)原理及其參數

支持向量回歸機(SVR)首先通過核函數[14]將低維的非線性回歸問題映射至高維的空間，在高維的空間計算回歸函數，基于結構風險最小化原則[15]，有效避免因數據量不足而引起的預測精度過低、泛化能力差等問題。具體的計算過程為：

給定訓練樣本數據S={(xi,yi)|xi∈R,yi∈R,i=1,2,…,m}，其中xi為模型輸入，yi為模型輸出，m為樣本數量。則SVR的回歸函數為：

y=wT·φ(x)+b

(7)

式中：wT為權重向量；b為偏置向量；φ(x)用于將輸入向量映射至高維空間。

假設我們能容忍f(x)與yi之間最多有ε的偏差，即僅當f(x)與yi之間差值的絕對值大于ε時才計算損失，于是SVR問題可以轉化為：

(8)

其中：

f(x)-yi≤ε+ξi

運用拉格朗日乘子求解式(8)得：

(9)

由Karush-Kuhn-Tucker[16]條件解得的SVR模型可以表示為：

(10)

式中：K(xi,x)=φ(xi)Tφ(xj)為核函數。

支持向量回歸機中常用的核函數包括線性核函數、多項式核函數和高斯徑向基核函數[17](RBF)。針對徑流預測問題，RBF核函數具有超參數少、映射維度高與決策邊界多樣的優點，可以更好地適應非線性的徑流預測問題。

RBF核函數的主要超參數有懲罰系數C、不敏感損失函數ε。除此之外，RBF核函數還有一個獨有參數gamma，隱含地決定了數據映射到新的特征空間后的分布，gamma值越大，支持向量越少，gamma值越小，支持向量越多，支持向量的個數會影響訓練與預測的速度。

1.4 引力搜索算法(GSA)

引力搜索算法(GSA)是一種新的啟發式搜索算法[18]，它受萬有引力定律的啟發，適用于模型的參數尋優，具有較好的全局搜索能力。GSA可以描述為一個n維的空間中有s個粒子，則第i個粒子的位置定義為：

Xi=(xi1,…,xid,…,xin)

(11)

式中：xid表示第i個粒子在第d維的位置；n為搜索空間的維度。

在開始搜索前，所有粒子的位置都是隨機的。在某一時刻t，粒子i與j之間的引力為：

(12)

式中：Mpi(t)為施力物體j的慣性質量；Mai(t)為受力物體i的慣性質量；G(t)為引力常數，其值隨時間的變化而變化，如式(13)所示(通常設置G0為100，α為10)；Rij(t)為i與j之間的歐式距離，如式(14)所示；ε為一個小的常數，用來防止i與j之間的歐式距離為0。

(13)

(14)

粒子i在t時刻受到的其它粒子的引力為：

(15)

式中：randj為[0,1]范圍內的隨機數。

在第d維，粒子i的加速度為：

(16)

式中：Mii(t)為粒子i的慣性質量。

假設重力質量和慣性質量相等，則重力質量和慣性質量可更新為：

Mai=Mpi=Mii=Mi,i=1,2,…,n

(17)

(18)

(19)

式中：fiti(t)為粒子i在t時刻的適應度值。worst(t)與best(t)在求解最小值問題時定義為：

(20)

(21)

在求解最大值問題時，則定義為：

(22)

(23)

將粒子i的加速度加到當前的速度上，就可以得到現在的粒子速度和位置：

(24)

(25)

每次迭代后，同步更新粒子的位置和速度，直至達到最大迭代次數或達到要求的精度。

1.5 評價指標

為了更好地反映模型的預測效果。本實驗選取均方誤差(MSE)和納什系數(NSE)對測試集的預測結果進行評價。

(26)

(27)

式中：yi為i時刻的預測值；y0為i時刻的實測值；y為實測值的均值；n為測試集樣本數量。

1.6 VMD-GSA-SVR預測模型

VMD-GSA-SVR預測模型是一種新型的集成預測模型，其具體的預測流程如圖1所示。

圖1 VMD-GSA-SVR模型預測流程圖Fig.1 VMD-GSA-SVR model prediction flow char

1) 對實測月徑流數據進行M-K突變點檢測，突變點前的數據作為訓練集，突變點后的數據作為測試集。

2) 對訓練集數據進行VMD分解，確定模態數K，再將測試集數據分解成同樣的K個模態分量。

3) 將訓練集輸入SVR模型中，運用GSA對SVR的三個參數C、gamma、ε進行全局尋優，用MSE作為GSA的適應度值公式。

4) 將測試集代入經過參數優化的GSA-SVR模型進行預測，并用評價指標進行評價。

2 實例分析

2.1 資料來源與實驗環境搭建

咸陽水文站位于陜西省咸陽市東關，始建于1931年6月10日。1957年6月20日，其基本水尺斷面遷至秦皇路咸陽橋上游一側(上游2.6 km處)，更名為咸陽(二)站。該站東經108°42′，北緯34°19′，控制流域面積46 827 km2，距河源606.9 km，距河口211.1 km。臨潼水文站距咸陽站下游53.7 km，位于臨潼區行者鄉，東經109°12′，北緯34°26′，集水面積97 299 km2。本文選取咸陽站1956—2010年660個月的實測月徑流數據，選取臨潼站1960—2006年564個月的實測月徑流數據，如圖2和圖3所示。

圖2 咸陽站月徑流序列Fig.2 Monthly runoff series of Xianyang Station

圖3 臨潼站月徑流序列Fig.3 Monthly runoff series of Lintong Station

本實驗編程使用python 3.7，其中VMD分解使用第三方庫vmdpy；SVR使用sklearn庫；GSA與對比優化算法使用第三方庫optimal。

2.2 訓練集與測試集的劃分

圖4 咸陽站月徑流M-K突變點檢測結果Fig.4 Monthly runoff M-K mutation point detection results at Xianyang Station

圖5 臨潼站月徑流M-K突變點檢測結果Fig.5 Monthly runoff M-K mutation point detection results at Lintong Station

由圖4與圖5可知，咸陽站的突變點在1985年左右，因此選擇1985年前348個月的實測月徑流作為訓練集，1985年后312個月的數據作為測試集；臨潼站的突變點在1989年左右，因此選擇1989年前300個月的實測月徑流作為訓練集，1989年后264個月的數據作為測試集。

2.3 月徑流的VMD分解

多數信號分解方法不能手動選擇模態數K，因此，很多研究在實驗開始時就將整個徑流序列進行分解，再劃分訓練集與測試集，這樣做實際上將測試集的未知信息帶入到了模型的訓練過程中。鑒于此，可先對訓練集進行分解，再選擇和訓練集相同的模態個數K，對測試集進行分解。

確定訓練集VMD的分解模態數K時，采用枚舉法，兩個站的訓練集數據在K=8時分解效果最好，且無模態混疊現象[19]。其余初始參數γ、α分別設置為0、2 000。測試集數據則使用與訓練集相同的參數進行分解。

如圖6與圖7所示，兩個站點的訓練集數據和測試集數據分別被分解為7個周期性分量(IMF)和1個殘差(R)。IMF的頻率由大到小排列，R表示序列的趨勢走向。每一個模態都表現出原始序列的特征，使模型能更加準確地學習徑流序列的周期性與規律性特征。根據月徑流的年際變化規律，將每一個模態分量滯后12個月作為模型的輸入，原序列的第13個月作為輸出。在考慮滯后的情況下，咸陽站的訓練集樣本實際為336個月，測試集為300個月；臨潼站的訓練集樣本實際為288個月，測試集為252個月。

圖6 咸陽站月徑流VMD分解Fig.6 VMD decomposition of monthly runoff at Xianyang Station

圖7 臨潼站月徑流VMD分解Fig.7 VMD decomposition of monthly runoff at Lintong Station

2.4 模型預測及結果評價

為了便于比較，選擇無優化調參的SVR和PSO優化調參的SVR作為對比。SVR的核函數選擇RBF核函數，其三個參數C、gamma、ε的搜索區間分別為[0.01,100]、[0.000 001,1]、[0.000 001,1]，迭代次數n為100次，種群規模N=40。其中GSA的初始參數G0=100，α為10；PSO的初始參數c1=1.5、c2=1.7。

輸入訓練集進行模型訓練，得到優化后的SVR參數，如表1所示。

表1 優化后的SVR參數結果Tab.1 Optimized SVR parameter results

圖8和圖9分別為咸陽站、臨潼站測試集預測結果。表2為各站預測評價結果。

圖8 咸陽站測試集預測結果Fig.8 Prediction results of test set at Xianyang Station

表2 各站預測評價結果Tab.2 Prediction and evaluation results of each station

利用優化好的參數建立PSO-SVR、GSA-SVR模型，將測試集輸入模型得出預測結果并進行評價。由圖8和圖9可知，在未使用優化算法時，VMD-SVR模型的擬合效果較VMD-PSO-SVR模型、VMD-GSA-SVR模型明顯偏差，這證明了優化算法對SVR學習過程的重要性。相較于VMD-PSO-SVR模型, VMD-GSA-SVR模型對峰值與谷值的預測更為精確，表明GSA算法對SVR參數的尋優能力明顯優于PSO算法。引力搜索具有良好的全局搜索能力，而粒子群算法易陷入局部最優，因此，通過引力搜索優化的支持向量回歸機取得了更優的預測效果。

由表2可知，VMD-GSA-SVR模型對咸陽站月徑流預測的MSE為0.306 3，較VMD-SVR、VMD-PSO-SVR模型分別降低了84%、35%；VMD-GSA-SVR模型對咸陽站月徑流預測的NSE為0.946 2，較VMD-SVR、VMD-PSO-SVR模型分別提高了45%、2%。VMD-GSA-SVR模型對臨潼站月徑流預測的MSE為0.712 1，較VMD-SVR、VMD-PSO-SVR模型分別降低了了66%、9%；VMD-GSA-SVR模型對臨潼站月徑流預測的NSE為0.952 0，較VMD-SVR、VMD-PSO-SVR模型分別提高了61%、2%。由此可知，VMD-GSA-SVR模型在兩個站的預測結果評價中均達到最優，這說明該模型具有良好的擬合精度以及泛化能力。

3 結 論

本文針對月徑流小樣本數據，建立了VMD分解、GSA與SVR組合的VMD-GSA-SVR模型，并將其應用于渭河流域咸陽站和臨潼站的實測月徑流預測中。

1) 利用M-K突變點檢測劃分訓練集和測試集，能更好地測試模型對未知數據的擬合程度。

2) 先對訓練集進行VMD分解，再將訓練集的VMD參數代入測試集進行分解。這樣在分解過程中就不會將測試集的信息帶入模型的訓練過程，更加符合預測實際。

3) 利用GSA對SVR的三個參數進行尋優。相較于PSO-SVR模型，GSA-SVR模型在預測中取得了更低的誤差以及更高的精度，證明GSA是一種可靠有效的優化算法。

4) 將VMD-GSA-SVR模型應用于渭河流域咸陽站和臨潼站的實測月徑流預測中。相較于VMD-SVR模型與VMD-PSO-SVR模型，VMD-GSA-SVR模型均取得了最優的預測結果。該模型為渭河流域的月徑流預測提供了一條新的途徑。