機載LiDAR測量復雜地形中濾波算法的比較研究

李 豆, 李朋飛, 穆興民, 姚頑強, 湯伏全, 李 婷

(1.西安科技大學 測繪科學與技術學院, 西安 710054;2.西北農林科技大學 黃土高原土壤侵蝕與旱地農業國家重點實驗室,陜西 楊凌 712100; 3.中國科學院 水利部 水土保持研究所, 陜西 楊凌 712100)

黃土高原地形破碎復雜，高精度地形信息的高效獲取是區域微地形變化監測的關鍵，可為土壤侵蝕、地表沉降等過程的研究奠定基礎。常用的地形獲取方法有實時動態載波相位差分技術、攝影測量、干涉雷達，但消除植被影響的濾波算法尚不成熟，影響其地形獲取的精度[1-3]。激光雷達(Light Detection And Ranging，LiDAR)是新興的主動遙感技術，精度可達mm級，且穿透植被能力強[4-7]，其通過發射激光束照射物體表面并接收返回信號，測得激光束發射至返回的時間，推算物體的距離，獲取每一個光斑的三維坐標，形成高密度三維點云[8-10]。機載LiDAR是LiDAR發展的重要方向，由LiDAR、全球定位導航系統、慣性導航系統等組成[11]，可快速、高效地獲取大范圍區域的高精度地形信息。在國外，機載LiDAR獲取地形信息的研究已有報道[12]，然而其在國內尤其黃土高原的應用鮮見報道。

機載LiDAR獲取的原始點云包含植被、地形、建筑等信息，需濾波以準確提取地面點云，方可實現地形信息的高精度表達[13]。國內外學者已開發眾多點云濾波算法[14]，可分為基于坡度濾波算法[15]、基于形態學濾波算法[16]、基于曲面擬合濾波算法[17]、基于不規則三角網的濾波算法等[18]類型。不同濾波算法已應用于世界不同地區，但在黃土高原地區的應用鮮有報道。地形、植被和點云密度是濾波精度的主要影響因素[19]。黃土高原地形復雜，陡峭區域的點云獲取較為困難，加之近年來植被恢復[20]，部分區域植被覆蓋度較高，為已有濾波算法的應用帶來挑戰。因此，需明確各影響因素對不同濾波算法的影響，為算法的選取和改進提供依據。

本文以黃土高原溝壑區董莊溝為研究區，選用多尺度曲率濾波算法(Multi-scale curvature classification,MCC)、擴展窗口高程閾值算法(Elevation threshold with expand window,ETEW)、適應三角網濾波算法(Adaptive Triangle Irregular Network,ATIN)、漸進形態學濾波算法(Progressive Morphological Filtering,PM)、基于坡度的濾波算法(Slope Based Filtering,SBF)5種算法濾波獲取的無人機LiDAR點云。分析不同算法的濾波精度及其與影響因素(坡度、植被覆蓋、點云密度)的關系，為黃土高原高精度地形信息的高效獲取提供參考。

1 研究區概況

董莊溝位于甘肅省慶陽市西峰區境內(107°30′—107°37′E，35°41′—35°44N)，是黃土高原溝壑區典型小流域之一(圖1)。董莊溝溝長1 600 m，溝道比降8.93%，流域面積1.15 km2，海拔1 135～1 350 m。年均氣溫9.3℃，年均降水量546.9 mm，其中6—9月降水量占年降水量的69.9%。流域為第四紀黃土所覆蓋，總厚度約250 m。地貌類型主要包括塬面、梁峁坡、溝谷。塬面地形平坦，坡度在5°以下，多為農耕地；梁峁坡為連接塬面的緩坡帶，一般為10°～20°；梁峁坡以下為溝谷，呈V形，溝谷坡一般為40°～60°的陡坡和大于60°的立壁。流域內以荒草地為主，除村莊、道路旁和部分溝頭有小型林帶外，無整塊大片林帶。

2 數據源與研究方法

2.1 數據獲取

利用北科天繪機載LiDAR系統Sky-Lark于2018年7月采集獲取點云數據。數據采集中，無人機LiDAR采用矩形航線沿董莊溝流域邊界飛行，掃描視場角70°×40°，回波次數4次，激光采樣頻率300 kHz，掃描轉速3 000 r/min，航高190 m，飛行速度17 km/h。數據采集完成后進行飛行軌跡和姿態解算、點云配準和精化、坐標轉換等預處理，生成.las格式點云數據。依據所獲取點云，在董莊溝流域不同位置選取12個樣本區域(圖1，表1)，用于本次研究。樣本選取確保涵蓋流域內的典型地貌并保證地貌的完整性。

圖1 董莊溝及12個樣本區域的位置

表1 各樣本區域信息

2.2 數據處理與分析

2.2.1 濾波算法 (1) 多尺度曲率濾波算法(MCC)。MCC定義一個包含所有點x，y，z坐標的向量Z(s)，x和y為平面坐標，z為點云高程。利用薄板樣條插值法Z(s)得到連續柵格表面。在柵格表面傳遞一個3×3的格網定義一個新的向量X(s)，其z為3×3柵格的平均高程；比較Z(s)中各點的高程與X(s)高程和曲率容差t之和的大小，當Z(s)

(2) 擴展窗口高程閾值算法(ETEW)。ETEW根據相鄰地面點的高程差和地面點與非地面點高程差的差異進行點云分類，通過對比預設半徑內點的高程差與預先設置的閾值，剔除非地面點[22]。首先將點云數據細分成方形單元，保留單元內最低高程的點；其次計算單元中高程最低點和其他點的高程差，高程差小于閾值時視為地面點，大于閾值則為非地面點。隨著單元格和閾值的增大，重復此過程，直到不再剔除非地面點。ETEW分類的效果和單元格大小、坡度因子、迭代次數有密切關系[22]。

(3) 適應三角網濾波算法(ATIN)。ATIN以點到三角形不規則網絡(TIN)表面的距離為主要約束條件，進行濾波計算。首先將整體區域劃分成小的單元格網，每個單元格只保留高程最低的點(種子點)構建初始稀疏的TIN；然后計算每個點到最近三角形表面的距離以及點到三角形頂點的直線與三角形表面的夾角。當距離和角度都小于預先設定的閾值時，則將該點歸類為地面點；后續迭代中使用這些基本點來構造新的TIN，計算其余點到表面的距離和角度，并與閾值進行比較，直到未分類的點都不可添加到地面點集時，迭代終止[18]。該算法分類結果和單元格的選取、初始三角網的大小關系最密切。

(4) 漸進形態學濾波算法(PM)。數學形態學濾波算法基于集合理論從圖像中提取特征，采用腐蝕和膨脹兩個基礎算子組成“開”運算對LiDAR點云數據分類。Zhang等[16]在此基礎上開發了PM，其核心原理是逐漸增大窗口尺寸及對應的高差閾值，在保留地面數據的同時，去除不同尺寸的非地面目標的測量值。首先將矩形格網覆蓋在點云數據集上，每個單元格網包含一個測量點，其為單元格內的最低高程點。單元格中測量點的高程構成一個初始近似曲面；然后進行“開”運算，獲得次生曲面，將第i次的曲面與第(i-1)次的曲面之間高程差值和閾值進行比較，判斷單元格網中的點是否為非地面點[23]。單元格和起始窗口大小對PM結果影響較大。

(5) 基于坡度的濾波算法(SBF)。SBF通過坡度差來分離地面點和非地面點。算法原理是將一個矩形網格覆蓋在點云數據集上，每個激光點分配到一個單元網格，如果同一單元格中有多個點下落，則選擇高程最低的點作為數組元素；LiDAR點在給定半徑范圍內，該點與任意其他點之間斜率的最大值小于預先設定的閾值時，被劃分為地面點[15]。單元格大小、搜索半徑、最大坡度對該算法影響較大。

2.2.2 精度評估 (1) 參考數據獲取。通過對比參考點云數據與不同算法的分類結果，確定濾波精度。參考數據可通過純人工濾波、GNSS/全站儀測量、自動+手動濾波等方法獲取[22]。純人工濾波是作業人員根據對研究區的先驗知識，通過手動分類實現地面點與非地面點的區分，分類精度受作業人員知識、經驗等影響較大。GNSS/全站儀測量所得地面點數據密度較低，且野外工作量大，不易實現。“自動+手動”方法可兼顧點云精度和工作量，已被廣泛用于獲取參考點云數據[18]。因此，本研究選用“自動+手動”方式，在Terrasolid 2016軟件中處理得到參考點云。首先，通過Terrasolid 2016內置濾波算法剔除雜點和噪聲點，初步分離地面點和非地面點。然后，結合無人機LiDAR獲取的董莊溝高分辨率照片，開展目視判讀，修正植被覆蓋度、地形陡峭等敏感區域的分類結果，獲取高精度的點云分類數據。

(2) 精度評價。點云分類精度評價采用國際攝影測量與遙感學會推薦的方法。該方法計算3類誤差，分別為Ⅰ類誤差、Ⅱ類誤差和總誤差[17]。其中Ⅰ類誤差表示地面點云錯分為非地面點云的比例，Ⅱ類誤差表示非地面點云錯分為地面點的比例，總誤差代表總體點云錯分比例。

(1)

(2)

(3)

式中:a為正確分類的地面點數量;b為地面點誤分為非地面點的數量；c為非地面點誤分為地面點的數量;d為正確分類非地面點的數量。

2.2.3 濾波算法參數確定 高精度的地形獲取需要完整的地面點云，以便有效監測微地形變化。Ⅰ類誤差導致的地面點云缺失難以在后續數據處理中消除，是地形信息準確表達的主要難點，而Ⅱ類誤差(非地面點云誤分為地面點云)則相對容易去除[17]。因此，本文濾波中，在盡可能保證地形信息完整的同時，降低各類誤差，各參數初始取值為允許范圍的下限，以固定步長逐步增加參數取值，直至參數取值上限，對比不同參數取值下的濾波誤差，獲取參數的最優取值(表2)。

表2 各濾波算法主要參數的最優取值

3 結果與分析

3.1 不同濾波算法的精度對比

各濾波算法關鍵參數的取值見表2，在最優參數下，各樣本區域Ⅰ類誤差中，MCC和PM結果通常最低，ETEW和SBF結果最高，ATIN結果處于中間，但偶爾低于MCC結果。Ⅱ類誤差中，SBF和ETEW結果最低，MCC和PM結果最高，ATIN介于中間。12個樣本的不同濾波方法誤差見圖2。分析表明：樣本5，6及8—12中，Ⅰ類誤差大小順序為MCC

圖2 5種濾波算法的Ⅰ類誤差、Ⅱ類誤差和總誤差

利用ArcGIS中的地形轉柵格法(Topo to Raster)插值點云生成連續表面，獲取地形陰影圖(圖3)，對比分析不同算法濾波效果的空間分布。根據地形陰影圖，PM結果能夠完整保留地形陡峭區域(圖3中的框線區域)，ATIN能夠部分保留地形陡峭區域，而MCC，ETEW，SBF剔除了陡峭區域，導致陡峭區域地形模型細節缺失。這說明PM能較好適用于陡峭區域的地面點云濾波，ATIN部分適用于陡峭區域濾波，而MCC，ETEW和SBF難以用于陡峭區域。綜合考慮，PM能夠較好地適應黃土高原的復雜地形地貌特征。

圖3 5種濾波算法在樣本1區域濾波結果

3.2 影響濾波精度的因素

3.2.1 坡 度 基于各樣本數據的地面點云參考數據，生成0.5 m分辨率的DEM，計算不同樣本區域的平均坡度(表1)，分析坡度對不同濾波算法精度的影響。結果表明，各算法的Ⅰ類誤差與平均坡度呈顯著正相關(圖4)(p<0.05)，而Ⅱ類誤差和總誤差與坡度的相關性不顯著(p>0.05)。MCC和PM的Ⅰ類誤差隨坡度上升速率最小，ETEW的Ⅰ類誤差隨坡度上升速率最大，SBF和ATIN的Ⅰ類誤差隨坡度上升速率居中。坡度小于45°的9個樣本區域中，6個樣本區域的MCC的Ⅰ類誤差最小，坡度最高的3個樣本區域(48.54°，49.72°，54.95°)PM的Ⅰ類誤差最低。

圖4 坡度與濾波算法誤差的關系

Ⅰ類誤差較Ⅱ類誤差和總誤差對坡度變化更敏感。原因在于各算法傾向于將陡峭區域整體誤分為非地面點，造成Ⅰ類誤差的迅速上升，而陡峭區域點云中包含的部分真實非地面點被正確分為非地面點云，緩慢降低了Ⅱ類誤差。各算法Ⅰ類誤差對坡度變化敏感性差別的原因可能在于，ETEW僅依據高差閾值判斷歸類點云，不能有效考慮陡峭區域的影響，SBF和ATIN通過選擇合適的坡度(或夾角)閾值一定程度上降低了坡度對濾波結果的影響，而MCC利用曲率容差有效考慮了坡度的影響，因此Ⅰ類誤差隨坡度變化較小[21]。PM通過能夠有效考慮了陡峭區域的地形變化(圖1)，因此對坡度變化不敏感。

3.2.2 植被覆蓋度 基于參考數據中植被點云與原始點云點數量的比值，計算各樣本區域植被覆蓋度[19]，分析植被覆蓋對濾波結果的影響。結果表明，Ⅰ類誤差與植被覆蓋度呈顯著負相關(p<0.05)，而Ⅱ類誤差和總誤差與植被覆蓋度相關性不顯著(p>0.05)(圖5)。MCC和PM的Ⅰ類誤差隨植被覆蓋度上升下降速率最低，ETEW的Ⅰ類誤差隨植被覆蓋度增加降低速率最高，SBF和ATIN居中。植被覆蓋度較低85%的區域，MCC和PM的Ⅰ類誤差較ATIN小，而在植被覆蓋度大于85%時，ATIN的Ⅰ類誤差小于PM，但植被覆蓋度達到90%時，ATIN的Ⅰ類誤差與MCC接近。

圖5 植被覆蓋度與濾波算法誤差關系與坡度和植被覆蓋度的關系

Sithole等[17]等研究認為Ⅰ類誤差隨植被覆蓋度上升而增大，而本文得出與之相反的結論。本研究區域地形復雜，地形不僅直接影響濾波結果，而且顯著影響植被分布，植被覆蓋度高的區域比較平緩，Ⅰ類誤差較小。而植被覆蓋度低的區域，地形陡峭，坡度因素占主導，導致Ⅰ類誤差比較高。同時，地形對植被覆蓋的影響解釋了各算法的Ⅰ類誤差對植被覆蓋度變化的敏感性差異，也解釋了Ⅰ類誤差較Ⅱ類誤差對植被覆蓋度的變化更為敏感。

3.2.3 點云密度 對各樣本區域原始的點云數據進行抽稀，獲得不同抽稀密度(5%，10%，25%，50%，100%)的樣本數據，利用濾波算法對不同點云密度數據進行濾波，分析點云密度對不同濾波算法精度的影響。結果表明，隨著點云密度增加，MCC，ATIN，PM的Ⅰ類誤差無明顯變化，Ⅱ類誤差和總誤差緩慢下降，在5%和100%點云密度下的Ⅰ類誤差略高于10%～50%點云密度；SBF和ETEW的Ⅰ類誤差增加，Ⅱ類誤差和總誤差下降(圖6)。

圖6 濾波算法誤差隨點云密度的變化

通常情況下，低點云密度下裸地與植被不易分離，導致Ⅱ類誤差增大[22-23]，因此隨著點云密度增加各算法的Ⅱ類誤差都呈下降趨勢。SBF和ETEW的Ⅰ類誤差與點云密度正相關，原因是隨著點云密度的升高，固定鄰域內的點數增加，大于高程或坡度閾值而被歸于非地面點的數量增加，這部分點中包含真實地面點云，它們的剔除提升了Ⅰ類誤差。其次MCC，ATIN，PM誤差隨點云密度變化較為緩慢，原因可能是此3種濾波算法均需利用點云構建連續表面，連續表面的質量影響算法的精度[15]，本文獲取的原始的點云密度較高(平均密度26～56個/m2)，在其基礎上抽稀對連續表面質量影響不明顯。

4 結 論

本文以黃土高原溝壑區董莊溝為研究區，在流域內選取12個樣本區域，選用MCC，ETEW，ATIN，PM，SBF共5種濾波算法對獲取的無人機LiDAR點云數據進行濾波，對比分析不同算法的濾波精度及其與影響因素(坡度、植被覆蓋、點云密度)的關系。結果表明，MCC和PM的Ⅰ類誤差通常最低，ETEW和SBF最高，ATIN處于中間(偶爾小于MCC)，各算法的Ⅱ類誤差大小順序與Ⅰ類誤差大致相反。PM在陡峭區域的濾波效果最優，能夠較好地適應黃土高原的復雜地形地貌特征。坡度和植被覆蓋度主要影響各算法的Ⅰ類誤差，坡度與植被覆蓋度呈顯著負相關，各算法的Ⅰ類誤差隨坡度顯著上升，隨植被覆蓋度顯著下降。MCC和PM對坡度不敏感，Ⅰ類誤差隨坡度上升速率最小；ETEW受坡度因素影響最大，Ⅰ類誤差隨坡度上升速率最快；SBF和ATIN可一定程度緩解坡度對濾波結果的影響，Ⅰ類誤差上升速率居中。隨著點云密度增加，MCC，ATIN，PM的Ⅰ類誤差無明顯變化，Ⅱ類誤差和總誤差緩慢下降，SBF和ETEW的Ⅰ類誤差增加，Ⅱ類誤差和總誤差下降。