隧道中轴线的提取对于分析隧道几何信息有重要的作用。中轴线的提取方法主要有截面圆提取法、基于投影的中线提取法、基于法矢的提取方法等。截面圆提取法根据被测物体的长度，按一定的间隔获取若干个截面圆或椭圆，然后将截取的圆心或椭圆的圆心进行直线拟合，得到中轴线，但该方法提取精度受人为因素影响较大，可靠性不足。基于投影的中轴线提取法将点云分别投影至两个互相垂直的平面上，然后提取两个平面的边线点并求取平均值拟合出中轴线^[1-3]，但该方法易受到粗差的影响，而且中轴线方向提取精度低。基于法矢的提取方法以法矢为基础，构建法矢与中轴线间的关系进行中轴线的提取，对圆柱状的点云具有很好的适用性，其中法矢的计算精度对中轴线的提取结果有重要的影响。蓝秋萍等^[4]利用隧道表面点云法矢与中轴线之间的垂直关系，截取不同区段的点云采用高斯球映射的方法拟合出该区段的中轴线，由于受到法矢精度的影响，中轴线提取精度较低。汪子豪^[5]为提高中轴线的提取精度，设定一阈值将法矢拟合精度较低的点滤除，得到可靠性比较高的点，再参与中轴线的拟合，以达到提高中轴线拟合精度的目的。为了得到准确的中轴线结果，本文通过质量控制的方法对隧道点云法矢进行精确计算，并利用圆柱多段拟合方法进行隧道中轴线的提取。

1 隧道中轴线提取方法

图 1给出了隧道中轴线提取流程。

1.1 数据预处理

预处理是数据分析的前提，为得到完整、可靠的隧道几何信息需要对原始的点云数据进行拼接、去噪和按距离采样等处理。处理后的隧道表面点云数据没有明显的噪声点，分布均匀，比较平顺，有利于后期中轴线的提取。

1.2 隧道的分割

对于直线型隧道，由于其直线度较好，可直接进行中轴线的提取。而对于曲线型的隧道数据，由于其走向弯曲不能直接提取中轴线，可将隧道沿其走向分成几个小的直线段，然后对每一段分割区域用圆柱多段拟合方法提取中轴线。将曲线型的隧道分割为一段段的直线段处理(图 2，方框区域即为沿隧道走向分割的直线段)，有利于把复杂的弯曲段隧道中轴线的提取简化为直线型隧道中轴线的提取。

1.3 法矢精确计算

点云表面法矢计算结果存在二值性和分布变化不连续性，而且也受点分布疏密和邻域选用范围等的影响，因此对法矢的精确计算是后期数据质量的重要保证。采取3个步骤对法矢的计算进行质量保证和质量控制^[6]。

1) 质量保证。采样的目的是使点云分布均匀，采样的意义和效果见文献[6]。

2) 法矢计算。在每个点的邻域范围内选择k个点，采用主成分分析的方法得到局部拟合平面的参数值a₀、b₀、c₀、d₀，进而得到该点处的平面法矢。误差方程如下：

$ \mathit{\boldsymbol{V = A\hat X}}-\mathit{\boldsymbol{Z}} $

(1)

式中，$ \mathit{\boldsymbol{A = }}\left[{\begin{array}{*{20}{c}} {{x_1}}&{{y_1}}&1\\ {{x_2}}&{{y_2}}&1\\ \vdots & \vdots & \vdots \\ {{x_k}}&{{y_k}}&1 \end{array}} \right], \mathit{\boldsymbol{\hat X = }}\left[{\hat a\;\;\;\;\hat b\;\;\;\;\;\hat d} \right], \mathit{\boldsymbol{Z = }}{\left[{{z_1}\;\;\;{z_2}\; \cdots \;{z_k}} \right]^{\rm{T}}}$。拟合初值${{\mathit{\boldsymbol{\hat X}}}_0} = \left[{\frac{{{a_0}}}{{{c_0}}}\;\;\frac{{{b_0}}}{{{c_0}}}\;\;\frac{{{d_0}}}{{{c_0}}}} \right] $。

3) 质量控制。利用邻域传递方式解决法矢计算结果的二值性问题，根据各点的邻域法矢分布，采用双边滤波器对法矢分布及变化的连续性进行调整，构造双边滤波器^[7-8]：

$ {n_p} = \frac{{\sum\limits_{i \in {\rm{KNN}}} {n{p_i} \cdot \left( {\left\| {{P_i}-P} \right\|} \right){G_P}\left( {\left\| {{{\tilde P}_i}-P} \right\|} \right)} }}{{\sum\limits_{i \in {\rm{KNN}}} {{G_s}\left( {\left\| {{P_i}-P} \right\|} \right){G_P}\left( {\left\| {{{\tilde P}_i} - P} \right\|} \right)} }} $

(2)

其中，G_s(x)和G_P(x)分别为空间距离与预测距离的高斯核函数。

经过以上步骤，可以获取高精度的法矢(图 3)，为中轴线的提取提供高精度的数据保证。

1.4 中轴线粗提取

隧道中某一分割的直线段内的中轴线可以用一个定位点加一个方向矢量来表示。假设被分割的某直线段数据有n个点云数据{(x_i, y_i, z_i)，i=0, 1, 2…n}，隧道中心点(中轴线定位点)坐标为o(x₀, y₀, z₀)= $\sum\limits_{i = 1}^n {\left( {{x_i}, {y_i}, {z_i}} \right)} /n $。将计算完的法矢平移至隧道中心点处并单位化，构成一圆环数据的点集。对这些点集进行最小二乘拟合^[9]，可以得到平面方程aX +bY +cZ +d=0，其中X ={x₁, x₂, …x_n}^T，Y ={y₁, y₂, …y_n}^T，Z ={z₁, z₂, …z_n}^T。将此平面的法矢n₀(a₀, b₀, c₀)作为中轴线方向的初始计算结果。

理论上隧道表面点的法矢{n_i(a_i, b_i, c_i)，i=0, 1, 2…n}与中轴线方向n₀(a₀, b₀, c₀)夹角为90°，但实际中却很难满足这样的条件。因此，为了得到中轴线更优的结果，对于与中轴线夹角跟理论值偏差较大的点(即满足n₀(a₀, b₀, c₀)· n_i(a_i, b_i, c_i)/| n₀|| n_i|－π/2>ε点)予以删除，阈值的确定可以参考文献[5]，再利用上面介绍的方法重新拟合得到较可靠的中轴线数据。

1.5 中轴线精提取

上述方法得到的中轴线精度受到中轴线定位点选取和方向线偏差的影响，只能得到粗略的结果。由于外界压力等影响，隧道会呈现扁率很小的椭圆，可近似视为圆形。为了能够得到精度较高的中轴线，首先在分割的每一个大区段内再进行次小区域的分割，对次小区域进行最小二乘圆柱面拟合，将得到的中轴线的改正值用来修正其位置点和方向线，以得到较精确的值；然后将该区段中m段次小区段分别进行最小二乘圆柱拟合；最后将所求的每个次小区段中轴线的值求平均作为该大区段中轴线的最终值。

在理想状态下隧道表面点到旋转轴的距离应该相等，但事实上存在少许的偏差(d₁≠d₂≠d₃≠d₄), 如图 4所示。为了得到中轴线的最优解，对该直线段中的中轴线位置进行调整，以满足隧道表面各点到中轴线的距离的偏差最小。圆柱面可认为是由到中轴线上一定距离R的点的集合，因此圆柱面确定可以用7个参数来表示，分别为中轴线定位点o(x₀, y₀, z₀)、中轴线方向向量n₀(a₀, b₀, c₀)以及圆柱半径R，测量点与轴线的关系如图 5所示。

图 5中，p_i(x_i, y_i, z_i)为测量点，测量点到中轴线的投影与测量点的连线为R (即圆柱的拟合半径)，o(x₀, y₀, z₀)为中轴线定位点，n₀(a₀, b₀, c₀)为中轴线的方向向量，p′_i为p_i在中轴线上的投影。那么圆柱面的数学表达式为^[10]：

$ \begin{array}{l} R_0^2 = {\left( {{x_i} - {x_0}} \right)^2} + {\left( {{y_i} - {y_0}} \right)^2} + {\left( {{z_i} - {z_0}} \right)^2}\\ - \frac{{{{\left[ {{a_0}\left( {{x_i} - {x_0}} \right) + {b_0}\left( {{y_i} - {y_0}} \right) + {c_0}\left( {{z_i} - {z_0}} \right)} \right]}^2}}}{{a_0^2 + b_0^2 + c_0^2}} = \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;op_i^2 - {{op'_i}^2} \end{array} $

(3)

$ \begin{array}{l} A = {\left( {{x_i}- {x_0}} \right)^2} + {\left( {{y_i}- {y_0}} \right)^2} + {\left( {{z_i}- {z_0}} \right)^2} - \\ \frac{{{{\left[{{a_0}\left( {{x_i}-{x_0}} \right) + {b_0}\left( {{y_i}-{y_0}} \right) + {c_0}\left( {{z_i}-{z_0}} \right)} \right]}^2}}}{{a_0^2 + b_0^2 + c_0^2}} -R_0^2 \end{array} $

(4)

由于式(4)为关于(x₀, y₀, z₀, a₀, b₀, c₀, R₀)的非线性函数，对其进行线性化得到误差方程：

$ \begin{array}{l} v = \frac{{\partial A}}{{\partial {x_0}}}\delta {x_0} + \frac{{\partial A}}{{\partial {y_0}}}\delta {y_0} + \frac{{\partial A}}{{\partial {z_0}}}\delta {z_0} + \frac{{\partial A}}{{\partial {a_0}}}\delta {a_0} + \\ \frac{{\partial A}}{{\partial {b_0}}}\delta {b_0} + \frac{{\partial A}}{{\partial {c_0}}}\delta {c_0} + \frac{{\partial A}}{{\partial R}}\delta {R_0}-{A_0} \end{array} $

(5)

式中，(x₀, y₀, z₀, a₀, b₀, c₀)初值可由§1.4得到，R₀初值可通过将(x₀, y₀, z₀, a₀, b₀, c₀)代入式(3)得到。将式(5)扩展为整个点集，则有：

$ \begin{array}{l} {\mathit{\boldsymbol{V}}_S} = \left[{\frac{{\partial S}}{{\partial {x_0}}}\;\;\frac{{\partial S}}{{\partial {y_0}}}\;\;\;\frac{{\partial S}}{{\partial {z_0}}}\;\;\;\frac{{\partial S}}{{\partial {a_0}}}\;\;\;\frac{{\partial S}}{{\partial {b_0}}}\;\;\frac{{\partial S}}{{\partial {c_0}}}\;\;\;\frac{{\partial S}}{{\partial {R_0}}}} \right] \cdot \\ {\left[{\delta {x_0}\;\;\;\delta {y_0}\;\;\;\delta {z_0}\;\;\;\delta {a_0}\;\;\;\delta {b_0}\;\;\;\delta {c_0}\;\;\;\delta {R_0}} \right]^{\rm{T}}} -{\mathit{\boldsymbol{S}}_0} = \mathit{\boldsymbol{B\hat X}} -\mathit{\boldsymbol{l}} \end{array} $

(6)

依据间接平差理论解算出(δx₀, δy₀, δz₀, δa₀, δb₀, δc₀, δR₀)，对中轴线的初值进行改正，进而得到更为准确的提取结果。为了得到更稳定的结果，可以将式(6)计算出的结果作为初值再次代入式(6)中解算，直至两次提取出的旋转轴的方向矢量夹角小于设定的阈值。最终将m个区段所得的中轴线的值求平均，即(x₀, y₀, z₀, a₀, b₀, c₀)= $ \left( {\frac{{\sum {{x_0}} }}{m}, \frac{{\sum {{y_0}} }}{m}, \frac{{\sum {{z_0}} }}{m}, \frac{{\sum {{a_0}} }}{m}, \frac{{\sum {{b_0}} }}{m}, \frac{{\sum {{c_0}} }}{m}} \right)$为最终结果。

2 实验与分析 2.1 最适采样间隔分析

采集的隧道点云具有近密远疏的特征，点云分布不均对中轴线的提取有一定的影响，点云密度过大，数据冗余度高，中轴线不稳定；点云密度过小，隧道的一些特征会被忽略，中轴线提取不准确。为得到中轴线稳定的结果，需要对点云进行采样，找到合适的采样区间。本文采用Regil VZ-400三维激光扫描仪以较高分辨率获取的一段隧道实测数据进行实验。截取一段点云，并进行不同间隔的采样，表 1为采样间隔从4~100 mm、由密到疏提取的27组隧道中轴线的位置和方向数据，中轴线与各坐标轴之间的夹角(弧度)如图 6所示。从图 6可以看出，中轴线在采样间隔为0.026~0.035 mm时与各坐标轴间保持较稳定的角度值，定位点的坐标值也在亚mm的范围内微小变动(图 6(d)中实心点为采样间隔为0.026~0.035 mm的定位点)，说明本方法能够得到较为稳定的结果。对于隧道这样的大型构筑物，点云之间太密，数据冗余度高，会增加计算量；点云之间太疏，又不能真实地反映构筑物的表面信息。图 6、表 1信息有助于选择合适的采样间隔，得到稳定的提取结果。

为验证本文方法是否满足重复提取精度的要求，将采样间隔为0.026~0.035 mm的数据求取标准差。中轴线方向与坐标轴夹角的标准差为1.202″，定位点的标准差为0.024 mm，能够满足工程需要。因此，点云数据采样间隔在0.026~0.035 mm之间，中轴线提取的精度能够满足重复提取的要求且稳定性好。

2.2 多段拟合分析

隧道实际施工跟理论设计间存在稍许偏差(在限差内)。为了得到隧道在直线段和曲线段的精确中轴线数据，沿隧道走向将隧道分成若干段，然后对每一段拟合获取精确的中轴线信息，具体步骤见上文，这样便可以将长隧道进行分段化处理，保证结果的可靠性。实验所用直线段数据约为6 m，分为两段，曲线段数据约15 m，分为4段，隧道弯曲程度不大，中轴线粗提取结果和精提取结果如图 7所示。

表 2为中轴线精提取数据。为了更好地对隧道走向的变化进行定量分析，将隧道中轴线与z坐标轴的夹角用角度值表示。隧道直线段两段间的轴向变化约6′，在分段拟合后可以发现轴向的微小变化，但此微小变化在工程建设的误差要求范围内，可以忽略。从图 7(d)中可以明显看出，第1段与第2段存在稍许的弯曲变化。从表 2定量分析可得，两段轴向偏差达到近1°56′，而后3段则变化不大，这与隧道真实的走向相一致。由此可知，分段拟合方法可以用于隧道直线段和曲线段中轴线的提取工作，精度满足隧道中轴线提取的要求。

3 结语

本文以隧道的三维点云数据为研究对象，提出了一种基于圆柱多段拟合提取隧道中轴线的方法。首先对隧道数据进行拼接、去噪、采样；然后为方便数据处理，将隧道截取为多个小区段；再通过质量控制方法获取高质量的法矢数据进而得到中轴线初值；最后基于圆柱多段拟合提取精度较高的中轴线结果。由实验分析可知，隧道最适采样间隔在0.026~0.035 mm之间，此方法不仅在隧道直线段有良好的适用性，在隧道弯曲处也有同样的适用性，为隧道中轴线的提取提供了一种可参考的方案。但该方法对于大数据量的隧道点云处理还有一定局限性。