ABSTRACT

Terrain Laser Scanning (TLS) plays an important role in engineering survey, 3D reconstruction and cultural relics preservation. The way laser scanner acquiring point cloud data must be several stations attaching because of views interruption, therefore, registration of point cloud is necessary. This article put forward an automatic registration method for TLS, extracts patches in point cloud firstly, counting transform parameters using patches combination and voting, and output the transformed point cloud automatically in the end. Experiments shows that our method can finish the registration of multiple-station city environment TLS robustly and efficiently, which makes the method useful in practice.

Keywords:Laser Scanning, Registration, Patches, Combination, Voting

基于面组合投票的地面站点云配准方法

周桐，梁福逊，杨必胜

武汉大学测绘遥感信息工程国家重点实验室，湖北 武汉

收稿日期：2017年6月10日；录用日期：2017年7月2日；发布日期：2017年7月5日

摘 要

地面站三维激光扫描在工程测量，三维重建，文物保护等领域发挥越来越重要的作用。地面激光扫描仪由于视角和遮挡的问题，必须通过多站采集拼接的方式获取数据，故配准是必不可少的步骤。文章提出了一种自动化激光点云配准方案，从点云中提取面片，通过面组合和最佳投票的方法计算转换参数，并自动完成多站点云的坐标转换与输出。经实验证实，本文的配准方法能够稳健且高效地完成城市环境点云的配准，有很好的实用价值。

关键词 :激光扫描，配准，面片，组合，投票

Copyright © 2017 by author and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

地面站LiDAR由于其高精度、高密度和非接触式采集的特点，能够在较短时间内获得场景表面详细的三维信息，适合多种环境下的工程测量与地理信息分析，在城市区域测图、地下管线测量，古文物重建及保护等多种场合发挥越来越重要的作用。由于站载激光点云的扫描特点，单站点云不可避免存在视角遮挡的问题，必须通过多站采集和拼接的方式获取一片区域的完整数据，所以，配准是激光点云数据处理中必不可少的步骤。

目前工程较多采用人工交互的方法完成地面站LiDAR配准，包括靶标方案和手工选取同名特征的方案。靶标方案需要在数据采集过程中架设多个靶标并保证通视，通过同名靶标对应来求解转换参数以实现配准，费时费力，影响了点云配准的效率，且不能处理多期数据的配准问题。手工选取同名特征方案则是在数据后处理过程中，人工寻找并指定重叠区域的同名点(或同名直线、平面)，求解几何转换参数，完成配准，灵活性高，但需要大量人工操作，影响效率，文献 [1] [2] [3] 等均采用该类方案。

自动化地面站点云配准需要解决两个问题：任意站两两配准以及多站转换参数的优化及输出。三维点云配准最常见的算法是迭代最近点算法(ICP)及其延伸算法，该方法对每个点统计最近邻域内两站点之间的偏移，整体纠正偏移量，迭代计算转换参数，全过程自动且求解结果精度高，但算法对迭代初值有要求，容易陷入局部不收敛的情况。所以，位置姿态信息未知情况下的粗配准，才是点云自动化配准的真正难点。粗配准与精配准的结合，是点云自动化配准最常见的方案。

关于两两自动化粗配准的方法国内外有不少的研究，大致可以分为三类。第一类方法是借助图像配准的手段完成点云配准：将三维点云转换成二维全景深度图，进行图像配准，再将二维转换关系反对应到三维空间，实现三维点云的配准。文献 [4] 和文献 [5] 都采用此类方法，该类方法对深度图像的视角选取有一定的要求，且基于深度图像的特征描述子获得的匹配关系往往不够稳定，对于任意架站的激光点云的配准不能让人满意。第二类方法是基于点特征的配准：基于点特征描述子，通过适当的同名点匹配策略实现配准。Aiger D等 [6] 提出了4PCS (4 points congruent set)的方法，通过创新的搜索策略，实现点云中4对同名点对应关系的遍历查找，Theiler P W等 [7] [8] 提出了基于特征点的K-4PCS的方案，提高了点云配准的适用性。另外，Li N, Cheng P [9] 等提出了在地形表面上提取特征描述完成配准的算法，彭博 [10] 提出了基于曲率图的配准方法。该类方法能自动完成两站点云的配准，但由于激光点云不同的场景环境和剧烈变化的点密度，很难找到稳定的描述子保证点云配准的鲁棒性。第三类方法是基于几何基元的配准：选择各种几何元素(直线、曲线、平面等)作为基元，计算匹配关系和转换参数以完成配准。文献 [11] 提出了一种基于建筑物立面组合关系的机载点云和地面点云匹配的方法；文献 [12] 提出了一种带有语义的杆特征元素实现两站点云的配准，并自动分析多站匹配关系，完成多站配准；文献 [13] 提出了曲线特征的数学模型描述，通过曲线的局部匹配实现点云的配准，适用于佛像等特殊场景的点云数据；文献 [14] 也提出了类似的几何基元配准方法。相比于点特征，几何基元(杆、线、面)在配准流程中更加稳定可靠，成功率更高，适用场合更广泛。如何自动且高效地提取几何基元，寻找对应关系，以及对不同类型几何基元的转换模型求解，是该类方法的难点所在。

文章提出了一种主要适用于城市环境的地面站激光点云的自动配准方法。该方法从原始点云中自动提取面片特征，通过遍历组合和最佳投票的策略，自动求解两站点云刚性变换参数，最后根据两两转换关系一次性输出全部站的点云。该方法结合了激光点云面片提取算法，用于配准基元的获取，创新性地利用投票策略，求得契合关系最佳的刚性变换参数，且能够自动处理多站点云的复杂匹配关系。多组数据实验表明，本文方法能够实现包含人工建筑物场景的点云配准，算法表现稳定，自动化程度高。

2. 点云配准方法

2.1. 综述

本文提出了一种基于平面特征的配准方法，分为两个部分：两两配准(流程图如图1所示)和多站点云输出。首先，通过区域生长的算法，自动提取点云中的全部平面，记录平面参数；其次，求解基于平面元素的转换模型，搜索候选组合，投票得到最优刚性变换参数；最后，分析站与站之间的匹配关系，完成全部点云的坐标转换及输出。下面对方法流程展开介绍。

图1. 两站点云配准方法流程图

2.2. 区域生长法面片提取

本文采用区域生长的方法从原始点云中自动提取面片，以面片作为基本几何基元，实现两站点云的配准。面片提取详细方法介绍如下。

首先，对大数据量的原始点云进行格网抽稀，简化点云，保持密度稳定同时降低计算量。然后，对整个场景的点云进行邻域内曲率和法向量分析，即计算邻域点的坐标分布进行主成分分析，过程如公式 (1)和(2)所示：

(1)

其中， (3*3)代表邻域内点坐标的方差协方差矩阵，n为邻域内点的个数，代表邻域内任意一点，代表邻域点的平均坐标。对其求解特征值和特征向量，如公式2所示：

(2)

其中，分别代表协方差矩阵的三个特征值，;和分别对应的三个特征向量。根据主成分分析理论，显然，最小主分量对应的特征向量即为平面的法向量。另外，邻域内曲率计算如公式3所示：

(3)

根据上述计算的曲率和法向量，按照图2所示的区域生长算法完成面片的提取。

这里，生长条件描述为：“该点相对于种子点的法向量夹角小于阈值，且点到平面距离小于阈值”。上述区域生长算法对输入的离散点进行分割，获得平面点的聚类，图3所示的某一城市环境点云数据面片提取结果中，不同平面的点用不同的颜色标示。

对每个面片，根据点集的坐标求解出三维空间内的平面方程：

(4)

这里，，如公式(2)所示， (为平面点的平均坐标)。记录每个平面的参数A，B，C，D。

2.3. 候选面组合的筛选

由于站载LiDAR的水平架站的特点，同一场景内多站之间的转换关系主要表现在xoy二维平面内，只需通过水平平面的平移和旋转，即可实现站载激光点云的配准。因此我们只选取垂直于水平面的立面，将其转换为二维平面内的线段，作为配准的基本基元。该基元在二维平面的数学方程可表示为：

(5)

这里已知(如公式4)， (为面片点在二维空间上的平均坐标)，记录下每一个立面基元对应的二维参数A，B，d以及。

二维平面内确定一组转换参数，需要至少两条直线。因此，只需从两站待配准的点云中分别选择两个立面元素(对应在二维平面是两条直线)作为候选组合，即可完成转换参数的计算。

记两站待配准的点云分别为和，其中为参考，两站的配准问题即可表示为求解到的刚性变换参数。记中的候选组合为和,的候选组合和，则候选组合的选取需要满足以下条件：

图2. 区域生长法提取面片流程

图3. 某场景原始点云及面片提取结果

1) 与不平行，与不平行；

2) 和之间的夹角与和之间的夹角近似相等。

只有满足上述两个条件的候选组合对，才能保证到求解得到的转换参数满足刚性变换条件。以上述条件为指导，遍历和两站中所有的立面组合，将所有符合条件的立面组合对记录下来，保存为候选组合。

2.4. 刚性变换参数计算

对于1.2所得的每一个候选组合，都能有效地计算一组刚性变换参数来描述两站点云的坐标转换，本节从刚性变换模型来阐述转换参数的计算，具体分为两个步骤：平移量计算和旋转角度的确定。

平移量计算相对简单，每一个立面元素在xoy二维平面内都可以用公式(5)表示成一条直线，故同一站的两个立面元素和易求得二维交点，平移量可表示为：

(6)

其中表示基准站的立面元素在二维平面的交点，表示站对应的交点。

对于二维平面内的旋转角度的求解，需要确定候选组合的对应关系。记和两个立面的候选组合分别为和,和，候选组合对的对应关系有两种(对应,对应；或对应,对应)，如图4所示。

通过以下方式确定方向向量：记两个组合对应的二维平面的交点分别为和，四个立面元素的重心点在二维平面的坐标分别为, , ,，则对于四个立面元素，可以唯一确定其在二维空间内的方向向量(如图4中直线箭头所示)：

(7)

分别检验组合对的两种不同对应关系，判断其方向向量夹角是否相同，如果夹角相同，则证明对应关系正确，否则选择另外一种对应关系即可。确定了候选组合对中四个立面元素正确的对应关系，易求得夹角，进而旋转参数可唯一表示，如公式8所示：

(8)

对于任意一组候选组合对，根据公式6和公式8，可获得唯一确定的转换参数，表示两站点云之间的刚性变换。

图4. 候选组合的对应关系分析

2.5. 最佳转换参数的选择

由于每一组候选组合对都可以求得一组刚性变换的转换参数，都能独立表示一次刚性变换。故如何从诸多候选的刚性变换参数中评选出最优的一组转换参数，是本节讨论的内容。本文提出了投票的方法对每一个转换参数进行评价，从中求得一组最佳的转换参数。

对于每一个立面元素，在二维空间内表示为公式(5)所示的直线公式，通过转换参数易求得该直线经过坐标转换后的公式。对于中的每一个立面元素，经过刚性变换，求得其在参考坐标下新的坐标表示。将中每一个立面元素与进行对比，如果公式(5)对应的参数完全符合，则认为来自两站的立面元素在新的坐标系统中契合，记录为1个投票数。遍历判断每一对立面元素是否契合，求得最终投票数总和。投票数代表了经过坐标转换后两组点云对应的立面元素的契合程度，显然，投票数越高，契合程度越好，该转换参数越可信。

每一组候选组合对应的转换参数，都可以按照上述步骤记录投票，求得转换参数对应投票数最高的一组，即为最佳候选组合对对应的最优转换参数。

以此最优转换参数为初始值，对两站点云执行迭代最邻近点(ICP)算法，对粗配准获得的刚性变换参数进行精确纠正。至此，全自动两两配准的过程全部完成。

2.6. 多站点云自动输出

有了站与站之间的两两配准的方法，还需要进行以下工作，才能完成多站点云配准结果的输出：

1) 根据通视关系，指定哪些站之间需要执行两两配准；

2) 指定基准站，将每站的转换参数统一转化到基准站坐标系下；

3) 点云坐标转换和输出。

本文提出了基于深度优先搜索(DFS)的方法，自动推荐最佳基准站，并计算每站相对于基准站的累积转换参数，最终自动输出全部点云。

2.6.1. 连通关系有向图构建

站与站之间的两两匹配关系，在扫描仪架站的时候就已指定，并通过草图或其他方式记录下来。在两两配准步骤中，就是按照站与站之间的连接关系指定哪些站需要执行配准。根据这样的连接关系，自动构建连通关系有向图：顶点代表每个点云架站，有向边代表两站之间的匹配关系，包括是否执行过配准，以及哪一站是参考站。

连通关系有向图是站与站之间的连接关系的数学描述，以此为基础可进行最佳基准站的分析，及转换路径的计算。

2.6.2. 基准站自动选取

假设任意一站为基准站，从该顶点出发，执行深度优先搜索(DFS)算法，遍历到全部顶点，记录下每一个顶点从基准站出发的最短配准路径，同时记录该路径的长度，作为该路径的代价。根据误差传播理论，路径长度越小，误差传播越小。所以依据最小生成树原理，遍历寻找一个最佳的顶点，该顶点作为起始点，使得所有站的路径代价总和达到最小值。该顶点对应的站点即为推荐最佳基准站。

2.6.3. 转换参数的确定及点云输出

基准站可以按照2.6.2节的方法自动选择，也可以手动指定。以基准站坐标系为统一参考坐标系，求解每一站点云相对于基准站的刚性变换参数，即可获得每站转换后的结果。

基于2.6.1节获得的连通关系有向图，根据深度优先搜索(DFS)算法，获得每一站点云到基准站的最短路径，路径中的每一条边代表了一组两两配准的参数，计算累积转换参数，可以获得每一站点云转换到基准站坐标系的刚性变换参数，进而可以自动完成每一站点云的坐标转换及输出。

3. 点云配准实验

3.1. 实验数据及运行环境

本文对三组城市环境的固定站激光点云数据进行配准测试。三组数据都采用Regel VZ-400扫描仪进行数据采集。数据分别来自不同类型的城市环境，采集过程中记录架站的相对位置关系草图，无靶标，不记录架站的位置和姿态信息。所有数据的原始坐标都是扫描仪局部坐标：笛卡尔右手坐标系，坐标原点为扫描仪中心，Z轴与水平面垂直，X轴方向任意。三组数据具体信息如表1所示：

本文全部实验算法采用C++编写，Visual Studio 2013编译环境下运行，硬件采用Intel Xeon E3处理器，3.5GHz，8G内存，单线程运行，无GPU加速。在此环境下对配准算法的运算效率进行分析。

3.2. 两两配准表现

根据图1所示的两两配准的流程，对原始数据进行配准。图5是场景B中某两站的两两配准过程，从上到下依次是原始数据、面片提取结果、最优同名面片(左站右站各两个，分别用绿色和蓝色标示)以及最终的点云拼接结果。

本文对三个场景一共进行了51次两两配准，均成功完成粗配准，精度足够满足ICP的收敛条件。通过ICP算法精化转换参数。

效率方面，三个场景的配准算法表现及效率如表2所示(不包括数据读入和写入耗时)。由表2可知，本文采用的区域生长面片提取算法，不仅面片提取效果稳定可靠，而且对于大数据量的点云数据有很高的计算效率，每站点云在面片提取算法上的耗时不超过5秒，两站点云从原始数据输入，到配准完成，整个流程时间消耗在10秒以内。对于数十站的大场景激光扫描数据，完成配准不过几分钟的时间，充分证明了本文方法效率方面具有很好的实用性。

表1. 三个场景信息

表2. 两两配准算法效率

图5. 两两配准过程示意图

3.3. 多站输出及配准结果

本文为多站点云自动化输出编写独立的交互模块，如图6所示。该模块根据用户指定的两两匹配关系(来自数据采集时的架站记录)，自动分析站与站的连通关系，通过2.6.2节的方法，推荐最佳基准站(也可手工指定基准站)，根据基准站和连通关系图，自动实现每站点云配准参数的换算以及配准后点云文件的输出。

三个场景的点云配准结果如图7所示。由配准结果图可知，本文算法在各种类型城市环境的表现都比较稳定，场景B(第二幅图)所示的环境建筑物密集，遮挡严重，仍能自动实现配准的完整过程。图8展示了场景C中一个建筑物的局部匹配细节效果，不同站的点用不同颜色表示。

图6. 多站输出界面功能

4. 结论与展望

本文提出了一种地面站激光点云全自动配准方法，创新性利用面片的相对位置组合关系解决配准问题，完整实现了从两两配准，到多站点云输出的过程。面片作为一种良好的几何基元，对点云密度变化有很强适应性，面片位置关系组合起来用于配准，能很大程度克服视角遮挡的难题，确保算法的稳健性。

实验证明该方法有以下优势，第一，配准算法对遮挡有较强的适应性，对重叠度没有过高的要求，比起基于点特征自动配准算法有很大的进步，大大提高了数据获取时架站的灵活性；第二，算法效率高，从大文件原始点云中进行平面分割，到站间相对位置匹配，整个计算过程只需数秒时间，有很强的实用性；第三，本文提供了大场景多站点云配准的完整解决方案，自动化程度高，大大减少了人工工作量。

然而，本文方法过于依赖平面元素，仅适用于城市环境的配准工作，对于林区或自然地形等不存在人工建筑物的场景，很难提取到稳定可靠的平面元素，配准效果不佳。多类几何元素(如杆状物、三维曲线、脊线等)共同组合，能够适用于更多场景的配准问题，是可行的改进方向。

综上，本文方法能够很好解决城市环境下站载LiDAR的配准问题，算法稳健，自动化程度高，在激光点云数据处理过程中有很高的实用价值。

图7. 三个场景架站及配准结果示意图

图8. 某建筑物配准的局部效果

致谢

本文工作在杨必胜教授指导和梁福逊硕士帮助下完成，感谢广东省国土资源测绘院提供数据采集的便利，感谢国家科技支撑项目“村镇区域综合防灾减灾信息系统研究及示范”(编号2014BAL05B07)为本文工作的支持，特此致谢。

项目基金

国家科技支撑项目，编号2014BAL05B07。

文章引用

周桐,梁福逊,杨必胜. 基于面组合投票的地面站点云配准方法
A Registration Method for TLS Based on Patches Combination and Voting[J]. 测绘科学技术, 2017, 05(03): 106-118. http://dx.doi.org/10.12677/GST.2017.53013

参考文献 (References)