利用隧道激光点云提取中轴线及进行整体变形分析

林景峰，俞家勇，田茂义，徐飞，周茂伦

(1.山东科技大学 测绘科学与工程学院，山东 青岛 266590；2.安徽建筑大学 土木工程学院，合肥 230601；3.青岛秀山移动测量有限公司，山东 青岛 266590)

0 引言

对隧道工程进行变形监测是保障施工及运营安全必不可少的工作之一。目前隧道监测方法使用收敛仪、断面仪、全站仪等进行变形分析，这些传统测量方法可以采集到高精度的隧道断面数据，但采集监控点数据的工作效率低、控制点易被破坏等突出问题对隧道在施工与运营过程中造成了大量的经济损失。三维激光扫描技术具有不接触、高密度、高精度、快速性、自动化等优点，随着三维激光扫描技术在国内外的高速发展，在隧道工程领域三维激光扫描技术逐渐替代了传统测量方法[1-2]。

Han等[3]利用最小距离投影法直接估计三维离散点云，识别出整个隧道表面任何点位的变形位移量，提出了一种严格的协方差传播方法，对隧道数据进行降噪，并通过对公路隧道的仿真实验和实例分析，验证了该方法在变形监测中的实用性。王令文等[4]以隧道设计中轴线为基准提取隧道横断面，通过建立横断面上的局部坐标系，使坐标系原点与横断面所在的设计中轴线点重合，分析理论断面和实际断面的径向差值，得到隧道收敛变形情况。徐光华[5]根据中轴线提取横断面并结合粗滤噪和精滤噪对点云数据进行去噪，通过叠加2期点云数据建立的数字表面模型分析隧道变形。李珵等[6]基于双向投影法获取隧道中轴线并提取隧道横断面，采用极坐标法进行隧道局部变形分析。卢小平等[7]结合双向投影法和二次曲线拟合隧道中轴线，通过构建极坐标叠加分析多期数据，以获取隧道在各个方向上的位移变化量。许磊等[8]提出一种基于射线相交的隧道断面自动测量方法，根据设计线位提取断面线，并采用RANSAC算法进行断面滤波，滤波后的断面数据与全站仪数据对比误差可用于隧道变形分析。纪思源等[9]基于隧道点云双向投影提取中轴线，根据中轴线进行点云分割，并对截取的横断面进行椭圆拟合，精度满足隧道变形分析。Chen等[10]基于控制网络和最小二乘优化方法使点云配准误差最小化，通过将隧道三维点云投影到水平面上，转换成二维平面图像，然后采用骨架化方法估计隧道的三维中轴线，并以此为基准截取隧道横断面。通过实例验证，该方法能够连续提取精度较高的正交横断面进行变形分析。杜黎明等[11]将隧道三维数据投影至二维平面，生成二维图像提取中轴线，之后以k近邻法计算各断面的缓冲区，并采用投影法构建断面点集，然后使用迭代椭圆拟合方法拟合出轮廓线，通过与全站仪数据比较，误差在隧道变形范围内。杜盼晓等[12]证实了竣工后隧道并不是严格的椭圆形，对隧道竣工后与设计模型之间的变形进行分析，结果表明径向方向上可以监测大于15 mm的变形，验证了三维激光扫描技术在隧道变形分析中的可行性。

以上学者证实了三维激光扫描技术在隧道变形监测中的实用性与可靠性。本文针对双向投影方法难以提取全局中轴线的问题，基于单向投影提取水平中轴线并以此截取横断面，采用迭代拟合空间圆获取隧道全局三维中轴线。为了综合分析隧道整体与局部变形的情况，本文根据隧道横断面点到三维中轴线与设计半径的径向变形量渲染点云数据，使隧道整体变形可视化，依据可视化效果分析隧道整体以及局部变形情况。

1 隧道整体变形分析方法

首先，将隧道三维点云数据投影于XOY水平面，通过构建Delaunay三角网提取出边界点集；其次，基于水平边界点与隧道中心的几何关系提取水平中轴线；然后，根据水平中轴线将整条隧道切割成若干段点云数据，对每段点云使用投影法得到横断面点云，并迭代拟合空间圆进行降噪获取隧道中轴线；最后，通过获取横断面点到拟合圆心的距离与隧道设计半径的相对变形量，绘制隧道整体变形图。算法流程图如图1所示。

图1 隧道整体变形分析算法流程图

1.1 水平中轴线提取

隧道中轴线能够准确地提取横断面及表达隧道整体走势，目前中轴线提取有基于双向投影与单向投影的方法。双向投影方法在直线型隧道能够提取精度较高的中轴线，但是在弯曲型隧道中难以提取完整的中轴线。相比于双向投影方法，单向投影方法在弯曲型隧道中能够提取整条隧道的中轴线，且精度满足隧道工程误差容忍度[13-14]。因此，本文首先将隧道三维激光点云投影至XOY水平面以便提取点云边界；其次基于下采样边界特征提取算法获取边界点，利用垂线公理及中点公式提取初始中轴线；最后使用距离阈值算法剔除粗差点，并根据工程需求提取一定间距的水平中轴线用于截取横断面。

1) 基于平面几何关系提取初始中轴线。目前提取平面边界算法较多，如经纬线扫描法、网格划分法以及三角网法等，其中Delaunay三角网算法适用于带状数据的边界提取。Delaunay三角网算法的边界提取原则：Delaunay三角网在生成三角网过程中，以最近的三点形成三角形，且各三角形的边皆不相交，每条边都使用所在的三角形的顶点验算并记录相应的2个顶点坐标。整个验算过程中，若一条边只被一个三角形使用过，表明这2个顶点所确定的边是边界边，存储这2个顶点坐标；若一条边被2个三角形使用过，表明这2个顶点确定的边在整个三角网内部，将其剔除[15]。因此，一般采用Delaunay三角网算法提取边界点，是对全部平面点云进行构网，计算时间较长。

为了提高边界提取的效率，本文提出一种基于下采样的边界特征提取算法。该方法通过点云分块和三角化进而提取出精确的边界特征，其具体步骤如下。

(1)如图2(a)所示，在XOY平面内，沿点云的X轴方向设置间距为ΔX将平面点云分成若干块点云数据，并提取每一块点云中xmin、xmax、ymin、ymax对应的点组成边界点集 Ⅰ。为了避免当隧道垂直于X轴时提取边界点集较少情况，沿Y轴方向以同样的方法提取边界点集 Ⅱ。其中，ΔX取值根据工程实际需要进行选取。

(2)利用这2个点集组成隧道初始边界点集。

(3)根据初始边界点云的密度自适应设置三角形边长，以e≤2ΔX为限制条件构建Delaunay三角网，最后提取三角网中只被一个三角形使用过的边界边对应的2个顶点构成隧道边界点集，如图2(b)所示。

图2 水平边界点集提取

图3 左、右边界线提取

如图3所示，基于下采样边界特征提取算法获取的边界点集，根据选取隧道点云数据的长度s，区分隧道边界与隧道断面分割线a与b的几何关系，并依据隧道前进方向提取左、右边界线。

图4 水平中轴线提取示意图

对于平行曲线型的2组数据，使用KD树算法能够快速准确地在另一组数据中找到与之对应的最近k个点。如图4所示，以左边界点集Pi(xp i，yp i，zp i)构建KD树，沿前进方向搜索右边界最近的2个点Qi与Qj，利用垂线公理获取在右边界垂足点Hi。根据每个左边界点Pi与其对应的垂足点Hi,求得中轴线点集Gi(i=1，2，…，n)(式(1))。

(1)

为了获取更加准确的中轴线，对右边界以同样的方法获取一列中点Mi，构成初始水平中轴线点集。

2) 基于距离阈值优化水平中轴线。若水平中轴线存在粗差点，会导致提取的横断面与隧道不垂直。本文根据曲线长度和控制点密度，对初始水平中轴线进行初步插值，并拟合中轴线剔除粗差点，但是在狭长的空间内，若使用多项式整体拟合中轴线会存在偏差。因此,本文提出利用距离阈值的方法，判断初始水平中轴线点与加密点是否为控制点，分段拟合水平中轴线。如图5所示，根据点到直线距离d，剔除超出阈值εd的中轴线上的粗差点。该方法的步骤如下。

(1)根据设计的距离阈值εd，在图5中以第一个控制点a1为起点，计算控制点a2到直线a1a3的垂直距离d1(式(2))。

(2)

(2)若a2是控制点，则以其为起点计算点a3到直线a2a4的距离d2；否则以下一个点a3作为循环的起始点，计算点a4到直线a3a5的距离。

(3)按此步骤，重复计算至中轴线的末端，将大于距离阈值的粗差点剔除。

图5 距离阈值方法原理

其中，ai(i=1，2，…，n)是初始水平中轴线经过加密后的控制点，按照隧道前进方向排序组成的点集。距离阈值εd大小决定剔除控制点数量多少，太大会导致部分粗差点未被剔除，反之保留控制点数量太少，因此，本文以每个点ai到其相应直线的垂直距离di(i=1，2，…，n)的平均值作为设计的距离阈值εd(式(3))。

(3)

因为通过均匀分布的横断面与设计模型的相对变形量能够较好地分析隧道整体变形情况，所以对剔除粗差点后的初始水平中轴线，根据工程需求提取间距为g的水平中轴线。

1.2 横断面提取

利用横断面相对变形量能够获取隧道变形情况。对于直线度较好的直线型隧道，可以直接基于中轴线提取与隧道正交的横断面，但是对于复杂的弯曲型隧道，需要将其分割成若干直线段再进行数据处理[16]。因此，本文首先基于水平中轴线切割隧道点云，其次采用投影法获取横断面点云，最后通过迭代拟合空间圆方法进行降噪，以提取精确的隧道中轴线与横断面点云数据。

根据中点和法向量，则过中点处的法平面方程如式(4)所示。

x-xg0+k(y-yg0)=0

(4)

根据中点所在的法平面方程，提取距离此法平面g/2范围内的点云[17]，将点云使用集合Ψ(xi，yi，zi)表示，如式(5)所示。

(5)

2)基于投影法获取横断面点云。虽然隧道点云密度很大，但是落到法平面上的点较少。如图6所示，设过中点G0的法平面为M，将分割的点云至法平面距离小于g/2的点均视为法平面上点，并将这些点投影到法平面上，投影变换如式(6)所示。

(6)

式中：T′i(x′i，y′i，z′i)为距离法平面小于g/2的点云投影到法平面上对应的坐标，使用Ti(xi，yi，zi)表示。在空间坐标系下，T′i投影前、后在3个坐标轴方向上的变化量为Δxi、Δyi和Δzi，在法平面上相应的坐标Ti的表达如式(7)所示。

(7)

图6 横断面点投影示意图

3)横断面降噪。隧道是一种在封闭空间内狭长的建筑物，其点云内部噪点影响点云精度分析[18]。原始点云经过预处理后，还存在部分噪点。为了精确分析隧道的变形，采用迭代拟合空间圆的方法降噪，并通过半径内符合精度评价降噪精度。

(1)拟合空间圆。如图7所示，空间圆是球体S与平面P的交线。球面方程如式(8)所示。

(x-x0)2+(y-y0)2+(z-z0)2=R2

(8)

式中：(x0，y0，z0)为球心O；R为球体半径。

根据法平面点Ti到球心的距离的关系|TiO|=R，将法平面点代入式(8)得到式(9)。

(xi-x0)2+(yi-y0)2+(zi-z0)2=R2

(9)

根据间接平差原理建立误差方程(式(10))。

(10)

式(10)是非线性的误差方程，将误差方程线性化得到式(11)[19]。

(11)

根据式(4)所求法平面如式(12)所示。

Ax+By+Cz+D=0

(12)

式中：A=1；B=ki；C=0；D=-xi-kiyi。

联合式(6)、式(7)及式(12)，将式(11)所求的球心坐标O(x0，y0，z0)投影至法平面上，得到空间圆的圆心O1(x1，y1，z1)，并求得球心至法平面的垂直距离D1与空间圆半径R1，如式(13)所示。

(13)

图7 拟合空间圆示意图

(2)迭代拟合空间圆降噪。隧道内表面粗糙度、连接螺栓和螺帽等因素会影响隧道变形分析的精度。本文结合隧道设计半径与隧道变形以及噪点，通过设定径向阈值δ，采用迭代拟合空间圆的方法逐步剔除噪点，以获取精确的轮廓线，并提取降噪后的圆心与半径。圆心作为隧道中轴线绘制隧道纵断面图，半径作为点云内符合精度评定的基准。

首先，获取横断面上任意一点Ti到拟合圆心O1的距离λi与拟合圆半径R1的径向差值Δλi，通过比较径向差值Δλi与径向阈值δ的大小，对大于径向阈值的横断面点予以剔除，保留小于径向阈值的横断面点。

(14)

Δλi=λi-R1

(15)

然后，通过设置递减阈值Δδ和迭代次数k,对保留的横断面点云数据进行迭代空间圆拟合，逐步剔除噪点，提取降噪后横断面数据拟合的圆心作为中轴线点。

最后，按此步骤对每个横断面依次降噪，提取拟合圆心构成隧道中轴线，并采用点云内符合精度的中误差评估降噪精度(式(16))。

(16)

式中：Δ为横断面点到相应拟合圆心的距离与拟合圆半径R1的差值；n为横断面点云数量。

1.3 整体变形可视化

一般将扫描点云经过配准、分割、降噪等过程，提取实测断面数据，通过实测断面与设计断面对比，获取隧道整体变形。但是实测断面以一定间隔提取，断面提取位置存在随机性，导致隧道的局部相对变形被忽略。同时，通过拟合得到的近似断面与实际断面的点云数据存在一定的偏差，导致隧道变形分析不准确。因此，本文根据点云数据全覆盖的特点，通过横断面上的点Ti到拟合圆心O1的距离λi与隧道设计半径R2的径向差值Δd渲染点云，然后将附上颜色属性的点云投影反算至初始位置，使隧道整体变形可视化，同时直观地表达各区域的局部变形情况。隧道的整体变形情况如图8、式(17)所示。

(17)

图8 整体变形分析图

2 实验分析

本文实验数据是利用FARO Focus 3D X130采集的一段隧道数据，该段弯曲型的圆形隧道全长为1 000 m，起点里程为K20+100，数据采集于运营阶段，半径设计值为R，最后通过Faro Scene 7.0进行数据预处理并导出扫描点间隔为2 cm的隧道点云数据文件。通过数据预处理与VC++编程，得到隧道整体变形，并根据隧道整体变形的可视化效果分析隧道局部变形情况。

2.1 基于边界点提取水平中轴线

本文采用下采样边界特征提取算法提取边界点，通过实验分析，ΔX为原始点云密度的2.5倍时能够有效保留边界特征，同时保证边界提取的效率。首先，沿X轴和Y轴方向设置间距Δx=Δy=0.05 m的分割平面点云；然后，提取对应的点组成初始边界点集；最后，限定三角形边长e≤10 cm，构建Delaunay三角网，获取隧道精确边界点集。使用本次实验数据(一共2 278 015个点)与重采样数据(一共1 141 315个点)验证下采样边界特征提取算法的适用性，对这2组数据分别采用传统Delaunay三角网算法与下采样边界特征提取算法获取隧道边界，提取边界点集数量与计算时间见表1。由表1可知，传统Delaunay三角网算法随着点云数量的增加构网时间明显增加。以传统Delaunay三角网算法为基准，利用下采样边界特征提取算法粗提取全部实验数据的边界点减少平面点云数据达到97%。通过下采样边界特征提取算法精确提取2组数据的边界点集，数量均约占传统Delaunay三角网算法精确提取边界点数量的97%，说明使用下采样边界特征提取算法并不会损失隧道边界特征。对于数据量为1 141 315的重采样数据集，下采样边界特征提取算法计算效率相比于传统Delaunay三角网算法提高并不明显，但对于数据量为2 278 015的全部实验数据集，计算时间比传统Delaunay三角网算法缩短了5 146 s。

表1 传统Delaunay三角网算法与下采样边界特征提取算法比较

基于边界点集提取左右边界线，使用式(1)获取隧道初始水平中轴线控制点，采用距离阈值方法优化水平中轴线，并提取间隔为0.2 m的水平中轴线作为提取横断面的基准及表达隧道整体走势(图9)。

图9 隧道整体走势图

2.2 拟合横断面及精度分析

基于提取的水平中轴线控制点将点云分割为4 908个直线段，并计算每相邻2个中轴线控制点在中点处的法平面，实现对点云数据提取。本文采用迭代拟合空间圆的方法在每个法平面上进行降噪，通过实验分析，设置初始径向阈值δ≥±0.1 m，递减间距Δδ=0.005 m，迭代次数k=10时，能有效地剔除隧道壁以外的噪点，同时不破坏隧道管壁点云。为了更好地展示隧道整体横断面提取效果，提取间距为5 m的横断面和中轴线点进行整体显示，如图10(a)与图10(b)所示，并提取了一段间距为0.2 m的5 m局部数据表达本文算法拟合中轴线以及横断面降噪效果，如图10(c)所示。

图10 降噪隧道横断面与中轴线

通过分析横断面拟合的精度，检验点云的内符合精度及横断面降噪情况。统计所有降噪后横断面点到相应拟合圆心的距离与拟合圆半径的径向差值，结果如图11所示。图11中，点云的偏移量主要集中在±2 cm范围内，最大不超过4 cm；经平差计算，降噪后的点云平均偏差为0.5 cm，中误差为1 cm，说明降噪效果较为理想，点云的内符合精度较高。

图11 点云的内符合精度分析

2.3 隧道整体收敛变形分析

通过横断面点云的径向差值渲染，得到隧道整体变形的可视化图形。因整段隧道数据量较大，不便整体显示，选取里程K20+245～K20+520的隧道数据分析隧道整体变形情况，如图12所示。

图12 隧道整体变形分析图

根据图12(a)的可视化效果，分析隧道的整体变形情况。图中隧道顶部为黄色，两侧呈现青色和绿色，符合隧道在运营阶段的变形趋势[20]。根据可视化图形的颜色变化情况，该段隧道可分为3种主要变形，分别提取相应位置的断面进行变形分析。为了更好地表达横断面相对变形，将变形量放大20倍(图12(b))。横断面505邻近的两侧区域呈现青色，顶部区域呈现黄色，表明隧道在该局部区域内变形为扁平的圆形。横断面948邻近的两侧区域呈现绿色，顶部区域呈现黄色，表明隧道该区域内横向均匀收敛和纵向均匀沉降。横断面1114邻近的两侧区域呈现绿色，顶部局部区域呈现橙色，表明隧道在该区域内纵向沉降变形较大。此外从图12中可以观察出，该段隧道局部变形不连续。初步推断该段隧道受外部因素，导致局部变形不连续，但是变形量在隧道变形的可控范围内。

通过本文算法对整段隧道进行整体变形分析，得知该段隧道在安全变形范围内，但受外部及其他因素的影响，导致隧道局部存在微小变形。在后期进行隧道安全检测等工作时，为了防止隧道变形由量变引起质变，应着重关注隧道整体变形，分析图12中的颜色变化区域，具体分析其变化因素，保障隧道安全运营及人员安全。

3 结束语

本文基于三维激光扫描技术，研究了利用隧道三维点云提取隧道中轴线及整体变形分析方法，并在某段中长距离的弯曲型隧道中进行实验，提取了精度较高的水平中轴线用于截取横断面，通过横断面点云的径向差值得到隧道整体变形分析图。在中轴线提取过程中，采用下采样边界特征提取算法提高了边界点云提取的效率，通过距离阈值方法剔除粗差点得到精度较高的水平中轴线。在横断面提取过程中，采用迭代拟合空间圆的方法获取轮廓线，减少了传统轮廓线提取旋转点云的过程。最后，通过点云相对变形量获取隧道整体变形分析图，直观地描述隧道的整体变形趋势。本研究为井筒等圆型建筑的动态变形监测分析提供了数据基础。