Song's Blog

外观

球面单元映射的 Jacobian 连续性分析

约 2629 字大约 9 分钟

球面网格

2026-04-28

做球面网格时，一个很自然的问题是：单元是怎么贴到球面上的？这个问题看起来只是几何构造，但一旦进入有限元，尤其是要用到 Piola 变换的 Hdiv 空间，映射的微分结构就会直接影响离散误差。

这篇文章想比较两种常见做法：

  • 先在线性三角形（四边形）里插值，再做径向投影；
  • 先给球面做全局参数化，再把参考单元映到参数域里。

我们真正想问的不是“点能不能落到球面上”，而是更细一点的问题：相邻两个单元在公共边上，映射的 Jacobian 到底是不是连续的。

我们考虑一个从参考单元到球面的映射

F:(ξ,η)↦x(ξ,η)∈R3,∣x∣=RF : (\xi,\eta) \mapsto \mathbf x(\xi,\eta)\in\mathbb R^3, \qquad |\mathbf x|=R

对应的微分写成

DF=∂x∂(ξ,η)DF=\frac{\partial \mathbf x}{\partial(\xi,\eta)}

当两个单元共享一条边时，我们关心三件事：

  • 边上位置是否连续；
  • 沿公共边的切向导数是否连续；
  • 横跨公共边的导数（法向导数）是否连续。

前两件事解决的是“能不能拼起来”，第三件事解决的是“拼起来以后是不是光滑”。如果横跨公共边的导数不连续，那么这个映射通常只有 C0C^0 连续性，而不是 C1C^1 连续性。

线性插值加径向投影

先看最常见的一种构造。给定三角形三个顶点

Xi=(xi,yi,zi),i=1,2,3\mathbf X_i = (x_i,y_i,z_i), \quad i=1,2,3

先在平面三角形上做线性插值：

X(ξ,η)=X1(1−ξ−η)+X2ξ+X3η\mathbf X(\xi,\eta) = \mathbf X_1(1-\xi-\eta)+\mathbf X_2\xi+\mathbf X_3\eta

再做径向投影到球面：

x=RX∣X∣\mathbf x = R \frac{\mathbf X}{|\mathbf X|}

这个映射的 Jacobian 是

J=R∣X∣(I−XXT∣X∣2)(X2−X1,X3−X1)J = \frac{R}{|\mathbf X|} \left( I - \frac{\mathbf X\mathbf X^T}{|\mathbf X|^2} \right) \begin{pmatrix} \mathbf X_2-\mathbf X_1, & \mathbf X_3-\mathbf X_1 \end{pmatrix}

公式本身并不复杂，但连续性问题要在单元交界处看。设两个相邻三角形共享边 X1–X2X_1\text{--}X_2：

  • 左单元：(X1,X2,X3)(X_1,X_2,X_3)
  • 右单元：(X1,X4,X2)(X_1,X_4,X_2)

相邻三角形单元的邻接关系

它们各自的线性映射部分写成

XA(ξ,η)=X1(1−ξ−η)+X2ξ+X3η\mathbf X_A(\xi,\eta) = \mathbf X_1(1-\xi-\eta)+\mathbf X_2\xi+\mathbf X_3\eta

XB(ξ,η)=X1(1−ξ−η)+X4ξ+X2η\mathbf X_B(\xi,\eta) = \mathbf X_1(1-\xi-\eta)+\mathbf X_4\xi+\mathbf X_2\eta

这里有一个小细节：同一条公共边，在左单元里对应 η=0\eta=0，在右单元里对应 ξ=0\xi=0。

边上的位置与切向导数

先看最容易的一步。在公共边上，左右两侧其实都退化成同一条线段参数化

Xe(t)=X1(1−t)+X2t\mathbf X_e(t)=\mathbf X_1(1-t)+\mathbf X_2 t

所以投影以后，边上的位置完全一样：

xA(t,0)=xB(0,t)=RXe(t)∣Xe(t)∣\mathbf x_A(t,0)=\mathbf x_B(0,t) = R\frac{\mathbf X_e(t)}{|\mathbf X_e(t)|}

再看沿边方向的切向导数。由于这条边在左单元中由 ξ\xi 参数化，在右单元中由 η\eta 参数化，所以要比较的是

∂ξxA和∂ηxB\partial_\xi \mathbf x_A \quad \text{和} \quad \partial_\eta \mathbf x_B

而线性部分满足

∂ξXA=X2−X1,∂ηXB=X2−X1\partial_\xi \mathbf X_A = \mathbf X_2-\mathbf X_1, \qquad \partial_\eta \mathbf X_B = \mathbf X_2-\mathbf X_1

于是

∂ξxA=DF(Xe)(X2−X1)=∂ηxB\partial_\xi \mathbf x_A = DF(\mathbf X_e)(\mathbf X_2-\mathbf X_1) = \partial_\eta \mathbf x_B

也就是说，这种构造至少能保证两件事：

  • 边上的点能接上；
  • 沿边方向的切向导数也能接上。

为什么法向导数一般不连续

真正的问题出在“跨过这条边往单元内部走”时的导数。为此，把径向投影单独记成

x=F(X)=RX∣X∣\mathbf x = F(\mathbf X) = R\frac{\mathbf X}{|\mathbf X|}

它的微分对任意方向向量 V\mathbf V 都有

DF(X) V=R∣X∣(I−XXT∣X∣2)VDF(\mathbf X)\,\mathbf V = \frac{R}{|\mathbf X|} \left( I-\frac{\mathbf X\mathbf X^T}{|\mathbf X|^2} \right)\mathbf V

现在看横跨公共边的方向。在左单元里，这个方向对应 ∂η\partial_\eta；在右单元里，对应 ∂ξ\partial_\xi。而

∂ηXA=X3−X1,∂ξXB=X4−X1\partial_\eta \mathbf X_A = \mathbf X_3-\mathbf X_1, \qquad \partial_\xi \mathbf X_B = \mathbf X_4-\mathbf X_1

所以边上两侧的法向导数分别是

∂ηxA=DF(Xe)(X3−X1)\partial_\eta \mathbf x_A = DF(\mathbf X_e)(\mathbf X_3-\mathbf X_1)

∂ξxB=DF(Xe)(X4−X1)\partial_\xi \mathbf x_B = DF(\mathbf X_e)(\mathbf X_4-\mathbf X_1)

它们一减，得到

∂ηxA−∂ξxB=DF(Xe)(X3−X4)\partial_\eta \mathbf x_A-\partial_\xi \mathbf x_B = DF(\mathbf X_e)(\mathbf X_3-\mathbf X_4)

进一步展开得

∂ηxA−∂ξxB=R∣Xe∣(I−XeXeT∣Xe∣2)(X3−X4)\partial_\eta \mathbf x_A-\partial_\xi \mathbf x_B = \frac{R}{|\mathbf X_e|} \left( I-\frac{\mathbf X_e\mathbf X_e^T}{|\mathbf X_e|^2} \right) (\mathbf X_3-\mathbf X_4)

这个式子很说明问题。若要法向导数连续，必须有

(I−XeXeT∣Xe∣2)(X3−X4)=0\left( I-\frac{\mathbf X_e\mathbf X_e^T}{|\mathbf X_e|^2} \right) (\mathbf X_3-\mathbf X_4)=0

也就是说，X3−X4\mathbf X_3-\mathbf X_4 在球面切平面里的分量必须正好消失。换句话说，它必须平行于径向方向。这个条件在一般网格里几乎不会碰巧成立，所以通常都会得到

∂ηxA≠∂ξxB\partial_\eta \mathbf x_A \neq \partial_\xi \mathbf x_B

这就是这类映射最关键的特征：位置是连续的，沿边的切向导数也是连续的，但横跨边界的导数通常会跳。

所以它更像是“分片拼起来的球面”，而不是一个在边附近真正光滑的统一曲面。用一句话概括就是

映射本身在单元边界是连续的，但 Jacobian 一般不连续；因此整体通常只有 C0 连续性\boxed{ \text{映射本身在单元边界是连续的，但 Jacobian 一般不连续；因此整体通常只有 } C^0 \text{ 连续性} }

经纬参数化映射

再看另一种思路。不是每个单元自己拼几何，而是先给整个球面建立一套统一参数，比如经纬度 (lon,lat)(\mathrm{lon},\mathrm{lat})，然后把参考单元映到参数域里的小区域。定义

lat=lat1+η(lat2−lat1)lon=lon1+ξ(lon2−lon1)\begin{aligned} \mathrm{lat} &= \mathrm{lat}_1 + \eta(\mathrm{lat}_2-\mathrm{lat}_1) \\ \mathrm{lon} &= \mathrm{lon}_1 + \xi(\mathrm{lon}_2-\mathrm{lon}_1) \end{aligned}

然后令

x=R(cos⁡(lat)cos⁡(lon)cos⁡(lat)sin⁡(lon)sin⁡(lat))\mathbf x = R \begin{pmatrix} \cos(\mathrm{lat})\cos(\mathrm{lon})\\ \cos(\mathrm{lat})\sin(\mathrm{lon})\\ \sin(\mathrm{lat}) \end{pmatrix}

对应 Jacobian 写成

J=(∂ξx,∂ηx)J = \begin{pmatrix} \partial_\xi \mathbf x, & \partial_\eta \mathbf x \end{pmatrix}

其中

∂ξx=RΔlon(−cos⁡(lat)sin⁡(lon)cos⁡(lat)cos⁡(lon)0)\partial_\xi \mathbf x = R\Delta\mathrm{lon} \begin{pmatrix} -\cos(\mathrm{lat})\sin(\mathrm{lon})\\ \cos(\mathrm{lat})\cos(\mathrm{lon})\\ 0 \end{pmatrix}

∂ηx=RΔlat(−sin⁡(lat)cos⁡(lon)−sin⁡(lat)sin⁡(lon)cos⁡(lat))\partial_\eta \mathbf x = R\Delta\mathrm{lat} \begin{pmatrix} -\sin(\mathrm{lat})\cos(\mathrm{lon})\\ -\sin(\mathrm{lat})\sin(\mathrm{lon})\\ \cos(\mathrm{lat}) \end{pmatrix}

现在考虑两个相邻经纬单元：

  • 单元 1：[lon1,lon2][\mathrm{lon}_1,\mathrm{lon}_2]
  • 单元 2：[lon2,lon3][\mathrm{lon}_2,\mathrm{lon}_3]

它们的公共边满足

lon=lon2\mathrm{lon} = \mathrm{lon}_2

和前一种做法最不一样的地方在于：这里两个单元本来就是同一个全局映射

x(lon,lat)=R(cos⁡(lat)cos⁡(lon)cos⁡(lat)sin⁡(lon)sin⁡(lat))\mathbf x(\mathrm{lon},\mathrm{lat}) = R \begin{pmatrix} \cos(\mathrm{lat})\cos(\mathrm{lon})\\ \cos(\mathrm{lat})\sin(\mathrm{lon})\\ \sin(\mathrm{lat}) \end{pmatrix}

在不同区域上的局部限制。所以只要回到几何参数 (lon,lat)(\mathrm{lon},\mathrm{lat}) 来看，

∂x∂lon,∂x∂lat\frac{\partial \mathbf x}{\partial \mathrm{lon}}, \quad \frac{\partial \mathbf x}{\partial \mathrm{lat}}

它们天然就是同一个全局微分结构的一部分，因此在公共边上不会无缘无故跳掉。

当然，如果只盯着各自参考单元里的局部 Jacobian，也会看到一个细节：

∂ξx∝Δlon\partial_\xi \mathbf x \propto \Delta\mathrm{lon}

这里记

Δlon1=lon2−lon1,Δlon2=lon3−lon2\Delta \mathrm{lon}_1=\mathrm{lon}_2-\mathrm{lon}_1, \quad \Delta \mathrm{lon}_2=\mathrm{lon}_3-\mathrm{lon}_2

所以还得分两种情况来看。若网格均匀，也就是相邻两个经度区间等宽，

lon3−lon2=lon2−lon1⟺Δlon2=Δlon1\mathrm{lon}_3-\mathrm{lon}_2=\mathrm{lon}_2-\mathrm{lon}_1 \quad\Longleftrightarrow\quad \Delta\mathrm{lon}_2 = \Delta\mathrm{lon}_1

那么两侧以参考坐标写出的 Jacobian 确实完全一样：

J1=J2J_1 = J_2

但若网格非均匀，

Δlon1≠Δlon2\Delta\mathrm{lon}_1 \neq \Delta\mathrm{lon}_2

那么局部 Jacobian 会带上不同的尺度因子。这里要小心：这不表示几何映射本身不连续，只表示你拿来描述它的参考坐标缩放不一样。

所以这一类映射的本质不是“每个单元都恰好对上了”，而是“它们本来就是同一个全局参数化的局部片段”。因此可以概括为

几何映射由全局参数化给出，因此 Jacobian 在几何意义下连续\boxed{ \text{几何映射由全局参数化给出，因此 Jacobian 在几何意义下连续} }

两种构造的本质区别

线性插值加投影经纬参数化
是否来自全局参数化否是
单元映射来源每个单元独立构造全局映射的局部限制
边界位置连续连续
切向导数连续连续
法向导数一般不连续几何上连续
Jacobian边界跳跃几何上一致

如果把差别压缩成一句话，那就是：

  • 第一种方法是“每个单元自己造，再想办法投到球面上”；
  • 第二种方法是“先有整个球面，再把单元切出来”。

前者的问题不在位置，而在微分结构不是全局统一的；后者的优势也不只是公式好看，而是它确实来自一个统一的几何对象。

对有限元的影响

对 RT 元这类 H(div)H(\mathrm{div}) 有限元，Piola 变换写作

u=1JDFu^\mathbf u = \frac{1}{J} DF \hat{\mathbf u}

这时几何映射就不只是背景了，因为它的微分会直接进入离散空间的定义。

若采用第一种映射，则

  • DFDF 在边界上不连续
  • 法向导数出现跳变
  • Jacobian 在相邻单元之间发生跃迁

这会使通量自由度匹配更敏感，也更容易放大插值误差。

若采用第二种映射，则

  • FF 来自统一的全局参数化
  • DFDF 在几何上连续变化
  • Piola 变换更稳定

所以在相同离散阶数下，第二种做法通常会带来更好的几何一致性，也更不容易把几何误差放大到插值误差里。

总结

回头看，这两类构造的区别其实很清楚。

第一种映射里，相邻单元在公共边上能把点接起来，切向方向也能对上，但法向导数通常会跳，所以整体更像一个 C0C^0 的分片球面。

第二种映射里，各个单元本来就来自同一个全局参数化，因此共享同一套几何微分结构，边界上的导数行为也更自然。

第一种是“拼出来的球面”，第二种是“参数化出来的球面”。

实际计算建议

  • 对经纬网格，优先使用全局参数化映射。
  • 对非结构球面网格，线性插值加投影仍然可用，但要预期边界微分结构不连续带来的误差放大。
  • 对 RT 元或更一般的 H(div)H(\mathrm{div}) 方法，Jacobian 与导数连续性尤其重要。

版权归属：Guisong Wu

许可证：署名 4.0 国际 (CC-BY-4.0)