1. 引言
我国经济水平发展迅速,城镇化水平不断提升,人们生活水平日益提高,包括SO2、CO2、粉尘等污染物的排放也逐渐增多,加重了我国大气生态环境污染情况,导致霾污染天气频发,雾和霾均属于大气气溶胶系统的组成部分,但是霾的核心是一些粉尘微粒,是空气动力学直径为2.5微米或更小的微粒(PM2.5)大量聚积的过程[1][2]。由于PM2.5直径极其微小,导致其可以通过人的呼吸进入人们的呼吸道以及肺泡中,在空气中携带的有毒物质和重金属污染物也会附着在这些细小颗粒物上随之进入人体参与血液循环,对人体造成严重损伤[3][4]。对PM2.5的浓度进行估算是有重要意义的研究。
PM2.5主要通过地面空气质量监测网站获取,由于我国监测网络运行较晚,自2013年才开始运行,且地面监测站点分布不均,在一定程度上影响了对PM2.5污染影响程度评估的能力[5],无法支撑我国大范围的PM2.5浓度监测需求,所以越来越多的学者选择使用气溶胶光学厚度来进行PM2.5浓度估算,能有效弥补地面站点不足的缺点[6][7],气溶胶光学厚度(Aerosol Optical Depth, AOD)为介质的消光系数在垂直方向上的积分,用来描述气溶胶对光的削减作用,是气溶胶最重要的参数之一,AOD是一个无量纲正值。卫星遥感AOD数据可通过不同的数据模型进行PM2.5的浓度估算,证明了PM2.5与AOD之间会随着时间和空间的变化而相应产生变化[8],卫星遥感气溶胶光学厚度已经被广泛应用于PM2.5的浓度估算中[9][10],在吴宇宏等[11]的研究中,AOD因子与PM2.5的相关因系数高达0.65,在吴迪[12]的研究中,AOD因子与PM2.5浓度的相关系数为0.50。通过卫星获取的AOD有着大范围高精度的优点,但是当遇到多云、多雨或者雪覆盖导致地面反射率增高时,卫星遥感获取的AOD数据会存在大范围的缺失。缺失的AOD数据对PM2.5浓度估算精度有一定的影响,所以在对PM2.5进行估算时,对AOD数据进行补值是必要的处理。
传统AOD补值方法包括克里金插值、泛克里金插值和反距离加权等[13],在研究PM2.5-AOD的时空关系中,常用的统计方法有线性回归模型、多元线性回归模型[14][15]、线性混合效应(LME)模型[16]、广义加法模型(GAM)[17]和地理加权回归模型(GWR)[18][19]等,KIM[15]等人利用MODIS AOD数据和部分ERA5气象再分析数据建立了一个多元线性回归(MLR)模型对PM2.5浓度进行估算,总体相关系数大于0.8。这些实验中使用的AOD补值方法均为传统补值方法,对AOD的时空分布会产生一定的影响,进而影响PM2.5浓度估算的精度。为进一步提高精度,有学者在混合效应模型中加入与PM2.5相关的时间和空间辅助变量,包括气温、降水、REA5气象再分析资料、土地利用数据和人口数据等变量,相较于普通的线性回归模型精度有所提高。也有学者用不同的方法对AOD进行补值得到更好的PM2.5浓度估算结果,Liu[20]等人将MERRA-2 AOD进行重采样处理为Himawari-8 AOD进行数据填补,并将填补好的Himawari-8 AOD和NDVI、路网、人口密度和坐标等作为预测因子使用随机森林模型对地面PM2.5进行浓度估算,结果显示R2为0.88。机器学习也常用于PM2.5的浓度估算中,Li[21]等人利用MAIAC AOD数据,结合其他时空预测因子和空间自相关性构建了一个每日混合效应空间模型,利用嵌入的结构化和非结构化空间随机效应来解释空间自相关性。然后基于自举聚合(Catboost)对基本模型的点估计进行平均,以减少预测中的方差。然后,开发了约束优化,对PM2.5浓度的完整时间序列进行了预测。
在PM2.5浓度估算的众多常用变量中,AOD与PM2.5浓度的相关性较强,但卫星遥感AOD数据存在大面积缺失,传统的AOD补值方法只顾及了该数据的空间相关性或者时空相关性,没有综合考虑两者的互相关关系,对补值结果会产生一定的影响,从而影响PM2.5浓度估算的精度。因此,本文提出一种新的AOD补值方法对MAIAC AOD进行补值,并将AOD补值结果结合ERA5气象再分析数据使用Catboost算法进行PM2.5浓度进行估算。
2. 研究概况
2.1. 研究区概况
本次研究区域为中华人民共和国的陆地部分,中国陆地总面积约为960万平方千米。目前中国有34个省级行政区,包括23个省、5个自治区、4个直辖市和两个特别行政区(https://www.gov.cn/guoqing/index.htm)。按行政区主要划分为七个大区,分别是东北地区、华北地区、华中地区、华南地区、华东地区、西北地区和西南地区,研究区域如图1所示。
Figure 1.Overview of the study area and distribution of PM2.5stations
图1.研究区概况及PM2.5站点分布
2.2. 数据来源及处理
本研究使用数据包括MAIAC AOD数据、AERONET AOD数据、ERA5气象再分析数据和PM2.5站点数据。其详细信息及来源介绍如下。
2.2.1. MAIAC AOD数据
美国国家宇航局(NASA)分别于1999年和2002年发射了Terra卫星和Aqua卫星,搭载于两个卫星上的MODIS传感器提供了AOD数据[22],多角度大气校正算法(Multiangle Implementation of Atmospheric Correct, MAIAC) AOD是Lyapustin等人发布的全球覆盖高空间分辨率(1 km)每日AOD数据集(MCD19A2),MAIAC AOD属于MODIS的L2A级产品[23]。MCD19A2 Version 6的数据在2023年7月31日之后弃用,官方更新了MCD19A2 Version 6.1数据[24],该数据将AOD值的有效范围从第六版本的0-3更新到了现在的0-6,数据可以从NASA官网网站进行获取(https://ladsweb.modaps.eosdis.nasa.gov/)。根据最新的MAIAC数据操作指南对该数据进行下载中国陆地区域2016年3月~2021年2月的AOD数据文件进行预处理并质量控制,选取AOD_QA为“0000”代表Best quality的数据并将MCD19A2中Terra和Aqua卫星的Optical_Depth_055波段数据融合作为实验数据。
2.2.2. AERONET AOD数据
AERONET AOD常作为卫星遥感AOD的验证数据,可以从全球布站的气溶胶特性地基观测网(http://aeronet.gsfc.nasa.gov/)进行获取。MODIS过境时间为上午10点30分(Terra卫星过境时间)和下午13点30分(Aqua卫星过境时间),以AERONET站点为中心采用50 km × 50 km的空间窗口计算卫星的AOD空间平均值,匹配地面观测点时间 ± 30 min的AERONET AOD平均值[25],以此进行验证分析。
MAIAC AOD数据是550 nm波段处的数据,AERONET站点提供的AOD数据并不包括550 nm波段的AOD数据,因此,需要通过AERONET AOD 440 nm和675 nm两个波段的AOD插值出550 nm波段的AOD值[26],波长计算公式如下:
(2.1)
(2.2)
式中
和
为波长,
为波长
和
之间的波长指数,
、
和
分别对应波长为
、
和
时的AOD数据。
2.2.3. ERA5气象再分析数据
本次研究所用的ERA5再分析数据可以从欧洲中期天气预报中心(ECMWF)获取(https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era5),下载2016~2020年中国陆地区域的气象数据,包括温度(TEMP)、风速(WS)、边界层高度(BLH)和相对湿度(RH),气象数据的空间分辨率为25 km。
现有研究表明,温度、风速、边界层高度和相对湿度等气象数据与PM2.5有高相关性,作为预测PM2.5的辅助变量可以有效提高预测精度[27]。地面气象监测点虽然可以提供高精度的气象数据,但是由于地面站点的空间分布有限,难以提供大面积的气象监测数据,且站点可供下载的气象参数较少加上时间跨度较大,相对而言,选择气象再分析数据能提供更多的参数选择和较高的时空分辨率和时间跨度[28]。
ERA5气象再分析数据的风速(WS)变量可以由10 m风速u变量和10 m风速v分量计算得到。计算公式如下:
(2.3)
式中WS为最终使用的风速变量;u10为10 m风速u变量;v10为10 m风速v变量。
2.2.4. PM2.5站点监测数据
PM2.5的数据可以从中国环境监测总站的全国城市空气质量实时发布平台(https://air.cnemc.cn:18007/)进行获取,本次研究下载了2016~2020年的PM2.5地面监测数据。在研究时间范围内的PM2.5监测站点有1500个左右,相对于本文的研究区域来说,站点数据冗余过多,经过筛选以及剔除(数据缺失和异常)多余站点后剩余138个站点。PM2.5监测站点分布见图1。
3. 分析方法
3.1. MAIAC AOD时空补值模型
AOD数据是具有时空属性的时空数据,在进行AOD数据的补值时,需要同时兼顾时间和空间的影响。本文对MAIAC AOD数据分别使用Prophet-LSTM时序组合模型进行时间维度补值和P-Bshade模型进行空间维度补值,然后将两个维度的补值结果进行线性融合,最终得到顾及时空影响的AOD补值结果。总流程如图2所示。
Figure 2.Flow chart of the spatiotemporal compensation model
图2.时空补值模型流程图
按下式将时间维和空间维补值结果进行线性融合。
式中,A表示空间维度权重;B表示时间维度权重;
表示最终补值结果;
表示空间维补值结果;
表示时间维补值结果。本文所用的时间维和空间维的插值方法中,涉及到计算空间或时间协方差的影响,顾及到时空平衡性,本文的A、B取值均为0.5。
3.1.1. Prophet-LSTM时间序列补值组合模型
Prophet是一种基于加法模型的时间序列数据预测,有具体的数学模型,能快速地进行时间预测,在建模过程中考虑了趋势线、季节性、周期性,以及外生变量等因素的影响,预测效果好,相对于传统时序模型有很大优势。Prophet对于异常值、丢失的数据具有健壮性,可以对杂乱的数据进行合理的预测[29]。LSTM模型是一种特殊的循环神经网络(RNN)模型,LSTM模型的“门”机制使得信息可以在时间序列中正确地流动,“门”机制可以限制流量信息,从而使得梯度在反向传播中不会消失或者爆炸。LSTM模型对不同类型的数据有很强的适应性,可以捕获长时间序列数据的非线性关系[30][31],LSTM让循环神经网络具备更强更好的记忆性能,可以有效地处理较长的时间序列数据。
Prophet-LSTM时间序列补值组合模型由Prophet模型和LSTM模型两个部分组成,Prophet模型负责为LSTM模型提供完整的AOD时间序列,LSTM模型负责对AOD数据进行时间维度补值。具体模型构建如图3。
Figure3.Prophet-LSTM combined timing model
图3.Prophet-LSTM组合时序模型
Prophet-LSTM算法流程描述:
① 数据预处理:选择用Prophet补齐MAIAC AOD数据的时间序列,并用mask标记缺失值。对于数据中的缺失值,将其替换为NaN (Not a Number),以便在后续步骤中进行处理。同时,创建掩码来标记原始序列中的缺失值位置,将缺失值位置的掩码设置为1,非缺失值位置的掩码设置为0。
② 为了进行训练,将完整的序列按照长度为64的窗口进行循环拆分,并将每个窗口作为一个样本输入到模型中,保证模型可以对时间序列的不同部分进行学习和预测,同时也便于训练和批量处理。在拆分序列的过程中,只使用已有的数据作为输入,补齐的数据仅用于填充缺失部分,不参与损失函数的计算。这样可以确保模型在训练时只利用真实的数据进行学习,而不会受到补齐数据的影响。
③ 在LSTM模型中,为了提高模型表达能力,在LSTM网络中堆叠了4个LSTM层,每个LSTM层都具有相同的隐藏状态大小(hidden_size),以确保信息的传递和记忆能力。并且将LSTM中的普通卷积换成了CSPConv并添加空间通道注意力(Spatial Channel Attention),空间通道注意力是一种自适应地调整通道权重的方法,可以使模型更关注重要的特征通道。在CSPConv Block的最后一个残差块作为空间通道注意力模块,以增强模型的特征表达能力。
④ 为了避免过拟合问题,在每个LSTM层之间添加了一个Dropout层,Dropout层可以在训练过程中随机丢弃一部分神经元的输出,这样可以减轻网络对某些局部特征的依赖,减少过拟合风险。
⑤ 在计算损失函数时,使用带有掩码(mask)的损失函数,只计算非缺失部分的损失,对于每个时间切片,只计算掩码为0部分的损失。对于每个时间切片,根据掩码来选择是否计算该时间切片的损失。以预测序列与真实序列之间的均方差(Mean Squared Error, MSE)作为损失函数。使用带有掩码的损失函数时,只计算非缺失部分的预测值和真实值之间的方差。
⑥ 在训练阶段,通过反向传播和优化算法对模型进行训练,直到模型收敛并达到最佳性能。训练收敛后,使用构建的网络,为序列中缺失的部分预测其数值。对于连续缺失的情况,我们逐步迭代地预测序列的缺失部分,先预测第一个缺失值,然后将其用于下一个缺失值的预测,以此类推。
⑦ 模型在训练的时候同时训练正反方向的序列,避免序列头的缺失;对于中段缺失的部分,使用正反模型预测值的均值;对于头端的缺失,使用对侧方向的预测值。训练完成后,我们可以使用已经训练好的网络来预测序列中缺失部分的数值。对于连续缺失的情况,我们逐步迭代地预测序列的缺失部分,先预测第一个缺失值,然后将其用于下一个缺失值的预测,以此类推。
3.1.2. P-Bshade模型空间维度补值
P-Bshade方法是在空间维度进行的插值方法,计算原理如下:
(3.1)
式中,
表示第i个空间周围采样数据的观测值;
表示第i个空间周围采样数据对缺失数据的空间贡献权重;
可以用下式计算求得:
(3.2)
式中,方程中间的矩阵为待求矩阵;
为拉格朗日系数。方程左边的矩阵中
为第i个空间附近采样点的时间序列与第
个空间附近采样点的时间序列的协方差,
为第i个空间附近采样点的时间序列与缺失数据点的时间序列的期望比。方程右边的矩阵中
为第i个空间周围采样点的时间序列与缺失数据点的时间序列的协方差,并满足
,
。
首先选取每天的MAIAC AOD缺失数据附近n个相关性最大的空间采样数据进行插值计算,采用相关系数R来说明相关性的强弱,对于一个AOD缺失序列,计算其附近空间点的AOD数据序列和缺失AOD数据点的数据序列的相关系数R,R越大则表示相关性越强,反之则越弱。之后在缺失点附近找到非空AOD序列且相关系数R最大的十组序列,构建拉格朗日方程组之后求解权重,最后得出缺失值。
3.2. Catboost模型构建
Catboost模型[32]是一种梯度提升术算法,是基于一种对称二叉树为基学习器的GBDT框架,专门用于处理类别型特征的机器学习问题,还解决了模型训练和预测过程中的梯度偏移和预测偏移问题,从而达到减少过拟合风险的目的,进而提高模型算法的准确性和泛化能力。该模型对PM2.5浓度显示出较好的估算能力[33][34]。在进行类别型特征处理的时候,常用的方法是用整个数据集的均值代替类别特征的值,并进行平滑处理,这样的方法对于单个样本的估计量是有偏的,Catboost算法通过使用除去该样本的数据子集而不是整个数据集来进行估计,公式如下:
(3.3)
式中
为第k个样本的第i个样本特征,
为第j个样本的标签特征值,
为第k个样本前的第j个样本的第i个类别特征,
为随机序列中在第k个样本前的数据集,a通常为大于0的参数,p为先验项。
本文使用Catboost模型来建立MAIAC AOD-PM2.5之间的关系,并加入温度(TEMP)、风速(WS)、边界层高度(BLH)和相对湿度(RH)作为相关预测因子,与AOD一起建立模型对PM2.5浓度进行估算,模型建立流程如图4:
Figure4.Flow chart of PM2.5estimation by Catboost
图4.Catboost估算PM2.5浓度流程图
3.3. 评价指标
AOD补值和PM2.5的浓度估算都采用平均绝对误差(MAE)、均方根误差(RMSE)和皮尔逊相关系数(R)作为评价指标来进行精度验证,AOD补值精度验证使用地面站点AERONET AOD数据与补值结果进行验证,PM2.5的估算精度验证使用PM2.5站点与估算结果进行验证。各评价指标计算见下式:
(3.4)
式中
和
分别是模型计算和站点监测的AOD值(PM2.5值),n是样本的个数,
为模型计算值和站点监测AOD值(PM2.5浓度)的协方差,
为模型计算结果的标准差,
为站点检测值的标准差。
4. 结果分析
4.1. 精度指标评价
据气象学标准定义的实际季节来定义季节:春季为3、4和5月,夏季为6、7和8月,秋季为9、10和11月,冬季为12月、次年的1月和2月。本文对2020年3月~2021年2月的MAIAC AOD进行时空补值并对PM2.5浓度进行估算。图5为AOD时空补值数据与地面站点数据的线性拟合结果对比图,图6为PM2.5浓度估算结果与地面站点数据的线性拟合结果图。
Figure5.Fitting diagram of AOD top-up results
图5.AOD补值结果拟合图
本文所用的Prophet-LSTM + P-Bshade时空补值模型相对于常用的经典克里金模型和时空克里金模型有较高的拟合精度,R、MASE和MAE分别为0.891、0.275和0.183,经典克里金补值模型的R、MASE和MAE分别为0.751、0.518和0.337空克里金补值模型的R、MASE和MAE分别为0.834、0.521和0.341。从各模型补值的结果和站点监测值的收敛来看,本文提出的补值方法收敛性最好。经典克里金模型是基于数据空间相关性进行的补值方法,未考虑数据的时间相关性,因而模型补值结果较差,时空克里金模型是经典克里金模型在空间维度上的延伸,在补值时需要多个时间点的空间截面数据进行时空连续性插值,在进行大范围的补值研究时,计算量十分庞大,且AOD的补值中存在低估现象。本文提出的Prophet-LSTM + P-Bshade时空补值模型基于AOD数据的时空相关性对其进行了补值,且相对于对比模型来说有最好的补值效果。
将ERA5气象数据和经过Prophet-LSTM + P-Bshade补值后的AOD数据作为估算变量输入到Catboost模型中,即可得到PM2.5浓度的空间分布数据,图6为模型估算结果与PM2.5站点拟合结果图。
本文选择LightGBM等八种常用机器学习、随机森林方法与Catboost模型进行对比,通过图6可知,Catboost模型的拟合精度最高,R、MASE和MAE分别为0.93、15.89 μg∙m−3和10.54 μg∙m−3,这些指数均优于其他八个对比模型。
其中较为传统的Bagging模型和KNN模型相对来说拟合精度较差,LightGBM、XGBoost和Catboost模型是GBDT的三大主流模型,都是在GBDT算法框架下进行了改进。这三种模型在PM2.5浓度估算中都显示有较高的拟合精度。其中Catboost模型主要使用了Ordered Target Statistics方法将类别特征转化为数值特征、基于贪心策略的特征组合方法、使用Orcdered boostng避免梯度偏移问题和使用对称二叉树作为基模型,其拟合效果最好。
Figure 6.PM2.5concentration fitting results
图6.PM2.5浓度拟合结果图
4.2. PM2.5时空分布分析
将补值后的AOD数据与ERA5气象数据输入Catboost模型中可以得到2020年3月至2021年2月的每日PM2.5浓度估算空间分布图,再按4.1节方法划分季节,按季节取均值可以得到PM2.5季均空间分布图,如图7。
从季节分布上来看,2020年四个季节PM2.5浓度分布特征较为明显,整体呈现冬季(51.37 μg∙m−3) > 春季(37.40 μg∙m−3) > 秋季(24.86 μg∙m−3) > 夏季(21.95 μg∙m−3)的季节分布特点。从空间分布上来看,我国陆地区域PM2.5浓度呈现东高西低的特点。东部地区经济发达,我国三大经济圈(环渤海经济圈、长江三角洲经济圈和珠江三角洲经济圈)主要城市基本都在中国东部,高速城市化、工业化带来经济发展的同时,城市生态环境空间被大量蚕食、大量的流动人口朝着经济发达的地方聚集,污染排放的强度和密度剧增,使得这些经济发达的地方成为大气环境污染的重灾区[35]。西部地区经济相对落后,PM2.5浓度整体浓度偏低,但是塔里木盆地区域PM2.5浓度局部较高,中国最大的沙漠——塔克拉玛干沙漠位于塔里木盆地附近,长期的沙尘暴以及风蚀现象,会将高浓度的气溶胶颗粒物带到不同的地方从而影响当地气候,包括整个西北地区和华北地区等,不同程度上都会受到沙尘的影响,沙漠活动等自然因素影响是导致塔里木盆地区域的PM2.5浓度较高的主要原因[36]。
Figure 7.Spatial distribution of seasonal average PM2.5concentrations in 2020
图7.2020年PM2.5浓度季均空间分布图
5. 结论
本文使用建立的Prophet-LSTM + P-Bshade时空补值模型对MAIAC AOD数据进行时空补值,并将其与ERA5气象再分析数据通过Catboost模型进行PM2.5浓度估算,为我国大气环境污染治理提供了数据支撑。对2020年3月~2021年2月的全国陆地区域进行PM2.5浓度估算,分析得到以下结论。
(1) Prophet-LSTM + P-Bshade时空补值模型精度明显优于传统补值方法,R、MASE和MAE分别为0.891、0.275和0.183。
(2) Catboost模型在PM2.5浓度估算中相较于其他八个常用模型显示更高的估算精度,R、MASE和MAE分别为0.93、15.89 μg∙m−3和10.54 μg∙m−3。
(3) 中国陆地区域2020年的PM2.5浓度在季节尺度分布上明显,整体呈现冬季 > 春季 > 秋季 > 夏季的季节分布特点。在空间分布上,PM2.5浓度整体呈现东部地区较高,塔里木盆地区域局部较高的特点。
(4) 提出的Prophet-LSTM + P-Bshade时空补值模型在将AOD数据进行时间维和空间维补值结果线性融合时,本文顾及时空平稳性将两个维度的权重都取值为0.5,未来可以尝试其他的方法来设定时空维度权重,以此达到更好的补值效果。