基于AGRI 颗粒物浓度遥感反演及季节变化分析

2021-03-17鲍艳松许梦婕顾英杰南京信息工程大学气象灾害预报预警与评估协同创新中心气象环境卫星工程与应用联合实验室中国气象局气溶胶与云降水重点开放实验室江苏南京210044南京信息工程大学大气物理学院江苏南京210044

中国环境科学 2021年2期

董焱,许丹*,鲍艳松,许梦婕,顾英杰 (1.南京信息工程大学,气象灾害预报预警与评估协同创新中心,气象环境卫星工程与应用联合实验室,中国气象局气溶胶与云降水重点开放实验室,江苏南京 210044；2.南京信息工程大学大气物理学院,江苏南京 210044)

大气颗粒物浓度PM2.5是对人类健康影响较大的大气污染指标,因为其颗粒直径较小,粒子可以从肺部进入血管,同时有研究表明PM2.5微粒上附着的有害的物质及重金属可能进一步危害人类健康,目前PM2.5浓度的研究广泛受到关注[1-2].气溶胶光学厚度(AOD)是卫星和地面遥感大气气溶胶容易获得的光学参数,通过卫星还能得到区域气溶胶分布特征.该参数表示整层大气(从地面到大气层顶)气溶胶的削光量,是气溶胶消光系数在垂直方向上的积分.由于气溶胶本身对光路存在丁达尔效应,因此对大气辐射和能量的收支平衡起着重要的作用,从而导致气候变化的影响[3-4].因此大气气溶胶的定量是卫星遥感的重要参数之一[5-6].

卫星仪器探测的每一个像元都是反映整层大气柱的气溶胶光学特性,而大气边界层是人类活动的主要范围.通过卫星仪器(一般是可见光波段)估算近地面PM2.5浓度时,当高层气溶胶含量较高时,两者的相关性会受到影响.粒子的吸湿作用如华东地区工业排放导致的粉煤灰以及未燃烧的碳在大气中凝结水汽[7],导致AOD 受相对湿度影响较大.而颗粒物浓度结果不受相对湿度的影响.因此,AOD 与PM2.5在通常情况下存在非线性相关性[8].

利用卫星传感器在不同波段下观测气溶胶的光学特性,来估算观测区域的PM2.5[9].其观测结果具有高时效性、监测成本低、且监测范围广等优势[10-11].目前,关于此方面的工作,主要分为数学统计物理模型、人工智能,以及化学模式模拟等方法.其中数学统计的方法开展的最多,包括回归模型[12]、多元回归模型[13-14]、地理加权回归[15]等,在此基础上建立了气溶胶细模态光学厚度与大气颗粒物PM2.5质量浓度之间的回归关系[16];以及气溶胶反射率与大气颗粒物PM2.5质量浓度的线性回归关系等[17].目前对于人工智能/机器学习法开展的研究不多,基于神经网络方法对PM2.5质量浓度进行估算[18].非线性的机器学习的方法也进行了很多不同算法之间的尝试,但是和数学统计的方法一样都需要较多的数据进行支撑,对数据的质量的要求较高.而利用化学模式进行模拟,其依赖于模式本身.而物理机制建立模型的方法具有物理意义、通用型强、业务维护成本较低等特点.

由于我国国土面积辽阔、区域人口密度不协调、部分地区自然条件恶劣等一系列外在因素,想实现地面检测网成为了一项具有挑战性的任务[19].如今,华北地区形成了较为密集的观测网,但地面监测网需要协调较多仪器,单个仪器之间的校准定标都存在一定差异.对于PM2.5的浓度局地趋势延判具有较高难度.而通过卫星观测反演可以更好的解决上述问题.

利用FY4A/AGRI 仪器反演得到的AOD 数据,以及对应的地面相对湿度数据,估算近地面大气颗粒物PM2.5浓度.根据大气气溶胶消光的物理机制和垂直分布特征,运用大气辐射传输模型6S(Second Simulation of the Satellite Signal in the Solar Spectrum)能见度经验转换公式,计算得到气溶胶标高并对卫星探测的整层大气消光系数进行垂直订正,估算近地面PM2.5质量浓度数.

1 数据来源和研究区域

1.1 FY4A/AGRI AOD 数据

FY4A 是中国发射的第二代静止气象卫星,2016年12 月11 日发射升空,并在2017 年9 月25 日正式投入用户使用[20].其搭载的扫描辐射成像仪跻身于世界静止轨道成像仪最先进行列[21](Advanced Geostationary Radiation Imager,AGRI).AGRI 每15min 生成一副全圆盘影像观测,且一共拥有14 个通道,其中两个可见光通道(红光),蓝光以及近红外、短波红外、中波红外和热红外通道等[22].除可见光通道外, AGRI 的空间分辨率在4km,这有利于多通道数据的计算,为利用数据进行同化和反演提供了较大的便利[23].

本文基于团队自主研发的FY4A/AGRI 气溶胶光学厚度数据集进行分析研究.该AOD 数据集采用暗像元法,利用AGRI 前6 个通道,基于6S 模型建立查算表反演得到,详细原理和验证过程另有文章阐明,在此不再赘述.本文利用此数据集,对PM2.5浓度进行估算.

1.2 NPP/VIIRS AOD 数据

搭载在NPP 卫星上共有五个传感器,其中包括VIIRS(Visible infrared Imaging Radiometer)即可见光红外成像辐射仪.在2011 年10 月28 日启动,星下点空间分辨率为400m[24],VIIRS 相对于MODIS 具有更高的分辨率[25].

1.3 PM2.5 数据和相对湿度数据

图1 华北东部地区环境监测站点分布Fig.1 Map showing the location of the 180air quality monitoring stations

采用地面的相对湿度数据参与近地面颗粒物PM2.5的估算.同时,采用了地面颗粒物PM2.5数据来进行PM2.5估算的检验对比.利用美国METONE 公司所生产的手持式空气尘埃粒子计数器进行近地面PM2.5测量.其数据覆盖了京津冀和山东地区.数据来源于各个省市的环境监测中心180 个地面站点(图1)信息,包括小时平均相对湿度数据和对应的颗粒物数据,如图一所示.采用2018 年4 月～2019 年1 月不同季节高污染天气前后5d 的数据样本进行与地面数据的精度检验,并提取2018 年4 月18 日重污染天气个例经预处理后与地面172 个站点进行个例分析对比.

2 区域PM2.5 浓度反演原理

2.1 辐射传输模型

在原理上,由于卫星所得到的AOD 数据,每一个像元所对应的是卫星探测到的整层气溶胶光学厚度即大气柱上的总消光系数,但是近地面的颗粒物PM2.5浓度所指的是对流层以下甚至是干粒子的质量浓度,并且大气气溶胶的吸湿性对其消光特性具有相当大的影响[26].所以AOD和PM2.5两者的相关性受大气分布和大气状态等一系列因素的影响[27].本文通过对AOD-PM2.5建模的分析,将影响较大因子放入模型中参与计算,具体计算流程如下图2 所示.

图2 实验流程Fig.2 Experimental flow chart

2.2 垂直订正与湿度订正

2.2.1 垂直订正假设不同的大气层各边界之间相互平行,则AOD 是对整个大气柱的气溶胶消光系数积分,可以通过公式(1)进行表示:

式中:τa(λ)为AOD、Ka为气溶胶消光系数在波长为λ 高度为z 时的值.可以通过估算的方式得到Ka(λ,z)的函数表达式[28],具体如下:

式中:Ka,o(λ)是表示地表气溶胶的消光系数即z=0; HA为气溶胶标高,m;气溶胶标高可以看成大气边界层高度(Boundary Layer Height,简称BLH)[29].因此,式(1)、式(2)结合可以得到:

由式(3)中可以发现,当Ka,o(λ)、τa(λ)和HA如果有两个是已知的,则可以计算到第三个.由于本文基于业务化系统进行考虑,对于大气边界层高度数据的获取利用再分析资料的时效性都无法实现实时的业务化系统监测,因此对于如何利用已知的AOD 数据来估算气溶胶标高变成了难题.通过在大气辐射传输模型6S 中有AOD 与能见度的经验转换模型:

式中:V 为近地面能见度m[30].增加在水平方向上的瑞利散射和臭氧吸收所带来的影响[31-32],将0.02 的对比阈值更改到了0.0416,得到公式(5),具体如下:

根据式(5)估算出气溶胶标高HA.再将HA带入式(4)中,得到垂直订正后的消光系数Ka,o(λ).

2.2.2 湿度订正由于PM2.5指的是干粒子的质量浓度,需要对AOD 数据进行湿度订正[33].因此,吸湿性增长因子可以表示为[34]:

式中:RH 为地表相对湿度; f(RH)为吸湿性增长因子;g 为经验拟合系数.华北东部地区的经验拟合系数g一般取值为0.38[35].

因此,湿度订正后的消光系数Ka,o(λ)可以表示为[36]:

因此,将上述式(8)与式(5)结合,得到近地面干气溶胶消光系数Ka,0,Dry(λ),公式如下:

2.2.3 气溶胶类型参数化根据粒子的米散射理论,大气消光系数Ka定义为[37]:

式中:Qext(m, r ,λ )为大气气溶胶消光效率,m、r 分别表示单个粒子的复折射指数、粒子尺度,n(r)表示数密度谱分布.粒子尺度的取值范围为(0,x/2)µm.

颗粒物PMx在大气中的质量浓度可以表示为:

式中:ρ为对应尺度的气溶胶质量密度,单位为g/m3,北京附近地区,ρ取1.5g/m3[38].Hansen等[39]研究发现,消光效率Qext和有效半径reff可以表示为:

表1 MODIS 在不同气溶胶类型的消光特性表Table 1 Extinction characteristics of MODIS in different aerosol types

因此,将式(9)、(10)、(11)、(12)整理:

式中:Ka为湿度订正和垂直订正之后的大气气溶胶消光系数,X 取2.5.在MODIS 的二级数据MOD04产品中,将不同气溶胶类型(Aerosol Type 即Aerosol_Type_Land)进行了对消光效率Qext和有效半径reff的分类,如表1 所示,选取其中大陆型的消光效率与有效半径作为本次实验的参数.

2.3 星-地数据时空匹配

将不同卫星的AOD 数据进行预处理,包含对个别数据文件的边角不整齐所进行的统一切割;对数据的缺省进行相应的质量控制;对数据在程序中占用空间太大所进行的数据格式转化以及数据保留小数点后几位的选择等.

将每一个卫星像素点匹配上与其像素点距离最近的地面信息数据,如式(14)所示:

式中: Dmin表示为两点的最小距离;lonsat表示卫星某个像素点的经度数值;latsat表示卫星某个像素点的纬度数值.同样,longro表示地面测站某个像素点的经度数值,latgro表示地面测站某个像素点的纬度数值.

在时间匹配上,根据不同卫星数据的过境时间信息寻找就近的相对湿度数据来进行估算以及对应地面颗粒物浓度数据进行结果的检验.数据总结发现:MODIS 与VIIRS 都以16d 为一个飞行周期,且每天飞过实验地区为10:00～11:00.因为本文所使用的地面信息站点的数据为小时平均,与卫星数据进行时间匹配时间窗小于±1h.

3 结果与讨论

3.1 AGRI 个例结果展示

图3 为FY-4A/AGRI PM2.5 区域分布图,取自2019 年10 月28 日11:00 数据个例.在山东省的西南地区以及东部的半岛沿海地区存在PM2.5浓度较大的区域.同时,北京和天津的PM2.5浓度同样较高.

图3 FY-4A/AGRI PM2.5 区域分布Fig.3 FY-4A/AGRI PM2.5 regional distribution map

图4 是基于AGRI 计算出的PM2.5结果与地面站点匹配的数据折线图,取自2018年4月18日11:00数据.利用AGRI 得到的PM2.5与地面观测值整体趋势一致.在110号站点附件AGRI都存在高值,其原因为卷云的干扰或是云边界像素,导致云识别没有将此像素剔除,由于卷云反射率与地面不同,因此卫星计算值偏高.

图4 AGRI 反演结果与地面结果对比Fig.4 Comparison of AGRI inversion results and ground results

3.2 不同季节各卫星的对比

由于AOD 具有光学特性,其散射和吸收也会因不同季节太阳直射点的变化而受到影响.对颗粒物浓度夏季＜秋季/春季＜冬季,春季和秋季在PM2.5污染程度大致相同[40].本文随机选取了夏季2019 年9 月18 日、秋季2018 年10 月13 日、春季2018 年4 月18 日和冬季2019 年1 月10 日进行对比,如图5 所示.四个季节相比地面站点颗粒物PM2.5的高值区与低值区两个卫星都有较好的表现,其次相比于地面颗粒物站点数值,卫星数据估算的结果使其PM2.5的整体趋势显得更加直观.在秋季(Ⅰ),无论是地面站点还是利用卫星资料估算的颗粒物PM2.5都在天津地区存在一个高值区,在河北和站东北部沿海地区具有较高的PM2.5质量浓度.其质量浓度数趋势沿着海湾向四周递减.利用卫星估算的PM2.5结果与地面站点探测数据结果吻合.在春季(Ⅱ),地面站点PM2.5数据与卫星估算PM2.5结果在北京地区附近都存在高值.两个卫星的结果具有较高的一致性,与地面站点探测数据相统一.在冬季(Ⅲ)MODIS 与VIIRS 对比差异相对较大.因为算法不同,因此导致舍去的像元不同.在MODIS 数据结果中,河北大部分地区包括北京和天津地区都存在缺省值,主要是这些地区没有达到MODIS 官方暗像元的级别要求从而在反演中将其剔除.MODIS 数据删去了大量的像素点,但是其保留下来的像素点与VIIRS 数据依然存在较高的一致性.相比于地面站点河北南部以及山东北部的颗粒物高值地带,这两颗卫星都有较好的表现.在夏季(Ⅳ),山东地区从沿海到内陆PM2.5质量浓度都有逐步变大的趋势,两颗卫星结果都与地面站点大致相符,趋势一致,但是VIIRS 不如MODIS 结果明显.其原因主要是基于仪器标定的差异以及VIIRS AOD 网格数据精度较低所导致.在地面站点图中,北京地区具有一个PM2.5的低值区.这两个卫星相对于地面结果同样具有相同的趋势.可以初步判断,此方法可以用在多种卫星估算颗粒物PM2.5质量浓度上,受季节的变化影响较小,模型具有较高的稳定性.

3.3 反演结果精度检验

根据上述时空匹配的方法,对AGRI AOD 数据进行了多个个例结果的运行,将地面颗粒物PM2.5站点数据与AGRI、MODIS 和VIIRS PM2.5数据进行统计和评价.其中包括对散点对比、样本个数(N)、相关系数(R)、均方根误差(RMSE)、平均偏差(ME)以及相对误差(MRE)等几个方面的精度评价.在时间匹配上,采用的是地面站点PM2.5质量浓度小时平均数据;空间匹配上,地面PM2.5站点数据与卫星像素点值的距离范围小于0.01 经纬度.若地面站点或卫星像素点存在缺省则此匹配点不参与验证.在质量控制中,剔除了卫星数据中的缺测值(-9999)以及地面站点的漏报值(-999);同时,若卫星估算值和地面站点值差的绝对值大于其三倍的均方根误差(RMSE)进行剔除,以保证其它外在因素的干扰.

将地面测站结果与卫星数据结果进行上述预处理后,结果如图6 和表2 所示.选取春季2018 年4 月14～18 日的10:00 数据;夏季为2018 年7 月13～18 日10:00;秋季和冬季分别选取2018 年10月13～17 日10:00 和2019 年1 月10～14 日10:00数据,分别对不同卫星仪器的数据进行了精度检验分析.

图5 不同季节的多颗卫星反演结果与地面结果比较Fig.5 Comparison of satellite retrieval results and ground results in different seasons

图6 四个季节不同AOD 数据估算近地面PM2.5 结果与地面测站PM2.5 结果散点对比Fig.6 Scatter comparison of PM2.5 results from different satellite instruments in four seasons and PM2.5 results from ground stations

如图7 所示,采用中值滤波的方法对平均值进行平滑处理.其中,夏季(b)PM2.5平均浓度在京津冀地区达到了60μg/m3.其他3 个季节的平均浓度最大值约在80～100μg/m3之间.因此,夏季平均浓度相较于其他三个季节较低.相反,秋冬两季,北京和天津地区的PM2.5平均浓度值为90～100μg/m3,比其他地区高,冬季更为明显.并且除极值区外,冬季的低值区也有近50μg/m3,相较于其他3 个季节的低值区也相对偏高.进一步说明,夏季的颗粒物浓度普遍低,而冬季的颗粒物浓度普遍较高.

在以下4 个季节中AGRI 与地面数据相比都有相对较好的结果,相关系数R 都在0.8 以上,进一步反映出AGRI 的稳定性.AGRI 在四个季节相对误差(MRE)都小于20%,其中夏季达到17.24%.由表2 可知,AGRI 在冬季和春季的精度检验结果相对较好,夏季的精度普遍偏低[41].原因主要有两点,其一是因为夏季的颗粒物浓度相对于其他三个季节浓度相对较低,卫星仪器的灵敏程度对测量浓度较低的PM2.5偏差较大.其二是因为夏季的太阳高度角的变化使太阳直射点北移,导致暗像元算法受到季节的影响从而影响了暗像元识别的等级.但上述两点仍需要进一步验证.AGRI 的结果明显优于VIIRS 和MODIS 的结果.AGRI 在4 个季节中的相对误差都小于VIIRS 和MODIS,AGRI 与地面的数据更加接近,AGRI 的稳定性更好.在颗粒物浓度较大的春秋两季,AGRI 的结果总体看也略优.夏季颗粒物浓度低于其他3 个季节.因此夏季的散点在低值区普遍较多,在相关性上也略微偏低.但其结果仍不亚于MODIS 和VIIRS 的夏季精度检验结果.