基于椭圆柱模型的隧道断面拟合算法研究

陈潇，王鑫森，俞雪薇

(天津市勘察院，天津 300191)

1 引 言

隧道断面测量是调线调坡设计和结构限界检测中的一个关键环节，可用于检验施工过程中产生的施工误差以及结构变形等不利因素，作为后续工程的参考依据，因此其测量的精度与数据的准确性对于整条地铁线路的施工建设至关重要。传统的全站仪测量方法采集数据速度慢、效率低，且仅能反映不连续的结构特征，难以作为长距离、大范围的数据采样手段。而地面三维激光扫描技术的出现为解决该问题带来了革命性突破，它通过高速激光测距迅速获取大量的目标对象表面数据点，可提供全面的点云信息，具有非接触性、高精度、高密度等优点，克服了传统测量技术的局限性和片面性[1，2，3，7]。

利用三维激光扫描仪在盾构隧道内部施测往往需面对较为复杂的现场环境，施工器具、金属支架、临时轨道均会对扫描激光准确捕捉隧道内壁点产生不利影响[4，5]。因此，如何通过滤波有效剔除多余的噪声点，同时快速、准确地求解隧道断面参数、提取隧道中心线都是工程中亟须解决的技术难题。虽然国外的一些厂商(如瑞士Amberg)针对此类问题都提供过专业的解决方案，但其硬件设备和软件系统的价格十分昂贵，难以在隧道工程中普及应用。为此本文提出了一套可行的滤波及拟合算法，并已在工程实践中进行验证。

2 主要技术流程

受施工工艺影响，地铁盾构隧道的断面形状一般皆非标准圆，而在数学模型中，正圆可视为椭圆的一个特例，故本文以通用的椭圆柱体作为地铁区间隧道的数学模型[8]。本文的主要技术流程为：

(1)对扫描点云数据进行配准，并手动删除明显、易剔的杂点；

(2)利用天宝公司Realworks软件对隧道点云进行横断面切片处理[6]，该步骤需提供隧道中心线的初始值，一般经两至三次迭代切片便可批量获取高质量的隧道横断面切片；

(3)建立实地测量坐标系与断面坐标系之间的转换关系，利用最小二乘原理拟合模型参数和隧道中心坐标，其计算结果作为初始值进入下一轮迭代，同时每次迭代之前需先完成一次点云滤波(滤波阈值逐次减小)，直至迭代收敛；

(4)记录、保存各断面的拟合参数，包括椭圆模型参数、旋转参数、平移参数及拟合中心点三维坐标，以固定格式输出成果。其流程图如图1所示：

图1 总体技术流程图

3 数学模型及算法关键点

3.1 椭圆柱模型建立

建立隧道模型的关键在于求解模型的几何参数以及确定实地测量坐标系与模型坐标系(即断面坐标系)之间的转换关系。如图2所示，其中O-XYZ为测量坐标系，O′-X′Y′Z′为断面坐标系，O′点为断面(即椭圆)中心，X′轴与椭圆柱中轴线重合，Y′、Z′轴在隧道断面内。

图2 测量坐标系与断面坐标系的示意图

以上两种坐标系均为左手系，二者间的坐标转换关系为：

(1)

式中，△x、△y、△z为坐标系平移参数，R为旋转矩阵，包含9个旋转参数：

隧道断面在O′-X′Y′Z′系下的方程为：

y′2/a2+z′2/b2=1

(2)

写成一般式为：

λy′2+μz′2=1

(3)

将式(1)代入式(3)，按泰勒公式展开后，可构建下列误差方程：

v=a1dλ+a2dμ+a3dr21+a4dr22+a5dr23+a6dr31+a7dr32+a8dr33+a9△y+a10△z+l0

(4)

式中：

l0=λ(△y+r21x+r22y+r23z)2

+μ(△z+r31x+r32y+r33z)2-1

a1=r21x2+r22y2+r23z2+2r21r22xy+2r21r23xz+2r22r23yz

+2r21x△y+2r22y△y+2r23z△y+△y2

a2=r31x2+r32y2+r33z2+2r31r32xy+2r31r32xz+2r32r33yz

+2r31x△z+2r32y△z+2r33z△z+△z2

a3=2λr21x2+2λr22xy+2λr23xz+2λx△y

a4=2λr22y2+2λr21xy+2λr23yz+2λy△y

a5=2λr23z2+2λr21xy+2λr22yz+2λz△y

a6=2μr31x2+2μr32xy+2μr33xz+2μx△z

a7=2μr32y2+2μr31xy+2μr33yz+2μy△z

a8=2μr33z2+2μr31xz+2μr32xz+2μz△z

a9=2λ△y+2λr21x+2λr22y+2λr23z

a10=2μ△z+2μr31x+2μr32y+2μr33z

以上各系数均以未知参数的近似值计算。

由于椭圆方程(2)不含x项，式(4)中不含△x和r11、r12、r13。旋转矩阵为正交矩阵，故式中的6个元素必须满足如下限制条件：

(5)

将式(5)线性化展开后，与式(4)共同构建附有限制条件的间接平差[9]。解算该方程组，应至少选取10个隧道断面点的数据。

3.2 坐标重心化

测量坐标系一般采用国家坐标系或地方坐标系，其平面坐标值远大于高程值，亦远大于断面坐标系下的坐标值。为避免由此造成的数值运算精度损失，同时为增强两种坐标间的可对比性，应对点云的测量坐标实施坐标重心化处理，即用各断面点的三维坐标分别减去同一断面内所有点的对应坐标分量均值，然后以重心化坐标代替测量坐标进行平差计算。

3.3 参数初始化

上述附有限制条件的间接平差方程包含10个未知参数，分别为2个椭圆模型参数λ、μ， 2个平移参数△y、△z， 6个旋转参数r21～r33。另外为求取隧道中心线，X方向上的平移参数△x应一并求解，所以待求的参数初始值共计11个。

(1)椭圆模型参数

椭圆模型参数的初始化一般采用隧道的设计参数，即令椭圆的长轴a与短轴b皆等于设计半径，在初次平差时视其为标准圆进行处理。

(2)旋转参数

旋转矩阵的9个参数实际是由3个独立的旋转角所确定，但因缺少公共点，无法利用式(1)反算旋转参数(或旋转角)的初始值。若以单位矩阵作为其初始值，即假定旋转角的初始值皆为零，则试验证明：当测量坐标系与断面坐标系的x轴(或y轴)夹角小于30°时，迭代平差最终趋于收敛；否则，结果将无法收敛[8，10]。因此该方法不具备普适性。

为实现旋转矩阵的初始化，本文需基于下列两点假设：

①假设两坐标系的Z轴平行。考虑在竖直方向上，Z轴与Z′轴的夹角往往较小，故暂且视其为平行轴，如此可将三轴旋转问题简化为单轴旋转问题(绕Z轴旋转)，仅剩求解1个旋转角θ。

②假设断面切片上所有点的x′坐标值均为0。截取隧道断面时严格控制切片厚度(约 2 cm)，计算时仅考虑y′、z′分量变化，而令x′分量始终为0。

如此可利用式(1)中的第一个方程建立唯一的函数关系式：

x′=△x+cosθ·x+sinθ·y=0

(6)

值得注意的是，式中x、y均为重心化测量坐标。因重心化测量坐标系与断面坐标系的原点在同一平面内，故△x=0，式(6)可变换为：

sinθ/cosθ=-x/y

(7)

正弦与余弦的比值易受误差影响而产生较大偏差，对此应通过多点求解再取均值，然后代入下列限制条件方程求解：

sin2θ+cos2θ=1

(8)

该二次方程具有正、负两根，取正或取负的差别仅取决于断面坐标系X轴的不同指向(大里程或小里程)，不会影响拟合结果，因此计算时可任选其一，进而求出全部旋转参数的初始值。

(3)平移参数

本文采用的滤波方法是基于断面几何中心的迭代筛选，平移参数初始值的偏差过大很可能导致大量有效点云被自动剔除，选择合适的初始值对于后续的点云滤波具有十分重要的意义。

由上文可知平移参数△x始终为零，其余的△y、△z可通过式(4)和式(5)平差求解。考虑到首次平差时点云数据未经滤波，若噪声点在局部较为集中(如铺轨后的轨道及道床扫描点均视为噪声点)，易引起隧道模型的长轴、短轴发生较大的形变，极大地影响拟合精度。因此，本文在初始化平移参数时，省略式(4)中dλ、dμ两项改正，即使用固定的椭圆模型参数计算平移参数△y、△z，得到更准确的初始值。

3.4 迭代滤波

扫描点云中存在大量的噪声点必然会严重影响模型参数和中心坐标的拟合精度，因此本文针对隧道噪声点的分布特点，并结合隧道设计参数(半径)，代入逐级缩小的阈值进行迭代滤波，直至收敛稳定为止。

根据平面椭圆的数学特性，可定义如下滤波函数：f(y′，z′)=(λy′2+μz′2)1/2。对于椭圆平面内的任意一点，若该点恰恰处在椭圆曲线上，则f(y′，z′)=1；当该点落在椭圆内部时，有0≤f(y′，z′)<1；而当该点位于椭圆外部时，有f(y′，z′)≥1。滤波阈值的设置应按从大到小的次序逐级递减(如 30 cm、 20 cm、10 cm、5 cm)，具体数值视实际需求而定。由于滤波函数为无量纲表达式，其阈值也应改用无量纲的数值，即：(1±p/r)，其中p代表某一阈值，r为设计半径值。

点云滤波的基本方法步骤为：

(1)利用既有的平移参数和旋转参数将点云的重心化测量坐标转换为断面坐标；

(2)将各点一一代入滤波函数进行检验，剔除超限点，再以过滤后的点云数据重新平差，计算全部未知参数；

(3)重复上述步骤(1)和(2)，直到相邻两次迭代的改正数小于0.001，停止迭代；

(4)采用更小的阈值，重复上述步骤(1)～(3)，直到最小阈值条件下的收敛完毕。

4 工程试验

本文选择天津地铁3号线某段隧道作为研究对象，该区间位于市内繁华地带，采用盾构法施工，隧道断面的设计内径为 5.5 m。本次试验使用瑞士Amberg的隧道状态激光全息成像扫描系统进行点云数据采集，该系统集成最新的Profiler 6012型360°激光扫描仪与GRP 5000隧道检测小车，其采样频率可达100万点/秒，测距精度优于 ±1.0 mm(相对仪器中心 5 m半径)。为验证本文算法的拟合效果，扫描得到的原始点云采用两种方式处理：一是用系统配套的专业限界及结构变形分析软件，二是用本文阐述的方法及编制的C#软件。

4.1 滤波分析

本次滤波采用的阈值依次为R±20 cm、R±10 cm、R±5 cm，其中R为设计半径 2.75 m。以上阈值分别除以设计半径可得到对应的比例阈值，即：(1±0.07)、(1±0.04)、(1±0.02)。图3所示为某断面作逐次滤波时先后得到效果示意图：

图3逐次滤波效果图

从图3中可以看出，隧道的原始点云中混杂着钢支架、电缆、道床等大量噪声点，利用本文的方法经三次滤波，已将主要的噪声点滤除，剩余的点云数据(除道床区域外)皆均匀分布于隧道内壁，这说明本方法具有较好的去噪效果。

4.2 拟合分析

本次试验扫描的隧道长度约400 m，分别采用TMS软件与Realworks软件基于相同的设计中心线、在相同里程处进行断面切片，相邻切片的间距约为 8 m(弯道及特征点处需要加密)。限于篇幅，仅选取部分数据作对比分析。

图4 拟合中心对比

如图4所示，图中横轴表示隧道断面中心点号C1～C15，竖轴表示由TMS软件或自编程序拟合的实际中心点与设计中心点的点位偏差△P。不难发现二者的拟合结果存在一定的差别：TMS软件解算结果在数值上与采用的设计中心线更为吻合，其△P均值为 1.9 mm，最大值也仅为 5.8 mm；而自编软件的对比结果多介于4 mm～8 mm，仅一点的偏差超过 8 mm(10.8 mm)。该差异产生的原因主要在于：两种软件在指定位置截取的切片会存在细小的差别，进而影响其中心点的拟合坐标。对比而言，作为专业的隧道点云处理软件，TMS明显更具有精度优势及自动化优势，但以本文的方法作为一种更经济的替代方案，也基本满足工程需求。

各断面的椭圆度 表1

表1为通过本方法平差计算出的上述15个断面的长、短半轴及椭圆度。各断面的长半轴与短半轴的较差约 1 cm，其椭圆度值基本处于3.0～5.0的区间浮动，这充分说明较之标准圆柱体，本文选用椭圆柱体作为盾构隧道的模型是更符合实际的。另外由于长、短轴的差别较小，点云滤波环节若作为正圆模型进行处理，会令运算更为简洁高效。

5 结 语

针对地铁隧道的形态及其扫描点云的分布特点，本文提出了采用椭圆柱体进行建模的方法，并推导参数拟合及点云滤波的计算公式，编制相关的计算程序。经对比试验证明：本方法的拟合结果与预期效果相符。在顾及成本的前提下，以本方法替代专业厂商系统软件的部分关键功能具有可行性，可满足工程需求。此外，本方法中的隧道断面切片仍然依靠Realworks等处理软件以半自动化的方式提取，因此探索基于点云的隧道断面自动化截取将是下一步研究工作的重点。