基于VMD和LSTM方法的北京市PM2.5短期預(yù)測

秦喜文,王強進(jìn),王新民,郭佳靜,初 曉

1.長春工業(yè)大學(xué)數(shù)學(xué)與統(tǒng)計學(xué)院，長春 130012 2.長春工業(yè)大學(xué)研究生院，長春 130012 3.長春財經(jīng)學(xué)院信息工程學(xué)院，長春 130122

0 引言

大氣污染所帶來的危害與人類自身的健康和安全息息相關(guān)，近年來，隨著工業(yè)生產(chǎn)規(guī)模的不斷擴大和城鎮(zhèn)化進(jìn)程的快速發(fā)展，大氣環(huán)境的污染問題變得越來越嚴(yán)重，因此，人們必須對大氣污染問題予以重視[1]?！办F霾”主要由氣態(tài)污染物(二氧化硫、氮氧化物)和可吸入顆粒物(PM2.5)等組成，其中PM2.5是加重霧霾天氣污染的主要原因[2]。PM2.5又稱細(xì)顆粒物，是指環(huán)境空氣中直徑≤2.5 μm的懸浮顆粒，其成分構(gòu)成復(fù)雜，可以由硫和氮的氧化物轉(zhuǎn)化而成，還包括氨氣以及建筑灰塵等[3]。研究顯示，PM2.5具有較強的地域性和季節(jié)性差異[4]。比如我國北方地區(qū)，由于在冬天受到供暖期污染排放的影響較大，所以在冬季時空氣污染更為嚴(yán)重。人們在霧霾環(huán)境中，不論是日常出行還是身體健康狀況都受到了很大的影響[5]。因此建立準(zhǔn)確、可靠、有效的建模方法，向公眾公布未來幾天內(nèi)預(yù)期污染物的濃度是非常重要的，其可以指導(dǎo)預(yù)警決策活動。

近年來，針對PM2.5的預(yù)測問題，研究人員給出了各種數(shù)據(jù)驅(qū)動的模型，其中包括線性和非線性模型。傳統(tǒng)的線性模型，如整合移動平均自回歸模型[6-7]、多元線性回歸模型[8]均基于所研究序列是線性的假設(shè)，此外，分位數(shù)回歸[9]等模型也被廣泛應(yīng)用于PM2.5的預(yù)測領(lǐng)域。然而，在實際生活中，PM2.5時間序列的特性呈現(xiàn)出非線性和時變復(fù)雜性。因此，非線性模型的應(yīng)用也越來越多，例如，基于人工神經(jīng)網(wǎng)絡(luò)[10]、支持向量機[11]、隨機森林[12]和其他混合方法[13-14]的學(xué)習(xí)模型，可以根據(jù)非線性映射和自學(xué)習(xí)的能力進(jìn)行高精度預(yù)測。

變分模態(tài)分解(variational mode decomposition，VMD)是一種新型的自適應(yīng)分解方法[15]，能夠很好地處理非線性和非平穩(wěn)序列，自提出以來就被廣泛應(yīng)用到故障診斷與時間序列預(yù)測中[16-19]。長短時記憶模型(long short term memory neural network，LSTM)早在1997年就被提出[20]，是深度學(xué)習(xí)的經(jīng)典方法之一，由于其能更好地發(fā)現(xiàn)長期依賴關(guān)系而被廣泛用于處理序列信息，被應(yīng)用到多種研究領(lǐng)域，如時間序列預(yù)測[21-23]、文本挖掘[24]、分類[25]等方面。本文提出基于VMD和LSTM相結(jié)合的混合預(yù)測模型(VMD-LSTM)，并與其他9種PM2.5預(yù)測模型的短期預(yù)測結(jié)果進(jìn)行比較，以期為有效預(yù)測PM2.5濃度提供新思路和新方法。

1 改進(jìn)的變分模態(tài)分解

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

VMD是一種新的自適應(yīng)信號處理算法，其在處理非線性、非平穩(wěn)信號時有效解決了經(jīng)驗?zāi)B(tài)分解(empirical mode decomposition，EMD)算法中存在的模態(tài)混疊現(xiàn)象和端點效應(yīng)問題；此外，由于VMD實質(zhì)上是由多個自適應(yīng)維納濾波組組成，所以它在應(yīng)對噪聲敏感問題上也更具有優(yōu)勢[10,15]。VMD算法為了實現(xiàn)原始信號的有效分解，需要通過迭代搜索一個受約束的變分模型的最優(yōu)解，從而確定每個調(diào)幅調(diào)頻子信號的帶寬和中心頻率。具體來說，VMD算法主要可以分為構(gòu)造變分問題和求解變分問題兩個部分。

首先引入“模態(tài)”的概念。VMD算法可以把原始信號序列分解成幾個固有模態(tài)函數(shù)(intrinsic mode function，IMF)，即調(diào)幅調(diào)頻子信號uk(t)：

uk(t)=Ak(t)cos[φk(t)]。

(1)

式中：k為固有模態(tài)函數(shù)的數(shù)目；t為時間；Ak(t)為瞬時幅值，且滿足Ak(t)≥0；cos[φk(t)]為瞬時頻率；φk(t)為非遞減函數(shù)。

為了保證稀疏性，可以利用VMD算法將原始輸入信號f分解為一系列調(diào)幅調(diào)頻子信號uk，分解后得到的各模態(tài)應(yīng)該滿足重構(gòu)以后近似等于原始輸入序列的約束條件，且應(yīng)該保證每個模態(tài)的估計帶寬之和最小。

構(gòu)造變分問題的過程需要通過以下3個步驟實現(xiàn)：1)對模態(tài)函數(shù)uk進(jìn)行Hilbert變換得到其對應(yīng)的解析信號，進(jìn)而得到單邊頻譜；2)為了將模態(tài)函數(shù)的中心頻帶調(diào)整到基頻帶上，將中心頻率ωk的指數(shù)函數(shù)e-jωkt與單邊頻譜相乘；3)對解調(diào)信號進(jìn)行高斯平滑處理，以此獲得每段帶寬。此時要求解的帶約束變分問題的目標(biāo)函數(shù)為

使得∑kuk=f。

(2)

式中：?t為偏導(dǎo)函數(shù)；δ(t)為單位脈沖函數(shù)；j為虛數(shù)單位；*為卷積。

在求解受約束變分模型的最優(yōu)解時，需要將式(2)中受到約束的變分問題轉(zhuǎn)換成無約束變分問題，通過引入懲罰因子α和拉格朗日算子λ來實現(xiàn)，此時構(gòu)造的增廣Lagrange函數(shù)表達(dá)式為

L({uk},{ωk},λ)=

(3)

式中，〈·〉表示做乘積。

1.2 閾值法確定K的數(shù)目

VMD在非線性非平穩(wěn)數(shù)據(jù)分解中具有較高的精確性，并且噪聲魯棒性更強，K值的選取是否合理將會直接影響最終的預(yù)測效果[26]。通常情況下，VMD算法在初始化分解時需要人為假定分解的數(shù)目K值，此時就可以實現(xiàn)將復(fù)雜的數(shù)據(jù)適當(dāng)?shù)胤殖蒏個調(diào)幅調(diào)頻分量。若K值小于實際待分解數(shù)目，則信號分解不徹底，無法完全提取信號中蘊含的有用信息；若K值大于實際待分解數(shù)目，則會產(chǎn)生過度分解現(xiàn)象，分解結(jié)果中就會產(chǎn)生虛假分量。大量實驗結(jié)果表明，原始信號與分解后信號的誤差平方和隨著K的增加而呈現(xiàn)遞減趨勢，而迭代終止條件ε則呈現(xiàn)出先減小再增大的變化過程。因此本文提出了閾值法確定K的取值。改進(jìn)VMD算法的計算步驟如下所示。

1)輸入分解數(shù)目K=2，3，…，N，N<15；

2)在每一個分解數(shù)目下，得到分解后的固有模態(tài)函數(shù)與原始序列之間的誤差平方和，記為eK，且e2>e3>…>eN。

3)在每一個分解數(shù)目下，得到每一次分解所對應(yīng)的迭代終止條件εK；

4)給定閾值p，記錄下使得eK滿足eK≤p的K值；

5)在滿足4)的分解數(shù)目下，將滿足min(εK)所對應(yīng)的K值確定為最終的分解數(shù)目。

2 LSTM神經(jīng)網(wǎng)絡(luò)

長短時記憶模型是一種時間遞歸神經(jīng)網(wǎng)絡(luò)，是循環(huán)神經(jīng)網(wǎng)絡(luò)(recurrent neural network，RNN)的一種變體，其對RNN進(jìn)行了改進(jìn)并有效避免了常規(guī)RNN網(wǎng)絡(luò)的梯度消失問題，對于時間序列數(shù)據(jù)的預(yù)測具有非常好的表現(xiàn)能力。LSTM模型建立了一個記憶存儲結(jié)構(gòu)，用來存儲網(wǎng)絡(luò)循環(huán)過程中產(chǎn)生的一系列狀態(tài)。該模型的記憶結(jié)構(gòu)如圖1所示。

圖1 LSTM的記憶結(jié)構(gòu)

LSTM模型的門有遺忘門、輸入門和輸出門3種。其中，在時刻t，xt表示記憶單元的輸入值，ht表示記憶單元隱藏層的當(dāng)前值，因此設(shè)置輸入門、遺忘門和輸出門的初始值分別為it、ft和ot：

it=σ(Wixt+Uiht-1+bi);

ft=σ(Wfxt+Ufht-1+bf);

ot=σ(Woxt+Uoht-1+Voct+bo)。

(4)

式中：σ為sigmoid激活函數(shù)，其值域為[0,1]；W為權(quán)重矩陣；U為輸入層到隱藏層的參數(shù)矩陣；V為隱藏層到輸出層的參數(shù)矩陣；b為偏置項；下標(biāo)i、f、o分別代表輸入門、遺忘門、輸出門的簡寫標(biāo)識；ct為記憶單元的候選值。

ct=tanh(Wcxt+Ucht-1+bc)。

式中,下標(biāo)c代表記憶細(xì)胞。

從圖1可以看出該記憶存儲結(jié)構(gòu)的核心部分就是細(xì)胞的自循環(huán)，其權(quán)重值是1，并且與時間的變化情況無關(guān)。其中輸入單元可以決定信息的狀態(tài)或決定有多少信息可以輸入，而輸出單元則可以調(diào)整記憶結(jié)構(gòu)決定輸出多少信息以及輸出何種信息。此外，還可以通過遺忘單元對記憶結(jié)構(gòu)中的信息連接進(jìn)行調(diào)整，決定遺忘哪些信息，防止內(nèi)部單元信息值無限制增長。因此，LSTM模型可以較好地處理時間序列數(shù)據(jù)，消除冗余噪聲，改善RNN模型的缺陷。

3 基于VMD和LSTM方法的北京市PM2.5短期預(yù)測

3.1 數(shù)據(jù)來源及描述

本文數(shù)據(jù)來自于空氣質(zhì)量歷史數(shù)據(jù)查詢網(wǎng)站(https://www.aqistudy.cn/historydata/)，該網(wǎng)站囊括了全國各個城市的空氣質(zhì)量數(shù)據(jù)，最早可追溯到2013-12，空氣質(zhì)量數(shù)據(jù)包括9個指標(biāo)：日期、AQI(空氣質(zhì)量指數(shù))、質(zhì)量等級，以及PM2.5、PM10、SO2、CO、NO2和O3的質(zhì)量濃度。本文選取北京市PM2.5的歷史數(shù)據(jù)作為研究對象，圖2為該數(shù)據(jù)集的時序圖。

從圖2可以發(fā)現(xiàn)該數(shù)據(jù)周期性明顯，并且從2017-06起開始數(shù)據(jù)峰值明顯下降。該樣本的時間區(qū)間為2013-12-01—2019-08-26，樣本數(shù)為2 093。本文旨在建立PM2.5的短期預(yù)測模型，因此將樣本劃分為訓(xùn)練集與測試集，其中訓(xùn)練集包含2 083個樣本，測試集包含10個樣本。

圖2 北京市PM2.5數(shù)據(jù)時序圖

3.2 模型評價指標(biāo)

本文通過與主成分回歸(principal component regression, PCR)、整合移動平均自回歸(autore-gressive integrated moving average, ARIMA)、隨機森林回歸(random forest regression, RFR)、支持向量回歸(support vector regression, SVR)、最小二乘支持向量回歸(least squares support vector regression,LS-SVR)、長短時記憶網(wǎng)絡(luò)這6種基準(zhǔn)模型，以及VMD-RFR、VMD-SVR和VMD-LS-SVR 3種組合模型的預(yù)測結(jié)果進(jìn)行比較，對所提出的混合預(yù)測模型VMD-LSTM的性能進(jìn)行評價。評價指標(biāo)包括均方誤差(EMSE)、平均絕對誤差(EMAE)、平均絕對百分比誤差(EMAPE)、均方根誤差(ERMSE)、異方差調(diào)整的均方誤差(EHMSE)、異方差調(diào)整的平均絕對誤差(EHMAE)和協(xié)議指數(shù)(IA)。各指標(biāo)計算公式如下：

(5)

(6)

(7)

(8)

(9)

(10)

(11)

通常來說，IA值介于0到1之間，其越接近于1，模型的預(yù)測精度越高；而其他評價指標(biāo)的值則越接近于0，模型的預(yù)測精度越高。

3.3 預(yù)測結(jié)果比較分析

由于PM2.5數(shù)據(jù)序列的波動程度較大，呈現(xiàn)出非線性、非平穩(wěn)的復(fù)雜特性，尖峰現(xiàn)象嚴(yán)重。為了更好地擬合PM2.5預(yù)測模型，更為全面地研究其波動特征，將具有高強度降噪特點的VMD分解方法應(yīng)用于PM2.5的預(yù)測中，以提升預(yù)測精度，分解結(jié)果如圖3所示。

從圖3可以看出，PM2.5數(shù)據(jù)序列具有明顯的周期波動特征，且頻率分布由低到高，共分解為10個固有模態(tài)函數(shù)，運用VMD分解方法可以增加重構(gòu)精度和提升數(shù)據(jù)的有效信息含量。

a.IMF1；b.IMF2；c.IMF3；d.IMF4；e.IMF5；f.IMF6；g.IMF7；h.IMF8；i.IMF9；j.IMF10。

本文用樣本集的前2 083個數(shù)據(jù)作為訓(xùn)練集，分別建立了PCR、ARIMA、RFR、SVR、LS-SVR、LSTM、VMD-RFR、VMD-SVR、VMD-LS-SVR，以及本文提出的VMD-LSTM模型。PCR模型的建立，是先以數(shù)據(jù)集中的PM2.5序列作為因變量，PM10、SO2、CO、NO2和O3的質(zhì)量濃度作為自變量，運用SAS9.4版本軟件編程實現(xiàn)該過程；然后通過比較累計貢獻(xiàn)率和特征值選取了前2個主成分進(jìn)行回歸并預(yù)測未來10期的PM2.5值。ARIMA模型的建立也是使用SAS軟件完成，先通過圖2(時序圖)得到PM2.5序列有以年為單位的周期性；之后通過ADF(augmented dickey-fuller)等平穩(wěn)性檢驗也證實了該序列為非平穩(wěn)非白噪聲序列，進(jìn)行1階12步差分后使其平穩(wěn)化；最后建立ARIMA模型。RFR、SVR、LS-SVR、VMD-RFR、VMD-SVR、和VMD-LS-SVR建模過程都是通過R軟件和MATLAB軟件實現(xiàn)，建模形式與PCR相同。LSTM通過Python實現(xiàn)，VMD-LSTM則是同時應(yīng)用MATLAB和Python實現(xiàn)的。

經(jīng)過上述過程，每個模型都得到了相應(yīng)的10個預(yù)測值。將它們與原始序列以折線圖的形式呈現(xiàn)在圖4中。

圖4 后10期PM2.5預(yù)測值

表1為各個模型的評價指標(biāo)匯總表，從表1可以發(fā)現(xiàn)：VMD-LSTM模型具有最小的EMSE、EMAE、EMAPE、ERMSE和EHMAE，僅在EHMSE指標(biāo)中，VMD-LSTM模型的EHMSE稍大于RFR模型；同時，對于指標(biāo)IA而言，VMD-LSTM模型的IA值最接近于1。因此可以認(rèn)為VMD-LSTM模型對PM2.5序列的預(yù)測精度最高，是準(zhǔn)確有效的預(yù)測模型。此外，由于PM2.5數(shù)據(jù)序列的周期波動性，VMD分解方法并非對所有機器學(xué)習(xí)方法都有效，而LSTM模型特有的門記憶結(jié)構(gòu)可以較好地處理時間序列數(shù)據(jù)，因而具有更好的性能表現(xiàn)。

表1 評價指標(biāo)

4 結(jié)論與建議

1)將改進(jìn)后的VMD和LSTM相結(jié)合，提出了一種新型的時間序列預(yù)測模型VMD-LSTM，將其應(yīng)用到PM2.5序列的預(yù)測分析，同時與其他9種經(jīng)典對照模型(PCR、ARIMA、RFR、SVR、LS-SVR、LSTM、VMD-RFR、VMD-SVR、VMD-LS-SVR)的短期預(yù)測結(jié)果進(jìn)行比較。

2)基于VMD-LSTM的PM2.5短期預(yù)測結(jié)果具有更小的誤差，其預(yù)測精度最高；RFR、SVR模型在PM2.5的預(yù)測方面仍然具有較好的準(zhǔn)確性。各評價指標(biāo)的結(jié)果均顯示，VMD-LSTM模型的預(yù)測效果優(yōu)于單獨的LSTM模型的預(yù)測效果。

3)LSTM作為一種經(jīng)典的深度學(xué)習(xí)方法，為PM2.5短期時間序列預(yù)測提供了新的研究思路和途徑。與此同時，將VMD引入到PM2.5序列預(yù)測中，能夠?qū)⒎蔷€性非平穩(wěn)序列的不同尺度特征完整地分解出來，從而提高預(yù)測精度。此外，由于VMD和LSTM都是較為新型的方法，其可改進(jìn)之處還有很多，所以在時間序列預(yù)測的發(fā)展中具有遠(yuǎn)大的前景。后續(xù)將深入研究并改進(jìn)該模型，將運用更多的數(shù)據(jù)來驗證該模型的適用性。