基于變分模態(tài)分解和引力搜索算法的徑流預(yù)測模型研究

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

(1.陜西省江河水庫工作中心， 陜西 西安 710018；2.西安理工大學(xué) 西北旱區(qū)生態(tài)水利國家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710048)

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

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

1 研究方法

1.1 Mann-Kendall(M-K)檢驗(yàn)

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

M-K檢驗(yàn)是一種檢驗(yàn)水文時(shí)間序列的趨勢以及突變點(diǎn)的非參數(shù)統(tǒng)計(jì)檢驗(yàn)[10]。給定月徑流時(shí)間序列變量x={x1,x2,…,xn}，則統(tǒng)計(jì)量Sk為：

(1)

其中：

(2)

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

E(Sk)=0

(3)

(4)

x的檢驗(yàn)統(tǒng)計(jì)量定義為：

(5)

在顯著性水平1-α(Z1-α)下，通過比較UFn與標(biāo)準(zhǔn)正態(tài)變量來檢測x的趨勢，α(0<α<0.5)是M-K檢驗(yàn)錯(cuò)誤拒絕原假設(shè)H0的容忍概率。當(dāng)|UFn|≥Z1-α?xí)r，拒絕H0；否則，接受H0。在拒絕H0的條件下，如果|UFn|大于(小于)零，x隨時(shí)間呈單調(diào)增加(減少)的趨勢。

1.2 變分模態(tài)分解(VMD)

鑒于月徑流時(shí)間序列的高度不穩(wěn)定性，引入電力系統(tǒng)的信號(hào)分解模式[11]，將月徑流序列分解為多個(gè)分量，以便模型更好地學(xué)習(xí)月徑流的變化規(guī)律。

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

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

(6)

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

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

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

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

y=wT·φ(x)+b

(7)

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

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

(8)

其中：

f(x)-yi≤ε+ξi

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

(9)

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

(10)

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

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

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

1.4 引力搜索算法(GSA)

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

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

(11)

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

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

(12)

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

(13)

(14)

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

(15)

式中：randj為[0,1]范圍內(nèi)的隨機(jī)數(shù)。

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

(16)

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

假設(shè)重力質(zhì)量和慣性質(zhì)量相等，則重力質(zhì)量和慣性質(zhì)量可更新為：

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

(17)

(18)

(19)

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

(20)

(21)

在求解最大值問題時(shí)，則定義為：

(22)

(23)

將粒子i的加速度加到當(dāng)前的速度上，就可以得到現(xiàn)在的粒子速度和位置：

(24)

(25)

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

1.5 評(píng)價(jià)指標(biāo)

為了更好地反映模型的預(yù)測效果。本實(shí)驗(yàn)選取均方誤差(MSE)和納什系數(shù)(NSE)對測試集的預(yù)測結(jié)果進(jìn)行評(píng)價(jià)。

(26)

(27)

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

1.6 VMD-GSA-SVR預(yù)測模型

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

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

1) 對實(shí)測月徑流數(shù)據(jù)進(jìn)行M-K突變點(diǎn)檢測，突變點(diǎn)前的數(shù)據(jù)作為訓(xùn)練集，突變點(diǎn)后的數(shù)據(jù)作為測試集。

2) 對訓(xùn)練集數(shù)據(jù)進(jìn)行VMD分解，確定模態(tài)數(shù)K，再將測試集數(shù)據(jù)分解成同樣的K個(gè)模態(tài)分量。

3) 將訓(xùn)練集輸入SVR模型中，運(yùn)用GSA對SVR的三個(gè)參數(shù)C、gamma、ε進(jìn)行全局尋優(yōu)，用MSE作為GSA的適應(yīng)度值公式。

4) 將測試集代入經(jīng)過參數(shù)優(yōu)化的GSA-SVR模型進(jìn)行預(yù)測，并用評(píng)價(jià)指標(biāo)進(jìn)行評(píng)價(jià)。

2 實(shí)例分析

2.1 資料來源與實(shí)驗(yàn)環(huán)境搭建

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

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

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

本實(shí)驗(yàn)編程使用python 3.7，其中VMD分解使用第三方庫vmdpy；SVR使用sklearn庫；GSA與對比優(yōu)化算法使用第三方庫optimal。

2.2 訓(xùn)練集與測試集的劃分

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

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

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

2.3 月徑流的VMD分解

多數(shù)信號(hào)分解方法不能手動(dòng)選擇模態(tài)數(shù)K，因此，很多研究在實(shí)驗(yàn)開始時(shí)就將整個(gè)徑流序列進(jìn)行分解，再劃分訓(xùn)練集與測試集，這樣做實(shí)際上將測試集的未知信息帶入到了模型的訓(xùn)練過程中。鑒于此，可先對訓(xùn)練集進(jìn)行分解，再選擇和訓(xùn)練集相同的模態(tài)個(gè)數(shù)K，對測試集進(jìn)行分解。

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

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

圖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 模型預(yù)測及結(jié)果評(píng)價(jià)

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

輸入訓(xùn)練集進(jìn)行模型訓(xùn)練，得到優(yōu)化后的SVR參數(shù)，如表1所示。

表1 優(yōu)化后的SVR參數(shù)結(jié)果Tab.1 Optimized SVR parameter results

圖8和圖9分別為咸陽站、臨潼站測試集預(yù)測結(jié)果。表2為各站預(yù)測評(píng)價(jià)結(jié)果。

圖8 咸陽站測試集預(yù)測結(jié)果Fig.8 Prediction results of test set at Xianyang Station

表2 各站預(yù)測評(píng)價(jià)結(jié)果Tab.2 Prediction and evaluation results of each station

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

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

3 結(jié) 論

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

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

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

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

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