王逸超，吳小忠，陳仲偉，彭康博

(國(guó)網(wǎng)湖南省電力有限公司經(jīng)濟(jì)技術(shù)研究院，長(zhǎng)沙 410003)

0 引言

機(jī)載激光雷達(dá)(light detection and ranging，LiDAR)作為一種主動(dòng)遙感測(cè)量技術(shù)，它可以不受光照等條件的影響獲取高精度的密集點(diǎn)云，近年來(lái)得到了快速發(fā)展，已廣泛應(yīng)用于DTM生產(chǎn)、城市建模、電力線檢測(cè)等各種領(lǐng)域[1]。系統(tǒng)獲取的三維點(diǎn)云數(shù)據(jù)中除了地面點(diǎn)通常還包括建筑物、樹(shù)木、低矮植被、橋梁等非地面點(diǎn)，分離地面點(diǎn)和非地面點(diǎn)的過(guò)程稱為點(diǎn)云濾波，對(duì)于大多數(shù)應(yīng)用，點(diǎn)云濾波是機(jī)載LiDAR點(diǎn)云數(shù)據(jù)處理的核心步驟，會(huì)直接影響建模和特征提取的精度。點(diǎn)云濾波通常通過(guò)設(shè)定閾值來(lái)進(jìn)行地面點(diǎn)和非地面點(diǎn)的判斷，大多都假設(shè)地形坡度變化在一定范圍內(nèi)[2]。濾波閾值選擇的標(biāo)準(zhǔn)一般為坡度、高差距離、點(diǎn)到面的距離，數(shù)據(jù)組織形式可以基于單點(diǎn)的、格網(wǎng)、三角網(wǎng)等。

目前的濾波方法按照概念大致可分為3種，基于坡度的方法[3-4]、基于形態(tài)學(xué)的方法[5-6]和基于參數(shù)曲面的方法[7-9]。Vosselman提出了一種基于坡度的濾波方法，通過(guò)比較激光點(diǎn)和周圍鄰近點(diǎn)的坡度，若小于一定的坡度閾值，則將其當(dāng)做地面點(diǎn)[3]。Susaki提出一種自適應(yīng)坡度的濾波方法[4]。但基于坡度的方法隨著坡度的增大濾波的效果也會(huì)變差。Zhang等提出了一種漸進(jìn)形態(tài)學(xué)濾波方法，根據(jù)地形坡度在不同窗口尺寸內(nèi)設(shè)定不同的高差閾值[5]，但其閾值由固定的坡度計(jì)算得到，在地形突出的地方容易引起誤判。Hui等提出了基于多級(jí)克里金插值的形態(tài)學(xué)濾波方法[6]。Axelsson提出了漸進(jìn)加密三角網(wǎng)濾波方法，并應(yīng)用于TerraSolid軟件[7]。Zhao等提出了一種改進(jìn)的漸進(jìn)加密三角網(wǎng)濾波算法用于植被地區(qū)的濾波[8]。加密三角網(wǎng)濾波方法適用范圍較廣，但其實(shí)現(xiàn)復(fù)雜，參數(shù)設(shè)置多，時(shí)間復(fù)雜度較高。Kraus和Pfeifer提出一種迭代線性最小二乘插值算法用于植被地區(qū)的點(diǎn)云濾波，隨后擴(kuò)展了該算法用于城市地區(qū)建筑物和樹(shù)木的過(guò)濾[9]。此外，目前已有一些利用多回波和波形信息進(jìn)行濾波的一些研究，不同地物的反射回波具有不同的回波寬度，因此可利用回波寬度作為閾值或權(quán)重輔助進(jìn)行濾波[10-11]，但其在實(shí)際生產(chǎn)過(guò)程中的作用還待進(jìn)一步驗(yàn)證。

可以看出，大多數(shù)算法并非一次完成濾波操作，大多通過(guò)迭代的方式進(jìn)行。而相比只利用點(diǎn)云高差的算法，利用局部平面(曲面)的算法更加魯棒，濾波效果也更好，但是算法的效果在不同的地形上差別很大，通常都不太能夠處理復(fù)雜城市地區(qū)和地形不連續(xù)的地方[12]。另外，大多數(shù)算法都假設(shè)最低點(diǎn)為地面點(diǎn)，地面點(diǎn)若包含低位粗差點(diǎn)則會(huì)導(dǎo)致生成的DTM在局部下降。由于低位粗差點(diǎn)分布不規(guī)則，在無(wú)人工干預(yù)的情況下很難精確剔除所有的粗差點(diǎn)[6]。

針對(duì)上述問(wèn)題，本文提出了一種改進(jìn)的多級(jí)曲面濾波方法，該方法結(jié)合了形態(tài)學(xué)濾波和多級(jí)插值濾波的思想。在粗差點(diǎn)去除階段采用一種簡(jiǎn)單快速有效的基于坐標(biāo)直方圖統(tǒng)計(jì)的粗差點(diǎn)檢測(cè)方法。在種子點(diǎn)選擇階段，對(duì)初始種子點(diǎn)進(jìn)行形態(tài)學(xué)運(yùn)算進(jìn)一步篩選得到可靠的擬合點(diǎn)。通過(guò)建立不同分辨率的數(shù)據(jù)集，由粗到細(xì)逐級(jí)迭代得到最終的地面模型。

1 改進(jìn)的多級(jí)移動(dòng)曲面濾波算法

算法的總體流程圖如圖1所示。包括粗差點(diǎn)檢測(cè)、網(wǎng)格劃分、曲面擬合，自適應(yīng)閾值設(shè)置等步驟。在粗差點(diǎn)檢測(cè)階段本文提出了一種新的粗差點(diǎn)檢測(cè)方法，可以在保留正常點(diǎn)的同時(shí)快速有效的去除低位粗差點(diǎn)。在曲面擬合濾波階段，通過(guò)增加邊界和空白區(qū)域的插值處理以處理邊界和數(shù)據(jù)缺失地區(qū)的濾波運(yùn)算，對(duì)初始選擇種子點(diǎn)進(jìn)行一次形態(tài)學(xué)開(kāi)運(yùn)算得到更為可靠的地面點(diǎn)種子點(diǎn)。

圖1 總體流程圖

1.1 粗差點(diǎn)檢測(cè)

機(jī)載LiDAR測(cè)量得到的點(diǎn)云數(shù)據(jù)通常會(huì)包含粗差點(diǎn)，一般為高程異常點(diǎn)。高位粗差點(diǎn)通?？稍跒V波過(guò)程中去掉，而低位粗差點(diǎn)對(duì)濾波效果影響很大。本文采用一種統(tǒng)計(jì)高程直方圖的方式來(lái)快速有效去除低位粗差點(diǎn)，低位粗差點(diǎn)通常數(shù)量較少，且高程上不連續(xù)，通過(guò)設(shè)定一定的高程間距統(tǒng)計(jì)所有點(diǎn)的高程分布，可以計(jì)算得到點(diǎn)云的高程直方圖。直方圖橫坐標(biāo)為當(dāng)前塊內(nèi)點(diǎn)云的高程范圍，縱坐標(biāo)為落在各個(gè)區(qū)間上的點(diǎn)云個(gè)數(shù)。根據(jù)數(shù)據(jù)區(qū)域的地形起伏情況可以對(duì)數(shù)據(jù)進(jìn)行分塊處理，分別統(tǒng)計(jì)每塊內(nèi)點(diǎn)云的高程值直方圖，若數(shù)據(jù)區(qū)域內(nèi)有異常點(diǎn)，則其對(duì)應(yīng)的直方圖區(qū)間的點(diǎn)云個(gè)數(shù)數(shù)量較少，與其他正常點(diǎn)的區(qū)間分布存在明顯差異，設(shè)定低位粗差點(diǎn)的判別條件如下：若第i個(gè)間隔的點(diǎn)云個(gè)數(shù)小于一定閾值，且位置i小于最大點(diǎn)云個(gè)數(shù)所在間隔位置，則該間隔內(nèi)的點(diǎn)為低位粗差點(diǎn)。

圖2顯示了低位粗差點(diǎn)對(duì)高程直方圖分布的影響，圖2(a)顯示了無(wú)低位粗差點(diǎn)的高程直方圖分布，圖2(b)顯示了有粗差點(diǎn)的高程直方圖分布?？梢钥闯龅臀淮植铧c(diǎn)會(huì)使直方圖分布不連續(xù)，最大值位置偏移，通過(guò)設(shè)置合適的閾值可快速檢測(cè)去除低位粗差點(diǎn)。圖3顯示了使用本文的粗差點(diǎn)檢測(cè)方法與統(tǒng)計(jì)平均點(diǎn)間距法去除粗差點(diǎn)的效果對(duì)比。本文方法中高程間距值為0.5m，最小點(diǎn)云個(gè)數(shù)為10，統(tǒng)計(jì)平均點(diǎn)間距方法中臨近點(diǎn)搜索個(gè)數(shù)為6，標(biāo)準(zhǔn)差系數(shù)值為1.0。可以看出本文基于統(tǒng)計(jì)高程直方圖的粗差點(diǎn)檢測(cè)算法可以有效地去除低位粗差點(diǎn)，而統(tǒng)計(jì)平均點(diǎn)間距的方法在去除低位粗差點(diǎn)的同時(shí)也會(huì)去除部分有效的激光點(diǎn)。

圖2 不同區(qū)域點(diǎn)云高程直方圖分布

圖3 粗差點(diǎn)檢測(cè)方法對(duì)比

1.2 曲面擬合濾波

去除粗差點(diǎn)后，將點(diǎn)云數(shù)據(jù)按照網(wǎng)格大小劃分網(wǎng)格，每個(gè)網(wǎng)格內(nèi)存儲(chǔ)落在網(wǎng)格內(nèi)的點(diǎn)的索引。將網(wǎng)格內(nèi)的最低點(diǎn)作為初始種子點(diǎn)，通過(guò)形態(tài)學(xué)開(kāi)運(yùn)算和空白區(qū)域插值得到可靠種子點(diǎn)。算法首先設(shè)置較大的網(wǎng)格大小，得到一個(gè)粗略的地面模型，然后逐漸縮小網(wǎng)格大小不斷細(xì)化，得到最終的地面模型。曲面擬合濾波的步驟如下：

①建立網(wǎng)格：根據(jù)當(dāng)前網(wǎng)格大小建立網(wǎng)格，計(jì)算每個(gè)點(diǎn)的格網(wǎng)號(hào)并在格網(wǎng)內(nèi)存儲(chǔ)該點(diǎn)的索引。對(duì)每個(gè)當(dāng)前網(wǎng)格，尋找3×3的臨近窗口范圍內(nèi)的最低點(diǎn)，對(duì)每個(gè)網(wǎng)格內(nèi)的最低點(diǎn)進(jìn)行一次形態(tài)學(xué)開(kāi)運(yùn)算以保證其盡可能的為地面點(diǎn)。在邊界和空白窗口處進(jìn)行插值處理保證臨近窗口的完整性，以更好地?cái)M合地形，同時(shí)避免種子點(diǎn)個(gè)數(shù)太少導(dǎo)致曲面擬合失敗的情況。其中，高程值插值為周圍窗口內(nèi)的最低點(diǎn)，坐標(biāo)插值為該網(wǎng)格中心點(diǎn)坐標(biāo)位置。

②通過(guò)步驟1得到的9個(gè)地面點(diǎn)，建立二次多項(xiàng)式曲面方程，擬合方程如下：

z=a0+a1y+a2y2+a3x+a4xy+a5x2

(1)

二次多項(xiàng)式曲面有6個(gè)參數(shù)，觀測(cè)量通常多余6個(gè)，可通過(guò)建立誤差方程：

v=Bx-l

(2)

最小二乘可以求解得到曲面參數(shù)向量。

③計(jì)算高差閾值，地面點(diǎn)和非地面點(diǎn)的劃分是由激光點(diǎn)到擬合曲面的距離決定的，不同窗口尺寸和不同的地形變化需要設(shè)置不同的高差閾值。通常，地形起伏較大的區(qū)域，擬合曲面和真實(shí)地形有一定的偏差，因此需要設(shè)置較大的閾值，反之，則需要較小的閾值。另外窗口的尺寸越大，地形的起伏的可能性就越大，因此本文考慮窗口尺寸和曲面擬合均方根誤差的影響，窗口尺寸和均方根誤差越小，擬合的曲面越接近于真實(shí)地形，高差閾值就需要設(shè)置得相對(duì)較小，反之需要設(shè)置較為寬松的高差閾值。閾值設(shè)置如下：

T=Tmin+Wratio*(Tmax-Tmin)

(3)

式中：Wratio為窗口尺寸和高程方差歸一化加權(quán)之后的值；Tmin、Tmax為設(shè)置高差閾值的最小、最大值。

④更新網(wǎng)格大小，窗口尺寸的選擇對(duì)于濾波效果的影響同樣很大。窗口尺寸的選擇，一般有線性增長(zhǎng)和指數(shù)增長(zhǎng)2種方式，線性增長(zhǎng)方式可以很好地保持地形特征，但是會(huì)增加濾波的迭代次數(shù)?？紤]到點(diǎn)云一般數(shù)據(jù)量較大，以指數(shù)方式增長(zhǎng)窗口大小可以降低大范圍地區(qū)點(diǎn)云濾波的時(shí)間。通常，點(diǎn)云的平均點(diǎn)間距可以作為窗口尺寸的最小值，而窗口尺寸的最大值，則應(yīng)接近于建筑物的大小，為了保證窗口尺寸最大時(shí)，窗口內(nèi)的最低點(diǎn)為地面點(diǎn)，從而使大型建筑物在濾波過(guò)程中去除。濾波過(guò)程中每一級(jí)窗口的大小可以由初始值和增長(zhǎng)方式計(jì)算得到。

2 數(shù)據(jù)及實(shí)驗(yàn)分析

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

ISPRS委員會(huì)提供了用于濾波實(shí)驗(yàn)的基準(zhǔn)數(shù)據(jù)集，參考數(shù)據(jù)由人工濾波操作生成。數(shù)據(jù)集由Optech ALTM掃描儀在Vaihingen實(shí)驗(yàn)區(qū)和Stuttgart市區(qū)的8個(gè)地區(qū)采集得到。數(shù)據(jù)集共包含15個(gè)樣區(qū)，其中城市地區(qū)點(diǎn)云平均間距為1.0～1.5 m，農(nóng)村地區(qū)平均點(diǎn)間距為2.0～3.5 m。樣區(qū)地物類型豐富，包括大型建筑物、低矮植被、陡坡等，本文選擇其中的4個(gè)典型樣區(qū)進(jìn)行實(shí)驗(yàn)驗(yàn)證。樣區(qū)12帶有小型目標(biāo)物(汽車等)；樣區(qū)21具有狹窄的橋梁；樣區(qū)31地形不連續(xù)且?guī)в械臀淮植铧c(diǎn)；樣區(qū)52為帶有低矮植被，不連續(xù)且坡度較大的山地。

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

點(diǎn)云濾波效果的定量評(píng)價(jià)指標(biāo)通常是計(jì)算濾波的三類誤差，I類誤差、II類誤差和總誤差。其中I類誤差為所有地面點(diǎn)中地面點(diǎn)被當(dāng)做非地面點(diǎn)的比例，II類誤差為所有非地面點(diǎn)中被當(dāng)做地面點(diǎn)的比例，總誤差為分類錯(cuò)誤的點(diǎn)在總點(diǎn)數(shù)中的比例?？傉`差反映了濾波算法的總體性能。我們?cè)赩S2015平臺(tái)上利用c++語(yǔ)言實(shí)現(xiàn)了本文算法，并計(jì)算I類誤差、II類誤差和總誤差作為濾波精度的指標(biāo)，與其他方法進(jìn)行分析比較。

表1顯示了本文方法對(duì)不同測(cè)區(qū)數(shù)據(jù)濾波結(jié)果與其他8種經(jīng)典的濾波方法的對(duì)比結(jié)果。其他濾波方法結(jié)果數(shù)據(jù)來(lái)源于文獻(xiàn)[12]及其附錄。從表中可以看出，本文方法的濾波的總誤差整體較小，且優(yōu)于大部分方法，尤其在地形平坦的地區(qū)。不同地區(qū)的濾波精度差別較大，這跟測(cè)區(qū)地形的復(fù)雜度有關(guān)。

表1 濾波結(jié)果比較

其中測(cè)區(qū)12、21、31地形相對(duì)平坦，坡度變化不大，濾波效果較好，測(cè)區(qū)21、31的總誤差在1.5%左右，在所有濾波算法中誤差最低。隨著坡度增大，濾波誤差增大，這與現(xiàn)有的濾波方法的結(jié)論保持一致。測(cè)區(qū)52的總誤差為9.23%，相比于其他地形平坦的測(cè)區(qū)，總誤差增大了不少，但仍然好于大多數(shù)算法?？赡艿脑蚴怯捎谠诘匦纹鸱^大的地方，擬合曲面與實(shí)際的地形表面存在差異，且受窗口大小的影響較大，相對(duì)來(lái)說(shuō)，基于三角網(wǎng)的濾波效果要好一些，可能的原因是三角形比格網(wǎng)更適合表達(dá)起伏地形。

值得注意的是，沒(méi)有一種濾波算法在所有地形情況下濾波效果都達(dá)到最佳，用戶需要根據(jù)地形情況選擇合適的算法。另外參數(shù)的設(shè)置也會(huì)影響最終的濾波精度，實(shí)驗(yàn)表明，高差閾值的最小值通常在0.3 m到0.5 m之間，最大值在1.5 m到2 m之間。經(jīng)過(guò)兩到三次迭代通常會(huì)達(dá)到比較好的濾波效果。但是通過(guò)設(shè)置參數(shù)不可能完全消除誤差，I類誤差和II類誤差具有相關(guān)性，I類誤差減少，II類誤差就會(huì)增大，反之亦然。因此，用戶必須在減少I類誤差和II類誤差之間做出選擇，對(duì)于在后續(xù)的質(zhì)量控制過(guò)程中，II類錯(cuò)誤要比I類錯(cuò)誤更容易發(fā)現(xiàn)和被人工修復(fù)。但是對(duì)于濾波算法來(lái)說(shuō)，大多都會(huì)盡量減少II類誤差，因?yàn)橐恍┩蛊鸹蛘咻^低的非地面點(diǎn)對(duì)于自動(dòng)化地形重建的效果影響很大。

圖4顯示了參考數(shù)據(jù)和本文算法進(jìn)行點(diǎn)云濾波的實(shí)驗(yàn)結(jié)果及地形重建的效果圖，圖4(a)為標(biāo)準(zhǔn)的濾波參考數(shù)據(jù)集，其中紅色點(diǎn)為非地面點(diǎn)，藍(lán)色點(diǎn)為地面點(diǎn)。圖4(b)顯示了濾波后的地面點(diǎn)，對(duì)比圖4(a)可以看出濾波算法在去除了非地面點(diǎn)的同時(shí)保留了絕大部分的非地面點(diǎn)。圖4(c)顯示了濾波后的地面點(diǎn)通過(guò)反距離權(quán)重法利用Surf軟件進(jìn)行插值重建的地形效果，可以看出重建的地形效果接近于真實(shí)的地表形狀。

圖4 點(diǎn)云曲面擬合濾波

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

針對(duì)機(jī)載LiDAR點(diǎn)云濾波處理中存在的一些問(wèn)題，提出并實(shí)現(xiàn)了一種改進(jìn)的多級(jí)移動(dòng)曲面點(diǎn)云濾波方法，該方法結(jié)合了形態(tài)學(xué)運(yùn)算和曲面插值方法的優(yōu)點(diǎn)，兼顧了算法運(yùn)行速度和精度。實(shí)驗(yàn)表明，該方法能夠處理各種地形的數(shù)據(jù)，在平坦和地形起伏區(qū)域均能一定程度上提高點(diǎn)云的濾波精度，具有良好的魯棒性和適用性，可適用于大多數(shù)地形復(fù)雜的地區(qū)點(diǎn)云的濾波處理。本文的貢獻(xiàn)在于在數(shù)據(jù)預(yù)處理階段采用了一種基于統(tǒng)計(jì)直方圖的方法可以快速有效地去除低位粗差點(diǎn)，在地面種子點(diǎn)選取的過(guò)程中結(jié)合形態(tài)學(xué)開(kāi)運(yùn)算以獲取更加可靠的地面點(diǎn)，同時(shí)增加了對(duì)邊界和空白區(qū)域的插值處理以及自適應(yīng)窗口閾值的設(shè)置。如何進(jìn)一步提高算法在地形坡度劇烈變化地區(qū)的濾波精度是本文下一步研究的重點(diǎn)，利用多源數(shù)據(jù)和額外的輔助信息進(jìn)行濾波是進(jìn)一步提高濾波精度的一個(gè)研究方向。