CN109919955A - 地基式激光雷达点云的隧道轴线提取和分割方法 - Google Patents

Abstract

本发明提出的是一种地基式激光雷达点云的隧道轴线提取和分割方法，包括如下步骤：（1）基于2D投影法提取隧道轴线点；（2）构建轴线MST；（3）基于SG滤波进行隧道轴线平滑；（4）基于改进RANSAC算法进行轴线分割；（5）总体最小二乘轴线拟合。本发明采取先分割后拟合的策略，实现隧道轴线方程的分段表达；降低分割的复杂度，实现对隧道要素的精确分割；使用总体最小二乘方法拟合隧道轴线，降低对噪声的敏感程度。

Description

地基式激光雷达点云的隧道轴线提取和分割方法

技术领域

本发明涉及的是一种地基式激光雷达点云的隧道轴线提取和分割方法，属于隧道断面测量技术领域。

背景技术

使用地基式三维激光扫描对地铁隧道进行变形监测是当前地铁隧道监测的一种趋势，与传统隧道变形监测方式相比，地基式三维激光具有无接触、无损害、效率高、精度高、全天时、全天候等优点。现有技术针对利用地基式三维激光进行隧道变形监测的方法进行了大量研究，现有检测方法中，应用最为广泛的方法是提取隧道断面，然后通过对比不同时期的隧道断面来分析隧道的健康状况。在提取隧道断面时，最关键的步骤就是提取隧道轴线来确定隧道在空间中的姿态，以保证提取的隧道断面和隧道走向严格正交。

隧道轴线提取方法较多，如文献“托雷. 基于三维激光扫描数据的地铁隧道变形监测[D]. 中国地质大学(北京).”提出了一种“基于投影提取隧道轴线”的方法，该方法将隧道点云投影到XOY面、XOZ面内，在两个二维平面内，提取隧道的边界，然后基于边界提取投影后的隧道轴线方程，将两个二维轴线方程联合即得到隧道轴线在空间中的方程。这方法可以快速提取隧道轴线点，但该方法对于非线性隧道、线性和非线性混合隧道无法通过拟合获取完整的隧道轴线方程。

文献“Xie X , Lu X . Development of a 3D Modeling Algorithm for TunnelDeformation Monitoring Based on Terrestrial Laser Scanning[J]. UndergroundSpace, 2017:S2467967417300120.”中也提出了“基于隧道圆柱面拟合的方法”来提取隧道轴线。该方法事先将隧道人工进行分段，然后对分割后的隧道点云使用圆柱面进行拟合，拟合得到的圆柱面的轴线作为隧道轴线。由于采用的是逐段拟合的方法，所以该方法提取的隧道轴线不连续。此外，由于该方法采用的是圆柱面拟合求取隧道轴线，所以对于非线性隧道、线性和非线性隧道连接处，轴线提取结果偏差较大。综上所述，当前存在的隧道轴线提取方法尚不能精确提取、非线性轴线、线性和非线性混合隧道轴线。此外，传统方法尚不能完全克服隧道内外噪声点、离值点对轴线提取的干扰。

发明内容

本发明的目的在于解决现有各种隧道轴线提取方法存在的上述缺陷，提出一种从地基式三维隧道点云中提取隧道轴线，并对隧道轴线进行分割的方法，采取先分割后拟合的策略，实现隧道轴线方程的分段表达；降低分割的复杂度，实现对隧道要素的精确分割；使用总体最小二乘方法拟合隧道轴线，降低对噪声的敏感程度。

本发明的技术解决方案：地基式激光雷达点云的隧道轴线提取和分割方法，包括如下步骤：

（1）基于2D投影法提取隧道轴线点；

（2）构建轴线MST（Minimum Spanning Tree，最小生成树）；MST的构建首先将轴线点以邻接多重表的方式进行存储，并生成无向图，然后以距离最为权重，对无向图进行最大深度遍历获取轴线的MST；

（3）基于SG滤波的隧道轴线平滑；SG（Savitzky-Golay）滤波是一种基于局部多项式拟合的滤波方法。该滤波方法从轴线点的第i个点开始，创建一个的一个窗口，对该窗口的轴线点使用多项式方程进行拟合，然后插值得到第i个点的坐标。这样从第i点起，对轴线点逐个进行上述处理得到最终滤波后的点；

（4）基于改进RANSAC（Random Sample Consensus，随机抽样一致）算法的轴线分割；在轴线分割时如果使用经典RANSAC算法，会出现分割不准确或者过分割的情况，因此需要对RASNSAC算法进行改进。经典RANSAC算法实现步骤如下：

①首先输入数据，在该数据中随机抽取若干个点，求解这若干个点所属的参数模型；

②使用求解出的参数模型测试所有的输入数据，如果某个点适用于估计的模型，认为该点为局内点；

③如果如果有足够多的点被归类为假设的局内点，那么估计的模型就足够合理；

④然后，用所有假设的局内点去重新估计模型；

⑤最后，通过估计局内点与模型的错误率来评估模型：这个过程被重复执行固定的次数，每次产生的模型要么因为局内点太少而被舍弃，要么因为比现有的模型更好而被选用；

（5）总体最小二乘轴线拟合。

本发明的优点：

（1）采用先分割后拟合的策略，提取连续完整的隧道轴线；将轴线方程用分段方程进行表达，在解决隧道轴线连续性的同时提升了轴线拟合的精度；

（2）将轴线数据转化至正切空间，将数据空间中的线性和非线性分割问题变为正切空间中的线性识别问题，同时结合改进的RANSAC算法识别线性元素，增加了轴线分割的准确性；

（3）采用了总体最下二乘方法，确保轴线拟合估计更加精确，同时克服隧道轴线的噪声的影响，在轴线包含噪声的情况下依然能够准确拟合隧道轴线方程。

附图说明

附图1是本发明地基式激光雷达点云的隧道轴线提取和分割方法的总体步骤流程图。

附图2是隧道轴线点提取投影示意图，其中（a）是XOY面投影，（b）是XOZ面投影。

附图3是轴线点恢复拓扑示意图，其中（a）是轴线点构建MST示意图，（b）是最大深度遍历结果示意图。

附图4是SG算法滤波前和滤波后的轴线点在正切空间线性分布比较图，其中（a）是无噪声的轴线点转化至正切空间的轴线点分布图，（b）是增加噪声后的轴线点转化正切空间的轴线点分布图，（c）是滤波后的轴线点转化至正切空间的轴线点分布图。

具体实施方式

隧道轴线反映了隧道点云在空间中分布，也同时反映了局部范围内隧道的姿态和走向，依据轴线上指定的某一点，可方便提取出在该点垂直于隧道轴线的隧道断面。本发明将三维空间轴线表达问题转化为对应两个二维投影平面的共同表达问题。为此首先判断隧道延伸方向，并将点云投影到两个垂直的二维平面内。例如，若隧道延X轴方向延伸，则把隧道三维点云分别投影到XOY面和XOZ面；若延Y轴方向延伸则把隧道三维点云分别投影到XOY，YOZ面，然后在两个二维影平面内分别提取隧道轴线，最后联合这两个投影平面内的轴线方程，共同表达隧道三维轴线方程。本发明以XOY面为例介绍隧道轴线点的提取方法。

如图1所示，本发明提出的地基式激光雷达点云的隧道轴线提取和分割方法，具体包括：输入隧道点云数据，根据输入的点云判断隧道的大致走向，然后依据隧道走向将点云投影到两个垂直的平面内，例如XOY面和XOZ面内。在两个平面内分别进行以下步骤提取隧道轴线：首先使用切片法提取隧道轴线点，对提取的轴线点生成无向图并使用欧式距离作为权重，进行最大深度遍历，构建轴线点的最小生成树MST，在此基础上使用SG滤波算法对轴线点进行滤波，然后把滤波后的轴线点转化至正切空间，在正切空间内进行线性识别。最后把正切空间识别结果映射到数据空间，并对分割结果采用多项式进行拟合以获取最终的隧道轴线。

地基式激光雷达点云的隧道轴线提取和分割方法的具体操作步骤如下：

（1）基于2D投影法提取隧道轴线点

为求取XOY面轴线，在XOY面内，以横坐标最小值为起点，以给定阈值为间隔，在X轴的铅垂线方向进行竖直切割，直到达到点云横坐标的最大值为止；

。在该切片内搜索点云纵坐标的最大值和最小，从而计算出对应的中点坐标，如图2所示，依据上述过程，求取所有切片的中点集合。类似的方法可以用于求取另一个投影平面XOZ面内隧道的轴线方程。

（2）构建轴线MST

由于步骤（1）中所述提取的投影后的隧道轴线是一系列无序的点，为了后续处理方便需要恢复这些点的拓扑关系，即点序的顺序或者逆序拓扑关系。为此，我们首先对这些无序点生成无向图，而后使用点间的欧式距离度量作为权重遍历最大深度子树得到组成隧道轴线点的MST以获取轴线点顺序和逆序拓扑关系。如图3所示，图3中的（a）是隧道轴线点生成无向图，从图中可以看出轴线上包含少量的噪声。图3中的（b）表示了遍历最大深度子树后的结果。从图中可以看出，基于隧道轴线MST的最大深度子树遍历除了能够恢复这些点序列的拓扑外，MST还具有滤除部分初始隧道轴线点噪声的作用，起到了对初始轴线点进一步优化的效果。

（3）基于SG滤波的隧道轴线平滑

构建轴线MST后虽然可以去除一部分轴线噪声，但是轴线仍包含大量噪声，所以从步骤（2）中获得的具有拓扑关系的轴线数据需要进一步地平滑来抑制噪声。SG滤波算法是一种理想的滤波方法，它可以在去除噪声的同时保留轴线的细节特征。经过SG滤波后，轴线的线性要素和非线性要素特征更加明显，这将有助于后续轴线分割算法快速识别出线性要素和非线性要素。图4比较了滤波前与滤波后轴线点转换至正切空间的分布情况。从图中可以看出，理想情况下轴线的线性元素在正切空间被转化为水平直线，而非线性元素被转化为具有斜率的直线。那么轴线在数据空间的分割的问题就被转化为正切空间的线性识别问题。

（4）基于改进RANSAC算法的轴线分割

对于经过滤波后的轴线点，本发明使用改进的RANSAC算法对其实现了线性要素和非线性要素的分割。经典的RANSAC算法直接对隧道轴线进行分割，会出现过分割的情况。例如，对于曲率较小的轴线段将错误的分割为多个直线段。为此，本发明提出了使用改进的RANSAC算法来解决上述问题。改进的RANSAC算法步骤如下：

①把轴线点转化至正切空间，其中为当前点与起始点之间的路径累计长度，为当前点的切线方向与横坐标的正方向之间的夹角；

②使用RANSAC算法拟合面内的线性要素，此时轴线点被分割为局内点和局外点，求出局内点中点作为拟合直线的位置属性；

③确定局内点的坐标范围，若局外点在局内点范围内，则把该点从局外点中删除；

④剩余的局外点作为输入，循环步骤②和③，直到处理完所有的局外点；

⑤根据正切空间中提取的线性要素的位置属性从左到右对线性元素进行排列，然后求出相邻两个线性元素的交点，即为分割位置；

⑥把交点映射到数据空间，对轴线进行分割。

（5）总体最小二乘轴线拟合

对提取的隧道轴线经过构建MST、SG滤波、轴线分割处理后，需要对分割后的轴线使用多项式方程进行拟合，以获取轴线的方程。本发明采用总体最小二乘方法拟合多项式方程，这样使得估计方法更为合理，也保证了即使轴线在包含噪声的情况下依然可以拟合出隧道轴线。

一般多项式方程为：

（1）

令，其中变量x为观测值，那么上式非线性方程可以变化成多元线性方程：

（2）

假设有m个曲线上的点坐标，这m个点坐标构成如下方程组：

将其写成矩阵的形式为：

（3）

其中，B=，。

Y和B均为观测值，为了顾及到两者的误差，采用总体最小二乘法进行平差，式（3）的误差方程形式为：

（4）

其中，，分别为误差矩阵。数学上通常采用矩阵的奇异值分解（SVD）来求参数的总体最小二乘解。

本发明基于隧道点云数据，提出了基于切片法提取隧道轴线，并使用改进的RANSAC算法对隧道轴线分割，对提取既有直线又曲线的隧道轴线具有很好的适用性。采用本发明所提出的方法可以对隧道的曲直线段进行准确的分割，然后获取连续的隧道轴线方程。本发明把数据空间中的曲直线分割问题转化为正切空间的线性识别问题，增加识别的准确性的同时也简化了计算过程，使用总体最小二乘方法拟合隧道轴线，进一步增加了拟合准确性。

Claims (6)

1.地基式激光雷达点云的隧道轴线提取和分割方法，其特征是包括如下步骤：

（1）基于2D投影法提取隧道轴线点；

（2）构建轴线MST；

（3）基于SG滤波的隧道轴线平滑；

（4）基于改进RANSAC算法的轴线分割；

（5）总体最小二乘轴线拟合。

2.根据权利要求1所述的地基式激光雷达点云的隧道轴线提取和分割方法，其特征是所述步骤（1）基于2D投影法提取隧道轴线点具体包括：

输入隧道点云数据，然后将点云投影到互相垂直的XOY面和XOZ面内；在XOY面内，以横坐标最小值为起点，以给定阈值为间隔，在X轴的铅垂线方向进行竖直切割，直到达到点云横坐标的最大值为止；

，在该切片内搜索点云纵坐标的最大值和最小，从而计算出对应的中点坐标，求取所有切片的中点集合，即隧道轴线点；由此类推，求取另一个投影平面XOZ面内隧道的轴线方程。

3.根据权利要求1所述的地基式激光雷达点云的隧道轴线提取和分割方法，其特征是所述步骤（2）构建轴线MST具体包括：

首先对步骤（1）所得到的无序的隧道轴线点生成无向图，而后使用点间的欧式距离度量作为权重遍历最大深度子树，得到组成隧道轴线点的MST，以获取轴线点顺序和逆序拓扑关系数据。

4.根据权利要求1所述的地基式激光雷达点云的隧道轴线提取和分割方法，其特征是所述步骤（3）基于SG滤波的隧道轴线平滑具体包括：

将步骤（2）中获得的轴线点顺序和逆序拓扑关系数据经过SG滤波算法处理，将轴线的线性元素在正切空间转化为水平直线，非线性元素转化为具有斜率的直线。从而将轴线在数据空间的分割的问题转化为正切空间的线性识别问题。

5.根据权利要求1所述的地基式激光雷达点云的隧道轴线提取和分割方法，其特征是所述步骤（4）基于改进RANSAC算法的轴线分割具体包括：

①把轴线点转化至正切空间，其中为当前点与起始点之间的路径累计长度，为当前点的切线方向与横坐标的正方向之间的夹角；

②使用RANSAC算法拟合面内的线性要素，此时轴线点被分割为局内点和局外点，求出局内点中点作为拟合直线的位置属性；

③确定局内点的坐标范围，若局外点在局内点范围内，则把该点从局外点中删除；

④剩余的局外点作为输入，循环步骤②和③，直到处理完所有的局外点；

⑤根据正切空间中提取的线性要素的位置属性从左到右对线性元素进行排列，然后求出相邻两个线性元素的交点，即为分割位置；

⑥把交点映射到数据空间，对轴线进行分割。

6.根据权利要求1所述的地基式激光雷达点云的隧道轴线提取和分割方法，其特征是所述步骤（5）总体最小二乘轴线拟合具体包括：

一般多项式方程为：

（1）

令，其中变量x为观测值，将上式非线性方程变化成多元线性方程：

（2）

假设有m个曲线上的点坐标，这m个点坐标构成如下方程组：

将其写成矩阵的形式为：

（3）

其中，B=，。

Y和B均为观测值，为了顾及到两者的误差，采用总体最小二乘法进行平差，式（3）的误差方程形式为：

（4）

其中，，分别为误差矩阵；采用矩阵的奇异值分解来求参数的总体最小二乘解。