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

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

0 引 言

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

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

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

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

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

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

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

1.1 多尺度網(wǎng)格構建與孤立點去除

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

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

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

圖2 多尺度網(wǎng)格示意圖

點與網(wǎng)格間的索引關系計算公式為:

(1)

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

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

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

1.2.1 網(wǎng)格屬性劃分和種子點選取

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

(2)

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

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

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

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

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

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

(3)

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

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

1.2.3 點云多尺度濾波

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

1.3 地面點去噪與數(shù)據(jù)輸出

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

1.4 濾波精度評價參數(shù)

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

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

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

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

2 算法實驗與結果分析

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

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

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

圖3 單一場景實驗點云數(shù)據(jù)及濾波結果

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

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

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

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

2.2.1 多地物的區(qū)域濾波實驗

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

圖4 區(qū)塊1的點云數(shù)據(jù)及濾波結果

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

2.2.2 植被覆蓋度高的區(qū)域濾波實驗

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

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

2.2.3 附含建筑物群的區(qū)域濾波實驗

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

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

圖6 區(qū)塊3點云數(shù)據(jù)及濾波結果

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

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

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

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

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

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

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

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

3 結 論

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

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

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

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

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