中國(guó)區(qū)域Landsat 地表反射率產(chǎn)品

彭燕，何國(guó)金，張兆明，尹然宇

1. 中國(guó)科學(xué)院空天信息創(chuàng)新研究院，北京 100094

2. 中國(guó)科學(xué)院大學(xué)，北京 100049

3. 海南省地球觀測(cè)重點(diǎn)實(shí)驗(yàn)室，海南三亞 572029

引 言

Landsat 系列衛(wèi)星自1972 年發(fā)射以來(lái)，已獲取了大量中空間分辨率的衛(wèi)星影像，這些影像記錄著人類(lèi)活動(dòng)和自然變化，成為最長(zhǎng)時(shí)間序列的星載陸地觀測(cè)數(shù)據(jù)集[1]。地表反射率常用于土地覆蓋及變化研究，是很多地表地球物理參數(shù)（如葉面積指數(shù)、葉綠素和生物量）反演的關(guān)鍵輸入?yún)?shù)。同時(shí)，隨著對(duì)地觀測(cè)領(lǐng)域進(jìn)入大數(shù)據(jù)時(shí)代，衛(wèi)星遙感數(shù)據(jù)應(yīng)用用戶(hù)和研究者不僅希望遙感數(shù)據(jù)在幾何位置上具有一致性，同時(shí)也對(duì)遙感數(shù)據(jù)的輻射一致性提出了更高的要求，以更好地進(jìn)行遙感數(shù)據(jù)應(yīng)用和信息挖掘分析。由于衛(wèi)星遙感數(shù)據(jù)的大氣校正是一項(xiàng)繁瑣而專(zhuān)業(yè)的工作，因此，研究生產(chǎn)高質(zhì)量中國(guó)區(qū)域的Landsat 地表反射率產(chǎn)品并公開(kāi)共享，具有重要的現(xiàn)實(shí)意義。

本文提供30 m 空間分辨率的中國(guó)區(qū)域Landsat 系列的地表反射率產(chǎn)品，給出了中國(guó)區(qū)域Landsat數(shù)據(jù)地表反射率產(chǎn)品生產(chǎn)的技術(shù)流程：首先采集覆蓋中國(guó)區(qū)域的Landsat 正射影像，然后經(jīng)過(guò)輻射定標(biāo)，獲取6S（Second Simulation of Satellite Signal in the Solar Spectrum）輻射傳輸模型所需要的幾何參數(shù)和大氣參數(shù)并進(jìn)行大氣校正，生成中國(guó)區(qū)域Landsat 地表反射率產(chǎn)品。

1 數(shù)據(jù)采集和處理方法

1.1 數(shù)據(jù)采集方法

本文所生產(chǎn)的中國(guó)區(qū)域地表反射率產(chǎn)品主要是在Landsat 系列正射影像數(shù)據(jù)的基礎(chǔ)上進(jìn)行生產(chǎn)的，1980s-2012 年為L(zhǎng)andsat 5 數(shù)據(jù)，2000-2003 年期間為無(wú)條帶的Landsat 7 數(shù)據(jù)，2013-2019 年采用Landsat 8 數(shù)據(jù)，這些數(shù)據(jù)主要由中國(guó)遙感衛(wèi)星地面站接收。對(duì)于太陽(yáng)天頂角大于76°的Landsat 數(shù)據(jù)不進(jìn)行大氣校正處理。

1.2 數(shù)據(jù)處理方法

中國(guó)區(qū)域地表反射率產(chǎn)品基于6S 輻射傳輸模型進(jìn)行大氣校正獲得。6S 模型是目前比較完善的大氣校正模型之一，由Vermote 等人在5S 模型的基礎(chǔ)上改進(jìn)而來(lái)[2]，適用于0.25-4 μm 波長(zhǎng)范圍內(nèi)電磁波的大氣輻射傳輸?shù)哪M。假設(shè)氣溶膠類(lèi)型為大陸型，并將地表視為平面朗伯體，則其地表反射率可根據(jù)表觀反射率計(jì)算而得到。其公式為[3]：

式中，ρs為地表反射率；ρTOA為表觀反射率，通過(guò)輻射定標(biāo)可計(jì)算得到；Tg(OG)為O2、O3、CO2、NO2、CH4氣體總的透過(guò)率，Tg(OG)=Tg(O2)×Tg(O3)×Tg(CO2)×Tg(NO2)×Tg(CH4)；ρR+A為瑞利和氣溶膠的反射率；TR+A為瑞利和氣溶膠透過(guò)率；Tg(H2O)為水汽的透過(guò)率；SR+A為瑞利和氣溶膠球面反照率。而Tg(OG)、ρR+A、TR+A、Tg(H2O)、SR+A等大氣相關(guān)系數(shù)均可通過(guò)調(diào)用6S 模型計(jì)算得到。而基于6S 模型進(jìn)行大氣校正的關(guān)鍵和難點(diǎn)在于大氣參數(shù)的獲取。本數(shù)據(jù)集是基于公式（1）的計(jì)算原理和6S 輻射傳輸模型[1,3-4]進(jìn)行生產(chǎn)的，首先需要進(jìn)行輻射定標(biāo)，然后獲取6S 模型所需要的大氣參數(shù)、DEM 以及幾何參數(shù)進(jìn)行大氣校正，其技術(shù)路線(xiàn)如圖1 所示。

1.2.1 輻射定標(biāo)

表觀反射率與大氣頂層進(jìn)入衛(wèi)星傳感器的光譜輻射亮度、日地間距離、大氣頂層的平均太陽(yáng)光譜輻照度，以及太陽(yáng)的天頂角有關(guān)。Landsat 5/7 的表觀反射率計(jì)算公式為：

式中，ρ為表觀反射率，d為日地距離，θz為太陽(yáng)天頂角（單位為度，與元數(shù)據(jù)文件中給出的太陽(yáng)高度角互為余角），ESUN（Mean solar exoatmospheric spectral irradiances）是大氣層頂平均太陽(yáng)光譜輻照度，Lλ為光譜輻射亮度，單位為w/(m2?μm?sr)，利用定標(biāo)系數(shù)進(jìn)行線(xiàn)性計(jì)算得到。

圖1 地表反射率產(chǎn)品生產(chǎn)技術(shù)路線(xiàn)圖

Landsat 8 的表觀反射率計(jì)算公式為：

式中，Mρ和Aρ可從元數(shù)據(jù)文件中獲取，分別為Mρ= REFLECTANCE_MULT_BAND_x，Aρ= REFLECTANCE_ADD_BAND_x，x 為波段號(hào)。

1.2.2 模型參數(shù)獲取

6S 模型的輸入?yún)?shù)主要包括太陽(yáng)天頂角、方位角和觀測(cè)天頂角、觀測(cè)方位角等觀測(cè)幾何參數(shù)以及水汽、臭氧、氣溶膠、氣壓等大氣參數(shù)和海拔高度信息。在本文中，觀測(cè)幾何參數(shù)可根據(jù)元數(shù)據(jù)文件獲取，海拔高度信息則是利用0.05°空間分辨率的GCM DEM（Global Climate Model Digital Elevation Model）[5]得到的（https://search.earthdata.nasa.gov/search?q=ASTER&ok=ASTER）。對(duì)于Landsat 5/7 數(shù)據(jù)而言，大氣校正所需要的水汽數(shù)據(jù)以及氣壓數(shù)據(jù)均來(lái)自于NCEP（NOAA National Centers for Environmental Prediction）再分析數(shù)據(jù)（http://dss.ucar.edu/），空間分辨率為2.5°×2.5°，臭氧 數(shù) 據(jù) 來(lái) 自 TOMS （ Total Ozone Mapping-Spectrometer ） 數(shù) 據(jù)（https://ozoneaq.gsfc.nasa.gov/data/ozone/#），空間分辨率為1.25°×1.00°，氣溶膠光學(xué)厚度則是利用暗目標(biāo)法反演得到的[3,6]。而對(duì)于Landsat 8 數(shù)據(jù)而言，大氣校正所需要的水汽數(shù)據(jù)、臭氧數(shù)據(jù)、氣溶膠光學(xué)厚度（AOT）來(lái)自于0.05°空間分辨率的MODIS09CMA（MODIS Aerosol Optical Thickness Climate Modeling Grid）和MODIS09CMG（MODIS surface reflectance Climate Modeling Grid）數(shù)據(jù)[7]（https://earthdata.nasa.gov/earth-observation-data/near-real-time/download-nrt-data/modis-nrt）。

2 數(shù)據(jù)樣本描述

本數(shù)據(jù)集是以景為單位存放文件夾，由每個(gè)波段的地表反射率產(chǎn)品、質(zhì)量文件（Quality Assessment, QA）、元數(shù)據(jù)、縮略圖組成。

（1）文件夾的命名規(guī)則為衛(wèi)星-傳感器-path-row-成像日期-LSR，如L5-TM-115-026-19840418-LSR。

（2）每個(gè)波段的地表反射率產(chǎn)品：Landsat 5/7 包括波段1，2，3，4，5，7（即分別為藍(lán)、綠、紅、近紅、中紅外1、中紅外2），Landsat 8 包括波段1，2，3，4，5，6，7（即分別為深藍(lán)、藍(lán)、綠、紅、近紅、中紅外1、中紅外2）。命名規(guī)則為衛(wèi)星-傳感器-path-row-成像日期-LSR-BX（X 表示第幾個(gè)波段）.TIF，如L5-TM-115-026-19840418-LSR-B1.TIF。影像的空間分辨率為30 m，投影坐標(biāo)系是WGS84 UTM。為了降低存儲(chǔ)空間，地表反射率結(jié)果由原本的0-1 范圍內(nèi)的浮點(diǎn)型均乘以10000變成16 位整型，背景填充值為-9999，并進(jìn)行了“LZW”的無(wú)損壓縮。

（3）質(zhì)量文件（QA）：包括在原始數(shù)據(jù)的基礎(chǔ)上生成的QA，為PIXEL-QA，主要是對(duì)填充值（Fill）、晴空（Clear）、云（Cloud）、云置信度（Cloud Confidence）、云陰影（Cloud Shadow）、冰雪（Snow/Ice）以及水（Water）等信息進(jìn)行標(biāo)識(shí)，命名規(guī)則為：衛(wèi)星-傳感器-path-row-成像日期-PIXELQA.TIF。和在地表反射率的基礎(chǔ)上生成的QA，主要是對(duì)氣溶膠相關(guān)的信息進(jìn)行標(biāo)識(shí)，Landsat 5/7 為L(zhǎng)SR-CLOUD-QA，命名規(guī)則為衛(wèi)星-傳感器-path-row-成像日期-LSR-CLOUD-QA.TIF；Landsat 8 為L(zhǎng)SR-AEROSOL，命名規(guī)則為衛(wèi)星-傳感器-path-row-成像日期-LSR-AEROSOL.TIF。Landsat 5/7 和Landsat 8 的QA 的屬性分別如表1-4 所示。

表1 Landsat 5/7 的PIXEL-QA 屬性表

表2 Landsat 5/7 的LSR-CLOUD-QA 的屬性表

表3 Landsat 8 的PIXEL-QA 屬性表

表4 Landsat 8 的LSR-AEROSOL 屬性表

（4）元數(shù)據(jù)文件：包括原始數(shù)據(jù)的元數(shù)據(jù)文件（MTL.txt 文件），命名規(guī)則為衛(wèi)星-傳感器-pathrow-成像日期-MTL.txt；和地表反射率產(chǎn)品的XML 元數(shù)據(jù)文件，命名規(guī)則為衛(wèi)星-傳感器-path-row-成像日期-LSR.xml。XML 元數(shù)據(jù)文件描述了地表反射率產(chǎn)品的基本信息以及各個(gè)波段產(chǎn)品的相關(guān)信息。

（5）縮略圖：包括512 像素大小和1024 像素大小的縮略圖，命名規(guī)則分別為衛(wèi)星-傳感器-pathrow-成像日期-LSR-THUMB.JPG 和衛(wèi)星-傳感器-path-row-成像日期-LSR-BROWSER.JPG。

圖2 為2005 年江西省的地表反射率結(jié)果示意圖，其中圖2（a）為大氣校正前的影像，圖2（b）為地表反射率結(jié)果，圖2（c）為每景數(shù)據(jù)的成像時(shí)間，如0420 表示2005 年4 月20 日的Landsat5 數(shù)據(jù)。從圖2 中可以看出，經(jīng)過(guò)大氣校正后能較大地消除大氣的影響，具有較好的輻射一致性，尤其是2005 年6 月23 日這一景數(shù)據(jù)，其地表反射率結(jié)果明顯消除了大部分薄云薄霧的影響。

圖2 江西省地表反射率產(chǎn)品示意圖

3 數(shù)據(jù)質(zhì)量控制和評(píng)估

為了評(píng)估本數(shù)據(jù)集的可靠性，在全國(guó)范圍內(nèi)均勻選取了12 個(gè)軌道號(hào)如圖3 所示，每個(gè)軌道號(hào)分別選取了Landsat5/7/8 地表反射率（具體的時(shí)間見(jiàn)圖3 的列表），同時(shí)每景影像均隨機(jī)選取50 個(gè)樣本點(diǎn)，與USGS 發(fā)布的地表反射率產(chǎn)品進(jìn)行對(duì)比驗(yàn)證。結(jié)果顯示本數(shù)據(jù)與USGS 所發(fā)布的地表反射率產(chǎn)品的相關(guān)性系數(shù)R2均在0.99 以上，均方根誤差（RMSE）均小于1%，這表明本數(shù)據(jù)集與USGS所發(fā)布的地表反射率產(chǎn)品具有很好的一致性，如圖4 所示。Feng 等人將USGS 發(fā)布的全球2000 年和2005 年的Landsat 5/7 地表反射率產(chǎn)品分別與相應(yīng)的MODIS 地表反射率產(chǎn)品進(jìn)行比較分析，結(jié)果表明Landsat 5 TM 的均方根差（RMSD）約為2.2%-3.5%，Landsat 7 ETM+的均方根差約為1.3%-2.8%[8]。同時(shí)利用2014 年6 月11 日獲取的南京地區(qū)的實(shí)測(cè)光譜數(shù)據(jù)對(duì)Landsat 8 數(shù)據(jù)地表反射率產(chǎn)品進(jìn)行驗(yàn)證，結(jié)果顯示其地表反射率產(chǎn)品的RMSD 約為2%-4%之間。并對(duì)USGS 發(fā)布的Landsat8地表反射率產(chǎn)品進(jìn)行驗(yàn)證，結(jié)果顯示USGS 發(fā)布的Landsat8地表反射率結(jié)果的RMSD分別為4.07%、4.37%、3.49%和5.3%，而本產(chǎn)品的RMSD 分別為2.71%、3.1%、2.71%、2.08%，本產(chǎn)品的RMSD均低于USGS 發(fā)布的地表反射率產(chǎn)品[1]。

4 數(shù)據(jù)價(jià)值

本文推出的1986-2019 年中國(guó)區(qū)域Landsat 系列衛(wèi)星長(zhǎng)時(shí)序地表反射率產(chǎn)品，采用目前國(guó)際上較為成熟的地表反射率算法，精度較高，未來(lái)將持續(xù)更新，在森林、水資源、氣候變化等領(lǐng)域長(zhǎng)時(shí)序信息挖掘分析方面具有重要的應(yīng)用價(jià)值。

圖3 對(duì)比驗(yàn)證的軌道號(hào)分布示意圖

5 數(shù)據(jù)使用方法和建議

本數(shù)據(jù)集可通過(guò)地球大數(shù)據(jù)科學(xué)工程（ CASEarth ） Databank 在線(xiàn)服務(wù)網(wǎng)址（http://databank.casearth.cn）獲取數(shù)據(jù)。用戶(hù)注冊(cè)成功并登錄系統(tǒng)后，進(jìn)入平臺(tái)產(chǎn)品查詢(xún)界面，產(chǎn)品類(lèi)型選擇“LSR”，然后在所需要的數(shù)據(jù)下點(diǎn)擊下載，隨后會(huì)彈出一個(gè)下載對(duì)話(huà)框，該框中列出了各波段的地表反射率結(jié)果、QA 文件和元數(shù)據(jù)文件，根據(jù)需要下載數(shù)據(jù)即可。用戶(hù)可根據(jù)行政區(qū)、地圖選擇以及行列號(hào)等方式查詢(xún)所需要的數(shù)據(jù)。目前共享的產(chǎn)品主要包括1980s-2012 年云量小于50%的Landsat 5、2000-2003 年云量小于20%的Landsat 7 以及2018-2019 年的全部Landsat 8 的地表反射率產(chǎn)品，后續(xù)作者將持續(xù)生產(chǎn)我國(guó)區(qū)域云量大于50%的Landsat 5、云量大于20%的Landsat 7 和2019年以后的Landsat 8 地表反射率產(chǎn)品，以提供更好的、持續(xù)的數(shù)據(jù)共享服務(wù)。若平臺(tái)系統(tǒng)中暫時(shí)缺少用戶(hù)所需的數(shù)據(jù)或者有與本數(shù)據(jù)相關(guān)的其他數(shù)據(jù)需求，可通過(guò)咨詢(xún)本文作者進(jìn)行申請(qǐng)。

圖4 與USGS 地表反射率產(chǎn)品對(duì)比驗(yàn)證結(jié)果

致 謝

衷心感謝王桂周和龍騰飛在產(chǎn)品生產(chǎn)時(shí)遇到的存儲(chǔ)與效率問(wèn)題給予的建設(shè)性意見(jiàn)。

數(shù)據(jù)作者分工職責(zé)

彭燕（1988—），女，湖南省郴州市人，在讀博士，工程師，研究方向?yàn)檫b感圖像智能處理。主要承擔(dān)工作：算法集成程序編寫(xiě)，數(shù)據(jù)生產(chǎn)流程設(shè)計(jì)，論文撰寫(xiě)。

何國(guó)金（1968—），男，福建省龍巖市人，博士，研究員，研究方向?yàn)檫b感數(shù)據(jù)智能處理與信息挖掘。主要承擔(dān)工作：總體思路與方案設(shè)計(jì)，論文修改。

張兆明（1980—），男，河南省鄭州市人，博士，副研究員，研究方向?yàn)檫b感數(shù)據(jù)智能處理與信息提取。主要承擔(dān)工作：技術(shù)指導(dǎo)，論文修改。

尹然宇（1996—），男，山東省臨沂市人，在讀博士，研究方向?yàn)檫b感圖像智能處理。主要承擔(dān)工作：數(shù)據(jù)挑選與整合。