2013-2020年京津冀地区PM2.5浓度时空变化模拟及趋势分析

杨 晓 辉，肖 登 攀,*，柏 会 子，唐 建 昭，王 卫，郭 风 华，刘 剑 锋

(1.河北省科学院地理科学研究所/河北省地理信息开发应用技术创新中心，河北 石家庄 050011；2.河北师范大学地理科学学院，河北 石家庄 050024；3.河北省环境演变与生态建设实验室，河北 石家庄 050024)

0 引言

PM2.5(空气动力学粒径<2.5 μm的细颗粒物)可以携带有毒有害物质[1]，长期暴露在PM2.5下会对人类健康产生不良影响[2,3]。中国快速的城市化和工业化导致PM2.5成为空气污染的主导因素[4]，全面了解PM2.5的空间分布具有重要意义[5]。然而，由于地面监测站点数量稀疏，可获得的监测数据有限，严重阻碍PM2.5的相关研究[6,7]。卫星遥感可为地面PM2.5污染监测提供有价值的信息，特别是对于无PM2.5监测站点的地区[8,9],卫星反演的气溶胶光学厚度(AOD)与PM2.5浓度高度相关[8,10,11]，高空间分辨率的AOD产品能以更精细的尺度对PM2.5浓度进行更精确的评估[12-14]。2018年美国国家航空航天局(NASA)基于多角度大气校正(MAIAC)算法[15]生成1 km高空间分辨率AOD产品(MCD19A2)，在预测PM2.5浓度方面显现出优越性能[16-18]。近年来，许多学者通过比例因子法、半经验公式和统计回归模型等探索PM2.5浓度和卫星AOD数据之间的定量关系，其中，统计回归模型以快速、简单和准确等特点得到广泛应用。从早期研究的简单线性回归[10]到先进的高级统计模型，如线性混合效应(LME)模型[19]、广义加和(GAM)模型[20]、地理加权回归(GWR)模型[21]、时空线性混合效应(STLME)模型[22,23]和地理时间加权回归(GTWR)模型[24]等，之后通过组合两个模型(如LME+GWR、LME+GAM和TEFR+GWR等)[25-27]或三阶段统计模型(如IPW+GAMM+KED)[28]减少单一模型预测带来的偏差[29]。此外，部分学者虽然使用随机森林(RF)模型[30]、自适应深度神经网络(SADNN)模型[6]和支持向量机(SVM)[31]等机器学习模型预测PM2.5浓度，但不能给出特定的模型参数以解释PM2.5-AOD的时空关系。

研究表明，在模型中添加交互项(AOD的二次项)可以更好地描述非线性效应[32]，在PM2.5污染较严重区域能纠正PM2.5浓度预测偏差[32，33]。因此，本研究选择LME+GWR两阶段统计回归模型以精确捕捉PM2.5-AOD的时空变化，以MAIAC AOD为主要变量，气象和土地利用因子为辅助变量，预测每日地面PM2.5浓度，同时模型中加入AOD的二次项以及行星边界层高度(PBLH)和AOD的交互项，以解释PM2.5-AOD的非线性关系，并使用十折交叉验证评估模型预测的准确性，最后利用模型建立长时间序列高空间分辨率的PM2.5浓度数据集，分析PM2.5地基监测网络建立以来京津冀地区2013―2020年PM2.5浓度时空变化趋势，以期为京津冀地区大气污染防治工作提供参考和建议。

1 研究区与数据

1.1 研究区概况

京津冀地区是我国北方重要的行政、经济和文化中心，包括北京市、天津市和河北省的11个地级市(图1)，该地区人口稠密，以煤炭为主要能源的第二产业排放各种大气污染物，尤其在内陆平原地区，大气污染物难以扩散，空气污染严重。2013―2020年《中国环境公报》(http://www.cnemc.cn/jcbg/zghjzkgb/)指出：中国空气质量较差的10个城市中，京津冀地区历年分别占7、8、7、6、6、5、4和1个。虽然近年来该地区空气质量有所改善，但仍需重视PM2.5浓度的时空分布特征和变化趋势研究。

图1 研究区PM2.5监测站点分布Fig.1 Study area with PM2.5 monitoring stations

1.2 数据来源和处理

基于AOD的PM2.5浓度预测模型加入气象和土地利用等因子，可显著提高模型预测精度[22,34]。本研究采用LME+GWR两阶段统计回归模型，除主要预测因子AOD外，其他辅助预测因子包括行星边界层高度(PBLH)、近地面温度(TEMP)、近地面风速(WS)、近地面相对湿度(RH)、近地面气压(PRS)、降水量(PRCP)、森林覆盖率(FC)和建筑用地覆盖率(UC)，通过预测变量选择和多重共线性分析筛选得以利用。数据集涵盖时间为2013年1月1日至2020年12月31日，详细信息见表1。

表1 数据来源与时空分辨率Table 1 Data source and spatial-temporal resolution

(1)PM2.5监测数据。京津冀80个监测站点PM2.5浓度的小时尺度数据均来自全国城市空气质量实时发布平台。为确保PM2.5数据有效性，在日均PM2.5浓度拟合过程中，剔除不在国家环境空气质量标准(NAAQS)[35]监测范围内的数值(即PM2.5浓度<2 μg/m3或PM2.5浓度>500 μg/m3)。

(2)MAIAC AOD数据。MAIAC算法在反演过程中利用时间序列分析和空间分析增强了对比局部尺度气溶胶变化的能力，提高了云雪检测和气溶胶反演的质量[36]。从NASA获得MODIS C6 MAIAC Terra/Aqua AOD(0.55 μm)产品，利用气溶胶地基观测网(http://aeronet.gsfc.nasa.gov/)中Beijing、Beijing-CAMS和Xianghe的 AERONET AOD 数据(版本 3，级别 2)的验证结果(图2)表明，决定系数R2为0.81～0.91，均方根预测误差RMSPE为0.10～0.25 μg/m3，73.15%～83.74%的样本落在期望误差区间(±(0.05+20%×AERONET AOD))内，满足验证精度要求。

(3)气象数据。2013-2020年PBLH数据来源于GEOS5-FP的网格化数据产品；2013-2018年其余气象数据(TEMP、WS、RH、PRS和PRCP)来源于国家青藏高原数据中心[37]，由于数据集只涵盖1979-2018年，从中国气象数据网下载分辨率相近的2019-2020年的气象数据。

(4)土地利用数据。土地利用数据来源于地理国情监测云平台，本研究选取2015年土地利用数据代表2013-2020年的土地利用状况，提取计算整个研究区建筑用地覆盖率和森林覆盖率加入模型。

为使各类数据的时空分辨率与每日PM2.5浓度相匹配，日均PBLH数据由MODIS卫星过境期间两次观测均值获得；所有日均气象变量使用双线性内插法重新采样到与MAIAC AOD相同空间范围(1 km×1 km)；森林和建筑用地覆盖率在以每个PM2.5监测站点为中心的1 km×1 km缓冲区内进行平均。将上述MAIAC AOD、气象和土地利用变量按经纬度提取到80个监测站点，并根据站点和日期将提取结果与PM2.5浓度数据进行匹配，最后，删除MAIAC AOD和PM2.5的非随机缺失数据，2013-2020年分别获得13 272(300 d)、15 413(325 d)、15 657(316 d)、21 882(330 d)、16 961(323 d)、10 942(300 d)、10 987(281 d)和11 698(299 d)条数据记录。

注：红色虚线是回归线,黑线是1∶1 线,灰线代表期望误差区间。

2 研究方法

2.1 共线性诊断

为提高预测模型的稳定性，必须考虑自变量的共线性。本研究选择方差膨胀因子(VIF)(式(1))和容差(TV)(式(2))诊断所选变量之间的共线性[38]。VIF越大，说明自变量之间存在共线性的可能性越大，当自变量的TV>0.1且VIF<10时，表明自变量之间无共线性问题存在。由表2可知，参与模型的所有自变量每年的VIF和TV均满足VIF< 10和TV> 0.1的条件，表明自变量之间不存在共线性问题，均可考虑用于模型拟合。

表2 变量共线性分析中方差膨胀因子(VIF)和容差(TV)范围Table 2 Range of variance inflation factor (VIF) and tolerance value (TV) in the analysis of variable collinearity

(1)

(2)

式中：Ri为自变量Xi对其余自变量作回归分析的负相关系数。

2.2 标准差椭圆(SDE)

标准差椭圆(SDE)是基于一组离散点的平均中心概括地理要素(本文为PM2.5浓度)空间特征(集中趋势、离散度和方向趋势)[39]的方法，输出参数包括标准差椭圆、平均中心、长短轴和方位角[40]。其中，标准差椭圆表示要素空间分布的主体区域，其变化可以描述整个空间动态过程的特征；平均中心的移动可揭示要素的整体演化轨迹；椭圆长轴和短轴的变化表示特定空间方向的扩展或收缩；方位角的变化表示整体要素在特定空间方向上的变化[41]。

2.3 模型构建与验证

本研究采用由线性混合效应(LME)模型和地理加权回归(GWR)模型组成的两阶段统计回归模型对PM2.5-AOD关系的时空变化进行模拟，第一阶段LME模型用以校正PM2.5-AOD的时间变化关系，模型中添加了AOD的二次项以及AOD和PBLH之间的相互作用，以解释 AOD和PM2.5之间的非线性关系。计算公式为：

(3)

第二阶段GWR模型用以校正PM2.5-AOD的空间异质性，通过对LME模型的残差进行建模，每天拟合一次以说明时间可变性，具体计算公式为：

PM2.5_resist=β0(us,vs)+β1(us,vs)AODst+εst

(4)

式中：PM2.5_resist为LME模型中地面监测站点s第t天的残差值；AODst为监测站点s第t天的MAIAC AOD值；β0(us,vs)和β1(us,vs)分别为空间坐标为(us,vs)监测站点s处的回归截距和回归斜率。

本文采用十折交叉验证[29]检测模型的过拟合程度，首先，整个模型拟合数据集被随机分成10个子集，每个子集包含大约10%的数据集，在每次交叉验证中选择一个子集作为测试样本，并使用剩余的9个子集拟合模型以对测试样本进行预测，此过程重复10次以确保所有子集都被预测。此外，对实测和预测的PM2.5浓度进行线性回归拟合，利用拟合后R2、斜率、RMSPE和相对预测误差(RPE)评估模型性能。计算公式为：

(5)

(6)

3 结果与分析

3.1 统计分析

京津冀地区地面监测站点观测PM2.5年均浓度范围为40.97～91.27 μg/m3，2013年最高，2020年最低，表明过去8年间PM2.5污染程度呈下降趋势，但仍超过国家环境空气质量二级标准(GB3095-2012)的限值(35 μg/m3)。同期年平均AOD范围为0.37～0.69，森林和建筑用地年均覆盖率分别为3%～5%和55%～78%，2013―2020年各气象变量的范围见表3。

表3 建模变量统计指标Table 3 Statistical indicators of modeling variables

3.2 模型拟合与验证

2013-2020年LME+GWR模型拟合和交叉验证结果对比(图3)显示，在捕捉每日PM2.5浓度方面，模型表现性能优异。模型拟合时，数据分布向回归线集中，模型拟合R2范围为0.89～0.97，表明两阶段模型可以解释地面89%～97%的PM2.5浓度变化；斜率为0.89～1.04，表明模型的预测偏差较小，能较好地克服高值低估问题；RMSPE和RPE分别为6.85～24.60 μg/m3和16.67%～26.95%。与模型拟合相比，交叉验证后模型R2为0.85～0.95，RMSPE为7.87～29.90 μg/m3，RPE为19.19%～32.72%，可以看出，交叉验证后模型R2减小，而RMSPE和RPE增加，表明模型存在轻微的过度拟合。

图3 2013-2020年模型拟合和交叉验证结果对比Fig.3 Comparison of model fitting and cross-validation results from 2013 to 2020

不同年份中，2020年模型性能最好，交叉验证后模型R2最高，为0.95，RMSPE和RPE值最小，分别为7.87 μg/m3和19.19%，这主要是由于随着政府对“大气污染防治行动计划”的开展，PM2.5浓度出现下降。2020年大约93.72%的数据样本PM2.5浓度小于100 μg/m3。相比之下，2013年模型性能最差，交叉验证后模型R2最低，预测不确定性最大，主要原因是超过32.38%的数据样本PM2.5浓度大于100 μg/m3，相对离散的数据样本可能增加了模型拟合的难度。综上，LME+GWR模型稳健，使用1 km MAIAC AOD产品结合LME+GWR模型在京津冀地区可以很好地模拟每日地面PM2.5浓度。

3.3 京津冀地区PM2.5浓度时空分布

基于标准差椭圆分析2013-2020年京津冀地区PM2.5浓度空间格局的逐年变化(图4)，可以看出，PM2.5主要集中在京津冀平原区，说明在整个研究期内，京津冀平原区是全域PM2.5总量的主要贡献者。PM2.5浓度重心逐渐向东北方向移动，仍处于沧州市内。椭圆长轴从2013年的242.17 km减至2017年的239.03 km，表明该时段京津冀地区PM2.5年均浓度在主要方向上呈现空间集聚趋势，但集聚程度不大，之后又轻微增至2020年的242.18 km，表明该时段PM2.5浓度在主要方向上出现分散趋势；短轴从2013年的101.40 km增至2020年的105.11 km，表明京津冀地区PM2.5在地理空间上呈扩张分散化趋势；从方位角变化范围看，PM2.5浓度总体呈“西南—东北”方向空间分布格局，该地理空间走向与京津冀地形密切相关；转角由2013年的40.53°波动缩小至2020年的39.56°，表明PM2.5空间分布格局存在沿西北方向微弱偏移趋势。

图4 2013-2020年PM2.5浓度重心和标准差椭圆的空间变化Fig.4 Spatial changes of PM2.5 concentrations gravity center and standard deviational ellipse from 2013 to 2020

基于两阶段统计回归模型预测的2013-2020年京津冀地区PM2.5年均浓度分布状况(图5)可知，PM2.5浓度空间分布趋势大致相同，南北差异显著。PM2.5浓度低值区位于西部和北部山区，高值区位于平原，与标准差椭圆分析结果一致。总体而言，PM2.5浓度呈现出“北部山区低，南部平原高(平原西南高浓度、东北中等浓度)”的空间分布格局。2013-2020年京津冀地区PM2.5年均浓度分别为69.67、65.31、49.26、51.17、44.96、43.11、34.54和32.02 μg/m3，可以看出，2013-2014年PM2.5年均浓度波动较小，2015年出现大幅下降，2017年PM2.5年均浓度下降35.47%，达到“大气污染防治行动计划”中“到2017年下降25%”的要求，主要原因在于2017年是《大气十条》实施效果的验收年，政府防治力度加大[42]。此外，2019年PM2.5年均浓度下降更明显，2020年下降54.04%，相对2017年下降28.78%，PM2.5污染程度总体下降显著。从区域看，高浓度区明显缩小，污染城市主要集中在邯郸、邢台、石家庄、保定和衡水。2013年高浓度区范围最大，2017-2020年低浓度区域逐渐扩大，2019和2020年京津冀全域PM2.5浓度降至55 μg/m3以下。

图5 2013-2020年PM2.5年均浓度空间分布Fig.5 Spatial distribution of annual mean PM2.5 concentrations from 2013 to 2020

京津冀地区PM2.5污染也表现出较强的季节性变化。整体看，PM2.5浓度呈现“冬季高，夏季低，春秋过渡”的季节变化特点(图6)。2013年和2014年全年PM2.5浓度较高，冀中南地区春季PM2.5浓度大于85 μg/m3，冬季大于115 μg/m3；2015-2020年京津冀地区PM2.5浓度逐渐降低，夏季均小于55 μg/m3，2019-2020年甚至低于30 μg/m3；2013-2018年秋季、冬季和春季污染范围逐渐缩小，污染严重地区主要集中在邯郸、邢台、保定和石家庄等市；2019-2020年春夏秋三季均小于55 μg/m3。2013-2020年冬季PM2.5浓度均值依次为117.46 μg/m3、84.24 μg/m3、75.30 μg/m3、72.72 μg/m3、55.97 μg/m3、44.63 μg/m3、51.48 μg/m3和51.42 μg/m3，逐年下降并在2017年后趋于稳定，2018-2020年冬季重污染区范围显著减小，轻污染区范围扩大，导致相对2017年PM2.5浓度变化不大；2017年和2020年相比2013年冬季PM2.5浓度分别下降了61 μg/m3(52%)和66 μg/m3(56%)。

图6 2013-2020年PM2.5季均浓度空间分布Fig.6 Spatial distribution of seasonal mean PM2.5 concentrations from 2013 to 2020

为了解京津冀地区各市域之间的差异，利用预测结果计算各城市PM2.5年均浓度(图7)。2013-2020年承德和张家口PM2.5浓度较低，8年内均低于35 μg/m3；邯郸、邢台和衡水PM2.5浓度居前三位，2013年PM2.5浓度均超过100 μg/m3，到2020年PM2.5浓度分别降为47.17 μg/m3、47.28 μg/m3和50.05 μg/m3，与沧州和廊坊接近。此外，2013-2017年北京PM2.5污染在13个城市中排名第10位，2018年之后降至第11位，属中低污染水平，可能缘于重污染行业的搬迁。

图7 2013-2020年各城市PM2.5年均浓度Fig.7 Annual mean PM2.5 concentrations in each city from 2013 to 2020

4 结论与讨论

本文基于京津冀地区2013-2020年的日均PM2.5浓度数据，结合MAIAC AOD、气象和土地利用等因子，利用两阶段统计回归模型模拟京津冀地区每日PM2.5浓度的时空分布，其中第一阶段采用LME模型反映PM2.5-AOD关系随时间变化规律，第二阶段加入GWR模型反映PM2.5-AOD关系的空间变异，主要研究结论如下：1)十折交叉验证结果表明LME+GWR模型稳健，使用1 km MAIAC AOD产品结合LME+GWR模型对京津冀地区PM2.5-AOD之间的时空变化模拟性能优异；2)PM2.5浓度总体呈现出“西南—东北”方向空间分布格局，且存在沿西北方向微弱偏移趋势；PM2.5浓度的重心向东北轻微移动，结合椭圆的长、短半轴可知，京津冀PM2.5浓度呈现轻微扩张分散化趋势；3)由LME+GWR模型模拟得到的京津冀地区2013-2020年地面PM2.5浓度整体呈现“南部平原高、北部山区低”和“冬高夏低、春秋过渡”的时空特征，2018-2020年冬季低污染区域有扩大趋势，需格外重视。

与已有研究相比，本研究构建的模型(LME+GWR)性能良好。对于低空间分辨率(3～10 km)的 PM2.5模拟研究[26,43,44]，该模型具有更高的R2和更低的RMSPE；与使用高空间分辨率(1 km)AOD预测PM2.5浓度的模型[12,45]性能相似。尽管本研究使用了Terra和Aqua融合后的MAIAC AOD数据，但仍没有达到区域内全覆盖。AOD数据缺失会对模型预测精度产生一定影响，主要表现为“高值低估”[37]，这是后续研究的方向之一。如Chen等[46]开发了混合效应模型和逆距离加权(IDW)两步插值技术填补AOD缺失值，整体上AOD缺失率从87.91%(Auqa)和85.47%(Terra)降低至融合后的13.83%且R2达0.76。 此外，未来可以在更大的研究区域内考虑更多相关时空因素，以扩大研究区域并提高模型模拟精度。