单光子激光雷达(single photon LiDAR，SPL)可将激光脉冲测距系统的探测灵敏度逼近理论极限值，为实现激光雷达远程探测提供了有效的技术手段^[1]。随着多波束分光技术、单光子阵列探测技术的发展，以多波束、高重频、微脉冲并行发射，阵列单光子探测为特征的新一代三维成像激光雷达可实现远距离、宽幅、高分辨、高帧频的三维立体成像，在冰盖高程测量、海冰厚度反演、水位监测、森林生物量估算、地形测绘等方面^[2-5]拥有广阔的应用前景。

由于探测系统的灵敏度高，太阳背景辐射和暗电流噪声极易触发探测电路，SPL与传统的激光雷达点云数据相比具有背景噪声大特点，在点云去噪方面具有较大的难度；且不同探测环境、不同系统参数也会导致SPL点云空间分布的差异，因此需要研究针对性的点云去噪算法，为后续数据应用奠定基础。

目前常用的单光子点云去噪方法可分为两类：一类是基于局部统计的方法，通过统计中心点与邻域点的距离以滤除噪声点；另一类是基于密度聚类的方法，以中心点一定邻域范围内点的数量判断噪声点。Herzfeld等^[6]提出一种基于空间分析和离散数学的处理方法，由高程-频率直方图划分地面与树冠数据，再由密度-频率直方图确定地面与树冠高度。夏少波等^[7]针对森林区域数据提出一种基于局部距离统计的去噪算法，计算点与其k近邻的距离之和，统计距离分布直方图，将局部距离之和大于阈值的点视为噪声点。Zhang和Kerekes^[8]介绍了一种基于密度的空间聚类算法——应用程序空间聚类(density-based spatial clustering of application with noise, DBSCAN)算法，根据整体点云的密度自适应地确定椭圆搜索邻域密度阈值，提取符合阈值条件的树冠与地面点云。Ma等^[9]提出一种基于改进DBSCAN的自适应光子信号提取方法，将点云在垂直方向上划分为多个区间，判断噪声区间与信号区间并分别计算密度，以此确定DBSCAN算法密度阈值，完成陆地与森林点云的提取。Popescu等^[10]提出一种结合了局部统计和聚类的方法，首先统计格网内点云识别可能的信号网格；然后利用聚类分析的方法，以正态分布作为光子的密度函数，通过最大化似然函数提取信号光子。

基于统计的方法对于密度变化具有较好的适应性，且在被测物体结构复杂时能保留细节信息，但无法消除与主体结构连接的成群点云；基于密度的方法则在大量离散点云远离主体点云，位置分散时能取得比较好的效果，适用于单光子数据背景噪声大的特点，但是该方法对于密度十分敏感^[11]。上述基于局部统计或基于密度聚类的单一方法大多用于植被和自然地表，而本文的研究区域地表类型丰富，包括树林、水面、房屋、裸地等多种类型，单一去噪方法存在普适性较差与精度低的问题。因此本文通过高程直方图确定可能的信号高程区间后，将基于局部统计与基于密度聚类的方法相结合，进行两次去噪：首先将基于密度的DBSCAN算法进行改进，使得密度参数由点云密度自适应确定，使用改进的DBSCAN(modified DBSCAN, MDBCAN)算法进行第一次去噪，去除大部分背景噪声；再使用一种基于局部统计的统计移除离群点(statistical outlier removal, SOR)算法去除地物附近的残余噪声，以获得更好的点云去噪效果。

1 SPL系统简介

本文涉及的机载多波束单光子三维成像激光雷达由激光器子系统、光学收发子系统、信号采集处理与综合控制子系统、机载高精度组合导航、一维扫描光电吊舱组成。系统主要参数见表 1。

其中，光纤激光器发射脉宽2 ns(FWHM)，重频40 kHz，单脉冲能量200 μJ。1 064 nm波长的高重频激光，经光学系统扩束后，通过衍射分光元件(diffractive optical element, DOE)分为64束，单波束发散角为50 μrad，相邻波束间隔为330 μrad，分束均匀度为92 %。波束在接收望远镜采用透射式镜头，在望远镜焦面处通过阵列光纤实现焦面阵列耦合接收，并将单波束接收视场限制在144 μrad以内。当机载平台飞行高度6 km时，垂直于机下点的单波束足印直径为0.3 m，接收视场地面直径为0.864 m，相邻波束中心点间隔为2 m，多波束形成的地面足印与接收视场如图 1(a)所示。阵列光纤为1分4光纤，望远镜端排布与发射光束在地面的足印对应，另一端为4束4×4排布光纤，每束4×4光纤阵列通过双远心镜头、窄带滤光片，将散射光窄带滤波后耦合至单光子探测器每一个像元。单光子探测器选用4个4×4铟镓砷盖格APD(InGaAs Gm-APD)实现激光回波微弱光子信号的接收探测。阵列光纤分束示意如图 1(b)所示。信号采集处理与控制单元控制各单机按流程与时序实现既定的工作模式，并实现单光子探测器输出的64路脉冲的飞行时间测量，测量精度优于100 ps，对应距离精度1.5 cm。

机载高精度组合导航获取雷达位置与姿态信息，并提供秒脉冲信号实现系统时间同步。一维扫描光电吊舱采用单轴高速扫描方式，垂直于飞行方向实现±60°的大范围指向调节和1 km范围内的高速扫描，并通过内部光栅实现高精度扫描角度的实时测量。信号采集处理与控制单元将测量数据与机载惯导数据、扫描角度打包存储。此外，系统还配备一台可见光相机，可实现激光三维点云与可见光复合成像。

根据上述系统参数，针对不同地物目标，在能见度为23 km条件下，系统的回波光子数预估结果如表 2所示。

2 数据处理

本文的单光子激光雷达共有64通道，单通道数据可投影为二维剖面点云，横坐标为点云在飞行方向上的位置，纵坐标为点云高程。以单通道数据为例说明点云去噪方法与流程：首先将数据沿飞行方向分段，统计每段点云数据的高程直方图，以此确定信号点云的高度区间，再对可能的信号点云分别采用基于密度的MDBSCAN算法与基于局部统计的SOR算法对点云进行粗去噪和精去噪。对64个单通道数据进行去噪处理后，再根据探测器的通道排布及POS数据生成三维点云。点云去噪的技术路线如图 2所示。

为减少MDBSCAN与SOR去噪算法的计算量，使用高程直方图确定信号点云高度区间，初始垂直方向距离窗长度100 m，将点云数据沿飞行距离分段后，统计段内光子的高程直方图，将分布曲线视为混合高斯函数，计算高斯分量中中心的最大值μ₁和最小值μ₂，以及对应的标准差σ₁和σ₂，仅保留高度区间在[μ₁-6σ₁, μ₂+6σ₂]内的光子点云。

2.1 MDBSCAN算法点云粗去噪

单光子激光雷达的噪声点云在空间中呈现随机分布特点，与地物点云相比密度更低，因此采用基于密度的方法可以有效区别单光子数据中的噪声点与地物点。DBSCAN算法是一种基于密度的聚类方法^[12-13]，在机载与星载单光子激光雷达数据处理中已得到广泛应用^[8-9]。该算法将具有足够高密度的区域划分为簇(有效信号)，并且能在噪声的空间中发现任意形状的簇。DBSCAN有2个重要参数，分别为半径参数R和密度参数MinPts，往往需要经过多次实验确定，且全局参数可能存在局部不适用问题。同时，DBSCAN的搜索邻域为圆形，不符合实际地物分布情况。因此本文从搜索邻域形状与密度参数两方面对DBSCAN算法进行改进。

参考文献[8]的思路，首先将DBSCAN算法的搜索邻域由圆形改进为椭圆形，以适应点云在水平方向密度大于竖直方向密度的特点。则点p、q的距离计算方法如下

$\operatorname{dist}(p, q)=\left[\frac{\left(x_p-x_q\right)^2}{a^2}+\frac{\left(h_p-h_q\right)^2}{b^2}\right]^{\frac{1}{2}} \text {, }$

(1)

其中：x是光子在飞机飞行方向上的位置，h为光子对应的高程，a、b分别为椭圆的长短半轴，且a>b。若dist(p, q)≤R，则q是p的R邻域内一点。

不同的下垫面会带来不同点云密度与背景噪声率^[14]，因此飞行航线上的点云不可避免地会产生密度差异。为降低算法的复杂程度，a、b设为固定值，MinPts采用一种自适应的方法计算产生。具体方法如下：

1) 对于每一个数据段，计算点云密度

$\rho=N \cdot S /(H \cdot L),$

(2)

其中：ρ为粗略估计的椭圆邻域内点云数量，N为该段数据内点云数量，S为搜索邻域的面积，S=πab，H为该段数据的垂直方向距离窗长度，L为飞行方向距离。

2) 计算密度阈值

$\text { MinPts }=n \cdot \rho,$

(3)

其中：n为阈值。沿飞行距离分段的数据段内地物类型基本一致，点云密度差异较小，由点云密度计算出每段的MinPts，以适应段间点云密度差异，代入搜索邻域为椭圆的MDBSCAN算法中，完成数据的第一次去噪。

2.2 SOR算法点云精去噪

MDBSCAN算法能滤除大部分的背景噪声，但是对于密度较大的噪声及地物附近的噪声滤除效果较差，因此本文又引入一种基于局部统计的去噪算法——SOR算法^[15]以进一步获得精细的去噪结果。

SOR算法假设噪声点之间的距离更远，统计点云中每个点的邻域距离均值，并假设所有结果应满足高斯分布，且高斯分布的形状由均值和标准差决定，则邻域距离均值在标准范围之外的点为离群点。具体方法如下：

1) 假设点云数量为n，计算点云中每个点与最近邻K个点的距离均值r_{ave, i}，i = 1, 2, ⋯, n，得到集合r_ave={r_{ave, 1}, r_{ave, 2}, ⋯, r_{ave, n}}；

2) 计算r_ave的均值mean与标准差std，再次遍历点云中的所有点，若r_{ave, i}>mean+m·std，则第i个点为噪声点，从点云中删除，否则保留。

文中的相应参数都是参考相关文献和经过多次实验确定的，具体数值见3.1节。

3 实验结果与对比分析 3.1 实验过程

飞行实验地点为湖北荆门，时间为2021年8月某日下午，飞机飞行速度约50 m/s，高度约1 km，飞行航线及实验区影像见图 3，图中红线为飞机航线，飞行方向从左至右。以推扫模式下单通道的数据为例，使用本文方法对光子点云数据进行去噪：将点云在剖面投影后按照200 m的飞行距离划分数据段，对于每段数据仅保留可能的信号高程区间内的点云；采用2.1节的方法进行粗去噪，椭圆邻域长半轴a取2 m，短半轴b取1 m，密度阈值中MinPts的n取1.8；采用2.2节方法进行精去噪，k值设为9，阈值m设为0.9。为证明该方法的有效性，将本文基于密度与局部统计的去噪方法(MDBSCAN+SOR)与MDBSCAN、原始DBSCAN+SOR进行对比。原始DBSCAN+SOR方法中DBSCAN算法搜索半径R为2 m，密度参数MinPts设为20，不同方法之间其余共有的参数取值相同。

3.2 结果与分析

选取典型地物的点云去噪结果进行分析，不同方法之间的对比如图 4、图 5所示，分别对应图 3①处厂房和②处树木，黑色点为噪声点，红色点为信号点。图 4(a)、5(a)中都存在大片的背景噪声，且有部分信号点云被判为噪声点云，说明原始DBSCAN算法不能很好地适用于单光子点云去噪。MDBSCAN方法对于背景噪声的去除效果较好，但是图 4(b)、5(b)的屋顶与树木冠层周围都围绕着明显的噪声点云，证明单纯的密度聚类方法仅能有效识别远离地物的背景噪声点云。图 4(c)、5(c)显示本文方法对背景噪声有良好的去除效果，地物点云完整，结构清晰，地物周围的噪声点相比于前2种方法都有所减少，在3种方法中去噪效果最好。

为定量地对比3种方法的去噪结果，对图 3中标注③处的路面点云进行直线拟合，计算路面点云与拟合直线的均方根误差，结果如表 3所示。从表中可以看出，原始DBSCAN+SOR的均方根误差最大，去噪效果最差，MDBSCAN方法次之，而本文方法的去噪效果最好，均方根误差仅为0.27 m。对3种方法在整条航线上单通道数据的计算时间进行比较，结果如表 3所示，MDBSCAN方法因只使用了单一算法，耗时仅为38.85 s；原始DBSCAN+SOR方法时间最长，达到73.53 s；而本文方法耗时61.63 s，相比于MDBSCAN方法增加22.78 s，但少于原始DBSCAN+SOR方法。同时以人工判别的信号点云作为标准，进一步计算3种方法在图 4、图 5中的准确率Acc、召回率R、精确率P，结果如表 4所示。3种方法按照各项指标从低到高排序大致为原始DBSCAN+SOR、MDBSCAN与本文方法，与路面点云直线拟合的结果相同，仅有MDBSCAN方法的召回率R高于本文方法，因为该方法少了1次去噪，使得被误判为噪声点的信号点减少，但是总体判断正确的点云少于本文方法。

4 结论

本文针对自研的机载单光子激光雷达系统提出一种基于密度与局部统计的二维剖面单光子点云去噪方法：在使用高程直方图确定信号点云的大致高程区间后，提出一种基于密度的MDBSCAN算法，将算法密度参数改进为由点云密度自适应决定，完成点云的粗去噪；再使用基于局部统计的SOR算法进一步去除密度较高的背景噪声与地物附近噪声。实验证明该方法对于不同场景和不同密度的点云均具有良好的适应性，去噪效果优于常规的单光子点云去噪方法，能够满足国产单光子激光雷达获取高精度地表三维轮廓的技术需求。