基于漸進(jìn)形態(tài)學(xué)開運(yùn)算和偏度平衡法的LiDAR數(shù)據(jù)濾波方法

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

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

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

1 算法原理及流程

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

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

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

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

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

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

1.2 算法流程

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

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

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

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

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

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

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

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

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

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

2 算法驗(yàn)證

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

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

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

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

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

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

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

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

表2 實(shí)驗(yàn)結(jié)果

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

3 結(jié)論

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

［1］張小紅.機(jī)載激光雷達(dá)測量技術(shù)理論與方法［M］.武漢:武漢大學(xué)出版社，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ì)平，林祥國.多核處理器的機(jī)載激光雷達(dá)點(diǎn)云并行三角網(wǎng)漸進(jìn)加密濾波方法［J］.測繪學(xué)報，2013，42(3):331-336.

［8］隋立春，張熠斌，張碩，等.基于漸進(jìn)三角網(wǎng)的機(jī)載LiDAR點(diǎn)云數(shù)據(jù)濾波［J］.武漢大學(xué)學(xué)報:信息科學(xué)版，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］沈晶，劉紀(jì)平，林祥國.用形態(tài)學(xué)重建進(jìn)行機(jī)載激光雷達(dá)數(shù)據(jù)濾波的方法研究［J］.武漢大學(xué)學(xué)報:信息科學(xué)版，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］董保根，秦志遠(yuǎn)，陳靜，等.無需閾值支持的機(jī)載LiDAR點(diǎn)云數(shù)據(jù)濾波方法［J］.計算機(jī)工程與應(yīng)用，2012，49(15):219-223.

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

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

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

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

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