乘法的误差传递公式 - CSDN
文章推薦指數: 80 %
加法中的误差传递: X=u±v 则X的均方差为: σX=sqrt(σu^2+σv^2)...乘法中的误差传递: 除法中的误差传递: 有限次幂的误差的传播: 可以使用蒙特卡罗法来验证其误差: 如下 ...
精华内容
下载资源
问答
我要提问
误差传递公式
万次阅读
2020-12-3010:50:33
1.加法中的误差传递:即:若有X=u±v则X的均方差为:σX^2=sqrt(σu^2+σv^2);2.乘法中的误差传递:3.除法中的误差传递:4.有限次幂的误差的传播:
1.加法中的误差传递:
即:若有X=u±v
则X的均方差为:
σX^2=sqrt(σu^2+σv^2);
2.乘法中的误差传递:
3.除法中的误差传递:
4.有限次幂的误差的传播:
转载自其他博客,侵删。
收起
展开全文
数学
协方差、方差、误差传递公式
万次阅读
2018-10-0120:06:51
加法中的误差传递:X=u±v则X的均方差为:σX=sqrt(σu^2+σv^2)...乘法中的误差传递: 除法中的误差传递: 有限次幂的误差的传播: 可以使用蒙特卡罗法来验证其误差:如下面的程序...
转载自:https://www.cnblogs.com/heaventian/archive/2012/11/24/2786241.html
加法中的误差传递:
X=u±v
则X的均方差为:
σX=sqrt(σu^2+σv^2);
乘法中的误差传递:
除法中的误差传递:
有限次幂的误差的传播:
可以使用蒙特卡罗法来验证其误差:
如下面的程序用来验证出发的误差:
N=1e6;x=10+randn(N,1);y=5+randn(N,1)*2;std(x./y)mean(x./y)
收起
展开全文
协方差
VINS-Mono之IMU预积分,预积分误差、协方差及误差对状态量雅克比矩阵的递推方程的推导
千次阅读
2020-02-0711:44:14
关于公式(1)(1)(1)的推导,这里首先引入四元素左乘右乘及导数定理:根据《视觉SLAM14讲》3.4.2的四元数乘法,我们引入左乘和右乘符号如下:qa⊗qb=R(qb)qa=[sbzb−ybxb−zbsbxbybyb−xbsbzb−xb−yb−zbsb]...
文章目录
1.前言2.IMU模型3.基于世界坐标系下的IMU运动模型3.1连续形式下的IMU运动模型3.2离散形式下的IMU运动模型3.2.1欧拉法离散形式3.2.2中值法离散形式
4.IMU预积分(基于第K帧IMUbody坐标系下的运动模型)4.1连续形式下的IMU运动模型4.2离散形式下的IMU运动模型4.2.1两帧之间PVQ增量的欧拉法离散形式4.2.2两帧之间PVQ增量的中值法离散形式
5.非线性方程的误差递推方程5.1一段时间内多个IMU数据积分形成的预积分量的协方差计算5.2非线性方程的误差递推方程推导5.2.1基于一阶泰勒展开的误差递推方程5.2.2基于误差随时间变化(误差导数)的递推方程5.2.3基于泰勒展开和误差随时间变化的方法对比
6.PVQ增量误差、协方差及雅克比矩阵的递推方程6.1连续形式下的误差、协方差及雅克比矩阵的递推方程6.2离散形式下(中值法)增量误差、协方差及雅克比矩阵的递推方程
1.前言
本博客借鉴了崔华坤的《VINS论文推导及代码解析》和VINS-Mono理论学习——IMU预积分Pre-integration(Jacobian协方差)的内容,因为确实写得太好了,然后有些地方加入自己一些理解。
VINS-MONO论文中的IV-B.IMUPre-integration介绍了IMU预积分模型,Foster的倆篇论文对IMU预积分理论进行详细分析。
值得注意的是,本文使用了四元素中Hamilton(可参考《QuaternionkinematicsfortheerrorstateKalmanfilter》,右手系)的表示和JPL(可参考《IndirectKalmanFilterfor3DAttitudeEstimation》,左手系)表示,因为Vins-mono里面也是用了俩种表示混用,俩者的具体区别可以看《QuaternionkinematicsfortheerrorstateKalmanfilter》的第17页,且论文中用到的JPL的右乘左乘定义还与《IndirectKalmanFilterfor3DAttitudeEstimation》不同,这种问题估计让无数人心塞。
/(ㄒoㄒ)/~~
为什么需要对IMU进行预积分?传统捷联惯性导航算法,在已知
k
k
k时刻下的IMU状态量(姿态、速度和位移)情况下,利用IMU测量的线加速度和角速度,通过积分预算得到
k
+
1
k+1
k+1时刻下的状态量。
然而在非线性优化的VIO中,各个节点的状态量都是估计值,当算法对这些状态量优化时,每次调整都需要在它们之间重新积分,导致绝对位姿被优化时对状态量进行重复积分。
IMU预积分的提出使得优化算法可对IMU的相对测量进行处理,使它与绝对位姿解耦,或者只要线性运算就可以进行矫正。
2.IMU模型
IMU测量值包括加速度计得到的(测量值)线加速度
a
t
^
\hat{a_{t}}
at^和陀螺仪得到的角加速度
w
^
t
\hat{w}_{t}
w^t【论文式(1)】
a
^
t
=
a
t
+
b
a
t
+
R
w
t
g
w
+
n
a
\hat{a}_{t}=a_{t}+b_{at}+R_{w}^{t}g^{w}+n_{a}
a^t=at+bat+Rwtgw+na
w
^
t
=
w
t
+
b
w
t
+
n
w
\hat{w}_{t}=w_{t}+b_{wt}+n_{w}
w^t=wt+bwt+nw其中
t
t
t下标表示在IMU的体(body)坐标系下,
a
t
a_{t}
at、
w
t
w_{t}
wt分别表示IMU真实的线加速度和角速度,并受到加速度偏置(bias)
b
a
t
b_{at}
bat、陀螺仪偏置
b
w
t
b_{wt}
bwt和附加噪声
n
a
n_a
na、
n
n
w
n_{n_{w}}
nnw的影响。
计算得到的线加速度
a
t
^
\hat{a_{t}}
at^是重力加速度和物体加速度的合矢量。
假设附加噪声为高斯噪声:
n
a
∼
(
0
,
σ
a
2
)
,
n
a
∼
(
0
,
σ
w
2
)
n_{a}\sim(0,\sigma_{a}^{2}),\n_{a}\sim(0,\sigma_{w}^{2})
na∼(0,σa2), na∼(0,σw2)加速度计偏置和陀螺仪偏置被建模为随机游走,其导数为高斯分布:【论文式(2)】
b
˙
a
t
=
n
b
a
,
b
˙
w
t
=
n
b
w
\dot{b}_{at}=n_{ba},\\dot{b}_{wt}=n_{bw}
b˙at=nba, b˙wt=nbw
n
b
a
∼
N
(
0
,
σ
b
a
2
)
,
n
b
w
∼
N
(
0
,
σ
b
w
2
)
n_{ba}\simN(0,\sigma_{ba}^{2}),\n_{bw}\simN(0,\sigma_{bw}^{2})
nba∼N(0,σba2), nbw∼N(0,σbw2)
3.基于世界坐标系下的IMU运动模型
3.1连续形式下的IMU运动模型
对于图像帧
k
k
k和
k
+
1
k+1
k+1,IMUbody坐标系对应为
b
k
b_{k}
bk和
b
k
+
1
b_{k+1}
bk+1,位置、速度和姿态状态值PVQ(Pose、Velocity、Quaternion)可以根据
[
t
k
,
t
k
+
1
]
[t_{k},t_{k+1}]
[tk,tk+1]时间间隔内的IMU测量值,在世界坐标系下进行传递:【论文式(3)(4)】
p
b
k
+
1
w
=
p
b
k
w
+
v
b
k
w
Δ
t
k
+
∫
∫
t
∈
[
t
k
,
t
k
+
1
]
(
R
t
w
(
a
^
t
−
b
a
t
−
n
a
)
−
g
w
)
d
t
2
p^{w}_{b_{k+1}}=p^{w}_{bk}+v_{b_{k}}^{w}\Deltat_{k}+\int\int_{t\in[t_{k},t_{k+1}]}(R_{t}^{w}(\hat{a}_{t}-b_{at}-n_{a})-g^{w})dt^{2}
pbk+1w=pbkw+vbkwΔtk+∫∫t∈[tk,tk+1](Rtw(a^t−bat−na)−gw)dt2
v
b
k
+
1
w
=
v
b
k
w
+
∫
t
∈
[
t
k
,
t
k
+
1
]
(
R
t
w
(
a
^
t
−
b
a
t
−
n
a
)
−
g
w
)
d
t
v_{b_{k+1}}^{w}=v_{bk}^{w}+\int_{t\in[t_{k},t_{k+1}]}(R^{w}_{t}(\hat{a}_{t}-b_{at}-n_{a})-g^{w})dt
vbk+1w=vbkw+∫t∈[tk,tk+1](Rtw(a^t−bat−na)−gw)dt
q
b
k
+
1
w
=
q
b
k
w
⊗
∫
t
∈
[
t
k
,
t
k
+
1
]
1
2
q
t
b
k
⊗
[
(
w
^
t
−
b
w
t
−
n
w
)
0
]
d
t
=
q
b
k
⊗
∫
t
∈
[
t
k
,
t
k
+
1
]
1
2
Ω
(
w
^
t
−
b
w
t
−
n
w
)
q
t
b
k
d
t
(1)
q_{b_{k+1}}^{w}=q_{b_{k}}^{w}\otimes\int_{t\in[t_{k},t_{k+1}]}\frac{1}{2}q_{t}^{b_{k}}\otimes\begin{bmatrix}(\hat{w}_{t}-b_{wt}-n_{w})\\0\end{bmatrix}dt\\=q_{bk}\otimes\int_{t\in[t_{k},t_{k+1}]}\frac{1}{2}\Omega(\hat{w}_{t}-b_{wt}-n_{w})q_{t}^{bk}dt\tag{1}
qbk+1w=qbkw⊗∫t∈[tk,tk+1]21qtbk⊗[(w^t−bwt−nw)0]dt=qbk⊗∫t∈[tk,tk+1]21Ω(w^t−bwt−nw)qtbkdt(1)
Ω
(
w
)
=
[
−
[
w
]
×
w
−
w
T
0
]
,
[
w
]
×
=
[
0
−
w
z
w
y
w
z
0
−
w
x
−
w
y
w
x
0
]
\Omega(w)=\begin{bmatrix}-[w]_{\times}&w\\-w^{T}&0\end{bmatrix},\[w]_{\times}=\begin{bmatrix}0&-w_{z}&w_{y}\\w_{z}&0&-w_{x}\\-w_{y}&w_{x}&0\end{bmatrix}
Ω(w)=[−[w]×−wTw0], [w]×=⎣⎡0wz−wy−wz0wxwy−wx0⎦⎤其中
Δ
t
k
\Deltat_{k}
Δtk是
[
t
k
,
t
k
+
1
]
[t_{k},t_{k+1}]
[tk,tk+1]之间的时间间隔,
R
t
w
R_{t}^{w}
Rtw为t时刻IMUbody坐标系到世界坐标系的旋转矩阵,
q
t
b
k
q_{t}^{bk}
qtbk为用四元素表示的
t
t
t时刻从IMUbody坐标系到
k
k
k时刻IMUbody坐标系的旋转,这里的四元素实部在后,虚部在前,为了与论文保持一致,即四元素JPL表达。
这里的
Ω
(
w
)
\Omega(w)
Ω(w)表示四元素右乘。
关于公式
(
1
)
(1)
(1)的推导,这里首先引入四元素左乘右乘及导数定理:这里与《IndirectKalmanFilterfor3DAttitudeEstimation》中公式(10)不同,主要是因为左右手系不同,我们引入左乘和右乘符号如下:
q
a
⊗
q
b
=
R
(
q
b
)
q
a
=
[
s
b
z
b
−
y
b
x
b
−
z
b
s
b
x
b
y
b
y
b
−
x
b
s
b
z
b
−
x
b
−
y
b
−
z
b
s
b
]
[
x
a
y
a
z
a
s
a
]
=
L
(
q
a
)
q
b
=
[
s
a
−
z
a
y
a
x
a
z
a
s
a
−
x
a
y
a
−
y
a
x
a
s
a
z
a
−
x
a
−
y
a
−
z
a
s
a
]
[
x
b
y
b
z
b
s
b
]
q_{a}\otimesq_{b}=R(q_{b})q_{a}=\begin{bmatrix}s_{b}&z_{b}&-y_{b}&x_{b}\\-z_{b}&s_{b}&x_{b}&y_{b}\\y_{b}&-x_{b}&s_{b}&z_{b}\\-x_{b}&-y_{b}&-z_{b}&s_{b}\end{bmatrix}\begin{bmatrix}x_{a}\\y_{a}\\z_{a}\\s_{a}\end{bmatrix}\\=L(q_{a})q_{b}=\begin{bmatrix}s_{a}&-z_{a}&y_{a}&x_{a}\\z_{a}&s_{a}&-x_{a}&y_{a}\\-y_{a}&x_{a}&s_{a}&z_{a}\\-x_{a}&-y_{a}&-z_{a}&s_{a}\end{bmatrix}\begin{bmatrix}x_{b}\\y_{b}\\z_{b}\\s_{b}\end{bmatrix}
qa⊗qb=R(qb)qa=⎣⎢⎢⎡sb−zbyb−xbzbsb−xb−yb−ybxbsb−zbxbybzbsb⎦⎥⎥⎤⎣⎢⎢⎡xayazasa⎦⎥⎥⎤=L(qa)qb=⎣⎢⎢⎡saza−ya−xa−zasaxa−yaya−xasa−zaxayazasa⎦⎥⎥⎤⎣⎢⎢⎡xbybzbsb⎦⎥⎥⎤为了简化,令
q
=
[
x
y
z
s
]
=
[
w
s
]
q=[x\y\z\s]=[w\s]
q=[x y z s]=[w s],则有:
R
(
q
)
=
Ω
(
w
)
+
s
I
4
×
4
=
[
−
[
w
]
×
w
−
w
T
0
]
+
s
I
4
×
4
R(q)=\Omega(w)+sI_{4\times4}=\begin{bmatrix}-[w]_{\times}&w\\-w^{T}&0\end{bmatrix}+sI_{4\times4}
R(q)=Ω(w)+sI4×4=[−[w]×−wTw0]+sI4×4
L
(
q
)
=
Ψ
(
w
)
+
s
I
4
×
4
=
[
[
w
]
×
w
−
w
T
0
]
+
s
I
4
×
4
L(q)=\Psi(w)+sI_{4\times4}=\begin{bmatrix}[w]_{\times}&w\\-w^{T}&0\end{bmatrix}+sI_{4\times4}
L(q)=Ψ(w)+sI4×4=[[w]×−wTw0]+sI4×4对于四元素的求导,我们定义
q
t
q_{t}
qt为
t
t
t时刻下的单位四元素(这里的四元素实部在后,虚部在前),
w
w
w为
q
t
q_{t}
qt确定的角速度,则关于
q
t
q_{t}
qt的导数为:
q
˙
t
=
1
2
[
−
[
w
]
×
w
−
w
T
0
]
q
t
=
1
2
Ω
(
w
)
q
t
=
1
2
R
(
[
w
0
]
)
q
t
=
1
2
q
t
⊗
[
w
0
]
\dot{q}_{t}=\frac{1}{2}\begin{bmatrix}-[w]_{\times}&w\\-w^{T}&0\end{bmatrix}q_{t}=\frac{1}{2}\Omega(w)q_{t}=\frac{1}{2}R(\begin{bmatrix}w\\0\end{bmatrix})q_{t}=\frac{1}{2}q_{t}\otimes\begin{bmatrix}w\\0\end{bmatrix}
q˙t=21[−[w]×−wTw0]qt=21Ω(w)qt=21R([w0])qt=21qt⊗[w0]
当看到
q
˙
=
1
2
q
t
⊗
[
0
w
]
\dot{q}=\frac{1}{2}q_{t}\otimes\begin{bmatrix}0\\w\end{bmatrix}
q˙=21qt⊗[0w]说明该四元素实部在前,虚部在后(Hamilton表示),只是表示不同而已。
因此对于IMU连续形式下的旋转状态(用四元素表示)推导,我们有:
q
b
k
+
1
w
=
q
b
k
w
⊗
q
b
k
+
1
b
k
=
q
b
k
w
⊗
∫
t
∈
[
t
k
,
t
k
+
1
]
q
t
˙
b
k
d
t
=
q
b
k
w
⊗
∫
t
∈
[
t
k
,
t
k
+
1
]
1
2
q
t
b
k
⊗
[
w
t
b
k
0
]
d
t
=
q
b
k
w
⊗
∫
t
∈
[
t
k
,
t
k
+
1
]
1
2
q
t
b
k
⊗
[
w
t
^
−
b
w
t
−
n
w
0
]
d
t
=
=
q
b
k
⊗
∫
t
∈
[
t
k
,
t
k
+
1
]
1
2
Ω
(
w
^
t
−
b
w
t
−
n
w
)
q
t
b
k
d
t
q_{b_{k+1}^{w}}=q_{b_{k}}^{w}\otimesq_{b_{k+1}}^{b_{k}}=q_{b_{k}}^{w}\otimes\int_{t\in[t_{k},t_{k+1}]}\dot{q_{t}}^{b_{k}}dt=q_{bk}^{w}\otimes\int_{t\in[t_{k},t_{k+1}]}\frac{1}{2}q_{t}^{b_{k}}\otimes\begin{bmatrix}w_{t}^{b_{k}}\\0\end{bmatrix}dt\\=q_{b_{k}}^{w}\otimes\int_{t\in[t_{k},t_{k+1}]}\frac{1}{2}q_{t}^{b_{k}}\otimes\begin{bmatrix}\hat{w_{t}}-b_{wt}-n_{w}\\0\end{bmatrix}dt\\==q_{bk}\otimes\int_{t\in[t_{k},t_{k+1}]}\frac{1}{2}\Omega(\hat{w}_{t}-b_{wt}-n_{w})q_{t}^{bk}dt
qbk+1w=qbkw⊗qbk+1bk=qbkw⊗∫t∈[tk,tk+1]qt˙bkdt=qbkw⊗∫t∈[tk,tk+1]21qtbk⊗[wtbk0]dt=qbkw⊗∫t∈[tk,tk+1]21qtbk⊗[wt^−bwt−nw0]dt==qbk⊗∫t∈[tk,tk+1]21Ω(w^t−bwt−nw)qtbkdt
3.2离散形式下的IMU运动模型
3.2.1欧拉法离散形式
使用欧拉法,即
k
+
1
k+1
k+1时刻的位姿是用第
k
k
k时刻的测量值
a
^
b
k
\hat{a}_{b_{k}}
a^bk,
w
^
b
k
\hat{w}_{b_{k}}
w^bk来计算的:
p
b
k
+
1
w
=
p
b
k
w
+
v
b
k
w
Δ
t
k
+
1
2
a
^
b
k
δ
t
2
p^{w}_{b_{k+1}}=p_{b_{k}}^{w}+v_{b_{k}}^{w}\Deltat_{k}+\frac{1}{2}\hat{a}_{b_{k}}\deltat^{2}
pbk+1w=pbkw+vbkwΔtk+21a^bkδt2
v
b
k
+
1
w
=
v
b
k
w
+
a
^
b
k
δ
t
v_{b_{k+1}}^{w}=v_{b_{k}}^{w}+\hat{a}_{b_{k}}\deltat
vbk+1w=vbkw+a^bkδt
q
b
k
+
1
w
=
q
b
k
w
⊗
[
1
1
2
w
^
b
k
δ
t
]
q^{w}_{b_{k+1}}=q^{w}_{b_{k}}\otimes\begin{bmatrix}1\\\frac{1}{2}\hat{w}_{b_{k}}\deltat\end{bmatrix}
qbk+1w=qbkw⊗[121w^bkδt]其中
a
^
b
k
=
q
b
k
w
(
a
b
k
−
b
a
k
)
−
g
w
\hat{a}_{b_{k}}=q_{b_{k}}^{w}(a_{b_{k}}-b_{ak})-g^{w}
a^bk=qbkw(abk−bak)−gw
w
^
b
k
=
w
b
k
−
b
w
k
\hat{w}_{b_{k}}=w_{b_{k}}-b_{wk}
w^bk=wbk−bwk
3.2.2中值法离散形式
使用中值法,即
k
+
1
k+1
k+1时刻的位姿是用俩个时刻
k
k
k和
k
+
1
k+1
k+1测量值
a
a
a,
w
w
w的平均值来计算的:
p
b
k
+
1
w
=
p
b
k
w
+
v
b
k
w
Δ
t
k
+
1
2
a
^
ˉ
t
δ
t
2
p^{w}_{b_{k+1}}=p_{b_{k}}^{w}+v_{b_{k}}^{w}\Deltat_{k}+\frac{1}{2}\bar{\hat{a}}_{t}\deltat^{2}
pbk+1w=pbkw+vbkwΔtk+21a^ˉtδt2
v
b
k
+
1
w
=
v
b
k
w
+
a
^
ˉ
t
δ
t
v_{b_{k+1}}^{w}=v_{b_{k}}^{w}+\bar{\hat{a}}_{t}\deltat
vbk+1w=vbkw+a^ˉtδt
q
b
k
+
1
w
=
q
b
k
w
⊗
[
1
1
2
w
^
ˉ
t
δ
t
]
q^{w}_{b_{k+1}}=q^{w}_{b_{k}}\otimes\begin{bmatrix}1\\\frac{1}{2}\bar{\hat{w}}_{t}\deltat\end{bmatrix}
qbk+1w=qbkw⊗[121w^ˉtδt]其中
a
^
ˉ
t
=
1
2
[
q
b
k
w
(
a
b
k
−
b
a
k
)
−
g
w
+
q
b
k
+
1
w
(
a
b
k
+
1
−
b
a
k
+
1
)
−
g
w
]
\bar{\hat{a}}_{t}=\frac{1}{2}[q_{b_{k}}^{w}(a_{b_{k}}-b_{ak})-g^{w}+q_{b_{k+1}}^{w}(a_{b_{k+1}}-b_{ak+1})-g^{w}]
a^ˉt=21[qbkw(abk−bak)−gw+qbk+1w(abk+1−bak+1)−gw]
w
^
ˉ
t
=
1
2
(
w
b
k
−
b
w
k
+
w
b
k
+
1
−
b
w
k
+
1
)
=
1
2
(
w
b
k
+
w
b
k
+
1
)
−
b
w
k
\bar{\hat{w}}_{t}=\frac{1}{2}(w_{b_{k}}-b_{wk}+w_{b_{k+1}}-b_{wk+1})\\=\frac{1}{2}(w_{b_{k}}+w_{b_{k+1}})-b_{w{k}}
w^ˉt=21(wbk−bwk+wbk+1−bwk+1)=21(wbk+wbk+1)−bwk假设在短时间内加速度计和陀螺仪的偏置不变,则有:
b
a
k
=
b
a
k
+
1
,
b
w
k
=
b
w
k
+
1
b_{ak}=b_{ak+1},\b_{wk}=b_{wk+1}
bak=bak+1, bwk=bwk+1
4.IMU预积分(基于第K帧IMUbody坐标系下的运动模型)
通过公式
(
1
)
(1)
(1)可以看到,IMU的积分需要依赖与第
k
k
k帧的
v
v
v和
R
R
R(基于世界坐标系下的),当我们在后端进行非线性优化时,需要迭代更新第
k
k
k帧的
v
v
v和
R
R
R,这将导致我们需要根据每次迭代后的值重新进行积分,这将非常耗时。
我们考虑将优化变量从第
k
k
k帧到第
k
+
1
k+1
k+1帧的IMU积分项中分离开来。
4.1连续形式下的IMU运动模型
IMU预积分的思想就是将参考坐标系从世界坐标系
w
w
w调整为第
k
k
k帧的IMUbody坐标系
b
k
b_{k}
bk下,可通过在等式俩端同时乘以
R
w
b
k
R^{b_{k}}_{w}
Rwbk得到:【论文式[5][6]】
R
w
b
k
p
b
k
+
1
w
=
R
w
b
k
(
p
b
k
w
+
v
b
k
w
Δ
t
k
−
1
2
g
w
Δ
t
k
2
)
+
α
b
k
+
1
b
k
R^{b_{k}}_{w}p^{w}_{b_{k+1}}=R^{b_{k}}_{w}(p_{b_{k}}^{w}+v^{w}_{b_{k}}\Deltat_{k}-\frac{1}{2}g^{w}\Deltat_{k}^{2})+\alpha^{b_{k}}_{b_{k+1}}
Rwbkpbk+1w=Rwbk(pbkw+vbkwΔtk−21gwΔtk2)+αbk+1bk
R
w
b
k
v
b
k
+
1
w
=
R
w
b
k
(
v
b
k
w
−
g
w
Δ
t
k
)
+
β
b
k
+
1
b
k
R^{b_{k}}_{w}v_{b_{k+1}}^{w}=R_{w}^{b_{k}}(v_{b_{k}}^{w}-g^{w}\Deltat_{k})+\beta^{b_{k}}_{b_{k+1}}
Rwbkvbk+1w=Rwbk(vbkw−gwΔtk)+βbk+1bk
q
w
b
k
⊗
q
b
k
+
1
w
=
γ
b
k
+
1
b
k
(2)
q_{w}^{b_{k}}\otimesq_{b_{k+1}}^{w}=\gamma_{b_{k+1}}^{b_{k}}\tag{2}
qwbk⊗qbk+1w=γbk+1bk(2)其中
α
b
k
+
1
b
k
=
∫
∫
t
∈
[
t
k
,
t
k
+
1
]
R
t
b
k
(
a
^
t
−
b
a
t
−
n
a
)
d
t
2
\alpha^{b_{k}}_{b_{k+1}}=\int\int_{t\in[t_{k},t_{k+1}]}R_{t}^{b_{k}}(\hat{a}_{t}-b_{at}-n_{a})dt^{2}
αbk+1bk=∫∫t∈[tk,tk+1]Rtbk(a^t−bat−na)dt2
β
b
k
+
1
b
k
=
∫
t
∈
[
t
k
,
t
k
+
1
]
R
t
b
k
(
a
^
t
−
b
a
t
−
n
a
)
d
t
\beta_{b_{k+1}}^{b_{k}}=\int_{t\in[t_{k},t_{k+1}]}R_{t}^{b_{k}}(\hat{a}_{t}-b_{at}-n_{a})dt
βbk+1bk=∫t∈[tk,tk+1]Rtbk(a^t−bat−na)dt
γ
b
k
+
1
=
∫
t
∈
[
t
k
,
t
k
+
1
]
1
2
Ω
(
w
^
t
−
b
w
t
−
n
w
)
γ
t
b
k
d
t
\gamma_{b_{k+1}}=\int_{t\in[t_{k},t_{k+1}]}\frac{1}{2}\Omega(\hat{w}_{t}-b_{wt}-n_{w})\gamma_{t}^{b_{k}}dt
γbk+1=∫t∈[tk,tk+1]21Ω(w^t−bwt−nw)γtbkdt
此时的积分结果
α
b
k
+
1
b
k
\alpha^{b_{k}}_{b_{k+1}}
αbk+1bk、
β
b
k
+
1
b
k
\beta^{b_{k}}_{b_{k+1}}
βbk+1bk、
γ
b
k
+
1
b
k
\gamma^{b_{k}}_{b_{k+1}}
γbk+1bk可以理解为
b
k
+
1
b_{k+1}
bk+1对
b
k
b_{k}
bk的相对运动量,
b
k
b_{k}
bk的状态并不会对其产生影响,因此将其作为非线性优化变量,可以避免状态的重复传递。
注意,这是在假设IMU偏置
b
a
b_{a}
ba、
b
w
b_{w}
bw已经确定的情况下,实际上偏置也是需要优化的变量,那么每次迭代时,
b
a
b_{a}
ba、
b
w
b_{w}
bw发生改变,得重新根据公式求得所有帧之间的IMU预积分。
当偏置变换很小时,可以将
α
b
k
+
1
b
k
\alpha_{b_{k+1}}^{b_{k}}
αbk+1bk、
β
b
k
+
1
b
k
\beta_{b_{k+1}}^{b_{k}}
βbk+1bk、
γ
b
k
+
1
b
k
\gamma_{b_{k+1}}^{b_{k}}
γbk+1bk按其偏置的一阶近似来调整,否则就进行重新传递。
【论文式[12]】(这部分只是抛出一个概念,后面会讲为什么这样写)
α
b
k
+
1
b
k
≈
α
^
b
k
+
1
b
k
+
J
b
a
α
δ
b
a
+
J
b
w
α
δ
b
w
\alpha_{b_{k+1}}^{b_{k}}\approx\hat{\alpha}_{b_{k+1}}^{b_{k}}+J_{b_{a}}^{\alpha}\deltab_{a}+J_{b_{w}}^{\alpha}\deltab_{w}
αbk+1bk≈α^bk+1bk+Jbaαδba+Jbwαδbw
β
b
k
+
1
b
k
≈
β
^
b
k
+
1
b
k
+
J
b
a
β
δ
b
a
+
J
b
w
β
δ
b
w
\beta_{b_{k+1}}^{b_{k}}\approx\hat{\beta}_{b_{k+1}}^{b_{k}}+J_{b_{a}}^{\beta}\deltab_{a}+J_{b_{w}}^{\beta}\deltab_{w}
βbk+1bk≈β^bk+1bk+Jbaβδba+Jbwβδbw
γ
b
k
+
1
b
k
≈
γ
^
b
k
+
1
b
k
[
1
1
2
J
b
w
γ
δ
b
w
]
\gamma_{b_{k+1}}^{b_{k}}\approx\hat{\gamma}_{b_{k+1}}^{b_{k}}\begin{bmatrix}1\\\frac{1}{2}J_{b_{w}}^{\gamma}\deltab_{w}\end{bmatrix}
γbk+1bk≈γ^bk+1bk[121Jbwγδbw]
4.2离散形式下的IMU运动模型
4.2.1两帧之间PVQ增量的欧拉法离散形式
欧拉法给出第i时刻与第i+1时刻的预积分量估计值的关系:【论文式[7]】
α
^
i
+
1
b
k
=
α
^
i
b
k
+
β
^
i
b
k
δ
t
+
1
2
R
(
γ
^
i
b
k
)
(
a
^
i
−
b
a
i
)
δ
t
2
\hat{\alpha}^{b_{k}}_{i+1}=\hat{\alpha}^{b_{k}}_{i}+\hat{\beta}_{i}^{b_{k}}\deltat+\frac{1}{2}R(\hat{\gamma}_{i}^{b_{k}})(\hat{a}_{i}-b_{ai})\deltat^{2}
α^i+1bk=α^ibk+β^ibkδt+21R(γ^ibk)(a^i−bai)δt2
β
^
i
+
1
b
k
=
β
^
i
b
k
+
1
2
R
(
γ
^
i
b
k
)
(
a
^
i
−
b
a
i
)
δ
t
\hat{\beta}^{b_{k}}_{i+1}=\hat{\beta}^{b_{k}}_{i}+\frac{1}{2}R(\hat{\gamma}_{i}^{b_{k}})(\hat{a}_{i}-b_{ai})\deltat
β^i+1bk=β^ibk+21R(γ^ibk)(a^i−bai)δt
γ
^
i
+
1
b
k
=
γ
^
i
b
k
⊗
γ
^
i
+
1
i
=
γ
^
i
b
k
⊗
[
1
1
2
(
w
^
i
−
b
w
i
)
δ
t
]
\hat{\gamma}^{b_{k}}_{i+1}=\hat{\gamma}^{b_{k}}_{i}\otimes\hat{\gamma}^{i}_{i+1}=\hat{\gamma}^{b_{k}}_{i}\otimes\begin{bmatrix}1\\\frac{1}{2}(\hat{w}_{i}-b_{wi})\deltat\end{bmatrix}
γ^i+1bk=γ^ibk⊗γ^i+1i=γ^ibk⊗[121(w^i−bwi)δt]其中
i
i
i是
[
t
k
,
t
k
+
1
]
[t_{k},t_{k+1}]
[tk,tk+1]之间的离散时间
4.2.2两帧之间PVQ增量的中值法离散形式
代码中采用的基于中值法的IMU预积分公式,这在Estimator::processIMU()函数中的IntegrationBase::push_back()函数中得以实现,注意这里积分出来的是前后两帧之间的IMU增量信息。
α
^
i
+
1
b
k
=
α
^
i
b
k
+
β
^
i
b
k
δ
t
+
1
2
a
^
ˉ
i
δ
t
2
\hat{\alpha}^{b_{k}}_{i+1}=\hat{\alpha}^{b_{k}}_{i}+\hat{\beta}_{i}^{b_{k}}\deltat+\frac{1}{2}\bar{\hat{a}}_{i}\deltat^{2}
α^i+1bk=α^ibk+β^ibkδt+21a^ˉiδt2
β
^
i
+
1
b
k
=
β
^
i
b
k
+
a
^
ˉ
i
δ
t
\hat{\beta}^{b_{k}}_{i+1}=\hat{\beta}^{b_{k}}_{i}+\bar{\hat{a}}_{i}\deltat
β^i+1bk=β^ibk+a^ˉiδt
γ
^
i
+
1
b
k
=
γ
^
i
b
k
⊗
γ
^
i
+
1
i
=
γ
^
i
b
k
⊗
[
1
1
2
w
^
ˉ
i
δ
t
]
\hat{\gamma}^{b_{k}}_{i+1}=\hat{\gamma}^{b_{k}}_{i}\otimes\hat{\gamma}^{i}_{i+1}=\hat{\gamma}^{b_{k}}_{i}\otimes\begin{bmatrix}1\\\frac{1}{2}\bar{\hat{w}}_{i}\deltat\end{bmatrix}
γ^i+1bk=γ^ibk⊗γ^i+1i=γ^ibk⊗[121w^ˉiδt]其中
a
^
ˉ
t
=
1
2
[
γ
^
i
b
k
(
a
^
i
−
b
i
)
+
γ
^
i
+
1
b
k
(
a
^
i
+
1
−
b
a
i
)
]
\bar{\hat{a}}_{t}=\frac{1}{2}[\hat{\gamma}^{b_{k}}_{i}(\hat{a}_{i}-b_{i})+\hat{\gamma}^{b_{k}}_{i+1}(\hat{a}_{i+1}-b_{ai})]
a^ˉt=21[γ^ibk(a^i−bi)+γ^i+1bk(a^i+1−bai)]
w
^
ˉ
i
=
1
2
(
w
^
i
+
w
^
i
+
1
)
−
b
w
i
\bar{\hat{w}}_{i}=\frac{1}{2}(\hatw_{i}+\hatw_{i+1})-b_{wi}
w^ˉi=21(w^i+w^i+1)−bwi初始状态下
α
b
k
b
k
\alpha_{b_{k}}^{b_{k}}
αbkbk、
β
b
k
b
k
\beta_{b_{k}}^{b_{k}}
βbkbk为0,
γ
b
k
b
k
\gamma_{b_{k}}^{b_{k}}
γbkbk为单位四元素,
n
a
n_{a}
na、
n
w
n_{w}
nw被视为0,
i
i
i为在
[
k
,
k
+
1
]
[k,k+1]
[k,k+1]中IMU测量值的某一时刻,
δ
t
\deltat
δt为IMU测量值
i
i
i和
i
+
1
i+1
i+1之间的时间间隔。
5.非线性方程的误差递推方程
5.1一段时间内多个IMU数据积分形成的预积分量的协方差计算
一个IMU数据作为测量值的噪声方差我们能够标定。
现在,一段时间内多个IMU数据积分形成的预积分量的方差呢?
要推导预积分量的协方差,需要知道IMU噪声和预积分量之间的线性递推关系。
假设已知了相邻时刻误差的线性传递方程:
η
i
k
+
1
=
F
k
η
i
k
+
G
k
n
k
\eta_{ik+1}=F_{k}\eta_{ik}+G_{k}n_{k}
ηik+1=Fkηik+Gknk其中
η
i
k
\eta_{ik}
ηik为状态量误差且
η
i
k
=
[
δ
θ
i
k
,
δ
v
i
k
,
δ
p
i
k
]
\eta_{ik}=[\delta\theta_{ik},\deltav_{ik},\deltap_{ik}]
ηik=[δθik,δvik,δpik],
n
k
n_{k}
nk为测量噪声且
n
k
=
[
n
k
g
,
n
k
g
]
n_{k}=[n_{k}^{g},n_{k}^{g}]
nk=[nkg,nkg]。
可以看出误差的传递由倆部分组成:当前时刻的误差传递给下一时刻,当前时刻测量噪声传递给下一时刻。
5.2非线性方程的误差递推方程推导
通常对于状态量之间的递推关系是非线性的方程如
x
k
+
1
=
f
(
x
k
,
u
k
)
x_{k+1}=f(x_{k},u_{k})
xk+1=f(xk,uk),其中状态量
x
x
x、
u
u
u为系统的输入量。
可以用俩种方法来推导状态误差传递的线性递推关系:
一种是基于一阶泰勒展开的误差递推方程一种是基于误差随时间变化的递推方程(论文是基于误差随时间变化来推导的)
5.2.1基于一阶泰勒展开的误差递推方程
令状态量为
x
^
=
x
+
δ
x
\hat{x}=x+\deltax
x^=x+δx,其中
x
^
\hat{x}
x^表示计算得到的值,带有误差,真值为
x
x
x,误差为
δ
x
\deltax
δx。
另外,输入量
u
u
u的噪声为
n
n
n。
应用一阶泰勒展开(EKF的协方差预测也是用了泰勒展开),非线性系统
x
^
k
+
1
=
f
(
x
^
k
,
u
^
k
)
\hat{x}_{k+1}=f(\hat{x}_{k},\hat{u}_{k})
x^k+1=f(x^k,u^k)的状态误差的线性递推关系为:
δ
x
k
+
1
=
F
δ
x
k
+
G
n
k
\deltax_{k+1}=F\deltax_{k}+Gn_{k}
δxk+1=Fδxk+Gnk其中,F是状态量
x
k
+
1
x_{k+1}
xk+1对状态量
x
k
x_{k}
xk的雅克比矩阵,G是状态量
x
k
+
1
x_{k+1}
xk+1对输入量
u
k
u_{k}
uk的雅克比矩阵。
证明:对非线性状态方程进行一阶泰勒展开有:
x
^
k
+
1
=
f
(
x
^
k
,
u
^
k
)
\hat{x}_{k+1}=f(\hat{x}_{k},\hat{u}_{k})
x^k+1=f(x^k,u^k)
x
k
+
1
+
δ
x
k
+
1
=
f
(
x
k
+
δ
x
k
,
u
k
+
n
k
)
x_{k+1}+\deltax_{k+1}=f(x_{k}+\deltax_{k},u_{k}+n_{k})
xk+1+δxk+1=f(xk+δxk,uk+nk)
x
k
+
1
‾
+
δ
x
k
+
1
=
f
(
x
k
,
u
k
)
‾
+
F
δ
x
k
+
G
n
k
(3)
\underline{x_{k+1}}+\deltax_{k+1}=\underline{f(x_{k},u_{k})}+F\deltax_{k}+Gn_{k}\tag{3}
xk+1+δxk+1=f(xk,uk)+Fδxk+Gnk(3)
5.2.2基于误差随时间变化(误差导数)的递推方程
如果我们能够推导状态误差随时间变化的导数关系,如:
δ
˙
x
=
A
δ
x
+
B
n
\dot{\delta}_{x}=A\deltax+Bn
δ˙x=Aδx+Bn则误差状态的传递方程为:
δ
x
k
+
1
=
δ
x
k
+
δ
˙
x
k
Δ
t
\deltax_{k+1}=\deltax_{k}+\dot{\delta}x_{k}\Deltat
δxk+1=δxk+δ˙xkΔt
→
δ
x
k
+
1
=
(
I
+
A
Δ
t
)
δ
x
k
+
B
Δ
t
n
k
(4)
\rightarrow\deltax_{k+1}=(I+A\Deltat)\deltax_{k}+B\Deltatn_{k}\tag{4}
→δxk+1=(I+AΔt)δxk+BΔtnk(4)
由公式
(
3
)
(
4
)
(3)(4)
(3)(4)得
F
=
I
+
A
Δ
t
F=I+A\Deltat
F=I+AΔt,
G
=
B
Δ
t
G=B\Deltat
G=BΔt
5.2.3基于泰勒展开和误差随时间变化的方法对比
基于一阶泰勒展开的误差递推方程不香吗,为什么会用误差随时间的变化来进行误差递推方程的构建呢?
在VIO系统中,我们已经知道了状态的导数和状态之间的转移矩阵,如我们已经知道速度和状态量之间的关系:
v
˙
w
=
R
b
w
a
b
+
g
w
\dot{v}^{w}=R^{w}_{b}a^{b}+g^{w}
v˙w=Rbwab+gw那我们就可以推导速度的误差导数和状态误差之间的关系,在每一项上都加上各自的误差就有:
v
˙
w
+
δ
˙
v
w
=
R
b
w
e
x
p
(
[
δ
θ
]
×
)
(
a
b
+
δ
a
b
)
+
g
+
δ
g
\dot{v}^{w}+\dot{\delta}v^{w}=R_{b}^{w}exp([\delta\theta]_{\times})(a^{b}+\deltaa^{b})+g+\deltag
v˙w+δ˙vw=Rbwexp([δθ]×)(ab+δab)+g+δg
v
˙
w
+
δ
˙
v
w
=
R
b
w
(
I
+
[
δ
θ
]
×
)
(
a
b
+
δ
a
b
)
+
g
+
δ
g
\dot{v}^{w}+\dot{\delta}v^{w}=R_{b}^{w}(I+[\delta\theta]_{\times})(a^{b}+\deltaa^{b})+g+\deltag
v˙w+δ˙vw=Rbw(I+[δθ]×)(ab+δab)+g+δg
v
˙
‾
w
+
δ
˙
v
w
=
R
b
w
a
b
+
g
‾
+
R
b
w
δ
a
b
+
R
b
w
[
δ
θ
]
×
(
a
b
+
δ
a
b
)
+
δ
g
\underline{\dot{v}}^{w}+\dot{\delta}v^{w}=\underline{R_{b}^{w}a^{b}+g}+R_{b}^{w}\deltaa^{b}+R_{b}^{w}[\delta\theta]_{\times}(a^{b}+\deltaa^{b})+\deltag
v˙w+δ˙vw=Rbwab+g+Rbwδab+Rbw[δθ]×(ab+δab)+δg
δ
˙
v
w
=
R
b
w
δ
a
b
+
R
b
w
[
δ
θ
]
×
(
a
b
+
δ
a
b
)
+
δ
g
\dot{\delta}{v}^{w}=R_{b}^{w}\deltaa^{b}+R_{b}^{w}[\delta\theta]_{\times}(a^{b}+\deltaa^{b})+\deltag
δ˙vw=Rbwδab+Rbw[δθ]×(ab+δab)+δg
δ
˙
v
w
=
R
b
w
δ
a
b
−
R
b
w
[
a
b
]
×
δ
θ
+
δ
g
(5)
\dot{\delta}v^{w}=R_{b}^{w}\deltaa^{b}-R_{b}^{w}[a^{b}]_{\times}\delta\theta+\deltag\tag{5}
δ˙vw=Rbwδab−Rbw[ab]×δθ+δg(5)由此就能以此类推,轻易写出整个A和B其他方程了。
考虑到公式的编辑篇幅,为了对一些求导公式进行简化,对求导公式进行了简化(下同):
∂
a
w
∂
δ
θ
=
l
i
m
δ
θ
→
0
R
b
w
e
x
p
(
[
δ
θ
]
×
)
a
b
−
R
b
w
a
b
δ
θ
→
∂
a
w
∂
δ
θ
=
∂
R
b
w
e
x
p
(
[
δ
θ
]
×
)
a
b
∂
δ
θ
\frac{\partiala^{w}}{\partial\delta\theta}=\underset{\delta\theta\rightarrow0}{lim}\frac{R_{b}^{w}exp([\delta\theta]_{\times})a^{b}-R_{b}^{w}a^{b}}{\delta\theta}\rightarrow\frac{\partiala^{w}}{\partial\delta\theta}=\frac{\partialR_{b}^{w}exp([\delta\theta]_{\times})a^{b}}{\partial\delta\theta}
∂δθ∂aw=δθ→0limδθRbwexp([δθ]×)ab−Rbwab→∂δθ∂aw=∂δθ∂Rbwexp([δθ]×)ab
6.PVQ增量误差、协方差及雅克比矩阵的递推方程
6.1连续形式下的误差、协方差及雅克比矩阵的递推方程
在第4部分我们已经完成了IMU预积分测量值的推导,而要将IMU预积分运用到非线性优化中,需要建立线性高斯误差状态递推方程(可参考第5部分),由线性高斯系统的协方差,推到方程协方差矩阵,并纠结对应的雅克比矩阵。
首先由于四元素
γ
t
b
k
\gamma_{t}^{b_{k}}
γtbk被参数化过,我们将其误差定义为围绕其均值的扰动:【论文式(8)】
γ
t
b
k
≈
γ
^
t
b
k
⊗
[
1
1
2
δ
θ
t
b
k
]
\gamma_{t}^{b_{k}}\approx\hat{\gamma}_{t}^{b_{k}}\otimes\begin{bmatrix}1\\\frac{1}{2}\delta\theta_{t}^{b_{k}}\end{bmatrix}
γtbk≈γ^tbk⊗[121δθtbk]由预积分的连续形式建立的运动模型,由公式(2)可得一段时间内IMU构建的预积分量作为测量值,对俩时刻之间的状态量进行约束,可得到其在
k
k
k和
k
+
1
k+1
k+1时刻下的误差项为:【论文公式(24)】
[
δ
α
b
k
+
1
b
k
δ
β
b
k
+
1
b
k
δ
θ
b
k
+
1
b
k
δ
b
a
δ
b
w
]
=
[
R
w
b
k
(
p
b
k
+
1
w
−
p
b
k
w
−
v
k
w
Δ
t
+
1
2
g
w
Δ
t
2
)
−
α
b
k
+
1
b
k
R
w
b
k
(
v
b
k
+
1
w
−
v
b
k
w
+
g
w
Δ
t
)
−
β
b
k
+
1
b
k
2
[
q
w
b
k
⊗
q
b
k
+
1
w
⊗
γ
b
k
b
k
+
1
]
x
y
z
b
a
b
k
+
1
−
b
a
b
k
b
w
b
k
+
1
−
b
w
b
k
]
\begin{bmatrix}\delta\alpha_{b_{k+1}}^{b_{k}}\\\delta\beta_{b_{k+1}}^{b_{k}}\\\delta\theta_{b_{k+1}}^{b_{k}}\\\deltab_{a}\\\deltab_{w}\end{bmatrix}=\begin{bmatrix}R^{b_{k}}_{w}(p^{w}_{b_{k+1}}-p_{b_{k}}^{w}-v_{k}^{w}\Deltat+\frac{1}{2}g^{w}\Deltat^{2})-\alpha_{b_{k+1}}^{b_{k}}\\R^{b_{k}}_{w}(v_{b_{k+1}}^{w}-v_{b_{k}}^{w}+g^{w}\Deltat)-\beta^{b_{k}}_{b_{k+1}}\\2[q_{w}^{b_{k}}\otimesq^{w}_{b_{k+1}}\otimes\gamma_{b_{k}}^{b_{k+1}}]_{xyz}\\b_{ab_{k+1}}-b_{ab_{k}}\\b_{wb_{k+1}}-b_{wb_{k}}\end{bmatrix}
⎣⎢⎢⎢⎢⎢⎡δαbk+1bkδβbk+1bkδθbk+1bkδbaδbw⎦⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎡Rwbk(pbk+1w−pbkw−vkwΔt+21gwΔt2)−αbk+1bkRwbk(vbk+1w−vbkw+gwΔt)−βbk+1bk2[qwbk⊗qbk+1w⊗γbkbk+1]xyzbabk+1−babkbwbk+1−bwbk⎦⎥⎥⎥⎥⎥⎤
我们给出在
t
t
t时刻误差项的线性化递推方程为:【论文式(9)】:
[
δ
α
˙
t
b
k
δ
β
˙
t
b
k
δ
θ
˙
t
b
k
δ
b
˙
a
t
δ
b
˙
w
t
]
=
[
0
I
0
0
0
0
0
−
R
t
b
k
[
a
^
t
−
b
a
t
]
×
−
R
t
b
k
0
0
0
−
[
w
^
−
b
w
t
]
×
0
−
I
0
0
0
0
0
0
0
0
0
0
]
[
δ
α
t
b
k
δ
β
t
b
k
δ
θ
t
b
k
δ
b
a
t
δ
b
w
t
]
+
[
0
0
0
0
−
R
t
b
k
0
0
0
0
−
I
0
0
0
0
I
0
0
0
0
I
]
[
n
a
n
w
n
b
a
n
b
w
]
=
F
t
δ
z
t
b
k
+
G
t
n
t
(6)
\begin{bmatrix}\delta\dot{\alpha}_{t}^{b_{k}}\\\delta\dot{\beta}_{t}^{b_{k}}\\\delta\dot{\theta}_{t}^{b_{k}}\\\delta\dot{b}_{at}\\\delta\dot{b}_{wt}\end{bmatrix}=\begin{bmatrix}0&I&0&0&0\\0&0&-R_{t}^{b_{k}}[\hat{a}_{t}-b_{at}]_{\times}&-R_{t}^{b_{k}}&0\\0&0&-[\hat{w}-b_{wt}]_{\times}&0&-I\\0&0&0&0&0\\0&0&0&0&0\end{bmatrix}\begin{bmatrix}\delta\alpha_{t}^{b_{k}}\\\delta\beta_{t}^{b_{k}}\\\delta\theta_{t}^{b_{k}}\\\deltab_{at}\\\deltab_{wt}\end{bmatrix}+\begin{bmatrix}0&0&0&0\\-R_{t}^{b_{k}}&0&0&0\\0&-I&0&0\\0&0&I&0\\0&0&0&I\end{bmatrix}\begin{bmatrix}n_{a}\\n_{w}\\n_{ba}\\n_{bw}\end{bmatrix}\tag{6}\\=F_{t}\deltaz_{t}^{b_{k}}+G_{t}n_{t}
⎣⎢⎢⎢⎢⎡δα˙tbkδβ˙tbkδθ˙tbkδb˙atδb˙wt⎦⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎡00000I00000−Rtbk[a^t−bat]×−[w^−bwt]×000−Rtbk00000−I00⎦⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎡δαtbkδβtbkδθtbkδbatδbwt⎦⎥⎥⎥⎥⎤+⎣⎢⎢⎢⎢⎡0−Rtbk00000−I00000I00000I⎦⎥⎥⎥⎥⎤⎣⎢⎢⎡nanwnbanbw⎦⎥⎥⎤=Ftδztbk+Gtnt(6)可以简写为:
δ
z
˙
t
b
k
=
F
t
δ
z
t
b
k
+
G
t
n
t
\delta\dot{z}_{t}^{b_{k}}=F_{t}\deltaz_{t}^{b_{k}}+G_{t}n_{t}
δz˙tbk=Ftδztbk+Gtnt其中
F
t
F_{t}
Ft是15x15,
G
t
G_{t}
Gt是12x12,
δ
z
t
b
k
\deltaz_{t}^{b_{k}}
δztbk是15x1,
n
t
n_{t}
nt是12x1。
根据导数的定义有:
δ
z
˙
t
b
k
=
l
i
m
δ
t
→
0
δ
z
t
+
δ
t
b
k
−
δ
z
t
b
k
δ
t
\delta\dot{z}_{t}^{b_{k}}=\underset{\deltat\rightarrow0}{lim}\frac{\deltaz_{t+\deltat}^{b_{k}}-\deltaz_{t}^{b_{k}}}{\deltat}
δz˙tbk=δt→0limδtδzt+δtbk−δztbk即:
δ
z
t
+
δ
t
b
k
=
δ
z
t
b
k
+
δ
z
˙
t
b
k
δ
t
=
δ
z
t
b
k
+
(
F
t
δ
z
t
b
k
+
G
t
n
t
)
δ
t
=
(
I
+
F
t
δ
t
)
δ
z
t
b
k
+
G
t
δ
t
n
t
\deltaz_{t+\deltat}^{b_{k}}=\deltaz_{t}^{b_{k}}+\delta\dot{z}_{t}^{b_{k}}\deltat=\deltaz_{t}^{b_{k}}+(F_{t}\deltaz_{t}^{b_{k}}+G_{t}n_{t})\deltat=(I+F_{t}\deltat)\deltaz_{t}^{b_{k}}+G_{t}\deltatn_{t}
δzt+δtbk=δztbk+δz˙tbkδt=δztbk+(Ftδztbk+Gtnt)δt=(I+Ftδt)δztbk+Gtδtnt这个形式与公式(4)不谋而合,一般连续形式下的微小时间用
δ
t
\deltat
δt表示,离散形式下的时间区间用
Δ
t
\Deltat
Δt表示,其实都是一个意思。
进一步的,令
F
=
(
I
+
F
t
δ
t
)
,
V
=
G
t
δ
t
F=(I+F_{t}\deltat),V=G_{t}\deltat
F=(I+Ftδt),V=Gtδt,有:
δ
z
t
+
δ
t
b
k
=
F
δ
z
t
b
k
+
V
n
t
(7)
\deltaz_{t+\deltat}^{b_{k}}=F\deltaz_{t}^{b_{k}}+Vn_{t}\tag{7}
δzt+δtbk=Fδztbk+Vnt(7)由此我们已经知道了该非线性系统的线性化递推方程。
至于下一时刻的协方差,首先对于一个高斯分布的变量线性变换后得到新的变量的协方差计算,我们有以下推导:已知一个变量
y
=
A
x
y=Ax
y=Ax,
x
∼
(
0
,
Σ
x
)
x\sim(0,\Sigma_{x})
x∼(0,Σx),则有
Σ
y
=
A
Σ
x
A
T
\Sigma_{y}=A\Sigma_{x}A^{T}
Σy=AΣxAT,即:
Σ
y
=
E
(
(
A
x
)
(
A
x
)
T
)
=
E
(
A
x
x
T
A
T
)
=
A
Σ
x
A
T
\Sigma_{y}=E((Ax)(Ax)^{T})=E(Axx^{T}A^{T})=A\Sigma_{x}A^{T}
Σy=E((Ax)(Ax)T)=E(AxxTAT)=AΣxAT
得到了递推方程
(
7
)
(7)
(7)之后,根据高斯分布线性变换协方差的传递规律,我们可以预测下一时刻的协方差:【论文式(10)】
P
t
+
δ
t
b
k
=
(
I
+
F
t
δ
t
)
P
t
b
k
(
I
+
F
t
δ
t
)
T
+
(
G
δ
t
)
Q
(
G
t
δ
t
)
T
P_{t+\deltat}^{b_{k}}=(I+F_{t}\deltat)P_{t}^{b_{k}}(I+F_{t}\deltat)^{T}+(G\deltat)Q(G_{t}\deltat)^{T}
Pt+δtbk=(I+Ftδt)Ptbk(I+Ftδt)T+(Gδt)Q(Gtδt)T其中初始协方差
P
b
k
b
k
=
0
P_{b_{k}}^{b_{k}}=0
Pbkbk=0,
Q
Q
Q代表噪声项的对角协方差矩阵:
Q
12
×
12
=
[
σ
a
2
0
0
0
0
σ
w
2
0
0
0
0
σ
b
a
2
0
0
0
0
σ
b
w
2
]
Q^{12\times12}=\begin{bmatrix}\sigma_{a}^{2}&0&0&0\\0&\sigma_{w}^{2}&0&0\\0&0&\sigma_{ba}^{2}&0\\0&0&0&\sigma_{bw}^{2}\end{bmatrix}
Q12×12=⎣⎢⎢⎡σa20000σw20000σba20000σbw2⎦⎥⎥⎤
同时也可以根据递推方程得到对应的雅克比矩阵迭代公式:【论文式[11]】
J
t
+
δ
t
=
(
I
+
F
t
δ
t
)
J
t
J_{t+\deltat}=(I+F_{t}\deltat)J_{t}
Jt+δt=(I+Ftδt)Jt其中初始雅克比矩阵为:
J
b
k
=
I
J_{b_{k}}=I
Jbk=I。
关于公式
(
6
)
(6)
(6)的推导:先把之前的预积分再写一遍:
α
b
k
+
1
b
k
=
∫
∫
t
∈
[
t
k
,
t
k
+
1
]
R
t
b
k
(
a
^
t
−
b
a
t
−
n
a
)
d
t
2
\alpha^{b_{k}}_{b_{k+1}}=\int\int_{t\in[t_{k},t_{k+1}]}R_{t}^{b_{k}}(\hat{a}_{t}-b_{at}-n_{a})dt^{2}
αbk+1bk=∫∫t∈[tk,tk+1]Rtbk(a^t−bat−na)dt2
β
b
k
+
1
b
k
=
∫
t
∈
[
t
k
,
t
k
+
1
]
R
t
b
k
(
a
^
t
−
b
a
t
−
n
a
)
d
t
\beta_{b_{k+1}}^{b_{k}}=\int_{t\in[t_{k},t_{k+1}]}R_{t}^{b_{k}}(\hat{a}_{t}-b_{at}-n_{a})dt
βbk+1bk=∫t∈[tk,tk+1]Rtbk(a^t−bat−na)dt
γ
b
k
+
1
=
∫
t
∈
[
t
k
,
t
k
+
1
]
1
2
Ω
(
w
^
t
−
b
w
t
−
n
w
)
γ
t
b
k
d
t
\gamma_{b_{k+1}}=\int_{t\in[t_{k},t_{k+1}]}\frac{1}{2}\Omega(\hat{w}_{t}-b_{wt}-n_{w})\gamma_{t}^{b_{k}}dt
γbk+1=∫t∈[tk,tk+1]21Ω(w^t−bwt−nw)γtbkdt1)首先对于
δ
α
˙
t
b
k
\delta\dot{\alpha}_{t}^{b_{k}}
δα˙tbk、
δ
b
˙
a
t
\delta\dot{b}_{at}
δb˙at、
δ
b
˙
w
t
\delta\dot{b}_{wt}
δb˙wt,显然根据定义有:
δ
α
˙
t
b
k
=
α
˙
^
t
b
k
−
α
˙
t
b
k
=
β
^
t
b
k
−
β
t
b
k
=
δ
β
t
b
k
\delta\dot{\alpha}_{t}^{b_{k}}=\hat{\dot{\alpha}}_{t}^{b_{k}}-\dot{\alpha}_{t}^{b_{k}}=\hat{\beta}_{t}^{b_{k}}-\beta_{t}^{b_{k}}=\delta\beta_{t}^{b_{k}}
δα˙tbk=α˙^tbk−α˙tbk=β^tbk−βtbk=δβtbk
δ
b
˙
a
t
=
b
˙
a
t
−
0
=
n
b
a
\delta\dot{b}_{at}=\dot{b}_{at}-0=n_{b_{a}}
δb˙at=b˙at−0=nba
δ
b
˙
w
t
=
b
˙
w
t
−
0
=
n
b
w
\delta\dot{b}_{wt}=\dot{b}_{wt}-0=n_{b_{w}}
δb˙wt=b˙wt−0=nbw其中上标"^"(如
α
˙
^
t
b
k
\hat{\dot{\alpha}}_{t}^{b_{k}}
α˙^tbk)代表真实测量值,包含了噪声,
α
˙
t
b
k
\dot{\alpha}_{t}^{b_{k}}
α˙tbk代表无噪声的理论值,偏置被建立成随机游走模型,其导数服从0均值的高斯分布。
2)然后是
δ
β
˙
t
b
k
\delta\dot{\beta}_{t}^{b_{k}}
δβ˙tbk:
β
˙
t
b
k
\dot{\beta}_{t}^{b_{k}}
β˙tbk为理论值,即不考虑存在噪声
n
a
n_{a}
na和
b
b
a
b_{ba}
bba时应为:
β
˙
t
b
k
=
R
t
b
k
(
a
^
t
−
b
a
t
)
\dot{\beta}_{t}^{b_{k}}=R_{t}^{b_{k}}(\hata_{t}-b_{at})
β˙tbk=Rtbk(a^t−bat)
β
˙
t
b
k
\dot\beta_{t}^{b_{k}}
β˙tbk的实际测量值,即考虑噪声时应为:
β
˙
^
t
b
k
=
R
^
t
b
k
(
a
^
t
−
b
^
a
t
−
n
a
)
=
R
t
b
k
e
x
p
(
[
δ
θ
]
×
)
(
a
^
t
−
n
a
−
b
a
t
−
δ
b
a
t
)
=
R
t
b
k
(
1
+
[
δ
θ
]
×
)
(
a
^
t
−
n
a
−
b
a
t
−
δ
b
a
t
)
=
R
t
b
k
(
a
^
t
−
n
a
−
b
a
t
−
δ
b
a
t
+
[
δ
θ
]
×
(
a
^
t
−
b
a
t
)
)
=
R
t
b
k
(
a
^
t
−
n
a
−
b
a
t
−
δ
b
a
t
−
[
a
^
t
−
b
a
t
]
×
δ
θ
)
(8)
\hat{\dot{\beta}}_{t}^{b_{k}}=\hat{R}_{t}^{b_{k}}(\hata_{t}-\hatb_{at}-n_{a})=R_{t}^{b_{k}}exp([\delta\theta]_{\times})(\hata_{t}-n_{a}-b_{at}-\deltab_{at})\\=R_{t}^{b_{k}}(1+[\delta\theta]_{\times})(\hata_{t}-n_{a}-b_{at}-\deltab_{at})\\=R_{t}^{b_{k}}(\hata_{t}-n_{a}-b_{at}-\deltab_{at}+[\delta\theta]_{\times}(\hata_{t}-b_{at}))\\=R_{t}^{b_{k}}(\hata_{t}-n_{a}-b_{at}-\deltab_{at}-[\hata_{t}-b_{at}]_{\times}\delta\theta)\tag{8}
β˙^tbk=R^tbk(a^t−b^at−na)=Rtbkexp([δθ]×)(a^t−na−bat−δbat)=Rtbk(1+[δθ]×)(a^t−na−bat−δbat)=Rtbk(a^t−na−bat−δbat+[δθ]×(a^t−bat))=Rtbk(a^t−na−bat−δbat−[a^t−bat]×δθ)(8)则有
δ
β
˙
t
b
k
=
β
˙
^
t
b
k
−
β
˙
t
b
k
=
−
R
t
b
k
[
a
^
t
−
b
a
t
]
×
δ
θ
−
R
t
b
k
δ
b
a
t
−
R
t
b
k
n
a
\delta\dot{\beta}_{t}^{b_{k}}=\hat{\dot{\beta}}_{t}^{b_{k}}-\dot{\beta}_{t}^{b_{k}}=-R_{t}^{b_{k}}[\hata_{t}-b_{at}]_{\times}\delta\theta-R_{t}^{b_{k}}\deltab_{at}-R_{t}^{b_{k}}n_{a}
δβ˙tbk=β˙^tbk−β˙tbk=−Rtbk[a^t−bat]×δθ−Rtbkδbat−Rtbkna其中公式
(
8
)
(8)
(8)的推导原因可借鉴公式
(
5
)
(5)
(5),证毕!3)接下来是
δ
θ
˙
t
b
k
\delta\dot{\theta}_{t}^{b_{k}}
δθ˙tbk,先推个
δ
q
˙
t
b
k
\delta\dot{q}_{t}^{b_{k}}
δq˙tbk
q
˙
t
b
k
\dotq_{t}^{b_{k}}
q˙tbk的理论值,即考虑噪声
n
w
n_{w}
nw和
n
b
w
n_{b_{w}}
nbw:
q
˙
t
b
k
=
1
2
q
t
b
k
⊗
[
(
w
^
t
−
b
w
t
)
0
]
=
1
2
Ω
(
w
^
t
−
b
w
t
)
q
t
b
k
\dotq_{t}^{b_{k}}=\frac{1}{2}q_{t}^{b_{k}}\otimes\begin{bmatrix}(\hat{w}_{t}-b_{wt})\\0\end{bmatrix}=\frac{1}{2}\Omega(\hat{w}_{t}-b_{wt})q_{t}^{b_{k}}
q˙tbk=21qtbk⊗[(w^t−bwt)0]=21Ω(w^t−bwt)qtbk
q
˙
t
b
k
\dotq_{t}^{b_{k}}
q˙tbk的真实测量值:
q
˙
t
^
b
k
=
1
2
q
t
^
b
k
⊗
[
(
w
^
t
−
n
w
−
b
w
t
−
δ
b
w
t
)
0
]
=
1
2
q
t
b
k
⊗
δ
q
t
b
k
⊗
[
(
w
^
t
−
n
w
−
b
w
t
−
δ
b
w
t
)
0
]
\hat{\dotq_{t}}^{b_{k}}=\frac{1}{2}\hat{q_{t}}^{b_{k}}\otimes\begin{bmatrix}(\hat{w}_{t}-n_{w}-b_{wt}-\deltab_{wt})\\0\end{bmatrix}\\=\frac{1}{2}q_{t}^{b_{k}}\otimes\deltaq_{t}^{b_{k}}\otimes\begin{bmatrix}(\hat{w}_{t}-n_{w}-b_{wt}-\deltab_{wt})\\0\end{bmatrix}
q˙t^bk=21qt^bk⊗[(w^t−nw−bwt−δbwt)0]=21qtbk⊗δqtbk⊗[(w^t−nw−bwt−δbwt)0]根据导数性质,有:
q
˙
t
^
b
k
=
(
q
t
b
k
⊗
˙
δ
q
t
b
k
)
=
q
˙
t
b
k
⊗
δ
q
t
b
k
+
q
t
b
k
⊗
δ
q
˙
t
b
k
=
1
2
q
t
b
k
⊗
[
(
w
^
t
−
b
w
t
)
0
]
⊗
δ
q
t
b
k
+
q
t
b
k
⊗
δ
q
˙
t
b
k
\hat{\dotq_{t}}^{b_{k}}=(q_{t}^{b_{k}}\dot\otimes\deltaq_{t}^{b_{k}})=\dotq_{t}^{b_{k}}\otimes\deltaq_{t}^{b_{k}}+q_{t}^{b_{k}}\otimes\delta\dotq_{t}^{b_{k}}=\frac{1}{2}q_{t}^{b_{k}}\otimes\begin{bmatrix}(\hat{w}_{t}-b_{wt})\\0\end{bmatrix}\otimes\deltaq_{t}^{b_{k}}+q_{t}^{b_{k}}\otimes\delta\dotq_{t}^{b_{k}}
q˙t^bk=(qtbk⊗˙δqtbk)=q˙tbk⊗δqtbk+qtbk⊗δq˙tbk=21qtbk⊗[(w^t−bwt)0]⊗δqtbk+qtbk⊗δq˙tbk即得到等式:
1
2
q
t
b
k
⊗
δ
q
t
b
k
⊗
[
(
w
^
t
−
n
w
−
b
w
t
−
δ
b
w
t
)
0
]
=
1
2
q
t
b
k
⊗
[
(
w
^
t
−
b
w
t
)
0
]
⊗
δ
q
t
b
k
+
q
t
b
k
⊗
δ
q
˙
t
b
k
\frac{1}{2}q_{t}^{b_{k}}\otimes\deltaq_{t}^{b_{k}}\otimes\begin{bmatrix}(\hat{w}_{t}-n_{w}-b_{wt}-\deltab_{wt})\\0\end{bmatrix}=\frac{1}{2}q_{t}^{b_{k}}\otimes\begin{bmatrix}(\hat{w}_{t}-b_{wt})\\0\end{bmatrix}\otimes\deltaq_{t}^{b_{k}}+q_{t}^{b_{k}}\otimes\delta\dotq_{t}^{b_{k}}
21qtbk⊗δqtbk⊗[(w^t−nw−bwt−δbwt)0]=21qtbk⊗[(w^t−bwt)0]⊗δqtbk+qtbk⊗δq˙tbk
2
δ
q
˙
t
b
k
=
q
t
b
k
⊗
[
(
w
^
t
−
n
w
−
b
w
t
−
δ
b
w
t
)
0
]
−
[
(
w
^
t
−
b
w
t
)
0
]
⊗
δ
q
t
b
k
2\delta\dotq_{t}^{b_{k}}=q_{t}^{b_{k}}\otimes\begin{bmatrix}(\hat{w}_{t}-n_{w}-b_{wt}-\deltab_{wt})\\0\end{bmatrix}-\begin{bmatrix}(\hat{w}_{t}-b_{wt})\\0\end{bmatrix}\otimes\deltaq_{t}^{b_{k}}
2δq˙tbk=qtbk⊗[(w^t−nw−bwt−δbwt)0]−[(w^t−bwt)0]⊗δqtbk
2
δ
q
˙
t
b
k
=
R
(
[
(
w
^
t
−
n
w
−
b
w
t
−
δ
b
w
t
)
0
]
)
δ
q
t
b
k
−
L
(
[
(
w
^
t
−
b
w
t
)
0
]
)
δ
q
t
b
k
2\delta\dotq_{t}^{b_{k}}=R(\begin{bmatrix}(\hat{w}_{t}-n_{w}-b_{wt}-\deltab_{wt})\\0\end{bmatrix})\deltaq_{t}^{b_{k}}-L(\begin{bmatrix}(\hat{w}_{t}-b_{wt})\\0\end{bmatrix})\deltaq_{t}^{b_{k}}
2δq˙tbk=R([(w^t−nw−bwt−δbwt)0])δqtbk−L([(w^t−bwt)0])δqtbk
2
δ
q
˙
t
b
k
=
[
δ
θ
˙
t
b
k
0
]
=
[
−
[
2
w
^
t
−
2
b
w
t
−
n
w
−
δ
b
w
t
]
×
−
n
w
−
δ
b
w
t
(
n
w
+
δ
b
w
t
)
T
0
]
[
1
2
δ
θ
t
b
k
1
]
2\delta\dotq_{t}^{b_{k}}=\begin{bmatrix}\delta\dot{\theta}_{t}^{b_{k}}\\0\end{bmatrix}=\begin{bmatrix}-[2\hat{w}_{t}-2b_{wt}-n_{w}-\deltab_{wt}]_{\times}&-n_{w}-\deltab_{wt}\\(n_{w}+\deltab_{wt})^{T}&0\end{bmatrix}\begin{bmatrix}\frac{1}{2}\delta\theta_{t}^{b_{k}}\\1\end{bmatrix}
2δq˙tbk=[δθ˙tbk0]=[−[2w^t−2bwt−nw−δbwt]×(nw+δbwt)T−nw−δbwt0][21δθtbk1]
δ
θ
˙
t
b
k
=
−
[
2
w
^
t
−
2
b
w
t
−
n
w
−
δ
b
w
t
]
×
1
2
δ
θ
t
b
k
−
n
w
−
δ
b
w
t
≈
−
[
w
^
t
−
b
w
t
]
×
δ
θ
t
b
k
−
n
w
−
δ
b
w
t
\delta\dot{\theta}_{t}^{b_{k}}=-[2\hat{w}_{t}-2b_{wt}-n_{w}-\deltab_{wt}]_{\times}\frac{1}{2}\delta\theta_{t}^{b_{k}}-n_{w}-\deltab_{wt}\\\approx-[\hat{w}_{t}-b_{wt}]_{\times}\delta\theta_{t}^{b_{k}}-n_{w}-\deltab_{wt}
δθ˙tbk=−[2w^t−2bwt−nw−δbwt]×21δθtbk−nw−δbwt≈−[w^t−bwt]×δθtbk−nw−δbwt证毕!
6.2离散形式下(中值法)增量误差、协方差及雅克比矩阵的递推方程
根据连续形式的增量误差递推方程,同理我们可以推导出离散形式的方程,首先直接给出PVQ增量误差在离散形式下的矩阵形式,这在vins_estimator/src/factor/integration_base.h中的voidintegrationBase::midPointIntegration()函数中得以实现(这里的矩阵V几个变量与VINS中的代码差个负号,体现在前四个噪声项,应该是推导过程中与VINS添加噪声的方式不同,即是减去一个噪声还是加上一个噪声,问题不是特别大)
[
δ
α
k
+
1
δ
θ
k
+
1
δ
β
k
+
1
δ
b
a
k
+
1
δ
b
w
k
+
1
]
=
[
I
f
01
δ
t
f
03
f
04
0
f
11
0
0
−
δ
t
0
f
21
I
f
23
f
24
0
0
0
I
0
0
0
0
0
I
]
[
δ
α
k
δ
θ
k
δ
β
k
δ
b
a
k
δ
b
w
k
]
+
[
v
00
v
01
v
02
v
03
0
0
0
−
1
2
δ
t
0
−
1
2
δ
t
0
0
−
1
2
R
k
δ
t
v
21
−
1
2
R
k
+
1
δ
t
v
23
0
0
0
0
0
0
δ
t
0
0
0
0
0
0
δ
t
]
[
n
a
k
n
w
k
n
a
k
+
1
n
w
k
+
1
n
b
a
n
b
w
]
(9)
\begin{bmatrix}\delta\alpha_{k+1}\\\delta\theta_{k+1}\\\delta\beta_{k+1}\\\deltab_{ak+1}\\\deltab_{wk+1}\end{bmatrix}=\begin{bmatrix}I&f_{01}&\deltat&f_{03}&f_{04}\\0&f_{11}&0&0&-\deltat\\0&f_{21}&I&f_{23}&f_{24}\\0&0&0&I&0\\0&0&0&0&I\end{bmatrix}\begin{bmatrix}\delta\alpha_{k}\\\delta\theta_{k}\\\delta\beta_{k}\\\deltab_{ak}\\\deltab_{wk}\end{bmatrix}+\begin{bmatrix}v_{00}&v_{01}&v_{02}&v_{03}&0&0\\0&-\frac{1}{2}\deltat&0&-\frac{1}{2}\deltat&0&0\\-\frac{1}{2}R_{k}\deltat&v_{21}&-\frac{1}{2}R_{k+1}\deltat&v_{23}&0&0\\0&0&0&0&\deltat&0\\0&0&0&0&0&\deltat\end{bmatrix}\begin{bmatrix}n_{ak}\\n_{wk}\\n_{ak+1}\\n_{wk+1}\\n_{ba}\\n_{bw}\end{bmatrix}\tag{9}
⎣⎢⎢⎢⎢⎡δαk+1δθk+1δβk+1δbak+1δbwk+1⎦⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎡I0000f01f11f2100δt0I00f030f23I0f04−δtf240I⎦⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎡δαkδθkδβkδbakδbwk⎦⎥⎥⎥⎥⎤+⎣⎢⎢⎢⎢⎡v000−21Rkδt00v01−21δtv2100v020−21Rk+1δt00v03−21δtv2300000δt00000δt⎦⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎢⎡naknwknak+1nwk+1nbanbw⎦⎥⎥⎥⎥⎥⎥⎤(9)其中:
f
01
=
1
2
f
21
δ
t
,
f
03
=
−
1
4
(
R
k
+
R
k
+
1
)
δ
t
2
,
f
04
=
1
2
f
24
δ
t
f_{01}=\frac{1}{2}f_{21}\deltat,\f_{03}=-\frac{1}{4}(R_{k}+R_{k+1})\deltat^{2},\f_{04}=\frac{1}{2}f_{24}\delta{t}
f01=21f21δt, f03=−41(Rk+Rk+1)δt2, f04=21f24δt
f
11
=
I
−
[
w
^
k
+
w
^
k
+
1
2
−
b
w
k
]
×
δ
t
f_{11}=I-[\frac{\hat{w}_{k}+\hat{w}_{k+1}}{2}-b_{wk}]_{\times}\deltat
f11=I−[2w^k+w^k+1−bwk]×δt
f
21
=
−
1
2
R
k
[
a
^
k
−
b
a
k
]
×
δ
t
−
1
2
R
k
+
1
[
a
^
k
+
1
−
b
a
k
]
×
(
I
−
[
w
^
k
+
w
^
k
+
1
2
−
b
w
k
]
×
δ
t
)
δ
t
f_{21}=-\frac{1}{2}R_{k}[\hata_{k}-b_{ak}]_{\times}\deltat-\frac{1}{2}R_{k+1}[\hata_{k+1}-b_{ak}]_{\times}(I-[\frac{\hat{w}_{k}+\hat{w}_{k+1}}{2}-b_{wk}]_{\times}\deltat)\deltat
f21=−21Rk[a^k−bak]×δt−21Rk+1[a^k+1−bak]×(I−[2w^k+w^k+1−bwk]×δt)δt
f
23
=
−
1
2
(
R
k
+
R
k
+
1
)
δ
t
,
f
24
=
1
2
R
k
+
1
[
a
^
k
+
1
−
b
a
k
]
×
δ
t
2
f_{23}=-\frac{1}{2}(R_{k}+R_{k+1})\deltat,\f_{24}=\frac{1}{2}R_{k+1}[\hat{a}_{k+1}-b_{a{k}}]_{\times}\deltat^{2}
f23=−21(Rk+Rk+1)δt, f24=21Rk+1[a^k+1−bak]×δt2
v
00
=
−
1
4
R
k
δ
t
2
,
v
01
=
v
03
=
1
2
v
21
δ
t
,
v
02
=
−
1
4
R
k
+
1
δ
t
2
v_{00}=-\frac{1}{4}R_{k}\deltat^{2},\v_{01}=v_{03}=\frac{1}{2}v_{21}\deltat,\v_{02}=-\frac{1}{4}R_{k+1}\deltat^{2}
v00=−41Rkδt2, v01=v03=21v21δt, v02=−41Rk+1δt2
v
21
=
v
23
=
1
4
R
k
+
1
[
a
^
k
+
1
−
b
a
k
]
×
δ
t
2
v_{21}=v_{23}=\frac{1}{4}R_{k+1}[\hat{a}_{k+1}-b_{ak}]_{\times}\deltat^{2}
v21=v23=41Rk+1[a^k+1−bak]×δt2————————————其中上述符号中,为了与预积分(4.2.2)的推导对应,应有:
R
k
=
=
γ
^
i
b
k
.
t
o
R
o
t
a
t
i
o
n
M
a
t
r
i
x
(
)
R_{k}==\hat{\gamma}_{i}^{b_{k}}.toRotationMatrix()
Rk==γ^ibk.toRotationMatrix(),
a
^
k
=
=
a
i
^
\hat{a}_{k}==\hat{a_{i}}
a^k==ai^,
b
a
k
=
=
b
i
b_{ak}==b_{i}
bak==bi,
R
k
+
1
=
γ
^
i
+
1
b
k
.
t
o
R
o
t
a
t
i
o
n
M
a
t
r
i
x
(
)
R_{k+1}=\hat{\gamma}_{i+1}^{b_{k}}.toRotationMatrix()
Rk+1=γ^i+1bk.toRotationMatrix(),
a
^
k
+
1
=
=
a
^
i
+
1
\hat{a}_{k+1}==\hat{a}_{i+1}
a^k+1==a^i+1,
b
a
k
=
b
a
i
b_{ak}=b_{ai}
bak=bai,
w
^
k
=
=
w
^
i
\hat{w}_{k}==\hat{w}_{i}
w^k==w^i,
w
^
k
+
1
=
=
w
^
i
+
1
\hat{w}_{k+1}==\hat{w}_{i+1}
w^k+1==w^i+1,
b
w
k
=
=
b
w
i
b_{wk}==b_{wi}
bwk==bwi
.
t
o
R
o
t
a
t
i
o
n
M
a
t
r
i
x
(
)
.toRotationMatrix()
.toRotationMatrix()表示将单位四元素转成旋转矩阵————————————
简写为:
δ
z
k
+
1
15
×
1
=
F
15
×
15
δ
z
k
15
×
1
+
V
15
×
18
n
18
×
1
(10)
\deltaz_{k+1}^{15\times1}=F^{15\times15}\deltaz_{k}^{15\times1}+V^{15\times18}n^{18\times1}\tag{10}
δzk+115×1=F15×15δzk15×1+V15×18n18×1(10)协方差的迭代公式为(协方差的初始值
P
0
=
0
P_{0}=0
P0=0):
P
k
+
1
=
F
P
k
F
T
+
V
Q
V
T
P_{k+1}=FP_{k}F^{T}+VQV^{T}
Pk+1=FPkFT+VQVT其中噪声项的对角协方差矩阵
Q
Q
Q为:
Q
18
×
18
=
(
σ
a
2
,
σ
w
2
,
σ
a
2
,
σ
w
2
,
σ
b
a
2
,
σ
b
w
2
)
Q^{18\times18}=(\sigma_{a}^{2},\sigma_{w}^{2},\sigma_{a}^{2},\sigma_{w}^{2},\sigma_{ba}^{2},\sigma_{bw}^{2})
Q18×18=(σa2,σw2,σa2,σw2,σba2,σbw2)。
雅克比矩阵的迭代公式为(雅克比矩阵的初始值为
J
0
=
I
J_{0}=I
J0=I):
J
k
+
1
=
F
J
k
J_{k+1}=FJ_{k}
Jk+1=FJk这里计算出来的
J
k
+
1
J_{k+1}
Jk+1只是为了给后面提供对bias的雅克比矩阵。
关于公式(9)的推导:1)对于
δ
b
a
k
+
1
\deltab_{ak+1}
δbak+1、
δ
b
w
k
+
1
\deltab_{wk+1}
δbwk+1,根据偏置模型定义显然有:
δ
b
a
k
+
1
=
δ
b
a
k
+
n
b
a
δ
t
,
δ
b
w
k
+
1
=
δ
b
w
k
+
n
b
w
δ
t
\deltab_{ak+1}=\deltab_{ak}+n_{b_{a}}\deltat,\\deltab_{wk+1}=\deltab_{wk}+n_{b_{w}}\deltat
δbak+1=δbak+nbaδt, δbwk+1=δbwk+nbwδt2)对于
δ
θ
k
+
1
\delta\theta_{k+1}
δθk+1:之前已经得到了角度误差导数的连续形式:
δ
θ
˙
t
b
k
=
−
[
w
^
t
−
b
w
t
]
×
δ
θ
t
b
k
−
n
w
−
δ
b
w
t
\delta\dot{\theta}_{t}^{b_{k}}=-[\hat{w}_{t}-b_{wt}]_{\times}\delta\theta_{t}^{b_{k}}-n_{w}-\deltab_{wt}
δθ˙tbk=−[w^t−bwt]×δθtbk−nw−δbwt那么其中值积分的离散形式为:
δ
θ
˙
k
=
−
[
w
^
k
+
w
^
k
+
1
2
−
b
w
k
]
×
δ
θ
k
b
k
−
n
w
k
+
n
w
k
+
1
2
−
δ
b
w
k
\delta\dot{\theta}_{k}=-[\frac{\hat{w}_{k}+\hat{w}_{k+1}}{2}-b_{wk}]_{\times}\delta\theta_{k}^{b_{k}}-\frac{n_{wk}+n_{wk+1}}{2}-\deltab_{wk}
δθ˙k=−[2w^k+w^k+1−bwk]×δθkbk−2nwk+nwk+1−δbwk根据导数定义可得下一时刻的误差为:
δ
θ
k
+
1
=
δ
θ
k
+
δ
θ
˙
k
δ
t
=
(
I
−
[
w
^
k
+
w
^
k
+
1
2
−
b
w
k
]
×
δ
t
)
δ
θ
k
−
n
w
k
+
n
w
k
+
1
2
δ
t
−
δ
b
w
k
δ
t
\delta\theta_{k+1}=\delta\theta_{k}+\delta\dot{\theta}_{k}\deltat=(I-[\frac{\hat{w}_{k}+\hat{w}_{k+1}}{2}-b_{wk}]_{\times}\deltat)\delta\theta_{k}-\frac{n_{wk}+n_{wk+1}}{2}\deltat-\deltab_{wk}\deltat
δθk+1=δθk+δθ˙kδt=(I−[2w^k+w^k+1−bwk]×δt)δθk−2nwk+nwk+1δt−δbwkδt整理得
δ
θ
k
+
1
=
(
I
−
[
w
^
k
+
w
^
k
+
1
2
−
b
w
k
]
×
δ
t
)
δ
θ
k
−
δ
t
δ
b
w
k
−
δ
t
2
n
w
k
−
δ
t
2
n
w
k
+
1
=
f
11
δ
θ
k
−
δ
t
δ
b
w
k
−
δ
t
2
n
w
k
−
δ
t
2
n
w
k
+
1
\delta\theta_{k+1}=(I-[\frac{\hat{w}_{k}+\hat{w}_{k+1}}{2}-b_{wk}]_{\times}\deltat)\delta\theta_{k}-\deltat\deltab_{wk}-\frac{\deltat}{2}n_{wk}-\frac{\deltat}{2}n_{wk+1}=f_{11}\delta\theta_{k}-\deltat\deltab_{wk}-\frac{\deltat}{2}n_{wk}-\frac{\deltat}{2}n_{wk+1}
δθk+1=(I−[2w^k+w^k+1−bwk]×δt)δθk−δtδbwk−2δtnwk−2δtnwk+1=f11δθk−δtδbwk−2δtnwk−2δtnwk+1故有:
f
11
=
I
−
[
w
^
k
+
w
^
k
+
1
2
−
b
w
k
]
×
δ
t
f_{11}=I-[\frac{\hat{w}_{k}+\hat{w}_{k+1}}{2}-b_{wk}]_{\times}\deltat
f11=I−[2w^k+w^k+1−bwk]×δt3)对于
δ
β
k
+
1
:
\delta\beta_{k+1}:
δβk+1:之前已经得到了速度误差导数的连续形式:
δ
β
˙
t
b
k
=
−
R
t
b
k
[
a
^
t
−
b
a
t
]
×
δ
θ
−
R
t
b
k
δ
b
a
t
−
R
t
b
k
n
a
\delta\dot{\beta}_{t}^{b_{k}}=-R_{t}^{b_{k}}[\hata_{t}-b_{at}]_{\times}\delta\theta-R_{t}^{b_{k}}\deltab_{at}-R_{t}^{b_{k}}n_{a}
δβ˙tbk=−Rtbk[a^t−bat]×δθ−Rtbkδbat−Rtbkna那么其中值积分的离散形式为:
δ
β
˙
k
=
−
1
2
R
k
[
a
^
k
−
b
a
k
]
×
δ
θ
k
−
1
2
R
k
+
1
[
a
^
k
+
1
−
b
a
k
+
1
]
×
δ
θ
k
+
1
−
1
2
(
R
k
+
R
k
+
1
)
δ
b
a
k
−
1
2
R
k
n
a
k
−
1
2
R
k
+
1
n
a
k
+
1
\delta\dot{\beta}_{k}=-\frac{1}{2}R_{k}[\hata_{k}-b_{ak}]_{\times}\delta\theta_{k}-\frac{1}{2}R_{k+1}[\hata_{k+1}-b_{ak+1}]_{\times}\delta\theta_{k+1}-\frac{1}{2}(R_{k}+R_{k+1})\deltab_{ak}-\frac{1}{2}R_{k}n_{ak}-\frac{1}{2}R_{k+1}n_{ak+1}
δβ˙k=−21Rk[a^k−bak]×δθk−21Rk+1[a^k+1−bak+1]×δθk+1−21(Rk+Rk+1)δbak−21Rknak−21Rk+1nak+1其中假设在
k
k
k时刻和
k
+
1
k+1
k+1时刻下的加速度计偏置
b
a
b_{a}
ba相等(下同),把之前的
δ
θ
k
+
1
\delta\theta_{k+1}
δθk+1代进去有:
δ
β
˙
k
=
{
−
1
2
R
k
[
a
^
k
−
b
a
k
]
×
−
1
2
R
k
+
1
[
a
^
k
+
1
−
b
a
k
]
×
(
I
−
[
w
^
k
+
w
^
k
+
1
2
−
b
w
k
]
×
δ
t
)
}
δ
θ
k
+
1
4
R
k
+
1
[
a
^
k
+
1
−
b
a
k
]
×
δ
t
(
n
w
k
+
n
w
k
+
1
)
+
1
2
R
k
+
1
[
a
^
k
+
1
−
b
a
k
]
×
δ
t
δ
b
w
k
−
1
2
(
R
k
+
R
k
+
1
)
δ
b
a
k
−
1
2
R
k
n
a
k
−
1
2
R
k
+
1
n
a
k
+
1
\delta\dot{\beta}_{k}=\begin{Bmatrix}-\frac{1}{2}R_{k}[\hata_{k}-b_{ak}]_{\times}-\frac{1}{2}R_{k+1}[\hata_{k+1}-b_{ak}]_{\times}(I-[\frac{\hat{w}_{k}+\hat{w}_{k+1}}{2}-b_{wk}]_{\times}\deltat)\end{Bmatrix}\delta\theta_{k}\\+\frac{1}{4}R_{k+1}[\hata_{k+1}-b_{ak}]_{\times}\deltat(n_{wk}+n_{wk+1})+\frac{1}{2}R_{k+1}[\hata_{k+1}-b_{ak}]_{\times}\deltat\deltab_{wk}\\-\frac{1}{2}(R_{k}+R_{k+1})\deltab_{ak}-\frac{1}{2}R_{k}{n_{ak}}-\frac{1}{2}R_{k+1}{n_{ak+1}}
δβ˙k={−21Rk[a^k−bak]×−21Rk+1[a^k+1−bak]×(I−[2w^k+w^k+1−bwk]×δt)}δθk+41Rk+1[a^k+1−bak]×δt(nwk+nwk+1)+21Rk+1[a^k+1−bak]×δtδbwk−21(Rk+Rk+1)δbak−21Rknak−21Rk+1nak+1根据导数定义可得到下一时刻的误差为:
δ
β
k
+
1
=
δ
β
k
+
δ
β
˙
k
δ
t
=
{
−
1
2
R
k
[
a
^
k
−
b
a
k
]
×
−
1
2
R
k
+
1
[
a
^
k
+
1
−
b
a
k
]
×
(
I
−
[
w
^
k
+
w
^
k
+
1
2
−
b
w
k
]
×
δ
t
)
}
δ
t
δ
θ
k
+
δ
β
k
+
{
1
4
R
k
+
1
[
a
^
k
+
1
−
b
a
k
]
×
δ
t
(
n
w
k
+
n
w
k
+
1
)
+
1
2
R
k
+
1
[
a
^
k
+
1
−
b
a
k
]
×
δ
t
δ
b
w
k
−
1
2
(
R
k
+
R
k
+
1
)
δ
b
a
k
−
1
2
R
k
n
a
k
−
1
2
R
k
+
1
n
a
k
+
1
}
δ
t
\delta\beta_{k+1}=\delta\beta_{k}+\delta\dot{\beta}_{k}\deltat\\=\begin{Bmatrix}-\frac{1}{2}R_{k}[\hata_{k}-b_{ak}]_{\times}-\frac{1}{2}R_{k+1}[\hata_{k+1}-b_{ak}]_{\times}(I-[\frac{\hat{w}_{k}+\hat{w}_{k+1}}{2}-b_{wk}]_{\times}\deltat)\end{Bmatrix}\deltat\delta\theta_{k}+\delta\beta_{k}+\\\begin{Bmatrix}\frac{1}{4}R_{k+1}[\hata_{k+1}-b_{ak}]_{\times}\deltat(n_{wk}+n_{wk+1})+\frac{1}{2}R_{k+1}[\hata_{k+1}-b_{ak}]_{\times}\deltat\deltab_{wk}-\frac{1}{2}(R_{k}+R_{k+1})\deltab_{ak}-\frac{1}{2}R_{k}{n_{ak}}-\frac{1}{2}R_{k+1}{n_{ak+1}}\end{Bmatrix}\deltat
δβk+1=δβk+δβ˙kδt={−21Rk[a^k−bak]×−21Rk+1[a^k+1−bak]×(I−[2w^k+w^k+1−bwk]×δt)}δtδθk+δβk+{41Rk+1[a^k+1−bak]×δt(nwk+nwk+1)+21Rk+1[a^k+1−bak]×δtδbwk−21(Rk+Rk+1)δbak−21Rknak−21Rk+1nak+1}δt整理后得到:
δ
β
k
+
1
=
f
21
δ
θ
k
+
δ
β
k
+
f
23
δ
b
a
k
+
f
24
δ
b
w
k
−
1
2
R
k
δ
t
n
a
k
−
1
2
R
k
+
1
δ
t
n
a
k
+
1
+
v
21
n
w
k
+
v
23
n
w
k
+
1
\delta\beta_{k+1}=f_{21}\delta\theta_{k}+\delta\beta_{k}+f_{23}\deltab_{ak}+f_{24}\deltab_{wk}-\frac{1}{2}R_{k}\deltatn_{ak}\\-\frac{1}{2}R_{k+1}\deltatn_{ak+1}+v_{21}n_{wk}+v_{23}n_{wk+1}
δβk+1=f21δθk+δβk+f23δbak+f24δbwk−21Rkδtnak−21Rk+1δtnak+1+v21nwk+v23nwk+1其中
f
21
=
−
1
2
R
k
[
a
^
k
−
b
a
k
]
×
δ
t
−
1
2
R
k
+
1
[
a
^
k
+
1
−
b
a
k
]
×
(
I
−
[
w
^
k
+
w
^
k
+
1
2
−
b
w
k
]
×
δ
t
)
δ
t
f_{21}=-\frac{1}{2}R_{k}[\hata_{k}-b_{ak}]_{\times}\deltat-\frac{1}{2}R_{k+1}[\hata_{k+1}-b_{ak}]_{\times}(I-[\frac{\hat{w}_{k}+\hat{w}_{k+1}}{2}-b_{wk}]_{\times}\deltat)\deltat
f21=−21Rk[a^k−bak]×δt−21Rk+1[a^k+1−bak]×(I−[2w^k+w^k+1−bwk]×δt)δt
f
23
=
−
1
2
(
R
k
+
R
k
+
1
)
δ
t
f_{23}=-\frac{1}{2}(R_{k}+R_{k+1})\deltat
f23=−21(Rk+Rk+1)δt
f
24
=
1
2
R
k
+
1
[
a
^
k
+
1
−
b
a
k
]
×
δ
t
2
f_{24}=\frac{1}{2}R_{k+1}[\hat{a}_{k+1}-b_{a_{k}}]_{\times}\deltat^{2}
f24=21Rk+1[a^k+1−bak]×δt2
v
21
=
v
23
=
1
4
R
k
+
1
[
a
^
k
+
1
−
b
a
k
]
×
δ
t
2
v_{21}=v_{23}=\frac{1}{4}R_{k+1}[\hat{a}_{k+1}-b_{ak}]_{\times}\deltat^{2}
v21=v23=41Rk+1[a^k+1−bak]×δt2证毕!4)对于
δ
α
k
+
1
:
\delta\alpha_{k+1}:
δαk+1:之前已经得到了速度误差导数的连续形式:
δ
α
˙
t
b
k
=
α
˙
^
t
b
k
−
α
˙
t
b
k
=
β
^
t
b
k
−
β
t
b
k
=
δ
β
t
b
k
\delta\dot{\alpha}_{t}^{b_{k}}=\hat{\dot{\alpha}}_{t}^{b_{k}}-\dot{\alpha}_{t}^{b_{k}}=\hat{\beta}_{t}^{b_{k}}-\beta_{t}^{b_{k}}=\delta\beta_{t}^{b_{k}}
δα˙tbk=α˙^tbk−α˙tbk=β^tbk−βtbk=δβtbk则离散形式为:
δ
α
˙
t
b
k
=
1
2
δ
β
k
+
1
2
δ
β
k
+
1
=
δ
β
k
+
1
2
f
21
δ
θ
k
−
1
4
(
R
k
+
R
k
+
1
)
δ
t
δ
b
a
k
+
1
2
f
24
δ
b
w
k
−
1
4
R
k
δ
t
n
a
k
−
1
4
R
k
+
1
δ
t
n
a
k
+
1
+
1
2
v
21
n
w
k
+
1
2
v
23
n
w
k
+
1
\delta\dot{\alpha}_{t}^{b_{k}}=\frac{1}{2}\delta\beta_{k}+\frac{1}{2}\delta\beta_{k+1}\\=\delta\beta_{k}+\frac{1}{2}f_{21}\delta\theta_{k}-\frac{1}{4}(R_{k}+R_{k+1})\deltat\deltab_{ak}+\frac{1}{2}f_{24}\deltab_{wk}\\-\frac{1}{4}R_{k}\deltatn_{ak}-\frac{1}{4}R_{k+1}\deltatn_{ak+1}+\frac{1}{2}v_{21}n_{wk}+\frac{1}{2}v_{23}n_{wk+1}
δα˙tbk=21δβk+21δβk+1=δβk+21f21δθk−41(Rk+Rk+1)δtδbak+21f24δbwk−41Rkδtnak−41Rk+1δtnak+1+21v21nwk+21v23nwk+1根据导数定义可得下一时刻的误差为:
δ
α
k
+
1
=
δ
α
k
+
δ
α
˙
δ
t
=
δ
α
k
+
f
01
δ
θ
k
+
δ
t
δ
β
k
+
f
03
δ
b
a
k
+
f
04
δ
b
w
k
+
v
00
n
a
k
+
v
01
n
w
k
+
v
02
n
a
k
+
1
+
v
03
n
w
k
+
1
\delta\alpha_{k+1}=\delta{\alpha}_{k}+\delta\dot{\alpha}\deltat\\=\delta\alpha_{k}+f_{01}\delta\theta_{k}+\deltat\delta\beta_{k}+f_{03}\deltab_{ak}+f_{04}\deltab_{wk}+v_{00}n_{ak}+v_{01}n_{wk}+v_{02}n_{ak+1}+v_{03}n_{wk+1}
δαk+1=δαk+δα˙δt=δαk+f01δθk+δtδβk+f03δbak+f04δbwk+v00nak+v01nwk+v02nak+1+v03nwk+1其中
f
01
=
1
2
f
21
δ
t
,
f
03
=
−
1
4
(
R
k
+
R
k
+
1
)
δ
t
2
,
f
04
=
1
2
f
24
δ
t
f_{01}=\frac{1}{2}f_{21}\deltat,\f_{03}=-\frac{1}{4}(R_{k}+R_{k+1})\deltat^{2},\f_{04}=\frac{1}{2}f_{24}\delta{t}
f01=21f21δt, f03=−41(Rk+Rk+1)δt2, f04=21f24δt
v
00
=
−
1
4
R
k
δ
t
2
,
v
01
=
v
03
=
1
2
v
21
δ
t
,
v
02
=
−
1
4
R
k
+
1
δ
t
2
v_{00}=-\frac{1}{4}R_{k}\deltat^{2},\v_{01}=v_{03}=\frac{1}{2}v_{21}\deltat,\v_{02}=-\frac{1}{4}R_{k+1}\deltat^{2}
v00=−41Rkδt2, v01=v03=21v21δt, v02=−41Rk+1δt2
收起
展开全文
IMU预积分
神经网络学习4【误差传递与权重更新】
2021-08-1516:23:52
1.4.误差反馈1.4.1误差反馈校正权重矩阵可以理解,输出和误差都是多个节点共同作用的结果,那么该如何更新链接权重?思考一下,得到误差后,该怎么分配?平均分的话是否会有失公平?毕竟我们在之前的学习中...
1.误差反馈
1.1误差反馈校正权重矩阵
可以理解,输出和误差都是多个节点共同作用的结果,那么该如何更新链接权重?思考一下,得到误差后,该怎么分配?平均分的话是否会有失公平?毕竟我们在之前的学习中了解到,前一层每个节点的贡献都是不一样的。
考虑极端情况,当某权重为0时,它对下一个节点的贡献为0;这时如果误差仍然平均分配显然是不那么合适的。
但我们很容易想到这个标准:较大链接权重的连接分配更多的误差。
将同样的思想扩展到多个节点。
如果我们拥有100个节点链接到输出节点,那么我们要在这100条链接之间,按照每条链接对误差所做贡献的比例(由链接权重的大小表示),分割误差。
我们使用权重,将误差从输出向后传播到网络中。
我们称这种方法为反向传播。
1.2多个输出节点反向传播误差
第一个输出节点的误差标记为e1=(t1-o1)。
按照所连接链接的比例,也就是权重w1,1和w2,1,对误差e1进行分割。
同理,按照权重w1,2和w2,2的比例分割e2。
1.3反向传播误差到更多层中
将输出误差标记为eoutput,将在输出层和隐藏层之间的链接权重标记为who。
通过将误差值按权重的比例进行分割,我们计算出与每条链接相关的特定误差值。
采用与隐藏层节点相关联的这些误差ehidden,再次将这些误差按照输入层和隐藏层之间的链接权重wih进行分割。
如果神经网络具有多个层,那么我们就从最终输出层往回工作,对每一层重复应用相同的思路。
具有实际数字的3层网络中,误差如何向后传播:
进一步向后工作:
1.神经网络通过调整链接权重进行学习。
这种方法由误差引导,误差就是训练数据所给出正确答案和实际输出之间的差值。
2.简单地说,在输出节点处的误差等于所需值与实际值之间的差值。
3.与内部节点相关联的误差并不显而易见。
一种方法是按照链路权重的比例来分割输出层的误差,然后在每个内部节点处重组这些误差。
1.4使用矩阵乘法进行反向传播误差
计算的起始点是在神经网络最终输出层中出现的误差。
为隐藏层的误差构建矩阵。
可以观察到,最重要的事情是输出误差与链接权重wij的乘法。
这些分数的分母是一种归一化因子。
如果我们忽略了这个因子,那么我们仅仅失去后馈误差的大小。
也就是说,可以使用简单得多的e1*w1,1来代替e1*w1,1/(w1,1+w2,1)。
因此,我们得到所希望的矩阵,使用矩阵的方法来向后传播误差:
反向传播误差可以表示为矩阵乘法。
无论网络规模大小,这使我们能够简洁地表达反向传播误差。
前向馈送信号和反向传播误差都可以使用矩阵计算而变得高效。
2.更新权重
2.1梯度下降法
目前我们已经完成了让误差反向传播到网络的每一层,我们的目的是要使用误差来指导如何调整链接权重,从而使得神经网络更适合这个样本集,输出更准确的答案。
如果我们将复杂困难的函数当作网络误差,那么找到最小值就意味着最小化误差。
这样我们就可以改进网络输出。
如此我们就可以理解要用梯度下降法来更新权重的意义所在了。
可以从梯度的物理意义上直观理解,就像是在三维地形中,在曲面的任意一点放置一个静止的小球,在重力作用下小球会向下滚动,那么初始时刻,小球滚动的方向一定是该点周围下降最明显的地方,即负梯度方向。
但最终小球会在摩擦力和重力作用下停下,可以想象,这必是一个“凹谷”,也就是这个三维曲面的(局部)极小值。
但如果这个地形非常复杂,我们掉入的仅仅是这个点周围某个“凹谷”,而非整个地形中最低(即我们所要求的最小值)呢?该如何解决这个问题?首先,采用梯度下降法,我们必须选择一个起点。
这一次,我们所在之处的斜率为正,因此我们向左移动。
为了避免超调,避免在最小值的地方来回反弹,就要改变步子大小,这是一个必要的优化。
回想刚刚的问题,如果恰好掉入了错误的山谷怎么办?记得我们刚刚讲的随机选择起点吗?我们可以随机取,并且多次取,尽可能地,把错误的结果过滤掉。
多次训练神经网络,确保并不总是终止于错误的山谷。
不同的起始点意味着选择不同的起始参数,在神经网络的情况下,这意味着选择不同的起始链接权重。
像这样的尝试,终究会带领我们走进正确的山谷。
2.2误差函数
那么,梯度下降法所作用的误差函数应该如何构造呢?像我们之前那样,目标值-真实值=误差?但这种情况下的误差会有正负,一旦我们求和,就有可能得到非常趋近于0的结果,这会让我们误以为误差很小。
误差函数如何构造?看下表:
第三种选择是差的平方,即(目标值-实际值)^2。
我们更喜欢使用第三种误差函数,而不喜欢使用第二种误差函数,原因有以下几点:
使用误差的平方,我们可以很容易使用代数计算出梯度下降的斜率。
误差函数平滑连续,这使得梯度下降法很好地发挥作用——没有间断,也没有突然的跳跃。
越接近最小值,梯度越小,这意味着,如果我们使用这个函数调节步长,超调的风险就会变得较小。
下图显示了两个链接权重,这次,误差函数是三维曲面,这个曲面随着两个链接权重的变化而变化。
这个表达式表示了当权重wj,k改变时,误差E是如何改变的。
这是误差函数的斜率,也就是我们希望使用梯度下降的方法到达最小值的方向。
下为推导过程:
第一部分就是(目标值-实际值)。
在sigmoid中的求和表达式也很简单,就是进入最后一层节点的信号,我们可以称之为ik,这样它看起来比较简单。
这是应用激活函数之前,进入节点的信号。
最后一部分是前一隐藏层节点j的输出。
同理,得到输入层和隐藏层之间权重调整:权重改变的方向与梯度方向相反,我们用数学的形式来表达这个因子。
(和之前在线性分类器里学到的是不是原理相同?newA=oldA+ΔA)
更新后的权重wj,k是由刚刚得到误差斜率取反来调整旧的权重而得到的。
正如我们先前所看到的,如果斜率为正,我们希望减小权重,如果斜率为负,我们希望增加权重,因此,我们要对斜率取反。
符号α是一个因子,这个因子可以调节这些变化的强度,确保不会超调。
我们通常称这个因子为学习率。
同理,这个表达式不仅适用于隐藏层和输出层之间的权重,而且适用于输入层和隐藏层之间的权重。
由于学习率只是一个常数,并没有真正改变如何组织矩阵乘法,因此我们省略了学习率α。
1.神经网络的误差是内部链接权重的函数。
2.改进神经网络,意味着通过改变权重减少这种误差。
3.直接选择合适的权重太难了。
另一种方法是,通过误差函数的梯度下降,采取小步长,迭代地改进权重。
所迈出的每一步的方向都是在当前位置向下斜率最大的方向,这就是所谓的梯度下降。
4.使用微积分可以很容易地计算出误差斜率。
2.3误差更新成功范例
更新隐藏层和输出层之间的权重w1,1。
当前,这个值为2.0。
第一项(tk-ok)得到误差e1=0.8。
S函数内的求和Σjwj,koj为(2.0×0.4)+(3.0*0.5)=2.3。
sigmoid1/(1+e-2.3)为0.909。
中间的表达式为0.909*(1-0.909)=0.083。
由于我们感兴趣的是权重w1,1,其中j=1,因此最后一项oj也很简单,也就是oj=1。
此处,oj值就是0.4。
将这三项相乘,最后我们得到-0.0265。
如果学习率为0.1,那么得出的改变量为-(0.1*-0.02650)=+0.002650。
因此,新的w1,1就是原来的2.0加上0.00265等于2.00265。
收起
展开全文
神经网络
python
GRU网络
千次阅读
2019-09-2910:24:52
采用后向误差传播算法来学习网络,所以先得求损失函数对各参数的偏导(总共有7个): 其中各中间参数为: 在算出了对各参数的偏导...
收起
gru
推导和实现:全面解析高斯过程中的函数最优化(附代码&公式)
千次阅读
2018-04-0100:00:00
在上式的右边,分子中的第一项需要我们对测量过程中的误差来源做一些假设,分子中的第二项是先验概率,在这里我们必须采用最合理的假设。
例如,我们将在下面看到先验概率有效地决定了f函数在给定平滑度的概率。
在...
收起
神经网络基本原理、误差逆传播BP算法公式推导与多层神经网络的Python实现
2021-12-0501:36:13
学习西瓜书神经网络章节的笔记,对神经网络的相关公式进行了推导,并用python进行实现,欢迎交流指正~
收起
神经网络
python
机器学习
人工智能
MakeYourOwnNeuralNetwork(十一)-----多层神经网络层反向传输误差
千次阅读
2018-05-0320:20:26
利用矩阵乘法计算反向传输误差实际上是如何更新权重(一)实际上是如何更新权重(二)权重更新实例多层神经网络层反向传输误差下图所示的是包含一个输出层、一个隐藏层和一个输出层的简单的神经网络。
从上图最...
收起
神经网络
Java浮点数float和double精确计算的精度误差问题总结
千次阅读
2017-07-0515:21:45
构造函数不恰当,在传递给JDBC setBigDecimal() 方法时,会造成似乎很奇怪的JDBC驱动程序中的异常。
例如,考虑以下JDBC代码,该代码希望将数字 0.01 存储到小数字段:PreparedStatement...
收起
MLP公式推导及numpy实现
2021-06-0217:10:36
文章目录全连接神经网络公式推导及numpy实现1.预备知识1.1链导法则1.2Sigmoid激活函数1.3Softmax激活函数1.4交叉熵损失2.前向传播3.梯度反向传播4.训练流程5.代码实现全连接神经网络公式推导及numpy实现...
收起
mlp
神经网络
机器学习
计算机组成原理复习笔记-2
千次阅读
多人点赞
2018-05-2710:37:04
创建海明编码的方法:首先根据公式确定编码所需的校验位数目r,算出编码字的位长度n=m+r,从右向左从1开始编号。
位数是2的指数幂的位设置为奇偶校验位,其他位为数据位。
对于各个编码位置,第b位编码由满足b=...
收起
计算机组成原理
10-误差反向传播法(一)——计算图
2020-07-2523:58:42
误差反向传播法上一篇文章介绍的是用数值微分法计算梯度,但这种方法比较耗时间,接下来介绍新的梯度计算法:误差反向传播法。
在此之前,先介绍计算图。
一、计算图计算图用节点和箭头表示,节点表示某种运算...
收起
深度学习
深度学习-误差反向传播直观理解
2019-11-2621:17:31
刚开始接触机器学习或者深度学习的时候,学到反向传播的时候总是感觉很抽象,看了知乎和一些教材后还是感觉不够直观,后来看到一处介绍后豁然开朗,在此做下笔记方便查阅,全文不涉及公式推导,仅作为直观理解,待...
收起
反向传播
深度学习
误差反向传播算法浅解
万次阅读
2017-09-2809:24:48
反向传播(英语:Backpropagation,缩写为BP)是“误差反向传播”的简称。
由于多层前馈神经网络的训练经常采用误差反向传播算法,人们也常把多层前馈神经网络称为BP网络。
收起
BP算法
反向传播
反向传播算法(过程及公式推导)
万次阅读
多人点赞
2016-04-0121:19:56
表示Hadamard乘积,用于矩阵或向量之间点对点的乘法运算。
公式1的推导过程如下:公式2(由后往前,计算每一层神经网络产生的错误): 推导过程:...
收起
反向传播
算法步骤
并行处理mpi矩阵乘法
千次阅读
2014-01-1121:53:55
基于MPI并行方法实现矩阵乘法目录1. 实验目的32. 实验环境43. 实验内容43.1. 实验题目43.2. 实验过程53.2.1. 集群使用53.2.2. 源码及解析53.3. 执行时间截图123.3.1. 基准程序...
收起
python深度学习入门-误差反向传播法
2021-05-0118:08:10
深度学习入门-误差反向传播法摘要通过使用计算图,可以直观地把握计算过程。
计算图的节点是由局部计算构成的。
局部计算构成全局计算。
计算图的正向传播进行一般的计算。
通过计算图的反向传播,可以计算...
收起
第五章误差反向传播算法——基于numpy的代码详解
千次阅读
2020-04-1512:43:44
softmax层的反向传播比较复杂,我们将分6个步骤进行:步骤一:cross entropy error层的反向传播的值传递过来:步骤二:“×”节点将正向传播的值翻转后相乘,需要用到的公式是:步骤三:正向传播时若有分支流出...
收起
神经网络
python
深度学习
numpy
到底什么是最小二乘法
2018-08-0210:00:10
最小二乘法,又是一个即熟悉又陌生的名字。
对于学工科的我,简直就是听着最小二乘长大的(汗。
。
。
...),并没有系统的了解一下什么是最小二乘法。
包括最小二乘这个叫法,也从来都不理解(一直以...
收起
最小二乘法
误差反向传播和深度学习相关技巧总结
2019-10-2311:20:36
误差反向传播和深度学习相关技巧总结一、误差反向传播法1.几个问题误差反向传播的目的是什么?为了能够更高效的计算权重...2.简单层(加法、乘法层)、激活函数层、Affine/softmax层的实现说明:通过计算图,...
收起
神经网络学习
Matlab嵌套传递函数简化_MATLAB的数据处理方法及图形绘制详解
2020-11-2013:26:38
矩阵的左除和右除:注意只有当两个矩阵中前一个矩阵的列数和后一个矩阵的行数相同时,才可以进行乘法运算。
a\b运算等效于求a*x=b的解;而a/b等效于求x*b=a的解。
只有方阵才可以求幂。
2.逆矩阵与行列式计算求逆:inv...
收起
r语言中正定矩阵由于误差不正定_Kalman滤波在MOT中的应用(一)——理论篇
2020-10-2618:14:52
可利用他们的均方误差来计算:由此可将均方误差不断地传递下去,从而估计估算出最优的温度值。
为了方便解释卡尔曼滤波方程,下面以一辆小车的运动为例,假设我们已知上一时刻小车的状态,现在要估计当前时刻的状态...
收起
深度学习入门笔记(六):误差反向传播算法
千次阅读
多人点赞
2020-07-1521:37:04
误差反向传播法——能够高效计算权重参数的梯度的方法。
要正确理解误差反向传播法,有两种方法:一种是基于数学式;另一种是基于计算图(computationalgraph)。
前者是比较常见的方法,机器学习相关的图书中多数都是...
收起
神经网络
深度学习
机器学习
人工智能
深度学习入门学习笔记之——误差反向传播法
2021-01-1321:57:57
误差反向传播法上一章中,我们介绍了神经网络的学习,并通过数值微分计算了神经网络的权重参数的梯度(严格来说,是损失函数关于权重参数的梯度)。
数值微分虽然简单,也容易实现,但缺点是计算上比较费时间。
本章...
收起
神经网络
深度学习
深度学习入门之4--误差反向传播法
2020-01-2916:49:15
目录1计算图1.1计算图求解...3.2乘法节点反向传播4简单层的实现4.1乘法层的实现4.2加法层的实现5激活函数模块化5.1Relu层5.2Sigmoid层5.3Affine层5.3.1Affine层5.3.2批版本...
收起
深度学习之神经网络传递流程
2020-05-1412:30:02
核心思路:传递的过程其实就是矩阵的乘法以及计算加权平均后经过激活函数就完成了一次传递,以此类推进行传递1.数据输入网络输入层:输入层中输入的是一些矩阵,比如可以输入(1,2)(1,2)(1,2)这个点坐标作为一个...
收起
深度学习
神经网络
空空如也
空空如也
1
2
3
4
5
...
20
收藏数
2,908
精华内容
1,163
热门标签
九九乘法表
矩阵乘法
99乘法表
打印九九乘法表
python九九乘法表
九九乘法表c语言编程
大整数乘法
python矩阵乘法
java九九乘法表
关键字:乘法的误差传递公式
延伸文章資訊
- 1實驗數據的處理與分析 - 台灣師範大學物理系
誤差傳遞:. 經常一個物理量是經由測量數個物理量,再藉由關係式計算而得出。 例如:動量是由測量值質量 ...
- 2誤差傳遞乘法完整相關資訊
提示:乘法的誤差傳遞公式:(. . ̅. ).[PDF] 誤差分析簡介2012年9月27日· 如無法準確估計誤差值時,該怎麼辦?例如: 單擺的擺動... Error propagation. 為...
- 3乘法的误差传递公式 - CSDN
加法中的误差传递: X=u±v 则X的均方差为: σX=sqrt(σu^2+σv^2)...乘法中的误差传递: 除法中的误差传递: 有限次幂的误差的传播: 可以使用蒙特卡罗法来验证其误差: 如下...
- 4實驗一數據處理與分析
乘法. →有效數字為9.10. • 除法. → 有效數字為18 ... 實驗誤差. • 系統誤差. ─ 設備系統誤差. ─ 環境系統誤差. ─ 人為誤差 ... 誤差傳遞. • 加減的誤差傳遞.
- 5誤差傳遞加減乘除在PTT/Dcard完整相關資訊 - 小文青生活
... 當n 個量相加減時,其誤差傳遞的一般性公式為:. ∑. = = + ... (2) 乘除的誤差傳遞: yx xy =.[PDF] 基于误差传递模型的精密装配几何误差灵敏度分析- 北京理工...