论文信息：Zhang, Y., Mooney, W. D., Chen, C., & Du, J. (2019). Interface inversion of gravitational data using spherical triangular tessellation: an application for the estimation of the Moon’s crustal thickness. Geophysical Journal International, 217(1), 703–713. https://doi.org/10.1093/gji/ggz026

从重力数据里反演密度界面的起伏（例如盆地的基底界面、Moho 起伏、月球地壳厚度），是重力勘探与行星地球物理中最经典的一类问题。规模一大，就必须考虑球面曲率，问题也就自然落到球坐标系下的正演与反演上来。本文以一作张壹为代表，提出用球面三角剖分（Spherical Triangular Tessellation, STT）作为 3D 密度模型的基本表示，在球坐标系下完成密度界面的正演与反演，并将其运用于月球地壳厚度的估计。下面围绕这篇 GJI 论文，梳理研究动机、方法要点与主要结果。

为什么需要"球面三角剖分"这类表示

在大区域尺度上用重力数据做界面反演，现有方法大致分两类：一类在球谐（SH）域开展工作（如 Wieczorek & Phillips, 1998 将 Parker 公式推广到球谐域并用于月球地壳厚度）；另一类在空间域把球面划分成矩形区域，用**球面棱柱单元（tesseroid）**逼近界面起伏（如 Heck & Seitz, 2007；Grombein 等, 2013；Uieda & Barbosa, 2017，后者用于南美洲 Moho）。这两类做法都很有价值，但在一些场景下存在不便：

  • 不规则的数据覆盖与高纬度区域：球谐方法在球面上是全局、等分辨率的，难以针对稀疏或非规则观测做局部加密；
  • 变分辨率 / 多分辨率模型：tesseroid 的矩形剖分在球面上做局部分辨率调整并不自然，相邻网格大小与形状难以平滑过渡；
  • 全球等分辨率模型：矩形剖分在极区会出现网格严重变形、单元形状不均匀的问题。

相比之下，STT——即对球面做三角剖分（可由 Delaunay 三角剖分、递归细分，或直接连接矩形格点对角线得到）——用连续三角面片逼近曲面，单元之间没有缝隙，且能自然地实现局部分辨率的加密。这正是本文选择 STT 作为模型构建基础的出发点。

方法：用三角棱柱搭起 3D 密度模型

作者沿球向"拉伸"STT 结构，构造出由上、下两个边界曲面夹着一个质量层的 3D 模型：把两个边界曲面上对应三角面片的侧边相连，就得到一个个三角棱柱体单元（triangular prism），棱柱的侧棱正好沿径向。模型的密度属性已知、参考界面固定后，重力场就完全由目标界面各个顶点的地心半径决定；反演目标界面，即是通过拟合重力观测来求解这些顶点的半径。

正演上，一个三角棱柱的引力场可表示为其两个"共原点四面体"引力场之差（g_prism = g(A'B'C'O) − g(ABCO)），而四面体外部引力场可采用 Werner & Scheeres (1996) 的多面体解析公式严格计算。相较于球谐级数在限制球内截断带来的距离依赖误差，解析方法在不同观测距离上具有一致的精度——这已在 Zhang 等 (2018) 中验证。

反演上，顶点半径与模型重力场之间是非线性关系。作者通过非线性共轭梯度法最小化目标函数。目标函数以"模型数据与观测差在观测标准差 σi 之内"为目标（使归一化残差趋近 1）。为此需要计算目标函数对各顶点半径的偏导数，而某顶点的偏导可分解为与其相邻的所有三角棱柱单元偏导之和；每个四面体内的偏导公式由本文附录给出。

合成算例：区域与全球两关都过了

作者先用两个合成算例检验方法。区域算例中，模型范围 32°–58°E、22°–48°N，地形与密度分布给定（图4），用约 1.0° 分辨率、1493 个三角面片、802 个顶点的 STT 构建模型，并在参考球上方 100 m 处正演重力与径向重力梯度、加入高斯噪声后反演。结果显示反演地形与真实地形高度吻合（图6），重力数据与径向梯度数据给出的平均误差分别约为 0.017 m（标准差 0.89 m）与 0.006 m（标准差 0.59 m），最大差异约 ±4 m（图7）——反演地形略偏光滑，这正是最小二乘意义下目标函数对非光滑起伏的天然抑制。

全球算例中，用一个球内接二十面体作为目标界面（其顶点半径为 1e+2 km，参考界面位于其下方 25 km，密度 1000 kg/m³），采用约 4° 分辨率、5120 个面片、2562 个顶点的 STT，在 10 km 高度正演并反演。反演结果在很大程度上恢复了二十面体的形状，重力与径向梯度数据对应界面与理论值的平均差异分别约为 2.89 m 与 −8.51 m。两个算例共同说明：该方法能够在球坐标下反演出常密度或变密度界面，反演的精细程度主要由数据覆盖与观测精度决定。

主要结果：月球地壳厚度

在月球应用中，作者采用月球绝对重力场模型 GRGM660PRIM 与地形模型 LRO LTM05 (LOLA)。构造一个四层密度模型（月海玄武岩、地壳、地幔、核），先扣除月核与月海盆地的重力效应，再将剩余重力异常归因于一个以地幔—地壳密度差为横向密度变化的单一质量层，其上下界面分别为月球的 Moho 与核幔边界（CMB）——这样参考界面取 CMB，目标界面即 Moho 起伏。模型参数见表1：月核半径 370 km、密度 6600 kg/m³，地幔密度取 3356 kg/m³，地壳平均密度 2794 kg/m³（取 7% 孔隙度），月海玄武岩密度 3159 kg/m³。

反演采用约 1.0° 分辨率、81,920 个三角面片、40,962 个顶点的 STT，参考球半径 1738 km，重力数据在参考球上方 10 km 计算。反演得到的平均 Moho 深度约为 43.1 km（图15），最大/最小深度分别为 −3.8 km 与 −76.2 km，与 Hikida & Wieczorek (2007) 及 Wieczorek 等 (2013) 的结果一致。

反演的月球Moho深度
图1　反演的月球 Moho 深度分布（引自论文 Figure 15，为论文核心结果之一）

平均 Moho 深度约 43.1 km，最大/最小分别为 −3.8 与 −76.2 km。月海盆地区域呈现较浅的 Moho（对应大型撞击抬升的地幔），远侧山区则较深，反映月球地壳在正反两面与撞击盆地区域的显著差异。

作者进一步用 LOLA 地形叠加该 Moho 模型得到月球地壳厚度（图16）：平均地壳厚度约 42.0 km，最大值 82.2 km（集中于远侧山区），最小值约 0.29 km（出现在莫斯科海 Moscoviense）；此外 Crisium、Humboldtianum、Grimaldi 与 Apollo 等盆地的地壳厚度不足 5 km。这一模型与 Wieczorek 等 (2013) 的 GRAIL 结果高度一致，两者的 Moho 深度差异服从平均值为 −27.9 m、标准差 2.17 km 的近正态分布（图17），验证了结果的可靠性。

月球地壳厚度
图2　基于 LOLA 地形与上图中 Moho 模型得到的月球地壳厚度（引自论文 Figure 16，为论文核心结果之一）

平均地壳厚度约 42.0 km，远侧山区最厚（约 82 km），大型月海盆地地壳显著减薄，莫斯科海下方最薄处仅约 0.29 km。结果与 Wieczorek 等 (2013) 的 GRAIL 模型高度一致。

讨论：哪些参数在"卡脖子"

作者特别强调两个关键参数对结果的影响：地幔密度与地壳/月海密度。由于残差重力异常不直接依赖地幔密度，若地幔密度取更大值，反演的 Moho 起伏会更大；考虑到平均 Moho 约 43 km 且至少有一个大型撞击盆地已穿透整层地壳，作者折中取 3356 kg/m³。同样，月海玄武岩与地壳密度的不确定性直接影响 Moho 起伏的估计——更好的密度约束将带来更可靠的月球 Moho/地壳厚度模型。

结论与展望

本文提出的方法把**球面三角剖分（STT）**作为球坐标系下构建 3D 密度模型的基本表示——它支持构造完全三维、可局部加密（变分辨率）的模型，能灵活适应不规则数据覆盖与特定研究区域；在前人工作中（Zhang & Chen, 2018；Zhang 等, 2018）的正演基础上，本文给出了完整的面反演算法。除重力数据及其径向梯度外，正演与反演算法同样适用于重力场的其他分量，且通过改写目标函数与偏导数，多种分量的联合反演也很容易实现。该方法不仅复原了月球 Moho 起伏，还揭示出大型撞击盆地与正反两面之间显著的地壳厚度差异，凸显了月球地质演化的鲜明特征。由于框架的一般性，它同样适用于地球乃至其他行星的诸多区域性与全球性界面反演问题。作为开放可复现的技术路线，这一思路后续亦可与更高分辨率观测和更细的界面模型结合，获得更精细的结果。

延伸阅读与引用

  • Wieczorek, M. A., & Phillips, R. J. (1998). Potential anomalies on a sphere: applications to the thickness of the lunar crust. JGR, 103(E1).
  • Hikida, H., & Wieczorek, M. A. (2007). Crustal thickness of the Moon: new constraints from gravity inversions using polyhedral shape models. Icarus, 192(1).
  • Wieczorek, M. A., et al. (2013). The crust of the Moon as seen by GRAIL. Science, 339.
  • Uieda, L., & Barbosa, V. C. (2017). Fast nonlinear gravity inversion in spherical coordinates with application to the South American Moho. GJI, 208(1).
  • 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. GJI, 215(1).

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