利用層次聚類對(duì)移動(dòng)曲面擬合濾波算法快速分類的研究

唐 菓，邢承濱，朱 磊，鄧興升，丁美青

(長沙理工大學(xué) 交通運(yùn)輸工程學(xué)院，湖南 長沙 410004)

激光雷達(dá)是利用紅外光技術(shù)對(duì)地物進(jìn)行探測。相較于數(shù)字?jǐn)z影測量的被動(dòng)式掃描，激光雷達(dá)發(fā)射脈沖是一種新型的主動(dòng)式遙感方式。激光雷達(dá)具有全天候觀測，受外界干擾少等優(yōu)點(diǎn)。并且激光雷達(dá)利用GPS和IMU可以直接獲取物體的三維坐標(biāo)，這不同于攝影測量繁瑣的數(shù)據(jù)處理過程，但是激光雷達(dá)缺少紋理，不能像數(shù)字?jǐn)z影測量那樣直觀辨別測區(qū)地物。激光雷達(dá)現(xiàn)在應(yīng)用于各個(gè)行業(yè)，如電力、林業(yè)、古建筑保護(hù)、精密儀器組裝等行業(yè)[1-2]。其高密度、高精度數(shù)據(jù)在各個(gè)行業(yè)發(fā)揮著越來越重要的作用。

如今激光雷達(dá)技術(shù)硬件已經(jīng)趨于完善。但是數(shù)據(jù)處理仍然是眾多學(xué)者研究的難點(diǎn)。其中濾波技術(shù)(地面點(diǎn)和非地面點(diǎn)的分離稱為濾波)更是重點(diǎn)。多年來，眾多學(xué)者提出了多種濾波算法：形態(tài)學(xué)濾波算法、三角網(wǎng)濾波算法、坡度算法、局部分割算法等。雖然各種算法對(duì)于部分地形取得較好濾波效果，但是仍然存在較多問題需要解決[3-4]。

張小紅首先提出移動(dòng)曲面濾波算法，該算法實(shí)現(xiàn)簡單，能夠適應(yīng)多種地形點(diǎn)云過濾[5]；蘇偉提出一種多級(jí)改進(jìn)型移動(dòng)曲面濾波算法[6]；朱笑笑提出多級(jí)移動(dòng)曲面擬合的自適應(yīng)閾值點(diǎn)云濾波方法，用于嘗試點(diǎn)云閾值自適應(yīng)的問題[7]；邢承濱提出置信區(qū)間估計(jì)理論，旨在檢驗(yàn)種子點(diǎn)選取中存在的粗差問題[8]。本文基于移動(dòng)曲面算法高差閾值難以確定的問題，提出一種層次聚類算法，利用相鄰數(shù)據(jù)點(diǎn)的水平距離與高差進(jìn)行聚類，判斷地面點(diǎn)與地物點(diǎn)的分類。

1 建立初始地面

改進(jìn)型濾波算法首先將測區(qū)劃分為大型規(guī)則格網(wǎng)，每一個(gè)格網(wǎng)中包含一定數(shù)量的數(shù)據(jù)點(diǎn)。通常格網(wǎng)大小由測區(qū)內(nèi)最大建筑物的大小決定，邢承濱提出等值線方法確定格網(wǎng)內(nèi)最大建筑物大小[9]。對(duì)于每一個(gè)網(wǎng)格建立索引，保證每一個(gè)點(diǎn)都有相應(yīng)格網(wǎng)索引與其對(duì)應(yīng)。i,j為格網(wǎng)索引號(hào)，位于同一格網(wǎng)內(nèi)數(shù)據(jù)點(diǎn)索引號(hào)相同。

i=floor[max(y)-min(y)/l]+1，

(1)

j=floor[max(x)-min(x)/l]+1.

(2)

如圖1所示，圖1(a)表示測區(qū)內(nèi)數(shù)據(jù)點(diǎn)的分布,利用xmin,ymin,xmax,ymax確定數(shù)據(jù)分布；圖1(b)為將測區(qū)分割為大型格網(wǎng),其中i為大格網(wǎng)的行號(hào)，j為大格網(wǎng)的列號(hào)，為了格網(wǎng)劃分均勻并且易于建立格網(wǎng)索引，使得格網(wǎng)為正方形，步長dx=dy。

圖1 數(shù)據(jù)點(diǎn)建立格網(wǎng)索引

對(duì)于劃分的格網(wǎng)，每一個(gè)格網(wǎng)都是一個(gè)獨(dú)立的測區(qū)，判斷每個(gè)格網(wǎng)內(nèi)的數(shù)據(jù)點(diǎn)的類型需要借助格網(wǎng)的初始曲面。擬合初始曲面通常采用選取地面種子點(diǎn)擬合的方式，由于每個(gè)格網(wǎng)大小在一定范圍內(nèi)，通常格網(wǎng)內(nèi)區(qū)域地形較為簡單，因此可以采用二次曲面擬合初始格網(wǎng)內(nèi)的地面。式(3)為二次曲面擬合公式。

f(x,y)=a1x2+a2y2+a3xy+a4x+a5y+a6.

(3)

對(duì)于格網(wǎng)內(nèi)初始地面的建立，需要借助格網(wǎng)中的初始種子點(diǎn)。為了更好選擇格網(wǎng)中的初始種子點(diǎn)，通常建立一個(gè)4×4的小型格網(wǎng)，如圖2所示將每個(gè)格網(wǎng)劃分為16個(gè)均分的小型格網(wǎng)，選取小型格網(wǎng)中最低點(diǎn)作為種子點(diǎn)，并將其擬合曲面作為真實(shí)參考地面。

如圖3所示對(duì)于某一區(qū)域內(nèi)，利用種子點(diǎn)擬合的初始地面曲面與原始地形曲面相對(duì)比，顯示地面點(diǎn)和非地面點(diǎn)的分布。

對(duì)于測區(qū)內(nèi)所有數(shù)據(jù)點(diǎn)，都會(huì)存在兩個(gè)數(shù)據(jù)值，相同的平面坐標(biāo)(x,y),一個(gè)真實(shí)的高程值HT，一個(gè)擬合地面高程Hn。真實(shí)高程值HT表示數(shù)據(jù)點(diǎn)的真實(shí)高程，Hn表示數(shù)據(jù)點(diǎn)擬合出的地面高程，如果數(shù)據(jù)點(diǎn)為地面點(diǎn)，該點(diǎn)的真實(shí)高程值與擬合高程值近似，h趨近于0；如果該數(shù)據(jù)點(diǎn)為地物點(diǎn)(非地面點(diǎn))，真實(shí)高程值與擬合高程的高差h表示該地物點(diǎn)的高程。例如，此數(shù)據(jù)點(diǎn)為樹木，h表示樹木的高度。

h=HT-Hn.

(4)

圖2 格網(wǎng)內(nèi)種子點(diǎn)的選擇

圖3 擬合地形與原始地形曲面比較

2 層次聚類自適應(yīng)閾值分類方法

2.1 三維數(shù)據(jù)的降維處理

利用相鄰點(diǎn)的水平距離和高程距離進(jìn)行分類。數(shù)據(jù)點(diǎn)三維空間坐標(biāo)轉(zhuǎn)換為二維平面坐標(biāo)Δx，dh，將水平距離和高差作為新的變量對(duì)每一個(gè)格網(wǎng)內(nèi)數(shù)據(jù)點(diǎn)進(jìn)行濾波處理。

如圖4所示，曲線表示種子點(diǎn)擬合曲面的剖面。白色圓形點(diǎn)表示真實(shí)地面，紅色圓形點(diǎn)表示地物點(diǎn)(非地面點(diǎn))。A，B為相鄰兩個(gè)數(shù)據(jù)點(diǎn)，A為地面點(diǎn)，B為建筑物點(diǎn)，hA為地面點(diǎn)高程與擬合高程之間高差，hB為B點(diǎn)地物點(diǎn)高程與擬合地面點(diǎn)之間的高差；ΔAB為AB兩點(diǎn)之間水平距離，dhAB為相鄰點(diǎn)AB的高差。

dhAB=hB-hA.

(5)

圖4 相鄰點(diǎn)高差分布圖

2.2 層次聚類算法

將上述降維處理后的Δx與dh利用層次聚類的方式進(jìn)行分類處理，以此作為濾波結(jié)果的重要判斷標(biāo)準(zhǔn)。聚類是一個(gè)將數(shù)據(jù)集中在某些方面相似的數(shù)據(jù)成員進(jìn)行分類組織的過程，聚類技術(shù)經(jīng)常被稱為無監(jiān)督學(xué)習(xí)，廣泛應(yīng)用于遙感影像的分類處理中。本文提出分層聚類閾值確定算法，通過計(jì)算樣本間或者簇間的距離進(jìn)行樣本合并，最終合并為幾個(gè)大類。本算法研究的重點(diǎn)是建筑地區(qū)，通常會(huì)將數(shù)據(jù)集合分為3類：高大地物區(qū)域、低矮地物區(qū)域、地面。

2.3 分層聚類的并項(xiàng)原理

對(duì)于層次聚類的數(shù)據(jù)都會(huì)增加新的屬性，如圖5所示，圓形、米狀、星狀表示3種不同的地物分類。高大地物區(qū)域和低矮地物區(qū)域統(tǒng)稱為地物點(diǎn)區(qū)域(非地面點(diǎn))，濾波的完成需要合并高大地物和低矮地物。層次聚類的結(jié)果只會(huì)描述點(diǎn)的分類屬性，例如一類、二類、三類，但是并不能確定每一類的性質(zhì)。算法的難點(diǎn)在于3類地形的屬性未知，層次聚類本身是一種非監(jiān)督分類，缺乏先驗(yàn)條件。本文統(tǒng)計(jì)三類地形點(diǎn)云的個(gè)數(shù)，將個(gè)數(shù)最多的數(shù)據(jù)點(diǎn)定為地面點(diǎn)，另外兩種地形合并為地物點(diǎn)。

2.4 算法過程

具體算法如圖6所示。

圖5 分層聚類并項(xiàng)

圖6 層次聚類移動(dòng)曲面算法流程

1)將測區(qū)數(shù)據(jù)劃分為規(guī)則大格網(wǎng)，并建立索引；

2)將大格網(wǎng)劃分16個(gè)小格網(wǎng)，選取小格網(wǎng)中最低點(diǎn)作為種子點(diǎn)擬合為二次曲面；

3)計(jì)算大格網(wǎng)內(nèi)擬合曲面與真實(shí)地面之間的高差作為地物點(diǎn)的高差；

4)利用k-d樹搜索格網(wǎng)內(nèi)每個(gè)種子點(diǎn)平面內(nèi)最鄰近的數(shù)據(jù)點(diǎn)，計(jì)算水平相鄰點(diǎn)水平距離；

5)利用最鄰近點(diǎn)的索引，索引出相鄰點(diǎn)的高差，計(jì)算相鄰點(diǎn)之間的相鄰高差；

6)將相鄰點(diǎn)水平距離和對(duì)應(yīng)高差距離建立層次聚類數(shù)據(jù)集，利用層次聚類算法將數(shù)據(jù)集分為三類；

7)將地物地面邊界點(diǎn)與低矮地物點(diǎn)合并為地物點(diǎn)，集中于0點(diǎn)附近數(shù)據(jù)點(diǎn)歸為地面點(diǎn)；

8)遍歷所有大格網(wǎng)，統(tǒng)計(jì)格網(wǎng)內(nèi)數(shù)據(jù)點(diǎn)的分類；

9)合并所有地面點(diǎn)，合并所有地物點(diǎn)完成測區(qū)濾波。

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

本文通過3種方式驗(yàn)證本文算法的科學(xué)性：將算法應(yīng)用于不同地區(qū)檢驗(yàn)算法的可行性，每一個(gè)測區(qū)包含整體濾波效果和格網(wǎng)內(nèi)濾波效果；將本文層次聚類算法對(duì)城市地區(qū)和村鎮(zhèn)地區(qū)進(jìn)行對(duì)比分析，驗(yàn)證本文算法在何種地形更優(yōu)越；與其他多種濾波算法進(jìn)行對(duì)比，分析本文算法的優(yōu)缺點(diǎn)和適用范圍。

3.1 實(shí)驗(yàn)數(shù)據(jù)

本文數(shù)據(jù)采用2003年攝影測量與遙感協(xié)會(huì)(ISPRS)公布的雷達(dá)點(diǎn)云數(shù)據(jù)。該數(shù)據(jù)為Optech ALTM機(jī)載激光掃描儀獲取的Vaihingen/Enz和Stuttgart市的點(diǎn)云數(shù)據(jù)。其中4個(gè)城市數(shù)據(jù)：Site1～Site4，點(diǎn)間距1.0～1.5 m；村莊數(shù)據(jù)：Site5～Site8，點(diǎn)間距2.0～3.5 m。在8個(gè)數(shù)據(jù)集中選取15組樣本，表1對(duì)每個(gè)樣本的屬性進(jìn)行詳細(xì)解釋[10]。

表1 測區(qū)地形特征

誤差的判別標(biāo)準(zhǔn)依據(jù)2003年ISPRS小組的評(píng)估報(bào)告。表2將濾波誤差分為3種：I類誤差(Type I Error)：將地面點(diǎn)誤分為非地面點(diǎn)；II類誤差(Type II Error):將非地面點(diǎn)誤分為地面點(diǎn)；總誤差(Total Error)：所有誤分點(diǎn)數(shù)/數(shù)據(jù)點(diǎn)總數(shù)[11]。

表2 濾波誤差定義

3.2 層次聚類算法濾波精度對(duì)比分析

3.2.1 層次聚類算法實(shí)驗(yàn)結(jié)果分析

本文算法統(tǒng)計(jì)15個(gè)樣本數(shù)據(jù)進(jìn)行計(jì)算，見表3。不僅計(jì)算整個(gè)樣本的三類誤差，同時(shí)隨機(jī)抽取部分格網(wǎng)樣本中數(shù)據(jù)進(jìn)行單獨(dú)格網(wǎng)的計(jì)算。

測區(qū)內(nèi)數(shù)據(jù)進(jìn)行統(tǒng)計(jì)計(jì)算，得到整體樣本，I類誤差獲得較好的擬合效果。對(duì)于地面的分離具有較高精度，均保持在5%以內(nèi)，但是II類誤差，城市測區(qū)Samp11～Samp42均能取得較好的分類。以樣本Samp41為例，樣本城區(qū)聚類低值點(diǎn)。樣本包含14個(gè)大型格網(wǎng)，統(tǒng)計(jì)計(jì)算每一個(gè)格網(wǎng)的精度。

由表4可知，對(duì)于城市地區(qū)總體效果較好，格網(wǎng)15存在NAN現(xiàn)象為本測區(qū)15格網(wǎng)內(nèi)不存在地面點(diǎn)，而分層聚類算法判斷出2個(gè)地面點(diǎn)，按照標(biāo)準(zhǔn)無法計(jì)算出I類誤差。

圖7 Samp41數(shù)據(jù)分析描述

表3 15組樣本三類誤差統(tǒng)計(jì)

表4 Samp41全部格網(wǎng)數(shù)據(jù)統(tǒng)計(jì)

圖7為樣本Samp41的數(shù)據(jù)描述。由圖7(a)和圖7(b)可以看出，地形為不連續(xù)地形。部分地形間斷明顯。由圖7(c)顯示，紅藍(lán)顏色為I，II類誤差分布，誤差主要分布在間斷地形。地形高程陡變造成分層聚類算法的部分失真，存在較為明顯的錯(cuò)誤分類。由圖7(d)地面點(diǎn)分布與圖7(b)對(duì)比，間斷地形內(nèi)利用聚類算法得到很好的消除。

從數(shù)據(jù)結(jié)論判斷層次聚類移動(dòng)曲面擬合算法可較精準(zhǔn)的對(duì)地面點(diǎn)保留和對(duì)地物點(diǎn)剔除。

3.2.2 層次聚類算法與其它濾波算法對(duì)比分析

自適應(yīng)移動(dòng)曲面算法為朱笑笑提出的新型移動(dòng)曲面濾波算法[7]，PTD(三角網(wǎng)致密化算法)濾波算法為Axelssion提出的三角網(wǎng)漸進(jìn)加密算法[3]。兩種算法都在同一組樣本下取得較好的濾波效果，其中PTD濾波算法應(yīng)用在專業(yè)軟件。表5通過3種算法對(duì)比，說明本文算法對(duì)于城區(qū)樣本濾波效果較好，I類誤差的計(jì)算結(jié)果要略優(yōu)于自適應(yīng)多級(jí)移動(dòng)曲面。且對(duì)樣本Samp11有較大程度提高，由21.52%提高到2.36%，說明利用水平距離和相鄰點(diǎn)高差聚類對(duì)于地面點(diǎn)的判斷能取得較好效果。

但是本文算法也存在一些問題。對(duì)于村鎮(zhèn)測區(qū)，Samp51—Samp71中，I類誤差精度保持較高，誤差控制在6%以內(nèi)，但是II類誤差較大，普遍在15%左右。對(duì)上述數(shù)據(jù)進(jìn)行統(tǒng)計(jì)計(jì)算，對(duì)于不同地形結(jié)構(gòu)進(jìn)行統(tǒng)計(jì)分析。

以樣本Samp53間斷地形樣本為例。將最初格網(wǎng)閾值選擇為40 m，將Samp53劃分為132個(gè)格網(wǎng)。每一個(gè)格網(wǎng)選擇最低點(diǎn)擬合為初始地面。為了體現(xiàn)本文算法對(duì)于樣本的普遍適用性，隨機(jī)選取10個(gè)格網(wǎng)。統(tǒng)計(jì)格網(wǎng)中數(shù)據(jù)點(diǎn)的I，II類誤差，如表6所示。

表5 3類算法15組樣本精度對(duì)比

表6 Samp53隨機(jī)格網(wǎng)I，II類誤差統(tǒng)計(jì)

以第3個(gè)格網(wǎng)為例，格網(wǎng)中包含297個(gè)數(shù)據(jù)點(diǎn)，利用初始種子點(diǎn)擬合為初始地面模型，圖8(a)為初始地面種子點(diǎn)的分布，圖8(b)為初始地面模型。利用初始地面模型，計(jì)算每一個(gè)數(shù)據(jù)點(diǎn)的擬合高程。利用k-d樹搜索數(shù)據(jù)點(diǎn)水平距離最鄰近點(diǎn)，計(jì)算相鄰點(diǎn)高差。將水平距離和相鄰高差作為層次聚類的依據(jù)，如圖8(c)表示格網(wǎng)聚類成果。通常地形分類為3類：地面點(diǎn)、低矮地物、大型地物。地面點(diǎn)的特點(diǎn)：水平距離和高差距離集中在0附近。大型地物出現(xiàn)在高差較大地區(qū)，小型地物會(huì)出現(xiàn)在高差較小地區(qū)。具體分布需利用分層聚類算法計(jì)算，圖8(d)描述測區(qū)誤差點(diǎn)的分布。

同時(shí)圖9(a)選擇Samp53樣本中第一個(gè)格網(wǎng)，格網(wǎng)中有數(shù)據(jù)點(diǎn)248個(gè)。又聚類算法呈現(xiàn)出相鄰點(diǎn)高差存在巨大差異但是不呈現(xiàn)集聚性，由此判斷為山體斷裂線。由初始數(shù)據(jù)點(diǎn)呈現(xiàn)測區(qū)為山體地形。同時(shí)由二類誤差分布顯示。絕大多數(shù)誤差集中在斷裂線處。斷裂線處誤差主要因?yàn)閿嗔丫€處地面點(diǎn)誤分為地物點(diǎn)，造成I類誤差偏高。斷裂線處個(gè)別山體陡峭，真實(shí)高程遠(yuǎn)遠(yuǎn)大于擬合高程，造成層次聚類中存在錯(cuò)誤分類。

圖8 Samp53第3格網(wǎng)數(shù)據(jù)描述

圖9 Samp53第1格網(wǎng)數(shù)據(jù)描述

圖9(b)為第1格網(wǎng)中數(shù)據(jù)點(diǎn)聚類為3類，星號(hào)為I類誤差，點(diǎn)號(hào)為II類誤差，米狀為地面點(diǎn)，米點(diǎn)具有集聚性，并且集中在0軸附近。圖9(c)所示I，II類都有明顯的特征，高差較大，但是I類誤差水平距離較小，II類誤差水平距離較大，高差距離較小。但是從圖9(d)和圖9(e)，過濾后的地面區(qū)域較為平滑，剔除絕大多數(shù)地物點(diǎn)。整體分析，層次聚類對(duì)于山體地區(qū)還有較大需要改進(jìn)的地方，對(duì)于城鎮(zhèn)地區(qū)濾波效果良好。

4 結(jié)束語

本文提出的聚類算法利用擬合地面和真實(shí)地形做高差，利用水平距離和相鄰點(diǎn)高差層次聚類，在城區(qū)取得較好效果，但是在山體地區(qū)，由于山體斷裂線和陡變的地形，使得擬合地面與真實(shí)地形高差較大，造成II類誤差較大，總體而言濾波精度較好。由于實(shí)驗(yàn)數(shù)據(jù)有限，本文算法并沒有進(jìn)行對(duì)于大規(guī)模地形的檢校，對(duì)于大規(guī)模地形的適應(yīng)性處理還有待研究。