基于空间插值算法的2014—2018年河南省PM2.5时空分布及演变特征研究

张馨文

摘 要：本文选取河南省2014—2018年监测网络中269个监测站点的空气质量资料，在地理统计学和GIS的支持下，分别用四种不同的空间插值方法对各年PM2.5年均浓度进行插值，通过模拟精度的对比以及制图分析发现：泛克里金和普通克里金插值精度较高;而可视化方面普通克里金具有最佳效果。将普通克里金插值结果与前期探索性数据进行整合，得出结论：PM2.5浓度大致从由中北部向四周递减变为由北部向南部递减，在东西方向上呈现中部高两侧低的特点，2014年以来河南省大气质量整体在改善。

关键词：PM2.5;空间插值;时空分布

Abstract： In this paper， selected the air quality data of 269 monitoring stations in the monitoring network of Henan Province from 2014 to 2018， and the annual average concentration of PM2.5 was interpolated by four different spatial interpolation methods with the support of geostatistics and GIS. Through the comparison of simulation precision and cartographic analysis， it was found that Ordinary Kriging and Universal Kriging have the highest accuracy， while the Ordinary Kriging has the best effect in the visual. Combining the Ordinary Kriging interpolation with the ESDA， the result showed that the concentration of PM2.5 is changing from the pattern of decline from north-central to the sidesinto the north to the south， whilein the east-west direction the middle is higher than the sides， and the overall air quality has improved in Henan Province since 2014.

Keywords： PM2.5;spatial interpolation;spatial and temporal distribution

空氣质量与每个人的生活息息相关。但是，随着我国城市化和工业化进程不断加快，大范围雾霾天气频频发生，严重影响了居民的健康水平、经济与环境的可持续发展[1]。作为重点大气污染源[2]，PM2.5一方面粒径细小，能较长时间悬浮于空气中，显著减弱大气能见度;另一方面，其化学组成复杂，富含大量有害物质，极易损害人体健康。党的十九大报告中指出“污染防治”是我国全面建成小康社会实现第一个一百年奋斗目标必须打好的“三大攻坚战”之一，以防治PM2.5为重点的大气污染防治工作在我国必然被着力推进。然而，近年来，河南省空气质量形势不容乐观，尤其到18年末，多地持续出现雾霾红色预警，17个地级市中就有12个频频登上全国城市空气质量排行的末位。由此，结合区域特征，利用更适宜的方法探讨河南省PM2.5的时空分布规律，成为优化大气防治措施的关键。

国内有关PM2.5时空分布的研究区域主体往往以全国为重点[3-4]，而区域尺度上，多数研究基本关注京津冀、长三角等污染核心地带[5-6]，尚未发现以河南省为区域主体的研究。同时，研究数据主要来自遥感反演和地面监测[7]。由于反演所用的光学气溶胶厚度值（AOD）时间序列上不稳定，只能通过分辨率较粗的气象数据修正[8]，反演模型精度难以提升。针对时间尺度不断增长的需求，越来越多的学者选择将研究数据主体转向具有稳定时序的监测数据。

本文选择河南省具有完整监测数据的2014—2018年为研究时段，基于392个监测站实时监测的空气质量数据，利用4种不同空间插值方法探索出PM2.5模拟分布的最佳结果，进而获得2014—2018年来河南省PM2.5质量浓度时空分布格局和演变规律，为相关研究和河南省大气污染防治政策优化提供科学的参考依据。

1 研究数据与方法

1.1 研究区概况

河南省位于中国中部，界于北纬31°23'—36°22'，东经10°21'—116°39'，东接安徽、山东，北接河北、山西，西连陕西，南临湖北，总面积16.7万km2，大部分地区属北亚热带向暖温带过渡的大陆性季风气候，南北过渡性明显。经济方面，河南省地处沿海开放地区与中西部地区的结合部，是中国经济活动承东启西、通南达北的中间地带。截至2018年，河南省区域总人口和年GDP总量分别位于全国省份排行的第三和第四名。

1.2 数据来源与处理

本实验的数据主要包括：①从真气网（https：//www.zq12369.com/）上爬取的河南省17个地级市共包括国控点和省控点在内的392个监测站2014年1月1日至2018年12月31日监测的空气质量数据;②来自2014—2018年《河南统计年鉴》的各县（市）年末常住人口、GDP数据;③由Esri公司提供的ChinaData矢量数据集（包括省级行政区、地级市、县城等）。

针对各站点监测的PM2.5逐时浓度数据，剔除了数据不完整的123个站点，余下269个监测站作为后续插值的采样点（如图1所示），逐步求得各采样点24h日平均质量浓度和年平均质量浓度。

1.3 探索性空间数据分析

探索性空间数据分析（Exploratory Spatial Data Analysis，ESDA）不仅能显示空间数据的内在规律，而且是后续插值方法相关参数和模型选择的依据。

1.3.1 空间自相关分析。空间插值建立在空间相关性分析的基础上，若分析检测的是同一属性变量，称为空间自相关，反映的是一个区域单元上的某种地理现象或某一属性值与邻近区域单元上同一现象或属性值的相关程度[9]。常用Moran’s I系数来描述空间邻近区域单元的相似程度，计算结果在-1～1，大于0呈正相关关系，属性值高或低的地区都有聚集现象;若结果小于0，呈负相关关系，数据空间分布呈现高低间隔分布的状态;绝对值越接近1说明相关性越强，若结果趋近0，则相邻区域单元相关性低，在研究区呈随机分布状态。

1.3.2 同步统计检验。同步统计检验可以完成检验数据分布、寻找离群值以及进行全局趋势分析等一系列任务。直方图可以对大量采样点进行整理加工，找出其数据分布的规律和形态，以便对其总体分布特征进行推断。QQ图是一种用于用于检验数据分布的相似性的统计图，被检验的采样点数据越符合所指定的分布，越聚集在一条直线上。趋势分析图用X、Y平面表示采样点位置，用Z轴表示数据属性值，通过观察Z轴数据值投影到X、Y平面散点曲线的趋势，判断数据是否存在全局趋势以及其趋势形态，若存在，可使用某种确定性内插方法生成表面，或通过范克里金或协同克里金模型移除全局趋势。

1.4 插值方法

采样点的数量有限，空间分布也不够均衡，要获得区域内整体数据的分布认知，必然要进行空间插值。空间插值的方法多样，每种方法都有其特点和适用性。用多种方法对PM2.5浓度插值，在保证预测结果适用性的同时还可对比出其误差规律[10]，并探索与PM2.5有较强相关性的因素，为河南省PM2.5浓度预测和调控工作提供参考。

空间插值的公式为：

1.4.1 确定性空间数据插值。反距离加权法（Inverse Distanceweight，IDW）认为，预测点的数据值是其邻近采样点综合贡献的结果，贡献程度与距离成反比。确定权重的公式为：

1.4.2 克里金（Kriging）插值。克里金插值法不仅考虑了预测点与邻近采样点的空间距离关系，还考虑了其位置关系。该方法建立在区域化变量理论分析的基础上，以半变异函数作为工具，对有限区域内的区域化变量取值进行线性无偏最优估计。

半变异函数是一个关于数据点的半变异值与数据点间距离的函数，计

1.4.2.2 泛克里金（Universal Kriging，UK）。对不满足平稳性假设的数据，UK可以对空间坐标表示的变量按照自相关进行回归分析，即线性无偏最优估计时考虑其趋势，其模型可表示为：

1.4.2.3 协同克里金法（Collaborative Kriging，CK）。CK是OK区域化变量单个向多个的拓展，通过计算属性间的半方差函数和交叉半方差函数利用变量间的空间连续性、相关关系[11]，较为复杂。

1.5 插值检验

1.5.1 交叉验证。本实验选择ArcMap10.4软件Geostatistical Analyst模块中的交叉验证工具调整插值方法的参数，根据平均误差最接近于零，均方根预测误差最小，平均标准误差最接近于均方根预测误差，标准均方根预测误差最接近1的标准判断出克里金插值最合适的半变异函数模型。

1.5.2 子集验证。实验通过子集验证的方式比较四种不同插值方法的精度。先从研究区的269个采样点中随机选择20个作为检验站点，剩余249个采样点作为分析插值数据重新进行插值，各点分布如图1所示。对不同插值结果的精度選择绝对平均绝对误差（Mean Absolute Error，MAE）和均方根误差（Root-Mean-Square Error，RMSE）进行比较。二者越小插值精度越高;MAE相等时，RMSE越小越好。

2 插值应用过程

首先观察ESDA结果（见图2），直方图显示年均PM2.5浓度的中位数基本都与平均值较为接近，说明数据接近正态分布;而正态QQ图则可以显示出数据总体是倾向于正态分布的，插值前无需进行相应的数据转换。通过趋势分析可以发现，河南省各市的年均PM2.5浓度具有一定的空间相关性，其基本趋势为北高南低，东西向为中部高两侧低。

OK不考虑移除趋势，其精度主要受半变异函数的影响。本文在保证变异函数的“步长个数[×]步长大小=研究区域长度的一半”的情况下，分别选择球面、高斯、指数、三角函数四种半变异函数模型进行插值，对比交叉验证结果后，除2015年选择三角函数外，其余各年都选择高斯函数做OK的半变异函数模型。UK综合了探索性数据分析结果去除数据的趋势特征，选择趋势的移除阶数为1，确定高斯函数作为适合表面的核函数，其半变异函数建模过程与普通克里金保持一致。CK考虑雾霾与经济发展的紧密关联，分析二者的正交协方差函数云图后，选择各市县的年末GDP作为增加的区域化变量，半变异函数建模过程依旧与OK保持一致。IDW需要确定的参数极少，这里经调试后统一设置指数为2，权重取自15个相邻要素。

3 实验结果分析

3.1 2014—2018年河南PM2.5浓度的空间集聚趋势

实验利用GeoDa软件中的Arc Distance距离算法生成空间权重矩阵计算各年的Moran’s I，结果如图3所示。

从图3可以看出，PM2.5的分布在2014和2017年具有一定的空间自相关性，存在集聚现象;在2016和2018年存在空间相关性，但相关性不高;而2015年几乎没有自相关性。观测Moran散点图可以看出浓度值的聚集类型。相关性最强的是2014年，其散点图中显示，大部分测点在一、三象限，说明在2014年河南省内PM2.5多呈冷、热点分明的分布模式。

3.2 插值预测结果分析

3.2.1 不同方法预测结果比较。对插值结果制图得到不同方法模拟的PM2.5浓度空间分布图像（图略）。IDW预测结果非常平滑，但能明显看出PM2.5浓度的极值只出现在采样点处，受采样点分布的影响较大，尤其是在河南省中北部地区产生很多孤立异常的图斑，与实际分布情况有差距。OK所得浓度结果过渡性好，比较符合实际，可以体现出空间变化趋势，制图的视觉效果也很理想，但其存在预测的最高值偏低，而最低值偏高的问题。UK也可以反映出PM2.5的空间变化趋势，但成图存在块状效果。CK明显受到GDP变量的影响，使得模拟的PM2.5浓度在部分地区空间变异过大。GDP极值处常有孤立异常的图斑存在，与实际情况有区别，这也说明实际的PM2.5浓度分布受众多因素影响，仅考虑GDP作为协同变量还不够。可视化上OK最优，可以准确地表现PM2.5的空间分布特征与变化规律。

3.2.2 2014—2018年河南PM2.5的空间分布及动态变化。4种插值方法均大致显示出河南省的PM2.5浓度分布格局与变化规律，综合趋势面分析结果：近年来，PM2.5浓度大致从由中北部向四周递减变为由北部向南部递减，在东西方向上依然呈现中部高两侧低的特点。此外，插值预测到的PM2.5浓度的年均值的域值逐年降低，说明2014年以来河南省大气质量整体上是逐年改善的。

3.3 插值精度分析

由于IDW没有对预测误差产生估计，并不适用于交叉验证，实验转而利用子集验证的方式进行误差比对，对比结果如表1所示。首先确定实测值与各种插值方法的估算值无显著差异。

综合两指标的计算值及其排序情况得出结论：河南省2014和2016年的年均PM2.5浓度用UK插值所得结果精度最高，2015、2017年和2018年用OK插值所得结果精度最高。虽然不同插值方法插值精度没有特别显著的差别，但克里金法要普遍优于确定性空间插值方法。

4 讨论与总结

本文选取了2014—2018年河南省17个地级市269个监测站点的空气质量资料，以此为采样点检验数据分布特征，利用IDW、OK、UK、CK共4种方法分别进行PM2.5年均浓度的插值，制作了各方法模拟各年PM2.5空间分布的图像，通过交叉验证确定高斯函数为OK拟合区域化变量的最优函数，接着进行子集验证，结合ESDA与OK的结果，最后得出结论：①2014和2016年UK的精度最高，2015、2017年和2018年OK精度最高，视觉上OK最能反映空间变化规律;②PM2.5年均浓度大致从由中北部向四周递减变为由北部向南部递减，在东西方向上依然呈现出中部高两侧低的特点;③插值预测到的PM2.5浓度的年均值的域值在逐年降低，说明2014年以来河南省大气质量整体上是逐年改善的。

此外，CK的应用过程表明，近年来，经济的发展与空气质量的冲突在逐渐减弱，体现了河南省经济结构正逐步向与环境协调化发展转型的特征。深入研究PM2.5的影响因素，如风力风向、湿度、汽车尾气排放，从而确定最适宜研究区的多协同变量的克里金模型也是今后研究的主要方向。

参考文献：

[1]Lelieveld J， Evans J S， Fnais M， et al. The contribution of outdoor air pollution sources to premature mortality on a global scale[J]. Nature，2015（7569）：367-371.

[2]盛黎，繆育聪，徐权周，等.对话黄顺祥：走近中国大气污染防治[J].科学通报，2018（16）：1515-1518.

[3]王振波，方创琳，许光，等.2014年中国城市PM2.5浓度的时空变化规律[J].地理学报，2015（11）：1720-1734.

[4]李名升，任晓霞，于洋，等.中国大陆城市PM污染时空分布规律[J].中国环境科学，2016（3）：641-650.

[5]杨兴川，赵文吉，熊秋林，等.2016年京津冀地区PM2.5时空分布特征及其与气象因素的关系[J].生态环境学报，2017（10）：1747-1754.

[6]戴昭鑫，张云芝，胡云锋，等.基于地面监测数据的2013—2015年长三角地区PM2.5时空特征[J].长江流域资源与环境，2016（5）：1004-8227.

[7]熊欢欢，梁龙武，曾赠，等.中国城市PM2.5时空分布的动态比较分析[J].资源科学，2017（1）：136-146.

[8]张云芝.基于地面监测和遥感反演的京津冀地区PM2.5时空分布研究[D].呼和浩特：内蒙古师范大学，2016.

[9]刘湘南，王平，关丽，卢浩，张春晓.GIS空间分析[M].3版.北京：科学出版社，2008.

[10]朱求安，张万昌，余钧辉.基于GIS的空间插值方法研究[J].江西师范大学学报（自然科学版），2004（2）：183-188.

[11]章清，张海涛，郭龙，等.基于主成分分析的协同克里格插值模型对土壤铜含量的空间分布预测[J].华中农业大学学报，2016（1）：60-68.