馬俊文，嚴(yán)京海，孫瑞雯，劉保獻(xiàn)

北京市生態(tài)環(huán)境監(jiān)測(cè)中心，大氣顆粒物監(jiān)測(cè)技術(shù)北京市重點(diǎn)實(shí)驗(yàn)室，北京 100048

在環(huán)境監(jiān)測(cè)部門開展空氣質(zhì)量監(jiān)測(cè)預(yù)報(bào)分析業(yè)務(wù)中，目前采用的主流空氣預(yù)報(bào)分析模型仍然主要是空氣數(shù)值預(yù)報(bào)模型。從研究大氣流體運(yùn)動(dòng)學(xué)和應(yīng)用大氣內(nèi)部流體物理學(xué)的規(guī)律關(guān)系出發(fā)，例如利用大氣流體動(dòng)力學(xué)、熱力學(xué)等，用傳統(tǒng)數(shù)理分析方法建立相應(yīng)的具體數(shù)學(xué)大氣物理規(guī)律模型，然后可以借助大型高性能計(jì)算機(jī)的計(jì)算能力，建立空氣污染物濃度的動(dòng)態(tài)分布運(yùn)輸和擴(kuò)散數(shù)值預(yù)測(cè)模型。例如：XU等[1]基于尺度空氣質(zhì)量模式系統(tǒng)(CMAQ)提出了CMAQMOS模型，很好糾正了CMAQ中平均污染物排放清單所導(dǎo)致的系統(tǒng)性預(yù)測(cè)誤差；CHUANG等[2]應(yīng)用結(jié)合化學(xué)的天氣研究與預(yù)測(cè)模型(WRF-Chem-MADRID)對(duì)美國(guó)東南部地區(qū)空氣質(zhì)量進(jìn)行預(yù)測(cè),結(jié)果顯示在臭氧和可吸入顆粒物預(yù)測(cè)中表現(xiàn)良好。

近年來(lái)，研究人員嘗試基于機(jī)器學(xué)習(xí)的深度學(xué)習(xí)算法預(yù)測(cè)PM2.5濃度。例如，在分析時(shí)間特征維度，王瑋等[3]開展了GRU神經(jīng)網(wǎng)絡(luò)在可吸入顆粒物預(yù)測(cè)中的應(yīng)用研究，白盛楠等[4]開展了LSTM神經(jīng)網(wǎng)絡(luò)在可吸入顆粒物濃度預(yù)測(cè)中的應(yīng)用研究。然而，由于空氣質(zhì)量在時(shí)間維度具有很強(qiáng)的季節(jié)性，因此，研究人員開始將季節(jié)性分解算法和深度學(xué)習(xí)算法相結(jié)合，開展了PM2.5濃度預(yù)測(cè)模型的研究。例如，劉銘等[5]提出EMD分解算法和LSTM算法的組合算法EMD-LSTM，開展了PM2.5濃度預(yù)測(cè)模型研究；劉炳春等[6]提出小波分解變換算法Wavelet和LSTM算法的組合算法Wavelet-LSTM，開展了空氣污染物濃度預(yù)測(cè)模型研究；王曉飛等[7]提出Prophet分解方法和LSTM的組合算法Prophet-LSTM，開展了PM2.5濃度預(yù)測(cè)模型研究。

但是，空氣質(zhì)量監(jiān)測(cè)數(shù)據(jù)除了具有時(shí)間特征外，還具有空間特征，隨著污染物傳輸過(guò)程以及氣象因素的影響，各監(jiān)測(cè)站的監(jiān)測(cè)數(shù)據(jù)具有很強(qiáng)的空間傳輸演化特征。因此，研究人員開始關(guān)注基于時(shí)空特征的PM2.5預(yù)測(cè)模型。例如，宋飛揚(yáng)等[8]利用KNN算法提取目標(biāo)站點(diǎn)空間因素，提出KNN算法和LSTM算法的組合算法KNN-LSTM，開展了PM2.5濃度預(yù)測(cè)模型研究；劉旭林等[9]將地理空間網(wǎng)格化，利用CNN提取空間特征，提出CNN和Seq2seq的組合算法CNN-Seq2seq，開展了未來(lái)1 h PM2.5濃度的預(yù)測(cè)研究。

現(xiàn)有研究中時(shí)間特征提取基本是使用LSTM模型實(shí)現(xiàn)，空間特征提取基本是將非歐氏距離空間的監(jiān)測(cè)站分布數(shù)據(jù)通過(guò)地理空間網(wǎng)格化后轉(zhuǎn)換到歐氏距離空間，然后利用面向歐氏距離空間的二維卷積神經(jīng)網(wǎng)絡(luò)或KNN算法提取空間特征。2017年，Thomas Kpif提出了專門面向圖提取圖結(jié)構(gòu)特征的圖卷積神經(jīng)網(wǎng)絡(luò)(GCN)，可以直接用于處理非歐氏距離空間監(jiān)測(cè)站地理分布數(shù)據(jù)。黃偉建等[10]利用圖卷積神經(jīng)網(wǎng)絡(luò)改進(jìn)門控循環(huán)單元網(wǎng)絡(luò)的輸入部分，提出面向空氣質(zhì)量的時(shí)空混合預(yù)測(cè)模型。但是，該時(shí)空混合預(yù)測(cè)模型并未分別提取時(shí)間特征和空間特征。因此，本文基于LSTM算法構(gòu)建時(shí)間特征提取組件，基于GCN算法構(gòu)建空間特征提取組件，然后將時(shí)間特征和空間特征進(jìn)行聯(lián)合，提出預(yù)測(cè)未來(lái)1 h PM2.5濃度的LSTM-GCN組合模型。

1 LSTM-GCN組合模型構(gòu)建

1.1 LSTM模型

機(jī)器學(xué)習(xí)領(lǐng)域的深度學(xué)習(xí)網(wǎng)絡(luò)算法主要有3類：第一類是循環(huán)神經(jīng)網(wǎng)絡(luò)，適用于處理序列數(shù)據(jù)；第二類是卷積神經(jīng)網(wǎng)絡(luò)，適用于處理圖形圖像數(shù)據(jù)；第三類是全連接神經(jīng)網(wǎng)絡(luò)，適用于解決分類問(wèn)題。其中，LSTM網(wǎng)絡(luò)是循環(huán)神經(jīng)網(wǎng)絡(luò)的主要變體之一,它在充分繼承了循環(huán)神經(jīng)網(wǎng)絡(luò)遞歸運(yùn)算優(yōu)點(diǎn)的同時(shí),能夠很好地解決目前循環(huán)神經(jīng)網(wǎng)絡(luò)在應(yīng)用深度機(jī)器學(xué)習(xí)算法過(guò)程中,產(chǎn)生的梯度遞歸消失和梯度遞歸爆炸問(wèn)題。

LSTM結(jié)構(gòu)的核心思想是在循環(huán)神經(jīng)網(wǎng)絡(luò)基礎(chǔ)上，引入一個(gè)叫單元狀態(tài)的連接，最后的狀態(tài)不再簡(jiǎn)單地存儲(chǔ)，而是通過(guò)LSTM神經(jīng)網(wǎng)絡(luò)的訓(xùn)練機(jī)制來(lái)選擇更新狀態(tài)。LSTM 的訓(xùn)練機(jī)制是通過(guò)遺忘門、輸入門和輸出門來(lái)控制信息的刪減或增加[11]。

1.2 GCN模型

GCN是一種基于圖數(shù)據(jù)的深度學(xué)習(xí)方法，其核心思想是學(xué)習(xí)一個(gè)映射函數(shù)，可以使圖中的節(jié)點(diǎn)能夠聚合自身和鄰居節(jié)點(diǎn)的特征，生成節(jié)點(diǎn)的新表示。卷積神經(jīng)網(wǎng)絡(luò)可以很好地處理歐氏空間數(shù)據(jù)以提取空間特征，對(duì)于非歐氏空間里的圖數(shù)據(jù)，可以使用GCN提取空間特征。

(1)

第二步是提取特征矩陣X的空間信息,具體換算方式見式(2)。

(2)

式中：σ(·)為激活函數(shù)；W(i)為第i層權(quán)值矩陣；H(i)為第i層的激活值；H(0)=X。

1.3 LSTM-GCN組合模型

設(shè)計(jì)LSTM-GCN組合模型結(jié)構(gòu)如圖1所示。

圖1 LSTM-GCN模型結(jié)構(gòu)Fig.1 LSTM-GCN model structure

圖1中，LSTM-GCN模型包含時(shí)間特征組件和空間特征組件。時(shí)間特征組件的輸入是多條時(shí)序數(shù)據(jù)，針對(duì)每條時(shí)序數(shù)據(jù)應(yīng)用LSTM模型，然后將其輸出合并于特征聯(lián)合層，由多維輸入一維化層輸出所提取的時(shí)間特征；空間特征組件的輸入是鄰接矩陣A和特征矩陣X，然后經(jīng)過(guò)GCN模型、全連接層后，由多維輸入一維化層輸出所提取的空間特征。最后，將時(shí)間特征和空間特征進(jìn)行特征聯(lián)合后，經(jīng)全連接層輸出預(yù)測(cè)值。

選擇均方根誤差(RMSE)和平均絕對(duì)誤差(MAE)為模型預(yù)測(cè)效果評(píng)價(jià)指標(biāo)，具體換算見式(3)和式(4)。

(3)

(4)

2 時(shí)空趨勢(shì)分析和數(shù)據(jù)預(yù)處理

以北京市35個(gè)空氣質(zhì)量監(jiān)測(cè)站2018—2020年P(guān)M2.5小時(shí)監(jiān)測(cè)數(shù)據(jù)、氣象監(jiān)測(cè)數(shù)據(jù)、站點(diǎn)分布經(jīng)緯度數(shù)據(jù)為樣本，首先分析PM2.5時(shí)空變化趨勢(shì)，然后構(gòu)建LSTM-GCN模型輸入數(shù)據(jù)和預(yù)測(cè)數(shù)據(jù)。

2.1 時(shí)空趨勢(shì)分析

Prophet時(shí)序模型是Facebook公司于2017年發(fā)布的，它可以分解出時(shí)序數(shù)據(jù)的趨勢(shì)，且彌補(bǔ)了傳統(tǒng)時(shí)序模型對(duì)時(shí)序數(shù)據(jù)過(guò)于局限、需要填充缺失值、模型不靈活等不足。因此，基于該模型分解35個(gè)監(jiān)測(cè)站PM2.5濃度數(shù)據(jù)的趨勢(shì)分量，然后按年計(jì)算趨勢(shì)分量的均值，繪制用于表征時(shí)空趨勢(shì)的雷達(dá)圖，如圖2所示。

圖2 時(shí)空趨勢(shì)圖Fig.2 Temporal and spatial trends

從圖2可知，2018—2020年，在時(shí)間趨勢(shì)方面，除了第15#站點(diǎn)的趨勢(shì)分量年均值在2018年較低以外，其余34個(gè)站點(diǎn)的趨勢(shì)分量年均值均隨時(shí)間逐年下降；在空間趨勢(shì)方面，站點(diǎn)之間的趨勢(shì)分量年均值相對(duì)變化趨勢(shì)各年基本一致。其中，第15#站點(diǎn)在2018年進(jìn)行了停站改造，影響了數(shù)據(jù)完整性，導(dǎo)致趨勢(shì)分量年均值偏低。

2.2 時(shí)間特征組件輸入數(shù)據(jù)和預(yù)測(cè)數(shù)據(jù)構(gòu)建

假設(shè)基于t時(shí)刻之前的連續(xù)24 h PM2.5濃度數(shù)據(jù)，預(yù)測(cè)t時(shí)刻PM2.5濃度，則設(shè)計(jì)構(gòu)建時(shí)間特征組件輸入數(shù)據(jù)以及對(duì)應(yīng)的模型預(yù)測(cè)數(shù)據(jù)的程序流程分6個(gè)步驟。

1)對(duì)PM2.5小時(shí)數(shù)據(jù)按時(shí)間順序排序，數(shù)據(jù)有35列，分別對(duì)應(yīng)35個(gè)監(jiān)測(cè)站。

2)剔除極大和極小異常值，按列執(zhí)行歸一化，即對(duì)每個(gè)監(jiān)測(cè)站的PM2.5小時(shí)數(shù)據(jù)分別進(jìn)行歸一化。

3)定義存儲(chǔ)時(shí)間特征組件輸入數(shù)據(jù)的變量M，存儲(chǔ)LSTM-GCN模型預(yù)測(cè)數(shù)據(jù)的變量Y。

4)建立滑動(dòng)窗口，窗口大小為25，并置于數(shù)據(jù)頂端。

5)提取滑動(dòng)窗口內(nèi)的數(shù)據(jù)，前24條數(shù)據(jù)追加到變量M，第25條數(shù)據(jù)追加到變量Y。

6)如果窗口滑動(dòng)到達(dá)數(shù)據(jù)底端，則終止程序，否則，滑動(dòng)窗口向下一小時(shí)滑動(dòng)，返回第5步。

時(shí)間特征組件輸入數(shù)據(jù)是35個(gè)站點(diǎn)t時(shí)刻之前的連續(xù)24 h PM2.5濃度數(shù)據(jù)，用Mk表示，則其數(shù)據(jù)格式如式(5)所示。

(5)

式中：k=1,2,…,35。

LSTM-GCN的預(yù)測(cè)數(shù)據(jù)是35個(gè)站點(diǎn)t時(shí)刻PM2.5濃度預(yù)測(cè)值,用Yt表示，則其數(shù)據(jù)格式如式(6)所示。

(6)

2.3 空間特征組件輸入數(shù)據(jù)構(gòu)建

空間特征組件用于提取35個(gè)監(jiān)測(cè)站在空間分布上的演化趨勢(shì)特征。其中，監(jiān)測(cè)站在北京市的空間分布情況如表1所示。

表1 35個(gè)站點(diǎn)空間分布Table 1 Spatial distribution of 35 stations

空間特征組件的輸入數(shù)據(jù)有站點(diǎn)分布圖和站點(diǎn)特征數(shù)據(jù)。站點(diǎn)分布圖是以站點(diǎn)為頂點(diǎn)，以站點(diǎn)之間的距離倒數(shù)為邊權(quán)構(gòu)建的加權(quán)無(wú)向圖，站點(diǎn)分布圖用鄰接矩陣表示，其數(shù)據(jù)格式如式(7)所示。

(7)

式中：yij為站點(diǎn)i和站點(diǎn)j之間地理距離的倒數(shù)，yii=0；i,j=1,2,…,n；n=35。

站點(diǎn)之間的地理距離可以由站點(diǎn)經(jīng)緯度數(shù)據(jù)計(jì)算得到，假設(shè)站點(diǎn)P的經(jīng)緯度為(P1，P2)，站點(diǎn)Q的經(jīng)緯度為(Q1，Q2)，則計(jì)算站點(diǎn)P和站點(diǎn)Q之間的地理距離D，具體換算方式見式(8)。

D=R×arccos[sin(P2)×sin(Q2)+

cos(P2)×cos(Q2)×cos(P1-Q1)]×π/180

(8)

式中：R為赤道半徑，取值為6 371.004 km。

站點(diǎn)特征數(shù)據(jù)是t-1時(shí)刻每個(gè)站點(diǎn)所具有的特征，假設(shè)每個(gè)站點(diǎn)有m個(gè)特征，則形成的特征矩陣如式(9)所示。

(9)

式中：zij是站點(diǎn)i的第j個(gè)特征；n=35。

因?yàn)槊總€(gè)監(jiān)測(cè)站點(diǎn)除了監(jiān)測(cè)PM2.5外，還監(jiān)測(cè)相關(guān)氣象數(shù)據(jù)，且氣象條件的變化對(duì)PM2.5有影響[12]，因此，使用PM2.5、風(fēng)速、風(fēng)向、溫度、濕度構(gòu)建監(jiān)測(cè)站點(diǎn)的特征數(shù)據(jù)。數(shù)據(jù)預(yù)處理方法：對(duì)PM2.5、風(fēng)速、溫度、濕度4列數(shù)值型數(shù)據(jù)，剔除異常數(shù)據(jù)后分別執(zhí)行歸一化；對(duì)風(fēng)向數(shù)據(jù)，按照東、南、西、北、西北、東北、東南、西南8個(gè)方向，使用獨(dú)熱編碼(One-Hot Encoding)方法，編碼為由0和1組成的8位數(shù)據(jù)，每個(gè)數(shù)據(jù)位代表1個(gè)方向。例如，若8位數(shù)據(jù)中代表方向東的數(shù)據(jù)位為1，代表其他方向的數(shù)據(jù)位為0，則表示風(fēng)向?yàn)闁|風(fēng)。以上2類數(shù)據(jù)按列合并后，形成包含12列的特征數(shù)據(jù)。

按照時(shí)間維度依次提取35個(gè)站點(diǎn)同一時(shí)刻的特征數(shù)據(jù)，形成特征矩陣，則其階數(shù)為35×12階，行代表35個(gè)監(jiān)測(cè)站點(diǎn)，列代表監(jiān)測(cè)站點(diǎn)特征。

LSTM-GCN模型輸入和輸出矩陣的階數(shù)如圖3所示。

圖3 模型輸入/輸出矩陣階數(shù)Fig.3 Input/output matrix order of model

圖3中，LSTM-GCN模型輸入有3個(gè)矩陣，分別是24×35階的PM2.5時(shí)間序列矩陣、35×35階的站點(diǎn)分布圖鄰接矩陣、35×12階的監(jiān)測(cè)站特征矩陣，輸出預(yù)測(cè)值是1×35階的矩陣。

3 仿真實(shí)驗(yàn)

本文實(shí)驗(yàn)基于TensorFlow開發(fā)，它是一個(gè)基于Python語(yǔ)言的較為成熟的開放源代碼的深度學(xué)習(xí)框架，應(yīng)用較為廣泛，尤其在自然語(yǔ)言處理、推薦系統(tǒng)、圖形分類、音頻處理等領(lǐng)域應(yīng)用較為豐富[13]。

模型化的開發(fā)工具可以選擇Keras來(lái)實(shí)現(xiàn),它是一個(gè)可以在Tensorflow里被導(dǎo)入、能夠被用于開發(fā)深度學(xué)習(xí)模型的高級(jí)API,在極大程度簡(jiǎn)化模型的基本構(gòu)建設(shè)計(jì)過(guò)程和模型訓(xùn)練方法的同時(shí),還可以完全兼容在Tensorflow里已有的模型開發(fā)相關(guān)的絕大部分功能,具有易于模塊化、可彈性組合、用戶友好、易于進(jìn)行擴(kuò)展的四大優(yōu)點(diǎn)[14]。此外，針對(duì)復(fù)雜模型，可以使用Keras函數(shù)式API進(jìn)行定義，然后把模型看作層，使用傳遞張量的方法來(lái)調(diào)用模型。

本實(shí)驗(yàn)除了實(shí)現(xiàn)LSTM-GCN模型外，還分離時(shí)間特征組件和空間特征組件，形成2個(gè)獨(dú)立的LSTM模型和GCN模型，作為評(píng)估LSTM-GCN模型的對(duì)比模型。

模型訓(xùn)練和預(yù)測(cè)程序流程圖如圖4所示。

圖4 模型訓(xùn)練和預(yù)測(cè)程序流程圖Fig.4 Flow chart of modeltraining and prediction

從圖4可見，樣本數(shù)據(jù)預(yù)處理后，按照時(shí)間順序排序，劃分前80%為訓(xùn)練數(shù)據(jù)，后20%為測(cè)試數(shù)據(jù)，設(shè)定模型激活函數(shù)為ReLU[15]，損失函數(shù)為均方誤差(MSE)，優(yōu)化方法為自適應(yīng)運(yùn)動(dòng)估計(jì)算法(Adam)。

同時(shí)，將時(shí)空地理加權(quán)回歸模型(GTWR)作為對(duì)比模型，GTWR是一種通過(guò)構(gòu)建時(shí)空非平衡關(guān)系的局部線性回歸模型，設(shè)置模型解釋變量為t時(shí)刻PM2.5濃度、風(fēng)速、溫度、濕度，因變量為t+1時(shí)刻PM2.5濃度。

根據(jù)每個(gè)目標(biāo)站點(diǎn)上實(shí)際值和4個(gè)目標(biāo)模型的預(yù)測(cè)值,分別計(jì)算所對(duì)應(yīng)的4個(gè)模型的評(píng)估指標(biāo)RMSE和MAE，計(jì)算結(jié)果如圖5和圖6所示。

圖5 各監(jiān)測(cè)站RMSE對(duì)比圖Fig.5 RMSE comparison of monitoring stations

從圖5和圖6可知，4個(gè)模型對(duì)35個(gè)監(jiān)測(cè)站的預(yù)測(cè)精度具有較為一致的波動(dòng)趨勢(shì)，且LSTM-GCN模型預(yù)測(cè)精度最高，GCN模型預(yù)測(cè)精度最差。

圖6 各監(jiān)測(cè)站MAE對(duì)比圖Fig.6 MAE comparison of monitoring stations

計(jì)算35個(gè)站點(diǎn)RMSE和MAE的平均值，結(jié)果如表2所示。

表2 模型評(píng)估結(jié)果Table 2 Model evaluation results

從表2可見，LSTM-GCN模型評(píng)估結(jié)果優(yōu)于LSTM模型、GCN模型、GTWR模型，相較于LSTM模型RMSE、MAE分別降低了11.68%、7.34%；相較于GCN模型RMSE、MAE分別降低了40.22%、36.37%；相較于GTWR模型RMSE、MAE分別降低了17.52%、23.69%。

以1#站點(diǎn)為例，抽取同一時(shí)間的真實(shí)值和4個(gè)模型的預(yù)測(cè)值，繪制曲線如圖7所示，以展示時(shí)間維度模型預(yù)測(cè)效果。

從圖7可知，在時(shí)間維度，4個(gè)模型的預(yù)測(cè)值基本隨著實(shí)際值波動(dòng)，且趨勢(shì)較為一致，預(yù)測(cè)效果較好。

圖7 1#站PM2.5預(yù)測(cè)效果Fig.7 PM2.5 prediction effect of 1# station

隨機(jī)選擇某時(shí)刻，獲取35個(gè)監(jiān)測(cè)站實(shí)際值和預(yù)測(cè)值，繪制曲線如圖8所示，以展示空間維度模型預(yù)測(cè)效果。

從圖8可知，在空間維度，4個(gè)模型預(yù)測(cè)值的波動(dòng)幅度和演化趨勢(shì)同實(shí)際值較為一致，預(yù)測(cè)效果較好。

圖8 35個(gè)站PM2.5預(yù)測(cè)效果Fig.8 PM2.5 prediction effect of 35 stations

應(yīng)用LSTM-GCN模型預(yù)測(cè)北京市35個(gè)監(jiān)測(cè)站2021年1—7月的PM2.5數(shù)據(jù)，隨機(jī)抽取部分時(shí)間的預(yù)測(cè)值和實(shí)際值，繪制對(duì)比圖，其中6個(gè)站的對(duì)比圖如圖9所示。

圖9 6個(gè)監(jiān)測(cè)站的LSTM-GCN預(yù)測(cè)值和實(shí)際值對(duì)比圖Fig.9 Comparison between LSTM-GCN predicted valueand actual value of six monitoring stations

從圖9可知，6個(gè)監(jiān)測(cè)站的預(yù)測(cè)值很好地?cái)M合了實(shí)際值，說(shuō)明預(yù)測(cè)效果較好。

4 結(jié)論

本文提出的LSTM-GCN模型同時(shí)考慮空氣質(zhì)量數(shù)據(jù)時(shí)間特征和空間特征，從模型預(yù)測(cè)和評(píng)估結(jié)果來(lái)看，能夠進(jìn)行PM2.5數(shù)據(jù)的預(yù)測(cè)，且預(yù)測(cè)精度比LSTM模型、GCN模型、GTWR模型均有所提升。LSTM-GCN模型在2021年1—7月的PM2.5濃度預(yù)測(cè)效果表明，該模型具有一定的可用性。

但是，LSTM-GCN模型沒有分解PM2.5的季節(jié)性，后續(xù)可以再研究疊加一種適用的季節(jié)性分解算法，并增加LSTM和GCN的層數(shù)，進(jìn)一步優(yōu)化模型。

LSTM-GCN模型相比傳統(tǒng)數(shù)值模型具有數(shù)據(jù)需求少、運(yùn)行速度快、計(jì)算資源需求少等優(yōu)勢(shì)，有望在短時(shí)預(yù)測(cè)、數(shù)據(jù)審核、異常數(shù)據(jù)判定等業(yè)務(wù)中成為一種技術(shù)支撐手段。