陳斐然，李 浩，李語旻，陳仁喜

（1.河海大學(xué) 地球科學(xué)與工程學(xué)院，江蘇 南京 210098）

基于改進坡度的自適應(yīng)數(shù)學(xué)形態(tài)學(xué)點云濾波

陳斐然1，李 浩1，李語旻1，陳仁喜1

（1.河海大學(xué) 地球科學(xué)與工程學(xué)院，江蘇 南京 210098）

首先介紹了傳統(tǒng)形態(tài)學(xué)算法，并分析了其應(yīng)用方面的不足；然后提出了一種基于改進坡度參數(shù)的自適應(yīng)形態(tài)學(xué)濾波方法；最后利用實際點云數(shù)據(jù)進行了實驗。該方法通過劃分格網(wǎng)自動計算當前區(qū)域坡度的平均值，從而得到適當?shù)钠露葏?shù)用于后續(xù)高度閾值的計算。實驗結(jié)果表明，該方法能更有效地識別地面點和地物點，對被建筑物包圍的植被點濾除效果明顯。

LiDAR點云；數(shù)學(xué)形態(tài)學(xué)；濾波；分類

LiDAR不僅能快速獲取大面積點云數(shù)據(jù)，而且能從中得到高精確DEM，因此受到各測繪行業(yè)的高度重視。點云濾波是LiDAR數(shù)據(jù)處理的重要內(nèi)容之一，也是最先需要解決的問題。濾波分類精度的高低直接影響點云數(shù)據(jù)的使用質(zhì)量。點云濾波即對在地形表面隨機分布的離散點進行識別與處理，區(qū)分落在地面上的點（地面點）和落在地面物體上的點（植被﹑建筑物﹑車輛等地物點），分別形成地面點與地物點數(shù)據(jù)集合。

目前國內(nèi)外已有許多專家學(xué)者對點云濾波方法進行了詳細研究[1-2]。張子林[3]等使用基于TIN的濾波算法對高程突變地物進行了濾波，但對低矮的地面物體進行濾波時，誤差會顯著增大。參考文獻[4]～[6]提出一種基于坡度變化的濾波算法，但該算法需提前知道地形坡度值，并確定窗口大小，區(qū)域內(nèi)每個點都需進行坡度計算，并同其他點進行比較，以確定該點是否為地面點。同時，如果僅根據(jù)坡度設(shè)置統(tǒng)一的閾值，則會濾除掉一些真實的地形細節(jié)特征，造成分類誤差，所以必須根據(jù)作業(yè)區(qū)域的實際地形狀況設(shè)置不同的過濾參數(shù)值。另外一種常用的濾波算法是基于數(shù)學(xué)形態(tài)學(xué)原理的濾波算法[7-11]。數(shù)學(xué)形態(tài)學(xué)算法不僅可用于點云數(shù)據(jù)濾波，也可用于建筑物或道路點云數(shù)據(jù)的提取。然而現(xiàn)有的基于經(jīng)典數(shù)學(xué)形態(tài)學(xué)算法的LiDAR點云濾波方法存在以下不足：①無法根據(jù)地形特點動態(tài)調(diào)整窗口。利用一個固定窗口大小的結(jié)構(gòu)元素很難過濾各種尺寸的地物點，若窗口過小則只有較小的地物點被去除，如汽車或植被，但此時窗口無法完全包住大型建筑物，使得大型建筑物不能被有效濾除；而較大的濾波窗口對于含有山坡和低矮山包的地區(qū)而言，很容易將山包和山坡當成非地面點而錯誤濾除。②無法顧及窗口內(nèi)地形的細節(jié)特征。使用一個固定的“高度閾值”作為判斷是否為地面點的唯一標準顯然是不合理的。③坡度參數(shù)沒有考慮實際地形情況。坡度參數(shù)的取值由用戶自己輸入，并沒有真正考慮地形起伏的實際情況，對后期濾波的準確程度產(chǎn)生影響。

本文對傳統(tǒng)數(shù)學(xué)形態(tài)學(xué)濾波算法進行了改進，提出了一種基于改進坡度參數(shù)的自適應(yīng)數(shù)學(xué)形態(tài)學(xué)濾波算法，不僅能自動調(diào)整窗口大小，而且可以通過劃分格網(wǎng)來自動計算當前區(qū)域坡度的平均值，得到適當?shù)钠露葏?shù)用于后續(xù)高度閾值的計算，使得對LiDAR點云的識別更加準確。

1 數(shù)學(xué)形態(tài)學(xué)濾波方法

1.1 經(jīng)典數(shù)學(xué)形態(tài)學(xué)濾波

數(shù)學(xué)形態(tài)學(xué)是基于集合運算原理進行特征提取的，其中集合表示圖像中不同的對象。腐蝕與膨脹分別用于減少或增大圖像中特征形狀的尺寸。腐蝕運算即以窗口內(nèi)高程的最小值代替窗口中心的高程值，而膨脹運算即以窗口內(nèi)高程的最大值代替窗口中心的高程值。

在采用數(shù)學(xué)形態(tài)學(xué)方法對LiDAR點云進行濾波時，一般使用這兩個基礎(chǔ)算子的組合，主要是開運算和閉運算，其表達式為：

式中，Θ為腐蝕運算符號；⊕為膨脹運算符號；○為開運算符號；●為閉運算符號。

開運算是先腐蝕后膨脹的運算過程，可先將比結(jié)構(gòu)尺寸小的地物點去除，再膨脹恢復(fù)被腐蝕的建筑物邊緣形狀。這種去除比結(jié)構(gòu)元素窗口尺寸小的物體同時保留比結(jié)構(gòu)元素窗口大的物體特性，常被用于提取LiDAR點云數(shù)據(jù)的地面點，因此在LiDAR點云濾波中得到廣泛應(yīng)用。

1.2 漸進式數(shù)學(xué)形態(tài)學(xué)濾波

針對單一窗口和單一閾值數(shù)學(xué)形態(tài)學(xué)濾波方法存在的問題，ZHANG K[7]等提出的漸進式數(shù)學(xué)形態(tài)學(xué)濾波方法可從LiDAR點云中濾除大部分的地物點（見圖1）。

圖1 漸進式數(shù)學(xué)形態(tài)學(xué)濾波方法

由圖1可知，在第一次濾波處理后，比窗口尺寸L1小的樹木點會被濾除，而大型建筑物則被保留；同樣尺寸小于L1的地面特征被削平，山坡最高點的高程與第一次濾波后的高程差值為dhmax(t),1，而山坡上隨機點P的高程與第一濾波后的高程差值為dhp,1。在隨后的迭代過程中，窗口大小增加到L2，此時大型建筑物尺寸小于當前窗口尺寸，大型建筑物點被濾除，而地形特征點前后兩次濾波的高程差值為dhmax(t),2。

引入坡度參數(shù)和高度閾值可避免濾波后起伏地形的高點被濾除，但無法保留地形的細節(jié)特征，將地形坡度參數(shù)s和每次迭代高程差閾值dhT,k分別定義為：

式中，dhmax(t),k為第k次濾波后與前一次濾波的最大高程差；wk為第k次濾波的窗口大??；c為格網(wǎng)間距；dh0為初始高差閾值；dhmax為最大高差閾值。如果某個點前后兩次濾波的高程差值小于閾值，則分類成地面點，否則為地物點。

對于城市區(qū)域的LiDAR數(shù)據(jù)而言，若dhT,k取值適當，則可確保比s小的地形特征表面被保留下來。城市中主要的非地面點對象為建筑物﹑汽車﹑植被等，通過逐漸增大的濾波窗口重復(fù)使用開運算，汽車與植被點在首次運算時就會被濾除，而高大的房屋會在最后幾次運算中被過濾掉。一般情況下，將最低房屋的高度設(shè)為最大高差閾值dhmax，以保證復(fù)雜建筑物被完全識別濾除。

2 自適應(yīng)數(shù)學(xué)形態(tài)學(xué)濾波方法

[7]中，s的取值是由用戶根據(jù)地形實際情況自己輸入的，城市與山區(qū)地形的s取值有較大差別，且城市內(nèi)部中的s也是不同的。本文采用首先規(guī)則劃分格網(wǎng)并進行高程賦值，然后利用坡度公式求解當前區(qū)域各格網(wǎng)的坡度值，最終取其平均值作為區(qū)域近似坡度值的改進坡度參數(shù)形態(tài)學(xué)濾波方法，以提高算法在不同地形環(huán)境下的自適應(yīng)能力和完善性。改進算法流程如圖2所示。

圖2 改進算法流程圖

格網(wǎng)尺寸過大會降低LiDAR點云高程的精度，而網(wǎng)格尺寸過小則會造成數(shù)據(jù)計算量過大﹑運行緩慢。將區(qū)域總面積除以區(qū)域內(nèi)點云總數(shù)得到的平均點間距作為規(guī)則化點云格網(wǎng)的大小c，這樣做的目的是使得大部分格網(wǎng)中只包含1個LiDAR點，滿足一個激光點對應(yīng)一個格網(wǎng)的要求。在此基礎(chǔ)上，由已知區(qū)域內(nèi)最大和最小的x﹑y值與某一個格網(wǎng)中激光點坐標計算得到該格網(wǎng)的行列號，并依次遍歷全部點云數(shù)據(jù)。

格網(wǎng)規(guī)則劃分完成后，則需對每個格網(wǎng)進行高程賦值。這里分為兩種情況，一種是當有多個點落在同一個格網(wǎng)中時，則取其中最低點的高程作為該格網(wǎng)的高程；另一種是當格網(wǎng)中沒有任何激光點時，則將最鄰近點的高程值賦給該格網(wǎng)。所有格網(wǎng)賦值完成后，再利用新生成的規(guī)則格網(wǎng)高程數(shù)據(jù)建立初始LiDAR點云的數(shù)字表面模型（DSM）用于下一步坡度計算。

對于初始DSM，為了自動計算地表面的坡度值，本文使用二階坡度計算公式來計算每個格網(wǎng)所對應(yīng)的坡度值：

式中，Zi,j+1﹑Zi,j-1﹑Zi+1,j﹑Zi-1,j分別對應(yīng)相應(yīng)格網(wǎng)行列號的格網(wǎng)高程值；c為格網(wǎng)大小。

將式（5）﹑（6）代入式（7）中就可計算得到各格網(wǎng)的坡度S(i, j)，但該坡度并不是所需要的實際地形坡度，而是地表坡度。在得到的地表坡度中，除了一小部分地面坡度外，其余絕大多數(shù)是植被﹑建筑物等地物點產(chǎn)生的坡度。為了能獲得實際地形的近似坡度值，必須進行預(yù)處理。本文首先利用合適窗口大小的開算子遍歷全部地表坡度值，以剔除建筑物和部分植被對坡度的影響。由于坡度平均值代表區(qū)域內(nèi)地形的普遍情形，為了進一步消除建筑物周圍植被對坡度的干擾，再對當前區(qū)域的全部格網(wǎng)坡度值進行相應(yīng)運算，取其平均值作為該區(qū)域的近似坡度，并用于后續(xù)高度閾值的計算，最大限度地保留實驗區(qū)域的地形細節(jié)特征。

3 實驗結(jié)果與分析

3.1 實驗數(shù)據(jù)與結(jié)果

實驗數(shù)據(jù)采用ISPRS網(wǎng)站提供的城市點云數(shù)據(jù)，主要是由大小不一的建筑物及不同高度的植被所組成。實驗區(qū)域大小為1 450 m×505 m，點云高差為50 m，點云間距為1.5 m。原始點云如圖3所示，可以看出實驗區(qū)分布著大﹑中﹑小不規(guī)則形狀的建筑物，同時有道路﹑汽車﹑樹木等地物。圖4為采用基于改進坡度參數(shù)的自適應(yīng)數(shù)字形態(tài)學(xué)濾波方法濾波后的結(jié)果。

在原始點云中選取一處復(fù)雜城市環(huán)境作為樣本，共有26 984個離散激光點，圖5b為采用漸進式數(shù)學(xué)形態(tài)學(xué)濾波算法的結(jié)果，圖5c為采用本文算法的結(jié)果。

圖3 原始數(shù)據(jù)

圖4 改進算法濾波結(jié)果

圖5 復(fù)雜區(qū)域濾波結(jié)果

3.2 實驗結(jié)果分析

實驗結(jié)果表明，基于改進坡度參數(shù)的數(shù)學(xué)形態(tài)學(xué)濾波方法能更有效地區(qū)分地面點與地物點。在圖4中，絕大多數(shù)植被與建筑物都被濾除掉，同時大部分地面細節(jié)特征得以保留下來。從圖5b可以看出，漸進式數(shù)學(xué)形態(tài)學(xué)濾波對被建筑物包圍的植被點濾除效果不佳，而本文提出的改進算法能有效去除建筑物周圍的植被點，保留地面真實信息。

為了能更好地定量分析實驗結(jié)果，將樣本區(qū)域的人工判讀分類結(jié)果作為準確分類結(jié)果（地面點14 023 個，非地面點12 961個），統(tǒng)計參考文獻[7]的算法和本文提出的改進算法的誤差。第Ⅰ類誤差﹑第Ⅱ類誤差和總誤差的定義如表1所示。表2為在樣本區(qū)域濾波結(jié)果中兩種算法的誤差對比。

表1 濾波誤差的定義

表2 濾波結(jié)果定量分析

從表2可以看出，改進算法的第Ⅰ類誤差較參考文獻[7]算法有明顯降低，說明地面點識別分類更加精確；而第Ⅱ類誤差沒有太大的變化，總誤差則明顯低于參考文獻[7]算法。

4 結(jié) 語

本文所提出的基于改進坡度參數(shù)的自適應(yīng)數(shù)學(xué)形態(tài)學(xué)濾波方法對現(xiàn)有的數(shù)學(xué)形態(tài)學(xué)濾波算法進行了改進，采用了自動劃分格網(wǎng)求解坡度平均值作為坡度參數(shù)，并用于后續(xù)高度閾值計算的方法，有效解決了被建筑物包圍的植被點濾除問題。實驗結(jié)果顯示，該方法能有效區(qū)分地面點與非地面點，并進一步保留了地形的細節(jié)特征，降低了濾波誤差。

參考文獻

[1] 王金亮,陳聯(lián)君.激光雷達點云數(shù)據(jù)的濾波算法述評[J].遙感技術(shù)與應(yīng)用,2010,25(5):632-638

[2] ZHANG K, Whitman D. Comparison of Three Algorithms for Filtering Airbone LiDAR Data[J].Photogrammetric Engineering & Remote Sensing,2005,71(3):313-324

[3] 張子林,門林杰.基于TIN改進的LiDAR濾波算法 [J].工程地球物理學(xué)報,2012,9(3):362-364

[4] 許可,王運巧,胡少興.基于坡度與迭代運算的LiDAR點云數(shù)據(jù)濾波[J].機械工程與自動化,2016(2):71-73

[5] 張寧寧,楊英寶,于雙.基于坡度和區(qū)域生長的城市LiDAR點云濾波方法[J].地理空間信息,2016,14(3):30-32

[6] Sithole G. Filtering of Laser Altimetry Data Using a Slope Adaptive Filter[J]. International Archives of Photogrammetry and Remote Sensing,2001,34(3/W4):22-24

[7] ZHANG K,CHEN S C, Whitman D, et al. A Progressive Morphological Filter for Removing Nonground Measurement from LiDAR Data[J]. IEEE Transactions on Geoscience and Remote Sensing,2003,41(4):872-882

[8] 隋立春,張熠斌,柳艷,等.基于改進的數(shù)學(xué)形態(tài)學(xué)算法的LiDAR點云數(shù)據(jù)濾波[J].測繪學(xué)報,2010,39(4):390-396

[9] CHEN Q, GONG P , Baldocchi D ,et al. Filtering Airborne Laser Scanning Data with Morphological Methods[J].Photogrammetry Engineering and Remote Sensing,2007,73(2):175-185

[10] 羅伊萍,姜挺,龔志輝,等.基于自適應(yīng)和多尺度數(shù)學(xué)形態(tài)學(xué)的點云數(shù)據(jù)濾波方法[J].測繪科學(xué)技術(shù)學(xué)報,2009,26(6):426-429

[11] 李峰,崔希民,袁德寶,等.改進坡度的LiDAR點云形態(tài)學(xué)濾波算法[J].大地測量與地球動力學(xué),2012,32(5):128-132

P228

B

1672-4623（2017）09-0028-04

10.3969/j.issn.1672-4623.2017.09.010

2016-09-18。

項目來源：國家自然科學(xué)基金資助項目（41471276）。

陳斐然，碩士研究生，研究方向為LiDAR數(shù)據(jù)處理與地圖制圖。