采用PPI算法改進(jìn)的一種數(shù)學(xué)形態(tài)學(xué)端元提取方法

徐 君,王彩玲,王 麗

1.西安航空學(xué)院電子工程學(xué)院，陜西 西安 710077；2.西安石油大學(xué)計算機學(xué)院，陜西 西安 710065

現(xiàn)有的端元提取算法大多只是利用高光譜影像的光譜信息，把高光譜數(shù)據(jù)當(dāng)作空間無序的光譜矢量來處理，對圖像的空間信息利用不足。實際上，地物在圖像空間上的分布是有規(guī)律的，是具有一定的幾何形狀的，充分利用圖像的空間信息有益于提高端元提取的效率和精度[1-4]。文獻(xiàn)[5]提出的自動形態(tài)學(xué)端元提取(automated morphological endmember extraction，AMEE)算法同時利用了圖像的光譜信息和空間信息，在試驗中取得了不錯的效果。

數(shù)學(xué)形態(tài)學(xué)的兩個基本操作是膨脹和腐蝕，灰度圖像的膨脹和腐蝕操作是在結(jié)構(gòu)元素內(nèi)尋找像元的灰度最大值和最小值來取代目標(biāo)像元的灰度值[6]。而高光譜圖像是一個數(shù)據(jù)立方體，采用尋找灰度最大值和最小值進(jìn)行膨脹和腐蝕操作顯然是不合理的。AMEE算法對高光譜圖像的膨脹和腐蝕操作進(jìn)行了拓展性定義，高光譜圖像的膨脹和腐蝕不再是尋找像元的灰度最大值和最小值，而是尋找結(jié)構(gòu)元素內(nèi)的最純像元和混合度最大像元。AMEE算法通過計算和比較結(jié)構(gòu)元素內(nèi)像元之間或者像元與凸面單形體中心之間的光譜角距離來定量化衡量像元的相對純度，它將最純像元與混合程度最大像元的光譜角距離定義為離心率指數(shù)(morphological eccentricity index，MEI)，并賦值給結(jié)構(gòu)元素內(nèi)的最純像元，擴大結(jié)構(gòu)元素并不斷迭代更新，最后得到一幅MEI圖像，根據(jù)MEI值由高到低在像元中選取端元。

然而AMEE算法也存在一些缺陷。首先，AMEE算法將結(jié)構(gòu)元素內(nèi)所有像元的均值光譜作為混合程度最大的像元，這在結(jié)構(gòu)元素內(nèi)混合像元占大多數(shù)時是合理的，但結(jié)構(gòu)元素內(nèi)的純凈像元較多時，像元的均值光譜將更接近純像元，此時距離均值光譜較遠(yuǎn)的是混合程度較大的像元，而不是純像元，此時的膨脹和腐蝕操作就是不正確的。其次，AMEE算法將結(jié)構(gòu)元素內(nèi)最純像元與混合程度最大像元之間的光譜角距離定義為MEI來表示像元純度，但不同結(jié)構(gòu)元素內(nèi)混合程度最大的像元也是不一樣的，因此計算MEI值時的參考標(biāo)準(zhǔn)是不同的，MEI值的大小也就無法真實有效地表示像元的純度[7]。

另外，從現(xiàn)有的純像元指數(shù)(pure pixel index，PPI)端元提取算法的角度來看，PPI算法將特征空間中所有的樣本點在隨機生成的直線上進(jìn)行投影，然后統(tǒng)計每個像元投影落在直線兩端的次數(shù)來選取端元[8-11]。這種方法只是把高光譜圖像看成是無序光譜向量的集合，僅僅利用了圖像的光譜信息，而忽略了圖像中地物分布的空間結(jié)構(gòu)信息。

受AMEE算法啟發(fā)，如果能在結(jié)構(gòu)元素內(nèi)運用PPI算法尋找最純凈像元，并給像元賦值相應(yīng)的純凈指數(shù)，拓展定義形態(tài)學(xué)中的膨脹操作和腐蝕操作對整幅圖像進(jìn)行處理，擴大結(jié)構(gòu)元素不斷迭代運算并累計更新像元的PPI值，最終可得到一幅PPI值圖像，那么端元就可以在PPI累計值較大的像元中選取。這種方法可以彌補上述AMEE算法存在的不足，改善AMEE算法的性能。

1 現(xiàn)有的AMEE算法與PPI算法

1.1 AMEE算法

AMEE算法將傳統(tǒng)數(shù)學(xué)形態(tài)學(xué)理論進(jìn)行了拓展性定義，給出了高光譜圖像的形態(tài)學(xué)膨脹和腐蝕操作方法。AMEE算法中定義了兩種腐蝕和膨脹操作[5]。

(1)將結(jié)構(gòu)元素內(nèi)某像元到其他像元的距離累計之和作為純度指標(biāo)。結(jié)構(gòu)元素K內(nèi)某個像元f(x,y)到其他像元的距離累計之和用式(1)表示

?(s,t)∈K

(1)

此時，膨脹和腐蝕操作的拓展定義如下

d(x,y)=(f?K)(x,y)=

(2)

e(x,y)=(f⊕K)(x,y)=

(3)

式中，d(x,y)、e(x,y)是結(jié)構(gòu)元素K內(nèi)的像元，d(x,y)與K內(nèi)其他像元的距離累計之和最大，是K中最純的像元；e(x,y)與K內(nèi)其他像元的距離累計之和最小，是K中混合度最大的像元。

(2)將結(jié)構(gòu)元素內(nèi)像元到結(jié)構(gòu)元素中心的距離作為純度指標(biāo)。假設(shè)結(jié)構(gòu)元素K像元數(shù)量為M，那么K的中心可定義為

(4)

結(jié)構(gòu)元素K內(nèi)像元f(x,y)到K的中心的距離可表示為

D′(f(x,y),K)=dist(f(x,y),cK)

(5)

此時膨脹和腐蝕操作可分別拓展定義為

d′(x,y)=(f?K)(x,y)=

(6)

e′(x,y)=(f⊕K)(x,y)=

(7)

式中，d′(x,y)是距離K的中心最遠(yuǎn)的像元，即相對最純的像元；e′(x,y)是距離K的中心最近的像元，即混合度最大的像元。

AMEE算法將d(x,y)與e(x,y)之間的光譜夾角距離定義為MEI，MEI用來定量表示像元的相對純度，并賦值給結(jié)構(gòu)元素中最純的像元，計算公式如式(8)所示

MEI(n,m)=dist(d(x,y),e(x,y))

(8)

AMEE算法可以描述如下：從最小結(jié)構(gòu)元素Kmin開始，在整幅影像中移動，每移動一次，都會在結(jié)構(gòu)元素內(nèi)找到一個最純凈像元和一個混合度最大的像元，計算MEI值并將其賦給最純凈像元。然后逐漸增大結(jié)構(gòu)元素的尺寸，并重復(fù)執(zhí)行上述操作，達(dá)到預(yù)設(shè)的最大迭代次數(shù)Imax后終止，最終獲得一幅MEI圖像，參考MEI的值由高到低選取端元。

1.2 PPI算法

PPI算法的思想是基于凸面幾何學(xué)理論。凸面幾何學(xué)理論認(rèn)為高光譜影像的所有像元在高維光譜特征空間中對應(yīng)的樣本點呈散點圖分布，所有樣本點包含在一個凸面單形體內(nèi)部，那些純凈像元(端元)位于凸面單形體的頂點處[12-17]。

凸面單形體的特點決定了位于單形體頂點的樣本點在特征空間中任意直線上的投影必直線的兩端。以二維空間為例，如圖1所示，端元e1、e2、e3位于三角形的頂點，混合像元是由端元混合而成，都位于三角形的內(nèi)部。在這個二維空間中隨機生成一條直線，將所有樣本點投影到直線上，三角形的頂點必將投影到線段的兩端，三角形內(nèi)部的點會投影到線段的內(nèi)部[18]。

圖1 PPI算法中像元投影原理Fig.1 Schematic diagram of pixel projection principle in PPI algorithm

在特征空間中生成大量的隨機直線進(jìn)行投影，統(tǒng)計每個像元落在線段兩個端點的次數(shù)作為像元純凈指數(shù)，用來定量表示該像元的純凈度。顯然，某個像元對應(yīng)的像元純凈指數(shù)越大，它是純像元的可能性就越高[19-22]。

PPI算法雖然原理比較簡單直觀，但有PPI算法沒有利用圖像的空間信息，需要利用可視化工具手動選擇端元，自動化程度不高，而且投影直線都是隨機生成的，導(dǎo)致每次算法運行所提取的端元并不一定相同。

2 PPI-AMEE算法

膨脹和腐蝕操作的目的是找出結(jié)構(gòu)元素內(nèi)最純像元和混合度最大的像元。AMEE算法通過計算結(jié)構(gòu)元素內(nèi)的最純像元與混合度最大像元之間的光譜角距離作為MEI表示像元的純度。而PPI算法利用端元位于凸面單形體頂點的特性，將所有樣本點投影到在特征空間中隨機生成的直線上，統(tǒng)計樣本點投影在直線兩端的次數(shù)作為純像元指數(shù)來表示像元的純度。鑒于此，如果在結(jié)構(gòu)元素內(nèi)部用PPI代替MEI來尋找最純凈像元和混合程度最嚴(yán)重的像元，進(jìn)而拓展性地定義形態(tài)學(xué)中的膨脹和腐蝕操作，也能達(dá)到同樣甚至更優(yōu)的端元提取效果。本文結(jié)合PPI算法和AMEE算法提出一種PPI-AMEE端元提取方法。

以二維空間為例，圖2所示為像元在特征空間中的點云分布，結(jié)構(gòu)元素覆蓋范圍內(nèi)的像元只是其中的一小部分。圖2中三角形框abc和aef分別代表不同的結(jié)構(gòu)元素所覆蓋的像元。根據(jù)AMEE方法，需要在結(jié)構(gòu)元素范圍內(nèi)尋找距離中心最遠(yuǎn)的像元作為最純像元，然而觀察三角形框abc和aef，距離中心最遠(yuǎn)的像元分別是偽端元c和e，并不是實際的端元a，因此AMEE方法在局部范圍內(nèi)膨脹所得的最純像元可能并不準(zhǔn)確。

圖2 不同結(jié)構(gòu)元素覆蓋像元的空間分布Fig.2 Spatial distribution of pixels covered by different structural elements

如果在結(jié)果元素范圍內(nèi)采用PPI算法尋找最純像元，三角形框abc和aef的3個頂點均能被投影到隨機生成的直線的兩端，但每個頂點投影到直線兩端的次數(shù)并不相同，真實端元a投影到直線兩端的次數(shù)不一定是最多。但如果在圖像上移動結(jié)構(gòu)元素，比如從三角形框abc變化到三角形框aef，原來的頂點b、c變?yōu)榱薳、f，也就是說隨著結(jié)構(gòu)元素的變化，偽端元也是變化的，但真實端元a卻依舊能夠投影到線段兩端。也就說只要結(jié)構(gòu)元素內(nèi)存在真實的純像元，就一定能被PPI算法投影到直線的兩端，而偽端元則隨著結(jié)構(gòu)元素的改變可能不會再次被PPI算法找到。

如果對結(jié)構(gòu)元素內(nèi)的像元采用PPI算法，記錄每個像元的純度指數(shù)，那么隨著結(jié)果元素的增大和改變，偽純像元的純度指數(shù)不會累計增大，而實際純像元的純度指數(shù)則不斷累計增大。因此，可以利用PPI指數(shù)代替AMEE算法中的MEI指數(shù)來對形態(tài)學(xué)中的膨脹和腐蝕操作重新進(jìn)行拓展定義。

膨脹操作：將結(jié)構(gòu)元素覆蓋范圍內(nèi)的像元在隨機生成的直線上投影，統(tǒng)計每個像元落在線段兩端的次數(shù)，次數(shù)最多者被認(rèn)為純度最高，并把這個次數(shù)作為純度指標(biāo)賦給對應(yīng)的像元。

腐蝕操作：統(tǒng)計每次投影距離線段中點最近的像元，并記錄次數(shù)，次數(shù)最多者被認(rèn)為混合程度最嚴(yán)重。

本文提出的方法不需要進(jìn)行腐蝕操作，只需要進(jìn)行膨脹操作。記錄結(jié)構(gòu)元素內(nèi)每一個像元投影到直線兩端的次數(shù)，并賦值給對應(yīng)像元作為純度指數(shù)。如前所述，隨著結(jié)構(gòu)元素的增大和變化，只有真實的純像元的PPI值在不斷累計增大，偽端元的PPI值不會持續(xù)增大，最終形成一幅類似于AMEE方法中的MEI圖像的PPI圖像，端元就從具有較大PPI值的像元中選取。

傳統(tǒng)的PPI算法需要隨機性地生成投影直線，算法執(zhí)行過程中需要用戶主觀介入觀察，具體生成多少投影向量，需要由用戶主觀設(shè)定，沒有統(tǒng)一的標(biāo)準(zhǔn)。本文提出的PPI-AMEE算法只在結(jié)構(gòu)元素內(nèi)部運行PPI算法，將結(jié)構(gòu)元素的像元在特征空間中兩兩連線生成投影向量，這樣使得投影向量盡可能地分布在不同方向上，而不會像傳統(tǒng)PPI算法那樣在對整幅圖像的所有像元投影時，由于生成的投影向量是完全隨機的，有可能使得生成的投影向量主要集中在某些特定方向上，而使得最終結(jié)果并不準(zhǔn)確。而且由于結(jié)構(gòu)元素所包含的像元數(shù)量是確定的，因此兩兩連線生成投影向量的數(shù)量是確定的，而且隨著結(jié)構(gòu)元素的擴大，像元數(shù)目增多，生成的投影向量也相應(yīng)地自動增多。這樣以來，算法運行過程中不再需要用戶主觀設(shè)定投影向量的數(shù)目，在一定程度上提高了算法的自動程度。

這里定義結(jié)構(gòu)元素K內(nèi)某個像元f(x,y)投影后所得PPI值的累計之和為

(9)

相應(yīng)的膨脹和腐蝕操作可以分別拓展性地定義為

dPPI(x,y)=(f?K)(x,y)=

(10)

ePPI(x,y)=(f⊕K)(x,y)=

(11)

PPI-AMEE算法的具體步驟如下：

(1)先對整幅圖像進(jìn)行最小噪聲分離變換(minimum noise fraction，MNF)變換進(jìn)行降維和去噪處理，估算端元數(shù)目m；

(2)設(shè)定結(jié)構(gòu)元素的最小尺寸Kmin和最大尺寸Kmax，得出最大迭代次數(shù)Imax；

(3)令i=1，所有像元的PPI初值P(f(x,y),K)=0，從最小結(jié)構(gòu)元素Kmin開始執(zhí)行；

(4)按照PPI算法拓展定義的形態(tài)學(xué)結(jié)構(gòu)算子進(jìn)行膨脹操作，得到結(jié)構(gòu)元素內(nèi)各像元的純度指數(shù)；

(5)i=i+1，如果i=Imax，順序執(zhí)行步驟(6)，否則增大結(jié)構(gòu)元素K，跳回步驟(4)執(zhí)行；

(6)輸出PPI圖像，參考PPI值較大的m個像元確定為端元。

算法的流程如圖3所示。

3 試驗結(jié)果與分析

首先采用模擬數(shù)據(jù)比較PPI、AMEE、PPI-AMEE這3種算法的性能，選取美國地質(zhì)勘探局(United States Geological Survey，USGS)光譜庫中的4種礦物光譜，如圖4所示，按照Dirichlet分布生成一個隨機數(shù)字序列來構(gòu)造混合像元的豐度矩陣，生成大小為200×200的模擬圖像，并在模擬圖像中分別添加信噪比(SNR)分別為10 db、20 db、30 db、40 db、50 db的高斯噪聲。表1給對比分析了不同信噪比情況下，PPI、AMEE、PPI-AMEE所提取的4個端元與USGS光譜庫中對應(yīng)光譜之間的光譜角距離(spectral angle distance,SAD)的均值、光譜信息散度(spectral information divergence,SID)的均值。

圖3 PPI-AMEE算法流程Fig.3 Flow chart of PPI-AMEE algorithm

圖4 構(gòu)建模擬數(shù)據(jù)的4種光譜Fig.4 Four spectra of simulated data

本文的PPI算法進(jìn)行了1000次隨機投影，由于PPI算法運行時間與隨機投影的次數(shù)有關(guān)，投影次數(shù)越多，運行時間越長，因此比較PPI算法與其他算法的運行效率不具有意義。試驗中，AMEE算法和PPI-AMEE算法迭代時，設(shè)置它們的最小結(jié)構(gòu)元素和最大結(jié)構(gòu)元素均為3×3和15×15，比較它們不同信噪比情況下的運行時間，試驗結(jié)果見表1。

表1顯示，隨著信噪比的提高，3種算法所提取的4個端元與USGS光譜庫對應(yīng)光譜之間的光譜角距離均值、光譜信息散度均值逐漸變小，也就是說端元提取精度逐漸變高。從算法運行時間來看，在信噪比較低時PPI-AMEE的運行效率比AMEE稍低，但隨著信噪比的提高PPI-AMEE的運行效率逐漸優(yōu)于AMEE。

表1 不同信噪比情況下PPI、AMEE、PPI-AMEE的模擬數(shù)據(jù)試驗結(jié)果Tab.1 Experimental results of simulated data for PPI,AMEE and PPI-AMEE under different SNR

真實高光譜數(shù)據(jù)試驗中采用的是美國內(nèi)華達(dá)州Cuprite地區(qū)的AVIRIS高光譜數(shù)據(jù)，該地區(qū)有詳細(xì)的地物調(diào)查報告，是檢驗各種分類算法和混合像元分解算法的典型數(shù)據(jù)[23-25]。該數(shù)據(jù)包含224個波段，波長范圍為400～2500 nm，去除其中受水蒸氣影響較大和信噪比較低的波段，剩余188個波段。選擇15、69、157波段作為紅、綠、藍(lán)通道合成假彩色圖像，如圖5所示。

圖5 Cuprite地區(qū)的AVIRIS數(shù)據(jù)假彩色合成圖Fig.5 False-color composition image of AVIRIS data in Cuprite region

利用Hysime算法[26]估計端元的數(shù)目為9。

這里重點研究對比4種礦物：明礬石(alunite)、鈣芒硝(buddingtonite)、方解石(calcite)和高嶺石(kaolinite)，因為這4種礦物都有比較明顯的純像元分布區(qū)域。

試驗中分別采用AMEE算法、PPI算法和PPI-AMEE算法提取9個端元，然后計算其與USGS光譜庫中的4種礦物alunite、buddingtonite、calcite、kaolinite標(biāo)準(zhǔn)光譜的光譜角。AMEE算法和PPI-AMEE算法中最小結(jié)構(gòu)元素和最大結(jié)構(gòu)元素分別設(shè)置為3×3和15×15。用PA1、PA2，…，PA9表示PPI-AMEE算法提取的9個端元，表2所示為PPI-AMEE算法提取端元與USGS標(biāo)準(zhǔn)光譜之間的光譜角混淆矩陣。

從表2中看出，PA3與alunite標(biāo)準(zhǔn)光譜的光譜角最小，因此PA3即被確認(rèn)為PPI-AMEE算法提取的對應(yīng)alunite的端元。同理，PA5、PA4、PA9分別對應(yīng)buddingtonite、calcite、kaolinite。

然后，分別單獨采用AMEE算法和PPI算法提取Cuprite地區(qū)alunite、buddingtonite、calcite、kaolinite的端元光譜，連同上述PPI-AMEE算法所提取的端元光譜曲線一起示于圖6中，并與USGS光譜庫中標(biāo)準(zhǔn)光譜曲線進(jìn)行對比。

表2 PPI-AMEE提取端元與USGS標(biāo)準(zhǔn)光譜之間的光譜角混淆矩陣Tab.2 Confusion matrix of spectral angle distances between PPI-AMEE endmembers and USGS reference spectra

表3中的試驗數(shù)據(jù)比較了PPI、AMEE、PPI-AMEE這3種算法所提取的4個端元光譜與各自對應(yīng)的標(biāo)準(zhǔn)光譜之間的光譜角距離、光譜信息散度。

圖6 PPI、AMEE、PPI-AMEE提取的端元與USGS標(biāo)準(zhǔn)光譜的對比Fig.6 Comparison of enmembers extracted by PPI,AMEE,PPI-AMEE with USGS standard spectra

表3 不同信噪比情況下PPI、AMEE、PPI-AMEE的真實高光譜數(shù)據(jù)試驗結(jié)果Tab.3 Experimental results of real hyperspectral data for PPI,AMEE and PPI-AMEE

真實數(shù)據(jù)的試驗結(jié)果表明，與AMEE算法和PPI算法相比，PPI-AMEE算法端元與USGS標(biāo)準(zhǔn)光譜之間的光譜角距離和光譜信息散度最小，其端元提取精度總體上優(yōu)于其他兩種算法。

4 結(jié) 論

本文提出的PPI-AMEE算法融合了AMEE算法和PPI算法的設(shè)計思想，采用PPI算法對數(shù)學(xué)形態(tài)學(xué)中的膨脹和腐蝕操作作了拓展性的定義，在結(jié)構(gòu)元素內(nèi)通過累計各像元的PPI指數(shù)尋找最純像元。PPI算法只是把整幅高光譜影像當(dāng)作無序向量進(jìn)行投影的，僅僅利用了圖像的光譜信息。PPI-AMEE算法在結(jié)構(gòu)元素局部范圍內(nèi)運行PPI算法，然后結(jié)合數(shù)學(xué)形態(tài)學(xué)的膨脹操作對整幅圖像進(jìn)行處理，同時兼顧了圖像的空間信息和光譜信息。試驗結(jié)果表明，本文提出的PPI-AMEE算法的端元提取精度優(yōu)于AMEE算法和PPI算法。