论文信息：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(1), 93–111. https://doi.org/10.1007/s00190-017-1057-3

计算起伏密度界面（例如地球波动的地表）产生的重力及其梯度，是地球物理正演与重力数据处理的基础问题。区域乃至全球尺度下，传统的局部平面地形改正不再适用，需要真正考虑地球曲率的球面或椭球地形（连同船测水深）影响。本文由一作张壹与中国地质大学（武汉）的陈超合作完成，提出了一种以多面体来表示全球密度界面的新方法（Polyhedral Global Interface Model，PGIM），并配套设计了**近区替换（Near Zone Replacement，NZR）**的自适应加速策略，以解析解一次性给出位、重力与重力梯度张量的正演结果。下面围绕这篇 Journal of Geodesy 论文梳理其动机、思路与主要结论。

研究背景：球面/椭球地形重力效应为什么难算

把起伏地形与海底的作用从重力观测中扣除，即所谓地形改正（terrain correction），是重力测量被用来研究地球外部形状与壳内密度结构的前提。对局部小范围，Bullard-B 改正足以处理地球曲率；但当研究跨越较大地理区域甚至全球、且涉及高山地区剧烈的高程变化时，就必须做球面或椭球地形（含水深）改正（Lafehr, 1998；An et al., 2015）。

已有正演方法大体分两类，各有其短板：

  • 球谐域方法：从 Parker 公式的笛卡尔频域推广到球谐域（Wieczorek & Phillips, 1998）或直接用于卫星重力及其梯度（Balmino et al., 2011；Hirt et al., 2012；Bouman et al., 2013）。但如 Eshagh (2013) 所指出的，球谐表达中含 $(r/R)^n$ 项，当观测点地心距 $r$ 大于平均半径 $R$ 时，高阶球谐会数值失稳。
  • 空域方法：把地表分成规则小块（球面棱柱、tesseroid 等）再累加（Heck & Seitz, 2007；Wild-Pfeiffer, 2008；Asgharzadeh et al., 2007；Grombein et al., 2013；Roussel et al., 2015）。这类方法精度随源与观测点距离变化、计算耗时，且观测点贴近源时可能奇异；更重要的是，它们大多建立在地理坐标（经度—纬度）网格上，而经纬度网格在从赤道向两极移动时会严重变形，导致极区分辨率退化（Sahr et al., 2003）。

这些缺陷构成了本文的出发点：能否构建一种全球分辨率近似均匀、且无需频繁处理奇异的密度界面模型，并同时、精确地给出重力与梯度？

多面体表示的全球密度界面模型
图1　本级 8 的 PGIM 地球表面模型（引自论文 Figure 4）

以二十面体为基、经多级递归细分后映射到地表得到的多面体界面模型。因采用 geodesic 剖分而非经纬度网格，它在全球各处保持近乎一致的分辨率与形状。上方放大图展示了模型的实际三角面片结构。

方法：用二十面体剖分构建全球多面体界面模型（PGIM）

本文的方法建立在测地离散全球网格（Geodesic DGGS）思想之上。与依赖地理坐标的网格不同，DGGS 用球面多边形（如球面三角形、六边形）作为地球表面的拓扑等价替代，其核心优势是在全球保持几乎恒定的分辨率（Sahr et al., 2003）。PGIM 的构建分五步：

  1. 选取规则基底多面体：本文选用**正二十面体（icosahedron）**作为基底。与其它柏拉图立体相比，二十面体面小、向球面映射时畸变最小，因而三角面片的大小和形状差异也被压到最小。
  2. 确定方位：最常用的是把一个顶点置于极点、一条边与本初子午线对齐（图1a）。
  3. 层次化剖分：采用 Class I（交替）递归细分，把每个二十面体面一分为四（连接各边中点）。细分 $n$ 次后，面片数为 $20 \times 4^n$；$n=0$ 即基底二十面体，$n$ 越大剖分越细。
  4. 确定顶点地心距：由于 PGIM 顶点坐标大多不与等角高程网格（如 ETOPO1）的节点重合，先用等角网格计算界面地心距，再用球面双线性插值（Wang et al., 2006）求各顶点处的地心距。
  5. 逆映射构造 3-D 模型：把剖分顶点按地心距映射回界面表面，即得球面或椭球密度界面模型（图1）。

不同界面可方便地通过更换数据（ETOPO1 地形水深、WGS84 参考椭球、EGM96 大地水准面）定义。PGIM 的分辨率以等效角衡量：从第 6 级的 1.07° 到第 10 级的 0.07°，每升一级分辨率近似翻倍；第 8 级即有 1,310,720 个三角面片，平均面积约 389 km²、平均边长约 30 km。全体面片面积、边长分布集中，最小内角大于 54°，近似等边，故全球分辨率一致。

正演方面，PGIM 外部重力场（位、梯度、梯度张量）采用 Werner & Scheeres (1996) 的闭式解析解，把多面体场表示为所有面与棱贡献之和，理论上正演过程不引入数值误差、也无须处理奇异点。此正演算法对位、梯度、梯度张量使用同一套核矩阵，因而可"一劳永逸"地同时给出多个分量，且算法全球统一、无需拼接，适用于陆地、海洋乃至海底、航空与卫星数据。

方法：近区替换（NZR）如何省下四倍计算量

球面/椭球地形重力效应有个显著特征：观测点附近的起伏对重力与梯度影响大，而远处起伏的影响迅速衰减（Mikuška et al., 2006；Lafehr, 1991）。据此可适当降低远处界面的分辨率而不引入明显误差，同时大幅削减计算量。NZR 的实现是：对观测点 $P$，先用整体（低分辨率）模型算出重力效应 $g_{glo}$（图6a）；再按阈值角 $\theta_P$ 圈定"近区"，把近区替换成更细（$n+1$ 或 $n+2$ 级）的局部细分模型（图6b）；于是 $P$ 点的总效应为

$$g_{NZR} = g_{glo} + g_{imp} - g_{ori},$$

即"整体值 + 细化局部改进值 − 局部原值"（图6c）。细化的局部模型边缘会留下缝隙（缺口），但只要阈值角足够大（远处地形的贡献本就可忽略），缺口带来的误差即可忽略。

误差分析表明，PGIM 与参考球面的差异在二十面体原顶点附近更小、在面心附近相对较大；级数越高差异越小——第 10 级的 PGIM 相对理论固体球解（979,828.494 mGal）仅约 0.5 mGal 误差，足以满足全球研究精度。

应用：青藏高原—印度洋区域的椭球 Bouguer 壳层改正

作者把 PGIM 用于 60°E–120°E、0°N–60°N（同时覆盖青藏高原与部分印度洋）区域、基于 WGS84 参考椭球的 Bouguer 壳层改正。该区域地形起伏近乎全球最剧烈，是检验方法的好场所。为做改正，用地形水深面、大地水准面等构造了三个 PGIM（M1：陆地 2.67 g/cm³；M2：陆地—海洋差分层 1.64 g/cm³；M3：海洋 1.03 g/cm³），壳层改正取 $V_{tc}=v_{M1}-v_{M2}-v_{M3}$。

沿 85°E、15°N–30°N 的典型剖面先选定参数（NZR 阈值 20°、全球/区域模型取第 8 级、局部模型第 10 级）。与不用 NZR 相比，带 NZR 的运行时间从约 3046 s 降到 752 s，相差约 4 倍，而由此带来的重力/梯度误差几乎不可察觉。与 Du et al. (2012) 的空间方法对比，两者差异在印度洋与印度半岛地区为 −4 至 1 mGal、在青藏高原增大到 −11 至 5 mGal，垂直梯度差异在低海拔区约 ±1E、在青藏高原为 −25E 至 18E——总体差异微小，精度足以满足区域/全球研究。

重力扰动与布格异常对比
图2　应用区域的重力改正结果（引自论文 Figure 13）

(a) 由 EIGEN-6C4 全球重力模型算得的重力扰动；(b) 用 (a) 扣除地形重力效应（图12a）得到的 Bouguer 重力异常；(c) ICGEM 在线服务（Bouguer 平板法）得到的 Bouguer 异常。本文方法在海洋/陆地之间幅度更均衡，范围约 −496 至 +502 mGal，并因扣除了地球曲率而比平板法更贴近真实的重力效应。

同一参数下的区域应用结果（图2b）表明，用 EIGEN-6C4 重力扰动扣除 PGIM 地形效应后，所得的 Bouguer 异常与 ICGEM 在线服务（采用 Bouguer 平板法、忽略地球曲率）相比有两点不同：ICGEM 结果在陆地上显示极强的负异常（喜马拉雅一带最大约 −724.58 mGal）、海洋上则是中等正异常（最大约 +365 mGal）；而本文方法给出的幅度在海洋正异常与陆地负异常之间更为均衡，总体范围约 −496.47 至 +502.49 mGal，因为它在扣除时正确计入了椭球曲率与曲面地形效应，从而更真实地反映了地球内部密度对比对重力的影响。

结论与展望

文章的核心贡献可概括为四点：(1) 界面用三角面片逼近，分辨率在全球近似恒定，且细分过程层次化、易于获得多分辨率模型；(2) 正演从位到梯度张量全部采用无奇异的精确解析解，且各分量的核矩阵统一、计算开销低，算法全球统一无需拼接，可应用于陆地、海洋、海底、航空及卫星数据等多种场景；(3) PGIM 的总误差主要来自原始高程数据，而用平面三角片近似球面曲率与界面起伏引入的误差可忽略，对需要超高精度的场合建议不使用近区替换；(4) 除地形效应外，该方法也适合重建任意近似等深（isometric）体的密度界面并正演其外部场，可推广到均衡改正等应用。

这项工作的意义在于，为区域/全球尺度的球面或椭球地形重力效应提供了一个分辨率均匀、精度高且计算高效的新工具——它在概念上把地球表面的处理和局部平板近似提升到了全球统一的多面体框架。作者也指出，方法源码可按索取提供，方便后续研究复现与拓展。

延伸阅读与引用

  • 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.
  • Sahr, K., White, D., & Kimerling, A. J. (2003). Geodesic discrete global grid systems. Cartography and Geographic Information Science, 30(2), 121–134.
  • Grombein, T., Seitz, K., & Heck, B. (2013). Optimized formulas for the gravitational field of a tesseroid. Journal of Geodesy, 87(7), 645–660.
  • Tsoulis, D. (2012). Analytical computation of the full gravity tensor of a homogeneous arbitrarily shaped polyhedral source using line integrals. Geophysics, 77(2), F1–F14.
  • Du, J., Chen, C., Liang, Q., Wang, L., Zhang, Y., & Wang, Q. (2012). Gravity anomaly calculation based on volume integral in spherical cap and comparison with the Tesseroid–Taylor series expansion approach. Acta Geodaetica et Cartographica Sinica, 41(3), 339–346 (in Chinese).

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