从局部导数看场 · 第 2 章
散度与旋度:从雅可比中提取什么
向量场在一点附近的行为由雅可比完整描述,散度和旋度只是从中提取的两部分:散度是迹,读出局部体积变化;旋度来自反对称部分,等于局部角速度的两倍。本章推导这两点、它们的积分形式与零恒等式,并用光流散度估计碰撞时间。
系列目录第 2 / 3 章
问题是什么
一个向量场在某一点附近”在做什么”?散度 (Divergence) 和旋度 (Curl) 各自回答了其中哪一部分,又漏掉了什么?
教材通常直接给出 和 的分量公式,再配上”源”和”涡旋”的图。这样容易记住公式,却看不出它们为什么长这个样子。本章从第 1 章的雅可比出发,说明散度和旋度都是从雅可比中提取出来的信息。
直觉怎么理解
把向量场 想成流体的速度场,在某一点放一个很小的流体团,观察它在一小段时间内的变化。这个变化可以拆成四种:
- 平移:整体随速度 移动;
- 均匀膨胀或压缩:体积变大或变小,形状不变;
- 旋转:像刚体一样转动;
- 保体积形变:一个方向拉长、另一个方向压扁,体积不变。
平移由场的值决定,后三种由雅可比决定。散度读出的是第二种,旋度读出的是第三种。第四种两者都读不到。
下图是四个典型的二维场。前两个分别是纯膨胀和纯旋转;第三个是简单剪切,它同时包含旋转和形变;第四个是纯形变,散度和旋度都为零,但它显然不是一个平凡的场。
数学 / 几何表达
记 ,下标表示偏导,例如 。雅可比为
散度是雅可比的迹
为什么迹对应体积变化。 速度场在时间 内把点 移到 。相对于流体团中心,小位移 被映射为
小体积元的体积缩放比例是这个线性映射的行列式。按 Leibniz 公式, 是所有置换对应乘积之和:
- 恒等置换给出对角线乘积 ;
- 其余每个置换至少包含两个非对角元素,每个都带一个 ,所以贡献为 。
于是
散度是局部体积的瞬时相对变化率:正值表示膨胀,负值表示压缩,零表示体积一阶不变。散度为零既不意味着速度为零,也不意味着形状不变。
注意区分两个容易混淆的量: 描述一个空间映射本身的局部体积缩放, 描述速度场造成的体积变化率。两者通过上面的短时间映射 联系起来。
旋度来自雅可比的反对称部分
每个分量都是雅可比中一对对称位置元素之差,所以它只依赖反对称部分 。
旋度的一半是局部角速度。 记 ,叉乘矩阵为
逐项写出 :
与 逐个比较: 位置 , 位置 , 位置 ,完全吻合。因此
刚体以角速度 转动时,速度为 。比较可知,流体团的局部刚体旋转角速度为
三维反对称矩阵恰好有三个独立分量,所以能用一个轴向量表示,这也是旋度只在三维时才是向量的原因。二维情形下,反对称部分只有一个分量,旋度写成标量 ,表示垂直于平面的转动。
例子。 是绕 轴、角速度为 2 的刚体旋转。它的散度为 0,旋度为 ,正好是角速度的两倍。
雅可比的完整分解
把雅可比分成对称部分和反对称部分,再从对称部分中分出迹:
三部分的自由度分别为 ,合起来正好是雅可比的九个元素。散度只读出第一部分(1 个数),旋度只读出第三部分(3 个数),中间 5 个自由度的保体积形变两者都读不到。
反例。 的雅可比为 ,散度和旋度都为零,但它沿 方向拉伸、沿 方向压缩。所以仅凭一点的散度和旋度,无法恢复该点的雅可比。
简单剪切。 图中的 的雅可比为 。它的对称部分 是形变,反对称部分 是角速度为 (顺时针)的旋转,二维旋度 。看起来只是”一层层平移”的流动,其实同时包含旋转和形变。
积分形式:通量与环流
散度和旋度还有不依赖坐标的定义。
散度是单位体积的净流出通量。 取以 为中心、边长为 的小立方体。垂直于 轴的两个面上,净流出通量为
另外两对面同理。三项相加,得到净流出通量 。所以
旋度分量是单位面积的环流量。 对过 、法向为 的小曲面 ,
其中边界方向与 满足右手定则。
把这些局部关系在整个区域上累加,内部相邻小块的贡献互相抵消,只剩下边界,于是得到两个积分定理:
| 定理 | 公式 | 含义 |
|---|---|---|
| 散度定理 (Divergence Theorem) | 内部净源的总和等于边界的净通量 | |
| Stokes 定理 | 面内局部环流的总和等于边界环流 |
两个定理都要求场足够光滑、区域足够规则,并且定向一致。第 3 章推导泊松方程时会用到散度定理。
两个零恒等式
梯度场无旋:。以第一分量为例,它等于 。 二阶连续可微时混合偏导可交换,所以为零;另外两个分量同理。
旋度场无散:。展开得
六项两两抵消,原因同样是混合偏导可交换。
两个恒等式都要求场足够光滑,在奇点处不能直接套用。
反过来,零旋度能否推出存在势函数? 只有在单连通区域上才行。考虑定义在 上的场
直接计算可得 ,所以处处 。但沿单位圆一周,,环流量为 。如果 ,任意闭合回路上的环流都应为零,所以这个场没有全局势函数。原因是区域在原点处有一个洞。
使用约定
- 是微分算符,不是数值向量。“点乘""叉乘”的写法便于记忆,但作用于乘积时仍需使用乘积法则,例如 。
- 本章公式只适用于笛卡尔坐标。柱坐标和球坐标中,基向量随位置变化,散度和旋度的分量公式会多出尺度因子项,不能照搬。
- 有些中文材料把 称为”哈密顿算子”。为避免与量子力学中的哈密顿算符 混淆,本站称它为 nabla 或 del。
机器人中的例子
光流的散度与碰撞时间
无人机或机器人沿光轴方向接近一个正对着的平面时,图像上的光流 (Optical Flow) 从一点向外辐射。这一点称为扩展焦点 (Focus of Expansion, FOE)。光流的散度直接给出碰撞时间 (Time to Contact, TTC)。
设相机以速度 沿光轴前进,平面在深度 处。在相机坐标系中,场景点的速度为 。归一化图像坐标为 ,,对时间求导:
光流场 正是图中的”源”场,FOE 在主点。它的散度为
这个结论有三个值得注意的地方:
- 只要测出光流的散度,就能得到碰撞时间 ,不需要知道深度 和速度 各自是多少;
- 换成像素坐标 时,光流和坐标同时乘以 ,所以散度不变,也不需要知道焦距;
- 散度反映的是”膨胀”,而第 1 章投影雅可比的第三列告诉我们,沿光轴的运动正是使图像点沿径向移动的那部分。
正因为这些性质,光流散度被广泛用于微型飞行器的避障和降落控制,例如保持散度恒定以实现平滑减速着陆。上面的推导只适用于纯前向运动和正对的平面;存在横向平移、平面倾斜或旋转时,散度中还会出现其他项,需要先做补偿。
光流的旋度则对应绕光轴的转动。若相机以角速度 绕光轴旋转,场景点在相机坐标系中的速度为 ,图像光流为 ,其二维旋度为 :大小是角速度的两倍,方向与相机转动相反。这与上一节”旋度的一半是角速度”一致。
可视化 / 代码
第一段代码在网格上数值计算四个典型场的散度和旋度,对应图中的四个面板,并给出画图的代码。
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(-1, 1, 41)
y = np.linspace(-1, 1, 41)
X, Y = np.meshgrid(x, y) # 形状 (ny, nx),X 沿 axis=1 变化
fields = {
"source": (X, Y), # 源
"vortex": (-Y, X), # 涡旋
"shear": (Y, 0 * X), # 简单剪切
"strain": (X, -Y), # 纯形变
}
def div_curl(P, Q):
dP_dy, dP_dx = np.gradient(P, y, x)
dQ_dy, dQ_dx = np.gradient(Q, y, x)
return dP_dx + dQ_dy, dQ_dx - dP_dy
for name, (P, Q) in fields.items():
d, c = div_curl(P, Q)
print(f"{name:7s} div = {round(d.mean(), 2) + 0.0:+.2f} curl = {round(c.mean(), 2) + 0.0:+.2f}")
# source div = +2.00 curl = +0.00
# vortex div = +0.00 curl = +2.00
# shear div = +0.00 curl = -1.00
# strain div = +0.00 curl = +0.00
fig, axes = plt.subplots(1, 4, figsize=(12, 3.4))
s = slice(None, None, 4)
for ax, (name, (P, Q)) in zip(axes, fields.items()):
ax.quiver(X[s, s], Y[s, s], P[s, s], Q[s, s])
ax.set_title(name)
ax.set_aspect("equal")
plt.show()
第二段代码合成接近平面时的光流,验证散度等于 ,并验证绕光轴旋转时旋度为 。
import numpy as np
W, Z = 2.0, 10.0 # 接近速度 2 m/s,深度 10 m,tau = 5 s
f = 400.0 # 焦距(像素)
px = np.linspace(-200, 200, 81) # 以主点为原点的像素坐标
U, V = np.meshgrid(px, px)
xn, yn = U / f, V / f # 归一化坐标
# 纯前向运动:像素光流 = f * (W/Z) * (x, y)
flow_u, flow_v = f * xn * W / Z, f * yn * W / Z
div = np.gradient(flow_u, px, axis=1) + np.gradient(flow_v, px, axis=0)
print(div.mean(), 2 * W / Z, "tau =", 2 / div.mean()) # 约 0.4 0.4 tau = 5.0(浮点舍入)
# 绕光轴旋转:像素光流 = f * omega * (y, -x)
omega = 0.3
rot_u, rot_v = f * omega * yn, -f * omega * xn
curl = np.gradient(rot_v, px, axis=1) - np.gradient(rot_u, px, axis=0)
print(curl.mean(), -2 * omega) # 约 -0.6 -0.6
参考资料
- Koenderink, J. J., & van Doorn, A. J. (1976). Local structure of movement parallax of the plane. Journal of the Optical Society of America, 66(7), 717–723.
- Lee, D. N. (1976). A theory of visual control of braking based on information about time-to-collision. Perception, 5(4), 437–459.
- Herissé, B., Hamel, T., Mahony, R., & Russotto, F.-X. (2012). Landing a VTOL unmanned aerial vehicle on a moving platform using optical flow. IEEE Transactions on Robotics, 28(1), 77–89.
- Schey, H. M. (2005). Div, Grad, Curl, and All That (4th ed.). W. W. Norton.
- OpenStax. Calculus Volume 3, Section 6.5: Divergence and Curl. https://openstax.org/books/calculus-volume-3/pages/6-5-divergence-and-curl
- 3Blue1Brown. Divergence and curl: The language of Maxwell’s equations, fluid flow, and more. https://www.3blue1brown.com/lessons/divergence-and-curl/
- MIT OpenCourseWare. 18.02SC Multivariable Calculus (Fall 2010),第四单元:三维空间中的三重积分与曲面积分(含散度定理与 Stokes 定理)。https://ocw.mit.edu/courses/18-02sc-multivariable-calculus-fall-2010/