基于点云数据的隧道变形分析系统研究

冯存定，谭志杰，姚竟成，郭 飞，毛朝阳，张瑞鹏

(1.中铁七局集团路桥工程有限公司，陕西 宝鸡 721000；2.西安铁路职业技术学院，西安 710000)

随着西部交通事业的迅速发展，隧道工程的数量呈持续上升状态。在隧道工程施工过程中，隧道变形监测是确保隧道施工质量的重要手段。三维激光扫描技术近几年发展迅速，其所具有的测量速度快、精度高、可三维表达被测物体、抗干扰强等特点，被广泛应用于实际工程中，包括建筑物变形[1]、桥梁工程、隧道工程[2-3]、边坡工程等。托雷[4]运用RANSAC方法对提取的隧道断面进行曲面拟合；程效军等[5]运用三维激光扫描仪进行隧道断面曲线拟合，并分析断面形变、隧道开裂、隧道渗水等；李祖来[6]运用三维激光扫描仪获取隧道点云数据，通过三次样条曲线拟合隧道断面轮廓，用最小距离投影(MDP)分析隧道断面形变量；于宝兴[7]通过三维激光扫描仪分别对地铁、公路隧道进行监测，提取隧道连续断面进行隧道断面变形分析；陈国栋等[8-9]将三维激光扫描技术应用于既有铁路限界检测中。

文中利用三维激光扫描仪对中兰客专后沟梁隧道进行断面扫描获得点云数据，根据线路设计里程及变形监测相关技术规范要求绘制隧道断面图、拟合隧道变形分析点位坐标，基于既有平台实现坐标转换、断面分析、监测报表生成等模块设计，开发基于点云数据的隧道净空三维量测及变形分析系统，从而提高变形监测效率和施工质量，并节约工程投资。

1 系统设计

1.1 系统框架

根据隧道施工监测的内容，在设计系统功能时从隧道点云数据转换和隧道断面分析两个方面入手，对隧道点云数据进行整理和分析。通过VBA语言进行编译，在CAD菜单中形成如图1所示的菜单选项。

图1 CAD菜单功能

根据监测报表生成的需要，在Excel中通过VBA编译，形成如图2所示的菜单选项。

图2 Excel菜单功能

1.2 数据转换模块

1.2.1 隧道中轴线拟合

隧道中轴线是隧道施工过程中保障隧道按设计方向掘进与隧道横断面垂直的虚拟轴线，是隧道施工、监测的基准线，其空间信息是表达隧道走向和空间姿态的重要参考[10]。因此，利用隧道点云数据拟合隧道中轴线是提取隧道断面的重要一步，进而才可以进行断面拟合，分析断面变形信息，完成隧道监测。目前，并没有专业的软件可以直接提取隧道中轴线，通过点云数据快速准确地拟合中轴线成为众多学者研究的重点，也取得了一定的成果，并提出很多中轴线拟合的算法。

虽然从点云数据中提取隧道中轴线的方法有很多，但基本思路都一致，首先提取隧道边界点拟合隧道轮廓，然后根据边界拟合中轴线[6]。虽然获取隧道边界点的方法有很多，为了提高中轴线拟合的精度，文中采用算法编程进行隧道边界点提取。点云边界点提取的算法有极坐标法拟合边界线、方位角法拟合边界线、微分法拟合边界线等。文中采用改进微分法进行拟合。具体步骤如下:

1)将获取的隧道点云数据投影至扫描坐标系的XOY面内，然后指定Y坐标值，寻找该坐标值下X的最大最小值。依次等间隔循环找遍所有Y值，即可得到点云数据投影与XOY平面的边界。

2)指定Y值，并给定微量t值(t=0.25)，通过冒泡法获得区间内X最大最小值，计算X的平均值，即可将此次坐标值作为中轴线上点的坐标，依次寻找到每个Y值对应的中轴线上的点。

3)根据拟合的中轴线点，便可拟合得到中轴线。

CAD中实现隧道点云边界点拟合的主要代码如下：

ReDim maxs(n)

ReDim mins(n)

For j = 1 To m

mx = 0

mi = 90000000

For i = 1 To n

If ZB(i, 1)> min + j * t And ZB(i,1)<=min+(j + 1)* t Then

If ZB(i, 2)>= mx Then

mx = ZB(i, 2)

End If

If ZB(i, 2)<= mi Then

mi = ZB(i, 2)

End If

End If

Next i

maxs(j)= mx

mins(j)= mi

Next j

1.2.2 隧道坐标旋转

点云数据是分布在三维空间的整体数据模型，需要将点云数据进行细化，提取有用的点云数据。因此，需要将点云数据进行坐标旋转，使隧道中轴线与Y轴平行，即可根据Y轴坐标的增加提取与Y轴正交的断面，并且将有用的点云数据保存至Excel数据库中，便于后期断面数据处理时调用。

根据布尔沙模型和隧道坐标转换关系，文中目的是将得到的点云测量坐标转化为点云数据垂直于空间坐标系XOZ面，使点云中轴线平行于Y轴，所以可以简化布尔沙模型[11]，公式中的缩放系数、平移参数、绕Y轴旋转参数都可以不用考虑。简化后的模型公式为

(1)

(2)

(3)

X′=XcosϖZ+YsinϖZcosϖx+ZsinϖZsinϖX,

(4)

Y′=-XsinϖZ+YcosϖZcosϖX+ZsinϖXcosϖZ,

(5)

Z′=-YsinϖX+ZcosϖX.

(6)

CAD中实现坐标转换的主要程序：

For i = 1 To n

Input JHJ3, x, y, z

ZB(i, 1)= y * Cos(angx)+ x * Sin(angx)* Cos(angy)+ z * Sin(angx)* Sin(angy)

ZB(i, 2)=-y * Sin(angx)+ x * Cos(angx)* Cos(angy)+ z * Sin(angy)* Cos(angx)

ZB(i, 3)=-x * Sin(angy)+ z * Cos(angy)

对于任意里程的断面，点云是由离散的一系列点组成的，这时候需要确定阈值d，即断面切片厚度。切片厚度d的断面包含Ya±(d/2)范围内，将这些点投影到Y=Ya的XOZ面上即可对断面进行拟合分析。值得注意的是断面截取厚度d的选取，如果取值比较小，则会造成断面包含的数据量较小，对断面信息提取不全；如果取值过大，运算量会变大，噪点增多，去噪效果不理想，影响断面信息的提取，所以断面厚度d的确定非常关键。文中选择d为5 cm。

1.2.3 方位角排序

隧道断面一般呈现椭圆形或者圆形，三维激光扫描的点云数据在进行隧道点位排序时操作难度大，因此，利用断面点和断面中心点进行坐标反算，计算出方位角，根据方位角大小对断面点云数据进行排序。

已知两个点的坐标，进而计算两点之间的距离和坐标方位角的方法称为坐标反算[12]，计算简单，规律性强，在VBA中实现简便。

CAD中实现方位角计算的主要程序如下：

DltX = X1-X0

DltZ = Z1-Z0 + 1E-20

方法具有较强通用性，对不同型号的零损耗电抗器，只需给定电抗器的特性参数，就可以利用该方法确定最优方案。

A_tmp = Pi *(1-Sgn(DltZ)/2)-Atn(DltX/DltZ)

1.3 断面变形分析模块

1.3.1 断面曲线拟合

对于隧道的变形情况，隧道的曲线部分一般都设计成圆形。采用圆曲线和椭圆曲线拟合隧道的每个断面点云数据相对比较合理，模型简单且计算量小，能很好地反应隧道的整体变形信息。但是从具体的断面来看，隧道在各方向上是不均匀变形的，这时如果用圆曲线和椭圆曲线拟合不能分析出断面上具体的变形位置，也不能体现出局部的变形差异，因此，采用其他曲线模型来拟合断面数据是有必要的。

CAD中有三次样条曲线拟合命令[13]，可以直接调用CAD三次样条曲线命令进行曲线拟合。在调用命令之前需要提取断面数据，选择点作为曲线拟合基础数据。基础点的多少决定了曲线拟合的精度，如果基础点少，曲线拟合效果不佳。如果基础点多，拟合曲线会出现凸点，影响曲线光滑性[14]。断面图为变形监测提供了依据，观察不同时期的隧道断面变形比单一点的位移特征更加直观。

CAD中实现断面曲线拟合主要程序如下：

Dim splineObj As AcadSpline

Dim startTan(0 To 2)As Double

Dim endTan(0 To 2)As Double

startTan(0)= 0: startTan(1)= 0: startTan(2)= 0

endTan(0)= 0: endTan(1)= 0: endTan(2)= 0

Set splineObj = ThisDrawing.ModelSpace.AddSpline(Fitpoints, startTan, endTan)

splineObj.Closed2 = True

splineObj.color = acRed

mmm = splineObj.Area

1.3.2 隧道超欠挖计算

隧道超欠挖计算是隧道开挖掘进施工过程中十分重要的监测工序。超欠挖控制需要实测隧道实际断面，并将其与设计断面进行比较，即准确地测出隧道实际断面轮廓线，并将其与设计断面归入同一坐标系内，可清晰直观地观测到超欠挖情况。文中在拟合隧道断面曲线基础上，将设计断面加载于图形之中，编写程序计算隧道超欠挖情况，用样条曲线进行断面拟合，其结果更加接近于隧道开挖真实情况，计算结果更准确、更精细。

1.3.3 特征点信息提取

在数据转换完成后，数据库中不仅存放着隧道断面的点云坐标数据，还有计算得到的各点与拟合断面中心点的方位角数据。指定方位角可获得相应位置的点坐标，以此获得隧道断面的各特征点坐标数据，为监测报表提供数据支持。每期监测点都以固定方位角进行选择，保证每期数据提取的监测点位置相同。

1.4 监测报表生成模块

在Excel中通过VBA进行编辑，使Excel拥有自动提取特征点数据信息、汇总每期数据、生成监测报表等功能。为了便于数据的管理和调取，每期有效点云数据存储于一个Excel工作簿中，分为两个工作表：一个工作表存放原始点云数据，另一个工作表单独存放提取的各断面特征点数据。同时另建一个汇总表存放每期特征点数据，特征点数据的汇总经过VBA编写，由Excel自动进行数据汇总和整理。新建“监测报表”Excel工作簿，并在第一个工作表中设计监测报表模板，编写VBA程序，自动调取相关数据，计算对应的变形量，对变形量设置阈值，判断变形量大小，进而对超过阈值的变形值进行特殊标记，以起到警示作用。

2 工程应用

2.1 工程概况

研究区域为中兰客运专线后沟梁隧道(见图3)，该隧道位于甘肃省兰州市兰州新区中川镇境内，是中兰客专兰州新区联络线右线关键工程。隧道起讫里程为LDyK241+056～LDyK241+350，全长294 m，为单线隧道，设计时速160 km·h-1，为满足过梁需要，采用时速250 km·h-1标准设计。全隧道位于R=1 600 m的右转曲线上，洞内纵坡为-25‰/294 m的单面下坡。

图3 中兰客专后沟梁隧道

2.2 点云数据采集及预处理

2.2.1 数据采集

采用Z+F 5010C型三维激光扫描仪作为隧道点云数据获取的仪器。测站间隔为25 m，设站位置位于隧道中线，保证成“一”字型排列，且随着隧道施工的进度合理地安排测站数。三维激光扫描仪采集的点云数据是在以测站为原点建立的坐标系下，为了使点云数据规划到统一的隧道施工坐标系下，需要在隧道外或隧道内固定位置布置标靶，并采集基准标靶在隧道施工坐标系下的坐标。隧道内采用全站仪进行基准标靶坐标采集，隧道外可用全站仪或RTK进行采集。在三维激光扫描仪进行隧道数据采集时，可以自动识别标准靶纸或靶球，为后期点云数据坐标转换提供高精度的基准数据。基准标靶的布设应该遵循分布均匀、避免共线、基础稳定可长期保存的原则。隧道洞内靶球是点云数据拼接的重要依据，因此，在布设时应固定在相对稳定的地方，防止靶球移动影响数据采集。隧道洞内进行靶球布设时，在测站两侧大概12.5 m处布设3个靶球，要求靶球不共线、分布均匀、与测站通视没有障碍物遮挡。

2.2.2 隧道点云预处理

三维激光扫描仪的扫描距离可以达到几百米，但是在隧道这种特殊的线状建筑物中，扫描仪采集的点云密度和点位精度会随着距离的增加大幅度下降。为了获得符合精度要求的高密度隧道点云数据，需要将扫描仪的扫描范围控制在几十米以内，因此，全隧道点云数据的获取必须经过多次扫描才能获得。

在点云数据采集过程中，每设站1次都会以本站为坐标原点建立测站坐标系获得点云数据，而且存在大量的冗余点云和噪声点云，所以采集的点云数据需要经过预处理才可以使用。隧道点云数据需要进行点云拼接、点云去躁、点云压缩等工作。将经过预处理的点云数据导出为XYZ格式的点云文本数据便可在系统中使用。

2.3 隧道断面分析

2.3.1 数据转换

选择功能菜单中数据转换模块的中轴线拟合命令，导入经过预处理的隧道点云数据，按照系统提示输入对应指令，系统将筛选隧道点云投影在坐标面的边界点，并拟合隧道在该投影面的中轴线。通过中轴线计算隧道与坐标轴的夹角。图4为隧道投影在XOY面的图形及计算得到的隧道与中轴线夹角。

图4 中轴线及夹角

已知隧道与坐标轴的夹角，可选择数据转换菜单中的坐标转换命令，根据系统提示输入相应的夹角开始坐标转换，坐标转换原理在上节中已经进行详细说明。完成坐标转换后的点云数据存放在新建的Excel数据库中[14-17]。数据转换完成后可根据系统提示选择是否将转换完成的隧道点云数据展绘在CAD中。图5为存有点云数据的Excel数据库和展绘的隧道点云数据。

图5 数据库及隧道点云数据

得到坐标转换后的隧道点云数据可选择方位角排序选项，系统按照上节方位角排序的原理进行后台计算，并将方位角值写入Excel数据库中。图6所示为方位角排序后的结果，图中第4列为方位角值。在方位角排序时，系统提供操作者两种方式：一种是对所有断面进行排序，缺点是排序过程耗时较长，优点是1次排序即可；另一种是指定断面排序，优点是仅排序需要断面，耗时短，缺点是调用后续断面时还需进行排序。两种模式各有利弊，根据实际需要进行选择。

图6 排序后的点云数据

2.3.2 断面分析

对点云数据进行数据转换后开始进行断面分析。选择断面分析菜单中的断面提取命令，输入数据库文件，系统将链接到相应数据库文件，进行断面提取，系统设置了全隧道断面提取和指定隧道断面提取，方便操作者根据不同需要进行提取、展示和分析。图7所示为提取的连续隧道断面。

图7 连续隧道断面

连续断面的提取更多的是为了展示，单一断面的提取是为了进行分析。提取单一断面后进行断面曲线的拟合，选择菜单中的断面曲线拟合命令，系统可按方位角筛选自动拟合隧道断面曲线，但其拟合效果取决于断面点筛选的精细程度；也可通过手动选点进行断面曲线拟合。图8所示为手动选择20个断面点进行断面拟合的断面曲线。

图8 拟合隧道断面曲线

完成隧道断面曲线拟合工作，加载设计断面曲线。调用隧道超欠挖计算程序，输入超欠挖计算指令“CC”，根据系统提示，先后选择拟合的隧道断面曲线和设计断面曲线，系统进行超欠挖计算，并将超欠挖结果显示在CAD图形中。图9所示为超欠挖计算示例。

图9 隧道超欠挖计算

超欠挖计算结束，可选择菜单中的监测点提取命令，进行断面监测点的提取。系统根据设置好的阈值以方位角大小进行特征点提取，将提取的数据保存在Excel数据库中，为监测报表的生成提供数据支撑。

2.3.3 监测报表

系统设计的变形监测报表自动生成部分在Excel数据库中实现。首先在按监测日期保存的数据库文件中，选择监测报表下的“提取特征点”选项，系统将会自动提取特征点数据存放在特征点数据模板中，图10所示为提取的特征点数据。

图10 特征点数据提取

特征点数据的汇总：打开存放在监测报表文件夹中的特征点数据汇总文件，点击 “监测报表”中的“特征点数据汇总”选项。系统将调取存放在“监测报表原始文件”文件夹中所有期数据中的“特征点数据”工作表，将其复制于特征点数据汇总工作簿中。图11所示为汇总的3期特征点数据，工作表以监测日期命名。

图11 特征点数据汇总

监测报表自动生成：打开存放在监测报表文件夹中的监测报表文件，点击“监测报表”中的“监测报表生成”选项。系统会根据输入的日期调取特征点数据汇总文件中相应数据，并计算与前一期数据的变形值与初期数据的变形值，如图12所示。

图12 变形监测报表

3 结 论

1)研究隧道点云中轴线提取方法，并对此进行编译，形成子菜单选项，可在软件中直接操作，提高隧道点云数据处理的效率。根据得到的隧道点云中轴线，对点云数据进行坐标转换，使点云数据Y坐标与隧道中轴线大致平行，根据隧道里程进行连续断面提取，并将断面点云数据保存在Excel数据库中，进而通过提取断面点云进行变形分析。

2)调用CAD三次样条曲线命令，根据断面点拟合隧道断面曲线，并计算断面曲线所围闭合图形的面积；加载设计断面，平移后计算隧道超欠挖情况，将超欠挖信息标注在CAD图形中。

3)以AutoCAD软件为平台、VBA为二次开发工具、Excel作为数据容器，结合ActiveX技术，建立AutoCAD环境中的对象和VBA开发工具之间的联系，使它们之间能够快速地进行数据交换和绘图，形成以处理隧道断面数据为主要功能的变形分析系统。

文中数据分析系统在实际工程中进行了检验，能够快速提取断面进行分析，满足施工精度要求，但是系统功能还不全面，可以在此基础上继续进行二次开发，使系统的应用更为广泛。