基于贝叶斯分层自回归时空模型的北京PM 2.5预测

王静 曹春正

摘要为解决PM 2.5的多站点同步预测问题，提出一种贝叶斯框架下的分层自回归时空模型．将PM 2.5日均浓度真实值视为潜在时空过程，利用一阶自回归过程刻画时间相关性，并基于Matérn过程捕获空间相关性，极大程度地提高了降维和同步预测的效率．此外，还将日最高温度、相对湿度和风速等气象因素作为解释变量，用于提升PM 2.5的预测效果．借助模型的分层结构，通过贝叶斯方法结合马尔可夫链蒙特卡罗（MCMC）算法实现参数估计和预测过程．对北京市日均PM 2.5浓度的实证分析表明，模型在空间和时间维度上均有良好的插值或预测效果．

关键词贝叶斯;分层模型;自回归;时空模型;PM 2.5预测;马尔可夫链蒙特卡罗（MCMC）

中图分类号

X513

文献标志码

A

收稿日期

2021-11-18

资助项目

江苏省自然科学基金（BK20191394）;国家社会科学基金重大项目（16ZDA047）

作者简介

王静，女，硕士生，研究方向为PM 2.5等时空数据建模及分析.jwang_jane@163.com

曹春正（通信作者），男，教授，主要研究方向为函数型数据分析、气象中的复杂数据建模及分析.caochunzheng@163.com

0 引言

PM 2.5作为主要空气污染物之一，由于粒径小，可以直接被人体吸入，且在大气中的停留时间长、输送距离远，因而对人体健康和大气环境质量的影响很大．医学研究表明，过高浓度的PM 2.5不仅会导致心肺疾病发病率和死亡率升高［1］，而且还会对人体的心血管系统、神经系统和免疫系统造成影响［2-3］，甚至会对染色体和 DNA 等不同水平的遗传物质产生毒性作用，引起癌症和出生缺陷的发生［4-5］．

PM 2.5的研究工作包括数据采集方法、机理成因和影响因素等［6-7］．从统计学的角度，一个地区一段时间的PM 2.5浓度作为典型的时空数据集，相关研究重点是空间插值以及时间上的短期或长期预测等．PM 2.5的空间插值较流行的是时空克里金方法［8-9］，该方法能够基于时空数据时空位置关系与时空变异特征，实现对未观测位置的线性无偏最优估计，而PM 2.5在时间维度上的预测可使用机理分析或统计建模方法等．机理分析方法主要对大气污染物的产生、转换和扩散的物理化学过程进行建模，如CMAQ模型［10］等;统计建模方法即通过数据捕捉特征，得到污染物浓度变化规律，包括多元线性回归（Multivariable Linear Regression，MLR）［11］、广义加性模型（Generalized Additive Model，GAM）［12-13］，以及统计学习模型，如BP神经网络［14］和长短期记忆网络模型（Long Short-Term Memory，LSTM）［15-16］等的各种扩展模型．相比机理分析方法，统计方法较少依赖于污染源数据、传输模式和物理机理，更侧重数据本身的规律，定量分析的精确度更具优势，是处理复杂数据的强有力工具．

近年来，许多研究关注于PM 2.5浓度的时空特征和统计推断．例如：Cheam等［17］将EM算法应用于参数时空混合模型的推断中，用于对空气质量数据进行聚类;Clifford等［18］在半参数时空模型的基础上，利用高斯马尔可夫随机场对空间随机效应和非参数时间趋势进行近似，对大气颗粒物浓度进行预测等贝叶斯推断．这些研究更多关注模型和计算的灵活性，没有考虑在触发空气污染方面起重要作用的气象变量．也有一些研究开发了包含气象变量的时空模型，并将其用于时空预测［19-20］，包括Wan等［21］通过建立精细的参数统计模型对北京市PM 2.5浓度进行综合研究，分析PM 2.5浓度的时空依赖结构并进行了预测．但是在应对大规模数据，特别是多站点同步预测时，此类时空模型将面临过高的计算复杂度难题．

本文在贝叶斯框架下，结合分层模型理论，同时考虑气象因素，对北京市35个空气质量监测点的日均PM 2.5浓度整体建立了贝叶斯分层自回归（Bayesian Hierarchical Autoregression，BHAR）时空模型．模型有三个优点：一是利用分层结构描述出清晰的变量关联性和时空结构关系;二是利用贝叶斯方法可同时达到参数估计和多站点同步预测的目的;三是可以对除现有空气质量监测点之外，有气象信息的地点进行预测，解决部分地区空气质量监测点分布稀疏的问题．BHAR时空模型在潜在时空过程中同时对时间和空间相关性进行拟合，实现降维，解决了传统时空模型计算复杂度较高的问题．进一步借助R软件中的spTimer包［22］，使用马尔可夫链蒙特卡罗（Markov Chain Monte Carlo，MCMC）算法对模型进行参数估计和预测．

1 数据初步分析

本文的研究区域北京位于115.7°～117.4°E，39.4°～41.6°N，地势西北高、东南低．西部、北部和东北部三面环山，东南部平原逐渐向渤海倾斜，见图1．

北京市原有27个空气质量监测点，2012年新增了8个PM 2.5监测点．由于从2021年1月23日起，北京市生态环境监测中心的空气质量发布平台（http：／／zx.bjmemc.com.cn）更新了北京市的空气质量监测点名称，按照城六区、东南部、东北部、西南部、西北部五个区域对监测点重新进行了划分（表1），同时考虑到PM 2.5污染多发生在冬季［23］，因此本文收集了北京市35个监测点2021年2月1日至2021年3月31日每天24 h的PM 2.5质量浓度（μg／m3）数据用于后续分析．本文将基于采集数据对PM 2.5每日24 h的平均质量浓度（日均质量浓度）进行建模和预测．

为初步探索北京市PM 2.5浓度的时空分布特征，首先分别按小时计算所有站点一天中每一小时的PM 2.5平均质量浓度，按天计算一周中每一天的平均质量浓度，然后绘制箱线图．由图2a可以观察到，PM 2.5质量浓度在一天中呈现出先下降再上升的趋势，自凌晨起逐渐下降，直至14时开始出现明显上升趋势，这可能与地形和冬季的气候条件有关，这些条件导致了热逆温和早晨、晚上的弱风，阻止了污染物扩散［24］;同时，还可以发现在早高峰和晚高峰时间段内，PM 2.5质量浓度有小幅度上升，可能是受到汽车尾气的影响．图2b显示PM 2.5在一周中的变化也具有规律性，随着一周的人類活动水平不断提高，PM 2.5质量浓度也逐渐升高，而周一和休息日，PM 2.5质量浓度明显较低．

接下来分析北京PM 2.5浓度的空间分布特征.首先绘制2021年2月和3月北京市的平均PM 2.5质量浓度空间分布图．由图3可以观察到北部山区的PM 2.5质量浓度较低，其中延庆县的浓度最低，可能是受到了山谷地形的影响，而中心城区的PM 2.5质量浓度较高显然也是合理的．

然后借助全局莫兰指数分析PM 2.5质量浓度分布的空间自相关性．全局莫兰指数公式为

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

式中w ij是距离权重，代表站点i和j之间的加权距离，y i代表PM 2.5质量浓度值．经过计算得到PM 2.5质量浓度的全局莫兰指数I=0.2，为正值，且显著性检验p值等于3.9×10-6，说明北京市35个空气质量监测点2—3月PM 2.5质量浓度分布存在显著的正相关性及空间聚集性．

进一步收集北京地区延庆、密云、北京三个气象台站2021年2月1日到2021年3月31日每日的风速（m／s）、相对湿度（％）和最高温度（℃）数据，数据来自中国气象数据网（http：／／data.cma.cn）．将三个气象变量分别简记为WS、RH和MT，详细汇总信息见表2．为了分析各气象变量和PM 2.5质量浓度的相关性，首先为三个气象站匹配了距离最近的空气质量监测点．与延庆、密云和北京三个气象站相匹配的空气质量监测点分别为延庆夏都、怀柔新城以及大兴黄村三个监测点．分别计算PM 2.5质量浓度与三个气象变量的Spearman相关系数，结果见表2最后一列．结果表明，相对湿度与PM 2.5呈现出较强的正相关．2—3月北京的相对湿度比较低，此时大气污染物的化学聚合作用使得PM 2.5增加的速率，高于沉降作用使得PM 2.5降低的速率．相对湿度对PM 2.5影响呈正向，即PM 2.5质量浓度随着相对湿度的上升而增加．风速与PM 2.5呈现负相关，而温度与PM 2.5相关性较弱．

2 BHAR时空模型

2.1 模型建立

假设Z（s，t）代表站点s在时间t的PM 2.5质量浓度实际观测值，其对应的真实浓度值通过潜在的随机过程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;β为回归系数;ε（s，t）为误差项，通常假定是白噪声过程，即ε（s，t）～GP（0，σ2 ε）．在空间统计中σ2 ε 通常被称为块金值，当采样点的距离为0时，半方差函数值也应为0，但由于存在测量误差和空间变异，使得两采样点非常接近时，对应半方差函数值不为0，即存在块金值．

对潜在的污染物排放水平Y（s，t）建立一阶自回归模型［22］：

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

其中η（s，t）为残差随机项，用来刻画潜在污染物排放水平的時空随机效应．假定η（s，t）在时间上独立而在空间上满足高斯过程GP（0，Σ η），其中Σ η=σ2 η S η，σ2 η 是不随空间变化的方差，S η代表与空间相关的协方差矩阵，通常用Matérn族相关函数来刻画［25］．此时η（s，t）的协方差矩阵是n×n维的，而不是nT×nT维的，实现了降维，简化了计算．Matérn族相关函数的一般形式为

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

其中：u=‖s i-s j‖表示监测点s i和s j之间的距离，在这里选择的是欧氏距离;用来控制空间相关性的衰减速度，即距离u越大，衰减速度越快;ν为控制平滑程度的参数;K ν为ν阶的第二类贝塞尔函数．当ν=0.5时，Matérn族相关函数退化为指数相关函数，即κ（u;）=exp（-u）;当ν=3／2时，κ（u;）=（1+u）exp（-u）;当ν→∞时，Matérn族相关函数退化为高斯过程函数，即κ（u;）=exp（-2u2）．

综上，对于实测数据，BHAR时空模型结构如下：

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 η）．根据分层模型结构，BHAR时空模型可分为三层：第一层表示在给定潜在时空过程和参数的条件下原始数据的分布;第二层表示参数给定的条件下潜在过程的分布Y t|Θ;第三层代表引入的参数或超参数的先验分布．其中第二层中的过程可以添加不同水平的解释［26］，第一水平表示真正的潜在过程，第二水平刻画了潜在过程的时空随机效应．

2.2 参数估计与预测

BHAR时空模型中的待估参数为Θ=｛β，ρ，σ2 ε，σ2 η，，ν｝，利用MCMC方法进行参数估计．除了和ν以外的其他参数都存在共轭先验分布，即在给定先验分布后β，ρ，σ2 ε，σ2 η 都可以求出标准的满条件后验分布，进一步采用Gibbs抽样法进行参数估计．固定ν=0.5，对于参数采用Metropolis-Hastings（MH）算法进行估计．

对Z（s，t）的预测可以分为三类：一是预测未知监测点s 0在已知时间t的值;二是预测已知监测点s在未知时间点t 0的值;三是预测未知监测点s 0在未知时间点t 0的值．将第一类预测为空间插值，第二、三类预测为在时间上的预测．

首先介紹空间插值，在未知监测点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）．可以根据初始条件Y 0的先验分布来计算Y（s 0，0）．当然如果指定Y 0为一固定常数，那么Y（s 0，0）也可以被认为是同样的常数［19］，因此为了简便，Y 0通常选择固定值．

Z（s 0，t）的预测一般基于其后验分布π（Z（s 0，t）|Z），该后验分布可以通过对联合后验分布积分得到：

π（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和监测点s的值．预测值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）．

在时间维度上的预测过程与空间插值过程类似，对某一站点s（包含现有监测点或任意指定位置作为监测点），做向前一步时间预测，同样可以基于Z（s，T+1）的后验分布根据类似空间插值的MCMC抽样方法实现．与空间插值的不同主要体现在时间维度上的预测需要从具有零均值、方差为σ2 η S η的边缘分布出发来模拟Y（s，T+1），而不是条件分布．因为从0时刻一直到时刻T，观测点的信息已经全部被用来获得Y（s，T），在未来的T+1时刻，除了回归项X（s，T+1）的新值外，已经没有可用于条件分布的新信息．

3 实例分析

结合北京市35个空气质量监测点的空间分布情况，选取其中9个监测点作为空间验证集，同时选取2021年3月30日和3月31日两天作为时间验证集，剩下的26个监测点数据作为训练集用来拟合模型．利用R软件包spTimer同时实现以下参数估计和预测的计算过程．由参数估计表（表3）可以看到，WS和MT两个变量的回归系数β 1和β 3的估计的95％ 置信区间包含零点，因此是不显著的．气象变量中只有相对湿度RH的回归系数β 2显著且为正值，这与数据初步分析的结果一致，进一步说明了2—3月北京的相对湿度RH对PM 2.5有显著的正向影响．

为了评估BHAR时空模型的预测性能，采用均方根误差RMSE和平均绝对误差MAE两个测量指标对预测数据和原始数据进行误差对比分析．RMSE和MAE的公式如下：

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

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

首先对空间验证集的9个监测点进行空间插值，BHAR时空模型可以实现对9个监测点的同步预测．为了进行对比，进一步利用R软件的mgcv包中的gam函数对训练集数据进行GAM模型拟合，预测空间验证集中各站点的PM 2.5质量浓度，同时计算以上两个测量指标．Sorek-Hamer等［12］ 利用GAM能够拟合解释变量与被解释变量间非线性关系的优势，提高了PM 2.5质量浓度的预测效果．对比结果见表4，可以看到BHAR时空模型的RMSE和MAE都约为GAM的1／3，证明本文提出的BHAR时空模型的空间插值效果一致优于GAM．

进一步，分别对空间验证集和训练集中的监测点做同步时间预测，预测3月30日和31日两天的日均PM 2.5质量浓度值．选取LSTM作为主要对比模型，选取常规的ARIMA模型作为基准对比模型．将得到的测量指标列在表5中，可以看到ARIMA模型的预测效果最差，其次是LSTM模型，BHAR模型的时间预测效果最好．本文提出的BHAR时空模型对所有监测点进行整体建模，充分考虑了空间和时间相关性，因此得到了准确度较高的预测结果．

4 总结

本文建立的BHAR时空模型是将一个区域的PM 2.5数据看作时间序列的空间过程，从整体上拟合PM 2.5质量浓度的时间和空间相关性特征，实现了对特定站点的PM 2.5空间插值和时间上的短期预测功能，预测效果优于GAM和LSTM．该模型不仅适用于PM 2.5质量浓度预测，还可以推广到其他空气质量地面监测数据，例如PM 10和O 3浓度等．在贝叶斯框架下建模，对模型的不确定性更具有包容性，而且提前给定先验分布可以融合专家知识，提高预测精度．进一步，建立分层模型可以更加清晰地刻画出数据内部的潜在时空过程，增强模型的可解释性．同时，模型的分层结构也使得参数估计和预测等推断过程更加便捷．

本文选取了风速、湿度和温度三个气象因素作为解释变量，用以提高模型的实际预测效果．为了简化模型，本文采用的是固定系数，但气象变量对PM 2.5的影响也可能会随着时间和空间的变化而存在差异，因此在后续研究中将考虑拟合变系数模型．本文采用一阶自回归来刻画时间维度上的相关性，达到了较好的数值效果，同时也降低了计算复杂度．在实际应用中，可以针对数据特点，结合模型选择方法选取合适的自回归阶数．此外，本文模型中用于刻画空间相关性的Matérn核函数采用了齐次的欧氏距离，在有更多的地理细节特征下，可以考虑非齐次的欧氏距离或其他非欧距离，捕捉更为真实的空间相关性．

参考文献

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］ 郭玉明，刘利群，陈建民，等.大气可吸入颗粒物与心脑血管疾病急诊关系的病例交叉研究［J］.中华流行病学杂志，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，Fiotakis 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］ 卢月明，王亮，仇阿根，等.局部加权线性回归模型的PM 2.5空间插值方法［J］.测绘科学，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期间北京空气质量的影响［J］.环境科学，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循环神经网络的PM 2.5预测［J］.计算机应用与软件，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年冬季京津冀地区PM 2.5质量浓度时空分布特征［J］.桂林理工大学学报，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］ 王晨，时悦，景悦，等.基于遥感数据的京津冀地区PM 2.5时空分布特征［J］.环境监测管理与技术，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）