1 引 言

地面激光雷达扫描(terrestrial laser scanning)是一种近些年来飞速发展的主动式空间数据获取技术，可用于生成数字高程模型、数字线划图、三维城市模型等多种数字空间产品，在工程测量、数字城市建设、文物保护等领域具有重要的应用。为了完成对被扫描物体或区域的完整覆盖，往往需要在不同位置架设多个扫描站点。而每一站点的激光点云均处于彼此独立的三维坐标系下，无法直接用于后期处理，因此需要通过将不同站点的点云整合到统一的全局坐标系，即三维点云的拼接。点云拼接是点云处理流程中的重要步骤，点云拼接的速度也直接影响点云数据处理的整体效率。

为了更好地解决点云拼接问题，学术界已提出了许多种点云拼接方法。从用于计算转换矩阵的几何特征模型分析，现有的拼接方法可大致分为基于点、基于线、基于面、基于影像4类。

文献^[1]借鉴经典图像匹配SURF^[2]的方法计算提取点云中的特征角点及其法向量，从而形成每一点云独有的点特征签名(point feature histogram signature)。通过对相邻站数据的特征签名进行RANSAC(random sample consensus^[3])统计匹配出两个坐标系的同名点对，并进一步使用奇异值分解方法解算出6个未知参数。文献^{[4, 5]}同样使用类似的点特征方法进行拼接。

文献^[6]通过人工选定两个点云中的同名直线段进行拼接。文献^[4]使用自动匹配同名直线段拼接点云与影像。文献^[7]提出了融合面和曲线的最小二乘拼接方法。

文献^[8]将传统ICP中最小化对应点总距离的方式变化为最小化点到对应面的总距离，可达到比传统ICP更好的收敛率。文献^[9]将点到面ICP中非线性最小二乘解算的过程简化为线性最小二乘解算，在基本保证精度的情况下提高了算法效率。文献^[10]将点云中的特定特征拟合为圆柱、圆锥、圆球等常见几何形状，并利用同名物体的几何参数变化计算坐标系的6参数变换矩阵。文献^{[11, 12]}中使用的几何特征为点云中的平面特征。

文献^[13]提出了一种基于共线方程的改进丹麦法选权迭代法，可实现任意角度光学影像和点云的精确配准。文献^[14]使用与激光扫描仪同步拍摄的光学影像作为桥梁，借助计算机图形领域相邻影像匹配计算相对转换矩阵的方式计算相邻两站的拼接参数。文献^[15]也使用了类似的拼接方式，但使用的是激光点云反射强度所生成的灰度影像。

以上现有方法大多需要基准站和拼接站点云具有较多的同名特征，而全自动化的特征匹配对重叠区域的面积和特征类型具有严格要求，在实际作业中难以满足。目前商业化点云拼接软件如RiScan、Cyclone^{[16, 17]}，所提供的拼接方法主要是通过同名点选取，即根据人工观察判断出每两个相邻站点云的同名点，并且需满足至少4对和不可共面这两个条件，过程十分漫长。本文针对现有拼接方法同名点选取效率低下的固有问题，提出一种融合点云语义特征与测站中心GPS位置的地面激光点云拼接算法。该方法对于角点特征不明显的点云数据均有良好的适应性，可大幅降低同名面特征的标定时间，有效提高现有拼接过程的效率。

激光扫描可快速获取空间场景的海量信息，但由于激光点其实是对现实世界的采样，无法控制具体某个激光点落在某个精确的空间位置(如墙角)上。因此，即使对于两站相邻激光测站中的同名物体，也无法找到真正的“同名点”。例如在点间距均为5 cm的情况下，某一墙角的坐标只能假定为其范围在2.5 cm内的某一激光点坐标。此外，人工选点的操作方式也很难保证选择的是最佳点。因此，基于同名点计算出的转换矩阵难免带有较大误差，只能通过选择多对同名点以进行最小二乘平差。然而，目前主流激光扫描仪均已配备高精度的GPS设备，提供的扫描中心精确三维坐标可以作为计算拼接转换矩阵的一项特征。此外，激光扫描对于现实场景的采样率非常高，足以记录并恢复场景中的各种面。人类世界中包含地面、墙等大量平面特征。通过同一面上多个激光点拟合出的平面参数相对精确。使用同名“面特征”代替同名“点特征”，并结合测站中心的GPS位置，可以显著提高拼接算法的精度及效率。文献^{[18, 19]}提出了一种利用语义知识自动化识别地面激光点云中建筑物特征并进行三维重建的方法。该方法同样可以应用在点云拼接中以提高算法的自动化程度。通过对点云分割片的语义化解读，可以只保留对拼接有用的地面和墙体特征，从而将原始点云缩减为专门用于拼接的“特征子集”点云。对两个“特征子集”点云进行同名面片标定的效率，将远远高于从所有点中进行解读及标定的效率。

本文提出的点云拼接方法分为特征识别、初始拼接、精确拼接3个阶段(见图 1)。首先，通过分割及语义识别判断并提取出场景中的地面和墙体，并通过人工指定的方式标定两站间的同名面片；其次，结合测站中心的GPS位置估算出相邻测站点云的初始平移向量及旋转矩阵；最后，使用点到面迭代最近点(iterative closest point，ICP)^{[9, 10]}进一步修正平移及旋转参数以最小化两站数据之间的误差，从而达到站点的精确拼接。

图 1 初始姿态计算流程图 Fig. 1 Data flow of initial pose estimation

图选项

地面激光点云特征识别目的是从点云中自动识别提取出适合作为同名靶标平面的墙体及地面部分激光点簇，分为点云分割与语义特征识别两个步骤。首先，使用平面拓展法^[20](见图 2)将点云按照空间共性将点云进行平面分组，已形成可供进一步信息提取的平面分割。随后，将点云所获得的场景特征按照类型分类为：地面、建筑物、植被、临时物体(如车辆、行人)、杆状物体(如交通杆、电线杆、路灯)。以上特征中，植被、临时物体、杆状物体均曲面较多且面积较小，不适合作为计算初始姿态的平面特征。地面是地面激光点云中的主要部分，建筑物的墙面也是点云数据中的常见特征。地面和墙面均为点云中的显著平面特征，因此主要从此两类面片中选择用于计算初始姿态的同名平面。文献^[11]中给出了一种基于语义的特征识别方法。首先，通过知识总结出特征的语义化约束。随后，计算出每一分割片可用于语义推理的特征描述符，即大小、位置、方向、拓扑关系。通过对特征描述符及特征约束的概率化对比，即可判断出每一分割片所对应的特征类型。图 2给出了一幅点云数据的分割效果及自动提取的墙体特征。

图 2 点云的分割与特征识别 Fig. 2 Segmentation and feature recognition from a laser point cloud

图选项

点云经过特征识别处理后被简化为只包含适合于作为同名特征的点云分割片。本方法首先借助人工标定一到两对同名面特征，并结合通过RTK测量的基准站和拼接站中心GPS位置，计算出两个测站的初始转换矩阵，随后使用点到面最小距离的ICP方法进一步计算出精确拼接矩阵。

(1) 旋转。设基准站一靶标分割片的平面法向量为[P_A1 P_B1 P_C1]，其在拼接站同名靶标分割片的平面法向量为[P_A2 P_B2 P_C2]。拼接站对基准站的旋转方式为将拼接站坐标系绕轴L旋转α度，其中L为两平面法向的叉乘向量，α满足

(2) 平移。设旋转后基准站内两靶标分割片所在的平面的法向量为[A₁ B₁ C₁]、[A₂ B₂ C₂]，则两平面方程为

计算拼接站中心平移所依据的假设为：拼接站中心在基准站的坐标(X，Y，Z)距离基准站中心(0，0，0)的距离约等于GPS位置的距离R；在基准站坐标系下(X，Y，Z)距第一个基准靶标平面的距离，等于在拼接站坐标系中心距对应同名标靶平面的距离d₁；在基准站坐标系下(X，Y，Z)距第二个基准靶标平面的距离，等于在拼接站坐标系中心距对应同名标靶平面的距离d₂。

将以上假设转化为下列三元二次方程组

点击浏览公式

点击浏览公式

点击浏览公式

为计算方便，引入临时变量[a b c]、[A B C]、[a′ b′ c′]

平面(2)与平面(3)的交线方程为

点击浏览公式

解方程组可得到3种情况：

当b′²-4a′c′＞0时，方程组有两个解

用两个解分别计算两站同名平面的距离，距离较近的为正确值，距离较远的为镜像值。

当b′²-4a′c′=0时

当b′²-4a′c′＜0时，方程无解，测站中心GPS位置测量错误时会发生该状况，此时(X，Y，Z)应满足以下假设条件：①在线(4)上；②在过圆心0，0，0、垂直且与线(4)相交的直线上。

求解步骤为：

求过圆心且垂直于交线的平面

点击浏览公式

求平面(5)与线(4)的交点，即拼接站中心在基准站坐标系下的坐标

经过初始拼接后的拼接站坐标系与基准站坐标系大致吻合，但由于受到平面曲率、平面拟合精度、RTK精度的影响，两站数据之间仍存在一定误差，须经过迭代最近点ICP算法以进一步精确拼接。本文采用的点到面ICP^[6]是经典ICP的一种变种方法，其递归收敛的最小距离为拼接站点到基准站临近三角格网的三维距离。

本试验两组数据所采用的激光扫描仪为国产第一代地面激光扫描仪Stonex X300^[21]。该扫描仪具有300 m测距范围，360°×90°可视范围，扫描频率最高达40 kHz。除了配备激光扫描功能外，X300还配备有两个1070万像素的高清摄像头(本次试验中未使用相片)，并可外置连接GPS/RTK设备。扫描测站的中心坐标可通过置于扫描仪上的GPS坐标减去连接杆高度获得。算法由C++语言开发实现，运行平台为Windows 7 32位操作系统，机器配置为Intel i5 750 (2.67 GHz)，内存3 GB。

本组测试数据的作业区域为河南省南乐县王洪店村居民区。为了提高点云渲染及算法运行效率，原始激光点云数据均已通过k-d树抽稀处理至平均点间距5 cm。抽稀后基准站点个数为506 083，拼接站点个数505 436。本测区由于处于村庄居民区内，平面特征较为明显。本组测试数据初始拼接所选用的同名面为一对同名地面和一对同名墙面。两站重叠区域距离大约10 m，在此重叠区域内选取同名地面选取较为容易。但由于村庄建筑结构较为单一，同名墙面的选取往往需要视觉对比周围地物(如电线杆、大门形状)以作出正确判断。此外，此测区中狭长街道较多、路口较少，因此墙体特征一般只能选择到平行于测站中心连线方向的墙体，而很少能够找到垂直于测站中心连线方向的墙体。表 1列出了使用本文方法以及Riscan软件拼接模块(tie point selection方法)进行数据拼接的拼接速度以及中误差对比。图 3(a)给出了本组测试数据的拼接结果，由于数据场景的复杂性，在使用Riscan的拼接过程中我们很难寻找到重叠区域的对应物体，因此同名点对的选择过程较长，且无法确保理想的同名点对选择。

图 3 基于GPS位置和同名平面特征的拼接 Fig. 3 Registration from GPS positions and planar features

图选项

本组测试数据的作业区域为广东省广州市某写字楼区域。原始激光点云数据同样已通过k-d树抽稀处理至平均点间距3 cm。抽稀后基准站点个数为228 590，拼接站点个数276 935。此测区由于建筑物结构较复杂且地物特征多样，因此同名平面特征的识别较5.1节容易，甚至某些大型货车的侧面也可作为同名面特征。由于每站扫描之前扫描仪均一经过调平，因此RTK获得测站中心GPS坐标之后，只需选择一对不处于水平方向的同名面特征即可。本组测试数据初始拼接所选用的同名平面为某一建筑物的墙面。表 2列出了使用本文方法以及Riscan软件拼接模块(tie point selection方法)进行数据拼接的拼接速度以及中误差对比。图 3(b)给出了本组测试数据的拼接结果。本组测试数据场景较5.1节的数据简单，因此选择面状特征和点状特征均比较容易。与5.1节试验相比，本方法对于拼接速度与精度的提升并不十分明显。

以上测试数据的结果表明，首先，最终拼接质量与选取一对或者两对同名平面没有太大关系。只要初始拼接能够将数据大致逼近至ICP算法所能处理的范围以内，则即使初始拼接步骤之后存在一定误差也可以完成拼接。以本两组数据为例，只要选取的同名面正确且两组同名面不平行，则均可完成初始拼接并且误差足够使点到ICP达到正确收敛。但对于山区、农田等平面特征稀缺的情况，可以预见到本方法所采取的语义特征识别将难以开展。其次，拼接中误差与物体表面平滑程度有较大关系，这是因为在精确拼接阶段算法收敛的约束为点到面的总距离最小，而凹凸不平的区域所拟合的平面距点本身就具有较大距离。再次，拼接效率主要取决于人工解读场景的速度，而算法自动运行的时间一般均在数秒之内。在操作员熟练度相当的情况下，同名面判读的效率主要取决于场景是否易于视觉判读。如南乐数据由于对称及相似结构很多，造成一定程度上的人工判读困难，而广州数据由于建筑物外形独特，且路面有汽车、装饰等物体，操作员选定同名面片或者同名点均相对容易。如能将使用数码相片对点云着色，相信可大幅提高点云特征的人工判读速度。

本文给出了一种基于同名语义特征及测站GPS位置的拼接算法，可显著提高含有规则平面结构地面激光点云的拼接效率。与基于点特征的拼接方法相比，本文通过语义知识自动提取出地面特征和墙体特征作为面拼接所需的同名平面，从而克服了传统方法中人工匹配同名特征效率低下甚至无法找到同名点的实际难题，以及自动匹配同名特征对数据重叠区域要求较高的固有弊端。通过对两组地面激光点云数据拼接的试验，表明本方法具有良好的拼接精度以及显著的效率提升。为了进一步提高拼接过程的速度，下一步将在以下几个方面进行深入研究：①丰富语义知识库，使本方法可自动识别水塔、烟囱等曲面特征，从而利用此类同名特征进行自动化拼接；②将指北针、倾斜补偿器两类芯片植入X300扫描仪，并通过算法实现基于GPS、指北针数据和倾角数据的全自动化点云数据拼接。