Single photon denoising algorithm combined with improved DBSCAN and statistical filtering

近些年, 单光子雷达技术逐渐发展成熟^[1-3]，作为一种新型激光3维雷达，光子计数激光雷达拥有自己显著的特点，如回波信号强度低、噪声大于激光回波、光电探测系统处于光子水平等^[4]。不同于传统雷达探测回波时需要较多的光子形成弱光信号，单光子雷达使用单光子敏感探测器，仅需要检测到若干光子的返回信息便可得到目标位置信息，使得发射的光子更有效被利用，从而减少激光发射脉冲所需能量，因此单光子探测技术在采集数据方面具有较大优势，在当下以及未来，光子计数激光雷达在侦察、搜索等方面会发挥越来越重要的作用。

光子计数型激光雷达引入单光子探测器件和时间相关单光子计数(time-correlated single-photon counting，TCSPC)技术，使系统具有探测极微弱光信号的能力和皮秒量级的时间分辨率^[5]，即便在系统设计中采用超窄带滤波，但由于系统本身的灵敏度极高，探测目标背景仍会产生大量的噪声，因此，如何有效地去除噪声是发挥光子计数雷达探测能力的关键。目前常用的噪声去除和数据提取算法主要有残差比较算法、泊松滤波算法和距离相关检测算法等^[6-8]。

HORAN等人^[9]根据扫描数据中信号光子和噪声光子的分布特征，结合地形设定阈值，从而分离出目标点。HERZFELD等人^[10]为模拟的ICESat-2数据分析而开发的计算算法，主要是借助特征向量计算出各项异性的密度中心，再根据地面和树冠设定密度阈值，达到识别地物的目的。LI等人^[11]对初步去噪后的点云使用改进后的统计滤波方法，设定阈值进行去噪，但参量的设定比较依赖经验，对结果也有着较大的影响。LI等人^[12]依靠计算点云密度进行去噪，首先进行图像分割进行粗去噪，再定义一种线性径向基函数进行精去噪，得到结果与设定的阈值比较，得出较好的结果。XU等人^[13]在对单光子激光数据特点精确分析的基础上，提出一种基于地形相关和最小二乘曲线拟合的数据处理方法，获得较为精确的结果。上述算法要求的参量较多，且参量数值选取有较高的要求，作者在研究单光子数据空间分布特征的基础上，提出一种基于密度、参量自适应较强的去噪算法，并通过美国国家航空航天局(National Aeronautics and Space Administration，NASA)提供的多波束试验激光雷达(multiple altimeter beam experimental lidar, MABEL)实际飞行数据^[14-15]进行验证，得到较好的结果。

光子计数激光雷达^[16]的系统框图如图 1所示。

Figure 1.  Photon counting lidar system

如图 1所示，光子计数激光雷达的工作原理是：脉冲激光器发射的脉冲经分光片后分成2束，一束充当计时脉冲照射到光电探测器，触发计时信号；其余部分照射到物体，经反射后产生回波被单光子探测器接收，并在计时电路中产生停止信号。TCSPC模块计算2次信号的时间差，得到本次脉冲的来回飞行时间。随着探测时间和探测次数的增加，得到扫描光子随时间分布的数据，传递给后续模块通过时间获取算法得到具体目标距离数据及成像。

由于系统本身的灵敏度很高，在探测过程中目标背景会产生大量的噪声，激光雷达的噪声大致分为3类：探测器噪声、放大器噪声、背景的辐射噪声^[17]。随着制造工艺和致冷方式的进步，探测器噪声和放大器噪声得到有效抑制，为抑制背景噪声，除了在接收光学系统添加窄带滤波片外，还会采用距离门技术。

距离门技术是指给探测器添加一个时间区间，只有在这个时间区间内返回的回波光子才可以被探测到，在这个时间段外的回波光子不会被探测到。虽然如此，但大量的背景噪声仍会被引入到探测数据中，这就需要更近一步的去噪^[18-19]。

常规的点云密度定义D_i是指在规定的半径eps内包含点云的数目，如下式所示：

式中, p_i是中心点坐标，p_j是临近点坐标，i, j=1, 2, 3…; eps是指基于密度的空间聚类噪声应用(density-based spatial clustering of applications with noise，DBSCAN)算法中的半径参量。

与常规的点云密度定义相反，本文中规定点密度x_i表示单个点云的密度，如(2)式所示。由k维树可以得到k个临近点到中心点的距离数组{x₁, x₂, …, x_k}，由于k值是输入参量，因此点密度x_i表示包含k个点云的最小距离。

分割阈值x_m根据点密度x_i求出，其大小为整体点密度的均值，计算下式所：

由光子计数激光雷达的特性以及对实际数据集的研究可知，噪声点云的数量远远大于信号点云，因此分割阈值的范围在0.8x_m~x_m之间均可以满足实验要求，经多次实验可知，分割阈值的取值为0.9x_m时结果最佳。

DBSCAN算法根据人为输入参量eps以及minpts(DBSCAN算法中的固定变量符号，表示eps距离内包含的个数) 衡量点云密度，遍历所有点云，将点云划分成一个个密度不等的簇。由于输入参量有较大的主观性，如果和数据集密度特征不匹配容易造成过分割或欠分割。

为了弥补原DBSCAN算法的缺点，本文中将点密度x_i用于衡量单个点云密度，由于点密度x_i和密度阈值x_m均由数据集本身密度特征自动计算，因此在核心点“生长”为簇的过程中，避免了人工带来的误差。

改进后的聚类算法由于具有密度参量自适应的特性，可以对绝大多数的点云集进行合理分割。

本文中去噪方法可以分为2步：粗去噪和精去噪。首先对数据进行粗去噪，去除大部分明显的噪声点，可以极大地简化后续算法的计算量；然后根据剩余点云的密度特征和空间分布特征，执行精去噪，获得目标点云数据。算法流程如图 2所示。

Figure 2.  The flow chart of noise filter method

(1) 坐标转换。实际实验数据是大地坐标系(L, M, H)，即由经度L、维度M和高程H组成，利用坐标转换^[20-21]，将大地坐标系数据转化成空间直角坐标系(x, y, z)数据；(2)建立空间拓扑关系。对直角坐标系下的点云数据通过k维树建立空间拓扑关系，由近邻索引计算点云之间密度；(3)设定阈值去除明显噪声。选取合适的阈值进行点云分割(本文中选取0.9倍密度均值作为阈值)，滤除明显的离散噪声点。

(1) 计算剩余数据的点密度x_i。对剩下的数据建立k维树，由最近邻索引得到距离数组{x₁, x₂, …, x_k}，求出点密度x_i；(2)根据点密度进行点云分簇。随机选取一个点，当该点点密度满足阈值x_m时进行区域增长，逐步遍历其周围的点，将满足条件的点归为一簇，重复此过程，直至所有点云均被处理；(3)去除团状噪声簇。数据分簇后，检查每个点云簇包含的点云数量，对于点云数小于一定值的点云簇视为噪声簇去除，具体分析如下：由于本文中选取的研究对象是山脉，在噪声影响的下，目标也可能被分成几段，所以经过试验选取阈值为0.5T，既可以达到去噪目的也尽可能保留信号点，其中，阈值T=数据集点云数÷聚类后簇数，查看每个簇点云数目，如果小于0.5T，则作为密集噪声簇去除; (4)去除目标很近的离散噪声。遍历点云得到临近的k个点, 并求出距离均值构成1维数组D，然后计算数组D的均值d_i和标准差σ，设定阈值滤除离散点。

定义d_i为点p_i(x_i, y_i, z_i)到其k个临近点p_ij(x_ij, y_ij, z_ij)的平均距离，d_i为d_i的均值，σ为d_i的标准差，j=1, 2, …, k。

本文中的算法，需要设定的参量只有k，k表示选取临近点的数目，k值大小的意义在于点密度的衡量标准。k值很小, 话点密度由最近的较少点决定，增加了特殊值出现的情况，如某个信号点恰好周围最近的几个点较稀疏，或者某个噪声点最近的几个点距离很近等；k值很大, 会使信号点与噪声点的密度分割线较为模糊，虽然可以额外保留分布特殊的信号点，但也容易保留较多的噪声点。

本文中的数据来源是美国NASA提供的MABEL实际飞行数据，目标区域是位于美国加州Sierras-Forest地区，采用2010年12月飞行试验的532nm波段第50通道的观测数据进行实验。坐标转化后的数据总共有169511个点，其中噪声数据点有123280个，目标点云有46231个。

本文中选取不同k值代入算法进行运算，去噪后的点云数据如表 1所示。

此外还借助2个参考标准衡量k值对去噪结果的影响。一个是滤波精度N，其公式如下：

式中，N_c表示正确区分的信号点数目，N_t表示信号点总数。滤波精度N主要描述了点云正确识别率，但没有考虑噪声的影响。基于此，本文中还采用交并比U来描述信号点与噪声点的关系，公式如下：

式中, A表示实际信号点云数据集，B表示去噪后点云数据集，(A∩B)表示的是去噪后保留的点云数，(A∪B)表示的是实际信号点云数据和去噪后保留的噪声数据之和。一般来说，在计算机检测任务中，U>0.5的结果就是可以接受的，U越大，则点云图的质量也就越好。具体数据如表 2所示。

从表 1可知，k值和去噪后的信号点和噪声点数目均呈正相关，这和上面对k值的分析相符。由表 2可知，随着k值的增长，N的增长率逐渐接近于0，而U却逐渐降低，表明噪声点的增长率大于信号点的增长率，继续增大k值得不偿失。总体上来看，U保持在75%~80%左右，表明k值的适当选取对点云图质量影响不是很大，符合本文中算法参量自适应的要求。

由第3节中分析可知，本文中的算法参量具有参量自适应强的特点，k值对最终结果的影响并不是很大，因此令k=23，首先将数据由大地坐标系转换到空间直角坐标系，并将坐标系原点平移到数据最小处，然后应用本文中的算法进行去噪。去噪后，剩余噪声点为6290个，保留的目标点为41733个，滤波精度N=90.27%，交并比U=79.5%。

图 3是坐标转换且进行平移后的原始数据。其中黑色是噪声，细线是山脉。从图中可以看出, 噪声点的数量远大于信号点数量。图 4是粗去噪的结果。从图中可以看出，大量的噪声点已经被去掉，可以近似看出山脉的轮廓。图 5和图 6是精去噪的结果图。其中图 5是利用改进DBSCAN算法进行二次去噪后的结果，图中极大部分明显的团状噪声点已经被去除。图 6是应用统计滤波处理后的最终去噪结果，目标主体附近的离散噪声点云也有了明显的去除。图 7是部分放大后的场景图。从中可以看出，目标细节有了更多的保留，山坡上的树木也可以分辨。

Figure 3.  3-D scene graph

Figure 4.  Rough denoising result graph

Figure 5.  Second denoising result graph

Figure 6.  Second denoising result graph

Figure 7.  Partially enlarged view

半径滤波算法也是常用的滤波算法，其原理是假定数据集中每一个点在给定的半径r内至少存在k个点，符合假定条件的作为信号点保存下来，不符合条件的作为噪声点除去，其中半径r以及数值k由人工指定。为了优化结果，事先选取数值k=30，半径r在密度阈值附近多个范围内取值，观察不同参量下的半径滤波结果，并与本文中去噪算法结果进行比较。

在k=30时，计算出的r=21。由于r值越大，保留的噪声点就越多，因此对于参量r主要选取小于21的值。令r为3, 6, 9, 12, 15, 18, 21和24，去噪后数据如表 3所示，去噪效果较好的结果如图 8所示。

Figure 8.  Filtering results under different r values

结合图 8和表 3可以发现，在不同的参量下，半径滤波的结果起伏非常大。半径滤波的准确性建立在较为苛刻的先验知识上，需要人为地设定参量，不具备普遍适用性，最好的参量情况下，其结果达到78.39%。

参考文献[12]中提出的去噪算法也采用滤波精度作为衡量算法的指标，该文献中滤波精度分别达到90.26%以及96.11%。本文中算法在输入参量k>20后，其滤波精度均达到90%左右，点云去噪结果较好，且可选参量范围较广，突出本文中算法参量自适应能力强的特点。

DBSCAN算法分割结果的好坏取决于参量设定，需要eps, minpts和数据集密度特征相匹配，因此，DBSCAN算法在处理不同数据集时需要多次实验来确定最佳参量，而本文中改进后的DBSCAN算法仅需要设定临近点数目，相应的距离参量可以由k维树根据数据集密度特征自动得出，因此改进的DBSCAN算法在适应性上大大提高。其次，本文中对改进后DBSCAN算法分割后的结果进行处理，去除了目标附近大部分团状噪声点云，再根据剩余点云特征，设定统计滤波法的阈值对点云进行最终处理，得到结果。

半径滤波算法是应用极为广泛的密度滤波算法，本文中利用半径滤波算法处理数据集，将得到结果与本文中算法相对比，可以看出, 本文中算法几乎不需要先验知识，在参量自适应性上远远优于半径滤波算法。此外, 本文中的去噪算法在参量适中的情况下, U普遍在75%~80%(半径滤波算法的U在20%~79%), 滤波精度在85%~93%，去噪结果也好于半径滤波算法。两种算法对比，更加突出本文中算法的优越性。

但本文中的算法也存在一定的缺点，对于距离目标较近且密度极大的团状点云无法有效的去除，只能尽量减少其大小与数量，此外，如果数据集很大的情况下，k维树的搜索速度也逐渐下降。由图 5粗去噪结果图可以清晰地看出扫描线结构，在粗去噪的基础上，结合扫描线间距以及同一扫描线上点间距进行精确去噪，这是作者下一步进行改进的方向。