在地铁施工运营的过程中，地铁隧道的结构安全极其重要，地铁隧道结构竣工测量和安全监测都需要地铁隧道断面测量，地铁隧道断面监测是地铁隧道结构安全监测中最重要的部分，近年来主要使用测量机器人和三维激光扫描仪进行监测，然而测量机器人只能反映单点的位移变形情况，不能反映整个断面的变形特征，而传统的三维激光扫描仪单点测量精度不高^[1]，并且需要进行点云配准拼接，徕卡MS50是全站扫描仪，结合了测量机器人的高精度点测量与海量点云数据，并且点云坐标是全站仪坐标系统，无需安置标靶，可以反映整体的隧道断面变形情况^[2]。本文使用MS50对广州某地铁隧道进行扫描，隧道断面为盾构机开挖形成的圆断面，利用直接点云分析方法提取隧道断面中轴线与隧道断面，椭圆拟合之后与测量机器人数据进行比较，结果表明本文提出的方法效果良好^[3]。

1 基于点云的隧道断面提取方法

MS50采集的点云位置信息是全站仪坐标，仪器整平后才能作业，补偿器打开后，Z轴方向为严格天顶方向，指向仪器的竖向扫描面，X、Y为平面坐标。为方便数据处理，在野外扫描时，使Y方向指向隧道前进方向，X方向为与ZOY平面构成左手坐标系的方向。

1.1 隧道中轴线提取

为了保证截取的断面是垂直于隧道的前进方向，所以在截取隧道断面前需要知道截取断面处的隧道走向及姿态。隧道中轴线是用来表达隧道的姿态与走势的虚拟的空间曲线，可得中轴线上任意一点的切线，垂直于切线的平面即为隧道断面所在平面，故必须提取中轴线。本文采用的隧道中轴线拟合方法如下：

1) 将隧道点云数据投影至XOY平面与ZOY平面上，从Y坐标的最小值开始到最大值，以δy为固定间隔进行取值，沿着Y坐标方向分别在XOY与ZOY平面上搜索X坐标与Z坐标的最大值与最小值。

2) 计算任意一点y_i(i=1, 2, …, n)对应于X坐标与Z坐标的均值：

$ \left\{ \begin{align} &{{x}_{i}}=(\max {{x}_{i}}+\min {{x}_{i}})/2 \\ &{{y}_{i}}=(\max {{z}_{i}}+\min {{z}_{i}})/2 \\ \end{align} \right. $

(1)

式中，min x_i与max x_i为对应于y_i的X坐标的最小值与最大值；min z_i与max z_i为对应于y_i的Z坐标的最小值与最大值。理论上min x_i与max x_i的差值为隧道宽度，min z_i与max z_i的差值为隧道高度^[4]。

3) 把y_i分别和求出的x_i、y_i进行三次多项式拟合，有些y_i对应的max y_i和min y_i可能有粗差，两者大小很接近，完全达不到隧道的最大宽度，为了剔除粗差点，本文使用RANSAC(random sample consensus)算法进行参数估计^[5]，得到的隧道中轴线表示为：

$ \left\{ \begin{align} &x=f(y)={{a}_{0}}{{y}^{3}}+{{a}_{1}}{{y}^{2}}+{{a}_{2}}y+{{a}_{3}} \\ &z=g(y)={{b}_{0}}{{y}^{3}}+{{b}_{1}}{{y}^{2}}+{{b}_{2}}y+{{b}_{3}} \\ \end{align} \right. $

(2)

得到的两条中轴线XOY平面内的即是隧道平曲线，YOZ平面内的即是隧道竖曲线。

1.2 隧道断面生成

通过空间曲线方程可以求出方程上任意一点的切线向量，此切向量即为对应于这点的断面法向量，这样就可以得到这个点的断面平面方程。断面的中轴线方程可以表示为：

$ x=f(y), y=y, z=g(y) $

(3)

用空间曲线参数方程求导，中轴线上的某点为M(x₀, y₀, z₀)，切线向量为$ \overset{\to }{\mathop{T}}\, =({f}'({{y}_{0}}), 1, {g}'({{y}_{0}})) $，法平面方程即M点对应的隧道断面方程为:

$ {f}'({{y}_{0}})(x-{{x}_{0}})+(y-{{y}_{0}})+{g}'({{y}_{0}})(z-{{z}_{0}})=0 $

(4)

1.3 隧道断面附近点集投影至断面平面

虽得到隧道的断面平面方程，但点云数据不可能都在平面内，因此需确定一邻域，邻域内的离散点构成的点集内的所有点视为断面上的点，并把这些点投影至平面，本文中若点云中的点到平面的距离小于d，则判断此点在点集内。阈值d的选取根据中轴线的曲率决定，中轴线在该点曲率越大，则阈值d越大，反之越小。如图 1所示，中轴线上任意一点M(x₀，y₀，z₀)，MV是其切线，L是其法平面，距离法平面为d的两个平面内的点即是要进行投影的点集^[6]。

点的投影坐标通过下面的方法求解：已知待投影区域一点的坐标为Q_i(x_i，y_i，z_i)，投影到隧道断面平面内的坐标为Q_i′(x_i′，y_i′，z_i′)，联立方程组进行求解^[7]：

$ \left\{ \begin{align} &{f}'({{y}_{0}})(x-{{x}_{0}})+(y-{{y}_{0}})+{g}'({{y}_{0}})(z-{{z}_{0}})=0 \\ &{{{{x}'}}_{i}}={{x}_{i}}+k{f}'({{y}_{0}}), {{{{y}'}}_{i}}={{y}_{i}}+k, {{{{z}'}}_{i}}={{z}_{i}}+k{g}'({{y}_{0}}) \\ \end{align} \right. $

(5)

可得$ k=[{f}'({{y}_{0}})({{x}_{0}}-{{x}_{i}})+({{y}_{0}}-{{y}_{i}})+{g}'({{y}_{0}})({{z}_{0}}-{{z}_{i}})]/[{f}'{{({{y}_{0}})}^{2}}+{{1}^{2}}+{g}'{{({{y}_{0}})}^{2]}} $, 从而可以求出Q_i′(x_i′，y_i′，z_i′)。这样就求出了中轴线上某一点M(x₀，y₀，z₀)附近邻域的断面点云集。

1.4 隧道断面空间平面旋转

由于隧道断面所在平面为空间平面，而曲线拟合均在二维平面进行，故拟合前需将空间平面上所有点旋转至二维平面，本文中把断面所在平面的法向量旋转至垂直于坐标轴Y轴，这样椭圆所在平面就为平行于XOZ的平面。

需经两次旋转，第一次旋转使法向量X值变为零，第二次使法向量Z值变为零，平面法向量由$ [{f}'({{y}_{0}}), 1, {g}'({{y}_{0}})] $变为$ [0, \sqrt{{f}'{{({{y}_{0}})}^{2}}+1+{g}'{{({{y}_{0}})}^{2}}}, 0] $，根据空间7参数旋转变换理论^{[8, 9]}，可得旋转矩阵为：

$ {{R}_{0}}=\left[ \begin{matrix} 1/\sqrt{{{a}^{2}}+1} & -a/\sqrt{{{a}^{2}}+1} & 0 \\ a/\sqrt{{{a}^{2}}+1+{{b}^{2}}} & 1/\sqrt{{{a}^{2}}+1+{{b}^{2}}} & b/\sqrt{{{a}^{2}}+1+{{b}^{2}}} \\ -ab/(\sqrt{{{a}^{2}}+1}\sqrt{{{a}^{2}}+1+{{b}^{2}}}) & -b/(\sqrt{{{a}^{2}}+1}\sqrt{{{a}^{2}}+1+{{b}^{2}}}) & \sqrt{{{a}^{2}}+1}/\sqrt{{{a}^{2}}+1+{{b}^{2}}} \\ \end{matrix} \right.\left. \begin{align} & \\ & \\ & \\ \end{align} \right] $

(6)

式中，α=f′(y₀)；b=g′(y₀)。旋转为刚性旋转，点与点之间的位置关系不变。

1.5 隧道断面轮廓线拟合

由于地铁隧道断面形状大部分都是标准圆，受压之后变形成椭圆，故本文研究椭圆拟合断面轮廓，即把点云数据以椭圆为模型进行拟合，计算椭圆方程的各个参数。

椭圆属于圆锥曲线，圆锥曲线的平面二次曲线方程的一般形式为：

$ \begin{align} & \ \ \ \ \ \ \ \ \ \ \ F(\mathit{\boldsymbol{a}}, \mathit{\boldsymbol{x}})=\mathit{\boldsymbol{a}}\centerdot \mathit{\boldsymbol{x}}= \\ & A{{x}^{2}}+Bxy+c{{y}^{2}}+Dx+Ey+1=0 \\ \end{align} $

在椭圆拟合中，本文采用直接最小二乘法进行拟合^[10]。由椭圆一般方程和椭圆标准方程的参数对应关系，可以计算椭圆参数(长半轴α、短半轴b、离心率e)分别如下：

$ a=\sqrt{\frac{2(A{{X}_{c}}^{2}+C{{Y}_{c}}^{2}+B{{X}_{c}}{{Y}_{c}}-1)}{A+C+\sqrt{{{(A-C)}^{2}}+{{B}^{2}}}}} $

(7)

$ b=\sqrt{\frac{2(A{{X}_{c}}^{2}+C{{Y}_{c}}^{2}+B{{X}_{c}}{{Y}_{c}}-1)}{A+C-\sqrt{{{(A-C)}^{2}}+{{B}^{2}}}}} $

(8)

$ e=c/a=\sqrt{{{a}^{2}}-{{b}^{2}}}/a $

(9)

在分析地铁隧道一个断面的形变中，主要提取的是椭圆几何中心与长短半轴、离心率。这样，根据隧道内的任意里程即Y坐标，可以计算该里程对应的隧道断面平面以及分析断面变形情况。

2 工程应用案例 2.1 点云数据采集和预处理

本次实验中隧道内径D=2.7×2=5.4 m，隧道未铺轨。测站间距S=1.5D=8.1 m，MS50基本架设于隧道中线，测站间距大概为9 m，一共扫描4站，扫描长度约35 m^[11]。隧道原始数据如图 2所示。

野外获取的点云数据需使用Cyclone软件进行格式转换，使用Geomagic进行点云去噪与点云重采样。经预处理后，截取其中一段点云数据(-8 m＜Y＜8 m)进行分析。

2.2 中轴线生成

以投影到YOX平面的点集求中轴线，同理可得YOZ平面上的点集中轴线。

1) 把点云数据投影至YOX平面，直接取点云的Y坐标和X坐标即可；

2) 在平面点云中求离散的二维中轴线点：①把Y值按一定步长进行等分，本文把Y值从最小到最大分为步长为0.03 m的等差数列；②对每个Y值，找到对应的X值的均值。笔者按照如下方法进行寻找：设当前的Y值为-5 m，设置搜索阈值d为步长的一半，即d=0.015 m，则从-4.985到-5.015区间内的Y值都认为其值为-5 m，在这个区间内，从二维点集(x, y)中遍历每个Y值对应的X值，从这些X值中找到最大与最小值，根据最大与最小值求出对应于当前Y=Y₀的X值的均值; ③剔除偏差大的点。对应于每个Y₀的X值的最大与最小值，加入条件：X_max-X_min>0.8D，否则剔除这组X值数据。

拟合的隧道中轴线方程为：

$ \begin{align} &x=f(y)=7.9354\times {{10}^{-5}}{{y}^{2}}+2.9221\times {{10}^{-4}}{{y}^{2}}-0.0037y-0.0161 \\ &z=g(y)=2.5982\times {{10}^{-5}}{{y}^{2}}-8.6315\times {{10}^{-4}}{{y}^{2}}+0.0037y-0.0162 \\ \end{align} $

(10)

2.3 断面生成与拟合

1) 生成断面所在平面

截取Y=-5 m处的隧道断面，将Y值代入中轴线方程(10)可得隧道中轴线点坐标为M(0.009 8, -5, -0.023 9)，断面所在平面法向量为$ \overset{\to }{\mathop{T}}\, =(-0.0066, 1, 0.0123) $，断面所在平面方程为：

$ -0.0066(x-0.0098)+(y+5)+0.0123(z+0.0239)=0 $

(11)

2) 断面附近邻域点集投影至断面所在平面

Y=-5 m点集投影至该平面。若计算隧道内所有点的距离则耗时较大，故设置一范围-4.8 m＜Y＜-5.2 m，点集Q中所有的点都在此范围内，只需在此范围搜寻。求得的邻域点集Q (部分)见表 1。

求出Q后，按照点到平面的投影公式求出Q在断面所在平面的投影点集Q₁，见表 1(部分)。隧道中轴线、M点对应的隧道断面所在平面、投影到断面平面的点集Q₁如图 3所示。

3) 断面投影点集旋转至二维平面

空间平面点集Q₁左乘旋转矩阵R₀可得点集Q₂，则Q₂=R₀ Q₁，其中∣R₀ ∣=1，旋转后点之间的几何位置关系不变，故旋转后得出的椭圆参数即为原椭圆参数。由上述旋转矩阵公式可得：α=f′(y₀)=-0.006 6, b=g′(y₀)=0.012 3，则旋转矩阵为：

$ {{\mathit{\boldsymbol{R}}}_{0}}=\left[ \begin{matrix} 0.9999&0.0066&0 \\ -0.0066&0.9999&0.0123 \\ 0&-0.0123&0.9999 \\ \end{matrix} \right] $

(12)

可得旋转至二维平面后的点集Q ₂(部分)见表 1。

使用Q ₂中点的X值与Z值进行直接最小二乘椭圆拟合，椭圆拟合效果如图 4所示。拟合的椭圆方程为：

$ -0.13311{{x}^{2}}+0.00089xy-0.13600{{y}^{2}}-\\0.00766x+0.01029+0.98164=0 $

(13)

由上述椭圆参数提取公式可得α=2.716 7 m，b=2.686 3 m，e=0.147 8。

2.4 测量机器人数据采集与分析

在隧道Y=-5 m处安置6个反射片，分别在东、西、东南、西南、西北、西南方向，安置时注意6个反射片严格垂直于中轴线，都在同一个断面平面内。使用MS50对这6个点进行方向观测法观测4个测回，观测完毕后，使用GNPS软件计算坐标，分别为A₁(-2.037，-4.996，1.771) m、A₂(1.700，-4.992，2.099) m、A₃(1.825，-5.012，-1.996) m、A₄ (-1.409，-4.983，-2.292) m、A₅(-2.709，-4.857，0.001) m、A₆(2.726，-5.016，0.002) m。使用这6个断面点采用直接椭圆拟合方法进行拟合。拟合的椭圆方程为：

$ -0.13294{{x}^{2}}-0.00018xy-0.13616{{y}^{2}}+\\0.00198x-0.00004y+0.98173=0 $

由椭圆参数提取公式可得α=2.717 5 m，b=2.685 2 m，e=0.153 8。

该地铁隧道断面设计为一标准圆，其半径设计值为2.7 m，可见测量机器人得到的结果与激光点云得到的结果基本一致，两者数据差别在1 mm左右，说明两者结合可以更好地分析隧道形变，更有效地掌握断面的整体变形信息。

3 结束语

徕卡MS50进行点云扫描时可以高效获取高精度的点云数据，并且无需点云拼接、配准，结合了传统的测量机器人高精度点测量与三维激光扫描仪海量点云扫描的优势，本文研究内容对于全站扫描仪在地铁监测中可提供一定的经验指导。