Sentinel-2A與Landsat 8 OLI逐像元輻射歸一化方法研究

徐玉雯，張 浩，陳正超，景海濤

(1.河南理工大學(xué)測(cè)繪與國(guó)土信息工程學(xué)院，河南 焦作 454003；2.中國(guó)科學(xué)院遙感與數(shù)字地球研究所數(shù)字地球重點(diǎn)實(shí)驗(yàn)室，北京 100094)

隨著在軌運(yùn)行衛(wèi)星數(shù)量和種類(lèi)的增加，越來(lái)越多的遙感應(yīng)用傾向于綜合應(yīng)用多源遙感數(shù)據(jù)開(kāi)展高頻次的應(yīng)用監(jiān)測(cè)，這就要求不同傳感器數(shù)據(jù)解譯具有一致性，通過(guò)輻射歸一化使數(shù)據(jù)在空間和時(shí)間上保持連續(xù)性，使定量研究的重點(diǎn)聚焦在地表的變化而不過(guò)分考慮傳感器的差異。目前常用的輻射歸一化方法可分為絕對(duì)輻射歸一化(absolute radiometric normalization,ARN)及相對(duì)輻射歸一化(relative radiometric normalization,RRN)[1]。絕對(duì)輻射歸一化是將圖像的DN值轉(zhuǎn)換成地表反射率。相對(duì)輻射歸一化是以同區(qū)域多時(shí)相相似波段灰度值之間滿足線性關(guān)系為前提，主要分為基于分布的相對(duì)輻射歸一化方法和基于像元對(duì)的輻射歸一化[2]?；诜植嫉南鄬?duì)輻射歸一化方法簡(jiǎn)單、客觀，考慮全局分布特征，歸一化結(jié)果受主觀因素影響較少，在多源遙感器之間的輻射歸一化更顯優(yōu)勢(shì)[3]。基于像元對(duì)相對(duì)輻射歸一化方法主要原理是從待歸一化的影像對(duì)中選取部分或全部的像元對(duì)作為樣本點(diǎn)回歸得到校正系數(shù)，而樣本點(diǎn)質(zhì)量的好壞直接影響歸一化校正的結(jié)果[4-9]。上述相對(duì)輻射校正方法，不同影像相似波段之間的變換均屬于線性變換，即同一波段所有像元具有統(tǒng)一的變換系數(shù)，不同的是選取偽不變特征點(diǎn)的方式及回歸方程的擬合方法有所差別，并且這些方法多應(yīng)用于同類(lèi)型傳感器多時(shí)相影像之間的歸一化，沒(méi)有考慮不同傳感器對(duì)相同地物的成像差異及不同地物之間的光譜差異。文獻(xiàn)[10]考慮地物差異利用表觀反射率光譜庫(kù)，通過(guò)回歸分析與線性組合變換實(shí)現(xiàn)了SPOT 5與THEOS之間的輻射歸一化，但是該歸一化是基于表觀輻射亮度進(jìn)行的，沒(méi)有在地表反射率層面上實(shí)現(xiàn)輻射歸一化[10]。本文考慮不同傳感器光譜響應(yīng)差異對(duì)地物反射率的影響，以兩景幾乎同時(shí)成像的Sentinel-2A和Landsat 8 OLI數(shù)據(jù)為例，基于光譜庫(kù)構(gòu)建兩傳感器VNIR相似波段之間的匹配轉(zhuǎn)換模型，以實(shí)現(xiàn)逐像元?dú)w一化。

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

本文采用準(zhǔn)同步過(guò)境Sentinel-2A及Landsat 8 OLI為數(shù)據(jù)源，以內(nèi)蒙古中東部達(dá)里諾爾湖附近地區(qū)為研究區(qū)，兩景影像在2017年7月17日成像時(shí)間分別為3:05:41、2:58:44，相差7 min。Sentinel-2是全球環(huán)境與安全監(jiān)測(cè)計(jì)劃(GMES)下的高空間分辨率衛(wèi)星，包括A、B兩星影像，單景成像幅寬與用戶分發(fā)幅寬分別為290、109.8 km，光譜范圍為0.4～2.4 μm，包含13個(gè)波段，空間分辨率分為10、20、60 m 3個(gè)等級(jí)[11]，其Level-1C數(shù)據(jù)是經(jīng)正射校正和亞像元級(jí)幾何精校正后的表觀反射率產(chǎn)品。Landsat 8 OLI陸地成像儀成像寬幅為185 km，包括9個(gè)波段，空間分辨率為30 m，其中全色波段15 m，L1T產(chǎn)品是使用地面控制點(diǎn)和數(shù)字高程模型數(shù)據(jù)進(jìn)行精確校正后的DN值數(shù)據(jù)[12]。同作為提供全球免費(fèi)觀測(cè)服務(wù)的新型中高分辨率傳感器，兩者在定量遙感方面應(yīng)用廣泛，波段設(shè)置也較為相似，VNIR波段光譜響應(yīng)函數(shù)如圖1所示，波段參數(shù)對(duì)比見(jiàn)表1。

表1 Sentitle-2A與Landsat 8 OLI可見(jiàn)—近紅外波段對(duì)比

數(shù)據(jù)預(yù)處理主要包括大氣校正、無(wú)效值掩模、空間重采樣、裁剪等。Sentinel-2A影像采用歐空局官方提供的L2A_AtmCorr(Sentinel-2 atmospheric correction)算法反演地表反射率，其中氣溶膠和水汽的反演分別基于改進(jìn)的DDV(dense dark vegetation)算法及APDA(atmospheric pre-corrected differential absorption)算法，校正得到的地表反射率產(chǎn)品以整型存儲(chǔ)，放大系數(shù)為10 000，有效值范圍0～10 000[13]。Landsat 8 OLI影像采用NASA官方提供大氣校正算法進(jìn)行大氣校正，該算法是基于MODIS影像獲取相應(yīng)大氣參數(shù)經(jīng)反演得到地表反射率[14]，放大系數(shù)為10 000，地表反射率可能值為-2000～12 000，有效值范圍0～10 000。無(wú)效值對(duì)影像歸一化及后續(xù)分析并無(wú)意義，需要對(duì)Sentinel-2A與Landsat 8圖像上無(wú)效反射率值進(jìn)行掩模處理。然后將Sentinel-2A的藍(lán)、綠、紅波段10 m與近紅外20 m的空間分辨率采用像素聚合方法重采樣到與Landsat 8相同的空間分辨率，根據(jù)二者同時(shí)刻成像重疊情況，從配準(zhǔn)好的兩影像中裁剪出600×600像元大小的圖像作為研究區(qū)，如圖2所示。研究區(qū)中包含了濃密植被、草地、高亮白沙地、水體等主要地物。

2 逐像元?dú)w一化方法

2.1 光譜匹配因子庫(kù)建立

光譜庫(kù)是由光譜成像儀在實(shí)驗(yàn)室或野外對(duì)不同地物按照嚴(yán)格規(guī)范的流程取的各類(lèi)地物反射光譜數(shù)據(jù)的集合，具有多樣性、穩(wěn)定性的特點(diǎn)，對(duì)遙感圖像的解譯、影像地物的匹配識(shí)別具有重要的作用。ENVI中內(nèi)嵌的光譜庫(kù)包含了國(guó)際上常用的5種光譜庫(kù)：USGS、JPL、JHU、IGCP-264、ASTER，近2000條光譜，但并不是所有的地物光譜都涵蓋VNIR波長(zhǎng)范圍的反射率值。對(duì)經(jīng)過(guò)篩選的反射率光譜整合、重采樣，通過(guò)與Sentinel-2A和Landsat 8 OLI的VNIR波段的光譜響應(yīng)函數(shù)按照式(1)進(jìn)行等效計(jì)算，形成不同地物類(lèi)別的等效光譜庫(kù)。同一地物因傳感器光譜特性差異反映出的波段等效值也不同，兩傳感器等效光譜庫(kù)對(duì)應(yīng)地物之間的比值則組成了光譜匹配因子庫(kù)。

(1)

式中，ρj為某一地物波段j的等效反射率；λ1、λ2為該波段的起止波長(zhǎng)；Rj(λ)為波段j的響應(yīng)函數(shù)，ρ對(duì)應(yīng)該波段光譜庫(kù)實(shí)際反射率值，Rj(λ)和ρ是間隔1 nm的不連續(xù)的值，對(duì)其積分通常用累加值替代。

另外，由于ENVI自帶光譜庫(kù)中缺少清澈水體的標(biāo)準(zhǔn)光譜，試驗(yàn)中又無(wú)實(shí)測(cè)類(lèi)似的光譜數(shù)據(jù)。又因達(dá)里諾爾湖地處內(nèi)蒙古中東地區(qū)，水體受污染和破壞的程度較小，查詢往年同地區(qū)相同季節(jié)的Hyperion高光譜數(shù)據(jù)，下載得到2012年7月27日一景影像。經(jīng)過(guò)壞波段剔除、壞線條紋去除、減弱smlie效應(yīng)、大氣校正等處理后獲取可見(jiàn)近紅外的48個(gè)有效波段，在達(dá)里諾爾湖上抽樣100個(gè)水體的像元點(diǎn)平均以獲取連續(xù)的水體光譜，錄入光譜庫(kù)，與光譜庫(kù)其他光譜同作處理。

2.2 逐像元?dú)w一化模型

傳感器成像過(guò)程中，每個(gè)像元對(duì)應(yīng)的實(shí)際地物往往是由一種或多種地物混合而成的，反映出的地物光譜反射率值也是一種或多種地物光譜的組合。通過(guò)完全約束最小二乘約束對(duì)Sentinel-2A圖像進(jìn)行線性分解獲取典型地物的豐度圖像，豐度圖像客觀地表示了每種地物在每個(gè)像元中所占的比例，而每個(gè)混合像元的光譜反映值則是端元光譜與豐度的線性組合[14-17]。

最小二乘約束的原則是像元內(nèi)每種端元的比例值范圍在0～1之間，并且和為1，即

∑wi=1wi∈[0,1]

(2)

基于光譜解混原理，利用光譜庫(kù)中匹配到的標(biāo)準(zhǔn)光譜去表示影像像元地物的光譜，通過(guò)統(tǒng)計(jì)分析和多元線性回歸分析，實(shí)現(xiàn)影像實(shí)際光譜到光譜庫(kù)光譜的轉(zhuǎn)換。假設(shè)ρn{ρn1,ρn2,…,ρni}為歸一化時(shí)所用目標(biāo)影像Sentinel-2A各波段各像元點(diǎn)的光譜集合，sn{s1n,s2n,…,skn}表示從光譜庫(kù)中經(jīng)匹配分析獲取的像元代表光譜在波段n處的等效值，設(shè)目標(biāo)影像像元多元線性回歸表達(dá)式為

ρni=a0i+w1is1ni+w2is2ni+…+wkiskni

(3)

式中，i為像元點(diǎn)數(shù)；n為波段數(shù)，參考影像與目標(biāo)影像的相似波段一一對(duì)應(yīng)；wki為權(quán)重，表示光譜庫(kù)中某一匹配光譜占目標(biāo)影像實(shí)際像元光譜的比重；k表示典型光譜的數(shù)量。各代表光譜經(jīng)匹配因子gkn調(diào)整之后，按像元順序組合，則得到歸一化之后各波段各像元點(diǎn)的光譜集合ρN{ρN1,ρN2,…,ρNi}，即

lkNi=gknskni

(4)

ρNi=b0i+w1il1Ni+w2nl2Ni+…+wkilkNi

(5)

3 數(shù)據(jù)結(jié)果處理與分析

3.1 典型地物解混

根據(jù)影像中地物類(lèi)別，結(jié)合幾種經(jīng)典的選取端元的方法及目視判斷，選取植被、草地、高亮白沙地、水體4種典型地物作為合適的待選端元，綜合采用最小距離及光譜角法在等效光譜庫(kù)中匹配和典型地物光譜最相似的標(biāo)準(zhǔn)地物光譜，并將其等效值作為輸入對(duì)Sentient-2A影像進(jìn)行解混，得到的豐度圖像如圖3所示。從解混結(jié)果上看，各典型地物的代表性可以較好地體現(xiàn)出來(lái)。

3.2 歸一化結(jié)果對(duì)比

基于4種典型地物的比例圖像及依據(jù)光譜庫(kù)計(jì)算的匹配因子庫(kù)，按照式(5)逐像元校正組合之后獲取歸一化之后的Sentinel-2A影像，同原始影像及采用HLS[17]線性參數(shù)歸一化后影像，分別與同時(shí)刻的Landsat 8參考影像進(jìn)行散點(diǎn)回歸分析，如圖4所示。散點(diǎn)圖中實(shí)線為線性擬合線，虛線為y=x,擬合方程、R2、RMSE統(tǒng)計(jì)見(jiàn)表2。

表2 相關(guān)性分析統(tǒng)計(jì)

從圖表中可以看出，雖然Sentinel-2A與Landsat 8 OLI的波段設(shè)置相似，選取的影像對(duì)成像時(shí)間僅相差7 min，大氣條件接近，地物類(lèi)型也相同，但是經(jīng)過(guò)各自官方大氣校正算法反演的反射率影像卻具有差異。對(duì)比原始的Sentinel-2A反射率影像，發(fā)現(xiàn)未歸一化的Sentinel-2A VNIR波段反射率比Landsat 8 OLI偏高，藍(lán)波段最為明顯，且從藍(lán)波段到近紅外波段這種偏高的趨勢(shì)依次遞減，相關(guān)性逐漸增加，各波段的R2依次為0.805 9、0.935 2、0.961 1、0.991 2。即隨著波長(zhǎng)的增加，兩傳感器影像反演得到的反射率更加相近，受大氣校正方法差異及輻射響應(yīng)差異的影響越來(lái)越小。但是從影像中各地物類(lèi)型分開(kāi)來(lái)看，并不是所有地物的各波段的Sentinel-2A反射率值均比Landsat 8 OLI高，比如植被區(qū)，在紅波段Sentinel-2A的反射率值比Landsat 8 OLI低，而其他3個(gè)波段則相反，4個(gè)波段在不同類(lèi)型的地物區(qū)呈現(xiàn)不同的高低值。

HLS計(jì)劃下的光譜通道歸一化系數(shù)是通過(guò)擬合模擬的Landsat 8和Sentinel-2地表反射率值得到的。從中間的散點(diǎn)圖可以看出，同原始影像一樣，經(jīng)HLS校正系數(shù)歸一化的影像R2隨波長(zhǎng)的增加而增大，相比于原始影像，藍(lán)波段相關(guān)性有細(xì)微的提高，紅、近紅外波段保持不變，綠波段降低。其次，就RMSE而言，藍(lán)、紅波段的RMSE相對(duì)于原始影像各減少了1.1%、2.7%，綠、近紅外波段保持不變。綜合來(lái)看，經(jīng)HLS校正系數(shù)歸一化之后，藍(lán)、紅波段效果有所改善，但是并不明顯，近紅外波段無(wú)變化，綠波段引入了額外的誤差。這可能由HLS系數(shù)針對(duì)的是相對(duì)一致的大氣校正方法，而本文采用的是各自官方的大氣校正算法所致。

逐像元?dú)w一化后的Sentinel-2A VNIR 4個(gè)波段從R2、RMSE、擬合系數(shù)3個(gè)方面看，都更接近于參考影像，輻射一致性明顯增強(qiáng)。相比原始影像R2，除紅波段稍差之外，其他3個(gè)波段均有所提高，藍(lán)波段提高的最為明顯，直接從0.805 9增至0.813 2。從RMSE分析，藍(lán)、綠、近紅波段RMSE各降低9.7%、1.5%、1.4%，近紅外保持不變。從擬合方程上看，經(jīng)逐像元?dú)w一化后各波段的slope更接近于1，藍(lán)、綠、紅、近紅外波段的slope分別由原來(lái)的1.131 1、1.182 6、1.181 5、1.055 2降低到1.040 0、1.178 9、1.149 0、1.054 1。相對(duì)于一個(gè)波段所有像元均采用同樣的歸一化校正系數(shù)的方法，本文提出的通過(guò)匹配因子調(diào)整的逐像元?dú)w一化方法可以一定程度上減少傳感器之間的差異，增強(qiáng)多源傳感器之間的輻射一致性，有利于遙感影像的綜合利用。

4 結(jié) 語(yǔ)

Sentinel-2A與Landsat 8 OLI均是新一代中高分辨率傳感器，提供全球觀測(cè)服務(wù)，波段設(shè)置相似，數(shù)據(jù)也可免費(fèi)獲取，且均具有大幅寬、高精度、高空間分辨率的特點(diǎn)，在定量遙感方面應(yīng)用廣泛。為了方便用戶更加高效合理地利用影像數(shù)據(jù)，兩者在眾多大氣校正算法的基礎(chǔ)上發(fā)展了各自的反射率反演算法，但即便是同時(shí)刻成像的影像，得到的反射率數(shù)據(jù)差異也比較大。因此，對(duì)其綜合應(yīng)用時(shí)，需要考慮兩者大氣校正算法的特點(diǎn)及光譜響應(yīng)的差異，將其歸一化至同一輻射基準(zhǔn)。

基于光譜庫(kù)的逐像元?dú)w一化方法可以有效地消除部分光譜響應(yīng)帶來(lái)的差異，相對(duì)于傳統(tǒng)簡(jiǎn)單的相對(duì)輻射歸一化所有波段均采用一個(gè)套轉(zhuǎn)換系數(shù)，本文方法考慮地物的類(lèi)型及傳感器波段設(shè)置特點(diǎn)，經(jīng)歸一化后的影像與目標(biāo)影像各波段之間的相關(guān)性有所提高，RMSE降低，輻射一致性增強(qiáng)。更重要的是，與傳統(tǒng)方法相比，該方法為多源中高分辨率遙感圖像高精度輻射歸一化提供了新思路。

然而，在試驗(yàn)過(guò)程中也發(fā)現(xiàn)，當(dāng)?shù)匚镱?lèi)型復(fù)雜時(shí)，典型地物的選擇、組合、光譜匹配精度等因素在一定程度影響歸一化效果。今后研究工作考慮在傳統(tǒng)輻射歸一化的基礎(chǔ)之上再應(yīng)用本文方法，結(jié)合兩者的優(yōu)勢(shì)，在其他多源中高分辨率反射率影像上驗(yàn)證。