☰
机器人位置与姿态描述:旋转矩阵、欧拉角与四元数详解
2026/9/30 1:27:15 网站建设 项目流程

调机械臂的时候,最让人抓狂的往往不是电机不转,而是位置到了、姿态却歪了十几度。我第一次被"位置与姿态描述"狠狠教育,是在一次抓取实验里:末端执行器停在了目标上方偏左大约三厘米的位置,姿态还整体绕自身轴多转了小半圈,夹爪对着空气开合。复盘下来,问题不在控制器参数,而在我对"姿态"的理解——我把它当成了三个可以像位置一样做加减的数。位置与姿态描述看上去只是机器人学导论里的第一章内容,实际上它是后面运动学、轨迹规划、手眼标定、力控这些东西的公共地基,地基歪一厘米,上层就歪一米。

这篇内容适合三类人看:刚上机器人学导论、被旋转矩阵和欧拉角绕晕的学生;手头要把机械臂、云台、无人机、移动平台的状态量对接起来的工程实现者;以及想把书本符号落成能跑代码的实践者。我会按"为什么难、怎么表示、怎么算、怎么排错"的顺序,把位置与姿态描述的几套主流表示法讲清楚,中间给出可以直接抄下来跑的 Python 片段,也会说说我踩过的那些坑。

1. 位置好写,姿态难缠:一次抓取失败拆出来的差别

1.1 位置是标准向量,三个数可以放心加减

位置这件事之所以简单,是因为它天然活在三维欧氏空间里。一个点 P 在参考系 A 下的位置,就是三个实数组成的有序三元组,写成列向量形式 p^A = [x, y, z]^T。这个上标 A 不是什么装饰,它表示"这组数字是在 A 坐标系里量出来的"。同一个点,在 A 系里是 [1, 0, 0],把参考系整体转 90 度之后,数字可能变成 [0, 1, 0],但点是同一个点。很多初学者的混乱就是从忽略这个上标开始的——手里拿着三个数,却说不清它们是相对谁的。

位置向量满足我们熟悉的全部代数性质:p1 + p2 有意义(位移叠加),k·p 有意义(缩放),p1 - p2 表示从 p2 指向 p1 的位移向量。这意味着位置可以用线性代数那一整套工具处理,求平均、做最小二乘、插值,都不会出问题。这也是为什么位置的代码通常写起来很顺手:三行数组,加减乘除,一目了然。

我在项目里有一个硬性习惯:任何一个位置变量命名里必须带参考系信息,比如p_base_tcp、p_cam_target。这不是洁癖。做手眼标定时,相机系下的目标点和基座系下的目标点数值上可能只差一个变换矩阵,但一旦命名混了,你就会拿着相机系的 z 去和基座系的 z 做比较,然后花半天时间怀疑标定结果。位置的坑基本都在"忘了它是相对谁"这一件事上。

1.2 姿态没有减法,因为旋转不满足交换律

姿态的麻烦在于,它不是向量空间里的元素,而是三维旋转群 SO(3) 里的元素。SO(3) 是一个群,有乘法(旋转的复合)和单位元(不转),但它不是线性空间——两个旋转相加没有意义,你也不能把一个旋转乘以 2 表示"转两倍那么多"。想转两倍,得用轴角形式把角度翻倍,而不是把旋转矩阵乘二。

更反直觉的是旋转的复合不可交换。你可以做个小实验:拿一本书放在桌上,先绕世界坐标系的 Z 轴转 90 度,再绕世界坐标系的 X 轴转 90 度;记下书的最终朝向。然后把书复位,把顺序反过来,先绕 X 转 90 度再绕 Z 转 90 度。两次结果完全不同,误差可能到 90 度量级。这就是为什么姿态描述里"先转哪个后转哪个"是必须写进文档的信息,而不能靠默认。

这里有个特别容易被混淆的点:角速度是可以相加的,有限转动不行。原因在于角速度是瞬时量,它活在姿态流形的切空间里,切空间是线性空间,所以能加起来;而有限转动是流形上的"点",两点之间没有加法。我见过不少人在多传感器融合里把 IMU 的三个角速度积分出来的角度直接相加,然后奇怪为什么姿态会漂。实际上那三个角速度积分出来的是一个旋转,得按旋转的复合规则乘起来,而不是向量加。

1.3 开工之前先把符号约定钉死

机器人学的符号约定是个重灾区,不同教材的同一个符号含义可能正好相反。我建议在动手写第一行代码之前,先在自己的项目里写一份一页纸的约定文档,把下面这些事定死:坐标系一律用右手系;位置和旋转都用列向量相乘的约定,也就是 p^A = R_AB · p^B;下标写成"目标系_参考系"还是"参考系_目标系"要统一;上标表示被描述量所在(或所表达于)的坐标系。

以 R_AB 为例,我采用的是"把 B 系中的坐标转换到 A 系"这个语义,也就是它同时是 B 系的三个基轴在 A 系中的排列。这个约定下 R_BA = R_AB^T,链式关系是 R_AC = R_AB · R_BC,非常好记。但如果你看的教材采用相反约定,所有公式都得跟着翻,抄公式时不能只抄一半。

另一个必须提前定的是角度单位。ROS 体系里很多接口用弧度,但示教器上显示的是度;配置文件里写的是度,代码里忘了转,末端姿态就会以 57.3 倍的错误比例歪掉。我的做法是变量名带单位后缀,theta_rad、theta_deg,转换只在输入输出边界做一次,中间层永远只认弧度。这条规则帮我省掉的调试时间,比任何一个算法优化都多。

2. 旋转矩阵:九个数字里的六条硬约束

2.1 方向余弦:把 B 的三根轴搬到 A 里量一遍

旋转矩阵最朴素的理解方式,就是"把 B 坐标系的三根单位轴,拿到 A 坐标系里逐个描述一遍,然后按列拼起来"。B 系的 X 轴在 A 系里是一个单位向量,写成三分量;Y 轴、Z 轴同理。把这三个列向量并排放成一个 3×3 矩阵,得到的就正好是 R_AB。

这样构造出来的矩阵,每个元素都是两个单位向量夹角的余弦,所以叫方向余弦矩阵。它有一个非常直观的语义:某一列就是对应坐标轴在参考系里的方向。调试的时候我经常直接打印旋转矩阵,然后看一眼第三列——它就是末端 Z 轴(通常对应工具朝向)在世界系里的方向,如果这一列和期望的接近垂直方向差别很大,那问题基本就在姿态上而不是位置上。

这个"按列读"的习惯在排查问题时效费比极高。比如机械臂末端本应该竖直向下对着桌面,那么工具 Z 轴在世界系里应该接近 [0, 0, -1](取决于你工具系的 Z 指向内还是外)。你只要看 R 的第三列是不是接近这个值,就知道姿态对不对,不需要去反解欧拉角。

2.2 正交与行列式:九个数字只有三个自由度

旋转矩阵必须满足两个条件:R^T R = I,且 det(R) = +1。第一条是正交性,保证三根轴互相垂直且长度为一;第二条排除了镜像反射(那种变换行列式是 -1,它把左手系变成右手系,不是刚体旋转)。

九个元素加上六个独立约束(正交性给出六个独立方程),剩下来正好三个自由度,和"姿态有三个自由度"这个基本事实吻合。这也解释了为什么用旋转矩阵做优化时特别难受:你要么把九个数当自由变量然后每步做正交化投影,要么干脆换参数化。工程上后者更常见——用四元数或轴角去做优化和插值,只在需要和硬件接口对接时才转成旋转矩阵。

数值上还有一件事必须留意:旋转矩阵连乘多了以后会漂移。浮点误差累积几千步,R^T R 与单位阵的偏差可能到 1e-6 量级甚至更大,行列式也会从 1 慢慢偏出去。这时候如果拿去求逆,误差会被放大。我的处理方式是在关键节点做一次正交化,最简单的是 SVD 法:

import numpy as np def orthonormalize(R): U, _, Vt = np.linalg.svd(R) Rn = U @ Vt if np.linalg.det(Rn) < 0: U[:, -1] *= -1 Rn = U @ Vt return Rn

如果只是在做视觉显示或者低频计算,这一步可以省;但如果是几百赫兹的迭代估计环路,建议每步或者每几十步做一次,成本很低,稳定性提升明显。

2.3 复合顺序:左乘和右乘差的是一整个坐标系

旋转的复合是姿态描述里最容易出错的地方。R_AC = R_AB · R_BC 这个式子要读成"先把 C 系里的点转到 B 系,再转到 A 系"。矩阵乘法从右往左作用在向量上,这个顺序不能乱。

左乘和右乘的差别,用一句话概括:左乘相当于绕参考系(固定系)旋转,右乘相当于绕当前坐标系(动系)旋转。举个例子,一个物体当前姿态是 R,你想让它绕自身 Z 轴转 30 度,写 R · Rz(30°);想让它绕世界 Z 轴转 30 度,写 Rz(30°) · R。这两个结果一般不同,只有在这两种 Z 轴恰好平行时才相等。

我见过的最典型的翻车场景是云台控制。操作手想让云台先俯仰再偏航,代码里按偏航矩阵乘俯仰矩阵的顺序写了,结果每次俯仰之后偏航的方向都跟着变了,操作手感极其别扭。改成左乘固定轴顺序之后立刻就对了。所以当有人问你"姿态控制不听话"的时候,第一个该问的问题往往不是算法,而是"你矩阵乘的顺序是什么"。

def rot_x(t): c, s = np.cos(t), np.sin(t) return np.array([[1, 0, 0], [0, c, -s], [0, s, c]]) def rot_y(t): c, s = np.cos(t), np.sin(t) return np.array([[c, 0, s], [0, 1, 0], [-s, 0, c]]) def rot_z(t): c, s = np.cos(t), np.sin(t) return np.array([[c, -s, 0], [s, c, 0], [0, 0, 1]]) # 绕世界 Z 转 90°,再绕世界 X 转 90° R1 = rot_x(np.pi/2) @ rot_z(np.pi/2) # 交换顺序 R2 = rot_z(np.pi/2) @ rot_x(np.pi/2) print(np.round(R1, 3)) print(np.round(R2, 3))

跑一下这两行,你会看到两个完全不同的矩阵,这是理解复合顺序最快的办法,比看十页推导都管用。

3. 欧拉角与 RPY:说人话的代价是死锁

3.1 十二种约定与"内旋外旋"这组绕人的词

欧拉角的核心思想是把一个任意姿态拆成三次绕坐标轴的旋转,每次一个角度。因为每次可以选择 X、Y、Z 三根轴中的任意一根,且相邻两次不能选同一根,组合起来一共有 12 种约定。这个数字本身就说明了问题的复杂度:如果你只写"欧拉角是 (30°, 45°, 60°)",这句话信息量约等于零,因为不知道约定就无法还原姿态。

这 12 种约定又被分成两大类。第一类是"绕动轴旋转",也叫内旋,每次旋转都绕着上一次旋转之后的新坐标轴;第二类是"绕定轴旋转",也叫外旋,每次旋转都绕原始参考系的固定轴。有个非常实用的等价关系:绕动轴按 X→Y→Z 顺序转,结果等于绕定轴按 Z→Y→X 顺序转。记住这条,看到不同教材的描述就能互相翻译。

我的建议是:在代码注释和接口文档里,永远把约定写全,格式写成"内旋 ZYZ"或者"定轴 XYZ(等价内旋 ZYX)",不要只写"欧拉角"。团队协作里因为省略约定造成的返工,我见过的次数两只手数不过来。

3.2 ZYZ 和 RPY 各自的地盘

机器人领域最常见的两种约定,一种用在机械臂腕部,一种用在飞行器和移动平台。

机械臂的球腕(三个旋转轴交于一点的那种手腕)几乎都用 ZYZ 欧拉角。原因是它的三个角有明确的物理含义:最外面那个角对应整条臂的方位,中间那个角对应腕部的倾斜(通常限制在 0 到 180 度之间),最里面那个角是末端绕自身轴的滚转。这种参数化的好处是中间角接近 0 或 180 度时腕部会出现奇异,而这个奇异恰好和机械臂运动学奇异对应得上,分析起来很整齐。

飞行器、移动机器人、视觉标定里则普遍用 RPY,也就是绕固定轴依次 X(roll)、Y(pitch)、Z(yaw),数学上写成 R = Rz(γ) · Ry(β) · Rx(α),注意乘法顺序是从右往左对应 X、Y、Z,写的时候别把顺序弄反。它的好处是三个角有很强的直觉:roll 是左右倾斜,pitch 是抬头低头,yaw 是左右转向。IMU 输出的姿态角基本都是这个格式。

从旋转矩阵反解 RPY 的公式值得记一下,因为实际项目里几乎天天用:

def R_to_rpy(R): # R = Rz(yaw) @ Ry(pitch) @ Rx(roll) sy = -R[2, 0] sy = np.clip(sy, -1.0, 1.0) pitch = np.arcsin(sy) if abs(abs(sy) - 1.0) < 1e-8: # 奇异位形,roll 与 yaw 不独立,约定 roll = 0 roll = 0.0 yaw = np.arctan2(-R[0, 1], R[1, 1]) else: roll = np.arctan2(R[2, 1], R[2, 2]) yaw = np.arctan2(R[1, 0], R[0, 0]) return roll, pitch, yaw

注意那个np.clip。如果旋转矩阵带一点数值误差,-R[2,0]可能变成 1.0000000002,arcsin会直接返回 nan,然后整个姿态估计链路就崩了。这个小防护我在每个反解函数里都加。

3.3 万向节死锁:几何图像加上数值验证

死锁这个词被讲得特别玄,其实几何图像很清楚。以 RPY 为例,roll 绕 X 转,pitch 绕 Y 转,yaw 绕 Z 转。当 pitch 等于正负 90 度时,Y 轴和 Z 轴发生了关系——第一次 roll 旋转会把 Z 轴转到原来的 X 位置附近,结果就是 roll 和 yaw 变成了绕同一根物理轴旋转,两个角度不再独立。三个自由度退化成了两个,你无法单独调整 roll 和 yaw。

注意死锁不是数值问题,也不是实现 bug,它是这套参数化方法固有的拓扑缺陷。欧拉角本质上是用三个数去覆盖一个拓扑上不是环面、也不是球面的流形,任何最小参数化都必然存在奇异,只是位置不同。所以"换一种欧拉角约定就能避开死锁"这种想法是不成立的,只能把奇点挪到你不常用的区域去。

数值上怎么验证自己的代码有没有正确处理死锁?我的做法很土但很有效:构造一个 pitch 恰好等于 90 度的矩阵,转成 RPY 再转回矩阵,对比两者。如果回来的是 nan 或者完全不对,说明缺了上面的分支处理。

def check_gimbal_lock(): R = rot_z(0.7) @ rot_y(np.pi/2) @ rot_x(0.3) r, p, y = R_to_rpy(R) R2 = rot_z(y) @ rot_y(p) @ rot_x(r) print("pitch =", p, " 重建误差 =", np.max(np.abs(R - R2)))

在死锁位形上,重建出来的矩阵和原始矩阵应该完全一致(差异在 1e-15 量级),但反解出的 roll 和 yaw 可能和输入值不一样——因为不唯一。如果你的单元测试断言"反解出的角度必须等于输入角度",它在死锁点必然会失败。这个测试用例的写法本身就是个坑,我建议改成断言"重建矩阵一致",鲁棒得多。

3.4 反解分支与角度连续性

除了死锁,欧拉角还有两个实际使用中的麻烦。

第一个是多解。同一个旋转矩阵,通常有两组欧拉角解,比如 RPY 里 (roll, pitch, yaw) 和 (roll+180°, 180°-pitch, yaw+180°) 都是合法解。反解函数默认返回哪一组,取决于 atan2 的分支和 pitch 的取值范围。做轨迹规划时,如果相邻两个路点的解不连续,中间就会插出一条莫名其妙的翻转路径,机械臂会突然抡半圈。解决办法是在反解之后做一次分支选择:把候选解都算出来,选和上一个时刻角度差最小的那一组。

第二个是 ±180 度附近的跳变。角度从 179 度走到 -179 度,数值上跳了 358 度,但物理上只转了 2 度。对角度做求导、滤波、PID 控制时,如果直接处理原始数值,就会在跳变点产生巨大的冲击。标准做法是做角度归一化,把差值折到 (-180°, 180°] 区间:

def wrap_pi(a): return (a + np.pi) % (2 * np.pi) - np.pi

这个函数短到不值得称为算法,但我在几乎所有角度相关的代码里都会定义它,而且只在"计算两个角度的差"时使用。越过一次 180 度导致的控制器震荡,排查起来非常费时间,因为现象是偶发的,数据看起来也"没错"。

从工程实践角度看,我更倾向的结论是:欧拉角和 RPY 适合做人的接口和局部表示,不适合做内部的连续状态量。内部状态一律用旋转矩阵或四元数,只在显示、日志、给操作员看的时候转成角度。这个原则帮我避免了大量和分支、跳变、死锁相关的杂事。

4. 轴角、四元数与旋转的最小表示

4.1 罗德里格斯公式:从一个轴和一个角度长出整张旋转矩阵

欧拉定理告诉我们,任何一个三维旋转都可以表示成绕某一根固定轴转某一个角度。于是姿态可以用四元数之外的另一种方式描述:一个单位向量 k 加上一个角度 θ,总共四个数、一个单位长度的约束,正好三个自由度。这个表示叫轴角,物理意义极其清晰,做视觉伺服和手眼标定时特别直观——误差就是"绕某个方向转多少度"。

从轴角到旋转矩阵的转换由罗德里格斯公式给出:

R = I + sinθ · [k]× + (1 - cosθ) · [k]ײ

其中 [k]× 是由向量 k 构成的反对称矩阵,也就是叉乘的矩阵形式。这个公式值得抄下来存着,因为在力控和阻抗控制里,力矩与姿态误差的关系直接和这个公式相关。它的一个直观解释是:第一项是原始方向,第二项控制小角度时的线性响应,第三项补足旋转的非线性部分。

def skew(k): return np.array([[0, -k[2], k[1]], [k[2], 0, -k[0]], [-k[1], k[0], 0]]) def axis_angle_to_R(k, theta): k = np.asarray(k, dtype=float) k = k / np.linalg.norm(k) K = skew(k) return np.eye(3) + np.sin(theta) * K + (1 - np.cos(theta)) * (K @ K) def R_to_axis_angle(R): theta = np.arccos(np.clip((np.trace(R) - 1) / 2, -1.0, 1.0)) if theta < 1e-8: return np.array([1.0, 0.0, 0.0]), 0.0 if abs(theta - np.pi) < 1e-6: # 180 度附近,反对称部分退化,改用对角线求解 A = (R + np.eye(3)) / 2 k = np.sqrt(np.clip(np.diag(A), 0, None)) idx = int(np.argmax(k)) if k[idx] > 1e-8: k = A[:, idx] / k[idx] return k / np.linalg.norm(k), theta k = np.array([R[2, 1] - R[1, 2], R[0, 2] - R[2, 0], R[1, 0] - R[0, 1]]) k = k / (2 * np.sin(theta)) return k, theta

需要重点提醒的是 180 度附近的分支。当 θ 接近 π 时,反对称部分 R - R^T 趋近于零矩阵,用 k = (R - R^T) / (2 sinθ) 求出来的轴会被噪声彻底淹没。这是我在做半个圆周翻转的关节标定时踩过的坑,最后的做法是切到对角线的平方根方法,就是上面代码里那一段。这类"在某一点附近特殊处理"的分支,每个姿态表示法里都有,写通用函数时必须考虑。

4.2 四元数的半角和它的双覆盖

四元数用四个数表示旋转,q = [w, x, y, z] 或者写成标量加向量的形式。单位四元数对应旋转,表示方式是 q = [cos(θ/2), sin(θ/2) · k],其中 k 和 θ 就是轴角表示里的轴和角度。注意那个二分之一角,这是四元数最关键也最容易被忽略的细节。

为什么是半角?因为四元数是通过旋转的"双覆盖"来作用的:单位四元数构成的空间 S³ 到 SO(3) 是一个二对一的映射,q 和 -q 表示同一个旋转。这个性质带来两个直接后果。第一,两个姿态之间的"角度距离"在四元数空间里是实际转角的一半,做误差计算时不要忘了乘 2。第二,四元数插值走的路程是实际旋转的一半弧长,所以天然的插值结果比线性插值平滑得多。

四元数的乘法在工程代码里天天用,公式不难但容易写错某一项的符号。我习惯直接用固定写法:

def qmul(q1, q2): 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, ])

写完一定要做回归测试:用 qmul 把同一个旋转连续乘两次,结果应该等于角度翻倍的同一根轴的旋转。这个小测试能在一分钟内抓出符号错误,比对着公式检查半天快得多。

另外,从四元数转旋转矩阵的时候,如果 q 没有归一化(比如多次乘法累积了浮点误差),转出来的矩阵会带一个整体缩放,把位置一起放大,现象是"机械臂好像变长了一点"。我一般在 qmul 之后立刻归一化,代价极小。

4.3 slerp:为什么姿态插值必须走这条路

两个姿态之间怎么平滑过渡,是最能体现表示法差异的地方。对欧拉角做线性插值,看起来简单,实际会出大问题:当两个姿态的欧拉角相差较大,或者经过死锁区域时,插值出来的路径会和"最短旋转路径"偏离很远,甚至中途出现奇怪的翻转。原因还是那个老问题——欧拉角是三个独立的数,线性插值等于在三个坐标轴上分别走直线,这和旋转流形上的测地线不是一回事。

四元数的球面线性插值(slerp)解决的正是这个问题。它的思路是:在单位球面上沿着两个四元数之间的最短大圆弧走,而不是在四维空间里走直线。公式不复杂,工程实现上要注意两个细节。

def slerp(q0, q1, t): q0 = q0 / np.linalg.norm(q0) q1 = q1 / np.linalg.norm(q1) d = float(np.dot(q0, q1)) if d < 0.0: q0 = -q0 d = -d if d > 0.9995: q = q0 + t * (q1 - q0) return q / np.linalg.norm(q) th0 = np.arccos(np.clip(d, -1.0, 1.0)) th = th0 * t q2 = q1 - q0 * d q2 = q2 / np.linalg.norm(q2) return q0 * np.cos(th) + q2 * np.sin(th)

第一个细节是if d < 0: q0 = -q0。因为 q 和 -q 表示同一个旋转,如果不做这一步,插值可能沿着"绕远路"的方向走,实际效果是末端突然转一个大圈。第二个细节是d > 0.9995时切到线性插值,因为此时 sinθ 接近零,除法会数值不稳定。这两个分支我在每一份姿态插值代码里都写,属于标准防护。

还有一个常被忽略的问题:slerp 给出的是姿态插值,但通常我们还要同时插值位置,两者用同一个 t 参数。这在大部分场景没问题,但如果是抓取这类高精度任务,建议对位置用较高次的多项式或者样条,对姿态用 slerp,然后把两者在时间上对齐,而不是简单地都用线性。

4.4 四种表示法到底该怎么选

我把常用的几种姿态表示放在一起做个对照,这张表我在项目里给新人讲过很多次。

表示法参数个数冗余优点缺点典型用途
旋转矩阵96 个约束复合简单、作用直观冗余大、迭代会漂移数值计算中间量、链式变换
欧拉角 / RPY3无人易读、直观有死锁、多解、跳变显示、日志、人机接口
轴角41 个约束物理意义清晰0 度和 180 度需分支误差定义、力控
四元数41 个约束插值好、复合高效不直观、双覆盖内部状态量、姿态估计

选型的基本逻辑是:内部计算用旋转矩阵和四元数,对外接口用人能看懂的角度。具体一点,运动学链式求变换用 4×4 齐次矩阵最方便;姿态估计滤波器的状态量用四元数;姿态误差和期望力矩用轴角;示教器显示和配置文件用 RPY。中间每一步转换都要写单元测试,因为这些转换函数一旦写错,错误会以"看起来对了但就是不精确"的形式潜伏很久。

5. 齐次变换:把位置和姿态打包进一次乘法

5.1 4×4 分块和齐次坐标存在的理由

有了旋转矩阵和位置向量,一个坐标系相对另一个坐标系的完整描述就是一个 4×4 矩阵:

T_AB = [[R_AB, p_AB], [0, 0, 0, 1]]

左上 3×3 是姿态,右上 3×1 是 B 系原点在 A 系中的位置,最后一行固定是 [0, 0, 0, 1]。这个结构叫齐次变换矩阵,它把"旋转加平移"这两件事统一成了一次矩阵乘法。

齐次坐标的引入乍看像是为了凑维度,实际上解决了两个根本问题。第一,平移在三维里是加法,旋转是乘法,两者本来不能统一成一次线性运算,扩到四维之后都能写成乘法。第二,它顺便表达了"点和向量的区别"——点的第四分量是 1,向量的第四分量是 0。这个区别非常实用:一个向量经过变换矩阵时不应该被平移,而一个点应该被平移。

我强烈建议在代码里显式区分这两类数据。比如定义一个函数transform_point和一个transform_vector,内部用不同的齐次分量,不要图省事都用 1。把方向向量当点来变换,是所有变换库里最隐蔽的错误之一,因为它在平移量很小的时候几乎看不出来。

5.2 链式相乘和求逆的解析式

齐次变换最大的价值在于链式关系:T_AC = T_AB · T_BC。机械臂的正运动学就是把从基座到末端的每一段变换乘起来;手眼标定就是把相机到末端的变换和基座到相机的变换串起来。整个过程是纯矩阵乘法,写起来干净利落,这也是为什么几乎所有机器人库都用齐次变换作为标准接口。

更值一提的是求逆。数值求 4×4 矩阵的逆是通用操作,但齐次变换的逆有解析形式,而且比数值求逆快得多、稳得多:

T_BA = T_AB^(-1) = [[R_AB^T, -R_AB^T · p_AB], [0, 0, 0, 1]]

这个公式的推导很简单:逆变换的姿态就是转置,位置则是把原点变换回去。它的实现只要几次矩阵乘法,而且天然保持正交性,不会像数值求逆那样引入误差。

def inv_transform(T): R = T[:3, :3] p = T[:3, 3] Ti = np.eye(4) Ti[:3, :3] = R.T Ti[:3, 3] = -R.T @ p return Ti

这段代码我在每个项目里都会写一遍。用np.linalg.inv(T)也能跑,但在几百赫兹的环路里,解析形式的性能优势和数值稳定性优势都很实在。而且如果 T 因为浮点误差而不再是严格的齐次变换,通用求逆会把这个误差原封不动放大,解析式则不会。

5.3 三个高频混用的错误

第一个错误是变换顺序反了。T_AB · T_BC 和 T_BC · T_AB 完全不是一回事。判断方法很简单:看下标能不能"约掉",能约掉的就是正确顺序,就像单位换算一样。这个"下标相消"的检查法,我在写多级变换链的时候每次都会用,只要下标接不上,基本就能断定顺序错了。

第二个错误是齐次矩阵没有保持结构。比如做插值时对四个块独立插值,插完之后左上角不再是正交矩阵,导致缩放变形。正确做法是把旋转和平移分开插值(旋转用 slerp),插完再拼回齐次矩阵。

第三个错误是混用坐标系。手眼标定里最容易出现,因为涉及基座、末端、相机、标定板四个坐标系。我见过一个案例,工程师把标定板的位姿从相机系直接当成了基座系下的目标位姿,误差看起来很小(因为相机离标定板近),但机械臂怎么调都差几毫米。最后是在每个变换上强制加命名,用静态检查工具禁止无命名的变换相乘才解决的。

6. DH 参数:让"描述"变成机械臂真能到的坐标

6.1 四个参数各自管什么

前面讲的都是"两个坐标系之间怎么描述",而机械臂正运动学的任务是"给一组关节角,末端在哪"。中间需要的是一座桥,这座桥就是 Denavit-Hartenberg 参数。

标准 DH 约定下,相邻两个连杆坐标系之间的变换由四个参数决定:连杆长度 a_i(沿前一关节轴方向的偏移)、连杆扭转角 alpha_i(两根相邻关节轴的夹角)、连杆偏距 d_i(沿关节轴方向的偏移)、关节角 theta_i(绕关节轴的转角)。对于旋转关节,theta_i 是变量,其他三个是常值;对于移动关节,d_i 是变量。这个分工很干净,所以 DH 参数一直是教科书的标配。

单个连杆的变换可以写成四个基本变换的连乘,顺序固定:

A_i = Rz(theta_i) · Tz(d_i) · Tx(a_i) · Rx(alpha_i)

展开之后就是那个非常标准的 4×4 矩阵:

def dh_transform(a, alpha, d, theta): ct, st = np.cos(theta), np.sin(theta) ca, sa = np.cos(alpha), np.sin(alpha) return np.array([ [ct, -st*ca, st*sa, a*ct], [st, ct*ca, -ct*sa, a*st], [0.0, sa, ca, d], [0.0, 0.0, 0.0, 1.0], ])

这段代码短到可以背下来,但它的正确性完全依赖于前面的约定。一旦换到别的约定,矩阵形式就变了。

6.2 标准 DH 与改进 DH 的账要算清

这里必须说清楚一件事:DH 参数有两套主流约定,标准 DH 和改进 DH(也叫 Craig 约定),它们把坐标系固连的位置不一样。标准 DH 把坐标系固连在连杆的远端(靠近下一个关节),改进 DH 固连在近端(靠近上一个关节)。结果是同一个机械臂,两套约定给出的参数表完全不同,变换矩阵的乘法顺序也不同。

改进 DH 的顺序是:A_i = Tx(a_{i-1}) · Rx(alpha_{i-1}) · Rz(theta_i) · Tz(d_i)。

def mdh_transform(a_prev, alpha_prev, d, theta): ct, st = np.cos(theta), np.sin(theta) ca, sa = np.cos(alpha_prev), np.sin(alpha_prev) return np.array([ [ct, -st, 0.0, a_prev], [st*ca, ct*ca, -sa, -d*sa], [st*sa, ct*sa, ca, d*ca], [0.0, 0.0, 0.0, 1.0], ])

我踩过的坑是:从厂商手册抄参数表,手册用的是改进 DH,代码里写的是标准 DH 的变换函数,结果末端位置差了十几厘米。查了半天坐标变换,最后发现是两套约定混用。从那以后我定了一条规矩:任何一个 DH 相关的文件,第一行注释必须写明用的是哪套约定,且不允许在一个项目里同时出现两套,除非有明确的转换说明。

6.3 手写正运动学并做闭环自检

把参数表和变换函数凑在一起,正运动学就是一行循环:

def forward_kinematics(dh_params, joint_angles): T = np.eye(4) for (a, alpha, d, theta_offset), q in zip(dh_params, joint_angles): T = T @ dh_transform(a, alpha, d, q + theta_offset) return T

跑通之后最重要的一件事是自检。我常用的自检方法有三种。

第一种是零位检查。把所有关节角置零,看末端位姿是否和机械臂手册上标注的零位一致。这一步能抓出大部分符号错误和偏置遗漏,因为零位是已知的。

第二种是单轴旋转检查。只让第一个关节转 90 度,观察末端位置的变化是否符合几何直觉。比如基座关节转 90 度,末端应该绕基座 Z 轴画弧,高度不变。如果高度变了,说明有轴的对应关系搞错了。

第三种是数值微分检查。用有限差分算末端的线速度和角速度,和雅可比矩阵的解析结果对比。这两者在不奇异的位置应该吻合到 1e-6 量级。这个方法稍复杂,但它能抓出那种"位置看着对、姿态旋向反了"的隐蔽错误,我建议在正式使用前至少做一次。

提示:正运动学自检最好写成固定的测试脚本,每次改完参数就重跑一遍。我见过太多因为改了 DH 表里一个小数点导致整台设备行为异常的情况,而这种问题靠肉眼翻参数表基本不可能发现。

7. 姿态出错时的排查链路

7.1 先查约定,再查数据

姿态问题排查有个原则:不要在还不确定约定是否统一的时候去分析数值。所以第一步永远是核对约定链条。我会按顺序问自己四个问题:输入数据的姿态用的是什么表示法?它的旋转顺序是内旋还是外旋?角度单位是度还是弧度?它描述的是"坐标系到坐标系"还是"物体到世界"?

这四个问题里任何一个不清楚,后面所有的计算都建立在不牢靠的假设上。实际项目里最省事的做法是找一组已知答案的数据做端到端验证。比如把相机正对一块标定板,此时标定板相对相机应该有一个接近单位阵的旋转矩阵。如果算出来不是,问题就在链条中的某个环节,而不是在某个公式上。

第二步才是查数值。查数值也有顺序:先看旋转矩阵是否正交(把 R^T R 打印出来,非对角元应该在 1e-6 以内),再看行列式是否为 +1,再看是否有明显的数量级异常(比如某一行全是零,那多半是数据没读进来)。这三步一分钟就能做完,却能筛掉一大半低级问题。

7.2 可视化:别用肉眼看九宫格

九个数一屏,人眼是分辨不出对错的。姿态调试一定要可视化。我常用的三种手段,成本从低到高。

第一种是打印三根轴的方向向量。旋转矩阵的每一列就是一根轴,直接看第三列是否指向期望方向,比看整张矩阵快得多。如果工具 Z 轴应该朝下而第三列接近 [0, 0, 1],那就是差了 180 度,问题可能出在轴的定义上。

第二种是在三维里画坐标轴。用一个简单的坐标系可视化工具,把每个关键坐标系的三根轴画出来,箭头方向一目了然。多关节情况下配合滑块调角度,很快就能定位到是哪一段变换错了。这个手段在我定位手眼标定问题时效率最高。

第三种是把姿态角随时间的变化画成曲线。如果在某个时刻曲线上出现尖刺,那基本就是 ±180 度跳变或者分支选择出问题。这类问题只看单帧数据永远发现不了。

我更倾向的做法是:姿态相关的代码一旦改过,先把曲线画出来看一遍再上真机。哪怕只是改了一个符号,也值得多花两分钟。

7.3 一张高频坑位对照表

现象最可能的原因快速验证方式
末端位置对、姿态整体差 180 度轴的定义相反(例如 Z 指向内/外)检查旋转矩阵第三列符号
姿态在小角度对,大角度偏欧拉角顺序或内旋外旋搞反用 90 度级别的组合旋转验证
姿态角曲线出现尖刺±180 度跳变未归一化打印连续两帧的角度差
反解出现 nanarcsin/arccos 输入越界加 clip,检查死锁位形
姿态缓慢漂移迭代连乘未做正交化打印 R^T R - I 的范数
插值过程突然翻大圈四元数符号未统一检查相邻帧点积是否为负
末端位置整体偏一点向量被当成点做了平移检查齐次分量用的是 0 还是 1

这张表里的每一条我都在真实项目里遇到过,其中"姿态缓慢漂移"和"插值翻大圈"这两条花的时间最长,因为它们的现象不具备偶发性之外的特征,数据看起来都没问题。

8. 我自己用下来比较稳的几个习惯

最后说说几个长期用下来觉得值得坚持的做法。

第一个是"边界转换"原则:整个系统内部只保留一种姿态表示。我的默认组合是齐次矩阵做链式变换、四元数做状态量、RPY 只用于显示。所有转换函数集中在同一个模块里,每个函数都配单元测试。这样一来,当出现问题时,可怀疑的范围就压缩到了少数几个函数上,而不是散落各处。

第二个是"数据自带坐标系"原则。任何位置和姿态变量,要么名字里带坐标系,要么封装成带参考系标签的结构体。这条规则在只有一个人的项目里显得啰嗦,但只要超过两个人协作,回报立刻显现。我现在的习惯是连日志打印都带坐标系名,回看几十条日志时能省下大量猜测。

第三个是"用重建法验证反解"。任何一个从姿态反解角度的函数,都要配一个"反解再重建"的测试,断言重建出的旋转矩阵与原矩阵一致,而不是断言角度相等。这个小改动让我的测试用例在死锁位形和 180 度附近都能稳定通过,不再需要为特殊情况写例外。

第四个是保留一组"已知答案"的基准数据。我手上有一份从实际设备采出来的姿态序列,包含小角度、大角度、接近死锁、跨越 ±180 度这几类。每次改动姿态相关的代码,先拿这份数据跑一遍,看输出是否和上次一致。这比重新设计测试用例省事得多,也更容易发现"改了一个地方、影响了另一个地方"的回归问题。

如果后面还要往下扩展,我建议的下一步是看看姿态误差在优化问题里怎么定义。因为一旦涉及轨迹优化或者标定求解,误差函数写成 R1 - R2 的范数是不对的,正确的写法是 log(R1^T R2) 这类基于李代数的方式。那个话题比姿态描述本身要绕一些,但理解了今天这些表示法的来龙去脉之后,再看它就不会觉得突兀了。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询