基于不同濾波方法的機載激光雷達數據處理及表面模型構建1)

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

(東北林業大學，哈爾濱，150040)

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

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

1 研究區概況及數據獲取

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

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

2 數據處理

選擇了3種最為普遍有效的濾波方法對研究區進行LIDAR數據濾波處理。

2.1 濾波原理

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

圖1 不規則三角網加密算法

圖2 回波垂直剖面示意圖

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

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

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

圖3 定權示意圖

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

2.2 濾波試驗

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

圖5 消除噪聲點后按高程顯示的數據

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

圖6 原始點云

圖7 濾波后獲得的地面點

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

圖8 首次回波點云

圖9 尾次回波點云

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

圖10 迭代線性最小二乘內插法分類示意圖

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

3 模型生成

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

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

3.1 由TIN濾波結果構建

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

圖12 TIN三角網構建的DEM

圖13 平坦與陡峭地區DEM對比

3.2 由多重回波濾波結果構建

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

圖14 首回波生成的DSM

圖15 尾回波生成的DSM

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

3.3 由迭代線性最小二乘內插結果構建

圖17和圖18分別為利用迭代線性最小二乘法過濾出的地面點與非地面點生成的DEM及DSM?？芍瑢τ谥脖簧L狀況良好的林區，迭代線性預測算法的效果很好，幾乎所有地面點都能被分離出來，故生成的DEM能很好地反應微小的地形起伏。圖19反映了該算法在提取大型建筑物及低矮灌木植被的不足。

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

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

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

4 效果評價

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

5 結論

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

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

［2］曾齊紅.機載激光雷達點云數據處理與建筑物三維重建［D］.上海:上海大學，2009.

［3］張小紅.利用機載LIDAR雙次回波高程之差分類激光腳點［J］.測繪科學，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.