趙　巖,　張春晶,　王東輝

(1.黑龍江科技大學(xué) 電氣與控制工程學(xué)院, 哈爾濱 150022; 2.黑龍江科技大學(xué) 電子與信息工程學(xué)院,哈爾濱 150022; 3.哈爾濱工程大學(xué) 信息與通信工程學(xué)院, 哈爾濱 150001)

?

基于廣義形態(tài)學(xué)的高光譜圖像解混算法

趙巖1,張春晶2,王東輝3

(1.黑龍江科技大學(xué) 電氣與控制工程學(xué)院, 哈爾濱 150022; 2.黑龍江科技大學(xué) 電子與信息工程學(xué)院,哈爾濱 150022; 3.哈爾濱工程大學(xué) 信息與通信工程學(xué)院, 哈爾濱 150001)

綜合利用像元的光譜和空間信息,提出一種基于廣義形態(tài)學(xué)的混合像元分解算法。引入基準(zhǔn)像元,避免形態(tài)學(xué)算子在像元排序規(guī)則和替換準(zhǔn)則上存在的局限性。針對(duì)圖像中不同類別交界處存在的交叉替換現(xiàn)象,廣義開-閉算子采用修正能量函數(shù)作為距離測(cè)度進(jìn)行計(jì)算。端元提取之后利用全約束最小二乘法進(jìn)行豐度反演。利用礦物光譜的模擬數(shù)據(jù)和真實(shí)數(shù)據(jù)進(jìn)行實(shí)驗(yàn),結(jié)果表明,該算法無(wú)需先驗(yàn)信息便可自動(dòng)進(jìn)行混合像元分解,分解精度較高。

高光譜混合像元分解; 廣義形態(tài)學(xué); 修正能量函數(shù)

地球表面很少由均一物質(zhì)組成,當(dāng)具有不同波譜屬性的物質(zhì)出現(xiàn)在同一個(gè)像元內(nèi)時(shí),就會(huì)出現(xiàn)混合現(xiàn)象,即混合像元(Mixed Pixel)[1-4]。高光譜圖像混合像元分解技術(shù)可應(yīng)用于礦產(chǎn)資源勘探等領(lǐng)域。在線性混合模型[5]進(jìn)行光譜解混操作之前,需要進(jìn)行光譜端元提取,為光譜解混獲取必要的先驗(yàn)知識(shí)。光譜端元的選擇應(yīng)當(dāng)具有代表性,成為圖像內(nèi)大多數(shù)像元的類別成分集合[6]。端元波譜作為混合像元分解過程中的參考波譜,影響混合像元分解結(jié)果的精度[7]。常用的獲取端元波譜的方法有純凈像元指數(shù)(Pixel Purity Index, PPI)、N-FINDR算法、頂點(diǎn)成分分析(Vector Component Analysis, VCA)、迭代誤差分析(Iterative Error Analysis, IEA)、獨(dú)立成分分析(Independent Component Analysis, ICA)等。這些方法只是考慮光譜信息,而忽略了像素間的空間相關(guān)性[8]。自動(dòng)形態(tài)學(xué)端元提取(Automated Morphological Endmember Extraction, AMEE)[9],在高光譜圖像處理中引入反映像素間空間相關(guān)性信息的形態(tài)學(xué)算子。AMEE算法膨脹運(yùn)算應(yīng)用到圖像中不同地物類別的邊界處時(shí),出現(xiàn)交叉替換現(xiàn)象。筆者引入正則化的修正能量函數(shù)作為距離測(cè)度,同時(shí)給出關(guān)鍵參量——基準(zhǔn)像元的計(jì)算方法,通過計(jì)算廣義形態(tài)開-閉算子進(jìn)行端元提取,之后進(jìn)行全約束最小二乘法[10]豐度反演,名為廣義形態(tài)學(xué)端元提取-全約束最小二乘算法(Generalized Morphological Endmember Extraction-Fully Constrainded Least Squares, GMEE-FCLS)。

1　廣義形態(tài)學(xué)算子

1.1線性光譜混合模型

線性光譜混合模型(Linear Spectral Mixture Model,LSMM)計(jì)算簡(jiǎn)單,應(yīng)用廣泛,是混合像元分解的常用方法。設(shè)X為像元光譜向量,S=[s1,s2,…,sN]為端元光譜矩陣,A=[a1,a2,…,aN]T是N維向量,其各分量元素為端元豐度,N為隨機(jī)噪聲。線性光譜混合模型為:

X=SA+N,

(1)

ai≥0,

(2)

(3)

1.2能量函數(shù)

能量函數(shù)(Energy Function)為相鄰像元的能量測(cè)度[11],它包含圖像的空間信息?？紤]N維光譜向量si=[si1,si2,…,siN]T,sj=[sj1,sj2,…,sjN]T,則si和sj的EF為

(4)

1.3修正能量函數(shù)

常用的表征光譜間相似性的距離測(cè)度有光譜角(Spectral Angle Distance, SAD)、歐式距離(Euclidean Minimum Distance, EMD)、光譜信息散度(Spectral Information Divergence, SID)、馬氏距離(Mahalanobis Distance, MD)等。文中引入修正能量函數(shù)(Modified Energy Function, MEF)計(jì)算該距離。將EF中的歐式距離d2替換為加入正則化項(xiàng)的馬氏距離d3。將正則化項(xiàng)加入?yún)f(xié)方差矩陣中,可以提高求取協(xié)方差矩陣的逆陣的穩(wěn)定性,使其適用于稍小的樣本[12]。馬氏距離可以排除光譜波段之間的相關(guān)性的干擾。修正能量函數(shù)(ECF)為

(5)

結(jié)構(gòu)元素的尺寸和結(jié)構(gòu)形狀影響形態(tài)學(xué)圖像空間信息的提取[13]。廣義形態(tài)學(xué)算子采用兩個(gè)不同尺寸結(jié)構(gòu)元素,輸出統(tǒng)計(jì)偏倚小于經(jīng)典形態(tài)學(xué)算子[14],端元提取性能得到改善。基于導(dǎo)出序的廣義形態(tài)學(xué)算子將光譜向量作為整體處理[15],廣義開-閉算子(GOC)和廣義閉-開算子(GCO)[16]為

GOC(f(x,y))=(f°B1·B2)(x,y),

(6)

GCO(f(x,y))=(f·B1°B2)(x,y),

(7)

通常情況下,純像元的GOC值最大,GCO值約為零,而混合像元的GCO值最大,GOC值約為零。

1.5基準(zhǔn)像元

(8)

其中,度量算子d(f(x,y),Bi)=ECF(f(x,y),e(j,k)),e(j,k)為基準(zhǔn)像元,

e(j,k)=argmax(GCO(f(x,y))。

(9)

2　端元提取算法與豐度反演

2.1基于廣義形態(tài)學(xué)算子的端元提取算法

基于廣義形態(tài)學(xué)算子的端元提取步驟:

步驟1將利用虛擬維度法(Virtual Dimensionality, VD)[18]估計(jì)出的端元數(shù)作為算法的終止條件。當(dāng)已提取的端元數(shù)達(dá)到VD估計(jì)的端元數(shù)時(shí),結(jié)束運(yùn)算。

步驟2設(shè)定結(jié)構(gòu)元素B1、B2的空間尺寸。

步驟3遍歷整幅圖像計(jì)算基準(zhǔn)像元e(j,k)。

步驟7若未滿足步驟1的終止條件,重復(fù)步驟4。

2.2全約束最小二乘法豐度反演

采用滿足“非負(fù)”和“和為1”的全約束最小二乘法進(jìn)行豐度反演。

(10)

將式(10)代入“非負(fù)”約束最小二乘法[8](Nonegatively Constrainded Least Squares, NCLS)得到全約束最小二乘法的豐度解。

3　實(shí)驗(yàn)分析

3.1性能指標(biāo)

解混精度評(píng)價(jià)采用光譜角距離(dSAD)和均方根誤差(Re)這兩個(gè)常用的指標(biāo)。

對(duì)于第i個(gè)端元,dSAD定義為

(11)

Re定義為

(12)

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

模擬實(shí)驗(yàn)數(shù)據(jù)由USGS礦物光譜庫(kù)中選取出五種線性獨(dú)立的端元光譜,按照Dirichlet分布進(jìn)行混合。端元豐度之和進(jìn)行歸一化處理,加入不同強(qiáng)度的白噪聲,波段數(shù)為50,圖像大小為256×256像元。

圖1　不同噪聲強(qiáng)度時(shí)不同算法性能的比較

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

真實(shí)圖像采用美國(guó)內(nèi)華達(dá)州Cuprite地區(qū)的AVIRIS(AirborneVisible/InfraredImagingSpectrometer)高光譜數(shù)據(jù)(400×350像元,波長(zhǎng)范圍1.99～2.48μm,成像于1995年),共計(jì)172～221波段間的50個(gè)波段數(shù)據(jù),圖2為第196波段。實(shí)驗(yàn)中去掉大氣吸收的221波段,實(shí)際采用49個(gè)波段。

圖2　內(nèi)華達(dá)州Cuprite地區(qū)的AVIRIS圖像

表1　Cuprite地區(qū)數(shù)據(jù)的光譜角距離比較

圖3　基于GMEE-FCLS算法的解混結(jié)果

4　結(jié)束語(yǔ)

基于廣義形態(tài)學(xué)算子的高光譜混合像元分解方法,從結(jié)構(gòu)元素中像元的排序規(guī)則和替換準(zhǔn)則入手,引入基準(zhǔn)像元,并以修正能量函數(shù)為距離測(cè)度。該GMEE-FCLS算法結(jié)合光譜和空間信息,考慮光譜波段的相關(guān)性。礦物光譜的模擬和真實(shí)數(shù)據(jù)實(shí)驗(yàn)分析表明，文中提出的方法能夠提高混合像元分解的精度。

[1]GREER J B. Sparse demixing of hyperspectral images[J]. IEEE Transactions on Image Processing, 2012, 21(1): 219-228.

[2]YOKOYA N, CHANUSSOT J, IWASAKI A. Nonlinear unmixing of hyperspectral data using semi-nonnegative matrix factorization[J]. IEEE Transactions on Geoscience and Remote Sensing, 2014 52(2): 1430-1437.

[3]趙春暉, 崔士玲, 劉務(wù). 基于分層的多端元光譜解混算法[J]. 光電子·激光, 2014, 25(9): 1830-1836.

[4]楊可明, 王林偉, 劉士文, 等. 壓縮感知和萬(wàn)有引力模型相結(jié)合的高光譜混合像元分解[J]. 測(cè)繪學(xué)報(bào), 2014, 43(10): 1068-1074.

[5]NEVILLE R A, STAENZ K, SZEREDI T, et al. Automatic endmember extraction from hyperspectral data for mineral exploration[C]//Fourth International Airborne Remote Sensing Conference and Exhibition. Canada： [s.n.], 1999: 21-24.

[6]王立國(guó), 趙春暉. 高光譜圖像處理技術(shù)[M]. 北京: 國(guó)防工業(yè)出版社, 2013: 11.

[7]鄧書斌. ENVI遙感圖像處理方法[M]. 北京: 科學(xué)出版社, 2010: 334.

[8]JI R R, GAO Y, HONG R C, et al. Spectral-spatial constraint hyperspectral image classification[J]. IEEE Transactions Geoscience and Remote Sensing, 2014, 52(3): 1811-1824.

[9]杜培軍, 譚琨, 夏俊士. 高光譜遙感影像分類與支持向量機(jī)應(yīng)用研究[M]. 北京: 科學(xué)出版社, 2012.

[10]張兵, 高連如. 高光譜圖像分類與目標(biāo)識(shí)別[M]. 北京: 科學(xué)出版社, 2011: 147-148.

[11]孫顯, 付琨, 王宏琦. 高分辨率遙感圖形理解[M]. 北京: 科學(xué)出版社, 2011: 145.

[12]NASRABADI N M. Regularization for spectral matched filter and RX anomaly detector[C]//SPIE, Defense and Security Symposium. [S.l.]: International Society for Optics and Photonics, 2008, 6966: 1-12.

[13]劉生, 王曉宇, 邱新法. 一種高分辨率遙感影像的數(shù)學(xué)形態(tài)學(xué)濾波算法[J]. 氣象與減災(zāi)研究, 2008, 31(4): 48-51.

[14]段汕. 形態(tài)學(xué)及其在遙感影像處理中的應(yīng)用研究[D]. 武漢: 武漢大學(xué), 2004: 115-124.

[15]PLAZA A, MARTINEZ P, PLAZA J, et al. Dimensionality reduction and classification of hyperspectral image data using sequences of extended morphological transformations[J]. IEEE Transactions on Geoscience and Remote Sensing, 2005, 43(3): 466-479.

[16]趙春暉. 數(shù)字形態(tài)濾波器理論及其算法研究[D]. 哈爾濱: 哈爾濱工業(yè)大學(xué), 1998: 59.

[17]王瀛, 梁楠, 郭雷. 一種基于修正擴(kuò)展形態(tài)學(xué)算子的高光譜遙感圖像端元提取算法[J]. 光子學(xué)報(bào), 2012, 41(6): 672-677.

[18]CHANG CHEIN-I, DU QIAN. Estimation of number of spectrally distinct signal sources in hyperspectral imagery[J]. IEEE Transactions on Geoscience and Remote Sensing, 2004, 42(3): 608-619.

(編輯徐巖)

Generalized morphology-based hyperspectral image unmixing algorithm

ZHAOYan1,ZHANGChunjing2,WANGDonghui3

(1.School of Electric & Control Engineering, Heilongjiang University of Science & Technology, Harbin 150022, China; 2.School of Electronic & Information Engineering, Heilongjiang University of Science & Technology, Harbin 150022, China; 3.College of Information & Communication Engineering, Harbin Engineering University, Harbin 150001, China)

This paper introduces an unmixing algorithm based on generalized morphology, thanks to comprehensive utilization of pixel spectral and spatial information. The algorithm study consists of introducing reference pixel and avoiding morphological operator which suffers from limitations in the pixel arrangement rule and replacement criteria; calculating the generalized open-closing operator using the modified energy function as a distance measure in response to cross substitution occurring in different types of junction in the image; and performing abundance inversion using fully constrainded least squares, after the endmember extraction. The experiment using simulated data and real data of mineral spectrum shows that proposed algorithm has a demonstrated ability to automatically decompose mixed pixel without prior information and exhibits a higher decomposition precision.

hyperspectral unmixing; generalized morphology; modified energy function

2015-03-23

黑龍江省教育廳科學(xué)技術(shù)研究面上項(xiàng)目(12541734)

趙巖(1976-),女,黑龍江省齊齊哈爾人,講師,碩士,研究方向:圖像處理、智能控制,E-mail:zh-ao-yan@sina.com。

10.3969/j.issn.2095-7262.2015.03.023

TN911.73

2095-7262(2015)03-0343-05

A