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

LIDAR點(diǎn)云數(shù)據(jù)全自動(dòng)濾波算法研究

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

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

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

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

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

0引言

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

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

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

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

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

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

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

圖1　點(diǎn)云濾波處理流程圖

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

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

1.3精簡(jiǎn)移動(dòng)最小二乘法擬合DEM

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

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

(1)

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

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

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

(2)

(3)

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

(4)

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

式 (4) 可寫為

(5)

式中：

(6)

(7)

(8)

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

(9)

式中：

A(x)=PTW(x)P.

(10)

B(x)=PTW(x).

(11)

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

(12)

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

(13)

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

(14)

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

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

圖2　MLS算法擬合DEM

2實(shí)驗(yàn)與應(yīng)用

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

圖3　精簡(jiǎn)MLS算法擬合DEM

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

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

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

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

圖6　激光點(diǎn)云高程與GPS高程數(shù)據(jù)對(duì)比

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

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

3結(jié)論

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

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

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

參考文獻(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]隋立春,張熠斌.基于改進(jìn)的數(shù)學(xué)形態(tài)學(xué)算法的LiDAR點(diǎn)云數(shù)據(jù)濾波[J].測(cè)繪學(xué)報(bào),2010, 39(4): 390-396.

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

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

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

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

[7]左志權(quán),張祖勛,張劍清.知識(shí)引導(dǎo)下的城區(qū)LiDAR點(diǎn)云高精度三角網(wǎng)漸進(jìn)濾波方法[J].測(cè)繪學(xué)報(bào),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]李樂(lè)林. 基于等高線族的機(jī)載LiDAR數(shù)據(jù)建筑物三維模型重建方法研究[D].武漢:武漢大學(xué)測(cè)繪遙感信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,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]趙波,邊馥苓. 面向移動(dòng)GIS的動(dòng)態(tài)四叉樹空間索引算法[J].計(jì)算機(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ìn)的移動(dòng)最小二乘法[J].力學(xué)季刊,2003,24(2):266-272.

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

[18]任紅萍, 程玉民,張武.改進(jìn)的移動(dòng)最小二乘插值法研究[J].工程數(shù)學(xué)學(xué)報(bào),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

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

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

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

中圖分類號(hào)：P237

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

doi:10.3969/j.issn.1671-6833.201504004

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