基于均值漂移聚類的端元束提取

陳立偉，邱艷芳，朱海峰，王立國

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

由于空間分辨率的限制以及地物分布的復(fù)雜多樣性，混合像元普遍存在于高光譜遙感圖像中。因此，實(shí)際應(yīng)用中必須解決光譜解混的問題。光譜解混主要分為端元提取和豐度反演2個(gè)步驟。近年來，有多種端元提取方法相繼提出，最具有代表性的為像元純度指數(shù)算法(pixel purity index，PPI)[1]，內(nèi)部最大體積算法 (N-FINDR)[2]，單形體增長算法(the simplex growing algorithm，SGA)[3]，和頂點(diǎn)成分分析(vertex component analysis，VCA）[4]等。這些端元提取方法只能從高光譜數(shù)據(jù)中提取單一端元代表一類地物。由于光照輻射條件的變化、時(shí)空等自然因素變化的不可測(cè)影響，高光譜遙感圖像中地物普遍存在光譜變異性[5]。此時(shí)，單一端元對(duì)地物類別的代表性有限，影響對(duì)真實(shí)地物豐度反演的精度。因此，從高光譜遙感圖像中為每一種地物類別提取多個(gè)代表端元，并形成端元束（endmember bundle，EB）[6-7]，是處理光譜變異性的一種簡單有效方法。

端元束提取分為2個(gè)步驟：獲取候選端元，對(duì)候選端元進(jìn)行聚類。聚類的目的是形成每種地物的代表性端元光譜集合。在聚類階段，目前常用的算法就是K 均值（K-means）聚類[8]，K-means算法對(duì)于大型的數(shù)據(jù)計(jì)算簡單高效，且時(shí)間和空間復(fù)雜度都較低，在高光譜圖像分類中也有應(yīng)用[9-10]。但是K均值聚類算法需要預(yù)先知道類簇的數(shù)目，在光譜聚類中表現(xiàn)為需要預(yù)先知道高光譜遙感圖像中地物的類別數(shù)目。但是地物類別數(shù)目的估計(jì)通常是很困難的[11]。地物類別數(shù)目的估計(jì)誤差會(huì)對(duì)聚類結(jié)果產(chǎn)生直接的影響，從而增加多端元光譜混合分析的計(jì)算量[12]。

1 候選端元提取

端元提取的方法有很多，大致可以分為基于幾何學(xué)和基于統(tǒng)計(jì)學(xué)的算法。這些算法中大多數(shù)都是利用光譜的信息進(jìn)行端元提取，而目前已有研究者對(duì)空間信息加以利用，與光譜信息結(jié)合提取圖像中的端元，將此類方法運(yùn)用于端元束提取的潛力已經(jīng)被證明[13]，因此文中候選端元的提取采取空間信息和光譜信息相結(jié)合的方法[13]，主要分為3個(gè)步驟。

1）PPI（pixel purity index）預(yù)處理。運(yùn)用 PPI算法對(duì)候選端元進(jìn)行粗略的篩選，此處閾值T設(shè)置為0，生成隨機(jī)向量的數(shù)目設(shè)置為K（地物種類數(shù)目為10～15時(shí)，K設(shè)置為10 000）。經(jīng)過PPI算法的計(jì)算會(huì)提取出大量的候選端元，此步驟也是為了減少下一步驟的計(jì)算量和復(fù)雜程度。

2）同質(zhì)性指數(shù)（homogeneity index，HI）的計(jì)算。文獻(xiàn)[14-15]中指出，端元一般位于空間同質(zhì)區(qū)域內(nèi)，混合像元一般位于過渡區(qū)域或者不均勻的區(qū)域。在此基礎(chǔ)上，文獻(xiàn)[16]將相鄰像素的相似度定義為HI，通過對(duì)像素點(diǎn)與鄰域內(nèi)其他像素的光譜信息散度（spectral information divergence，SID）來定量測(cè)量HI。HI值越小，像素點(diǎn)為純像元的概率越大。

3）基于區(qū)域的端元選擇。在該步驟中，根據(jù)HI的閾值來選擇端元。在不同的環(huán)境，對(duì)于不同的材料，端元的光譜變異程度不同。為此，對(duì)整個(gè)圖像進(jìn)行分塊，自適應(yīng)地為不同的塊選擇不同的閾值，從不同塊中獲得端元集合，其中初始閾值由HI值的統(tǒng)計(jì)直方圖確定。根據(jù)選候端元的比例自適應(yīng)地調(diào)整該閾值，圖像塊的大小是依據(jù)圖像的復(fù)雜程度來劃分的。如果地物分布廣泛且均勻，則塊的尺寸可以稍大[17]，反之亦然。

2 基于均值漂移聚類的端元束提取

2.1 均值漂移聚類

均值漂移（mean shift）算法是一種無參密度估計(jì)算法或稱核密度估計(jì)算法[18]，在聚類、圖像分割、跟蹤等方面有著廣泛的應(yīng)用。Mean shift是一個(gè)向量，它的方向指向當(dāng)前點(diǎn)上概率密度梯度的方向。在聚類中，該算法完全依靠特征空間中的樣本點(diǎn)進(jìn)行分析，不需要任何先驗(yàn)知識(shí)，數(shù)據(jù)集中的每一點(diǎn)都可以作為初始點(diǎn)，它能對(duì)任何維度、任何分布的采樣點(diǎn)進(jìn)行快速聚類，迭代效率高。

均值漂移指的就是聚類中心根據(jù)偏移向量的均值進(jìn)行移動(dòng)，偏移均值為式中：Sh指的是以數(shù)據(jù)x為中心，半徑為h的圓形區(qū)域；k是指區(qū)域Sh內(nèi)所包含數(shù)據(jù)點(diǎn)的個(gè)數(shù)；xi是指圓形區(qū)域范圍內(nèi)的第i個(gè)數(shù)據(jù)點(diǎn)。的方向就是聚類中心移動(dòng)的方向，聚類中心的移動(dòng)距離就是的模。

某一時(shí)刻聚類中心的計(jì)算為

式中：Mt是指在t狀態(tài)下的偏移均值；xt、xt+1分別指的是狀態(tài)t和t+1時(shí)的聚類中心。

2.2 基于光譜角距離準(zhǔn)則的均值漂移聚類

傳統(tǒng)的均值漂移聚類算法通常使用歐氏距離作為聚類的判決準(zhǔn)則，同一地物的變異光譜在空間中的位置比較靠近，但是當(dāng)變異較大時(shí)，變異端元的位置也可能距離較遠(yuǎn)，傳統(tǒng)的均值漂移算法就不再適用。同一物質(zhì)的變異光譜之間在多維空間中的光譜角較小，故本文在均值漂移算法采用光譜角距離（spectral angle distance，SAD）作為聚類的判決準(zhǔn)則。2條光譜 xa和 xb之間的SAD計(jì)算為

2）標(biāo)記以該中心光譜為中心，光譜角距離為半徑，半徑為h的圓形范圍內(nèi)所有的光譜數(shù)據(jù)，記做集合，將這些點(diǎn)劃分為一個(gè)簇。同時(shí)，記這些光譜的標(biāo)記次數(shù)為1。

6）重復(fù)步驟1）～5）直到所有的樣本都被標(biāo)記訪問。根據(jù)每個(gè)類，對(duì)每個(gè)樣本的訪問頻率，取訪問頻率最大的那個(gè)類，作為當(dāng)前點(diǎn)集的所屬類。

均值漂移聚類方法不需要預(yù)先知道地物的種類數(shù)目，完全依靠現(xiàn)有的光譜數(shù)據(jù)進(jìn)行分析，盡可能避免由預(yù)估地物類別數(shù)目算法帶來的誤差，更準(zhǔn)確地形成端元束。下面通過實(shí)驗(yàn)對(duì)改進(jìn)前后的均值漂移聚類方法以及K均值聚類算法在端元束聚類方面進(jìn)行了比較。

3 實(shí)驗(yàn)結(jié)果與討論

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

3.1.1 變異光譜的合成

要合成混合高光譜圖像，首先需要合成變異光譜。本文從USGS光譜庫中選取4種地物光譜作為原始光譜，如表1所示?？紤]光照引起的變化，本文基于端元擾動(dòng)理論去合成光譜[19]。文中將各地物光譜的變異系數(shù)分別為設(shè)置為：0.1、0.35、0.25和0.25；使用的USGS光譜庫中光譜波長在380～2 500 nm，共224個(gè)波段，合成后光譜如圖1所示。

表1 選擇的地物名稱

圖1 合成的變異光譜

3.1.2 混合高光譜圖像的合成

實(shí)驗(yàn)中使用的是大小為1 2 8×128×224，由上述變異光譜線性混合而成的圖像，并在圖像中加入40 dB的高斯白噪聲。合成的混合高光譜圖像中存在純像元，將變異光譜和原始光譜均視為端元，序號(hào)為1～3的端元豐度大小均以各自圓心向外發(fā)散，圓心在高光譜圖像中的位置分別為[50,30]、[100,64]以及[50,100]，在距離圓心越遠(yuǎn)位置的混合像元中，該端元所占的比例越小。序號(hào)為4的端元?jiǎng)t與之相反，豐度分布情況如圖2。

圖2 豐度分布

3.1.3 實(shí)驗(yàn)結(jié)果及分析

在端元集提取部分，有幾點(diǎn)需要指出：1）PPI預(yù)處理算法中的閾值設(shè)置為T=0，可以確保提取到的光譜中包含了所有的4類地物；2）PPI預(yù)處理只是做一個(gè)粗略的篩選，其中隨機(jī)向量的數(shù)量遠(yuǎn)遠(yuǎn)小于傳統(tǒng)PPI算法中隨機(jī)向量的數(shù)量，本次實(shí)驗(yàn)中將K設(shè)置為50 000；3）在實(shí)驗(yàn)中，根據(jù)HI直方圖將數(shù)據(jù)塊的大小設(shè)置為20×20，并將閾值設(shè)置為0.013（如圖3所示）。

在光譜聚類部分，1）由于混合圖像中選取了4類地物，故K-means聚類算法中聚類中心的個(gè)數(shù)設(shè)置為4；2）均值漂移聚類算法中的半徑設(shè)置為0.9，中心點(diǎn)移動(dòng)距離誤差設(shè)置為0.14，與聚類中心之間的歐氏距離誤差設(shè)置為0.2；3）改進(jìn)后均值漂移聚類算法中的半徑設(shè)置為2.26，中心點(diǎn)移動(dòng)距離誤差設(shè)置為0.14，與聚類中心之間的歐氏距離誤差設(shè)置為0.16；4）聚類后的光譜束與USGS光譜庫中的光譜進(jìn)行對(duì)比，選擇光譜角距離最小的地物作為光譜束所代表的地物類別。

圖3 合成數(shù)據(jù)的HI直方圖

按照改進(jìn)后的均值漂移算法進(jìn)行聚類，其結(jié)果如圖4所示，圖中依次為Actinolite、Brucite、Corundum以及Goethite。

圖4 改進(jìn)后的均值漂移聚類結(jié)果

本文將通過光譜類內(nèi)部的分散程度和光譜類之間的分散程度來度量聚類性能的優(yōu)劣。光譜類內(nèi)部分散程度越小且光譜類之間的分散程度越大則聚類效果越好。實(shí)驗(yàn)中以各類光譜中光譜的平均值作為類中心光譜，光譜內(nèi)部分散程度采用的是類內(nèi)各光譜到類中心光譜的SAD值來表示，光譜類間的分散程度采用的是各類的中心光譜之間的SAD來度量，SAD值越小說明越聚集。因?yàn)橐唤M光譜代表的是一類地物，所以每一類地物都會(huì)有相應(yīng)的一組類內(nèi)SAD值以及類間SAD值。

分別使用K-means和未改進(jìn)的meanshift的算法進(jìn)行對(duì)比，每組類內(nèi)SAD值的統(tǒng)計(jì)信息結(jié)果如圖5所示，其中的數(shù)字表示的是該類光譜束中的所擁有的光譜的數(shù)目。

各類光譜束的中心光譜間的SAD值如表2所示。序號(hào)所對(duì)應(yīng)的端元與表1中一致。

圖5 地物類內(nèi)SAD值的統(tǒng)計(jì)信息

表2 各類的中心光譜間的SAD

實(shí)驗(yàn)結(jié)果表明，在未知端元數(shù)目的情況下，均值聚類算法也能聚類出正確類別數(shù)目的地物；由圖5可以看出，Brucite使用3種方法的聚類效果完全相同；Corundum使用改進(jìn)后的均值漂移算法時(shí)的最大觀測(cè)點(diǎn)較小，離群點(diǎn)較少，整體數(shù)據(jù)比較集中；Actinolite使用改進(jìn)后的均值漂移算法聚類式的箱圖更加的矮胖，數(shù)據(jù)集中效果更好。Goethite光譜束中，改進(jìn)后的meanshift聚類較K-means聚類箱圖更加矮胖，數(shù)據(jù)整體比較集中；較未改進(jìn)的meanshift算法整體觀測(cè)點(diǎn)更小，表明光譜束之間的光譜角更小。表2中顯示了各類中心光譜之間的SAD值，使用K均值算法時(shí)，Actinolite與Goethite、Corundum與Goethite之間中心光譜的SAD值都較小，分別為0.177 5和0.087 6，類間的分離程度不高；改進(jìn)后的均值漂移算法整體上的類間SAD值略高于K均值算法以及傳統(tǒng)均值漂移算法，類與類之間的差異性更大，分離程度更高。

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

為了評(píng)估均值漂移聚類算法在不同的分析場(chǎng)景中的有效性，接下來用真實(shí)圖像進(jìn)行實(shí)驗(yàn)。實(shí)驗(yàn)所采用的覆蓋Cuprite場(chǎng)景的真實(shí)數(shù)據(jù)集于1997年被機(jī)載可見紅外成像光譜儀（AVIRIS）傳感器捕獲。Cuprite原始數(shù)據(jù)有224個(gè)波段，范圍覆蓋370～2 480 nm。此處選擇了Cuprite數(shù)據(jù)中的50個(gè)波段，大小為190×250的數(shù)據(jù)。由于USGS光譜庫有大量的可用于現(xiàn)場(chǎng)地面真實(shí)數(shù)據(jù)，該數(shù)據(jù)集已成為用于算法評(píng)估的常用基準(zhǔn)數(shù)據(jù)集。Swayze和Clark[20]也制作了一份關(guān)于該地區(qū)基本事實(shí)的報(bào)告。端元集提取部分需要指出的是：1）實(shí)驗(yàn)中的PPI預(yù)處理算法將閾值設(shè)置為0，生成的隨機(jī)向量個(gè)數(shù)為12 000；2）Cuprite數(shù)據(jù)的HI直方圖如圖6所示，由圖6將數(shù)據(jù)塊的大小設(shè)置為20×20，并將閾值設(shè)置為0.014；3）在聚類部分，地物類別的數(shù)目為12[20]，故K-means聚類算法中中心點(diǎn)個(gè)數(shù)設(shè)置為12；均值漂移聚類算法的半徑設(shè)置為5；中心點(diǎn)移動(dòng)距離誤差設(shè)置為3.2；聚類距離誤差設(shè)置為13；以上數(shù)值在改進(jìn)后的均值漂移聚類算法中分別設(shè)置為 0.03、0.002、0.037。

通過2種均值漂移算法聚類之后均可得到11種地物類型，如圖7和8所示。

圖6 真實(shí)圖像數(shù)據(jù)的HI直方圖

圖7 均值漂移聚類結(jié)果

圖8 改進(jìn)后的均值漂移聚類結(jié)果

K-means算法聚類之后共得到12個(gè)光譜束，將這些光譜束與USGS光譜庫中的地物光譜進(jìn)行對(duì)比，同樣得到11種地物，其中與改進(jìn)后的均值漂移算法所得的地物類型相同的地物有7種，如圖9所示，依次為：Oligoclase HS110.3B、Mizzonite HS350.3B HLSep、 Lazurite HS418.3B、Desert_Varnish ANP90-14、 Mizzonite NMNH113775-1、Kaolin/Smect KLF506 95%K、Montmorillonite CM26。

根據(jù)結(jié)果圖7可明顯看出，在未經(jīng)改進(jìn)的meanshift算法聚類中，有明顯的聚類錯(cuò)誤情況出現(xiàn)，聚類效果不如K-means算法以及改進(jìn)后的meanshift算法，因此接下來只將K-means以及改進(jìn)后的meanshift算法進(jìn)行比較。將2種算法所得的這7種相同類型的地物的SAD統(tǒng)計(jì)信息進(jìn)行比較，各端元束內(nèi)部的SAD信息如圖10所示（除去類內(nèi)只有一條光譜的地物，Mizzonite NMNH 113775-1和 Montmorillonite CM26）；各類光譜束的中心光譜間的SAD的信息如圖11所示。

圖9 K-means聚類結(jié)果

由圖10可以看出，對(duì)meanshift算法進(jìn)行改進(jìn)后，聚類結(jié)果中所有地物均無離散點(diǎn)的出現(xiàn)，整體而言，聚合程度較好；其中Desert_Varnish ANP90-14的類內(nèi)SAD值在最大值、最小值以及中值等方面小于K-means算法，具有明顯更優(yōu)的聚類效果；Oligoclase HS110.3B和Lazurite HS418.3B 2種地物的聚類情況略遜于K-means算法。由圖11可以看出，2種算法得出的各個(gè)端元束的中心光譜之間SAD值的分布情況大致相同，在預(yù)先未知地物種類的情況下，改進(jìn)后的均值漂移算法所得到的結(jié)果較K-means算法無明顯缺點(diǎn)存在。

圖10 2種聚類方法的5種相同地物的類內(nèi)SAD統(tǒng)計(jì)信息

圖11 各類的中心光譜間的SAD統(tǒng)計(jì)信息

4 結(jié)論

端元束提取中的聚類影響多端元光譜混合分析的效率。本文對(duì)傳統(tǒng)的均值漂移聚類算法進(jìn)行改進(jìn)，將改進(jìn)后的方法應(yīng)用于光譜聚類中，并將其結(jié)果與K-means聚類方法以及傳統(tǒng)的均值漂移算法結(jié)果進(jìn)行比較。

1)在不需要預(yù)先知道地物種類數(shù)的前提條件下，均值漂移聚類算法能夠依靠數(shù)據(jù)自身的特點(diǎn)對(duì)數(shù)據(jù)進(jìn)行聚類。

2)通過模擬數(shù)據(jù)和真實(shí)數(shù)據(jù)實(shí)驗(yàn)的驗(yàn)證，基于光譜角距離的均值漂移算法較原方法更加適合光譜聚類，且在整體效果上較K-means算法有一定的優(yōu)勢(shì)，聚類后的數(shù)據(jù)集類內(nèi)的離散程度更低，類之間的分離程度更高。

端元束的有效聚類，有利于提高多端元光譜混合分析的效率，對(duì)于高光譜遙感定量解譯具有重要意義。