传统的点云配准方法^[1]是基于离散点来实施的, 虽然取得了一定的成效, 但其配准精度易受点特征提取精度的影响^[2], 由于被扫描的物体之间相互遮挡, 寻找完全对应的同名点特征十分困难。相对而言, 直线特征则具有丰富和稳定的几何特性, 广泛存在于被扫描的对象上, 即使存在遮挡情况, 也只需提取点特征所在的直线特征, 既能降低实施的复杂度, 又能保证相邻测站点云之间的高精度配准。近年来, 有学者对线特征点云配准方法进行了研究, 提出基于四元数表示的点云配准直接解法^[3-4], 该方法虽无需迭代, 却不能同时求解所有的配准参数; 王永波等^[5]提出一种基于Plücker直线描述的点云配准非迭代方法, 提高了配准精度; 盛庆红等^[6-7]提出Plücker直线表达的配准模型, 充分利用了直线的几何约束能力, 使配准模型形式更加简洁。本文使用直线特征代替传统的点特征配准基元, 根据点在配准平面上的几何关系, 采用单位四元数表征点云配准模型, 从而实现配准参数的同步求解, 增强了理论的严密性。该方法不仅适用于遮挡造成的点云配准精度低的问题, 也适用于大旋转角的点云配准。

1 单位四元数表示的旋转运动

四元数的基本定义为^[8]：

$ \mathit{\boldsymbol{\widehat h}} = {h_0} + {h_x}{\rm{i}} + {h_y}{\rm{j}} + {h_z}{\rm{k}} $

(1)

式中, i、j、k为虚数单位, h₀和h_x、h_y、h_z分别为实部和虚部。

若四元数$\mathit{\boldsymbol{\hat h}}$的模‖$\mathit{\boldsymbol{\hat h}}$‖²=1, 则四元数$\mathit{\boldsymbol{\hat h}}$称为单位四元数, 而旋转矩阵R与单位四元数间的关系为^[9]：

$ \mathit{\boldsymbol{R}} = \left[ {\begin{array}{*{20}{c}} {h_0^2 + h_x^2 - h_y^2 - h_z^2}&{2{h_x}{h_y} - 2{h_0}{h_z}}&{2{h_x}{h_z} + 2{h_0}{h_y}}\\ {2{h_x}{h_y} + 2{h_0}{h_z}}&{h_0^2 - h_x^2 + h_y^2 - h_z^2}&{2{h_y}{h_z} - 2{h_0}{h_x}}\\ {2{h_x}{h_z} - 2{h_0}{h_y}}&{2{h_y}{k_z} + 2{h_0}{h_x}}&{h_0^2 - h_x^2 - h_y^2 + h_z^2} \end{array}} \right] $

(2)

则旋转参数可表示为：

$ \left\{ {\begin{array}{*{20}{l}} {\tan \alpha = - {\mathit{\boldsymbol{R}}_{13}}/{\mathit{\boldsymbol{R}}_{33}}}\\ {\sin \beta = - {\mathit{\boldsymbol{R}}_{23}}}\\ {\tan \gamma = {\mathit{\boldsymbol{R}}_{21}}/{\mathit{\boldsymbol{R}}_{22}}} \end{array}} \right. $

(3)

2 依据点面关系构建点云配准模型 2.1 点特征点云配准数学模型

点云配准是一种刚体变换^[10], 其形式为：

$ {\left[ {\begin{array}{*{20}{l}} X\\ Y\\ Z \end{array}} \right]_Q} = \mathit{\boldsymbol{R}}{\left[ {\begin{array}{*{20}{l}} X\\ Y\\ Z \end{array}} \right]_P} + \left[ {\begin{array}{*{20}{l}} {\Delta X}\\ {\Delta Y}\\ {\Delta Z} \end{array}} \right] $

(4)

式中, X_Q、Y_Q、Z_Q为基准测站Q坐标系下的坐标, X_P、Y_P、Z_P为待配准测站P坐标系下的坐标, R为a_i、b_i、c_i(i= 1, 2, 3)等元素组成的旋转矩阵, ΔX、ΔY、ΔZ为待配准测站P坐标系转换到基准测站Q坐标系下的平移参数。

2.2 线基元点云配准模型

如图 1所示, 从待配准测站P坐标系和基准测站Q坐标系下的点云数据上提取直线段AB和CD(AB、CD为同名直线段, 但端点不一定为同名点), 其中点A和点B在P坐标系下的坐标为(X_PA, Y_PA, Z_PA)和(X_PB, Y_PB, Z_PB), 点A和点B在Q坐标系下的坐标(均按式(4)表示)为(X_QA, Y_QA, Z_QA)和(X_QB, Y_QB, Z_QB), 点P在Q坐标系下的坐标为(ΔX, ΔY, ΔZ), 点C和点D在Q坐标系下的坐标为(X_QC, Y_QC, Z_QC)和(X_QD, Y_QD, Z_QD)。若点云数据完成配准, 则平面PAB与平面QCD完全重合, 即点A和点B到平面PCD的距离应为0。基于此, 利用单位四元数来描述基于线基元的点云配准模型。

平面QCD的法向量n(n_Qx, n_Qy, n_Qz)(用Q坐标系下的坐标表示)为：

$ \mathit{\boldsymbol{n}} = \mathit{\boldsymbol{PC}} \times \mathit{\boldsymbol{PD}} $

(5)

式中, 向量PC=(X_QC-ΔX, Y_QC-ΔY, Z_QC-ΔZ), 向量PD=(X_QD-ΔX, Y_QD-ΔY, Z_QD-ΔZ)。

设平面PCD的一般方程为：

$ Ax + By + Cz + D = 0 $

(6)

式中, A=n_Qx, B=n_Qy, C=n_Qz。

由于点P(ΔX, ΔY, ΔZ)在平面PCD上, 故可得D=－n_QxΔX－n_QyΔY－n_QzΔZ, 且点A和点B也在平面PCD上, 即满足式(6), 整理得：

$ \left\{ {\begin{array}{*{20}{l}} {{F_A} = {n_{Qx}}{{\bar X}_{QA}} + {n_{Qy}}{{\bar Y}_{QA}} + {n_{Qz}}{{\bar Z}_{QA}} = 0}\\ {{F_B} = {n_{Qx}}{{\bar X}_{QB}} + {n_{Qy}}{{\bar Y}_{QB}} + {n_{Qz}}{{\bar Z}_{QB}} = 0} \end{array}} \right. $

(7)

式中, X_QA=a₁X_PA+a₂Y_PA+a₃Z_PA, Y_QA=b₁X_PA+b₂Y_PA+b₃Z_PA, Z_QA=c₁X_PA+c₂Y_PA+c₃Z_PA, X_QB=a₁X_PB+a₂Y_PB+a₃Z_PB, Y_QB=b₁X_PB+b₂Y_PB+b₃Z_PB, Z_QB=c₁X_PB+c₂Y_PB+c₃Z_PB。

2.3 平差模型解算

将式(7)按Taylor展开, 取一次项进行线性化, 可得关于点A和点B到平面PCD的距离误差方程式：

$ \left\{\begin{array}{c} v_{A}=b_{11} \mathrm{d} h_{0}+b_{12} \mathrm{d} h_{x}+b_{13} \mathrm{d} h_{y}+b_{14} \mathrm{d} h_{z}+ \\ b_{15} \mathrm{d} \Delta X+b_{16} \mathrm{d} \Delta Y+b_{17} \mathrm{d} \Delta Z-l_{A} \\ v_{B}=b_{21} \mathrm{d} h_{0}+b_{22} \mathrm{d} h_{x}+b_{23} \mathrm{d} h_{y}+b_{24} \mathrm{d} h_{z}+ \\ b_{25} \mathrm{d} \Delta X+b_{26} \mathrm{d} \Delta Y+b_{27} \mathrm{d} \Delta Z-l_{B} \end{array}\right. $

(8)

式中, ${b_{11}} = \frac{{\partial {F_A}}}{{\partial {h_0}}}, {b_{12}} = \frac{{\partial {F_A}}}{{\partial {h_x}}}, {b_{13}} = \frac{{\partial {F_A}}}{{\partial {h_y}}}, {b_{14}} = \frac{{\partial {F_A}}}{{\partial {h_z}}}, {b_{15}} = \frac{{\partial {F_A}}}{{\partial \Delta X}}, {b_{16}} = \frac{{\partial {F_A}}}{{\partial \Delta Y}}, {b_{17}} = \frac{{\partial {F_A}}}{{\partial \Delta Z}}, {b_{21}} = \frac{{\partial {F_B}}}{{\partial {h_0}}}, {b_{22}} = \frac{{\partial {F_B}}}{{\partial {h_x}}}, {b_{23}} = \frac{{\partial {F_B}}}{{\partial {h_y}}}, {b_{24}} = \frac{{\partial {F_B}}}{{\partial {h_z}}}, {b_{25}} = \frac{{\partial {F_B}}}{{\partial \Delta X}}, {b_{26}} = \frac{{\partial {F_B}}}{{\partial \Delta Y}}, {b_{27}} = \frac{{\partial {F_B}}}{{\partial \Delta Z}}, {l_A} = - F_A^0, {l_B} = - F_B^0$。

若从点云数据中提取n对直线, 则根据式(8)可列出2n个误差方程式, 并表示为矩阵形式：

$ \mathit{\boldsymbol{V}} = \mathit{\boldsymbol{B\hat X}} - \mathit{\boldsymbol{L}} $

(9)

式中, V为2n×1的残差向量, B为2n×7的系数矩阵, L为2n×1的常数项矩阵。其法方程为：

$ {\mathit{\boldsymbol{N}}_{bb}}\mathit{\boldsymbol{\hat X}} - \mathit{\boldsymbol{W}} = {\rm{0}} $

(10)

式中, ${\mathit{\boldsymbol{N}}_{bb}} = {\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{PB}}, \mathit{\boldsymbol{W}} = {\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{PL}}$。

由于单位四元数$\mathit{\boldsymbol{\hat h}}$还满足‖$\mathit{\boldsymbol{\hat h}}$‖²=1, 对其线性化得：

$ {\mathit{\boldsymbol{C}}_x}\mathit{\boldsymbol{\hat X}} + {\mathit{\boldsymbol{W}}_x} = 0 $

(11)

式中, C_x=[2h₀ 2h_x 2h_y 2h_z 0 0 0], W_x=h₀²+h_x²+h_y²+h_z²－1。

联立式(10)和式(11), 至少需要从待配准测站和基准测站的点云数据中提取4对直线段, 依据附有限制条件的间接平差原理^[11]求解旋转参数和平移参数。具体实施过程为：

1) 确定参数的初始值：h₀ = 1, h_x=h_y=h_z=0, 平移参数ΔX = ΔY= ΔZ =0。

2) 将h₀、h_x、h_y、h_z的值代入式(2), 求出旋转矩阵R。

3) 计算参数改正数$\mathit{\boldsymbol{\hat X}}$的值, 更新h₀、h_x、h_y、h_z的值：h₀=h₀+dh₀, h_x=h_x+dh_x, h_y=h_y+dh_y, h_z=h_z+dh_z, ΔX=ΔX+dΔX, ΔY=ΔY+dΔY, ΔZ =ΔZ+dΔZ。

4) 判断参数改正数是否均小于限定的阈值：若max{dh₀, dh_x, dh_y, dh_z}>ε₁=10^－6或者max{dΔX, dΔY, dΔZ}>ε₂=10^－3, 则需要迭代, 重复步骤2)~4), 直至满足要求; 否则, 迭代结束。

5) 迭代终止后, 得到最终的h₀、h_x、h_y、h_z代入式(2), 计算旋转矩阵R, 由式(3)确定旋转角度α、β、γ。

3 实验结果及分析

本文设计3个实验以验证点云配准算法的正确性及实用性。

3.1 实验1

采用文献[12]中的房屋结构点云数据作为数据源, 挑选不同顶点的坐标数据构成直线, 且待配准测站(P坐标系)中的直线段端点与基准测站(Q坐标系)的两端点并不一定完全对应。P坐标系的直线特征的坐标数据见表 1, 本实验模拟了小转角和大转角2种实验方案的配准参数, 具体见表 2, 表 3为本文算法与文献[4]算法的解算结果对比。

由表 3可知, 对于小转角(< 3°)和大转角的算例, 本文算法的解算结果与模拟值较为吻合, 证明本文算法是正确有效的。文献[4]算法计算的α最大偏差为0.000 8°, β最大偏差为0.000 2°, γ最大偏差为0.000 7°, ΔX最大偏差为0.000 4 m, ΔY最大偏差为0.000 2 m, ΔZ最大偏差为0.000 1 m; 而本文算法计算的α最大偏差为0.000 2°, β最大偏差为0.000 1°, γ最大偏差为0.000 2°, ΔX最大偏差为0.000 1 m, ΔY最大偏差为0.000 1 m, ΔZ最大偏差为0.000 0 m。由此可知, 本文算法对于模拟数据的计算精度高于文献[4]算法。虽然本文为迭代算法, 但迭代次数均在10次以内, 特别是对于大转角, 并未出现迭代不收敛的现象, 因为本文的点云配准算法是基于单位四元数进行的, 增强了算法的稳定性和适应能力。

在图 2和3中, 虚线代表待配准测站上的直线, 实线代表基准测站上的直线。由图可知, 在配准前, 小转角数据和大转角数据下待配准测站和基准测站上的同名直线并不重合, 而使用本文基于直线特征的点云配准算法后, 实线和虚线的重合度较高, 说明本文配准算法是可行的。

3.2 实验2

为进一步验证本文算法的可靠性及实用性, 与文献[6]算法进行对比分析。实验采用文献[5]中实测的建筑物立面的LiDAR点云数据, 其中从待配准测站和基准测站中共提取7对同名直线数据, 限于篇幅数据并未完全列出, 配准前后直线的空间分布见图 4。

本文算法与文献[6]算法的配准结果见表 4, 从解算的旋转矩阵和平移参数来看, 结果比较接近。表 5为配准精度对比, 由表可知, 配准后得到的直线方向向量偏差的中误差相差不大, 本文算法的矩偏差中误差均优于文献[6]算法。结合图 4可以看出, 配准后基准测站直线和待配准测站直线重合, 证明本文算法能够正确实现配准。

3.3 实验3

实验3利用三维激光扫描仪获取太原理工大学虎峪校区训练馆前的一个标志物的点云数据, 共扫描4站点云, 对点云数据进行精简去噪后, 选择第1测站数据作为基准测站数据(图 5(a)), 第4测站的点云数据作为待配准测站数据(图 5(b)), 共选取4对同名直线, 空间分布见图 6。

由图 6可以看出, 从基准测站和待配准测站提取的同名直线端点并不是同名点。对比本文算法与文献[4]算法, 配准结果和配准精度见表 6和7。由表 6、7可以看出, 2种算法的配准参数比较接近, 文献[4]算法与本文算法计算的单位方向向量偏差中误差都为0.000 1 m, 而矩向量偏差分别为0.015 8 m和0.013 1 m, 本文算法稍优于文献[4]算法。主要是因为文献[4]算法是先计算旋转矩阵, 后计算平移向量, 平移向量的计算依赖于旋转矩阵; 而本文算法实现了旋转矩阵和平移向量的同步求解, 从而提高了计算精度。由图 7可以看出, 依据本文算法计算的配准参数, 基准测站和待配准测站的点云数据能较好地拼接在一起。

4 结语

从应用方面来看, 本文推导的线特征点云配准模型能够克服传统的基于点特征的配准算法不能解决缺少完全对应的同名特征点不足的问题, 且不要求同名直线段的端点具有一致性。同时, 由于采用单位四元数进行描述, 线性化程度较高, 对于大转角也能较快收敛, 对初值的依赖性也降低。从实验结果来看, 该模型具有一定的实用价值。