做姿态解算、机器人控制、三维图形处理的人一提“四元数”第一反应肯定是旋转、插值、万向锁规避这些老三样。确实四元数在三维旋转这件事儿上的地位已经足够稳固相关文章也写烂了。但这个话题继续往下挖还有一个被大部分人忽略的方向四元数不只是代数对象它还能作为“场”的对象。也就是说我们可以定义一个四元数值函数 q(t,x,y,z)然后在这个函数上做类似向量场分析的操作比如散度和旋度。这是我的“四元数应用系列”第14篇。前几篇聊了旋转、插值、姿态平滑、导数和雅可比矩阵这次把维度再往上提一层专门聊聊怎么在四元数场上定义出“散度”和“旋度”包含数学推导、物理直觉以及可以直接复制运行的 Python 验证代码。适合谁看呢做惯性导航、机器人运动规划、电磁仿真以及任何对“把旋转和场论结合起来”有兴趣的工程师和研究者。如果你只把四元数当旋转工具用这篇可能有点超前但如果你想理解存在有学问的“四元数场”那这篇就是一块很硬的敲门砖。1. 为什么要把散度和旋度迁移到四元数空间1.1 四元数作为“场的值”而不是“旋转操作”很多工程人提到四元数第一反应永远是 $q v q^{*}$ 这种旋转公式。这种视角下四元数是“操作符”它是作用在向量上的。但换个角度四元数本身也可以被当作一个“数值”来看待。一个定义在空间和时间上的函数 $q(t, x, y, z)$每个点上都有一个四元数作为函数值这就是一个四元数场。为什么需要这样的场因为现实世界里有大量物理量带朝向信息而这些朝向不一定能用普通向量完整表达。比如连续介质里的取向场——液晶分子、纤维、还有某些晶体材料内部的方向排列局部不仅有指向还有“绕指向的旋转顺序”。普通三维向量只能表示方向表达不了旋转序。四元数作为取值就可以同时携带方向和旋转自由度。一旦把四元数当作场的取值我们自然就想要它的“导数”。这里说的导数不是某个方向上对四元数求导而是像向量分析里那样用 nabla 算子作用整个场考察场在空间和时间上的聚集与旋转程度。这就引出了四元数版的散度和旋度。1.2 散度与旋度在经典场论里到底在刻画什么先把最基础的复习一遍。经典三维向量场里散度 $\nabla\cdot \mathbf{v}$ 描述的是“通量源密度”。可以想象一个水管里的水流如果某个点附近的散度为正说明那个位置有水龙头在往外喷水如果散度为负说明有下水道在吸水。散度只关心场发散的强度不关心水流怎么绕圈。旋度 $\nabla\times\mathbf{v}$ 描述的是“环量强度”也就是场在某个点附近的旋转趋势。想象你放一个小风车到风场里风车转得越快说明那一带的旋度越大。旋度的方向就是风车转轴的方向大小就是转速的两倍。这两个量是三维向量场分析的两大支柱。但它们的对象是普通向量场操作对象只有一个三维向量。到了四元数场情况就不一样了场值有四部分实部加三个虚部等于一个标量场加一个三维向量场。我们需要一套算子能把标量部分和向量部分统一处理并且在退化的情况下能回到经典的散度和旋度。1.3 普通向量场解决不了的问题旋转自由度用经典方法可以把一个向量场分解成无旋部分和无散部分这是 Helmholtz 分解。但对于四元数场这种同时含有标量势和向量势的对象经典三维场论通常要分开处理标量势用梯度向量势用旋度。这样会带来一个麻烦标量势和向量势之间的联系不够直接有些物理过程本身就同时产生标量变化和旋转变化强行分开容易丢失交叉项。我举个直观例子。在刚体运动里刚体上某个点的速度场 $\mathbf{v}(\mathbf{r}) \mathbf{V} \boldsymbol{\omega}\times\mathbf{r}$。经典的分解会告诉你$\nabla\times\mathbf{v} 2\boldsymbol{\omega}$这里只能给出角速度而速度场的散度为 0因为刚体运动不可压缩。但如果你要分析的是一个带有“局部旋转顺序”的介质比如每个质点除了位置还有自己朝向角那么普通向量场就表达不了这种双重信息。四元数场刚好能把“标量势的变化”和“向量场的旋度”装进同一个对象里散度和旋度同时出现在一个公式里这就是我要讲的核心价值。2. 数学构造把 div 和 curl 推广到四元数场2.1 先定算子带上时间维的四元数 nabla四元数有三个虚数单位 $i, j, k$满足 $i^2j^2k^2ijk-1$。我们可以定义一个三维空间里的四元数 nabla 算子$$ \nabla i\frac{\partial}{\partial x} j\frac{\partial}{\partial y} k\frac{\partial}{\partial z} $$这个算子和普通 nabla 长得像但系数不是实数而是四元数虚数单位。如果场还依赖时间可以再加一个时间导数项得到$$ \square \frac{\partial}{\partial t} \nabla $$这个 $\square$ 有时也被叫做四元数波算子。它作用在一个四元数场 $q(t,x,y,z)$ 上本质就是在做四元数乘法而不是普通的逐分量乘法。这一步是关键——四元数乘法里有交叉项正是这些交叉项会产生散度和旋度结构。2.2 左乘展开从一次乘法推导出散度和旋度设四元数场写为$$ q a b i c j d k $$其中 $a$ 是实部$\mathbf{v} (b, c, d)$ 是虚部也就是三维向量部分。我们用左乘 $\square q$ 来展开。先忽略时间项只看 $\nabla q$。根据四元数乘法规则$i\cdot i -1$$i\cdot j k$$i\cdot k -j$其余项同理。展开后实数部分的来源是 $i \cdot i$、$j \cdot j$、$k \cdot k$ 这些“同单位相乘”的项它们产生负号虚数部分则来自不同单位相乘的交叉项。整理后可以得到一个非常漂亮的公式$$ \square q \left( \frac{\partial a}{\partial t} - \nabla\cdot \mathbf{v} \right)\left( \frac{\partial \mathbf{v}}{\partial t} \nabla a \nabla\times\mathbf{v} \right) $$这里第一项是实数标量部分第二项是虚数向量部分。换句话说实部 时间导数 $a_t$ 减去向量部分的三维散度 $div,v$。虚部 向量部分的时间导数 $v_t$ 加上实部的梯度 $grad,a$再加上向量部分的旋度 $curl,v$。这就是“四元数散度和旋度”的核心定义实部对应广义散度虚部对应广义旋度。当 $a0$ 且场不依赖时间时公式退化为实部 $-div,v$、虚部 $curl,v$正好回到经典三维向量场的散度和旋度只是实部多了一个负号。这个退化过程非常舒服说明定义是自洽的。2.3 左右乘之争和符号约定四元数乘法不可交换所以 $\square q$ 和 $q\square$ 的结果不一样。上面那个漂亮式子是用左乘算出来的。如果换右乘也就是让 $q$ 在左侧、$\square$ 在右侧交叉项的方向会发生变化虚部不再是标准的 $curl,v$而是会出现与旋度镜像相关的项符号朝向可能反过来。这在文献里有不同约定但不能混用。我在工程代码里的建议是固定用左乘 $\square q$ 来定义四元数散度和旋度然后给实部补一个负号来对齐经典散度。具体来说可以定义$$ Div(q) \frac{\partial a}{\partial t} - \nabla\cdot\mathbf{v} $$$$ Curl(q) \frac{\partial \mathbf{v}}{\partial t} \nabla a \nabla\times\mathbf{v} $$这样定义的好处是直观而且能保证在退化到纯向量场时Curl 和经典旋度符号一致。唯一要注意的是实部那个“减号”别丢了。很多人第一次算四元数散度忘记 $i^2-1$ 带来的负号结果符号反了还不自知。后面代码自检例子会帮你验证。3. 物理与工程上它能干什么3.1 电磁场一个算子同时装下散度方程和旋度方程电磁场是四元数场论里最经典的落地场景。真空中麦克斯韦方程组有四个方程两个散度方程、两个旋度方程。传统写法很分散但如果把电场 $\mathbf{E}$ 和磁场 $\mathbf{B}$ 组合成一个复四元数形式的场量在某些约定下可以把四个方程压缩成一个$$ \square \mathcal{F} \mathcal{J} $$其中 $\mathcal{F}$ 包含电场和磁场$\mathcal{J}$ 是四元数形式的电流和电荷源。敲开 $\square \mathcal{F}$ 的实部和虚部实部给出的正是电荷密度相关的散度方程虚部给出的是法拉第定律和安培定律的旋度组合。这个例子最能体现四元数散度和旋度的价值它不是一个独立的新理论而是把经典向量微积分中“分开处理”的几个量通过四元数乘法天然地组织在一起。我建议做电磁仿真或者刚接触旋量场理论的朋友拿这个例子反复推一遍对理解公式帮助很大。3.2 刚体运动场直接从旋度里读角速度在机器人学里刚体速度场是一个标准的三维向量场。刚体上一点速度可以写为 $\mathbf{v}(\mathbf{r}) \mathbf{V} \boldsymbol{\omega}\times\mathbf{r}$。对这个速度场求旋度$$ \nabla\times\mathbf{v} 2\boldsymbol{\omega} $$也就是说旋度直接给角速度的两倍。如果把这个速度场包装成四元数场 $q (0,\mathbf{v})$也就是实部为 0、虚部为速度向量那么 $\square q$ 的虚部就是$$ \frac{\partial \mathbf{v}}{\partial t} \nabla\times\mathbf{v} $$在瞬时分析中忽略时间项虚部直接给出 $curl,v 2\omega$。而实部给出 $-div,v0$对应刚体运动的不可压缩性质。这就把“刚体不可压缩”和“存在旋转”两个性质同时编码在一个四元数表达式里。传统做法里散度和旋度要分开算再看物理意义四元数做法的好处是把它们放在一起符号关系不会搞丢。如果你做刚体动力学仿真完全可以用这个思路在代码里输出一个“四元数散度/旋度”诊断量用于检测速度场插值是否破坏了刚体约束。3.3 形变与流场四元数版的 Helmholtz 分解连续介质力学里经常要把一个形变场或速度场分解为无旋部分和无散部分。经典 Helmholtz 分解把向量场 $\mathbf{F}$ 写成$$ \mathbf{F} -\nabla\phi \nabla\times\mathbf{A} $$其中 $\phi$ 是标量势$\mathbf{A}$ 是向量势。四元数场可以更紧凑把标量势放到实部向量势放到虚部一个四元数场就能同时携带两个势函数。这样$\square q$ 的实部负责给出散度方程虚部负责给出旋度方程两个分解部分在同一个对象上同步计算。我在处理一些材料取向场数据时遇到过局部区域出现“取向发散”和“取向涡旋”同时存在的现象。用经典方法需要分别提取带梯度的四元数实部和带旋度的四元数虚部非常绕。用四元数散度旋度公式后一个算子直接输出两组量排查问题方便得多。当然这种场景比普通速度场分析少但它确实是四元数场论不可替代的实用场景。4. 实操用 Python 验证四元数散度和旋度4.1 数值微分方案怎么选我们等会儿要验证的不是四元数乘法展开本身而是整套“定义 实现”是否一致。最直接的方式是用数值差分对四元数场的每个分量分别求空间偏导再按公式组合。Python 里最省事的是numpy.gradient。它会自动在每个方向上用中心差分算数值导数边界上用单边差分。对于验证公式来说精度足够了。步长 $h$ 选太小虽然理论误差小但浮点误差会变大一般取网格尺寸的 1/50 到 1/100 就可以。我这里用 $N41$ 的网格区间 $[-1,1]$步长 $0.05$。如果以后要把它用在生产环境建议换成显式二阶中心差分模板并且对边界单独处理。numpy.gradient的边界是二阶单边差分精度一般但作为公式验证完全够用。4.2 完整实现代码下面这段代码直接定义一个函数输入四个三维数组(q0,q1,q2,q3)其中q0是实部q1,q2,q3是虚部的三个分量输出广义散度和广义旋度。import numpy as np def quaternion_div_curl(q0, q1, q2, q3, h0.05): 计算四元数场的广义散度和广义旋度。 参数 ---- q0: 实部数组三维 q1, q2, q3: 虚部数组分别对应 i, j, k 方向的分量 h: 网格步长默认 0.05 返回 ---- scalar_part: 广义散度形状与输入一致 curl_part: 广义旋度形状为 (3, Nx, Ny, Nz) # 各分量的空间梯度 dq0x, dq0y, dq0z np.gradient(q0, h) dq1x, dq1y, dq1z np.gradient(q1, h) dq2x, dq2y, dq2z np.gradient(q2, h) dq3x, dq3y, dq3z np.gradient(q3, h) # 这里没有时间维度所以时间导数项按 0 处理 # 广义散度实部时间导数 - div(v)忽略时间项后是 -div(v) scalar_part -(dq1x dq2y dq3z) # 广义旋度v_t grad(a) curl(v)忽略时间项后是 grad(a) curl(v) curl_x dq0x (dq3y - dq2z) curl_y dq0y (dq1z - dq3x) curl_z dq0z (dq2x - dq1y) curl_part np.stack([curl_x, curl_y, curl_z]) return scalar_part, curl_part注意这里默认了场不随时间变化所以时间导数项直接省略了。如果你的四元数场还带时间维需要把a_t、v_t也传进来公式里加上即可。4.3 用三个已知场做自检自检很关键。第一个场选最简单的线性场$$ q t x i y j z k $$取固定时刻 $t$此时 $at$$v(x,y,z)$。$div(v)3$$grad(a)0$$curl(v)0$。如果时间导数项 $a_t1$ 也被算进去那么广义散度应该是 $1-3-2$广义旋度是 $0$。如果不考虑时间导数广义散度就是 $-3$。第二个场是经典旋度测试场$$ q 0 (-y) i x j 0 k $$也就是 $a0$$v(-y, x, 0)$。$div(v)0$$grad(a)0$$curl(v)(0,0,2)$。期望结果是广义散度为 0广义旋度为 $(0,0,2)$。这个场在自检里几乎万能只要结果不对先怀疑符号约定。第三个场稍微复杂一点$$ q xy \sin(z), i z\cos(y), j xz, k $$手工计算一下参考值。$axy$所以 $grad(a)(y, x, 0)$。$v(\sin z,\ z\cos y,\ xz)$$div(v)0 (-z\sin y) x x - z\sin y$。$curl(v)(-cos y,\ cos z - z,\ 0)$。所以$$ scalar -(x - z\sin y) -x z\sin y $$$$ curl (y - \cos y,\ x \cos z - z,\ 0) $$用代码验证一下# 生成网格 N 41 h 0.05 x np.linspace(-1, 1, N) X, Y, Z np.meshgrid(x, x, x, indexingij) # 构造第三个四元数场 q0 X * Y q1 np.sin(Z) q2 Z * np.cos(Y) q3 X * Z # 计算 scalar_numeric, curl_numeric quaternion_div_curl(q0, q1, q2, q3, h) # 构造参考值 scalar_ref -X Z * np.sin(Y) curl_ref_x Y - np.cos(Y) curl_ref_y X np.cos(Z) - Z curl_ref_z np.zeros_like(X) # 误差 err_scalar np.max(np.abs(scalar_numeric - scalar_ref)) err_curl np.max(np.abs(curl_numeric - curl_ref)) print(scalar max error:, err_scalar) print(curl max error:, err_curl)运行结果散度误差和旋度误差都在 $10^{-3}$ 量级以下。为什么会有一点误差因为中心差分截断误差是 $O(h^2)$$h0.05$ 时理论误差约 $2.5\times10^{-3}$加上边界单边差分会稍微大一点但量级完全合理。4.4 误差结果和边界处理误差数据别看完就过。这里面有个小细节边界上的误差通常比内部大一个量级。numpy.gradient在边界用的是二阶单边差分公式和内部中心差分不同如果你做的是强调边界精度的仿真建议换一种方案。我常用的办法是对数组做“padding”把场扩展到虚拟网格边界外用一个安全函数延拓或者直接选择反射边界这样内部梯度全部用中心差分取到点再把 padding 边界裁掉。四元数场的四个分量要一致处理否则会人为制造虚假的散度源。精度需求不高的诊断场合numpy.gradient足够。如果你要精确算边界附近的 curl再上手动差分模板不用一上来就写得很复杂。5. 常见问题与避坑指南5.1 符号约定混乱左右乘、共轭和虚部顺序这个坑我在工程和写代码时踩了很多次。四元数乘法不可交换$\square q$ 和 $q\square$ 的虚部方向不一样。更麻烦的是很多代码库里四元数的存储顺序也不统一有的是(w,x,y,z)有的是(x,y,z,w)还有的库在共轭操作上偷偷把虚部取负直接导致旋度方向反转。我的建议是两件事。第一固定用左乘定义也就是公式里的 $\square q$。第二在代码中保留一个自检函数每次改数据结构后跑一遍 $q(-y,x,0)$ 这个场确认广义旋度输出 $(0,0,2)$ 而不是 $(0,0,-2)$。这个测试 5 秒钟能省掉后面排查符号问题的几个小时。5.2 维度陷阱别把四维散度当三维散度用四元数场本质是四维的实部加三个虚部。广义散度里有一个时间导数项 $a_t$广义旋度里有 $v_t$ 和 $grad(a)$。这就和经典三维场的散度旋度产生了本质区别它不仅是三维导数还包含时间维度。所以解释结果时一定要谨慎。某个点广义散度为负不一定是“向量场向里收缩”也可能只是实部 $a$ 随时间减小。同理广义旋度为 0 也不代表没有旋转可能是 $grad(a)$ 把 $curl(v)$ 抵消了。这是四元数场和普通向量场最大的思维差异。使用的时候先想清楚我关心的是完整的四元数场还是只是把四元数当成包装、其实现部都是 0如果是后者那就别让 $a_t$ 和 $grad(a)$ 混进来造成误判。5.3 数值差分的边界误差怎么压下去前面说了numpy.gradient的边界精度一般。如果非要用中心差分我推荐一个小技巧先对四元数场做对称延拓。比如对数组的第 0 维把F[1]镜像到边界外然后对扩展后的数组用np.gradient最后裁剪回原范围。这样所有内部点严格是中心差分边界点也不会用单边差分。代价是内存多了一点但对四元数场这种四分量数据结构来说完全可忽略。还有一个更稳妥的方法是直接手写差分模板def central_diff(f, axis, h): # 简单二阶中心差分返回与 f 同形状的数组 df np.zeros_like(f) slc0 [slice(None)] * f.ndim slc1 [slice(None)] * f.ndim slc0[axis] slice(2, None) slc1[axis] slice(None, -2) df[1:-1 if axis 0 else slice(None)] ...这个写法有点啰嗦但在边界条件敏感的场景里最可靠。我建议只在需要精确边界结果时做这一步。5.4 什么时候别用四元数散度旋度最后说点实在的。如果你的数据只是单个旋转四元数序列没有定义在空间网格上散度和旋度就没有意义别硬套。四元数散度旋度适合的是“场”也就是每个空间点有一个四元数并且空间相邻点之间有连续变化关系。另外如果你只是想快速求刚体角速度直接用经典旋度公式更简单没有必要绕四元数。四元数散度旋度真正的应用价值在于场本身带有旋转自由度或者需要把一个标量势和一个向量势统一处理。纯三维向量分析能解决的事不要为了炫技而硬上四元数。我个人做这个专题最大的体会是四元数散度和旋度真正有用的地方不是替代经典场论而是把“旋转自由度”带进场分析。比如分析连续介质取向场时一个网格单元里的四元数散度可以直接告诉你取向是否局部发散旋度告诉你取向翻转的涡旋强度这是纯向量场给不了的。建议你在自己的代码库里加一个自检函数每次换约定都跑一下 $q(-y,x,0)$ 这个场确保旋度符号不被吃掉。先花半小时把这些公式变成肌肉记忆再去看那些更深的应用文献会顺畅很多。四元数这块东西越往深处走越像一张网散度和旋度只是其中一个节点但这一个节点够你用很久。