張 鑫，趙 祥，d，劉素紅

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

1 研究背景

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

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

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

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

2 反演方法

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

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

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

2.1 篩選反演像元

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

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

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

對(duì)比使用NDVI和ARVI計(jì)算在不同能見度下計(jì)算的Cv，使用模擬數(shù)據(jù)進(jìn)行計(jì)算Cv為0～1時(shí)計(jì)算得到的植被混合比和真實(shí)值的誤差和，得到表1。

表1 不同能見度下分別使用NDVI和ARVI計(jì)算不同Cv的均方根誤差(RMSE)Table 1 The RMSerrors of Cv on different VIScalculated through NDVI and ARVI

由表可知，在能見度VIS＜23km時(shí)，NDVI由于受到大氣影響，計(jì)算出的植被混合比誤差均大于ARVI，能見度越小誤差越大。通過使用ARVI計(jì)算植被混合比精度得到了提高。

2.2 反演氣溶膠光學(xué)厚度

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

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

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

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

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

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

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

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

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

3 結(jié)果與分析

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

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

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

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

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

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

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

圖3 3組Cv反演結(jié)果與真實(shí)值對(duì)比圖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)時(shí)，同樣比較模擬值和真實(shí)值，得到4圖。

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

3.2 真實(shí)遙感影像分析

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

圖4 4組和5組Cv反演結(jié)果與真實(shí)值對(duì)比圖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é)厚度的空間變化。

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

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

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

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

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

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

4 結(jié)論

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

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

［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ìn)展［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àn)公式及其應(yīng)用［J］.遙感學(xué)報(bào)，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))