0 概述

近年来, 随着传感器技术及相关算法的快速发展, 激光雷达测量技术在不同领域得到广泛应用, 如激光扫描仪以自动高效的方式提供精确的测量服务, 激光摄像机收集大量的光谱信息数据。激光扫描仪数据和数码相机数据的特性互补, 这对于集成数据处理方案的设计具有重要价值。

目前, 国内外学者提出了不同的方法来配准扫描仪和相机中的2种异源数据。文献[1]利用图像建模的方法, 基于单个图像进行激光扫描仪和未校准数码相机的自动配准, 其中, 通过2D图像和3D点云对应自动解决配准的问题, 2D和3D点云通过距离和角度的空间不变量搜索进行匹配。但是, 文献[1]方法必须考虑若干输入因素, 如特殊目标、消失点和几何约束等, 导致其算法复杂度高, 时间效率低。文献[2]针对研究对象的编码目标提出一种配准策略, 该目标由校准的数码相机记录, 刚性连接到激光扫描仪中, 使用全自动的处理方法解决激光扫描仪坐标系内相机的空间位置和方向问题, 被识别的编码目标则被应用到3D相似性变换中。然而, 文献[2]方法需要一台连接到激光扫描仪的相机和放置在研究对象上的一些编码目标, 导致算法的鲁棒性较差。文献[3]提出一种激光雷达和相机系统的自动、无目标、外部校准算法, 其适用于现场校准, 并且使用基于互信息(Mutual Information, MI)的框架来交叉记录由相机和激光模态测量的强度和反射率信息, 但是, 该算法对硬件的要求高, 时间效率低。文献[4]提出一种基于MI的车载激光雷达点云数据与全景图配准算法, 通过使用统一的球面全景成像模型, 将车载激光雷达点云数据生成的深度图与MI提取的全景图进行匹配。但是, 该方法采用球面全景成像模式, 导致其配准精度较低。

上述针对2种异源数据配准的研究并未对多幅局部点云与多张纹理影像的自动配准进行深入探讨。为了解决地面激光扫描仪在大场景扫描中无法一次性采集到完整点云数据, 以及点云数据与单张纹理影像配准难以得到完整RGB信息的问题, 本文提出一种多幅点云数据与纹理序列间的自动配准方法。对局部点云进行3D-3D配准以得到完整的点云数据, 采用中心投影方式生成强度图像, 利用SIFT算法将强度图像与矫正过的纹理影像进行匹配, 在此基础上, 对纹理影像进行融合预处理操作, 以解决配准过程中的纹理重叠问题并保证算法的稳健性和可靠性。

1 研究方法 1.1 方法步骤

如图 1所示, 本文方法的步骤如下:

1) 进行多幅局部点云数据的预处理及局部点云数据的配准。

2) 采用中心投影的方法将配准后的点云投影生成强度图像并对强度图像进行处理。

3) 对纹理序列进行畸变矫正并实现强度图像与纹理序列的特征匹配与优化。

4) 采用融合预处理方法解决多张纹理影像与强度图像匹配过程中的纹理重叠问题, 通过共线方程得到带有纹理的点云数据, 最终实现多幅点云数据与纹理影像序列间的自动配准。

1.2 多幅点云数据的预处理及配准

为了得到完整的点云, 对采集获得的旋转错位、平移错位以及不完整的局部点云进行配准。因为激光扫描仪采集数据量大、扫描范围广, 为了加快实验中的数据处理速度, 在点云数据配准之前, 需要对局部点云进行预处理操作以获得更有效的数据, 该过程的主要方法有点云数据的去噪、抽稀和滤波等。

1.2.1 点云数据的预处理

本文利用PCL(Point Cloud Library)中的体素滤波器^[5]对每一幅点云数据进行滤波处理, 即对实验输入的三维局部点云数据创建一个三维的体素栅格, 然后在每个体素栅格内用体素中可以近似表示其他点的重心点来描述该体素内的所有点。相对于其他方法, 该方法不仅减少了点云数据, 还保持了点云的形状特性, 可以更准确地表示采样点对应的曲面。本文实验对采集到的3幅局部点云进行研究, 其中, 第1幅局部点云在经过体素滤波前后的效果如图 2所示。

1.2.2 点云配准

多幅点云数据配准是计算机技术中一个重要问题, 迭代最近点(Iterative Closest Point, ICP)是图形图像领域的经典算法, 在科学研究中使用较为广泛。但是, 经过滤波去噪后的点云数据仍然十分复杂, 仅使用ICP算法会陷入局部最优。因此, 本文采用NDT(Normal Distributions Transform)与ICP相结合的方法进行多幅点云匹配。

NDT算法^[6]将参考点云所占空间划分成指定大小的网格或者体素, 并计算每个网格的多维正态分布参数, 利用这些参数计算出配准点云变换到参考点云网格中的变换矩阵, 然后将变换矩阵作为一个良好的迭代初值提供给ICP。ICP算法^[7]则利用NDT算法得到的初值, 经过多次迭代得到最优的匹配结果。配准后的效果如图 3所示。

1.3 强度图像的生成及处理

点云数据与纹理影像的自动配准主要通过地面激光扫描仪的反射强度图像与纹理的匹配来完成。虽然地面激光扫描仪大多采用球面投影^[8]方法获取三维点云数据, 但是, 直接利用球面投影方法生成强度图像, 会在之后的特征匹配中造成较大的形变, 最终导致2种异源数据配准失败。而中心投影^[9]方法和数码相机的投影方法相同, 减少了2种非同源数据间的差别, 提高了配准成功概率, 因此, 本文采用中心投影的方法生成强度图像。

1.3.1 基于中心投影的强度图像生成

为了实现3D点云的2D表示, 采用中心投影的方式将每个激光扫描仪点的笛卡尔坐标投影到虚拟图像平面。虚拟摄像机的坐标系定义如图 4所示。

3D相机的坐标系由帧{i, j, k, O}定义, 投影中心则依赖于激光扫描仪坐标系{X, Y, Z, O}的原点, 虚拟图像平面由向量i和k决定。通过具有最小和最大水平角度的激光扫描仪的点来计算摄像机坐标系的Z轴(矢量k), 使得它指向点云的水平中心。计算过程如式(1)所示:

$ \mathit{\boldsymbol{k}} = \frac{{\mathit{\boldsymbol{p}}_\partial ^{\min }}}{{\left| {\mathit{\boldsymbol{p}}_\partial ^{\min }} \right|}} + \frac{{\mathit{\boldsymbol{p}}_\partial ^{\max }}}{{\left| {\mathit{\boldsymbol{p}}_\partial ^{\max }} \right|}}, {\mathit{\boldsymbol{k}}_\mathit{\boldsymbol{Z}}} = 0, \mathit{\boldsymbol{k}} = \frac{{ - \mathit{\boldsymbol{k}}}}{{|\mathit{\boldsymbol{k}}|}} $

(1)

在通常情况下, Y轴的方向被定义为摄像机坐标系(矢量j)的向上方向, 可以将近似向上方向设置为激光扫描仪坐标系的Z轴, 摄像机坐标系的其余帧矢量如式(2)所示:

$ \begin{array}{*{20}{l}} {\mathit{\boldsymbol{i}} = {{\left( {\begin{array}{*{20}{l}} 0&0&1 \end{array}} \right)}^{\rm{T}}} \times \mathit{\boldsymbol{k}}, \mathit{\boldsymbol{i}} = \frac{\mathit{\boldsymbol{i}}}{{|\mathit{\boldsymbol{i}}|}}}\\ {\mathit{\boldsymbol{j}} = \mathit{\boldsymbol{k}} \times \mathit{\boldsymbol{i}}, \mathit{\boldsymbol{j}} = \frac{\mathit{\boldsymbol{j}}}{{|\mathit{\boldsymbol{j}}|}}} \end{array} $

(2)

利用应用矩阵V使点p=(X, Y, Z)^T从激光扫描仪坐标系变换到3D相机坐标系, 如式(3)所示:

$ \mathit{\boldsymbol{V}} = {\left[ {\begin{array}{*{20}{l}} \mathit{\boldsymbol{i}}&\mathit{\boldsymbol{j}}&\mathit{\boldsymbol{k}}&0\\ 0&0&0&0 \end{array}} \right]^{\rm{T}}}, \mathit{\boldsymbol{\widehat p}} = \mathit{\boldsymbol{V}}{\mathit{\boldsymbol{p}}_h} $

(3)

其中, p_h=(w·X w·Y w·Z w)^T, w=1。

虚拟平面到投影中心点的距离由用户设定的理论距离c决定。在3D相机坐标系中, 通过投影矩阵M来描述点$\mathit{\boldsymbol{\widehat p}}$对于平面的投影, 求解过程如式(4)所示:

$ \mathit{\boldsymbol{\tilde p}} = \left( {\begin{array}{*{20}{c}} {\tilde x}\\ {\tilde y}\\ {\tilde z}\\ {\tilde w} \end{array}} \right) = \mathit{\boldsymbol{MV}}{\mathit{\boldsymbol{p}}_h} = \mathit{\boldsymbol{M}}\mathit{\boldsymbol{\widehat p}} = \left[ {\begin{array}{*{20}{c}} 1&0&0&0\\ 0&1&0&0\\ 0&0&1&0\\ 0&0&{\frac{{ - 1}}{c}}&0 \end{array}} \right] \cdot \left( {\begin{array}{*{20}{c}} {\hat x}\\ {\hat y}\\ {\hat z}\\ {\hat w} \end{array}} \right) $

(4)

欧几里德图像坐标的结果p′来自归一化齐次坐标向量${\mathit{\boldsymbol{\tilde p}}}$, 如式(5)所示:

$ {\mathit{\boldsymbol{p}}^\prime } = {\left( {\begin{array}{*{20}{c}} {\tilde x}&{\tilde y}&{\tilde z}&{\tilde w} \end{array}} \right)^{\rm{T}}} \cdot \frac{1}{{\tilde w}} = {\left( {\begin{array}{*{20}{c}} {{x^\prime }}&{{y^\prime }}&{ - c}&1 \end{array}} \right)^{\rm{T}}} $

(5)

通过对投影坐标的最小值和最大值的计算得到图像的大小, 或者使用用户自定义的像素大小Δpx计算图像高度I_h和宽度I_w, 结果如式(6)所示:

$ \begin{array}{*{20}{l}} {{I_w} = \left( {{{x'}_{\max }} - {{x'}_{\min }}} \right)/\Delta px}\\ {{I_h} = \left( {{{y'}_{\max }} - {{y'}_{\min }}} \right)/\Delta px} \end{array} $

(6)

最后, 计算配准后点云数据左上角的点的投影坐标, 如式(7)所示:

$ \begin{array}{*{20}{l}} {{x_p} = \left( {x' - {{x'}_{\min }}} \right)/\Delta px}\\ {{y_p} = {I_k} - \left( {y' - {{y'}_{\min }}} \right)/\Delta px} \end{array} $

(7)

1.3.2 强度图像处理

本文方法在生成强度图像后, 采用消隐算法解决同一束光线上存在2个或多个三维激光点投影在同一个图像格网上的问题, 即遍历生成的强度图像, 距离扫描仪扫描中心偏近的原始激光点作为该图像格网的投影点, 并且记录每一个格网点对应的三维激光点坐标。另外, 本文采用反权距离插值(Inverse Distance Weighted, IDW)^[10]方法来解决生成的强度图像中有像素值为0的问题。

1.4 强度图像和纹理序列的特征匹配及优化评估

在图像配准前, 需要对纹理影像进行畸变矫正处理^[11-12]。畸变矫正处理本质上是对一幅退化图像进行恢复, 包括空间坐标变换、像素灰度级确定2个主要部分。空间坐标变换是使像素落在正确的位置, 灰度级确定则是对纹理影像进行视觉上的修正。在两者的特征匹配过程中, 为了降低错误匹配对对配准结果的影响, 本文对匹配结果进行优化与评估处理。

1.4.1 特征匹配

本文采用SIFT算法^[13-14]实现强度图像和纹理影像的匹配。首先, 实现SIFT特征检测要通过连续应用高斯低通滤波器构建图像尺度空间金字塔, 在被检测的图像中, 金字塔局部极值被检测为关键点, 使用泰勒展开来计算它们的准确位置; 随后, 分析所选关键点的2个稳定性标准, 如果校正位置处的泰勒展开小于峰值阈值, 则丢弃具有低对比度的关键点, 相反, 则丢弃位于边缘上的极值。最后, 通过手动优化这些阈值来检索最佳候选对应集。部分特征检测结果如图 5所示。

1.4.2 匹配优化与评估

本文引入极线几何^[15-16]并采用RANSAC算法对SIFT提供的候选对应点进行评估与优化。从不同摄像机和视角拍摄的相同场景的图像与它们的极线几何相关, 3×3的同质矩阵F是2个图像中对应点之间映射的代数表示。考虑到一对同质图像的坐标p₁=(x₁, y₁, l)^T∈I₁和p₂=(x₂, y₂, l)^T∈I₂, 两者都表示同质世界坐标下的点Q=(X, Y, Z, l)^T的投影, 如式(8)所示:

$ \begin{array}{l} {\left( {x_1^\prime , y_1^\prime , w_1^\prime } \right)^{\rm{T}}} = {\mathit{\boldsymbol{P}}_1}\mathit{\boldsymbol{Q}}, {x_1} = x_1^\prime /w_1^\prime , {y_1} = y_1^\prime /w_1^\prime \\ {\left( {x_2^\prime , y_2^\prime , w_2^\prime } \right)^{\rm{T}}} = {\mathit{\boldsymbol{P}}_2}\mathit{\boldsymbol{Q}}, {x_2} = x_2^\prime /w_2^\prime , {y_2} = y_2^\prime /w_2^\prime \end{array} $

(8)

其中, P_i=K_i[R_i|t_i], P_i是点p_i的投影矩阵, K_i是相机的内参, R_i、t_i分别为点p_i的旋转、平移矩阵。在世界坐标系中的摄像机位置为C₁, t₁=－R₁C₁将点从世界坐标系变换到摄像机坐标系, 如果p₁和p₂确实显示相同的场景点P, 则p₂位于表示p₁的极线的同质线I₂上, 如图 6所示, 其中, I₂由式(9)表示:

$ {I_2} = {C^\prime } \times p_2^\prime $

(9)

其中, p′₂=P₂P⁺p₁, C′=P₂C₁∈I₂, P⁺是P₁的广义逆矩阵。相反, 如果p₁和p₂显示的是不同的场景点, 式(10)则表示以|I₂|为单位的p₂与I₂的距离, 其中, 基本矩阵F具有7个自由度, 且det(F)=0。

$ \mathit{\boldsymbol{p}}_2^{\rm{T}}{I_2} = \mathit{\boldsymbol{p}}_2^{\rm{T}}{\mathit{\boldsymbol{P}}_2}{C_1} \times {\mathit{\boldsymbol{P}}_2}{\mathit{\boldsymbol{P}}^ + }{\mathit{\boldsymbol{p}}_1} = \mathit{\boldsymbol{p}}_2^{\rm{T}}\mathit{\boldsymbol{F}}{\mathit{\boldsymbol{p}}_1} = 0 $

(10)

图 6(a)、图 6(b)中的直线表示极线, 白色粗正交线表示在假定的F下产生的误差e₅₆。

在RANSAC算法^[17]中, 将F作为模型参数, 目的是识别正确的匹配子集及其相关的F。RANSAC依次选择7个匹配(共识集), 并使用式(10)中的条件拟合模型矩阵F对其余部分进行评估, RANSAC是一种迭代的概率优化算法, 该算法会在完成预定义的迭代次数或者模型很好地拟合预定数量的匹配之后终止。误差函数E(F)通过模型F来描述匹配对p₁, p₂)的准确率, 如式(11)所示:

$ E(\mathit{\boldsymbol{F}}) = \mathit{\boldsymbol{p}}_2^{\rm{T}}\mathit{\boldsymbol{F}}{\mathit{\boldsymbol{p}}_1} $

(11)

图 6(c)中的直方图表示候选对应集中的误差值分布。

1.5 纹理序列的融合预处理及自动配准

采用多张纹理影像的主要目的是使点云数据获得完整的RGB信息, 解决单张纹理影像不能满足大场景点云数据对RGB信息的需求问题。然而, 在纹理影像序列与强度图像的配准过程中存在纹理重叠问题, 因此, 本文采用融合算法^[18]对纹理影像序列进行预处理操作。采用SIFT算法对纹理序列进行特征检测, 得到对应的不变矩阵, 利用加权融合方法得到加权系数, 再通过加权系数对待处理的纹理序列进行融合预处理。本文的融合预处理并非文献[18]中真正意义上实现2张纹理影像的拼接, 而是极大地减少纹理影像序列之间的重叠, 进而解决在纹理影像序列与强度图像配准过程中存在的纹理重叠问题^[19]。

将特征匹配得到的参数值代入共线方程^[20]中, 即将纹理影像序列的纹理信息赋给点云数据, 实现点云数据与纹理序列的配准。共线方程如式(12)所示:

$ \begin{array}{*{20}{l}} {x - {x_0} = - f\frac{{{\mathit{\boldsymbol{r}}_{11}}\left( {X - {X_0}} \right) + {\mathit{\boldsymbol{r}}_{21}}\left( {Y - {Y_0}} \right) + {\mathit{\boldsymbol{r}}_{31}}\left( {Z - {Z_0}} \right)}}{{{\mathit{\boldsymbol{r}}_{13}}\left( {X - {X_0}} \right) + {\mathit{\boldsymbol{r}}_{23}}\left( {Y - {Y_0}} \right) + {\mathit{\boldsymbol{r}}_{33}}\left( {Z - {Z_0}} \right)}}}\\ {y - {y_0} = - f\frac{{{\mathit{\boldsymbol{r}}_{12}}\left( {X - {X_0}} \right) + {\mathit{\boldsymbol{r}}_2}\left( {Y - {Y_0}} \right) + {\mathit{\boldsymbol{r}}_{32}}\left( {Z - {Z_0}} \right)}}{{{\mathit{\boldsymbol{r}}_{13}}\left( {X - {X_0}} \right) + {\mathit{\boldsymbol{r}}_{23}}\left( {Y - {Y_0}} \right) + {\mathit{\boldsymbol{r}}_{33}}\left( {Z - {Z_0}} \right)}}} \end{array} $

(12)

其中, x、y是已知的图像坐标, X、Y、Z表示构成对象的相应的已知激光扫描仪坐标, r_ij是组成未知的3×3旋转矩阵的元素, X₀、Y₀、Z₀表示未知的相机位置, f表示主焦距。

2 实验结果与分析

本文利用VS2015开发平台并结合PCL类库、OPENCV类库编程进行实验。采用Trimble公司生产的TX6三维地面激光扫描仪获取点云数据, 实验场景是中北大学教学楼主楼。实验总共采集到3幅局部点云, 总计7 306 818个点, 局部点云数据效果如图 7所示, 点云配准结果见前文图 3。

为了验证采用中心投影的方法生成的强度图像与纹理影像匹配点对的准确性与可靠性, 本文还对利用球面投影生成的强度图像和纹理影像的匹配进行了实验。2种投影方法生成的强度图像对比如图 8所示。

第1次实验将2种投影方法生成的强度图像分别直接和纹理序列1进行匹配, 匹配结果如图 9所示。第2次实验将2种投影方法生成的强度图像分别直接和纹理序列2进行匹配, 匹配结果如图 10所示。表 1所示为2种投影方法特征匹配的结果对比, 从表 1可以看出, 与球面投影方法相比, 采用中心投影方法得到的强度图像与纹理序列匹配的可靠性和稳定性更高, 中心投影方法在很大程度上减少了2种异源数据的差异性, 提高了2种数据自动配准的准确率。

球面投影生成的强度图像与纹理影像匹配效果较差, 将其得到的参数代入共线方程后得不到完整的带RGB信息的点云数据, 而采用本文方法则可以获得很好的配准效果, 如图 11所示。

如表 2所示, 将本文方法与文献[4, 20]方法的配准结果进行对比。从表 2可以看出, 本文方法缩短了计算时间, 提高了算法的执行效率, 均方根误差为3.4个像素, 配准精度得到提升。

综合表 1、表 2及图 11可以看出, 在多幅点云数据与纹理序列间的配准问题中, 本文方法可以获得较好的配准结果。

3 结束语

本文提出一种点云数据与纹理序列的自动配准方法。为解决ICP配准方法容易陷入局部最优的问题, 结合NDT与ICP算法来配准局部点云, 利用SIFT算法对纹理影像和配准后的点云进行特征匹配, 同时采用融合预处理方法解决纹理重叠问题。在此基础上, 通过共线方程实现多幅点云数据与纹理序列的自动配准。实验结果验证了该方法的高效性。在点云数据转换成强度图像的过程中存在数据丢失的现象, 对本文方法进行改进以解决该问题将是下一步的研究方向。