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

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

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

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

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

表1 Sentitle-2A與Landsat 8 OLI可見—近紅外波段對比

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

2 逐像元歸一化方法

2.1 光譜匹配因子庫建立

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

(1)

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

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

2.2 逐像元歸一化模型

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

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

∑wi=1wi∈[0,1]

(2)

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

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

(3)

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

lkNi=gknskni

(4)

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

(5)

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

3.1 典型地物解混

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

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

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

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

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

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

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

4 結(jié) 語

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

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

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