多传感器融合定位十一-基于滤波的融合方法Ⅱ
- 1. 编码器运动模型及标定
- 1.1 编码器基础知识
- 1.2 编码器运动模型
- 1.2.1 旋转半径求解
- 1.2.2 角速度求解
- 1.2.3 线速度求解
- 1.2.4 位姿求解
- 1.3 编码器的标定
- 1.3.1 轮子半径标定
- 1.3.2 轮子与底盘中心距离标定
- 2. 融合编码器的滤波方法
- 2.1 核心思路
- 2.2 观测量定义
- 2.3 观测方程推导
- 3. 融合运动约束的滤波方法
- 4. 融合点云特征的滤波方法
- 4.1 整体思路
- 4.2 滤波模型
- 4.3 位姿更新
- 4.4 相似工作
Reference:
- 深蓝学院-多传感器融合
- 多传感器融合定位理论基础
文章跳转:
- 多传感器融合定位一-3D激光里程计其一:ICP
- 多传感器融合定位二-3D激光里程计其二:NDT
- 多传感器融合定位三-3D激光里程计其三:点云畸变补偿
- 多传感器融合定位四-3D激光里程计其四:点云线面特征提取
- 多传感器融合定位五-点云地图构建及定位
- 多传感器融合定位六-惯性导航原理及误差分析
- 多传感器融合定位七-惯性导航解算及误差分析其一
- 多传感器融合定位八-惯性导航解算及误差分析其二
- 多传感器融合定位九-基于滤波的融合方法Ⅰ其一
- 多传感器融合定位十-基于滤波的融合方法Ⅰ其二
- 多传感器融合定位十一-基于滤波的融合方法Ⅱ
- 多传感器融合定位十二-基于图优化的建图方法其一
- 多传感器融合定位十三-基于图优化的建图方法其二
- 多传感器融合定位十四-基于图优化的定位方法
- 多传感器融合定位十五-多传感器时空标定(综述)
1. 编码器运动模型及标定
1.1 编码器基础知识
编码器
感应轮子的旋转,并在旋转时输出脉冲,脉冲数与转过的角度呈线性比例关系。
脉冲对应的是角度增量,有时也用增量除以时间,转成轮子转动的角速度输出。
需要注意的是,编码器只是各种转角测量方式中的一种,其他还有轮速计、霍尔传感器等,本课程以编码器为例子讲解模型,但同样适用于其他形式的传感器。
编码器安装方式有单轮、双轮、三轮,本课程推导仅围绕双轮差分模型进行展开。
该模型中,需要用到以下变量:
- r L \boldsymbol{r}_L rL:左轮半径
- r R \boldsymbol{r}_R rR:右轮半径
- d \boldsymbol{d} d:轮子离底盘中心的距离
- ω L \boldsymbol{\omega}_L ωL:左轮自转角速度
- ω R \boldsymbol{\omega_R} ωR:右轮自转角速度
- v L \boldsymbol{v}_L vL:左轮线速度
- v R \boldsymbol{v}_R vR:右轮线速度
实际使用时,标定完成后,
r
L
\boldsymbol{r}_L
rL、
r
R
\boldsymbol{r}_R
rR、
d
\boldsymbol{d}
d 为固定参数,
ω
L
\boldsymbol{\omega}_L
ωL 和
ω
R
\boldsymbol{\omega}_R
ωR 为测量值,而
v
L
\boldsymbol{v}_L
vL 和
v
R
\boldsymbol{v}_R
vR 可以通过下式计算得到:
v
L
=
ω
L
r
L
v
R
=
ω
R
r
R
\begin{aligned} \boldsymbol{v}_L & =\boldsymbol{\omega}_L \boldsymbol{r}_L \\ \boldsymbol{v}_R & =\boldsymbol{\omega}_R \boldsymbol{r}_R \end{aligned}
vLvR=ωLrL=ωRrR
1.2 编码器运动模型
运动模型的作用是,使用前述已知量,求解以下变量:
ω
\boldsymbol{\omega}
ω : 底盘中心的角速度
v
\boldsymbol{v}
v : 底盘中心的线速度
r
\boldsymbol{r}
r : 底盘中心圆弧运动旋转半径
1.2.1 旋转半径求解
双轮差分模型下,左右轮圆弧运动的角速度相等,且等于底盘中心圆弧运动的角速度(两个轮子的自转角速度是相同的),因此有:
ω
=
v
L
r
−
d
=
v
R
r
+
d
\boldsymbol{\omega}=\frac{\boldsymbol{v}_L}{\boldsymbol{r}-\boldsymbol{d}}=\frac{\boldsymbol{v}_R}{\boldsymbol{r}+\boldsymbol{d}}
ω=r−dvL=r+dvR由此可以得出:
v
L
(
r
+
d
)
=
v
R
(
r
−
d
)
\boldsymbol{v}_L(\boldsymbol{r}+\boldsymbol{d})=\boldsymbol{v}_R(\boldsymbol{r}-\boldsymbol{d})
vL(r+d)=vR(r−d)移项可得:
(
v
R
−
v
L
)
r
=
(
v
R
+
v
L
)
d
\left(\boldsymbol{v}_R-\boldsymbol{v}_L\right) \boldsymbol{r}=\left(\boldsymbol{v}_R+\boldsymbol{v}_L\right) \boldsymbol{d}
(vR−vL)r=(vR+vL)d从而可以得到:
r
=
v
R
+
v
L
v
R
−
v
L
d
\boldsymbol{r}=\frac{\boldsymbol{v}_R+\boldsymbol{v}_L}{\boldsymbol{v}_R-\boldsymbol{v}_L} \boldsymbol{d}
r=vR−vLvR+vLd
1.2.2 角速度求解
把旋转半径的求解结果,代入角速度公式,即可得到:
ω
=
v
L
v
R
+
v
L
v
R
−
v
L
d
−
d
=
v
R
−
v
L
2
d
\boldsymbol{\omega}=\frac{\boldsymbol{v}_L}{\frac{\boldsymbol{v}_R+\boldsymbol{v}_L}{\boldsymbol{v}_R-\boldsymbol{v}_L} \boldsymbol{d}-\boldsymbol{d}}=\frac{\boldsymbol{v}_R-\boldsymbol{v}_L}{2 \boldsymbol{d}}
ω=vR−vLvR+vLd−dvL=2dvR−vL
1.2.3 线速度求解
利用旋转角速度和旋转半径的结果,可以直接得到线速度:
v
=
ω
r
=
v
R
−
v
L
2
d
v
R
+
v
L
v
R
−
v
L
d
=
v
R
+
v
L
2
\boldsymbol{v}=\boldsymbol{\omega} \boldsymbol{r}=\frac{\boldsymbol{v}_R-\boldsymbol{v}_L}{2 \boldsymbol{d}} \frac{\boldsymbol{v}_R+\boldsymbol{v}_L}{\boldsymbol{v}_R-\boldsymbol{v}_L} \boldsymbol{d}=\frac{\boldsymbol{v}_R+\boldsymbol{v}_L}{2}
v=ωr=2dvR−vLvR−vLvR+vLd=2vR+vL
1.2.4 位姿求解
假设
x
k
,
y
k
,
θ
k
\boldsymbol{x}_k, \boldsymbol{y}_k, \boldsymbol{\theta}_k
xk,yk,θk 为当前时刻位姿,
x
k
−
1
,
y
k
−
1
,
θ
k
−
1
\boldsymbol{x}_{k-1}, \boldsymbol{y}_{k-1}, \boldsymbol{\theta}_{k-1}
xk−1,yk−1,θk−1 为上一时刻的位姿,则有:
θ
k
=
θ
k
−
1
+
ω
Δ
t
x
k
=
x
k
−
1
+
v
Δ
t
cos
(
θ
k
−
1
)
y
k
=
y
k
−
1
+
v
Δ
t
sin
(
θ
k
−
1
)
\begin{aligned} & \boldsymbol{\theta}_k=\boldsymbol{\theta}_{k-1}+\boldsymbol{\omega} \Delta t \\ & \boldsymbol{x}_k=\boldsymbol{x}_{k-1}+\boldsymbol{v} \Delta t \cos \left(\boldsymbol{\theta}_{\boldsymbol{k - 1}}\right) \\ & \boldsymbol{y}_k=\boldsymbol{y}_{k-1}+\boldsymbol{v} \Delta t \sin \left(\boldsymbol{\theta}_{\boldsymbol{k}-1}\right) \end{aligned}
θk=θk−1+ωΔtxk=xk−1+vΔtcos(θk−1)yk=yk−1+vΔtsin(θk−1)其中:
Δ
t
=
t
k
−
t
k
−
1
\Delta t=t_k-t_{k-1}
Δt=tk−tk−1
1.3 编码器的标定
标定可以理解为运动模型求解过程的反向过程,具体是指在已知底盘中心线速度、角速度的情况下,求解轮子半径、轮子离底盘中心距离等。
已知量:
v
v
v:底盘中心的线速度
ω
\omega
ω:底盘中心的角速度
待求解量:
r
L
\boldsymbol{r}_L
rL:左轮半径
r
R
\boldsymbol{r}_R
rR:右轮半径
d
\boldsymbol{d}
d:轮子离底盘中心的距离
实际标定时,线速度、角速度由其他传感器提供(比如雷达点云和地图匹配),且为了简化模型,认为雷达装在底盘中心正上方。(这里雷达最好的方法是先建好一个点云地图,然后在点云地图里面匹配,然后拿它做观测,去提供线速度和角速度的结果,而不要用激光里程计-----里程计本身就是有累计误差的)
1.3.1 轮子半径标定
由于速度的求解公式为:
v
=
v
R
+
v
L
2
=
ω
R
r
R
+
ω
L
r
L
2
\boldsymbol{v}=\frac{\boldsymbol{v}_R+\boldsymbol{v}_L}{2}=\frac{\boldsymbol{\omega}_R \boldsymbol{r}_R+\boldsymbol{\omega}_L \boldsymbol{r}_L}{2}
v=2vR+vL=2ωRrR+ωLrL它可以重新写为:
[
ω
R
ω
L
]
[
r
R
r
L
]
=
2
v
\left[\begin{array}{ll} \boldsymbol{\omega}_R & \boldsymbol{\omega}_L \end{array}\right]\left[\begin{array}{l} \boldsymbol{r}_R \\ \boldsymbol{r}_L \end{array}\right]=2 \boldsymbol{v}
[ωRωL][rRrL]=2v当有多组测量值时,可以构成如下方程组:
[
ω
R
0
ω
L
0
ω
R
1
ω
L
1
⋮
⋮
ω
R
N
ω
L
N
]
[
r
R
r
L
]
=
[
2
v
0
2
v
1
⋮
2
v
N
]
\left[\begin{array}{cc} \boldsymbol{\omega}_{R 0} & \boldsymbol{\omega}_{L 0} \\ \boldsymbol{\omega}_{R 1} & \boldsymbol{\omega}_{L 1} \\ \vdots & \vdots \\ \boldsymbol{\omega}_{R N} & \boldsymbol{\omega}_{L N} \end{array}\right]\left[\begin{array}{l} \boldsymbol{r}_R \\ \boldsymbol{r}_L \end{array}\right]=\left[\begin{array}{c} 2 \boldsymbol{v}_0 \\ 2 \boldsymbol{v}_1 \\ \vdots \\ 2 \boldsymbol{v}_N \end{array}\right]
ωR0ωR1⋮ωRNωL0ωL1⋮ωLN
[rRrL]=
2v02v1⋮2vN
这是典型的最小二乘问题,可用最小二乘标准形式计算。
1.3.2 轮子与底盘中心距离标定
由于角速度的求解公式为:
ω
=
v
R
−
v
L
2
d
\boldsymbol{\omega}=\frac{\boldsymbol{v}_R-\boldsymbol{v_L}}{2 \boldsymbol{d}}
ω=2dvR−vL在经过轮子半径标定之后,分子上的两项可认为是已知量,因此可以得到:
d
=
v
R
−
v
L
2
ω
\boldsymbol{d}=\frac{\boldsymbol{v}_R-\boldsymbol{v}_L}{2 \boldsymbol{\omega}}
d=2ωvR−vL虽然可直接求解,但是为了抑制噪声带来的影响,因多次采样计算取平均。
2. 融合编码器的滤波方法
2.1 核心思路
在上一节课滤波模型的基础上增加编码器进行融合,有一种非常简单的方法,即使用编码器解算的速度作为观测量,加入原来模型的观测方程中,而其他环节保持不变。
2.2 观测量定义
编码器提供的是载体系下的速度观测,在前
(
x
)
(\mathrm{x})
(x)-左
(
y
)
(\mathrm{y})
(y)-上(z)坐标系的定义下,
x
\mathrm{x}
x 方向的速度分量是已知的
v
x
b
=
v
m
\boldsymbol{v}_x^b=\boldsymbol{v}_m
vxb=vm。 另外,在以车作为载体的情况下,由于车的侧向和天向没有运动,因此又有
v
y
b
=
0
\boldsymbol{v}_y^b=0
vyb=0,
v
z
b
=
0
\boldsymbol{v}_z^b=0
vzb=0。
基于此,我们可以认为
b
b
b 系 3 个维度的速度分量都是可观测的。
2.3 观测方程推导
由于导航解算得到的是
w
w
w 系下得速度,而速度观测是
b
b
b 系下得,因此需要推导二者之间的误差关系,才能得到相应的观测方程。
推导方法仍按照第6讲的固定套路进行。
- 写出不考虑误差时的方程:
v b = R b w v w \boldsymbol{v}^b=\boldsymbol{R}_{b w} \boldsymbol{v}^w vb=Rbwvw - 写出考虑误差时的方程:
v ~ b = R ~ b w v ~ w \tilde{\boldsymbol{v}}^b=\tilde{\boldsymbol{R}}_{b w} \tilde{\boldsymbol{v}}^w v~b=R~bwv~w - 写出真实值与理想值之间的关系:
v ~ b = v b + δ v b v ~ w = v w + δ v w R ~ b w = R ~ w b T = ( R w b ( I + [ δ θ ] × ) ) T = ( I − [ δ θ ] × ) R b w \begin{aligned} & \tilde{\boldsymbol{v}}^b=\boldsymbol{v}^b+\delta \boldsymbol{v}^b \\ & \tilde{\boldsymbol{v}}^w=\boldsymbol{v}^w+\delta \boldsymbol{v}^w \\ & \tilde{\boldsymbol{R}}_{b w}=\tilde{\boldsymbol{R}}_{w b}^T=\left(\boldsymbol{R}_{w b}\left(\boldsymbol{I}+[\delta \boldsymbol{\theta}]_{\times}\right)\right)^T \\ & =\left(\boldsymbol{I}-[\delta \boldsymbol{\theta}]_{\times}\right) \boldsymbol{R}_{b w} \end{aligned} v~b=vb+δvbv~w=vw+δvwR~bw=R~wbT=(Rwb(I+[δθ]×))T=(I−[δθ]×)Rbw - 把3)中的关系带入2)式:
v b + δ v b = ( I − [ δ θ ] × ) R b w ( v w + δ v w ) \boldsymbol{v}^b+\delta \boldsymbol{v}^b=\left(\boldsymbol{I}-[\delta \boldsymbol{\theta}]_{\times}\right) \boldsymbol{R}_{b w}\left(\boldsymbol{v}^w+\delta \boldsymbol{v}^w\right) vb+δvb=(I−[δθ]×)Rbw(vw+δvw) - 把1)中的关系带入4)式:
R b w v w + δ v b = ( I − [ δ θ ] × ) R b w ( v w + δ v w ) \boldsymbol{R}_{b w} \boldsymbol{v}^w+\delta \boldsymbol{v}^b=\left(\boldsymbol{I}-[\delta \boldsymbol{\theta}]_{\times}\right) \boldsymbol{R}_{b w}\left(\boldsymbol{v}^w+\delta \boldsymbol{v}^w\right) Rbwvw+δvb=(I−[δθ]×)Rbw(vw+δvw) - 化简方程:
δ v b = R b w δ v w − [ δ θ ] × R b w v w = R b w δ v w − [ δ θ ] × v b = R b w δ v w + [ v b ] × δ θ \begin{aligned} \delta \boldsymbol{v}^b & =\boldsymbol{R}_{b w} \delta \boldsymbol{v}^w-[\delta \boldsymbol{\theta}]_{\times} \boldsymbol{R}_{b w} \boldsymbol{v}^w \\ & =\boldsymbol{R}_{b w} \delta \boldsymbol{v}^w-[\delta \boldsymbol{\theta}]_{\times} \boldsymbol{v}^b \\ & =\boldsymbol{R}_{b w} \delta \boldsymbol{v}^w+\left[\boldsymbol{v}^b\right]_{\times} \delta \boldsymbol{\theta} \end{aligned} δvb=Rbwδvw−[δθ]×Rbwvw=Rbwδvw−[δθ]×vb=Rbwδvw+[vb]×δθ
根据前一章内容,状态量为
δ
x
=
[
δ
p
δ
v
δ
θ
δ
b
a
δ
b
ω
]
\delta \boldsymbol{x}=\left[\begin{array}{c} \delta \boldsymbol{p} \\ \delta \boldsymbol{v} \\ \delta \boldsymbol{\theta} \\ \delta \boldsymbol{b}_a \\ \delta \boldsymbol{b}_\omega \end{array}\right]
δx=
δpδvδθδbaδbω
而融合编码器以后,观测量变为
y
=
[
δ
p
‾
δ
v
‾
b
δ
θ
‾
]
\boldsymbol{y}=\left[\begin{array}{c} \delta \overline{\boldsymbol{p}} \\ \delta \overline{\boldsymbol{v}}^b \\ \delta \overline{\boldsymbol{\theta}} \end{array}\right]
y=
δpδvbδθ
其中
δ
v
‾
b
\delta \overline{\boldsymbol{v}}^b
δvb 的观测值可以通过下式获得
δ
v
‾
b
=
v
~
b
−
v
b
=
R
~
b
w
v
~
w
−
[
v
m
0
0
]
\delta \overline{\boldsymbol{v}}_b=\tilde{\boldsymbol{v}}^b-\boldsymbol{v}^b=\tilde{\boldsymbol{R}}_{b w} \tilde{\boldsymbol{v}}^w-\left[\begin{array}{c} \boldsymbol{v}_m \\ 0 \\ 0 \end{array}\right]
δvb=v~b−vb=R~bwv~w−
vm00
此时的观测方程
y
=
G
t
δ
x
+
C
t
n
\boldsymbol{y}=\boldsymbol{G}_t \delta \boldsymbol{x}+\boldsymbol{C}_t \boldsymbol{n}
y=Gtδx+Ctn 中的各变量应重新写为
G
t
=
[
I
3
0
0
0
0
0
R
b
w
[
v
b
]
×
0
0
0
0
I
3
0
0
]
C
t
=
[
I
3
0
0
0
I
3
0
0
0
I
3
]
n
=
[
n
δ
p
ˉ
x
n
δ
p
ˉ
y
n
δ
p
ˉ
z
n
δ
v
ˉ
x
b
n
δ
v
ˉ
y
b
n
δ
v
ˉ
z
b
n
δ
θ
ˉ
x
n
δ
θ
ˉ
y
n
δ
θ
ˉ
z
]
T
\begin{aligned} & \boldsymbol{G}_t=\left[\begin{array}{ccccc} \boldsymbol{I}_3 & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \boldsymbol{R}_{b w} & {\left[\boldsymbol{v}^b\right]_{\times}} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \boldsymbol{I}_3 & \mathbf{0} & \mathbf{0} \end{array}\right] \\ & \boldsymbol{C}_t=\left[\begin{array}{ccc} \boldsymbol{I}_3 & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \boldsymbol{I}_3 & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \boldsymbol{I}_3 \end{array}\right] \\ & \boldsymbol{n}=\left[\begin{array}{lllllllll} n_{\delta \bar{p}_x} & n_{\delta \bar{p}_y} & n_{\delta \bar{p}_z} & n_{\delta \bar{v}_x^b} & n_{\delta \bar{v}_y^b} & n_{\delta \bar{v}_z^b} & n_{\delta \bar{\theta}_x} & n_{\delta \bar{\theta}_y} & n_{\delta \bar{\theta}_z} \end{array}\right]^T \end{aligned}
Gt=
I3000Rbw00[vb]×I3000000
Ct=
I3000I3000I3
n=[nδpˉxnδpˉynδpˉznδvˉxbnδvˉybnδvˉzbnδθˉxnδθˉynδθˉz]T随后,便可以使用新的观测方程,不改变其他方程,直接按照原有流程进行Kalman滤波融合。
3. 融合运动约束的滤波方法
很多时候,硬件平台并没有编码器,不能直接使用上一小节的模型,但是车本身的运动特性(即侧向速度和天向速度为0)仍然可以使用。
它对观测量带来的改变仅仅是少了一个维度( x x x 方向),而推导方法并没有改变,因此此处直接给出该融合模式下的推导结果。
新的观测量为
y
=
[
δ
p
‾
[
δ
v
‾
b
]
y
z
δ
θ
‾
]
\boldsymbol{y}=\left[\begin{array}{c} \delta \overline{\boldsymbol{p}} \\ {\left[\delta \overline{\boldsymbol{v}}^b\right]_{y z}} \\ \delta \overline{\boldsymbol{\theta}} \end{array}\right]
y=
δp[δvb]yzδθ
[
∙
]
y
z
[\bullet]_{y z}
[∙]yz 表示只取三维向量或矩阵的后2行
此时的观测方程
y
=
G
t
δ
x
+
C
t
n
\boldsymbol{y}=\boldsymbol{G}_t \delta \boldsymbol{x}+\boldsymbol{C}_t \boldsymbol{n}
y=Gtδx+Ctn 中的各变量应重新写为
G
t
=
[
I
3
0
0
0
0
0
[
R
b
w
]
y
z
[
[
v
b
]
×
]
y
z
0
0
0
0
I
3
0
0
]
C
t
=
[
I
3
0
0
0
I
2
0
0
0
I
3
]
n
=
[
n
δ
p
ˉ
x
n
δ
p
ˉ
y
n
δ
p
ˉ
z
n
δ
v
ˉ
y
b
n
δ
v
ˉ
z
b
n
δ
θ
ˉ
x
n
δ
θ
ˉ
y
n
δ
θ
ˉ
z
]
T
\begin{aligned} & \boldsymbol{G}_t=\left[\begin{array}{ccccc} \boldsymbol{I}_3 & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & {\left[\boldsymbol{R}_{b w}\right]_{y z}} & {\left[\left[\boldsymbol{v}^b\right]_{\times}\right]_{y z}} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \boldsymbol{I}_3 & \mathbf{0} & \mathbf{0} \end{array}\right] \\ & \boldsymbol{C}_t=\left[\begin{array}{ccc} \boldsymbol{I}_3 & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \boldsymbol{I}_2 & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \boldsymbol{I}_3 \end{array}\right] \\ & \boldsymbol{n}=\left[\begin{array}{llllllll} n_{\delta \bar{p}_x} & n_{\delta \bar{p}_y} & n_{\delta \bar{p}_z} & n_{\delta \bar{v}_y^b} & n_{\delta \bar{v}_z^b} & n_{\delta \bar{\theta}_x} & n_{\delta \bar{\theta}_y} & n_{\delta \bar{\theta}_z} \end{array}\right]^T \end{aligned}
Gt=
I3000[Rbw]yz00[[vb]×]yzI3000000
Ct=
I3000I2000I3
n=[nδpˉxnδpˉynδpˉznδvˉybnδvˉzbnδθˉxnδθˉynδθˉz]T随后的Kalman流程仍然与之前保持一致。
4. 融合点云特征的滤波方法
4.1 整体思路
以IMU做状态预测,以特征中的点-面距离、点-线距离为约束(观测),修正误差。
论文题目:LINS: A Lidar-Inertial State Estimator for Robust and Efficient Navigation
4.2 滤波模型
-
状态定义
位姿定义:
x w b k : = [ p w b k , q w b k ] \mathbf{x}_w^{b_k}:=\left[\mathbf{p}_w^{b_k}, \mathbf{q}_w^{b_k}\right] xwbk:=[pwbk,qwbk]相对位姿相关:
x b k + 1 b k : = [ p b k + 1 b k , v b k + 1 b k , q b k + 1 b k , b a , b g , g b k ] \mathbf{x}_{b_{k+1}}^{b_k}:=\left[\mathbf{p}_{b_{k+1}}^{b_k}, \mathbf{v}_{b_{k+1}}^{b_k}, \mathbf{q}_{b_{k+1}}^{b_k}, \mathbf{b}_a, \mathbf{b}_g, \mathbf{g}^{b_k}\right] xbk+1bk:=[pbk+1bk,vbk+1bk,qbk+1bk,ba,bg,gbk]状态量:
δ x : = [ δ p , δ v , δ θ , δ b a , δ b g , δ g ] \delta \mathbf{x}:=\left[\delta \mathbf{p}, \delta \mathbf{v}, \delta \boldsymbol{\theta}, \delta \mathbf{b}_a, \delta \mathbf{b}_g, \delta \mathbf{g}\right] δx:=[δp,δv,δθ,δba,δbg,δg]状态量修正:
x b k + 1 b k = − x b k + 1 b k ⊞ δ x = [ − p b k + 1 b k + δ p − v k b k + δ v − q b k + 1 b k ⊗ exp ( δ θ ) − b a + δ b a − b g + δ b g − g b k + δ g ] \mathbf{x}_{b_{k+1}}^{b_k}={ }^{-} \mathbf{x}_{b_{k+1}}^{b_k} \boxplus \boldsymbol{\delta} \mathbf{x}=\left[\begin{array}{c} -\mathbf{p}_{b_{k+1}}^{b_k}+\boldsymbol{\delta} \mathbf{p} \\ -\mathbf{v}_k^{b_k}+\boldsymbol{\delta} \mathbf{v} \\ -\mathbf{q}_{b_{k+1}}^{b_k} \otimes \exp (\boldsymbol{\delta} \boldsymbol{\theta}) \\ -\mathbf{b}_a+\boldsymbol{\delta} \mathbf{b}_a \\ -\mathbf{b}_g+\delta \mathbf{b}_g \\ -\mathbf{g}^{b_k}+\boldsymbol{\delta} \mathbf{g} \end{array}\right] xbk+1bk=−xbk+1bk⊞δx= −pbk+1bk+δp−vkbk+δv−qbk+1bk⊗exp(δθ)−ba+δba−bg+δbg−gbk+δg 其中 − x b k + 1 b k { }^{-} \mathbf{x}_{b_{k+1}}^{b_k} −xbk+1bk 表示预测值。 -
状态方程
δ x ˙ ( t ) = F t δ x ( t ) + G t w \delta \dot{\mathbf{x}}(t)=\mathbf{F}_t \delta \mathbf{x}(t)+\mathbf{G}_t \mathbf{w} δx˙(t)=Ftδx(t)+Gtw其中
F t = [ 0 I 0 0 0 0 0 0 − R t b k [ a ^ t ] × − R t b k 0 − I 3 0 0 − [ ω ^ t ] × 0 − I 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ] G t = [ 0 0 0 0 − R t b k 0 0 0 0 − I 3 0 0 0 0 I 3 0 0 0 0 I 3 0 0 0 0 ] w = [ n a T , n g T , n b a T , n b g T ] T \begin{aligned} \mathbf{F}_t & =\left[\begin{array}{cccccc} 0 & \mathbf{I} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & -\mathbf{R}_t^{b_k}\left[\hat{\mathbf{a}}_t\right]_{\times} & -\mathbf{R}_t^{b_k} & \mathbf{0} & -\mathbf{I}_3 \\ 0 & 0 & -\left[\hat{\omega}_t\right]_{\times} & \mathbf{0} & -\mathbf{I}_3 & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \end{array}\right] \quad \mathbf{G}_t=\left[\begin{array}{cccc} \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ -\mathbf{R}_t^{b_k} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & -\mathbf{I}_3 & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{I}_3 & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{I}_3 \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \end{array}\right] \\ \mathbf{w} & =\left[\mathbf{n}_a^T, \mathbf{n}_g^T, \mathbf{n}_{b_a}^T, \mathbf{n}_{b_g}^T\right]^T \end{aligned} Ftw= 000000I000000−Rtbk[a^t]×−[ω^t]×0000−Rtbk000000−I30000−I30000 Gt= 0−Rtbk000000−I3000000I3000000I30 =[naT,ngT,nbaT,nbgT]T -
观测方程
观测的计算与loam中前后帧匹配的思想一致,都是计算 点-面、点-线的残差
其中
p ^ i l k = ( R l b ) T ( R b k + 1 b k ( R l b p i l k + 1 + p l b ) + p b k + 1 b k − p l b ) \hat{\mathbf{p}}_i^{l_k}=\left(\mathbf{R}_l^b\right)^T\left(\mathbf{R}_{b_{k+1}}^{b_k}\left(\mathbf{R}_l^b \mathbf{p}_i^{l_{k+1}}+\mathbf{p}_l^b\right)+\mathbf{p}_{b_{k+1}}^{b_k}-\mathbf{p}_l^b\right) p^ilk=(Rlb)T(Rbk+1bk(Rlbpilk+1+plb)+pbk+1bk−plb)为了计算观测方程,需要计算残差对状态量的雅可比
H k = ∂ f ∂ p ^ i l k ⋅ ∂ p ^ i l k ∂ δ x \mathbf{H}_k=\frac{\partial \mathbf{f}}{\partial \hat{\mathbf{p}}_i^{l_k}} \cdot \frac{\partial \hat{\mathbf{p}}_i^{l_k}}{\partial \delta \mathbf{x}} Hk=∂p^ilk∂f⋅∂δx∂p^ilk上式中包含两部分,第一部分已经讲过,此处只推导第 二部分。
除了旋转和平移外,残差项对其它量的导数均为 0 。
a.对平移求导
∂ p ^ i l k ∂ δ p = ∂ [ ( R l b ) ⊤ ( R b k + 1 b k ( R l b p i l k + 1 + p l b ) + p b k + 1 b k − p l b ) ] ∂ δ p = ∂ [ ( R l b ) ⊤ p b k + 1 b k ] ∂ δ p = ( R l b ) ⊤ \begin{aligned} & \frac{\partial \hat{\mathbf{p}}_i^{l_k}}{\partial \delta \mathbf{p}} \\ = & \frac{\partial\left[\left(\mathbf{R}_l^b\right)^{\top}\left(\mathbf{R}_{b_{k+1}}^{b_k}\left(\mathbf{R}_l^b \mathbf{p}_i^{l_{k+1}}+\mathbf{p}_l^b\right)+\mathbf{p}_{b_{k+1}}^{b_k}-\mathbf{p}_l^b\right)\right]}{\partial \delta \mathbf{p}} \\ = & \frac{\partial\left[\left(\mathbf{R}_l^b\right)^{\top} \mathbf{p}_{b_{k+1}}^{b_k}\right]}{\partial \delta \mathbf{p}} \\ = & \left(\mathbf{R}_l^b\right)^{\top} \end{aligned} ===∂δp∂p^ilk∂δp∂[(Rlb)⊤(Rbk+1bk(Rlbpilk+1+plb)+pbk+1bk−plb)]∂δp∂[(Rlb)⊤pbk+1bk](Rlb)⊤b.对旋转求导
∂ p ^ i l k ∂ δ x = ∂ [ ( R l b ) T ( R b k + 1 b k ( R l b p i l k + 1 + p l b ) + p b k + 1 b k − p l b ) ] ∂ δ θ = ( R l b ) T ∂ [ R b k + 1 b k ( R l b p i l k + 1 + p l b ) ] ∂ R b k + 1 b k ∂ R b k + 1 b k ∂ δ θ = − ( R l b ) T [ R b k + 1 b k ( R l b p i l k + 1 + p l b ) ] × R b k + 1 b k J r − 1 ( θ ) = − ( R l b ) T R b k + 1 b k [ R l b p i l k + 1 + p l b ] × R b k b k + 1 R b k + 1 b k J r − 1 ( θ ) = − ( R l b ) T R b k + 1 b k [ R l b p i l k + 1 + p l b ] × J r − 1 ( θ ) \begin{aligned} & \frac{\partial \hat{\mathbf{p}}_i^{l_k}}{\partial \delta \mathbf{x}} \\ = & \frac{\partial\left[\left(\mathbf{R}_l^b\right)^T\left(\mathbf{R}_{b_{k+1}}^{b_k}\left(\mathbf{R}_l^b \mathbf{p}_i^{l_{k+1}}+\mathbf{p}_l^b\right)+\mathbf{p}_{b_{k+1}}^{b_k}-\mathbf{p}_l^b\right)\right]}{\partial \delta \boldsymbol{\theta}} \\ = & \left(\mathbf{R}_l^b\right)^T \frac{\partial\left[\mathbf{R}_{b_{k+1}}^{b_k}\left(\mathbf{R}_l^b \mathbf{p}_i^{l_{k+1}}+\mathbf{p}_l^b\right)\right]}{\partial \mathbf{R}_{b_{k+1}}^{b_k}} \frac{\partial \mathbf{R}_{b_{k+1}}^{b_k}}{\partial \delta \boldsymbol{\theta}} \\ = & -\left(\mathbf{R}_l^b\right)^T\left[\mathbf{R}_{b_{k+1}}^{b_k}\left(\mathbf{R}_l^b \mathbf{p}_i^{l_{k+1}}+\mathbf{p}_l^b\right)\right]_{\times} \mathbf{R}_{b_{k+1}}^{b_k} \boldsymbol{J}_r^{-1}(\boldsymbol{\theta}) \\ = & -\left(\mathbf{R}_l^b\right)^T \mathbf{R}_{b_{k+1}}^{b_k}\left[\mathbf{R}_l^b \mathbf{p}_i^{l_{k+1}}+\mathbf{p}_l^b\right]_{\times} \mathbf{R}_{b_k}^{b_{k+1}} \mathbf{R}_{b_{k+1}}^{b_k} \boldsymbol{J}_r^{-1}(\boldsymbol{\theta}) \\ = & -\left(\mathbf{R}_l^b\right)^T \mathbf{R}_{b_{k+1}}^{b_k}\left[\mathbf{R}_l^b \mathbf{p}_i^{l_{k+1}}+\mathbf{p}_l^b\right]_{\times} \boldsymbol{J}_r^{-1}(\boldsymbol{\theta}) \end{aligned} =====∂δx∂p^ilk∂δθ∂[(Rlb)T(Rbk+1bk(Rlbpilk+1+plb)+pbk+1bk−plb)](Rlb)T∂Rbk+1bk∂[Rbk+1bk(Rlbpilk+1+plb)]∂δθ∂Rbk+1bk−(Rlb)T[Rbk+1bk(Rlbpilk+1+plb)]×Rbk+1bkJr−1(θ)−(Rlb)TRbk+1bk[Rlbpilk+1+plb]×Rbkbk+1Rbk+1bkJr−1(θ)−(Rlb)TRbk+1bk[Rlbpilk+1+plb]×Jr−1(θ)预测 :
δ x t τ = ( I + F t τ Δ t ) δ x t τ − 1 P t τ = ( I + F t τ Δ t ) P t τ − 1 ( I + F t τ Δ t ) T + ( G i τ Δ t ) Q ( G t t Δ t ) T K k , j = P k H k , j T ( H k , j P k H k , j T + J k , j M k J k , j T ) − 1 \begin{aligned} & \delta \mathbf{x}_{t_\tau}=\left(\mathbf{I}+\mathbf{F}_{t_\tau} \Delta t\right) \delta \mathbf{x}_{t_{\tau-1}} \\ & \mathbf{P}_{t_\tau}=\left(\mathbf{I}+\mathbf{F}_{t_\tau} \Delta t\right) \mathbf{P}_{t_{\tau-1}}\left(\mathbf{I}+\mathbf{F}_{t_\tau} \Delta t\right)^T+\left(\mathbf{G}_{i_\tau} \Delta t\right) \mathbf{Q}\left(\mathbf{G}_{t_t} \Delta t\right)^T \\ & \mathbf{K}_{k, j}=\mathbf{P}_k \mathbf{H}_{k, j}^T\left(\mathbf{H}_{k, j} \mathbf{P}_k \mathbf{H}_{k, j}^T+\mathbf{J}_{k, j} \mathbf{M}_k \mathbf{J}_{k, j}^T\right)^{-1} \end{aligned} δxtτ=(I+FtτΔt)δxtτ−1Ptτ=(I+FtτΔt)Ptτ−1(I+FtτΔt)T+(GiτΔt)Q(GttΔt)TKk,j=PkHk,jT(Hk,jPkHk,jT+Jk,jMkJk,jT)−1迭代观测: Δ x j = K k , j ( H k , j δ x j − f ( − x b k + 1 b k ⊞ δ x j ) ) \Delta \mathbf{x}_j=\mathbf{K}_{k, j}\left(\mathbf{H}_{k, j} \delta \mathbf{x}_j-f\left({ }^{-} \mathbf{x}_{b_{k+1}}^{b_k} \boxplus \boldsymbol{\delta} \mathbf{x}_j\right)\right) Δxj=Kk,j(Hk,jδxj−f(−xbk+1bk⊞δxj))
δ x j + 1 = δ x j + Δ x j \delta \mathbf{x}_{j+1}=\delta \mathbf{x}_j+\Delta \mathbf{x}_j δxj+1=δxj+Δxj后验方差: P k + 1 = ( I − K k , n H k , n ) P k ( I − K k , n H k , n ) T + K k , n M k K k , n T \mathbf{P}_{k+1}=\left(\mathbf{I}-\mathbf{K}_{k, n} \mathbf{H}_{k, n}\right) \mathbf{P}_k\left(\mathbf{I}-\mathbf{K}_{k, n} \mathbf{H}_{k, n}\right)^T+\mathbf{K}_{k, n} \mathbf{M}_k \mathbf{K}_{k, n}^T Pk+1=(I−Kk,nHk,n)Pk(I−Kk,nHk,n)T+Kk,nMkKk,nT
4.3 位姿更新
x w b k + 1 = [ p w b k + 1 q w b k + 1 ] = [ R b k b k + 1 ( p w b k − p b k + 1 b k ) q b k b k + 1 ⊗ q w b k ] \mathbf{x}_w^{b_{k+1}}=\left[\begin{array}{c} \mathbf{p}_w^{b_{k+1}} \\ \mathbf{q}_w^{b_{k+1}} \end{array}\right]=\left[\begin{array}{c} \mathbf{R}_{b_k}^{b_{k+1}}\left(\mathbf{p}_w^{b_k}-\mathbf{p}_{b_{k+1}}^{b_k}\right) \\ \mathbf{q}_{b_k}^{b_{k+1}} \otimes \mathbf{q}_w^{b_k} \end{array}\right] xwbk+1=[pwbk+1qwbk+1]=[Rbkbk+1(pwbk−pbk+1bk)qbkbk+1⊗qwbk]
4.4 相似工作
论文题目:FAST-LIO: A Fast, Robust LiDAR-inertial Odometry Package by Tightly-Coupled Iterated Kalman Filter