單光子激光雷達數(shù)據(jù)去噪與濾波算法

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

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

0 引言

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

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

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

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

1.1 實驗數(shù)據(jù)介紹

（1）MABEL 數(shù)據(jù)

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

圖1 MABEL剖面點云數(shù)據(jù)Fig.1 Diagram of MABEL point cloud data

（2）傳統(tǒng)機載激光雷達數(shù)據(jù)

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

圖2 DTM數(shù)據(jù)與MABEL數(shù)據(jù)的疊加顯示Fig.2 Overlay diagram of DTM data and MABEL data

1.2 MABEL點云去噪和濾波算法

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

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

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

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

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

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

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

步驟3：以0.1m為距離間隔統(tǒng)計平均距離Di的頻數(shù)，生成頻數(shù)直方圖。

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

（2）基于統(tǒng)計分析的點云濾波算法

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

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

圖3 使用統(tǒng)計分析方法提取的點云結(jié)果Fig.3 Point cloud extraction results using statistical analysis method

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

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

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

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

步驟4：設(shè)置自適應(yīng)閾值［17］以區(qū)分地面點和非地面點

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

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

2 數(shù)據(jù)處理與結(jié)果分析

2.1 點云初步去噪實驗

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

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

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

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

2.2 點云濾波處理實驗

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

圖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

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

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

圖9 濾波后的結(jié)果Fig.9 Diagram of filtered results

2.3 實驗結(jié)果分析

（1）點云去噪結(jié)果分析

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

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

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

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

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

（2）點云濾波結(jié)果分析

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

圖11 MABEL地面點云與DTM內(nèi)插數(shù)據(jù)對比Fig.11 Comparison between MABEL ground point cloud and DTM interpolation data

圖12 MABEL地面點云生成DEM與DTM內(nèi)插數(shù)據(jù)對比Fig.12 Comparison between DEM generated by MABEL ground point cloud and DTM interpolation data

3 結(jié)論

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

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