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

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

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

3.加拿大滑鐵盧大學(xué)地理與環(huán)境學(xué)院，加拿大滑鐵盧，N2L2R7

4.武漢大學(xué)遙感信息工程學(xué)院，武漢430079

1 前言

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

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

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

2 復(fù)雜城區(qū)激光點云濾波算法

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

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

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

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

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

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

圖1 離散度計算

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

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

2.2 基于區(qū)域的精分類

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

圖2 基于區(qū)域分類濾波流程圖

（1）建立Delaunay三角網(wǎng)數(shù)據(jù)結(jié)構(gòu)。其目的是用來尋找臨近點，形成一個個聯(lián)通區(qū)域。

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

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

圖3 全局高程閾值的確定

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

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

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

2.3 后向地形點加密

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

其中，如果待定激光點Lp同時滿足這兩個判據(jù)，則Lp被認(rèn)定為地形點。

圖4 判據(jù)計算示意圖

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

3 實驗結(jié)果與分析

實驗1 ISPRS測試數(shù)據(jù)

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

圖5 ISPRS sample數(shù)據(jù)

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

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

實驗2 三個不同國家地區(qū)數(shù)據(jù)實驗

圖6 三個不同國家地區(qū)數(shù)據(jù)實驗

表1 濾波數(shù)據(jù)的I型誤差和II型誤差分析表

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

4 結(jié)語

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

[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.