一種高光譜遙感圖像快速譜聚類算法

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

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

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

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

1 譜聚類算法

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

1.1 最小割集規則

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

(1)

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

1.2 規范割集規則

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

(2)

該規則不僅可以量度簇內樣本間的相似程度,還可以量度簇間樣本間的相似程度。

1.3 平均割集規則

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

(3)

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

1.4 比率割集規則

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

(4)

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

1.5 譜聚類算法圖像矩陣

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

1.5.1 相似矩陣

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

(5)

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

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

(6)

1.5.2 度矩陣

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

(7)

1.5.3 拉普拉斯矩陣

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

L=D-A

(8)

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

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

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

正則化譜聚類算法框架：

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

(1) 由鄰接權重矩陣構成相似圖。

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

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

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

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

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

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

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

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

(1) 由鄰接權重矩陣構成相似圖。

(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空間中的數據劃分為C1,C2,…,Ck個簇。

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

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

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

(9)

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

(10)

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

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

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

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

3 試驗驗證

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

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

4 結 語

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