论文信息：Zhang, Y., Mooney, W. D., & Chen, C. (2018). Forward calculation of gravitational fields with variable resolution 3D density models using spherical triangular tessellation: Theory and Applications. Geophysical Journal International, 215(1), 702–712. https://doi.org/10.1093/gji/ggy278

在全球或区域尺度上，重力正演是约束地球及其他行星密度结构的基本手段。然而，如何在球坐标系下既精确又高效地构造三维密度模型并进行正演，长期存在一个"分辨率"难题：传统的经纬度网格在接近两极时会严重畸变，无法保证全球均匀的分辨率。本文以第一作者张壹（本文作者之一）为代表的工作，提出了一套基于球面三角剖分（Spherical Triangular Tessellation，STT）的完整方案——从构造可变分辨率的三维密度模型，到实现从地表到卫星高度几乎恒定精度的重力场正演。下面围绕这篇 GJI 论文梳理其动机、方法要点与两个应用实例。

为什么要做"变分辨率"的球面正演

在地球物理中，重力正演的现有方法可分为空间域与球谐域两类。空间域常用做法是把球面划分为小块（如球面棱柱 tesseroid），通过求各块效应的叠加得到重力场；由于 tesseroid 外场没有解析解，通常要借助泰勒展开或 Gauss–Legendre 数值积分（Heck & Seitz, 2007；Asgharzadeh et al., 2007 等）。另一个普遍问题是这些方法大多以**地理坐标（经纬度）**为底层网格，而经纬度网格单元的实际面积与形状随纬度剧烈变化，越靠近两极畸变越严重。

为此，本文采用了测地离散全球网格系统（DGGS）的思想，用球面三角形作为覆盖地球表面的基本单元——即球面三角剖分。STT 的首要优点是能在全球提供几乎恒定的分辨率，且由于只有三角形面片，它在几何上非常适合刻画起伏的界面；对相同的空间分辨率，STT 的顶点数相比球坐标下的矩形网格约少 25%，因而更节省计算量，尤其适合大规模正演与反演。

方法：从二十面体出发的球面三角剖分建模

STT 的构建从一个正二十面体开始：把它的每个三角形面不断细分（连接各边中点），再把这些细分面片映射到其外接球面上，就得到层层加密的球面三角网格（图1）。一个 n 级 STT 由 20 × 4ⁿ 个三角形面片组成，其分辨率约为 60°/2ⁿ，每升一级、边长减半、分辨率提高一倍。

球面三角剖分构建过程示意
图1　球面三角剖分（STT）的构建过程（引自论文 Figure 1）

(a) 正二十面体示意；(b) 二十面体各面的细分；(c) 把细分面片映射到外接球面得到的 STT。图片来源张壹与陈超（2018）。

为了使分辨率能够按需连续变化，本文引入了三类约束条件：点约束（在指定位置细化）、弧约束（沿大圆弧轨迹细化，复杂路径可拆成多条弧）、以及球面多边形约束（在指定区域内细化，同时可用来裁剪出区域模型）。对每条约束设定终止条件（最大细分级数或最大分辨率）来控制加密过程；距离目标区越远，细分级数越低，从而实现"近处加密、远处稀疏"的自适应网格。

在 STT 之上，先用球面双线性插值把网格顶点映射到起伏界面上，构造二维界面模型，再以该界面为模板、向不同深度"挤出"形成多层三维密度模型。相邻层连接侧面棱后，模型的基本单元是球面三角棱柱（图2）。由于各层界面起伏不同、层与层可能相交，棱柱会退化为金字塔、四面体等四种基本形态。所幸重力场是势场，一个三角棱柱的外场总可等价地表示为两个"一个顶点位于原点"的四面体的外场之差，因此各种退化形态都无需特殊处理，统一用 Werner & Scheeres（1996）的闭式解析公式计算即可。

三维密度模型几何单元示意
图2　三维密度模型的基本几何单元（引自论文 Figure 6）

蓝色与红色三角网格分别表示模型的顶面与底面。当顶、底面分离时为截断三角棱柱（Type 1）；共顶点或共边时退化为金字塔或四面体（Type 2、3）；顶底面相交时又成另一种形态（Type 4）。

为了检验精度，作者以地球参考球（半径 6371.0088 km、密度 1.0 g/cm³）为对象，用 6–9 级 STT 构造模型，并在 1 km、10 km、25 km、50 km、100 km、150 km 与 250 km 等多个观测高度上正演重力与径向重力梯度，与理论值对比（表1、表2）。结果表明：正演误差随模型分辨率快速下降——1 km 高度处 6 级模型的差异约 23.7 mGal，9 级模型则降至约 0.374 mGal（标准差 0.04 mGal）。更重要的是，随观测高度上升误差几乎不增加：9 级模型在 1 km 与 250 km 高度之差仅约 7.8% 的相对变化，说明该方法能从地形表面一直到卫星高度都保持近乎恒定的计算精度。

应用一：西南美国沿海的地形重力效应

第一个算例是正演地形重力效应。作者用 ETOPO1 模型的地形与水深构建密度模型，并采用多分辨率 STT来提升计算效率：对观测点而言，远处地形的分辨率与精度没有近处重要，因此细分级数由目标区的约 0.0625°（10 级）向外递减到 8 级。正演结果与 Uieda 等（2016）采用 0.01°×0.01° tesseroid 的方法对比，最大正演重力约 189.43 mGal、最小约 −262.785 mGal，两种方法在地形起伏剧烈处的差异约为 ±1 mGal。作者认为其余差异主要来自 tesseroid 用断续球面棱柱近似地表面的误差，而本方法因模型分辨率足够高、正演误差很小，在地形陡变区具有更好的精度。

应用二：北美大陆的残差地幔重力异常

更综合的算例是正演北美大陆的残差地幔重力异常。作者将 NACr14 模型（Tesauro et al., 2014）与 CRUST1.0 模型（Laske et al., 2013）结合，构建了一个包含水、冰、三层沉积、上中下地壳与莫霍面的九层一体化地壳密度模型；NACr14 提供北美的大陆地壳分层与界面深度（其密度由 P 波速度—密度经验关系得到），其余区域采用 CRUST1.0。三维模型采用多分辨率 STT：在北美大陆边界用 7 级细分（约 0.5°分辨率），全球其余区域用 6 级，以提高计算效率与精度。

正演针对三层一维参考模型（0–15 km 密度 2.7 g/cm³、至 40 km 为 2.94 g/cm³、至 75 km 为 3.35 g/cm³）进行改正，再结合 EIGEN-6c4 重力异常扣除各层效应，得到残差地幔重力异常（图3）。

北美大陆残差地幔重力异常
图3　应用算例：残差地幔重力异常（引自论文 Figure 15）

北美大陆的残差地幔重力异常 (a) 引自 Kaban 等（2014）与 (b) 本文用 EIGEN-6c4 重力扰动扣除三维密度模型正演改正后得到的结果。残差异常大多在 ±400 mGal 以内，反映上地幔密度变化。

结果显示，残差异常大多在 ±400 mGal 以内，最显著的特征是稳定古老的中东部与相对活动的西部新生代—中生代省之间的差异，与 Mooney & Kaban（2010）及 Kaban 等（2014）的先前结果高度一致。最大的负异常（< −300 mGal）位于美国西南部，对应低 Pn 波速与低上地幔速度，指示热成因；而在元古宙与太古宙地体上方异常为正且较强（> 100 mGal），与上地幔正地震波速异常一致。

结语与展望

这项工作提出了一套以球面三角剖分为基础的完整空间域方法：既能构建具有可变分辨率、起伏界面的三维密度模型，又能实现从地面到卫星高度的恒定精度正演。在模型构造上，它结合了 Ballard 等（2009）的变分辨率 STT 细分与张壹和陈超（2018）的界面建模；在正演上，用两个四面体的差来统一处理各种退化棱柱单元，兼得闭式解析的精度与自适应网格的高效。作者指出，该方法适用于对精度和建模灵活性要求较高的区域与全球研究，适合作为后续重力反演的正演引擎。

作为后续延伸，读者可以进一步关注这一团队基于同样思路在球面网格上开展的地球物理反演工作，以及 DGGS 类网格在重力、磁力等多源数据联合建模中的应用趋势。

延伸阅读与引用

  • Zhang, Y., & Chen, C. (2018). Forward calculation of gravity and its gradient using polyhedral representation of density interfaces: An application of spherical or ellipsoidal topographic gravity effect. Journal of Geodesy, 92(2), 205–218.
  • Ballard, S., Hipp, J. R., & Young, C. J. (2009). Efficient and accurate calculation of ray theory seismic travel time through variable resolution 3D earth models. Seismological Research Letters, 80(6), 989–999.
  • Werner, R. A., & Scheeres, D. J. (1996). Exterior gravitation of a polyhedron derived and compared with harmonic and mascon gravitation representations of asteroid 4769 Castalia. Celestial Mechanics and Dynamical Astronomy, 65(3), 313–344.
  • Kaban, M. K., Tesauro, M., Mooney, W. D., & Cloetingh, S. A. (2014). Density, temperature, and composition of the North American lithosphere—New insights from a joint analysis of seismic, gravity, and mineral physics data: 1. Density structure of the crust and upper mantle. Geochemistry, Geophysics, Geosystems, 15(12), 4781–4807.
  • Tesauro, M., Kaban, M. K., Mooney, W. D., & Cloetingh, S. (2014). NACr14: A 3D model for the crustal structure of the North American Continent. Tectonophysics, 631, 65–86.

本文为对该论文的中文解读，图片均引自论文原文。