胡永杰，姚曉偉，胡倩偉,韓鳳納

(1.新疆兵團勘測設(shè)計院(集團)有限責(zé)任公司，新疆 烏魯木齊 830002；2.東華理工大學(xué) 測繪工程學(xué)院，江西 南昌 330013；3.北京四維遠見信息技術(shù)有限公司，北京 100070)

?

基于點云空洞修復(fù)和TPS模型的數(shù)學(xué)形態(tài)學(xué)機載LIDAR濾波

胡永杰1,2，姚曉偉1，胡倩偉3,韓鳳納1

(1.新疆兵團勘測設(shè)計院(集團)有限責(zé)任公司，新疆 烏魯木齊 830002；2.東華理工大學(xué) 測繪工程學(xué)院，江西 南昌 330013；3.北京四維遠見信息技術(shù)有限公司，北京 100070)

摘要：在分析現(xiàn)有濾波算法的基礎(chǔ)上，結(jié)合數(shù)據(jù)驅(qū)動和模型驅(qū)動式算子各自的優(yōu)點，提出基于點云空洞修復(fù)和TPS變形模型的數(shù)學(xué)形態(tài)學(xué)機載LIDAR點云濾波，該方法首先提取和修復(fù)由水域造成的大面積點云空洞，采用多尺度形態(tài)學(xué)開算子作用于修復(fù)的數(shù)據(jù)，得到近似裸露地表面；然后利用2D空間的TPS變形模型，以近似地表面為基礎(chǔ)，插值原始點云，根據(jù)插值與原始點云高程的差值大小去識別地面點和非地面點。通過定量分析，驗證該方法不僅有較高的濾波精度，而且也能較好的保留裸露地表的細節(jié)特征，同時該方法有助于輔助人工處理，提高數(shù)據(jù)處理的質(zhì)量。

關(guān)鍵詞：空洞修復(fù)；多尺度形態(tài)學(xué)開算子；2D空間TPS模型；機載LIDAR濾波

機載LIDAR技術(shù)是一種主動式遙感技術(shù)，它能夠快速獲取大面積區(qū)域內(nèi)高精度的空間三維地形信息，已經(jīng)在地球空間信息學(xué)科的許多領(lǐng)域得到廣泛應(yīng)用[1-2]。機載LIDAR獲取的原始數(shù)據(jù)，其結(jié)構(gòu)是離散的，以不規(guī)則方式分布在地形表面，因此稱為點云。以點云為基礎(chǔ)進行濾波分類，是該系統(tǒng)數(shù)據(jù)后處理重要的組成部分，也是進行后續(xù)高層次數(shù)據(jù)處理的基礎(chǔ)。而點云數(shù)據(jù)濾波及其配套的質(zhì)量控制是生成DEM最關(guān)鍵也是最耗時的一步，約占后處理流程60%～80%的時間,因此找到一種濾波精度較高的方法對提高DEM的質(zhì)量是至關(guān)重要的，同時近幾年DEM一個重要的應(yīng)用領(lǐng)域就是虛擬地理環(huán)境(Geographic Virtual Environments)，要提高人在虛擬環(huán)境的沉浸感，對地形起伏的仿真度提出更高要求[3]。這就要求濾波算法不僅要盡量減少分類誤差，同時也要盡可能保留地形特征的細節(jié)信息[4]。

目前濾波算法大致可以分為兩類：一類是數(shù)據(jù)驅(qū)動式的濾波算法，例如基于數(shù)學(xué)形態(tài)學(xué)的方法[5]、基于TIN結(jié)構(gòu)的分層方法[6]；另一類是模型驅(qū)動的濾波算法，例如基于分層穩(wěn)健線性估計的濾波算法[7]、基于樣條插值和區(qū)域增長的濾波算法[8]。這兩類算法大多都從單一的高差或坡度作為分類的依據(jù)，對地形情況復(fù)雜的場景很難獲得理想的結(jié)果，另一方面這兩類算法幾乎都沒有考慮到大面積的數(shù)據(jù)空洞對濾波保留地形特征的影響。

針對這兩類濾波算法所存在的這些不足，本文提出基于形態(tài)學(xué)理論和TPS變形模型的機載LIDAR點云濾波，該算法對點云的濾波不僅綜合考慮高差和坡度的影響，而且也減少大面積的空洞對濾波保留地形特征的影響。通過試驗證明該方法不僅有較高的濾波精度，而且也能較好地保留地形的細節(jié)特征信息,同時該方法也有助于輔助人工處理，提高數(shù)據(jù)處理的質(zhì)量。

1算法原理及流程

1.1近似裸露地表的獲取

1.1.1水體邊界的提取及填充

首先需要將離散的點云柵格化，根據(jù)三維坐標(biāo)的最大最小值，確定柵格圖的范圍，每個柵格值的大小取1 m范圍內(nèi)最近點的高程值。由于水體對激光的吸收，導(dǎo)致傳感器無法接收到回波信號，從而使部分區(qū)域點云過于稀疏,這可能會導(dǎo)致有一部分柵格無值。通常水體的高程在一個相鄰的區(qū)域內(nèi)高程最低，因此可以用水體邊界高程值最低的點填充水體空洞區(qū)域[9]。首先將柵格圖生成二值圖f(x)，這里引入形態(tài)學(xué)梯度算子檢測水體空洞區(qū)域的外邊界，算子如式(1)所示。

(1)

其中g(shù)稱為結(jié)構(gòu)窗口，窗口大小根據(jù)式(2)確定。

(2)

圖1　水域邊界提取

盲目的用高程最低值填充點云空洞會造成局部地形區(qū)域的不連續(xù)性，破壞地形特征[10]，因此水體造成的大面積空洞，采用水體邊界高程最低值的點填充,最終得到以點云高程值為基礎(chǔ)的柵格表面SF。

1.1.2獲取近似裸露地表面

離散點云形成柵格圖的時候，采用最近點填充柵格的方式，這必然會導(dǎo)致柵格圖像中有植被、房屋等非地面柵格單元，因此必須剔除這些非地面柵格，這里采用多尺度形態(tài)學(xué)開算子，該方法能夠較好的保留地形起伏的裸露地表特征[11]。其過程如下：

1)將柵格表面SF賦給SL，形成多尺度形態(tài)學(xué)窗口向量w={w1,w2,…,|wi≤wmax}，其中

wi=2bik+1.

(3)

式中:k=0.5,bi=i(i=0,1,2…)。

2)計算柵格濾波閾值et。從w中取出一個值

et=swic.

(4)

式中：s為測區(qū)的最大坡度估值；c為柵格大小。

3)對表面SL執(zhí)行形態(tài)學(xué)開算子得到表面ST。開算子定義為先對目標(biāo)圖像f通過結(jié)構(gòu)窗口g被腐蝕，然后在腐蝕的目標(biāo)圖像上進行膨脹操作，其數(shù)學(xué)過程為

(5)

其結(jié)構(gòu)窗口g大小為wi。

4)標(biāo)示地面柵格。當(dāng)ST-SL>et,在SL上對應(yīng)柵格標(biāo)示為地面柵格。

5)執(zhí)行SL=ST,重復(fù)2)到5)的步驟直到wi等于wmax為止。

通過以上的步驟，最終可以得到機載LIDAR掃描區(qū)的一個近似的表面S，該表面位于該區(qū)域DEM與DSM表面之間，因此可以作為點云濾波的參考表面。

1.2點云濾波

1.2.12D空間的TPS變形模型

薄板樣條函數(shù)TPS模型模擬一個薄金屬板在點約束下的彎曲變形，用一個簡單的代數(shù)式表示彎曲的能量變化，屬于非線性變換方法。假定有n觀測值z(xi,yi)，則觀測值可由式(6)表示[12]。

(6)

其中i=1,2,…,n，f(xi,yi)是觀測值z(xi,yi)的估值函數(shù)，ε(xi,yi)是估計誤差。

若使估計函數(shù)最大限度的接近觀測值，則需要使用最小化函數(shù)實現(xiàn)，算式為

(7)

其中平滑參數(shù)λ>0,本文中λ=1，其中If為估計函數(shù)的平滑項，定義為

(8)

f(xi,yi)定義為

(9)

(10)

1.2.2濾波過程

TPS相對于其它的插值方法具有反映高程異常變化的物理特性，具有光滑、連續(xù)、彈性好的特點，而且該方法不需要測量點呈規(guī)則分布[13],因此該方法較適合離散點云的插值，它能插值出光滑的表面，有利于提高插值精度[14]。濾波的過程如下：

2)求每個離散點Pi(xi,yi,zi)對應(yīng)地表面S上的近似坡度，計算濾波閾值hyi。坡度的算式為

(11)

其中，F(xiàn)x(i,j)，F(xiàn)y(i,j)為水平垂直梯度，水平梯度算式為

(12)

式中：i,j為第i個點對應(yīng)柵格S的行列號，需要注意的是：第1列上的水平梯度計算方法為第2列與第1列的差值，最后一列為最后兩列的差值。垂直梯度與水平梯度計算方法類似。

(13)

式中：hy為固定的一個閾值；sc為縮放因子(一般取值為：0.0到2.5)。

2試驗結(jié)果與分析

實驗數(shù)據(jù)是由國際攝影測量與遙感協(xié)會網(wǎng)(http://www.itc.nl/isprswgIII-3/filte-rtest)提供的8種場景LIDAR數(shù)據(jù)，專門用于測試濾波算法，其中城市和鄉(xiāng)村地區(qū)的數(shù)據(jù)各占 4 景;這些數(shù)據(jù)包含平原、植被、建筑物、公路、鐵路、橋梁、電線、水域、陡坡等不同的地形地貌特征。

本文算法用MATLAB實現(xiàn)，定量分析采用ISPRS提出的I類誤差、II類誤差和總誤差作為算法精度的評定準(zhǔn)則，參數(shù)和實驗結(jié)果如表1所示。

Terrasolid開發(fā)人員將上述實驗數(shù)據(jù)在Terracan中半自動化地進行了處理，濾波結(jié)果達到該處理的最佳狀態(tài)。這里將該軟件的處理結(jié)果與本方法進行比較，比較結(jié)果如圖2所示。

圖2　My Method 與TerraScan的比較

I類誤差和II類誤差反映算法的適應(yīng)性，總誤差反眏算法的可行性[15],圖2中該方法與Terracan軟件中的算法濾波結(jié)果比較，一類誤差和總誤差明顯優(yōu)于Terracan軟件的濾波，這反映出該算法能夠較好的保留裸露地表上的點，從機載LIDAR點云中能夠提取質(zhì)量較高的DEM；總誤差較低能總體反映算法在實際生產(chǎn)應(yīng)用中具有較高的應(yīng)用價值和優(yōu)越性，II類誤差雖然較Terracan軟件的濾波結(jié)果高，但是平均誤差也能控制在7.6%以下，這也表明該方法能夠適應(yīng)不同復(fù)雜程度的地形環(huán)境。表1中該方法II類誤差較I類誤差稍大，這一特點與當(dāng)下大部分濾波算法剛好相反(即I類誤差大于II類誤差)，經(jīng)驗表明，與Ⅰ類誤差相比，Ⅱ類誤差的點云大部分有著高于地面點云的高程，所以在可視化環(huán)境下，人工干預(yù)時能更容易地識別和修復(fù)Ⅱ類誤差帶來的分類錯誤，這能較好地提高數(shù)據(jù)處理效率和數(shù)據(jù)質(zhì)量。因此該算法無論是自動化處理還是輔助人工處理都能較好提高濾波精度和數(shù)據(jù)質(zhì)量。

如圖3所示，S12和S52分別為城市和山區(qū)的浮雕圖，圖3(b)為人工逐點分類的結(jié)果圖，可認為是真實的裸露地表。對比圖3(c)、(b)和(a)可以明顯地看到，本文方法無論在城區(qū)還是山區(qū)都能夠較好地保持原始裸露地貌，而且對于地形起伏較大的山區(qū)也能保留地形起伏的輪廓細節(jié)特征；圖中紅色箭頭所指示的陡坎、斜坡區(qū)域內(nèi)的植被幾乎都被剔除，較好地保持了原始的地形地貌，沒有出現(xiàn)過度剔除地形特征的現(xiàn)象。這種情況在虛擬現(xiàn)實中能夠改善場景的真實性，為客戶提供較滿意的虛擬沉浸感。

3結(jié)論

本文結(jié)合數(shù)據(jù)驅(qū)動和模型驅(qū)動式算子各自的優(yōu)點，提出基于點云空洞修復(fù)和TPS變形模型的數(shù)學(xué)形態(tài)學(xué)機載LIDAR點云濾波，通過試驗驗證該方法相對于商用軟件Terracan，濾波精度較高。I類誤差平均為2.19%，II類誤差平均為7.54%，總誤差平均為3.06%，因此該算法不僅提高濾波精度，而且該算法能夠適應(yīng)不同復(fù)雜程度的地形環(huán)境，適應(yīng)性強;I類誤差小于II類誤差，說明該方法能較好地保存原始裸露地表上的點，用該方法提取的DEM用于三維可視化和虛擬現(xiàn)實中，可以明顯提高場景真實感和沉浸感，所以在可視化環(huán)境下，人工干預(yù)時能更容易地識別和修復(fù)Ⅱ類誤差帶來的分類錯誤，這有助于輔助人工處理提高工作效率，提升數(shù)據(jù)處理的質(zhì)量。當(dāng)然該方法有些不足，比如在試驗區(qū)S53中二類誤差較大，S53為人為改造出現(xiàn)的間斷地形，因此該方法應(yīng)進一步改進，以降低在間斷地形上的二類誤差。

圖3　S12和S52浮雕圖

參考文獻：

[1]吳芳,張宗貴,郭兆成,等.基于機載LiDAR點云濾波的礦區(qū)DEM 構(gòu)建方法[J].國土資源遙感,2015,27(1):62-67.

[2]肖春蕾,郭兆成,張宗貴,等.利用機載LiDAR數(shù)據(jù)提取與分析地裂縫[J].國土資源遙感,2014,26(4):111-118.

[3]PINGEL T,CLARKE K C.An improved simple morphological filter for the terrain classification of airborne LIDAR data[J].ISPRS Journal of Photogrammetry and Remote Sensing,2013,77:21-30.

[4]張小紅.機載激光雷達測量技術(shù)理論與方法[M].武漢:武漢大學(xué)出版社，2007:93-95.

[5]沈晶,劉紀(jì)平,林祥國. 用形態(tài)學(xué)重建進行機載激光雷達數(shù)據(jù)濾波的方法研究 [J]. 武漢大學(xué)學(xué)報(信息科學(xué)版), 2011, 36(2)：167-170.

[6]亢曉琛，劉紀(jì)平，林祥國．多核處理器的機載激光雷達點云并行三角網(wǎng)漸進加密濾波方法[J]．測繪學(xué)報,2013,42(3):331-336.

[7]PFEIFER N, STADLER P, BRIESE C. Derivation of Digital Terrain Models in the SCOP Environment [C]. OEEPE Workshop on Airborne Laser Scanning and Interferometric SAR for Detailed Digital Elevation Models. Stockholm:Sweden, 2001:1-12.

[8]成曉倩，趙紅強.基于區(qū)域生長的LIDAR點云數(shù)據(jù)濾波[J].國土資源遙感，2008，20(4)：6-8.

[9]QI Chen,PENG Gong,DENNLS B,et al.Filtering Airborne Laser Scanning Data With Morphological Methods[J].Photogrammetric Engineering & Remote Sensing,2007(73):175-185.

[10] 劉文.基于數(shù)學(xué)形態(tài)學(xué)原理和TerraScan的Lidar點云數(shù)據(jù)分類方法研究[D].青島：國家海洋局第一研究所，2008:28-45.

[11] 隋立春，張熠斌，柳艷，等. 基于改進的數(shù)學(xué)形態(tài)學(xué)算法的LiDAR點云數(shù)據(jù)濾波[J].測繪學(xué)報，2010,39(4)：390-396.

[12] BOER E P J, MDEBEURS K,HARTKAMP A D.Kriging and thin plate splines for mapping climate variables[J].International Journal of Applied Earth Observation and Geoinformation,2001,3(2):146-154.

[13] 杜國軍，賈良文. 薄板樣條函數(shù)在空間數(shù)據(jù)插值中的應(yīng)用[J]. 計算機工程與應(yīng)用,2009,45(36):238-240.

[14] MONGUS D, ZLIK B. Parameter-free ground filtering of LiDAR data for automatic DTM generation[J]. ISPRS Journal of Photogrammetry and Remote Sensing,2011,67:1-12.

[15] 吳叢叢,盧小平. 基于TIN的LiDAR數(shù)據(jù)濾波算法研究[J].測繪通報，2013(3):32-34.

[責(zé)任編輯：張德福]

Filtering of airborne LiDAR based on mathematical morphological of the repairing point cloud hole and TPS model

HU Yongjie1,2,YAO Xiaowei1,HU Qianwei3,HAN Fengna1

(1.XPCC Surveying & Designing Institute(Group) Co.,Ltd., Urumqi 830002,China;2.Institute of Surveying and Mapping, East China Institute of Technology, Nanchang 330013,China;3.Beijing Geo-VisionTech.Co.,Ltd., Beijing 100070,China)

Abstract：Based on the basic analysis of existing filtering algorithm,combiningrespective advantages of data-driven and model-driven,a filtering method is proposed based on mathematical morphological of the repairing point cloud hole and TPS model.This method first extracts and repaires the point cloud,using multi-scale morphological opening operator acting on the restoration of data to form an approximate bare ground surface.Then by the TPS deformation model of 2D space based on approximate surface,interpolation of the original point cloud,the difference between the original point cloud elevation interpolation,ground point and non-ground points are identified.In this paper,through quantitative analysis shows that the method not only has high filtering precision but also better retains the minutiae of bare surface.Meanwhile,the method can also help auxiliary artificial processing and improve the quality of data processing.

Key words：repairing holes;multi-scale morphological opening operator;2D space TPS model;airborne LiDAR filtering

DOI:10.19349/j.cnki.issn1006-7949.2016.07.006

收稿日期：2015-08-20

作者簡介：胡永杰(1986-)，男，助理工程師.

中圖分類號：P237.3

文獻標(biāo)識碼：A

文章編號：1006-7949(2016)07-0028-05