基于光譜相似度量的高光譜圖像異常檢測(cè)算法

王玉磊，趙春暉，齊 濱

(哈爾濱工程大學(xué)信息與通信工程學(xué)院，哈爾濱150001)

高光譜遙感圖像具有很高的光譜分辨率，是一種圖譜合一的新型遙感數(shù)據(jù)。借助其豐富的光譜信息，可以反映目標(biāo)間的細(xì)微差異，使人們能夠發(fā)現(xiàn)許多在常規(guī)遙感中無(wú)法或難以探測(cè)的地面目標(biāo)，這為目標(biāo)檢測(cè)提供了有利的支持，因而越來(lái)越受到人們的重視。異常檢測(cè)是高光譜圖像處理的一個(gè)重要分支，它能夠在沒(méi)有先驗(yàn)光譜信息的情況下檢測(cè)到與周圍環(huán)境存在光譜差異的目標(biāo)，在軍事防御、農(nóng)業(yè)、地質(zhì)探測(cè)、醫(yī)療診斷以及環(huán)境監(jiān)控等方面都有重要的應(yīng)用，成為了高光譜目標(biāo)檢測(cè)領(lǐng)域的一個(gè)研究熱點(diǎn)。從已有文獻(xiàn)中可知，異常目標(biāo)具有以下特性:(1)它是unexpected;(2)與背景相比，它的樣本空間非常小;(3)它的發(fā)生概率非常低;(4)異常目標(biāo)與周圍背景的光譜信息有很大的區(qū)別。

文獻(xiàn)［1］提出了一種基于廣義似然比檢測(cè)(GLRT)的異常算子，并稱之為RX異常算子。目前RX算子及其改進(jìn)模型已成功應(yīng)用于高光譜圖像異常探測(cè)［2－6］，并已成為高光譜異常檢測(cè)的一個(gè)基準(zhǔn)算法。該算法是一種局部異常檢測(cè)算法，用于在高斯背景統(tǒng)計(jì)特性和空間白化的條件下檢測(cè)空間模式己知而光譜特性未知的目標(biāo)物，其漸進(jìn)意義下的檢測(cè)器將原始算法簡(jiǎn)化為求待檢測(cè)點(diǎn)到周圍背景均值的馬氏距離。RX算法是在一些簡(jiǎn)化的假設(shè)條件下構(gòu)造的似然比檢測(cè)算子，直接利用RX算法對(duì)高光譜圖像進(jìn)行處理將會(huì)產(chǎn)生較高的虛警概率。這主要是由兩方面原因造成的: (1)RX算法中所采用的局部統(tǒng)計(jì)模型假定數(shù)據(jù)是空間不相關(guān)的，或者是空間白化的，且數(shù)據(jù)要求服從局部正態(tài)分布，然而，由于實(shí)際地物分布是復(fù)雜多變的，這種假設(shè)不能全面地描述真實(shí)場(chǎng)景的情況;(2)當(dāng)RX算法直接用于高光譜圖像處理時(shí)，需要計(jì)算樣本的協(xié)方差矩陣及其逆矩陣，而協(xié)方差矩陣的維數(shù)會(huì)隨著波段數(shù)目的增加而迅速增加，這將會(huì)帶來(lái)巨大的計(jì)算量。

近年來(lái)，隨著統(tǒng)計(jì)學(xué)習(xí)理論的快速發(fā)展和逐步完善，應(yīng)用核機(jī)器學(xué)習(xí)算法對(duì)高光譜圖像進(jìn)行處理成為可能。將高光譜數(shù)據(jù)映射到高維的特征空間進(jìn)行處理，充分挖掘其隱含的非線性信息，進(jìn)一步提高高光譜小目標(biāo)與背景之間的分離性能。目前已有學(xué)者Kwon對(duì)傳統(tǒng)高光譜小目標(biāo)檢測(cè)算法，RX算子進(jìn)行了改進(jìn)［2］，利用核函數(shù)將原始空間光譜信號(hào)映射到高維特征空間，進(jìn)而實(shí)現(xiàn)了小目標(biāo)的異常檢測(cè);通過(guò)與傳統(tǒng)高光譜小目標(biāo)檢測(cè)算法檢測(cè)性能的比較，得出核算法在提高目標(biāo)檢測(cè)概率和降低虛警方面有顯著提高，然而該算法仍然是基于RX的基本假設(shè)基礎(chǔ)的前提下，并且仍然需要巨大的計(jì)算量。

為了去除對(duì)背景模型的假設(shè)以及提高異常檢測(cè)算法的實(shí)時(shí)性，本文提出了一種新型無(wú)假設(shè)模型的快速異常檢測(cè)算法。首先討論了光譜相似度量的方法及其效果，選用光譜角余弦作為光譜相似度量準(zhǔn)則;然后從端元提取入手，采用迭代誤差分析方法(IEA)對(duì)高光譜數(shù)據(jù)進(jìn)行異常端元提取;最后采用核函數(shù)方法，利用得到的異常端元作為參考數(shù)據(jù)，以核光譜角余弦(KSAC)作為相似度量準(zhǔn)則，對(duì)所有像素點(diǎn)進(jìn)行判別，得到異常檢測(cè)結(jié)果。

1 光譜相似度量

很多光譜圖像處理算法都是基于光譜相似或者光譜距離的，從而估計(jì)兩個(gè)光譜向量的相似程度。

歐氏距離是傳統(tǒng)的基于幾何學(xué)的向量距離，其計(jì)算公式如下

歐氏距離是向量距離最簡(jiǎn)單的計(jì)算方法，但是卻不是精確的方法。

圖1和圖2分別給出了兩組模擬數(shù)據(jù)。圖1給出了兩個(gè)弱相關(guān)的向量(r2≈0.000 6)，而圖2則給出了兩個(gè)強(qiáng)正相關(guān)的向量(r2≈0.999 3)。

通過(guò)式(1)計(jì)算兩組數(shù)據(jù)的歐氏距離均為0.481 18。若以歐氏距離作為相似度標(biāo)準(zhǔn)，兩組數(shù)據(jù)的歐氏距離相同，應(yīng)得到相同的判別結(jié)果。然而，從圖1和圖2中很明顯可以看出，第一組數(shù)據(jù)是弱相關(guān)甚至不相關(guān)的，而第二組數(shù)據(jù)具有強(qiáng)相關(guān)性。因此歐氏距離并不能很好地反應(yīng)向量的相似程度。

圖1 兩個(gè)弱相關(guān)向量(歐氏距離為0.481 18)Fig.1 Two weakly correlated vectors

圖2 兩個(gè)強(qiáng)正相關(guān)向量(歐氏距離為0.481 18)Fig.2 Two strong positive correlated vectors

光譜角余弦(Spectral Angel Cosine，SAC)是一種計(jì)算光譜相似度的方法［7］，它同時(shí)反映了光譜向量在數(shù)值上的差異和光譜曲線之間、形狀之間的差異。相比馬氏距離，用光譜角余弦作為光譜相似度測(cè)量具有明顯的優(yōu)點(diǎn)。光譜角余弦的計(jì)算公式如下

式中:x和y為像元的光譜向量;〈·，·〉為內(nèi)積運(yùn)算。光譜角作為兩個(gè)光譜向量的夾角，光譜角越小，兩個(gè)光譜向量的相似度就越高，SAC就越接近1。

同樣以圖1和圖2兩組模擬數(shù)據(jù)為例，以光譜角度作為光譜相似度量的標(biāo)準(zhǔn)，計(jì)算兩組數(shù)據(jù)的相似程度。圖1通過(guò)式(2)計(jì)算得兩個(gè)光譜向量的光譜角為32.662 6°;而圖2通過(guò)式(2)計(jì)算得兩個(gè)光譜向量的光譜角為3.799 2°。兩組數(shù)據(jù)得到的光譜角度差異較大，通過(guò)設(shè)定適當(dāng)?shù)慕嵌乳撝担瑘D1兩個(gè)光譜向量可認(rèn)為是不相似的，而圖2兩個(gè)光譜向量可認(rèn)為是相似的，用此標(biāo)準(zhǔn)得到的結(jié)果是比較符合實(shí)際情況的。

2 基于光譜相似度量的異常檢測(cè)算法

2.1 端元提取算法

光譜相似度量方法需要有特定的參考光譜向量，在此可用端元提取的方法來(lái)獲得。

高光譜數(shù)據(jù)端元提取是理解高光譜數(shù)據(jù)，進(jìn)而對(duì)高光譜數(shù)據(jù)進(jìn)行下一步處理(如混合像元分解、分類、地物識(shí)別)的重要前提。目前，高光譜數(shù)據(jù)端元提取的主要方法有純像元指數(shù)法(Pure Pixel Index，PPI)、N－FINDR、迭代誤差分析法(Iterative Error Analysis，IEA)、頂點(diǎn)成分分析(Vertex Component Analysis，VCA)、最大距離法單行體體積法等［8－10］。在此選用了迭代誤差分析法對(duì)高光譜數(shù)據(jù)進(jìn)行端元提取。

迭代誤差分析(IEA)是一種不需要對(duì)原始數(shù)據(jù)降維或去冗余而直接對(duì)數(shù)據(jù)進(jìn)行處理的端元提取方法。迭代誤差分析是建立在線性混合模型的基礎(chǔ)上，線性混合模型可表示為:

在式(1)中，假設(shè)端元矩陣S已知，則x的最小均方估計(jì)為:

均方根誤差表示為模長(zhǎng):

迭代誤差分析算法的步驟如下:

首先，給定一個(gè)初始向量(一般為圖像中所有光譜的均值向量)，用它來(lái)代替式(2)中的S;

第二，用得到的S計(jì)算均方根誤差，對(duì)圖像數(shù)據(jù)中每個(gè)向量通過(guò)式(3)求取均方根誤差，得到一幅誤差圖像;

第三，選擇誤差圖像中較大的點(diǎn)，并以此點(diǎn)為中心，以一定角度為門限，選擇與此點(diǎn)角度較小的點(diǎn)進(jìn)行平均，將得到的平均波譜作為第一個(gè)端元。

重復(fù)上面步驟，直至找到規(guī)定數(shù)量的端元或者均方根誤差小于一個(gè)特定門限。

2.2 異常檢測(cè)算法

在采用IEA算法獲取異常端元后，就可用異常端元作為參考端元，對(duì)所有像元進(jìn)行相似度量，從而實(shí)現(xiàn)異常檢測(cè)。相似度量準(zhǔn)則采用核光譜角余弦準(zhǔn)則［11］。

光譜角余弦如式(2)所示，將原始輸入數(shù)據(jù)空間通過(guò)非線性映射函數(shù)Φ映射到高維特征空間中，此時(shí)式(2)變?yōu)?/p>

由于特征空間的維數(shù)很高(甚至可能趨于無(wú)窮)，同時(shí)由于非線性映射函數(shù)未知，這使得直接在特征空間進(jìn)行計(jì)算變得很困難甚至不可能實(shí)現(xiàn)。核函數(shù)方法可以將特征空間的點(diǎn)積轉(zhuǎn)化為輸入空間數(shù)據(jù)的核函數(shù)實(shí)現(xiàn)，從而巧妙地避開(kāi)了在高維的特征空間計(jì)算點(diǎn)積。

核函數(shù)k(x，y)表示如下:

常用的核函數(shù)有高斯徑向基核函數(shù)和多項(xiàng)式核函數(shù)，數(shù)學(xué)模型分別如下:

高斯徑向基核函數(shù):

多項(xiàng)式核函數(shù):

在此選用高斯徑向基核函數(shù)(RBF)，這是因?yàn)楫?dāng)選用高斯徑向基核函數(shù)時(shí)，若x=y，則有k(x，y)=1，因而式(6)的分母部分為1，可大大降低式(6)的復(fù)雜度。

將式(7)和式(8)代入式(6)并化簡(jiǎn)得KSAC的數(shù)學(xué)表達(dá)式為:

式中:x是待檢測(cè)像元的光譜向量;y是參考像元的光譜向量，即端元光譜向量。

根據(jù)式(10)，將IEA算法獲取的異常端元代入y，x為待檢測(cè)像元，設(shè)定合適的閾值即能得到異常檢測(cè)結(jié)果。算法流程如3所示。

圖3 算法流程Fig.3 Algorithm flowchart

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

為了驗(yàn)證算法的有效性，利用真實(shí)的高光譜圖像進(jìn)行仿真實(shí)驗(yàn)，實(shí)驗(yàn)中采用的高光譜圖像數(shù)據(jù)是由AVIRIS傳感器獲取的美國(guó)圣地亞哥機(jī)場(chǎng)圖像，如圖4所示。該數(shù)據(jù)已經(jīng)過(guò)大氣校正和幾何校正等預(yù)處理，是以反射率數(shù)據(jù)形式進(jìn)行存儲(chǔ)的，關(guān)于該高光譜圖像數(shù)據(jù)的具體的參數(shù)如表1所示。

從原始圖像中截取60×60的一段含有四個(gè)飛機(jī)目標(biāo)的子圖，圖5為截取的高光譜子圖的第一個(gè)波段示意圖及地面目標(biāo)分布。本文分別以光譜角余弦(IEA_SAC)與核光譜角余弦(IEA_ KSAC)作為光譜相似度量準(zhǔn)則，設(shè)計(jì)了高光譜異常檢測(cè)算法。

圖4 高光譜圖像數(shù)據(jù)源Fig.4 Hyperspectral data

表1 實(shí)驗(yàn)用高光譜圖像數(shù)據(jù)源參數(shù)Table 1 The parameters of hyperspectral imagery for numerical experiment

圖5 第1波段及地面目標(biāo)分布圖Fig.5 The 1st band and real target distribution

4.1 端元提取

實(shí)驗(yàn)中，在對(duì)高光譜數(shù)據(jù)進(jìn)行歸一化處理后，用IEA算法對(duì)高光譜數(shù)據(jù)進(jìn)行端元提取。其中，IEA算法有兩個(gè)重要的參數(shù)，即參與平均的波譜數(shù)目R和角度門限θ。參與平均的波譜數(shù)目R決定了所提取出的端元是背景端元還是異常端元，這是因?yàn)楸尘皵?shù)目要遠(yuǎn)遠(yuǎn)多于異常點(diǎn)的數(shù)目。因此，在提取出的端元中，如果參與平均的波譜數(shù)目較多，一般為背景端元;反之，參與平均的波譜數(shù)目較少，應(yīng)為異常端元。角度門限θ決定了光譜是否參與平均，一般來(lái)說(shuō)，角度門限θ較小。

為了驗(yàn)證參數(shù)范圍對(duì)算法魯棒性的影響，我們討論參與平均的波譜數(shù)目R和角度門限θ對(duì)端元提取的影響。實(shí)驗(yàn)中對(duì)IEA算法分別設(shè)置不同的參數(shù)進(jìn)行端元提取，得到的結(jié)果如表2和表3所示。

表2 θ=8°時(shí)IEA算法中參與平均的波譜數(shù)目Table2 The average spectrum number in IEA at θ=8°

表3 R=150時(shí)IEA算法中參與平均的波譜數(shù)目Table3 The average spectrum number in IEA at R=150

通過(guò)表2可以看出，參與平均的像素?cái)?shù)R對(duì)IEA算法沒(méi)有太大影響。但是當(dāng)R取值很小，各個(gè)端元參與平均的波譜數(shù)目比較接近，此時(shí)無(wú)法判斷該端元是異常端元還是背景端元;當(dāng)R到達(dá)一定值后，此參數(shù)對(duì)端元提取的結(jié)果就沒(méi)有太大的影響。此時(shí)，為了節(jié)省計(jì)算量，應(yīng)盡量選取較小的R值，故需根據(jù)與圖像總像素?cái)?shù)的比例選取合適的R值。由表3分析可得，當(dāng)角度門限θ為2°～5°，由于角度過(guò)小，得到的結(jié)論并不完全準(zhǔn)確。隨著角度門限的增大，參與端元1平均的像素?cái)?shù)始終較少，但是光譜之間的差異會(huì)隨之增大，此時(shí)不能保證參與平均的像素均屬于同一端元。因此，角度門限應(yīng)盡量取較小值，但并不是越小越好。

本文中設(shè)置參數(shù)R=400，θ=8°，在IEA算法的每次循環(huán)中，由R和θ所確定的參與平均波譜的數(shù)目如表4所示。

表4 IEA算法中參與平均的波譜數(shù)目Table 4 The average spectrum number in IEA algorithm

從表4中可以看出，參與端元1平均的波譜數(shù)目只有36個(gè)，說(shuō)明此端元為異常，即目標(biāo)端元。

4.2 異常檢測(cè)

將得到的目標(biāo)端元作為參考端元，可以進(jìn)行相似度量從而完成異常檢測(cè)。首先采用式(2)的光譜角余弦表達(dá)式作為相似度量判別式，異常檢測(cè)結(jié)果如圖6所示，圖6(a)為原始高光譜圖像第一波段，圖6(b)為IEA_SAC算法檢測(cè)結(jié)果二值化圖。

圖6 IEA_SAC異常檢測(cè)結(jié)果Fig.6 IEA_SAC anomaly detection result

采用式(10)的核光譜角余弦表達(dá)式作為相似度量判別式。檢測(cè)結(jié)果如圖7所示。其中圖7 (a)為原始高光譜圖像第一波段，圖7(b)為IEA_ KSAC算法檢測(cè)結(jié)果二值化圖。

圖7 IEA_KSAC異常檢測(cè)結(jié)果Fig.7 IEA_KSAC anomaly detection result

IEA_SAC算法與IEA_KSAC算法的閾值取值區(qū)間如表5所示。閾值取值區(qū)間是指得到余弦最大值與最小值之間的差異。從表5中可以看出，KSAC的取值區(qū)間明顯大于SAC取值區(qū)間，可明顯降低算法對(duì)閾值選擇的敏感度。同時(shí)，由于引入高斯徑向基核函數(shù)對(duì)KSAC進(jìn)行化簡(jiǎn)，從而無(wú)需多次計(jì)算內(nèi)積，使式(10)的計(jì)算復(fù)雜度低于式(5)。從表2算法運(yùn)行時(shí)間可以看出，KSAC算法的運(yùn)行時(shí)間還不到SAC運(yùn)行時(shí)間的一半，當(dāng)數(shù)據(jù)量較大時(shí)，KSAC的效率要明顯高于SAC。

表5 閾值取值區(qū)間與運(yùn)行時(shí)間Table5 Threshold value interval and running time

從檢測(cè)結(jié)果二值化圖(見(jiàn)圖7和圖8)中可以看出，基于核光譜角余弦的異常檢測(cè)算法的檢測(cè)效果比較理想，飛機(jī)目標(biāo)被很好地檢測(cè)出來(lái)。同時(shí)，由于該算法沒(méi)有大量的矩陣運(yùn)算及內(nèi)積運(yùn)算，因此算法復(fù)雜度較低，運(yùn)算速度快。

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

本文提出了一種新型的高光譜異常檢測(cè)算法。該算法利用光譜相似度量的方法，通過(guò)計(jì)算當(dāng)前像素點(diǎn)與參考異常端元的相似性進(jìn)行異常判斷。文中給出了基于高斯徑向基的核光譜角余弦(KSAC)的定義及表達(dá)式，并設(shè)計(jì)了基于KSAC的高光譜圖像異常檢測(cè)算法。該方法適用于實(shí)際的高光譜圖像異常檢測(cè)，并具有檢測(cè)效果好、算法復(fù)雜度低、運(yùn)算速度快的特點(diǎn)。然而，核函數(shù)中參數(shù)優(yōu)化是一個(gè)重要的研究方向，文中并沒(méi)有給出核函數(shù)及參數(shù)的優(yōu)化，因此在參數(shù)選擇上徑向基核函數(shù)的寬度是靠大量仿真實(shí)驗(yàn)比較得出的，沒(méi)有實(shí)現(xiàn)自動(dòng)選取，該項(xiàng)工作有待進(jìn)一步進(jìn)行研究。

［1］Irving S Reed，Yu Xiao－li.Adaptive multiple－band cfar detection of an optical pattern with unknown spectral distribution［J］.EEE Transactions on Acoustics Speech and Signal Processing，1990，38(10):1760－1769.

［2］Heesung K，Nasrabadi N M.Kernel RX－algorithm:A nonlinear anomaly detector for hyperspectral imagery［J］.IEEE Transactions on Geoscience and Remote Sensing，2005，43(2):388－397.

［3］梅鋒，趙春暉.基于空域?yàn)V波的核RX高光譜圖像異常檢測(cè)算法［J］.哈爾濱工程大學(xué)學(xué)報(bào)，2009，30 (6):697－702.

Mei Feng，Zhao Chun－h(huán)ui.Spatial filter based anomaly detection algorithm for hyperspectral imagery kernel RX detectors［J］.Journal of Harbin Engineering University，2009，30(6):697－702.

［4］Zhao Chun－h(huán)ui，Wang Yu－lei，Mei Feng.Kernel ICA feature extraction for anomaly detection in hyperspectral imagery［J］.Chinese Journal of Electronics，2012，21 (2):265－269.

［5］趙春暉，胡春梅.基于目標(biāo)正交子空間投影加權(quán)的高光譜圖像異常檢測(cè)算法［J］.吉林大學(xué)學(xué)報(bào):工學(xué)版，2011，41(5):1468－1474.

Zhao Chun－h(huán)ui，Hu Chun－mei.Weighted anomaly detection algorithm for hyperspectral image based on target orthogonal subspace projection［J］.Journal of Jilin U－niversity(Engineering and Technology Edition)，2011，41(5):1468－1474.

［6］王斯博，楊春玲.基于RX算法的高光譜紅外弱小目標(biāo)檢測(cè)［J］.紅外技術(shù)，2010，32(4)::204－208.

Wang Si－bo，Yang Chun－ling.Detection of dim small targetbased oil rx algorithm for hyperspectral ir imagery［J］.Infrared Technology，2010，32(4):204－208.

［7］James Norman Sweet.The spectral similarity scale and its application to the classification of hyperspectral remote sensing data［D］.New York:State university of New York，College of Environmental Science and Forestry，2002.

［8］Antonio Plaza，Pablo Martínez，Rosa Pérez，et al.A quantitative and comparative analysis of endmember extraction algorithms from hyperspectral data［J］.IEEE Transactions on Geoscience and Remote Sensing，2004，42(3):650－663.

［9］Sun Li－xin，Zhang Ying，Bert Guindon.Improved iterative error analysis for endmember extraction from hyperspectral imagery［J］.Processing of SPIE，2008，7086: 11－18.

［10］李智勇，郁文賢，趙和鵬.迭代誤差分析方法在高光譜異常檢測(cè)中的應(yīng)用［J］.系統(tǒng)工程與電子技術(shù)，2008，30(12):2340－2344.

Li Zhi－yong，Yu Wen－xian，Zhao He－peng.Application of iterative error analysis in hyperspectral anomaly detection［J］.Systems Engineering and Electronics，2008，30(12):2340－2344.

［11］諶德榮，孫波，陶鵬，等.基于核光譜角余弦的高光譜圖像空間鄰域聚類方法［J］.電子學(xué)報(bào)，2008，36 (10):1992－1995.

Chen De－rong，Sun Bo，Tao Peng，et al.Spatial neighboring clustering method for hyperspectral imagery based on kernel spectral angel cosine［J］.Acta Electronica Sinia，2008，36(10):1992－1995.