基于多尺度的城區機載點云濾波方法研究

鐘良，湯璇，管海燕，鄔建偉

1.長江委長江空間信息技術工程有限公司，武漢430010

2.武漢軍械士官學校計算機系，武漢430075

3.加拿大滑鐵盧大學地理與環境學院，加拿大滑鐵盧，N2L2R7

4.武漢大學遙感信息工程學院，武漢430079

1 前言

機載Lidar點云數據為三維空間中呈現隨機分布的點云。在點云中，有些點是真實地形點，有些則是人工地物（比如建筑物、橋梁、塔、車輛等）或者是自然地物（如樹木、灌木等）。從激光掃描點云中區分地形點子集與非地形點子集（包含人工地物與自然地物）稱之為濾波。目前，自動區別地形點和地物點一直是研究的熱點，特別是針對場景中含有很多各種復雜地形地貌特征情況。目前學術界已經提出了很多濾波方法：數學形態學[1-4]、線性預測[5-6]、漸進加密[7-8]和分割[9-11]，其他濾波方法還有全波形[12]、機器學習[13-14]、掃描線[15-16]。

以上算法都是基于對地面點和非地面點的理解概念或者思路來設計的，從算法設計角度來說可分為直接應用于不規則分布的點集，或者將原始點集重采樣為規則格網再利用成熟的圖像處理技術進行處理；從算法對激光點的幾何關系的分析上可分為考察單點與單點之間的關系，單點與點集或者點集與點集之間的關系。從算法對分析點云數據特征中某種不連續測度上可分為高程差、坡度、到結構表面的最小距離或到參數平面的最小距離等。總的來說上述的各種濾波算法都基于各種不同的假設前提，如坡度、平面、表面和分割等；有些算法的處理方式是一次單步完成，有些是通過迭代處理，這也意味著它們分別在準確性與速度方面的存在著或多或少的劣勢[17]。

考慮到對復雜城區點云數據進行濾波、生成DTM，針對城區建筑物偏多，局部地形較為平坦的特性，本文提出了一種以離散度分析為基礎的平面點提取與基于強度方差的區域增長相結合的組合濾波算法。算法分為三步：首先計算每個點的局部空間分布特性，將激光點大致分為平面點，邊緣點和離散點；其次將平面點作為分析對象跟蹤出一組聯通的區域，利用開運算和閉運算去除零散的小區域，基于強度方差測度將其分類為地面，非地面以及未確定；最后后向搜索從邊緣點以及未確定點中加密地形點。

2 復雜城區激光點云濾波算法

本算法分為三個步驟：基于離散度的初始激光點粗分類，基于區域的精分類，以及后向地形點加密。算法首先分析并計算每個點與其周圍一定鄰域內激光點之間的幾何特征值關系，將點云粗分類為平面點、邊緣點和離散點；然后對平面點進行區域跟蹤，利用高程以及強度方差等將平面點分類為地面點、建筑物點以及為確定類別點；最后對地面點構建不規則三角網，反向分析確定類別點以及邊緣點用以加密地面點集。

2.1 基于離散度的初始激光點粗分類

激光點云空間離散度可以通過計算每個點與其周圍一定鄰域范圍內所有點的空間幾何分布來獲取。一般樹木、小型地物（如汽車等）等在空間各個方向上都較為分散，而裸露地面以及建筑物在空間的離散程度可以認為沿二維表面（不一定水平）分布，這種離散性可以通過離散矩陣的特征值來分析。具體方法：首先構建激光點云的二維K-D樹，搜索激光點云中每個點鄰域內的全部相鄰點，建立該激光點在空間上的3×3離散矩陣，見公式（1）：

其中，Sj為第j個點的3×3離散矩陣，n為第j個點鄰域內相鄰點數目，vi為第j個點的第i相鄰點的空間坐標vi=(xi，yi，zi)，為vi的轉置矩陣，M為激光點云點數。

將離散矩陣Sj作奇異值分解，可獲取該點矩陣的三個特征值λ0、λ1、λ2，并將特征值從大到小排列，e0、e1、e2為三個特征值對應的特征矢量[18]。則公式（1）可以表示為：

根據特征值，一般設定三個類別：（1）λ2＜＜λ1≈λ0，則該點被標記為平面類（圖1（a））；（2）λ2≈λ1＜＜λ0，則該點被標記為邊緣類（圖1（b））；（3）λ0≈λ1≈λ2，并且都足夠大，則標記為空間離散類（圖1（c））。

圖1 離散度計算

根據以上的描述，無法確定閾值來判定這三個類別。因此公式（3）可表示為：

在(Slinear，Splanar，SSphere)中，最大的值就為該Sj所屬的類別。一般大部分房屋頂由平面組成，因此呈現出平面特征，但是其中部分由于屋頂結構復雜、形狀多變，因此會有些屋頂出現少量離散的點，且地面上小物體（如小汽車等）也會造成少量離散點；基本上裸露地面會出現大量的平面點特性；而由于城市中樹木稀疏，因此各向離散性能夠很好地表現出來，這樣就可以去除大部分植被點以及建筑物邊緣點等。經過離散度計算，平面點作為下步處理的輸入數據來獲取地形點。

2.2 基于區域的精分類

經過離散度計算后，Lidar點集中保留有地面點和建筑物點，以及其他零散的激光點，為此需要對區域進行精分類以將地面點從建筑物點中分離出來。由于地面是由連續聯通的面片構成，因此本文不是處理單個的點，而是將區域作為分析單元。如圖2所示為本階段的流程圖。

圖2 基于區域分類濾波流程圖

（1）建立Delaunay三角網數據結構。其目的是用來尋找臨近點，形成一個個聯通區域。

（2）開運算與閉運算的運用。由于地形和建筑物也是具有一定面積的連續平面，開運算和閉運算可用于去除那些激光點很少的小區域。一般這個點數量閾值Tφ需要根據點密度來參考。在本研究中，點數量閾值Tφ設為20，如果區域內點的數量小于20個點，這個區域將被歸并到邊緣點中。

（3）計算全局高程閾值和局部高程閾值。全局高程閾值確定HG是根據輸入數據的高程直方圖分析獲取的。如圖3所示，對建筑物和裸露地面，在高程直方圖中會出現多個波峰，一般第一個波谷就是區別二者的全局高程閾值（圖中大約為150m）。局部高程閾值的確定HL是通過計算每個待分類區域的最大范圍W以及中心點坐標。然后基于中心點，在原始點云獲取一個2W的區域。則2W區域內高程平均值被作為局部高程閾值HL。

圖3 全局高程閾值的確定

（4）如果待判定區域，其平局高程均大于全局閾值和局部閾值，則該區域被認為是建筑物平面；區域平均高程均小于全局和局部閾值，則被分類為地面點；剩下的不符合這兩個判斷條件的，將進行進一步的分析。

（5）根據上一個步驟所獲取的地面點和建筑物類別點，再分別計算已分類地面點和建筑物點的強度方差，其方差計算公式如下：

通過以上三個步驟，平面點基本分類為地面點，建筑物點和待定點。地面點則被用于建立Delaunay三角網作為初始地形。待定點和前面分類為邊緣點將被帶入下一個步驟，進行后向地面點加密獲取。因為邊緣點一般是處于斷裂線附近，其中有部分地面點。

2.3 后向地形點加密

對于任意一個Lidar點如果同時滿足兩個地形判定條件，則將其作為地形點。假設N1～N3為已判定的地形點，Lp為待判定激光點，d為待判定激光點Lp到N1～N3所構成的三角面片的距離，(α，β，λ)分別為Lp到N1～N3構成的三角面片的夾角。給定兩個閾值：距離閾值dmax和角度閾值θmax，則兩類判據分別定義為[7]：

其中，如果待定激光點Lp同時滿足這兩個判據，則Lp被認定為地形點。

圖4 判據計算示意圖

將初始地面點集構建不規則三角網，從邊緣點集中需找潛在的地面點加入到地面點集中。遍歷不規則三角網，計算落在三角形中每個激光點到三角形的距離以及三個夾角，如果滿足給定的高程差閾值以及角度閾值則認為該點是地面點。在TIN中所有三角形遍歷結束后，將獲得的地面點重新加入不規則三角網。重復迭代執行，直至滿足給定的迭代最大次數或者迭代搜索獲取的地面點小于給定最少地面點，迭代結束。

3 實驗結果與分析

實驗1 ISPRS測試數據

本文從ISPRS網站提供的實驗數據中選擇了四個城區測試數據集中的幾個典型場景（數據來自http：//www.commission3.isprs.org/wg3/index.html），如圖5所示。其中點云的水平分辨率均為0.67個點/m2（點間距為1.0～1.5m）。

圖5 ISPRS sample數據

對激光掃描數據進行濾波的主要目標是生成DTM，因此濾波結果中存在的兩種誤差將影響DTM的生成質量。第一種誤差稱為I型誤差，指的是將地形點錯誤的分類為非地形點；第二種誤差稱為I型誤差，指的是將非地形點錯誤的分類為地形點。I型誤差造成精度上的損失，II型誤差將導致DTM質量的惡化。因此I型誤差和II型誤差也成為評價濾波算法的主要指標。根據ISPRS提供的測試數據，表1為對本文濾波結果的I型誤差和II型誤差。

根據表1，其濾波結果還是相當滿意的，尤其是Type II誤差較小，基本上保持了地形的特征。數據Sample11、Sam ple22和Sample24的I型誤差較高，產生的主要原因是地形的坡度變化比較大。由于算法設計主要是針對城市地區，所有濾波的假設條件就使得算法對這類地形具有一定的缺陷。Sam ple22、Sample23的II型誤差相對較高，主要原因是一些建筑物和植被小于算法設定的閾值，因此將一部分點錯誤的分類為地面點了。

實驗2 三個不同國家地區數據實驗

圖6 三個不同國家地區數據實驗

表1 濾波數據的I型誤差和II型誤差分析表

實驗2選取了三塊數據進行實驗，數據一是山西太原部分城區激光點云，如圖6（a）所示；數據二是加拿大多倫多城區數據，如圖6（b）所示；數據三是德國Toposys公司提供的Mannheim城區數據（點平均密度為4 points/m2），如圖6（c）所示。這三塊原始點云數據是根據高程進行分色顯示的，因此從圖中可以看出這三塊地區地勢較為平坦，是典型的城區地貌。圖6（d）、（e）和（f）是三角網加密迭代濾波實驗結果。由于沒有實地的DEM數據做參考，難以定量地衡量其提取的DTM的質量，但是總體感覺上提取的質量還是不錯的，基本上去除了非地形點，保留了地形特征。圖6（f）顯示的是將Mannheim城區激光點云濾波實驗結果投影到其配準好的航空影像上，圖中綠色像素是投到二維影像上的三維激光點。從實驗結果可以看出，基本上將非地形點去除掉，保留了基本地形點。

4 結語

通過上述的實驗比較與分析可知，本文提出的多尺度從粗到細的濾波算法不僅考慮了地形的連續性，也最大顧及了整個掃描區域地形的特征，其算法性能也能保證DTM的精度和質量。這種基于點與其周圍領域內其他點空間關系來計算離散度能夠很好地反映點的特性，可以去除大量的離散點和邊緣點；同樣通過基于區域的濾波精分類，很大程度上減少了候選地形點中的非地形點，提高DTM的質量，減少了II型誤差。由于邊緣點中包含地面點和非地面點，后向加密地形點可以發現更多的地面點，從而提高濾波的質量。實驗證明利用多尺度的濾波算法能夠在大面積、大數據量、復雜城區得到運用，為下一步Lidar數據的應用打下了很好的基礎。

[1]Lindenberger J.Laser-Profilmenssungen zur topographischen Gelandeaufnahme[D].Munich：Deutsche Geod?tische Kommission，1993.

[2]Kilian J，Haala N.Capture and evaluation of airborne laser scanner data[J].International Archives of Photogrammetry and Remote Sensing，1996，31：383-388.

[3]Vosselman G.Slope based filtering of laser altimetry data[C]//Proceedings of International A rchives of Photogrammetry，Remote Sensing and Spatial Information Sciences，WG III/3，Amsterdam，2000.

[4]Silván-Cárdenas J L，Wang L.A multi-resolution approach for filtering Lidar altimetry data[J].ISPRS Journal of Photogrammetry and Remote Sensing，2000，61（1）：11-22.

[5]K raus K，Pfeifer N.Determination of terrain models in wooded areas with airborne laserscanner data[J].ISPRS Journal of Photogrammetry and Remote Sensing，1998，53：193-203.

[6]Kobler A，Pfeifer N.Repetitive interpolation：a robust algorithm for DTM generation from Aerial laser scanner data in forested terrain[J].Remote Sensing of Environment，2007，108：9-23.

[7]Axelsson P.DEM generation from laser scanner data using adaptive TIN models[C]//Proceedings of International Archives of the Photogrammetry，Remote Sensing and Spatial Information Sciences XXXIII（B4/1），1999：110-117.

[8]Sohn G，Dowman I.Terrain surface reconstruction by the use of tetrahedron model with the MDL criterion[C]//Proceedings of the Photogrammetric Computer Vision Symposium，2002.

[9]Voegtle T，Steinle E.On the quality of object classification and automated building modelling based on laserscanning data[C]//Proceedings of IAPRSIS XXXIV 3W 13，2003：149-155.

[10]Tóvári D，Pfeifer N.Segmentation based robust interpolation-a new approach to laser data filtering[C]//Proceedings of ISPRSWG III/3，III/4，V/3 Workshop“Laser Scanning 2005”，Enschede，the Netherlands，2005.

[11]Forlani G.Complete classification of raw LIDAR data and 3D reconstruction of building[J].Pattern Analysis and Applications，2006，8（4）：357-374.

[12]Doneus M，Briese C.Digital terrain model ling for archaeological interpretation within forested areas using full waveform laserscanning[C]//Proceeding of the 7th International Symposium on Virtual Reality，A rchaeology and Cultural Heritage VAST，2006.

[13]Lodha S，K reps E J，Helmbold D，et al.Aerial lidar data classification using Support Vector Machines（SVM）[C]//Proceedings of the 3rd International Symposium on 3D Data Processing，Visualization，and Transmission（3DPVT’06），2006：567-574.

[14]Lu W L，Murphy K P，Little J J，et al.A hybrid conditional random field for estimating the underlying ground surface from airborne lidar data[J].IEEE Transactions on Geoscience and Remote Sensing，2009，47（8）：2913-2922.

[15]Shan J，Sampath A.Urban DEM generation from raw Lidar data：a labeling algorithm and its performance[J].Photogrammetric Engineering&Remote Sensing，2005，71（2）：217-226.

[16]M eng X L，Wang L.A multi-directional ground filtering algorithm for airborne lidar[J].ISPRS Journal of Photogrammetry and Remote Sensing，2009，64：117-124.

[17]Sithole G.Experimental comparison of filtering algorithm s for bare-earth extraction from airborne laser scanning point clouds[J].ISPRS Journal of Photogrammetry and Remote Sensing，2004，59（1/2）：85-101.

[18]Poullis C，You S.Delineation and geometric modeling of road networks[J].ISPRS Journal of Photogrammetry and Remote Sensing，2010，65：165-181.