基于譜線形狀與信息量差異的高光譜解混NMF初始化方法

袁德有， 袁 林

(1.南陽理工學院數學與統計學院,南陽 473004； 2.南陽理工學院經濟管理學院,南陽 473004)

基于譜線形狀與信息量差異的高光譜解混NMF初始化方法

袁德有1， 袁 林2

(1.南陽理工學院數學與統計學院,南陽 473004； 2.南陽理工學院經濟管理學院,南陽 473004)

在高光譜像元解混應用中，好的端元光譜矩陣初始化方法對于提高盲信號分解精度具有重要意義。針對空間分辨率較高的高光譜數據，提出了一種新的面向非負矩陣分解(non-negative matrix factorization，NMF)的初始化方法。該方法通過計算像元在譜線形狀和信息量差異等方面的參數，利用像元譜線峭度、KL散度和光譜角等參量，從眾多混合像元中識別出純像元； 并分辨出不同類型純像元(或類純像元)之間的差別，從中選擇最適合代表每一類型端元的純像元(或類純像元)作為算法的初值像元，完成端元矩陣的初始化。將此方法分別用于模擬數據和真實數據的實驗結果表明，該方法能夠明顯提高高光譜混合數據的NMF精度，相比其他常用初始化方法具有更好的效果。

初始化； 盲信號分解； 非負矩陣分解(NMF)； 譜線形狀； 信息量差異

0 引言

非負矩陣分解(non-negative matrix factorization，NMF)技術，由于其“非負性約束”和“局部構成整體”的概念符合人們的日常感知習慣，同時還具有物理意義明確、實現簡便、占用存儲空間少等優點，因此得到了廣泛應用[1]。近年來，研究人員開始將NMF引入到高光譜像元解混應用中[2]； 但在取得一定研究成果的同時，也帶來了不少問題，初始化就是其中之一。NMF是一個迭代過程，需要確定初始值； 好的初始值可以加快算法收斂速度，提高分解精度。但目前的NMF算法多采用取隨機值、主成分分析(principle component analysis, PCA)、奇異值分解[3](singular value decomposition，SVD)和模糊C均值聚類[4](fuzzy C-mean cluster, FCM)等方法進行初始化。隨機初始化最簡單，但效果往往也最差； 其他方法多是從文本分類[5]、語音處理等領域移植過來的，不能很好地符合高光譜混合數據的特點，實際應用效果不甚理想[6]。為此，本文主要針對空間分辨率相對較高、地物分布不太復雜的高光譜混合數據，在端元類型和混合特性未知的前提下，為解決利用NMF進行像元盲分解時的初始化問題，提出一種新的綜合像元譜線形狀特征及信息量差異特征等信息的初始化方法。該方法基于如下假設： 對于場景中每一類型的端元，高光譜數據中至少含有一個對應的純像元或類純像元(部分混合像元中，光譜和反射率特性等與占主導地位的某一種端元類似的像元，本文稱之為“類純像元”)。

1 初始化的目標與關鍵

用于高光譜像元解混的NMF算法，其初始化分為2部分： ①端元光譜矩陣W; ②豐度矩陣H。由于大多數真實高光譜數據的端元混合情況未知，因此H只需要采取隨機生成的方式，同時滿足非負性約束和全加性約束即可。端元提取是高光譜像元解混的關鍵步驟，也是地物類型識別和豐度反演的基礎和前提，目前尚缺少高效的W初始化方法，所以本文主要研究W的初始化方法。

W初始化的目標，即從原始高光譜數據的眾多像元中，選擇每類端元對應的一個純像元(或類純像元)，將其對應光譜值作為NMF算法初始值。難點在于，要選出最能代表每類端元的一個像元，并保證沒有遺漏和重復，否則將對結果造成不利影響。

選擇各類型端元對應的純像元(或類純像元)關鍵在于： ①從眾多混合像元中識別出純像元； ②分辨出不同類型純像元之間的差別，保證選擇結果在端元類別上沒有重復。這樣，每類端元都有且僅有一個純像元(或類純像元)作為初始值參與NMF運算，從而實現對真實W的最優近似。

2 初始化方法

首先需要確定端元數量，進而確定要選擇的波段數量。可利用PCA、最小噪聲分離(minimum noise fraction，MNF)和SVD等方法估計高光譜數據的端元數量[7-8]。本文選擇PCA變換確定波段數，設x={x1,x2,…,xn-1,xn}為一個像元所有n個波段光譜響應信號組成的向量，則X={x1,x2,…,xN-1,xN}就是由圖像中所有N個像元對應向量組成的矩陣。用PCA方法對X進行處理，得到一系列主成分分量，以及它們所對應的特征值，特征值的大小與各主成分分量的信息含量相對應。所以，觀察特征值的分布情況，確定大特征值的數目k，就可估計出該高光譜圖像覆蓋區域的端元數目。

k的確定可采用以下方法： 設定一個閾值(如99.73%)，假設前k個主成分分量對應的特征值累加占總特征值λ的百分比為p，即

(1)

若達到閾值要求，則k即為所求端元數目。確定端元數量后，主要基于像元的譜線形狀和光譜信息量差異2個方面進行初始化。首先研究像元的譜線形狀，包括像元本身譜線形狀的特點，以及像元之間譜線形狀的差異。為了利用像元本身的譜線形狀從眾多混合像元中識別出純像元，本文參考了獨立成分分析(independent component analysis，ICA)理論和中心極限定理?？梢酝浦?，純像元(端元)或類純像元的譜線形狀趨于非高斯分布，混合程度高的像元譜線趨于高斯分布。因此，可通過計算像元譜線的非高斯化程度來標識像元的純度。

峭度(kurtosis)是曲線非高斯性的自然度量指標[9]，其計算公式為

(2)

需要特別注意的是，K的閾值設置要適度，過小會使混合程度較高的像元無法剔除； 而過大則會將部分譜線形狀非高斯性相對較弱的純像元排除。多次實驗結果表明，K的閾值取K平均值的0.6倍以下時，基本不會發生誤排除純像元的情況。本文試驗中K的閾值取K平均值的0.5倍。利用K選出第一個初值像元，并留下絕大多數純像元(或類純像元)。

接下來分辨不同類型純像元之間的差別。本文利用KL散度(Kullback-Leibler divergence)和光譜角(spectral angle, SA)[7]2種參數，通過綜合像元間的信息量差異和譜線形狀差異，得出像元間光譜差異的定量指標； 然后選擇與已選初值像元之間光譜差異最大的像元，使得選出的初值像元是能夠分別代表所有種類端元的純像元(或類純像元)。

一般而言，同種端元對應的純像元(或類純像元)之間光譜差異很小，反之異種端元對應的純像元(或類純像元)之間光譜差異則較大[8]。通過計算待選像元與已選像元之間光譜信息的KL散度，選出最大KL散度對應的待選像元，即可選出與已選像元不屬于同一端元種類的純像元或純度很高的類純像元。

對于離散隨機變量，其概率分布P和Q的KL散度DKL(P‖Q)定義為

(3)

KL散度僅當概率P和Q各自總和均為1，且對于任何i，均滿足P(i)>0及Q(i)>0時才有意義。

NMF分解結果中還包含了各類端元的譜線形狀信息。因此，初始化時也應加入對像元間譜線形狀差異的度量。本文利用光譜角來度量像元間的譜線形狀差異。

2個像元t和r之間SA的計算公式為

(4)

式中i為波段序號。

通過計算KL散度和SA參數，并將二者作加權和(本文實驗中，根據多次實驗結果，權重系數比值定為0.6)，可更加全面地度量像元間的光譜特性差異。在初值選擇的每次迭代過程中，選擇具有最大加權和的像元，就能夠有效地保證新選像元與所有已選像元均分屬于異類端元。

3 初始化流程

W初始化方法流程如圖1所示。

圖1 W初始化流程Fig.1 Flowchart of initialization for W

初始化方法步驟如下： ①讀取數據，轉換數據維度； ②估算需要選擇的初值像元數量(PCA等)，并計算所有像元譜線的K，選出第一個初值像元； ③完成混合矩陣(H)的初始化； ④計算候選像元與已選像元集合之間KL散度和SA的加權和，確定新的初值像元； ⑤重復步驟④，直到選擇的初值像元數量滿足要求； ⑥根據選擇的所有初值像元，完成W的初始化； ⑦按照“非負”且“和為1”的條件約束初始化混合系數矩陣。

4 實驗與分析

4.1 模擬數據實驗

模擬數據為長和寬各10個像元、波段數為224的高光譜圖像，波長范圍約為0.4～2.5 μm，組成100行、224列的混合像元矩陣。該矩陣由端元矩陣和混合矩陣相乘得到，其中，端元矩陣由5種ENVI標準光譜庫中的端元數據(黑色涂料、瀝青、煤渣、松木和水泥)組成，共5行、224列； 混合矩陣為所有元素同時滿足“和為1”且“非負”2個約束條件的隨機數矩陣，共100行、5列。將混合像元矩陣的第15,25,35，45和75個像元重新賦值為純像元，分別對應5種端元。5種端元的光譜曲線圖如圖2。

圖2 5種端元初始光譜曲線Fig.2 Spectral curves of five end-members

對模擬數據的100個混合像元按1～100依次編號，則本文方法選中的初值像元序號(按照選擇的順序)依次為： 35，25，75，45和39，按照序號大小順序依次為： 25，35，39，45和75。其中，序號為25，35，45和75的4種像元都被準確選出，序號為15的像元則被誤選為序號為39的像元。究其原因，可能是由于直接將15號像元賦值為純像元時，較嚴重地偏離了“和為1”的約束。

利用本文方法得到初始像元光譜(圖2)，以該結果對NMF基本算法進行初始化，得到的NMF分解結果中端元波形如圖3所示。

圖3 NMF分解結果Fig.3 NMF decomposition results

從圖3可以看出，NMF分解結果中的端元波形與初始端元光譜十分相似。進一步應用本文初始化方法的NMF端元估計結果，以其與真實端元之間SA的均方根誤差(root mean square error，RMSE)為度量參數，進行精度分析和對比。首先，計算NMF基本算法(應用本文方法進行初始化)端元估計結果和真實端元之間的SA，得到矩陣

(5)

每行的最小值為5種真實端元光譜與對應的NMF端元光譜分解結果的SA，分別為0.033 18，0.192 33，0.018 62，0.016 31和0.004 22，RMSE為0.088 0。

為了更加直觀，圖4分別列出了5種真實端元光譜和NMF結果中對應波形估計的細節對比。

(a) 端元1 (b) 端元2 (c) 端元3

(d) 端元4(e) 端元5

圖4真實端元和對應NMF分解結果的光譜細節對比

Fig.4Detailedcomparisonbetweenspectralcurvesofrealend-membersandcorrespondingNMFdecompositionresults

從圖4可知，除圖4(b)中所示的端元2的波形估計誤差(與真實端元光譜間SA為0.192 33)偏大以外，其余4種端元估計的統計誤差和實際光譜形狀誤差均非常小，只有較明顯的幅度誤差，這4種端元都被很好地識別和估計出來。

表1列出了本文方法與其余幾種常用初始化方法的分解精度對比結果； 表2列出了在不考慮分解效果最差的一個端元(圖4(b)所示端元2)時，本文方法與原理較為類似的頂點成分分析(vertex component analysis，VCA)方法的精度比較。

表1 模擬數據實驗中本文方法與其他幾種常用初始化方法的分解精度對比Tab.1 Comparison between decomposition accuracies of method proposed in this paper and several other initialization methods in experiment of simulated data

表2 剔除一個精度最差的端元 (端元2)后本文方法與VCA的精度對比Tab.2 Comparison between accuracies of method proposed in this paper and VCA method without considering the worst estimated end-member (end-member 2)

從表1可以看出，在利用NMF算法進行盲分解時，本文方法相比其他大多數初始化方法，可以獲得更精確的分解結果； 但由于其中一個端元(端元2)的分解誤差較大，導致其整體精度略低于VCA。

從表2可以看出，本文方法比VCA的局部精度更高。具體而言，就是在不考慮誤差最大的一個端元(端元2)時，本文方法相對VCA精度更高。若假設誤差最大的端元為非感興趣信息，將其視為背景噪聲，只要求準確分解出數據中部分感興趣端元的光譜信息，則本文方法相對VCA將更具有優勢。

4.2 真實數據實驗

選取的真實實驗數據為美國圣地亞哥市的AVIRIS高光譜數據。該數據的光譜范圍為0.4～2.5 μm，光譜分辨率為10 nm，波段數為224個(剔除無效或噪聲較大波段以及水汽吸收波段后，實際有效波段數為189個)； 空間分辨率為3.5 m，圖像行、列數均為40。

該數據中含有感興趣地物類型主要包括4種地物，其余類型的小目標作為背景噪聲加以忽略。4種地物的光譜真值如圖5(a)所示，圖5(b)為應用本文方法進行初始化后的NMF分解結果。

(a) 4種地物光譜真值 (b) NMF分解結果

圖54種地物光譜真值與NMF分解結果對比

Fig.5ComparisonbetweenoriginalspectralandNMFdecompositionresultsforfourkindsofgroundobjects

表3列出本文方法與其余幾種常用初始化方法的分解精度對比結果。

表3 真實數據實驗中本文方法與其他幾種常用初始化方法的分解精度對比Tab.3 Comparison between decomposition accuracies of method proposed in this paper and several other initialization methods in experiment of real data

從表3可以看出，本文方法應用于真實數據時的精度略低于模擬數據，但相對其他常用方法仍具有更高的精度。與模擬數據實驗結果相反，本文方法應用于真實數據的精度略微高于VCA，與FCM精度也非常接近。首先，上述結果證明了本文方法在應用于真實數據時的可行性； 其次，也說明本方法的性能尚不夠穩定，處理的數據不同時，相對其他方法的精度也有所不同。

5 結論

本文針對缺少像元混合信息，同時含有純像元(或類純像元)的高光譜數據盲分解過程，提出一種新的非負矩陣分解(NMF)初始化方法。該方法利用像元譜線峭度、KL散度和光譜角等參量，從眾多混合像元中識別出純像元，并分辨出不同類型純像元之間的差別，從中選出最有資格代表各類型端元的純像元或類純像元作為初值像元，完成端元光譜矩陣的初始化。實驗結果說明，該方法能有效提高NMF算法的精度，比其他大多數常用初始化方法具有更好的效果。

需要強調的是，本文方法基于選擇已有像元作為初值，并要求高光譜數據中具有每種端元對應的純像元(或類純像元)，這一前提對于部分真實高光譜數據(例如空間分辨率很低、場景內地物分布情況和紋理復雜的高光譜數據)，是難以滿足或難以確定是否滿足的，這就限制了該方法的適用范圍。另外，該方法的性能尚不夠穩定； 在峭度閾值判斷中，權重系數的選擇主要依靠經驗和實驗結果來確定，尚缺乏具體流程和理論依據。因此，下一步的工作主要包括如何擴展本文方法的適用范圍及如何確定權重系數選擇的流程和理論依據。

[1] Berry M W,Browne M,Langville A N,et al.Algorithms and applications for approximate nonnegative matrix factorization[J].Computational Statistics and Data Analysis,2007,52(1):155-173.

[2] 李二森,張保明,楊 娜,等.非負矩陣分解在高光譜圖像解混中的應用探討[J].測繪通報,2011(3):7-10.

Li E S,Zhang B M,Yang N,et al.Discussion of the NMF’s application for hyperspectral imagery unmixing[J].Bulletin of Surveying and Mapping,2011(3):7-10.

[3] Boutsidis C,Gallopoulos E.SVD based initialization:A head start for nonnegative matrix factorization[J].Pattern Recognition,2008,41(4):1350-1362.

[4] Alshabrawy O S,Ghoneim M E,Awad W A,et al.Underdetermined blind source separation based on fuzzy C-means and semi-nonnegative matrix factorization[C]//Proceedings of 2012 Federated Conference on Computer Science and Information Systems.Wroclaw,Poland:IEEE,2012:695-700.

[5] 翟亞利,吳 翊.NMF初始化研究及其在文本分類中的應用[J].計算機工程,2008,34(16):191-193,197.

Zhai Y L,Wu Y.Study of non-negative matrix factorization initialization and its application to text classification[J].Computer Engineering,2008,34(16):191-193,197.

[6] Lee D D,Seung H S.Algorithms for non-negative matrix factorization[C]//Proceedings of the 13th International Conference on Neural Information Processing Systems.Denver,CO:ACM,2000:535-541.

[7] 余先川,安衛杰,呂中華,等.一種基于光譜角和光譜距離自動加權融合分類方法[J].地質學刊,2012,36(1):33-36.

Yu X C,An W J,Lyu Z H,et al.Automatic weighting fusion classification method based on spectral angle and spectral distance[J].Journal of Geology,2012,36(1):33-36.

[8] Wang R J,Zhan Y J,Zhou H F.RETRACTED:A method of underdetermined blind source separation with an unknown number of sources[J].Engineering Applications of Artificial Intelligence,2011:556-562.doi:10.1016/j.engappai.2011.06.003.

[9] 毋文峰,陳小虎,蘇勛家,等.基于峭度的ICA特征提取和齒輪泵故障診斷[J].機械科學與技術,2011,30(9):1583-1587.

Wu W F,Chen X H,Su X J,et al.ICA feature extraction and fault diagnosis based on Kurtosis for a gear pump[J].Mechanical Science and Technology for Aerospace Engineering,2011,30(9):1583-1587.

Aninitializationmethodofnon-negativematrixfactorizationforhyperspectraldataunmixingbasedonspectralshapeandinformationdissimilarity

YUAN Deyou1， YUAN Lin2

(1.SchoolofMathematicsdissimilarityandStatistics，NanyangInstituteofTechnology，Nanyang473004，China；2.SchoolofEconomicsandManagement，NanyangInstituteofTechnology，Nanyang473004，China)

When blind signal separation technique is applied to unmixing hyperspectral data, a good initialization is vital for improving separating precision. Aimed at the hyperspectral data with relatively high spatial resolution and simple surface features, the authors put forward a reasonable hypothesis that the data contain pure pixel or approximate pure pixel corresponding to the each type of end-members, and proposed a new initialization method of non-negative matrix factorization(NMF), which has great potential in pixel unmixing. By calculating parameters to quantify the spectral shape and information difference among candidate pixels, this method extracts pure pixels from mixed pixels, recognizes the information dissimilarity among different types of pure pixels and choose the existing pixels that are most suitable for representing each type of end-members as NMF’s initial values. The experimental results show that the method proposed in this paper can improve NMF’s decomposition accuracy of hyperspectral data significantly, and its performance is better than that of other NMF initialization methods.

initialization； blind signal separation； non-negative matrix factorization(NMF)； spectral shape； information dissimilarity

10.6046/gtzyyg.2017.04.17

袁德有,袁林.基于譜線形狀與信息量差異的高光譜解混NMF初始化方法[J].國土資源遙感,2017,29(4)：114-119.(Yuan D Y,Yuan L.An initialization method of non-negative matrix factorization for hyperspectral data unmixing based on spectral shape and information dissimilarity[J].Remote Sensing for Land and Resources,2017,29(4)：114-119.)

TP 751

A

1001-070X(2017)04-0114-06

2016-12-16；

2017-02-27

河南省高等學校重點科研項目“Smith正規型在有限域上有理點個數中的應用”(編號： 17A110010)資助。

袁德有(1960-),男,教授,主要從事小波分析等方面的研究。Email： yuandeyou1960@163.com。

(責任編輯：張仙)