一種正交子空間投影高光譜圖像端元提取算法

趙　巖,王東輝，張春晶，黃奕程，于明華

(1. 黑龍江科技大學(xué) 電氣與控制工程學(xué)院，哈爾濱 150022；2. 哈爾濱工程大學(xué) 信息與通信工程學(xué)院，哈爾濱 150001；3. 齊齊哈爾車(chē)輛集團(tuán)，黑龍江 齊齊哈爾 161000)

?

一種正交子空間投影高光譜圖像端元提取算法

趙巖1,王東輝2，張春晶1，黃奕程3，于明華3

(1. 黑龍江科技大學(xué) 電氣與控制工程學(xué)院，哈爾濱 150022；2. 哈爾濱工程大學(xué) 信息與通信工程學(xué)院，哈爾濱 150001；3. 齊齊哈爾車(chē)輛集團(tuán)，黑龍江 齊齊哈爾 161000)

端元提取是高光譜圖像混合像元分解的關(guān)鍵問(wèn)題。針對(duì)正交子空間投影方法進(jìn)行端元提取需要端元先驗(yàn)知識(shí)的問(wèn)題，提出一種基于光譜最小信息熵的正交子空間投影高光譜圖像端元提取方法。以光譜最小信息熵判定最優(yōu)端元子集，同時(shí)將正交投影散度作為不同地物光譜向量之間相似性的測(cè)度指標(biāo)用于判別端元。利用模擬數(shù)據(jù)和真實(shí)數(shù)據(jù)對(duì)算法進(jìn)行驗(yàn)證，結(jié)果表明：算法不需先驗(yàn)知識(shí)，能夠自動(dòng)進(jìn)行端元提取，且精度較高。

高光譜圖像；正交子空間投影；端元提??；光譜最小信息熵；正交投影散度

高光譜圖像以像元為單位獲取地物信息，所記錄的是像元內(nèi)全部地物光譜的綜合。由于空間分辨率的約束和地物的復(fù)雜多樣性的影響，一個(gè)像元內(nèi)經(jīng)常包含有多種地物，稱為混合像元(Mixed Pixel)[1-2]，作為分解基礎(chǔ)的基本地物成分被稱為端元(Endmember)[3]。為了更準(zhǔn)確地描述真實(shí)地表的覆蓋情況，需要對(duì)混合像元進(jìn)行解混，計(jì)算出每一種地物成分在該像元中的所占比重。使用最廣泛的的線性光譜解混(Linear spectral unmixing)是利用線性光譜混合模型(Linear Spectral Mixing Model, LSMM)把圖像X中的每個(gè)混合像元分解成端元及其對(duì)應(yīng)豐度，進(jìn)而端元矩陣S和豐度矩陣A的過(guò)程[4]。端元提取的好壞是混合像元光譜分解效果的關(guān)鍵[6]。端元提取的主要算法有基于凸面幾何學(xué)的Pixel Purity Index (PPI)、N-FINDR、Vertex Component Analysis (VCA)[7]，基于數(shù)學(xué)形態(tài)學(xué)的Automated Morphological Endmember Extraction (AMEE)，可以同時(shí)進(jìn)行端元提取和豐度反演的Independent Component Correlation Algorithm (ICA)等。正交子空間投影算法(Orthogonal Subspace Projection, OSP)是在已知豐度分布的情況下根據(jù)信噪比最小推導(dǎo)出的，Harsanyi和Chang將OSP用于目標(biāo)檢測(cè)，通過(guò)信號(hào)光譜逐步提取出目標(biāo)信號(hào)，但需要提供端元的先驗(yàn)知識(shí)[8-9]。在OSP的基礎(chǔ)上，近年來(lái)研究者提出基于后驗(yàn)知識(shí)的信號(hào)空間正交投影(Signature Space orthogonal Projection, SSP)、目標(biāo)空間正交投影(Target Signature space orthogonal Projection, TSP)及斜子空間投影(Oblique Subspace Projection, OBSP)[10]，這些方法都屬于監(jiān)督分解。通常情況下獲取端元的先驗(yàn)知識(shí)是比較困難的。本文以光譜最小信息熵[11]判定最優(yōu)端元子集，并利用正交投影散度作為相似性測(cè)度來(lái)判定端元，提出一種非監(jiān)督的改進(jìn)正交子空間投影方法(Improved Orthogonal Subspace Projection, IOSP)，用以自動(dòng)提取端元。

1　正交子空間投影端元提取算法

設(shè)X為L(zhǎng)×1維的像元光譜向量，S=[s1,s2,…,sp]為L(zhǎng)×P維的端元光譜矩陣，A=[a1,a2,…,ap]T是P維向量，其中ai為對(duì)應(yīng)端元的豐度，N為L(zhǎng)×1維的隨機(jī)噪聲。像元的線性光譜混合模型表示為:

X=SA+N

(1)

(2)

(3)

OSP端元提取的步驟：

1) 根據(jù)凸面幾何學(xué)理論，通常認(rèn)為端元位于特征空間的角點(diǎn)上，即所謂最大光譜向量法。滿足式(4)，得到一個(gè)待選端元，它是圖像中亮度最大的像元。

(4)

2) 判斷待選端元是否為噪聲，如果是噪聲，除去；如果不是噪聲，則為端元。

3) 消除已經(jīng)提取端元在圖像中的影響，生成新的圖像數(shù)據(jù)。

OSP認(rèn)為像元是由目標(biāo)信號(hào)d和背景端元U兩部分組成。設(shè)d=si,1≤i≤P，則背景端元向量空間U=(s1…si-1si+1,…,sp)。因此，式(4)可以寫(xiě)為:

(5)

式中ai為si的豐度；b=(a1…ai-1ai+1…ap)T為對(duì)應(yīng)背景端元的豐度向量。

通過(guò)正交子空間投影算子U⊥對(duì)背景及噪聲進(jìn)行抑制，由最小二乘估計(jì)得到U⊥=I-UU#，其中，式中⊥為正交補(bǔ)空間，#為Moore-Penrose廣義偽逆矩陣，U#=(UTU)-1UT，-1為矩陣逆運(yùn)算，I為P×P維的單位陣。在假設(shè)噪聲與信號(hào)不相關(guān)的條件下，提取信號(hào)d的算子為：

(6)

其中，λ=(dTU⊥d)-1為一個(gè)表示分解精度的常數(shù)項(xiàng)[12]。

從式(6)可知，OSP算法在進(jìn)行端元提取時(shí)，需要圖像端元的先驗(yàn)知識(shí)。

4) 判斷是否滿足結(jié)束條件(已提取的端元數(shù)達(dá)到指定的數(shù)目)，若滿足，則結(jié)束計(jì)算；若不滿足，則循環(huán)步驟1)～3)。

2　改進(jìn)的正交子空間投影端元提取算法

2.1光譜最小信息熵

光譜最小信息熵(SMSE)算法利用高斯分布概率函數(shù)得到高光譜圖像像元光譜的概率，用來(lái)計(jì)算光譜的信息熵，之后通過(guò)設(shè)置閾值篩選出小于閾值的信息熵的光譜，并累加光譜曲線概率，經(jīng)過(guò)運(yùn)算獲得端元。 SMSE算法如下式：

(7)

其中，

式中xi為(j,k)像元在第i個(gè)波段的灰度值;μi和σi為第i個(gè)波段灰度值的均值和方差;P(j,k)為光譜曲線概率;a為常數(shù)。

2.2基于正交投影散度的端元判定

正交投影散度(orthogonal projection divergence, OPD)來(lái)源于OSP的原理[13]，其目的是最大程度地分離目標(biāo)信息和背景噪聲信息，計(jì)算得到的正交投影后的殘差，作為衡量不同地物光譜向量之間相似性大小的指標(biāo)。

考慮兩個(gè)M維光譜信號(hào)ri=[ri1,ri2,…,riM]T，rj=[rj1,rj2,…,rjM]T，則ri和rj之間的OPD為：

(8)

IOSP端元提取的步驟：

1) 按SMSE算法計(jì)算獲得待選端元。

2) 按基于正交投影散度的端元判定法判斷待選端元是否為噪聲，如果是噪聲，做出標(biāo)記，排除；如果不是噪聲，則為端元。

3) 基于正交投影散度的端元判定法獲得1個(gè)端元信號(hào)d=s1，利用d的正交補(bǔ)空間d⊥消除端元的d=s1影響，其中d⊥=I-d(dTd)-1dT。將d⊥作用于OSP模型中，d⊥X=d⊥Ub+d⊥N，噪聲有效地壓制為d⊥N。

將d⊥X作為新的圖像數(shù)據(jù)，提取端元e2，獲得d=[e1,e2]，循環(huán)運(yùn)算，直至提取P個(gè)端元。

3　仿真實(shí)驗(yàn)

3.1性能指標(biāo)

采用光譜角度距離(Spectral Angle Distance, SAD)作為指標(biāo)進(jìn)行精度評(píng)價(jià)。

對(duì)于第i個(gè)端元，SAD定義為：

(9)

圖1　Alunite, Buddingtonite, Chalcedony, Kaolinite的光譜Fig.1　Reflectances of Alunite, Buddingtonite, Chalcedony, Kaolinite

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

從USGS礦物光譜庫(kù)中選取4種相互之間線性獨(dú)立的端元光譜(Alunite, Buddingtonite, Chalcedony, Kaolinite)，見(jiàn)圖1。反射率的波長(zhǎng)為0.38～2.5 μm，按照Dirichlet分布進(jìn)行混合，先對(duì)端元豐度之和進(jìn)行歸一化，然后加上不同的高斯白噪聲，波段數(shù)是224，圖像大小是256×256像元，形成模擬實(shí)驗(yàn)數(shù)據(jù)。

抗噪性能實(shí)驗(yàn)：在混合光譜中加入隨機(jī)噪聲。

(10)

圖2　噪聲強(qiáng)度不同下算法性能的比較Fig.2　Comparison of Algorithm performance of different oise intensity

圖3　Nevada州Cuprite地區(qū)的AVIRIS圖像(193波段)Fig. 3　AVIRIS image data in cuprite region (193 band)

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

采用美國(guó)Cuprite Nevada地區(qū)的AVIRIS高光譜數(shù)據(jù)(圖像大小400×350像元，波長(zhǎng)為1.99～2.48 μm，于1995年成像)，172～221波段間的共50個(gè)波段數(shù)據(jù)，見(jiàn)圖3。去掉大氣吸收的221波段，實(shí)際采用49個(gè)波段進(jìn)行實(shí)驗(yàn)，圖像大小為13.3 MB。

IOSP、VCA、OSP算法所輸出Cuprite地區(qū)5種常見(jiàn)礦物端元與相應(yīng)實(shí)際地物參考光譜的相似度見(jiàn)表1。采用虛擬維度(VD)方法[14]獲得的圖像端元數(shù)目?？梢?jiàn)，IOSP考慮端元光譜相關(guān)性影響，進(jìn)一步提高了提取端元的純度，表現(xiàn)較為穩(wěn)定。

表1　3種算法的SAD比較

4　結(jié)　論

本文在研究OSP算法的基礎(chǔ)上，以光譜最小信息熵判定最優(yōu)端元子集，將正交投影散度用于不同地物光譜向量之間相似性的測(cè)度，利用正交投影散度均值判定候選端元是端元，還是噪聲，逐個(gè)獲得高光譜圖像中的所有端元。筆者提出的IOSP算法克服了OSP需要先驗(yàn)知識(shí)的不足，通過(guò)模擬、真實(shí)數(shù)據(jù)實(shí)驗(yàn)分析，表明IOSP算法能夠自動(dòng)進(jìn)行端元提取，且提取結(jié)果較好。

[1]Thouvenin P A, Dobigeon N, Tourneret J Y. Hyperspectral unmixing with spectral variability using a perturbed linear mixing model[J]. IEEE Transactions on Signal Processing, 2016, 64(2):525-538.

[2]Lin C H, Chi C Y, Wang Y H, et al. A fast hyperplane-based minimum-volume enclosing simplex algorithm for blind hyperspectral unmixing[J]. IEEE Transactions on Signal Processing, 2016, 64(8):1946-1961.

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

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

[5]劉國(guó)華.基于圖像處理的鉆頭測(cè)量系統(tǒng)[J].黑龍江大學(xué)工程學(xué)報(bào)，2016，7(1)：92-96.

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

[7]Chang C I, Chen S Y, Li H C, et al. Comparative study and analysis among ATGP, VCA, and SGA for finding endmembers in hyperspectral umagery[J]. IEEE Journal of Selected Topics in Applied Earth Observations & Remote Sensing, 2016:1-27.

[8]Harsanyi J C, Chang C I. Hyperspectral image classification and dimensionality reduction: an orthogonal subspace projection approach[J]. IEEE Transactions on Geoscience and Remote Sensing, 1994, 32(4): 779-785.

[9]Chang C I. Further results on relationship between spectral unmixing and subspace projection[J]. IEEE Transactions on Geoscience and Remote Sensing, 1998, 36(3): 1030-1032.

[10] Tu T M, Chen C H, Chang C I. A posteriori least squares orthogonal subspace projection approach to desired signature extraction and detection[J]. IEEE Transactions on Geoscience and Remote Sensing, 1997, 35(1):127-139.

[11] 楊可明, 劉士文, 王林偉,等. 光譜最小信息熵的高光譜影像端元提取算法[J]. 光譜學(xué)與光譜分析, 2014(8):2229-2233.

[12] Chang C I, Zhao X L, Mark L G, et al. Least squares subspace projection spproach to mixed pixel classification for hyperspectral images[J]. IEEE Transactions on Geoscience and Remote Sensing, 1998, 36(3): 898-912.

[13] Keshava N. Distance metrics and band selection in hyperspectral processing with applications to material identification and spectral libraries[J]. IEEE Trans. Geoscience Remote Sensing, 2004, 42(7): 1552-1565.

[14] Chang C I, Du Q. Estimation of number of spectrally distinct signal sources in hyperspectral imagery[J]. IEEE Transactions on Geoscienee and Remote Sensing, 2004, 42(3): 608-619.

An Orthogonal Subspace Projection for endmember extraction algorithm in hyperspectral images

ZHAO Yan1，WANG Dong-Hui2，ZHANG Chun-Jing1，HUANG Yi-Cheng3，YU Ming-Hua3

(1.SchoolofElectricalandControlEngineering,HeilongjiangUniversityofScienceandTechnology,Harbin150022,China;2.CollegeofInformationandCommunicationEngineering,HarbinEngineeringUniversity,Harbin150001,China;3.QiqiharVehicleGroup,Qiqihar161000,Heilongjiang,China)

Endmember extraction is the key problem of hyperspectral images mixed pixel decomposition. To address Orthogonal Subspace Projection method for endmember extraction needing to solve the problem of endmember prior knowledge, a hyperspectral image endmember extraction method of Orthogonal Subspace Projection based on spectral minimum shannon entropy was proposed. It determined the optimal endmembers subset using the criteria of spectral minimum shannon entropy, and the Orthogonal Projection Divergence as a measure index of similarity between different spectral vector for judging the endmember. The algorithm was verified by the simulated hyperspectral images and real hyperspectral images. Experimental results show that the proposed algorithm can extract endmember automatically without prior knowledge and achieve relatively high extraction accuracy.

hyperspectral images; Orthogonal Subspace Projection(OSP); endmember extraction; Spectral Minimum Shannon Entropy(SMSE); Orthogonal Projection Divergence(OPD)

10.13524/j.2095-008x.2016.03.046

2016-07-30

黑龍江省教育廳科學(xué)技術(shù)研究項(xiàng)目(12541734)；國(guó)家自然科學(xué)基金資助項(xiàng)目(61275010)

趙巖(1976-)，女，黑龍江齊齊哈爾人，講師，博士研究生，研究方向：圖像處理、智能控制，E-mail:zh_ao_yan@sina.com。

TP751

A

2095-008X(2016)03-0082-05