基于点云数据的隧道中轴线自动提取方法研究

2016-03-16StudyonPointCloudDatabasedAutomaticExtractionMethodforTunnelCentralAxis

重庆建筑 2016年1期

关键词：最小二乘法边界点中轴线

Study on Point Cloud Data-based Automatic Extraction Method for Tunnel Central Axis

李兵，杨勇，吴英勇（西南大学资源环境学院，重庆　400715）

基于点云数据的隧道中轴线自动提取方法研究

Study on Point Cloud Data-based Automatic Extraction Method for Tunnel Central Axis

李兵，杨勇，吴英勇
（西南大学资源环境学院，重庆400715）

摘要：隧道中轴线可以反映隧道的整体走势，而中轴线的准确性影响着隧道点云数据的后续处理。对此，该文提出一种自动提取中轴线的方法，将点云数据投影后利用单位圆极坐标相邻角度差来确定边界点，根据RanSAC算法对不同参数模型进行统计，对直线段、缓和曲线及曲线的曲段进行区分，然后利用RanSAC算法对不同区段拟合，提取各区段的中轴线，对各区段中轴线进行加权整体最小二乘法拟合，最终得到整段中轴线。该文采用上海隧道点云数据进行了中轴线自动提取实验。实验证明了该方法能有效提取中轴线。

关键词：中轴线；边界点；RanSAC算法；最小二乘法

Abstract:The tunnel central axis can reflect the overall trend of the tunnel，and its accuracy affects the processing of the tunnel point cloud data. In this paper, the automatic extraction method for central axis is presented, with which the boundary points can be determined. Different sections are fitted by RanSAC algorithm，and the central axis can be obtained through that of each section. The automatic extraction experiment is done with the point cloud data of Shanghai tunnel, and it shows that this method is proven to be effective.

Keywords:central axis; boundary points; RanSAC algorithm; least square method

引言

隧道中轴线的检测不仅在施工过程中可以控制保证施工仪器沿着设计的轴线掘进，还可以随时提供掘进仪器的位置姿态，保证隧道贯穿精度，如图1所示。再者，在隧道保养维护期间，可以通过中轴线的测量来检测隧道的变形情况。

图1　隧道中轴线偏差（水平）

隧道中轴线是由各自独立的平曲线和竖曲线组合而成。平曲线包括直线、缓和曲线、曲线等3种；竖曲线包括直线、凸曲线、凹曲线等3种[1]。

通过中轴线，很多有关方面的问题都能很好的控制，如人体肠道的自动虚拟导航就是使用肠道的中心线来控制虚拟摄像头的运动和方向的[2]。Ingmar Bitter等人在一种光滑、准确的中轴线提取算法一文中提出通过离散中心线使用差值及逼近曲线方式来创建一个连续的中心线的方法[3]。AlexandruTelea等（2002）提出一种简单的计算任意3维物体的中轴线方法，通过计算水平方向边界像素位置的参数的快速匹配方法[4]。Kimmel等人提出一种基于距离转换的提取中轴线的方法，主要是根据分段边界最大曲率来计算中轴线[5]，然而这种方法的准确性依赖于在存在噪声的边界上准确地检测出曲率。曹斌（2008）[6]和王荣（2011）[7]提出一种计算隧道轴线三维坐标偏差的方法，根据隧道轴线的特点对不同线型进行分析，实现了隧道轴线三维坐标的批量快速计算，为盾构推进过程中偏差的确定以及隧道环片拼装后中心坐标偏移的监测以及环片后期监测提供了一种实现方法。陆步云采用数理统计随机抽样的方法[8]，对隧道轴线数据进行分析、计算从而达到对隧道整体质量的控制。托雷等人[9]（2013）通过投影的方式，利用Ransac算法及最小二乘法提取中轴线，之后王保前[10]（2013）提出利用法向式中轴线提取的方法，以投影式中轴线提取为基础，确定轴线提取的初始位置，其次沿垂直于轴线切线的方向截取隧道断面点，并对断面点截取区域进行局部曲面拟合，确定该断面点坐标及其对应的局部曲面法向量，然后，依据隧道设计的几何结构，该法向量应该处处垂直于隧道轴线切线方向，依此为约束条件调整隧道轴线方向，最后，沿初始轴线方向设定步长依此截取断面，实现隧道中轴线的动态提取。

本文根据隧道中轴线的特点，结合点云数据投影XOY面及YOZ面及利用单位圆极坐标相邻角度差来提取边界点，应用RanSAC算法进行各种参数模型的统计自动识别不同的区段，然后对整个应用RanSAC进行模型粗差点的剔除，拟合模型隧道点云数据的不同区段进行加权整体最小二乘法拟合，最终得到整体的中轴线。

1　基于RanSAC算法的中轴线自动提取

隧道中轴线中平曲线包括直线、缓和曲线、曲线3种，而竖曲线包括直线、凸曲线、凹曲线，但是竖曲线的曲率变化一般不明显，为了实现中轴线的自动提取，首先要对投影的边界进行区段的分类，在此基础上对中轴线提取，最终完成整体的拟合。本文的提取的流程为图2所示。

图2　中轴线提取流程

1.1移动区域法边界点自动提取

将隧道点云数据投影到XOY面及YOZ面上，中轴线可以根据XOY面的两边界的中心线及YOZ面的两边界的中心线得到，因此，点云数据的边界点提取对中轴线有重要作用。隧道是由直线、缓和曲线及曲线组成，而且提取的边界点中存在部分的粗差点，通过RanSAC算法剔除粗差点，根据不同区段的模型自动提取直线、缓和曲线及曲线，最后对重合部位进行加权整体最小二乘法拟合。

本文通过点云数据在XOY面及YOZ面上的投影，根据图3中所示，假设一半径为R的圆，如图3中红色圆所示，以一点P为圆心，提取在圆内的点作为该点的邻近点Pi（i =1，2，...），然后都用极坐标表示，L为极坐标的起始方向。预先设定一个阈值角度β（假设β=150O），如果相邻的极坐标角度差大于阈值β，则P点为边界上的点，否则不是。对所有点分别进行上述方法的处理，则提取出位于边界上的点。

图3　移动区域法提取边界点

提取得到的边界点中会存在部分的粗差点，通过RanSAC算法[11]对边界点进行粗差的剔除。

1.2基于模型参数迭代统计法提取边界线

边界线的直线、缓和曲线及曲线部分分别建立数学模型。

假设直线模型（XOY面）：

缓和曲线模型（XOY面）：

曲线模型（XOY面）：

其中，a,b;c,d,e,f;g,h,k三组分别为直线模型、缓和曲模型及曲线模型的模型参数。

由于边界是由三种模型组合而成，只需要确定出组合模型的最优参数即可确定出唯一的边界曲线。本文通过RanSAC算法对模型参数迭代次数统计的方法，具有统计值最大的模型参数作为模型的最优参数，如图4所示。

图4　模型参数迭代次数统计

对投影后的两边界都应用RanSAC算法迭代方法统计模型的最优参数，然后根据最优参数对边界进行拟合。

1.3模型中轴线的自动提取

边界点拟合后，在边界1上重新等间隔采集点集1，并分别计算出每个点在其法线方向上与边界2的交点（点集2），然后利用点集1、2求得一列中点；再对边界2以同样方法求得一列中点，根据这两列中点的均值拟合曲线，得到最终的中轴线如图5所示。Pl点的法线方向上与边界2的交点为Pl’，计算PI与PI’的中点Ml’；再根据PI’在边界2上的法线方向上与边界I的交点为PI’’，计算PI’与PI’’的中点MI’’，最终以MI’与MI’’的均值作为中轴线上的点MI。逐次计算边界点，将得到最终的中轴线。

图5　中轴线提取

由于模型的不同，此中轴线的提取是两边界模型对应模型之间的提取，又因为RanSAC算法提取边界点时，模型之间提取的边界点可能存在重合的现象（如图6所示），所以需要对重叠但不吻合的部位进行处理。

图6　区段自动提取

1.4基于加权最小二乘法的中轴线整体拟合

第1.1.2节得到中轴线，为了解决重合部位不吻合的现象，本文采用加权最小二乘法整体对不同中轴线之间进行平差拟合[13]，减小重合部位不吻合的影响。方法如下：如果两种不同类型线重合部位吻合的话，他们的对应点坐标值是一样，因此，根据Xtc-Xt=0,Xc-Xtc=0建立模型，设重合部位调整模型为：

其中a,b;c,d,e,f;g,h,k分别为区段i，j的模型参数。

不同模型之间的重合可能是直线与缓和曲线、缓和曲线与曲线、曲线与缓和曲线等方式，在本文中依然是根据不同模型参数的个数来判定重合的区段，然后对采用相应的调整模型，主要是提高重合部位点的权重。

根据间接平差的原理，设：

ii=1,2,…n1;jj=1,2,…n2，kk=1,2,…n3，yii,yjj,ykk重合部位对应点的y坐标值。

B1表示直线模型的隧道中轴线，Bc表示曲线模型的隧道中轴线，Btc表示缓和曲线模型的隧道中轴线；C表示重合部位点的集合；i表示模型个数，m表示模型重采样的个数。

2　实验与分析

本文采用RIEGL-LMS VZ-400对上海地铁隧道进行点云数据的采集（图7），表1为仪器的参数及点云数据的点数。

表1　仪器参数及点数

图7　实验数据（LMS VZ-400）

2.1边界点提取实验

将点云数据投影到XOY面及YOZ面，根据扫描线小角度的移动不断的提取扫描线X坐标值的最大及最小，然后反求得三维坐标，进而提取边界点。图8显示只是部分点云数据。

图8　边界点

2.2边界自动提取

根据投影得到的边界点，RanSAC算法剔除粗差，同时进行模型参数的迭代统计，选择最优的模型参数，作为边界线模型的参数。以一边为例，如图9。

图9　边界线

根据RanSAC的迭代模型参数统计出三种不同的区段的峰值（如图10）。分别为直线段、缓和曲线段及曲线段（图11）。

图10　RanSAC迭代模型参数统计

图11　边界线分区段

2.3中轴线提取

两条边界线通过自动提取方式确定并分别拟合重新等间隔采样，依照第2.1.2节提取的方法对XOY面及YOZ面的中轴线提取。图12为部分点云数据的中轴线提取示意图。

图12　中轴线提取

隧道中轴线的提取是各分段独自提取，不同区段之间的重合性，及模型的不同，必然对出现重合部位不吻合的现象，图13为各区段中轴线提取后的组合在一起的中轴线。白框中是重合部位的局部放大图，分别显示曲线与缓和曲线区段、缓和曲线与直线区段之间的偏差。

图13　中轴线提取

为了消除上述重合部位不吻合的问题，本文采用加权整体最小二乘法，重合部位点的权值根据多次实验表明p=5时，整体拟合的效果较佳。图14为整体拟合后的效果。

图14　提取效果

在重合部位分别随机抽取12对对应点，分别计算点对点之间的偏差，根据偏差来衡量拟合的精度（表2）。

表2　距离精度

表2显示整体拟合后点位偏差变化最大为0.007m，中误差为0.002m；而未拟合前点位偏差最大为0.059m，中误差为0.026m，因此，自动提取中轴线的方法是有效的。

3　结论

本文针对隧道点云数据中轴线提取问题，提出点云数据投影后利用等间距移动扫描线上X轴最大及最小值来求得边界点，通过RanSAC算法剔除粗差，并且结合中轴线的特点，建立不同区段模型，RanSAC算法进行模型参数迭代运算，统计模型参数迭代次数，统计值最大的作为最优模型参数，根据最优模型参数拟合边界线，利用边界点求得中轴线，对区段中轴线之间重合不吻合的问题，提出加权整体最小二乘拟合的方法，最终实现隧道点云数据的中轴线的自动提取。实验结果表明该方法对隧道点云数据中轴线的自动提取有一定的借鉴作用。

参考文献：

[1]徐伟忠，岳秀平.盾构法隧道轴线测量及控制技术[A].大直径隧道与城市轨道交通工程技术——2005上海国际隧道工程研讨会文集[C].2005:272-278.

[2] L.Hong，S.Muraki，A.Kaufman，etc.Virtual Voyage: Interactive Navigation in the Human Colon [C]. Proceedings of Siggraph，1997（27，34）.

[3] Ingmar Bitter，Mie Sato，Bender，M，etc.CEASAR: a smooth，accurate and robust centerline extraction algorithm[C]. Proceedings Visualization 2000：45-52.

[4] Alexandru Telea，Jarke J，van Wijk.An Augmented Fast Marching Method for Computing Skeletons and Centerlines [C]. Joint EUROGRAPHICS-IEEE TCVG Symposium on Visualization（2002）：251-ff.

[5] R.KIMMEL，D.SHAKED，N.KIRYATI，etc.Skeletonization vis Distance Maps and Level Sets [J].Computer Vision and Image Understanding，1995（3）382-391.

[6]曹斌，鲍峰.一种隧道轴线偏差的快速计算方法[A].第八届全国工程地质大会论文集[C].2008:121-123.

[7]王荣，杨雪.计算隧道轴线三维坐标偏差方法的介绍及探讨[J].测绘技术装备，2011（02）:37-39.

[8]陆步云.上海地铁2号线采用数理统计控制区间隧道轴线质量[A].地下工程施工与风险防范技术——2007第三届上海国际隧道工程研讨会文集[C]，2007.

[9]托雷，康志忠，谢远成，等.利用三维点云数据的地铁隧道断面连续截取方法研究[J].武汉大学学报：信息科学版，2013 （02）:171-175+185.

[10]王保前.融合反射值影像和卡尔曼滤波的三维点云全局拼接方法[D].北京：中国地质大学，2013.

[11] FISCHLER MA，BOLLES RC.Random Sample Consensus：A Paradigm for Model Fitting with Applications to Image analysis and automated Cartography [J].Computer Vision，1987：726-740.

[12]贾丰曼，康志忠.基于SIFT算法与贝叶斯抽样一致性检验的影像同名点匹配[J].测绘学报，2013.

[13]侯超钧，曾艳姗，吴东庆，等.全局连续的分段最小二乘曲线拟合方法[J].重庆师范大学学报：自然科学版，2011，（06）:44-48.

责任编辑：孙苏

施工经验

采暖系统施工要点

1干管支管有坡度

管内存有空气会引起系统不热。干管支管均应设置坡度，干管坡向末端，支管坡向暖气片，且有可靠的排气措施，避免永久性空气随水流一起流动，以保证系统正常运行。

2预防回水走近路

对于下供下回系统中的单立管连接，不要一下穿到底，就是说不能用三通分别接暖气片，即上层暖气的回水接下层的供水，对回水走近路引起上层暖气片不热有很好的效果。

3须设置泄水装置

考虑到检修方便，应在系统最低点安装泄水装置，并排入到卫生间下水口或不碍事的地方。

4暖气片须做复验