基于点云的劲性骨架拱桥自动化线形监测方法

白祖应， 卯申殷， 胡开心， 王义成， 周银， 王玥， 应春莉， 韩达光

(1. 云南交投普澜高速公路有限公司， 普洱 665000； 2. 重庆鲁汶智慧城市与可持续发展研究院， 重庆 401135；3. 重庆交通大学土木工程学院， 重庆 400074； 4. 奥斯陆城市大学技术&艺术与设计学院， 奥斯陆 0130)

近二十年来,劲性骨架拱桥因其受力性能优越、跨越能力强、施工便捷，以及抗震性能好等诸多优点。大跨度劲性骨架拱桥在中国得到了迅猛发展，尤其是在山区桥梁中占有重要位置。根据相关统计表明，随着拱桥跨径的增大，钢管混凝土劲性骨架拱桥在钢筋混凝土拱桥中占比呈现出递增趋势[1]。尤其是400 m以上跨径拱桥中，占有绝对优势[2]。随着云桂铁路南盘江特大桥(主跨416 m)、沪昆高铁北盘江特大桥(主跨445 m)的相继建成，劲性骨架拱桥跨径不断突破。目前在建的天峨龙滩特大桥将达到600 m[3]。类似其他大跨径桥梁，拱桥的拱圈线形测量对分析桥梁的结构受力状况是十分必要的。拱桥的拱圈线形是其在长期的动、静荷载和自然老化作用下，被认为桥梁几何形态变化及发展趋势是桥梁结构安全预警的关键指标与依据[4]。然而，对于这类大跨径桥梁来说，使用接触测定的方法进行测量操作难度大、成本较高，传统的全站仪、水准仪等测定仪，由于人工操作因素较多、定位点有限，导致在桥梁检测时精度难以控制[5]。

根据测量手段，可以将当前测量分为两大类：一种为单点式测量，另一种为全覆盖式测量。单点式测量即传统的全站仪、激光干涉仪、全球定位系统[6](global positioning system,GPS)等技术。该类方法可实现单点的精准测量，操作简易，其中某些全站仪可达毫米级别精度。但由于该类测量的测点有限，人工影响因素大，故难以保证结构的整体高精度线形。全覆盖式测量包括摄影测量、三维激光扫描等方式。由于近景摄影测量综合精度较低，较难达到工程高精度测量要求[7]。三维激光扫描通过获得物体表面的三维坐标和密集信息，被誉为继GPS技术后测绘领域的又一次技术革命。近年来，越来越多的学者将三维激光扫描应用于工程测量。Armesto等[8]在古建筑拱桥的变形监测中，利用三维激光扫描技术获取点云，并以统计非参数的方法处理得到拱的精确几何尺寸，以此来分析拱的变形。Park等[9]对钢梁的变形监测做出了研究，通过提取扫描点云的单点数据变形获取。中国在三维激光扫描技术应用上起步较晚。徐进军等[10]根据斜拉桥桥塔的特点，通过提取桥塔棱线处的激光点云数据，以获取不同荷载情况下整个桥塔的挠度变化。任杰[11]对某高铁桥采用三维激光扫描技术进行拱桥基础支座、桥梁结构拱顶进行连续定点定时扫描，达到灌注混凝土砂浆时线形监测的目的，并与全站仪进行对比验证了准确性。邓晓隆等[12]提出了一种基于标靶点云的桥梁变形监测方法，有效地避免了单点测量难以获取结构全面变形的缺点。

目前，前人研究多是通过扫描拱桥或其他桥梁有限单点的关键特征方式，并配合相关算法完成定位点或线形的提取。然而，对于拱桥拱圈分环分段浇筑及拱上立柱等复杂工序来说，拱圈线形连续、高精度的自动化提取是拱桥施工监测的关键，目前相关的研究尚少，也是目前行业急需解决的难题之一。基于此，采用三维激光扫描仪获取多期该拱桥拱圈底部点云数据，并提出一套拱桥自动化高精度线形提取算法；进一步的，将多期线形做基于曲率的径向滑窗算法叠差变形提取，并与现场全站仪实测数据进行对比验证。该变形将对劲性骨架拱桥拱圈分环分段浇筑及拱上立柱等重要施工工序的施工，以及拱桥的后期运维也具有重要潜在意义。

1 拱圈变形提取

1.1 算法概述

由于施工时间的不同，现场的点云采集工作会分多期进行，导致每一站的数据处于不同坐标系下。故首先将采集的多期点云进行拼站处理。然后对拱桥点云进行降噪，降噪的主要对象为第三类噪点。考虑到拱圈底部点云为圆弧状分布的特点，选用三棱锥法去噪。再定义最佳拟合表面对拱圈底部点云进行分区域拟合，一般来说，考虑到拱桥的空间变形，将其划分为上中下3块等分区域；进一步的，基于拟合多项式曲率的径向滑窗算法叠差变形提取。点云处理流程如图1所示。

n为循环的次数图1 点云处理流程Fig.1 Point cloud processing flow

1.2 点云预处理

点云拼站是不同期点云数据对比的基础，该步骤实质是将位于不同坐标系的点云统一坐标系。通过迭代的方式使得对应点间的平方距离之和最小从而确定两片点云之间的旋转平移矩阵，ICP(iterative closest point)算法是应用最广泛且影响最深的配准手段。

(1)

式(1)中：si为源点云中任意点i；mc(i)为源点云s中任意点i在目标点云M中的对应点；ICP算法采用迭代求解旋转矩阵Rrot；t为平移向量。

点云噪点是影响点云精度的关键因素，故本文对点云去噪进行重点研究。除仪器本身及目标物特性所产生的第一、二类噪点外，第三类噪点是主要需要处理的噪点。第三类噪点包括：漂移点、孤立点、冗余点、以及混杂点。前3种都可以采用人机交互的方式进行可靠的去除，第4种小尺度噪声是算法去噪的重点研究对象。

常见的去噪算法为中值滤波算法、均值滤波算法、高斯滤波算法等，对于曲面点云，采用适应性较强的三棱锥法去噪，如图2所示。假设任一点P作为三棱锥的顶点，另外任意3个连续点P1、P2和P3连接成底面。由于空间坐标点可知，该三棱锥3个侧面的顶部夹角α1、α2和α3。该算法将通过顶点3个夹角的大小，来判断点P是否为噪点。首先排除当P点为噪声点的条件是，其中任意两个夹角均属于[90°，180°) 范围。接着，根据点云噪点程度不同，分为如下两级噪点。当α1、α2和α3同属于(60°，90°]，为一级噪声点;当α1、α2和α3同属于(30°，60°]，为二级噪声点;当α1、α2和α3同属于(0°，30°]，为体外孤点。将评估判断后点移至统计的正确位置或者剔除弥补扫描误差，最终得到去噪后的点云。

P为噪点顶点；P1、P2、P3为与噪点P形成棱锥的底面顶点；C1、C2、C3为下底面三角形边长；L1、L2、L3为棱锥棱长；α1、α2、α3 为三棱锥三个侧面的顶部夹角。图2 三棱锥法去噪原理Fig.2 Denoising principle of triangular pyramid method

1.3 多项式曲面拟合

多项式曲面拟合是求得最佳拟合表面算法的基础，故先对其原始做详细的介绍。设拱圈底部点云实际高成为Zi、曲面拟合高程为zi、拟合差为ξi，则三者的关系可表示为

zi=Zi-ξi

(2)

多项式曲面拟合模型认为ξi与平面坐标(xi，yi)之间存在着相关关系，即

ξi=f(xi,yi)+εi

(3)

式(3)中：εi为随机误差；f(xi,yi)为多项式曲面函数关系式，可表示为

(4)

当有n(n≥m)个已知Zi、zi点和xi、yi点(以下称为已知点)时，由式(3)在∑ε2=min的条件下，其中，ε为随机误差，min为随机误差平方和的最小值，∑为对样本所有点的随机误差平方值求和。运用最小二乘法求得aj=(j=0,1,2,…,m-1)的最小二乘估值，然后将aj，和待定高程的点(以下称为待定点)的坐标同时代入ξ=f(x,y)，进而由式(2)求得该点的zi。

为区别已知点的ξ和待定点的ξ，将已知点中的ξ用Li表示。将式(3)写成矩阵形式为

L=AX+ε

(5)

式(5)中：

其中，n≥m时以ε的数学期望E(ε)=0为条件，求得其最小二乘解，可表示为

X=(AΤPknowA)-1(AΤPknowL)

(6)

式(6)中：Pknow为已知点观测值的权阵。

多项式曲面拟合是一种精度较高，且鲁棒性好的曲面拟合算法。但由于多项式拟合曲面次数未知，不同阶次的曲面拟合结果精度有差异。因此，为了达到最佳拟合效果，现定义最佳拟合表面。最佳拟合表面应满足两个要求：一是标准偏差应小于给定阈值；二是选定表面的标准偏差应小于高阶表面。

1.4 叠差变形提取

滑窗算法的主要原理是通过设定固定的窗口沿着一定的方向滑动，窗口可以是圆形或者矩形，并给窗口设置一定的阈值，进而达到对数据的计算限制在各个窗口内进行单独分析的目的。为保证变形的连续性，采取窗口中心间距小于窗口阈值(即相邻插值点所取区域部分重合)的滑窗模式即b模式进行拱圈线形提取。具体步骤如下。

步骤1将1.3节求得的f(xi,yi)及点云导入本算法。

步骤2首先规定起始端为横坐标最小值xmin，以及末端为最大值xmax。

步骤3根据拱圈线形精度要求选取x轴间距Δx、y轴间距Δy及插值点的分布位置(xi,yi)，这些参数满足式(7)。

d为滑窗算法窗口长度；h为滑窗算法窗口宽度图3 拱圈叠差变形提取简图Fig.3 Schematic diagram of arch ring overlap deformation extraction

(7)

式(7)中：ymin为纵坐标最小值。

步骤4对最佳拟合表面提取的插值点进行同次数拟合，得到f(xi,zi)，并绕y轴进行坐标旋转变化。

(8)

式(8)中：夹角φi=arctanf′(xi,zi)；x′、y′、z′为旋转变化后的坐标值。

步骤5视拱圈底部点云数据为总点集C，在总点集C中搜寻满足响应窗口阈值的点集Ci，在算法程序中将点集Ci设定为元胞数组，点集Ci满足式(9)。

(9)

步骤6通过算法搜寻后能够得到各插值点对应的元胞数组Ci，能够针对性地对各插值点的数据进行分析。如果采用该元胞数组的平均点值代表其高程，则当点云缺失时会产生较大的偏差，如图3所示。采用加权质心点表示该元胞数组的高程，则有

(10)

式(10)中：mk为窗口阈值的点集Ci对应点的权重，用该点到点集Ci平均点距离的倒数来计算；k为该区域内点的个数；n为点集Ci中点的个数；Δi为插值点的间隔步长与间隔个数之积；zk为该区域内各点的z轴方向数值；z′i为加权质心点表示该元胞数组的的高程。

然后是对窗口滑动的路径连接所提取的若干个离散的插值点进行加权质心点的差值计算，即完成叠差分析。最后进行逆变换，得到拱圈整体精确线形变形节点数据，并进行3次样条曲线插值获得拱桥变形曲线。

2 实验验证

2.1 拱圈模型点云获取及预处理

对某钢筋混凝土拱圈模型进行单点加载实验，其中加载的具体位置为模型的拱顶中心处，以探讨三维激光扫描仪在桥梁变形监测方面中的应用效果。该实验系统由模型、三维激光扫描仪(3D laser scanner，3DLS)、百分表(dial indicator，DI)等组成。该模型是由某拱桥拱圈按35∶1比例缩小进行制作，值得注意的是，受试验场地的限制，拱圈底座无法实验固结。该实验方案最终确定，采用横系梁连接的方式对两侧底座完成固定。如图4所示，将该模型完成百分表位移测点(①～⑨)的布置。

该拱圈模型计算跨径L为350 cm，计算矢高S为90 cm，矢跨比S/L为0.173。底部横系梁横截面为30 cm×30 cm的正方形。采用徕卡P50扫描仪获取不同工况下的拱圈点云数据，如图5所示。由于同一时期能完成加载，故同时期不动站获取的点云数据处于同一坐标系，无需拼站。如果拼站，则需采用ICP原理进行两期或多期的点云拼接(具体原理见1.2节)，且拼站精度，即均方误差εmin需达到一定值。采用三棱锥法进一步对点云进行去噪，如图6所示。并将拱圈底部点云沿着横桥向划分成3个等分区域，目标为提取拱圈的上、中、下3条线形。

①～⑨为百分表测点图4 百分表测点布置Fig.4 Measuring point arrangement of dial indicator

图5 拱圈点云Fig.5 Arch point cloud

图6 拱圈底部点云降噪Fig.6 Noise reduction of point cloud at the bottom of arch ring

2.2 实验数据处理

2.2.1 拱圈多项式曲面拟合

通过编程的方式实现点云数据多项式曲面拟合的处理。经试验得知，将多项式曲面拟合标准差设置为0.15 mm。结果表明：在对区域1、区域2、区域3拟合曲面时，5次多项次是最优选择，具体如表1所示。图7为区域1的多项式曲面拟合效果，整体上拟合曲面与点云贴合良好，且局部也能保证足够的准确度。后续将基于该次数的多项式曲面曲率信息进行径向滑窗算法叠差变形提取。

表1 拟合标准差Table 1 Fitting standard deviation

2.2.2 拱圈叠差变形提取

将步骤1～步骤6进行对应的编程处理，其中，期望值之差为每个插值点的变形值。对于遮挡区域，如立柱等部分，则采用3次样条曲线对遮挡区域前后点进行插值连接。分别将3个区域的变形进行提取，并与百分表所测数据进行对比分析。

2.3 实验结果分析

如图8所示，由于百分表位置与拱圈点云上边线一侧，故二者的变形十分接近。整体来看，拱圈点云的上、中、下3条线的变形趋势与百分表基本一致。由于百分表分布间距，定位点之间采用三次样条曲线拟合，结果表明：百分表线形与点云3条线形在较多区域偏差在为0.05～0.1 mm，尤其拱圈的反弯点位置。点云3条变形之间的相互差值绝大多在0.05 mm内，且差值并不等同于变形提取的实际偏差，其中还包括拱不均匀变形、扭转等因素。

图7 区域1最佳拟合表面Fig.7 Region one best fit surface

图8 实验室拱圈点云与百分表变形曲线图Fig.8 Deformation curve of point cloud of laboratory arch ring and dial indicator

值得注意的是，鉴于实验室良好的环境，以及距离、点云拼站等对结果影响可忽略不计的情况下，得到上述结果。在实际工程中需要给定刚体变换Rk和tk下源点云与目标点云均方根误差εk(具体见1.2节)，保证毫米级的误差范围的基础上，才能进一步进行相关结果分析。

3 实桥验证

3.1 拱桥点云获取及处理

某山区大跨径劲性骨架拱桥跨径为155 m，拱圈混凝土采用分环分段浇筑工艺，共分为三环十一段，施工工序繁多、复杂。后续在拱圈上浇筑拱上立柱，并搭设预应力简支梁作为拱桥桥面。拱圈的线形作为拱圈混凝土及拱上立柱浇筑的关键特征，直接影响该桥的施工安全及成桥质量，故采用基于点云的劲性骨架拱桥自动化线形监测方法对该拱桥进行研究。

在拱圈浇筑成型后，拱上立柱搭建完成、以及桥面系简支梁安装完成都是拱桥施工中关键的施工状态。故对在3次不同时期(拱圈浇筑成型、拱上立柱搭建完成、桥面系简支梁安装完成)分别对拱桥进行点云数据获取，并以拱圈浇筑成型时为参考，另两期数据分别与其采用ICP算法配准。配准结果为均方误差εmin。拱桥旁边的引桥桥墩在拱桥施工中可认为是稳定的，将此部分共同点云作为点云拼站的基础数据(图9)。其中，拱上立柱搭建完成时与拱圈浇筑成型点云配准的均方误差εmin=0.023 mm，桥面系简支梁安装完成与拱圈浇筑成型点云配准的均方误差εmin=0.025 mm。配准的点云大面积高度重合，说明点云配准精度较高，降低了配准误差对桥梁挠度分析的影响。

对不同期拱圈底部点云沿着横桥向进行等间距的区域划分，分别为上中下三部分区域点云。后续的点云数据处理步骤同内容1.3节、1.4节，在此不再赘述。

图9 拱圈成型与桥面系完成两期点云拼站Fig.9 Arch ring forming and bridge deck system completion of two phases of point cloud splicing station

3.2 结果与分析

由于现场施工质量等因素影响，难以保证与目标线形完全一致，故实际线形走势更为复杂。如图10所示，拱上立柱搭建完成、桥面系简支梁安装完成两个重要阶段的拱圈底部点云3个区域提取的变形趋势基本吻合，全局精度可保证2 mm范围内。另外，该精度还包括拱桥施工发生不均匀变形、扭转等因素影响。相对于本文算法提取的点云变形曲线，由全站仪测量拱圈表贴反光片所定位点拟合而成的三次样条变形曲线局部精度稍差，尤其在拱圈反弯点、以及拱顶等较多区域可达3 mm及以上，这将拱桥的施工质量及后续的变形监测带来挑战。

鉴于实桥线形复杂的特点，采用点云变形色谱图进行算法精度的验证。重点选取桥面系简支梁安装完成阶段的点云变形曲线和全站仪定位点插值曲线差别较大的区域进行分析，如图拱圈A、B、C 3处横坐标分别为-36.0、3.5、40 m对应的上中下线形变形值。在点云处理平台Geomagic Control中完成两期局部变形分析，如图10所示。可以看出，变形值都更接近于前者，进一步说明所提算法的精度及可靠性。

A1、A2、A3分别为上、中、下部分在A点处对应的变形值；B1、B2、B3分别为上、中、下部分在B点处对应的变形值；C1、C2、C3分别为 上、中、下部分在C点处对应的变形值图10 现场拱圈点云与全站仪变形曲线Fig.10 Deformation curve of on-site arch circle point cloud and total station

总体来看，点云变形曲线和全站仪定位点插值曲线都能对拱圈的变形提供基本的保障，考虑到表贴反光片有效寿命以及更高精度变形曲线对于其表贴密度的要求，建议将三维激光扫描获取的拱桥点云数据作为主要参考对象。

4 结论

提出一套基于点云的拱桥自动化高精度线形监测算法，并将该算法在实验室的拱圈模型中进行了精度验证。进一步的，将该算法应用于实桥，同步采用全站仪测量的数据，以及点云变形色谱图进一步保证该方法的精度及可靠性。提取的多条变形曲线对于拱桥尤其是大跨度拱桥的分环分段浇筑及拱上立柱等关键施工工序，以及拱桥后期运维具有潜在的重要意义。

随着工程领域数字技术的不断突破和发展，全面、高精度的数据作为工程的重要保障将变得越来越不可或缺，对应结构的关键特征会更丰富、精确地被识别与提取。同时，该类数据的处理对行业人员提出了更高的要求，在数字化、智能化的今天，技术提升及转型成为必然趋势。先进数据获取技术及自动化算法的研究与应用将为智慧交通带来更多助力。