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/m³）数据用于后续分析.本文将基于采集数据对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质量浓度分布的空间自相关性.全局莫兰指数公式为

式中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相关性较弱.

注:***表示p＜0.001.

2 BHAR时空模型

2.1 模型建立

假设Z（s，t）代表站点s在时间t的PM_2.5质量浓度实际观测值，其对应的真实浓度值通过潜在的随机过程Y（s，t）刻画，二者满足如下测量误差模型:

其中:s =s₁，s₂，···，s_n为n个站点的地理位置; t =1，2，···，T为时间（d）; X（s，t）代表p维气象变量，即$\mathit{X}（s，t）={\left(x_1（s，t），x_2（s，t），\cdots ，x_p（s，t）\right)}^{\mathrm{T}}$; β为回归系数; ε（s，t）为误差项，通常假定是白噪声过程，即ε（s，t）~GP（0，${\sigma }_{\epsilon }^2$）.在空间统计中${\sigma }_{\epsilon }^2$通常被称为块金值，当采样点的距离为0时，半方差函数值也应为0，但由于存在测量误差和空间变异，使得两采样点非常接近时，对应半方差函数值不为0，即存在块金值.

对潜在的污染物排放水平Y（s，t）建立一阶自回归模型^[22]:

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

其中:u =‖s_i-s_j‖表示监测点s_i和s_j之间的距离，在这里选择的是欧氏距离; $\varphi$用来控制空间相关性的衰减速度，即距离u越大，衰减速度越快; ν为控制平滑程度的参数; K_ν为ν阶的第二类贝塞尔函数.当ν=0.5时，Matérn族相关函数退化为指数相关函数，即$\kappa （u;\varphi ）=\mathrm{e}\mathrm{x}\mathrm{p}（-\varphi u）$; 当ν=3/2时，$\kappa （u;\varphi ）=（1+\varphi u）\mathrm{e}\mathrm{x}\mathrm{p}（-\varphi u）$; 当ν→∞时，Matérn族相关函数退化为高斯过程函数，即$\kappa （u;\varphi ）=\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\varphi }^2{u}^2\right)$.

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

其中${\mathit{Z}}_t={\left(Z\left(s_1，t\right)，\cdots ，Z\left(s_n，t\right)\right)}^{\mathrm{T}}，{\mathit{Y}}_t=\left(Y\left(s_1，t\right)\right.，{\left.\cdots ，Y\left(s_n，t\right)\right)}^{\mathrm{T}}，{\mathit{\epsilon }}_t={\left(\epsilon \left(s_1，t\right)，\cdots ，\epsilon \left(s_n，t\right)\right)}^{\mathrm{T}}，{\mathit{\eta }}_{\mathit{t}}={\left(\eta \left(s_1，t\right)，\cdots ，\eta \left(s_n，t\right)\right)}^{\mathrm{T}}$，且${\mathit{\epsilon }}_t\sim N\left(0，{\sigma }_{\epsilon }^2{\mathit{I}}_n\right)，{\mathit{\eta }}_t\sim N\left(0，{\sigma }_{\eta }^2{\mathit{S}}_{\eta }\right)$.根据分层模型结构，BHAR时空模型可分为三层:第一层表示在给定潜在时空过程和参数的条件下原始数据的分布; 第二层表示参数给定的条件下潜在过程的分布Y_t|Θ; 第三层代表引入的参数或超参数的先验分布.其中第二层中的过程可以添加不同水平的解释^[26]，第一水平表示真正的潜在过程，第二水平刻画了潜在过程的时空随机效应.

2.2 参数估计与预测

BHAR时空模型中的待估参数为$\mathit{\Theta }=\{\mathit{\beta }，\rho ，\left.{\sigma }_{\epsilon }^2，{\sigma }_{\eta }^2，\varphi ，\nu \right\}$，利用MCMC方法进行参数估计.除了$\varphi$和ν以外的其他参数都存在共轭先验分布，即在给定先验分布后β，ρ，${\sigma }_{\epsilon }^2$，${\sigma }_{\eta }^2$都可以求出标准的满条件后验分布，进一步采用Gibbs抽样法进行参数估计.固定ν=0.5，对于参数$\varphi$采用Metropolis-Hastings（MH）算法进行估计.

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

首先介绍空间插值，在未知监测点s₀，由方程（1）可以得到:

其中$\mathit{Y}\left(s_0，t\right)=\rho \mathit{Y}\left(s_0，t-1\right)+\mathit{\eta }\left(s_0，t\right)$.显然，$\mathit{Y}\left(s_0，t\right)$只能由t之前所有时间点的$\mathit{Y}\left(s_0，·\right)$顺序确定，且包括$\mathit{Y}\left(s_0，0\right)$.可以根据初始条件Y₀的先验分布来计算$\mathit{Y}\left(s_0，0\right)$.当然如果指定Y₀为一固定常数，那么$\mathit{Y}\left(s_0，0\right)$也可以被认为是同样的常数^[19]，因此为了简便，Y₀通常选择固定值.

Z（s₀，t）的预测一般基于其后验分布$\pi （\mathit{Z}\left.\left(s_0，t\right)\mid \mathit{Z}\right)$，该后验分布可以通过对联合后验分布积分得到:

上式中Z，Y分别代表已知时间t和监测点s的值.预测值Z（s₀，t）的估计通过MCMC按成分抽样法获得，具体抽样方法如下:

1）从后验分布$\pi （\mathbf{\Theta }，\mathit{Y}\mid \mathit{Z}）$中抽取随机样本${\mathbf{\Theta }}^{（j）}，{\mathit{Y}}^{（j）}$;

2）从后验分布$\pi \left(\mathit{Y}\left(s_0，t\right)\mid {\mathbf{\Theta }}^{（j）}，{\mathit{Y}}^{（j）}\right)$中抽取样本${\mathit{Y}}^{（j）}\left(s_0，t\right)$;

3）从后验分布$\pi \left(\mathit{Z}\left(s_0，t\right)\mid {\mathit{Y}}^{（j）}\left(s_0，t\right)，{\sigma }_{\epsilon }^{2（j）}\right)$中抽取样本${\mathit{Z}}^{（j）}\left(s_0，t\right)$.

在时间维度上的预测过程与空间插值过程类似，对某一站点s（包含现有监测点或任意指定位置作为监测点），做向前一步时间预测，同样可以基于Z（s，T +1）的后验分布根据类似空间插值的MCMC抽样方法实现.与空间插值的不同主要体现在时间维度上的预测需要从具有零均值、方差为σ²η 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两个变量的回归系数β₁和β₃的估计的95% 置信区间包含零点，因此是不显著的.气象变量中只有相对湿度RH的回归系数β₂显著且为正值，这与数据初步分析的结果一致，进一步说明了2—3月北京的相对湿度RH对PM_2.5有显著的正向影响.

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

首先对空间验证集的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₁₀和O₃浓度等.在贝叶斯框架下建模，对模型的不确定性更具有包容性，而且提前给定先验分布可以融合专家知识，提高预测精度.进一步，建立分层模型可以更加清晰地刻画出数据内部的潜在时空过程，增强模型的可解释性.同时，模型的分层结构也使得参数估计和预测等推断过程更加便捷.

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

参考文献