1. 为什么我们要死磕四元数转旋转矩阵
做三维开发的人,迟早会撞上四元数这堵墙。你可能是做游戏引擎的,可能是搞机器人运动学的,也可能是在做飞行器姿态解算,甚至只是想在三维软件里写个插件让模型转得自然一点。不管哪条路,只要涉及姿态表示,四元数和旋转矩阵这两个东西就绕不开。四元数紧凑、插值平滑、不会出现万向节死锁,适合做运算和存储;旋转矩阵直观、易于和向量相乘、方便与其它线性代数工具衔接,适合做变换和渲染。所以实际工程里最常见的操作就是:内部用四元数维护姿态,需要作用于向量或传给渲染管线时,转成旋转矩阵。
这个转换本身并不复杂,公式网上一搜一大把。但我见过太多人直接复制一段代码,跑起来发现旋转方向反了,或者转出来的矩阵行列式是负的,又或者归一化没做导致缩放被悄悄引入。问题不在于公式难,而在于大多数人没有真正推导过一遍,不知道每个元素是怎么来的,也就无法在出错时定位。这篇文章我就从一个从业者的角度,把四元数转旋转矩阵这件事从头推一遍,把每一步的来龙去脉讲清楚,再补上实操中真正会踩的坑。目标很明确:看完之后你不光能默写出公式,还能在别人代码出错时一眼看出问题在哪。
适合的读者包括三维图形开发、机器人姿态控制、惯性导航算法、物理引擎相关方向的从业者和学生。只要你会基本的向量运算和矩阵乘法,剩下的推导我尽量用生活化的方式讲明白。
2. 先把基础概念捋直:四元数到底在表达什么
2.1 四元数的结构与单位约束
四元数可以写成 q = w + xi + yj + zk,通常记作 q = (w, x, y, z),其中 w 是实部,x、y、z 是虚部。在姿态表示的场景里,我们用的都是单位四元数,也就是满足 w² + x² + y² + z² = 1。这个约束极其重要,后面推导和实操都依赖它。
单位四元数和旋转的对应关系是这样的:它表示绕某个单位轴 n = (nx, ny, nz) 旋转角度 θ,那么
- w = cos(θ/2)
- x = nx · sin(θ/2)
- y = ny · sin(θ/2)
- z = nz · sin(θ/2)
注意这里是半角。很多人第一次接触会疑惑为什么是 θ/2 而不是 θ,这跟四元数的代数结构有关,简单说是因为四元数乘法作用在向量上时会出现两次旋转的复合效果,用半角才能让最终结果对应一次 θ 的旋转。这个半角特性也是四元数能平滑插值(比如球面线性插值)的根本原因。
2.2 四元数旋转向量的核心公式
四元数旋转一个三维向量 v 的标准做法是构造一个纯虚四元数 p = (0, v),然后计算
p' = q · p · q*
其中 q* 是 q 的共轭,对单位四元数来说共轭等于逆,即 q* = (w, -x, -y, -z)。算出来的 p' 依然是纯虚四元数,它的虚部就是旋转后的向量 v'。
这个公式是整个推导的起点。它看起来简洁,但直接拿去做工程有两个问题:一是每次旋转向量都要做两次四元数乘法,计算量偏大;二是很多系统(渲染管线、物理引擎)只接受矩阵形式的变换。所以我们希望把 q · p · q* 这个操作等价地写成一个 3×3 矩阵 R 乘以 v,也就是 v' = R · v。接下来的推导,本质上就是把四元数乘法展开,整理成矩阵形式。
2.3 为什么不用欧拉角直接转矩阵
有人会问,既然欧拉角转矩阵那么直观,为什么不干脆用欧拉角?原因在于欧拉角有三个致命问题:万向节死锁、插值不自然、依赖旋转顺序约定。四元数没有这些问题,所以工业界普遍用四元数做内部表示。而四元数转矩阵这一步,就是把“好用的运算表示”翻译成“好用的变换表示”,是管线里必不可少的一环。理解了这一点,你就明白为什么这个转换值得认真推导,而不是抄个公式了事。
3. 推导全过程:从四元数乘法到旋转矩阵
3.1 把四元数乘法展开成代数式
先定义两个四元数的乘法。设 a = (aw, ax, ay, az),b = (bw, bx, by, bz),则
a · b = (aw·bw - ax·bx - ay·by - az·bz, aw·bx + ax·bw + ay·bz - az·by, aw·by - ax·bz + ay·bw + az·bx, aw·bz + ax·by - ay·bx + az·bw)
这个乘法规则来自虚数单位的关系 i² = j² = k² = ijk = -1。记住这个展开式,后面全靠它。
现在令 q = (w, x, y, z),p = (0, vx, vy, vz),q* = (w, -x, -y, -z)。我们要算 t = q · p,再算 p' = t · q*。
先算 t = q · p。因为 p 的实部为 0,代入乘法公式:
- tw = -x·vx - y·vy - z·vz
- tx = w·vx + y·vz - z·vy
- ty = w·vy - x·vz + z·vx
- tz = w·vz + x·vy - y·vx
这一步就是把 p 的分量代进去,注意 p 的实部 bw = 0,所以很多项直接消失了。
3.2 完成第二次乘法并提取虚部
接着算 p' = t · q*。q* = (w, -x, -y, -z),代入乘法公式,我们只关心 p' 的虚部,因为实部理论上应该是 0(旋转后仍是纯虚四元数)。
p'x = tw·(-x) + tx·w + ty·(-z) - tz·(-y) = -tw·x + tx·w - ty·z + tz·y
p'y = tw·(-y) - tx·(-z) + ty·w + tz·(-x) = -tw·y + tx·z + ty·w - tz·x
p'z = tw·(-z) + tx·(-y) - ty·(-x) + tz·w = -tw·z - tx·y + ty·x + tz·w
现在把前面 tw、tx、ty、tz 的表达式代进来,展开整理。这一步是纯代数,比较繁琐,但没有任何技巧,就是耐心。我建议你拿张纸自己推一遍,推过一次之后印象会非常深。
以 p'x 为例,代入后:
p'x = -(-x·vx - y·vy - z·vz)·x + (w·vx + y·vz - z·vy)·w - (w·vy - x·vz + z·vx)·z + (w·vz + x·vy - y·vx)·y
展开各项:
= x²·vx + xy·vy + xz·vz
- w²·vx + wy·vz - wz·vy
- wz·vy + xz·vz - z²·vx
- wy·vz + xy·vy - y²·vx
按 vx、vy、vz 合并同类项:
- vx 的系数:x² + w² - z² - y²
- vy 的系数:xy - wz - wz + xy = 2xy - 2wz
- vz 的系数:xz + wy + wy + xz = 2xz + 2wy
所以 p'x = (w² + x² - y² - z²)·vx + (2xy - 2wz)·vy + (2xz + 2wy)·vz。
3.3 利用单位约束化简对角线元素
上面 vx 的系数是 w² + x² - y² - z²。利用单位四元数约束 w² + x² + y² + z² = 1,可以把它改写成 1 - 2(y² + z²),或者等价地 2(w² + x²) - 1。这两种写法在不同场景下各有便利,代码里常见的是 1 - 2(y² + z²) 这种形式,因为它只需要算两个平方。
同理对 p'y 和 p'z 做完全一样的展开整理,最终得到:
- p'y = (2xy + 2wz)·vx + (w² - x² + y² - z²)·vy + (2yz - 2wx)·vz
- p'z = (2xz - 2wy)·vx + (2yz + 2wx)·vy + (w² - x² - y² + z²)·vz
把这三个式子写成矩阵形式 v' = R · v,就得到经典的旋转矩阵:
R = [ 1-2(y²+z²) 2(xy-wz) 2(xz+wy) ] [ 2(xy+wz) 1-2(x²+z²) 2(yz-wx) ] [ 2(xz-wy) 2(yz+wx) 1-2(x²+y²) ]
这就是最终结果。注意矩阵里每一项的符号,尤其是第二、三项里 wz、wx、wy 的正负号,这是最容易抄错的地方。
3.4 用旋转轴角验证公式正确性
推导完不能就这么信了,得验证。最简单的验证方式是取一个绕 z 轴旋转 θ 的四元数:q = (cos(θ/2), 0, 0, sin(θ/2)),即 w = cos(θ/2),z = sin(θ/2),x = y = 0。
代入矩阵:
- R11 = 1 - 2z² = 1 - 2sin²(θ/2) = cos θ
- R12 = 2(xy - wz) = -2wz = -2cos(θ/2)sin(θ/2) = -sin θ
- R21 = 2(xy + wz) = 2wz = sin θ
- R22 = 1 - 2(x² + z²) = 1 - 2sin²(θ/2) = cos θ
- R33 = 1 - 2(x² + y²) = 1
得到的矩阵是
[ cosθ -sinθ 0 ] [ sinθ cosθ 0 ] [ 0 0 1 ]
这正是绕 z 轴旋转 θ 的标准旋转矩阵。验证通过。同样的方法你可以取绕 x 轴、y 轴的特例再验一遍,确保三个方向都对。我强烈建议每个第一次推导的人都做这个验证,它能帮你确认符号约定没有搞反。
4. 代码实现与工程落地细节
4.1 一份可直接用的转换实现
把上面的矩阵翻译成代码,核心就是一个函数。下面用 Python 写一版,逻辑清晰,方便对照公式:
def quaternion_to_matrix(w, x, y, z): # 建议先归一化,防止输入不是单位四元数 norm = (w*w + x*x + y*y + z*z) ** 0.5 w, x, y, z = w/norm, x/norm, y/norm, z/norm xx, yy, zz = x*x, y*y, z*z xy, xz, yz = x*y, x*z, y*z wx, wy, wz = w*x, w*y, w*z return [ [1 - 2*(yy + zz), 2*(xy - wz), 2*(xz + wy)], [2*(xy + wz), 1 - 2*(xx + zz), 2*(yz - wx)], [2*(xz - wy), 2*(yz + wx), 1 - 2*(xx + yy)] ]这段代码有几个细节值得说。第一,我提前算了 xx、yy 这些中间量,避免重复乘法,这在被频繁调用的渲染循环里能省下可观的浮点运算。第二,归一化放在函数内部做,虽然多了一次开方,但能防止上游传进来的四元数因为累积误差而偏离单位长度。第三,矩阵按行返回,如果你的图形库用的是列主序(比如某些 OpenGL 约定),记得转置。
4.2 归一化到底该不该做
这是个经常被争论的点。理论上如果四元数始终是单位长度,归一化就是多余的。但实际工程里,四元数经过多次乘法、插值、积分之后,长度会慢慢漂移。漂移带来的后果是矩阵不再正交,旋转里混入了缩放和剪切,模型会看起来“发胖”或者“扭曲”。
我的经验是:在转换函数内部做一次归一化,成本极低,收益极大。除非你在做极致的性能优化,且能保证上游每一步都严格归一化,否则不要省这一步。曾经有个项目因为省了归一化,跑了几个小时后场景里的物体开始肉眼可见地变形,排查了半天才定位到是四元数漂移。这个坑我替你踩过了。
4.3 性能优化:什么时候值得展开
如果你的转换每秒只调用几次,怎么写都无所谓。但在粒子系统、骨骼动画、大规模刚体模拟里,这个函数可能每秒被调用几十万次。这时候可以考虑:
- 把中间量提取出来复用,减少乘法次数
- 如果上游保证单位四元数,去掉归一化
- 用 SIMD 指令并行处理多个四元数
- 在 GPU 上做转换,把矩阵直接算在着色器里
不过要提醒一句,优化之前先测。我见过有人花大力气优化这个函数,结果 profile 一看它只占总耗时的 0.5%,纯属白费功夫。先定位瓶颈,再动手。
5. 常见问题与排查技巧实录
5.1 旋转方向反了怎么办
这是最高频的问题。转出来的结果和预期差一个反向,通常有三个原因:
第一,四元数的共轭用错了。旋转公式是 q · p · q*,如果你写成了 q* · p · q,得到的就是反向旋转。检查你的乘法顺序。
第二,矩阵转置了。行主序和列主序搞混,表现就是旋转方向相反。确认你的图形 API 用的是哪种约定。
第三,四元数本身的符号约定。q 和 -q 表示同一个旋转,但如果你在构造四元数时半角符号搞错,也会导致方向问题。
排查方法很简单:取一个绕单轴旋转 90 度的四元数,手动算一遍矩阵,和代码输出对比。差异一眼就能看出来。
5.2 矩阵行列式不为 1
正常的旋转矩阵行列式必须等于 1。如果你算出来是 -1 或者别的值,说明有问题。行列式为负通常是符号错误或者矩阵构造错误;行列式不为 1 但为正,多半是四元数没有归一化导致的缩放。
快速检查:把矩阵三列分别求模,正常应该都等于 1,且两两正交。如果模不等于 1,先查归一化。
5.3 常见问题速查表
| 现象 | 可能原因 | 排查方向 |
|---|---|---|
| 旋转方向相反 | 共轭顺序错、矩阵转置 | 检查 q·p·q* 顺序和存储约定 |
| 物体变形缩放 | 四元数未归一化 | 在转换前做归一化 |
| 特定角度出错 | 半角用错、符号抄错 | 用轴角特例验证 |
| 数值不稳定 | 累积误差漂移 | 定期重新归一化 |
| 性能不达标 | 重复计算、未复用中间量 | profile 后针对性优化 |
5.4 几个容易被忽略的实操心得
第一,写单元测试。取若干已知的轴角组合,构造四元数,转矩阵,再和直接用轴角公式算的矩阵对比。这个测试能挡住 90% 的低级错误。
第二,注意浮点精度。单精度浮点在极端角度下可能出现微小误差,如果对精度要求高,考虑用双精度做中间计算。
第三,文档里写清楚符号约定。团队协作时,四元数的分量顺序(wxyz 还是 xyzw)和矩阵的主序约定一定要统一,否则对接时必然出问题。我见过两个模块因为分量顺序不一致,联调了一整天才发现。
6. 从公式到直觉:怎么真正理解这个转换
推导做完、代码写完,其实还可以再往深走一步:建立对这个转换的直觉。旋转矩阵的三列,本质上就是旋转后的三个基向量。也就是说,R 的第一列是原 x 轴 (1,0,0) 旋转后的结果,第二列是原 y 轴旋转后的结果,第三列是原 z 轴旋转后的结果。你可以把 v = (1,0,0) 代入公式验证一下,得到的正好是矩阵第一列。
这个视角很有用。当你不确定矩阵某个元素该是什么符号时,就想象一个基向量转过去应该落在哪里,符号自然就出来了。比如绕 z 轴转 90 度,x 轴应该转到 y 轴方向,那么矩阵第一列应该是 (0,1,0),对照公式就能确认符号。
另一个直觉是:四元数的四个分量里,w 控制旋转角度,x、y、z 控制旋转轴。当 w = 1 时没有旋转,矩阵退化成单位阵;当 w = 0 时旋转 180 度,矩阵里只剩下轴相关的项。理解了这一点,你就能大致预判不同四元数会产生什么样的矩阵,调试时心里有底。
我个人在实际操作中的体会是,这个转换看起来是个纯数学问题,但真正决定成败的是工程细节:归一化做没做、符号约定统不统一、测试覆盖够不够。公式本身半小时就能推完,但把这些细节都处理妥当,往往需要在实际项目里摔打几轮。所以如果你现在正在做相关开发,别急着抄公式,先自己推一遍,再写几个测试用例,后面能省下大量排查时间。