一種DEM 輔助下的LiDAR 點(diǎn)云PTD 濾波改進(jìn)算法

鄭 斌，鄒學(xué)忠，李小昱

（1.江蘇省測(cè)繪工程院，江蘇 南京 210013）

LiDAR作為一種主動(dòng)式、高精度、快速、智能化實(shí)現(xiàn)地物三維形態(tài)特性的測(cè)繪高新技術(shù)，已廣泛應(yīng)用于基礎(chǔ)測(cè)繪、電力、林業(yè)和城市三維建模等領(lǐng)域，其中點(diǎn)云濾波是點(diǎn)云數(shù)據(jù)后處理的關(guān)鍵步驟[1]。點(diǎn)云濾波技術(shù)是從LiDAR 數(shù)據(jù)點(diǎn)中獲取地形表面點(diǎn)的過程，大致分成基于坡度、基于數(shù)學(xué)形態(tài)學(xué)、基于TIN、基于聚類分割、基于曲面擬合以及基于機(jī)器學(xué)習(xí)的濾波6 類[2-4]。近年來，一些研究者通過實(shí)驗(yàn)測(cè)試對(duì)上述主流濾波算法進(jìn)行了精度對(duì)比，其中基于機(jī)器學(xué)習(xí)的濾波算法在訓(xùn)練樣本足夠的情況下能獲得較好的濾波效果，但數(shù)據(jù)處理時(shí)間較長，對(duì)計(jì)算機(jī)硬件的要求較高[5-6]；基于聚類分割的濾波算法因能獲得更多聚類后的語義信息，可達(dá)到較高的濾波精度[7-8]；基于TIN和曲面擬合的濾波算法采用漸進(jìn)迭代的方式也能獲得較好的濾波精度[9-12]。復(fù)雜環(huán)境下的點(diǎn)云濾波技術(shù)仍是當(dāng)前LiDAR點(diǎn)云研究的熱點(diǎn)方向。本文在經(jīng)典漸進(jìn)加密不規(guī)則三角網(wǎng)（PTD）濾波算法的基礎(chǔ)上，提出了一種DEM輔助下的PTD濾波改進(jìn)算法。該算法將往期DEM數(shù)據(jù)提取的地形高程信息和地形梯度信息作為輔助，改進(jìn)PTD中初始地面種子點(diǎn)的選取方法，優(yōu)化地面點(diǎn)判斷參數(shù)，并對(duì)往期DEM數(shù)據(jù)和現(xiàn)勢(shì)LiDAR點(diǎn)云數(shù)據(jù)之間的地形變化進(jìn)行檢測(cè)和處理，以改善濾波效果。

1 算法原理

傳統(tǒng)PTD 濾波算法的穩(wěn)定性和濾波效果均較好，能適應(yīng)相對(duì)復(fù)雜的地形，但初始種子點(diǎn)的選取對(duì)濾波結(jié)果影響較大，在不同坡度地形條件下需輸入不同的地面點(diǎn)判斷參數(shù)才能獲得較好的結(jié)果[13]。因此，本文提出了一種DEM輔助下的PTD濾波改進(jìn)算法，旨在優(yōu)化初始種子點(diǎn)的選取，降低不同地形對(duì)修改參數(shù)的依賴，最終提升點(diǎn)云處理的自動(dòng)化程度、減少后期人工作業(yè)的工作量。

1）通過統(tǒng)一坐標(biāo)系統(tǒng)、組織點(diǎn)云數(shù)據(jù)和剔除點(diǎn)云噪聲點(diǎn)對(duì)LiDAR點(diǎn)云數(shù)據(jù)進(jìn)行預(yù)處理，再按照DEM地形面約束將預(yù)處理后的點(diǎn)云數(shù)據(jù)分為非地面點(diǎn)和待分類點(diǎn)，分別加入非地面點(diǎn)集合VNon-Ground和待分類點(diǎn)集合VNon-Classfication。具體處理流程為：①遍歷點(diǎn)云數(shù)據(jù)中的每個(gè)點(diǎn)，若不是僅一次回波點(diǎn)或多次回波中的末次回波點(diǎn)，則該點(diǎn)為非地面點(diǎn)，加入VNon-Ground，否則保留為待分類點(diǎn)，加入VNon-Classfication；②設(shè)定高差閾值HThreshold，并判斷HNon-Classfication中的每個(gè)點(diǎn)與相同平面坐標(biāo)位置DEM 的高程差值是否超過HThreshold，若超過，則判斷為待定非地面點(diǎn)，加入VNon-Ground，否則保留為待分類點(diǎn)，加入VNon-Classfication（圖1）。

圖1 非地面點(diǎn)剔除示意圖

2）利用DEM 數(shù)據(jù)的高程信息，從點(diǎn)云數(shù)據(jù)中選取可信的地面點(diǎn)，加入地面點(diǎn)集合VGround。首先對(duì)Li?DAR點(diǎn)云數(shù)據(jù)建立規(guī)則格網(wǎng)索引；然后在每個(gè)格網(wǎng)內(nèi)選擇高程最小點(diǎn)作為種子點(diǎn)加入VGround，格網(wǎng)尺寸應(yīng)不小于車輛、低矮植被等近地面地物的尺寸，實(shí)驗(yàn)中設(shè)置為5～10 m；再從DEM數(shù)據(jù)中提取斷裂線區(qū)域，將落入斷裂線區(qū)域的格網(wǎng)劃分為更小的格網(wǎng)，并分別在每個(gè)小格網(wǎng)內(nèi)選擇高程最小點(diǎn)作為種子點(diǎn)加入VGround（圖2）。

圖2 種子點(diǎn)選取示意圖

3）將地面點(diǎn)作為地面種子點(diǎn)構(gòu)建初始地形TIN模型，并利用DEM數(shù)據(jù)的坡度信息對(duì)地面點(diǎn)判斷參數(shù)進(jìn)行優(yōu)化；然后在初始地形TIN模型的基礎(chǔ)上對(duì)VNon-Classification進(jìn)行地面點(diǎn)判別，并加入VGround。具體處理流程為：①遍歷VNon-Classification中的所有點(diǎn)，按照式（1）計(jì)算該點(diǎn)處的相對(duì)地面高度閾值Td、地形角閾值Tθ、d和θ；②判斷是否滿足d＜Td，且θ＜Tθ，若滿足，則該點(diǎn)標(biāo)記為地面點(diǎn)，加入VGround。

式中，a0、a1、b0、b1分別為常數(shù)系數(shù)；為某激光腳點(diǎn) 所 處 位 置 的DEM 地 形 梯 度，=(f(x)，f(y)) ；f(x)、f(y)分別為x和y方向的坡度值；為地形梯度向量的模。

4）利用TIN 模型對(duì)非地面點(diǎn)進(jìn)行分析處理，將VNon-Ground中地形凸起變化區(qū)域的點(diǎn)重新分類到地面點(diǎn)，并加入VGround?；趨^(qū)域生長的變化區(qū)域地面點(diǎn)探測(cè)原理見圖3。具體處理流程為：①遍歷VNon-Ground中的所有點(diǎn)，計(jì)算該點(diǎn)與TIN模型中所在三角形的相對(duì)地面高度d、地形角θ和最近地面點(diǎn)距離L；②判斷是否同時(shí)滿足d＜Td、θ＜Tθ、L＜TL，若滿足，則該點(diǎn)標(biāo)記為地面點(diǎn)，加入VGround，插入到初始地形TIN模型中并更新，其中xp、yp為待判斷點(diǎn)P的x、y坐標(biāo)，xp'、yp'為距離待判點(diǎn)最近的地面點(diǎn)P′的x、y坐標(biāo)；③重復(fù)步驟②，直至不再有新的點(diǎn)滿足條件為止。

圖3 基于區(qū)域生長的變化區(qū)域地面點(diǎn)探測(cè)原理

由于在地形變化區(qū)域，地形情況已發(fā)生明顯變化，往期DEM數(shù)據(jù)不能提供可靠的參考信息。實(shí)驗(yàn)中Td小于非地面點(diǎn)粗剔除中的HThreshold，Tθ應(yīng)設(shè)置較大，TL為格網(wǎng)索引尺寸的兩倍。

2 算法實(shí)現(xiàn)

2.1 DEM輔助下的PTD濾波實(shí)驗(yàn)

本文以6 km 的DEM 標(biāo)準(zhǔn)圖幅作為實(shí)驗(yàn)樣本（圖4），樣本區(qū)內(nèi)包含平原、丘陵、水域、山地等地物類型，可充分驗(yàn)證概率波策略對(duì)于復(fù)雜地形的適應(yīng)性和對(duì)于不同地形的普適性。樣本區(qū)LiDAR點(diǎn)云獲取時(shí)間為2017 年3 月，平均點(diǎn)密度為1 m；用于輔助濾波的DEM 為航攝法獲取并制作的往期基礎(chǔ)測(cè)繪成果，數(shù)據(jù)現(xiàn)勢(shì)性為2007年10月，格網(wǎng)間距為5 m。輔助濾波完成后生成的DEM模型見圖5，可以看出，實(shí)驗(yàn)樣本的濾波效果良好，尤其在城區(qū)和山區(qū)，無需人工編輯；但由于水面會(huì)吸收或反射激光雷達(dá)信號(hào)，使得水面上的激光腳點(diǎn)較少，導(dǎo)致河流或水塘等區(qū)域的濾波效果不佳（圖6），因此仍需單獨(dú)進(jìn)行補(bǔ)點(diǎn)和添加特征線等處理，以保證地表模型的真實(shí)與自然。

圖4 實(shí)驗(yàn)樣本區(qū)域影像與往期DEM示意圖

圖5 實(shí)驗(yàn)樣本輔助分類DEM效果圖

2.2 實(shí)驗(yàn)結(jié)果評(píng)價(jià)與分析

本文采用激光點(diǎn)云數(shù)據(jù)的3 類誤差評(píng)價(jià)對(duì)本文的濾波效果進(jìn)行定量評(píng)價(jià)[14-15]。I 類誤差又稱去真誤差，即將地面反射的激光腳點(diǎn)數(shù)據(jù)誤歸為地物反射的激光腳點(diǎn)并濾波；II類誤差又稱納偽誤差，即將地物反射的激光腳點(diǎn)數(shù)據(jù)誤歸為地面反射的激光腳點(diǎn)數(shù)據(jù)并保留；對(duì)兩種誤差進(jìn)行加權(quán)求和，可有效判斷激光點(diǎn)云數(shù)據(jù)濾波分類的結(jié)果，反映濾波方法的可行性，濾波總體精度的計(jì)算公式為：

式中，a為將地面激光點(diǎn)誤認(rèn)為地物激光點(diǎn)的數(shù)目；b為正確分類的地面激光點(diǎn)數(shù)目；c為將地物反射的激光點(diǎn)誤認(rèn)為地面激光點(diǎn)的數(shù)目；d為正確分類的地物激光點(diǎn)數(shù)目。

Kappa 系數(shù)常用于一致性檢驗(yàn)，基于混淆矩陣計(jì)算得到所需參數(shù)的精度指標(biāo)，范圍在0～1 之間，其中0.61～0.80 表示高度一致性，0.81～1 表示幾乎完全一致，計(jì)算公式為：

本文將經(jīng)典PTD濾波算法結(jié)果與改進(jìn)的PTD濾波算法結(jié)果進(jìn)行對(duì)比，結(jié)果見表1，可以看出，改進(jìn)算法的I 類誤差從6.12%降至3.3%，II 類誤差從25.59%降至15.45%，總體精度從83.57%提升至90.27%；且Kappa系數(shù)為0.81，表示分類結(jié)果幾乎完全一致。

表1 實(shí)驗(yàn)樣本濾波分類精度統(tǒng)計(jì)表

在江蘇省“十三五”基礎(chǔ)測(cè)繪機(jī)載LiDAR點(diǎn)云數(shù)據(jù)獲取與處理項(xiàng)目中，利用該算法對(duì)江蘇省全域10.26萬km2范圍內(nèi)的點(diǎn)云數(shù)據(jù)進(jìn)行濾波精細(xì)化分類，達(dá)到了良好的效果。江蘇省測(cè)繪產(chǎn)品質(zhì)量監(jiān)督檢驗(yàn)站對(duì)項(xiàng)目生產(chǎn)的機(jī)載激光LiDAR 數(shù)據(jù)和DEM 進(jìn)行了抽樣檢測(cè)，結(jié)果表明各樣本的LiDAR 激光點(diǎn)云高程中誤差、DEM 高程中誤差質(zhì)量均合格，符合GB/T 18316-2008《數(shù)字測(cè)繪成果質(zhì)量檢查與驗(yàn)收》的指標(biāo)要求[16-17]。

3 結(jié) 語

本文提出的DEM輔助下的PTD濾波算法，對(duì)經(jīng)典PTD濾波進(jìn)行了改進(jìn)和優(yōu)化。經(jīng)過實(shí)驗(yàn)驗(yàn)證，改進(jìn)算法提高了實(shí)驗(yàn)樣本的總體精度和Kappa系數(shù)，改善了濾波效果。在江蘇省“十三五”省級(jí)基礎(chǔ)測(cè)繪機(jī)載LiDAR項(xiàng)目中，利用往期5 m格網(wǎng)DEM輔助完成了全省1 m點(diǎn)密度的點(diǎn)云濾波分類，經(jīng)檢驗(yàn)成果質(zhì)量良好。江蘇省大范圍應(yīng)用效果表明，該濾波算法可在已有較低密度格網(wǎng)DEM 輔助下，對(duì)高密度待濾波Li?DAR點(diǎn)云進(jìn)行自動(dòng)化處理，尤其是在植被覆蓋密集的山區(qū)、建筑物間距較小的城區(qū)能發(fā)揮更大優(yōu)勢(shì)，降低點(diǎn)云濾波的錯(cuò)誤率，減少后期人工干預(yù)的工作量。當(dāng)然，對(duì)于地形已發(fā)生明顯變化的區(qū)域，已有DEM數(shù)據(jù)不能提供可靠的參考信息，不適用于該方法，需將變化區(qū)域單獨(dú)劃分出來處理。