王靜 曹春正

摘要為解決PM 2.5的多站點同步預(yù)測問題，提出一種貝葉斯框架下的分層自回歸時空模型．將PM 2.5日均濃度真實值視為潛在時空過程，利用一階自回歸過程刻畫時間相關(guān)性，并基于Matérn過程捕獲空間相關(guān)性，極大程度地提高了降維和同步預(yù)測的效率．此外，還將日最高溫度、相對濕度和風(fēng)速等氣象因素作為解釋變量，用于提升PM 2.5的預(yù)測效果．借助模型的分層結(jié)構(gòu)，通過貝葉斯方法結(jié)合馬爾可夫鏈蒙特卡羅（MCMC）算法實現(xiàn)參數(shù)估計和預(yù)測過程．對北京市日均PM 2.5濃度的實證分析表明，模型在空間和時間維度上均有良好的插值或預(yù)測效果．

關(guān)鍵詞貝葉斯;分層模型;自回歸;時空模型;PM 2.5預(yù)測;馬爾可夫鏈蒙特卡羅（MCMC）

中圖分類號

X513

文獻(xiàn)標(biāo)志碼

A

收稿日期

2021-11-18

資助項目

江蘇省自然科學(xué)基金（BK20191394）;國家社會科學(xué)基金重大項目（16ZDA047）

作者簡介

王靜，女，碩士生，研究方向為PM 2.5等時空數(shù)據(jù)建模及分析.jwang_jane@163.com

曹春正（通信作者），男，教授，主要研究方向為函數(shù)型數(shù)據(jù)分析、氣象中的復(fù)雜數(shù)據(jù)建模及分析.caochunzheng@163.com

0 引言

PM 2.5作為主要空氣污染物之一，由于粒徑小，可以直接被人體吸入，且在大氣中的停留時間長、輸送距離遠(yuǎn)，因而對人體健康和大氣環(huán)境質(zhì)量的影響很大．醫(yī)學(xué)研究表明，過高濃度的PM 2.5不僅會導(dǎo)致心肺疾病發(fā)病率和死亡率升高［1］，而且還會對人體的心血管系統(tǒng)、神經(jīng)系統(tǒng)和免疫系統(tǒng)造成影響［2-3］，甚至?xí)θ旧w和 DNA 等不同水平的遺傳物質(zhì)產(chǎn)生毒性作用，引起癌癥和出生缺陷的發(fā)生［4-5］．

PM 2.5的研究工作包括數(shù)據(jù)采集方法、機理成因和影響因素等［6-7］．從統(tǒng)計學(xué)的角度，一個地區(qū)一段時間的PM 2.5濃度作為典型的時空數(shù)據(jù)集，相關(guān)研究重點是空間插值以及時間上的短期或長期預(yù)測等．PM 2.5的空間插值較流行的是時空克里金方法［8-9］，該方法能夠基于時空數(shù)據(jù)時空位置關(guān)系與時空變異特征，實現(xiàn)對未觀測位置的線性無偏最優(yōu)估計，而PM 2.5在時間維度上的預(yù)測可使用機理分析或統(tǒng)計建模方法等．機理分析方法主要對大氣污染物的產(chǎn)生、轉(zhuǎn)換和擴散的物理化學(xué)過程進(jìn)行建模，如CMAQ模型［10］等;統(tǒng)計建模方法即通過數(shù)據(jù)捕捉特征，得到污染物濃度變化規(guī)律，包括多元線性回歸（Multivariable Linear Regression，MLR）［11］、廣義加性模型（Generalized Additive Model，GAM）［12-13］，以及統(tǒng)計學(xué)習(xí)模型，如BP神經(jīng)網(wǎng)絡(luò)［14］和長短期記憶網(wǎng)絡(luò)模型（Long Short-Term Memory，LSTM）［15-16］等的各種擴展模型．相比機理分析方法，統(tǒng)計方法較少依賴于污染源數(shù)據(jù)、傳輸模式和物理機理，更側(cè)重數(shù)據(jù)本身的規(guī)律，定量分析的精確度更具優(yōu)勢，是處理復(fù)雜數(shù)據(jù)的強有力工具．

近年來，許多研究關(guān)注于PM 2.5濃度的時空特征和統(tǒng)計推斷．例如：Cheam等［17］將EM算法應(yīng)用于參數(shù)時空混合模型的推斷中，用于對空氣質(zhì)量數(shù)據(jù)進(jìn)行聚類;Clifford等［18］在半?yún)?shù)時空模型的基礎(chǔ)上，利用高斯馬爾可夫隨機場對空間隨機效應(yīng)和非參數(shù)時間趨勢進(jìn)行近似，對大氣顆粒物濃度進(jìn)行預(yù)測等貝葉斯推斷．這些研究更多關(guān)注模型和計算的靈活性，沒有考慮在觸發(fā)空氣污染方面起重要作用的氣象變量．也有一些研究開發(fā)了包含氣象變量的時空模型，并將其用于時空預(yù)測［19-20］，包括Wan等［21］通過建立精細(xì)的參數(shù)統(tǒng)計模型對北京市PM 2.5濃度進(jìn)行綜合研究，分析PM 2.5濃度的時空依賴結(jié)構(gòu)并進(jìn)行了預(yù)測．但是在應(yīng)對大規(guī)模數(shù)據(jù)，特別是多站點同步預(yù)測時，此類時空模型將面臨過高的計算復(fù)雜度難題．

本文在貝葉斯框架下，結(jié)合分層模型理論，同時考慮氣象因素，對北京市35個空氣質(zhì)量監(jiān)測點的日均PM 2.5濃度整體建立了貝葉斯分層自回歸（Bayesian Hierarchical Autoregression，BHAR）時空模型．模型有三個優(yōu)點：一是利用分層結(jié)構(gòu)描述出清晰的變量關(guān)聯(lián)性和時空結(jié)構(gòu)關(guān)系;二是利用貝葉斯方法可同時達(dá)到參數(shù)估計和多站點同步預(yù)測的目的;三是可以對除現(xiàn)有空氣質(zhì)量監(jiān)測點之外，有氣象信息的地點進(jìn)行預(yù)測，解決部分地區(qū)空氣質(zhì)量監(jiān)測點分布稀疏的問題．BHAR時空模型在潛在時空過程中同時對時間和空間相關(guān)性進(jìn)行擬合，實現(xiàn)降維，解決了傳統(tǒng)時空模型計算復(fù)雜度較高的問題．進(jìn)一步借助R軟件中的spTimer包［22］，使用馬爾可夫鏈蒙特卡羅（Markov Chain Monte Carlo，MCMC）算法對模型進(jìn)行參數(shù)估計和預(yù)測．

1 數(shù)據(jù)初步分析

本文的研究區(qū)域北京位于115.7°～117.4°E，39.4°～41.6°N，地勢西北高、東南低．西部、北部和東北部三面環(huán)山，東南部平原逐漸向渤海傾斜，見圖1．

北京市原有27個空氣質(zhì)量監(jiān)測點，2012年新增了8個PM 2.5監(jiān)測點．由于從2021年1月23日起，北京市生態(tài)環(huán)境監(jiān)測中心的空氣質(zhì)量發(fā)布平臺（http：／／zx.bjmemc.com.cn）更新了北京市的空氣質(zhì)量監(jiān)測點名稱，按照城六區(qū)、東南部、東北部、西南部、西北部五個區(qū)域?qū)ΡO(jiān)測點重新進(jìn)行了劃分（表1），同時考慮到PM 2.5污染多發(fā)生在冬季［23］，因此本文收集了北京市35個監(jiān)測點2021年2月1日至2021年3月31日每天24 h的PM 2.5質(zhì)量濃度（μg／m3）數(shù)據(jù)用于后續(xù)分析．本文將基于采集數(shù)據(jù)對PM 2.5每日24 h的平均質(zhì)量濃度（日均質(zhì)量濃度）進(jìn)行建模和預(yù)測．

為初步探索北京市PM 2.5濃度的時空分布特征，首先分別按小時計算所有站點一天中每一小時的PM 2.5平均質(zhì)量濃度，按天計算一周中每一天的平均質(zhì)量濃度，然后繪制箱線圖．由圖2a可以觀察到，PM 2.5質(zhì)量濃度在一天中呈現(xiàn)出先下降再上升的趨勢，自凌晨起逐漸下降，直至14時開始出現(xiàn)明顯上升趨勢，這可能與地形和冬季的氣候條件有關(guān)，這些條件導(dǎo)致了熱逆溫和早晨、晚上的弱風(fēng)，阻止了污染物擴散［24］;同時，還可以發(fā)現(xiàn)在早高峰和晚高峰時間段內(nèi)，PM 2.5質(zhì)量濃度有小幅度上升，可能是受到汽車尾氣的影響．圖2b顯示PM 2.5在一周中的變化也具有規(guī)律性，隨著一周的人類活動水平不斷提高，PM 2.5質(zhì)量濃度也逐漸升高，而周一和休息日，PM 2.5質(zhì)量濃度明顯較低．

接下來分析北京PM 2.5濃度的空間分布特征.首先繪制2021年2月和3月北京市的平均PM 2.5質(zhì)量濃度空間分布圖．由圖3可以觀察到北部山區(qū)的PM 2.5質(zhì)量濃度較低，其中延慶縣的濃度最低，可能是受到了山谷地形的影響，而中心城區(qū)的PM 2.5質(zhì)量濃度較高顯然也是合理的．

然后借助全局莫蘭指數(shù)分析PM 2.5質(zhì)量濃度分布的空間自相關(guān)性．全局莫蘭指數(shù)公式為

I=n∑ni=1∑nj=1w ij（y i-）（y j-）∑ni=1∑nj=1w ij∑ni=1（y i-）2，? （1）

式中w ij是距離權(quán)重，代表站點i和j之間的加權(quán)距離，y i代表PM 2.5質(zhì)量濃度值．經(jīng)過計算得到PM 2.5質(zhì)量濃度的全局莫蘭指數(shù)I=0.2，為正值，且顯著性檢驗p值等于3.9×10-6，說明北京市35個空氣質(zhì)量監(jiān)測點2—3月PM 2.5質(zhì)量濃度分布存在顯著的正相關(guān)性及空間聚集性．

進(jìn)一步收集北京地區(qū)延慶、密云、北京三個氣象臺站2021年2月1日到2021年3月31日每日的風(fēng)速（m／s）、相對濕度（％）和最高溫度（℃）數(shù)據(jù)，數(shù)據(jù)來自中國氣象數(shù)據(jù)網(wǎng)（http：／／data.cma.cn）．將三個氣象變量分別簡記為WS、RH和MT，詳細(xì)匯總信息見表2．為了分析各氣象變量和PM 2.5質(zhì)量濃度的相關(guān)性，首先為三個氣象站匹配了距離最近的空氣質(zhì)量監(jiān)測點．與延慶、密云和北京三個氣象站相匹配的空氣質(zhì)量監(jiān)測點分別為延慶夏都、懷柔新城以及大興黃村三個監(jiān)測點．分別計算PM 2.5質(zhì)量濃度與三個氣象變量的Spearman相關(guān)系數(shù)，結(jié)果見表2最后一列．結(jié)果表明，相對濕度與PM 2.5呈現(xiàn)出較強的正相關(guān)．2—3月北京的相對濕度比較低，此時大氣污染物的化學(xué)聚合作用使得PM 2.5增加的速率，高于沉降作用使得PM 2.5降低的速率．相對濕度對PM 2.5影響呈正向，即PM 2.5質(zhì)量濃度隨著相對濕度的上升而增加．風(fēng)速與PM 2.5呈現(xiàn)負(fù)相關(guān)，而溫度與PM 2.5相關(guān)性較弱．

2 BHAR時空模型

2.1 模型建立

假設(shè)Z（s，t）代表站點s在時間t的PM 2.5質(zhì)量濃度實際觀測值，其對應(yīng)的真實濃度值通過潛在的隨機過程Y（s，t）刻畫，二者滿足如下測量誤差模型：

Z（s，t）=XT（s，t）β+Y（s，t）+ε（s，t），? （2）

其中：s=s 1，s 2，…，s n為n個站點的地理位置;t=1，2，…，T為時間（d）;X（s，t）代表p維氣象變量，即X（s，t）=（x 1（s，t），x 2（s，t），…，x p（s，t））T;β為回歸系數(shù);ε（s，t）為誤差項，通常假定是白噪聲過程，即ε（s，t）～GP（0，σ2 ε）．在空間統(tǒng)計中σ2 ε 通常被稱為塊金值，當(dāng)采樣點的距離為0時，半方差函數(shù)值也應(yīng)為0，但由于存在測量誤差和空間變異，使得兩采樣點非常接近時，對應(yīng)半方差函數(shù)值不為0，即存在塊金值．

對潛在的污染物排放水平Y(jié)（s，t）建立一階自回歸模型［22］：

Y（s，t）=ρY（s，t-1）+η（s，t），（3）

其中η（s，t）為殘差隨機項，用來刻畫潛在污染物排放水平的時空隨機效應(yīng)．假定η（s，t）在時間上獨立而在空間上滿足高斯過程GP（0，Σ η），其中Σ η=σ2 η S η，σ2 η 是不隨空間變化的方差，S η代表與空間相關(guān)的協(xié)方差矩陣，通常用Matérn族相關(guān)函數(shù)來刻畫［25］．此時η（s，t）的協(xié)方差矩陣是n×n維的，而不是nT×nT維的，實現(xiàn)了降維，簡化了計算．Matérn族相關(guān)函數(shù)的一般形式為

κ（u;，ν）=12ν-1Γ（ν）（2νu）νK ν（2νu），? >0，ν>0，? （4）

其中：u=‖s i-s j‖表示監(jiān)測點s i和s j之間的距離，在這里選擇的是歐氏距離;用來控制空間相關(guān)性的衰減速度，即距離u越大，衰減速度越快;ν為控制平滑程度的參數(shù);K ν為ν階的第二類貝塞爾函數(shù)．當(dāng)ν=0.5時，Matérn族相關(guān)函數(shù)退化為指數(shù)相關(guān)函數(shù)，即κ（u;）=exp（-u）;當(dāng)ν=3／2時，κ（u;）=（1+u）exp（-u）;當(dāng)ν→∞時，Matérn族相關(guān)函數(shù)退化為高斯過程函數(shù)，即κ（u;）=exp（-2u2）．

綜上，對于實測數(shù)據(jù)，BHAR時空模型結(jié)構(gòu)如下：

Z t=X tβ+Y t+ε t， （5）

Y t=ρY t-1+η t，? （6）

其中Z t=（Z（s 1，t），…，Z（s n，t））T，Y t=（Y（s 1，t），…，Y（s n，t））T，ε t=（ε（s 1，t），…，ε（s n，t））T，η t=（η（s 1，t），…，η（s n，t））T，且ε t～N（0，σ2 ε I n），η t～N（0，σ2 η S η）．根據(jù)分層模型結(jié)構(gòu)，BHAR時空模型可分為三層：第一層表示在給定潛在時空過程和參數(shù)的條件下原始數(shù)據(jù)的分布;第二層表示參數(shù)給定的條件下潛在過程的分布Y t|Θ;第三層代表引入的參數(shù)或超參數(shù)的先驗分布．其中第二層中的過程可以添加不同水平的解釋［26］，第一水平表示真正的潛在過程，第二水平刻畫了潛在過程的時空隨機效應(yīng)．

2.2 參數(shù)估計與預(yù)測

BHAR時空模型中的待估參數(shù)為Θ=｛β，ρ，σ2 ε，σ2 η，，ν｝，利用MCMC方法進(jìn)行參數(shù)估計．除了和ν以外的其他參數(shù)都存在共軛先驗分布，即在給定先驗分布后β，ρ，σ2 ε，σ2 η 都可以求出標(biāo)準(zhǔn)的滿條件后驗分布，進(jìn)一步采用Gibbs抽樣法進(jìn)行參數(shù)估計．固定ν=0.5，對于參數(shù)采用Metropolis-Hastings（MH）算法進(jìn)行估計．

對Z（s，t）的預(yù)測可以分為三類：一是預(yù)測未知監(jiān)測點s 0在已知時間t的值;二是預(yù)測已知監(jiān)測點s在未知時間點t 0的值;三是預(yù)測未知監(jiān)測點s 0在未知時間點t 0的值．將第一類預(yù)測為空間插值，第二、三類預(yù)測為在時間上的預(yù)測．

首先介紹空間插值，在未知監(jiān)測點s 0，由方程（1）可以得到：

Z（s 0，t）=XT（s 0，t）β+Y（s 0，t）+ε（s 0，t），? （7）

其中Y（s 0，t）=ρY（s 0，t-1）+η（s 0，t）．顯然，Y（s 0，t）只能由t之前所有時間點的Y（s 0，·）順序確定，且包括Y（s 0，0）．可以根據(jù)初始條件Y 0的先驗分布來計算Y（s 0，0）．當(dāng)然如果指定Y 0為一固定常數(shù)，那么Y（s 0，0）也可以被認(rèn)為是同樣的常數(shù)［19］，因此為了簡便，Y 0通常選擇固定值．

Z（s 0，t）的預(yù)測一般基于其后驗分布π（Z（s 0，t）|Z），該后驗分布可以通過對聯(lián)合后驗分布積分得到：

π（Z（s 0，t）|Z）=∫π（Z（s 0，t）|Y（s 0，t），σ2 ε）×

π（Y（s 0，t）|Θ，Y）×

π（Θ，Y|Z）dY（s 0，t）dYdΘ，? （8）

上式中Z，Y分別代表已知時間t和監(jiān)測點s的值．預(yù)測值Z（s 0，t）的估計通過MCMC按成分抽樣法獲得，具體抽樣方法如下：

1）從后驗分布π（Θ，Y|Z）中抽取隨機樣本Θ（j），Y（j）;

2）從后驗分布π（Y（s 0，t）|Θ（j），Y（j））中抽取樣本Y（j）（s 0，t）;

3）從后驗分布π（Z（s 0，t）|Y（j）（s 0，t），σ2（j） ε）中抽取樣本Z（j）（s 0，t）．

在時間維度上的預(yù)測過程與空間插值過程類似，對某一站點s（包含現(xiàn)有監(jiān)測點或任意指定位置作為監(jiān)測點），做向前一步時間預(yù)測，同樣可以基于Z（s，T+1）的后驗分布根據(jù)類似空間插值的MCMC抽樣方法實現(xiàn)．與空間插值的不同主要體現(xiàn)在時間維度上的預(yù)測需要從具有零均值、方差為σ2 η S η的邊緣分布出發(fā)來模擬Y（s，T+1），而不是條件分布．因為從0時刻一直到時刻T，觀測點的信息已經(jīng)全部被用來獲得Y（s，T），在未來的T+1時刻，除了回歸項X（s，T+1）的新值外，已經(jīng)沒有可用于條件分布的新信息．

3 實例分析

結(jié)合北京市35個空氣質(zhì)量監(jiān)測點的空間分布情況，選取其中9個監(jiān)測點作為空間驗證集，同時選取2021年3月30日和3月31日兩天作為時間驗證集，剩下的26個監(jiān)測點數(shù)據(jù)作為訓(xùn)練集用來擬合模型．利用R軟件包spTimer同時實現(xiàn)以下參數(shù)估計和預(yù)測的計算過程．由參數(shù)估計表（表3）可以看到，WS和MT兩個變量的回歸系數(shù)β 1和β 3的估計的95％ 置信區(qū)間包含零點，因此是不顯著的．氣象變量中只有相對濕度RH的回歸系數(shù)β 2顯著且為正值，這與數(shù)據(jù)初步分析的結(jié)果一致，進(jìn)一步說明了2—3月北京的相對濕度RH對PM 2.5有顯著的正向影響．

為了評估BHAR時空模型的預(yù)測性能，采用均方根誤差RMSE和平均絕對誤差MAE兩個測量指標(biāo)對預(yù)測數(shù)據(jù)和原始數(shù)據(jù)進(jìn)行誤差對比分析．RMSE和MAE的公式如下：

RMSE=1n∑ni=1（ i-Z i）2， （9）

MAE=1n∑ni=1| i-Z i|. （10）

首先對空間驗證集的9個監(jiān)測點進(jìn)行空間插值，BHAR時空模型可以實現(xiàn)對9個監(jiān)測點的同步預(yù)測．為了進(jìn)行對比，進(jìn)一步利用R軟件的mgcv包中的gam函數(shù)對訓(xùn)練集數(shù)據(jù)進(jìn)行GAM模型擬合，預(yù)測空間驗證集中各站點的PM 2.5質(zhì)量濃度，同時計算以上兩個測量指標(biāo)．Sorek-Hamer等［12］ 利用GAM能夠擬合解釋變量與被解釋變量間非線性關(guān)系的優(yōu)勢，提高了PM 2.5質(zhì)量濃度的預(yù)測效果．對比結(jié)果見表4，可以看到BHAR時空模型的RMSE和MAE都約為GAM的1／3，證明本文提出的BHAR時空模型的空間插值效果一致優(yōu)于GAM．

進(jìn)一步，分別對空間驗證集和訓(xùn)練集中的監(jiān)測點做同步時間預(yù)測，預(yù)測3月30日和31日兩天的日均PM 2.5質(zhì)量濃度值．選取LSTM作為主要對比模型，選取常規(guī)的ARIMA模型作為基準(zhǔn)對比模型．將得到的測量指標(biāo)列在表5中，可以看到ARIMA模型的預(yù)測效果最差，其次是LSTM模型，BHAR模型的時間預(yù)測效果最好．本文提出的BHAR時空模型對所有監(jiān)測點進(jìn)行整體建模，充分考慮了空間和時間相關(guān)性，因此得到了準(zhǔn)確度較高的預(yù)測結(jié)果．

4 總結(jié)

本文建立的BHAR時空模型是將一個區(qū)域的PM 2.5數(shù)據(jù)看作時間序列的空間過程，從整體上擬合PM 2.5質(zhì)量濃度的時間和空間相關(guān)性特征，實現(xiàn)了對特定站點的PM 2.5空間插值和時間上的短期預(yù)測功能，預(yù)測效果優(yōu)于GAM和LSTM．該模型不僅適用于PM 2.5質(zhì)量濃度預(yù)測，還可以推廣到其他空氣質(zhì)量地面監(jiān)測數(shù)據(jù)，例如PM 10和O 3濃度等．在貝葉斯框架下建模，對模型的不確定性更具有包容性，而且提前給定先驗分布可以融合專家知識，提高預(yù)測精度．進(jìn)一步，建立分層模型可以更加清晰地刻畫出數(shù)據(jù)內(nèi)部的潛在時空過程，增強模型的可解釋性．同時，模型的分層結(jié)構(gòu)也使得參數(shù)估計和預(yù)測等推斷過程更加便捷．

本文選取了風(fēng)速、濕度和溫度三個氣象因素作為解釋變量，用以提高模型的實際預(yù)測效果．為了簡化模型，本文采用的是固定系數(shù)，但氣象變量對PM 2.5的影響也可能會隨著時間和空間的變化而存在差異，因此在后續(xù)研究中將考慮擬合變系數(shù)模型．本文采用一階自回歸來刻畫時間維度上的相關(guān)性，達(dá)到了較好的數(shù)值效果，同時也降低了計算復(fù)雜度．在實際應(yīng)用中，可以針對數(shù)據(jù)特點，結(jié)合模型選擇方法選取合適的自回歸階數(shù)．此外，本文模型中用于刻畫空間相關(guān)性的Matérn核函數(shù)采用了齊次的歐氏距離，在有更多的地理細(xì)節(jié)特征下，可以考慮非齊次的歐氏距離或其他非歐距離，捕捉更為真實的空間相關(guān)性．

參考文獻(xiàn)

References

［1］ Pope C A III，Burnett R T，Thun M J，et al.Lung cancer，cardiopulmonary mortality，and long-term exposure to fine particulate air pollution［J］.JAMA，2002，287（9）：1132-1141

［2］ Guo Y M，Jia Y P，Pan X C，et al.The association between fine particulate air pollution and hospital emergency room visits for cardiovascular diseases in Beijing，China［J］.Science of the Total Environment，2009，407（17）：4826-4830

［3］ 郭玉明，劉利群，陳建民，等.大氣可吸入顆粒物與心腦血管疾病急診關(guān)系的病例交叉研究［J］.中華流行病學(xué)雜志，2008，29（11）：1064-1068

GUO Yuming，LIU Liqun，CHEN Jianmin，et al.Association between the concentration of particulate matters and the hospital emergency room visits for circulatory diseases：a case-crossover study［J］.Chinese Journal of Epidemiology，2008，29（11）：1064-1068

［4］ Valavanidis A，F(xiàn)iotakis K，Vlachogianni T.Airborne particulate matter and human health：toxicological assessment and importance of size and composition of particles for oxidative damage and carcinogenic mechanisms［J］.Journal of Environmental Science and Health，Part C：Environmental Carcinogenesis ＆ Ecotoxicology Reviews，2008，26（4）：339-362

［5］ Samet J M，DeMarini D M，Malling H V.Do airborne particles induce heritable mutations？ ［J］.Science，2004，304（5673）：971-972

［6］ Guo S，Hu M，Zamora M L，et al.Elucidating severe urban haze formation in China［J］.Proceedings of the National Academy of Sciences of the United States of America，2014，111（49）：17373-17378

［7］ Sun Y L，Wang Z F，Du W，et al.Long-term real-time measurements of aerosol particle composition in Beijing，China：seasonal variations，meteorological effects，and source analysis［J］.Atmospheric Chemistry and Physics，2015，15（17）：10149-10165

［8］ 盧月明，王亮，仇阿根，等.局部加權(quán)線性回歸模型的PM 2.5空間插值方法［J］.測繪科學(xué)，2018，43（11）：79-84，91

LU Yueming，WANG Liang，QIU Agen，et al.PM 2.5 spatial interpolation method based on local weighted linear regression model［J］.Science of Surveying and Mapping，2018，43（11）：79-84，91

［9］ Gething P，Atkinson P，Noor A，et al.A local space-time Kriging approach applied to a national outpatient malaria dataset［J］.Computers ＆ Geosciences，2007，33（10）：1337-1350

［10］ Byun D，Schere K L.Review of the governing equations，computational algorithms，and other components of the models-3 community multiscale air quality （CMAQ） modeling system［J］.Applied Mechanics Reviews，2006，59（2）：51-77

［11］ 李穎若，汪君霞，韓婷婷，等.利用多元線性回歸方法評估氣象條件和控制措施對APEC期間北京空氣質(zhì)量的影響［J］.環(huán)境科學(xué)，2019，40（3）：1024-1034

LI Yingruo，WANG Junxia，HAN Tingting，et al.Using multiple linear regression method to evaluate the impact of meteorological conditions and control measures on air quality in Beijing during APEC 2014［J］.Environmental Science，2019，40（3）：1024-1034

［12］ Sorek-Hamer M，Strawa A W，Chatfield R B，et al.Improved retrieval of PM 2.5 from satellite data products using non-linear methods［J］.Environmental Pollution，2013，182C：417-423

［13］ Yu S，Wang G N，Wang L，et al.Estimation and inference for generalized geoadditive models［J］.Journal of the American Statistical Association，2020，115（530）：761-774

［14］ Bai Y，Li Y，Wang X X，et al.Air pollutants concentrations forecasting using back propagation［J］.Atmospheric Pollution Research，2016，7（3）：557-566

［15］ Zhou Q P，Jiang H Y，Wang J Z，et al.A hybrid model for PM 2.5 forecasting based on ensemble empirical mode decomposition and a general regression neural network［J］.Science of the Total Environment，2014，496：264-274

［16］ 白盛楠，申曉留.基于LSTM循環(huán)神經(jīng)網(wǎng)絡(luò)的PM 2.5預(yù)測［J］.計算機應(yīng)用與軟件，2019，36（1）：67-70，104

BAI Shengnan，SHEN Xiaoliu.PM 2.5 prediction based on LSTM recurrent neural network［J］.Computer Applications and Software，2019，36（1）：67-70，104

［17］ Cheam A S M，Marbac M，McNicholas P D.Model-based clustering for spatiotemporal data on air quality monitoring［J］.Environmetrics，2017，28（3）：e2437

［18］ Clifford S，Low-Choy S，Mazaheri M，et al.A Bayesian spatiotemporal model of panel design data：airborne particle number concentration in Brisbane，Australia［J］.Environmetrics，2019，30（7）：e2597

［19］ Nicolis O，Díaz M，Sahu S K，et al.Bayesian spatiotemporal modeling for estimating short-term exposure to air pollution in Santiago de Chile［J］.Environmetrics，2019，30（7）：e2574

［20］ Padilla L，Lagos-lvarez B，Mateu J，et al.Space-time autoregressive estimation and prediction with missing data based on Kalman filtering［J］.Environmetrics，2020，31（7）：e2627

［21］ Wan Y T，Xu M Y，Huang H，et al.A spatio-temporal model for the analysis and prediction of fine particulate matter concentration in Beijing［J］.Environmetrics，2021，32（1）：e2648

［22］ Bakar K S，Sahu S K.spTimer：spatio-temporal Bayesian modeling using R［J］.Journal of Statistical Software，2015，63（15）：1-32

［23］ 梁麗思，靖娟利，王安娜，等.2014—2019年冬季京津冀地區(qū)PM 2.5質(zhì)量濃度時空分布特征［J］.桂林理工大學(xué)學(xué)報，2020，40（4）：788-797

LIANG Lisi，JING Juanli，WANG Anna，et al.Spatial-temporal distribution characteristics of PM 2.5 concentrations in Beijing-Tianjin-Hebei region in winter from 2014 to 2019［J］.Journal of Guilin University of Technology，2020，40（4）：788-797

［24］ 王晨，時悅，景悅，等.基于遙感數(shù)據(jù)的京津冀地區(qū)PM 2.5時空分布特征［J］.環(huán)境監(jiān)測管理與技術(shù)，2020（1）：37-41

WANG Chen，SHI Yue，JING Yue，et al.Spatial and temporal distribution characteristics of PM 2.5 in Beijing-Tianjin-Hebei region based on remote sensing data［J］.The Administration and Technique of Environmental Monitoring，2020（1）：37-41

［25］ Handcock M S，Stein M L.A Bayesian analysis of Kriging［J］.Technometrics，1993，35（4）：403-410

［26］ Gelfand A E.Hierarchical modeling for spatial data problems［J］.Spatial Statistics，2012，1：30-39

Prediction of PM 2.5 concentration in Beijing based on Bayesian

hierarchical autoregressive spatio-temporal model

WANG Jing1 CAO Chunzheng1

1School of Mathematics and Statistics，Nanjing University of Information Science ＆ Technology，Nanjing 210044

Abstract Here，a hierarchical autoregressive spatio-temporal model under the Bayesian framework is proposed to address the simultaneous multi-site PM 2.5 prediction.The true daily average concentration of PM 2.5 is regarded as a potential spatio-temporal process，then the temporal correlation is described by the first-order autoregressive process and the spatial correlation is captured based on the Matérn process，which greatly improves the efficiency in dimension reduction and synchronous prediction.In addition，meteorological factors such as daily maximum temperature，relative humidity and wind speed are used as explanatory variables to improve the prediction accuracy.The combination of Bayesian method and MCMC can realize parameter estimation and prediction process due to the models hierarchical structure.The empirical analysis of daily PM 2.5 concentration in Beijing shows that the proposed model has good interpolation or prediction performance in both spatial and temporal dimensions.

Key words Bayesian method;hierarchical model;autoregressive;spatio-temporal model;PM 2.5 prediction;Markov Chain Monte Carlo （MCMC）