基于卫星多角度观测的气溶胶遥感反演算法研究进展

吴孔逸,侯伟真,史正,许华,温亚南

(1自然资源部信息中心,北京 100812;2中国科学院空天信息创新研究院,国家环境保护卫星遥感重点实验室,北京 100101)

0 引言

气溶胶是大气中重要的组成部分,对全球辐射平衡、气候变化和人类健康等起着直接或间接的作用[1-6]。气溶胶光学厚度(Aerosol optical depth,AOD)是描述气溶胶光学特性的重要参数,也是地基和卫星遥感需要反演的首要光学参数。地基观测可以获得某一位置高时间分辨率的气溶胶特性,却无法对较大的范围进行监测[7,8]。相比于地基观测,星载监测气溶胶具有范围广、成本低、不受地理条件约束等特点[9,10]。针对气溶胶卫星遥感反演,由于要用少量观测来估计复杂的地气系统,本质上属于一个“病态”(Ill-posed)反演问题[11],而对应地表反射率的准确估计和地气解耦是其中的关键[12]。

搭载在美国EOS(Earth observation system)系列卫星上的中分辨率成像光谱仪(Moderate resolution imaging spectroradiometer,MODIS)是国际上具有代表性的单角度多光谱传感器,其气溶胶反演方法主要包括基于地表反射率波段传递的暗目标(Dark target,DT)算法[13-15]和基于全球地表反射率库的深蓝(Deep blue,DB)算法[16,17]。搭载在风云三号(FY-3)卫星上的中分辨率光谱成像仪(Medium-resolution spectral imager,MERSI),是我国发展的与MODIS类似的传感器;杨磊库等[18]通过借鉴MODIS的DT算法思想,针对升级版的MERSI-II传感器,进一步构建了全球适用的陆地气溶胶反演算法。

随着新型传感器的不断发展,卫星遥感对地观测模式也从最初比较典型的单角度观测,发展到现在流行的双角度和多角度观测。对于城市、沙漠等亮地表情况,气溶胶的标量贡献相对于地表反射贡献一般要小很多,属于弱信息,如果缺少地表先验知识很难进行有效的地气解耦。由于大气和地表对大气层顶卫星观测贡献的比率随不同观测角度而有所不同,多角度气溶胶遥感针对这一特性将二者进行分离[19-21]。此外,最新发展的偏振探测模式通过与多光谱和多角度的观测方式相结合,使得对气溶胶的观测信息得以丰富,可以获取更多气溶胶微物理和光学参数信息[22-24],有助于气溶胶特性的提取,为大气气溶胶定量遥感反演提供了新的途径。

双角度是多角度遥感的一种特殊情况,欧洲空间局(European Space Agency,ESA)研制的沿轨迹扫描辐射计(Along-track scanning radiometer,ATSR)系列传感器是经典的双角度观测传感器,对应的ATSR气溶胶反演方法是国际上双角度遥感气溶胶中最具代表性的经典方法[25-28],其核心思想也同样适用于两个传感器在几乎同时对同一地区进行观测的情况[29,30]。AASTR(Advanced along-track scanning radiometer)的后续传感器(Sea and land surface temperature radiometer,SLST,也称SLSTR)已和OLCI(Ocean and land color instrument)一同搭载在Sentineil-3卫星上,可为OLCI传感器的多光谱观测数据进行同步大气校正。

目前,我国也在发展自己的星载大气同步校正仪,如搭载在高分多模卫星上的SMAC(Synchronization monitoring atmospheric corrector)和环境减灾二号A/B卫星上的PSAC(Polarized scanning atmospheric corrector)等[31-34],通过多光谱和偏振观测模式相结合来同步探测和反演关键大气参数,为搭载在同一卫星平台上的主载荷传感器提供大气同步校正,以提升卫星遥感图像大气校正的质量。ATSR系列传感器双角度观测的通道设计及其双角度气溶胶遥感的反演方法,可为我国星载大气同步校正仪的研制和应用提供很好的借鉴。本文将针对国际上双角度、多角度及偏振卫星遥感反演算法的最新研究进展进行总结和综述,并为我国新发展的多角度观测卫星系列载荷气溶胶反演算法的研究提供有力的支撑。

1 ATSR双角度气溶胶反演

1.1 ATSR系列传感器

迄今为止,搭载ATSR系列传感器的卫星已经发射了5颗。其中,ATSR-1搭载在ERS-1卫星上,于1991年7月发射;ATSR-2搭载在ERS-2卫星上,于1995年4月发射;AATSR搭载在ENVISAT(Environment satellite)卫星上,于2002年3月发射[35]。2012年4月AATSR传感器停止工作,AATSR的后续传感器为搭载在哨兵3A/3B(Sentinel-3A/3B)上的SLSTR,分别于2016年2月16日和2019年8月25日发射升空。其中,ATSR-1传感器的波段设计主要用于红外探测获取海面温度,从ATSR-2开始设置了可见光和近红外波段的通道。ATSR-2和AATSR,统称(A)ATSR,通过采用底向(0°)和前向(55°)两个角度在短时间联合进行对地观测,进而获取对同一区域两个角度的多光谱观测信息,对应的7个波段分别为0.56、0.66、0.87、1.61、3.70、10.80、12.00 μm,其中前4个波段可用于大气气溶胶的研究[36]。

SLSTR继承了AATSR的双角度观测设计,但与AATSR不同的是它的双角度分别为底向和后向,并在AASTR原有波段的基础上添加了新的波段,与(A)ATSR相比,SLSTR增加了1.38 μm和2.25 μm两个红外通道,以及用于火点探测的3.74 μm和10.95 μm通道,共11个波段,可以更好地用于大气气溶胶的研究[37]。SLSTR底向的扫描幅宽为1675 km,前向的扫描幅宽为750 km。搭载在同一卫星平台上的OLCI传感器继承了MERIS(Medium resolution imaging spectrometer)传感器的设计,并将光谱通道在0.40～0.90 μm的范围内从15个增加到了21个,5个相机呈扇形分布进行穿轨(Cross-track)观测,幅宽约为1269 km[38,39]。在可见光和近红外波段,SLSTR的像元分辨率为500 m,双角度观测的幅宽覆盖OLCI,可以对其进行同步大气校正[35,40]。

1.2 ADV算法

ATSR-DV(ATSR dual view)算法简称ADV算法,其核心是通过大气层顶观测的表观反射率与由一个辐射传输模型建立的查找表(Look-up table,LUT)进行结果对比,从无云的像元中反演得到大气气溶胶的光学特性[41]。该算法同时利用双角度和多波段的信息进行气溶胶反演,而不需要事先确定地表反射率的信息。

ADV气溶胶反演算法基于如下两个假设:

1)不同类型气溶胶混合体在大气层顶的表观反射率可由其中每种气溶胶的表观反射率加权平均得到[42]。

2)在一定的光学厚度范围内,给定波段λ对应程辐射的表观反射率可近似表示为一个与气溶胶光学厚度τ有关的线性形式[43,44],即

式中:ρa(λ)为程辐射的表观反射率,ρ0(λ)为无气溶胶大气τ=0时程辐射的表观反射率,c(λ)为斜率。

对于一个朗伯(Lambert)地表,卫星传感器在大气层顶(Top of the atmosphere,TOA)接收到的表观反射率可以表示为

式中:ρTOA(λ)为表观反射率,ρa(λ)为程辐射的反射率,ρs(λ)为地表反射率,T(λ)为总透过率 (太阳-地表-卫星),S(λ)为半球反照率。

在双角度反射算法中,设两个方向的地表反射率之间的关系[45]为

式中:ρs,f(λ)为前向(Forward)的地表反射率;k为与波段无关的系数,主要与太阳和卫星的观测角度有关;ρs,n(λ)为底向(Nadir)的地表反射率。

通过结合两个角度的观测信息,联合式(2)和(3),并假设ρs(λ)S(λ)≪1可得到

这样可以将未知的地表反射率ρs(λ)消去,系数k可由1.6 μm波段的两个角度的表观反射率来确定(受大气影响较小)。将式(1)带入式(4)中,可推导得到气溶胶光学厚度(AOD)的迭代公式为

式中:ρTOA,f(1.6 μm)和ρTOA,n(1.6 μm)分别为卫星遥感前向和底向实际观测的表观反射率值。

利用式(5)和(6),从τ=0开始进行迭代,计算得到1.6 μm波段对应程辐射反射率的值,初始的k可直接由 ρTOA,f(1.6 μm)和 ρTOA,n(1.6 μm)取比值得到,对应参数的都可以利用查找表或辐射传输模式模拟计算得到,进而可以迭代得到对应波段的τ值。一般情况下,要使AOD达到1%的收敛精度,迭代几次即可满足对应的收敛条件[46,47]。

此外,k比值法为ADV迭代法的进一步改进,由式(4)可得到

这样每个波段(0.87 μm除外,植被在该波段有较强的反射特性)的k值可以通过模式模拟或直接查表计算得到(AOD都归结到同一波段上),而1.6 μm对应的初始k为

进而可引入反演所需的代价函数

ADV算法,包括迭代反演方法和k比值反演方法,都是基于同一像元地表的两个观测方向地表反射率的比值不变这个基础假设。Flowerdew和Haigh[45]已验证了这种假设关系是成立的,比值k在大部分区域主要与太阳和卫星的观测角度有关,与观测波段基本无关。其中,迭代方法只能利用单个波段的观测信息进行迭代反演;而k比值反演方法可以结合多个波段的观测信息进行联合反演。

1.3 SU算法

SU-AATSR算法以斯旺西大学(Swansea University)命名,也是一种针对(A)ATSR卫星气溶胶参数反演的常用方法,它主要是以North提出的适用于(A)ATSR卫星观测的地表反射率经验模型为基础[48]。

Grey等[49,50]基于(A)ATSR观测研究了地表二向反射,给出了与波段和观测几何相关的地表反射率经验模型

式中:Ω为观测几何(ATSR-2和AATSR的垂直和斜向观测);D(λ)为漫射辐照度的分量;v(Ω)只与观测几何有关;w(λ)只与波长有关;γ为多次散射系数,一般取0.3;g的表示形式为

在反演过程中,利用模型校正后的地表反射率ρs与模拟得到的ρmod值误差最小,对应的误差函数可表示为

这样就可以反演得到AOD和对应的气溶胶模型参数。

在AOD和大气气溶胶模型反演过程中,式(12)中需要的地表反射率ρs(λ,Ω)可通过对表观反射率进行大气校正得到[51],即

其中

式中:ρTOA(λ)、ρa(λ)、S(λ)和T(λ)的定义与式 (2)中的符号定义一致。

1.4 ORAC算法

ORAC(Oxford-RAL aerosol and cloud)算法是由牛津大学(University of Oxford)和卢瑟福·阿普尔顿实验室(Rutherford Appleton Laboratory,RAL)共同开发完成,并由早期的ATSR-2算法改进而来,可以反演气溶胶光学厚度、粒子有效半径和地表反射率等参数[52,53]。该算法首先基于最优估计(Optimal estimation,OE)理论建立耦合先验知识约束的价值函数,然后通过优化迭代方法获得价值函数的最优解,进而实现气溶胶和地表多参数的联合反演。

最优估计反演一般通过定义一个前向模型F来描述卫星观测数据对气溶胶和地表的依赖,其表达式为

式中:y为对应的观测向量,x为包含需要反演参数的状态向量,ϵ为误差向量。基于贝叶斯原理和高斯估计,y对于x的概率分布函数可以表示为

式中:Sϵ为测量误差协方差矩阵,xa为对应状态向量的先验参数向量,Sa为先验误差协方差矩阵。反演问题就是寻找价值函数式(16)的极小值,即x相对于y的概率最大化。

ORAC反演算法利用列文伯格-马夸尔特(Levenberg-Marquardt,LM)数值优化方法对价值函数进行极小化处理,即为一个优化迭代的过程。若y是m维的,x是n维的,则x的迭代表达式为

式中:K为雅可比(Jacobian)矩阵,下标k为当前迭代的次数,γ为可变参数,D为一个n×n的对角尺度矩阵。K为一个m×n的矩阵,其中的元素形式为

γ是关系到LM算法效率和鲁棒性的一个关键参数,在迭代过程中,ORAC采用因子10来增大或减小γ的大小。矩阵D是用来保证状态参数在相同的量级范围内,以保证数值计算的稳定性。

在LM优化算法中,当x满足一个收敛标准或当迭代次数已超出一个设置的最大值时(即迭代不收敛),迭代过程停止。与传统反演方法相比,最优估计方法具有两个明显的优势:一是将先验知识和观测信息一起包含在反演过程中,将反演问题转换为一个最优化问题,利用尽可能多的信息来约束反演;二是将前向模型和参数的多源误差一起纳入到反演过程中,通过为反演误差提供约束和对误差传递的严格控制,在算法上保证了反演的精度。

2 多角度气溶胶反演

2.1 星载多角度标量反演算法

MISR传感器的全称为多角度成像光谱仪(Multiangle imaging spectroradiometer),与MODIS传感器一同搭载在Terra卫星上,于1999年12月升空。MISR是可提供多角度、连续、高空间分辨率的传感器,其主要功能是全球多角度地形地物数据的获取。MISR 采用 0°、±26.1°、±45.6°、±60.0°和 ±70.5°(符号“±”表示前向和后向)共9个不同的角度进行对地观测,每个角度都设有4个通道(446.4、557.5、671.7、866.4 nm),共36个通道,可获得275、550、1100 m空间分辨率的观测[19]。

MISR气溶胶反演算法采用的气溶胶模型由一系列不同的气溶胶组成,每一种气溶胶混合物最多可以包括3种不同的气溶胶组分,粒子谱分布采用对数谱分布。MISR多角度气溶胶反演算法假定地球表面每个17.6 km×17.6 km的区域内大气气溶胶在水平方向上成分均一,并包含256个1.1 km×1.1 km的子区域。将一个像元表观辐亮度可以表示为程辐射和地表贡献两部分,选取区域中所有样本的最小值作为偏移等效反射率,然后定义散射矩阵(协方差矩阵),将协方差矩阵分解成特征值和特征向量的形式,进而利用经验正交函数分解(Empirical orthogonal function,EOF)方法实现地气的分离[54-57]。

对于EOF方法,假设一个像元的观测辐亮度信息Lx,y,λ[58]可以表示为

式中i和j表示MISR卫星传感器9个不同的观测角度的编号。通过EOF处理,协方差矩阵Cλ可表示为特征值和特征向量组合的形式,则地表贡献部分可以表示为

式中Ax,y,λ为对应的系数。对应地,

因此,在17.6 km×17.6 km的区域范围的等效辐亮度可以表示为

其中的系数

在此基础上,通过对设定的代价函数取极小,即可实现气溶胶参数的反演。

在卫星星下点的观测角度上,由于亮地表等原因,陆地上空的轻度雾霾很难被探测到,然而MISR通过多角度观测可以有效探测轻度灰霾[20]。这是因为大气程辐射的光程距离越大,灰霾对程辐射的影响就越大,观测结果也就越明显。在ATSR的气溶胶参数反演工作中,可考虑借用MISR所采用的混合气溶胶模型,但是由于角度太少,只有两个特征值,经验正交函数方法不适用。同样,MISR由于缺少1.6 μm的波段设置,也不利用ADV方法进行气溶胶参数的反演。

Zhang等[59]将EOF方法推广到了具有多角度探测能力的法国“太阳伞”(Polarization and anisotropy of reflectances for atmospheric science coupled with observations from a lidar,PARASOL)卫星POLDER(Polarization and directionality of the earth′s reflectances)传感器观测数据的反演处理上,实现了POLDER多角度地表反射贡献的估算。同时结合6SV辐射传输模式建立气溶胶反演查找表,最终实现了多光谱和多角度约束条件的AOD的反演。此外,与从传感器设置方面设计多角度气溶胶反演算法的思路有所不同,还有一些基于卫星单角度观测合成多角度观测的气溶胶反演算法方面的研究,如Lyapustin等[60-62]所发展的MAIAC(Multiangle implementation of atmospheric correction)算法和相关评估应用[63,64]以及Xue等[65]所提出的SRAP(Synergetic retrieval of aerosol properties)算法。

2.2 星载和机载多角度偏振反演算法

法国的POLDER系列传感器可以获取多光谱、多角度的辐射(标量)和偏振观测信息,在提高大气参数反演精度方面取得了较好效果。一般情况下,卫星在大气层顶观测的偏振信息主要来自大气,地表的偏振反射贡献相比大气程辐射的偏振贡献要小很多[66,67];通过对地表偏振反射及瑞利(Rayleigh)散射偏振贡献部分的进一步订正,可实现对陆地上空气溶胶参数的反演[68-75]。其中,Dubovik等所发展的GRASP(Generalized retrieval of aerosol and surface properties)方法是国际上广泛认可的、同最具代表的最优化反演模型和算法[11,76-80]。在多角度偏振卫星遥感反演过程中,代价函数一般采用多角度偏振表观反射率的模拟值与真实观测值的残差取极小的形式,由于残差在多角度观测时的累积效应导致最优模型选取出现较大偏差,会进一步影响气溶胶的反演精度。Zhang等[72]提出了一种基于残差分类排序(Grouped residual error sorting,GRES)的细颗粒物反演新策略,通过对模拟残差、光学厚度分类排序和查找,有效减少了各类误差的多角度累计效应,实现最接近真实大气细颗粒物模型和参数的优选。Fan等[81]基于多角度偏振观测和神经网络模型对海洋上空的气溶胶参数进行了反演处理和结果验证,阐明了神经网络可用于海洋上空气溶胶的业务化反演。

继法国POLDER系列传感器之后,我国高分五号(GF-5)卫星搭载的大气气溶胶多角度偏振探测仪(Directional polarimetric camera,DPC)已成为全球气候变化、污染传输和气溶胶特性研究领域重要的数据来源[82-85]。Li等[86]针对DPC提出了一套基于最优化的反演理论框架,综合利用DPC多角度的强度和偏振观测量进行气溶胶和地表多个参数的联合反演,并利用UNL-VRTM(Unified linearized vector radiative transfer model)系统地进行了信息量分析方面的研究[87,88]。这里,最优化反演建模所采用的代价函数[86]为

对应的梯度向量为

式中:γ1和γ2表示对应观测部分和先验部分的系数,用来进一步调整两部分对代价函数的贡献比率,其他的参数定义与式(17)一致。要进行多参数联合反演,等价于如下形式的最优化问题[89]

式中:l和u为对应的边界约束条件,需要由先验知识来确定,并采用含有约束的拟牛顿算法L-BFGS-B进行优化模型的迭代求解[90]。

为了减少每个波段需要反演的地表核函数的系数,对核驱动的地表双向反射分布函数(Bidirectional reflectance distribution function,BRDF)形式进行改进和简化[91],即

式中:fgeom和fvol为与观测几何有关几何光学和体积表面散射的核函数,系数k(λ)与波长有关,而k1和k2与波长无关,这样针对DPC传感器443、490、565、670、865 nm的5个波段的观测就需要反演7个BRDF核参数,需要反演的BRDF参数数量可大大减少。针对地表偏振反射贡献的描述,引入单参数的地表偏振双向反射分布函数(Bidirectional polarized reflectance distribution function,BPDF)模型[92]。因此,包含反演参数的状态向量可以表示为

式中:τa(λ)表示可进一步反演得到的AOD,Qext为对应的消光效率因子,可由UNL-VRTM模式所包含的Mie散射模块计算得到。

在此最优化反演框架的基础上,Zheng等[94]利用POLDER多角度偏振观测数据反演了细、粗模态气溶胶体积柱浓度和光谱的气溶胶光学厚度,反演结果与POLDER官方产品具有较好的一致性。Ge等[95]利用地表偏振反射率基本不随波长变化特性,发展了一种新的细模态气溶胶光学厚度反演算法,有效地避免使用BPDF模型带来的估算误差,实现了气溶胶参数和地表偏振反射率的同时反演。此外,基于SONET(Sun-sky radiometer observation network)16个站点7年连续的地基观测[8],Li等[96]建立了10种“基因”模态气溶胶模型,通过组合后可表述中国区域94%以上的气溶胶特征,显著优于国外模型,为反演中所需的普适气溶胶模型提供了关键参考。针对机载偏振传感器的观测,基于传统的查找表方法,Waquet等[97,98]开展了针对MicroPOL(Micro polarimeter)气溶胶反演研究,Wang等[99,100]和Qie等[101]分别开展了针对我国AMPR(Advanced atmosphere multi-angle polarization radiometer)传感器观测的气溶胶光学厚度反演研究。基于最优估计方法,Waquet等[102]和Wu等[103]分别对RSP(Research scanning polarimeter)进行了气溶胶多参数的反演,Xu等[104-106]对AirMSPI(Airborne multi-angle spectropolarimetric imager)进行了反演,并分别利用地基观测数据对气溶胶反演参数进行了精度验证。之后,Puthukkudy等[107]对AirHARP(Airborne hyper angular rainbow polarimeter)进行了气溶胶多参数的反演,Fu等[108]对SPEX(Spectropolarimeter for planetary exploration)机载观测进行反演,进一步阐述了多角度偏振观测在气溶胶多参数遥感反演上的独特优势。另外,针对美国TMEPO(Tropospheric emissions:monitoring of pollution)静止卫星计划的GEO-TASO(Geostationary trace gas and aerosol sensor optimization)机载高光谱传感器,Hou等[109-111]基于主成分分析(Principal component analysis,PCA)进行光谱重建的方法,综合利用尽可能多的高光谱观测信息,建立了面向静止高光谱卫星遥感的气溶胶和地表多参数最优化反演理论框架和反演算法,并利用地面实测数据对光谱的AOD和细粒子AOD反演结果进行了验证。

3 结论

针对国际上多角度(含双角度)及偏振卫星遥感反演算法的最新研究进展进行总结和综述,可为我国多角度及偏振遥感反演算法的设计和发展提供有效的参考和关键的支撑。ATSR系列传感器采用的ADV算法利用两个观测角度地表反射率的比值关系进行地气解耦合,避免了对先验地表反射率库的依赖,实现了植被等暗地表条件下的气溶胶参数反演。SU算法是一种基于地表经验模型的方法,通过引入双角度观测条件下地表反射率模型的先验约束,进而实现气溶胶参数的反演。ORAC算法将先验知识和观测信息一起包含在反演过程中,将反演问题转换为一个最优化问题,利用优化方法进行迭代求解,保证了反演的精度。这三种典型的双角度气溶胶反演算法的不同在于气溶胶模型的选取和地气解耦合方法上的处理,从而为双角度气溶胶遥感反演提供了不同的解决方案。MISR多角度气溶胶反演算法结合9个角度的观测信息,通过设定不同的混合气溶胶模型,利用经验正交函数的方法实现地气的分离,进而反演得到AOD等气溶胶参数。偏振遥感通过与多光谱和多角度的观测方式相结合,使得对气溶胶的观测信息更加丰富,更有助于气溶胶光学和微物理参数的提取。基于国际上新发展的最优估计反演框架,多角度偏振反演算法可充分利用多光谱、多角度、标量和偏振的观测信息来约束反演,实现气溶胶和地表关键参数的联合反演,为大气气溶胶定量遥感反演提供了新的途径。