利用層次聚類對移動曲面擬合濾波算法快速分類的研究

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

(長沙理工大學 交通運輸工程學院，湖南 長沙 410004)

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

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

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

1 建立初始地面

改進型濾波算法首先將測區劃分為大型規則格網，每一個格網中包含一定數量的數據點。通常格網大小由測區內最大建筑物的大小決定，邢承濱提出等值線方法確定格網內最大建筑物大小[9]。對于每一個網格建立索引，保證每一個點都有相應格網索引與其對應。i,j為格網索引號，位于同一格網內數據點索引號相同。

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

(1)

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

(2)

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

圖1 數據點建立格網索引

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

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

(3)

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

如圖3所示對于某一區域內，利用種子點擬合的初始地面曲面與原始地形曲面相對比，顯示地面點和非地面點的分布。

對于測區內所有數據點，都會存在兩個數據值，相同的平面坐標(x,y),一個真實的高程值HT，一個擬合地面高程Hn。真實高程值HT表示數據點的真實高程，Hn表示數據點擬合出的地面高程，如果數據點為地面點，該點的真實高程值與擬合高程值近似，h趨近于0；如果該數據點為地物點(非地面點)，真實高程值與擬合高程的高差h表示該地物點的高程。例如，此數據點為樹木，h表示樹木的高度。

h=HT-Hn.

(4)

圖2 格網內種子點的選擇

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

2 層次聚類自適應閾值分類方法

2.1 三維數據的降維處理

利用相鄰點的水平距離和高程距離進行分類。數據點三維空間坐標轉換為二維平面坐標Δx，dh，將水平距離和高差作為新的變量對每一個格網內數據點進行濾波處理。

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

dhAB=hB-hA.

(5)

圖4 相鄰點高差分布圖

2.2 層次聚類算法

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

2.3 分層聚類的并項原理

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

2.4 算法過程

具體算法如圖6所示。

圖5 分層聚類并項

圖6 層次聚類移動曲面算法流程

1)將測區數據劃分為規則大格網，并建立索引；

2)將大格網劃分16個小格網，選取小格網中最低點作為種子點擬合為二次曲面；

3)計算大格網內擬合曲面與真實地面之間的高差作為地物點的高差；

4)利用k-d樹搜索格網內每個種子點平面內最鄰近的數據點，計算水平相鄰點水平距離；

5)利用最鄰近點的索引，索引出相鄰點的高差，計算相鄰點之間的相鄰高差；

6)將相鄰點水平距離和對應高差距離建立層次聚類數據集，利用層次聚類算法將數據集分為三類；

7)將地物地面邊界點與低矮地物點合并為地物點，集中于0點附近數據點歸為地面點；

8)遍歷所有大格網，統計格網內數據點的分類；

9)合并所有地面點，合并所有地物點完成測區濾波。

3 實驗分析

本文通過3種方式驗證本文算法的科學性：將算法應用于不同地區檢驗算法的可行性，每一個測區包含整體濾波效果和格網內濾波效果；將本文層次聚類算法對城市地區和村鎮地區進行對比分析，驗證本文算法在何種地形更優越；與其他多種濾波算法進行對比，分析本文算法的優缺點和適用范圍。

3.1 實驗數據

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

表1 測區地形特征

誤差的判別標準依據2003年ISPRS小組的評估報告。表2將濾波誤差分為3種：I類誤差(Type I Error)：將地面點誤分為非地面點；II類誤差(Type II Error):將非地面點誤分為地面點；總誤差(Total Error)：所有誤分點數/數據點總數[11]。

表2 濾波誤差定義

3.2 層次聚類算法濾波精度對比分析

3.2.1 層次聚類算法實驗結果分析

本文算法統計15個樣本數據進行計算，見表3。不僅計算整個樣本的三類誤差，同時隨機抽取部分格網樣本中數據進行單獨格網的計算。

測區內數據進行統計計算，得到整體樣本，I類誤差獲得較好的擬合效果。對于地面的分離具有較高精度，均保持在5%以內，但是II類誤差，城市測區Samp11～Samp42均能取得較好的分類。以樣本Samp41為例，樣本城區聚類低值點。樣本包含14個大型格網，統計計算每一個格網的精度。

由表4可知，對于城市地區總體效果較好，格網15存在NAN現象為本測區15格網內不存在地面點，而分層聚類算法判斷出2個地面點，按照標準無法計算出I類誤差。

圖7 Samp41數據分析描述

表3 15組樣本三類誤差統計

表4 Samp41全部格網數據統計

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

從數據結論判斷層次聚類移動曲面擬合算法可較精準的對地面點保留和對地物點剔除。

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

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

但是本文算法也存在一些問題。對于村鎮測區，Samp51—Samp71中，I類誤差精度保持較高，誤差控制在6%以內，但是II類誤差較大，普遍在15%左右。對上述數據進行統計計算，對于不同地形結構進行統計分析。

以樣本Samp53間斷地形樣本為例。將最初格網閾值選擇為40 m，將Samp53劃分為132個格網。每一個格網選擇最低點擬合為初始地面。為了體現本文算法對于樣本的普遍適用性，隨機選取10個格網。統計格網中數據點的I，II類誤差，如表6所示。

表5 3類算法15組樣本精度對比

表6 Samp53隨機格網I，II類誤差統計

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

同時圖9(a)選擇Samp53樣本中第一個格網，格網中有數據點248個。又聚類算法呈現出相鄰點高差存在巨大差異但是不呈現集聚性，由此判斷為山體斷裂線。由初始數據點呈現測區為山體地形。同時由二類誤差分布顯示。絕大多數誤差集中在斷裂線處。斷裂線處誤差主要因為斷裂線處地面點誤分為地物點，造成I類誤差偏高。斷裂線處個別山體陡峭，真實高程遠遠大于擬合高程，造成層次聚類中存在錯誤分類。

圖8 Samp53第3格網數據描述

圖9 Samp53第1格網數據描述

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

4 結束語

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