图优化相关知识
1.位姿图
1.1 概念
- 用一个图(Graph)来表示SLAM问题
- 图中的节点来表示机器人的位姿 ( x , y , y a w ) (x,y,yaw) (x,y,yaw)
- 两个节点之间的边表示两个位姿的空间约束(相对位姿关系以及对应方差)
- 一旦形成回环即可进行优化消除误差
- 里程积分的相对位姿视为预测值
- 回环计算的相对位姿视为观测值
- Graph-based SLAM:构建图并调整各节点的位姿,让预测与观测的误差最小
1.2 构建
帧间边
-
里程计测量:
-
相邻节点之间的相对位姿关系,可以由里程计、IMU、帧间匹配计算得到
回环边
-
通过回环检测得到:
-
节点 i i i和节点 j j j在空间上相邻(观测到同样的数据),但是时间上不相邻
-
用帧间匹配算法计算一个相对位姿
2.回环检测方法
一种简单的回环检测方法
⇓
\Darr
⇓
图示对应步骤
- 把节点分为active和inactive两部分
- 跟当前节点(红色节点)时间相近的节点称为active node(黄色节点),其他的称为inactive node
⇓ \Darr ⇓ - 找到当前节点周围一定范围内所有inactive节点,作为回环候选帧(绿色节点)
- 当前节点和回环候选帧进行匹配,根据得分判断是否形成回环
3.非线性最小二乘原理
3.1 要解决的问题
给定一个系统,其状态方程为: f ( x ) = z f(x)=z f(x)=z
- x x x表示系统的状态向量—即需要估计的值
- z z z表示系统的观测值,可以通过传感器进行直接观测
- f ( x ) f(x) f(x)表示一个非线性的映射函数,状态向量 x x x可以通过非线性函数 f ( x ) f(x) f(x)映射得到 z z z
目的就是给定该系统的
n
n
n个混有噪声的观测值
(
z
1
,
…
,
z
n
)
(z_1,\dots,z_n)
(z1,…,zn),估计状态向量
x
x
x,使得其经过
f
(
x
)
f(x)
f(x)映射之后的预测值和观测值的误差最小。
原理跟线性最小二乘基本相同,只是状态方程
f
(
x
)
f(x)
f(x)是一个非线性函数。
3.2 示意图
- x x x表示机器人的位置
- f ( x ) f(x) f(x)为观测模型,节点之间相对位姿计算函数
- z z z为帧间匹配或者回环检测计算出来的相对位姿
- 找到最优的 x x x,让预测和观测的误差最小
3.3 误差函数
目标为最小化预测和观测的差,因此误差即为预测和观测的差:
e
i
(
x
)
=
f
i
(
x
)
−
z
i
′
e_i(x)=f_i(x)-z_i'
ei(x)=fi(x)−zi′
假设误差服从高斯分布,即
e
i
(
x
)
∼
N
(
0
,
Ω
i
)
e_i(x)\sim N(0,\Omega_i)
ei(x)∼N(0,Ωi),
Ω
i
\Omega_i
Ωi为对应的信息矩阵。
我们定义误差的联合概率分布为:
G
(
e
i
(
x
)
)
=
∏
i
1
(
2
π
)
D
/
2
∣
Ω
i
∣
1
/
2
exp
[
−
1
2
e
i
(
x
)
T
Ω
i
e
i
(
x
)
]
G(e_i(x))=\prod_i \frac{1}{(2\pi)^{D/2}|\Omega_i|^{1/2}} \exp[-\frac{1}{2}e_i(x)^T\Omega_ie_i(x)]
G(ei(x))=i∏(2π)D/2∣Ωi∣1/21exp[−21ei(x)TΩiei(x)]
最终目标是使得误差尽可能趋近于0(均值),等价于每个高斯分布取得最大值。
因此误差的联合概率分布
G
(
e
i
(
x
)
)
G(e_i(x))
G(ei(x))取得最大值。
对
G
(
e
i
(
x
)
)
G(e_i(x))
G(ei(x))取对数:
ln
(
G
(
e
i
(
x
)
)
)
=
∑
1
(
2
π
)
D
/
2
∣
Ω
i
∣
1
/
2
−
1
2
∑
e
i
(
x
)
T
Ω
i
e
i
(
x
)
\ln(G(e_i(x)))=\sum \frac{1}{(2\pi)^{D/2}|\Omega_i|^{1/2}} - \frac{1}{2} \sum e_i(x)^T\Omega_ie_i(x)
ln(G(ei(x)))=∑(2π)D/2∣Ωi∣1/21−21∑ei(x)TΩiei(x)
即若让
G
(
e
i
(
x
)
)
G(e_i(x))
G(ei(x))取最大值,就是让
∑
e
i
(
x
)
T
Ω
i
e
i
(
x
)
\sum e_i(x)^T\Omega_ie_i(x)
∑ei(x)TΩiei(x)取最小值。
令非线性最小二乘的目标函数为:
min
F
(
x
)
=
min
∑
e
i
(
x
)
T
Ω
i
e
i
(
x
)
\min F(x)=\min \sum e_i(x)^T\Omega_ie_i(x)
minF(x)=min∑ei(x)TΩiei(x)
3.4 求解
目标函数:
min
x
F
(
x
)
\min_x F(x)
xminF(x)
直接想法:
求
F
(
x
)
F(x)
F(x)关于变量
x
x
x的导数,令其等于0,求解方程即可。
对于凸函数来说,上述想法是可行的,但对于非凸函数,通常采用基于梯度的优化方法。
3.4.1 线性化
F
(
x
)
F(x)
F(x)为关于
x
x
x的非线性方程,将其化为关于
x
x
x的线性方程。
F
(
x
)
=
∑
e
i
(
x
)
T
Ω
i
e
i
(
x
)
F(x)=\sum e_i(x)^T\Omega_ie_i(x)
F(x)=∑ei(x)TΩiei(x)
误差函数
e
i
(
x
)
e_i(x)
ei(x)是非线性函数,因此
F
(
x
)
F(x)
F(x)是关于
x
x
x的非线性函数。对误差函数
e
i
(
x
)
e_i(x)
ei(x)进行线性化得:
e
i
(
x
+
Δ
x
)
=
e
i
(
x
)
+
J
i
(
x
)
Δ
x
e_i(x+\Delta x)=e_i(x)+J_i(x)\Delta x
ei(x+Δx)=ei(x)+Ji(x)Δx
其中,
J
J
J为映射函数
F
(
x
)
F(x)
F(x)对状态向量
x
x
x的导数,称之为Jacobian矩阵。
J
i
(
x
)
=
(
∂
f
i
(
x
)
∂
x
1
,
∂
f
i
(
x
)
∂
x
2
,
…
,
∂
f
i
(
x
)
∂
x
n
)
J_i(x)=(\frac{\partial f_i(x)}{\partial x_1},\frac{\partial f_i(x)}{\partial x_2},\dots,\frac{\partial f_i(x)}{\partial x_n})
Ji(x)=(∂x1∂fi(x),∂x2∂fi(x),…,∂xn∂fi(x))
因此,函数
F
(
x
)
F(x)
F(x)可线性化为:
F
(
x
+
Δ
x
)
=
∑
e
i
T
(
x
+
Δ
x
)
Ω
i
e
i
(
x
+
Δ
x
)
=
∑
(
e
i
(
x
)
+
J
i
Δ
x
)
T
Ω
i
(
e
i
(
x
)
+
J
i
Δ
x
)
=
∑
(
e
i
T
Ω
i
e
i
+
e
i
T
Ω
i
J
i
Δ
x
+
Δ
x
T
J
i
T
Ω
i
e
i
+
Δ
x
T
J
i
T
Ω
i
J
i
Δ
x
)
=
∑
(
e
i
T
Ω
i
e
i
+
2
e
i
T
Ω
i
J
i
Δ
x
+
Δ
x
T
J
i
T
Ω
i
J
i
Δ
x
)
=
∑
c
i
+
∑
(
2
b
i
T
Δ
x
+
Δ
x
T
H
i
Δ
x
)
\begin{aligned} F(x+\Delta x)&=\sum e_i^T(x+\Delta x)\Omega_ie_i(x+\Delta x) \\ &=\sum (e_i(x)+J_i\Delta x)^T\Omega_i(e_i(x)+J_i\Delta x) \\ &=\sum (e_i^T\Omega_ie_i+e_i^T\Omega_iJ_i\Delta x+\Delta x^TJ_i^T\Omega_ie_i+\Delta x^TJ_i^T\Omega_iJ_i\Delta x) \\ &=\sum (e_i^T\Omega_ie_i+2e_i^T\Omega_iJ_i\Delta x+\Delta x^TJ_i^T\Omega_iJ_i\Delta x) \\ &=\sum c_i + \sum (2b_i^T\Delta x+\Delta x^T H_i \Delta x) \end{aligned}
F(x+Δx)=∑eiT(x+Δx)Ωiei(x+Δx)=∑(ei(x)+JiΔx)TΩi(ei(x)+JiΔx)=∑(eiTΩiei+eiTΩiJiΔx+ΔxTJiTΩiei+ΔxTJiTΩiJiΔx)=∑(eiTΩiei+2eiTΩiJiΔx+ΔxTJiTΩiJiΔx)=∑ci+∑(2biTΔx+ΔxTHiΔx)
其中,
b
i
T
=
e
i
T
Ω
i
J
i
,
H
i
=
J
i
T
Ω
i
J
i
b_i^T=e_i^T\Omega_iJ_i,H_i=J_i^T\Omega_iJ_i
biT=eiTΩiJi,Hi=JiTΩiJi。
F
(
x
+
Δ
x
)
=
∑
c
i
+
∑
(
2
b
i
T
Δ
x
+
Δ
x
T
H
i
Δ
x
)
=
∑
c
i
+
∑
2
b
i
T
Δ
x
+
Δ
x
T
∑
H
i
Δ
x
=
∑
c
i
+
2
b
T
Δ
x
+
Δ
x
T
H
Δ
x
\begin{aligned} F(x+\Delta x)&=\sum c_i + \sum (2b_i^T\Delta x+\Delta x^T H_i \Delta x) \\ &=\sum c_i + \sum 2b_i^T\Delta x+\Delta x^T \sum H_i \Delta x \\ &=\sum c_i + 2b^T\Delta x+\Delta x^T H \Delta x \end{aligned}
F(x+Δx)=∑ci+∑(2biTΔx+ΔxTHiΔx)=∑ci+∑2biTΔx+ΔxT∑HiΔx=∑ci+2bTΔx+ΔxTHΔx
F
(
x
+
Δ
x
)
F(x+\Delta x)
F(x+Δx)为关于变量
Δ
x
\Delta x
Δx的二次函数,令其关于
Δ
x
\Delta x
Δx的导数等于0,可求解得到
F
(
x
+
Δ
x
)
F(x+\Delta x)
F(x+Δx)的极值,即
∂
F
(
x
+
Δ
x
)
∂
Δ
x
=
2
b
+
2
H
Δ
x
\frac{\partial F(x+\Delta x)}{\partial \Delta x}=2b+2H\Delta x
∂Δx∂F(x+Δx)=2b+2HΔx
令
∂
F
(
x
+
Δ
x
)
∂
Δ
x
=
0
\frac{\partial F(x+\Delta x)}{\partial \Delta x}=0
∂Δx∂F(x+Δx)=0,得
H
Δ
x
=
−
b
⟹
Δ
x
∗
=
−
H
−
1
b
H\Delta x=-b \implies \Delta x^*=-H^{-1}b
HΔx=−b⟹Δx∗=−H−1b
令
x
=
x
+
Δ
x
∗
x=x+\Delta x^*
x=x+Δx∗,然后不断迭代,直至收敛即可。
3.5 流程总结
- 线性化误差函数:
e i ( x + Δ x ) = e i ( x ) + J i ( x ) Δ x e_i(x+\Delta x)=e_i(x)+J_i(x)\Delta x ei(x+Δx)=ei(x)+Ji(x)Δx - 构建线性系统:
b T = ∑ e i T Ω i J i H = ∑ J i T Ω i J i H Δ x = b \begin{aligned} &b^T=\sum e_i^T\Omega_iJ_i \\ &H=\sum J_i^T\Omega_iJ_i \\ &H\Delta x=b \end{aligned} bT=∑eiTΩiJiH=∑JiTΩiJiHΔx=b - 求解线性系统:
Δ x ∗ = − H − 1 b \Delta x^*=-H^{-1}b Δx∗=−H−1b - 更新解,并不断迭代直至收敛:
x = x + Δ x ∗ x=x+\Delta x^* x=x+Δx∗
4.非线性最小二乘原理在SLAM中的应用
4.1 误差函数
4.1.1 定义
- 观测值为匹配计算得到的节点
i
i
i和节点
j
j
j的相对位姿
z i j ′ = ( t i j , θ i j ) Z i j ′ = V 2 T ( z i j ′ ) \begin{aligned} z_{ij}'&=(t_{ij},\theta_{ij}) \\ Z_{ij}'&=V2T(z_{ij}') \end{aligned} zij′Zij′=(tij,θij)=V2T(zij′) - 预测值为里程积分得到的当前节点
i
i
i和节点
j
j
j的相对位姿
Z i j = f ( x i , x j ) = X i − 1 X j X i = V 2 T ( x i ) X j = V 2 T ( x j ) \begin{aligned} Z_{ij}&=f(x_i,x_j)=X_i^{-1}X_j \\ X_i&=V2T(x_i) \\ X_j&=V2T(x_j) \end{aligned} ZijXiXj=f(xi,xj)=Xi−1Xj=V2T(xi)=V2T(xj)
已知:
X i = [ R i t i 0 1 ] ⟹ X i − 1 = [ R i T − R i T t i 0 1 ] X j = [ R j t j 0 1 ] \begin{aligned} X_i&= \begin{bmatrix} R_i & t_i \\ 0 & 1 \end{bmatrix} \implies X_i^{-1}= \begin{bmatrix} R_i^T & -R_i^Tt_i \\ 0 & 1 \end{bmatrix} \\ X_j&= \begin{bmatrix} R_j & t_j \\ 0 & 1 \end{bmatrix} \end{aligned} XiXj=[Ri0ti1]⟹Xi−1=[RiT0−RiTti1]=[Rj0tj1]
则预测值:
Z i j = X i − 1 X j = [ R i T − R i T t i 0 1 ] [ R j t j 0 1 ] = [ R i T R j R i T ( t j − t i ) 0 1 ] z i j = T 2 V ( Z i j ) = [ R i T ( t j − t i ) θ j − θ i ] \begin{aligned} Z_{ij}&=X_i^{-1}X_j= \begin{bmatrix} R_i^T & -R_i^Tt_i \\ 0 & 1 \end{bmatrix} \begin{bmatrix} R_j & t_j \\ 0 & 1 \end{bmatrix}= \begin{bmatrix} R_i^TR_j & R_i^T(t_j-t_i) \\ 0 & 1 \end{bmatrix} \\ z_{ij}&=T2V(Z_{ij})= \begin{bmatrix} R_i^T(t_j-t_i) \\ \theta_j-\theta_i \end{bmatrix} \end{aligned} Zijzij=Xi−1Xj=[RiT0−RiTti1][Rj0tj1]=[RiTRj0RiT(tj−ti)1]=T2V(Zij)=[RiT(tj−ti)θj−θi] - 误差函数的定义
e i j ( x ) = T 2 V ( Z i j ′ − 1 Z i j ) e_{ij}(x)=T2V(Z_{ij}'^{-1}Z_{ij}) eij(x)=T2V(Zij′−1Zij) - 误差函数的矩阵形式
e i j ( x ) = [ R i j T ( R i T ( t j − t i ) − t i j ) θ j − θ i − θ i j ] e_{ij}(x)= \begin{bmatrix} R_{ij}^T(R_i^T(t_j-t_i)-t_{ij}) \\ \theta_j-\theta_i-\theta_{ij} \end{bmatrix} eij(x)=[RijT(RiT(tj−ti)−tij)θj−θi−θij] - 对应的Jacobian矩阵
∂ e i j ( x ) ∂ x i = [ − R i j T R i T R i j T ∂ R i T ∂ θ ( t j − t i ) 0 − 1 ] ∂ e i j ( x ) ∂ x j = [ R i j T R i T 0 0 1 ] \begin{aligned} \frac{\partial e_{ij}(x)}{\partial x_i}&= \begin{bmatrix} -R_{ij}^TR_i^T & R_{ij}^T \frac{\partial R_i^T}{\partial \theta}(t_j-t_i)\\ 0 & -1 \end{bmatrix} \\ \frac{\partial e_{ij}(x)}{\partial x_j}&= \begin{bmatrix} R_{ij}^TR_i^T & 0 \\ 0 & 1 \end{bmatrix} \end{aligned} ∂xi∂eij(x)∂xj∂eij(x)=[−RijTRiT0RijT∂θ∂RiT(tj−ti)−1]=[RijTRiT001]
4.1.2 线性化
误差函数:
e
i
j
(
x
+
Δ
x
)
=
e
i
j
(
x
)
+
J
i
j
(
x
)
Δ
x
J
i
j
=
∂
e
i
j
(
x
)
∂
x
e_{ij}(x+\Delta x)=e_{ij}(x)+J_{ij}(x)\Delta x \\ J_{ij}=\frac{\partial e_{ij}(x)}{\partial x}
eij(x+Δx)=eij(x)+Jij(x)ΔxJij=∂x∂eij(x)
因为误差函数只跟
x
i
x_i
xi和
x
j
x_j
xj有关,因此具有下列性质:
∂
e
i
j
(
x
)
∂
x
=
(
0
,
…
,
∂
e
i
j
(
x
)
∂
x
i
,
…
,
∂
e
i
j
(
x
)
∂
x
j
,
…
,
0
)
J
i
j
=
(
0
,
…
,
A
i
j
,
…
,
B
i
j
,
…
,
0
)
\begin{aligned} \frac{\partial e_{ij}(x)}{\partial x}&=(0,\dots,\frac{\partial e_{ij}(x)}{\partial x_i},\dots,\frac{\partial e_{ij}(x)}{\partial x_j},\dots,0) \\ J_{ij}&=(0,\dots,A_{ij},\dots,B_{ij},\dots,0) \end{aligned}
∂x∂eij(x)Jij=(0,…,∂xi∂eij(x),…,∂xj∂eij(x),…,0)=(0,…,Aij,…,Bij,…,0)
4.2 固定坐标系
- 观测值观测到的两个位姿之间的相对位姿
- 满足相对位姿约束的解有无穷多组
- 为了让解唯一,必须加入一个约束条件让某一个位姿固定,一般选择第一个位姿,即:
Δ x 1 = 0 \Delta x_1=0 Δx1=0
等价于:
[ 1 0 0 0 1 0 0 0 1 ] Δ x 1 = [ 0 0 0 ] \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \Delta x_1= \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix} 100010001 Δx1= 000 - 加入的约束为:
[ 1 0 0 0 1 0 0 0 1 ] Δ x 1 = [ 0 0 0 ] \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \Delta x_1= \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix} 100010001 Δx1= 000 - 求解的线性系统为:
H Δ x = − b H\Delta x=-b HΔx=−b - 因此等价于:
H 11 + = [ 1 0 0 0 1 0 0 0 1 ] H_{11}+= \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} H11+= 100010001
4.3 构建线性系统
- 已知误差项和Jacobian矩阵 A i j A_{ij} Aij和 B i j B_{ij} Bij
- 向量
b
b
b的更新为:
b i T + = e i j T Ω i j A i j b j T + = e i j T Ω i j B i j \begin{aligned} b_i^T&+=e_{ij}^T\Omega_{ij}A_{ij} \\ b_j^T&+=e_{ij}^T\Omega_{ij}B_{ij} \\ \end{aligned} biTbjT+=eijTΩijAij+=eijTΩijBij - 矩阵
H
H
H的更新为:
H i i + = A i j T Ω i j A i j H i j + = A i j T Ω i j B i j H j i + = B i j T Ω i j A i j H j j + = B i j T Ω i j B i j \begin{aligned} H_{ii}&+=A_{ij}^T\Omega_{ij}A_{ij} \\ H_{ij}&+=A_{ij}^T\Omega_{ij}B_{ij} \\ H_{ji}&+=B_{ij}^T\Omega_{ij}A_{ij} \\ H_{jj}&+=B_{ij}^T\Omega_{ij}B_{ij} \\ \end{aligned} HiiHijHjiHjj+=AijTΩijAij+=AijTΩijBij+=BijTΩijAij+=BijTΩijBij
4.4 求解
- 已知矩阵 H H H和向量 b b b
- 求解线性方程组:
Δ x = − H − 1 b \Delta x=-H^{-1}b Δx=−H−1b - 不断进行迭代,直至收敛:
x = x + Δ x x=x+\Delta x x=x+Δx