殷 飛 齊 華,2 薛曉濱 盧建康

(1.西南交通大學(xué)信息科學(xué)與技術(shù)學(xué)院,四川成都 610031；2.西南交通大學(xué)土木工程學(xué)院,四川成都 610031；

3.中鐵二院工程集團有限責(zé)任公司信息中心,四川成都 610031；

4.中鐵二院工程集團有限責(zé)任公司測繪分院,四川成都 610031)

機載激光雷達(dá)測量(Light Detection and Ranging)技術(shù)是近些年快速發(fā)展起來的一種空間信息獲取技術(shù),能夠快速高效地獲取數(shù)字地表三維模型。目前,LiDAR系統(tǒng)硬件技術(shù)發(fā)展得比較成熟,大部分屬于硬件及系統(tǒng)集成方面的問題已得到解決。但數(shù)據(jù)后處理技術(shù)卻發(fā)展相對滯后[1],主要包括濾波分類、地物提取、三維重建及海量點云數(shù)據(jù)簡化?！笆晃濉逼陂g,國家863計劃地球觀測與導(dǎo)航技術(shù)領(lǐng)域也開始進行LiDAR數(shù)據(jù)后處理方面的研究[7]。

機載激光雷達(dá)系統(tǒng)采集到的原始數(shù)據(jù)點離散分布于三維空間,這些點可能打到真實地表面,也有可能打到植被或人工建筑物等地物表面。為了從點云數(shù)據(jù)中得到數(shù)字高程模型,需要分離地面點和地物點,即Li-DAR數(shù)據(jù)濾波。作為LiDAR數(shù)據(jù)后處理最為重要的組成部分之一,濾波結(jié)果直接影響所生成DEM(Digital Elevation Model)的質(zhì)量,也是進行地物提取與三維重建的必要前提。

如何實現(xiàn)由LiDAR數(shù)據(jù)快速、準(zhǔn)確地提取出地面點,一直是國內(nèi)外專家學(xué)者們研究的重點、難點和熱點之一。出于商業(yè)方面的考慮,現(xiàn)有的幾種商業(yè)軟件其數(shù)據(jù)處理算法都不對用戶公開。

1 濾波原理及方法回顧

LiDAR數(shù)據(jù)濾波的原理基于高程突變,即鄰近的兩個激光點間高程變化較大時,一般認(rèn)為不是由于地形的陡然起伏所造成的,更為可能是由于地面存在地物而引起,高程較大的點很可能是地物點。若兩點間的距離越近,高差越大,較高的點更有可能為地物點。故判斷某點是否地面點時,要考慮該點到參考地面點的距離。隨著兩點間距離的增加,其高差也會變大[6]。

迄今為止,國內(nèi)外的專家學(xué)者對濾波算法進行了深入研究。最簡單的濾波算法認(rèn)為一定區(qū)域內(nèi)高程最低的點即為地面點；該算法基于平坦地形,誤差較大。德國斯圖加特大學(xué)的Lindenberger提出了數(shù)學(xué)形態(tài)學(xué)濾波算法,對LiDAR數(shù)據(jù)進行形態(tài)開運算并利用自回歸過程改善開運算的結(jié)果[4]。奧地利維也納大學(xué)的Kraus和Pfeifer提出了迭代線性最小二乘內(nèi)插濾波算法,根據(jù)計算得到的殘差值給每一個激光點賦權(quán)值,由此分離出地面點[3]。Vosselman提出了基于坡度的濾波算法,根據(jù)最大坡度 γ和垂直測量精度 σz確定高差閾值函數(shù)[6]。Axelsson提出了三角網(wǎng)迭代加密濾波算法,通過不斷向三角面加入滿足高差閾值和角度閾值的激光點,達(dá)到重建地表面的目的[2]。武漢大學(xué)張小紅提出了移動曲面擬合濾波算法,用二次曲面去擬合空間曲面,擬合值與原始值間高差小于閾值則視為地面點[8]。

目前還沒有哪一種濾波算法對任意地形的處理結(jié)果都很完美。Sithole通過實驗分析了8種典型濾波算法,指出各種算法都存在一定的問題[5]。本文通過研究三角網(wǎng)迭代加密濾波算法,提出了一種改進方法。

2 濾波算法

數(shù)據(jù)處理過程中,若將原始離散激光點云內(nèi)插為規(guī)則格網(wǎng)DSM(Digital Surface Model)格式,雖然可以引入成熟的圖像處理算法,但勢必會引入誤差,另外數(shù)據(jù)的精度也會降低。在原始LiDAR數(shù)據(jù)中,相同的平面坐標(biāo)可以對應(yīng)多個不同高程值的激光點。以規(guī)則格網(wǎng)DSM表示時,一對X、Y坐標(biāo)只能對應(yīng)一個高程值。因此,認(rèn)為規(guī)則格網(wǎng)表示方式會在一定程度上丟失數(shù)據(jù)原有的某些屬性,本文選擇對基于離散點的三角網(wǎng)迭代加密法進行改進。

Axelsson提出的三角網(wǎng)迭代加密濾波算法思想：首先對測區(qū)進行分塊處理(分塊大小以測區(qū)地物最大長度為準(zhǔn)),初始點選擇地面上一些高程較低的點,然后選取每塊中的最低點作為初始地面點并建立不規(guī)則三角網(wǎng)模型。定位待判斷點P對應(yīng)的三角形,P及其在三角面上的投影點與三角形頂點間的角度分別記為α、β、γ,P到三角面的垂直距離記為d。若 α、β、γ和 d都小于設(shè)定的閾值,則接受該點為地面點并重新建立三角網(wǎng)(見圖1)。接著判斷下一個點,直到所有的點都被分為地面點或非地面點。通過不斷地加入滿足條件的激光點來擴大地面模型,每加入一個新點都更新不規(guī)則三角網(wǎng),最終得到非常接近地表面的三角網(wǎng)模型[2]。

圖1 濾波算法參數(shù)示意

經(jīng)研究,該算法存在的主要缺點有：

①該算法假定一定區(qū)域內(nèi)高程最低的點為地面點,此時如果存在低于地面的點,則會錯誤分類到地面點類,因此需要先分離低點。此外,沒有考慮空中噪點對濾波算法的影響。

②參數(shù)值不能隨著地形自適應(yīng)改變。在絕對理想的平地區(qū)域,各激光點的高差為0,若此時的高差閾值還取一固定值,不可避免錯誤地將地物點分類為地面點。此外,陡坡邊緣的點可能被錯誤分類為地物點。

針對上述問題,本文對該算法作了以下改進,算法流程見圖2。

圖2 算法流程

(1)對原始激光點云數(shù)據(jù)進行預(yù)處理,剔除錯誤點,包括低點、噪點等。

剔除單個低點：

若某一個點為低點,則其高程值與領(lǐng)域內(nèi)其他點的高差都應(yīng)滿足一定的高差閾值。設(shè)待判斷點為Pa,以Pa為中心點,搜索一定鄰域范圍內(nèi)的其他激光點,若對鄰域內(nèi)任意激光點P都有ZPa-ZP＜-H,H為設(shè)置的高差閾值,則可判斷Pa為低點。用數(shù)學(xué)形式表達(dá)為：對于激光點云集P,若其鄰域內(nèi)任一點滿足：PT={Pa∈P|?Pr∈Pn：ZPa-ZP＜-H},則判斷 Pa為低點 。

剔除低點簇：

若某點的鄰域內(nèi)幾個連續(xù)點均為低點,則其高程與內(nèi)插高程之差應(yīng)滿足一定的高差閾值。設(shè)低點簇中的某個待判斷點為Pa,以Pa為中心點,搜索一定鄰域范圍內(nèi)的其他激光點并進行插值,其內(nèi)插高程為Pi,若ZPa-ZPi＜-H,H為設(shè)置的高差閾值,則可判斷Pa為低點(見圖2)。該方法也可用于剔除單個低點。

本文提到的噪點專指空中點,噪點不僅會影響數(shù)據(jù)瀏覽,還會影響濾波算法的使用,所以應(yīng)先剔除這些點。由于這些點的高程與周圍點高程相比突變較大,一般可以通過設(shè)置絕對高程閾值去除。

(2)對測區(qū)分塊并建立虛擬格網(wǎng)索引,取每一塊中高程最低的點PL作為種子點。查找塊的凸包并將PL的高程賦給邊界點,由凸包與PL建立初始三角網(wǎng),此時所有三角形位于同一平面。根據(jù)測區(qū)情況設(shè)置初始高差閾值。

(3)分塊時保證各塊間有一定程度的重疊。

(4)高差閾值設(shè)置為自適應(yīng)函數(shù)。

式中 α——待判斷點對應(yīng)三角形的坡度；

d——激光點的平均間距；

σz——輸入的常數(shù)。

3 濾波實驗及分析

3.1 濾波實驗

為了有針對性地比較濾波算法,特地選擇了兩塊典型區(qū)域的數(shù)據(jù)進行濾波實驗。

實驗數(shù)據(jù)1為某城區(qū)實測LiDAR數(shù)據(jù),數(shù)據(jù)區(qū)域為長方形,測區(qū)面積為250 m ×250 m,原始文件5.68 M,點數(shù)為213093個,平均間距約0.77 m,密度約3.4個/m2。原始數(shù)據(jù)中包含有低點,該區(qū)域地形平坦,建筑物較多且復(fù)雜,地物主要為植被、汽車、房屋、橋梁等,見圖 3(a)。

實驗數(shù)據(jù)2為某森林區(qū)域?qū)崪yLiDAR數(shù)據(jù),測區(qū)面積為146 m ×80 m,原始文件1.57M,點數(shù)為58 894個,平均間距約0.46 m,密度約5個/m2。該區(qū)域地形平坦,地物主要為植被,見圖3(b)。

由圖4(a)、圖4(b)可看出,改進后算法剔除了錯誤點、小汽車等近地點。同時有少數(shù)面積較小的房屋未過濾,形成小的突起,房屋邊緣的過濾效果不是很理想。此外有部分道路點被錯誤分類為地物。

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

圖4 城區(qū)濾波結(jié)果暈渲圖

圖5 森林區(qū)域濾波結(jié)果暈渲圖

由圖5(a)、圖5(b)可看出,森林區(qū)域采用改進后算法處理的結(jié)果過濾了大部分植被,但仍有少數(shù)極低的灌木點錯誤分類為地面點,同時將部分地面點錯誤分類了地物點。

3.2 濾波難點分析

機載激光雷達(dá)數(shù)據(jù)濾波結(jié)果與激光點的密度、測區(qū)地形和地物的復(fù)雜程度都有很大關(guān)系。一般情況下,大多數(shù)濾波算法都能夠處理地形、地物不復(fù)雜的區(qū)域,如平原、城市、稀疏的樹林等。從工程實踐可知：困難區(qū)域的處理是一個難點。因此,能否準(zhǔn)確對任意復(fù)雜地形、地物進行濾波處理是評價濾波算法優(yōu)劣的一個重要指標(biāo)。為了定性分析濾波算法,特選取幾處濾波難點,包括錯誤點、復(fù)雜地物、附著地物、植被(低矮植被和坡面上的植被)等(如表1)。

3.3 定量分析

ISPRS組織第Ⅲ委員會第三工作組于2004年發(fā)布了濾波方法比較報告,文章中將濾波誤差分為兩類。第Ⅰ類誤差：將地面點錯誤分類為地物點；第Ⅱ類誤差：錯誤接收地物點并歸入地面點[5]。

表1 濾波實驗定性比較

第Ⅰ類誤差：

第Ⅱ類誤差：

總誤差：

式中 a——正確標(biāo)識為地面點的數(shù)目；

b——錯誤分類為地物點的數(shù)目(由于第Ⅰ類誤差產(chǎn)生)；

c——錯誤分類為地面點的數(shù)目(由于第Ⅱ類錯誤產(chǎn)生)；

d——正確分類為地面點的數(shù)目。濾波算法考慮最小化第Ⅱ類誤差,即應(yīng)剔除盡量多的地物點,當(dāng)?shù)匚锖苄』蚝苜N近地表面時,甚至可以剔除部分真實地面點。其原因是第Ⅱ類誤差對DEM質(zhì)量的影響比第Ⅰ類誤差要大(如表2、表3所示)。

表2 數(shù)據(jù)1濾波結(jié)果定量比較

三角網(wǎng)迭代加密算法將部分低矮植被點和近地點錯誤分類為地面點,導(dǎo)致第Ⅱ類誤差較大。改進后采用坡度閾值函數(shù),有效剔除了低植被點,處理結(jié)果的第Ⅱ類誤差減小,總誤差有一定程度減小。同時由于剔除了部分真實地面點,導(dǎo)致第Ⅰ類誤差變大。

表3 數(shù)據(jù)2濾波結(jié)果定量比較

4 結(jié)論

通過研究三角網(wǎng)迭代加密濾波算法,提出了一種改進方法。在濾波處理前先剔除錯誤點,以降低其對后續(xù)算法的影響。針對不同地形情況的參數(shù)設(shè)置是一個難點,提出了根據(jù)三角形坡度建立自適應(yīng)閾值函數(shù)的方法,并選擇兩塊典型區(qū)域數(shù)據(jù)驗證了改進后濾波算法的可行性。實驗證明：該方法能夠有效處理Li-DAR數(shù)據(jù),剔除錯誤點和近地點并得到DEM數(shù)據(jù),且第Ⅱ類誤差、總誤差相對減小。但是過濾小面積房屋的效果不是很好,可以考慮結(jié)合回波信息進行處理。此外,由于LiDAR數(shù)據(jù)是海量數(shù)據(jù),本算法中又涉及動態(tài)構(gòu)建不規(guī)則三角網(wǎng),其時間復(fù)雜度較高。因此,選擇速度較快的構(gòu)網(wǎng)方法是降低濾波時間復(fù)雜度的關(guān)鍵。

[1]Axelsson P.Processing of laser scanner data-algorithmsand applications.ISPRSJournal of Photogrammetry and Remote Sensing,1999,54：138-147

[2]Axelsson.DEM generation from laser scanner data using adaptive TIN models.International Archives of Photogrammetry and Remote Sensing,2000,33(B4)：110-117

[3]Kraus K,Pfeifer N.Determination of terrain modelsin wooded areas with airborne laser scanner data.ISPRS Journal of Photogrammetry and Remote Sensing,1998,53：193-203

[4]Lindenberger J.Laser-profilmessungen zur topographischen gelandeaufnahme. Verlag der Bayerischen Akademie der Wissenschaften,1993

[5]Sithole G,Vosselman,G.Experimental comparison of filter algorithms for bare-Earth extraction from airborne laser scanning point clouds.ISPRS Journal of Photogrammetry and Remote Sensing,2004,59：85-101

[6]Vosselman G.Slope based filtering of laser altimetry data.IAPRS,2000,33(B3/2)：935-942

[7]科技部.863計劃地球觀測與導(dǎo)航技術(shù)領(lǐng)域2006年度專題申請指南[R].北京：科技部,2006

[8]張小紅.機載激光雷達(dá)測量技術(shù)理論與方法[M].武漢：武漢大學(xué)出版社,2007