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

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

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

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

1 譜聚類算法

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

1.1 最小割集規(guī)則

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

(1)

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

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

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

(2)

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

1.3 平均割集規(guī)則

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

(3)

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

1.4 比率割集規(guī)則

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

(4)

式中，|A|表示A子圖的頂點(diǎn)個(gè)數(shù)；|B|表示B子圖中的頂點(diǎn)個(gè)數(shù)。當(dāng)目標(biāo)函數(shù)最小時(shí)，簇間樣本數(shù)據(jù)的相似性最小。

1.5 譜聚類算法圖像矩陣

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

1.5.1 相似矩陣

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

(5)

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

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

(6)

1.5.2 度矩陣

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

(7)

1.5.3 拉普拉斯矩陣

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

L=D-A

(8)

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

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

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

正則化譜聚類算法框架：

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

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

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

(3) 由拉普拉斯矩陣計(jì)算K個(gè)特征向量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個(gè)簇。

輸出：A1，A2，…，Ak個(gè)類別。

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

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

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

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

(3) 由廣義特征值公式Lu=λDu計(jì)算K個(gè)特征向量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個(gè)簇。

輸出：A1,A2,…,Ak個(gè)類別。

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

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

(9)

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

(10)

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

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

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

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

3 試驗(yàn)驗(yàn)證

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

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

4 結(jié) 語

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