CN115690184A - 一种基于三维激光扫描的隧道掌子面位移量测方法 - Google Patents

Abstract

本发明公开了一种基于三维激光扫描的隧道掌子面位移量测方法，包括将全站式扫描仪布置于新开挖掌子面前，获取不同时期掌子面点云数据；对获取的掌子面点云数据进行处理，分离掌子面点云数据，根据基准点将不同时期掌子面点云数据转换到同一个坐标系中；设置体素块，将不同时期分离的掌子面点云数据分别分割为方形点云数据；采用点云特征描述子对方形点云数据进行配准，得到不同时期点云坐标对应关系；根据点云坐标对应关系，计算掌子面变形位移。本发明主要是实现了掌子面整***移变形测量，无需设置靶点，可实现掌子面上任意多个点的测量，并且避免了工作人员在施工初期与掌子面的直接接触，极大地保障了相关工作人员的人身安全。

Description

一种基于三维激光扫描的隧道掌子面位移量测方法

技术领域

本发明涉及掌子面位移监测技术领域，具体涉及一种基于三维激光扫描的隧道掌子面位移量测方法。

背景技术

在新建隧道开挖过程中，隧道掌子面作为直接工作区域，其稳定性应作为重点关注对象，掌子面附近往往集中了大量人员与机械设备，一旦掌子面失稳，将可能带来人员伤亡、机械设备受损、工期滞后等严重后果，因此有必要对掌子面变形进行实时监测。隧道掌子面变形能够直接反应隧道掌子面稳定状态，因此，实时获取掌子面变形具有重要意义。

目前，由于对掌子面的位移监测难度较大，实际较少进行该监测。对掌子面位移的监测，可采用在掌子面上设置监测靶标，用全站仪进行监测，这种方法虽然监测精度较高，但需在被测处放置特定靶标，新开挖掌子面是否稳定具有不确定性，安装测点靶标时不能确保人员安全；同时，相关测量都为点测量，无法获取掌子面整体变形情况。

三维激光扫描技术在一定条件下可用于作隧道变形非接触式测量，相较于一般的隧道变形量测，可实现对任意表面位置的位移测量，具有一定的技术优势。三维激光扫描可以获得整个观测物的表面空间信息，通过对点云数据的分析可以得到结构的变形情况，其变形分析将传统的固定单点分析拓展为整体分析，可以获取掌子面整体变形。

由于三维激光扫描的诸多优点，相关工作人员围绕三维激光扫描监测开展了大量的工作，其中，现有文献1公开了一种基于移动三维激光扫描的隧道掌子面变形报警方法，实现了基于特定监测点的变形监测；现有文献2公开了一种隧道局部变形识别方法及装置，通过统一大地坐标获取监测结构的相对变形；现有文献3公开了一种圆形隧道收敛变形的测量方法，将不同时期圆形隧道点云归于统一坐标系中，通过拟合算法分析获取圆形隧道收敛变形数据；

上述方法或针对掌子面进行特定点测量，或用于分析结构变形，而目前对掌子面任意点测量方法还较为缺乏。

发明内容

针对现有技术中的上述不足，本发明提供了一种基于三维激光扫描的隧道掌子面位移量测方法，实现准确、实时、全面的检测掌子面围岩动态变形的自动化测量。

为了达到上述发明目的，本发明采用的技术方案为：

一种基于三维激光扫描的隧道掌子面位移量测方法，包括以下步骤：

将全站式扫描仪布置于新开挖掌子面前，获取不同时期掌子面点云数据；

对获取的掌子面点云数据进行处理，分离掌子面点云数据，根据基准点将不同时期掌子面点云数据转换到同一个坐标系中；

设置体素块，将不同时期分离的掌子面点云数据分别分割为方形点云数据；

采用点云特征描述子对方形点云数据进行配准，得到不同时期点云坐标对应关系；

根据点云坐标对应关系，计算掌子面变形位移。

可选地，所述对获取的掌子面点云数据进行处理，具体包括：

根据隧道设计的标准中轴线数据和全站式扫描仪所处的里程和坐标数据，计算掌子面处的里程，然后对所有点云数据进行分割，得到包含掌子面的点云数据；

采用统计去噪及半径去噪的方式对包含掌子面的点云数据进行降噪处理；

采用基于距体素中心距离的降采样方法对降噪后的点云数据进行降采样。

可选地，所述采用基于距体素中心距离的降采样方法对包含掌子面的点云数据进行降采样，具体包括：

在三维空间中构建若干个设定边长的立方体栅格，将包含掌子面的点云数据通过栅格划分进行体素化；

从所有划分的体素栅格中剔除空的体素栅格；

以体素栅格为单位，对体素栅格内的点云数据进行降采样。

可选地，所述分离掌子面点云数据，具体包括：

从掌子面点云数据中任意选取一个点，通过宽度优先搜索所有与该点密度可达的点；如果该点是核心对象，则将该点和搜索的点记为同一个集合；如果该点是边界对象，则将搜索的点记为噪点，直到找到一个完整的族；

随机从掌子面点云数据中重新选取一个点，重复上述操作，得到下一个族，直到所有的点都被标记过为止。

可选地，所述利用点云特征描述子对分块的掌子面点云数据进行配准，得到不同时期点云坐标对应关系，具体包括：

根据邻域内的任意两点p_i、p_j及其对应的法线n_i、n_j，将两点构成的直线与法线的夹角较小的点作为源点p_s，另一点为p_t，其法线分别为n_s、n_t，并以源点建立坐标轴；

根据建立的坐标系计算以下参数：

根据上式参数计算两点的

特征，其中α＝v·n_t、

θ＝arctan(w·u，u·n_t)，并以直方图的形式进行统计，形成简化的点直方图SPFH，作为点云特征描述子；

利用点云特征描述子确定源点云与目标点云中的一致性匹配点对，对点云数据进行全局配准；

采用Point-to-Plane ICP算法对点云数据进行精配准；

根据全局配准的目标变换矩阵和精配准的旋转变换矩阵得到最终变换矩阵，根据最终变换矩阵计算源点云与目标点云的最佳配准状态，得到源点云与目标点云的对应关系。

可选地，所述利用点云特征描述子确定源点云与目标点云中的一致性匹配点对，对点云数据进行全局配准，具体包括：

利用点云特征描述子确定源点云P＝{p_i}，i＝0，1，2…，n与目标点云Q＝{q_i}，i＝0，1，2…，m中的一致性匹配点对K＝{(p_k，q_k)}；

以匹配点对间的空间距离最小作为目标函数，构建目标变换矩阵的估计函数；

其中，E(T)为目标变化矩阵估计函数；p(x)为稳健优化器；μ为优化器参数；x为控制残差参数；Tq为变换后的目标对应点；

采用Black-Rangarajan对偶和线性分解方法对目标变换矩阵的估计函数进行优化求解，得到全局配准的目标变换矩阵。

可选地，所述采用Black-Rangarajan对偶和线性分解方法对目标变换矩阵的估计函数进行优化求解，得到全局配准的目标变换矩阵，具体包括：

根据一致性点对间的线形关系，对目标变换矩阵的估计函数进行优化，表示为

其中，E(T，L)为目标变换矩阵T的估计函数，T为求解的目标变换矩阵，L＝{l_p，q}为一致性点对间的线形关系，l_p，q为原点云与目标点云对应点之间的距离，(p，q)为原点云与目标点云对应点，p为原点云，q为目标点云，K为对应点云集合，ψ(l_p，q)为限制函数，μ为；

将目标变换矩阵T中的旋转分量w和平移分量t线形分解为6个向量，表示为

ξ＝(w，t)＝(α，β，γ，t_x，t_y，t_z)

其中，α为沿X轴旋转分量，β为沿Y轴旋转分量，γ为沿Z轴旋转分量，t_x为X方向位移分量，t_y为Y方向位移分量，t_z为Z方向位移分量；

则变换矩阵T可由ξ通过下式表示：

其中，T^k为最后一次迭代的估计变换矩阵；

采用高斯-牛顿法求解分解向量，表示为：

其中，

为分解向量转置，J_r为其对应的雅各比矩阵，r为剩余向量；

根据分解向量ξ求解目标变换矩阵T，再使用求解得到目标变换矩阵的估计函数E(T，L)，即可获取粗配准变换矩阵。

可选地，所述采用Point-to-Plane ICP算法对点云数据进行精配准具体包括：

以源点云与目标点云经过全局配准的对应点对所在切平面的距离的平方和最小化作为目标，建立目标损失函数为

M＝T(t_x，t_y，t_z)·R(α，β，γ)

r₁₁＝cosγcosβ

r₁₂＝-sinγcosα+cosγsinβsinα，

r₁₃＝sinγsinα+cosγsinβcosα

r₂₁＝sinγcosβ

r₂₂＝cosγcosα+sinγsinβsinα

r₂₃＝-cosγsinα+sinγsinβcosα，

r₃₁＝-sinβ

r₃₂＝cosβsinα

r₃₃＝cosβcosα

其中，M_opt为目标损失函数，M为源点云和目标点云之间的刚性变换矩阵，i为一致对应点序号，T(t_x，t_y，t_z)为矩阵变换平移分量，R(α，β，γ)为矩阵变化旋转分量，R_x(α)、R_y(β)、R_z(γ)为旋转矩阵；

对旋转变换矩阵进行近似，表示为

其中，

为近似的旋转变换矩阵；

根据近似的旋转变换矩阵对目标损失函数进行转换，表示为

其中，

为转换的目标损失函数，

为变换后的刚性变换矩阵，p_i为一致对应点源点，q_i为一致对应点目标点，n_i为目标点切平面法向量；

将转换后的目标损失函数进行展开，表示为：

其中，p_ix为对应点原点x坐标值，p_iy为对应点原点y坐标值，p_iz为对应点原点z坐标值，q_ix为对应点目标点x坐标值，q_iy为对应点目标点y坐标值，q_iz为对应点目标点z坐标值，n_ix为目标点切平面法向量的x分量，n_iy为目标点切平面法向量的y分量，n_iz为目标点切平面法向量的z分量，α，β，γ为旋转分量，t_x、t_y、t_z为平移分量；

将所有点对的目标损失函数采用矩阵表示为

Ax-b

a_i1＝n_izp_iy-n_iyp_iz

a_i2＝n_ixp_iz-n_izp_ix

a_i3＝n_iyp_ix-n_ixp_iy

x＝(α β γ t_x t_y t_z)^T

其中，n_Nx为第N个对应目标点切平面法向量的x分量对应点，n_Ny为第N个对应目标点切平面法向量的y分量对应点，n_Nz为第N个对应目标点切平面法向量的z分量对应点，p_Nx为第N个对应源点x分量，p_Ny为第N个对应源点y分量，p_Nz为第N个对应源点z分量，q_Nx为第N个对应目标点x分量，q_Ny为第N个对应目标点y分量，q_Nz为第N个对应目标点z分量；

根据目标损失函数的矩阵表示，将目标损失函数转换为

x_opt＝argmin_x|Ax-b|²

其中，x_opt为转换的目标损失函数；

对转换的目标损失函数采用SVD分解进行求解，得到旋转变换矩阵。

可选地，还包括对掌子面变形进行可视化处理以及信息整理。

可选地，所述对掌子面变形进行可视化处理以及信息整理具体包括：

通过对变形值进行等比划分，以结构设计允许变形值为限值，对不同变形值点云进行染色，变形小于设计变形允许值，颜色由蓝到红，超过设计允许值使用灰色作为标注，确定变形值与颜色之间对应关系，使得变形直观显示；

统计大于某一变形值的变形区域数量，并与变形限制值进行对比，显示掌子面安全状态信息。

本发明具有以下有益效果：

本发明仅通过掌子面固有特征，可以实现对掌子面变形的监控量测，能够显著提高检测范围，实现掌子面的全面测量；并且本发明以掌子面的空间拓扑作为配准目标，不需要在隧道主体结构上选择或埋设标志物作为测点，且避免了工作人员在施工初期与掌子面的直接接触，极大地保障了相关工作人员的人身安全。

附图说明

图1为本发明实施例中一种基于三维激光扫描的隧道掌子面位移量测方法的流程示意图。

图2为本发明实施例中定位棱镜的安装示意图；

图3为本发明实施例中点云配准示意图；

图4为本发明实施例中点云坐标轴构建示意图；

图5为本发明实施例中Point-to-Plane ICP算法示意图。

具体实施方式

下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

如图1所示，本发明实施例提供了一种基于三维激光扫描的隧道掌子面位移量测方法，包括以下步骤S1至S5：

S1、将全站式扫描仪布置于新开挖掌子面前，根据基准点获取不同时期掌子面点云数据；

在本发明的一个可选实施例中，本发明将三套棱镜分别安装于隧道二次衬砌之上，利用全站式扫描仪中全站功能获取其相对坐标，建立基准点；

其中，定位棱镜布置于二次衬砌之上，应当尽可能靠近掌子面，错位布置，装置为可拆卸结构，以其中两点作为基准点，另一点作为复核点。

本发明不需要在掌子面上设置监测点，且全站式扫描仪安装位置距掌子面有一定距离，提高了监测技术人员的安全。

本发明将全站式扫描仪布置于新开挖掌子面前，并对基准点进行扫描，获取全站式扫描仪定位信息，然后使用扫描功能对掌子面进行扫描，获取不同时期掌子面点云数据；

其中，全站式扫描仪安装位置距离掌子面控制在20-40m之间。在进行掌子面扫描前应当先对已知点进行全站扫描，获取全站式扫描仪基础状态；隧道开挖及找顶完成后，应尽量提早对隧道掌子面进行扫描，并记录其扫描时间，以初次扫描点云作为源点云，作为隧道掌子面变形对比分析的基准。

随后，结合现场情况以及监测要求，重复上述步骤，获取多期点云数据。

本发明对掌子面扫描数据采集时间短，扫描数据量大，提高了掌子面数据采集效率。

S2、对获取的掌子面点云数据进行处理，分离掌子面点云数据，根据基准点将不同时期掌子面点云数据转换到同一个坐标系中；

在本发明的一个可选实施例中，本发明对获取的掌子面点云数据进行处理，具体包括：

根据隧道设计的标准中轴线数据和全站式扫描仪所处的里程和坐标数据，计算掌子面处的里程，然后对所有点云数据进行分割，得到包含掌子面的点云数据；

采用统计去噪及半径去噪的方式对包含掌子面的点云数据进行降噪处理；

采用基于距体素中心距离的降采样方法对降噪后的点云数据进行降采样。

具体而言，本发明对点云数据进行处理使用包括软件Cloud-compare，python中open3D扩展包等程序；

首先，依据隧道设计的标准中轴线数据和全站式扫描仪所处的里程、坐标，计算出掌子面处的里程，使用Cloud-compare对点云进行分割，获取包含掌子面的点云数据

其次，用统计去噪及半径去噪的方式对点云进行降噪处理；

然后，采用基于距体素中心距离的降采样方法对点云进行降采样。

本发明采用的统计去噪是通过统计邻域内的点云距离，计算其标准差，利用高斯分布来剔除邻域内平均距离大于阈值的点，具体方法为：

对于点云中任一点p到其邻域内的点p_i(i＝1，2，...，n)的平均距离，设点p的坐标为{x，y，z}，点p_i的坐标为{x_i，y_i，z_i}，点p到点p_i的距离为d_i，其平均距离为

则

由于点云扫描的离散型，可假定点云中的点与其临域内的点的平均距离

在统计意义上满足正态分布，即N(μ，σ²)的分布，则点云距离的数学期望μ与标准差σ的计算公式如下：

其中，

为点云中第i个点到其临域内的点的平均距离；m为点云中的点数。

考虑正态分布的原理设置距离阈值ε，并将点云内到邻域内最近临近点距离大于ε的点噪声剔除：

ε＝θσ+μ

其中，θ为比例系数，与临域内的点数有关。

本发明采用的半径去噪利用距离阈值将点云数据分簇，将簇中点数低于数量阈值的点云作为离群点去除，具体方法为：

对于点云中任一点p，查找其邻域内欧式距离小于距离阈值的所有点，标记为近邻n_i(n＝1，2，...，n)，对于n_i仍按照距离阈值获取近邻，直至没有新的近邻出现，则将上述点划分为一簇。通过数量阈值将簇中点数低于阈值的簇作为离群点去除，完成去噪。

本发明采用基于距体素中心距离的降采样方法对包含掌子面的点云数据进行降采样，具体包括：

在三维空间中构建若干个设定边长的立方体栅格，将包含掌子面的点云数据通过栅格划分进行体素化；

从所有划分的体素栅格中剔除空的体素栅格；

以体素栅格为单位，对体素栅格内的点云数据进行降采样。

在本发明的一个可选实施例中，本发明分离掌子面点云数据，具体包括：

从掌子面点云数据中任意选取一个点，通过宽度优先搜索所有与该点密度可达的点；如果该点是核心对象，则将该点和搜索的点记为同一个集合；如果该点是边界对象，则将搜索的点记为噪点，直到找到一个完整的族；

随机从掌子面点云数据中重新选取一个点，重复上述操作，得到下一个族，直到所有的点都被标记过为止。

具体而言，掌子面与隧道中轴线形成的夹角较大，而掌子面周边区域曲面近似与隧道中轴线平行，夹角较小，将点云投影到与隧道中轴线垂直的平面。掌子面投影区域大量重合，导致隧道轮廓处的点云非常密集，轮廓内部掌子面区域的点云较为稀疏，呈现明显的特征；因此本发明使用基于密度聚类算法对点云进行分离，具体方法为：

步骤一：从任意个体i开始，通过宽度优先搜索所有与i密度可达的个体，如果个体i是核心对象，将它们记为同一个集合；如果个体i是边界对象，就将之记为噪点，直到找到一个完整的族：

步骤二：随机选取一个新的个体i，进行处理，得到下一个族。算法一直进行下去，直到所有的个体都被标记过为止。

直接密度可达：如果i是一个核心对象，且j在i的邻域内，即点j在以某i为圆心、R为半径的圆范围，则i到j直接密度可达。

密度可达：给定个体P₁，P₂，…，P_n，任意个体P_i，，其中i＝l,2,…,n，存在P_i到P_i+1直接密度可达，则称P₁到P_n密度可达。

最后，本发明根据建站信息将两期点云通过坐标统一的方式转换至同一坐标系下。

S3、设置体素块，将不同时期分离的掌子面点云数据分别分割为方形点云数据；

在本发明的一个可选实施例中，本发明设置体素块大小为0.15*0.15*1m的矩形块，将不同时期分离的掌子面点云数据分割为0.15*0.15的方形点云数据。

S4、采用点云特征描述子对方形点云数据进行配准，得到不同时期点云坐标对应关系；

在本发明的一个可选实施例中，本发明利用点云特征描述子对分块的掌子面点云数据进行配准，得到不同时期点云坐标对应关系，具体包括：

根据邻域内的任意两点及其对应的法线，将两点构成的直线与法线的夹角较小的点作为原点，从原点建立坐标轴；

计算两点的

特征，并以直方图的形式进行统计，形成简化的点直方图SPFH，作为点云特征描述子；

利用点云特征描述子确定源点云与目标点云中的一致性匹配点对，对点云数据进行全局配准；

采用Point-to-Plane ICP算法对点云数据进行精配准；

根据全局配准的目标变换矩阵和精配准的旋转变换矩阵得到最终变换矩阵，根据最终变换矩阵计算源点云与目标点云的最佳配准状态，得到源点云与目标点云的对应关系。

其中利用点云特征描述子确定源点云与目标点云中的一致性匹配点对，对点云数据进行全局配准，具体包括：

利用点云特征描述子确定源点云P＝{p_i}，i＝0，1，2…，n与目标点云Q＝{q_i}，i＝0，1，2…，m中的一致性匹配点对K＝{(p_k，q_k)}；

以匹配点对间的空间距离最小作为目标函数，构建目标变换矩阵的估计函数；

其中，(T)为目标变化矩阵估计函数；ρ(x)为稳健优化器；μ为优化器参数；x为控制残差参数；Tq为变换后的目标对应点；

采用Black-Rangarajan对偶和线性分解方法对目标变换矩阵的估计函数进行优化求解，得到全局配准的目标变换矩阵，具体包括：

根据一致性点对间的线形关系，对目标变换矩阵的估计函数进行优化，表示为

其中，E(T，L)为目标变换T的估计函数，T为求解目标变换矩阵，L＝{l_p，q}为一致性点对间的线形关系，l_p，q为原点云与目标点云对应点之间的距离，(p，q)为原点云与目标点云对应点，p为原点云，q为目标点云，K为对应点云集合，ψ(l_p，q)为限制函数；

将目标变换矩阵T中的旋转分量w和平移分量t线形分解为6个向量，表示为：

ξ＝(w，t)＝(α，β，γ，t_x，t_y，t_z)

其中，α为沿x轴旋转分量，β为沿y轴旋转分量，γ为沿z轴旋转分量，t_x为x方向位移分量，t_y为y方向位移分量，t_z为z方向位移分量；

则变换矩阵T可由ξ通讨下式表示：

其中，T^k为最后一次迭代的估计变换矩阵；

采用高斯-牛顿法求解分解向量，表示为

其中，

为分解向量转置，J_r为其对应的雅各比矩阵(Jacobian)，r为剩余向量；

根据分解向量求解目标变换矩阵，再求解得到目标变换矩阵的估计函数。

其中采用Point-to-Plane ICP算法对点云数据进行精配准具体包括：

以源点云与目标点云经过全局配准的对应点对所在切平面的距离的平方和最小化作为目标，建立目标损失函数为

M＝T(t_x，t_y，t_z)·R(α，β，γ)

r₁₁＝cosγcosβ

r₁₂＝-sinγcosα+cosγsinβsinα，

r₁₃＝sinγsinα+cosγsinβcosα

r₂₁＝sinγcosβ

r₂₂＝cosγcosα+sinγsinβsinα

r₂₃＝-cosγsinα+sinγsinβcosα，

r₃₁＝-sinβ

r₃₂＝cosβsinα

r₃₃＝cosβcosα

其中，M_opt为目标损失函数，M为源点云和目标点云之间的刚性变换矩阵，i为一致对应点序号，T(t_x，t_y，t_z)为矩阵变换平移分量，R(α，β，γ)为矩阵变化旋转分量，R_x(α)、R_y(β)、R_z(γ)为旋转矩阵；

对旋转变换矩阵进行近似，表示为

其中，

为近似的旋转变换矩阵；

根据近似的旋转变换矩阵对目标损失函数进行转换，表示为

其中，

为转换的目标损失函数，

变换后的刚性变换矩阵，p_i为一致对用点源点，q_i为一致对应点目标点，n_i为目标点切平面法向量；

将转换后的目标损失函数进行展开，表示为：

其中，p_ix为对应点原点x坐标值，p_iy为对应点原点y坐标值，p_iz为对应点原点z坐标值，q_ix为对应点目标点x坐标值，q_iy为对应点目标点y坐标值，q_iz为对应点目标点z坐标值，n_ix为目标点切平面法向量的x分量，n_iy为目标点切平面法向量的y分量，n_iz为目标点切平面法向量的z分量，α，β，γ为旋转分量，t_x、t_y、t_z为平移分量；

假设共有N对对应点对(p_i，q_i)，1≤i≤N；通过上述的线性近似处理，将所有点对的目标损失函数采用矩阵表示为

Ax-b

a_i1＝n_izp_iy-n_iyp_iz

a_i2＝n_ixp_iz-n_izp_ix

a_i3＝n_iyp_ix-n_ixp_iy

x＝(α β γ t_x t_y t_z)^T

其中，n_Nx为第N个对应目标点切平面法向量的x分量对应点，n_Ny为第N个对应目标点切平面法向量的y分量对应点，n_Nz为第N个对应目标点切平面法向量的z分量对应点，p_Nx为第N个对应源点x分量，p_Ny为第N个对应源点y分量，p_Nz为第N个对应源点z分量，q_Nx为第N个对应目标点x分量，q_Ny为第N个对应目标点y分量，q_Nz为第N个对应目标点z分量；

根据目标损失函数的矩阵表示，将目标损失函数转换为

x_opt＝argmin_x|Ax-b|²

其中，x_opt为转换的目标损失函数；

对转换的目标损失函数采用SVD分解进行求解，得到旋转变换矩阵。

具体而言，如图3所示，点云配准过程即寻求一个包含平移和旋转的刚性变换(Rigid transformation)，使得不同时期的点云数据的空间位置差异最小化。

首先，特征匹配为表现为识别不同局部点云之间特征，求得最相似区域，以其中两片点云质心作为对应同名点，本发明使用特征描述为FPFH(Fast Point FeatureHistogram)描述子，所述点云描述字表现为计算该点及与其直接相连的近邻点的(α，

θ)3个特征，将上述特征以直方图的形式进行统计，形成简化的点直方图SPFH(Simple PointFeature Histograms)。

特征描述子计算如下所示：

对于邻域内任两点p_i、p_j，及其对应的法线n_i、n_j，选定两点p_s与p_t，将上述两点构成的直线与法线的夹角较小的点作为原点，则从原点建立坐标轴如图4所示。

其次，已知源点云S＝{A_i}，i＝0，1，2…，n和目标点云T＝{A′_i}，i＝0，1，2…，m；存在一个初始刚性变换矩阵M，使得S，T在三维空间中对齐。通过空间点云FPFH描述子初步确定源点云与目标点云中的一致性匹配点对(Correspondence sets)K＝{(p_k，q_k)}。则所求目标变换T能使得配准后的匹配点对间的空间距离最小。定义针对目标变换矩阵的估计函数：

其中，E(T)为目标变化矩阵估计函数；ρ(x)为稳健性优化器(Robust penalty)，用于提高计算效率，避免在优化目标变换的过程中采样、验证、删减或重新计算对应关系。

由于上式难以直接进行优化求解，故算法在稳健估计和线性过程之间使用了Black-Rangarajan对偶，具体如下：

设L＝{l_p，q}为一致性点对间的线形关系，则对目标变换T的估计函数进行优化如下：

E(T，L)最小时，函数偏导满足下式：

求解得：

将上式代入目标变换矩阵的估计函数即可求解目标变换矩阵T。

考虑计算复杂度，使用上述计算步骤的基础上再此对算法进行优化：

1)将目标变换矩阵T中的旋转分量w和平移分量t线形分解为6个向量：

ξ＝(w，t)＝(α，β，γ，a，b，c)

则变换矩阵T可由ξ通过下式表示：

上式中的T^k为最后一次迭代的估计变换矩阵；

2)使用高斯-牛顿法求解ξ：

上式中，r为剩余向量，J_r为其对应的雅各比矩阵(Jacobian)；

3)由上式，将ξ代入T^k中，求解T，再将T回代入式(8)计算E(T，L)；

然后在全局配准的基础上使用Point-to-Plane ICP算法对点云进行精配准，具体方法为：

已知源点云P和目标点云Q，将经过初始变换的源点云中的点与目标点云建立对应关系，设对应点对为(p_i，q_i)。Point-to-Plane ICP以p_i到q_i所在切平面的距离的平方和作为评价标准，如图5所示，Point-to-Plane ICP的算法目标即是使此距离最小化。

首先定义目标损失函数：

其中，p_i表示源点云上的点，q_i表示目标点云上p_i的对应点，n_i为q_i所在切平面的法向量，M为源点云和目标点云之间的刚性变换。M包含了点云在空间中的平移和旋转，以矩阵形式进行表示：

M＝T(t_x，t_y，t_z)·R(α，β，γ)

r₁₁＝cosγcosβ

r₁₂＝-sinγcosα+cosγsinβsinα，

r₁₃＝sinγsinα+cosγsinβcosα

r₂₁＝sinγcosβ

r₂₂＝cosγcosα+sinγsinβsinα

r₂₃＝-cosγsinα+sinγsinβcosα，

r₃₁＝-sinβ

r₃₂＝cosβsinα

r₃₃＝cosβcosα

为降低求解复杂度，将非线性求解问题进行线性近似，方法如下：

当角度θ≈0时，有sinθ≈θ，cosθ≈1，因此当θ，β，γ≈0时，可对旋转变换进行如下近似：

则所求最优变换可近似为：

将上式代入目标损失函数，得：

将上式展开即完成了求解非线性最小二乘问题的线性近似：

假设共有N对对应点对(p_i，q_i)，1≤i≤N；通过上述的线性近似处理，将所有点对的目标损失函数采用矩阵表示为

Ax-b

其中：

a_i1＝n_izp_iy-n_iyp_iz

a_i2＝n_ixp_iz-n_izp_ix

a_i3＝n_iyp_ix-n_ixp_iy

x＝(α β γ t_x t_y t_z)^T

由此，对于目标损失函数的求解即可转化为如下形式：

x_opt＝argmin_x|Ax-b|²

上式通过SVD分解进行求解：

x_opt＝A*b

最后，根据全局配准模型坐标变换矩阵以及精配准获得最终变换矩阵，通过获取两期点云的最佳配准状态，即可获得两片点云对应关系。

S6、根据点云坐标对应关系，计算掌子面变形位移。

在本发明的一个可选实施例中，本发明根据点云对应关系，并计算两片点云之间欧氏距离，该距离即为掌子面在该范围内变形值；根据配准结果可知不同时期点云间对应关系，即：点云A(x_i，y_i，z_i)和点云A′(x_j，y_j，z_j)之间欧氏距离平均值计算有

该距离就是掌子面在该区域的变形值。

在本发明的一个可选实施例中，本发明还包括对掌子面变形进行可视化处理以及信息整理，具体包括：

通过对变形值进行等比划分，以结构设计允许变形值为限值，对不同变形值点云进行染色，变形小于设计变形允许值，颜色由蓝到红，超过设计允许值使用灰色作为标注，确定变形值与颜色之间对应关系，使得变形直观显示；

统计大于某一变形值的变形区域数量，并与变形限制值进行对比，显示掌子面安全状态信息。

本发明可对较多的掌子面特征点进行位移分析，实现了对掌子面的整***移分析，增强了掌子面位移分析效果。

综上所述，通过本发明方案，可以实现对掌子面变形的监控量测，在保障相关工作人员的安全的同时，测量精度也满足规范要求。本发明以掌子面的空间拓扑作为配准目标，不需要在隧道主体结构上选择或埋设标志物作为测点，且避免了工作人员在施工初期与掌子面的直接接触，极大地保障了相关工作人员的人身安全。

本发明是参照根据本发明实施例的方法、设备(***)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

本发明中应用了具体实施例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (10)

1.一种基于三维激光扫描的隧道掌子面位移量测方法，其特征在于，包括以下步骤：

将全站式扫描仪布置于新开挖掌子面前，获取不同时期掌子面点云数据；

对获取的掌子面点云数据进行处理，分离掌子面点云数据，根据基准点将不同时期掌子面点云数据转换到同一个坐标系中；

设置体素块，将不同时期分离的掌子面点云数据分别分割为方形点云数据；

采用点云特征描述子对方形点云数据进行配准，得到不同时期点云坐标对应关系；

根据点云坐标对应关系，计算掌子面变形位移。

2.根据权利要求1所述的一种基于三维激光扫描的隧道掌子面位移量测方法，其特征在于，所述对获取的掌子面点云数据进行处理，具体包括：

根据隧道设计的标准中轴线数据和全站式扫描仪所处的里程和坐标数据，计算掌子面处的里程，然后对所有点云数据进行分割，得到包含掌子面的点云数据；

采用统计去噪及半径去噪的方式对包含掌子面的点云数据进行降噪处理；

采用基于距体素中心距离的降采样方法对降噪后的点云数据进行降采样。

3.根据权利要求2所述的一种基于三维激光扫描的隧道掌子面位移量测方法，其特征在于，所述采用基于距体素中心距离的降采样方法对包含掌子面的点云数据进行降采样，具体包括：

在三维空间中构建若干个设定边长的立方体栅格，将包含掌子面的点云数据通过栅格划分进行体素化；

从所有划分的体素栅格中剔除空的体素栅格；

以体素栅格为单位，对体素栅格内的点云数据进行降采样。

4.根据权利要求1所述的一种基于三维激光扫描的隧道掌子面位移量测方法，其特征在于，所述分离掌子面点云数据，具体包括：

从掌子面点云数据中任意选取一个点，通过宽度优先搜索所有与该点密度可达的点；如果该点是核心对象，则将该点和搜索的点记为同一个集合；如果该点是边界对象，则将搜索的点记为噪点，直到找到一个完整的族；

随机从掌子面点云数据中重新选取一个点，重复上述操作，得到下一个族，直到所有的点都被标记过为止。

5.根据权利要求1所述的一种基于三维激光扫描的隧道掌子面位移量测方法，其特征在于，所述利用点云特征描述子对分块的掌子面点云数据进行配准，得到不同时期点云坐标对应关系，具体包括：

根据邻域内的任意两点p_i、p_j及其对应的法线n_i、n_j，将两点构成的直线与法线的夹角较小的点作为源点p_s，另一点为p_t，其法线分别为n_s、n_t，并以源点建立坐标轴；

根据建立的坐标系计算以下参数：

根据上式参数计算两点的

特征，其中α＝v·n_t、

θ＝arctan(w·u，u·n_t)，并以直方图的形式进行统计，形成简化的点直方图SPFH，作为点云特征描述子；

利用点云特征描述子确定源点云与目标点云中的一致性匹配点对，对点云数据进行全局配准；

采用Point-to-Plane ICP算法对点云数据进行精配准；

根据全局配准的目标变换矩阵和精配准的旋转变换矩阵得到最终变换矩阵，根据最终变换矩阵计算源点云与目标点云的最佳配准状态，得到源点云与目标点云的对应关系。

6.根据权利要求5所述的一种基于三维激光扫描的隧道掌子面位移量测方法，其特征在于，所述利用点云特征描述子确定源点云与目标点云中的一致性匹配点对，对点云数据进行全局配准，具体包括：

利用点云特征描述子确定源点云P＝{p_i}，i＝0，1，2...，n与目标点云Q＝{q_i}，i＝0，1，2...，m中的一致性匹配点对K＝{(p_k，q_k)}；

以匹配点对间的空间距离最小作为目标函数，构建目标变换矩阵的估计函数；

其中，E(T)为目标变化矩阵估计函数；ρ(x)为稳健优化器；μ为优化器参数，x为控制残差参数；Tq为变换后的目标对应点；

采用Black-Rangarajan对偶和线性分解方法对目标变换矩阵的估计函数进行优化求解，得到全局配准的目标变换矩阵。

7.根据权利要求6所述的一种基于三维激光扫描的隧道掌子面位移量测方法，其特征在于，所述采用Black-Rangarajan对偶和线性分解方法对目标变换矩阵的估计函数进行优化求解，得到全局配准的目标变换矩阵，具体包括：

根据一致性点对间的线形关系，对目标变换矩阵的估计函数进行优化，表示为

其中，E(T，L)为目标变换矩阵T的估计函数，T为求解的目标变换矩阵，L＝{l_p，q}为一致性点对间的线形关系，l_p，q为原点云与目标点云对应点之间的距离，(p，q)为原点云与目标点云对应点，p为原点云，q为目标点云，K为对应点云集合，ψ(l_p，q)为限制函数；

将目标变换矩阵T中的旋转分量w和平移分量t线形分解为6个向量，表示为

ξ＝(w，t)＝(α，β，γ，t_x，t_y，t_z)

其中，α为沿X轴旋转分量，β为沿Y轴旋转分量，γ为沿Z轴旋转分量，t_x为X方向位移分量，t_y为Y方向位移分量，t_z为Z方向位移分量；

则变换矩阵T可由ξ通过下式表示：

其中，T^k为最后一次迭代的估计变换矩阵；

采用高斯-牛顿法求解分解向量，表示为：

其中，

为分解向量转置，J_r为其对应的雅各比矩阵，r为剩余向量；

根据分解向量ξ求解目标变换矩阵T，再使用求解得到目标变换矩阵的估计函数E(T，L)，即可获取粗配准变换矩阵。

8.根据权利要求5所述的一种基于三维激光扫描的隧道掌子面位移量测方法，其特征在于，所述采用Point-to-Plane ICP算法对点云数据进行精配准具体包括：

以源点云与目标点云经过全局配准的对应点对所在切平面的距离的平方和最小化作为目标，建立目标损失函数为

M＝T(t_x，t_y，t_z)·R(α，β，γ)

r₁₁＝cosγcosβ

r₁₂＝-sinγcosα+cosγsinβsinα，

r₁₃＝sinγsinα+cosγsinβcosα

r₂₁＝sinγcosβ

r₂₂＝cosγcosα+sinγsinβsinα

r₂₃＝-cosγsinα+sinγsinβcosα，

r₃₁＝-sinβ

r₃₂＝cosβsinα

r₃₃＝cosβcosα

其中，M_opt为目标损失函数，M为源点云和目标点云之间的刚性变换矩阵，i为一致对应点序号，T(t_x，t_y，t_z)为矩阵变换平移分量，R(α，β，γ)为矩阵变化旋转分量，R_x(α)、R_y(β)、R_z(γ)为旋转矩阵；

对旋转变换矩阵进行近似，表示为

其中，

为近似的旋转变换矩阵；

根据近似的旋转变换矩阵对目标损失函数进行转换，表示为

其中，

为转换的目标损失函数，

为变换后的刚性变换矩阵，p_i为一致对应点源点，q_i为一致对应点目标点，n_i为目标点切平面法向量；

将转换后的目标损失函数进行展开，表示为：

其中，p_ix为对应点原点x坐标值，p_iy为对应点原点y坐标值，p_iz为对应点原点z坐标值，q_ix为对应点目标点x坐标值，q_iy为对应点目标点y坐标值，q_iz为对应点目标点z坐标值，n_ix为目标点切平面法向量的x分量，n_iy为目标点切平面法向量的y分量，n_iz为目标点切平面法向量的z分量，α，β，γ为旋转分量，t_x、t_y、t_z为平移分量；

将所有点对的目标损失函数采用矩阵表示为

Ax-b

a_i1＝n_izp_iy-n_iyp_iz

a_i2＝n_ixp_iz-n_izp_ix

a_i3＝n_iyp_ix-n_ixp_iy

x＝(α β γ t_x t_y t_z)^T

其中，n_Nx为第N个对应目标点切平面法向量的x分量对应点，n_Ny为第N个对应目标点切平面法向量的y分量对应点，n_Nz为第N个对应目标点切平面法向量的z分量对应点，p_Nx为第N个对应源点x分量，p_Ny为第N个对应源点y分量，p_Nz为第N个对应源点z分量，q_Nx为第N个对应目标点x分量，q_Ny为第N个对应目标点y分量，q_Nz为第N个对应目标点z分量；

根据目标损失函数的矩阵表示，将目标损失函数转换为

x_opt＝argmin_x|Ax-b|²

其中，x_opt为转换的目标损失函数；

对转换的目标损失函数采用SVD分解进行求解，得到旋转变换矩阵。

9.根据权利要求1所述的一种基于三维激光扫描的隧道掌子面位移量测方法，其特征在于，还包括对掌子面变形进行可视化处理以及信息整理。

10.根据权利要求9所述的一种基于三维激光扫描的隧道掌子面位移量测方法，其特征在于，所述对掌子面变形进行可视化处理以及信息整理具体包括：

通过对变形值进行等比划分，以结构设计允许变形值为限值，对不同变形值点云进行染色，变形小于设计变形允许值，颜色由蓝到红，超过设计允许值使用灰色作为标注，确定变形值与颜色之间对应关系，使得变形直观显示；

统计大于某一变形值的变形区域数量，并与变形限制值进行对比，显示掌子面安全状态信息。

Images