一种不同植被覆盖条件下的气溶胶光学厚度反演方法

张 鑫，赵 祥，d，刘素红

(北京师范大学 a.地理学与遥感科学学院;b.遥感科学国家重点实验室;c.环境遥感与数字城市北京市重点实验室;d.北京师范大学全球变化与地球系统科学研究院，北京 100875)

1 研究背景

由于大气的存在，太阳辐射经过气体分子的吸收和气溶胶粒子的散射，得到减弱，同时部分散射信号直接或间接经过地物反射进入到传感器，又得到增强。反映在实际处理中，大气影响降低了图像的反差比，使图像可读性降低，增加了解译的困难。气溶胶光学厚度(AOD)是形容大气效应中气溶胶散射对大气影响的参数，是气溶胶粒子在整层大气中消光系数的积分。也是大气校正中最为重要的参数［1］。如何反演气溶胶光学厚度，对影像进行大气校正是至关重要的［2－3］。

到目前为止，遥感影像的大气校正方法很多。如Richter等人提出的直方图匹配法［4－5］，这种算法假定清晰区域和模糊区域的地表反射率的直方图是相同的，在一景影像中辨认出清晰和模糊的区域，匹配其直方图得到气溶胶光学厚度。Liang等人提出了聚类匹配法［6－7］，这种做法类似于直方图匹配法，假定每一种地物在不同大气条件下的平均反射率是相同的，对地物的划分更加精细。还有基于统计的不变目标法［8－9］，这种做法假定图像上存在一些反射率稳定的像元，这些像元在不同时相的影像中存在线性关系，当确定这个线性关系的时候，即可估算出气溶胶光学厚度。Kaufman等人提出的暗目标法是一种常用的大气纠正方法［10］。该方法的基本原理是假定需要校正的图像中存在黑暗像元区域(一般为浓密植被或者干净的水体)，假定其地表朗伯面反射，大气性质均一。在忽略大气多次散射辐照作用和邻近像元漫反射作用的前提下，反射率很小的黑暗像元受到大气的影响，使得其反射率相对增加，这个增加的部分可以认为是大气所产生的影响。当前大多传感器(MODIS［11］，AVHRR，ETM+等)使用暗目标法进行大气校正，这种校正方法简单、易行、自动化程度高，且精度较高。然而，暗目标法需要能够选取合适的暗目标像元，在校正中需要引入短波红外的数据(SWIR)来确定暗目标的位置［12－13］。

当前存在很多对地观测卫星，如SPOT，CBERS01/02 CCD，HJ－A/B CCD等数据，这些数据往往只提供可见光至近红外的光谱信息，这给暗目标法的大气校正带来了困难［14］，与此同时需要选取浓密植被的必要条件也在一定程度上限制了暗目标方法的应用。Guanter［15］等人提出一种针对ENVISAT/MERIS的大气校正方法，这种方法利用植被的光谱特性，针对不同植被覆盖度的像元在各个波段间的关系使用最优化算法耦合出气溶胶光学厚度，这给大气校正提出了新的思路。本方法同样基于不同植被覆盖条件下的遥感影像，通过建立方程式的方式求得气溶胶光学厚度。

2003年10月21日中巴地球资源卫星02星(CBERS02)成功发射，其CCD相机提供450nm至730nm的可见/近红外波段的影像，本文方法基于可见光至近红外波段的遥感信息，针对CBERS02星CCD影像进行大气校正，验证其校正结果并进行分析。

2 反演方法

根据空气的污染程度，可以设定在一定范围内气溶胶光学厚度变化相对保持不变，如在干净的影像中，可以设定10km×10km范围内保持不变，如果影像受气溶胶影像比较严重，可以设定1km×1km范围内保持不变。同时假定在这个范围内，植被与裸土类型是相同的。根据这个假设，在反演的过程中对遥感影像划分网格，在每个网格内，每个像元具有相同的气溶胶光学厚度，且纯净植被和裸土具有相同的光谱信息。对每个网格分别进行反演，得到光学厚度分布图。具体的操作流程如图1。

图1 反演操作流程Fig.1 Flow chart of the retrieval

如图1所示，首先，对遥感影像进行定标处理，得到每个像元的表观辐亮度。第二部对影像进行网格划分。在每个网格内，选取反演所需的像元。

2.1 筛选反演像元

本方法通过拟合不同植被混合比的像元间的关系来反演气溶胶光学厚度，因此，为了得到更好的结果，在每个网格中选取植被混合比(Cv)相差最大的3个像元作为反演的像元。本文假定像元中植被覆盖部分的植被类型一致且密度相同，使用等密度模型(dense vegetation model)来计算，如式(1)。

式中NDVI为归一化植被指数，其中NDVIsoil和NDVIveg分别为对应于裸土和高垂直密度植被的植被指数。Gutman文中使用NDVI来计算混合像元的植被混合比，然而大气中的水汽、臭氧、气溶胶、锐利散射等增加或减少红光和近红外波段的反射，这影响了NDVI在受到大气干扰条件下的精度，于是本文引入了抗大气植被指数(ARVI)来替代NDVI计算植被混合比。如式(2):

ARVI使用蓝光，红光和近红外3个波段的通道组合进行计算，可以在一定程度上消除气溶胶对植被指数的干扰［17－18］，ARVI的计算公式如下:

对比使用NDVI和ARVI计算在不同能见度下计算的Cv，使用模拟数据进行计算Cv为0～1时计算得到的植被混合比和真实值的误差和，得到表1。

表1 不同能见度下分别使用NDVI和ARVI计算不同Cv的均方根误差(RMSE)Table 1 The RMSerrors of Cv on different VIScalculated through NDVI and ARVI

由表可知，在能见度VIS＜23km时，NDVI由于受到大气影响，计算出的植被混合比误差均大于ARVI，能见度越小误差越大。通过使用ARVI计算植被混合比精度得到了提高。

2.2 反演气溶胶光学厚度

根据计算网格内所有像元的植被混合比，筛选出Cv＞0.5的像元。在其中选出3个植被混合比相差最大的像元(最大、最小、和中间值)作为反演气溶胶指数的像元。本方法认定其地表反射率是植被混合比的线性加成。如式(4):

式中ρveg和ρsoil为当前网格内植被像元与裸土的地表反射率。根据式(4)得到以下反演关系式:

式中ρ1，ρ2，ρ3分别为这3 个像元的地表反射率，由辐射传输方程得到，在假定地表朗伯的情况下，有式(6)［15］:

式中Lp，EdTv，S使用MODTRAN计算得到。为了简化计算过程，我们建立一个反射率计算查找表(LUTS)。查找表各个参数的设置如下:20个大气能见度(2～60km);9个太阳天顶角(10°～78°)。根据查找表可以很方便和快速地根据给定的气溶胶光学厚度计算出像元的地表反射率。其中气溶胶光学厚度和大气能见度是可以根据线性公式转换得到的。有以下MODTRAN模拟的经验公式［19］:

式中a，b在水汽含量为3g/cm2时的取值。

表2 气溶胶光学厚度和大气能见度转换系数Table 2 The conversion coefficients between AOD and VIS

由于地表反射率是由查找表计算得到的，我们不能简单的以解三元一次方程的方式对方程求解，于是使用一个一维迭代寻求最优值的方式，遍历所有AOD范围内的值，每引入一个AOD，根据式(5)中前2个方程得到ρveg和ρsoil，代入方程3中，求出方程左右两边之差δ:

当δ最小时，求得方程的最优解AOD。

最后，得到每个网格的气溶胶光学厚度分布图，使用一个3×3的空间滤波对结果进行平滑，得到最终的气溶胶光学厚度图。

3 结果与分析

本文分为2个部分对本文的算法结果进行分析:前半部分使用模拟数据来验证算法的精确性，后半部分使用本文算法对CBERS卫星数据的影像进行大气校正，检验其校正结果。

3.1 模拟数据分析

第一部分选择模拟数据对算法进行验证，我们获取植被和裸土的真实地表反射率(如图2)，由式(4)计算得到不同植被混合比下的地表反射率。使用查找表逆向计算这些地表反射率在不同气溶胶光学厚度下的表观辐亮度，通过得到的表观辐亮度使用本文的算法计算出气溶胶光学厚度。最后与真实的结果进行对比，以此来效验算法的效果。

图2 纯净植被和裸土的地表反射率光谱信息Fig.2 The spectral information of surface reflectance on pure vegetation and bare soil

首先选择3组Cv(0.2，0.5，0.8)的表观辐亮度进行计算，计算出的光学厚度与真实值进行比较，得到图3(a)。

如图3(a)所示，反演结果与真实值的均方根误差(RMSE)达到了0.14，拟合系数(R2)达到0.94，误差较大。为了提高反演的精度，我们另外使用了3组大于0.5的Cv(0.5，0.7，0.9)进行计算，得到的气溶胶光学厚度和真实值的对比如图3(b)。拟合系数(R2)提高到0.9818，均方根误差(RMSE)达到0.055，反演精度得到了提高，说明当植被混合比较高时，本文方法可以获得更好的精度。

同时，由于本文使用了一种一维搜索算法解方程组，当引入更多组数据时可以获得更好的结果，于是引入更多组数据进行分析，计算方式如式(9):

图3 3组Cv反演结果与真实值对比图Fig.3 Comparison between the retrieved values and the observed values with 3 groups of Cv

当引入4组Cv(0.5，0.6，0.7，0.9)，5组Cv(0.5，0.6，0.7，0.8、0.9)时，同样比较模拟值和真实值，得到4图。

由图4可知，当引入4～5组Cv进行计算时，拟合系数和均方根误差的变化并不大，在选取4组Cv时均方根误差甚至出现了增加，说明3组反演数据已经可以反演出较好的结果，考虑到增加了反演像元的个数和计算量，本文仍然选取3个像元作为反演像元。

3.2 真实遥感影像分析

为了进一步检验方法的有效性，本文选取了一景遥感影像进行验证。本文选取了2007年6月24号(Path 8/Row 61)中国陕西省的一景CBERS 0 2 BCCD影像(如图5(a))。校正CBERS影像数据，通过分析反演得到的地表反射率来验证本算法校正结果。使用本算法对此景影像进行计算，得到气溶胶光学厚度分布图如图5(b)。

图4 4组和5组Cv反演结果与真实值对比图Fig.4 Comparison of the retrieved values and the observed values with 4 groups of Cv and 5 groups of Cv

图5 气溶胶光学厚度分布图Fig.5 The original image and the VISmap made by this algorithm

如图5所示，根据原始影像显示，影像东边中部植被显示较暗;西南边建筑部分也显得模糊，可见能见度较低;西边中央与西南角植被地区红色鲜艳清晰，能见度较高。由气溶胶光学厚度分布图可知，这些肉眼可见的信息均显现了出来。由此可见，使用本文方法反演的CBERS影像的能见度结果能够比较合理地反映出气溶胶光学厚度的空间变化。

同时，使用本文方法，对影像进行大气校正，得到影像地表反射率，与校正前的表观反射率和真实地物的反射率进行对比，观察大气校正的效果。

如图6所示，通过大气校正前后的反射率对比，大气校正前图像整体偏暗，对比度不高，经过校正后，影像对比度有所提高，植被部分红色更加鲜艳，校正后模糊部分有所减弱，地物显得更加清晰可辨。出现这种差别的原因一方面在于大气的吸收和散射消弱了地表反射能量，消减了影像信号，导致气溶胶光学厚度较高的地方影像出现模糊不清;另一方面大气在近红外波段对植被能量的影响，导致校正前植被显得较为暗淡，校正后植被在近红外波段能量明显提高，使之在校正后显得更加鲜艳。

图6 校正前后效果图(左边为校正前效果，右边为校正后效果)Fig.6 Images before and after atmospheric correction

除了通过视觉比较校正前后的差异，反射率光谱特征也是分辨校正效果的重要手段，本文在研究区选取特征水体、植被像元，对其光谱信息进行平均，比较校正前后表观反射率和地表反射率差异(如图 7)。

图7 光谱信息及水体对比图Fig.7 Comparison of the spectral information of vegetation(left)and water(right)before and after atmospheric correction

如图7所示，植被反射率光谱信息在未经大气校正前，其表观反射率在蓝光、绿光和红光波段反射率较高且数值相差不大，没有明显的绿峰;近红外波段与可见光数值差距较小。经过大气校正后。蓝光波段的反射率值降低，近红外升高，与地物的真实地表反射率更为接近。水体反射率光谱在经过大气校正后在可见光波段反射率都有不同程度的降低，更加接近水体的真实测量值。

4 结论

本文构建了一种精确的针对不同植被覆盖度条件下的气溶胶光学厚度反演方法。通过使用可见光到近红外波段(400～900nm)影像信息中不同植被混合度像元的关系构建气溶胶光学厚度的反演关系式，最终使用一种一维迭代的方法解得关系式，得到反演的气溶胶光学厚度分布图。

方法的验证分为2个部分，首先使用模拟数据进行反演，反演结果和真实值的拟合系数(R2)达到0.981 8，均方根误差(RMSE)达到0.055。说明本方法在能找到合适的像元的情况下可以得到较好的反演结果，第2部分对CBERS影像进行气溶胶光学厚度的反演和大气校正，通过视觉比较和光谱信息分析，本文方法在视觉上恢复了地物的明暗差别，增强了影像对比度，同时可以反映出影像能见度的情况分布，经过大气校正后，可以使得地物更加接近真实地表反射率，有利于影像的判读与分析。

［1］KAUFMAN Y J，TANRE D，BOUCHER O.A Satellite View of Aerosols in the Climate System［J］.Nature，2002，419(6903):215－223.

［2］亓雪勇，田庆久.光学遥感大气校正研究进展［J］.国土资源遥感，2005，(4):1－6.(QI Xue-yong，TIAN Qing-jiu.The Advances in the Study of Atmospheric Correction for Optical Remote Sensing［J］.Remote Sensing for Land ＆ Resources，2005，(4):1－6.(in Chinese))

［3］LIANG S L.Recent Developments in Estimating Land Surface Biogeophysical Variables from Optical Remote Sensing［J］.Progress in Physical Geography，2007，31(5):501－516.

［4］RICHTER R.A Spatially Adaptive Fast Atmospheric Correction Algorithm［J］.International Journal of Remote Sensing，1996，17(6):1202－1211.

［5］RICHTER R.Atmospheric Correction of Satellite Data with Haze Removal Including a Haze/Clear Transition Region［J］.Computers ＆ Geosciences，1996，22(6):675－681.

［6］LIANG SL，FANG H L，CHEN M Z.Atmospheric Cor-rection of Landsat ETM+Land Surface Imagery-Part I:Methods［J］.IEEE Transactions on Geoscience and Remote Sensing，2001，39(11):2490－2498.

［7］LIANG S L，FANG H L，MORISETTE J T，etal.Atmospheric Correction of Landsat ETM Plus Land Surface Imagery-Part II:Validation and Applications［J］.IEEE Transactions on Geoscience and Remote Sensing，2002，40(12):2736－2746.

［8］SCHOTT JR，SALVAGGIO C，VOLCHOK W J.Radiometric Scene Normalization Using Pseudoinvariant Features［J］.Remote Sensing of Environment，1988，26(1):1－14.

［9］HALL F G，BOTKIN D B，STREBEL D E，etal.Large-Scale Patterns of Forest Succession as Determined by Remote Sensing［J］.Ecology，1991，72(2):628－640.

［10］梁顺林.定量遥感［M］.北京:科学出版社，2009.(LIANG Shun-lin.Quantitative Remote Sensing［M］.Beijing:Science Press，2009.(in Chinese))

［11］VERMOTE E F，El SALEOUSN，JUSTICE C O，etal.Atmospheric Correction of Visible to Middle-Infrared EOS-MODISData over Land Surfaces:Background，Operational Algorithm and Validation［J］.Journal of Geophysical Research，1997，102(D14):17131－17141.

［12］VERMOTE E F，El SALEOUSN，JUSTICE CO.Atmospheric Correction of MODIS Data in the Visible to Middle-Infrared:First Results［J］.Remote Sensing of Environment，2002，83(1/2):97－111.

［13］KAUFMAN Y J，WALD A E，REMER L A，etal.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.

［14］RICHTER R，SCHLAPFER D，MULLER A.An Automatic Atmospheric Correction Algorithm for Visible/NIR Imagery［J］.International Journal of Remote Sensing，2006，27(9/10):2077－2085.

［15］GUANTER L，GONZALEZ-SANPEDRO M D，MORENO J.A Method for the Atmospheric Correction of ENVISAT/MERIS Data over Land Targets［J］.International Journal of Remote Sensing，2007，28(3/4):709－728.

［16］GUTMAN G，IGNATOV A.The Derivation of the Green Vegetation Fraction from NOAA/AVHRR Data for Use in Numerical Weather Prediction Models［J］.International Journal of Remote Sensing，1998，19(8):1533－1543.

［17］KAUFMAN Y J，TANRE D，HOLBEN B N，etal.Atmospheric Effects on the NDVI——Strategies for Its Removal［C］∥Proceedings of the Geoscience and Remote Sensing Symposium，IGARSS，1992，2:1238－1241.

［18］KAUFMAN Y J，TANRE D.Atmospherically Resistant Vegetation Index(ARVI)for EOS MODIS［J］.IEEE Transactions on Geoscience and Remote Sensing，1992，30(2):261－270.

［19］何立明.气溶胶光学厚度与水平气象视距相互转换的经验公式及其应用［J］.遥感学报，2003，7(5):373－378.(HE Li-ming.Analysis and Application for the Empirical Relative between Aerosol Optical Depth and Horizontal Meteorological Range［J］.Journal of Remote Sensing，2003，7(5):373－378.(in Chinese))