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

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

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

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

1 候選端元提取

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

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

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

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

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

2.1 均值漂移聚類

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

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

某一時刻聚類中心的計算為

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

2.2 基于光譜角距離準則的均值漂移聚類

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

2）標記以該中心光譜為中心，光譜角距離為半徑，半徑為h的圓形范圍內所有的光譜數據，記做集合，將這些點劃分為一個簇。同時，記這些光譜的標記次數為1。

6）重復步驟1）～5）直到所有的樣本都被標記訪問。根據每個類，對每個樣本的訪問頻率，取訪問頻率最大的那個類，作為當前點集的所屬類。

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

3 實驗結果與討論

3.1 模擬數據實驗

3.1.1 變異光譜的合成

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

表1 選擇的地物名稱

圖1 合成的變異光譜

3.1.2 混合高光譜圖像的合成

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

圖2 豐度分布

3.1.3 實驗結果及分析

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

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

圖3 合成數據的HI直方圖

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

圖4 改進后的均值漂移聚類結果

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

分別使用K-means和未改進的meanshift的算法進行對比，每組類內SAD值的統計信息結果如圖5所示，其中的數字表示的是該類光譜束中的所擁有的光譜的數目。

各類光譜束的中心光譜間的SAD值如表2所示。序號所對應的端元與表1中一致。

圖5 地物類內SAD值的統計信息

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

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

3.2 真實數據實驗

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

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

圖6 真實圖像數據的HI直方圖

圖7 均值漂移聚類結果

圖8 改進后的均值漂移聚類結果

K-means算法聚類之后共得到12個光譜束，將這些光譜束與USGS光譜庫中的地物光譜進行對比，同樣得到11種地物，其中與改進后的均值漂移算法所得的地物類型相同的地物有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。

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

圖9 K-means聚類結果

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

圖10 2種聚類方法的5種相同地物的類內SAD統計信息

圖11 各類的中心光譜間的SAD統計信息

4 結論

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

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

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

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