摘要
针对三维激光扫描密集点云提取洞室表面变形信息的问题,本文提出一种基于改进的Alpha Shapes算法识别洞室轮廓点云和多尺度模型到模型的点云比对(Multiscale Model-to-Model Cloud Comparison,M3C2)的洞室表面变形监测方法.首先对获取到的两期洞室表面点云数据进行配准,采用改进的Alpha Shapes算法识别洞室表面外轮廓点云.获得的两期洞室表面外轮廓点云经精配准后,再采用M3C2算法进行各点变形值计算,最后进行距离聚类提取连续形变区域.实验结果表明:该方法能够有效剔除点云中细小沟壑处的点及受到混合像元影响的点,在洞室截面到扫描仪距离10 m的范围内,两期点云剔除率分别为14.17%及13.52%,在70 m范围内,分别为6.25%及6.42%;该方法能够准确高效地提取出2倍配准误差以上的洞室表面形变区域.
关键词
Abstract
Aiming at the extraction of cavern surface deformation from three-dimensional laser scanning dense point clouds,we propose a method integrating the Multiscale Model-to-Model Cloud Comparison (M3C2) with an improved Alpha Shapes algorithm.First,the two-phase surface point cloud data are registered,and the improved Alpha Shapes algorithm is used to identify the outer contour point clouds.After the fine registration of these two-phase outer contour point clouds,the M3C2 algorithm calculates the deformation value of each point,and finally the continuous deformation regions are extracted through distance clustering.Experimental results show that the proposed method effectively eliminates the points at small furrows as well as those affected by mixed pixels.Specifically,the removal rates of point clouds in the two phases within 10 m from the scanner to the cavern section are 14.17% and 13.52%,respectively,which are 6.25% and 6.42% within 70 m.This method accurately and efficiently extracts the cavern surface deformation regions with more than twice the registration error.
0 引言
三维激光扫描技术凭借其数据采集速度快、精度高、非接触式以及受环境因素影响小等特点,已被广泛应用于各个领域的变形监测中[1-4].地下洞室工程地质复杂、施工环境恶劣,传统方法依靠经纬仪、水准仪、GPS 接收机等设备获取数据,不仅安全风险高、效率低,且获取的数据为离散信息,存在一定局限性[5-6].因此,将三维激光扫描技术应用于铁路隧道、水工隧洞、人防洞室等地下洞室的变形监测中具有极大优势[7-8].
目前,已有学者在基于三维激光扫描技术的地下洞室变形监测方面做了大量研究工作.研究洞室变形的方法主要有以下两类:第一类是通过研究洞室断面点云来分析洞室整体变形;第二类是通过对洞室点云进行特征提取或三维建模来分析洞室整体变形.第一类方法如:Sun等[9]将隧道点云数据通过圆柱投影扩展为二维图像,从隧道断面的实测数据和设计数据中获得变形值;Yasuda等[10]针对圆形隧道断面,根据变形前后隧道轮廓傅里叶级数表示的差异,来估计隧道变形;王永锋等[11]将配准后的隧道点云切片得到二维断面数据,利用断面的k近邻点拟合曲线提取隧道断面变形值,并采用弧线投影的方法表达变形;周人飞等[12]采用RandLANet分割隧洞点云断面,再分别利用RANSAC和LSM对隧洞中轴线和各断面椭圆曲线进行拟合,计算隧洞的纵向变形量.第二类方法如:Pan等[13]构建了一个综合大型溶洞地上和地下特征的三维模型,以该模型作为研究对象,分析了结构在罕遇地震中的变形响应;Jiang等[14]根据拱的几何特征进行隧道表面点云间的坐标变换,再基于多尺度模型云比较算法完成隧道地表序列变形计算;韦征等[15]对点云数据采用贪婪三角法进行曲面拟合得到法向量,形成法向量概率密度函数以判断隧道的变形情况.
综上所述,现有的基于三维激光扫描技术的地下洞室形变监测方法已较为成熟,但针对局部微量的形变仍缺乏精细提取能力,主要体现在数据处理中忽略了施工期洞室表面平整度低、细小沟壑众多,且密集点云存在“厚度”特征等特点,难以提取洞室表面的真实微小形变.此外,洞室拱顶的大量支护锚杆和水泥柱块,导致点云数据中存在“拖尾噪声”,会引起变形值的误判.因此,在精细的洞室表面形变区域提取中需要消除上述影响.本文引入改进的Alpha Shapes算法识别洞室表面外轮廓点云和多尺度模型到模型的点云比对(Multiscale Model-to-Model Cloud Comparison,M3C2)算法提取洞室表面微小形变量方法,从而实现更精确、完整地反映地下洞室表面的形变提取.
1 基于改进Alpha Shapes算法的轮廓点云识别方法
1.1 改进的Alpha Shapes算法
Alpha Shapes算法于1994年由Edelsbrunner等[16]提出,被广泛应用于三维表面重建.本文采用Alpha Shapes算法识别轮廓点云.为了区分具有“厚度”点云的内外侧,通过限制Alpha Shapes算法中α球在点云表面“滚动”的位置设置方向,保证改进后的算法提取具有“厚度”点云的轮廓均为“外侧”轮廓的点.通过上述具有“方向性”的限制,Alpha Shapes算法可排除密集点云“厚度”及点云中噪声数据的影响,增强变形提取的准确度.
改进算法的主要思想是用一个半径为α的球在点集上“滚动”,若点集不能使α球通过,那么使α球无法通过的3个点构成该点集的一个轮廓三角形,如图1所示.当α球“滚动”至所有可能的轮廓三角形后,即可将所有轮廓三角形的集合确定为整个点集的轮廓面[17-18].

图1Alpha Shapes算法示意
Fig.1Schematic of Alpha Shapes algorithm
设洞室中轴线方向与X轴同向,三维激光扫描仪安置于坐标原点处,将点云数据距离扫描仪近的一侧视为“外侧”,距离洞室墙体近的一侧为“内侧”,如图2所示.将入射激光方向在YOZ面上的投影向量记,轮廓三角形指向内侧的法线记作,如图3所示.令与它的夹角θ小于设定阈值β(此处为90°),即可限制α球滚动位置为“外侧”,获取洞室点云外侧的轮廓三角形,进而确定洞室表面外轮廓面.

图2识别洞室表面外轮廓点云
Fig.2Recognition of point clouds for outer contour of cavern surface

图3α球的位置确定
Fig.3Localization of α-ball
由于洞室拱顶表面裸露的支护锚杆和水泥柱一般呈现垂直于洞室拱顶的局部表面形态,在点云中呈现为长短不一的锥状突起点云簇,如图4所示.各点云簇中还含有一定数量的混合像元噪声点.为消除这些因素对洞室表面变形值提取的影响,进一步将θ角的阈值β设置为60°,使点云簇不被识别为外轮廓三角形的顶点.
计算α球的球心坐标的方法如下:
已知3个点P1(x1,y1,z1)、P2(x2,y2,z2)、P3(x3,y3,z3)坐标及α值,确定过此3个点且半径为α球的球心O.

图4洞室拱顶支护锚杆点云
Fig.4Point cloud of bolt supporting vault of cavern
首先确定过P1、P2、P3三点的平面П方程:
(1)
再求解球心O在平面П上的投影点P(xp,yp,zp)的坐标,列矩阵方程:
(2)
根据式(2),得到P点坐标(xp,yp,zp),计算OP点之间的距离d:
(3)
式中:r为平面П截α球的截面圆的半径.
可求得球心坐标O(xo,yo,zo):
(4)
式中:为轮廓三角形的法向量,.
每组数据均可确定2个球心坐标,分别判断其θ角是否满足阈值,若有一个满足,则P1、P2、P3可确定一个外轮廓三角形.
由于锚杆和水泥柱底部易被识别为外轮廓三角形的顶点,从而出现零星、离散的轮廓点,故需对这些点进行剔除.遍历所有外轮廓三角形,若某外轮廓三角形的3个顶点均不与其他2个及以上三角形共用,则将该轮廓三角形判断为离群三角形并剔除.提取剩余外轮廓三角形的所有顶点,对这些顶点使用半径滤波滤除稀疏离群点,滤波后输出剩余点云,即为最终的洞室表面外轮廓点.
1.2 基于改进Alpha Shapes算法的洞室表面外轮廓点云识别方法
设扫描点云S含有n个三维空间点,利用改进的Alpha Shapes算法识别洞室表面外轮廓点云的步骤如下:
1)设置参数α以及角度阈值β.
2)从点集S中选取一点P1,将与之距离小于2α的所有点构成子集S1.
3)在S1中任取两点P2和P3,求同时过P1、P2、P3三点且半径为α,并且满足角度阈值β的球是否存在.若存在,记球心的坐标为O;若不存在,则重新选取点P2、P3.
4)遍历点集S1,计算球心O到S1中除P1、P2、P3外所有点的距离,记为集合d.若d中任意距离值都大于半径α的值,则将P1、P2、P3构成的三角形记作边界三角形;否则认为P1、P2、P3构成的三角形不是边界三角形.
5)选取S1中的下一组点作为P2、P3,重复步骤3)和4),直至遍历完S1中所有可能的P2、P3组合.
6)选取S中的下一个点作为P1,重复步骤2)—5),直至遍历完S中的所有点.
7)遍历所有边界三角形,记录各个顶点参与构成的轮廓三角形个数m,若某边界三角形3个顶点的m值均小于3,则删除该边界三角形.
8)提取所有边界三角形的顶点构成点集S2,输出点集S2即为洞室表面外轮廓点云.
最后,对得到的外轮廓点云进行半径滤波处理,剔除离散点,得到最终的洞室表面外轮廓点云.
2 基于Alpha Shapes轮廓点云识别的洞室变形监测方法
2.1 基本思路
首先,对于两期洞室表面密集点云数据,利用球形标靶拟合得到的球心作为特征点,进行特征点匹配完成粗配准,对配准后的两期数据,使用改进的Alpha Shapes算法识别洞室表面外轮廓点云.在对两期洞室表面外轮廓点云进行精配准后,使用M3C2算法进行各点变形值计算,最后对形变点通过聚类获得形变区域.变形监测的具体流程如图5所示.
2.2 洞室点云预处理方法
洞室点云预处理可分为两步:
1)基于特征点匹配的粗配准.在三维激光扫描数据采集过程中设置球形标靶,根据靶球表面点云数据,采用基于入射角定权的球形标靶点云拟合方法[19]求取靶球球心坐标,以此作为同名特征点.根据同名点之间的关系求得点云的旋转、平移矩阵,根据矩阵完成点云变换.

图5洞室变形监测方法流程
Fig.5Flow of the proposed cavern deformation monitoring method
2)基于Alpha Shapes算法的洞室表面外轮廓点云识别.Alpha Shapes算法需要反复进行近邻点查找与球心计算,为减少计算量和计算复杂度,先对点云数据进行抽稀以及切片处理.对点云切片采用1.2中的算法识别两期洞室外轮廓点云,外轮廓点云已消除了大部分噪声,比原始数据更准确地反映了洞室表面的真实起伏.
2.3 洞室表面变形监测方法
洞室表面变形监测可分为三步:
1)基于ICP算法的洞室表面外轮廓点云精配准.为保证后续变形监测的精度,采用迭代最近点算法(Iterative Closest Point,ICP)[20]对两期洞室表面外轮廓点云进行精配准.将第1、2期数据分别作为目标点集和待配准点集,计算旋转及平移矩阵,根据矩阵完成点云变换.
2)基于M3C2算法的洞室表面变形值计算.M3C2点云比对算法可直接在点云上计算点对距离,以检测复杂的形变[21].本文采用 M3C2算法计算两期外轮廓点云的距离,主要步骤为选取特征点云、拟合法线、计算距离变化.
首先,将使用改进的Alpha Shapes算法提取的两期洞室表面外轮廓点云作为特征点云.遍历第1期点云,以第1期点云中的每个点作为核心点Pcore,如图6a所示.
然后进行法线拟合.设置搜索半径R,对Pcore半径R内的邻域点,进行最小二乘平面拟合,得到平面法向量.以通过Pcore的法向量为轴线,设置底面半径r构建圆柱体.圆柱体分别与两期点云相交,形成点云柱,将点集A、B中的点分别投影至轴线,得到投影点点集,如图6b所示.
最后进行距离计算.分别计算点集A′、B′的重心,计算M1与M2之间的欧氏距离,作为两期数据在Pcore处的变形值:
(5)

图6M3C2算法示意
Fig.6Schematic of M3C2 algorithm
3)基于距离聚类的洞室表面形变区域提取.为定位和量化形变区域,在计算出各个核心点的变形值DM3C2之后,对发生变形的核心点进行距离聚类,以达到提取洞室表面连续形变区域,并确定其位置以及形变量的目的.
3 实验处理与分析
3.1 洞室表面变形提取模拟实验
为了验证引入改进的Alpha Shapes算法识别洞室表面外轮廓点云和M3C2算法提取洞室表面微小形变量方法的正确性和有效性,设计模拟变形提取实验.选择洞室拱顶某处长宽均约0.5 m的点云数据作为实验区域,区域中包含一处锚杆(图7a).红色点云为第1期数据,蓝色点云为第2期数据,平均点距为5 mm.两期点云精配准均方根误差为2.57 mm.
在第2期数据Z坐标正向添加一个旋转抛物面的模拟变形,形成新的点云数据(下称“模拟数据”),如图7b所示.设变形中心的(x,y)坐标为(0,0),模拟变形范围(单位:m)为<0.1,模拟变形值为
(6)
分别对两期数据分别采用改进的Alpha Shapes算法识别洞室表面外轮廓点,α的取值为2倍平均点距,得到洞室外轮廓点云,如图8所示.
由图8可知,洞室表面外轮廓点云中,锚杆处点云簇和细小沟壑处点云已被剔除,剔除前后点云数量对比如表1所示.
由表1可知,第1期、第2期以及模拟数据中分别有19.11%、19.01%以及18.35%的噪声点被剔除.
表1Alpha Shapes算法前后点云数量变化
Table1Changes in number of points before and after the Alpha Shapes algorithm

将两期洞室表面外轮廓点云进行ICP配准,对配准后的第1期数据分别与第2期数据、模拟数据采用M3C2算法对每个点进行变形值计算,R值设定为3倍平均点距,r值为2倍平均点距,结果如图9
图7实验区数据
Fig.7Experimental area point cloud data

图8基于Alpha形状算法的外轮廓点云识别
Fig.8Recognition of outer contour point clouds via Alpha Shapes algorithm
图9变形值计算结果
Fig.9Deformation values
(7)
均方根误差为2.2 mm.实验显示,模拟变形值能够被较为完整、准确地反映出来.
3.2 实验概况
本文以某大型水力电站地下主副厂房洞作为研究对象设计了实验,主副厂房洞开挖长度为179.0 m,宽25.0 m,高57.0 m.实验数据收集采用Z+F IMAGER@5016三维激光扫描仪,获取了两期洞室点云数据.将易发生表面损伤的洞室拱顶部分作为实验区域,裁切点云获得洞室拱顶点云数据,如图10所示.

图10洞室拱顶点云数据
Fig.10Point cloud data of cavern vault
3.3 洞室点云预处理
实验数据采集中使用了4个相同的球形标靶,半径为72.5 mm,安置在洞室拱顶边缘.三维激光扫描仪架设在坐标原点O处,洞室中轴线方向与X轴同向,1、2号靶球安置在中轴线正向10 m处,3、4号靶球安置在中轴线负向10 m处,如图11所示.

图11三维激光扫描仪及靶球位置
Fig.11Position of 3D laser scanner and target balls
利用基于入射角加权法分别对两期球形标靶进行球心拟合,得到的球心坐标作为同名特征点,完成两期数据的粗配准,均方根误差为2.9 mm.得到的变换矩阵如表2所示,配准结果如图12所示.
表2粗配准参数
Table2Coarse registration parameters


图12两期点云粗配准
Fig.12Point clouds
a.before coarse registration; b.after coarse registration
取位于X轴正向一侧长度70 m的洞室表面点云数据,以0.4 m的固定厚度沿垂直于洞室X轴的方向切片,各切片x坐标区间为[0,0.4),···,[69.6,70.0],两期数据各产生175个点云切片.使用改进的Alpha Shapes轮廓点云识别算法处理各切片,α选取为2倍平均点距.实验区域两期点云数据各切片及其外轮廓点云的数据点个数如图13所示.
由图13可知,两期数据总体相似:x坐标距离扫描仪10 m以内的切片,点的剔除率约10%~20%;x坐标距离扫描仪10 m以上的切片,点的剔除率在10%以下.实验区域第1期数据总剔除率为6.25%,其中x坐标10 m以内剔除率为13.52%;第2期数据总剔除率为6.42%,其中x坐标10 m以内剔除率为14.17%.数据中的锥状凸起点云簇也已被剔除,如图14所示.
3.4 洞室表面变形监测
对得到的两期洞室外轮廓点云进行ICP精配准,对配准后的两期外轮廓点云采用M3C2算法计算变形值,R选取为3倍平均点距,r选取为2倍平均点距.计算得到21 648 255个核心点处的变形值,其平均值为0.5 mm,均方根误差为0.6 mm.核心点中变形值大于5 mm的形变点共44 103个,占所有核心点百分比为0.20%,其位置分布如图15所示.

图13各切片原点云及其外轮廓点云点数
Fig.13Point count of each slice origin point cloud and its outer contour point cloud

图14Alpha Shapes算法剔除点云簇
Fig.14Elimination of point cloud clusters via Alpha Shapes algorithm
利用同样的数据进行对比实验,不利用改进的Alpha Shapes算法处理洞室点云而直接提取形变点,形变点数量分布如表3所示.

图15识别出的形变点
Fig.15Identifying deformation points
表3形变点的变形值分布
Table3Distribution of deformation values for core points with 5 mm or more deformation values

由表3可知,未经改进的Alpha Shapes算法处理的点云,提取的形变点中有29.46%的点变形值在10 mm以上.而经改进的Alpha Shapes算法处理后,形变点中仅有8.26%的点变形值在10 mm以上,且高变形值区间内的形变点数量明显减少.两次实验中对同一区域提取的形变点如图16所示.
由图16b可知,未经Alpha Shapes算法处理,锚杆处锥状点云簇均被误判为形变点,而在图16a中,该部位噪声点已经过Alpha Shapes算法剔除,未造成误判.结果表明,识别出“混合像元”点云,可减少形变提取中噪声点被误提取为形变点.
将形变点进行距离聚类,共得到包含36 639个形变点在内的210个形变区域.目视判读形变区域,其中35个位于锚杆、钢筋、水泥柱周边,由于噪声点欠剔除而被误判为存在形变.其余175个形变区域平均变形值的频率分布如图17所示.
由图17可知,175个形变区域中,各区域平均变形值最大为18.4 mm,有147个形变区域的平均变形值在7.5 mm以下.形变区域的位置分布如图18所示.

图16形变点的比较
Fig.16Comparison of local deformation points identified from data

图17各形变区域平均变形值频率分布直方图
Fig.17Frequency distribution histogram of average deformation value in each deformation region
由图18可知,聚类得到的175个形变区域中,有140个区域分布于x坐标值在0~6 m范围内的形变区A中,可判断此处发生了大面积连续形变,其平均变形值为6.5 mm.其余35个形变区域离散分布于洞室拱顶各处,例如形变区B,定位其质心坐标为(34.95,4.29,17.87)(m),区域平均变形值为11.7 mm.验证了本文提出的方法剔除洞室表面点云噪声点,并定位和量化局部微小形变区域的可行性.
3.5 结果分析
改进的Alpha Shapes轮廓点云识别算法在点云密度较大、与洞室表面入射角近直角的区域,识别洞室表面轮廓点云的表现最好,噪声点剔除率达到15%,有效解决了由于洞室点云表面起伏引起的难以找到正确匹配点的问题.在锚杆、水泥柱和由于混合像元造成的噪声点的剔除方面,Alpha Shapes算法结合半径滤波,呈现了良好的噪声剔除功能.然而,在面对锚杆点云时,仍存在过剔除或欠剔除的现象,此部分的点云易在后续对比中被识别为形变,从而造成了16.7%的形变区域误判.

图18形变区域聚类
Fig.18Clustering of deformation areas
M3C2算法的洞室表面变形监测能够反映两期点云数据之间在各核心点法线方向上的变形,在经聚类提取形变区域的过程中,存在16.9%的零星粗差点未被聚类.
4 结论
本文针对三维激光扫描技术监测与提取洞室表面变形信息的问题,提出一种引入改进的Alpha Shapes算法的基于M3C2的洞室表面变形监测方法,通过实验验证得到以下结论:
1)针对洞室表面点云存在锚杆和水泥柱等突起、混合像元噪声以及表面起伏难以找到匹配点的问题,本文将改进的Alpha Shapes算法用于洞室表面外轮廓点云的识别,可有效排除以上噪声点云对变形值计算的影响.
2)本文设计的M3C2算法提取核心点变形值,通过距离聚类提取变形区域的方法,能够有效估计变形值并提取形变区域.