题目:流形上的运动学 - Notes for “Kalman Filter on Differentiable Manifolds” - IKFoM (I)
文章目录
- 前言
- I. About eq (12) in "Example 2: Attitude kinematics in a global reference frame"
- II. About eq (13) in "Example 3: Attitude kinematics in body frame"
- 题外话一: right- ⊕ \oplus ⊕ 和 left- ⊕ \oplus ⊕
- III. About eq (15) in "Example 5: Vectors of known magnitude in body frame"
- IV. About eqs (16)~(17) in "Example 6: Bearing-distance parameterization of visual landmarks"
- 题外话二: S 2 \mathbb{S}^2 S2 is not a Lie group 三维空间的球面 (2-Sphere) 不是李群
- 总结
- 参考文献
前言
Dongjiao He, Wei Xu, Fu Zhang 的论文 “Kalman Filter on Differentiable Manifolds”[1] 中构建的 IKFoM (“流形上的迭代 Kalman 滤波”) 的系统架构比较复杂. 另外, 为了通用化的缘故, 同时在欧几里得空间、 S O ( 3 ) SO(3) SO(3)、 S 2 \mathbb{S}^{2} S2 上推导, 增加了论文阅读的难度. 尝试在论文原作者简练的推导过程上, 以自己的理解进行详细和啰嗦地推导, 作为自己学习理解的记录.
本文是涉及到论文 “Kalman Filter on Differentiable Manifolds”[1] 的 “IV. CANONICAL REPRESENTATION OF ON-MANIFOLD SYSTEMS” 部分的阅读笔记. 主要关于
- S O ( 3 ) SO(3) SO(3) 上的姿态运动学
- S 2 \mathbb{S}^2 S2 上重力向量运动学
- S 2 \mathbb{S}^2 S2 上方位向量运动学
本文是相关主题的第一篇, 后一篇请参考
流形上的复合函数偏导数 - Notes for “Kalman Filter on Differentiable Manifolds“ - IKFoM (II)
I. About eq (12) in “Example 2: Attitude kinematics in a global reference frame”
下面是关于流形上的姿态运动学. 这里对符号进行了简单替换, 以使得表示更直观和符合习惯.
1. 姿态矩阵的导数 Time Derivative of an Attitude Matrix
文献 [2] 中对姿态矩阵的导数有进行推导, 是本小节的参考来源.
运动物体相对于全局坐标系 {G} 的姿态 —— B G R {_B^{G}\mathbf{R}} BGR (原文 x \mathbf{x} x) —— 运动物体的姿态在全局坐标系下的描述
运动物体在全局坐标系下的角速度 —— G ω ^G\boldsymbol{\omega} Gω —— 在采样期间 Δ t \Delta t Δt 内保持不变
假设运动物体上任一固定点
P
P
P 在运动物体局部坐标系 {B} 上的坐标为
B
P
^B\mathbf{P}
BP, 则在全局坐标系下的坐标为
G
P
=
B
G
R
B
P
(1-1)
{^G\mathbf{P}} = {_B^G\mathbf{R}}{^B\mathbf{P}} \tag{1-1}
GP=BGRBP(1-1)
因为
B
P
^B\mathbf{P}
BP 在局部坐标系 {B} 上固定, 故
B
P
˙
=
0
(1-2)
{^B\dot{\mathbf{P}}}=\mathbf{0} \tag{1-2}
BP˙=0(1-2)
对式 (1-1) 两边求导得到
G
P
˙
=
B
G
R
˙
B
P
+
B
G
R
B
P
˙
=
B
G
R
˙
B
P
(1-3)
{^G\dot{\mathbf{P}}} = {_B^G\dot{\mathbf{R}}}{^B\mathbf{P}} + {_B^G{\mathbf{R}}}{^B\dot{\mathbf{P}}} = {_B^G\dot{\mathbf{R}}}{^B\mathbf{P}} \tag{1-3}
GP˙=BGR˙BP+BGRBP˙=BGR˙BP(1-3)
又由全局坐标系下角速度和线速度的关系
G
P
˙
=
G
ω
×
G
P
=
[
G
ω
]
×
G
P
=
(
1
−
1
)
[
G
ω
]
×
B
G
R
B
P
(1-4)
{^G\dot{\mathbf{P}}} = {^G\boldsymbol{\omega}} \times {^G\mathbf{P}} = [{^G\boldsymbol{\omega}}]_\times {^G\mathbf{P}} \overset{\rm{(1-1)}}= [{^G\boldsymbol{\omega}}]_\times {_B^G\mathbf{R}}{^B\mathbf{P}} \tag{1-4}
GP˙=Gω×GP=[Gω]×GP=(1−1)[Gω]×BGRBP(1-4)
结合式 (I-3) 和 (I-4) 得
B
G
R
˙
B
P
=
[
G
ω
]
×
B
G
R
B
P
(1-5)
{_B^G\dot{\mathbf{R}}}{^B\mathbf{P}} = [{^G\boldsymbol{\omega}}]_\times {_B^G\mathbf{R}}{^B\mathbf{P}} \tag{1-5}
BGR˙BP=[Gω]×BGRBP(1-5)
因为
B
P
^{B}\mathbf{P}
BP 是任意的, 故要使式 (I-5) 成立需要满足
B
G
R
˙
=
[
G
ω
]
×
B
G
R
(1-6)
{_B^G\dot{\mathbf{R}}} = [{^G\boldsymbol{\omega}}]_\times {_B^G\mathbf{R}} \tag{1-6}
BGR˙=[Gω]×BGR(1-6)
即为原文中
x
˙
=
⌊
G
ω
⌋
⋅
x
\dot{\mathbf{x}} = \lfloor^G\boldsymbol{\omega}\rfloor\cdot \mathbf{x}
x˙=⌊Gω⌋⋅x
2. 离散化 Zero-Order Hold Discretization
由式 (1-6) 离散化得到
B
G
R
k
+
1
−
B
G
R
k
Δ
t
=
[
G
ω
k
]
×
B
G
R
k
⇒
B
G
R
k
+
1
=
B
G
R
k
+
Δ
t
[
G
ω
k
]
×
B
G
R
k
⇒
B
G
R
k
+
1
=
(
I
+
Δ
t
[
G
ω
k
]
×
)
B
G
R
k
(1-7)
\begin{aligned} \frac{{_B^G{\mathbf{R}}_{k+1}} - {_B^G{\mathbf{R}}_{k}}}{\Delta t} & = [{^G\boldsymbol{\omega}}_k]_\times {_B^G\mathbf{R}}_{k}\\ \Rightarrow \qquad {_B^G{\mathbf{R}}_{k+1}} & = {_B^G{\mathbf{R}}_{k}} + {\Delta t} [{^G\boldsymbol{\omega}}_k]_\times {_B^G\mathbf{R}}_{k}\\ \Rightarrow \qquad {_B^G{\mathbf{R}}_{k+1}} & =\left( \mathbf{I} + {\Delta t} [{^G\boldsymbol{\omega}}_k]_\times \right){_B^G\mathbf{R}}_{k}\\ \end{aligned} \tag{1-7}
ΔtBGRk+1−BGRk⇒BGRk+1⇒BGRk+1=[Gωk]×BGRk=BGRk+Δt[Gωk]×BGRk=(I+Δt[Gωk]×)BGRk(1-7)
由
E
x
p
(
⋅
)
{\rm{Exp}}(\cdot)
Exp(⋅) 的泰勒展开
E
x
p
(
Δ
t
G
ω
k
)
=
e
x
p
(
Δ
t
[
G
ω
k
]
×
)
=
∑
k
=
0
∞
Δ
t
k
k
!
[
G
ω
k
]
×
k
(1-8)
{\rm{Exp}}(\Delta t {^G\boldsymbol{\omega}_k}) = {\rm{exp}}(\Delta t [{^G\boldsymbol{\omega}_k}]_{\times}) = \sum_{k=0}^{\infty} \frac{{\Delta t}^{k}}{k!}[{^G\boldsymbol{\omega}_k}]_{\times}^k \tag{1-8}
Exp(ΔtGωk)=exp(Δt[Gωk]×)=k=0∑∞k!Δtk[Gωk]×k(1-8)
取一阶泰勒展开为
E
x
p
(
Δ
t
G
ω
k
)
≈
I
+
Δ
t
[
G
ω
k
]
×
(1-9)
{\rm{Exp}}(\Delta t {^G\boldsymbol{\omega}_k}) \approx \mathbf{I} + \Delta t [{^G\boldsymbol{\omega}_k}]_{\times} \tag{1-9}
Exp(ΔtGωk)≈I+Δt[Gωk]×(1-9)
则式 (1-7) 继续变换为
B
G
R
k
+
1
=
E
x
p
(
Δ
t
G
ω
k
)
B
G
R
k
(1-10)
{_B^G{\mathbf{R}}_{k+1}} ={\rm{Exp}}(\Delta t \,{^G\boldsymbol{\omega}_k})\,{_B^G\mathbf{R}}_{k} \tag{1-10}
BGRk+1=Exp(ΔtGωk)BGRk(1-10)
即原文中
x
k
+
1
=
E
x
p
(
Δ
t
G
ω
k
)
⋅
x
k
{\mathbf{x}}_{k+1} = {\rm{Exp}}(\Delta t \,{^G\boldsymbol{\omega}_k})\cdot{\mathbf{x}}_{k}
xk+1=Exp(ΔtGωk)⋅xk
3. 指数映射的变换 Transformation of a Exponential Map
由 “机器人中常用矩阵等式-I (汇总)” 中 Identity 11[3], 可知
[
R
k
T
G
ω
k
]
×
=
R
k
T
[
G
ω
k
]
×
R
k
(1-11)
[\mathbf{R}_{k}^{\rm\small T} \, {^{G}{\boldsymbol{\omega}}_k}]_{\times} = \mathbf{R}_{k}^{\rm\small T}\, [{^{G}{\boldsymbol{\omega}}_k}]_{\times}\, \mathbf{R}_{k} \tag{1-11}
[RkTGωk]×=RkT[Gωk]×Rk(1-11)
由旋转矩阵的正交特性可知
[
R
k
T
G
ω
k
]
×
m
=
R
k
T
[
G
ω
k
]
×
m
R
k
(
m
=
0
,
1
,
2
,
⋯
)
(1-12)
[\mathbf{R}_{k}^{\rm\small T} \, {^{G}{\boldsymbol{\omega}}_k}]_{\times}^{m} = \mathbf{R}_{k}^{\rm\small T}\, [{^{G}{\boldsymbol{\omega}}_k}]_{\times}^{m}\, \mathbf{R}_{k} \qquad (m = 0, 1, 2, \cdots) \tag{1-12}
[RkTGωk]×m=RkT[Gωk]×mRk(m=0,1,2,⋯)(1-12)
则指数映射有
E
x
p
(
Δ
t
(
R
k
T
G
ω
k
)
)
=
e
x
p
(
Δ
t
[
R
k
T
G
ω
k
]
×
)
=
∑
k
=
0
∞
Δ
t
k
k
!
[
R
k
T
G
ω
k
]
×
k
=
R
k
T
(
∑
k
=
0
∞
Δ
t
k
k
!
[
G
ω
k
]
×
k
)
R
k
=
R
k
T
E
x
p
(
Δ
t
G
ω
k
)
R
k
(1-13)
\begin{aligned} {\rm Exp}\left( \Delta t\,(\mathbf{R}_{k}^{\rm\small T}\, {^{G}{\boldsymbol{\omega}}_k})\right) &= {\rm exp}\left( \Delta t\,[\mathbf{R}_{k}^{\rm\small T}\, {^{G}{\boldsymbol{\omega}}_k}]_{\times}\right)\\ &= \sum_{k=0}^{\infty} \frac{{\Delta t}^{k}}{k!}[\mathbf{R}_{k}^{\rm\small T}\, {^G\boldsymbol{\omega}_k}]_{\times}^k\\ &= \mathbf{R}_{k}^{\rm\small T}\left( \sum_{k=0}^{\infty} \frac{{\Delta t}^{k}}{k!}[ {^G\boldsymbol{\omega}_k}]_{\times}^k\right) \mathbf{R}_{k}\\ &= \mathbf{R}_{k}^{\rm\small T}\, {\rm Exp}(\Delta t {^G\boldsymbol{\omega}_k})\, \mathbf{R}_{k} \end{aligned} \tag{1-13}
Exp(Δt(RkTGωk))=exp(Δt[RkTGωk]×)=k=0∑∞k!Δtk[RkTGωk]×k=RkT(k=0∑∞k!Δtk[Gωk]×k)Rk=RkTExp(ΔtGωk)Rk(1-13)
两边同时左乘以
R
k
R_{k}
Rk 得到
R
k
E
x
p
(
Δ
t
(
R
k
T
G
ω
k
)
)
=
E
x
p
(
Δ
t
G
ω
k
)
R
k
(1-14)
\mathbf{R}_{k}\, {\rm Exp}\left(\Delta t\,(\mathbf{R}_{k}^{\rm\small T}\, {^{G}{\boldsymbol{\omega}}_k})\right) = {\rm Exp}(\Delta t {^G\boldsymbol{\omega}_k})\, \mathbf{R}_{k} \tag{1-14}
RkExp(Δt(RkTGωk))=Exp(ΔtGωk)Rk(1-14)
即原文中
x
k
+
1
=
E
x
p
(
Δ
t
G
ω
k
)
⋅
x
k
=
x
k
E
x
p
(
Δ
t
(
x
k
T
⋅
G
ω
k
)
)
\mathbf{x}_{k+1}={\rm Exp}(\Delta t {^G\boldsymbol{\omega}_k})\cdot \mathbf{x}_{k} = \mathbf{x}_{k}\, {\rm Exp}\left( \Delta t\,(\mathbf{x}_{k}^{\rm\small T}\cdot {^{G}{\boldsymbol{\omega}}_k})\right)
xk+1=Exp(ΔtGωk)⋅xk=xkExp(Δt(xkT⋅Gωk))
这样完整得到了原文中式 (12).
II. About eq (13) in “Example 3: Attitude kinematics in body frame”
物体相对于全局坐标系运动的角速度在全局坐标系 {G} 中的描述记为 G ω ^G\boldsymbol{\omega} Gω;
物体相对于全局坐标系运动的角速度在运动物体局部坐标系 {B} 中的描述记为 B ω ^B\boldsymbol{\omega} Bω.
两者相互之间的转换关系为
G
ω
=
B
G
R
B
ω
(2-1)
{^G\boldsymbol{\omega}} = {_B^{G}\mathbf{R}}\, {^B\boldsymbol{\omega}} \tag{2-1}
Gω=BGRBω(2-1)
1. 姿态矩阵的导数 Time Derivative of an Attitude Matrix
将式 (2-1) 代入式 (1-6) , 再结合式 (1-11) 得到
B
G
R
˙
=
[
G
ω
]
×
B
G
R
(
2
−
1
)
=
[
B
G
R
B
ω
]
×
B
G
R
(
1
−
11
)
=
B
G
R
[
B
ω
]
×
B
G
R
T
B
G
R
=
B
G
R
[
B
ω
]
×
(2-2)
\begin{aligned} {_B^G\dot{\mathbf{R}}} & = [{^G\boldsymbol{\omega}}]_\times\, {_B^G\mathbf{R}} \\ {\small (2-1)} \qquad & = [{_B^{G}\mathbf{R}}\, {^B\boldsymbol{\omega}}]_{\times}\, {_B^G\mathbf{R}}\\ {\small (1-11)} \qquad& = {_B^{G}\mathbf{R}}\, [{^B\boldsymbol{\omega}}]_{\times} \, {_B^{G}\mathbf{R}}^{\rm \small T}\, {_B^G\mathbf{R}} \\ &= {_B^{G}\mathbf{R}}\, [{^B\boldsymbol{\omega}}]_{\times} \end{aligned}\tag{2-2}
BGR˙(2−1)(1−11)=[Gω]×BGR=[BGRBω]×BGR=BGR[Bω]×BGRTBGR=BGR[Bω]×(2-2)
即原文中的
x
˙
=
x
⋅
⌊
B
ω
⌋
\dot{\mathbf{x}} = \mathbf{x}\cdot \lfloor {^B\boldsymbol{\omega}}\rfloor
x˙=x⋅⌊Bω⌋
2. 离散化 Zero-Order Hold Discretization
由式 (2-2) 做离散化可知
B
G
R
k
+
1
−
B
G
R
k
Δ
t
=
B
G
R
k
[
B
ω
k
]
×
⇒
B
G
R
k
+
1
=
B
G
R
k
+
B
G
R
k
Δ
t
[
B
ω
k
]
×
⇒
B
G
R
k
+
1
=
B
G
R
k
(
I
+
Δ
t
[
B
ω
k
]
×
)
(
1st-Order Taylor Approximation
)
⇒
B
G
R
k
+
1
=
B
G
R
k
e
x
p
(
Δ
t
[
B
ω
k
]
×
)
⇒
B
G
R
k
+
1
=
B
G
R
k
E
x
p
(
Δ
t
B
ω
k
)
(2-3)
\begin{aligned} \frac{{_B^G{\mathbf{R}}_{k+1}} -{_B^G{\mathbf{R}}_{k}}} {\Delta t} & = {_B^{G}\mathbf{R}_k}\, [{^B\boldsymbol{\omega}_k}]_{\times}\\ \Rightarrow\qquad{_B^G{\mathbf{R}}}_{k+1} &= {_B^G{\mathbf{R}}_{k}} + {_B^{G}\mathbf{R}_k}{\Delta t} [{^B\boldsymbol{\omega}_k}]_{\times}\\ \Rightarrow\qquad{_B^G{\mathbf{R}}}_{k+1} &= {_B^G{\mathbf{R}}_{k}} \left( \mathbf{I}+ {\Delta t} [{^B\boldsymbol{\omega}}_k]_{\times}\right)\\ ({\small\text{1st-Order Taylor Approximation}})\Rightarrow\qquad {_B^G{\mathbf{R}}}_{k+1} &= {_B^G{\mathbf{R}}_{k}}\,{\rm exp}({\Delta t} [{^B\boldsymbol{\omega}}_k]_{\times})\\ \Rightarrow \qquad {_B^G{\mathbf{R}}}_{k+1} &= {_B^G{\mathbf{R}}_{k}}\,{\rm Exp}({\Delta t} \, {^B\boldsymbol{\omega}}_k)\\ \end{aligned} \tag{2-3}
ΔtBGRk+1−BGRk⇒BGRk+1⇒BGRk+1(1st-Order Taylor Approximation)⇒BGRk+1⇒BGRk+1=BGRk[Bωk]×=BGRk+BGRkΔt[Bωk]×=BGRk(I+Δt[Bωk]×)=BGRkexp(Δt[Bωk]×)=BGRkExp(ΔtBωk)(2-3)
即原文中
x
k
+
1
=
x
k
⋅
E
x
p
(
Δ
t
B
ω
k
)
\mathbf{x}_{k+1} = \mathbf{x}_{k}\cdot {\rm Exp}({\Delta t} \, {^B\boldsymbol{\omega}}_k)
xk+1=xk⋅Exp(ΔtBωk)
这样就得到了原文中式 (13).
题外话一: right- ⊕ \oplus ⊕ 和 left- ⊕ \oplus ⊕
文献 [4] 分别针对相对坐标系下的切向量扰动和全局坐标系下的切向量扰动定义了 right- ⊕ \oplus ⊕、right- ⊖ \ominus ⊖、left- ⊕ \oplus ⊕、left- ⊖ \ominus ⊖.
The right operators are
right- ⊕ : Y = X ⊕ X τ ≜ X ∘ E x p ( X τ ) ∈ M right- ⊖ : X τ = Y ⊖ X ≜ L o g ( X − 1 ∘ Y ) ∈ T X M \begin{aligned} \text{right-} \oplus &: \quad \mathcal{Y} = \mathcal{X} \oplus {^{\mathcal{X}}\tau} \triangleq \mathcal{X}\circ \rm{Exp}({^\mathcal{X}}\tau) \in \mathcal{M}\\ \text{right-} \ominus &: \quad {^{\mathcal{X}}\tau} = \mathcal{Y} \ominus \mathcal{X} \triangleq \rm{Log}( \mathcal{X}^{-1}\circ {\mathcal{Y}} ) \in T_{\mathcal{X}}\mathcal{M} \end{aligned} right-⊕right-⊖:Y=X⊕Xτ≜X∘Exp(Xτ)∈M:Xτ=Y⊖X≜Log(X−1∘Y)∈TXM
where X τ {^\mathcal{X}\tau} Xτ belongs to the tangent space at X \mathcal{X} X. That is, X τ {^\mathcal{X}\tau} Xτ is expressed in the local frame at X \mathcal{X} X.
The left operators are
left- ⊕ : Y = ε τ ⊕ X ≜ E x p ( ε τ ) ∘ X ∈ M left- ⊖ : ϵ τ = Y ⊖ X ≜ L o g ( Y ∘ X − 1 ) ∈ T ε M \begin{aligned} \text{left-} \oplus &: \quad \mathcal{Y} = {^{\boldsymbol{\varepsilon}}\tau} \oplus\mathcal{X} \triangleq \rm{Exp}({^{\boldsymbol{\varepsilon}}}\tau) \circ \mathcal{X} \in \mathcal{M}\\ \text{left-} \ominus &: \quad {^{\boldsymbol{\epsilon}}\tau} = \mathcal{Y} \ominus \mathcal{X} \triangleq \rm{Log}( {\mathcal{Y}} \circ \mathcal{X}^{-1} ) \in T_{\boldsymbol{\varepsilon}}\mathcal{M} \end{aligned} left-⊕left-⊖:Y=ετ⊕X≜Exp(ετ)∘X∈M:ϵτ=Y⊖X≜Log(Y∘X−1)∈TεM
We say that ε τ {^{\boldsymbol{\varepsilon}}\tau} ετ is expressed in the global frame.
式 (1-10) 中的 G ω k {^G\boldsymbol{\omega}_k} Gωk 是全局坐标系下的角速度 (角速度向量相对于姿态矩阵, 视作为扰动), 故而是 left- ⊕ \oplus ⊕.
式 (2-3) 中的 B ω k {^B\boldsymbol{\omega}}_k Bωk 是运动物体局部坐标系下的角速度, 故而是 right- ⊕ \oplus ⊕.
翻来覆去, 这些定义是自洽的.
III. About eq (15) in “Example 5: Vectors of known magnitude in body frame”
1. 姿态矩阵的导数 Time Derivative of an Attitude Matrix
重力向量在全局坐标系 {G} 下的描述记为
G
g
∈
S
2
(
∥
g
∥
)
{^G\mathbf{g}}\in \mathbb{S}^2(\parallel\mathbf{g}\parallel)
Gg∈S2(∥g∥). 重力向量在局部坐标系 {B} 下的描述记为
B
g
^B\mathbf{g}
Bg (原文中记为
x
\mathbf{x}
x). 两者之间的坐标转换关系
G
g
=
B
G
R
B
g
(3-1)
^G\mathbf{g} = {_B^G\mathbf{R}}\,{^B\mathbf{g}} \tag{3-1}
Gg=BGRBg(3-1)
重力向量在全局坐标系下是不变的, 即
G
g
˙
=
0
(3-2)
^G\dot{\mathbf{g}}=\mathbf{0} \tag{3-2}
Gg˙=0(3-2)
对式 (3-1) 两边求导可得
B
G
R
˙
B
g
+
B
G
R
B
g
˙
=
0
(3-3)
{_B^G\dot{\mathbf{R}}}\,{^B\mathbf{g}} + {_B^G\mathbf{R}}\,{^B\dot{\mathbf{g}}} =\mathbf{0} \tag{3-3}
BGR˙Bg+BGRBg˙=0(3-3)
式 (2-2) 代入式 (3-3) 得
B
G
R
[
B
ω
]
×
B
g
+
B
G
R
B
g
˙
=
0
(3-4)
{_B^{G}\mathbf{R}}\, [{^B\boldsymbol{\omega}}]_{\times} \,{^B\mathbf{g}} + {_B^G\mathbf{R}}\,{^B\dot{\mathbf{g}}} =\mathbf{0} \tag{3-4}
BGR[Bω]×Bg+BGRBg˙=0(3-4)
两边同乘以
B
G
R
T
{^G_B\mathbf{R}^{\rm\small T}}
BGRT 可得
B
g
˙
=
−
[
B
ω
]
×
B
g
(3-5)
{^B\dot{\mathbf{g}}} = - [{^B\boldsymbol{\omega}}]_{\times} \,{^B\mathbf{g}} \tag{3-5}
Bg˙=−[Bω]×Bg(3-5)
即原文中的
x
˙
=
−
⌊
B
ω
⌋
x
\dot{\mathbf{x}} = -\lfloor ^B\boldsymbol{\omega}\rfloor \mathbf{x}
x˙=−⌊Bω⌋x
2. 离散化 Zero-Order Hold Discretization
对式 (3-5) 离散化得到
B
g
k
+
1
−
B
g
k
Δ
t
=
−
[
B
ω
k
]
×
B
g
k
⇒
B
g
k
+
1
=
B
g
k
−
Δ
t
[
B
ω
k
]
×
B
g
k
⇒
B
g
k
+
1
=
(
I
−
Δ
t
[
B
ω
k
]
×
)
B
g
k
(
1st-Order Taylor Approximation
)
⇒
B
g
k
+
1
=
E
x
p
(
−
Δ
t
B
ω
k
)
B
g
k
(3-6)
\begin{aligned} \frac{{^B\mathbf{g}_{k+1}} - {^B\mathbf{g}_{k}}}{\Delta t} &= - [{^B{\boldsymbol{\omega}}_k}]_{\times} \,{^B\mathbf{g}_k}\\ \Rightarrow \qquad {^B\mathbf{g}_{k+1}} & = {^B\mathbf{g}_{k}} - {\Delta t} \,[{^B{\boldsymbol{\omega}}_k}]_{\times} \,{^B\mathbf{g}_k}\\ \Rightarrow \qquad {^B\mathbf{g}_{k+1}} & = \left(\mathbf{I}- {\Delta t} \,[{^B{\boldsymbol{\omega}}_k}]_{\times} \right){^B\mathbf{g}_k}\\ ({\small\text{1st-Order Taylor Approximation}}) \quad \Rightarrow\qquad {^B\mathbf{g}_{k+1}} &= {\rm Exp}(- {\Delta t} \,{^B{\boldsymbol{\omega}}_k} ){^B\mathbf{g}_k} \end{aligned}\tag{3-6}
ΔtBgk+1−Bgk⇒Bgk+1⇒Bgk+1(1st-Order Taylor Approximation)⇒Bgk+1=−[Bωk]×Bgk=Bgk−Δt[Bωk]×Bgk=(I−Δt[Bωk]×)Bgk=Exp(−ΔtBωk)Bgk(3-6)
即原文中的
x
k
+
1
=
E
x
p
(
−
Δ
t
B
ω
k
)
x
k
{\mathbf{x}_{k+1}} = {\rm Exp}(- {\Delta t} \,{^B{\boldsymbol{\omega}}_k} ){\mathbf{x}_k}
xk+1=Exp(−ΔtBωk)xk
这样就得到了原文中式 (15).
IV. About eqs (16)~(17) in “Example 6: Bearing-distance parameterization of visual landmarks”
视觉特征点 (Visual landmark) 在相机坐标系 {C} 中的方位向量 (Bearing vector) 记作 x ∈ S 2 ( 1 ) \mathbf{x}\in \mathbb{S}^2(1) x∈S2(1) (即 x ∈ R 3 \mathbf{x}\in \mathbb{R}^3 x∈R3 且 ∣ x ∣ = 1 |\mathbf{x}|=1 ∣x∣=1);
特征点在相机坐标系 {C} 中的深度 (Depth) 记作 d ( ρ ) ∈ R d(\rho)\in \mathbb{R} d(ρ)∈R, 为基于参数 ρ ∈ R \rho \in \mathbb{R} ρ∈R 的标量.
相机在全局坐标系 {G} 中的姿态矩阵记作 C G R {_{C}^{G}{\mathbf{R}}} CGR (原文中记作 G R C {^{G}{\mathbf{R}}_{C}} GRC);
相机在全局坐标系 {G} 中的位置向量记作
C
G
P
{_{C}^{G} \mathbf{P}}
CGP (原文中记作
G
P
C
{^{G}\mathbf{P}_{C}}
GPC).
1. 方位向量的性质 Properties of a Bearing Vector
根据 Fig 1 所示, 特征点在相机坐标系 {C} 中的方位向量
x
=
[
s
i
n
ϕ
c
o
s
θ
s
i
n
ϕ
s
i
n
θ
c
o
s
ϕ
]
(4-1)
\mathbf{x} = \begin{bmatrix} {\rm sin}\phi \,{\rm cos}\theta \\ {\rm sin}\phi \,{\rm sin}\theta \\ {\rm cos}\phi\end{bmatrix} \tag{4-1}
x=
sinϕcosθsinϕsinθcosϕ
(4-1)
对方位向量求时间导数, 得到
x
˙
=
[
c
o
s
ϕ
c
o
s
θ
ϕ
˙
−
s
i
n
ϕ
s
i
n
θ
θ
˙
c
o
s
ϕ
s
i
n
θ
ϕ
˙
+
s
i
n
ϕ
c
o
s
θ
θ
˙
−
s
i
n
ϕ
ϕ
˙
]
(4-2)
\dot{\mathbf{x}} = \begin{bmatrix} {\rm cos}\phi \,{\rm cos}\theta \, \dot{\phi} - {\rm sin}\phi \,{\rm sin}\theta \, \dot{\theta}\\ {\rm cos}\phi \,{\rm sin}\theta \, \dot{\phi} + {\rm sin}\phi \,{\rm cos}\theta \,\dot{\theta}\\ -{\rm sin}\phi\, \dot{\phi} \end{bmatrix} \tag{4-2}
x˙=
cosϕcosθϕ˙−sinϕsinθθ˙cosϕsinθϕ˙+sinϕcosθθ˙−sinϕϕ˙
(4-2)
计算
x
T
x
˙
=
[
s
i
n
ϕ
c
o
s
θ
s
i
n
ϕ
s
i
n
θ
c
o
s
ϕ
]
[
c
o
s
ϕ
c
o
s
θ
ϕ
˙
−
s
i
n
ϕ
s
i
n
θ
θ
˙
c
o
s
ϕ
s
i
n
θ
ϕ
˙
+
s
i
n
ϕ
c
o
s
θ
θ
˙
−
s
i
n
ϕ
ϕ
˙
]
=
0
(4-3)
\mathbf{x}^{\small\rm T} \dot{\mathbf{x}} = \begin{bmatrix} {\rm sin}\phi \,{\rm cos}\theta & {\rm sin}\phi \,{\rm sin}\theta & {\rm cos}\phi\end{bmatrix} \begin{bmatrix} {\rm cos}\phi \,{\rm cos}\theta \, \dot{\phi} - {\rm sin}\phi \,{\rm sin}\theta \, \dot{\theta}\\ {\rm cos}\phi \,{\rm sin}\theta \, \dot{\phi} + {\rm sin}\phi \,{\rm cos}\theta \,\dot{\theta}\\ -{\rm sin}\phi\, \dot{\phi} \end{bmatrix} = 0 \tag{4-3}
xTx˙=[sinϕcosθsinϕsinθcosϕ]
cosϕcosθϕ˙−sinϕsinθθ˙cosϕsinθϕ˙+sinϕcosθθ˙−sinϕϕ˙
=0(4-3)
另外根据 机器人中常用矩阵等式-I (汇总) 中 “[Identity 12] Identity About Skew-Symmetric Matrix” 可知
x
T
[
C
ω
]
×
x
=
0
(4-4)
\mathbf{x}^{\small\rm T} [{^C\boldsymbol{\omega}}]_{\times}\mathbf{x} = 0 \tag{4-4}
xT[Cω]×x=0(4-4)
2. 深度参数的导数 Time Derivative of a Depth Parameter
特征点在全局坐标系中的位置为
C
G
R
(
x
d
(
ρ
)
)
+
C
G
P
{_C^G\mathbf{R}} (\mathbf{x}d(\rho))+{_{C}^G\mathbf{P}}
CGR(xd(ρ))+CGP. 为固定于全局坐标系中的几何点, 所以位置坐标对时间求导为零, 即
d
(
C
G
R
(
x
d
(
ρ
)
)
+
C
G
P
)
d
t
=
0
⇒
C
G
R
˙
(
x
d
(
ρ
)
)
+
C
G
R
(
x
˙
d
(
ρ
)
)
+
C
G
R
(
x
d
′
(
ρ
)
ρ
˙
)
+
C
G
P
˙
=
0
(4-5)
\begin{aligned} \frac{d({{_C^G\mathbf{R}} (\mathbf{x}d(\rho))+{_{C}^G\mathbf{P}}})}{dt} &= \mathbf{0}\\ \Rightarrow \qquad {{_C^G\dot{\mathbf{R}}} (\mathbf{x}d(\rho))+{{_C^G\mathbf{R}} (\dot{\mathbf{x}}d(\rho))}+{{_C^G\mathbf{R}} (\mathbf{x}} d'(\rho)\dot{\rho}) + {_{C}^G\dot{\mathbf{P}}}} &= \mathbf{0} \end{aligned}\tag{4-5}
dtd(CGR(xd(ρ))+CGP)⇒CGR˙(xd(ρ))+CGR(x˙d(ρ))+CGR(xd′(ρ)ρ˙)+CGP˙=0=0(4-5)
由式 (2-2) 可知
C
G
R
˙
=
C
G
R
[
C
ω
]
×
(4-6)
{_C^G\dot{\mathbf{R}}} = {_C^{G}\mathbf{R}}\, [{^C\boldsymbol{\omega}}]_{\times} \tag{4-6}
CGR˙=CGR[Cω]×(4-6)
式 (4-6) 代入式 (4-5) 得到
C
G
R
[
C
ω
]
×
‾
x
d
(
ρ
)
+
C
G
R
x
˙
d
(
ρ
)
+
C
G
R
x
d
′
(
ρ
)
ρ
˙
+
C
G
P
˙
=
0
(4-7)
\underline{{_C^{G}\mathbf{R}} \, [{^C\boldsymbol{\omega}}]_{\times}}\, \mathbf{x}\,d(\rho) + {{_C^G\mathbf{R}}\, \dot{\mathbf{x}} \,d(\rho)}+{{_C^G\mathbf{R}} \,\mathbf{x}}\, d'(\rho)\,\dot{\rho} + {_{C}^G\dot{\mathbf{P}}} = \mathbf{0} \tag{4-7}
CGR[Cω]×xd(ρ)+CGRx˙d(ρ)+CGRxd′(ρ)ρ˙+CGP˙=0(4-7)
左右两边都左乘
C
G
R
T
{_C^{G}\mathbf{R}}^{\small\rm T}
CGRT, 并利用姿态矩阵的正交性, 得到
[
C
ω
]
×
x
d
(
ρ
)
+
x
˙
d
(
ρ
)
+
x
d
′
(
ρ
)
ρ
˙
+
C
G
R
T
C
G
P
˙
=
0
(4-8)
[{^C\boldsymbol{\omega}}]_{\times}\, \mathbf{x}\,d(\rho) + \dot{\mathbf{x}} \,d(\rho) + \mathbf{x}\, d'(\rho)\,\dot{\rho} + {{_C^G{\mathbf{R}}}^{\small\rm T}} \, {_{C}^G\dot{\mathbf{P}}} = \mathbf{0} \tag{4-8}
[Cω]×xd(ρ)+x˙d(ρ)+xd′(ρ)ρ˙+CGRTCGP˙=0(4-8)
定义相机相对于全局坐标系 {G} 的运动速度向量在相机坐标系 {C} 中进行描述所得的速度向量记为
C
v
≜
G
C
R
C
G
P
˙
=
C
G
R
T
C
G
P
˙
(4-9)
^C{\mathbf{v}} \triangleq {_G^C{\mathbf{R}}} \,{_{C}^G\dot{\mathbf{P}}} = {_C^G{\mathbf{R}}}^{\small\rm T}\, {_{C}^G\dot{\mathbf{P}}} \tag{4-9}
Cv≜GCRCGP˙=CGRTCGP˙(4-9)
式 (4-9) 代入式 (4-8)
[
C
ω
]
×
x
d
(
ρ
)
+
x
˙
d
(
ρ
)
+
x
d
′
(
ρ
)
ρ
˙
+
C
v
=
0
(4-10)
[{^C\boldsymbol{\omega}}]_{\times}\, \mathbf{x}\,d(\rho) + \dot{\mathbf{x}} \,d(\rho) + \mathbf{x}\, d'(\rho)\,\dot{\rho} + {^C{\mathbf{v}}} = \mathbf{0} \tag{4-10}
[Cω]×xd(ρ)+x˙d(ρ)+xd′(ρ)ρ˙+Cv=0(4-10)
式 (4-10) 两边左乘
x
T
\mathbf{x}^{\small\rm T}
xT 得到标量方程
x
T
[
C
ω
]
×
x
d
(
ρ
)
+
x
T
x
˙
d
(
ρ
)
+
x
T
x
d
′
(
ρ
)
ρ
˙
+
x
T
C
v
=
0
(4-11)
\mathbf{x}^{\small\rm T}\,[{^C\boldsymbol{\omega}}]_{\times}\, \mathbf{x}\,d(\rho) + \mathbf{x}^{\small\rm T}\,\dot{\mathbf{x}} \,d(\rho) + \mathbf{x}^{\small\rm T}\,\mathbf{x}\, d'(\rho)\,\dot{\rho} + \mathbf{x}^{\small\rm T}\, {^C{\mathbf{v}}} = \mathbf{0} \tag{4-11}
xT[Cω]×xd(ρ)+xTx˙d(ρ)+xTxd′(ρ)ρ˙+xTCv=0(4-11)
由式 (4-3) 及式 (4-4) 可知式 (4-11) 的第一、二项为零, 以及
x
T
x
=
1
\mathbf{x}^{\small\rm T} \mathbf{x} = 1
xTx=1 可知
d
′
(
ρ
)
ρ
˙
+
x
T
C
v
=
0
(4-12)
d'(\rho)\,\dot{\rho} + \mathbf{x}^{\small\rm T}\, {^C{\mathbf{v}}} = \mathbf{0} \tag{4-12}
d′(ρ)ρ˙+xTCv=0(4-12)
最后得到
ρ
˙
=
−
x
T
C
v
/
d
′
(
ρ
)
(4-13)
\dot{\rho} = - \mathbf{x}^{\small\rm T} \, {^C{\mathbf{v}}}/d'(\rho) \tag{4-13}
ρ˙=−xTCv/d′(ρ)(4-13)
3. 方位向量的时间导数 Time Derivative of a Bearing Vector
将式 (4-13) 代入式 (4-10) 得到
[
C
ω
]
×
x
d
(
ρ
)
+
x
˙
d
(
ρ
)
−
x
x
T
C
v
+
C
v
=
0
⇒
x
˙
=
−
[
C
ω
]
×
x
+
1
d
(
ρ
)
(
x
x
T
−
I
)
C
v
(4-14)
\begin{aligned} \left[{^C\boldsymbol{\omega}}\right]_{\times}\, \mathbf{x}\,d(\rho) + \dot{\mathbf{x}} \,d(\rho) - \mathbf{x}\, \mathbf{x}^{\small\rm T}\, {^C{\mathbf{v}}} + {^C{\mathbf{v}}} = \mathbf{0}\\ \Rightarrow\qquad \dot{\mathbf{x}} = -\left[{^C\boldsymbol{\omega}}\right]_{\times}\, \mathbf{x} + \frac{1}{d(\rho)}\,\left(\mathbf{x}\, \mathbf{x}^{\small\rm T} - \mathbf{I}\right)\,{^C{\mathbf{v}}} \end{aligned} \tag{4-14}
[Cω]×xd(ρ)+x˙d(ρ)−xxTCv+Cv=0⇒x˙=−[Cω]×x+d(ρ)1(xxT−I)Cv(4-14)
其中
(
x
x
T
−
I
)
C
v
(
[
x
]
×
2
=
x
x
T
−
I
)
=
[
x
]
×
2
C
v
Associative law of multiplication of matrices
=
[
x
]
×
(
[
x
]
×
C
v
)
(
[
x
]
×
y
=
−
[
y
]
×
x
)
=
−
[
[
x
]
×
C
v
]
×
x
(4-15)
\begin{aligned} &\; \left(\mathbf{x}\, \mathbf{x}^{\small\rm T} - \mathbf{I}\right)\,{^C{\mathbf{v}}}\\ {\small ([\mathbf{x}]_{\times}^{2} = \mathbf{x}\, \mathbf{x}^{\small\rm T} - \mathbf{I})} \qquad =& \; [\mathbf{x}]_{\times}^{2} {^C\mathbf{v}} \\ {\small\text{Associative law of multiplication of matrices}}\qquad = &\; [\mathbf{x}]_{\times} \left([\mathbf{x}]_{\times}\, {^C\mathbf{v}}\right) \\ {\small ([\mathbf{x}]_{\times} \mathbf{y} = -[\mathbf{y}]_{\times} \mathbf{x} )} \qquad =& \; - \left[[\mathbf{x}]_{\times}\, {^C\mathbf{v}}\right]_\times \mathbf{x}\\ \end{aligned} \tag{4-15}
([x]×2=xxT−I)=Associative law of multiplication of matrices=([x]×y=−[y]×x)=(xxT−I)Cv[x]×2Cv[x]×([x]×Cv)−[[x]×Cv]×x(4-15)
式 (4-15) 代入式 (4-14) 得到
x
˙
=
−
[
C
ω
]
×
x
−
1
d
(
ρ
)
[
[
x
]
×
C
v
]
×
x
=
−
(
[
C
ω
]
×
+
1
d
(
ρ
)
[
[
x
]
×
C
v
]
×
)
x
=
−
[
C
ω
+
1
d
(
ρ
)
[
x
]
×
C
v
]
×
x
(4-16)
\begin{aligned} \dot{\mathbf{x}} &= -\left[{^C\boldsymbol{\omega}}\right]_{\times}\, \mathbf{x} - \frac{1}{d(\rho)} \left[[\mathbf{x}]_{\times}\, {^C\mathbf{v}}\right]_\times \mathbf{x}\\ &= -\left(\left[{^C\boldsymbol{\omega}}\right]_{\times} + \frac{1}{d(\rho)} \left[[\mathbf{x}]_{\times}\, {^C\mathbf{v}}\right]_\times\right) \mathbf{x}\\ &= - \left[{^C\boldsymbol{\omega}} + \frac{1}{d(\rho)} [\mathbf{x}]_{\times}\, {^C\mathbf{v}}\right]_\times\, \mathbf{x} \end{aligned} \tag{4-16}
x˙=−[Cω]×x−d(ρ)1[[x]×Cv]×x=−([Cω]×+d(ρ)1[[x]×Cv]×)x=−[Cω+d(ρ)1[x]×Cv]×x(4-16)
推到得到和原文中 (17) 一致的式子.
4. 离散化 Zero-Order Hold Discretization
先简化式 (4-16)
Ω
≜
C
ω
+
1
d
(
ρ
)
[
x
]
×
C
v
(4-17)
\boldsymbol{\Omega} \triangleq {^C\boldsymbol{\omega}} + \frac{1}{d(\rho)} [\mathbf{x}]_{\times}\, {^C\mathbf{v}} \tag{4-17}
Ω≜Cω+d(ρ)1[x]×Cv(4-17)
则式 (4-16) 简写为
x
˙
=
−
[
Ω
]
×
x
(4-18)
\dot{\mathbf{x}} = - \left[\boldsymbol{\Omega} \right]_\times\, \mathbf{x} \tag{4-18}
x˙=−[Ω]×x(4-18)
对式 (4-18) 离散化得到
x
k
+
1
−
x
k
Δ
t
=
−
[
Ω
k
]
×
x
k
⇒
x
k
+
1
=
x
k
−
Δ
t
[
Ω
k
]
×
x
k
⇒
x
k
+
1
=
(
I
−
Δ
t
[
Ω
k
]
×
)
x
k
(
1st-Order Taylor Approximation
)
⇒
x
k
+
1
=
E
x
p
(
−
Δ
t
Ω
k
)
x
k
(4-19)
\begin{aligned} \frac{{\mathbf{x}_{k+1}} - {\mathbf{x}_{k}}}{\Delta t} &= - [{{\boldsymbol{\Omega}}_k}]_{\times} \,{\mathbf{x}_k}\\ \Rightarrow \qquad {\mathbf{x}_{k+1}} &= {\mathbf{x}_{k}} - {\Delta t} [{{\boldsymbol{\Omega}}_k}]_{\times} \,{\mathbf{x}_k}\\ \Rightarrow \qquad {\mathbf{x}_{k+1}} &= \left( {\mathbf{I}} - {\Delta t} [{{\boldsymbol{\Omega}}_k}]_{\times} \right) {\mathbf{x}_k}\\ ({\small\text{1st-Order Taylor Approximation}}) \quad \Rightarrow\qquad {\mathbf{x}_{k+1}} &= {\rm Exp}( - {\Delta t} {{\boldsymbol{\Omega}}_k} )\, {\mathbf{x}_k}\\ \end{aligned}\tag{4-19}
Δtxk+1−xk⇒xk+1⇒xk+1(1st-Order Taylor Approximation)⇒xk+1=−[Ωk]×xk=xk−Δt[Ωk]×xk=(I−Δt[Ωk]×)xk=Exp(−ΔtΩk)xk(4-19)
式 (4-17) 代入式 (4-19) 得到
x
k
+
1
=
E
x
p
(
−
Δ
t
(
C
ω
k
+
1
d
(
ρ
)
[
x
k
]
×
C
v
k
)
)
x
k
(4-20)
{\mathbf{x}_{k+1}} = {\rm Exp}\left( - {\Delta t} \left({^C\boldsymbol{\omega}_k} + \frac{1}{d(\rho)} [\mathbf{x}_k]_{\times}\, {^C\mathbf{v}_k}\right) \right)\, {\mathbf{x}_k} \tag{4-20}
xk+1=Exp(−Δt(Cωk+d(ρ)1[xk]×Cvk))xk(4-20)
和原文中 (17) 一致.
题外话二: S 2 \mathbb{S}^2 S2 is not a Lie group 三维空间的球面 (2-Sphere) 不是李群
文献 [5] 中说明了 S 2 \mathbb{S}^2 S2 不是李群, 也就是上面提到的重力向量、方位向量都不是李群. 这里鹦鹉学舌讲一下大概原理.
毛球定理 (Hairy ball theorem): 在三维空间中球面 S 2 \mathbb{S}^2 S2 上不存在处处连续的非零并且与球面相切的向量场. 比如一个长满毛圆球, 把毛都贴在圆球表面理顺的话, 必然有 “漩”. 再比如, 地球表面永远存在气旋和风眼, 在风眼处风平浪静, 但四周都有风环绕.
而李群上都存在一个处处连续且非零的切向量场. 这就和毛球定理矛盾, 故而三维空间的球面 (2-Sphere) 不是李群.
因为不是李群, 所以重力向量相关的式 (3-6) 和方位向量相关的式 (4-20) 与标准 right- ⊕ \oplus ⊕ 和 left- ⊕ \oplus ⊕ 形式有所不同.
总结
本文主要是阅读论文 “Kalman Filters on Differentiable Manifolds”[1] 过程中, 对
- S O ( 3 ) SO(3) SO(3) 上的姿态运动学
- S 2 \mathbb{S}^2 S2 上重力向量的运动学
- S 2 \mathbb{S}^2 S2 上方位向量的运动学
进行的详细啰嗦版推导记录.
参考文献
[1] Dongjiao He, Wei Xu, Fu Zhang, “Kalman Filters on Differentiable Manifolds”, https://arxiv.org/abs/2102.03804, 2021, Version V3
[2] Shiyu Zhao, “Time Derivative of Rotation Matrices: A Tutorial”, https://arxiv.org/abs/1609.06088, 2016
[3] wzf@robotics_notes, “机器人中常用矩阵等式-I (汇总)”, https://blog.csdn.net/woyaomaishu2/article/details/132351402?spm=1001.2014.3001.5502
[4] Joan Solà, Jeremie Deray, Dinesh Atchuthan, “A micro Lie theory for state estimation in robotics”, https://arxiv.org/abs/1812.01537, 2021, Version v9
[5] Seewoo Lee, “ S 2 S^2 S2 is not a Lie group”, https://math.berkeley.edu/~seewoo5/S2notLiegrp.pdf#:~:text=Since%20%28S2%29%20%3D%202%2C%20it,can%27t%20admit%20a%20Lie%20group%20structure, 2018