一種改進(jìn)的形態(tài)學(xué)機(jī)載LiDAR濾波算法

朱 磊

（浙江中海達(dá)空間信息技術(shù)有限公司，浙江 湖州 313299）

機(jī)載激光雷達(dá)（LiDAR）以飛機(jī)為平臺，在空中實(shí)現(xiàn)主動對地觀測，可以準(zhǔn)確獲取測區(qū)內(nèi)目標(biāo)地物的幾何位置信息。它將激光測距技術(shù)、差分GNSS技術(shù)以及慣性測量技術(shù)進(jìn)行高度融合，具有自動化度高、精度高、全天候和受外部環(huán)境影像小的特點(diǎn)。因此，機(jī)載激光雷達(dá)測量技術(shù)被廣泛應(yīng)用于數(shù)字地面模型（DSM）提取、電力巡檢、災(zāi)害監(jiān)測預(yù)防以及林業(yè)調(diào)查管理等領(lǐng)域。

通過無人機(jī)搭載機(jī)載LiDAR可以獲取大量的地表點(diǎn)云數(shù)據(jù)，還需要對這些地表點(diǎn)云數(shù)據(jù)進(jìn)行濾波才能準(zhǔn)確提取地物的坐標(biāo)和高程信息，從而高精度地表達(dá)測區(qū)地形信息[1-2]。研究者往往會根據(jù)區(qū)域的不同特點(diǎn)選擇不同原理的濾波算法進(jìn)行分析，但是在植被覆蓋多的地形復(fù)雜區(qū)域采用已有的濾波算法，其地形的適應(yīng)性及濾波的準(zhǔn)確度等方面存在明顯的不足?；谏鲜鰡栴}，該文將薄板樣條插值（TPS）理論引入點(diǎn)云濾波算法中，提出了一種基于TPS插值理論改進(jìn)的數(shù)字形態(tài)學(xué)點(diǎn)云濾波算法，對LiDAR數(shù)據(jù)進(jìn)行濾波測試，驗(yàn)證了改進(jìn)算法的適應(yīng)性和可行性，提高了機(jī)載LiDAR在復(fù)雜區(qū)域應(yīng)用時所獲取的點(diǎn)云數(shù)據(jù)的精度和效率。

1 濾波算法

1.1 經(jīng)典的形態(tài)學(xué)濾波算法

經(jīng)典的形態(tài)學(xué)濾波算法流程是先使用柵格化的激光點(diǎn)云數(shù)據(jù)，取每個柵格內(nèi)最低點(diǎn)進(jìn)行開運(yùn)算，計算高差閾值，從而計算開運(yùn)算前后柵格化點(diǎn)云數(shù)據(jù)的高程差，再對高程差大于閾值的點(diǎn)云數(shù)據(jù)進(jìn)行判別和過濾，進(jìn)而獲得濾波結(jié)果[3]。經(jīng)典的形態(tài)學(xué)濾波算法中提出的高差閾值是通過結(jié)合實(shí)際測區(qū)地形的坡度值及其濾波窗口2個參數(shù)確定的，利用線性方程逐漸改變窗口的大小，其線性方程如公式（1）所示。

式中：wk為濾波窗口大??；k為系數(shù)，k=0，1，2，3，…，m；b為窗口的初始大小，最大濾波窗口為2mb+1。

選用2kb+1作為窗口尺寸，可以滿足濾波窗口以中心點(diǎn)對稱的需求，這樣可以對濾波的開運(yùn)算過程進(jìn)行優(yōu)化。對測區(qū)內(nèi)高程的閾值來說，該文根據(jù)測區(qū)內(nèi)實(shí)際地形的坡度計算相關(guān)數(shù)據(jù)，該方法的優(yōu)點(diǎn)是原理比較簡單，算法比較容易實(shí)現(xiàn)。但是當(dāng)采用該方法對較大窗口進(jìn)行濾波處理時，不可避免地會出現(xiàn)將地形起伏點(diǎn)判別為非地面點(diǎn)的問題，在復(fù)雜的城市、丘陵等地區(qū)濾波效果較差。

1.2 薄板樣條插值

薄板樣條插值是基于樣條線的數(shù)據(jù)插值和平滑技術(shù)，通過測區(qū)內(nèi)給定的不規(guī)則分布的多個控制點(diǎn)來生成一個平滑的曲面，在生成的曲面上要保證所有控制點(diǎn)的坡度值是最小的[4-7]。

與其他插值相比，TPS插值能很好地對地形進(jìn)行模擬，更真實(shí)地反映地形坡度變化的特征，TPS插值中多個控制點(diǎn)可以不規(guī)則分布，比較適用于離散點(diǎn)插值的情況，當(dāng)應(yīng)用于復(fù)雜條件下時，同樣可以得到很好的插值結(jié)果。

1.3 改進(jìn)的基于TPS插值的形態(tài)學(xué)算法

以經(jīng)典數(shù)學(xué)形態(tài)學(xué)濾波算法為基礎(chǔ)，結(jié)合TPS插值，根據(jù)濾波窗口大小的不同進(jìn)行TPS插值，然后對地形起伏度進(jìn)行計算。地形起伏度是基于TPS插值后的DEM經(jīng)膨脹運(yùn)算的結(jié)果和初始擬合的DEM進(jìn)行差值運(yùn)算的，這種算法的優(yōu)勢是可以在一定程度上保留地形起伏的細(xì)節(jié)，提高單獨(dú)使用形態(tài)學(xué)進(jìn)行濾波時數(shù)據(jù)的精度?；赥PS差值的形態(tài)學(xué)算法的具體實(shí)施步驟如下。

1.3.1 點(diǎn)云去噪

激光點(diǎn)云中的噪聲點(diǎn)主要是因個別地物高程值出現(xiàn)異常而產(chǎn)生的，其與周圍地物點(diǎn)的點(diǎn)云數(shù)據(jù)存在一定的差值，可以根據(jù)噪聲點(diǎn)與周圍點(diǎn)的高程差的不同，將噪聲點(diǎn)分為2個類別，即高位噪聲點(diǎn)和低位噪聲點(diǎn)，再分別對高位和低位噪聲點(diǎn)進(jìn)行去除。首先，采用基于密度的去噪算法，通過輸入密度參數(shù)對某點(diǎn)一定距離內(nèi)的LiDAR點(diǎn)云數(shù)據(jù)進(jìn)行計算，以刪除小于閾值的噪聲點(diǎn)。其次，在基于密度去噪算法的基礎(chǔ)上進(jìn)行人工去噪，通過人工操作去除剩余噪聲點(diǎn)。

1.3.2 建立格網(wǎng)

在上述去噪的基礎(chǔ)上，對點(diǎn)云數(shù)據(jù)進(jìn)行格網(wǎng)化處理和格網(wǎng)劃分，建立格網(wǎng)索引并確定格網(wǎng)的邊界，格網(wǎng)索引的行號和列號分別用i和j來表示。將分化后每個格網(wǎng)內(nèi)的點(diǎn)云數(shù)據(jù)的最小值賦值給該格網(wǎng)，但是在點(diǎn)云數(shù)據(jù)格網(wǎng)化后存在一部分格網(wǎng)中不包括點(diǎn)云數(shù)據(jù)的情況，針對這種情況下的格網(wǎng)，可以采用最鄰近法賦值法，即將格網(wǎng)周圍距離最近的點(diǎn)賦值于該格網(wǎng)，這樣避免因點(diǎn)云數(shù)據(jù)柵格化的缺失對結(jié)果精度造成影響。

1.3.3 形態(tài)學(xué)開運(yùn)算

繼續(xù)對點(diǎn)云數(shù)據(jù)進(jìn)行形態(tài)學(xué)開運(yùn)算，當(dāng)進(jìn)行開運(yùn)算時，要考慮當(dāng)前窗口尺寸的大小，通過運(yùn)算得出該層數(shù)據(jù)前、后的高差變化。

1.3.4 TPS插值計算地形起伏度

從格網(wǎng)化的LiDAR點(diǎn)云數(shù)據(jù)中挑選格網(wǎng)高程最小值并賦值給該格網(wǎng)，對當(dāng)前窗口下所有的最小高程點(diǎn)進(jìn)行TPS插值，從而得到平滑的擬合曲面，再對曲面進(jìn)行膨脹運(yùn)算，進(jìn)而得到地形的坡度值。

1.3.5 濾波判斷

對初始點(diǎn)云進(jìn)行TPS插值處理后，得到對應(yīng)柵格的地形起伏度。地形坡度S為常量，為了提高適用性，在得到地形起伏度后，要計算不同層級的地形坡度，從而求出高差閾值，再將其與該點(diǎn)窗口下形態(tài)學(xué)開運(yùn)算前后的高差作差，將差值與閾值進(jìn)行比較，如果差值大于閾值，就被視為地物點(diǎn)剔除；如果差值小于閾值，就視為地面點(diǎn)，并代入下一級重新建立格網(wǎng)索引，重復(fù)該過程直到窗口尺寸小于最小窗口尺寸。不同層級的地形坡度如公式（2）所示。

式中：S（i，j）為坡度；FX為水平梯度；FY為垂直梯度；nc為常量；z為高程值；C為柵格大小。

2 試驗(yàn)與分析

2.1 試驗(yàn)數(shù)據(jù)

為了驗(yàn)證該文提出的改進(jìn)濾波算法的可行性，選取飛馬D2000搭載LiDAR2000測量系統(tǒng)對試驗(yàn)區(qū)進(jìn)行航攝并測量激光雷達(dá)，試驗(yàn)區(qū)域?yàn)榍鹆甑貛В匚镱愋桶ǚ课?、植被、道路以及溝壑等，面積約為5.56 km2，最大高程為268.480 m，最小高程值為200.589 m，點(diǎn)云密度為1.13 點(diǎn)/m2。根據(jù)所獲取的點(diǎn)云（如圖1所示），在試驗(yàn)區(qū)內(nèi)選取8個樣本區(qū)域，包括橋梁、規(guī)則建筑物、植被、陡坡、不連續(xù)陡坡、植被、橋梁以及不規(guī)則建筑物用于該文的試驗(yàn)研究。

圖1 研究區(qū)原始點(diǎn)云示意圖

2.2 精度評估與評價

隨機(jī)抽取8個樣本區(qū)域點(diǎn)云數(shù)據(jù)中的800個點(diǎn)作為參考數(shù)據(jù)，參考數(shù)據(jù)作為對濾波精度評價的依據(jù)，該試驗(yàn)參考數(shù)據(jù)采用“TerraSolid軟件+人工判讀”的方式區(qū)分地面點(diǎn)和非地面點(diǎn)，其精度滿足試驗(yàn)參考數(shù)據(jù)的要求。

該試驗(yàn)點(diǎn)云分類精度采用國際通用的由混淆矩陣推導(dǎo)的Ⅰ類誤差、Ⅱ類誤差和總誤差進(jìn)行評價，并將其作為檢驗(yàn)濾波效果的精度評價指標(biāo)。

2.3 數(shù)據(jù)處理與分析

分別使用經(jīng)典形態(tài)學(xué)濾波算法和該文提出的改進(jìn)算法對8組樣本數(shù)據(jù)進(jìn)行濾波處理，得到的濾波結(jié)果如圖2所示，2種濾波算法總誤差對比如圖3所示，評價標(biāo)準(zhǔn)見表1，Ⅰ類誤差、Ⅱ類誤差和總誤差精度見表2。

表1 評價標(biāo)準(zhǔn)表

表2 精度評價指標(biāo)

圖2 濾波結(jié)果

圖3 2種濾波算法的總誤差對比圖

由圖3可知，該文提出的基于TPS插值的改進(jìn)的數(shù)學(xué)形態(tài)學(xué)濾波算法中的Ⅰ類誤差、Ⅱ類誤差比經(jīng)典形態(tài)學(xué)濾波算法的誤差小。

由表2和圖3可知，對經(jīng)典算法中Ⅰ類誤差進(jìn)行統(tǒng)計發(fā)現(xiàn)，只有樣本三和樣本六有2組樣本數(shù)據(jù)的誤差小于10%，誤差分別為6.34%、5.43%，其余6組數(shù)據(jù)的誤差為10%～20%；對Ⅱ類誤差進(jìn)行統(tǒng)計發(fā)現(xiàn)，只有數(shù)據(jù)樣本六的誤差小于20%，其余樣本數(shù)據(jù)的誤差都大于20%，樣本二最大，達(dá)到52.38%；對總誤差進(jìn)行統(tǒng)計發(fā)現(xiàn)，只有樣本三和樣本六有2組樣本數(shù)據(jù)的誤差小于10%，誤差值分別為9.69%、7.76%，其余樣本的誤差為7.76%～26.79%；在基于TPS插值的濾波改進(jìn)算法所得的試驗(yàn)結(jié)果中，8組數(shù)據(jù)中Ⅰ類誤差的統(tǒng)計結(jié)果是所有樣本的誤差都小于7%，樣本一誤差為1.35%；Ⅱ類誤差的統(tǒng)計結(jié)果有4組樣本數(shù)據(jù)的誤差小于10%，分別為樣本一、樣本二、樣本七和樣本八，其中樣本一的誤差最低，為4.39%；總誤差的統(tǒng)計結(jié)果顯示8組數(shù)據(jù)的誤差為1.96%～9.69%，最小值為樣本一，誤差為1.96%。對這8組數(shù)據(jù)進(jìn)行分析可以說明，基于TPS插值的改進(jìn)算法的各類誤差均比經(jīng)典算法低。

3 結(jié)語

該文基于2種濾波算法計算試驗(yàn)區(qū)內(nèi)8個點(diǎn)云數(shù)據(jù)特征樣本，通過處理LiDAR點(diǎn)云數(shù)據(jù)得出以下2個結(jié)論：1）該文提出的基于TPS插值形態(tài)學(xué)濾波算法濾波效果較好，尤其是在不規(guī)則建筑物和陡坡的情況中。2）該文提出的改進(jìn)的形態(tài)學(xué)機(jī)載激光LiDAR算法能大幅度地降低Ⅰ類誤差、Ⅱ類誤差和總誤差，經(jīng)過改進(jìn)的基于TPS的形態(tài)學(xué)濾波算法生成的DEM數(shù)據(jù)更貼近真實(shí)地形。

綜上所述，基于TPS插值的漸進(jìn)形態(tài)學(xué)濾波算法可以比較完整地體現(xiàn)、保留地形特征，提高了LiDAR點(diǎn)云濾波的精度和效率，在植被茂密和建筑物較多、較復(fù)雜的區(qū)域的應(yīng)用效果也比經(jīng)典算法好。