羅伊萍，姜 挺，王 鑫，張 銳，羅 勝

(信息工程大學(xué)測(cè)繪學(xué)院，河南鄭州 450052)

基于數(shù)學(xué)形態(tài)學(xué)的LiDAR數(shù)據(jù)濾波新方法

羅伊萍，姜 挺，王 鑫，張 銳，羅 勝

(信息工程大學(xué)測(cè)繪學(xué)院，河南鄭州 450052)

數(shù)學(xué)形態(tài)學(xué)濾波是從激光雷達(dá)點(diǎn)云數(shù)據(jù)中識(shí)別地面點(diǎn)、創(chuàng)建數(shù)字高程模型的一種重要方法。在分析現(xiàn)有濾波方法的優(yōu)劣性以及數(shù)學(xué)形態(tài)學(xué)濾波方法存在的問題基礎(chǔ)上，提出一種新的具有一定自適應(yīng)性的數(shù)學(xué)形態(tài)學(xué)濾波算法。該方法通過分析LiDAR數(shù)據(jù)的特點(diǎn)，利用形態(tài)學(xué)算子提取原始點(diǎn)云數(shù)據(jù)的空白區(qū)域，進(jìn)行精確的重采樣生成DSM;然后將DSM多尺度形態(tài)學(xué)濾波的結(jié)果作為初始的DEM，分析由于局部地形突變導(dǎo)致可能存在誤分類區(qū)域;最后提取該區(qū)域邊緣外的地面點(diǎn)作為新的種子點(diǎn)，利用最小二乘平面擬合的方法進(jìn)行區(qū)域生長搜索誤分類地面點(diǎn)，從而提高分類精度。試驗(yàn)結(jié)果表明該方法能夠有效識(shí)別地面點(diǎn)和地物點(diǎn)，并且保留地形的細(xì)節(jié)信息。

激光雷達(dá);濾波;數(shù)學(xué)形態(tài)學(xué)濾波;數(shù)字高程模型

一、引 言

機(jī)載激光雷達(dá)(light deteation and ranging，LiDAR)技術(shù)自20世紀(jì)80年代引入攝影測(cè)量領(lǐng)域以來，受到廣大研究者的極大關(guān)注。LiDAR系統(tǒng)利用機(jī)載激光雷達(dá)測(cè)距系統(tǒng)和GPS/IMU能直接獲取地面點(diǎn)的三維坐標(biāo)，形成離散的、不規(guī)則的三維點(diǎn)云數(shù)據(jù)。為了從離散的點(diǎn)云數(shù)據(jù)中獲取數(shù)字高程模型(DEM)，必須對(duì)點(diǎn)的屬性進(jìn)行區(qū)分，將原始的點(diǎn)云數(shù)據(jù)分為地面點(diǎn)集和非地面點(diǎn)集，這個(gè)過程被稱為濾波。

目前，已有不少國內(nèi)外學(xué)者對(duì)激光雷達(dá)點(diǎn)云數(shù)據(jù)的濾波方法進(jìn)行了深入的研究，主要的濾波算法可分為以下三類:基于表面的內(nèi)插濾波算法［1-2］、基于區(qū)域的濾波算法［3-5］和基于約束曲面的濾波算法［6-7］。基于表面的內(nèi)插濾波算法的核心思想是通過一個(gè)較粗的起始DEM，逐步從備選數(shù)據(jù)點(diǎn)篩選并內(nèi)插加密DEM達(dá)到分類的目的。因此，濾波結(jié)果受到初始DEM影響較大并且誤差會(huì)隨著迭代的過程積累。而基于Snake樣條曲面等約束曲面的濾波方法認(rèn)為地面是一個(gè)連續(xù)且平緩變化的表面，這樣就過于強(qiáng)調(diào)地形的平緩變化而忽略了地形的復(fù)雜性［8］。

數(shù)學(xué)形態(tài)學(xué)濾波方法屬于基于區(qū)域的濾波算法，與基于表面的內(nèi)插濾波算法相反，它是一種自下而上、從局部出發(fā)擴(kuò)展到全局的濾波方法。借鑒柵格圖像的處理思想，根據(jù)激光點(diǎn)云數(shù)據(jù)生成的深度圖像，利用形態(tài)學(xué)開運(yùn)算剔除高于地面的點(diǎn)，通過逐步增大濾波窗口，得到逼近地形的一個(gè)表面。然而在濾波過程中，僅考慮窗口內(nèi)的地形特征，容易受到局部地形的影響，例如道路旁存在溝渠，當(dāng)濾波窗口尺寸大于道路寬度時(shí)會(huì)導(dǎo)致路面點(diǎn)被分為非地面點(diǎn)。另外，濾波過程中的參數(shù)需要用戶根據(jù)不同地形設(shè)置，自適應(yīng)性不強(qiáng)［3-4］。針對(duì)數(shù)學(xué)形態(tài)學(xué)濾波算法存在的問題，本文提出一種新的多尺度數(shù)學(xué)形態(tài)學(xué)濾波策略用于LiDAR數(shù)據(jù)濾波。

二、多尺度數(shù)學(xué)形態(tài)學(xué)濾波方法

數(shù)學(xué)形態(tài)學(xué)的語言和理論是基于集合運(yùn)算原理提取圖像中的特征。腐蝕和膨脹運(yùn)算是形態(tài)學(xué)圖像處理的基礎(chǔ)，通常用于“減少”(腐蝕)或“增大”(膨脹)圖像中特征形狀的尺寸?；诨叶葓D像的腐蝕和膨脹運(yùn)算是在結(jié)構(gòu)元素定義的鄰域內(nèi)選擇圖像像素值和結(jié)構(gòu)元素相作用后的最小或最大像素值。

在LiDAR點(diǎn)云數(shù)據(jù)生成的規(guī)則化數(shù)字表面模型(DSM)中，腐蝕和膨脹運(yùn)算可以定義為:

腐蝕

膨脹

其中，f為規(guī)則化的DSM;g為結(jié)構(gòu)元素;Z(i，j)為腐蝕或膨脹運(yùn)算后規(guī)則化的DSM中第i行第j列的高程值;w為結(jié)構(gòu)元素的窗口。

腐蝕和膨脹運(yùn)算組合后，形成開運(yùn)算和閉運(yùn)算用于LiDAR點(diǎn)云數(shù)據(jù)濾波，開運(yùn)算和閉運(yùn)算可以定義為:

開運(yùn)算

閉運(yùn)算

式中，開運(yùn)算通過先腐蝕后膨脹的方法，先將比結(jié)構(gòu)元素尺寸小的樹木點(diǎn)等非地面點(diǎn)從原始點(diǎn)云數(shù)據(jù)中移除掉，然后通過膨脹運(yùn)算恢復(fù)被腐蝕掉的建筑物邊緣形狀。然而，利用一個(gè)固定窗口大小的結(jié)構(gòu)元素很難移除各種尺寸的非地面物體，如果窗口過小，則大部分地面點(diǎn)被保留，只有小的非地面目標(biāo)，比如車或者樹木被移除，而城區(qū)的大型建筑物不能被移除。因此必須利用多尺度的濾波窗口進(jìn)行迭代計(jì)算，在迭代過程中逐步增大濾波窗口的尺寸。

每次迭代開運(yùn)算時(shí)，濾波的窗口尺寸wk按指數(shù)增長，wk=2bk+1，式中 k為迭代的次數(shù)，k=1，2，…，M;b為初始的窗口尺寸大小，一般為2。為了保證建筑物都被去除掉，最后一次迭代的窗口尺寸wk必須大于試驗(yàn)區(qū)域中最大建筑物的面積，一般迭代次數(shù)為5或6次。為了確保在濾波的過程中保留細(xì)節(jié)的地形特征，設(shè)置每次迭代高差閾值dhT，k

式中，wk為第k次濾波的窗口大小;c為DSM的格網(wǎng)間距;dh0和dhmax分別為最小和最大高差閾值;s為地形坡度參數(shù)。從式(5)可以看出，當(dāng)濾波窗口增大時(shí)，高差閾值dhT，k隨之增大，增幅大小由地形坡度s決定。如果DSM中某個(gè)格網(wǎng)點(diǎn)開運(yùn)算前后的高程值之差小于本次迭代的高差閾值dhT，k，則認(rèn)為是地面點(diǎn)，否則認(rèn)為是非地面點(diǎn)。

三、DEM提取過程

1.點(diǎn)云數(shù)據(jù)格網(wǎng)化

1)在點(diǎn)云數(shù)據(jù)的xy平面范圍內(nèi)，選擇合適的格網(wǎng)間距構(gòu)成m×n平面格網(wǎng)。對(duì)于每個(gè)離散的激光腳點(diǎn)，根據(jù)其平面坐標(biāo)將其分配到相應(yīng)的格網(wǎng)中。為了使盡可能多的激光點(diǎn)用于初始的DSM，一般選取比點(diǎn)間距稍小的格網(wǎng)間距值。由于點(diǎn)云數(shù)據(jù)分布的不均勻性，如圖1(a)所示，一個(gè)格網(wǎng)中可能會(huì)沒有或者有一個(gè)或多個(gè)激光點(diǎn)。如果一個(gè)格網(wǎng)內(nèi)有多個(gè)激光點(diǎn)落入，高程值越低的點(diǎn)作為地面點(diǎn)的概率越大，因此選取最低高程值作為格網(wǎng)的取值。記錄下每個(gè)格網(wǎng)中激光點(diǎn)的原始索引號(hào)，用于最后的精度評(píng)定。

2)在初始的DSM中，由于選取了格網(wǎng)內(nèi)所有激光腳點(diǎn)的最低值，因此由激光掃描系統(tǒng)誤差和隨機(jī)誤差引起的少量無意義低值粗差點(diǎn)也被保留下來。考慮到粗差點(diǎn)與相鄰點(diǎn)之間的高差形成脈沖波峰，在初始表面模型中的5×5鄰域內(nèi)，如果中心點(diǎn)高程值與24個(gè)鄰域點(diǎn)高程值的差值均大于高差閾值3 m，那么判定該點(diǎn)為噪聲點(diǎn)并且該格網(wǎng)中沒有激光點(diǎn)。

以m×n的二值標(biāo)記矩陣記錄下對(duì)應(yīng)格網(wǎng)中是否有激光點(diǎn)，矩陣值為0則沒有激光點(diǎn)，如圖1(c)所示。

2.DSM重采樣

DSM重采樣是針對(duì)格網(wǎng)化后的DSM中的空白格網(wǎng)進(jìn)行數(shù)據(jù)內(nèi)插。由于激光雷達(dá)的視場(chǎng)角比傳統(tǒng)攝影測(cè)量相機(jī)視場(chǎng)角小，飛行時(shí)容易造成掃描漏洞，如圖1(c)的頂部有細(xì)長的掃描漏洞。另外，水體對(duì)1 064～1 600 nm波段的激光有吸收效應(yīng)，也會(huì)造成數(shù)據(jù)空白區(qū)域，如圖1(c)的左側(cè)中部區(qū)域。如果存在高層建筑物，由于建造物的遮擋也會(huì)產(chǎn)生數(shù)據(jù)空白區(qū)域。由于這種數(shù)據(jù)的不均衡性，導(dǎo)致數(shù)據(jù)內(nèi)插的搜索半徑的大小很難確定，尤其是存在大量湖泊河流的區(qū)域，并且對(duì)建筑物遮擋的空白區(qū)域內(nèi)插生成的DSM擴(kuò)大了建筑物的面積。因此，要先針對(duì)這種大面積的數(shù)據(jù)空白區(qū)域進(jìn)行數(shù)據(jù)重采樣，然后以較小的搜索半徑對(duì)剩余的空白格網(wǎng)內(nèi)插，以提高重采樣的效率和精度。

首先，利用3×3的形態(tài)學(xué)窗口對(duì)二值標(biāo)記矩陣作閉運(yùn)算，去掉由于激光點(diǎn)分布不均勻造成的稀疏空白格網(wǎng)。其次，為了消除水域上的零散激光點(diǎn)，采用小面積消除法，求得標(biāo)記矩陣內(nèi)所有值為1的四連通區(qū)域，其個(gè)數(shù)小于40的連通區(qū)域認(rèn)為是由水體反射的激光點(diǎn)造成的，將它們賦值為0。然后，二值矩陣取反，再次使用小面積消除法消除零散的無激光點(diǎn)格網(wǎng)，得到圖1(d)。此時(shí)僅保留了面積大于40個(gè)格網(wǎng)數(shù)的水域和因飛行掃描漏洞造成的數(shù)據(jù)空白區(qū)域，為了求得這些區(qū)域的邊緣，將標(biāo)記矩陣用3×3的形態(tài)學(xué)窗口作膨脹運(yùn)算，再減去原標(biāo)記矩陣即得到如圖1(e)所示的邊緣。一般認(rèn)為水域是比較平緩的高程一致性區(qū)域，因此將水域邊緣的最低值作為水域的高程取值。掃描漏洞和建筑物遮擋區(qū)域也作相同處理。最后，由合適的搜索半徑對(duì)其余零散的空白格網(wǎng)進(jìn)行最鄰近值內(nèi)插，得到重采樣的DSM，如圖1(b)所示。

圖1 點(diǎn)云數(shù)據(jù)重采樣

3.多尺度數(shù)學(xué)形態(tài)學(xué)濾波

在點(diǎn)云數(shù)據(jù)規(guī)則化和重采樣后，本文采用上述的基于坡度和高差閾值的多尺度數(shù)學(xué)形態(tài)學(xué)濾波方法。針對(duì)城市區(qū)域地形起伏較小以及建筑物的最大面積一般不會(huì)超過120 m×120 m的實(shí)際情況，濾波的參數(shù)選擇如表1所示。

表1 濾波參數(shù)設(shè)置

4.質(zhì)量控制

由于數(shù)學(xué)形態(tài)學(xué)濾波僅考慮窗口內(nèi)的地形特征，在局部地形有突變的情況下，會(huì)導(dǎo)致相當(dāng)數(shù)目的地面點(diǎn)分類為非地面點(diǎn)，從而影響分類的精度。在這種情況下由于地面的連續(xù)性，誤分類區(qū)域在一定范圍內(nèi)與另外的地面區(qū)域相連接，通過提取相連接地面區(qū)域的邊緣點(diǎn)作為種子點(diǎn)，通過區(qū)域生長法搜索誤分類區(qū)域中的地面點(diǎn)，從而改進(jìn)分類精度。

將多尺度數(shù)學(xué)形態(tài)學(xué)濾波后的結(jié)果作為初始的DEM，并且從中提取非地面點(diǎn)集的二值化標(biāo)記圖像。用小面積消除法剔除標(biāo)記圖像中小面積的樹木點(diǎn)，再用3×3的窗口對(duì)標(biāo)記圖像進(jìn)行膨脹運(yùn)算再減去標(biāo)記圖像，得到大面積非地面點(diǎn)區(qū)域的外邊緣點(diǎn)，即相連接地面區(qū)域的邊緣點(diǎn)。從這些邊緣點(diǎn)中，根據(jù)它們3×3鄰域的最大和最小高程值之差，剔除建筑物的外邊緣點(diǎn)，保留有可能存在誤分類地面點(diǎn)集的邊緣點(diǎn)作為種子點(diǎn)。取種子點(diǎn)5×5鄰域中的地面點(diǎn)作最小二乘平面擬合，以擬合平面作為初始地面，搜索標(biāo)記圖像中與初始地面的高差小于高差閾值的點(diǎn)，作為誤分類的地面點(diǎn)。

四、試驗(yàn)與分析

1.試驗(yàn)結(jié)果

試驗(yàn)采用ISPRS網(wǎng)站提供的點(diǎn)云數(shù)據(jù)，這些點(diǎn)云數(shù)據(jù)是在2003年由ISPRS Commission III提供給廣大學(xué)者的，用于比較不同濾波算法的試驗(yàn)。試驗(yàn)選取城市區(qū)域的四塊測(cè)試數(shù)據(jù)Site1-Site4的末次回波激光點(diǎn)作為原始數(shù)據(jù)，點(diǎn)間距為1.5 m，重采樣的格網(wǎng)間距為1 m。試驗(yàn)區(qū)域散布著大、中、小型的不規(guī)則形狀的建筑物，同時(shí)有橋梁、隧道、鐵路、汽車以及與建筑物相鄰的樹木。在四塊測(cè)試數(shù)據(jù)中選取了九個(gè)樣本數(shù)據(jù)，對(duì)樣本數(shù)據(jù)進(jìn)行了手工分類，將激光腳點(diǎn)精確分類為地面點(diǎn)和非地面點(diǎn)兩類。利用分類后的樣本數(shù)據(jù)，可以對(duì)濾波算法進(jìn)行分類誤差的定量分析。

圖2(a)給出了應(yīng)用多尺度數(shù)學(xué)形態(tài)學(xué)算法對(duì)重采樣后Site2的DSM進(jìn)行濾波的結(jié)果，其中空白地區(qū)為分類后的非地面點(diǎn)。從圖2(a)可以看出，在A、C處由于地下通道與旁邊的路面存在高程突變，而在B處也由于地面不連貫變化，導(dǎo)致部分地面點(diǎn)明顯被誤分為非地面點(diǎn)。對(duì)濾波結(jié)果進(jìn)行質(zhì)量控制，搜索誤分類的地面點(diǎn)后得到最終分類結(jié)果，如圖2(c)所示。

圖3分別給出了Site1、Site3和Site4的分類后地面點(diǎn)重采樣生成的DEM。試驗(yàn)結(jié)果表明該算法能夠有效地識(shí)別地面和地物點(diǎn)，地面細(xì)節(jié)信息得到保留;樹木被有效識(shí)別，所有建筑物被濾除掉。

圖2 Site2試驗(yàn)結(jié)果

圖3 地面點(diǎn)重采樣的DEM效果圖

2.試驗(yàn)分析

文獻(xiàn)［9］對(duì)不同濾波方法從算法上進(jìn)行了試驗(yàn)比較，表2給出了本文濾波方法與文獻(xiàn)［9］和文獻(xiàn)［4］中所涉及算法的總誤差對(duì)比。

表2 樣本的總誤差對(duì)比

由表2可以看出，在陡坡地帶的樣本11、復(fù)雜建筑物密集的地形高程變化較大的樣本23以及鐵路等高頻率的地形起伏樣本42地區(qū)，Axelsson的內(nèi)插濾波方法的精度稍微高于本文和文獻(xiàn)［4］的形態(tài)學(xué)方法，因?yàn)閺拇值郊?xì)的處理方式能避免大地形起伏導(dǎo)致距離計(jì)算值較大而引起的誤判。而在地形平緩、建筑物密集的樣本12、樣本21、樣本31以及樣本41地區(qū)，本文方法明顯優(yōu)于內(nèi)插方法和文獻(xiàn)［4］的形態(tài)學(xué)方法，尤其是數(shù)據(jù)樣本21和樣本41，這是由于多尺度形態(tài)學(xué)濾波方法采取由局部出發(fā)擴(kuò)展到全局的策略，在平緩地區(qū)不容易導(dǎo)致誤判，尤其是在有數(shù)據(jù)空白和一個(gè)大型建筑物的樣本41，Axelsson的內(nèi)插濾波方法對(duì)數(shù)據(jù)缺失比較敏感，而文獻(xiàn)［4］的形態(tài)學(xué)方法以高程較低值填補(bǔ)數(shù)據(jù)空白后，則會(huì)將高于數(shù)據(jù)空白地區(qū)的周邊地面點(diǎn)分類為非地面點(diǎn)，同時(shí)該方法的濾波窗口最大值小于建筑物的面積，導(dǎo)致部分建筑物點(diǎn)被分類成地面點(diǎn)，而本文以形態(tài)學(xué)濾波的結(jié)果作為初始DEM，通過進(jìn)一步搜索誤分類的地面點(diǎn)明顯提高了分類的精度。

五、結(jié)束語

本文提出一種新的數(shù)學(xué)形態(tài)學(xué)濾波策略，該方法對(duì)點(diǎn)云數(shù)據(jù)進(jìn)行合理的重采樣，并用多尺度數(shù)學(xué)形態(tài)學(xué)濾波后作為初始的DEM，通過搜索誤分類地面點(diǎn)的方式提高分類精度。該方法解決了數(shù)學(xué)形態(tài)學(xué)濾波算法中存在的兩個(gè)問題:局部地形不連續(xù)對(duì)分類的影響以及坡度或高差閾值的人工選取。試驗(yàn)結(jié)果表明該方法能夠有效識(shí)別地面點(diǎn)和地物點(diǎn)，并且保留地形的細(xì)節(jié)信息。

［1］KRAUS K，PFEIFER N.Determination of Terrain Models in Wooded Areas with Airborne Laser Scanner Data［J］.ISPRS Journal of Photogrammetry and Remote Sensing，1998，53(4):193-203.

［2］PFEIFER N，STADLER P.Derivation of Digital Terrain Models in the SCOP++Environment［C］∥Proceedings of OEEPE Workshop on Airborne Laser Scanning and Interferometric SAR for Detailed Digital Elevation Models.Stockholm:［s.n.］，2001.

［3］ZHANG K，CHEN S C，WHITMAN D，et al.A Progressive Morphological Filter for Removing Nonground Measurement from LiDAR Data［J］.IEEE Transactions on Geoscience and Remote Sensing，2003，41(4):872-882.

［4］CHEN QGONG P，BALDOCCHI D，et al.Filtering Airborne Laser Scanning Data with Morphological Methods［J］.Photogrammetric Engineering and Remote Sensing，2007，73(2):175-185.

［5］張熠斌，隋立春，曲佳，等.基于數(shù)學(xué)形態(tài)學(xué)算法的機(jī)載 LiDAR點(diǎn)云數(shù)據(jù)快速濾波［J］.測(cè)繪通報(bào)，2009(5):16-18.

［6］ELMQVIST M JUNGERT E，LANTZ F，et al.Terrain Modelling and Analysis Using Laser Scanner Data［J］.International Archives of Photogrammetry and Remote Sensing，2001，34(3/W4):219-224.

［7］BROVELLI M A，CANNATA M.Digital Terrain Model Reconstruction in Urban Areas from Airborne Laser Scanning Data:the Method and an Example of the Town of Pavia(Northern Italy)［J］.Computers ＆ Geosciences，2004，30(4):325-331.

［8］黃先鋒，李卉，王瀟，等.機(jī)載LiDAR數(shù)據(jù)濾波方法評(píng)述［J］.測(cè)繪學(xué)報(bào)，2009，38(5):465-469.

［9］SITHOLE G，VOSSELMAN G.Experimental Comparison of Filter Algorithms for Bare-Earth Extraction from Airborne Laser Scanning Point Clouds［J］.ISPRS Journal of Photogrammetry and Remote Sensing， 2004，59(1-2):85-101.

［10］穆超，余潔，許磊，等.基于高分辨率遙感影像的DSM建筑物點(diǎn)的提取研究［J］.武漢大學(xué)學(xué)報(bào):信息科學(xué)版，2009，34(4):414-417.

A New Filtering Method for LiDAR Data Based on Mathematic Morphological Approach

LUO Yiping，JIANG Ting，WANG Xin，ZHANG Rui，LUO Sheng

0494-0911(2011)03-0015-05

P237

B

2010-02-23

羅伊萍(1982—)，女，湖南永州人，博士生，研究方向?yàn)長iDAR數(shù)據(jù)處理和航天遙感工程。