受扫描仪自身分辨率的限制及系统噪声的影响，待配准的多个测站点云间并不存在严格意义上的同名点对，利用同名点方法难以准确地获取测站间的匹配关系^[15]。利用建筑物场景中普遍存在大量平面特征的特点，将海量点云的配准转换为少量同名平面的配准，是简化计算、提高配准准确性的有利途径。

本文采用一种高效RANSAC方法^[16]从点云中提取平面，然后利用MCMD_Z方法^[17]对提取的平面进行粗差剔除，拟合高精度平面基元。对属于同一个平面的点集进行主成分分析，将第一和第二主成分方向分别作为长和宽的方向，构建点集平面投影的外接矩形作为平面基元，矩形的几何中心即为平面基元的几何中心。计算平面基元的几何中心c_p、单位法向量n_p。采用何文峰的方法^[15]计算点云在平面上的投影面积S_p。按照不同视角中相同平面法向量方向一致的原则，利用Pathak的方法^[18]对平面基元的法向量进行调整，使场景中平面的法向量统一朝向外侧。由于平面基元的面积和二面角在旋转平移变换中具有不变性，故本文采用面积和二面角作为相似性测度。

建筑物场景通常包含大量对称或相似的平面。如果直接在目标和参考平面集中搜索同名平面，计算量较大且部分相似平面难以区分。针对此，本文通过对平面基元分组以简化场景。

平面基元分组即依据平行平面法向量方向一致原则，将场景划分为内部存在平行关系的平面基元组。借助平面基元组求解粗略旋转参数，在粗略旋转参数约束下建立平面匹配，减少投票方法的穷举次数，简化计算。

在建筑物场景点云配准中，通常难以事先确定平面基元组的数目。因此，本文采用层次聚类方法，根据平面的方向进行分组。可认为法向量夹角小于夹角阈值t_θ的平面方向相同，可将其归为同一平面基元组。假设有一个球心位于原点，半径为单位长度的球，则可以将场景中所有平面(图 1(a))的单位法向量看作一个分布在该球表面的点集(图 1(b))，利用平均距离对点集进行层次聚类(图 1(c))，对聚类后的类别求类别中心；选取距离类别中心最近的球面样本点对应的法向量为该类别的组法向量；根据分类树结果选取合适的夹角阈值t_θ，并得到距离阈值t_d=r·t_θ，其中r=1，为球面半径；按照距离阈值选取合理的分类作为分组结果(图 1(d))。

图  1  利用层次聚类对平面基元分组示意图

Figure 1.  Schematic Diagram of Planar Primitives Grouping Using Hierarchical Cluster

计算从目标点云到参考点云的变换参数，需在两个集合中找到3对或3对以上不平行的同名平面。根据同名平面基元组的法向量计算旋转变换参数，利用法向量和原点-平面距离获取平移矩阵。若平面基元正确匹配，则同名平面的法向量方向一致，此时两视角中方向一致的平面对数目最多。在此约束下，对相同面积的平面基元建立匹配，并利用单位四元数法^[19]计算变换参数。

同一坐标系内的各个平面间的二面角不随系统的旋转平移变换而改变，二面角的这一特性决定了它可以作为用于计算旋转变换参数的几何不变量。由于所有平面基元都已按照是否平行分为若干组，各平面基元组之间均互不平行，则3个平面基元组A、B、C的组法向量n_A、n_B、n_C之间存在3个夹角θ_AB、θ_BC、θ_AC，这3个夹角分别对应3个平面基元组之间的姿态关系。本文采用3个平面基元组间的法向量夹角θ_AB、θ_BC、θ_AC和平面基元的面积作为不变量，结合RANSAC方法建立同名平面匹配，计算从目标数据集P到参考数据集Q的旋转矩阵R和平移矩阵t，步骤如下。

(1) 将3个平面基元组的组法向量n₁、n₂、n₃作为确定旋转变换模型的最小观测值集合，计算参考平面集合中任意两个平面基元组的组法向量夹角，夹角关系保存为LIST_Q。

(2) 随机抽取目标平面集合中的3个平面基元组G_Pⁱ(i=1，2，3)，计算组间代表法向夹角，得到3个角度值θ₁₂^P，θ₂₃^P，θ₁₃^P。在LIST_Q中搜索与这3个夹角差异小于t_θ的对应角，并结合平面基元组内基元的面积值S₁₃ⁱ，判断G_Pⁱ(i=1，2，3)和参考平面集合中的3个平面基元组G_Qⁱ(i=1，2，3)中是否存在同名平面(图 2)。如果平面基元组夹角存在对应关系，且对应平面基元组内存在面积差别小于阈值t_s的平面基元(图 3)，则认为可能存在同名平面。若不存在夹角和面积对应关系，则重新抽取新的平面基元组。

图  2  与目标平面基元组对应的参考平面基元组搜索过程

Figure 2.  Searching the Reference Primtive Groups Corresponding to Object Primitive Groups

图  3  目标平面与参考平面相同的情况

Figure 3.  Same Plane in Object and Reference Plane Sets

(3) 对可能存在同名关系的情况，根据平面基元组对应的3对法向量n_P1、n_P2、n_P3和n_Q1、n_Q2、n_Q3计算旋转变换R_i。使用该旋转变换R_i对目标平面集合中的所有平面进行旋转，计算旋转后的目标平面集合与参考平面集合中平面基元组的法向量夹角，统计法向量方向一致的平面基元组对的数目n_c。

(4) 重复步骤(2)~(3)K次，抽样次数K＞log(1－P_e)/log(1－P_o)。其中P_e为期望随机选择的3个平面基元组在重叠区域的概率，P_o为估计的两视角数据间的重叠区域的比例。

(5) 选取具有最大平行平面基元组对数目n_cmax的G_Pi^j和G_Qi^j(i=1，2，3)作为初步最佳映射。由于具有最大n_cmax的平面基元组对应可能不止一种情况，因此j≥1，该组合对应的R_m^j即为初步最佳组旋转矩阵。按照R_m^j对G_Pi^j(i=1，2，3；j≥1)作旋转变换得到变换后的平面基元组G_Pi^m_j。

(6) 在对应组G_Pi^mj和G_Qi^j中找到面积相同(面积差别小于阈值t_s且最接近)的平面基元Plane_Pi^j和Plane_Qi^j(i=1，2，3，j≥1)建立匹配关系。检查Plane_Pi^j和Plane_Qi^j对应的3对几何中心组成的3个矢量$\overrightarrow{P{{O}^{j}}_{i}}$(起点为Plane_Pi^j的中心，终点为Plane_Qi^j的中心)是否存在方向和长度一致关系(图 4)，以免场景中存在具有相同面积且平行的平面出现误匹配。若存在一致，计算对应的四元数q=q₀q₁q₂q₃^T和平移向量t。

图  4  平面基元几何中心矢量约束

Figure 4.  Constraint with Vectors from Reference to Target Planar Primitive Centers

(7) 计算所有G_Pi^m_j和G_Qi^j(i=1，2，3； j≥1)的组合对应的q和t，对计算结果中具有近似相等q和t的对应平面重新建立平面映射关系。按照最小二乘原则对匹配的平面特征计算最佳旋转矩阵R和平移矩阵，此变换参数也是目标点云到参考点云的变换参数。

为验证算法自动配准的精度，在点云的公共区域选取同名平面，利用三个平面求交点，将交点作为控制点，计算控制点转换后的精度^[20]：

式中，i=1，2，3，…，n，n为特征点对的数目；X_i、Y_i、Z_i为控制点在标点云中控制点转换到参考坐标系中的坐标值，ΔX_i=X_i－X′_i，ΔY_i=Y_i－Y′_i，ΔZ_i=Z_i－Z′_i。

本文采用Reigl VZ-400地面三维激光扫描仪，工作时，扫描仪沿水平面顺时针方向匀速旋转，激光束竖向逐列扫描被测物体表面，并自动解析记录激光的回波信号。本文进行了多组实验，现举两例予以说明。

实验一 扫描场景为一片建筑群，点云区域为200×100 m。扫描内容主要包含建筑物立面、行人、车辆、道路和植被。两测站重叠区域约70%。图 5为未经配准的两测站点云，图 5中黑色为目标点云P，点数量约2.43×10⁶个，灰色为参考点云Q，点数量约3.61×10⁶个。忽略面积过小或所含点数过少的平面，从目标点云P中提取平面基元32个，从参考点云Q中提取平面基元43个。按照0.1°夹角阈值，P、Q分别被分为9组和15组(图 6)，图 6中竖线为分组阈值。从P和Q中寻找7对存在公共平面的平面基元组，根据面积相等原则，从公共平面基元组中选出7对面积差别最小的平面基元建立匹配，计算变换参数。图 7为配准后的两测站点云。选取7个控制点检查，变换精度σ_P=0.015 m。

图  5  未经配准的两站点云

Figure 5.  Unregistered Point Clouds

图  6  P与Q中平面基元分组结果

Figure 6.  Result of Planar Primitives Grouping in P and Q

图  7  配准后的两站点云

Figure 7.  Registered Point Clouds

实验二 点云区域为500×300 m。扫描内容主要包含建筑物墙面、道路和植被。两测站重叠区域约40%。图 8为未经配准的两视角点云，黑色为目标点云P，点数量约1.40×10⁶个，灰色为参考点云Q，点数量约8.37×10⁵个。忽略面积过小或所含点数过少的平面，从目标点云P中提取平面基元26个，从参考点云Q中提取平面基元38个。按照0.5°夹角阈值分别被分为7组和10组(图 9)。从中选出5对匹配基元，计算得到变换参数。图 10为配准后的两测站点云。选取5个控制点检查，变换精度σ_P=0.062 m。

图  8  未经配准的两站点云

Figure 8.  Unregistered Point Clouds

图  9  中平面基元分组结果

Figure 9.  Result of Planar Primitives Grouping

图  10  配准后的两站点云

Figure 10.  Registered Point Clouds

实验一的变换精度高于实验二，主要是由于匹配的同名平面基元越多，求解的变换参数越精确。此外，当场景较大时，距离较远的建筑物点云较为稀疏，该区域噪声所占比例较大，造成平面基元的拟合误差较大，利用这些平面基元的配准结果也出现相应的偏差。当扫描场景较小时，一般只针对目标区域精细扫描，噪声所占比例较少，配准结果也较为理想。在实验二中，受场景重叠度限制，可用于匹配的同名平面基元数量较少。为了满足三对以上不平行平面的配准要求，选用了较远区域的建筑物立面。由于实验二中扫描仪设置的采样间隔较大，平面点云获取完整度降低，使同名平面基元的面积偏差增大，导致可匹配平面较少。且路面存在起伏，地表平面拟合误差较大，提取的平面基元包含噪声较多。多方面因素造成了实验二的配准精度降低。实际应用中，为了获取建筑物完整表面，会从多个角度对建筑物整体进行扫描，多个测站形成一个闭合环，通过分配闭合差的方法进一步减小相邻两站配准带来的累积误差^[21]。

与基于ICP的配准方法相比，本文方法只需要两次扫描获取的数据具有足够的重叠度，就可以求出用于粗配准的变换参数，更加方便灵活，无需标靶辅助和人工选取同名特征；相对于直接从点云中提取同名点的配准方法，本文方法在扫描线较为稀疏的点云配准时具有更好的鲁棒性；与其他利用平面特征的配准方法相比，本文方法采用平面分组的策略将场景中的大量平面简化为几个方向相同的平面集合，有利于降低搜索同名平面的复杂度。本文方法的配准精度对平面基元的参数精度较为依赖，要求建筑物场景点云相邻两站的重叠区域内至少包含3对或3对以上互不平行的同名平面，且平面的整体形状获取较为完整。

本文提出了一种基于平面基元组的点云自动配准方法。该方法将平面基元按法向量方向分组，采用RANSAC方法匹配平面基元，实现建筑物点云的自动配准。实验结果表明，该方法具有自动化程度高，无需人工选取初值的优点，可用于建筑物场景点云粗配准。