LIDAR點云數(shù)據(jù)全自動濾波算法研究

李　健, 方宏遠(yuǎn), 崔雅博, 范　濤

LIDAR點云數(shù)據(jù)全自動濾波算法研究

李健1, 方宏遠(yuǎn)1, 崔雅博2, 范濤3

(1.鄭州大學(xué) 水利與環(huán)境學(xué)院，河南 鄭州 450001；2.開封大學(xué) 實驗實訓(xùn)中心，河南 開封 475004；

3.河南省地質(zhì)環(huán)境監(jiān)測院，河南 鄭州 450001)

摘要：提出了一種基于移動最小二乘法的點云數(shù)據(jù)全自動濾波算法，該方法首先對LIDAR點云數(shù)據(jù)進行合理分塊，并建立分塊網(wǎng)格的動態(tài)四叉樹空間索引，便于數(shù)據(jù)操作和管理.對分塊網(wǎng)格中的點云數(shù)據(jù)利用精簡移動最小二乘法擬合出參考地形，將擬合得到的參考地形用于LIDAR點云高程閾值的迭代計算，將每次迭代前后高差小于閾值的點劃為地面點，其余點劃分為非地面點，迭代運算直至閾值滿足要求為止.實驗表明，精簡移動二乘法效率高，計算量小，并且精度高，適合點云數(shù)據(jù)DEM(digital elevation model)擬合，利用該算法對LIDAR點云數(shù)據(jù)進行濾波的速度快、精度高，能夠有效地識別地面點和非地面點，并保留地形的細(xì)節(jié)信息.

關(guān)鍵詞：點云數(shù)據(jù)；數(shù)字地面模型；濾波算法；動態(tài)四叉樹；移動最小二乘法

0引言

隨著激光技術(shù)的快速發(fā)展和完善，激光數(shù)據(jù)在眾多領(lǐng)域得到了廣泛的應(yīng)用.由于激光能在短時間內(nèi)獲得地物三維坐標(biāo)信息，并且數(shù)據(jù)量極大，故而如何快速從海量LIDAR點云數(shù)據(jù)中提取有用的信息是目前研究的熱點和難點[1].國內(nèi)外許多學(xué)者都對點云濾波進行了討論和研究，并且提出了許多濾波算法，包括基于數(shù)學(xué)形態(tài)學(xué)的濾波算法[2-3]、基于坡度的濾波算法[4-6]、基于TIN的漸進加密算法[7-8]等，都取得了一些研究成果，但其中還存在一些問題尚未解決.不管是機載LIDAR數(shù)據(jù)還是地面LIDAR數(shù)據(jù)大部分是基于激光點云中高程突變信息進行濾波，假定點云中高程低的點為地面點，高程較高的點為非地面點，由于系統(tǒng)誤差的存在，這種情況未必完全正確.另外一些濾波算法適用范圍有限.從上述問題可看出，提出一種簡單、快速、適用范圍廣、效率高的點云濾波算法是非常必要的[9].

由于激光點數(shù)據(jù)量大，并且點云數(shù)據(jù)的不規(guī)則、散亂復(fù)雜等性質(zhì)決定了點云數(shù)據(jù)處理工作的復(fù)雜困難[10-11].針對LIDAR點云數(shù)據(jù)的特點，筆者提出了先將點云數(shù)據(jù)進行網(wǎng)格分塊，保證點云數(shù)據(jù)的原始性，減少單次數(shù)據(jù)處理量.對分塊數(shù)據(jù)建立空間索引，提高點云數(shù)據(jù)處理的效率.

1關(guān)鍵技術(shù)與算法

1.1LIDAR點云數(shù)據(jù)的濾波流程

將海量激光點云分塊并建立相應(yīng)的空間索引關(guān)系后，進行地物的自動過濾處理，濾波要考慮當(dāng)前點所在的網(wǎng)格，并對其進行計算，每次計算的結(jié)果再以索引的方式動態(tài)存儲，作為下次迭代計算的基礎(chǔ)數(shù)據(jù)，具體濾波流程如圖1所示.

1.2點云數(shù)據(jù)的網(wǎng)格分塊與動態(tài)四叉樹空間索引

為了進行激光點云的海量數(shù)據(jù)管理、處理與顯示，對激光點云分塊處理顯得尤為必要.分塊的大小直接影響到數(shù)據(jù)處理層次及深度，相應(yīng)地影響算法的效率.分塊越小，分割越細(xì)，效率就越低，其合并的區(qū)域相對增大，數(shù)據(jù)的壓縮比就越高；反之，效率就越高，而壓縮比相對降低.最小格網(wǎng)大小的選擇應(yīng)是最小采樣間距的整數(shù)倍，具體數(shù)值的確定取決于被測對象的復(fù)雜度、儀器的最小采樣間距以及期望的數(shù)據(jù)壓縮比.

圖1　點云濾波處理流程圖

為了高效地管理和存儲分塊網(wǎng)格及網(wǎng)格內(nèi)的激光點云數(shù)據(jù)，需要對所有存在點的網(wǎng)格與點云之間建立索引關(guān)系，在對點云進行處理和過濾時，只需要考慮點所在的網(wǎng)格即可，這里采用四叉樹結(jié)構(gòu)[12].

由于激光點云的分布不均勻和邊界形狀極其不規(guī)則，為了克服常規(guī)四叉樹空間索引結(jié)構(gòu)中的問題，筆者采用空間動態(tài)四叉樹的方法對分塊網(wǎng)格建立索引關(guān)系.其算法要點為：在開始建立四叉樹時，不需要事先確定工作區(qū)域的范圍，只須把要插入空間數(shù)據(jù)庫的第一個對象的MBR(minimum bounding rectangle)中心作為四叉樹的頂點，隨著數(shù)據(jù)處理工作的進行，對作業(yè)空間進行分解.

1.3精簡移動最小二乘法擬合DEM

移動最小二乘法提供一種較高次數(shù)的多項式逼近方式對散亂點云進行曲面擬合[13]，并且要求擬合函數(shù)在各個節(jié)點處的誤差的平方和最小，能夠保證比較高的精度.用該算法擬合的曲面比較平滑，與實際曲面近似[14-15].但是移動最小二乘算法比較復(fù)雜，運算效率不高，如果點云數(shù)據(jù)量比較大，則處理比較困難，因此需要對該方法進行精簡，在保證精度前提下提高其運算效率.精簡移動最小二乘法采用帶權(quán)的正交函數(shù)作為基函數(shù)[16-18]，在求解系數(shù)矩陣時可只考慮對角線元素，不用求逆矩陣，減少了運算量、提高了運算效率，同時也提高了精度，適合于數(shù)據(jù)量比較大的點云數(shù)據(jù)擬合.

最小二乘擬合函數(shù)：

(1)

式中：αi(X)為待求系數(shù)(i=1,2,...，n)，是坐標(biāo)X的函數(shù)；pi(X)稱為基函數(shù)，它是一個k階完備的多項式；n是基函數(shù)的項數(shù).

對于二維空間擬合，基函數(shù)PT(X)的形式為

線性基：PT=(1,x1,x2).

(2)

(3)

為了使得曲面在點x局部擬合最優(yōu)，定義誤差方程：

(4)

式中：w(X-Xi)為權(quán)重函數(shù);m為計算域內(nèi)節(jié)點總數(shù).

式 (4) 可寫為

(5)

式中：

(6)

(7)

(8)

待求系數(shù)α(x)通過求取目前函數(shù)J的極值獲取(見式(9)).

(9)

式中：

A(x)=PTW(x)P.

(10)

B(x)=PTW(x).

(11)

為了避免式(9)在求解過程中出現(xiàn)病態(tài)和奇異，這時如果假定對于點集{x}和權(quán){wiJB(i=1,2,...n)},若存在一組函數(shù)pi(x)(i=1,2,...,n)滿足：

(12)

則稱pi(x)(i=1,2,…n)為點集{x}和權(quán){wi}的正交函數(shù)集，那么移動最小二乘可變?yōu)椋?/p>

(13)

那么很容易就解出αi(x)：

(14)

待定系數(shù)αi(x)可以通過求解式(14)得到，并且αi(x)的求解過程避免了求矩陣的逆，避免了求解病態(tài)方程，計算效率得到了提高.

為了驗證該算法，利用掃描點云數(shù)據(jù)中2 658個點來檢測對比精簡移動最小二乘算法與移動最小二乘算法擬合的DEM，圖2為在MATLAB中使用移動最小二乘擬合的DEM，耗時156s，從擬合DEM數(shù)據(jù)可以看出構(gòu)建的DEM精度較低，而使用IDL(interactivedatalanguage)編程開發(fā)的精簡移動最小二乘法擬合的DEM(如圖3所示)，只用了30s，時間減少了126s，效率得到提高，從效果來看,擬合DEM比較平滑，與實際地形較接近.

圖2　MLS算法擬合DEM

2實驗與應(yīng)用

筆者主要通過IDL編程實現(xiàn)上述算法， IDL是面向矩陣、語法簡單的第四代可視化語言，可以應(yīng)用于任何領(lǐng)域的三維數(shù)據(jù)可視化、數(shù)值計算、三維圖形建模、科學(xué)數(shù)據(jù)讀取等功能中.IDL是完全面向矩陣的，它具有快速分析超大規(guī)模數(shù)據(jù)的能力，這在海量點云處理方面具有獨特的優(yōu)勢[19].

圖3　精簡MLS算法擬合DEM

為了驗證本算法的精度，將濾波算法得到的地面點云數(shù)據(jù)進行抽稀，與本區(qū)域內(nèi)采集GPS測量數(shù)據(jù)進行對比，來評定該算法的精度.通過十站掃描采集約一平方公里的點云，原始點云數(shù)據(jù)約2.5G(如圖4)，獲得本區(qū)域內(nèi)200個GPS點數(shù)據(jù).將掃描點云數(shù)據(jù)拼接后使用筆者提出的算法對點云數(shù)據(jù)進行濾波，得到地面點云數(shù)據(jù)并進行抽稀，如圖5所示.將抽稀后的地面點云數(shù)據(jù)在CAD中進行高程展點，并將200個GPS點高程數(shù)據(jù)加入CAD中，將兩組數(shù)據(jù)進行統(tǒng)計對比，如圖6所示.淺色為激光點云高程數(shù)據(jù)，深色為GPS高程數(shù)據(jù).在這里由于獲得的地面點云數(shù)據(jù)與采集的GPS點坐標(biāo)不是同名點，因此只能根據(jù)離GPS點最近或者附近的地面激光點的高程進行統(tǒng)計，存在一些高程異常點，高程點誤差范圍見表1.

圖4　原始掃描激光點云數(shù)據(jù)

利用筆者提出的濾波算法對十站掃描點云數(shù)據(jù)濾波，過濾后地面激光點云個數(shù)約為8 145 631個.在濾波算法中分塊網(wǎng)格大小為2 m，迭代3次，普通筆記本電腦過濾時間約1 min將過濾后激光點云數(shù)據(jù)構(gòu)建間距2 m的DEM網(wǎng)格如圖7所示.

圖5　濾波抽稀10 m點云數(shù)據(jù)

圖6　激光點云高程與GPS高程數(shù)據(jù)對比

絕對值高程誤差/cm激光點云數(shù)目<1542～5776～102811～152216～2014高程異常點7

圖7　過濾后的激光點云數(shù)據(jù)生成DEM數(shù)據(jù)

3結(jié)論

(1)根據(jù)LIDAR點云數(shù)據(jù)的特點，提出了基于網(wǎng)格分塊的過濾方法，并建立了分塊網(wǎng)格的動態(tài)四叉樹的空間索引方法，解決了海量激光點云數(shù)據(jù)的存儲和管理，進而實現(xiàn)了海量激光點云數(shù)據(jù)的快速處理.并且使用動態(tài)最小二乘法擬合DEM對點云數(shù)據(jù)進行濾波，擬合平面比較平滑，與實際地形接近.

(2)從分類得到點云數(shù)據(jù)與GPS高程點精度對比可以看出，筆者所采用的濾波算法精度較高，能夠快速有效地從激光點云數(shù)據(jù)中識別地面點和非地面點，實用性強.

(3)從文中兩個案例可以看出本算法適用于機載LIDAR和地面LIDAR的濾波處理，且能夠高效地對海量點云數(shù)據(jù)進行濾波處理.但對一些特殊地形(如懸崖)仍然存在一定的問題，需要采用轉(zhuǎn)換坐標(biāo)系或進行局部過濾等處理方法.

參考文獻(xiàn)：

[1]AXELSSON P. Processing of laser scanner data algorithms and applications [J]. ISPRS international journal of photogrammetry and remote sensing, 1999,54(2/3):138-147.

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

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

[4]楊應(yīng),蘇國中,周梅.影像分類信息支持的LiDAR點云數(shù)據(jù)濾波方法研究[J].武漢大學(xué)學(xué)報(信息科學(xué)版),2010, 35(12): 1453-1456.

[5]袁楓,張繼賢,張力,等.結(jié)合強度信息的LiDAR數(shù)據(jù)濾波方法[J].測繪科學(xué),2010,35(5): 39-41.

[6]劉凱斯.機載激光LiDAR點云數(shù)據(jù)濾波和分類算法研究[D].北京：首都師范大學(xué)資源環(huán)境與旅游學(xué)院,2014.

[7]左志權(quán),張祖勛,張劍清.知識引導(dǎo)下的城區(qū)LiDAR點云高精度三角網(wǎng)漸進濾波方法[J].測繪學(xué)報,2012, 41(2): 246-251.

[8]AXELSSON P. DEM Generation from laser scanner data using adaptive TIN models [J]. International archives of photogrammetry and remote sensing, 2000,XXXIII (B4):110-117.

[9]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, 5(4):193-203.

[10]李樂林. 基于等高線族的機載LiDAR數(shù)據(jù)建筑物三維模型重建方法研究[D].武漢:武漢大學(xué)測繪遙感信息工程國家重點實驗室,2012.

[11]KHAMARRUL A R, MICHELE S C J， VAN W., et al. Generating an optimal DTM from airborne laser scanning data for landslide mapping in a tropical forest environment [J]. Geomorphology,2013,190:112-125.

[12]趙波,邊馥苓. 面向移動GIS的動態(tài)四叉樹空間索引算法[J].計算機工程,2007,33(15):86-88.

[13]LANCASTER P. SALKAUSKAS K. Surfaces generated by moving least squares methods [J]. Mathematics of computation,1981,37:141-158.

[14]TSCHKO T B, Lu Y Y, GU L. Elements-free galerkin methods [J]. International journal for numerical methods in engineering,1994,37:229-256.

[15]LEVIN D. The Approximation power of moving least squares [J]. Mathematics of computation,1998, 67:1517～1531.

[16]陳美娟,程玉民.改進的移動最小二乘法[J].力學(xué)季刊,2003,24(2):266-272.

[17]倪慧,李重,宋紅星,等.帶插值條件的最小移動二乘曲線擬合[J].浙江理工大學(xué)學(xué)報,2011,28(1):135-139.

[18]任紅萍, 程玉民,張武.改進的移動最小二乘插值法研究[J].工程數(shù)學(xué)學(xué)報,2010,27(6):1021-1029.

[19]韓培友. IDL可視化分析與應(yīng)用[M].西安:西北工業(yè)大學(xué)出版社,2006.

An Automatic Point Clouds Filtering Algorithm Based on Grid Partition and Simplified Moving Least Squares

LI Jian1, FANG Hongyuan1, CUI Yabo2, FAN Tao3

(1 College of Water Conservancy & Environmental Engineering, Zhengzhou University, Zhengzhou 450001， China; 2 Training Center，Kaifeng University, Kaifeng 475004， China; 3 Henan Geological Environment Information, Zhengzhou 450001， China)

Abstract:An automatic point clouds filtering algorithm is presented on the basis of Grid Partition using Dynamic Quad Trees and reference surface fitted by Moving Least Squares. The filtering processing contains three major steps: Firstly, it gives the LIDAR point clouds reasonable grid partitions and establishes the corresponding dynamic quad trees spatial indices. Secondly, the points in the partitioned grids are utilized to fit a DEM reference plane using moving least squares technology. Finally, the elevation threshold is setup to separate ground points from those non-ground ones who are positioned above the reference plane and have a distance exceeding the threshold value to the plane. The aforementioned steps have to be repeated on the obtained ground points with gradually decreasing thresholds and grid size until desired precision is achieved. The experiments show that simplified moving least squares is high efficiency, small amount of calculation and high precision DEM data for point cloud fitting, and the filtering algorithm has high precision and can effectively identify ground points and non-ground ones without losing the detailed information of topography.

Key words:point clouds data; DEM; filtering algorithm; dynamic quad trees; moving least square

收稿日期:2015-04-02；

修訂日期：2015-10-28

基金項目：國家自然科學(xué)基金青年基金資助項目(41404096)；河南省教育廳基金資助項目(14A420002,15A420002)

作者簡介：李健(1983—)，男，河南孟州人，鄭州大學(xué)講師，博士，主要從事點云數(shù)據(jù)處理，E-mail:jianli@zzu.edu.cn.

文章編號:1671-6833(2016)01-0092-05

中圖分類號：P237

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

doi:10.3969/j.issn.1671-6833.201504004

引用本文：李健, 方宏遠(yuǎn), 崔雅博，等.LIDAR點云數(shù)據(jù)全自動濾波算法研究[J].鄭州大學(xué)學(xué)報(工學(xué)版)，2016,37(1):92-96.