郭 强，唐 家 奎*，何 文 通，田 媛，于 新 菊,2

(1.中国科学院大学,北京 100049;2.中国科学院遥感与数字地球研究所,北京 100101)



利用MODIS可见光波段反演陆地气溶胶光学厚度

郭 强1，唐 家 奎1*，何 文 通1，田 媛1，于 新 菊1,2

(1.中国科学院大学,北京 100049;2.中国科学院遥感与数字地球研究所,北京 100101)

利用Herold等建立的地表反射率库及MODIS遥感影像研究城市区和非城市区典型地物在可见光红蓝波段地表反射率的比值特性。在此基础上，利用MODIS 1 km分辨率遥感影像红蓝可见光波段实现了气溶胶光学厚度的反演,采用卫星过境时间前后半小时北京和香河AERONET站的气溶胶光学厚度观测平均值作为验证参考。结果显示，66.67%的反演结果处于±0.05±0.15 τ的误差界限内,反演算法不受地表反射率的限制，而且只利用了可见光红蓝波段，避免缺少近红外波段数据的限制。

气溶胶光学厚度；遥感；MODIS；地表反射率

0 引言

大气气溶胶作为大气中最活跃的3种因素之一，它的辐射强迫效应深刻影响着全球的辐射能量平衡及气候变化[1]。气溶胶颗粒在大气中的停留时间，一般为几天到几周，而且其特征随空间和时间都有显著的变化。目前研究确定大气气溶胶对气候系统的精确辐射强迫作用仍十分困难,而准确的气溶胶监测对于研究全球辐射能量平衡及环境监测保护具有重要意义。

卫星遥感具有大尺度同步观测、高时效性等特点，使之成为监测大气气溶胶的有效手段。卫星传感器接收的信号来源于地表和大气两部分。利用卫星影像反演气溶胶特性首先需要从卫星接收信号中区分地表贡献和大气贡献，即实现地气解耦。不同卫星传感器利用卫星观测值的光谱、角度及极化等特性实现地气解耦。Kaufman等[2,3]基于地表反射率特征提出暗目标法，其基本原理是在浓密植被区及暗色土壤区，可见光的红、蓝波段地表反射率较低且与短波红外波段存在线性相关关系。在短波红外波段气溶胶的散射吸收作用很小，假设气溶胶光学厚度为0，卫星观测到的大气层顶表观反射率就是地表真实反射率。通过短波红外波段的反射率推断可见光波段的地表反射率，实现地气解耦。该方法局限于暗目标区域，对于城市、干旱和半干旱等高地表反射率区域无法获得气溶胶光学厚度(Aerosol Optical Depth,AOD)。此外，算法需要短波红外通道，而目前在轨的遥感卫星大多没有短波红外通道，尤其是高分卫星，限制了该算法的应用。Tanré等[4,5]假定一定范围内复杂地表上空大气透过率相同，提出了结构函数法以解决高地表反射率对于气溶胶监测的限制。结构函数法较少受到地表反射率大小的限制，但由于对几何校正的要求较高，业务化运行比较困难。Deuzé等[6]从POLDER卫星观测获得的大气后向散射偏振信息中分离出气溶胶贡献，进而反演得到气溶胶光学厚度,但偏振的方法只能反演细粒子气溶胶。Tang等[7]基于AQUA和TERRA数据提出了双星协同反演算法，实现了气溶胶光学厚度和地表反射率的同时反演，部分解决了亮地表区域反演无效的难题，但使用两颗卫星的数据反演的光学厚度会受到传感器的性能退化、不同传感器影像之间图像匹配精度等影响。Knapp等[8]假设地表反射率特征在一段时间内保持不变，而且至少有一天是不受气溶胶影响的清洁日。考虑到残云的影响，选取每个像元位置第二最低的反射率值合成一幅参考影像。假设一个背景AOD并用其校正参考影像，得到参考的地表反射率用于反演AOD，算法的误差主要来源于地表反射率和背景AOD的假设。Hsu等[9,10]根据在红、蓝波段AOD对天顶辐亮度有显著的贡献提出了基于地表反射率库的深蓝(Deep Blue)算法，并成功应用于撒哈拉沙漠、阿拉伯半岛等干旱半干旱地区，但该方法需要事先建立地表反射率库作为先验知识。Diner等[11]假设地表辐射亮度的角度分布与波长无关，利用多角度观测数据消除地表反射率的影响，反演气溶胶光学厚度。Wong等[12]利用模型模拟将卫星观测到的表观反射率解析为气溶胶反射率和大气瑞利散射引起的路径反射率，在此基础上利用最小反射率技术(MRT)确定季节性的地表反射率,然后利用SBDART辐射传输模型构建查找表LUT，反演得到香港地区500 m分辨率的AOD。该反演算法有希望应用于其他区域的城市/局部尺度高分辨率AOD反演。Remer等[13]提出了3 km分辨率的AOD反演算法，并将加入MODIS C006版本的AOD反演算法中。MODIS 3 km反演算法相比10 km反演算法缩小了反演窗口，精度也略差于10 km分辨率算法。但提高分辨率的AOD产品可以为相关领域的研究提供重要的工具。

本文基于MODIS数据探讨北京及周边不同地表覆盖类型区域的可见光红蓝波段地表反射率比值特征，并尝试利用可见光红蓝波段卫星观测值反演气溶胶光学厚度，从而摆脱反演算法对近红外波段的依赖。

1 数据与算法

1.1 MODIS数据

采用的数据包括卫星数据和地面站点数据。卫星数据包括MODIS Level 1B地表反射率数据、Level 2气溶胶光学厚度数据、Level 3合成的地表反射率数据和Level 3土地覆盖类型产品。MODIS是搭载于EOS/TERRA和EOS/AQUA太阳同步轨道系列卫星的主要传感器。TERRA的过境时间为上午10∶30，AQUA的过境时间为下午13∶30。MODIS具有从可见光到远红外的36个光谱通道，地面分辨率分别为250 m、500 m和1 000 m，扫描幅宽2 330 km。反演采用1 km分辨率的地表反射率数据，时间范围为2009年4-9月。Kaufman等[2]指出，遥感影像的空间分辨率为250 m及更高时，大气点扩散函数引起的临近像元影像显著。当空间分辨率为1 km甚至更低时，临近像元效应很小，在大气校正过程中临近像元效应带来的误差远小于气溶胶特征的不确定性带来的误差。本文选取1 km分辨率的数据，因此可以不考虑临近像元效应的影响。

1.2 AERONET数据

地面站点数据来源于AERONET站点AOD数据。AERONET是NASA建立的覆盖全球的地基气溶胶监测网络。站点观测设备主要采用法国CIMEL公司的偏振太阳光度计CE318-II和标准太阳光度计CE318-I，观测波段主要有8个：0.340 μm、0.380 μm、0.440 μm、0.500 μm、0.675 μm、0.870 μm、0.936 μm、1.020 μm。AERONET气溶胶光学厚度数据分为3个等级：1级数据未做云去除；1.5级数据仅做了云去除；2级数据既经过云去除又经过人工检查。本文采用其官网(http://aeronet.gsfc.nasa.gov)发布的2级AOD数据对卫星反演的气溶胶光学厚度进行验证。AERONET 2级AOD数据的精度为±0.02[14]，因此可以认为是AOD的真值。AERONET站点在0.55 μm波长没有AOD观测值，本文采用Angstrom公式拟合得到。

τλ=β·λ-α

(1)

式中:β为大气浑浊度指数，α为Angstrom指数，λ为波长，τλ为λ对应的气溶胶光学厚度。

采用北京站(位于城市区域)、香河站(位于非城市区域)2个AERONET站点验证反演的结果。

1.3 反演原理与算法

Flowerdew等[15]研究指出，地物对于不同波长的反射特征取决于地物的微观结构，即地物的分子、原子与电磁波的作用机理。因此对于不同的波长，地物BRDF的大小不同，但地物反射率的角度分布特征即BRDF的形状特征主要取决于地物的宏观结构，与波长的相关性很小，几乎可以忽略。因此有：

(2)

式中:rm、rn分别表示两个不同观测几何条件的地表反射率，m、n分别表示两个不同的观测几何条件；b1、b2表示两个不同的波段。

由上式可以导出：

(3)

式中:比例系数k为常数。

通过以上分析可知，对于同一种地物，在两个观测几何条件下两个波段反射率的比值相同。Levy等[16]通过对全球AERONET站点地表反射率的研究发现，在暗像元区可见光红蓝波段反射率的相关性比可见光蓝波段与短波红外波段的相关性高，而且更加稳定。在MODIS 官方C005及以后版本的反演算法中用以下两个等式实现地气解耦：

(4)

Levy等[16]对AERONET站点的研究发现，全球范围内红蓝波段地表反射率比值恒定在1.923附近，夏季城市区比较低(1.305左右)。因此考虑利用这个特性实现AOD的反演。

根据Vermote等[17]的理论，在地表朗伯体、大气水平均一的假设条件下，卫星平台观测到的大气顶部反射率可以表示为：

(5)

式中:ρ0是大气路径辐射项等效反射率，ρs是地表反射率，T(μs)、T(μv)分别是入射方向和观测方向的大气透过率，S为大气下界面的半球反射率。T(μs)、T(μv)总是同时出现，因此将T(μs)T(μv)作为一个参数考虑;未知数ρ0、S、T(μs)T(μv)与AOD、气溶胶粒子单次散射率、相函数、不对称指数等有关。

在实际反演过程中，首先针对试验区域确定一个气溶胶模型，然后用6s、MODTRAN等辐射传输模型模拟在不同观测几何情况下AOD与上述3个参数之间的对应关系，据此建立查找表得到AOD。MODIS红、蓝波段地表反射率之间存在如下关系：

(6)

(7)

不同区域地表覆盖类型不同，假设可见光红波段与蓝波段地表反射率比值在城市区、植被区固定。采取3种方法验证：1)根据Herold等[18]提出的地表反射率库，选取包括植被区、高亮地表、城市区等典型地物，用MODIS波段响应函数拟合地表反射率，考察红、蓝波段的反射率比值；2)选取一系列无云的MODIS影像，提取香河和北京站点位置的表观反射率，根据AERONET的观测值插值得到0.55μmAOD，然后用6s校正，得到香河和北京站点的地表反射率，最后考察反射率比值；3)选取春、夏、秋、冬4幅无云影像的MODIS地表反射率产品MOD09，考察地表反射率比值。具体反演流程见图1。

2 气溶胶光学厚度反演

2.1 红、蓝波段地表反射率比值估计

图1 算法反演流程

Fig.1Flowchartofalgorithmprocessing

Herold等[18]利用AVIRIS数据研究了城市地区典型地物的光谱特征，并建立了城市地区包括植被、土壤、沥青道路、水泥道路、屋顶等多种典型地物的光谱反射率数据库。AVIRIS数据具有0.4～2.5 μm范围内224个光谱波段，空间分辨率为4 m，光谱分辨率为10 nm。从上述光谱数据库中选取了9种典型地物类型，用MODIS红蓝波段响应函数拟合得到反射率，计算红、蓝波段地表反射率比值与蓝波段地表反射率的关系，结果如图2所示(见封3)。

从图2中蓝色椭圆内的地物可以看出，裸土、绿色植被、不进行光合作用的枯萎植被大多集中在1.4～2.1的范围内。城市周边主要地物是植被和高亮的裸土，因此可以推断在城市周边区域像元的红、蓝波段反射率比值处于1.4～2.1。城市区域的典型地物选取了3种典型的屋顶、3种道路。城市地区地物在蓝波段的地表反射率分散很大，例如水泥道路的反射率分散在0.13～0.46的范围内。红、蓝波段地表反射率比值分散也较大，分布在1.0～3.2之间。城市地物反射率比值的差异性反映了城市区地物的复杂性。由于选取的MODIS影像分辨率为1 km，在城市区多为混合像元，且受到临近像元的影响，因此需针对MODIS数据特点考察不同地物地表反射率特征。

为了考察MODIS影像上地物地表反射率特征，选取2009年北京周边61景无云的MODIS 1 km分辨率影像。为了降低地表BRDF的影响，选取的影像在AERONET北京站和香河站的太阳天顶角需小于30°。对于每一景影像，提取AERONET香河站和北京站的表观反射率，用卫星过境时间AERONET站点观测的气溶胶光学厚度进行大气校正，得到红、蓝波段地表反射率，进而得到地表反射率的比值(图3)，可以看出，AERONET站点多数比值处在1.65～2.25之间。

图3 AERONET站点红蓝波段地表反射率的比值

Fig.3 The ratio of red/blue band at the AERONET sites

为了进一步考察研究区域地表反射率比值特性，选取了四季的MOD09地表反射率产品，共4景影像。影像的成像时间分别为春分、夏至、秋分和冬至日，空间范围为北京周边400×400 km。用MODIS Level 3土地覆盖类型数据MOD12将研究区分为城市区和非城市区，分别考察不同区域地表反射率特征。对于每种区域，在影像上随机选取2 000个样点，考察地表反射率比值，结果如图4所示。

图4a-图4d分别为城市区四季红蓝波段地表反射率，图4e-图4h分别为非城市区四季红蓝波段地表反射率。可以看到，城市区和非城市区可见光红蓝波段地表反射率都具有很强的相关性，相关系数都在0.86以上。图4中的直线是回归线，颜色条显示的是散点的密度。可以看到绝大多数点分布在回归直线两侧很小的范围，显示了回归关系的稳定性。此外，除了春季非城市区相关性都大于城市区，这与城市区地表覆盖复杂有关。四季城市区和非城市区的回归直线的斜率如表1。

图4 MOD04红、蓝波段地表反射率关系

Fig.4 The relationship of red/blue band from MOD04

表1 红蓝波段地表反射率回归直线的斜率

Table 1 Slope of the regression line between red and blue band reflectance

春夏秋冬城市区1．93741．90921．94581．9105非城市区1．98871．92011．95571．9520

2.2 气溶胶光学厚度反演实例

为了验证算法的可行性，选取2009年北京周边4-9月AERONET香河站和北京站上空无云的53景MODIS 1 km影像进行反演实验。2009年5月3日和6月20日两景遥感影像中研究区域上空云量很少。图5(见封3)是两个日期研究区域影像真彩色合成。

反演前需对遥感影像进行辐射定标、重投影、云去除、裁剪等操作。云像元的去除采用红波段阈值法。Ackerman等[19]通过研究发现，在MODIS影像红波段云像元的表观反射率一般大于0.18，并将该值作为MODIS云检测产品MOD06的陆地区域默认阈值。因此采用0.18作为云像元去除的阈值。查找表是基于6s辐射传输模型设定不同的观测几何条件计算得到的。6s辐射传输模型可以模拟0.25～4.0 μm波长范围内卫星接收的信号。它采用连续阶散射(Successive Orders of Scattering)方法求解辐射传输方程，考虑了吸收气体(水汽、二氧化碳、甲烷、臭氧、二氧化氮等)的吸收、散射及气溶胶粒子的多次散射。设定条件包括几何条件、光谱波段、大气模式、气溶胶模式。选择数据的时间都是4-9月之间，因此大气模式选取中纬度夏季，气溶胶模式选择大陆性气溶胶。

图6(见封3)是两个日期反演得到的1 km分辨率气溶胶光学厚度。可以看出，两个日期北京城区上空的气溶胶光学厚度都明显大于周边地区，另外，天津也是气溶胶光学厚度的高值带。从气溶胶光学厚度的区域分布看，两个日期都呈现西北低、东南高的分布趋势。北京西北部山区植被茂密、人口稀少，气溶胶光学厚度小，大气污染相对较轻；北京城区及其东南方向平原地区人口密集，气溶胶光学厚度相对较大，反映了这些地区大气污染严重。6月20日北京城区东北部的山区上空有一自东向西的气溶胶光学厚度高值分布带，可能是由周边高污染地区的污染物通过风力的作用漂移、聚集形成。研究区域地表覆盖多样，包含了城市、浓密植被、高反射裸土等地物，反演结果较好地反映了不同类型区域上空气溶胶粒子的空间分布特征。

为了验证反演结果，获取两景影像卫星过境时间北京和香河AERONET站点观测的气溶胶光学厚度，与反演结果对比，结果如表2所示。5月3日的反演结果与AERONET观测值很接近，两个站点最大的误差为0.088。6月20日的反演结果与AERONET结果有一定的偏差，两个站点最大的误差达0.173。

表2 反演结果与AERONET观测结果对比

Table 2 Retrieved AOD and observed value from AERONET sites

北京站香河站反演AERONET误差反演AERONET误差2009-05-030．2370．2630．0260．2020．2900．0882009-06-200．3600．4980．1380．2560．4290．173

3 反演结果验证

为了进一步验证算法的精度，选取2009年4-9月香河和北京AERONET站点上空无云的75景影像进行了反演试验。同时选取卫星过境时间前后半小时之内AERONET站点的观测值取平均，得到光学厚度值。两者的对比如图7(见封3)所示。

图7中橙色实线为回归线，黄色虚线为误差界限±0.05±0.15τ，灰色虚线为1∶1线。总的来说，算法反演结果与AERONET观测结果有较好的相关性，75个有效反演值中，有66.67%的反演结果处于±0.05±0.15 τ的误差界限内。图8为反演值和AERONET站点观测值的时间序列散点图,可以看到，AERONET北京站和香河站上空反演结果与观测值的时间序列比较吻合。

图8 反演结果与AERONET观测值对比散点图

Fig.8 Scatter plot of AOD retrieved and observed value from AERONET

4 结论

本文提出了一种利用MODIS影像可见光波段数据反演高分辨率气溶胶光学厚度的算法，利用地表反射率库与MODIS遥感影像分析了北京周边不同典型地物类型可见光红、蓝波段地表反射率的比值特征，实现了气溶胶光学厚度的反演。选取了2009年4-9月北京和香河AERONET站点上空无云的多景遥感影像进行了反演实验，并将反演结果与AERONET站点观测值进行了对比验证。结果显示，有66.67%的反演结果处于±0.05±0.15 τ的误差界限内。

本文反演算法不受地表反射率的限制，可应用于包含城市区、高亮裸土等典型地物的广大区域。另外，反演算法只利用了可见光红蓝波段，有潜力应用于包含可见光红蓝波段的高分辨率遥感卫星数据。但需要指出的是，算法的反演需要提前估计大气模型和气溶胶模型，这是算法的误差来源之一。

[1] SOLOMON S,QIN D,MANNING M,et al.Climate Change 2007——The Physical Science Basis.Contribution of Working Group I to the Fourth Assessment Report of the IPCC[M].London：Cambridge University Press,2007.

[2] KAUFMAN Y J,WALD A E,REMER L A,et al.The MODIS 2.1 μm channel-correlation with visible reflectance for use in remote sensing of aerosol[J].Geoscience and Remote Sensing,1997,35(5)：1286-1298.

[3] KAUFMAN Y J,TANRÉ D,GORDON H R,et al.Passive remote sensing of tropospheric aerosol and atmospheric correction for the aerosol effect[J].J.Geophys.Res.,1997,102(D14)：16815-16830.

[4] TANRÉ D,DESCHAMPS P Y,DEVAUX C,et al.Estimation of Saharan aerosol optical thickness from blurring effects in thematic mapper data[J].J.Geophys.Res,1988,93(D12)：15955-15964.

[5] HOLBEN B N,VERMOTE E,KAUFMAN Y J,et al.Aerosol retrieval over land from AVHRR data——Application for atmosphere correction[J].IEEE Transactions on Geoscience and Remote Sensing,1992,30(2)：212-222.

[6] DEUZÉ J L,BRÉON F M,DEVAUX C,et al.Remote sensing of aerosols over land surfaces from POLDER/ADEOS-1 polarized measurments[J].Journal of Geophysical Research,2001,106(D5)：4913-4926.

[7] TANG J,XUE Y,YU T,et al.Aerosol optical thickness determination by exploiting the synergy of TERRA and AQUA MODIS[J].Remote Sens.Environ.,2005,94(3)：327-334.

[8] KNAPP K R,FROUIN R,KONDRAGUNTA S,et al.Toward aerosol optical depth retrievals over land from GOES visible radiances：Determining surface reflectance[J].International Journal of Remote Sensing,2005,26(18)：4097-4116.

[9] HSU N C,TSAY S C,KING M D,et al.Aerosol properties over bright-reflecting source regions[J].IEEE Trans.Geosci.Remote Sens.,2004,42(3):557-569.

[10] HSU N C,TSAY S C,KING M D,et al.Deep blue retrievals of Asian aerosol properties during ACE-Asia[J].IEEE Trans.Geosci.Remote Sens.,2006,44(11)：3180-3195.

[11] DINER D J,MARTONCHIK J V,KAHN R A.Using angular and spectral shape similarity constrains to improve MISR aerosol and surface retrievals over land[J].Remote Sensing of Environment,2005,94(2)：155-171.

[12] WONG M S,NICHOL J E,LEE K H.An operational MODIS aerosol retrieval algorithm at high spatial resolution,and its application over a complex urban region[J].Atmospheric Research,2011,99(3)：579-589.

[13] REMER L A,MATTOO S,LEVY R C,et al.MODIS 3 km aerosol product：Algorithm and global perspective[J].Atmos.Meas.Tech.Discuss.,2013,6(1)：69-112.

[14] HOLBEN B N,ECK T F,SLUTSKER I,et al.AERONET——A federated instrument network and data archive for aerosol characterization[J].Remote Sens.Environ.,1998,66(1)：1-16.

[15] FLOWERDEW R J,HAIGH J D.An approximation to improve accuracy in the derivation of surface reectances from multi-look satellite radiometers[J].Geophys.Research Letters,1995,22(13)：1693-1696.

[16] LEVY R,REMER L,TANRE D,et al.Algorithm for remote sensing of tropospheric aerosol over dark targets from MODIS：Collections 005 and 051：Revision 2[DB/OL].MODIS Algorithm Theoretical Basis Document,2009.

[17] VERMOTE E,TANRE D,DEUZJ L,et al.Second simulation of the satellite signal in the solar spectrum：An overview[J].IEEE Trans.Geosci.Remote Sensing,1997,35(3)：675-686.

[18] HEROLD M,GARDNER M,ROBERTS D A.Spectral resolution requirements for mapping urban areas[J].IEEE Transactions on Geoscience and Remote Sensing,2003,41(9)：1907-1919.

[19] ACKERMAN S A,STRABALA K I,MENZEL W P,et al.Discriminating clear-sky from clouds with MODIS[J].J.Geophys.Res,1998,103(D24):32141-32157.

Aerosol Optical Depth Retrieval over Land Using MODIS Visible Bands Imagery

GUO Qiang1,TANG Jia-kui1,HE Wen-tong1,TIAN Yuan1,YU Xin-ju1,2

(1.UniversityofChineseAcademyofSciences，Beijing100049; 2.InstituteofRemoteSensingandDigitalEarthofChineseAcademyofSciences，Beijing100101，China)

A surface reflectance library proposed by Herold et.al.(2003)and MODIS imagery are used to estimate the ratio of surface reflectance between visible red and blue bands in both urban and non-urban areas.On this basis,blue and red bands of MODIS 1 km resolution imagery are used to retrieve aerosol optical depth.The observations of Beijing and Xianghe AERONET sites are used to validate the results.In order to minimize the error,the averaged measurements within half an hour before and after the MODIS overpass are used in the validation.The result shows that 66.67% of the retrievals are within the expected error envelope of ±0.05±0.15 τ compared to the AERONET observed AOD.The application of this algorithm is not subject to any limitations of surface reflectance.Moreover,the algorithm only utilizes the visible blue and red bands,without using NIR band.

aerosol optical depth;remote sensing;MODIS;surface reflectance

2014-04-17；

2014-09-08

国家973支撑项目(2013CB733402);中国科学院数字地球重点实验室开放基金项目(2011LDE015);环保公益行业科研专项(201309011);中国科学院大学校长基金项目；浙江省重大科技专项(2012C13011-2)

郭强(1989-),男,硕士研究生,从事大气遥感研究。 *通讯作者E-mail：jktang@ucas.ac.cn

10.3969/j.issn.1672-0504.2015.02.009

P407

A

1672-0504(2015)02-0038-06