基于卫星遥感AOD的华北地区2003—2014年PM2.5浓度 时空分布特征

陈优芳，周一敏，赵昕奕

(1. 北京大学城市与环境学院，北京100871；2. 地表过程分析与模拟教育部重点实验室，北京100871)

1 引 言

气溶胶细颗粒物(空气动力学直径小于2.5 μm，PM2.5) 在霾的形成和发展中扮演主要角色[1]，由此造成的空气污染导致了许多严重的环境问题，如城市能见度下降、影响交通安全、粮食产量下降等[2-4]。更严重的是，由于体积很小，它们可渗入人体支气管和肺部，进而引起严重的心血管疾病、呼吸道疾病甚至肺癌[5]。最新的研究已经证实PM2.5会影响人体发病率和死亡率水平，对大众健康有重要影响[6]，并在所有健康风险影响因素中被列为第九大影响因素[7]。由于工业的快速发展、发达的交通网络和密度较大的人口，华北地区是世界上大气气溶胶浓度最高的地区之一[8-9]。自1980年代以来，PM2.5污染天数显著增加[10-12]，更重要的是，长期的重污染事件越来越频繁。根据中华人民共和国环境保护部(MEP)发布的2013年度74个主要城市的年度空气质量数据，京津冀地区所有城市的PM2.5年平均浓度均超标，其中7 个城市排名前10 位。一些城市的空气污染天数占全年总天数的40%以上。

重污染事件给人类健康、气候、食品生产、交通安全等带来了一系列不利影响。因此，研究PM2.5的时空变化特征对进一步预测PM2.5浓度和防治PM2.5污染具有重要意义。基于地面站点的观测可提供准确、及时的PM2.5浓度。但是，因为运营成本和技术要求较高，地面站点的空间和时间覆盖范围有限，无法捕捉到高时空分辨率的PM2.5浓度。而卫星遥感可提供长时间序列、大范围的观测数据而被广泛应用于各领域研究中[13-16]。在PM2.5浓度研究中，卫星遥感可提供各种尺度的气溶胶光学厚度(AOD) 产品，从而反演出未覆盖地面PM2.5观测区域的PM2.5浓度[17-21]。AOD 是气溶胶最基本的光学特性之一，定义为介质的消光系数在垂直方向上的积分，可用来描述大气气溶胶对光的衰减作用，推算大气气溶胶的含量、评估大气污染的程度与确定气溶胶的气候效应[22-23]。卫星观测AOD 始于1970s，目前已发展出包括MODIS、MISR、GEOS、POLDER、OMI、SeaWiFS、CALIOP等在内的多数据来源、多种分辨率的卫星AOD 数据，许多研究使用这些数据对气溶胶浓度进行反演[24]。在众多卫星AOD 产品中，中分辨率成像光谱仪(MODIS)AOD 产品由于其高分辨率和准确度而被广泛应用于地面PM2.5浓度的反演[25]。反演效果被许多研究认可，通过与地面观测系统的比较，一般认为MODIS AOD 的系统误差为±0.05，反演算法带来的误差约为15%，达到设计目标[9,26-27]。焦利民等[28]发现京津冀地区AOD 与PM2.5浓度具有很高的空间相关性，二者年均值相关系数约为0.6，但相关性存在时空差异，秦玮等[29]对南京地区、范萌等[30]对珠三角地区的研究也有类似结论。可见，利用MODIS_AOD 对某区域PM2.5演变情况进行研究是可行的。

反演方法主要包括简单线性模型，土地利用回归模型，地理加权回归(GWR)和具有附加变量的复杂模型等[31-33]。这些方法中PM2.5与AOD 之间的关系是固定的，但是由于PM2.5浓度受排放条件、气象条件等影响，而这些影响因素具有时空变化，因此PM2.5-AOD 关系随着时空变化而变化[34]。因此，我们在反演时使用了混合效应模型(Mixed effect Model)。混合效应模型是包含固定效应和随机效应的统计模型。这些模型在对相同统计单位进行重复测量的情况下特别有用。在1950年代，亨德森提供了固定效应的最佳线性无偏估计(BLUE) 和随机效应的最佳线性无偏预测(BLUP)[35]。从那时起，混合建模已成为统计研究的一种主要手段。Lee 等[36]通过开发线性混合效应(LME)模型将这种方法引入到PM2.5反演研究中。以前的研究表明，LME 模型的表现优于其他统计模型，如地理加权模型[36]和线性回归模型[37]。

目前，现有的研究主要是结合特定的污染过程中气象要素和环流对PM2.5浓度的影响进行研究，即对某次短时间内的霾事件过程进行研究[38-39]，但在长时间尺度的研究上，对于华北地区PM2.5浓度的时空分布特征研究较少，难以为PM2.5污染的预测和防治提供更有效的信息。因此，本研究利用了2013—2015年MODIS 的AOD 3 km 分辨率产品(MYDO4_3k)和环保部公布的2014年4月—2015年3月各个城市的PM2.5浓度数据对华北地区长时间序列的PM2.5浓度的时空分布特征进行了分析。

2 数据来源与研究方法

2.1 研究区概况

本研究的研究区域为华北地区及其周边小部分地区，包括北京、天津、河北、山东、陕西、河南、内蒙古、辽宁等八个省区市，位于113～123 °E，34～43 °N 之间(图1)。整个研究区域的地形呈西部高、东南低的分布，是大气气溶胶分布的重要成因之一。由于丰富的矿产资源，华北平原是我国的重要工业基地，其产业结构也是形成空气污染的重要原因之一。

图1 研究区及测站城市示意图

2.2 研究数据

2.2.1 卫星遥感数据

卫星遥感数据为NASA 发布的2003年1月—2015年3月的MODIS 传感器MYD04(http://ladsweb.nascom.nasa.gov/data/search.html) 的 3 km 分辨率AOD 逐日数据产品。该卫星每日经过华北地区的时间为北京时间13 时，即该卫星每日对华北地区的观测时间为每日13 时。MYD04_3k数据包括了气溶胶的光学性质(例如光学厚度和粒径分布)、质量浓度等特征，具有高时空分辨率的特点。在MODIS 系列5 数据之前，气溶胶数据一般使用MYD04 产品，但是其空间分辨率只有10 km，导致目前空气污染的相关研究精度不够。因此，在系列6 产品中，NASA 使用暗目标气溶胶算法创建了3 km 气溶胶产品，MYD04_3k 与MYD04 两个产品的算法一致，但是MYD04_3k空间分辨率更高。本研究将使用与地面观测数据相匹配的2014年4—12月的数据进行模型拟合，并用2015年1—3月的卫星数据对模型效果进行检验。

2.2.2 地面观测数据

地面观测数据为环保部公布的各地(城市)实时空气质量PM2.5浓度数据(http://air.ePMap.org/)。本研究使用了上述八个省区市中地处113～123°E，34～43 °N 范围内的41 个城市总计215 个站点、2014年4月—2015年3月的逐小时PM2.5数据。为了匹配卫星观测时间，取11：00—15：00 的小时数据的平均数据作为日值数据。根据《中国国家环境空气质量标准(GB3095-2012)》，本研究中的PM2.5污染指的是PM2.5质量浓度大于75 μg/m3时的浓度，而污染日数比则是指在研究时段内的PM2.5质量浓度值大于75 μg/m3的天数总和与有PM2.5质量浓度数据的天数总和的比值。

2.3 研究方法

2.3.1 聚类分析

MYD04_3k 的AOD 数据为日值数据，观测时间为北京时间每天13 时左右。由于卫星观测的原因，MYD04_3k 的AOD 数据存在大量缺测值并且缺失情况每日不同，所以需要对缺失的观测结果进行补充。考虑到气溶胶在大气中的空间一致性和连续性，本研究通过空间聚类分析将研究区域分为几个不同区域，并且使用每一类中的有观测数据的点的AOD 大小的平均值作为该区域无观测点的替代；若是一个区域内所有点在当天均无观测值，则视为该区域在当天没有观测结果。

2.3.2 Mixed Effect Model 一元线性回归算法反演PM2.5浓度

对华北地区进行聚类分区后，本研究使用Mixed Effect Model 一元线性回归算法，通过MYD04_3k 的AOD 数据反演了2003—2014年的华北地区各个小区域PM2.5浓度值。

首先，利用地面观测的PM2.5浓度数据校准MYD04_3k 的AOD 数据。AOD 数据取该城市中当天所有点的AOD 值的平均值作为该城市当天的AOD 值。对每个城市逐日的PM2.5浓度数据和AOD 数据进行一元线性回归。本研究既不是对每个城市分别进行线性回归，也不是对所有城市的所有天数统一进行线性回归，而是参考Lee 等[36]和Xie 等[40]的方法利用Mixed Effect Model 进行线性回归。Mixed effect Model 主要是在一元线性回归中的固定参数的基础上，对不同站点和不同时间的数据添加随机参数，从而更精确地模拟PM2.5的浓度。本研究在反演过程中先将反演日期内研究区中地面观测数据和MODIS AOD 数据都有观测数据的天数挑选出来，而后在SPSS 中通过分析-混合模型-线性的工具对数据进行混合效应模型拟合，最终获得各个城市的参数。在Mixed Effect Model 中：

其中β 和α 是修正后的斜率和截距，εi为不同城市的修正值，对于不同城市来说，ε 不同，且∑εi=0。

在获得了各个城市的一元线性回归参数后，对ε 在整个研究区域进行数值插值，并利用插值的结果反演出整个研究区域的2003—2014年的PM2.5浓度数值大小，结合聚类分析的结果最终得到了华北地区2003—2014年的PM2.5浓度。

2.3.3 EOF 分解

本文使用经验正交函数(EOF)对反演的PM2.5浓度时间和空间特征分析。EOF，也称为特征向量分析或主成分分析(PCA)，是通过矩阵数据中的主要成分对数据结构进行分析的方法。 Lorenz[41]在1950年代首次将其引入气象学和气候学研究，现在广泛应用于地球科学和水电学等领域。EOF 分析能够将时变信号或数据集分解为空间函数部分和时间函数部分。空间函数部分不随时间变化并且可总结数据的地理分布特征，而仅依赖于时间变化的时间函数部分由变量的线性组合组成，因此可展现数据的时间变化。前几个分解成分占原始信号总方差的很大一部分，可看作原始信号的主要信息集中在这些成分上。该方法快速收敛，并且容易集中大量数据信息。因此，我们通过EOF分解可了解PM2.5的时间序列和空间模式。

EOF 分解过程主要如下。

(1) 将PM2.5反演数据预处理成距平形式，得到数据矩阵Xm×n，即m 个变量、n 次观测。

(2) 计算X 矩阵与其转置矩阵XT的交叉积，得到协方差矩阵，

(3) 计算方阵C 的特征根(λ1，征根λm)和特征向量Vm×m，二者满足：

其中E 是m×m 维对角阵，即，

一般将特征根从大到小顺序排列，即λ1＞λ2＞……＞λm。因为数据X 是PM2.5浓度，所以λ 应该大于或者等于0。每个非0 的特征根对应一列特征向量值，也称为EOF。如λ1对应的特征向量值称为第一个EOF 模态，也就是V 的第一列，它反映的是PM2.5的空间分布特征。

(4) 计算主成分。将EOF 投影到原始资料矩阵上，就得到所有空间特征向量对应的时间系数(即主成分)，即，

其中PC 中每行数据就是对应每个特征向量的时间系数，也就是PM2.5的时间变化特征。

3 研究结果

3.1 对研究区进行空间聚类分析

影响AOD 反演PM2.5效果的重要因素是AOD 数据的缺失。由于卫星观测数据大量缺失，必须对缺失的数据进行补充，而补充数据的前提则是对研究区进行空间聚类。空间聚类分析结果如图2，表1 对9 个区域所在位置进行了具体分析，结合EOF 分析的结果可看到，区域7 为污染情况最严重的区域，区域5、区域8 和区域6 的污染情况仅次于区域7；区域4 为污染情况快速减少的区域；区域9 是污染较少的区域；区域1、区域2、区域3 则是污染最少甚至可说几乎没有污染的区域。9 个小区域内的各个观察点PM2.5浓度的相关系数显著高于区域外的各个点，且在该区域与其他区域的分界处相关系数快速下降，说明这几个区域内部一致性较高，区域聚类效果较好。

图2 空间聚类分析结果行政区切割

表1 9 个区域所在位置

3.2 分区反演浓度

为了评估模型反演的拟合效果，计算了反演PM2.5浓度和地面观测PM2.5浓度之间的确定系数(R2)，平均误差(MPE)和均方根误差(RMSE)。

R2如图3 所示，大多数城市通过0.01 显著性检验(P＜0.01)，威海、青岛和德州通过0.10(P＜0.10)的显著性检验，秦皇岛、菏泽、大同和葫芦岛通过0.05 显著性检验(P＜0.05)，但枣庄和郑州未通过0.10 显著性检验(P＞0.10)。 因此，我们使用了通过0.10 显著性检验的39 个城市的数据来构建混合效应模型。总体而言，北方城市的相关系数略高于南方城市，尤其是河北南部、河南北部和山东南部城市的相关系数。这些城市通常是PM2.5污染严重的地区。也就是说，对于PM2.5污染严重的地区，模型拟合结果较差，而污染较少的地区则效果较好。

图3 混合效应模型中每个城市的拟合系数(R2)点的颜色越深，R2 越大。

最后对相关系数通过0.10 的显著性检验的39个城市进行Mixed Effect 一元线性回归，使用2014年4—12月的地 面PM2.5和MODIS 卫星AOD 作为模型的训练样本，2015年1—3月的数据作为检测样本，对比发现经过Mixed Effect Model 处理后，整个区域的地面观测与卫星观测的数据的相关性显著提高(拟合优度提高50%左右)。2014年4—12月共有1 973 个数据，MPE 和RMSE 分别为21.6 和30.1。此外我们利用2015年1—3月8 个城市共141 个数据检验了模型效果，MPE 和RMSE 分别为25.6 和38.5(图4)。

图5 为各个城市经过Mixed Effect 一元线性回归算法后得到的截距的空间插值结果。可看到截距在整个研究区内呈西南高东北低的趋势，特别是河北南部和山东西部地区的截距明显高于其他地区，表明在卫星观测到相同AOD 的情况下，河北南部和山东西部地区的PM2.5浓度要高于其他地区，这可能是由于不同地区的污染物化学类型不同造成的。

图6 为对研究区反演后的PM2.5浓度进行经验正交(EOF)分析得到的第一空间模态，其对整个PM2.5浓度场的方差贡献率为73%，且第二～第六空间模态的方差贡献率相当并且都较少，故认为第一空间模态可代表整个研究区PM2.5浓度的空间分布情况：河北南部、河南北部和山东西部地区为PM2.5污染(下称污染)空间高值区域，从其向周围扩散一直到太行山脉和山东半岛，整个区域的污染空间分布较均匀，梯度较小；其东部与山东半岛交界的地区污染梯度较大，随着经度增加快速减少。河北中部到北京南部、天津北部，污染的等值线与地形海拔高度的等高线几乎重合，在此地区污染情况随海拔高度的增加迅速减少。在40 °N以北地区，污染情况总体较少，特别是张家口及其周边地区几乎没有污染。

从整个EOF 分析的结果来看，地形对污染的影响非常明显，整个山脉将污染分为两个部分，山脉的东南平原地区为污染较严重的地区，而山脉西北则污染较少。其原因主要有两方面，一是平原地区为主要的工业产地，特别是污染较大的重工业；二是由于地形的阻挡，平原地区产生的污染无法扩散到山脉的西北地区，故而山脉西北地区的污染较少。

图4 通过0.10 显著性检验的39 个城市混合效应模型模拟效果 a. 2014年4—12月的反演PM2.5 浓度和实测PM2.5 浓度拟合情况；b. 2015年1—3月反演PM2.5 浓度和实测PM2.5 浓度拟合情况，灰线表示拟合的线性回归线。

图5 研究区各个城市在Mixed Effect 算法中截距空间插值结果

对数据进行反演后，利用空间聚类的结果对9个研究分区的卫星观测数据进行补充，使用每一类中有观测数据的格点的反演PM2.5浓度的平均值作为该区域无观测数据格点的替代浓度；若是一个区域内所有点在当天均无观测值，则视为该区域在当天没有观测结果。在经过数据补充之后，在研究的总共4 808日中，研究区内9 个区域数据齐全的日数共有2 017 天，8 个区域有数据的有740 天。3.3 研究区PM2.5浓度在2003—2014年的时空分布

图6 研究区EOF 第一空间模态

在补充了部分缺失数据后，通过修正过的卫星反演的PM2.5浓度得到华北地区PM2.5浓度在2003—2014年的变化情况。用每个区域中的平均值代表这个区域的污染情况得到不同区域污染日数百分比的年际变化(图7)，通过该图也可看出研究区污染情况的时空分布特征。从空间分布看，总体情况与EOF 第一主成分的空间格局类似，太行山脉是污染强弱的明显的分界线，山脉东南部的污染显著高于西北部；整个研究区北部和山东半岛地区污染情况较少；河北南部和山东西部地区是污染最严重的区域，其周围区域污染程度仅次于这一区域，南北方向上的污染程度稍高于东西方向。从区域来看，北方4 个区域(区域1～区域4)污染情况相似，几乎无污染，其中区域1 最好，区域4 最差；区域7 为污染情况最严重的地区，区域5 和8 次之，区域6 再次之。区域与区域间的相关系数也可以反映出这种空间联系性，从表2 可发现，区域1、区域2、区域3 和区域4 之间都有较高的相关系数，区域5、区域6 和区域7 之间的相关系数较高(下划线数据)，区域6、区域7 和区域8之间相关系数较高(下划线数据)；区域9 污染较少，但与其他8个区域的相关性均较低，可能与其独特的地理位置有关。

图7 整个研究区和9 个区域污染日数百分比的年际变化 蓝色系曲线为污染较少的区域，与左侧坐标相对应；红色系曲线为污染较多的区域，与右侧坐标相对应；黑色曲线为整个研究的污染日数比，与左侧坐标相对应。

从时间分布来看，2004年和2005年两年为污染情况较少的年份，从2006年开始，污染程度急剧加强，并且在随后两年保持高值，但是2009年稍稍缓和，可能与保障2008年奥运会顺利召开所采取的大气污染治理措施有关。之后污染再次加剧，不过2013年后，污染呈现下降趋势。

表2 不同区域污染情况的相关系数 均通过0.01 显著性检验。

4 讨论与展望

本研究使用了MODIS 的AOD 数据和PM2.5的地面观测数据进行统计分析。此前，大部分研究主要使用单站或者多个站点的观测结果进行研究，导致研究结果在空间和时间上的局限性。本研究使用了上述两套数据，可以提供大范围和长时间的气溶胶观测数据，大大提高了研究的可靠性。但是由于卫星遥感的AOD 有数据缺失的问题，本研究中采用了聚类分区，然后用区域平均的结果对区域中的缺失数据进行填补的方法，引入了误差。虽然两者的相关性较高，但需要探讨是否有其他更好的方法可对缺失数据进行填补。

此外，从长期或平均的状态来说，城市污染程度取决于城市的能源结构、交通和工业排放污染物的多少[42-43]，但从短期或实时状态而言，却主要与当地、当时的气象条件有关。大尺度环流形势和局地气象条件不仅可影响污染物的生成、积累和清除，同时也是影响区域输送的重要外部条件[44-45]。在华北地区，污染物浓度的变化很大程度上取决于局地的气象条件，如风速、边界层高度、相对湿度、降水等[46-47]，因而在未来需要着重对气象要素如何影响PM2.5浓度时空分布特征进行研究。

5 结 论

本研究利用了卫星遥感气溶胶光学厚度资料、地面观测的PM2.5浓度数据，采用了混合效应模型、聚类分析、EOF 分解等统计分析方法，研究了2003—2014年华北及周边地区PM2.5浓度的时空分布特征。

(1) 卫星遥感的AOD 与地面观测PM2.5浓度值有较高的相关系数，可利用MODIS 卫星遥感的AOD 数据对地面观测的PM2.5浓度进行反演。

(2) 在华北和周边地区，PM2.5浓度呈现出明显的空间分布特征，40 °N 以北的区域和山东半岛地区污染情况较少，太行山脉是污染强弱的明显的分界线，山脉东南部的污染显著高于西北部，且在地势变化的地方出现明显的突变；河北南部、河南北部和山东西北小部分区域是污染最严重的地区，其周围区域污染程度仅次于这一区域，南北方向的污染程度稍高于东西方向。

(3) 从时间分布来看，2004 和2005年两年为污染情况较少的年份，从2006年开始，污染程度急剧加强，2009年稍稍缓和，但从2010年开始又逐渐加重，2013年后呈现下降趋势。河北南部和山东西北小部分区域的污染情况与北京南部、天津和河北中部的污染情况相关性较高。