李 豆, 李朋飛, 穆興民, 姚頑強(qiáng), 湯伏全, 李 婷

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

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

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

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

1 研究區(qū)概況

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

2 數(shù)據(jù)源與研究方法

2.1 數(shù)據(jù)獲取

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

圖1 董莊溝及12個樣本區(qū)域的位置

表1 各樣本區(qū)域信息

2.2 數(shù)據(jù)處理與分析

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

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

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

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

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

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

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

(1)

(2)

(3)

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

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

表2 各濾波算法主要參數(shù)的最優(yōu)取值

3 結(jié)果與分析

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

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

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

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

圖3 5種濾波算法在樣本1區(qū)域濾波結(jié)果

3.2 影響濾波精度的因素

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

圖4 坡度與濾波算法誤差的關(guān)系

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

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

圖5 植被覆蓋度與濾波算法誤差關(guān)系與坡度和植被覆蓋度的關(guān)系

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

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

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

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

4 結(jié) 論

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