基于稀疏和正交約束非負矩陣分解的高光譜解混

陳善學 儲成泉

摘 要：針對基于非負矩陣分解（NMF）的高光譜解混存在的容易陷入局部極小值和受初始值影響較大的問題，提出一種稀疏和正交約束相結合的NMF的線性解混算法SONMF。首先，從傳統的基于NMF的高光譜線性解混方法出發，分析高光譜數據本身的理化特性;然后，結合豐度的稀疏性和端元的獨立性兩個方面，將稀疏非負矩陣分解（SNMF）和正交非負矩陣分解（ONMF）兩種方法結合應用到高光譜解混當中。模擬數據和真實數據實驗表明，相比頂點成分分析法（VCA）、SNMF和ONMF這三種參考解混算法，所提算法提高了線性解混的性能;其中，評價指標光譜角距離（SAD）降低了0.012～0.145。SONMF能夠結合兩種約束條件的優勢，彌補傳統基于NMF線性解混方法對高光譜數據表達的不足，取得較好的效果。

關鍵詞：非負矩陣分解; 高光譜解混; 稀疏; 正交;獨立性

中圖分類號：?TN79

文獻標志碼：A

Hyperspectral unmixing based on sparse and orthogonal constrained non-negative matrix factorization

CHEN Shanxue， CHU Chengquan*

School of Communication and Information Engineering， Chongqing University of Posts and Telecommunications， Chongqing 400065， China

Abstract：?Aiming at the problem that hyperspectral unmixing based on Non-negative Matrix Factorization （NMF） is easy to fall into local minimum and greatly affected by initial value， a linear unmixing algorithm based on Sparse and Orthogonal constrained Non-negative Matrix Factorization （SONMF） was proposed. Firstly， based on the traditional NMF hyperspectral linear unmixing method， the physical and chemical properties of the hyperspectral data was analyzed. Then the sparsity of the abundance and the independence of the endmember were combined together， two methods of Sparse Non-negative Matrix Factorization （SNMF） and Orthogonal Non-negative Matrix Factorization （ONMF） were combined and applied into hyperspectral unmixing. The experiments on simulation data and real data show that， compared with the three reference unmixing algorithms of Vertex Component Analysis （VCA）， SNMF and ONMF， the proposed algorithm has improved the performance of linear unmixing， in which the Spectral Angle Distance （SAD） is reduced by 0.012 to 0.145. SONMF can combine the advantages of the two constraints to make up for the lack the expression of hyperspectral data by traditional NMF based linear unmixing methods， and achieve good results.

Key words：?Non-negative Matrix Factorization （NMF）; hyperspectral unmixing; sparsity; orthogonality; independence

0 引言

近年來，傳感器技術的快速發展大大加快了高光譜遙感圖像處理技術在遙感領域的實際應用，高光譜遙感技術的快速發展使得遙感領域的許多方向發生了翻天覆地的變化。特別是在對生活環境特別重視的現在，高光譜技術在探測地表環境、地質災害等領域發揮著越來越關鍵的作用[1]。但是，由于遙感傳感器的空間分辨率較低，相鄰物質可能混合在一起形成混合像元。高光譜解混[2]就是將混合像元分解為各種地物的端元光譜以及對應的豐度系數。由于在混合像元分解時有關端元和豐度的任何信息都是未知的，而且受到復雜的自然條件和地物環境等因素的影響，混合像元分解問題變得更加困難，所以在解決混合像元分解問題之前，首先要確定光譜混合模型。

線性光譜混合模型（Linear Spectral Mixed Model， LSMM）[3]假定端元之間互相獨立，是一個描述混合像元現象的有效模型。非負矩陣分解（Non-negative Matrix Factorization， NMF）[4]在基于LSMM的高光譜解混中得到廣泛的關注，它采用乘子算法將數據分解為兩個非負矩陣，不需要每個端元都存在純像素，非負的約束自然保證了LSMM中豐度的非負約束（Abundance Non-negativity Constraint， ANC）;然而，由于函數本身的非凸性，NMF具有大量的局部極小值，這將影響最終的解混結果。為了解決這個問題，需要在NMF模型中引入額外的結構性約束。例如，利用高光譜的空間特性[5-6]，將空間信息的約束引入到基于NMF的方法中，或者充分利用高光譜固有的稀疏性特性[7-8]，將稀疏性約束引入到基于NMF的解混方法中。近年來，一種新的稀疏約束NMF （L1/2-NMF）解混方法[9]在高光譜解混中得到了廣泛的應用。L1/2-NMF已經被證明是相對較優的方法。然而，與其他稀疏約束NMF一樣，L1/2-NMF的解混性能仍然受到初始設置和噪聲干擾的影響[10]。

考慮到稀疏約束NMF的高光譜解混方法存在的問題，本文采用經典的頂點成分分析法（Vertex Component Analysis，VCA）[11]得到端元，通過全約束最小二乘法（Fully Constrained Least Square， FCLS）[12]得到豐度，將該算法得到的端元和豐度矩陣作為算法輸入的初始值，同時在NMF中加入了正交性約束即正交非負矩陣分解（Orthogonal Non-negative Matrix Factorization， ONMF）[13-14]。通過進行模擬數據和真實數據實驗分析，稀疏和正交約束非負矩陣分解（Sparse and Orthogonal constrained Non-negative Matrix Factorization， SONMF）在對高光譜數據進行解混時能夠充分利用數據本身的特殊性，在保證端元相互獨立同時減小噪聲對解混結果的影響。本文的技術路線如圖1所示。

1 線性光譜混合模型

在線性光譜混合模型（LSMM）中，每一個端元的光譜特性是其端元光譜特性和豐度比例的線性組合，用矩陣表示為：

R = EA + N? （1）

s.t.? A ≥0， 1 Tp A = 1 Tn

式中： R =[ r i，…， r n]∈ R d×n表示高光譜圖像的n個像元以及d個波段; E =[ e 1， e 2，…， e p]∈ R d×p表示混合矩陣包含p個端元; A =[ a 1， a 2，…， a n]∈ R p×n表示端元的豐度矩陣; N 表示誤差矩陣。豐度的物理意義明確了豐度矩陣中的每個元素非負，且每一列的和為1。這兩個約束條件可以用以下的公式表示：

a i≥0; i=1，2，…，n

（2）

∑ n i=1? a i=1

（3）

2 非負矩陣分解與線性解混的結合

非負矩陣分解的數學模型是將高維的原矩陣分解為兩個低維非負矩陣，分解的過程實際上是一個不斷迭代優化的過程，建立一個合理的目標函數，再交替迭代求得分解后的非負矩陣，形如：

V ≈ WH

（4）

矩陣 V ∈ R m×n分解成兩個矩陣 W 和 H ，這里的 W ∈ R m×t， H ∈ R t×n，并且滿足t

非負矩陣分解的表達形式十分契合LSMM，所以可以將非負矩陣分解的方法應用于高光譜解混的條件，由于非負矩陣分解是個NP-hard問題，解決這類優化問題需要引入目標函數，文獻[4]中提出了一種基于歐氏距離的目標函數

J（ W ， H ）= 1 2 ‖ Y - WH ‖2F

（5）

其中，‖·‖2F表示Frobenius范數，因為目標函數是非凸的，所以尋找全局的最優解很困難的。因此，使用迭代方法如乘法更新規則被認為是解決這類基于NMF的問題的有效工具之一。式（5）的乘法更新法則可以表示為：

W ← W .*? VH T? WHH T? H ← H .*? W T V?? W T WH

（6）

2.1 基于稀疏約束的非負矩陣分解

在高光譜解混中，其中經常被用到的約束條件是基于豐度稀疏性的約束，稀疏性約束指的是每個混合像元中包含的端元數量遠小于總的端元數量，在傳統的稀疏約束的NMF（Sparse Non-negative Matrix Factorization， SNMF）解混方法中，L1 / 2正則化[9]作為稀疏約束項被引入NMF以追求更稀疏的豐度表示，由于L1/2-NMF能夠有效利用數據本身的稀疏性，與其他稀疏NMF方法相比，它顯示出較大的優勢，基于L1/2-NMF的目標函數可以表示為：

J（ W ， H ）= 1 2 ‖ Y - WH ‖2F+α‖ H ‖1/2

（7）

式中：‖ H ‖1/2=∑ L，N l，n=1 （ h l×n）1/2 表示正則化，α∈ R +表示正則化參數。

L1/2-NMF算法解混的效果受到端元以及豐度初始設置值的影響較大，另外，噪聲對解混結果也有很大的影響，所以為了取得更好解混效果，需要將更優的約束條件引入到NMF目標函數中去。

2.2 基于稀疏和正交約束的非負矩陣分解

由于高光譜解混提取的端元是相互獨立的，鑒于這樣的特殊性，本文將正交約束[13]應用到非負矩陣分解中，該目標函數可以表示為：

J（ W ， H ）= 1 2 ‖ Y - WH ‖2F （8）

s.t.? W ≥0， H ≥0， WW T= I

正交約束被證明等價于K-mean聚類等價，式（8）中對端元矩陣 W 加入正交約束，對提取出來的端元而言，保證了獨立性。考慮到高光譜數據集高維度的特點，通常基于正交非負矩陣分解算法復雜度為O（m2n），每當m增加時，相應的算法也會變得更加復雜，直接引入正交約束之后需要一定的先驗知識并且需要很大計算量，迭代求解的復雜度也隨之增加，所以將正交非負矩陣分解應用到高光譜解混中時需要降低算法復雜度。將正交約束直接表示在目標函數中，將正交約束作為目標函數的一部分，不需要借助額外的約束條件，這種方法可以一定程度上減少計算量，改進后的基于正交約束的NMF的目標函數可以表示為：

J（ W ， H ）= 1 2 ‖ Y - WH ‖2F+β‖ W T W - I ‖22 （9）

s.t.? W ≥0， H ≥0

式中：β表示正交回歸參數，‖ W T W - I ‖22表示2范數約束。

通過對式（5）引入端元正交性約束，可將豐度稀疏性和端元正交性約束非負矩陣分解（SONMF）的目標函數表示為：

J（ W ， H ）=? 1 2 ‖ Y - WH ‖2F+α‖ H ‖1/2+

β‖ W T W - I ‖22

（10）

s.t.? W ≥0， H ≥0， 1 Tp H = 1 Tn

對于式（10）的求解，這里使用乘法更新法則對矩陣 W 和 H 進行更新。首先對目標函數分別求關于 W 和 H 的偏導數，因為矩陣 W 和 H 是非負的，記

ξ W =? W?? WHH T+β WW T W

ξ H =? H?? W T WH + α 2? H - 1 2

（11）

將ξ W ，ξ H 代入：

W = W -ξ W?J?W

H = H -ξ H?J?H

（12）

使用最小二乘迭代規則，得到 W 和 H 的更新規則：

W ← W .*? VH T+2β W?? WHH T+2β WW T W

H ← H .*? W T V?? W T WH + α 2? H - 1 2

（13）

2.3 預處理環節

將SONMF算法應用于高光譜解混時，需要考慮一些預處理的環節。首先是端元數目的確定，本文采用最小誤差信號子空間識別（Hyperspectral Signal Identification by Minimum Error，HySime）[15]方法確定端元的數目;接著是端元和豐度初始化的問題，考慮到基于范數約束NMF受初始值影響較大，

用頂點成分分析法得到端元矩陣，通過全約束最小二乘法得到豐度矩陣，將得到的端元和豐度矩陣作為算法輸入的初始值，這樣不僅可以降低初始值對解混的影響，還可以加快算法的執行速度，減少算法的執行時間。關于參數設置的問題，稀疏正則化參數α主要根據矩陣豐度稀疏度來設置，具體的表達式如下：

λ= 1? L? ∑ l?? N -‖ x l‖1/‖ x l‖2? N -1

（14）

其中 x l表示高光譜數據第l維的光譜。正交化參數的設置根據文獻[14]設為0.05。

最后是停止條件設置的問題，本文停止算法的條件是迭代次數達到500時停止。

SONMF算法概括如下：

1）輸入圖像矩陣 Y ，用HySime算法估計端元的個數，采用VCA初始化端元矩陣 W ，利用FCLS初始化豐度矩陣 H ，設定參數α，β。

2）將初始化的矩陣 W 和 H 代入乘法更新式（13）進行更新。

3）如果目標函數式達到最大的迭代次數Tmax，算法停止運行，否則返回2）。

2.4 收斂性的證明

正交約束目標函數式（9）收斂性的證明參照了文獻[4]的思路，證明目標函數不遞增，考慮矩陣 W 中的某個元素wij，對式（9）求 J 對wij的一階偏導數和二階偏導數可得：

J′ij（w）= WHH T+2β WW T W - YH T-2β W

（15）

J″ij（w）=（ HH T）ij+2β（（ W T W ）jj+w2ii-1）

（16）

構造關于 J 的輔助函數，記

F（w，wtij）=Jij（wtij）+J′ij（wtij）（w-wtij）+

（ WHH T）ij+2β（ WW T W ）ij wtij ?（w-wtij）2

（17）

證明收斂性只要證明F（w，wtij）≥Jij（w）即可。

F（w，wtij）=Jij（wtij）+J′ij（wtij）（w-wtij）+

（ WHH T）ij+2β（ WW T W ）ij wtij ?（w-wtij）2

（18）

根據條件（ W T W ）jj=1， Jij（w） 可以表示為：

Jij（w）= Jij（wtij）+J′ij（w-wtij）+

[（ HH T）ij+2βw2ii]（w-wtij）2

（19）

同時可以得到：

（ WHH T）ij=∑ k p=1 wtip（ HH T）pj≥wtij（ HH T）jj

（20）

（ WW TW）ij=∑ k p=1 wip（ W T W ）pj≥wii（ W T W ）≥???? wii （ ∑ k p=1 wpiwpj ） ≥wiiwiiwij=w2iiwij

（21）

通過比較式（17）和（19）得到F（w，wtij）≥Jij（w），所以原目標函數式是收斂的。

3 實驗分析

3.1 精度評價指標

衡量高光譜解混的精度通常通過兩種評價指標：光譜角距離（Spectral Angle Distance，SAD）和均方根誤差（Root Mean Square Error，RMSE）。

SAD的計算公式如下：

RSAD（ a ， b ）=arccos ??a T b? ‖ a ‖2·‖ b ‖2

（22）

其中： a 表示提取的端元光譜矢量， b 表示參考光譜矢量即光譜庫中的已知光譜矢量。

RMSE的計算公式如下：

RRMSE= 1 n ∑ n i=1?? 1 L ‖ r i-?? i‖22

（23）

其中：{ ??i}ni=1表示重構的圖像。RMSE的值越小表示原始圖像與重構圖像間的誤差越小，解混的精度越高。

3.2 合成數據分析

本實驗在處理器i5，CPU頻率3.5GHz，內存8GB的計算機上使用Matlab2014a操作平臺進行實驗，合成數據是從美國地質勘探局USGS光譜中選取五個端元Alunite、Buddingtonit、Calcite、Kaolinite和Muscovite，它們的光譜包含224個波段，波長范圍為0.38～2.5μm，光譜圖見圖2。

將這五個端元以Dirichlet分布的形式進行混合，同時對端元豐度的和進行歸一化操作，為了模擬數據采集過程，按照指定的信噪比（Signal-to-Noise Ratio，SNR）在模擬圖像中加入不同的噪聲等級，其中SNR可表示為：

SNR=10 lg E[ x T x ] E[ n T n ]

（24）

式中： x 和 n 表示原始信號和相應的噪聲，E[·]表示期望值。

實驗比較了VCA[11]、ONMF[13]、SNMF[9]以及本文采用的SONMF算法的性能。

實驗1

不同噪聲條件下算法魯棒性分析。本實驗研究高斯噪聲對實驗的影響，所以保持圖像大小為100×100，端元數目為5，SNR從15 dB到35dB，間隔為5dB。

實驗2

不同像元數目條件下對實驗的影響。本實驗研究像元數目對實驗的影響，所以保持端元數目為5，SNR為20dB的條件不變，圖像大小分別為40×40，60×60，80×80，100×100。

實驗結果如圖3、4所示，可以看出，SONMF解混效果整體要優于SNMF和ONMF的解混效果。

本文采用的數據集是1997年機載可見光及紅外光譜成像儀采集到的Urban地區高光譜圖像（307×307像素），總共有210個波段。除去低信噪比和水汽吸收的波段（1-4，76，87，101-111，136-153，198-210），只有162個有效波段，該地區包含了屋頂、樹木、草地、道路等多種地物的混合，圖5顯示了四種地物光譜的真實豐度圖。

首先通過HySime估計端元數目為4，再使用VCA算法提取的端元及FCLS估計的豐度作為初始值，分別使用SNMF、ONMF和SONMF算法對該地區的高光譜圖像進行處理，表1給出了不同算法最終的SAD的標準值和均值，加粗的字體是提取較好的結果，本文提出的SONMF算法的解混結果是比較理想的，均值SAD降低0.012～0.145。

圖6給出了這四種算法對Urban地區高光譜數據進行解混得到的端元對應的豐度圖。

將圖6中提取的光譜豐度圖與圖5中地物光譜的真實端元豐度圖進行對照可以看出，除了VCA算法分解出的端元豐度圖效果較差之外，其他三種算法分解出的端元豐度圖分布在道路和草地部分與真實端元豐度圖分布比較接近;然而在屋頂和樹木部分，SONMF算法分解出端元豐度圖明顯要優于其他的算法。

4 結語

混合像元處理是高光譜遙感領域的一個熱門研究，本文采用線性光譜混合模型來解決光譜分解的問題，重點討論了非負矩陣分解的方法，非負矩陣分解沒有獨立的解決方案，需要在目標函數中結合其他的約束條件求解，本文將稀疏性與端元正交性約束結合作為非負矩陣分解的約束條件，同時在迭代初始值方面作了一些改進。本文算法在保證提取端元的獨立性和稀疏端元本身豐度的同時，減小了噪聲和局部極小值對解混解混結果的影響，取得了良好的解混效果。但本文算法相對比較復雜，如何進一步通過預處理的手段降低算法復雜度是下一步的研究工作。

參考文獻

[1]?童慶禧，張兵，鄭蘭芬.高光譜遙感——原理、技術與應用[M].北京：高等教育出版社， 2006：15-30. （TONG Q X， ZAHNG B， ZHENG L F. Hyperspectral Remote Sensing——Principle， Technology and Application [M]. Beijing： Higher Education Press， 2006： 15-30.）

[2]?BIOUCAS-DIAS J M， PLAZA A， DOBIGEON N， et al. Hyperspectral unmixing overview： geometrical， statistical， and sparse regression-based approaches [C]// Proceedings of the 2011 IEEE International Geoscience and Remote Sensing Symposium. Piscataway， NJ： IEEE， 2011： 1135-1138.

[3]?張兵，孫旭.高光譜圖像混合像元分解[M].北京：科學出版社， 2015：6-22. （ZAHNG B， SUN X. Hyperspectral Image Unmixing [M]. Beijing： Science Press， 2015： 6-22.）

[4]??LEE D D， SEUNG S S. Learning the parts of objects by non-negative matrix factorization [J]. Nature， 1999， 401： 788-791.

[5]?袁博.馬爾可夫隨機場的空間相關模型在非負矩陣分解線性解混中的應用[J].計算機應用，2017，37（12）：3563-3568. （YUAN B. Application of MRFs spatial correlation model in NMF-based linear unmixing [J]. Journal of Computer Applications， 2017， 37（12）： 3563-3568.）

[6]?WANG X， ZHONG Y， ZHANG L， et al. Spatial group sparsity regularized nonnegative matrix factorization for hyperspectral unmixing [J]. IEEE Transactions on Geoscience and Remote Sensing， 2017， 55（11）：6287-6304.

[7]?RAJABI R， GHASSEMIAN H. Spectral unmixing of hyperspectral imagery using multilayer NMF [J]. IEEE Geoscience and Remote Sensing Letters， 2014， 12（1）： 38-42.

[8]?ZHAO Y， ZHOU Z， WANG D， et al. Hyperspectral image unmixing algorithm based on endmember-constrained nonnegative matrix factorization [J]. Frontiers of Optoelectronics， 2016， 9（4）： 627-632.

[9]??QIAN Y， JIA S， ZHOU J， et al. Hyperspectral unmixing via L1/2 sparsity-constrained nonnegative matrix factorization [J]. IEEE Transactions on Geoscience & Remote Sensing， 2011， 49（11）： 4282-4297.

[10]?余肖玲.非負矩陣分解理論及其在高光譜解混中的應用[D]. 成都：成都理工大學， 2015：45-60. （YU X L. The theory of Nonnegative matrix factorization and its application in the hyperspectral unmixing [D]. Chengdu： Chengdu University of Technology， 2015： 45-60.）

[11]?NASCIMENTO J M P，BIOOUCASDIAS J M B. Vertex component analysis： a fast algorithm to unmix hyperspectral data [J]. IEEE Transactions on Geoscience and Remote Sensing， 2005， 43（4）： 898-910.

[12]?HEINZ D C， CHANG C. Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery [J]. IEEE Transactions on Geoscience and Remote Sensing， 2002， 39（3）： 529-545.

[13]?POMPILI F， GILLIS N， ABSIL P， et al. Two algorithms for orthogonal nonnegative matrix factorization with application to clustering [J]. Neurocomputing， 2014， 141：15-25.

[14]?李孟杰，謝強，丁秋林.基于正交非負矩陣分解的K-means聚類算法研究[J].計算機科學，2016，43（5）：204-208. （LI M J， XIE Q， DING Q L. Orthogonal Non-negative Matrix Factorization for K-means Clustering [J]. Computer Science， 2016， 43（5）： 204-208.）

[15]?BIOUCASDIAS J M， NASCIMENTO J M P. Hyperspectral subspace identification [J]. IEEE Transactions on Geoscience and Remote Sensing， 2008， 46（8）：2435-2445.