曹彬才，王建榮，胡燕，呂源，楊秀策

(1.地理信息工程國家重點實驗室，西安 710054；2.西安測繪研究所，西安 710054)

0 引言

美國2018年9月發(fā)射了第二代對地觀測激光雷達衛(wèi)星ICESat-2(ice，cloud and land elevation satellite-2)，旨在利用激光雷達高精度特點在全球范圍內(nèi)開展冰蓋變化監(jiān)測、海面高度測量、植被覆蓋反演等，為研究碳循環(huán)和全球變暖等科學問題提供技術支撐[1]。與此同時，ICESat-2高精度點云數(shù)據(jù)在高精度地形測量、激光點云輔助遙感影像平差、淺海水深測繪等[2-4]傳統(tǒng)測繪鄰域也得到了廣泛應用，ICESat-2點云數(shù)據(jù)處理成為測繪行業(yè)的研究熱點。

ICESat-2搭載的激光載荷ATLAS(advanced topographic laser altimeter system)采用了光子計數(shù)探測體制，這種新型激光雷達具有功耗小、重量輕、靈敏度高等特點[5]，比傳統(tǒng)線性探測激光雷達更適合用作星載平臺。高靈敏度探測帶來的缺點是數(shù)據(jù)噪聲多，ATLAS接收器會將大量的太陽輻射光子和大氣散射光子記錄為光子事件，導致數(shù)據(jù)信噪比差，因此點云去噪對光子計數(shù)激光雷達至關重要。

當前光子點云去噪總體上有兩種思路[6]：一是基于數(shù)字圖像處理技術，將剖面點云轉換為影像，采用邊緣檢測[7]、區(qū)域檢測等識別噪聲；二是逐點計算某個局部統(tǒng)計量，利用其分布特征計算全局閾值并區(qū)分信號和噪聲，如直方圖法[8]、點云聚類法[9-10]、隨機森林法[11]等。其中DBSCAN(density-based spatial clustering of applications with noise)是一種較為成熟的聚類算法，能將半徑Eps范圍內(nèi)點個數(shù)超過最小值MinPts的點歸為一類，廣泛用于激光點云去噪。Zhang等[12]將密度均值作為MinPts，并將圓形搜索核改進為平行地面的橢圓，處理有植被區(qū)域的MABEL(multiple altimeter beam experimental LiDAR)單光子數(shù)據(jù)時效果較好；Ma等[13]在處理淺海測深數(shù)據(jù)時使用DBSCAN進行去噪，給出了海洋場景下MinPts經(jīng)驗公式，半徑Eps則采用經(jīng)驗值，在淺海區(qū)域也能夠有效識別噪聲；李文杰等[14]通過數(shù)據(jù)集本身分布特征生成Eps候選集，通過窮舉方式確定最佳參數(shù)；魏碩等[15]通過k維樹求取點云密度進行粗去噪，然后運用改進DBSCAN算法和統(tǒng)計濾波算法進行精去噪。本文詳細論述了當前光子計數(shù)激光雷達去噪中DBSCAN具體應用方法，探討了關鍵參數(shù)MinPts、Eps自適應確定的可行性，并使用不同地物類型ICESat-2數(shù)據(jù)開展精度驗證。

1 DBSCAN去噪方法原理

1.1 DBSCAN算法

DBSCAN算法是一種基于密度的空間聚類算法，通過尋找密度相連的點的最大集合來分離信號點和噪聲點[16]。對于集群的每個點，給定半徑Eps鄰域必須至少包含最少數(shù)量的點，即鄰域中的密度必須超過閾值MinPts。DBSCAN能自動將密度足夠大的點區(qū)域劃分為簇，數(shù)據(jù)集中不屬于任何簇的點則被視為噪聲。該算法原理簡單、運算速度快，且不需要預先指定簇的個數(shù)，核心問題是要輸入半徑Eps和鄰域最小點數(shù)MinPts兩個基本參數(shù)。

如圖1所示，ICESat-2光子點云常被投影到沿飛行方向和高程方向組成的二維平面，X軸表示沿飛行方向距離，可通過數(shù)據(jù)文件中的速度與時間參量相乘得到，Y軸表示高程。圖中可明顯觀察到位于地形線附近點云更加密集，地形線上下方噪聲點在空間分布上更加稀疏，這也是基于密度算法成功識別信號和噪聲的必要條件。

圖1 ICESat-2典型光子計數(shù)點云剖面圖

1.2 參數(shù)確定方法

1)鄰域最小點數(shù)MinPts。經(jīng)驗公式一：文獻[12]通過定義目標點范圍內(nèi)密度值來確定算法MinPts。

(1)

MinPts≥4·ρ

(2)

式中：ρ表示密度值；N是數(shù)據(jù)集點總數(shù)；S為數(shù)據(jù)集二維剖面面積；s1表示局部搜索范圍面積，s1=π·Eps2。文獻[10]在處理MABEL數(shù)據(jù)時將Eps設為固定值2，并根據(jù)式(1)、式(2)計算MinPts。

經(jīng)驗公式二：文獻[13]在處理海洋場景淺海測深點云去噪時，給出以下MinPts經(jīng)驗公式。

(3)

(4)

(5)

式中：S′、N′分別是某背景范圍面積和對應的點總數(shù)。選擇高程Y軸最低的5 m范圍作為背景，沿軌方向X保持不變，仍為S對應的沿軌長度。

2)半徑Eps?，F(xiàn)有去噪算法處理ICESat-2數(shù)據(jù)時搜索半徑基本都采用經(jīng)驗參數(shù)[17-18]，文獻[14]提出了一種基于K平均最鄰近法尋找最優(yōu)Eps參數(shù)的思路，步驟如下。

步驟1：根據(jù)K平均最鄰近法求出數(shù)據(jù)集D的候選Eps參數(shù)集合DEps。首先計算每個點的K最鄰近點對應的距離(K=1，2，3，…，N)，形成距離矩陣DN×N，該矩陣中每一行表示一個點，每一列表示從小到大排列的K最鄰近距離，每一列平均值即為K平均最鄰近距離。

步驟2：在已知MinPts下，將DEps中Eps參數(shù)逐個帶入DBSCAN算法中進行聚類運算，當生成的簇數(shù)連續(xù)3次相同時認為結果趨于穩(wěn)定，簇數(shù)M為最優(yōu)簇數(shù)。

步驟3：繼續(xù)執(zhí)行步驟2直到簇數(shù)不等于M，選用簇數(shù)為M時的最大K平均最鄰近距離作為最優(yōu)Eps參數(shù)。

3)本文改進方法。圖2展示了圖1數(shù)據(jù)Eps列表與K值的關系，隨著K值增大，Eps參數(shù)隨之增大。按照上述Eps尋優(yōu)計算思路，本文通過對多類型ICESat-2數(shù)據(jù)的實驗與分析，總結出DBSCAN聚類個數(shù)與K值關系(圖3)具有下述特點。

圖2 參數(shù)Eps列表與K的關系示意圖

圖3 典型場景下聚類數(shù)與K的關系

①植被、城市等陸地場景下，聚類個數(shù)與K值關系線呈現(xiàn)“雙峰”分布式，即K值從小到大，聚類數(shù)先變大，再減小并趨于平穩(wěn)，隨后再次變大，最后減小直到穩(wěn)定不發(fā)生變化。

②冰蓋海洋等場景下，聚類個數(shù)與K值關系線呈現(xiàn)“單峰”分布式，即聚類數(shù)隨著K值的增大而變大，隨即減小趨于平穩(wěn)。

同時，本文通過試錯觀察等方法表明：最佳K值如圖3中紅點所示位置，第一種情況位于兩波峰之間，第二種情況位于波峰右側。因此本文改進Eps自動確定的步驟如下。

步驟1：利用MinPts值限制最大K值范圍，Kmax=D·MinPts，Kmax表示最大K值，D=15，此范圍內(nèi)足以使得聚類數(shù)趨于穩(wěn)定，同時避免計算所有K值，大幅減小計算量。

步驟2：生成聚類數(shù)與K值關系曲線，利用一階求導等方法判斷波峰個數(shù)。當曲線呈現(xiàn)雙峰分布時，取波峰之間聚類數(shù)穩(wěn)定的最小K值為最佳K，對應Eps為最佳搜索半徑；當曲線呈單峰分布時，取波峰右側聚類數(shù)穩(wěn)定的最小K值為最佳K。

2 實驗和分析

2.1 ICESat-2單光子數(shù)據(jù)

ICESat-2衛(wèi)星軌道高度約500 km，受傳播過程中的大氣散射、目標漫反射等影響，一束激光脈沖僅有數(shù)個或數(shù)十個光子能返回接收器。不同拍攝條件、地物類型對應點云的信噪比和點云疏密程度有較大差異，ICESat-2在官方文檔泊松去噪算法里將地物類型劃分為陸地冰、海冰、海洋、陸地、內(nèi)陸水5種類型，分別設計經(jīng)驗參數(shù)[19]。

本文選擇了植被、城市、冰蓋、海洋4種不同類別的光子數(shù)據(jù)開展實驗，數(shù)據(jù)詳細信息見表1。Data1至Data8的拍攝地點分別為美國新罕布什爾州、陜西城固市、拉斯維加斯、西安市、格陵蘭島、南極大陸、南海珊瑚礁、夏威夷海岸，均選擇了強波束開展實驗。

表1 實驗數(shù)據(jù)詳細信息

2.2 實驗過程及結果分析

實驗中利用經(jīng)驗公式二計算MinPts，由于文獻[14]提出的計算全部K均值思路計算量過大，非常耗時，因此使用改進后的Eps確定方法。為便于分析，這里對比固定參數(shù)DBSCAN(記為DB1)、本文自適應DBSCAN(記為DB2)算法以及官方推薦算法DRAGANN，實驗主要對算法參數(shù)的自適應確定、計算耗時以及去噪精度等方面開展評價。

1)參數(shù)自適應確定。表2為算法參數(shù)自適應確定及耗時情況，表中MinPts和Eps為自適應算法DB2計算的最優(yōu)參數(shù)。本文程序采用VS2010 C++語言編寫，并使用激光點云處理庫PCL用作構建KDtree并尋找最鄰近K值，自適應參數(shù)算法中采用CPU多核并行運算。

表2 參數(shù)確定及耗時

從表2可知，Data1、Data2兩處植被區(qū)的MinPts和Eps區(qū)域對比而言差異不大，分別為11、11.5和13、12.5，固定參數(shù)時DB1計算時間僅僅為42″和55″，而自適應算法DB2時間達到35′54″和68′32″，效率很低。

Data3、Data4兩處城市區(qū)的核心參數(shù)區(qū)域對比而言差異較大，Data3的MinPts和Eps分別為12、7.3，Data4為6、14.9，這是由于Data3的點密度更大，并且美國該城市區(qū)域為低矮建筑，Data4西安城市點密度稍小，且建筑高度較大，因此需要更大的Eps才能保證將地形和建筑物這種在剖面上斷開的對象聚類為一類。

Data5、Data6是冰蓋高反射區(qū)域，平均一束激光返回的光子數(shù)要遠大于植被，因此MinPts比前幾組數(shù)據(jù)都要大，分別為58和30。Data5位于格陵蘭島平地區(qū)域，點密度較大，在搜索半徑Eps為1.97 m時即可有效區(qū)分信號和噪聲。Data5位于南極某山地，冰雪覆蓋略少，且存在起伏斷裂，因此自適應Eps稍大，為11.2 m。

Data7、Data8是海面區(qū)域，反射光子低于冰蓋區(qū)，與陸地相差不大，兩組數(shù)據(jù)的MinPts和Eps分別為14、7.4和6、6.1。

2)計算耗時。光子點云投影到飛行方向剖面后，實際上成為二維數(shù)據(jù)，DBSCAN算法的空間復雜度為O(n2)，n表示點個數(shù)，循環(huán)生成Eps列表空間復雜度為O(n)，因此總的空間復雜度為O(n2)+O(n)。

2.3 去噪精度

本文自適應DBSCAN算法去噪精度如表3所示，8組數(shù)據(jù)的整體去噪精度分別為97.3%、97.6%、99.2%、98.3%、97.2%、98.6%、97.9%和99.1%，均表現(xiàn)較為優(yōu)異，整體水平與泊松去噪、DRAGANN去噪等幾乎一致。圖4、圖5、圖6分別為Data1、Data3、Data5去噪效果，從目視觀察來看，本文DBSCAN算法均有效識別出信號，Data1的植被和地形線被正確識別，Data3城市區(qū)域沒有因為地形不連續(xù)而出現(xiàn)識別錯誤，Data5的冰面點云密集度高，也被正確識別為信號點。不足之處在于地形線下方有明顯的噪聲點被識別為信號，這是空間密度類算法的固有缺陷。以Data5為例，最佳Eps=2.9 m時，地形線下方2.9 m以內(nèi)的噪聲點會被識別為信號，而冰面處地形平坦，點密度高，地形線下方的噪聲點在視覺上很容易發(fā)現(xiàn)，但自動算法卻分類錯誤。

表3 實驗數(shù)據(jù)自適應DBSCAN去噪混淆矩陣

圖4 Data1植被DBSCAN去噪效果

圖5 Data3城市DBSCAN去噪效果

圖6 Data5冰蓋DBSCAN去噪效果

3 結束語

DBSCAN作為典型的基于空間密度的聚類算法，具有實現(xiàn)簡單、分類效果好、能實現(xiàn)任意形狀的點聚類等特點。缺點是DBSCAN算法對核心參數(shù)MinPts和Eps非常敏感，兩個參數(shù)直接決定聚類效果。本文研究了多場景下ICESat-2光子點云DBSCAN去噪效果。實驗表明，不同場景下的光子點云特征差異較大，不宜使用一個固定參數(shù)處理多種場景。針對不同目標類型，以最終聚類數(shù)為參考雖然可以自適應確定DBSCAN算法的Eps參數(shù)，實現(xiàn)全自動化去噪，但這種思路缺乏光子對光子密度的針對性設計，計算復雜度高，耗時較長，與現(xiàn)有的主流算法相比效率較低。下一步應針對不同地物類型在不同高度的點云密度變化特點，設計自適應濾波模型，力求在處理效率上得到大幅改進。