基于高程統計的機載LiDAR點云三角網漸進濾波方法

陳琳，范湘濤，杜小平

（1.中國科學院遙感與數字地球研究所 數字地球重點實驗室，北京 100094；2.中國科學院大學，北京 100049）

1 引 言

激光掃描測距技術（Light Detection And Ranging，LiDAR）是一種快速直接獲取地表模型的技術［1］。機載LiDAR點云數據包括地面點之外，還有樹木、建筑物、橋梁、車輛等信息，獲取真實的數字高程模型需要剔除這些非地面點，即點云數據濾波［2］。在機載LiDAR數據處理過程中，濾波是最關鍵的部分之一，是進行后續高層次數據處理的基礎［3］。準確的濾波結果才能保證各種應用的有效性。

經過十幾年的發展，濾波算法已經取得一定的研究成果。現有濾波方法根據技術路線差異，主要分為3類：形態學方法、基于內插的方法和基于曲面約束的方法［4］。基于形態學的點云濾波方法使用與被測地物相似的結構元，進行腐蝕、膨脹等形態學運算過濾地物點數據。眾多學者開展了此方面的研究［5－8］，這類方法比較適用于城市區域，結果依賴移動窗口尺寸，對非平坦地區效果并不理想。基于內插的點云濾波方法首先建立粗糙的初始DEM，然后逐漸內插加密DEM，去除地物點。其中Axelsson提出的不規則三角網加密算法是比較有代表性的方法［9－11］，該方法在城市和森林地區都有較好的濾波效果，并且已經應用到商業軟件TerraScan中。基于曲面約束的方法假設復雜地形在局部地區可以用簡單的曲面表達，進而排除地物點。曲面的種類有樣條曲面［12］、正交多項式曲面［13］和移動二次曲面［14－15］等。此類方法計算相對復雜，擬合過程會造成誤差累積，較適合平緩地形。由于地形地物的復雜性，沒有一種算法可以完全適用。

2 方 法

2.1 基于高程統計的閾值分割

不同地物要素會導致局部地區點云的高程分布差異，可通過高程統計得到高程閾值進行數據分割。大部分數理統計量與數據的空間位置和數據格式無關［17］，基于統計的處理方法普適性強。

高程直方圖反映點云數據的高程分布。對比高程直方圖，迭代計算點云數據的高程閾值，得到的最佳閾值將點云數據分割為兩個高程差異明顯的數據集。當研究區地形較為平緩、地物與地面差異明顯時，單次分割即可將地物與地面點有效分離；反之，需要對已經分割得到的數據集進行多次分割，逐級分離地物與地面點。其具體處理流程為：

①繪制高程直方圖。點云數據的高程直方圖以高程為橫坐標，頻數為縱坐標，反映各高程與其出現頻數之間的關系。

②參數閾值計算。計算原理根據文獻［18］：查找待分類點云的最大高程值gmax和最小高程值gmin，取兩者平均值為初始閾值g0；g0把點云分成大于g0和小于g0的兩部分，g1為兩部分期望值的平均值；如果g0與g1之差的絕對值大于自定義的極小值，則將g1賦值給g0，重復上一步的閾值計算，反之，則g1是分割的最佳閾值，退出計算。其中，極小值是一盡可能小的實數，以保證最終得到的最佳閾值趨于穩定（在本實驗中自定義極小值取為0.01）。

③多閾值分割。如果待分類點云中有多個不同高度的目標，可對已分割結果進行多次分割，逐步將地物點剔除。

2.2 三角網漸進濾波

三角網漸進濾波算法默認局部地區的最低點為地面點，由這些最低點創建初始不規則三角網（Irregular Triangle Network，TIN）地形，通過計算點到初始地形的距離和角度判定點的類別，其具體流程為：

①參數設置。參數包括：研究區最大建筑物尺寸、點到三角面的最大高程閾值和點到三角面節點的最大角度閾值。

1)小范圍內精細化導航，避免游客迷路。AR導航更直觀、準確，能讓步行的用戶知曉自己周邊位置場景的更多細節。

②選擇種子點。以最大建筑物尺寸劃分格網后，查找每個格網中的最小高程點，由這些點創建初始TIN。

③加密三角網。查找激光點，如果點到它落入的三角面片的距離和角度滿足閾值參數，則判定該點為地面點，并加密到三角網中，反復迭代直到沒有新的點被判定為地面點為止。

2.3 基于高程統計的三角網漸進濾波

基于高程統計的閾值分割方法簡單易行，但由于缺少空間約束條件，濾波得到的地面點集會混雜較多的低矮地物點。而三角網漸進濾波方法缺少先驗知識，在選擇種子點時格網尺寸受到建筑物尺寸限制。本文方法引入高程統計分割結果作為先驗知識，在此基礎上劃分格網時，格網尺寸不再受建筑物尺寸約束，可以定義小尺寸的格網以提取更多初始地面點，同時保證初始值的準確度，得到更加細致的初始地形。

由于先驗知識的準確度會影響后期三角網濾波的精度，對點云數據進行高程統計分析時有三個注意點。第一，研究區范圍較大或地形劇烈變化時，低地處具有較高高程的地物點與高地處的地面點可能有相同的高程值（“同高異類”），需要對研究區按照高程和平面位置進行分區統計。具體操作時，將原始數據按照橫向和縱向分別二等分，得到4個次級尺寸的區域，如果小尺寸區域內仍存在“同高異類”現象，繼續對小尺寸區域按照橫向和縱向進行二等分，直到分區內不存在“同高異類”為止；同時，如果相鄰的小區域內點云高程分布相近，地形變化相對平緩，也可將相鄰小區域進行合并。第二，為了避免地物邊緣點混入分割結果中造成初始點選擇誤差，在進行分割時，可以適當調低分割閾值，舍棄部分地面點，以保證統計分割結果的準確率。第三，在微小的地形變化中，單純依靠高程統計分割方法很難剔除低矮植被點，但由于植被是略高于局部地表的，按照選擇極低點的假設條件，這些點不會影響種子點的精度。

由于高程統計分割方法在建筑物密集的城區具有較高的準確性，而在地形起伏較大的林區和山區準確性相對較低，所以基于高程統計的三角網漸進濾波方法更適合于城區點云數據。

該方法濾波主要流程：

①剔除系統粗差點。依次查找每個點一定范圍內（以5m～8m為半徑的圓周內）的其他激光點，如果該點比任意其他激光點都低1m以上，則將該點判定為粗差點，予以剔除。此處的1m為經驗值，可以根據地形起伏情況，手動調整。

②繪制點云數據的高程直方圖，統計點云高程得到閾值參數，對照直方圖進行高程分割。當研究區范圍過大或地形變化較劇烈時，按照地形特征將點云數據分成范圍較小的多個區域，對每個子區域的點云數據執行高程統計和閾值分割后，合并地面點分割結果。

③對分割得到的備選地面點劃分格網，提取每個格網內的最低點作為種子點，創建初始三角網地形。

④計算粗差處理得到的所有激光點到初始地形的距離和角度，滿足閾值參數的點被加密到三角網中進行下一次迭代過程，當沒有新的地面點加入到三角網時，迭代計算結束，得到濾波結果。濾波流程如圖1所示。

圖1 基于高程統計的三角網漸進濾波流程圖

3 數據與實驗

為了檢驗此算法，使用ISPRS公布的濾波實驗數據，其中Site2數據覆蓋范圍廣，用以解釋分區原理，其他3組數據有半自動半手工分類結果作為真值，與原始三角網漸進濾波方法相比進行精度驗證。

3.1 實驗數據

所有數據都為城區數據，點密度為0.67pts／m2，點間距為1.0m～1.5m。Samp21中高程為288.48m～320.28m，研究區內有橋梁、公路、小隧道和植被等。samp23中高程為262.27m～348.29m，有較明顯的地形起伏，研究區內有大型建筑物、不規則形狀的建筑物以及汽車等人工地物。samp31中高程為226.94m～343.95m，數據中存在部分極低點，研究區內分布有高低起伏的建筑物、零散汽車以及分布于建筑物間的各種植被。Site2為大范圍數據集，覆蓋區域包含了samp21和samp23，高程范圍為248.1m～482.63m，地形變化較大，地物復雜。

3.2 實驗結果

Site2數據集包含了408921個激光點，覆蓋范圍廣，研究區內地形變化較大，在進行高程統計前，需對該數據集進行分區以保證子區域內沒有“同高異類”現象，然后對分區內數據進行獨立的統計分割，再將分割結果合并，提取地面點，建立初始地形。圖2分別為原始數據集、數據分區結構以及基于本文算法的濾波結果。圖2中的（a）圖是原始數據集按照高程大小進行渲染的顯示結果。從圖中可以看出，數據集的左上部分地形較平坦，建筑物與植被高差較明顯，而數據集的右下部分地形變化劇烈，植被以及建筑物分布較復雜。因此，在圖2（b）的數據分區結構中，右下側被劃分得更為詳細，以確保高程統計方法可以準確地識別地物與地形間的高程差異。圖2（c）為基于高程分割數據集進行三角網漸進濾波后得到的最終結果。該方法有效地剔除了研究區內的建筑物和植被，同時，在高程渲染圖的左上角和左下角顯示不明顯的橋梁也被有效地剔除了。在剔除地物點的同時，此方法有效保留了研究區最右側的高地地形。

圖2 大數據集分區

圖3是3個數據集的ISPRS公布的半自動半手工分類結果、經原始三角網方法和本文方法濾波得到的對比結果圖。

圖3中原始算法與本文算法的濾波效果都比較理想，在3組數據中橋梁、大部分建筑物和零碎植被點都被較好地剔除了。不同之處在于，本文算法在3組數據的邊界建筑物濾波效果優于原始算法。

在samp21中，對比右下角的建筑物，原始三角網算法錯誤地保留了部分建筑物邊界，并且濾波結果中建筑物的邊緣也相對粗糙，本文方法較好地剔除了這部分地物點；同時，對比左下角處的地物部分，原始算法相比較本文算法，過度地剔除了一部分地面點。在samp23中，本文方法較原始三角網算法，更好地剔除了右邊界處的建筑物，建筑物邊沿與真值數據也很匹配，但是左邊界下半部分的少量地物點被錯誤地保留了。在samp31中，本文方法正確地剔除了下邊界左部的植被、下邊界中部的建筑物以及中部的植被建筑物交錯地物點。

圖3 3組數據濾波結果對比圖

3.3 誤差分析

根據ISPRS公布的誤差評定方法，對3組數據集進行誤差統計。表1和表2分別是原始三角網漸進濾波方法與本文方法得到的3個數據集的統計結果。其中，以ISPRS公布的半自動半手動分類結果作為真值判斷依據，a為被正確分類的地面點，b為被錯分類為地物點的地面點，c為被錯誤分類為地面點的地物點，d為被正確分類的地物點；一類誤差為＊100%，二類誤差為＊100%。

在原始算法精度較高的情況下，3組數據經過本文方法濾波后，提取了更多地面點的同時，有效剔除了地物點的影響，精度均得到一定的提高。在samp21和samp31數據集中，研究區地形起伏較小，進行高程統計分割后，先驗知識的準確度也相對較高，兩個數據集的濾波精度提升幅度，相比較地形起伏更加劇烈的samp23數據集更大，所以研究區的地形條件是影響誤差的一個重要因素。所有數據集在先驗知識的基礎上提取了更多的初始地面點，待分類激光點相對減少，同時建立了更加細致的初始地形，所以在三角網漸進濾波時，待分類激光點到地形的參數值更加準確，濾波精度得到提高。

表1 原始三角網漸進算法精度評價

表2 本文算法濾波結果

原始三角網漸進濾波算法中，格網尺寸的劃分會受到研究區內最大建筑物尺寸影響，以免劃分的格網中出現沒有地面點的情況。在本文方法中，經過高程統計分割后，在很大程度上解決了格網尺寸的約束，格網可以按照用戶自定義的尺寸進行劃分。為了測試格網尺寸對濾波結果的影響，本文從5m～40m劃分了8個不同的格網尺寸，然后統計了3組數據在不同尺寸下得到的濾波結果的兩類誤差，圖4為3組數據的誤差變化情況。研究發現，本方法得到的濾波結果，二類誤差相對平穩，一類誤差隨著格網尺寸變化相對劇烈，并且隨著格網尺寸變大，一類誤差總體呈上升趨勢。當格網尺寸設定在10m～20m之間時，一類誤差和二類誤差較小且相對平衡，可以得到比較理想的濾波結果。

圖4 本文算法中兩類誤差隨格網尺寸變化

4 結束語

本文在三角網漸進濾波方法及其改進方案［19－21］基礎上，提出基于高程統計的三角網漸進濾波方法，該方法先從點云數據中提取先驗知識，基于先驗知識再進行三角網濾波，經過3組數據驗證，均得到較好的濾波效果，比較適合于城區點云數據處理。在進行高程統計分割時，當研究區面積較大或者地形起伏較劇烈時，需要對研究區進行分區統計，以保證統計分割的精度。由于有高程統計分割得到的先驗知識，在進行格網尺寸劃分時，該方法可以不受原始三角網漸進濾波算法中最大建筑物尺寸的約束，自由選擇尺寸，經過實驗驗證，此方法以10m～20m的格網尺寸得到的濾波結果更優。與原始三角網漸進濾波方法相比，此方法提高了濾波精度，但對數據進行分區時，本文方法未引入有效的評估參數，自動化水平有待進一步提高。

［1］劉春，陳華云，吳杭斌.激光三維遙感的數據處理與特征提取［M］.北京：科學出版社，2010.

［2］周曉明，馬秋禾，許曉亮，等.LIDAR點云濾波算法分析——以ISPRS測試實驗為參考［J］.測繪工程，2011，20（5）：36－39.

［3］劉經南，許曉東，張曉紅，等.機載激光掃描測高數據分層迭代選權濾波方法及其質量評價［J］.武漢大學學報：信息科學版，2008，30（6）：551－555.

［4］黃先鋒，李卉，王瀟，等.機載LiDAR數據濾波方法評述［J］.測繪學報，2009，38（5）：466－469.

［5］KILIAN J，HAALA N，ENGLICH M.Capture and evaluation of airborne laser scanner data［J］.ISPRS，1996，31（B3）：383－388.

［6］CHANG K J，TABOADA A，CHAN Y C.Geological and morphological study of the Jiufengershan landslide triggered by the Chi－Chi Taiwan earthquake［J］.Geomorphology，2005，71（3－4）：293－309.

［7］梁欣廉，張繼賢，李海濤.一種應用于城市區域的自適應形態學濾波方法［J］.遙感學報，2007，11（2）：276－281.

［8］隋立春，張熠斌，柳艷，等.基于改進的數學形態學算法的LiDAR點云數據濾波［J］.測繪學報，2010，39（4）：390－396.

［9］AXELSSON P.Processing of laser scanner data－algorithms and applications［J］.ISPRS Journal of Photogrammetry and Remote Sensing，1999，54（2－3）：138－147.

［10］AXELSSON P.DEM generation from laser scanner data using adaptive TIN models［J］.International Archives of Photogrammetry and Remote Sensing，2000，33（Part4－B4／1）：110－117.

［11］AXELSSON P.Ground estimation of laser data using adaptive TIN－models［C］.Proceedings of OEEPE Workshop on Airborne Laser Scanning and Interferometric SAR for Detailed Digital Elevation Models，Stockholm，Sweden，March，2001：185－208.

［12］ELMQVIST M.Terrain modeling ＆analysis using laser scanner data［J］.ISPRS，2001，34（3／W4）：219－224.

［13］BROVELLI M A，CANNATA M.Digital terrain model reconstruction in urban areas from airborne laser scanning data：the method and the example of the town of Pavia（Northern Italy）［J］.ISPRS，2002，34（2）：43－48.

［14］蘇偉，孫中平，趙冬玲，等.多級移動曲面擬合LIDAR數據濾波算法［J］.遙感學報，2009，13（5）：833－839.

［15］尚大帥，馬東洋，趙羲，等.一種基于移動曲面擬合的機載LIDAR點云數據濾波方法［J］.測繪技術裝備，2012，14（2）：23－26.

［16］SITHOLE G，VOSSELMAN G.Experimental comparison of filter algorithms for bare－Earth extraction from airborne laser scanning point clouds［J］.ISPRS Journal of Photogrammetry and Remote Sensing，2004，59（1－2）：85－101.

［17］楊曉云，梁鑫.基于參數統計的LiDAR數據分割算法［J］.廣西民族師范學院學報，2012，29（3）：33－35.

［18］莊軍，李弼程，陳剛.一種有效的文本圖像二值化方法［J］.微計算機信息，2005，21（06X）：56－57.

［19］李卉，李德仁，黃先鋒，等.一種漸進加密三角網LIDAR點云濾波的改進算法［J］.測繪科學，2009，34（3）：39－41.

［20］隋立春，張熠斌，張碩，等.基于漸進三角網的機載LiDAR點云數據濾波［J］.武漢大學學報：信息科學版，2011，36（10）：1159－1163.

［21］左志權，張祖勛，張劍清.知識引導下的城區LiDAR點云高精度三角網漸進濾波方法［J］.測繪學報，2012，41（2）：246－251.