非負(fù)矩陣分解與光譜解混

孫莉,于瑞林,吳杰芳

非負(fù)矩陣分解與光譜解混

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

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

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

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

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

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

1.1 NMF的基本思想

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

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

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

1.2 NMF問題模型

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

1.3 求解NMF問題的有效算法

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

ANLS算法框架

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

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

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

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

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

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

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

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

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

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

（3）投影梯度法 (PGM)

（4）投影牛頓法 (PNM)

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

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

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

2 光譜解混中的NMF

2.1 遙感圖像的矩陣表示

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

2.2 遙感圖像的矩陣表示

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

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

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

2.3 遙感圖像的矩陣表示

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

3 數(shù)值比較

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

圖1 Jasper Ridge的AVIRIS偽彩色圖

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

圖2 Jasper Ridge真實(shí)地物豐度譜圖

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

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

下標(biāo)表示的第個(gè)分量。

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

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

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

表4 USGS譜庫(kù)的真實(shí)光譜和提取的端元光譜的SAD比較

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

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

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

4 小結(jié)

NMF應(yīng)用于光譜解混工作，在給定端元個(gè)數(shù)的條件下，可同時(shí)提取出端元光譜和每個(gè)端元的豐度譜圖。本文總結(jié)了5種求解NMF的有效方法，并具體應(yīng)用于高光譜圖像的解混。結(jié)果表明，有效集型算法的分解結(jié)果與真實(shí)地物的吻合度較高，且運(yùn)行時(shí)間有明顯優(yōu)勢(shì)。我們將進(jìn)一步將有效集識(shí)別方法[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] 孫莉,趙庚星.基于正交非負(fù)矩陣分解的高光譜遙感圖像混合像元分解[J].山東省農(nóng)業(yè)大學(xué)學(xué)報(bào)(自然科學(xué)版),2017,48(2):264-267

[9] 孫莉,賀國(guó)平,房亮.基于求解大規(guī)模界約束問題的三種有效集識(shí)別策略的比較[J].數(shù)值計(jì)算與計(jì)算機(jī)應(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

國(guó)家自然科學(xué)基金項(xiàng)目(11701337,10901094,11301307);山東省自然科學(xué)基金項(xiàng)目(ZR2016DM03);山東省高等學(xué)?？萍加?jì)劃項(xiàng)目(J16LI16);大學(xué)生創(chuàng)新創(chuàng)業(yè)訓(xùn)練計(jì)劃項(xiàng)目(201410434097)

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