密度聚类算法在光子点云去噪中的应用与评价

曹彬才，王建荣，胡燕，吕源，杨秀策

(1.地理信息工程国家重点实验室，西安 710054；2.西安测绘研究所，西安 710054)

0 引言

美国2018年9月发射了第二代对地观测激光雷达卫星ICESat-2(ice，cloud and land elevation satellite-2)，旨在利用激光雷达高精度特点在全球范围内开展冰盖变化监测、海面高度测量、植被覆盖反演等，为研究碳循环和全球变暖等科学问题提供技术支撑[1]。与此同时，ICESat-2高精度点云数据在高精度地形测量、激光点云辅助遥感影像平差、浅海水深测绘等[2-4]传统测绘邻域也得到了广泛应用，ICESat-2点云数据处理成为测绘行业的研究热点。

ICESat-2搭载的激光载荷ATLAS(advanced topographic laser altimeter system)采用了光子计数探测体制，这种新型激光雷达具有功耗小、重量轻、灵敏度高等特点[5]，比传统线性探测激光雷达更适合用作星载平台。高灵敏度探测带来的缺点是数据噪声多，ATLAS接收器会将大量的太阳辐射光子和大气散射光子记录为光子事件，导致数据信噪比差，因此点云去噪对光子计数激光雷达至关重要。

当前光子点云去噪总体上有两种思路[6]：一是基于数字图像处理技术，将剖面点云转换为影像，采用边缘检测[7]、区域检测等识别噪声；二是逐点计算某个局部统计量，利用其分布特征计算全局阈值并区分信号和噪声，如直方图法[8]、点云聚类法[9-10]、随机森林法[11]等。其中DBSCAN(density-based spatial clustering of applications with noise)是一种较为成熟的聚类算法，能将半径Eps范围内点个数超过最小值MinPts的点归为一类，广泛用于激光点云去噪。Zhang等[12]将密度均值作为MinPts，并将圆形搜索核改进为平行地面的椭圆，处理有植被区域的MABEL(multiple altimeter beam experimental LiDAR)单光子数据时效果较好；Ma等[13]在处理浅海测深数据时使用DBSCAN进行去噪，给出了海洋场景下MinPts经验公式，半径Eps则采用经验值，在浅海区域也能够有效识别噪声；李文杰等[14]通过数据集本身分布特征生成Eps候选集，通过穷举方式确定最佳参数；魏硕等[15]通过k维树求取点云密度进行粗去噪，然后运用改进DBSCAN算法和统计滤波算法进行精去噪。本文详细论述了当前光子计数激光雷达去噪中DBSCAN具体应用方法，探讨了关键参数MinPts、Eps自适应确定的可行性，并使用不同地物类型ICESat-2数据开展精度验证。

1 DBSCAN去噪方法原理

1.1 DBSCAN算法

DBSCAN算法是一种基于密度的空间聚类算法，通过寻找密度相连的点的最大集合来分离信号点和噪声点[16]。对于集群的每个点，给定半径Eps邻域必须至少包含最少数量的点，即邻域中的密度必须超过阈值MinPts。DBSCAN能自动将密度足够大的点区域划分为簇，数据集中不属于任何簇的点则被视为噪声。该算法原理简单、运算速度快，且不需要预先指定簇的个数，核心问题是要输入半径Eps和邻域最小点数MinPts两个基本参数。

如图1所示，ICESat-2光子点云常被投影到沿飞行方向和高程方向组成的二维平面，X轴表示沿飞行方向距离，可通过数据文件中的速度与时间参量相乘得到，Y轴表示高程。图中可明显观察到位于地形线附近点云更加密集，地形线上下方噪声点在空间分布上更加稀疏，这也是基于密度算法成功识别信号和噪声的必要条件。

图1 ICESat-2典型光子计数点云剖面图

1.2 参数确定方法

1)邻域最小点数MinPts。经验公式一：文献[12]通过定义目标点范围内密度值来确定算法MinPts。

(1)

MinPts≥4·ρ

(2)

式中：ρ表示密度值；N是数据集点总数；S为数据集二维剖面面积；s1表示局部搜索范围面积，s1=π·Eps2。文献[10]在处理MABEL数据时将Eps设为固定值2，并根据式(1)、式(2)计算MinPts。

经验公式二：文献[13]在处理海洋场景浅海测深点云去噪时，给出以下MinPts经验公式。

(3)

(4)

(5)

式中：S′、N′分别是某背景范围面积和对应的点总数。选择高程Y轴最低的5 m范围作为背景，沿轨方向X保持不变，仍为S对应的沿轨长度。

2)半径Eps。现有去噪算法处理ICESat-2数据时搜索半径基本都采用经验参数[17-18]，文献[14]提出了一种基于K平均最邻近法寻找最优Eps参数的思路，步骤如下。

步骤1：根据K平均最邻近法求出数据集D的候选Eps参数集合DEps。首先计算每个点的K最邻近点对应的距离(K=1，2，3，…，N)，形成距离矩阵DN×N，该矩阵中每一行表示一个点，每一列表示从小到大排列的K最邻近距离，每一列平均值即为K平均最邻近距离。

步骤2：在已知MinPts下，将DEps中Eps参数逐个带入DBSCAN算法中进行聚类运算，当生成的簇数连续3次相同时认为结果趋于稳定，簇数M为最优簇数。

步骤3：继续执行步骤2直到簇数不等于M，选用簇数为M时的最大K平均最邻近距离作为最优Eps参数。

3)本文改进方法。图2展示了图1数据Eps列表与K值的关系，随着K值增大，Eps参数随之增大。按照上述Eps寻优计算思路，本文通过对多类型ICESat-2数据的实验与分析，总结出DBSCAN聚类个数与K值关系(图3)具有下述特点。

图2 参数Eps列表与K的关系示意图

图3 典型场景下聚类数与K的关系

①植被、城市等陆地场景下，聚类个数与K值关系线呈现“双峰”分布式，即K值从小到大，聚类数先变大，再减小并趋于平稳，随后再次变大，最后减小直到稳定不发生变化。

②冰盖海洋等场景下，聚类个数与K值关系线呈现“单峰”分布式，即聚类数随着K值的增大而变大，随即减小趋于平稳。

同时，本文通过试错观察等方法表明：最佳K值如图3中红点所示位置，第一种情况位于两波峰之间，第二种情况位于波峰右侧。因此本文改进Eps自动确定的步骤如下。

步骤1：利用MinPts值限制最大K值范围，Kmax=D·MinPts，Kmax表示最大K值，D=15，此范围内足以使得聚类数趋于稳定，同时避免计算所有K值，大幅减小计算量。

步骤2：生成聚类数与K值关系曲线，利用一阶求导等方法判断波峰个数。当曲线呈现双峰分布时，取波峰之间聚类数稳定的最小K值为最佳K，对应Eps为最佳搜索半径；当曲线呈单峰分布时，取波峰右侧聚类数稳定的最小K值为最佳K。

2 实验和分析

2.1 ICESat-2单光子数据

ICESat-2卫星轨道高度约500 km，受传播过程中的大气散射、目标漫反射等影响，一束激光脉冲仅有数个或数十个光子能返回接收器。不同拍摄条件、地物类型对应点云的信噪比和点云疏密程度有较大差异，ICESat-2在官方文档泊松去噪算法里将地物类型划分为陆地冰、海冰、海洋、陆地、内陆水5种类型，分别设计经验参数[19]。

本文选择了植被、城市、冰盖、海洋4种不同类别的光子数据开展实验，数据详细信息见表1。Data1至Data8的拍摄地点分别为美国新罕布什尔州、陕西城固市、拉斯维加斯、西安市、格陵兰岛、南极大陆、南海珊瑚礁、夏威夷海岸，均选择了强波束开展实验。

表1 实验数据详细信息

2.2 实验过程及结果分析

实验中利用经验公式二计算MinPts，由于文献[14]提出的计算全部K均值思路计算量过大，非常耗时，因此使用改进后的Eps确定方法。为便于分析，这里对比固定参数DBSCAN(记为DB1)、本文自适应DBSCAN(记为DB2)算法以及官方推荐算法DRAGANN，实验主要对算法参数的自适应确定、计算耗时以及去噪精度等方面开展评价。

1)参数自适应确定。表2为算法参数自适应确定及耗时情况，表中MinPts和Eps为自适应算法DB2计算的最优参数。本文程序采用VS2010 C++语言编写，并使用激光点云处理库PCL用作构建KDtree并寻找最邻近K值，自适应参数算法中采用CPU多核并行运算。

表2 参数确定及耗时

从表2可知，Data1、Data2两处植被区的MinPts和Eps区域对比而言差异不大，分别为11、11.5和13、12.5，固定参数时DB1计算时间仅仅为42″和55″，而自适应算法DB2时间达到35′54″和68′32″，效率很低。

Data3、Data4两处城市区的核心参数区域对比而言差异较大，Data3的MinPts和Eps分别为12、7.3，Data4为6、14.9，这是由于Data3的点密度更大，并且美国该城市区域为低矮建筑，Data4西安城市点密度稍小，且建筑高度较大，因此需要更大的Eps才能保证将地形和建筑物这种在剖面上断开的对象聚类为一类。

Data5、Data6是冰盖高反射区域，平均一束激光返回的光子数要远大于植被，因此MinPts比前几组数据都要大，分别为58和30。Data5位于格陵兰岛平地区域，点密度较大，在搜索半径Eps为1.97 m时即可有效区分信号和噪声。Data5位于南极某山地，冰雪覆盖略少，且存在起伏断裂，因此自适应Eps稍大，为11.2 m。

Data7、Data8是海面区域，反射光子低于冰盖区，与陆地相差不大，两组数据的MinPts和Eps分别为14、7.4和6、6.1。

2)计算耗时。光子点云投影到飞行方向剖面后，实际上成为二维数据，DBSCAN算法的空间复杂度为O(n2)，n表示点个数，循环生成Eps列表空间复杂度为O(n)，因此总的空间复杂度为O(n2)+O(n)。

2.3 去噪精度

本文自适应DBSCAN算法去噪精度如表3所示，8组数据的整体去噪精度分别为97.3%、97.6%、99.2%、98.3%、97.2%、98.6%、97.9%和99.1%，均表现较为优异，整体水平与泊松去噪、DRAGANN去噪等几乎一致。图4、图5、图6分别为Data1、Data3、Data5去噪效果，从目视观察来看，本文DBSCAN算法均有效识别出信号，Data1的植被和地形线被正确识别，Data3城市区域没有因为地形不连续而出现识别错误，Data5的冰面点云密集度高，也被正确识别为信号点。不足之处在于地形线下方有明显的噪声点被识别为信号，这是空间密度类算法的固有缺陷。以Data5为例，最佳Eps=2.9 m时，地形线下方2.9 m以内的噪声点会被识别为信号，而冰面处地形平坦，点密度高，地形线下方的噪声点在视觉上很容易发现，但自动算法却分类错误。

表3 实验数据自适应DBSCAN去噪混淆矩阵

图4 Data1植被DBSCAN去噪效果

图5 Data3城市DBSCAN去噪效果

图6 Data5冰盖DBSCAN去噪效果

3 结束语

DBSCAN作为典型的基于空间密度的聚类算法，具有实现简单、分类效果好、能实现任意形状的点聚类等特点。缺点是DBSCAN算法对核心参数MinPts和Eps非常敏感，两个参数直接决定聚类效果。本文研究了多场景下ICESat-2光子点云DBSCAN去噪效果。实验表明，不同场景下的光子点云特征差异较大，不宜使用一个固定参数处理多种场景。针对不同目标类型，以最终聚类数为参考虽然可以自适应确定DBSCAN算法的Eps参数，实现全自动化去噪，但这种思路缺乏光子对光子密度的针对性设计，计算复杂度高，耗时较长，与现有的主流算法相比效率较低。下一步应针对不同地物类型在不同高度的点云密度变化特点，设计自适应滤波模型，力求在处理效率上得到大幅改进。