一種不同植被覆蓋條件下的氣溶膠光學(xué)厚度反演方法

張 鑫，趙 祥，d，劉素紅

(北京師范大學(xué) a.地理學(xué)與遙感科學(xué)學(xué)院;b.遙感科學(xué)國家重點實驗室;c.環(huán)境遙感與數(shù)字城市北京市重點實驗室;d.北京師范大學(xué)全球變化與地球系統(tǒng)科學(xué)研究院，北京 100875)

1 研究背景

由于大氣的存在，太陽輻射經(jīng)過氣體分子的吸收和氣溶膠粒子的散射，得到減弱，同時部分散射信號直接或間接經(jīng)過地物反射進入到傳感器，又得到增強。反映在實際處理中，大氣影響降低了圖像的反差比，使圖像可讀性降低，增加了解譯的困難。氣溶膠光學(xué)厚度(AOD)是形容大氣效應(yīng)中氣溶膠散射對大氣影響的參數(shù)，是氣溶膠粒子在整層大氣中消光系數(shù)的積分。也是大氣校正中最為重要的參數(shù)［1］。如何反演氣溶膠光學(xué)厚度，對影像進行大氣校正是至關(guān)重要的［2－3］。

到目前為止，遙感影像的大氣校正方法很多。如Richter等人提出的直方圖匹配法［4－5］，這種算法假定清晰區(qū)域和模糊區(qū)域的地表反射率的直方圖是相同的，在一景影像中辨認(rèn)出清晰和模糊的區(qū)域，匹配其直方圖得到氣溶膠光學(xué)厚度。Liang等人提出了聚類匹配法［6－7］，這種做法類似于直方圖匹配法，假定每一種地物在不同大氣條件下的平均反射率是相同的，對地物的劃分更加精細(xì)。還有基于統(tǒng)計的不變目標(biāo)法［8－9］，這種做法假定圖像上存在一些反射率穩(wěn)定的像元，這些像元在不同時相的影像中存在線性關(guān)系，當(dāng)確定這個線性關(guān)系的時候，即可估算出氣溶膠光學(xué)厚度。Kaufman等人提出的暗目標(biāo)法是一種常用的大氣糾正方法［10］。該方法的基本原理是假定需要校正的圖像中存在黑暗像元區(qū)域(一般為濃密植被或者干凈的水體)，假定其地表朗伯面反射，大氣性質(zhì)均一。在忽略大氣多次散射輻照作用和鄰近像元漫反射作用的前提下，反射率很小的黑暗像元受到大氣的影響，使得其反射率相對增加，這個增加的部分可以認(rèn)為是大氣所產(chǎn)生的影響。當(dāng)前大多傳感器(MODIS［11］，AVHRR，ETM+等)使用暗目標(biāo)法進行大氣校正，這種校正方法簡單、易行、自動化程度高，且精度較高。然而，暗目標(biāo)法需要能夠選取合適的暗目標(biāo)像元，在校正中需要引入短波紅外的數(shù)據(jù)(SWIR)來確定暗目標(biāo)的位置［12－13］。

當(dāng)前存在很多對地觀測衛(wèi)星，如SPOT，CBERS01/02 CCD，HJ－A/B CCD等數(shù)據(jù)，這些數(shù)據(jù)往往只提供可見光至近紅外的光譜信息，這給暗目標(biāo)法的大氣校正帶來了困難［14］，與此同時需要選取濃密植被的必要條件也在一定程度上限制了暗目標(biāo)方法的應(yīng)用。Guanter［15］等人提出一種針對ENVISAT/MERIS的大氣校正方法，這種方法利用植被的光譜特性，針對不同植被覆蓋度的像元在各個波段間的關(guān)系使用最優(yōu)化算法耦合出氣溶膠光學(xué)厚度，這給大氣校正提出了新的思路。本方法同樣基于不同植被覆蓋條件下的遙感影像，通過建立方程式的方式求得氣溶膠光學(xué)厚度。

2003年10月21日中巴地球資源衛(wèi)星02星(CBERS02)成功發(fā)射，其CCD相機提供450nm至730nm的可見/近紅外波段的影像，本文方法基于可見光至近紅外波段的遙感信息，針對CBERS02星CCD影像進行大氣校正，驗證其校正結(jié)果并進行分析。

2 反演方法

根據(jù)空氣的污染程度，可以設(shè)定在一定范圍內(nèi)氣溶膠光學(xué)厚度變化相對保持不變，如在干凈的影像中，可以設(shè)定10km×10km范圍內(nèi)保持不變，如果影像受氣溶膠影像比較嚴(yán)重，可以設(shè)定1km×1km范圍內(nèi)保持不變。同時假定在這個范圍內(nèi)，植被與裸土類型是相同的。根據(jù)這個假設(shè)，在反演的過程中對遙感影像劃分網(wǎng)格，在每個網(wǎng)格內(nèi)，每個像元具有相同的氣溶膠光學(xué)厚度，且純凈植被和裸土具有相同的光譜信息。對每個網(wǎng)格分別進行反演，得到光學(xué)厚度分布圖。具體的操作流程如圖1。

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

如圖1所示，首先，對遙感影像進行定標(biāo)處理，得到每個像元的表觀輻亮度。第二部對影像進行網(wǎng)格劃分。在每個網(wǎng)格內(nèi)，選取反演所需的像元。

2.1 篩選反演像元

本方法通過擬合不同植被混合比的像元間的關(guān)系來反演氣溶膠光學(xué)厚度，因此，為了得到更好的結(jié)果，在每個網(wǎng)格中選取植被混合比(Cv)相差最大的3個像元作為反演的像元。本文假定像元中植被覆蓋部分的植被類型一致且密度相同，使用等密度模型(dense vegetation model)來計算，如式(1)。

式中NDVI為歸一化植被指數(shù)，其中NDVIsoil和NDVIveg分別為對應(yīng)于裸土和高垂直密度植被的植被指數(shù)。Gutman文中使用NDVI來計算混合像元的植被混合比，然而大氣中的水汽、臭氧、氣溶膠、銳利散射等增加或減少紅光和近紅外波段的反射，這影響了NDVI在受到大氣干擾條件下的精度，于是本文引入了抗大氣植被指數(shù)(ARVI)來替代NDVI計算植被混合比。如式(2):

ARVI使用藍光，紅光和近紅外3個波段的通道組合進行計算，可以在一定程度上消除氣溶膠對植被指數(shù)的干擾［17－18］，ARVI的計算公式如下:

對比使用NDVI和ARVI計算在不同能見度下計算的Cv，使用模擬數(shù)據(jù)進行計算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 反演氣溶膠光學(xué)厚度

根據(jù)計算網(wǎng)格內(nèi)所有像元的植被混合比，篩選出Cv＞0.5的像元。在其中選出3個植被混合比相差最大的像元(最大、最小、和中間值)作為反演氣溶膠指數(shù)的像元。本方法認(rèn)定其地表反射率是植被混合比的線性加成。如式(4):

式中ρveg和ρsoil為當(dāng)前網(wǎng)格內(nèi)植被像元與裸土的地表反射率。根據(jù)式(4)得到以下反演關(guān)系式:

式中ρ1，ρ2，ρ3分別為這3 個像元的地表反射率，由輻射傳輸方程得到，在假定地表朗伯的情況下，有式(6)［15］:

式中Lp，EdTv，S使用MODTRAN計算得到。為了簡化計算過程，我們建立一個反射率計算查找表(LUTS)。查找表各個參數(shù)的設(shè)置如下:20個大氣能見度(2～60km);9個太陽天頂角(10°～78°)。根據(jù)查找表可以很方便和快速地根據(jù)給定的氣溶膠光學(xué)厚度計算出像元的地表反射率。其中氣溶膠光學(xué)厚度和大氣能見度是可以根據(jù)線性公式轉(zhuǎn)換得到的。有以下MODTRAN模擬的經(jīng)驗公式［19］:

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

表2 氣溶膠光學(xué)厚度和大氣能見度轉(zhuǎn)換系數(shù)Table 2 The conversion coefficients between AOD and VIS

由于地表反射率是由查找表計算得到的，我們不能簡單的以解三元一次方程的方式對方程求解，于是使用一個一維迭代尋求最優(yōu)值的方式，遍歷所有AOD范圍內(nèi)的值，每引入一個AOD，根據(jù)式(5)中前2個方程得到ρveg和ρsoil，代入方程3中，求出方程左右兩邊之差δ:

當(dāng)δ最小時，求得方程的最優(yōu)解AOD。

最后，得到每個網(wǎng)格的氣溶膠光學(xué)厚度分布圖，使用一個3×3的空間濾波對結(jié)果進行平滑，得到最終的氣溶膠光學(xué)厚度圖。

3 結(jié)果與分析

本文分為2個部分對本文的算法結(jié)果進行分析:前半部分使用模擬數(shù)據(jù)來驗證算法的精確性，后半部分使用本文算法對CBERS衛(wèi)星數(shù)據(jù)的影像進行大氣校正，檢驗其校正結(jié)果。

3.1 模擬數(shù)據(jù)分析

第一部分選擇模擬數(shù)據(jù)對算法進行驗證，我們獲取植被和裸土的真實地表反射率(如圖2)，由式(4)計算得到不同植被混合比下的地表反射率。使用查找表逆向計算這些地表反射率在不同氣溶膠光學(xué)厚度下的表觀輻亮度，通過得到的表觀輻亮度使用本文的算法計算出氣溶膠光學(xué)厚度。最后與真實的結(jié)果進行對比，以此來效驗算法的效果。

圖2 純凈植被和裸土的地表反射率光譜信息Fig.2 The spectral information of surface reflectance on pure vegetation and bare soil

首先選擇3組Cv(0.2，0.5，0.8)的表觀輻亮度進行計算，計算出的光學(xué)厚度與真實值進行比較，得到圖3(a)。

如圖3(a)所示，反演結(jié)果與真實值的均方根誤差(RMSE)達到了0.14，擬合系數(shù)(R2)達到0.94，誤差較大。為了提高反演的精度，我們另外使用了3組大于0.5的Cv(0.5，0.7，0.9)進行計算，得到的氣溶膠光學(xué)厚度和真實值的對比如圖3(b)。擬合系數(shù)(R2)提高到0.9818，均方根誤差(RMSE)達到0.055，反演精度得到了提高，說明當(dāng)植被混合比較高時，本文方法可以獲得更好的精度。

同時，由于本文使用了一種一維搜索算法解方程組，當(dāng)引入更多組數(shù)據(jù)時可以獲得更好的結(jié)果，于是引入更多組數(shù)據(jù)進行分析，計算方式如式(9):

圖3 3組Cv反演結(jié)果與真實值對比圖Fig.3 Comparison between the retrieved values and the observed values with 3 groups of Cv

當(dāng)引入4組Cv(0.5，0.6，0.7，0.9)，5組Cv(0.5，0.6，0.7，0.8、0.9)時，同樣比較模擬值和真實值，得到4圖。

由圖4可知，當(dāng)引入4～5組Cv進行計算時，擬合系數(shù)和均方根誤差的變化并不大，在選取4組Cv時均方根誤差甚至出現(xiàn)了增加，說明3組反演數(shù)據(jù)已經(jīng)可以反演出較好的結(jié)果，考慮到增加了反演像元的個數(shù)和計算量，本文仍然選取3個像元作為反演像元。

3.2 真實遙感影像分析

為了進一步檢驗方法的有效性，本文選取了一景遙感影像進行驗證。本文選取了2007年6月24號(Path 8/Row 61)中國陜西省的一景CBERS 0 2 BCCD影像(如圖5(a))。校正CBERS影像數(shù)據(jù)，通過分析反演得到的地表反射率來驗證本算法校正結(jié)果。使用本算法對此景影像進行計算，得到氣溶膠光學(xué)厚度分布圖如圖5(b)。

圖4 4組和5組Cv反演結(jié)果與真實值對比圖Fig.4 Comparison of the retrieved values and the observed values with 4 groups of Cv and 5 groups of Cv

圖5 氣溶膠光學(xué)厚度分布圖Fig.5 The original image and the VISmap made by this algorithm

如圖5所示，根據(jù)原始影像顯示，影像東邊中部植被顯示較暗;西南邊建筑部分也顯得模糊，可見能見度較低;西邊中央與西南角植被地區(qū)紅色鮮艷清晰，能見度較高。由氣溶膠光學(xué)厚度分布圖可知，這些肉眼可見的信息均顯現(xiàn)了出來。由此可見，使用本文方法反演的CBERS影像的能見度結(jié)果能夠比較合理地反映出氣溶膠光學(xué)厚度的空間變化。

同時，使用本文方法，對影像進行大氣校正，得到影像地表反射率，與校正前的表觀反射率和真實地物的反射率進行對比，觀察大氣校正的效果。

如圖6所示，通過大氣校正前后的反射率對比，大氣校正前圖像整體偏暗，對比度不高，經(jīng)過校正后，影像對比度有所提高，植被部分紅色更加鮮艷，校正后模糊部分有所減弱，地物顯得更加清晰可辨。出現(xiàn)這種差別的原因一方面在于大氣的吸收和散射消弱了地表反射能量，消減了影像信號，導(dǎo)致氣溶膠光學(xué)厚度較高的地方影像出現(xiàn)模糊不清;另一方面大氣在近紅外波段對植被能量的影響，導(dǎo)致校正前植被顯得較為暗淡，校正后植被在近紅外波段能量明顯提高，使之在校正后顯得更加鮮艷。

圖6 校正前后效果圖(左邊為校正前效果，右邊為校正后效果)Fig.6 Images before and after atmospheric correction

除了通過視覺比較校正前后的差異，反射率光譜特征也是分辨校正效果的重要手段，本文在研究區(qū)選取特征水體、植被像元，對其光譜信息進行平均，比較校正前后表觀反射率和地表反射率差異(如圖 7)。

圖7 光譜信息及水體對比圖Fig.7 Comparison of the spectral information of vegetation(left)and water(right)before and after atmospheric correction

如圖7所示，植被反射率光譜信息在未經(jīng)大氣校正前，其表觀反射率在藍光、綠光和紅光波段反射率較高且數(shù)值相差不大，沒有明顯的綠峰;近紅外波段與可見光數(shù)值差距較小。經(jīng)過大氣校正后。藍光波段的反射率值降低，近紅外升高，與地物的真實地表反射率更為接近。水體反射率光譜在經(jīng)過大氣校正后在可見光波段反射率都有不同程度的降低，更加接近水體的真實測量值。

4 結(jié)論

本文構(gòu)建了一種精確的針對不同植被覆蓋度條件下的氣溶膠光學(xué)厚度反演方法。通過使用可見光到近紅外波段(400～900nm)影像信息中不同植被混合度像元的關(guān)系構(gòu)建氣溶膠光學(xué)厚度的反演關(guān)系式，最終使用一種一維迭代的方法解得關(guān)系式，得到反演的氣溶膠光學(xué)厚度分布圖。

方法的驗證分為2個部分，首先使用模擬數(shù)據(jù)進行反演，反演結(jié)果和真實值的擬合系數(shù)(R2)達到0.981 8，均方根誤差(RMSE)達到0.055。說明本方法在能找到合適的像元的情況下可以得到較好的反演結(jié)果，第2部分對CBERS影像進行氣溶膠光學(xué)厚度的反演和大氣校正，通過視覺比較和光譜信息分析，本文方法在視覺上恢復(fù)了地物的明暗差別，增強了影像對比度，同時可以反映出影像能見度的情況分布，經(jīng)過大氣校正后，可以使得地物更加接近真實地表反射率，有利于影像的判讀與分析。

［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］亓雪勇，田慶久.光學(xué)遙感大氣校正研究進展［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，F(xiàn)ANG 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，F(xiàn)ANG 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］.北京:科學(xué)出版社，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］何立明.氣溶膠光學(xué)厚度與水平氣象視距相互轉(zhuǎn)換的經(jīng)驗公式及其應(yīng)用［J］.遙感學(xué)報，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))