基于改进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 聚类算法的搜索形状，将圆形搜索区域换成椭圆形区域沿水平方向搜索，以便更加精准地获取浅海区域的海面及海底光子，并剔除其噪声影响。

致谢：武汉大学纪雪，山东科技大学石通为文章修改提供了帮助，在此一并致谢。