袁德有， 袁 林

(1.南陽理工學(xué)院數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,南陽 473004； 2.南陽理工學(xué)院經(jīng)濟(jì)管理學(xué)院,南陽 473004)

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

袁德有1， 袁 林2

(1.南陽理工學(xué)院數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,南陽 473004； 2.南陽理工學(xué)院經(jīng)濟(jì)管理學(xué)院,南陽 473004)

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

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

0 引言

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

1 初始化的目標(biāo)與關(guān)鍵

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

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

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

2 初始化方法

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

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

(1)

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

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

(2)

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

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

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

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

(3)

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

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

2個(gè)像元t和r之間SA的計(jì)算公式為

(4)

式中i為波段序號。

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

3 初始化流程

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

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

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

4 實(shí)驗(yàn)與分析

4.1 模擬數(shù)據(jù)實(shí)驗(yàn)

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

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

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

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

圖3 NMF分解結(jié)果Fig.3 NMF decomposition results

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

(5)

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

為了更加直觀，圖4分別列出了5種真實(shí)端元光譜和NMF結(jié)果中對應(yīng)波形估計(jì)的細(xì)節(jié)對比。

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

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

圖4真實(shí)端元和對應(yīng)NMF分解結(jié)果的光譜細(xì)節(jié)對比

Fig.4Detailedcomparisonbetweenspectralcurvesofrealend-membersandcorrespondingNMFdecompositionresults

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

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

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

表2 剔除一個(gè)精度最差的端元 (端元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算法進(jìn)行盲分解時(shí)，本文方法相比其他大多數(shù)初始化方法，可以獲得更精確的分解結(jié)果； 但由于其中一個(gè)端元(端元2)的分解誤差較大，導(dǎo)致其整體精度略低于VCA。

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

4.2 真實(shí)數(shù)據(jù)實(shí)驗(yàn)

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

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

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

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

Fig.5ComparisonbetweenoriginalspectralandNMFdecompositionresultsforfourkindsofgroundobjects

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

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

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

5 結(jié)論

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

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

[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] 李二森,張保明,楊 娜,等.非負(fù)矩陣分解在高光譜圖像解混中的應(yīng)用探討[J].測繪通報(bào),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初始化研究及其在文本分類中的應(yīng)用[J].計(jì)算機(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] 余先川,安衛(wèi)杰,呂中華,等.一種基于光譜角和光譜距離自動(dòng)加權(quán)融合分類方法[J].地質(zhì)學(xué)刊,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].機(jī)械科學(xué)與技術(shù),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

河南省高等學(xué)校重點(diǎn)科研項(xiàng)目“Smith正規(guī)型在有限域上有理點(diǎn)個(gè)數(shù)中的應(yīng)用”(編號： 17A110010)資助。

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

(責(zé)任編輯：張仙)