基于改進(jìn)坡度的自適應(yīng)數(shù)學(xué)形態(tài)學(xué)點(diǎn)云濾波

陳斐然，李 浩，李語(yǔ)旻，陳仁喜

（1.河海大學(xué) 地球科學(xué)與工程學(xué)院，江蘇 南京 210098）

基于改進(jìn)坡度的自適應(yīng)數(shù)學(xué)形態(tài)學(xué)點(diǎn)云濾波

陳斐然1，李 浩1，李語(yǔ)旻1，陳仁喜1

（1.河海大學(xué) 地球科學(xué)與工程學(xué)院，江蘇 南京 210098）

首先介紹了傳統(tǒng)形態(tài)學(xué)算法，并分析了其應(yīng)用方面的不足；然后提出了一種基于改進(jìn)坡度參數(shù)的自適應(yīng)形態(tài)學(xué)濾波方法；最后利用實(shí)際點(diǎn)云數(shù)據(jù)進(jìn)行了實(shí)驗(yàn)。該方法通過劃分格網(wǎng)自動(dòng)計(jì)算當(dāng)前區(qū)域坡度的平均值，從而得到適當(dāng)?shù)钠露葏?shù)用于后續(xù)高度閾值的計(jì)算。實(shí)驗(yàn)結(jié)果表明，該方法能更有效地識(shí)別地面點(diǎn)和地物點(diǎn)，對(duì)被建筑物包圍的植被點(diǎn)濾除效果明顯。

LiDAR點(diǎn)云；數(shù)學(xué)形態(tài)學(xué)；濾波；分類

LiDAR不僅能快速獲取大面積點(diǎn)云數(shù)據(jù)，而且能從中得到高精確DEM，因此受到各測(cè)繪行業(yè)的高度重視。點(diǎn)云濾波是LiDAR數(shù)據(jù)處理的重要內(nèi)容之一，也是最先需要解決的問題。濾波分類精度的高低直接影響點(diǎn)云數(shù)據(jù)的使用質(zhì)量。點(diǎn)云濾波即對(duì)在地形表面隨機(jī)分布的離散點(diǎn)進(jìn)行識(shí)別與處理，區(qū)分落在地面上的點(diǎn)（地面點(diǎn)）和落在地面物體上的點(diǎn)（植被﹑建筑物﹑車輛等地物點(diǎn)），分別形成地面點(diǎn)與地物點(diǎn)數(shù)據(jù)集合。

目前國(guó)內(nèi)外已有許多專家學(xué)者對(duì)點(diǎn)云濾波方法進(jìn)行了詳細(xì)研究[1-2]。張子林[3]等使用基于TIN的濾波算法對(duì)高程突變地物進(jìn)行了濾波，但對(duì)低矮的地面物體進(jìn)行濾波時(shí)，誤差會(huì)顯著增大。參考文獻(xiàn)[4]～[6]提出一種基于坡度變化的濾波算法，但該算法需提前知道地形坡度值，并確定窗口大小，區(qū)域內(nèi)每個(gè)點(diǎn)都需進(jìn)行坡度計(jì)算，并同其他點(diǎn)進(jìn)行比較，以確定該點(diǎn)是否為地面點(diǎn)。同時(shí)，如果僅根據(jù)坡度設(shè)置統(tǒng)一的閾值，則會(huì)濾除掉一些真實(shí)的地形細(xì)節(jié)特征，造成分類誤差，所以必須根據(jù)作業(yè)區(qū)域的實(shí)際地形狀況設(shè)置不同的過濾參數(shù)值。另外一種常用的濾波算法是基于數(shù)學(xué)形態(tài)學(xué)原理的濾波算法[7-11]。數(shù)學(xué)形態(tài)學(xué)算法不僅可用于點(diǎn)云數(shù)據(jù)濾波，也可用于建筑物或道路點(diǎn)云數(shù)據(jù)的提取。然而現(xiàn)有的基于經(jīng)典數(shù)學(xué)形態(tài)學(xué)算法的LiDAR點(diǎn)云濾波方法存在以下不足：①無法根據(jù)地形特點(diǎn)動(dòng)態(tài)調(diào)整窗口。利用一個(gè)固定窗口大小的結(jié)構(gòu)元素很難過濾各種尺寸的地物點(diǎn)，若窗口過小則只有較小的地物點(diǎn)被去除，如汽車或植被，但此時(shí)窗口無法完全包住大型建筑物，使得大型建筑物不能被有效濾除；而較大的濾波窗口對(duì)于含有山坡和低矮山包的地區(qū)而言，很容易將山包和山坡當(dāng)成非地面點(diǎn)而錯(cuò)誤濾除。②無法顧及窗口內(nèi)地形的細(xì)節(jié)特征。使用一個(gè)固定的“高度閾值”作為判斷是否為地面點(diǎn)的唯一標(biāo)準(zhǔn)顯然是不合理的。③坡度參數(shù)沒有考慮實(shí)際地形情況。坡度參數(shù)的取值由用戶自己輸入，并沒有真正考慮地形起伏的實(shí)際情況，對(duì)后期濾波的準(zhǔn)確程度產(chǎn)生影響。

本文對(duì)傳統(tǒng)數(shù)學(xué)形態(tài)學(xué)濾波算法進(jìn)行了改進(jìn)，提出了一種基于改進(jìn)坡度參數(shù)的自適應(yīng)數(shù)學(xué)形態(tài)學(xué)濾波算法，不僅能自動(dòng)調(diào)整窗口大小，而且可以通過劃分格網(wǎng)來自動(dòng)計(jì)算當(dāng)前區(qū)域坡度的平均值，得到適當(dāng)?shù)钠露葏?shù)用于后續(xù)高度閾值的計(jì)算，使得對(duì)LiDAR點(diǎn)云的識(shí)別更加準(zhǔn)確。

1 數(shù)學(xué)形態(tài)學(xué)濾波方法

1.1 經(jīng)典數(shù)學(xué)形態(tài)學(xué)濾波

數(shù)學(xué)形態(tài)學(xué)是基于集合運(yùn)算原理進(jìn)行特征提取的，其中集合表示圖像中不同的對(duì)象。腐蝕與膨脹分別用于減少或增大圖像中特征形狀的尺寸。腐蝕運(yùn)算即以窗口內(nèi)高程的最小值代替窗口中心的高程值，而膨脹運(yùn)算即以窗口內(nèi)高程的最大值代替窗口中心的高程值。

在采用數(shù)學(xué)形態(tài)學(xué)方法對(duì)LiDAR點(diǎn)云進(jìn)行濾波時(shí)，一般使用這兩個(gè)基礎(chǔ)算子的組合，主要是開運(yùn)算和閉運(yùn)算，其表達(dá)式為：

式中，Θ為腐蝕運(yùn)算符號(hào)；⊕為膨脹運(yùn)算符號(hào)；○為開運(yùn)算符號(hào)；●為閉運(yùn)算符號(hào)。

開運(yùn)算是先腐蝕后膨脹的運(yùn)算過程，可先將比結(jié)構(gòu)尺寸小的地物點(diǎn)去除，再膨脹恢復(fù)被腐蝕的建筑物邊緣形狀。這種去除比結(jié)構(gòu)元素窗口尺寸小的物體同時(shí)保留比結(jié)構(gòu)元素窗口大的物體特性，常被用于提取LiDAR點(diǎn)云數(shù)據(jù)的地面點(diǎn)，因此在LiDAR點(diǎn)云濾波中得到廣泛應(yīng)用。

1.2 漸進(jìn)式數(shù)學(xué)形態(tài)學(xué)濾波

針對(duì)單一窗口和單一閾值數(shù)學(xué)形態(tài)學(xué)濾波方法存在的問題，ZHANG K[7]等提出的漸進(jìn)式數(shù)學(xué)形態(tài)學(xué)濾波方法可從LiDAR點(diǎn)云中濾除大部分的地物點(diǎn)（見圖1）。

圖1 漸進(jìn)式數(shù)學(xué)形態(tài)學(xué)濾波方法

由圖1可知，在第一次濾波處理后，比窗口尺寸L1小的樹木點(diǎn)會(huì)被濾除，而大型建筑物則被保留；同樣尺寸小于L1的地面特征被削平，山坡最高點(diǎn)的高程與第一次濾波后的高程差值為dhmax(t),1，而山坡上隨機(jī)點(diǎn)P的高程與第一濾波后的高程差值為dhp,1。在隨后的迭代過程中，窗口大小增加到L2，此時(shí)大型建筑物尺寸小于當(dāng)前窗口尺寸，大型建筑物點(diǎn)被濾除，而地形特征點(diǎn)前后兩次濾波的高程差值為dhmax(t),2。

引入坡度參數(shù)和高度閾值可避免濾波后起伏地形的高點(diǎn)被濾除，但無法保留地形的細(xì)節(jié)特征，將地形坡度參數(shù)s和每次迭代高程差閾值dhT,k分別定義為：

式中，dhmax(t),k為第k次濾波后與前一次濾波的最大高程差；wk為第k次濾波的窗口大小；c為格網(wǎng)間距；dh0為初始高差閾值；dhmax為最大高差閾值。如果某個(gè)點(diǎn)前后兩次濾波的高程差值小于閾值，則分類成地面點(diǎn)，否則為地物點(diǎn)。

對(duì)于城市區(qū)域的LiDAR數(shù)據(jù)而言，若dhT,k取值適當(dāng)，則可確保比s小的地形特征表面被保留下來。城市中主要的非地面點(diǎn)對(duì)象為建筑物﹑汽車﹑植被等，通過逐漸增大的濾波窗口重復(fù)使用開運(yùn)算，汽車與植被點(diǎn)在首次運(yùn)算時(shí)就會(huì)被濾除，而高大的房屋會(huì)在最后幾次運(yùn)算中被過濾掉。一般情況下，將最低房屋的高度設(shè)為最大高差閾值dhmax，以保證復(fù)雜建筑物被完全識(shí)別濾除。

2 自適應(yīng)數(shù)學(xué)形態(tài)學(xué)濾波方法

[7]中，s的取值是由用戶根據(jù)地形實(shí)際情況自己輸入的，城市與山區(qū)地形的s取值有較大差別，且城市內(nèi)部中的s也是不同的。本文采用首先規(guī)則劃分格網(wǎng)并進(jìn)行高程賦值，然后利用坡度公式求解當(dāng)前區(qū)域各格網(wǎng)的坡度值，最終取其平均值作為區(qū)域近似坡度值的改進(jìn)坡度參數(shù)形態(tài)學(xué)濾波方法，以提高算法在不同地形環(huán)境下的自適應(yīng)能力和完善性。改進(jìn)算法流程如圖2所示。

圖2 改進(jìn)算法流程圖

格網(wǎng)尺寸過大會(huì)降低LiDAR點(diǎn)云高程的精度，而網(wǎng)格尺寸過小則會(huì)造成數(shù)據(jù)計(jì)算量過大﹑運(yùn)行緩慢。將區(qū)域總面積除以區(qū)域內(nèi)點(diǎn)云總數(shù)得到的平均點(diǎn)間距作為規(guī)則化點(diǎn)云格網(wǎng)的大小c，這樣做的目的是使得大部分格網(wǎng)中只包含1個(gè)LiDAR點(diǎn)，滿足一個(gè)激光點(diǎn)對(duì)應(yīng)一個(gè)格網(wǎng)的要求。在此基礎(chǔ)上，由已知區(qū)域內(nèi)最大和最小的x﹑y值與某一個(gè)格網(wǎng)中激光點(diǎn)坐標(biāo)計(jì)算得到該格網(wǎng)的行列號(hào)，并依次遍歷全部點(diǎn)云數(shù)據(jù)。

格網(wǎng)規(guī)則劃分完成后，則需對(duì)每個(gè)格網(wǎng)進(jìn)行高程賦值。這里分為兩種情況，一種是當(dāng)有多個(gè)點(diǎn)落在同一個(gè)格網(wǎng)中時(shí)，則取其中最低點(diǎn)的高程作為該格網(wǎng)的高程；另一種是當(dāng)格網(wǎng)中沒有任何激光點(diǎn)時(shí)，則將最鄰近點(diǎn)的高程值賦給該格網(wǎng)。所有格網(wǎng)賦值完成后，再利用新生成的規(guī)則格網(wǎng)高程數(shù)據(jù)建立初始LiDAR點(diǎn)云的數(shù)字表面模型（DSM）用于下一步坡度計(jì)算。

對(duì)于初始DSM，為了自動(dòng)計(jì)算地表面的坡度值，本文使用二階坡度計(jì)算公式來計(jì)算每個(gè)格網(wǎng)所對(duì)應(yīng)的坡度值：

式中，Zi,j+1﹑Zi,j-1﹑Zi+1,j﹑Zi-1,j分別對(duì)應(yīng)相應(yīng)格網(wǎng)行列號(hào)的格網(wǎng)高程值；c為格網(wǎng)大小。

將式（5）﹑（6）代入式（7）中就可計(jì)算得到各格網(wǎng)的坡度S(i, j)，但該坡度并不是所需要的實(shí)際地形坡度，而是地表坡度。在得到的地表坡度中，除了一小部分地面坡度外，其余絕大多數(shù)是植被﹑建筑物等地物點(diǎn)產(chǎn)生的坡度。為了能獲得實(shí)際地形的近似坡度值，必須進(jìn)行預(yù)處理。本文首先利用合適窗口大小的開算子遍歷全部地表坡度值，以剔除建筑物和部分植被對(duì)坡度的影響。由于坡度平均值代表區(qū)域內(nèi)地形的普遍情形，為了進(jìn)一步消除建筑物周圍植被對(duì)坡度的干擾，再對(duì)當(dāng)前區(qū)域的全部格網(wǎng)坡度值進(jìn)行相應(yīng)運(yùn)算，取其平均值作為該區(qū)域的近似坡度，并用于后續(xù)高度閾值的計(jì)算，最大限度地保留實(shí)驗(yàn)區(qū)域的地形細(xì)節(jié)特征。

3 實(shí)驗(yàn)結(jié)果與分析

3.1 實(shí)驗(yàn)數(shù)據(jù)與結(jié)果

實(shí)驗(yàn)數(shù)據(jù)采用ISPRS網(wǎng)站提供的城市點(diǎn)云數(shù)據(jù)，主要是由大小不一的建筑物及不同高度的植被所組成。實(shí)驗(yàn)區(qū)域大小為1 450 m×505 m，點(diǎn)云高差為50 m，點(diǎn)云間距為1.5 m。原始點(diǎn)云如圖3所示，可以看出實(shí)驗(yàn)區(qū)分布著大﹑中﹑小不規(guī)則形狀的建筑物，同時(shí)有道路﹑汽車﹑樹木等地物。圖4為采用基于改進(jìn)坡度參數(shù)的自適應(yīng)數(shù)字形態(tài)學(xué)濾波方法濾波后的結(jié)果。

在原始點(diǎn)云中選取一處復(fù)雜城市環(huán)境作為樣本，共有26 984個(gè)離散激光點(diǎn)，圖5b為采用漸進(jìn)式數(shù)學(xué)形態(tài)學(xué)濾波算法的結(jié)果，圖5c為采用本文算法的結(jié)果。

圖3 原始數(shù)據(jù)

圖4 改進(jìn)算法濾波結(jié)果

圖5 復(fù)雜區(qū)域?yàn)V波結(jié)果

3.2 實(shí)驗(yàn)結(jié)果分析

實(shí)驗(yàn)結(jié)果表明，基于改進(jìn)坡度參數(shù)的數(shù)學(xué)形態(tài)學(xué)濾波方法能更有效地區(qū)分地面點(diǎn)與地物點(diǎn)。在圖4中，絕大多數(shù)植被與建筑物都被濾除掉，同時(shí)大部分地面細(xì)節(jié)特征得以保留下來。從圖5b可以看出，漸進(jìn)式數(shù)學(xué)形態(tài)學(xué)濾波對(duì)被建筑物包圍的植被點(diǎn)濾除效果不佳，而本文提出的改進(jìn)算法能有效去除建筑物周圍的植被點(diǎn)，保留地面真實(shí)信息。

為了能更好地定量分析實(shí)驗(yàn)結(jié)果，將樣本區(qū)域的人工判讀分類結(jié)果作為準(zhǔn)確分類結(jié)果（地面點(diǎn)14 023 個(gè)，非地面點(diǎn)12 961個(gè)），統(tǒng)計(jì)參考文獻(xiàn)[7]的算法和本文提出的改進(jìn)算法的誤差。第Ⅰ類誤差﹑第Ⅱ類誤差和總誤差的定義如表1所示。表2為在樣本區(qū)域?yàn)V波結(jié)果中兩種算法的誤差對(duì)比。

表1 濾波誤差的定義

表2 濾波結(jié)果定量分析

從表2可以看出，改進(jìn)算法的第Ⅰ類誤差較參考文獻(xiàn)[7]算法有明顯降低，說明地面點(diǎn)識(shí)別分類更加精確；而第Ⅱ類誤差沒有太大的變化，總誤差則明顯低于參考文獻(xiàn)[7]算法。

4 結(jié) 語(yǔ)

本文所提出的基于改進(jìn)坡度參數(shù)的自適應(yīng)數(shù)學(xué)形態(tài)學(xué)濾波方法對(duì)現(xiàn)有的數(shù)學(xué)形態(tài)學(xué)濾波算法進(jìn)行了改進(jìn)，采用了自動(dòng)劃分格網(wǎng)求解坡度平均值作為坡度參數(shù)，并用于后續(xù)高度閾值計(jì)算的方法，有效解決了被建筑物包圍的植被點(diǎn)濾除問題。實(shí)驗(yàn)結(jié)果顯示，該方法能有效區(qū)分地面點(diǎn)與非地面點(diǎn)，并進(jìn)一步保留了地形的細(xì)節(jié)特征，降低了濾波誤差。

參考文獻(xiàn)

[1] 王金亮,陳聯(lián)君.激光雷達(dá)點(diǎn)云數(shù)據(jù)的濾波算法述評(píng)[J].遙感技術(shù)與應(yīng)用,2010,25(5):632-638

[2] ZHANG K, Whitman D. Comparison of Three Algorithms for Filtering Airbone LiDAR Data[J].Photogrammetric Engineering & Remote Sensing,2005,71(3):313-324

[3] 張子林,門林杰.基于TIN改進(jìn)的LiDAR濾波算法 [J].工程地球物理學(xué)報(bào),2012,9(3):362-364

[4] 許可,王運(yùn)巧,胡少興.基于坡度與迭代運(yùn)算的LiDAR點(diǎn)云數(shù)據(jù)濾波[J].機(jī)械工程與自動(dòng)化,2016(2):71-73

[5] 張寧寧,楊英寶,于雙.基于坡度和區(qū)域生長(zhǎng)的城市LiDAR點(diǎn)云濾波方法[J].地理空間信息,2016,14(3):30-32

[6] Sithole G. Filtering of Laser Altimetry Data Using a Slope Adaptive Filter[J]. International Archives of Photogrammetry and Remote Sensing,2001,34(3/W4):22-24

[7] 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

[8] 隋立春,張熠斌,柳艷,等.基于改進(jìn)的數(shù)學(xué)形態(tài)學(xué)算法的LiDAR點(diǎn)云數(shù)據(jù)濾波[J].測(cè)繪學(xué)報(bào),2010,39(4):390-396

[9] CHEN Q, GONG P , Baldocchi D ,et al. Filtering Airborne Laser Scanning Data with Morphological Methods[J].Photogrammetry Engineering and Remote Sensing,2007,73(2):175-185

[10] 羅伊萍,姜挺,龔志輝,等.基于自適應(yīng)和多尺度數(shù)學(xué)形態(tài)學(xué)的點(diǎn)云數(shù)據(jù)濾波方法[J].測(cè)繪科學(xué)技術(shù)學(xué)報(bào),2009,26(6):426-429

[11] 李峰,崔希民,袁德寶,等.改進(jìn)坡度的LiDAR點(diǎn)云形態(tài)學(xué)濾波算法[J].大地測(cè)量與地球動(dòng)力學(xué),2012,32(5):128-132

P228

B

1672-4623（2017）09-0028-04

10.3969/j.issn.1672-4623.2017.09.010

2016-09-18。

項(xiàng)目來源：國(guó)家自然科學(xué)基金資助項(xiàng)目（41471276）。

陳斐然，碩士研究生，研究方向?yàn)長(zhǎng)iDAR數(shù)據(jù)處理與地圖制圖。