SLAM学习
首发于SLAM学习

第三讲:三维空间的刚体运动

题图来自monicore | Pixabay

三维空间的刚体运动

本讲主要主要关注一个问题:一个刚体在三维空间中的运动是如何描述的

为什么要研究这个问题?因为当我们描述机器人的位姿时,就是在描述一个刚体在三维空间中的运动。而且我们还要对这个估计出来的位姿进行进一步的优化。如何机器人的位姿的表述就决定了机器人的运动方程和优化方程的表述。

三维空间中,刚体的运动可以用两个概念来表示:旋转平移。平移比较简单一些,一般用一个表示位移的向量来表示。而旋转则有多种表示方法,例如旋转矩阵、旋转向量等等,不同的表示方法各有优劣。

向量

在描述旋转矩阵前我们先明确向量这个概念。向量是空间中的一个具体实物且不和任何实数相关联。为了描述向量,应该先确定一个具体的坐标系,明确该坐标系的线性基 [\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3] 后才能够确定一个向量  \mathbf{a} 在该坐标系下的坐标:

\mathbf{a}=[\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3] \begin{bmatrix} a_1 \\ a_2 \\ a_3 \end{bmatrix}

因此,向量坐标的具体取值,和向量本身选取的坐标系相关。

明确的向量的坐标后,就可以谈论向量间的运算,主要有两种:内积外积(又叫叉积):

\mathbf{a} \cdot \mathbf{b} =\mathbf{a}^T\mathbf{b}=|\mathbf{a}||\mathbf{b}|\cdot \cos\langle\mathbf{a}, \mathbf{b}\rangle=\sum\limits^3_{i=1}a_ib_i

\mathbf{a}\times\mathbf{b}= \begin{Vmatrix} \mathbf{i} & \mathbf{j} & \mathbf{k} \\ a_1 & a_2 &a_3 \\ b_1 & b_2 & b_3 \\ \end{Vmatrix} = \begin{bmatrix} a_2 b_3 - a_3 b_2 \\ a_3 b_1 - a_1 b_3 \\ a_1 b_2 - a_2 b_1 \\ \end{bmatrix} = \begin{bmatrix} 0 & -a_3 & a_2 \\ a_3 & 0 & -a_1 \\ -a_2 & a_1 & 0 \end{bmatrix} \mathbf{b} =\mathbf{a} ^\wedge \mathbf{b}

其中,符号 ^\wedge 表示取该向量的反对称矩阵。

外积的计算结果仍然是一个向量,方向垂直于 \mathbf{a}\mathbf{b} , 大小为 |\mathbf{a}||\mathbf{b}|\cdot \sin\langle\mathbf{a}, \mathbf{b}\rangle . 值得注意的是,外积还可以表示向量的旋转。因为任何旋转都可以描述为:绕某个旋转轴顺时针/ 逆时针旋转多少角度。后面会在详细描述。

欧式变换

假设存在两个坐标系:一个世界坐标系,定义为某个墙角和它的三条边,并是一个惯性系;一个机器人坐标系,是一个随机器人移动的坐标系。假设机器人观察到了某个向量 \mathbf{p} , 它在这两个坐标系中分别有一套坐标。前面说了,向量是一个客观存在的实体,那么必然有一个关系能够将这两套坐标联系起来。

这个关系就是欧式变换。因为机器人的运动是一个刚体运动,所以同一个向量在不同坐标系下的模长和方向都不会发生变化。这样一个欧式变换就是由一个旋转和一个平移两部分组成。

旋转的描述

旋转矩阵

我们先考虑旋转。对于向量 \mathbf{p} ,它在两个坐标系下的表示应该是相等的:

[\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3] \begin{bmatrix} a_1 \\ a_2 \\ a_3 \end{bmatrix} = [\mathbf{e}_1', \mathbf{e}_2', \mathbf{e}_3'] \begin{bmatrix} a_1' \\ a_2' \\ a_3' \end{bmatrix}

等式右边同时左乘 [\mathbf{e}_1^T, \mathbf{e}_2^T, \mathbf{e}_3^T]^T 可得

\begin{bmatrix} a_1 \\ a_2 \\ a_3 \end{bmatrix} = \begin{bmatrix} e_1^T e_1' & e_1^T e_2' & e_1^T e_3' \\ e_2^T e_1' & e_2^T e_2' & e_2^T e_3' \\ e_3^T e_1' & e_3^T e_2' & e_3^T e_3' \\ \end{bmatrix} \begin{bmatrix} a_1' \\ a_2' \\ a_3; \end{bmatrix} = \mathbf{R} \mathbf{a'}

我们把中间的矩阵拿出来,定义成一个矩阵 \mathbf{R} ,它描述了旋转本身,又称为旋转矩阵,是由两组基之间的内积组成的。而且它是一个行列式为1正交矩阵。它的逆(亦即转置)描述了一个相反的的旋转。

上述定义是一个充分必要条件,所以,我们可以将旋转矩阵的集合定义如下,称为特殊正交群

SO(n)=\{\mathbf{R}\in R^{n\times n} | \mathbf{R}\mathbf{R}^T=\mathbf{I},\ det(\mathbf{R})=1\}

旋转向量

旋转矩阵可以很好的描述相机的旋转了,但也有两个问题:

  • 一次旋转只有3个自由度,但旋转矩阵却有9个量。显然这冗余了;
  • 旋转矩阵自身带有两个约束。在进行估计和优化时,我们显然不想要这些附加的约束。

所以我们想要一个更紧凑的表达,最好就是3个量。而旋转向量就可以。

前面提到,任意旋转都可以用一个旋转轴和一个旋转角来描述。那么,我们可以用一个向量,其方向为旋转轴 \mathbf{n} 的方向,大小则为旋转角 \theta ,这种向量就是旋转向量,表示为 \theta\mathbf{n} . 前面说到,外积可以用来表示旋转就是因为外积可以用来表示旋转向量:考虑两个不平行的向量 \mathbf{a}\mathbf{b} ,在右手法则下,用右手的4个指头从 \mathbf{a} 转向 \mathbf{b} ,大拇指的朝向就是旋转向量的方向,即 \mathbf{a} \times \mathbf{b} 的方向,而它的大小则由 \mathbf{a}\mathbf{b} 的夹角决定。

旋转向量和旋转矩阵的转换可以用罗德里格斯公式表示:

\mathbf{R}=\cos \theta \mathbf{I} + (1-\cos\theta)\mathbf{n}\mathbf{n}^T+\sin\theta \ \mathbf{n}^\wedge

反之则有

tr(\mathbf{R})=\cos\theta \ tr(\mathbf{I})+(1-\cos\theta)\ tr(\mathbf{n}\mathbf{n}^T)+\sin\theta\ tr(\mathbf{n}^\wedge)=1+2\cos\theta

所以

\theta=arccos(\frac{tr(\mathbf{R})-1}{2})

而旋转轴上的向量在旋转后不发生改变,所以:

\mathbf{R}\mathbf{n}=\mathbf{n}

说明 \mathbf{n} 是矩阵 \mathbf{R} 的特征值 1 所对应的特征向量。求解后并归一化则可得到 \mathbf{n} .

欧拉角

除了旋转向量,我们也可以用欧拉角来紧凑地描述旋转。一个旋转可以分解成3次分别绕X,Y,Z 轴的旋转来表示,在航空摄影测量中,一般用“翻转 - 航偏 - 俯仰”(roll - yaw - pitch),也即XZY来表示。先绕X 轴旋转roll 角度,再绕Y 轴旋转yaw 角度,最后按照Z 轴旋转pitch 角度。这三个旋转矩阵相乘就得到了总的旋转矩阵。

此时,就可以用 [r, p, y]^T 这三个量来描述任意旋转,这种表示方式会比其他方式更直观、更易理解。

四元数

既然我们已经有了旋转向量和欧拉角,为什么还有个四元数(Quaternion)?因为欧拉角和旋转向量具有奇异性(万向锁问题)。不存在不带有奇异性的三维向量描述方式。因此我们需要用到四元数,它既是紧凑的,也没有奇异性

定义:一个四元数包含一个实部和三个虚部。

\mathbf{q}=q_0+q_1\mit{i}+q_2\mit{j}+q_3\mit{k} = [s,\mathbf{v}]

其中,后面的等式将四元数表达成一个标量和一个向量, \mit{i,j,k} 表示四元数的三个虚部,满足:

\begin{align} \begin{cases} i^2=j^2=k^2=-1 \\ ij=k,\ ji=-k \\ jk=i,\ kj=-i \\ ki=j,\ ji=-j \end{cases} \end{align}

若一个四元数虚部全为0,则它是一个实四元数;如其实部为零,则称它为虚四元素。而且,一个虚四元数对应一个空间点

我们能用单位四元数来表示三维空间中的任意一个旋转。我们先考虑下复数。在负数中,乘以i 表示在复平面内旋转90度。但在四元数中,情形却有微妙的不同:乘以i 表示旋转180度,这样才能保证 ij=-k 的性质。而 i^2=-1 ,说明绕i 轴旋转360度后得到一个相反的东西,而要旋转720度(两周)才能得到它原先的样子。(?)

假设某个旋转的旋转向量为 \theta\mathbf{n} , 则

\mathbf{q}= [ \cos\frac{\theta}{2},\ n_x\sin\frac{\theta}{2},\ n_y\sin\frac{\theta}{2},\ n_z\sin\frac{\theta}{2} ]^T

反之则有

\theta=2\cdot \arccos q_0 \\ [n_x,\ n_y,\ n_z]^T=[q_1,\ q_2,\ q_3]^T/ \sin\frac{\theta}{2}

上式给人一种“转了一半”的感觉。将上式中的 \theta 加上 2\pi 后得到一个相同的旋转,但是对应的四元数却变成了 -\mathbf{q} . 所以,在四元数中,任意的旋转都可以由两个互为相反数的四元数表示

而四元数和旋转矩阵的关系为:

\mathbf{R}= \begin{bmatrix} 1 - 2q_2^2 - 2q_3^2 & 2q_1q_2 - 2q_0q_3 & 2q_1q_3+2q_0q_2 \\ 2q_1q_2 + 2q_0q_3 & 1 - 2q_1^2 - 2q_3^2 & 2q_2q_3 - 2q_0q_1 \\ 2q_1q_3 - 2q_0q_2 & 2q_2q_3 + 2q_0q_1 & 1 - 2q_1^2 - 2q_2^2 \end{bmatrix}

设矩阵 \mathbf{R}=\{m_{ij}\},\ i,j \in [1, 2, 3] , 则由上式可以推得:

q_0 = \frac{\sqrt{tr(\mathbf{R}) + 1}}{2},\ q_1=\frac{m_{23}-m_{32}}{4q_0}, \ q_2=\frac{m_{31}-m_{13}}{4q_0},\ q3=\frac{m_{12}-m_{21}}{4q_0}

NOTE: 由于 \mathbf{q}-\mathbf{q} 表示同一个旋转,所以一个旋转矩阵对应的四元数表示并不惟一且存在其他转换公式。在实际中,如果 q_0 接近于0,会造成其他3个数的解不稳定,应采用其他公式。

四元数的运算

这里只列出关于四元数的简单运算。

现有两个四元数 \mathbf{q}_a,\ \mathbf{q}_b , 向量表示分别为 [s_a, \mathbf{v}_a], \ [s_b, \mathbf{v}_b]

  • 加减法

\mathbf{q}_a \pm \mathbf{q}_b=[s_a \pm s_b,\ \mathbf{v}_a \pm \mathbf{v}_b]

  • 乘法

乘法将两个四元数的每一项两两相乘然后相加。注意其满足之前所述的四元数的性质。

\mathbf{q}_a \pm \mathbf{q}_b=[s_a \pm s_b,\ \mathbf{v}_a \pm \mathbf{v}_b]

该乘法定义下,两个实四元数的乘积仍然是实四元数。由于最后一项外积的存在,四元数的乘法不满足交换律,除非 \mathbf{v}_a\mathbf{v}_bR^ 3 中共线(此时外积为零)。

  • 共轭

共轭即将虚部取为相反数:

\mathbf{q}_a^*=[s_a,\ -\mathbf{v}_a]

四元数共轭与其自身相乘可得一个实四元数,实部为模长的平方。

\mathbf{q}^*\mathbf{q}=\mathbf{q}\mathbf{q}^*=[s_a^2+\mathbf{v}^T\mathbf{v},\ 0]

  • 模长

\begin{Vmatrix}\mathbf{q}_a\end{Vmatrix}=\sqrt{s_a^2+x_a^2+y_a^2+z_a^2}

可知,两个四元数乘积的模即为二者模的乘积。所以,单位四元数相乘后仍然得到单位四元数。

\begin{Vmatrix}\mathbf{q}_a\mathbf{q}_b\end{Vmatrix}=\begin{Vmatrix}\mathbf{q}_a\end{Vmatrix}\begin{Vmatrix}\mathbf{q}_b\end{Vmatrix}

\mathbf{q}^{-1}=\frac{\mathbf{q}^*}{\begin{Vmatrix}\mathbf{q}\end{Vmatrix}}

可见,对于单位四元数,其逆和共轭相同。按照定义,四元数和自己逆的乘积为实四元数 \mathbf{1} . 最后,四元数乘积的逆有和矩阵相乘相似的性质:

(\mathbf{q}_a\mathbf{q}_b)^{-1}=\mathbf{q}_b^{-1}\mathbf{q}_a^{-1}

  • 数乘和点乘

四元数的数乘与点乘与向量相同:

k\mathbf{q}=[ks,\ k\mathbf{v}]

\mathbf{q}_a\mathbf{q}_b=s_as_b+x_ax_b \mit{i}+y_ay_b\mit{j}+z_az_b\mit{k}

用四元数表示旋转

要用四元数来计算旋转,首先,用一个虚四元数来表示三维空间点:

\mathbf{p}=[0,\ x, \ y, \ z]^T=[0,\mathbf{v}]

则旋转后的点 p' 则为

\mathbf{p}'=\mathbf{q}\mathbf{p}\mathbf{q}^{-1}

变换矩阵

解决了旋转问题,就可以考虑平移了。平移可以简单地用一个平移向量 \mathbf{t} 来表示。则变换后的坐标为

\mathbf{a}'=\mathbf{R}\mathbf{a}+\mathbf{t}

但是,上面这个式子并不是线性的。所以我们引入了齐次坐标。在一个三维向量末尾加上一个1,就得到一个四维向量,成为齐次坐标。上式就可以写成

\begin{bmatrix} \mathbf{a}' \\ 1 \end{bmatrix}= \begin{bmatrix} \mathbf{R} & \mathbf{t} \\ \mathbf{0}^T & 1 \end{bmatrix} \begin{bmatrix} \mathbf{a}'\\ 1 \end{bmatrix}= \mathbf{T} \begin{bmatrix} \mathbf{a} \\ 1 \end{bmatrix}

类似旋转矩阵 \mathbf{R} ,变换矩阵 \mathbf{T} 可以组成一个集合,叫特殊欧式群

SE(3)= \begin{Bmatrix} T= \begin{bmatrix} \mathbf{R} & \mathbf{t} \\ \mathbf{0}^T & 1 \end{bmatrix} \in R^{4 \times 4} |\ \mathbf{R} \in SO(3), \ t \in R^3 \end{Bmatrix}


最后,连发了两片之后才发现公式没有附上来。只好又在修改了一番。修改完后又发现知乎的编辑器的公式只能偶尔居中,然后还没有编号。后面的文章常常会写到 式(10) 之类的。只能在想办法了。

编辑于 2018-03-28

文章被以下专栏收录