知星悟地

基于地震/岩石学先验与卫星重磁观测的球谐域联合反演，重建全球岩石圈密度与磁化强度结构。

本分类将分享教学相关的博文：课程资料、教学笔记与技术教程等。 内容整理中，敬请期待。

本分类将分享与个人兴趣爱好相关的博文。 内容整理中，敬请期待。

论文信息：Zhang, Y., Xu, Y., & Yang, B. (2021). Lévy Gradient Descent: Augmented Random Search for Geophysical Inverse Problems. Surveys in Geophysics, 42, 1109–1132. https://doi.org/10.1007/s10712-021-09644-6 地球物理反演的难点在于：一方面，位场、电磁、地震波等观测的解析都依赖对地下模型的"正演"，而反演往往是非线性、欠定（ill-posed）的；另一方面，如何评估反演结果的可信度——即给出模型参数的不确定度——始终是实际工作中的难题。本文以一作张壹为代表提出了一种新的随机搜索算法 Lévy 梯度下降（L-GD），它把自然界中常见的**莱维飞行（Lévy flights）**思想引入梯度下降框架，试图同时获得传统两类反演方法的优点：梯度法的效率，以及随机/贝叶斯方法的全局收敛性与不确定度估计。下面围绕这篇 Surveys in Geophysics 论文，梳理其动机、方法、数值实验与结论。 研究背景：确定性反演与概率反演的取舍 把地下空间离散成网格，构造"数据拟合 + 模型约束"的目标函数，反演问题就转化为一个最优化问题。按求解算法，主流方法大体可分两类： 确定性反演：基于目标函数对模型参数的梯度的算法，如经典梯度下降、共轭梯度（CG）、L-BFGS 等拟牛顿类方法。它们对连续凸问题极其高效，广泛用于大规模 2D/3D 反演（重力、磁法、MT、地震走时等）。但面对非线性、病态的非凸问题，这类算法容易收敛到局部极小值（或鞍点），结果强烈依赖初始模型；更重要的是，它们通常不提供模型参数的不确定度。若要估计误差，只能在线性化、局部子空间或修改正则化约束等近似下进行，统计有效性常受到质疑。 概率反演（贝叶斯反演）：把模型参数视为概率分布，用 MCMC 等采样技术探索后验分布。它能天然给出误差分布，但计算量往往比梯度法大一到两个数量级，且随参数个数急剧增长，难以直接用于大规模 3D 反演；通常也不直接给出一个"最优模型"。 于是自然产生一个设想：能否设计一种方法，既像梯度法那样高效、可直接用于大规模反演，又具备随机/贝叶斯方法的全局收敛性与不确定度估计能力？本文提出的 L-GD 正是对这一问题的尝试。 方法：用 Lévy 飞行驱动梯度下降 L-GD 的基本思想非常直观：与经典梯度下降一样，每次迭代都沿目标函数梯度给出的下降方向搜索；唯一的区别在于，它不使用固定（或逐渐减小的）步长，而是从 Lévy 分布中随机抽取步长。 莱维飞行是一类特殊的随机游走：路径由大量小步长与少量大步长组成（图1）。这种运动模式广泛存在于动物觅食、流体动力学、光的输运、人类出行乃至地震时空分布等领域。与著名的布朗运动（步长服从正态分布、覆盖面积小）相比，莱维飞行在步数相同时能覆盖大得多的区域，因而对稀疏目标的搜索效率更高。 关键的性质在于：莱维分布的最大值是无穷大的。这意味着即使偶然陷入局部极小值，L-GD 在有限次搜索内也有概率通过一个很大的步长"跳出"局部极值——这正是其全局收敛性的来源。同时，由于搜索过程本质上是一个马尔可夫过程，随着搜索步增长，它会收敛到与初始模型无关的某种后验分布，因此搜索路径的统计量可以用来评估反演模型的不确定度，而不需要额外的后验采样。 算法对外输出三类结果：平均模型 m_mean、标准差 m_sd（即不确定度）以及历史最优模型 m_best（使目标函数最小的那个）；并设有停滞检测（stagnation test）以避免梯度消失导致算法卡在鞍点。为保证不同量纲、不同尺度的参数能同时反演，算法对参数空间做了归一化（尺度因子 p）。 莱维步长采用Mantegna 算法生成，其中指数因子 𝛽 ∈ (1,2) 控制分布形态（𝛽 越大分布越快趋近高斯过程）。数值统计实验表明：𝛽 越大，最大和平均步长越小、路径越稳定；𝛽 越小，越容易产生超大步长，对参数空间的穿越效率越高（表1中最大步长的均值跨了约三个数量级）。因此 𝛽 就像一把"旋钮"，用来在路径稳定与搜索强度之间权衡。加上一个步长尺度因子 𝛼（默认约 0.01–0.1），用于整体缩放步长、调节收敛速度。 ...

论文信息：Li, D., Liang, Q., Du, J., Sun, S., Zhang, Y., & Chen, C. (2020). Transforming Total‐Field Magnetic Anomalies Into Three Components Using Dual‐Layer Equivalent Sources. Geophysical Research Letters, 47, e2019GL084607. https://doi.org/10.1029/2019GL084607 磁法勘探中大都是单分量、标量的观测，而地下磁性体的信息其实藏在更完整的矢量场里。本文由中国地质大学（武汉）与华中科技大学、浙江大学等单位的同行完成，本文作者之一张壹参与了其中浙江大学一方的合作。我们提出的双层等效源方法，用"浅层+深层"两层等效源结合预条件反演，把日常测量得到的总场磁异常 ΔB 稳定地转换为三分量磁场（Bx、By、Bz），尤其在长波长场上精度得到明显提升。下面围绕这篇 GRL 研究短文，梳理它的动机、方法与主要结论。 背景与动机：标量 ΔB 的局限 野外磁测（地面、航空、海洋）最常见到的输出是总场磁异常 ΔB——扣除背景地磁场 B0 后，反映地下磁性体分布的那个标量。ΔB 虽然对磁性体敏感，却同时强烈依赖于岩石磁化强度、磁性体的产状方向以及所处的纬度。换句话说，同一个地下模型在不同位置、不同方向的测量下，ΔB 的形态会有很大差异，这给数据处理与地质解释带来麻烦。 一个行之有效的破解办法，是把标量的 ΔB 转换为磁场的三个方向分量 Bx、By、Bz。矢量化的磁场一方面能压低"磁化方向对异常形态的影响"，一方面能压低反演与解释中的非唯一性，为后续正演与反演提供更干净的输入。 问题在于：传统转换手段（如快速傅里叶变换 FFT）要求观测位于水平平面且数据点在规则网格上。可实际野外数据往往落在起伏地形、不规则分布的测点上，FFT 类方法就"水土不服"了。而等效源方法正好能在任意曲面、任意分布上处理位场数据——这是它被广泛采用的根本原因。 方法：从单层等效源到双层等效源 等效源的思想很简单：在观测面之下布设一层假想的离散小体（等效源单元），通过反演确定每个单元的物性参数（这里为磁化强度），使得这一层假想源正演出的场能拟合观测数据；此后在这个等效源上做任意正演，即可得到上延、下延、分量转换等结果。 然而传统的单层等效源有一个致命短板：经验上（Dampney, 1969）等效源层应放在观测点网格间距的 2–6 倍深度处，它主要用来刻画短波长信息；而位场异常（含磁异常）总是同时包含长、短波长成分。一个贴近地面的单层等效源，无法充分捕获由深部源产生的长波长异常，于是当把 ΔB 转换为三分量、尤其做向上延拓或提取长波长分量时，误差就会显著放大，边界效应也难以压制。 为此本文引入双层等效源（图1）： 浅层：紧贴观测面下方、随地形起伏的曲面/层面，由小而密的矩形棱柱单元组成，负责恢复短波长场； 深层：位于地壳较深位置、置于一个平面上，单元尺度更大，负责捕获长波长信号。深层深度是一个关键参数，可通过磁数据的功率谱分析估计出一个平均深度（本此研究实测中由谱分析定在 5 km）。 图1　合成数据上三种方法重建三分量的绝对误差对比（引自论文 Figure 4，为论文核心结果图）第1–4列分别为 ΔB、By、Bx、Bz 的重建误差分布；(a)–(d) 单层等效源、(e)–(h) 双层等效源（无 PCG）、(i)–(l) 双层等效源（含 PCG）。可见单层源在三分量上误差显著，加入深层源后明显改善，再引入预条件（PCG）后边界效应进一步被压制。 ...

论文信息：Xu, Y., Zhang, Y., Yang, B., & Bao, X. (2022). Phanerozoic Evolution of Lithospheric Structures of the North China Craton. Geophysical Research Letters, 49, e2022GL098341. https://doi.org/10.1029/2022GL098341 克拉通岩石圈如何被改造，是固体地球科学长期争论的话题。华北克拉通（NCC）凭借几十年密集的地震学、年代学与岩石学投入，被公认为研究克拉通岩石圈改造机制的最佳天然实验室。本文由浙江大学徐义贤团队完成，本文作者之一张壹（第二作者）参与本文的数据整理。作者提出了一种从现今岩石圈结构反演古生代岩石圈结构的均衡方法，给出了华北克拉通显生宙岩石圈演化的非均匀图景，并揭示了由深部热流驱动的岩石圈"冷却—加厚 / 加热—减薄"自激振荡行为。下面围绕这篇 GRL 论文，梳理其动机、方法与主要结论。 研究背景：克拉通岩石圈减薄机制之争 克拉通是长期构造稳定、自约 25 亿年以来几乎未变形的古老大陆。华北克拉通在约 480 Ma 时仍保持稳定，其基于金伯利岩中深源捕虏体确定的岩石圈厚达约 200 km；然而大量研究表明，其岩石圈已遭受显著改造——尤其是东部陆块（EB）的岩石圈已减薄至仅约 70 km，远薄于典型克拉通岩石圈的厚度。 对克拉通岩石圈改造的驱动机制目前并无共识：一种观点强调周缘板片俯冲（如华北周缘俯冲）与地幔柱上涌引发的热机械侵蚀（thermomechanical erosion）是主因；另一种则认为，仅由地球稳定的传导热流出地即可驱动岩石圈经历冷却/加厚与加热/减薄的周期性演化（Houseman & Houseman, 2010 数值试验）。长期以来缺少在改造发生前重建岩石圈结构的有效方法，导致各机制在"被摧毁"的克拉通中所起的作用难以厘清。本文正是针对这一空缺，首次尝试把评估造山带塌陷历史的均衡方法推广到可被视为古老造山带的克拉通上。 方法：由现今结构反演古生代岩石圈结构 论文的核心是一套基于均衡与稳态热传导的分析方法。研究区选在华北克拉通核心区（108°–118°E，34°–40°N），现今地壳厚度与岩石圈厚度由 1732 个台站的 P、S 波接收函数成像提供，现今莫霍温度取自已有热流—岩石圈结构研究。 作者沿用 Sandiford & Powell（1990）的 fc − f l 参数化：令 fc、f l 分别为变形前后地壳、岩石圈厚度之比，则岩石圈的演化可由地温场与重力场所控制的高程变化、形变后各深度的温度来刻画。在此基础上，作者扩展了 Nelson（1992） 由均衡理论建立的"造山带塌陷前后结构联系"公式（本文还修正了其附录中缺失的一项），从而在重力均衡与稳态热传导同时满足的假设下，把地壳加厚/减薄与岩石圈加厚/减薄关联起来。当岩石圈改造前的壳-岩比 ψ0 已知（本文取 0.267，接近稳定北美克拉通平均值）时，即可解出改造前的古生代地壳厚度与岩石圈厚度，以及恢复后的莫霍温度。 ...

论文信息：Zhang, Y., Sun, S., Mooney, W. D., & Xu, Y. (2026). Lithospheric magnetization derived from petrological and satellite constraints. Journal of Geophysical Research: Solid Earth, 131, e2025JB032111. https://doi.org/10.1029/2025JB032111 岩石圈磁化强度是理解地球磁场的岩石圈分量、反演地下物质组成与热状态的关键参数。然而，由于观测手段的固有局限，我们长期只能"看到"它的一小部分。本文以一作张壹带领的研究为代表，提出了一种岩石学先验 + 卫星磁测数据联合反演的新框架，重建了覆盖全波段的全球岩石圈感应磁化强度（以垂直积分磁化率 VIS 表征）分布。下面围绕这篇 JGR: Solid Earth 论文，梳理研究动机、方法要点与主要结论。 为什么长波长磁化结构一直难以获得 卫星磁测（CHAMP、Swarm 等）为我们提供了全球尺度的岩石圈磁场观测，如 LCS-1、CHAOS、MF7、CM6 等球谐（SH）模型。但这些模型的功率谱（Lowes, 1974）清楚地表明：岩石圈磁场只有在球谐阶数 16 阶及以上才主导观测量；在更低的阶数（长波长部分），信号被强大得多的地核磁场完全淹没（图1）。这意味着，仅凭卫星数据，我们无法可靠地分离出长波长的岩石圈磁场。 图1　各磁场模型的功率谱（SH 1–80 阶，400 km 高度，引自论文 Figure 1）岩石圈磁场的信号在 SH 16 阶及以上才占主导，更低阶的信号被地核磁场掩盖。这决定了单纯卫星反演无法恢复长波长磁化结构。 此外，还存在一类被称为磁湮灭子（magnetic annihilator）的问题——某些磁化强度分布（例如以固定磁化率磁化的球壳）不产生任何可观测磁场。它使得岩石圈磁化的全球平均值及其长波长特征本质上无法仅靠磁场观测确定。因此，非磁学手段（即地质—岩石学约束）对于确定全球平均磁化强度和长波长结构是不可或缺的。 方法：如何用"岩石学 + 卫星"补齐波长缺口 针对上述问题，本文采用球谐域反演策略，把长波长与短波长分别交给两类不同的约束： 长波长（SH 0–16 阶）：由岩石学先验模型 SM3-SI（Hemant & Maus, 2005）提供约束。该模型基于世界地质构造图、主要岩石类型的实验室磁化率测量以及地壳地震厚度建立，可以给出可靠的全球平均值与长波长磁化变化，从而规避磁湮灭子问题。 短波长（SH 16–80 阶）：由卫星磁数据（采用 CHAOS-8 模型）约束，卫星数据在这些阶数上分辨率高、可信度强。 反演在单个球面等效源层上进行（球面三角剖分生成约 81,920 个三角棱柱单元，垂向积分磁化率 VIS = 磁化率 × 层厚）。目标函数由数据拟合项与模型约束项之和构成，采用 Lévy-Gradient Descent（L-GD） 随机搜索算法迭代求解——该算法不仅能收敛，还能顺带给出反演模型的不确定度估计。合成数据实验表明，恢复的磁场在各波段误差均小于 1 nT，验证了反演框架的有效性。 ...

论文信息：Zhang, Y., & Xu, Y. (2025). Stochastic multi-observables inversion for the 3D thermochemical structure of lithosphere in spherical coordinates: Theory and synthetic examinations. Journal of Geophysical Research: Solid Earth, 130, e2024JB029717. https://doi.org/10.1029/2024JB029717 岩石圈与上地幔中岩石的物性，由它的组成（化学成分）与所处环境的温度、压力状态共同决定，合称为热化学结构。它是认识地球深部物质组成、动力学状态及演化历史的核心依据。然而，温度与组成在控制物性上存在强烈的"等效性"——不同温度—组成组合可产生几乎相同的密度、波速等物性，因此单靠某一种地球物理观测很难把两者区分开。针对这一问题，本文以第一作者张壹等为代表，提出了一套在球坐标下用多类地球物理观测联合反演 3D 热化学结构并同步给出不确定度的框架：基于四面体自适应网格，正演重力、大地水准面、高程、波速、密度与静岩压力等多种观测量，再用随机优化算法求解并得到模型与其估计误差。下面围绕这篇 JGR: Solid Earth 论文，梳理其动机、方法要点与合成检验结果。 为什么需要在球坐标下做热化学反演 上地幔热化学性质的研究区域常常横跨数个经度、纬度，直达大陆尺度。在这种尺度上，地球的曲率变得不可忽略，直接决定了几何建模与地球物理观测正演的精度。然而，此前广泛采用的多观测概率反演方法（如 Afonso 等基于 Fullea 等开发的框架）虽然能同时利用面波频散、体波到时、重力、大地水准面与大地电磁等观测，却大多在笛卡尔坐标下进行，且作为贝叶斯框架内的随机方法计算代价极高——随着未知量维数增长，所需采样时间呈指数膨胀。当 3D 反演的未知量动辄超过数百万个时，这类方法便难以支撑大尺度应用。作者由此提出要在球坐标球壳空间内，构建一套兼顾几何精度、灵活性与计算效率的可扩展反演方案。 方法：四面体自适应网格与技术要点 整个方案围绕三条主线展开： 自适应四面体网格：建模空间采用球壳，网格用四面体而非规则的球面网格剖分。理由很清楚——随半径减小球壳体积收缩，规则网格要么因浅部需要细网格而在深部堆积大量不必要的单元，要么为了效率牺牲浅部分辨率；而非结构化四面体网格可以自适应改变单元大小，在保证效率的同时兼顾分辨率，并能高保真地刻画起伏界面（如 Moho、岩石圈—软流圈边界 LAB 以及板内异常体），还避免了球坐标两极点处的单元畸变，也便于与许多基于非结构化网格的地球动力学软件协同。 岩性热物性计算：采用成熟的 CFMAS（CaO–FeO–MgO–Al₂O₃–SiO₂）矿物体系（约占地壳与上地幔 98 wt%），以镁数 Mg# = MgO/(MgO+FeO) 刻画组成变化（肥沃地幔 Mg#≈89、亏损地幔 Mg#≈94）。针对抽样 Mg#，利用矿物数据库的氧化物统计关系组装代表性全岩组成，在给定温度—压力（400–2,200 K、1–15 GPa）下通过 Gibbs 自由能最小化（Perple_X 实现）确定平衡矿物组合，再依 Stixrude & Lithgow-Bertelloni 的热物性公式与 Voigt-Reuss-Hill 平均计算波速、密度等物性，并建成热物性参考网格（T–P–Mg# 三维插值），从而快速查算物性及其对热化学条件的偏导。合成结果显示密度约 3.0–4.0 g/cm³、VS 约 3.8–6.0 km/s、VP 约 6.0–10.6 km/s，与前人结果接近。 密度—压力耦合与正演：由于密度依赖于压力，作者用有限元方法把静岩压力作为密度—压力耦合问题的解求出，确保模型内部自洽；温度场在传导主导区用有限元求解三维稳态热传导方程、在次岩石圈对流区用线性插值；重力与大地水准面用四面体的多面体重力解析式（Werner & Scheeres 方法）正演，高程按均衡（等静压）思想计算。 随机优化：同时给出模型与不确定度 反演被当作一个多任务评估问题：对地震、重力、大地水准面等每种观测定义 L2 范数的数据失配函数，叠加平滑约束与参考模型约束后构成总体目标函数。求解采用作者此前提出的 Lévy 梯度下降（L-GD）随机优化算法——它结合梯度下降与 Lévy 飞行（多数短步、偶发长跳），对非线性、非凸问题有良好的全局收敛特性，且相比常用随机优化搜索效率更高。更重要的是，L-GD 在寻找最优解的同时，会在解附近随机采样，从而顺带给出反演模型的不确定度估计。为平衡多种观测的贡献，各数据失配函数的权重由深度学习中的 GradNorm 算法迭代确定，使各观测项保持接近的收敛速率。 ...

论文信息：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）。 ...

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