曹彬才，王建荣，胡燕，吕源，杨秀策
(1.地理信息工程国家重点实验室，西安 710054；
2.西安测绘研究所，西安 710054)
美国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.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.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去噪效果
DBSCAN作为典型的基于空间密度的聚类算法，具有实现简单、分类效果好、能实现任意形状的点聚类等特点。缺点是DBSCAN算法对核心参数MinPts和Eps非常敏感，两个参数直接决定聚类效果。本文研究了多场景下ICESat-2光子点云DBSCAN去噪效果。实验表明，不同场景下的光子点云特征差异较大，不宜使用一个固定参数处理多种场景。针对不同目标类型，以最终聚类数为参考虽然可以自适应确定DBSCAN算法的Eps参数，实现全自动化去噪，但这种思路缺乏光子对光子密度的针对性设计，计算复杂度高，耗时较长，与现有的主流算法相比效率较低。下一步应针对不同地物类型在不同高度的点云密度变化特点，设计自适应滤波模型，力求在处理效率上得到大幅改进。
猜你喜欢光子聚类噪声噪声可退化且依赖于状态和分布的平均场博弈数学年刊A辑(中文版)(2020年3期)2020-10-27基于K-means聚类的车-地无线通信场强研究铁道通信信号(2019年6期)2019-10-08偏振纠缠双光子态的纠缠特性分析电子制作(2019年12期)2019-07-16控制噪声有妙法中学生数理化·八年级物理人教版(2017年9期)2017-12-20基于高斯混合聚类的阵列干涉SAR三维成像雷达学报(2017年6期)2017-03-26基于Spark平台的K-means聚类算法改进及并行化实现互联网天地(2016年1期)2016-05-04基于改进的遗传算法的模糊聚类算法智能系统学报(2015年4期)2015-12-27光子嫩肤在黄褐斑中的应用中国医疗美容(2015年2期)2015-07-19一种基于白噪声响应的随机载荷谱识别方法噪声与振动控制(2015年4期)2015-01-01车内噪声传递率建模及计算振动、测试与诊断(2014年4期)2014-03-01