多圖正則多核非負矩陣分解高光譜圖像解混

劉 敬，李康欣，張 悠，劉 逸

（1. 西安郵電大學 電子工程學院，陜西 西安 710121；2. 西安電子科技大學 電子工程學院，陜西 西安 710071）

1 引言

高光譜圖像［1］有高的光譜分辨率，包含豐富的圖像及光譜信息，而由于高光譜傳感器的低空間分辨率和不同純物質波譜的混合，導致混合像元［2］的產生，極大地影響了高光譜遙感圖像的應用。為改善高光譜圖像分解精度，高光譜解混［3-5］已成為熱點，可用線性或非線性方式將混合像元分解，同時提取端元與豐度。端元是混合像元分解出的純物質光譜，而豐度［6］則是每個像元中端元的貢獻。

早期高光譜解混算法主要采用線性混合模型（Linear Mixture Model，LMM），如基于幾何的頂點成分分析法（Vertex Component Analysis，VCA）［7］、最小封閉體積的單純形法（Minimum Volume Enclosing Simplex，MVES）［8］，基于統計的貝葉斯方法、獨立成分分析（Independent Component Analysis，ICA）［9］和非負矩陣分解（Nonnegative Matrix Factorization，NMF）［10］算法等。在NMF 算法中，先找到一組非負基，然后將原始數據映射到這組基上，且數據在每個基上的表達非負。NMF 十分適合應用于高光譜解混，得到端元和豐度矩陣。為保持數據空間固有的流形結構，Yang 等人在NMF 中加入了圖正則算法，稱為圖正則非負矩陣分解［11］，該算法利用一種圖相似性描述樣本之間的關系，充分考慮了局部不變性。由于定義樣本之間的關系及其權重矩陣的方法很多，單圖的選擇至關重要。為解決圖選擇的問題，Shu 等人提出多圖正則非負矩陣分解（Multi-graph Regularized Nonnegative Matrix Factorization，MGNMF）［12］，采用一組已知的多個圖拉普拉斯矩陣，通過學習得到的加權參數組合，去逼近原始數據。多圖可更準確地刻畫樣本的相似性，進而更好地表達原始數據的結構。雖然多圖解決了圖正則非負矩陣分解（Graph Regularized Nonnegative Matrix Factorization，GNMF）算法中圖選擇困難的問題，但其屬于線性解混算法［13-14］，難以適應真實場景中復雜的非線性光譜混合結構，所以，這一問題也促進了非線性解混算法研究。

NMF 是一種典型的線性解混算法，不適合提取數據的非線性混合結構［15］，而核方法可以解決這個問題。核方法［16］是將非線性關系轉變為線性關系的一種過程方法，將低維非線性混合結構的數據映射到高維核空間，在核空間中數據混合結構呈現線性，因此在核空間可實現高光譜圖像的非線性解混。Yan 等人提出的包含純像元的核非負矩陣分解（Pure-pixels Kernel Nonnegative Matrix Factorization，pKNMF）與不含純像元的核非負矩陣分解（Non-pure-pixel Kernel Nonnegtive Matrix Factorization，npKNMF），是將KNMF［17］算法直接應用于高光譜數據，得到比NMF 好的分類效果。但KNMF 的性能很大程度上取決于核函數的選擇。多核非負矩陣分解（Multi-kernel Nonnegtive Matrix Factorization，MKNMF）算法采用多個核函數的組合，并為每個核函數設置適當的權重參數。相比KNMF 算法，多核NMF 可自適應地選擇核函數并加權，有更好的映射能力。Yao 等人將MKNMF 引入圖正則NMF 算法中，在高光譜數據集中得到了更好的驗證［18］。

在許多復雜自然場景中存在大量的非線性混合現象，如沙地和礦物混合區的密集混合現象，以及植被和建筑物覆蓋區的多級混合現象。基于線性混合模型的線性解混算法不適合于非線性混合情況，所以，本文提出了一種非線性解混算法——多圖正則多核非負矩陣分解（Multigraph Regularized Multi-kernel Nonnegative Matrix Factorization，MGMKNMF），先使用多核函數構造適合于高光譜數據的核空間，然后在多核空間的基礎上為目標函數添加多圖正則項。本文提出的算法有以下兩個優點：（1）與KNMF 算法相比，MGMKNMF 算法采用核函數權重將多個不同參數的高斯核函數聯合起來，并在學習過程中不斷更新核函數權重，避免了單核的唯一性，使構造的核空間更合適，也解決了多個高斯核函數權重選擇困難的問題。（2）與GNMF 和MGNMF 算 法 相 比，MGMKNMF 算 法 是 非 線 性方法，在多核空間構造多圖，圖權重將多個圖拉普拉斯矩陣線性組合，并與豐度矩陣最終構成多核空間的多圖正則項，且在學習過程中不斷更新圖權重。相比原空間的單圖和多圖，多核空間中的多圖可更精確地刻畫原始數據的非線性流形結構，更適合對真實場景中復雜的非線性光譜混合結構進行非線性解混。本文采用兩個真實高光譜數據集Jasper Ridge 和Cuprite，并采用廣義雙線性模型（GBM）［19］和Hapke［20］非線性模型分別生成兩個模擬數據集，將所提MGMKNMF 算法與GNMF、npKNMF、核稀疏非負矩陣分解（Kernel Sparse Nonnegative Matrix Factorization，KSNMF）、基于核的字典剪枝非線性光譜解混（Kernel-based Nonlinear Spectral Unmixing with Dictionary Pruning，KDP）、多圖正則核非負矩陣（Multi-graph Regularized Kernel Nonnegative Matrix Factorization， MGKNMF）算法比較，實驗結果表明，MGMKNMF的光譜角距離（Spectral Angel Distance，SAD）和 均 方 根 誤 差（Root Mean Square Error，RMSE）相比其他算法均有較為顯著的下降。

2 相關工作

2.1 非負矩陣分解NMF

NMF 可用于高光譜遙感影像的無監督解混。給定一個數據矩陣X=[x1，…，xN]∈RM×N，其中X的每一列都是樣本向量。NMF 通過矩陣分解將原始高秩矩陣分解為兩個低秩矩陣A∈RM×K和S∈RK×N的相乘，并加上非負的約束。

非負矩陣分解的目標函數：

其中：X是原始數據矩陣，在高光譜數據中，A和S分別代表端元矩陣與豐度矩陣。

2.2 圖正則非負矩陣分解GNMF

NMF 將非負矩陣X分解為基矩陣A和編碼矩陣S的乘積(X≈AS)，用于高光譜圖像的無監督解混時，即將高光譜數據集X分解為端元矩陣A與豐度矩陣S的乘積。GNMF 將圖正則項添加到NMF 的目標函數中，改善了未考慮數據流形結構所帶來的問題。對有N個樣本的數據集X構造一個K近鄰圖，圖中頂點為X中各像素點。Nn為X中樣本xn的K近鄰集，將每個頂點xn與屬于它的Nn連接，并定義圖權重矩陣W∈RN×N，原空間中頂點xn和xm間的權重Wnm越大，在子空間中的sn和sm距離也越近。通過權重矩陣W構造圖正則項：

其中：sn和sm為像元xn和xm在端元基向量上的編碼向量，即豐度；Wnm為xn和xm間的權重；S為豐度矩陣；D為對角矩陣，D-W是圖拉普拉斯矩陣；“tr（·）”表示矩陣的對角線元素之和。常用定義權重矩陣W的方法有：0-1 加權、熱核加權和點積加權等。

將式（3）與（2）結合，得到GNMF 的目標函數：

其中：α為權重參數，GNMF 求解約束最小化問題minOGNMF(A，S，W) s.t.A≥0，S≥0。

2.3 多圖正則非負矩陣分解MGNMF

根據不同定義權重矩陣的模型，可計算相應的圖權重矩陣和圖拉普拉斯矩陣。MGNMF 采用不同數量的最近鄰構建圖并進行加權，經過算法自動選擇，得到最優多圖正則項。若已知一組M種模型的圖權重矩陣｛W1，W2，…，WM｝和相應的圖拉普拉斯矩陣｛L1，L2，…，LM｝，將這M個圖權重矩陣線性組合，則相應的圖拉普拉斯矩陣也進行相同的線性組合：

其中，γm是第M個圖權重矩陣和圖拉普拉斯矩陣的權重。MGNMF 為一組預先計算得到的候選圖確定最佳的圖權重，而不是先選最佳圖矩陣模型并估計參數。MGNMF 的多圖正則項為：

多圖正則項比單圖正則項更精準，GNMF 的唯一權重不可靠。MGNMF 無需選擇唯一的圖權重矩陣模型，且通過學習所得圖權重向量γ=[γ1，…，γL]T對M個圖拉普拉斯矩陣進行最優線性組合。

2.4 多核非負矩陣分解MKNMF

NMF 是線性方法，不能很好地處理數據中的非線性結構，而KNMF 通過核方法，將原始數據 映 射 到 高 維 核 空 間 ：X→Φ(X)=[Φ(x1)，…，Φ(xn)]，可以解決數據的非線性問題。核矩陣K=Φ(X)TΦ(X)，K∈RN×N。則在高維核空間，NMF 可以表示為Φ(X)≈AS。其中A為端元矩陣，S為豐度矩陣。以核空間中樣本Φ(X) 作為基向量，得到端元矩陣A=Φ(X)F，A=[a1，a2，…，aP]，P為端元個數，F矩陣包含核空間中所有樣本對構造各端元的貢獻，F∈RN×P且為矩陣F中 第n行p列個元素，F矩陣中的第p列為核空間中所有樣本對構造第p個端元的貢獻。當核函數確定后，KNMF 有唯一的核空間，這種選擇核空間的方式并不準確。

MKNMF 算法將L個不同核函數對應的核空間聯合起來，以構造一個更合適的希爾伯特空間，這L個不同核空間對應的核矩陣為核函數權重向量τ=[τ1，…，τL]T將這L個不同的核空間線性組合，組合后的多核核空間的核矩陣為：

通過將學習好的參數τ帶入上式，避免了不同核函數權重分配的問題；而多個核函數可構造出更適合原數據的核空間，比單核更可靠。將式（4）代入KNMF 目標函數中，得MKNMF 的目標函數：

3 MGMKNMF

本文提出的MGMKNMF 在多核空間構造多圖，為更新多圖，給定參數τ，用歐幾里德距離的平方重新定義多核空間中xn的K近鄰集Nτ n：

多核空間中的多圖正則化項為：

最終將式（8）與式（10）結合，得到MGMKNMF 的目標函數：

多核空間可更好地挖掘數據間的非線性關系，在多核空間嵌入多圖能更好地表達數據的非線性流形結構。式（11）中的‖τ‖2可防止參數τ過度偏向到一個核函數中；‖γ‖2項可避免γ偏向到一個權重構造函數中；α，β，μ均是權衡上式所用的權衡參數，其值均為非負。因α是約束多圖正則項的參數，與約束核權重和圖權重的β，μ相比，α應大于β，μ。根據參考文獻［21］和多次實驗的結果，本文實驗中：Cuprite 數據α=100，Jasper Ridge 數據α=20，HAPKE 模擬數據α=20，GBM 模擬數據α=20；所有數據的β均為10，μ均為10。

根據參考文獻［21］和［22］，式（11）中的MGMKNMF 目標函數是非凸的，無法得到全局最小值。MGMKNMF 采用分步迭代策略優化目標函數，可得到局部最小值［21-22］，具體如下：

（1）固定τ和γ，更新F，S。式（11）可寫為：

假設φlp和ψpn分別是flp和spn的拉格朗日乘子，令φlp=φ，ψpn=ψ，則式（12）的拉格朗日函數為：

式（13）分別對F和S求偏導得：

結合卡羅需-庫恩-塔克（KKT）條件：φlp flp=0，ψpn spn=0 有：

得到F和S的乘法更新法則分別為：

（2）固定F，S和γ，更新τ，式（11）可寫為：

其中，Z=I-FS，gl=tr[Kl ZZT]。

（3）固 定F和S和τ，更 新γ，式（11）可 以寫為：

式（18）和式（19）的約束二次規劃問題可根據文獻［23］中的方法來解決。

MGMKNMF 解混算法總結如下：

MGMKNMF 解混算法輸入：原始高光譜數據X，L 個核矩陣{ K1，…，Kl，…，KL}，M 個圖拉普拉斯矩陣{ L1，…，Lm，…，LM}，最大迭代次數T Step1. 初始化矩陣F0和S0，Step2. 初始化核權重變量τ0 l =1 L，l=1，…，L，初始化圖權重γ0m=1 M，m=1，…，M for t=1 to T do通過式（7）更新圖Gτt 和相應的拉普拉斯矩陣Lτt；通過式（11）和式（12）更新Ft 和St；通過式（13）和式（14）更新核權重τt 和圖權重γt；end輸出：F=Ft-1，S=St-1，τ=τt-1，γ=γt-1

4 實驗結果與分析

將 所 提 MGMKNMF 算 法 與 GNMF、npKNMF、KSNMF、KDP 和MGKNMF 解 混 算法進行對比，使用SAD 和RMSE 作為評估指標，采用HAPKE 和GBM 模擬數 據，以及Cuprite 和Jasper Ridge 真實數據驗證該算法的有效性。所有實驗中，我們選擇0-1 加權圖、熱核加權圖和點積加權圖的圖權重矩陣模型構成最終圖正則項。多核函數選擇高斯核函數，所選取的核參數為1/32，1/16，1/8，1/4，1/2，1，2，4，8，16，32 以構成不同的核函數。各算法的迭代次數T均設置為200。

4.1 評價標準

SAD值反映了解混所得端元光譜與原端元光譜之間相似性，定義為VSAD=arccos(aTb/||a||||b||)，其中，a，b是兩個端元光譜。

RMSE 值反映解混所得豐度與實驗室測量的實際豐度間的差別，定義為VRMSEsi=(1/N|si-s^i|2)1/2，其中，si和s^i分別是實際豐度和解混所得豐度。

4.2 模擬數據實驗與分析

本文采用HAPKE 和GBM 兩種非線性模型生成模擬數據集。HAPKE 模型是一種緊密混合模型；GBM 是一種雙線性混合模型，是LMM 線性結構與端元間的二次散射項的加權組合。本文從美國光譜庫（USGS）隨機選擇6 種地物光譜作為端元，如圖1。并且豐度矢量滿足豐度非負和豐度和為一的約束。最后模擬生成的高光譜每個大小為20×20，每個像素波段為224，并加入不同信噪比（Signal to Noise Ratio，SNR）的零均值高斯噪聲來更好地接近真實數據。

圖1 光譜庫中隨機生成的端元光譜Fig.1 Endmember spectra randomly generated by spectrum library

4.3 模擬數據實驗結果及分析

表1～4 分別顯示了端元數目為6，SNR 不同時，各算法在HAPKE 和GBM 模擬數據上的SAD、RMSE 值。從表中可以看出，隨著SNR 值的增加，SAD 和RMSE 的值均呈下降趨勢，而MGMKNMF 算法與其他6 種算法相比具有較好的準確度。NMF 和GNMF 因為屬于線性解混算法，對非線性數據解混精度都很差，但GNMF算法在NMF 的基礎上增加了圖正則約束，揭示了數據內在固有的流形結構，相比NMF 算法有進一步提升；剩下的5 種算法均利用了核函數的概念，更適合于非線性數據。npKNMF 算法在核空間中應用了NMF，相比NMF 準確度有所提高；而KSNMF 在豐度中加了L1范數使豐度更加稀疏，得到了更好的準確率；KDP 在進行端元選擇時使用了大型光譜庫，相比其余兩者更加準確，MGKNMF 在單核空間中構造多圖正則項。而提出的MGMKNMF 算法不僅用多圖來刻畫數據內在流形結構，更使用多核學習找到合適的核參數與核函數，進一步增加了算法的準確度。

表5～8 分別顯示了SNR 為40 dB，端元數目P不同時，各算法在這2 個模型的模擬數據上的SAD 與RMSE 值。因每次生成模擬數據時是在USGS 庫隨機選取端元，導致表5～8 的結果在P=6 時，與表1～4 中SNR 為40 dB 時的結果不同。可以看出，因為模擬數據非線性程度高，端元數的增加會導致算法解混能力的下降。相比NMF 算 法，MGMKNMF 算 法 在HAPKE 數 據中，端元為6 時，SAD 值和RMSE 值分別減少了約0.17 和0.13；在GBM 數據中，分別減少了約0.12 和0.17。和其余算法相比，MGMKNMF 算法也基本保持著最優結果。

表1 不同SNR 下HAPKE 模型各算法的SAD 值Tab.1 SAD value of each algorithm of HAPKE model under different SNR

表2 不同SNR 下GBM 模型各算法的SAD 值Tab.2 SAD value of each algorithm of GBM model under different SNR

表3 不同SNR 下HAPKE 模型各算法的RMSE 值Tab.3 RMSE value of each algorithm of HAPKE model under different SNR

表4 不同SNR 下GBM 模型各算法的RMSE 值Tab.4 RMSE value of each algorithm of GBM model under different SNR

表5 不同端元數目下HAPKE 模型各算法的SAD 值Tab.5 SAD value of each algorithm of HAPKE model under different number of endmembers

表6 不同端元數目下GBM 模型各算法的SAD 值Tab.6 SAD value of each algorithm of GBM model under different number of endmembers

表7 不同端元數目下HAPKE 模型各算法的RMSE 值Tab.7 RMSE value of each algorithm of HAPKE model under different number of endmembers

表8 不同端元數目下GBM 模型各算法的RMSE 值Tab.8 RMSE value of each algorithm of GBM model under different number of endmembers

4.4 真實數據實驗結果及分析

本文采用真實地物Cuprite 和Jasper Ridge 數據集對MGMKNMF 的有效性進行驗證，兩個數據 集 均 可 以 在https：//rslab. ut. ac. ir/data 中下載。

Cuprite 是高光譜解混研究常用的數據集，包含美國內華州Cuprite 礦區。在除去低信噪比和吸水通道后，留有188 個通道可以使用。每張圖像大小為250×191，共有12 種類別。

表9 總結了在Cuprite 數據集各類算法的SAD 值。可以看出，MGMKNMF 算法的平均SAD 值是最優的。

圖2 為MGMKNMF 算法在Cuprite 數據集上的豐度圖。結合表9，在此算法下，12 種端元對應的豐度圖應是最優的，可清晰地看出各類別的劃分情況。

圖2 MGMKNMF 算法在Cuprite 數據上的豐度圖Fig.2 Abundance graphs of MGMKNMF algorithm on Cuprite

表9 不同算法在Cuprite 數據的SAD 值Tab.9 SAD values of Cuprite data by different algorithms

由于Jasper Ridge 數據集太復雜無法處理，我們僅考慮100×100 像素的子圖像，子圖像的第一個像素從原始圖像中的第（105，269）像素開始，且去除了低噪聲和水吸收波段，保留198 個通道。所以，該高光譜遙感數據中只有樹、水、土壤和道路這4 類端元。

從表10 可以看到，提出的MGMKNMF 算法在Jasper Ridge 數據集上仍然有效，平均SAD 值依舊是最優的。KDP 算法在KSNMF算法的基礎上利用光譜庫挑選端元，而MGKNMF 在單核空間應用多圖，精度僅次于MGMKNMF。

表10 不同算法在Jasper Ridge 數據的SAD 值Tab.10 SAD values of Jasper Ridge data by different algorithms

圖3 為Jasper Ridge 數據集各算法解混出的豐度圖，從左至右分別對應樹、水、土壤和道路這4 個端元。

圖3 各算法在Jasper Ridge 的豐度圖Fig.3 Abundance graphs of each algorithm on Jasper Ridge

從圖3 可看出，NMF 算法與真實地物豐度相差很大，尤其在水這一端元上；GNMF 算法比NMF 算法的豐度圖清晰。而其余的5種算法均在核空間進行，豐度圖明顯優于NMF 和GNMF，其中水的豐度圖更接近真實地物。相比其他算法，在核空間進行的算法盡管其豐度圖的區別甚微，但結合表10 的數據，可 以 看 到，MKNGNMF 算 法 的SAD 值 是最優的。

5 結論

本文提出了一種MGMKNMF 高光譜非線性解混算法。該算法用多個核函數構造出了多核空間，且在學習過程中不斷更新核函數權重，更有利于揭示原始數據的非線性結構；并在多核空間用圖權重向量將多個圖拉普拉斯矩陣線性組合，與豐度矩陣最終構成多核空間的多圖正則項。相比單圖正則項，多核空間的多圖正則項將更加逼近原始數據的非線性流形結構。基于2 個真實數據集和2 個模擬數據集的實驗結果表明，相比GNMF、npKNMF、KSNMF、MGKNMF 和KDP 算法，MGMKNMF 算法確實是最優的，它更合適于復雜的高光譜數據。