單光子激光雷達數據去噪與濾波算法

李 凱，張永生，童曉沖，楊偉銘，董 鵬

（1.軍事科學院系統工程研究院后勤科學與技術研究所，北京100071；2.戰略支援部隊信息工程大學地理空間信息學院，鄭州450001；3.北京航天控制儀器研究所，北京100039）

0 引言

20世紀90年代早期，第一臺商用機載激光雷達系統TopScan研制成功。經過20余年的發展，激光雷達技術已經被廣泛應用于測繪領域，成為地理空間信息獲取的一種重要手段。然而，隨著傳統線性激光雷達技術的成熟，其測繪能力的提升空間也越來越小。線性激光雷達一般采用較大的發射脈沖能量以及較低的飛行高度（較低的條帶寬度）以獲取足夠的激光脈沖回波能量，但這都導致獲取激光雷達數據的成本和時間代價增加［1］。因此，當使用線性激光雷達針對較大面積區域進行測繪任務或對同一區域進行長期動態監測時，時間和成本代價會限制傳統激光雷達的應用能力。

近些年，單光子激光雷達技術開始漸漸成熟［2-4］。不同于線性激光雷達在探測回波信號時需要數以百計的光子形成弱光返回信號，單光子激光雷達采用單光子敏感探測器，也就是說，每發射一個激光脈沖，新型的激光雷達只需要探測到若干光子的返回信息，就可以獲取被探測物體的位置信息。這使得激光光源產生的光子得到更有效地利用，大大減少了激光發射脈沖所需能量，同時使得激光雷達的作業高度大幅提高。這些技術優勢也使單光子激光雷達更適用于衛星激光測高任務，美國2018年9月15日發射的新一代監測冰蓋和海平面變化的ICESat-2衛星所搭載的激光高度計ATLAS就放棄了第一代傳感器GLAS所采用的全波形記錄技術，而采用了單光子探測技術。GLAS激光發射頻率僅為40Hz，而ATLAS激光脈沖發射頻率設計為10kHz，這使沿軌方向地面激光足印的間距由172m左右減小至約70cm，極大地提高了ATLAS沿軌方向的數據采集密度［5］。為了驗證ICESat-2搭載傳感器的可用性并開發適用于單光子激光雷達數據的處理算法，NASA研制了機載單光子激光雷達系統MABEL并向全球用戶免費提供其數據產品［6］。因此，通過對MABEL數據進行研究，能夠為機載和星載單光子激光雷達技術的應用提供理論依據。

單光子探測技術在采集數據時具有一定的優勢，但是其極高的探測靈敏度同時也給點云數據帶來了較大的背景噪聲［7］。因此，單光子激光雷達的數據去噪問題是獲取精確地形首先需要解決的問題。本文提出一種改進的基于局部距離統計的點云去噪算法初步去除了單光子激光雷達噪聲信息，從剩余的點云數據中進一步利用統計分析方法提取地面點，同時進一步去除地面以下噪聲點，并以傳統機載激光雷達的數字地形模型（Digital Terrain Model，DTM）作為參考數據，對從單光子激光雷達數據中提取的地面點云和生成的數字高程模型（Digital Elevation Model，DEM）數據進行對比驗證。

1 單光子激光雷達點云去噪和濾波模型

1.1 實驗數據介紹

（1）MABEL 數據

MABEL是搭載在NASA ER-2飛機上的單光子激光雷達設備，MABEL的飛行高度約為20km，飛行速度在200m／s左右，MABEL的激光脈沖重復頻率可在5kHz～25kHz范圍內調整。因此，在標稱飛行條件下，若采用10kHz的脈沖重復頻率，沿軌方向激光腳點間距為2cm。此外，MABEL利用發射器光纖分束器將一個發射脈沖分離為8個近紅外（1064nm）波束和 16個綠色（532nm）波束，且可以根據測量需求在-3°～＋3°范圍內調整發射波束的下視角［8］。結合飛機上搭載的GNSS和慣性導航系統，MABEL在工作狀態下能夠記錄每個返回光子的往返時間和三維位置信息。MABEL點云數據處理結果分段存儲于HDF格式文件中，每個文件包含60s飛行數據。其中，由每個子波束返回信息得到的點云數據存儲在一個通道（channal）中。本文選擇MABEL在North Carolina中北部區域的一個完整數據段進行處理，該數據在位于近紅外波段的44通道返回的光子數量最多，而噪聲信息也相對較多。因此，實驗中以44通道的點云數據為例檢驗本文提出的算法并利用參考數據進行驗證。MABEL點云數據為剖面點云，在進行去噪處理時通常在二維空間內進行，即將平面坐標轉換為沿軌距離，縱坐標為點云高程信息。本文使用的MABEL點云數據如圖1所示。

圖1 MABEL剖面點云數據Fig.1 Diagram of MABEL point cloud data

（2）傳統機載激光雷達數據

傳統機載激光雷達數據來源于G-LiHT（Goddard's LiDAR， Hyperspectral ＆ Thermal Imager）系統。單個的傳感器獲得的數據在與其他數據融合時，由于坐標基準等不一致，融合存在較大困難，而G-LiHT則將多種傳感器集成在一起，減輕了數據融合處理的難度。該系統集激光雷達、成像光譜儀、熱成像儀、GPS-INS于一體，以同時獲取植被結構、地物光譜、表面溫度等數據，并通過這些數據之間的融合處理為大尺度生態系統的科學研究提供高分辨率（＜1m）觀測數據［9］。 本文使用的數據獲取時間為2011年8月，其中激光雷達點云由Riegl VQ-480獲取。G-LiHT還提供了由點云數據處理得到的1m分辨率DTM數據，本文即利用該數據與MABEL數據處理結果進行對比。在UTM（Universal Transverse Mercator）17N投影帶中，DTM數據與MABEL數據的覆蓋范圍如圖2所示，該圖的左上角坐標為（650360.0000E， 3997988.0000N）， 右下角坐標為（652768.0000E， 3997988.0000N）。 圖 2中， 白色區域為DTM數據，里面的黑色實線為MABEL點云數據運行軌跡的投影。與MABEL數據段作對比的DTM數據由5個文件拼接而成，故相鄰文件之間存在少量的數據缺失。

圖2 DTM數據與MABEL數據的疊加顯示Fig.2 Overlay diagram of DTM data and MABEL data

1.2 MABEL點云去噪和濾波算法

（1）改進的基于局部距離統計的點云初步去噪算法

當前，針對MABEL點云數據的去噪算法主要可以分為兩類。一類將MABEL點云數據柵格化為數字影像，再通過邊緣檢測等數字影像處理技術分離噪聲點［10-11］。這種算法在點云數據柵格化過程中造成信息損失，因此使去噪結果精度降低［12］。另外，點云柵格化形成的圖像尺寸很大，在地面上一小段距離就能生成上百萬像素的圖像，這使得數據處理的計算量和時間增加［13］。另一類則根據信號和噪聲光子的分布特征，直接對點云數據進行處理。這種算法利用局部點云之間的統計量作為判斷閾值，較常用的統計量有局部點云的密度、模、均值、方差等［13-15］。該算法能夠剔除大部分的背景噪聲，保留較多的點云信息。但是，隨著統計范圍的增加，不同區域地形坡度、點云密度、地物形態的變化也越顯著，全局閾值的確定也越來越困難，噪聲光子難以完全消除，尤其是誤分類的真實地表下的噪聲點云使提取的地形產生較大誤差。

文獻［12］依據噪聲點和非噪聲點的局部密度差異，提出了一種基于局部距離統計的點云去噪算法。該算法計算點云中每個點到最臨近k（k＝50）個點的局部距離之和，并根據局部距離之和的頻數直方圖擬合Gauss函數，計算Gauss函數的均值和方差，以均值與t（t＝2）倍標準差之和作為閾值去除噪聲點。將該算法應用于小段MABEL實驗數據，實驗結果較好，算法有效剔除了大部分噪聲光子，同時能夠處理具有一定地形坡度的點云信息，且該算法計算較為簡單。但是，該算法在實際實驗中保留了少量地面以下的噪聲信息，在濾波處理時，這些點都是局部最低點，將這些點全部作為地面種子點擬合地形曲線，會給濾波結果帶來較大誤差。此外，在進行實驗時僅選取MABEL在532nm通道的一小段完全被植被覆蓋的區域，該段區域地表下的噪聲光子數較少，且地形特征較為單一。實際情況下，在完整的MABEL數據段中既有植被覆蓋區域，也有裸露的地表。而在不同的區域，非噪聲點云的密度也不相同，在點云稀疏的裸露地表區域、植被冠層頂部以及植被冠層內部非噪聲點云的密度甚至會小于噪聲點云。這就會使得在選擇較大閾值以保留所有非噪聲點情況下，更多的噪聲點被分類為非噪聲點，選擇較小閾值又會使非噪聲點被錯誤分類為噪聲點。因此，本文使用噪聲信息更加豐富的完整MABEL數據段檢驗該算法和本文提出的改進算法。

一般來說，在MABEL點云信息中，非噪聲點云在水平方向（沿軌方向）的聚集程度明顯高于垂直方向（如圖1所示）。因此，在使用局部距離統計作為閾值判定某點是否為噪聲點時，本文給水平方向的距離增加一個適當的權值因子ρ，使水平方向點云信號聚合度更高。這樣，在原始距離統計算法中使用的Euclidean距離表達式就變為

式（1）中，P、Q為點云中任意兩點， （xP，yp）、（xQ，yQ）分別為P、Q在水平方向沿軌距離和垂直方向上的高程。在本文實驗中，ρ取0.1。通過加入權值因子，使沿水平方向或接近水平方向的高密度點云更可能被分類為非噪聲點云。依據式（1），MABEL點云初步去噪算法的步驟如下：

步驟1：建立 MABEL點云數據的 K-D樹索引，搜索距離每個點最近的k個點。

步驟2： 計算第i點Pi到距離最近k個點Qj（j＝1，2，…，k）的“距離”之和的平均值Di， 距離計算公式如下

步驟3：以0.1m為距離間隔統計平均距離Di的頻數，生成頻數直方圖。

步驟4：以頻數直方圖峰值位置對應的距離作為均值，以峰值位置對應距離與最小平均距離的差值作為標準差，用Gauss函數擬合頻數直方圖，閾值設置為均值與t倍標準差之和，保留平均距離小于閾值的點為初選信號點，其余點為噪聲點。

（2）基于統計分析的點云濾波算法

提取去噪后點云中地面點的方法通常是在一定的窗口范圍內查找點云高程的最低點，以最低點作為地面種子點，再通過曲線擬合方法選取地面點［12，16］。在實際情況下，點云經過去噪處理后，地面以下通常還會含有一定量的噪聲信息，而選取地面種子點時由于這些點低于地面，往往被誤選為地面種子點。利用這些點進行曲面擬合會造成植被和地面點云的錯誤分類，且使提取的地形低于真實地形。因此，需要在濾波的過程中進一步剔除地面以下的噪聲信息。迭代曲線擬合是一種可選的方案，然而地面種子點分布本來就較為稀疏，迭代過程又進一步造成數據的抽稀，這會造成擬合結果精度降低，同樣影響分類結果。此外，該算法在迭代過程中閾值的選取也較為困難。

Horan等［13］引入統計分析的方法剔除了Greenland島冰川區域采集的單光子激光雷達數據中的噪聲信息，該方法認為單光子激光雷達點云的分布近似于Dirac delta函數，而該函數中的極大值即對應地面點云的分布。實現該方法的基本思路為：在沿軌方向上以一定的間隔將點云分割，對于每一小段點云，將點云的高程四舍五入至整數米，并計算取整后該段點云的模（頻數最大的高程值），則與模的高差在設定經驗閾值內的點被選取為地面點。然而，將該方法直接應用于植被覆蓋區域的點云，存在以下問題：1）在某一分段區域可能不存在信號點云，則計算得到該段點云的模值為某一噪聲點的高程，這可能嚴重偏離實際值，造成較大粗差。例如以20m為間隔分段時，在圖1中沿軌方向760m～780m不存在信號點云，該段點云模值為1134m，遠遠超出了實際高程。2）在樹木茂密區域，林下地面返回的信號點云稀疏，而樹冠冠層部分點云密度較高，此時計算得到該段點云模值位于植被冠層區域，也造成較大粗差，這種情況如圖3所示。圖3為從圖1中截取的從5980m～6200m的點云數據，圖中綠色點為利用文獻［13］中統計分析方法提取出的點云，圖中縱坐標顯示范圍僅包含了非噪聲點云部分，且點云均放大顯示。從圖3中明顯可以看出，這些被選出的點云雖然高程分布的聚集度較大，但并非地面點，而是樹冠部分點云。而地面被樹冠遮擋，獲得的點云信息很少。

圖3 使用統計分析方法提取的點云結果Fig.3 Point cloud extraction results using statistical analysis method

為解決上面提到的兩個問題，本文針對點云初步去噪算法處理過的點云，使用統計分析的方法對點云進行濾波，提取地面點，并進一步剔除地面以下噪聲點云。經過對點云的初步去噪，距離非噪聲點云較遠的點已被完全剔除，在提取地面點云時不會出現嚴重偏離實際值的情況。針對統計分析方法提取的植被冠層區域點云，考慮到地面以下距離地面較遠的點云大部分已被剔除，本文采用以下方法進行判定：若某段點云的模與該段點云最低點高程差值大于經驗閾值，認為該段提取點云為植被冠層區域點云，否則為地面點云。經驗閾值根據初步去噪結果和平均樹高進行選擇，針對本文實驗數據，經驗閾值設置為8m。對于被提取點云判定為地面點云的分段點云，選取高程值與模差值絕對值最小的點作為該段點云的地面種子點；對于被提取點云判定為植被冠層點云的分段點云，對該段點云兩側分段點云的地面種子點進行插值，將距離插值點最近的點作為該段點云的地面種子點。得到整個數據段的地面種子點后，即可利用最小二乘局部曲線擬合方法［12］進行點云濾波。因此，本文點云濾波過程如下：

步驟1：在沿軌方向以一定大小的窗口將經過去噪處理的點云分段，并將點云的高程值四舍五入至整數米。

步驟2：計算窗口內點云數據高程值的模值和最低點的高程，并計算模值與窗口內最低點高程的差值。若差值小于經驗閾值Δh， 則取窗口內與模值差值絕對值最小的點為該窗口內的地面種子點；若差值大于Δh， 則選取距離該窗口左右兩側窗口的地面種子點的連線垂直距離最近的點為該窗口內的地面種子點。

步驟3：對于每個窗口內的地面種子點，搜索在沿軌方向距離該點最近的9個種子點，以這10個點擬合二次曲線。

步驟4：設置自適應閾值［17］以區分地面點和非地面點

式（3）中，hmax為窗口內最高點高程；hmin為窗口內最低點高程；χ根據地形確定，一般設為0.1。計算窗口內每個點到步驟3得到的二次曲線在垂直方向的距離，若該距離小于或等于閾值，判定該點為地面點，否則判定為非地面點。

在對全部點云進行濾波處理后，點云可分為3部分，即：地面點云、植被冠層點云（地面點云之上的點云）和噪聲點云（地面點云之下的點云）。因此，該濾波算法進一步剔除了部分噪聲點云。獲取地面點云后，即可進一步提取DEM，并與該區域傳統激光雷達DEM數據進行對比。

2 數據處理與結果分析

2.1 點云初步去噪實驗

按照上文提出的點云初步去噪算法對MABEL點云實驗數據進行去噪處理，k設置為50，平均距離Di的頻數直方圖如圖4所示。圖4中對平均距離大于10m的頻數縮小了縱坐標范圍重新顯示，以便于查看。由圖4可以看出，點云平均距離主要集中在10m以內，且這個范圍內頻數隨點云平均距離增加而減少，直方圖峰值位置對應的平均距離為0.45m，最小平均距離為0.11m。因此，擬合Gauss函數均值為0.45m，標準差為0.34m。經驗值t選取的原則為：使去噪結果中有盡量多的噪聲點被剔除，而盡量多的非噪聲點被保留。根據本文實驗數據，經驗值t取16時去噪效果較好，此時閾值為5.89m，初步去噪結果如圖5所示。為了進行對比，本文使用文獻［12］中的去噪算法對該數據進行了處理，k仍設置為50，擬合Gauss函數均值為3.15m，方差為2.58m，經驗值t取5，去噪結果如圖6所示。圖5和圖6中，綠色點云為提取的非噪聲點。從圖中可以看出，兩種算法都較好地實現了噪聲點與非噪聲點的分離。還可以看出的是，圖5中更多地面下噪聲點被探測出來。實際上該數據共有點數51986個，圖5中被分類為噪聲的點數為1792個，圖6中被分類為噪聲的點數為1603個。

圖4 各點平均距離頻數直方圖Fig.4 Frequency histogram of each point average distance

圖5 本文算法初步去噪結果Fig.5 Preliminary de-noising results using the algorithm in this paper

圖6 文獻［12］算法去噪結果Fig.6 De-noising results using the algorithm proposed by reference ［12］

2.2 點云濾波處理實驗

將點云初步去噪結果以20m大小的窗口進行分段，將每段窗口內點云高程值四舍五入后計算窗口內點云的模值，將窗口內與模值高程差絕對值最小的點作為候選地面種子點，圖7顯示了經2.1節初步去噪之后保留的點云結果與候選地面種子點分布情況。從候選地面點的沿軌方向的分布中依稀可以看出地面地形，但是其中有大量的植被點云被選取出來。如果直接利用這些點恢復地形，將造成很大誤差。由圖7還可以看出，真實地面一般都在高程最低點之上，若直接采用窗口內高程最低點作為地面種子點，同樣會造成較大誤差。為選取合適的地面種子點，本文計算模值與該窗口最低點高程的差值。當差值小于8m時，保留該點為地面種子點，否則剔除該點，并選擇距離該點兩側窗口內地面種子點連線垂直距離最近的點為該窗口內的地面種子點。經過處理后，得到最終的地面種子點，如圖8所示。由圖8可以看出，雖然有少量地面種子點偏離真實地面，但是沿軌方向地面種子點的分布已經比較接近真實地面地形。

圖7 初步去噪后保留的點云與候選地面種子點Fig.7 Remained point cloud and candidate ground seed points after preliminary de-noising

圖8 初步去噪后保留的點云與地面種子點Fig.8 Remained point cloud and ground seed points after preliminary de-noising

為了進一步提取地面點云，本文聯合與每個窗口地面種子點距離最近的9個地面種子點擬合二次曲線

式（4）中， （x， y）為地面種子點坐標，p1、p2、p3為擬合參數。本文選擇較多的種子點來擬合二次曲線，這能減小偏離真實地面的種子點對真實地形的影響。對于窗口內每個點，計算x取該點橫坐標時二次曲線函數值與該點高程的差值Δh， 將差值絕對值與由式（3）得到的自適應閾值T進行比較，當≤T時，判定該點為地面點；當Δh＞0、＞T時，則判定該點為植被點；當Δh＜0、＞T時，則判定該點為噪聲點。照此方法對所有窗口內點云進行分類，結果如圖9所示。濾波后，點云被分為植被、地面、噪聲3類。其中，植被點為18988個，地面點為23421個，噪聲點為7785個。從濾波結果可以看出，地面點、植被點和噪聲點得到了較好的區分。

圖9 濾波后的結果Fig.9 Diagram of filtered results

2.3 實驗結果分析

（1）點云去噪結果分析

為了分析本文所采用的兩種去噪算法之間的區別，分別從圖5和圖6中截取兩段點云數據放大顯示， 如圖10所示。 其中， 圖10（a）和圖10（c）分別為截取自圖5沿軌方向上600m～1800m、6000m～7200m的點云， 圖10（b）和圖10（d）分別為截取自圖6沿軌方向上600m～1800m、6000m～7200m的點云。對比圖10（a）和圖 10（b）可以看出， 本文去噪算法和文獻［12］去噪算法在植被以上區域剔除的噪聲光子相同；而在地面點云和臨近地面以下的點云，二者處理結果有較大差別。這些差別主要表現在兩方面：1）文獻［12］去噪算法中有一部分地面點云被誤分類為噪聲點，如沿軌方向680m～780m的裸地區域以及780m～1100m的被植被覆蓋地面。這兩類點的共同特征是分布比較稀疏，因此較難與非噪聲點云進行區分。而本文去噪算法較好地將這些地面點正確提取出來，只有在740m處兩個地面點被誤判為噪聲點。2）在臨近地表以下區域的噪聲點云中，本文算法識別了更多的噪聲點。在圖10（a）中，去噪后地表以下信號點距離地面最大距離約15m，而圖10（b）中這個距離在8m以內。這非常有利于在下一步的濾波處理中設置較小閾值，以區分地面點云和植被點云。 圖10（c）和圖 10（d）除了在地面點云和臨近地面以下點云處理結果上與圖10（a）和圖10（b）表現出一致的差別外，二者在植被以上的噪聲光子上的表現也不盡相同。其中，圖10（a）（本文算法）成功剔除了6622m和7175m處兩個孤立的噪聲點，而圖10（b）結果中將其判斷為非噪聲點。

（a）截取自圖5沿軌方向上600m～1800m的點云

（b）截取自圖6沿軌方向上600m～1800m的點云

總體來講，在整個數據處理結果中，文獻［12］去噪算法將部分非噪聲點云分類為噪聲點，同時將部分噪聲點分類為地面點，尤其是臨近地面以下的噪聲點，而這兩種錯誤分類的情況在本文算法結果中都相對較少。例如，經人工判讀方法簡單統計，文獻［12］算法的去噪結果將地面點分類為噪聲點的個數為75個，而本文算法僅將6個地面點錯誤分類為噪聲點。由于被錯誤分類為噪聲點的地面點多處于點云稀疏的區域，這些點被剔除后相當于進一步抽稀了此處點云密度，甚至造成點云漏洞，這將會給下一步的濾波帶來較大的誤差。而臨近地面以下噪聲點的去除則對下一步濾波過程中的閾值設置非常重要，當最低點與實際地面高程差越小時，閾值設置越能有效將植被與地面分離。當高程差接近或達到樹冠高度時，則很難找到合適的閾值，甚至導致濾波算法失效。因此，經過本文初步去噪處理，噪聲點云和非噪聲點云得到了較好的分離，同時為進一步的濾波提供了較好的實驗條件。

圖10 截取的兩段點云數據放大圖Fig.10 Enlarged images of the two segmented point cloud data

（2）點云濾波結果分析

本文采用定量的方法分析濾波算法提取的地面點云數據質量，將實驗結果與該區域采用傳統激光雷達得到的DTM數據進行對比。由于所采用的DTM數據高程采用EGM96大地高，而MABEL點云高程為WGS84橢球高，因此需要對DTM進行數據轉換，這可以利用美國大地測量局提供的高程轉換工具［18］完成。對于濾波后地面點云中的每一點，本文利用雙線性內插法計算對應該點平面坐標處的DTM高程，最終得到MABEL點云和對應DTM點高程在沿軌方向上的分布，如圖11（a）所示。計算可得到MABEL地面點云相對于傳統激光雷達高程的均方根誤差RMSE為2.98m，相關系數R2為0.9938。如果對提取MABEL地面點云進一步規則化處理，可以得到線狀DEM，將DEM分辨率設置為3m，并將三維點云二次曲面內插法簡化為曲線形式，即采用二次曲線內插方法內插得到待求點高程，同樣利用雙線性內插法計算對應該點平面坐標處DTM高程，得到MABEL點云生成的DEM和對應DTM點高程在沿軌方向上的分布，如圖12（a）所示。可計算得到MABEL點云生成的DEM相對于傳統激光雷達高程的均方根誤差RMSE為2.85m，相關系數R2為0.9931。MABEL點云生成的DEM與傳統激光雷達相關性較MABEL原始地面點云更好，均方根誤差也較小，這是由于在內插獲得DEM時減弱了誤差較大點的影響。在圖 11（a）和圖 12（a）中， 沿軌方向 8000m～10000m的部分區域沒有數據，這是由于部分DTM數據缺失造成的（如圖2所示）。在這兩幅圖中，也可以看出MABEL數據與DTM數據之間的高相關性。本文使用初步去噪后點云高程的最低點作為地面種子點進行濾波實驗，得到MABEL點云、MABEL點云生成的DEM和對應DTM點高程在沿軌方向上分布，分別如圖11（b）、 圖 12（b）所示。 計算得到 MABEL地面點云、MABEL點云生成的DEM相對于傳統激光雷達高程的均方根誤差RMSE分別為6.50m和8.02m，相關系數R2分別為0.9728和0.9702。從圖11（b）、 圖12（b）中明顯可以看出， 該方法提取的地面點云高程于傳統激光雷達高程。而在本文濾波實驗中，通過采用模值和窗口最低點高程的差值作為閾值，合理的確定了地面種子點，避免了選擇地面以下噪聲點作為地面種子點帶來的較大誤差。因此，利用本文統計分析濾波提取得到的DEM和DTM與傳統采用高程最低點作為種子點得到的DEM和DTM相比更接近于參考數據。

圖11 MABEL地面點云與DTM內插數據對比Fig.11 Comparison between MABEL ground point cloud and DTM interpolation data

圖12 MABEL地面點云生成DEM與DTM內插數據對比Fig.12 Comparison between DEM generated by MABEL ground point cloud and DTM interpolation data

3 結論

本文提出了改進的基于局部距離統計的點云去噪算法，利用完整的MABEL數據段進行實驗，得到了較好的去噪效果。在MABEL數據去噪實驗中，傳統算法將部分非噪聲點云分類為噪聲點，同時將部分噪聲點分類地面點，尤其是臨近地面以下的噪聲點，而這兩種錯誤分類的情況在本文算法結果中都相對較少。傳統算法去噪結果將地面點分類為噪聲點的個數為75個，而本文算法僅將6個地面點錯誤分類為噪聲點。點云初步去噪后，仍然存在較多的噪聲信號，這可以通過濾波結果看出，臨近地面以下噪聲點占濾波前點云總數的比例達到15.5%。這些點也是傳統的去噪方法難以剔除的，因此本文從初步去噪后的點云中利用統計分析方法提取地面點，間接的剔除地面以下噪聲點，同時完成了點云分類。在濾波實驗過程中，采用了一些經驗值，這些數值的設置需要結合實驗數據特點進行合理選擇。最后，本文還利用傳統激光雷達獲取DTM驗證了MABEL地面點云提取結果。實驗結果顯示，利用本文算法得到的MABEL地面點云相對于傳統激光雷達高程的均方根誤差RMSE為2.98m，相關系數R2為0.9938。MABEL點云生成的DEM相對于傳統激光雷達高程的均方根誤差RMSE為2.85m，相關系數R2為0.9931。使用初步去噪后點云高程的最低點作為地面種子點得到MABEL地面點云、MABEL點云生成的DEM相對于傳統激光雷達高程的均方根誤差RMSE分別為6.50m和8.02m，相關系數R2分別為0.9728和0.9702。對兩種方法實驗結果對比可以看出，利用本文統計分析濾波提取得到的DEM和DTM的精度優于傳統采用高程最低點作為種子點得到的DEM和DTM。

本文主要關注MABEL數據中地面點云的提取與處理，雖然同時獲得了植被區域的點云，但并未對其質量進行評價。這主要是由于植被附近區域點云分布狀況復雜，一些與植被接近的噪聲點由于點云分布較為分散，人工尚且無法準確判讀。因此，植被區域點云噪聲去除還需要進行更深入的研究。