時培強,江 虹
(西南科技大學 信息工程學院,四川 綿陽 621010)
激光雷達(Light Detection and Ranging,LiDAR)是一種應用全球定位系統(GPS)﹑慣性測量單元(IMU)﹑激光測距儀和計算機技術于一體的系統[1]。與傳統航測相比,機載LiDAR具有精度高﹑速度快﹑受外界因素影響小﹑自動化程度高等技術優點,且能提供高精度﹑高密度的三維信息[2]。LiDAR在三維數字城市建模﹑城市規劃﹑海岸線檢測﹑電力和公路選線巡線﹑森林資源的管理和評估等領域,具有廣闊的應用前景[3]。
激光雷達數據由地面點和地物點組合而成。濾波是指從離散的點云數據中將地面點和非地面點進行分離[4]。由于點云數據呈現出的無規律性和不確定性,濾波一直是激光雷達數據處理的首要問題,且其精度直接影響后續工作的準確性。
典型算法主要分為:三角網濾波算法[5-6]﹑曲面擬合濾波算法[7-10]﹑數學形態學濾波算法[11]﹑移動窗口法和基于坡度的濾波算法[12]等。ISPRS第Ⅲ委員會對不同的濾波算法進行了全面分析,結果表明:多數濾波算法對特定地形可以得到理想的效果,但對于多種復雜地形濾波效果不佳[13-14]。三角網濾波算法整體適應性較好,但計算量和復雜度相對較大,且會丟棄部分地面點;曲面擬合算法對地形特征保護較好,算法相對簡單,但種子點的不當選取會影響濾波精度;數學形態學算法簡單易實現,計算效率高,需要根據實際數據確定窗口大小;移動窗口算法的濾波適用性相對較差;基于坡度的濾波算法主要依賴初始坡度的設定。實際研究中發現,可以將均值限差法﹑曲面擬合算法與三角形角度限制法結合,通過偏度和峰度的變化對閾值進行相應調整,從而提高濾波效果。該方法以均值限差法獲取可靠性較高的地面點為基礎進行曲面擬合,通過對不斷更新的區域進行曲面擬合,完成“粗分類”“細分類”濾波,然后將剩余點與鄰近已確定的地面點通過三角形角度限制法判斷,利用各方法的優點進行逐級濾波,最終使濾波結果更加合理。
本文的濾波算法是一個逐級濾波處理過程,思路為:區域廣且點云密度大時,如果對整體進行處理會影響濾波的效率和準確性,所以本文采用數據格網化,并進行區域分塊。首先對點云數據進行格網化分,其次對整個區域進行分塊,分塊后的各區域應包含地面點,最后將分塊后的區域再分成均等的若干小區域,在小區域內選擇最低點,并利用均值限差法判斷當前最低點是否為地面點。如果當前區域選擇的地面點數大于二次曲面擬合方程求解系數個數,則對所選的地面點進行曲面擬合,并利用擬合曲面對此區域的點云數據實施Ⅱ類誤差最小化的“粗分類”過程。通過對加入點前后兩次偏度和峰度的計算,根據偏度和峰度的變化對該區域中的閾值進行動態調整。當滿足地面閾值條件時,將點加入到地面點集合;反之,當滿足地物閾值條件時,將點加入地物點集合。
經過濾波后,大部分建筑物﹑植被和高程極值點被濾除。然后,分別對小區域以“粗分類”的地面點結果為基礎選擇下一級地面點。在小區域內進行曲面擬合,利用擬合曲面方程對小區域進行“細分類”。根據偏度和峰度的變化對閾值進行相應調整。當滿足地面閾值條件時,將點加入地面點集合;反之,當滿足地物閾值條件時,將點加入到地物點集合。遍歷所有分塊,直到整個區域被覆蓋,從而實現點云數據的初次分類。最后,對部分剩余的點通過與其鄰近且已確定的地面點利用三角形角度限制判別法再次進行分類。角度閾值的設定需要根據峰度和偏度的變化作相應調整,將判斷點之間形成的角度與設定的角度閾值進行比較。當滿足地面閾值條件時,將其加入地面點集合;反之,加入到地物點集合,從而實現點云數據的整體分類。濾波流程如圖1所示。

圖1 濾波算法流程
統計學中,從樣本的概率密度曲線直觀看來,偏度(Skewness)就是其尾部的相對長度[15],峰度(Kurtosis)就是其尾部的厚度。通常,偏度和峰度分別用sk和ku表示。若某一隨機變量X的三階矩和四階矩同時存在,則計算公式為:

其中N表示樣本總數,xi為隨機樣本點,σ和μa分別表示樣本的標準差和算術平均值,定義為:


偏度和峰度均是為無量綱的量。若sk>0,則稱該分布為正偏態或者右偏態;若sk<0,則稱該分布為負偏態或者左偏態。|sk|越大,表示其偏離程度越大。類似地,若ku>0,則該分布比較陡峭;若ku<0,則該分布比較平坦。偏度和峰度均是相對于正態分布來比較的。正態分布的偏度和峰度均為0。偏度和正態分布如圖2所示。本文利用偏度﹑峰度的變化對閾值進行相應調整,通過比較新增加點云數據對偏度和峰度的影響,定位到新增點云的類型。如果前后偏度差小于0,表示新增點接近均值,更偏向于是地面點。如果偏度差大于0,表示新增點偏離均值較大,更偏向于地物點。如果前后峰度差大于0,表示新增點云與當前點云數據偏差不大,偏向于是地面點。如果峰度差小于0,表示新增點云與當前數據偏差較大,偏向于地物點。通過偏度﹑峰度的變化量對閾值進行調整,將高程差值與閾值進行比較,決定新增點云是地面點還是地物點,以提高濾波的精度。


圖2 不對稱分布及正態分布
數據格網化是將所測區域的點云數據用規則的格網全部覆蓋,如圖3所示。

圖3 格網化分
首先統計所測區域點云總數N,并得到X方向的最小值XMin﹑最大值YMax以及Y方向的最小值YMin﹑最大值YMax。所需矩形區域的面積大小為Area=(XMax-XMin)(YMax-YMin)。由點云總數N和區域面積Area計算點云密度,根據點云密度確定格網間隔,使點云數據分布均勻。假設根據以上信息得到規則格網的大小為m×n,則某點的坐標信息(Xk,Yk,Zk)所對應的行列號(i, j)為:

如圖4所示,根據某點(圖中標注點)的坐標信息可以得到所在格網位置。每個格網需要記錄該格網內包含的點數和當前點所對應的序列號,通過序列號可以快速索引到該點包含的信息。

圖4 格網索引
均值限差法是以選中點為中心,并設定窗口大小,計算窗口內所有點的均值與中心點之間的差值,通過比較差值與設定的閾值判斷中心點的取舍。本文根據非地形點的高程大于地形點的高程且非地形點的高程一般偏離均值較大的特性,計算此區域內的點云高程平均值,然后計算選中點與平均值之差,通過與設定的閾值進行比較,與平均值高差大于設定閾值的點被視為非地面點,如圖5所示。由于點云數據包含一些高程極高點或極低點,而這些極點并不是真實數據,會影響濾波精度,需剔除這些誤差點。本文采用此方法通過剔除極低點來選擇初始地面種子點。首先選擇局部區域的最低點為中心點,以此點為中心計算其周圍若干格內的點云高程均值,然后得到中心點高程與均值的差值。通過比較差值與設定閾值確定此點是否為地面點,從而排除極低點的干擾。此方法可以有效得到最初地面種子點,從而為曲面擬合提供準確﹑可靠的地面點。

圖5 均值限差法
由于地形表面比較復雜,可以看作是一個有多個獨立的曲面組合而成的空間曲面[16]。二次曲面擬合對地形特征保護較好,通過部分地面點可得到局部地形的總體特征,利用地形特征能夠有效辨別鄰域內的點是否屬于地面點。每一個獨立的曲面可以用一個二次曲面進行擬合,其擬合方程為:

式(8)中,(xn, yn, zn)為點云n 的3維坐標值,式(9)中X為二次多項式的系數,即待求參數。式(7)擬合方程用矩陣表示為Ax=b,其一般包含m個等式和n個未知數,且m>n。該方程組一般無解,為了得到盡可能合適的X且使等式盡量成立,引入殘差平方和函數:

式(11)中s(x)取最小值:

通過對式(11)中的s(x)進行微分求最值,得到:

當方程個數大于方程系數個數時,由此求出方程系數,將相應區域的點云數據代入此方程等到相應的高程值。同時,更新所在區域的偏度和峰度,通過偏度和峰度的變化量相應調整閾值。原始點云高程和擬合曲面高程的差值小于地面限差值的點被視為地面點;反之,差值大于地物限差值的點被視為非地面點。
當高程差值不符合限差條件時,當前點需要進行下一步的細分類,如圖6所示。當所測區域的窗口較大時,此時偏度和峰度的變化量相對較大,對閾值的設定權重需要降低,設置的限差值應較??;反之,當所測區域的窗口較小時,偏度和峰度的變化量相對較小,對閾值的設定權重較大,設定的限差值應較大。閾值的動態設置可以提高算法的適應性。

圖6 曲面擬合
當水平距離較近時,垂直高度差越小,形成的空間立體角越小;反之,水平距離較近,垂直高度差越大,形成的空間立體角越大,如圖7所示。

圖7 角度限制法
文中通過偏度和峰度的變化量對角度閾值進行相應調整。當前判斷點與鄰近已確定地面點構成三角形。當二者構成的角度大于設定的角度閾值時,則當前判斷點被視為非地面點;反之,當角度小于設定的角度閾值時,則當前點被視為地面點。空間立體角的計算公式為:

采用ISPRS第Ⅲ委員會在2003年提供的數據對本文算法的適應性﹑可行性進行驗證。本文選取包含高大建筑物﹑地物﹑低矮植被等分布的地區(Sample12),包含激光點云總數52 119個;另一區域包含高大建筑物﹑植被﹑建筑物與植被混合區(Sample31),包含激光點云總數28 862個。
濾波評價標準包含Ⅰ類誤差﹑Ⅱ類誤差以及總誤差三種。Ⅰ類誤差是指將地面點誤判為非地面點的比率,Ⅱ類誤差是指將非地面點誤判為地面點的比率,總誤差是由Ⅰ類誤差與Ⅱ類誤差加權求和得到的[17]。其中,Ⅰ﹑Ⅱ類誤差表現算法的適應性,總誤差反映算法的可行性。
利用本文所提的算法分別對Sample12區域和Sample31區域進行濾波處理??梢园l現,在城市區域濾波效果較好,在地形起伏較大的山區或地物點與地面點高程相近處出現部分分類誤差。如圖8﹑圖9所示,通過對濾波效果進行分析,可以得出以下結論:
(1)對于平坦的道路,本文算法能夠比較準確地從點云數據中提取出道路附近的地面點。由于道路起伏較小,曲面擬合可以獲取較為準確的原始地形信息。當道路旁低矮的植被與地面點接近時,會將部分低矮植被點誤分為地面點。在保證Ⅱ類誤差最小化的前提下,通過偏度和峰度的變化對閾值進行相應調整,能夠比較準確地分離出地面點。
(2)對于高大建筑物和植被,本文算法能夠從點云數據中將絕大部分建筑物和植被點提取出來,并可對建筑物和植被點進行準確分類。對于建筑物的起始邊緣或低矮植被與地面接近的情況,此時會影響部分建筑物邊緣點和低矮植被的誤分類,如圖10所示。隨著相差距離的增大,分類結果會越來越準確。
(3)對于建筑物與植被混合區,本文算法能夠將大部分點云正確分類,但仍存在部分誤分類情況,如圖11所示。在建筑物與植被之間會存在少量的地面點或低矮植被。對于此部分點云,在經過粗分類和細分類后不能準確分類。根據此區域的偏度和峰度的變化對角度閾值進行調整,對此部分點云使用角度限制法可以有效提高分類效果。和峰度的變化不明顯,會將部分低矮植被點誤分為地面點,從而影響了Ⅱ類誤差。為了降低Ⅱ類誤差,本文將地面點閾值區間設置較小。當地形出現明顯突變時,不可避免會將此部分地面點誤分為地物點,從而增加了Ⅰ類誤差率。

圖8 Sample12區域濾波結果對比


圖9 Sample31區域濾波結果對比

圖10 Sample12區域局部誤差

表3是本文算法的誤差與部分經典算法的誤差之間的對比??梢园l現,本文算法在測試區域總誤差處于中等以上水平,可知所提算法具有一定的可行性,并能夠得到可靠的濾波效果。

圖11 Sample31區域局部誤差
從表1和表2可以看出,通過對濾波誤差進行分析,可以得出以下結論:
(1)不同的測試區域雖然三類誤差有所區別,但三類誤差中,Ⅱ類誤差在不同測試區域都是最低的。這驗證了本文算法是在降低Ⅱ類誤差的前提下,進一步提升了整體的濾波精度。
(2)Ⅱ類誤差出現的主要原因是地面上有低矮植被,當低矮植被與地面點高程差很小時,偏度

表1 Sample12區域濾波結果統計

表2 Sample31區域濾波結果統計

表3 濾波結果比較
本文提出一種多級相結合的LiDAR點云數據濾波算法。主要工作和結論如下:
(1)以均值限差法獲取可靠性較高的初始地面種子點為基礎,提高曲面擬合的準確性。以曲面擬合算法對數據進行濾波,能有效剔除地物點,同時保留原始地貌,且算法簡單。
(2)以三角形角度限制法對剩余點進行判斷,可以很好地利用坡度信息和已經確定為地面點的數據進行處理,進一步提高濾波的精度。
(3)通過實驗驗證,該濾波算法通過偏度和峰度的變化對所在區域的閾值進行相應調整,對含有高大建筑物和植被的城市區域有較好的適應性,且取得了理想的濾波效果。
盡管本文算法的可行性和實用性較高,偏度和峰度的變化在城市區域有著較高的適應性,能夠對閾值進行相應調整,但對于地形起伏較大的山區效果不顯著,濾波精度較低。此外,區域分塊需要根據建筑物的大小劃分,自適應性也需要進一步完善。以后的研究需要進一步對算法進行改進,并對已確定的地物點進行插值處理,以便更好地保持地形原貌,也可以利用回波強度信息進一步提高濾波精度。
[1] 賴旭東.機載激光雷達基礎原理與應用[M].北京:電子工業出版社,2010.
LAI Xu-dong.Basic Principles and Applications of Airborne Lidar[M].Beijing:Publishing House of Electronics Industry,2010.
[2] 張小紅.機載激光雷達測量技術理論與方法[M].武漢:武漢大學出版社,2007.
ZHANG Xiao-hong.Airborne Lidar Measurement Technology Theory and Method[M].Wuhan:Wuhan University Press,2007.
[3] 張玉方,程新文,歐陽平等.機載LIDAR數據處理及其應用綜述[J].工程地球物理學報,2008,5(01):119-124.
ZHANG Yu-fang,CHENG Xin-wen,OUYANG Ping,et al.The Data Processing Technology and Application of Airborne LIDAR[J].Chinese Journal of Engineering Geophysics,2008,5(01):119-124.
[4] 劉經南,許曉東,張小紅等.機載激光掃描測高數據分層迭代選權濾波方法及其質量評價[J].武漢大學學報:信息科學版,2008,33(06):551-555.
LIU Jing-nan,XU Xiao-dong,ZHANG Xiao-hong,et al.Adaptive Hierarchical and Weighted Iterative Filtering of Airborne LIDAR Data and Its Quality Assessment[J].Geomatics and Information Science of Wuhan University,2008,33(06):551-555.
[5] Axelsson P.DEM Generation from Laser Scanner Data Using Adaptive TIN Models[J].International Archives of Photogrammetry and Remote Sensing,2000,33(B4/1;PART 4):111-118.
[6] 曾妮紅,岳迎春,魏占營等.車載LiDAR點云濾波的改進不規則三角網加密方法[J].測繪科學,2016,41(09):136-139.
ZENG Ni-hong,YUE Ying-chun,WEI Zhan-ying,et al.An Improved Irregular Triangular Network Encryption Method of Vehicle-borne LiDAR Point Clouds[J].Science of Surveying and Mapping,2016,41(09):136-139.
[7] 張小紅.機載激光掃描測高數據濾波及地物提取[D].武漢:武漢大學,2002.
ZHANG Xiao-hong.Airborne Laser Scanning Altimetry Data Filtering and Features Extraction[D].Wuhan:Wuhan University,2002.
[8] 熊俊華,方源敏,鄧德標.最小二乘曲面擬合的LiDAR數據濾波方法[J].測繪科學,2013,38(04):74-76.
XIONG Jun-hua,FANG Yuan-min,DENG De-biao,et al.Surface Fitting Filtering based on Least Square Method from LiDAR Data[J].Science of Surveying and Mapping,2013,38(04):74-76.
[9] 孫崇利,蘇偉,武紅敢等.改進的多級移動曲面擬合激光雷達數據濾波方法[J].紅外與激光工程,2013,42(02):349-354.
SUN Chong-li,SU Wei,WU Hong-gan,et al.Improved Hierarchical Moving Curved Filtering Method of LiDAR Data[J].Infrared and Laser Engineering,2013,42(02):349-354
[10] 張皓,張永生,劉軍等.一種基于平面擬合的LIDAR點云濾波方法[J].測繪科學,2009,34(04):141-143.
ZHANG Hao,ZHANG Yong-sheng,LIU Jun,et al.A Method for Filtering LIDAR Points Cloud based on Planar Fitting[J].Science of Surveying and Mapping,2009,34(04):141-143.
[11] 董保根,秦志遠,朱傳新等.關于機載LiDAR點云數據形態學濾波的幾點思考[J].測繪科學,2013,38(04):23-25.
DONG Bao-gen,QIN Zhi-yuan,ZHU Chuan-xin,et al.Perspectives of Morphological Filtering for Airborne LiDAR Point Clouds Data[J].Science of Surveying and Mapping,2013,38(04):23-25.
[12] Vosselman G.Slope based Filtering of Laser Altimetry Data[J].International Archives of Photogrammetry and Remote Sensing,2000,33(B3/2;PART 3):935-942.
[13] Sithole G,Vosselman G.Experimental Comparison of Filter Algorithms for Bare-Earth Extraction from Airborne Laser Scanning Point Clouds[J].Isprs Journal of Photogrammetry& Remote Sensing,2004,59(01-02):85-101.
[14] Sithole G,Vosselman G.Filtering of Airborne Laser Scanner Data based on Segmented Point Clouds[C].Laser Scanning,2005.
[15] 董保根,秦志遠,陳靜等.無需閾值支持的機載LiDAR點云數據濾波方法[J].計算機工程與應用,2013,49(15):219-223.
DONG Bao-gen,QIN Zhi-yuan,CHEN Jing,et al.Threshold-free Method for Airborne LiDAR Point Clouds Data Filtering[J].Computer Engineering and Appl ications,2013,49(15):219-223.
[16] 劉志青,李鵬程,郭海濤等.融合強閾值三角網與總體最小二乘曲面擬合濾波[J].紅外與激光工程,2016,45(04):128-135.
LIU Zhi-qing,LI Peng-cheng,GUO Hai-tao,et al.Integrating Strict Threshold Triangular Irregular Networks and Curved Fitting based on Total Least Squares for Filtering Method[J].Infrared and Laser Engin eering,2016,45(04):128-135.
[17] 黃先鋒,李卉,王瀟等.機載LiDAR數據濾波方法評述[J].測繪學報,2009,38(05):466-469.
HUANG Xian-feng,LI Hui,WANG Xiao,et al.Filter Algorithms of Airborne LiDAR Data:Review and Prospects[J].Acta Geodaetica et Cartographica Sinica,2009,38(05):466-469.