結(jié)合CSF和TIN的機(jī)載LiDAR點(diǎn)云濾波算法

崔 浩,高 飛,余 敏,葉周潤

(合肥工業(yè)大學(xué) 土木與水利工程學(xué)院,安徽 合肥 230009)

在工程建設(shè)、道路設(shè)計(jì)、滑坡監(jiān)測(cè)、土地覆蓋分類、森林管理和道路設(shè)計(jì)中,高分辨率數(shù)字地形模型(digital terrain model,DTM)有著至關(guān)重要的作用。光探測(cè)與測(cè)距技術(shù)(light detection and ranging,LiDAR)又稱激光雷達(dá)，是一種有效的收集上空三維點(diǎn)云的方法,并廣泛用于生成DTM。為了獲取高精度DTM,需要將地面和非地面測(cè)量數(shù)據(jù)進(jìn)行分離,這是一個(gè)濾波過程。目前地面濾波算法主要有基于坡度的方法、基于數(shù)學(xué)形態(tài)的方法和基于地表的方法等。

基于坡度的方法通常假設(shè)在一個(gè)鄰近區(qū)域內(nèi),地形的坡度變化通常是漸進(jìn)的,而建筑物或樹木與地面之間的坡度變化很大。基于這個(gè)假設(shè)文獻(xiàn) [1]開發(fā)一種基于坡度的濾波算法,利用LiDAR點(diǎn)與其相鄰點(diǎn)之間坡度的差值來進(jìn)行濾波。該方法在復(fù)雜地形如陡坡、斷裂地形處容易出現(xiàn)濾波錯(cuò)誤,具有一定的局限性[2]。

基于數(shù)學(xué)形態(tài)的方法是通過數(shù)學(xué)形態(tài)學(xué)來剔除非地面LiDAR點(diǎn)。對(duì)于這種濾波方法,選擇最佳窗口大小是至關(guān)重要的[3-6]。較小尺寸的窗口可以有效地濾除較小的物體,但會(huì)保留較大的建筑物在地面點(diǎn)中,而大的窗口尺寸往往會(huì)濾除平滑的地形細(xì)節(jié),如山峰、山脊和懸崖,故通常需要利用研究區(qū)域的額外先驗(yàn)知識(shí)來確定合適的窗口大小。

基于地表的方法主要是建立地形模型,然后通過迭代運(yùn)算近似得到原始地形,其中最經(jīng)典的就是文獻(xiàn)[7]提出的漸進(jìn)式不規(guī)則三角網(wǎng)(triangular irregular network,TIN)算法。基于TIN濾波的可靠性和準(zhǔn)確性已經(jīng)被大量的研究所證實(shí),但其在處理不連續(xù)地形時(shí)濾波效果較差[8-10]。

針對(duì)上述問題,許多研究者提出改進(jìn)的方法。文獻(xiàn)[11]利用數(shù)字形態(tài)學(xué)算法提取種子點(diǎn);文獻(xiàn)[12]利用均方差、點(diǎn)密度和小格網(wǎng)高程統(tǒng)計(jì)數(shù)據(jù)來選取種子點(diǎn);文獻(xiàn)[13]采用鄰域卷積和進(jìn)行種子點(diǎn)的判斷與虛擬種子點(diǎn)的構(gòu)造;文獻(xiàn)[14]采用多核并行處理技術(shù)加快三角網(wǎng)漸進(jìn)加密濾波方法的算法效率。

單一的算法往往都具有一定的局限性,而絕大多數(shù)的改進(jìn)算法都是在TIN算法基礎(chǔ)上進(jìn)行的,因此本文提出一種結(jié)合布料模擬濾波(cloth simulation filtering,CSF)[15]與TIN的點(diǎn)云濾波算法。CSF算法的主要優(yōu)點(diǎn)在于:① 可以應(yīng)用于陡坡、間斷區(qū)域和狹長地形;② 相比于其他濾波算法要設(shè)置復(fù)雜的濾波參數(shù),CSF算法參數(shù)設(shè)置較為簡便;③ 其他算法針對(duì)不同的地形要對(duì)算法參數(shù)進(jìn)行調(diào)試,找出最優(yōu)的地形值,而CSF算法的布料模擬率可以針對(duì)不同地形進(jìn)行選擇,自動(dòng)進(jìn)行處理[16-17]。但是CSF算法最后的閾值設(shè)置僅僅依賴于點(diǎn)云數(shù)據(jù)與格網(wǎng)數(shù)據(jù)之間的距離,因此存在一定的局限性。本文將CSF算法與TIN優(yōu)勢(shì)互補(bǔ),通過CSF獲取初始種子點(diǎn),然后利用初始種子點(diǎn)整體構(gòu)建格網(wǎng),通過判斷該點(diǎn)與三角網(wǎng)之間的角度和距離是否小于一定的閾值來獲取地面點(diǎn),將該地面點(diǎn)與之前的三角網(wǎng)再次整體構(gòu)網(wǎng),最終獲取完整地面點(diǎn)。

1 本文算法原理

本文算法主要包括點(diǎn)云粗差剔除、CSF粗濾波和改進(jìn)的TIN精濾波3個(gè)步驟。

1.1 點(diǎn)云的粗差提取

為減弱粗差對(duì)地面點(diǎn)提取的影響,本文利用CloudCompare軟件對(duì)機(jī)載LiDAR點(diǎn)云進(jìn)行粗差探測(cè)和剔除。

1.2 布料模擬粗濾波

傳統(tǒng)的CSF中,往往是通過判斷點(diǎn)云數(shù)據(jù)中的點(diǎn)與格網(wǎng)粒子之間的距離是否小于閾值,從而確定其是否為地面點(diǎn),但是這種方法通常會(huì)導(dǎo)致誤判或者漏判,并不能作為判斷地面點(diǎn)的準(zhǔn)確依據(jù)。本文首先采用CSF對(duì)粗差剔除后的點(diǎn)云進(jìn)行濾波,濾除建筑物、構(gòu)筑物等較大的地物,然后進(jìn)行種子點(diǎn)的選取,將其作為三角網(wǎng)濾波的初始種子點(diǎn),解決TIN初始種子點(diǎn)選取的難題,提高算法的精度和自適應(yīng)性。

CSF是基于表面的點(diǎn)云濾波算法,主要是對(duì)于物理過程的模擬。假設(shè)一塊足夠柔軟的布被放在地面上,這塊布因重力作用而下落,可以黏著在物體表面上,地形表面空洞的點(diǎn)最終會(huì)由于黏著在空洞周圍的非空洞點(diǎn)而拉回到地面,最終,布料表面的形狀就是數(shù)字地表模型(digital surface model,DSM)。但是,如果將地形上下翻轉(zhuǎn),并且該布料具有一定的硬度,那么最終布料的形狀就是DTM,以此為基礎(chǔ)將原始點(diǎn)劃分為地面和非地面部分。

本文中CSF的主要步驟如下:

(1) 將機(jī)載點(diǎn)云高程數(shù)據(jù)上下顛倒,如圖1所示。

圖1 點(diǎn)云數(shù)據(jù)倒置

(2) 設(shè)置布料網(wǎng)格分辨率,從而得到格網(wǎng)粒子,但是格網(wǎng)中所有點(diǎn)的高程坐標(biāo)在反轉(zhuǎn)點(diǎn)云的最高點(diǎn)之上。

(3) 把所有點(diǎn)云數(shù)據(jù)與格網(wǎng)數(shù)據(jù)投影到一個(gè)水平面,并為每個(gè)粒子找到其最鄰近點(diǎn)(corresponding point,CP),記錄其投影前的高程(intersection height value,IHV)。

(4) 對(duì)于所有可移動(dòng)的格網(wǎng)粒子,粒子的位置移動(dòng)取決于外部和內(nèi)部2個(gè)方面的驅(qū)動(dòng)因素,表達(dá)式為:

(1)

其中:m為粒子的質(zhì)量,其值通常設(shè)置為1;X為時(shí)刻t粒子的位置;Fext(X,t)為格網(wǎng)粒子受到的外力(重力或者遇到障礙物的碰撞力);Fint(X,t)為粒子在位置X和時(shí)刻t的內(nèi)力(相互連接的力)。內(nèi)部力和外部力隨著時(shí)間t變化而變化。

粒子受到外部驅(qū)動(dòng)因素影響而產(chǎn)生的位移表達(dá)式為:

(2)

其中:Δt為時(shí)間步長;G為常數(shù);X(t+Δt)為下一步點(diǎn)的位置;X(t)為當(dāng)前點(diǎn)位置;X(t-Δt)為上一步點(diǎn)的位置。若時(shí)間步長和初始位置都確定,則可以計(jì)算粒子的當(dāng)前位置。

將粒子當(dāng)前位置與該格網(wǎng)粒子對(duì)應(yīng)CP點(diǎn)的IHV進(jìn)行對(duì)比,若格網(wǎng)粒子高程等于或者低于IHV,則該粒子的高程被設(shè)置為IHV并且為不可移點(diǎn)。

(5) 對(duì)于格網(wǎng)粒子,計(jì)算每個(gè)粒子受到內(nèi)部驅(qū)動(dòng)因素所產(chǎn)生的位移,若2個(gè)相鄰點(diǎn)是可移動(dòng)點(diǎn),則其中一點(diǎn)向上移動(dòng),另一點(diǎn)向下移動(dòng);若2個(gè)相鄰點(diǎn)中的一點(diǎn)是不可移動(dòng)點(diǎn),另一點(diǎn)是可移動(dòng)點(diǎn),則可移動(dòng)點(diǎn)的移動(dòng)距離與布料的硬度相關(guān),布料硬度越大,移動(dòng)的距離越小,如圖2所示。

圖2 CSF內(nèi)力作用

不可移動(dòng)點(diǎn)與可移動(dòng)點(diǎn)之間的距離為d,當(dāng)布料硬度為1時(shí),可移動(dòng)點(diǎn)的距離上移d/2;當(dāng)布料硬度為2時(shí),可移動(dòng)點(diǎn)的距離上移3d/4;當(dāng)布料硬度為3時(shí),可移動(dòng)點(diǎn)的距離上移7d/8。

位移量公式如下:

(3)

其中:di為粒子的位移量;當(dāng)粒子為可移動(dòng)時(shí),b=1,不可移動(dòng)時(shí)b=0;P0為待定的移動(dòng)點(diǎn);Pi為P0的相鄰粒子;n為點(diǎn)在垂直方向上的單位向量[0 0 1]T。

(6) 重復(fù)進(jìn)行步驟(4)、步驟(5)計(jì)算,當(dāng)?shù)螖?shù)達(dá)到設(shè)置的最大迭代次數(shù)或者所有粒子的移動(dòng)距離都小于閾值時(shí),停止迭代。

(7) 辨別地面點(diǎn)與非地面點(diǎn)的方法為:若LiDAR點(diǎn)云與格網(wǎng)粒子之間的距離小于預(yù)先設(shè)置的閾值,則認(rèn)為其是地面點(diǎn);反之則為地物點(diǎn)。閾值可以選擇一個(gè)較小的值,從而為以后的TIN濾波提供一定數(shù)量的種子點(diǎn)。

1.3 改進(jìn)的TIN精濾波

傳統(tǒng)的TIN濾波算法中,初始種子點(diǎn)的選取主要依靠最大建筑物的尺寸及其范圍,通常會(huì)導(dǎo)致種子點(diǎn)的密度較小,初始種子點(diǎn)構(gòu)成的TIN不能充分表示地貌,因此,本文采用CSF獲取格網(wǎng)粒子作為構(gòu)建初始TIN的種子點(diǎn),解決傳統(tǒng)TIN點(diǎn)云濾波算法中初始種子點(diǎn)密度過小的問題,增加算法的魯棒性和自適應(yīng)性。

本文改進(jìn)的TIN精濾波主要步驟如下:

(1) 首先對(duì)CSF后的點(diǎn)云進(jìn)行分析,設(shè)置網(wǎng)格分辨率,對(duì)其進(jìn)行格網(wǎng)劃分,并且選取每個(gè)格網(wǎng)中的最低點(diǎn)為初始種子點(diǎn)。

(2) 利用點(diǎn)云整體四周的邊界點(diǎn)為輔助點(diǎn),邊界點(diǎn)的高程一般為最近的種子點(diǎn)的高程,然后用種子點(diǎn)與輔助點(diǎn)構(gòu)造初始的Delaunay三角網(wǎng)。

(3) 遍歷所有三角形和點(diǎn),確定每個(gè)點(diǎn)在哪個(gè)三角形的投影中,并且該三角形與該點(diǎn)構(gòu)成一個(gè)四面體,如圖3所示。

圖3 TIN濾波原理示意圖

其中點(diǎn)P為三角形上方的任意離散點(diǎn),點(diǎn)P到點(diǎn)A的連線與三角形平面構(gòu)成∠a,點(diǎn)P到點(diǎn)B的連線與三角形平面構(gòu)成∠b, 點(diǎn)P到點(diǎn)C的連線與三角形平面構(gòu)成∠c,點(diǎn)P到三角形平面的距離為d′,如果d′、∠a、∠b、∠c都小于某一設(shè)定的閾值,那么該點(diǎn)為地面點(diǎn),反之為地物點(diǎn),角度閾值一般設(shè)置為6°～50°,每次迭代過程中可以適當(dāng)減小閾值,最終閾值一般為2°～10°。

(4) 取濾波后的地面點(diǎn)和種子點(diǎn)重新構(gòu)建TIN。

(5) 重復(fù)步驟(3)、步驟(4),進(jìn)行迭代,直到?jīng)]有新的地面點(diǎn)加入或達(dá)到最大迭代次數(shù)。

2 實(shí)驗(yàn)結(jié)果和分析

為了驗(yàn)證該算法的有效性,選取國際攝影測(cè)量與遙感學(xué)會(huì)(International Society for Photogrammetry and Remote Sensing,ISPRS)網(wǎng)站(https://www.isprs.org/)的3組濾波測(cè)試樣本數(shù)據(jù)進(jìn)行實(shí)驗(yàn),原始點(diǎn)云數(shù)據(jù)如圖4所示。第1組為Samp11點(diǎn)云數(shù)據(jù)(圖4a),覆蓋的區(qū)域地形起伏較大,包含較多復(fù)雜建筑物,還有路燈、汽車等小型移動(dòng)地物,左上方有一個(gè)山坡,整個(gè)區(qū)域包括37 922個(gè)離散點(diǎn),區(qū)域高差為107.31 m;第2組為Samp31點(diǎn)云數(shù)據(jù)(圖4b),覆蓋區(qū)域?yàn)槌鞘袇^(qū)域,地勢(shì)較為平緩,包含大量建筑物,區(qū)域內(nèi)包含28 836個(gè)離散點(diǎn),區(qū)域高差為35.49 m;第3組為Samp54點(diǎn)云數(shù)據(jù)(圖4c),地勢(shì)較為平緩,包含大量植被及房屋,區(qū)域內(nèi)包含8 596個(gè)離散點(diǎn),區(qū)域高差為32.99 m。

圖4 3組原始點(diǎn)云數(shù)據(jù)

本文算法采用點(diǎn)云庫(Point Cloud Library,PCL)C++語言編程對(duì)上述3組點(diǎn)云數(shù)據(jù)進(jìn)行濾波實(shí)驗(yàn),并且將濾波后的點(diǎn)云數(shù)據(jù)與ISPRS官網(wǎng)生成的點(diǎn)云數(shù)據(jù)(參考數(shù)據(jù))進(jìn)行對(duì)比,如圖5～圖7所示。從圖5～圖7可以看出,本文算法對(duì)各種場(chǎng)景的點(diǎn)云數(shù)據(jù)有較好的濾波效果,可較好地濾除建筑物、構(gòu)筑物和樹木等。

圖5 Samp11濾波結(jié)果

圖6 Samp31濾波結(jié)果

圖7 Samp54濾波結(jié)果

根據(jù)ISPRS規(guī)定的Ⅰ類誤差、Ⅱ類誤差評(píng)價(jià)標(biāo)準(zhǔn)對(duì)實(shí)驗(yàn)結(jié)果進(jìn)行精度評(píng)定,為了便于對(duì)比分析,同時(shí)采用CSF算法、傳統(tǒng)TIN算法對(duì)3組實(shí)驗(yàn)數(shù)據(jù)進(jìn)行濾波實(shí)驗(yàn)并進(jìn)行精度評(píng)價(jià),Ⅰ類誤差、Ⅱ類誤差、總誤差、E1和E2結(jié)果見表1所列。表1中:Ⅰ類誤差為將原本是地面點(diǎn)的點(diǎn)劃分入地物點(diǎn)集;Ⅱ類誤差為將非地面點(diǎn)錯(cuò)誤地劃分入地面點(diǎn)集;總誤差為Ⅰ類誤差和Ⅱ類誤差占總點(diǎn)數(shù)的比例;E1為本文算法與CSF算法相比的誤差變化量;E2為本文算法與TIN算法相比的誤差變化量。

從Ⅰ類誤差和Ⅱ類誤差看,本文算法在Samp11中的Ⅰ類誤差較低,說明該算法對(duì)于坡度較大的地區(qū),地物點(diǎn)的精度較高,而在Samp31和Samp54中,本文算法Ⅰ類誤差與CSF算法相比偏高，但與TIN算法相比較小,說明其在地勢(shì)坡度較低區(qū)域的地物點(diǎn)精度較低,因此本文算法在平坦地勢(shì)上低矮地物的激光腳點(diǎn)存在欠濾波的現(xiàn)象,對(duì)于此類地物該算法的濾波效果有待提高。

與傳統(tǒng)TIN算法相比,本文Ⅰ類誤差降低幅度最大,在坡度較大的區(qū)域Samp11中降低50%以上,在建筑物密集區(qū)域Samp31、Samp54中降低幅度較小;與CSF算法相比,本文算法在Samp11中Ⅰ類誤差降低7.35%。

對(duì)Ⅱ類誤差分析可知,相比于CSF算法,本文算法在降低Ⅱ類誤差方面具有一定的優(yōu)勢(shì),在Samp31中降低1.45%;相比于TIN算法,本文算法Ⅱ類誤差值增大,濾波精度降低。

從總誤差看,本文算法在Samp11中較傳統(tǒng)TIN算法最大降低26.51%,說明該算法與傳統(tǒng)TIN算法相比具有更高的精度;與CSF算法相比,本文算法總誤差在Samp11中增加0.47%,在Samp31中增加1.68%,在Samp54中增加5.68%。

綜上所述,本文算法在坡度較大的區(qū)域可以顯著地降低Ⅰ類誤差,并將Ⅱ類誤差控制在一定的范圍內(nèi),驗(yàn)證了該算法具有較好的普適性與有效性,能夠獲得更接近真實(shí)地形的高精度數(shù)字高程模型(digital elevation model,DEM)。

表1 3種算法濾波結(jié)果誤差對(duì)比 單位:%

3 結(jié) 論

本文提出一種CSF 與TIN相結(jié)合的地面點(diǎn)濾波算法,該算法在坡度較大的區(qū)域中可以提高地面點(diǎn)提取精度。與傳統(tǒng)TIN濾波方法相比,本文算法在CSF的結(jié)果上適當(dāng)獲取種子點(diǎn),不僅可以避免地形因素的干擾,而且在種子點(diǎn)選取上可以不受建筑物尺寸影響,解決了種子點(diǎn)選取敏感問題和閾值選取問題。與CSF算法相比,本文算法可以在較大坡度的區(qū)域提高地面點(diǎn)提取精度,但是在坡度較低的區(qū)域內(nèi),提取點(diǎn)云點(diǎn)的精度相對(duì)降低,但可以保持在一定的范圍內(nèi)。本文算法需針對(duì)坡度較低的區(qū)域進(jìn)行研究,進(jìn)而保證其在坡度較低區(qū)域中地面點(diǎn)提取精度。