基于漸進形態學開運算和偏度平衡法的LiDAR數據濾波方法

林金彥，鄒時林* ，胡永杰，2

(1.東華理工大學測繪工程學院，江西南昌330013;2.新疆生產建設兵團勘測規劃設計研究院，新疆烏魯木齊830002)

機載LiDAR(Light Detection And Ranging)技術是融合了激光雷達技術，全球定位系統(GPS)和慣性導航系統(INS)為一體的現代化對地觀測的最新技術，具有受天氣影響小、分辨率高、探測范圍大、能夠部分穿透樹林遮擋，直接獲取真實地表的高精度三維信息，是快速獲取高精度地形信息的全新手段［1］，因而在DEM生產方面被廣泛使用。但是LiDAR采集的點云數據既包括地面點，又包括非地面點(房屋、植被等)，因此在生產DEM之前，有必要對地面點和非地面點進行分離，即點云濾波。從原理來看，目前點云濾波大致可以分為4 類［2］:基于形態學的方法［3-4］、基于表面的方法［5-6］、漸近加密的方法［7-8］、基于分割的方法［9-10］。從目前來看這些方法幾乎都要進行參數閾值的設置，例如Zhang等提出的漸進形態學濾波算法，需要輸入5個參數閾值，對于有一些算法雖然參數輸入較少，但參數大小設置門檻較高，難以適合不同的地形環境［11］;沈晶等用形態學重建的方法進行機載LiDAR點云濾波，參數輸入很少，僅輸入一個高差閾值參數，而且實驗結果也比較滿意，但在實際生產中，參數的輸入大小難以把握［12］;Bartels等提出了基于偏度平衡法的點云濾波算法，該算法是典型的非監督分類，無需閾值和任何先驗知識，效率高，能夠自動的完成地面點和非地面的分離［13-14］;楊曉云等提出的基于參數統計的LiDAR數據分割［15］，實質也屬于偏度平衡法，但是基于這種偏度參數統計的方法有一定的局限性，僅適用于較為平坦的區域，且假設建筑物等非地面點均高于地面點，因此對有起伏的城區，濾波會出現失效的情況。針對這種情況，筆者對該算法進行了改進，提出了結合漸進形態學開運算和偏度平衡法的LiDAR數據濾波方法，實驗表明改進的濾波算法對地形起伏的城區濾波適應性更強，濾波更為合理和穩定，能夠有效地提高濾波精度。

1 算法原理及流程

1.1 偏度系數 根據中心極限定理可知，自然狀態下量測的樣本數據服從正態分布［16］，根據這一原理可以假設:采集的點云數據中，地面點的概率密度服從正態分布，而非地面點對這種正態分布進行了干擾。因此要使點云中地面點分布達到正態分布的狀態，需要從LiDAR數據樣本中去除干擾地面點分布的非地面點，從而“校正”數據的整體概率密度分布，這種“校正”過程會持續進行直至等于0為止，此時地面點的分布即接近于正態分布。偏度的數學表達式為:

式中，sk為偏度;N為樣本總數;xi為任意樣本點，i∈{1，2，…，N}。σ為樣本標準方差，ua為算術平均值，它們由公式(2)和(3)定義:

用偏度SK描述隨機變量分布:若sk＞0，則呈正偏態分布;若sk＜0，則呈負偏態分布;若sk=0，則呈正態分布。

根據以上的思想，可以把采集的點云高程作為隨機變量，通過統計高程的偏度系數，通過迭代，不斷剔除高程最高的點，直到sk小于等于0為止，剩余的點就為裸露地表上的點。這種方法作為一種全局的濾波器，對點云的最大數量，數據格式和分辨率都沒有限制，但該算法的前提假設條件是地面點服從正態分布，因而必須要限制點云中地面點的最小數據量，如果不滿足這個條件，這個假設條件就不成立。根據文獻［13］、［14］提出的最小樣本尺寸原理，最小點數為:

式中，zα/2在α下的水平值;σ0為成功濾波后的點云的標準方差;E為所有點云數據與濾波后點云數據之間最小距離所容許的最小幅度。σ0由于事先未知，所以根據經驗值可以取0.75，假設允許的最小幅度值E為0.05，則最小點數大概為866，顯然對于城區來講，大面積的掃描，這是滿足要求的［17］。

在滿足假設的條件下，該算法在實際實驗中會出現過分割現象或者分割失敗的情況，尤其對一些有起伏的城區，濾波精度明顯下降［18-19］。產生這種情況的原因在于該算法僅考慮了完全是由像房子、橋梁、植被等地物引起的高程突變，也就是假設建筑物等非地面點均高于地面點，并沒有考慮像陡坎、斜坡等地形起伏對高程突變的影響。針對該算法所存在的這種不足，該研究采用一種轉換的方式，即把由地形引起的高程突變轉換成平地，轉換方法是:首先通過漸進形態學開運算的方法得到近似DEM，然后求得每個離散點到近似DEM的高差，這樣就使由地形引起的高程突變轉換成平地了，這樣就滿足了偏度平衡法的假設條件——建筑物等非地面點均高于地面點，即滿足了高程突變不是由地形引起的條件。

1.2 算法流程

1.2.1 數據預處理。數據預處理的目的在于剔除由于系統原因產生的高程過高過低點，該算法進行這一步的目的在于LiDAR數據虛擬格網化時，提高每個格網最低點為地面點的概率。基于偏度參數統計的方法，假設地面點近似服從正態分布，由于局外點的高程值要么遠遠大于周圍的激光腳點高程，要么遠遠小于周圍的激光腳點高程，因此根據高程均值和方差來估計和判斷激光腳點是否屬于局外點［20］。局外點的判斷閾值Zmin和Zmax計算公式如下:

式中，Zmean是高程均值;σ2是方差。對于高程小于Zmin或大于Zmax的點作為過低過高的點從數據集中剔除，為下一步的數據處理做準備。

1.2.2 近似DEM的獲取。這里采用漸進形態學開算子用于裸露地表的獲取，這是在于漸進形態學開運算通過不斷地調整窗口大小，不僅能剔除較小的地物點，也能較好地剔除較大的地物點［11］，同時形態學開算子運算效率高，過程如下:

(1)載入初始點云數據，進行柵格化，生成規則分布初始網格DSM，格網的大小取1 m2。

(2)對DSM不斷的增大窗口執行開運算，直到窗口尺寸達到最大Wmax，最大窗口一般取最大建筑物的寬度，而現實生活中絕大部分房屋寬度都小50 m［21］，最大建筑物長度一般取50～100 m較為合適［22］，最小窗口尺寸取1 m，窗口的迭代根據下式:

式中，k為迭代次數，k=1，2，…，M;wk即第k次濾波的次數;w0為初始窗口的尺寸大小;wk為第k次窗口的尺寸大小。開運算的根據下式獲得:

式中，Z為DSM柵格值的大小;g為窗口大小;(Zog)為開運算，即對原始離散點云插值出的數字表面模型先執行腐蝕操作，接下來進行膨脹操作。先運用腐蝕運算剔除尺寸小于窗口的小的地物點(離散的樹木、車等)，再用膨脹運算把誤作地物點剔除的大物體的邊界進行恢復。

1.2.3 濾波的實現。步驟如下:

(1)根據三次樣曲線和近似DEM，對每個點插值求得每個點的高差為V。

(2)對V求偏度值sk，當sk＞0時，剔出高差最大的那個值，循環迭代，直到sk≤0，那么剩余的高差值對應的點為地面點。

2 算法驗證

算法驗證采用MATLAB實現，可視化利用了Quick_Terrain_Reader和ArcGIS10.0軟件。為了充分驗證算法的可行性和合理性，這里采用了2個實驗，實驗一以城市區域簡單房屋為例;實驗二以復雜城區為例。

2.1 改進算法合理性驗證 這里采用ISPRS(國際攝影測量與遙感協會網)提供的某城區的一塊以房屋為主的區域為例，該區域的一個最重要的特點就是存在復雜四合院式的房屋，房屋周圍植被多，同時較四合院高，地形相對比較平坦，如圖1a所示，該數據針對每個點進行了精確的人工分離，人工提取的地面點如圖1b所示。

將該數據用該研究提出的改進算法處理，得到的結果如圖2b所示。

圖2a為偏度平衡法直接對原始點云處理的結果，比較圖2a跟圖1b，明顯可以看出，在紅色圈里的房屋并沒有被剔除掉，這主要是由于在圖1a中箭頭所指的房屋較紅色圈里的房屋高，這導致了偏度系數過早收斂，這導致較低的地物不能被有效的剔除，如圖2a中箭頭所指的較低矮的房屋任仍然沒有被去除掉。但對該算法進行改進，實驗效果明顯，如圖2b所示，比較圖2b、圖2a、圖1b發現，改進的算法能夠有效剔除非地面點建筑、植被等，而且反映了改進的算法能夠剔除低矮物體，同時減緩了偏度平衡法偏度系統的收斂，有利于提高濾波精度。

通過該實驗也反映了偏度平衡及其在平坦的城市區域濾波也有可能出現失效或者達不到預期的效果，從側面反映了偏度平衡不具有廣泛的適應性。

2.2 改進算法適應性驗證 這里選取了ISPRS提供的3組復雜城市區域用于改進算法適應性驗證，實驗數據描述見表1。實驗區A、B、C屬于城市區域，整體地形有起伏，不存在絕對水平的區域，實驗區內幾乎包括城市區域的所有特征地物，如橋梁、植被、房屋等，這有利于驗證改進算法的適應性。

對3組數據分別用偏度平衡法和改進的算法處理，處理的結果采用ISPRS提出的3類誤差的評價方法:Ⅰ類誤差，地面點被錯分為非地面點的概率;Ⅱ類誤差，非地面點錯分為地面點的概率;總誤差，反映了分類結果與參考數據不一致的概率。處理結果見表2。

由表2可知，改進后的偏度平衡算法總誤差明顯降低，而且較改進前算法更為穩定。較低的總誤差，可以反映該算法有較好的適應性，II類誤差也有較大的降低，實驗區B降低了61.13%，實驗C區降低了30.47%，這可以證明該研究改進后的算法能夠較好地剔除非地面點(植被、房屋等)，從另一面可以說明該算法能夠較好地適應有地形起伏的城市區域的濾波。

表2 實驗結果

圖3為濾波前后的可視化圖，從整體上看，濾波處理后能剔除幾乎所有的房屋、植被，橋梁等非地面點，從濾波后B地區的浮雕圖中，可以清晰地看出地形起伏的輪廓，這說明改進后的算法可以較好地保留地形細節信息。

3 結論

該研究針對偏度平衡法不適應地形起伏的城區濾波的缺陷，提出了改進的方法，結合漸進形態學開運算和偏度平衡法的LiDAR數據濾波方法，并通過2個實驗進行了算法驗證。驗證結果表明，改進的算法在降低了參數輸入門檻的情況下，同時輸入較少的參數，就能夠快速地對城市區域濾波，獲得較高精度的DEM，同時還能夠保留詳細的地形起伏信息，這對于城市區域的快速DEM提取有著重要的參考價值和實際意義。

［1］張小紅.機載激光雷達測量技術理論與方法［M］.武漢:武漢大學出版社，2007:8-11.

［2］SHAN J，TOTH C K.Topographic laser ranging and scanning:Principles and processing［M］.Publisher:CRC Press，2008:312-322.

［3］PINGEL T J，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］CHEN Q，GONG 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］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，54(4):193-203.

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

［7］亢曉琛，劉紀平，林祥國.多核處理器的機載激光雷達點云并行三角網漸進加密濾波方法［J］.測繪學報，2013，42(3):331-336.

［8］隋立春，張熠斌，張碩，等.基于漸進三角網的機載LiDAR點云數據濾波［J］.武漢大學學報:信息科學版，2011，36(10):167-170.

［9］SITHOLE G.Segmentation and classification of airborne laser scanner data［D］.TU Delft:Aerospace Engineering，2005:35-68.

［10］NARDINOCCHIC C，FORLANI G，ZINGARETTI P.Classification and Filtering of laser data［C］//MIAPRS.Germany Dresden，2003:245-251.

［11］ZHANG K，CHEN S C，WHITMAN D，et al.Aprogressive morphological filter for removing nonground measurements from airborne LiDAR data［J］.IEEE Trans.Geoscience and Remote Sensing，2003，41(4):872-882.

［12］沈晶，劉紀平，林祥國.用形態學重建進行機載激光雷達數據濾波的方法研究［J］.武漢大學學報:信息科學版，2011，36(2):167-170.

［13］BARTELS M，WEI H，MASON D C.DTM generation from LiDAR data using Skewness Balancing［C］//Proceedings of the 18th International Conference on Pattern Recognition.Hong Kong:Pattern Recognition，2006:566-569.

［14］BARTELS M，WEI H.Threshold-free object and ground point separation in LiDAR data［J］.Pattern Recognition Letters，2010，31(10):1089-1099.

［15］楊曉云，岑敏儀，賈洪果.基于參數統計的LiDAR數據分割算法［J］.測繪通報，2013(10):20-22.

［16］DUDA R O，HART P E，STORK D G.Pattern classification［M］.New York:Wiley，2001.

［17］董保根，秦志遠，陳靜，等.無需閾值支持的機載LiDAR點云數據濾波方法［J］.計算機工程與應用，2012，49(15):219-223.

［18］龔亮，張永生，施群山，等.基于高程統計方法的機載LiDAR點云數據濾波［J］.測繪與空間地理信息，2012，35(2):42-45.

［19］萬劍華，黃榮剛，周行，等.基于曲率統計的LiDAR點云二次濾波方法［J］.中國石油大學學報:自然科學版，2013，37(1):56-60.

［20］賴旭東.記載激光雷達基礎原理與應用［M］.北京:電子工業出版社，2010:172-175.

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

［22］張褶斌.機載LiDAR點云數據處理理論及技術研究［D］.西安:長安大學，2010:33-34.