徐 建 諸葛淑敏 朱潔瓊
(浙江省測繪科學(xué)技術(shù)研究院,浙江 杭州 310001)
機(jī)載激光雷達(dá)(light detection and ranging,LiDAR)技術(shù)能夠在短時間內(nèi)獲取海量地表點云數(shù)據(jù),可對高精度三維地形信息進(jìn)行準(zhǔn)確還原,目前已廣泛應(yīng)用于城市三維建模、地面特征與土地提取覆蓋、高精度數(shù)字高程模型(digital elevation model,DEM)生產(chǎn)等領(lǐng)域[1-3]。機(jī)載激光LiDAR技術(shù)獲取地表點云數(shù)據(jù)主要由非地面點以及地面點組成,點云數(shù)據(jù)得以廣泛應(yīng)用的前提是對地面點與非地面點進(jìn)行分類,又稱之為點云濾波。目前點云濾波算法多是基于地面點與非地面點的高程差,主要有移動曲面擬合濾波算法、基于地形坡度的濾波算法以及數(shù)學(xué)形態(tài)學(xué)濾波算法等[4-5]。移動曲面擬合濾波算法由張小紅等人[6]提出,該算法主要是通過二次曲面對點云表面進(jìn)行過濾;文獻(xiàn)[7]提出的基于地形坡度的濾波算法通過對于相鄰點高差與設(shè)置高差閾值實現(xiàn)點云屬性的判斷;文獻(xiàn)[8]提出的數(shù)學(xué)形態(tài)學(xué)濾波算法通過對固定窗口內(nèi)點云的開運(yùn)算,根據(jù)高差變化實現(xiàn)點云屬性的判斷。經(jīng)典的數(shù)學(xué)形態(tài)學(xué)濾波算法是基于固定窗口,該算法能夠?qū)崿F(xiàn)固定大小地物的分類,為了提升該算法精度,在經(jīng)典數(shù)學(xué)形態(tài)學(xué)濾波算法的基礎(chǔ)上,文獻(xiàn)[9]通過改變窗口大小并設(shè)置高差閾值進(jìn)行算法改進(jìn),提出了漸進(jìn)形態(tài)學(xué)濾波算法。但是該算法假設(shè)地形坡度是連續(xù)的,對于復(fù)雜地形的點云濾波效果不佳,文獻(xiàn)[10]在漸進(jìn)形態(tài)學(xué)濾波算法的基礎(chǔ)上,將算法從一維擴(kuò)展至二維,較原有算法提升了地形的適應(yīng)性,但是需要設(shè)置大量參數(shù)。
針對現(xiàn)有漸進(jìn)形態(tài)學(xué)濾波算法對于地形復(fù)雜區(qū)域濾波時細(xì)節(jié)特征保留不完整、誤差較大等問題,本文引入薄板樣條插值(thin plate spline,TPS)算法,建立組合濾波算法,通過TPS對不同大小窗口進(jìn)行處理,對膨脹運(yùn)算后的地形起伏度以及開運(yùn)算后的高程差進(jìn)行計算,將地形起伏度與高程差作差,通過對比差值與設(shè)置閾值大小實現(xiàn)點云濾波。為檢驗本文提出算法的有效性與優(yōu)越性,使用兩組不同地形條件的機(jī)載LiDAR點云數(shù)據(jù)進(jìn)行實驗,并與其他濾波算法的濾波結(jié)果進(jìn)行對比。
數(shù)學(xué)形態(tài)學(xué)濾波算法由Lindenberger于1993年首次提出,目前主要用于密集匹配點云與點云濾波,數(shù)學(xué)形態(tài)學(xué)濾波算法的原理是:非地面點在經(jīng)過形態(tài)學(xué)計算以后,高程會產(chǎn)生較大變化,設(shè)置高差閾值,將高差大于高差閾值的點劃分為非地面點,實現(xiàn)地面點、非地面點的分離,該算法濾波結(jié)果受濾波窗口的影響較大。數(shù)學(xué)形態(tài)學(xué)濾波包括4種運(yùn)算,分別為腐蝕運(yùn)算、膨脹運(yùn)算、開運(yùn)算、閉運(yùn)算,其中膨脹運(yùn)算與腐蝕運(yùn)算是兩個相反的過程,膨脹運(yùn)算是計算窗口范圍內(nèi)最大高程值,腐蝕運(yùn)算是計算窗口范圍內(nèi)最小高程值,完成膨脹運(yùn)算與腐蝕運(yùn)算后,將兩個基礎(chǔ)算子組合成開運(yùn)算與閉運(yùn)算。窗口大小設(shè)置的研究成為數(shù)學(xué)形態(tài)學(xué)濾波算法改進(jìn)的研究方向,如文獻(xiàn)[9]借助一個不斷改變窗口大小的開算子,通過逐格網(wǎng)進(jìn)行非地形點濾除獲取地面點,該方法中閾值的設(shè)置與窗口大小相對應(yīng),當(dāng)窗口大小大于最大建筑物尺寸時停止迭代。但是該方法使用于連續(xù)坡度區(qū)域,在復(fù)雜城市區(qū)域的地形特征濾波效果較差。
本文考慮優(yōu)化改進(jìn)經(jīng)典的漸進(jìn)式形態(tài),首先根據(jù)不同地物在高程上的突變,使用TPS對不同窗口大小的最低點進(jìn)行插值處理,得到近地表面DEM,其次計算DEM地形坡度并對其進(jìn)行膨脹運(yùn)算,最后計算點的高程與地形起伏度差值并與閾值對比得到最終濾波結(jié)果。
薄板樣條插值是一種具有高魯棒性的數(shù)據(jù)插值方法,該插值算法是通過分布不規(guī)則的控制點構(gòu)建光滑表面,能夠滿足最小曲率面要求[11-12]。彎曲能量的表達(dá)式為
(1)
式中,S(x,y)為TPS近似表面函數(shù),表示為
(2)
式中,n為控制點數(shù)量;wi為控制點權(quán)重;a0、a1、a2是趨勢函數(shù)系數(shù);R為核函數(shù),表示為
R(r)=r2logr
(3)
式中,r為兩點的歐式距離;S(x,y)若要能夠進(jìn)行二次積分需滿足條件
(4)
式中,(xi,yi)為控制點坐標(biāo);根據(jù)式(4)得到線性方程
(5)
(6)
相比其他插值算法,TPS插值算法能夠更好地反映出地形高程變化,具有連續(xù)性、光滑性特征。此外,TPS插值算法具有更好的地形適應(yīng)性,對控制點分布的要求低,無論是簡單地形,還是復(fù)雜地形均有較好的插值效果。
結(jié)合多級形態(tài)學(xué)濾波算法與薄板樣條插值算法的優(yōu)勢,建立一種新的組合濾波算法,使用薄板樣條插值算法對不同濾波窗口在不同層級上進(jìn)行插值處理,對插值后的結(jié)果進(jìn)行膨脹運(yùn)算并減去原始擬合面,得到地形起伏度,從而實現(xiàn)濾波精度的有效提升,地形起伏度的計算公式為[13]
σDEM=δDEM-ODEM
(7)
式中,σDEM為地形起伏度;δDEM為原始ODEM膨脹運(yùn)算結(jié)果。窗口尺寸的設(shè)置原則為從上到下依次減小,直至與最大建筑物尺寸近乎相等,本文提出的組合濾波算法技術(shù)路線如圖1所示。
圖1 本文組合濾波算法技術(shù)流程
本文提出的組合濾波算法的具體實施步驟為:
1)點云去噪。噪聲點主要為高程異常點,根據(jù)高程可將噪聲點分為高位噪聲點與低位噪聲點,主要通過兩個步驟進(jìn)行點云去噪,首先設(shè)置數(shù)量閾值,計算點云閾值距離內(nèi)的點云數(shù)量,剔除低于數(shù)量閾值噪聲點,其次進(jìn)行人工校正[14-15]。
2)建立格網(wǎng)。對點云數(shù)據(jù)進(jìn)行格網(wǎng)劃分,其中索引行號為i,索引列號為j,格網(wǎng)內(nèi)高程最低值為該格網(wǎng)高程值。
3)形態(tài)學(xué)開運(yùn)算。在當(dāng)前窗口尺寸下對柵格化點云數(shù)據(jù)進(jìn)行形態(tài)學(xué)開運(yùn)算,記開運(yùn)算前后的高差為H1。
4)地形起伏度計算。將柵格內(nèi)高程最小點指定給該柵格,使用TPS對所有高程最小點進(jìn)行插值,對插值后得到的擬合參考面進(jìn)行形態(tài)學(xué)膨脹運(yùn)算,將地形起伏度記為H2。
5)濾波判斷。獲取地形起伏度后,計算不同層級的地形坡度,同時計算處高差閾值DH,以增強(qiáng)在實際地形中的適應(yīng)性。計算H2與H1的差值,對比差值與DH大小,如果大于DH,將初始點視為地物點并刪剔除;如果小于DH,將初始點視為地物點并重復(fù)進(jìn)行上述步驟,直至窗口大小小于最小窗口尺寸。不同層級地形坡度公式為[16]
(8)
式中,FX、FY分別為水平、垂直梯度;nc為坡度變化常量;C為柵格尺寸;z為對應(yīng)參考高程值。
選取兩處具有不同地形特征的機(jī)載LiDAR點云數(shù)據(jù)進(jìn)行實驗,如圖2所示,其中實驗區(qū)域1位于城市建筑物居民區(qū),區(qū)域內(nèi)建筑物較為密集,地形起伏較小,地形地物有建筑物、道路、陡坎、植被等,共含激光點1 135 868個;區(qū)域2位于城市城市郊區(qū),區(qū)域內(nèi)建筑物較為稀疏,地形起伏較大,同樣有建筑物、道路、陡坎、植被等地物,共含激光點746 982個。
(a)區(qū)域1
使用國際攝影測量與遙感學(xué)會(International Society for Photogrammetry and Remote Sensing,ISPRS)在2003年提出的交叉表格形式對點云濾波結(jié)果進(jìn)行量化,如表1所示。
表1 交叉表格
表1中,a為正確分類的地面點個數(shù);b為錯誤分類的非地面點個數(shù);c為錯誤分類的地面點個數(shù);d為正確分類的非地面點個數(shù)[17]。
將Ⅰ類誤差、Ⅱ類誤差以及總誤差作為衡量點云濾波效果的評價指標(biāo)[18-20]
式中,n為激光點總數(shù);e為真實地面點個數(shù);f為真實非地面點個數(shù)。同時可通過Kappa系數(shù)[21]對濾波結(jié)果與參考數(shù)據(jù)的一致性進(jìn)行表征
(12)
式中,h為分類為非地面點個數(shù);g為分類為地面點的個數(shù)。
使用點云庫PCL結(jié)合C++語言實現(xiàn)本文機(jī)載LiDAR點云濾波算法,處理時設(shè)置好相關(guān)參數(shù),將點云濾波結(jié)果與參考數(shù)據(jù)對比,實驗區(qū)域1點云濾波結(jié)果與參考數(shù)據(jù)如圖3所示,實驗區(qū)域2點云濾波結(jié)果與參考數(shù)據(jù)如圖4所示。
(a)濾波結(jié)果
(a)濾波結(jié)果
由圖3、圖4可知,本文提出機(jī)載點云濾波算法能夠有效濾除建筑物、植被等非地面點,但是相較于參考數(shù)據(jù),本文算法沒有將非地面點完整濾除干凈,存在將少量非地面點分類為地面點,地面點分類為非地面點的情況。
為對比得出本文提出算法的高效性與優(yōu)越性,使用經(jīng)典形態(tài)學(xué)濾波算法與經(jīng)典不規(guī)則三角網(wǎng)濾波算法進(jìn)行相同實驗[22],并對比不同濾波算法實驗結(jié)果。使用Ⅰ類誤差、Ⅱ類誤差以及總誤差定量評價不同算法濾波結(jié)果,精度統(tǒng)計結(jié)果如表2所示。
表2 不同點云濾波方法濾波結(jié)果精度統(tǒng)計
由表2可知,3種算法點云濾波結(jié)果的Ⅱ類誤差要小于Ⅰ類誤差,由于Ⅰ類誤差表示將地面點分類為非地面點的誤差,多發(fā)生于斜坡點云中,更易將該類點云誤分為非地面點。同時可以看到,無論是對于實驗區(qū)域1還是實驗區(qū)域2,本文提出濾波算法濾波結(jié)果Ⅰ類誤差、Ⅱ類誤差以及總誤差均不同程度減少,在實驗區(qū)域1,本文算法濾波結(jié)果較形態(tài)學(xué)濾波算法濾波結(jié)果的Ⅰ類誤差、Ⅱ類誤差以及總誤差分別減少了10.19%、7.98%、9.51%,較不規(guī)則TIN濾波算法濾波結(jié)果的Ⅰ類誤差、Ⅱ類誤差以及總誤差分別減少了9.35%、7.41%、7.36%;在實驗區(qū)域2,本文算法濾波結(jié)果較形態(tài)學(xué)濾波算法濾波結(jié)果的Ⅰ類誤差、Ⅱ類誤差以及總誤差分別減少了7.44%、6.40%、6.85%,較不規(guī)則TIN濾波算法濾波結(jié)果的Ⅰ類誤差、Ⅱ類誤差以及總誤差分別減少了6.31%、5.77%、6.19%。同時可以看到實驗區(qū)域1的地物更多,環(huán)境更為復(fù)雜,本文濾波算法仍然具有較好的濾波效果,本文算法對于實驗區(qū)域1的誤差改善效果要優(yōu)于實驗區(qū)域2,表明本文濾波算法受地形因素影響小。綜上所述,無論是在建筑物復(fù)雜區(qū)域還是地形復(fù)雜區(qū)域,本文算法均能夠有效改進(jìn)Ⅰ類誤差、Ⅱ類誤差以及總誤差,提升點云濾波精度,具有較高的地形適應(yīng)度,通過本文算法能夠為獲取與實際地形貼合度更高的DEM做數(shù)據(jù)準(zhǔn)備。
點云濾波是激光LiDAR點云數(shù)據(jù)處理過程中的關(guān)鍵步驟,一直以來成為點云數(shù)據(jù)處理領(lǐng)域的研究熱點,作為一種常用的激光LiDAR點云濾波算法,漸進(jìn)形態(tài)學(xué)濾波算法已經(jīng)逐漸無法滿足日益提升的點云數(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)勢的基礎(chǔ)上,結(jié)合TPS多級插值曲面擬合算法,使用TPS在漸進(jìn)形態(tài)學(xué)濾波算法不同尺寸的窗口下進(jìn)行處理,由上至下進(jìn)行迭代直至窗口尺寸小于預(yù)設(shè)尺寸。使用兩組不同地形條件的機(jī)載LiDAR點云數(shù)據(jù)進(jìn)行濾波實驗,結(jié)果表明較經(jīng)典形態(tài)學(xué)濾波算法與經(jīng)典不規(guī)則三角網(wǎng)濾波算法,本文算法能夠有效降低Ⅰ類誤差、Ⅱ類誤差以及總誤差,具有良好的可靠性與地形適用性。下一步將重點研究算法效率的提升以及參數(shù)設(shè)置的自適應(yīng)性,進(jìn)一步降低點云濾波的誤差。