简介:这份PDF文档围绕四元数与三维旋转之间的关系展开,面向计算机图形学、机器人学、计算机视觉方向的学习者与开发者,帮助读者从几何与线性代数角度理解四元数如何表示三维空间中的旋转。内容从复数与2D旋转的类比切入,逐步推导四元数的定义、加减法与乘法运算规则,并给出四元数到旋转矩阵的对应公式,最后延伸至实际应用场景。资源包内共1个PDF文件,大小约478KB,篇幅紧凑,适合作为姿态表示与旋转计算的入门参考或复习材料。目前已有182人学习浏览。读者可借此理清四元数乘法与旋转复合的对应关系,掌握旋转矩阵的构造方式,并理解其在动画、机器人姿态解算中的使用逻辑,为后续工程实践打下数学基础。
1. 四元数 PDF 到底在讲什么:从姿态解算翻车说起
如果你写过飞控、做过机械臂正逆解,或者只是想在 Unity 里让一个物体平滑转向,你大概率绕不开四元数。而当你搜到一份名为quaternion.pdf的资料时,真正想知道的往往不是“四元数是什么”,而是:它能不能帮我解决手头这个姿态解算的玄学问题?这份 PDF 里那套公式,我能不能直接抄进代码里跑通?
四元数最核心的价值,是用四个数表达三维旋转,避免欧拉角的万向锁,也比旋转矩阵更省内存、更适合插值。一份典型的quaternion.pdf会覆盖定义、乘法、共轭、归一化、转旋转矩阵、球面线性插值(Slerp)这几块。它面向的是做 IMU 姿态融合、机器人运动学、图形学动画的工程师。我见过太多人把四元数当成黑匣子,乘法顺序写反、共轭当逆用、插值不归一化,结果姿态输出抖得像筛糠。这篇笔记就顺着这份 PDF 的知识脉络,把每个公式落到可复现的代码和参数上,让你少走我当年走过的弯路。
2. 四元数基础运算:把 PDF 里的公式翻译成能跑的代码
一份quaternion.pdf的前半部分,基本都在定义和推导基本运算。很多人看公式觉得懂了,一写代码就错,问题出在没把“约定”钉死。四元数有实部在前的(w, x, y, z),也有虚部在前的(x, y, z, w),PDF 里用哪种,你的代码就必须跟哪种,混用是血泪级的翻车源头。
2.1 定义、单位四元数与三个必须固定的约定
四元数写成 $q = w + xi + yj + zk$,其中 $i^2=j^2=k^2=ijk=-1$。表示旋转时,单位四元数满足 $w^2+x^2+y^2+z^2=1$,对应绕单位轴 $\hat{u}$ 转 $\theta$ 角:
$$q = \left[\cos\frac{\theta}{2},\ \hat{u}\sin\frac{\theta}{2}\right]$$
这里第一个坑就是半角。为什么是 $\theta/2$ 不是 $\theta$?因为四元数乘法对旋转是“双覆盖”,$q$ 和 $-q$ 表示同一旋转。PDF 里通常会推导这一点,但实操时你只需要记住:角度进去先除以 2。
三个必须固定的约定:第一,分量顺序,我一般统一用(w, x, y, z);第二,乘法是左乘还是右乘,对应旋转是作用在向量上还是坐标系上;第三,旋转方向是主动旋转(转向量)还是被动旋转(转坐标系)。这三个约定不统一,后面所有公式都会错。
import numpy as np def quat_normalize(q): """归一化四元数,输入输出均为 (w, x, y, z)""" q = np.asarray(q, dtype=float) n = np.linalg.norm(q) if n < 1e-12: raise ValueError("零四元数无法归一化") return q / n def quat_from_axis_angle(axis, angle_rad): """绕单位轴 axis 旋转 angle_rad,返回 (w, x, y, z)""" axis = np.asarray(axis, dtype=float) axis = axis / np.linalg.norm(axis) # 轴必须先单位化 half = angle_rad * 0.5 # 半角,别漏 w = np.cos(half) xyz = axis * np.sin(half) return np.array([w, xyz[0], xyz[1], xyz[2]])逻辑说明:quat_normalize里设了1e-12的阈值,因为浮点累乘后模长会漂移,直接除零会得到 NaN,姿态就彻底崩了。quat_from_axis_angle强制先单位化轴,如果你传进来的轴本身带模长,结果就不是单位四元数,后续转矩阵会引入缩放。参数上,angle_rad用弧度,别传角度,我见过有人传 90 进去,结果转了 90 弧度,物体直接飞出屏幕。
2.2 四元数乘法与共轭:顺序错了姿态就反了
四元数乘法是这套体系里最容易写错的地方。给定 $q_1=(w_1,x_1,y_1,z_1)$、$q_2=(w_2,x_2,y_2,z_2)$,Hamilton 积为:
$$q_1 q_2 = \begin{bmatrix} w_1w_2 - x_1x_2 - y_1y_2 - z_1z_2 \ w_1x_2 + x_1w_2 + y_1z_2 - z_1y_2 \ w_1y_2 - x_1z_2 + y_1w_2 + z_1x_2 \ w_1z_2 + x_1y_2 - y_1x_2 + z_1w_2 \end{bmatrix}$$
注意这里没有交换律,$q_1q_2 \neq q_2q_1$。物理含义是:先做 $q_2$ 的旋转再做 $q_1$ 的旋转,等价于 $q_1q_2$(在主动旋转约定下)。如果你把顺序写反,姿态就会朝反方向转,而且在小角度下不明显,大角度下直接翻车。
共轭 $q^* = (w, -x, -y, -z)$,对单位四元数来说,共轭等于逆,$q^{-1} = q^*$。这一点非常重要:只有单位四元数才能这么偷懒,非单位四元数的逆要除以模长平方。
def quat_multiply(q1, q2): """Hamilton 积,输入输出 (w, x, y, z)""" w1, x1, y1, z1 = q1 w2, x2, y2, z2 = q2 return np.array([ w1*w2 - x1*x2 - y1*y2 - z1*z2, w1*x2 + x1*w2 + y1*z2 - z1*y2, w1*y2 - x1*z2 + y1*w2 + z1*x2, w1*z2 + x1*y2 - y1*x2 + z1*w2, ]) def quat_conjugate(q): """共轭,单位四元数下等价于逆""" w, x, y, z = q return np.array([w, -x, -y, -z]) def quat_rotate_vector(q, v): """用单位四元数 q 旋转三维向量 v""" q = quat_normalize(q) v_quat = np.array([0.0, v[0], v[1], v[2]]) q_conj = quat_conjugate(q) # v' = q * v * q^-1 return quat_multiply(quat_multiply(q, v_quat), q_conj)[1:]逻辑说明:quat_rotate_vector里先把向量补成实部为 0 的纯四元数,再做 $qvq^{-1}$。这里调用了两次quat_multiply,顺序不能换。参数上,q进来先归一化,防止累积误差导致旋转带缩放。如果你发现旋转后的向量模长变了,八成是q没归一化,或者共轭用成了非单位四元数的逆。
2.3 转旋转矩阵与欧拉角:对接现有管线的桥梁
很多现有系统只吃旋转矩阵或欧拉角,所以quaternion.pdf一定会讲转换。单位四元数转旋转矩阵:
$$R = \begin{bmatrix} 1-2(y^2+z^2) & 2(xy-wz) & 2(xz+wy) \ 2(xy+wz) & 1-2(x^2+z^2) & 2(yz-wx) \ 2(xz-wy) & 2(yz+wx) & 1-2(x^2+y^2) \end{bmatrix}$$
转欧拉角时要注意旋转顺序,常见的有 ZYX(偏航-俯仰-滚转)和 XYZ。不同顺序结果完全不同,PDF 里如果没写清楚,你得从上下文推断。
def quat_to_matrix(q): """单位四元数转 3x3 旋转矩阵""" w, x, y, z = quat_normalize(q) return np.array([ [1-2*(y*y+z*z), 2*(x*y-w*z), 2*(x*z+w*y)], [2*(x*y+w*z), 1-2*(x*x+z*z), 2*(y*z-w*x)], [2*(x*z-w*y), 2*(y*z+w*x), 1-2*(x*x+y*y)], ]) def quat_to_euler_zyx(q): """转 ZYX 欧拉角,返回 (roll, pitch, yaw),单位弧度""" w, x, y, z = quat_normalize(q) # roll (绕 X) sinr = 2*(w*x + y*z) cosr = 1 - 2*(x*x + y*y) roll = np.arctan2(sinr, cosr) # pitch (绕 Y) sinp = 2*(w*y - z*x) sinp = np.clip(sinp, -1.0, 1.0) # 防浮点越界 pitch = np.arcsin(sinp) # yaw (绕 Z) siny = 2*(w*z + x*y) cosy = 1 - 2*(y*y + z*z) yaw = np.arctan2(siny, cosy) return roll, pitch, yaw逻辑说明:quat_to_matrix直接套公式,注意它假设q是单位四元数,非单位会引入缩放。quat_to_euler_zyx里对sinp做了clip,因为浮点误差可能让它略大于 1,arcsin会返回 NaN,这是姿态显示突然变 NaN 的常见原因。参数上,返回的是弧度,要显示成角度自己乘180/pi。欧拉角在俯仰接近 ±90° 时会万向锁,这也是为什么内部解算要用四元数,只在最后显示时才转欧拉角。
3. 姿态解算实战:从 IMU 数据到稳定四元数
quaternion.pdf讲完基础运算,通常会进入应用部分,而最典型的应用就是 IMU 姿态解算。你手上有陀螺仪和加速度计,陀螺仪积分得到角度但会漂移,加速度计测重力方向但噪声大。四元数在这里的角色,是把两者融合成一个稳定的姿态表示。
3.1 陀螺仪积分:四元数微分方程怎么离散
陀螺仪输出角速度 $\omega = (\omega_x, \omega_y, \omega_z)$,四元数随时间演化为:
$$\dot{q} = \frac{1}{2} q \otimes \begin{bmatrix}0 \ \omega\end{bmatrix}$$
离散化后,常用一阶近似:
$$q_{k+1} = q_k + \frac{\Delta t}{2} q_k \otimes \begin{bmatrix}0 \ \omega\end{bmatrix}$$
然后归一化。这里 $\Delta t$ 是采样周期,单位秒。如果你用 100Hz 采样,$\Delta t = 0.01$。
def quat_integrate_gyro(q, omega, dt): """陀螺仪积分更新四元数 q: 当前姿态 (w,x,y,z) omega: 角速度 (rad/s) 三元组 dt: 采样周期 秒 """ q = quat_normalize(q) omega_quat = np.array([0.0, omega[0], omega[1], omega[2]]) dq = 0.5 * quat_multiply(q, omega_quat) * dt q_new = q + dq return quat_normalize(q_new) # 必须归一化逻辑说明:quat_integrate_gyro里先构造纯四元数omega_quat,再算 $q \otimes \omega$,乘 0.5 和dt得到增量。关键在最后一步归一化,因为一阶近似会引入模长误差,不归一化的话积分几分钟后四元数模长就飘了,姿态跟着缩放。参数上,omega必须是弧度每秒,很多 IMU 原始输出是度每秒,记得先转。dt要尽量准,用实际时间戳差值比用固定周期更稳。
3.2 加速度计修正:用重力方向拉回漂移
纯陀螺仪积分会漂,需要用加速度计观测的重力方向来修正。思路是:把当前四元数估计的重力方向 $\hat{g}_{est} = R(q) \cdot [0,0,1]^T$ 和加速度计实测方向 $\hat{a}$ 做叉积,得到误差向量,再用它去修正角速度。
def quat_correct_accel(q, accel, gain=0.02): """用加速度计修正姿态漂移 q: 当前四元数 accel: 加速度计读数 (g),三元组 gain: 修正增益,典型 0.01~0.05 """ q = quat_normalize(q) a = np.asarray(accel, dtype=float) n = np.linalg.norm(a) if n < 1e-6: return q # 自由落体或无效数据,跳过 a = a / n # 单位化实测重力方向 # 估计重力方向:把机体 Z 轴转到世界系 g_est = quat_rotate_vector(q, np.array([0.0, 0.0, 1.0])) # 叉积得到误差 error = np.cross(a, g_est) # 用误差修正角速度(这里简化成直接修正四元数) if np.linalg.norm(error) < 1e-8: return q correction = gain * error # 把修正量转成小角度四元数增量 angle = np.linalg.norm(correction) axis = correction / angle dq = quat_from_axis_angle(axis, angle) return quat_normalize(quat_multiply(dq, q))逻辑说明:quat_correct_accel先判断加速度模长,太小说明在自由落体或数据无效,直接跳过,否则会引入错误修正。gain是修正强度,太大姿态会跟着加速度噪声抖,太小漂移压不住,我一般从 0.02 起调。叉积np.cross(a, g_est)的方向决定了修正方向,如果你发现越修越偏,把叉积顺序反过来试试,这就是约定问题。参数上,accel单位是 g 还是 m/s² 不影响,因为做了归一化。
3.3 融合流程与参数整定:一个能跑的最小闭环
把积分和修正串起来,就是一个最小姿态解算闭环。每来一帧 IMU 数据:先用陀螺仪积分更新四元数,再用加速度计修正,最后归一化输出。
class AttitudeEstimator: def __init__(self, gain=0.02): self.q = np.array([1.0, 0.0, 0.0, 0.0]) # 初始水平 self.gain = gain def update(self, gyro, accel, dt): """gyro: rad/s, accel: 任意单位, dt: 秒""" # 1. 陀螺仪积分 self.q = quat_integrate_gyro(self.q, gyro, dt) # 2. 加速度计修正 self.q = quat_correct_accel(self.q, accel, self.gain) # 3. 输出 return self.q.copy() # 使用示例 est = AttitudeEstimator(gain=0.02) # 假设 100Hz,陀螺仪静止但有零偏,加速度计指向重力 for _ in range(1000): q = est.update(gyro=[0.01, 0.0, 0.0], accel=[0.0, 0.0, 1.0], dt=0.01) roll, pitch, yaw = quat_to_euler_zyx(q)逻辑说明:AttitudeEstimator把状态封在类里,初始四元数[1,0,0,0]表示无旋转。update里先积分后修正,顺序不能反,否则修正量会被积分带偏。参数上,gain是唯一需要整定的,静止时看 roll/pitch 是否收敛到 0,如果缓慢漂就加大,如果抖动就减小。陀螺仪零偏最好先标定,静止时取几百帧平均,从gyro里减掉,否则gain再大也压不住持续漂移。
4. 避坑与排查:四元数代码里那些让你怀疑人生的错误
四元数的问题有个特点:错了不一定报错,但结果就是不对,排查起来很折磨。下面这几条是我和身边人真实踩过的,按“现象 → 原因 → 解决”整理。
现象:姿态缓慢旋转,静止时也漂。原因:陀螺仪零偏没标定,或者dt用了名义值而非实际时间戳。解决:静止采集 500~1000 帧求陀螺仪均值,从每帧减去;dt用相邻时间戳差值,别用固定 0.01。
现象:姿态突然翻转 180° 或跳变。原因:四元数没归一化,模长漂移后转矩阵出现缩放;或者q和-q表示同一旋转但插值时走了远路。解决:每次更新后强制归一化;插值前判断点积,为负则取反其中一个。
现象:加速度计修正越修越偏。原因:叉积顺序反了,或者gain太大把噪声当信号。解决:确认np.cross(a, g_est)方向,反了就对调;gain从 0.01 开始,逐步加到能压住漂移又不抖。
现象:转欧拉角时 pitch 输出 NaN。原因:arcsin输入因浮点误差略超 ±1。解决:np.clip(sinp, -1.0, 1.0),这是后悔药级别的操作。
现象:Slerp 插值结果不平滑,中间突然加速。原因:没处理四元数双覆盖,两个四元数点积为负时走了长弧。解决:插值前若dot(q1,q2) < 0,令q2 = -q2,再插值。
注意:四元数乘法和共轭的约定一旦定下,全项目统一。我见过有人在姿态解算里用
(w,x,y,z),在图形显示里用(x,y,z,w),对接时没转换,物体姿态直接镜像,查了一整天才发现。
5. 进阶技巧:Slerp 插值与四元数双覆盖的工程处理
当你把基础解算跑通,quaternion.pdf后半部分通常会讲插值,而 Slerp 是四元数最漂亮的特性之一。它能在两个姿态之间沿最短弧匀速旋转,这是欧拉角做不到的。但工程上有个细节:四元数双覆盖意味着 $q$ 和 $-q$ 是同一姿态,如果直接对两个“数值上”的四元数插值,可能走了长弧,表现为旋转突然反向加速。
def quat_slerp(q1, q2, t): """球面线性插值,t in [0,1]""" q1 = quat_normalize(q1) q2 = quat_normalize(q2) dot = np.dot(q1, q2) # 双覆盖处理:点积为负说明走了长弧,取反 if dot < 0.0: q2 = -q2 dot = -dot # 点积接近 1 时退化为线性插值,避免除零 if dot > 0.9995: result = q1 + t * (q2 - q1) return quat_normalize(result) theta_0 = np.arccos(dot) # 夹角 theta = theta_0 * t sin_theta = np.sin(theta) sin_theta_0 = np.sin(theta_0) s1 = np.cos(theta) - dot * sin_theta / sin_theta_0 s2 = sin_theta / sin_theta_0 return quat_normalize(s1 * q1 + s2 * q2)逻辑说明:quat_slerp先归一化输入,算点积。dot < 0时把q2取反,这是处理双覆盖的关键,不做这一步插值会走长弧。dot > 0.9995时两个四元数几乎重合,sin_theta_0接近 0 会除零,所以退化成线性插值再归一化。参数上,t从 0 到 1,返回对应姿态。这个函数在动画、云台平滑转向里非常实用。
验证 Slerp 是否正确,有个简单方法:固定q1和q2,让t从 0 到 1 走 100 步,把每步转成欧拉角,看角度变化是否单调平滑。如果中间有跳变,八成是双覆盖没处理。另一个技巧是,如果你只需要插值一小段,可以用nlerp(线性插值后归一化)近似,计算量小很多,在dot > 0.9时视觉上几乎看不出差别,这是性能敏感场景的常用取舍。
我自己现在的习惯是:任何四元数进函数先归一化,任何插值先查点积符号,任何转欧拉角先 clip。这三条帮我省了无数个调试的夜晚。四元数不难,难的是约定和边界,把这两样钉死,它就是你手里最稳的姿态工具。希望帮到你。
本文还有配套的精品资源,点击获取