孫莉,于瑞林,吳杰芳

非負矩陣分解與光譜解混

孫莉1,于瑞林1,吳杰芳2

1. 山東農(nóng)業(yè)大學 信息科學與工程學院, 山東 泰安 271018 2. 泰山學院 數(shù)學與統(tǒng)計學院, 山東 泰安 271000

非負矩陣分解(NMF)用兩個非負矩陣的乘積近似原始數(shù)據(jù)對應(yīng)的非負矩陣，它為基于線性光譜混合模型的光譜解混提供了新途徑。給出NMF在光譜解混中三個矩陣的具體含義后，用五種求解NMF的有效算法，對Jasper Ridge的高光譜遙感圖像進行解混。討論了五種算法的迭代方式以及收斂性質(zhì)。實驗結(jié)果表明，五種算法能成功分離出4種端元光譜以及相應(yīng)的豐度譜圖，其中有效集型算法表現(xiàn)突出。

非負矩陣分解; 光譜解混; 界約束優(yōu)化; 有效集

非負矩陣分解(Nonnegative matrix factorization, NMF)的思想可以追溯到1994年P(guān)attero和Tapper的文章[1]，隨后1999年，Lee DD和Seung HS獨立的提出了NMF的概念后[2]，該方法迅速引起了眾多學者的關(guān)注，蓬勃發(fā)展起來。用兩個非負矩陣的乘積近似原始數(shù)據(jù)對應(yīng)的非負矩陣，非負性的要求保證了通過部分表示原始數(shù)據(jù)的形式，即樣本數(shù)據(jù)只允許加性組合，這使得分解結(jié)果對原始數(shù)據(jù)的特征及結(jié)構(gòu)具有相當?shù)谋磉_能力，有利于發(fā)現(xiàn)數(shù)據(jù)的局部特征。

1 NMF的基本思想與主要求解算法

1.1 NMF的基本思想

若的各列表示每一樣本的相關(guān)屬性，則經(jīng)NMF分解得到的矩陣稱為基矩陣，稱為系數(shù)矩陣。將和按列分塊，即=(1,2,…,w)，=(1,2,…,v)，由(1)可得個等式，

w≈Uv，=1,2,…,(2)

(2)式說明，每一個樣本w可近似的看作非負矩陣的列向量的非負線性組合，組合系數(shù)是v的分量。于是，矩陣可以看作是對原始矩陣進行線性逼近的一組基，而則是樣本集在基上的非負線性投影系數(shù)。尋找合適的基向量組，用其表征大量的數(shù)據(jù)向量，一方面可以給出數(shù)據(jù)之間潛在的結(jié)構(gòu)關(guān)系，另一方面實現(xiàn)了數(shù)據(jù)壓縮，用少量的數(shù)據(jù)實現(xiàn)對原始數(shù)據(jù)很好的逼近。通過NMF對矩陣進行分解后，所得矩陣和在不同的應(yīng)用領(lǐng)域有不同的含義，它在圖像分析、文本聚類、語音處理、生物醫(yī)學等領(lǐng)域應(yīng)用廣泛。本文第二部分，我們會給出NMF應(yīng)用于光譜解混時，和的具體含義。

1.2 NMF問題模型

求解矩陣和的常用方法是通過最小化和的差距，即

1.3 求解NMF問題的有效算法

針對NMF問題模型，已經(jīng)提出了多種行之有效的求解方法，現(xiàn)有算法可分為兩類：第一類，基于梯度下降的交替迭代公式，通過公式交替實現(xiàn)和的更新。目前求解NMF的常用方法——乘法迭代公式便屬于這一類。算法的特點是結(jié)構(gòu)簡單，易于實現(xiàn)，但其在運行效率與收斂性分析上有所欠缺。第二類，交替非負最小二乘(Alternating Nonnegative Least Squares, ANLS)框架下的有效算法，ANLS算法框架如下：

ANLS算法框架

步驟 1 給定初始矩陣0和0，:=0。

步驟 2 依次實現(xiàn)和的更新如下:

基于ANLS算法框架下的優(yōu)化方法有投影梯度法、投影(擬)牛頓法、有效集型方法和塊轉(zhuǎn)軸方法等。2007年Lin CJ的文獻[3]是利用傳統(tǒng)界約束優(yōu)化算法求解NMF問題的開端。以優(yōu)化理論為基礎(chǔ)，這些算法具有良好的收斂性質(zhì)。

下面給出5種求解NMF的算法。

（1）乘法迭代公式(MU)

該方法一經(jīng)提出，便成為求解NMF問題的標準算法[4]，其迭代規(guī)則為：

（2）交替最小二乘算法(ALS)

通過交替求解最小二乘問題實現(xiàn)和的更新，一旦矩陣元素不滿足非負性要求，直接將該元素定義為0。忽略非負性約束，無法保證目標函數(shù)的單調(diào)下降，因此無法獲得算法的收斂性結(jié)論，但這一算法運行速度快，其迭代規(guī)則如下：

和已知，通過UUV=UW，求出，將中的負數(shù)直接賦值為0；

和已知，通過VVU=VW，求出，將中的負數(shù)直接賦值為0。

（3）投影梯度法 (PGM)

（4）投影牛頓法 (PNM)

Gong PH和Zhang CS以Lin CJ的投影梯度法為基礎(chǔ)給出了投影牛頓法[5]，它通過乘子信息同樣將變量分為有效變量和非有效變量，以投影牛頓方向為搜索方向?qū)崿F(xiàn)非有效變量的更新。這一方法的優(yōu)點在于一步迭代允許有多個指標進出有效集，且搜索方向的定義借助于目標函數(shù)的二階信息，它獲得了二階收斂的速度。

（5）有效集型算法 (AS)

Kim H和Park H提出了求解NMF的有效集型算法[6]，他們將待求變量分為兩類：有效變量(非負約束中等式成立的變量)和非有效變量，對非有效變量通過求解無約束的最小二乘問題實現(xiàn)更新。隨后以目標函數(shù)下降為前提，重新制定有效變量和非有效變量。算法是否停止迭代，通過檢測最優(yōu)化條件判斷，即非有效變量的非負性和有效變量對應(yīng)乘子的非負性是否滿足。這一算法的主要問題在于每步迭代僅允許一個指標進出有效集，若初始矩陣的選取不合適，算法的計算復(fù)雜度很高。

2 光譜解混中的NMF

2.1 遙感圖像的矩陣表示

高光譜圖像可以看成是三維數(shù)據(jù)構(gòu)成的立方體，它的和維表現(xiàn)的是地球表面的坐標，第三維是波段，它由光波的頻率所決定。例如，將一幅AVIRIS圖像讀取為××的矩陣，其中×對應(yīng)圖像的寬度和高度，對應(yīng)波段數(shù)，于是一幅AVIRIS圖像對應(yīng)了幅×的二維圖像，每幅圖像對應(yīng)地球上相同位置在不同波長下的影像。不妨設(shè)圖像的像素點有個，=×，于是將上述矩陣重新排列可得到×的矩陣。的行向量代表某一像素點在不同波段的灰度值(DN值)，而的列向量對應(yīng)某一波段下，個像素點的DN值。

2.2 遙感圖像的矩陣表示

通常遙感影像以像元為單位獲取地物信息，所記錄的是像元內(nèi)所有地物光譜的綜合。在實際中，由于空間分辨率的限制以及地物的復(fù)雜多樣性，一個像元內(nèi)往往會包含多種地物類型，稱為混合像元。為提高分類結(jié)果對真實地表覆蓋的描述準確性，需要進入像元內(nèi)部，將混合像元分解為不同的“基本組成單元”或“端元”，并求得這些基本組分所占的比例，這就是所謂的“光譜解混”過程。

在線性光譜混合模型中，混合光譜是端元光譜與其比例的線性組合，其數(shù)學表達式如下：

其中，表示端元的個數(shù)，為維光譜向量對應(yīng)圖像中某一像元在個波段的DN值，1,2,…,e分別對應(yīng)個端元對應(yīng)的光譜向量，c代表像元中端元e所占的比例，為誤差項。

2.3 遙感圖像的矩陣表示

非負矩陣分解與光譜解混中線性光譜混合模型相吻合，為光譜解混提供了一種新的途徑。利用NMF實現(xiàn)的分解后，得到的矩陣為豐度矩陣，它的行給出某一像元中各端元所占比例，列向量則對應(yīng)某個端元的空間分布。為端元矩陣，它的行向量為某個端元的光譜向量[7,8]。

3 數(shù)值比較

本文所采用的數(shù)據(jù)來自于2006年機載可見光及紅外成像光譜儀(AVIRIS)采集到的Jasper Ridge高光譜圖像，圖像大小為512×614像素，自380 nm至2500 nm，有224個波段，光譜分辨率達到9.46 nm。我們截取了100×100像素的子圖，該區(qū)域包含了樹木、道路、裸土等多種地物的混合，去除信噪比低和水蒸氣吸收波段(1~3，108~112，154~166以及220~224)，余下198個有效波段。圖1為第35，20，7波段圖像合成的偽彩色圖像。

圖1 Jasper Ridge的AVIRIS偽彩色圖

Jasper Ridge地區(qū)的樹木、水體、裸土以及道路，4種地物的真實豐度譜圖，如圖2所示。

圖2 Jasper Ridge真實地物豐度譜圖

設(shè)置端元數(shù)目為4，即=4，利用NMF實現(xiàn)上述遙感圖像的光譜解混，算法程序由Matlab 7.0編寫。我們主要從解混后的豐度譜圖，分解所得的端元光譜與真實端元光譜的余弦值以及運行時間3個方面比較5種算法的性能。求解NMF的5種算法分別為乘法迭代公式(MU)，交替最小二乘算法(ALS)，投影梯度法(PGM)，投影牛頓法(PNM)以及有效集型算法(AS)。

NMF問題(2)非凸，導(dǎo)致算法的運行效率與初始點的選取密切相關(guān)。數(shù)值測試中，利用(,)，(,)隨機生成初始矩陣0和0，隨后5種算法以相同的0，0開始運行。算法的最高迭代次數(shù)設(shè)置為1000，運行時間上限為104s，PGM和PNM增加了下述終止準則，即：

下標表示的第個分量。

圖3 采用5種算法獲得的Jasper Ridge四地物豐度譜圖比較

圖3為采用5種算法對圖像數(shù)據(jù)矩陣進行NMF分解，所得豐度矩陣對應(yīng)的豐度譜圖，圖中第1至4列分別對應(yīng)樹木、水體、裸土以及道路。通過對比發(fā)現(xiàn)，5種算法均成功獲得樹木與裸土的豐度譜圖。其中MU，PGM，PNM，AS所得樹木的豐度譜圖與真實地物吻合度較高，同時，裸土豐度譜圖以ALS和AS算法所得為佳。水體與道路的分離不顯著，如圖3所示，第二列圖像對應(yīng)的水體豐度譜圖中，水體與道路的信息一并顯示出來。通過NMF分解所得的端元矩陣，提供了四種端元：樹木，水體，裸土以及道路的端元光譜信息。下面計算分離的光譜信息與美國地質(zhì)調(diào)查局(USGS)庫中的參考光譜間的光譜角距離。光譜角距離(Spectral Angel Distance, SAD)，即兩個光譜向量之間的夾角，計算公式如下：

其中V，V表示不同端元的光譜向量，通過SAD計算出的角度越小代表與參考光譜越匹配。

表4 USGS譜庫的真實光譜和提取的端元光譜的SAD比較

由表1看出AS算法提取出的4個端元光譜與真實光譜最接近，PGM和PNM算法所得的豐度譜圖以及4個端元光譜信息十分接近。事實上，PGM和PNM算法的終止準則以及算法框架大致相同，主要區(qū)別在于搜索方向的定義，PGM借助于投影梯度方向，而PNM則建立在投影牛頓方向的基礎(chǔ)上。最后，利用5種算法完成光譜解混工作所需的計算時間如表2所示，時間單位：秒。

表5 不同算法完成光譜解混的時間比較

PNM采用了借助二階信息的搜索方向，計算時間較之PGM明顯減少。5種算法中AS算法的運行時間最少，經(jīng)分析，AS算法在兩種情形會終止迭代，一是超出最大迭代次數(shù)；二是滿足一階最優(yōu)性條件，即非有效變量非負性以及有效變量對應(yīng)乘子的非負性滿足。本測試在第二種情形終止迭代。

4 小結(jié)

NMF應(yīng)用于光譜解混工作，在給定端元個數(shù)的條件下，可同時提取出端元光譜和每個端元的豐度譜圖。本文總結(jié)了5種求解NMF的有效方法，并具體應(yīng)用于高光譜圖像的解混。結(jié)果表明，有效集型算法的分解結(jié)果與真實地物的吻合度較高，且運行時間有明顯優(yōu)勢。我們將進一步將有效集識別方法[9]推廣至正則化NMF問題[7,8]的求解。

[1] Paatero P, Tapper U. Positive matrix factorization: a non-negative factor model with optimal utilization of error estimates of data values[J]. Environmetricsa, 1994,5(2):111-126

[2] Lee DD, Seung HS. Learning the parts of objects by non-negative matrix factorization[J]. Nature, 1999,401(6755):788-791

[3] Lin CJ. Projected gradient methods for non negative matrix factorization[J]. Neural Comput.,2007,19(10):2756-2779

[4] Lee DD, Seung HS. Algorithms for non-negative matrix factorization[C]//Neural Information Processing Systems Conference. Massachusetts USA: MIT Press, 2001:556-562

[5] Gong PH, Zhang CS. Efficient nonnegative matrix factorization via projected Newton method[J]. Pattern Recognition, 2012,45(9):3557-3565

[6] Kim H, Park H. Nonnegative matrix factorization based on alternating nonnegativity constrained least squares and active set method[J]. SIAM Journal on Matrix Analysis and Applications, 2008,30(2):713-730

[7] Gillis N, Plemmons RJ. Sparse nonnegative matrix underapproximation and its application to hyperspectral image analysis[J]. Linear Algebra and its Applications, 2013,438(10):3991-4007

[8] 孫莉,趙庚星.基于正交非負矩陣分解的高光譜遙感圖像混合像元分解[J].山東省農(nóng)業(yè)大學學報(自然科學版),2017,48(2):264-267

[9] 孫莉,賀國平,房亮.基于求解大規(guī)模界約束問題的三種有效集識別策略的比較[J].數(shù)值計算與計算機應(yīng)用,2009,30(1):41-47

Non-negative Matrix Factorization and Hyperspectral Unmixing

SUN Li1, YU Rui-lin1, WU Jie-fang2

1.271018,2.271000,

Nonnegative matrix factorization approximates the data matrix with the production of two low-rank matrices. It provides a new approach to hyperspectral image analysis, which is based on the standard linear mixing model. After introducing the implications of the three matrices in NMF, we summarize five different computational algorithms which solve NMF successfully, and their iteration procedures and convergence properties are analyzed. In the experimental test, five algorithms are employed to unmix the hyperspectral image Jasper Ridge. Numerical results show that, they can detect the endmembers spectral signatures and the corresponding fractional abundance. A good choice is the active set type method which performs noticeably well.

Nonnegative matrix factorization; hyperspectral unmixing; bound constrained optimization; active set

O221;TP751.1

A

1000-2324(2019)05-0908-05

10.3969/j.issn.1000-2324.2019.05.038

2015-09-05

2015-10-29

國家自然科學基金項目(11701337,10901094,11301307);山東省自然科學基金項目(ZR2016DM03);山東省高等學校科技計劃項目(J16LI16);大學生創(chuàng)新創(chuàng)業(yè)訓練計劃項目(201410434097)

孫莉(1980-),女,博士,副教授,研究方向為最優(yōu)化算法與理論. E-mail:sunlishi@hotmail.com