点云对齐/轨迹对齐方法及论文讲解

轨迹对齐方法及论文讲解"/>

点云对齐/轨迹对齐方法及论文讲解

文章目录

• • 点云对齐/轨迹对齐方法及论文讲解
  • • 1 Umeyama算法
    • 2 Horn四元数闭合求解算法
    • Marta的迭代最小二乘解法

点云对齐/轨迹对齐方法及论文讲解

已知两组对应点云，求解两组点云之间的坐标系变换。或者计算里程计的精度问题，得到里程计估计的时间离散位姿后以及轨迹位姿真值后，如何评判里程计的精度。
本文介绍三种方法：
1.Umeyama算法
参考文献：
《Least-Squares Fitting of Two 3-D Point Sets》
2.Horn四元数闭合求解算法
参考文献：
《Closed-form solution of absolute orientation using unit quaternions》
3.Marta的迭代最小二乘解法
参考文献：
《Trajectory Alignment and Evaluation in SLAM: Horn’s Method vs Alignment on the Manifold》

1 Umeyama算法

原理说明
代码实现及工具使用

2 Horn四元数闭合求解算法

使用单位四元数实现求解旋转闭式解。
设两组点集分别为 r l \bm r_l rl​和 r r \bm r_r rr​，则其误差方程构成如下：
e i = r ( r , i ) − s R ( r ( l , i ) ) − r 0 \bm e_i=\bm r_(r,i)-sR(\bm r_(l,i) )-\bm r_0 ei​=r(​r,i)−sR(r(​l,i))−r0​
其中待求解的变量为旋转矩阵R，平移向量 r 0 \bm r_0 r0​及尺度因子 s s s

（1）求解最少满足的点数
对齐两个轨迹或者对齐两组点云所要求解的自由度包含两组坐标系的3个平移，3个旋转及一个放缩尺度共7个自由度。3个点能建立9个约束(一个三维点对应3个约束)。因此3个点则能够满足求解约束。

（2）点集质心约束及平移求解
设两组点集质心对应计算如下：

误差构成的最小二乘问题：

通过对 r 0 \bm r_0 r0​求导可得:
r 0 = r ‾ r − s R r ‾ l \bm r_0=\bm {\overline{r}}_r-sR\bm {\overline{r}}_l r0​=rr​−sRrl​
在两组点集的质心处建立新的坐标系，轴向与原坐标系一致，从而得到点集在新坐标系下的坐标：
r l , i ′ = r l , i − r ‾ l , r r , i ′ = r r , i − r ‾ r \bm r'_{l,i}=\bm r_{l,i}-\bm {\overline{r}}_l,\bm r'_{r,i}=\bm r_{r,i}-\bm {\overline{r}}_r rl,i′​=rl,i​−rl​,rr,i′​=rr,i​−rr​
将质心代入到误差方程，可消去 r 0 r_0 r0​，误差方程变为：
e i ′ = r r , i ′ − s R r l , i ′ \bm e'_i=\bm r'_{r,i}-sR\bm r'_{l,i} ei′​=rr,i′​−sRrl,i′​

（3）尺度求解
1）由新的误差方程构成最小二乘问题：

得到关于 s s s的一元二次方程，解得：

2）未知R阵求解s
将误差方程进行如下调整：

由于 R T R = I R^T R=I RTR=I将关于 s s s与 R R R分离，构成关于 s s s的最小二乘问题，从而在未知 R R R阵的情况下得到 s s s值，即

另外，同时需满足：
min ⁡ R ∑ i = 1 n r r , i ′ T R r l , i ′ \mathop{\min}\limits_{R} \sum_{i=1}^{n} \bm r_{r,i}^{'T}R\bm r'_{l,i} Rmin​i=1∑n​rr,i′T​Rrl,i′​

（4）三点求解旋转
在已知两组坐标系之间的旋转变换后，求解平移及尺度因子是相对比较简单的。因此求解 R R R阵是该问题的关键。
设3个不共线点在两组坐标系下的表示为： r l , 1 , r l , 2 , r l , 3 ∈ r l \bm r_{l,1}, \bm r_{l,2}, \bm r_{l,3}∈\bm r_l rl,1​,rl,2​,rl,3​∈rl​及 r r , 1 , r r , 2 , r r , 3 ∈ r r \bm r_{r,1}, \bm r_{r,2}, \bm r_{r,3}∈r_r rr,1​,rr,2​,rr,3​∈rr​

利用3个点分别建立新的坐标系，设 r l , 1 r_{l,1} rl,1​为坐标系l的原点， r l , 1 → r l , 2 r_{l,1}→r_{l,2} rl,1​→rl,2​为x轴方向，y轴位于三个点构成的平面上，垂直于x轴，z轴符合右手定则正交与x轴与z轴。如下示意图：

令： x l = r l , 2 − r l , 1 ， x ^ l = x l / ‖ x l ‖ \bm x_l=\bm r_{l,2}-\bm r_{l,1}，\hat {\bm x}_l=\bm x_l/‖\bm x_l‖ xl​=rl,2​−rl,1​，x^l​=xl​/‖xl​‖，则 x ^ l \hat {\bm x}_l x^l​为指向x轴的单位向量
令： y l = ( r l , 3 − r l , 1 ) − [ ( r l , 3 − r l , 1 ) ⋅ x ^ l ] x ^ l ， y ^ l = y l / ‖ y l ‖ \bm y_l=(\bm r_{l,3}-\bm r_{l,1})-[(\bm r_{l,3}-\bm r_{l,1}) \cdot \bm {\hat x}_l ] \bm {\hat x}_l，\bm {\hat y}_l=\bm y_l/‖\bm y_l ‖ yl​=(rl,3​−rl,1​)−[(rl,3​−rl,1​)⋅x^l​]x^l​，y^​l​=yl​/‖yl​‖
[ ( r l , 3 − r l , 1 ) ⋅ x ^ l ] x ^ l [(\bm r_{l,3}-\bm r_{l,1})\cdot\bm {\hat x}_l ] \bm {\hat x}_l [(rl,3​−rl,1​)⋅x^l​]x^l​为 r l , 3 − r l , 1 \bm r_{l,3}-\bm r_{l,1} rl,3​−rl,1​向量在x轴方向上的投影向量，即为坐标系原点到第三点 r l , 3 \bm r_{l,3} rl,3​在x轴的投影点的向量，因此 ( r l , 3 − r l , 1 ) − [ ( r l , 3 − r l , 1 ) ⋅ x ^ l ] x ^ l (\bm r_{l,3}-\bm r_{l,1})-[(\bm r_{l,3}-\bm r_{l,1}) \cdot \bm {\hat x}_l ] \bm {\hat x}_l (rl,3​−rl,1​)−[(rl,3​−rl,1​)⋅x^l​]x^l​为 r l , 3 r_{l,3} rl,3​在y轴的投影点的向量。 y ^ l \bm {\hat y}_l y^​l​为其单位向量，即指向y轴的单位向量。
令： z ^ l = x ^ l × y ^ l \bm {\hat z}_l=\bm {\hat x}_l×\bm {\hat y}_l z^l​=x^l​×y^​l​为对应z轴的单位向量，得到新坐标系单位矩阵构成的基向量 M l = ∣ x ^ l , y ^ l , z ^ l ∣ \bm M_l=|\bm {\hat x}_l,\bm {\hat y}_l,\bm {\hat z}_l | Ml​=∣x^l​,y^​l​,z^l​∣

同理使用第二组点集构成新的坐标系单位向量 M r = ∣ x ^ r , y ^ r , z ^ r ∣ \bm M_r=|\bm {\hat x}_r,\bm {\hat y}_r,\bm {\hat z}_r| Mr​=∣x^r​,y^​r​,z^r​∣
M l \bm M_l Ml​与 M r \bm M_r Mr​分别为三个点构成的坐标系的基础向量在两个原始坐标系的表示。如下图示意图所示：

因此：
M r = R M l \bm M_r=R\bm M_l Mr​=RMl​
从而可知：
R = M r M l T R=\bm M_r \bm M_l^T R=Mr​MlT​
（5）四元数求解旋转矩阵
使用四元数 q \bm q q表示旋转：
q = q 0 + i q x + j q y + k q z \bm q=q_0+iq_x+jq_y+kq_z q=q0​+iqx​+jqy​+kqz​
r ′ = R r r'=Rr r′=Rr
将空间点表示成思维虚四元数
r ‾ = [ 0 , r ] T = [ 0 , x , y , z ] T = 0 + i x + j y + k z \bm{\overline r}=[0,\bm r]^{\rm T}=[0,x,y,z]^{\rm T}=0+ix+jy+kz r=[0,r]T=[0,x,y,z]T=0+ix+jy+kz
r ′ = q r ‾ q − 1 \bm r'=\bm {q\overline rq}^{-1} r′=qrq−1
其中， q − 1 = q ∗ / q 2 \bm q^{-1}=\bm q^*/\bm q^2 q−1=q∗/q2, q ∗ \bm q^* q∗为 q q q的共轭四元数， q ∗ = q 0 − i q x − j q y − k q z \bm q^*=q_0-iq_x-jq_y-kq_z q∗=q0​−iqx​−jqy​−kqz​
文中使用单位四元数，因此 q 2 = q ⋅ q = 1 q^2=q\cdot q=1 q2=q⋅q=1，因此， q − 1 = q ∗ \bm q^{-1}=\bm q^* q−1=q∗

定义四元数的两种乘法
1）直接

2）点乘

要求解
max ⁡ R ∑ i = 1 n ( r r , i ′ T R r l , i ′ ) \mathop{\max}\limits_R \sum_{i=1}^n(\bm r_{r,i}'^T R\bm r'_{l,i}) Rmax​i=1∑n​(rr,i′T​Rrl,i′​)
使用四元数表示：
max ⁡ q ∑ i = 1 n ( r ‾ r , i ′ ( q r ‾ l , i q − 1 ) ) \mathop{\max}\limits_{\bm q} \sum_{i=1}^n(\bm{\overline r}'_{r,i} (\bm{q\overline r}_{l,i}\bm q^{-1} )) qmax​i=1∑n​(rr,i′​(qrl,i​q−1))
向量乘法交换律：

由四元数的性质:

上式可化为：

其中： r ‾ = [ 0 , r ] T = [ 0 , x , y , z ] T \bm{\overline r}=[0,\bm r]^{\rm T}=[0,x,y,z]^{\rm T} r=[0,r]T=[0,x,y,z]T


以此得：

令:

上式的解：q为N对应最大特征值对应的特征向量，证明过程见参考论文附文A3
（6）总结求解流程：
1）整理残差方程

2）求解尺度因子

3）求解 R R R阵对应四元数表示 q \bm q q

求解其最大特征值对应的特征向量
4）求解偏移量
r 0 = r ‾ r − s R r ‾ l r_0=\bm{\overline r}_r-sR\bm{\overline r}_l r0​=rr​−sRrl​

Marta的迭代最小二乘解法

采用最小二乘迭代求解轨迹对齐变换矩阵，除考虑平移，也考虑了姿态角构成误差方程的约束：
设里程计估计得到的位姿组为 P 1 , P 2 , ⋯ P N {P_1,P_2,⋯P_N} P1​,P2​,⋯PN​，真实轨迹对应的 Q 1 , Q 2 , ⋯ Q N {Q_1,Q_2,⋯Q_N} Q1​,Q2​,⋯QN​，定义两组位姿 P i P_i Pi​和 G i G_i Gi​之间的误差为：
E i = Q i − 1 T P i E_i =Q_i^{-1} TP_i Ei​=Qi−1​TPi​
每个位姿由3轴旋转 R R R及3轴平移构成的4*4矩阵，该矩阵符合李群性质。
对误差进行log-v操作，将误差变换到切空间：该操作与将旋转矩阵转换为李代数(轴角)一致。

定义两组轨迹的总误差约束方程为：

使用LM方法迭代求解 T T T阵

注：