基于Himawari-8卫星的气溶胶光学厚度反演

陈 瑞 游浩妍 万 翔

(1. 自然资源部第一地理信息制图院, 陕西 西安 710054;2. 自然资源部陕西基础地理信息中心, 陕西 西安 710054)

0 引言

气溶胶是大气系统重要构成,对气候变化、人类健康有重要影响,在全球和区域气候变化中扮演十分重要的角色[1-2]。气溶胶光学厚度是描述光的衰减作用的无量纲参数,能表征大气浑浊度,是大气气溶胶基本光学特性之一[3]。卫星遥感可以开展气溶胶光学特性的宏观监测,为获取气溶胶时空分布特征、变化趋势及溯源信息提供有效监测手段。

高时间分辨率是地球静止卫星的优势,能实现分钟级的对地观测,对气溶胶的日变化情况研究具有巨大的应用潜力[4-6]。葵花8号 (Himawari-8,H8)是地球静止卫星的代表,于2015年7月开始向全球实时广播数据[7]。目前基于H8数据开展的气溶胶监测已有相关研究:文献[8]利用暗像元法开展了H8气溶胶反演,并对结果进行验证,精度可行;文献[9]利用新研究算法对H8卫星气溶胶光学厚度进行反演,并通过全球气溶胶地基观测网(aerosol robotic network，AERONET)站点进行误差分析,说明通过H8反演的结果在雾霾传输、成因分析等方面有着广泛应用价值;文献[10]构建了基于最优估计技术的H8卫星气溶胶反演算法,并验证了可行性。本文参考中分辨率成像光谱仪(moderate-resolution imaging spectroradiometer，MODIS)暗目标法,构建H8的气溶胶光学厚度(aerosol optical depth，AOD)反演算法,并将结果与AERONET、MODIS AOD进行对比,对反演结果进行验证分析。

1 资料与方法

1.1 资料

1.1.1Himawari-8数据

H8数据来源于日本气象厅,于2015年7月正式投入使用,用户通过网站注册后可以免费获取数据(https://www.eorc.jaxa.jp/ptree/)。搭载的高级成像仪(Advanced Himawari Imager,AHI)传感器,包括16个通道,覆盖可见光、近红外、红外波段,星下点空间分辨率第3波段为0.5 km,1、2、4波段为1 km,5～16波段为2 km[11]。先进葵花成像仪(advanced himawari imager，AHI)全圆盘图完成周期是10 min,特定地区编码扫描可达2.5 min。本文利用H8的1～6波段进行AOD反演,其中云检测利用第1(0.46 μm)、2(0.51 μm)、4(0.86 μm)、5(1.6 μm)波段进行计算;AOD利用1(0.46 μm)、3(0.64 μm)、6(2.3 μm)波段进行反演。

1.1.2MOD04数据

MOD04数据是美国国家航空航天局(National Aeronautics and Space Administration，NASA)在国际上公开发布的包含各类气溶胶特性参数的AOD产品,数据在官网(https://ladsweb.modaps.eosdis.nasa.gov/search/)免费获得。产品经过大量地面实测站点验证,具有较高的精度,空间分辨率3 km,一天过境2次,过境时间为地方时10 ：30和13 ：30左右,难以满足小区域气溶胶的日变化监测分析,常作为对比产品来检验反演结果的精度。本文利用MOD04产品中暗像元法反演的气溶胶产品,通过时空匹配作为H8 AOD反演结果的对比验证数据,空间分辨率为3 km。

1.1.3AERONET数据

AERONET可以采集更为精确的气溶胶特性数据,是NASA联合其他科研机构等设立的自动观测网,提供了全球范围内多个站点的气溶胶光学厚度数据[12],利用CE-318型太阳辐射计以及统一的算法获取AOD信息,数据可从官网免费获得(https://aeronet.gsfc.nasa.gov/)。本文利用北京区域AERONET长期观测站点:Beijing、Beijing_CAMS、XiangHe,采用经人工筛选质量级别高的Level2.0产品作为验证数据,利用440～675nm波长处的AOD,插值得到550nm波段的AOD。

1.2 反演原理

卫星获取的表观反射率变化能够体现气溶胶消光作用,从而进行AOD反演。根据大气辐射传输方程,在地表郎伯体和大气水平均一的前提下,表观反射率ρTOA可以表示为式(1)。

(1)

式中,μS,μV,φ分别为太阳天顶角余弦、观测天顶角余弦、相对方位角;ρ0是大气中分子和气溶胶散射的反射率;ρS是地表反射率；T(μS)是入射方向的大气透过率；T(μV)是观测方向的大气透过率;S为大气下界面的半球反射率。对于单次散射,ρ0与气溶胶光学厚度τ、气溶胶散射相函数Pa、单次散射反照率ω0以及分子散射造成的路径辐射ρm有关[13],公式表示为式(2)。

(2)

由式(2)可见,当ρS很小时,ρTOA主要取决于大气贡献项;当ρS很大时,地表成为主要贡献项。卫星数据可以提供角度信息及表观反射率信息,气溶胶散射相函数Pa、单次散射反照率ω0以及路径辐射ρm、大气透过率T(μS)T(μV)、半球反射率S等参数我们通过6S模型计算获取。因此,求得地表反射率后,卫星观测的表观反射率除掉地表反射率部分即得到AOD值。

MODIS暗目标法是在清洁大气情况下,暗目标地物(植被、水体、黑土等)在近红外通道地表反射率很小,且近红外波段与红蓝波段地表反射率具有较高的相关性,因此利用近红外波段估算红蓝波段的地表反射率,通过6S模型建立查找表,实现AOD反演[14]。

1.3 反演方法

利用交互式编程语言实现AOD的反演,工艺如下:将红外波段中数值小于0.25的像元选择为暗像元,并得到红蓝波段相同位置的暗像元;根据MODIS C6版本中暗目标法的地表反射率关系计算红蓝通道地表反射率;读取每个像元角度信息,并寻找查找表中对应的角度值上下限,对上下限角度进行插值,得到大气透过率T(μS)T(μV)、半球反射率S、路径辐射ρm等大气参数;通过3个大气参数、红、蓝波段地表反射率结合式(1),得出红蓝波段的理论表观反射率。最接近实际表观反射率的查找表AOD为反演结果,最后对红蓝波段的反演结果进行平均。

1.3.1云掩膜

统计H8数据云覆盖区域的表观反射率情况,设置阈值进行去云处理,当3×3窗口的表观反射率ρ0.47>0.25,或者表观反射率ρ0.51标准差>0.006,则窗口内的9个像元均认为是云。

式中,σ0.51为标准差;ρ0.51为H8数据波长0.51 μm波段;ρ0.47为H8数据波长0.47 μm波段;μ0.51为窗口内H8波长0.51 μm波段的平均值。

1.3.2地表波段关系

MODIS C6版本中暗目标法的地表反射率关系为式(5)和式(6)。

式中

H8波段设置同MODIS相似,利用MODIS波段间地表反射率关系进行H8卫星AOD反演。由于H8 AHI传感器没有1.24 μm通道,统计0.86 μm波段与1.24 μm波段关系,以替代1.24 μm计算的NDVI[15]。最终适用于H8的波段地表反射率关系见式(8)和式(9)。

(9)

式中,

(10)

1.3.3查找表建立

利用6S模型建立北京区域5月1日查找表,输入参数包括:平均海拔0.47 km;太阳天顶角为24°～68°,间隔1°;观测天顶角为51°～55°,间隔1°;相对方位角为0°～180°,间隔12°;550 nmAOD为0～2.5,间隔0.1;选择大陆型气溶胶;根据上述参数以及6S模型模拟的大气参数建立査找表,如表1所示。

表1 查找表示例

1.3.4查找表插值

查找表包含的角度值不完全等于像元读取的角度,故对在查找表中查找到对应角度的上下限进行插值得到前面3个参数的插值结果,以此作为参数计算表观反射率。分析查找表发现:相同AOD的半球反照率相同,太阳天顶角(Solar zenith angle, SOZ)和卫星天顶角(satellite zenith angle，SAZ)同大气总通过率成反比,与大气中气溶胶反射率成正比;相对方位角与大气中气溶胶反射率成反比,与大气总透过率无关。插值方法见式(13)～式(16)。

式中,C为插值系数;Z为太阳天顶角;ZA为卫星天顶角,φ为此像元读取的角度数据;S为半球反照率;T为大气总透过率;P为大气中气溶胶反射率,通过查找表1得出,Z1、ZA1为查找出的角度数据上限,S1、T1、P1是上限对应的半球反照率、大气总透过率、大气中气溶胶反射率,Z2、ZA2、S2、T2、P2为查找出的角度数据下限。

2 结果与分析

本文以北京为研究区域,选择2020年5月1日UTC时间00 ：00—08 ：00逐小时数据进行反演,当日植被情况较好、云量较少、有地面验证点,可保证足够的反演结果。当日空气质量指数(air quality index，AQI)为204,属于严重污染天气。

2.1 与MODIS结果对比

H8卫星反演AOD结果为2 km分辨率,重采样为3 km分辨率后与MODIS暗像元的AOD产品结果进行对比分析,从空间分布图可见:H8和MODIS的AOD结果空间分布一致,但存在细微差别：①来自云掩膜方法不同;②分辨率不同重采样后进行对比。

从回归分析结果可见:H8反演AOD结果与MODIS AOD结果相关系数为0.854,相关性较好,表明算法具有可行性,如图1、图2所示。

(a)UTC 03 ：50 H8反演

AOD结果 (b)UTC 03 ：45 MODIS AOD

(c)UTC 05 ：30 H8反演

AOD结果 (d)UTC 05 ：25 MODIS AOD

图1 2020年5月1日H8反演AOD结果与MODIS AOD结果分布对比

图2 2020年5月1日UTC 03 ：45及05 ：25 H8 AOD反演结果与MODIS AOD结果对比

2.2 与AERONET站点实测结果对比

图3 2020年5月1日H8 AOD结果与AERONET站点对比

H8卫星反演结果为550 nm处的AOD,而AERONET没有直接观测的数据,根据其440 nm、500 nm和670 nm处的观测结果,利用二次多项式拟合方法,计算得到550 nm处AOD直接观测结果。结果对比时,时间上将H8观测时间前后5 min的地测结果平均值作为该时刻观测值;空间上,以站点为中心,半径10 km区域内H8 AOD像元均值为该站点处AOD反演结果。验证结果如图3所示,可见AOD反演结果和AERONET站点结果相关系数为0.809,相关性较高,与MODIS的对比分析情况相似,具有低值区AOD略高、高值区略低的趋势。

3 结束语

本文构建了H8的AOD反演算法,并以北京地区2020年5月1日为对象反演了AOD结果,通过与AERONET 3个站点Beijing-RADI、Beijing-CAMS、Xianghe的观测数据(MODIS AOD产品)进行对比,结果具有较好的相关性和验证精度,可见利用暗目标法对H8卫星反演AOD具有较大的应用潜力。相比于MODIS AOD产品,该算法时间分辨率高,能得到10 min级、空间分辨率2 km的AOD日变化结果,对于动态监测气溶胶情况具有较好的潜力。存在的问题及下一步工作包括:①目前查找表构建时采用大陆型气溶胶,下一步通过研究区域气溶胶实际构成情况,采用自定义气溶胶模式,提升AOD反演结果精度;②因地面观测站点有限,将MODIS地表反射率关系转化为H8卫星,后续可进一步优化波段地表反射率关系,实现地表反射率的精确估计。