地球静止卫星和网格化站点支持下的PM2.5精细制图

范东浩，秦凯，杜娟，何秦，辛世纪，刘鼎医

1.中国矿业大学环境与测绘学院,徐州221116;

2.中国科学院空天信息创新研究院,北京100094;

3.江苏省徐州环境监测中心,徐州221002

1 引言

PM2.5因其粒径小、面积大、活性强、易附毒害物质且能在大气中停留较长时间、输送较远距离等特点（Saffari 等，2014），对人体健康和大气环境质量有很大危害，是中国主要的空气污染源之一（Cao等，2012）。

基于卫星遥感数据估算近地面颗粒物的方法有比例因子法或化学传输模式法（Cameron 等，2017）、物理经验法（张莹和李正强，2013）及统计模型法等（Zhang 等，2021）。统计模型法又包括Wang 和Christopher（2003）采用的简单线性回归方法、经研究发现预测水平显著高于简单线性回归的多元线性回归模型方法（Koelemeijer 等，2006；Wu 等，2012；潘锦秀等，2019），地理加权回归模型（GWR）（Hu，2009），地理时空加权回归方法（GTWR）（Huang 等，2010），神经网络方法（Ordieres等，2005；郭建平等，2013；Chen等，2020）等方法。在卫星选用上，早期多选用MODIS 传感器数据，目前日本葵花8 号卫星（Himawari-8/AHI）及韩国千里眼卫星（COMS/GOCI）的卫星遥感数据正广泛被应用在近地面颗粒物遥感估算（Xu 等，2015；Wang 等，2017；Tang等，2019）。通常卫星AOD产品是经过卫星表观反射率（TOA）反演而来（Hsu 等，2004；Guo等，2009）。基于此，可以直接建立TOA 产品和地面站点监测的PM2.5浓度的反演模型（Shen 等，2018；Fan等，2021；Yin等，2021）。对于PM2.5选取目前多采用中国环境监测总站CNEMC（China Environmental Monitoring Center）的国控站点数据（Zou 等，2017；Wang 等，2018；Tang 等，2019；Bai 等，2019a）。当前，许多地方政府完成了网格化空气质量监测网络建设，为精细化监管城市空气质量奠定了基础。本文旨在检验网格化地面监测站点支持下的卫星遥感PM2.5精细化制图能力，以徐州为例，以网格化地面监测站点和Himawari-8/AHI 及COMS/GOCI 卫星数据为主，辅以气象、人口、道路网密度等数据，开展对比实验研究。

2 研究区域和数据

2.1 研究区域概况

本文研究区域为江苏省徐州市（33°43′N—34°58′N、116°22′E—118°40′E），是淮海经济区中心城市，下辖2 市（新沂市、邳州市）、3 县（丰县、沛县、睢宁县）、5 区（云龙区、鼓楼区、泉山区、铜山区、贾汪区）。徐州周边有大型煤炭或煤电基地及大量化工企业，主要工业区位于贾汪区内（樊文智等，2018；张宇静等，2019）。徐州市市区（铜山区、泉山区、云龙区和鼓楼区）是徐州市整体各类大气污染物浓度较高的区域。除城市排放因素外，还与徐州市区四面环山的盆地效应有关。

2.2 Himawari-8/AHI

日本宇航局JAXA（Japan Aerospace Exploration Agency）于2014-10-07 发射了新一代静止卫星“葵花8 号”（Himawari-8），并于2015-07 开始业务化运行。Himawari-8 卫星具有良好的时间分辨率，可以实现10 min/次的高频观测，其上搭载的AHI 传感器拥有16 个波段，可见光波段范围内空间分辨率最高可达0.005°。

本文使用了逐小时的Himawari-8/AHI L3 级TOA（表观反射率）数据，时间阶段为2019-07—2020-06 期间每天的8：00—15：00（为方便数据处理，选取与COMS/GOCI一致的时间序列），空间分辨率最高为0.005°（Bessho 等，2016），包括蓝光（1 波段，470 nm，0.01°分辨率）、绿光（2 波段，510 nm，0.01°分辨率）、红光（3 波段，640 nm，0.005°分辨率）及短波红外波段（6波段，2300 nm，0.02°分辨率）。此外，使用红光和近红外波段（4波段，860 nm，0.01°分辨率）计算了归一化植被指数NDVI（Normalized Difference Vegetation Index），还用到了相关角度数据等。

对于表观反射率的处理，参照Yan 等（2020）中将AHI 传感器的反照率数据转为TOA 数据的方法，如式（1）所示：

式中，ρ*为表观反射率，θ0为太阳天顶角，λalbdeo为反照率。

2.3 COMS/GOCI卫星数据

韩国的通信海洋气象卫星，或称千里眼卫星COMS（Communication Ocean Meteorological Satellite）于2010年6月发射，定轨于128.2°E，是一颗静止轨道多用途卫星。COMS 卫星主要搭载3 种有效载荷：Ka 频段通信仪、地球同步海洋水色成像仪GOCI（Geostationary Ocean Color Imager）和气象成像仪。

本文所用为COMS 卫星的GOCI 传感器数据，以（36°N，130°E），为观测中心，陆地范围包括朝鲜半岛、中国东部和日本，空间分辨率最高可达0.005°，时间分辨率为1 h。GOCI 每天可传输北京时间8：00—15：00共8幅图像（Ryu等，2012）。

本文采用GOCI 的L1B 波段数据作为建模数据，共使用8 个光谱通道进行气溶胶反演，前6 个为可见光波段，而后2 个为近红外波段，并参照Choi 等（2016）所给出的如式（2），对L1B 数据进行预处理。

式中，λ为上述的8个光谱波段波长，L(λ)为GOCI的辐射高度，μ0为观测时太阳天顶角，E0(λ)为太阳辐射能流。

此外还采用YEAR-V2 版本（Choi 等，2016）的GOCI-AOD 数据，其空间分辨率较粗仅为6 km，因此只作为掩膜的判别数据，不作为建模的特征参数。

2.4 PM2.5浓度地面观测数据

徐州于2018年建成了覆盖全市的空气质量网格化监测网络，本文使用其中的172个站点（图1）的逐小时观测数据。考虑到逐小时PM2，5浓度变化较大，可能存在随机偏差，参照原环境保护部发布的《灰霾污染日判别标准》，筛选PM2.5浓度在0—500 μg/m3内的数据以减小站点数据漂移带来的影响。为了与卫星遥感数据进行时空匹配，以0.005°空间分辨率构建网格，同一网格内的不同站点求均值。最终，形成了覆盖徐州市的165 个网格，所用质控合格率为99%，数据有效率为91%（数据有效率=采用的数据数/实际上传的数据数；数据质控合格率=采用的数据数/应上传的总数）。此外，本文还使用了7 个国控站点的逐小时PM2.5数据（http：//106.37.208.233：20035/［2021-07-19］）作为对比，分别是黄河新村、淮塔、新城区、桃园路、农科院、鼓楼区政府和铜山区招生办站点。

图1 徐州地区CNMEC国控站点及网格化监测站点的空间分布Fig.1 Spatial distribution of CNMEC’s national control stations and grid monitoring stations in Xuzhou

2.5 其他数据集

为了获得大范围的气象数据，本文使用ECMWF（欧洲中期天气预报中心）的ERA5 气象数据（Fifth Generation Atmospheric Reanalysis）（Hersbach等，2020），其时间分辨率为逐小时，空间分辨率为0.25°。本文主要使用了9 个气象参数，即地表气压、总降雨量、2 m 露点温度、2 m 气温、10 m经向风、10 m 纬向风、边界层高度、高植被叶面积指数和低植被叶面积指数。

人口数据来自于世界人口网站WorldPop（https：//www.worldpop.org［2021-07-19］），其通过人口普查及top-down的估计方法得到100 m空间分辨率的网格化数据（Stevens 等，2015；Lloyd 等，2019），本文选用2019年及2020年年均数据。土地利用数据来源于清华大学发布的FROM-GLC 分类模型（Finer Resolution Observation and Monitoring of Global Land Cover）（Gong 等，2013）。道路网密度数据来源于OpenStreetMap 网站（https：//www.openstreetmap.org［2021-07-19］）。

2.6 研究方法

本文使用了两种集成式机器学习方法，分别为Boosting的极端梯度提升XGBoost和Bagging的随机森林，及基于核密度估计和加权最小二乘法的一种局部回归方法—地理时空加权回归（GTWR）算法。

XGBoost 在梯度提升决策树GBDT（Gradient Boosting Decision Tree）的基础上改进得到，GBDT的核心在于每一次迭代都是向残差降低的梯度方向继续建立新的基础学习器，以此将弱分类器转为强分类器。XGBoost 在GBDT 基础上对损失函数进行二阶泰勒展开，通过正则化处理和并行计算，来避免过拟合（Chen和Guestrin，2016）。

随机森林RF（Random Forest），已广泛用于解决回归和分类问题（Bai等，2019b；Wu等，2020），其特点是自展法随机采样（Bootstrap Sampling）。通过不断采样得到采样集，基于每个采样集训练出一个基础学习器，再将这些基础学习器通过结合形成强学习器。

时空地理加权回归GTMR（Geographically and Temporally Weight Regression）（Huang 等，2010）在地理加权回归模型（GWR）的基础上又加入了时间的影响，可以解决时空非平稳性问题，能有效地解决回归分析中无法同时考虑空间和时间特征的问题，可用于PM2.5浓度预测中（Li等，2020）。

技术路线如图2所示，将卫星的TOA 数据、气象数据、土地利用类型数据与PM2.5站点数据等进行时空匹配后（空间匹配选用0.005°×0.005°空间分辨率网格）进行数据清洗，再划分训练集和测试集通过XGBoost、RF和GTWR 这3种方法进行模型训练，按交叉验证结果择优选为模型预测的主要方法，并将估算数据通过GOCI-AOD 掩膜得到最终的估算结果（即将AOD 数据缺失地区的PM2.5数据设零值从而减弱因云雾引起的干扰），并与CHAP和TAP数据集进行对比分析。

图2 技术路线图Fig.2 Flowchart of the methodology

3 结果与分析

3.1 模型选取和参数优化

根据徐州地区的大致经纬度可分为280×480个空间分辨率为0.005°×0.005°的网格，并将卫星遥感数据等通过最临近插值方法填入预定义网格，经过时空匹配后，最终得到了37万行样本数据。

由图3 可知，XGBoost 模型和RF 模型的R2及RMSE 要优于GTWR 模型，因为GTWR 方法中随多特征参数误差的累积，整体验证结果会快速变差；XGBoost 模型的R2及RMSE 在测试集和训练集都具有一致性，且离散点分布较为均匀；而RF模型的训练集拟合程度明显比测试集好，可能存在过拟合情况。因此，选用XGBoost 模型开展后续研究。

图3 XGBoost、RF和GTWR 3种方法的测试集及训练集结果对比Fig.3 Comparison of test set and train set results of XGBoost，RF and GTWR

考虑到各个参数在模型中重要程度不同且可能相互影响，选用不同的参数组合进行XGBoost模型优化，按照数据来源分为11 种组合（表1），其中H8 代表选用葵花8 号（Himawari-8）卫星数据，GOCI 代表选用GOCI 载荷数据，GLC 代表选用FROM-GLC 土地利用数据，ERA5 代表选用ERA5气象数据，下同。结果表明（表1），H8 与ERA5结合的模型综合性能最优，而不包含GLC 的组合优于包含GLC 的组合，如H8+ERA5 组合优于H8+ERA5+GOCI 组合，其原因可能是GLC 产品分辨率较粗，无法反映PM2.5精细化制图所需要的土地利用类型细节变化。

表1 不同特征参数组合的PM2.5浓度反演结果比较Table 1 Statistics of PM2.5 retrieval uncertainties based on various combination of input data

进一步，选取性能靠前的3 组参数组合，即H8+ERA5，GOCI+ERA5，H8+GOCI+ERA5，与国控站点开展时间序列对比分析。考虑到GOCI 载荷仅观测获得北京时间8：00—15：00的表观反射率数据，所以采用该时段的均值进行对比。图4 中，XGB_HGE 代表使用XGBoost 算法并选用H8+GOCI+ERA5 数据，XGB_HE 代表使用XGBoost 算法并选用H8+ERA5 数据，XGB_GE 代表使用XGBoost 算法并选用GOCI+ERA5 数据，下同。由图4 可知，H8+GOCI+ERA5 比另外两者的误差更小，因此作为最终的制图参数组合。

图4 3种特征参数组合PM2.5预测结果与国控站点对比图Fig.4 Comparison diagram of PM2.5 prediction results of three characteristic parameters combined with national station measurements

3.2 制图结果分析

通过与网格化站点及国控站点进行对比，如图5 和图6所示，在年均和月均时间尺度的逐小时制图中，本文采用模型的精细化成图结果与站点数据基本匹配，能较好地反映出不同时间不同地区的PM2.5浓度变化趋势。

图5 网格化监测站点和国控站点数据与本文PM2.5预测结果9：00、12：00和15：00等3个时次的对比（年均）Fig.5 Spatial distribution comparison between the annual averaged model derived PM2.5 concentrations and those that are measured from national station measurements

图6 网格化站点与本文PM2.5预测结果逐小时对比（月均）Fig.6 Comparison between the prediction Results of PM2.5and data from National Station Measurements（monthly average）

3.3 与TAP和CHAP数据集对比

为了检验本文研究中采用了网格化站点数据后的PM2.5精细化制图能力，进一步与仅使用国控站点的的清华大学TAP（Xiao等，2021）和马里兰大学的CHAP 数据集（Wei 等，2021）进行对比分析。TAP 全称为Tracking Air Pollution in China，即中国大气成分近实时追踪数据集，该数据集融合地面观测、卫星遥感、排放清单和模式模拟等多源数据，构建多尺度、近实时的中国大气气溶胶和气态污染物浓度数据集。目前TAP 数据集提供2000年至今的中国PM2.5浓度数据，空间分辨率为10 km，可下载数据的时间分辨率为每日（http：//tapdata.org/［2020-07-19］）。CHAP 全称为China HighAirPollutants，基于多源卫星遥感和人工智能技术，考虑空气污染的时空异质性，结合地基测量、遥感产品、大气再分析和模式模拟等资料生产得到。数据集包含PM1、PM2.5、PM10等7 种主要空气污染物。本文选取由MODIS 卫星AOD 数据估算的空间分辨率为1 km 的每日PM2.5数据产品（备注：CHAP的逐小时数据仅提供2018年全年）。

考虑到年均数据无法体现制图细节，逐小时数据可能存在缺失，主要从月均和日均两个尺度检验采用网格化站点数据后对于PM2.5精细化制图的贡献。图7 为月均值尺度下的本文结果和TAP、CHAP 数据集与地面检测站点的对比。其中XGB_HGE为11种特征参数组合选出的性能最优的XGBoost 模型。因为TAP 和CHAP 所用月均值数据均为每日24 h 的月均，而本文模型所采用的为每日8—15 h 的月均，国控站点数据做相应两种月均值对照。可以看到除部分月份的部分站点外，本文所采用模型基本与TAP 及CHAP 数据保持一致，在部分月份与地面观测拟合程度更好。

图7 不同模型预测PM2.5结果与国控站点数据对比Fig.7 Comparison of prediction results of PM2.5 by different models with data from national station measurements

如图8所示，2019-09-09，本文结果很好地表征了市区（铜山区、泉山区、云龙区和鼓楼区）和工业区（贾汪区）的PM2.5污染情况（＞100 μg/m3），与地面观测一致；TAP数据虽然显示了PM2.5空间分布的东西差异性，但是PM2.5浓度都低于100 μg/m3，无法表征市区和工业区的高值；而CHAP数据因为使用MODIS 数据，每天过境只有两次，数据缺失较多。如图9所示，2020-03-18，本文结果与地面观测一致，而TAP和CHAP数据均无法表征市区和工业的污染情况。随机抽取的两天数据表明，本文结果空间覆盖度更好，且更能反映出徐州市区和周边区县的PM2.5空间分布差异性。

图8 本研究的单日PM2.5空间分布图与仅使用国控站点的结果对比（2019-09-09）Fig.8 Comparison between the Daily Spatial Distribution of PM2.5 Prediction in this Research and the Results Only Using National Station Measurements（2019-09-09）

图9 本研究的单日PM2.5空间分布图与仅使用国控站点的结果对比（2020-03-18）Fig.9 Comparison between the Daily Spatial Distribution of PM2.5Prediction in this Research and the Results Only Using National Station Measurements（2020-03-18）

4 结论

本文主要使用Himawari-8/AHI 和COMS/GOCI地球静止卫星的表观反射率和徐州地区的空气质量网格化监测数据，通过机器学习和时空地理加权回归方法开展了0.005°分辨率的PM2.5浓度精细化制图研究。通过与国控站点和现有数据集的对比分析，得出以下结论：

（1）地球静止卫星遥感数据因其具有更高的时间分辨率，在模型预测中可以更好地契合逐小时的城市网格化监测站点数据，比极轨卫星数据更适合城市PM2.5精细化制图。

（2）基于城市网格化监测站点及卫星遥感数据的估算结果，在一定程度上可以弥补国控站点稀疏的不足，能够更好地反映出城市内不同区域的PM2.5浓度分布差异性，可以更好地服务于城市空气质量精准管控。

志 谢本研究使用了由日本JAXA 提供的Himawari-8/AHI（https：//www.eorc.jaxa.jp/ptree/）以及韩国KOSC（韩国海洋科学技术研究院）提供的COMS/GOCI（http：//kosc.kiost.ac.kr/）卫星数据。此外，本文研究还获得美国马里兰大学韦晶博士（https：//weijing-rs.github.io）和韩国延世大学Kim Jhoon 教授团队（http：//atrad.yonsei.ac.krl）所提供的相关数据，在此一并表示感谢。