陳傳法，王夢櫻，楊 帥，王 珍

(1.山東科技大學(xué) 測繪與空間信息學(xué)院，山東 青島 266590；2.北京市勘察設(shè)計研究院有限公司，北京 100038)

隨著激光雷達(dá)技術(shù)的發(fā)展，機(jī)載激光雷達(dá)(light detection and ranging，LiDAR)已成為采集地面點信息的主流工具[1]。與傳統(tǒng)攝影測量技術(shù)相比，LiDAR能夠快速獲取高精度、高分辨率的真實地面三維坐標(biāo)數(shù)據(jù)，已廣泛應(yīng)用于數(shù)字高程模型(digital elevation model，DEM)構(gòu)建、水文分析[2]、森林清查[3]、城市三維可視化[4-5]、土地覆蓋和土地利用分類等[6]。其中，LiDAR點云使用前一般需要將地面點從原始點云數(shù)據(jù)中分離出來，即濾波[7]。然而，原始LiDAR點云數(shù)據(jù)中通常包含許多建筑物、植被等非地面點，特別是對于地形起伏變化較大的林區(qū)，分類地面點與非地面點仍是一項挑戰(zhàn)。

目前，許多學(xué)者對LiDAR數(shù)據(jù)濾波展開了一系列研究，并提出了相應(yīng)的濾波算法。根據(jù)濾波原理，現(xiàn)有濾波算法可大致分為：形態(tài)學(xué)濾波[8-9]、坡度濾波[10]、分割濾波[11-12]和內(nèi)插濾波[13]。實際應(yīng)用中各種濾波算法各有優(yōu)缺點，其中插值濾波方法因原理清晰、能更好地模擬真實地表以及濾波可靠性高等優(yōu)勢被廣泛關(guān)注。例如，Axelsson[14]提出漸進(jìn)不規(guī)則三角網(wǎng)(triangulated irregular network，TIN)加密濾波方法(progressive TIN densification，PTD)，根據(jù)局部低點構(gòu)造初始的稀疏TIN，并根據(jù)每個點與TIN的距離和角度判斷地面點。Haugerud等[15]提出一種虛擬森林砍伐算法，用原始點云構(gòu)建TIN，然后計算采樣點及其對應(yīng)表面的曲率。Evans等[16]提出一種多尺度曲率分類方法，該方法采用正規(guī)化樣條代替TIN插值并在不同分辨率下構(gòu)建柵格曲面。Chen等[17]提出多分辨率分層點云濾波算法(multi-resolution hierarchical classification，MHC)，用局部最低點作為地面種子點，通過薄板樣條函數(shù)(thin plate spline，TPS)構(gòu)造地面參考面。分析表明，基于插值的濾波在地形相對平坦、地形較為簡單的地區(qū)表現(xiàn)良好，但在地形復(fù)雜林區(qū)由于獲取的初始地面點少，插值迭代計算受初始DEM影響增大且誤差累積增大，其濾波效果并不理想。為此，提出一種適用于林區(qū)機(jī)載LiDAR點云的多分辨率層次插值濾波方法，該方法首先借助形態(tài)學(xué)濾波和穩(wěn)健的標(biāo)準(zhǔn)分?jǐn)?shù)(z-score)方法提高地面種子點數(shù)量和質(zhì)量，以提高初始地面參考面的精度；然后從低層到高層濾波過程中，借助自適應(yīng)坡度閾值選擇地面點，以保證濾波在地形崎嶇林區(qū)的適應(yīng)性。

圖1 多分辨率層次插值濾波算法流程圖

1 算法原理及流程

首先格網(wǎng)化原始點云，并使用形態(tài)學(xué)迭代開運算獲取潛在地面種子點；然后借助穩(wěn)健“z-score”準(zhǔn)則剔除潛在地面種子點中的非地面點；最后分三層濾波，每一層用薄板樣條內(nèi)插構(gòu)建一定網(wǎng)格分辨率的地面參考面，并通過待分類點與參考面的距離判斷該點是否為地面點；濾波出的地面點用于更新參考面，直到該層沒有地面點加入。算法的流程如圖1所示，其中涉及到的重點環(huán)節(jié)介紹如下。

1.1 初始地面種子點的選取

1) 點云格網(wǎng)化。以c作為格網(wǎng)大小(例如c=1 m)，落入每個網(wǎng)格內(nèi)的最低點作為對應(yīng)的格網(wǎng)值(圖2)。格網(wǎng)大小取決于點云密度，一般取平均點間距。對于剩余的未被填充的格網(wǎng)，利用彈簧質(zhì)子圖像修復(fù)技術(shù)填補(bǔ)，修復(fù)后形成的最小值表面作為迭代開運算的初始表面。

圖2 點云格網(wǎng)化

2) 迭代開運算。開運算為先腐蝕后膨脹。開運算結(jié)構(gòu)元素的選取直接影響運算的質(zhì)量，因此選擇各向同性并且不會對運算結(jié)果產(chǎn)生偏移的圓盤形作為開運算結(jié)構(gòu)元素的形狀，同時采用多尺度結(jié)構(gòu)元素以避免結(jié)構(gòu)元素過大或過小造成的誤分類。為此，定義兩個參數(shù)：最大窗口半徑(Wmax) 和坡度值(s)。其中，由Wmax形成一個動態(tài)變化窗口的開運算器，并計算用于分類格網(wǎng)的高程閾值：

ET=s×Wi。

(1)

其中：Wi=1,2,…,Wmax，表示每次迭代開運算的窗口大小；ET表示每次迭代開運算前后兩個表面之間對應(yīng)格網(wǎng)的差值閾值。如果開運算前后格網(wǎng)高差閾值大于ET，則標(biāo)記該格網(wǎng)為非地面格網(wǎng)。每次開運算后得到的表面被當(dāng)作下一次開運算的初始表面(圖3)，迭代全部完成后形成一個新的最小值表面。由于選擇最低點時可能會選中極低的無意義的點，這會對開運算結(jié)果產(chǎn)生影響。因此，要翻轉(zhuǎn)最后得到的表面并用單一的窗口半徑再進(jìn)行一次開運算。最后，將未標(biāo)記為非地面格網(wǎng)中的點作為潛在的地面種子點。

圖3 迭代開運算示意圖

1.2 種子點粗差剔除

種子點粗差剔除旨在剔除潛在地面種子點集中的非地面點，以提高初始構(gòu)建的地面參考面精度。異常值檢測的方法有許多，本研究擬基于z-scores來檢驗異常值，即

(2)

(3)

式中，

其中，常數(shù)1.483是使MAD正態(tài)分布無偏的校正因子；median為取一組數(shù)據(jù)的中位數(shù)。

1.3 地面點選擇

在獲取地面種子點之后，通過多分辨率插值濾波選擇剩余采樣點中的地面點。該方法包括三個層次，其中參考DEM分辨率和殘差閾值從低到高逐漸增大。在第i層(i= 1,2,3…)地面點選取步驟如下:

總輻射年最大值為1 112.I7 W/m2：觀測期間凈輻射年平均 48.24 W/m2，春季 65.01 W/m2，夏季93.79 W/m2，秋季 30.39 W/m2，冬季 2.42 W/m2，最大 807.98 W/m2，最小 186.99 W/m2。

1) 利用TPS對地面種子點插值生成分辨率為h的參考DEM。

2) 計算每個采樣點與其對應(yīng)DEM網(wǎng)格點所在3×3移動網(wǎng)格中9個點的殘差。如果9個殘差中至少有5個小于設(shè)定閾值，則該點判斷為地面點。經(jīng)典MHC方法僅考慮高程差，因此在陡坡和地形變化劇烈的區(qū)域容易導(dǎo)致點的錯分。為此，本研究在殘差閾值中引入坡度因子，以提高濾波方法對各種地形區(qū)域的自適應(yīng)性。每個網(wǎng)格的殘差閾值T計算如下：

T=ti+e×SF，

(4)

其中，ti為每一層中固定的閾值(i= 1,2,3)；e為尺度縮放因子；SF為每個點云對應(yīng)的DEM上的近似坡度，計算公式為：

(5)

其中，F(xiàn)是X、Y和Z之間的函數(shù)，Z=F(X,Y)。

如圖4所示，在陡坡上固定閾值易將地面點誤分類為非地面點，而自適應(yīng)坡度閾值可以減少誤分類。

圖4 適應(yīng)坡度閾值和固定閾值分類點示意圖

3) 更新地面點和地面種子點。地面種子點是分辨率為h的網(wǎng)格單元中局部最低點。

4) 重復(fù)步驟1)～4)直到?jīng)]有地面點加入。

2 實驗及結(jié)果分析

分別采用基準(zhǔn)數(shù)據(jù)和林區(qū)數(shù)據(jù)驗證新方法的濾波精度，并將結(jié)果與常用濾波方法比較。采用的濾波精度評價指標(biāo)為I類誤差(T.I)、II類誤差(T.II)、總誤差(T.E)和Kappa系數(shù)，計算公式如下：

T.I=b/(a+b)×100%，T.II=c/(c+d)×100%，T.E=(b+c)/e×100%，κ=(p0-pc)/(1-pc)×100%。

(6)

其中：a表示正確分類的地面點數(shù)；b表示誤分為非地面點的地面點數(shù)；c表示誤分為地面點的非地面點數(shù)；d表示正確分類的非地面點數(shù)；e=a+b+c+d；p0=(a+d)/e；pc=((a+b)×(a+c)+(c+d)×(b+d))/e2。

I類誤差表示地面點錯分為非地面點比例，II類誤差表示非地面點錯分為地面點比例，總誤差為衡量兩類誤差的綜合指標(biāo)，I類、II類誤差和總誤差越小以及Kappa系數(shù)越大表示濾波精度越高。

2.1 ISPRS基準(zhǔn)數(shù)據(jù)

采用國際攝影測量與遙感協(xié)會(International Society for Photogrammetry and Remote Sensing，ISPRS)提供的6個山地區(qū)域(samp51-samp71)樣本來測試新算法性能。其中，每個樣本的點云數(shù)據(jù)都通過半自動過濾和人工編輯分類出地面點和非地面點。山地區(qū)域一般地勢復(fù)雜且有植被覆蓋，因此新算法將其視為林區(qū)進(jìn)行濾波處理。為了測試算法的最佳性能，通過重復(fù)試驗確定新方法的濾波參數(shù)如表1所示，其中初始分辨率h均為2 m。為了增加算法在各類地形的適應(yīng)性，規(guī)定最優(yōu)參數(shù)設(shè)置準(zhǔn)則：最大窗口半徑(Wmax)應(yīng)稍大于數(shù)據(jù)處理區(qū)域最大非地面物尺寸；初始高差閾值(t)應(yīng)稍小于低矮植被高度；坡度(s)為0.2；坡度系數(shù)(e)可設(shè)置為1。根據(jù)此準(zhǔn)則給出單參數(shù)集(Wmax=13 m，s=0.2，t=0.29 m，e=1)及其濾波結(jié)果如表1所示。

表1 新方法濾波誤差及處理時間

由表1可得，單參數(shù)集取得的樣本平均Kappa系數(shù)和總誤差分別為82.31%和3.17%，最優(yōu)參數(shù)取得的樣本平均Kappa系數(shù)和總誤差分別為87.88%和1.89%。單參數(shù)與最優(yōu)參數(shù)的總誤差僅相差約2%，說明該方法對參數(shù)設(shè)置不敏感。由處理時間可以看出，該方法運算效率整體較好，運算效率與點數(shù)相關(guān)，點數(shù)越少，處理時間越短。表1中samp51的Kappa系數(shù)最大，高達(dá)95.73%，samp61的總誤差最小，僅為0.74%。samp52的總誤差最大，samp53的Kappa系數(shù)最小。圖5和圖6展示了這兩個樣本濾波前后DEM的對比以及I類和II誤差的空間分布情況。由圖5(a)可見，samp52中斷的陡坡、沿河低矮植被及河岸高低變化的地形造成了II類誤差增大；分布在河道兩側(cè)的錯分的地面點(圖5(c))造成了濾波后DEM中河道兩側(cè)凸起問題。由圖6(a)可見，samp53具有不連續(xù)和陡峭的梯田，其中西南角地形變化明顯，導(dǎo)致II類誤差較大(圖6(c))。

(空心圓點表示Ⅰ類誤差；實心圓點表示Ⅱ類誤差)

將新方法與近5年提出的10種濾波算法[17-26]進(jìn)行比較見表2，可見：新方法在6個樣本中有4個總誤差最小，樣本samp52和samp71的總誤差與最佳結(jié)果接近，僅比相應(yīng)最佳濾波方法誤差高0.07%和0.30%。新方法的平均總誤差為1.89%，明顯優(yōu)于其他濾波方法。

表2 新方法與10種濾波算法總誤差對比

2.2 林地區(qū)域數(shù)據(jù)

本研究從美國國家科學(xué)基金會提供的開放地形數(shù)據(jù)網(wǎng)站(http://www.opentopography.org/)下載了6個具有不同地形特征和植被覆蓋的樣本(圖7)，以進(jìn)一步驗證新方法在不同景觀下的適應(yīng)性。這6組樣本均按照1∶1 000標(biāo)準(zhǔn)圖幅大小取樣。數(shù)據(jù)1為美國猶他州中南部垂直懸崖樣本，地形陡峭有較多陡坡及懸崖，覆蓋植被低矮且稀疏(圖7(a))。數(shù)據(jù)2(圖7(b))和3(圖7(c))為美國科羅拉多州緩慢移動的滑坡樣本；數(shù)據(jù)2地形平緩，覆蓋植被高且分布較密集，數(shù)據(jù)3地形較平緩，包含較多褶皺，覆蓋植被高且分布不均勻。數(shù)據(jù)4(圖7(d))為新西蘭惠靈頓森林樣本，地形陡峭，覆蓋植被高且密集。數(shù)據(jù)5(圖7(e))和數(shù)據(jù)6(圖7(f))為加利福尼亞州猛犸湖熔巖表面樣本，包含高低起伏較大的山包，覆蓋植被較高且分布較密集。

將新方法與形態(tài)濾波算法(morphological filtering algorithm，MF)[27]和漸近不規(guī)則三角網(wǎng)加密濾波算法(PTD)[14]濾波結(jié)果比較(表3)可知：新方法的總誤差低于其他方法，Kappa系數(shù)最大。新方法的濾波精度最高，平均總誤差為6.82%，Kappa系數(shù)為85.61%，其次是PTD，MF濾波結(jié)果最差。因此，與經(jīng)典濾波方法相比，新方法更適合于森林地區(qū)。

圖7 6組樣本數(shù)據(jù)衛(wèi)星影像

表3 新方法與兩種常用濾波算法誤差對比

三種方法構(gòu)建的DEM精度表明(表4)：新方法生成的DEM精度整體較好，6組數(shù)據(jù)中有5組數(shù)據(jù)誤差最小。數(shù)據(jù)2的DEM中誤差大于其他方法，這是由于新方法II類誤差較大(表3)，即濾波結(jié)果中非地面點較多引起的。數(shù)據(jù)集中各方法對數(shù)據(jù)2構(gòu)建的DEM中誤差最小，這是由于在平緩和較光滑的地形上各類濾波算法運行良好；各方法在數(shù)據(jù)4構(gòu)建的DEM中誤差最大，這是由于太過密集的植被易被認(rèn)為是地面點造成的。由表4可知，本方法獲得的DEM平均精度優(yōu)于其他方法。

表4 各個區(qū)域DEM中誤差對比

以數(shù)據(jù)1為例，圖8展示了各個方法構(gòu)建的DEM。結(jié)果顯示：新方法生成的DEM與標(biāo)準(zhǔn)DEM相似，而MF和PTD在地形復(fù)雜的區(qū)域會保留較多植被點造成DEM表面粗糙。

圖8 數(shù)據(jù)1標(biāo)準(zhǔn)DEM以及新方法、MF和PTD獲取的DEM

3 結(jié)論

為了提升林區(qū)點云濾波算法精度，本研究提出一種多分辨率層次插值濾波方法。該方法借助形態(tài)學(xué)開運算和穩(wěn)健z-score方法來精確地獲取大量地面種子點，并提出了一種考慮地形坡度的自適應(yīng)殘差閾值。借助ISPRS山區(qū)數(shù)據(jù)分析表明，該方法濾波平均總誤差和Kappa系數(shù)分別為1.89%和87.88%；對森林?jǐn)?shù)據(jù)集結(jié)果分析表明，該方法平均總誤差為6.80%，Kappa系數(shù)為85.61%。而且新方法精度均明顯優(yōu)于常用的插值濾波方法。然而，本研究構(gòu)建的濾波方法濾波結(jié)果中II誤差明顯大于I類誤差，這將導(dǎo)致濾波后的地面點中含有大量非地面點，進(jìn)而影響DEM建模精度。因此后續(xù)研究將借助對II類誤差有抵抗力的濾波算法進(jìn)一步提升本文算法的精度。