前言

四元数quaternion,一个熟悉又陌生的东西,我曾经短暂地接触过,但是又倏忽地将它抛弃。如今,我再次拜倒在它的脚下,只因为我现在不得不使用它,希望它能够不计前嫌 ,留我一条生路。

四元数基础

参考自Understanding Quaternions,用于记录学习四元数的过程。

定义

四元数是最简单的一种超复数,可以看作是复数在三维空间中的推广,其共有四个分量,定义为:

q=a+bi+cj+dkq=a+b\vec{i}+c\vec{j}+d\vec{k}

其中aa称为实部,b,c,db,c,d称为虚部,i,j,k\vec{i},\vec{j},\vec{k}是四元数的基本单位,满足:

i2=j2=k2=ijk=1\vec{i}^2=\vec{j}^2=\vec{k}^2=\vec{i}\cdot\vec{j}\cdot\vec{k}=-1\\

即三个基本单位是完备的正交基。特别的,在复数域下具有重要意义的欧拉公式,在四元数空间中也有以下的形式:

eiθ=cosθ2+sinθ2ni=q0+qie^{i\theta}=\cos\frac{\theta}{2}+\sin\frac{\theta}{2}\vec{n}\cdot\vec{i}=q_0+\vec{q}\cdot\vec{i}

四元数可以看作是超球面上的一点,即a2+b2+c2+d2=1a^2+b^2+c^2+d^2=1,而纯四元数(类比于纯虚数),是在a=0a=0的子空间中的一点。

共轭运算

四元数的共轭运算类似于复数,即将全部的虚部都取反,得到是共轭四元数:

q=abicjdkq^*=a-b\vec{i}-c\vec{j}-d\vec{k}

根据共轭运算可以定义四元数的模长,与复数相同:

q=qq=a2+b2+c2+d2|q|=\sqrt{q^*q}=\sqrt{a^2+b^2+c^2+d^2}

特别地,可以定义四元数的逆,称满足下式的另一个目标四元数为当前四元数的逆:

qq1=q1q=1q1=qq2\begin{aligned} qq^{-1}&=q^{-1}q=1\\\\ q^{-1}&=\frac{q^*}{|q|^2} \end{aligned}

四则运算

四元数的运行遵从与复数相似的法则,满足线性运算,四元数的乘除运算也与复数保持一致,即有:

q1q2=(a1a2b1b2c1c2d1d2)+(a1b2+b1a2+c1d2d1c2)i+(a1c2b1d2+c1a2+d1b2)j+(a1d2+b1c2c1b2+d1a2)k=(a1b1q1q2,a1q2+b1q1+q1×q2)\begin{aligned} q_1q_2&=(a_1a_2-b_1b_2-c_1c_2-d_1d_2)+(a_1b_2+b_1a_2+c_1d_2-d_1c_2)\vec{i}\\ &+(a_1c_2-b_1d_2+c_1a_2+d_1b_2)\vec{j}+(a_1d_2+b_1c_2-c_1b_2+d_1a_2)\vec{k}\\ &=(a_1b_1-\vec{q}_1\cdot\vec{q}_2,a_1\vec{q}_2+b_1\vec{q}_1+\vec{q}_1\times\vec{q}_2) \end{aligned}

其中最后一行所写的公式是乘法运算的快捷公式。

四元数的导数

将四元数看作是一个函数,当其虚部实部满足导数的条件时,可以认为它也存在相应的导数。仍然以q=a+bi+cj+dkq=a+b\vec{i}+c\vec{j}+d\vec{k}为例,不妨设qq是一个单位四元数,则就有:

a2+b2+c2+d2=1a^2+b^2+c^2+d^2=1

对其求导,得到:

aa˙+bb˙+cc˙+dd˙=0a\dot{a}+b\dot{b}+c\dot{c}+d\dot{d}=0

单位四元数多用于计算旋转,因此本式实际指出了一种最基本的动力学方程形式。

四元数与旋转

四元数现在基本上只用于计算坐标系旋转变换,与其对应 的另一个著名表达式是欧拉轴角公式。以书中作者给出的例子为参考,记偏航角为ψ\psi,俯仰角为θ\theta,滚转角为γ\gamma,分别对应绕Z0Z_0轴,YY'轴和X1X_1轴旋转,如下图所示:
欧拉角的坐标系旋转

三轴欧拉角

绕三个轴的欧拉角旋转矩阵可以写为:

R=RZ(ψ)RY(θ)RX(γ)R=R_Z(\psi)R_Y(\theta)R_X(\gamma)

如以逆时针转角为正,则可以求得其中绕Z0Z_0轴旋转的矩阵为:

RZ(ψ)=[cosψsinψ0sinψcosψ0001]R_Z(\psi)=\begin{bmatrix} \cos\psi&\sin\psi&0\\ -\sin\psi&\cos\psi&0\\ 0&0&1 \end{bmatrix}

YY'轴旋转的矩阵为:

RY(θ)=[cosθ0sinθ010sinθ0cosθ]R_Y(\theta)=\begin{bmatrix} \cos\theta&0&-\sin\theta\\ 0&1&0\\ \sin\theta&0&\cos\theta \end{bmatrix}

X1X_1轴旋转的矩阵为:

RX(γ)=[1000cosγsinγ0sinγcosγ]R_X(\gamma)=\begin{bmatrix} 1&0&0\\ 0&\cos\gamma&\sin\gamma\\ 0&-\sin\gamma&\cos\gamma \end{bmatrix}

四元数表示旋转

使用四元数表示旋转,首先需要将空间中的三维向量使用四元数表示出来,这里需要考虑四维超球的投影。空间中任意一个三维向量都可以表示为四维超球与a=0a=0交平面上的一点,即:

p=0+pxi+pyj+pzkp=0+p_x\vec{i}+p_y\vec{j}+p_z\vec{k}

设某三维向量p0p_0在三维空间中发生一次旋转qq,到达p1p_1的位置,旋转过程可以使用四元数q=q0+q1i+q2j+q3kq=q_0+q_1\vec{i}+q_2\vec{j}+q_3\vec{k}表示,则有:

p0=qp1qp_0=q\cdot p_1\cdot q^*

一定注意这里是结果在与四元数旋转在相乘!

计算该式,得到:

p1=(q0+q1i+q2j+q3k)(0+px1i+py1j+pz1k)(q0q1iq2jq3k)=[px1(q02+q12q22q32)+2py1(q1q2q0q3)+2pz1(q1q3+q0q2)]i+[2px1(q2q1+q0q3)+py1(q02q12+q22q32)+2pz1(q2q3q0q1)]j+[2px1(q3q1q0q2)+2py1(q3q2+q0q1)+pz1(q02q12q22+q32)]k    {px0=px1(q02+q12q22q32)+2py1(q1q2q0q3)+2pz1(q1q3+q0q2)py0=2px1(q2q1+q0q3)+py1(q02q12+q22q32)+2pz1(q2q3q0q1)pz0=2px1(q3q1q0q2)+2py1(q3q2+q0q1)+pz1(q02q12q22+q32)\begin{aligned} p_1&=(q_0+q_1\vec{i}+q_2\vec{j}+q_3\vec{k})\cdot(0+p_{x1}\vec{i}+p_{y1}\vec{j}+p_{z1}\vec{k})\cdot(q_0-q_1\vec{i}-q_2\vec{j}-q_3\vec{k})\\ &=[p_{x1}(q_0^2+q_1^2-q_2^2-q_3^2)+2p_{y1}(q_1q_2-q_0q_3)+2p_{z1}(q_1q_3+q_0q_2)]\vec{i}\\ &+[2p_{x1}(q_2q_1+q_0q_3)+p_{y1}(q_0^2-q_1^2+q_2^2-q_3^2)+2p_{z1}(q_2q_3-q_0q_1)]\vec{j}\\ &+[2p_{x1}(q_3q_1-q_0q_2)+2p_{y1}(q_3q_2+q_0q_1)+p_{z1}(q_0^2-q_1^2-q_2^2+q_3^2)]\vec{k}\\ \implies&\left\{ \begin{aligned} p_{x0}&=p_{x1}(q_0^2+q_1^2-q_2^2-q_3^2)+2p_{y1}(q_1q_2-q_0q_3)+2p_{z1}(q_1q_3+q_0q_2)\\ p_{y0}&=2p_{x1}(q_2q_1+q_0q_3)+p_{y1}(q_0^2-q_1^2+q_2^2-q_3^2)+2p_{z1}(q_2q_3-q_0q_1)\\ p_{z0}&=2p_{x1}(q_3q_1-q_0q_2)+2p_{y1}(q_3q_2+q_0q_1)+p_{z1}(q_0^2-q_1^2-q_2^2+q_3^2) \end{aligned} \right. \end{aligned}

请记住这个公式,它是后面一切转换的来源

特别地,对于单位四元数,(即逆等于其共轭的四元数),有:

p1=qp0q=q1p0q\begin{aligned} p_1&=q^*\cdot p_0\cdot q\\ &=q^{-1}\cdot p_0\cdot q \end{aligned}

千万不要记反了顺序!

欧拉角与四元数的关系

根据上节公式,现对三轴欧拉角逐个讨论其四元数等价式。首先对于绕Z0Z_0轴旋转的情况,其旋转矩阵将向量p0=px0i+py0j+pz0k\vec{p}_0=p_{x0}\vec{i}+p_{y0}\vec{j}+p_{z0}\vec{k}旋转至p1=px1i+py1j+pz1k\vec{p}_1=p_{x1}\vec{i}+p_{y1}\vec{j}+p_{z1}\vec{k},则其反过程是:

[px0py0pz0]=[cosψsinψ0sinψcosψ0001][px1py1pz1]    {px0=px1cosψpy1sinψpy0=px1sinψ+py1cosψpz0=pz1\begin{bmatrix} p_{x0}\\ p_{y0}\\ p_{z0} \end{bmatrix}=\begin{bmatrix} \cos\psi&-\sin\psi&0\\ \sin\psi&\cos\psi&0\\ 0&0&1 \end{bmatrix}\begin{bmatrix} p_{x1}\\ p_{y1}\\ p_{z1} \end{bmatrix}\implies\left\{ \begin{aligned} p_{x0}&=p_{x1}\cos\psi-p_{y1}\sin\psi\\ p_{y0}&=p_{x1}\sin\psi+p_{y1}\cos\psi\\ p_{z0}&=p_{z1} \end{aligned} \right.

将上式与四元数等价式进行比较,可以得到以下九项式子:

{q02+q12q22q32=cosψq1q2q0q3=12sinψq1q3+q0q2=0q2q1+q0q3=12sinψq02q12+q22q32=cosψq2q3q0q1=0q3q1q0q2=0q3q2+q0q1=0q02q12q22+q32=1\left\{ \begin{aligned} q_0^2+q_1^2-q_2^2-q_3^2&=\cos\psi\\ q_1q_2-q_0q_3&=-\frac{1}{2}\sin\psi\\ q_1q_3+q_0q_2&=0\\ q_2q_1+q_0q_3&=\frac{1}{2}\sin\psi\\ q_0^2-q_1^2+q_2^2-q_3^2&=\cos\psi\\ q_2q_3-q_0q_1&=0\\ q_3q_1-q_0q_2&=0\\ q_3q_2+q_0q_1&=0\\ q_0^2-q_1^2-q_2^2+q_3^2&=1\\ \end{aligned} \right.

本式实际上也给出了已知四元数,求解旋转矩阵的方法。直接求解这个方程,显然可知q1=q2=0q_1=q_2=0,然后进一步代回原式,得到:

{q0=cosψ2q1=0q2=0q3=sinψ2\left\{ \begin{aligned} q_0&=\cos\frac{\psi}{2}\\ q_1&=0\\ q_2&=0\\ q_3&=\sin\frac{\psi}{2} \end{aligned} \right.

因此,绕Z0Z_0轴旋转的四元数为:

qZ=cosψ2+sinψ2kq_Z=\cos\frac{\psi}{2}+\sin\frac{\psi}{2}\vec{k}

同理,可以得到绕YY'轴旋转的四元数为:

qY=cosθ2+sinθ2jq_Y=\cos\frac{\theta}{2}+\sin\frac{\theta}{2}\vec{j}

X1X_1轴旋转的四元数为:

qX=cosγ2+sinγ2iq_X=\cos\frac{\gamma}{2}+\sin\frac{\gamma}{2}\vec{i}

牢记定义的逆时针为正

旋转矩阵与四元数的关系

欧拉角与四元数有着显著和优美的联系,而实际中的旋转矩阵就没有如此好看的表达式,一般人估计硬记也记不住,求解其的方法仍然需要从四元数的等价式出发,只不过这里不再是欧拉轴角旋转矩阵,而是一般意义上的旋转矩阵。仍然有:

{px0=r11px1+r12py1+r13pz1py0=r21px1+r22py1+r23pz1pz0=r31px1+r32py1+r33pz1\left\{ \begin{aligned} p_{x0}&=r_{11}p_{x1}+r_{12}p_{y1}+r_{13}p_{z1}\\ p_{y0}&=r_{21}p_{x1}+r_{22}p_{y1}+r_{23}p_{z1}\\ p_{z0}&=r_{31}p_{x1}+r_{32}p_{y1}+r_{33}p_{z1} \end{aligned} \right.

对比四元数等价式,可以得到:

{q02+q12q22q32=r11q1q2q0q3=12r12q1q3+q0q2=12r13q2q1+q0q3=12r21q02q12+q22q32=r22q2q3q0q1=12r23q3q1q0q2=12r31q3q2+q0q1=12r32q02q12q22+q32=r33\left\{ \begin{aligned} q_0^2+q_1^2-q_2^2-q_3^2&=r_{11}\\ q_1q_2-q_0q_3&=\frac{1}{2}r_{12}\\ q_1q_3+q_0q_2&=\frac{1}{2}r_{13}\\ q_2q_1+q_0q_3&=\frac{1}{2}r_{21}\\ q_0^2-q_1^2+q_2^2-q_3^2&=r_{22}\\ q_2q_3-q_0q_1&=\frac{1}{2}r_{23}\\ q_3q_1-q_0q_2&=\frac{1}{2}r_{31}\\ q_3q_2+q_0q_1&=\frac{1}{2}r_{32}\\ q_0^2-q_1^2-q_2^2+q_3^2&=r_{33} \end{aligned} \right.

本方程更加复杂,不过也是能求解的,首先加入先验的单位四元数条件,即q02+q12+q22+q32=1q_0^2+q_1^2+q_2^2+q_3^2=1,然后将第一式、第四式和第九式相加,得到:

4q021=r11+r22+r334q_0^2-1=r_{11}+r_{22}+r_{33}

由此可以解出q0q_0,但是需要注意q_0应该取正号还是负号,实际上,只要保证q1,q2,q3q_1,q_2,q_3均与q0q_0耦合,就不用担心符号问题了,单位四元数本身符号不影响旋转结果。因此解得:

{q0=12r11+r22+r33+1q1=r32r234q0q2=r13r314q0q3=r21r124q0\left\{ \begin{aligned} q_0&=\frac{1}{2}\sqrt{r_{11}+r_{22}+r_{33}+1}\\ q_1&=\frac{r_{32}-r_{23}}{4q_0}\\ q_2&=\frac{r_{13}-r_{31}}{4q_0}\\ q_3&=\frac{r_{21}-r_{12}}{4q_0} \end{aligned} \right.

上式也给出了已知四元数求解旋转矩阵的方法。

四元数表示旋转

所谓表示旋转,就是将欧拉角的三轴旋转过程用四元数表达即可。仍然按照Z0YX1Z_0-Y'-X_1的顺序,旋转矩阵如上式可写为:

R=RZ(ψ)RY(θ)RX(γ)R=R_Z(\psi)R_Y(\theta)R_X(\gamma)

而复合旋转的四元数将与上式有相同的结构,如下所示:

q=qZqYqX=(cosψ2+sinψ2k)(cosθ2+sinθ2j)(cosγ2+sinγ2i)q=q_Zq_Yq_X=\left(\cos\frac{\psi}{2}+\sin\frac{\psi}{2}\vec{k}\right)\left(\cos\frac{\theta}{2}+\sin\frac{\theta}{2}\vec{j}\right)\left(\cos\frac{\gamma}{2}+\sin\frac{\gamma}{2}\vec{i}\right)

可将复合后的四元数直接视为一个四元数,即可完成旋转后向量的表示计算。

参考文献

  1. Understanding Quaternions
  2. Du Peng, Hu Haibao, Ding Dong, and Li Zhuoyue. Understanding Quaternions. 1st ed. New York, USA: Nova Science Publishers, Inc., 2020. https://lccn.loc.gov/2020037162.