四元数与旋转——学习笔记(三)

来源：互联网发布：手机免费视频软件编辑：程序博客网时间：2024/06/09 22:39

联系方式：860122112@qq.com

四元数与旋转——学习笔记(一)
四元数与旋转——学习笔记(二)
四元数与旋转——学习笔记(三)

五、四元数的微分和积分

假设前面公式(1)给出的四元数是时间 t 的函数。则对 q(t) 求导有：

q = q 0 ˙ + q ˙ = q 0 ˙ + q 1 ˙ i^+ q 2 ˙ j^+ q 3 ˙ k^

设

p(t) 是另一个四元数，则对前面公式(3)关于时间

t 求导：

d d t (p q) = p ˙ q + p q ˙

四元数函数 q(t) 描述的是某移动物体方向的变化，具体一点就是物体的自身坐标系相对于固定的世界坐标系的变化（是坐标系旋转！）。令 ω(t) 是物体坐标系相对于世界坐标系变化的角速度。角速度可以由牛顿动力学方程确定。

定理3
令 q(t) 是一个单位四元数函数，ω(t) 是由 q(t) 确定的角速度。则 q(t) 的导数为：
(14)

q ˙ = 1 2 ω q

证明
在

t+Δt 时刻，旋转可以描述为

q(t+Δt)。在

Δt 过程中，物体坐标系在经过了

q(t) 旋转的前提下，又经过了额外的微小旋转。这个额外的微小旋转的瞬时旋转轴为

ω^=ω/∥ω∥，旋转角度为

Δθ=∥ω∥Δt，可以用一个单位四元数描述：

Δ q = = cos Δ θ 2 + ω^sin Δ θ 2 cos ∥ ω ∥ Δ t 2 + ω^sin ∥ ω ∥ Δ t 2

t+Δt 时刻的旋转可以描述成四元数序列

q(t),Δq，由第四节的内容有：

q (t + Δ t) = Δ q q (t)

（注意先旋转的四元数在右边。有童鞋可能会疑惑这是坐标系旋转，为什么不是左边？因为这里的坐标系旋转其实就是将物体坐标系各个基向量进行旋转，与前面提到的指定向量不动坐标系旋转不同）
接下来推导

q˙，首先求差值：

q (t + Δ t) - q (t) = = (cos ∥ ω ∥ Δ t 2 + ω^sin ∥ ω ∥ Δ t 2) q (t) - q (t) - 2 sin 2 ∥ ω ∥ Δ t 2 q (t) + ω^sin ∥ ω ∥ Δ t 2 q (t)

等式右边第一项是高阶项趋近于0，可省略。因此

q ˙ (t) = = = = = lim Δ t \to 0 q ( t + Δ t ) - q ( t ) Δ t ω^lim Δ t \to 0 sin ( ∥ ω ∥ Δ t / 2 ) Δ t q (t) ω^d d t sin (∥ ω ∥ t 2) | t = 0 q (t) ω^∥ ω ∥ 2 q (t) 1 2 ω (t) q (t)

在实际中，能得到的角速度其实是相对于物体坐标系下的角速度

ω′（例如通过IMU获得）。世界坐标系下的角速度

ω 不动，物体坐标系旋转，则角速度

ω 在物体坐标系下的表示为

ω′=q∗ωq。将

ω=qω′q∗ 代入(14)，得到：

q ˙ = 1 2 q ω'

若

q˙ 已知，(14)两边同时右乘

q∗，可以求得角速度：

ω = 2 q ˙ q *

根据(14)，四元数的二阶导为：

q ¨ = = = 1 2 (ω ˙ q + ω q ˙) 1 2 ω ˙ q + 1 4 ω ω q (- 1 4 ∥ ω ∥ 2 + 1 2 ω ˙) q

当 q 的一阶导和二阶导都是已知时，由上式两边同时右乘 q∗ 可以求出角加速度：

ω ˙ = = = 2 q ¨ q * - ω q ˙ q * 2 q ¨ q * - 2 q ˙ q * q ˙ q * 2 (q ¨ q * - (q ˙ q *) 2)

微分方程(14)可以用数值积分求解。令

h 为时间步长，

qk 和

ωk 是

kh 时刻的四元数和角速度。则在下一时刻的四元数近似为：

q k + 1 = q k + 1 2 h ω k q k

由于增量计算，

qk+1不一定是单位四元数，所以需要规范化为单位四元数：

q k + 1 \leftarrow q k + 1 ∥ q k + 1 ∥

六、四元数插值

四元数插值意思就是给定物体起始方位和终点方位，需要生成从起始到终点的平滑旋转运动（图形和动画经常用到）。令起始和最终的旋转表示为以下两个单位四元数：

r 1 = cos θ 1 2 + u^1 sin θ 1 2

r 2 = cos θ 2 2 + u^2 sin θ 2 2

要使插值有意义，

r1≠r2 必须成立。

1、旋转轴和角度的恒定变化率
旋转角的线性插值为
(15)

θ (τ) = (1 - τ) θ 1 + τ θ 2

且

τ∈[0,1]。但是旋转轴的线性插值

v=(1−τ)u^1+τu^2 得到的

v 不是单位四元数，规范化后得

w=v/∥v∥，

w(τ) 的曲线不是以恒定的速率变化，这会造成由

u^1到

u^2旋转时不平稳。

由于 u^1 和 u^2 可看做都处于一个单位球面上，则可以在它们球面的最短路径（即最短的弧）上做插值。如下图：
这里写图片描述
令一个点 u^(τ) 以恒定速率在最短的弧上从 u^1 到 u^2，即 τ由 0 到 1 增长。本质上， u^(τ) 是弧在[0,1]上恒定速率的参数化。要推导 u^(τ)，首先构造 u^1 和 u^2 所在平面的法向量：

n^= u ^ 1 \times u ^ 2 ∥ u ^ 1 \times u ^ 2 ∥

n^ 确定后，关于

n^ 由

u^1 旋转到

u^2 的角度为

ϕ∈(0,π]，且

ϕ = arccos (u^1 \cdot u^2)

因此

u^(τ) 由

u^1 绕轴

n^ 旋转角度

τϕ 后确定，该四元数运算为：

u^(τ) = (cos τ ϕ 2 + n^sin τ ϕ 2) u^1 (cos τ ϕ 2 - n^sin τ ϕ 2)

事实上，

u^ 有不包含四元数乘法的更简洁形式：
(16)

u^(τ) = sin ( ( 1 - τ ) ϕ ) sin ϕ u^1 + sin ( τ ϕ ) ϕ u^2

证明上式成立，首先考虑

u^1 和

u^2 在

xy 平面上，且

n^ 沿

z 轴方向。令

u^1=(cosθ1,sinθ1) ，

u^2=(cos(θ1+ϕ),sin(θ1+ϕ))，则：

= = = = sin ( ( 1 - τ ) ϕ ) sin ϕ u^1 + sin ( τ ϕ ) ϕ u^2 (sin ( ( 1 - τ ) ϕ ) cos θ 1 + sin ( τ ϕ ) cos ( θ 1 + ϕ ) sin ϕ, sin ( ( 1 - τ ) ϕ ) sin θ 1 + sin ( τ ϕ ) sin ( θ 1 + ϕ ) sin ϕ) (sin ϕ ( cos ( τ ϕ ) cos θ 1 - sin ( τ ϕ ) sin θ 1 ) sin ϕ, sin ϕ ( cos ( τ ϕ ) sin θ 1 + sin ( τ ϕ ) cos θ 1 ) sin ϕ) (cos (θ 1 + τ ϕ), sin (θ 1 + τ ϕ) u^(τ)

若

u^1 和

u^2 不在

xy 平面上，先通过旋转让

n^ 和

z 轴方向一致。令

R 为对应的旋转矩阵，则有：

u^(τ) = = = R - 1 (R u^) R - 1 (sin ( ( 1 - τ ) ϕ ) sin ϕ (R u^1) + sin ( τ ϕ ) ϕ (R u^2)) sin ( ( 1 - τ ) ϕ ) sin ϕ u^1 + sin ( τ ϕ ) ϕ u^2

最后，由(15)(16)得到起始和最终单位四元数

r1,r2 的插值公式：

r (τ) = = cos θ ( τ ) 2 + u^1 (τ) sin θ ( τ ) 2 cos ( 1 - τ ) θ 1 + τ θ 2 2 + (sin ( ( 1 - τ ) ϕ ) sin ϕ u^1 + sin ( τ ϕ ) ϕ u^2) sin ( 1 - τ ) θ 1 + τ θ 2 2

此插值公式满足旋转角和旋转轴变化速率恒定。

2、球面线性插值
球面线性插值Slerp(spherical linear interpolation)的形式：

r (τ) = r 1 (r * 1 r 2) τ

τ∈[0,1]，单位四元数的幂由前面给出的公式(6)得到。

七、关于四元数的微分

对关于四元数 q=q0+q 求微分，其中 q=(q1,q2,q3) 。可以把四元数看成一个列向量 (q0,q1,q2,q3)T 。
给定一个向量 u=(u1,u2,u3)T ，记向量叉乘 u× 是 3×3 的反对称阵：

u \times = ⎡ ⎣ ⎢ 0 u 3 - u 2 - u 3 0 u 1 u 2 - u 1 0 ⎤ ⎦ ⎥

此外，有：

\partial ( u \times v ) \partial v = \partial ( u \times v ) \partial u = = u \times - \partial ( v \times u ) \partial u - v \times

给定两个四元数

p=p0+p 和

q=q0+q ，列向量形式分别是

(p0,p1,p2,p3)T 和

(q0,q1,q2,q3)T ，则以下的求微分得到

4×4 矩阵：

\partial ( p q ) \partial p = = \partial \partial p [p 0 q 0 - p \cdot q p 0 q + q 0 p + p \times q] [q 0 q - q T q 0 I 3 \times 3 - q \times]

\partial ( p q * ) \partial p = [q 0 - q q T q 0 I 3 \times 3 + q \times]

\partial ( p q ) \partial q = [p 0 p - p T p 0 I 3 \times 3 + p \times]

\partial ( p q * ) \partial q = [p 0 p p T - p 0 I 3 \times 3 - p \times]

令 v 是纯四元数，则有：

\partial ( q v ) \partial q = [0 v - v T - v \times]

\partial ( v q ) \partial q = [0 v - v T v \times]

\partial ( q v ) \partial v = [- q T q 0 I 3 \times 3 + q \times]

\partial ( v q ) \partial v = [- q T q 0 I 3 \times 3 - q \times]

前面的公式(10)关于 v 的微分为：

\partial ( q v q * ) \partial v = (q 20 - ∥ q ∥ 2) I 3 \times 3 + 2 q q T + 2 q 0 q \times

因为

\partial ( ∥ q ∥ 2 v ) \partial q = = \partial ( v q T q ) \partial q 2 v q T

\partial ( ( q \cdot v ) q ) \partial q = = \partial ( ( v T q ) q ) \partial q v T q I 3 \times 3 + q v T

所以公式(10)关于四元数 q 的微分为：

\partial ( q v q * ) \partial q = 2 (q 0 v + q \times v, - v q T + (v \cdot q) I 3 \times 3 + q v T - q 0 v \times)

参考文献
Yan-Bin Jia. Quaternion and Rotation

阅读全文

0 0