王藍(lán)星，王群明,2，童小華

1.同濟(jì)大學(xué)測繪與地理信息學(xué)院，上海 200092; 2.上海市數(shù)字光學(xué)前沿科學(xué)研究基地，上海 200092

遙感影像對于全球地表監(jiān)測具有重要的作用。其中，高光譜影像能提供地物豐富的光譜信息，被廣泛應(yīng)用于環(huán)境監(jiān)測、精細(xì)農(nóng)業(yè)、異常目標(biāo)檢測、海洋監(jiān)測及地表礦物制圖等領(lǐng)域[1-5]，具有不可替代的地位。GF-5衛(wèi)星于2018年5月發(fā)射，其作為全譜段的高光譜衛(wèi)星，首次同時(shí)實(shí)現(xiàn)了對大氣和陸地的綜合觀測，是中國高分專項(xiàng)中一顆重要的科研衛(wèi)星，具有重大的應(yīng)用價(jià)值。然而，云層(尤其是厚云)遮擋對高光譜影像的應(yīng)用造成了較大的影響。研究高光譜影像的云去除問題具有重要意義。

現(xiàn)有的厚云去除(以下簡稱為“云去除”)方法大多針對于多光譜影像，主要分為3類：①基于空間信息的方法。此類方法僅利用云影像上非云區(qū)域中的有效信息(一般為云覆蓋區(qū)域的鄰近有效區(qū)域)對云覆蓋區(qū)域進(jìn)行重建[6-7]。文獻(xiàn)[8]通過分析影像結(jié)構(gòu)和紋理的空間模式，計(jì)算待重建影像上已知的最相似圖塊對之間的偏移量，在能量函數(shù)最小化的前提下，通過復(fù)制已知像元的方式填補(bǔ)缺失值。然而，由于此類方法利用的已知信息較為有限，預(yù)測結(jié)果的不確定性較大(尤其在云覆蓋區(qū)域較大且紋理較為復(fù)雜時(shí))。因此，此類方法更適用于小范圍的云層去除。②基于時(shí)間信息的方法[9-11]。此類方法利用同一地區(qū)其他時(shí)間點(diǎn)獲取的無云影像(即輔助影像)提供輔助信息，其關(guān)鍵在于利用輔助影像和云影像中共同無云區(qū)域的有效信息對兩者之間的關(guān)系進(jìn)行準(zhǔn)確建模。然而，由于不同的地物往往隨時(shí)間發(fā)生不同程度的變化，輔助影像與目標(biāo)云影像之間的關(guān)系會(huì)隨時(shí)間間隔增大變得愈加復(fù)雜，這給兩者之間關(guān)系模型的構(gòu)建造成了困難。因此，與目標(biāo)云影像在時(shí)間上最接近的無云影像是輔助影像的最佳選擇。能否獲得時(shí)間上與目標(biāo)云影像足夠接近的輔助影像是此類方法能否獲得令人滿意的重建結(jié)果的關(guān)鍵因素。全局直方圖匹配(global linear histogram match，GLHM)方法[12]基于輔助影像與云影像之間共同的無云區(qū)域數(shù)據(jù)構(gòu)建兩者之間的線性關(guān)系，從而對剩余缺失區(qū)域進(jìn)行預(yù)測。GLHM方法最初為解決Landsat 7 ETM+條帶問題而提出，但對于云去除問題同樣適用，該方法簡單易實(shí)現(xiàn)，能夠快速對大范圍缺失數(shù)據(jù)進(jìn)行重建。③基于時(shí)空信息的方法。此類方法綜合利用時(shí)間和空間信息，相較于前兩種方法，具有更好的穩(wěn)定性，能夠應(yīng)對不同情況下云下缺失信息的重建任務(wù)[13-16]。例如，增強(qiáng)的近鄰相似像元插值法(modified neighborhood similar pixel interpolator,MNSPI)[17]在預(yù)測每個(gè)云像元時(shí)，利用云影像本身的無云信息獲得一個(gè)空間預(yù)測項(xiàng)，同時(shí)利用輔助影像上對應(yīng)位置的有效信息估計(jì)一個(gè)時(shí)間預(yù)測項(xiàng)，最后根據(jù)影像中的景觀同質(zhì)程度和像元信息隨時(shí)間變化程度為空間項(xiàng)和時(shí)間項(xiàng)賦予不同的權(quán)重，得到最終結(jié)果。MNSPI方法計(jì)算效率高且具有較好的準(zhǔn)確性，是一種常用的云去除方法。此外，基于機(jī)器學(xué)習(xí)(如深度學(xué)習(xí))的方法也得到了越來越廣泛的應(yīng)用[18-19]，但其一般對訓(xùn)練樣本的數(shù)量有較高的要求，如通常需要大量待預(yù)測區(qū)域之外的數(shù)據(jù)用于模型訓(xùn)練。文獻(xiàn)[20]提出一種自適應(yīng)回歸方法，綜合利用長時(shí)間序列數(shù)據(jù)(部分影像包含局部云污染)進(jìn)行云去除。文獻(xiàn)[21]基于光譜時(shí)間度量，使用K-近鄰回歸方法利用一年的無云影像預(yù)測缺失值。

與上述針對多光譜影像云去除的研究相比，高光譜影像的云去除面臨著更多的挑戰(zhàn)，其研究整體上遠(yuǎn)落后于多元譜影像云去除。分析上述云去除方法可以發(fā)現(xiàn)，輔助影像的利用是云去除的必要手段。然而，現(xiàn)有的高光譜數(shù)據(jù)時(shí)間分辨率往往較低。具體來說，GF-5衛(wèi)星的重訪周期為51 d[22]，EO-1高光譜數(shù)據(jù)的時(shí)間分辨率被設(shè)計(jì)為200 d[23]。因此，時(shí)域上最鄰近的同源輔助影像與目標(biāo)高光譜影像之間通常存在著較大的時(shí)間間隔，其間地物覆蓋可能發(fā)生了較大的變化(例如植被的自然枯萎，農(nóng)作物的輪作，城區(qū)的擴(kuò)張等)，使得輔助影像的參考價(jià)值大大降低。因而，對于高光譜影像(如GF-5和EO-1高光譜影像)的云去除，尋求其他具有更高時(shí)間分辨率的輔助影像(如多光譜影像)十分必要。

Landsat系列數(shù)據(jù)是至今應(yīng)用最廣泛的多光譜數(shù)據(jù)之一[24-30]。目前，Landsat 8衛(wèi)星仍在穩(wěn)定運(yùn)行，其全球周期性覆蓋，較高的時(shí)間分辨率(16 d)[31]使得其更有可能提供時(shí)間上更接近于目標(biāo)高光譜云影像的輔助影像。同時(shí)，其30 m空間分辨率和GF-5及EO-1高光譜影像一致。此外，Landsat 8衛(wèi)星數(shù)據(jù)可供用戶免費(fèi)下載，是一種易獲取的輔助影像。然而，現(xiàn)有云去除方法通?；谳o助影像與云影像波段一一區(qū)間對應(yīng)的假設(shè)，即對任一含云波段，選取輔助影像中同區(qū)間的波段進(jìn)行重建。但是，由于光譜分辨率不一致，多光譜影像與高光譜影像各波段的光譜區(qū)間(即光譜分辨率)設(shè)置存在較大的差異?，F(xiàn)有方法無法用于基于多光譜輔助影像的高光譜云去除。為此，必須尋求更有效的波段映射模型以充分利用多光譜數(shù)據(jù)進(jìn)行高光譜數(shù)據(jù)的云去除。

文獻(xiàn)[32]提出一種基于空譜信息的隨機(jī)森林(spatial-spectral-based random forest,SSRF)方法，利用同源多光譜輔助影像對多光譜(即Landsat和Sentinel-2)遙感影像進(jìn)行云去除。SSRF方法繼承了經(jīng)典的隨機(jī)森林(RF)方法[33]在描述因變量與自變量之間復(fù)雜非線性關(guān)系時(shí)的優(yōu)勢，其在對多光譜云覆蓋影像的各個(gè)波段分別進(jìn)行云去除時(shí)，可同時(shí)利用時(shí)域近鄰影像的多個(gè)波段構(gòu)建擬合模型，無須假設(shè)輔助影像與云影像波段一一對應(yīng)。這些特點(diǎn)為異源數(shù)據(jù)的利用提供了可能。本文基于已有的SSRF方法，研究將其擴(kuò)展至異源影像間的云去除情形，提出一種利用Landsat 8 OLI影像作為輔助影像對高光譜影像(如GF-5和EO-1高光譜影像)進(jìn)行云去除的方法(記為SSRF_M方法)，以解決同源高光譜輔助影像因時(shí)間間隔較長而存在較大地物變化的難題。

1 本文方法

SSRF_M方法同時(shí)利用Landsat 8 OLI數(shù)據(jù)的6個(gè)波段(blue、green、red、NIR、SWIR 1和SWIR 2波段)對GF-5或EO-1高光譜影像的每個(gè)云覆蓋波段單獨(dú)進(jìn)行重建。具體來說，SSRF_M方法在構(gòu)建樣本時(shí)，對于輔助數(shù)據(jù)中的有效像元，利用以其為中心的3×3像元大小的圖像塊在6個(gè)波段的Landsat反射率值作為樣本的輸入(即自變量)，對應(yīng)輸出(即因變量)為待預(yù)測數(shù)據(jù)中同一位置下的有效高光譜像元的反射率值。本文對每個(gè)波段缺失像元的預(yù)測如下

(1)

圖1 SSRF_M流程(以高光譜單波段為例)

SSRF_M的實(shí)施主要包括兩個(gè)階段：SSRF模型訓(xùn)練和云下信息重建(即預(yù)測階段)。其對高光譜影像中每個(gè)云覆蓋波段單獨(dú)處理，具體如下。

1.1 訓(xùn)練階段

1.1.1 構(gòu)建訓(xùn)練樣本

(1)獲取待重建高光譜影像的云掩膜，用于區(qū)分云覆蓋區(qū)域和無云區(qū)域。原始云掩膜中通常賦予云及其陰影不同取值以便區(qū)分，在本文中將云和陰影均視為缺失區(qū)域生成新的掩膜進(jìn)行預(yù)測。

(2)獲取與待重建高光譜影像在時(shí)間上最接近的覆蓋同一地區(qū)的無云Landsat 8 OLI影像，作為輔助影像。

(3)利用云掩膜提取輔助影像與云影像上共同的無云區(qū)域數(shù)據(jù)。同時(shí)，利用云掩膜在輔助影像上提取云影像中云覆蓋區(qū)域?qū)?yīng)位置的有效數(shù)據(jù)。

(4)利用上一步中的無云區(qū)域數(shù)據(jù)構(gòu)建訓(xùn)練數(shù)據(jù)。具體來說，自變量(即輸入)為輔助的Landsat數(shù)據(jù)中的圖像塊，因變量(即輸出)為目標(biāo)云覆蓋高光譜波段中對應(yīng)的無云數(shù)據(jù)。根據(jù)常用經(jīng)驗(yàn)值，本文隨機(jī)選擇30%的樣本進(jìn)行訓(xùn)練，剩下的70%的樣本用于評價(jià)訓(xùn)練模型的準(zhǔn)確性。同時(shí)，將樹的數(shù)量設(shè)置為100。

1.1.2 SSRF模型訓(xùn)練

將得到的訓(xùn)練數(shù)據(jù)輸入SSRF，獲得自變量(即6個(gè)Landsat波段下的3×3像元圖像塊)和因變量(即高光譜單個(gè)波段中的缺失像元)之間的非線性關(guān)系模型(即式(1)中的f)。本文基于R軟件運(yùn)行SSRF程序。

1.2 預(yù)測階段

1.2.1 構(gòu)建預(yù)測樣本

利用構(gòu)建訓(xùn)練樣本階段中步驟(3)中得到的云覆蓋區(qū)對應(yīng)Landsat輔助影像上的有效數(shù)據(jù)，構(gòu)建測試數(shù)據(jù)(即圖塊自變量)，作為預(yù)測數(shù)據(jù)的輸入。

1.2.2 SSRF模型預(yù)測

將上一步中的測試數(shù)據(jù)輸入通過步驟1.2訓(xùn)練好的SSRF模型，得到云像元的預(yù)測數(shù)據(jù)。以上步驟依次用于高光譜影像各波段進(jìn)行云去除。

2 試 驗(yàn)

2.1 試驗(yàn)數(shù)據(jù)

本文利用GF-5和EO-1高光譜影像，選取3個(gè)不同的區(qū)域進(jìn)行模擬的厚云去除試驗(yàn)。每個(gè)區(qū)域均包含一景用于生成模擬云影像的無云高光譜影像，還包括用于提供時(shí)間輔助信息的一景同源高光譜影像以及一景異源(Landsat 8 OLI)多光譜影像。試驗(yàn)中使用的3類數(shù)據(jù)的空間分辨率均為30 m(圖2)。使用數(shù)據(jù)的具體信息見表1。如圖2所示，對于每個(gè)區(qū)域，本文基于無云參考影像生成模擬云影像。特別地，區(qū)域3的云掩膜基于其他區(qū)域真實(shí)的云生成，最大限度地接近于真實(shí)的云覆蓋情況，而其余兩個(gè)區(qū)域中，云均為人工隨機(jī)生成。各區(qū)域中高光譜影像和多光譜影像顯示時(shí)RGB波段均一一對應(yīng)。區(qū)域1和區(qū)域2均選自中國北京，為典型的城市區(qū)域，異質(zhì)性較強(qiáng)，紋理復(fù)雜。其中，區(qū)域1中同源輔助影像和異源輔助影像均與模擬的云覆蓋影像較為相似，其原因可能在于該地區(qū)種植的植被四季常青。然而，在區(qū)域2中，同源的輔助影像與云覆蓋影像之間存在明顯的色調(diào)差異，其原因在于兩者之間時(shí)間間隔較大，地物發(fā)生了一定程度的季節(jié)性變化。相比之下，與云覆蓋影像之間時(shí)間間隔更小的異源輔助影像與前者在視覺上十分接近。區(qū)域3選自美國法戈，其作為典型的農(nóng)業(yè)區(qū)域，勻質(zhì)性較好，然而對比同源輔助影像和云覆蓋影像看可以看出，該區(qū)域地物發(fā)生了顯著的變化，整體色彩差異較大，給同源高光譜輔助信息的利用造成了一定的困難。反之，Landsat 8 OLI影像的獲取時(shí)間更接近于云覆蓋影像，二者在視覺上更為接近。試驗(yàn)中使用的高光譜影像均剔除了質(zhì)量較差的波段。此后，本文分別利用與參考影像同源的無云影像以及Landsat 8 OLI無云影像作為輔助影像，重建模擬云影像的云下缺失信息。

注：3個(gè)區(qū)域高光譜影像分別以波段150、10、3，163、111、65和161、104、58為RGB組合；區(qū)域1的Landsat影像以red、blue、green，區(qū)域2和區(qū)域3的Landsat影像均以SWIR1、NIR、blue為RGB組合。

表1 厚云去除試驗(yàn)數(shù)據(jù)

2.2 試驗(yàn)內(nèi)容

為論證SSRF_M的有效性，本文將其與經(jīng)典的MNSPI方法進(jìn)行了對比。MNSPI方法作為一種典型的主流方法，利用同源數(shù)據(jù)提供輔助信息，要求輔助影像與云影像波段之間存在一一對應(yīng)的關(guān)系，不適用于利用異源輔助影像的情況?？紤]到MNSPI方法僅能利用單個(gè)波段的輔助數(shù)據(jù)參與對應(yīng)波段的云去除，為進(jìn)一步在方法層面驗(yàn)證SSRF_M的有效性，本文測試了同樣計(jì)算簡單，且能同時(shí)利用多個(gè)波段進(jìn)行預(yù)測的GLHM方法。為使GLHM方法在數(shù)據(jù)利用層面上與SSRF_M更接近，同樣基于3×3局部圖像塊進(jìn)行運(yùn)算，故將其稱為利用多光譜影像的空譜GLHM(spatial-spectral GLHM，SSGLHM_M)方法。總之，對于每個(gè)區(qū)域的模擬云影像，同時(shí)測試了3種云去除的方法：①SSRF_M：基于SSRF方法，利用Landsat 8 OLI作為輔助影像；②MNSPI_H：基于MNSPI方法，利用與云影像同源的其他時(shí)間點(diǎn)獲取的無云高光譜影像作為輔助影像。③SSGLHM_M：基于GLHM方法，利用Landsat 8 OLI輔助影像所提供的空間和光譜信息。需要說明的是，本文提出的SSRF_M僅需利用單景輔助影像，具有簡單易實(shí)現(xiàn)的優(yōu)勢，故并未將其與需要大量輔助影像或訓(xùn)練數(shù)據(jù)(例如需多景或一年內(nèi)的所有可用輔助影像)的方法[18-21]進(jìn)行對比。

對3種方法的預(yù)測結(jié)果進(jìn)行了定性和定量評價(jià)，前者包括兩種不同波段組合的視覺效果展示，后者包括4種評價(jià)指標(biāo)，即均方根誤差(RMSE)，相關(guān)系數(shù)(correlation coefficient,CC)，通用圖像質(zhì)量指數(shù)(universal image quality index,UIQI)[34]和光譜角距離(spectral angle mapper,SAM)。評價(jià)指標(biāo)均基于云區(qū)域的所有像元計(jì)算。值得注意的是，前3項(xiàng)指標(biāo)均基于單個(gè)波段進(jìn)行計(jì)算，而SAM先基于單個(gè)像元運(yùn)算然后取所有像元結(jié)果的均值。4種評價(jià)指標(biāo)的計(jì)算公式如下。

(1)RMSE。RMSE能夠衡量預(yù)測影像與參考影像之間的差異性，其值越小表明預(yù)測結(jié)果越理想，理想值為0。RMSE定義如下

(2)

式中，Pb(x,y)和Rb(x,y)分別表示大小為m×n像元的影像P和R在位置(x,y)處波段b的像元值。

(2)CC。CC反映預(yù)測影像與參考影像之間的相關(guān)性，其值越接近于1表明預(yù)測影像與參考影像越接近。CC定義如下

(3)

(3)UIQI。UIQI用于評價(jià)預(yù)測影像與參考影像之間空間結(jié)構(gòu)的相似程度，其值越接近于1表明預(yù)測影像的重建質(zhì)量越好。UIQI定義為

(4)

式中，σPbRb為影像P和R在波段b的協(xié)方差；σPb和σRb分別為影像P和R在b波段的標(biāo)準(zhǔn)差。

(4)SAM。SAM用于衡量預(yù)測影像與參考影像在光譜維度上的接近程度

SAM=arccos

(5)

2.3 試驗(yàn)結(jié)果及分析

圖3中以假彩色(同一區(qū)域用同種波段組合)顯示出3種云層去除方法的結(jié)果。由圖3可以看出，3種方法在城市區(qū)(區(qū)域1和區(qū)域2)都以較好的視覺連續(xù)性重建了云下區(qū)域的信息。對于地物發(fā)生更為復(fù)雜變化的農(nóng)業(yè)區(qū)(區(qū)域3)，SSGLHM_M的結(jié)果中非云區(qū)域與云區(qū)域之間存在較為明顯的接縫，而MNSPI_H的結(jié)果中存在著明顯的色調(diào)異常。例如，在圖3中，MNSPI_H將區(qū)域3中部紫色的地物錯(cuò)誤地預(yù)測為粉色和藍(lán)色。相反地，本文提出的SSRF_M在區(qū)域3中同樣表現(xiàn)出令人滿意的性能，其重建結(jié)果在視覺上具備較好的連續(xù)性，且無明顯的色調(diào)異常。為更清楚地對比3種方法的重建結(jié)果，圖3對每個(gè)區(qū)域放大展示了一個(gè)子區(qū)域?？梢钥闯?，MNSPI_H方法的重建結(jié)果均存在著明顯的噪聲，且在區(qū)域1和區(qū)域3中存在著明顯的色調(diào)異常。相較之下，兩種利用多光譜影像作為輔助影像的方法(即SSRF_M和SSGLHM_M)能得到更接近于參考影像的重建結(jié)果。這表示與獲取時(shí)間距目標(biāo)影像更遠(yuǎn)的同源高光譜影像相比，在時(shí)間上更接近的多光譜影像能提供更準(zhǔn)確的輔助信息。進(jìn)一步地，與SSGLHM_M相比，SSRF_M的重建結(jié)果在色調(diào)上與參考影像更為接近。

注：3個(gè)區(qū)域分別以波段150、10、3，163、111、65和161、104、58為RGB組合；白色線框內(nèi)為模擬的云區(qū)域。

為定量評估3種方法，將4種評價(jià)指標(biāo)對應(yīng)的值繪制成條形圖，如圖4所示，其中RMSE、CC和UIQI為其在所有波段的平均值。由圖4可知，在3個(gè)試驗(yàn)區(qū)域中，SSRF_M重建結(jié)果的RMSE最小，CC和UIQI最大，而MNSPI_H方法的RMSE最大，CC與UIQI最小。對于SAM，SSRF_M在區(qū)域1的結(jié)果與GLHM_M幾乎相同，而在區(qū)域2和區(qū)域3中，前者的結(jié)果具有較為明顯的優(yōu)勢。此外，圖5展示了3個(gè)區(qū)域數(shù)據(jù)中每個(gè)波段的RMSE值。由圖5可知，在3種方法的預(yù)測結(jié)果中，SSRF_M在各個(gè)波段下的RMSE幾乎始終處于最下方(即RMSE值最小)。由此可以說明，本文提出的SSRF_M能以更高的精度恢復(fù)云下缺失信息。

注：RMSE、CC和UIQI為各波段的平均值，SAM先基于單個(gè)像元計(jì)算然后取所有像元結(jié)果的均值。

圖5 區(qū)域1—區(qū)域3模擬云層去除結(jié)果各波段RMSE折線圖

此外，為進(jìn)一步檢驗(yàn)各方法預(yù)測光譜的準(zhǔn)確性，在每個(gè)區(qū)域中隨機(jī)選取了兩個(gè)像元，并展示其在參考影像及不同方法預(yù)測結(jié)果中的光譜曲線，如圖6所示?？梢园l(fā)現(xiàn)，SSRF_M預(yù)測的光譜曲線普遍與參考影像更為接近，表明其能夠更大程度地恢復(fù)缺失像元的光譜信息。

圖6 區(qū)域1—區(qū)域3中3種方法預(yù)測的像元光譜曲線

3 結(jié) 論

本文提出一種利用時(shí)間分辨率較高的多光譜影像(Landsat 8 OLI影像)對時(shí)間分辨率較低的高光譜影像(GF-5和EO-1高光譜影像)進(jìn)行云去除的方法(SSRF_M)，一定程度上解決了同源高光譜輔助影像因時(shí)間間隔過大而可用價(jià)值偏低的問題。SSRF_M在描述因變量與自變量之間復(fù)雜非線性關(guān)系時(shí)的優(yōu)勢使其能夠同時(shí)利用輔助多光譜影像的多個(gè)波段數(shù)據(jù)對高光譜各個(gè)云覆蓋波段進(jìn)行重建，無須輔助波段和云覆蓋波段區(qū)間一一對應(yīng)。試驗(yàn)表明，與利用同源高光譜數(shù)據(jù)作為輔助影像的方法(MNSPI_H)及利用Landsat 8 OLI數(shù)據(jù)作為輔助影像的線性擬合方法(SSGLHM_M)相比，SSRF_M能得到更準(zhǔn)確的結(jié)果。值得注意的是，SSRF方法對空間信息的挖掘能力仍有待進(jìn)一步提高，在更大的空間范圍內(nèi)篩選一定數(shù)量的相似像元是一種可能的有效方法。此外，也可考慮發(fā)展其他有效的學(xué)習(xí)模型，在預(yù)測時(shí)充分利用空間鄰域信息。