數(shù)學形態(tài)學與布料模擬濾波算法適用性分析

王 凱，魯鐵定*，魯春陽，劉 盈

(1.東華理工大學測繪工程學院，330013，南昌；2.河南城建學院測繪與城市空間信息學院，467036，河南，平頂山；3.井岡山大學電子與信息工程學院，343009，江西，吉安)

0 引言

機載激光探測和測距(LiDAR)作為一種高效、快速的大面積采集三維(3D)點云的遙感技術，已經在各個領域得到了廣泛的應用，例如生成數(shù)字地形模型(DTM)[1]、森林生態(tài)系統(tǒng)調查[2]和三維建筑建模[3]。大多數(shù)點云數(shù)據應用中，區(qū)分為地面點和非地面點是一個必不可少的步驟，該過程也稱為點云濾波。在過去的幾十年中，眾多學者根據不同的理論基礎，提出了許多自動過濾地面點的點云濾波算法，這些算法主要可分為6類[4]。

基于數(shù)學形態(tài)學的點云濾波算法具有原理簡單、實現(xiàn)效率高的特點，被眾多學者研究，并進行相關改進。形態(tài)學濾波(Morphological filters)首先應用于數(shù)字圖像處理，其能夠去除圖像中的細小部分，從而使圖像更加平滑。Lindenberger[4]等首先將數(shù)學形態(tài)學方法引入到機載LiDAR點云濾波中。Zhang[5]等提出了一種經典的漸進形態(tài)學濾波方法。在該方法中，濾波窗口按一定的規(guī)則從小到大，不同的濾波窗口對應不同的高程閾值。此后，許多學者大都基于此方法進行形態(tài)學濾波算法改進。改進主要從以下3個方面入手[6]：1)進行格網誤差的改進；2)對地形坡度、最大濾波窗口、高差閾值等參數(shù)設置的改進；3)地形坡度常量假設的改進。盡管已經產生了許多不同的改進算法，但是基于數(shù)學形態(tài)學的濾波算法對濾波窗口大小的選取敏感，以及很容易對陡峭地形中的點進行錯誤分類，如何提高基于形態(tài)學的濾波算法在地形起伏較大的區(qū)域中的魯棒性仍是眾多學者研究的重點與難點。

Zhang[7]等2016年提出的布料模擬濾波算法(CSF)，因其用戶定義的參數(shù)少且易于設置，開始受到廣泛關注。張凡[8]等將布料模擬算法在LiDAR點云中進行生成DEM實驗，結果表明該算法生成的DEM能夠很好地表現(xiàn)連續(xù)起伏區(qū)域的地形和細節(jié)特征，與實際地表地形能夠很好地擬合。李雅盟[9]等提出一種顧及地形特征的布料模擬濾波改進方法，采用基于坡度的動態(tài)格網分割方法，對粗差剔除后的點云建立格網索引；然后利用每個格網的鄰域格網中的最低點建立曲面方程來擬合高程值，最后計算真實高程與擬合高程差值從而實現(xiàn)高程歸一化；最后使用布料模擬算法模擬布料下降過程得到地形布料的最終形態(tài)，進而通過限定閾值實現(xiàn)地面點提取。結果表明，該改進算法在保證大范圍復雜場景區(qū)域濾波正確率的基礎上，對不同地形具有較強的自適應性，且提高了濾波計算效率。Cai[10]等提出了一種將CSF和漸進不規(guī)則三角網加密(PTDF)相結合的濾波算法。在所提出的算法中，利用CS得到了一種高質量的初始臨時數(shù)字地形模型(DTM)，并根據統(tǒng)計分析，從初始臨時DTM中估計PTD的參數(shù)閾值。最后，利用具有自適應參數(shù)閾值的PTD對初始臨時DTM進行了細化。這些細化實現(xiàn)了濾波精度的提高以及相關參數(shù)的優(yōu)化。該算法的濾波精度優(yōu)于直接使用PTDF算法。此外，與公開的改進PTDF算法相比，該算法不僅精度高，且更實用。針對布料模擬的相關特點，在點云分類方面，佟國峰[11]等針對室外大場景提出了一種基于聚類分割的點云分類算法，該方法通過濾除法向量差大的點改善地面濾波效果，為后面的點云分割提供了更好的非地面點數(shù)據。在點云精簡算法方面，李繞波[12]等提出一種基于布料模擬濾波、方法庫和曲率分級等綜合算法的點云精簡優(yōu)化方法。由于布料模擬濾波算法(CSF)提出時間較晚，對布料模擬的其他相關應用，眾多學者還在不斷開發(fā)研究中。

本文通過對提出時間早、做過大量相關改進的經典的數(shù)學形態(tài)學濾波與最近提出的相關應用和改進都做的較少的布料模擬濾波算法進行濾波實驗并進行精度評定。比較、分析和總結2種濾波算法的優(yōu)缺點，對其適應性進行分析，希望能夠對以后相關算法的研究有參考借鑒意義。

1 點云濾波原理

1.1 數(shù)學形態(tài)學濾波

數(shù)學形態(tài)學濾波算法原理：點云在進行形態(tài)學運算后非地面點的高程會發(fā)生很大的變化。將高度差大于閾值的點分類為非地面點并進行濾除。濾波窗口的大小對濾波結果的影響很大。

膨脹、腐蝕，閉運算和開運算是構成數(shù)學形態(tài)學運算的4個公式。f(x,y)為輸入圖像，b(x,y)為“結構元素”，定義如下。

1.1.1 膨脹 結構元素b對f進行灰度膨脹記為(f⊕b)，如式(1)：

(f⊕b)(s,t)=max{f(s-t,t-y)+b(x,y)|(s-x),(t-y)∈Df;(x,y)∈Db}

(1)

其中Df和Db分別是f和b的定義域。膨脹運算是在由結構元素確定的領域中，選取(f⊕b)的最大值。

1.1.2 腐蝕 結構元素b對f進行灰度腐蝕記為(fΘb)，如式(2)：

(fΘb)(s,t)=min{f(s+t,t+y)+b(x,y)|(s+x),(t+y)∈Df;(x,y)∈Db}

(2)

其中Df和Db分別是f和b的定義域。腐蝕是在由結構元素確定的領域中選取(fΘb)的最小值。

1.1.3 閉運算和開運算 由膨脹和腐蝕進行組合。

閉運算，用b閉運算f，記為f·b，定義如式(3)：

f·b=(f。b)⊕b

(3)

將結構元素內的所有數(shù)據填充為區(qū)域最大值，然后再填充為該區(qū)域最小值。

開運算，用b開運算f，記為f。b，定義如式(4)：

f。b=(f。b)⊕b

(4)

將結構元素內的所有數(shù)據填充為區(qū)域最小值，然后再填充為該區(qū)域最大值。

在漸進形態(tài)學濾波時，濾波結果的好壞對濾波窗口的大小選取十分敏感。對于濾波窗口的增長方式，本文選擇方法是通過線性方式增加窗口大小，如式(5)：

wk=2kb+1

(5)

其中k=0,1,2,3,…,m；b是初始濾波窗口，最大濾波窗為2mb+1；采用2kb+1作為窗口大小可保證濾波窗口繞中心點對稱，簡化了開運算。線性增加窗口大小的優(yōu)點是可以有效地保留漸進變化的地形特征。但是，對于具有大型非地面物體的區(qū)域，則需消耗大量的計算時間。

高程閾值可基于研究區(qū)域的地形坡度決定，對于最大高程差與地形dhmax(t)，k，窗口大小wk和地形坡度s的關系如公式(6)，得出坡度常數(shù)。

(6)

因此，高程閾值dhT,k可由公式(7)得：

(7)

其中：dh0是原始高程閾值；s是坡度；c是單元格網大小；dhmax是最大高程閾值。

以指數(shù)方式增長的濾波窗口優(yōu)點在于增長速度快、效率高，但缺點在于濾波精度與線性增長方式相比效果不是很好。因此，本文采取濾波窗口線性增長方式。窗口大小wk和高程閾值dhT,k上增長方式如式(6)、式(7)所示，使用C++語言與點云庫PCL相結合，并將所得的點云數(shù)據在VS2013平臺上進行數(shù)學形態(tài)學濾波實驗。

1.2 布料模擬濾波

CSF算法原理：首先使用布料模擬(CS)估算初始地形，然后基于實際地形應在某個局部區(qū)域內接近水平面的假設，僅通過使用高程信息將未過濾點的地面點添加到初始地形。布料模擬濾波算法(CSF)，首先將LiDAR點云翻轉倒置，即將點云坐標中的z變?yōu)?z，此后布料粒子因受重力影響而下降，通過分析布料節(jié)點與相對應的LiDAR點云之間的作用關系，從而確定布料粒子最終所處的位置，以此來確定布料最后所形成的形狀，從而達到點云濾波的目的，布料模擬過程如圖1所示[6]。

圖1 布料模擬過程圖

由于限制了布料粒子的移動方向，因此2個具有不同高度的粒子將試圖移動到同一高度。若2個相關聯(lián)的粒子都為可移動的，則沿相反方向將它們移動相同的位移。如果其中一個不可移動，則另一個將被移動。否則，若這2個粒子有同一的高度值，則不移動。因此，通過式(8)計算每個粒子的位移(矢量)：

(8)

圖2 剛性參數(shù)化

CSF的主要過程如下：先將布料粒子和LiDAR點云投影到同一水平面，其次在此2D平面中為每個布料粒子找到與其相臨近的LiDAR點(稱為對應點，CP)。定義交叉點高度值(IHV)以記錄CP的高度值(投影前，該值表示粒子可以到達的最低位置)。每次迭代過程中，將該粒子的當前高度值(CHV)與IHV進行比較；若CHV≦ IHV，則粒子移回IHV的位置，并設為不可移動。

CSF主要由4個用戶定義的參數(shù)組成：網格分辨率(GR)，表示相鄰粒子之間的水平距離；時間步長(DT)，表示每次迭代過程中粒子受重力影響而產生的位移；剛性度(RI)，其控制布料的剛度；距離閾值(hcc)，該距離閾值根據與布料網格的距離決定將LiDAR點的最終分類為BE(地面點)和OBJ(地物點)。

CSF算法的實現(xiàn)過程如下。

1)使用第三方軟件(如Cloudcompare)自動或手動去除躁。

2)使原始LiDAR點云倒置。

3)布料格網的生成，根據網格分辨率(GR)確定粒子數(shù),布料的起始位置一般設置在最高點之上。

4)將所有LiDAR點云和網格粒子投影到一個水平面，并找到該平面中每個網格粒子的對應點(CP)，并記錄IHV。

5)對每個布料粒子，若該粒子為可移動，計算其受重力影響產生的位移，并將該粒子的高度與IHV進行比較。若粒子的高度≦IHV，則將該粒子置于IHV的高度，并設為“不可移動”粒子。

6)計算每個布料粒子受內力影響產生的位移距離。

7)重復5)～6)至所有布料粒子的最大變化高度足夠小或超設定的最大迭代次數(shù)時，則布料模擬停止。

8)計算布料粒子與其相對應點(CP)之間的距離。

9)區(qū)分地面點和地物點。若每個LiDAR點到模擬粒子的距離小于距離閾值(hcc)，分類為BE(地面點)，否則分類為OBJ(地物點)。

除了這些用戶定義的參數(shù)外，根據CSF算法原理，當將原始LiDAR點云倒置時，地面上方的物體將出現(xiàn)在地面測量點以下。因此，根據實驗樣本地形特征的不同，在視覺上將實驗樣本分為不同的組別，以此依據來設定濾波參數(shù)。例如地形非常平坦且沒有陡峭或有坡度的梯田，則將RI設置為相對較大的值(RI= 3)，且不需要坡度處理(ST=false)；如存在陡峭的斜坡(例如河岸、溝渠和梯田)，則需要中等軟布(RI=2)且需要進行坡度處理(ST =true)；當處理非常陡峭的斜坡時，需要非常柔軟的布料(RI=1)且需要進行坡度處理(ST=true)。

2 實驗結果與分析

由于Zhang[5]等提出的漸進形態(tài)學濾波方法最為經典，較具有代表性。此后的相關數(shù)學形態(tài)學濾波方法的改進大都是基于這種經典算法的改進，本文選用Zhang[5]等提出的漸進形態(tài)學進行點云濾波實驗，與Zhang[7]等2016年提出的布料模擬的濾波算法進行實驗，最后分別對2種濾波算法進行精度評定并進行適應性分析。實驗數(shù)據采用ISPRS(國際攝影測量和遙感學會)數(shù)據，對2種算法的性能進行了測試，該實驗數(shù)據由Optech ALTM掃描儀收集，包括由不同地形環(huán)境組成的8個站點(命名為站點1～8)：4個城市站點和4個農村站點。從站點1～7中選擇了15個具有不同地形特征和土地覆蓋類型的參考樣本，進行濾波精度評定。由于缺乏參考數(shù)據，因此站點8被排除在外。對于每個實驗樣本，參考數(shù)據是用先驗知識和可用的航空圖像手工識別生成的[13]，15組樣本數(shù)據地形特征如表1所示。

表1 15組樣本的地形特征

本文實驗數(shù)據使用C++語言與點云庫PCL相結合，并將所得的點云數(shù)據在VS2013、點云處理軟件Cloudcompare進行數(shù)據處理。為盡量避免不同噪聲點對濾波結果的影響，在進行濾波實驗前，15組樣本數(shù)據都通過手工進行了去噪處理。此外，為定量分析2種算法的濾波性能，根據ISPRS規(guī)定的Ⅰ類誤差(將地面點錯分為地物點的誤差)、Ⅱ類誤差(將地物點錯分為地面點的誤差)和總誤差對實驗結果進行定量精度評定分析。15組樣本數(shù)據濾波實驗結果Ⅰ類誤差、Ⅱ類誤差和總誤差如表2所示。

表2 2種算法濾波結果

從表2中可以看出，2種算法在地形平坦區(qū)域(S21、S31、S42、S51、S54)都取得了良好的濾波效果，Ⅰ類誤差、Ⅱ類誤差和總誤差都處于較低水平，較具有代表性的有S21、S31，如圖3、圖4所示。

(a)S21 DSM (b)參考DEM

(a)S21 DSM (b)參考DEM

2種算法在城市等綜合復雜地形區(qū)域(S11、S22、S23、S24、S41)濾波效果一般，由于數(shù)學形態(tài)學濾波算法需要一定的先驗知識去預估最大建筑物尺寸，以設定最大濾波窗口大小，而CSF的局限性在于已將粒子運動的物理運動過程改為2個離散的步驟，當在城市等綜合地形復雜區(qū)域處理超大型低矮建筑物時，布料粒子可能會粘在屋頂或低矮植被上，從而可能會將某些地物點錯誤地分類為地面點。所以2種算法在城市等綜合復雜地形區(qū)域的濾波效果一般，優(yōu)缺點各異。較具代表性的有S22、S24，如圖5、圖6所示。

(a)S22 DSM (b)參考DEM

(a)S24 DSM (b)參考DEM

在山區(qū)等陡峭地形區(qū)域(S52、S53、S61、S71)2種算法都取得了較差的濾波效果。此外，這4組樣本數(shù)據都在鄉(xiāng)村區(qū)域，地面點數(shù)量遠大于非地面數(shù)量，地面上覆蓋的大多數(shù)是植被及低矮灌木。在用數(shù)學形態(tài)學濾波算法時，選取相對小的濾波窗口就能獲得相對好的濾波效果。而布料模擬濾波算法(CSF)，斜率是影響該濾波算法準確性的重要因素。斜坡上，模擬的布料會躺在斜坡上，但無法完美貼合粘附在地面上，即使進行后處理后，濾波精度的提升也相對有限。較具有代表性的有S52、S53，如圖7、圖8所示。

(a)S52 DSM (b)參考DEM

(a)S53 DSM (b)參考DEM

3 總結與展望

基于數(shù)學形態(tài)學的濾波算法原理簡單、效率高，其缺點在于濾波窗口大小的選取，這需要對地形有一定的先驗知識。布料模擬濾波算法(CSF)的優(yōu)點在于用戶定義的參數(shù)少且易于設置。但局限性在于將粒子運動的物理運動過程改為2個離散的步驟，當在城市等綜合地形復雜區(qū)域處理超大型低矮建筑物時，布料粒子可能會粘在屋頂或低矮植被上，從而可能會將某些非地面點錯誤地分類為地面點。2種算法優(yōu)缺點各異，在地形平坦區(qū)域都取得了良好的濾波效果，在城市等綜合復雜地形區(qū)域取得了相對較好的濾波效果，但從便利性來看，由于形態(tài)學濾波需要一定的先驗知識估算最大建筑物大小，從而確定濾波窗口大小的選取，布料模擬濾波在處理各類地形上更具便利性，設置參數(shù)少。2種算法在山區(qū)等陡峭地形區(qū)域取得了較差的濾波效果。這2種濾波算法和其他的濾波算法一樣，在陡峭地形區(qū)域的濾波精度低，濾波效果差。如何有效提升地形陡峭區(qū)域的濾波精度，這也是機載LiDAR點云濾波算法一直研究的熱點與難點。