劉萬軍　楊秀紅　曲海成,2　孟　煜　姜慶玲

1(遼寧工程技術(shù)大學軟件學院　遼寧 葫蘆島 125105)2(哈爾濱工業(yè)大學電子與信息工程學院　黑龍江 哈爾濱 150006)3(鐵嶺師范高等?？茖W?！∵|寧 鐵嶺 112000)

?

基于端元子集優(yōu)選的高光譜解混算法研究

劉萬軍1楊秀紅1曲海成1,2孟煜1姜慶玲3

1(遼寧工程技術(shù)大學軟件學院遼寧 葫蘆島 125105)2(哈爾濱工業(yè)大學電子與信息工程學院黑龍江 哈爾濱 150006)3(鐵嶺師范高等?？茖W校遼寧 鐵嶺 112000)

摘要針對采用最大體積單體MVS(Maximization Volume Simplex)端元提取算法進行端元初選時存在相似端元光譜問題，提出一種光譜信息散度SID(Spectral Information Divergence)和光譜梯度角SGA(Spectral Gradient Angle)相結(jié)合以區(qū)分兩個相似端元光譜的方法。該方法對經(jīng)過端元初選之后的端元子集進行端元的二次選擇，采用以SID_SG作為最相似端元選擇的判據(jù)，除去相似端元，降低相似端元對解混精度的影響，利用全約束最小二乘法進行豐度估計。實驗結(jié)果表明，提出的優(yōu)化方法與傳統(tǒng)方法相比，提高了端元的選擇精度，重構(gòu)影像與原始影像之間的均方根誤差RMSE(Root Mean Square Error)也有所降低，分布更加均勻。該方法對高光譜遙感影像進行深度解譯具有十分重要的意義。

關(guān)鍵詞光譜解混端元初選二次提取除去端元解混算法

0引言

高光譜遙感技術(shù)已經(jīng)廣泛應(yīng)用于軍事和民用等各個領(lǐng)域。然而，阻礙高光譜遙感技術(shù)向定量化方向深入發(fā)展的主要因素是混合像元的普遍存在[1]。因為高光譜傳感器有限的空間分辨率，每個像元對應(yīng)著較大的區(qū)域，所以每個像元一般包含多種物質(zhì)，多種物質(zhì)光譜混合的混合光譜就是傳感器探測的光譜[2]。解決混合像元問題最有效的方法就是混合像元分解[3]。高光譜混合像元分解一般包括兩個方面，分別是端元提取[4]和豐度估計[5]。當前端元提取方法有N-FINDR[6]、頂點成分分析[7]VCA(Vertex Component Analysis)、迭代誤差分析IEA(Iterative Error Analysis)[8]和純像元指數(shù)PPI(Pure Pixel Index)等。傳統(tǒng)的豐度估計算法主要有非限制性最小二乘法，非負限制性最小二乘法，和為一限制性最小二乘法和全限制性最小二乘法[9]。在進行高光譜解混時需要盡可能地選擇出最精確的端元，降低解混時出現(xiàn)的誤差。

本文結(jié)合兩種常用的光譜區(qū)分方法，提出了一種將光譜信息散度和光譜梯度角相結(jié)合的方法用以區(qū)分兩個相似端元光譜。該方法不僅能夠從形狀上反映兩條光譜曲線的匹配程度，而且能夠?qū)蓷l有較小差異的光譜曲線做出較大的響應(yīng)。首先利用最大體積單體端元提取算法進行端元的初次提取，然后再利用光譜信息散度和光譜梯度角SID_SG(Spectral Information Divergence_ Spectral Gradient)相結(jié)合的算法找出端元子集中存在的相似端元，去除相似端元，最后利用全約束最小二乘法FCLS(Full Constraint Least Square)進行端元的豐度估計，從而實現(xiàn)高光譜的混合像元分解。

1高光譜解混算法

1.1最大體積單體MVS

一般而言，高光譜圖像中的每個像元都能看作是圖像中各端元的線性混合的像元:

(1)

(2)

0≤ki≤1

(3)

式中，x是高光譜圖像中任一P維的光譜向量(P為圖像的波段數(shù)目)，E=[e1,e2,…,en]是P×N的矩陣，N是端元的數(shù)目，E中的每一列都是端元向量。c=(c1,c2,…,cn)T是系數(shù)向量，ci是像元x中端元ei所占比例，n表示誤差項。

圖1　兩個波段三個端元

在誤差項n非常小時，只有滿足上面3個式子的點的集合，才能夠構(gòu)成高維空間的凸集，這些端元就坐落在凸面的單形體的各個頂點[10]。這里采用兩個波段、三個端元作為例子，從而清晰地表明它們之間存在的幾何關(guān)系，如圖1所示。從圖1中看出,三角形體的各個頂點分別是端元a、b和c，那么三角形內(nèi)部的點就對應(yīng)著圖像中混合像元。所以提取高光譜圖像的端元問題就轉(zhuǎn)化成求單形體的頂點問題[11]。

從圖1可以看出，圖1中的端元a，端元b以及端元c所構(gòu)成的三角形面積是圖像中的任意三個像元所構(gòu)成的面積集合中的最大的一個，表示如下:

S(a,b,c)=max{S(i,j,k)}

(4)

式中，S(i,j,k)表示的是由像元i、j和k所構(gòu)成的三角形面積。所以尋找單形體的頂點的問題就轉(zhuǎn)化為求最大體積單體的問題，體積公式如下:

(5)

(6)

式中，ei表示第i個端元的列向量，V表示的是由e1,e2,…,en這N個端元所構(gòu)成的單形體體積，|·|表示行列式的運算符。因為這里用到了求行列式的運算，所以要求E必須為方陣。

1.2光譜信息散度SID

Ii(x)=-lgai

(7)

Ii(y)=-lgbi

(8)

通過式(7)和式(8)，可得y關(guān)于x的相對熵：

(9)

同樣可以得到x關(guān)于y的相對熵：

(10)

因而x和y的光譜信息散度為：

SID(x，y)=D(x‖y)+D(y‖x)

(11)

光譜信息散度法的優(yōu)點是可以對兩條光譜進行整體上的比較。

1.3光譜梯度角SGA

光譜梯度角是在對兩條光譜向量一階求導的基礎(chǔ)上求解它們的梯度向量，然后再得出兩個梯度向量間的廣義夾角大小[13]。

光譜曲線x的梯度向量為：

sGA(x)=[x2-x1,x3-x2,…,xn-xn-1]

(12)

光譜曲線y的梯度向量為：

sGA(y)=[y2-y1,y3-y2,…,yn-yn-1]

(13)

則梯度向量的廣義夾角大小為：

(14)

該方法的優(yōu)勢在于能夠反映光譜曲線的斜率變化，能反映出光譜曲線的局部特征的變化。

1.4信息散度與梯度角相結(jié)合

兩條光譜曲線的銳角夾角范圍是(0,π/2)，并且梯度向量的銳角夾角范圍同樣是(0,π/2)。因為正切函數(shù)的特性是當夾角是π/4時值為1；當夾角大于π/4時，其值隨著夾角的增大迅速變大。根據(jù)正切的函數(shù)特征，將光譜梯度夾角控制在[π/4，π/2]范圍內(nèi)，當夾角值改變較小時，其結(jié)果就能有明顯變化。因此為了不僅可以在整體上比較兩條光譜曲線的匹配程度，而且可以從局部上區(qū)分兩條光譜曲線的差異，提出將SID與SGA相結(jié)合的公式如下：

(15)

1.5均方根誤差RMSE

RMSE[14]是均方根誤差，是一個度量準則，用來評價兩幅圖像的相似性，均方根誤差越小，說明二者越相近。

(16)

式中，假設(shè)端元豐度矩陣為S=(s11,s12,…,s1n;…;sp1,sp2,…,spn),真實豐度矩陣為Z=(z11,z12,…,z1n;…;zp1,zp2,…,zpn),sij和zij分別為矩陣S和Z的第i行第j列元素。

1.6豐度估計誤差

(17)

式中，m是圖像中混合光譜條數(shù)；n是端元組中端元個數(shù)；qij是第i條光譜第j個端元的計算豐度值；pij是第i條光譜第j個端元的假定豐度值；fa是豐度估計誤差。

2基于端元子集優(yōu)選的高光譜解混

利用最大體積單體端元提取算法進行端元的初次提取，先去除部分端元，再進行端元的二次提取迭代循環(huán)，去掉相似的端元，剩下的端元為最佳端元，最后再利用全約束最小二乘法進行豐度的估計，最終實現(xiàn)高光譜的混合像元分解。端元初次提取的目的是為了去除部分多余的端元，端元二次選擇的目的是進一步從原始高光譜影像中找到純像元，從而準確求解各純物質(zhì)在混合像元中所占的比例。以SID_SG作為最相似的端元選擇準則，能夠選擇出最相似的端元組，去除均方根誤差較大的端元，保留均方根誤差較小的端元，從而降低相似端元對解混精度的影響。

算法流程如圖2所示。

圖2　基于端元子集優(yōu)選的高光譜解混算法流程圖

端元二次選擇優(yōu)選算法步驟：

步驟1找出一組相似端元。將端元集合利用光譜匹配算法SID_SG進行第一次循環(huán)，找出最相似的一組端元，分別是端元A和端元B。

步驟2計算出去除端元A后的均方根誤差。去除端元集合中的端元A，利用剩下的端元計算出均方根誤差。

步驟3計算出去除端元B后的均方根誤差。去除端元集合中的端元B，利用剩下的端元計算出均方根誤差。

步驟4比較步驟2和步驟3中的均方根誤差。將步驟2和步驟3中計算出的均方根誤差進行對比，留下均方根誤差比較小的端元。

步驟5重復步驟1～4，直至選出均方根誤差最小并且求解豐度所運行的時間是最小的端元組，最終形成新的端元集合。

選出優(yōu)化后的端元組之后，再采用全約束最小二乘法求解高光譜最終的豐度。

3實驗結(jié)果分析

3.1端元數(shù)目不同時的對比實驗

該實驗采用的模擬圖像大小是64×64，端元數(shù)目是3至8個，并且添加的信噪比SNR(Signal to Noise Ratio)為30 dB。在端元數(shù)目不同的情況下對比N-FINDR、單體增長算法SGA(Simplex Growing Algorithm)和最大體積單體(MVS)這三種端元提取算法的優(yōu)劣，見表1所示。

表1　端元數(shù)不同時的實驗對比

從表1中可以看出，N-FINDR、單體增長算法(SGA)、最大體積單體(MVS)這三種算法在端元數(shù)目不同的情況下對端元提取的結(jié)果進行相似性度量的對比。隨著端元數(shù)目的增多， MVS算法的效果是最好的，SGA算法要優(yōu)于N-FINDR算法，所以采用MVS算法進行端元初次選擇的效果要優(yōu)于SGA算法和N-FINDR算法。

3.2信噪比不同時的對比實驗

該實驗采用的模擬圖像是64×64，端元數(shù)目是4條，并且添加的信噪比為10至50 dB。在不同的信噪比下對比N-FINDR，單體增長算法(SGA)，最大體積單體(MVS)這三種端元提取算法的魯棒性，見表2所示。

表2　信噪比不同時的實驗對比

從表2中可以看出，N-FINDR、單體增長算法(SGA)、最大體積單體(MVS)這三種算法在信噪比不同的情況下對端元提取的結(jié)果進行相似性度量的對比。對比表2中數(shù)據(jù)發(fā)現(xiàn)N-FINDR、SGA和MVS這三種算法的魯棒性從高到低依次是MVS算法、SGA算法、N-FINDR算法，所以采用MVS算法進行端元的初次選擇的穩(wěn)定性要優(yōu)于SGA算法和N-FINDR算法。

3.3分類精度對比實驗

實驗數(shù)據(jù)集為高光譜研究中經(jīng)常采用的Pavia數(shù)據(jù)[15]，該研究區(qū)域位于意大利北部的帕維亞，該地區(qū)地形復雜，主要有灌木、水體、森林、農(nóng)田、草地和城鎮(zhèn)等地物類型。

為了更好地證明SID_SG方法的優(yōu)越性，利用實驗數(shù)據(jù)將SID_SG方法與光譜角匹配SAM(Spectral Angle Mapping)法，光譜信息散度(SID)法，光譜梯度角(SGA)法這三種方法做對比實驗，將它們在分類精度上進行對比驗證，見表3所示。

表3　不同分類方法的精度分析

從表3中可以看出，分類精度最高的是SID_SG，達到78.96%，它對所調(diào)查的樣本點正確分類的個數(shù)是30個。分類精度最低的是SAM，它的精度是60.51%，它對所調(diào)查的樣本點正確分類的個數(shù)是24個。SID的分類精度低于SID_SG的分類精度，高于SAM的分類精度，它對所調(diào)查的樣本點正確分類的個數(shù)是27個，其精度是69.75%。SGA的分類精度也低于SID_SG的分類精度，高于SAM的分類精度，它對所調(diào)查的樣本點正確分類的個數(shù)是25個，其精度是62.81%。從這些數(shù)據(jù)可以看出，SID_SG的分類精度要高于SAM、SGA和SID，所以采用SID_SG算法提取相似端元的效果要優(yōu)于SAM、SGA和SID。

3.4均方根誤差及豐度估計時間對比

實驗數(shù)據(jù)集為高光譜解混研究中經(jīng)常采用的美國 Nevada 州南部沙漠地區(qū)的Cuprite數(shù)據(jù)[16]。該數(shù)據(jù)中共有224個波段，波長范圍為0.37到2.48 μm，光譜分辨率是10 nm。

RMSE是均方根誤差，是一個度量準則，用來評價兩幅圖像的相似性，均方根誤差越小，說明二者越相近。

利用實驗數(shù)據(jù)，根據(jù)本文優(yōu)選算法分別提取出不同的端元數(shù)，得到的RMSE隨著端元個數(shù)增加的變化情況，如圖3所示。利用全約束最小二乘法解混求豐度運行的時間，見表4所示。根據(jù)圖3中不同端元數(shù)目得出的均方根誤差值以及表4中給出的在不同端元數(shù)目下利用全約束最小二乘法解混求解豐度所運行的時間值可以看出去除十個相似的端元是最適宜的，所以在端元二次提取時提取出12個端元是最合適的。

圖3　端元數(shù)目和RMSE值的關(guān)系

端元數(shù)目求解豐度所運行時間(s)221156211097208681969318668176341660315578145421352112503

3.5整體豐度估計誤差對比

為了研究噪聲對本文算法性能的影響, 本文在混合光譜中加入了隨機噪聲。加入噪聲的公式為：

(18)

式中，r′是加噪聲之后的光譜；r是無噪聲的光譜；u(0,1)是均值為0，方差為1的隨機噪聲；SNR是信噪比；X是噪聲的假想發(fā)射率，假定為0.5。

本文實驗中將本文算法、MVS算法以及FCLS算法進行對比，這三種算法的豐度估計誤差隨信噪比的變化情況，如圖4所示。

圖4　三種算法的豐度估計誤差隨信噪比變化情況

從圖4中可以看出, 當信噪比SNR的范圍是[12，100]時, 本文算法豐度估計誤差fa是最小的，并且當信噪比達到100時，本文算法的豐度估計誤差小于0.1，綜上所述，本文算法要優(yōu)于最大體積單體算法和全約束最小二乘算法。

3.6優(yōu)選前后端元及誤差分布對比

利用Cuprite數(shù)據(jù)進行優(yōu)選前后的誤差分布圖對比。在端元初選時，利用最大體積單體端元提取算法進行端元的一次提取處理，提取出22個端元，得到其誤差分布圖，如圖5所示。對22個端元進行端元的二次提取處理后，得到12個端元的誤差分布圖，如圖6所示。

圖5　初次提取誤差分布圖　　圖6　二次提取誤差分布圖

對比圖5和圖6，可以看出優(yōu)選后誤差分布更加均勻，降低了誤差，提高了端元選擇的精度，減少了相似端元對解混精度的影響。

4結(jié)語

本文提出了一種將光譜信息散度和光譜梯度角相結(jié)合的方法用以區(qū)分兩個相似端元光譜。在利用最大體積單體端元提取算法進行端元初選之后，再采用以SID_SG作為最相似端元選擇的判據(jù)進行端元的二次選擇，除去相似的端元，降低相似端元對解混精度的影響。從實驗的結(jié)果可以看出, 提出的基于端元子集優(yōu)選的高光譜解混算法比傳統(tǒng)算法提高了端元的選擇精度，分布也更加均勻，同時也減少了重構(gòu)影像與原始影像之間的均方根誤差(RMSE)。該方法對高光譜遙感影像進行深度解譯具有十分重要的意義。

參考文獻

[1] Chen G,Song F,Wang J,et al.FRET spectral unmixing:a ratiometric fluorescent nanoprobe for hypochlorite[J].Chemical Communications,2011,24(4):2949-2951.

[2] Roth K L,Dennison P E,Roberts D A.Comparing endmember selection techniques for accurate mapping of plant species and land cover using imaging spectrometer data[J].Remote Sensing of Environment,2012,32(8):139-152.

[3] Gu Y,Wang S,Jia X.Spectral unmixing in multiple-kernel hilbert space for hyperspectral imagery[J].Gon and Rmo Nng Ranaon on,2013,51(10):3968-3981.

[4] Martin G,Plaza A.Spatial-spectral preprocessing prior to endmember identification and unmixing of remotely sensed hyperspectral data[J].Ld O N Ald Arh Obrvaon and Rmo Nng Jornal of,2012,5(12):380-395.

[5] Gonzalez C,Resano J,Plaza A,et al.FPGA implementation of abundance estimation for spectral unmixing of hyperspectral data using the image space reconstruction algorithm[J].Ld O N Ald Arh Obrvaon and Rmo Nng Jornal of,2012,5(3):248-261.

[6] Xiong W,Chang C,Wu C,et al.Fast algorithms to implement N-FINDR for hyperspectral endmember extraction[J].Ld O N Ald Arh Obrvaon and Rmo Nng Jornal of,2011,4(3):545-564.

[7] Krafft C,Diderhoshan M A,Recknagel P,et al.Crisp and soft multivariate methods visualize individual cell nuclei in Raman images of liver tissue sections[J].Vibrational Spectroscopy,2011,55(1):90-100.

[8] Fang Y,Yan H.Error analysis for remote nonlinear iterative learning control system with wireless channel noise[J].Journal of Shanghai University (English Edition),2011,15(3):7-11.

[9] Chang C I.An information-theoretic approach to spectral variability,similarity,and discrimination for hyperspectral image analysis[J].IEEE Transactionson InformationTheory,2000,45(5):1927-1932.

[10] 耿修瑞,張兵,張霞,等.一種基于高維空間凸面單行體體積的高光譜圖像解混算法[J].自然科學進展,2004,7(14):810-814.

[11] 耿修瑞,趙永超,周冠華.一種利用單行體體積自動提取高光譜圖像端元的算法[J].自然科學進展,2006,9(16):1196-1200.

[12] Jarosik N,Bennett C L,Dunkley J,et al.Seven-year wilkinson microwave anisotropy probe (WMAP) observations: sky maps,systematic errors,and basic results[J].The Astrophysical Journal Supplement Series,2011,192(2):14.

[13] Liu J,Zhang J.A new maximum simplex volume method based on householder transformation for endmember extraction[J].Ranaon on Gon and Rmo Nng,2012,50(6):104-118.

[14] Chai T,Draxler R R.Root mean square error (RMSE) or mean absolute error (MAE)?[J].Geoscientific Model Development Discussions,2014,7(1):1525-1534.

[15] Brasco L,Pratelli A.Sharp stability of some spectral inequalities[J].Geometric and Functional Analysis,2012,22(1):107-135.

[16] Chen B,Vodacek A,Cahill N D.A novel adaptive scheme for evaluating spectral similarity in high-resolution urban scenes[J].Selected Topics in Applied Earth Observations and Remote Sensing,IEEE Journal of,2013,9(3):1376-1385.

收稿日期：2014-09-28。國家高技術(shù)研究發(fā)展計劃項目(2012AA 12A405)；國家自然科學基金項目(61172144)。劉萬軍，教授，主研領(lǐng)域：數(shù)字圖像處理，運動目標檢測與跟蹤。楊秀紅，碩士。曲海成，博士。孟煜，碩士。姜慶玲，博士。

中圖分類號TP75

文獻標識碼A

DOI:10.3969/j.issn.1000-386x.2016.07.058

HYPERSPECTRAL UNMIXING ALGORITHM BASED ON ENDMEMBER SUBSET OPTIMISATION

Liu Wanjun1Yang Xiuhong1Qu Haicheng1,2Meng Yu1Jiang Qingling3

1(SchoolofSoftware,LiaoningTechnicalUniversity,Huludao125105,Liaoning,China)2(SchoolofElectronicsandInformationEngineering,HarbinInstituteofTechnology,Harbin150006,Heilongjiang,China)3(TielingNormalCollege,Tieling112000,Liaoning,China)

AbstractFor the problem of maximisation volume simplex (MVS)endmember extraction algorithm that when applying in primary endmembers selection it would have similar endmembers spectra, we proposed an algorithm which distinguishes two similar endmember spectra by combining the spectral information divergence (SID) and the spectral gradient angle (SGA). This algorithm carries out secondary selection on the endmember subset derived from primary endmembesr selection, and adopts SID_SG rule as the criteria for selecting the most similar endmembers to remove the similar endmembers and to reduce the influence of similar endmembers on unmixing accuracy, and uses the full constraint least square for abundance estimation. Experimental results showed that proposed optimisation algorithm improved the accuracy of endmember selection than the traditional method, the root mean square error (RMSE) between reconstruction images and original image was rather reduced, and the distribution was more evenly as well. The algorithm is of very great significance to the deep interpretation on hyperspectral remote sensing image.

KeywordsSpectral unmixingPrimary endmembers selectionSecondary extractionRemoval of endmemberUnmixing algorithm