一种改进的基于聚类的单光子点云数据去噪方法与流程

本发明属于单光子探测的信号处理领域，具体涉及一种改进的基于聚类的单光子点云数据去噪方法。

背景技术：

单光子探测技术是近几年发展起来的新型探测技术，以目前已得到应用的计数星载激光测高仪为例，其采用单光子探测器件作为接收器件，灵敏度比传统激光测高仪提高了2-3个数量级，更易实现微脉冲、高重频、多波束的直接三维成像，具有很大的技术优势和发展前景。光子计数星载激光测高仪测量频率高，数据量大，探测地表类型多种多样，因而研究快速、有效、对多种目标都具有良好鲁棒性的点云数据去噪方法对光子计数激光测高仪数据处理具有重要意义。现有的去噪办法包含基于栅格，基于局部统计信息和基于聚类的三类方法。其中，基于栅格的方法在目标起伏比较厉害或坡度较陡峭时容易出现误判；基于部统计信息的方法时间复杂度高，运算时间久，且由于实际地形多变，密度阈值选取困难；基于聚类的方法受点云分布不均匀的影响较大，固定的阈值难以同时对反射率和环境差异较大的目标正确提取。

技术实现要素：

针对背景技术中存在的问题，本发明提供了一种改进的基于聚类的单光子点云数据去噪方法。该方法通过粗去噪简化点云数据，减少计算量；通过点云均匀化平均点云密度，防止因反射率和环境差异造成的局部点云分布不均，噪声密度大于或接近信号密度的现象；通过精去噪中的自适应dbscan算法对均匀化后的点云进行聚类和目标信号提取。本技术方案有运算速度快，去噪效果好，信号提取准确等突出特点。

本发明采用的具体技术方案如下：

一种改进的基于聚类的单光子点云数据去噪方法，包含：

步骤1：将光子计数激光测高仪的三维点云数据转换为二维点云数据，二维点云数据通过粗去噪方法得到粗去噪后二维点云数据；

步骤2：将粗去噪后二维点云数据通过均匀化处理得到均匀化后二维点云数据；

步骤3：将均匀化后二维点云数据通过精去噪方法得到精去噪后二维点云数据；

作为优选，步骤1中所述三维点云数据转换为二维点云数据为：

将光子计数激光测高仪的三维点云数据p＝(mk,nk,yk)用二维点云数据表示：

p＝(xk,yk)，k∈[1,nt]

其中，nt为三维点云离散点的数量即二维点云离散点的数量，mk表示第k个三维点云离散点的经度，nk表示第k个三维点云离散点的纬度，xk表示第k个二维点云离散点的水平距离值即第k个二维点云离散点的沿轨距离值，且yk表示第k个二维点云离散点的垂直距离值即第k个二维点云离散点的高程；

步骤1中所述粗去噪方法为：

将二维点云数据按照每lv米为一个间隔，划分成m＝h/lv个高程片，m为高程片的数量，高程的范围为：

h＝max(yk)-min(yk)，k∈[1,nt]

其中，nt为二维点云离散点的数量，max(yk)为yk，k∈[1,nt]的最大值，max(yk)为yk，k∈[1,nt]的最小值；

每个高程片内二维点云离散点为：

pi＝(xi,j,yi,j)，i∈[1,m]，j∈[1,ei]

lv×(i-1)＜yi,j＜lv×i

其中，xi,j为第i个高程片内第j个二维点云离散点沿轨距离值，yi,j为第i个高程片内第j个二维点云离散点高程，ei为第i个高程片内二维点云离散点的数量；

以高程片i为横坐标，第i个高程片内二维点云离散点的数量ei为纵坐标建立高程统计分布直方图，并用巴特沃斯滤波器对高程统计分布直方图进行五点加权滤波：

ei,filt＝a1ei-2+a2ei-1+a3ei+a4ei+1+a5ei+2i∈[3,m]

a＝(a1,a2,a3,a4,a5)

其中，ei,filt为滤波后第i个高程片内二维点云离散点的数量，a1为滤波器第一加权系数，a2为滤波器第二加权系数，a3为滤波器第三加权系数，a4为滤波器第四加权系数，a5为滤波器第五加权系数；

根据滤波后高程片内二维点云离散点的数量得到高程阈值：

ed,filt＝max{e1,filt,e2,filt,...,em,filt}

et＝ed,filt/2

其中，max为求最大值，et为高程阈值；

当高程片(i，ei,filt)i∈[1,m]满足以下条件：

ei<et,i∈[3,m]；

(i，ei,filt)为非极小值点，即ei,filt至少要大于ei-1,filt，ei+1,filt中的任意一个；

lv×|i-d|>＝λ，d为ed,filt对应的高程片；

满足以上条件在|i-d|，i∈[1,m]中搜索最小值：

当il<d时，下限高程阈值tstart＝lv×il，|il-d|在|i-d|，i∈[1,m]中为最小值；

当ih>d时，上限高程阈值tend＝lv×ih，|ih-d|在|i-d|，i∈[1,m]中为最小值；

步骤1中所述粗去噪后二维点云数据为：

pfilt,

其中pfilt＝(xk,yk)，yk∈(tstart,tend),k∈[1,np]，np为粗去噪后二维点云离散点的数量；

作为优选，步骤2中所述均匀化处理为：

在沿着搭载激光测高仪的卫星或飞机的运行轨道方向，将粗去噪后二维点云数据pfilt按照每lh米为一个间隔，划分成k＝l/lh个沿轨距离片，距离的范围为：

l＝max(xk)-min(xk)，k∈[1,np]

其中，max(xk)为xk，k∈[1,np]的最大值，min(xk)为xk，k∈[1,np]的最小值；

每个沿轨距离片内粗去噪后二维点云离散点为：

pi＝(xu,v,yu,v)，u∈[1,k]，v∈[1,hu]

lh×(u-1)＜xu,v＜lh×u

其中，xu,v为第u个沿轨距离片内第v个粗去噪后二维点云离散点的沿轨距离值，yu,v为第u个沿轨距离片内第v个粗去噪后二维点云离散点的高程，hu为第u个沿轨距离片内粗去噪后二维点云离散点的数量；

以沿轨距离片u为横坐标，第u个沿轨距离片内粗去噪后二维点云离散点的数量hu为纵坐标建立高程统计分布直方图，然后以沿轨距离统计分布直方图为基础，计算粗去噪后二维点云离散点落在每一个沿轨距离片的概率；

其中，pdf(u)表示粗去噪后二维点云离散点落在第u个沿轨距离片的概率，其中u∈{1,2,…，k}；

然后对每个点云离散点的沿轨距离xu,v，进行调整得到均匀后的沿轨距离xu,v,2：

xu,v,2＝k×lhcdf(u)+k×xu,v,rempdf(u+1)

其中，xu,v,2代表均匀化后第u个沿轨距离片中第v个离散点的沿轨距离值，xu,v,rem表示xu,v除以lh的余数，通过点云均匀化，用xu,v,2的值替换xu,v可以得到步骤2中均匀化后二维点云数据paver＝(xk,2,yk),k∈[1,np]；

作为优选，步骤3中所述精去噪方法为：

根据邻域半径eps、核心点阈值minpts以及均匀化后二维点云数据paver通过dbscan聚类算法挑选均匀化后二维点云数据paver中的信号点即精去噪后二维点云数据psignal；

核心点阈值minpts为：

其中，nt为步骤1中所述二维点云离散点的数量，m为步骤1中所述高程片的数量，m1为步骤1中所述高程片内二维点云离散点的数量ei大于nt/m的高程片的数量，n1为m1个高程片中所有二维点云离散点的数量，m2为步骤1中所述高程片内二维点云离散点的数量ei小于等于nt/m的高程片的数量，n2为m2个高程片中所有二维点云离散点的数量；

根据邻域半径eps，核心点阈值minpt和均匀化后二维点云数据paver通过dbscan算法，即可得到精去噪后二维点云数据psignal。

本发明具有以下优点和增益效果：

本发明通过粗去噪预处理了单光子点云数据，粗去噪通过点云数据的垂直距离统计分布直方图来粗略剔除较明显的噪声点，删除了冗余数据，提高了后续算法的效率。

本发明的点云均匀化方法避免了因局部点云分布不均，噪声密度大于或接近于信号密度的现象。实验表明均匀化后的数据，在使用相同的dbscan聚类算法时，提高了查全率和查准率的调和平均值。

本发明的自适应dbscan算法相比于原有的dbscan算法在选取核心点时更具优势，结合点云均匀化方法，提高了f值平均。

附图说明

图1：本发明方法流程图；

图2：点云粗去噪原理；

图3：点云均匀化前后对比；

图4：核心点阈值minpts计算流程；

图5：原始的模拟点云数据；

图6：局部信息统计效果；

图7：原始dbscan聚类效果；

图8：改进的基于聚类的单光子点云数据去噪方法效果；

图9：经过和不经过点云均匀化处理的点云数据去噪方法的f值对比；

图10：三种去噪方法结果对比；

图11：三种去噪方法中目标点的平均高程值与真实目标点平均高程值的差值对比。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。

本发明根据美国icesat-2卫星搭载的光子计数激光测高仪作为本发明实施方式的光子计数激光测高仪，以icesat-2模拟三维点云数据为实施方式数据源，进行本发明的实施方式。

下面结合图1至图11介绍本发明的实施方式，具体步骤如下：

步骤1：将光子计数激光测高仪的三维点云数据转换为二维点云数据，二维点云数据通过粗去噪方法得到粗去噪后二维点云数据；

步骤1中所述三维点云数据转换为二维点云数据为：

将光子计数激光测高仪的三维点云数据p＝(mk,nk,yk)用二维点云数据表示：

p＝(xk,yk)，k∈[1,nt]

其中，nt为三维点云离散点的数量即二维点云离散点的数量，mk表示第k个三维点云离散点的经度，nk表示第k个三维点云离散点的纬度，xk表示第k个二维点云离散点的水平距离值即第k个二维点云离散点的沿轨距离值，且yk表示第k个二维点云离散点的垂直距离值即第k个二维点云离散点的高程；

步骤1中所述粗去噪方法为：

将二维点云数据按照每lv＝1米为一个间隔，划分成m＝h/lv个高程片，m为高程片的数量，高程的范围为：

h＝max(yk)-min(yk)，k∈[1,nt]

其中，nt为二维点云离散点的数量，max(yk)为yk，k∈[1,nt]的最大值，max(yk)为yk，k∈[1,nt]的最小值；

每个高程片内二维点云离散点为：

pi＝(xi,j,yi,j)，i∈[1,m]，j∈[1,ei]

lv×(i-1)＜yi,j＜lv×i

其中，xi,j为第i个高程片内第j个二维点云离散点沿轨距离值，yi,j为第i个高程片内第j个二维点云离散点高程，ei为第i个高程片内二维点云离散点的数量；

以高程片i为横坐标，第i个高程片内二维点云离散点的数量ei为纵坐标建立高程统计分布直方图，并用巴特沃斯滤波器对高程统计分布直方图进行五点加权滤波：

ei,filt＝a1ei-2+a2ei-1+a3ei+a4ei+1+a5ei+2i∈[3,m]

a＝(a1,a2,a3,a4,a5)＝(0.0625,0.25,0.375,0.25,0.0625)

其中，ei,filt为滤波后第i个高程片内二维点云离散点的数量，a1为滤波器第一加权系数，a2为滤波器第二加权系数，a3为滤波器第三加权系数，a4为滤波器第四加权系数，a5为滤波器第五加权系数；

根据滤波后高程片内二维点云离散点的数量得到高程阈值：

ed,filt＝max{e1,filt,e2,filt,...,em,filt}

et＝ed,filt/2

其中，max为求最大值，et为高程阈值；

当高程片(i，ei,filt)i∈[1,m]满足以下条件：

ei<et,i∈[3,m]；

(i，ei,filt)为非极小值点，即ei,filt至少要大于ei-1,filt，ei+1,filt中的任意一个；

lv×|i-d|>＝λ＝80m，d为ed,filt对应的高程片；

满足以上条件在|i-d|，i∈[1,m]中搜索最小值：

当il<d时，下限高程阈值tstart＝lv×il，|il-d|在|i-d|，i∈[1,m]中为最小值；

当ih>d时，上限高程阈值tend＝lv×ih，|ih-d|在|i-d|，i∈[1,m]中为最小值；

步骤1中所述粗去噪后二维点云数据为：

pfilt,

其中pfilt＝(xk,yk)，yk∈(tstart,tend),k∈[1,np]，np为粗去噪后二维点云离散点的数量；

步骤2：将粗去噪后二维点云数据通过均匀化处理得到均匀化后二维点云数据；

步骤2中所述均匀化处理为：

在沿着搭载激光测高仪的卫星或飞机的运行轨道方向，将粗去噪后二维点云数据pfilt按照每lh＝200米为一个间隔，划分成k＝l/lh个沿轨距离片，距离的范围为：

l＝max(xk)-min(xk)，k∈[1,np]

其中，max(xk)为xk，k∈[1,np]的最大值，min(xk)为xk，k∈[1,np]的最小值；

每个沿轨距离片内粗去噪后二维点云离散点为：

pi＝(xu,v,yu,v)，u∈[1,k]，v∈[1,hu]

lh×(u-1)＜xu,v＜lh×u

其中，xu,v为第u个沿轨距离片内第v个粗去噪后二维点云离散点的沿轨距离值，yu,v为第u个沿轨距离片内第v个粗去噪后二维点云离散点的高程，hu为第u个沿轨距离片内粗去噪后二维点云离散点的数量；

以沿轨距离片u为横坐标，第u个沿轨距离片内粗去噪后二维点云离散点的数量hu为纵坐标建立高程统计分布直方图，然后以沿轨距离统计分布直方图为基础，计算粗去噪后二维点云离散点落在每一个沿轨距离片的概率；

其中，pdf(u)表示粗去噪后二维点云离散点落在第u个沿轨距离片的概率，其中u∈{1,2,…，k}；

然后对每个点云离散点的沿轨距离xu,v，进行调整得到均匀后的沿轨距离xu,v,2：

xu,v,2＝k×lhcdf(u)+k×xu,v,rempdf(u+1)

其中，xu,v,2代表均匀化后第u个沿轨距离片中第v个离散点的沿轨距离值，xu,v,rem表示xu,v除以lh的余数，通过点云均匀化，用xu,v,2的值替换xu,v可以得到步骤2中均匀化后二维点云数据paver＝(xk,2,yk),k∈[1,np]；

步骤3：将均匀化后二维点云数据通过精去噪方法得到精去噪后二维点云数据；

步骤3中所述精去噪方法为：

根据邻域半径eps、核心点阈值minpts以及均匀化后二维点云数据paver通过dbscan聚类算法挑选均匀化后二维点云数据paver中的信号点即精去噪后二维点云数据psignal；

邻域半径为eps＝9m；

核心点阈值minpts为：

其中，nt为步骤1中所述二维点云离散点的数量，m为步骤1中所述高程片的数量，m1为步骤1中所述高程片内二维点云离散点的数量ei大于nt/m的高程片的数量，n1为m1个高程片中所有二维点云离散点的数量，m2为步骤1中所述高程片内二维点云离散点的数量ei小于等于nt/m的高程片的数量，n2为m2个高程片中所有二维点云离散点的数量；

根据邻域半径eps，核心点阈值minpt和均匀化后二维点云数据paver通过dbscan算法，即可得到精去噪后二维点云数据psignal。

由图2的原理框图可知，粗去噪统计了二维点云离散点事件在各个高程中出现的概率，并将二维点云离散点事件发生概率较高的数百米高程作为信号所在的高程提取出来，滤掉了大量的噪声，为后续的信号处理算法提供了便利。

由图3可知，由于地形目标反射率和坡度等参数的差异性，二维点云数据在有的地方稀疏，有的地方密集，会出现部分噪声点的密度大于或接近目标点的密度的情况，影响目标点的提取，提取效果的影响如图9所示。二维点云数据均匀化效果如图3所示。

核心点阈值minpts的选取规则如图4的流程所示。icesat-2模拟的三维点云数据，使用局部信息统计方法的效果，使用原始dbscan聚类方法的效果和使用本发明方法的效果分别如图5、图6、图7和图8所示。

由图10可知，本发明方法对点云去噪和信号的提取有着较大的提升。相比于原有的dbscan算法，本发明在所有环境中都能够获取更高的f值。在低噪声情况下，本发明和基于局部统计信息的方法效果相差不大，但是在高噪声情况下本发明能更加有效的提取目标点。并且本发明的计算复杂度小于基于局部统计的方法，能够更加快速的处理点云数据。

由图11可知本文发明获取的平均高程差值波动较小。由统计理论可以得到本发明、基于局部统计方法和基于dbscan方法获取的平均高程差值的均方误差分别为0.052cm、0.077cm和0.153cm，由此可见本发明能更加有效的提取目标点，减少所获平均高程值的随机误差。

本文中所描述的具体实施案例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施案例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。