气溶胶光学厚度备考复习

第2章 气溶胶光学厚度反演的原理和方法

气溶胶光学厚度（Aerosol Optical Depth ）简称AOD ，定义为介质的消光系数在垂直方向上的积分，描述的是气溶胶对光的消减作用[7]。它是气溶胶最重要的参数之一，表征大气浑浊程度的关键物理量，也是确定气溶胶气候效应的重要因素。。通常高的AOD 值预示着气溶胶纵向积累的增长，因此导致了大气能见度的降低。现阶段对于AOD 的监测主要有地基遥感和卫星遥感两种方法。其中地基遥感又有多种形式：多波段光度计遥感、全波段太阳直接辐射遥感、激光雷达遥感等。其中多波段光度计遥感是目前地基遥感研究中采用的最广泛的方法。美国NASA 和法国LOA-PHOTONS 联合建立的全球地基气溶胶遥感观测网AERONET 所使用的就是多波段太阳光度计（Sun/Sky Photomerers ），在全球共布设1217个站点长期观测全球气溶胶的光学特性，积累了大量的AOD 数据，并用作检测气溶胶光学厚度反演精度的标准。而近年来卫星遥感技术的快速发展，多种传感器被用来研究气溶胶特性，加上经济发展带来的大气污染问题使得利用卫星遥感资料反演AOD 成为热门课题。 2.1 气溶胶光学厚度反演的基本原理

大气光学厚度是指沿辐射传输路径单位截面上气体吸收和粒子散射产生的总消弱，是无纲量值。在可见光和近红外波段，它可以由下列公式计算得出：

)()()()()()(a 21m λτ+λτ+λτ+λτ+λτ=λτμωω (2-1)

其中)(λτ表示大气总的光学厚度，)(m λτ表示整层大气的分子散射光学厚度，)(1λτω表示氧气的吸收光学厚度，)(2λτω表示臭氧的吸收光学厚度，)(λτμ表示

水汽的吸收光学厚度，)(a λτ表示气溶胶光学厚度[21; 22]。

卫星遥感反演大气气溶胶是利用卫星传感器探测到的大气顶部的反射率，也称为表观反射率，可以表示为[23]：

F /L s s *

μπ=ρ (2-2) 其中，L 表示卫星传感器探测到的辐射值，F s 表示大气上界太阳辐射通量，

μs 表示太阳天顶角θs 的余弦值。ρ*

与地表二项反射率之间的关系可以表达为：

ρ

-θθφθθρ+

φθθρ=φθθρνννν*

s d s s a s *

s 1)

(T )(F ),,(),,(),,( (2-3)

其中，θν表示传感器天顶角，θs 表示太阳天顶角，φ表示太阳方位角和卫星方位角确定的相对方位角；),,(s a φθθρν表示由大气分子和气溶胶散射造成的路径辐射，它与地表状况无关；)(F s d θ表示地表反射率归一化为零时总的向下辐射通量，也可以称为总的向下透过率，由于气溶胶粒子对太阳光的吸收和散射作用，它的值小于1.0；)(T θν是向上进入卫星传感器视场方向的总透过率，S 是大气后向散射比。

在单次散射近似中，路径辐射),,(s a φθθρν与气溶胶光学厚度τa 和单次散射反射率ω0之间的关系如下[24]：

),,(s a φθθρν=μμφθθτω+

φθθρνννS

s a a 0s m 4)

,,(P ),,(

(2-4)

其中),,(s m φθθρν是分子散射造成的路径辐射，它取决与大气模式，μν表示传感器天顶角的余弦值，μS

表示太阳天顶角的余弦值。

在式（2-3）中，)(F s d θ、)(T θν和S 取决于ω0，τa 和),,(P s a φθθν。

假设地表是均匀朗伯表面，大气垂直均匀变化，将式（2-4）代入（2-3）得：

ρ

τω-τωθτωθφθθρ+

μμφθθτω+

φθθρ=φθθρννννν*

a a 0a a 0S a a 0S d s S

s a a 0s m s *

)P ,,(s 1)

P ,,,(T )P ,,,(F ),,(4)

,,(P ),,(),,( (2-5)

上式中),,(s φθθρν为假设的朗伯体特性的地表反射率，卫星传感器接收到的表观发射率),,(s *

φθθρν既是地表反射率),,(s φθθρν的函数，又是气溶胶光学厚度τa 的函数[25]。表观反射率),,(s *

φθθρν以及太阳和传感器的几何参数),,(s φθθν可以从卫星遥感资料中获取，

假如可以得知地表反射率),,(s *

φθθρν，并用气溶胶类型和大气模式来确定ω0和),,(P s a φθθν的相关参数，理论上就可以计算得出地面上空的气溶胶光学厚度τa 。反之，若已知地面上空气溶胶光学厚度τa 、气溶胶类型以及大气模式，也可以反演出地表反射率),,(s *

φθθρν[26]。

由式（2-3）可以得出反演气溶胶的最优条件是地表反射率低且光谱波段波波长较短。在地表反射率角度的情况下，气溶胶散射引起的路径辐射

)

,,(s a φθθρν（与地表状况无关）对表观反射率),,(s *

φθθρν起主要作用，此时反演气溶胶光学厚度误差较小；而在地表反射率较大的情况下，地表贡献项

ρ

-θθφθθρνν*

s d s s 1)(T )(F ),,(对表观反射率),,(s *

φθθρν影响较大，此时反演精度较低[27]。

2.2 气溶胶光学厚度反演的主要方法

自20世纪70年代中期开始，利用卫星数据反演气溶胶光学厚度的研究已经

有40年的历史，反演的方法有单通道算法、多通道算法、暗像元法、结构函数法、深蓝算法、多星协同反演法、海陆对比法、多角度偏振法、热辐射对比等[28]。目前有代表性的常用气溶胶光学厚度反演算法有两种：一种是通过路径辐射项求取气溶胶光学厚度的暗像元法，另一种是通过透过率求取气溶胶光学厚度的对比法。

暗像元算法英文全称为Dense Dark Vegetation ，简称DDV ，它是通过路径辐射项来计算气溶胶光学厚度。由于地表物体的复杂多样性造成反射率变化范围很大，很难从辐射值中分理处辐射项，如果想通过辐射项来获取气溶胶信息，就必须使地表辐射值较小且能确定其精确值，这样就能够最大限度的消除地表反射率的不确定性带来的影响。在卫星影像中，大量浓密植被区由于在可见光波段反射率极低(约为0.01~0.02)，它们被称作暗像元。1988年Kaufman 等利用大多数地物在红（0.60~0.68μm ）蓝(0.40~0.48μm)波段反射率低的特性,根据归一化植被指数（NDVI ）或近红外通道（2.1μm ）的表观反射率进行暗像元的识别，并假定已知这些暗像元在红蓝通道的地表反射率，依据一定的关系反演气溶胶光学厚度

[29]

。通过大量的卫星影像资料，考虑到多种地表覆盖物，拟合得到红（0.66μm ）

蓝（0.47μm ）和中红外通道（2.1μm ）地表反射率的关系[29]：

4

2*1.2blue

*

1

.2red ρ=

ρρ

=ρ (2-6)

对于2.1μm 通道卫星观测表观反射率几乎不受气溶胶影响，其值接近地表反射率，因此可以用2.1μm 通道的表观反射率代替地表反射率，根据式（2-6）计算出红蓝通道的地表反射率。再假定合适的气溶胶模型，就可以计算出气溶胶光学厚度。

然而对于干旱、半干旱以及冬季城市等高发射率地区，用暗像元法来反演气溶胶光学厚度还存在不少困难，在以上地区由于地表的非均一性使确定地表反射率的精确值十分不易[30]。当地表反射率升高时，气溶胶的指示作用降低[31]。一般情况下，当地表反射率较低时，传感器接收到的辐射值随着气溶胶的增多而迅速变大，暗像元法就是利用浓密植被在红、蓝波段的地表反射了和气溶胶的这种关系来反演光学厚度；随着地表反射率的不断增大，接收到的辐射值随气溶胶的增多而增大速度减缓，当地表反射率增大到某一程度时，辐射值将不随气溶胶的增多而增大，甚至会出现降低的趋势。为了使暗像元方法的应用更为广泛，2002年Kaufman 等通过大量的数据验证，对以前的方法进行了扩展：对于星下点，暗像元法适用于中红外通道表观反射率小于0.4的区域；对于非星下点，需要考虑太阳和卫星的几何参数，适用范围可以扩展到中红外通道表观反射率小于

)1

1(0.1250μ+μ的地区（其中μ为卫星天顶角的余弦值，0μ为太阳天顶角的

余弦值），红蓝通道的地表反射率依然遵循式（2-6）的函数关系。

扩展后的暗像元法是由Levy 和Remer 等人提出的称为V5.2的气溶胶反演算法[32; 33]。与原来的算法相比，它主要有两方面的改进：首先提出VIS/SWIR （可见光和中红外通道的反射率比值）不再是一个常量，而是一个函数关系；其次考虑了NDVI （植被指数）对于地表反射率的影响，红蓝通道和中红外通道的地表反射率关系如式（2-7）所示[34]：

)

(g )(f s

66.0s 47.0s 1.2s 66.0ρ=ρρ=ρ (2-7)

该方法中的VIS/SWIR 的变化跟不但跟地表植被指数有关，根据Remer 和Gatebe 等的研究表明它还与太阳天顶角、传感器方位角、散射角有关，尤以散射角的影响最大[35; 36]。它可以表示为：

)cos sin sin cos cos (cos 001φθθ+θθ-=Θ-

(2-8)

其中θ表示观测天顶角， θ0表示太阳天顶角，φ表示太阳方位角和卫星方位角的相对方位角。根据式（2-8）可将式（2-7）扩展为：

int y slope )(g int y slope )(f 66

.047.066.047.0s

66.0s 66.0s 47.01

.266.01.266.0s

12.2s

1.2s

66.0+ρ=ρ=ρ+ρ=ρ=ρ (2-9)

其中

005

.0y i n t 49.0slope 003.000025.0yint 27.0002.0slope slope 66.047.066.047.01.266.0NDVI 1.266.01.66.0R SWI ==+Θ=-Θ+= (2-10)

NDVI SWIR 可表示为：

)()(NDVI m

1.2m

2.1m 1.2m 2.1SWIR

ρ+ρρ-ρ=

(2-11)

其中，ρm

2.1和ρm

1.2分别表示MODIS 第5波段和第7波段接收到的辐射值，当NDVI SWIR 大于0.6时，表示该地区的植被状况较好；当NDVI SWIR 小于0.2时，表示该地区植被稀疏[37]。

当25.0NDVI SWIR <时，48.0slope NDVI 1.266.0R

SWI =； 当0.75NDVI SWIR >时，58

.0slope NDVI 1.266.0R SWI =； 当0.75NDVI 25.0SWIR ≤≤时，)

25.0NDVI (*2.048.0slope SWIR NDVI 1.266.0R SWI -+= 其中slope 和yint 为统计方程的经验系数，slope 为偏移量，yint 为截距。 将（2-10）代入（2-9）可得出： a ． 植被指数小于0.25的情况：

005

.049.0033

.000025.0)27.0002.048.0(s

66.0s 47.0s

1.2s

66.0+?ρ=ρ+Θ+-Θ+?ρ=ρ (2-12)

b ．植被指数大于0.75的情况：

005

.049.0033

.000025.0)27.0002.058.0(s

66.0s 47.0s

1.2s

66.0+?ρ=ρ+Θ+-Θ+?ρ=ρ (2-13)

c ． 植被指数介于0.25和0.75之间的情况：

005

.049.0033.000025.0]27.0002.0)25.0NDVI (2

.048.0[s

66.0s

47.0SWIR s

1.2s

66.0+?ρ=ρ+Θ+-Θ+-+?ρ=ρ (2-14)

对于旱季和中高纬度地区的冬季，植被覆盖少，地物反射率高，此时应用暗目标法会造成一定的误差，因此，针对陆地亮地表的情况，发展了结构函数法来反演气溶胶光学厚度。

结构函数法也称为对比法，是早期用来研究陆地污染气溶胶采用的卫星遥感算法。它跟暗像元法一样采用的是红、蓝通道数据。对比法需要假定在同一地区同一时间段内地表反射率不变[38]。然后找出污染极小或无污染的一天作为“清洁日”，以它的气溶胶数据作为参考反演“污染日”的大气气溶胶光学厚度。Tanré和Holben 在采用对比法进行气溶胶反演的时候，假定在同一卫星观测时刻，对于两个相邻的像元大气是均一的，并在反演的过程中引入了结构函数的概念，使得地表反射率的确定变得容易[39; 40]。该算法主要利用的是表观反射率的地表贡献

项来反演气溶胶光学厚度。在暗像元方法不适用的地区，结构函数法为气溶胶光学厚度的反演提供了一条新的途径。

在反演过程中，需要通过一段时间内卫星观测数据的分析，选定其中无污染（气溶胶光学厚度极小）的一天作为“清洁日”，然后通过地面观测或者其他途径确定这一天的气溶胶光学厚度值，作为背景的气溶胶信息。假定地表目标无变化，通过透射函数的变化就能获取其它“污染日”的气溶胶光学厚度。由于总透射函数对于相函数细节不敏感，结构函数法对散射相函数有较大的独立性，获取气溶胶光学厚度主要根据单次散射反照率和不对称因子。 由式（2-3）得出相邻两点的表观反射率差值为：

)θ)(θ(T V S )j i,(*

)j (i,ρ?=ρ? (2-15)

地表结构函数为[40]：

)d j i,j i,(d)n(m 1

)d (M 2n 1i d m 1

j 2ρ+-ρ∑∑-=

=-= (2-16) 改进后的地表结构函数为[41]：

∑∑ρ-ρ+ρ-ρ+ρ+-ρ-=

=-=++n 1i d m 1

j d j i,j i,2

j d,i j i,222

])()()d j i,j i,([

d)d)(m -3(n 1)d (M (2-17)

在卫星高度上的结构函数为：

)(T )(T )d (M )d (M V 2

S 222*θθ=

(2-18)

假定“清洁日”t 1和“污染日”t 2地表特征保持不变，则)t d,(M 12=)t d,(M 22，观测结构函数关系根据式（2-12）得：

)]

t (),t (),t ([T )]t (),t (),t ([T )t ,d (M )t ,d (M 2a 2v 2s 1a 1v 1s 2*

1*

τθθτθθ= (2-19)

由上式可得在已知“清洁日”的气溶胶光学厚度的情况下，“污染日”的气溶胶光学厚度就可以由卫星观测值获取。

近年来随着卫星遥感技术的不断发展，对地观测卫星携带的传感器种类不断增多，这对于不同平台之间的结合研究起到了推动作用。为了解决气溶胶光学厚度反演中未知

参数过多导致气溶胶光学厚度信息难以提取的问题，唐家奎等提出了基于Terra 和Aqua 双星的MODIS 遥感数据协同反演算法[17]。该算法具有以下几种优点：（1）不需要事先假设气溶胶的类型等参数。（2）不需要假设或利用其它统计资料预先估计地表反射率，双星协同算法是将真实的表观反射率作为反演气溶胶光

学厚度的变量之一与其同时反演。（3）不收地表反射类型的影响，可应用高反射率的陆地气溶胶反演。（4）可以同时获取两次过境时间的气溶胶信息，有利于分析气溶胶的动态变化规律[17]。基本原理如下：

在整个大气层范围内对微元体积大气的辐射性质进行积分称为地球—大气系统的辐射传输模型。根据以前学者的研究分析[42]，地球—大气系统之间的辐射传输方程可以表达为：

?λλλ

λσ+-ω''γλ'πσ=δρθ')r ,z (I )k (d )r ,r ,z (),z (I 4dz I cos （2-20）

其中θ'表示传感器天顶角，ρ表示空气密度，σ表示散射系数，z 表示高度，

)r ,r ,z ('γλ

表示散射相函数，ω'表示固体角，r 表示方位（天顶角，方位角），

)r ,z (I λ表示在z 高度r 方向的辐射强度。方程（2-20）中不能够得到)r ,z (I λ

的

解析解，只能够获得方程的近似解，Xue 等通过用多个不同变量将方程简化后得到地表反射率L 和表观反射率L '之间的关系如下[43]：

θ'

ξ-θ'ξ

-τ'-+-'τ'-+-'=

λ

λ

sec )b a (sec )b a (0

0e

)L 1(b )a b L (e )L 1(a )a b L (L （2-21）

其中ξ表示后向散射系数（通常取值为0.1），θ'=sec a ，b=2，τλ

0表示表

示气溶胶光学厚度。在双星协同反演模型中，仅考虑了空气分子和气溶胶粒子散射两种影响因素，将大气气溶胶光学厚度分为大气分子的瑞利散射和大气气溶胶粒子的散射两部分：

)()(L M 0∞τ+∞τ=τλ

λλ

（2-22）

根据Linke 等的研究大气分子的瑞丽散射近似解为[44]：

λ

=∞τ-λ00879.0)(09

.4M （2-23）

利用?ngstr?m 浊度公式可以求出气溶胶粒子散射光学厚度为：

λβ=∞τα

-λ)(L

（2-24）

将式（2-24）和（2-23）以及（2-22）代入（2-21）中，可以得出地表反射率L 以及?ngstr?m 浊度系数β和波长指数α之间的关系方程。

首先假设在两次卫星过境时间间隔内，地表反射率不变，且气溶胶类型保持不变，这样波长指数α不变，只有气溶胶粒子浓度变化，即浊度系数β是变化的。然后利用0.47,0.55以及0.66m

μ三个通道实现双星协同反演算法。

但该算法在实际应用中存在一些制约因素：（1）反演过程中很难获得方程的解析解，只能通过迭代法获取其近似解，计算量大，很难将其业务化。（2）有些像元出现不收敛的情况误差较大.（2）不同卫星传感器之间会存在图像的配准误差，反演过程中影响精度，导致部分反演结果误差较大。（3）不同卫星过境时间有一定间隔，气溶胶以及大气中其它成分的性质数量会发生变化，导致部分像元反演结果无法收敛。

20世纪70年代Kondratyev等提出了利用高反差地表发反演晴空条件下陆地上空的气溶胶光学厚度[45; 46]。利用该方法反演陆陆地上空气溶胶光学厚度需要在两个空间位置相近的区域，假设大气光学特性不变，选择明暗两种对比明显的像元进行反演。这里要求的明暗像元是相对而言，并不局限于浓密植被覆盖的特定区域，在一定程度上解决了亮地表上空反演的问题。但实际应用中，该方法要求传感器误差极小，且地表光谱反射值的差异要足够大，实现业务化反演比较困难。

2.3 本章小结

本章节主要介绍了利用遥感影像数据反演气溶胶光学厚度的基本原理，并介绍暗像元法、V5.2算法、结构函数法、双星协同反演算法以及深蓝算法和高反差地法。详细讲述了各自算法的原理，并描述了它们的应用范围、反演效果以及优缺点。总体来说，暗像元法应用最为广泛，它利用红、蓝以及中红外通道信息，根据它们之间地表反射率的关系进行气溶胶光学厚度的反演，但仅适用于地表反射率较低时的情况，相对来说结构函数法在针对城市及其周边地区的气溶胶污染监测有一定的优势。双星协同反演算法虽然不受地表类型限制，适用于城市等亮地表的气溶胶光学厚度反演，但由于该算法计算量大而且容易受到外界环境（如传感器间隔时间内温度、湿度）等的影响部分像元存在较大误差。高反差地法获取的反演结果只是相对气溶胶光学厚度，很难实现业务化。而V5.2算法即具有暗像元算法的优点，又在一定程度上减小了“亮”地表带来的影响。

本文的研究区域—北京市位于华北大平原北部，四季分明，在春、夏、秋季节城市及周边植被覆盖良好，冬季干旱城区地表反射率高，若利用暗像元法会存在一定误差，而扩展的暗像元法在一定程度上可以解决此问题。