<img src="https://www.zjujournals.com/eng/images/1008-973X/images/logo.png" class="img-responsive">

图 1 正态参数求解的方法

Fig.1 Methods for solving normal distribution parameters

令$ {\varTheta _X} = \{ {{\boldsymbol{\mu}} _i},\boldsymbol{\varSigma}_i\} _{i = 1}^{{K_X}} $和$ {\varTheta _Y} = \{ {{\boldsymbol{\mu }}_j},{{\boldsymbol{\varSigma}} _j}\} _{j = 1}^{{K_Y}} $分别表示源点云和目标点云各自空间所对应的正态分布参数集，其中$ {K_X} $、$ {K_Y} $表示源点云和目标点云具体包含的正态分布个数. 规定$ {{\boldsymbol{R}}^0}、{{\boldsymbol{t}}^0}$表示当前迭代过程的旋转矩阵和平移矩阵. 以Stoyanov等^[23]的研究内容为基础，对式（1）进行改进优化，以兼顾配准的精度和效率为目标，提出隶属于“分布到分布”配准模式的新的目标函数：

$ \left. \begin{split} \underset{\boldsymbol{R}, {\boldsymbol{t}}}{\arg \min } &\;{\sum}_{i=1}^{K_X} w_{i j} \left(\boldsymbol{\mu}_{i j}^{\mathrm{T}}\left(\boldsymbol{R}^0 {\boldsymbol{\varSigma}}_i\left(\boldsymbol{R}^0\right)^{\mathrm{T}}+{\boldsymbol{\varSigma}}_j\right)^{-1} \boldsymbol{\mu}_{i j}\right) ; \\&\text { s.t. } \boldsymbol{R}^{\mathrm{T}} \boldsymbol{R}={{\boldsymbol{I}}}_3, \operatorname{det}\;(\boldsymbol{R})=1 .\end{split}\right\} $

(2)

式中：$ ({\boldsymbol{R}},{\text{ }}{\boldsymbol{t}}) $表示待优化的变换矩阵，$ N({{\boldsymbol{\mu}} _j},{{\boldsymbol{\varSigma}} _j}) $表示源点云所处空间正态分布$ N({{\boldsymbol{\mu}} _i},{{\boldsymbol{\varSigma}} _i}) $经当前刚性变换$ ({{\boldsymbol{R}}^0},{{\boldsymbol{t}}^0}) $后的最近邻匹配结果，在此基础上可以定义$ {{\boldsymbol{\mu}} _{ij}} = {\boldsymbol{R}}{{\boldsymbol{\mu}} _i}+{\boldsymbol{t}} - {{\boldsymbol{\mu}} _j} $；$ {w_{ij}} $表示权重系数，能够指明$ N({{\boldsymbol{\mu}} _i},{{\boldsymbol{\varSigma}} _i}) $和$ N({{\boldsymbol{\mu}} _j},{{\boldsymbol{\varSigma}} _j}) $之间的相似度，使算法能够识别潜在的非重叠区域.

为了优化求解式（2），采用如图2所示的配准策略. 图中，初始点云当中的蓝色部分表示源点云，红色部分表示目标点云. 整个算法框架共分为3层，其中K-means聚类层并不参与算法整体的大循环，但会通过内部迭代来生成稳定的聚簇结果和对应的正态分布参数集；KL散度计算层和李代数求解器层则作为整体共同参与算法的整个迭代过程，直至达到大循环次数阈值上限或相邻2次配准结果没有明显区别，再输出最终的刚性变换矩阵.

图 2

图 2 基于正态分布相似性的双视角点云配准方法示意图

Fig.2 Diagram of pair-wise point cloud registration method based on normal distribution similarity

2.2. 任务求解

由式（2）可知，在利用李代数求解器优化配准结果前，需要2类关键信息：第1类信息是正态分布参数集合$ {\varTheta _X} $和$ {\varTheta _Y} $，一经K-means聚类层获取后无须频繁求解；第2类信息是最近邻匹配正态分布的相似性权重$ {w_{ij}} $，须利用KL散度计算层进行计算，且在算法整体大循环的每一次迭代中须重新计算.

给定初始刚性变换以及2帧点云各自的初始聚簇中心，规定h表示K-means聚类层的内部迭代次数，H表示KL散度计算层和李代数求解器层的整体迭代次数. 基于正态分布相似性的双视角点云配准方法步骤如下.

2.2.1. K-means聚类层

利用K-means聚类算法为源点云和目标点云分别生成若干正态聚簇，并利用其中的数据计算对应的均值和协方差矩阵.

1）为每个源数据$ {{\boldsymbol{x}}_i} \in {{X}} $搜索其在源点云当中的最近邻正态分布聚簇：

$ {c_h}(i) = \mathop {\arg \min }\limits_{k \in \{ 1,2, \cdots ,{K_X}\} } \left\| {{{\boldsymbol{x}}_i} - {\boldsymbol{\mu}} _k^{h - 1}} \right\|_2^2. $

(3)

2）利用新一轮的聚簇结果，为每个正态分布更新其均值信息：

$ {\boldsymbol{\mu}} _k^h = \frac{{\sum\nolimits_{i = 1}^M {I\left[ {{c_h}(i) = k} \right]} {{\boldsymbol{x}}_i}}}{{\sum\nolimits_{i = 1}^M {I\left[ {{c_h}(i) = k} \right]} }}. $

(4)

式中：$ I\left[ \cdot \right] $表示判断内部条件是否成立的指示函数，若内部条件成立，值为1，否则为0.

3）重复上述环节直至h达到阈值上限或聚簇均值稳定，略去上下标h，利用最终的聚簇信息为每个正态分布求解其对应的协方差矩阵：

$ {{\boldsymbol{\varSigma}} _k} = \frac{{\sum\nolimits_{i = 1}^M {I\left[ {c(i) = k} \right]} {\text{ }}({{\boldsymbol{x}}_i} - {{\boldsymbol{\mu}} _k}){{({{\boldsymbol{x}}_i} - {{\boldsymbol{\mu}} _k})}^{\mathrm{T}}}}}{{\sum\nolimits_{i = 1}^M {I\left[ {c(i) = k} \right]} }}. $

(5)

注意，K-means聚类层的上述步骤须分别应用在源点云和目标点云，在此之后，即可得到它们各自的正态分布组件参数$ {\varTheta _X} $和$ {\varTheta _Y} $.

2.2.2. KL散度计算层

利用聚簇均值之间的空间距离来搜索最近邻匹配正态分布，并计算它们的KL散度，以此作为识别非重叠区域的依据.

1）为每个源正态分布$ N({{\boldsymbol{\mu}} _i},{{\boldsymbol{\varSigma}} _i}) $搜索其在目标点云空间的最近邻匹配正态分布：

$ {d_H}(i) = \mathop {\arg \min }\limits_{k \in \{ 1,2, \cdots ,{K_Y}\} } {\left\| {{{\boldsymbol{R}}_{H - 1}}{{\boldsymbol{\mu}} _i}+{{\boldsymbol{t}}_{H - 1}} - {{\boldsymbol{\mu}} _k}} \right\|_2}. $

(6)

2）在获取最近邻匹配分布的基础上，求解两者的KL散度，以描述它们的数据相似性：

$ \begin{split} {\mathrm{k}}{{\mathrm{l}}_{i,{d_H}(i)}} =& {\mathrm{KL}}(N({{\boldsymbol{\mu}} _{{d_H}(i)}},{{\boldsymbol{\varSigma}} _{{d_H}(i)}})\left\| {N({{\boldsymbol{\mu}} _i'},{{\boldsymbol{\varSigma}} _i'})} \right.) = \\&\frac{1}{2}\Bigg[\log\; \frac{{\det \;({{\boldsymbol{\varSigma}} _i'})}}{{\det\; ({{\boldsymbol{\varSigma}} _{{d_H}(i)}})}} - D+{\mathrm{tr}}\;({({\varSigma _i'})^{ - 1}}{{\boldsymbol{\varSigma}} _{{d_H}(i)}}) + \\& {{\text{(}}{{\boldsymbol{\mu}} _{{d_H}(i)}} - {{\boldsymbol{\mu}} _i'}{\text{)}}^{\mathrm{T}}}{({{\boldsymbol{\varSigma}} _i'})^{ - 1}}{\text{(}}{{\boldsymbol{\mu}} _{{d_H}(i)}} - {{\boldsymbol{\mu}} _i'}{\text{)}}\Bigg]. \end{split}$

(7)

式中：D表示正态分布维度，等于点云的坐标维度3；$ \det\; ( \cdot ) $和$ {\mathrm{tr}}\;( \cdot ) $分别表示求解矩阵的行列式和矩阵的迹的操作；正态分布$ N({{\boldsymbol{\mu}} _i'},{\boldsymbol{\varSigma} _i'}) $表示$ N({{\boldsymbol{\mu}} _i},{\boldsymbol{\varSigma} _i}) $经过当前刚性变换$ ({{\boldsymbol{R}} _{H - {\text{1}}}},{{\boldsymbol{t}} _{H - 1}}) $后形成的正态参数，有$ {{\boldsymbol{\mu}} _i'} = {{\boldsymbol{R}} _{H - {\text{1}}}}{{\boldsymbol{\mu}} _i}+{{\boldsymbol{t}} _{H - 1}}, $$ {\boldsymbol{\varSigma} _i'} = {{\boldsymbol{R}} _{H - {\text{1}}}}{\boldsymbol{\varSigma} _i}{\boldsymbol{R}} _{_{H - {\text{1}}}}^{\mathrm{T}} $.

3）计算最近邻正态分布匹配对的权重系数，以判断其属于真实重叠的可靠性：

$ {w_{i,{d_H}(i)}} = \frac{{\min\; ({\mathrm{k}}{{\mathrm{l}}_{m,{d_H}(m)}})}}{{{\mathrm{k}}{{\mathrm{l}}_{i,{d_H}(i)}}}};\;\;\; m = 1,2, \cdots ,{K_X}. $

(8)

式（2）所需的正态分布参数集信息以及相似性权重系数信息全部求解完毕.

2.2.3. 李代数求解器层

将最终的求解过程映射至李代数空间，获取新一轮刚性变换结果.

为了便于书写，后续求解过程略去下标$ {d_H}(i) $等变量信息，以目标正态分布$ N({{\boldsymbol{\mu}} _j},{\boldsymbol{\varSigma} _j}) $表示源正态分布$ N({{\boldsymbol{\mu}} _i},{\boldsymbol{\varSigma} _i}) $的最近邻匹配，以$ {w_i} $表示两者之间的相似性权重系数. 对式（2）进行优化，可以形成如下新的目标函数：

$ \min\; L({{\boldsymbol{R}} _H},{{\boldsymbol{t}} _H}) = \sum\nolimits_{i = 1}^{{K_X}} {{w_i} ({\boldsymbol{r}}_i^{\mathrm{T}}{{\boldsymbol{\varOmega}} _i}{{\boldsymbol{r}}_i})} . $

(9)

式中：$ \left(\boldsymbol{R}_H, \boldsymbol{t}_H\right) $表示新一轮待优化的刚性变换配准结果；$ {{\boldsymbol{r}}_i} $和$ {{\boldsymbol{\varOmega}} _i} $分别表示残差和信息矩阵，分别满足$ {{\boldsymbol{r}}_i} = {{\boldsymbol{R}} _H}{{\boldsymbol{\mu}} _i}+{{\boldsymbol{t}} _H} - {{\boldsymbol{\mu}} _j} $，$ {{\boldsymbol{\varOmega}} _i} = {({{\boldsymbol{\varSigma}} _i'}+{{\boldsymbol{\varSigma}} _j})^{ - 1}} $. 由Zhu等^[31]的研究工作，可直接给出对应于李代数空间的新一轮的刚性变换配准结果：

$ {{\boldsymbol{\xi}} _H} = - {({\boldsymbol{G}})^\dagger }{\boldsymbol{b}}. $

(10)

式中：$ {( \cdot )^\dagger } $表示对矩阵求解其伪逆矩阵的操作. 除此之外，G和b的求解过程表示如下：

$ \left. \begin{gathered} {\boldsymbol{G}} = \sum\nolimits_{i = 1}^{{K_X}} {{w_i}} {\boldsymbol{G}}_i^{\mathrm{T}}{{\boldsymbol{\varOmega}} _i} {{\boldsymbol{G}}_i} , \\ {\boldsymbol{b}} = \sum\nolimits_{i = 1}^{{K_X}} {{w_i}{\boldsymbol{G}}_i^{\mathrm{T}}{{\boldsymbol{\varOmega}} _i} {\boldsymbol{r}}_i^0 }, \\ {{\boldsymbol{G}}_i} = \left[ {\begin{array}{*{20}{c}} { - {{({{\boldsymbol{\mu}} _i'})}^ \wedge }},&{{{\boldsymbol{I}}_3}} \end{array}} \right] , \\ {\boldsymbol{r}}_i^0 = {{\boldsymbol{\mu}} _i'} - {{\boldsymbol{\mu}} _j} . \\ \end{gathered} \right\} $

(11)

式中：${\boldsymbol{G}} \in {\bf{R}}^{6 \times 6} $，${\boldsymbol{b}} \in {\bf{R}}^{6 } $，${\boldsymbol{G}}_i \in {\bf{R}}^{3 \times 6} $，${\boldsymbol{r}}_i^0 \in {\bf{R}}^{3} $，$ {( \cdot )^ \wedge } $表示对三维向量求解其对应的反对称矩阵. 在式（11）的基础上，将式（10）的结果重新映射至李群空间，可以得到

$ \left. \begin{gathered} {{\boldsymbol{R}} _H} = \exp \;({({{\boldsymbol{\xi}} _{{H_R}}})^ \wedge }){{\boldsymbol{R}} _{H - 1}}, \\ {{\boldsymbol{t}} _H} = \exp \;({({{\boldsymbol{\xi}} _{{H_R}}})^ \wedge }){{\boldsymbol{t}} _{H - 1}}+{{\boldsymbol{\xi}} _{{H_t}}}. \\ \end{gathered} \right\} $

(12)

式中：$ {{\boldsymbol{\xi}} _{{H_R}}} $和$ {{\boldsymbol{\xi}} _{{H_t}}} $分别表示处于李代数空间中的新一轮刚性变换配准结果$ {{\boldsymbol{\xi}} _H}\; ({\boldsymbol{\xi}}_H \in {\bf{R}}^6) $所对应的旋转矩阵分量和平移矩阵分量；$ \exp\; ( \cdot ) $表示由李代数空间到李群空间的指数映射操作，基于微小扰动假设可以定义为$ \exp \;({({{\boldsymbol{\xi}} _{{H_R}}})^ \wedge }) = {{\boldsymbol{I}}_3}+{({{\boldsymbol{\xi}} _{{H_R}}})^ \wedge } $. 至此，李代数求解器层处理完毕.

综上，在K-means聚类层获取正态分布参数集的基础上，重复迭代KL散度计算层和李代数求解器层，直至整体大循环次数H达到其阈值上限或相邻2次求解得到的刚性变换矩阵没有明显差异，即可输出最终结果.

2.3. 算法实现和复杂度分析

规定h和H的阈值上限都为200. 在此基础上，本研究所提的基于正态分布相似性的双视角点云配准方法的具体实现如算法1所示.

算法1　基于正态分布相似性的双视角点云配准方法

输入：源点云X，目标点云Y，初始旋转矩阵和平移矩阵，

X和Y各自的初始聚簇中心

1）初始化h=1，H=1；

Repeat：

2）令h=h+1；

3）根据式（3）~（5），分别为源点云和目标点云生成正态分布参数集；

Until：h大于200或前后2次聚簇均值不再发生变化；

Repeat：

4）令H=H+1；

5）根据式（6）搜索最近邻匹配正态分布；

6）根据式（7）、（8）计算最近邻正态分布的相似性权重；

7）根据式（11）计算李代数求解器相关参数；

8）根据式（10）、（12）优化新一轮刚性变换配准结果；

Until：H大于200或前后2次配准结果没有明显变化；

输出：最优旋转矩阵R和平移矩阵t

在对所提方法进行复杂度分析之前，作出以下假设：1）待配准的2帧点云所包含的数据量不会相差太大，统一用N表示；2）待配准的2帧点云数据分布密度相对均匀，因此其对应的正态分布聚簇数量不会相差太大，统一用K表示.

2.3.1. 算法的时间复杂度分析

基于上述假设，所提方法的时间复杂度如表1所示，具体分析如下.

表 1 双视角点云配准方法的时间复杂度

Tab.1 Time complexity analysis of pair-wise point cloud registration

所属模块	具体操作	时间复杂度
K-means聚类	建立k-d树	$ 2O(hK \times \lg\; K) $
K-means聚类	更新聚簇并计算正态参数	$ 2O(hN\times \lg \;K) $
KL散度计算	搜索最近邻匹配正态分布	$ O(HK\times \lg \;K) $
KL散度计算	计算KL散度并生成权重	$ O(HK) $
李代数求解器	优化新一轮刚性变换	$ O(HK) $

1） K-means聚类层. 为源点云当中的K个聚簇中心建立k-d树，复杂度为$ O(K\times \lg \;K) $；在此基础上，为每个数据搜索对应的最近邻簇中心，复杂度为$ O(N\times \lg K) $. 正态参数的计算可以一并包含在上述过程中，不占用额外时间. 对目标点云的处理思路完全一致，因此这一层的整体时间复杂度可以记作$ 2(O(hK\times \lg K)+O(hN\times \lg K)) $.

2） KL散度计算层. 为源点云的每一个正态分布搜索其在目标点云空间的最近邻匹配，复杂度为$ O(K\times \lg K) $；遍历每一对最近邻正态分布，计算其KL散度，复杂度记作$ O(K) $. 由于整体迭代次数为H，这一层的时间复杂度最终可以表示为$ O(HK\times \lg K)+O(HK) $.

3）李代数求解器层. 求解过程相当于遍历所有正态组件参数，整体复杂度为$ O(HK) $.

2.3.2. 算法的空间复杂度分析

令D表示输入数据的维度，结合上述假设，可知2帧点云各自占用的空间大小为$ O(DN) $，后续环节分析如下.

1） K-means聚类层. 该层的目标是更新簇中心和协方差矩阵，两者所占用的空间大小直接关系于正态聚簇个数K，因此复杂度可以分别表示为$ O(DK) $和$ O({D^2}K) $. 另外，由于更新过程依赖于数据间的最近邻匹配搜索，须将簇中心存储到k-d树结构，因此这一过程会额外占用$ O(DK) $大小的空间. 注意，源点云和目标点云要分别进行K-means聚类操作，两者所耗用的空间资源是相似的，都适用于上述分析.

2） KL散度计算层. 该层的任务是求解最近邻匹配正态分布聚簇之间的相似性权重，对应的空间复杂度与正态分布数量呈线性相关，因此可以直接表示为$ O(K) $.

3）李代数求解器层. 该层直接利用上述环节所获取的正态参数集信息和相似性权重信息来优化刚性变换配准结果，期间没有大量额外空间须被耗用，因此复杂度用$ O(1) $表示.

基于正态分布相似性的双视角点云配准方法空间复杂度如表2所示. 本研究输入点云的维度D=3，所提方法的空间复杂度与K可以看作线性关系，而K一般又远小于点云数据量N，因此所提方法空间复杂度整体较低.

表 2 双视角点云配准方法的空间复杂度

Tab.2 Spatial complexity analysis of pair-wise point cloud registration

所属模块	具体操作	空间复杂度
K-means聚类	建立k-d树	$ 2O(DK) $
K-means聚类	求解正态分布聚簇中心	$ 2O(DK) $
K-means聚类	求解正态分布协方差矩阵	$ 2O({D^2}K) $
KL散度计算	计算KL散度并生成权重	$ O(K) $
李代数求解器	优化新一轮刚性变换	$ O(1) $

3. 实验结果

为了验证本研究所提出的基于正态分布相似性的双视角点云配准方法的性能，使用4个斯坦福数据集^[32]（Bunny、Armadillo、Buddha、Dragon）进行消融实验、精度实验、效率实验等. 4个数据集包含的点云帧数分别为10、12、15和15，每帧点云平均拥有的数据量分别为36227、25635、73267和31280. 须说明，同一数据集相邻2帧点云的数据量一般相差不大，但其数据重叠百分比为30%~100%，因此对于同一数据集相邻2帧两两配准的平均性能结果将在一定程度上反映算法能否有效处理部分重叠的点云配准任务.

将本研究所提方法与其他8种较流行的算法进行对比，具体包括TrICP^[10]、NDT^[22]、SpICP^[11]、SmICP^[20]、AAICP^[13]、刚性变换一致（rigid transformation consensus, RTC）算法^[16]、镜像迭代（mirrored iterative closest point, MICP）算法^[17]、RSICP^[21]，所有测试均在Matlab R2021a上实现. 如无特殊说明，在实验过程中，将对同一数据集的相邻2帧点云两两配准，最终求解指标的平均值以说明不同算法的性能. 所用的精度指标包括平均均方根误差RMSE和总体配准成功率SR，定义如下：

$ \left. \begin{gathered} {\mathrm{RMSE}} = \frac{1}{F}\sum\nolimits_{f = 1}^F {\sqrt {\frac{1}{{{N_f}}}\sum\nolimits_{i = 1}^{{N_f}} {\left\| {T_f^{\mathrm{E}}{{\boldsymbol{x}}_i} - T_f^{\mathrm{G}}{{\boldsymbol{x}}_i}} \right\|} _2^2} } , \\ {\mathrm{SR}} = \frac{1}{F}\sum\nolimits_{f = 1}^F {I\left[ {{{\mathrm{RMSE}}_f} < \varepsilon \times {{\mathrm{RMSE}}_{{\mathrm{fini}}}}} \right]} {\text{ }}. \\ \end{gathered} \right\} $

(13)

式中：F表示数据集两两配准的总次数，具体为该数据集的点云帧数减去1；$ {N_f} $表示第$ f $次配准即数据集第$ f $和第$ {f+1} $帧点云配准时源点云所包含的数据量大小；$ {\boldsymbol{T}}_f^{\mathrm{E}} $和$ {\boldsymbol{T}}_f^{\mathrm{G}} $分别表示第f次配准过程中刚性变换矩阵的估计结果和真值结果；$ \varepsilon $为阈值系数，设置为0.15；$ {{\mathrm{RMSE}}_f} $表示第f次配准后的均方根误差结果；$ {\mathrm{RMSE}}_{{\mathrm{fini}}} $表示第f次配准前的初始误差结果.

对于待配准的某特定数据集，指标RMSE能够反映方法对不同重叠百分比下的各相邻帧进行配准时的平均精确度，而指标SR能够反映方法在配准过程中保持良好性能的稳定性.

3.1. 参数选择

本研究所提方法的参数K，表示源点云或目标点云的正态分布聚簇个数. 一般而言，这个数量难以直接估计，但可以通过确定每个分布平均包含的数据量进行反推. 令$ {N_{\mathrm{P}}} $表示每个正态分布聚簇平均包含的数据量，有$ K \approx N/{N_{\mathrm{P}}} $. 由该式不难推断出，当$ {N_{\mathrm{P}}} $过小时，K过大，算法依然相当于传统“点到点”配准，可能会耗用大量时间资源；当$ {N_{\mathrm{P}}} $过大时，K过小，原始点云当中的大量局部细节会丢失，无法保证配准精度的稳定性. 当$ {N_{\mathrm{P}}} $取不同值时各数据集第1、2帧的配准精度和用时情况如图3所示. 为了兼顾精度和效率，本研究所提方法在所有配准过程中设置$ {N_{\mathrm{P}}} = 36 $，即将每个正态分布聚簇当中平均包含的数据量视作固定常量.

图 3

图 3 正态分布聚簇内平均数据量不同时的配准结果

Fig.3 Registration results with different mean data volume within normal distribution clusters

3.2. 精度和效率实验

在刚性变换的真值结果中加入随机扰动，以生成初始配准参数. 在斯坦福的不同数据集上对各种算法进行测试，所用指标包括平均均方根误差RMSE、平均配准成功率SR、平均配准用时t. 配准结果如表3所示. 表中，粗体表示配准中的最优值，下划线表示次优值，Initial表示未使用任何配准方法. 为了便于观察差异，如图4所示给出了用不同方法配准Bunny第2、3帧和Dragon第12、13帧的结果示意图，这些点云数据之间的重叠率约为65%. 图中，RP表示经配准后2帧点云的拟合度，颜色越偏蓝色表示算法越能处理部分重叠的点云配准任务，越偏红色表示处理能力越差. 此外，为了进一步说明不同方法的运行效率，如图5所示给出了它们在不同数据集上配准的平均用时示意.

表 3 不同方法在4个数据集上的精度结果对比

Tab.3 Accuracy comparison of different methods on four datasets

配准方法	Bunny		Armadillo		Buddha		Dragon
配准方法	RMSE	SR	RMSE	SR	RMSE	SR	RMSE	SR
Initial	13.4487	—	19.3405	—	16.9126	—	17.9725	—
TrICP	2.759 4	0.78	0.1728	1.00	2.5540	0.64	0.9472	0.93
NDT	12.7811	0.00	17.4936	0.00	6.3208	0.50	17.6523	0.00
SpICP	5.7842	0.56	1.7818	0.91	5.1119	0.21	4.4801	0.57
SmICP	24.8743	0.22	0.2554	1.00	0.3101	1.00	1.1178	0.86
AAICP	16.2111	0.11	0.2774	1.00	0.2705	1.00	0.9205	0.93
RTC	3.8347	0.56	2.5884	0.73	1.7657	0.86	1.6472	0.79
MICP	5.4027	0.44	0.2120	1.00	0.249 6	1.00	0.9594	0.93
RSICP	27.9910	0.56	0.1591	1.00	0.2510	1.00	0.335 2	1.00
本研究方法	0.4249	1.00	0.161 9	1.00	0.2472	1.00	0.3289	1.00

图 4

图 4 不同方法配准结果的颜色渐变示意图

Fig.4 Color gradient illustration of registration results of different methods

图 5

图 5 不同双视角点云配准方法的运行时间对比

Fig.5 Comparison of relative runtime of different pair-wise point cloud registration methods

在配准精度上，如表3和图4所示，本研究所提方法基本上可以获得最优的配准结果，且其稳定性较高，在4个数据集均取得了100%的配准成功率. 在用时效率上，如图5所示（仅截取50 s以内部分，耗时最长的SpICP在4个数据集上平均配准时间分别为266、147、645、152 s），本研究所提方法在整体上处于中上游水平，虽然略慢于NDT和MICP算法，但在精度上却远超它们.

TrICP、SpICP、AAICP和MICP都是对传统ICP的改进，隶属于“点到点”配准模式. TrICP主动引入重叠百分比参数来预测重叠部分，虽然有效，但依然仅仅利用点对之间的空间距离来完成配准任务，一旦非重叠区域面积增大或更加精细复杂，很可能因陷入大量迭代而致使效率和精度同时下降. SpICP基于稀疏表示理论来增强算法整体的鲁棒性，但由于引入了增广拉格朗日函数，算法效率处于对比当中的最差. AAICP将安德鲁加速算法引入传统ICP框架，效率得到有效提升；此外，该算法基于威尔斯函数设计了鲁棒性误差度量，来避免噪声值和异常值的不利影响，但其准确性同样依赖于真实重叠百分比. MICP在配准过程中须为数据建立可靠的邻域空间，算法本身非常高效，但最终的配准结果并不稳定，尤其是在真实重叠百分比较低时，其搜索的邻域空间往往并不可靠甚至完全错误，致使最终的配准精度大幅下降.

NDT和RTC都是基于概率的配准算法. NDT隶属于“点到分布”配准，非常高效，但须对点云空间进行单元体素划分，提前确定合适的网格大小，否则数据分布极不均匀，导致正态参数求解的偏差过大. RTC会为最近邻匹配点对分配一个后验概率，以指示其属于真实匹配的可能性大小，在本质上仍属于“点到点”配准，虽能在一定程度上解决部分重叠问题，但其本身须固定若干高斯参数，参数选取不当容易导致配准失败.

SmICP和RSICP虽然也是对ICP的改进，但将传统“点到点”转化为“点到面”配准，使用对称目标函数来尽可能优化点面距离，从而能够提升配准的精度. 此外，相较于SmICP，RSICP由于引入鲁棒性对称度量，既能拓宽收敛域又能提升收敛速度，因此整体性能会更强. 然而这类“点到面”的算法都须估计法向量，在迭代过程中可能会耗用更多时长，且在处理像Bunny数据集这样的低重叠配准任务时，其度量结果可能会出现较大的偏差，致使配准完全失败.

综上所述，相较于其他方法，本研究所提方法有性能上的优势，能兼顾效率和精度. 具体来说，算法采用“分布到分布”配准技巧，能够避免频繁对所有数据进行遍历迭代，耗时大幅降低；算法利用最近邻正态分布之间的相似性来识别潜在的非重叠区域，从而能有效完成部分重叠配准任务，提升配准精度.

3.3. 鲁棒性实验

为了验证本研究所提方法对于数据异常值的鲁棒性，为斯坦福各个数据集分别添加信噪比SNR为50、25 dB的高斯噪声. 此外，为了保证公平性，所有方法在实验过程中须对同一数据集的相邻2帧进行两两配准，最终求解平均均方根误差RMSE作为该数据集的一次实验数据，在此基础上对该数据集独立重复进行共计30次实验后求取实验数据均值. 最终的实验结果如表4所示，其中指标RMSE以“平均值±标准差”的形式进行表示，粗体部分表示对比中的最优值，下划线部分表示对比中的次优值. 由实验结果可知，在加入高斯噪声后，不同方法的配准精度都会出现一定程度的扰动，且噪声越剧烈，扰动的幅度可能越大. 在相同的噪声水平下，本研究所提方法能保持相对更稳定的配准效果. 一方面，该算法利用K-means进行正态聚类和参数求解，灵活性和可靠性较高，即使出现噪声值或异常值，其正态分布拟合也不会出现较大波动，配准所用的均值和协方差矩阵信息相对稳定；另一方面，该算法使用基于KL散度的权重系数来识别非重叠的数据区域，能够将配准的注意力集中在真实重叠部分，从而规避异常数据的影响，进一步提升配准过程的鲁棒性.

表 4 不同高斯噪声（50、25 dB）下的配准结果

Tab.4 Comparison of registration results under Gaussian noise of 50 dB and 25 dB

配准方法	SNR/dB	RMSE
配准方法	SNR/dB	Bunny	Armadillo	Buddha	Dragon
TrICP	50	2.7612±0.0017	0.1727±0.0003	2.5530±0.0009	0.9465±0.0004
TrICP	25	2.7563±0.0081	0.1733±0.0014	2.4407±0.0199	0.9474±0.0028
NDT	50	12.7890±0.0062	22.2961±0.0181	6.6362±0.4974	24.0140±0.0787
NDT	25	12.7977±0.0232	17.3343±0.3254	6.8398±0.5018	17.5718±0.0923
SpICP	50	5.7519±0.0089	1.6955±0.0031	5.2335±0.0019	4.6071±0.0038
SpICP	25	5.8516±0.0171	1.7163±0.0020	5.1914±0.0025	4.5882±0.0046
SmICP	50	20.3221±2.2295	0.2554±0.0002	0.3103±0.0003	1.1177±0.0005
SmICP	25	20.4048±1.6212	0.2548±0.0009	0.3108±0.0003	1.1218±0.0020
AAICP	50	16.2103±0.0013	0.2774±0.0002	0.2706±0.0003	0.9206±0.0003
AAICP	25	16.2135±0.0051	0.2778±0.0009	0.2709±0.0006	0.9204±0.0010
RTC	50	3.8398±0.0047	2.5661±0.0134	1.7657±0.0004	1.6704±0.0372
RTC	25	3.8460±0.0246	2.6874±0.1911	1.7591±0.0027	1.4611±0.1315
MICP	50	5.3758±0.0735	0.2111±0.0014	0.2535±0.0010	0.9608±0.0035
MICP	25	5.4674±0.0644	0.2111±0.0022	0.2553±0.0019	0.9907±0.0473
RSICP	50	27.9803±0.0200	0.1592±0.0002	0.2510±0.0003	0.3353±0.0005
RSICP	25	27.7358±0.5571	0.1583±0.0006	0.2499±0.0007	0.3354±0.0012
本研究方法	50	0.4265±0.0078	0.1611±0.0017	0.2470±0.0009	0.3287±0.0012
本研究方法	25	0.4294±0.0185	0.1615±0.0015	0.2478±0.0012	0.3293±0.0028

NDT和RTC均以高斯概率模型为基础进行配准. NDT使用规则单元来划分正态分布，一旦噪声数据增加，固定网格的数据信息变化相对剧烈，因此正态参数求解误差增大，方法的鲁棒性下降. RTC为每一对最近邻匹配点计算其是否属于“真实匹配”的后验概率，而在噪声数据较多时，点对的概率评估可靠性势必会受到影响，致使最终的配准结果波动幅度较大.

TrICP、SpICP、SmICP、AAICP、MICP和RSICP均为基于ICP的改进算法. 其中，TrICP和MICP能够将注意力放在数据的重叠区，对于异常信息有一定的过滤和清除作用，因此具有一定的鲁棒性，但当高斯噪声加剧时，这种识别能力可能会受到一定影响. SpICP和AAICP都在传统“点到点”配准目标函数的基础上进行改进和优化，前者使用l_p范数替换原先的l₂范数，后者引入基于威尔斯函数的鲁棒性误差度量，都能在一定程度上增强配准结果的稳定性. SmICP和RSICP均采用“点到面”的方式对点云数据进行配准，能够有效规避“点到点”配准的信息度量误差，且两者都使用对称型目标函数来降低算法对于噪声数据的敏感程度. 然而，这些算法虽然在其他数据集上有较好的配准表现，但在处理像Bunny数据集这样重叠百分比极低的数据时依然容易出现剧烈的精度偏移.

3.4. 消融实验

本研究所提方法因隶属“分布到分布”配准而拥有较高的效率，因引入由KL散度推导的相似性权重系数而拥有较好的精度. 从这2个要素出发，可以设计以下实验：1）“点到分布”+“无相似性权重”（P2DNDT），即将K-means聚类层只应用在目标点云，将配准问题转化为源数据的概率生成和最大问题，中间不涉及KL散度计算，但利用李代数求解器对结果进行求解；2）“分布到分布”+“无相似性权重”（D2DNDT），即配准任务不包含KL散度计算层，除此之外的过程皆与本研究所提方法一致. 须注意，由于KL散度计算须面向“分布到分布”，不能设计“点到分布”+“有相似性权重”的实验. 最终的消融实验结果如表5所示，用于比较的实验指标包括平均均方根误差RMSE、平均配准用时t和整体配准成功率SR. 表中，粗体表示结果中的最优值. 由实验结果可知，D2DNDT的效率略高于所提方法的，这是因为其无须经过KL散度计算层处理，但其也因此无法对非重叠区域进行识别，配准精度大幅下降. 与之相比，P2DNDT须遍历所有源数据点，信息较为全面可靠，因此在精度方面会有一定保障，但依然低于本研究所提方法的，且该算法须重复遍历源点云当中的所有数据，会消耗大量的时间.

表 5 在不同数据集上消融实验的结果

Tab.5 Results of ablation experiments on different datasets

数据集	Initial	P2DNDT			D2DNDT			本研究方法
数据集	RMSE	RMSE	t /s	SR	RMSE	t /s	SR	RMSE	t /s	SR
Bunny	13.4487	4.4592	19.75	0.33	13.7552	3.99	0.33	0.4249	3.74	1.00
Armadillo	19.3405	0.6452	14.58	0.91	0.7699	1.64	0.82	0.1619	2.33	1.00
Buddha	16.9126	0.4295	42.88	1.00	0.4184	4.35	1.00	0.2472	9.47	1.00
Dragon	17.9725	0.8273	15.99	1.00	3.6687	1.93	0.71	0.3289	3.56	1.00

3.5. 真实场景实验

为了进一步论证所提算法的性能，在真实数据集GazeboS和GazeboW^[33]上测试，两者各包含约30帧的点云，每帧点云含3~5万个数据，这些数据由配备扫描仪的机器人采集得到. 如表6所示为不同方法对2个数据集相邻2帧两两配准后的实验结果. 表中，粗体表示结果中的最优值，下划线表示结果中的次优值. 可以看出，除了NDT之外，多数算法都能得到精度相对良好的配准结果，这是因为点云数据是由机器人沿二维闭合路径按照较小的视角变化采集得到的，相邻2帧的重叠百分比较高，因此它们的配准精度均较高. 整体来看，本研究所提方法在精度上相较于TrICP、SmICP、AAICP等算法具有明显的优势；在2个数据集上精度分别排名第2的RSICP算法和SpICP算法，其平均配准效率均落后于本研究所提方法，其中，RSICP的效率不足本研究所提方法的60%，而SpICP的效率甚至不足1%. 综上所述，本研究所提方法在真实环境场景配准方面也具有较好的应用潜力.

表 6 在真实场景数据集上的配准结果

Tab.6 Registration results of different real scene datasets

配准方法	GazeboS			GazeboW
配准方法	RMSE	t /s	SR	RMSE	t /s	SR
Initial	0.2689	—	—	0.3685	—	—
TrICP	0.0279	13.88	0.84	0.0204	2.83	1.00
NDT	0.1489	2.85	0.26	0.0784	1.97	0.63
SpICP	0.0302	1242.40	0.77	0.019 6	402.32	1.00
SmICP	0.0709	34.42	0.74	0.0255	17.37	0.97
AAICP	0.0299	5.77	0.71	0.0228	2.29	1.00
RTC	0.0261	5.66	0.74	0.0229	1.94	1.00
MICP	0.0277	2.24	0.71	0.0575	1.54	0.93
RSICP	0.0211	26.04	0.93	0.0197	7.88	0.97
本研究方法	0.0208	9.48	0.90	0.0192	4.25	1.00

4. 结　论

（1）所提算法通过引入K-means聚类策略，将传统“点到点”配准模式转化为“分布到分布”的配准问题，不仅规避了传统正态分布变换规则划分体素单元致使数据不均的弊端，也有效保证了配准任务的整体效率.

（2）为了解决部分重叠配准任务，所提算法引入相似性权重对非重叠区域进行识别. 方法利用均值之间的距离来搜索最近邻匹配分布，并在此基础上计算KL散度，散度越小，意味2个分布越相似，为它们赋予权重则越高. 经过这一环节处理，可以有效提升配准结果的精度.

（3）在公开数据集上的结果表明，相较于目前主流的基于ICP或“点到点”配准的算法，所提方法能够在兼顾效率的同时提升配准结果的精度，且可以有效容忍噪声的存在，在真实数据集上也表现出了良好的适应能力和应用能力.

（4）尽管所提方法在性能上有一定优势，但也存在待改进的方面. 由于方法更适合用于庞杂精细模型的双视角配准任务，因此如何将其应用至深度网络进行嵌入特征分析，或如何将其应用至多视角点云配准，将成为后续工作的重点.

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

GOJCIC Z, ZHOU C, WEGNER J D, et al. Learning multiview 3D point cloud registration [C]// IEEE/CVF Conference on Computer Vision and Pattern Recognition. Seattle: IEEE, 2020: 1759–1769.

[2]

LI P, WANG R, WANG Y, et al

Evaluation of the ICP algorithm in 3D point cloud registration

[J]. IEEE Access, 2020, 8: 68030- 68048

DOI:10.1109/ACCESS.2020.2986470 [本文引用: 1]

[3]

CHEN X, MA H, WAN J, et al. Multi-view 3D object detection network for autonomous driving [C]// IEEE Conference on Computer Vision and Pattern Recognition. Honolulu: IEEE, 2017: 6526–6534.

[4]

ZHENG Y, LI Y, YANG S, et al

Global-PBNet: a novel point cloud registration for autonomous driving

[J]. IEEE Transactions on Intelligent Transportation Systems, 2022, 23 (11): 22312- 22319

DOI:10.1109/TITS.2022.3153133 [本文引用: 1]

[5]

范光宇, 宫宇宸, 饶蕾, 等

基于灰度相似性的激光点云与全景影像配准

[J]. 浙江大学学报: 工学版, 2022, 56 (8): 1633- 1639

FAN Guangyu, GONG Yuchen, RAO Lei, et al

Registration of laser point cloud and panoramic image based on gray similarity

[J]. Journal of Zhejiang University: Engineering Science, 2022, 56 (8): 1633- 1639

[6]

龚靖渝, 楼雨京, 柳奉奇, 等

三维场景点云理解与重建技术

[J]. 中国图象图形学报, 2023, 28 (6): 1741- 1766

DOI:10.11834/jig.230004 [本文引用: 1]

GONG Jingyu, LOU Yujing, LIU Fengqi, et al

Scene point cloud understanding and reconstruction technologies in 3D space

[J]. Journal of Image and Graphics, 2023, 28 (6): 1741- 1766

DOI:10.11834/jig.230004 [本文引用: 1]

[7]

BESL P J, MCKAY N D

A method for registration of 3-D shapes

[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1992, 14 (2): 239- 256

DOI:10.1109/34.121791 [本文引用: 1]

[8]

BIBER P, STRASSER W. The normal distributions transform: a new approach to laser scan matching [C]// IEEE International Conference on Intelligent Robots and Systems. Las Vegas: IEEE, 2003: 2743–2748.

[9]

CHARLES R Q, HAO S, MO K, et al. PointNet: deep learning on point sets for 3D classification and segmentation [C]// IEEE Conference on Computer Vision and Pattern Recognition. Honolulu: IEEE, 2017: 77–85.

[10]

CHETVERIKOV D, SVIRKO D, STEPANOV D, et al. The trimmed iterative closest point algorithm [C]// International Conference on Pattern Recognition. Quebec City: IEEE, 2002: 545–548.

[11]

BOUAZIZ S, TAGLIASACCHI A, PAULY M

Sparse iterative closest point

[J]. Computer Graphics Forum, 2013, 32 (5): 113- 123

DOI:10.1111/cgf.12178 [本文引用: 2]

[12]

PAVLOV A L, OVCHINNIKOV G W, DERBYSHEV D Y, et al. AA-ICP: iterative closest point with Anderson acceleration [C]// IEEE International Conference on Robotics and Automation. Brisbane: IEEE, 2018: 3407–3412.

[13]

ZHANG J, YAO Y, DENG B

Fast and robust iterative closest point

[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2022, 44 (7): 3450- 3466

DOI:10.1016/j.compeleceng.2022.108098 [本文引用: 1]

[14]

YANG J, LI H, CAMPBELL D, et al

Go-ICP: a globally optimal solution to 3D ICP point-set registration

[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2016, 38 (11): 2241- 2254

DOI:10.1109/TPAMI.2015.2513405 [本文引用: 1]

[15]

SHI X, LIU T, HAN X

Improved iterative closest point (ICP) 3D point cloud registration algorithm based on point cloud filtering and adaptive fireworks for coarse registration

[J]. International Journal of Remote Sensing, 2020, 41 (8): 3197- 3220

DOI:10.1080/01431161.2019.1701211

[16]

LI Z, WANG C, MA J, et al

Point set registration via rigid transformation consensus

[J]. Computers and Electrical Engineering, 2022, 101: 108098

[17]

徐文菲, 金莉, 韩旭, 等

面向缺失点云配准的镜像迭代最近点算法

[J]. 西安交通大学学报, 2023, 57 (7): 201- 212,220

DOI:10.7652/xjtuxb202307019 [本文引用: 2]

XU Wenfei, JIN Li, HAN Xu, et al

Mirrored iterative closest point algorithm for missing point cloud registration

[J]. Journal of Xi’an Jiaotong University, 2023, 57 (7): 201- 212,220

DOI:10.7652/xjtuxb202307019 [本文引用: 2]

[18]

CHEN Y, MEDIONI G

Object modelling by registration of multiple range images

[J]. Image and Vision Computing, 1992, 10 (3): 145- 155

DOI:10.1016/0262-8856(92)90066-C [本文引用: 1]

[19]

RUSINKIEWICZ S, LEVOY M. Efficient variants of the ICP algorithm [C]// 3rd International Conference on 3-D Digital Imaging and Modeling. Quebec City: IEEE, 2001: 145–152.

[20]

RUSINKIEWICZ S

A symmetric objective function for ICP

[J]. ACM Transactions on Graphics, 2019, 38 (4): 1- 7

DOI:10.1016/j.isprsjprs.2022.01.019 [本文引用: 2]

[21]

LI J, HU Q, ZHANG Y, et al

Robust symmetric iterative closest point

[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2022, 185: 219- 231

[22]

MAGNUSSON M, LILIENTHAL A, DUCKETT T

Scan registration for autonomous mining vehicles using 3D-NDT

[J]. Journal of Field Robotics, 2007, 24 (10): 803- 827

DOI:10.1002/rob.20204 [本文引用: 2]

[23]

STOYANOV T, MAGNUSSON M, LILIENTHAL A J. Point set registration through minimization of the L2 distance between 3D-NDT models [C]// IEEE International Conference on Robotics and Automation. Saint Paul: IEEE, 2012: 5196–5201.

[24]

AOKI Y, GOFORTH H, SRIVATSAN R A, et al. PointNetLK: robust and efficient point cloud registration using PointNet [C]// IEEE/CVF Conference on Computer Vision and Pattern Recognition. Long Beach: IEEE, 2019: 7156–7165.

[25]

WANG Y, SOLOMON J. Deep closest point: learning representations for point cloud registration [C]// IEEE/CVF International Conference on Computer Vision. Seoul: IEEE, 2019: 3522–3531.

[26]

WANG Y, SUN Y, LIU Z, et al

Dynamic graph CNN for learning on point clouds

[J]. ACM Transactions on Graphics, 2019, 38 (5): 1- 12

[27]

YUAN W, ECKART B, KIM K, et al. DeepGMR: learning latent Gaussian mixture models for registration [C]// European Conference on Computer Vision. Glasgow: Springer, 2020: 733–750.

[28]

XU H, LIU S, WANG G, et al. OMNet: learning overlapping mask for partial-to-partial point cloud registration [C]// IEEE/CVF International Conference on Computer Vision. Montreal: IEEE, 2021: 3112–3121.

[29]

MEI G, TANG H, HUANG X, et al. Unsupervised deep probabilistic approach for partial point cloud registration [C]// IEEE/CVF Conference on Computer Vision and Pattern Recognition. Vancouver: IEEE, 2023: 13611–13620.

[30]

ZHU J, JIANG Z, EVANGELIDIS G D, et al

Efficient registration of multi-view point sets by K-means clustering

[J]. Information Sciences, 2019, 488: 205- 218

DOI:10.1016/j.ins.2019.03.024 [本文引用: 1]

[31]

ZHU J, MU J, YAN C B, et al

3DMNDT: 3D multi-view registration method based on the normal distributions transform

[J]. IEEE Transactions on Automation Science and Engineering, 2024, 21 (1): 488- 501

DOI:10.1109/TASE.2022.3225679 [本文引用: 1]

[32]

Stanford University. The Stanford 3D scanning repository [EB/OL]. (2014-08-19) [2021-10-01]. https://graphics.stanford.edu/data/3Dscanrep/.