孟文君，李杰，張凱，2，劉長達，唐秋華，

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

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

針對光子去噪問題，目前主要有以下三類方法：（1）基于柵格圖像處理的去噪算法（Chen et al，2015； Magruder et al，2012），該算法原理簡單，但在點云和柵格圖像的轉(zhuǎn)換中會損失部分精度；（2）基于局部統(tǒng)計參數(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），該算法通過利用噪聲光子在空間中分布離散的特點進行聚類分析，從而剔除噪聲光子。主要的密度聚類方法有貝葉斯（李明等，2018）、DBSCAN（Density-Based Spatial Clustering of Applications with Noise）和OPTICS（Ordering Points to Identify the Clustering Structure）。在這三種方法中，DBSCAN 算法原理較為簡單，也是目前針對光子去噪最常用的手段之一。

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

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

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

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

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

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

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

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

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

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

2 研究方法

2.1 DBSCAN 算法原理

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

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

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

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

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

最后，就可以通過公式（3）計算MinPts 值

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

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

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

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

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

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

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

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

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

2.3 精度評價指標(biāo)

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

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

3 結(jié)果與分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

表3 去噪精度評價結(jié)果

4 結(jié)論

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

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

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

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