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

趙春暉，成寶芝，楊偉超

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

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

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

1 線性光譜混合模型

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

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

2 約束非負矩陣分解的解混方法

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

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

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

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

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

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

因此，以J1(M)和J2(M)為約束項的非負矩陣分解的目標函數為

若將式(7)作為改進的MCMDNMF的目標函數，沒有考慮高光譜豐度總和為1這一約束條件，而這個條件作為軟約束可以使解混的結果更精確.本文采用文獻[13]中所使用的表示形式，即式中:δ為權值，控制總和為1這一約束條件對目標函數的影響，1T表示一個向量的所有分量都為1，即當迭代求解α時，可用代替U代替M.

MCMDNMF采用迭代算法計算W和H的最優解:

依據式(9)，需要對J1(M)和J2(M)進行求導，可得

因此，目標函數的梯度為

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

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

3 仿真實驗與分析

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

3.1 性能評價指標

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

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

3.2 合成高光譜數據

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

圖1 合成高光譜數據的端元光譜圖Fig.1 Endmember spectral of synthetic data

3.2.1 算法的抗噪聲能力

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

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

3.2.2 空間維數不同時算法的解混能力

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

通過比較可以發現，從整體上來說，MCMDNMF算法和MVCNMF算法性能相近，空間維數越高，MCMDNMF性能越優;MCMDNMF算法性能明顯要優于CNMF和MDNMF算法.

3.3 真實的高光譜圖像

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

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

圖4 Cuprite區域的AVIRIS數據偽彩色圖Fig.4 AVIRIS data false-color image of cuprite

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

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

圖5 Cuprite區域的AVIRIS數據含有的礦物基于MCMDNMF算法的解混結果Fig.5 Unmixing results for AVIRIS data of cuprite region using MCMDNMF

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

4 結束語

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

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

[1]張偉，杜培軍，張華鵬.基于神經網絡的高光譜混合像元分解方法研究[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]吳波，趙銀娣，周小成.端元約束下的高光譜混合像元非負矩陣分解[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].遙感學報，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].應用科技，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.