基于自適應二級網(wǎng)格的二次曲面點云濾波方法

皇永波 盧小平＊ 周雨石 劉雪晴 朱夢豪

（河南理工大學測繪與國土信息工程學院，河南焦作 454003）

三維激光掃描技術(shù)在大范圍數(shù)字高程模型的高精度實時獲取、城市三維模型重建、局部區(qū)域的地理信息獲取等方面表現(xiàn)出強大優(yōu)勢[1]。原始的點云數(shù)據(jù)包含大量地物信息，應用于建設(shè)數(shù)字高程模型時，需要對點云數(shù)據(jù)進行濾波得到地面點云。目前，點云數(shù)據(jù)的地面點濾波方法主要包括坡度濾波算法[2]、數(shù)學形態(tài)學濾波算法[3]、最小二乘曲面擬合濾波算法等[4]。述濾波算法均存在各自的問題，例如坡度濾波閾值的自適應性低，且容易錯誤的選擇地面種子；數(shù)學形態(tài)學濾波算法存在窗口尺寸不易確定和細節(jié)地形方塊效應問題；最小二乘曲面擬合濾波算法則對擬合半徑依賴性大。本文提出了一種自適應大小的二級網(wǎng)格與曲面結(jié)合的濾波方法，根據(jù)兩種地形點云數(shù)據(jù)進行試驗，對比常用的坡度濾波和形態(tài)學濾波，分析試驗結(jié)果的1類誤差、2類誤差和總誤差。

1 適應大小的二級網(wǎng)格構(gòu)建方法

1.1 虛擬網(wǎng)格坡度濾波

虛擬網(wǎng)格是對點云劃區(qū)分塊，將點云區(qū)域分成M×N個大小相等的矩形網(wǎng)格，使每個網(wǎng)格都包含離散的激光腳點，網(wǎng)格大小M和N的計算公式：

式中：xmax和xmin——點云數(shù)據(jù)中x坐標的最大值和最小值；ymax和ymin——y坐標的最大值和最小值。

確定M和N的值后，對每個激光腳點建立索引（X，Y），建立索引的具體計算方法：

式中：INT——取整數(shù)；x和y——激光腳點的坐標；xmin和ymin——點云數(shù)據(jù)二維平面坐標的最小值。

構(gòu)建虛擬網(wǎng)格后遍歷網(wǎng)格內(nèi)所有的激光腳點，將每個網(wǎng)格內(nèi)高程最低的激光腳點作為網(wǎng)格內(nèi)的地面種子點。使用虛擬網(wǎng)格和地面種子點實現(xiàn)坡度濾波算法，設(shè)置坡度閾值T，計算虛擬網(wǎng)格內(nèi)初始地面種子點P0和激光腳點Pi之間的坡度值，將坡度值小于閾值T的點分類為地面點，否則分為非地面點。Pi和P0兩點之間的坡度值Gi：

式中：li0——Pi和P0之間的歐式距離；hi0——Pi和P0兩點之間的高差。

1.2 建立自適應大小的二級網(wǎng)格

統(tǒng)計坡度濾波后的地面點，根據(jù)地面點重新構(gòu)建虛擬網(wǎng)格。將每個虛擬網(wǎng)格中高程最低點選為初始地面種子點p0，再以p0點作為節(jié)點，對虛擬網(wǎng)格進行分割，將虛擬網(wǎng)格初步分割成4個不相等的二級網(wǎng)格。保證每一網(wǎng)格內(nèi)激光腳點的個數(shù)保持在10個左右，因此合并小于10個激光腳點的二級網(wǎng)格，再對二級網(wǎng)格內(nèi)的腳點建立索引[5]。構(gòu)建二級網(wǎng)格：

（1）選擇坡度濾波后的地面點重新構(gòu)建虛擬網(wǎng)格，網(wǎng)格邊長S的大小由以下公式計算：

式中：m——粗提取的地面點云區(qū)域的面積；n——地面點云激光腳點的數(shù)量。

（2）在每個網(wǎng)格內(nèi)選擇初始地面種子點，以種子點為中心將虛擬網(wǎng)格初步分割成四個二級網(wǎng)格。

（3）根據(jù)虛擬網(wǎng)格內(nèi)激光腳點與初始地面種子點之間的位置關(guān)系建立二級網(wǎng)格的索引，數(shù)學關(guān)系為：

式中：（x0，y0）——地面種子點坐標；（xi，yi）——網(wǎng)格內(nèi)激光腳點坐標；（a，b，c，d）——虛擬網(wǎng)格內(nèi)激光腳點所對應的二級網(wǎng)格的索引編號。

（4）統(tǒng)計每個二級網(wǎng)格內(nèi)的點數(shù)，4個二級網(wǎng)格中有3個以上網(wǎng)格的點數(shù)小于10個，重新合并為一級網(wǎng)格，否則按順時針方向?qū)⑿∮?0個激光腳點的二級網(wǎng)格與相鄰的二級網(wǎng)格進行合并，并建立合并后的網(wǎng)格索引，計算公式為：

式中：A，B，C，D——合并后二級網(wǎng)格的編號。

1.3 基于二級網(wǎng)格擬合二次曲面

在每個二級網(wǎng)格中選取6個高程最低的地面點，采用二次曲面擬合的方法擬合曲面[4]：

式中：（x，y）——地面點坐標；z——擬合曲面后的高程；(a0,a1,a2,a3,a4,a5)——多項式系數(shù)。

多項式的系數(shù)可由求解多項式的最小殘差平方和得到：

式中：Emin——最小殘差平方和；zi——選取二級網(wǎng)格內(nèi)的激光腳點的高程。

E的最小值是由ai控制，ai對求偏導，將求解公式轉(zhuǎn)換為求極值的過程：

式中：A——待定系數(shù)為[a0,a1,a2,a3,a4,a5]T的矩陣；B和C——包含（xi，yi，zi）的矩陣。

將每個虛擬網(wǎng)格內(nèi)粗提取的地面點，帶入到公式(10)中計算[a0,a1,a2,a3,a4,a5]T的值，即求解出二次曲面的多項式的系數(shù)。

1.4 本文濾波方法

本文濾波方法流程如圖1所示。

圖1 算法流程

二次曲面由粗提取的地面點擬合而成，必須剔除粗差點減少誤差。

計算激光腳點擬合前后均方根誤差剔除粗差點：

式中：νi——激光腳點擬合前后差值；——差值平均值；n——二級網(wǎng)格內(nèi)的點數(shù)；σ——差值的均方根誤差。

將σ作為誤差閾值，若地面激光腳點擬合前后的差值大于2倍的均方根誤差，將此腳點判為粗差點進行剔除。剔除粗差點后重新計算二次曲面多項式，并且重新計算曲面點的均方根誤差2σ。將坡度濾波得到的非地面點坐標代入對應的二次曲面多項式計算高程值，重新執(zhí)行式（11），將νi小于2σ的點歸類為地面點，否則歸類為非地面點。

2 試驗分析

本文針對不同的地形區(qū)域進行試驗，試驗1使用平緩的城市區(qū)域點云，試驗2使用地形起伏區(qū)域的點云。濾波結(jié)果評價標準使用ISPRS于2003年提出的濾波誤差的評判標準，如表1所示。

表1 點云濾波誤差評判標準

地勢較為平坦的地區(qū)，地形坡度的值為0°～10°，地勢起伏較大的地區(qū)的坡度閾值為10°～30°，試驗1的坡度閾值取為10°，坡度濾波網(wǎng)格大小為3 m，試驗2的坡度閾值取為15°，網(wǎng)格大小為5 m。

2.1 地勢平緩區(qū)域試驗

試驗數(shù)據(jù)是DublinCity數(shù)據(jù)集中的一部分，數(shù)據(jù)編號為T_315500_234500_NE。原始數(shù)據(jù)集如圖2所示。

圖2 點云數(shù)據(jù)集

點云總數(shù)12 198 716個，地面點5 957 008個，非地面點個數(shù)為6 241 708個。

本文濾波方法與虛擬網(wǎng)格坡度濾波和數(shù)學形態(tài)學濾波相比1類誤差減少19.589%、21.293%，2類誤差減少15.573%、17.257%，總誤差減少17.534%、19.202%。

試驗1的誤差統(tǒng)計如表2所示。

表2 試驗1的誤差統(tǒng)計 單位：%

試驗結(jié)果如圖3所示。

圖3 試驗1結(jié)果

2.2 地勢起伏區(qū)域試驗

原始點云數(shù)據(jù)如圖4所示。

圖4 原始點云數(shù)據(jù)圖

總點云數(shù)3 123 944個。地面點云和非地面點云通過人工進行分類，地面激光腳點1 225 880個，非地面激光腳點1 898 064個。本文濾波方法與虛擬網(wǎng)格坡度濾波和數(shù)學形態(tài)學濾波相比1類誤差減少22.861%、28.916%，2類誤差減少26.942%、23.047%，總誤差減少25.351%、25.080%。地面點云結(jié)果如圖5所示。

圖5 地面點云結(jié)果

試驗2誤差統(tǒng)計如表3所示。

表3 試驗2誤差統(tǒng)計 單位：%

3 結(jié)語

本文介紹了自適應大小的二級網(wǎng)格與二次曲面相結(jié)合的點云的濾波方法，采用擬合二次曲面能夠減少地面種子點選取造成的誤差。由于三種濾波方法均不能刪除建筑物的所有邊緣輪廓點，后續(xù)可針對此問題進行深入研究。