鲍艳磊 田 冰# 胡引翠 张文静 王晨旭 王天宇

(1.河北师范大学资源与环境科学学院，河北 石家庄 050024；2.河北省环境演变与生态建设实验室，河北 石家庄 050024)

自20世纪90年代以来，雾霾逐渐成为中国严重的环境问题之一[1]。PM2.5是雾霾的主要污染物，许多流行病学研究表明，PM2.5对人体健康产生不良影响，造成呼吸、心血管、免疫和生殖系统出现功能障碍，严重甚至会导致死亡[2-4]。近年来，其引起了国家和社会的广泛关注。目前，PM2.5自动监测依赖于地面监测格网，由于目前中国监测格网数量还有所欠缺，对于研究区域的PM2.5浓度获取存在困难，而利用遥感卫星获得的气溶胶光学厚度(AOD)数据反映近地面污染情况具有连续性好、传输速度快等难以比拟的优势。以AOD为介质的消光系数在垂直方向上的积分，是表征大气浑浊程度的关键物理量，也是确定气溶胶气候效应的重要因素。通常，高的AOD预示着气溶胶纵向积累的增长，直接导致大气能见度的降低。国内外许多学者开展了利用AOD估算PM2.5的研究。研究表明，AOD和PM2.5之间存在明显的正相关关系。近年来，不同学者基于多元线性模型[5]、土地利用回归模型[6-9]、协同时空地理加权回归模型[10]等估算PM2.5，基本得到了很好的结果。由此可知，利用AOD来估算地面PM2.5是可行的，对于了解PM2.5的时空分布及弥补没有监测站点地区的数据缺失是具有参考价值的。

2018年，北京市、天津市、河北全省PM2.5年平均值分别约为51、52、56 μg/m3，京津冀地区平均值超过了《环境空气质量标准》(GB 3095—2012)一级标准(35 μg/m3)。针对该地区大气质量严重污染问题，《大气污染防治行动计划》提出建立京津冀等三大区域大气污染联防联控机制，开启京津冀区域协作治理的新进程。

郝静等[11]使用空间分辨率为10 km的AOD数据，利用混合效应模型对京津冀内陆平原PM2.5浓度进行时空变化定量模拟，得到R2=0.78。景悦等[12]2890利用基于每日的混合效应模型探究AOD与PM2.5的相关关系，R2=0.92。由此说明，高精度的遥感数据和混合效应模型为预测PM2.5浓度提供了更好的数据源和数据模型。但PM2.5与AOD的关系会受到气象因素和其他因素的影响而存在不稳定性。国内外众多学者利用AOD来估算PM2.5浓度时加入气象要素，结果发现，这些因素会在一定程度提高两者相关性[13-15]。薛文博等[16]研究表明，PM2.5与订正过的AOD的R2为0.77；同时引入相对湿度、风速，近地面PM2.5与反演值的R2升至0.79。目前，在混合效应模型中引入气象因子作为自变量对京津冀地区PM2.5与AOD相关关系的研究较少，本研究尝试以AOD数据和气象要素为预测因子，得到最优的混合效应模型，以期达到更好预测京津冀地区PM2.5浓度的效果。

1 研究区概况

京津冀是中国的“首都经济圈”，包括北京市、天津市及河北省11个地级市。京津冀地区包括多种类型的地貌，包括西部的坝上高原、燕山、太行山脉和中部、东南部的平原地区。京津冀地区属于温带季风性气候，具有明显的季节变化特征。燃煤、机动车、工业废气排放多种污染物共生局面加之西、北面环山的地形不利于雾霾的扩散情况导致近年来京津冀地区大气污染超标频度较高。为更好监测PM2.5对大气环境影响，目前京津冀地区有PM2.5监测站点80个、气象站点26个，具体位置如图1所示。

2 数据与方法

2.1 数据来源

2.1.1 中分辨率成像光谱仪(MODIS)AOD数据

MODIS数据通过将MODIS搭载在Terra卫星上获得，光谱范围0.4～14.4 μm。本研究采用美国国家航空航天局(NASA)网站(https：//ladsweb.nascom.nasa.gov)发布的以日为时间分辨率、3 km为空间分辨率的MOD04 L2 C6产品数据，时间为2017年1月1日至12月31日。

首先在ENVI5.1软件中对下载的原始数据进行几何校正、拼接和裁剪等预处理，之后在ArcGIS 10.2软件中获取PM2.5监测站点周围3×3个栅格像元的AOD平均值与对应PM2.5监测站点相匹配。

2.1.2 地面PM2.5数据

京津冀地区的80个PM2.5监测站点的数据来自中国环境监测中心官方网站(http：//106.37.208.233：20035/)。Terra卫星在中国的过境时间是10：30，因此本研究采用卫星过境时段(即10：00—11：00)地面监测站点的PM2.5平均值。在以3 km为单元划分的京津冀地区监测格网中可能会存在一个或多个PM2.5监测站点，此时PM2.5取监测格网中所有监测站点的平均值。

2.1.3 气象和边界数据

京津冀地区26个气象站点的气象数据可从中国气象科学数据共享服务网(http：//data.cma.cn/)获得，包括相对湿度、风速、气温和气压；因气象站点与PM2.5监测站点经纬度不同，因此需要将气象站点的数据进行插值得到与PM2.5监测站点相匹配的气象数据。

省、地、市行政边界数据来源于国家基础地理信息系统网站(http：//ngcc.sbsm.gov.cn/)。

2.2 研究方法

2.2.1 混合效应模型

混合效应模型是指既包含固定效应又包括随机效应的模型。固定效应类似于标准回归系数，直接估计得到。随机效应从它们的方差和协方差估计值中总结而来。随机效应用于解释AOD数据、气象因子与PM2.5之间的日变化关系，以随机截距或随机系数的形式呈现。

景悦等[12]2892-2894指出，基于日变化的混合效应模型能在时间上解释PM2.5和AOD的关系变化，基于站点的混合效应模型也能在一定程度上解释空间上PM2.5和AOD的关系变化，而且两种模型均优于简单的线性回归。式(1)表示在自变量只有AOD的情况下建立的3种模型：(1)站点变化影响下的混合效应模型(基于空间的模型)；(2)时间变化影响下的混合效应模型(基于时间的模型)；(3)PM2.5和AOD关系在时间和站点共同影响下的混合效应模型(基于时空的模型)。3种模型使用相同的数据，基于不同的变量字段来实现。式(2)表示增添了气象因子作为影响因子的3种混合效应模型，反映PM2.5和AOD关系在不同情况下的变化。

图1 京津冀PM2.5监测站点和气象站点分布Fig.1 Distribution of Beijing-Tianjin-Hebei PM2.5 monitoring site and meteorological site

Mij=(α+uj)+(β1+νj)×Aij+Si+εij

(1)

Mij=(α+uj)+(β1+νj)×Aij+(β2+κj)×Tij+(β3+lj)×Wij+(β4+mj)×Rij+β5×Pij+Si+εij

(2)

式中：Mij为第i个监测格网第j天的PM2.5，μg/m3；α为固定截距；uj为第j天随机截距；β1、β2、β3、β4、β5为与AOD、气温、风速、相对湿度、气压有关的固定斜率；vj、κj、lj、mj为第j天与AOD、气温、风速、相对湿度、气压有关的随机效应斜率；Aij为第i个监测格网第j天的AOD；Si为第i个监测格网的随机效应；εij为第i个监测格网第j天的随机误差项；Tij为第i个监测格网第j天的气温，℃；Wij为第i个监测格网第j天的风速，m/s；Rij为第i个监测格网第j天的相对湿度，%；Pij为第i个监测格网第j天的气压，Pa。

2.2.2 模型验证方法

采用十折交叉验证法验证混合效应模型的准确度。具体实现过程：将数据集随机分成10等份，轮流将9份参与建模，剩余一份作为测试集，重复10次直至所有数据都被预测出来，然后将10次验证结果的均值作为对预测模型精度的估计。

模型拟合效果的判断以PM2.5监测值与模拟值的R2、十折交叉验证后相关系数(CVR2)、均方根误差(RMSE，μg/m3)及平均绝对误差(MAE，μg/m3)表示。

3 结果与讨论

3.1 模型拟合

采用基于不同类型的混合效应模型，并尝试在固定效应和随机效应中增添、改变气象因子的类型和数量，来进行模型的拟合和验证，共尝试进行了78种类型的含有不同自变量的混合效应模型。

在固定效应和随机效应中加入不同气象因子后，R2、CVR2、RMSE、MAE的取值见表1。R2和CVR2取值越大，两者的相关性越强；RMSE和MAE取值越小，模型的误差越小。综合分析，基于空间的模型的R2、CVR2比另外两种模型小，基于时间的模型的R2、CVR2基本都比基于时空的模型小；基于空间的模型的RMSE、MAE明显大于另外两种模型，后两种模型数值有交叉，但总体基于时空的模型的RMSE、MAE小于基于时间的模型。

在78种模型中，分别选择基于空间、时间和时空的最优模型进行分析。基于空间的模型，在固定效应和随机效应中同时含有AOD、风速和相对湿度的模型得到的结果最优；基于时间的模型，在固定效应和随机效应中同时含有AOD、风速、相对湿度、气温和气压的模型的结果最优；基于时空模型，在固定效应和随机效应中同时含有AOD、风速、相对湿度和气温，且仅在固定效应中含有气压的模型的结果最优。将这3种最优模型经十折交叉模型校正，结果见图2。模型R2=0.90，校正后基于时空最优模型拟合程度最好，CVR2=0.81，RMSE、MAE分别为13.44、10.12 μg/m3。表1中，基于时空的模型的CVR2比R2小，说明模型存在部分过度拟合，但总体来说模型的显著性明显，说明气象因子可提高模型的拟合精度和PM2.5的预测精度。

表1 R2、CVR2、RMSE、MAE的取值

图2 模型模拟值与监测值间校正后得到的相关关系Fig.2 Correlation diagram between predicted and actual values of the model

基于时空的模型中引入气温、相对湿度、风速和气压4种气象因子，此模型得到的固定效应(包括固定截距和固定斜率)的解见表2。其中，混合效应模型的固定效应类似于标准回归系数，表示全年内自变量对于PM2.5的固定影响，AOD固定斜率估计值为21.861 6，说明PM2.5与AOD呈现正相关；相对湿度固定斜率估计值0.371 0，说明相对湿度越大越有利于大气颗粒物质的形成，PM2.5浓度越大；气温固定斜率斜率估计值为0.064 2，说明PM2.5与气温呈现正相关，气温通过控制边界层高度影响大气颗粒物质的扩散；气压固定斜率估计值为0.000 1，说明气压增大，PM2.5浓度也提高，但影响程度较小。风速固定斜率估计值为-0.453 2，表明PM2.5浓度与风速呈反比，风有利于大气颗粒物的扩散，风速越大越有利于空气中的PM2.5的输送与扩散。风速与相对湿度的系数较高，是相对主要的气象影响因子。

由表3可知，各城市的CVR2为0.71～0.88，RMSE为9.79～14.92 μg/m3，MAE为8.19～11.80 μg/m3。根据CVR2，北京市、邢台市、秦皇岛市、张家口市、廊坊市和保定市的拟合效果较好，石家庄市、承德市的拟合效果较差，但总体能在一定程度上表达PM2.5与AOD的关系，表现为较高的拟合效果、较低的误差。

3.2 京津冀地区PM2.5浓度模型适应性与时空分布分析

2017年，京津冀地区PM2.5监测值为1～232 μg/m3，平均为45.23 μg/m3；PM2.5模拟值为0～229 μg/m3，平均为44.60 μg/m3。分别将各监测站点的PM2.5监测、模拟年均值插值到整个研究区，结果见图3。PM2.5浓度空间差异性和数值吻合度均较高，说明模型的模拟效果较好，因此可利用该模型对缺乏监测数据地区以及大范围研究区进行PM2.5的模拟与预测。

为方便分析，本研究将京津冀分为4个地区，分别为北京市、天津市、河北平原区(廊坊市、保定市、石家庄市、衡水市、邢台市和邯郸市)和河北非平原区(唐山市、秦皇岛市、沧州市、张家口市和承德市)。2017年，北京市、天津市、河北非平原区、河北平原区的PM2.5模拟年均值分别为42.51、51.38、38.32、46.72 μg/m3，整体上北部低于南部、东部沿海低于西部地区。

本研究还将一年分为冷暖两季(冷季为10月至次年3月，暖季为4—9月)，对PM2.5浓度及模型适应性进行分析，结果见图4。京津冀全区冷、暖季PM2.5模拟平均值分别为44.99、44.49 μg/m3。全区暖、冷季PM2.5模拟值与监测值分别相差0.31、0.36 μg/m3，其中北京市分别相差0、1.76 μg/m3，天津市分别相差0.17、2.00 μg/m3，河北平原区分别相差0.35、1.59 μg/m3，河北非平原区分别相差0.83、1.10 μg/m3。全区平均值相对误差值不超过0.008，分区平均值相对误差值为0～0.037。总体来说，暖季相较于冷季估计均值更精确，模型对PM2.5浓度估计精度较高。

表2 固定效应的解

表3 基于城市统计的混合效应模型拟合结果及校正结果

图3 2017年PM2.5监测与模拟年均值分布对比Fig.3 Comparison of the PM2.5 distribution of monitoring values and simulation values in 2017

图4 基于不同地区的PM2.5监测值与模拟值对比Fig.4 Comparison of PM2.5 monitoring values and simulated values in different regions

研究区4个地区PM2.5模拟年均值均大于35 μg/m3，但是不同地区也存在差别。河北非平原区中的张家口市、承德市PM2.5年均值较低，分析原因可能是因为当地工厂相对较少，空气流通状况好；秦皇岛市为沿海城市，常年有海风的作用，因此PM2.5年均值相对较低，空气质量较好。河北平原区工厂多、人口密集，产生PM2.5的源头多；冷季受西北风影响，随风飘荡的污染物在遇到太行山时无法翻越，造成累积，然后便有可能形成气团污染到其他周边城市，造成二次及三次污染，导致内部平原地区PM2.5年均值明显偏高。北京市人口密集，但工厂外迁，冬季盛行西北风，湿度较小，较有利于PM2.5的扩散，所以PM2.5年均值相对不是很高。整体来说，京津冀地区呈现南高北低、平原地区东南高西北低、非平原地区高纬度区域相对低的空间格局，京津冀地区空气质量亟待解决。

4 结 语

(1) MODIS 3 km AOD的数据全球开放、获取方便，因此本研究使用此数据来预测地面的PM2.5浓度。目前，3 km的AOD产品可通过MODIS的暗像元算法(DT)和深蓝算法(DB)得到，但通过DB得到的产品没有发布，应用DT在冷季云、雾及重污染天气反演得到的AOD产品质量较差，导致部分冬季天气恶劣地区参与建模的数据较少。在以后的研究中可尝试将3、10 km数据进行融合或采用SARA算法反演出500 m分辨率的AOD产品填补AOD的缺失，以提高研究的精度。另外，以后的研究中还可加入县级监测站点，增加站点数量与均匀程度进一步提高模型的精度。

(2) 局地气象是影响PM2.5浓度的重要因子，每天每个地方的气象因子存在差异，因此可通过最优插值得到监测站点的各种气象因子的数值。在构建模型时加入不同的气象因子，所得的CVR2存在不同。其中，风速和相对湿度的影响系数较大。

(3) 基于时空的模型拟合结果比基于空间、时间的模型优。

(4) 整体来说，京津冀地区呈现南高北低、平原地区东南高西北低、非平原地区高纬度区域相对低的空间格局，京津冀地区空气质量亟待解决。