密度聚類算法在光子點(diǎn)云去噪中的應(yīng)用與評(píng)價(jià)

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

(1.地理信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，西安 710054；2.西安測(cè)繪研究所，西安 710054)

0 引言

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

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

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

1 DBSCAN去噪方法原理

1.1 DBSCAN算法

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

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

圖1 ICESat-2典型光子計(jì)數(shù)點(diǎn)云剖面圖

1.2 參數(shù)確定方法

1)鄰域最小點(diǎn)數(shù)MinPts。經(jīng)驗(yàn)公式一：文獻(xiàn)[12]通過(guò)定義目標(biāo)點(diǎn)范圍內(nèi)密度值來(lái)確定算法MinPts。

(1)

MinPts≥4·ρ

(2)

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

經(jīng)驗(yàn)公式二：文獻(xiàn)[13]在處理海洋場(chǎng)景淺海測(cè)深點(diǎn)云去噪時(shí)，給出以下MinPts經(jīng)驗(yàn)公式。

(3)

(4)

(5)

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

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

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

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

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

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

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

圖3 典型場(chǎng)景下聚類數(shù)與K的關(guān)系

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

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

同時(shí)，本文通過(guò)試錯(cuò)觀察等方法表明：最佳K值如圖3中紅點(diǎn)所示位置，第一種情況位于兩波峰之間，第二種情況位于波峰右側(cè)。因此本文改進(jìn)Eps自動(dòng)確定的步驟如下。

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

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

2 實(shí)驗(yàn)和分析

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

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

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

表1 實(shí)驗(yàn)數(shù)據(jù)詳細(xì)信息

2.2 實(shí)驗(yàn)過(guò)程及結(jié)果分析

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

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

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

從表2可知，Data1、Data2兩處植被區(qū)的MinPts和Eps區(qū)域?qū)Ρ榷圆町惒淮?，分別為11、11.5和13、12.5，固定參數(shù)時(shí)DB1計(jì)算時(shí)間僅僅為42″和55″，而自適應(yīng)算法DB2時(shí)間達(dá)到35′54″和68′32″，效率很低。

Data3、Data4兩處城市區(qū)的核心參數(shù)區(qū)域?qū)Ρ榷圆町愝^大，Data3的MinPts和Eps分別為12、7.3，Data4為6、14.9，這是由于Data3的點(diǎn)密度更大，并且美國(guó)該城市區(qū)域?yàn)榈桶ㄖ珼ata4西安城市點(diǎn)密度稍小，且建筑高度較大，因此需要更大的Eps才能保證將地形和建筑物這種在剖面上斷開(kāi)的對(duì)象聚類為一類。

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

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

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

2.3 去噪精度

本文自適應(yīng)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去噪效果，從目視觀察來(lái)看，本文DBSCAN算法均有效識(shí)別出信號(hào)，Data1的植被和地形線被正確識(shí)別，Data3城市區(qū)域沒(méi)有因?yàn)榈匦尾贿B續(xù)而出現(xiàn)識(shí)別錯(cuò)誤，Data5的冰面點(diǎn)云密集度高，也被正確識(shí)別為信號(hào)點(diǎn)。不足之處在于地形線下方有明顯的噪聲點(diǎn)被識(shí)別為信號(hào)，這是空間密度類算法的固有缺陷。以Data5為例，最佳Eps=2.9 m時(shí)，地形線下方2.9 m以內(nèi)的噪聲點(diǎn)會(huì)被識(shí)別為信號(hào)，而冰面處地形平坦，點(diǎn)密度高，地形線下方的噪聲點(diǎn)在視覺(jué)上很容易發(fā)現(xiàn)，但自動(dòng)算法卻分類錯(cuò)誤。

表3 實(shí)驗(yàn)數(shù)據(jù)自適應(yīng)DBSCAN去噪混淆矩陣

圖4 Data1植被DBSCAN去噪效果

圖5 Data3城市DBSCAN去噪效果

圖6 Data5冰蓋DBSCAN去噪效果

3 結(jié)束語(yǔ)

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