張亞平，張 宇，楊 楠,2，羅 曉，羅 謙

(1. 哈爾濱工業(yè)大學(xué)交通科學(xué)與工程學(xué)院，黑龍江 哈爾濱 150090; 2. 國土資源部城市土地資源監(jiān)測與仿真重點實驗室，廣東 深圳 518034; 3. 中國民用航空局第二研究所，四川 成都 610041)

譜聚類算法作為聚類分析中的一個全新分支，無需對數(shù)據(jù)的全局結(jié)構(gòu)作出假設(shè)，適用于任意空間分布數(shù)據(jù)的聚類，具有識別非凸分布數(shù)據(jù)、高效聚類的能力，適合于解決諸多實際問題，因而在短短的幾年時間里就引起了中外學(xué)術(shù)界的廣泛關(guān)注。近年來，國內(nèi)外學(xué)者在譜聚類算法理論和應(yīng)用研究方面展開了大量研究工作，取得了諸多理論和應(yīng)用方面的研究成果，促進了譜聚類算法體系及其應(yīng)用技術(shù)的發(fā)展。目前，譜聚類算法已經(jīng)成功應(yīng)用于人臉識別[1]、圖像分割[2]、醫(yī)學(xué)圖像分析、信息檢索[3]、電力系統(tǒng)建模[4]、蛋白質(zhì)數(shù)據(jù)分析[5]等領(lǐng)域。

本文的待分類圖像為高光譜圖像，文獻[6]對高光譜圖像的分類算法進行了系統(tǒng)性的介紹。當前，關(guān)于遙感圖像的機器學(xué)習(xí)算法的研究很多，如文獻[7]利用SVM研究了幾何特征、圖論特征對土地利用分類的影響，文獻[8]利用卷積神經(jīng)網(wǎng)絡(luò)對地表覆蓋類型的分類精度進行評價，文獻[9]將FSVM算法與ISODATA算法相結(jié)合，新的算法更適用于高分辨率遙感圖像，分類精度也得到了較大提高。但是關(guān)于遙感圖像的譜聚類算法的研究卻很少。本文將非監(jiān)督分類算法在機器學(xué)習(xí)領(lǐng)域中的研究成果引入遙感圖像數(shù)據(jù)處理領(lǐng)域，嘗試應(yīng)用集成在Sklearn模塊中的快速解求大規(guī)模矩陣端元奇異值的Lanczos算法[10]來求解拉普拉斯矩陣的少數(shù)最小特征值及其對應(yīng)的特征向量，以解決譜聚類算法計算效率過低的問題，實現(xiàn)能夠區(qū)分遙感圖像像素空間分布形態(tài)復(fù)雜的點群的非監(jiān)督分類算法；然后與傳統(tǒng)的K-均值算法數(shù)據(jù)進行對比，發(fā)現(xiàn)譜聚類算法易于識別線性地物，說明其應(yīng)用于遙感圖像分類的可行性。

1 譜聚類算法

聚類算法主要包含構(gòu)成圖和裁剪圖兩個步驟。先由眾多數(shù)據(jù)點構(gòu)成一張圖(graph)，然后再按照一定的準則進行切圖(graph)。主要的切圖規(guī)則有最小割集規(guī)則(minimum cut)、平均割集規(guī)則(average cut)、規(guī)范化割集規(guī)則(normalized cut)、比例割集規(guī)則(ratio cut)等[11]。

1.1 最小割集規(guī)則

最小割集規(guī)則較適宜于將譜圖G分割為A、B兩個子圖的情況。A、B兩個子圖合并起來即為G，且A、B之間沒有重疊。即A∪B=G，A∩B=?。定義一個不同邊權(quán)重之和的損失函數(shù)為

(1)

最小割集規(guī)則為cut(A,B)等于最小值時的分割規(guī)則。使用最小割集規(guī)則對譜圖進行分割可以取到很好的聚類結(jié)果，但是這種分割規(guī)則只考慮了兩個子圖之間的權(quán)值，沒有考慮子圖內(nèi)部的權(quán)值問題，致使較易發(fā)生譜圖歪斜偏向小區(qū)域分割的現(xiàn)象。因此，學(xué)者們還提出了規(guī)范割集規(guī)則與比例割集規(guī)則，以防止譜圖歪斜分割的問題。

1.2 規(guī)范割集規(guī)則

定義一個目標函數(shù)，目標函數(shù)為Ncut(A,B)，當致使目標函數(shù)Ncut(A,B)最小的時候的一種譜圖分割規(guī)則稱為規(guī)范割集規(guī)則。目標函數(shù)Ncut(A,B)為

(2)

該規(guī)則不僅可以量度簇內(nèi)樣本間的相似程度,還可以量度簇間樣本間的相似程度。

1.3 平均割集規(guī)則

平均割集規(guī)則目標系數(shù)Avcut(A,B)為

(3)

式中,|A|、|B|分別代表給A、B子圖各自的頂點數(shù)目。該規(guī)則的目標函數(shù)表示A、B子圖各自與損失函數(shù)的比值之和。

1.4 比率割集規(guī)則

比率割集規(guī)則目標函數(shù)Rcut(A,B)為

(4)

式中，|A|表示A子圖的頂點個數(shù)；|B|表示B子圖中的頂點個數(shù)。當目標函數(shù)最小時，簇間樣本數(shù)據(jù)的相似性最小。

1.5 譜聚類算法圖像矩陣

譜聚類算法的流程如圖1所示，可以劃分為以下幾個主要步驟：①相似度矩陣計算；②度矩陣的計算；③拉普拉斯矩陣的計算；④特征值與特征向量的計算；⑤K-均值聚類。主要包括相似度矩陣、度矩陣、拉普拉斯矩陣3個重要矩陣。

1.5.1 相似矩陣

相似矩陣A為一個對稱矩陣，它的每個元素是由不同樣本之間的相似度組成的，相似度一般由樣本之間的距離來度量，也有用高斯核函數(shù)與余弦相似度來度量的。如果用距離來度量樣本相似度，通常矩陣中只保留距離小于給定閾值的相似度值，此時的矩陣A稱為K-鄰域矩陣，是一個稀疏對稱矩陣。有許多表示樣本點之間距離的方式，最常見的是歐氏距離。歐氏距離的公式為

(5)

式中,向量xi和xl(i,l=1,2,…,n)為遙感圖像中的兩個像元；n代表遙感圖像像元總數(shù)。

當使用高斯核函數(shù)計算相似矩陣A時，核函數(shù)中的參數(shù)σ代表核函數(shù)的響應(yīng)寬度，與K-鄰域矩陣中的閾值K具有相同的作用，此時的相似矩陣A實質(zhì)上稀疏核矩陣。給定遙感圖像中的兩個像素向量xi和xl(i,l=1,2,…,n)，n代表遙感圖像像素總數(shù)，則高斯核函數(shù)可表示為[12]

(6)

1.5.2 度矩陣

度矩陣D(degree matrix)是一個對角矩陣，是在矩陣A的基礎(chǔ)上得到的派生矩陣。其對角線元素為矩陣A中相應(yīng)的行或列的所有元素之和。矩陣D可以表示為[13]

(7)

1.5.3 拉普拉斯矩陣

拉普拉斯矩陣L也稱為導(dǎo)納矩陣、基爾霍夫矩陣，是圖論中一個圖(graph)直觀的矩陣表示。拉普拉斯矩陣主要有以下幾個特征：①拉普拉斯矩陣為對稱陣，且特征值非負，即為半正定矩陣；②拉普拉斯矩陣的零特征值的個數(shù)即為圖連通區(qū)域的個數(shù)；③拉普拉斯矩陣最小的非零特征值與圖的代數(shù)連通度相等；④拉普拉斯矩陣每一行之和均為零。由相似矩陣A和度矩陣D得到拉普拉斯矩陣的數(shù)學(xué)基礎(chǔ)公式為[14]

L=D-A

(8)

譜聚類算法的核心部分是解求拉普拉斯矩陣L的特征值和特征向量，而一幅遙感圖像的像素數(shù)目一般在數(shù)萬以上，如何快速解求這一超高階拉普拉斯矩陣的特征值與特征向量問題，從而使譜聚類算法能夠應(yīng)用于遙感圖像的非監(jiān)督分類研究，是需要妥善解決的關(guān)鍵科學(xué)問題之一。這一尚待解決的關(guān)鍵科學(xué)問題的存在，也正是導(dǎo)致目前遙感圖像處理領(lǐng)域尚無研究者應(yīng)用譜聚類算法進行遙感圖像非監(jiān)督分類研究的原因之一。

2 高光譜圖像譜聚類算法框架

高光譜圖像的數(shù)據(jù)量大，常含有上百個波段。因此對算法的運算是一個挑戰(zhàn)，在面向高光譜圖像的譜聚類算法中，為了提升算法的運算速度集成了Lanczos算法，使得譜聚類算法更加快速。本文提出的高光譜圖像譜聚類算法框架如圖2所示。傳統(tǒng)意義上的譜聚類算法主要分為非正則化譜聚類算法(Unnormalized spectral clustering)和正則化譜聚類算法(Unnormalized spectral clustering)。正則化譜聚類算法與非正則化譜聚類算法框架如下[2]。

正則化譜聚類算法框架：

輸入：相似矩陣A∈Rn×n,聚類數(shù)目K。

(1) 由鄰接權(quán)重矩陣構(gòu)成相似圖。

(2) 計算非正則化拉普拉斯矩陣L。

(3) 由拉普拉斯矩陣計算K個特征向量u1,u2,…,uk。

(4) 由特征向量u1,u2,…,uk組成特征矩陣U∈Rn×k的每一列，得到特征矩陣。

(5) 也可由yi∈Rk(i=1,2，…,n)組成特征矩陣的每一行，得到特征矩陣。

(6) 用K-均值算法將Rk空間中的數(shù)據(jù)劃分為C1,C2,…,Ck個簇。

輸出：A1，A2，…，Ak個類別。

非正則化譜聚類算法框架：

輸入：相似矩陣A∈Rn×n,聚類數(shù)目K。

(1) 由鄰接權(quán)重矩陣構(gòu)成相似圖。

(2) 計算非正則化拉普拉斯矩陣L。

(3) 由廣義特征值公式Lu=λDu計算K個特征向量u1,u2,…,uk。

(4) 由特征向量u1,u2,…,uk組成特征矩陣U∈Rn×k的每一列，得到特征矩陣。

(5) 也可由yi∈Rk(i=1,2,…,n)組成特征矩陣的每一行，得到特征矩陣。

(6) 用K-均值算法將Rk空間中的數(shù)據(jù)劃分為C1,C2,…,Ck個簇。

輸出：A1,A2,…,Ak個類別。

遙感圖像譜聚類算法是在Python集成開發(fā)環(huán)境(Python-IDEL)中，通過無縫集成Pyhthon-GDAL函數(shù)庫和sklearn開源代碼庫實現(xiàn)的。遙感圖像譜聚類算法需依次完成下列操作過程：遙感圖像的輸入操作，遙感圖像的數(shù)據(jù)預(yù)處理，K-鄰域矩陣或高斯核矩陣的計算，拉普拉斯矩陣構(gòu)造，拉普拉斯矩陣特征值和特征向量求解，K-均值聚類算法輸入數(shù)據(jù)構(gòu)造(用拉普拉斯矩陣的數(shù)個最小特征值對應(yīng)的特征向量構(gòu)造新的待分類數(shù)據(jù))，待分類數(shù)據(jù)的K-均值聚類分析，分類結(jié)果圖像的輸出，其算法框架可由圖2表示。

由于遙感圖像的類型較多，數(shù)據(jù)觀測尺度和記錄方式不盡相同。如多光譜遙感圖像像元每個波段的DN值(Digital Number)為0～255的整數(shù)值，在進行譜聚類之前，可將DN值轉(zhuǎn)化成實數(shù)值，公式為

(9)

本文待分類的遙感數(shù)據(jù)為高光譜遙感數(shù)據(jù)，盡管各波段的光譜值為實數(shù)值，但是不同波段的光譜數(shù)據(jù)量綱可能不同，因此，需在分類前進行數(shù)據(jù)預(yù)處理以統(tǒng)一高光譜數(shù)據(jù)各波段的數(shù)據(jù)量綱。可運用數(shù)據(jù)正規(guī)化變換方法對高光譜數(shù)據(jù)進行預(yù)處理，公式為

(10)

式中,xmaxj和xminj分別為高光譜遙感圖像第j波段光譜觀測值的最大和最小值。正規(guī)化變換后的高光譜遙感圖像各波段光譜取值為0～1之間的實數(shù)。

預(yù)處理后的遙感圖像數(shù)據(jù)可以直接應(yīng)用式(5)計算K-鄰域矩陣或應(yīng)用式(6)計算高斯核函數(shù)矩陣。在計算K-鄰域矩陣時，距離閾值需要根據(jù)遙感圖像數(shù)據(jù)的具體特征來確定一個適當值，如果該閾值過大，將導(dǎo)致相似矩陣A的非零元素過多，導(dǎo)致A不滿足稀疏性；如果距離閾值過小，會丟失一些必要的像素相似性信息，使矩陣A不能正確反映遙感圖像數(shù)據(jù)類群結(jié)構(gòu)信息。在應(yīng)用高斯核函數(shù)計算高斯核矩陣時，核函數(shù)響應(yīng)寬度需要根據(jù)遙感圖像的特征選取，過大或過小的響應(yīng)寬度都會影響遙感圖像分類效果。

計算出相似矩陣A之后，可以很容易地在矩陣A基礎(chǔ)上構(gòu)造拉普拉斯矩陣L，將矩陣A每一行(或列)全部元素求和，即可得到矩陣L的對角線元素。然后，將矩陣A的非對角線元素取相反數(shù)，即可獲得矩陣L的非對角線元素。此時，拉普拉斯矩陣的構(gòu)造過程結(jié)束。

當拉普拉斯矩陣L滿足稀疏性時，可以用于Lanczos算法快速求解矩陣L的少數(shù)最大(或最小)特征值及其對應(yīng)的特征向量。在譜聚類算法中，需要求出矩陣L的前G(G<

3 試驗驗證

為了驗證算法的可行性，采用美國圣地亞哥市機場400×400×189的高光譜數(shù)據(jù)立方體(AVIRIS)進行分類試驗，為使遙感圖像的像元灰度值更加真實地反映地物的反射率，需要對遙感圖像做大氣校正。該研究區(qū)黑暗像元較多，因此采用忽略大氣散射作用與相鄰像元漫反射作用的簡化黑暗像元大氣校正法，采取的確定黑暗像素值的方法為波段最小值法。圖3左、中、右分別為研究區(qū)域原圖、譜聚類算法分類結(jié)果和K-均值算法分類結(jié)果。

由圖3的譜聚類算法分類結(jié)果和K-均值算法分類結(jié)果對比分析可得，在譜聚類算法分類結(jié)果中停機坪上的兩架飛機被很好地識別出來，對道路的區(qū)分也很好，但是在K-均值算法分類結(jié)果中卻將飛機與地面劃分為一類。對于圖中右下角4棟房屋的頂面，因為日照角度的原因?qū)е乱粋?cè)上有陰影，反映在圖像上房屋兩側(cè)的像元灰度值有差異，譜聚類算法據(jù)此將同一房屋的兩側(cè)分為不同的類，K-均值算法分類結(jié)果則并未出現(xiàn)這樣的現(xiàn)象。由此可見，譜聚類算法易于識別K-均值算法不易識別的地物類別，即譜聚類算法對像元灰度值差異比較敏感。對于線性地物、兩條不同材質(zhì)道路的分界線譜聚類算法能較好地識別，而K-均值算法往往對道路的邊界線(即灰度值突變處)反應(yīng)不敏感。譜聚類算法的時間復(fù)雜度比K-均值算法高，因此運行速度比較慢，經(jīng)過與Lanczos算法集成之后的譜聚類算法在運行速度上已經(jīng)得到了比較大的提升，通過對譜聚類算法和K-均值算法分類,數(shù)目均設(shè)定為5個。兩次的試驗研究表明：K-均值算法的運算時間為45.9 s，譜聚類算法的時間為40多分鐘，雖然集成Lanczos的譜聚類算法的運行效率仍低于K-均值算法，但卻可以保證更高的分類精度。

4 結(jié) 語

本文提出了一種基于高光譜遙感圖像的譜聚類算法，應(yīng)用快速解求大規(guī)模矩陣端元奇異值的Lanczos算法求解拉普拉斯矩陣的最小特征值及其對應(yīng)的特征向量，以此提高譜聚類算法的運算速度，初步解決了將譜聚類算法應(yīng)用于遙感圖像處理中運算速度慢的問題，并發(fā)現(xiàn)譜聚類算法易于識別線性地物，驗證了譜聚類算法應(yīng)用于遙感圖像分類的可行性。