基于多尺度網格的點云自適應坡度濾波算法

趙明君, 劉 超, 高 翔, 李慧慧

(1.安徽理工大學 空間信息與測繪工程學院,安徽 淮南 232001; 2.安徽理工大學 礦山采動災害空天地協同監測與預警安徽普通高校重點實驗室,安徽 淮南 232001)

0 引 言

三維激光掃描技術是一種新型的高時空分辨率、高精度的地理信息數據觀測手段[1],所獲取的點云數據包括測量區域內所有的地表及地物的三維信息,因此,點云數據后處理算法研究中的關鍵問題,就是如何從大量的散亂點云數據中快速準確地將地表點云數據與其他地物點云數據區分[2]。點云濾波就是指區分點云中地面點和非地面點的過程[3]。

在點云數據濾波方面,常用的濾波算法有4種,分別為基于數學形態學、基于不規則三角網(triangulated irregular network,TIN)、基于曲面擬合、基于坡度[4-7]。文獻[8]提出了基于數學形態學濾波算法,該算法的關鍵步驟是如何選擇合適的濾波窗口;文獻[9]提出了一種漸進式形態學濾波方法,雖然取得了較好的點云濾波效果,但是在采用大窗口對點云數據濾波時,仍容易將地形起伏劇烈處的點誤判為非地面點,降低了分類精度;文獻[10]給出了一種基于TIN的點云濾波算法,該算法是基于二維層次的搜索方法,其計算量和算法的復雜度相對較大;基于曲面擬合算法的核心是通過足夠的點擬合出較準確的地形曲面,再通過高差值判斷緩沖區內的點是否為地面點[11];文獻[12]提出了一種基于坡度的濾波算法,并將其應用到點云濾波過程中,該算法的基本思路是利用地面點與其鄰近的非地面點之間的坡度值大于設置閾值來區分非地面點與地面點。

上述4種點云濾波方法中,基于坡度濾波算法的濾波原理較簡單,可以適應多種類型的點云數據,但在實際復雜多樣的地形中,該算法對整個點云區域按單一的坡度閾值進行濾波是不完全合理的。文獻[13]針對該缺點進行了改進,提出了一個基于核函數的濾波方法,增強了坡度濾波算法在陡峭地形區域中的濾波精度;文獻[14]中利用數字地面模型(digital terrain model,DTM)進行分層迭代更新得出坡度閾值,得到了坡度閾值隨著地形坡度變化而變化的濾波效果;文獻[15]提出在陡坎、斜坡等地形變化較大處通過設置坡度增量來區分非地面點與地面點,以減小坡度濾波時的分類誤差;文獻[16]通過增加不斷變化的窗口,利用DTM不斷更新來實現坡度閾值的適應性設置,可以顧及點云區域的整體地形起伏和局部細節。這些基于傳統濾波算法的改進方法都是基于整體地形變化而更新坡度閾值設置的,沒有考慮到局部細節的地形結構;同時多次的迭代計算也增加了算法的復雜度,降低了算法的效率。因此,上述改進方法并未完全實現坡度閾值的自適應設置。

坡度濾波算法需要解決以下2個難點:① 如何根據地表地形變化得出最優坡度濾波閾值;② 在地形變化劇烈的邊緣帶過度濾波造成的地表點云數據損失問題。為此,本文提出了一種基于多尺度網格的點云自適應坡度濾波算法(簡稱“本文算法”),通過統計分析局部地形坡度值得出最優坡度閾值,實現了真實地形的坡度閾值自適應設置,從而分離點云數據中的地面點與非地面點。

1 點云自適應坡度濾波算法

本文算法主要內容包括:① 點云數據組織,即進行孤立點去除與網格構建;② 坡度閾值計算與自適應設置,即進行網格屬性劃分與種子點選取、點云數據的自適應坡度濾波、多尺度網格的點云重復濾波;③ 地面點云數據去噪與數據輸出。該算法對復雜地形自適應設置坡度閾值,并通過不同尺度網格重復濾波,盡可能剔除非地面點,保留真實地形。點云自適應坡度濾波算法流程圖如圖1所示。

圖1 點云自適應坡度濾波算法流程圖

1.1 多尺度網格構建與孤立點去除

(1) 多尺度網格構建。原始點云數據基本上是散亂無序的,點間無拓撲關系。為了克服算法運算效率低和信息損失的問題,在光探測和測距(Light Detection and Ranging,LiDAR)(簡稱“激光雷達”)數據處理中引入虛擬網格概念。多尺度網格的示意圖如圖2所示。

圖2a中,圓點表示點云數據,長方體表示相應尺度的網格;圖2b所示為在虛擬網格內將三維坐標點投影到X-Y平面,方格的顏色由深到淺表示虛擬網格尺度由小到大;圖2c所示為基于點的高程對網格分層時點云的存在形式;圖2d所示為包含地面點與植被點時,點云數據在網格中的空間位置。

由圖2b可知,可以將點云數據按照其X-Y坐標與網格建立起索引機制,每個點所在的網格或每個網格包含的點均可以快速查詢。

圖2 多尺度網格示意圖

點與網格間的索引關系計算公式為:

(1)

其中:(X,Y)為網格號;(x,y)為點云的平面坐標;(xmin,ymin)為整個數據集的最小平面坐標;m為網格單元的尺度;INT表示對計算結果向下取整。

(2) 孤立點去除。從圖2c可以看出,通過在網格內進行分層高程統計,可以得出在不同高程區間內的點云個數,此時若某一層中點云個數不足平坦區域內點云密度的10%,則認為該層的點云為孤立點,并進行刪除。此步驟主要是減少孤立點,為后續運算得到較純凈的點云數據。

1.2 坡度閾值自適應設置與點云濾波

1.2.1 網格屬性劃分和種子點選取

(1) 網格屬性劃分。從圖2d可以看出,該網格包含了地面點云數據和植被點云數據。當網格平面尺度大小一定時,給定網格高度,可以定義空間占比F為:

(2)

其中:hi為網格內第i點高程;hmin′為網格內次低點高程;H為網格高度;n為網格內點云個數。F反映了點云數據在空間存在形式的差異,根據F值可以定義網格的不同語義屬性,同時根據網格屬性,可以在后續計算中給出不同網格的地形因子,為坡度閾值的計算分配參數。選取1個實驗區,網格平面尺度為5 m,H=25 m,計算不同屬性網格的F值,結果見表1所列。

表1 實驗區不同語義屬性的空間占比

(2) 初始地面種子點選取。初始地面種子點的準確選取是基于坡度的濾波算法有效運行的關鍵因素。原始點云數據經過孤立點去除和點云格網構建后,通常把每個網格內的高程最低點作為地面種子點,但該方法選取的種子點不是最優的地面種子點。為提高所選取的地面種子點的可靠性,對種子點選取方法做如下改進:① 初始地面種子點必須在最大尺度的網格內選取,為保證每個網格內包含地面種子點,網格的最大尺度應盡量大于點云區域內建筑物或樹冠的寬度;② 為避免所選取的最低點是局部微地形中位于地面以下的噪聲點,對每個網格先按照高程由高到低排序,然后選取高程值從最小值開始位置排列在點數5%的激光點作為地面種子點;③ 如果少量網格包含的點全是非地面點,那么無法選取地面點,此時計算該網格內待擬定的種子點與周圍8個鄰格內種子點的高差值Δhi(i=1,2,3,…,8),若Δhi均大于網格尺度值的2倍,則認為該種子點為非地面點,予以去除,并刪除該網格內的點云數據。上述的3個約束條件的目的是最大概率保證選取的地面種子點為真實地面點。

1.2.2 坡度閾值求解與自適應設置

根據F值劃分出網格屬性并選出單一網格內的地面種子點后,對于不同的網格屬性需要給出相應的地形計算因子σ,σ為網格內點云與種子點間高差的倍數。網格的語義屬性僅能粗略地描述地形與地物的空間關系。對于地形復雜的區域,網格內局部地形的變化仍然存在。為更加真實反演網格內的地形變化,在網格內計算出激光點與種子點的平面距離,以種子點為圓心,按平面距離從小到大、等點數間隔逐步分割的原則,將每個網格分割成若干個輻射圈。在每個輻射圈選取高差值為hpt時所對應的激光點為參考點pt,多個參考點pt即為參考點集Pt。

hpt=max{ht≤σumin,umin=minHt}，hpt∈{Ht}

(3)

其中:Ht為輻射圈內點云數據中第t點與該輻射圈內高程最低激光點間的高差;σ為該網格對應的地形計算因子。計算參考點集與種子點間的坡度值si,則其中的最大值為坡度閾值Si。在不同的網格內得出不同的坡度閾值,以解決坡度閾值隨地形變化的自適應設置問題。實驗區不同語義屬性下的σ見表2所列,該σ是通過地面點樣本實驗得出的。

表2 實驗區不同語義屬性的地形計算因子

1.2.3 點云多尺度濾波

通過σ得出單個網格的Si后,計算網格內所有點與種子點的坡度值,并與Si進行比較。當坡度值大于Si時,將該點分為非地面點,予以去除。進行濾波操作時,算法按照網格平面尺度由大到小重復濾波。在較大尺度網格進行濾波時,可以去除較大的植被與建筑物點云；在小尺度網格進行濾波時, 可以濾除與地面點接近的地物點,如低矮植被、草叢、低矮地物等。按固定尺度間隔由大到小重復濾波使得較多的地面點數據得以保留。

1.3 地面點去噪與數據輸出

通過自適應坡度濾波后的點云中仍然可能存在少量噪聲點。這是由于少部分區域的網格屬性或σ值很難完全符合真實地形的表達。此類的噪聲點多為孤立點,因此設立較小尺度的網格,求取待定點的網格序號及其鄰域網格,當鄰域網格內不含點云的網格個數大于4時,則認為該點為噪聲點,去除該點。最后去噪后的點云數據即為提取出的地面點。

1.4 濾波精度評價參數

為了量化表達基于多尺度網格的點云自適應濾波算法在不同場景中的濾波效果,利用點云數據分類中常用的混淆矩陣來定量評價算法的濾波精度。令T表示正確分類,F表示誤分類,P表示地面點,N表示非地面點,則分類結果有以下4種:① TP,即濾波結果中分類正確的地面點數(TP);② FP,即非地面點被誤分類為地面點的個數(FP);③ TN,即分類正確的非地面點數(TN);④ FN,即地面點被誤分類為非地面點地的個數(FN)。濾波精度的評價參數有Ⅰ類誤差e1、Ⅱ類誤差e2、總體誤差e及Kappa系數。e1、e2、e的計算公式為：

e1=FN/(TP+FN),e2=FP/(FP+TN),

e=(FN+FP)/N總,N總=TP+FP+TN+FN。

Ⅰ類誤差反映了提取地面點的完整度,Ⅱ類誤差反映了提取地面點的純度,總體誤差反映了算法的可行性。

2 算法實驗與結果分析

本文在PCL點云庫的基礎上采用Matlab進行計算和點云數據濾波成果展示。 在安徽理工大學校內選取多處復雜的地形場景,利用地基站掃描儀和無人機攝影測量獲取實驗數據,并從定性和定量計算2個角度進行點云濾波實驗結果分析。

2.1 單一場景下的點云濾波實驗

實驗區域為安徽理工大學校園內一處緩坡,所得數據為范圍較小的邊坡點云數據。該組數據采用海達云HS650脈沖式三維激光掃描儀獲得,數據為三維坐標值,無顏色及紋理信息。實驗區域的點云平均密度為11 200 點/m2,區域的長度約為18.0 m,寬度約為7.5 m,最大的高程差為5.0 m。實驗區主要為緩坡地形,無較大的起伏,最大坡度約為30°,包括的地物為樟樹、低矮的闊葉樹及少量雜草。 單一場景實驗區域的點云數據及濾波結果如圖3所示。

圖3 單一場景實驗點云數據及濾波結果

從圖3b可以看出，實驗區內的植被點被有效分離,包括與地面較接近的低矮植被也得到濾除,同時整個地面的點云數據得到保留,經過去噪處理后,實現了地面地形的真實還原。 由于采用的數據是由地基站掃描儀得到的,在圖3c中出現的裂縫是在掃描過程中被植被遮擋而形成的空洞。從圖3d可以看出，在濾除非地面點的同時,盡可能地保證了分離出的非地面點是真實地表中的非地面點,有效減少了地面點的誤分類。

2.2 復雜度場景下的點云濾波實驗

在安徽理工大學校園內選取3塊實驗區。實驗區的選取標準是盡量表達復雜的自然場景, 包括平地、緩坡、陡坡、河流、道路、建筑群、大小不同的植被、灌木從等,以驗證本文算法在復雜地形場景中的濾波效果和算法的實用性。

該3塊實驗區域采用大疆精靈4RTK無人機飛行器拍攝影像數據,再對影像數據進行三維實景模型構建,最后輸出數據格式為LAS的點云數據。點云數據包括點位三維坐標值、點位顏色屬性值及點位強度信息等。

2.2.1 多地物的區域濾波實驗

該區域為區塊1，平均點云密度為260點/m2,區域的長度約為100.0 m,寬度約為56.0 m,最大的高程差為14.0 m。整個區域包含復雜地形和大量植被,是典型的濾波難度大的實驗區域。該區域的南邊為平地和斜坡,斜坡中部高程值較小,兩邊地形坡度較大;區域中間部分為寬度約為12 m的河溝,北邊為坡度變化小的緩坡;區域內地物種類繁多,有涼亭、假山石、護欄、較高的梧桐樹、中等高度的樟樹、地表灌木叢、河中的蘆葦叢等,同時測區內含有2條完整道路線。區塊1的點云數據及濾波結果如圖4所示。

圖4 區塊1的點云數據及濾波結果

算法的初始網格寬度為8.0 m,最低網格寬度為0.5 m。從圖4b可以看出，與原始數據對比,地表上的地物基本上得到濾除，包括涼亭、梧桐樹、樟樹及地表的灌木叢；但是在河中的蘆葦叢處仍然存在較多的點云數據,通過實際考察可知，蘆葦叢高度與河岸邊高度基本一致,在濾波閾值計算時不能很好地分離數據,造成了點云的殘留。從圖4b還可以看出河溝的邊岸線和部分邊坡都得到保留,整個地面和坡面數據保留較為完整,無噪聲點。圖4b中出現的點云空洞是由于數據是由影像數據生成的點云數據, 存在較多地面點被植被遮擋的情況,但是濾波結果對于地形的表達較完整。

2.2.2 植被覆蓋度高的區域濾波實驗

該區域為區塊2，主要為坡度變化明顯的斜坡和大量植被,植被有樟樹、針葉松、闊葉樹、地表灌木、地表草叢等。該區域地形兩邊為緩坡,中間為凸出的坡地；平均點云密度為960點/m2,區域的長度約為50.0 m,寬度約為40.0 m,最大的高程差為7.8 m。區塊2的點云數據及濾波結果如圖5所示。

因為區塊內的非地面點僅包含植被點,所以以最大樹冠直徑5.0 m作為初始網格寬度,最低網格寬度為0.3 m,以達到非地面點徹底濾除的效果。 對比圖5c與圖5d可知,該區域內的非地面點濾除效果良好,并且通過局部地形自適應坡度閾值的設置可以很好地在顧及地形變化的同時剔除非地面點。不同尺度的網格設置可以去除不同大小的植被。點云中的空洞是濾波算法將非地面點剔除和數據源2個方面造成的。

2.2.3 附含建筑物群的區域濾波實驗

該區域為區塊3，平均點云密度為570點/m2,區域的長度約為80.0 m,寬度約為40.0 m,最大的高程差為22.0 m。該區域地形較簡單,主要是緩坡和平坦地面;主要的地物類型是建筑物和地面植被。 區塊3的點云數據及濾波結果如圖6所示。

從圖6b可以看出，點云數據經過濾波后可以得到很完整的地面點云數據,區域中的建筑物和植被均已濾除。在實驗中，初始網格寬度設置時需要考慮網格內包含的真實地面點數據,因此初始寬度設置為10.0 m,最低寬度為0.5 m。實驗結果表明,建筑物點云很難一次徹底濾除, 對于這種混合地形的點云數據,需要多次計算才能得到真實地面點云。

圖6 區塊3點云數據及濾波結果

2.3 復雜場景下的點云濾波精度評定

基于TIN的點云濾波、基于二次曲面擬合濾波、傳統坡度濾波是常用的3種濾波算法，將本文算法與上述3種濾波算法進行對比。3個區塊4種濾波算法的分類結果、濾波誤差計算結果和算法耗時對比見表3所列。3個區塊的樣本點是利用Terrasoild軟件和手工交互的方式進行分類后的地面點和非地面點。

表3 3個區塊4種濾波算法的濾波誤差計算結果和算法耗時對比

點云濾波算法一般要求在控制Ⅱ類誤差的基礎上盡量降低Ⅰ類誤差,即在保證分離出的地面點最大程度上是地面點的同時,盡可能分離出較多的地面點。Kappa系數的大小反映了分離出的地面點位置的準確性程度。總體精度和Kappa系數保證了濾波得到的點云數據的數量和質量。

通過計算結果可以分析得出:

(1) 本文算法與傳統坡度濾波算法相比,在總體精度和Kappa系數上都有所提高,表明本文算法分離出的地面點更多,并且與實際地表地形更加符合;但是在計算效率上小于傳統的坡度濾波算法,這是由于多尺度網格的設立使得本文算法進行了多次運算,造成了本文算法計算效率降低。同其他2種算法相比,本文算法計算效率較高。

(2) 在同一區塊，本文算法總體精度較高,表明本文算法與其他3種濾波算法相比同樣具有很好的應用價值,甚至優于其他3種算法;本文算法Kappa系數值較大或者接近于其他3種算法,尤其在區塊1的復雜地形中,本文算法的Kappa系數值最高,表明其在地形模擬和地面點提取方面具有很好的效果。

(3) 區塊1內4種算法的分類精度和Kappa系數都低于區塊2和區塊3,這主要是由于區塊1的實際地形復雜,地面點和非地面點的分離界限不明確,濾波算法很難從無序點云中提取出地面點;但是Kappa系數值在0.7附近,說明可以獲得大體的地形數據,但是局部細節仍有損失。區塊2內的地形包括平地、緩坡及斜坡,但是整體地形變化較均勻,含有的低矮植被掩蓋了真實地表信息,本文算法的總體精度和計算效率雖然不是最優,但是仍具有較好的濾波效果;區塊3內4種算法分類精度和Kappa系數值優于區塊1和區塊2,但實際實驗中該區塊由于含有高大建筑物,需要經過多次濾波后才能得出最佳分類效果。

3 結 論

為實現海量LiDAR點云數據中地面點與非地面點的有效分離,本文在傳統坡度濾波算法的基礎上,提出基于多尺度網格的點云自適應坡度濾波算法。將本文算法在3個復雜場景下進行應用驗證，并與3種常用濾波算法的分類精度進行對比,得到以下結論:

(1) 本文提出在構建的虛擬網格內計算點云空間占比,并以此劃分網格語義屬性,通過地形計算因子計算不同網格的坡度閾值,達到坡度濾波閾值隨地形變化自適應設置的目的。

(2) 多級尺度網格構建和坡度閾值的自適應設置在一定程度上解決了坡度濾波算法在邊緣和陡坡區域的過濾波問題,有效降低了坡度濾波算法的Ⅱ類誤差。

(3) 本文算法總體精度大于85%,Kappa系數值大于0.7,表明該算法在大多數地形中有較好的適用性和分類的準確性。

(4) 本文算法雖然可以自適應設置坡度閾值,但是前期設置的參數需要選取樣本值進行實驗來獲取,增加了計算難度,該問題需進一步研究;對于混合地形區域,本文算法仍有局限性,需要多次實驗才能取得良好的濾波效果。