姚杰,侯至群,朱大明,熊建華,付志濤,杜思雨,袁付松
(1.昆明理工大學國土資源工程學院,云南 昆明 650093; 2.昆明市測繪研究院,云南 昆明 650051; 3.昆明市子午環(huán)測繪咨詢服務有限責任公司,云南 昆明 650000)
機載LiDAR技術已廣泛應用于地面數字高程模型(DEM)獲取[1]、道路提取[2]、電力線提取[3]、森林參數估計[4]、三維城市模型建立[5]等領域,濾波是這些應用中一項必要的處理環(huán)節(jié),是一個過濾地物點,保留地面點的過程,目的是為了突出重要信息。近些年來,許多專家學者對點云濾波算法進行了針對性的研究,主要有基于鄰域點三維關系的濾波算法,如漸進加密三角網濾波(PTD)[6]、坡度濾波[7]、形態(tài)學濾波[8];基于平面擬合濾波算法,如最小二乘擬合濾波[9]、RANSC濾波[10];基于物理過程模擬和濾波,如布料濾波等[11]。各種濾波算法經過不斷改進,在各自應用領域都實現了較好的濾波效果。但是在面對大面積綜合復雜地形區(qū)域時濾波效果仍不理想。如在地形起伏較大區(qū)域,PTD濾波、CSF濾波、形態(tài)學濾波難以保留陡坡和山頂等地形特征信息,坡度濾波能夠較好地保留地形特征信息但又難以過濾植被和房屋分布密集區(qū)域。針對這一問題,本文提出一種組合濾波方法,首先利用改進漸進加密三角網濾波算法對LiDAR點云數據進行整體處理,再通過鄰域高差濾波對缺失的陡坡和山頂進行補全,最后通過平面擬合濾波去除近地點噪聲,取得了較好的濾波效果,得到了完整且光滑的地面點云數據;并利用同名點較差的方法對地面點的整體精度進行了評定。
組合濾波由三種濾波算法組合而成,本文重點介紹三種濾波算法的原理和處理結果,預處理、裁剪等基礎步驟不再詳述,算法流程如圖1所示。
圖1 組合濾波算法流程圖
漸進加密三角網算法(PTD)[12]最初是Axelsson提出的,算法通過預先生成一個稀疏的三角格網,再通過逐級迭代加密的方式進行分類,直至所有點云分類完畢。Zhao等人對PTD算法進行了改進[13],提出了IPTD濾波算法(Improved Progressive TIN Densification,IPTD)。該算法包括三個步驟:①利用形態(tài)學方法獲取潛在的地面種子點;②獲得準確的地面種子點;③建立基于TIN的模型,對TIN進行迭代致密化。量取最大建筑物尺寸作為初始格網大小,對點云數據進行格網化處理,對于不含建筑物的點云數據,以默認值作為初始格網大小。取格網內的高程最低點作為起始種子點。遍歷所有待分類點,查詢各點水平面投影所落入的三角形,計算點到三角形底面的距離及點到三角形三個頂點與三角形所在平面形成夾角的最大值,將其分別與迭代距離與迭代角度進行比較,如果小于對應閾值,則將此點判定為地面點,并加入三角網中;如果大于對應閾值,則去除。重復此過程,直至所有地面點分類完畢。
但是,按照格網最低點先向下逐級致密化再向上逐級致密化的算法原理,如圖2所示,山頂點和坡頂點需事先已知。如果初始網格中沒有包含這些點,那么陡坡點就不會被保留下來。盡管可以通過調整迭代角度和停止迭代三角網邊長來逼近坡頂,但即使調整參數也不能完全還原陡坡點和山頂點,且過度的調整參數會使得其他類別的誤差增多。
圖2 IPTD向下致密化原理
鄰域高差濾波是利用種子點與鄰域點高差閾值來獲取小于閾值目標點的原理進行濾波的算法,通過控制高差閾值來控制點云的坡度大小。從坡度計算公式(式1)可知,坡度是通過高差與水平距離的比值計算所得,那么通過控制一定半徑范圍內高差閾值,就可實現保留陡坡點云、過濾非地面點云的濾波效果,算法原理如圖3所示。根據鄰域高差濾波原理可知,在鄰域搜索半徑大于最小點間距的情況下,搜索半徑越小,則對非地面點的過濾效果越好。但是如果鄰域搜索半徑小于最小點間距,那么該算法會失去濾波功能,所以濾波結果會存在搜索半徑距離的近地點噪聲,并且大面積的建筑物屋頂點云因坡度滿足閾值也不會被過濾。綜上所述,鄰域高差濾波作用的區(qū)域有限,無法大面積應用,但是其保留陡坡和山頂處點云的特點卻能夠很好地補充IPTD濾波算法的不足。
圖3 鄰域高差濾波原理
坡度α=arctan(高差/水平距離)
(1)
平面擬合濾波算法的主要功能是對近地點噪聲進行過濾,并使點云表面光滑。其基本步驟包括:①以種子點為參考原點,在原點周圍一定搜索半徑范圍內搜索鄰域點,將鄰域點利用最小二乘方法進行平面擬合;②計算點云到擬合平面的投影距離;③根據實際需求設置距離閾值,該閾值代表將點云向擬合平面調整的最大距離。平面擬合濾波可以讓點云表面更加光滑、與地面更加趨近一致;尤其是對表面粗糙的地面點處理后,點云排列會更加規(guī)則。因此,該方法能夠有效過濾IPTD濾波和鄰域高差濾波獲取結果中的近地點噪聲。
實驗區(qū)位于昆明市長水機場附近,面積約為 3.75 km2,如圖4、圖5所示,點云數據由大棕熊100固定翼飛機搭載RIEGLVQ-1560i機載激光雷達設備于2020年采集所得。總點數4.9億個,最大點密度 2 616 per/m2,最小點密度 1 per/m2,平均點密度為 131 per/m2。本文使用LiDAR360、ArcGIS等軟件對點云數據進行相關處理和分析。
圖4 研究區(qū)全景
圖5 局部細節(jié)展示
(1)IPTD濾波
首先采用IPTD算法對原始LiDAR數據進行初次濾波,分離地面點與非地面點。經過反復實驗,確定參數設置為:最大建筑物尺寸 400 m,最大值迭代角度10°,地形坡度角為88°。當邊長小于 5 m時減小迭代角度,小于 1 m時停止構建三角形。經過濾波后,得到地面點云如圖6所示。
圖6 IPRD濾波后得到的地面點云
從結果看出,IPTD點云濾波后陡坡和山頂處點云容易被當作非地面點錯誤過濾掉,如圖7所示。為了正確描述缺失處的地形特征,在小比例尺地形測繪時,目前通常做法是利用插值[14]方法進行近似補全,但是,在大比例尺地形測繪時,插值方法會導致特征地形處出現失真。采用本文提出的方法,可以有效減少這類錯誤的出現。
圖7 部分陡坡等地形變化特征點
(2)補全缺失點云并進行鄰域高差濾波
通過對初次濾波數據的檢查,矢量裁剪了15塊缺失區(qū)域(如圖8所示),一類為自然區(qū)域,如自然土坡、山體陡坎邊緣、山頂等區(qū)域,如圖9中的2、15部分;另一類為人工區(qū)域,如道路護坡、擋土墻等區(qū)域,如圖9中的7、9部分。這些區(qū)域在構建DEM時非常重要,因此有必要對這些區(qū)域進行補全并進行第二次濾波。
圖8 15塊缺失區(qū)域
圖9 部分缺失區(qū)域補全示意圖
根據每塊缺失地形的最大坡度計算高差閾值,經鄰域高差濾波后不僅過濾了大部分非地面點數據,還完整的保留山坡的全部信息。但是由于大量的鄰域搜索計算,整體濾波計算時間較長,且近地表面搜索半徑范圍內存在一定的噪聲,如圖10所示。上文2.2中提到,這是由于鄰域高差濾波算法對小于搜索半徑的點不敏感所致。針對這一問題,提出利用平面擬合濾波的方法對近地點噪聲進行過濾。
圖10 鄰域高差濾波結果表面噪聲
(3)平面擬合濾波
在前兩輪濾波之后,地面上大部分的地物點都已經被過濾去除。針對近地點噪聲問題,提出利用平面擬合濾波算法,首先通過一定范圍內的鄰域點擬合平面,并計算鄰域內所有點到擬合平面的距離;再根據實際噪聲高度,設置相應的距離閾值,修改噪聲點高程;最終得到了光滑的地面點數據,如圖11所示。
圖11 平面擬合濾波平滑表面
經過上述三輪濾波處理后,得到了平滑的地面點數據。但是,補全后的地面數據是分散的,故需將15處缺失區(qū)域地面點和初次濾波得到的地面點數據進行融合,得到完整的地面點數據。因與初次濾波得到地面數據運用的算法不同,融合后會出現點云密度不均勻現象,如圖12所示,再利用重采樣的方法即可獲得均勻且完整的地面點云數據。
圖12 補全區(qū)域與初始地面點密度不一致
因LiDAR激光穿透力有限,建筑物和茂密植被區(qū)域仍然會出現空洞,按照《激光雷達數據處理技術規(guī)范》CH/T 8023-2011中 1∶500地形圖的DEM格網間距為 0.5m×0.5 m的要求設置,利用反距離權重插值[15]的方法內插出空洞區(qū)高程值,格式為TIFF文件。
點云濾波評價方法有定量分析、定性分析等評定方法,也有通過構建DEM后選點與實測點進行交叉驗證、同名點較差等評定方法。因本實驗數據為實測數據,故采用構建DEM后采取同名點較差的方法更加精確可靠。
為了對整體測區(qū)精度進行準確評定,將測區(qū)劃分為平地、丘陵、山地等區(qū)域,在不同區(qū)域內均勻選取一定數量的高程點。通過計算選取的DEM高程點與實測檢核點的高程殘差△h的中誤差來評定DEM模型的精度,從而達到評定組合濾波結果精度的目的。實驗區(qū)檢核點由GNSS系統(tǒng)采集,經過平面、高程基準轉換統(tǒng)一后,在DEM相同位置選取同名點共189個,如圖13所示,按照中誤差計算公式(式2)計算。
圖13 同名點分布
式中,mh為公共點高程中誤差,單位為(m);△hi為檢核點高程較差,單位為米(m);n為參與精度評定的點數為參與精度評定的點數。經過計算,不同區(qū)域同名點較差的中誤差如表1所示,其中平地、丘陵地和山地的精度都遠高于標準限差值,所得DEM數據準確。
表1 同名點較差中誤差
從誤差結果分布圖14可以看出:①總體誤差值均優(yōu)于規(guī)范限差并呈正態(tài)分布;②平地和山地中誤差較小,平地精度又優(yōu)于山地,平地主要為農田和空地,山地主要為植被;丘陵地中誤差較大,丘陵地中主要地物為廠房和居民房。丘陵地區(qū)存在較大的誤差值,主要是本文提出的組合濾波效果在房屋分布密集區(qū)域并沒有獨特的優(yōu)勢和針對性改進,并且構建居民區(qū)和廠區(qū)DEM時運用了內插算法,且內插面積較大,也是導致誤差值較大的原因之一。因此,本文提出的組合濾波算法在較多建筑物區(qū)域處理效果方面還有待進一步提高。
圖14 誤差分布
點云濾波是LiDAR點云數據處理過程中最重要的一環(huán),從地物影響因素的角度出發(fā),不僅要考慮如地形、建筑物、植被等因素影響,還要考慮多種因素的聯合影響,在單一濾波算法中加強對一種因素影響的控制同時會減弱對另一種影響因素的控制,難以適應大面積地形復雜的區(qū)域。本文提出一種組合濾波的方案,首先對初始LiDAR點云數據利用IPTD濾波算法去除了大部分的地物點,通過鄰域高差濾波的方法對IPTD濾波在陡坡、山頂等區(qū)域的缺失點進行了補全,最后通過平面擬合濾波的方法對數據進行了平滑,消除了近地點噪聲;經過融合、重采樣、內插等處理之后,制作了實驗區(qū) 1∶500地面DEM,通過同名點較差的方法進行精度評定,精度均滿足規(guī)范要求,證明組合濾波是一種有效可靠的濾波算法,并取得了如下結論:
(1)相比以往單一濾波方法,組合濾波方法的適應性更強,可以做到各濾波算法之間取長補短,有效降低了單一濾波算法的不足對結果的影響。
(2)組合濾波方法適用于大面積綜合場景下的LiDAR數據,對自然和人工陡坡、山頂等地形的濾波效果好;但是也存在處理步驟多、智能化程度低、對建筑物密集地區(qū)的過濾能力不足等問題。