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

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

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

?

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

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

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

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

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

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

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

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

1算法原理及流程

1.1近似裸露地表的獲取

1.1.1水體邊界的提取及填充

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

(1)

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

(2)

圖1　水域邊界提取

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

1.1.2獲取近似裸露地表面

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

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)計(jì)算柵格濾波閾值et。從w中取出一個(gè)值

et=swic.

(4)

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

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

(5)

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

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

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

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

1.2點(diǎn)云濾波

1.2.12D空間的TPS變形模型

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

(6)

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

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

(7)

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

(8)

f(xi,yi)定義為

(9)

(10)

1.2.2濾波過(guò)程

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

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

(11)

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

(12)

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

(13)

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

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

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

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

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

圖2　My Method 與TerraScan的比較

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

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

3結(jié)論

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

圖3　S12和S52浮雕圖

參考文獻(xiàn)：

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

[2]肖春蕾,郭兆成,張宗貴,等.利用機(jī)載LiDAR數(shù)據(jù)提取與分析地裂縫[J].國(guó)土資源遙感,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]張小紅.機(jī)載激光雷達(dá)測(cè)量技術(shù)理論與方法[M].武漢:武漢大學(xué)出版社，2007:93-95.

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

[6]亢曉琛，劉紀(jì)平，林祥國(guó)．多核處理器的機(jī)載激光雷達(dá)點(diǎn)云并行三角網(wǎng)漸進(jìn)加密濾波方法[J]．測(cè)繪學(xué)報(bào),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]成曉倩，趙紅強(qiáng).基于區(qū)域生長(zhǎng)的LIDAR點(diǎn)云數(shù)據(jù)濾波[J].國(guó)土資源遙感，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點(diǎn)云數(shù)據(jù)分類方法研究[D].青島：國(guó)家海洋局第一研究所，2008:28-45.

[11] 隋立春，張熠斌，柳艷，等. 基于改進(jìn)的數(shù)學(xué)形態(tài)學(xué)算法的LiDAR點(diǎn)云數(shù)據(jù)濾波[J].測(cè)繪學(xué)報(bào)，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] 杜國(guó)軍，賈良文. 薄板樣條函數(shù)在空間數(shù)據(jù)插值中的應(yīng)用[J]. 計(jì)算機(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].測(cè)繪通報(bào)，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

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

中圖分類號(hào)：P237.3

文獻(xiàn)標(biāo)識(shí)碼：A

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