孟文君，李杰，張凱，2，劉長(zhǎng)達(dá)，唐秋華，

（1.山東科技大學(xué) 測(cè)繪與空間信息學(xué)院，山東 青島 266590；2.自然資源部海洋測(cè)繪重點(diǎn)實(shí)驗(yàn)室，山東 青島 266590；3.自然資源部第一海洋研究所，山東 青島 266061）

第二代星載激光雷達(dá)冰、云和陸地測(cè)高衛(wèi)星ICESat-2（The Ice，Cloud，and Land Elevation Satellite-2）是美國(guó)國(guó)家航空航天局（National Aero原nautics and Space Administration，NASA）繼ICESat失效后，于2018 年發(fā)射的又一星載激光雷達(dá)衛(wèi)星。不同于ICESat 搭載的全波形激光雷達(dá)系統(tǒng)（Geo原science Laser Altimeter System，GLAS），ICESat-2搭載了光子計(jì)數(shù)激光測(cè)高儀（Advanced Topograph原ic Laser Altimeter System，ATLAS），使用全新的微脈沖多波束光子計(jì)數(shù)式激光雷達(dá)，這是該技術(shù)首次應(yīng)用于衛(wèi)星平臺(tái)（Markus et al，2017；Neumann et al，2019）。光子計(jì)數(shù)激光雷達(dá)系統(tǒng)采用更加靈敏的單光子探測(cè)器，具有更高的脈沖重復(fù)頻率，更高的精度和更小的激光足?。∟eumann et al，2019），主要用于冰蓋高程測(cè)量及變化監(jiān)測(cè)、海冰高程測(cè)量及厚度反演、陸地高程測(cè)量和森林探測(cè)、湖泊水位和海平面變化監(jiān)測(cè)等（葛莉等，2017；Wang et al，2011；楊帆 等，2011；胡國(guó)軍 等，2015；胥喆等，2017）。但同時(shí)由于發(fā)射的激光脈沖能量小，記錄的光子事件中包括各類噪聲（大氣散射、太陽(yáng)輻射及儀器噪聲），信噪比低，給光子數(shù)據(jù)處理與應(yīng)用研究帶來(lái)了巨大挑戰(zhàn)。

針對(duì)光子去噪問(wèn)題，目前主要有以下三類方法：（1）基于柵格圖像處理的去噪算法（Chen et al，2015； Magruder et al，2012），該算法原理簡(jiǎn)單，但在點(diǎn)云和柵格圖像的轉(zhuǎn)換中會(huì)損失部分精度；（2）基于局部統(tǒng)計(jì)參數(shù)的去噪算法（Gwenzi et al，2016；謝峰 等，2017；Herzfeld et al，2019；夏少波等，2014），該算法有較好的適應(yīng)性，但受閾值影響較大，如何選取有效閾值仍有待研究；（3）基于密度空間聚類的去噪算法（Zhang et al，2014；Zhu et al，2020；Wang et al，2016），該算法通過(guò)利用噪聲光子在空間中分布離散的特點(diǎn)進(jìn)行聚類分析，從而剔除噪聲光子。主要的密度聚類方法有貝葉斯（李明等，2018）、DBSCAN（Density-Based Spatial Clustering of Applications with Noise）和OPTICS（Ordering Points to Identify the Clustering Structure）。在這三種方法中，DBSCAN 算法原理較為簡(jiǎn)單，也是目前針對(duì)光子去噪最常用的手段之一。

由于DBSCAN 算法對(duì)輸入?yún)?shù)較為敏感，如何選擇合適的參數(shù)成為影響光子去噪效果的關(guān)鍵因素。秦磊等（2020）通過(guò)計(jì)算Eps 取0～30 內(nèi)的整數(shù)時(shí)的F1原score，得到最佳Eps 值，從而進(jìn)行DBSCAN 聚類，但這種方法具有較大的偶然性；馬躍等（2020） 提出了一種計(jì)算淺海中MinPts（Minimum points） 參數(shù)的改進(jìn)公式，聚類效果良好，但并未說(shuō)明Eps 的計(jì)算方法。

針對(duì)上述問(wèn)題，本文通過(guò)改進(jìn)運(yùn)-平均最近鄰法確定最優(yōu)Eps 參數(shù)，結(jié)合馬躍提出的計(jì)算公式，得到最佳聚類參數(shù)進(jìn)行DBSCAN 聚類，從而區(qū)分出噪聲光子和信號(hào)光子，并對(duì)結(jié)果進(jìn)行了精度驗(yàn)證。

1 研究區(qū)域與數(shù)據(jù)

ICESat-2 是NASA 于2018 年9 月15 日發(fā)射的一顆非太陽(yáng)同步軌道衛(wèi)星，軌道高度約500 km，以91 天為一個(gè)重訪周期，可覆蓋全球88毅N—88毅S的區(qū)域范圍，將全球分為1387 個(gè)地面參考航跡。其搭載的光子計(jì)數(shù)激光測(cè)高儀ATLAS 采用全新的光子計(jì)數(shù)模式，這種測(cè)量方式以10 Hz 的高重復(fù)頻率發(fā)射低脈沖能量（48～170 滋J）的單脈沖藍(lán)綠激光（532 nm），并且具有1.5 ns 的較短脈沖寬度，可沿軌連續(xù)獲取足印間距約0.7 m、足印大小17 m的重疊光斑。ATLAS 同時(shí)發(fā)射3 對(duì)激光，相鄰每對(duì)激光的間距為3.3 km，每對(duì)激光包含一強(qiáng)一弱兩束激光，組間強(qiáng)弱激光的間距為90 m，其能量之比約4 頤1，以便實(shí)現(xiàn)對(duì)坡度地區(qū)的檢測(cè)。

ICESat-2 共有Level-1、Level-2、Level-3A、Level-3B 四級(jí)，以及ATL00—ATL21 等21 種標(biāo)準(zhǔn)數(shù)據(jù)產(chǎn)品。其中Level-2 級(jí)產(chǎn)品ATL03 全球定位光子數(shù)據(jù)，記錄了光子數(shù)據(jù)的緯度、經(jīng)度、高度、光子往返時(shí)間、激光器位置和姿態(tài)角等信息，并提供陸地、海洋、海冰、陸地冰、內(nèi)陸水5 種地表類型的“置信度”參數(shù)來(lái)分離噪聲點(diǎn)和信號(hào)點(diǎn)。本文使用ATL03 數(shù)據(jù)進(jìn)行去噪處理。

為了驗(yàn)證去噪算法在淺海區(qū)域的有效性，本文選取了ICESat-2 衛(wèi)星于2019 年3 月23 日和2019年7 月21 日（Neumann et al，2019）經(jīng)過(guò)南海某兩處島礁的測(cè)量數(shù)據(jù)進(jìn)行實(shí)驗(yàn)（圖2）。

圖1 ICESat-2 波束分布示意圖

兩處研究區(qū)域均位于中國(guó)南海，區(qū)域A 所在的礁盤東西長(zhǎng)約6.8 km，南北寬約3.5 km，面積約0.4 km2，臺(tái)風(fēng)大潮時(shí)礁盤常被海水淹沒(méi)。區(qū)域B所在的礁盤東西長(zhǎng)約29.632 km，南北寬約9.26 km，低潮時(shí)礁盤可露出海面。兩區(qū)域所在的礁盤相距約12.96 km。

為了驗(yàn)證本文去噪方法在淺海區(qū)域的有效性，利用藍(lán)綠激光對(duì)水的穿透性，分別截取ICESat-2軌跡中經(jīng)過(guò)兩處礁盤區(qū)域的一段作為實(shí)驗(yàn)樣本進(jìn)行驗(yàn)證（圖2）。這兩處區(qū)域包含水深約1.5～3 m 的較淺區(qū)域，藍(lán)綠激光可以探測(cè)到水底，對(duì)驗(yàn)證淺海區(qū)域的數(shù)據(jù)去噪效果具有一定的代表性。

圖2 ICESat-2 激光光束在南海某島礁區(qū)域附近的軌跡

2 研究方法

2.1 DBSCAN 算法原理

ATLAS 激光器發(fā)射和接收的信號(hào)均為弱信號(hào)，大量噪聲混雜在信號(hào)光子之間，但是ATL03 數(shù)據(jù)在空間分布上信號(hào)比噪聲更加密集，可以利用噪聲比信號(hào)更加離散的空間分布特點(diǎn)，使用DBSCAN算法剔除離群點(diǎn)。

DBSCAN 算法由MartinEster、Hans-PeterKriegel等人于1996 年提出，該算法通過(guò)尋找一個(gè)個(gè)密度相連的點(diǎn)的最大集合（即簇）來(lái)分離信號(hào)點(diǎn)和噪聲點(diǎn)?；贒BSCAN 算法，對(duì)于聚類中的所有點(diǎn)，當(dāng)其在給定半徑鄰域內(nèi)的點(diǎn)密度超過(guò)給定閾值時(shí)，該點(diǎn)將被識(shí)別為“信號(hào)”。鄰域由給定兩點(diǎn)的加權(quán)距離函數(shù)定義。DBSCAN 算法簡(jiǎn)單且運(yùn)算速度較快，算法不需要預(yù)先指定簇的個(gè)數(shù)，但需要輸入兩個(gè)重要的參數(shù)：Eps （鄰域半徑）和MinPts（鄰域最小點(diǎn)數(shù)）。

由于算法存在對(duì)輸入的參數(shù)極為敏感的特點(diǎn)，因此選擇合適的MinPts 和Eps 參數(shù)尤為重要。本文的參數(shù)確定過(guò)程如下。

2.1.1 確定參數(shù)MinPts首先將飛躍研究區(qū)域沿軌道方向的每連續(xù)10 000 個(gè)原始光子設(shè)為一組，分組計(jì)算MinPts 參數(shù)（Ma et al，2019，2020）。然后，按公式（1）計(jì)算給定半徑Ra內(nèi)包括噪聲光子和信號(hào)光子的預(yù)期總光子數(shù)雜暈員。

式中，N1是總光子數(shù)，h 代表每10 000 個(gè)點(diǎn)中最高點(diǎn)與最低點(diǎn)之差的垂直范圍，l 代表每10 000 個(gè)點(diǎn)中最遠(yuǎn)點(diǎn)與最近點(diǎn)之差的沿軌范圍。接下來(lái)，按公式（2） 計(jì)算給定半徑Ra內(nèi)的預(yù)期噪聲光子數(shù)雜暈圓。假設(shè)h2=5，相應(yīng)的N2是對(duì)應(yīng)于最低5 m 處的光子數(shù)。

最后，就可以通過(guò)公式（3）計(jì)算MinPts 值

2.1.2 生成候選Eps 列表 對(duì)于判定參數(shù)Eps，將使用K-平均最近鄰法和數(shù)學(xué)期望法生成候選Eps 列表（李文杰等，2019）。

為了生成候選Eps 列表，需要假設(shè)一個(gè)大小為1000 個(gè)點(diǎn)的窗口，按沿軌飛行方向移動(dòng)窗口，逐個(gè)計(jì)算每個(gè)窗口內(nèi)的Eps 值。在每個(gè)窗口內(nèi)計(jì)算所有點(diǎn)之間的歐氏距離，并按距離由近及遠(yuǎn)的順序進(jìn)行排序，生成一個(gè)距離矩陣。距離矩陣展示了窗口內(nèi)所有點(diǎn)之間排序后的距離情況，行號(hào)代表窗口內(nèi)的所有點(diǎn)，列號(hào)代表窗口內(nèi)的該點(diǎn)到其他點(diǎn)排序后的距離。為了獲取K-平均最近鄰距離，將距離矩陣按列取平均，此時(shí)行號(hào)為K 值，得到所有K 值對(duì)應(yīng)的平均距離，這些K 值和平均距離就構(gòu)成了候選Eps 列表。

以往的實(shí)驗(yàn)證明使用沿軌道方向每連續(xù)的10 000 個(gè)點(diǎn)算得的MinPts 值在該區(qū)間內(nèi)穩(wěn)定。為了更好地適應(yīng)地形起伏，取數(shù)據(jù)中連續(xù)的10 000 個(gè)點(diǎn)為一組，分組計(jì)算參數(shù)MinPts，MinPts 參數(shù)固定后，用移動(dòng)窗口的方式，逐個(gè)計(jì)算每個(gè)窗口內(nèi)的候選Eps 列表。

2.2 自適應(yīng)確定最優(yōu)參數(shù)

自適應(yīng)確定最優(yōu)參數(shù)的步驟如下：

步驟一：將所有數(shù)據(jù)按飛躍研究區(qū)域的沿軌方向進(jìn)行分組，每10 000 個(gè)點(diǎn)一組，按公式（3）計(jì)算一次MinPts。

步驟二：同一組內(nèi)的MinPts 參數(shù)相同，移動(dòng)窗口，根據(jù)K-平均最近鄰法依次生成每個(gè)窗口內(nèi)的候選Eps 列表。

步驟三：將候選Eps 列表和由公式得到的MinPts參數(shù)，輸入DBSCAN 算法對(duì)數(shù)據(jù)集進(jìn)行聚類分析，分別得到不同K 值下的聚類簇?cái)?shù)。當(dāng)生成的簇?cái)?shù)連續(xù)6 次相同且聚類數(shù)為2 時(shí)，認(rèn)為聚類結(jié)果趨于穩(wěn)定，數(shù)據(jù)集內(nèi)的光子點(diǎn)云被分成噪聲光子和信號(hào)光子兩類。

步驟四： 繼續(xù)執(zhí)行步驟三，直到循環(huán)結(jié)束，并選用當(dāng)簇?cái)?shù)為2 時(shí)所對(duì)應(yīng)的最小K 值作為最優(yōu)K 值。最優(yōu)K 值對(duì)應(yīng)的平均距離則為最優(yōu)Eps 參數(shù)。

2.3 精度評(píng)價(jià)指標(biāo)

ICESat-2 數(shù)據(jù)受背景噪聲，尤其是太陽(yáng)背景噪聲的影響很大，噪聲離散的分布于整個(gè)研究區(qū)域。通過(guò)定義精確度Precision、召回率Recall 和F1 值評(píng)定指標(biāo)評(píng)價(jià)精度。設(shè)精確度Precision 為判斷正確的數(shù)據(jù)個(gè)數(shù)占識(shí)別出的數(shù)據(jù)個(gè)數(shù)的百分比；召回率Recall 為判斷正確的數(shù)據(jù)個(gè)數(shù)占實(shí)際的數(shù)據(jù)總數(shù)的百分比。F1 值是為了平衡召回率和精確度，假定二者同等重要的一個(gè)綜合評(píng)價(jià)指標(biāo)。

式中，TP 是被正確分類為信號(hào)點(diǎn)的光子數(shù)，F(xiàn)P 是實(shí)際為噪聲但被錯(cuò)分為信號(hào)的光子個(gè)數(shù)，F(xiàn)N 是實(shí)際為信號(hào)點(diǎn)卻被錯(cuò)誤的地分為噪聲的光子數(shù)。

3 結(jié)果與分析

下文將以區(qū)域A 為例，詳細(xì)說(shuō)明實(shí)驗(yàn)的具體步驟。區(qū)域A 參與實(shí)驗(yàn)的光子共有4000 個(gè)，原始數(shù)據(jù)如圖3 所示，噪聲光子雜亂無(wú)章，信號(hào)光子和噪聲光子難以分開，給后續(xù)數(shù)據(jù)處理帶來(lái)困難。

圖3 原始數(shù)據(jù)所有光子（區(qū)域A）

該實(shí)驗(yàn)區(qū)域的全部光子按公式（3）參與計(jì)算MinPts 值，結(jié)果為MinPts =18。按K-平均最近鄰法分別計(jì)算4000 個(gè)光子在每個(gè)窗口內(nèi)的Eps 值。得到4 個(gè)聚類簇?cái)?shù)與K 值的關(guān)系圖，如圖4 所示?？梢钥闯龃?cái)?shù)穩(wěn)定在2 時(shí)的最小K 值分別為50、68、56、51。圖5 展示了K 值與Eps 的對(duì)應(yīng)關(guān)系，最小K 值對(duì)應(yīng)的最佳Eps 分別為1.67、1.89、1.53、1.43。用實(shí)驗(yàn)獲得的參數(shù)分別對(duì)每個(gè)窗口內(nèi)的點(diǎn)聚類，聚類結(jié)果如圖6 所示。

圖4 K 值與分類個(gè)數(shù)（區(qū)域A）

圖5 K 值與Eps（區(qū)域A）

圖6 各組分類情況（區(qū)域A）

為了更加直觀地表現(xiàn)本文方法的去噪效果，將本實(shí)驗(yàn)結(jié)果與ATL03 算法提取的官方去噪結(jié)果進(jìn)行對(duì)比分析。選擇ATL03 數(shù)據(jù)中提供的置信度為4的高置信度光子作為信號(hào)光子，其余光子均為噪聲光子進(jìn)行顯示。區(qū)域A 的具體去噪結(jié)果如圖7、圖8所示。用同樣的方法對(duì)區(qū)域B 進(jìn)行處理，得到去噪結(jié)果對(duì)比如圖9、圖10 所示。

圖7 本文方法去噪結(jié)果（區(qū)域A）

圖8 ATL03 算法提取的官方去噪結(jié)果（區(qū)域A）

圖9 本文方法去噪結(jié)果（區(qū)域B）

圖10 ATL03 算法提取的官方去噪結(jié)果（區(qū)域B）

可以看出，在區(qū)域A，官方提供的置信度分離噪聲光子的方法只去除了少部分噪聲，對(duì)于大多數(shù)的離群噪聲點(diǎn)，官方提供的分類方法沒(méi)有將其區(qū)分出來(lái)。本文的去噪方法，可將大部分噪聲點(diǎn)去除，但將橫軸距離0～5 m 處的部分信號(hào)光子誤判為噪聲，這是由于DBSCAN 聚類方法自身的缺陷導(dǎo)致。針對(duì)區(qū)域B，使用官方提供的置信度分離噪聲光子，同樣未將大部分噪聲光子識(shí)別出來(lái)。本文的去噪方法，只有少量殘留噪聲點(diǎn)未被去除，這些未被去除的光子位置靠近海面，其密度與海面信號(hào)光子密度相近，這是由于DBSCAN 算法難以區(qū)分密度相近的點(diǎn)云數(shù)據(jù)導(dǎo)致。

由于缺少有效的驗(yàn)證數(shù)據(jù)，本文采用人工判讀的方法對(duì)原始光子進(jìn)行標(biāo)注，將標(biāo)注后的結(jié)果與實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比分析，并用混淆矩陣對(duì)分類結(jié)果進(jìn)行定量分析。表1、表2 為用本文方法得到的去噪點(diǎn)云與真實(shí)點(diǎn)云對(duì)比混淆矩陣，表3 列出了試驗(yàn)數(shù)據(jù)的去噪精度。結(jié)果表明，本文提出的算法去噪精度高，去噪效果良好。

表1 點(diǎn)云去噪混淆矩陣（區(qū)域A）

表2 點(diǎn)云去噪混淆矩陣（區(qū)域B）

表3 去噪精度評(píng)價(jià)結(jié)果

4 結(jié)論

本文根據(jù)DBSCAN 算法公式計(jì)算MinPts 值，再結(jié)合K-平均最近鄰法選擇Eps，分離出淺海區(qū)域的信號(hào)光子，通過(guò)南海某島礁的ICESat-2 數(shù)據(jù)驗(yàn)證，得到以下結(jié)論：

（1）本文提出的算法能應(yīng)用于地勢(shì)平緩的淺海區(qū)域海面激光光子數(shù)據(jù)去噪，去噪精度達(dá)到98%，優(yōu)于傳統(tǒng)去噪算法精度，可為ICESat-2 海面數(shù)據(jù)的去噪工作提供參考；

（2） 受DBSCAN 聚類算法自身限制的影響，在光子數(shù)據(jù)稀疏及地勢(shì)變化劇烈的局部區(qū)域，存在錯(cuò)判的情況，從而影響去噪精度。針對(duì)此類情況，將在后續(xù)的研究中改進(jìn)DBSCAN 聚類算法的搜索形狀，將圓形搜索區(qū)域換成橢圓形區(qū)域沿水平方向搜索，以便更加精準(zhǔn)地獲取淺海區(qū)域的海面及海底光子，并剔除其噪聲影響。

致謝：武漢大學(xué)紀(jì)雪，山東科技大學(xué)石通為文章修改提供了幫助，在此一并致謝。