Sentinel-2A衛(wèi)星數(shù)據(jù)處理分析及在干旱河谷提取中的應(yīng)用

楊 斌， 李 丹， 高桂勝， 陳 財， 王 磊

(西南科技大學(xué)環(huán)境與資源學(xué)院，綿陽 621010)

0 引言

哨兵-2A(Sentinel-2A)衛(wèi)星作為歐洲空間局(European Space Agency，ESA)“哥白尼”計劃系列中的一顆重要的光學(xué)遙感衛(wèi)星，于2015年6月23日發(fā)射升空，衛(wèi)星運行在高度為786 km、傾角為98.5°的太陽同步軌道上。Sentinel-2系列計劃發(fā)射2顆衛(wèi)星，構(gòu)建雙星對地觀測平臺，主要用于全球生態(tài)環(huán)境的監(jiān)測[1]，并且該系列衛(wèi)星與Landsat系統(tǒng)在對地觀測領(lǐng)域互為補充，具有較高空間分辨率、免費共享使用等特點，因而對該衛(wèi)星的深入研究與應(yīng)用具有重要意義。

Sentinel-2A衛(wèi)星具有13個光譜波段，涵蓋可見光、近紅外和短波紅外波段，擁有10 m，20 m和60 m3個不同尺度的空間分辨率，重訪周期為10 d，Sentinel -2B衛(wèi)星發(fā)射升空后，Sentinel-2這2顆衛(wèi)星每隔5 d就可對全球N84°～S56°范圍內(nèi)的所有地物覆蓋一遍[2]。Sentinel-2A衛(wèi)星自發(fā)射升空后經(jīng)過近半年的運行和調(diào)試，從2015年12月3日起，正式向全球用戶提供免費下載服務(wù)[3]，國外研究人員前期則曾模擬Sentinel-2數(shù)據(jù)波段進行了小麥土豆的葉面積指數(shù)提取[4]及植被葉綠素含量的評估[5-6]，而后又利用Sentinel-2A衛(wèi)星數(shù)據(jù)進行了溫室效應(yīng)變化監(jiān)測[7]、紅邊光譜指數(shù)適應(yīng)性分析[8]和云檢測分析[9]，而國內(nèi)研究人員對于Sentinel-2A衛(wèi)星數(shù)據(jù)的應(yīng)用研究還不多見。為進一步了解認識Sentinel-2A衛(wèi)星影像數(shù)據(jù)的基本特點，通過對該影像進行數(shù)據(jù)處理與分析，選取研究區(qū)域嘗試運用該數(shù)據(jù)相關(guān)指標進行地學(xué)環(huán)境領(lǐng)域綜合應(yīng)用，為了解該遙感衛(wèi)星影像數(shù)據(jù)的應(yīng)用潛力提供較為全面的基礎(chǔ)信息。

1 Sentinel-2A衛(wèi)星概述

1.1 衛(wèi)星數(shù)據(jù)基本特征

Sentinel-2A衛(wèi)星屬于全新的高空間分辨率多光譜成像衛(wèi)星，主要用于包括陸地植被[11]、土壤以及水資源、內(nèi)河水道和沿海區(qū)在內(nèi)的全球陸地觀測以及緊急救援服務(wù)[12-14]。選取目前應(yīng)用較為廣泛且參數(shù)相近的陸地資源觀測遙感衛(wèi)星Landsat8，分別從波段設(shè)置、波長范圍、幅寬大小、空間分辨率、時間分辨率和輻射分辨率等基本特征參數(shù)對這2種衛(wèi)星數(shù)據(jù)進行對比，如表1所示。

表1 2種衛(wèi)星數(shù)據(jù)主要參數(shù)對比Tab.1 Comparison of spectral bands data between two satellites

從表1中可以得出Sentinel-2A數(shù)據(jù)除在空間分辨率、時間分辨率及幅寬方面具有一定的優(yōu)勢之外，在紅光邊緣波段、近紅外波段均比Landsat8數(shù)據(jù)具有更細致、更全面的波譜分辨率[15-16]，其中，Sentinel-2A波段中的B1是氣溶膠波段，B5，B6，B7和B8a是植被紅邊波段，B9是水汽波段，B10是卷云波段，B11和B12適合于識別提取雪、云、冰信息(圖1)，由此可見，利用Sentinel-2A數(shù)據(jù)在全球植被生態(tài)特征監(jiān)測與分析方面具有更強的應(yīng)用潛力。

圖1 Sentinel-2A衛(wèi)星數(shù)據(jù)波段分布Fig.1 Band distribution of Sentinel-2A satellite

1.2 數(shù)據(jù)組織形式及產(chǎn)品格式

Sentinel-2A衛(wèi)星數(shù)據(jù)的分發(fā)處理模式與Sentinel系列衛(wèi)星相一致，分別通過地面控制中心、地面數(shù)據(jù)處理中心及數(shù)據(jù)分發(fā)中心3個環(huán)節(jié)進行。按照ESA發(fā)布的Sentinel-2A數(shù)據(jù)產(chǎn)品等級依次劃分為L0，L1A，L1B，L1C和L2A這5種不同層次的數(shù)據(jù)格式[17]，其中L0級別為衛(wèi)星獲取的初始壓縮數(shù)據(jù)，并提供生產(chǎn)L1A級別數(shù)據(jù)所需的所有信息，L1A級別是在L0級別基礎(chǔ)上，對所有波段進行解壓，并對短波紅外波段進行重新排列處理； L1B級別數(shù)據(jù)是在L1A級別數(shù)據(jù)的基礎(chǔ)上，經(jīng)過輻射校正，并基于有效的地面控制點構(gòu)建物理幾何精校正模型，并附帶到該級別產(chǎn)品數(shù)據(jù)中； L1C級別數(shù)據(jù)則是在L1B級別數(shù)據(jù)的基礎(chǔ)上，通過重采樣、正射校正和大氣表觀反射率計算處理等操作，生產(chǎn)云、陸地/水體掩模數(shù)據(jù)。其中，L0，L1A和L1B級別的數(shù)據(jù)不對用戶發(fā)布，L1C級別數(shù)據(jù)則通過ESA SciHub網(wǎng)站對用戶免費公開發(fā)布。

ESA對外共享發(fā)布的L1C級別數(shù)據(jù)都是按照100 km2作為最小瓦片單元進行存儲，該數(shù)據(jù)經(jīng)過UTM/WGS84投影正射校正處理，且最基本的文件單元包括影像數(shù)據(jù)文件(IMG_ DATA )、影像數(shù)據(jù)輔助文件(AUX_DATA)、影像質(zhì)量指標文件(QI_DATA)以及元數(shù)據(jù)信息(Metadata)，對L1C級別基本數(shù)據(jù)的說明分析如表2所示。

表2 L1C級別數(shù)據(jù)說明Tab.2 Data illustrate of L1C basic unit

根據(jù)ESA發(fā)布的Sentinel-2A用戶手冊指出，L1C產(chǎn)品的獲取、處理、歸檔及發(fā)布均由Sentinel-2衛(wèi)星數(shù)據(jù)地面數(shù)據(jù)分發(fā)中心提供，而 L2A級別數(shù)據(jù)的生成處理與應(yīng)用分析則需用戶自行利用L1C級別數(shù)據(jù)和Sentinel-2 Toolbox(SNAP軟件工具箱)處理完成。

1.3 L2A級別數(shù)據(jù)生成

通過對Sentinel-2A用戶手冊分析得出，L1C級別數(shù)據(jù)雖然經(jīng)過了正射校正、大氣表觀反射率處理等操作，但在光譜特征、地形變換等方面仍然具有一定的缺陷，而L2A級別數(shù)據(jù)處理則基于查找表法(LUT)利用已經(jīng)計算出的不同復(fù)雜大氣條件下的氣溶膠光學(xué)厚度、反射率和相函數(shù)等指標，進行大氣下墊面發(fā)射率校正處理，最大程度考慮了大氣對遙感數(shù)據(jù)的影響(即大氣精校正)。由于 L2A級別數(shù)據(jù)是在L1C級別數(shù)據(jù)基礎(chǔ)上用戶自行生成，因此需對L2A級別數(shù)據(jù)生成進行探討。

L2A級別數(shù)據(jù)生成是基于L1C級別數(shù)據(jù)調(diào)用SNAP軟件Sentinel-2 Toolbox中的sen2cor處理模塊，進行大氣校正和幾何精校正，生成精度更高的大氣下墊面反射率影像數(shù)據(jù)，還附加生成了氣溶膠厚度分布圖、水汽分布圖和場景分類圖等環(huán)境質(zhì)量分析數(shù)據(jù)[18]。Sentinel-2 Toolbox主要針對于Sentinel-2系列光學(xué)影像數(shù)據(jù)開發(fā)的可視化遙感處理工具箱軟件，其處理操作可通過桌面可視化和命令行運行界面2種方式調(diào)用，并提供了基于Java和 Python語言的開發(fā)接口與環(huán)境，便于用戶自行開發(fā)定制。根據(jù)ESA網(wǎng)站信息及發(fā)布的Sentinel-2衛(wèi)星數(shù)據(jù)用戶指南[19]，Sentinel-2 Toolbox集成在SNAP軟件中，而在處理生成L2A級別數(shù)據(jù)之前，需先安裝SNAP軟件再運用Python語言命令行安裝sen2cor處理模塊，然后在SNAP軟件中調(diào)用sen2cor進行L2A級別數(shù)據(jù)的生成。

生成的L2A級別數(shù)據(jù)除完成大氣下墊面反射率校正(bottom-of-atmosphere，BOA)、卷云計算(cirrus)和地形校正(terrain)等操作外，還將波段保留為10個(去除B1，B9和B10)，且空間分辨率可以根據(jù)需求選擇10 m，20 m和60 m的任何一種類型，數(shù)據(jù)的精度和應(yīng)用價值均得到了顯著提升[19]。此外，SNAP軟件還可利用L2A級別數(shù)據(jù)生成葉面積指數(shù)(leaf area index，LAI)、吸收光合有效輻射比例(fraction of absorbed photosynthetically active radiation，F(xiàn)APAR)、葉綠素含量(chlorophyll content in the leaf，CAB)、冠層含水量(canopy water content，CWC)和植被覆蓋度(fraction of vegetation cover，F(xiàn)VC)等生物物理參數(shù)，可為大面積進行植被生態(tài)環(huán)境的潛力研究與評價提供重要數(shù)據(jù)。

2 研究區(qū)概況及數(shù)據(jù)預(yù)處理

2.1 研究區(qū)概況及數(shù)據(jù)獲取

本文選取四川省茂縣與黑水縣交界區(qū)域為研究區(qū)。該區(qū)經(jīng)緯度坐標范圍為E103°17′～103°44′，N31°47′～31°58′，長41.61 km，寬19.89 km，總面積約848 km2。該區(qū)地處岷江上游黑水河支流區(qū)域，具有山地地形和青藏高原季風(fēng)氣候特征，干旱河谷比較發(fā)育。

遙感數(shù)據(jù)采用2016年5月5日Sentinel-2A L1C級別數(shù)據(jù)，投影/坐標系為UTM/WGS84，整景數(shù)據(jù)的含云量為14.84%，但研究區(qū)范圍內(nèi)云量較少，符合研究需求。

研究區(qū)范圍及其影像如圖2所示。

(a) Sentinel-2A B4(R)，B3(G)，B2(B)真彩色合成影像

(a) 研究區(qū)DEM

圖2研究區(qū)Sentinel-2A數(shù)據(jù)范圍及其遙感影像

Fig.2DistributionofSentinel-2Adataintheresearcharea

2.2 數(shù)據(jù)預(yù)處理

利用裁剪后的L1C級別數(shù)據(jù)，將其導(dǎo)入SNAP軟件中，調(diào)用SNAP軟件中sen2cor處理模塊，選取輸出數(shù)據(jù)精度為10 m，分別對研究區(qū)數(shù)據(jù)進行大氣下墊面發(fā)射率校正、卷云和地形校正處理，生成研究區(qū)L2A級別數(shù)據(jù)及氣溶膠光學(xué)厚度、水汽分布和場景分類信息，如圖3所示。對處理后的L2A影像數(shù)據(jù)各波段基本特征值統(tǒng)計結(jié)果如表3所示。

(a) Sentinel-2A B8(R)，B4(G)，B3(B)假彩色合成影像 (b) 氣溶膠光學(xué)厚度

(c) 水汽分布 (d) 場景分類

從表3可以看出，B2，B3，B4和B8的標準差較大，在第一主成分分量中所占比例也較高，說明其信息量較為豐富，數(shù)據(jù)較為穩(wěn)定，從而可得出紅邊波段在第一主成分分量中所占比例較為穩(wěn)定，對于生態(tài)環(huán)境監(jiān)測具有較好的應(yīng)用價值。

2.3 數(shù)據(jù)產(chǎn)品分析

分別從數(shù)據(jù)生成機理、流程與應(yīng)用效果等方面對L2A級別數(shù)據(jù)生成的氣溶膠光學(xué)厚度、水汽分布和場景分類數(shù)據(jù)類型與影像數(shù)據(jù)的光譜特征進行重點分析[20]。

1)氣溶膠光學(xué)厚度圖層是利用濃密植被算法(dense dark vegetation)，使用短波紅外波段(B12)、紅光波段(B4)和藍光波段(B2)生成，該圖層能有效反演大氣透明度。圖3(b)中可反映出研究區(qū)內(nèi)高山無云區(qū)的氣溶膠光學(xué)厚度較低。

2)水汽分布圖層是利用大氣預(yù)校正微分吸收算法(atmospheric pre-corrected differential absorption，APDA)，使用Sentinel-2A數(shù)據(jù)中的窄近紅外波段(B8a)和水汽波段(B9)生成。其中，將B8a設(shè)置為大氣窗口區(qū)參考通道，將B9 設(shè)置為水汽吸收過程中的測量通道，吸收深度的計算是在假設(shè)測量通道的表面反射率與參考通道相同的情況下獲得。圖3(c)中可反映出研究區(qū)河谷區(qū)內(nèi)水汽含量較大。

3)場景分類圖層是利用閾值測試法(threshold test)和光譜大氣表觀反射率特性，檢測出研究區(qū)云、雪、陰影和地物類型。其中算法可獲得4種不同類型的云類型(低概率云層、中概率云層、高概率云層和薄卷云)，并可生成云層陰影、植被、裸地、水體、雪和黑暗區(qū)域(盲區(qū))等類型。圖3(d)中顯示出了研究區(qū)內(nèi)10種不同的場景類型。

4)獲得的L2A級別影像數(shù)據(jù)波段將原有的13個波段變?yōu)?0個波段。其光譜特征也經(jīng)過了大氣下墊面發(fā)射率校正、卷云和地形校正處理，從生成的20 m空間分辨率L2A級別數(shù)據(jù)[21]中分別獲取研究區(qū)內(nèi)同一地點的水體、植被和裸地3種不同地物光譜特征曲線。圖4中反映出L1C和L2A級別數(shù)據(jù)在光譜特征中的差異性。根據(jù)經(jīng)驗得出，L2A級別數(shù)據(jù)的光譜特性與真實地物光譜特性更為接近。

(a) L1C級別數(shù)據(jù)(b) L2A級別數(shù)據(jù)

圖4典型地物光譜特征分析

Fig.4Spectralcharacteristicsanalysisoftypicalobject

3 生物量提取分析

基于L2A級別數(shù)據(jù)，利用SNAP軟件中生物物理量處理器模塊[22](biophysical processor)生成LAI，F(xiàn)APAR，CAB，CWC和FVC這5種因子，并利用波段運算方法提取出研究區(qū)NDVI(圖5)。該模塊是將事先處理好的區(qū)域植被特征與Sentinel-2數(shù)據(jù)冠層表觀(top of canopy)反射率相關(guān)聯(lián)構(gòu)建一個綜合數(shù)據(jù)庫，并通過訓(xùn)練定標好的神經(jīng)網(wǎng)絡(luò)去評估遙感圖像中每個冠層的特性，而模型中每個神經(jīng)網(wǎng)絡(luò)元的構(gòu)成則涉及到L2A級別數(shù)據(jù)中的B3，B4，B5，B6，B7，B8a，B11，B12和天頂角余弦值、太陽高度角余弦值、相對方位角余弦值等11個輸入?yún)?shù)和5個具有正切S形曲線(sigmoid)傳遞函數(shù)的神經(jīng)元輔助參數(shù)。此模型算法最大程度精細計算出各植被特征與冠層表觀發(fā)射率之間的關(guān)系，可為生態(tài)環(huán)境定量遙感分析提供強有力的數(shù)據(jù)支持。

(a) LAI (b) FAPAR

(c) CAB (d) CWC

(e) FVC (f) NDVI

圖5研究區(qū)植被生物量因子數(shù)據(jù)分布

Fig.5Distributionofvegetable-biophysicalindexdataintheresearcharea

得出的5種生物物理量因子中，LAI是眾多植被-大氣相互作用模型(特別涉及到碳循環(huán)和水循環(huán)模型)中的一個關(guān)鍵參數(shù)[23]，圖5(a)反映出河谷地帶的LAI顯著低于海拔較高的山坡及山頂?shù)貐^(qū)； FAPAR是描述植被結(jié)構(gòu)以及與之相關(guān)的物質(zhì)與能量交換過程的基本生理變量，也是遙感估算陸地生態(tài)系統(tǒng)植被第一性生產(chǎn)力(net primary productivity,NPP)的重要參數(shù)，該指標也能有效表征植被生長狀態(tài)和演化過程，圖5(b)中也反映出云雪與植被的顯著反差，還能有效顯示出河谷地帶植被長勢較差； CAB能夠真實有效地指示出地表植被的生長和營養(yǎng)狀況(受氮素脅迫程度)、光合作用能力、所處的生長發(fā)育階段及是否受到自然環(huán)境破壞或病蟲害影響等狀態(tài)信息，圖5(c)中CAB對植被特征區(qū)分精細度較高，并反映出黑水河上游區(qū)域(左下方)植被長勢較好； CWC能夠有效評估植被的生長狀況，對干旱監(jiān)測和生態(tài)環(huán)境改善等具有現(xiàn)實意義，但圖5(d)中CWC對于區(qū)域植被特征顯示不明顯； 而FVC直接表征地表植被覆蓋的多少，能大致反映出地表資源與生態(tài)環(huán)境的好壞程度，圖5(e)中FVC能顯著區(qū)分出裸露地區(qū)與植被覆蓋度較高地區(qū)的差異性。綜上所述，這5種因子對于區(qū)域植被情況特征分布均具有一定的指示作用，但具體的應(yīng)用效果還有待進一步探討。

4 干旱河谷提取與分析

由于該研究區(qū)位于青藏高原東緣地帶的岷江上游流域，長年受西南橫斷山區(qū)季風(fēng)氣候的影響，氣候溫暖、干燥度與輻射度強、雨熱同步等特點，且河谷區(qū)內(nèi)植被發(fā)育情況較高海拔山坡差，形成了山地典型垂直分異植被特性。從遙感假彩色圖像上可明顯看出干旱河谷區(qū)內(nèi)影像呈現(xiàn)出暗灰色調(diào)，其邊界在圖3(a)中的顯示較為明顯。

通過Sentinel-2A提取出的5種生物物理量因子和NDVI值能有效反映研究區(qū)內(nèi)植被發(fā)育長勢情況，對干旱河谷邊界的定量分析與提取具有較好的指示作用[24]。選取基于知識經(jīng)驗的決策樹分類法，將上述5種生物物理量因子和NDVI參數(shù)導(dǎo)入到ENVI軟件中進行圖像特征分析，如表4所示。

表4 各植被生物量圖像特征分析Tab.4 Analysis of image feature based on every vegetable-biophysical data

結(jié)合各指標影像特征反映得出標準差較高且信息量較大的因子為CAB。并通過查閱該區(qū)域干旱河谷前期研究資料和野外實際調(diào)研[25]，得出研究區(qū)范圍內(nèi)干旱河谷分布與海拔高度分布具有極其緊密的聯(lián)系，該區(qū)干旱河谷海拔主要集中分布在600～2 400 m區(qū)間的河谷范圍內(nèi)[25]。利用研究區(qū)30 m空間分辨率的DEM(從美國USGS網(wǎng)站下載)和CAB因子，分別對比假彩色影像選擇該因子數(shù)據(jù)的合理閾值(本文閾值選為32.368 9)，在ENVI軟件中調(diào)用決策樹分類法，按圖6所示流程定量提取出研究區(qū)的干旱河谷分布邊界與面積，其結(jié)果如圖7所示。

圖6 決策樹分類流程Fig.6 Flowchart of decision tree classification

圖7 干旱河谷提取分布Fig.7 Distribution of dry river valleys

通過上述分析與提取，將決策樹分類結(jié)果導(dǎo)入到ArcGIS 軟件中，運用ArcGIS空間分析功能，計算出該研究區(qū)內(nèi)干旱河谷的面積為114.147 7 km2，占研究區(qū)總面積的13.46%。對干旱河谷空間分布進行的分析表明，干旱河谷區(qū)主要集中分布在黑水河河谷坡度小于25°的范圍內(nèi)。對照該區(qū)域的地質(zhì)圖，可以看出干旱河谷區(qū)內(nèi)灰?guī)r、千枚巖和石英巖較為發(fā)育。結(jié)合包維楷等[26]、楊兆平等[27]前期研究成果與作者以往的研究積累及野外實地考察，選取6個野外觀測樣點(圖7)，進行實地觀測驗證，提取精度達到90%以上。研究結(jié)果表明，干旱河谷提取結(jié)果能較好表征干旱河谷分布情況，為岷江上游干旱河谷大面積國土資源普查和生態(tài)環(huán)境變遷分析研究提供了方法支持。

5 結(jié)論

Sentinel-2A光學(xué)遙感衛(wèi)星作為對地觀測陸地資源衛(wèi)星系列大家庭的一員，以其涵蓋面積大、時間和空間分辨率高、數(shù)據(jù)免費共享等優(yōu)勢受到廣大遙感應(yīng)用工作者的廣泛關(guān)注。本次遙感應(yīng)用的主要結(jié)論如下：

1)從Sentinel-2A遙感數(shù)據(jù)的基本屬性、光譜特征、產(chǎn)品特征和應(yīng)用特征4個方面進行系統(tǒng)分析，利用SNAP軟件sen2cor模塊進行數(shù)據(jù)處理，處理結(jié)果可以較好滿足后續(xù)應(yīng)用需要。

2)以四川省岷江上游流域黑水河區(qū)域為研究區(qū)進行了遙感生物量因子提取分析，并利用分析后的生物量因子和DEM數(shù)據(jù)，采用決策樹分類法模型進行了干旱河谷區(qū)域定量提取。結(jié)果表明，Sentinel-2A系列衛(wèi)星數(shù)據(jù)光譜特征豐富、應(yīng)用前景廣闊、數(shù)據(jù)可擴展性較強，可為植被生態(tài)環(huán)境定量分析和地學(xué)綜合應(yīng)用提供更多的遙感數(shù)據(jù)資源。

志謝： 感謝ESA提供的Sentinel-2A數(shù)據(jù)資源及學(xué)習(xí)文檔，感謝http: //forum.step. esa.int/c/s2tbx/sen2cor論壇網(wǎng)站上國內(nèi)外同行給予的幫助！