一種改進(jìn)的漸進(jìn)形態(tài)學(xué)濾波算法及應(yīng)用

徐 建 諸葛淑敏 朱潔瓊

(浙江省測(cè)繪科學(xué)技術(shù)研究院,浙江 杭州 310001)

0 引言

機(jī)載激光雷達(dá)(light detection and ranging,LiDAR)技術(shù)能夠在短時(shí)間內(nèi)獲取海量地表點(diǎn)云數(shù)據(jù),可對(duì)高精度三維地形信息進(jìn)行準(zhǔn)確還原,目前已廣泛應(yīng)用于城市三維建模、地面特征與土地提取覆蓋、高精度數(shù)字高程模型(digital elevation model,DEM)生產(chǎn)等領(lǐng)域[1-3]。機(jī)載激光LiDAR技術(shù)獲取地表點(diǎn)云數(shù)據(jù)主要由非地面點(diǎn)以及地面點(diǎn)組成,點(diǎn)云數(shù)據(jù)得以廣泛應(yīng)用的前提是對(duì)地面點(diǎn)與非地面點(diǎn)進(jìn)行分類(lèi),又稱(chēng)之為點(diǎn)云濾波。目前點(diǎn)云濾波算法多是基于地面點(diǎn)與非地面點(diǎn)的高程差,主要有移動(dòng)曲面擬合濾波算法、基于地形坡度的濾波算法以及數(shù)學(xué)形態(tài)學(xué)濾波算法等[4-5]。移動(dòng)曲面擬合濾波算法由張小紅等人[6]提出,該算法主要是通過(guò)二次曲面對(duì)點(diǎn)云表面進(jìn)行過(guò)濾;文獻(xiàn)[7]提出的基于地形坡度的濾波算法通過(guò)對(duì)于相鄰點(diǎn)高差與設(shè)置高差閾值實(shí)現(xiàn)點(diǎn)云屬性的判斷;文獻(xiàn)[8]提出的數(shù)學(xué)形態(tài)學(xué)濾波算法通過(guò)對(duì)固定窗口內(nèi)點(diǎn)云的開(kāi)運(yùn)算,根據(jù)高差變化實(shí)現(xiàn)點(diǎn)云屬性的判斷。經(jīng)典的數(shù)學(xué)形態(tài)學(xué)濾波算法是基于固定窗口,該算法能夠?qū)崿F(xiàn)固定大小地物的分類(lèi),為了提升該算法精度,在經(jīng)典數(shù)學(xué)形態(tài)學(xué)濾波算法的基礎(chǔ)上,文獻(xiàn)[9]通過(guò)改變窗口大小并設(shè)置高差閾值進(jìn)行算法改進(jìn),提出了漸進(jìn)形態(tài)學(xué)濾波算法。但是該算法假設(shè)地形坡度是連續(xù)的,對(duì)于復(fù)雜地形的點(diǎn)云濾波效果不佳,文獻(xiàn)[10]在漸進(jìn)形態(tài)學(xué)濾波算法的基礎(chǔ)上,將算法從一維擴(kuò)展至二維,較原有算法提升了地形的適應(yīng)性,但是需要設(shè)置大量參數(shù)。

針對(duì)現(xiàn)有漸進(jìn)形態(tài)學(xué)濾波算法對(duì)于地形復(fù)雜區(qū)域?yàn)V波時(shí)細(xì)節(jié)特征保留不完整、誤差較大等問(wèn)題,本文引入薄板樣條插值(thin plate spline,TPS)算法,建立組合濾波算法,通過(guò)TPS對(duì)不同大小窗口進(jìn)行處理,對(duì)膨脹運(yùn)算后的地形起伏度以及開(kāi)運(yùn)算后的高程差進(jìn)行計(jì)算,將地形起伏度與高程差作差,通過(guò)對(duì)比差值與設(shè)置閾值大小實(shí)現(xiàn)點(diǎn)云濾波。為檢驗(yàn)本文提出算法的有效性與優(yōu)越性,使用兩組不同地形條件的機(jī)載LiDAR點(diǎn)云數(shù)據(jù)進(jìn)行實(shí)驗(yàn),并與其他濾波算法的濾波結(jié)果進(jìn)行對(duì)比。

1 濾波原理

1.1 漸進(jìn)形態(tài)學(xué)濾波

數(shù)學(xué)形態(tài)學(xué)濾波算法由Lindenberger于1993年首次提出,目前主要用于密集匹配點(diǎn)云與點(diǎn)云濾波,數(shù)學(xué)形態(tài)學(xué)濾波算法的原理是:非地面點(diǎn)在經(jīng)過(guò)形態(tài)學(xué)計(jì)算以后,高程會(huì)產(chǎn)生較大變化,設(shè)置高差閾值,將高差大于高差閾值的點(diǎn)劃分為非地面點(diǎn),實(shí)現(xiàn)地面點(diǎn)、非地面點(diǎn)的分離,該算法濾波結(jié)果受濾波窗口的影響較大。數(shù)學(xué)形態(tài)學(xué)濾波包括4種運(yùn)算,分別為腐蝕運(yùn)算、膨脹運(yùn)算、開(kāi)運(yùn)算、閉運(yùn)算,其中膨脹運(yùn)算與腐蝕運(yùn)算是兩個(gè)相反的過(guò)程,膨脹運(yùn)算是計(jì)算窗口范圍內(nèi)最大高程值,腐蝕運(yùn)算是計(jì)算窗口范圍內(nèi)最小高程值,完成膨脹運(yùn)算與腐蝕運(yùn)算后,將兩個(gè)基礎(chǔ)算子組合成開(kāi)運(yùn)算與閉運(yùn)算。窗口大小設(shè)置的研究成為數(shù)學(xué)形態(tài)學(xué)濾波算法改進(jìn)的研究方向,如文獻(xiàn)[9]借助一個(gè)不斷改變窗口大小的開(kāi)算子,通過(guò)逐格網(wǎng)進(jìn)行非地形點(diǎn)濾除獲取地面點(diǎn),該方法中閾值的設(shè)置與窗口大小相對(duì)應(yīng),當(dāng)窗口大小大于最大建筑物尺寸時(shí)停止迭代。但是該方法使用于連續(xù)坡度區(qū)域,在復(fù)雜城市區(qū)域的地形特征濾波效果較差。

本文考慮優(yōu)化改進(jìn)經(jīng)典的漸進(jìn)式形態(tài),首先根據(jù)不同地物在高程上的突變,使用TPS對(duì)不同窗口大小的最低點(diǎn)進(jìn)行插值處理,得到近地表面DEM,其次計(jì)算DEM地形坡度并對(duì)其進(jìn)行膨脹運(yùn)算,最后計(jì)算點(diǎn)的高程與地形起伏度差值并與閾值對(duì)比得到最終濾波結(jié)果。

1.2 薄板樣條插值

薄板樣條插值是一種具有高魯棒性的數(shù)據(jù)插值方法,該插值算法是通過(guò)分布不規(guī)則的控制點(diǎn)構(gòu)建光滑表面,能夠滿足最小曲率面要求[11-12]。彎曲能量的表達(dá)式為

(1)

式中,S(x,y)為T(mén)PS近似表面函數(shù),表示為

(2)

式中,n為控制點(diǎn)數(shù)量;wi為控制點(diǎn)權(quán)重;a0、a1、a2是趨勢(shì)函數(shù)系數(shù);R為核函數(shù),表示為

R(r)=r2logr

(3)

式中,r為兩點(diǎn)的歐式距離;S(x,y)若要能夠進(jìn)行二次積分需滿足條件

(4)

式中,(xi,yi)為控制點(diǎn)坐標(biāo);根據(jù)式(4)得到線性方程

(5)

(6)

相比其他插值算法,TPS插值算法能夠更好地反映出地形高程變化,具有連續(xù)性、光滑性特征。此外,TPS插值算法具有更好的地形適應(yīng)性,對(duì)控制點(diǎn)分布的要求低,無(wú)論是簡(jiǎn)單地形,還是復(fù)雜地形均有較好的插值效果。

1.3 組合算法流程

結(jié)合多級(jí)形態(tài)學(xué)濾波算法與薄板樣條插值算法的優(yōu)勢(shì),建立一種新的組合濾波算法,使用薄板樣條插值算法對(duì)不同濾波窗口在不同層級(jí)上進(jìn)行插值處理,對(duì)插值后的結(jié)果進(jìn)行膨脹運(yùn)算并減去原始擬合面,得到地形起伏度,從而實(shí)現(xiàn)濾波精度的有效提升,地形起伏度的計(jì)算公式為[13]

σDEM=δDEM-ODEM

(7)

式中,σDEM為地形起伏度;δDEM為原始ODEM膨脹運(yùn)算結(jié)果。窗口尺寸的設(shè)置原則為從上到下依次減小,直至與最大建筑物尺寸近乎相等,本文提出的組合濾波算法技術(shù)路線如圖1所示。

圖1 本文組合濾波算法技術(shù)流程

本文提出的組合濾波算法的具體實(shí)施步驟為:

1)點(diǎn)云去噪。噪聲點(diǎn)主要為高程異常點(diǎn),根據(jù)高程可將噪聲點(diǎn)分為高位噪聲點(diǎn)與低位噪聲點(diǎn),主要通過(guò)兩個(gè)步驟進(jìn)行點(diǎn)云去噪,首先設(shè)置數(shù)量閾值,計(jì)算點(diǎn)云閾值距離內(nèi)的點(diǎn)云數(shù)量,剔除低于數(shù)量閾值噪聲點(diǎn),其次進(jìn)行人工校正[14-15]。

2)建立格網(wǎng)。對(duì)點(diǎn)云數(shù)據(jù)進(jìn)行格網(wǎng)劃分,其中索引行號(hào)為i,索引列號(hào)為j,格網(wǎng)內(nèi)高程最低值為該格網(wǎng)高程值。

3)形態(tài)學(xué)開(kāi)運(yùn)算。在當(dāng)前窗口尺寸下對(duì)柵格化點(diǎn)云數(shù)據(jù)進(jìn)行形態(tài)學(xué)開(kāi)運(yùn)算,記開(kāi)運(yùn)算前后的高差為H1。

4)地形起伏度計(jì)算。將柵格內(nèi)高程最小點(diǎn)指定給該柵格,使用TPS對(duì)所有高程最小點(diǎn)進(jìn)行插值,對(duì)插值后得到的擬合參考面進(jìn)行形態(tài)學(xué)膨脹運(yùn)算,將地形起伏度記為H2。

5)濾波判斷。獲取地形起伏度后,計(jì)算不同層級(jí)的地形坡度,同時(shí)計(jì)算處高差閾值DH,以增強(qiáng)在實(shí)際地形中的適應(yīng)性。計(jì)算H2與H1的差值,對(duì)比差值與DH大小,如果大于DH,將初始點(diǎn)視為地物點(diǎn)并刪剔除;如果小于DH,將初始點(diǎn)視為地物點(diǎn)并重復(fù)進(jìn)行上述步驟,直至窗口大小小于最小窗口尺寸。不同層級(jí)地形坡度公式為[16]

(8)

式中,FX、FY分別為水平、垂直梯度;nc為坡度變化常量;C為柵格尺寸;z為對(duì)應(yīng)參考高程值。

2 實(shí)驗(yàn)數(shù)據(jù)與精度評(píng)價(jià)指標(biāo)

2.1 實(shí)驗(yàn)數(shù)據(jù)

選取兩處具有不同地形特征的機(jī)載LiDAR點(diǎn)云數(shù)據(jù)進(jìn)行實(shí)驗(yàn),如圖2所示,其中實(shí)驗(yàn)區(qū)域1位于城市建筑物居民區(qū),區(qū)域內(nèi)建筑物較為密集,地形起伏較小,地形地物有建筑物、道路、陡坎、植被等,共含激光點(diǎn)1 135 868個(gè);區(qū)域2位于城市城市郊區(qū),區(qū)域內(nèi)建筑物較為稀疏,地形起伏較大,同樣有建筑物、道路、陡坎、植被等地物,共含激光點(diǎn)746 982個(gè)。

(a)區(qū)域1

2.2 精度評(píng)價(jià)指標(biāo)

使用國(guó)際攝影測(cè)量與遙感學(xué)會(huì)(International Society for Photogrammetry and Remote Sensing,ISPRS)在2003年提出的交叉表格形式對(duì)點(diǎn)云濾波結(jié)果進(jìn)行量化,如表1所示。

表1 交叉表格

表1中,a為正確分類(lèi)的地面點(diǎn)個(gè)數(shù);b為錯(cuò)誤分類(lèi)的非地面點(diǎn)個(gè)數(shù);c為錯(cuò)誤分類(lèi)的地面點(diǎn)個(gè)數(shù);d為正確分類(lèi)的非地面點(diǎn)個(gè)數(shù)[17]。

將Ⅰ類(lèi)誤差、Ⅱ類(lèi)誤差以及總誤差作為衡量點(diǎn)云濾波效果的評(píng)價(jià)指標(biāo)[18-20]

式中,n為激光點(diǎn)總數(shù);e為真實(shí)地面點(diǎn)個(gè)數(shù);f為真實(shí)非地面點(diǎn)個(gè)數(shù)。同時(shí)可通過(guò)Kappa系數(shù)[21]對(duì)濾波結(jié)果與參考數(shù)據(jù)的一致性進(jìn)行表征

(12)

式中,h為分類(lèi)為非地面點(diǎn)個(gè)數(shù);g為分類(lèi)為地面點(diǎn)的個(gè)數(shù)。

3 實(shí)驗(yàn)結(jié)果分析

使用點(diǎn)云庫(kù)PCL結(jié)合C++語(yǔ)言實(shí)現(xiàn)本文機(jī)載LiDAR點(diǎn)云濾波算法,處理時(shí)設(shè)置好相關(guān)參數(shù),將點(diǎn)云濾波結(jié)果與參考數(shù)據(jù)對(duì)比,實(shí)驗(yàn)區(qū)域1點(diǎn)云濾波結(jié)果與參考數(shù)據(jù)如圖3所示,實(shí)驗(yàn)區(qū)域2點(diǎn)云濾波結(jié)果與參考數(shù)據(jù)如圖4所示。

(a)濾波結(jié)果

(a)濾波結(jié)果

由圖3、圖4可知,本文提出機(jī)載點(diǎn)云濾波算法能夠有效濾除建筑物、植被等非地面點(diǎn),但是相較于參考數(shù)據(jù),本文算法沒(méi)有將非地面點(diǎn)完整濾除干凈,存在將少量非地面點(diǎn)分類(lèi)為地面點(diǎn),地面點(diǎn)分類(lèi)為非地面點(diǎn)的情況。

為對(duì)比得出本文提出算法的高效性與優(yōu)越性,使用經(jīng)典形態(tài)學(xué)濾波算法與經(jīng)典不規(guī)則三角網(wǎng)濾波算法進(jìn)行相同實(shí)驗(yàn)[22],并對(duì)比不同濾波算法實(shí)驗(yàn)結(jié)果。使用Ⅰ類(lèi)誤差、Ⅱ類(lèi)誤差以及總誤差定量評(píng)價(jià)不同算法濾波結(jié)果,精度統(tǒng)計(jì)結(jié)果如表2所示。

表2 不同點(diǎn)云濾波方法濾波結(jié)果精度統(tǒng)計(jì)

由表2可知,3種算法點(diǎn)云濾波結(jié)果的Ⅱ類(lèi)誤差要小于Ⅰ類(lèi)誤差,由于Ⅰ類(lèi)誤差表示將地面點(diǎn)分類(lèi)為非地面點(diǎn)的誤差,多發(fā)生于斜坡點(diǎn)云中,更易將該類(lèi)點(diǎn)云誤分為非地面點(diǎn)。同時(shí)可以看到,無(wú)論是對(duì)于實(shí)驗(yàn)區(qū)域1還是實(shí)驗(yàn)區(qū)域2,本文提出濾波算法濾波結(jié)果Ⅰ類(lèi)誤差、Ⅱ類(lèi)誤差以及總誤差均不同程度減少,在實(shí)驗(yàn)區(qū)域1,本文算法濾波結(jié)果較形態(tài)學(xué)濾波算法濾波結(jié)果的Ⅰ類(lèi)誤差、Ⅱ類(lèi)誤差以及總誤差分別減少了10.19%、7.98%、9.51%,較不規(guī)則TIN濾波算法濾波結(jié)果的Ⅰ類(lèi)誤差、Ⅱ類(lèi)誤差以及總誤差分別減少了9.35%、7.41%、7.36%;在實(shí)驗(yàn)區(qū)域2,本文算法濾波結(jié)果較形態(tài)學(xué)濾波算法濾波結(jié)果的Ⅰ類(lèi)誤差、Ⅱ類(lèi)誤差以及總誤差分別減少了7.44%、6.40%、6.85%,較不規(guī)則TIN濾波算法濾波結(jié)果的Ⅰ類(lèi)誤差、Ⅱ類(lèi)誤差以及總誤差分別減少了6.31%、5.77%、6.19%。同時(shí)可以看到實(shí)驗(yàn)區(qū)域1的地物更多,環(huán)境更為復(fù)雜,本文濾波算法仍然具有較好的濾波效果,本文算法對(duì)于實(shí)驗(yàn)區(qū)域1的誤差改善效果要優(yōu)于實(shí)驗(yàn)區(qū)域2,表明本文濾波算法受地形因素影響小。綜上所述,無(wú)論是在建筑物復(fù)雜區(qū)域還是地形復(fù)雜區(qū)域,本文算法均能夠有效改進(jìn)Ⅰ類(lèi)誤差、Ⅱ類(lèi)誤差以及總誤差,提升點(diǎn)云濾波精度,具有較高的地形適應(yīng)度,通過(guò)本文算法能夠?yàn)楂@取與實(shí)際地形貼合度更高的DEM做數(shù)據(jù)準(zhǔn)備。

4 結(jié)束語(yǔ)

點(diǎn)云濾波是激光LiDAR點(diǎn)云數(shù)據(jù)處理過(guò)程中的關(guān)鍵步驟,一直以來(lái)成為點(diǎn)云數(shù)據(jù)處理領(lǐng)域的研究熱點(diǎn),作為一種常用的激光LiDAR點(diǎn)云濾波算法,漸進(jìn)形態(tài)學(xué)濾波算法已經(jīng)逐漸無(wú)法滿足日益提升的點(diǎn)云數(shù)據(jù)處理要求,基于此,本文在漸進(jìn)形態(tài)學(xué)濾波算法的基礎(chǔ)上進(jìn)行算法改進(jìn),提出一種基于TPS插值的漸進(jìn)形態(tài)學(xué)濾波算法。該改進(jìn)算法在形態(tài)學(xué)算法優(yōu)勢(shì)的基礎(chǔ)上,結(jié)合TPS多級(jí)插值曲面擬合算法,使用TPS在漸進(jìn)形態(tài)學(xué)濾波算法不同尺寸的窗口下進(jìn)行處理,由上至下進(jìn)行迭代直至窗口尺寸小于預(yù)設(shè)尺寸。使用兩組不同地形條件的機(jī)載LiDAR點(diǎn)云數(shù)據(jù)進(jìn)行濾波實(shí)驗(yàn),結(jié)果表明較經(jīng)典形態(tài)學(xué)濾波算法與經(jīng)典不規(guī)則三角網(wǎng)濾波算法,本文算法能夠有效降低Ⅰ類(lèi)誤差、Ⅱ類(lèi)誤差以及總誤差,具有良好的可靠性與地形適用性。下一步將重點(diǎn)研究算法效率的提升以及參數(shù)設(shè)置的自適應(yīng)性,進(jìn)一步降低點(diǎn)云濾波的誤差。