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

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

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

1 算法原理及流程

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

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

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

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

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

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

1.2 算法流程

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

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

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

(1)載入初始點云數(shù)據(jù)，進行柵格化，生成規(guī)則分布初始網(wǎng)格DSM，格網(wǎng)的大小取1 m2。

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

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

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

1.2.3 濾波的實現(xiàn)。步驟如下:

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

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

2 算法驗證

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

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

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

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

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

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

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

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

表2 實驗結果

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

3 結論

該研究針對偏度平衡法不適應地形起伏的城區(qū)濾波的缺陷，提出了改進的方法，結合漸進形態(tài)學開運算和偏度平衡法的LiDAR數(shù)據(jù)濾波方法，并通過2個實驗進行了算法驗證。驗證結果表明，改進的算法在降低了參數(shù)輸入門檻的情況下，同時輸入較少的參數(shù)，就能夠快速地對城市區(qū)域濾波，獲得較高精度的DEM，同時還能夠保留詳細的地形起伏信息，這對于城市區(qū)域的快速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］亢曉琛，劉紀平，林祥國.多核處理器的機載激光雷達點云并行三角網(wǎng)漸進加密濾波方法［J］.測繪學報，2013，42(3):331-336.

［8］隋立春，張熠斌，張碩，等.基于漸進三角網(wǎng)的機載LiDAR點云數(shù)據(jù)濾波［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，F(xiàn)ORLANI 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］沈晶，劉紀平，林祥國.用形態(tài)學重建進行機載激光雷達數(shù)據(jù)濾波的方法研究［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］楊曉云，岑敏儀，賈洪果.基于參數(shù)統(tǒng)計的LiDAR數(shù)據(jù)分割算法［J］.測繪通報，2013(10):20-22.

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

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

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

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

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

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

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