本文将介绍一篇关于 四元数运动学的误差卡尔曼滤波
经典论文。本文结构如下:
- 第1章四元数定义和性质介绍,包括:
加法、减法、乘法(矩阵表示)、模、幂数、指数运算
等。 - 第2章旋转群定义和性质介绍,包括:
旋转矩阵表示旋转运算、四元数表示旋转运算、旋转矩阵和四元数转换
。 - 第3章李群定义和性质介绍,包括:
李群加减法运算、李群四种导数定义、常用雅可比矩阵、四元数导数
。 - 第4章IMU运动学方程介绍,包括:
局部坐标系统微分方程、全局坐标系统微分方程、全局坐标系统离散时刻运动学方程
。 - 第5章IMU/GPS融合定位介绍,包括:
融合定位预测、融合定位更新
。
论文英文版链接:https://hal.archives-ouvertes.fr/hal-01122406/document (October 12, 2017)
论文中文版链接:https://blog.csdn.net/sunqin_csdn/article/details/108560582
1. 四元数定义和性质
1.1 四元数定义
四元数可以这样进行定义,假设有两个复数
A
=
a
+
b
i
A=a+b i
A=a+bi 和
C
=
c
+
d
i
C=c+di
C=c+di,则
Q
Q
Q 可以写成
Q
=
A
+
C
j
Q=A+Cj
Q=A+Cj,并且规定
k
≜
i
j
k \triangleq i j
k≜ij,那么则得到四元数
Q
Q
Q 为:
Q
=
a
+
b
i
+
c
j
+
d
k
∈
H
Q=a+bi+cj+dk \in{\mathbb{H}}
Q=a+bi+cj+dk∈H
其中
{
a
,
b
,
c
,
d
}
∈
R
\{a, b, c, d\} \in \mathbb{R}
{a,b,c,d}∈R,
{
i
,
j
,
k
}
\{i,j,k\}
{i,j,k} 为三个虚数单位,三者并且满足以下性质:
i
2
=
j
2
=
k
2
=
i
j
k
=
−
1
i^2=j^2=k^2=ijk=-1
i2=j2=k2=ijk=−1
进一步,从中可以得出:
i
j
=
−
j
i
=
k
,
j
k
=
−
k
j
=
i
,
k
i
=
−
i
k
=
j
ij=-ji=k, \quad jk=-kj=i, \quad ki=-ik=j
ij=−ji=k,jk=−kj=i,ki=−ik=j
四元数由实部
和虚部
组成,如果实部为零,则得到纯四元数,记为
H
p
\mathbb{H}_p
Hp:
Q
=
b
i
+
c
j
+
d
k
∈
H
p
Q=bi+cj+dk \in \mathbb{H}_p
Q=bi+cj+dk∈Hp
为了计算方便,常把四元数写成标量
和向量
组合的形式,例如
Q
=
q
w
+
q
v
Q=q_w+\mathbf{q}_v
Q=qw+qv,其中
q
w
q_w
qw 是实部,
q
v
=
q
x
i
+
q
y
j
+
q
z
k
\mathbf{q}_v=q_xi+q_yj+q_zk
qv=qxi+qyj+qzk 为虚部,四元数也可以写成4维向量形式,即:
q
≜
[
q
w
q
v
]
=
[
q
w
q
x
q
y
q
z
]
\mathbf{q} \triangleq\left[\begin{array}{l} q_{w} \\ \mathbf{q}_{v} \end{array}\right]=\left[\begin{array}{l} q_{w} \\ q_{x} \\ q_{y} \\ q_{z} \end{array}\right]
q≜[qwqv]=⎣⎢⎢⎡qwqxqyqz⎦⎥⎥⎤
1.2 四元数主要性质
(1)加法:
p
±
q
=
[
p
w
p
v
]
±
[
q
w
q
v
]
=
[
p
w
±
q
w
p
v
±
q
v
]
\mathbf{p} \pm \mathbf{q}=\left[\begin{array}{l} p_{w} \\ \mathbf{p}_{v} \end{array}\right] \pm\left[\begin{array}{l} q_{w} \\ \mathbf{q}_{v} \end{array}\right]=\left[\begin{array}{l} p_{w} \pm q_{w} \\ \mathbf{p}_{v} \pm \mathbf{q}_{v} \end{array}\right]
p±q=[pwpv]±[qwqv]=[pw±qwpv±qv]
(2)乘积:
p
⊗
q
=
[
p
w
q
w
−
p
x
q
x
−
p
y
q
y
−
p
z
q
z
p
w
q
x
+
p
x
q
w
+
p
y
q
z
−
p
z
q
y
p
w
q
y
−
p
x
q
z
+
p
y
q
w
+
p
z
q
x
p
w
q
z
+
p
x
q
y
−
p
y
q
x
+
p
z
q
w
]
\mathbf{p} \otimes \mathbf{q}=\left[\begin{array}{l} p_{w} q_{w}-p_{x} q_{x}-p_{y} q_{y}-p_{z} q_{z} \\ p_{w} q_{x}+p_{x} q_{w}+p_{y} q_{z}-p_{z} q_{y} \\ p_{w} q_{y}-p_{x} q_{z}+p_{y} q_{w}+p_{z} q_{x} \\ p_{w} q_{z}+p_{x} q_{y}-p_{y} q_{x}+p_{z} q_{w} \end{array}\right]
p⊗q=⎣⎢⎢⎡pwqw−pxqx−pyqy−pzqzpwqx+pxqw+pyqz−pzqypwqy−pxqz+pyqw+pzqxpwqz+pxqy−pyqx+pzqw⎦⎥⎥⎤
上式也可以写成标量
和向量
相乘的形式:
p
⊗
q
=
[
p
w
q
w
−
p
v
⊤
q
v
p
w
q
v
+
q
w
p
v
+
p
v
×
q
v
]
\mathbf{p} \otimes \mathbf{q}=\left[\begin{array}{c} p_{w} q_{w}-\mathbf{p}_{v}^{\top} \mathbf{q}_{v} \\ p_{w} \mathbf{q}_{v}+q_{w} \mathbf{p}_{v}+\mathbf{p}_{v} \times \mathbf{q}_{v} \end{array}\right]
p⊗q=[pwqw−pv⊤qvpwqv+qwpv+pv×qv]
这里需要注意的是,四元数相乘没有交换律,
即
p
⊗
q
≠
q
⊗
p
\mathbf{p} \otimes \mathbf{q} \neq \mathbf{q} \otimes \mathbf{p}
p⊗q=q⊗p。四元数相乘也可以写成矩阵和向量相乘的形式
,即:
q
1
⊗
q
2
=
[
q
1
]
L
q
2
=
[
q
2
]
R
q
1
\mathbf{q}_1 \otimes \mathbf{q}_2 = [\mathbf{q}_1]_L\mathbf{q}_2 = [\mathbf{q}_2]_R\mathbf{q}_1
q1⊗q2=[q1]Lq2=[q2]Rq1
其中:
[
q
]
L
=
[
q
w
−
q
x
−
q
y
−
q
z
q
x
q
w
−
q
z
q
y
q
y
q
z
q
w
−
q
x
q
z
−
q
y
q
x
q
w
]
,
[
q
]
R
=
[
q
w
−
q
x
−
q
y
−
q
z
q
x
q
w
q
z
−
q
y
q
y
−
q
z
q
w
q
x
q
z
q
y
−
q
x
q
w
]
[\mathbf{q}]_{L}=\left[\begin{array}{cccc} q_{w} & -q_{x} & -q_{y} & -q_{z} \\ q_{x} & q_{w} & -q_{z} & q_{y} \\ q_{y} & q_{z} & q_{w} & -q_{x} \\ q_{z} & -q_{y} & q_{x} & q_{w} \end{array}\right], \quad[\mathbf{q}]_{R}=\left[\begin{array}{cccc} q_{w} & -q_{x} & -q_{y} & -q_{z} \\ q_{x} & q_{w} & q_{z} & -q_{y} \\ q_{y} & -q_{z} & q_{w} & q_{x} \\ q_{z} & q_{y} & -q_{x} & q_{w} \end{array}\right]
[q]L=⎣⎢⎢⎡qwqxqyqz−qxqwqz−qy−qy−qzqwqx−qzqy−qxqw⎦⎥⎥⎤,[q]R=⎣⎢⎢⎡qwqxqyqz−qxqw−qzqy−qyqzqw−qx−qz−qyqxqw⎦⎥⎥⎤
或者也可以写成下面的形式:
[
q
]
L
=
q
w
I
+
[
0
−
q
v
⊤
q
v
[
q
v
]
×
]
,
[
q
]
R
=
q
w
I
+
[
0
−
q
v
⊤
q
v
−
[
q
v
]
×
]
[\mathbf{q}]_{L}=q_{w} \mathbf{I}+\left[\begin{array}{cc} 0 & -\mathbf{q}_{v}^{\top} \\ \mathbf{q}_{v} & {\left[\mathbf{q}_{v}\right]_{\times}} \end{array}\right], \quad[\mathbf{q}]_{R}=q_{w} \mathbf{I}+\left[\begin{array}{cc} 0 & -\mathbf{q}_{v}^{\top} \\ \mathbf{q}_{v} & -\left[\mathbf{q}_{v}\right]_{\times} \end{array}\right]
[q]L=qwI+[0qv−qv⊤[qv]×],[q]R=qwI+[0qv−qv⊤−[qv]×]
反对称矩阵
[
∙
]
×
[\bullet]_{\times}
[∙]× 定义如下,
[
a
×
]
T
=
−
[
a
]
×
[\mathbf{a}_{\times}]^T=-[\mathbf{a}]_{\times}
[a×]T=−[a]×,并且有
[
a
]
×
b
=
a
×
b
[\mathbf{a}]_{\times}\mathbf{b}=\mathbf{a}\times \mathbf{b}
[a]×b=a×b。不过有时反对称矩阵符号也会写成
a
∧
\mathbf{a}^{\wedge}
a∧ 的形式。
[
a
]
×
≜
[
0
−
a
z
a
y
a
z
0
−
a
x
−
a
y
a
x
0
]
[\mathbf{a}]_{\times} \triangleq\left[\begin{array}{ccc} 0 & -a_{z} & a_{y} \\ a_{z} & 0 & -a_{x} \\ -a_{y} & a_{x} & 0 \end{array}\right]
[a]×≜⎣⎡0az−ay−az0axay−ax0⎦⎤
因此,四元数相乘可以写成:
q
⊗
x
⊗
p
=
(
q
⊗
x
)
⊗
p
=
[
p
]
R
[
q
]
L
x
=
q
⊗
(
x
⊗
p
)
=
[
q
]
L
[
p
]
R
x
\begin{aligned} \mathbf{q} \otimes \mathbf{x} \otimes \mathbf{p} &=(\mathbf{q} \otimes \mathbf{x}) \otimes \mathbf{p}=[\mathbf{p}]_{R}[\mathbf{q}]_{L} \mathbf{x} \\ &=\mathbf{q} \otimes(\mathbf{x} \otimes \mathbf{p})=[\mathbf{q}]_{L}[\mathbf{p}]_{R} \mathbf{x} \end{aligned}
q⊗x⊗p=(q⊗x)⊗p=[p]R[q]Lx=q⊗(x⊗p)=[q]L[p]Rx
可以得到;
[
p
]
R
[
q
]
L
=
[
q
]
L
[
p
]
R
[\mathbf{p}]_{R}[\mathbf{q}]_{L}=[\mathbf{q}]_{L}[\mathbf{p}]_{R}
[p]R[q]L=[q]L[p]R
(3)同一性:
q
1
=
1
=
[
1
0
v
]
\mathbf{q}_{1}=1=\left[\begin{array}{c} 1 \\ 0_{v} \end{array}\right]
q1=1=[10v]
(4)共轭四元数:
q
∗
≜
q
w
−
q
v
=
[
q
w
−
q
v
]
\mathbf{q}^{*} \triangleq q_{w}-\mathbf{q}_{v}=\left[\begin{array}{c} q_{w} \\ -\mathbf{q}_{v} \end{array}\right]
q∗≜qw−qv=[qw−qv]
可以得到:
q
⊗
q
∗
=
q
∗
⊗
q
=
q
w
2
+
q
x
2
+
q
y
2
+
q
z
2
=
[
q
w
2
+
q
x
2
+
q
y
2
+
q
z
2
0
v
]
\mathbf{q} \otimes \mathbf{q}^{*}=\mathbf{q}^{*} \otimes \mathbf{q}=q_{w}^{2}+q_{x}^{2}+q_{y}^{2}+q_{z}^{2}=\left[\begin{array}{c} q_{w}^{2}+q_{x}^{2}+q_{y}^{2}+q_{z}^{2} \\ \mathbf{0}_{v} \end{array}\right]
q⊗q∗=q∗⊗q=qw2+qx2+qy2+qz2=[qw2+qx2+qy2+qz20v]
同时有 ( p ⊗ q ) ∗ = q ∗ ⊗ p ∗ (\mathbf{p} \otimes \mathbf{q})^* = \mathbf{q}^* \otimes \mathbf{p}^* (p⊗q)∗=q∗⊗p∗。
(5)模:
∥ q ∥ ≜ q ⊗ q ∗ = q ∗ ⊗ q = q w 2 + q x 2 + q y 2 + q z 2 ∈ R \|\mathbf{q}\| \triangleq \sqrt{\mathbf{q} \otimes \mathbf{q}^{*}}=\sqrt{\mathbf{q}^{*} \otimes \mathbf{q}}=\sqrt{q_{w}^{2}+q_{x}^{2}+q_{y}^{2}+q_{z}^{2}} \in \mathbb{R} ∥q∥≜q⊗q∗=q∗⊗q=qw2+qx2+qy2+qz2∈R
可以得到: ∥ p ⊗ q ∥ = ∥ q ⊗ p ∥ = ∥ p ∥ ∥ q ∥ \|\mathbf{p} \otimes \mathbf{q}\|=\|\mathbf{q} \otimes \mathbf{p}\|=\|\mathbf{p}\|\|\mathbf{q}\| ∥p⊗q∥=∥q⊗p∥=∥p∥∥q∥
(6)逆:
q
−
1
=
q
∗
/
∥
q
∥
2
\mathbf{q}^{-1}=\mathbf{q}^{*} /\|\mathbf{q}\|^{2}
q−1=q∗/∥q∥2
1.3 四元数其余性质
(1)从上面的公式可以得到:
p
v
⊗
q
v
−
q
v
⊗
p
v
=
2
p
v
×
q
v
\mathbf{p}_{v} \otimes \mathbf{q}_{v}-\mathbf{q}_{v} \otimes \mathbf{p}_{v}=2 \mathbf{p}_{v} \times \mathbf{q}_{v}
pv⊗qv−qv⊗pv=2pv×qv
(2)纯四元数相乘:
p
v
⊗
q
v
=
−
p
v
⊤
q
v
+
p
v
×
q
v
=
[
−
p
v
⊤
q
v
p
v
×
q
v
]
\mathbf{p}_{v} \otimes \mathbf{q}_{v}=-\mathbf{p}_{v}^{\top} \mathbf{q}_{v}+\mathbf{p}_{v} \times \mathbf{q}_{v}=\left[\begin{array}{c} -\mathbf{p}_{v}^{\top} \mathbf{q}_{v} \\ \mathbf{p}_{v} \times \mathbf{q}_{v} \end{array}\right]
pv⊗qv=−pv⊤qv+pv×qv=[−pv⊤qvpv×qv]
可以得到:
q
v
⊗
q
v
=
−
q
v
⊤
q
v
=
−
∥
q
v
∥
2
\mathbf{q}_{v} \otimes \mathbf{q}_{v}=-\mathbf{q}_{v}^{\top} \mathbf{q}_{v}=-\left\|\mathbf{q}_{v}\right\|^{2}
qv⊗qv=−qv⊤qv=−∥qv∥2
如果
u
\mathbf{u}
u 是单位四元数,则有;
u
⊗
u
=
−
1
\mathbf{u} \otimes \mathbf{u}=-1
u⊗u=−1
(3)纯四元数幂次方:
如果
v
=
u
θ
\mathbf{v}=\mathbf{u}\theta
v=uθ 是纯四元数,且
u
\mathbf{u}
u 是单位四元数,
θ
=
∣
∣
v
∣
∣
\theta=||\mathbf{v}||
θ=∣∣v∣∣,则
v
\mathbf{v}
v 的幂次方形式为:
v
2
=
−
θ
2
,
v
3
=
−
u
θ
3
,
v
4
=
θ
4
,
v
5
=
u
θ
5
,
v
6
=
−
θ
6
\mathbf{v}^{2}=-\theta^{2} \quad, \quad \mathbf{v}^{3}=-\mathbf{u} \theta^{3} \quad, \quad \mathbf{v}^{4}=\theta^{4} \quad, \quad \mathbf{v}^{5}=\mathbf{u} \theta^{5} \quad, \quad \mathbf{v}^{6}=-\theta^{6}
v2=−θ2,v3=−uθ3,v4=θ4,v5=uθ5,v6=−θ6
对于纯单位四元数
u
\mathbf{u}
u,可以得到:
u
2
=
−
1
,
u
3
=
−
u
,
u
4
=
1
,
u
5
=
u
,
u
6
=
−
1
\mathbf{u}^{2}=-1 \quad, \quad \mathbf{u}^{3}=-\mathbf{u} \quad, \quad \mathbf{u}^{4}=1 \quad, \quad \mathbf{u}^{5}=\mathbf{u} \quad, \quad \mathbf{u}^{6}=-1
u2=−1,u3=−u,u4=1,u5=u,u6=−1
(4) 纯四元数指数:
四元数指数形式为:
e
q
≜
∑
k
=
0
∞
1
k
!
q
k
∈
H
e^{\mathbf{q}} \triangleq \sum_{k=0}^{\infty} \frac{1}{k !} \mathbf{q}^{k} \in \mathbb{H}
eq≜k=0∑∞k!1qk∈H
对于纯四元数
v
=
u
θ
\mathbf{v}=\mathbf{u}\theta
v=uθ ,
θ
=
∣
∣
v
∣
∣
\theta=||\mathbf{v}||
θ=∣∣v∣∣。根据幂次方形式可以得到其指数形式为:
e
u
θ
=
(
1
−
θ
2
2
!
+
θ
4
4
!
+
⋯
)
+
(
u
θ
−
u
θ
3
3
!
+
u
θ
5
5
!
+
⋯
)
e^{\mathbf{u} \theta}=\left(1-\frac{\theta^{2}}{2 !}+\frac{\theta^{4}}{4 !}+\cdots\right)+\left(\mathbf{u} \theta-\frac{\mathbf{u} \theta^{3}}{3 !}+\frac{\mathbf{u} \theta^{5}}{5 !}+\cdots\right)
euθ=(1−2!θ2+4!θ4+⋯)+(uθ−3!uθ3+5!uθ5+⋯)
最终,可以得到其表达形式为:
e
v
=
e
u
θ
=
cos
θ
+
u
sin
θ
=
[
cos
θ
u
sin
θ
]
e^{\mathbf{v}}=e^{\mathbf{u} \theta}=\cos \theta+\mathbf{u} \sin \theta=\left[\begin{array}{c} \cos \theta \\ \mathbf{u} \sin \theta \end{array}\right]
ev=euθ=cosθ+usinθ=[cosθusinθ]
这就是四元数指数形式对应的欧拉方程
。同时也可以得到:
e
−
v
=
(
e
v
)
∗
e^{-\mathbf{v}}=\left(e^{\mathbf{v}}\right)^{*}
e−v=(ev)∗
2. 旋转群、旋转矩阵、四元数
2.1 旋转方程与旋转群 S O ( 3 ) SO(3) SO(3)
假设有向量
x
\mathbf{x}
x 绕着轴
u
\mathbf{u}
u 旋转,角度为
ϕ
\phi
ϕ。则可以得到旋转后的向量
x
′
\mathbf{x}^{\prime}
x′为:
x
′
=
x
∥
+
x
⊥
cos
ϕ
+
(
u
×
x
)
sin
ϕ
\mathbf{x}^{\prime}=\mathbf{x}_{\|}+\mathbf{x}_{\perp} \cos \phi+(\mathbf{u} \times \mathbf{x}) \sin \phi
x′=x∥+x⊥cosϕ+(u×x)sinϕ
这就是旋转方程
,后面我们会进行证明。
下面给出旋转群
S
O
(
3
)
SO(3)
SO(3) 的定义:
S
O
(
3
)
:
{
r
:
R
3
→
R
3
/
∀
v
,
w
∈
R
3
,
∥
r
(
v
)
∥
=
∥
v
∥
,
r
(
v
)
×
r
(
w
)
=
r
(
v
×
w
)
}
S O(3):\left\{r: \mathbb{R}^{3} \rightarrow \mathbb{R}^{3} / \forall \mathbf{v}, \mathbf{w} \in \mathbb{R}^{3},\|r(\mathbf{v})\|=\|\mathbf{v}\|, r(\mathbf{v}) \times r(\mathbf{w})=r(\mathbf{v} \times \mathbf{w})\right\}
SO(3):{r:R3→R3/∀v,w∈R3,∥r(v)∥=∥v∥,r(v)×r(w)=r(v×w)}
可以看到,一个向量
v
\mathbf{v}
v 旋转后,其模并不改变,向量间的角度也不会改变
。
2.2 旋转群和旋转矩阵
(1)旋转矩阵基本性质
已知一个旋转矩阵
R
∈
R
3
×
3
\mathbf{R}\in\mathbb{R}^{3\times3}
R∈R3×3,定义旋转操作为
r
(
v
)
=
R
v
r(\mathbf{v})=\mathbf{R}\mathbf{v}
r(v)=Rv。因为旋转不改变向量的模,可以得到:
(
R
v
)
⊤
(
R
v
)
=
v
⊤
R
⊤
R
v
=
v
⊤
v
(\mathbf{R} \mathbf{v})^{\top}(\mathbf{R} \mathbf{v})=\mathbf{v}^{\top} \mathbf{R}^{\top} \mathbf{R} \mathbf{v}=\mathbf{v}^{\top} \mathbf{v}
(Rv)⊤(Rv)=v⊤R⊤Rv=v⊤v
得到旋转矩阵的正交性
:
R
⊤
R
=
I
=
R
R
⊤
\mathbf{R}^{\top} \mathbf{R}=\mathbf{I}=\mathbf{R} \mathbf{R}^{\top}
R⊤R=I=RR⊤
从上式可以得到,旋转矩阵的逆是其的转置
R
−
1
=
R
⊤
\mathbf{R}^{-1}=\mathbf{R}^{\top}
R−1=R⊤,旋转矩阵行列式
为
det
(
R
)
=
1
\operatorname{det}(\mathbf{R})=1
det(R)=1。
(2)旋转矩阵指数映射
上式两边求导:
d d t ( R ⊤ R ) = R ˙ ⊤ R + R ⊤ R ˙ = 0 \frac{d}{d t}\left(\mathbf{R}^{\top} \mathbf{R}\right)=\dot{\mathbf{R}}^{\top} \mathbf{R}+\mathbf{R}^{\top} \dot{\mathbf{R}}=0 dtd(R⊤R)=R˙⊤R+R⊤R˙=0
可以得到;
R ⊤ R ˙ = − ( R ⊤ R ˙ ) ⊤ \mathbf{R}^{\top} \dot{\mathbf{R}}=-\left(\mathbf{R}^{\top} \dot{\mathbf{R}}\right)^{\top} R⊤R˙=−(R⊤R˙)⊤
可以看出
R
⊤
R
˙
\mathbf{R}^{\top}\dot{\mathbf{R}}
R⊤R˙ 是一个反对称矩阵
。反对称的
3
×
3
3\times3
3×3 矩阵可以用符号
s
o
(
3
)
\mathfrak{so}(3)
so(3)表示,也称为群
S
O
3
SO{3}
SO3 的李代数
。
定义向量
ω
∈
R
3
↔
[
ω
]
×
∈
s
o
(
3
)
\boldsymbol{\omega} \in \mathbb{R}^{3} \leftrightarrow[\boldsymbol{\omega}]_{\times} \in \mathfrak{s o}(3)
ω∈R3↔[ω]×∈so(3),则上式可以写为:
R
˙
=
R
[
ω
]
×
\dot{\mathbf{R}}=\mathbf{R}[\boldsymbol{\omega}]_{\times}
R˙=R[ω]×
在原点,旋转矩阵
R
=
I
\mathbf{R}=\mathbf{I}
R=I,上式简化为
R
˙
=
[
ω
]
×
\dot{\mathbf{R}}=[\boldsymbol{\omega}]_{\times}
R˙=[ω]×。表示李代数为旋转矩阵在原点的导数
。
如果
ω
\boldsymbol{\omega}
ω 是常数,上式微分方程的积分形式可以写为:
R
(
t
)
=
R
(
0
)
e
[
ω
]
×
t
=
R
(
0
)
e
[
ω
t
]
×
\mathbf{R}(t)=\mathbf{R}(0) e^{[\boldsymbol{\omega}]_{\times} t}=\mathbf{R}(0) e^{[\omega t]_{\times}}
R(t)=R(0)e[ω]×t=R(0)e[ωt]×
定义向量
ϕ
≜
ω
Δ
t
\phi \triangleq \omega \Delta t
ϕ≜ωΔt,则旋转矩阵
R
\mathbf{R}
R 可以写成:
R
=
e
[
ϕ
]
×
\mathbf{R}=e^{[\phi]_{\times}}
R=e[ϕ]×
这就是旋转矩阵的指数映射
,从李代数
s
o
(
3
)
\mathfrak{so}(3)
so(3) 到群
S
O
(
3
)
SO(3)
SO(3):
exp : s o ( 3 ) → S O ( 3 ) ; [ ϕ ] × ↦ exp ( [ ϕ ] × ) = e [ ϕ ] × \exp : \mathfrak{s o}(3) \rightarrow S O(3) ;[\boldsymbol{\phi}]_{\times} \mapsto \exp \left([\boldsymbol{\phi}]_{\times}\right)=e^{[\boldsymbol{\phi}]_{\times}} exp:so(3)→SO(3);[ϕ]×↦exp([ϕ]×)=e[ϕ]×
也可以写成从实数
R
3
\mathbb{R}^3
R3 到群
S
O
(
3
)
SO(3)
SO(3)的映射。
Exp
:
R
3
→
S
O
(
3
)
;
ϕ
↦
Exp
(
ϕ
)
=
e
[
ϕ
]
×
\operatorname{Exp}: \mathbb{R}^{3} \rightarrow S O(3) ; \phi \mapsto \operatorname{Exp}(\phi)=e^{[\phi]_{\times}}
Exp:R3→SO(3);ϕ↦Exp(ϕ)=e[ϕ]×
实数、李代数数、旋转群三者之间关系如下图所示
:
(3) 旋转矩阵与旋转向量:罗德里格斯公式
已知旋转
R
\mathbf{R}
R 为:
且旋转轴
u
\mathbf{u}
u 满足
[
u
]
×
2
=
u
u
⊤
−
I
,
[
u
]
×
3
=
−
u
×
[\mathbf{u}]_{\times}^2=\mathbf{u}\mathbf{u}^{\top}-\mathbf{I},[\mathbf{u}]_{\times}^3=-\mathbf{u}_{\times}
[u]×2=uu⊤−I,[u]×3=−u×。最终上式可以写成如下形式,也被称为罗德里格斯公式
。
R
=
I
+
sin
ϕ
[
u
]
×
+
(
1
−
cos
ϕ
)
[
u
]
×
2
\mathbf{R} = \mathbf{I} + \sin\phi [\mathbf{u}]_{\times} + (1 − \cos\phi) [\mathbf{u}]_{\times}^2
R=I+sinϕ[u]×+(1−cosϕ)[u]×2
指数映射的逆运算为对数映射,定义为:
log
:
S
O
(
3
)
→
s
o
(
3
)
;
R
↦
log
(
R
)
=
[
u
ϕ
]
×
\log : SO(3) \rightarrow \mathfrak{s o}(3) ; \mathbf{R} \mapsto \log(\mathbf{R}) =[\mathbf{u}\phi]_{\times}
log:SO(3)→so(3);R↦log(R)=[uϕ]×
其中
:
ϕ
=
arccos
(
t
r
a
c
e
(
R
)
−
1
2
)
\phi=\arccos(\frac{trace(\mathbf{R})-1}{2})
ϕ=arccos(2trace(R)−1)
u = ( R − R ⊤ ) ∨ 2 sin ϕ \mathbf{u}=\frac{(\mathbf{R}-\mathbf{R}^{\top})^∨}{2\sin\phi} u=2sinϕ(R−R⊤)∨
其中 [ ∙ ] ∨ [\bullet]^∨ [∙]∨ 是反对称矩阵的逆操作。
(4) 旋转运算
回到最开始的问题, 向量
x
\mathbf{x}
x 绕轴
u
\mathbf{u}
u 转动角度
ϕ
\phi
ϕ 之后的向量
x
′
\mathbf{x^{\prime}}
x′ 为:
x
′
=
R
x
=
(
I
+
sin
ϕ
[
u
]
×
+
(
1
−
cos
ϕ
)
[
u
]
×
2
)
x
=
x
+
sin
ϕ
[
u
]
×
x
+
(
1
−
cos
ϕ
)
[
u
]
×
2
x
=
x
+
sin
ϕ
(
u
×
x
)
+
(
1
−
cos
ϕ
)
(
u
u
⊤
−
I
)
x
=
x
∥
+
x
⊥
+
sin
ϕ
(
u
×
x
)
−
(
1
−
cos
ϕ
)
x
⊥
=
x
∥
+
(
u
×
x
)
sin
ϕ
+
x
⊥
cos
ϕ
\begin{aligned} \mathbf{x}^{\prime} &=\mathbf{R} \mathbf{x} \\ &=\left(\mathbf{I}+\sin \phi[\mathbf{u}]_{\times}+(1-\cos \phi)[\mathbf{u}]_{\times}^{2}\right) \mathbf{x} \\ &=\mathbf{x}+\sin \phi[\mathbf{u}]_{\times} \mathbf{x}+(1-\cos \phi)[\mathbf{u}]_{\times}^{2} \mathbf{x} \\ &=\mathbf{x}+\sin \phi(\mathbf{u} \times \mathbf{x})+(1-\cos \phi)\left(\mathbf{u} \mathbf{u}^{\top}-\mathbf{I}\right) \mathbf{x} \\ &=\mathbf{x}_{\|}+\mathbf{x}_{\perp}+\sin \phi(\mathbf{u} \times \mathbf{x})-(1-\cos \phi) \mathbf{x}_{\perp} \\ &=\mathbf{x}_{\|}+(\mathbf{u} \times \mathbf{x}) \sin \phi+\mathbf{x}_{\perp} \cos \phi \end{aligned}
x′=Rx=(I+sinϕ[u]×+(1−cosϕ)[u]×2)x=x+sinϕ[u]×x+(1−cosϕ)[u]×2x=x+sinϕ(u×x)+(1−cosϕ)(uu⊤−I)x=x∥+x⊥+sinϕ(u×x)−(1−cosϕ)x⊥=x∥+(u×x)sinϕ+x⊥cosϕ
2.3 旋转群与四元数
(1)四元数旋转基本性质
四元数旋转公式为:
r
(
v
)
=
q
⊗
v
⊗
q
∗
r(\mathbf{v})=\mathbf{q}\otimes\mathbf{v}\otimes\mathbf{q}^*
r(v)=q⊗v⊗q∗
旋转并不改变向量的长度,两边取模可得:
∣
∣
q
⊗
v
⊗
q
∗
∣
∣
=
∣
∣
q
∣
∣
2
∣
∣
v
∣
∣
=
∣
∣
v
∣
∣
||\mathbf{q}\otimes\mathbf{v}\otimes\mathbf{q}^*||=||\mathbf{q}||^2||\mathbf{v}||=||\mathbf{v}||
∣∣q⊗v⊗q∗∣∣=∣∣q∣∣2∣∣v∣∣=∣∣v∣∣
最终可以得到:
q
⊗
q
∗
=
1
=
q
∗
⊗
q
\mathbf{q}\otimes\mathbf{q}^*=1=\mathbf{q}^*\otimes\mathbf{q}
q⊗q∗=1=q∗⊗q
与上一节关于旋转矩阵性质十分相似,即 R ⊤ R = I = R R ⊤ \mathbf{R}^{\top}\mathbf{R}=\mathbf{I}=\mathbf{R}\mathbf{R}^{\top} R⊤R=I=RR⊤。同样地,单位四元数也构成了一个群,用 S 3 S^3 S3表示。
(2)四元数指数映射
假设有一个四元数
q
∈
S
3
\mathbf{q}\in S^3
q∈S3,满足
q
∗
⊗
q
=
1
\mathbf{q}^*\otimes\mathbf{q}=1
q∗⊗q=1,两边求导可得:
d
(
q
∗
⊗
q
)
d
t
=
q
˙
∗
⊗
q
+
q
∗
⊗
q
˙
=
0
\frac{d(\mathbf{q}^*\otimes\mathbf{q})}{d t}=\dot{\mathbf{q}}^{*} \otimes\mathbf{q}+\mathbf{q}^{*} \otimes\dot{\mathbf{q}}=0
dtd(q∗⊗q)=q˙∗⊗q+q∗⊗q˙=0
可以得到:
q
∗
⊗
q
˙
=
−
(
q
˙
∗
⊗
q
)
=
−
(
q
∗
⊗
q
˙
)
∗
\mathbf{q}^{*} \otimes\dot{\mathbf{q}}=-(\dot{\mathbf{q}}^{*} \otimes\mathbf{q})=-(\mathbf{q}^{*} \otimes\dot{\mathbf{q}})^*
q∗⊗q˙=−(q˙∗⊗q)=−(q∗⊗q˙)∗
可以证明,
q
∗
⊗
q
˙
\mathbf{q}^{*} \otimes\dot{\mathbf{q}}
q∗⊗q˙ 是一个纯四元数,使用
Ω
\mathbf{\Omega}
Ω 表示纯四元数,则有:
q
∗
⊗
q
˙
=
Ω
=
[
0
Ω
]
\mathbf{q}^{*} \otimes\dot{\mathbf{q}}=\mathbf{\Omega}=\left[\begin{array}{c} 0 \\ \mathbf{\Omega} \end{array}\right]
q∗⊗q˙=Ω=[0Ω]
最终可以得到,四元数导数为:
q
˙
=
q
⊗
Ω
\dot\mathbf{q}=\mathbf{q}\otimes\mathbf{\Omega}
q˙=q⊗Ω
在原点,四元数
q
=
1
\mathbf{q}=1
q=1。此时得到李代数
q
˙
=
Ω
\dot\mathbf{q}=\mathbf{\Omega}
q˙=Ω。积分可得:
q
(
t
)
=
q
(
0
)
⊗
e
Ω
t
\mathbf{q}(t)=\mathbf{q}(0)\otimes e^{\mathbf{\Omega} t}
q(t)=q(0)⊗eΩt
此时,可以得到四元数的指数映射方程
为:
exp
:
H
p
→
S
3
;
V
↦
exp
(
V
)
=
e
V
\exp : \mathbb{H}_p \rightarrow S^3 ; \boldsymbol{V} \mapsto \exp \left(\boldsymbol{V}\right)=e^{\boldsymbol{V}}
exp:Hp→S3;V↦exp(V)=eV
同样地,也可以写成实数
R
\mathbb{R}
R 到 群
S
3
S^3
S3 的映射。
Exp
:
R
3
→
S
3
;
ϕ
↦
Exp
(
ϕ
)
=
e
ϕ
/
2
\operatorname{Exp}: \mathbb{R}^3 \rightarrow S^3 ; \boldsymbol{\phi} \mapsto \operatorname{Exp} \left(\boldsymbol{\phi}\right)=e^{\boldsymbol{ \boldsymbol{\phi}/2}}
Exp:R3→S3;ϕ↦Exp(ϕ)=eϕ/2
实数、四元数、旋转群
三者之间关系为;
有时,为了方便,会引入角速度向量
ω
=
2
Ω
∈
R
3
\mathbf{\omega}=2\mathbf{\Omega}\in\mathbb{R}^3
ω=2Ω∈R3,则四元数及其导数
可以写为:
q
˙
=
1
2
q
⊗
ω
\dot{\mathbf{q}}=\frac{1}{2}\mathbf{q}\otimes\mathbf{\omega}
q˙=21q⊗ω
q = e ω t / 2 \mathbf{q}=e^{\mathbf{\omega}t/2} q=eωt/2
(3) 四元数与旋转向量
令
ϕ
=
ϕ
u
\boldsymbol{\phi}=\phi\mathbf{u}
ϕ=ϕu,其中
u
\mathbf{u}
u为旋转向量,
ϕ
\phi
ϕ 为旋转角度。则四元数与旋转向量方程
为:
q
=
Exp
(
ϕ
u
)
=
e
ϕ
u
/
2
=
[
cos
(
ϕ
/
2
)
u
sin
(
ϕ
/
2
)
]
\mathbf{q}=\operatorname{Exp}(\phi\mathbf{u}) =e^{\phi\mathbf{u}/2}=\left[\begin{array}{c} \cos(\phi/2)\\ \mathbf{u}\sin(\phi/2) \end{array}\right]
q=Exp(ϕu)=eϕu/2=[cos(ϕ/2)usin(ϕ/2)]
从四元数求旋转向量和角度
公式为:
ϕ = 2 arctan ( ∣ ∣ q v ∣ ∣ , q w ) \phi=2\arctan(||\mathbf{q}_v||,q_w) ϕ=2arctan(∣∣qv∣∣,qw)
u = q v / ∣ ∣ q v ∣ ∣ \mathbf{u}=\mathbf{q}_v/||\mathbf{q}_v|| u=qv/∣∣qv∣∣
(4) 旋转运算
回到最开始的问题,向量
x
\mathbf{x}
x 经过四元数旋转后的向量
x
′
\mathbf{x}^{\prime}
x′为:
x
′
=
q
⊗
x
⊗
q
∗
=
(
cos
ϕ
2
+
u
sin
ϕ
2
)
⊗
(
0
+
x
)
⊗
(
cos
ϕ
2
−
u
sin
ϕ
2
)
=
x
cos
2
ϕ
2
+
(
u
⊗
x
−
x
⊗
u
)
sin
ϕ
2
cos
ϕ
2
−
u
⊗
x
⊗
u
sin
2
ϕ
2
=
x
cos
2
ϕ
2
+
2
(
u
×
x
)
sin
ϕ
2
cos
ϕ
2
−
(
x
(
u
⊤
u
)
−
2
u
(
u
⊤
x
)
)
sin
2
ϕ
2
=
x
(
cos
2
ϕ
2
−
sin
2
ϕ
2
)
+
(
u
×
x
)
(
2
sin
ϕ
2
cos
ϕ
2
)
+
u
(
u
⊤
x
)
(
2
sin
2
ϕ
2
)
=
x
cos
ϕ
+
(
u
×
x
)
sin
ϕ
+
u
(
u
⊤
x
)
(
1
−
cos
ϕ
)
=
(
x
−
u
u
⊤
x
)
cos
ϕ
+
(
u
×
x
)
sin
ϕ
+
u
u
⊤
x
=
x
⊥
cos
ϕ
+
(
u
×
x
)
sin
ϕ
+
x
∥
\begin{aligned} \mathbf{x}^{\prime} &=\mathbf{q} \otimes \mathbf{x} \otimes \mathbf{q}^{*} \\ &=\left(\cos \frac{\phi}{2}+\mathbf{u} \sin \frac{\phi}{2}\right) \otimes(0+\mathbf{x}) \otimes\left(\cos \frac{\phi}{2}-\mathbf{u} \sin \frac{\phi}{2}\right) \\ &=\mathbf{x} \cos ^{2} \frac{\phi}{2}+(\mathbf{u} \otimes \mathbf{x}-\mathbf{x} \otimes \mathbf{u}) \sin \frac{\phi}{2} \cos \frac{\phi}{2}-\mathbf{u} \otimes \mathbf{x} \otimes \mathbf{u} \sin ^{2} \frac{\phi}{2} \\ &=\mathbf{x} \cos ^{2} \frac{\phi}{2}+2(\mathbf{u} \times \mathbf{x}) \sin \frac{\phi}{2} \cos \frac{\phi}{2}-\left(\mathbf{x}\left(\mathbf{u}^{\top} \mathbf{u}\right)-2 \mathbf{u}\left(\mathbf{u}^{\top} \mathbf{x}\right)\right) \sin ^{2} \frac{\phi}{2} \\ &=\mathbf{x}\left(\cos ^{2} \frac{\phi}{2}-\sin ^{2} \frac{\phi}{2}\right)+(\mathbf{u} \times \mathbf{x})\left(2 \sin \frac{\phi}{2} \cos \frac{\phi}{2}\right)+\mathbf{u}\left(\mathbf{u}^{\top} \mathbf{x}\right)\left(2 \sin ^{2} \frac{\phi}{2}\right) \\ &=\mathbf{x} \cos \phi+(\mathbf{u} \times \mathbf{x}) \sin \phi+\mathbf{u}\left(\mathbf{u}^{\top} \mathbf{x}\right)(1-\cos \phi) \\ &=\left(\mathbf{x}-\mathbf{u} \mathbf{u}^{\top} \mathbf{x}\right) \cos \phi+(\mathbf{u} \times \mathbf{x}) \sin \phi+\mathbf{u} \mathbf{u}^{\top} \mathbf{x} \\ &=\mathbf{x}_{\perp} \cos \phi+(\mathbf{u} \times \mathbf{x}) \sin \phi+\mathbf{x}_{\|} \end{aligned}
x′=q⊗x⊗q∗=(cos2ϕ+usin2ϕ)⊗(0+x)⊗(cos2ϕ−usin2ϕ)=xcos22ϕ+(u⊗x−x⊗u)sin2ϕcos2ϕ−u⊗x⊗usin22ϕ=xcos22ϕ+2(u×x)sin2ϕcos2ϕ−(x(u⊤u)−2u(u⊤x))sin22ϕ=x(cos22ϕ−sin22ϕ)+(u×x)(2sin2ϕcos2ϕ)+u(u⊤x)(2sin22ϕ)=xcosϕ+(u×x)sinϕ+u(u⊤x)(1−cosϕ)=(x−uu⊤x)cosϕ+(u×x)sinϕ+uu⊤x=x⊥cosϕ+(u×x)sinϕ+x∥
2.4 旋转矩阵与四元数
无论向量
x
\mathbf{x}
x 是用四元数进行旋转,还是旋转矩阵进行旋转,旋转后的结果都是一致的,则有:
q
⊗
x
⊗
q
∗
=
R
x
\mathbf{q}\otimes\mathbf{x}\otimes\mathbf{q}^*=\mathbf{R}\mathbf{x}
q⊗x⊗q∗=Rx
因此可以得到四元数到旋转矩阵的
转换关系为:
R
=
[
q
w
2
+
q
x
2
−
q
y
2
−
q
z
2
2
(
q
x
q
y
−
q
w
q
z
)
2
(
q
x
q
z
+
q
w
q
y
)
2
(
q
x
q
y
+
q
w
q
z
)
q
w
2
−
q
x
2
+
q
y
2
−
q
z
2
2
(
q
y
q
z
−
q
w
q
x
)
2
(
q
x
q
z
−
q
w
q
y
)
2
(
q
y
q
z
+
q
w
q
x
)
q
w
2
−
q
x
2
−
q
y
2
+
q
z
2
]
\mathbf{R}=\left[\begin{array}{ccc} q_{w}^{2}+q_{x}^{2}-q_{y}^{2}-q_{z}^{2} & 2\left(q_{x} q_{y}-q_{w} q_{z}\right) & 2\left(q_{x} q_{z}+q_{w} q_{y}\right) \\ 2\left(q_{x} q_{y}+q_{w} q_{z}\right) & q_{w}^{2}-q_{x}^{2}+q_{y}^{2}-q_{z}^{2} & 2\left(q_{y} q_{z}-q_{w} q_{x}\right) \\ 2\left(q_{x} q_{z}-q_{w} q_{y}\right) & 2\left(q_{y} q_{z}+q_{w} q_{x}\right) & q_{w}^{2}-q_{x}^{2}-q_{y}^{2}+q_{z}^{2} \end{array}\right]
R=⎣⎡qw2+qx2−qy2−qz22(qxqy+qwqz)2(qxqz−qwqy)2(qxqy−qwqz)qw2−qx2+qy2−qz22(qyqz+qwqx)2(qxqz+qwqy)2(qyqz−qwqx)qw2−qx2−qy2+qz2⎦⎤
为了书写方便,旋转矩阵也可以写成
:
R
=
(
q
w
2
−
q
v
⊤
q
v
)
I
+
2
q
v
q
v
⊤
+
2
q
w
[
q
v
]
×
\mathbf{R}=(q_w^2-\mathbf{q}_v^{\top}\mathbf{q}_v)\mathbf{I}+2\mathbf{q}_v\mathbf{q}_v^{\top}+2\mathbf{q}_w[\mathbf{q}_v]_{\times}
R=(qw2−qv⊤qv)I+2qvqv⊤+2qw[qv]×
旋转矩阵
R
\mathbf{R}
R 和四元数
q
\mathbf{q}
q 总结如下表所示。
3. 李群扰动导数
3.1 李群加减运算符
在实数空间 R n \mathbb{R}^{n} Rn 中,加减运算符使用 + , − +,- +,− 表示。但是在李群 S O ( 3 ) SO(3) SO(3) 上,定义的加减运算符为 ⊕ , ⊖ \oplus, \ominus ⊕,⊖。
加法运算为: S = R ⊕ θ ≜ R ∘ Exp ( θ ) R , S ∈ S O ( 3 ) , θ ∈ R 3 \mathrm{S}=\mathrm{R} \oplus \boldsymbol{\theta} \triangleq \mathrm{R} \circ \operatorname{Exp}(\boldsymbol{\theta}) \quad \mathrm{R}, \mathrm{S} \in S O(3), \boldsymbol{\theta} \in \mathbb{R}^{3} S=R⊕θ≜R∘Exp(θ)R,S∈SO(3),θ∈R3
可以发现运算符满足任意
S
O
(
3
)
SO(3)
SO(3) 的表示,对于四元数可旋转矩阵,可以写成:
q
S
=
q
R
⊕
θ
=
q
R
⊗
Exp
(
θ
)
R
S
=
R
R
⊕
θ
=
R
R
⋅
Exp
(
θ
)
\begin{aligned} \mathbf{q}_{\mathrm{S}} &=\mathbf{q}_{\mathrm{R}} \oplus \boldsymbol{\theta}=\mathbf{q}_{\mathrm{R}} \otimes \operatorname{Exp}(\boldsymbol{\theta}) \\ \mathbf{R}_{\mathrm{S}} &=\mathbf{R}_{\mathrm{R}} \oplus \boldsymbol{\theta}=\mathbf{R}_{\mathrm{R}} \cdot \operatorname{Exp}(\boldsymbol{\theta}) \end{aligned}
qSRS=qR⊕θ=qR⊗Exp(θ)=RR⊕θ=RR⋅Exp(θ)
减法运算为:
θ
=
S
⊖
R
≜
log
(
R
−
1
∘
S
)
R
,
S
∈
S
O
(
3
)
,
θ
∈
R
3
\boldsymbol{\theta}=\mathrm{S} \ominus \mathrm{R} \triangleq \log \left(\mathrm{R}^{-1} \circ \mathrm{S}\right) \quad \mathrm{R}, \mathrm{S} \in S O(3), \boldsymbol{\theta} \in \mathbb{R}^{3}
θ=S⊖R≜log(R−1∘S)R,S∈SO(3),θ∈R3
对于旋转矩阵和四元数,其减法运算为:
θ
=
q
s
⊖
q
R
=
log
(
q
R
∗
⊗
q
S
)
θ
=
R
s
⊖
R
R
=
log
(
R
R
⊤
R
S
)
\begin{array}{l} \boldsymbol{\theta}=\mathbf{q}_{\mathrm{s}} \ominus \mathbf{q}_{\mathrm{R}}=\log \left(\mathbf{q}_{R}^{*} \otimes \mathbf{q}_{\mathrm{S}}\right) \\ \boldsymbol{\theta}=\mathbf{R}_{\mathrm{s}} \ominus \mathbf{R}_{\mathrm{R}}=\log \left(\mathbf{R}_{\mathrm{R}}^{\top} \mathbf{R}_{\mathrm{S}}\right) \end{array}
θ=qs⊖qR=log(qR∗⊗qS)θ=Rs⊖RR=log(RR⊤RS)
3.2 李群四种导数定义
(1)向量函数与向量的导数
已知函数
f
:
R
m
→
R
n
f: \mathbb{R}^{m} \rightarrow \mathbb{R}^{n}
f:Rm→Rn,则函数导数为:
∂
f
(
x
)
∂
x
≜
lim
δ
x
→
0
f
(
x
+
δ
x
)
−
f
(
x
)
δ
x
∈
R
n
×
m
\frac{\partial f(\mathbf{x})}{\partial \mathbf{x}} \triangleq \lim _{\delta \mathbf{x} \rightarrow 0} \frac{f(\mathbf{x}+\delta \mathbf{x})-f(\mathbf{x})}{\delta \mathbf{x}} \quad \in \mathbb{R}^{n \times m}
∂x∂f(x)≜δx→0limδxf(x+δx)−f(x)∈Rn×m
则欧拉积分公式为:
f
(
x
+
Δ
x
)
≈
f
(
x
)
+
∂
f
(
x
)
∂
x
Δ
x
∈
R
n
f(\mathbf{x}+\Delta \mathbf{x}) \approx f(\mathbf{x})+\frac{\partial f(\mathbf{x})}{\partial \mathbf{x}} \Delta \mathbf{x} \quad \in \mathbb{R}^{n}
f(x+Δx)≈f(x)+∂x∂f(x)Δx∈Rn
(2)李群函数与李群的导数
已知函数
f
:
S
O
(
3
)
→
S
O
(
3
)
f: SO(3) \rightarrow SO(3)
f:SO(3)→SO(3),则函数导数为:
∂
f
(
R
)
∂
θ
≜
lim
δ
θ
→
0
f
(
R
⊕
δ
θ
)
⊖
f
(
R
)
δ
θ
∈
R
3
×
3
=
lim
δ
θ
→
0
log
(
f
−
1
(
R
)
f
(
R
Exp
(
δ
θ
)
)
)
δ
θ
\begin{aligned} \frac{\partial f(\mathrm{R})}{\partial \boldsymbol{\theta}} & \triangleq \lim _{\delta \boldsymbol{\theta} \rightarrow 0} \frac{f(\mathrm{R} \oplus \delta \boldsymbol{\theta}) \ominus f(\mathrm{R})}{\delta \boldsymbol{\theta}} & \in \mathbb{R}^{3 \times 3} \\ &=\lim _{\delta \boldsymbol{\theta} \rightarrow 0} \frac{\log \left(f^{-1}(\mathrm{R}) f(\mathrm{R} \operatorname{Exp}(\delta \boldsymbol{\theta}))\right)}{\delta \boldsymbol{\theta}} \end{aligned}
∂θ∂f(R)≜δθ→0limδθf(R⊕δθ)⊖f(R)=δθ→0limδθlog(f−1(R)f(RExp(δθ)))∈R3×3
则欧拉积分公式为:
f
(
R
⊕
Δ
θ
)
≈
f
(
R
)
⊕
∂
f
(
R
)
∂
θ
Δ
θ
≜
f
(
R
)
Exp
(
∂
f
(
R
)
∂
θ
Δ
θ
)
∈
S
O
(
3
)
f(\mathrm{R} \oplus \Delta \boldsymbol{\theta}) \approx f(\mathrm{R}) \oplus \frac{\partial f(\mathrm{R})}{\partial \boldsymbol{\theta}} \Delta \boldsymbol{\theta} \triangleq f(\mathrm{R}) \operatorname{Exp}\left(\frac{\partial f(\mathrm{R})}{\partial \boldsymbol{\theta}} \Delta \boldsymbol{\theta}\right) \quad \in S O(3)
f(R⊕Δθ)≈f(R)⊕∂θ∂f(R)Δθ≜f(R)Exp(∂θ∂f(R)Δθ)∈SO(3)
(3)李群函数与向量的导数
已知函数
f
:
R
m
→
S
O
(
3
)
f: \mathbb{R}^{m} \rightarrow SO(3)
f:Rm→SO(3),则函数导数为:
∂
f
(
x
)
∂
x
≜
lim
δ
x
→
0
f
(
x
+
δ
x
)
⊖
f
(
x
)
δ
x
∈
R
3
×
m
=
lim
δ
x
→
0
log
(
f
−
1
(
x
)
f
(
x
+
δ
x
)
)
δ
x
\begin{aligned} \frac{\partial f(\mathbf{x})}{\partial \mathbf{x}} & \triangleq \lim _{\delta \mathbf{x} \rightarrow 0} \frac{f(\mathbf{x}+\delta \mathbf{x}) \ominus f(\mathbf{x})}{\delta \mathbf{x}} \quad \in \mathbb{R}^{3 \times m} \\ &=\lim _{\delta \mathbf{x} \rightarrow 0} \frac{\log \left(f^{-1}(\mathbf{x}) f(\mathbf{x}+\delta \mathbf{x})\right)}{\delta \mathbf{x}} \end{aligned}
∂x∂f(x)≜δx→0limδxf(x+δx)⊖f(x)∈R3×m=δx→0limδxlog(f−1(x)f(x+δx))
则欧拉积分公式为:
f
(
x
+
Δ
x
)
≈
f
(
x
)
⊕
∂
f
(
x
)
∂
x
Δ
x
≜
f
(
x
)
Exp
(
∂
f
(
x
)
∂
x
Δ
x
)
∈
S
O
(
3
)
f(\mathbf{x}+\Delta \mathbf{x}) \approx f(\mathbf{x}) \oplus \frac{\partial f(\mathbf{x})}{\partial \mathbf{x}} \Delta \mathbf{x} \triangleq f(\mathbf{x}) \operatorname{Exp}\left(\frac{\partial f(\mathbf{x})}{\partial \mathbf{x}} \Delta \mathbf{x}\right) \quad \in S O(3)
f(x+Δx)≈f(x)⊕∂x∂f(x)Δx≜f(x)Exp(∂x∂f(x)Δx)∈SO(3)
(4)向量函数与李群的导数
已知函数
f
:
S
O
(
3
)
→
R
3
f: SO(3) \rightarrow \mathbb{R}^3
f:SO(3)→R3,则函数导数为:
∂
f
(
R
)
∂
θ
≜
lim
δ
θ
→
0
f
(
R
⊕
δ
θ
)
−
f
(
R
)
δ
θ
=
lim
δ
θ
→
0
f
(
R
Exp
(
δ
θ
)
)
−
f
(
R
)
δ
θ
\begin{aligned} \frac{\partial f(\mathrm{R})}{\partial \boldsymbol{\theta}} & \triangleq \lim _{\delta \theta \rightarrow 0} \frac{f(\mathrm{R} \oplus \delta \boldsymbol{\theta})-f(\mathrm{R})}{\delta \boldsymbol{\theta}} \\ &=\lim _{\delta \boldsymbol{\theta} \rightarrow 0} \frac{f(\mathrm{R} \operatorname{Exp}(\delta \boldsymbol{\theta}))-f(\mathrm{R})}{\delta \boldsymbol{\theta}} \end{aligned}
∂θ∂f(R)≜δθ→0limδθf(R⊕δθ)−f(R)=δθ→0limδθf(RExp(δθ))−f(R)
则欧拉积分公式为:
f
(
R
⊕
δ
θ
)
≈
f
(
R
)
+
∂
f
(
R
)
∂
θ
Δ
θ
≜
f
(
R
)
+
Exp
(
∂
f
(
R
)
∂
θ
Δ
θ
)
∈
R
3
f(\mathrm{R} \oplus \delta \boldsymbol{\theta}) \approx f(\mathrm{R})+\frac{\partial f(\mathrm{R})}{\partial \boldsymbol{\theta}} \Delta \boldsymbol{\theta} \triangleq f(\mathrm{R})+\operatorname{Exp}\left(\frac{\partial f(\mathrm{R})}{\partial \boldsymbol{\theta}} \Delta \boldsymbol{\theta}\right) \quad \in \mathbb{R}^3
f(R⊕δθ)≈f(R)+∂θ∂f(R)Δθ≜f(R)+Exp(∂θ∂f(R)Δθ)∈R3
3.3 常用旋转雅可比矩阵 (重点)
已知向量 a \mathbf{a} a,旋转角度 θ \theta θ,旋转轴 u \mathbf{u} u,可以有三种旋转表示: θ = θ u , q = q { θ } , R = R { θ } \boldsymbol{\theta}=\theta\mathbf{u},\mathbf{q}=\mathbf{q}\{\boldsymbol{\theta}\},\mathbf{R}=\mathbf{R}\{\boldsymbol{\theta}\} θ=θu,q=q{θ},R=R{θ}。下面求三者各自对应的雅可比矩阵。
(1)对向量求雅可比
∂
(
q
⊗
a
⊗
q
∗
)
∂
a
=
∂
(
R
a
)
∂
a
=
R
\frac{\partial(\mathbf{q} \otimes \mathbf{a} \otimes \mathbf{q} *)}{\partial \mathbf{a}}=\frac{\partial(\mathbf{R} \mathbf{a})}{\partial \mathbf{a}}=\mathbf{R}
∂a∂(q⊗a⊗q∗)=∂a∂(Ra)=R
(2)对四元数求雅可比
四元数 q = w + v \mathbf{q}=w+\mathbf{v} q=w+v,则旋转后向量 a ′ \mathbf{a}^{\prime} a′为:
a ′ = q ⊗ a ⊗ q ∗ = ( w + v ) ⊗ a ⊗ ( w − v ) = w 2 a + w ( v ⊗ a − a ⊗ v ) − v ⊗ a ⊗ v = w 2 a + 2 w ( v × a ) − [ ( − v ⊤ a + v × a ) ⊗ v ] = w 2 a + 2 w ( v × a ) − [ ( − v ⊤ a ) v + ( v × a ) ⊗ v ] = w 2 a + 2 w ( v × a ) − [ ( − v ⊤ a ) v − ( v × a ) ⊤ v + ( v × a ) × v ] = w 2 a + 2 w ( v × a ) − [ ( − v ⊤ a ) v + ( v ⊤ v ) a − ( v ⊤ a ) v ] = w 2 a + 2 w ( v × a ) + 2 ( v ⊤ a ) v − ( v ⊤ v ) a \begin{aligned} \mathbf{a}^{\prime} &=\mathbf{q} \otimes \mathbf{a} \otimes \mathbf{q}^{*} \\ &=(w+\mathbf{v}) \otimes \mathbf{a} \otimes(w-\mathbf{v}) \\ &=w^{2} \mathbf{a}+w(\mathbf{v} \otimes \mathbf{a}-\mathbf{a} \otimes \mathbf{v})-\mathbf{v} \otimes \mathbf{a} \otimes \mathbf{v} \\ &=w^{2} \mathbf{a}+2 w(\mathbf{v} \times \mathbf{a})-\left[\left(-\mathbf{v}^{\top} \mathbf{a}+\mathbf{v} \times \mathbf{a}\right) \otimes \mathbf{v}\right] \\ &=w^{2} \mathbf{a}+2 w(\mathbf{v} \times \mathbf{a})-\left[\left(-\mathbf{v}^{\top} \mathbf{a}\right) \mathbf{v}+(\mathbf{v} \times \mathbf{a}) \otimes \mathbf{v}\right] \\ &=w^{2} \mathbf{a}+2 w(\mathbf{v} \times \mathbf{a})-\left[\left(-\mathbf{v}^{\top} \mathbf{a}\right) \mathbf{v}-(\mathbf{v} \times \mathbf{a})^{\top} \mathbf{v}+(\mathbf{v} \times \mathbf{a}) \times \mathbf{v}\right] \\ &=w^{2} \mathbf{a}+2 w(\mathbf{v} \times \mathbf{a})-\left[\left(-\mathbf{v}^{\top} \mathbf{a}\right) \mathbf{v}+\left(\mathbf{v}^{\top} \mathbf{v}\right) \mathbf{a}-\left(\mathbf{v}^{\top} \mathbf{a}\right) \mathbf{v}\right] \\ &=w^{2} \mathbf{a}+2 w(\mathbf{v} \times \mathbf{a})+2\left(\mathbf{v}^{\top} \mathbf{a}\right) \mathbf{v}-\left(\mathbf{v}^{\top} \mathbf{v}\right) \mathbf{a} \end{aligned} a′=q⊗a⊗q∗=(w+v)⊗a⊗(w−v)=w2a+w(v⊗a−a⊗v)−v⊗a⊗v=w2a+2w(v×a)−[(−v⊤a+v×a)⊗v]=w2a+2w(v×a)−[(−v⊤a)v+(v×a)⊗v]=w2a+2w(v×a)−[(−v⊤a)v−(v×a)⊤v+(v×a)×v]=w2a+2w(v×a)−[(−v⊤a)v+(v⊤v)a−(v⊤a)v]=w2a+2w(v×a)+2(v⊤a)v−(v⊤v)a
对实部和虚部依次求偏导为:
∂
a
′
∂
w
=
2
(
w
a
+
v
×
a
)
∂
a
′
∂
v
=
−
2
w
[
a
]
×
+
2
(
v
⊤
a
I
+
v
a
⊤
)
−
2
a
v
⊤
=
2
(
v
⊤
a
I
+
v
a
⊤
−
a
v
⊤
−
w
[
a
]
×
)
\begin{aligned} \frac{\partial \mathbf{a}^{\prime}}{\partial w} &=2(w \mathbf{a}+\mathbf{v} \times \mathbf{a}) \\ \frac{\partial \mathbf{a}^{\prime}}{\partial \mathbf{v}} &=-2 w[\mathbf{a}]_{\times}+2\left(\mathbf{v}^{\top} \mathbf{a} \mathbf{I}+\mathbf{v} \mathbf{a}^{\top}\right)-2 \mathbf{a} \mathbf{v}^{\top} \\ &=2\left(\mathbf{v}^{\top} \mathbf{a} \mathbf{I}+\mathbf{v} \mathbf{a}^{\top}-\mathbf{a} \mathbf{v}^{\top}-w[\mathbf{a}]_{\times}\right) \end{aligned}
∂w∂a′∂v∂a′=2(wa+v×a)=−2w[a]×+2(v⊤aI+va⊤)−2av⊤=2(v⊤aI+va⊤−av⊤−w[a]×)
最终雅可比矩阵为:
∂
(
q
⊗
a
⊗
q
∗
)
∂
q
=
2
[
w
a
+
v
×
a
∣
v
⊤
a
I
3
+
v
a
⊤
−
a
v
⊤
−
w
[
a
]
×
]
∈
R
3
×
4
\frac{\partial(\mathbf{q} \otimes \mathbf{a} \otimes \mathbf{q} *)}{\partial \mathbf{q}}=2\left[w \mathbf{a}+\mathbf{v} \times \mathbf{a} \mid \mathbf{v}^{\top} \mathbf{a} \mathbf{I}_{3}+\mathbf{v} \mathbf{a}^{\top}-\mathbf{a} \mathbf{v}^{\top}-w[\mathbf{a}]_{\times}\right] \in \mathbb{R}^{3 \times 4}
∂q∂(q⊗a⊗q∗)=2[wa+v×a∣v⊤aI3+va⊤−av⊤−w[a]×]∈R3×4
(3)李群的右雅可比
右雅可比的矩阵可以参考下图。
已知: Exp ( θ ) ⊕ δ ϕ = Exp ( θ + δ θ ) \operatorname{Exp}(\boldsymbol{\theta}) \oplus \delta \boldsymbol{\phi}=\operatorname{Exp}(\boldsymbol{\theta}+\delta \boldsymbol{\theta}) Exp(θ)⊕δϕ=Exp(θ+δθ),则 δ ϕ = log ( Exp ( θ ) − 1 ∘ Exp ( θ + δ θ ) ) = Exp ( θ + δ θ ) ⊖ Exp ( θ ) \delta \boldsymbol{\phi}=\log \left(\operatorname{Exp}(\boldsymbol{\theta})^{-1} \circ \operatorname{Exp}(\boldsymbol{\theta}+\delta \boldsymbol{\theta})\right)=\operatorname{Exp}(\boldsymbol{\theta}+\delta \boldsymbol{\theta}) \ominus \operatorname{Exp}(\boldsymbol{\theta}) δϕ=log(Exp(θ)−1∘Exp(θ+δθ))=Exp(θ+δθ)⊖Exp(θ)。
右雅可比矩阵计算公式为:
J
r
(
θ
)
=
lim
δ
θ
→
0
Exp
(
θ
+
δ
θ
)
⊖
Exp
(
θ
)
δ
θ
=
lim
δ
θ
→
0
log
(
Exp
(
θ
)
⊤
Exp
(
θ
+
δ
θ
)
)
δ
θ
if using
R
=
lim
δ
θ
→
0
log
(
Exp
(
θ
)
∗
⊗
Exp
(
θ
+
δ
θ
)
)
δ
θ
if using
q
.
\begin{aligned} \mathbf{J}_{r}(\boldsymbol{\theta}) &=\lim _{\delta \boldsymbol{\theta} \rightarrow 0} \frac{\operatorname{Exp}(\boldsymbol{\theta}+\delta \boldsymbol{\theta}) \ominus \operatorname{Exp}(\boldsymbol{\theta})}{\delta \boldsymbol{\theta}} \\ &=\lim _{\delta \boldsymbol{\theta} \rightarrow 0} \frac{\log \left(\operatorname{Exp}(\boldsymbol{\theta})^{\top} \operatorname{Exp}(\boldsymbol{\theta}+\delta \boldsymbol{\theta})\right)}{\delta \boldsymbol{\theta}} & \text { if using } \mathbf{R} \\ &=\lim _{\delta \boldsymbol{\theta} \rightarrow 0} \frac{\log \left(\operatorname{Exp}(\boldsymbol{\theta})^{*} \otimes \operatorname{Exp}(\boldsymbol{\theta}+\delta \boldsymbol{\theta})\right)}{\delta \boldsymbol{\theta}} & \text { if using } \mathbf{q} . \end{aligned}
Jr(θ)=δθ→0limδθExp(θ+δθ)⊖Exp(θ)=δθ→0limδθlog(Exp(θ)⊤Exp(θ+δθ))=δθ→0limδθlog(Exp(θ)∗⊗Exp(θ+δθ)) if using R if using q.
具体形式为:
J
r
(
θ
)
=
I
−
1
−
cos
∥
θ
∥
∥
θ
∥
2
[
θ
]
×
+
∥
θ
∥
−
sin
∥
θ
∥
∥
θ
∥
3
[
θ
]
×
2
J
r
−
1
(
θ
)
=
I
+
1
2
[
θ
]
×
+
(
1
∥
θ
∥
2
−
1
+
cos
∥
θ
∥
2
∥
θ
∥
sin
∥
θ
∥
)
[
θ
]
×
2
\begin{aligned} \mathbf{J}_{r}(\boldsymbol{\theta}) &=\mathbf{I}-\frac{1-\cos \|\boldsymbol{\theta}\|}{\|\boldsymbol{\theta}\|^{2}}[\boldsymbol{\theta}]_{\times}+\frac{\|\boldsymbol{\theta}\|-\sin \|\boldsymbol{\theta}\|}{\|\boldsymbol{\theta}\|^{3}}[\boldsymbol{\theta}]_{\times}^{2} \\ \mathbf{J}_{r}^{-1}(\boldsymbol{\theta}) &=\mathbf{I}+\frac{1}{2}[\boldsymbol{\theta}]_{\times}+\left(\frac{1}{\|\boldsymbol{\theta}\|^{2}}-\frac{1+\cos \|\boldsymbol{\theta}\|}{2\|\boldsymbol{\theta}\| \sin \|\boldsymbol{\theta}\|}\right)[\boldsymbol{\theta}]_{\times}^{2} \end{aligned}
Jr(θ)Jr−1(θ)=I−∥θ∥21−cos∥θ∥[θ]×+∥θ∥3∥θ∥−sin∥θ∥[θ]×2=I+21[θ]×+(∥θ∥21−2∥θ∥sin∥θ∥1+cos∥θ∥)[θ]×2
其中,右雅可比矩阵有如下性质:
Exp
(
θ
+
δ
θ
)
≈
Exp
(
θ
)
Exp
(
J
r
(
θ
)
δ
θ
)
Exp
(
θ
)
Exp
(
δ
θ
)
≈
Exp
(
θ
+
J
r
−
1
(
θ
)
δ
θ
)
log
(
Exp
(
θ
)
Exp
(
δ
θ
)
)
≈
θ
+
J
r
−
1
(
θ
)
δ
θ
\begin{aligned} \operatorname{Exp}(\boldsymbol{\theta}+\delta \boldsymbol{\theta}) & \approx \operatorname{Exp}(\boldsymbol{\theta}) \operatorname{Exp}\left(\mathbf{J}_{r}(\boldsymbol{\theta}) \delta \boldsymbol{\theta}\right) \\ \operatorname{Exp}(\boldsymbol{\theta}) \operatorname{Exp}(\delta \boldsymbol{\theta}) & \approx \operatorname{Exp}\left(\boldsymbol{\theta}+\mathbf{J}_{r}^{-1}(\boldsymbol{\theta}) \delta \boldsymbol{\theta}\right) \\ \log (\operatorname{Exp}(\boldsymbol{\theta}) \operatorname{Exp}(\delta \boldsymbol{\theta})) & \approx \boldsymbol{\theta}+\mathbf{J}_{r}^{-1}(\boldsymbol{\theta}) \delta \boldsymbol{\theta} \end{aligned}
Exp(θ+δθ)Exp(θ)Exp(δθ)log(Exp(θ)Exp(δθ))≈Exp(θ)Exp(Jr(θ)δθ)≈Exp(θ+Jr−1(θ)δθ)≈θ+Jr−1(θ)δθ
(4)对旋转向量的雅可比
∂
(
q
⊗
a
⊗
q
∗
)
∂
δ
θ
=
∂
(
R
a
)
∂
δ
θ
=
lim
δ
θ
→
0
R
{
θ
+
δ
θ
}
a
−
R
{
θ
}
a
δ
θ
=
lim
δ
θ
→
0
(
R
{
θ
}
Exp
(
J
r
(
θ
)
δ
θ
)
−
R
{
θ
}
)
a
δ
θ
=
lim
δ
θ
→
0
(
R
{
θ
}
(
I
+
[
J
r
(
θ
)
δ
θ
]
×
)
−
R
{
θ
}
)
a
δ
θ
=
lim
δ
θ
→
0
R
{
θ
}
[
J
r
(
θ
)
δ
θ
]
×
a
δ
θ
=
lim
δ
θ
→
0
−
R
{
θ
}
[
a
]
×
J
r
(
θ
)
δ
θ
δ
θ
=
−
R
{
θ
}
[
a
]
×
J
r
(
θ
)
\begin{aligned} \frac{\partial\left(\mathbf{q} \otimes \mathbf{a} \otimes \mathbf{q}^{*}\right)}{\partial \delta \boldsymbol{\theta}}=\frac{\partial(\mathbf{R} \mathbf{a})}{\partial \delta \boldsymbol{\theta}} &=\lim _{\delta \boldsymbol{\theta} \rightarrow 0} \frac{\mathbf{R}\{\boldsymbol{\theta}+\delta \boldsymbol{\theta}\} \mathbf{a}-\mathbf{R}\{\boldsymbol{\theta}\} \mathbf{a}}{\delta \boldsymbol{\theta}} \\ &=\lim _{\delta \boldsymbol{\theta} \rightarrow 0} \frac{\left(\mathbf{R}\{\boldsymbol{\theta}\} \operatorname{Exp}\left(\mathbf{J}_{r}(\boldsymbol{\theta}) \delta \boldsymbol{\theta}\right)-\mathbf{R}\{\boldsymbol{\theta}\}\right) \mathbf{a}}{\delta \boldsymbol{\theta}} \\ &=\lim _{\delta \boldsymbol{\theta} \rightarrow 0} \frac{\left(\mathbf{R}\{\boldsymbol{\theta}\}\left(\mathbf{I}+\left[\mathbf{J}_{r}(\boldsymbol{\theta}) \delta \boldsymbol{\theta}\right]_{\times}\right)-\mathbf{R}\{\boldsymbol{\theta}\}\right) \mathbf{a}}{\delta \boldsymbol{\theta}} \\ &=\lim _{\delta \boldsymbol{\theta} \rightarrow 0} \frac{\mathbf{R}\{\boldsymbol{\theta}\}\left[\mathbf{J}_{r}(\boldsymbol{\theta}) \delta \boldsymbol{\theta}\right]_{\times} \mathbf{a}}{\delta \boldsymbol{\theta}} \\ &=\lim _{\delta \boldsymbol{\theta} \rightarrow 0}-\frac{\mathbf{R}\{\boldsymbol{\theta}\}[\mathbf{a}]_{\times} \mathbf{J}_{r}(\boldsymbol{\theta}) \delta \boldsymbol{\theta}}{\delta \boldsymbol{\theta}} \\ &=-\mathbf{R}\{\boldsymbol{\theta}\}[\mathbf{a}]_{\times} \mathbf{J}_{r}(\boldsymbol{\theta}) \end{aligned}
∂δθ∂(q⊗a⊗q∗)=∂δθ∂(Ra)=δθ→0limδθR{θ+δθ}a−R{θ}a=δθ→0limδθ(R{θ}Exp(Jr(θ)δθ)−R{θ})a=δθ→0limδθ(R{θ}(I+[Jr(θ)δθ]×)−R{θ})a=δθ→0limδθR{θ}[Jr(θ)δθ]×a=δθ→0lim−δθR{θ}[a]×Jr(θ)δθ=−R{θ}[a]×Jr(θ)
3.4 四元数导数
四元数导数为:
q
˙
≜
lim
Δ
t
→
0
q
(
t
+
Δ
t
)
−
q
(
t
)
Δ
t
=
lim
Δ
t
→
0
q
⊗
Δ
q
L
−
q
Δ
t
=
lim
Δ
t
→
0
q
⊗
(
[
1
1
2
Δ
ϕ
L
]
−
[
1
0
]
)
Δ
t
=
lim
Δ
t
→
0
q
⊗
(
[
0
1
2
Δ
ϕ
L
]
Δ
t
=
1
2
q
⊗
[
0
ω
L
]
\begin{aligned} \dot{\mathbf{q}} & \triangleq \lim _{\Delta t \rightarrow 0} \frac{\mathbf{q}(t+\Delta t)-\mathbf{q}(t)}{\Delta t} \\ &=\lim _{\Delta t \rightarrow 0} \frac{\mathbf{q} \otimes \Delta \mathbf{q}_{\mathcal{L}}-\mathbf{q}}{\Delta t} \\ &=\lim _{\Delta t \rightarrow 0} \frac{\mathbf{q} \otimes\left(\left[\begin{array}{c} 1 \\ \frac{1}{2} \Delta \boldsymbol{\phi}_{\mathcal{L}} \end{array}\right]-\left[\begin{array}{l} 1 \\ \mathbf{0} \end{array}\right]\right)}{\Delta t} \\ &=\lim _{\Delta t \rightarrow 0} \frac{\mathbf{q} \otimes\left(\left[\begin{array}{c} 0 \\ \frac{1}{2} \Delta \boldsymbol{\phi}_{\mathcal{L}} \end{array}\right]\right.}{\Delta t} \\ &=\frac{1}{2} \mathbf{q} \otimes\left[\begin{array}{c} 0 \\ \boldsymbol{\omega}_{\mathcal{L}} \end{array}\right] \end{aligned}
q˙≜Δt→0limΔtq(t+Δt)−q(t)=Δt→0limΔtq⊗ΔqL−q=Δt→0limΔtq⊗([121ΔϕL]−[10])=Δt→0limΔtq⊗([021ΔϕL]=21q⊗[0ωL]
定义:
Ω
(
ω
)
≜
[
ω
]
R
=
[
0
−
ω
⊤
ω
−
[
ω
]
×
]
=
[
0
−
ω
x
−
ω
y
−
ω
z
ω
x
0
ω
z
−
ω
y
ω
y
−
ω
z
0
ω
x
ω
z
ω
y
−
ω
x
0
]
\boldsymbol{\Omega}(\boldsymbol{\omega}) \triangleq[\boldsymbol{\omega}]_{R}=\left[\begin{array}{cc} 0 & -\boldsymbol{\omega}^{\top} \\ \boldsymbol{\omega} & -[\boldsymbol{\omega}]_{\times} \end{array}\right]=\left[\begin{array}{cccc} 0 & -\omega_{x} & -\omega_{y} & -\omega_{z} \\ \omega_{x} & 0 & \omega_{z} & -\omega_{y} \\ \omega_{y} & -\omega_{z} & 0 & \omega_{x} \\ \omega_{z} & \omega_{y} & -\omega_{x} & 0 \end{array}\right]
Ω(ω)≜[ω]R=[0ω−ω⊤−[ω]×]=⎣⎢⎢⎡0ωxωyωz−ωx0−ωzωy−ωyωz0−ωx−ωz−ωyωx0⎦⎥⎥⎤
最终得到四元数和旋转矩阵的导数公式(局部坐标):
q
˙
=
1
2
Ω
(
ω
L
)
q
=
1
2
q
⊗
ω
L
,
R
˙
=
R
[
ω
L
]
×
\dot{\mathbf{q}}=\frac{1}{2} \boldsymbol{\Omega}\left(\boldsymbol{\omega}_{\mathcal{L}}\right) \mathbf{q}=\frac{1}{2} \mathbf{q} \otimes \boldsymbol{\omega}_{\mathcal{L}}, \quad \dot{\mathbf{R}}=\mathbf{R}\left[\boldsymbol{\omega}_{\mathcal{L}}\right]_{\times}
q˙=21Ω(ωL)q=21q⊗ωL,R˙=R[ωL]×
全局坐标下导数为:
q
˙
=
1
2
ω
G
⊗
q
,
R
˙
=
[
ω
G
]
×
R
\dot{\mathbf{q}}=\frac{1}{2} \boldsymbol{\omega}_{\mathcal{G}}\otimes\mathbf{q}, \quad \dot{\mathbf{R}}=[\boldsymbol{\omega}_{\mathcal{G}}]_{\times}\mathbf{R}
q˙=21ωG⊗q,R˙=[ωG]×R
4. IMU运动学方程
在文中作者列出了使用误差卡尔曼滤波的原因
:
方向误差状态极小
(和自由度有相同的参数数量),能够避免过度参数化
以及由此产生的协方差矩阵奇异
的风险。- 误差状态系统始终在接近原点的位置运行,因此
远离可能的参数奇异性,万向节锁定问题
等,从而保证了线性化有效性始终不变
。 - 误差状态量始终很小,这意味着所有
二阶项都可以忽略不计
。 这使得 Jacobian 的计算非常容易和快速。 - 误差状态量动态变化缓慢,因为所有的
large-signal
动态变化已经被包含在名义状态里了。因此,在卡尔曼滤波过程中,相较于预测过程,可以以更低的频率来执行更新过程
。
在文中作者对误差状态卡尔曼滤波还做了如下解释:
- 在误差状态卡尔曼滤波中,会讨论到三种状态量:
真值状态,名义状态和误差状态
,真值状态由名义状态和误差状态组合表示(如线性组合,四元数乘积或矩阵乘积方式等)。其思想是将名义状态视为 large-signal(以非线性方式可积分)
,将误差状态视为 small-signal(线性可积分并适用于线性高斯滤波)
。 - 在使用 IMU 进行定位时,高频 IMU 数据
u
m
\mathbf{u_m}
um 被积分到名义状态
x
\mathbf{x}
x 中。名义状态未考虑噪声
w
\mathbf{w}
w,结果它会积累误差。误差状态
δ
x
\delta \mathbf{x}
δx 通过误差状态卡尔曼滤波进行估计,
包含了所有的噪声和扰动
。误差状态由小幅度的信号组成,其更新函数由线性动态系统定义,其动态、控制和测量矩阵均由名义状态的值计算得出。误差状态卡尔曼修正是在其它传感器信息(例如GPS,视觉等)到达时执行的。
4.1 局部坐标系统连续时刻微分方程
下表是误差状态卡尔曼滤波中使用到的所有变量,需要注意的时:这里的方向误差
δ
q
\delta\mathbf{q}
δq 是定义在局部坐标
下的。后面会再给出全局坐标
下的方向计算公式。
4.1.1 真值状态微分方程
首先介绍真值状态微分方程
:
p
˙
t
=
v
t
v
˙
t
=
a
t
q
˙
t
=
1
2
q
t
⊗
ω
t
a
˙
b
t
=
a
w
ω
˙
b
t
=
ω
w
g
˙
t
=
0
\begin{aligned} \dot{\mathrm{p}}_{t} &=\mathbf{v}_{t} \\ \dot{\mathbf{v}}_{t} &=\mathbf{a}_{t} \\ \dot{\mathbf{q}}_{t} &=\frac{1}{2} \mathbf{q}_{t} \otimes \boldsymbol{\omega}_{t} \\ \dot{\mathbf{a}}_{b t} &=\mathbf{a}_{w} \\ \dot{\boldsymbol{\omega}}_{b t} &=\boldsymbol{\omega}_{w} \\ \dot{\mathbf{g}}_{t} &=0 \end{aligned}
p˙tv˙tq˙ta˙btω˙btg˙t=vt=at=21qt⊗ωt=aw=ωw=0
其中,加速度真值
a
t
\mathbf{a}_t
at,角速度真值
ω
t
\mathbf{\omega_t}
ωt 可以从 IMU 中获得,IMU测量值
a
m
\mathbf{a}_m
am 和
ω
m
\mathbf{\omega}_m
ωm 分别为:
a
m
=
R
t
⊤
(
a
t
−
g
t
)
+
a
b
t
+
a
n
ω
m
=
ω
t
+
ω
b
t
+
ω
n
\begin{aligned} \mathbf{a}_{m} &=\mathbf{R}_{t}^{\top}\left(\mathbf{a}_{t}-\mathbf{g}_{t}\right)+\mathbf{a}_{b t}+\mathbf{a}_{n} \\ \boldsymbol{\omega}_{m} &=\boldsymbol{\omega}_{t}+\boldsymbol{\omega}_{b t}+\boldsymbol{\omega}_{n} \end{aligned}
amωm=Rt⊤(at−gt)+abt+an=ωt+ωbt+ωn
对上面两个方程进行移项变换,可以得到真值
a
t
\mathbf{a}_t
at 和
ω
t
\mathbf{\omega_t}
ωt 为:
a
t
=
R
t
(
a
m
−
a
b
t
−
a
n
)
+
g
t
ω
t
=
ω
m
−
ω
b
t
−
ω
n
\begin{aligned} \mathbf{a}_{t} &=\mathbf{R}_{t}\left(\mathbf{a}_{m}-\mathbf{a}_{b t}-\mathbf{a}_{n}\right)+\mathbf{g}_{t} \\ \boldsymbol{\omega}_{t} &=\boldsymbol{\omega}_{m}-\boldsymbol{\omega}_{b t}-\boldsymbol{\omega}_{n} \end{aligned}
atωt=Rt(am−abt−an)+gt=ωm−ωbt−ωn
最终,真值状态微分方程
为:
p
˙
t
=
v
t
v
˙
t
=
R
t
(
a
m
−
a
b
t
−
a
n
)
+
g
t
q
˙
t
=
1
2
q
t
⊗
(
ω
m
−
ω
b
t
−
ω
n
)
a
˙
b
t
=
a
w
ω
˙
b
t
=
ω
w
g
˙
t
=
0
\begin{aligned} \dot{\mathbf{p}}_{t} &=\mathbf{v}_{t} \\ \dot{\mathbf{v}}_{t} &=\mathbf{R}_{t}\left(\mathbf{a}_{m}-\mathbf{a}_{b t}-\mathbf{a}_{n}\right)+\mathbf{g}_{t} \\ \dot{\mathbf{q}}_{t} &=\frac{1}{2} \mathbf{q}_{t} \otimes\left(\boldsymbol{\omega}_{m}-\boldsymbol{\omega}_{b t}-\boldsymbol{\omega}_{n}\right) \\ \dot{\mathbf{a}}_{b t} &=\mathbf{a}_{w} \\ \dot{\boldsymbol{\omega}}_{b t} &=\boldsymbol{\omega}_{w} \\ \dot{\mathbf{g}}_{t} &=0 \end{aligned}
p˙tv˙tq˙ta˙btω˙btg˙t=vt=Rt(am−abt−an)+gt=21qt⊗(ωm−ωbt−ωn)=aw=ωw=0
4.1.2 名义状态微分方程
由于名义状态不包含误差和扰动项,可以直接得到其微分方程
为:
p
˙
=
v
v
˙
=
R
(
a
m
−
a
b
)
+
g
q
˙
=
1
2
q
⊗
(
ω
m
−
ω
b
)
a
˙
b
=
0
ω
˙
b
=
0
g
˙
=
0
\begin{aligned} \dot{\mathbf{p}} &=\mathbf{v} \\ \dot{\mathbf{v}} &=\mathbf{R}\left(\mathbf{a}_{m}-\mathbf{a}_{b}\right)+\mathbf{g} \\ \dot{\mathbf{q}} &=\frac{1}{2} \mathbf{q} \otimes\left(\boldsymbol{\omega}_{m}-\boldsymbol{\omega}_{b}\right) \\ \dot{\mathbf{a}}_{b} &=0 \\ \dot{\boldsymbol{\omega}}_{b} &=0 \\ \dot{\mathrm{g}} &=0 \end{aligned}
p˙v˙q˙a˙bω˙bg˙=v=R(am−ab)+g=21q⊗(ωm−ωb)=0=0=0
4.1.3 误差状态微分方程
在忽略二阶项后,这里先给出误差状态的线性微分方程
。这里重点是速度误差状态量
δ
v
{\delta\mathbf{v}}
δv 和方向误差状态量
δ
θ
\delta\boldsymbol{\theta}
δθ 微分方程的推导。
δ
p
˙
=
δ
v
δ
v
˙
=
−
R
[
a
m
−
a
b
]
×
δ
θ
−
R
δ
a
b
+
δ
g
−
R
a
n
δ
θ
˙
=
−
[
ω
m
−
ω
b
]
×
δ
θ
−
δ
ω
b
−
ω
n
δ
a
˙
b
=
a
w
δ
ω
˙
b
=
ω
w
δ
g
˙
=
0
\begin{aligned} \dot{\delta\mathbf{p}} &=\delta \mathbf{v} \\ \dot{\delta \mathbf{v}} &=-\mathbf{R}\left[\mathbf{a}_{m}-\mathbf{a}_{b}\right]_{\times} \delta \boldsymbol{\theta}-\mathbf{R} \delta \mathbf{a}_{b}+\delta \mathbf{g}-\mathbf{R a}_{n} \\ \dot{\delta\boldsymbol{\theta}} &=-\left[\boldsymbol{\omega}_{m}-\boldsymbol{\omega}_{b}\right]_{\times} \delta \boldsymbol{\theta}-\delta \boldsymbol{\omega}_{b}-\boldsymbol{\omega}_{n} \\ \dot{\delta\mathbf{a}}_{b} &=\mathbf{a}_{w} \\ \dot{\delta\boldsymbol{\omega}}_{b} &=\boldsymbol{\omega}_{w} \\ \dot{\delta \mathbf{g}} &=0 \\ \end{aligned}
δp˙δv˙δθ˙δa˙bδω˙bδg˙=δv=−R[am−ab]×δθ−Rδab+δg−Ran=−[ωm−ωb]×δθ−δωb−ωn=aw=ωw=0
(1)速度误差状态量微分方程
在推导 速度误差状态微分方程
时先给出两个基本方程:
R
t
=
R
(
I
+
[
δ
θ
]
×
)
+
O
(
∥
δ
θ
∥
2
)
v
˙
=
R
a
B
+
g
\begin{aligned} \mathbf{R}_{t} &=\mathbf{R}\left(\mathbf{I}+[\delta \boldsymbol{\theta}]_{\times}\right)+O\left(\|\delta \boldsymbol{\theta}\|^{2}\right) \\ \dot{\mathbf{v}} &=\mathbf{R a}_{\mathcal{B}}+\mathbf{g} \end{aligned}
Rtv˙=R(I+[δθ]×)+O(∥δθ∥2)=RaB+g
由表3可知:
R
t
=
R
δ
R
\mathbf{R}_{t}=\mathbf{R} \delta \mathbf{R}
Rt=RδR,同时
δ
R
=
e
[
δ
θ
]
×
\delta \mathbf{R}=e^{[\delta \boldsymbol{\theta}]_{\times}}
δR=e[δθ]×,泰勒展开保留二阶项即可得到
R
t
\mathbf{R}_t
Rt。这里引入两个变量:
a
B
\operatorname{a}_{\mathcal{B}}
aB 是载体坐标下加速度large-signal值
,
δ
a
B
\delta \mathbf{a}_{\mathcal{B}}
δaB 为samll-signal值
。
a
B
\operatorname{a}_{\mathcal{B}}
aB 和
δ
a
B
\delta \mathbf{a}_{\mathcal{B}}
δaB方程为:
a
B
≜
a
m
−
a
b
δ
a
B
≜
−
δ
a
b
−
a
n
\begin{array}{c} \mathbf{a}_{\mathcal{B}} \triangleq \mathbf{a}_{m}-\mathbf{a}_{b} \\ \delta \mathbf{a}_{\mathcal{B}} \triangleq-\delta \mathbf{a}_{b}-\mathbf{a}_{n} \end{array}
aB≜am−abδaB≜−δab−an
代入上式,可得真值状态
a
t
\mathbf{a}_t
at 为:
a
t
=
R
t
(
a
m
−
a
b
t
−
a
n
)
+
g
t
=
R
t
(
a
m
−
a
n
−
(
a
b
+
δ
a
b
)
)
+
g
t
=
R
t
(
a
B
+
δ
a
B
)
+
g
+
δ
g
\begin{aligned} \mathbf{a}_{t} &=\mathbf{R}_{t}\left(\mathbf{a}_{m}-\mathbf{a}_{b t}-\mathbf{a}_{n}\right)+\mathbf{g}_{t} \\ &= \mathbf{R}_{t}(\mathbf{a}_{m} -\mathbf{a}_{n} - (\mathbf{a}_{b} + \delta \mathbf{a}_{\mathcal{b}})) + \mathbf{g}_{t} \\ &= \mathbf{R}_{t}\left(\mathbf{a}_{\mathcal{B}}+\delta \mathbf{a}_{\mathcal{B}}\right)+\mathbf{g}+\delta \mathbf{g}\end{aligned}
at=Rt(am−abt−an)+gt=Rt(am−an−(ab+δab))+gt=Rt(aB+δaB)+g+δg
现在用两种形式来表示
v
˙
t
\dot \mathbf{v}_t
v˙t:
v
˙
+
δ
v
˙
=
v
˙
t
=
R
t
(
a
m
−
a
b
t
−
a
n
)
+
g
t
\dot{\mathbf{v}}+\dot{\delta \mathbf{v}}=\dot{\mathbf{v}}_{t}=\mathbf{R}_{t}\left(\mathbf{a}_{m}-\mathbf{a}_{b t}-\mathbf{a}_{n}\right)+\mathbf{g}_{t}
v˙+δv˙=v˙t=Rt(am−abt−an)+gt
可得:
v
˙
+
δ
v
˙
=
v
˙
t
=
R
(
I
+
[
δ
θ
]
×
)
(
a
B
+
δ
a
B
)
+
g
+
δ
g
\dot{\mathbf{v}}+\dot{\delta \mathbf{v}}=\dot{\mathbf{v}}_{t}=\mathbf{R}\left(\mathbf{I}+[\delta \boldsymbol{\theta}]_{\times}\right)\left(\mathbf{a}_{\mathcal{B}}+\delta \mathbf{a}_{\mathcal{B}}\right)+\mathbf{g}+\delta \mathbf{g}
v˙+δv˙=v˙t=R(I+[δθ]×)(aB+δaB)+g+δg
对上式右边进行展开可得:
R
a
B
+
g
+
δ
v
˙
=
v
˙
t
=
R
a
B
+
R
δ
a
B
+
R
[
δ
θ
]
×
a
B
+
R
[
δ
θ
]
×
δ
a
B
+
g
+
δ
g
\mathbf{R a}_{\mathcal{B}}+\mathbf{g}+\dot{\delta \mathbf{v}}=\dot{\mathbf{v}}_{t}=\mathbf{R a}_{\mathcal{B}}+\mathbf{R} \delta \mathbf{a}_{\mathcal{B}}+\mathbf{R}[\delta \boldsymbol{\theta}]_{\times} \mathbf{a}_{\mathcal{B}}+\mathbf{R}[\delta \boldsymbol{\theta}]_{\times} \delta \mathbf{a}_{\mathcal{B}}+\mathbf{g}+\delta \mathbf{g}
RaB+g+δv˙=v˙t=RaB+RδaB+R[δθ]×aB+R[δθ]×δaB+g+δg
两边消除相同项
R
a
B
+
g
\mathbf{R a}_{\mathcal{B}}+\mathbf{g}
RaB+g,可得
δ
v
˙
\dot {\delta \mathbf{v}}
δv˙ 为:
δ
v
˙
=
R
(
δ
a
B
+
[
δ
θ
]
×
a
B
)
+
R
[
δ
θ
]
×
δ
a
B
+
δ
g
\dot{\delta \mathbf{v}}=\mathbf{R}\left(\delta \mathbf{a}_{\mathcal{B}}+[\delta \boldsymbol{\theta}]_{\times} \mathbf{a}_{\mathcal{B}}\right)+\mathbf{R}[\delta \boldsymbol{\theta}]_{\times} \delta \mathbf{a}_{\mathcal{B}}+\delta \mathbf{g}
δv˙=R(δaB+[δθ]×aB)+R[δθ]×δaB+δg
忽略上式中的二阶项
,且已知
[
a
]
×
b
=
−
[
b
]
×
a
[\mathbf{a}]_{\times} \mathbf{b}=-[\mathbf{b}]_{\times} \mathbf{a}
[a]×b=−[b]×a,可得:
δ
v
˙
=
R
(
δ
a
B
−
[
a
B
]
×
δ
θ
)
+
δ
g
\dot{\delta \mathbf{v}}=\mathbf{R}\left(\delta \mathbf{a}_{\mathcal{B}}-\left[\mathbf{a}_{\mathcal{B}}\right]_{\times} \delta \boldsymbol{\theta}\right)+\delta \mathbf{g}
δv˙=R(δaB−[aB]×δθ)+δg
最终可得:
δ
v
˙
=
R
(
−
[
a
m
−
a
b
]
×
δ
θ
−
δ
a
b
−
a
n
)
+
δ
g
\dot{\delta \mathbf{v}}=\mathbf{R}\left(-\left[\mathbf{a}_{m}-\mathbf{a}_{b}\right]_{\times} \delta \boldsymbol{\theta}-\delta \mathbf{a}_{b}-\mathbf{a}_{n}\right)+\delta \mathbf{g}
δv˙=R(−[am−ab]×δθ−δab−an)+δg
对上式进行整理,最终可得速度误差状态量的微分方程为
:
δ
v
˙
=
−
R
[
a
m
−
a
b
]
×
δ
θ
−
R
δ
a
b
+
δ
g
−
R
a
n
\dot{\delta \mathbf{v}}=-\mathbf{R}\left[\mathbf{a}_{m}-\mathbf{a}_{b}\right]_{\times} \delta \boldsymbol{\theta}-\mathbf{R} \delta \mathbf{a}_{b}+\delta \mathbf{g}-\mathbf{R} \mathbf{a}_{n}
δv˙=−R[am−ab]×δθ−Rδab+δg−Ran
(2)方向误差状态量微分方程
同样地在推导 方向误差状态量微分方程
时也先给出两个基本方程,下面是真值状态
和名义状态
四元数的微分方程:
q
˙
t
=
1
2
q
t
⊗
ω
t
q
˙
=
1
2
q
⊗
ω
\begin{aligned} \dot{\mathbf{q}}_{t} &=\frac{1}{2} \mathbf{q}_{t} \otimes \boldsymbol{\omega}_{t} \\ \dot{\mathbf{q}} &=\frac{1}{2} \mathbf{q} \otimes \boldsymbol{\omega} \end{aligned}
q˙tq˙=21qt⊗ωt=21q⊗ω
这里也引入角速度的 large-signal
和 small-signal
值,其表达式为:
ω
≜
ω
m
−
ω
b
δ
ω
≜
−
δ
ω
b
−
ω
n
\begin{array}{c} \boldsymbol{\omega} \triangleq \boldsymbol{\omega}_{m}-\boldsymbol{\omega}_{b} \\ \delta\boldsymbol{\omega} \triangleq -\delta\boldsymbol{\omega}_{b}-\boldsymbol{\omega}_{n} \end{array}
ω≜ωm−ωbδω≜−δωb−ωn
因此,角速度真值
可以写为:
ω
t
=
ω
+
δ
ω
\boldsymbol{\omega_{t}}=\boldsymbol{\omega}+\delta \boldsymbol{\omega}
ωt=ω+δω
也使用两种形式来表示
q
˙
t
\dot \mathbf{q}_t
q˙t:
q
⊗
˙
δ
q
=
q
˙
t
=
1
2
q
t
⊗
ω
t
\mathbf{q}\dot{\otimes} \delta \mathbf{q}=\dot{\mathbf{q}}_{t}=\frac{1}{2} \mathbf{q}_{t} \otimes \omega_{t}
q⊗˙δq=q˙t=21qt⊗ωt
上式两边同时展开可得:
q
˙
⊗
δ
q
+
q
⊗
δ
q
˙
=
q
˙
t
=
1
2
q
⊗
δ
q
⊗
ω
t
\dot{\mathbf{q}} \otimes \delta \mathbf{q}+\mathbf{q} \otimes \dot{\delta \mathbf{q}}=\dot \mathbf{q}_t=\frac{1}{2} \mathbf{q} \otimes \delta \mathbf{q} \otimes \boldsymbol{\omega}_{t}
q˙⊗δq+q⊗δq˙=q˙t=21q⊗δq⊗ωt
等式左边进一步展开可得:
1
2
q
⊗
ω
⊗
δ
q
+
q
⊗
δ
q
˙
=
q
˙
t
=
1
2
q
⊗
δ
q
⊗
ω
t
\frac{1}{2} \mathbf{q} \otimes \boldsymbol{\omega} \otimes \delta \mathbf{q}+\mathbf{q} \otimes \dot{\delta \mathbf{q}}=\dot \mathbf{q}_t=\frac{1}{2} \mathbf{q} \otimes \delta \mathbf{q} \otimes \boldsymbol{\omega}_{t}
21q⊗ω⊗δq+q⊗δq˙=q˙t=21q⊗δq⊗ωt
上式移项合并可得:
q
⊗
δ
q
˙
=
1
2
q
⊗
δ
q
⊗
ω
t
−
1
2
q
⊗
ω
⊗
δ
q
=
1
2
q
⊗
(
δ
q
⊗
ω
t
−
ω
⊗
δ
q
)
\mathbf{q} \otimes \dot{\delta \mathbf{q}}=\frac{1}{2} \mathbf{q} \otimes \delta \mathbf{q} \otimes \boldsymbol{\omega}_{t} - \frac{1}{2} \mathbf{q} \otimes \boldsymbol{\omega} \otimes \delta \mathbf{q} =\frac{1}{2} \mathbf{q} \otimes( \delta \mathbf{q} \otimes \boldsymbol{\omega}_{t} - \boldsymbol{\omega} \otimes \delta \mathbf{q})
q⊗δq˙=21q⊗δq⊗ωt−21q⊗ω⊗δq=21q⊗(δq⊗ωt−ω⊗δq)
两边消除同类项可得:
2
δ
q
˙
=
δ
q
⊗
ω
t
−
ω
⊗
δ
q
2 \dot{\delta \mathbf{q}}= \delta \mathbf{q} \otimes \boldsymbol{\omega}_{t} - \boldsymbol{\omega} \otimes \delta \mathbf{q}
2δq˙=δq⊗ωt−ω⊗δq
且已知
δ
q
=
e
δ
θ
/
2
\delta \mathbf{q}=e^{\delta \boldsymbol{\theta} / 2}
δq=eδθ/2,则
[
0
δ
˙
θ
]
=
2
δ
q
˙
\left[\begin{array}{c}0 \\\dot{\delta} \boldsymbol{\theta}\end{array}\right]=2 \dot{\delta \mathbf{q}}
[0δ˙θ]=2δq˙,代入上式可得:
[
0
δ
˙
θ
]
=
2
δ
q
˙
=
δ
q
⊗
ω
t
−
ω
⊗
δ
q
\left[\begin{array}{c}0 \\\dot{\delta} \boldsymbol{\theta}\end{array}\right] = 2 \dot{\delta \mathbf{q}} = \delta \mathbf{q} \otimes \boldsymbol{\omega}_{t} - \boldsymbol{\omega} \otimes \delta \mathbf{q}
[0δ˙θ]=2δq˙=δq⊗ωt−ω⊗δq
根据四元数矩阵乘法,可得:
[
0
δ
˙
θ
]
=
[
q
]
R
(
ω
t
)
δ
q
−
[
q
]
L
(
ω
)
δ
q
=
(
[
q
]
R
(
ω
t
)
−
[
q
]
L
(
ω
)
)
δ
q
\left[\begin{array}{c}0 \\\dot{\delta} \boldsymbol{\theta}\end{array}\right] = [\mathbf{q}]_{R}\left(\boldsymbol{\omega}_{t}\right) \delta \mathbf{q}-[\mathbf{q}]_{L}(\boldsymbol{\omega}) \delta \mathbf{q}= ( [\mathbf{q}]_{R}\left(\boldsymbol{\omega}_{t}\right)-[\mathbf{q}]_{L}(\boldsymbol{\omega}))\delta \mathbf{q}
[0δ˙θ]=[q]R(ωt)δq−[q]L(ω)δq=([q]R(ωt)−[q]L(ω))δq
将四元数转换为矩阵表示,可知:
[
q
]
R
(
ω
t
)
=
[
0
−
(
ω
t
)
⊤
(
ω
t
)
−
[
ω
t
]
×
]
[\mathbf{q}]_{R}\left(\boldsymbol{\omega}_{t}\right) = \left[\begin{array}{cc}0 & -\left(\boldsymbol{\omega}_{t}\right)^{\top} \\(\boldsymbol{\omega}_{t}) & -\left[\boldsymbol{\omega}_{t}\right]_{\times}\end{array}\right]
[q]R(ωt)=[0(ωt)−(ωt)⊤−[ωt]×]
[
q
]
L
(
ω
)
=
[
0
−
(
ω
)
⊤
(
ω
)
[
ω
]
×
]
[\mathbf{q}]_{L}\left(\boldsymbol{\omega}\right) = \left[\begin{array}{cc}0 & -\left(\boldsymbol{\omega}\right)^{\top} \\\left(\boldsymbol{\omega}\right) & \left[\boldsymbol{\omega}\right]_{\times}\end{array}\right]
[q]L(ω)=[0(ω)−(ω)⊤[ω]×]
带入上式可得:
[
0
δ
˙
θ
]
=
[
0
−
(
ω
t
−
ω
)
⊤
(
ω
t
−
ω
)
−
[
ω
t
+
ω
]
×
]
[
1
δ
θ
/
2
]
+
O
(
∥
δ
θ
∥
2
)
\left[\begin{array}{c}0 \\\dot{\delta} \boldsymbol{\theta}\end{array}\right] =\left[\begin{array}{cc}0 & -\left(\boldsymbol{\omega}_{t}-\boldsymbol{\omega}\right)^{\top} \\\left(\boldsymbol{\omega}_{t}-\boldsymbol{\omega}\right) & -\left[\boldsymbol{\omega}_{t}+\boldsymbol{\omega}\right]_{\times}\end{array}\right]\left[\begin{array}{c}1 \\\delta \boldsymbol{\theta} / 2\end{array}\right]+O\left(\|\delta \boldsymbol{\theta}\|^{2}\right)
[0δ˙θ]=[0(ωt−ω)−(ωt−ω)⊤−[ωt+ω]×][1δθ/2]+O(∥δθ∥2)
根据 ω t = ω + δ ω \omega_{t}=\boldsymbol{\omega}+\delta \boldsymbol{\omega} ωt=ω+δω,上式进一步处理为:
[ 0 δ ˙ θ ] = [ 0 − δ ω ⊤ δ ω − [ 2 ω + δ ω ] × ] [ 1 δ θ / 2 ] + O ( ∥ δ θ ∥ 2 ) \left[\begin{array}{c}0 \\\dot{\delta} \boldsymbol{\theta}\end{array}\right] =\left[\begin{array}{cc}0 & -\delta \boldsymbol{\omega}^{\top} \\\delta \boldsymbol{\omega} & -[2 \boldsymbol{\omega}+\delta \boldsymbol{\omega}]_{\times}\end{array}\right]\left[\begin{array}{c}1 \\\delta \boldsymbol{\theta} / 2\end{array}\right]+O\left(\|\delta \boldsymbol{\theta}\|^{2}\right) [0δ˙θ]=[0δω−δω⊤−[2ω+δω]×][1δθ/2]+O(∥δθ∥2)
最后可得 δ θ ˙ = δ ω − [ ω ] × δ θ − 1 2 [ δ ω ] × δ θ + O ( ∥ δ θ ∥ 2 ) \dot{\delta \boldsymbol{\theta}}=\delta \boldsymbol{\omega}-[\boldsymbol{\omega}]_{\times} \delta \boldsymbol{\theta}-\frac{1}{2}[\delta \boldsymbol{\omega}]_{\times} \delta \boldsymbol{\theta}+O\left(\|\delta \boldsymbol{\theta}\|^{2}\right) δθ˙=δω−[ω]×δθ−21[δω]×δθ+O(∥δθ∥2),忽略掉二阶项,可得: δ ˙ θ = − [ ω ] × δ θ + δ ω \dot{\delta} \boldsymbol{\theta}=-[\boldsymbol{\omega}]_{\times} \delta \boldsymbol{\theta}+\delta \boldsymbol{\omega} δ˙θ=−[ω]×δθ+δω。
最终可得方向误差状态微分方程为:
δ
˙
θ
=
−
[
ω
m
−
ω
b
]
×
δ
θ
−
δ
ω
b
−
ω
n
\dot{\delta} \boldsymbol{\theta}=-\left[\boldsymbol{\omega}_{m}-\boldsymbol{\omega}_{b}\right]_{\times} \delta \boldsymbol{\theta}-\delta \boldsymbol{\omega}_{b}-\boldsymbol{\omega}_{n}
δ˙θ=−[ωm−ωb]×δθ−δωb−ωn
4.2 全局坐标系统连续时微分方程
上一节中 方向误差状态量是在局部坐标系下定义的
,本节将给出全局坐标系
下误差状态微分方程。在全局坐标系下,方向四元数真值
为:
q
t
=
δ
q
⊗
q
\mathbf{q}_{t}=\delta \mathbf{q} \otimes \mathbf{q}
qt=δq⊗q
上式证明过程如下:设在局部坐标系下方向误差项为
δ
q
L
\delta \mathbf{q}_{L}
δqL,
q
\mathbf{q}
q 为局部坐标到全局坐标的旋转四元数,则在全局坐标下的方向误差为:
δ
q
G
=
q
⊗
δ
q
L
⊗
q
∗
\delta\mathbf{q}_{G} = \mathbf{q} \otimes\delta \mathbf{q}_L \otimes \mathbf{q}^*
δqG=q⊗δqL⊗q∗已知在局部坐标下方向误差真值为:
q
t
=
q
⊗
δ
q
L
\mathbf{q}_{t}=\mathbf{q} \otimes \delta\mathbf{q}_L
qt=q⊗δqL,代入上式可得:
q
t
=
q
⊗
q
∗
⊗
δ
q
G
⊗
q
=
δ
q
G
⊗
q
\mathbf{q}_{t}=\mathbf{q} \otimes \mathbf{q}^*\otimes\delta\mathbf{q}_G\otimes\mathbf{q}=\delta\mathbf{q}_G\otimes\mathbf{q}
qt=q⊗q∗⊗δqG⊗q=δqG⊗q,证明完毕。
与局部坐标下误差状态量微分方程相比,只有速度误差状态量和方向误差状态量微分方程
不同,下面将给出推导过程。
(1)速度误差状态量微分方程
在推导速度误差状态微分方程
时这里也给出两个基本方程:
R t = ( I + [ δ θ ] × ) R + O ( ∥ δ θ ∥ 2 ) v ˙ = R a B + g \begin{aligned} \mathbf{R}_{t} &=(\mathbf{I}+[\delta \boldsymbol{\theta}]_{\times})\mathbf{R}+O\left(\|\delta \boldsymbol{\theta}\|^{2}\right) \\ \dot{\mathbf{v}} &=\mathbf{R a}_{\mathcal{B}}+\mathbf{g} \end{aligned} Rtv˙=(I+[δθ]×)R+O(∥δθ∥2)=RaB+g
全局坐标下 R t = δ R ⋅ R \mathbf{R}_{t} = \delta \mathbf{R} \cdot \mathbf{R} Rt=δR⋅R,同时 δ R = e [ δ θ ] × \delta \mathbf{R}=e^{[\delta \boldsymbol{\theta}]_{\times}} δR=e[δθ]×。
这里同样引入两个变量:
a
B
\mathbf{a}_{\mathcal{B}}
aB 是载体坐标下加速度large-signal值
,
δ
a
B
\delta \mathbf{a}_{\mathcal{B}}
δaB为samll-signal值
。
同样的,也用两种形式来表示
v
˙
t
\dot \mathbf{v}_t
v˙t:
v
˙
+
δ
v
˙
=
v
˙
t
=
(
I
+
[
δ
θ
]
×
)
R
(
a
B
+
δ
a
B
)
+
g
+
δ
g
\dot{\mathbf{v}}+\dot{\delta \mathbf{v}}=\dot{\mathbf{v}}_{t}=\left(\mathbf{I}+[\delta \boldsymbol{\theta}]_{\times}\right) \mathbf{R}\left(\mathbf{a}_{\mathcal{B}}+\delta \mathbf{a}_{\mathcal{B}}\right)+\mathbf{g}+\delta \mathbf{g}
v˙+δv˙=v˙t=(I+[δθ]×)R(aB+δaB)+g+δg
同时对上式进行展开可得:
R
a
B
+
g
+
δ
v
˙
=
R
a
B
+
R
δ
a
B
+
[
δ
θ
]
×
R
a
B
+
[
δ
θ
]
×
R
δ
a
B
+
g
+
δ
g
\mathbf{R a}_{\mathcal{B}}+\mathbf{g}+\dot{\delta \mathbf{v}}=\mathbf{R} \mathbf{a}_{\mathcal{B}}+\mathbf{R} \delta \mathbf{a}_{\mathcal{B}}+[\delta \boldsymbol{\theta}]_{\times} \mathbf{R} \mathbf{a}_{\mathcal{B}}+[\delta \boldsymbol{\theta}]_{\times} \mathbf{R} \delta \mathbf{a}_{\mathcal{B}}+\mathbf{g}+\delta \mathbf{g}
RaB+g+δv˙=RaB+RδaB+[δθ]×RaB+[δθ]×RδaB+g+δg
两边消除相同项
R
a
B
+
g
\mathbf{R a}_{\mathcal{B}}+\mathbf{g}
RaB+g,可得
δ
v
˙
\dot {\delta \mathbf{v}}
δv˙ 为:
δ
v
˙
=
R
δ
a
B
+
[
δ
θ
]
×
R
a
B
+
[
δ
θ
]
×
R
δ
a
B
+
δ
g
\dot{\delta \mathbf{v}}=\mathbf{R}\delta \mathbf{a}_{\mathcal{B}} + [\delta \boldsymbol{\theta}]_{\times} \mathbf{R}\mathbf{a}_{\mathcal{B}}+[\delta \boldsymbol{\theta}]_{\times}\mathbf{R} \delta \mathbf{a}_{\mathcal{B}}+\delta \mathbf{g}
δv˙=RδaB+[δθ]×RaB+[δθ]×RδaB+δg
消除上式中的二阶项
,且已知
[
a
]
×
b
=
−
[
b
]
×
a
[\mathbf{a}]_{\times} \mathbf{b}=-[\mathbf{b}]_{\times} \mathbf{a}
[a]×b=−[b]×a,可得:
δ
v
˙
=
R
δ
a
B
−
[
R
a
B
]
×
δ
θ
+
δ
g
\dot{\delta \mathbf{v}}=\mathbf{R} \delta \mathbf{a}_{\mathcal{B}}-\left[\mathbf{R} \mathbf{a}_{\mathcal{B}}\right]_{\times} \delta \boldsymbol{\theta}+\delta \mathbf{g}
δv˙=RδaB−[RaB]×δθ+δg
最终可得全局坐标系下 速度误差状态量的微分方程
:
δ
v
˙
=
−
[
R
(
a
m
−
a
b
)
]
×
δ
θ
−
R
δ
a
b
+
δ
g
−
R
a
n
\dot{\delta \mathbf{v}}=-\left[\mathbf{R}\left(\mathbf{a}_{m}-\mathbf{a}_{b}\right)\right]_{\times} \delta \boldsymbol{\theta}-\mathbf{R} \delta \mathbf{a}_{b}+\delta \mathbf{g}-\mathbf{R} \mathbf{a}_{n}
δv˙=−[R(am−ab)]×δθ−Rδab+δg−Ran
(2)方向误差状态量微分方程
和局部坐标下推导一样,真值状态和名义状态四元数微分方程为:
q
˙
t
=
1
2
q
t
⊗
ω
t
q
˙
=
1
2
q
⊗
ω
\begin{aligned} \dot{\mathbf{q}}_{t} &=\frac{1}{2} \mathbf{q}_{t} \otimes \boldsymbol{\omega}_{t} \\ \dot{\mathbf{q}} &=\frac{1}{2} \mathbf{q} \otimes \boldsymbol{\omega} \end{aligned}
q˙tq˙=21qt⊗ωt=21q⊗ω
现在也使用两种形式来表示
q
˙
t
\dot \mathbf{q}_t
q˙t,但记住这里使用的是全局坐标下的方向误差
:
( δ q ⊗ q ) ˙ = q ˙ t = 1 2 q t ⊗ ω t \dot{(\delta \mathbf{q} \otimes \mathbf{q})} = \dot{\mathbf{q}}_{t}=\frac{1}{2} \mathbf{q}_{t} \otimes \omega_{t} (δq⊗q)˙=q˙t=21qt⊗ωt
上式两边同时展开可得:
δ
q
˙
⊗
q
+
δ
q
⊗
q
˙
=
1
2
δ
q
⊗
q
⊗
ω
t
\dot{\delta \mathbf{q}} \otimes \mathbf{q}+\delta \mathbf{q} \otimes \dot{\mathbf{q}}=\frac{1}{2} \delta \mathbf{q} \otimes \mathbf{q} \otimes \boldsymbol{\omega}_{t}
δq˙⊗q+δq⊗q˙=21δq⊗q⊗ωt
进一步展开可得:
δ q ˙ ⊗ q + 1 2 δ q ⊗ q ⊗ ω = 1 2 δ q ⊗ q ⊗ ω t \dot{\delta \mathbf{q}} \otimes \mathbf{q}+\frac{1}{2} \delta \mathbf{q} \otimes \mathbf{q} \otimes \boldsymbol{\omega}=\frac{1}{2} \delta \mathbf{q} \otimes \mathbf{q} \otimes \boldsymbol{\omega}_{t} δq˙⊗q+21δq⊗q⊗ω=21δq⊗q⊗ωt
根据
ω
t
=
ω
+
δ
ω
\boldsymbol{\omega}_{t}=\boldsymbol{\omega}+\delta \boldsymbol{\omega}
ωt=ω+δω,移项可得:
δ
q
˙
⊗
q
=
1
2
δ
q
⊗
q
⊗
δ
ω
\dot{\delta \mathbf{q}} \otimes \mathbf{q}=\frac{1}{2} \delta \mathbf{q} \otimes \mathbf{q} \otimes \delta \boldsymbol{\omega}
δq˙⊗q=21δq⊗q⊗δω
上式两边右乘
q
∗
\mathbf{q}^{*}
q∗,可得:
δ
q
˙
=
1
2
δ
q
⊗
q
⊗
δ
ω
⊗
q
∗
=
1
2
δ
q
⊗
(
R
δ
ω
)
=
1
2
δ
q
⊗
δ
ω
G
\begin{aligned}\dot{\delta \mathbf{q}} &=\frac{1}{2} \delta \mathbf{q} \otimes \mathbf{q} \otimes \delta \boldsymbol{\omega} \otimes \mathbf{q}^{*} \\&=\frac{1}{2} \delta \mathbf{q} \otimes(\mathbf{R} \delta \boldsymbol{\omega}) \\&=\frac{1}{2} \delta \mathbf{q} \otimes \delta \boldsymbol{\omega}_{G}\end{aligned}
δq˙=21δq⊗q⊗δω⊗q∗=21δq⊗(Rδω)=21δq⊗δωG
且已知
δ
q
=
e
δ
θ
/
2
\delta \mathbf{q}=e^{\delta \boldsymbol{\theta} / 2}
δq=eδθ/2,则
[
0
δ
˙
θ
]
=
2
δ
q
˙
\left[\begin{array}{c}0 \\\dot{\delta} \boldsymbol{\theta}\end{array}\right]=2 \dot{\delta \mathbf{q}}
[0δ˙θ]=2δq˙,代入上式可得:
[
0
δ
θ
˙
]
=
2
δ
q
˙
=
δ
q
⊗
δ
ω
G
=
Ω
(
δ
ω
G
)
δ
q
=
[
0
−
δ
ω
G
⊤
δ
ω
G
−
[
δ
ω
G
]
×
]
[
1
δ
θ
/
2
]
+
O
(
∥
δ
θ
∥
2
)
\begin{aligned}\left[\begin{array}{c}0 \\\dot{\delta \boldsymbol{\theta}}\end{array}\right] = 2 \dot{\delta\mathbf{q}}&=\delta \mathbf{q} \otimes \delta \boldsymbol{\omega}_{G} =\Omega\left(\delta \boldsymbol{\omega}_{G}\right) \delta \mathbf{q} =\left[\begin{array}{cc}0 & -\delta \boldsymbol{\omega}_{G}^{\top} \\\delta \boldsymbol{\omega}_{G}&-\left[\delta \boldsymbol{\omega}_{G}\right]_{\times}\end{array}\right]\left[\begin{array}{c}1 \\\delta \boldsymbol{\theta} / 2\end{array}\right]+O\left(\|\delta \boldsymbol{\theta}\|^{2}\right)\end{aligned}
[0δθ˙]=2δq˙=δq⊗δωG=Ω(δωG)δq=[0δωG−δωG⊤−[δωG]×][1δθ/2]+O(∥δθ∥2)
最后可得
δ
θ
˙
=
δ
ω
G
−
1
2
[
δ
ω
G
]
×
δ
θ
+
O
(
∥
δ
θ
∥
2
)
\dot{\delta \boldsymbol{\theta}}=\delta \boldsymbol{\omega}_{G}-\frac{1}{2}\left[\delta \boldsymbol{\omega}_{G}\right]_{\times} \delta \boldsymbol{\theta}+O\left(\|\delta \boldsymbol{\theta}\|^{2}\right)
δθ˙=δωG−21[δωG]×δθ+O(∥δθ∥2),忽略掉二阶项,可得:
δ
˙
θ
=
δ
ω
G
=
R
δ
ω
\dot{\delta} \boldsymbol{\theta}=\delta \boldsymbol{\omega}_{G}=\mathbf{R} \delta \boldsymbol{\omega}
δ˙θ=δωG=Rδω。最终全局坐标系下方向误差状态微分方程为:
δ
θ
˙
=
−
R
δ
ω
b
−
R
ω
n
\dot{\delta \boldsymbol{\theta}}=-\mathbf{R} \delta \boldsymbol{\omega}_{b}-\mathbf{R} \boldsymbol{\omega}_{n}
δθ˙=−Rδωb−Rωn
4.3 全局坐标系统离散时刻运动学
(1)名义状态量运动学方程
根据名义状态量微分方程,可以得到名义状态运动学方程:
p
←
p
+
v
Δ
t
+
1
2
(
R
(
a
m
−
a
b
)
+
g
)
Δ
t
2
v
←
v
+
(
R
(
a
m
−
a
b
)
+
g
)
Δ
t
q
←
q
⊗
{
(
ω
m
−
ω
b
)
Δ
t
}
a
b
←
a
b
ω
b
←
ω
b
g
←
g
\begin{aligned} {\mathbf{p}} &\leftarrow\mathbf{p} +\mathbf{v}\Delta t+\frac{1}{2}(\mathbf{R}(\mathbf{a}_m-\mathbf{a}_b)+\mathbf{g})\Delta t^2\\ {\mathbf{v}} &\leftarrow\mathbf{v} +(\mathbf{R}(\mathbf{a}_m-\mathbf{a}_b)+\mathbf{g})\Delta t \\ {\mathbf{q}} &\leftarrow\mathbf{q} \otimes \{(\boldsymbol{\omega}_m-\boldsymbol{\omega}_b)\Delta t\} \\ \mathbf{a}_b &\leftarrow \mathbf{a}_b \\ \boldsymbol{\omega}_b &\leftarrow \boldsymbol{\omega}_b \\ \mathbf{g}&\leftarrow\mathbf{g} \end{aligned}
pvqabωbg←p+vΔt+21(R(am−ab)+g)Δt2←v+(R(am−ab)+g)Δt←q⊗{(ωm−ωb)Δt}←ab←ωb←g
(2)误差状态量运动学方程
根据前面两小节的推导,误差状态量运动学方程为:
δ
p
←
p
+
δ
v
Δ
t
δ
v
←
v
+
(
−
[
R
(
a
m
−
a
b
)
]
×
δ
θ
−
R
δ
a
b
+
δ
g
)
Δ
t
+
v
i
δ
θ
←
δ
θ
−
R
δ
ω
b
Δ
t
+
θ
i
δ
a
b
←
δ
a
b
+
a
i
δ
ω
b
←
δ
ω
b
+
ω
i
δ
g
←
δ
g
\begin{aligned} \delta{\mathbf{p}} &\leftarrow\mathbf{p} +\delta\mathbf{v}\Delta t \\ \delta{\mathbf{v}} &\leftarrow\mathbf{v} +(-\left[\mathbf{R}\left(\mathbf{a}_{m}-\mathbf{a}_{b}\right)\right]_{\times} \delta \boldsymbol{\theta}-\mathbf{R} \delta \mathbf{a}_{b}+\delta \mathbf{g})\Delta t +\mathbf{v}_{\mathbf{i}}\\ \delta{\boldsymbol{\theta}} &\leftarrow\delta{\boldsymbol{\theta}} -\mathbf{R}\delta\boldsymbol{\omega}_b\Delta t+\boldsymbol{\theta}_{\mathbf{i}} \\ \delta\mathbf{a}_b &\leftarrow \delta\mathbf{a}_b +\mathbf{a}_{\mathbf{i}}\\ \delta\boldsymbol{\omega}_b &\leftarrow \delta\boldsymbol{\omega}_b+\boldsymbol{\omega}_{\mathbf{i}} \\ \delta\mathbf{g}&\leftarrow\delta\mathbf{g} \end{aligned}
δpδvδθδabδωbδg←p+δvΔt←v+(−[R(am−ab)]×δθ−Rδab+δg)Δt+vi←δθ−RδωbΔt+θi←δab+ai←δωb+ωi←δg
其中,
v
i
,
θ
i
,
a
i
,
ω
i
\mathbf{v}_{\mathbf{i}},\boldsymbol{\theta}_{\mathbf{i}},\mathbf{a}_{\mathbf{i}},\boldsymbol{\omega}_{\mathbf{i}}
vi,θi,ai,ωi 是输入高斯白噪声:
V
i
=
σ
a
~
n
2
Δ
t
2
I
[
m
2
/
s
2
]
Θ
i
=
σ
ω
~
n
2
Δ
t
2
I
[
r
a
d
2
]
A
i
=
σ
a
w
2
Δ
t
I
[
m
2
/
s
4
]
Ω
i
=
σ
ω
w
2
Δ
t
I
[
r
a
d
2
/
s
2
]
\begin{aligned} \mathbf{V}_{\mathbf{i}} &=\sigma_{\tilde{\mathbf{a}}_{n}}^{2} \Delta t^{2} \mathbf{I} &\left[m^{2} / s^{2}\right] \\ \Theta_{\mathbf{i}} &=\sigma_{\tilde{\omega}_{n}}^{2} \Delta t^{2} \mathbf{I} &\left[\mathrm{rad}^{2}\right] \\ \mathbf{A}_{\mathbf{i}} &=\sigma_{\mathbf{a}_{w}}^{2} \Delta t \mathbf{I} &\left[\mathrm{m}^{2} / \mathrm{s}^{4}\right] \\ \Omega_{\mathbf{i}} &=\sigma_{\boldsymbol{\omega}_{w}}^{2} \Delta t \mathbf{I} &\left[\mathrm{rad}^{2} / \mathrm{s}^{2}\right] \end{aligned}
ViΘiAiΩi=σa~n2Δt2I=σω~n2Δt2I=σaw2ΔtI=σωw2ΔtI[m2/s2][rad2][m2/s4][rad2/s2]
5. IMU/GPS融合定位
5.1 融合定位预测
现在来看一个实例,即IMU/GPS融合定位。关于IMU/GPS融合定位的工作原理可以查看博文:动手学无人驾驶(6):基于IMU和GPS数据融合的自车定位。
卡尔曼滤波中,首先给出系统的名义状态变量
x
\mathbf{x}
x,误差状态量
δ
x
\delta\mathbf{x}
δx,输入为
u
m
\mathbf{u}_m
um,噪声为
i
\mathbf{i}
i:
x
=
[
p
v
q
a
b
ω
b
g
]
,
δ
x
=
[
δ
p
δ
v
δ
θ
δ
a
b
δ
ω
b
δ
g
]
,
u
m
=
[
a
m
ω
m
]
,
i
=
[
v
i
θ
i
a
i
ω
i
]
\mathbf{x}=\left[\begin{array}{c} \mathbf{p} \\ \mathbf{v} \\ \mathbf{q} \\ \mathbf{a}_{b} \\ \boldsymbol{\omega}_{b} \\ \mathbf{g} \end{array}\right], \quad \delta \mathbf{x}=\left[\begin{array}{c} \delta \mathbf{p} \\ \delta \mathbf{v} \\ \delta \boldsymbol{\theta} \\ \delta \mathbf{a}_{b} \\ \delta \boldsymbol{\omega}_{b} \\ \delta \mathbf{g} \end{array}\right], \quad \mathbf{u}_{m}=\left[\begin{array}{c} \mathbf{a}_{m} \\ \boldsymbol{\omega}_{m} \end{array}\right], \quad \mathbf{i}=\left[\begin{array}{c} \mathbf{v}_{\mathbf{i}} \\ \boldsymbol{\theta}_{\mathbf{i}} \\ \mathbf{a}_{\mathbf{i}} \\ \boldsymbol{\omega}_{\mathbf{i}} \end{array}\right]
x=⎣⎢⎢⎢⎢⎢⎢⎡pvqabωbg⎦⎥⎥⎥⎥⎥⎥⎤,δx=⎣⎢⎢⎢⎢⎢⎢⎡δpδvδθδabδωbδg⎦⎥⎥⎥⎥⎥⎥⎤,um=[amωm],i=⎣⎢⎢⎡viθiaiωi⎦⎥⎥⎤
则ESKF 预测方程
可以写为:
δ
^
x
←
F
x
(
x
,
u
m
)
⋅
δ
^
x
P
←
F
x
P
F
x
⊤
+
F
i
Q
i
F
i
⊤
\begin{aligned} \hat{\delta} \mathbf{x} & \leftarrow \mathbf{F}_{\mathbf{x}}\left(\mathbf{x}, \mathbf{u}_{m}\right) \cdot \hat{\delta} \mathbf{x} \\ \mathbf{P} & \leftarrow \mathbf{F}_{\mathbf{x}} \mathbf{P} \mathbf{F}_{\mathbf{x}}^{\top}+\mathbf{F}_{\mathbf{i}} \mathbf{Q}_{\mathbf{i}} \mathbf{F}_{\mathbf{i}}^{\top} \end{aligned}
δ^xP←Fx(x,um)⋅δ^x←FxPFx⊤+FiQiFi⊤
其中,状态转移矩阵
F
x
\mathbf{F}_{\mathbf{x}}
Fx为:
F
x
=
[
I
I
Δ
t
0
0
0
0
0
I
−
[
R
(
a
m
−
a
b
)
]
×
Δ
t
−
R
Δ
t
0
I
Δ
t
0
0
I
0
−
R
Δ
t
0
0
0
0
I
0
0
0
0
0
0
I
0
0
0
0
0
0
I
]
.
\mathbf{F}_{\mathbf{x}}=\left[\begin{array}{cccccc} \mathbf{I} & \mathbf{I} \Delta t & 0 & 0 & 0 & 0 \\ 0 & \mathbf{I} & -\left[\mathbf{R}\left(\mathbf{a}_{m}-\mathbf{a}_{b}\right)\right]_{\times} \Delta t & -\mathbf{R} \Delta t & 0 & \mathbf{I} \Delta t \\ 0 & 0 & \mathbf{I} & 0 & -\mathbf{R} \Delta t & 0 \\ 0 & 0 & 0 & \mathbf{I} & 0 & 0 \\ 0 & 0 & 0 & 0 & \mathbf{I} & 0 \\ 0 & 0 & 0 & 0 & 0 & \mathbf{I} \end{array}\right] .
Fx=⎣⎢⎢⎢⎢⎢⎢⎡I00000IΔtI00000−[R(am−ab)]×ΔtI0000−RΔt0I0000−RΔt0I00IΔt000I⎦⎥⎥⎥⎥⎥⎥⎤.
系统输入雅可比矩阵
F
i
\mathbf{F}_{\mathbf{i}}
Fi 为:
F
i
=
[
0
0
0
0
I
0
0
0
0
I
0
0
0
0
I
0
0
0
0
I
0
0
0
0
]
\mathbf{F}_{\mathbf{i}}=\left[\begin{array}{llll} 0 & 0 & 0 & 0 \\ \mathbf{I} & 0 & 0 & 0 \\ 0 & \mathbf{I} & 0 & 0 \\ 0 & 0 & \mathbf{I} & 0 \\ 0 & 0 & 0 & \mathbf{I} \\ 0 & 0 & 0 & 0 \end{array}\right]
Fi=⎣⎢⎢⎢⎢⎢⎢⎡0I000000I000000I000000I0⎦⎥⎥⎥⎥⎥⎥⎤
状态转移协方差矩阵
Q
i
\mathbf{Q}_{\mathbf{i}}
Qi 为:
Q
i
=
[
V
i
0
0
0
0
Θ
i
0
0
0
0
A
i
0
0
0
0
Ω
i
]
\mathbf{Q}_{\mathbf{i}}=\left[\begin{array}{cccc} \mathbf{V}_{\mathbf{i}} & 0 & 0 & 0 \\ 0 & \Theta_{\mathbf{i}} & 0 & 0 \\ 0 & 0 & \mathbf{A}_{\mathbf{i}} & 0 \\ 0 & 0 & 0 & \boldsymbol{\Omega}_{\mathbf{i}} \end{array}\right]
Qi=⎣⎢⎢⎡Vi0000Θi0000Ai0000Ωi⎦⎥⎥⎤
需要注意的是
,由于
δ
x
\delta\mathbf{x}
δx 初始化为0
,其在预测阶段一直为0
,因此使用名义状态量作为真值状态量
。
5.2 融合定位更新
回顾误差卡尔曼滤波更新步骤:
K
=
P
H
⊤
(
H
P
H
⊤
+
V
)
−
1
δ
x
^
←
K
(
y
−
h
(
x
^
t
)
)
x
t
^
←
x
⊕
δ
x
^
P
←
(
I
−
K
H
)
P
(
I
−
K
H
)
⊤
+
K
V
K
⊤
\begin{aligned} \mathbf{K} &=\mathbf{P} \mathbf{H}^{\top}\left(\mathbf{H P H}^{\top}+\mathbf{V}\right)^{-1} \\ \hat{\delta \mathbf{x}} & \leftarrow \mathbf{K}\left(\mathbf{y}-h\left(\hat{\mathbf{x}}_{t}\right)\right) \\ \hat{\mathbf{x}_t} & \leftarrow \mathbf{x} \oplus\hat{\delta\mathbf{x}}\\ \mathbf{P} & \leftarrow(\mathbf{I}-\mathbf{K H}) \mathbf{P} (\mathbf{I}-\mathbf{K H})^{\top}+\mathbf{K}\mathbf{V}\mathbf{K}^{\top} \end{aligned}
Kδx^xt^P=PH⊤(HPH⊤+V)−1←K(y−h(x^t))←x⊕δx^←(I−KH)P(I−KH)⊤+KVK⊤
其中, y \mathbf{y} y 是测量值, H \mathbf{H} H 是对误差状态量 δ x \delta\mathbf{x} δx 的测量矩阵,由于误差状态量为0,使用名义状态量作为真值状态量。
测量矩阵
H
\mathbf{H}
H 为:
H
≜
∂
h
∂
δ
x
∣
x
=
∂
h
∂
x
t
∣
x
∂
x
t
∂
δ
x
∣
x
=
H
x
X
δ
x
\left.\mathbf{H} \triangleq \frac{\partial h}{\partial \delta \mathbf{x}}\right|_{\mathbf{x}}=\left.\left.\frac{\partial h}{\partial \mathbf{x}_{t}}\right|_{\mathbf{x}} \frac{\partial \mathbf{x}_{t}}{\partial \delta \mathbf{x}}\right|_{\mathbf{x}}=\mathbf{H}_{\mathbf{x}} \mathbf{X}_{\delta \mathbf{x}}
H≜∂δx∂h∣∣∣∣x=∂xt∂h∣∣∣∣x∂δx∂xt∣∣∣∣x=HxXδx
矩阵
H
x
\mathbf{H}_\mathbf{x}
Hx为 测量函数对真值状态量的雅可比矩阵,根据使用传感器的不同,其形式也不同。真值状态对误差状态矩阵
X
δ
x
\mathbf{X}_{\delta\mathbf{x}}
Xδx 为:
X
δ
x
=
[
∂
(
p
+
δ
p
)
∂
δ
p
∂
(
v
+
δ
v
)
∂
δ
v
0
∂
(
q
⊗
δ
q
)
∂
δ
θ
∂
(
a
b
+
δ
a
b
)
∂
δ
a
b
0
∂
(
ω
b
+
δ
ω
b
)
∂
δ
ω
b
∂
(
g
+
δ
g
)
∂
δ
g
]
\mathbf{X}_{\delta \mathbf{x}}=\left[\begin{array}{cccc}\frac{\partial(\mathbf{p}+\delta \mathbf{p})}{\partial \delta \mathbf{p}} & & & & & \\& \frac{\partial(\mathbf{v}+\delta \mathbf{v})}{\partial \delta \mathbf{v}} & & & 0& \\& & \frac{\partial(\mathbf{q} \otimes \delta \mathbf{q})}{\partial \delta \boldsymbol{\theta}} & & & \\& & & \frac{\partial\left(\mathbf{a}_{b}+\delta \mathbf{a}_{b}\right)}{\partial \delta \mathbf{a}_{b}} & & \\& 0 & && \frac{\partial\left(\boldsymbol{\omega}_{b}+\delta \boldsymbol{\omega}_{b}\right)}{\partial \delta \boldsymbol{\omega}_{b}} & \\& & & & &\frac{\partial(\mathbf{g}+\delta \mathbf{g})}{\partial \delta \mathbf{g}}\end{array}\right]
Xδx=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡∂δp∂(p+δp)∂δv∂(v+δv)0∂δθ∂(q⊗δq)∂δab∂(ab+δab)0∂δωb∂(ωb+δωb)∂δg∂(g+δg)⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤
进一步简化,可得:
X
δ
x
=
[
I
6
0
0
0
Q
δ
θ
0
0
0
I
9
]
\mathbf{X}_{\delta \mathbf{x}}=\left[\begin{array}{cccc} \mathbf{I}_6 & 0& 0 \\0 &\mathbf{Q}_{\delta\boldsymbol\theta} &0 \\0&0 & \mathbf{I}_9 \end{array}\right]
Xδx=⎣⎡I6000Qδθ000I9⎦⎤
其中,
Q
δ
θ
=
1
2
[
−
q
x
−
q
y
−
q
z
q
w
q
z
−
q
y
−
q
z
q
w
q
z
q
y
−
q
x
q
w
]
\mathbf{Q}_{\delta\boldsymbol\theta}=\frac{1}{2}\left[\begin{array}{cccc}-q_x &-q_y & -q_z \\ q_w &q_z&-q_y \\-q_z &q_w &q_z \\ q_y &-q_x & q_w\end{array}\right]
Qδθ=21⎣⎢⎢⎡−qxqw−qzqy−qyqzqw−qx−qz−qyqzqw⎦⎥⎥⎤
接下来是更新状态量:
p
←
p
+
δ
p
v
←
v
+
δ
v
q
←
q
{
δ
θ
^
}
⊗
q
a
b
←
a
b
+
δ
a
b
ω
b
←
ω
b
+
δ
ω
b
g
←
g
+
δ
g
\begin{aligned} {\mathbf{p}} &\leftarrow\mathbf{p} + \delta\mathbf{p}\\ {\mathbf{v}} &\leftarrow\mathbf{v} + \delta\mathbf{v} \\ {\mathbf{q}} &\leftarrow\mathbf{q}\{\hat{\delta\boldsymbol{\theta}}\} \otimes \mathbf{q} \\ \mathbf{a}_b &\leftarrow \mathbf{a}_b + \delta\mathbf{a}_b \\ \boldsymbol{\omega}_b &\leftarrow \boldsymbol{\omega}_b +\delta\boldsymbol{\omega}_b\\ \mathbf{g}&\leftarrow\mathbf{g} + \delta\mathbf{g} \end{aligned}
pvqabωbg←p+δp←v+δv←q{δθ^}⊗q←ab+δab←ωb+δωb←g+δg
在误差卡尔曼滤波中,还需要对协方差矩阵
P
\mathbf{P}
P 进行重置操作:
P
←
G
P
G
⊤
\mathbf{P}\leftarrow\mathbf{G}\mathbf{P}\mathbf{G}^{\top}
P←GPG⊤
其中,矩阵
G
\mathbf{G}
G:
G
=
[
I
6
0
0
0
I
+
[
1
2
δ
θ
^
]
×
0
0
0
I
9
]
\mathbf{G}=\left[\begin{array}{cccc} \mathbf{I}_6 & 0& 0 \\0 &\mathbf{I}+[\hat{\frac{1}{2}\delta\boldsymbol{\theta}}]_{\times} &0 \\0&0 & \mathbf{I}_9 \end{array}\right]
G=⎣⎡I6000I+[21δθ^]×000I9⎦⎤
局部坐标和全局坐标方向误差总结如下表4所示:
至此,本文介绍完毕:)。