孫　蒙　顧和和

1　中國礦業(yè)大學(xué)環(huán)境與測繪學(xué)院，徐州市大學(xué)路1號，221116

基于微分形態(tài)學(xué)斷面的機(jī)載LiDAR數(shù)據(jù)濾波新方法

孫蒙1顧和和1

1中國礦業(yè)大學(xué)環(huán)境與測繪學(xué)院，徐州市大學(xué)路1號，221116

摘要：提出一種基于微分形態(tài)學(xué)斷面的機(jī)載LiDAR點云數(shù)據(jù)濾波新方法。該方法由點云數(shù)據(jù)構(gòu)建規(guī)則格網(wǎng)，去除粗差點；對構(gòu)建的每一個格網(wǎng)進(jìn)行多尺度分解，獲取初始地面點及地物點；分別利用曲面逼近及微分形態(tài)學(xué)斷面構(gòu)建DTM(digital terrain model)，通過閾值函數(shù)判別二者之間的殘差，確定最終地面點。使用國際攝影測量與遙感學(xué)會提供的測試數(shù)據(jù)進(jìn)行實驗，并與8種經(jīng)典濾波算法比較分析，表明該方法能夠有效去除地物點和保留地面點，并降低總誤差。

關(guān)鍵詞：微分形態(tài)學(xué)斷面；多尺度分解；數(shù)學(xué)形態(tài)學(xué)；濾波；精度評定

LiDAR可低成本、快速、大范圍、實時獲取高精度三維點云[1]，被廣泛應(yīng)用到構(gòu)建數(shù)字高程模(digital elevation model,DEM)及獲取各種圖像產(chǎn)品[2-3]。由于受建筑物、植被和車輛等地物的影響，原始三維點云中含有非地面點，為此需對其進(jìn)行濾波處理[4]。經(jīng)典濾波算法有基于內(nèi)插[5-6]、坡度[7-8]和形態(tài)學(xué)的濾波[9]等?；趦?nèi)插的濾波方法認(rèn)定地球表面為連續(xù)變化的曲面，通過一個起始DEM逐步從備選點篩選并內(nèi)插加密實現(xiàn)分類目的，易受初始DEM的影響且誤差會隨迭代過程積累。基于坡度的濾波方法需量測兩點之間的高差或坡度，通過預(yù)設(shè)的坡度閾值完成分類，易受地形變化及閾值選取的影響，濾波方法自適應(yīng)性不強(qiáng)?；谛螒B(tài)學(xué)的濾波方法僅考慮窗口內(nèi)的地形特征，易受局部地形的影響，且濾波過程中參數(shù)需要用戶根據(jù)不同地形設(shè)置，自適應(yīng)性不強(qiáng)。針對數(shù)學(xué)形態(tài)學(xué)濾波算法存在的問題，本文提出一種基于微分形態(tài)學(xué)斷面(differential morphological profiles, DMPs)的機(jī)載LiDAR數(shù)據(jù)濾波新方法。使用國際攝影測量與遙感學(xué)會(ISPRS)提供的樣本數(shù)據(jù)進(jìn)行實驗表明，該算法可以得到較好的濾波結(jié)果并使總誤差較小。

1DMPs濾波原理和步驟

為克服地形起伏的影響，考慮到濾波區(qū)域目標(biāo)尺寸形狀差異，本文在預(yù)估的粗糙近似曲面上基于微分形態(tài)學(xué)斷面進(jìn)行多尺度分解。算法主要流程分為點云格網(wǎng)化、DMPs構(gòu)建和點云濾波3步(圖1)。

圖1　微分形態(tài)學(xué)斷面濾波流程圖Fig.1　The framework of the DMPs filter

1.1點云格網(wǎng)化

設(shè)原始LiDAR點云集合為L={li}，則預(yù)處理后構(gòu)建的格網(wǎng)表示為E={gi}。其中E的范圍由L的邊界確定；gi不僅包含格網(wǎng)內(nèi)所有的原始點云信息，同時特指該格網(wǎng)內(nèi)最低點的三維坐標(biāo)信息。規(guī)則格網(wǎng)化的步驟如下。

1)確定格網(wǎng)分辨率Rg。Rg取原始LiDAR點云密度DL的倒數(shù)，即

(1)

2)確定單元格值。若格網(wǎng)化后單元格內(nèi)存在的原始最低點為p，則該單元格的值為p點高程值g[p]；對于空單元格，其值g[p*]由鄰近的單元格內(nèi)插得到。本文采用反距離加權(quán)法[10]，即

(2)

由于多路徑效應(yīng)及儀器誤差的影響，原始點云中存在粗差點，本文通過構(gòu)建格網(wǎng)剔除低位粗差點(遠(yuǎn)低于地面的點)。

1.2DMPs構(gòu)建

通過多尺度格網(wǎng)分解方法對格網(wǎng)g中含有的特征進(jìn)行映射分析，其決定了點云濾波后生成的DTM的精度。首先將預(yù)處理后獲得的格網(wǎng)g分解為低頻帶g↓和高頻帶g↑，g↓通過高斯核函數(shù)G對格網(wǎng)g進(jìn)行卷積運(yùn)算得到，該卷積運(yùn)算可剔除格網(wǎng)g中大部分的非地面特征；g↑為獲取g↓后格網(wǎng)g的剩余部分，即

(3)

式中，*為卷積運(yùn)算，G為高斯核函數(shù)。由于g↑為濾除大部分非地面特征的近似地面，接下來的處理主要是針對g↑。

設(shè)rωs(g↑)為對g↑進(jìn)行數(shù)學(xué)形態(tài)學(xué)的開運(yùn)算，ωs是大小為s的圓盤結(jié)構(gòu)元素。通過一系列呈漸進(jìn)變化的濾波窗口，采用形態(tài)學(xué)開運(yùn)算粒度測試技術(shù)對g↑進(jìn)行多尺度分解。這些濾波窗口集合為s={0,1,2,…,S}，其中S為所包含目標(biāo)的最大尺寸。故格網(wǎng)g的微分形態(tài)學(xué)斷面表示為[11]：

(4)

式中，i∈[1,I]，I為集合s的勢減1。如圖2所示，對于每一個特定點p，Δ(g↑)為g↑的分解，Δ(g↑)[p]為該處長度為I的響應(yīng)矢量，該矢量包含與每一個rωsi有關(guān)的高差值。

圖2　格網(wǎng)分解過程Fig.2　Decomposition of a grid

對Δ(g↑)進(jìn)行特征映射，可獲得g′和g*。g′是表征每一點p處所獲得的Δ(g↑)的最大相應(yīng)函數(shù)；g*是在p點得到最大響應(yīng)值時求取最小濾波窗口的函數(shù)。定義函數(shù)θws(g↑)用于求取p點處的g′和g*，即

(5)

式中，∨為上確界，∨為下確界。g′描述了目標(biāo)特征與周邊地形之間的高差，g*則是對目標(biāo)特征大小的估計。

1.3點云濾波

盡管通過設(shè)置g′的閾值可以識別和移除不滿足假設(shè)的目標(biāo)，但是一些小塊低矮地物會對g′的判別有影響，故應(yīng)考慮將g*也作為判別準(zhǔn)則。本文采用線性閾值函數(shù)T(g*)[p]，定義為：

(6)

式中，k和n為用戶自定義常量，描述了特征尺寸與殘差的相關(guān)系數(shù)，與地形環(huán)境及地物復(fù)雜程度有關(guān)。地面點集SG和非地面點集SNG表示為：

(7)

(8)

對于獲取的地面點集SG，采用式(2)內(nèi)插生成DTM。為此，可在閾值函數(shù)T(g*[p])基礎(chǔ)上，采用外置形態(tài)學(xué)梯度閾值函數(shù)?DTM,即?DTM=δωi(DTM)-DTM。其中，δωi(DTM)為用最小尺寸的結(jié)構(gòu)元素(如3×3)進(jìn)行形態(tài)學(xué)膨脹運(yùn)算得到的結(jié)果。則一系列地面點集LG可通過下式獲?。?/p>

(9)

式中，dh(li)為原始點云L中某一點li與相應(yīng)DTM之間的絕對高程差，b為用于定義DTM與地面點之間的最大高差。

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

2.1實驗數(shù)據(jù)

本文用來測試的數(shù)據(jù)為ISPRS第Ⅲ小組2003年采集的一組數(shù)據(jù)集。該數(shù)據(jù)集位于德國南部的Vaihingen/Enz和Stuttgart市中心，包含多種特征目標(biāo)(植被、建筑物、道路、鐵路、河流、橋梁、電力線、水面等)，由Optech公司的ALTM機(jī)載LiDAR測量系統(tǒng)獲取，其首末次回波皆可用。數(shù)據(jù)集由15個樣本組成，其中包含平均點間距為1～1.5 m的城區(qū)(samp11～42)及平均點間距為2～3.5 m的山區(qū)(samp51～71)。每個樣本都由測試數(shù)據(jù)與參考數(shù)據(jù)2個部分組成，參考數(shù)據(jù)通過半自動濾波及結(jié)合影像的人工識別進(jìn)行地面點與非地面點的精確區(qū)分。

2.2實驗參數(shù)選取

本文方法需要用戶自定義4個常量，而其中常量S一般指濾波區(qū)域內(nèi)的最大地物尺寸，決定著算法的運(yùn)行效率；常量b為高差閾值，與初始地面點生成的DTM精度有關(guān)；常量k、n的取值受地形環(huán)境的復(fù)雜程度影響，本文采取單一變量法，k與n在0～1之間呈線性遞增取值，通過比較濾波后生成的DTM確定最優(yōu)值。

2.3實驗結(jié)果分析

通過C++編程實現(xiàn)該算法，并對這15個樣本進(jìn)行濾波處理。實驗所用最優(yōu)參數(shù)、濾波所用時間及結(jié)果見表1。可以看出，濾波過程所耗時間皆在0.3 s以內(nèi)，算法每s可處理百余萬點云，說明該算法運(yùn)行效率高。由濾波結(jié)果中存在的Ⅱ類誤差可以看出，無論在山區(qū)還是城區(qū)，其結(jié)果都能滿足實際需求，證明該算法濾波自適應(yīng)強(qiáng)。

表1　實驗參數(shù)值

2.3.1定性比較

通過表1，分別從城區(qū)和山區(qū)樣本中選擇總誤差最大的樣本samp11及內(nèi)插生成的DTM平均誤差最大的樣本samp53進(jìn)行定性比較，圖3顯示了每一個樣本的原始DSM、濾波后DTM以及濾波后Ⅰ類和Ⅱ類誤差的空間分布。從圖3(a)可以看出，由于陡坡上建筑物多、目標(biāo)復(fù)雜以及裸露梯田的存在，導(dǎo)致建筑多的地方濾波后Ⅰ類誤差明顯增大。Ⅱ類誤差主要受低矮植被的影響(圖3(c))。從圖3(d)可以看出，由于存在粗糙狹窄的褶皺與陡崖，導(dǎo)致濾波后Ⅰ類誤差主要分布在陡崖處及褶皺地帶(圖3(f))。

圖3　樣本濾波前后的DTM及Ⅰ類、Ⅱ類誤差空間分布Fig.3　The DTM before and after filtering and type Ⅰ errors and type Ⅱ errors of samples

2.3.2定量分析

將15個樣本經(jīng)DMPs濾波后與ISPRS提供的8種經(jīng)典濾波算法進(jìn)行對比分析，結(jié)果如表2所示。表中黑色加粗字體即為每一行的最小值，表示相對應(yīng)的濾波方法在該樣本中濾波效果最好。通過表2可以看出，本文提出的算法在樣本samp11、samp12、samp21、samp22、samp24、samp31、samp41、samp42、samp53、samp54、samp61及總平均誤差上優(yōu)于8種經(jīng)典算法。因為以上這些樣本區(qū)域內(nèi)無復(fù)雜建筑物、地表連續(xù)變化，故通過DMPs構(gòu)建生成的初始DTM精度高，能濾除大部分的非地面點，保證后續(xù)濾波精度。在樣本samp11、samp41及samp53中雖然總誤差比較小，但是由于非地面點數(shù)遠(yuǎn)小于地面點數(shù)，錯分類的非地面點會使Ⅱ類誤差明顯增大。在剩余樣本中本文算法濾波結(jié)果與最優(yōu)濾波算法Axelsson結(jié)果相差不大。剩余的樣本特征為復(fù)雜建筑物或低矮植被，如在樣本samp51與samp52中，由于區(qū)域內(nèi)低矮的植被點較多，根據(jù)生成的形態(tài)學(xué)梯度閾值函數(shù)進(jìn)行判別時，結(jié)構(gòu)元素內(nèi)最低的植被點會被保留下來，使總誤差及Ⅱ類誤差增大。而Axelsson算法假設(shè)地形局部區(qū)域是平坦的，然后選擇局部區(qū)域內(nèi)的最低點作為種子點來構(gòu)建初始稀疏三角網(wǎng)，并依據(jù)三角網(wǎng)內(nèi)隨機(jī)點距離三角形最近頂點的高度夾角和垂距來判斷該點是否為地面點。如果是地面點，則加密TIN模型，依次迭代直到?jīng)]有新的地面點為止。該算法比較適用于建筑物密集的城區(qū)以及地形不連續(xù)的區(qū)域，故在剩余樣本中該方法的濾波結(jié)果較好。

表2　MHC與8種經(jīng)典濾波算法總誤差比較

3結(jié)語

本文借鑒微分形態(tài)學(xué)斷面在遙感影像處理方面的應(yīng)用，提出了一種新的形態(tài)學(xué)濾波策略。該方法對點云數(shù)據(jù)進(jìn)行規(guī)則格網(wǎng)化、濾除粗差點；對構(gòu)建的每一個格網(wǎng)進(jìn)行多尺度分解，獲取初始地面點及地物點；分別利用曲面逼近及微分形態(tài)學(xué)斷面構(gòu)建DTM，通過閾值函數(shù)判別二者之間的殘差，確定最終地面點。該方法解決了數(shù)學(xué)形態(tài)學(xué)濾波算法中存在的兩個問題:局部地形不連續(xù)對分類的影響以及坡度或高差閾值的人工選取。但是在建筑物復(fù)雜的城區(qū)以及低矮植被較多的區(qū)域，該方法不能很好地克服因形態(tài)學(xué)梯度誤差引起的誤分類。比較實驗結(jié)果可知，該方法能快速有效地識別地面點和地物點，并且保留地形的細(xì)節(jié)信息，滿足精度需求。

參考文獻(xiàn)

[1]Yun J, Gao C X, Zhu S L, et al. High-Peak-Power,Single-Mode,Nanosecond Pulsed,All-Fiber Laser for High Resolution 3D Imaging LiDAR System[J]. Chinese Optics Letters, 2012, 10(12):39-41

[2]Lai X D,Han Z J. Hige Fidelity DEM Generation Based on LiDAR Data[J].Journal of Computers,2012,7(8):2 071-2 077

[3]范士俊，張愛武，胡少興，等．基于隨機(jī)森林的機(jī)載激光全波形點云數(shù)據(jù)分類方法[J].中國激光，2013，40(9)：216-222(Fan Shijun,Zhang Aiwu,Hu Shaoxing, et al. A Method of Classification for Aribone Full Waveform LiDAR Data Based on Random Forest[J].Chinese Journal of Lasers, 2013,40(9)：216-222)

[4]劉凱斯．機(jī)載激光LiDAR點云數(shù)據(jù)濾波和分類算法研究[D]．北京：首都師范大學(xué)， 2014(Liu Kaisi. The Research of Airborne Laser LiDAR Point Cloud Data Filtering and Classification Algorithm[D]. Beijing: Capital Normal University, 2014)

[5]Vosselman G. Slope Based Filtering of Laser Altimetry Data[J]. International Archives of Photogrammetry and Remote Sensing and Spatial Information Sciences, 2000,33(B3): 935-942

[6]Mongus D, Zalik B.Parameter-Free Ground Filtering of LiDAR Data for Automatic DTM Generation[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2012,67(1):1-12

[7]Filin S, Pfeifer N. Segmentation of Airborne Laser Scanning Data Using a Slope Adaptive Neighborhood[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2006,60(2),71-80

[8]Tovari D,Pfeifer N. Segmentation Based Robust Interpolation a New Approach to Laser Data Filtering[J].International Archives of Photogrammetry and Remote Sensing and Spatial Information Sciences,2005,36(19),79-84

[9]羅伊萍,姜挺,龔志輝. 基于自適應(yīng)和多尺度數(shù)學(xué)形態(tài)學(xué)的點云數(shù)據(jù)濾波方法[J]. 測繪科學(xué)技術(shù)學(xué)報,2009,26(6):426-429(Luo Yiping, Jiang Ting ,Gong Zhihui.An Adaptive and Muli-Scale Mathematic Morph-ological Filter for Point Cloud Data Filtering[J]. Journal of Geomatics Science and Technology,2009,26(6):426-429)

[10]Chaplot V, Darboux F, Bourennane H, et al. Accuracy of Interpolation Techniques for the Derivation of Digital Elevation Models in Relation to Landform Types and Data Density[J].Geomorphology,2006,77(1-2):126-141

[11]Pesaresi M, Benediktsson J A. A New Approach for the Morphological Segmentation of High-Resolution SatelliteImagery[J]. IEEE Transactions on Geoscience and Remote Sensing,2001, 39(2): 309-320

Foundation support:Commonweal Research Project on Land and Resources,No.201111013; Priority Academic Program Development of Jiangsu Higher Education Institutions,No.SZBF2011-6-B35.

About the first author:SUN Meng, postgraduate, majors in airborne LiDAR point cloud data processing，E-mail：519282720@qq.com.

A New Filtering Method for Aribone LiDAR Data Based on Differential Morphological Profiles

SUNMeng1GUHehe1

1School of Environment Science and Spatial Informatics, China University of Mining and Technology,1 Daxue Road, Xuzhou 221116, China

Abstract:This paper proposes a new filtering method for LiDAR Data. The proposed approach constructs the connectivity of a grid over the LiDAR point-cloud in order to perform multi-scale data decomposition. This is realized by forming a top-hat scale-space using differential morphological profiles (DMPs) on points’ residuals from the approximated surface. The International Society for Photogrammetry and Remote Sensing (ISPRS) reference dataset is used to test the method. The experimental results show that the proposed method can effectively remove non-ground points, keep the ground points, and is effective at minimizing total error rates.

Key words:differential morphological profiles；multi-scale data decomposition；mathematical morphology；filtering；accuracy assessment

收稿日期：2015-07-08

第一作者簡介：孫蒙，碩士生，主要研究方向為機(jī)載LiDAR點云數(shù)據(jù)處理，E-mail：519282720@qq.com。 通訊作者：顧和和，教授，主要研究方向為測量數(shù)據(jù)處理、地籍測量，E-mail：guhh@263.net。

DOI：10.14075/j.jgg.2016.07.007

文章編號：1671-5942(2016)07-0591-04

中圖分類號：P237

文獻(xiàn)標(biāo)識碼：A

Corresponding author:GU Hehe, professor, majors in surveying data processing and cadastral surveying, E-mail：guhh@263.net.

項目來源：國土資源公益性行業(yè)科研專項(201111013)；江蘇省高校優(yōu)勢學(xué)科建設(shè)工程(SZBF2011-6-B35)。