陳靜杰 孟 琦

(中國民航大學(xué)電子信息與自動化學(xué)院 天津 300300)

0 引 言

民航節(jié)能減排“十三五”規(guī)劃對民航業(yè)溫室氣體排放的要求以及國際上即將實行的碳排放交易體系，使我國機(jī)場節(jié)能減排壓力日益增加[1]。因此，準(zhǔn)確預(yù)測機(jī)場能耗意義重大。而機(jī)場非飽和狀態(tài)與飽和狀態(tài)相互轉(zhuǎn)換的生命周期特點導(dǎo)致非飽和機(jī)場能耗數(shù)據(jù)有效樣本量較少且具有非線性和非平穩(wěn)特性，預(yù)測難度較大。

能耗預(yù)測方法的研究大致分為單一模型法和組合模型法。常見單一模型有時間序列法[2]、支持向量回歸法[3-4]和神經(jīng)網(wǎng)絡(luò)法[5]等，但以上單一模型在處理非線性且非平穩(wěn)的時間序列時均存在一定的局限性。因此，能夠結(jié)合單一模型優(yōu)點的組合模型研究開始受到越來越多關(guān)注。組合模型主要有：1) 并聯(lián)式組合法,該方法利用組合權(quán)系數(shù)將單一模型預(yù)測結(jié)果進(jìn)行加權(quán)相加，組合模型是單一模型的線性疊加[6-7]。該方式忽略了模型間的非線性關(guān)系且組合權(quán)重的確定決定了模型精度的高低。2) 串聯(lián)式組合法，該方法利用前一模型的輸出作為后一模型的輸入，后一模型的輸出為最終結(jié)果，該方式的組合模型精度受前一模型輸出影響且易導(dǎo)致誤差疊加[8]。3) 基于序列分解的組合預(yù)測法，該方法將原始序列分解為一系列子序列，對子序列分別建模預(yù)測并疊加得組合預(yù)測值。如小波變換法[9]、經(jīng)驗?zāi)B(tài)分解EMD(Empirical Mode Decomposition)方法[10-11]、聚類經(jīng)驗?zāi)Ｊ椒纸釫EMD(Ensemble Empirical Mode Decomposition)方法[12]和自適應(yīng)噪聲完整集成經(jīng)驗?zāi)B(tài)分解CEEMDAN(Complete Ensemble Empirical Mode Decomposition with Adaptive Noise)方法[13]等。以上研究通過分解技術(shù)降低了原始序列的非線性與非平穩(wěn)特征。但由于各種分解算法的固有缺陷，使得受外界隨機(jī)因素影響的時間序列經(jīng)分解所得的高頻分量的非線性和非平穩(wěn)性依然較高[14]。

綜上可知，對高頻分量進(jìn)一步弱線性化和平穩(wěn)化是基于序列分解的組合預(yù)測方法精度提高的關(guān)鍵，因此，提出一種基于兩步分解法和季節(jié)差分自回歸滑動平均模型SARIMA相結(jié)合的組合預(yù)測方法。首先，利用CEEMDAN算法將原始序列分解為一系列從高頻到低頻的分量，依據(jù)各分量的樣本熵SE(Sample Entropy)識別其復(fù)雜度進(jìn)行模式重組。然后，利用變分模態(tài)分解VMD(Variational Mode Decomposition)算法把重組分量中的高頻復(fù)雜分量再次分解，得到一系列呈現(xiàn)弱非線性且相對平穩(wěn)的子序列。最后，采用SARIMA對子序列分別進(jìn)行建模，將各子序列預(yù)測值疊加得機(jī)場能耗預(yù)測值。以某非飽和機(jī)場能耗數(shù)據(jù)作為樣本進(jìn)行實驗，驗證所提方法具有較高的預(yù)測精度。

1 兩步分解方法

兩步分解法的分解過程見圖1。第一步分解，利用CEEMDAN算法將原始序列分解為一系列本征模函數(shù)IMF和一個剩余趨勢分量RES。利用SE值識別各分量的復(fù)雜度進(jìn)行模式重組，降低預(yù)測模型的運(yùn)算量。第二步分解，利用VMD算法對高頻復(fù)雜分量再次分解，得到一系列子序列VM。經(jīng)過以上兩步分解，原始時間序列被分解為一系列呈現(xiàn)弱非線性且相對平穩(wěn)的子序列。

圖1 兩步分解流程圖

其中，IMF1,IMF2,…,IMF(n-1)和RES為CEEMDAN分解所得子序列，IMF1′,IMF2′,…,IMFm′為經(jīng)樣本熵識別重組后的子序列，VM1,…,VMu為VMD分解所得子序列。

1.1 CEEMDAN-SE一步分解

定義EMD算法分解所得分量為Ek(·)，CEEMDAN算法分解所得分量為FIM，k(t)，k為分解階段。CEEMDAN-SE分解具體步驟如下：

(1) 向原始時間序列加入N次白噪聲構(gòu)造新序列:

xi(t)=x(t)+ε0wi(t)

(1)

式中：x(t)為實際時間序列，wi(t)為服從N(0,1)分布的白噪聲，i=1,2,…,N，ε0為自適應(yīng)系數(shù)。

(2)

(3) CEEMDAN的第1個分量FIM，1(t)和第1階段(k=1)的余量信號r1(t)為:

(3)

r1(t)=x(t)-FIM,1(t)

(4)

(4) 對信號r1(t)+ε1E1(wi(t))進(jìn)行EMD分解，則FIM，2(t)為:

(5)

(5) 對其余階段，即k=2,3,…,K，與步驟(3)相同，計算第k個余量信號為:

rk(t)=rk-1(t)-FIM,k(t)

(6)

計算第k+1個模態(tài)分量為:

(7)

(6) 重復(fù)步驟(5)，直至所獲余量信號不能再分解為止，即余量信號極值點數(shù)不超過兩個。最終序列分解為K個IMF分量和一個余量R(t):

(8)

(7) 計算各FIM，k(t)和R(t)的SE值，將SE值相差較小的FIM，k(t)合并為新分量。SE的具體算法見文獻(xiàn)[15]，設(shè)時間序列為{u(t),t=1,2,…,N},u(t)為FIM，k(t)或R(t)，m為嵌入維數(shù)，r為相似容限,定義X(i)=[u(1),u(2),…,u(i+m-1)],i=1,2,…,N-m+1,矢量X(i)與X(j)之間的距離為:

(9)

式中：j=1,2,…,N-m+1且i≠j，則SE表示為:

(10)

經(jīng)以上步驟,原始時間序列分解為一系列具有復(fù)雜度差異的子序列，設(shè)各子序列的SE均值γ為門限值，對SE大于γ的子序列進(jìn)行第二步分解。

1.2 VMD二步分解

VMD算法實現(xiàn)過程如下：

(2) 根據(jù)下式更新參數(shù)yk和wk。

(11)

(12)

(3) 更新參數(shù)λ:

(13)

式中：τ為噪聲容限，為達(dá)良好去噪效果可設(shè)為0。

(4) 迭代停止判斷條件如下式，若不滿足，則m加1，返回步驟(2),否則停止迭代，得k個模態(tài)分量。

(14)

經(jīng)過以上兩步分解，原始時間序列被分解為一系列呈現(xiàn)弱非線性且相對平穩(wěn)的子序列。

2 預(yù)測模型

2.1 SARIMA預(yù)測模型

季節(jié)ARIMA模型記作SARIMA(p,d,q)(P,D,Q)S，公式如下:

(15)

2.2 基于兩步分解法和SARIMA的組合預(yù)測法

CEEMDAN-SE-VMD-SARIMA組合預(yù)測方法具體過程如下：

(1) 將原始時間序列進(jìn)行兩步分解，得到一系列子序列。

(16)

(17)

式中：yi(t)為殘差修正后各子序列預(yù)測結(jié)果；y(t)為子序列預(yù)測結(jié)果相加所得的最終預(yù)測結(jié)果。

3 實例分析

3.1 數(shù)據(jù)選取及誤差指標(biāo)

實驗程序在Windows 7 64位系統(tǒng)Intel-I7 CPU、32 GB內(nèi)存的計算機(jī)上進(jìn)行，在MATLAB 2009a軟件平臺上進(jìn)行仿真。實驗數(shù)據(jù)為某非飽和機(jī)場2013年-2016年的機(jī)場綜合能耗數(shù)據(jù)，數(shù)據(jù)采樣周期為1個月。利用2013年-2015年月度能耗數(shù)據(jù)對2016年機(jī)場月度能耗進(jìn)行預(yù)測，用2016年的能耗實際數(shù)據(jù)對預(yù)測方法的預(yù)測精度進(jìn)行評價。采用平均相對誤差(RME)、均方根誤差(RMSE)和判定系數(shù)(R2)對模型預(yù)測效果進(jìn)行比較。

(18)

(19)

(20)

3.2 預(yù)測過程

3.2.1 兩步分解

首先利用CEEMDAN算法分解原始能耗時間序列，分解結(jié)果如圖2所示。

圖2 能耗原始序列和CEEMDAN分解所得子序列

由圖2可以看出，原始時間序列被分解為5個IMF分量和1個剩余趨勢分量RES。由于分解后的分量較多，若對每一個分量都進(jìn)行SARIMA建模并進(jìn)行殘差修正，計算規(guī)模較大，且會造成誤差疊加。為降低上述影響，利用SE識別各子序列的復(fù)雜度進(jìn)行模式重組。由式(10)可知，SampEn值與m和r有關(guān)，由于SE具有較好的一致性，其值的變化趨勢不受m和r影響，一般情況下取m=1或2，r=0.1～0.25SD，SD為時間序列標(biāo)準(zhǔn)差，本文選取m=2，r=0.2SD。各子序列的SE值結(jié)果如圖3所示。

圖3 各IMF及RES的樣本熵

由圖3可知，IMF3和IMF4的樣本熵值差異較小，可將其合并為一個子序列。經(jīng)過重組后子序列變?yōu)?個。結(jié)合圖2和圖3可知，IMF1和IMF2的復(fù)雜度較高依然具有較高的非線性和非平穩(wěn)性特征。

計算得IMF1和IMF2的SE值均大于γ，采用VDM算法再次分解IMF1和IMF2，結(jié)果如圖4所示。

(a) IMF1經(jīng)VMD分解所得子序列

(b) IMF2的VMD分解所得子序列圖4 VDM分解

由圖4可以看出，IMF1和IMF2經(jīng)再次分解后分別得到4個充分降低了非線性和非平穩(wěn)性的子序列。至此，經(jīng)兩步分解得到11個子序列。

3.2.2 模型比較分析

在經(jīng)過兩步分解的基礎(chǔ)上，對分解所得子序列建立不同的SARIMA模型，并對擬合殘差建模，修正預(yù)測結(jié)果，將各子序列的預(yù)測值相加得機(jī)場能耗預(yù)測值。為驗證CEEMDAN-SE-VMD-SARIMA組合預(yù)測方法對提高機(jī)場能耗預(yù)測精度的有效性，分別利用SARIMA模型、CEEMDAN-SE-ARIMA法和VMD-SARIMA法對相同的能耗時間序列進(jìn)行預(yù)測，對比分析各個方法得到的預(yù)測值。機(jī)場能耗預(yù)測結(jié)果對比如圖5所示，各模型的預(yù)測誤差指標(biāo)對比見表1。

圖5 能耗實際值和各模型預(yù)測值

表1 四種方法誤差指標(biāo)比較

根據(jù)圖5和表1可知，CEEMDAN-SE-VMD-SARIMA組合預(yù)測方法的擬合效果最好，其RME和RMSE均低于其他方法，其擬合優(yōu)度R2高于其他方法，說明本文提出的CEEMDAN-SE-VMD-SARIMA組合預(yù)測方法的預(yù)測精度最高。CEEMDAN-SE-SARIMA方法和VMD-SARIMA方法的誤差小于SARIMA模型，擬合優(yōu)度R2高于SARIMA模型，說明通過分解技術(shù)可以提高模型的預(yù)測能力。而相比于單步分解的組合預(yù)測法，CEEMDAN-SE-VMD-SARIMA組合預(yù)測方法的誤差更小，擬合優(yōu)度R2更高。這說明在處理非線性非平穩(wěn)的機(jī)場能耗時間序列時，對序列進(jìn)行弱非線性化和平穩(wěn)化預(yù)處理對提高預(yù)測精度有很大的幫助，而一步分解得到的子序列中的高頻分量依然存在著較高的非線性和非平穩(wěn)性，直接將其舍棄或?qū)⑵漕A(yù)測結(jié)果加入最終結(jié)果均會導(dǎo)致預(yù)測精度的下降。相比于單步分解法，兩步分解法能夠更好地降低各個子序列的非線性和非平穩(wěn)性，從而有效地提高模型的預(yù)測精度。

4 結(jié) 語

針對非飽和機(jī)場能耗時間序列具有的非線性和非平穩(wěn)性特點，提出一種基于兩步分解法和SARIMA的組合預(yù)測方法。

(1) CEEMDAN-SE方法優(yōu)于一般的序列分解方法，能夠有效地對原始時間序列信號進(jìn)行分解并減少建模運(yùn)算量。

(2) 采用CEEMDAN-SE-VMD兩步分解法來處理原始能耗時間序列的非線性和非平穩(wěn)性，分解效果明顯優(yōu)于一步分解法。

(3) 對子序列建立SARIMA預(yù)測模型，有效提取了序列的周期特征，并且利用擬合殘差預(yù)測值修正最初的預(yù)測值，提高了預(yù)測準(zhǔn)確性。

(4) 通過對幾種方法預(yù)測結(jié)果的對比分析，可以得出，本文方法大大提升了機(jī)場能耗預(yù)測的精度，也適用于其他非飽和型機(jī)場，可為相關(guān)機(jī)場制定能耗有關(guān)決策、能耗綜合管理和節(jié)能績效評價提供一定的參考。