吳翔，王鳳艷，林楠，王明常

1.吉林大學(xué) 地球探測(cè)科學(xué)與技術(shù)學(xué)院，長(zhǎng)春 130026；2.吉林建筑大學(xué) 測(cè)繪與勘查工程學(xué)院，長(zhǎng)春 130018

0 引言

三維激光掃描技術(shù)是通過發(fā)射激光來獲取物體表面的三維空間信息和物理信息，實(shí)現(xiàn)了測(cè)繪領(lǐng)域的又一次技術(shù)革新[1]。與一些傳統(tǒng)的全站儀測(cè)量技術(shù)、GPS測(cè)量技術(shù)相比，三維激光掃描技術(shù)有著快速、準(zhǔn)確和非接觸性等特點(diǎn)，廣泛應(yīng)用于城市建設(shè)、工程設(shè)計(jì)和地貌分析等領(lǐng)域[2]。但三維激光掃描技術(shù)獲取的點(diǎn)云數(shù)據(jù)量大且雜亂無序，不利于后續(xù)的三維場(chǎng)景識(shí)別和地物分類等分析與應(yīng)用。因此如何有效、快速和準(zhǔn)確的對(duì)點(diǎn)云數(shù)據(jù)進(jìn)行分類成了眾多學(xué)者研究的重點(diǎn)。近年來，針對(duì)三維激光點(diǎn)云數(shù)據(jù)的分類問題，國(guó)內(nèi)外眾多學(xué)者展開了相關(guān)研究。例如：李真等[3]利用K-means算法對(duì)不同種類的樹進(jìn)行試驗(yàn)，通過計(jì)算高度和強(qiáng)度直方圖確定K值和初始聚類中心點(diǎn)來做聚類分析，實(shí)現(xiàn)了單株樹木原始點(diǎn)云數(shù)據(jù)分類；胡健等[4]利用地面三維激光掃描儀采集雪地的點(diǎn)云數(shù)據(jù)，采用相關(guān)分析獲取點(diǎn)云回波強(qiáng)度與入射角度、距離之間的修正函數(shù)，最終實(shí)現(xiàn)不同類型雪及地物的分類；Zhang et al.[5]基于自然指數(shù)函數(shù)閾值將點(diǎn)云劃分為不同大小的層次簇，并對(duì)其形狀特征進(jìn)行提取和編碼，實(shí)現(xiàn)了城市環(huán)境中的地物分類；Lai和Fox[6]先利用Meanshift聚類方法對(duì)環(huán)境點(diǎn)云進(jìn)行分割，然后再利用構(gòu)成物體點(diǎn)云的宏觀特征將各物體分類；王書民等[7]提出了結(jié)合點(diǎn)云數(shù)據(jù)的波形信息，利用模糊C均值的分類方法，將研究區(qū)域中樹木、地面和建筑進(jìn)行分類；袁夏等[8]提出利用改進(jìn)的高斯混合模型對(duì)地面不同地物進(jìn)行分類；倪明等[9]以地面拍攝影像為輔助條件，結(jié)合高程、紋理等信息進(jìn)行地面點(diǎn)云分類。根據(jù)上述研究可以看出，學(xué)者們主要利用K-means算法、Meanshift算法等傳統(tǒng)聚類算法對(duì)點(diǎn)云進(jìn)行分類。但是這些算法都需建立在樣本空間為凸的條件下進(jìn)行分類，當(dāng)樣本空間不為凸時(shí)會(huì)陷入局部最優(yōu)，導(dǎo)致分類效果不佳[10]。因此，為了能在任意的樣本空間中進(jìn)行聚類且收斂于最優(yōu)解，學(xué)者們開始采用譜聚類算法進(jìn)行聚類[11]。譜聚類算法的核心思想來源于譜圖理論，由于其具有在任意形狀的樣本空間上進(jìn)行聚類并收斂于全局最優(yōu)解等傳統(tǒng)聚類算法難以實(shí)現(xiàn)的特點(diǎn)，因此被廣泛應(yīng)用于人臉識(shí)別、文本分類和計(jì)算機(jī)視覺等多個(gè)領(lǐng)域[12-14]。

筆者基于Python平臺(tái)，以三維激光點(diǎn)云數(shù)據(jù)為數(shù)據(jù)源，利用ARCGIS軟件將點(diǎn)云數(shù)據(jù)轉(zhuǎn)換為柵格數(shù)據(jù)，并結(jié)合紋理以及形狀信息構(gòu)建譜聚類算法模型，對(duì)實(shí)驗(yàn)區(qū)內(nèi)的典型地物進(jìn)行分類，與傳統(tǒng)的K-means算法和高斯混合模型進(jìn)行精度對(duì)比分析，探討譜聚類算法用于典型地物的分類可行性。

1 方法

1.1 譜聚類算法原理

譜聚類算法是聚類分析中的一種新型聚類算法，相較于傳統(tǒng)的聚類算法，它建立在譜圖理論基礎(chǔ)之上[15]。其基本思想：首先利用高斯核函數(shù)計(jì)算數(shù)據(jù)之間的歐式距離構(gòu)建相似度矩陣，將聚類問題轉(zhuǎn)換為圖分割問題；之后根據(jù)相似度矩陣來獲取拉普拉斯矩陣，并利用一種連續(xù)的松弛形式將圖分割問題轉(zhuǎn)化為拉普拉斯矩陣的譜分解問題，即計(jì)算拉普拉斯矩陣前K個(gè)最小特征值對(duì)應(yīng)的特征向量；最后利用經(jīng)典的聚類算法(如K-means算法)對(duì)其進(jìn)行聚類。

(1)相似度矩陣、度矩陣、拉普拉斯矩陣

譜聚類算法的相似度矩陣是通過數(shù)據(jù)之間的相似度來構(gòu)建的，其步驟：根據(jù)給定的數(shù)據(jù)集X，X={x1,x2,…,xn}，將每個(gè)數(shù)據(jù)點(diǎn)xi看成無向圖G=(V,E)的一個(gè)頂點(diǎn)vi，并根據(jù)vi和vj連成的邊Eij賦權(quán)值wij，可以得到無向加權(quán)圖G(V,E,W)。通過該方法可以把數(shù)據(jù)集X轉(zhuǎn)化為一個(gè)帶有權(quán)值的無向圖G(V,E,W)，并得到相似度矩陣W。相似度是利用高斯核函數(shù)來計(jì)算的，公式為：

(1)

式中：xi表示第i個(gè)點(diǎn);‖xi-xj‖2表示為xi、xj之間的歐式距離;σ為高斯核參數(shù)，控制樣本點(diǎn)的鄰域?qū)挾?，σ越大表示樣本點(diǎn)與距離較遠(yuǎn)的樣本點(diǎn)的相似度越大，反之亦然。

將相似度矩陣W的每行元素相加可以得到度矩陣D，公式為：

(2)

式中：D為di組成的n×n對(duì)角矩陣，D=diag(d1,d2, …,dn)。

雖然相似度矩陣包含了聚類所需要的數(shù)據(jù)信息，然而在實(shí)際過程中，圖分割的過程都是通過求解拉普拉斯矩陣的特征值與特征向量來完成的。拉普拉斯矩陣L可用下式求得：

L=D-W

(3)

(2)圖的分割

通過上述分析，利用高斯核函數(shù)構(gòu)造的相似度矩陣完成了從聚類到圖分割的過程。那么，怎么對(duì)圖進(jìn)行最優(yōu)分割是譜聚類算法的一個(gè)難點(diǎn)。學(xué)者們研究出比較好的處理方法是將圖分割問題轉(zhuǎn)化為求解連續(xù)的拉普拉斯矩陣特征向量問題。即根據(jù)拉普拉斯矩陣求解的前K個(gè)最小特征值對(duì)應(yīng)的特征向量，可以構(gòu)成相應(yīng)的向量空間，其中不同的向量空間對(duì)應(yīng)著不同的分割，即可以達(dá)到圖分割的目的。目前常見的圖分割準(zhǔn)則主要有最小割集準(zhǔn)則、比例割集準(zhǔn)則、規(guī)范割集準(zhǔn)則和平均割集準(zhǔn)則等。本文所使用的圖分割準(zhǔn)則為比例割集準(zhǔn)則，原理：

令圖G=(V，E)，圖G可分為A、B兩部分，并有A∩B=φ,A∪B=V，兩個(gè)節(jié)點(diǎn)間邊的權(quán)值為W，那么圖G劃分為A、B兩部分的代價(jià)函數(shù)為：

(4)

其中|A|、|B|表示子圖A、B中節(jié)點(diǎn)的個(gè)數(shù)。

(3)算法步驟

①根據(jù)給定的數(shù)據(jù)集X，X={x1,x2,…,xm}，利用高斯核函數(shù)計(jì)算相似度矩陣W，并構(gòu)造拉普拉斯矩陣L；

②求取拉普拉斯矩陣L的前K個(gè)最小特征值，并計(jì)算對(duì)應(yīng)的特征向量u1,u2,…,uk，構(gòu)成矩陣U={u1,u2,…,uk}；

③取yi為矩陣U的第i行向量，并利用其他經(jīng)典聚類算法(如K-means算法)，將新樣本點(diǎn)Y={y1,y2,…,yn}聚類成簇C1,C2,…,Ck。

1.2 譜聚類算法建模過程

筆者利用譜聚類算法處理高維空間數(shù)據(jù)的優(yōu)勢(shì)，將其應(yīng)用到三維激光點(diǎn)云數(shù)據(jù)中，構(gòu)建基于譜聚類算法的點(diǎn)云數(shù)據(jù)分類方法。具體建模過程：利用sklearn模塊(scikit-learn是python中的第三方庫(kù)，庫(kù)中包括分類、回歸、聚類和降維等相關(guān)算法)中的譜聚類代碼集成到遙感影像譜聚類算法中，分別將僅含反射強(qiáng)度特征和加入了紋理、形狀信息的多源特征建立譜聚類模型，利用混淆矩陣來評(píng)價(jià)譜聚類算法模型的分類結(jié)果，并根據(jù)總體分類精度和Kappa系數(shù)來判斷其分類效果(圖1)。

2 數(shù)據(jù)采集與處理

2.1 數(shù)據(jù)采集

以吉林省長(zhǎng)春市某一場(chǎng)地為實(shí)驗(yàn)區(qū)，實(shí)驗(yàn)區(qū)面積為100 m×100 m，利用Z+F IMAGER 5010C掃描儀對(duì)實(shí)驗(yàn)區(qū)進(jìn)行數(shù)據(jù)采集，采集方式為多個(gè)測(cè)站進(jìn)行獨(dú)立掃描，共掃描了14站，實(shí)驗(yàn)區(qū)包括建筑用地、林地、裸地和草地等地物(圖2)。

2.2 特征分析

利用獲取的實(shí)驗(yàn)區(qū)點(diǎn)云數(shù)據(jù)，導(dǎo)出地物的反射強(qiáng)度值，并隨機(jī)選取同一地物樣本各20個(gè)(表1)，分析不同地物反射強(qiáng)度值差異性。從圖3中可以看出，不同地物的反射強(qiáng)度值不同且其反射強(qiáng)度值都在一定范圍上下波動(dòng)，驗(yàn)證了反射強(qiáng)度可以作為分類信息將草地、建筑用地、林地和裸地等不同地物區(qū)分。

2.3 點(diǎn)云柵格化

在特征分析過程中得出了反射強(qiáng)度可以作為分類信息區(qū)分草地、建筑用地、林地和裸地等不同地物，故采用點(diǎn)云數(shù)據(jù)中的反射強(qiáng)度信息轉(zhuǎn)換為柵格數(shù)據(jù)。因此利用柵格數(shù)據(jù)來做后續(xù)研究，將點(diǎn)云數(shù)據(jù)典型地物的分類問題轉(zhuǎn)換為柵格數(shù)據(jù)的分類問題，利用相關(guān)算法對(duì)柵格影像進(jìn)行分類研究。圖4為經(jīng)過拼接、去噪后的點(diǎn)云數(shù)據(jù)，圖5為轉(zhuǎn)換后的柵格數(shù)據(jù)。

圖2 實(shí)驗(yàn)區(qū)示意圖Fig.2 Experimental area diagram

表1 典型地物反射強(qiáng)度值

2.4 特征信息提取

(1)紋理信息提取

在利用柵格影像進(jìn)行分類時(shí)，由于物體的邊緣呈現(xiàn)灰度的不連續(xù)性和混合像元的存在，僅依靠反射強(qiáng)度信息通常存在丟失像元的現(xiàn)象，而將紋理和反射強(qiáng)度信息結(jié)合起來進(jìn)行分析可以更好地體現(xiàn)地物的宏觀性質(zhì)和細(xì)微信息。本文采用統(tǒng)計(jì)方法中的灰度共生矩陣(gray level co-occurrence matrix, GLCM) 方式進(jìn)行分析。GLCM的紋理分析方法由西班牙數(shù)學(xué)家Haralick于1973年首先提出[16]，GLCM為方陣，維數(shù)即灰度級(jí)數(shù)，通過計(jì)算各灰度級(jí)之間的聯(lián)合條件概率密度P(i,j|d,θ)來表示紋理，公式：

圖3 實(shí)驗(yàn)區(qū)地物反射強(qiáng)度折線圖Fig.3 Line graph of ground object reflectance in experimental area

圖4 點(diǎn)云數(shù)據(jù)Fig.4 Point cloud data

圖5 柵格數(shù)據(jù)Fig.5 Raster data

P(i,j|d,θ)=φ{(diào)(x,y)|f(x,y)=i,f(x+dx,y+dy)=j}

x,y=0, 1, …,N-1

(5)

本文選用的是GLCM紋理分析方法中的方差(variance)和平均值(mean)來進(jìn)行圖像紋理特征提取，提取結(jié)果如圖6、7所示。

(2)形狀信息提取

在柵格影像中，地物的內(nèi)部信息豐富，使得在統(tǒng)計(jì)反射強(qiáng)度和紋理信息時(shí)，由于細(xì)節(jié)過多會(huì)出現(xiàn)大量的“鹽分”現(xiàn)象，進(jìn)而導(dǎo)致圖像中出現(xiàn)許多空白點(diǎn)和大量的噪聲點(diǎn)。本次實(shí)驗(yàn)采用高斯濾波方法對(duì)影像進(jìn)行濾波處理并生成0～255的灰度圖像，通過直方圖分析設(shè)置閾值進(jìn)行邊緣信息提取，生成邊緣密度二值圖像(圖8)。

圖6 方差Fig.6 Variance

圖7 平均值Fig.7 Mean

圖8 形狀信息Fig.8 Shape information

為了實(shí)驗(yàn)更充分，將提取的紋理、形狀和反射強(qiáng)度信息進(jìn)行合成，分別對(duì)只包含反射強(qiáng)度信息和多源信息融合后的圖像進(jìn)行分類，對(duì)比分類結(jié)果。

2.5 分類結(jié)果

本文利用譜聚類算法對(duì)實(shí)驗(yàn)區(qū)數(shù)據(jù)進(jìn)行分類，為了對(duì)比，分別運(yùn)用了K-means算法和高斯混合模型。試驗(yàn)中為了測(cè)試邊緣形狀、圖像紋理等信息對(duì)點(diǎn)云數(shù)據(jù)分類的影響，分別加入了形狀信息和紋理信息(圖9)。

對(duì)比不同分類結(jié)果圖可以看出，傳統(tǒng)的K-means算法和高斯混合模型都很難分出裸地和草地兩種相近的地物，并且存在漏分和錯(cuò)分現(xiàn)象，而譜聚類算法表現(xiàn)的更加優(yōu)越，對(duì)于線性地物、不同紋理的分界線都實(shí)現(xiàn)了較高的識(shí)別；同時(shí)，對(duì)于反射強(qiáng)度和紋理相近的地物類型，如林地和草地，譜聚類算法分類的效果更好，漏分、錯(cuò)分現(xiàn)象明顯減少，譜聚類對(duì)數(shù)據(jù)分布的適應(yīng)性更強(qiáng)，聚類效果更好。

a.K-means(反射強(qiáng)度信息);b.高斯混合(反射強(qiáng)度信息);c.譜聚類(反射強(qiáng)度信息);d.K-means(多源信息); e.高斯混合(多源信息);(f)譜聚類(多源信息)。圖9 分類結(jié)果對(duì)比圖Fig.9 Comparison of classification results

3 精度驗(yàn)證

本次實(shí)驗(yàn)進(jìn)行分類結(jié)果精度評(píng)價(jià)的方法是混淆矩陣[17]，在圖像上對(duì)應(yīng)不同地物隨機(jī)選取1 682個(gè)樣本進(jìn)行實(shí)驗(yàn)，對(duì)只含有反射強(qiáng)度信息和加入了紋理信息、形狀信息的分類結(jié)果分別計(jì)算混淆矩陣中的用戶精度、總體分類精度和Kappa系數(shù)等指標(biāo)(表2、圖10)。

表2 分類精度指標(biāo)對(duì)比表

圖10 分類精度對(duì)比圖Fig.10 Comparison chart of classification accuracy

從表2可以看出：譜聚類算法對(duì)不同地物的分類精度均高于K-means算法和高斯混合模型，其中譜聚類算法對(duì)只含有反射強(qiáng)度信息的分類結(jié)果的總體分類精度為80.61%，Kappa系數(shù)為0.698 2，在結(jié)合紋理和形狀信息后總體分類精度和Kappa系數(shù)分別達(dá)到81.36%和0.713 8，驗(yàn)證了譜聚類算法用于三維激光點(diǎn)云數(shù)據(jù)分類的高效性和適用性，且結(jié)合形狀和紋理信息后的多源信息圖像分類精度會(huì)高于僅含反射強(qiáng)度信息數(shù)據(jù)的分類精度；從地物類型來看，建筑用地和草地的分類精度較高，林地和裸地的分類精度一般，通過結(jié)合形狀信息和紋理信息后分類精度會(huì)有一定程度的提高，其中總體分類精度和Kappa系數(shù)分別提高了0.75%和1.56。綜上，基于遙感圖像的地物分類算法同樣適用于點(diǎn)云數(shù)據(jù)轉(zhuǎn)換成的柵格影像，并且譜聚類算法對(duì)點(diǎn)云數(shù)據(jù)柵格影像的分類精度較高，可行性和適用性都較高，該方法為典型地物分類提供了一個(gè)新途徑。

4 結(jié)論

(1)對(duì)于實(shí)驗(yàn)區(qū)內(nèi)的草地、建筑用地、林地和裸地等地物，基于反射強(qiáng)度、紋理和形狀的多源信息相較于僅基于反射強(qiáng)度信息的點(diǎn)云分類，可以更有效地進(jìn)行地類識(shí)別。

(2)譜聚類算法對(duì)點(diǎn)云數(shù)據(jù)的分類效果較好，其總體分類精度和Kappa系數(shù)都高于傳統(tǒng)的K-means算法和高斯混合模型。但譜聚類算法耗時(shí)較長(zhǎng)，后續(xù)的研究需對(duì)其參數(shù)進(jìn)行一定的優(yōu)化。

(3)將點(diǎn)云數(shù)據(jù)轉(zhuǎn)換為柵格數(shù)據(jù)后再進(jìn)行分類是一種可行的方法，對(duì)地面激光點(diǎn)云數(shù)據(jù)分類有一定的參考價(jià)值。