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

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

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

0 引言

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

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

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

1 單光子激光雷達(dá)點(diǎn)云去噪和濾波模型

1.1 實(shí)驗(yàn)數(shù)據(jù)介紹

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

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

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

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

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

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

1.2 MABEL點(diǎn)云去噪和濾波算法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.1 點(diǎn)云初步去噪實(shí)驗(yàn)

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

圖4 各點(diǎn)平均距離頻數(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 文獻(xiàn)［12］算法去噪結(jié)果Fig.6 De-noising results using the algorithm proposed by reference ［12］

2.2 點(diǎn)云濾波處理實(shí)驗(yàn)

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

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

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

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

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

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

2.3 實(shí)驗(yàn)結(jié)果分析

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

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

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

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

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

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

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

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

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

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

3 結(jié)論

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

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