利用約束非負(fù)矩陣分解的高光譜解混算法

趙春暉，成寶芝，楊偉超

(哈爾濱工程大學(xué)信息與通信工程學(xué)院，黑龍江哈爾濱150001)

高光譜遙感數(shù)據(jù)是通過高光譜遙感測量儀器定量獲取的，是由數(shù)十至數(shù)百個波長在0.3～2.5μm之間的窄波段(波段寬度小于10 nm)形成的像元組成，可以提供幾乎連續(xù)的地物光譜曲線[1]，這些地物光譜曲線相比于多光譜來說更全面反映了成像的目標(biāo)特征.但是，由于遙感儀器中使用的傳感器空間分辨率受到技術(shù)條件的限制和觀測的地物情況的復(fù)雜性，使得高光譜圖像存在著混合像元的情況.如果把混合像元作為純像元進(jìn)行分類、目標(biāo)探測等應(yīng)用研究，結(jié)果會有很大的誤差.這使高光譜解混問題成為近年來遙感領(lǐng)域的一個研究熱點(diǎn).目前，利用非負(fù)矩陣分解進(jìn)行高光譜解混得到充分關(guān)注，Miao等[2]提出的最小體積約束的非負(fù)矩陣分解(minimum volume constrained nonnegative matrix factorization， MVC-NMF)算法;Alexis Huck 等［3］提出的結(jié)合最小離差約束的非負(fù)矩陣分解(minimum dispersion constrained nonnegative matrix factorization，MIDINMF)算法;Jia S等［4］提出的基于稀疏性和分段平滑性約束的非負(fù)矩陣分解;YU Yue等［5］提出的采用最小距離約束的非負(fù)矩陣分解(minimum distance constrained nonnegative matrix factorization，MDNMF)算法，吳波等［6］提出的基于端元約束的非負(fù)矩陣分解(constrained nonnegative matrix factorization，CNMF)算法.由于高光譜數(shù)據(jù)的復(fù)雜性，這些改進(jìn)的約束非負(fù)矩陣分解算法都有一定的局限性.

本文通過對高光譜混合像元含有的端元光譜的光譜特性和豐度分布特性的深入分析，提出了一個新的約束非負(fù)矩陣分解算法，即以最小估計豐度協(xié)方差和單形體各頂點(diǎn)到中心點(diǎn)均方距離總和最小為約束條件的非負(fù)矩陣分解(minimum covariance and minimum distances nonnegative matrix factorization，MCMDNMF)算法.這個新的算法充分考慮了混合像元的光譜特性和空間特性，沒有對端元構(gòu)成的單形體體積、空間分布的稀疏性和平滑性進(jìn)行約束，避免了上述算法的缺陷，也不需要有純像元存在這一先驗(yàn)條件.

1 線性光譜混合模型

在分析高光譜含有的混合像元時，一般應(yīng)用線性光譜混合模型［7］進(jìn)行分析，混合像元是由各類端元和其對應(yīng)的豐度線性混合組成.假設(shè)U是波段為的高光譜圖像，M是一個L×P的光譜特征矩陣，α=［α1，…，αP］T為端元列向量對應(yīng)的豐度向量，M=［m1，…，mP］為 P 個端元向量，線性混合模型可以寫成

式中:n為一個L維的噪聲或者誤差;P個端元向量和對應(yīng)的豐度都是未知量，它們需要滿足2個約束條件:1)端元光譜及其豐度是非負(fù)的，即 mi≥0，0≤αi≤1;2)豐度總和為 1，即…，P.

2 約束非負(fù)矩陣分解的解混方法

非負(fù)矩陣分解［8］是由Lee等提出的，分解后的矩陣所有分量均為非負(fù)值，是一種新的矩陣分解方法.非負(fù)矩陣分解能使原始數(shù)據(jù)結(jié)構(gòu)得到清晰的展示，使高維的原始數(shù)據(jù)得到一定程度的維數(shù)約減.非負(fù)矩陣分解模型和基于線性光譜混合模型的混合像元分解相符合，非常適合于高光譜數(shù)據(jù)解混問題.但是由于原始非負(fù)矩陣分解方法目標(biāo)函數(shù)的非凸性，導(dǎo)致非負(fù)矩陣分解在某些情況下僅能得到局部最優(yōu)解，所以目標(biāo)函數(shù)需要增加約束條件.

本文采用一種新的約束非負(fù)矩陣分解方法，即MCMDNMF算法進(jìn)行背景信息端元的提取和混合像元分解.以最小估計豐度協(xié)方差和單形體各頂點(diǎn)到中心點(diǎn)均方距離總和最小這2個條件作為約束，下面分析一下這2個約束條件.

根據(jù)混合像元的光譜特性和空間分布結(jié)構(gòu)特點(diǎn)，由式(1)可知，利用最小二乘法可以得到豐度向量α無約束解的估計為

式中:φ2為誤差或者噪聲n的方差，一般設(shè)定n為零均值的高斯白噪聲.

式中:tr(·)表示矩陣的跡.如果將式(4)作為MCMDNMF的目標(biāo)函數(shù)，那么在迭代規(guī)則計算過程中存在微分過程求解復(fù)雜的問題.根據(jù)文獻(xiàn)［6］的分析，對式(4)進(jìn)行修正后作為約束項(xiàng):

僅式(5)作為非負(fù)矩陣分解的約束條件還不能得到理想的解混結(jié)果，因此，需要再增加約束條件提高解混精度.通過對端元光譜分布的單形體進(jìn)行分析，引入單形體各頂點(diǎn)到中心點(diǎn)距離總和作為另一個約束條件[5，12]，即

因此，以J1(M)和J2(M)為約束項(xiàng)的非負(fù)矩陣分解的目標(biāo)函數(shù)為

若將式(7)作為改進(jìn)的MCMDNMF的目標(biāo)函數(shù)，沒有考慮高光譜豐度總和為1這一約束條件，而這個條件作為軟約束可以使解混的結(jié)果更精確.本文采用文獻(xiàn)[13]中所使用的表示形式，即式中:δ為權(quán)值，控制總和為1這一約束條件對目標(biāo)函數(shù)的影響，1T表示一個向量的所有分量都為1，即當(dāng)?shù)蠼猞習(xí)r，可用代替U代替M.

MCMDNMF采用迭代算法計算W和H的最優(yōu)解:

依據(jù)式(9)，需要對J1(M)和J2(M)進(jìn)行求導(dǎo)，可得

因此，目標(biāo)函數(shù)的梯度為

由式(9)、(12)、(13)可得

根據(jù)式(14)和(15)可以求得M和α的最優(yōu)解，也就是解混所要得到的最優(yōu)端元和對應(yīng)的豐度.

3 仿真實(shí)驗(yàn)與分析

利用合成的高光譜數(shù)據(jù)和真實(shí)的高光譜圖像進(jìn)行仿真實(shí)驗(yàn)和分析.仿真實(shí)驗(yàn)是為了驗(yàn)證提出的MCMDNMF算法的有效性，并與前文提到的MVCNMF、CNMF和MDNMF幾種典型算法進(jìn)行性能比較.

3.1 性能評價指標(biāo)

光譜角距離(spectral angel distance，SAD)[14]和均方根誤差(root mean square error，RMSE)[15]是評價高光譜數(shù)據(jù)解混結(jié)果性能好壞的常用指標(biāo).光譜角距離用來比較真實(shí)端元光譜M及其估計值M^的相似性，如下:

均方根誤差用來度量端元對應(yīng)的真實(shí)豐度α和估計的豐度α^之間的差異，如下:

3.2 合成高光譜數(shù)據(jù)

從USGS光譜數(shù)據(jù)庫[16]提取了6種線性獨(dú)立的端元光譜，如圖1所示，這些端元光譜具有224個光譜波段，反射率的波長范圍在 0.38～2.5μm 之間，6個端元光譜按Dirichlet分布進(jìn)行混合，形成相應(yīng)的豐度，保證豐度非負(fù)且總和為1的約束條件.對MCMDNMF算法的抗噪聲能力和空間維數(shù)不同時算法的解混性能進(jìn)行了仿真實(shí)驗(yàn)，并與其他幾種算法在相同仿真環(huán)境下的性能情況進(jìn)行了比較.

圖1 合成高光譜數(shù)據(jù)的端元光譜圖Fig.1 Endmember spectral of synthetic data

3.2.1 算法的抗噪聲能力

驗(yàn)證在相同的混合光譜，不同信噪比的情況下MCMDNMF算法和其他幾種算法的抗噪聲能力.仿真設(shè)定的像元數(shù)目為10 000，分別對信噪比為15、20、25、30、35、∞ dB(假設(shè)合成光譜數(shù)據(jù)不含噪聲)進(jìn)行仿真實(shí)驗(yàn)，使用SAD和RMSE指標(biāo)進(jìn)行性能評價.仿真時用的像元是通過Dirichlet分布混合而成的，并且每次生成的噪聲信號也是隨機(jī)產(chǎn)生的.所以，SAD和RMSE求取平均值評估算法的性能，故圖2中的評價指標(biāo)值都為平均值，分別用和RMSE表示.通過比較可以發(fā)現(xiàn)，MCMDNMF算法和MVCNMF算法性能相近，在低信噪比時，性能優(yōu)于MVCNMF;高信噪比時，性能比MVCNMF稍差.MCMDNMF算法性能明顯要優(yōu)于CNMF和MDNMF算法.

圖2 信噪比不同時的算法性能比較Fig.2 Comparison of the algorithms with different SNR levels

3.2.2 空間維數(shù)不同時算法的解混能力

驗(yàn)證在相同的混合光譜，空間維數(shù)不同時MCMDNMF算法和其他幾種算法的解混性能和魯棒性，空間維數(shù)不同也就是混合像元數(shù)目發(fā)生變化.因此，在像元數(shù)目為500、1 000、5 000、8 000、10 000，信噪比為 35 dB的情況下，對算法進(jìn)行了仿真比較，如圖3所示.

通過比較可以發(fā)現(xiàn)，從整體上來說，MCMDNMF算法和MVCNMF算法性能相近，空間維數(shù)越高，MCMDNMF性能越優(yōu);MCMDNMF算法性能明顯要優(yōu)于CNMF和MDNMF算法.

3.3 真實(shí)的高光譜圖像

仿真用的數(shù)據(jù)來自于1997年AVIRIS(airborne visible/infrared imaging spectrometer)在美國內(nèi)華達(dá)州采集到的Cuprite區(qū)域的一個實(shí)際的高光譜圖像，圖像包含著多種礦物[17-18].去除了低信噪比和水蒸氣吸收波段(為1～2，104～113，148～167和221～224波段)，剩余的波段為188個，以 12、70、159波段作為RGB通道的偽彩色圖像如圖4.為了提高混合像元的分解精度，MCMDNMF算法的初始值由頂點(diǎn)成分分析(vertex component analysis，VCA)算法得到.

圖3 空間維數(shù)不同時算法解混性能比較Fig.3 Comparison of the algorithms with different spatial dimension

圖4 Cuprite區(qū)域的AVIRIS數(shù)據(jù)偽彩色圖Fig.4 AVIRIS data false-color image of cuprite

圖5為MCMDNMF算法對圖4進(jìn)行混合像元分解得到的9個端元豐度譜圖，解混得到的真實(shí)端元光譜的確定可以根據(jù)文獻(xiàn)[17-18]提供的真實(shí)地物情況進(jìn)行比對分析和篩選，本文參考了文獻(xiàn)[2]中采用的方法.其中，Kaolinite被分為2個不同端元，因?yàn)檫@種礦物在整個豐度分布區(qū)域上，不同部分成分有差別，導(dǎo)致了其端元光譜在不同區(qū)域發(fā)生變化.

為了定量評估MCMDNMF和其他幾種算法在真實(shí)高光譜圖像情況下的解混性能，使用SAD指標(biāo)進(jìn)行性能評價比較，結(jié)果如表1所示.MCMDNMF和MVCNMF算法能提取出9個端元光譜，但MVCNMF算法提取的端元譜精度低;CNMF有3種端元不能正確的提取出來，分別為Kaolinite、Alunite和Chalcedony端元;MDNMF算法有4種端元不能正確的提取出來，分別為 Kaolinite、Alunite、Chalcedony和Jarosite端元;ICA算法有3種端元不能提取，提取出來的端元的性能評價指標(biāo)很差.總體來說，MCMDNMF算法性能優(yōu)于其他幾種進(jìn)行對比的算法性能.

圖5 Cuprite區(qū)域的AVIRIS數(shù)據(jù)含有的礦物基于MCMDNMF算法的解混結(jié)果Fig.5 Unmixing results for AVIRIS data of cuprite region using MCMDNMF

表1 AVIRIS數(shù)據(jù)USGS譜庫的真實(shí)光譜和估計的端元光譜的SAD比較Table 1 SAD comparison for the real AVIRIS data in USGS and estimated endmember spectral

4 結(jié)束語

本文提出的改進(jìn)的非負(fù)矩陣分解算法，在標(biāo)準(zhǔn)非負(fù)矩陣分解算法的基礎(chǔ)上，增加了2個約束條件，即最小估計豐度協(xié)方差和單形體各頂點(diǎn)到中心點(diǎn)均方距離總和最小.這2個約束條件既符合高光譜圖像混合像元中含有的端元光譜的分布特性，又使非負(fù)矩陣分解可以獲取全局最優(yōu)解.使用合成的高光譜數(shù)據(jù)和真實(shí)的高光譜圖像進(jìn)行仿真實(shí)驗(yàn)，利用SAD和RMSE評價指標(biāo)進(jìn)行比較，MCMDNMF算法在性能和魯棒性方面優(yōu)于其他進(jìn)行比較的算法.

但是，MCMDNMF算法和其他基于非負(fù)矩陣分解的高光譜解混算法一樣，由于收斂過程是基于梯度下降方法，存在著收斂速度慢、計算時間長、算法的實(shí)時性不好等問題，這都需要在以后的研究中進(jìn)一步改進(jìn).

[1]張偉，杜培軍，張華鵬.基于神經(jīng)網(wǎng)絡(luò)的高光譜混合像元分解方法研究[J].測繪通報，2007(7):23-26.ZHANG Wei，DU Peijun，ZHANG Huapeng.Mixed pixel decomposition in hyperspectral remote sensing image based on neural network[J].Bulletin of Surveying and Mapping，2007(7):23-26.

[2]MIAO L，QI H.Endmember extraction from highly mixed data using minimum volume constrained nonnegative matrix factorization[J].IEEE Geoscience and Remote Sensing Letters，2007，45(3):765-777.

[3]ALEXIS H，MIREILLE G，JACQUES B T.Minimum dispersion constrained nonnegative matrix factorization to unmix hyperspectral data[J].IEEE Geoscience and Remote Sensing Letters，2008，48(6):2590-2602.

[4]JIA S，QIAN Y.Constrained nonnegative matrix factorization for hyperspectral unmixing[J].IEEE Geoscience and Remote Sensing Letters，2009，47(1):161-173.

[5]YU Yue，GUO Shan，SUN Weidong.Minimum distance constrained nonnegative matrix factorization for the endmember extraction of hyperspectral images[C]//Proceedings of Remote Sensing and GIS Data Processing and Applications.Wuhan，2007:6790151-6790159.

[6]吳波，趙銀娣，周小成.端元約束下的高光譜混合像元非負(fù)矩陣分解[J].計算機(jī)工程，2008，34(22):229-230.WU Bo，ZHAO Yindi，ZHOU Xiaocheng.Unmixing mixture pixels of hyperspectral imagery using endmember constrained nonnegative matrix factorization[J].Computer Engineering，2008，34(22):229-230.

[7]KESHAVA N.A survey of spectral unmixing algorithms[J].Lincoln Laboratory Journal，2003，14(1):55-78.

[8]LEE D D，SEUNG H S.Learning the parts of objects by of nonneative matrix factorizationn[J].Nature，1999，401(6755):788-791.

[9]JOHNSON R A，WIECHERN D W.Applied multivariate statistical analysis[M].6thed Englewood Cliffs:Prentice Hall Incorporated，2007:389-392.

[10]YANG He，DU Qian，SU Hongjun，et al.An efficient method for supervised hyperspectral band selection[J].IEEE Geoscience and Remote Sensing Letters，2011，8(1):138-142.

[11]吳波，張良培，李平湘.高光譜端元自動提取的迭代分解方法[J].遙感學(xué)報，2005，9(3):287-292.WU Bo，ZHANG Liangpei，LI Pingxiang.Automatic extraction of endmember from hyperspectral imagery by iterative unmixing[J].Journal of Remote Sensing，2005，9(3):287-292.

[12]ARNGREN M，SCHMIDT M N，LARSEN J.Unmixing of hyperspectral images using Bayesian nonnegative matrix factorization with volume prior[J].Journal of Signal Processing Systems，2011，65(3):479-496.

[13]HEINZ D C，CHANG C I.Fully constrained least squares linear spectral mixture analysis method for material quantification in Hyperspectral imagery[J].IEEE Geoscience and Remote Sensing Letters，2001，39(3):529-545.

[14]KESHAVA N，MUSTARD J F.Spectral unmixing[J].IEEE Signal Process Mag，2002，19(1):44-57.

[15]PLAZA A，MARTINEZ P，PEREZ R，et al.A quantitative and comparative analysis of endmember extraction algorithms from hyperspectral data[J].IEEE Geoscience and Remote Sensing Letters，2004，42(3):650-663.

[16]CLARK R N，SWAYZ A G，WISEE R，et al.USGS digital spectral library splib05a[M].[s.n.]:USGS Open File Report，2003:3-395.

[17]SWAYZE G.The hydrothermal and structural history of the cuprite mining district，southwestern Nevada:an integrated geological and geophysical approach[D].Boulder:University of Colorado，1997:399.

[18]CLARK R N，SWAYZE G A.Evolution in imaging spectroscopy analysis and sensor signal-to-noise:an e-amination of how far we have come[C]//Proc 6th Annu JPL Airborne Earth Sci Workshop.Colorado，1996:49-53.

[19]王立國，趙妍，王群明.基于POCS的高光譜圖像超分辨率方法[J].應(yīng)用科技，2010，37(10):26-30.WANG Liguo， ZHAO Yan， WANG Qunming. POCS based super-resolution method for hyperspectral imagery[J].Applied Science and Technology，2010，37(10):26-30.