基于不同濾波方法的機載激光雷達(dá)數(shù)據(jù)處理及表面模型構(gòu)建1)

劉 芳 張 瓊 范文義 陳 成 楊樹文 李 典

(東北林業(yè)大學(xué)，哈爾濱，150040)

激光雷達(dá)(LIDAR)是一種通過位置、距離、角度等觀測數(shù)據(jù)直接獲取對象表面點三維坐標(biāo)、回波和強度信息，實現(xiàn)地表信息提取和三維場景重建的對地觀測技術(shù)［1］。激光雷達(dá)技術(shù)在對地探測方面有強大的優(yōu)勢，具有空間與時間分辨率高、動態(tài)探測范圍大、能夠部分穿越樹林遮擋、直接獲取真實地表的高精度三維信息等特點。這使得激光雷達(dá)可以快速高精度地獲取地表物體的水平和垂直結(jié)構(gòu)信息。機載激光雷達(dá)作為一種新興的獲取信息的技術(shù)，因其獨特的優(yōu)勢而被廣泛應(yīng)用于各個不同的領(lǐng)域。

激光雷達(dá)數(shù)據(jù)處理是獲取信息的前提，而激光雷達(dá)點云數(shù)據(jù)的濾波和分類則是數(shù)據(jù)后處理的主要任務(wù)之一。利用激光雷達(dá)點云數(shù)據(jù)可以獲取反映地形起伏變化的數(shù)字高程模型(DEM)和目標(biāo)區(qū)域的表面三維信息——數(shù)字表面模型(DSM)以及數(shù)字城市模型等。從原始點云數(shù)據(jù)中獲取DEM和DSM這兩種表面模型最關(guān)鍵的步驟是如何有效地、可靠地濾除噪聲點，并分別提取出僅反映地形表面信息的地面點和只反映地物表面信息的地物點，這個過程即濾波分類。筆者綜合考慮研究區(qū)特有的地形、植被等因素，對涼水國家級自然保護(hù)區(qū)不同類型的區(qū)域采用不同的算法進(jìn)行濾波試驗，取得了較好的效果。

1 研究區(qū)概況及數(shù)據(jù)獲取

研究區(qū)域為黑龍江省涼水國家級自然保護(hù)區(qū)，位于北緯 47°6'49″～ 47°16'10″，東經(jīng) 128°47'8″～128°57'19″。保護(hù)區(qū)全境為山地，總面積6394 hm2，屬小興安嶺南部達(dá)里帶嶺支脈東坡。

LIDAR數(shù)據(jù)于2009年9月采用LiteMapper 5600激光雷達(dá)系統(tǒng)獲得。激光掃描儀為Riegl LMSQ560，其波長為1550 nm，回波寬度分辨率0.15 m，垂直精度可達(dá)0.15 m。采樣點云密度大于2點/m2，飛行高度達(dá)800 m，平均對地飛行速度是180 km/h。此外同步獲取了空間分辨率為0.2 m的CCD影像共189幅，總覆蓋面積約為1.2 km×1.6 km。LIDAR數(shù)據(jù)保存成 LAS格式，投影方式為UTM，參考橢球為WGS84，每個激光點包含了激光點三維坐標(biāo)值、強度值、回波類型等信息。

2 數(shù)據(jù)處理

選擇了3種最為普遍有效的濾波方法對研究區(qū)進(jìn)行LIDAR數(shù)據(jù)濾波處理。

2.1 濾波原理

①不規(guī)則三角網(wǎng)算法。Axelsson提出了一種漸進(jìn)式不規(guī)則三角網(wǎng)加密算法。其基本原理是:首先選取區(qū)域內(nèi)的少量低點作為三角網(wǎng)的種子點，然后將點云中的點逐步加入到三角網(wǎng)中，如果要加入的點到三角形的距離(d)以及由于加入該點所形成的新三角形的內(nèi)角(α、β、γ)均小于所設(shè)定的閾值，則該點為地面點。在每輪加入點運算前，都須通過運算重新生成閾值。當(dāng)沒有點可再加入三角網(wǎng)時(即沒有點符合閾值標(biāo)準(zhǔn)時)，運算結(jié)束，形成最終的三角網(wǎng)［2］。這種方法通過反復(fù)迭代不斷向模型中添加新點，每個新加的點都使模型更加接近實際地表面。算法如圖1所示。

圖1 不規(guī)則三角網(wǎng)加密算法

圖2 回波垂直剖面示意圖

②多重回波算法。大部分激光雷達(dá)系統(tǒng)除了記錄高度信息外也記錄了回波及回波的強度信息，打在不同地物上的激光脈沖回波能夠反映不同的信息。對于植被，一個激光脈沖可以穿過樹葉、枝干到達(dá)地面，得到不同的回波信號，即返回不同的三維點信息［3］。激光脈沖經(jīng)發(fā)射后到達(dá)地面的過程中被不同高度的植被返回，首次回波一般是植被冠層表面返回的信息，而中間表面則包括第二次和中間次回波，地面表面則往往是由最后一次回波返回的點構(gòu)成的。多重回波過濾就是利用回波數(shù)目及回波次數(shù)的特性對點云進(jìn)行過濾，由最后一次回波可以得到DEM，而由首次和第二次回波能得到DSM。圖2是由不同高度的植被冠層反射回來的LIDAR激光脈沖信息［4］。

③迭代線性最小二乘內(nèi)插法。此法最初由奧地利維也納大學(xué)的Krasu和Pfeifer等人提出［5］。其核心思想是基于地物點的高程比對應(yīng)區(qū)域地形表面激光腳點的高程大，經(jīng)過線性最小二乘內(nèi)插后，激光腳點高程擬合殘差相對于擬合后的地形參考面不服從正態(tài)分布(見圖3)，高出地面的地物點高程擬合殘差都為正值且數(shù)值較大，而地面點的擬合殘差較小且可能為負(fù)值［6］。該方法首先用所有激光腳點的高程觀測值按等權(quán)計算出初步曲面模型(該曲面介于真實地面和植被覆蓋面之間，其結(jié)果使擬合后真實地面腳點的殘差為負(fù)值的概率變大)。然后依據(jù)每一個激光點到該新生成的表面的距離和方向利用穩(wěn)健估計的權(quán)函數(shù)關(guān)系式(Ⅰ)來計算其權(quán)值p:

式中:v為高程值與擬合面的殘差;p表示內(nèi)插時使用的權(quán)值;a和b決定內(nèi)插權(quán)函數(shù)的陡峭程度;g決定著哪種點的權(quán)值為1;地上偏移參數(shù)w指決定某一點對中間表面是否有影響的上限值。經(jīng)過幾次迭代后最終的裸露地表點從中間表面中抽取出來，所有高程值滿足權(quán)函數(shù)關(guān)系式(1)中第一二種情況的點都被作為裸露地面點。然后用這些點進(jìn)行內(nèi)插處理生成最終的DEM。

圖3 定權(quán)示意圖

圖4 去除明顯高點及孤立點

2.2 濾波試驗

為了說明這些方法的可行性，分別選擇了研究區(qū)內(nèi)不同區(qū)塊的點云數(shù)據(jù)進(jìn)行處理，這些區(qū)域既包含裸露地、水域，也有成片的樹林及房屋等地物。為了有效地過濾原始數(shù)據(jù)，使所獲取的DEM更加符合實際區(qū)域地形，首先對原始數(shù)據(jù)進(jìn)行粗差點剔除，包括剔除明顯的低點和孤立點以及空中懸浮點，再進(jìn)行濾波分類試驗。

圖5 消除噪聲點后按高程顯示的數(shù)據(jù)

①TIN三角網(wǎng)算法濾波效果。圖6是用于檢驗TIN三角網(wǎng)算法濾波效果的原始點云。該區(qū)域包含植被、河流湖泊、房屋及道路。首先根據(jù)試驗區(qū)地形選擇不同的最陡坡度、最大內(nèi)插角度和最大內(nèi)插距離。由于該區(qū)域比較平坦，上述參數(shù)分別設(shè)置為45°、55°、1 m時效果最好，設(shè)置最大建筑物為30 m，分離出地面點后便得到圖7的結(jié)果。對比原始點云及同步獲取的影像發(fā)現(xiàn)幾乎所有的地面點被有效地保留下來。然后利用Terrascan中的其他分類工具分出植被和建筑物點，這些過程都需要根據(jù)該區(qū)域的特點選擇合適的參數(shù)，才能得到比較理想的結(jié)果。

圖6 原始點云

圖7 濾波后獲得的地面點

②多重回波算法濾波效果。筆者選取試驗區(qū)中包含植被和人工建筑物的一片區(qū)域進(jìn)行濾波分類。圖8為根據(jù)不同回波類型記錄值分離出來的首次激光脈沖腳點，對比首尾回波可以發(fā)現(xiàn)直接或穿透植被打在地面上的點比較少，不利于構(gòu)建地面模型;圖9為尾次激光脈沖腳點;對比兩圖可以得出如下結(jié)論:對于森林地區(qū)，大部分打在植被冠層上的激光脈沖束不是直接返回而是向下繼續(xù)穿透直至抵達(dá)地面。但是能夠到地面的激光腳點數(shù)量還是不多。此外對于人工建筑物特別是屋頂，其上的激光腳點則往往只有一次回波，所以這部分腳點既屬于首次回波也屬于尾次回波。

圖8 首次回波點云

圖9 尾次回波點云

③迭代線性最小二乘內(nèi)插法效果。筆者選取了圖10左側(cè)區(qū)域(均為植被)進(jìn)行濾波實驗，設(shè)置參數(shù) a、b、g、w 分別為0.5、2、0、0.5，迭代 8 次，并設(shè)置濾波的格網(wǎng)大小為3 m，得到的結(jié)果如圖10右側(cè)所示，其中右上部分為植被部分的點云，而右下部分為過濾出的地面點云。

圖10 迭代線性最小二乘內(nèi)插法分類示意圖

圖10從激光點云角度反映了這一算法在提取大型建筑物及低矮灌木植被的缺陷。借助激光雷達(dá)同步獲取的影像發(fā)現(xiàn)圖11對應(yīng)的試驗區(qū)有大面積建筑和低矮灌叢，對于該試驗區(qū)，筆者試了不同的a、b、g和w參數(shù)，但是效果一直不理想，總是有殘存的植被點和建筑物。

3 模型生成

采用上述不同算法濾波后筆者將各自的濾波結(jié)果保存成ASCII文件，再生成相應(yīng)的數(shù)字高程模型與數(shù)字表面模型。

圖11 迭代線性最小二乘內(nèi)插算法錯分的地面點

3.1 由TIN濾波結(jié)果構(gòu)建

圖12即為用TIN濾波所得地面點構(gòu)建的DEM，與傳統(tǒng)的DEM相比，由該方法構(gòu)建的地面模型更能反映地表的實際情況，微小的地形起伏很好地被展現(xiàn)出來。而圖13為DEM渲染及等高線疊加顯示圖，充分反映了TIN三角網(wǎng)濾波算法的適用范圍，此圖右半部分地形比較平坦，而左半部分是起伏坡度較大的山，可以看出平坦地區(qū)的等高線特別細(xì)碎，構(gòu)建的DEM也很破碎;而用左邊陡峭地區(qū)濾出的地面點生成的等高線平緩又光滑，證明該部分DEM質(zhì)量很好。

圖12 TIN三角網(wǎng)構(gòu)建的DEM

圖13 平坦與陡峭地區(qū)DEM對比

3.2 由多重回波濾波結(jié)果構(gòu)建

圖14和圖15分別是首回波和尾回波生成的部分DSM，對比可看出首回波蘊含的植被信息。圖16為首回波脈沖腳點的高程同尾次回波脈沖腳點的高程作差后形成的正規(guī)化數(shù)字表面模型(nDSM);在森林地區(qū)nDSM被稱作冠層高度模型(CHM)，它反映的是冠層的高度信息。可看出這個模型幾乎呈平面，沒有山體那種大的起伏或坡度，因此正規(guī)化的表面模型，消除了傳統(tǒng)表面模型中地形起伏變化對地物高程及形狀的干擾，能獲取更準(zhǔn)確的地物形態(tài)和高度信息。

圖14 首回波生成的DSM

圖15 尾回波生成的DSM

圖16 首尾兩次回波差生成的nDSM

3.3 由迭代線性最小二乘內(nèi)插結(jié)果構(gòu)建

圖17和圖18分別為利用迭代線性最小二乘法過濾出的地面點與非地面點生成的DEM及DSM。可知，對于植被生長狀況良好的林區(qū)，迭代線性預(yù)測算法的效果很好，幾乎所有地面點都能被分離出來，故生成的DEM能很好地反應(yīng)微小的地形起伏。圖19反映了該算法在提取大型建筑物及低矮灌木植被的不足。

圖17 迭代線性法分離出的地面點生成的DEM

圖18 迭代線性法分離出的植被點生成的DSM

圖19 迭代線性法失效生成的DEM

4 效果評價

激光雷達(dá)點云濾波分類之后還需要進(jìn)行濾波效果評價，評價分為定性評價與定量評價兩種。定量分析是統(tǒng)計濾波結(jié)果中有多少非地面點被誤分為地面點，有多少地面點被誤分為非地面點，從而得出TypeⅠ型錯誤TypeⅡ型錯誤的百分比［7］，然而對不同地區(qū)的海量數(shù)據(jù)進(jìn)行統(tǒng)計往往是很困難的。因此本研究僅從定性的角度評價3種算法的濾波質(zhì)量，具體為:可視化分析(利用分出的地面點在ArcGIS中生成渲染圖，通過縮放觀測表面模型的構(gòu)建狀況)、等高線分析(根據(jù)等高線的疏密、平緩、光滑程度等特性來判斷濾波結(jié)果好壞)、剖面圖分析(直接針對濾波后的點云進(jìn)行分析，從立面評價濾波效果，判斷點是否錯分)。

5 結(jié)論

采取這些手段評價分析濾波分類效果后，得出或驗證了各方法的適用范圍及優(yōu)缺點，結(jié)論如下:①對高程突變地物，TIN不規(guī)則三角網(wǎng)迭代算法的過濾效果較好。這是由于高大建筑物和植被與其鄰近地面點之間形成了明顯的高程突變，但在過濾灌叢或低矮的地面物體時，產(chǎn)生過大誤差。②利用激光雷達(dá)的不同地物對應(yīng)回波的高程差不同進(jìn)行濾波的多重回波算法簡單有效，根據(jù)植被點和與其他點的回波類型情況不同能夠很快地把植被點與其他激光腳點分開，特別適用于分類植被激光腳點;但由于穿透率的問題，末回波點也即地面點稀少，因此利用末回波構(gòu)建的DEM不能很好地反映出精細(xì)地形。③迭代線性最小二乘內(nèi)插算法借助鄰近激光腳點間的高程突變來區(qū)分地面點與非地面點，能很好地獲得地形趨勢面。但是該方法在有大型建筑物存在的區(qū)域往往不適用，它僅能削去房屋的棱角，而無法完全濾掉建筑物。對地形起伏變化不大的森林地區(qū)，其效果則是其他方法無法比擬的。④筆者根據(jù)涼水自然保護(hù)區(qū)特有的地形、植被等因素進(jìn)行分塊，然后對每一分塊采用適合其地形特點的濾波算法進(jìn)行處理，所得效果較佳。對平坦地區(qū)設(shè)置不同的參數(shù)采用迭代線性最小二乘內(nèi)插法進(jìn)行處理，對有許多高大植被和建筑物的地區(qū)采用不規(guī)則三角網(wǎng)算法濾波，對森林地區(qū)采用多重回波算法分離出地面點集，再利用各小分塊間重疊區(qū)域的特征進(jìn)行拼接得到?jīng)鏊貐^(qū)完整的地面點集并生成DEM。

［1］周淑芳，李增元，等.基于機載激光雷達(dá)數(shù)據(jù)的DEM獲取及應(yīng)用［J］.遙感技術(shù)與應(yīng)用，2007，22(3):356－360.

［2］曾齊紅.機載激光雷達(dá)點云數(shù)據(jù)處理與建筑物三維重建［D］.上海:上海大學(xué)，2009.

［3］張小紅.利用機載LIDAR雙次回波高程之差分類激光腳點［J］.測繪科學(xué)，2006，31(4):48－50.

［4］Akay A E，Oˇguz H，Karas I R，et al.Using LiDAR Technology in Forestry Activities［J］.Environ Monit Assess，2009，151:117－125.

［5］Krasu K，Pfeifer N.Determination of Terrain Models in Wooded Areas with Airborne Laser Scanner Data［J］.ISPRS Journal of Photogrammertry and Remote Sensing，1998，53(4):193－203.

［6］Pfeifer N，Reiter T，Briese C，et al.Interpolation of High Quality Ground Models from Laser Scanner Data in Forested Areas［J］.International Archives of Photogrammertry and Remote Sensing，1999，32(3/W14):31－36.

［7］Sithole G，Vosselman G.Experimental Comparison of Filter Algorithms for Bare－Earth Extraction from Aiborne Laser Scanning Point Clouds［J］.ISPRS Journal of Photogrammertry ＆ Remote Sensing，2004，59:85－101.