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

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

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

0 引言

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

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

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

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

1 方法

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

圖1 本文方法流程圖

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

1.1 確定未變化集

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

圖2 從散點圖確定NC

圖中HVW是NC選擇的限制條件[4]，由HPW控制，HPW是不變區域一側的垂直寬度。HVW和HPW的關系如式(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中，對不變特征點的選擇做進一步的優化。采用以下規則來尋找可以使用的擬合樣本點數據[6]：對影像的每個像元設置初始權值為1，每一次迭代過程中計算均值向量與方差矩陣，用典型相關分析的方法計算每個MAD變量，根據新計算的MAD變量更新權值，值域區間為[0，1]；未發生變化的像元具有較大的權值，經過若干次迭代計算之后，每個像元的權值會趨于穩定，通過權值與閾值的比較，便可判定每個像元點是否屬于不變特征點。更新權值的公式如式(3)～式(5)所示。

(3)

(4)

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

(5)

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

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

Var11=∑11+λ1Ω1

(6)

Var22=∑22+λ2Ω2

(7)

式中：λi是規則化參數；Ωi是N×N的對角線矩陣。

2 數據

2.1 實驗數據

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

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

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

本文分別選取分辨率相近的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顆衛星影像的對應波段信息對比如表1及圖3所示。由表1及圖3可知，對于GF-1 PMS2和Sentinel-2A MSI而言，除了藍波段，GF-1 PMS2的波譜帶寬明顯比Sentinel-2A寬，但中心波長差距較小；GF-1 WFV4和Landsat-8 OLI 2類傳感器在可見光波段波譜帶寬設置較相近，在近紅外波段波譜帶寬設置差異顯著，GF-1在近紅外波段處有較寬的波譜通道，并且在近紅外波段二者中心波長相距較大。

表1 所選影像波段信息

圖3 傳感器光譜響應曲線

2.2 數據處理

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

3 實驗與分析

3.1 不變特征點提取

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

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

(8)

(9)

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

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

由于閾值t的穩定性，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)是由正則化的迭代多元變化檢測規則提取出來的不變特征點。這些散點圖有一些共同的特征：不變特征點逐漸收斂，最后一次散點圖中的點數量減少，沿直線分布；主軸的斜率從第2次到第3次有一定的變化；散點的變化范圍逐漸變小，且各點分布較均勻。這些特征表明了該方法在不變特征點選擇中的有效性和穩定性。

傳感器波段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不變特征點集作為驗證數據用來進行光譜輻射歸一化結果的驗證。根據式(8)計算回歸判定系數R2，式(10)計算相關系數r，式(11)計算平均值差異Δmean，計算結果如表4所示。

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

(10)

(11)

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

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

3.3 結果分析

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

圖6 輻射歸一化前后對比

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

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

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

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

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

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

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

4 結束語

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

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