袁 博

(南陽(yáng)理工學(xué)院 計(jì)算機(jī)與信息工程學(xué)院，河南 南陽(yáng) 473004)

馬爾可夫隨機(jī)場(chǎng)的空間相關(guān)模型在非負(fù)矩陣分解線性解混中的應(yīng)用

袁 博*

(南陽(yáng)理工學(xué)院 計(jì)算機(jī)與信息工程學(xué)院，河南 南陽(yáng) 473004)

針對(duì)基于非負(fù)矩陣分解(NMF)的高光譜解混存在的初始化與“局部極小”等問(wèn)題，提出一種基于馬爾可夫隨機(jī)場(chǎng)(MRF)的空間相關(guān)約束NMF線性解混算法(MRF-NMF)。首先，通過(guò)基于最小誤差的高光譜信號(hào)識(shí)別(HySime)法估算端元數(shù)量，同時(shí)利用頂點(diǎn)成分分析(VCA)和全約束最小二乘法(FCLS)初始化端元矩陣與豐度矩陣；其次，利用MRF模型建立描述地物空間分布規(guī)律的能量函數(shù)，以此描述地物分布的空間相關(guān)特征；最后，將基于MRF的空間相關(guān)約束函數(shù)與NMF標(biāo)準(zhǔn)目標(biāo)函數(shù)以交替迭代的形式參與解混，得出高光譜數(shù)據(jù)的端元信息與豐度分解結(jié)果。理論分析和真實(shí)數(shù)據(jù)實(shí)驗(yàn)結(jié)果表明，在高光譜數(shù)據(jù)空間相關(guān)程度較低的情況下，相比最小體積約束的NMF(MVC-NMF)、分段平滑和稀疏約束的NMF(PSNMFSC)和交互投影子梯度非負(fù)矩陣分解(APS-NMF)三種參考算法，所提算法的端元分解精度仍分別提高了7.82%、12.4%和10.1%，其豐度分解精度仍分別提高了8.34%、12.6%和9.87%。MRF-NMF能夠彌補(bǔ)NMF對(duì)于空間相關(guān)特征描述能力的不足，減小解混結(jié)果中地物的空間能量分布誤差。

非負(fù)矩陣分解；高光譜線性解混；空間相關(guān)；馬爾可夫隨機(jī)場(chǎng); 交替迭代；空間能量

0 引言

“混合像元”是影響高光譜定量遙感應(yīng)用精度的關(guān)鍵問(wèn)題之一。隨著高光譜數(shù)據(jù)應(yīng)用領(lǐng)域的不斷拓展及應(yīng)用精度要求的不斷提高，各種解混方法層出不窮。其中，文獻(xiàn)[1]提出的非負(fù)矩陣分解(Non-negative Matrix Factorization, NMF)技術(shù)由于適合處理高維數(shù)據(jù)，且排除了無(wú)意義的負(fù)值元素，形式上類似于線性光譜混合模型，使得基于NMF的解混算法成為高光譜線性解混研究中的熱點(diǎn)。

由于存在初始化、“非凸”導(dǎo)致“局部極小”等問(wèn)題，標(biāo)準(zhǔn)NMF的解混效果并不理想，大量用于高光譜解混的NMF擴(kuò)展方法不斷涌現(xiàn)。這些方法的共同特點(diǎn)在于：通過(guò)在標(biāo)準(zhǔn)NMF目標(biāo)函數(shù)內(nèi)部加入描述某種高光譜圖像典型特征的約束項(xiàng)，達(dá)到縮小解空間、加快迭代并提高解混精度等目的。例如，文獻(xiàn)[2]提出了一種基于NMF的光譜解混算法，向標(biāo)準(zhǔn)NMF算法中加入了端元光譜的平滑性和豐度的稀疏性約束；文獻(xiàn)[3]提出了最小體積約束的NMF(Minimum Volume Constrained NMF, MVC-NMF)算法；文獻(xiàn)[4]提出了交互投影子梯度非負(fù)矩陣分解(NMF with Alternating Projected Subgradients, APS-NMF)；文獻(xiàn)[5]也提出了分段平滑和稀疏約束的NMF(Piecewise Smoothness NMF with Sparseness Constraints, PSNMFSC)，但需要預(yù)先指定稀疏度。

上述NMF擴(kuò)展方法的共同問(wèn)題在于：首先，僅將高光譜數(shù)據(jù)視為各類光譜信息或統(tǒng)計(jì)信息的集合，忽略了空間屬性，特別是各地物(端元)類型之間顯著的空間相關(guān)特征；其次，擴(kuò)展后的目標(biāo)函數(shù)往往包含多個(gè)內(nèi)部函數(shù)項(xiàng)，容易產(chǎn)生相互干擾。上述兩點(diǎn)在不同程度上制約了高光譜解混精度的進(jìn)一步提高。

近年來(lái)，基于空間相關(guān)特征的高光譜解混研究也取得了大量學(xué)術(shù)成果。例如，文獻(xiàn)[6]提出了一種基于圖正則化的半監(jiān)督NMF算法(Graph regularized-based Semi-supervised NMF, GSNMF)，克服了NMF忽略樣本數(shù)據(jù)局部幾何結(jié)構(gòu)的缺陷；文獻(xiàn)[7]利用訓(xùn)練數(shù)據(jù)構(gòu)造結(jié)構(gòu)化字典，添加空間相關(guān)性約束項(xiàng)和訓(xùn)練數(shù)據(jù)的空間信息，提出了一種新的高光譜解混方法；文獻(xiàn)[8]利用先驗(yàn)概率密度函數(shù)表達(dá)兩個(gè)相鄰區(qū)域的空間相關(guān)程度，提出了一種區(qū)域相關(guān)的NMF解混算法；文獻(xiàn)[9]提出了一種基于稀疏約束和圖正則化的半監(jiān)督NMF方法，在進(jìn)行低維非負(fù)分解時(shí)，能夠保持?jǐn)?shù)據(jù)幾何結(jié)構(gòu)。

上述算法需要大量先驗(yàn)知識(shí)，對(duì)于部分先驗(yàn)知識(shí)不足或難以實(shí)地獲取的真實(shí)高光譜數(shù)據(jù)來(lái)說(shuō)，算法精度與應(yīng)用效果受到一定制約。針對(duì)上述問(wèn)題，本文提出一種基于馬爾可夫隨機(jī)場(chǎng) (Markov Random Field, MRF)模型的NMF線性解混算法(NMF linear unmixing algorithm based on MRF, MRF-NMF)。MRF-NMF可歸類為非監(jiān)督算法，對(duì)先驗(yàn)知識(shí)依賴程度低，利用NMF實(shí)現(xiàn)高光譜數(shù)據(jù)的線性解混并保證基本的解混精度，通過(guò)引入MRF模型描述高光譜圖像的地物空間相關(guān)特征，修正NMF標(biāo)準(zhǔn)目標(biāo)函數(shù)產(chǎn)生的解混誤差，進(jìn)一步改善解混精度。

1 基于MRF模型的NMF解混

由于具有二維空間圖像屬性，高光譜遙感圖像在包含大量光譜信息和統(tǒng)計(jì)信息的同時(shí)，也含有豐富的空間信息。例如，自然地物的空間分布具有成片、連續(xù)的特點(diǎn)，導(dǎo)致高光譜圖像中鄰近像元取值相近，且端元的空間分布具有連續(xù)性和相關(guān)性。

高光譜圖像的空間能量大小由地物的空間變化頻率和幅度決定，頻率越低、幅度越小，則空間能量越小。由于端元的空間相關(guān)性，決定了局部鄰域內(nèi)端元變化地頻率低、幅度小，即理想的NMF分解結(jié)果中，端元的“空間能量”應(yīng)盡可能小。因此，一般情況下度量空間能量大小的能量函數(shù)取值越小，越接近地面真實(shí)情況。

文獻(xiàn)[10]指出，如果NMF解混方法沒(méi)有考慮鄰近地物間的相似關(guān)系，忽略了端元分布的空間相關(guān)性，解混結(jié)果中端元分布受到其他類型端元或噪聲的干擾程度就高，容易導(dǎo)致端元“空間能量”偏大，亦即端元在各混合像元中的分布比例失真程度偏大，影響解混精度。

為解決上述問(wèn)題，需要構(gòu)建一種描述端元分布空間相關(guān)性的模型，使解混過(guò)程中端元分布“空間能量”不斷向理想狀態(tài)(最小)靠近，提高結(jié)果中端元空間分布與地面真實(shí)情況的符合程度。

1.1 MRF模型

MRF模型是根據(jù)像元局部空間相關(guān)性建立其聯(lián)合概率分布的一種有力工具，其最大優(yōu)點(diǎn)在于既能利用像元本身的特征信息，又能利用相鄰像元之間的相關(guān)信息。設(shè)x={x1,x2,…,xB}是隨機(jī)變量簇X={X1,X2,…,XB}的一個(gè)實(shí)現(xiàn)，所有可能的實(shí)現(xiàn)組成一個(gè)空間T，N={Nt,t∈T}為定義在T上的鄰域系統(tǒng)，若隨機(jī)場(chǎng)X滿足條件(1)，則稱X是以N為鄰域系統(tǒng)的MRF：

(1)

MRF模型能夠提供關(guān)于高光譜圖像的一種統(tǒng)計(jì)描述，很好地體現(xiàn)出每個(gè)像元關(guān)于其鄰近像元的條件分布(亦即端元分布的空間相關(guān)程度)，而該條件分布即表現(xiàn)為MRF的局部特性。當(dāng)鄰域系統(tǒng)階數(shù)足夠大，任何定義在T上的高光譜圖像均可看成是MRF的一個(gè)實(shí)現(xiàn)。

1.2 描述地物空間相關(guān)特征的MRF模型

由于MRF的定義沒(méi)有明確給出聯(lián)合概率的具體形式，使用起來(lái)很不方便。根據(jù)文獻(xiàn)[11], Hammersley-Clifford定理建立了MRF和Gibbs隨機(jī)場(chǎng)之間的等效性，將MRF模型的最優(yōu)解問(wèn)題轉(zhuǎn)化為求解Gibbs隨機(jī)場(chǎng)能量函數(shù)的最小值問(wèn)題。

不同的能量函數(shù)形式可構(gòu)造出不同的MRF模型，較常見的MRF模型有自生模型(Auto Model)和多級(jí)邏輯模型(Multi-Level Logistic, MLL)等。本文采用Deng等[12]提出的MRF模型來(lái)描述端元空間相關(guān)性，該模型將能量函數(shù)分為兩部分，如式(2)所示：

E=ER+EF

(2)

式中：EF是表示圖像本身特征的能量函數(shù)，只和觀測(cè)值分布有關(guān)，如式(3)所示；ER是表示像元與鄰域相關(guān)性的能量函數(shù)。

(3)

(4)

式中：P為端元數(shù)；當(dāng)x=xi+r時(shí)，δ(x,xi+r)=-1，否則δ(x,xi+r)=0；c表示鄰域模型的階數(shù)。

(5)

將式(5)代入式(2)～(4)，設(shè)W為全部像元對(duì)應(yīng)特征向量的均值矩陣，可推出以U為自變量的MRF模型能量函數(shù)E(U)，如式(6)所示：

E(U)=ER(U)+EF(U)=

(6)

如前文所述，由于能量函數(shù)E(U)取最小值時(shí)， MRF模型描述的“空間能量”最小，空間分布情況最理想。因此，對(duì)式(6)關(guān)于U求導(dǎo)，經(jīng)過(guò)化簡(jiǎn)得到最終結(jié)果如式(7)所示：

U=(RT·R)-1·(W·R)T

(7)

1.3 MRF-NMF模型

MRF-NMF通過(guò)NMF和MRF交替運(yùn)行的方式計(jì)算分離矩陣U。其中，NMF步驟選擇基于歐氏距離的標(biāo)準(zhǔn)目標(biāo)函數(shù)(8)作為保證解混精度的基礎(chǔ)；迭代規(guī)則選擇最大似然法，得到的端元光譜矩陣和豐度矩陣的乘性迭代公式分別為式(9)和式(10)。

(8)

(9)

(10)

MRF步驟中，利用1.2節(jié)中的能量函數(shù)模型及相關(guān)公式，對(duì)NMF解混結(jié)果中的端元分布情況進(jìn)行修正，使其逐步向“空間能量”最小的理想分布狀態(tài)逼近。

綜上所述，MRF-NMF算法步驟可簡(jiǎn)單歸納如下：

整地，首先對(duì)地面雜物清理，清除雜草，保證苗圃地的清凈；其次對(duì)高洼地面采用推土機(jī)進(jìn)行平整，將土層較厚的土壤推到低洼的地方。

1)利用基于最小誤差的高光譜信號(hào)識(shí)別(Hyperspectral Signal identification by minimum error, HySime)法估算端元數(shù)量P；

2)基于估算出的端元數(shù)量P，結(jié)合頂點(diǎn)成分分析(Vertex Component Analysis, VCA)和全約束最小二乘法(Fully Constrained Least Square, FCLS)，初始化端元矩陣M和豐度矩陣S；

3)NMF步驟：利用標(biāo)準(zhǔn)NMF方法進(jìn)行迭代運(yùn)算，得到分離矩陣U的迭代結(jié)果；

4)MRF步驟：歸一化U的每一列，同時(shí)估計(jì)像元對(duì)應(yīng)特征向量的均值矩陣W，然后利用式(7)計(jì)算U；

5)重復(fù)步驟3)～4)，直到各自的停止準(zhǔn)則同時(shí)滿足，得到一個(gè)估計(jì)的成分，將其轉(zhuǎn)換為矩陣，即獲得一個(gè)端元分布。繼續(xù)迭代，直至滿足閾值條件，得到端元分布的估計(jì)結(jié)果。

步驟1)中采用的端元數(shù)量估計(jì)算法(HySime)，是一種估計(jì)高光譜信號(hào)子空間的方法，雖然計(jì)算過(guò)程復(fù)雜，但不需任何參數(shù)，具有自適應(yīng)性，估計(jì)準(zhǔn)確度較高。文獻(xiàn)[13]給出了該算法的原理與詳細(xì)實(shí)現(xiàn)過(guò)程。

步驟2)中用于初始化端元矩陣的VCA是一種端元識(shí)別算法，優(yōu)點(diǎn)是運(yùn)行時(shí)間短、效率高、精度較高，缺點(diǎn)是對(duì)純像元比較依賴。但文獻(xiàn)[14]實(shí)驗(yàn)結(jié)果表明，即使純像元不存在，VCA作為NMF的端元初始化算法的精度與綜合性能也是比較理想的。

從步驟3)和4)可以看出，MRF-NMF摒棄了在NMF目標(biāo)函數(shù)內(nèi)部添加新約束項(xiàng)的傳統(tǒng)做法，通過(guò)MRF構(gòu)建與NMF目標(biāo)函數(shù)交替運(yùn)行的獨(dú)立輔助函數(shù)，避免了內(nèi)部函數(shù)項(xiàng)之間的相互干擾。

下面討論該算法的收斂性。

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

實(shí)驗(yàn)運(yùn)行的軟硬件環(huán)境：操作系統(tǒng)Windows 7 Service Pack 1 64 b；處理器Intel Xeon E3- 1230 v3 3.30 GHz；內(nèi)存8.00 GB；硬盤Seagate 1 TB 7 200 rpm。

本文根據(jù)文獻(xiàn)[17]計(jì)算全局Moran’I指數(shù)，定量描述實(shí)驗(yàn)數(shù)據(jù)的空間相關(guān)程度，采用了兩組空間相關(guān)程度不同的真實(shí)高光譜數(shù)據(jù)進(jìn)行解混實(shí)驗(yàn)。全局Moran’I指數(shù)越高，表示圖像空間相關(guān)程度越高，反之則越低。

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

實(shí)驗(yàn)數(shù)據(jù)1選擇美國(guó)華盛頓特區(qū)的機(jī)載高光譜數(shù)字圖像收集實(shí)驗(yàn)儀器(HYperspectral Digital Imagery Collection Experiment, HYDICE)獲取的高光譜圖像，如圖1所示。該數(shù)據(jù)的全局Moran’I指數(shù)為0.498 6，說(shuō)明其空間相關(guān)程度較高。圖像行列數(shù)均為400，具有明顯的地物分布空間相關(guān)性；原始波段數(shù)量210，波段范圍0.4～2.4 μm，波段寬度為10 nm，涵蓋了可見光和近紅外譜段范圍；在去除了0.9～1.4 μm范圍內(nèi)的大氣吸收波段后，剩余191波段。

圖1 美國(guó)華盛頓特區(qū)HYDICE數(shù)據(jù)偽彩色圖像(R:60，G:25，B:15)Fig. 1 Pseudo-color image of HYDICE data in Washington, DC, USA (R:60，G:25，B:15)

通過(guò)目視解譯，圖像中主要含有植被、裸土、水體3種地物，其他類型地物還包括線狀道路、小型建筑物等。為簡(jiǎn)化描述，實(shí)驗(yàn)中不考慮面積比例很小的地物類型，即認(rèn)為所有像元中都只包含植被、裸土、水體3種主要地物。

由于HYDICE為機(jī)載光譜成像儀，空間分辨率高，可認(rèn)為圖像中每種端元都含有大量純像元。因此，本文通過(guò)在原圖像中人工選擇參考點(diǎn)的方式收集每種端元(地物)的光譜作為參考值，利用全約束最小二乘法計(jì)算端元的豐度參考值。其中，得到的端元參考光譜如圖2所示。由圖2可以看出，相對(duì)其余兩種地物，水體的光譜曲線差異最為顯著。由于實(shí)驗(yàn)數(shù)據(jù)為遙感影像像元亮度(Digital Number, DN)值，且缺少進(jìn)行反射率反演的必要參數(shù)，圖2中光譜曲線的y軸數(shù)值為DN值，而非地表反射率。

接下來(lái)，通過(guò)均值法重采樣使圖像的空間分辨率降低為原來(lái)的0.1倍，由于空間分辨率的大幅降低，形成了大量混合像元。然后，對(duì)生成的含有大量混合像元的新圖像進(jìn)行解混實(shí)驗(yàn)，以驗(yàn)證本文算法的效果和精度。由于在地物范圍確定的前提下，豐度分解結(jié)果的尺度是和空間分辨率大小嚴(yán)格對(duì)應(yīng)的，這里同樣利用均值法重采樣技術(shù)把之前求出的端元豐度參考值進(jìn)行聚合，得到行、列元素?cái)?shù)量均為原始圖像0.1倍的新圖像的豐度參考值，作為解混結(jié)果精度分析中的近似真值。重采樣后對(duì)應(yīng)的豐度參考值如圖3所示。

需要說(shuō)明的是，圖3中純白色代表端元在該像元內(nèi)部面積比例為1(100%)，純黑色代表端元在該像元內(nèi)部面積比例為0，其余各階灰度分別對(duì)應(yīng)0～1的不同比例。下文其他豐度圖像中，不同灰度代表的含義與圖3相同。

圖2 HYDICE數(shù)據(jù)端元參考光譜Fig. 2 Endmember reference spectrum of HYDICE data

圖3 HYDICE數(shù)據(jù)豐度參考值Fig. 3 Abundance reference value of HYDICE data

為充分展示本文方法在解混真實(shí)高光譜圖像時(shí)的效果，選擇前文提到的三種比較有代表性的NMF解混算法MVC-NMF、PSNMFSC和APS-NMF作為參考算法，與MRF-NMF進(jìn)行實(shí)驗(yàn)對(duì)比。與本文算法類似，MVC-NMF等三種算法均以不同的原理與形式，在NMF解混過(guò)程中引入了高光譜圖像的空間信息，且文獻(xiàn)[3～5]表明，三種解混算法均有效且具有較為理想的精度。

MRF-NMF豐度分解結(jié)果如圖4所示，其中三幅子圖中的白色分別代表了某種地物類型在圖像中的分布位置。由圖4可以看出，三者與圖3中分別對(duì)應(yīng)的真值基本吻合，說(shuō)明MRF-NMF能夠?qū)⑷N不同光譜特征的地物有效分離。

圖4 MRF-NMF的豐度分解結(jié)果Fig. 4 Decomposition results of abundances with MRF-NMF

表1和表2分別列出了四種算法的端元光譜和豐度分布分解結(jié)果的精度分析。其中，利用光譜角距離(Spectral Angle Distance, SAD)計(jì)算端元光譜估計(jì)精度，利用均方根誤差(Root Mean Square Error, RMSE)計(jì)算端元豐度估計(jì)精度。

SAD的計(jì)算公式為：

(11)

整幅圖像的RMSE為：

(12)

從表1～2可以看出，無(wú)論是端元光譜，還是端元豐度，MRF-NMF在四種算法中的分解精度都最高。其中，MRF-NMF的端元光譜分解精度相比MVC-NMF提高了10.6%，比PSNMFSC提高了12.3%，比APS-NMF提高了14.1%；MRF-NMF的端元豐度分解精度相比MVC-NMF提高了14.4%，比PSNMFSC提高了15.9%，比APS-NMF提高了15.3%。

表1 端元光譜分解結(jié)果精度Tab. 1 Precision of decomposition results of endmembers’ spectrum

表2 豐度分布估計(jì)結(jié)果精度Tab. 2 Precision of estimation results of abundance distribution

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

實(shí)驗(yàn)數(shù)據(jù)1的自然場(chǎng)景以成片分布的水、植被和裸土等自然場(chǎng)景為主，空間相關(guān)特征顯著。為比較和驗(yàn)證算法對(duì)空間相關(guān)程度不同的場(chǎng)景的解混效果，選擇獲取于1995年7月美國(guó)Nevada州Cuprite采礦區(qū)的機(jī)載可見光/紅外成像光譜儀(Airborne Visible InfRared Imaging Spectrometer, AVIRIS)高光譜數(shù)據(jù)作為實(shí)驗(yàn)數(shù)據(jù)2，如圖5所示。實(shí)驗(yàn)數(shù)據(jù)2的全局Moran’I指數(shù)為0.224 7，相對(duì)于實(shí)驗(yàn)數(shù)據(jù)1，其空間相關(guān)程度明顯降低，說(shuō)明礦區(qū)內(nèi)各種礦物的分布較為分散和雜亂。

圖5 Cuprite礦區(qū)AVIRIS高光譜數(shù)據(jù)Fig. 5 AVIRIS hyperspectral data of Cuprite mining field

圖5中的A、B、C、K、M等字樣分別代表該地區(qū)5種廣泛分布的礦物Alunite、Buddingtointe、Calcite、Kaolinite和Muscovite的大致分布位置。

該圖像大小為400列、350行，空間分辨率為20 m，波長(zhǎng)范圍為1.99～2.48 μm ，光譜分辨率為10 nm，共50波段，依次為AVIRIS原始譜段中的第172～221波段。該地區(qū)位于美國(guó)Nevada州南部，地表多為裸露礦物，基本無(wú)植被覆蓋。相對(duì)于真實(shí)數(shù)據(jù)實(shí)驗(yàn)1中的美國(guó)華盛頓特區(qū)HYDICE高光譜數(shù)據(jù)，該數(shù)據(jù)的空間相關(guān)特征不明顯，幾種主要礦物空間分布的隨機(jī)性較大。

由于AVIRIS為機(jī)載光譜成像儀，空間分辨率(20 m)相對(duì)較高，可合理假設(shè)圖像中每種端元都含有一定數(shù)量的純像元。因此，本文通過(guò)在原圖像中人工選擇參考點(diǎn)的方式，收集每種端元(地物)的光譜作為參考值，利用全約束最小二乘法計(jì)算端元的豐度參考值。

下面重復(fù)實(shí)驗(yàn)數(shù)據(jù)1的實(shí)驗(yàn)環(huán)節(jié)，以此對(duì)比和驗(yàn)證算法針對(duì)不同場(chǎng)景時(shí)的性能表現(xiàn)。需要說(shuō)明的是，實(shí)驗(yàn)只針對(duì)上述五種主要礦物類型進(jìn)行解混實(shí)驗(yàn)，對(duì)于其余小目標(biāo)對(duì)應(yīng)的地物類型的處理與實(shí)驗(yàn)數(shù)據(jù)1完全相同，具體情況不再贅述。

圖6為本文算法MRF-NMF對(duì)應(yīng)的Cuprite采礦區(qū)AVIRIS數(shù)據(jù)豐度估計(jì)結(jié)果。由于篇幅限制，其他參考算法的結(jié)果圖不再一一列出。

將圖6中每種地物的MRF-NMF豐度估計(jì)值與前文中利用全約束最小二乘法計(jì)算出的豐度參考值進(jìn)行比較，相似度均在85%以上，說(shuō)明MRF-NMF算法能夠有效分離出五種不同地物類型及其大體分布位置。

圖6 Cuprite礦區(qū)AVIRIS高光譜數(shù)據(jù)MRF-NMF豐度估計(jì)結(jié)果Fig. 6 Abundance estimation results of AVIRIS hyperspectral data of Cuprite mining field with MRF-NMF

表3～4分別給出了四種解混算法的端元光譜和豐度解混的精度分析結(jié)果。

表3 Cuprite采礦區(qū)數(shù)據(jù)的端元光譜分解結(jié)果精度Tab. 3 Precision of endmembers’ spectrum decomposition results of Cuprite mining field

表4 Cuprite采礦區(qū)數(shù)據(jù)的豐度分布估計(jì)結(jié)果精度Tab. 4 Precision of abundance distribution estimation results of Cuprite mining field

從表3～4中可以看出，四種解混算法均能有效分解出五種主要礦物類型，本文提出的MRF-NMF算法性能最好。以端元光譜分解結(jié)果為準(zhǔn)，MRF-NMF比MVC-NMF、PSNMFSC和APS-NMF三種算法的分解精度分別提高了7.82%、12.4%和10.1%；以豐度估計(jì)結(jié)果為準(zhǔn)，MRF-NMF比MVC-NMF、PSNMFSC和APS-NMF三種算法的分解精度分別提高了8.34%、12.6%和9.87%。

上述結(jié)果表明，MRF-NMF在實(shí)驗(yàn)數(shù)據(jù)空間相關(guān)程度顯著減小的情況下，相對(duì)于其他三種NMF代表性算法，仍具有較為明顯的精度優(yōu)勢(shì)，只是比實(shí)驗(yàn)數(shù)據(jù)空間相關(guān)程度顯著時(shí)略有下降。

最后，為了更全面地分析算法解混性能，給出了上述四種算法在Matlab 7.0環(huán)境下分別解混上述兩組空間相關(guān)程度不同的實(shí)驗(yàn)數(shù)據(jù)時(shí)運(yùn)行時(shí)間的統(tǒng)計(jì)與比較結(jié)果，如表5所示。其中，為便于和Cuprite礦區(qū)數(shù)據(jù)進(jìn)行比較，截取了美國(guó)華盛頓特區(qū)HYDICE數(shù)據(jù)中相同行列數(shù)(400列、350行)的子區(qū)作為實(shí)驗(yàn)數(shù)據(jù)；波段數(shù)也都統(tǒng)一選定為50(從美國(guó)華盛頓特區(qū)HYDICE數(shù)據(jù)中隨機(jī)選取50個(gè)波段)，同時(shí)將兩組實(shí)驗(yàn)數(shù)據(jù)的感興趣端元數(shù)量均設(shè)定為5(HYDICE數(shù)據(jù)中增加屋頂、道路等兩類小面積感興趣地物)。

從表5中可以看出，本文算法MRF-NMF的運(yùn)行時(shí)間要明顯少于其他三種算法，運(yùn)算效率最高；同時(shí)該算法在處理空間相關(guān)程度較低的數(shù)據(jù)時(shí)耗時(shí)有所增加。分析其原因，主要是由于空間相關(guān)程度的降低使得本文算法迭代次數(shù)增多。

表5 空間相關(guān)程度不同時(shí)不同算法運(yùn)行時(shí)間比較 sTab. 5 Running time comparison of different algorithms with different spatial correlation s

3 結(jié)語(yǔ)

本文針對(duì)大多數(shù)NMF擴(kuò)展方法中容易忽略的端元分布空間相關(guān)特征，利用MRF模型對(duì)其進(jìn)行描述，提出了一種新的基于NMF的高光譜解混方法。該方法通過(guò)構(gòu)建端元空間分布的能量函數(shù)來(lái)表征空間相關(guān)特征，作為與標(biāo)準(zhǔn)NMF目標(biāo)函數(shù)相互獨(dú)立的輔助函數(shù)，舍棄了構(gòu)造NMF目標(biāo)函數(shù)內(nèi)部約束項(xiàng)的傳統(tǒng)做法，為基于NMF的高光譜解混研究提供了新的思路。

實(shí)驗(yàn)結(jié)果表明，對(duì)于大部分真實(shí)高光譜數(shù)據(jù)，所提方法的分解精度和計(jì)算效率均優(yōu)于其他幾種代表性對(duì)比算法。

盡管空間相關(guān)特征在高光譜圖像中普遍存在，但并非全部真實(shí)數(shù)據(jù)都同時(shí)具有顯著的(相鄰像元)空間相關(guān)特征。實(shí)驗(yàn)結(jié)果也表明了當(dāng)空間相關(guān)程度不高時(shí)，本文算法性能會(huì)出現(xiàn)小幅下降?？梢灶A(yù)見，如果出現(xiàn)實(shí)驗(yàn)圖像的像元之間空間相關(guān)程度接近于零或負(fù)相關(guān)等極端情況，算法性能可能會(huì)出現(xiàn)較為嚴(yán)重的退化。如何進(jìn)一步穩(wěn)定算法性能，將是下一步研究工作的主要內(nèi)容。

References)

[1] LEE D D, SEUNG H S. Learning the parts of objects by non-negative matrix factorization [J]. Nature, 1999, 401: 788-791.

[2] PAUCA V P, PIPER J, PLEMMONS R J. Nonnegative matrix factorization for spectral data analysis [J]. Linear Algebra and its Applications, 2006, 416(1): 29-47.

[3] MIAO L D, QI H R. 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.

[4] ZYMNIS A, KIM S J, SKAF J, et al. Hyperspectral image unmixing via alternating projected subgradients [C]// ACSSC 2007: Proceedings of the 2007 Conference Record of the Forty-First Asilomar Conference on Signals, Systems & Computers. Piscataway, NJ: IEEE, 2007: 1164-1168.

[5] JIA S, QIAN Y T. Constrained nonnegative matrix factorization for hyperspectral unmixing [J]. IEEE Transactions on Geoscience and Remote Sensing, 2009, 47(l): 161-173.

[6] 杜世強(qiáng),石玉清,王維蘭,等.基于圖正則化的半監(jiān)督非負(fù)矩陣分解[J].計(jì)算機(jī)工程與應(yīng)用,2012,48(36):194-200.(DU S Q, SHI Y Q, WANG W L, et al. Graph regularized-based semi-supervised non-negative matrix factorization [J]. Computer Engineering and Applications, 2012, 48(36): 194-200.)

[7] 劉建軍,吳澤彬,韋志輝,等.基于空間相關(guān)性約束稀疏表示的高光譜圖像分類[J].電子與信息學(xué)報(bào),2012,34(11):2666-2671.(LIU J J, WU Z B, WEI Z H, et al. Spatial correlation constrained sparse representation for hyperspectral image classification [J]. Journal of Electronics & Information Technology, 2012, 34(11): 2666-2671.)

[8] CHEN B L, LI M, WANG J X, et al. Disease gene identification by using graph kernels and Markov random fields [J]. Science China Life Sciences, 2014, 57(11): 1054-1063.

[9] 姜小燕,孫福明,李豪杰.基于圖正則化和稀疏約束的半監(jiān)督非負(fù)矩陣分解[J].計(jì)算機(jī)科學(xué),2016,43(7):77-82.(JIANG X Y, SUN F M, LI H J. Semi-supervised nonnegative matrix factorization based on graph regularization and sparseness constraints [J]. Computer Science, 2016, 43(7): 77-82.)

[10] DOBRUSHIN P L. The description of a random field by means of conditional probabilities and conditions of its regularity [J]. Theory of Probability and its Applications, 1968, 13 (2): 197-224.

[11] SPITZER F. Markov random fields and gibbs ensembles [J]. American Mathematical Monthly, 1971, 78 (2): 142-154.

[12] DENG H W, CLAUSI D A. Unsupervised image segmentation using a simple MRF model with a new implementation scheme [J]. Pattern Recognition, 2004, 37(12): 2323-2335.

[13] BIOUCAS-DIAS J M, NASCIMENTO J M P. Hyperspectral subspace identification [J]. IEEE Transactions on Geoscience and Remote Sensing, 2008, 46(8): 2435-2445.

[14] 袁博.高光譜像元解混盲處理方法研究——面向盲分解的相關(guān)性分析與應(yīng)用[D].北京:中國(guó)科學(xué)院大學(xué),2015:45-48.(YUAN B. Research on blind processing algorithms of hyperspectral unmixing — correlation analysis and application of blind decoposition [D]. Beijing: University of Chinese Academy of Sciences, 2015: 45-48.)

[15] LEE D D, SEUNG H S. Algorithms for nonnegative matrix factorization [J]. Advances in Neural Information Processing Systems, 2001, 13(6): 556-562.

[16] 吳金彪.D-N交替迭代法及其收斂性分析[J].數(shù)值計(jì)算與計(jì)算機(jī)應(yīng)用,2002,23(2):121-130.(WU J B. D-N commutative method and its convergence [J]. Journal on Numerical Methods and Computer Applications, 2002, 23(2): 121-130.)

[17] 朱子明,祁新華.基于Moran’I的閩南三角洲空間發(fā)展研究[J].經(jīng)濟(jì)地理,2009,29(12):1977-1980.(ZHU Z M, QI X H. Research on spatial development of minnan delta basing on Moran’I [J]. Economic Geography, 2009, 29(12): 1977-1980.)

YUANBo, born in 1982, Ph. D., lecturer. His research interests include hyperspectral data processing, Internet of things engineering.

ApplicationofMRF’sspatialcorrelationmodelinNMF-basedlinearunmixing

YUAN Bo*

(CollegeofComputerandInformationEngineering,NanyangInstituteofTechnology,NanyangHenan473004,China)

Aiming at the problems of initialization and “l(fā)ocal minima” of Non-negative Matrix Factorization (NMF) in hyperspectral unmixing, a spatial correlation constrained NMF linear unmixing algorithm based on Markov Random Field (MRF) (MRF-NMF) was proposed. Firstly, the number of endmembers was estimated by Hyperspectral Signal identification by minimum error (HySime) method, the endmember matrix and abundance matrix were initialized by Vertex Component Analysis (VCA) and Fully Constrained Least Squares (FCLS). Secondly, the energy function of depicting the spatial distribution characteristics of ground objects was established by using MRF to depict the spatial correlation distribution features of ground objects. Finally, the spatial correlation constraint function based on MRF and the NMF standard objective function were used for unmixing in the form of alternating iteration, and the endmember information and abundance decomposition results of hyperspectral data were obtained. The theoretical analysis and experimental results of real data show that, with hyperspectral data of low spatial correlation, compared with the three reference algorithms of Minimum Volume Constrained NMF(MVC-NMF), Piecewise Smoothness NMF with Sparseness Constraints (PSNMFSC) and NMF with Alternating Projected Subgradients (APS-NMF), the endmember decomposition precision of MRF-NMF increases by 7.82%, 12.4% and 10.1%, and the abundance decomposition precision of MRF-NMF increases by 8.34%, 12.6% and 9.87%. The proposed MRF-NMF can make up for NMF’s deficiency in depicting spatial correlation features, and reduce the spatial energy distribution error of ground objects.

Non-negative Matrix Factorization (NMF); hyperspectral linear unmixing; spatial correlation; Markov Random Field (MRF); alternative iteration; spatial energy

2017- 06- 12;

2017- 08- 29。

袁博(1982—)，男，河南南陽(yáng)人，講師，博士，主要研究方向：高光譜數(shù)據(jù)處理、物聯(lián)網(wǎng)工程。

1001- 9081(2017)12- 3563- 06

10.11772/j.issn.1001- 9081.2017.12.3563

(*通信作者電子郵箱nylgyb@163.com)

TP751.1

A