基于LiDAR點云濾波技術的風電場地形圖處理方法

張 兵，武殿文

（1.內蒙古電力勘測設計院有限責任公司，內蒙古 呼和浩特 010020；2.內蒙古能源規劃設計研究院，內蒙古 呼和浩特 012100）

0 引言

風電場地形圖測繪如今普遍采用無人機低空航測技術，以便獲取測區高精度數字正射影像圖(digital orthophoto map，DOM)和數字地表模型(digital surface model，DSM)。通過功能強大的無人機數據后處理軟件一般能夠得到較為理想的測區DOM，即使在非地面點復雜地區也不會有太多影響其精度的地方。但DSM對于非地面點的反應卻十分敏感。對于植被十分茂密的風電場測區，往往由于樹木高度的不確定性很難得到精確地形數據，進而影響風電場地形圖測繪的準確度。所以必須剔除非地面點對數字高程模型(digital elevation model，DEM)的影響。剔除非地面點的解決方法，傳統主要采用手工采集，然后減去樹木、房屋等地物點的高度，得到地面點的高程點繪制等高線，或在每個區域內部采取野外實測樹高，或者在航測內業采用立體測圖的辦法獲取樹高，對每個植被覆蓋區域采取減去固定高程差值的辦法來得到整個風電場區域的DEM。上述傳統方法采取人工干預的方式來剔除測區內植被或者房屋等非地面點，人工工作量很大、效率很低。另一種途徑是借鑒激光雷達(light detection and ranging，LiDAR)點云的處理方法對多視影像匹配點云處理，將已獲取的點云數據中地面點與非地面點分離，獲取剔除非地面點后的DEM，用于風電場矢量地形圖繪制。這種辦法人工參與少、作業效率高，在實際應用中反復驗證點云濾波后DEM精度情況顯得尤為重要。

1 LiDAR點云濾波技術的發展

點云數據濾波方法最早來自激光點云的數據處理，一直是激光點云應用最前沿和最難解決的棘手問題，機載LiDAR點云數據的濾波算法仍相對落后，許多學者對點云濾波算法進行了研究，并提出了多種濾波算法，可概括為2種：①利用點云數據的幾何特性(三維坐標)進行的濾波； ②利用點云數據的回波特性(回波強度、多回波信息) 進行的濾波[1-2]。LiDAR點云濾波技術通過近些年的發展，出現很多基于點云數據的幾何特性(數據本身)的濾波方法，有基于坡度、基于曲面擬合、基于數學形態學和基于布料模擬等。基于數據本身的濾波方法，這一過程是利用從數據內部提取的有用信息建立判別規則，實現過濾非地面點，獲取只包含地面點的點云數據的目標，這種濾波方法的對象不只是濾除植被，還包括了人工地物、移動物體等表面噪聲。基于數據本身的濾波目的是濾除點云數據或DSM的表面噪聲，在濾除高低不一的建筑物或孤立植被方面具有較好效果，但在大范圍植被覆蓋地區或成片森林，立體影像匹配生成的點云數據通常不包括地面點，因此不能滿足這類算法所要求的局部鄰域內最低點就是地面點的初始條件，無法有效濾除植被。但是針對多視影像匹配點云的濾波只能采用基于數據本身的濾波方法。目前基于數據本身的LiDAR點云濾波技術主要使用的方法包括基于坡度的濾波算法、基于不規則三角網濾波方法、基于移動曲面擬合的濾波算法、基于數學形態學濾波算法以及基于布料模擬的濾波算法[3-5]。本文將分別論述以上幾種濾波方法，并重點研究基于布料模擬的濾波算法。將基于布料模擬濾波算法研究成果應用于風電場地形圖等高線繪制。

2 多種基于數據本身LiDAR濾波技術的算法原理與分析

基于坡度的濾波算法由VOSSELMAN首先提出，其思想是依據坡度的差異區別出地面點和非地面點，但對整個點云區域采用固定坡度閾值的方式不夠合理[6]。基于坡度的濾波算法原理認定地形表面為平緩光滑曲面，局部區域地形急劇變化很大可能由非地面點引起，相鄰點間高差值和距離值的濾波函數被定義為地面點歸類函數，如式(1)所示。

式(1)中點pi與pj的差值不大于Δhmax(d)即可將pi判定為地面點，否則為非地面點。

圖1中基于不規則三角網濾波算法是通過初始地面點p1、p2、p3與待判定點p構建不規則三角網，分別計算d、a1、a2、a3與設定的閾值比較，如果小于閾值將p歸類為地面點并加入初始三角網，通過反復迭代計算直到沒有新的地面點結束。

圖1 不規則三角網略圖

基于移動曲面擬合算法強調地形的平滑性，將點云數據二位排序獲取地面點種子區域，計算種子點曲面方程并將臨近點代入計算擬合高程值，若不低于設定閾值，判定為地面點并吸收到二次曲面方程中，通過反復迭代計算直到所有點云判定結束，在改變網格尺寸和設置不同的窗口領域大小時，為了有效地將地面點和非地面點分開，需要設置不同的高差閾值[7]，如式(2)所示。

基于數學形態學的濾波算法對點云數據先進性腐蝕運算再進行膨脹運算，將區域內最小和最大高程值填充到結構元素內部，然后運用開閉運算反復迭代出點云數據內地面點和非地面點。形態學濾波主要依靠開運算實現點云類屬判斷[8-9]。主要運算式如式(3)～式(6)所示。

基于布料模擬濾波算法是根據牛頓第二定律原理在模擬過程中內部力和外部力同時作用且隨時間變化，假設布料足夠柔軟黏住地面則布料最終形狀會是DSM，考慮將點云倒置給布料賦予硬度因素則布料形狀貼近數字地面模型，物理表達式如式(7)所示。

在實際應用中，統計很多實際案例發現基于坡度的濾波算法對多數地形都具有適用性，需要根據不同地形設定坡度閾值，不適合復雜地形，坡度閾值的選定也較為困難。基于不規則三角網的濾波算法適用于地形起伏不大的地區，但是構建三角網運算量很大，距離和角度閾值設定很復雜且不易剔除植被。基于移動曲面擬合的濾波算法原理簡單算法速度快，但函數高程閾值很難判斷，窗口尺寸設置不好把握。基于數學形態學的濾波算法簡單易于實現，但在處理原始點云數據中容易丟失數據，固定的濾波窗口和閾值使濾波精度不高。而基于布料模擬濾波算法在實際應用中適用于大多數地貌類型，對于各種地物類型的濾波效果也更為優異，已被更多應用于DSM的點云濾波。在吸收很多人關于點云相關研究的基礎上人們建立起來一套大型跨平臺開源編程點云庫(point cloud library，PCL)，PCL實現了大量點云相關的通用算法，其中就包括各類點云濾波算法。利用PCL模塊化的現代編程模板庫實現點云濾波算法實踐顯得十分高效。

PCL可結合多種程序來編譯濾波算法，實現點云數據的自動濾波，提供了大量數據類型和數據結構，以及許多最先進的算法對3D數據進行處理，包括濾波、模型重估、表面重建等。PCL中的所有模塊和算法都是基于標準的第三方庫Boost、Eigen、FLANN、VTK、CUDA、OpenNI、Qhull來傳送處理數據的，避免多次復制系統中已存在的數據的需要。PCL利用OpenMP、GPU、CUDA等先進高性能計算技術，通過并行化提高程序實時性。從算法的角度，PCL是指納入了多種操作點云數據的三維處理算法。每一套算法都是通過基類進行劃分，把貫穿整個流水線處理技術的所有常見功能整合在一起，保持整個算法實現過程中的緊湊和結構清晰，提高代碼的重用性、簡潔性。

3 CSF算法原理與程序實現

基于布料模擬濾波(cloth simulation filter，CSF)算法相較于其他濾波算法具有算法精度高、適用面廣等優點，本文把基于布料模擬的濾波算法作為無人機多源數據匹配生成DSM的主要濾波方法。布料模擬是物理過程的模擬，符合牛頓第二定律原理。模擬布料落到倒置的點云表面上，分析布料的節點與相應的點之間的相互作用，以約定的時間內模擬布料的形狀，而布料節點在垂直方向上的位移依照式(8)所示。

算法的核心是根據PCL使用編程語言完成CSF算法編譯。其大致的實施過程是將原始點云數據倒置，設定布料模擬參數和初始位置。然后將點云數據與柵格質點投影到同一平面，結合式(6)和式(7)計算模擬布料所有質點的垂直位移，使用位移后的布料質點與初始點云之間進行距離計算，并與設定閾值比較，距離小于分類閾值則將此點歸為地面點，距離大于分類閾值則將此點歸為非地面點。如圖2所示是基于CSF算法簡易流程圖。

圖2 CSF簡易流程圖

具體在PCL中，根據提供的多種開發語言的點云處理技術算法源碼。首先在PCL中定義函數庫以及需要使用的規范文件，讀取和保存LAS格式點云數據。在進行CSF過程中關鍵在于濾波參數的設置，需要設置質點為不可移動，設定布料網格分辨率和布料硬度，設定布料模擬質點與原始點云的距離閾值以及重復計算最大迭代次數。以下是CSF基于LAS格式點云python簡易源碼。

# coding: utf-8

import laspy

import CSF

import numpy as np

import open3d as o3d

inFile = laspy.file.File(r” in.las”,mode=’r’)# read a las file

points = inFile.points

xyz=np.vstack((inFile.x,inFile.y,inFile.z)).transpose()# extract x,y,z and put into a list

csf = CSF.CSF()

# prameter settings

csf.params.bSloopSmooth =False

csf.params.cloth_resolution =0.2

csf.params.rigidness =5

csf.params.time_step =0.65

csf.params.class_threshold =0.5

csf.params.interations =200

#moredetailsaboutparameter:

csf.setPointCloud(xyz)

ground = CSF.VecInt()#

non_ground = CSF.VecInt()#

csf.do_filtering(ground,non_ground)#

outFile = laspy.file.File(r”ground.las”,

mode=’w’,header=inFile.header)

outFile.points = points[ground]# extract ground points,and save it to a las file.

outFile.close()#

4 基于DSM點云濾波方法的風電場地形圖應用

無人機低空攝影制圖技術用于提供風電場地形圖需要的DOM和DSM。多視影像匹配技術的發展帶來更為精細化的2D數據，針對無人機2D數據的風電場矢量地形圖繪制已經十分成熟。但是在內蒙古區域內依然存在不少地物植被眾多區域，因此在處理風電場矢量地形圖時不可避免對非地面點進行過濾和DEM重建。將LiDAR點云濾波技術引入到多視影像匹配DSM濾波，能夠提高后續制圖成果的效率和精度。如圖3所示是加入點云濾波步驟后的風電場地形圖制作流程圖。

圖3 風電場地形圖制作流程

數據濾波前后針對濾波區域進行了精度檢校。圖4為無人機低空攝影測量成果制作的DSM，圖5為采用布料濾波算法進行點云濾波后重新構建生成的DEM。圖6為采用原始DSM制作的等高線圖形，圖7為采用DEM制作的風電場等高線圖形。通過加入DSM濾波技術，能夠顯著提高風電場地形圖等高線制作的效率和精度。通過以上案例發現，在對原始DSM進行自動濾波后，大部分非地面點都被有效剔除，濾波后重構DEM能夠貼合真實地形(圖5)，但是在實際檢測中發現部分正常的地形點數據經過濾波后被帶偏移，部分地區高程點數值改變超過限差以上。針對以上問題在對現有數據進行濾波后再與原有DEM數據進行一次融合算法。該算法原理是首先對2次數據進行匹配，計算機自動鎖定區域內所有濾波區域然后計算出每個濾波區域的閉合界線。設定在閉合線里邊的數據采用濾波后數據，而在閉合線外部的數據替換成原有數據。在獲得更新數據之后，又進行了DEM精度檢校。采取以上辦法可以很好地解決部分正常數據被濾波數據帶偏的問題。最終再次重構DEM并生成等高線數據(圖7)，該等高線數據能夠真實體現風電場地形。

圖4 數字地表模型DSM

圖6 原始DSM生產等高線

圖7 重構高程模型生產等高線

針對點云濾波精度評定方法有目視定性評估法和統計學定量評價法：①目視定性評估法：是觀察濾波后地面點數據是否符合風電場地形圖制圖要求；②統計學定量評價法：點云濾波結果時定量分析第I類誤差、第II類誤差和總誤差。設定a為地面點誤分為非地面點的數量，b為非地面點誤分為地面點的數量，c為正確劃分地面點的數量，d為正確劃分非地面點的數量。誤差評價公式如式(9)～式(11)所示。

點云濾波算法一般要求在控制Ⅱ類誤差的基礎上盡量降低Ⅰ類誤差，即在保證分離出的地面點最大程度上是地面點的同時，盡可能分離出較多的地面點[10]。以上3類誤差比率系數值越小表明點云濾波效果越好。在多個工程實際應用中統計表明，能夠滿足風電場測圖需求的點云濾波3類誤差統計結果應滿足Ⅰ類誤差≤0.1，Ⅱ類誤差≤0.05，總誤差≤0.05。

5 結論

LiDAR點云濾波技術發展到今天，很多前人的經驗在經過發展之后能更好地對LiDAR點云數據進行濾波。將LiDAR點云濾波技術用于多視影像匹配獲取的DSM點云濾波，由于無人機航測數據的局限性只能在一定程度上解決植被等非地面點的干擾來獲取DEM，并且通過實際檢校發現很難達到更高精度。僅在風電場地形圖無人機航測實際應用中通過多次的濾波成果與工測數據對比，發現大多數能夠滿足風電場地形圖測繪精度要求。因此對于植被豐富地區的地形圖測繪，為了達到厘米級精度要求，想通過DSM濾波技術方法得到高精度DEM不現實：一方面需要廣泛運用無人機機載LiDAR技術采集測區激光點云，通過機載LiDAR多次回波技術直接穿透植被獲取真地表地形數據，才是一次性解決高精度真地形模型的辦法；另一方面由于單一一種LiDAR點云濾波技術都存在自身算法和閾值設定的局限性，因此很難確定到底哪種方法能夠完美適用于各種地形風電場地形圖的DSM濾波。對于植被豐富地區的地形圖測繪，無法滿足更高精度要求，需要更加深入地將多種LiDAR點云濾波技術合并使用，同時積極采納圖像學、人工智能學最新研究成果，提供適用于所有地形場景風電場DSM數據的點云濾波，以期達到更高精度和濾波閾值設定更智能化。