一种高程相关的机载LiDAR点云配准方法

吴明东，付建红，孟庆祥

(1.河南理工大学 测绘与国土信息工程学院，河南 焦作 454000；2.武汉大学 遥感信息工程学院，武汉 430079)

0 引言

机载激光扫描(light detection and ranging，LiDAR)集全球定位系统、惯性导航系统与激光测量3种技术于一体，能够快速高效地获取地物表面每个采样点的精确三维坐标，广泛应用于数字城市、灾害监测、海岸工程、林业等方面[1]。在对一个区域进行机载激光扫描时，往往需要多航带或多次扫描才能完成，而多次扫描的激光点云需要通过配准转换到统一的坐标系中[2-3]，这也是利用点云进行三维重建中最重要且最困难的步骤，直接影响着三维重建的精度与效率，是后续点云重建的技术基础[4]。

目前，点云配准的算法有很多，主要分为2类。一类是基于点云数据特性建立的特征匹配法，该方法利用物体表面上的几何特征信息，对特征进行提取和描述，通过确定特征之间的对应关系来实现点云配准。刘剑等[5]提出一种快速点特征直方图(FPFH)特征描述与Delaunay三角剖分相结合的特征点配准方法，有效地提升了点云配准的精度和效率。李健等[6]结合当下最热门的深度学习技术，将深度学习用于三维点云局部特征的提取，融合点云配准算法解决了多维度大角度变化的点云配准问题，使深度学习在点云配准领域实现了很好的应用。盛敏等[7]针对RGBD点云数据的配准问题，利用待配准的RGBD点云模型的曲率和颜色特征度定义了一种基于二者特征相似性的点云配准方法。但是这类方法需要点云包含丰富的颜色信息，同时对点云片的曲率有要求。机载LiDAR数据覆盖范围广、点云密度小，覆盖区内地物数据组成单一(主要为地表和建筑物屋顶)，同时缺少颜色信息，对曲率的计算和特征的提取相对困难。第二类是基于点云集之间相互迭代计算的方法。其中比较主流的是迭代最近点 (iterative closest point，ICP)算法，该算法是Besl等[8]于1992年提出的，其基本思想是通过旋转、平移使2个点集之间的欧氏距离最小。文献[9-10]分析了ICP算法的基本原理并用于对地面扫描数据的配准。乔世权等[11]通过定义源数据与目标数据的曲面距离，利用PCA算法计算旋转矩阵的初始值改进了传统的ICP算法，并用于点云配准。针对点云配准中的尺度和收敛速度问题，赵夫群等[12]通过在ICP算法中加入带边界的尺度矩阵，提出了一种改进的尺度迭代最近点算法，解决点云配准中尺度变换的问题。为了解决强噪声及密度不均匀点云高效配准的难题，彭真等[13]提出了一种基于关键点提取与优化的改进的 ICP点云配准算法。基于特征的匹配和基于ICP算法的匹配在地面激光扫描数据，特别是工业测量等方面得到了广泛的应用。对于机载LiDAR条带的配准问题，这些算法的研究和应用相对较少。

结合机载LiDAR点云数据特点，本文提出一种基于机载LiDAR点云高程相关的配准新方法。该方法首先在基准点云中选取某一目标点，通过目标点与周围若干点进行曲面拟合形成拟合曲面；然后在另一片点云建立的搜索区中遍历采样点，比较每个采样点和其邻近点组成的点集在拟合曲面上的拟合高程与点集的真实高程数据之间相关系数，选择相关系数最大的采样点作为同名点；最后通过2片点云中若干同名点对完成点云配准。

1 基于高程相关的机载点云配准

1.1 点云预处理

对于LiDAR点云数据而言，点云数据量越大将增加点云配准的时间、影响点云配准的效率；点云噪声会对点云同名点的寻找产生干扰，导致错误匹配点对[14]。因此，为保证点云配准的可靠性，在对点云进行配准前先对点云进行预处理。点云的预处理包括点云的去噪和精简。K-D树在点云去噪中有着广泛的运用，在进行点云去噪中，对全部点云进行遍历，并逐一查找各点的邻域，计算邻域内各点到中心点的平均距离，当平均距离过大时，可将该点视为噪声点。点云精简的过程中，在整体降低数据量的同时尽可能多地保留一些能够反映环境的表面特征的点位数据。

1.2 曲面拟合模型建立

对于任意机载LiDAR三维点集M中的一点p，定义p的一个邻域R，将邻域R内的所有点进行曲面拟合，拟合的曲面作为这些点所代表的局部真实场景的逼近，称为最邻近曲面。常见的曲面拟合算法有很多[15-19]，本文采用二次曲面拟合模型算法来拟合曲面，如式(1)所示。

(1)

式中：(Xi，Yi，Zi)为邻域R中第i个点的三维空间坐标值；a、b、c、d、e、f为相应的二次曲面拟合系数。

本文主要思想是在基准点云上选择局部点拟合曲面，计算出待配准点云上的局部点云在拟合曲面上的高程，求解该局部点云在拟合曲面上的高程和其真实高程之间的相关系数，并在搜索区域内移动拟合曲面，取相关系数最大时的局部点云作为基准点云的对应点云，进而确定同名点，完成点云配准。因此，如何快速、准确地选取合适的局部拟合曲面对于后续的配准工作有决定意义。选取的拟合曲面应遵循以下原则。

1)参与拟合的点应选择地形起伏较明显或与周围地面有明显区分的地物点。

2)曲面的范围选取不宜过大，必须保证所拟合的曲面为局部小范围曲面。

3)原始点云拟合后的曲面偏移距离要小。参与拟合的点云为局部真实的地面或地物点，若某一点或多点与拟合后的曲面偏移距离较大，则说明该局部点云不满足曲面模型，将重新选择点云。

图1 局部点云查找示意图

局部点云选择如图1所示。用于拟合曲面的局部点云选择过程为:假设点p为点云数据中满足曲面选取原则中的一个点，以该点为中心寻找在周围点云中距离点p一定阈值范围内的k个点，利用选取的k个点拟合二次曲面，并计算每个点的拟合残差，若拟合残差超限，则移动中心点位置，重新拟合曲面，直到生成合适的最邻近曲面。

1.3 高程相关系数求解

在本文中，求2组变量Z和Z′的相关系数ρ(Z，Z′)。相关系数的取值范围为[-1，1]，当ρ(Z，Z′)取正时，表示正相关；当ρ(Z，Z′)取负时表示负相关；ρ(Z，Z′)取0时，表示不相关。相关系数越大，表示2组变量相关性越大，当ρ(Z，Z′)=1时，说明2组变量完全相关。

1.4 点云配准

设M、N分别为存在交集的2个相邻机载LiDAR点云片S1和S2上的点集，S1、S2重合区域C=S1∩S2。记M落入C区域的采样点所构成的集合为CM={Mi|Mi∈R3}，N中落入C区域的采样点所构成的集合为CN={Ni|Ni∈R3}。具体关系如图2所示。

图2 点云数据集关系图

将点集M固定，点集N作为浮动点集参与配准。然后在CM中按照1.2节的方法选择一个初始目标点作为待定点p，以该点为中心搜索范围l以内的邻近点，建立拟合曲面ρ。为了在CN上搜索点p的同名点，须估计出该同名点可能存在的范围，建立一个搜索区域，如图3所示。在搜索区域中选择一点q作为局部采样点集的中心点，以l为半径查找q的邻接点组成采样点集。同时，为了方便求取搜索区中采样点集的平面坐标落在拟合曲面的坐标范围内，须对目标点集和采样点集进行坐标中心化处理，即将点集中心点平面坐标设置为(0，0)，点集中其余各点依次用相对于该中心点的坐标表示。通过遍历规则在搜索区域中选择不同的采样点，对比不同采样点集中各点在拟合曲面上的高程和各点真实的高程之间的相关系数，选择相关系数最大的一组采样点集合作为目标点集的对应点集，从而寻找出同名点对。

图3 搜索区域内查找点示意图

点云配准的具体过程如下。

1)2片相邻点云S1和S2，将S1固定作为基准点云，S2作为浮动点云与S1配准。

2)在CM中选择采样点p作为基准点云的目标点，并以p点为中心、l为半径搜索p点的邻近点，组成目标点集。以p点的平面坐标为(0，0)更改目标点集的平面坐标。

3)将p点及邻近点按照式(1)生成拟合曲面ρ，计算点集中每一点到拟合曲面的偏移距离，若偏移距离大于设定阈值则重新选取采样点，重复以上步骤，直到找到合适的曲面。

4)在CN中预估p的同名点可能出现的范围，建立搜索区。

5)在搜索区选择采样点q，并以q为中心、l为半径查找q点的邻近点，组成搜索点集。以q点的平面坐标为(0，0)更改搜索点集的平面坐标。

6)设搜索点集中每一个点的高程为Zi，并计算每一个点在拟合曲面上对应的高程Z′i。

7)计算Zi和Z′i的集合Z、Z′的相关系数ρ(Z，Z′)。

8)在搜索区中选择其他采样点，循环步骤5)、步骤6)、步骤7)，直到搜索区中的采样点遍历完。

9)计算相关系数ρmax=max{ρi}，对应搜索点集的采样点qi即为p点的同名点。

10)循环步骤2)至步骤9)，在CM中选择若干目标点，并在CN中寻找对应的同名点，生成同名点对(pk，qk)。

11)遍历所有的同名点对(pk，qk)，根据同名点之间的坐标关系计算出旋转和平移矩阵，反复迭代计算，直到达到最佳配准效果，完成配准。

点云配准的具体实施流程如图4所示。

图4 点云配准总流程图

2 实验及结果

机载LiDAR数据覆盖范围较广，且主要反映地表及地表各种建筑物信息。为验证本文方法的正确性和有效性，选取了4组机载LiDAR点云数据进行配准实验。实验系统环境为内存8 GB、Windows10操作系统，数据处理的开发平台为Visual Studio 2012，编程语言为C#语言，点云的可视化软件为Cloud Compare 3D点云处理软件。

第1组实验为某地表A的相邻2条航带的机载LiDAR数据，地面分辨率10 cm。图5(a)为该地表配准前的原始点云数据(红色部分表示基准点云，白色部分表示待配准点云)，可以看出2片点云在配准前存在明显的位置偏差。图5(b)为传统ICP算法配准结果，配准精度RMS(root mean square)为0.204 69 m。图5(c)为本文方法的配准结果，配准精度RMS为0.218 24 m。

图5 地表A配准图

第2组实验为某地表B的相邻2条航带的机载LiDAR数据，地面分辨率10 cm。图6(a)为该地表配准前的原始点云数据(红色部分表示基准点云，白色部分表示待配准点云)。第2组数据与第1组数据均为略有起伏的地表，与第1组数据的区别在于第2组数据地表B中存在一栋规则建筑物。图6(b)为传统ICP算法配准结果，配准精度RMS为0.412 20 m。图6(c)为本文方法的配准结果，配准精度RMS为0.398 42 m。

第3组实验为某建筑物的相邻2条航带的机载LiDAR数据，地面分辨率10 cm。图7(a)为该建筑物配准前的原始点云数据(红色部分表示基准点云，白色部分表示待配准点云)。图7(b)为传统ICP算法配准结果，配准精度RMS为0.426 53 m。图7(c)为本文方法的配准结果，配准精度RMS为0.385 73 m。

第4组实验为某建筑物屋顶的相邻2条航带的机载LiDAR数据。图8(a)为该建筑物屋顶配准前的原始点云数据(红色部分表示基准点云，白色部分表示待配准点云)。图8(b)为传统ICP算法配准结果，配准精度RMS为0.493 54 m。图8(c)为本文方法的配准结果，配准精度RMS为0.437 11 m。

图6 地表B配准图

图7 建筑物配准图

图8 建筑物屋顶配准图

在上述几组实验中，对点云预处理后，使用了传统ICP算法与本文所述算法进行了对比实验，实验相关参数及结果见表1所示。对比以上4组实验的结果以及表1数据中的配准精度可以看出：就机载LiDAR数据而言，对于地表或含有少量建筑物的地表点云数据的配准方面，传统ICP算法和本文提出的算法区别不大，均能达到较好的配准效果；对于单个地面建筑物的配准而言，本文提出的算法通过多次迭代后的配准精度与传统ICP算法的配准精度相当，甚至略优于ICP算法的配准精度。

表1 不同算法的配准结果

机载LiDAR数据为大范围内的条带点云数据，点云密度较小。以上几组实验数据均为局部小范围的数据，其实验结果存在一定局限。为验证本文算法对大范围地表点云的配准效果，进行了2组实验(实验5、实验6)。实验5为某城区相邻2条航带的机载LiDAR点云数据，覆盖范围约为600 m×200 m。区域内多为规则建筑物屋顶、马路及景观树。图9(a)为该城区2条相邻航带的原始LiDAR点云数据(红色部分、白色部分分别表示不同的条带数据)。图9(c)为该城区点云数据利用传统ICP算法配准结果。图9(e)为该城区使用本文方法的配准结果。实验6为某地表相邻2条航带的机载LiDAR点云数据，覆盖范围约为300 m×400 m，区域内为坡度较缓的地表(地面有少量树林、一条马路及路灯若干)。图9(b)为该地表2条相邻航带的原始LiDAR点云数据(红色部分、白色部分分别表示不同的条带数据)。图9(d)为该地表数据利用传统ICP算法配准结果。图9(f)为该地表使用本文方法的配准结果。

图9 实验5、实验6机载LiDAR点云配准图

实验5、实验6的配准精度结果如表2所示，表中还列出了2组实验的点云配准变换参数。从图9和表2中的实验结果可以看出，对于大范围的机载LiDAR条带点云数据的配准，本文所提出的基于机载LiDAR局部点云拟合高程相关的配准算法与经典的ICP算法精度相差无几，2组实验都达到了较高的精度水平。而本文提出的方法，通过在基准点云中选取特征点，在待配准点云中使用本文提出的算法计算特征点对应的同名点，通过所计算的若干同名点对反复迭代完成2片点云的配准，在2片待配准点云初始位置不是很好的情况下也能完成点云的高精度配准。

表2 实验5、实验6不同方法配准参数及精度表

3 结束语

针对机载LiDAR点云的特点，本文提出了一种基于机载LiDAR点云高程相关的配准方法，通过将基准点云片内局部范围内的点拟合曲面作为特征曲面，在待配准点云中对比求取在拟合曲面的拟合高程，计算拟合高程与实际高程的相关系数，选择待配准点云中相关系数最大位置作为配准的同名曲面。通过建立机载LiDAR三维点云的局部小范围拟合曲面，将离散机载三维点云配准问题简化为实际高程与拟合高程的相关性对比问题。实验表明，本文算法与传统ICP算法相比,精度相当。但是，在实验过程中也发现，对于拟合曲面点的选取、拟合半径的设定还有值得深入研究的地方，对于不同地形、地物、不同密度的点云设定不同的拟合半径和搜索半径往往得出不一样的结论。针对不同地形、地物，不同密度的点云该如何选择拟合半径这一问题将是后期研究工作的重点。