沈文娟， 李明詩

(南京林業(yè)大學林學院，南京 210037)

0 引言

Landsat衛(wèi)星已積累了超過40 a的連續(xù)對地觀測檔案數(shù)據(jù)，這些數(shù)據(jù)在區(qū)域生態(tài)環(huán)境監(jiān)測方面的應用優(yōu)勢明顯[1]。已有研究將其用于重構(gòu)過去幾十年的森林干擾及恢復歷史，并評價現(xiàn)存森林經(jīng)營方法的有效性[2-3]。陸地衛(wèi)星生態(tài)系統(tǒng)干擾自適應處理系統(tǒng)(landsat ecosystem disturbance adaptive processing system，LEDAPS)是美國航空航天局(NASA)的一個項目，主要研究長時間跨度的稠密時間序列Landsat GeoCover影像數(shù)據(jù)集的預處理算法。該影像數(shù)據(jù)集獲取的真實地表反射率能夠用于多種對象的長時間尺度評價。鑒于大區(qū)域范圍土地覆蓋變化的復雜性(包括火、森林采伐、城市化及農(nóng)業(yè)領(lǐng)域變化等)，干擾模式及其與氣候、生物多樣性和生態(tài)系統(tǒng)健康的定量化關(guān)系，實現(xiàn)長時間序列遙感數(shù)據(jù)的預處理和分析是科學分析生態(tài)系統(tǒng)變化的前提條件[4-5]。

現(xiàn)有用于大氣糾正的6S模型[6]，對云層較厚區(qū)域的反演精度較差[7]，應用時往往假設(shè)整個研究區(qū)的大氣條件均一，太陽和傳感器的觀測角度一致，地表具有相同的海拔高度[8]，但是對反演區(qū)的氣候和氣象資料以及地面狀況的了解信息不足，很多影像獲取時刻的氣象數(shù)據(jù)將很難準確獲得，所以常使用一些標準模式來估算相關(guān)參數(shù)，特別是水汽、O3含量[9]。本研究以ETM+和TM數(shù)據(jù)為例，實現(xiàn)了LEDAPS算法可處理數(shù)據(jù)的格式統(tǒng)一，并驗證長時間序列地面反射率處理產(chǎn)品的精度。在國內(nèi)推廣使用LEDAPS算法形成長時間序列可用的Landsat科學產(chǎn)品，對于推動多學科深入挖掘Landsat數(shù)據(jù)的科學價值具有重要的示范作用。

1 數(shù)據(jù)預處理

1.1 LEDAPS算法的原理

1.1.1 數(shù)據(jù)選擇與定標

LEDAPS的定標方法首先使用日平均定標燈亮度歷史數(shù)據(jù)移除先前錯誤的定標，再將修改后的定標系數(shù)用于相應日期的影像[4]。主要基于太陽方位、日地距離、TM或ETM+波段以及太陽輻照度，定標影像為表觀(top-of-otmosphere,TOA)反射率。 Landsat TOA反射產(chǎn)品的輸出格式為HDF-EOS，LEDAPS網(wǎng)站提供數(shù)據(jù)的下載。

1.1.2 大氣校正原理

大氣校正旨在補償大氣成分的散射和吸收效應，獲取精確的地表反射率值。表觀反射率ρTOA與地表反射率ρs的關(guān)系為

(1)

式中：Tg為其后面括號內(nèi)所含氣體的傳輸值；TR+A是瑞利和氣溶膠傳輸值；ρR+A是瑞利和氣溶膠大氣固有的反射值；SR+A是瑞利和氣溶膠球面反照率。傳輸、固有反射值和球面反照度由6S輻射傳輸模型計算所得[10]； O3含量來源于搭載在Nimbus-7，Meteor-3以及地球探測平臺上的O3總量繪圖系統(tǒng)(total ozone mapping spectrometer，TOMS)； 水蒸氣來源于NOAA國家環(huán)境預報中心(national centers for environmental prediction，NCEP)重分析數(shù)據(jù)，分辨率為2.5°×2.5°(http: //rda.ucar.edu/datasets/ds090.0/)。數(shù)字地形和NCEP表面壓力數(shù)據(jù)用于調(diào)整瑞利散射為實地的狀況。

Kaufman等[11]充分利用影像中濃密植被(DDV)信息直接提取氣溶膠光學厚度(aerosol optical thicrnecs,AOT)，在暗目標之間使用樣條算法空間插值獲取每個波段的氣溶膠光學厚度值[10,12-13]。將插值得到的AOT和通過相關(guān)資料得到的O3、大氣壓以及水汽值用于6S傳輸模型算法，獲取地表反射率。在處理過程中，還使用了Landsat7 自動云量評估算法進行云掩模[14]。在Landsat7項目中該方法提供以景為單位的云覆蓋元數(shù)據(jù)具有優(yōu)勢[4]。

1.2 地表反射率及其驗證

Masek等[4]總結(jié)了LEDAPS產(chǎn)品質(zhì)量評價的3種方法： ①將LEDAPS得到的地面反射數(shù)據(jù)與Terra MODIS日反射產(chǎn)品作比較； ②Landsat 影像氣溶膠光學厚度值與AERONET(aerosol robotic network)AOT作比較； ③Landsat地面反射值與高分辨率影像或者機載輻射數(shù)據(jù)的替代值作比較。本研究嘗試獲取LEDAPS處理前后的Landsat影像上典型地物反射率，并與標準波譜作比較，以此來判斷處理結(jié)果的精度。

1.3 LEDAPS模型實現(xiàn)

LEDAPS預處理代碼作為自動化的預處理和分析工具，是Landsat 1T數(shù)據(jù)定標、TOA反射率轉(zhuǎn)換、云掩模和大氣校正預處理的單獨版本，并且該軟件適用的數(shù)據(jù)為Landsat5和Landsat7標準EROS L1T產(chǎn)品，結(jié)果儲存在HDF-EOS格式中。HDF-EOS數(shù)據(jù)圖層包括所有反射波段的圖層(不含熱紅外波段)、藍光波段的氣溶膠光學厚度圖層以及一個“QA”圖層(含有云的分離標記，缺失數(shù)據(jù)和來源于Landsat4的數(shù)據(jù))。

LEDAPS模型(ftp://hrsl.arsusda.gov/pub/fgao/ Ledaps/)(圖1)含4個具體操作命令，分別為lndpm(提取標準Landsat產(chǎn)品元數(shù)據(jù)信息，生成后續(xù)處理使用文件)、lndcal(進行TM和ETM+數(shù)據(jù)的定標，獲取大氣上界反射率并得到大氣上界亮溫)、lndcsm(使用自動云量評估算法得到云掩模)、lndsr(使用MODIS/6S輻射傳輸模型和黑暗密集植被算法所得氣溶膠進行Landsat 反射波段的大氣校正)。模型輸出文件包含lndcal.*.hdf，lndth.*.hdf，lndcsm.*.hdf及l(fā)ndsr.*.hdf，每個文件均帶有地理信息的ENVI頭文件。

圖1 LEDAPS模型運行界面

2 數(shù)據(jù)格式的統(tǒng)一

在國內(nèi)推廣LEDAPS方法，可用的數(shù)據(jù)總結(jié)為以下3部分，即Landsat CDR(climate data records)產(chǎn)品、標準數(shù)據(jù)的LEDAPS處理數(shù)據(jù)以及非標準數(shù)據(jù)轉(zhuǎn)換后的LEDAPS處理數(shù)據(jù)。綜合以上數(shù)據(jù)，能夠形成完整的長時間序列數(shù)據(jù)產(chǎn)品。

2.1 標準數(shù)據(jù)準備

2.1.1 LEDAPS處理后的地表反射率數(shù)據(jù)

Landsat地表反射率CDR數(shù)據(jù)由美國地質(zhì)調(diào)查局(USGS)地球資源觀測和科學(earth resources observation and science，EROS)中心科學處理部門(ESPA)處理所得。ESPA提供了地表反射率數(shù)據(jù)以及相關(guān)的源數(shù)據(jù)和生成數(shù)據(jù)。數(shù)據(jù)由LEDAPS軟件處理所得。具體信息可參考http: //landsat.usgs.gov/documents/cdr_product_guide.pdf[15]。

2.1.2 LEDAPS可用的標準數(shù)據(jù)

標準數(shù)據(jù)的定義來源于USGS EDC或GeoCover的Landsat TM和ETM+數(shù)據(jù)，由LPGS或者NLAPS系統(tǒng)創(chuàng)建，具備有限的元數(shù)據(jù)格式(*_MTL.txt/*.met/*.H1/*.hdr)，其他數(shù)據(jù)格式或LPGS和NLAPS的其他類型均不能使用。含元數(shù)據(jù)(*_MTL.txt)格式的TM和ETM+數(shù)據(jù)部分可下載，數(shù)據(jù)包中含有Landsat單波段數(shù)據(jù)、GCP文件信息及元數(shù)據(jù)文件信息等。在LEDAPS處理中，主要使用數(shù)據(jù)包中的3種數(shù)據(jù)類型，其中元數(shù)據(jù)主要用于自動生成后續(xù)定標、掩模及大氣校正所需的信息。

2.2 非標準數(shù)據(jù)格式處理

LEDAPS只能運行標準數(shù)據(jù)。在國內(nèi)獲取到的非標準Landsat數(shù)據(jù)(相對于標準數(shù)據(jù)格式)需要統(tǒng)一成標準數(shù)據(jù)格式后才能運行LEDAPS，得到地表反射率系列產(chǎn)品。FAST B數(shù)據(jù)包含band1—7.dat和header.dat以及產(chǎn)品信息文件共9個文件。在數(shù)據(jù)格式統(tǒng)一之前需先將國內(nèi)地面站獲取的*.dat數(shù)據(jù)與標準數(shù)據(jù)進行配準和重采樣，轉(zhuǎn)換后投影為UTM，基準面為WGS-84，波段5和波段7的分辨率為30 m，波段6的分辨率為60 m。然后將單波段數(shù)據(jù)轉(zhuǎn)換成單波段GeoTiff格式，再模仿生成元數(shù)據(jù)信息。元數(shù)據(jù)信息主要參考頭文件和產(chǎn)品信息文件，影像對應日期、太陽方位角和高度角、定標參數(shù)、角點經(jīng)緯度坐標和投影坐標(分別在ENVI和ERDAS中讀取)等需設(shè)置，控制點文件信息可以省略。數(shù)據(jù)格式統(tǒng)一后，后續(xù)產(chǎn)品實現(xiàn)參照1.3節(jié)。

3 方法試驗與效果分析

LEDAPS在國內(nèi)的成功使用案例較少，可參照的相關(guān)信息不多。中國科學院凈月潭遙感實驗站有關(guān)于舊版本(ledapsSrc.20110204.tar.gz)源數(shù)據(jù)程序編譯安裝的信息。本文推廣標準數(shù)據(jù)和非標準數(shù)據(jù)的方法試驗用于Landsat影像p120r38不同年份的數(shù)據(jù)轉(zhuǎn)換和格式統(tǒng)一，完成了長時間序列(1987—2011)數(shù)據(jù)預處理及校正前后的評價。表1所列的影像信息除了包含多個時相的標準數(shù)據(jù)和非標準數(shù)據(jù)，還包含多個時相的Landsat CDR數(shù)據(jù)，含地表反射率數(shù)據(jù)，可供直接使用。

表1 長時間序列影像信息(1987—2011)

3.1 方法試驗

3.1.1 標準數(shù)據(jù)轉(zhuǎn)換方法

實現(xiàn)本試驗采用長時間序列影像數(shù)據(jù)中成像日期為2003-01-28的ETM+數(shù)據(jù)，數(shù)據(jù)包中內(nèi)容(參照2.1.2)有 L71120038_03820030128_B10-B80.TIF.L71120038_03820030128_GCP.txt及L71120038_03820030128_MTL.txt。經(jīng)過LEDAPS程序轉(zhuǎn)換得到的數(shù)據(jù)文件(參照1.3)有cloud_mask.log，L71120038_03820030128.carbon_met.txt,lndcal.L71120038_03820030128.hdf，lndcal.L71120038_03820030128.hdf.hdr，lndcal.L71120038_03820030128.txt，lndcsm.L71120038_03820030128.hdf，lndcsm.L71120038_03820030128.hdf.hdr，lndcsm.L71120038_03820030128.txt，lndsr.L71120038_03820030128.hdf，lndsr.L71120038_03820030128.hdf.hdr，lndsr.L71120038_03820030128.txt，lndth.L71120038_03820030128.hdf及l(fā)ndth.L71120038_03820030128.hdf.hdr。

3.1.2 非標準數(shù)據(jù)轉(zhuǎn)換和格式統(tǒng)一方法

首先將國內(nèi)地面站獲取的*.dat數(shù)據(jù)轉(zhuǎn)換成標準數(shù)據(jù)格式，之后采取LEDAPS處理。本試驗采用長時間序列影像數(shù)據(jù)中成像日期為2011-05-18的TM數(shù)據(jù)，數(shù)據(jù)包中內(nèi)容(參照2.2.2)有band1.dat，band2.dat，band3.dat，band4.dat，band5.dat，band6.dat，band7.dat，header.dat，及ProductDescription.self。配準和重采樣數(shù)據(jù)后得到新的數(shù)據(jù)包，投影為UTM，基準面為WGS-84。后續(xù)的處理參照3.1.1。

3.2 效果分析

LEDPAS處理的長時間序列地表反射率影像，其校正前后的影像視覺效果和反射輻射亮度有明顯變化，校正后的影像更加清楚。

限于篇幅，本文在標準數(shù)據(jù)和非標準數(shù)據(jù)中各選一個年份的影像LEDDAS校正和評價(圖2,3)。

圖2 LEDAPS校正前(左)后(右) ETM+3,2,1真彩色合成圖像(2003-01-28)

圖3 LEDAPS校正前(左)后(右) TM3,2,1真彩色合成圖像(2011-05-18)

從圖2(左)可以看出，校正前表觀反射率圖像上地物邊緣模糊，整幅圖像總體呈偏藍色調(diào)，這種輻射失真是由于在可見光波段，大氣窗口內(nèi)的太陽輻射受大氣分子吸收作用的影響較少，藍光波段受大氣瑞利散射影響較大造成的。校正后影像(圖2(右))氣溶膠基本被消除，地物輪廓清晰，地物可識別性高，圖像清晰度明顯提高，基本還原了地物真實反射率，因而圖像的亮度明顯提高，無論是視覺效果還是所反映物理信息的準確性都得到較大改善。

從圖3(左)可以看出，校正前表觀反射率影像整體顏色偏暗，色彩不分明，清晰度不高； 校正后影像(圖3(右))清晰度增加，能夠識別出的地物種類增多，層次更加豐富，影像質(zhì)量、視覺效果均得到較大改善，反演的物理信息更接近實況。特別是校正后地物邊緣變得清晰，影像對比度提高，其中植被的顏色更鮮明，可識別性更高。

將圖2，3中的局部區(qū)域(紅框)放大，能夠更好地比較校正前后的影像效果。

(a) 原始影像(b) TOA反射率影像(c) 地表反射率影像

(a) 原始影像(b) TOA反射率影像(c) 地表反射率影像

為了更好地理解本研究中標準數(shù)據(jù)和非標準數(shù)據(jù)經(jīng)LEDAPS校正前后各波段反射率的變化，在LEDAPS處理得到的TOA和地面反射率影像上選擇典型地物100個，將反射率取平均后與標準波譜曲線比較。圖4和圖5局部影像圖分別對應選取的典型地物水體和植被。由圖6(a)可知，LEDAPS處理后水體的波譜曲線形態(tài)與標準波譜曲線的相似，不同波段對應不同圖像波譜的具體趨勢比較與姚薇等[16]研究結(jié)果相似。由圖6(b)可知，經(jīng)過LEDAPS處理后的植被波譜曲線顯示了典型綠色植被波譜特征，而未處理之前的植被波譜曲線特征不明顯。此外，除處理后的近紅外波段的反射率值低于標準值外，其他各波段的反射率值與標準值基本相同，與姚薇等[16]研究結(jié)果基本相同。分析校正前后各波段的動態(tài)變換范圍，發(fā)現(xiàn)前3個波段的反射率值在校正后顯著減小，而后3個波段的反射率值在大氣校正后有所增大，這與大氣瑞利散射和氣溶膠的散射作用以及水汽的吸收密切關(guān)聯(lián)，與楊靜學等[9]的研究結(jié)果基本吻合。說明標準數(shù)據(jù)和非標準數(shù)據(jù)經(jīng)過LEDAPS處理后得到的地面反射率產(chǎn)品均能有效地降低大氣中臭氧、水汽及氣溶膠等對影像波段反射率的影響，從而更有助于地表真實光譜信息的提取與研究，為遙感分類及建模應用研究提供更精確的數(shù)據(jù)源。

(a) 水體(2003-01-28)(b) 植被(2011-05-18)

分別對比圖3和圖2以及圖5和圖4可以看出，非標準數(shù)據(jù)的圖像處理效果不及標準數(shù)據(jù)的處理效果。一方面，由于LEDAPS預處理方法的使用初期僅處理標準數(shù)據(jù)，本研究為了實現(xiàn)長時間序列反射率圖像統(tǒng)一轉(zhuǎn)換，將國內(nèi)地面站獲取的非標準數(shù)據(jù)統(tǒng)一至標準數(shù)據(jù)的格式，數(shù)據(jù)統(tǒng)一過程中不可避免存在誤差； 另一方面，被處理的Landsat數(shù)據(jù)來源于國內(nèi)外不同的生產(chǎn)機構(gòu)，其初級產(chǎn)品的處理級別有差異。雖存有以上不足，考慮到Landsat連續(xù)時相數(shù)據(jù)成功應用的案例，以及可用標準數(shù)據(jù)的有限性，為了滿足后續(xù)連續(xù)時相遙感影像相關(guān)的應用研究工作，現(xiàn)有的處理結(jié)果具有使用價值。并且在實現(xiàn)長時間序列的影像轉(zhuǎn)換中TM和ETM+影像產(chǎn)品對于真實地物的光譜反應符合波譜規(guī)律，表明兩者具備可比性。

4 結(jié)論

本研究介紹了一種長時間序列Landsat影像預處理程序LEDAPS，詳細描述了該程序的數(shù)據(jù)選擇、格式統(tǒng)一以及算法的實現(xiàn)過程，并在試驗中成功地運用LEDAPS實現(xiàn)了不同來源的Landsat影像數(shù)據(jù)格式的處理，為土地覆蓋變化和干擾等長時間序列的監(jiān)測、生物物理參數(shù)遙感反演提供了科學產(chǎn)品。本研究將有助于促進國內(nèi)形成處理長時間序列影像數(shù)據(jù)的準則。主要結(jié)論如下：

1)LEDAPS反射率轉(zhuǎn)換方法提供了較為全面的參數(shù)用于6S輻射傳輸模型生成地表反射率數(shù)據(jù)。

2)基于LEDAPS 的Landsat影像校正產(chǎn)品的生成比單一的6S模型校正結(jié)果使用更加完備的參數(shù)，結(jié)果精度更有說服力。前者能自動化生成長時間序列產(chǎn)品，后者模型方法不斷成熟并改進，為前者在方法的自動化、模塊化、使用化及工程化的發(fā)展提供詳盡的基礎(chǔ)。

3)結(jié)果評價指出標準和非標準Landsat影像數(shù)據(jù)經(jīng)LEDAPS處理前后反射波譜與標準波譜具有相似的形態(tài)，表明處理后的波譜更接近真實波譜，能有效地降低大氣中O3、水汽及氣溶膠等對影像波段反射率的影響。

Landsat影像以其能提供長時間序列產(chǎn)品的優(yōu)勢，為連續(xù)時相動態(tài)監(jiān)測提供了可能性。但該方法在非標準數(shù)據(jù)參數(shù)設(shè)置、算法驗證以及結(jié)果精度等方面尚有不足之處，是今后研究的重點。

參考文獻(References)：

[1] Cohen W B,Goward S N.Landsat’s role in ecological applications of remote sensing[J].Bioscience,2004,54(6):535-545.

[2] Kennedy R E,Cohen W B,Schroeder T A.Trajectory-based change detection for automated characterization of forest disturbance dynamics[J].Remote Sensing of Environment,2007,110(3):370-386.

[3] Huang C Q,Goward S N,Masek J G,et al.Development of time series stacks of Landsat images for reconstructing forest disturbance history[J].International Journal of Digital Earth,2009,2(3):195-218.

[4] Masek J G,Vermote E F,Saleous N,et al.A Landsat surface reflectance dataset for North America,1990-2000[J].Geoscience and Remote Sensing Letters,2006,3(1):68-72.

[5] Wolfe R,Masek J,Saleous N,et al.LEDAPS:Mapping North American disturbance from the Landsat record[C]//Proceedings of the 2004 IEEE International Geoscience and Remote Sensing Symposium.IEEE,2004.

[6] 亓雪勇,田慶久.光學遙感大氣校正研究進展[J].國土資源遙感,2005，17(4):1-6.

Qi X Y,Tian Q J.The advances in the study of atmospheric correction for optical remote sensing[J].Remote Sensing for Land and Resources,2005，17(4):1-6.

[7] 趙春江,宋曉宇,王紀華,等.基于6S模型的遙感影像逐像元大氣糾正算法[J].光學技術(shù),2007,33(1):11-15.

Zhao C J,Song X Y,Wang J H,et al.An algorithm based on 6S model for removing atmospheric effects form satellite imagery pixel by pixel[J].Optical Technique,2007,33(1):11-15.

[8] 徐永明,覃志豪,陳愛軍.基于查找表的MODIS逐像元大氣校正方法研究[J].武漢大學學報:信息科學版,2010,35(8):959-962.

Xu Y M,Qin Z H,Chen A J.A pixel-by-pixel atmospheric correction algorithm for MODIS data based on look-up table[J].Geomatics and Information Science of Wuhan University,2010,35(8):959-962.

[9] 楊靜學,王云鵬,楊勇.基于高程或氣溶膠厚度與6S模型校正參數(shù)回歸方程的遙感圖像大氣校正模型[J].遙感技術(shù)與應用,2009,24(3):331-340.

Yang J X,Wang Y P,Yang Y.Atmospheric correction model of remote sensing image by regression analysis between elevation/aerosol optical thickness and correction parameters of 6S model[J].Remote Sensing Technology and Application,2009,24(3):331-340.

[10]Vermote E F,Saleous N E,Justice C O,et al.Atmospheric correction of visible to middle-infrared EOS-MODIS data over land surfaces:Background,operational algorithm,and validation[J].Journal of Geophysical Research:Atmospheres(1984-2012),1997,102(D14):17131-17141.

[11]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].IEEE Transactions on Geoscience and Remote Sensing,1997,35(5):1286-1298.

[12]Vermote E F,Saleous N,Justice C O.Atmospheric correction of MODIS data in the visible to middle infrared:First results[J].Remote Sensing of Environment,2002,83:97-111.

[13]Ouaidrari H,Vermote E F.Operational atmospheric correction of Landsat TM data[J].Remote Sensing of Environment,1999,70(1):4-15.

[14]Irish R R.Landsat7 automatic cloud cover assessment:Algorithms for multispectral,hyperspectral,and ultraspectral imagery[J].Proceedings of SPIE,2000,4049:348-355.

[15]Department of the Interior U S.Geological Survey.Landsat climate data record(CDR)surface reflectance product guide version 2.0[EB/OL].[2013-03-27].http://landsat.usgs.gov/documents.

[16]姚薇,李志軍,姚珙,等.Landsat衛(wèi)星遙感影像的大氣校正方法研究[J].大氣科學學報,2011,34(2):251-256.

Yao W,Li Z J,Yao G,et al.Atmospheric correction model for Landsat images[J].Transactions of Atmospheric Sciences,2011,34(2):251-256.