基于形態(tài)學的機載LiDAR點云濾波研究

戴朦夢，趙俊三*，裴 旭(云南昆明理工大學國土資源工程學院，云南昆明650093)

機載LiDAR點云數(shù)據(jù)是當前獲得高精度DEM的重要手段之一，但激光雷達數(shù)據(jù)是分布在物體表面的點云，具有離散，不均勻分布，數(shù)據(jù)類型多等特點［1］，這種特點使得數(shù)據(jù)采集具有一定的盲目性，這種盲目性體現(xiàn)在反射激光信號數(shù)據(jù)可能來自裸露地表、建筑物和植被，還有可能來自天空飛行的物體以及由于其他誤差導致的極低無意義數(shù)據(jù)。因此要得到高精度的DEM，就必須要把裸露地表上的點和其他點分離開來，這個過程就叫濾波。

形態(tài)學濾波算法作為一種自下而上的基于局部高差突變思想的濾波器，其理論成熟，操作簡單，效率高［2］，因而其應用廣泛，但其在機載激光雷達點云濾波中的應用還處于不斷完善的階段。1996年Kilian用形態(tài)學獲取了DEM［3］，隨后Vosselman明確提出了“形態(tài)學濾波器”的概念［4］，此后的方法大都是基于該算法的擴展與改進的。Zhang提出了一種漸進多尺度形態(tài)學濾波，通過逐漸增大開運算的窗口尺寸和高程閾值的方法進行濾波，但是濾波尺寸和高程閾值只有對實驗區(qū)的地形和特征進行分析才能得到，在現(xiàn)實中卻很難做到，只有反復試驗才能得到［5］。Chen提出了機載激光雷達數(shù)據(jù)形態(tài)學濾波算法，該算法在漸進多尺度上，增加了一個條件，判斷最大地面物體邊緣高程是漸變還是突變，如果是漸變則保留，如果是高程差大于一定的值，則作為非地面點，去除掉，同時該方法采用圓形結(jié)構(gòu)作為開運算的單元［6］。Chen和Zhang的方法都會產(chǎn)生一個最小表面(Zimin)，并且對空洞采用一定的算法進行了修復，因此濾波精度也相對較好。Pingel［7］提出了改善的簡單形態(tài)學濾波算法(SMRF)，該算法采用了4個參數(shù)——坡度值(ED)、最大窗口半徑(Wmax)、高程閾值(EH)、高程縮放因子(SC)，同時采用D’Erric提出的基于彈簧質(zhì)子模型［8］用于空洞修復，漸進多尺度的濾波窗口半徑采用線性單位像素增加，濾波窗口采用圓形，點云濾波采用格網(wǎng)單位的地面點標記和原始點云的濾波，實驗表明這些措施能有效地改變?yōu)V波算法的濾波精度，具有較好的應用前景。

到目前為止，形態(tài)學濾波算法主要存在問題之一就是如何確定開算子鄰域大小，由于地形、地物分布的隨意性，固定尺寸很難滿足檢測要求，動態(tài)調(diào)整策略雖能取得不錯的效果，但需要人工指定一個最大尺寸［5］，但這在現(xiàn)實卻很難確定，需要實地調(diào)查和經(jīng)驗值，但鄰域過大或過小都會使算法濾波的精度降低，這不利于機載激光雷達數(shù)據(jù)的推廣和應用。筆者結(jié)合SMRF算法濾波精度高，效率高等特點，提出了基于多次曲線極值濾波窗口探測的形態(tài)學機載點云濾波算法，解決了需要人為確定最大窗口尺寸(即SMRF算法中濾波窗口半徑參數(shù))的問題。實驗表明該方法可行，而且濾波精度較高。

1 濾波算法描述

1.1 SMRF算法 SMRF算法是對Chen和Zhang的漸進多尺度算法的進一步改進，該算法需要輸入4個參數(shù):坡度值(ED)、最大窗口半徑(Wmax)、高程閾值(EH)、高程縮放因子(SC)，與其他的形態(tài)學濾波算法相比，其突出點主要表現(xiàn)在以下3個方面。

1.1.1 點云格網(wǎng)化。點云格網(wǎng)化以1 m作為格網(wǎng)的大小，每個格網(wǎng)用最近點最小高程值作為每個格網(wǎng)的值，SMRF算法要求最近點必須在1 m之內(nèi)的點，因而必然會出現(xiàn)一些格網(wǎng)無法填充高程值。因此采用D’Erric提出的基于彈簧質(zhì)子模型的圖形修復技術(shù)用于空洞修復［8］，這種方法相對于Laplacian技術(shù)、K鄰域的搜索的方法要好。這個過程就會形成最小的表面(Zimin)用于后續(xù)迭代開運算。

1.1.2 迭代開運算。迭代開運算的目的在于對每個格網(wǎng)進行地面格網(wǎng)和非地面格網(wǎng)的分類。定義2個變量LS和NS，分別代表開運算之前的地表面和開運算之后的地表面。由Wmax可以生成以一個像素為增量的線性多尺度形態(tài)學開運算迭代器，這相對于指數(shù)增加的開運算更能保持地面點的細節(jié)信息。在每次迭代開運算分類格網(wǎng)之前，要計算分類閾值，閾值的計算公式為:

式中，ET表示LS與NS之間對應格網(wǎng)的差值的范圍，對于差值大于ET的格網(wǎng)標記為非地面格網(wǎng)，差值小于ET的格網(wǎng)標記為地面點格網(wǎng);Wi為每次迭代的開運算半徑。迭代完成后會產(chǎn)生一個新的最小表面Zimin，產(chǎn)生新的Zimin，它代表了地面的近似表面，但是在產(chǎn)生最小表面，并沒有考慮的極低的無意義的點，這會對最終原始點云的濾波產(chǎn)生很大的影響。因此，對新的最小表面Zimin進行反轉(zhuǎn)，然后利用單一的窗口和高程閾值EH，對其進行開運算和判斷，剔除極小值的格網(wǎng)。

1.1.3 原始點云濾波。保留Zimin中地面格網(wǎng)，剔除非地面格網(wǎng)，對其進行空洞修復，在進行3次樣條插值產(chǎn)生一個臨時的DEM，原始點云的濾波閾值HY計算如下:

其中EP為每個點云對應的DEM上的近似坡度，計算原理為:

其中F是X、Y和Z之間的函數(shù)，Z=F(X，Y)，▽F為函數(shù)F的梯度值。計算每個點云的位置在DEM上的值與原始高程值的高程差絕對值看是否小于HY，小于HY的對應點為地面點，否則為非地面點。

1.2 基于高次曲線極值的濾波窗口探測的形態(tài)學機載點云濾波 在描述這個算法之前首先進行了這樣的假設:在形態(tài)學開算子用于點云濾波時，當濾波窗口半徑過大，處于地形起伏上的點就會被作為非地面點濾波掉，Ⅰ類誤差會變大，從而總誤差也會變大，當濾波窗口半徑過小時，城區(qū)的大部分建筑會被保留，從而Ⅱ類誤差會變大，總誤差也會變大。因此只有Ⅰ類誤差和Ⅱ類誤差同時較小時，總誤差也才會較小，也就說濾波窗口從小變大時，總誤差關于濾波窗口半徑的函數(shù)應該屬于一個拋物線開口向上的函數(shù)，即這個函數(shù)存在極值使得總誤差最小，這個最小誤差對應的窗口半徑就為最優(yōu)窗口半徑。但是在實際中對于像城市區(qū)域來說，大多地形比較平坦，因此當窗口過大時，對Ⅰ類誤差影響也較少，因而總誤差也是較小，如圖1所示，因而采用了高次多項式去探測最優(yōu)最大濾波窗口。改進算法后的濾波思路如圖2所示。

圖1 總誤差與最大濾波窗口的關系

在進行最大濾波窗口探測時，窗口半徑采樣必須要有一個采樣結(jié)束條件，由經(jīng)驗值可知，現(xiàn)實生活中的絕大數(shù)房屋的半徑都小于50 m［12］，因此在設計程序時可以設定一個結(jié)束閾值Hw≥50，這對結(jié)果并不會有太大的影響，因為通過實驗也發(fā)現(xiàn)通過探測優(yōu)化的最大濾波半徑都小于50。在進行探測的時候采用線性增加的方式進行，如下式所示:

式中，i=1，2，…，n。算法的具體實現(xiàn)流程見圖2。實現(xiàn)步驟如下:

(1)輸入原始數(shù)據(jù)，初始化開始濾波最大半徑Wi=1，flag=0，flag用于標記濾波是否結(jié)束。

(2)判斷Wi是否小于 Hw，如果小于，就將 Wmax=Wi，否則的話轉(zhuǎn)入第(4)步。

(3)進行SMRF濾波，結(jié)束后，判斷flag是否等于1，等于就結(jié)束濾波，否則，首先將Wi對應的總誤差Ti分別添加到向量W和T中，然后就按(5)式迭代Wi的值。

(4)判斷Wi是否大于Hw，如果是，就將W和T，按照(6)式進行最小二乘擬合，最后求得該多項式的極值，極值對應的窗口極值就為Woptim，就將Wmax=Woptim賦給，同時標記flag=1，然后轉(zhuǎn)到步驟(3)進行。

實驗表明，n取2或者13濾波結(jié)果較優(yōu)。

2 算法驗證

2.1 實驗數(shù)據(jù) 該研究實驗數(shù)據(jù)是在國際攝影測量與遙感協(xié)會網(wǎng)(ISPRS)提供8種場景LIDAR數(shù)據(jù)，專門用于測試濾波算法，其中城市和鄉(xiāng)村地區(qū)的數(shù)據(jù)各占4景;城市地區(qū)的點云間隔為1～1.5 m，而鄉(xiāng)村地區(qū)的點云間隔為2～3.5 m。這些數(shù)據(jù)包含平原、植被、建筑物、公路、鐵路、橋梁、電線、水域，陡坡等不同的地形地貌特征。從這8景數(shù)據(jù)中抽取代表不同地形特點的15個參考數(shù)據(jù)來測試算法的濾波精度，ISPRS已經(jīng)將這15個參考數(shù)據(jù)精確地分類為地面點類和非地面點類。

因此該研究對算法的測試就采用這15個人工精確分類的參考數(shù)據(jù)。

2.2 實驗結(jié)果 該研究算法用MATLAB實現(xiàn)，濾波參數(shù)參考文獻［7］，參數(shù)和實驗結(jié)果如表2所示。

圖2 改進算法流程

表2 實驗結(jié)果

圖3分別是S52、S232種樣例數(shù)據(jù)對應的原始數(shù)據(jù)、人工分離的參照數(shù)據(jù)和該算法濾波結(jié)果的浮雕圖。

2.2 結(jié)果分析 實驗結(jié)果表明，該算法不僅能免去了需要人工確定最大濾波窗口的困擾，而且也能有效的識別地面點和非地面點。如圖3所示，濾波后的地面點浮雕圖與參照圖相比，較好地保留了裸露地面原貌。在S32城區(qū)樹木被有效識別，處于大型建筑物包圍下的地面及建筑物和樹木包圍下的地面得到正確識別。在S52的山區(qū)，坡的起伏被很好的保留，說明該算法對于城區(qū)和山區(qū)都有較好的適應性。

為了更準確地評估過濾算法對結(jié)果DTM的影響，使用本文算法測試的15個參考數(shù)據(jù)的Ⅰ類、Ⅱ類和總誤差見表1;圖4是該算法與Pfeifer［9－11］和Axelsson［14］提出的分層穩(wěn)健線性迭代選權(quán)濾波和漸進三角網(wǎng)加密濾波算法以及SMRF算法的比較產(chǎn)生的Ⅰ類、Ⅱ類和總誤差圖，Pfeifer和Axelsson的算法已經(jīng)商用化，因此和它們比較分析具有重要的意義。

從Ⅰ類誤差來看，該算法相對于其他兩種算法較小，僅在S11、S52處波動較大，這說明該算法將地面點錯分為非地面點類的概率較小。在Ⅱ類誤差上，Pfeifer算法有較小的Ⅱ類誤差，而我們的算法產(chǎn)生了較大的Ⅱ類誤差，產(chǎn)生這種結(jié)果的原因是在于Pfeifer算法是自上而下的濾波，而形態(tài)學的濾波算法是自下而上的濾波，而且現(xiàn)在大多算法都具有較少的Ⅱ類誤差。Ⅰ類誤差和Ⅱ類誤差共同作用的結(jié)果一起影響著最終總誤差的精度，總誤差反映了算法的可行性，總體誤差越小表明算法的濾波效果越優(yōu)越［13］。從圖4可以看到，該研究改進后的算法在總誤差層面相對于Axelsson和Pfeifer的方法要小，僅在S23出稍微比Axelsson的算法高點。從3種誤差來看，改進的算法相比SMRF算法在某些點精度有所下降，但總體基本一致，而且在Ⅱ類誤差上，在像S52這種復雜地形地區(qū)卻要比SMRF精度要高。

經(jīng)驗表明，與Ⅰ類誤差相比，Ⅱ類誤差的點云大部分有著高于地面點云的高程，所以在人工干預時能更容易地修復Ⅱ類誤差帶來的分類錯誤。因此該研究提出的算法雖然相對于商業(yè)化的程序有較高的Ⅱ類誤差，但是Ⅰ類誤差和總誤差卻較少，從圖3也可以看出該算法也能較好地保存地形原貌。因此該算法在濾波方面精度較高，在點數(shù)據(jù)后處理中，使數(shù)據(jù)處理向著自動化處理發(fā)展，具有較好的應用前景和參考價值。

3 結(jié)論

筆者對SMRF算法進行了研究，該算法相比以往的算法具有最低的總誤差，因此機載點云濾波精度較高，但該算法需要輸入4個參數(shù)，這些參數(shù)通過優(yōu)化雖能得到較好的濾波結(jié)果，但參數(shù)輸入的不確定性，無形中提高了該算法推廣應用的門檻。因此，該研究提出了基于多次曲線極值的濾波窗口探測的形態(tài)學機載云濾波算法，該算法的主要貢獻在于結(jié)合SMRF算法效率高等特點，通過程序自動尋找最優(yōu)最大濾波窗口半徑，擺脫了以往需要人工確定多尺度形態(tài)學濾波開運算的最大濾波窗口問題。實驗表明該算法是可行的，同時在Ⅰ類誤差和總誤差層面都比商用化的算法較好且濾波精度較高，因此這對機載激光雷達點云濾波處理向著自動化發(fā)展具有重要的意義。

通過實驗，也發(fā)現(xiàn)該算法還存在的不足，主要是Ⅱ類誤差較大，因此未來的研究要對該算法進行改進，降低Ⅱ類誤差，同時還需要進一步研究降低算法輸入?yún)?shù)的門檻，主要是坡度值(ED)和高程閾值(EH)的確定。

圖3 S52、S23中改進后算法處理后的浮雕圖

圖4 幾種算法誤比較

［1］梁欣廉，張繼賢，李海濤，等.激光雷達數(shù)據(jù)特點［J］.遙感信息，2005(3):71－76.

［2］董保根，秦志遠，朱傳新，等.關于機載LiDAR點云數(shù)據(jù)形態(tài)學濾波的幾點思考［J］.測繪科學，2013，38(4):19 －21.

［3］KILIAN J，HAALA N，ENGLICH M.Capture and Evaluation of Airborne Laser Scanner Data［J］.IAPRS，1996，31(B3):383 －388.

［4］VOSSELMAN G.Slope Based Filtering of Laser Altimetry Data［J］.IAPRS，2000，33(B3/2):935 －942.

［5］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.

［6］CHEN Q.Filtering Airborne Laser Scanning Data with Morphological Methods［J］.Photogrammetric Engineering and Remote Sensing，2007，73(2):175 －185.

［7］PINGEL T J，CLARKE K C，MCBRIDE W M A.McBride.An improved simple morphological filter for the terrain classification of airborne LIDAR data［J］.ISPRS Journal of Photogrammetry and Remote Sensing，2013，77:21.

［8］D’ERRICO J.Inpaint_nans.m.［EB/OL］.http://www.mathworks.com/matlabcentral/FIleexchange/4551.(accessed 7.06.11).

［9］KRAUS K，PFEIFER N.Determination of Terrain Models in Wooded Areas with Airborne Laser ScannerData［J］.ISPRS Journal of Photogrammetry andRemote Sensing，1998，54(4):193 －203.

［10］PFEIFER N，SATDLER P，BRIESE C.Derivation of Digital Terrain Models in the SCOP Environment［C］.Stockholm，Sweden，OEEPE Workshop on Airborne Laser Scanning and Interferometric SAR for Detailed Digital Elevation Models，2001.

［11］KRAUS K，PFEIFER N.Advance DTM Generationfrom Lidar Data［J］.International Archives of Photogrammetry and Remote Sensing，2001，34(3/W4):23－35.

［12］李峰，崔希民，袁德寶.改進坡度的LiDAR點云形態(tài)學濾波算法［J］.大地測量與地球動力學，2012，32(5):128－132.

［13］吳叢叢，盧小平，李國利，等.基于TIN的LiDAR數(shù)據(jù)濾波算法研究［J］.測繪通報，2013(3):32 －35.

［14］AXELSSON P.DEM Generation from Laser Scanner Data Using Adaptive TIN Models［J］.InternationalArchives of Photogrammetry and Remote Sensing，2000，33(PartB4):110 －117.