点云配准是SLAM、计算机视觉、三维重建、多源三维数据融合等领域的关键技术，核心是求解待配准点云到目标点云的空间旋转平移矩阵。随着科学技术的发展，点云可通过激光雷达扫描仪、立体相机^[1]以及基于运动恢复结构^[2]（Structure From Motion，SFM）的无人机影像重建等方式获取。在实际应用中，传感器需要多次采集并进行点云配准才能获取目标或测区所需的精准完整点云。由于传感器采集数据对不同空间目标有不同视角、不同的传感方式，不同传感设备对同一传感目标，或同一目标的不同空间位置有不同的传感精度和传感能力，因此需要进行多源数据融合，比如激光点云与影像重建点云的配准工作，以获得更精准完整的三维空间信息。

在点云配准问题上，国内外学者已经做了许多相关工作。Besl^[3]等提出具有开创性的ICP算法，虽能提供收敛的精确结果，但应用条件苛刻，要求数据无噪声、完全重叠、初始位置良好^[4]。而实际的多源点云含有大量噪声、且部分重叠^[5]，为解决这一问题，业内学者提出了各种ICP变体算法。

针对含噪声和部分重叠点云的配准，Chetverikov^[6]提出了Trimmed ICP（TrICP）算法，利用最小截断二乘(Least Trimmed Squares, LTS)以点云重叠度作为比例来保留距离排序后距离较小的点对。戴静兰^[7]以点云的特征曲率提取特征点，通过特征点进行匹配提升ICP算法的效率。针对ICP对初始位置敏感的问题，钟莹^[8]使用主成分分析法（Principal Component Analysis, PCA）对点云进行初始配准，防止ICP陷入局部最优，并使用方向向量阈值剔除点对，提高了ICP算法的效率和精度。左超^[9]用信息论中的熵函数描述点云的空间分布状态，提出了ILSDE-ICP算法，为ICP精配准提供了良好的初值。在此基础上，陈杰^[10]提出结合最小空间分布熵与遗传算法的空间最优变换矩阵求解算法，有效提高了最小点云空间分布熵的求解效率。在激光与影像重建点云数据配准上，冯鸣^[11]等使用ICP算法对高陡崖的激光扫描点云与影像点云进行配准，获取了测区的完整点云。曹明兰^[12]以倾斜影像点云为主，激光点云为补充测量，使用ICP算法配准点云数据，构建了森林景区的三维模型。

除了ICP算法，也有学者提出了其他算法实现点云配准。Magnusson^[13]提出正态分布变换算法（The Normal Distributions Transform, NDT)，将三维目标点云表示成多个连续可微的概率密度函数，使用牛顿法优化得到待配准点云在目标点云上的最大得分以实现点云的配准。韩贤权^[14]提出了一种改进的粒子群优化法，在选取最佳配准点的基础上实现了散乱点云的精确配准。Aiger^[15]提出4PCS算法，用比例和距离约束以及RANSAC算法来查找两点云中的近似全等共面4点集，允许原始噪声、离群值，无需滤波去噪，无需初始配准，但需要点云间的重叠度信息。

以上方法能对普通点云进行配准并取得良好的配准效果，但是对于数据量大、非同源、含大量噪声且部分重叠的激光点云与影像重建点云，其稀疏程度、噪声程度等不同，非重叠区域的面积很大。对于上述基于特征的点云配准，无法提取出足够相同的特征；由于真实场景的点云尤其是影像重建点云噪声较多，提取的法向量误差也很大；对于3D-NDT和4PCS，对不同点云模型进行配准，需要多次实验以获取最优配准参数。因此，文中提出基于贡献因子的改进TrICP算法：CFB-ICP (Contribution Factor Based ICP)算法，使用改进体素降采样来得到降采样点云，再由贡献因子筛选出对配准贡献度更大的点对以加速配准过程，并根据距离曲线上若干点过原点的斜率（以下简称斜率）特性来自动计算并更新重叠度$\xi $，在保证配准结果正确的情况下，实现点云的高效、快速、全自动配准。

经典ICP算法的核心是查找最近邻点构成点对，通过匹配到的点对集进行误差函数最小化以求得旋转平移矩阵。实际应用中，待配准点云会包含大量的离群点。一方面由于点云的部分重叠，非重叠部分的点在参考点云中只能找到错误的匹配点对^[16]。另一方面，任何测量方式都会带来误差和噪声，这些点云与目标点云中的对应点存在偏差。

为克服上述问题，Chetverikov^[6]提出TrICP算法，使用LTS来剔除距离过大的点对，剔除的点对由点云重叠度$\xi $决定。对于重叠度，给出如下定义：设有待配准点集$P$，参考点集$Q$，${P_i}\{ i = 1,2,3, \cdots ,{N_P}\}$是$P$中的一点，${P_i}$在$Q$中的最近点是${Q_i}\{ i = 1,2,3, \cdots ,{N_Q}\}$，假设$P$中有一子集，子集中所有的点都能在参考点云$Q$中找到正确的对应点，设$A$是满足该条件的最大子集，即：

$\forall {A_i} \in \{ A\} ,\exists {Q_i} \in \;{\rm{\{ }}Q\} \;,\{ i = 0,1, \cdots ,{N_A}\}$，有$ \Vert {A}_{i}-{Q}_{i}\Vert < $$ \epsilon $，其中$ \epsilon $为距离阈值，则点云重叠度$\xi $为：

TrICP算法的步骤如下：

（1） 对$P$中的所有点${P_i}$查找其在$Q$中的最近邻点${Q_i}$，并计算其距离的平方$D_i^2$；

（2） 对所有$D_i^2$进行从小到大排序，保留前${N_{PO}}$个点对并计算它们的和为${S_{TS}}$，其中：

（3） 如果满足任意停止条件则停止算法，否则令$S_{TS}^\prime = {S_{TS}}$，继续下一步；

（4） 使用步骤（2）得到的${N_{PO}}$个点对代入目标函数并通过奇异值分解法^[3]（Singular Value Decomposition，SVD）求解旋转变换矩阵$\left( {R,T} \right)$：

（5） 用求得的$(R,T)$矩阵对$P$进行空间旋转平移变换：

转到步骤（1）。

其中，步骤（3）中满足以下任一情况即可以作为停止条件：

1) 迭代次数大于设定阈值${N_{iter}}$；

2) 点对距离均方误差(Mean Squared Error, MSE)值$e$小于给定阈值，其中：

3) 两次迭代的Trimmed MSE的绝对差值|e′−e|足够小，其中e′是上一次迭代的Trimmed MSE。

ICP算法是TrICP算法在重叠度等于1时的特殊情况^[6]。

对于重叠度，Chetverikov^[6]使用黄金搜索算法自动求得最佳值，但是这一过程需要多次配准以寻找最优的重叠度，需消耗巨大的计算资源。一般为手动设置，对于视觉上无法直接观察重叠度的点云模型，要多次手动设置重叠度才可得到最优结果。

为减少点云的噪声^[17]及数据量以提高配准精度与速度，在点云配准之前需要进行降采样。随机采样是常用的采样方法，Pomerleau^[18]使用固定数量点随机采样而不是固定比例，以此来保证不同规模的点云采样得到相同数量的点云，但是其结果并不均匀，采样结果的不均匀性和原点云一致^[13]。体素降采样是将点云划分至$n$个空间三维体素网格，公式(6)计算了每个体素的重心${P_j}\{ j = 1,2, \cdots ,n\}$，该体素内的所有点将用重心${P_j}$来表示，以达到均匀降采样目的：

式中：$\;n$为体素个数；$m$为第$j$个体素中的点云个数。

但该方法降采样出来的点云并不是原来的点，而是计算出来的重心，会增加点云配准的误差。文中提出改进体素降采样算法，利用公式(7)查找降采样后的每个点${P_j}$在原点云中的最近邻点${P_i}\{ i = 1,2,3, \cdots ,{N_P}\}$，并将其替换。这样采样出来的点依然属于原点云，并且避免了随机采样的不均匀。在此基础上再采用固定点云数量采样^[18]，使不同点云能采样出相同数量的点云以保证算法的速度以及算法间的可比性^{[16, 18]}。

在经典ICP与TrICP的算法中，当算法收敛到一定程度时，后续收敛变得相对很慢。以每次迭代重叠部分的MSE来绘制曲线图，如图1(c)所示，ICP算法与TrICP算法在一开始时收敛较快，后面部分处于缓慢收敛的状态。实际上，在利用最近点对求解的过程中，并非所有点对在配准上的贡献度都相同。当点云配准阶段还处于两点云距离较大时，参与求解的点对并不是实际理想的对应点，如图1(a)所示，而是最近距离点对，如图1(b)所示，并且还可能是噪声^[16]。经实验验证，此时距离相对大的点对对收敛到正确结果的贡献度高于距离小的点对，甚至距离小的点对是陷入局部最优的主要因素。

图  1  (a) 理想对应点对；(b) 最近距离对应点对；(c) 算法的收敛曲线；(d) 点对距离曲线；(e) 距离曲线上的点过原点的斜率；(f) 不同对应点对斜率变化曲线

Figure 1.  (a) Ideal corresponding points; (b) Nearest distance corresponding points; (c) Convergence curves of algorithms; (d) Distance curves of corresponding points; (e) Slopes between points on the distance curve and the original point; (f) Slope curves of different points

实验发现，点云配准过程中的点对距离${D_i}$的变化曲线随迭代次数增加趋于指数增长。使用上述采样方法使${N_P} = 2\;000$进行配准实验，在不同迭代次数下将点对距离${D_i}$由小到大排序后的曲线如图1(d)所示，可知随着迭代进行，重叠部分的点对的${D_i}$相对较小，非重叠部分的点与其最近点的${D_i}$远大于重叠部分的点对距离，在重叠到非重叠过渡段，距离${D_i}$开始快速增长。尤其对真实大场景的点云，重叠部分点对距离在0～2 m之间，而非重叠部分则是大于10 m甚至100 m。因此，重叠度$\xi $与距离曲线的快速增长区间有关。文中引入距离曲线上的若干点，对过原点斜率进行分析，如图1(e)所示，${k_i}$表示第$i$个点对的距离值过原点的斜率，当$i > 1\;750$，${D_i}$与${k_i}$迅速增大。对于一组确定的点云配准模型，${k_i}$总是随着迭代过程趋于稳定，如图1(f)所示。在点对距离较小的重叠部分，如${k_i { < 1\;750}}$集中在一个小范围，斜率曲线密集分布，而非重叠区域对应的${k_i { > 1\;750}}$明显变大，且斜率曲线分散。

基于上述分析，文中提出基于贡献因子$\alpha $的CFB-ICP算法，该算法的核心是，使用距离阈值$\delta $来得到贡献度大的点对（有效点对中距离更大的点对）来进行矩阵求解，同时根据距离曲线斜率值分布来自动快速更新重叠度$\xi $，以更新的$\xi $来参与下次迭代运算，并添加斜率变化约束，使算法在收敛到一定程度时恢复所有有效点对进行配准，以保证算法能收敛到一定精度，有效点对指重叠部分的点对。

改进算法对TrICP算法的明显提高主要有两方面：

（1）使用贡献度大的点对加速了配准的收敛过程；

（2）算法在迭代的同时自动计算重叠度$\xi $。

设初始重叠度$\xi = 1.0$，改进配准算法主要步骤如下：

（1）对$P$中的所有点${P_i}$查找其在$Q$中的最近邻点${Q_i}$，并计算其距离${D_i}\{ i = 1,2,3, \cdots ,{N_P}\}$；

（2）对所有${D_i}$进行从小到大排序，对由${D_i}$构成的距离曲线进行等间隔采样$j$个点，并求得采样点过原点的斜率${k_m},\;\;\{ m = 0,1, \cdots ,j\}$；

使用斜率变化约束以停止更新$\xi $：

式中：${\sigma _1}$初值为0；${k^\prime }_m$为上一次迭代的斜率，$m = $$ 0,1, \cdots,j$；若连续${n_{{\sigma _1}}}$次${k_m}$的变化率小于等于${R_c} = 0.01$，即${\sigma _1} > {n_{{\sigma _1}}}$，则停止更新$\xi $，${n_{{\sigma _1}}}$取经验值5。

公式(9)以斜率的绝对偏差来更新$\xi $，若连续${n_{{\sigma _2}}}$次${k_l}$的绝对偏差与整体的平均绝对偏差的比值大于$\beta $，即${\sigma _2} > {n_{{\sigma _2}}}$，则利用公式(10)和(11)更新$l$与$\xi $，其中${n_{{\sigma _2}}}$取经验值5，$\;\beta$取经验值2.5，$l$初始值为$l = j$；

公式(12)根据重叠度$\xi $计算${N_{PO}}$个点对的${D_i},\{ i = $$ 0,1, \cdots ,{N_{PO}}\}$之和记为${S_{TS}}$；

（3）如果满足任意停止条件则停止算法，否则令$S_{TS}^\prime = {S_{TS}}$；

（4）保留${N_c}$个点对，其点对距离${D_i}$满足阈值条件：

式中：${\delta _{max}} = {D_{{N_{PO}} }}$；

（5）更新阈值${\delta _{min}} = {D_d}$供下次迭代使用：

式中：${\delta _{min}}$初始值为0，$\alpha \in \left( {0,1} \right)$为贡献因子，经验值为0.5；

（6）使用步骤（4）得到的${N_c}$个点对代入目标函数，通过SVD求解最优旋转变换矩阵$\left( {R,T} \right)$：

（7）用求得的$(R,T)$矩阵对点云$P$进行空间旋转平移变换：

转到步骤（1）。

其中，步骤（3）中的停止条件与上述TrICP算法相同。为了点对数量不至于太少以保证结果的可靠性，增加了更新${\delta _{min}}$的约束，当${N_c} < {N_{PO}} \times 0.1$，即当前点对数小于有效点对数的十分之一时，更新${\delta _{min}} = {\delta _{min}} \times $$ 0.9$，直到${N_c} > {N_{PO}} \times 0.1$。

为检验改进算法性能，使用Pomerleau的模块化配准方法^[19]以及点云库PCL 1.10.0来实现ICP、TrICP、3D-NDT、4PCS和文中CFB-ICP算法，点云数据为斯坦福大学公开的0°与45°扫描Bunny点云组、TBunny (由0°Bunny裁剪得到)点云组以及Hippo点云组，在Bunny点云上添加噪声以检验算法鲁棒性，最后将算法用于滑坡体的激光与影像重建点云。实验平台为Windows 10+64位操作系统，Intel core i5-6500CPU，内存16 G。

CFB-ICP算法可直接对点云进行配准。其中$\alpha = 0.5$与$\;\beta = 2.5$由实验得出，重叠度精度与采样间隔$j$有关，文中取$j = 40$；为保证速度，需将点云降采样，但过少的点云不能得到正确的结果，经实验得出${N_P} = 2\;000$同时兼顾速度与结果的正确性；CFB-ICP能自动计算出点云的重叠度，Bunny为0.875，Hippo为0.725，TBunny为0.575，并将其应用于TrICP算法。点云初始状态及配准结果如图2~4所示，红色为目标点云，绿色为待配准点云，蓝色为配准结果点云。视觉上可以看出，对不同点云(Bunny与Hippo)以及不同重叠度点云(Bunny与TBunny)，CFB-ICP均取得较好的配准效果；对于ICP，配准后的点云并没有完全贴合；对于TrICP，使用由CFB-ICP得到的重叠度进行配准，Bunny与TBunny取得了良好的配准效果，而对于Hippo，得到了错误的局部最优结果，因为在初始Hippo的腹部位置的点云点对距离较小，这也证明上文所述，使用距离更近的点对进行矩阵求解反而可能使算法陷入局部最优；对于3D-NDT，三组点云也均得到较好效果；4PCS寻找共面4点集时为随机选点，导致每次结果不一致，若要得到相对更稳定的结果，需增大时间阈值，因此增加了时间消耗；3D-NDT与4PCS都对参数敏感^[20]，对不同的点云模型需要通过实验调整以获取最佳参数，相对来说使用不便且耗时。因此，CFB-ICP在点云的模型不同、重叠度更小的情况下均表现最优。

图  2  Bunny配准结果

Figure 2.  Registration results for Bunny

图  3  Hippo配准结果

Figure 3.  Registration results for Hippo

图  4  TBunny配准结果

Figure 4.  Registration results for TBunny

文中采用重叠部分点云的点对距离均方误差MSE以及运行时间作为评价指标，定量地对比算法结果。表1为各算法时间消耗以及MSE，算法时间取10次实验平均值；MSE均用同样的点对计算；由于4PCS算法每次结果不一致，其MSE也取10次实验的平均值。由表1可知，ICP与4PCS配准的MSE相对较大，且4PCS更耗时；TrICP对Hippo配准失败；3D-NDT也得到较小的MSE，但最耗时；TBunny的相对MSE大小结果与Bunny相似，均是CFB-ICP的MSE最小，虽然TrICP也成功配准且时间更少，但算法需要CFB-ICP提供重叠度信息。CFB-ICP配准的三种模型MSE最小，并且相比于TrICP算法，CFB-ICP对Bunny与Hippo的配准时间效率分别提升50%与64%，因此，CFB-ICP在效率和精度上均有优势。

为验证算法的抗噪声能力，同时在0°(${Q_0}$)与45°(${P_0}$)扫描的两个Bunny点云上添加高斯噪声。以对${P_0}$添加噪声为例，噪声点云的点数量为${N_G} = \gamma \times $$ {N_{{P_0}}}$，$\gamma \in \left[ {0,1} \right]$，即按原点云数量的0%～100%比例添加。从原点云中随机抽取上述${N_G}$个点，得到点云$G$，高斯噪声点由抽取点${G_i}$的xyz坐标分别加上一个高斯随机数得到，即：

${G_{iy}}$与${G_{iz}}$同理，$g$服从高斯分布$g\sim N\left( {u,{\sigma ^2}} \right)$，其中$u = 0,\;\;\sigma = 0.5$，均值$u$决定了噪声点云相对于原点云的整体偏移，设置为0；方差$\sigma $由具体使用的点云决定，该实验使用Bunny点云，设置$\sigma = 0.5$能很好地反映算法对大量噪声的鲁棒性。

将噪声点云$G$添加至原始待配准点云${P_0}$，得到最后用于实验的噪声点云，参考点云${Q_0}$同理：

此时，噪声点云总数为${N_{{P_0}}} = {N_{{P_0}}} + {N_G}$。

在采样点${N_P} = 2\;000$时，CFB-ICP算法与TrICP在原点云添加0%～60%比例的噪声下进行多次实验，均能得到正确的配准结果，其余算法效果较差；继续增加噪声至70%～100%，由于噪声的随机性导致点云的差异，出现部分配准失败，但若将采样点增至${N_P} = $$ 5\;000$，CFB-ICP算法与TrICP算法成功对点云配准，如图5为两点云同时添加100%噪声后的配准情况。表2为添加100%噪声后及${N_P} = 5\;000$时，各算法的时间消耗以及MSE，对比上文可知采样点增加后，所有算法配准时间都变长。TrICP与文中算法虽然都能对噪声点云进行配准，但是文中算法在时间上更有优势。所以可看出所提算法抗噪能力较强，且必要时可以增加采样点以牺牲时间来换取配准精度。

图  5  Bunny添加噪声后的配准结果

Figure 5.  Registration results after adding noise to Bunny

文中使用上述几种算法对“茂县624”滑坡体进行激光与影像重建点云的配准。该滑坡体垮塌体体积大约300万立方米，对基于多源数据融合的灾害监测研究具有代表性，影像重建点云由Pix4d对大疆如风8拍摄的现场影像进行重建获取，如图6(a)所示，该点云跨度最大方向为2 500 m，点数量超过200万点；激光扫描点云由Optech Polaris扫描获取，如图6(b)所示。因为激光扫描点云的精度高于影像重建点云，所以激光点云为配准输入的参考点云，影像重建点云为待配准源点云^[12]。配准结果如图7所示，文中CFB-ICP算法与TrICP算法均配准成功，但文中算法无需手动设置重叠度，且由表3可知，文中算法MSE更低且时间效率提升了67%。

图  6  茂县滑坡体点云。(a) 影像重建点云；(b) 激光扫描点云

Figure 6.  Point cloud of Maoxian landslide. (a) Image reconstruction point cloud; (b) Laser scanning point cloud

图  7  茂县滑坡体点云配准结果

Figure 7.  Registration results of Maoxian landslide point cloud

文中提出一种适应于多种点云模型并自动计算重叠度的自动配准算法，无需设置重叠度与复杂的参数即可实现点云配准。首先使用改进降采样方法，通过体素降采样点的最近点替换并随机采样固定数量点，解决了体素降采样引入误差以及随机采样不均匀的问题。针对多源点云、噪声、部分重叠以及配准参数确定问题，提出基于ICP的CFB-ICP算法，通过$\delta $阈值分离出对配准贡献度更大的点对，并采用距离曲线上的点过原点的斜率特性来自动计算重叠度，实现普通点云、激光与影像重建点云以及含大量噪声点云的高效配准。经对公开数据以及“茂县624”滑坡现场点云的实验验证，相比于ICP、TrICP、3D-NDT、4PCS算法，文中算法精度更高，时间更少，抗噪声能力更强且可以自动计算重叠度，无需根据配准点云模型变化而重新调节参数；即使存在大量噪声，也可通过增加采样点牺牲时间来换取精度。ICP对初始位置依赖很强，文中针对实际工程数据特点对其算法进行了改进，未来需探究改进算法对初始位置运算的鲁棒性，并深入探究非同源点云配准后数据误差分析以及去噪等问题。