利用Landsat-8 OLI反演大氣氣溶膠的可見光譜段地表反射率估算

呂春光1，2，田慶久1，2，王磊1，2，黃彥1，2，耿君1，2

(1.南京大學 國際地球系統(tǒng)科學研究所，南京 210023；2.江蘇省地理信息技術重點實驗室，南京 210023)

摘要：準確估算地表反射率的貢獻一直是遙感反演大氣氣溶膠光學厚度過程中的重點和難點。為了促進Landsat-8 OLI傳感器在地表參數定量化特別是大氣遙感領域的應用，本文提出一種利用OLI 1.6μm、2.2μm短波紅外譜段數據估算遙感影像可見光地表反射率的方法。該方法依托于MOD04產品地表反射率估算模式，通過光譜歸一化和構建新的短波紅外植被指數等過程，建立OLI地表反射率估算模式，通過誤差分析發(fā)現該模式能夠有效地降低由于傳感器光譜響應不同對估算結果的影響，對應用在OLI遙感影像的計算結果與同時間同區(qū)域MOD04產品地表反射率進行比較，表明其結果有較高的相關性和可靠性。

關鍵詞：地表反射率；氣溶膠；Landsat-8；光譜歸一化；植被指數

doi:10.3969/j.issn.1000-3177.2015.01.008

中圖分類號：TP79文獻標識碼：A

Estimation of Visible Surface Reflectance for Retrieving Aerosol

Optical Depth Using Landsat-8 OLI Data

LV Chun-guang1，2，TIAN Qing-jiu1，2，WANG Lei1，2，HUANG Yan1，2，GENG Jun1，2

(1.InternationalInstituteforEarthSystemScience，NanjingUniversity，Nanjing210023；

2.JiangsuProvincialKeyLaboratoryofGeographicInformationScienceand

Technology，NanjingUniversity，Nanjing210023)

Abstract：Accurate estimation to the contribution of the surface reflectance has been the focus and the difficult point in aerosol optical depth (AOD) retrieval from remote sensing data.In order to promote the applications of Landsat-8 OLI data in quantification of surface physical parameters particularly in the field of atmospheric remote sensing，a method to estimate the optical surface reflectance from OLI images was put forward using OLI 1.6μm and OLI 2.2μm short infrared bands.The method was proposed based on the surface reflectance estimation mode of MOD04 product，by means of processes such as spectral normalization，building a new short infrared vegetation index and so on.The error analysis shows the OLI estimation mode can effectively reduce the error effect caused by spectral response differences.The comparison of the surface reflectance estimation results between the OLI and MOD04 product in the same time and same area represents that the estimation results obtained by OLI estimation mode have high correlation and reliability.

Key words：surface reflectance；aerosol；Landsat-8；spectral normalization；vegetation index

1引言

大氣氣溶膠是大氣中重要的成分之一，它會影響大氣輻射狀況，對全球和區(qū)域氣候變化、大氣環(huán)境質量具有重要的作用和影響，是當前大氣遙感研究的重要領域之一[1]。遙感反演大氣氣溶膠光學厚度(Aerosol Optical Depth，AOD)具有重要的科學和現實意義[2]。有研究指出0.01的地表反射率估算誤差將會導致0.1的AOD反演誤差，同時將會顯著影響氣溶膠尺度等物理參數的估算[3-4]，因此，準確估算地表反射率的貢獻是AOD反演中的關鍵環(huán)節(jié)。

暗目標法是目前最常用的估算方法，它利用濃密植被區(qū)紅光和藍光波段地表反射率比較低，且與短波紅外地表反射率存在固定的線性關系，從而去除反演中地表對表觀反射率的貢獻[5-6]。更多的研究表明，在不同觀測條件下，可見光波段與短波紅外波段地表反射率的比值在一定的范圍內波動[7]。中分辨率成像光譜儀(Moderate Resolution Imaging Spectroradiometer，MODIS)C005算法是目前氣溶膠反演廣泛使用的地表反射率估算方法，它考慮了上述植被覆蓋情況和散射角對地表反射率估算的影響[8]，在植被覆蓋好的地區(qū)具有較高的精度[6]。

但是，由于MODIS氣溶膠產品(MOD 04)僅有10km的空間分辨率，從而限制了它的應用范圍[9]。2013年2月發(fā)射的Landsat-8衛(wèi)星上搭載的陸地成像儀(Operational Land Imager，OLI)具有30 m空間分辨率，覆蓋可見光到短波紅外的譜段設置，因此具備利用類似MODIS C005算法的可能性[10]。本文以暗目標法的MODIS地表反射率估算模式為依托，以實測地表反射率數據為基礎，考慮傳感器的光譜響應關系和譜段設置，嘗試建立針對OLI傳感器的可見光地表反射率估算模式，為OLI可見光地表反射率的應用提供支持。

2原理與方法

2.1原理

當遙感影像像元(水體像元除外)在2.1μm處的表觀反射率滿足0.01<ρm2.1<0.25時，該像元可以作為暗目標，在2.1μm處的地表反射率和0.47μm、0.66μm處的地表反射率之間的關系與散射角和植被的茂密程度有關，MODIS地表反射率標準估算模式如下[6-7]：

(1)

(2)

(3)

yint0.66/2.1=-0.00025Θ+0.033

(4)

(5)

(6)

式中，ρs0.66，ρs0.47，ρs2.1分別表示0.66μm，0.47μm，2.1μm波段的地表反射率；ρm1.24，ρm2.1分別表示1.24，2.1μm波段的表觀反射率；Θ表示散射角。

圖1　OLI和MODIS的不同波段的光譜響應函數與 不同濃密程度植被的光譜對應關系

OLI傳感器分別在藍光波段(0.48μm)、紅光波段(0.66μm)和短波紅外(2.1μm)3個波段范圍與MODIS傳感器存在對應關系，但是除紅光波段外，每個波段的響應明顯要寬于MODIS(圖1)。OLI傳感器在1.24μm處并沒有相應的波段設置，因而不能借助NDVISWIR來建立短波紅外與紅光反射率之間的關系，但OLI在1.6μm附近的波段位于大氣窗口，可考慮與短波紅外(2.12μm)結合建立地表反射率的關系。式(4)和式(5)表明散射角Θ在90°～180°變化將會導致MODIS短波紅外與紅光波段的斜率產生±0.09的差異，同時對截距產生-0.012～0.0105的差異。地表反射率隨散射角的變化關系也會隨傳感器的光譜響應不同而產生差異。

2.2實驗與方法

本研究于2013年3月～4月在安徽江蘇等地區(qū)采集了198組植被冠層光譜數據。實驗采用ASD野外光譜儀，其波長范圍為350nm～2500nm，其中350nm～1000nm的平均光譜分辨率為3nm，1000nm～2500nm的平均光譜分辨率為10nm，于北京時間11∶00～14∶00之間進行測量，目標光譜測定前、后均立即進行參考板校正，在每個測點采用5次平均值。為了得到OLI傳感器地表反射率，并對比與MODIS之間的差別，利用OLI和MODIS的光譜響應函數對實測光譜進行重采樣，分別得到MODIS和OLI傳感器短波紅外2.1μm、2.2μm、1.62μm、1.24μm(僅MODIS)、紅光波段和藍光波段的地表光譜反射率，光譜采樣定理如式(7)所示[11]：

(7)

式中，ρs(sensor)為光譜重采樣得到的傳感器地表反射率，λ為波長，ρs為實測地表反射率，f為傳感器光譜響應函數，λmin和λmax為波段起止波長。

3結果與分析

3.1建立OLI地表反射率估算模式

如圖2(a)，通過對短波紅外(2.1μm、2.2μm)、紅光波段和藍光波段MODIS與OLI地表反射率之間關系可以發(fā)現各波段之間分別存在以下關系：

(8)

式中，ρs(MODIS)2.1、ρs(MODIS)0.66、ρs(MODIS)0.47分別為模擬的MODIS地表反射率，ρs(OLI)2.2、ρs(OLI)0.65、ρs(OLI)0.48分別為模擬的OLI地表反射率。如圖2(a)所示，其在這兩個傳感器上有顯著的相關性(R2均高于0.99，如圖2(a))。將式(8)第2、3式帶入式(2)可得到OLI紅光波段和藍光波段的關系：

(9)

圖2

(10)

(11)

(12)

利用OLIslopeNDVI′SWIR0.65/2.2替換對應MODISslopeNDVISWIR0.66/2.12，對與每組得到的ρs(MODIS)0.66光譜歸一化后的ρs′(OIL)0.65進行比較，使得偏差最小，如式(13)所示：

(13)

由式(13)確定系數后，整理得到以下關系式：

(14)

yint0.65/2.2=-0.000284Θ+0.0352

(15)

將上式(14)、(15)帶入式(1)中，并與式(9)、(10)、(11)聯(lián)立即可得到OLI傳感器紅光和藍光波段地表反射率估算模式，Θ為散射角。

針對實測數據光譜重采樣得到的OLI 2.2μm反射率，利用上述推導得到的OLI地表反射率模式和直接使用MODIS地表反射率模式對紅光和藍光地表反射進行估算，并分別將估算得到的OLI地表反射率結果與MODIS標準估算模式得到的光譜歸一化后的地表反射率結果進行比較，為了表現大小差異，利用下列誤差公式進行評價：

(16)

式中，r為地表反射率估算模式誤差，ρ′可為利用OLI短波紅外波段、OLI估算模式或直接利用MODIS估算模式得到的反射率，ρ為利用MODIS短波紅外波段和估算模式計算，并進行OLI光譜歸一化后的反射率。圖3中，曲線(a)、(b)為對OLI短波紅外波段直接利用MODIS標準估算模式估算的誤差曲線，其平均誤差分別為22%和8.3%；而曲線(c)、(d)為利用建立的OLI估算模式的誤差曲線，平均相對誤差分別為小于1.5%和1%。這表明考慮傳感器光譜響應差異建立OLI估算模式可有效減小MODIS標準估算模式帶來的模式誤差，采用建立的OLI估算模式估算的地表反射率略低于標準值，而利用MODIS標準模式將會顯著高估濃密植被區(qū)的地表反射率。曲線(a)、(b)在短波紅外小于0.12的范圍內，其模式誤差r顯著大于曲線(c)、(d)，由于濃密植被的2.1μm短波紅外反射率主要集中在0.12以下，因此如不考慮傳感器光譜響應的差異，直接使用MODIS標準模式的反射率估算關系，其模式本身將會對估計值帶來較大誤差。

圖3　對OLI短波紅外數據使用不同地表反射率估算模式 產生的模式誤差隨2.1μm反射率的變化比較 ((a)直接利用MODIS標準估算模式估算紅光波段； (b)直接利用MODIS標準估算模式估算藍光波段； (c)利用建立的OLI估算模式估算紅光波段； (d)利用建立的OLI估算模式估算藍光波段)

3.2OLI影像地表反射率估算應用與分析

本文選取5景2013年9月1日中國東部地區(qū)的Landsat-8 OLI影像數據來進行紅光和藍光波段地表反射率估算。在這5景影像的獲取時間范圍內天氣晴好，各景影像的獲取時間、星下點位置和圖幅范圍如表1所示。選取對應上述OLI影像空間范圍的Terra MODIS數據，其時間范圍與OLI影像相差在40分鐘以內，表1中各序號對應影像的幾何位置關系如圖4所示。

首先根據影像頭文件信息逐像元計算太陽天頂角、方位角、觀測天頂角、方位角和散射角，并求得短波紅外植被指數；其次，根據植被指數和影像2.1μm波段設定閾值掩膜水體、城鎮(zhèn)、云等非暗像元信息，并對2.1μm波段大氣校正得到其地表反射率；最后，利用本文建立的OLI估算模式求得紅光和藍光波段地表反射率，結果如圖5(a)、(b)所示。

MOD04數據是全球范圍使用廣泛的MODIS氣溶膠產品，它對地表反射率的估算采用了MODIS標準反射率模式[12]。為了對結果進行比較，本文對OLI估算結果進行像元空間重采樣，與2013年9月1日MOD04產品的像元進行對應，并從中篩選出非暗像元比例小于10%的樣點作為純像元進行比較。樣點OLI和MODIS紅光波段地表反射率均值分別為0.0437和0.0485，藍光波段地表反射率均值分別為0.0279和0.0285。根據像元的經緯度位置進行編號，其中MOD04與OLI地表反射率對應樣點的情況如圖5(c)所示。從圖中可以發(fā)現，OLI估算結果在紅光和藍光波段均低于對應點位的MOD04結果，對應曲線的波動趨勢具有明顯的一致性。對應樣點OLI 和MOD04產品短波紅外(2.2μm、2.1μm)、紅光波段(0.6μm)和藍光波段(0.4μm)各自的相關關系如圖5(d)所示，相關系數R2分別為0.57、0.563和0.538。從圖5(d)中散點的分布情況看，MOD04數據某一值通常對應多個分布在周圍的不同OLI反射率值，這表明OLI數據對地表反射率的差異更敏感。

表1　所選衛(wèi)星影像信息

圖4　2013年9月1日Terra MODIS 1-2-1波段合成與 Landsat-8 OLI假彩色合成衛(wèi)星影像位置關系

4討論

圖5(c)表現了在同類型的反射率估算模式下，OLI紅光和藍光地表反射率整體比MODIS低10%和2.3%，導致這種差異的因素很復雜，其中包括傳感器的光譜響應、空間尺度效應、地物的反射特性和傳感器的觀測角度等。由于OLI的觀測天頂角很小(最大幅寬小于7°)，本文對OLI地表反射率，利用式(8)進行光譜歸一化，發(fā)現處理后的紅光波段和藍光波段地表反射率分別提高了3.5%和0.01%，但仍然比MODIS地表反射率低6.6%和2.1%。

圖5　OLI影像地表反射率估算結果及其與 MOD04數據的比較

理論上，OLI植被指數的相對增高，會使式(1)短波紅外波段2.2μm和紅光波段0.65μm反射率之間的斜率相對升高從而導致OLI可見光地表反射率高于MODIS，但實際過程中，散射角的相對大小對斜率有調節(jié)作用，并對截距的影響更加顯著，從而導致MODIS地表反射率的相對抬升。OLI中所選樣點的平均散射角為147.3°，而MODIS對應樣點的平均散射角僅為93.7°，根據斜率和截距公式得到各自斜率分別為0.539和0.502，截距分別為-0.005和0.002，而2.1μm平均反射率為0.09，這使得OLI與MODIS相比，受散射角影響的截距抵消了由植被指數和散射角帶來的斜率的升高，并進一步將反射率整體降低了0.0037。由此可見，散射角的差異對OLI與MODIS的反射率差異中起決定作用。

以上從地表反射率的影響因素和等式響應機制上對估算結果的差異進行了討論。此前已有相關的研究表明在傳感器對植被觀測方向與太陽同向時，其紅光反射率有可能較垂直觀測時升高[14，18]，這說明MODIS在較大觀測角度觀測同一目標時，其紅光地表反射率有可能比OLI要高。此外，幾何定位精度也可能導致兩者地表反射率出現差異，MODIS MOD04數據產品和Landsat OLI影像數據本身具備較高的幾何校正精度，對二者地表反射率相關性影響較小，由于篇幅和內容所限本文不再對其進行過多的論述。

隨著研究的深入，利用紅外對可見光地表反射率的估算模式從簡單的比值模型發(fā)展為引入觀測角和植被指數等具有物理意義的復雜模型，也出現了更多不同思路和針對不同傳感器估算可見光地表反射率的模式[7，19-23]。不同傳感器來探測地表反射率不但受到觀測幾何、目標特性、波譜響應、瞬時視場和輻射靈敏度的影響，還與傳感器的在軌時間、老化程度和衛(wèi)星的運行狀態(tài)有關，研究和發(fā)展針對不同傳感器，并引入更多物理參量提高結果精度和抗干擾因素能力在未來仍具有重要意義和實用價值。

5結束語

本文針對Landsat-8 OLI數據的譜段設置，在MODIS暗目標法估算可見光譜段地表反射率模式的基礎上，通過光譜歸一化以及構建新的短波紅外植被指數等過程，提出了一種新的用于Landsat-8 OLI 2.2μm、1.6μm短波紅外譜段數據估算0.65μm紅光和0.48μm藍光譜段地表反射率的方法，并利用該方法與同一時期同類型估算模式的MOD04產品地表反射率進行了比較。本文主要結論有：

(2)利用實測光譜數據的誤差分析表明，對OLI譜段直接使用MODIS標準估算模式進行地表反射率計算會產生較大的估算誤差，可使紅光和藍光估算結果平均偏高22%和8.3%，而使用建立的OLI估算模式所得紅光和藍光估算誤差則分別小于1.5%和1%。這說明考慮不同傳感器的光譜響應差異并對估算模式進行改進可有效減小估算模式本身所導致的誤差。

(3)對同一地區(qū)時間接近的OLI影像地表反射率的估算結果與MODIS產品進行比較發(fā)現，結果之間存在較高的線性相關性，而傳感器觀測角度的不同使得OLI地表反射率略低于MODIS，同時兩個傳感器地表反射率的差異也是由傳感器光譜響應差異、空間尺度效應、觀測幾何差異、地表植被結構以及衛(wèi)星運行狀態(tài)等諸多因素共同作用的結果。

致謝：本文感謝USGS earth explorer網站提供Landsat-8影像數據以及NASA LAADS提供的MODIS數據產品，并對各位審稿專家表示感謝。

參考文獻：

[1]DUBOVIK O，SMIRNOV A，HOLBEN B N，et al.Accuracy assessments of aerosol optical properties retrieved from aerosol robotic network (AONET) sun and sky radiance measurements[J].Journal of Geophysical Research-Atmospheres，2000，105(D8)：9791-9806.

[2]TORRE O，BHARTIA P K，HERMAN J R，et al.A long-term record of aerosol optical depth from TOMS observations and comparison to AERONET measurements[J].Journal of the Atmospheric Sciences，2002，59(3)：398-413.

[3]KAUFMAN Y J，TANRE D，GORDON H R，et al.Passive remote sensing of tropospheric aerosol and atmospheric correction for the aerosol effect[J].Journal of Geophysical Research-Atmospheres，1997，102(D14)：16815-16830.

[4]李莘莘，陳良富，陶金花，等.城市與冬季北方亮目標地區(qū)氣溶膠光學厚度反演[J].中國科學：地球科學，2012，(8)：1253-1263.

[5]KAUFMAN Y J，TANRE D，REMER L A，et al.Operational remote sensing of tropospheric aerosol over land from EOS moderate resolution imaging spectroradiometer[J].Journal of Geophysical Research-Atmospheres，1997，102(D14)：17051-17067.

[6]周春艷，柳欽火，唐勇，等.MODIS氣溶膠C004、C005產品的對比分析及其在中國北方地區(qū)的適用性評價[J].遙感學報，2009，(5)：854-872.

[7]LEVY R C，REMER L A，MATTOO，et al.Second-generation operational algorithm：Retrieval of aerosol properties over land from inversion of moderate resolution imaging spectroradiometer spectral reflectance[J].Journal of Geophysical Research-Atmospheres，2007，112(D13)：1-21.

[8]REMER L A，TANRE D，KAUFMAN Y J，et al.Algorithm for remote sensing of tropospheric aerosol from MODIS：Collection 005[J].National Aeronautics and Space Administration，2006：1-87.

[9]王玲，田慶久，李姍姍.利用MODIS資料反演杭州市500米分辨率氣溶膠光學厚度[J].遙感信息，2010，25 (3)：50-54.

[10]MARKHAM B L，STOREY J C，IRONS J R.Landsat data continuity mission-now Landsat-8：Six months on orbit[M].Bellingham：SPIE-Int SOC Optical Engineering，2013.

[11]梁順林，范聞捷.定量遙感[M].北京：科學出版社，2009.

[12]REMER L A，KAUFMAN Y J，TANRE D，et al.The MODIS aerosol algorithm，products，and validation[J].Journal of the Atmospheric Sciences，2005，62(4)：947-973.

[13]覃文漢，項月琴.植被結構及太陽/觀測角度對NDVI的影響[J].環(huán)境遙感，1996，(4)：285-290.

[14]張雪紅，田慶久，沈潤平.冬小麥冠層光譜的方向性特征分析[J].光譜學與光譜分析，2010，(6)：1600-1605.

[15]ZHANG X，YAN G，LI Q，et al.Evaluating the fraction of vegetation cover based on NDVI spatial scale correction model[J].Intnational Journal of Remote Sensing，2006，27(23-24)：5359-5372.

[16]欒海軍，田慶久，余濤，等.基于分形理論的NDVI連續(xù)空間尺度轉換模型研究[J].光譜學與光譜分析，2013，(7)：1857-1862.

[17]ZHANG X H.Study on spatial heterogeneity and scale effect of eucalyptus forest based on high resolution remote sensing[D].Nanjing：Nanjing University，2012.

[18]許民，宜樹華，葉柏生，等.植被蓋度及太陽/觀測角度對疏勒河上游NDVI和SAVI值的影響[J].干旱區(qū)資源與環(huán)境，2012，(5)：101-107.

[19]TANR E D，DESCHAMPS P Y，DEVAUX C，et al.Estimation of saharan aerosol optical thickness from blurring effects in thematic mapper data[J].Journal of Geophysical Research：Atmospheres (1984-2012)，1988，93 (D12)：15955-15964.

[20]盛莉，黃敬峰，郭瑞芳，等.綜合利用環(huán)境星CCD和紅外數據反演大氣氣溶膠光學厚度[J].紅外與毫米波學報，2013，(3)：277-282.

[21]王中挺，厲青，王橋，等.利用深藍算法從HJ-1數據反演陸地氣溶膠[J].遙感學報，2012，(3)：596-610.

[22]HAU N C，TSAY S C，KING M D，et al.Deep blue retrievals of asian aerosol properties during ACE-Asia[J].IEEE Transactions on Geoscience and Remote Sensing，2006，44(111)：3180-3195.

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