基于RNMU的多源星載SAR影像融合與土地覆蓋分類

李 俐 陳琦琦 張 超 尤淑撐 魏 海 付 雪

(1.中國農業大學土地科學與技術學院，北京 100083；2.農業農村部農業災害遙感重點實驗室，北京 100083；3.中國國土勘測規劃院，北京 100035；4.中國農業大學信息與電氣工程學院，北京 100083)

0 引言

及時準確地獲取土地覆蓋信息，能夠為土地資源的合理利用與監管提供有力的數據支持[1]。遙感技術由于具有覆蓋范圍廣、獲取頻率高的優勢，在大面積土地覆蓋信息提取方面得到了廣泛應用[2-3]。光學遙感影像憑借其與人眼視覺系統接近、直觀易解譯的特點首先被廣泛應用于土地覆蓋分類[4-6]。然而，光學數據的獲取受到光照條件的限制，在云、雨、雪、霧等天氣下，存在數據獲取困難、品質不佳等問題。以主動方式工作的合成孔徑雷達(Synthetic aperture radar，SAR)擺脫了光照的影響，能夠為土地覆蓋分類提供全天候、全天時微波遙感信息[7-9]。因此，SAR數據在土地覆蓋分類中的應用越來越受到重視。自2014年以來，相繼發射的Sentinel-1A/B和高分三號(GF-3)等SAR衛星，提供了越來越多的高分辨率、多極化、多波段的SAR數據，大大提高了SAR數據在土地覆蓋監測應用中的數據獲取能力。

然而，單星單極化SAR影像的分辨率與影像可解譯能力限制了SAR數據在土地覆蓋分類中的識別效果。多源遙感數據融合技術通過融合不同時相、不同來源的遙感數據，將目標的更多信息從不同角度展現出來，可有效提高土地覆蓋分類精度[10-12]。已有研究將光學影像與SAR影像融合用于土地覆蓋分類[13-15]。SAR與光學影像成像機理不同，簡單地選取不同波段或不同成分進行融合，可能造成信息損失和噪聲引入等問題，進而影響土地覆蓋分類精度[16]。多源SAR影像融合可以較好地規避以上問題，為及時獲取多波段、多極化、多角度的SAR數據提供了依據。目前，對于不同源SAR影像的融合，研究者進行了大量研究。針對多波段SAR影像，德國宇航局(Deutsches zentrum für luft- und raumfahrt，DLR)采用模糊專家系統法進行決策級融合[17]，該方法實時性好，并且具有一定的容錯能力，但預處理代價較高，影像中的原始信息損失較多。為了在盡可能多保留原始信息的同時降低信息冗余，文獻[18]將Gram-Schmidt(G-S)算法用于多時相多波段SAR圖像，有效提高了檢測性能。針對多極化SAR影像，文獻[19]采用非下采樣Contourlet變換(Nonsubsampled contourlet transform，NSCT)和脈沖耦合神經網絡(Pulse coupled neural network，PCNN)相結合的方法進行了像素級融合，保留了NSCT的多尺度性、多方向性、平移不變性和PCNN全局耦合性，但融合模型中參數設置復雜，不當的參數設置易導致影像融合效果不佳。針對多時相SAR影像，文獻[20]利用非負矩陣分解算法(Non-negative matrix factorization，NMF)進行特征級融合，有效保證了矩陣分解結果的可解釋性，實現簡單、且占用存儲空間較小，在SAR影像變化檢測應用中效果較好，但非負矩陣分解后，出現數據的稀疏性降低，因此可能導致融合后影像信噪比下降。

遞歸非負矩陣下近似(Recursive nonnegative matrix underapproximation，RNMU)算法在NMF算法的基礎上進行了改進，有效避免了數據稀疏性的降低，成功應用于高光譜影像的特征提取、分類與光譜混合分析[21]。針對不同源多極化多時相微波影像，本文嘗試將RNMU算法用于SAR影像融合，并將融合數據用于土地覆蓋分類，以有效融合多源信息，并降低數據冗余。

1 研究區概況及數據獲取

1.1 研究區概況

圖1 研究區及采樣點分布

選擇吉林省大安市為研究區。大安市位于123°8′45″～124°21′56″E、44°57′～45°45′51″N，屬于松嫩平原(圖1)，總面積約4 879 km2，海拔為120～160 m，地勢較為平坦。大安市屬中溫帶季風氣候，年平均降雨量413.7 mm[22]，降雨集中在6—9月，期間高品質、高分辨率的光學遙感數據獲取困難。根據大安市政府網提供的2019年土地覆蓋統計數據：現有耕地140 293 hm2；未利用地近217 607 hm2；大安市境內水資源豐富，水域面積占全省的1/7；濕地富集，種類繁多。大安市屬于農、牧生態系統的交界過渡區域，由于土地的不合理開發利用，近年來出現土地退化、濕地減少等生態問題[23]，因此及時了解大安市土地覆蓋現狀，可為土地利用規劃和生態保護提供基礎數據。

1.2 SAR數據源

選取GF-3 SAR數據和Sentinel-1 SAR數據進行多源SAR數據融合與土地覆蓋分類。GF-3衛星于2016年8月發射，是我國首顆最高分辨率可達1 m的C頻段多極化SAR衛星，GF-3 SAR數據有精細條帶2(FSⅡ)等12種成像模式，數據產品包括Level-0至Level-3級標準產品及Level-4級行業應用產品。Sentinel-1是由Sentinel-1A與Sentinel-1B兩顆衛星組成的星座，Sentinel-1(以下簡稱S-1)可提供重復觀測的 C 波段 SAR 數據[24]。S-1 SAR數據成像模式主要有干涉寬幅(IW)模式等4種，數據產品包括SAR Level-0 RAW數據、Level-1單視復數(Single look complex，SLC)、Level-1地距(Ground range detected，GRD)和Level-2海洋(Ocean，OCN)等產品。本文所用SAR數據具體參數見表1。

表1 SAR數據具體參數

1.3 樣本數據獲取

用于分類訓練與精度驗證的樣本數據，由2018年8月實地采集，并結合Google Earth影像人工辨讀選取得到。考慮到研究區為含有豐富濕地資源的糧食主產區，本文選取水體、建設用地、耕地、濕地、其他用地(未利用地、鹽堿地和草地等)共5類作為待分土地覆蓋類型。樣本(圖1)選取遵循平均分布于整個研究區域的原則，對于易分地類選取較少樣本，易混淆地類選用較多樣本，選取訓練樣本數為：水體27個、建設用地41個、耕地72個、濕地69個、其他用地70個。驗證樣本不同于訓練樣本，采用水體14個、建設用地21個、耕地46個、濕地35個、其他用地38個作為驗證樣本。

2 研究方法

2.1 技術流程

針對輸入的GF-3雙極化SAR數據和S-1雙極化SAR數據，在預處理的基礎上，進行影像配準、影像融合及土地覆蓋分類，處理流程如圖2所示。

圖2 技術流程圖

經預處理得到2顆衛星4種極化方式的后向散射系數圖像后，采用互相關[25](Cross-Correlation, CC)算法對圖像進行精確配準處理。然后，將RNMU方法引入多源SAR圖像融合并用于研究區土地覆蓋分類。作為對比，另一種同樣具有較小冗余度的融合算法G-S被用于融合和土地覆蓋分類。G-S融合方法基于多維正交線性變換, 通過正交化消除影像的冗余信息，并在第一正交分量中加入另一景影像信息，從而實現多源影像融合。

2.2 SAR數據預處理

SAR數據預處理主要包括輻射定標、斑點噪聲濾波、地形校正等。

相干斑噪聲是SAR遙感圖像獨有的乘性噪聲，會對影像解譯造成較大的干擾，斑點噪聲濾波是為了降低相干斑噪聲。由于GF-3和S-1衛星搭載不同的SAR傳感器且成像模式不同，導致不同源SAR影像相干斑噪聲不盡相同，因此本文在分別嘗試了若干種濾波器的基礎上，選擇窗口大小為5×5的Gamma濾波降低GF-3影像噪聲，在盡可能地保持影像空間細節的情況下提高信噪比[26]；選擇Refined Lee濾波降低S-1影像斑點噪聲，并將濾波窗口設為7×7[27]。

(1)

式中QV——該景影像量化前的最大值，可通過元數據獲取[28]

KdB——定標參數，通過元數據獲取

(2)

式中A——定標參數，通過元數據文件獲取[29]

2.3 RNMU多源SAR影像融合算法

2.3.1RNMU基本原理

RNMU算法核心思想是：對于任意l×k(l、k∈R+)非負實數的輸入矩陣W=[W1W2…Wk]，將W分解為非負的l×r的權重矩陣u和k×r的特征矩陣v，使得uvT≈W，求解u、v的優化過程即為W與uvT間距離的元素平方和最小化迭代求解的過程

uvT≤W

針對上述目標函數求解過程中可能出現的局部最優和病態優化問題，RNMU算法限制了權重矩陣和特征矩陣秩的大小，只保留r=1的矩陣，并增加約束條件uvT≤W。對于影像輸入來說，低秩性有效地去除了影像中的噪聲。每步僅計算秩為1的矩陣，最優解唯一，有效地避免了NMF算法應用于遙感影像時近似秩改變帶來的計算重復、計算效率低下以及融合后影像的品質易受初始化矩陣影響等問題。待融合的雷達影像中，每景影像對應特定的時相和極化方式，增加的約束條件保障了所提取的信息為原始數據的根本特征，使結果的稀疏性和局部特征的可分離性更優[30]。

(3)

式中m、n——輸入矩陣的行數和列數

式(3)的求解問題轉換為u、vT和Λ迭代優化的過程：

(1)初始化Λ為零矩陣O，初始化u、v為秩為1的向量，且滿足uvT≤W。

(2)更新u、vT

(3)用次梯度優化算法迭代更新拉格朗日乘子Λ

Λ=max(O,Λ+α(W-uvT)) (α→0)

(4)

(4)迭代步驟(2)、(3)，直到系統穩定，所得結果即為理想u、vT。

2.3.2基于RNMU的SAR影像融合

將GF-3單個極化影像和已配準的S-1單個極化影像分別作為輸入影像a、b，按列展開方式排列為兩個列向量Wa、Wb，構造出輸入矩陣W

[a11a21…am1…a1n…amn]T=Wa

(5)

b?[b11b21…bm1…b1n…bmn]T=Wb

(6)

W=[WaWb]

(7)

v=[v11v21…vm1…v1n…vmn]?

(8)

2.3.3融合結果評價指標

針對多源SAR影像融合，本文選用標準差(Standard deviation，SD)、平均絕對誤差(Mean absolute error，MAE)、信息熵(Entropy，H)[32]、信噪比(Signal-to-noise ratio，SNR)和峰值信噪比(Peak signal-to-noise ratio，PSNR)[33]等評價指標，分析評價RNMU融合影像的品質，其中，SD反映了融合前后影像的灰度分布情況；MAE反映融合后影像與原始影像之間的差異；H反映融合后影像信息豐富程度，SNR和PSNR反映影像噪聲情況。

3 實驗與結果分析

3.1 融合結果評價

為定量評價RNMU算法在多源SAR影像融合的效果，針對RNMU融合(R:HHVV, G:HVVH, B:HHVV-HVVH)影像與6月S-1(R:VV, G:VH, B:VV-VH)影像、6月GF-3(R:HH, G:HV, B:HH-HV)影像和G-S融合(R:HHVV, G:HVVH, B:HHVV-HVVH)影像，分別計算SD、MAE、H、SNR和PSNR等指標以比較融合前后和不同融合算法的性能差異。RNMU和G-S融合影像中，R：HHVV波段數據為GF-3 HH極化SAR數據和S-1 VV極化SAR數據分別經RNMU和G-S算法融合得到的雙極化SAR產品，G：HVVH波段數據為GF-3 HV極化SAR數據和S-1VH極化SAR數據分別經RNMU和G-S算法融合得到的雙極化SAR產品，融合后雙極化產品的差值HHVV-HVVH作為B波段。如表2所示，G-S融合和RNMU融合后影像與S-1和GF-3影像相比，SD較大表明融合后影像的灰度級分布較分散，包含的信息量較多；H較大表明兩種融合方法都增強了影像的細節表現能力。RNMU融合和G-S融合結果相比，RNMU融合后影像MAE較小，表明了RNMU融合算法具有更好的原始影像信息保留能力；而RNMU融合后的SNR和PSNR大于G-S融合后影像，體現了RNMU融合方法良好的影像噪聲抑制能力。GF-3的HH影像和S-1的VV影像經RNMU融合后SD較大，影像具有較大的對比度，影像灰度級分布較為分散，有利于地物分類。而對于兩個交叉極化的影像的融合，G-S算法得到的影像SD大于RNMU算法的結果。由于一般交叉極化信息弱于同極化，分類應用中同極化數據應用更為廣泛，因此，RNMU融合影像更適于土地覆蓋分類。

表2 不同融合算法定量評價

3.2 土地覆蓋分類結果與分析

3.2.1典型地物類型散射特性分析

地物的微波后向散射特性是利用SAR影像進行土地覆蓋分類的基礎[34]。圖3給出研究區6月GF-3雙極化(HH、HV)影像和S-1雙極化(VH、VV)影像以及10月S-1雙極化(VH、VV)影像中典型地物的后向散射系數分布圖。可以看出，無論是在兩時相S-1的VH和VV影像還是在單時相的GF-3 HH和HV影像中，水體均具有最低的后向散射系數，建設用地則具有較高的后向散射系數，而耕地、其他用地和濕地的后向散射系數較為接近，特別是在HH極化的GF-3數據中，水體與建設用地具有和其他的土地覆蓋類型良好的可分性。而耕地、其他用地和濕地的后向散射系數較為接近，耕地與濕地在S-1的VH極化下典型地物后向散射系數分離度略優于VV極化，但在GF-3的HH極化后向散射系數分離度優于HV極化，因此，單景影像中后向散射系數較為接近的耕地、濕地和其他用地需要利用SAR影像不同極化、不同時相的變化特性進行進一步區分。

圖3 不同時相不同極化影像典型地物后向散射系數分布圖

圖4 典型地物后向散射系數的均值和標準差的比較

為了更直觀看出不同地物的散射特性差異，圖4給出了融合前GF-3 SAR影像和S-1 SAR影像與RNMU融合后SAR圖像的均值和標準差。其中基于RNMU的融合，分別給出了研究區10月S-1 VV極化影像和6月GF-3 HH極化影像經RNMU算法融合得到的VVHH影像，10月S-1 VH極化影像和6月GF-3 HV極化影像經RNMU算法融合得到的VHHV影像(VVHH影像和VHHV影像疊加，簡寫為RNMU1)。6月S-1 VV極化影像、10月S-1 VV極化影像和6月GF-3 HH極化影像經RNMU算法融合得到VVVVHH影像，6月S-1 VH極化影像、10月S-1 VH極化影像和6月GF-3 HV極化影像經RNMU算法融合得到VHVHHV影像(VVVVHH影像和VHVHHV影像疊加，簡寫為RNMU2)。由圖4可以看出，6月的S-1影像中濕地和耕地上的植被長勢較好，因此比其他用地具有更高的后向散射系數，10月研究區的作物成熟或收獲，耕地的后向散射系數降低明顯，不同地物隨時間的變化特性差異相對明顯。而經RNMU融合后的影像，無論是RNMU1影像還是RNMU2影像各地類均值之間差異增大，因此增加了地類之間的可分性；而每種地類的標準差較S-1影像降低，進而提升了待融合影像(S-1)地類特征的一致性。因此利用不同時相、不同極化的多源SAR影像進行融合，能融合更多互補信息，得到更好的分類結果。

3.2.2分類結果分析

為分析融合影像在土地覆蓋分類應用中的性能改進，分別對研究區 RNMU1影像，RNMU2影像，6月、10月S-1影像和6月GF-3影像經G-S融合后影像(簡寫為G-S)，應用最大似然分類器(Maximum likelihood classifier，MLC)進行土地覆蓋分類，作為比較，同時給出6月單景S-1影像，6月單景GF-3影像，6月、10月S-1影像和6月GF-3經簡單波段疊加后影像的土地覆蓋分類結果。由圖5可以看出，利用G-S和RNMU2影像識別的水體基本相同，建設用地、濕地和其他用地存在不同程度差別。選取土地覆蓋類型多樣的區域(圖5矩形框內區域)作為驗證區域，利用3 m分辨率的Google高清圖像與S-1影像、G-S影像和RNMU2影像進行分類結果實驗驗證。如圖6所示，與S-1影像分類結果相比，G-S影像和RNMU2影像分出的地類邊界較清晰。RNMU2影像分出的建設用地(紅色橢圓框圈出)更接近真實值，G-S影像分類結果次之，S-1影像分出的建設用地較為破碎。

圖5 不同數據土地覆蓋分類圖

圖6 驗證區域的土地覆蓋分類比較

利用驗證樣本對分類結果進行精度評價，各數據的分類總體精度(OA)、制圖精度(PA)、用戶精度(UA)和Kappa系數(Kc)如表3所示。由表3可知，GF-3數據的土地覆蓋分類整體精度優于同時相S-1數據，不同時相不同源(S-1和GF-3)SAR數據經RNMU融合后數據分類OA分別由單時相SAR數據74.96%(S-1結果)和83.52%(GF-3結果)提高至86.08%(RMNU1影像結果)；Kc分別由0.55和0.71提高至0.74；耕地、建設用地和濕地的分類精度都有明顯提高。由于研究區濕地種類繁多，同時間不同類別的濕地變化情況不一。雖然加入其他時相數據可以提高其識別精度，但其分類精度與其他地物相比仍較低。使用相同數據源的情況下，G-S融合后分類結果要低于所有SAR數據直接疊加或者RNMU融合后的結果，可能是因為G-S融合圖像為了降低數據冗余，犧牲了一定的紋理信息，導致紋理特征明顯的耕地UA較低，而耕地在研究區面積占比較大，導致OA和Kc整體下降。基于RNMU融合后影像(RNMU2)地類識別效果，整體優于基于G-S算法融合結果和簡單波段疊加后影像，說明RNMU融合算法充分利用多時相、多極化SAR數據的互補信息，有效提高土地覆蓋分類精度，且融合后影像數據量減少了1/3，實現了數據壓縮，適用于較大區域的SAR數據。但對于易于分辨的水體地物，單時相S-1 SAR數據的水體分類PA高于RNMU2的結果，而融合后并沒有達到較好的效果，說明對于后向散射特征足夠明顯的地物，融合后數據可能出現過擬合現象，反而可能引入噪聲，降低分類精度。

表3 不同數據融合方式分類精度評價

注：表頭第2行表示數據的月份、來源以及分類所用的極化影像。

4 結論

(1)為充分利用多時相、多極化SAR互補數據，將RNMU算法引入多源SAR融合。RNMU融合后，影像SD、H均優于原始S-1、GF-3雙極化SAR影像；與G-S融合算法相比，MAE、H、SNR和PSNR均體現出優良的特性。說明RNMU融合算法有效降低了SAR數據噪聲，提高了融合后影像的信噪比，表明RNMU融合算法在降低數據冗余的同時有效融合了多源數據之間更多的互補信息。

(2)RNMU融合后，SAR影像綜合利用多時相、多極化SAR信息提高了土地覆蓋分類精度。融合兩時相SAR數據進行土地覆蓋分類的精度明顯高于融合前SAR數據，表明RNMU融合后多源SAR數據更適于土地覆蓋分類應用。

(3)對于后向散射特征明顯的水體地物，RNMU融合數據分類的PA較單時相S-1略有下降，這說明對于明顯特征地類可能存在過擬合現象。在進一步融合研究中，將不同區域、不同地類設定不同的權重，以期在增加互補信息的同時，降低噪聲干擾。