一种顾及地形特点的多波束点云数据去噪方法与流程

本发明涉及一种顾及地形特点的多波束点云数据去噪方法。
背景技术：
：船载多波束测深系统能够快速获取海底表面包含三维坐标和回波强度等信息的点云数据，为生成高精度DEM提供保障。在建立DEM之前，必须对点云数据进行去噪。由于海底地形的多样性和复杂性，噪声依附于地表而难以去除，手工处理海量数据消耗大量存储空间和计算时间，甚至有可能“误删”。目前，国内外关于点云数据去噪的研究较多，但对于噪声较多的多波束点云数据研究相对较少，现有技术方案难以去除近地表噪声，且无法在保证去噪效果的同时保证效率。下面给出几种与本发明相近的点云去噪的技术方案：第一种是PCL(PointCloudLibrary，点云库)在版本1.8.0中的统计滤波器，此种实现方案的原理为：计算每个点到近邻域各点的平均距离作为该点的值，然后计算平均距离的均值和标准差，判断该点是否在中误差范围内。该方法能够适应地形特点、执行效率高且去噪效果较优，但由于边缘数据点到k近邻点的平均距离值较大，往往处于标准差范围外，导致边缘信息难以较好保留，而且此种实现方案没有考虑到地形特征。第二种是基于邻域点特征选择的双边滤波，此种实现方案的原理为：首先通过散乱点云格网化建立拓扑关系，检索距离点云最近的k近邻点作为该点的邻域点，计算格网点云的平均欧式距离与某点邻域内平均欧式距离的比值T，若T小于给定阈值，采用k邻域内的点云计算滤波因子，否则采用整个网格内的点云。缺点：滤波因子计算复杂，阈值难以界定。第三种是高斯曲率滤波，此种实现方案的原理为：计算各处点云的曲率值，根据点云的变化情况采取不同的滤波方案，对于尖锐区域效果较好，但点云的曲率计算较为复杂，对于噪声较多的实测多波束点云数据，效率低且无法达到预期效果。第四种是建立格网索引滤波，此种实现方案的原理为：设置格网大小，对散乱点云建立格网索引，若格网内的点云数目小于给定的阈值，则删除格网内所有点云，其他基于格网的基本思想类似，不同的是对删除点做一些判断。缺点：此类方法难以区分近地表噪声数据，且对于阈值的选取要求较高，选择较为宽松则难以较好的去除噪声，设置较为苛刻，容易造成过度去除，且对于稀疏的非噪声点，难以判别。综上，现有技术中已有的点云去噪的技术方案具有如下缺点：1)效果局限：现有方案难以较好的区分近地表的噪声数据，或者无法保证在去除近地面噪声数据的同时较好的保留地形数据；2)适应地形能力差：现有方案分别适用于不同的点云数据，难以适应地形的变化。3)去噪效果与方法效率之间难以达到平衡：效率高的方法难以保证效果，效果较好的方法往往计算复杂，难以保证效率。技术实现要素：本发明的目的在于提出一种顾及地形特点的多波束点云数据去噪方法，该方法从顾及地形的角度出发，去除多波束点云数据中的近地表噪声和明显离群噪声数据，同时较好的保留边缘等信息，在保证以上效果的基础上优化设计方案，提高执行效率。为了实现上述目的，本发明采用如下技术方案：一种顾及地形特点的多波束点云数据去噪方法，包括如下步骤：1)输入点云集cloudPoint，建立KD索引树，遍历点云集cloudPoint中的每一个点searchPoint，执行步骤2)-步骤7)；遍历结束后，执行步骤8)；2)检查searchPoint是否是无效点，检查searchPoint是否已经被搜索过；若searchPoint是无效点或者已经被搜索过，则跳出本次循环，转到步骤1)；否则，转到步骤3)；3)设置搜索半径r，对searchPoint进行k近邻域半径搜索，搜索到点数记为K_Number，将搜索到的点和对应的索引号保存至向量容器pointVector中；4)保存向量容器pointVector中有效数据点的索引号和标记值0至哈希表1中，用于判断当前点searchPoint是否被搜索；5)设置阈值λ1，若K_Number小于λ1，则认定向量容器pointVector中的数据点为噪声数据点，标记为无效点并存储该点的索引号和距离值0至哈希表2中，执行步骤1)；否则，执行步骤6)；6)设置RANSAC平面拟合算法的阈值λ2，对向量容器pointVector中的点基于RANSAC算法拟合平面，平面记为A；7)计算向量容器pointVector中索引号对应的点到平面A的距离，并保存索引号和对应距离至哈希表2表中，转到步骤1)；8)计算哈希表2中不为0的距离值的平均值和中误差，分别记为u和σ；9)分别计算相邻平面法向量夹角|α|和投影距离|dij|，若同时不满足共向性和共面性条件，则认定当前拟合平面为离群面并标记相应的点为噪声数据点，若相邻平面不满足共向性但满足共面性，则认定当前拟合平面处地形为起伏较大的陡坡，标记相应点为非噪声数据点；10)设置阈值λ3，若哈希表2中值不为0的索引号对应的距离在[u-λ3·σ,μ+λ3·σ]之间，则认定该索引号对应cloudPoint中的数据点为非噪声数据点，否则认定为噪声数据点；11)分别保存非噪声数据点和噪声数据点数据至PCD格式文件中。优选地，所述步骤3)中，k近邻域半径搜索的过程为：对于KD索引树上任一点，搜索以该点为球心，以r为半径区域内的点云，若搜索到点的数目小于k，则认定该区域内所有的点作为该点的k近邻域数据点；否则，选择距离该点最近的k个点作为该点的k近邻域数据点。优选地，所述步骤3)中k的取值范围为：500<k≤1000。优选地，所述步骤3)中r值的计算公式为：S·k＝π·(r/3)2；其中，S表示各个数据点占据区域面积；根据设备的相关参数和实际外业情况分别计算采集点云数据的横向间距和纵向间距，则S＝横向间距×纵向间距。优选地，所述步骤5)中λ1值与k值大小关系为：λ1＝0.2×k。优选地，所述步骤6)中λ2值的确定步骤为：根据设备的相关参数和实际外业情况分别计算采集点云数据的横向间距和纵向间距，λ2取横向间距和纵向间距之中的最小值。优选地，所述步骤6)中平面A的拟合平面方程形式为：ax+by+cz+d＝0；其中，a,b,c,d为平面模型的拟合参数，保存参数a,b,c作为当前拟合平面的法向量。优选地，所述步骤7)中点到平面A的距离计算公式为：优选地，所述步骤9)中相邻平面法向量夹角|α|和投影距离|dij|的计算公式如下：相邻平面法向量夹角：投影距离：dij|＝max(|d·ni|,|d·nj|)＜dthreshold②其中，分别表示平面i和平面j的法矢量；dij表示两平面之间的投影距离；αthreshold表示角度阈值，dthreshold表示距离阈值；若相邻平面同时不满足公式①和②，则认定当前面为离群面；若相邻平面不满足公式①但满足公式②，则认定此处起伏较大，为地形变化剧烈的陡坡。优选地，所述步骤10)中λ3值与中误差σ值大小关系为：λ3＝2×σ。本发明具有如下优点：本发明基于KD索引树建立点云数据之间的拓扑关系，各点的近邻域数据基于RANSAC算法拟合局部平面，计算点云到各自局部拟合平面的距离，基于统计分析方法去噪，此外，去噪前根据相邻平面的法矢量特征作预判去除明显离群面，并保留陡坡处点云，防止过度去噪。本发明方法能够去除多波束点云数据中的近地表噪声和明显离群噪声数据，同时较好的保留边缘等信息，在保证以上效果的基础上优化设计方案，提高执行效率。附图说明图1为本发明中一种顾及地形特点的多波束点云数据去噪方法的流程示意图；图2为本发明中平面拟合示意图；图3为本发明中平面拓扑结构示意图；图4为本发明中去噪前的点云整体效果图；图5为本发明中去噪后的点云整体效果图；图6为本发明中近地面和明显噪声去噪前效果图；图7为本发明中近地面和明显噪声去噪后效果图；图8为本发明中不一致数据点的去噪前效果图；图9为本发明中不一致数据点的去噪后效果图；图10为利用本发明方法处理图4中区域3的地形边界去噪效果图；图11为利用PCL统计滤波器处理图4中区域3的地形边界去噪效果图。具体实施方式下面结合附图以及具体实施方式对本发明作进一步详细说明：首先给出与本发明相关的几个技术术语的名词解释：(1)点云数据：指使用测量仪器获得目标物体表面特性的海量点的集合，点云数据包括空间三维坐标(XYZ)和反射强度等信息。(2)多波束点云数据：指使用多波束测深系统采集的海底地形表面的点云数据。(3)点云数据一致性：指点云数据特征一致，符合相同的模型，如地形数据点云一致，噪声数据与地形数据不一致。(4)点云去噪：指去除点云数据中的噪声数据。下面将涉及本发明的技术方案拆分为五个技术环节进行分析：环节1基于KD索引树的k近邻域检索为提高点云数据的检索效率，去噪前需要建立点云间的拓扑关系，同时，当遍历下一个邻域时能够保证与当前邻域相邻，从而能够确定所有RANSAC拟合平面之间的拓扑关系。本方法基于KD索引树构建点云之间的拓扑关系，对各点进行近邻域检索。k近邻域检索方式如下：对于KD索引树上任一点，搜索以该点为球心，以r为半径区域内的点云，若搜索到点的数目小于k，则将该区域内所有的点作为该点的k近邻域数据点，否则，选择距离该点最近的k个点作为该点的k近邻域数据点。环节2随机采样一致性算法(RANSAC)拟合局部平面由于水下大部分地形坡度较缓，各点相对较大邻域范围内的非噪声点云可以近似看作平面，且在此范围内噪声较多，本发明基于RANSAC算法思想拟合局部平面。相比于最小二乘算法，RANSAC算法鲁棒性强，能够从包含大量噪声点的数据集中获得高精度的一致性模型参数。基于RANSAC思想拟合平面过程如下：1)从近邻域中随机选择三个点，计算平面方程ax+by+cz+d＝0；然后计算各点到平面的距离di；2)设置阈值t，若di≤t，则认为是非噪声点，统计内点个数；3)重复上述步骤1)和步骤2)m次，选择非噪声点最多的平面；4)特征值法拟合包含非噪声点最多的平面，获得拟合参数a,b,c,d。迭代次数m使用如下公式计算：其中p表示从数据集中随机选取的点均为非噪声点的概率，w表示每次从数据集中选取一个非噪声点的概率。通常p取99％，即以99％的概率选取的点均为非噪声点。环节3基于概率统计方法去噪本发明的技术方案是基于以下规律：计算点到各自RANSAC拟合平面的距离，对大量数据点对应的距离作统计分析，数据近似服从高斯分布。点到RANSAC拟合平面的距离计算公式：基于以上规律和统计分析理论，设定中误差阈值，认为处于给定中误差阈值之外的点云为噪声。环节4共面法向量特征去除离群面多波束点云去噪存在以下问题：1)噪声数据中存在能够拟合与非噪声点一致性良好的模型(本发明指离群面)时，算法失效；2)在全局统计下，陡坡处点云到拟合平面的距离往往处于给定中误差范围之外，导致过度去除。本发明方法在去噪前，结合共面法向量特征解决以上两个问题。由于两个相邻局部拟合平面共面，需要满足以下两个条件：1)共向性：平面法向量的夹角足够小，即朝向一致；2)共面性：平面上任意两点之间的距离在各自平面法向量的投影距离足够小。如图2所示，两个相邻局部拟合平面共面的数学描述如下：相邻平面法向量夹角：投影距离：|dij|＝max(|d·ni|,|d·nj|)＜dthreshold②其中，分别表示平面i和平面j的法矢量；dij表示两平面之间的投影距离；αthreshold表示角度阈值，dthreshold表示距离阈值。首先，确定所有RANSAC拟合平面之间的拓扑关系，如图3所示。然后，结合以上特征设定合理阈值αthreshold和dthreshold，若相邻平面同时不满足公式①和②，则认定当前面为离群面，若相邻平面不满足公式①但满足公式②，则认定此处起伏较大，为地形变化剧烈的陡坡，予以保留。环节5技术方案优化此环节涉及技术方案的一些优化，提高方法的执行效率。1)技术方案实现过程中的优化：①涉及点云数据的传递均使用该点的索引号，而非使用点云的三维坐标信息，降低存储空间；②算法中多次使用哈希表存储键-值，查询哈希表中键对应的值和哈希表中的键是否存在时，时间复杂度接近O(1)，如执行步骤中哈希表1中存储[点云索引号，标记值]键值对，检查键值即当前点是否存在；哈希表2中存储[点云索引号，距离值]键值对，获取索引号对应的键值，其中点云索引号唯一，距离值在下次遍历时更新。2)数据检索优化：①为提高RANSAC算法拟合平面的精度，RANSAC算法阈值之外的各点不再参与后续的平面拟合；②点云的遍历将使大量点云数据被重复检索，导致方法的整体时间复杂度较高。采用以下方案优化：由于地形数据点的k近邻域数据密度高，遍历过程中，任意一点的近邻域中各点对应各自近邻域RANSAC拟合平面基本相同(陡坡除外)，一次性保存最多k个数据点对应的距离，并将保存距离的点标记为已搜索点云，下次遍历时跳过已保存距离的点云数据。基于上述各技术环节，下面给出了本发明详细的技术方案流程，如图1所示。如图1所示，一种顾及地形特点的多波束点云数据去噪方法，包括如下步骤：1)输入点云集cloudPoint，建立KD索引树，遍历点云集cloudPoint中的每一个点searchPoint，执行步骤2)-步骤7)，遍历结束后，执行步骤8)；2)检查searchPoint是否已经被搜索过，检查searchPoint是否是无效点；若searchPoint已经被搜索过或者是无效点，则跳出本次循环，转到步骤1)；否则，转到步骤3)；3)设置搜索半径r，对searchPoint进行k近邻域半径搜索，搜索到点数记为K_Number，将搜索到的点和对应的索引号保存至向量容器pointVector中；4)保存向量容器pointVector中有效数据点的索引号和标记值0至哈希表1中，用于判断当前点searchPoint是否被搜索；5)设置阈值λ1，若K_Number小于λ1，则认定向量容器pointVector中的数据点为噪声数据点，标记为无效点并存储该点的索引号和距离值0至哈希表2中，执行步骤1)；否则，执行步骤6)；6)设置RANSAC平面拟合算法的阈值λ2，对定向量容器pointVector中的点基于RANSAC算法拟合平面，平面记为A；7)计算定向量容器pointVector中索引号对应的点到平面A的距离，并保存索引号和对应距离值至哈希表2表中，转到步骤1)；8)计算哈希表2中不为0的距离的平均值和中误差，分别记为u和σ；9)分别计算相邻平面法向量夹角|α|和投影距离|dij|，若同时不满足共向性和共面性条件，则认定当前拟合平面为离群面并标记相应的点为噪声数据点，若相邻平面不满足共向性但满足共面性，则认定当前拟合平面处地形为起伏较大的陡坡，标记相应点为非噪声数据点；10)设置阈值λ3，若哈希表2中值不为0的索引号对应的距离在[u-λ3·σ,μ+λ3·σ]之间，则认定该索引号对应cloudPoint中的数据点为非噪声数据点，否则认定为噪声数据点；11)分别保存非噪声数据点和噪声数据点数据至PCD格式文件中。下面对图1作如下说明：①规则1：将搜索到的有效点的距离值标记为0(也可以是其他值)，按照(索引号，距离值)的方式存入哈希表2中；②规则2：将有效点的索引号和计算的对应距离值，按照[索引号，距离值]的方式存入哈希表2中；③法矢量对拟合平面预判：详细内容见步骤9)；④外点是指噪声数据点，内点是指非噪声数据点。本发明中涉及的参数及含义和选取方式如下：(1)λ1值确定：由于半径搜索到的数据量相比于地形点云太少，说明此处均为密度低的噪声数据，这些数据使用RANSAC算法拟合平面达不到理想效果。且λ1需要滤掉密度明显低的离群噪声点同时满足RANSAC算法拟合平面模型(至少需要3个点)的要求，因此，阈值λ1的选取宽泛，通常选取为k值的20％左右。(2)λ2值的确定：该阈值确定RANSAC算法结束条件，根据设备的相关参数和实际外业情况大致计算采集点云间的横向和纵向距离，λ2取横向间距和纵向间距的最小值。例如：在多波束测量中，如果平均船速为6节(1海里/小时＝1.852公里/小时)，即3.08m/s。若行驶中每秒15ping，船行方向数据点间距大约0.2m。以多波束R2Sonic2024为例，1ping发送256个波束，若开角设为120°，平均水深为20m，可以得到扫宽大约为70m，船行垂直方向数据点间距大约0.27m。λ2取min{0.2,0.27}＝0.2。(3)k值确定：k值需要满足RANSAC算法计算拟合平面参数的要求，若点数较少，去除不一致噪声效果较差，若点数太多，算法迭代次数增加，导致算法时间复杂度增加。大量试验分析，多波束点云点数k>500效果较好，当然不要超过1000，否则拖慢速度。(4)去噪半径r值确定：为了减少对k值的限制和提高RANSAC拟合平面的精度，r取值较为宽泛，通常取3倍半径。根据λ2值确定中的例子，各个数据点占据区域面积S＝0.2*0.27＝0.054m。若取k＝700，由S·k＝π·(r/3)2，得去噪半径r为10.4m。(5)角度阈值αthreshold和距离阈值dthreshold确定：本方法认定35°以上为陡坡，角度阈值αthreshold选取35°，距离阈值dthreshold取去噪直径，即2r。下面对本发明方法进行实验测试分析：本发明中的编程实现平台基于表1，针对实测多波束数据进行去噪。参数设置如下：k＝700，去噪半径r＝11.5，λ2＝0.2，λ3＝2，λ1＝140，αthreshold＝45，dthreshold＝23。表1算法实验平台电脑(台式)CPU：Intel(R)Core(TM)i5-4590内存：8GB实验环境Qt4.8.7,VS2013编译器，PCLver1.8.0执行结果如表2所示。表2执行结果统计表算法去噪前点数确认非噪声数据点确认噪声数据点用时本发明方法275414273531188342s图4和图5分别为点云去噪前后的整体效果图，图6到图9为各个区域去噪的细节展示。其中，图6和图7为本发明针对近地表和明显离群噪声的去噪前后对比图；图8和图9为针对不一致噪声的去噪前后对比图。从点云效果来看，本发明能够去除多波束点云的近地表噪声和明显离群噪声。由于现有技术中PCL统计滤波器对点云数据去噪的处理能够适应地形特点、执行效率高且去噪效果较优的优点，因此，本发明还与此技术进行了对比分析。(1)原理分析：PCL统计滤波器的基本原理：计算每个点到k近邻点的平均距离作为该点的距离，然后计算所有点对应距离的均值和中误差，若该点在给定的中误差范围内，则认定为非噪声数据点，否则认定为噪声数据点。该方法存在如下问题：由于地形边缘区域的点云数据到k近邻点的平均距离值较大，往往处于中误差范围外，被认定为噪声数据点，这会导致边缘信息难以较好的保留。同时，算法没有考虑地形特征。本发明基本原理：依据地形特征变化进行去噪，结合RANSAC算法鲁棒性强和概率统计方法的优势，计算各个数据点到RANSAC拟合平面的距离作为该点的距离值，然后再计算均值和中误差，若该点在给定的中误差范围内，则认定为非噪声数据点，否则认定为噪声数据点。因而，本发明中的地形信息和边缘信息能够较好的保留。(2)实验分析PCL统计滤波器设置阈值k＝50和2倍中误差，对相同数据去噪。图10为本发明方法的边界去噪效果图，图11为PCL统计滤波器的边界去噪效果图。可以明显发现本发明相比PCL统计滤波器边界的去噪优势。(3)执行效率分析PCL统计滤波器算法对同样数据去噪的结果统计如表4所示：表3程序执行结果统计表算法去噪前点数确认内点确认外点用时PCL统计滤波器算法275414272873254155s本发明方法275414273531188342s由表3可以看出，本发明耗时42s，相比PCL统计滤波器算法，本发明执行效率更高。当然，以上说明仅仅为本发明的较佳实施例，本发明并不限于列举上述实施例，应当说明的是，任何熟悉本领域的技术人员在本说明书的教导下，所做出的所有等同替代、明显变形形式，均落在本说明书的实质范围之内，理应受到本发明的保护。当前第1页1 2 3