王 凱，魯鐵定*，魯春陽(yáng)，劉 盈

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

0 引言

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

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

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

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

1 點(diǎn)云濾波原理

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

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

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

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

1.1.2 腐蝕 結(jié)構(gòu)元素b對(duì)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 閉運(yùn)算和開(kāi)運(yùn)算 由膨脹和腐蝕進(jìn)行組合。

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

f·b=(f。b)⊕b

(3)

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

開(kāi)運(yùn)算，用b開(kāi)運(yùn)算f，記為f。b，定義如式(4)：

f。b=(f。b)⊕b

(4)

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

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

wk=2kb+1

(5)

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

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

(6)

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

(7)

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

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

1.2 布料模擬濾波

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

圖1 布料模擬過(guò)程圖

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

(8)

圖2 剛性參數(shù)化

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

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

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

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

2)使原始LiDAR點(diǎn)云倒置。

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

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

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

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

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

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

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

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

2 實(shí)驗(yàn)結(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)行點(diǎn)云濾波實(shí)驗(yàn)，與Zhang[7]等2016年提出的布料模擬的濾波算法進(jìn)行實(shí)驗(yàn)，最后分別對(duì)2種濾波算法進(jìn)行精度評(píng)定并進(jìn)行適應(yīng)性分析。實(shí)驗(yàn)數(shù)據(jù)采用ISPRS(國(guó)際攝影測(cè)量和遙感學(xué)會(huì))數(shù)據(jù)，對(duì)2種算法的性能進(jìn)行了測(cè)試，該實(shí)驗(yàn)數(shù)據(jù)由Optech ALTM掃描儀收集，包括由不同地形環(huán)境組成的8個(gè)站點(diǎn)(命名為站點(diǎn)1～8)：4個(gè)城市站點(diǎn)和4個(gè)農(nóng)村站點(diǎn)。從站點(diǎn)1～7中選擇了15個(gè)具有不同地形特征和土地覆蓋類型的參考樣本，進(jìn)行濾波精度評(píng)定。由于缺乏參考數(shù)據(jù)，因此站點(diǎn)8被排除在外。對(duì)于每個(gè)實(shí)驗(yàn)樣本，參考數(shù)據(jù)是用先驗(yàn)知識(shí)和可用的航空?qǐng)D像手工識(shí)別生成的[13]，15組樣本數(shù)據(jù)地形特征如表1所示。

表1 15組樣本的地形特征

本文實(shí)驗(yàn)數(shù)據(jù)使用C++語(yǔ)言與點(diǎn)云庫(kù)PCL相結(jié)合，并將所得的點(diǎn)云數(shù)據(jù)在VS2013、點(diǎn)云處理軟件Cloudcompare進(jìn)行數(shù)據(jù)處理。為盡量避免不同噪聲點(diǎn)對(duì)濾波結(jié)果的影響，在進(jìn)行濾波實(shí)驗(yàn)前，15組樣本數(shù)據(jù)都通過(guò)手工進(jìn)行了去噪處理。此外，為定量分析2種算法的濾波性能，根據(jù)ISPRS規(guī)定的Ⅰ類誤差(將地面點(diǎn)錯(cuò)分為地物點(diǎn)的誤差)、Ⅱ類誤差(將地物點(diǎn)錯(cuò)分為地面點(diǎn)的誤差)和總誤差對(duì)實(shí)驗(yàn)結(jié)果進(jìn)行定量精度評(píng)定分析。15組樣本數(shù)據(jù)濾波實(shí)驗(yàn)結(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àn)知識(shí)去預(yù)估最大建筑物尺寸，以設(shè)定最大濾波窗口大小，而CSF的局限性在于已將粒子運(yùn)動(dòng)的物理運(yùn)動(dòng)過(guò)程改為2個(gè)離散的步驟，當(dāng)在城市等綜合地形復(fù)雜區(qū)域處理超大型低矮建筑物時(shí)，布料粒子可能會(huì)粘在屋頂或低矮植被上，從而可能會(huì)將某些地物點(diǎn)錯(cuò)誤地分類為地面點(diǎn)。所以2種算法在城市等綜合復(fù)雜地形區(qū)域的濾波效果一般，優(yōu)缺點(diǎn)各異。較具代表性的有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ū)域，地面點(diǎn)數(shù)量遠(yuǎn)大于非地面數(shù)量，地面上覆蓋的大多數(shù)是植被及低矮灌木。在用數(shù)學(xué)形態(tài)學(xué)濾波算法時(shí)，選取相對(duì)小的濾波窗口就能獲得相對(duì)好的濾波效果。而布料模擬濾波算法(CSF)，斜率是影響該濾波算法準(zhǔn)確性的重要因素。斜坡上，模擬的布料會(huì)躺在斜坡上，但無(wú)法完美貼合粘附在地面上，即使進(jìn)行后處理后，濾波精度的提升也相對(duì)有限。較具有代表性的有S52、S53，如圖7、圖8所示。

(a)S52 DSM (b)參考DEM

(a)S53 DSM (b)參考DEM

3 總結(jié)與展望

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