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

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

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

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

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

1 前言

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

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

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

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

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

2.1 基于離散度的初始激光點(diǎn)粗分類

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

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

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

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

圖1 離散度計(jì)算

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

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

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

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

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

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

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

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

圖3 全局高程閾值的確定

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

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

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

2.3 后向地形點(diǎn)加密

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

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

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

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

3 實(shí)驗(yàn)結(jié)果與分析

實(shí)驗(yàn)1 ISPRS測(cè)試數(shù)據(jù)

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

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

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

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

實(shí)驗(yàn)2 三個(gè)不同國(guó)家地區(qū)數(shù)據(jù)實(shí)驗(yàn)

圖6 三個(gè)不同國(guó)家地區(qū)數(shù)據(jù)實(shí)驗(yàn)

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

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

4 結(jié)語(yǔ)

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