基于超像素的流形正則化稀疏約束NMF混合像元分解算法

李登剛 陳香香 李華麗 王忠美

摘 要：針對(duì)傳統(tǒng)非負(fù)矩陣分解（NMF）法用于高光譜圖像混合像元分解時(shí)產(chǎn)生的分解結(jié)果精度不高、對(duì)噪聲敏感等問題，提出一種基于超像素的流形正則化稀疏約束NMF混合像元分解算法——MRS-NMF。首先，通過基于熵率的超像素分割來構(gòu)造高光譜圖像的流形結(jié)構(gòu)，把原圖像分割為k個(gè)超像素塊并把每個(gè)超像素塊中具有相似性質(zhì)的數(shù)據(jù)點(diǎn)標(biāo)上相同的標(biāo)簽，定義像素塊內(nèi)有相同標(biāo)簽的任意兩個(gè)數(shù)據(jù)點(diǎn)之間的權(quán)重矩陣，然后將權(quán)重矩陣應(yīng)用于NMF的目標(biāo)函數(shù)中以構(gòu)造出流形正則化約束項(xiàng);第二，在目標(biāo)函數(shù)中添加二次拋物線函數(shù)以完成稀疏約束;最后，采用乘法迭代更新法則求解目標(biāo)函數(shù)以得到端元矩陣和豐度矩陣的求解公式，同時(shí)設(shè)置最大迭代次數(shù)和容忍誤差閾值，迭代運(yùn)算得到最終結(jié)果。該方法有效利用了高光譜圖像的光譜和空間信息。實(shí)驗(yàn)結(jié)果表明，在模擬的高光譜數(shù)據(jù)中，與傳統(tǒng)的流形稀疏約束的非負(fù)矩陣分解（GLNMF）、L1/2-NMF和頂點(diǎn)成分分析全約束最小二乘法（VCA-FCLS）等方法相比，MRS-NMF可以提高0.016～0.063的端元分解精度和0.01～0.05的豐度分解精度;而在真實(shí)的高光譜圖像中，MRS-NMF較傳統(tǒng)的GLNMF、頂點(diǎn)成分分析法（VCA）、最小體積約束的非負(fù)矩陣分解（MVCNMF）等方法可以平均提高0.001～0.0437的端元分解精度。所提MRS-NMF算法有效地提高了混合像元分解的精度，同時(shí)具有較好的抗噪性能。

關(guān)鍵詞： 混合像元分解;非負(fù)矩陣分解;超像素分割;流形正則化;稀疏性

中圖分類號(hào)：TP751.1

文獻(xiàn)標(biāo)志碼：A

Abstract：For the problems such as poor unmixing results and sensitivity to noise of traditional Nonnegative Matrix Factorization （NMF） applied to hyperspectral unmixing， a Manifold Regularized Sparse NMF with superpixel （MRS-NMF） algorithm for hyperspectral unmixing was proposed. Firstly， the manifold structure of hyperspectral image was constructed by superpixel segmentation based on entropy. The original image was divided into k-superpixel blocks， and the data points in each superpixel block with same property were labeled the same label. Weight matrices were defined between any two data points with the similar label in a superpixel block， and then the weight matrices were applied to the objective function of NMF to construct the manifold regularization constraint. Secondly， a quadratic parabola function was added to the objective function to complete the sparse constraint. Finally， the multiplicative iterative update rule was used to solve the objective function to obtain the solution formulas of endmember matrix and abundance matrix. At the same time， maximum iteration times and tolerate error threshold were set to get the final results by iterative operation. The proposed method makes full use of spectral and spatial information of hyperspectral images. Experimental results show that on synthetic data the unmixing accuracies of endmember and abundance based on proposed MRS-NMF are 0.016-0.063 and 0.01-0.05 respectively higher than those based on traditional methods like Graph-regularized L1/2-Nonnegative Matrix Factorization （GLNMF）， L1/2NMF and Vertex Component Analysis-Fully Constrained Least Squares （VCA-FCLS）; while on real hyperspectral images， the average unmixing accurary of endmember based on proposed MRS-NMF is 0.001-0.0437 higher than that of traditional GLNMF， Vertex Component Analysis （VCA） and Minimum Volume Constrained Nonnegative Matrix Factorization （MVCNMF）. This proposed algorithm improves the accuracy of unmixing effectively with good robustness to noise.Key words：?hyperspectral unmixing;Nonnegative Matrix Factorization （NMF）; superpixel segmentation; manifold regularization; sparseness

0 引言

隨著衛(wèi)星和高光譜成像技術(shù)的發(fā)展，高光譜圖像擁有較高的光譜分辨率以及豐富的光譜波段。然而由于空間分辨率的限制，高光譜圖像中常常會(huì)出現(xiàn)一個(gè)像素包含多種地物。地物是指地面上各種有形物（如山川、河流、森林等）和無形物（如省、縣界等）的總稱，從而形成了混合像元?；旌舷裨拇嬖冢沟脗鹘y(tǒng)基于像素級(jí)的高光譜圖像分類等圖像處理方法失去作用。為了解決這一問題，高光譜圖像混合像元分解技術(shù)隨之而生，混合像元分解是將混合的像元分解為各類物質(zhì)的組成部分（端元）及其相應(yīng)的組成比例（豐度）。作為一種新的有效的高光譜數(shù)據(jù)表述方式，混合像元分解技術(shù)近年來得到了廣泛的關(guān)注[1-2]，目前已在農(nóng)業(yè)調(diào)查、礦物挖掘、環(huán)境監(jiān)測、軍事監(jiān)督等領(lǐng)域得到了廣泛的應(yīng)用。

近年來，學(xué)者們提出了多種混合像元分解方法，總的來講，這些算法可以分成兩類：基于幾何學(xué)的方法和基于統(tǒng)計(jì)學(xué)的方法?；趲缀螌W(xué)的方法將混合像元構(gòu)成的空間看成為一個(gè)單形體，而各端元?jiǎng)t位于單形體的頂點(diǎn)，這類方法需要圖像中存在純像元，常見的方法有：頂點(diǎn)成分分析法（Vertex Component Analysis， VCA）[3]、純像素指數(shù)法（Pure Pixel Index， PPI）[4]、N-FINDR[5]、自動(dòng)形態(tài)學(xué)端元提取法（Automated Morphological Endmember Extraction， AMEE）[6]。得到端元后，通常通過全約束的最小二乘法（Fully Constrained Least Squares， FCLS）[7]來求解豐度?；诮y(tǒng)計(jì)學(xué)的方法則將混合像元分解看成盲信號(hào)分離問題，因?yàn)椴恍枰瓐D像中存在純凈像元，并且可以同時(shí)得到端元和豐度，故該方法非常適用于端元高度混合的復(fù)雜場景中。常見的基于統(tǒng)計(jì)學(xué)的方法是獨(dú)立成分分析法（Independent Component Analysis， ICA）[8]和非負(fù)矩陣分解法（Nonnegative Matrix Factorization， NMF）[9-14]。

NMF作為一種基于統(tǒng)計(jì)學(xué)的方法，不需要原圖像存在純像元，并且能夠同時(shí)得到端元和豐度，目前廣泛應(yīng)用于混合像元分解領(lǐng)域。然而由于目標(biāo)函數(shù)的非凸性，在求解過程中，NMF不能求得全局最優(yōu)解。因此，根據(jù)高光譜圖像端元和豐度的相關(guān)屬性，在NMF的目標(biāo)函數(shù)中加入了不同的約束項(xiàng)。

如利用端元幾何特性作為目標(biāo)函數(shù)約束項(xiàng)的方法：

最小體積約束的非負(fù)矩陣分解（Minimum Volume Constrained Nonnegative Matrix Factorization， MVCNMF）[9]、端元相異性約束的非負(fù)矩陣分解（Endmember Dissimilarity Constrained Non-negative Matrix Factorization， EDCNMF）[10]等。以及利用高光譜圖像豐度分布稀疏性作為約束項(xiàng)的方法，如L1/2稀疏約束的非負(fù)矩陣分解（L1/2 Sparsity-Constrained Nonnegative Matrix Factorization， L1/2NMF）[11]、稀疏性和平滑性約束的非負(fù)矩陣分解（Sparse and Smooth Nonnegative Matrix Factorization， SSNMF）等[12]。稀疏性是指在實(shí)際的高光譜圖像中，各類端元的分布往往會(huì)集中在局部區(qū)域，而不會(huì)充滿整個(gè)圖像空間[3]，是高光譜圖像的一種本質(zhì)屬性。還有利用流形正則化作為約束項(xiàng)的方法，如流形稀疏約束的非負(fù)矩陣分解（Graph-regularized L1/2-NMF， GLNMF）[13]、圖正則化稀疏約束的非負(fù)矩陣分解（Graph-regularized NMF， GNMF）[14]等。

研究[15]發(fā)現(xiàn)，不能通過圖像數(shù)據(jù)來填充滿高維歐氏空間，因此這些數(shù)據(jù)間存在一個(gè)潛在的流形結(jié)構(gòu)。如果原數(shù)據(jù)空間中的兩個(gè)數(shù)據(jù)點(diǎn)很相近，那么在新的空間中，它們?nèi)匀粫?huì)保持相似，這就是流形假設(shè)[16]。而對(duì)于高光譜圖像，如果兩個(gè)像元在空間結(jié)構(gòu)中相近，那么在豐度圖像中，它們也將保持這種相近的關(guān)系。

傳統(tǒng)的非負(fù)矩陣分解方法僅考慮了豐度的稀疏性或者流形結(jié)構(gòu)信息，導(dǎo)致分解精度不高、對(duì)噪聲敏感等問題。為此，本文充分利用了高光譜圖像的相關(guān)特性，將流形正則化約束和稀疏約束同時(shí)加入到了NMF模型中，提出了基于超像素的流形正則化稀疏約束NMF混合像元分解算法——MRS-NMF（Manifold Regularized Sparse NMF）。

與傳統(tǒng)算法相比，該算法既有效地利用了光譜信息，同時(shí)又充分挖掘了空間信息。此外，使用超像素分割來建立流形正則項(xiàng)，可以很好地描述像元之間的相似性。合成的模擬數(shù)據(jù)實(shí)驗(yàn)和真實(shí)的高光譜圖像數(shù)據(jù)實(shí)驗(yàn)結(jié)果表明，與其他傳統(tǒng)算法相比，該算法可以獲得更高的分解精度和更好的抗噪性能。

1 相關(guān)理論背景

1.1 線性混合模型

1.2 非負(fù)矩陣分解（NMF）

非負(fù)矩陣分解是由Lee等[17]提出，其核心思想是把一個(gè)非負(fù)矩陣分解為另外兩個(gè)非負(fù)矩陣的乘積，給定一個(gè)非負(fù)矩陣X∈Rn×m，NMF旨在尋找另外兩個(gè)非負(fù)矩陣A∈Rn×r和S∈Rr×m（rmin（m，n）），使得：

正如式（1）與式（4）所示，如果忽略噪聲和誤差，并且在NMF模型中加入豐度和為1約束：1TPS=1，那么NMF非常適合應(yīng)用于混合像元分解問題。通過計(jì)算X與AS的歐氏距離，從而得到NMF的目標(biāo)函數(shù)如下：

通過乘法迭代更新法則求解式（5），可以得到端元矩陣和豐度矩陣的解。但是在式（5）中，端元A和豐度S的解非唯一，對(duì)于任意一個(gè)可逆的矩陣D，如果A和S為目標(biāo)函數(shù)的解，那么AD、D-1S也將是目標(biāo)函數(shù)的解。因此，為了獲得全局最優(yōu)解，需要根據(jù)實(shí)際的物理意義在目標(biāo)函數(shù)中加入約束項(xiàng)以使端元和豐度的解唯一。

2 本文算法

本文提出的MRS-NMF算法旨在提高高光譜圖像混合像元分解精度及抗噪性。為了完成這一目標(biāo)，流形正則化約束和豐度稀疏性約束被添加到了目標(biāo)函數(shù)中。通過超像素分割來表示高光譜圖像的流形結(jié)構(gòu)信息，同時(shí)通過二次拋物曲線來模擬豐度圖像的稀疏性。最后通過最優(yōu)化算法來求解目標(biāo)函數(shù)，從而得到最終結(jié)果。MRS-NMF算法的模型定義如下：

2.1 構(gòu)造流形正則化約束

為了充分利用高光譜圖像的光譜和空間結(jié)構(gòu)信息，本文算法通過基于熵率的超像素分割（Entropy Rate Superpixel， ERS）[18]來構(gòu)造流形結(jié)構(gòu)。首先把原圖像分割為k個(gè)超像素塊，那么每個(gè)超像素塊中的各個(gè)數(shù)據(jù)點(diǎn)將具有相似的性質(zhì)，給它們標(biāo)上相同的標(biāo)簽。對(duì)于具有相同標(biāo)簽的像素塊內(nèi)任意兩個(gè)數(shù)據(jù)點(diǎn)xi和xj，定義它們之間的權(quán)重矩陣為權(quán)重矩陣應(yīng)設(shè)為0。而當(dāng)兩個(gè)數(shù)據(jù)點(diǎn)之間的空間距離越近時(shí)，它們之間的權(quán)值越大，空間距離越遠(yuǎn)時(shí)，權(quán)值越小。通過上述分析，將權(quán)重矩陣應(yīng)用于NMF的目標(biāo)函數(shù)中，構(gòu)造出流形正則化約束項(xiàng)如下：

其中：N表示像元的數(shù)量。由式（7）可知，當(dāng)兩個(gè)數(shù)據(jù)點(diǎn)xi和xj之間的距離越小時(shí)，Wij相對(duì)較大，使約束項(xiàng)f1（S）最小，從而豐度si和sj越相似。

2.2 構(gòu)造豐度稀疏約束

正如Qian等[11]所述，使用L1/2范數(shù)作為豐度的稀疏約束項(xiàng)，存在著對(duì)噪聲魯棒性能差等問題。

本文采用了文獻(xiàn)[5]和文獻(xiàn)[19]中使用的方法，通過向目標(biāo)函數(shù)添加二次拋物線函數(shù)來完成稀疏約束。二次拋物線定義為：

其中：P為端元個(gè)數(shù);ξ用來控制稀疏度。通過式（8）對(duì)整個(gè)豐度圖像進(jìn)行累加操作，從而得到了稀疏約束項(xiàng)如下：

根據(jù)上述分析，本文MRS-NMF算法的目標(biāo)函數(shù)如下：

2.3 目標(biāo)函數(shù)求解

得到了目標(biāo)函數(shù)后，最后的問題為目標(biāo)函數(shù)的求解，這里采用乘法迭代更新法則來求解式（10）所示的目標(biāo)函數(shù)。得到端元矩陣和豐度矩陣的求解公式如下：從式（12）可以看出，端元和豐度矩陣的計(jì)算方法中只存在加法和乘法運(yùn)算，從而保證了結(jié)果的非負(fù)性約束（ANC）。此外，為了使豐度滿足和為1約束（ASC），在端元和豐度矩陣每次迭代更新前按照式（13）進(jìn)行擴(kuò)展。

其中：δ為控制和為1約束（ASC）的權(quán)重參數(shù)，本文參考文獻(xiàn)[9，11]將δ設(shè)為15，1TP和1TN分別為P個(gè)元素和N個(gè)元素全為1的行向量。在每次迭代計(jì)算前，先通過代替A，代替X，計(jì)算完成后再從和中取出A和X。再進(jìn)行式（13）所示的擴(kuò)展，進(jìn)行下一次迭代計(jì)算，達(dá)到迭代停止條件時(shí)，停止迭代求出結(jié)果。

接下來的問題是迭代停止條件，采用設(shè)置最大迭代次數(shù)和容忍誤差。本文參考文獻(xiàn)[9-12]，將最大迭代次數(shù)設(shè)置為3000，容忍誤差設(shè)置為10-6。即當(dāng)連續(xù)10次的迭代誤差都小于10-6或者計(jì)算達(dá)到最大迭代次數(shù)時(shí)，停止迭代得到最終結(jié)果。

本文MRS-NMF算法步驟如下。

算法1 MRS-NMF算法。

輸入 高光譜圖像X，端元個(gè)數(shù)P。

輸出 端元矩陣A和豐度矩陣S。

步驟1 設(shè)置參數(shù)λ、 μ和K。

步驟2 使用ERS超像素分割構(gòu)造權(quán)重矩陣W。

步驟3 初始化端元矩陣A和豐度矩陣S，對(duì)豐度矩陣S的列作歸一化處理。

步驟4 按照式（11）迭代計(jì)算矩陣A。

步驟5 按照式（13）擴(kuò)展矩陣A和X分別為和。

步驟6 按照式（12）迭代計(jì)算矩陣S。

步驟7 返回步驟4繼續(xù)計(jì)算，直到滿足迭代停止條件時(shí)為止。

步驟8 達(dá)到迭代停止條件時(shí)，停止迭代，得到最終結(jié)果。

3 實(shí)驗(yàn)與結(jié)果分析

3.1 評(píng)價(jià)指標(biāo)

本章通過合成的模擬數(shù)據(jù)實(shí)驗(yàn)和真實(shí)的高光譜圖像數(shù)據(jù)實(shí)驗(yàn)來驗(yàn)證本文MRS-NMF算法的有效性。本文采用兩個(gè)評(píng)價(jià)指標(biāo)來衡量算法的性能，分別是用來評(píng)價(jià)端元的光譜角距離（Spectral Angle Distance， SAD）和用來評(píng)價(jià)豐度的均方根誤差（Root Mean Squared Error， RMSE）。SAD和RMSE的定義如下：

其中：a為參考的端元向量;為估計(jì)的端元向量;s為參考的豐度向量;為估計(jì)的豐度向量;N為豐度圖像像素的數(shù)量。從定義式可以看出，SAD和RMSE越小，估計(jì)值和真實(shí)值越相近，分解的結(jié)果精確度越高。

3.2 實(shí)驗(yàn)條件

本文所有的數(shù)據(jù)計(jì)算均在硬件環(huán)境為Intel Core i3-2330，內(nèi)存4GB，主頻2.2GHz，軟件環(huán)境為Matlab R2014b的計(jì)算機(jī)上進(jìn)行。

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

為了測試本文算法的性能以及確定算法中各個(gè)參數(shù)的最優(yōu)值，從美國地質(zhì)勘測局（United States Geological Survey， USGS）的光譜庫中選取了4種礦物光譜曲線來合成模擬圖像數(shù)據(jù)。這些光譜包含了224個(gè)光譜波段，波長范圍在0.38～2.5μm。圖1為這4種物質(zhì)的光譜曲線，它們分別是黃鉀鐵礬、葉蛇紋石、斧石和光鹵石。參考MVCNMF[9]中的方法，通過下述步驟來合成模擬圖像：1）將尺寸大小為64×64的圖像平分為8×8的小塊;2）從光譜庫中隨機(jī)選取任意一種純像元來填充上述各個(gè)小塊，這樣每個(gè)小塊的像素均為純像元;3）使用一個(gè)7×7的低通濾波器來產(chǎn)生混合像元，這樣可以得到每個(gè)像元的豐度比例;4）為了更進(jìn)一步去除圖像中的純像元，并且控制豐度的稀疏程度，對(duì)于豐度大于0.8的像素，將會(huì)被上述4種端元均勻填充，這樣各種端元的豐度均為1/4。為了更進(jìn)一步貼近真實(shí)的高光譜圖像，充分考慮可能出現(xiàn)的誤差和噪聲，通過向圖像添加不同信噪比（Signal-to-Noise Ratio， SNR）的零均值高斯白噪聲來實(shí)現(xiàn)。SNR定義為：

在本實(shí)驗(yàn)中，需要確定三個(gè)參數(shù)值：流形正則化約束參數(shù)μ、稀疏約束參數(shù)λ和超像素的個(gè)數(shù)K。其中：參數(shù)λ用來控制豐度的稀疏度，參數(shù)μ用來控制結(jié)構(gòu)信息所占的權(quán)重，而參數(shù)K則由圖像本身特點(diǎn)決定。將模擬數(shù)據(jù)中的信噪比SNR設(shè)為20dB，端元矩陣A和豐度矩陣S選擇隨機(jī)初始化，并在每次實(shí)驗(yàn)中保持不變。

參考文獻(xiàn)[12-13]中的參數(shù)設(shè)置，首先通過固定參數(shù)μ=0.1，K=400來確定參數(shù)λ的最優(yōu)值，將λ以0.025的步長從0增加到0.4。按照算法1中的步驟，首先求解出端元矩陣和豐度矩陣，然后根據(jù)式（14）和（15）分別計(jì)算出不λ同下SAD和RMSE的值。圖2（a）表示所提出MRS-NMF在不同λ值下SAD和RMSE的變化情況。從圖2可看出：在0～0.2的區(qū)間范圍內(nèi)，SAD和RMSE都隨著λ的增加而減小;而在λ大于0.2后，SAD和RMSE則隨著λ的增加而增加。所以，這里選擇λ的最優(yōu)值為0.2。

然后固定λ=0.2，K=400，改變?chǔ)桃?.025的步長從0增加到0.4。同樣按照算法1中的步驟，先求解出端元矩陣和豐度矩陣，根據(jù)式（14）和（15）計(jì)算出不同μ下，SAD和RMSE的值。不同μ值下，所提出方法求得的SAD和RMSE如圖2（b）所示。由圖2（b）可知：當(dāng)μ位于0～0.1的區(qū)間范圍內(nèi)時(shí)，SAD和RMSE都隨著μ的增加而減小;當(dāng)μ大于0.1后，SAD隨著的μ增加而增加，而RMSE則保持相對(duì)平穩(wěn)。所以綜合考慮SAD和RMSE的變化情況，這里選擇μ的最優(yōu)值為0.1。

最后固定λ和μ的值分別為0.2和0.1，改變K的值在60～2000內(nèi)變化。采用同樣的方法，求解出不同K值下SAD和RMSE的值。圖2（c）、（d）為在不同K值下，SAD和RMSE的變化情況。從中可看出：隨著K的增加，SAD和RMSE均呈現(xiàn)出一個(gè)先降后升的過程，且SAD和RMSE都在K=400時(shí)取得最優(yōu)值，由于K與圖形本身特征有關(guān)，所以在此選擇K的最優(yōu)值為400。

通過上述實(shí)驗(yàn)，從而確定了算法中各個(gè)參數(shù)的最優(yōu)值分別為：λ=0.2， μ=0.1，K=400。

3.3.2 不同算法比較

將本文MRS-NMF算法同其他流行的混合像元分解方法就解混的性能進(jìn)行了比較，對(duì)比算法包括：頂點(diǎn)成分分析全約束最小二乘法L1/2NMF[11]、GLNMF[13]。在本實(shí)驗(yàn)中，設(shè)置信噪比SNR=30dB，基于NMF的方法端元和豐度的初始化方法均采用隨機(jī)初始化，并且在各方法的對(duì)比實(shí)驗(yàn)中都采用相同的初始值，實(shí)驗(yàn)中的其他參數(shù)與上文一致。首先按照算法1中的步驟求解出端元矩陣和豐度矩陣，然后根據(jù)式（14）和（15）計(jì)算出SAD和RMSE。圖3為上述模擬數(shù)據(jù)重復(fù)20次實(shí)驗(yàn)時(shí)，不同方法所求得的SAD和RMSE的均值和方差。從圖3可看出：無論是SAD還是RMSE，VCA-FCLS算法分解性能最差，這是因?yàn)閂CA算法在純像元存在時(shí)，才能提取出端元，而本實(shí)驗(yàn)則不存在純像元。另一方面，VCA-FCLS算法的標(biāo)準(zhǔn)差比基于NMF的算法都要小，穩(wěn)定性相對(duì)較好，這是因?yàn)榛贜MF的算法本身不存在全部最優(yōu)解，且算法對(duì)初始值依賴性較大。最后比較所有方法，發(fā)現(xiàn)本文MRS-NMF算法分解精度最好，說明新方法有效地利用了高光譜圖像的空間信息和光譜信息以及豐度的稀疏性從而提高了分解效果。最后在所有基于NMF的方法中，本文MRS-NMF算法的標(biāo)準(zhǔn)差最小，這也說明本文算法比其他算法擁有較好的穩(wěn)定性。

3.3.3 抗噪性能比較

本部分通過向模擬數(shù)據(jù)添加不同程度的噪聲，以測試各方法的抗噪能力。實(shí)驗(yàn)中的其他參數(shù)與上文一致，僅改變所添加噪聲的信噪比SNR，以10dB為步長，將SNR從10dB增加到40dB以及考慮SNR為無窮大（無噪聲）的情況。采用同樣的方法，按照算法1中的步驟以及式（14）和（15）計(jì)算出SAD和RMSE。圖4為不同信噪比SNR下各算法的分解結(jié)果，其中基于NMF的方法，采用了隨機(jī)初始化方法，并且在各對(duì)比實(shí)驗(yàn)中初始值保持一致，所有的算法均取運(yùn)行20次后的平均值。從圖4可看出：VCA-FCLS算法的分解效果最差，同樣是因?yàn)槟M圖像中不存在純像元，而基于NMF的方法則不需要純像元的存在，效果較VCA-FCLS均要好。同時(shí)隨著信噪比的增加，各算法的解混性能均趨于好轉(zhuǎn)，且本文MRS-NMF算法比其他幾種方法要好，這也說明本文算法結(jié)合空間和光譜信息以及豐度的稀疏性后具有較好的抗噪性能。

采集的美國內(nèi)華達(dá)州Cuprite區(qū)域圖像來完成這一實(shí)驗(yàn)。該區(qū)域各類礦物分布錯(cuò)綜復(fù)雜，混合現(xiàn)象十分嚴(yán)重，被廣泛應(yīng)用于混合像元分解領(lǐng)域[9-14，19]。原始圖像包含了224個(gè)光譜波段，波長范圍為0.37～2.48μm，光譜分辨率為10nm，去除那些低信噪比以及被水蒸氣吸收的波段，在本實(shí)驗(yàn)中僅采用188個(gè)波段做后續(xù)的處理。同時(shí)，為了減輕計(jì)算負(fù)擔(dān)，提高運(yùn)算效率，本實(shí)驗(yàn)僅計(jì)算大小為250×191的子圖像區(qū)域。圖5為該子圖像的第80個(gè)波段。

根據(jù)虛擬維度算法[20]以及目前現(xiàn)有知識(shí)對(duì)該地區(qū)礦物的分析情況[9，19]，選擇端元個(gè)數(shù)為9。同時(shí)，在本實(shí)驗(yàn)中采用VCA-FCLS方法作為端元矩陣和豐度矩陣的初始化方法，以期提高計(jì)算速度。最大迭代次數(shù)設(shè)置為3000作為迭代停止條件，按照算法1中的步驟求解出端元矩陣和豐度矩陣。圖6列出了本文MRS-NMF算法所求解出來的端元曲線與光譜庫中參考端元曲線的對(duì)比情況以及分解后的豐度圖像。從圖6中可以看出，估計(jì)的端元曲線和參考曲線非常相近，另外和已經(jīng)公布的豐度參考圖像（http：//speclab.cr.usgs.gov/PAPERS/tetracorder/FIGURES/fig9b.cuprite95.tgif.2.2um_map.gif）相比，本文算法在豐度圖像上也獲得了較好的效果。求解的端元和參考端元之間的定量比較如表1所示。

由于MVCNMF只考慮了端元的幾何信息，而忽略了流形結(jié)構(gòu)信息，故也加入本次實(shí)驗(yàn)進(jìn)行對(duì)比。

從表1也可看出，本文MRS-NMF算法對(duì)大多數(shù)礦物都有較好的分解性能，同時(shí)誤差的平均值最小，精度最高。這也說明本文算法結(jié)合光譜和空間信息以及豐度稀疏性后能夠提高混合像元分解的性能，同時(shí)也說明了本文算法對(duì)真實(shí)的高光譜圖像也具有較好的適應(yīng)性。

4 結(jié)語

本文提出了一種基于超像素的流形正則化稀疏約束NMF混合像元分解算法——MRS-NMF。該算法基于NMF模型，不需要原圖像存在純凈像元，并且可以同時(shí)得到端元和豐度。MRS-NMF算法利用超像素分割來獲取高光譜圖像的流形結(jié)構(gòu)信息，同時(shí)在目標(biāo)函數(shù)中有效結(jié)合了流形結(jié)構(gòu)信息和豐度圖像稀疏特性。最后合成的模擬圖像和真實(shí)的高光譜圖像實(shí)驗(yàn)結(jié)果表明，本文算法較傳統(tǒng)混合像元分解具有更好的分解效果和抗噪性能。在下一步的研究工作中，將根據(jù)圖像的本身特性，設(shè)置算法自動(dòng)確定超像素的個(gè)數(shù)，從而提升混合像元分解的效率。

參考文獻(xiàn)（References）

[1] 徐晨光， 鄧承志， 朱華生. 近似稀疏約束的多層非負(fù)矩陣分解高光譜解混[J]. 紅外與激光工程， 2018， 47（11）： 257-265. （XU C G， DENG C Z， ZHU H S. Approximate sparse regularized multilayer NMF for hyperspectral unmixing[J]. Infrared and Laser Engineering， 2018， 47（11）： 257-265.）

[2] 王瀛， 何欣， 左方. 基于最大整體包容度約束非負(fù)矩陣分解的高光譜遙感圖像混合像元分析算法[J]. 光子學(xué)報(bào)， 2018， 47（3）： 136-144. （WANG Y， HE X， ZUO F. Mixed data analysis algorithm based on maximum overall coverage constraint nonnegative matrix factorization for hyperspectral image[J]. Acta Photonica Sinica， 2018， 47（3）： 136-144.）

[3] NASCIMENTO J M P， DIAS J M B. Vertex component analysis： a fast algorithm to unmix hyperspectral data[J]. IEEE Transactions on Geoscience and Remote Sensing， 2005， 43（4）： 898-910.

[4] BOARDMAN J W， KRUSE F A， GREEN R O. Mapping target signatures via partial unmixing of AVIRIS data [EB/OL]. [2019-01-10]. http：//www.hgimaging.com/PDF/boardman95.PDF.

[5] WINTER M E. N-FINDR： an algorithm for fast autonomous spectral end-member determination in hyperspectral data[J]. Proceedings of the SPIE， 1999， 3753： 266-275.

[6] PLAZA A， MARTINEZ P， PEREZ R， et al. Spatial/spectral endmember extraction by multidimensional morphological operations[J]. IEEE Transactions on Geoscience and Remote Sensing， 2002， 40（9）： 2025-2041.

[7] HEINZ D C， CHANG C. Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery[J]. IEEE Transactions on Geoscience and Remote Sensing， 2001， 39（3）： 529-545.

[8] WANG J， CHANG C. Applications of independent component analysis in endmember extraction and abundance quantification for hyperspectral imagery[J]. IEEE Transactions on Geoscience and Remote Sensing， 2006， 44（9）： 2601-2616.

[9] MIAO L， QI H. Endmember extraction from highly mixed data using minimum volume constrained nonnegative matrix factorization[J]. IEEE Transactions on Geoscience and Remote Sensing， 2007， 45（3）： 765-777.

[10] WANG N， DU B， ZHANG L. An endmember dissimilarity constrained non-negative matrix factorization method for hyperspectral unmixing[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing， 2013， 6（2）： 554-569.

[11] QIAN Y， JIA S， ZHOU J， et al. Hyperspectral unmixing via L1/2 sparsity constrained nonnegative matrix factorization[J]. IEEE Transactions on Geoscience and Remote Sensing， 2011， 49（11）： 4282-4297.

[12] WU C， SHEN C. Spectral unmixing using sparse and smooth nonnegative matrix factorization[C]// Proceedings of 21st International Conference on Geoinformatics. Piscataway： IEEE， 2013： 1-5.

[13] LU X， WU H. Manifold regularized sparse NMF for hyperspectral unmixing[J]. IEEE Transactions on Geoscience and Remote Sensing， 2013， 51（5）： 2815-2826.

[14] RAJABI R， GHASSEMIAN H. Sparsity constrained graph regularized NMF for spectral unmixing of hyperspectral data[J]. Journal of the Indian Society of Remote Sensing， 2015， 43（2）： 269-278.

[15] CAI D， HE X， HAN J， et al. Graph regularized nonnegative matrix factorization for data representation[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence， 2011， 33（8）： 1548-1560.

[16] BELKIN M， NIYOGI P. Laplacian eigenmaps and spectral techniques for embedding and clustering[C]// Proceedings of the 14th International Conference on Neural Information Processing Systems. Cambridge， MA： MIT Press， 2001： 585-591.

[17] LEE D D， SEUNG S H. Learning the parts of objects by non-negative matrix factorization[J]. Nature， 1999， 401（6755）： 788-791.

[18] LIU M Y， TUZEL O， RAMALINGAM S， et al. Entropy rate superpixel segmentation[C]// Proceedings of the 2011 IEEE Conference on Computer Vision and Pattern Recognition. Piscataway： IEEE， 2011： 2097-2104.

[19] LI D G， LI S T， LI H L. Hyperspectral image unmixing based on sparse and minimum volume constrained nonnegative matrix factorization[C]// Proceedings of the 6th Chinese Conference on Pattern Recognition. Berlin： Springer， 2014： 44-52.

[20] CHANG C， DU Q. Estimation of number of spectrally distinct signal sources in hyperspectral imagery[J]. IEEE Transactions on Geoscience and Remote Sensing， 2004， 42（3）： 608-619.