CN107123161A - 一种基于narf和fpfh的接触网零全网三维重建方法 - Google Patents

Abstract

本发明公开了一种基于NARF和FPFH的接触网零全网三维重建方法，包括以下步骤：通过安装有Kinect2.0深度摄像机的检测小车采集接触网初始点云；对初始点云进行预处理；针对接触网点云数据，提取两帧点云中对应关键点，采用归一化对齐径向特征NARF对点云图像边界进行检测，提取表面稳定但邻域深度信息发生变化的NARF关键点；采用快速点特征直方图FPFH算法对NARF关键点进行特征描述，确定两帧初始点云之间的关键点对应关系；利用关键点进行点云配准，通过SAC‑IA将误匹配点剔除；采用ICP算法进行精确配准，得到配准后的完全接触网三维点云模型。本发明能够良好建立可视化接触网全网模型，并有效检测其中缺损绝缘子存在的故障。

Description

一种基于NARF和FPFH的接触网零全网三维重建方法

技术领域

本发明涉及电气化铁路接触网三维模型重建以及检测领域，特别是一种基于NARF和FPFH的接触网零全网三维重建方法。

背景技术

接触网是电气化高速铁路供电系统的重要组成部分主要架构，由接触悬挂装置、定位装置、支持装置等部分组成，其性能直接影响高速列车的运行速度与安全，确保接触网参数符合设计规范要求对于高铁可靠的电能传输至关重要。接触网所处铁路沿线恶劣、复杂的工作环境，导致其故障频发，针对接触网设计、施工、状态的检测就显得至关重要。

三维检测目前已经被广泛应用于机械、航空、军工等各种领域，但针对于接触网故障方面的三维检测研究几乎没有。

一种适合在线三维检测的改进算法一文中提供了这样一种方法，即在传统相位测量轮廓术加入了Stoilov算法，使其应用于在线三维检测。传统的算法中含有开方和除法运算，对于CCD摄像头的非线性误差，投射光长，周围环境的光学干扰等比较敏感，容易出现较大误差。因此作者利用Stoilov算法，有效减弱了预案算法公式中的解相误差，极大地提高了在线三维检测的进度。在验证性的实物实验中，基于改进后的stoilov算法重建的三维物体具有非常好的保真度(详见钟立俊、曹益平：一种适合在线三维检测的改进算法，中国激光，2009)。

基于机器视觉的三角螺纹三维检测方法及实验研究一文中，作者根据三角形外螺纹的自己轮廓结构，首先通过CCD面阵摄像机对于螺纹进行表面信息的图像采集，然后对采集到的图像进行预处理，如图像增强处理、螺纹的边缘检测、提取螺纹边缘特征等手段；最后利用提取得到的螺纹特征点实现了外螺纹的三维重加，为外螺纹的三维重建提出了新的可行方案(详见万鹏：基于机器视觉的三角螺纹三维检测方法及实验研究，华南理工大学，2012)。

在三维检测技术在列车车轮检测与维修中的应用研究一文中，将三维重建技术应用于铁路系统中的列车车轮检车，其效率将远高于传统的人工检测手段。该方法源自逆向工程原理，结合计算机图像辅助技术来进行全参数的车轮三维检测。该方法在铁路车轮的检测规范前提下，论证了最优解决方法，将车轮的检测与维修紧密结合在一起，为铁路系统车轮检修提供了新的思路(详见程宏钊：三维检测技术在列车车轮检测与维修中的应用研究，西南交通大学，2013)。

基于三维模型的接触网检测技术研究一文中，利用光学扫描仪获得接触网各零部件的点云数据，并运用一定的点云配准算法将不用视角下测量得到的多片点云转换合并到同一坐标系下形成一个完整的数据点云，然后经曲面重建和渲染获得各部件三维模型(详见徐建芳：基于三维模型的接触网检测技术研究，西南交通大学，2014)。

总的来说，目前针对于接触网图像检测主要基于二维接触网图像，二维图像的检测技术蕴含的图像信息相对单一，存在图像检测的死角，无法快速、精确得检测接触网实时故障。

发明内容

本发明所要解决的技术问题是提供一种基于NARF和FPFH的接触网零全网三维重建方法，该方法采用三维点云处理技术对获取到的接触网全网点云数据进行预处理、配准、融合，获得三维模型，并提取其中缺损绝缘子，检测其故障；其能够良好建立可视化接触网全网模型，并有效检测其中缺损绝缘子存在的故障。

为解决上述技术问题，本发明采用的技术方案是：

一种基于NARF和FPFH的接触网零全网三维重建方法，包括以下步骤：

步骤1：通过安装有Kinect2.0深度摄像机的检测小车采集接触网初始点云，获得的初始点云数据是pcd格式文件，其包括点云位置、颜色和深度；

步骤2：对初始点云进行预处理，包括点云去噪、点云分割、点云精简、点云融合和点云派生；

步骤3：针对接触网点云数据，提取两帧点云中对应关键点，采用归一化对齐径向特征NARF对点云图像边界进行检测，提取表面稳定但邻域深度信息发生变化的NARF关键点，这些变化包括表面变化系数以及变化的主方向；

步骤4：采用快速点特征直方图FPFH算法对NARF关键点进行特征描述，确定两帧初始点云之间的关键点对应关系；

步骤5：利用关键点进行点云配准，通过SAC-IA将误匹配点剔除；

步骤6：采用ICP算法进行精确配准，得到配准后的完全接触网三维点云模型。

进一步的，还包括步骤7：

步骤7：利用点云分割手段，将点云配准后的接触网模型中的缺损腕臂绝缘子提取出来，估计提取得到的绝缘子表面法线，利用其表面的法线信息对缺损绝缘子进行检测。

进一步的，所述步骤3具体为：

步骤3.1：点云数据中某点曲率大的位置对应不同零件的结构分界点，边界兴趣值I(p)为：

I(p)＝I₁(p)·I₂(p)

式中：

I₁(p)为表面变化强度估计因子，I₂(p)为主方向辨别因子；p为边界点，n为点p的邻元素，σ为空间尺度，w_n为邻元素n对应的权重值，γ为点n对应一维角度；对于接触网曲面边缘点，w_n取值为1，而对于其它点，w_n取值为1-(1-ρ)³，ρ为p及其周边点进行主成分分析后得到的曲率值。

进一步的，所述步骤4具体为：

步骤4.1：对点云中任意点P_q查询其k领域内所有邻近点；

步骤4.2：对点P的k领域中任意一对点P_s和P_t估计法线n_s和n_t，其中s！＝t，在其中一个点上定义一个固定的u,v,w局部坐标系，计算n_s和n_t之间的偏差

步骤4.3：在步骤4.2的u,v,w局部坐标系中，用一组角度来估计法线n_s和n_t之间的偏差：

α＝v·n_t

式中，||P_t-P_s||＝d表示点P_s和P_t两点间的欧式距离，则P_q的所有k领域任意两点P_s和P_t的关系通过(α,β,θ,d)表示；k领域所有的四组值以统计的方式放入直方图，得到P_q的特征描述子；

步骤4.4：对于每一个查询点P_q，计算该点与它邻域点之间的一个元组α,θ；计算FPFH，重新确定每个点的k邻域，使用邻近SPFH计算P_q的最终的FPFH，其计算公式为式中，权重ω_k表示在给定的度量空间中，查询点P_q与其邻近点P_k之间的距离；

计算点云帧中全部点云关键点的FPFH特征描述子后，获得配准所需的关键点对应关系。

进一步的，所述步骤5中点云配准具体为：

两帧接触网点云对应关系确定后，对点云进行配准，其划分为粗配准和精配准；粗配准采用SAC-IA，将误匹配点剔除后，利用算法进行精确配准；经一系列点云配准获得完整的接触网零三维点云数据。

进一步的，所述步骤7具体为：

步骤7.1：分析一个协方差矩阵的特征矢量和特征值，而这个协方差矩阵C从每一个查询点P_i的邻近点中构造得到，具体的协方差矩阵如下：

上式中，k是查询点P_i邻近点的数目，是最邻近元素的三维质心，λ_j是协方差矩阵的第j个特征值，是第j个特征向量；

步骤7.2：确立切平面

对于一帧点云数据，其点云数量为n，若要计算其中某一个查询点P_i的法线，先设置平面方程：

ax+by+cz+d＝0

对于上式中进行条件约束，当满足a²+b²+c²＝1时，得到平面参数a、b、c、d；让点P_i周围的k个临近点到达该平面距离的平方和数值最小，即满足：

式中，d_i是点云数据集中任一查询点P_i＝(x_i,y_i,z_i)到对应平面的距离d_i＝|ax_i+by_i+cz_i-d|；利用拉格朗日乘子法求解极值，使得e→min，得到函数：

上式中等号两边对于d求偏导，并且令偏导数为零，得到：

式中，令质心为则

d_i＝|aΔx_i+bΔy_i+cΔz_i|

再对中等号两边对于a、b、c求偏导数，得到

求解平面参数a、b、c，转化为求解矩阵A的特征值和特征向量；矩阵A是一个三节对称矩阵，其特征值的求解公式为：

在a²+b²+c²＝1的约束条件下，求得利用质心求得d；

步骤7.3：调整法线方向

判断所有法线方向是否一致，假如实际视点为V_p，点云中所有任一法线n_i都指向该视点，则认为所有法线方向一致，即满足：

n_i·VP＞0

式中，VP是视点V_p到查询点P_i的向量；如果计算结果不满足式n_i·VP＞0的要求，说明点P_i的法向量n_i与点P_i到视点V_p的向量之间的向量夹角大与90°，点P_i的法向量n_i应该反向；根据此点云配准方法进行接触网三维重建过程中点云的配准，获得接触网全网三维模型，并对其中缺损的绝缘子进行故障检测。

与现有技术相比，本发明的有益效果是：

1、本发明采集接触网点云的方式是通过铁轨上架设的安装了Kinect2.0深度摄像机的检测小车来采集接触网初始点云，这种采集方式符合接触网实时检测要求，最后通过一系列处理获得的接触网全网三维点云模型能够高度还原实际接触网真实状态。

2、本发明利用绝缘子表面法线信息的三维检测方法能够直观、有效地检测出绝缘子缺损信息，满足当今对于接触网实时检测的客观需求。

附图说明

图1是接触网点云采集现场图。

图2是接触网点云数据采集原理图。

图3是初始Frame1点云可视化图。

图4是初始Frame2点云可视化图。

图5是Frame1和Frame2点云相对位置图。

图6是RadiusOutlierRemoval滤波原理示意图。

图7是Frame1下采样前的示意图。

图8是Frame1下采样后的示意图。

图9是Frame1的NARF关键点图。

图10是Frame2的NARF关键点图。

图11是Frame1某关键点的FPFH特征直方图。

图12是Frame2某关键点的FPFH特征直方图。

图13是接触网点云配准结果图。

图14是接触网点云中提取斜腕臂正常绝缘子点云示意图。

图15是接触网点云中提取斜腕臂缺损绝缘子点云示意图。

图16是正常棒式绝缘子的法线图。

图17是却损棒式绝缘子的法线图。

图18是正常绝缘子水平横截面的法线图。

图19是缺损绝缘子水平横截面的法线图。

图20是绝缘子相邻法线夹角值对比图。

具体实施方式

下面结合附图和具体实施方式对本发明作进一步详细的说明，本发明方法详述如下：

1、接触网三维点云数据获取及预处理

由于接触网特殊的存在形式，无法对其进行360°环绕拍摄，通过在铁轨上架设的安装了Kinect2.0深度摄像机的检测小车来采集接触网初始点云，获得初始点云数据是包涵点云位置、颜色、深度等信息的pcd格式文件；接触网初始点云数据包涵复杂的环境，以及很多噪点，需要进行初始点云进行预处理，包括点云去噪、点云分割、点云精简、点云融合、点云派生等多种处理方式。

2、对上述处理获得的单独、干净的接触网点云数据提取两帧点云中对应关键点，利用归一化对齐径向特征(Normal Aligned Radial Feature，NARF)通过点云图像边界检测结果，提取表面稳定但邻域实质变化极大的NARF兴趣点；接着，利用快速点特征直方图(Fast Point Feature Histograms，FPFH)算法对于NARF关键点进行特征描述，以此确定两帧初始点云之间的并确定匹配关系；通过SAC-IA(Sample Consensus Initial Alignment，SAC-IA)将误匹配点剔除后，最终利用ICP(Iterative Closest Point，ICP)算法进行精确配准，得到配准后的完全接触网三维点云模型。

3、利用点云分割手段从点云配准后的接触网模型中的缺损腕臂绝缘子提取出来，估计提取得到的绝缘子表面法线，利用其表面的法线信息，对缺损绝缘子进行检测。

具体如下：

1)NARF关键点检测

NARF算法通过点云图像边界检测结果，提取表面稳定但邻域实质变化极大的边缘兴趣点，重点分析对象结构，具有典型性和旋转不变形等优点，对于零件结构稀疏性较强的接触网点云数据较为适合。

接触网零件多为体表外形均匀变化，点云数据中某点曲率大的位置对应不同零件的结构分界点。边界兴趣值I(p)为：

I(p)＝I₁(p)·I₂(p)

式中：

I₁(p)为表面变化强度估计因子，I₂(p)为主方向辨别因子。p为边界点，n为点p的邻元素，σ为空间尺度，w_n为邻元素n对应的权重值，γ为点n对应一维角度。对于接触网曲面边缘点，w_n取值为1，而对于其它点，w_n取值为1-(1-ρ)³，ρ为p及其周边点进行主成分分析(Principal Component Analysis，PCA)后得到的曲率值。

2)快速点特征直方图FPFH特征描述

A、点特征直方描述子(Point Feature Histograms，PFH)，通过统计查询点与k领域之间以及k领域内所有点之间法线关系的变化情况，来描素对象的几何特征。一个查询点P_q的PFH计算影响区域是以它自身为圆心的范围，圆的半径为r，P_q的所有k领域元素全部相互连接在一个网络中，通过计算领域内任意两点以及它们的法线关系，最终获得PFH直方图。具体过程为：对点云中任意点P_q查询其k领域内所有邻近点。

B、对点P的k领域中任意一对点P_s和P_t(s！＝t)估计法线n_s和n_t，在其中一个点上定义一个固定的u,v,w局部坐标系，计算n_s和n_t之间的偏差，

C、在B中定义的u,v,w局部坐标系中，用一组角度来估计法线n_s和n_t之间的偏差：

α＝v·n_t

式中，||P_t-P_s||＝d表示点P_s和P_t两点间的欧式距离，那么P_q的所有k领域任意两点P_s和P_t的关系就可以通过(α,β,θ,d)表示。k领域所有的四组值以统计的方式放入直方图，就可以得到P_q的特征描述子。

D、快速点特征直方图(FPFH)由PFH改进而来，FPFH将算法的计算复杂度有效降低，但又保留了PFH的大部分识别特性，具有较好的鲁棒性和实时性。为了简化直方图的特征计算，其计算过程如下：简化点特征直方图(Simple Point Feature Histograms,SPFH)，即对于每一个查询点P_q，计算该点与它邻域点之间的一个元组α,θ；计算FPFH，重新确定每个点的k邻域，使用邻近SPFH来计算P_q的最终的FPFH，其计算公式：式中，权重ω_k在给定的度量空间中查询点P_q和其邻近点P_k之间的距离。计算点云帧中全部点云关键点的FPFH特征描述子后，可获得配准所需的关键点对应关系。

3)点云配准

两帧接触网点云对应关系确定后，对点云进行配准，可划分为粗配准和精配准两个过程；粗配准采用SAC-IA将误匹配点剔除后，最终利用算法进行精确配准；经一系列点云配准获得完整的接触网零三维点云数据。

4)基于绝缘子表面法线的三维检测

A、对于物体的表面法线估计的方法有很多，在此对于其中最简单的一种进行介绍，其核心思想：计算一个物体表面某一点的法线近似于估计该点相切面法线的问题，因此可以将此问题转化为一个最小二乘法平面拟合估计问题。基于此原理，就可以计算得到物体表面法线的解决方案，即分析一个协方差矩阵的特征矢量和特征值，而这个协方差矩阵C可以从每一个查询点P_i的邻近点中构造得到，具体的协方差矩阵如下：

上式中，k是查询点P_i邻近点的数目，是最邻近元素的三维质心，λ_j是协方差矩阵的第j个特征值，是第j个特征向量。下面物体法线估计的具体方法。

B、确立切平面

对于一帧点云数据，其点云数量为n，若要计算其中某一个查询点P_i的法线，那么先设置平面方程：

ax+by+cz+d＝0

一般来说点云数据集中的点在方向上存在偏差，为了尽量消除这些偏差，对于上式中进行条件约束，当满足a²+b²+c²＝1时，可以得到平面参数a、b、c、d。为了能够获得最佳的拟合平面，此时应该让点P_i周围的k个临近点到达该平面距离的平方和数值最小，即满足：

式中，d_i是点云数据集中任一查询点P_i＝(x_i,y_i,z_i)到对应平面的距离d_i＝|ax_i+by_i+cz_i-d|。利用拉格朗日乘子法求解极值，可以使得e→min，得到函数：

上式中等号两边对于d求偏导，并且令偏导数为零，得到：

式中，令质心为则

d_i＝|aΔx_i+bΔy_i+cΔz_i|

此时，再对中等号两边对于a、b、c求偏导数，得到

此时，求解平面参数a、b、c，就可以转化为求解矩阵A的特征值和特征向量。矩阵A式一个三节对称矩阵，其特征值的求解公式为：

在之前a²+b²+c²＝1的约束条件下，可求得因此，e的最小值就是矩阵A的最小特征值，它对应的特征向量为平面参数a、b、c，利用质心就可以求得d。

C、调整法线方向：利用上一小节中的原理可以计算得到点云中所有点的法向量，但是得到的法向量的方向往往无法一致。为了保持法向量方向的一致性，需要作进一步的处理。判断所有法线方向是否一致，假如实际视点为V_p，点云中所有任一法线n_i都指向该视点，则认为所有法线方向一致，即满足：

n_i·VP＞0

式中，VP是视点V_p到查询点P_i的向量。如果计算结果不满足式(4-10)的要求，说明点P_i的法向量n_i与点P_i到视点V_p的向量之间的向量夹角大与90°，点P_i的法向量n_i应该反向。

根据这一点云配准方法进行接触网三维重建过程中点云的配准，获得接触网全网三维模型，并对其中缺损的绝缘子进行故障检测。

图1为接触网三维点云采集现场

A、接触网点云数据获取点云数据采集包括光学法、计算机视觉法、断层成像法、时间飞行法、光编码技术等，本发明采用时间飞行法进行接触网点云采集。由于接触网特殊的存在形式，无法对其进行360°环绕拍摄，只能通过在铁轨上架设的安装了深度摄像机的检测小车来采集接触网初始点云，本发明中接触网点云采集原理图如图2所示。在检测小车行进的过程中，连续采样，由于小车在向前移动，因此不同帧点云深度信息不同，仰角在变化，配准不但涉及平移，还有旋转，难度更大。采样获得的其中两帧接触网初始点云如图3、图4所示，图3为点云Frame1，图4为点云Frame2。

B、接触网点云预处理

原始接触网点云数据量庞大，离散噪点对特征提取影响较大，需要对原始云进行预处理，包括去噪、背景分割、下采样等，经过预处理后的两帧点云在空间中的相对位置如图5所示。

a、点云去噪

受采集设备、环境等干扰，接触网原始点云包含测量噪声所引入的离群点。利用条件移除滤波器(Conditional Removal滤波器)定义接触网点云的密度阈值，当原始点云中点云密度小于阈值时为噪点。

b、点云分割

点云分割：接触网初始点云中包含自然背景、环境杂物等无用信息，利用半径滤波器(Radius Outlier滤波器)比接触网几何中心为球心设置一个半径和点云数量限值，并计算落在以给出半径画球内部的点云数量，当计算点云数量大于给出限值时保留改点，小于限值则剔除。其原理如图6，图中若限值为1个邻近点，图6中左侧圆的圆心处的点删除，若限制为2个邻近点，则左侧圆和右侧圆各自圆心处的点都删除。

c、点云下采样

原始采集点云较为密集，对后续法线计算和分析造成困难。利用体素滤波器(Voxel Grid滤波器)为原点云创建一个三维体素栅格(可以体素栅格理解为微小的三维立方体的集合)，每个体素栅格中，用体素中所有点的重心来近似体素中其它点，这样该体素就内所有点就用一个重心点最终表示，对于所有体素处理后得到过滤后点云。下采样前后Frame1、Frame2两帧点云数如表1所示，其中Frame1下采样前后的点云对比图如图7、图8所示，从图中可以看出，下采样后，同一点云数据同一位置的点数量明显减少，但其几何结构并未发生改变，这对于提高后续的点云处理的速率有很大的好处。

表1预处理前后点云数量变化图

C、接触网点云数据配准

通过NARF和FPFH算法得到Frame1和Frame2两帧点云的对应关键点对后，利用采样一致性初始配准算法(SAC-IA)对两帧点云进行初配准，利用ICP点云配准算法进行精确配准。

Frame1和Frame2两帧点云数据的NARF关键点，结果如图9、图10所示；Frame1和Frame2两帧的FPFH特征描述如图11、图12所示，图中，横坐标的数值是划分的从0-32的33个统计区间，纵坐标的数值是对应统计区间中包含的点(p₁,p₂,…,p_k)的个数。从图中看出Frame1和Frame2的一组对应关键点的FPFH特征直方图十分相似，虽然两张图的折线并不完全一致，但是可以看出折线的基本走势都一致的，纵坐标各个峰值所对应的横坐标也大致相同，这说明Frame1和Frame2两者选择的关键点的特征是极为相似的，这对于之后点云配准的精度十分有利。计算点云帧中全部点云关键点的FPFH特征描述子后，可获得配准所需的关键点对应关系。

a、SAC-IA初始配准

采样一致性初始配准算法主要分为两个部分：贪婪的初始配准对准方法，采样一致性算法。因为点云具有内部旋转不变形的特征，因此使用贪婪初始对准算法具有非常良好的鲁棒性。但是贪婪的初始对准算法计算复杂度较高，而且有可能只能得到局部的最优解，因此，还需要采用采样一致性方法，试图确定相同的对应几何关系，而不必计算有限个对应关系的所有组合。

SAC-IA基本过程如下：

(1)从点云Frame1中选取个样本点，同时计算样本点配对距离大于预设阈值的最小距离；

(2)对于个样本点，在点云Frame2中分别找到所有满足相似度条件的点云，并随机选择其中一些来计算采样点的对应关系；

(3)根据连个点云集合的对应关系，计算刚体变换矩阵，并通过计算度量错误来检验转换举证的质量。

错误度量可以由Huber评价公式检验：

(4)重复以上三个步骤，直到达到最佳误差度量。

根据上述计算原理，通过SAC-IA将误匹配点剔除，得到了SAC-IA初始转换矩阵如下：

b、ICP精确配准

通过SAC-IA算法去除了Frame1和Frame2两帧点云中存在的错误对应关系，并求得了粗配准转换矩阵后，为实现高精度点云配准，采用目前应用最广的ICP点云配准算法进行精确配准，并得到ICP精确配准的转化矩阵。

ICP精确配准的过程可以如下划分：

(1)接触网点集重采样：将Frame1和Frame2中所有点分别记为目标点云X＝{x₁,x₂,…,x_m}，参考点云Y＝{y₁,y₂,…,y_n}，其中X中点云数为m，Y中点云数为n，且m≤n。

(2)匹配点选择：利用四元素法设定初值，并通过Kd-Tree加速寻找最近点云。设旋转变换向量为单位四元数q_R＝[q_x,q_y,q_z,q_w]，其中q_x≥0，并且q_x+q_y+q_z+q_w＝1，平移变换向量q_T＝[t_x,t_y,t_z]。Frame1与Frame2的两片点云重心为

如此，可以计算获得两帧点云配准的最佳旋转矩阵R(q_R)和最佳平移向量q_T，

(3)赋予点对权值：Frame1与Frame2已匹配点对的权重通过它们的法向量来确定。设已匹配的对点的法向量分别为n₁和n₂，则权重为：

w＝n₁·n₂

(4)点对剔除：Frame1与Frame2并不是完全重合，所以一般将包含边缘的错误匹配关系点对剔除。

(5)选择合适的误差度量函数：采用基于最小二乘法的迭代方法实现坐标变换矩阵优化，使误差函数最小。误差函数定义为：

(6)最优化：重复上述步骤，实现误差最小化。利用最小误差获得的旋转矩阵R(q_R)和平移向量q_T将Frame1和Fame2配准，最终获得精确配准矩阵：

c、接触网点云配准结果

在经过关键点提取、描述，使用SAC-IA误匹配点剔除，ICP精确配准，最终配准效果如图13所示，图中相对深度位置靠后的为初始段云Frame1，目标点云Frame2与配准后的生成的新点云几乎重合在一起，但在边缘处可以看出配准后的点云更加完整。

D、基于接触网绝缘子表面法线的三维检测

a、提取棒式绝缘子

将接触网点云中的斜腕臂上的棒式绝缘子通过点云分割的方式提取出来。图14中展示了从接触网点云中提取正常绝缘子点云的过程；图15中展示了从接触网点云中提取缺损绝缘子点云的过程。

b、绝缘子法线计算

利用绝缘子表面法线信息对缺损的绝缘子进行检测，并将检测的信息量化。根据求解物体表面法线原理，分别对图14中正常的绝缘子和图15中破损的绝缘子求解它们的法线，其法线效果如图16、图17所示。

c、棒式绝缘子主体是由直径不同的两种圆柱形交替堆叠构成的长棒式物体，而且呈完全的中心对称，因此理论上其点云的表面法线一定指向其所在水平面圆的圆心。但实际采集获得的绝缘子点云，其表面在存在点云的缺失、重叠等情况，所以即使在滤波、去燥后也不可能完全的平顺、光滑，导致其表面法线也不可能完全指向圆心。图18和图19分别展示了正常的绝缘子和缺损的绝缘子在水平方向的横截面的法线图。

图18和图19中的右边是各自绝缘子水平横截面的法线图。图中两个横截面都不是完整的圆形，这是因为在如图2所示采集点云的过程中只拍摄了小车向前行进过程中的接触网点云数据，当小车经过接触网后并没有拍摄，所以采集到的点云不是360°的。

计算图18中正常绝缘子横截面中的所有相邻法线向量之间的夹角α_m，获得的所有夹角集{α₁,a₂…a_m}；同理，计算图19中缺损绝缘子横截面中的所有相邻法线向量之间的夹角β_n，获得的所有夹角集{β₁,β₂…β_n}。

将以上获得的角度集合{α₁,a₂…a_m}和{β₁,β₂…β_n}数据归纳整理后画出折线图，可以看出两者之间的差异，如图20所示。

图20的折线图中，横坐标表示的是两个绝缘子横截面表面相邻法线夹角个数，在此图中选择了{α₁,a₂…a_m}和{β₁,β₂…β_n}中各200个夹角值作为统计样本；纵坐标表示的是角度。

从图中可以看出，虚线折线代表的正常绝缘子横截面表面相邻法线之间的夹角角度，由于绝缘子本身点云表面的不完全光滑，这条折线在一定幅度内起伏，但是其纵坐标对应的各个角度值的方差较小，说明正常绝缘子表面向量法线之间的夹角值相差不多。图中实线折线代表的缺损绝缘子横截面表面相邻法线之间的夹角角度，这条折线基本上和虚线折线在同一水平线上起伏，除了存在一个远远超过其他角度值的非常规峰值，说明此处两条法线之间的夹角非常大，这个位置就是缺损绝缘子的缺损位置。从检测效果可知本发明能够获得较为理想的检测出接触网中零部件的缺损信息。

Claims (6)

1.一种基于NARF和FPFH的接触网零全网三维重建方法，其特征在于，包括以下步骤：

步骤1：通过安装有Kinect2.0深度摄像机的检测小车采集接触网初始点云，获得的初始点云数据是pcd格式文件，其包括点云位置、颜色和深度；

步骤2：对初始点云进行预处理，包括点云去噪、点云分割、点云精简、点云融合和点云派生；

步骤3：针对接触网点云数据，提取两帧点云中对应关键点，采用归一化对齐径向特征NARF对点云图像边界进行检测，提取表面稳定但邻域深度信息发生变化的NARF关键点，这些变化包括表面变化系数以及变化的主方向；

步骤4：采用快速点特征直方图FPFH算法对NARF关键点进行特征描述，确定两帧初始点云之间的关键点对应关系；

步骤5：利用关键点进行点云配准，通过SAC-IA将误匹配点剔除；

步骤6：采用ICP算法进行精确配准，得到配准后的完全接触网三维点云模型。

2.如权利要求1所述的一种基于NARF和FPFH的接触网零全网三维重建方法，其特征在于，还包括步骤7：

步骤7：利用点云分割手段，将点云配准后的接触网模型中的缺损腕臂绝缘子提取出来，估计提取得到的绝缘子表面法线，利用其表面的法线信息对缺损绝缘子进行检测。

3.如权利要求1所述的一种基于NARF和FPFH的接触网零全网三维重建方法，其特征在于，所述步骤3具体为：

步骤3.1：点云数据中某点曲率大的位置对应不同零件的结构分界点，边界兴趣值I(p)为：

I(p)＝I₁(p)·I₂(p)

式中：

<mrow> <msub> <mi>I</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>p</mi> <mo>)</mo> </mrow> <mo>=</mo> <munder> <mrow> <mi>m</mi> <mi>i</mi> <mi>n</mi> </mrow> <mi>i</mi> </munder> <mo>(</mo> <mn>1</mn> <mo>-</mo> <msub> <mi>w</mi> <msub> <mi>n</mi> <mi>i</mi> </msub> </msub> <mi>m</mi> <mi>a</mi> <mi>x</mi> <mo>(</mo> <mrow> <mn>1</mn> <mo>-</mo> <mfrac> <mrow> <mn>10</mn> <mo>&amp;CenterDot;</mo> <mo>|</mo> <mo>|</mo> <mi>p</mi> <mo>-</mo> <msub> <mi>n</mi> <mi>i</mi> </msub> <mo>|</mo> <mo>|</mo> </mrow> <mi>&amp;sigma;</mi> </mfrac> <mo>,</mo> <mn>0</mn> </mrow> <mo>)</mo> </mrow>

<mrow> <msub> <mi>I</mi> <mn>2</mn> </msub> <mo>=</mo> <munder> <mrow> <mi>m</mi> <mi>a</mi> <mi>x</mi> </mrow> <mrow> <mi>i</mi> <mo>-</mo> <mi>j</mi> </mrow> </munder> <mrow> <mo>(</mo> <mi>f</mi> <mo>(</mo> <msub> <mi>n</mi> <mi>i</mi> </msub> <mo>)</mo> <mi>f</mi> <mo>(</mo> <msub> <mi>n</mi> <mi>j</mi> </msub> <mo>)</mo> <mo>(</mo> <mrow> <mn>1</mn> <mo>-</mo> <mo>|</mo> <mi>c</mi> <mi>o</mi> <mi>s</mi> <mrow> <mo>(</mo> <mrow> <msubsup> <mi>&amp;gamma;</mi> <msub> <mi>n</mi> <mi>i</mi> </msub> <mo>&amp;prime;</mo> </msubsup> <mo>-</mo> <msubsup> <mi>&amp;gamma;</mi> <msub> <mi>n</mi> <mi>j</mi> </msub> <mo>&amp;prime;</mo> </msubsup> </mrow> <mo>)</mo> </mrow> <mo>|</mo> </mrow> <mo>)</mo> <mo>)</mo> </mrow> </mrow>

<mrow> <mi>f</mi> <mrow> <mo>(</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>=</mo> <msqrt> <mrow> <msub> <mi>w</mi> <mi>n</mi> </msub> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <mfrac> <mrow> <mn>2</mn> <mo>&amp;CenterDot;</mo> <mo>|</mo> <mo>|</mo> <mi>p</mi> <mo>-</mo> <mi>n</mi> <mo>|</mo> <mo>|</mo> </mrow> <mi>&amp;sigma;</mi> </mfrac> <mo>-</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> </mrow> </msqrt> </mrow>

I₁(p)为表面变化强度估计因子，I₂(p)为主方向辨别因子；p为边界点，n为点p的邻元素，σ为空间尺度，w_n为邻元素n对应的权重值，γ为点n对应一维角度；对于接触网曲面边缘点，w_n取值为1，而对于其它点，w_n取值为1-(1-ρ)³，ρ为p及其周边点进行主成分分析后得到的曲率值。

4.如权利要求1所述的一种基于NARF和FPFH的接触网零全网三维重建方法，其特征在于，所述步骤4具体为：

步骤4.1：对点云中任意点P_q查询其k领域内所有邻近点；

步骤4.2：对点P的k领域中任意一对点P_s和P_t估计法线n_s和n_t，其中s！＝t，在其中一个点上定义一个固定的u,v,w局部坐标系，计算n_s和n_t之间的偏差

步骤4.3：在步骤4.2的u,v,w局部坐标系中，用一组角度来估计法线n_s和n_t之间的偏差：

α＝v·n_t

<mrow> <mi>&amp;beta;</mi> <mo>=</mo> <mi>u</mi> <mo>&amp;CenterDot;</mo> <mfrac> <mrow> <mo>(</mo> <msub> <mi>P</mi> <mi>t</mi> </msub> <mo>-</mo> <msub> <mi>P</mi> <mi>s</mi> </msub> <mo>)</mo> </mrow> <mrow> <mo>|</mo> <mo>|</mo> <msub> <mi>P</mi> <mi>t</mi> </msub> <mo>-</mo> <msub> <mi>P</mi> <mi>s</mi> </msub> <mo>|</mo> <mo>|</mo> </mrow> </mfrac> </mrow>

<mrow> <mi>&amp;theta;</mi> <mo>=</mo> <mi>arctan</mi> <mrow> <mo>(</mo> <mi>w</mi> <mo>&amp;CenterDot;</mo> <mover> <msub> <mi>n</mi> <mi>t</mi> </msub> <mo>&amp;RightArrow;</mo> </mover> <mo>,</mo> <mi>u</mi> <mo>&amp;CenterDot;</mo> <mover> <msub> <mi>n</mi> <mi>t</mi> </msub> <mo>&amp;RightArrow;</mo> </mover> <mo>)</mo> </mrow> </mrow>

式中，||P_t-P_s||＝d表示点P_s和P_t两点间的欧式距离，则P_q的所有k领域任意两点P_s和P_t的关系通过(α,β,θ,d)表示；k领域所有的四组值以统计的方式放入直方图，得到P_q的特征描述子；

步骤4.4：对于每一个查询点P_q，计算该点与它邻域点之间的一个元组α,θ；计算FPFH，重新确定每个点的k邻域，使用邻近SPFH计算P_q的最终的FPFH，其计算公式为式中，权重ω_k表示在给定的度量空间中，查询点P_q与其邻近点P_k之间的距离；

计算点云帧中全部点云关键点的FPFH特征描述子后，获得配准所需的关键点对应关系。

5.如权利要求1所述的一种基于NARF和FPFH的接触网零全网三维重建方法，其特征在于，所述步骤5中点云配准具体为：

两帧接触网点云对应关系确定后，对点云进行配准，其划分为粗配准和精配准；粗配准采用SAC-IA，将误匹配点剔除后，利用算法进行精确配准；经一系列点云配准获得完整的接触网零三维点云数据。

6.如权利要求2所述的一种基于NARF和FPFH的接触网零全网三维重建方法，其特征在于，所述步骤7具体为：

步骤7.1：分析一个协方差矩阵的特征矢量和特征值，而这个协方差矩阵C从每一个查询点P_i的邻近点中构造得到，具体的协方差矩阵如下：

<mrow> <mi>C</mi> <mo>=</mo> <mfrac> <mn>1</mn> <mi>k</mi> </mfrac> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>k</mi> </munderover> <mo>&amp;CenterDot;</mo> <mrow> <mo>(</mo> <msub> <mi>P</mi> <mi>i</mi> </msub> <mo>-</mo> <mover> <mi>P</mi> <mo>&amp;OverBar;</mo> </mover> <mo>)</mo> </mrow> <mo>&amp;CenterDot;</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>P</mi> <mi>i</mi> </msub> <mo>-</mo> <mover> <mi>P</mi> <mo>&amp;OverBar;</mo> </mover> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mo>,</mo> <mi>C</mi> <mo>&amp;CenterDot;</mo> <mover> <msub> <mi>v</mi> <mi>j</mi> </msub> <mo>&amp;RightArrow;</mo> </mover> <mo>=</mo> <msub> <mi>&amp;lambda;</mi> <mi>j</mi> </msub> <mo>&amp;CenterDot;</mo> <mover> <msub> <mi>v</mi> <mi>j</mi> </msub> <mo>&amp;RightArrow;</mo> </mover> <mo>,</mo> <mi>j</mi> <mo>&amp;Element;</mo> <mo>{</mo> <mn>0</mn> <mo>,</mo> <mn>1</mn> <mo>,</mo> <mn>2</mn> <mo>}</mo> </mrow>

上式中，k是查询点P_i邻近点的数目，是最邻近元素的三维质心，λ_j是协方差矩阵的第j个特征值，是第j个特征向量；

步骤7.2：确立切平面

对于一帧点云数据，其点云数量为n，若要计算其中某一个查询点P_i的法线，先设置平面方程：

ax+by+cz+d＝0

对于上式中进行条件约束，当满足a²+b²+c²＝1时，得到平面参数a、b、c、d；让点P_i周围的k个临近点到达该平面距离的平方和数值最小，即满足：

<mrow> <mi>e</mi> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>n</mi> </munderover> <msubsup> <mi>d</mi> <mi>i</mi> <mn>2</mn> </msubsup> <mo>&amp;RightArrow;</mo> <mi>m</mi> <mi>i</mi> <mi>n</mi> </mrow> 2

式中，d_i是点云数据集中任一查询点P_i＝(x_i,y_i,z_i)到对应平面的距离d_i＝|ax_i+by_i+cz_i-d|；利用拉格朗日乘子法求解极值，使得e→min，得到函数：

<mrow> <mi>f</mi> <mo>=</mo> <mi>e</mi> <mo>-</mo> <mi>&amp;lambda;</mi> <mrow> <mo>(</mo> <msup> <mi>a</mi> <mn>2</mn> </msup> <mo>+</mo> <msup> <mi>b</mi> <mn>2</mn> </msup> <mo>+</mo> <msup> <mi>c</mi> <mn>2</mn> </msup> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mi>&amp;Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>n</mi> </munderover> <msup> <mi>d</mi> <mn>2</mn> </msup> <mo>-</mo> <mi>&amp;lambda;</mi> <mrow> <mo>(</mo> <msup> <mi>a</mi> <mn>2</mn> </msup> <mo>+</mo> <msup> <mi>b</mi> <mn>2</mn> </msup> <mo>+</mo> <msup> <mi>c</mi> <mn>2</mn> </msup> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>

上式中等号两边对于d求偏导，并且令偏导数为零，得到：

<mrow> <mi>d</mi> <mo>=</mo> <mfrac> <mn>1</mn> <mi>n</mi> </mfrac> <mrow> <mo>(</mo> <mi>a</mi> <munderover> <mi>&amp;Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>n</mi> </munderover> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>+</mo> <mi>b</mi> <munderover> <mi>&amp;Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>n</mi> </munderover> <msub> <mi>y</mi> <mi>i</mi> </msub> <mo>+</mo> <mi>c</mi> <munderover> <mi>&amp;Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>n</mi> </munderover> <msub> <mi>z</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> </mrow>

式中，令质心为则

d_i＝|aΔx_i+bΔy_i+cΔz_i|

再对中等号两边对于a、b、c求偏导数，得到

<mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mn>2</mn> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>n</mi> </munderover> <mrow> <mo>(</mo> <msub> <mi>a&amp;Delta;x</mi> <mi>i</mi> </msub> <mo>+</mo> <msub> <mi>b&amp;Delta;y</mi> <mi>i</mi> </msub> <mo>+</mo> <msub> <mi>c&amp;Delta;z</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <msub> <mi>&amp;Delta;x</mi> <mi>i</mi> </msub> <mo>-</mo> <mn>2</mn> <mi>&amp;lambda;</mi> <mi>a</mi> <mo>=</mo> <mn>0</mn> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mn>2</mn> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>n</mi> </munderover> <mrow> <mo>(</mo> <msub> <mi>a&amp;Delta;x</mi> <mi>i</mi> </msub> <mo>+</mo> <msub> <mi>b&amp;Delta;y</mi> <mi>i</mi> </msub> <mo>+</mo> <msub> <mi>c&amp;Delta;z</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <msub> <mi>&amp;Delta;y</mi> <mi>i</mi> </msub> <mo>-</mo> <mn>2</mn> <mi>&amp;lambda;</mi> <mi>b</mi> <mo>=</mo> <mn>0</mn> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mn>2</mn> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>n</mi> </munderover> <mrow> <mo>(</mo> <msub> <mi>a&amp;Delta;x</mi> <mi>i</mi> </msub> <mo>+</mo> <msub> <mi>b&amp;Delta;y</mi> <mi>i</mi> </msub> <mo>+</mo> <msub> <mi>c&amp;Delta;z</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <msub> <mi>&amp;Delta;z</mi> <mi>i</mi> </msub> <mo>-</mo> <mn>2</mn> <mi>&amp;lambda;</mi> <mi>c</mi> <mo>=</mo> <mn>0</mn> </mrow> </mtd> </mtr> </mtable> </mfenced>

求解平面参数a、b、c，转化为求解矩阵A的特征值和特征向量；矩阵A是一个三节对称矩阵，其特征值的求解公式为：

<mrow> <mi>&amp;lambda;</mi> <mo>=</mo> <mfrac> <mrow> <msup> <mrow> <mo>(</mo> <mi>A</mi> <mi>x</mi> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mo>&amp;CenterDot;</mo> <mi>x</mi> </mrow> <mrow> <msup> <mi>x</mi> <mi>T</mi> </msup> <mo>&amp;CenterDot;</mo> <mi>x</mi> </mrow> </mfrac> <mo>,</mo> <mrow> <mo>(</mo> <mi>x</mi> <mo>&amp;NotEqual;</mo> <mn>0</mn> <mo>)</mo> </mrow> </mrow>

在a²+b²+c²＝1的约束条件下，求得利用质心求得d；

步骤7.3：调整法线方向

判断所有法线方向是否一致，假如实际视点为V_p，点云中所有任一法线n_i都指向该视点，则认为所有法线方向一致，即满足：

n_i·VP＞0

式中，VP是视点V_p到查询点P_i的向量；如果计算结果不满足式n_i·VP＞0的要求，说明点P_i的法向量n_i与点P_i到视点V_p的向量之间的向量夹角大与90°，点P_i的法向量n_i应该反向；根据此点云配准方法进行接触网三维重建过程中点云的配准，获得接触网全网三维模型，并对其中缺损的绝缘子进行故障检测。