陳琳，范湘濤，杜小平

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

1 引 言

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

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

2 方 法

2.1 基于高程統(tǒng)計的閾值分割

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

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

①繪制高程直方圖。點云數(shù)據(jù)的高程直方圖以高程為橫坐標，頻數(shù)為縱坐標，反映各高程與其出現(xiàn)頻數(shù)之間的關(guān)系。

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

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

2.2 三角網(wǎng)漸進濾波

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

①參數(shù)設(shè)置。參數(shù)包括：研究區(qū)最大建筑物尺寸、點到三角面的最大高程閾值和點到三角面節(jié)點的最大角度閾值。

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

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

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

2.3 基于高程統(tǒng)計的三角網(wǎng)漸進濾波

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

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

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

該方法濾波主要流程：

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

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

③對分割得到的備選地面點劃分格網(wǎng)，提取每個格網(wǎng)內(nèi)的最低點作為種子點，創(chuàng)建初始三角網(wǎng)地形。

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

圖1 基于高程統(tǒng)計的三角網(wǎng)漸進濾波流程圖

3 數(shù)據(jù)與實驗

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

3.1 實驗數(shù)據(jù)

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

3.2 實驗結(jié)果

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

圖2 大數(shù)據(jù)集分區(qū)

圖3是3個數(shù)據(jù)集的ISPRS公布的半自動半手工分類結(jié)果、經(jīng)原始三角網(wǎng)方法和本文方法濾波得到的對比結(jié)果圖。

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

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

圖3 3組數(shù)據(jù)濾波結(jié)果對比圖

3.3 誤差分析

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

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

表1 原始三角網(wǎng)漸進算法精度評價

表2 本文算法濾波結(jié)果

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

圖4 本文算法中兩類誤差隨格網(wǎng)尺寸變化

4 結(jié)束語

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

［1］劉春，陳華云，吳杭斌.激光三維遙感的數(shù)據(jù)處理與特征提?。跰］.北京：科學出版社，2010.

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

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

［4］黃先鋒，李卉，王瀟，等.機載LiDAR數(shù)據(jù)濾波方法評述［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］梁欣廉，張繼賢，李海濤.一種應用于城市區(qū)域的自適應形態(tài)學濾波方法［J］.遙感學報，2007，11（2）：276－281.

［8］隋立春，張熠斌，柳艷，等.基于改進的數(shù)學形態(tài)學算法的LiDAR點云數(shù)據(jù)濾波［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數(shù)據(jù)濾波算法［J］.遙感學報，2009，13（5）：833－839.

［15］尚大帥，馬東洋，趙羲，等.一種基于移動曲面擬合的機載LIDAR點云數(shù)據(jù)濾波方法［J］.測繪技術(shù)裝備，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］楊曉云，梁鑫.基于參數(shù)統(tǒng)計的LiDAR數(shù)據(jù)分割算法［J］.廣西民族師范學院學報，2012，29（3）：33－35.

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

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

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

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