基于时空指示克里格的PM2.5不确定性分布

梅 杨,张文婷,杨 勇*,赵 玉,李露露 (.宁波市镇海规划勘测设计研究院,浙江 宁波 35200；2.华中农业大学资源与环境学院,湖北 武汉 430070；3.农业部长江中下游耕地保育重点实验室,湖北 武汉 430070)

近年来,随着国民经济的持续快速发展,我国PM2.5污染有常态化趋势[1],研究表明[2-5],若长时间暴露于高浓度 PM2.5中,不仅会引发各种呼吸系统病变,同时还会增加癌症的发病率,严重威胁人类身体健康.当前,对于 PM2.5的研究,主要有以下4个方面：(1)从时间角度出发,运用相关方法研究单个监测站点的PM2.5数据在时间维度上的变化特性,如李梓铭等[6]通过研究北京城区 PM2.5在不同时间尺度上的周期性得出,北京城区 PM2.5浓度存在多个明显的周期性变化;(2)从空间角度出发,运用相关方法研究不同监测站点的 PM2.5数据在空间上的变化特性,如 Gehrig等[7]研究了瑞士1998～2001年期间 PM2.5与 PM10变化特征,发现PM2.5和PM10存在高相关性,且冬季PM2.5浓度高于夏季;庄欣等[8]研究了珠三角大气污染的空间分布特征,结果表明在珠三角地区大气污染存在明显的区域性特征;(3)从 PM2.5自身角度出发,研究PM2.5的化学组成成分及来源,如周甜等[9]研究了华北平原夏季 PM2.5的化学组成成分及来源,结果指出PM2.5中二次无机离子的含量达到60%,其来源也受到工业源和尘源的影响;(4)从时空的角度出发,运用时空地统计学理论,研究PM2.5数据在时空维度的变化特征,如梅杨等[10]研究了山东省PM2.5的时空演变特征,结果发现PM2.5具有明显的时间污染特征和空间污染特征;Christakos等[11]在时空克里格的基础上,提出了一种新的时空降维预测方法,并以山东省 PM2.5为研究对象进行预测并对结果进行分析,证明了该方法优于时空克里格.此外,其他学者也从上述4个方面对PM2.5的时空浓度分布特征进行研究[12-15],但大多数均着力于对 PM2.5的浓度值进行分析,属于确定性分析范畴,而对于 PM2.5达到某种空气质量等级的概率性分析,国内外乏善可陈;同时,PM2.5作为大气污染物的一种,其变化受温度、湿度、风力等气象因素影响较大,若仅对PM2.5的浓度值做确定性分析,不仅会使得分析结果不够严谨,同时也不利于环保部门对空气污染来源的探索以及制定完善的空气污染防治措施.基于此,本研究拟以山东省 2014年PM2.5日均质量监测浓度为数据源,根据2012年我国环境保护部门颁布的《环境空气质量标准》[16]和环境保护环境监测司公布的空气质量指数对人体健康的影响状况[17]为阈值,扩展空间领域的指示克里格法,提出时空指示克里格法,并对山东省2014年PM2.5日均质量浓度进行不确定性分析.

1 材料与方法

1.1 数据来源及预处理

本研究 PM2.5数据来源于中国空气质量在线监测分析平台(http：//www.aqistu dy.cn/),该平台在山东省分布有99个监测点位(图1),每小时监测一次,并同步进行数据信息更新.数据监测时间为2014年1月1日0：00时至2014年12月31日23：00.对于获取的监测数据,由于数据更新周期短,数据量极大,不利于对数据进行时空建模分析,因此对各站点的监测数据按天取均值处理.图 2为各监测站点PM2.5日均质量浓度基本统计信息.

图1 山东省空气质量监测点位分布Fig.1 Spatial distribution of air quality monitoring stations in Shandong Province

图2 山东省PM2.5质量浓度基本统计特征Fig.2 Statistical characteristics of the PM2.5 concentrations for Shandong Province

1.2 多阈值时空指示克里格方法

作为一种大气污染物,PM2.5在不同尺度下其变量之间的自相关程度相差较大,且随着监测站点时间与空间距离的增大,其时空变异函数的随机成分也逐渐增大,若仅对PM2.5运用单一阈值的时空指示克里格(STIK)进行分析,即只分析 PM2.5在某一阈值下的时空结构特征,不仅掩盖了其在大尺度下的结构特征,同时也不利于分析其在小尺度下的宏观结构.因此,本研究采用多阈值时空指示克里格法进行分析.

1.2.1 多阈值时空指示变量与时空指示经验变异函数 多阈值时空指示克里格(MTSTIK),是指在多个时空阈值下,分别运用时空指示克里格法对时空变量进行时空指示预测分析[18].在MTSTIK中,定义D为空间域,T为时间域,若在研究区域存在指示值(阈值)Zc(c=1,2,…,n),则对于时空域D×T上的每一点Z(s,t)(s=(x,y)∈D,t∈T),其时空指示变量值I c (s , t; Z c )为：

对于各阈值下的时空指示变量,其时空指示经验变异函数与时空普通克里格经验变异函数

[10]类似,即在本征假设的条件下,定义时空指示变异函数为：

式中：hS,hT分别为空间和时间间隔变量,N(hS,hT)为样点中符合所定义空间和时间间隔的点对数.

1.2.2 时空理论变异函数模型 时空理论变异函数模型通常分为时空分离模型和时空非分离模型[19],两者的差别在于时空分离模型中存在独立的时间和空间部分,而时空非分离模型中则将时间与空间统一考虑.在本次研究中,PM2.5质量浓度变化受时间和地域的共同影响,若将 PM2.5质量浓度变化分别在时间维度和空间维度进行分解,不仅与实际情况不符,其分析结果也不具备科学性和严谨性.因此,本研究采用 Cressie等[20]提出的时空协方差函数模型(式 3)对各阈值下的时空经验指示变异函数进行拟合.

1.2.3 时空指示克里格预测 对于各阈值下的时空指示变量,其时空指示克里格插值理论基础同普通指示克里格插值算法类似,即构造克里格权重计算矩阵A⋅λ= b,然后通过时空变异函数计算A、b矩阵,求取时空权重,进而利用求取待估点小于阈值Zc的概率[19].

1.3 PM2.5时空分布的不确定性分析

根据各阈值对应的时空预测结果,绘制相应的时空立方体.并统计各空间点位对应的时空预测结果,分析各时空立方体结构特征以及计算各空间点位对应的年均概率值、月均概率值和部分城市城区范围内日均概率值,实现对山东省PM2.5时空分布的不确定性分析.

2 结果与分析

2.1 时空指示克里格阈值

表1 PM2.5浓度值对应的空气质量级别及对健康的影响状况Table 1 Air quality level and health effects to the concentration of PM2.5

MTSTIK关键在于如何确定多个既能反映研究区整体的分布特征,又能揭示局部时空结构特征的时空阈值.在本文中,根据我国环境保护部2012年颁布的《环境空气质量标准》中规定的日均 PM2.5浓度标准值和我国环境保护部环境监测司2012年公布的空气质量指数对人体健康的影响情况,依次选择 35,75,150和 250μg/m3作为本次研究的时空指示值(表1).

2.2 时空指示经验变异函数与模型拟合

按照公式(1)和对应的时空指示值,将山东省PM2.5质量浓度指示化.同时,设置空间步长为5km,最大空间计算距离为100km,时间步长为1d,最大时间计算距离为 14d,运用公式(2)计算各阈值下的时空指示经验变异函数,并根据公式(3)对各时空模型进行拟合(表2,图3).

表2 4种时空指示理论变异函数拟合参数Table 2 Parameters in 4 different spatio-temporal semivariogram models

由表2和图3可知,在空间距离大于100km时,各时空指示克里格变异函数值仍有增大趋势,说明山东省PM2.5浓度变化的空间自相关范围大于 100km;在时间距离在 3d时,各变异函数值已趋于稳定,说明山东省 PM2.5浓度变化的时间自相关范围为3d左右.

图3 4种时空指示理论模型拟合Fig.3 Spatio-temporal indicator semivariograms under 4thresholds

2.3 多阈值时空指示克里格预测立方体

设置MTSTIK预测格网大小为2km×2km×1d,运用时空指示克里格方法,对各阈值下的PM2.5指示值进行时空插值.并基于插值结果,在三维XYZ立方体中,以XOY平面为空间坐标,Z轴为时间坐标,绘制出各阈值对应的时空指示克里格预测立方体(图4).

图4 4种阈值下时空指示预测立方体Fig.4 The estimated cubes obtained by MTSTIK under 4thresholds

表3 四种时空立方体基本统计特征Table 3 Statistical characteristics of spatiotemporal indicators of PM2.5 estimated cubes under 4thresholds

对于各时空立方体,进行基本统计特征分析.由表3和图5可知,对于阈值Zc=35µg/m3,时空立方体数据值的平均值为 0.1309,四分位距为0.0593,数据值小于 0.2 的概率(ρ)为 81%,数据值大于0.8的概率为7%,说明其时空立方体数据值小,数据离散程度低,境内全年空气质量以超过0.8的概率达到空气质量优级别的时空占比为7%;对于阈值Zc=75µg/m3,时空立方体数据值的平均值为0.4945,四分位距为0.8782,数据值小于0.2的概率和大于0.8的概率均为34%,其余3个统计层次的概率值均为 10%左右,说明时空立方体数据值分布较为对称,离散程度高,境内全年空气质量以超过0.8的概率达到和超过轻度污染级别的时空占比均为 34%;对于阈值Zc=150µg/m3和Zc=250µg/m3,其时空立方体数据值的平均值分别为0.8111和0.9545,四分位距依次为0.2728和0,数据值大于0.8的概率均依次为74%和93%,说明其时空立方体数据值大,数据离散程度低,境内全年空气质量以超过0.8的概率超过严重污染级别的时空占比为1%.

图5 4种时空立方体日均概率百分比Fig.5 Percentages of daily average probability in estimated cubes under 4thresholds

2.4 PM2.5时空分布不确定性分析

2.4.1 山东省年均概率不确定性分析 根据各阈值下的时空预测立方体数据,获取各空间点位的日概率值,取其年平均值,绘制相应的山东省PM2.5年均概率分布(图6所示).

由图6可知,对于阈值Zc=35μg/m3,山东省境内东部沿海地域年均概率为 0.3～0.4,局部地域大于0.5,其余地域均小于 0.2,说明山东省大部分地域空气质量达到优级别的概率小于 0.2;对于阈值Zc=75μg/m3,东部沿海地域年均概率为0.6～0.9,其余地域为0.4～0.6,说明东部沿海空气质量达到轻度污染及以上级别的概率均大于0.6,其余地域大于0.4;对于阈值Zc=150μg/m3和Zc=250μg/m3,全境年均概率均大于0.8,部分地域接近于1,说明山东省境内空气质量超过严重污染级别的概率小于0.1.

此外,各尺度东部沿海空气质量达到各空气质量级别的概率均大于中部和西部,结合庄欣等[8]在对珠三角 PM2.5污染情况分析以及山东省地理位置和地理环境可知,受大气气流运输的影响,西部地区与河北、河南接壤,而陆雅静等[21]研究发现,河北省空气质量污染严重,从而使得山东省西部PM2.5污染较为严重;而东部则临海,距离河南河北等省较远,外省气流抵达东部区域时,空气中的颗粒物经过自然沉降后,气流携带颗粒物浓度降低,造成东部受到外省污染情况的影响远小于西部.

2.4.2 山东省月均概率不确定性分析 由上述分析可知,虽在不同阈值下山东省年均概率分布不确定性结构特征不尽相同,但在后续对各阈值下PM2.5月均浓度不确定性结构特征进行分析后发现,各阈值对应概率的时空分布模式大体相似,同时根据表 1空气质量分指数对健康的影响情况,在避免重复分析各阈值下时空不确定性特征的条件下,本文选取阈值Zc=75μg/m3的时空预测结果来分析山东省PM2.5月均概率和部分城市日均概率的时空分布特征(图7).

由图 7可知,从空间上,各月份达到 PM2.5≤75μg/m3的最小概率均位于西部区域(如菏泽、聊城、济宁等),尤其是在1月份,部分城市概率小于0.1;各月份概率最大位于东部沿海区域(青岛、烟台、威海等),绝大部分月份概率超过 0.6;中部为过渡区域,不同月份概率变化较大整体上呈现出从西至东,概率值逐渐增大,说明山东省空气质量从西至东,污染程度逐渐减轻,具有明显的空间分布特征,这与梅杨[10]在研究山东省 PM2.5质量浓度时空分布时得到的结论相同.

从时间上,各区域达到 PM2.5≤75μg/m3的概率值最大为7、8月份,大部分地域概率值大于0.6,部分地域概率超过 0.9;概率值最小为 1月份,大部分地域概率小于0.2,局部地域介于0.5～0.6之间.整体上,从 1～12 月,概率值先增大后减小,说明山东省空气污染程度从 1～12月,先减轻后加剧,具有明显的时间分布特征,这与康桂红等[22]对山东省2014～2015年间,不同月份PM2.5污染程度分析后得到的结果,以及成亚利对上海市 2014年PM2.5污染情况分析得到的结果相一致[23].

2.4.3 山东省部分城市单一阈值日均概率不确定性分析 根据山东省行政区划和城市经济排名,依次在山东省北部、西部、中部、南部、中东部以及东部沿海各选择一个城市作为PM2.5日均质量浓度不确定性分析观测点(北部选择东营市,西部选择济宁市,中部选择济南市,南部选择临沂市,中东部选择潍坊市,东部沿海选择烟台市),获取各城市城区范围在阈值 Zc=75µg/m3的日均概率值,并绘制相应的百分比(图8).

图6 多阈值时空指示克里格年均概率分布Fig.6 Spatial distribution of annual average probability by MTSTIK under various thresholds

图7 PM2.5≤75µg/m3下时空指示克里格月均概率分布特征Fig.7 Spatial distribution of monthly average probabilities for PM2.5≤75µg/m3

图8 山东省部分城市日均概率百分比Fig.8 Percentages of daily average probabilities for hot-spot cities in Shandong province

由图 8 可知,对于 PM2.5≤75μg/m3,6 个城市全年日均概率大于 0.75的占比中,最高为烟台市,超过70%,最低为济南市,仅为36%,说明烟台、东营、济宁、临沂、潍坊和济南6个城市全年的空气污染程度以0.75的概率达到轻度污染的天数最高为烟台市,260d左右,最低为济南市,仅有130d;全年日均概率小于0.25的占比中,除烟台市(13%)以外,其余五个城市比例均高于 30%(东营市最高,为 41%),说明全年空气污染程度以大于0.75的概率超过轻度污染的天数依次为48,150,139,130,105d.各城市空气污染程度依次为：烟台＜东营＜潍坊＜济南＜济宁≈临沂,大致符合从西至东,空气污染逐渐减轻的特点.此外,根据山东省《2014山东省环境状况公报》[24]显示,6个城市年均 PM2.5浓度大小顺序依次为：烟台＜潍坊＜东营＜济宁＜临沂＜济南[24],说明济南市PM2.5日均浓度相较于其他5个城市,其波动最为明显.

3 结论

3.1 山东省境内 PM2.5的空间自相关范围超过100km,时间自相关范围为3d左右.

3.2 对于阈值 Zc=35μg/m3,其时空预测立方体数据值偏小、离散程度低,各空间点位以超过0.8的概率达到空气质量优级别的时空占比为 7%;对于阈值 Zc=75μg/m3,其时空立方体数据值呈对称分布,离散程度高,其以超过 0.8的概率达到空气质量轻度污染级别的时空占比为 34%;对于阈值 Zc=150μg/m3和 Zc=250μg/m3,其时空立方体数据值较大、离散程度低,其以超过0.8的概率超过严重污染级别的时空占比为1%.

3.3 在年均概率分布上,山东省大部分地域空气质量达到优级别的概率介于 0～0.1;达到轻度污染级别的概率为东部大于 0.6,其余地区大于0.4;超过严重污染级别的概率均小于0.1.

3.4 从空间上,各月份达到 PM2.5≤75μg/m3的最小概率均位于西部区域,概率最大位于东部沿海区域,中部为过渡区域,.整体上呈现出从西至东,概率值逐渐增大;从时间上,月均概率值最大为7、8月份,概率值最小为1月份.整体上,从1～12月,概率值先增大后减小.

3.5 在济宁、济南、临沂、东营、潍坊和烟台6个城市中,全年空气污染程度以大于0.75的概率值达到和优于轻度污染的天数中,烟台市最高(超过260d),济南市最低(130d左右);以小于0.75的概率值超过轻度污染的天数中,最高为东营(150d左右),最低为烟台(48d).

[1]王庚辰,王普才.中国 PM2.5污染现状及其对人体健康的危害[J]. 科技导报, 2014,32(26)：72-78.

[2]戴海夏,宋伟民.大气 PM2.5的健康影响 [J]. 环境卫生学杂志,2001,28(5)：299-303.

[3]张文丽,崔九思,徐东群,等.大气中细颗粒物的污染特征及其生物效应 [J]. 中国环境卫生, 2003,(1)：90-102.

[4]覃知还.关于PM2.5的十个问答 [J]. 证券导刊, 2011,(46)：93-96.

[5]Weichenthal S, Villeneuve P J, Burnett R T, et al. Long-term exposure to fine particulate matter： association with nonaccidental and cardiovascular mortality in the agricultural health study cohort[J]. Environmental health perspectives, 2014,122(6)：609-15.

[6]李梓铭,孙兆彬,邵 勰,等.北京城区PM2.5不同时间尺度周期性研究 [J]. 中国环境科学, 2017,37(2)：407-415.

[7]Gehrig R, Buchmann B. Characterising seasonal variations and spatial distribution of ambient PM10and PM2.5concentrations based on long-term Swiss monitoring data [J]. Atmospheric Environment, 2003,37(19)：2571-2580.

[8]庄 欣,黄晓锋,陈多宏,等.基于日变化特征的珠江三角洲大气污染空间分布研究 [J]. 中国环境科学, 2017,37：2001-2006.

[9]周 甜,闫才青,李小滢,等.华北平原城乡夏季PM2.5组成特征及来源研究 [J]. 中国环境科学, 2017,37(9)：3227-3236.

[10]梅 杨,党丽娜,杨 勇,等.基于时空克里格的PM2.5时空预测及分析 [J]. 环境科学与技术, 2016,39(7)：157-163.

[11]Christakos G, Yang Y, Wu J P, et al. Improved space-time mapping of PM2.5distribution using a domain [J]. Ecological Indicators, 2018,85：1273-1279.

[12]Begum B A, Biswas S K, Hopke P K. Temporal variations and spatial distribution of ambient PM2.5and PM10concentrations in Dhaka, Bangladesh [J]. Scie of the Total Envi [J], 2006,358(1-3)：36—45.

[13]Yu H L, Chiting C, Lin S D, et al. Spatiotemporal analysis and mapping of oral cancer risk in Changhua County (Taiwan)： an application of generalized Bayesian maximum entropy method.[J]. Annals of Epidemiology, 2010,20(2)：99.

[14]Hu J, Wang Y, Ying Q, et al. Spatial and temporal variability of PM2.5, and PM10, over the North China Plain and the Yangtze River Delta, China [J]. Atmospheric Environment, 2014,95(1)：598-609.

[15]Kloog I, Chudnovsky A A, Just A C, et al. A New Hybrid Spatio-Temporal Model For Estimating Daily Multi-Year PM2.5Concentrations Across Northeastern USA Using High Resolution Aerosol Optical Depth Data [J]. Atmospheric Environment, 2014,95(1)：581-590.

[16]中国环境科学研究院.环境空气质量标准 [M]. 北京：中国环境科学出版社, 2012.

[17]环境保护部环境监测司,中国.“十二五”环境监测工作手册 [M].中国环境科学出版社, 2012.

[18]Yang Y, Christakos G. Uncertainty assessment of heavy metal soil contamination mapping using spatiotemporal sequential indicator simulation with multi-temporal sampling points [J].Environmental Monitoring & Assessment, 2015,187(9)：571.

[19]杨 勇,梅 杨,张楚天,等.基于时空克里格的土壤重金属时空建模与预测 [J]. 农业工程学报, 2014,30(21)：249-255.

[20]Cressie N, Huang H C. Classes of Nonseparable, Spatio-Temporal Stationary Covariance Functions [J]. Publications of the American Statistical Association, 1999,94(448)：1330-1339.

[21]陆雅静,王 辉,倪爽英,等.河北省 2013～2014年环境空气质量现状对比分析 [J]. 安徽农业科学, 2015,(14)：271-274.

[22]康桂红,孙兴池,韩永清,等.山东省大气污染时空分布特征分析[J]. 山东气象, 2016,36(1)：13-17.

[23]成亚利,王 波,上海市 PM2.5的时空分布特征及污染评估 [J].计算机与应用化学, 2014,31(10)：1189-1192.

[24]山东省环保保护厅.2014年山东省环境状况公报 [EB/OL].http：//www.zhb.gov.cn/gkml/hbb/qt/201506/t20150604_302942.html.