周堯民，黃恒君

（蘭州財經(jīng)大學(xué)統(tǒng)計學(xué)院，甘肅蘭州730020）

一、引言

改革開放以來，我國的工業(yè)、交通和服務(wù)業(yè)等發(fā)展迅速，其發(fā)展水平、實力和規(guī)模均達(dá)到世界先進(jìn)水平，人們的生活水平不斷提高，物質(zhì)需求得到極大的滿足。隨著工業(yè)發(fā)展和城市開發(fā)，能源消耗急劇增加，化石燃料、汽車尾氣等排放到大氣中，大氣環(huán)境受到嚴(yán)重污染，霧霾天氣頻發(fā)，且從局部污染向區(qū)域性污染擴(kuò)散。大氣污染已經(jīng)成為制約經(jīng)濟(jì)發(fā)展的重要因素之一，影響了人們的正常生活以及社會的安定（漆威，2015；宋凱藝和卞元超，2019）[1，2]。

近十年來，PM2.5（可吸入顆粒物，即大氣中直徑小于或等于2.5 微米的顆粒物）一直是我國大部分地區(qū)環(huán)境的首要污染物，其在大氣中的含量雖然很少，但對空氣質(zhì)量和能見度的影響很大（解堊，2011）[3]。更重要的是，PM2.5直徑小，成分復(fù)雜，含大量的有毒、有害物質(zhì)，且包含多種細(xì)小顆粒物。細(xì)小顆粒物進(jìn)入人體后會附著在呼吸道及肺葉上，對人體的一系列系統(tǒng)產(chǎn)生危害，而有害氣體、重金屬等會溶解于血液中，對人體健康的傷害更大（張義和王愛君，2020）[4]。PM2.5可以在大氣中長時間停留，長期高水平的PM2.5濃度會對生態(tài)環(huán)境、公共健康及社會經(jīng)濟(jì)構(gòu)成嚴(yán)重威脅（薛濤等，2020）[5]。全球疾病負(fù)擔(dān)研究中心（GBD）發(fā)布的報告顯示，2015 年全球約有400萬人因PM2.5污染而過早死亡，而我國就有100 多萬人，遠(yuǎn)高于歐洲和北美地區(qū)（Wang et al.，2020）[6]。PM2.5污染不僅會對健康造成直接損害，而且會增加健康支出，帶來經(jīng)濟(jì)損失。因此，PM2.5污染已經(jīng)引起研究者的廣泛關(guān)注（康曉明等，2015）[7]。

準(zhǔn)確預(yù)測PM2.5濃度的變化，從數(shù)據(jù)的角度來講，就是實時提供未來某時段PM2.5濃度的相關(guān)信息，使人們及時采取防護(hù)措施，從而在一定程度上降低大氣污染對人體的危害。基于此，本文利用“分解—聚類—集成”的學(xué)習(xí)范式，提取時間序列的各部分特征和長期趨勢，構(gòu)建先分解、再聚類、后集成的預(yù)測模型，并以北京市日均PM2.5濃度序列進(jìn)行實際預(yù)測研究。本文首先利用自適應(yīng)加噪聲完備集成經(jīng)驗?zāi)B(tài)分解（CEEMDAN）將原始序列進(jìn)行分解，并依據(jù)各分量的正則化長短時記憶神經(jīng)網(wǎng)絡(luò)（ELSTM）模型的預(yù)測效果，將預(yù)測效果不好的高頻分量利用變分模態(tài)分解（VMD）進(jìn)行二次分解，然后運(yùn)用基于形狀的時間序列聚類（K-shape）算法將子序列進(jìn)行聚類，將聚類結(jié)果作為預(yù)測模塊的輸入，并對各個預(yù)測結(jié)果進(jìn)行集成得到最終值。本文構(gòu)建的CEEMDANVMD-K-ELSTM 二次分解組合模型，可以提高模型的預(yù)測精度。

通過對城市空氣污染的研究文獻(xiàn)進(jìn)行梳理可以發(fā)現(xiàn)，PM2.5濃度預(yù)測屬于時間序列研究的一類問題，而時間序列預(yù)測問題一直受到眾多學(xué)者的關(guān)注，其采用的預(yù)測方法主要有經(jīng)典統(tǒng)計模型、支持向量機(jī)模型、神經(jīng)網(wǎng)絡(luò)模型等。經(jīng)典統(tǒng)計模型側(cè)重于線性回歸。Sun 等（2013）[8]為應(yīng)對某些關(guān)鍵氣象因素以及PM2.5濃度在先驗中所呈現(xiàn)出的非高斯分布，采用服從對數(shù)正態(tài)分布、Gamma 分布和廣義極值分布（GeV）的隱馬爾科夫模型，較為準(zhǔn)確地預(yù)測了北加州某地的PM2.5濃度超限天數(shù)，有效地減少了虛假警示。龔明等（2016）[9]建立了灰色馬爾科夫鏈模型，并在此基礎(chǔ)上對殘差進(jìn)行修正，融合了灰色模型和馬爾科夫模型的優(yōu)點(diǎn)，提高了預(yù)測精度以及預(yù)測值與實際值的吻合度。沈勁等（2020）[10]采用氣象因子聚類和多元回歸方法，基于廣東省的空氣質(zhì)量數(shù)據(jù)建立了統(tǒng)計預(yù)報模型，發(fā)現(xiàn)該模型能夠較好地模擬NO2、SO2、CO、PM10、PM2.5的日均濃度水平和變化趨勢。但這類方法容易受到各種空氣污染物濃度變化的影響，而且確切掌握其變化有著較大的難度，其預(yù)測精度往往有限。因此，經(jīng)典統(tǒng)計模型難以適應(yīng)時間序列內(nèi)含的非平穩(wěn)、含噪聲等特點(diǎn)，在預(yù)測中得不到較為準(zhǔn)確的結(jié)果。

由于經(jīng)典統(tǒng)計模型存在局限性，更多的學(xué)者開始采用時間序列法對PM2.5濃度進(jìn)行預(yù)測。自回歸移動平均模型（Autoregressive Integrated Moving Average Model，ARIMA）在時間序列的線性特征提取方面表現(xiàn)較好，Jian Le 等（2012）[11]利用 ARIMA 分析了氣壓、風(fēng)速、溫度和相對濕度在細(xì)微顆粒物濃度預(yù)測中的顯著影響作用。但這類模型在非線性特征的提取上表現(xiàn)不佳，時間序列的復(fù)雜性和非線性使其不能達(dá)到令人滿意的效果。為了解決以上問題，研究者開始采用機(jī)器學(xué)習(xí)方法進(jìn)行研究。其中，支持向量機(jī)（Support Vector Machine，SVM）在小樣本、非線性及高維模式識別中具有獨(dú)特優(yōu)勢。Sun 等（2017）[12]利用主成分分析和最小二乘支持向量機(jī)的混合模型對PM2.5濃度進(jìn)行了短期預(yù)測。李龍（2014）[13]則使用最小二乘SVM 模型結(jié)合氣象因素和污染物濃度特征預(yù)測PM2.5濃度，其與傳統(tǒng)的SVM 模型相比具有更好的預(yù)測精度和泛化能力。Zhou 等（2019）[14]將多任務(wù)算法（Multi-task Learning，MTL）與多輸入支持向量機(jī)（Multi-output Support Vector Machine，M-SVM）相結(jié)合，利用MTL 對M-SVM 模型進(jìn)行訓(xùn)練，以優(yōu)化模型參數(shù)提取非線性特征，并利用臺北市多個檢測站點(diǎn)的PM2.5濃度進(jìn)行了多步預(yù)測，以驗證模型的有效性。但是，PM2.5這類時間序列的樣本量往往很大，SVM 在處理較大的數(shù)據(jù)集時存在計算量大、處理時間長等問題。因此，部分學(xué)者利用神經(jīng)網(wǎng)絡(luò)配合其他算法處理較大規(guī)模的數(shù)據(jù)并進(jìn)行海量數(shù)據(jù)計算，取得了長足的進(jìn)步。在PM2.5濃度預(yù)測領(lǐng)域，石峰等（2017）[15]建立了基于灰狼群智能最優(yōu)化算法的神經(jīng)網(wǎng)絡(luò)預(yù)測模型，從非機(jī)理模型的角度結(jié)合氣象因素和空氣污染物對上海市的PM2.5濃度進(jìn)行了預(yù)測，其模型精度優(yōu)于BP 神經(jīng)網(wǎng)絡(luò)模型和支持向量回歸（Support Vector Regression，SVR）模型。周杉杉等（2018）[16]提出基于互信息最大相關(guān)和最小冗余準(zhǔn)則并結(jié)合粒子群優(yōu)化算法的混合特征選擇算法，利用遞歸模糊神經(jīng)網(wǎng)絡(luò)以最少的特征獲得最小的預(yù)測誤差，說明該方法適用于PM2.5濃度預(yù)測。

為了充分利用不同模型的優(yōu)勢，一些學(xué)者開始研究組合模型預(yù)測方法。Wang 等（2015）[17]利用泰勒展開修正模型誤差項將神經(jīng)網(wǎng)絡(luò)與支持向量機(jī)相結(jié)合，預(yù)測了太原市的SO2和PM10濃度。機(jī)器學(xué)習(xí)在時間序列預(yù)測中的運(yùn)用本質(zhì)上是對序列特征的監(jiān)督學(xué)習(xí)，只有最大限度地學(xué)習(xí)序列的時間窗特征，才能更好地提高機(jī)器學(xué)習(xí)類方法的預(yù)測效果，而這是一般的機(jī)器學(xué)習(xí)模型難以做到的。在深度學(xué)習(xí)中，長短期記憶神經(jīng)網(wǎng)絡(luò)（Long Short-Term Memory，LSTM）在空氣質(zhì)量預(yù)測中的應(yīng)用可以克服上述問題。Huang 等（2018）[18]將卷積神經(jīng)網(wǎng)絡(luò)（Convolutional Neural Networks，CNN）模型與 LSTM 模型相互融合，利用CNN 提取過去24 小時PM2.5濃度、風(fēng)力等信息，并將其序列輸入到LSTM 預(yù)測網(wǎng)絡(luò)中，其預(yù)測誤差小于 SVR、隨機(jī)森林（Random Forest，RF）等傳統(tǒng)模型以及單獨(dú)使用CNN 或LSTM 進(jìn)行預(yù)測的誤差。白盛楠等（2019）[19]采用灰色關(guān)聯(lián)度分析方法對多個氣象、大氣污染指標(biāo)進(jìn)行了關(guān)聯(lián)度分析，并通過搭建多變量的LSTM 循環(huán)神經(jīng)網(wǎng)絡(luò)PM2.5預(yù)測模型，實現(xiàn)了PM2.5日值濃度的準(zhǔn)確預(yù)測。蔣洪迅等（2021）[20]構(gòu)建了一種集成雙向長短期記憶網(wǎng)絡(luò)的神經(jīng)網(wǎng)絡(luò)預(yù)測模 型 DLENN （Double -LSTM Ensemble Neural Network），以內(nèi)含的兩個方向LSTM 分別刻畫PM2.5濃度變化的趨勢性和周期性，并利用線性回歸復(fù)合神經(jīng)網(wǎng)絡(luò)捕捉PM2.5濃度變化的隨機(jī)性，其結(jié)果證明了DLENN 預(yù)測模型穩(wěn)定優(yōu)于其他集成模型。

鑒于結(jié)合機(jī)器學(xué)習(xí)方法建立的組合預(yù)測模型存在非平穩(wěn)數(shù)據(jù)學(xué)習(xí)能力不足以及優(yōu)化過程容易過擬合等問題，一些學(xué)者提出了“分解—集成”的研究框架，將信號分解方法用于時間序列分析。Xiong 等（2019）[21]開展的時間序列預(yù)測工作就是從信號分解出發(fā)，有效降低了時間序列的非線性和非平穩(wěn)性。黃恒君和王偉科（2020）[22]將多模態(tài)分解與深度學(xué)習(xí)相結(jié)合，并利用多視角學(xué)習(xí)，提高了模型的預(yù)測精度。蔣峰等（2021）[23]利用變分模態(tài)分解（Variational Mode Decomposition，VMD）對 PM2.5濃度序列進(jìn)行分解并引入樣本熵對其進(jìn)行重構(gòu)，采用改進(jìn)的探路者算法優(yōu)化極限學(xué)習(xí)機(jī)（Extreme Learning Machine，ELM），最后利用極限學(xué)習(xí)機(jī)對每個重構(gòu)子序列進(jìn)行預(yù)測和集成，其預(yù)測的精度和穩(wěn)健性均有顯著的提升。為了更好地改進(jìn)預(yù)測效果，部分學(xué)者提出以二次分解的方式進(jìn)一步提取數(shù)據(jù)特征。Wang 等（2017）[24]研究了二次分解結(jié)構(gòu)，將經(jīng)驗?zāi)B(tài)分解（Empirical Mode Decomposition，EMD）所產(chǎn)生的分解信號利用小波分解（Wavelet Packet Decomposition，WPD）進(jìn)一步分解獲得最終的子序列?；诖?，本文在空氣污染研究中也采用二次分解結(jié)構(gòu)，以提高PM2.5濃度序列預(yù)測精度，并利用自適應(yīng)加噪聲的完備集成經(jīng)驗?zāi)B(tài)分解（CEEMDAN），改進(jìn)集成經(jīng)驗?zāi)B(tài)分解（EEMD）所缺失的完備性。

二、理論模型設(shè)計

雖然已有學(xué)者將二次分解技術(shù)運(yùn)用到時間序列預(yù)測中，并取得了一定的效果，但二次分解方法運(yùn)算時間過長以及運(yùn)算效率低等問題仍比較突出。因此，本文將二次分解與聚類相結(jié)合，合理減少子序列的數(shù)量并將其作為LSTM 神經(jīng)網(wǎng)絡(luò)的輸入，建立了多模態(tài)集成預(yù)測模型。首先，本文對數(shù)據(jù)進(jìn)行預(yù)先處理，運(yùn)用三次樣條插值法填補(bǔ)缺失值，并將無監(jiān)督序列轉(zhuǎn)化為有監(jiān)督序列，以配合LSTM 模型的輸入。其次，本文采用CEEMDAN 方法將原序列進(jìn)行分解，以避免模態(tài)混疊現(xiàn)象，改進(jìn)EEMD 對信號分解的不完整性，并進(jìn)一步將高頻信號利用VMD 進(jìn)行分解，以提取復(fù)雜分信號的潛在特征。再次，本文利用基于形狀相似度的時間序列聚類算法（K-shape）對分信號進(jìn)行聚類，用以區(qū)分所有成分之間的差異，并依據(jù)它們的特性將數(shù)據(jù)劃分為K 類。具體而言，各分信號和殘差成分的預(yù)測值可以分成不同的類別，每個類別中的序列具有相似的特征。第四，本文將LSTM 神經(jīng)網(wǎng)絡(luò)加入正則化項，在規(guī)避遞歸神經(jīng)網(wǎng)絡(luò)梯度爆炸問題的同時，提高模型的穩(wěn)定性及泛化能力。第五，本文采用“分解—聚類—集成”框架作為組合模型的運(yùn)行機(jī)制，以更好地適應(yīng)時間序列非平穩(wěn)、高波動、含噪聲的特點(diǎn)。具體而言，為了改善PM2.5濃度序列的預(yù)測精度，本文在“分解—聚類—集成”的研究范式下，從信號分解、信號預(yù)測以及對結(jié)果進(jìn)行集成等方面對已有模型進(jìn)行優(yōu)化和改進(jìn)，構(gòu)建了CEEMDAN-VMD-K-ELSTM 二次分解組合模型。模型的基本流程如圖1 所示。

在圖1 的框架中，本文采用如下方法構(gòu)建二次分解集成預(yù)測模型，即CEEMDAN-VMD-KELSTM：（1）利用自適應(yīng)白噪聲的完備集成經(jīng)驗?zāi)B(tài)分解（CEEMDAN），對PM2.5的原序列進(jìn)行信號分解；（2）將高頻信號通過變分模態(tài)分解（VMD）進(jìn)行二次分解，進(jìn)一步提取數(shù)據(jù)序列中的非線性和非平穩(wěn)性特征；（3）運(yùn)用基于形狀相似度的時間序列聚類算法（K-shape）將二次分解后的分信號進(jìn)行聚類，以減少預(yù)測模塊的計算量及運(yùn)行時間；（4）將聚類結(jié)果作為預(yù)測模塊的輸入，通過彈性正則化長短時神經(jīng)網(wǎng)絡(luò)（ELSTM）輸出各分信號的預(yù)測值；（5）對各分信號預(yù)測結(jié)果進(jìn)行集成，取其重復(fù)10 次的結(jié)果做平均，將其最終結(jié)果與基準(zhǔn)模型進(jìn)行比較分析，從而降低模型評價結(jié)果的隨機(jī)性。

圖1 CEEMDAN-VMD-K-ELSTM 模型流程

（一）集成經(jīng)驗?zāi)B(tài)分解（EEMD）

集成經(jīng)驗?zāi)B(tài)分解（EEMD）來源于經(jīng)驗?zāi)B(tài)分解（EMD），而EMD 是一種提取信號中非線性和非平穩(wěn)特征的技術(shù)，其本質(zhì)是將信號分解為具有不同頻率的本征模態(tài)分量（IMF）（陳仁祥等，2012）[25]。但分解過程中出現(xiàn)的跳躍式變化，會使EMD 的分解結(jié)果產(chǎn)生模態(tài)混疊現(xiàn)象，即一個IMF 中包含差異極大的特征時間尺度，或者相近的特征時間尺度分布在不同的IMF 中，出現(xiàn)相鄰兩個IMF 波形混疊難以分辨，從而使分解結(jié)果失去實際意義（張袁元等，2016）[26]。EEMD則可以有效改善EMD 所產(chǎn)生的模態(tài)混疊，其步驟可以簡述為：（1）s（t）為原始信號序列，vi（t）代表第i 次實驗中添加的白噪聲序列，其分布為標(biāo)準(zhǔn)正態(tài)分布，第i 次的信號序列可以表示為Si（t）=s（t）+vi（t），其中，i=1，…，i 代表實驗的次數(shù)；（2）將分信號序列 si（t）利用EMD 進(jìn)行分解，得到其中，k=1，…，K代表分解的模態(tài)個數(shù)；（3）s（t）的k 個模態(tài)分量為IMFk，對進(jìn)行平均可以得到

（二）自適應(yīng)完備集成經(jīng)驗?zāi)B(tài)分解（CEEMDAN）

EEMD 所添加的白噪聲序列會對原始信號產(chǎn)生一定程度的破壞并有殘余，而且求均值過程的處理較為復(fù)雜。針對以上問題，Torres 等（2011）[27]提出自適應(yīng)加噪聲的完備經(jīng)驗?zāi)B(tài)分解（CEEMDAN），即引入自適應(yīng)高斯白噪聲這一概念，通過在每個階段添加有限次的自適應(yīng)白噪聲，實現(xiàn)在較少的平均次數(shù)下，其重構(gòu)誤差接近于0。CEEMDAN 可以有效避免模態(tài)混疊問題的出現(xiàn)，并彌補(bǔ)EEMD 分解不完整的缺點(diǎn)，解決EEMD 計算效率低的問題（李峰等，2016）[28]。本文歸納出CEEMDAN 的算法步驟。

第一，利用EEMD 算法分解得到第一個模態(tài)分量：

第二，在第一階段（k=1）計算第一個余量：

第三，分解R1［n］+ε1E1（ωi［n］）（i=1，…，I）到第一個模態(tài)分量，則第二個模態(tài)分量可以表示為：

第四，對于k=2，…，K，我們計算第k 個余量：

第五，分解Rk［n］+εkEk（ωi［n］）（i=1，…，I）到第k個模態(tài)分量上，則第k+1 個模態(tài)分量可以表示為：

第六，重復(fù)第四和第五步，直到殘差分量不適合被分解時停止分解。最終的余量滿足：R［n］=X［n］-其中，K 表示分解得到的固有模態(tài)函數(shù)的數(shù)量，參數(shù)X［n］表示為

（三）變分模態(tài)分解（VMD）

變分模態(tài)分解（VMD）是Dragomiretskiy 和Zosso在2014 年提出的，它是一種新的混沌數(shù)據(jù)處理技術(shù)。VMD 分解是將信號分解為K 個本征模態(tài)函數(shù)（IMF），通過尋找一系列模態(tài)及各模態(tài)的中心頻率，重構(gòu)原始數(shù)據(jù)。VMD 分解的目的是使K 個本征模態(tài)函數(shù)的帶寬之和達(dá)到最小，進(jìn)而利用L2范數(shù)的平方最小達(dá)到上述要求。因此，分信號的瞬時頻譜具有一定的現(xiàn)實物理意義，其具體步驟分為兩步。

首先，構(gòu)造變分問題。假設(shè)原始信號f 被分解為k 個分量，為保證分解序列為具有中心頻率的有限帶寬的模態(tài)分量，同時確保各模態(tài)的估計帶寬之和最小，其約束條件為所有模態(tài)之和與原始信號f 相等。具體的變分約束表達(dá)式為：

式中，｛uk｝=｛u1，u2，…，uk｝是模態(tài)，｛ωk｝=｛ω1，ω2，…，ωk｝是模態(tài)的中心頻率，δ（t）為狄利克雷函數(shù)，*為卷積運(yùn)算。

其次，求解變分問題。我們引入拉格朗日乘子λ，轉(zhuǎn)變?yōu)闊o約束變分問題，得到如下增廣拉格朗日表達(dá)式：

式中，α 為平衡參數(shù)，其作用是減少高斯噪聲的干擾（朱敏等，2018）[29]。式（7）的優(yōu)化問題采用迭代方向乘子法（ADMM），即利用式（8）至式（10），迭代更新 u、ω、λ。

（四）基于形狀的時間序列聚類算法（K-shape）

由于時序數(shù)據(jù)具有特殊性，其聚類方法與截面數(shù)據(jù)聚類有所區(qū)別，因此，本文根據(jù)動態(tài)規(guī)劃原理，對時間序列進(jìn)行扭曲，進(jìn)行必要的錯位處理，以計算出最合適的距離，并依據(jù)時序數(shù)據(jù)的形狀相似性，將形狀相似的序列聚為一類，即K-shape 聚類算法（Gravano et al.，2016）[30]。

1.時間序列形狀相似度?；ハ嚓P(guān)測度是一種統(tǒng)計度量，由此可以確定x 和y 兩個序列的相似性。要實現(xiàn)平移不變性，計算互相關(guān)時應(yīng)保持y 序列不變，并將x 在y 上滑動，計算x 的每一個位移s 的內(nèi)積。若考慮所有的移動，CCw（x，y）=（c1，c2，…，cw），我們可以得到的互相關(guān)序列長度為2m-1，則有如下定義：

我們計算出使CCw（x，y）最大的ω，進(jìn)而得到x相對于y 的最佳移動s=w-m，進(jìn)而得到距離測度：

取值范圍是[0，2]，0 表示兩個序列最相似。

2.時間序列形狀提取。時間序列分析中的許多任務(wù)依賴于通過一個序列有效地總結(jié)一組時間序列的方法，這個摘要序列通常被稱為平均序列，其在聚類中則被稱為質(zhì)心。我們的目標(biāo)是找到與類內(nèi)所有其他時間序列之間距離平方和的最小值，這就變?yōu)橐粋€優(yōu)化問題：

該式需要對類內(nèi)所有的時間序列計算一個最佳的偏移。因為這里提到的方法是用在迭代聚類當(dāng)中，所以需要把前一次計算得到的聚類中心作為參考，并把所有的序列與這個參考的序列對齊。省略式（14）的分母可以得到：

為了簡單起見，我們用向量表示此方程，并假設(shè)序列已經(jīng)進(jìn)行了歸一化處理，得到下式：

歸一化數(shù)據(jù)，令 μk=μkQ，其中是單位矩陣，O 是全幺矩陣。用 S 代替我們得到：

其中，M=QT·S·Q，最大值 μk即為求瑞利商（Rayleigh quotient）最大化問題，同時最大值為矩陣M 對應(yīng)最大特征值的特征向量。

3.基于形狀的時間序列聚類。K-shape 算法需執(zhí)行兩個步驟：（1）在分配步驟中，算法通過將每個時間序列與所有摘要序列進(jìn)行比較，并將每個時間序列分配給最接近摘要序列的一類，以更新聚類中的成員關(guān)系；（2）在細(xì)化步驟中，通過更新聚類中心，反映前一步中聚類成員的變化。算法重復(fù)這兩個步驟，直到集群成員沒有變化，或者達(dá)到允許的最大迭代次數(shù)。在賦值步驟中，算法主要依賴時間序列形狀相似度中的距離測度，而在細(xì)化步驟中，算法主要依賴時間序列形狀提取中的聚類中心進(jìn)行計算。具體的算法流程如下：

（五）ELSTM 神經(jīng)網(wǎng)絡(luò)

長短時記憶神經(jīng)網(wǎng)絡(luò)最早是由Hochreiter 和Schmidhuber 提出的，它是循環(huán)神經(jīng)網(wǎng)絡(luò)（RNN）的改進(jìn)和發(fā)展（Dragomiretskiy，2014）。LSTM 因其獨(dú)特的單元結(jié)構(gòu)，在處理長期相關(guān)關(guān)系方面具有較明顯的優(yōu)勢，其結(jié)構(gòu)如圖2 所示。

從圖2 中可以看出，LSTM 包含一個或多個儲存器以及三個自適應(yīng)乘法門，分別為輸入門、輸出門和遺忘門。其中，輸入門是控制是否允許寫入，遺忘門是控制記憶單元的值是否需要更新，輸出門是控制是否允許輸出，通過這三個門就可以實現(xiàn)信息的保存和控制。

圖2 LSTM 神經(jīng)網(wǎng)絡(luò)的神經(jīng)元結(jié)構(gòu)

在時刻t，設(shè)xt代表PM2.5的時間序列，yt代表LSTM 的預(yù)測結(jié)果，ct和ht分別為神經(jīng)元狀態(tài)值和隱藏層狀態(tài)值，則LSTM 各單元的更新情況如下：

其中，Wih、Wfh、Wch、Woh分別是隱藏層狀態(tài)值 ht的權(quán)重矩陣，Wix、Wfx、Wcx、Wox分別是時間序列 xt的權(quán)重矩陣，Wic、Wfc、Woc分別是神經(jīng)元狀態(tài)值 ct與三個門函數(shù)的對角矩陣，bi、bf、bc、bo分別是偏置向量，Wyh和by是LSTM 網(wǎng)絡(luò)的輸出權(quán)重和偏置向量，σ（*）是 sigmoid 激活函數(shù)，g（*）和 h（*）是 tanh 激活函數(shù)，Φ 是softmax 激活函數(shù)。上述模型中所需學(xué)習(xí)訓(xùn)練的參數(shù)有：各個節(jié)點(diǎn)間的有偏連接權(quán)重、神經(jīng)元內(nèi)部的輸入連接權(quán)重、神經(jīng)元遞歸連接權(quán)重。我們對權(quán)重學(xué)習(xí)設(shè)置了不同系數(shù)的正則化項，進(jìn)而使用彈網(wǎng)懲罰項將其添加到目標(biāo)函數(shù)中，以防止模型學(xué)習(xí)過程中的過擬合，即有：

當(dāng) λ1≠0，λ2≠0 時，上式即為彈網(wǎng)懲罰，可以構(gòu)成ELSTM 模型，以提高模型的泛化能力。

三、實證分析

（一）數(shù)據(jù)描述

本文所用的北京市PM2.5濃度數(shù)據(jù)來自于中國空氣質(zhì)量在線監(jiān)測平臺（http：//www.cnemc.cn）。北京四季分明，冬天寒冷干燥，夏天炎熱少雨，再加上城市汽車總量和燃煤需求量大以及周邊城市的影響，北京的霧霾天氣頻發(fā)，一年中大多數(shù)時間的PM2.5濃度都高出正常值，呈現(xiàn)出明顯的非線性和非平穩(wěn)性特征（馬忠玉和肖宏偉，2017）。因此，選取北京市PM2.5濃度序列進(jìn)行預(yù)測，可以對模型的有效性進(jìn)行全面、系統(tǒng)的檢驗。本文以2014—2018 年的日均數(shù)據(jù)作為訓(xùn)練集，以2019 年一年的數(shù)據(jù)作為測試集，對預(yù)測模型進(jìn)行有效性檢驗。

圖3 PM2.5 濃度原始序列

從圖3 中可以看出，原始數(shù)據(jù)序列具有含噪音、高波動等特點(diǎn)。將月平均數(shù)據(jù)點(diǎn)標(biāo)注在圖中可以發(fā)現(xiàn)，每個年份的數(shù)據(jù)集均呈現(xiàn)較明顯的周期性和季節(jié)性。在長期趨勢中，北京市PM2.5濃度在2014—2017 年度水平較高、極值較多，在2018—2019 年度PM2.5濃度水平有所改善，突破 200（μg/m3）的天數(shù)明顯減少，總體呈現(xiàn)下降趨勢。

表1 PM2.5 數(shù)據(jù)的描述性統(tǒng)計分析

從表1 來看，數(shù)據(jù)集為非對稱分布，偏度較高，峰度值為6.83，數(shù)據(jù)較標(biāo)準(zhǔn)正態(tài)分布更為陡峭。在數(shù)據(jù)的前期處理中，本文運(yùn)用拉伊達(dá)法則（張德然，2003）對數(shù)據(jù)中的異常點(diǎn)進(jìn)行處理，對數(shù)據(jù)原有缺失點(diǎn)和被剔除的異常點(diǎn)利用三次樣條插值法進(jìn)行重新插值，得到處理后的數(shù)據(jù)集。以上的數(shù)據(jù)分布以及描述性統(tǒng)計分析充分說明，使用基于分布理論的傳統(tǒng)統(tǒng)計預(yù)測方法難以在PM2.5濃度序列預(yù)測中取得理想的效果。同時，數(shù)據(jù)的高波動、含噪聲、含缺失以及非線性等特點(diǎn)也證明了在進(jìn)行預(yù)測分析前進(jìn)行數(shù)據(jù)預(yù)處理具有合理性和必要性。

（二）評判標(biāo)準(zhǔn)

為了檢驗聚類任務(wù)的效果，本文利用輪廓系數(shù)評價聚類結(jié)果。輪廓系數(shù)可以同時計算類內(nèi)聚集度與類間分離度，檢驗簇內(nèi)樣本緊密程度以及簇間樣本遠(yuǎn)離程度。樣本i 的輪廓系數(shù)如下：

其中，a 為某個樣本與其所在簇內(nèi)其他樣本的平均距離，b 為某個樣本與其他簇樣本的平均距離。本文計算所有樣本對應(yīng)的輪廓系數(shù)并取均值作為該聚類結(jié)果的評價指標(biāo)，其取值范圍為[-1，1]，越接近1 說明聚類效果越好。

為了檢驗?zāi)Ｐ偷挠行裕疚牟捎昧巳N誤差分析方法，即平均絕對誤差（MAE）、均方根誤差（RMSE）和平均絕對誤差百分比（MAPE），這三種誤差分析方法的計算公式如下：

其中，yt和分別代表t 時刻的真實值和預(yù)測值，N 代表測試集中的時間點(diǎn)個數(shù)，即測試集的大小。

為了進(jìn)一步從統(tǒng)計學(xué)視角對不同預(yù)測模型的水平精度進(jìn)行比較分析，本文采用Diebold-Mariano 統(tǒng)計量測試不同模型的統(tǒng)計顯著性（孫少龍，2016）。DM 檢驗的原假設(shè)是測試模型與基準(zhǔn)模型的預(yù)測精度處于同一水準(zhǔn)，備擇假設(shè)是測試模型的預(yù)測精度顯著優(yōu)于基準(zhǔn)模型，損失函數(shù)選擇均方誤差（MSE），則DM 統(tǒng)計量的定義如下：

本文在集成預(yù)測時采用了滑動窗口前向滑動測試的方式，其目的是適應(yīng)時間序列數(shù)據(jù)的測試要求，即在測試集上統(tǒng)計真實值與預(yù)測值誤差的同時，隨測試的不斷進(jìn)行，將前一天的真實值納入模型的歷史數(shù)據(jù)集中，同時更新窗寬。該方式較固定分段的測試方式和不斷更新窗口的測試方式更符合預(yù)測時間序列的實際運(yùn)用，能夠在最近歷史數(shù)據(jù)和較遠(yuǎn)歷史數(shù)據(jù)的影響之間取得平衡。

（三）預(yù)測過程

本文將原始PM2.5濃度時間序列分解為多個子序列，其結(jié)果如圖4 所示。

圖4 PM2.5 原始序列CEEMDAN 分解

在圖4 所示的時域圖中，原始序列被分解為8條分量以及1 條趨勢項（Residual），且各個分量按照頻率從高到低依次排列，不同的分解信號反映了不同的信息。高頻數(shù)據(jù)分信號1 與分信號2 包含原始序列的震蕩信息，趨勢項則包含原始PM2.5濃度序列的趨勢信息。本文使用ELSTM 模型預(yù)測每一條分量與趨勢項，在ELSTM 模型的預(yù)測過程中，利用PM2.5濃度序列的連續(xù)7 個數(shù)據(jù)點(diǎn)預(yù)測第8 個數(shù)據(jù)點(diǎn)并依次向后滑動。由于本文采用的是日均數(shù)據(jù)，考慮到實際情況，每日的PM2.5濃度與附近7 日的PM2.5濃度應(yīng)有較明顯的相關(guān)關(guān)系，故按此設(shè)置，并對數(shù)據(jù)進(jìn)行一階差分，以緩解數(shù)據(jù)的不平穩(wěn)性。算法的迭代次數(shù)為100 次，每個小批量中包含的樣本數(shù)為50。為了確保對比的公平性和有效性，以上參數(shù)將用于本文所有的對比算法中。

基于上述設(shè)定，本文利用ELSTM 對每一條分量和趨勢項進(jìn)行預(yù)測，并對高頻分量的預(yù)測結(jié)果進(jìn)行展示。在圖5 的預(yù)測結(jié)果中，預(yù)測值與真實值較為接近，但頻率波動較大，包含的非線性特征明顯，預(yù)測效果也較差。因此，為了提高模型的預(yù)測能力，本文對難以準(zhǔn)確預(yù)測的高頻分量運(yùn)用VMD 方法進(jìn)行二次分解。VMD 方法要求預(yù)先設(shè)定分解模態(tài)數(shù)K，本文通過觀察各模態(tài)中心頻率的接近程度確定K 值，即各分量中心頻率出現(xiàn)相近的值時，就認(rèn)為出現(xiàn)VMD 過分解現(xiàn)象（陳東寧，2017）。高頻信號經(jīng)VMD分解后，不同K 值下各模態(tài)分量的中心頻率如表2所示。

圖5 高頻信號序列的ELSTM 預(yù)測

表2 不同K 值對應(yīng)的中心頻率

由表2 可知，模態(tài)數(shù)為10 時出現(xiàn)了中心頻率相近的模態(tài)分量，即出現(xiàn)過分解，故分信號數(shù)量確定為9。圖6 為高頻信號經(jīng)VMD 分解后的時域圖。

圖6 高頻信號VMD 分解時域

將高頻信號做二次分解可以獲取高頻信號中的非線性特征，但這同時會使子序列的數(shù)量成倍增加，加之ELSTM 神經(jīng)網(wǎng)絡(luò)的訓(xùn)練過程本來就很復(fù)雜，這樣會使預(yù)測模塊的訓(xùn)練時間過長。因此，本文采用時間序列聚類算法（K-shape），將具有相似特征的分信號進(jìn)行聚類，并將輪廓系數(shù)作為評價指標(biāo)。高頻信號1 的聚類結(jié)果如表3 所示。

表3 不同聚類中心數(shù)對應(yīng)的輪廓系數(shù)

由表3 可知，當(dāng)聚類中心數(shù)量為3 時，輪廓系數(shù)值最大，聚類效果最好。進(jìn)一步地，本文將所有高頻信號和低頻信號的分解結(jié)果進(jìn)行聚類，以提高預(yù)測模塊的運(yùn)行效率，降低預(yù)測模塊的運(yùn)行時間。

本文利用上述運(yùn)算得到的數(shù)據(jù)以及確定的模型參數(shù)作為預(yù)測的初始參數(shù)，對每個分信號進(jìn)行數(shù)據(jù)集分割、模型學(xué)習(xí)與預(yù)測。在訓(xùn)練模型之前，考慮時間序列的順序特點(diǎn)，本文將數(shù)據(jù)轉(zhuǎn)化為監(jiān)督問題的可訓(xùn)練形式，與ELSTM 神經(jīng)網(wǎng)絡(luò)的數(shù)據(jù)輸入要求進(jìn)行匹配，并將預(yù)測結(jié)果進(jìn)行逆差分轉(zhuǎn)換，還原為目標(biāo)預(yù)測值，最后進(jìn)行集成預(yù)測，將各個分信號的預(yù)測結(jié)果進(jìn)行疊加，作為最終的預(yù)測值。模型的預(yù)測效果如圖7 所示。

圖7 多模態(tài)集成預(yù)測效果

四、模型比較和魯棒性

為了說明基于“分解—聚類—集成”研究范式的二次分解組合模型的有效性，本文將CEEMDANVMD-K-ELSTM 模型與單一模型、一次分解集成模型、二次分解集成模型進(jìn)行對比，以探究模型復(fù)雜度的提升以及分解方法的組合對預(yù)測結(jié)果的影響，并利用DM 檢驗分析模型的預(yù)測精度是否存在顯著性差異。表4 為各個模型在北京市PM2.5濃度預(yù)測中的誤差值，本文分別從RMSE、MAE 和MAPE 標(biāo)準(zhǔn)方面對所有預(yù)測模型的性能進(jìn)行評價。

表4 不同模型的預(yù)測結(jié)果對比

由表4 的預(yù)測結(jié)果可知，ELSTM 較其他神經(jīng)網(wǎng)絡(luò)模型的精度更高，預(yù)測效果更好?；凇胺纸狻垲悺伞钡慕M合模型在預(yù)測性能上均優(yōu)于單一模型，表明“分解—聚類—集成”研究范式可以有效克服PM2.5濃度數(shù)據(jù)的高波動性、非線性特征對模型預(yù)測精度造成的影響，顯著提高模型的預(yù)測能力。本文所提出的二次分解模型在RMSE 和MAE 兩類評價中取得的效果最好，表明不同分解技術(shù)的組合對預(yù)測效果也有一定的影響，CEEMDAN 方法與VMD 方法組合的二次分解方法使得實驗結(jié)果達(dá)到最優(yōu)。

為了判斷CEEMDAN-VMD-K-ELSTM 模型的預(yù)測結(jié)果是否在統(tǒng)計學(xué)上顯著優(yōu)于基準(zhǔn)模型，本文利用DM 統(tǒng)計量進(jìn)行檢驗。DM 檢驗是根據(jù)DM 統(tǒng)計量的值判斷模型之間的預(yù)測精度是否在統(tǒng)計意義上具有顯著差異，表5 為DM 檢驗結(jié)果。

表5 DM 檢驗結(jié)果

（續(xù)表5）

由表5 可知，以本文提出的CEEMDAN-VMDK-ELSTM 作為測試模型時，其預(yù)測精度在0.01 的顯著性水平上顯著優(yōu)于其他基準(zhǔn)模型，二次分解組合模型顯著優(yōu)于一次分解模型，而單一模型中的ELSTM 神經(jīng)網(wǎng)絡(luò)顯著優(yōu)于其他網(wǎng)絡(luò)。

為了進(jìn)一步驗證模型的有效性，檢驗數(shù)據(jù)變化對模型的魯棒性影響，本文利用2015 年1 月1 日至2020 年6 月31 日的日均PM2.5濃度序列檢驗?zāi)Ｐ褪欠袢阅鼙３州^好的預(yù)測性能及穩(wěn)定性。窗寬設(shè)定為 30，驗證集數(shù)據(jù)為 2020 年 3 月 31 至 2020 年 6月 31 日。圖 8 為 CEEMDAN-VMD-K-ELSTM 模型與基準(zhǔn)模型在不同數(shù)據(jù)集上的預(yù)測比較結(jié)果。

圖8 不同數(shù)據(jù)集上的預(yù)測模型誤差對比

由圖 8 可知，CEEMDAN-VMD-K-ELSTM 組合模型在RMSE 和MAE 標(biāo)準(zhǔn)下仍具有最高的預(yù)測精度，即模型能夠?qū)Σ煌腜M2.5濃度時間序列進(jìn)行較為準(zhǔn)確的預(yù)測，具有良好的魯棒性。

五、結(jié)論與展望

空氣質(zhì)量研究一直是國內(nèi)重點(diǎn)關(guān)注的問題，污染物濃度預(yù)測更是空氣質(zhì)量研究的重中之重，其不僅可以有效預(yù)防嚴(yán)重空氣污染事件的發(fā)生，還可以幫助人們及時采取應(yīng)對措施。因此，構(gòu)建一個行之有效的空氣質(zhì)量預(yù)測模型具有重要的現(xiàn)實意義。

本文引入分解算法作為預(yù)處理工具，以提取輸入原始數(shù)據(jù)的內(nèi)在特征，而分解算法和深度學(xué)習(xí)在自然語言處理、計算機(jī)視覺等領(lǐng)域已經(jīng)取得很大的成就，尤其是在空氣質(zhì)量預(yù)測方面。由于數(shù)據(jù)存在非線性、非平穩(wěn)性及波動性的特性，以往的空氣質(zhì)量預(yù)測準(zhǔn)確率并不能令人滿意。因此，本文通過將不同的模態(tài)分解技術(shù)與常用機(jī)器學(xué)習(xí)模型、神經(jīng)網(wǎng)絡(luò)模型進(jìn)行組合對比分析，提出一種二層分解多模態(tài)集成預(yù)測方法，并在“分解—集成”的研究范式下進(jìn)行進(jìn)一步的拓展，將“分解—聚類—集成”的研究范式應(yīng)用于PM2.5濃度序列預(yù)測，采用基于時間序列形狀的聚類算法將分解后的時序數(shù)據(jù)進(jìn)行聚類，提高了模型整體的運(yùn)行效率。本文對北京市PM2.5濃度的實證分析，證明了二層分解與ELSTM 神經(jīng)網(wǎng)絡(luò)的組合在時間序列預(yù)測上可以獲得更為精準(zhǔn)的預(yù)測效果。

本文提出的基于聚類的二次分解集成模型還可以應(yīng)用于其他一些較為困難的預(yù)測任務(wù)，如金融時間序列預(yù)測、風(fēng)速預(yù)測、電力消耗預(yù)測等。就空氣質(zhì)量預(yù)測問題而言，本文依然存在諸多可拓展之處。首先，由實證分析可知，數(shù)據(jù)分解在混合集成學(xué)習(xí)范式中具有重要地位，今后應(yīng)探索更高效、可行的數(shù)據(jù)分解算法。其次，在單一預(yù)測中，為保證預(yù)測精度所使用的復(fù)雜模型，其時間復(fù)雜度大大提高，今后應(yīng)探究更高效的數(shù)據(jù)壓縮算法，以加快單一預(yù)測的速度。再次，影響PM2.5濃度的溫度、空氣、濕度等因素的相關(guān)數(shù)據(jù)量都較大，這使得精確預(yù)測PM2.5濃度存在困難，今后應(yīng)加入氣象數(shù)據(jù)、地理信息數(shù)據(jù)等影響空氣污染的其他因素，以提高預(yù)測精度，這有助于研究者提高對空氣質(zhì)量預(yù)測問題的認(rèn)知，構(gòu)建更廣義的數(shù)據(jù)融合預(yù)測模型。此外，人口規(guī)模、經(jīng)濟(jì)產(chǎn)業(yè)結(jié)構(gòu)與量級、政府空氣污染治理政策等因素也是除空氣和氣象地理數(shù)據(jù)之外應(yīng)考慮的重要方面，如何將這些因素與空氣質(zhì)量研究相結(jié)合也是一個需要解決的實際問題。