基于INLA-SPDE方法的区域污染物模拟与预测

袁 泽，陈 斌

基于INLA-SPDE方法的区域污染物模拟与预测

袁 泽，陈 斌

(北京大学地球与空间科学学院，北京 100871)

采用传统的空间插值方法对区域污染物进行模拟与预测，针对源数据分布不均，效果一般的问题，提出了采用INLA-SPDE模型来模拟与预测区域污染物的方法。模型的空间分量使用随机偏微分方程表达，时间分量则采用一阶时序自相关模型，同时还包含气象参数等10种协变量，以2019年度京津冀地区日均PM2.5浓度为例，逐月建立了时空模拟与预测模型。实验结果表明，与经典的克里金插值方法相比，在区域污染物分布的模拟上具有更好的效果，尤其在高值污染的预测上精度效果提升明显，同时可得到区域污染风险等级等多种结果。进一步基于模型的预测结果实现了京津冀地区日均PM2.5浓度时空可视化和虚拟仿真系统，为普通民众的出行或政府相关部门决策提供支持，验证了模型的实用性和价值。

细颗粒物PM2.5；贝叶斯时空建模；INLA算法；仿真系统；决策支持

空气污染一直是一个全球性问题[1]，随着我国高人口密度的增加和快速工业化的结合，与经济增长、工业化和城市化相关的严重空气污染事件频繁发生。当下大气颗粒物(particulate matter，PM)，即悬浮在空气当中的固体颗粒或液滴，是最主要的污染物之一。而细颗粒物(空气动力学直径小于等于2.5 μm的颗粒物称为PM2.5)因其粒径更为微小，且具有最为明显的污染效应。同时许多流行病学和毒理学研究表明，PM2.5暴露会导致许多不良健康影响[2]。基于此，政府已经将污染监测和预警工作提到了一个新的高度。

目前我国已经建立了一定规模的地面空气质量监测站点[3]，其监测结果可以有效地反映监测站周边的空气质量状况。但受限于建设成本以及空间占地等问题，地面监测站数量还比较有限，且多位于发达大城市周边，无法覆盖大陆全域。同时近年来随着监测网络的完善，产生了大量的时空数据，这些数据大多具有较高的时间维度分辨率，一般为小时级甚至分钟级，而在空间维度上则相对稀疏，空间数据大多会随着时间发生变化，对联合时间与空间的数据模型的研究日益增多。

在多种时空模型中，基于贝叶斯框架的时空建模方法近年来引起了越来越多学者的关注，其原因是时空数据具有内在的不确定性，而贝叶斯方法从概率的角度出发，能够对时空数据的不确定性进行建模，最终给出概率上最优的结果[4]。在研究空气质量时空分布模型时，SARTO等[5]使用分层贝叶斯模型对意大利中部佩鲁贾镇的PM2.5浓度分布进行了模拟，使用温度、湿度等气象因子作为协变量，同时添加空间效应和时间自回归项，获得了较好的拟合效果。

目前对于贝叶斯时空模型参数的后验分布求解，大都采用马尔可夫链蒙特卡洛方法(Markov chain Monto Carlo，MCMC)，但由于到达马尔可夫链的平稳态往往需要很长时间，在实际应用中，MCMC方法并非作为一种常规分析工具。

2009年，RUE等[6]针对一类特殊的时空隐高斯随机场的贝叶斯推断，提出了一种称为集成嵌套拉普拉斯(integrated nested Laplace approximations，INLA)的快速求解算法。该方法可以求解一类常见的结构化加法回归模型，这类模型符合被一些超参数约束的隐高斯场，同时具有非高斯响应变量。通过使用集成的INLA逼近及其简化版本，该算法可以直接计算出后验分布非常精确的逼近，同时在MCMC算法需要运行数小时或数天的情况下，该算法可在数秒或数分钟内提供更精确的估计。2011年，LINDGREN等[7]提出了一种随机偏微分方程(stochastic partial differential equations，SPDE)算法，其可有效地链接连续高斯随机场和离散的高斯马尔科夫随机场(Gaussian Markov random fields，GMRF)，从而实现对基于Matern相关系数的连续分布的空间效应的快速解算，以进一步加速INLA算法对贝叶斯时空模型后验分布的逼近。

目前，该方法在国外已经开始普及，在与空气质量有关的研究方面，CAMELETTI等[8]针对意大利皮埃蒙特地区的健康状况，采用INLA-SPDE方法获得了区域污染物浓度，并开展了对污染物浓度分布和区域健康风险关系的研究。FIORAVANTI等[9]基于意大利全国410个监测站点2015年的日均PM10浓度数据，使用INLA算法得到了意大利全域时空贝叶斯插值结果。但国内，该方法使用较少，几乎只用于流行病时空分布的分析与模拟[10]。

我国目前的空气质量地面监测站主要分布在城市中心地带，尚未形成空间分布均衡的监测网络，因此难以对大陆区域空气质量状况的时空分布进行合理地评估。因此本文以京津冀地区为例，综合利用研究区域内的气象数据和地理要素数据，建立了区域PM2.5浓度时空分布贝叶斯模型，并采取文献[6-7]提出的INLA-SPDE加速算法，实现了对区域日均PM2.5浓度的快速模拟与预测。并以该预测结果为基础，构建区域PM2.5浓度虚拟仿真决策支持系统，实现京津冀区域PM2.5浓度时空分布图、污染风险等级图的可视化，支持与行政区划图、高程图、人口图等进行叠加显示，从而验证模型的结果实用性。该系统还支持交互放置特定的高污染源于地图上，快速生成新情景下的空气质量空间分布图，对比观察污染物分布的变化，为污染物减排及控制等相关的空间问题提供决策支持。

1 原始数据处理与分析

本文使用的原始PM2.5浓度数据来自京津冀地区236个国控空气质量监测站，数据的原始时间分辨率为1 h，缺失值由站点当前时刻前后3 h污染物浓度数据取均值代替。同时为了保证数据分布的正态性，去除了大于650的极端污染值数据，日均浓度数据则由当日24 h污染物浓度数据求平均值获得，时间跨度从2019年1月到12月。最终得到月度PM2.5浓度日均值箱线图(图1)，图中的橘色虚线和红色虚线分别对应的PM2.5浓度为75 µg/m3和150 µg/m3，这是WHO标准中的第三级和第一级日均PM2.5浓度标准，也是我国分级标准中的污染和重度污染分界线。

图1 京津冀地区2019年月度PM2.5浓度日均值箱线图

在2019年中，PM2.5浓度日均值在0～626 µg/m3之间，年均值为49.9 µg/m3，中位数为37 µg/m3，上下四分位数分别为22.2 µg/m3和60.8 µg/m3。该箱线图与第2节研究一致，呈现出季节性规律：12月、1月、2月的PM2.5浓度日均值中位数大于50 µg/m3，夏季则小于30 µg/m3，其他月份在40 µg/m3。标准偏差也具有类似的规律，1月和2月的标准偏差大于60 µg/m3，夏季小于20 µg/m3，其他月份在35 µg/m3左右。

参与构建贝叶斯模型的协变量在经过反复试验，最终选定了10个时空变量作为模型的预报因子。具体的数据信息见表1。为了避免数据自身的量级对研究的影响，所有的数据均进行了标准化处理，即全部转换为了均值为0，标准差为1的正态分布。

表1 贝叶斯时空模型使用的协变量

使用的6种气象因子全部来自由美国国家环境中心(National Centers for Environment-Mental Prediction，NCEP)提供的6 h气象再分析资料。该数据提供了每日0点、6点、12点、18点的空间分辨率为0.25度的格网化气象要素，下载后以日为尺度重新计算了均值。气溶胶光学厚度(aerosol optical depth，AOD)数据则来自欧空局(Copernicus Atmosphere Monitoring Service，CAMS)气象再分析资料，是在空间分辨率为10 km的格网上对波长为550 nm的AOD数值模拟得到的结果。

京津冀地区各级道路数据获取来自OSM (OpenStreetMap Project)，从Geofabrik网站提供的服务下载得到。本文选取了OSM数据集中标记为highway的线状地物中的前四级标签作为基础道路数据，包括了京津冀地区的国道、省道以及县道。最终投入模型的路网密度数据则是通过将京津冀地区划分为10 km×10 km的格网，计算每个格网内道路数据的总长度(km)来获得的。海拔数据则来自USGS提供的GTOPO30全球数字高程模型，分辨率为30弧秒(约1 km)。人口数据采用WorldPop datasets提供的中国大陆1 km分辨率的人口密度数据。

2 京津冀地区PM2.5日均浓度时空贝叶斯建模

2.1 贝叶斯时空建模方案设计

对2019年全年的京津冀地区日均PM2.5浓度进行建模，按月份建立12个模型，具体模型如下：

图2 区域SPDE空间效应三角网示意图

2.2 预测结果精度评估

首先检验模型的站点预报精度，表2为12个月度模型在训练集和测试集上3次交叉验证后的平均结果，这些指标数值均在转换为原始数据尺度后计算获得。总体来说，12个模型在训练集和测试集上的预测均表现良好，而且效果相当，说明模型并不存在过拟合。RMSE和MAE的值在冬季要高于夏季，这与冬季PM2.5浓度值整体偏高有关。

表2 模型交叉验证结果(µg/m3)

观测值和预测值的关联性如图3所示，为了避免过多的散点图，本文将12个月度模型按季度展示其真实值与预测值的一致性。图中颜色越浅的部分代表处于此浓度的观测数据的数量越多，红色实线对应1﹕1参考线。无论是在训练集还是测试集上，点几乎都均匀地分布于参考线两侧，同时浅色区域几乎围绕着参考线两侧最近处分布，可以认为模型在4个季度对于大部分样本数据预测效果表现良好。不过从图3中也可以看到，在测试集上，当观测值大于250 µg/m3时，出现了较多偏小的预测值。

图3 模型预测值与真实值对比图((a)训练集；(b)测试集)

图4以污染最严重的1月和相对最轻的8月为例，绘制了区域内236个站点月度的日均观测值与预测置信区间的时序关联图，图中黄色区域代表预测浓度值的95%置信区间。该时序图表明京津冀地区PM2.5浓度日均值的时间变异性可以很好地被模型捕捉，观测值浓度越高，预测值的置信区间范围越大。但真实值总是落在预测的置信区间内部，证明了模型对于时序特征的复现能力。

2.3 预测结果对比分析与可视化

本文选取具有代表性的日期，将模型的预测结果与传统的空间插值通用方法——克里金插值方法进行比较。即在相同的训练集上建立INLA-SPDE模型和进行普通克里金空间插值，最终将测试集上的插值结果与INLA-SPDE模型在测试集上的预报结果进行对比，其中贝叶斯模型的预测结果均采用模型的预测后验均值数据。

图4 PM2.5日均浓度时序验证图((a)1月；(b)8月)

具体选取的实验数据从2019年1月7日—18日，这一时段内京津冀地区经历了一次污染从严重到消散的过程，具有一定的代表性。最终在完全相同的测试集上试验得到的RMSE数值见表3，可以发现使用的INLA-SPDE模型在空间上的预测效果要明显优于传统的克里金插值。根据表3的数据绘制得到图5，可以发现2种模型预测精度的差距在污染浓度越大时越明显，证明了本文模型在高污染值上的预测能力远强于传统的克里金插值方法。从PM2.5浓度预测的实际应用出发，对捕捉严重污染出现时机和位置的需求更大，因此本文模型具有较强的实用价值与现实意义。

图6以京津冀地区5 km×5 km的格网为基础，展示了克里金插值方法与本文模型的对比结果。选取了1月13日和15日分别为污染严重和开始消散的代表日期，值得注意的是左右2幅图颜色分级尺度不同。可以发现，本文模型对于空间效应的模拟与克里金插值相似，大尺度上的空间分布两者保持着高度一致性。该模型对于污染浓度偏高的值更为敏感，能更为准确地反映严重污染地带的分布，2幅图中极值亮点均出现在左图，相比与传统的克里金方法，INLA-SPDE模型对于小尺度上的空间变异性复现更为精确。

表3 INLA模型与克里金插值预测RMSE结果(µg/m3)

图5 INLA模型与克里金插值预测RMSE对比图

图6 INLA模型与克里金插值预测结果对比图((a)1月13日；(b)1月15日)

最后相较于传统的克里金空间插值，因为该模型的预测结果本质上是所有待预测格网点处污染浓度值的后验概率密度，因此还可以给出区域污染风险等级的空间分布，这对于人们日常的出行安排和政府部门空间决策提供了支持。京津冀地区典型日期的污染风险如图7所示。

图7 京津冀地区1月12日(左)和11月22日(右)污染风险等级空间分布图

3 区域PM2.5浓度虚拟仿真决策支持原型系统

本节基于第2节建立的模型构建京津冀地区PM2.5浓度虚拟仿真决策支持原型系统。主要目标是验证该模型的实用性，实现任意时间模型结果的区域可视化，呈现时空预测结果中蕴含的有效信息。并通过交互控制实现人工模拟大气环境的改变来观测区域效应，从而为普通民众的出行或政府部门相关人员的决策提供支持。

3.1 系统模块设计

原型系统包含3个基础模块，分别为模型层、数据层和交互层，如图8所示。

图8 原型系统工作模块图

将数据层中监测站的实时数据与有关地理要素数据合并成数据表，并存储到Postgres数据库中。气象数据为格网化数据，以GeoTIFF格式存储。

模型层输入为选定日期的测站历史数据，调用相关月度的INLA-SPDE模型，将结果数据存储到数据层后，再对区域PM2.5浓度预报的空间结果可视化。

交互模块通过Web前端界面加入相应交互按钮，支持人工模拟新污染源的添加与修改，求解得到交互后的模拟空气质量时空分布结果并可视化。

3.2 模型的原型系统验证

用于验证的原型系统基本界面主要包括，交互和结果展示2个区域。交互区如图9左侧边栏所示，包含日期选择、想要修改污染值的点位经纬度信息、污染辐射范围和浓度值的录入，其中位置的经纬度信息也可以通过地图区的绘制功能交互点选实现。

图9为主区域为京津冀地区PM2.5浓度时空分布结果地图，系统可根据所选可视化区域PM2.5浓度空间分布模拟结果。地图区通过左侧的功能栏，支持显示污染风险等级图和空间分布结果图与多种验证数据图的叠加。该结果展示区还包括若干功能，如上方显示的是PM2.5浓度超标的测站个数、预测PM2.5浓度超标的区域占京津冀地区的比例和重污染风险浓度在75%以上的区域比例，对于整个区域这一时刻的污染情况可以有直观且清晰的感知。右侧则显示了污染具体的行政区划分布情况，包括站点PM2.5浓度排序和以市为基本单位计算的高污染风险区域所占的比例，从而为特定城市的相关人群提供有价值的参考信息。

界面下方为交互模拟结果展示，根据交互区输入的相关信息得到新的空间分布结果并可视化，与主地图区对比，效果如图10所示。图中更改浓度值的点位为原结果图中的蓝色圆圈，将此处的污染值向下修正到了50 µg/m3，并设定辐射范围为100 km。在变动结果图中，可以清楚地看到蓝色圆圈周围的监测站值全部降低到75 µg/m3以下(图中的红色代表污染值在75 µg/m3以上的站点，绿色代表在75 µg/m3以下的站点)。从区域分布上可知，整个北京以北的地区污染状况显著改善，为了量化这种污染改善的程度，系统右下角还展示了不同污染等级的区域格网数目变动百分比，可以看到重度污染以上的区域数目大幅减少，而转为良好等级的区域数目则明显上升。

图9 原型系统基础界面图

(a) (b)

因此，通过该交互控制功能，可以清楚地展现在特定区域减排或设立新的污染源的影响程度，为相关政策的针对性制定和部分建筑的空间选址问题提供了重要的参考，真正实现了空间决策与支持功能，进一步验证了模型的实用价值。

4 结论与反思

本文提出了基于INLA-SPDE模型的区域污染物模拟与预测方法，实现了京津冀地区PM2.5浓度日均值的时空建模与预测，验证了模型预测精度，得到了京津冀地区5 km×5 km格网的日均PM2.5浓度时空分布预报图和污染风险等级图等。结果表明，INLA-SPDE方法在区域污染模拟与预测上相比传统空间插值方法具有更好的效果，尤其对于高污染值的预测更为精确，并能给出区域污染风险等级预测结果，在区域污染的预警和防控上具有重要的现实意义。基于该模型构建了京津冀地区PM2.5浓度预报与控制决策支持原型系统，解决了区域级别的PM2.5浓度时空可视化的问题。通过将人工模拟得到的仿真结果与原始结果对比，为人类的空间决策行为提供了有效地支持，验证了该模型的实用性与价值。

[1] Report from Intergovernmental Panel on Climate Change (IPCC). Climate change 2007: the physical science basis[EB/OL]. [2021-04-01]. http://www.doc88.com/p- 293361412559.html.

[2] APTE J S, BRAUER M, COHEN A J, et al. Ambient PM2.5reduces global and regional life expectancy[J]. Environmental Science & Technology Letters, 2018, 5(9): 546-551.

[3] CAI S Y, MA Q, WANG S X, et al. Impact of air pollution control policies on future PM2.5concentrations and their source contributions in China[J]. Journal of Environmental Management, 2018, 227: 124-133.

[4] GAUDARD M, KARSON M, LINDER E, et al. Bayesian spatial prediction[J]. Environmental and Ecological Statistics, 1999, 6(2): 147-171.

[5] SARTO S D, RANALLI M G, BAKAR K S, et al. Bayesian spatiotemporal modeling of urban air pollution dynamics[M]// Topics on Methodological and Applied Statistical Inference. Cham: Springer International Publishing, 2016: 95-103.

[6] RUE H, MARTINO S, CHOPIN N. Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations[J]. Journal of the Royal Statistical Society: Series B: Statistical Methodology, 2009, 71(2): 319-392.

[7] LINDGREN F, RUE H, LINDSTRÖM J. An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach[J]. Journal of the Royal Statistical Society: Series B: Statistical Methodology, 2011, 73(4): 423-498.

[8] CAMELETTI M, GÓMEZ-RUBIO V, BLANGIARDO M. Bayesian modelling for spatially misaligned health and air pollution data through the INLA-SPDE approach[J]. Spatial Statistics, 2019, 31: 100353.

[9] FIORAVANTI G, MARTINO S, CAMELETTI M, et al. Spatio-temporal modelling of PM10daily concentrations in Italy using the SPDE approach[J]. Atmospheric Environment, 2021, 248: 118192.

[10] 丁书珍, 张辉国, 胡锡健, 等. 利用R-INLA方法分析宏观因素对艾滋病疫情的影响[J]. 中国艾滋病性病, 2018, 24(12): 1192-1196.

DING S Z, ZHANG H G, HU X J, et al. Influence of macroscopic factors on HIV/AIDS epidemic based on R-INLA method[J]. Chinese Journal of AIDS & STD, 2018, 24(12): 1192-1196 (in Chinese).

[11] BLANGIARDO M, CAMELETTI M. Spatial and spatio-temporal Bayesian models with R-INLA[M]. Chichester: John Wiley & Sons, Ltd, 2015: 268-283.

[12] SIMPSON D, RUE H, RIEBLER A, et al. Penalising model component complexity: a principled, practical approach to constructing priors[EB/OL]. [2021-04-04]. https://arxiv.org/ pdf/1403.4630.pdf.

[13] SØRBYE S H, RUE H. Penalised complexity priors for stationary autoregressive processes[J]. Journal of Time Series Analysis, 2017, 38(6): 923-935.

Simulation and prediction of regional pollutants based on INLA-SPDE method

YUAN Ze, CHEN Bin

(School of Earth and Space Sciences, Peking University, Beijing 100871, China)

The simulation and prediction of regional pollutants generally use the traditional spatial interpolation method, which cannot obtain accurate results when the source data is not uniformly distributed. To address these problems, a method for simulation and prediction of regional pollutants based on the INLA-SPDE model was proposed. The interpolation model was based on a Bayesian hierarchical model where the spatial-component was represented through the stochastic partial differential equation (SPDE) approach, with a lag-1 temporal autoregressive component (AR1). In addition, the model included 10 spatial and spatio-temporal predictors such as meteorological variables. By building 12 models for each month with the integrated nested Laplace approximation (INLA), this research realized the spatio-temporal simulation and prediction of PM2.5concentration at daily resolution in the Beijing-Tianjin-Hebei region in 2019. Experiments show that compared with traditional Kriging interpolation methods, the proposed model can yield a better prediction of air pollutants at regional scale. Particularly, the prediction accuracy of high-value pollutants was improved significantly, and air pollutants exceedance probabilities can also be generated. Furthermore, a system for regional PM2.5concentration simulation and decision support was established, the system can provide support for the travel of ordinary people or the decision-making of government officials, and verify the practicability and value of the proposed model.

PM2.5; spatio-temporal Bayesian hierarchical model; integrated nested Laplace approximation; simulation system; decision support

18 May，2021；

TP 391

10.11996/JG.j.2095-302X.2022010125

A

2095-302X(2022)01-0125-08

2021-05-18；

2021-08-21

21 August，2021

袁 泽(1992–)，男，硕士研究生。主要研究方向为虚拟地理环境。E-mail：180121089@pku.edu.cn

YUAN Ze (1992–), master student. His main research interest covers virtual geographic environments. E-mail：180121089@pku.edu.cn

陈 斌(1973–)，男，教授，博士。主要研究方向为虚拟地理环境。E-mail：gischen@pku.edu.cn

CHEN Bin (1973–), professor, Ph.D. His main research interest covers virtual geographic environments.E-mail：gischen@pku.edu.cn