趙桂華，鄒曉亮，郭 麗

（1.西安測繪總站，陜西 西安 710054）

機載LiDAR點云數(shù)據(jù)自動生成DEM的方法與精度評價

趙桂華1，鄒曉亮1，郭 麗1

（1.西安測繪總站，陜西 西安 710054）

設計了一種點云數(shù)據(jù)快速處理自動生成DEM的算法，介紹了濾波結果的評價方法以及通過標準DEM評價內插的DEM整體精度的方法；并選用河北承德地區(qū)的機載激光LiDAR點云數(shù)據(jù)進行了實驗。結果表明，數(shù)據(jù)處理結果具有較高的精度，為激光點云數(shù)據(jù)自動生成DEM提供了一種有效的技術途徑。

LiDAR點云；DEM；濾波；DEM內插；濾波結果評價；精度評價

地球表面高低起伏呈現(xiàn)為一種連續(xù)變化， DEM[1]是定量描述地球表面地貌結構及空間變化的基礎數(shù)據(jù)，蘊含著大量豐富的地貌特征信息。從攝影測量技術的實現(xiàn)途徑來看，比較成熟的生成DEM 的方法包括：基于解析測圖儀的立體影像數(shù)字化掃描方法和基于全數(shù)字攝影測量立體環(huán)境下的人工采集方法。這兩 種方法生成的DEM，點位準確﹑內插精度高，能達到測繪產品需求的高精度要求，但成圖周期長，不利于大范圍﹑大面積的DEM 成果生產。近年來，出現(xiàn)了基于影像自動匹配生成數(shù)字表面模型（DSM）[1]，再通過濾波生成DEM的方法，如像素工廠[2-3]和全數(shù)字攝影測量系統(tǒng)，但這些設備硬件成本較高。隨著機載 LiDAR[4-5]技術的發(fā)展，大區(qū)域﹑大批量高效獲取激光點云數(shù)據(jù)成為可能，激光LiDAR技術成為了一種新的技術途徑。針對機載LiDAR點云數(shù)據(jù)，本文設計了一種從點云到DEM的數(shù)據(jù)快速處理算法，通過對點云數(shù)據(jù)進行濾波﹑內插自動生成DEM，再利用標準DEM數(shù)據(jù)對其進行整體精度評價，最后通過工程數(shù)據(jù)驗證了算法的可行性。

1 點云數(shù)據(jù)自動生成DEM算法的設計思路

機載LiDAR點云數(shù)據(jù)自動生成DEM算法的設計思想為：首先對LiDAR數(shù)據(jù)進行預處理剔除粗差點，再利用濾波算法對點云數(shù)據(jù)進行分類，獲取地面點和非地面點。本文選用TIN濾波﹑曲面濾波和坡度濾波3種方法對點云進行濾波分類，并對濾波結果進行評價。在計算第I類誤差和第II類誤差的基礎上，選擇濾波效果好的地面點數(shù)據(jù)進行內插計算，自動生成DEM。內插方法實現(xiàn)了反距離加權插值﹑趨勢面插值和Kriging插值的DEM自動生成，并將內插DEM與標準DEM進行了精度比對。該流程與常規(guī)DEM的成果采集相比，具有快速﹑高效的特點，如圖1所示。

圖1 機載LiDAR點云數(shù)據(jù)自動生成DEM算法流程圖

2 點云數(shù)據(jù)處理算法

2.1 點云數(shù)據(jù)預處理

激光LiDAR測量是一種無差別的測量方式，具有一定的盲目性和隨機性；對在其覆蓋區(qū)域內的所有目標（飛鳥﹑行人等）均進行測量，并盲目地記錄所有回波信號[5]，所以LiDAR系統(tǒng)獲取的點云數(shù)據(jù)存在一定量的粗差信息。數(shù)據(jù)預處理主要完成粗差剔除和地面種子點的選擇。粗差點通常與其鄰域數(shù)據(jù)點具有較大高差，一般可采用單一閾值法進行濾波處理，如選取地面高程的最大值作為粗差點剔除的閾值。為了精細化處理，可通過建立多級虛擬格網(wǎng)，結合種子點的選取，設定閾值和搜索窗口的半徑，對點云數(shù)據(jù)進行濾波。種子點的選擇是關鍵因素，因為初始地面種子點的選取直接影響數(shù)據(jù)點是地面點還是非地面點的判斷依據(jù)。結合多級虛擬格網(wǎng)技術，引入多尺度格網(wǎng)，篩選初始的地面種子點，為后續(xù)點云數(shù)據(jù)的濾波提供技術支撐。

2.2 點云數(shù)據(jù)濾波算法

點云數(shù)據(jù)濾波算法主要指基于相關信息建立若干判別規(guī)則，依據(jù)這些規(guī)則從點云數(shù)據(jù)中剔除非地面點，獲取地面點的過程。本文選擇了TIN濾波﹑曲面濾波和坡度濾波3種方法對點云數(shù)據(jù)進行濾波處理，并對濾波結果進行分析評價，選出濾波效果好的點云數(shù)據(jù)作為后續(xù)內插的原始數(shù)據(jù)。

1）TIN濾波。其算法原理為：先獲取一定的地面種子點組成初始的稀疏TIN；然后對各點進行判斷，若該點到三角面的垂直距離及角度小于設定的閾值，則將該點加入地面點集合，實現(xiàn)TIN的不斷加密，并重新計算TIN；再對非地面點集合內的點進行判斷；如此迭代，直到不再增加新的地面點或滿足給定條件為止。

從濾波原理上來看，每增加一個地面點，就要更新一次TIN；隨著地面點數(shù)量的增加，構建TIN 的計算量也迅速增加。為解決這一問題，本文采用虛擬格網(wǎng)鄰域搜索，對該算法進行改進：在獲得初始地面點的基礎上，對于非地面點集合中的某一點，搜索其周邊的10～20個地面點，構建局部TIN，并依據(jù)閾值對該數(shù)據(jù)點進行判斷，滿足條件就加入地面點集合，否則加入非地面點集合，直至完成所有非地面點的判斷；再基于新的地面點集合用相同的方法對新的非地面點集合再做一次篩選判斷，從而得到較準確的地面點集合和非地面點集合。

2）曲面濾波。其算法原理為：把地面看作是連續(xù)且平緩變化的表面，用帶有限制條件的參數(shù)曲面進行濾波約束。首先對離散的激光點云數(shù)據(jù)進行二維排序；再選取種子點區(qū)域，設定閾值，利用平面或曲面進行濾波；最后根據(jù)濾波結果不斷更新趨勢面，逐步完成整個測區(qū)的濾波過程。

3）坡度濾波。其算法原理為：地表的坡度變化比地物的坡度變化小，對于一個真實的地面點來說，與非地面點之間會有較大的坡度值。坡度閾值的設定是該算法的關鍵。本文結合多尺度虛擬格網(wǎng)技術和坡度的方向性，以一級格網(wǎng)為單元設定坡度閾值，待判斷點比地面點高的稱為上坡，比地面點低的稱為下坡。對于上坡的情況，由一級格網(wǎng)內的一級地面點與高程值最大的二級地面種子點之間的坡度值來確定；對于下坡的情況，一般取上坡坡度閾值的0.6倍作為下坡的坡度閾值；同時設定最大坡度閾值，防止地形突變引起閾值畸變。

2.3 濾波后點云數(shù)據(jù)的精度評價

根據(jù)國際攝影測量與遙感大會ISPRS小組在對經典濾波算法測評時提出的誤差評判標準[6]，將點云數(shù)據(jù)的誤分率分為3種：第I類誤差﹑第II類誤差和總誤差。第I類誤差是將地面點誤分為非地面點的誤差，又稱拒真誤差；第II類誤差是將非地面點誤分為地面點的誤差，又稱納偽誤差；第I類誤差和第II類誤差可以綜合得到總誤差。第I類誤差﹑第II類誤差反映了濾波算法的適應性，總誤差反映濾波算法的可行性。本文通過對點云濾波誤差的分析來評價點云濾波算法的性能。

2.4 濾波后點云數(shù)據(jù)的內插

DEM內插，即根據(jù)若干相鄰參考點的高程值求取待定點的高程值。任意一種內插方法都是基于原始地形起伏變化的連續(xù)光滑性，或者說鄰近的數(shù)據(jù)點間必須有很大的相關性，才可能由鄰近的數(shù)據(jù)點高程內插得到待定點的高程。本文采用反距離加權插值法﹑趨勢面插值法和Kriging插值法實現(xiàn)DEM的內插。

2.4.1 反距離加權插值法

反距離加權插值法[7]（IDW）基于相似相近原理，即兩個物體距離越近，其性質越相似。以待內插點為中心，在待內插的點X﹑Y 處內插變量Z 值時，在局部鄰域中計算各數(shù)據(jù)點的加權平均值，其公式為：

式中，Zi為估算點的數(shù)值；Z0為插值要素在第i 點的實測值；λi為預算過程中要使用的各種采樣點的權重；n為參與插值的實測點的數(shù)量。確定權重的公式為：

式中，P為冪指數(shù)；Di為第i點和各已知樣點之間的距離。采樣點在預測點值的計算過程中所占的權重受冪指數(shù)影響，即隨著采樣點與預測值之間距離的增加，采樣點對預測點影響的權重按指數(shù)規(guī)律減少。

2.4.2 趨勢面插值法

趨勢面插值法是利用多項式方程擬合已知點，用于估算其他點的高程值，屬于一種非精確的插值方法。該方法通過選擇一個二元函數(shù)來逼近采樣數(shù)據(jù)的整體變化趨勢。二元函數(shù)的一般形式為：

當p=0時，此二元函數(shù)表示水平面；當p=1時，此二元函數(shù)表示傾斜平面；當p=2時，此二元函數(shù)表示趨勢面，可用于模擬地形起伏。一次趨勢面需要估算3個系數(shù)，而3次趨勢面需要估算10個系數(shù)才能預測未知點的值，因此趨勢面模型次數(shù)越高，計算量越大。

2.4.3 Kriging插值法

Kriging[8-9]插值法也稱空間局部估計或空間局部插值，是地統(tǒng)計學中最經典的方法之一。在二階平穩(wěn)﹑本征假設等數(shù)學假設條件下，該方法是建立在變異函數(shù)理論及結構分析的基礎上，在有限區(qū)域內對區(qū)域化變量的取值進行無偏最優(yōu)估計的一種方法。本文實現(xiàn)了基于球狀模型的Kriging插值算法。球狀模型的一般公式為：

式中，C0為塊金常數(shù)；C0+C為基臺值；C為拱高； a為變程。當C0=0，C=1時，即標準球狀模型。

2.5 基于標準DEM的整體精度評價方法

為了對DEM 內插成果進行精度評價，本文提出了基于標準DEM的整體精度評價方法。整體精度評價是將標準DEM與內插的DEM結果進行整體比較，統(tǒng)計整體的誤差分布。標準DEM是利用 TerraSolid 軟件進行人工精細分類生成的DEM 數(shù)據(jù)。該評價方法的基本思路是從標準DEM 中隨機選取一定數(shù)量的點作為評價控制點，與內插的DEM進行比較，同時設定剔除誤差的閾值進行統(tǒng)計計算，統(tǒng)計點集的中誤差和平均誤差。分析濾波內插后點云數(shù)據(jù)中是否存在粗差點，若存在粗差點，則對其區(qū)域進行人工編輯，剔除內插有誤的點，再進行數(shù)據(jù)點集的誤差統(tǒng)計和分析。

3 實驗結果及精度分析

3.1 實驗數(shù)據(jù)

本文選取的實驗數(shù)據(jù)為河北省承德地區(qū)機載激光LiDAR ALS50傳感器獲取的離散激光點云，見圖2，激光點云間隔為1.2 m，選取677 m×686 m的區(qū)域，數(shù)據(jù)點共計28 880個；最小高程值為325.33 m，最大高程值為434.06 m；影像數(shù)據(jù)的分辨率為0.25 m，見圖3。

圖2 點云數(shù)據(jù)圖

圖3 影像數(shù)據(jù)圖

3.2 點云濾波

為了使濾波結果顯示得更直觀，采用濾波后的點云數(shù)據(jù)與對應影像數(shù)據(jù)疊加的方式來顯示濾波后的點云數(shù)據(jù)。利用TIN濾波﹑曲面濾波和坡度濾波3種方法濾波后所得到的地面點與對應影像的疊加情況，如圖 4～6所示。

圖4 TIN濾波

圖5 坡度濾波

圖6 曲面濾波

3.3 濾波結果評價

濾波的目的是為了最大限度地恢復真實地形，濾波后得到的地面點一般要求包含盡可能少的非地面點，即要盡可能地減小第II類誤差的發(fā)生。實驗中主要考慮第II類誤差的影響，并綜合考慮總誤差的影響。表 1反映了利用3種濾波算法對實驗數(shù)據(jù)濾波后的第I類誤差﹑第II類誤差和總誤差分布情況。

表1 不同濾波方法的誤差率/%

由表1可知，TIN濾波的第II類誤差率為1.42%，效果最好，但是總誤差率為20.05%；坡度濾波和一次曲面濾波的第II類誤差控制得較好，但坡度濾波的總誤差大于一次曲面濾波的總誤差；同時采用目視判讀的方法，將濾波后的點云數(shù)據(jù)與影像數(shù)據(jù)疊加進行判斷，目視結果也是一次曲面濾波優(yōu)于坡度濾波。

3.4 濾波后點云內插

根據(jù)濾波結果的誤差率分析，本次實驗內插前的點云數(shù)據(jù)選取曲面濾波算法所獲得的數(shù)據(jù)，如圖7所示。分別采用IDW插值﹑趨勢面插值和Kriging插值3種算法對內插前的地面點數(shù)據(jù)進行處理，其內插后的結果如圖8～10所示，形成內插DEM。

圖7 內插前的點云數(shù)據(jù)

3.5 整體精度評價

將濾波內插的DEM成果數(shù)據(jù)與標準DEM成果數(shù)據(jù)進行逐點比較，計算中誤差和平均誤差。在剔除大于7 m粗差點的基礎上，IDW內插的中誤差為0.61 m，平均誤差為0.30 m；趨勢面內插的中誤差為0.99 m，平均誤差為0.51 m；Kriging內插的中誤差為0.55 m，平均誤差為0.28 m。從整體評價的結果來看，IDW插值和Kriging插值方法的中誤差和平均誤差的趨勢較一致，而趨勢面插值法的誤差值相對較大?；贙riging插值的總體評價如圖11所示?？偟膩碚f，相較于3種方法，Kriging內插得到的DEM結果最優(yōu)。最終自動生成的DEM成果選取Kriging內插的DEM結果。

圖8 IDW內插

圖9 Kriging內插

圖10 趨勢面內插

圖11 Kriging內插總體評價

4 結 語

本文設計了機載 LiDAR點云數(shù)據(jù)快速處理生成DEM的流程，提出了通過標準DEM評價基于LiDAR點云濾波內插生成的DEM整體精度的方法，并利用機載LiDAR點云工程數(shù)據(jù)進行了實驗驗證。結果表明，Kriging插值算法在快速自動生成DEM方面精度高，處理流程總體設計合理，整體精度評價方法具有可操作性。該方法是行之有效的。

[1] 全曉萍,宋志勇.基于DEM 數(shù)據(jù)的地形分析研究與實現(xiàn)[J].科技信息,2007(28):345-349

[2] 曹敏,史照良.新一代海量影像自動處理系統(tǒng)“像素工廠”初探[J].測繪通報,2006(10):55-58

[3] 鄒曉亮,繆劍,張永生,等.基于像素工廠的無人機影像空三優(yōu)化技術[J]測繪科學技術學報,2012,29(5):362-367

[4] 余潔,張國寧,秦昆,等. LiDAR數(shù)據(jù)的過濾方法探討[J].地理空間信息,2006,4(4):8-10

[5] 賴旭東.機載激光雷達基礎原理與應用[M].北京:電子工業(yè)出版社,2010:37-99

[6] Sithole G, Vosselman G. Experimental Comparison of Filter Algorithms for Bare-earth Extraction from Airborne Laser Scanning Point Clouds[J]. ISPRS Journal of Photogrammetry & Remote Sensing,2004,59(1/2):85-101

[7] LIU X Y, GUO Q H, SU Y J, et al. Airborne LiDAR for DEM Generation: Some Critical Issues[J].Progress in Physical Geography,2008,32(1):31-49

[8] Vieira S R,Hatfield J I, Nielsen D R, et al. Geostatistical Theory and Application to Variability of Some Agronomical Properties[J]. Hilgardia,1983, doi:10.3733/hilg.v51n03p075

[9] 王政權.地統(tǒng)計學及在生態(tài)學中的應用[M].北京:科學出版社,1999:102-132

P23

B

1672-4623（2017）09-0009-04

10.3969/j.issn.1672-4623.2017.09.003

2016-07-15。

項目來源：西安科技創(chuàng)新工作站資助項目（2014-CX0452）。

趙桂華，碩士研究生，工程師，研究方向為激光LiDAR數(shù)據(jù)處理、數(shù)字攝影測量、影像分析與解譯。