一種基于譜聚類算法的高光譜遙感圖像分類方法*

楊隨心，耿修瑞，楊煒暾，趙永超，盧曉軍

(1 中國(guó)科學(xué)院電子學(xué)研究所 中國(guó)科學(xué)院空間信息處理與應(yīng)用系統(tǒng)技術(shù)重點(diǎn)實(shí)驗(yàn)室，北京 100190； 2 中國(guó)科學(xué)院大學(xué)，北京 100049； 3 中國(guó)國(guó)際工程咨詢公司，北京 100048)

高光譜圖像一般含有數(shù)十至數(shù)百個(gè)波段，其數(shù)據(jù)可以看作是一個(gè)三維立方體，其中每一個(gè)像素點(diǎn)抽取出來(lái)都是一條連續(xù)的光譜曲線，相鄰波段之間的波長(zhǎng)相差為納米級(jí)，其精細(xì)的光譜分辨率使得高光譜圖像應(yīng)用越來(lái)越廣泛[1]，例如在地質(zhì)制圖、礦石尋找、環(huán)境大氣監(jiān)測(cè)、農(nóng)林森林調(diào)查、海洋生物研究與海水分析等領(lǐng)域都有重要應(yīng)用[2]。隨著高光譜圖像應(yīng)用越來(lái)越廣泛，其圖像分析處理也變得越來(lái)越重要。高光譜圖像分類是其信息分析的重要手段之一，分類主要有兩種，有監(jiān)督分類和無(wú)監(jiān)督分類。有監(jiān)督分類方法眾多，如SAM[3]、最大似然[4]、SVM[5]和神經(jīng)網(wǎng)絡(luò)[6]等方法，分類精度隨著方法的發(fā)展越來(lái)越高，但是有監(jiān)督的分類必須以圖像存在樣本標(biāo)簽為前提。而在實(shí)際應(yīng)用中，獲得的高光譜圖像往往是沒(méi)有先驗(yàn)信息的，在這種情況下，無(wú)監(jiān)督聚類就顯得十分重要[7]。

在實(shí)際應(yīng)用中，常用的無(wú)監(jiān)督聚類方法有K-means算法[8-9]、均值漂移(mean-shift)算法[10]、迭代自組織數(shù)據(jù)分析(ISODATA)算法[11]和DBSCAN(density-based spatial clustering of applications with noise)方法[12]等。其中，K-means算法和ISODATA算法對(duì)于初始值的選取十分敏感，容易陷入局部最優(yōu)，且只有在數(shù)據(jù)為球狀或者接近球狀時(shí)聚類效果較好；Mean-Shift算法時(shí)間復(fù)雜度高，且半徑參數(shù)的選取會(huì)對(duì)聚類結(jié)果產(chǎn)生很大的影響；DBSCAN方法雖然不用指定最終聚類類別數(shù)目，但是當(dāng)空間聚類的密度不均勻、聚類間距差相差很大時(shí)，鄰域半徑和對(duì)象鄰域內(nèi)樣本數(shù)量閾值選取困難，聚類質(zhì)量較差。近年來(lái)，譜聚類方法[13-14]在眾多聚類算法中顯示出聚類精度高穩(wěn)定性好的優(yōu)勢(shì)，其應(yīng)用也越來(lái)越廣泛，但是在高光譜圖像中很少見(jiàn)到譜聚類算法的應(yīng)用。譜聚類方法進(jìn)行聚類的關(guān)鍵是數(shù)據(jù)間相似度矩陣的構(gòu)造以及特征值和特征向量的求解，相似度矩陣的大小取決于參與聚類的樣本個(gè)數(shù)，當(dāng)數(shù)據(jù)量非常大時(shí)，該矩陣占有很大的內(nèi)存，且特征值的求解面臨巨大壓力，所以該方法適用于樣本量較少的數(shù)據(jù)，而高光譜圖像的數(shù)據(jù)量非常大，致使譜聚類方法無(wú)法在高光譜圖像聚類中進(jìn)行應(yīng)用。針對(duì)這一問(wèn)題，本文提出結(jié)合K-Means的譜聚類方法(K-means spectral clustering，KSC)，有效結(jié)合K-means算法，將譜聚類應(yīng)用到高光譜圖像中，對(duì)其實(shí)現(xiàn)較高總體分類精度的聚類。

1 譜聚類方法

譜聚類方法是一種基于圖的聚類方法，能夠?qū)θ我庑螤畹臄?shù)據(jù)進(jìn)行最優(yōu)劃分。由于這一優(yōu)點(diǎn)，這些年譜聚類方法的應(yīng)用領(lǐng)域越來(lái)越廣泛。譜聚類算法以圖論為基礎(chǔ)，所有數(shù)據(jù)可以看成是無(wú)向帶權(quán)圖G=(V,E)，其中，V={v1,v2,…,vm}是頂點(diǎn)集合，E是邊的集合。任意2個(gè)頂點(diǎn)vi,vj之間的邊對(duì)應(yīng)的權(quán)值為wij，表示2個(gè)點(diǎn)之間的相似程度。對(duì)于所有數(shù)據(jù)點(diǎn)，任意2點(diǎn)之間的相似度構(gòu)成相似度矩陣W。譜聚類算法的分類思想為找到數(shù)據(jù)集中類內(nèi)相似度最大而類間相似度最小的劃分。計(jì)算圖劃分準(zhǔn)則的最優(yōu)解是一個(gè)NP難問(wèn)題，所以一般將此類問(wèn)題轉(zhuǎn)化為求解相似度矩陣的譜分解問(wèn)題，再利用譜分解得到合適的特征向量來(lái)描述數(shù)據(jù)的低維結(jié)構(gòu)，然后在低維空間再利用K-means等經(jīng)典方法得到最終的聚類結(jié)果。

對(duì)于輸入的具有m個(gè)數(shù)據(jù)點(diǎn)的數(shù)據(jù)集，X={x1,x2,…,xm}，可根據(jù)其中任意2點(diǎn)間的相似度建立相似度矩陣W∈m×m，任意2點(diǎn)間相似度定義如下

(1)

式中：σ為尺度參數(shù)，本文中σ的大小按照Z(yǔ)elnik-Manor和Perona在Self-Tuning譜聚類方法[15]中提出的“l(fā)ocal Scaling”思想進(jìn)行選取。

記矩陣D為度矩陣，度矩陣為對(duì)角矩陣，將W的每行元素相加得到度矩陣D。記L∈m×m為拉普拉斯矩陣：

L=D-W.

(2)

假定{x1,x2,…,xm}可以分為2類，分別記為A，B，設(shè)

(3)

比例割(Ratio Cut)函數(shù)[16]為

(4)

式中：|A|為A類數(shù)據(jù)點(diǎn)數(shù)目，|B|為B類數(shù)據(jù)點(diǎn)數(shù)目。顯然最小化比例割對(duì)應(yīng)著一個(gè)最佳的二分類問(wèn)題。假定有m維向量f={f1,f2,…,fm}，它的每個(gè)元素代表當(dāng)前數(shù)據(jù)點(diǎn)的類別歸屬，根據(jù)文獻(xiàn)[14]，向量f滿足下列等式

fTLf=|V|Rcut(A,B)，

(5)

式中：|V|表示數(shù)據(jù)樣本點(diǎn)數(shù)量。對(duì)于2類比例割譜聚類轉(zhuǎn)化為如下模型

(6)

其中Tr為矩陣的跡。對(duì)于k類問(wèn)題，比例割優(yōu)化模型為

(7)

模型(7)是一個(gè)標(biāo)準(zhǔn)的求解跡最小化問(wèn)題，根據(jù)拉普拉斯矩陣性質(zhì)，該優(yōu)化問(wèn)題最優(yōu)解可轉(zhuǎn)化為求L矩陣最小的k-1個(gè)特征值對(duì)應(yīng)的特征向量。

對(duì)于一幅M×N×L大小的高光譜圖像，建立的相似度矩陣大小為(M×N)2，如果矩陣中每個(gè)值按照1個(gè)字節(jié)來(lái)計(jì)算，那么該矩陣所需內(nèi)存約(M×N)2×10-9G。以實(shí)際高光譜圖像舉例，如果一幅高光譜圖像大小為500×500×200，那么相似度矩陣為25萬(wàn)×25萬(wàn)，所需計(jì)算機(jī)內(nèi)存約為62.5 G，這是一般PC機(jī)無(wú)法承受的。針對(duì)這一問(wèn)題，本文提出相應(yīng)的解決方法。

2 K均值-譜聚類方法

K均值-譜聚類(K-means spectral clustering，KSC)方法是結(jié)合K-means算法與譜聚類方法的優(yōu)點(diǎn)而實(shí)現(xiàn)聚類的一種方法。K-means算法對(duì)初始值選取敏感而易陷入局部最優(yōu)，致使當(dāng)聚類類別數(shù)較少時(shí)，整體的分類精度較低，而當(dāng)聚類類別數(shù)非常多時(shí)，局部的分類精度較高。如對(duì)于圖1(a)所示的數(shù)據(jù)分布示意圖進(jìn)行K-means算法聚類，當(dāng)聚類類別數(shù)較少，如圖1(b)所示分為2類時(shí)，由于K-means算法以距離為分類判別準(zhǔn)則，對(duì)于該形狀分布的數(shù)據(jù)，K-means算法會(huì)出現(xiàn)錯(cuò)分情況，如圖中虛線圓內(nèi)的右側(cè)紅色部分錯(cuò)分；當(dāng)聚類類別數(shù)較多，如圖1(c)所示分為很多類時(shí)，對(duì)于分成的每一類其類內(nèi)分類精度都比較高。如此，可以利用這一點(diǎn)，首先對(duì)高光譜圖像數(shù)據(jù)進(jìn)行K-means聚類，將其聚成很多類，從而可得到相應(yīng)類別數(shù)的聚類中心，該聚類中心可作為相應(yīng)類別的代表，接下來(lái)可以對(duì)其進(jìn)行譜聚類處理，這樣的處理在不損害聚類精度的同時(shí)大大降低了譜聚類的空間復(fù)雜度。

圖1 不同聚類類別數(shù)時(shí)K-means聚類效果示意圖Fig.1 Schematic diagram of K-means clustering results with different cluster numbers

此外，高光譜圖像一般包含數(shù)十甚至數(shù)百個(gè)波段，波段間相關(guān)性強(qiáng)、冗余度高，且部分波段含有大量噪聲信息。因此，需要先對(duì)高光譜圖像進(jìn)行降維處理，這樣不但可以提高信噪比，而且可以降低處理的時(shí)間復(fù)雜度。常用的特征降維方法有主成分分析[17](principal component analysis，PCA)、最小噪聲分離[18](minimum noise fraction rotation，MNF Rotation)以及快速獨(dú)立成分分析[19](fast independent component Analysis，F(xiàn)astICA)等，本文在對(duì)高光譜圖像進(jìn)行KSC之前采用MNF對(duì)數(shù)據(jù)進(jìn)行特征降維處理。MNF+KSC方法的主要算法流程如下：

1)首先采用MNF方法對(duì)原始高光譜數(shù)據(jù)進(jìn)行特征降維處理。

2)然后采用K-means方法對(duì)原始圖像進(jìn)行粗聚類，聚成C類，得到C個(gè)聚類中心點(diǎn)，記為{μ1,μ2,…,μC}。

3)對(duì)于上述聚類中心{μ1,μ2,…,μC}，根據(jù)任意2點(diǎn)間的相似度，構(gòu)建相似度矩陣W，W的大小為C×C。

4)然后根據(jù)W矩陣分別求得度矩陣D和拉普拉斯矩陣L。

5)求解拉普拉斯矩陣L的最小特征值對(duì)應(yīng)的特征矢量，進(jìn)行最終分類。

6)對(duì)應(yīng)到2)中每個(gè)聚類中心對(duì)應(yīng)的原始數(shù)據(jù)，進(jìn)行類別標(biāo)記。

針對(duì)上述提出的KSC方法以及MNF+KSC方法，我們使用不同的數(shù)據(jù)進(jìn)行實(shí)驗(yàn)來(lái)驗(yàn)證本文方法良好的聚類性能。

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

本文分別使用模擬數(shù)據(jù)和真實(shí)的高光譜圖像數(shù)據(jù)進(jìn)行聚類實(shí)驗(yàn)，通過(guò)與K-means算法的聚類結(jié)果進(jìn)行比較，證實(shí)本文方法的優(yōu)勢(shì)。

3.1 模擬數(shù)據(jù)

本文采用的模擬數(shù)據(jù)來(lái)源為Clustering benchmark datasets[20]，最早在2007年ACM TKDD(ACM Transactions on Knowledge Discovery from Data)上的Clustering Aggregation文章中使用，該數(shù)據(jù)是二維數(shù)據(jù)，共788個(gè)像元，根據(jù)歐氏距離劃分可分為7類，具體分布如圖2(a)所示，目前該數(shù)據(jù)被廣泛應(yīng)用于聚類實(shí)驗(yàn)[21-22]。對(duì)該數(shù)據(jù)分別采用K-Means算法、譜聚類算法和KSC算法進(jìn)行聚類實(shí)驗(yàn)，結(jié)果如下。

圖2 K-means算法、譜聚類算法以及KSC算法聚類效果對(duì)比圖Fig.2 Clustering result comparison among K-means, spectral clustering， and KSC

實(shí)驗(yàn)結(jié)果表明，對(duì)于該模擬數(shù)據(jù)，K-means算法聚類的結(jié)果受初始值選取的影響，且該數(shù)據(jù)形狀分布不規(guī)則，所以K-means聚類的效果并不能令人滿意，而采用譜聚類方法和KSC方法進(jìn)行聚類，可以得到較好的聚類效果，該聚類結(jié)果初步證實(shí)KSC方法的可行性。

3.2 真實(shí)高光譜圖像數(shù)據(jù)

為進(jìn)一步驗(yàn)證本文算法的聚類性能，采用兩幅真實(shí)的高光譜圖像進(jìn)行定量分析實(shí)驗(yàn)。采用KSC和MNF+KSC方法對(duì)上述2個(gè)真實(shí)的高光譜數(shù)據(jù)進(jìn)行聚類實(shí)驗(yàn)，并與K-Means算法的聚類結(jié)果進(jìn)行比較，分別從聚類結(jié)果的總體分類精度、Kappa系數(shù)以及計(jì)算耗時(shí)3方面進(jìn)行比較。本文中實(shí)驗(yàn)所用PC機(jī)信息如下：Windows 7系統(tǒng)，CPU為Intel Core i5-2430 2.40 Hz，內(nèi)存為4 G。

3.2.1 Salinas高光譜數(shù)據(jù)

數(shù)據(jù)一為覆蓋美國(guó)加利福尼亞州薩利納斯(Salinas)地區(qū)的AVIRIS數(shù)據(jù)，數(shù)據(jù)大小為512×217，空間分辨率為3.7 m。AVIRIS共有240個(gè)波段，其中波段108～112，154～167，224共20個(gè)波段為水吸收波段，在進(jìn)行聚類處理之前需要去掉這20個(gè)波段，剩余204個(gè)波段。該高光譜圖像包含的地物主要為植被和土壤，一共包含16種不同的地物，其中具體的每種地物樣本點(diǎn)數(shù)量如表1所示，該圖像地物對(duì)應(yīng)的真值圖如圖3所示。

表1 Salinas各地物樣本信息Table 1 Detailed information of various samples in Salinas data

圖3 Salinas地物真值圖Fig.3 Ground truth of Salinas data

為便于比較，對(duì)于K-means、KSC和MNF+KSC方法，我們?cè)O(shè)置最終的分類類別數(shù)均為16。對(duì)于KSC和MNF+KSC方法，首先利用K-means對(duì)該數(shù)據(jù)進(jìn)行粗聚類(人工設(shè)置類別數(shù)分別為600和300)，最終聚類結(jié)果如圖4所示，分別計(jì)算當(dāng)粗聚類類別數(shù)取600和300時(shí)，聚類結(jié)果的總體分類精度和Kappa系數(shù)，如表2所示。

Salinas數(shù)據(jù)的特點(diǎn)是包含的地物主要為植被和裸地，而它們各自間光譜相似性很高。圖5分別畫出3種不同植被和不同裸地的光譜曲線，從中可以看出它們各自之間光譜曲線非常相似，所以Salinas數(shù)據(jù)所包含的地物相對(duì)來(lái)說(shuō)較難區(qū)分。對(duì)于具有這樣特點(diǎn)的數(shù)據(jù)，當(dāng)粗聚類類別參數(shù)C分別取600和300時(shí)，從表2可以看出KSC與MNF+KSC相對(duì)K-Means算法來(lái)說(shuō)，聚類的總體分類精度及Kappa系數(shù)均有所提升。而與KSC的計(jì)算時(shí)間相比，MNF+KSC的計(jì)算效率和聚類精度均有所提高。

3.2.2 Pavia Center高光譜數(shù)據(jù)

數(shù)據(jù)二為覆蓋意大利北部帕維亞中心(Pavia Center)地區(qū)的ROSIS數(shù)據(jù)，去掉其中的無(wú)效數(shù)據(jù)后圖像大小為1 096×715，幾何分辨率為1.3 m，共102個(gè)波段。該圖像包含的地物有植被、建筑和水體等共9類地物，具體每種植被樣本點(diǎn)數(shù)量如表3所示，該圖像地物對(duì)應(yīng)的真值圖如圖6所示。

圖4 Salinas數(shù)據(jù)實(shí)驗(yàn)結(jié)果對(duì)比圖(粗聚類類別數(shù)為600)Fig.4 Comparison of clustering results in Salinas data (coarse clustering number is 600)

方法粗聚類類別數(shù)為600粗聚類類別數(shù)為300總體分類精度/%Kappa系數(shù)計(jì)算時(shí)間/s總體分類精度/%Kappa系數(shù)計(jì)算時(shí)間/sK-means62.990.589661.949262.990.589661.9 492KSC69.690.6327339.792067.700.6348210.6877MNF+KSC81.210.7880184.509378.320.755394.8923

圖5 Salinas數(shù)據(jù)中不同地物的光譜曲線Fig.5 Spectra of different objects in Salinas data

圖6 Pavia Center地物真值圖Fig.6 Ground truth of Pavia Center data

為便于比較，對(duì)于K-means算法、KSC算法和MNF+KSC方法，將最終的分類類別數(shù)目均設(shè)置為9。對(duì)于KSC和MNF+KSC方法，首先利用K-

表3 Pavia Center各地物樣本信息Table 3 Detailed information of various samples in Pavia Center data

means方法對(duì)該數(shù)據(jù)進(jìn)行粗聚類(人工設(shè)置類別數(shù)分別為600和300)，最終的聚類結(jié)果圖如圖7所示。計(jì)算當(dāng)粗聚類類別數(shù)分別取600和300時(shí)，聚類結(jié)果的總體分類精度和Kappa系數(shù)，如表4所示。

圖7 Pavia Center數(shù)據(jù)實(shí)驗(yàn)結(jié)果對(duì)比圖(粗聚類類別數(shù)為600)Fig.7 Comparison of clustering results in Pavia Center data (coarse clustering number is 600)

方法粗聚類類別數(shù)為600粗聚類類別數(shù)為300總體分類精度/%Kappa系數(shù)計(jì)算時(shí)間/s總體分類精度/%Kappa系數(shù)計(jì)算時(shí)間/sK-Means66.680.5755419.147666.680.5755419.1476KSC71.480.61142544.353172.860.63051517.0328MNF+KSC92.550.89971820.698890.410.8694961.2714

Pavia Center數(shù)據(jù)不同于Salinas數(shù)據(jù)，其特點(diǎn)是包含的不同地物之間光譜相似性都很低。圖8分別畫出3種地物的光譜曲線，從中可以看出它們之間的光譜曲線差異較大，所以Pavia Center數(shù)據(jù)所包含的地物相對(duì)來(lái)說(shuō)較易區(qū)分。對(duì)具有這樣特點(diǎn)的數(shù)據(jù)，當(dāng)粗聚類類別參數(shù)C分別取600和300時(shí)，從表4可以看出KSC與MNF+KSC相對(duì)K-means算法來(lái)說(shuō)，聚類結(jié)果的總體分類精度及Kappa系數(shù)均有較大幅度提升。而與KSC相比，MNF+KSC的計(jì)算效率和聚類精度均有所提高。

圖8 Pavia Center數(shù)據(jù)不同地物光譜曲線Fig.8 Spectra of different objects in Pavia Center data

由以上實(shí)驗(yàn)可知，KSC方法與K-means算法比較，聚類精度有所提升?？紤]高光譜圖像冗余性和信噪比后，采用的MNF+KSC的方法，在降低計(jì)算復(fù)雜度的同時(shí)，聚類精度又在KSC的基礎(chǔ)上有較大幅度提升。

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

本文針對(duì)譜聚類方法在眾多聚類方法中有明顯優(yōu)勢(shì)但是無(wú)法在高光譜圖像聚類中直接應(yīng)用的問(wèn)題，提出結(jié)合K-means的KSC方法，采用K-means首先對(duì)高光譜數(shù)據(jù)進(jìn)行粗聚類，然后再利用譜聚類方法對(duì)上述粗聚類中心做進(jìn)一步聚類處理，如此使得譜聚類方法在高光譜圖像中得以應(yīng)用。為進(jìn)一步降低計(jì)算復(fù)雜度和提高數(shù)據(jù)信噪比，采用MNF方法對(duì)高光譜數(shù)據(jù)進(jìn)行特征降維，然后再進(jìn)行KSC聚類處理。通過(guò)對(duì)真實(shí)高光譜圖像數(shù)據(jù)進(jìn)行實(shí)驗(yàn)，驗(yàn)證了KSC方法和MNF+KSC方法比K-means方法對(duì)高光譜圖像有更好的聚類性能。