一種機(jī)載LiDAR點(diǎn)云組合濾波方法的研究

姚杰，侯至群，朱大明，熊建華，付志濤，杜思雨，袁付松

(1.昆明理工大學(xué)國(guó)土資源工程學(xué)院，云南 昆明 650093； 2.昆明市測(cè)繪研究院，云南 昆明 650051； 3.昆明市子午環(huán)測(cè)繪咨詢服務(wù)有限責(zé)任公司，云南 昆明 650000)

1 引 言

機(jī)載LiDAR技術(shù)已廣泛應(yīng)用于地面數(shù)字高程模型(DEM)獲取[1]、道路提取[2]、電力線提取[3]、森林參數(shù)估計(jì)[4]、三維城市模型建立[5]等領(lǐng)域，濾波是這些應(yīng)用中一項(xiàng)必要的處理環(huán)節(jié)，是一個(gè)過濾地物點(diǎn)，保留地面點(diǎn)的過程，目的是為了突出重要信息。近些年來，許多專家學(xué)者對(duì)點(diǎn)云濾波算法進(jìn)行了針對(duì)性的研究，主要有基于鄰域點(diǎn)三維關(guān)系的濾波算法，如漸進(jìn)加密三角網(wǎng)濾波(PTD)[6]、坡度濾波[7]、形態(tài)學(xué)濾波[8]；基于平面擬合濾波算法，如最小二乘擬合濾波[9]、RANSC濾波[10]；基于物理過程模擬和濾波，如布料濾波等[11]。各種濾波算法經(jīng)過不斷改進(jìn)，在各自應(yīng)用領(lǐng)域都實(shí)現(xiàn)了較好的濾波效果。但是在面對(duì)大面積綜合復(fù)雜地形區(qū)域時(shí)濾波效果仍不理想。如在地形起伏較大區(qū)域，PTD濾波、CSF濾波、形態(tài)學(xué)濾波難以保留陡坡和山頂?shù)鹊匦翁卣餍畔?，坡度濾波能夠較好地保留地形特征信息但又難以過濾植被和房屋分布密集區(qū)域。針對(duì)這一問題，本文提出一種組合濾波方法，首先利用改進(jìn)漸進(jìn)加密三角網(wǎng)濾波算法對(duì)LiDAR點(diǎn)云數(shù)據(jù)進(jìn)行整體處理，再通過鄰域高差濾波對(duì)缺失的陡坡和山頂進(jìn)行補(bǔ)全，最后通過平面擬合濾波去除近地點(diǎn)噪聲，取得了較好的濾波效果，得到了完整且光滑的地面點(diǎn)云數(shù)據(jù)；并利用同名點(diǎn)較差的方法對(duì)地面點(diǎn)的整體精度進(jìn)行了評(píng)定。

2 組合濾波方法原理

組合濾波由三種濾波算法組合而成，本文重點(diǎn)介紹三種濾波算法的原理和處理結(jié)果，預(yù)處理、裁剪等基礎(chǔ)步驟不再詳述，算法流程如圖1所示。

圖1 組合濾波算法流程圖

2.1 漸進(jìn)加密三角網(wǎng)濾波算法

漸進(jìn)加密三角網(wǎng)算法(PTD)[12]最初是Axelsson提出的，算法通過預(yù)先生成一個(gè)稀疏的三角格網(wǎng)，再通過逐級(jí)迭代加密的方式進(jìn)行分類，直至所有點(diǎn)云分類完畢。Zhao等人對(duì)PTD算法進(jìn)行了改進(jìn)[13]，提出了IPTD濾波算法(Improved Progressive TIN Densification，IPTD)。該算法包括三個(gè)步驟：①利用形態(tài)學(xué)方法獲取潛在的地面種子點(diǎn)；②獲得準(zhǔn)確的地面種子點(diǎn)；③建立基于TIN的模型，對(duì)TIN進(jìn)行迭代致密化。量取最大建筑物尺寸作為初始格網(wǎng)大小，對(duì)點(diǎn)云數(shù)據(jù)進(jìn)行格網(wǎng)化處理，對(duì)于不含建筑物的點(diǎn)云數(shù)據(jù)，以默認(rèn)值作為初始格網(wǎng)大小。取格網(wǎng)內(nèi)的高程最低點(diǎn)作為起始種子點(diǎn)。遍歷所有待分類點(diǎn)，查詢各點(diǎn)水平面投影所落入的三角形，計(jì)算點(diǎn)到三角形底面的距離及點(diǎn)到三角形三個(gè)頂點(diǎn)與三角形所在平面形成夾角的最大值，將其分別與迭代距離與迭代角度進(jìn)行比較，如果小于對(duì)應(yīng)閾值，則將此點(diǎn)判定為地面點(diǎn)，并加入三角網(wǎng)中；如果大于對(duì)應(yīng)閾值，則去除。重復(fù)此過程，直至所有地面點(diǎn)分類完畢。

但是，按照格網(wǎng)最低點(diǎn)先向下逐級(jí)致密化再向上逐級(jí)致密化的算法原理，如圖2所示，山頂點(diǎn)和坡頂點(diǎn)需事先已知。如果初始網(wǎng)格中沒有包含這些點(diǎn)，那么陡坡點(diǎn)就不會(huì)被保留下來。盡管可以通過調(diào)整迭代角度和停止迭代三角網(wǎng)邊長(zhǎng)來逼近坡頂，但即使調(diào)整參數(shù)也不能完全還原陡坡點(diǎn)和山頂點(diǎn)，且過度的調(diào)整參數(shù)會(huì)使得其他類別的誤差增多。

圖2 IPTD向下致密化原理

2.2 鄰域高差濾波

鄰域高差濾波是利用種子點(diǎn)與鄰域點(diǎn)高差閾值來獲取小于閾值目標(biāo)點(diǎn)的原理進(jìn)行濾波的算法，通過控制高差閾值來控制點(diǎn)云的坡度大小。從坡度計(jì)算公式(式1)可知，坡度是通過高差與水平距離的比值計(jì)算所得，那么通過控制一定半徑范圍內(nèi)高差閾值，就可實(shí)現(xiàn)保留陡坡點(diǎn)云、過濾非地面點(diǎn)云的濾波效果，算法原理如圖3所示。根據(jù)鄰域高差濾波原理可知，在鄰域搜索半徑大于最小點(diǎn)間距的情況下，搜索半徑越小，則對(duì)非地面點(diǎn)的過濾效果越好。但是如果鄰域搜索半徑小于最小點(diǎn)間距，那么該算法會(huì)失去濾波功能，所以濾波結(jié)果會(huì)存在搜索半徑距離的近地點(diǎn)噪聲，并且大面積的建筑物屋頂點(diǎn)云因坡度滿足閾值也不會(huì)被過濾。綜上所述，鄰域高差濾波作用的區(qū)域有限，無法大面積應(yīng)用，但是其保留陡坡和山頂處點(diǎn)云的特點(diǎn)卻能夠很好地補(bǔ)充IPTD濾波算法的不足。

圖3 鄰域高差濾波原理

坡度α=arctan(高差/水平距離)

(1)

2.3 平面擬合濾波

平面擬合濾波算法的主要功能是對(duì)近地點(diǎn)噪聲進(jìn)行過濾，并使點(diǎn)云表面光滑。其基本步驟包括：①以種子點(diǎn)為參考原點(diǎn)，在原點(diǎn)周圍一定搜索半徑范圍內(nèi)搜索鄰域點(diǎn)，將鄰域點(diǎn)利用最小二乘方法進(jìn)行平面擬合；②計(jì)算點(diǎn)云到擬合平面的投影距離；③根據(jù)實(shí)際需求設(shè)置距離閾值，該閾值代表將點(diǎn)云向擬合平面調(diào)整的最大距離。平面擬合濾波可以讓點(diǎn)云表面更加光滑、與地面更加趨近一致；尤其是對(duì)表面粗糙的地面點(diǎn)處理后，點(diǎn)云排列會(huì)更加規(guī)則。因此，該方法能夠有效過濾IPTD濾波和鄰域高差濾波獲取結(jié)果中的近地點(diǎn)噪聲。

3 實(shí)驗(yàn)數(shù)據(jù)及實(shí)驗(yàn)內(nèi)容

3.1 研究區(qū)概況

實(shí)驗(yàn)區(qū)位于昆明市長(zhǎng)水機(jī)場(chǎng)附近，面積約為 3.75 km2，如圖4、圖5所示，點(diǎn)云數(shù)據(jù)由大棕熊100固定翼飛機(jī)搭載RIEGLVQ-1560i機(jī)載激光雷達(dá)設(shè)備于2020年采集所得?？傸c(diǎn)數(shù)4.9億個(gè)，最大點(diǎn)密度 2 616 per/m2，最小點(diǎn)密度 1 per/m2，平均點(diǎn)密度為 131 per/m2。本文使用LiDAR360、ArcGIS等軟件對(duì)點(diǎn)云數(shù)據(jù)進(jìn)行相關(guān)處理和分析。

圖4 研究區(qū)全景

圖5 局部細(xì)節(jié)展示

3.2 實(shí)驗(yàn)過程

(1)IPTD濾波

首先采用IPTD算法對(duì)原始LiDAR數(shù)據(jù)進(jìn)行初次濾波，分離地面點(diǎn)與非地面點(diǎn)。經(jīng)過反復(fù)實(shí)驗(yàn)，確定參數(shù)設(shè)置為：最大建筑物尺寸 400 m，最大值迭代角度10°，地形坡度角為88°。當(dāng)邊長(zhǎng)小于 5 m時(shí)減小迭代角度，小于 1 m時(shí)停止構(gòu)建三角形。經(jīng)過濾波后，得到地面點(diǎn)云如圖6所示。

圖6 IPRD濾波后得到的地面點(diǎn)云

從結(jié)果看出，IPTD點(diǎn)云濾波后陡坡和山頂處點(diǎn)云容易被當(dāng)作非地面點(diǎn)錯(cuò)誤過濾掉，如圖7所示。為了正確描述缺失處的地形特征，在小比例尺地形測(cè)繪時(shí)，目前通常做法是利用插值[14]方法進(jìn)行近似補(bǔ)全，但是，在大比例尺地形測(cè)繪時(shí)，插值方法會(huì)導(dǎo)致特征地形處出現(xiàn)失真。采用本文提出的方法，可以有效減少這類錯(cuò)誤的出現(xiàn)。

圖7 部分陡坡等地形變化特征點(diǎn)

(2)補(bǔ)全缺失點(diǎn)云并進(jìn)行鄰域高差濾波

通過對(duì)初次濾波數(shù)據(jù)的檢查，矢量裁剪了15塊缺失區(qū)域(如圖8所示)，一類為自然區(qū)域，如自然土坡、山體陡坎邊緣、山頂?shù)葏^(qū)域，如圖9中的2、15部分；另一類為人工區(qū)域，如道路護(hù)坡、擋土墻等區(qū)域，如圖9中的7、9部分。這些區(qū)域在構(gòu)建DEM時(shí)非常重要，因此有必要對(duì)這些區(qū)域進(jìn)行補(bǔ)全并進(jìn)行第二次濾波。

圖8 15塊缺失區(qū)域

圖9 部分缺失區(qū)域補(bǔ)全示意圖

根據(jù)每塊缺失地形的最大坡度計(jì)算高差閾值，經(jīng)鄰域高差濾波后不僅過濾了大部分非地面點(diǎn)數(shù)據(jù)，還完整的保留山坡的全部信息。但是由于大量的鄰域搜索計(jì)算，整體濾波計(jì)算時(shí)間較長(zhǎng)，且近地表面搜索半徑范圍內(nèi)存在一定的噪聲，如圖10所示。上文2.2中提到，這是由于鄰域高差濾波算法對(duì)小于搜索半徑的點(diǎn)不敏感所致。針對(duì)這一問題，提出利用平面擬合濾波的方法對(duì)近地點(diǎn)噪聲進(jìn)行過濾。

圖10 鄰域高差濾波結(jié)果表面噪聲

(3)平面擬合濾波

在前兩輪濾波之后，地面上大部分的地物點(diǎn)都已經(jīng)被過濾去除。針對(duì)近地點(diǎn)噪聲問題，提出利用平面擬合濾波算法，首先通過一定范圍內(nèi)的鄰域點(diǎn)擬合平面，并計(jì)算鄰域內(nèi)所有點(diǎn)到擬合平面的距離；再根據(jù)實(shí)際噪聲高度，設(shè)置相應(yīng)的距離閾值，修改噪聲點(diǎn)高程；最終得到了光滑的地面點(diǎn)數(shù)據(jù)，如圖11所示。

圖11 平面擬合濾波平滑表面

經(jīng)過上述三輪濾波處理后，得到了平滑的地面點(diǎn)數(shù)據(jù)。但是，補(bǔ)全后的地面數(shù)據(jù)是分散的，故需將15處缺失區(qū)域地面點(diǎn)和初次濾波得到的地面點(diǎn)數(shù)據(jù)進(jìn)行融合，得到完整的地面點(diǎn)數(shù)據(jù)。因與初次濾波得到地面數(shù)據(jù)運(yùn)用的算法不同，融合后會(huì)出現(xiàn)點(diǎn)云密度不均勻現(xiàn)象，如圖12所示，再利用重采樣的方法即可獲得均勻且完整的地面點(diǎn)云數(shù)據(jù)。

圖12 補(bǔ)全區(qū)域與初始地面點(diǎn)密度不一致

3.3 構(gòu)建DEM

因LiDAR激光穿透力有限，建筑物和茂密植被區(qū)域仍然會(huì)出現(xiàn)空洞，按照《激光雷達(dá)數(shù)據(jù)處理技術(shù)規(guī)范》CH/T 8023-2011中 1∶500地形圖的DEM格網(wǎng)間距為 0.5m×0.5 m的要求設(shè)置，利用反距離權(quán)重插值[15]的方法內(nèi)插出空洞區(qū)高程值，格式為TIFF文件。

4 精度評(píng)定

點(diǎn)云濾波評(píng)價(jià)方法有定量分析、定性分析等評(píng)定方法，也有通過構(gòu)建DEM后選點(diǎn)與實(shí)測(cè)點(diǎn)進(jìn)行交叉驗(yàn)證、同名點(diǎn)較差等評(píng)定方法。因本實(shí)驗(yàn)數(shù)據(jù)為實(shí)測(cè)數(shù)據(jù)，故采用構(gòu)建DEM后采取同名點(diǎn)較差的方法更加精確可靠。

為了對(duì)整體測(cè)區(qū)精度進(jìn)行準(zhǔn)確評(píng)定，將測(cè)區(qū)劃分為平地、丘陵、山地等區(qū)域，在不同區(qū)域內(nèi)均勻選取一定數(shù)量的高程點(diǎn)。通過計(jì)算選取的DEM高程點(diǎn)與實(shí)測(cè)檢核點(diǎn)的高程殘差△h的中誤差來評(píng)定DEM模型的精度，從而達(dá)到評(píng)定組合濾波結(jié)果精度的目的。實(shí)驗(yàn)區(qū)檢核點(diǎn)由GNSS系統(tǒng)采集，經(jīng)過平面、高程基準(zhǔn)轉(zhuǎn)換統(tǒng)一后，在DEM相同位置選取同名點(diǎn)共189個(gè)，如圖13所示，按照中誤差計(jì)算公式(式2)計(jì)算。

圖13 同名點(diǎn)分布

式中，mh為公共點(diǎn)高程中誤差，單位為(m)；△hi為檢核點(diǎn)高程較差，單位為米(m)；n為參與精度評(píng)定的點(diǎn)數(shù)為參與精度評(píng)定的點(diǎn)數(shù)。經(jīng)過計(jì)算，不同區(qū)域同名點(diǎn)較差的中誤差如表1所示，其中平地、丘陵地和山地的精度都遠(yuǎn)高于標(biāo)準(zhǔn)限差值，所得DEM數(shù)據(jù)準(zhǔn)確。

表1 同名點(diǎn)較差中誤差

從誤差結(jié)果分布圖14可以看出：①總體誤差值均優(yōu)于規(guī)范限差并呈正態(tài)分布；②平地和山地中誤差較小，平地精度又優(yōu)于山地，平地主要為農(nóng)田和空地，山地主要為植被；丘陵地中誤差較大，丘陵地中主要地物為廠房和居民房。丘陵地區(qū)存在較大的誤差值，主要是本文提出的組合濾波效果在房屋分布密集區(qū)域并沒有獨(dú)特的優(yōu)勢(shì)和針對(duì)性改進(jìn)，并且構(gòu)建居民區(qū)和廠區(qū)DEM時(shí)運(yùn)用了內(nèi)插算法，且內(nèi)插面積較大，也是導(dǎo)致誤差值較大的原因之一。因此，本文提出的組合濾波算法在較多建筑物區(qū)域處理效果方面還有待進(jìn)一步提高。

圖14 誤差分布

5 結(jié) 論

點(diǎn)云濾波是LiDAR點(diǎn)云數(shù)據(jù)處理過程中最重要的一環(huán)，從地物影響因素的角度出發(fā)，不僅要考慮如地形、建筑物、植被等因素影響，還要考慮多種因素的聯(lián)合影響，在單一濾波算法中加強(qiáng)對(duì)一種因素影響的控制同時(shí)會(huì)減弱對(duì)另一種影響因素的控制，難以適應(yīng)大面積地形復(fù)雜的區(qū)域。本文提出一種組合濾波的方案，首先對(duì)初始LiDAR點(diǎn)云數(shù)據(jù)利用IPTD濾波算法去除了大部分的地物點(diǎn)，通過鄰域高差濾波的方法對(duì)IPTD濾波在陡坡、山頂?shù)葏^(qū)域的缺失點(diǎn)進(jìn)行了補(bǔ)全，最后通過平面擬合濾波的方法對(duì)數(shù)據(jù)進(jìn)行了平滑，消除了近地點(diǎn)噪聲；經(jīng)過融合、重采樣、內(nèi)插等處理之后，制作了實(shí)驗(yàn)區(qū) 1∶500地面DEM，通過同名點(diǎn)較差的方法進(jìn)行精度評(píng)定，精度均滿足規(guī)范要求，證明組合濾波是一種有效可靠的濾波算法，并取得了如下結(jié)論：

(1)相比以往單一濾波方法，組合濾波方法的適應(yīng)性更強(qiáng)，可以做到各濾波算法之間取長(zhǎng)補(bǔ)短，有效降低了單一濾波算法的不足對(duì)結(jié)果的影響。

(2)組合濾波方法適用于大面積綜合場(chǎng)景下的LiDAR數(shù)據(jù)，對(duì)自然和人工陡坡、山頂?shù)鹊匦蔚臑V波效果好；但是也存在處理步驟多、智能化程度低、對(duì)建筑物密集地區(qū)的過濾能力不足等問題。