CN112862720A - 机载LiDAR点云中玻璃板漫反射噪声的去噪方法和系统 - Google Patents

Abstract

本发明公开了一种机载LiDAR点云中玻璃板漫反射噪声的去噪方法，包括：1、基于所有点云S_all估计点云中噪声点的高程范围I_noise并构建候选噪声点集合S_cand；S2、对候选噪声点集合S_cand内的点进行时间约束的三维欧几里得聚类，构成聚类集合S_allcluster；S3、判断聚类集合S_allcluster中每个聚类是否为噪声点聚类，并将噪声点聚类从点云S_all中去除；S4、对去除噪声点聚类后的点云采用点状噪声去噪法去噪。该方法可以将LiDAR在高海拔工作时密封玻璃板漫反射造成的大片、连续分布，甚至可能与地面存在交叉的噪声点与非噪声点区分开，能够最大限度地去除噪声点且保留非噪声点。

Description

机载LiDAR点云中玻璃板漫反射噪声的去噪方法和系统

技术领域

本发明属于机载激光雷达点云数据处理技术领域，具体涉及一种对在大飞行高度采集的LiDAR点云中，由玻璃板漫反射造成的噪声点进行去噪的方法。

背景技术

机载激光雷达(Light Detection And Ranging，LiDAR)已经成为陆地地形测绘的主要手段之一。通过在飞机平台上以极高频率发射激光脉冲并接收后向散射回波，配合高精度POS系统，可以快速获得大量散射面的三维坐标，得到的成果称为点云。

在对高原地区(例如青藏高原平均地形高程超过4000m)进行机载LiDAR测绘时，飞机的绝对飞行高度本身就比较高。如果设计的激光脉冲点密度较低，则航线与地面的相对高度可达4000m甚至5000m以上，飞机的绝对飞行高度可达7000m甚至8000m以上。有的LiDAR硬件可以在这么高的高度工作。但在这种高度空气稀薄、温度和气压偏低。为此，需要对机舱做密封处理，以防止舱内空气外泄，保证舱内气压温度满足机上人员生理需求。对LiDAR进行密封处理是机舱密封处理的重要一环。

在对LiDAR进行密封处理时，需要在激光发射器/接收器外安装透明平滑均匀的玻璃板用于密封。在LiDAR工作时，在玻璃板处，如果不考虑玻璃板对激光脉冲的吸收和散射，则激光脉冲可能发生透射、镜面反射和漫反射，如图1所示。对于玻璃板，要求其透射度高，即穿透玻璃板的激光脉冲功率占发射脉冲功率的比例越高越好，与此相对，要求镜面反射和漫反射的功率越低越好。要实现该目的，一来可以在玻璃板上涂抹防反射涂料，形成涂层，增强玻璃板对激光脉冲波长的透射能力；二来可以调整玻璃板的角度，避免激光脉冲以小入射角入射造成强烈的镜面反射或漫反射。但是，这两种做法都有问题。一般而言，某种涂料只能增强玻璃板对特定波长的透射能力，而不同LiDAR使用的波长不完全一样，这使得涂料的准备存在麻烦。涂料必须涂抹的非常均匀。如果不均匀，会导致透射出的激光脉冲功率不一致，影响回波信号的检测。此外，涂料的保存要求高。这使得目前在玻璃板上涂抹涂料的做法并非主流。调整玻璃板角度的方法不能完全消除玻璃板处的漫散射。玻璃板再平，相对于激光脉冲波长(微米级或亚微米级)来说，其粗糙度(普通玻璃以微米级居多)也不能完全忽略。玻璃板如果沾染灰尘，也可能导致漫反射增强。综上所述，如果使用了玻璃板，则很有可能在玻璃板处存在一定程度的漫反射，进一步地，如果入射激光脉冲在局部表面的入射角为0，则漫反射的激光脉冲也会进入激光接收器。

在激光接收器接收到回波后，需要判断接收到的是信号还是噪声。很明显，玻璃板漫反射的激光脉冲应被判断为噪声。但是，目前机载LiDAR广泛应用多脉冲技术，一个发射脉冲的回波信号还没被接收器接收到，后续的激光脉冲已被发射器发射出来，即同时在空中的激光脉冲有多个。在使用多脉冲技术工作时，通过玻璃板漫反射进入接收器的脉冲，可能与前期发射脉冲的回波在时间范围上发生重叠，即二者被激光接收器接收到的时间非常接近。如果二者的功率均超过背景噪声门限值，则LiDAR可能将漫反射造成的激光脉冲判断为正常信号，由此产生了噪声。

图2-(a)为机载LiDAR常见的点状噪声示意图(白色点为地面点，黑色点为点状噪声)；图2-(b)为玻璃板漫反射造成的噪声点分布示意图(白色点为地面点，黑色点为漫反射噪声点)。和常见的点状噪声相比，这些漫反射造成的噪声点不仅数量众多(可达几万甚至几十万)，而且分布连续，即一个噪声点和邻近噪声点的距离较小(可达米级)，构成噪点簇。有的噪点簇孤悬在空中，如图2-(b)中噪点簇C；有的噪点簇孤悬在地下，如噪点簇A；有的噪点簇与地面点有交叉，如噪点簇B。由于漫反射脉冲的功率较低，所以对应噪点的回波强度一般也较低。如果不对这种噪声点进行去噪，则很可能对后期提取地面点、构建数字高程模型(Digital Elevation Model，DEM)和数字表面模型(Digital Surface Model，DSM)等造成非常不利的影响。

这类噪声点固然可以通过人工手动去噪，但是效率低，容易存在误分、漏分。在自动点云去噪上，这类噪声点并不符合传统点状噪声去噪方法的假设，因此不能采用点状噪声去噪方法。采用限制高程范围的方法去噪并不完全可行，因为有的噪声点和非噪声点在高程分布范围上有重叠。对于噪点簇孤悬分布的情况，有可能基于地形连续分布假设进行去噪，但是对于噪点簇和地面有交叉的情况，无论是传统的点状噪声去噪法、基于曲率的去噪法、基于曲面拟合的去噪法，还是基于空间自相关的去噪方法均不能取得良好的效果。

发明内容

发明目的：针对现有技术中存在的问题，本发明提供一种机载LiDAR点云去噪方法，该方法可以将LiDAR在高海拔工作时密封玻璃板漫反射造成的大片、连续分布，甚至可能与地面存在交叉的噪声点与非噪声点区分开，最大限度地去除噪声点且保留非噪声点。

技术方案：本发明一方面公开了一种机载LiDAR点云中密封玻璃板漫反射噪声的去噪方法，包括：

S1、估计点云中噪声点的高程范围I_noise并构建候选噪声点集合S_cand，具体包括：

S11、在投影坐标系下将点云集合S_all内的点垂直投影到XY平面，计算S_all内点的X坐标最大值x_max和最小值x_min，Y坐标最大值y_max和最小值y_min，并基于XY坐标构建点云二维空间索引；

S12、对点云在XY平面的投影区域进行栅格化，每个栅格为边长是d_cell的正方形，

D为平均点密度；设X方向的栅格数为N_col，Y方向的栅格数为N_row；

为向上取整运算符；

每个栅格的中心点为二维噪声检测点，第i行第j列的二维噪声检测点cell_ij(x_dtc,y_dtc)坐标的计算方法是：

S13、构建可能噪点集合S_pn；利用点云二维空间索引在S_all中搜索以二维噪声检测点cell_ij为圆心，d_cell/2为半径的圆形邻域内的点，设搜索到的点构成集合s_ij，计算s_ij内点的最小高程z_{min ij}和最大高程z_{max ij}，得到高程差△z_ij：△z_ij＝z_{max ij}-z_{min ij}；

如果△z_ij大于等于高程差阈值△z_thre，则将s_ij内的点放入可能噪点集合S_pn；

将可能噪点集合S_pn中点的高程范围[z_{min pn},z_{max pn}]以粗粒度均匀划分为多个粗子区间，统计S_pn落入每个粗子区间的点数，取点数最多的L个粗子区间的并集作为噪声点的第一高程范围I₁；

将第一高程范围I₁以细粒度均匀划分为多个细子区间，统计S_pn落入每个细子区间的点数，取点数最多的K个细子区间的并集作为噪声点的高程范围I_noise；

S14、构建候选噪声点集合S_cand，所述候选噪声点集合S_cand由点云集合S_all中同时满足如下条件的点组成：

条件一：高程位于I_noise范围内；

条件二：回波强度位于预设的强度范围[it_min,it_max]内；

条件三：回波序号属于预设的回波序号S_echo；

对候选噪声点集合S_cand内的点建立点云三维空间索引；

S2、对候选噪声点集合S_cand内的点进行时间约束的三维欧几里得聚类，构建聚类集合S_allcluster，具体包括：

S21、为候选噪声点集合S_cand中的每个点添加状态标签，并初始化为未使用；建立聚类集合S_allcluster；

S22、遍历S_cand中的点，对于其中状态标签为未使用的点p₀构建聚类S_cluster，具体步骤包括：

S221、将p₀的状态标签设置为已使用，建立聚类S_cluster、种子点集合S_seed、扩张点集合S_new，将p₀加入S_cluster和S_seed；

S222、清空扩张点集合S_new，遍历S_seed内的点，对于其中点p_current，建立邻域点集合S_neighb，使用点云三维空间索引将S_cand中以p_current为中心，R为半径的球形邻域内除p_current外的点加入S_neighb中；

将S_neighb中状态标签为未使用，且与p_current的时间差小于时间差阈值△t_thre的点的状态标签设置为已使用，加入扩张点集合S_new中；

S223、在遍历完S_seed内的点后，如果S_new不为空，将S_new的点加入S_cluster，更新种子点集合S_seed为S_new；跳转至步骤S222进行新的扩张；

如果S_new为空，判断S_cluster内的点数是否位于[N_min,N_max]范围内，如是，则S_cluster为有效聚类，加入聚类集合S_allcluster中；N_min和N_max分别为预设的聚类最小点数阈值、聚类最大点数阈值；

跳转至步骤S221，为S_cand中下一个状态标签为未使用的点构建聚类；

S3、判断聚类集合S_allcluster中每个聚类是否为噪声点聚类，并将噪声点聚类从点云S_all中去除，具体包括：

S31、利用非聚类点集合S_nocluster内点的XY坐标重新构建点云二维索引，所述非聚类点集合S_nocluster由S_all去除聚类集合S_allcluster中的所有点得到，即：

S_nocluster＝S_all-S_allcluster；

S32、遍历S_allcluster内的每个聚类，判断每个聚类是否为噪声，包括：

S321、设当前聚类为S′_cluster，S′_cluster∈S_allcluster，将S′_cluster所有点垂直投影到XY平面，计算其二维边界，寻找位于二维边界上的点，设S′_cluster位于二维边界上的点构成集合S_bound，共N_bound个点；

S322、对S_bound内的点进行遍历，设其中第u个点为

利用点云二维空间索引在S_nocluster中搜索以

为圆心，d_cell/2为半径的圆形邻域内的邻域点，设这些邻域点构成集合

计算

内点的最大高程

和最小高程

以及

和

内点的最小高程差

其中min(·)为计算最小值函数；

S323、对S_bound内的所有N_bound个点，计算

的平均值

S324、如果

大于等于高程差阈值△z_thre，则认为S′_cluster内的点均为噪声，将S′_cluster内的点从点云S_all中去除；跳转至步骤S321处理S_allcluster内的下一个聚类；

S4、对去除噪声点聚类后的点云采用点状噪声去噪法进行去噪。

所述步骤S11中采用KD树或四叉树建立点云二维空间索引。

所述高程差阈值△z_thre的取值为：

如果测区60％以上面积为山区，则△z_thre需要满足

θ为地形坡度；

如果测区含有树木、建筑物，且树高、建筑物高的最大值约为z_{max hgt}，则△z_thre需要满足△z_thre>z_{max hgt}。

作为优选，所述步骤S13中还包括：

对可能噪点集合S_pn中点的高程范围[z_{min pn},z_{max pn}]进行扩展，得到扩展高程范围[z_{min pn}-α,z_{max pn}+β]，α,β分别为下、上界扩展常数，α≥0，β≥0，且α,β不同时为0；

对扩展高程范围[z_{min pn}-α,z_{max pn}+β]以粗粒度均匀划分为多个粗子区间，统计S_pn落入每个粗子区间的点数，取点数最多的L个粗子区间的并集作为噪声点的第二高程范围I₂；

计算第一高程范围I₁和第二高程范围I₂的并集，作为第三高程范围I₃；

将第三高程范围I₃以细粒度均匀划分为多个细子区间，统计S_pn落入每个细子区间的点数，取点数最多的K个细子区间的并集作为噪声点的高程范围I_noise。

作为优选，所述步骤S13还包括：

对高程范围I_noise的每个细子区间向两侧扩展，扩展后的细子区间并集为最终的噪声点的高程范围I_noise。

所述步骤S32中计算点云聚类二维边界的方法包括Alpha-shape法、栅格到矢量转换法、基于不规则三角网的边界追踪法。

作为优选，所述步骤S3之后还包括修改所述步骤S14中预设的回波序号S_echo，重复执行步骤S2和S3，多次去除噪声点聚类。

所述步骤S14中预设的回波序号S_echo的取值包括：唯一回波、多次回波的第一次回波、多次回波的非第一次回波。

所述点状噪声去噪法包括：统计离群点移除、半径离群点移除。

另一方面，本发明公开了实现上述去噪方法的机载LiDAR点云中玻璃板漫反射噪声的去噪系统，包括：

候选噪声点集合构建模块1，用于根据权利要求1中步骤S1估计点云中噪声点的高程范围I_noise并构建候选噪声点集合S_cand；

候选噪声点集合聚类模块2，用于根据权利要求1中步骤S2对候选噪声点集合S_cand内的点进行时间约束的三维欧几里得聚类，构建聚类集合S_allcluster；

去除噪声点聚类模块3，用于根据权利要求1中步骤S3判断聚类集合S_allcluster中每个聚类是否为噪声点聚类，并将噪声点聚类从点云S_all中去除；

点状噪声去噪模块4，用于对去除噪声点聚类后的点云采用点状噪声去噪法进行去噪。

有益效果：本发明公开的针对机载LiDAR点云中由于密封玻璃板漫反射造成的大片、连续分布的噪声点进行去噪的方法，先基于同一平面位置的高程差筛查可能存在噪声点的高程范围，再将点云按照时空属性进行分割，得到聚类，并基于地形连续分布假设判断每个聚类是否属于噪声。相比传统的仅利用坐标和分类信息进行去噪的方法，本发明提出的方法对点云分割的更为完全，可以在很大程度上将悬浮的噪点簇和穿过地面的噪点簇分割出来并判断为噪声点。

附图说明

图1为玻璃板处激光脉冲传输的示意图；

图2为噪声点垂直分布示意图；

图3为本发明公开的去噪方法流程图；

图4为本发明公开的去噪系统组成示意图。

具体实施方式

下面结合附图和具体实施方式，进一步阐明本发明。

对于机载LiDAR来说，原始点云数据一般经过坐标变换到投影坐标系下进行处理。投影坐标系是由大地基准面和地图投影方法共同确定的。在中国，机载LiDAR点云常用的高斯-克吕格投影、通用横轴墨卡托投影等，以中央经线和赤道的交点作为坐标原点，以中央经线的投影为纵坐标X轴，以赤道的投影为横坐标Y轴。点云高程(通常表达为Z坐标值)为点相对于大地基准面的高度。本实施例是在投影坐标系下对点云数据中的玻璃板漫反射噪声进行处理。

如图3所示，本发明公开了一种机载LiDAR点云中玻璃板漫反射噪声的去噪方法，包括：

S1、估计点云中噪声点的高程范围I_noise并构建候选噪声点集合S_cand，具体包括：

S11、在投影坐标系下将点云集合S_all内的点垂直投影到XY平面，计算S_all内点的X坐标最大值x_max和最小值x_min，Y坐标最大值y_max和最小值y_min，并基于XY坐标构建点云二维空间索引；此处可以采用KD树或四叉树建立点云二维空间索引；

S12、对点云在XY平面的投影区域进行栅格化，每个栅格为边长是d_cell的正方形，

D为平均点密度，单位一般为m^-2；设X方向的栅格数为N_col，Y方向的栅格数为N_row；

为向上取整运算符；

每个栅格的中心点为二维噪声检测点，第i行第j列的二维噪声检测点cell_ij(x_dtc,y_dtc)坐标的计算方法是：

S13、构建可能噪点集合S_pn；利用点云二维空间索引在S_all中搜索以二维噪声检测点cell_ij为圆心，d_cell/2为半径的圆形邻域内的点，设搜索到的点构成集合s_ij，计算s_ij内点的最小高程z_{min ij}和最大高程z_{max ij}，得到高程差△z_ij：△z_ij＝z_{max ij}-z_{min ij}；

如果△z_ij大于等于高程差阈值△z_thre，则将s_ij内的点放入可能噪点集合S_pn；

高程差阈值△z_thre的取值与测区具体的地形、地类有关；如果测区地形以山区居多，则△z_thre需要满足

θ为地形坡度；本实施例中设置当测区60％以上面积为山区时，认为测区地形以山区居多；如果测区含有树木、建筑物，且树高、建筑物高的最大值约为z_{max hgt}，则△z_thre需要满足△z_thre>z_{max hgt}。

将可能噪点集合S_pn中点的高程范围[z_{min pn},z_{max pn}]以粗粒度均匀划分为多个粗子区间，统计S_pn落入每个粗子区间的点数，取点数最多的L个粗子区间的并集作为噪声点的第一高程范围I₁；

本实施例中采用高程分布直方图来计算第一高程范围I₁，具体为：

计算S_pn内点的最小高程和最大高程z_{min pn},z_{max pn}，设直方图统计用的最小高程值z_low＝z_{min pn}；

对S_pn内所有点以较大组间距△z_L计算高程分布直方图，高程为z_pn的点所在的组号g为

其中

为向下取整符号，g的取值范围为

统计每组落入的S_pn内点的数量，则第g组的高程区间I_g,pn为I_g,pn＝[z_low+g△z_L,z_low+(g+1)△z_L]，称得到的直方图为第一直方图；

△z_L的取值应略微大于噪声点簇的高程差，可以通过点云处理软件目视识别出点云中的噪声点簇，获得单个噪点簇的最小高程z_noise,low和最大高程z_noise,high，计算高程差△z_noise：△z_noise＝z_noise,high-z_noise,low，对所有噪点簇，计算高程差△z_noise的最大值△z_noise,max，需要确保△z_L>△z_noise,max；

设漫反射噪声点高程主要集中在L个不重叠的区间(在边界处也无重叠)内，L≥1，在第一直方图中选择点数最多的前L个组，设对应的高程区间分别为I_1,1,I_1,2,…,I_1,L，则第一高程范围I₁为这些高程区间的并集：

I₁＝I_1,1∪I_1,2∪…∪I_1,L

L的取值可以通过使用点云处理软件查看点云获得，一般来说，对于小范围的点云，L不大于1；

为了避免由于区间端点设置导致同一噪点簇被分到多个组，导致每个组的点数均偏少，本实施例对可能噪点集合S_pn中点的高程范围[z_{min pn},z_{max pn}]进行扩展，对扩展后的高程范围进行第二次粗粒度划分，并获取第二高程范围，具体为：

对可能噪点集合S_pn中点的高程范围[z_{min pn},z_{max pn}]进行扩展，得到扩展高程范围[z_{min pn}-α,z_{max pn}+β]，α,β为下、上界扩展常数，α≥0，β≥0，且α,β不同时为0；本实施例中，令α＝△z_L/2，β＝0，即向下界扩展△z_L/2；

对扩展高程范围[z_{min pn}-α,z_{max pn}+β]以较大组间距△z_L进行粗粒度均匀划分，得到多个粗子区间，统计S_pn落入每个粗子区间的点数，取点数最多的L个粗子区间的并集作为噪声点的第二高程范围I₂；

计算第一高程范围I₁和第二高程范围I₂的并集，得到第三高程范围I₃：I₃＝I₁∪I₂；这样尽量不遗漏噪声点可能存在的高程范围。

为了尽量精确地确定噪点的高程范围，将第三高程范围I₃以细粒度均匀划分为多个细子区间，统计S_pn落入每个细子区间的点数，计算噪声点高程范围I_noise，具体步骤为：

以较小组间距△z_S将S_pn内高程位于I₃内的点依据高程再次计算高程分布直方图，设得到第三直方图；

△z_S的取值必须小于等于△z_L/2，并大于同一噪点簇不重复高程值中的相邻高程值的最大差△z_{max dif}，△z_{max dif}的计算方法是：对同一噪点簇内的点的高程值去重，从小到大排列，设为z₀,z₁,…,z_K，遍历z₀,z₁,…,z_K，设当前高程值为z_i，计算△z＝z_i-z_i-1，在各个△z中寻找最大值△z_{max dif}；一般，△z_S取1m或0.5m即可满足需求；

设第三直方图中点数大于0的组共构成K个不重叠的子区间(在边界处也无重叠)，即I_3,1,I_3,2,…,I_3,K。如果K≤L，此时I_noise＝I_3,1∪I_3,2∪…∪I_3,K；如果K>L，则从I_3,1,I_3,2,…,I_3,K中挑选落入点数最多的前L个区间I_4,1,I_4,2,…,I_4,L，此时I_noise＝I_4,1∪I_4,2∪…∪I_4,L。

为了尽可能全面地检测噪声点所在高程范围，避免前述步骤使用较大的组间距△z_L造成的噪声点高程范围不全，对高程范围I_noise的每个细子区间向两侧扩展，即对一个子区间[z_low,z_high]，向两侧扩展ξ的结果是变成[z_low-ξ,z_high+ξ]。扩展后的细子区间并集为最终的噪声点的高程范围I_noise；ξ>0，ξ的取值需要小于△z_S/2。

S14、构建候选噪声点集合S_cand，所述候选噪声点集合S_cand由点云集合S_all中同时满足如下条件的点组成：

条件一：高程位于I_noise范围内；

条件二：回波强度位于预设的强度范围[it_min,it_max]内；

条件三：回波序号属于预设的回波序号S_echo；

it_min可以简单取0，it_max可以使用点云处理软件，根据噪声点实际情况取值。

对候选噪声点集合S_cand内的点建立点云三维空间索引；此处索引可以是KD树或八叉树；

S2、对候选噪声点集合S_cand内的点进行时间约束的三维欧几里得聚类，构建聚类集合S_allcluster，具体包括：

S21、为候选噪声点集合S_cand中的每个点添加状态标签，并初始化为未使用；建立聚类集合S_allcluster；

S22、遍历S_cand中的点，对于其中状态标签为未使用的点p₀构建聚类S_cluster，具体步骤包括：

S221、将p₀的状态标签设置为已使用，建立聚类S_cluster、种子点集合S_seed、扩张点集合S_new，将p₀加入S_cluster和S_seed；

S222、清空扩张点集合S_new，遍历S_seed内的点，对于其中点p_current，建立邻域点集合S_neighb，使用点云三维空间索引将S_cand中以p_current为中心，R为半径的球形邻域内除p_current外的点加入S_neighb中；

R需要根据噪声点的实际情况取值，大于噪点簇中相邻噪点的最大距离；

将S_neighb中状态标签为未使用，且与p_current的时间差△t小于时间差阈值△t_thre的点的状态标签设置为已使用，加入扩张点集合S_new中；

设置时间差约束的原因是，和地面发生交叉的噪点簇和地面点在获取时间上可能有差异，种子点在扩张时只扩张到和自己获取时间接近的点，避免由噪声点扩张到地面点，或由地面点扩张到噪声点；△t_thre的取值与LiDAR的发射频率和点云的具体情况有关，一般可以先取0.02s至0.2s内的值，看是否能取得良好效果，如果效果不好再进行调整；

S223、在遍历完S_seed内的点后，如果S_new不为空，将S_new的点加入S_cluster，更新种子点集合S_seed为S_new；跳转至步骤S222进行新的扩张；

如果S_new为空，判断S_cluster内的点数是否位于[N_min,N_max]范围内，如是，则S_cluster为有效聚类，加入聚类集合S_allcluster中；N_min和N_max分别为预设的聚类最小点数阈值、聚类最大点数阈值，可以根据噪点簇的具体情况设置，其原则是尽量使得隶属于同一个噪点簇的点可以被分割到一个聚类，N_min需要大于点状噪声可能的最大点数，N_max需要小于连续分布的地面点的最小点数；

跳转至步骤S221，为S_cand中下一个状态标签为未使用的点构建聚类；

S3、判断聚类集合S_allcluster中每个聚类是否为噪声点聚类，并将噪声点聚类从点云S_all中去除，具体包括：

S31、对非聚类点集合S_nocluster内的点，依据其XY坐标重新构建点云二维索引，所述非聚类点集合S_nocluster由S_all去除聚类集合S_allcluster中的所有点得到，即：

S_nocluster＝S_all-S_allcluster；

此处的二维索引可以是KD树或四叉树；

S32、遍历S_allcluster内的每个聚类，判断每个聚类是否为噪声，包括：

S321、设当前聚类为S_c′_luster，S_c′_luster∈S_allcluster，将S_c′_luster所有点垂直投影到XY平面，计算其二维边界，寻找位于二维边界上的点，设S_c′_luster位于二维边界上的点构成集合S_bound，共N_bound个点；

计算点云聚类二维边界的方法包括Alpha-shape法、栅格到矢量转换法、基于不规则三角网的边界追踪法等。

S322、对S_bound内的点进行遍历，设其中第u个点为

利用点云二维空间索引在S_nocluster中搜索以

为圆心，d_cell/2为半径的圆形邻域内的邻域点，设这些邻域点构成集合

计算

内点的最大高程

和最小高程

以及

和

内点的最小高程差

其中min(·)为计算最小值函数；

S323、对S_bound内的所有N_bound个点，计算

的平均值

S324、如果

大于等于高程差阈值△z_thre，则认为S_c′_luster内的点均为噪声，将S_c′_luster内的点从点云S_all中去除；跳转至步骤S321处理S_allcluster内的下一个聚类；

该步骤基于地形分布连续性假设来判断聚类点是否噪声，认为噪声点聚类导致地形不连续，与聚类周边点间存在较大阶跃，而正常点与聚类周边点间没有较大阶跃。

S4、修改所述步骤S14中预设的回波序号S_echo，重复执行步骤S2和S3，多次去除噪声点聚类；

高于地面较多的噪点簇其回波一般是多次回波的第一次回波，低于地面较多的噪点簇一般是多次回波的非第一次回波(可能是第二次、第三次等)，而与正常地面有交叉的噪点簇在地面附近可能是唯一一次回波。如果在步骤S2对点云进行单次欧几里得聚类时S_echo包括所有回波序号，则有可能导致在聚类结果中地面点和噪声点混在一个聚类内，导致分割不完全。因此，设置S_echo共有3个可能的值，即唯一回波、多次回波的第一次回波、多次回波的非第一次回波。多次执行S2和S3，确保这3个S_echo值至少用到一次。

因为不能确保前面步骤去除所有噪声点，S_all中可能还有点状噪声存在，所以对去除噪声点聚类后的点云采用点状噪声去噪法进行去噪。可以采用的点状噪声去噪法包括：统计离群点移除、半径离群点移除等。

本实施例还公开了实现上述机载LiDAR点云中密封玻璃板漫反射噪声去噪方法的系统，如图4所示，包括：

候选噪声点集合构建模块1，用于根据1步骤S1估计点云中噪声点的高程范围I_noise并构建候选噪声点集合S_cand；

候选噪声点集合聚类模块2，用于根据步骤S2对候选噪声点集合S_cand内的点进行时间约束的三维欧几里得聚类，构建聚类集合S_allcluster；

去除噪声点聚类模块3，用于根据步骤S3判断聚类集合S_allcluster中每个聚类是否为噪声点聚类，并将噪声点聚类从点云S_all中去除；

点状噪声去噪模块4，用于根据步骤S4对去除噪声点聚类后的点云采用点状噪声去噪法去噪。

Claims (10)

1.一种机载LiDAR点云中玻璃板漫反射噪声的去噪方法，其特征在于，包括：

S1、估计点云中噪声点的高程范围I_noise并构建候选噪声点集合S_cand，具体包括：

S11、在投影坐标系下将点云集合S_all内的点垂直投影到XY平面，计算S_all内点的X坐标最大值x_max和最小值x_min，Y坐标最大值y_max和最小值y_min，并基于XY坐标构建点云二维空间索引；

S12、对点云在XY平面的投影区域进行栅格化，每个栅格为边长是d_cell的正方形，

D为平均点密度；设X方向的栅格数为N_col，Y方向的栅格数为N_row；

为向上取整运算符；

每个栅格的中心点为二维噪声检测点，第i行第j列的二维噪声检测点cell_ij(x_dtc,y_dtc)坐标的计算方法是：

S13、构建可能噪点集合S_pn；利用点云二维空间索引在S_all中搜索以二维噪声检测点cell_ij为圆心，d_cell/2为半径的圆形邻域内的点，设搜索到的点构成集合s_ij，计算s_ij内点的最小高程z_minij和最大高程z_maxij，得到高程差△z_ij：△z_ij＝z_maxij-z_minij；

如果△z_ij大于等于高程差阈值△z_thre，则将s_ij内的点放入可能噪点集合S_pn；

将可能噪点集合S_pn中点的高程范围[z_minpn,z_maxpn]以粗粒度均匀划分为多个粗子区间，统计S_pn落入每个粗子区间的点数，取点数最多的L个粗子区间的并集作为噪声点的第一高程范围I₁；

将第一高程范围I₁以细粒度均匀划分为多个细子区间，统计S_pn落入每个细子区间的点数，取点数最多的K个细子区间的并集作为噪声点的高程范围I_noise；

S14、构建候选噪声点集合S_cand，所述候选噪声点集合S_cand由点云集合S_all中同时满足如下条件的点组成：

条件一：高程位于I_noise范围内；

条件二：回波强度位于预设的强度范围[it_min,it_max]内；

条件三：回波序号属于预设的回波序号S_echo；

对候选噪声点集合S_cand内的点建立点云三维空间索引；

S2、对候选噪声点集合S_cand内的点进行时间约束的三维欧几里得聚类，构建聚类集合S_allcluster，具体包括：

S21、为候选噪声点集合S_cand中的每个点添加状态标签，并初始化为未使用；建立聚类集合S_allcluster；

S22、遍历S_cand中的点，对于其中状态标签为未使用的点p₀构建聚类S_cluster，具体步骤包括：

S221、将p₀的状态标签设置为已使用，建立聚类S_cluster、种子点集合S_seed、扩张点集合S_new，将p₀加入S_cluster和S_seed；

S222、清空扩张点集合S_new，遍历S_seed内的点，对于其中点p_current，建立邻域点集合S_neighb，使用点云三维空间索引将S_cand中以p_current为中心，R为半径的球形邻域内除p_current外的点加入S_neighb中；

将S_neighb中状态标签为未使用，且与p_current的时间差小于时间差阈值△t_thre的点的状态标签设置为已使用，加入扩张点集合S_new中；

S223、在遍历完S_seed内的点后，如果S_new不为空，将S_new的点加入S_cluster，更新种子点集合S_seed为S_new；跳转至步骤S222进行新的扩张；

如果S_new为空，判断S_cluster内的点数是否位于[N_min,N_max]范围内，如是，则S_cluster为有效聚类，加入聚类集合S_allcluster中；N_min和N_max分别为预设的聚类最小点数阈值、聚类最大点数阈值；

跳转至步骤S221，为S_cand中下一个状态标签为未使用的点构建聚类；

S3、判断聚类集合S_allcluster中每个聚类是否为噪声点聚类，并将噪声点聚类从点云S_all中去除，具体包括：

S31、对非聚类点集合S_nocluster内的点，依据其XY坐标重新构建点云二维索引，所述非聚类点集合S_nocluster由S_all去除聚类集合S_allcluster中的所有点得到，即：

S_nocluster＝S_all-S_allcluster；

S32、遍历S_allcluster内的每个聚类，判断每个聚类是否为噪声，包括：

S321、设当前聚类为S′_cluster，S′_cluster∈S_allcluster，将S_c′_luster所有点垂直投影到XY平面，计算其二维边界，寻找位于二维边界上的点，设S_c′_luster位于二维边界上的点构成集合S_bound，共N_bound个点；

S322、对S_bound内的点进行遍历，设其中第u个点为

利用点云二维空间索引在S_nocluster中搜索以

为圆心，d_cell2为半径的圆形邻域内的邻域点，设这些邻域点构成集合

计算

内点的最大高程

和最小高程

以及

和

内点的最小高程差

其中min(·)为计算最小值函数；

S323、对S_bound内的所有N_bound个点，计算

的平均值

S324、如果

大于等于高程差阈值△z_thre，则认为S′_cluster内的点均为噪声，将S′_cluster内的点从点云S_all中去除；跳转至步骤S321处理S_allcluster内的下一个聚类；

S4、对去除噪声点聚类后的点云采用点状噪声去噪法去噪。

2.根据权利要求1所述的机载LiDAR点云中玻璃板漫反射噪声的去噪方法，其特征在于，所述步骤S11中采用KD树或四叉树建立点云二维空间索引。

3.根据权利要求1所述的机载LiDAR点云中玻璃板漫反射噪声的去噪方法，其特征在于，所述高程差阈值△z_thre的取值为：

如果测区60％以上面积为山区，则△z_thre需要满足

θ为地形坡度；

如果测区含有树木、建筑物，且树高、建筑物高的最大值约为z_maxhgt，则△z_thre需要满足△z_thre>z_maxhgt。

4.根据权利要求1所述的机载LiDAR点云中玻璃板漫反射噪声的去噪方法，其特征在于，所述步骤S13中还包括：

对可能噪点集合S_pn中点的高程范围[z_minpn,z_maxpn]进行扩展，得到扩展高程范围[z_minpn-α,z_maxpn+β]，α,β为下、上界扩展常数，α≥0，β≥0，且α,β不同时为0；

对扩展高程范围[z_minpn-α,z_maxpn+β]以粗粒度均匀划分为多个粗子区间，统计S_pn落入每个粗子区间的点数，取点数最多的L个粗子区间的并集作为噪声点的第二高程范围I₂；

计算第一高程范围I₁和第二高程范围I₂的并集，作为第三高程范围I₃；

将第三高程范围I₃以细粒度均匀划分为多个细子区间，统计S_pn落入每个细子区间的点数，取点数最多的K个细子区间的并集作为噪声点的高程范围I_noise。

5.根据权利要求1所述的机载LiDAR点云中玻璃板漫反射噪声的去噪方法，其特征在于，所述步骤S13还包括：

对高程范围I_noise的每个细子区间向两侧扩展，扩展后的细子区间并集为最终的噪声点的高程范围I_noise。

6.根据权利要求1所述的机载LiDAR点云中玻璃板漫反射噪声的去噪方法，其特征在于，所述步骤S32中计算点云聚类二维边界的方法包括Alpha-shape法、栅格到矢量转换法、基于不规则三角网的边界追踪法。

7.根据权利要求1所述的机载LiDAR点云中玻璃板漫反射噪声的去噪方法，其特征在于，所述步骤S3之后还包括修改所述步骤S14中预设的回波序号S_echo，重复执行步骤S2和S3，多次去除噪声点聚类。

8.根据权利要求7所述的机载LiDAR点云中玻璃板漫反射噪声的去噪方法，其特征在于，所述步骤S14中预设的回波序号S_echo的取值包括：唯一回波、多次回波的第一次回波、多次回波的非第一次回波。

9.根据权利要求1所述的机载LiDAR点云中玻璃板漫反射噪声的去噪方法，其特征在于，所述点状噪声去噪法包括：统计离群点移除、半径离群点移除。

10.一种机载LiDAR点云中玻璃板漫反射噪声的去噪系统，其特征在于，包括：

候选噪声点集合构建模块(1)，用于根据权利要求1中步骤S1估计点云中噪声点的高程范围I_noise并构建候选噪声点集合S_cand；

候选噪声点集合聚类模块(2)，用于根据权利要求1中步骤S2对候选噪声点集合S_cand内的点进行时间约束的三维欧几里得聚类，构建聚类集合S_allcluster；

去除噪声点聚类模块(3)，用于根据权利要求1中步骤S3判断聚类集合S_allcluster中每个聚类是否为噪声点聚类，并将噪声点聚类从点云S_all中去除；

点状噪声去噪模块(4)，用于对去除噪声点聚类后的点云采用点状噪声去噪法去噪。

Images