基于數學形態(tài)學的LiDAR數據濾波新方法

羅伊萍，姜 挺，王 鑫，張 銳，羅 勝

(信息工程大學測繪學院，河南鄭州 450052)

基于數學形態(tài)學的LiDAR數據濾波新方法

羅伊萍，姜 挺，王 鑫，張 銳，羅 勝

(信息工程大學測繪學院，河南鄭州 450052)

數學形態(tài)學濾波是從激光雷達點云數據中識別地面點、創(chuàng)建數字高程模型的一種重要方法。在分析現(xiàn)有濾波方法的優(yōu)劣性以及數學形態(tài)學濾波方法存在的問題基礎上，提出一種新的具有一定自適應性的數學形態(tài)學濾波算法。該方法通過分析LiDAR數據的特點，利用形態(tài)學算子提取原始點云數據的空白區(qū)域，進行精確的重采樣生成DSM;然后將DSM多尺度形態(tài)學濾波的結果作為初始的DEM，分析由于局部地形突變導致可能存在誤分類區(qū)域;最后提取該區(qū)域邊緣外的地面點作為新的種子點，利用最小二乘平面擬合的方法進行區(qū)域生長搜索誤分類地面點，從而提高分類精度。試驗結果表明該方法能夠有效識別地面點和地物點，并且保留地形的細節(jié)信息。

激光雷達;濾波;數學形態(tài)學濾波;數字高程模型

一、引 言

機載激光雷達(light deteation and ranging，LiDAR)技術自20世紀80年代引入攝影測量領域以來，受到廣大研究者的極大關注。LiDAR系統(tǒng)利用機載激光雷達測距系統(tǒng)和GPS/IMU能直接獲取地面點的三維坐標，形成離散的、不規(guī)則的三維點云數據。為了從離散的點云數據中獲取數字高程模型(DEM)，必須對點的屬性進行區(qū)分，將原始的點云數據分為地面點集和非地面點集，這個過程被稱為濾波。

目前，已有不少國內外學者對激光雷達點云數據的濾波方法進行了深入的研究，主要的濾波算法可分為以下三類:基于表面的內插濾波算法［1-2］、基于區(qū)域的濾波算法［3-5］和基于約束曲面的濾波算法［6-7］。基于表面的內插濾波算法的核心思想是通過一個較粗的起始DEM，逐步從備選數據點篩選并內插加密DEM達到分類的目的。因此，濾波結果受到初始DEM影響較大并且誤差會隨著迭代的過程積累。而基于Snake樣條曲面等約束曲面的濾波方法認為地面是一個連續(xù)且平緩變化的表面，這樣就過于強調地形的平緩變化而忽略了地形的復雜性［8］。

數學形態(tài)學濾波方法屬于基于區(qū)域的濾波算法，與基于表面的內插濾波算法相反，它是一種自下而上、從局部出發(fā)擴展到全局的濾波方法。借鑒柵格圖像的處理思想，根據激光點云數據生成的深度圖像，利用形態(tài)學開運算剔除高于地面的點，通過逐步增大濾波窗口，得到逼近地形的一個表面。然而在濾波過程中，僅考慮窗口內的地形特征，容易受到局部地形的影響，例如道路旁存在溝渠，當濾波窗口尺寸大于道路寬度時會導致路面點被分為非地面點。另外，濾波過程中的參數需要用戶根據不同地形設置，自適應性不強［3-4］。針對數學形態(tài)學濾波算法存在的問題，本文提出一種新的多尺度數學形態(tài)學濾波策略用于LiDAR數據濾波。

二、多尺度數學形態(tài)學濾波方法

數學形態(tài)學的語言和理論是基于集合運算原理提取圖像中的特征。腐蝕和膨脹運算是形態(tài)學圖像處理的基礎，通常用于“減少”(腐蝕)或“增大”(膨脹)圖像中特征形狀的尺寸。基于灰度圖像的腐蝕和膨脹運算是在結構元素定義的鄰域內選擇圖像像素值和結構元素相作用后的最小或最大像素值。

在LiDAR點云數據生成的規(guī)則化數字表面模型(DSM)中，腐蝕和膨脹運算可以定義為:

腐蝕

膨脹

其中，f為規(guī)則化的DSM;g為結構元素;Z(i，j)為腐蝕或膨脹運算后規(guī)則化的DSM中第i行第j列的高程值;w為結構元素的窗口。

腐蝕和膨脹運算組合后，形成開運算和閉運算用于LiDAR點云數據濾波，開運算和閉運算可以定義為:

開運算

閉運算

式中，開運算通過先腐蝕后膨脹的方法，先將比結構元素尺寸小的樹木點等非地面點從原始點云數據中移除掉，然后通過膨脹運算恢復被腐蝕掉的建筑物邊緣形狀。然而，利用一個固定窗口大小的結構元素很難移除各種尺寸的非地面物體，如果窗口過小，則大部分地面點被保留，只有小的非地面目標，比如車或者樹木被移除，而城區(qū)的大型建筑物不能被移除。因此必須利用多尺度的濾波窗口進行迭代計算，在迭代過程中逐步增大濾波窗口的尺寸。

每次迭代開運算時，濾波的窗口尺寸wk按指數增長，wk=2bk+1，式中 k為迭代的次數，k=1，2，…，M;b為初始的窗口尺寸大小，一般為2。為了保證建筑物都被去除掉，最后一次迭代的窗口尺寸wk必須大于試驗區(qū)域中最大建筑物的面積，一般迭代次數為5或6次。為了確保在濾波的過程中保留細節(jié)的地形特征，設置每次迭代高差閾值dhT，k

式中，wk為第k次濾波的窗口大小;c為DSM的格網間距;dh0和dhmax分別為最小和最大高差閾值;s為地形坡度參數。從式(5)可以看出，當濾波窗口增大時，高差閾值dhT，k隨之增大，增幅大小由地形坡度s決定。如果DSM中某個格網點開運算前后的高程值之差小于本次迭代的高差閾值dhT，k，則認為是地面點，否則認為是非地面點。

三、DEM提取過程

1.點云數據格網化

1)在點云數據的xy平面范圍內，選擇合適的格網間距構成m×n平面格網。對于每個離散的激光腳點，根據其平面坐標將其分配到相應的格網中。為了使盡可能多的激光點用于初始的DSM，一般選取比點間距稍小的格網間距值。由于點云數據分布的不均勻性，如圖1(a)所示，一個格網中可能會沒有或者有一個或多個激光點。如果一個格網內有多個激光點落入，高程值越低的點作為地面點的概率越大，因此選取最低高程值作為格網的取值。記錄下每個格網中激光點的原始索引號，用于最后的精度評定。

2)在初始的DSM中，由于選取了格網內所有激光腳點的最低值，因此由激光掃描系統(tǒng)誤差和隨機誤差引起的少量無意義低值粗差點也被保留下來?？紤]到粗差點與相鄰點之間的高差形成脈沖波峰，在初始表面模型中的5×5鄰域內，如果中心點高程值與24個鄰域點高程值的差值均大于高差閾值3 m，那么判定該點為噪聲點并且該格網中沒有激光點。

以m×n的二值標記矩陣記錄下對應格網中是否有激光點，矩陣值為0則沒有激光點，如圖1(c)所示。

2.DSM重采樣

DSM重采樣是針對格網化后的DSM中的空白格網進行數據內插。由于激光雷達的視場角比傳統(tǒng)攝影測量相機視場角小，飛行時容易造成掃描漏洞，如圖1(c)的頂部有細長的掃描漏洞。另外，水體對1 064～1 600 nm波段的激光有吸收效應，也會造成數據空白區(qū)域，如圖1(c)的左側中部區(qū)域。如果存在高層建筑物，由于建造物的遮擋也會產生數據空白區(qū)域。由于這種數據的不均衡性，導致數據內插的搜索半徑的大小很難確定，尤其是存在大量湖泊河流的區(qū)域，并且對建筑物遮擋的空白區(qū)域內插生成的DSM擴大了建筑物的面積。因此，要先針對這種大面積的數據空白區(qū)域進行數據重采樣，然后以較小的搜索半徑對剩余的空白格網內插，以提高重采樣的效率和精度。

首先，利用3×3的形態(tài)學窗口對二值標記矩陣作閉運算，去掉由于激光點分布不均勻造成的稀疏空白格網。其次，為了消除水域上的零散激光點，采用小面積消除法，求得標記矩陣內所有值為1的四連通區(qū)域，其個數小于40的連通區(qū)域認為是由水體反射的激光點造成的，將它們賦值為0。然后，二值矩陣取反，再次使用小面積消除法消除零散的無激光點格網，得到圖1(d)。此時僅保留了面積大于40個格網數的水域和因飛行掃描漏洞造成的數據空白區(qū)域，為了求得這些區(qū)域的邊緣，將標記矩陣用3×3的形態(tài)學窗口作膨脹運算，再減去原標記矩陣即得到如圖1(e)所示的邊緣。一般認為水域是比較平緩的高程一致性區(qū)域，因此將水域邊緣的最低值作為水域的高程取值。掃描漏洞和建筑物遮擋區(qū)域也作相同處理。最后，由合適的搜索半徑對其余零散的空白格網進行最鄰近值內插，得到重采樣的DSM，如圖1(b)所示。

圖1 點云數據重采樣

3.多尺度數學形態(tài)學濾波

在點云數據規(guī)則化和重采樣后，本文采用上述的基于坡度和高差閾值的多尺度數學形態(tài)學濾波方法。針對城市區(qū)域地形起伏較小以及建筑物的最大面積一般不會超過120 m×120 m的實際情況，濾波的參數選擇如表1所示。

表1 濾波參數設置

4.質量控制

由于數學形態(tài)學濾波僅考慮窗口內的地形特征，在局部地形有突變的情況下，會導致相當數目的地面點分類為非地面點，從而影響分類的精度。在這種情況下由于地面的連續(xù)性，誤分類區(qū)域在一定范圍內與另外的地面區(qū)域相連接，通過提取相連接地面區(qū)域的邊緣點作為種子點，通過區(qū)域生長法搜索誤分類區(qū)域中的地面點，從而改進分類精度。

將多尺度數學形態(tài)學濾波后的結果作為初始的DEM，并且從中提取非地面點集的二值化標記圖像。用小面積消除法剔除標記圖像中小面積的樹木點，再用3×3的窗口對標記圖像進行膨脹運算再減去標記圖像，得到大面積非地面點區(qū)域的外邊緣點，即相連接地面區(qū)域的邊緣點。從這些邊緣點中，根據它們3×3鄰域的最大和最小高程值之差，剔除建筑物的外邊緣點，保留有可能存在誤分類地面點集的邊緣點作為種子點。取種子點5×5鄰域中的地面點作最小二乘平面擬合，以擬合平面作為初始地面，搜索標記圖像中與初始地面的高差小于高差閾值的點，作為誤分類的地面點。

四、試驗與分析

1.試驗結果

試驗采用ISPRS網站提供的點云數據，這些點云數據是在2003年由ISPRS Commission III提供給廣大學者的，用于比較不同濾波算法的試驗。試驗選取城市區(qū)域的四塊測試數據Site1-Site4的末次回波激光點作為原始數據，點間距為1.5 m，重采樣的格網間距為1 m。試驗區(qū)域散布著大、中、小型的不規(guī)則形狀的建筑物，同時有橋梁、隧道、鐵路、汽車以及與建筑物相鄰的樹木。在四塊測試數據中選取了九個樣本數據，對樣本數據進行了手工分類，將激光腳點精確分類為地面點和非地面點兩類。利用分類后的樣本數據，可以對濾波算法進行分類誤差的定量分析。

圖2(a)給出了應用多尺度數學形態(tài)學算法對重采樣后Site2的DSM進行濾波的結果，其中空白地區(qū)為分類后的非地面點。從圖2(a)可以看出，在A、C處由于地下通道與旁邊的路面存在高程突變，而在B處也由于地面不連貫變化，導致部分地面點明顯被誤分為非地面點。對濾波結果進行質量控制，搜索誤分類的地面點后得到最終分類結果，如圖2(c)所示。

圖3分別給出了Site1、Site3和Site4的分類后地面點重采樣生成的DEM。試驗結果表明該算法能夠有效地識別地面和地物點，地面細節(jié)信息得到保留;樹木被有效識別，所有建筑物被濾除掉。

圖2 Site2試驗結果

圖3 地面點重采樣的DEM效果圖

2.試驗分析

文獻［9］對不同濾波方法從算法上進行了試驗比較，表2給出了本文濾波方法與文獻［9］和文獻［4］中所涉及算法的總誤差對比。

表2 樣本的總誤差對比

由表2可以看出，在陡坡地帶的樣本11、復雜建筑物密集的地形高程變化較大的樣本23以及鐵路等高頻率的地形起伏樣本42地區(qū)，Axelsson的內插濾波方法的精度稍微高于本文和文獻［4］的形態(tài)學方法，因為從粗到細的處理方式能避免大地形起伏導致距離計算值較大而引起的誤判。而在地形平緩、建筑物密集的樣本12、樣本21、樣本31以及樣本41地區(qū)，本文方法明顯優(yōu)于內插方法和文獻［4］的形態(tài)學方法，尤其是數據樣本21和樣本41，這是由于多尺度形態(tài)學濾波方法采取由局部出發(fā)擴展到全局的策略，在平緩地區(qū)不容易導致誤判，尤其是在有數據空白和一個大型建筑物的樣本41，Axelsson的內插濾波方法對數據缺失比較敏感，而文獻［4］的形態(tài)學方法以高程較低值填補數據空白后，則會將高于數據空白地區(qū)的周邊地面點分類為非地面點，同時該方法的濾波窗口最大值小于建筑物的面積，導致部分建筑物點被分類成地面點，而本文以形態(tài)學濾波的結果作為初始DEM，通過進一步搜索誤分類的地面點明顯提高了分類的精度。

五、結束語

本文提出一種新的數學形態(tài)學濾波策略，該方法對點云數據進行合理的重采樣，并用多尺度數學形態(tài)學濾波后作為初始的DEM，通過搜索誤分類地面點的方式提高分類精度。該方法解決了數學形態(tài)學濾波算法中存在的兩個問題:局部地形不連續(xù)對分類的影響以及坡度或高差閾值的人工選取。試驗結果表明該方法能夠有效識別地面點和地物點，并且保留地形的細節(jié)信息。

［1］KRAUS K，PFEIFER N.Determination of Terrain Models in Wooded Areas with Airborne Laser Scanner Data［J］.ISPRS Journal of Photogrammetry and Remote Sensing，1998，53(4):193-203.

［2］PFEIFER N，STADLER P.Derivation of Digital Terrain Models in the SCOP++Environment［C］∥Proceedings of OEEPE Workshop on Airborne Laser Scanning and Interferometric SAR for Detailed Digital Elevation Models.Stockholm:［s.n.］，2001.

［3］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.

［4］CHEN QGONG P，BALDOCCHI D，et al.Filtering Airborne Laser Scanning Data with Morphological Methods［J］.Photogrammetric Engineering and Remote Sensing，2007，73(2):175-185.

［5］張熠斌，隋立春，曲佳，等.基于數學形態(tài)學算法的機載 LiDAR點云數據快速濾波［J］.測繪通報，2009(5):16-18.

［6］ELMQVIST M JUNGERT E，LANTZ F，et al.Terrain Modelling and Analysis Using Laser Scanner Data［J］.International Archives of Photogrammetry and Remote Sensing，2001，34(3/W4):219-224.

［7］BROVELLI M A，CANNATA M.Digital Terrain Model Reconstruction in Urban Areas from Airborne Laser Scanning Data:the Method and an Example of the Town of Pavia(Northern Italy)［J］.Computers ＆ Geosciences，2004，30(4):325-331.

［8］黃先鋒，李卉，王瀟，等.機載LiDAR數據濾波方法評述［J］.測繪學報，2009，38(5):465-469.

［9］SITHOLE G，VOSSELMAN G.Experimental Comparison of Filter Algorithms for Bare-Earth Extraction from Airborne Laser Scanning Point Clouds［J］.ISPRS Journal of Photogrammetry and Remote Sensing， 2004，59(1-2):85-101.

［10］穆超，余潔，許磊，等.基于高分辨率遙感影像的DSM建筑物點的提取研究［J］.武漢大學學報:信息科學版，2009，34(4):414-417.

A New Filtering Method for LiDAR Data Based on Mathematic Morphological Approach

LUO Yiping，JIANG Ting，WANG Xin，ZHANG Rui，LUO Sheng

0494-0911(2011)03-0015-05

P237

B

2010-02-23

羅伊萍(1982—)，女，湖南永州人，博士生，研究方向為LiDAR數據處理和航天遙感工程。