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

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

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

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

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

研究區(qū)域?yàn)楹邶埥鏊畤壹壸匀槐Ｗo(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。采樣點(diǎn)云密度大于2點(diǎn)/m2，飛行高度達(dá)800 m，平均對地飛行速度是180 km/h。此外同步獲取了空間分辨率為0.2 m的CCD影像共189幅，總覆蓋面積約為1.2 km×1.6 km。LIDAR數(shù)據(jù)保存成 LAS格式，投影方式為UTM，參考橢球?yàn)閃GS84，每個激光點(diǎn)包含了激光點(diǎn)三維坐標(biāo)值、強(qiáng)度值、回波類型等信息。

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

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

2.1 濾波原理

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

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

圖2 回波垂直剖面示意圖

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

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

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

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

圖4 去除明顯高點(diǎn)及孤立點(diǎn)

2.2 濾波試驗(yàn)

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

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

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

圖6 原始點(diǎn)云

圖7 濾波后獲得的地面點(diǎn)

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

圖8 首次回波點(diǎn)云

圖9 尾次回波點(diǎn)云

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

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

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

3 模型生成

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

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

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

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

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

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

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

圖14和圖15分別是首回波和尾回波生成的部分DSM，對比可看出首回波蘊(yùn)含的植被信息。圖16為首回波脈沖腳點(diǎn)的高程同尾次回波脈沖腳點(diǎn)的高程作差后形成的正規(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分別為利用迭代線性最小二乘法過濾出的地面點(diǎn)與非地面點(diǎn)生成的DEM及DSM?？芍瑢τ谥脖簧L狀況良好的林區(qū)，迭代線性預(yù)測算法的效果很好，幾乎所有地面點(diǎn)都能被分離出來，故生成的DEM能很好地反應(yīng)微小的地形起伏。圖19反映了該算法在提取大型建筑物及低矮灌木植被的不足。

圖17 迭代線性法分離出的地面點(diǎn)生成的DEM

圖18 迭代線性法分離出的植被點(diǎn)生成的DSM

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

4 效果評價

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

5 結(jié)論

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

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

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

［3］張小紅.利用機(jī)載LIDAR雙次回波高程之差分類激光腳點(diǎn)［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.