基于正則化IR-MAD的GF-1影像輻射歸一化

黃莉婷，焦偉利，龍騰飛，康傳利

(1.桂林理工大學 測繪地理信息學院，廣西 桂林 541006；2.中國科學院空天信息創(chuàng)新研究院，北京 100094)

0 引言

高分一號(GF-1)衛(wèi)星開啟了我國對地觀測的新時代，該衛(wèi)星具備高空間分辨率、高時間分辨率、寬覆蓋和多源遙感數(shù)據(jù)的綜合數(shù)據(jù)獲取能力，能滿足各行業(yè)的應用需求。同時，國內(nèi)外各類遙感衛(wèi)星的種類日益增多，使得遙感數(shù)據(jù)資源呈現(xiàn)多源性，降低了單顆衛(wèi)星在進行長時間序列研究中的困難。Landsat-8陸地成像儀(operational land imager，OLI)和Sentinel-2多光譜儀(multi spectral instrument，MSI)數(shù)據(jù)已成為多源影像來源的研究熱點，其數(shù)據(jù)的免費獲取極大地推進了中分辨率陸地表面成像的虛擬星座計劃。可見，GF-1、Landsat-8 OLI和Sentinel-2 MSI這3顆衛(wèi)星數(shù)據(jù)在當前以及今后相當長時間內(nèi)，是實現(xiàn)衛(wèi)星長時間序列數(shù)據(jù)集構(gòu)建的重要互補性數(shù)據(jù)源。

地表反射率(bottom-of-atmosphere corrected reflectance，BOA)是遙感定量反演的一個關(guān)鍵參數(shù)。目前可公開獲取的Landsat-8 OLI Level 2級BOA產(chǎn)品和Sentinel-2 L2A級BOA產(chǎn)品都具有較高的精度，但是GF-1地表反射率產(chǎn)品精度問題依然沒有得到足夠重視。對GF-1影像大氣校正的研究，目前主要是用基于大氣輻射傳輸模型的絕對輻射校正方法。該方法對輸入?yún)?shù)精度要求嚴格，并且精確的大氣參數(shù)難以獲取。而采用基于像元值經(jīng)驗模型的相對輻射校正方法來進行研究的還比較少，在現(xiàn)有的輻射歸一化研究中，也多集中應用在多時序遙感影像的變化檢測上。如果缺少相關(guān)大氣參數(shù)(氣溶膠、臭氧濃度等)，又要得到較為精確的BOA，那么相對輻射校正是一種較好的方法。此外，Landsat-8和Sentinel-2地表反射率產(chǎn)品為輻射校正提供了較好的參考數(shù)據(jù)。

相對輻射歸一化方法假設(shè)同一區(qū)域的多時相或同時相影像相同波段的灰度值存在線性關(guān)系[1]，目前常用的線性相對輻射歸一化方法包括暗集-亮集法(dark set-bright set normalization，DB)[2]、偽不變特征法(pseudo-invariant feature，PIF)[3]、自動散點控制回歸算法(automatic scattergram-controlled regression，ASCR)[4]、多元變化檢測法(multivariate alteration detection，MAD)[5]和迭代加權(quán)多元變化檢測法(iteratively reweighted MAD，IR-MAD)[6]。不少學者基于這些方法進行了相關(guān)的輻射校正研究[7-13]。在這些方法的基礎(chǔ)上，都采取了一些方法對2幅影像上沒有發(fā)生變化且光譜特性穩(wěn)定的同一地物點(即PIF點)進行高質(zhì)量的提取，實現(xiàn)了多時相遙感影像的輻射歸一化。但是這些研究大多數(shù)都是針對Landsat數(shù)據(jù)的，隨著越來越多的不同傳感器數(shù)據(jù)的綜合應用，現(xiàn)有的相對輻射歸一化方法均存在一定的限制，不能直接應用于GF-1數(shù)據(jù)[14]。如DB方法要求有足夠的暗點和亮點，這在有些情況下滿足不了。PIF方法選取的PIF點代表影像的大部分區(qū)域時，才能得到有效的輻射校正結(jié)果。通過水陸集群中心確定未變化集(no changed set，NC)可以排除掉一些輻射特性發(fā)生較大改變的像元，但是NC的確定需依據(jù)經(jīng)驗確定，當研究區(qū)水體較少時，不能確定NC。雖然IR-MAD方法能夠有效地集中變化信息，減少背景信息對變化信息的影響，對影像的要求比較少，但仍然存在多個閾值需要依據(jù)經(jīng)驗確定的缺點，且在多次迭代的過程中可能會出現(xiàn)近奇異矩陣。

為使GF-1影像具有較高精度的BOA，降低大氣校正的復雜度，本文提出一種簡便實用的方法，對大氣校正后的GF-1影像，以Landsat-8和Sentinel-2的BOA影像為基準，通過相對輻射歸一化提高GF-1 BOA精度，同時也可減少GF-1和Landsat-8、Sentinel-2之間的輻射差異，使得多源影像數(shù)據(jù)可以作為互補性數(shù)據(jù)聯(lián)合使用。

1 方法

本文提出一種基于正則化IR-MAD的GF-1影像輻射歸一化方法。該方法不要求影像同時具備暗點和亮點特征，并能消除大樣本集中樣本噪音對輻射校正結(jié)果的影響。其主要技術(shù)流程為：首先根據(jù)目標影像及基準影像的紅波段和近紅外波段散點圖確定NC；然后運用正則化IR-MAD規(guī)則對NC進一步篩選，得到不變特征點；最后由不變特征點建立線性回歸模型求解系數(shù)，得到輻射歸一化方程，進行輻射歸一化(圖1)。

圖1 本文方法流程圖

該輻射歸一化方法主要包含確定NC和提取不變特征點2個關(guān)鍵步驟，下面分別對其進行介紹。

1.1 確定未變化集

與ASCR方法不同的是，本文方法不用散點圖上的水陸集群中心確定NC，而是由正交回歸法確定初始回歸基線，再求解NC。首先，由目標影像及基準影像的紅波段和近紅外波段分別生成散點圖，散點圖橫坐標是目標影像，縱坐標是基準影像，如圖2所示。由線性變換確定的NC能初步表示“偽變化”的像元，即輻射特性幾乎沒有改變的像元。

圖2 從散點圖確定NC

圖中HVW是NC選擇的限制條件[4]，由HPW控制，HPW是不變區(qū)域一側(cè)的垂直寬度。HVW和HPW的關(guān)系如式(1)所示。

(1)

式中：a是由正交回歸法確定的初始回歸線斜率。分別確定了紅波段和近紅外波段散點圖的初始回歸線及HVW后，由公式(2)確定最終的NC[4]。

NC=(X，Y)：|y1-b1-a1x1|≤
HVWNC1&|y2-b2-a2x2|≤HVWNC2

(2)

式中：y1、y2分別是基準影像紅波段和近紅外波段的BOA值；x1、x2分別是目標影像紅波段和近紅外波段的BOA值。

1.2 提取不變特征點

在NC中，對不變特征點的選擇做進一步的優(yōu)化。采用以下規(guī)則來尋找可以使用的擬合樣本點數(shù)據(jù)[6]：對影像的每個像元設(shè)置初始權(quán)值為1，每一次迭代過程中計算均值向量與方差矩陣，用典型相關(guān)分析的方法計算每個MAD變量，根據(jù)新計算的MAD變量更新權(quán)值，值域區(qū)間為[0，1]；未發(fā)生變化的像元具有較大的權(quán)值，經(jīng)過若干次迭代計算之后，每個像元的權(quán)值會趨于穩(wěn)定，通過權(quán)值與閾值的比較，便可判定每個像元點是否屬于不變特征點。更新權(quán)值的公式如式(3)～式(5)所示。

(3)

(4)

Pr(nochange)=1-Pχ2；N(Z)

(5)

式中：γi是第i波段的方差；Z是迭代賦值引入的隨機變量，表示標準化MAD變量平方和，其服從自由度為N(波段數(shù))的卡方檢驗，通過卡方分布概率密度函數(shù)對Z進行加權(quán)；Pr(nochange)表示迭代權(quán)值，可以設(shè)定一個固定的閾值t(如t=0.95，便能得到很好的結(jié)果)，當Pr(nochange)>t時，判斷為不變特征點。

在多次迭代的過程中可能會出現(xiàn)近奇異矩陣，當矩陣中某一元素發(fā)生很小的變動就會產(chǎn)生很大的誤差。所以，Nielsen[15]在迭代的過程中引入正則化(也稱為“懲罰項”)，避免在計算中產(chǎn)生較大權(quán)重，影響不變特征點的判斷。本文在求解協(xié)方差矩陣時，也加入了正則化參數(shù)，如式(6)～式(7)所示。

Var11=∑11+λ1Ω1

(6)

Var22=∑22+λ2Ω2

(7)

式中：λi是規(guī)則化參數(shù)；Ωi是N×N的對角線矩陣。

2 數(shù)據(jù)

2.1 實驗數(shù)據(jù)

Sentinel-2A MSI其L2A級產(chǎn)品為經(jīng)過輻射定標、大氣校正的BOA反射率影像，可由歐空局(European Space Agency，ESA)發(fā)布的專門生產(chǎn)L2A級數(shù)據(jù)的插件Sen2cor根據(jù)需求自行生產(chǎn)[16]，或者在ESA官網(wǎng)(https：//scihub.copernicus.eu/dhus/#/home)進行下載，目前可下載到2018年12月后的中國區(qū)域L2A級數(shù)據(jù)。

Landsat-8 OLI的Level 2級BOA產(chǎn)品，由美國地質(zhì)調(diào)查局(United States Geological Survey，USGS)專門為Landsat-8數(shù)據(jù)設(shè)計的大氣校正程序LaSRC(landsat-8 surface reflectance code)生成。該程序利用海岸氣溶膠波段(band1)進行氣溶膠反演測試，運用MODIS輔助氣候數(shù)據(jù)并使用專門的輻射傳輸模型，被認為是最精確的Landsat-8 OLI大氣校正程序，其結(jié)果精度被廣泛認可和使用[17]。Level 2級數(shù)據(jù)可在USGS的網(wǎng)站(https：//earthexplorer.usgs.gov/)下載。

為保證傳感器在獲取遙感數(shù)據(jù)時受大氣影響的差異較小，應該盡量選取同一日期的太陽高度角和方位角接近的研究區(qū)，進行多源影像的相對輻射歸一化[18]。但是，由于不同衛(wèi)星的重訪周期和幅寬不相同，因此難以獲得同一地區(qū)的同一日期、同時過境的不同衛(wèi)星影像。本研究的影像選取遵循影像前后日期在7 d左右的規(guī)則。

本文分別選取分辨率相近的GF-1 PMS2和Sentinel-2A MSI，以及GF-1 WFV4和Landsat-8 OLI影像進行2組實驗。第1組選擇了2017年12月19日的GF-1 PMS2影像和2017年12月11日Sentinel-2A MSI影像，土地覆被類型以丘陵、水體為主；第2組選擇了2018年6月14日的GF-1 WFV4影像和2018年6月11日Landsat-8 OLI影像，土地覆被類型以耕地、裸土和不透水面(建筑物、道路等)為主。

3顆衛(wèi)星影像的對應波段信息對比如表1及圖3所示。由表1及圖3可知，對于GF-1 PMS2和Sentinel-2A MSI而言，除了藍波段，GF-1 PMS2的波譜帶寬明顯比Sentinel-2A寬，但中心波長差距較??；GF-1 WFV4和Landsat-8 OLI 2類傳感器在可見光波段波譜帶寬設(shè)置較相近，在近紅外波段波譜帶寬設(shè)置差異顯著，GF-1在近紅外波段處有較寬的波譜通道，并且在近紅外波段二者中心波長相距較大。

表1 所選影像波段信息

圖3 傳感器光譜響應曲線

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

獲取BOA需要對GF-1的PMS2和WFV4影像進行絕對輻射校正，在ENVI5.3中進行。首先把影像灰度值(digital number，DN)值轉(zhuǎn)換為輻射亮度，再將輻射亮度轉(zhuǎn)換成BOA。在本文實驗中具體處理步驟包括：輻射定標、大氣校正、正射校正和影像重采樣。Landsat-8 OLI和Sentinel-2A MSI影像經(jīng)過波段組合后，需要和處理過的GF-1 BOA影像進行精確的幾何配準，配準前GF-1需要重采樣，配準精度要求達到均方根誤差小于0.5個像元，影像配準的質(zhì)量會影響輻射歸一化的效果。

3 實驗與分析

3.1 不變特征點提取

相對輻射歸一化主要由不變特征點提取、求解回歸系數(shù)和輻射變換3個步驟組成，在本文中選取影像大小為1 000像素×1 000像素進行實驗。為減弱大氣輻射的影響，首先將目標影像和基準影像的近紅外和紅波段分別生成散點密度圖，確定初始回歸線；再由初始回歸線及式(1)、式(2)確定每個波段的NC。初始回歸線及NC如圖4(c2)、圖4(d2)、圖5(c2)和圖5(d2)所示。接著運用正則化IR-MAD規(guī)則對NC做進一步篩選。NC中每個像元根據(jù)自由度N為4(波段數(shù))以及影像差值的卡方變量Z，通過式(5)得到每個像元的權(quán)值。再將權(quán)值與閾值t進行比較，判斷該像元是否屬于變化和未變化像元，一般閾值t需要根據(jù)經(jīng)驗多次求取。為分析閾值t對求解輻射歸一化方程的影響，本文在t∈[0.91，0.95]等間隔取5個閾值：0.91、0.92、0.93、0.94、0.95，分別計算各波段的回歸判定系數(shù)R2(式(8))和均方根誤差RMSE(式(9))。

圖4 GF-1 PMS2和Sentinel-2A MSI 不變特征點的3次選擇過程

(8)

(9)

R2和RMSE計算結(jié)果如表2所示。分析表2可知，2組實驗中各波段在不同閾值t的情況下R2和RMSE均相差不大，但RMSE值會隨著閾值t的升高略有降低。分析其原因，是因為隨著閾值t的升高，樣本數(shù)量逐漸降低，導致RMSE值也逐漸減小。再分別對5個閾值輻射歸一化后的影像以相同的樣本數(shù)量計算R2和RMSE(表3)。由表3可知，2組實驗中各波段的RMSE的差距進一步減小，由于使用了相同的樣本集，所以各波段的R2相一致。該結(jié)果說明，所有閾值都能減小目標影像與基準影像之間的輻射差異，也進一步說明了該方法的有效性和穩(wěn)定性。

表3 相同樣本集不同閾值t的R2和RMSE

由于閾值t的穩(wěn)定性，2組實驗均取相同的閾值t=0.95進行下一步分析。

圖4和圖5分別顯示了2組實驗影像各波段不變特征點的3次選擇過程，圖中綠色直線代表初始回歸線，紅色直線為最終求解出來的擬合線。各實驗影像的BOA反射率值均被擴大了10 000倍。

圖5 GF-1 WFV4和Landsat-8 OLI 不變特征點的3次選擇過程

第1次圖(a1)、圖(b1)、圖(c1)、圖(d1)是目標影像和基準影像對應波段生成的散點圖；第2次圖(a2)、圖(b2)、圖(c2)、圖(d2)顯示的是由式(1)和式(2)確定的NC；第3次圖(a3)、圖(b3)、圖(c3)、圖(d3)是由正則化的迭代多元變化檢測規(guī)則提取出來的不變特征點。這些散點圖有一些共同的特征：不變特征點逐漸收斂，最后一次散點圖中的點數(shù)量減少，沿直線分布；主軸的斜率從第2次到第3次有一定的變化；散點的變化范圍逐漸變小，且各點分布較均勻。這些特征表明了該方法在不變特征點選擇中的有效性和穩(wěn)定性。

傳感器波段t0.91t0.92t0.93t0.94t0.95R2RMSER2RMSER2RMSER2RMSER2RMSE藍 0.9807.9160.9807.6480.9827.3380.9827.3490.9846.839PMS2和MSI綠 0.98610.7840.98810.3420.98810.1350.9889.3680.9908.736紅 0.97610.9940.97810.6120.97810.7460.97810.0840.9809.832近紅外0.91525.8500.93524.4750.99424.8560.92724.5230.92723.759藍 0.96614.3880.98314.8600.97013.6920.96814.0210.97014.071WFV4和OLI綠 0.97417.5430.98717.2460.97616.7540.97415.8920.97615.743紅 0.98422.7090.99122.6460.98421.5470.98420.9390.98419.911近紅外0.96042.2900.97842.6770.96641.2950.96640.7430.96438.917

3.2 求解輻射歸一化方程和輻射變換

將上述選出的不變特征點集運用正交回歸法分別計算得到PMS2和MSI以及WFV4和OLI在藍、綠、紅、近紅外波段的反射率歸一化方程(表4)。最后把各波段輻射歸一化方程應用在目標影像GF-1 PMS2和GF-1 WFV4相應波段中，得到輻射變換后的影像。為了提供能夠衡量總體差異性的指標，1/3不變特征點集作為驗證數(shù)據(jù)用來進行光譜輻射歸一化結(jié)果的驗證。根據(jù)式(8)計算回歸判定系數(shù)R2，式(10)計算相關(guān)系數(shù)r，式(11)計算平均值差異Δmean，計算結(jié)果如表4所示。

表4 輻射歸一化方程與R2、r、Δmean

(10)

(11)

在2組實驗中，各波段的R2均在0.92以上，說明所建立的回歸模型擬合效果較好。Δmean都較小，說明輻射變換后的GF影像和基準影像差異程度較小。

結(jié)合圖4、圖5(a1)、圖5(b1)、圖5(c1)、圖5(d1)分析，GF-1 PMS2和Sentinel-2A MSI以及GF-1 WFV4和Landsat-8 OLI的BOA值在可見光和近紅外波段都不同程度地存在一定偏差。對比表4中各傳感器對應波段的光譜值的相關(guān)系數(shù)r，均大于0.96，這表明對應波段反射率都具有很強的線性相關(guān)性，基準影像能夠顯著提升GF-1衛(wèi)星影像的輻射精度。

3.3 結(jié)果分析

1)目視效果分析。采用目視的方法對輻射歸一化后的影像進行定性評價，將歸一化后的影像與基準影像進行顯示對比，如果2個影像顏色、亮度等十分相似，則說明歸一化的效果較好[19]。圖6(a)和圖6(b)是GF-1 PMS2和Sentinel-2A MSI輻射歸一化前后結(jié)果的對比；圖6(c)和圖6(d)是GF-1 WFV4和Landsat-8 OLI輻射歸一化前后結(jié)果的對比，影像均在ENVI中選擇相同的拉伸方式顯示。觀察圖6(a)和圖6(b)，在PMS2(右)和MSI(左)影像接邊處，水體、植被、耕地等土地覆被類型的邊界鑲嵌處差異明顯減小，2個影像各土地覆被類型的顏色、亮度等十分相似。在圖6(c)和圖6(d)中，WFV4(右)和OLI(左)影像接邊處裸地、植被、耕地等土地覆被類型的邊界鑲嵌處差異較小，但在城市建筑用地區(qū)域，2個影像的接邊處還存在著一定亮度的差異。整體來看，在2組實驗中，目標影像與基準影像邊界鑲嵌處差異減小，可見數(shù)據(jù)輻射一致性明顯改善。經(jīng)過輻射歸一化后的影像，BOA值的分布會有所改變，從輻射歸一化前后4個波段的核密度估計(kernel density estimation，KDE)密度分布圖(圖7、圖8)可以看出，輻射歸一化后的影像BOA數(shù)值分布，更接近于基準影像BOA數(shù)值的分布特征。其中，圖7(d)有2個峰谷，分別代表的是影像上水、陸集群中心。圖8(d)中由于影像上沒有水像元，所以只有一個峰谷。本文方法不僅適用于有水陸分布的影像，對于缺少水體分布的影像也能很好地進行輻射歸一化。

圖6 輻射歸一化前后對比

圖7 PMS2影像、MSI影像和歸一化后影像BOA密度分布比較

圖8 WFV4影像、OLI影像和歸一化后影像BOA密度分布比較

2)典型地物波譜曲線分析。因受輻射畸變和不同傳感器的特性影響，導致地物波譜曲線有很大的差異，相對輻射歸一化的目的就是消除或減弱這些差異[20-21]。地物的波譜曲線越相似，說明地物的反射特性越相近，地物波譜曲線的變化在一定程度上可以說明相對輻射歸一化的效果[20]。本文分別從目標影像、輻射歸一化后影像與基準影像中選取幾種典型地物進行波譜曲線分析。從圖9和圖10中可以看出，經(jīng)輻射歸一化后的影像，其典型地物波譜曲線向基準影像拉近，說明輻射歸一化是有效果的。結(jié)合圖9、圖10及表4分析，在可見光波段(藍、綠、紅)對地物光譜特征的輻射歸一化效果相對較好，而近紅外波段輻射歸一化效果不如可見光波段。分析其原因，主要是由于PMS2和WFV4相較于MSI和OLI傳感器，在近紅外波段的光譜響應函數(shù)差異明顯大于其他3個可見光波段的差異，尤其是光譜范圍在770～900 nm之間的光譜響應曲線(圖3)。

3)與其他方法的比較。將本文方法和原ASCR方法、IR-MAD方法的歸一化效果進行比較，以均方根誤差RMSE作為歸一化影像與基準影像的相似性度量指標，RMSE的值越小，說明輻射歸一化效果越好。

以圖4和圖5最后一次篩選出來的不變特征點作為樣本集，分別計算目標影像和基準影像的RMSE，以及ASCR方法、IR-MAD方法和本文方法輻射歸一化前后影像的RMSE，計算結(jié)果如表5所示。由表5可知，經(jīng)輻射歸一化后的影像RMSE均小于未輻射校正前的，說明這3種方法都減小了目標影像和基準影像的輻射值差異，3種方法同樣在近紅外波段輻射歸一化效果不及可見光波段，其中本文方法與ASCR方法、IR-MAD方法相比，整體精度都較高。

表5 不同方法輻射歸一化的RMSE

ASCR方法采用整景影像比較多的像元信息進行輻射歸一化，因此易受到離初始回歸線較遠的像元影響。那些離初始回歸線較遠的像元不僅線性作用較差，反射率還發(fā)生了變化，會給輻射歸一化效果帶來一定的影響。IR-MAD方法基于輻射特征的變化進行不變特征點的選取，其MAD變量中的像元值反映了變化信息的情況，對未變化像元的選取更加準確。而在本文方法中，用初始回歸線先確定了NC，保證未變化像元的線性關(guān)系，再對NC采用IR-MAD的規(guī)則進一步篩選，再次控制了不變特征點的選取質(zhì)量，所建立的輻射歸一化方程的轉(zhuǎn)換效果較好。

4 結(jié)束語

在GF-1大氣校正后，本文基于IR-MAD，采用控制NC的線性關(guān)系、逐步篩選不變特征點的方法進行GF-1影像與其他影像的光譜輻射歸一化。實驗表明，GF-1 PMS2和Sentinel-2A MSI以及GF-1 WFV4和Landsat-8 OLI在BOA上有較好的線性關(guān)系，不變特征點的相關(guān)系數(shù)達到了0.96以上，說明通過輻射歸一化提高GF-1影像的輻射精度是可行的；與其他輻射歸一化方法相比，本文方法具有較好的精度。分別從目視效果、典型地物波譜曲線、均方根誤差等方面分析評價了輻射歸一化結(jié)果，取得了較為理想的效果，為提高GF-1影像BOA的輻射精度提供了一種解決思路。

本文主要針對大氣校正后的GF-1影像進行研究。由于地表具有二向反射分布特性，反射率隨太陽入射角及觀測角的變化而不同，導致不同時期和不同傳感器獲取的影像其成像幾何角度有所差異，所以應當進行BRDF校正以減小由成像幾何角度差異引起的反射率差異[22-23]，這在今后的研究中需要進一步分析。