1.1 参考系的旋转
我们使用下标表示特定参考系下的一个向量。惯性参考系(inertial frame)下的一个向量表示为 x ⃗ \vec x x,在本体参考系(body-reference)表示为 x ⃗ b \vec x_b xb。以此同时,我们假设上面这些坐标系的原点一致的(一样的),而且本体坐标系有不同的角度方向。该角度方向有几个众所周知的描述,包括欧拉角和欧拉参数(四元数)。前一种方法(指欧拉角)涉及围绕主轴的连续旋转,并且具有滚动(roll),俯仰(pitch)和偏航(yaw)的直观概念的牢固联系。欧拉角的问题之一是对于某些特定值,变换表现出不连续性。四元数呈现出更优雅和更健壮的方法,但具有更多的抽象。我们将使用欧拉角开发(develop)运动方程。
将三支铅笔粘在一起,形成右手三维坐标系。连续地围绕其自身的三个主轴旋转系统,很容易看出可以实现到达任何可能的方向。例如,考虑一下【yaw, pitch, roll】的顺序旋转:从与某些惯性坐标系一个方向开始,以yaw为轴旋转可移动的坐标系,然后绕着pitch轴旋转移动坐标系,接着以roll为轴旋转可移动坐标系(这里需要注意的是每次绕轴旋转旋都是以旋转后的可移动坐标系基础上旋转)。不用说,有许多有效的欧拉角旋转组可以达到给定的方向; 其中一些可能会使用相同的轴两次。
第一个问题是:在惯性空间中固定的点的坐标在可旋转的本体坐标系如何表示? 该变换使用3x3的矩阵表示,现在我们通过三个欧拉角的连续旋转推导。在第一次旋转之前,本体参考坐标与惯性系的坐标匹配(也就是开始旋转之前两个坐标系重合),即有
x
⃗
b
=
x
⃗
\vec x_b = \vec x
xb=x。现在让我们以yaw为轴(z轴)旋转可移动坐标系
ϕ
\phi
ϕ角度,我们有
x
⃗
b
1
=
[
cos
ϕ
sin
ϕ
0
−
sin
ϕ
cos
ϕ
0
0
0
1
]
x
⃗
b
0
=
R
(
ϕ
)
x
⃗
b
0
(1)
\vec x_b^1 = \begin{bmatrix} \cos \phi & \sin \phi & 0 \\ -\sin \phi & \cos \phi & 0 \\ 0 & 0 & 1 \end{bmatrix} \vec x^0_b = R(\phi)\vec x^0_b \tag{1}
xb1=⎣⎡cosϕ−sinϕ0sinϕcosϕ0001⎦⎤xb0=R(ϕ)xb0(1)
绕着z轴旋转不会改变该点的z坐标。其他轴根据基本三角学进行修改即可。现在应用第二次旋转,绕y轴(pitch)旋转
θ
\theta
θ角度
x
⃗
b
2
=
[
cos
θ
0
−
sin
θ
0
1
0
sin
θ
0
cos
θ
]
x
⃗
b
1
=
R
(
θ
)
x
⃗
b
1
(2)
\vec x_b^2 = \begin{bmatrix} \cos \theta & 0 & -\sin \theta \\ 0 & 1 & 0 \\ \sin \theta & 0 & \cos \theta \end{bmatrix} \vec x^1_b = R(\theta)\vec x^1_b \tag{2}
xb2=⎣⎡cosθ0sinθ010−sinθ0cosθ⎦⎤xb1=R(θ)xb1(2)
最后本体坐标系绕着x轴(roll)旋转
ψ
\psi
ψ角度
x
⃗
b
3
=
[
1
0
0
0
cos
ψ
sin
ψ
0
−
sin
ψ
cos
ψ
]
x
⃗
b
2
=
R
(
ψ
)
x
⃗
b
2
(3)
\vec x_b^3 = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos \psi & \sin \psi \\ 0 & -\sin \psi & \cos \psi \end{bmatrix} \vec x^2_b = R(\psi)\vec x^2_b \tag{3}
xb3=⎣⎡1000cosψ−sinψ0sinψcosψ⎦⎤xb2=R(ψ)xb2(3)
x
⃗
b
3
\vec x^3_b
xb3表示在完全变换本体坐标系下原始点在本体坐标系下的位置。这里,我们将要使用
x
⃗
b
\vec x_b
xb代替
x
⃗
b
3
\vec x^3_b
xb3来表示惯性坐标系下的一个点在本体坐标系下的表示。这三个无关的旋转可以通过矩阵相乘叠加在一起:
x
⃗
b
=
R
(
ψ
)
R
(
θ
)
R
(
ϕ
)
(4)
\vec x_b = R(\psi)R(\theta)R(\phi) \tag{4}
xb=R(ψ)R(θ)R(ϕ)(4)
=
[
c
θ
c
ϕ
c
θ
s
ϕ
−
s
θ
−
c
ψ
s
ϕ
+
s
ψ
s
θ
c
ϕ
c
ψ
c
ϕ
+
s
ψ
s
θ
s
ϕ
s
ψ
c
θ
s
ψ
s
ϕ
+
c
ψ
s
θ
c
ϕ
−
s
ψ
c
ϕ
+
c
ψ
s
θ
s
ϕ
c
ψ
c
θ
]
x
⃗
=
R
(
ϕ
,
θ
,
ψ
)
x
⃗
=\begin{bmatrix} c \theta c \phi & c \theta s \phi & - s \theta \\ -c \psi s \phi + s \psi s \theta c \phi & c \psi c \phi + s \psi s \theta s \phi & s \psi c \theta \\ s \psi s \phi + c \psi s \theta c \phi & - s \psi c \phi + c \psi s \theta s \phi & c \psi c \theta \end{bmatrix} \vec x \\ = R(\phi, \theta, \psi) \vec x
=⎣⎡cθcϕ−cψsϕ+sψsθcϕsψsϕ+cψsθcϕcθsϕcψcϕ+sψsθsϕ−sψcϕ+cψsθsϕ−sθsψcθcψcθ⎦⎤x=R(ϕ,θ,ψ)x
所有的变换矩阵(
R
(
ϕ
)
,
R
(
θ
)
,
R
(
ψ
)
R(\phi), R(\theta), R(\psi)
R(ϕ),R(θ),R(ψ),包括
R
(
ϕ
,
θ
,
ψ
)
\R(\phi, \theta, \psi)
R(ϕ,θ,ψ)都是正交矩阵:他们的逆等于他们的转秩。另外,我们应该注意到旋转矩阵
R
R
R对于所有方向的表示都是通用的,包括四元数。如所写的,三角函数以我们执行旋转的顺序作用于欧拉角。
在可移动(本体)参考系具有与惯性系不同的原点的情况下,我们有
x
⃗
=
x
⃗
0
+
R
T
x
⃗
b
(5)
\vec x = \vec x_0 + R^T \vec x_b \tag{5}
x=x0+RTxb(5)
x
⃗
0
\vec x_0
x0是以惯性坐标表示移动原点的位置。
1.2 微分旋转
现在考虑从一个坐标系系到另一个坐标系的旋转很小的角度,假设忽略高阶项给出
R
≃
[
1
δ
ϕ
−
δ
θ
−
δ
ϕ
1
δ
ψ
δ
θ
−
δ
ψ
1
]
=
[
0
δ
ϕ
−
δ
θ
−
δ
ϕ
0
δ
ψ
δ
θ
−
δ
ψ
0
]
+
I
3
×
3
(6)
R \simeq \begin{bmatrix} 1 & \delta \phi & - \delta \theta \\ -\delta \phi & 1 & \delta \psi \\ \delta \theta & -\delta \psi & 1 \end{bmatrix} \\ \tag{6} \\ =\begin{bmatrix} 0 & \delta \phi & - \delta \theta \\ -\delta \phi & 0 & \delta \psi \\ \delta \theta & -\delta \psi & 0 \end{bmatrix} + I_{3 \times 3}
R≃⎣⎡1−δϕδθδϕ1−δψ−δθδψ1⎦⎤=⎣⎡0−δϕδθδϕ0−δψ−δθδψ0⎦⎤+I3×3(6)
其中
I
3
×
3
I_{3 \times 3}
I3×3表示单位阵,
R
R
R有单位阵和
负
的
叉
乘
[
−
δ
E
⃗
×
]
负的叉乘[-\delta \vec E_{\times}]
负的叉乘[−δE×]组成,其中
δ
E
⃗
=
[
δ
ψ
,
δ
θ
,
δ
ϕ
]
\delta \vec E = [\delta \psi, \delta \theta, \delta \phi]
δE=[δψ,δθ,δϕ],欧拉角向量的排列顺序是轴
[
x
,
y
,
z
]
[x, y, z]
[x,y,z],对于很小的旋转可以完全解耦,也就是与排列顺序无关。由于
R
−
1
=
R
T
R^{-1} = R^T
R−1=RT,我们也有
R
−
1
=
I
3
×
3
+
δ
E
⃗
×
R^{-1} = I_{3 \times 3} + \delta \vec E_{\times}
R−1=I3×3+δE×,我们还有
x
⃗
b
=
x
⃗
−
δ
E
⃗
×
x
⃗
(7)
\vec x_b = \vec x - \delta \vec E \times \vec x \tag{7}
xb=x−δE×x(7)
x
⃗
=
x
⃗
b
+
δ
E
⃗
×
x
⃗
b
(8)
\vec x = \vec x_b + \delta \vec E \times \vec x_b \tag{8}
x=xb+δE×xb(8)
现在我们在本体坐标系中固定一个兴趣点,而不是在惯性系中,在本体坐标系中的位置用
r
⃗
\vec r
r表示(也称作半径)。微分旋转出现在很小的
δ
t
\delta t
δt的时间上,因此我们可以在该点旋转之前和之后写下它的位置,相对于第一个坐标系有:
x
⃗
(
t
)
=
r
⃗
x
⃗
(
t
+
δ
t
)
=
R
T
r
⃗
=
r
⃗
+
δ
E
⃗
×
r
⃗
(9)
\vec x(t) = \vec r \\ \vec x(t + \delta t) = R^T \vec r = \vec r + \delta \vec E \times \vec r \tag{9}
x(t)=rx(t+δt)=RTr=r+δE×r(9)
对
δ
t
\delta t
δt求导(即除以微分时间步长)有
δ
x
⃗
δ
t
=
δ
E
⃗
δ
t
×
r
⃗
=
ω
⃗
×
r
⃗
(10)
\frac{\delta \vec x}{\delta t} = \frac{\delta \vec E}{\delta t} \times \vec r \\ =\vec \omega \times \vec r \tag{10}
δtδx=δtδE×r=ω×r(10)
其中旋转速率向量
ω
≃
d
E
⃗
/
d
t
\omega \simeq d \vec E/ dt
ω≃dE/dt,因为这个无穷小旋转的欧拉角很小并且是解耦的。相同的叉积关系,相对于第二个坐标系有
x
⃗
b
(
t
)
=
R
r
⃗
=
r
⃗
−
δ
E
⃗
×
r
x
⃗
b
(
t
+
δ
t
)
=
r
⃗
(11)
\vec x_b(t) = R \vec r = \vec r - \delta \vec E \times r \\ \vec x_b(t + \delta t) = \vec r \tag{11}
xb(t)=Rr=r−δE×rxb(t+δt)=r(11)
可以得到
δ
x
⃗
b
δ
t
=
δ
E
⃗
δ
t
×
r
⃗
=
ω
⃗
×
r
⃗
(12)
\frac{\delta \vec x_b}{\delta t} = \frac{\delta \vec E}{\delta t} \times \vec r \\ =\vec \omega \times \vec r \tag{12}
δtδxb=δtδE×r=ω×r(12)
在原点被固定的旋转体中,恒定半径矢量的时间变化率是转速矢量
ω
⃗
\vec \omega
ω和半径矢量
r
⃗
\vec r
r本身的交叉积。得到的导数在移动本体坐标系中,在这种情况下,半径向量相对于本体坐标系的变化,我们需要增加一项:
d
x
⃗
b
d
t
=
δ
E
⃗
δ
t
×
r
⃗
+
∂
r
⃗
∂
t
=
ω
⃗
×
r
⃗
+
∂
r
⃗
∂
t
(13)
\frac{d\vec x_b}{d t} = \frac{\delta \vec E}{\delta t} \times \vec r + \frac{\partial \vec r}{\partial t}\\ =\vec \omega \times \vec r + \frac{\partial \vec r}{\partial t} \tag{13}
dtdxb=δtδE×r+∂t∂r=ω×r+∂t∂r(13)
最后,让原点也可以移动,有
d
x
⃗
b
d
t
=
ω
⃗
×
r
⃗
+
∂
r
⃗
∂
t
+
d
x
⃗
o
d
t
(14)
\frac{d \vec x_b}{dt} = \vec \omega \times \vec r + \frac{\partial \vec r}{\partial t} + \frac{d \vec x_o}{dt} \tag{14}
dtdxb=ω×r+∂t∂r+dtdxo(14)该项在本体参考系中经常被写作速度
v
⃗
\vec v
v:
v
⃗
=
ω
⃗
×
r
⃗
+
∂
r
⃗
∂
t
+
v
⃗
o
(15)
\vec v = \vec \omega \times \vec r + \frac{\partial \vec r}{\partial t} + \vec v_o \tag{15}
v=ω×r+∂t∂r+vo(15)
其中
v
⃗
o
\vec v_o
vo是原点的本体参考系速度。粒子的总速度等于参考帧原点的速度,加上由于该系的旋转而产生的分量。速度方程可以推广到任何本体参考矢量
f
⃗
\vec f
f:
d
f
⃗
d
t
=
∂
f
∂
t
+
ω
⃗
×
f
⃗
(16)
\frac{d \vec f}{d t} = \frac{\partial f}{\partial t} + \vec \omega \times \vec f \tag{16}
dtdf=∂t∂f+ω×f(16)
1.3 欧拉角的变化率
仅对于无穷小欧拉角的情况,欧拉角的时间变化率等于本体参考系旋转速率是正确的。例如,对于【yaw, pitchm=, roll】序列,定义的欧拉yaw角(首先旋转的那个轴)绝对不是最终本体的yaw轴,接下来的pitch和roll将要移动该轴。任何模拟的重要部分是对欧拉角的估计。因为物理确定的旋转率 ω ⃗ \vec \omega ω,我们找到一个映射 ω ⃗ → d E ⃗ / d t \vec \omega \to d\vec E / dt ω→dE/dt。
该想法意在考虑欧拉角的一个非常小的变化并确定该旋转向量的影响。第一个欧拉角暗含两个额外的旋转,第二个暗含一个旋转,最后一个没有额外的旋转:
ω
⃗
=
R
(
ψ
)
R
(
θ
)
{
0
0
d
ϕ
d
t
}
+
R
(
ψ
)
{
0
θ
d
t
0
}
+
{
d
ψ
d
t
0
0
}
(17)
\vec \omega = R(\psi)R(\theta) \begin{Bmatrix} 0 \\ 0 \\ \frac{d\phi}{dt} \end{Bmatrix} + R(\psi) \begin{Bmatrix} 0 \\ \frac{\theta}{dt} \\ 0 \end{Bmatrix} + \begin{Bmatrix} \frac{d\psi}{dt} \\ 0 \\ 0 \end{Bmatrix} \tag{17}
ω=R(ψ)R(θ)⎩⎨⎧00dtdϕ⎭⎬⎫+R(ψ)⎩⎨⎧0dtθ0⎭⎬⎫+⎩⎨⎧dtdψ00⎭⎬⎫(17)
=
[
1
0
−
sin
θ
0
cos
ψ
s
i
n
ψ
cos
θ
0
−
s
i
n
ψ
cos
ψ
cos
θ
]
{
d
ψ
d
t
d
θ
d
t
d
ϕ
d
t
}
=\begin{bmatrix} 1 & 0 & -\sin \theta \\ 0 & \cos \psi & sin \psi \cos \theta \\ 0 & -sin \psi & \cos \psi \cos \theta \end{bmatrix} \begin{Bmatrix} \frac{d\psi}{dt} \\ \frac{d\theta}{dt} \\ \frac{d\phi}{dt} \end{Bmatrix}
=⎣⎡1000cosψ−sinψ−sinθsinψcosθcosψcosθ⎦⎤⎩⎨⎧dtdψdtdθdtdϕ⎭⎬⎫
取逆有
E
⃗
d
t
=
[
1
sin
ψ
tan
θ
cos
ψ
tan
θ
0
cos
ψ
−
sin
ψ
0
sin
ψ
/
cos
θ
cos
ψ
/
cos
θ
]
ω
⃗
=
Γ
(
E
⃗
)
ω
⃗
(18)
\frac{\vec E}{dt} = \begin{bmatrix} 1 & \sin \psi \tan \theta & \cos \psi \tan \theta \\ 0 & \cos \psi & - \sin \psi \\ 0 & \sin \psi/\cos \theta & \cos \psi / \cos \theta \end{bmatrix} \vec \omega = \Gamma (\vec E) \vec \omega \tag{18}
dtE=⎣⎡100sinψtanθcosψsinψ/cosθcosψtanθ−sinψcosψ/cosθ⎦⎤ω=Γ(E)ω(18)
Γ
\Gamma
Γ在
θ
=
{
π
/
2
,
3
π
/
2
}
\theta = \{\pi/2, 3\pi/2\}
θ={π/2,3π/2}存在奇异性,因为除以
cos
θ
\cos \theta
cosθ,当车绕着y轴旋转90时,对于传播本体角度的方向来说,等式不成立(因为当
θ
=
9
0
∘
\theta = 90 ^\circ
θ=90∘时,
cos
θ
=
0
\cos \theta = 0
cosθ=0),在实际应用中是可能出现的,例如,在轨道卫星和机器人手臂中,四元数提供无缝映射。 对于大多数海洋船只来说,奇点是可以接受的,只要它不在偏航轴上!
reference
[1]MIT. Kinematics Of Moving Frames. url: https://ocw.mit.edu/courses/mechanical-engineering/2-154-maneuvering-and-control-of-surface-and-underwater-vehicles-13-49-fall-2004/lecture-notes/lec1.pdf (visited on 2004).