基于改進DBSCAN 算法的ICESat-2 海面點數據去噪處理及精度評估

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

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

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

針對光子去噪問題，目前主要有以下三類方法：（1）基于柵格圖像處理的去噪算法（Chen et al，2015； Magruder et al，2012），該算法原理簡單，但在點云和柵格圖像的轉換中會損失部分精度；（2）基于局部統計參數的去噪算法（Gwenzi et al，2016；謝峰 等，2017；Herzfeld et al，2019；夏少波等，2014），該算法有較好的適應性，但受閾值影響較大，如何選取有效閾值仍有待研究；（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 算法對輸入參數較為敏感，如何選擇合適的參數成為影響光子去噪效果的關鍵因素。秦磊等（2020）通過計算Eps 取0～30 內的整數時的F1原score，得到最佳Eps 值，從而進行DBSCAN 聚類，但這種方法具有較大的偶然性；馬躍等（2020） 提出了一種計算淺海中MinPts（Minimum points） 參數的改進公式，聚類效果良好，但并未說明Eps 的計算方法。

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

1 研究區域與數據

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

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

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

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

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

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

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

2 研究方法

2.1 DBSCAN 算法原理

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

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

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

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

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

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

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

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

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

2.2 自適應確定最優參數

自適應確定最優參數的步驟如下：

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

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

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

步驟四： 繼續執行步驟三，直到循環結束，并選用當簇數為2 時所對應的最小K 值作為最優K 值。最優K 值對應的平均距離則為最優Eps 參數。

2.3 精度評價指標

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

式中，TP 是被正確分類為信號點的光子數，FP 是實際為噪聲但被錯分為信號的光子個數，FN 是實際為信號點卻被錯誤的地分為噪聲的光子數。

3 結果與分析

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

圖3 原始數據所有光子（區域A）

該實驗區域的全部光子按公式（3）參與計算MinPts 值，結果為MinPts =18。按K-平均最近鄰法分別計算4000 個光子在每個窗口內的Eps 值。得到4 個聚類簇數與K 值的關系圖，如圖4 所示。可以看出簇數穩定在2 時的最小K 值分別為50、68、56、51。圖5 展示了K 值與Eps 的對應關系，最小K 值對應的最佳Eps 分別為1.67、1.89、1.53、1.43。用實驗獲得的參數分別對每個窗口內的點聚類，聚類結果如圖6 所示。

圖4 K 值與分類個數（區域A）

圖5 K 值與Eps（區域A）

圖6 各組分類情況（區域A）

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

圖7 本文方法去噪結果（區域A）

圖8 ATL03 算法提取的官方去噪結果（區域A）

圖9 本文方法去噪結果（區域B）

圖10 ATL03 算法提取的官方去噪結果（區域B）

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

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

表1 點云去噪混淆矩陣（區域A）

表2 點云去噪混淆矩陣（區域B）

表3 去噪精度評價結果

4 結論

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

（1）本文提出的算法能應用于地勢平緩的淺海區域海面激光光子數據去噪，去噪精度達到98%，優于傳統去噪算法精度，可為ICESat-2 海面數據的去噪工作提供參考；

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

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