Kriging点云滤波改进算法及监测试验研究

2022-11-19胡达黎永索张可能梁晓东梁小强吴有平

土木与环境工程学报 2022年6期

胡达,黎永索,张可能,梁晓东,梁小强,吴有平

（1. 湖南城市学院城市地下基础设施结构安全与防灾湖南省工程研究中心；土木工程学院，湖南益阳 413000；2. 中南大学 a.有色金属成矿预测与地质环境监测教育部重点实验室；b. 地球科学与信息物理学院，长沙 410083；3. 湖南联智科技股份有限公司，长沙 410200；4. 中国电建集团中南勘测设计研究院有限公司，长沙 410014）

对隧道围岩的稳定性进行监控及预报是确保隧道现场施工安全、提高施工信息化水平、优化围岩支护设计参数必不可少的一项重要工作。目前，隧道监控量测的主要方法仍然是采用精密水准仪、全站仪、智能收敛计等仪器设备进行测量，这些测量工作首先需要在隧道洞壁上及时布设监控测点（三角挂钩、反光片等）。这种传统监控量测方法实际操作起来非常麻烦，效率极低，且误差不容易控制，监测过程中极易造成对围岩稳定性评价的错误判断，使隧道施工安全性降低。

近年来，激光扫描技术逐渐兴起，该技术具有高精度、高效率、高分辨率、全自动数字化采集、数据量丰富等特点[1]。将三维激光扫描技术应用于隧道监控量测是未来发展的一个重要研究方向，学者们在该技术的探索和研究上取得了大量的成果。王涛等[2]基于点云数据预处理问题，采用小波分析格网化滤波方法，并结合南京地铁隧道进行了应用，取得了较好的滤波效果。许度等[3]提出了一种基于三维激光扫描技术的隧道变形非接触原位测试技术，建立了基于扫描信息的围岩结构面表征方法，并在中国锦屏地下深埋硐室群工程进行了应用。吴勇等[4]采用三维激光扫描技术进行运营期隧道检测，为隧道的病害检测及修复提供了一定参考。杜黎明等[5]提出了一种迭代椭圆拟合的方法，对隧道变形点云数据的提取进行了研究。张帆[6]提出了一种“聚合”算法，并应用于隧道管片结构变形可视化分析。徐卫东等[7]以北京地铁7号线为依托，运用三维激光扫描技术对隧道开挖的现状进行了测量。虞伟家[8]采用新型移动式三维激光测量技术，对盾构隧道进行了点云数据采集，并在苏州地铁进行了隧道监测试验。黄帆等[9]对盾构隧道管片变形进行了隧道监测试验，点云数据分析结果表明，三维激光扫描技术优势明显，对地铁施工具有一定的指导意义。上述研究成果表明，目前精度最高的三维激光扫描仪在很多领域的应用上已经可以替代全站仪，能够基本满足隧道监控量测的精度要求。

大量研究表明[10-11]，点云滤波算法成为三维激光扫描精度是否满足应用要求的关键。因此，对点云滤波进行针对性的研究十分必要。结合近年来各种滤波算法不同的理论背景，可以分成基于坡度、曲面拟合、分割、不规则三角网、形态学及机器学习的滤波算法6大类。基于坡度的滤波算法原理简单、易实现，但过度依赖阈值的设定，在地形起伏较大区域不适合海量数据的实时处理，滤波效果不佳[12-13]。基于曲面拟合的滤波算法有赖于插值方法的选择，所采用的多层级迭代方式会受到每一层级滤波结果的影响，容易出现误差传递与累积[14-16]。基于分割的滤波算法其滤波效果过分依赖聚类分割的结果，而点云分割方法的选择也对滤波结果有很大影响[17]。基于不规则三角网的滤波算法需要占用大量内存，并且对低位噪声敏感，极易误判低地势物点[18]。基于形态学的滤波算法原理简单、实现效率高，但此类算法在地形起伏较大区域的稳健性有待提高，如何提高其整体精度将是后续研究的重点[19]。基于机器学习的滤波算法需要大量的训练样本，其样本必须覆盖所有可能的地形特征，对计算机资源要求非常高，取得较好的滤波效果难度较大[20]。

综上所述，以上各种点云滤波算法均各有局限性，如何改进和完善滤波算法，提高算法精度成为监测试验能否成功的关键。由于Kriging插值算法在快速自动生成DEM方面精度高，整体精度评价方法具有可操作性[21]。笔者结合Kriging插值算法，通过现场试验段监测，对点云数据进行提取分析，将试验得出的变形数据与传统方法的量测数据进行对比，以验证其算法精度，为推进三维激光扫描技术在隧道监控量测中的应用提供理论依据。

1 Kriging滤波改进算法

Kriging插值法是一种运用结构分析与变异函数相关理论，在有限空间内针对区域化变量进行最优和无偏估值计算的预测方法，是研究空间变异和进行空间插值的一种线性无偏最优估计算法。Kriging插值法不仅考虑待估点位置与已知数据位置的相互关系，还考虑变量的空间相关性。假设区域化变量不是相互独立的，具有一定的随机性和结构特性，而且满足二阶平稳和内蕴平稳。同时，Kriging方法在数据网格化的过程中考虑了描述对象的空间相关性质，使插值结果更科学、更接近于实际情况，能给出插值的误差（Kriging方差），使插值的可靠程度一目了然[22]。考虑时间因素，将经典的Kriging空间插值进行拓展[23-24]，建立时空Kriging插值模型，并将其应用于点云格网化滤波分析，提高三维激光扫描数据获取的精度。

1.1 空间插值

Kriging插值法的基本数学模型[25]为

式中：Z′(X0)为预测点的估计值；Z(Xi)为预测点邻域内参与预测的参考点的值；λi为Kriging权系数，在无偏性和最小方差性条件下，其值依赖变异函数的计算结果来确定。

式（1）是n个数值的线性组合，Kriging方法的原则就是保证估计量无偏，且在估计方差σ2E最小的前提下，求出n个权值系数λi。

在无偏性条件下，为了使估计方差最小，用拉格朗日乘子法求条件极值，见式（3）。

F是n个权系数和μ的(n+1)元函数，求出F对λi和μ的偏导数，最后得到用半变异函数值表示的Kriging方程组。

式中：yi,j=y(xi,xj)=y(xi-xj)，为半变异函数值。

半变异函数或者半变差函数是从空间统计学中的方差概念演化而来的，将区域变量Z(x)在点x和（x+h）的值Z(x)和Z(x+h)之差的方差的1/2定义为区域变量Z(x)在x轴上的半变差函数，记为γ(h)，即

试验半变异函数计算公式为

式中：N(h)为被增量h分隔的试验数据对的数目。

为了估计区域化变量的未知值，需要将半变量试验函数拟合到相应的理论半变异函数模型中。典型Kriging算法提供的半变异函数模型包括高斯模型、线性模型、球面模型、阻尼正弦模型、指数模型等，其中，原点附近的球形变异函数模型是线性的，指数变异函数模型是直线，高斯变异函数模型是抛物线。考虑到本研究是将其应用到隧道变形监测的滤波算法中，选用球形模型，这样既考虑了储层参数的随机性，又考虑了储层参数的相关性，在满足插值方差最小的条件下，给出最佳线性无偏插值和方差。具体模型表达式为

式中：c0为块金值；c0+c为基台值；c为偏基台值；a为变程；h为滞后距。

滤波算法流程如图1所示。

图1 Kriging算法流程图Fig. 1 Kriging algorithm flow chart

1.2 时空插值

空间Kriging只能估计某一时间上未知区域的变形量，如果要估计任意时刻任意位置的变形量，需要进行时空Kriging插值。在隧道掘进过程中，监测点的布置、数据的采集以及位移值的变化不是空间和时间上的单一变化，而是时间与空间的统一。选用普通Kriging方法实现数据的时空插值[26-27]，计

算公式为

式中：Z*(s0,t0)为时空点(s0,t0)处的估计值；λi为临近观测值Z(si,ti)的加权系数，引入拉格朗日系数μ进行推导可得

式（9）中的加权系数λ的值可通过式（1）得到，继而代入式（8），可得研究区域内任意点的插值估计值。

2 现场监测试验方案

采用德国Z+F IMAGER 5010X三维激光扫描仪对长沙市城市固体废弃物处理场管理生活区隧道工程进行监控量测，并将所量测的数据与传统监控量测方法的数据进行对比研究。德国Z+F IMAGER 5010X三维激光扫描仪数据获取速率超过100万点每秒，50 m处精度最高可达0.8 mm，测程为0.3～187.3 m，是一台高性能、高精度的旗舰扫描仪，非常适合应用于高精度的工业测量和隧道工程测量。该仪器精度信息如表1所示。

表1 Z+F IMAGER 5010X三维激光扫描仪精度信息Table 1 Precision information of Z+F IMAGER 5010X 3D laser scanner

2.1 工程概况

所采用的试验数据来源于长沙市城市固体废弃物处理场管理生活区隧道工程监控量测项目现场。隧道场地位于长沙市望城区桥驿镇黑麋峰固体废弃物处理场西侧一带的山间空地，为望湘岩基的南部地区，以侵蚀构造地貌为主，属低山丘陵地貌，主要表现为低山—丘陵与冲沟形态，山脉脉络清晰，山顶一般较圆滑，个别较尖，山脊大多数开阔舒缓，局部地段狭窄，形成窄陡山脊。隧道起始桩号为K0+195～K0+285，长为90 m。根据勘察结果，该场地内地形起伏大，山坡较为陡峭，场地内分布的地层主要有第四系残坡积层，下伏基岩为燕山晚期花岗岩。

2.2 试验方案

地面三维激光扫描外业实施计划的制定需要根据测量任务、要求以及现场条件决定，主要包括坐标系、扫描仪及其配准靶标和扫描站的选择等。根据实施计划，提前布设好靶标，然后连接相关设备，设置扫描参数（如扫描范围、扫描距离、扫描间隔等），在不同的扫描站上架设扫描仪进行扫描。

2.2.1 控制点的布设与测量在隧道工程监测试验中，由于现场条件限制，控制点只能布设在一侧，如图2所示。因此，在隧道洞口布设了3个自制的平面靶标，如图3所示。作为控制坐标转换点，要求这些点的位置相对固定，并且不易被施工破坏。布设好控制转换点后，用全站仪测得靶标中心在施工控制坐标系下的坐标，用于后续转换。对于站间转换公共点，选择不受入射角影响、反射率高、方便布设的标准球形靶标，如图4所示。将3～4个球形靶标布设在两相邻测站之间，并在两测站上均对其进行精扫。

图2 平面控制点布设图Fig. 2 Layout of plane control points

图3 自制平面靶标Fig. 3 Homemade flat target

2.2.2 数据采集外业测量是实际获取数据的过程，扫描距离、物体表面材质、控制网、靶标测量精度、光斑大小、扫描点间距、点云拼接精度、全反射物质和外界环境等因素均对点云数据的质量产生直接影响。与其他光学电子仪器一样，外界环境温度、气压、空气质量等因素对激光回波信号的影响显著。因此，鉴于隧道内施工环境异常复杂，为避免数据采集质量受到影响，选取隧道内空气可见度较高的时段进行试验测试。

现场扫描试验前，先在隧道内确定好仪器和靶标的位置，确保在进行扫描时不会被打断或遮挡。布设好靶标球，如图4所示。架好扫描仪，设置好扫描参数，然后进行扫描作业。扫描作业时靶标球需要保持固定状态，避免工人触摸或误动，从而致使球的位置发生变化，影响后续点云配准精度。试验扫描目标为隧道衬砌外表面的目标点，所以不必进行全景粗扫，而是先直接进行全景目标扫描，再在此基础上对靶标球进行精扫（高精度高质量扫描）。试验时可按图5进行设站和靶标布设（站1～站3为测量站点），具体要求如下：两站之间间距约25～30 m，保证两站之间点云重合数不少于30%；靶标球或靶标板摆设与两站中间至两边站的距离比较均匀；靶标尽量上下左右起伏较大，保证空间位置形状良好；两站之间至少有4个靶标，便于后期拼站的精度。图6为三维激光扫描仪在隧道现场的扫描示意图。

图4 标准球形靶标Fig. 4 Standard spherical target

图5 扫描测站及球形靶标位置分布图Fig. 5 Location distribution of scanning station and spherical target

图6 隧道扫描示意图Fig. 6 Schematic diagram of tunnel scanning

2.3 数据处理

2.3.1 数据处理基本流程对于三维激光扫描仪获得的点云数据，并非每一个都能使用，过多的数据点会降低计算机运行的效率、增加存储空间。隧道点云整体图如图7所示。为避免上述问题，需要删除部分数据点，即对点云数据进行精简处理，点云数据处理的主要流程为：

图7 隧道点云整体图Fig. 7 Overall drawing of tunnel point cloud

1）点云编辑。裁减掉点云粗差以及与目标无关的点云。

2）点云配准与拼接。利用公共点坐标将不同测站测量的点云数据转换到同一坐标系中，实现点云的拼接，形成一个整体。目前，点云配准的方法主要有3种：人工靶标配准、点云自身配准和控制坐标系配准。配准后的点云图如图8所示。

图8 配准后的两个测站点云图Fig. 8 Cloud map of two stations after registration

3）建立拓扑关系。点云通常是孤立的点，每个点只与其一定范围内的周围点相关。空间点云拓扑关系的建立主要有八叉树法、网格法和k-d树法等。

4）点云数据的精简。点云数据精简算法大致可以分为4类：包围盒法、随机采样法、曲率采样法和均匀网格法。

5）点云数据的去噪与光顺。由于被测对象表面的粗糙度、波纹等缺陷和测量系统本身的影响，真实数据中混有噪声点。可根据点云质量和后续建模要求灵活选择合适的滤波算法。

6）孔洞的修补。激光扫描的过程中会因为各种原因（如局部遮挡）造成漏测，从而形成点云孔洞。修补算法可以分为2种：一种是在三角网格表面重建时进行三角网格面修补；另一种是先对散射点云数据孔洞进行修补，然后对三角网格面进行重建。

7）点云数据的分割。数据分割形成的不同曲面类型的子区域具有特征单一、凹凸一致的特点。对每一子区域进行单独重构有利于曲面拟合时减小误差和保持点云性质。目前，点云分割的算法主要有基于边的算法、基于面的算法和基于聚类的算法。

8）点云数据的三维建模。在三维模型的建构过程中，曲面重构是最关键也最复杂的一步。目前，主要有两类曲面重构方案：以三角网格面为基础的自由曲面重构方案和以样条曲面为基础的自由曲面重构方案。

2.3.2 点云数据处理点云滤波是点云数据处理中一个非常重要的环节，准确的滤波结果将有利于提高点云后处理结果的精度。根据普通Kriging插值的原理，采用Matlab程序实现Kriging滤波算法，并利用程序对点云数据进行处理。为减少滤波算法的工作量，选取最为接近垂直于靶标球表面的3个测点进行计算，取c0=0.01、c=0.2、a=20。取n=3，根据式（9），有μ=0.015 9、σ2E=0.054 7。

1）将外业扫描获取的点云数据导入到数据处理程序中，进行点云数据预处理操作（滤波、去噪），为后续操作做准备。

2）数据预处理完成后，在配准流程（Registering）中，分别用平面靶标拟合工具和球形靶标拟合工具对靶标进行拟合，提取出靶标中心在扫描仪坐标系中的三维坐标。

3）靶标中心坐标提取完成后，再进行点云配准，将全站仪测得的控制靶标的坐标保存在TXT文件中，并将点云数据转换到控制坐标系中。如果布设的公共点较多，配准时应该选择质量较好的点进行点云配准。最后将配准后的点云数据坐标信息导出到TXT文件中，完成数据处理流程。

3 数据提取及分析

3.1 监测信息提取

根据隧道表面的点云数据坐标，可以提取出任意桩号处的隧道断面数据，从而进行断面拟合，提取出中心坐标，然后拟合出中线。在这些数据上可以对超欠挖情况、拱顶沉降情况、收敛状态和轴线偏差等进行分析，对施工进行指导，及时发现不稳定及危险情况并上报，减少灾害发生及其不良后果的产生。

3.1.1 断面提取在隧道监控量测中，所有的测量工作都是在断面的基础上完成的，所以，三维激光扫描数据应用在隧道中的第一步就是断面提取，否则其他的监控量测项目将无法完成。可以根据隧道的设计资料来提取断面。由断面的定义可知，隧道中轴线上某点处的断面与该点处的切线垂直，即该点处的切线向量就是该断面的法向量。设该点在中轴线上的坐标为(x0,y0,z0)，该点处的切线向量为n=(A,B,C)，则该点处的断面方程为

中轴线上某点的三维坐标可以根据隧道段的线路设计参数确定，该点处的切线向量n=(1,k,i/cos(atank))，k为xoy平面中的切线斜率，i为坡度。

不管扫描仪的分辨率有多高，测量的点与点之间总会有一定的间隔，所以，在实际提取断面时，该断面并不是严格的平面，而是有一定厚度的。在实践应用中，一般提取1～2 cm厚度平面的数据作为断面数据。图9～图11为利用该方法提取的同一桩号处的隧道断面。图9为在施工坐标系中提取的断面；图10为施工坐标系中旋转到正面后的断面，图11为在独立坐标系中的断面，原点在轴线中心。

图9 三维空间坐标系中提取的断面Fig. 9 Extracted section in 3D coordinate system

图10 三维空间坐标系中提取的断面(横截面)Fig. 10 Extracted section in 3D coordinate system(cross section)

图11 平面坐标系中提取的断面Fig. 11 Extracted section in plane coordinate system

3.1.2 变形信息提取依据行业规范《公路隧道施工技术规范》（JTG/T 3660—2020），在隧道监测中，周边位移和拱顶沉降是必测项目，并且每5～50 m需要测量一个断面。三维激光扫描仪扫描的数据是连续、全面的，所以用扫描数据进行拱顶沉降和收敛变化的监测可以在任意断面处进行。而传统监控量测方法对隧道周边的收敛量测一般采用隧道净空变化测定计（简称收敛计）进行。拱顶沉降和收敛变化都是在断面上选择固定的点，对其进行连续的监测，收敛量测测点与拱顶下沉测点布置在同一断面。传统方法需要使用预埋件，而三维激光扫描方法不需要任何预埋件。根据隧道情况，拱顶沉降监测一般选取1～3个点，收敛变化一般选取2～3对点。

在断面中提取监测点数据时，一般以监测点为中心的面状数据为基础，然后对区域内的点进行距离加权平均求取监测数据。设选取的监测点设计坐标为(x0,y0,z0)，选取区域内的点坐标为(xi,yi,zi),i=1,2,3…n，则监测点的测量坐标为

获取了各个测点的监测数据后，就可以以第一期数据为基础，作出时间变化曲线，以反映各断面的收敛与沉降情况，还可以作出空间变化曲线，以反映隧道的整体变化情况。

3.2 数据分析

在隧道坐标文件的基础上进行数据挖掘与分析。根据隧道表面的点云数据坐标，可以提取出任意桩号处的隧道断面数据，从而进行断面拟合，提取出中心坐标，然后拟合出中线。选择该隧道工程具有代表性的3个断面进行对比分析，分别采用反距离倒权法（IDW）、不规则三角网法（TIN）与克里金法（Kriging）进行对比研究。3个断面在监测周期内的最终累计变形对比见表2、表3，图12～图17为所提取的K0+202、K0+207、K0+212断面的位移收敛示意图。

图12 K0+202断面(2016-09-30)Fig. 12 Section of K0+202 (2016-09-30)

图17 K0+212断面(2016-11-12)Fig. 17 Section of K0+212(2016-11-12)

表2 3种滤波方法对比Table 2 Comparison of three filtering methods

表3 试验结果对比Table 3 Comparison of test results

图13 K0+202断面(2016-10-29)Fig. 13 Section of K0+202(2016-10-29)

图14 K0+207断面(2016-10-07)Fig. 14 Section of K0+207(2016-10-07)

图15 K0+207断面(2016-11-06)Fig. 15 Section of K0+207(2016-11-06)

图16 K0+212断面(2016-10-13)Fig. 16 Section of K0+212(2016-10-13)

从表2中的对比结果来看，由于三维激光扫描获取了海量且较为均匀的监测数据，从而对IDW和Kriging方法十分有利，所得结果较相近；而TIN法的滤波结果误差相对较大。总体情况表明，相对于IDW法，Kriging改进算法的滤波更接近真实值，算法精度更高，且更适合数据较多的情况，同时也说明该方法是行之有效的。

由表3可知，三维激光扫描量测的累计变形总体相对较小，其中K0+207、K0+212拱顶沉降数据较为接近传统量测结果，周边收敛数据差异较大。结合现场情况并对设备精度进行分析可知，隧道内环境因素对量测精度影响极大。对于拱顶沉降，扫描仪只需提取Z轴方向的坐标点数据就可以得出沉降结果，因而比水平坐标的换算更为稳定；对于周边收敛，扫描仪必需提取X、Y、Z轴3个方向的坐标点数据才能有效得出收敛变形结果，由此造成的误差较大。除数据提取外，隧道内环境等因素也对监测结果造成较大的不确定性，进而造成误差累积。试验数据同时也说明，在隧道变形监测中，三维激光扫描技术的抗干扰能力和精度还有待进一步优化和提高。在特定的监测环境下，三维激光扫描技术所量测的隧道累计变形值在一定程度上能较好地反映隧道变形的真实情况。