尧图网站设计 尧图网站设计YAOTU DESIGN
ARTICLE DETAIL

资讯详情

深耕网站设计与一线实操的经验洞察。

四元数Hamilton与JPL约定:姿态解算避坑指南

四元数Hamilton与JPL约定:姿态解算避坑指南 在姿态解算这个圈子里四元数几乎是人人都用、但很少有人第一次就用对的东西。尤其是当你把一个模块从一个平台搬到另一个平台或者把一篇论文里的公式照抄进代码又或者从第三方设备手册里抠出一组数直接喂给自己的滤波器时十有八九会遇到同一个场景模长检查过了正交性检查过了数值也看不出异常可姿态就是反的或者转出来的角度差了 180 度。这个问题八成不在你的算法上而在你手上那组数到底是Hamilton notation还是JPL notation。这两个名字听起来像学术考据实际上它们直接决定了四元数的排列顺序、虚部符号、乘法展开、旋转矩阵形式甚至决定了角速度微分方程的矩阵长什么样。做机器人、惯性导航、视觉惯性融合、飞控、游戏引擎、仿真工具链的人只要碰姿态就绕不开这道坎。下面的内容我按为什么会分两套—差异到底在哪—手推一遍公式—姿态更新怎么做—移植时怎么三分钟定位错误—工具链怎么选这条线走公式和代码都给全目标是你看完能直接对着自己的工程改。1. 两套约定同时存在坑到底埋在哪1.1 用复数理解四元数一个不用背公式的入口先把记忆负担卸掉。四元数并不神秘它就是复数的二维推广。复数用一个实部加一个虚部表示平面上的旋转四元数用一个实部加三个虚部表示三维空间的旋转。你可以把三维旋转想象成绕某根轴转某个角度而四元数就是把这根轴和这个角度打包成四个数。理解到这一层你就能接受一件事既然打包方式是人定的那不同的人按不同顺序打包就必然产生不同约定。具体来说记 q w xi yj zk其中 w 叫标量部分x、y、z 叫矢量部分。三个虚单位的关系是 i² j² k² ijk -1。这组关系对应的是 Hamilton 的原始定义也是绝大多数机器人、SLAM、惯导代码默认采用的版本。它的几何意义很直观矢量部分的方向就是旋转轴方向矢量部分的长度是 sin(θ/2)标量部分是 cos(θ/2)其中 θ 是旋转角。这里有个初学者最容易忽略的点——为什么是半角。因为四元数对矢量做旋转时q 在公式里出现了两次一次左乘、一次右乘它的共轭两次叠加起来才凑成一个完整的 θ。这也顺带解释了四元数的双覆盖性质q 和 -q 表示同一个旋转因为两次取负在乘法里抵消了。这个性质在后面做插值和误差计算时会反复咬人先记住。1.2 三处差异顺序、乘法表、旋转语义Hamilton 和 JPL 的差异不是单点的而是三层叠加的这也是为什么很多人改了一处以为搞定了、结果还是错。第一层是存储顺序。Hamilton 习惯写成 (w, x, y, z)标量在前JPL 习惯写成 (x, y, z, w)标量在后。这一层最容易被发现也最容易被人当成就是换个顺序而已。但只改顺序往往不够。第二层是乘法表的方向。Hamilton 用 ij k交换两个虚单位会变号。JPL 采用的等价基可以理解成把三个虚单位整体取负会让乘法表的符号反过来。这一层决定了 q ⊗ p 和 p ⊗ q 哪个对应先转 p 再转 q。第三层是旋转语义。Hamilton 下的主动旋转写作 v q ⊗ v ⊗ q*q* 是共轭而 JPL 下同一个物理效果对应的公式方向相反等价的说法是它的旋转矩阵是 Hamilton 形式的转置。这一层决定了你的姿态是体系到惯性还是惯性到体系。三层叠加之后同一组四个数在两套约定下代表的物理姿态是不同的。更要命的是如果你只改了第一层剩下两层会让你得到一个看起来能跑、测出来是反的系统。1.3 为什么看起来全对、姿态全反我见过最典型的翻车现场是这样的从设备手册里抄来一组四元数标量在后于是你把它重排成标量在前喂给自己的 Hamilton 滤波器。滤波器不发散模长守恒输出的旋转矩阵行列式是 1正交性误差在 1e-12 量级。可你拿它去控制机械臂末端跑到了镜像位置。原因就在第二、三层。因为 JPL 的旋转矩阵与 Hamilton 的旋转矩阵互为转置转置对应的物理含义是旋转方向取反逆旋转。转置的旋转矩阵同样是正交的、行列式为 1所以所有合法性检查都能通过唯独方向错了。这种错误不会报异常只会让整个系统的工作方向反过来。还有一种更隐蔽的情况你的四元数是从某个中间件拿到的文档只说了标量在最后但没说旋转语义是哪个。这种情况下最省事的判断方式就是做一次单点测试——构造一个绕 z 轴转 90 度的姿态把自己那组数喂进旋转公式看矢量 (1, 0, 0) 被转到哪里。转到 (0, 1, 0) 就是你要的那种转到 (0, -1, 0) 就说明差了一次转置。这个测试具体怎么写我在第 5 章给完整脚本。2. 把 Hamilton 与 JPL 的差异逐条摊开2.1 存储顺序标量在前还是在后存储顺序的差异本身不涉及数学但它决定了你所有的索引访问是否写对。Hamilton 写法 q (w, x, y, z)代码里 q[0] 是标量。JPL 写法 q (x, y, z, w)代码里 q[3] 是标量。这个差异看起来微不足道可一旦你的代码里有几十处手写索引而中间又经过了几次重构和移植索引写错几乎是必然事件。我的建议是在任何跨模块的接口上禁止裸数组传递四元数。要么用一个带命名成员的类或者结构体把w、x、y、z做成具名字段要么在函数签名和注释里写清楚约定并且在入口处立刻做一次规范化转成你自己工程内部的统一表示。裸double q[4]在接口上传递是这个领域最常见的定时炸弹。另外要特别注意很多库的对外顺序和内部顺序是不一致的。典型例子是 Eigen构造的时候写Quaterniond(w, x, y, z)但coeffs()返回的是 (x, y, z, w)。也就是说它的构造接口是标量在前内部存储是标量在后而乘法语义又是 Hamilton 的。这种东西只看书是看不出来的必须动手试。2.2 基元关系ijk 等于 -1 还是 1Hamilton 的基元满足 i² j² k² ijk -1。如果你把三个虚单位整体取负即定义 I -i、J -j、K -k那么 I² J² K² -1 仍然成立但 IJK (-i)(-j)(-k) -ijk 1。这就是 JPL 那套基的来历本质上不是新的数学对象而是同一套结构换了一组符号。这个符号翻转的直接后果是乘法展开式里的交叉项全部反号。Hamilton 下两个纯四元数 p (0, u)、q (0, v) 相乘结果是 (-u·v, u × v)。换成反向基之后同样两个量相乘的矢量部分变成 -u × v。也就是说两套约定下乘法这个操作的物理含义是不同的。工程上你不需要记住两套乘法表只需要记住一句话JPL 的乘法与 Hamilton 的乘法相当于把两个因子的顺序颠倒过来。所以如果你的代码里用的是 Hamilton 实现而某个接口给的是 JPL 参数你不能只是重排数组还要把每次连乘的顺序反过来。这句话在实际调试时非常有用因为一旦你发现姿态是反的第一时间就可以去检查乘法顺序。2.3 旋转方向主动与被动的一线之隔主动旋转和被动旋转的区分是理解这两个约定最有效的一把钥匙。主动旋转是说物体转了坐标系不动被动旋转是说坐标系转了物体不动。从数学上这两者互为逆运算。Hamilton 的经典公式 v q ⊗ v ⊗ q* 描述的是主动旋转给定惯性系下的矢量 v乘完之后得到旋转后的新矢量。JPL 那套参数配合它的旋转矩阵使用时得到的是这个旋转的转置也就是逆旋转等价于被动旋转。在实际系统里这个差别对应的是一个非常具体的问题你的四元数描述的是体系到惯性系的旋转还是惯性系到体系的旋转。这两个恰好互为共轭。惯导里通常用 q 表示体系到惯性系也就是把体坐标系下的量转到导航系飞控里有的代码用反过来的一支。一旦搞混你的加速度补偿、重力投影全部反号表现出来就是静止时倾角算出来是负的或者飞机往左偏姿态显示往右。所以我一直强调一个习惯在每个四元数变量的命名里写清方向。比如q_b2n表示体系到导航系q_n2b表示导航系到体系。名字里带上方向比写十行注释都管用因为它会在你写错的时候当场露馅。2.4 差异对照总表把上面三层差异整理成表方便你贴在工位上对照对比项Hamilton notationJPL notation存储顺序(w, x, y, z)标量在前(x, y, z, w)标量在后基元关系i²j²k²ijk-1等价于取负基IJK1旋转公式v q ⊗ v ⊗ q*v q* ⊗ v ⊗ q旋转矩阵R(q)R(q)ᵀ乘法顺序q⊗p 表示先后 p 后 q与 Hamilton 相反输入角速度体系角速度右乘体系角速度左乘矩阵形式不同常见领域机器人、SLAM、惯导、游戏航天、轨道、部分仿真工具链注意这张表描述的是真正意义上的 JPL 约定。现实中大量的库自称标量在后但旋转语义仍然是 Hamilton 的也就是只改了顺序、没改语义。这种情况必须靠单点测试来区分不能靠名字。3. 从轴角一路推到旋转矩阵手算一遍3.1 轴角构造四元数与半角的来历给定单位旋转轴 u (ux, uy, uz) 和旋转角 θHamilton 约定下构造四元数的方法是标量部分 w cos(θ/2)矢量部分 (x, y, z) (ux, uy, uz) · sin(θ/2)最后做一次归一化理论上已经是单位四元数但浮点运算后建议再归一化一次举个能验证的例子绕 z 轴转 90 度u (0, 0, 1)θ 90°那么 w cos45° ≈ 0.7071z sin45° ≈ 0.7071x y 0。写成 Hamilton 顺序是 (0.7071, 0, 0, 0.7071)写成 JPL 顺序是 (0, 0, 0.7071, 0.7071)。到这里两种写法数值只差顺序。但如果按第 2 章的表真正的 JPL 表示要把矢量部分取负也就是 (0, 0, -0.7071, 0.7071)输出的旋转矩阵才是绕 z 转 90 度而不是 -90 度。半角为什么必须存在可以从罗德里格斯公式反推。把 q (s, u) 代入 v q ⊗ v ⊗ q*展开后矢量部分是v (s² - u·u) v 2 u (u·v) 2 s (u × v)代入 s cos(θ/2)、|u| sin(θ/2)用二倍角公式化简第一项系数变成 cos θ第二项系数变成 (1 - cos θ)第三项系数变成 sin θ。于是得到v cos θ · v (1 - cos θ) · u(u·v) sin θ · (u × v)这恰好就是经典的罗德里格斯旋转公式。也就是说四元数旋转和轴角旋转是同一件事半角是被两次乘法自然逼出来的不是谁拍脑袋定的。3.2 四元数乘法的完整展开与代码Hamilton 约定下的乘法展开写成代码比写成公式好记。设 p (w1, x1, y1, z1)q (w2, x2, y2, z2)import numpy as np def qmul_h(p, q): Hamilton 约定 (w, x, y, z) 的四元数乘法 p ⊗ q w1, x1, y1, z1 p w2, x2, y2, z2 q 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 qconj(q): w, x, y, z q return np.array([w, -x, -y, -z]) def qnorm(q): n np.linalg.norm(q) return q / n if n 1e-12 else np.array([1.0, 0.0, 0.0, 0.0])这段代码看着简单但它包含了两个容易出错的细节。第一个是虚部交叉项的符号x 分量里的y1*z2 - z1*y2对应的是 y × z 的方向关系写成-y1*z2 z1*y2就变成了逆序。第二个是当输入四元数不是单位长度时乘法结果的长度是两个输入长度之积如果你在链式乘法里忘了归一化误差会以乘积形式累积。验证乘法是否写对最快的办法是做一个双次旋转测试构造绕 x 轴转 90 度的 q1 和绕 x 轴再转 90 度的 q2q1 和 q2 数值相同。用你的乘法算 q1 ⊗ q2期望得到绕 x 轴转 180 度的四元数即 (0, 1, 0, 0)。这个测试能一次性抓出符号和顺序两类错误。3.3 DCM 的两种形式与互推把上一节的矢量展开写成矩阵形式Hamilton 约定下旋转矩阵DCM是def quat2dcm_h(q): Hamilton 约定 (w, x, y, z) 转旋转矩阵v_new R v_old w, x, y, z qnorm(q) return np.array([ [1-2*(y*yz*z), 2*(x*y-w*z), 2*(x*zw*y)], [2*(x*yw*z), 1-2*(x*xz*z), 2*(y*z-w*x)], [2*(x*z-w*y), 2*(y*zw*x), 1-2*(x*xy*y)], ])JPL 约定的旋转矩阵在同样的四个数排列下就是这个矩阵的转置。这个结论非常好用因为它给了你一个不依赖具体公式记忆的判据如果你手上有一个 JPL 参数的四元数想做 Hamilton 的主动旋转最省事的做法是把矩阵转置一次或者等价地把四元数取共轭再重排。在实际工程里我会在内部统一压成一种表示只要进到我的模块里一律先转成 Hamilton 的标量在前形式然后所有运算都用这一套实现。JPL 只作为进出口的序列化格式出现在边界上做一次转换函数就三个def ham_to_scalar_last(q_h): 仅换顺序旋转语义保持 Hamilton。Eigen/SciPy/ROS 风格 w, x, y, z q_h return np.array([x, y, z, w]) def ham_to_jpl_true(q_h): 真正 JPL 语义换顺序 矢量取负保证 R_J(jpl) R_H(ham) w, x, y, z q_h return np.array([-x, -y, -z, w]) def jpl_true_to_ham(q_j): x, y, z, w q_j return np.array([w, -x, -y, -z])这两个转换函数的区别就是我在 1.3 节说的那个致命歧义。你必须先确认对方到底是哪一种选错了姿态就直接反过来。3.4 从 DCM 反解四元数的分支处理把旋转矩阵反解成四元数比正推更容易踩坑因为它涉及除以接近零的数。教科书上的做法是用迹trace判断分支当 trace(R) 0 时用 w 做主元w sqrt(1 trace) / 2其余分量由矩阵元素差得到当 trace(R) ≤ 0 时取对角线最大的那一项做主元分别处理 x、y、z 三种情况。为什么必须分支理由很实在直接用 w sqrt(1trace)/2当旋转角接近 180 度时 trace 接近 -11trace 接近零开方后精度崩掉解出来的三个分量会被放大成噪声。用最大对角元做主元能保证分母始终远离零。这个技巧在姿态初始化、对准、以及长时间运行后的重同步场景里特别关键。还有一点反解出来的四元数与真值可能差一个整体符号因为双覆盖。如果你在做误差统计直接拿解出来的数和参考数相减会得到误差接近 2的荒谬结果。正确做法是先判断两者的点积符号如果为负就把其中一个整体取负再计算误差。4. 姿态更新微分方程、积分与归一化4.1 角速度耦合的两种写法姿态更新的核心是四元数微分方程。Hamilton 约定下如果 q 表示体系到惯性系的旋转而陀螺输出的角速度 ω 是体系下的量方程是q̇ 0.5 · q ⊗ ω其中 ω 作为纯四元数 (0, ωx, ωy, ωz) 参与乘法。展开成标量形式是ẇ -0.5 · (x·ωx y·ωy z·ωz)ẋ 0.5 · (w·ωx y·ωz - z·ωy)ẏ 0.5 · (w·ωy z·ωx - x·ωz)ż 0.5 · (w·ωz x·ωy - y·ωx)而 JPL 约定下标准的矩阵形式是 q̇ 0.5 · Ω(ω) · q其中 Ω 是一个 4×4 反对称矩阵其元素排布与上式的展开一一对应只是把标量放到了第四位。你可以自己推一遍验证把上面的 ẋ、ẏ、ż、ẇ 按 (x, y, z, w) 顺序重排成矩阵得到的正是那个经典的 Ω 矩阵。这也再次说明两套约定在物理层面是等价的差别只在书写和符号约定但写错一处就会让整个积分跑偏。这里有个实践中的判断准则你用什么坐标系下的角速度决定了是左乘还是右乘。体系角速度右乘惯性系角速度左乘。反过来用短时间看不出来长时间会积累成明显的姿态漂移。很多开源代码在注释里会写gyro in body frame这就是在提示你用右乘。4.2 一阶欧拉、闭式解与误差对比最简单的更新方式是一阶欧拉def integrate_euler(q, omega, dt): 一阶欧拉omega 为体系角速度Hamilton 约定 wq np.array([0.0, omega[0], omega[1], omega[2]]) dq 0.5 * qmul_h(q, wq) * dt return qnorm(q dq)它的截断误差是 O(dt²)在 dt 为 1 毫秒、角速度几百度每秒的场景下勉强够用但在高机动场景下会明显偏。想提升精度有两种常见做法。一种是用闭式解。假设在 dt 内角速度恒定模长 ω |ω|那么精确解是右乘一个增量四元数q_{k1} q_k ⊗ [cos(ω·dt/2), sin(ω·dt/2) · ω/ω]这个方法的好处是除了角速度恒定这个假设之外没有其他近似而且天然保持单位长度数值表现远好于一阶欧拉。代价是每次更新要算一次三角函数。实际测下来在 100 Hz 到 1 kHz 的更新率下这个开销完全可以接受。另一种是用中值或四阶龙格库塔。中值法的做法是用区间中点的角速度估计一次导数精度是二阶计算量只比欧拉多一点。四阶龙格库塔精度更高但每次更新要做四次四元数乘法和一次归一化在资源紧张的嵌入式平台上要权衡。下面这张表是我自己做过的一个参考对比条件是从静止开始给一个恒定角速度积分 60 秒用不同方法算末端姿态角误差积分方法阶数计算量相对60 秒累积误差量级一阶欧拉11.010⁻¹ 度中值法21.610⁻³ 度闭式解定角速度精确1.410⁻⁶ 度受浮点限制四阶龙格库塔44.510⁻⁶ 度注意表里的误差量级和采样率强相关。采样率越低一阶欧拉的劣势越明显。如果你只能跑欧拉法先把采样率提上去效果往往比换积分方法更立竿见影。4.3 归一化与零偏补偿的实操归一化这件事看起来是每次更新完除一下模长这么简单实际上有两个讲究。第一是不要每次都强行归一化。如果你用的是浮点运算归一化本身会引入舍入误差高频归一化反而会让低频分量被反复截断。比较折中的做法是每若干步归一化一次或者直接使用闭式解让更新本身保长。我一般会在每步更新后用平方和判断是否偏离 1 超过阈值比如 1e-6只有当超限的时候才做归一化。第二是零偏一定要在四元数层面之外处理。陀螺的零偏应该在角速度积分之前扣除而不是在姿态上打补丁。如果你把零偏修正塞进四元数更新里会让误差模型的物理意义混乱后期做标定和滤波的时候非常难拆。还有一个容易被忽略的细节是双覆盖在误差计算里的影响。假设你有一个估计四元数 q_est 和一个参考四元数 q_ref想算误差角正确做法是def angle_error_deg(q_est, q_ref): qe qnorm(q_est) qr qnorm(q_ref) if np.dot(qe, qr) 0: qr -qr # 消除双覆盖取最短弧 dq qmul_h(qconj(qe), qr) w min(1.0, abs(dq[0])) return 2.0 * np.degrees(np.arccos(w))少了if np.dot(...) 0这一句你的误差曲线会时不时冒出接近 360 度的尖峰很容易被误判成滤波器发散。这个坑我当年排查了整整一个下午最后发现只是符号问题。5. 移植排查实录怎么三分钟定位符号错5.1 三层自检法每次从外部引入一组四元数或者把代码搬到新平台我都会按固定顺序跑三层检查。第一层是结构检查看模长、看共轭关系、看两个合法旋转相乘后是否仍然是合法单位四元数。这一层能抓出索引写错、归一化漏做这类低级问题。第二层是矩阵一致性检查用同一组数分别算旋转矩阵然后检查 RᵀR 与单位阵的偏差、行列式是否接近 1。这层能抓出公式抄错的问题但它抓不出转置错误因为转置后的矩阵一样正交、一样行列式为 1。第三层是语义检查也就是单点已知角测试。构造一个已知的旋转看结果落在哪里。这一层是唯一能区分反了和对了的检查。前两层过了、第三层没过答案基本就锁定在约定问题上。顺序不能颠倒因为如果结构就有问题第三层的结论也靠不住。5.2 单点已知角测试脚本这个脚本我建议你不管现在有没有问题都先跑一遍把结果记下来当作基线。以后每次换库、升级依赖、换设备重新跑一次对比。import numpy as np def quat2dcm_h(q): w, x, y, z q / np.linalg.norm(q) return np.array([ [1-2*(y*yz*z), 2*(x*y-w*z), 2*(x*zw*y)], [2*(x*yw*z), 1-2*(x*xz*z), 2*(y*z-w*x)], [2*(x*z-w*y), 2*(y*zw*x), 1-2*(x*xy*y)], ]) def single_point_test(q, label): R quat2dcm_h(q) v_new R np.array([1.0, 0.0, 0.0]) print(f{label:28s} (1,0,0) - f({v_new[0]:.4f}, {v_new[1]:.4f}, {v_new[2]:.4f})) s np.sqrt(0.5) single_point_test([s, 0, 0, s], Hamilton 绕 z 90 度) single_point_test([s, 0, 0, -s], Hamilton 绕 z -90 度) single_point_test([0, 0, s, s], 标量在后、语义 Hamilton) single_point_test([0, 0, -s, s], 真 JPL 参数重排后喂 Hamilton)期望结果很明确绕 z 轴 90 度时(1, 0, 0) 应该被转到 (0, 1, 0) 附近绕 z 轴 -90 度时应该转到 (0, -1, 0) 附近。第三行是只换了顺序、语义没变的情况结果同样应该是 (0, 1, 0)。第四行是把真 JPL 参数直接重排后喂进 Hamilton 公式的典型错误结果会变成 (0, -1, 0)。把这个对比跑出来你以后一眼就能认出问题出在哪。5.3 常见问题速查表下面这个表是我这些年攒下来的问题清单基本都是检查全过但结果错这一类。建议收藏出问题的时候从第一行往下试。现象最可能的原因处理办法姿态方向恰好相反约定混淆导致的一次转置用 5.2 的脚本确认语义层绕某个轴转动时方向反了其他轴正常单个分量符号写错检查乘法展开式里的交叉项符号静止时倾角为负四元数方向定义反了体系到惯性 vs 惯性到体系取共轭或换用反向量测位置长时间运行后缓慢漂移角速度用了错误的坐标系确认是体系角速度右乘还是左乘误差曲线上出现接近 360 度的尖峰双覆盖没处理误差计算前统一点积符号两个模块数据对接后姿态跳变一个模块标量在前、一个在后在接口层做一次显式转换航向角在过 ±180 度时跳动四元数本身连续是转欧拉角时的分支问题在欧拉角转换处做角度解缠滤波收敛但姿态有小幅抖动归一化过于频繁引入截断噪声改为阈值触发式归一化6. 工具链的约定分布与接口设计6.1 常见库的倾向对照这部分信息变化很快不同版本之间还可能有差异所以下面只说倾向具体一定以你自己手上的版本加 5.2 的单点测试为准。表里不做绝对承诺。库或平台构造/输入顺序语义倾向Eigen构造按 (w, x, y, z)内部 coeffs 为 (x, y, z, w)Hamilton 主动旋转SciPy Rotationfrom_quat 接受 (x, y, z, w)Hamilton 主动旋转ROS geometry_msgs字段顺序 x, y, z, wHamiltonMATLAB Robotics Toolbox约定为 (s, v)标量在前HamiltonMATLAB Aerospace Toolbox以官方文档为准构造习惯标量在前Hamilton部分航天仿真工具链(x, y, z, w) 或 (s, v) 都有文档JPL矩阵为转置形式常见物理引擎(w, x, y, z) 或 (x, y, z, w) 并存多数 Hamilton注意左手系这张表最有价值的一行是 Eigen。因为它的构造和存储顺序不一致很多人第一次写都会写反而且写反之后程序不崩只是姿态不对。我的习惯是永远用Quaterniond q(w, x, y, z)构造需要取出数组的时候显式调用.coeffs()并立刻转换成我自己内部的标量在前表示中间不留下任何歧义空间。6.2 把语义写进接口名和文档真正能根治这类问题的办法不是记住差异而是让差异在你的代码里无法被忽略。第一招是命名带方向。q_b2n、q_n2b、q_enu_from_body这类名字比q、quat、quat_out强一百倍。它在编译期不会报错但在你读代码的时候会立刻提示方向是否搞反。第二招是约定带后缀。q_h表示 Hamilton 标量在前q_jpl表示 JPL 语义q_sl表示只换顺序、语义仍为 Hamilton。所有跨模块函数都带上这个后缀。这样即使半年后接手的人不懂这些名词看名字也能猜到需要转换。第三招是转换函数只有一个出口。不要在五个地方各写一遍重排代码那必然有一处写错。把所有转换收敛到一个模块里每个函数都配一个单点测试用起来才放心。第四招是在单元测试里固化语义。写一个测试用固定输入固定输出把绕 z 转 90 度这件事钉死。以后有人改动了乘法或者积分测试立刻报警。这类测试写起来三行代码挡住的却是最难查的那类问题。6.3 几条踩坑后总结的经验关于这件事我自己的体会是约定混乱带来的成本比任何算法本身的复杂度都高。四元数的公式早就写死了真正花时间的是搞清手上这组数到底代表什么。我在项目里做过一个决定后面再也没有被姿态方向的问题困扰过就是在系统边界上做一次强制规范化。所有外部输入——设备、文件、网络消息、其他模块——进来第一步就转成内部统一的 Hamilton 标量在前表示转换函数只有两个标量在后版、真 JPL 版每个都配单点测试。内部所有计算只用这一套包括乘法、转矩阵、积分、误差计算、插值。出口再按需要转回去。这样做的代价是每个接口多一次拷贝收益是整个系统只有一个地方可能出错而且那个地方有测试兜着。另一个经验是关于插值的。做姿态平滑或者传感器时间同步时线性插值四元数是不行的必须用球面插值。而球面插值的第一步就是判断两个四元数的点积符号取最短弧。这个判断在两套约定的代码里都要做不要以为换了一种约定就不需要。最后提醒一个很容易被忽略的场景从文件或者日志里回放数据做离线分析。这时候数据往往是几个月前录的录制时的代码版本和现在不同约定可能已经改过。回放之前先跑一次单点测试比你调半天滤波器参数管用得多。我自己就吃过这个亏当时以为是滤波器的过程噪声调得不合适折腾了两天最后发现是日志里的四元数在某个版本之后从标量在前改成了标量在后。如果后续要做的事更大一些比如把一个完整的姿态模块迁移到另一个平台我建议的推进顺序是先把两套表示之间的转换和单点测试建起来再把乘法、转矩阵、积分逐个对齐最后才去接传感器数据。中间任何一步对不上就停下来把差异定位清楚再往下走别指望后面能自动对上。
返回列表