基于椭圆约束的隧道点云滤波

吕亚磊，李永强，董亚涵，杨清科

(河南理工大学 测绘与国土信息工程学院，河南 焦作 454000)

随着大中城市对地下空间利用进程的加快，对于地下空间的监测成为一个严峻问题。移动LiDAR相对于传统监测方式和地面三维激光扫描仪，具有精度高、速度快、使用简单等优点在空间数据获取中具有广阔的应用前景[1-2]。当前研究中多通过地面三维激光扫描仪获取点云，进行隧道形变监测的研究。在隧道点云中，隧道壁上的附属管道、线路、检修平台等对隧道的监测、模型重建等带来不利影响，是隧道点云数据中的噪点[3]。谢雄耀等通过对隧道断面点云进行圆拟合滤除噪点，获得隧道壁点云[4-6]。然而隧道因施工、形变等因素可能会存在一定形变导致隧道断面不是圆形。琚俏俏等在以隧道截面为椭圆形、拟合出隧道轴线的基础上，在隧道轴线方向插入横截面，然后根据横截面法向量将点云投影到截面上，对截面投影点拟合椭圆，求出标准差，剔除误差较大的点[7-8]。李珵等通过投影法拟合隧道中心轴线，然后通过椭圆柱面拟合进行隧道内壁点云滤波[9-10]。地面三维激光扫描仪在隧道的应用极大增强隧道监测的效率，然而隧道环境一般为狭长结构，采用地面三维激光扫描仪进行获取数据，需要测站间有一定重叠度和设置大量标靶点进行拼接，再通过点云数据计算中心轴线，然后进行断面提取和拟合滤波。地面三维激光扫描仪获取隧道空间点云的方式不仅需要点云拼接等，而且获取断面切片时需要通过点云数据进行拟合中轴线，效率较低。根据隧道点云数据进行拟合中轴线，存在仅能对一段距离内的点云进行提取拟合，而无法对隧道的全局进行拟合的不利因素[11]。魏楚文通过移动测量小车获取隧道数据数据，验证移动LiDAR在隧道环境下的采集数据的可行性[12]。Boavida J等通过设置间隔25 m的平面标靶点校正航迹数据，对长距离隧道进行竣工测量[13]。Arastounia M等通过移动测量小车获取隧道点云数据，以隧道壁法向量指向隧道中轴线的特点拟合隧道中轴线，获取隧道断面并通过最小二乘椭圆拟合进行隧道噪点的滤除[14]。Zlot R等采用3D-SLAM算法对矿山巷道进行三维重建用于大型设备的运输安装[15]。上述算法仍然是通过点云数据拟合中轴线，然后进行隧道点云的滤波处理。

本文通过地下移动LiDAR系统采集地下空间的隧道数据，避免传统地面三维激光扫描仪需要频繁设站、点云拼接等采集步骤，提高隧道空间数据获取效率。以移动测量轨迹作为断面划分依据，避免点云数据计算中轴线的大量工作以及拟合误差。然后根据RANSAC算法对隧道点云数据进行粗滤波，最后通过稳健最小二乘椭圆拟合滤波进行精细滤波，获得准确的隧道壁，对滤波拟合结果进行滤波精度和拟合结果的鲁棒性进行分析，以便进行断面拟合、收敛分析、净空分析等应用。

1 隧道点云滤波方法

隧道为狭长结构，采用移动测量平台搭载LiDAR系统快速的获取隧道空间数据。但是隧道的设计中存在曲线，同时为满足点云拟合的需要，首先根据航迹对隧道点云进行断面切片划分，其次将点云转换到定义的载体坐标系中，获得点云相对于航迹点的相对坐标。然后在载体坐标系下采用RANSCA圆拟合滤波滤除距离隧道壁较远的粗差点云，最后采用稳健最小二乘拟合算法滤除距离隧道较近的粗差点并拟合隧道断面参数，如图1所示。

图1 隧道点云滤波流程

1.1 点云断面划分

隧道点云密度大，直接进行滤波处理较为困难，考虑隧道狭长结构，以及隧道拟合滤波的特点，在处理点云前将采集点云进行划分不同的断面。当前划分隧道点云的断面方法主要通过最小包围盒法、粗切点云拟合法、投影法等提取中轴线，然后进行断面划分。最小包围盒法在使用中主要使用紧密性较好的长方体盒模型进行确定中轴线，此方法在直线型隧道的中轴线提取效果好，但是对于曲线型隧道需要减小步长才能达到较好效果，且点云中存在噪声时需要先去除噪声才能取得较好效果。粗切点云拟合法主要采用手工选取点云边界点云进行拟合中轴线，此方法需要大量人工交互和多次重复才能获得较好结果。根据移动LiDAR系统数据特点，首先通过临近点距离最小原则对点云进行里程赋值，获得里程有序化点云，如图2所示。其次设置断面宽度，获取断面切片，如图3所示。

对于任意点P(xi,yi,zi),首先在航迹点中查找最邻近点Nj(xj,xj,xj)，使

(1)

将航迹点Nj(xj,xj,xj)的里程赋值给P(xi,yi,zi)。然后根据断面宽度和航迹里程便可确定断面切片中心航迹点里程为Sc，对于任意断面切片点集I，则有

I={Si|Sc-d/2≤Si≤Sc+d/2}.

(2)

图2 点云里程赋值示意图

图3 断面划分示意图

1.2 载体坐标系坐标获取

断面切片获取后的点云坐标仍为WGS-84坐标系，且断面切片可以为任意方向。对于隧道点云的断面方向任意问题，当前采用的处理方式有坐标变换法和柱面滤波法。文中以断面切片中心航迹点为原点，构建载体坐标系。载体坐标系中，Z轴指向平台前进方向，Y轴竖直向上，X轴与Y，Z轴构成右手坐标系，如图4所示。点云转换至载体坐标系后，点云向XOY平面投影，便可获得点云在二维平面的投影数据。

图4 载体坐标系定义

1.3 椭圆约束的隧道点云精细滤波

获得的隧道切片经RANSAC圆模型滤波后，距离隧道壁较远的粗差点得到滤除，对于和隧道壁贴合较为紧密的粗差点，采用椭圆约束的稳健特征值法进行粗差点的剔除。

在空间中，椭圆的一般方程可以定义为

F(x,y)=Ax2+Bxy+Cy2+Dx+Ey+1.

(3)

其中：A，B，C，D，E为空间曲线的参数。

对于隧道壁点，在最小二乘准则下满足点到椭圆的距离代数和最小，为计算方便，以点到椭圆曲线的不符值平方最小。因此，

Cy2+Dx+Ey+1)2=min.

(4)

(5)

可求得椭圆曲线方程A，B，C，D，E；对点云椭圆拟合后，引入稳健思想，进行稳健迭代，滤除和拟合结果偏差较大的点。

稳健特征值法椭圆拟合滤波具体算法如下：

1)以特征值法求得的曲线参数A0,B0,C0,D0,E0,F0。

2)利用参数A0,B0,C0,D0,E0,F0，计算点到曲线距离平方di。

3)计算点到曲线的距离di的标准差。

(6)

4)当di>2σ时，此点作为异常点剔除，否则保留当前点。

5)利用剩下的点重新计算椭圆曲线参数。

6)重复2)～5)步骤，直到所有点到曲线的距离di都满足给定阈值，即任意di<2σ。

7)计算最终的曲线参数A，B，C，D,E。

根据椭圆一般方程和标准方程之间的转换关系，可得到最终拟合椭圆的标准方程

(7)

2 应用案例及分析

为验证本文滤波算法的可行性以及滤波效果，选取一处隧道数据进行验证。在数据获取过程中，为保证精度可靠性，在进入隧道前的地面部分，采用GPS+IMU组合导航方式获取精确轨迹，在进入地下隧道中后，通过人工设置已知位置的标靶点参与组合导航解算，校正IMU漂移误差，提高数据位置精度。数据选取扫描隧道的19 m区域作为实验区域，该区域包含176万激光点，点间隔0.02 m。

图5 实验数据及粗差点说明

2.1 隧道断面切片生成

隧道点云数据获取的过程中，获取点云数据的同时获取到航迹数据。航迹数据实时获取扫描仪的位置和姿态。在获取的航迹数据基础上，通过航迹里程映射和设置断面切片厚度后获得隧道的断面切片。其中根据点云间隔为0.02 m设置隧道断面切片厚度为0.20 m。如图6所示。

图6 点云断面分割效果图

获取的隧道断面坐标仍然是WGS-84坐标系下的全局坐标，为能够使用平面拟合算法进行隧道粗差点滤波，将点云计算归化到以断面中心航迹点为原点的相对坐标。将相对坐标向XOY平面投影便得到投影在一个平面的点云数据。

2.2 隧道壁点云滤波

在点云切片获取后，隧道壁点云上的粗差点采用RANSAC算法进行初步的粗差点剔除。RANSAC算法为以迭代随机过程，参数选择对于滤波效果至关重要，根据隧道点云采集时选用的扫描仪精度为0.02 m, 本文考虑处理质量和效率，选择RANSAC算法的误差参数为0.10 m，滤波效果如图7所示。和原始断面点云相比，隧道断面上的接触网、线缆、底板轨道等附属设施大部分得到滤除，由此说明RANSAC算法能够有效去除距离隧道壁较远的粗差点云。

图7 RANSCA滤波前后效果对比图

RANSAC算法能够有效去除粗差点云，但是和隧道壁紧密连接的粗差点云数据并不能充分的去除，采用稳健最小二乘椭圆拟合算法进行粗差滤除过程中，采用k倍中误差作为滤除标准，结果如图8所示，其中图8a、8b、8c分别k=1、k=2、k=3时的滤波效果。从滤波效果图中可以看出，当k=3时，与隧道壁表面距离较近的粗差点云得到有效滤除，同时隧道壁点云结构较为完整。而对于k=2和k=2时粗差点云得到有效去除的同时出现隧道壁点云较多缺失。由此表明，隧道壁的细小附属物点云具有粗差点云的统计学特性，能够采用本文算法进行滤除。

图8 稳健最小二乘滤波效果图

2.3 滤波结果分析评价

对于数据滤波的结果分析，需要有一个标准数据集作为分析依据，本文通过人机交互对每一断面进行隧道壁点云提取，为隧道点云滤波的参考依据。与单独RANSAC算法进行隧道壁滤波方法进行滤波精度对比，并进行定量评价。

对于滤波精度评价，将滤波结果中存在的误差分为两类，I类误差是将隧道壁点云错误的划分为粗差点集中(弃真错误)，II类误差是将粗差点错误的划分为隧道壁点云(存伪错误)。对于任意一断面点云，作为集合P，经过手动提取得到隧道壁点云为集合A，粗差点为集合B。经算法滤波后，隧道壁点云为W，E1为第一类粗差点，E2为第二类粗差点，则：

集合关系

P=A∪B，

A∩B=∅，

E1=W∩B，

E2=W∩A.

I类误差e1和II类误差e2可以定义为

e1=E1/A，

e2=E2/B.

为验证算法的滤波精度，本文选取不同隧道壁完整度的断面A，B，C，在相同滤波效果的前提下进行误差分析，其中RANSAC算法的内点阈值为0.05 m,本文方法中，RANSAC阈值为0.1 m，稳健最小二乘滤波系数为3倍中误差。对应点云滤波前后效果如表1所示，滤波后误差结果如表2—表4所示。

表2、表3统计结果表明，在隧道壁完整或存在一半时，RANSAC算法I类误差较小，而II类误差较大，可见该方法保留隧道壁点被滤除较少，但是存在较多粗差点未被滤除。本文方法中，I类误差大于RANSAC算法，而II类误差小于RANSAC算法，可见在一部分隧道壁点被作为粗差点滤除的情况下，提高粗差点的滤除效果。表4表明，在隧道壁点云存在较少时，本文方法和RANSAC算法滤除效果近似，不具有显著优势。

表1 不同完整度断面滤波效果

表2 点云A精度评定

表3 点云B精度评定

表4 点云C精度评定

对于移动测量系统获取的隧道点云数据，由于地下空间无GNSS信号，绝对定位的位置精度不足，但是点云的相对位置关系仍然是高精度数据。文中采用隧道拟合后的椭圆长短轴、隧道断面的偏心率作为滤波拟合结果几何参数的评价方式。

根据本文算法对实验数据进行滤波拟合处理，最终获取断面切片94个，其中有效断面结果88个。为验证点云提取结果的可靠性，对隧道断面切片拟合结果的轴长、偏心率、轴长分布、轴长中误差进行分析。提取结果的长轴、短轴结果中，在隧道起始处，因隧道点云切片数量较少，提取结果出现较大偏离。后续点云切片结果中，长短轴稳定在2.66 m,如图9所示。拟合椭圆的偏心率在均值0.005浮动，且轴长主要集中在2.66 m,如图10所示。点云拟合结果的长短轴中误差分别为0.000 926 m和0.000 9 m，如表5所示。导致偏差的原因可以大致划分为2个原因，隧道的设计和施工结果并不能完全一致，可能存在施工误差。在隧道运行中存在一些微小形变。在隧道拟合过程中存在计算误差。

图9 隧道断面提取轴长结果

图10 断面提取结果稳定性分析

表5 提取结果误差分析 m

3 结束语

针对隧道点云中粗差点云滤除问题，提出采用RANSAC拟合滤波和基于稳健最小二乘拟合滤波相结合的滤波方法。主要包括移动LiDAR点云断面切片划分方法、RASCAC圆模型拟合滤波和稳健最小二乘椭圆拟合滤波，并通过具体工程实验数据对算法在隧道拟合滤波结果进行分析。对于滤波结果，对于断面完整度超过一半的点云，与RANSAC算法相比，能够更加有效的滤除非隧道壁点，而完整度较低的隧道点云，与RANSAC算法效果相近，并无显著优势。经过滤波后的隧道壁点云，拟合出隧道长短轴和变率的稳定性分析，证明RANSAC和稳健最小二乘椭圆拟合滤波的稳健性。随着地下移动LiDAR在隧道空间的应用，本文算法可以用于隧道模型重建等应用，但是此方法主要是针对圆形隧道进行，对于矩形、马蹄形等隧道内壁粗差点的滤除仍待进一步研究。