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

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

0 引言

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

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

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

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

1 點云濾波原理

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

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

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

1.1.1 膨脹 結(jié)構(gòu)元素b對f進(jìn)行灰度膨脹記為(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的定義域。膨脹運算是在由結(jié)構(gòu)元素確定的領(lǐng)域中，選取(f⊕b)的最大值。

1.1.2 腐蝕 結(jié)構(gòu)元素b對f進(jìn)行灰度腐蝕記為(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的定義域。腐蝕是在由結(jié)構(gòu)元素確定的領(lǐng)域中選取(fΘb)的最小值。

1.1.3 閉運算和開運算 由膨脹和腐蝕進(jìn)行組合。

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

f·b=(f。b)⊕b

(3)

將結(jié)構(gòu)元素內(nèi)的所有數(shù)據(jù)填充為區(qū)域最大值，然后再填充為該區(qū)域最小值。

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

f。b=(f。b)⊕b

(4)

將結(jié)構(gòu)元素內(nèi)的所有數(shù)據(jù)填充為區(qū)域最小值，然后再填充為該區(qū)域最大值。

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

wk=2kb+1

(5)

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

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

(6)

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

(7)

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

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

1.2 布料模擬濾波

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

圖1 布料模擬過程圖

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

(8)

圖2 剛性參數(shù)化

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

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

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

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

2)使原始LiDAR點云倒置。

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

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

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

6)計算每個布料粒子受內(nèi)力影響產(chǎn)生的位移距離。

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

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

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

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

2 實驗結(jié)果與分析

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

表1 15組樣本的地形特征

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

表2 2種算法濾波結(jié)果

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

(a)S21 DSM (b)參考DEM

(a)S21 DSM (b)參考DEM

2種算法在城市等綜合復(fù)雜地形區(qū)域(S11、S22、S23、S24、S41)濾波效果一般，由于數(shù)學(xué)形態(tài)學(xué)濾波算法需要一定的先驗知識去預(yù)估最大建筑物尺寸，以設(shè)定最大濾波窗口大小，而CSF的局限性在于已將粒子運動的物理運動過程改為2個離散的步驟，當(dāng)在城市等綜合地形復(fù)雜區(qū)域處理超大型低矮建筑物時，布料粒子可能會粘在屋頂或低矮植被上，從而可能會將某些地物點錯誤地分類為地面點。所以2種算法在城市等綜合復(fù)雜地形區(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ù)據(jù)都在鄉(xiāng)村區(qū)域，地面點數(shù)量遠(yuǎn)大于非地面數(shù)量，地面上覆蓋的大多數(shù)是植被及低矮灌木。在用數(shù)學(xué)形態(tài)學(xué)濾波算法時，選取相對小的濾波窗口就能獲得相對好的濾波效果。而布料模擬濾波算法(CSF)，斜率是影響該濾波算法準(zhǔn)確性的重要因素。斜坡上，模擬的布料會躺在斜坡上，但無法完美貼合粘附在地面上，即使進(jìn)行后處理后，濾波精度的提升也相對有限。較具有代表性的有S52、S53，如圖7、圖8所示。

(a)S52 DSM (b)參考DEM

(a)S53 DSM (b)參考DEM

3 總結(jié)與展望

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