塔形建筑物轮廓点Z坐标点云倾斜监测方法

屈乾龙，朱庆伟，艾卫涛，温 波，马肇祥

(1.西安科技大学 测绘科学与技术学院，陕西 西安 710054；2.河南城建学院 测绘与城市空间信息学院，河南 平顶山 467000)

0 引 言

随着时间的推移，古代高塔如今都是“十塔九斜”，塔体受基础和结构本身的影响及诸多外界因素而产生弯曲变形，综合表现为整个建筑体倾斜变形[1-3]。倾斜度是衡量施工技术质量和后期建筑物安全维护的重要指标，如何快速准确测量塔类建筑物倾斜度一直都是变形观测领域的研究热点。三维激光扫描技术以其无需直接接触即可快速获取物体表面高精度三维点云数据的特点被广泛应用于塔类建筑物的变形监测[4-7]。国内对塔形建筑物的变形监测研究从未停止过，梁华等基于重庆某铁塔点云数据，对铁塔结构底部特征进行圆柱拟合获取铁塔中心坐标，再分别提取关键部位点云，利用特征拟合获取几何参数，计算铁塔上、中、下端中心到铁塔轴线偏差距离[8]。辛星等基于三维激光扫描技术对塔形构筑物进行倾斜监测，利用Imageware软件对点云数据提取构筑物中轴线，计算出构筑物的倾斜值。最后与全站仪前方交会精度实验对比验证三维激光扫描技术在塔型构筑物倾斜监测中应用的可靠性[9]。王莉等以铜川市延昌塔为研究对象，采用三维激光扫描技术获取塔体点云数据，利用Geomagic软件进行建模，然后利用MATLAB软件对点云数据拟合计算，得到塔体中轴线，继而分析判断塔体的倾斜变化状况[10]。杨永林等以西安万寿寺塔为例，基于三维激光扫描技术对塔体进行变形监测，提出最小二乘塔体倾斜度计算方法。该方法通过提取每层塔体的轮廓点云，根据最小二乘原理拟合出每层塔体的圆心，对每层圆心进行线性拟合提取塔体中轴线，继而计算出塔体的倾斜量[11]。

以上学者在采用三维激光扫描技术做倾斜监测时，无论是横断面提取中心轴还是提取轮廓点拟合中心轴，均采用拟合中心轴的方法，数据处理方案较为单一，且均是将塔体看作一个刚体计算整体倾斜量，不能逐层计算。笔者提出一种点云轮廓点Z坐标倾斜监测数据处理方案，与提取中心轴法共同对西安大雁塔点云数据进行倾斜量计算，最后将计算的倾斜量进行对比。计算结果显示，点云数据处理方法对塔体建筑物倾斜监测具有可行性和工程实际意义。

1 研究对象概括

大雁塔位于中国陕西省西安市南郊大慈恩寺内，大雁塔是砖仿木结构的四方形楼阁式砖塔，由塔基、塔身、塔刹组成，现通高为64.517 m。塔基高4.200 m，南北约48.700 m，东西45.700 m；塔体呈方锥形，平面呈正方形，底边长为25.500 m，塔身高59.900 m，塔刹高4.870 m[12]。据记载，陕西测绘部门自1985年开始对其进行测量，其时大雁塔的倾斜速度处于加快过程中，到1996年左右塔倾斜达到最大程度，倾斜度达到1.011 m。如今，具有1 300多年历史的大雁塔整体属于动态平衡。大雁塔历史悠久且倾斜明显，故以大雁塔为例验证文中方法的可行性，对大雁塔的安全运营及维修加固具有一定的参考意义。

2 点云数据获取与处理

2.1 数据来源

三维激光扫描数据采集的误差来源大致分为3种：与目标物反射面有关的误差、仪器误差、外界环境条件。王俊杰等通过以Trimble GX三维激光扫描仪为例研究数据采集误差对数据精度的影响[13]。最终实验得知测距误差对点位精度影响微乎其微；由于扫描角的误差，点位误差随测量距离的增大而增大，可通过点位精度估计公式来计算出扫描最大点位误差，根据实际工作要求的点位精度来选取适当的扫描距离。文中数据是于2016年所测大雁塔点云数据。该数据是由Scan Station C10扫描仪采集的点云数据。设站情况如图1所示，站点分布在塔的四周。近距离站点共8站，主要是获取塔的中下部点云数据；远距离的站点共4站，主要是获取塔的顶部点云数据[12]。数据采集时，X轴在横向扫描面内，Y轴在横向扫描面内与X轴垂直，Z轴与横向扫描面垂直。文中所涉及的所有Z值均是指Z轴坐标值。

图1 设站分布

2.2 数据预处理

在获得原始点云数据后，第一步要对多站扫描数据进行拼接配准处理。由于三维激光扫描过程中难免会扫描到塔体周围各种地物而产生巨量的噪声点，为了防止拼接过程中因为海量点云数据而导致软件崩溃以及提高各站点云数据拼接速度，在原始数据导入Leica Cyclone软件后，对每站除待测塔体和标靶球区域外大部分噪声点云，直接在Cyclone软件中进行区域裁剪，快速除去大量噪声点云，这样后续点云数据拼接和融合速度大大提高[14-16]。然后逐站标记同名点和标靶球，通过标靶球和同名点混合拼接方法进行数据拼接。为了减少个别标靶球或同名点的误差传播，采用每2站数据拼接一次，然后再把拼接好的结果进行拼接[17]，严格把控各站拼接精度，对于拼接误差较大的测站，将该站标靶球和同名点删除后重新标记同名点再进行拼接，直至拼接误差在0.006 m左右。最终将各站基于仪器内部独立坐标系统的点云数据转化到同一坐标系中，形成一个整体。

拼接成一个整体后，仍存在大量的重复点云，所以需要通过设置平均采样间隔或其它参数设置来进行点云数据融合处理，去除大量重复冗余点云[18-20]。

3 倾斜监测方法

3.1 已有三维激光扫描倾斜监测方法

文献[9]提出基于中轴线上节点坐标偏移的方法，其思路是提取中轴线上节点的坐标[21-22]。基本流程如图2所示。

图2 切片中心轴倾斜监测基本流程

文献[10]利用Geomagic软件进行建模，构建OBJ三维格网模型；再通过软件提取每层塔体特征点坐标，计算出该层中心点坐标，利用 MATLAB软件拟合计算塔的中心轴线，最终计算分析塔体倾斜变化状况，基本原理如图3所示。

图3 拟合中心轴计算原理

文献[11]提出基于最小二乘的塔体中轴线提取方法，其思路是首先提取出每层塔体外轮廓点云的圆心坐标，然后基于最小二乘原理对所有层圆心坐标进行线性拟合获取塔体中轴线的方程，最后计算2期塔体的倾斜度，实现塔体的倾斜监测。

以上方法无论是对塔体进行区域分割还是基于塔体塔檐提取特征轮廓点，最后都是基于点云数据提取塔体中心轴的方法来确定塔体倾斜度，文中所提出的方法与该类方法有所不同，不用提取中心点即可计算出塔体倾斜度。

3.2 轮廓点Z坐标倾斜监测方法

笔者提出基于轮廓点云Z坐标提取计算法进行倾斜监测。基本思路是拟合提取塔体固定轮廓特征点坐标，比如每一层的塔檐。如图4所示，提取出每层塔檐外轮廓点云的坐标，A，B，C，D是塔檐的4个角，我们只需在提取出的外轮廓点云坐标中筛选出其Z坐标值最大点D和Z坐标值最小点B，此时LDO为D、B这2点Z坐标之差ΔZ，LDB可根据2点间距离公式求出，最后利用正弦定理公式(1)和反三角函数公式(2)即可算出该2点连线LDB与水平面夹角∠OBD，该角在数值上等于塔体倾斜角。

图4 基于轮廓点Z坐标倾斜监测法原理

(1)

又∠DOB=90°，故公式(1)可化成

(2)

如图5所示，这种是塔体倾斜特殊情况，塔体塔檐连线严格向另一侧倾斜，这时LAD，LBC线上各点的Z坐标总体相同，此时可直接选取同侧塔角点A，B，或者点D，C为Z坐标最大值和最小值点计算∠OEF。

图5 塔体倾斜特殊情况

本方法同样适用于圆形或六边形等规则塔体建筑，原理相同，在此不再赘述。

4 实验数据计算对比分析

4.1 拟合提取中心轴线倾斜监测法

拟合提取中心轴线倾斜监测法具体操作在此不再详细介绍，基本操作流程如图6所示。用Leica Cyclone软件将点云数据预处理后，将拼接后的点云数据导入Imageware软件中进行区域等距离切割点云，并将切割后的点云轮廓数据提取出来[23-25]，为了提高精度，准确计算出塔体各中心坐标，因此首先将轮廓点云拟合成等间距圆，如图7所示，共提取7条轮廓线。

图6 中心轴法实验流程

图7 切割轮廓点提取拟合圆

再计算出各层中心点圆心坐标并将计算得出的坐标和轮廓点坐标数据导入软件具体效果如图8所示。

图8 顶点及各层中心点效果

计算出各层中心点坐标导入MATLAB软件显示后明显发现各中心点与塔顶中心点并不在同一直线上，这是因为塔体建筑倾斜形变并不是线性形变的，因此中心轴法需要将塔体作为一个刚体整体分析。将中心点坐标和塔顶点坐标通过MATLAB编写程序，基于最小二乘拟合原理提取塔体中轴线，结果如图9所示。

图9 塔体中心圆点拟合直线(单位：m)

最终得到拟合后顶点和各层中心点坐标，见表1。

表1 拟合后顶点和各层中心点坐标

最终通过拟合中心坐标计算，该年大雁塔塔体大约整体倾斜0.960 m。

4.2 基于轮廓点Z坐标倾斜监测法

流程如图10所示。在点云拼接、去噪后，为了保证轮廓点提取精度，采用区域分割分层提取檐轮点云塔廓线，区域分割时，只需保证塔檐整体部分在分割区域内即可，如图11所示。

图10 轮廓点Z坐标提取法实验流程

图11 区域分割

因为在数据融合时为增加后期数据处理速度，故设置采样间隔进行点云抽稀，为防止顶点和塔檐四角点云误删而导致Z坐标值最大最小点提取不准确，故在Cloudcompare软件中选取各层塔檐四角和顶点坐标加入轮廓线提取坐标群，如图12所示。

图12 选取的塔角四点及塔顶

根据该方法进行轮廓点提取，最终各层塔檐共提取点云数量见表2。

表2 各层塔檐轮廓提取点云数量

提取各层轮廓点点云坐标后，导入表格从大到小排序，将Z坐标值最大最小点坐标单独导出，根据公式(1)、(2)，求得各层倾斜角度以及基于该层倾斜度乘以塔的通高所求的塔整体倾斜量xi，详细数据见表3。

表3 各层倾斜度及基于该层倾斜度的整体倾斜量

在求塔体整体倾斜量时，对表3中基于各层倾斜度整体倾斜量直接取算术平均值不能满足测量中的精度要求，为了提高塔体整体倾斜计算精度，引入权要素，对表3中整体倾斜量做加权处理。对于塔形建筑物而言，底层较顶层相对倾斜量较小，且对整体倾斜影响最大，由于激光扫描点精度相同，且轮廓点提取点云密度直接影响能否提取出真正Z坐标值最大最小点，每一个轮廓点都相当于同精度观测一次，因此用各层塔檐轮廓提取点数作为Ni次同精度观测值参与定权计算。7次相同精度观测塔体求解塔体倾斜量的过程，故可求同精度倾斜值算术平均值的权，具体公式见下式

(3)

(4)

(5)

式中Pi为第i层倾斜量的权；Ni为第i层轮廓点提取点云数，数据见表2；σ点为每一个点云的中误差，每层测量精度相同，其中误差为σ0，σi为第i层倾斜量中误差。

将公式(4)(5)导入公式(3)中得

(6)

式中C为单位权轮廓点的点数，个。

取第一层轮廓点提取数做单位权观测值，即C=2054。则根据公式(6)，可计算出每层倾斜量的权Pi。

(7)

由公式(7)可得倾斜量加权平均值即塔体整体倾斜量，结果为0.962 m。

4.3 数据分析与对比

通过提取的轮廓点最大值最小值点坐标数据对比后发现，文中方法所提取的Z坐标值最大、最小值点均为塔檐对角点，可根据文中方法原理可大致判断出塔体倾斜方向为东南向西北方向倾斜。基于文中方法，在能判断塔体整体大致倾斜走向情况下，可直接从塔檐各角坐标中提取Z坐标值最大、最小点计算，省去提取轮廓点坐标这一步骤，减少计算量，可大大提高数据处理效率。为了验证该规律的可靠性，在明确塔体大致倾斜方向后，随机提取某层塔柱沿该倾斜方向对角点来进行验证。最终，对塔体第3层和第4层塔柱的两对角点进行提取坐标计算塔体的倾斜量以此来验证该规律的可靠性，具体数据见表4。同样，和提取中心轴法相比文中方法原理易懂，数据计算处理更为简便。

根据表5各层倾斜度数据可以发现，每层倾斜角度并不相同，塔体底层倾斜量较小，随着塔体升高各层倾斜量不一，这也是提取中心轴法不便将塔体分层计算倾斜量原因，塔各层的变形并不是线性的，采用不同层数据进行计算整体倾斜量会有不同的结果。因此现有基于点云数据提取中心轴法进行塔形建筑物倾斜监测分析时，均是将塔体假设为一个刚体进行分析，利用塔体中心点数据拟合出塔体的中轴线数据再计算塔体倾斜量。而文中提出的方法是通过分层计算塔体倾斜角度，再基于该倾斜角度进行整体偏移量的加权计算，既可计算塔整体倾斜量，亦可计算各层倾斜量。通过表4可以看出，无论是从塔檐提取轮廓点还是从塔柱提取轮廓点，都可以满足文中轮廓点Z坐标倾斜监测方法。结合表5看出，2种方法计算塔体倾斜量尽管存在一些差异，所得出的大雁塔整体倾斜量相差0.002 m，约0.26%，但是在测量误差范围之内。总体来说，通过与拟合中心轴法整体倾斜量计算对比表明，文中所提的基于轮廓点Z坐标的塔形建筑物倾斜监测方法切实可行。

表4 数据特征验证

表5 两方法倾斜量计算结果对比

5 结 论

1)对大雁塔三维激光扫描数据进行研究计算，验证文中所提出的基于轮廓点Z坐标倾斜监测法的有效性，该方法可以对基于三维激光扫描塔形建筑物倾斜监测提供一种新的数据处理方案。

2)采用基于塔体各层倾斜角度加权计算塔形建筑物整体倾斜量的方法，因此可以明确获得塔体各层倾斜度，为塔体某层具体修缮维护提供数据参考。

3)在明确塔体基本倾斜方向后，可以直接提取该倾斜方向各层塔檐对角点云坐标数据，作为Z坐标最大最小值点参与倾斜量计算，减少点云计算量，可大大提高数据处理效率。

4)该方法适用性广泛，圆形、六边形等同类项目，可基于文中方法进行倾斜监测，对无中心轴不规则建筑物倾斜监测项目也具有一定的参考价值。