The plus operator
设 M \mathcal{M} M表示一个n维的流型,因为流型局部同胚与 R n \mathbb{R}^n Rn,所以我们可以通过定义符号 ⊞ \boxplus ⊞和 ⊟ \boxminus ⊟建立一个流型 M \mathcal{M} M的局部邻域和其正切空间的双射。
⊞ : M × R n → M ; ⊟ : M × M n → R n \boxplus:\mathcal{M}\times \mathbb{R}^n\to \mathcal{M};\quad\boxminus:\mathcal{M}\times \mathbb{M}^n\to \mathbb{R}^n ⊞:M×Rn→M;⊟:M×Mn→Rn
对于流型
S
O
(
3
)
SO(3)
SO(3),
⊞
:
S
O
(
3
)
×
R
n
→
S
O
(
3
)
\boxplus:SO(3)\times \mathbb{R}^n\to SO(3)
⊞:SO(3)×Rn→SO(3)生成一个
S
O
(
3
)
SO(3)
SO(3)中的元素,其结果是参考元素
R
∈
S
O
(
3
)
\bold{R}\in SO(3)
R∈SO(3)和一个(通常很小的)旋转的组合。这个旋转由流型在参考元素
R
\bold{R}
R处的正切空间中的一个向量
θ
∈
R
3
\theta\in\mathbb{R}^3
θ∈R3指定。
M
=
S
O
(
3
)
:
R
⊞
θ
=
R
E
x
p
(
θ
)
\mathcal{M}=SO(3):\bold{R}\boxplus\bm{\theta}=\bold{R}{\rm Exp}(\bm{\theta})
M=SO(3):R⊞θ=RExp(θ)
minus符号
⊟
:
S
O
(
3
)
×
S
O
(
3
)
→
R
3
\boxminus:SO(3)\times SO(3)\to \mathbb{R}^3
⊟:SO(3)×SO(3)→R3是plus符号的逆运算,它返回两个
S
O
(
3
)
SO(3)
SO(3)元素的角度差异向量
θ
∈
R
3
\bm{\theta}\in\mathbb{R}^3
θ∈R3,角度差异向量在参考元素
R
\bold{R}
R的正切空间中表示。
注意到这个符号可以定义于任意SO(3)的表示:
q
⊞
θ
=
q
⊗
E
x
p
(
θ
)
θ
=
q
S
⊟
q
R
=
L
o
g
(
q
R
∗
⊗
q
S
)
\bold{q}\boxplus \bm{\theta}=\bold{q}\otimes{\rm Exp}(\bm{\theta})\\ \bm{\theta}=\bold{q}_S\boxminus\bold{q}_R={\rm Log}(\bold{q}^{\ast}_R\otimes \bold{q}_S)
q⊞θ=q⊗Exp(θ)θ=qS⊟qR=Log(qR∗⊗qS)
在这两种情况下,请注意,即使向量差
θ
\bm{\theta}
θ通常被认为是很小的,上面的定义适用于
θ
\bm{\theta}
θ的任何值(直到SO(3)流形的第一个覆盖范围,也就是说,对于角
θ
<
π
\theta < \pi
θ<π
四种可能的导数定义
Functions from vector space to vector space
给定一个函数
f
:
R
m
→
R
n
f:\mathbb{R}^m\to\mathbb{R}^{n}
f:Rm→Rn
∂
f
(
x
)
∂
x
=
lim
δ
x
→
0
f
(
x
+
δ
x
)
−
f
(
x
)
δ
x
∈
R
n
×
m
f
(
x
+
Δ
x
)
≈
f
(
x
)
+
∂
f
(
x
)
∂
x
Δ
x
∈
R
n
\dfrac{\partial f\left( \bold{x}\right) }{\partial \bold{x}}=\lim _{\delta \bold{x}\rightarrow 0}\dfrac{f\left( \bold{x}+\delta \bold{x}\right) -f\left( \bold{x}\right) }{\delta \bold{x}} \in\mathbb{R}^{n\times m} \\ f\left( \bold{x}+\Delta \bold{x}\right) \approx f\left( \bold{x}\right) +\dfrac{\partial f\left( \bold{x}\right) }{\partial \bold{x}}\Delta \bold{x}\in \mathbb{R}^n
∂x∂f(x)=δx→0limδxf(x+δx)−f(x)∈Rn×mf(x+Δx)≈f(x)+∂x∂f(x)Δx∈Rn
Functions from S O ( 3 ) SO(3) SO(3) to S O ( 3 ) SO(3) SO(3)
给定函数
f
:
S
O
(
3
)
→
S
O
(
3
)
f:SO(3)\to SO(3)
f:SO(3)→SO(3),
R
∈
S
O
(
3
)
\bold{R} \in SO(3)
R∈SO(3)和一个局部的小角度变化
θ
∈
R
3
\bm{\theta}\in\mathbb{R}^3
θ∈R3,我们使用
{
⊞
,
⊟
}
\{\boxplus,\boxminus\}
{⊞,⊟}定义导数:
∂
f
(
R
)
∂
θ
=
lim
δ
θ
→
0
f
(
R
⊞
δ
θ
)
⊟
f
(
R
)
δ
θ
=
lim
δ
θ
→
0
L
o
g
(
f
−
1
(
R
)
f
(
R
E
x
p
(
δ
θ
)
)
)
δ
θ
∈
R
3
×
3
f
(
R
⊞
Δ
θ
)
≈
f
(
R
)
⊞
∂
f
(
R
)
∂
θ
Δ
θ
≈
f
(
R
)
E
x
p
(
∂
f
(
R
)
∂
θ
Δ
θ
)
∈
S
O
(
3
)
\begin{aligned} \dfrac{\partial f\left( \bold{R}\right) }{\partial \bm{\theta} }&=\lim _{\delta \bm{\theta} \rightarrow 0}\dfrac{f\left( \bold{R}\boxplus \delta \bm{\theta} \right) \boxminus f\left( \bold{R}\right) }{\delta \bm{\theta} }\\ &=\lim _{\delta \bm{\theta} \rightarrow 0}\dfrac{{\rm Log}\left( f^{-1}\left( \bold{R}\right) f\left( \bold{R} {\rm Exp}\left( \delta \bm{\theta} \right) \right) \right) }{\delta \bm{\theta} }\quad\in \mathbb{R}^{3\times3}\\ f\left( \bold{R}\boxplus \Delta \bm{\theta} \right) &\approx f\left( \bold{R}\right) \boxplus \dfrac{\partial f\left( \bold{R}\right) }{\partial \bm{\theta} }\Delta \bm{\theta} \approx f\left( \bold{R}\right) {\rm Exp}\left( \dfrac{\partial f\left( \bold{R} \right) }{\partial \bm{\theta} }\Delta \bm{\theta} \right) \quad\in SO(3) \end{aligned}
∂θ∂f(R)f(R⊞Δθ)=δθ→0limδθf(R⊞δθ)⊟f(R)=δθ→0limδθLog(f−1(R)f(RExp(δθ)))∈R3×3≈f(R)⊞∂θ∂f(R)Δθ≈f(R)Exp(∂θ∂f(R)Δθ)∈SO(3)
Functions from vector space to S O ( 3 ) SO(3) SO(3)
给定函数
f
:
R
m
→
S
O
(
3
)
f:\mathbb{R}^m\to SO(3)
f:Rm→SO(3),用
⊟
\boxminus
⊟进行
S
O
(
3
)
SO(3)
SO(3)差分,用-
进行向量差分。
∂
f
(
x
)
∂
x
=
lim
δ
x
→
0
f
(
x
+
δ
x
)
⊟
f
(
x
)
δ
x
=
lim
δ
x
→
0
L
o
g
(
f
−
1
(
x
)
f
(
x
+
δ
x
)
)
δ
x
∈
R
3
×
m
f
(
x
+
Δ
x
)
≈
f
(
x
)
⊞
∂
f
(
x
)
∂
x
Δ
x
≈
f
(
x
)
E
x
p
(
∂
f
(
x
)
∂
x
Δ
x
)
∈
S
O
(
3
)
\begin{aligned} \dfrac{\partial f\left( \bold{x}\right) }{\partial \bold{x}}&=\lim _{\delta \bold{x}\rightarrow 0}\dfrac{f\left( \bold{x}+\delta \bold{x}\right) \boxminus f\left( \bold{x}\right) }{\delta \bold{x}}\\ &=\lim _{\delta \bold{x}\rightarrow 0}\dfrac{{\rm Log}\left( f^{-1}\left( \bold{x}\right) f\left( \bold{x}+\delta \bold{x}\right) \right) }{\delta \bold{x}}\quad\in\mathbb{R}^{3\times m}\\ f\left( \bold{x}+\Delta \bold{x}\right) &\approx f\left( \bold{x}\right) \boxplus \dfrac{\partial f\left( \bold{x}\right) }{\partial \bold{x}}\Delta \bold{x}\approx f\left( \bold{x}\right) {\rm Exp}\left( \dfrac{\partial f\left( \bold{x}\right) }{\partial \bold{x}}\Delta \bold{x}\right) \quad \in SO(3) \end{aligned}
∂x∂f(x)f(x+Δx)=δx→0limδxf(x+δx)⊟f(x)=δx→0limδxLog(f−1(x)f(x+δx))∈R3×m≈f(x)⊞∂x∂f(x)Δx≈f(x)Exp(∂x∂f(x)Δx)∈SO(3)
Functions from S O ( 3 ) SO(3) SO(3) to vector space
对于函数
f
:
S
O
(
3
)
→
R
n
f:SO(3)\to \mathbb{R}^{n}
f:SO(3)→Rn,用
⊞
\boxplus
⊞进行
S
O
(
3
)
SO(3)
SO(3)组合,用-
进行向量差分。
∂
f
(
R
)
∂
θ
=
lim
δ
θ
→
0
f
(
R
⊞
δ
θ
)
−
f
(
R
)
δ
θ
=
lim
δ
θ
→
0
f
(
R
E
x
p
(
δ
θ
)
)
−
f
(
R
)
δ
θ
∈
R
n
×
3
f
(
R
⊞
δ
θ
)
≈
f
(
R
)
+
∂
f
(
R
)
∂
θ
Δ
θ
∈
S
O
(
3
)
\begin{aligned} \dfrac{\partial f\left( \bold{R}\right) }{\partial \bm{\theta} }&=\lim _{\delta \bm{\theta} \rightarrow 0}\dfrac{f\left( \bold{R}\boxplus \delta \bm{\theta} \right) -f\left( \bold{R}\right) }{\delta \bm{\theta} }\\ &=\lim _{\delta \bm{\theta} \rightarrow 0}\dfrac{f\left( \bold{R}{\rm Exp}\left( \delta \bm{\theta} \right) \right) -f\left( \bold{R}\right) }{\delta \bm{\theta} } \quad \in \mathbb{R}^{n\times 3}\\ f\left( \bold{R}\boxplus \delta \bm{\theta} \right) &\approx f\left( \bold{R}\right) +\dfrac{\partial f\left( \bold{R}\right) }{\partial \bm{\theta} }\Delta \bm{\theta} \quad \in SO(3) \end{aligned}
∂θ∂f(R)f(R⊞δθ)=δθ→0limδθf(R⊞δθ)−f(R)=δθ→0limδθf(RExp(δθ))−f(R)∈Rn×3≈f(R)+∂θ∂f(R)Δθ∈SO(3)
Jacobians of rotation
Jacobian with respect to the vector
向量旋转对于向量的求导是容易的。
∂
(
q
⊗
p
⊗
q
∗
)
∂
p
=
∂
(
R
p
)
∂
p
=
R
\dfrac{\partial \left( \bold{q}\otimes \bold{p} \otimes \bold{q}^{\ast }\right) }{\partial \bold{p}}=\dfrac{\partial \left( \bold{Rp}\right) }{\partial \bold{p}}=\bold{R}
∂p∂(q⊗p⊗q∗)=∂p∂(Rp)=R
Jacobian with respect to the quaternion
旋转对于四元数的导数是棘手的。为了方便起见,我们使用符号
q
=
[
w
v
]
=
w
+
v
\bold{q}=\begin{bmatrix}\bold{w} & \bold{v}\end{bmatrix}=\bold{w} + \bold{v}
q=[wv]=w+v
p
′
=
q
⊗
p
⊗
q
∗
=
(
w
+
v
)
⊗
p
⊗
(
w
−
v
)
=
w
2
p
+
w
(
v
⊗
p
−
p
⊗
v
)
−
v
⊗
p
⊗
v
=
w
2
p
+
2
w
(
v
×
p
)
−
[
(
−
v
T
p
+
v
×
p
)
⊗
v
]
=
w
2
p
+
2
w
(
v
×
p
)
−
[
(
−
v
T
p
)
v
−
(
v
×
p
)
T
v
+
(
v
×
p
)
×
v
]
=
w
2
p
+
2
w
(
v
×
p
)
−
[
(
−
v
T
p
)
v
+
(
v
T
v
)
p
−
(
v
T
p
)
v
]
=
w
2
p
+
2
w
(
v
×
p
)
+
2
(
v
T
p
)
v
−
(
v
T
v
)
p
\begin{aligned} \bold{p}'&=\bold{q}\otimes \bold{p}\otimes \bold{q}^{\ast }\\ &=\left( \bold{w}+\bold{v}\right) \otimes \bold{p}\otimes \left( \bold{w}-\bold{v}\right) \\ &=\bold{w}^{2}\bold{p}+\bold{w}\left( \bold{v}\otimes \bold{p}-\bold{p}\otimes \bold{v}\right) -\bold{v}\otimes \bold{p}\otimes \bold{v}\\ &=\bold{w}^{2}\bold{p}+2\bold{w}\left( \bold{v}\times \bold{p}\right) -\left[ \left( -\bold{v}^{T}\bold{p}+\bold{v}\times \bold{p}\right) \otimes \bold{v}\right] \\ &=\bold{w}^{2}\bold{p}+2\bold{w}\left( \bold{v}\times \bold{p}\right) -\left[ \left( -\bold{v}^{T}\bold{p}\right) \bold{v}-\bcancel{\left( \bold{v}\times \bold{p}\right) ^{T}\bold{v}}+\left( \bold{v}\times \bold{p}\right) \times \bold{v}\right] \\ &=\bold{w}^{2}\bold{p}+2\bold{w}\left( \bold{v}\times \bold{p}\right) -\left[ \left( -\bold{v}^{T}\bold{p}\right) \bold{v}+\left( \bold{v}^{T}\bold{v}\right) \bold{p}-\left( \bold{v}^{T}\bold{p}\right) \bold{v}\right] \\ &=\bold{w}^{2}\bold{p}+2\bold{w}\left( \bold{v}\times \bold{p}\right) +2\left( \bold{v}^{T}\bold{p}\right) \bold{v}-\left( \bold{v}^{T}\bold{v}\right) \bold{p} \end{aligned}
p′=q⊗p⊗q∗=(w+v)⊗p⊗(w−v)=w2p+w(v⊗p−p⊗v)−v⊗p⊗v=w2p+2w(v×p)−[(−vTp+v×p)⊗v]=w2p+2w(v×p)−[(−vTp)v−(v×p)Tv
+(v×p)×v]=w2p+2w(v×p)−[(−vTp)v+(vTv)p−(vTp)v]=w2p+2w(v×p)+2(vTp)v−(vTv)p
进而我们能够提取出导数
∂
p
′
∂
w
=
2
w
p
+
2
v
×
p
∂
p
′
∂
v
=
−
2
w
p
∧
+
2
(
v
p
T
+
v
T
p
I
)
−
2
p
v
T
=
2
(
v
T
p
I
+
v
p
T
−
p
v
T
−
w
p
∧
)
\begin{aligned} \dfrac{\partial \bold{p}'}{\partial \bold{w}}&=2\bold{wp}+2\bold{v}\times \bold{p}\\ \dfrac{\partial \bold{p}'}{\partial \bold{v}}&=-2\bold{wp}^{\wedge }+2\left( \bold{vp}^{T}+\bold{v}^{T}\bold{pI}\right) -2\bold{pv}^{T}\\ &=2\left( \bold{v}^{T}\bold{pI}+\bold{vp}^{T}-\bold{pv}^{T}-\bold{wp}^{\wedge}\right) \end{aligned}
∂w∂p′∂v∂p′=2wp+2v×p=−2wp∧+2(vpT+vTpI)−2pvT=2(vTpI+vpT−pvT−wp∧)
合并上式能够得出
∂
(
q
⊗
p
⊗
q
∗
)
∂
q
=
2
[
w
p
+
2
v
×
p
v
T
p
I
+
v
p
T
−
p
v
T
−
w
p
∧
]
∈
R
3
×
4
\dfrac{\partial \left( \bold{q}\otimes \bold{p}\otimes \bold{q}^{\ast }\right) }{\partial \bold{q}}=2\begin{bmatrix} \bold{wp}+2\bold{v}\times \bold{p} & \bold{v}^{T}\bold{pI}+\bold{vp}^{T}-\bold{pv}^{T}-\bold{wp}^{\wedge} \end{bmatrix}\in{\mathbb{R}^{3\times 4}}
∂q∂(q⊗p⊗q∗)=2[wp+2v×pvTpI+vpT−pvT−wp∧]∈R3×4
S O ( 3 ) SO(3) SO(3)的右Jacobian矩阵
E
x
p
(
θ
)
⊞
δ
ϕ
=
E
x
p
(
θ
+
δ
θ
)
⟺
E
x
p
(
θ
)
⋅
E
x
p
(
δ
ϕ
)
=
E
x
p
(
θ
+
δ
θ
)
⟺
δ
ϕ
=
L
o
g
(
E
x
p
(
θ
)
−
1
E
x
p
(
θ
+
δ
θ
)
)
=
E
x
p
(
θ
+
δ
θ
)
⊟
E
x
p
(
θ
)
\begin{aligned} &{\rm Exp}\left( \bm{\theta} \right) \boxplus \delta \bm{\phi} ={\rm Exp}\left( \bm{\theta} +\delta \bm{\theta} \right) \\ \iff &{\rm Exp}\left( \bm{\theta} \right) \cdot {\rm Exp}\left( \delta \bm{\phi} \right) ={\rm Exp}\left( \bm{\theta} +\delta \bm{\theta} \right) \\ \iff &\delta \bm{\phi} ={\rm Log}\left( {\rm Exp}\left( \bm{\theta} \right) ^{-1}{\rm Exp}\left( \bm{\theta} +\delta \bm{\theta} \right) \right) ={\rm Exp}\left( \bm{\theta} +\delta \bm{\theta} \right) \boxminus {\rm Exp}\left( \bm{\theta} \right) \end{aligned}
⟺⟺Exp(θ)⊞δϕ=Exp(θ+δθ)Exp(θ)⋅Exp(δϕ)=Exp(θ+δθ)δϕ=Log(Exp(θ)−1Exp(θ+δθ))=Exp(θ+δθ)⊟Exp(θ)
δ
ϕ
\delta \bm{\phi}
δϕ相对于
δ
θ
\delta \bm{\theta}
δθ的微分
∂
δ
ϕ
∂
δ
θ
=
lim
δ
θ
→
0
δ
ϕ
δ
θ
=
lim
δ
θ
→
0
E
x
p
(
θ
+
δ
θ
)
⊟
E
x
p
(
θ
)
δ
θ
\dfrac{\partial \delta \bm{\phi} }{\partial \delta \bm{\theta} }=\lim _{\delta \bm{\theta} \rightarrow 0}\dfrac{\delta \bm{\phi} }{\delta \bm{\theta} }=\lim _{\delta \bm{\theta} \rightarrow 0}\dfrac{{\rm Exp}\left( \bm{\theta} +\delta \bm{\theta} \right) \boxminus {\rm Exp}\left( \bm{\theta} \right) }{\delta \bm{\theta} }
∂δθ∂δϕ=δθ→0limδθδϕ=δθ→0limδθExp(θ+δθ)⊟Exp(θ)
这个Jacobian矩阵就是熟知的
S
O
(
3
)
SO(3)
SO(3)的右Jacobian矩阵,其定义为:
J
r
(
θ
)
=
∂
E
x
p
(
θ
)
∂
θ
=
lim
δ
θ
→
0
E
x
p
(
θ
+
δ
θ
)
⊗
E
x
p
(
θ
)
δ
θ
=
lim
δ
θ
→
0
L
o
g
(
E
x
p
(
θ
)
T
E
x
p
(
θ
+
δ
θ
)
)
δ
θ
if using
R
=
lim
δ
θ
→
0
L
o
g
(
E
x
p
(
θ
)
∗
⊗
E
x
p
(
θ
+
δ
θ
)
)
δ
θ
if using
q
\begin{aligned} \bold{J}_{r}\left( \bm{\theta} \right) &=\dfrac{\partial {\rm Exp}\left( \bm{\theta} \right) }{\partial \bm{\theta} } \\ &=\lim _{\delta \bm{\theta} \rightarrow 0}\dfrac{{\rm Exp}\left( \bm{\theta} +\delta \bm{\theta} \right) \otimes {\rm Exp}\left( \bm{\theta} \right) }{\delta \bm{\theta} } \\ &=\lim _{\delta \bm{\theta} \rightarrow 0}\dfrac{{\rm Log}\left( {\rm Exp}\left( \bm{\theta} \right) ^{T}{\rm Exp}\left( \bm{\theta} +\delta \bm{\theta} \right) \right) }{\delta \bm{\theta} }\qquad \text{if using }\bold{R} \\ &=\lim _{\delta \bm{\theta} \rightarrow 0}\dfrac{{\rm Log} \left( {\rm Exp}\left( \bm{\theta} \right) ^{\ast }\otimes {\rm Exp}\left( \bm{\theta} +\delta \bm{\theta} \right) \right) }{\delta \bm{\theta} } \qquad \text{if using }\bold{q} \end{aligned}
Jr(θ)=∂θ∂Exp(θ)=δθ→0limδθExp(θ+δθ)⊗Exp(θ)=δθ→0limδθLog(Exp(θ)TExp(θ+δθ))if using R=δθ→0limδθLog(Exp(θ)∗⊗Exp(θ+δθ))if using q
右Jacobian矩阵及其逆矩阵的封闭形式可以计算出:
J
r
(
θ
)
=
I
−
1
−
cos
∥
θ
∥
∥
θ
∥
2
θ
∧
+
∥
θ
∥
−
sin
∥
θ
∥
∥
θ
∥
3
θ
∧
2
J
r
−
1
(
θ
)
=
I
+
1
2
θ
∧
+
(
1
∥
θ
∥
2
−
1
+
cos
∥
θ
∥
2
∥
θ
∥
sin
∥
θ
∥
)
θ
∧
2
\bold{J}_{r}\left( \bm{\theta} \right) =I-\dfrac{1-\cos \left\| \bm{\theta} \right\| }{\left\| \bm{\theta} \right\| ^{2}}\bm{\theta} ^{\wedge }+\dfrac{\left\| \bm{\theta} \right\| -\sin \left\| \bm{\theta} \right\| }{\left\| \bm{\theta} \right\| ^{3}}\bm{\theta} ^{\wedge2}\\ \bold{J}_{r}^{-1}\left( \bm{\theta} \right) =I+\dfrac{1}{2}\bm{\theta} ^{\wedge }+\left( \dfrac{1}{\left\| \bm{\theta} \right\| ^{2}}-\dfrac{1+\cos \left\| \bm{\theta} \right\| }{2\left\| \bm{\theta} \right\| \sin \left\| \bm{\theta} \right\| }\right) \bm{\theta} ^{\wedge 2}
Jr(θ)=I−∥θ∥21−cos∥θ∥θ∧+∥θ∥3∥θ∥−sin∥θ∥θ∧2Jr−1(θ)=I+21θ∧+(∥θ∥21−2∥θ∥sin∥θ∥1+cos∥θ∥)θ∧2
SO(3)的右Jacobian矩阵对于任意
θ
\bm{\theta}
θ和小量
δ
θ
\delta \bm{\theta}
δθ具有以下性质:
E
x
p
(
θ
+
δ
θ
)
≈
E
x
p
(
θ
)
E
x
p
(
J
r
(
θ
)
δ
θ
)
E
x
p
(
θ
)
E
x
p
(
δ
θ
)
≈
E
x
p
(
θ
+
J
r
−
1
(
θ
)
δ
θ
)
\begin{aligned} {\rm Exp}\left( \bm{\theta} +\delta \bm{\theta} \right) &\approx {\rm Exp}\left( \bm{\theta} \right) {\rm Exp}\left( \bold{J}_{r}\left( \bm{\theta} \right) \delta \bm{\theta} \right) \\ {\rm Exp}\left( \bm{\theta} \right) {\rm Exp}\left( \delta \bm{\theta} \right) &\approx {\rm Exp}\left( \bm{\theta} +\bold{J}_{r}^{-1}\left( \bm{\theta} \right) \delta \bm{\theta} \right) \end{aligned}
Exp(θ+δθ)Exp(θ)Exp(δθ)≈Exp(θ)Exp(Jr(θ)δθ)≈Exp(θ+Jr−1(θ)δθ)
Jacobian with respect to the rotation vector
∂ ( q ⊗ p ⊗ q ∗ ) ∂ δ θ = ∂ ( R p ) ∂ δ θ = − R { θ } p ∧ J r ( θ ) \dfrac{\partial \left( \bold{q}\otimes \bold{p}\otimes \bold{q}^{\ast }\right) }{\partial \delta \bm{\theta} }=\dfrac{\partial \left( \bold{R} \bold{p}\right) }{\partial \delta \bm{\theta} }=-\bold{R}\left\{ \bm{\theta} \right\} \bold{p}^{\wedge} \bold{J}_{r}\left( \bm{\theta} \right) ∂δθ∂(q⊗p⊗q∗)=∂δθ∂(Rp)=−R{θ}p∧Jr(θ)
扰动
局部扰动
被扰动的朝向
q
~
\tilde{\bold{q}}
q~可表示为未经扰动的朝向
q
\bold{q}
q与局部小扰动
Δ
q
L
\Delta \bold{q}_{\mathcal{L}}
ΔqL的组合。由于Hamilton约定,这个局部扰动出现在复合积的右边,我们也可以给出等价的旋转矩阵形式
q
~
=
q
⊗
Δ
q
L
R
~
=
R
Δ
R
L
\tilde{\bold{q}}=\bold{q} \otimes \Delta \bold{q}_{\mathcal{L}}\qquad \tilde{\bold{R}}=\bold{R}\Delta \bold{R}_{\mathcal{L}}
q~=q⊗ΔqLR~=RΔRL
局部扰动很容易通过指数映射转换为切空间中定义的等效向量形式:
q
~
L
=
q
L
⊗
E
x
p
(
Δ
ϕ
L
)
R
~
L
=
R
L
E
x
p
(
Δ
ϕ
L
)
\tilde{\bold{q}}_{\mathcal{L}}=\bold{q}_{\mathcal{L}}\otimes {\rm Exp}\left( \Delta \bm{\phi} _{\mathcal{L}}\right) \qquad \tilde{\bold{R}}_{\mathcal{L}}=\bold{R}_{\mathcal{L}}{\rm Exp}\left( \Delta \bm{\phi} _{\mathcal{L}}\right)
q~L=qL⊗Exp(ΔϕL)R~L=RLExp(ΔϕL)
局部扰动的表示为:
Δ
ϕ
L
=
L
o
g
(
q
L
∗
⊗
q
~
L
)
=
L
o
g
(
R
L
T
⋅
R
~
L
)
\Delta \bm{\phi} _{\mathcal{L}}={\rm Log}\left(\bold{q}_{\mathcal{L}}^{\ast } \otimes \tilde{\bold{q}}_{\mathcal{L}}\right) ={\rm Log}\left( \bold{R}_{\mathcal{L}}^{T}\cdot\tilde{\bold{R}}_{\mathcal{L}} \right)
ΔϕL=Log(qL∗⊗q~L)=Log(RLT⋅R~L)
如果扰动角
Δ
ϕ
L
\Delta \bm{\phi}_{\mathcal{L}}
ΔϕL足够小,则四元数形式的扰动和旋转矩阵形式的扰动可以近似为泰勒展开直到线性项
Δ
q
L
=
E
x
p
(
Δ
ϕ
L
)
≈
1
+
1
2
Δ
ϕ
L
=
[
1
1
2
Δ
ϕ
L
]
Δ
R
L
=
E
x
p
(
Δ
ϕ
L
)
=
I
+
Δ
ϕ
L
∧
\begin{aligned} \Delta \bold{q}_{\mathcal{L}}&={\rm Exp}\left( \Delta \bm{\phi} _{\mathcal{L}}\right) \approx1+\dfrac{1}{2}\Delta \bm{\phi} _{\mathcal{L}}=\begin{bmatrix} 1 \\ \dfrac{1}{2}\Delta \bm{\phi} _{\mathcal{L}} \end{bmatrix} \\ \Delta \bold{R}_{\mathcal{L}}&={\rm Exp}\left( \Delta \bm{\phi} _{\mathcal{L}}\right) =I+\Delta \bm{\phi} _{\mathcal{L}}^{\wedge} \end{aligned}
ΔqLΔRL=Exp(ΔϕL)≈1+21ΔϕL=[121ΔϕL]=Exp(ΔϕL)=I+ΔϕL∧
因此,扰动可以在与
S
O
(
3
)
SO(3)
SO(3)流形相切的局部向量空间中指定。在这个向量空间中表示这些扰动的协方差矩阵是很方便的,即用一个3×3协方差矩阵表示。
全局扰动
全局扰动出现在复合积的左侧
q
~
G
=
E
x
p
(
Δ
ϕ
G
)
⊗
q
G
R
~
G
=
E
x
p
(
Δ
ϕ
G
)
⋅
R
G
\tilde{\bold{q}}_{\mathcal{G}}={\rm Exp}\left( \Delta \bm{\phi} _{\mathcal{G}}\right) \otimes \bold{q}_{\mathcal{G}} \qquad \tilde{\bold{R}}_{\mathcal{G}}={\rm Exp}\left( \Delta \bm{\phi} _{\mathcal{G}}\right) \cdot \bold{R}_{\mathcal{G}}
q~G=Exp(ΔϕG)⊗qGR~G=Exp(ΔϕG)⋅RG
全局扰动的表示为:
Δ
ϕ
G
=
L
o
g
(
q
~
G
⊗
q
G
∗
)
=
L
o
g
(
R
~
G
⋅
R
G
T
)
\Delta \bm{\phi} _{\mathcal{G}}={\rm Log}\left( \tilde{\bold{q}}_{\mathcal{G}}\otimes \bold{q}_{\mathcal{G}}^{\ast }\right) ={\rm Log}\left( \tilde{\bold{R}}_{\mathcal{G}}\cdot \bold{R}_{\mathcal{G}}^{T}\right)
ΔϕG=Log(q~G⊗qG∗)=Log(R~G⋅RGT)
全局扰动可以在与SO(3)流形原点处相切的向量空间中指定。
时间微分
在向量空间中表示局部扰动,我们可以很容易地得到时间导数的表达式。只要考虑
q
=
q
(
t
)
\bold{q} = \bold{q}(t)
q=q(t)为原始状态,
q
~
=
q
(
t
+
Δ
t
)
\tilde{\bold{q}} = \bold{q}(t +\Delta t)
q~=q(t+Δt)为扰动状态,并应用导数的定义
q
˙
=
lim
Δ
t
→
0
q
(
t
+
Δ
t
)
−
q
(
t
)
Δ
t
=
lim
Δ
t
→
0
q
⊗
Δ
q
L
−
q
Δ
t
=
lim
Δ
t
→
0
q
⊗
(
[
1
Δ
ϕ
L
/
2
]
−
[
1
0
]
)
Δ
t
=
lim
Δ
t
→
0
q
⊗
[
0
Δ
ϕ
L
/
2
]
Δ
t
=
1
2
q
⊗
[
0
ω
L
]
\begin{aligned} \dot{\bold{q}}&=\lim _{\Delta t\rightarrow 0}\dfrac{\bold{q}\left( t+\Delta t\right) -\bold{q}\left( t\right) }{\Delta t} \\ &=\lim _{\Delta t\rightarrow 0}\dfrac{\bold{q}\otimes \Delta \bold{q}_{\mathcal{L}}-\bold{q}}{\Delta t} \\ &=\lim _{\Delta t\rightarrow 0}\dfrac{\bold{q}\otimes \left( \begin{bmatrix} 1 \\ \Delta \bm{\phi} _{\mathcal{L}}/2 \end{bmatrix}-\begin{bmatrix} 1 \\ 0 \end{bmatrix}\right) }{\Delta t} \\ &=\lim _{\Delta t\rightarrow 0}\dfrac{\bold{q}\otimes \begin{bmatrix} 0 \\ \Delta \bm{\phi} _{\mathcal{L}}/2 \end{bmatrix}}{\Delta t} \\ &=\dfrac{1}{2}\bold{q}\otimes \begin{bmatrix} 0 \\ \bm{\omega} _{\mathcal{L}} \end{bmatrix} \end{aligned}
q˙=Δt→0limΔtq(t+Δt)−q(t)=Δt→0limΔtq⊗ΔqL−q=Δt→0limΔtq⊗([1ΔϕL/2]−[10])=Δt→0limΔtq⊗[0ΔϕL/2]=21q⊗[0ωL]
其中
Δ
ϕ
L
\Delta \bm{\phi}_{\mathcal{L}}
ΔϕL为局部扰动角,对应于由
q
\bold{q}
q定义的局部坐标系,
ω
L
\bm{\omega}_{\mathcal{L}}
ωL是对应的角度变化率
ω
L
(
t
)
=
d
ϕ
L
(
t
)
d
t
=
lim
Δ
t
→
0
Δ
ϕ
L
Δ
t
\bm{\omega} _{L}\left( t\right) =\dfrac{d\bm{\phi} _{\mathcal{L}}\left( t\right) }{dt}=\lim _{\Delta t\rightarrow 0}\dfrac{\Delta \bm{\phi} _{\mathcal{L}}}{\Delta t}
ωL(t)=dtdϕL(t)=Δt→0limΔtΔϕL
定义
Ω
(
ω
)
≜
[
ω
]
R
=
[
0
−
ω
T
ω
−
ω
∧
]
\bold{\Omega}(\bm{\omega})\triangleq \left[ \bm{\omega} \right]_R= \begin{bmatrix} 0&-\bm{\omega}^T\\ \bm{\omega} & -\bm{\omega}^{\wedge} \end{bmatrix}
Ω(ω)≜[ω]R=[0ω−ωT−ω∧]
可以得到
q
˙
=
1
2
Ω
(
ω
L
)
q
=
1
2
q
⊗
ω
L
,
R
˙
=
R
ω
L
∧
(1)
\dot{\bold{q}}=\dfrac{1}{2}\bold{\Omega}(\bm{\omega}_{\mathcal{L}})\bold{q}=\dfrac{1}{2}\bold{q}\otimes \bm{\omega}_{\mathcal{L}},\qquad \dot{\bold{R}}=\bold{R}\bm{\omega}_{\mathcal{L}}^{\wedge}\tag{1}
q˙=21Ω(ωL)q=21q⊗ωL,R˙=RωL∧(1)
这些表达式也可以在旋转群
S
O
(
3
)
SO(3)
SO(3)的框架中得出的.然而在微分框架下,我们能够清楚地将角速度
ω
L
\bm{\omega}_{\mathcal{L}}
ωL与一个特定的参考系联系起来.在上面这种情况下这个参考系是由方向
q
\bold{q}
q或
R
\bold{R}
R定义的局部坐标系.
与全局扰动相关的时间导数同理可得,其结果为
q
˙
=
1
2
ω
G
⊗
q
,
R
˙
=
ω
G
∧
R
(2)
\dot{\bold{q}}=\dfrac{1}{2}\bm{\omega}_{\mathcal{G}}\otimes \bold{q},\qquad \dot{\bold{R}}=\bm{\omega}_{\mathcal{G}}^{\wedge}\bold{R} \tag{2}
q˙=21ωG⊗q,R˙=ωG∧R(2)
其中
ω
G
(
t
)
≜
d
ϕ
G
d
t
\bm{\omega}_{\mathcal{G}}(t)\triangleq \dfrac{{\rm d}\bm{\phi}_{\mathcal{G}}}{{\rm d}t}
ωG(t)≜dtdϕG
是在全局坐标系中表示的角速度向量.
全局与局部的关系
1
2
ω
G
⊗
q
=
q
˙
=
1
2
q
⊗
ω
L
⟹
ω
G
=
q
⊗
ω
L
⊗
q
∗
=
R
ω
L
\dfrac{1}{2}\bm{\omega}_{\mathcal{G}}\otimes \bold{q}=\dot{\bold{q}}= \dfrac{1}{2}\bold{q}\otimes \bm{\omega}_{\mathcal{L}}\\ \implies \bm{\omega}_{\mathcal{G}}=\bold{q}\otimes \bm{\omega}_{\mathcal{L}}\otimes \bold{q}^{\ast}=\bold{R} \bm{\omega}_{\mathcal{L}}
21ωG⊗q=q˙=21q⊗ωL⟹ωG=q⊗ωL⊗q∗=RωL
类似的,考虑
Δ
ϕ
≈
ω
Δ
t
\Delta \bm{\phi} \approx \bm{\omega} \Delta t
Δϕ≈ωΔt
Δ
ϕ
G
=
q
⊗
Δ
ϕ
L
⊗
q
∗
=
R
Δ
ϕ
L
\Delta\bm{\phi}_{\mathcal{G}}=\bold{q}\otimes \Delta\bm{\phi}_{\mathcal{L}}\otimes \bold{q}^{\ast}=\bold{R} \Delta\bm{\phi}_{\mathcal{L}}
ΔϕG=q⊗ΔϕL⊗q∗=RΔϕL
也就是说,我们可以使用四元数或旋转矩阵将角速率向量
ω
\bm{\omega}
ω和小角扰动
Δ
ϕ
\Delta \bm{\phi}
Δϕ进行坐标系变换,就像它们是普通向量一样。
由
ω
=
u
ω
\bm{\omega} = \bold{u}\bm{\omega}
ω=uω,或
Δ
ϕ
=
u
Δ
ϕ
\Delta \bm{\phi} = \bold{u}\Delta\bm{\phi}
Δϕ=uΔϕ且旋转不改变向量长度,可以得到
u
G
=
q
⊗
u
L
⊗
q
∗
=
R
u
L
\bold{u}_{\mathcal{G}}=\bold{q}\otimes \bold{u}_{\mathcal{L}}\otimes \bold{q}^{\ast}=\bold{R} \bold{u}_{\mathcal{L}}
uG=q⊗uL⊗q∗=RuL
旋转速度的时间积分
对局部旋转速度定义的微分方程为(1),对全局旋转速率定义的微分方程为(2)。通过对旋转速度的微分方程积分来完成以四元数形式积分旋转随时间的变化。在通常的情况下,角速度由局部传感器测量,从而在离散时间
t
n
=
n
Δ
t
t_n = n\Delta t
tn=nΔt时提供局部测量
ω
(
t
n
)
\bm{\omega}(t_n)
ω(tn)。因为我们主要关注微分方程(1)。
q
˙
(
t
)
=
1
2
q
(
t
)
⊗
ω
(
t
)
\dot{\bold{q}}(t)=\dfrac{1}{2}\bold{q}(t)\otimes \bm{\omega}(t)
q˙(t)=21q(t)⊗ω(t)
定义
q
n
≜
q
(
t
n
)
,
ω
n
≜
ω
(
t
n
)
\bold{q}_n\triangleq \bold{q}(t_n),\bm{\omega}_n\triangleq \bm{\omega}(t_n)
qn≜q(tn),ωn≜ω(tn),由泰勒展开可得:
q
n
+
1
=
q
n
+
q
˙
n
Δ
t
+
1
2
!
q
¨
n
Δ
t
2
+
1
3
!
q
n
(
3
)
Δ
t
3
+
⋯
\bold{q}_{n+1}=\bold{q}_{n}+\dot{\bold{q}}_{n}\Delta t+\dfrac{1}{2!}\ddot{\bold{q}}_{n}\Delta t^{2}+\dfrac{1}{3!}{\bold{q}}^{(3)}_{n}\Delta t^{3}+\cdots
qn+1=qn+q˙nΔt+2!1q¨nΔt2+3!1qn(3)Δt3+⋯
反复运用四元数导数的表达式,假设
ω
¨
=
0
\ddot{\bm{\omega}}=0
ω¨=0,可以很容易的得到
q
n
\bold{q}_n
qn的逐次导数
q
˙
n
=
1
2
q
n
ω
n
q
¨
n
=
1
2
2
q
n
ω
n
2
+
1
2
q
n
ω
˙
q
¨
n
=
1
2
3
q
n
ω
n
3
+
1
2
2
q
n
ω
˙
ω
n
+
1
2
q
n
ω
n
ω
˙
q
n
(
i
⩾
4
)
=
1
2
i
q
n
ω
n
i
+
⋯
\begin{aligned} \dot{\bold{q}}_{n}&=\dfrac{1}{2}\bold{q}_{n}\bm{\omega}_{n}\\ \ddot{\bold{q}}_{n}&=\dfrac{1}{2^{2}}\bold{q}_{n}\bm{\omega}_{n}^{2}+\dfrac{1}{2}\bold{q}_{n}\dot{\bm{\omega} } \\ \ddot{\bold{q}}_{n}&=\dfrac{1}{2^{3}}\bold{q}_{n}\bm{\omega}_{n}^{3}+\dfrac{1}{2^2}\bold{q}_{n}\dot{\bm{\omega} }\bm{\omega}_n+\dfrac{1}{2}\bold{q}_n\bm{\omega}_n\dot{\bm{\omega}}\\ \bold{q}_n^{(i\geqslant 4)}&=\dfrac{1}{2^i}\bold{q}_n\bm{\omega}_n^i+\cdots \end{aligned}
q˙nq¨nq¨nqn(i⩾4)=21qnωn=221qnωn2+21qnω˙=231qnωn3+221qnω˙ωn+21qnωnω˙=2i1qnωni+⋯
零阶积分
前向积分
如果在
[
t
n
,
t
n
+
1
]
[t_n,t_{n+1}]
[tn,tn+1]期间是匀角速度运动,即
ω
=
0
˙
\dot{\bm{\omega}=0}
ω=0˙,则有
q
n
+
1
=
q
n
⊗
(
1
+
1
2
ω
n
Δ
t
+
1
2
!
(
1
2
ω
n
Δ
t
)
2
+
1
3
!
(
1
2
ω
n
Δ
t
)
3
+
⋯
)
\bold{q}_{n+1}=\bold{q}_{n}\otimes \left( 1+\dfrac{1}{2}\bm{\omega}_{n}\Delta t+\dfrac{1}{2!}\left( \dfrac{1}{2}\bm{\omega}_{n}\Delta t\right) ^{2}+\dfrac{1}{3!}\left( \dfrac{1}{2}\bm{\omega}_{n}\Delta t\right) ^{3}+\cdots \right)
qn+1=qn⊗(1+21ωnΔt+2!1(21ωnΔt)2+3!1(21ωnΔt)3+⋯)
这里我们很容易就可以识别出其中的泰勒展开级数就是
e
ω
n
Δ
t
/
2
e^{\bm{\omega}_n\Delta t/2}
eωnΔt/2,这个指数对应于增量旋转为
θ
=
ω
Δ
t
\bm{\theta}=\bm{\omega} \Delta t
θ=ωΔt
e
ω
Δ
t
/
2
=
ω
Δ
t
=
q
{
ω
Δ
t
}
=
[
cos
(
∥
ω
∥
Δ
t
/
2
)
ω
∥
ω
∥
sin
(
∥
ω
∥
Δ
t
/
2
)
]
e^{\bm{\omega} \Delta t/2}={\rm \bm{\omega}\Delta t}=\bold{q}\{\bm{\omega} \Delta t\}= \begin{bmatrix} \cos(\Vert\bm{\omega}\Vert\Delta t/2)\\ \frac{\bm{\omega}}{\Vert\bm{\omega}\Vert}\sin(\Vert\bm{\omega}\Vert\Delta t/2) \end{bmatrix}
eωΔt/2=ωΔt=q{ωΔt}=[cos(∥ω∥Δt/2)∥ω∥ωsin(∥ω∥Δt/2)]
因此,
q
n
+
1
=
q
n
⊗
q
{
ω
n
Δ
t
}
\bold{q}_{n+1}=\bold{q}_n\otimes \bold{q}\{\bm{\omega}_n \Delta t\}
qn+1=qn⊗q{ωnΔt}
反向积分
我们也可以认为周期
Δ
t
\Delta t
Δt上的恒定角速度为
ω
n
+
1
\bm{\omega}_{n+1}
ωn+1,即周期结束时测量的角速度。用类似的方式在将
q
n
\bold{q}_n
qn在
t
n
+
1
t_{n+1}
tn+1时刻泰勒展开,可以得到
q
n
=
q
n
+
1
+
q
˙
n
+
1
(
−
Δ
t
)
+
1
2
!
q
¨
n
+
1
(
−
Δ
t
)
2
+
1
3
!
q
n
+
1
(
3
)
(
−
Δ
t
)
3
+
⋯
\bold{q}_{n}=\bold{q}_{n+1}+\dot{\bold{q}}_{n+1}\left( -\Delta t\right) +\dfrac{1}{2!}\ddot{\bold{q}}_{n+1}\left( -\Delta t\right) ^{2}+\dfrac{1}{3!}\bold{q}_{n+1}^{\left( 3\right) }\left( -\Delta t\right) ^{3}+\cdots
qn=qn+1+q˙n+1(−Δt)+2!1q¨n+1(−Δt)2+3!1qn+1(3)(−Δt)3+⋯
类似的,假设
ω
˙
=
0
\dot{\bm{\omega}}=0
ω˙=0
q
n
+
1
(
i
)
=
q
n
+
1
⊗
(
1
2
ω
n
+
1
)
i
{\bold{q}}^{(i)}_{n+1}=\bold{q}_{n+1}\otimes(\frac{1}{2}\bm{\omega}_{n+1})^i
qn+1(i)=qn+1⊗(21ωn+1)i
于是
q
n
=
q
n
+
1
⊗
[
∑
k
=
0
∞
1
k
!
(
−
1
2
ω
n
+
1
Δ
t
)
k
]
=
q
n
+
1
⊗
E
x
p
(
−
ω
n
+
1
Δ
t
)
⟹
q
n
+
1
=
q
n
⊗
E
x
p
(
ω
n
+
1
Δ
t
)
\begin{aligned} \bold{q}_n&=\bold{q}_{n+1}\otimes\left[\sum_{k=0}^\infty \frac{1}{k!}\left(-\frac{1}{2}\bm{\omega}_{n+1}\Delta t\right)^k\right]\\ &=\bold{q}_{n+1}\otimes {\rm Exp}(-\bm{\omega}_{n+1}\Delta t) \\ \implies \bold{q}_{n+1}&=\bold{q}_n\otimes {\rm Exp}(\bm{\omega}_{n+1}\Delta t) \end{aligned}
qn⟹qn+1=qn+1⊗[k=0∑∞k!1(−21ωn+1Δt)k]=qn+1⊗Exp(−ωn+1Δt)=qn⊗Exp(ωn+1Δt)
所以最终可得
q
n
+
1
≈
q
n
⊗
q
{
ω
n
+
1
Δ
t
}
\bold{q}_{n+1}\approx \bold{q}_n\otimes \bold{q}\{\bm{\omega}_{n+1}\Delta t\}
qn+1≈qn⊗q{ωn+1Δt}
中点积分
类似的,认为周期
Δ
t
\Delta t
Δt上的恒定角速度为中间角速度(不一定是周期中点角速度)
ω
ˉ
=
ω
n
+
ω
n
+
1
2
\bar{\bm{\omega}}=\frac{\bm{\omega}_{n}+\bm{\omega}_{n+1}}{2}
ωˉ=2ωn+ωn+1
则有
q
n
+
1
≈
q
n
⊗
q
{
ω
ˉ
Δ
t
}
\bold{q}_{n+1}\approx \bold{q}_n\otimes \bold{q}\{\bar{\bm{\omega}}\Delta t\}
qn+1≈qn⊗q{ωˉΔt}