与
是否具有相关性。首先假设两者存在线性关系, 使用标准的高斯线性回归方法[14]为两者建模:述评
北京大学学报(自然科学版) 第61卷 第6期 2025年11月
Acta Scientiarum Naturalium Universitatis Pekinensis, Vol. 61, No. 6 (Nov. 2025)
doi: 10.13209/j.0479-8023.2025.021
摘要 从反演模型的角度, 阐述贝叶斯方法的基本原理、计算手段、发展历程和研究现状。详细介绍各种贝叶斯方法及其在地质学与地球化学各领域的应用, 并探讨未来发展方向以及需要解决的难点。
关键词 贝叶斯方法; 模型反演; 随机模拟; 先验概率; 后验概率; 似然函数
对地球科学问题的建模, 需要对所研究内容的物理、化学等过程进行数学描述, 从而实现问题的定量化理解。这样的建模通常分为模型正演和模型反演[1]两种策略。以地球动力学研究为例, 研究者期望对地幔对流的过程进行建模, 从而对地幔对流动力学过程(如流体动力学、地幔的流变学性质以及能量输运和相变等过程)进行数学描述, 获得对地球内部与外部物质交换及能量演化的理解[2]。建模完成后, 指定模型参数(如黏度等流变学参数以及热扩散系数、相变潜热等热力学参数), 将其输入演化模型, 可以得到与实测数据(如板块活动的速度以及俯冲带发生的样式、地表重力数据、地球温度场等)对比的结果, 这就是地球动力学的正演模型。抽象地说, 将科学问题建模并参数化→指定模型参数→计算模型推测结果, 这样的建模策略称为模型正演。
模型正演要求指定模型参数, 很多情况下, 这却是我们知之甚少的部分。仍然以地球动力学研究为例, 目前对地幔流变学性质的理论约束和实验约束较为缺乏[3], 自然而然的思路是, 将动力学模拟的实验结果与实测数据进行对比, 以两者的一致性作为探索地幔流变学参数的标准, 用模型推测结果来约束模型的输入[4–7]。这种已知模型及模型输出(或称为约束), 反过来寻求模型输入参数的过程, 称为模型反演[8]。对模型反演的需求普遍而迫切, 例如, 无论是从地球化学指标中还原大气氧化、地幔脱气以及大气–海洋相互作用等过程[9–10], 还是从造岩矿物的地球化学成分和模型年龄分布中还原大陆地壳生长过程[11–13], 这些地球科学问题的本质都是用模型推测结果来约束模型内部参数的模型反演问题。
求解一个反演问题, 往往需要通过特定的策略多次调整模型的输入参数, 多次运行正演模型, 直到模型推测结果与观测结果达到良好的匹配度。因此, 对于同样的模型, 反演模型的计算复杂度往往远高于正演模型[8]。通常, 传统的模型反演往往基于优化的方法, 而本文介绍的贝叶斯方法是一种不同思路的模型反演方法。我们以单变量的高斯线性回归为例, 简要介绍反演模型的流程, 以便引出贝叶斯反演方法。
假设我们需要探索某地质样本中两种地球化学指标(共 n组数据)与是否具有相关性。首先假设两者存在线性关系, 使用标准的高斯线性回归方法[14]为两者建模:
, (1)
其中, α 和 b 分别为线性回归的斜率与截距; σ2 为数据的方差; N(0, σ2)是均值为 0, 方差为 σ2 的正态分布。两组数据遵循的生成模型就是式(1)所示高斯线性模型。模型输出(即用于约束反演结果的观测值)是 n组样本{xi, yi}i=1,2,…,n, 模型参数(即反演目标)为{α, β, σ2}。
传统的反演方法往往采用优化的思路。在此思路下, 反演的目标是找到一组唯一的
。这里的
表示最优, 可以通过最小二乘法得到[8]。值得注意的是, 由于最小二乘法基于优化思路, 其背后隐藏着一条假设, 即存在一组这样的参数组合, 并且是唯一的(图 1(a))。
贝叶斯反演则不同, 它认为模型的输入参数(也就是这里的{α, β, σ2})有多种可能的取值组合。贝叶斯反演的目标不是找到唯一的一组最优值, 而是找到一系列可能的参数取值(图 1(b))及其概率分布(图 2)。本文中所有随机采样均通过随机数生成器获得。
这里的概率分布被称为后验概率(posterior pro-bability), 表示贝叶斯反演模型推测结果中, 不同参数取值出现的概率。事实上, 通过优化策略的最小二乘法得到的最优参数取值{α, β, σ2}, 也可以从数学的角度将它等价地看做一个概率分布, 只不过这种概率分布比较特殊, 仅仅这组最优取值有概率, 并且概率为 1[8], 即
, (2)
其中, p(θ)是优化方法的后验概率密度函数(posterior probability density function, posterior PDF); δ(·)为狄拉克 δ 函数(Dirac delta function[8]); 
表示优化模型反演得到的最优参数取值, 本例中
。因此, 可将贝叶斯反演视为优化反演的一种拓展。
在地球科学领域开展贝叶斯模型反演的应用研究时, 虽然不需要深入地了解其细节的数学原理, 但在解决科学问题时, 我们对贝叶斯反演模型的认识直接影响数据的解读方向是否合理, 因此相关理论与概念十分重要。对比优化反演结果(图 1(a))与贝叶斯反演结果(图 2)可见, 贝叶斯方法的最大优势在于信息量丰富。相比传统的优化思路中单一的最优参数取值, 贝叶斯方法利用后验分布, 可以得到一组构成特定概率分布的模型参数的可能取值, 我们可以从这个分布中获得参数取值的趋势和不确定度等多种信息。
此外, 基于优化的反演方法仅仅适用于诸如线性回归这样的简单模型, 但实际问题中的观测数据约束不足或模型具有强烈的非线性时, 反演结果往往具有多解性[8]。此时基于优化的反演方法常常遇到难以找到全局最优点的数学难题, 而贝叶斯方法可以天然地避免这样的问题, 这也是该方法的一大 优势。
总的来说, 贝叶斯方法在地球物理领域的应用较为成熟。近年来, 由于该方法的结果信息量丰富以及可避免优化方法遇到的多解性问题等优势, 地质与地球化学研究领域开始关注贝叶斯反演方法, 并且将其应用到诸如地球化学数据拐点检测[15–16]、热年代学分析[17–19]以及地球动力学反演[4–7]等分支领域。
贝叶斯反演方法完全基于贝叶斯公式[14]:
, (3)
其中, θ 为模型参数, Y 为观测样本或约束。在上述线性回归例子中, 观测样本 Y={xi, yi}i=1,2,…,n , 模型参数 θ={α, β, σ2}。
式(3)中 P(θ)称为先验概率(prior probability), 表示使用观测样本进行约束前, 我们对模型参数 θ={α, β, σ2}的认识。顾名思义, P(θ)描述对模型参数的先验认识。例如, 如果我们关心某个模型中某处的地表热流值 q(mW/m2), 基于现有的认识[2], 一种比较合适的做法是将其先验分布设定为在 30～120 范围内的均匀分布(uniform distribution): U(30, 120); 如果我们进一步了解到该处位于地表热流值更低的大陆地壳, 或许会基于这种更精确的信息, 将其先验分布范围调整为 U(30, 90)。
图1 基于最小二乘法优化反演(a)和贝叶斯反演(b)的线性回归分析结果
Fig. 1 Linear regression based on least squares optimization inversion (a) and Bayesian inversion (b)
例如, 在本文的贝叶斯线性回归实例中, 我们结合参数的意义, 将先验概率取为
(4)
其中, HalfNormal(σ2)是均值为 μ=0, 方差为 σ2 的正态分布的正半轴截断。
MAP(maximum a posteriori)是最大后验概率, 表示后验概率取最大值时的参数取值; 后验概率均值表示后验概率分布的期望值; 95%可信区间表示 2.5%～97.5%的可信区间
图2 基于贝叶斯反演的线性回归模型参数的后验分布以及相应的统计量
Fig. 2 Linear regression based on Bayesian inversion, showing the posterior distribution of model pa-rameters and their associated statistics
式(3)中的 P(Y|θ)称为似然函数(likelihood), 用来表示给定回归参数 θ={α, β, σ2}后, 获得样本 Y= {xi, yi}i=1, 2, …, n 的概率。可以直观地看到, P(Y|θ)与正演模型的形式直接相关。例如, 对于式(1)描述的高斯线性回归模型, 似然函数[14]为
(5)
值得注意的是, 似然函数不仅受到正演模型的影响, 还受到预测误差形式的影响。在本例中, 似然函数可以拆分为确定性的正演模型和概率性的预测误差两部分。
1)确定性的正演模型, 也称为正演模型预测值 d=f(θ), 即线性模型
, (6)
其中, f 表示正演模型描述的从模型参数 θ 到预测值d 的映射。
2)概率性的预测误差, 也称为正演模型预测误差
, 即正态误差
(7)
依据具体问题, 确定性的正演函数可能有各种复杂的形式。然而, 对于概率性预测误差的处理, 从下面大量例子中可以看到, 选择一个服从高斯分布的扰动项作为误差的描述是非常普遍的。推广到n 个约束的情形, 用数学语言描述就是
(8)
其中, f(θ)表示基于正演模型 f 和参数 θ, 对模型结果 Y 的正演预测函数; V 是样本协方差矩阵(covariance matrix), |V|表示其行列式的值; n 是样本数据(约束)的个数。
总之, 可以把这种贝叶斯框架下似然函数的定义视为两个步骤: 1)定义正演模型的预测函数 f(θ) (确定性); 2)定义该正演模型给出的预测与观测值的匹配度(概率性, 如高斯分布误差)。
式(3)中的 P(Y)称为归一化常数(evidence proba-bility), 可展开为
, (9)
即表示对所有可能的参数 θ 取值积分后, 获得样本Y 的概率。在涉及众多参数的复杂模型中, 这个积分的计算是非常困难的。事实上, 贝叶斯反演的计算难点几乎完全集中在这一项, 相应的各种算法都是为了解决该积分的计算问题[14]。另一方面, P(Y)这一项只是通过 P(Y|θ)与 P(θ)计算得到的积分常数, 不具有更多的物理含义以及讨论价值, 因此本文不赘述。
式(3)中 P(θ|Y)为前面讨论过的后验概率, 表示在获得样本 Y 的观测后, 对模型参数的新认识。后验概率即为贝叶斯反演的全部目标。
一方面完成对正演模型的定义, 从而确定似然函数 P(Y|θ), 另一方面取定对参数认识的先验分布 P(θ), 便可以进行取样, 获得观测样本 Y, 通过式(9), 积分得到经验概率 P(Y)。进一步, 便可以通过式(3), 计算得到后验概率 P(θ|Y)。例如, 对于本例的贝叶斯线性回归模型, 便可以按照上述方法得到关于斜率 α、截距 β 以及方差 σ2的后验分布(图 2)。
可以发现, 相比于只能得到一组最优选择
的传统反演策略, 贝叶斯方法可以获得更加丰富的信息。利用后验概率分布, 不仅可以进行类似传统优化策略反演的操作, 找到以最大后验概率作为标准的最优参数, 还可以采用更丰富的手段提取后验分布中的信息, 例如对后验概率求期望值, 得到参数的后验均值(这也是在贝叶斯反演中常用的指标, 见图 2):
, (10)
其中, E 表示求数学期望值。需要指出的是, 在很多情况下, 后验分布可能十分复杂, 此时仅仅采用最大后验概率或后验期望值等非常有限的信息作为模型推测结果是不合适的, 没有将贝叶斯反演模型的优势发挥出来。获得后验分布后, 应该综合地考虑整个分布的信息, 例如考察参数分布的更高阶统计量, 如计算参数的 n 阶中心矩:
。 (11)
例如,n 取 2 时, 表示参数的二阶中心矩即方差, 刻画数据后验分布的离散程度; n取 3 时, 表示参数的三阶中心矩, 刻画数据后验分布的偏度(skewness), 即分布的不对称程度。再例如, 通过参数的联合后验分布中彼此之间的相关系数, 定量地评价不同模型参数之间的相关程度[20]。
最后, 我们还可以非常方便地得到一种不确定性的度量——可信区间(confidence interval), 图 2展示 95%可信区间的模型推测结果。需要注意的是, 尽管形式上非常相似, 可信区间这一概念并不等同于经典频率学派统计方法中的置信区间(credible in-terval), 两者的区别与详细的数学讨论可以参考文献[14]。
贝叶斯后验概率的计算是一个复杂的数学问题[14,21–22], 其难点在于, 了解先验概率与似然函数之后, 还需要计算归一化常数 P(Y), 才能得到后验概率, 而 P(Y)往往需要通过复杂的高维积分来获得[14]。对于复杂问题, 需要根据科学目标来选择, 或者设计具体的算法。对于非常有限的一些简单问题, 可以得到 P(Y)的解析形式[14]。绝大多数情形下, 无法得到 P(Y)的解析形式, 此时主流的处理方式是基于采样的随机模拟方法[14], 即不计算后验概率值, 而是通过在这个后验分布中采样来获得信息。由于没有计算归一化常数 P(Y), 实际上, 采样是在
描述的分布中进行, 其中 Z 为常数[14]:
。 (12)
以使用随机模拟方法(Monte Carlo)计算 π 值为例, 引入拒绝–接受采样法(reject sampling)。取一个圆以及它的外切正方形, 在正方形中进行几何均匀分布随机地投点(图 3(a))。根据简单的几何概型, 某一次投点恰好在圆中的概率 P 为
。 (13)
根据大数定律, 进行足够多次投点后, 可以用投点在圆中的次数 N 圆 与总投点数 N 总数的比值作为 P 的近似值:
。 (14)
因此, 当 N 总数足够大时, 可以使用如下公式计算 π的近似值:
。 (15)
上例是典型的拒绝–接受采样, 下面将拒绝–接受采样法直接应用到贝叶斯后验分布采样的问题上。假设对待求解的贝叶斯反演问题已经完成采样, 获得样本(约束)Y以及对先验概率 P(θ)和似然函数 P(Y|θ)的定义, 于是, 得到两者的乘积 P(q)· P(Y|q), 即图 3(b)中斜线填充部分。根据式(3), 在确定样本 Y 后, 乘积 P(q)·P(Y|q)与后验概率 P(θ|Y)之间只相差一个常数因子。
选取一个适当的提议分布(proposal distribution) q(θ), 保证该概率分布的范围与先验分布相同, 并取一个常数的缩放因子 M, 保证
。 (16)
如图 3(b)所示, 式(16)的几何含义是对已知概率分布 q(θ)进行 M 倍缩放, 保证其取值始终大于 P(q) ·P(Y|q), 从而进行“覆盖”。这一步与计算 π 值的例子中用一个正方形来覆盖圆的想法是一致的。定义接受率:
。 (17)
单次采样按照如下步骤进行: 首先对我们给出的提议分布 q(θ)采样, 得到第 i 次采样值 θi。计算此次采样值对应的接受率:
。
接下来, 从一个均匀分布 U(0, 1)中采样, 得到采样值 ui。若 ui<αi , 则接受这次的采样值 θi , 否则将其舍去。经过对 q(θ)采样, 并基于均匀分布的拒绝–接受采样, 两者的复合等同于在后验分布
中采样。每一次采样是独立的, 经过足够多次采样后, 就得到服从后验分布
的参数样本。
图3 拒绝–接受采样法计算p值(a)和贝叶斯后验分布(b)示意图(修改自文献[14])
Fig. 3 Sketch maps of rejection-acceptance sampling method for calculating the value of p (a) and the Bayesian posterior distribution (b) (modified from Ref. [14])
需要指出的是, 此方法是一种随机模拟方法, 并没有真的计算出每一种参数取值的概率。我们所做的, 是通过一个复合的采样(先对已知分布 q(θ)采样, 再根据一个均匀分布决定是否接受)得到后验分布, 然后重复多次, 就可以用采样频数来逼近真实的后验分布。
拒绝–接受采样法的问题在于, 如果提议分布q(θ)选取不适当, 可能导致采样过程中接受率 α 非常低, 大量采样被拒绝, 从而造成采样的低效, 可见提议分布的选取对采样性能的影响很大[14]。无论如何, 作为随机模拟方法计算贝叶斯后验分布的基本方法, 拒绝–接受采样法为后续算法提供了对后验分布进行采样的基本思路。
马尔可夫链蒙特卡洛方法(Markov Chain Monta Carlo, MCMC)也是一种不计算归一化常数 P(Y)而直接对后验采样的方法。20 世纪 90 年代以来, 马尔可夫链蒙特卡洛方法已成为探索复杂贝叶斯模型的主要方法[23]。
与拒绝–接受采样法的主要区别在于, 马尔可夫链蒙特卡洛方法的每次采样不是独立的。通过构建一个马尔科夫链(Markov Chain) [24], 使得模型可以在参数空间随机游动, 从而更加精细地探索参数空间的局部结构, 实现更高的采样效率。在高维参数空间中, 马尔可夫链蒙特卡洛方法的表现通常优于拒绝–接受采样法。
马尔可夫链蒙特卡洛方法由 Metropolis 等[25]于1953 年提出, 并且由 Hastings[26]于 1970 年完善。早期的马尔可夫链蒙特卡洛方法被称为 Metropolis-Hastings(M-H)算法, 构建了后续各种马尔可夫链蒙特卡洛方法的基本框架。与拒绝–接受采样法对比可以发现, Metropolis-Hastings 算法同样具有提议分布和接受率的概念。另一方面, 以马尔可夫链为基础的马尔可夫链蒙特卡洛方法需要解决收敛性问题, Metropolis-Hastings 算法奠定了基本的理论框架, 引入细致平衡(detailed balance)的概念来为马尔可夫链收敛提供理论依据[14]。
之后, 马尔可夫链蒙特卡洛方法的算法探索与应用蓬勃发展。1984 年, Geman 等[27]首次提出适合高维采样的新算法——Gibbs 采样法, 并将其应用到具体的统计问题中。1990 年, Gelfand 等[23]进一步阐释了 Gibbs 采样法的思想内涵, 并提出简单随机效应模型, 使得马尔可夫链蒙特卡洛方法真正成为统计学界大规模使用的贝叶斯模型计算方法[21]。Green[28]于 1995 年提出逆跳马尔可夫链蒙特卡洛 方法(Reverse Jumping Markov Chain Monta Carlo, RJMCMC), 该方法使得贝叶斯模型可以处理变维(transdimensional)问题。
此外, 针对马尔可夫链蒙特卡洛方法的收敛加速以及优化方法也在持续发展中。应用广泛、效果优异的实例诸如引入辅助变量来加速收敛的哈密顿蒙特卡洛模拟(Hamilton Monte Carlo, HMC) [29]、基于多链采样避免探索陷入局部的平行退火(parallel tempering, PT)[30]以及基于积分近似方差减小技术的拉奥–布莱克威尔化(Rao-Blackwellisation)[31]等 算法。
变维, 指建模过程中, 将模型的变量数或者复杂度作为一个待求解的变量进行反演[32]。1995 年, Green[28]提出逆跳马尔可夫链蒙特卡洛方法, 使得贝叶斯模型可以处理变维贝叶斯问题, 极大地拓宽贝叶斯模型的应用领域。Robert 等[21]称这项工作为第二次马尔可夫链蒙特卡洛方法革命的开端(the start of the second MCMC revolution)。逆跳马尔可夫链蒙特卡洛算法的详细介绍可参见文献[28,33–34]。
将变维贝叶斯模型应用到实际科学问题中的动机在于, 在复杂的科学建模中, 我们往往不希望对数据的复杂度和规模信息做出先验的假设, 而是希望通过数据本身来确定。Malinverno 等[35]首次将变维贝叶斯模型应用到地球科学领域, 通过直流电阻率对地层厚度和地层层数信息进行反演推断。Sam-bridge 等[32]2013 年曾介绍变维贝叶斯模型在地球科学中的应用, 但集中在地球物理领域。
接下来, 我们列举一些变维贝叶斯模型在地质学与地球化学领域的应用实例。
Gallagher 等[15]2011 年将变维贝叶斯模型以及逆跳马尔可夫链蒙特卡洛方法应用到地球化学领域, 讨论了分段线性插值地球化学数据的时间演化算法: 假设待研究的地球化学数据(如 206Pb/207Pb 和Eu/Eu*等)在若干时间节点发生跳变, 在其他时刻则保持为常数。在这里, 发生跳变的时间节点个数描述了模型的复杂度。人为指定该“复杂度”显然不合理, 但也无法通过物理观测或地质观测得到。因此在他们的模型中, 用发生跳变的时间节点个数描述的“模型复杂度”被纳入未知变量中, 成为待反演的参数。他们称之为“拐点(changepoint)分析”。
之后, 拐点分析方法被进一步应用到地球科学中更广泛的领域, 如古地磁强度拐点分析[36]、热年代学与热演化历史分析[17–19]、古气候拐点分析[37]以及地球化学样本的进一步应用。例如, Spencer 等[16]利用地球化学指标的拐点判断陆生植物登陆时间。2016 年, Sambridge[38]从时间序列分析(time-series analysis)和回归(regression)分析的角度对拐点分析模型做了总结, 并列举大量研究实例。
除了拐点分析, 变维问题中的复杂度建模还可以有更多的诠释, 例如混合体系的混合组分数[39]、多项式拟合的阶数[40]以及空间几何参数拟合[41]等。这些应用集中在地球物理领域, 若将这种强有力的模型应用到地质与地球化学领域, 相信也会是非常有前景的研究方向。
地幔是地壳物质与能量交换的储库, 地壳与地幔的相互作用控制着地球的一阶演化。为此, 研究者从化学组成的角度对地幔进行研究[42–43], 例如对地幔物质端元的分类[43], 今天仍然被广泛应用与讨论。Liu 等[44]在对一个经典地球化学模型(地幔部分熔融模型)的研究中, 主要考虑地幔部分熔融过程中微量元素在熔体与残余部分之间的迁移, 地幔部分熔融结晶产生玄武岩, 残余体构成地幔岩, 微量元素在两者之间的分配受到初始状态、分异过程和分异程度影响。地球化学领域存在多种经典模型来描述这一过程, 包括存在解析解的简单模型(如批式熔融(batch melting)、增量批式熔融(incremental batch melting)、分离熔融(fractional melting)和动态熔融(dynamic melting)、不存在解析解的双孔隙熔融(double porosity melting)以及非平衡熔融(disequi-librium melting)等[45–47])。这些模型潜在的不同表现是该研究的动机之一。
Liu 等[44]考虑了两组正演模型: 1)非平衡分离熔融模型, 模型参数为
, 其中 F 为熔融分数, ϵ 为非平衡系数, ϵ=0 表示平衡情形; 2)非平衡动态熔融模型, 模型参数为
, 其中 α 表示熔体与残余固体的质量通量比值, α=0 表示分离熔融情形。事实上, 对于简单的非平衡分离熔融模型, 由于可获得解析解, 只通过传统的最小二乘反演便可以观测到模型参数 F 与 ϵ 的正相关关系。但是, 对于无显式解的复杂模型(如非平衡动态熔融模型), 优化最小二乘反演难以得到可靠的结果,F 与 ϵ 的正相关关系是否可以推广, 是一个待解决的问题。这也是本文论述贝叶斯反演的动机之一。
Liu 等[44]模型的约束 Y 是大洋中脊橄榄岩中单斜辉石微量元素(包括稀土元素和钇)的组成。简单模型(非平衡熔融模型)的正演预测存在解析解, 复杂情形(非平衡动态熔融模型)则需要求解一组常微分方程。我们将正演预测部分统一记为 d=f(θ), 预测误差依然采用高斯误差的形式:
。 (18)
值得注意的是, 
被定义为
, (19)
其中, 下角标 i 表示第 i 种组分;
表示基于正演模型 f 和参数值 θ, 对第 i 种组分的预测浓度; 
和
分别表示地幔源区某组分的浓度和浓度标准差, Ci 和 δCi 代表样本中的浓度和浓度标准差。相比于传统的高斯误差, 式(19)的分子部分采用了浓度的对数之差, 而非浓度之差; 在分母部分, 除了引入样本的误差, 还加入地幔源区的误差。这些调整引入更多的物理因素, 因此获得更好的反演效果。完成概率建模后, 对后验分布的采样使用 Me-tropolis-Hasting 算法。
Liu 等[44]的模型推测结果展示了在简单模型上与传统最小二乘优化反演结果的一致性。此外, 贝叶斯反演利用联合分布判断参数相关性的便捷性, 同样促进我们对地幔熔融过程的理解。Liu 等[44]发现, 非平衡系数 ϵ 与 143Nd/144Nd 之间无明显相关性, 指示地幔稀土元素分布的不均一性并非控制部分熔融平衡系数的重要因素; 贝叶斯模型还印证, 复杂熔融模型中的熔融分数 F 与非平衡系数 ϵ 依然存在正相关关系, 从而引出物理机制的进一步讨论。用贝叶斯反演的方式研究地幔化学成分还包括 Brown等[48]的工作。
除了地球化学观测的角度, 研究动力学过程也是研究地幔与地表物质及能量交换的重要手段。这方面的正演模型研究以大规模数值模拟为基础, 20世纪 60 年代早期开始兴起[48–50], 通过描述地幔的流体动力学过程、能量输运和热力学过程以及流变学和相变等过程, 模拟全球尺度的全地幔对流情形, 从而获得可与观测结果对比的结果。复杂的正演模型对反演带来极大的困难, 因此在地球动力学领域, 贝叶斯模型反演同样逐渐获得大家的关注。
本文主要关心地球动力学数值模拟中地幔流变学参数的反演问题。仍然从贝叶斯模型反演的框架来看, 正演预测模型 f 是全球的地幔对流动力学模型。这是一个数值模型, 建模过程包括全地幔的速度场、压强场、温度场以及流变学模型, 通过自适应网格求解模型化全球地幔对流, 即三维球壳的有限元离散化守恒方程[51]。需要进行反演的模型参数 θ 是流变学参数, 包括屈服应力 σy、应力应变幂次 n、应变活化能 Ea 以及上下地幔黏度等。模型约束 Y 则是可观测的现今全球板块运动速度分布。至此, 我们将模型的全部正演预测部分置于贝叶斯反演的框架下。
本例的特别之处在于正演模型庞大的计算量。我们知道, 每进行一次贝叶斯反演的后验采样, 都需要根据当前采样得到的参数值 θ 进行一次正演模型预测的计算, 而在本工作中, 这意味着一次复杂的全球尺度地幔对流演化数值模拟, 而单次数值模拟的计算复杂度可以达到千核并行的数量级。因此, 通过拒绝–接受采样或构建马尔科夫链的方式进行随机模拟的方法是不现实的。Hu 等[4]首先假设先验概率和似然函数都是高斯分布, 于是后验概率满足
, (20)
其中,
(21)
其中, Cdata 和 Cprior 分别是似然函数和先验分布的高斯协方差矩阵。进一步地, 他们对 J(θ)在最大后验概率处进行线性化近似。简单地说, 首先基于梯度下降的方法, 找到最大后验概率对应的参数取值θMAP (复杂的梯度信息通过伴随方法(adjoint method)得到[5,7]), 接着将后验分布近似为均值位于 θMAP 处的高斯分布, 则高斯分布的均值为 θMAP , 协方差矩阵为
是线性化的 f [7]。
通过上述方法, Hu 等[4]实现对后验分布的高斯分布近似, 避开了基于采样的贝叶斯后验推断方法的庞大计算量。作为对比, Holbach 等[6]同样对流变学参数进行模型反演, 但采用线性流变学模型, 并且模型的维度只有二维, 其正演模型的计算要求没有那么严格, 因此可以采用基于采样的马尔可夫链蒙特卡洛方法进行后验分布的推断。
最终, 反演模型得到与实验基本上一致的结果, 此外还带来新的物理学思考, 例如应力应变幂率指数 n 与应变活化能 Ea 之间的负相关关系。模型反演得到的不同俯冲带的逆冲断层强度明显不同, 增进了我们对不同地区驱动逆冲断层地震的应力以及控制有效黏度的各种因素及其相互作用的理解。
大陆地壳是地球独特的物质成分[52], 研究其形成及演化机制是地球科学的核心问题, 相关研究成果的综述可参阅文献[53–56]。大量大陆地壳研究相关的实地考察、样本收集与实验提供了海量数据, 为基于统计理论的贝叶斯反演方法提供了用武之地。
2023 年, Pease 等[57]在关于大陆地壳成分反演的研究中采用标准的贝叶斯反演模型框架, 通过观测到的大陆地壳波速等地球物理学信息, 反演大陆地壳不同深度的地球化学成分。其模型参数 θ 除包括大陆地壳中 SiO2, TiO2, Al2O3, FeO, MgO, CaO, Na2O, K2O, H2O 和 CO2共 10 种氧化物外, 还引入压强 p 和孔隙度 Φ 等影响大陆地壳地球物理学性质的参数, 观测约束 Y 则是大陆地壳的密度 ρ 以及大陆地壳的地震波速 VP 和 VP/VS等。因此, 正演模型f对应一个根据大陆地壳的化学成分预测其物理性质的过程。他们采用 Perple_X 开源代码库计算正演模型的预测部分, 这是一个基于岩石中矿物化学成分组合以及压强、孔隙度等参数计算岩石物理性质(包括地震波速等)的软件[58]。基于上述正演模型中预测部分的定义, 他们进一步引入预测误差, 并采用标准的高斯误差:
(22)
其中, Σ 为高斯分布的协方差矩阵。后验分布的推断则采用拒绝–接受采样法。至此, 他们将该工作完全置于贝叶斯反演模型的框架内。利用约束, 模型给出类似先前研究成果的下地壳化学成分, 但是给出的中–上地壳化学成分相对偏基性。相对于之前的工作, Pease 等[57]给出不确定性的后验分布, 并探究了不同的先验假设(无信息先验或采用过去的研究结论作为地壳成分的先验分布)对模拟结果的影响。
事实上, 以 Rudnick 等[52]为代表的大量早期工作给出不同的大陆地壳化学成分分布, 探究这些结论的可靠性是 Pease 等[57]的研究目的之一。Pease等[57]将这些研究结果以及数据库(例如 EarthChem, Crust1.0)提供的化学成分范围(如 SiO2含量)作为模型的先验知识, 探究后验分布对这些先验知识的敏感程度, 并且定量地描述先验知识的不确定如何传播到模型推测结果中, 得到的结论反而倾向于使用更少的先验假设。
值得一提的是, 文献[57]的研究团队与合作者将贝叶斯反演方法应用到大量地球科学相关领域, 包括盆地模型反演[59]、冰川侵蚀模型[60]、大气海洋碳循环模型[61]和小行星热演化模型[62]等。
大陆地壳成分的演化之所以重要, 是因为可以通过相关的质量守恒过程来约束大陆演化的地质过程。这方面的研究成果同样非常丰富[63–66], 其中不乏贝叶斯反演模型的应用。Ptáček 等[66]利用大陆上地壳沉积岩化学成分的演化以及不同大陆地壳地球化学成分端元(酸性岩、基性岩和科马提岩)的化学成分演化, 反演不同地质时期各大陆地壳地球化学成分端元对大陆地壳总体物质成分的贡献, 发现必须有超过 50%的酸性端元贡献, 才能解释 3.5Ga 的沉积岩微量元素化学组分, 指示 3.5Ga 的大陆地壳分异以及板块构造活动已经启动, 并且进入活跃状态。将沉积岩作为混合组分, 通过反演大陆地壳各个成分端元的贡献来理解大陆生长过程, 尽管富有争议[63–67], 依然是近来非常活跃的重要科学问题。无论如何, 贝叶斯反演为我们理解大陆生长这一重要问题提供了有力的数学工具。
另一个理解大陆地壳演化的重要角度是聚焦大陆地壳演化的物理过程, 此时贝叶斯模型反演同样可以发挥作用。Rosas 等[68] 2018 年的工作重点关注大陆地壳, 尤其是地球早期大陆地壳的生长和演化。他们将大陆地壳视为从地幔中抽取出来的一个物质储库, 两者之间存在物质的交换和化学元素的分配。地幔物质流入大陆地壳的过程称为大陆生长(continental crust generation), 反之, 大陆地壳物质流回地幔的过程称为大陆再循环(continental recyc-ling), 两者之间的差值则为大陆净增长(continental growth)[53]。Rosas 等[68]模型中的主要观测约束 Y 是147Sm-143Nd 和 146Sm-142Nd 两种放射性同位素体系, 这些同位素一方面遵循自身的衰变周期, 持续地转化, 另一方面遵循在地壳与地幔之间不同的配分系数, 随着两者的物质交换进行迁移。Rosas 等[68]的反演目标是通过对地质历史时期这两种同位素在地壳中的丰度, 反演大陆地壳与地幔的物质交换过程。因此, 在贝叶斯反演的框架下, 模型参数 θ 主要是大陆地壳与地幔物质交换过程的一些参数化指标, 例如大陆生长和消亡的启动时间 ts 和 tp、大陆体积的衰减速率 κg 以及大陆消亡速率的衰减速率 κr等, 共计 10 个参数, 观测值 Y 则为一组时间点{t}k及其对应的同位素观测值
和
。
值得注意的是, 尽管 Rosas 等[68]提到先验和后验这些概念, 显然运用了贝叶斯反演的思想, 但整个反演过程与贝叶斯的标准框架存在一些差别。对于所有的模型参数 θ, 他们将其指定为给定范围内的均匀分布, 这与贝叶斯反演别无二致。但是, 模型中没有给出似然函数, 而是定义了取定一组特定参数值后, 正演模型预测值与实测值的“全误差” (total misfit):
其中,
,
,
其中, 
(tk)以及
(tk)是 tk 时刻正演模型的预测值。根据同位素数据的特点, Δϵ 和 Δμ 的取值分别为1 和 2。这个“全误差”的引入并非为了给出似然函数和计算后验概率, 而是作为参数的筛选。Rosas等[68]的模型反演过程如下: 在均匀分布的模型参数空间中随机取样, 计算每一个样本的全误差 M, 只保留 M 小于给定标准的样本, 所有被保留的样本构成一个分布, 将其视为“后验分布”(这里使用带引号的“后验分布”, 以示其与贝叶斯意义上后验分布的区别)。
最终, 模型给出的“后验分布”指示大陆地壳的生长在相当早的时间(冥古宙末期, 4.0Ga 左右)已经启动, 大陆地壳的消亡以 2×1022～4×1022kg/Ga 的速率同时启动并逐渐衰减, 表明早期地球板块构造活动存在的可能性以及更加强烈的火山活动与大陆消亡过程。有趣的是, 同为贝叶斯反演方法在大陆动力学中的应用, Ptáček 等[66]的工作与 Rosas 等[68]的工作一样, 都得到类似的相对较早的大陆生长与板块构造启动时间。
从文献[68]中图 2 可以发现, 尽管以全误差 M的筛选作为生成“后验分布”的方式可能缺乏严格的数学基础, 但是其推测结果(强烈倾向于一个非常早期的板块构造启动时间)中蕴含的信息是显著的, 可能归因于模型的正演结构不那么复杂以及数据约束的充足。
Guo 等[69–70] 2020 和 2023 年的两篇文章延续了该思路, 用类似的随机模拟方法在参数空间中采样, 并使用特定的“全误差”进行筛选, 从而得到参数的“后验分布”, 进一步利用 Ar 同位素[69]和 Hf-Nd 同位素[70]对大陆地壳的演化过程进行约束。总的来看, 贝叶斯反演模型为大陆动力学的复杂建模和参数反演提供了强有力的数学工具, 打开了新的思路。
古地磁学涉及高度专业性的实验与样本分析技术和理论基础[71–72], 也极大地推动了地球科学领域重大科学问题的解决。早在 20 世纪 60 年代, 地球磁场古方向的测定就为板块构造理论的建立和发展提供了重要的依据和定量资料[73]。作为一个矢量, 磁场不仅具有方向(地磁古方向), 也具有强度(地磁古强度), 而无论是涉及的理论、实验技术和方法, 还是对样品的要求, 后者都复杂得多[72]。对数据质量和分析方法越来越高的要求, 使得贝叶斯反演方法在古地磁研究领域得到广泛的应用, 如一阶反转曲线(first order reversal curves, FORCs) [74]的计算以及古地磁强度或方向演化曲线的拟合[75–78]等。
Cych 等[75] 2021 的工作要解决的问题是地磁古强度 p(paleointensity)数据的筛选和计算。若要计算特定时间某个地点的地磁古强度数据, 一种主流的做法是对该地点多个样本进行 Arai 图[79]投图。在Arai 图中, 我们会对同一个样本的热退磁过程进行投图并进行线性拟合。理想情况下, 某个样本的所有数据点应该位于一条直线上, 可以利用该直线的斜率, 通过计算得到一个古强度数据 Bm (下角标 m表示第 m 个样本) [72]。
然而, 进行 Arai 图投图的假设可能被多种自然因素或材料性质干扰, 导致 Arai 图上的数据点偏离直线, 以致相应的斜率和Bm值计算结果不准确。为解决此问题, 传统的方法是通过一系列筛选, 将偏离直线较远的样本舍去, 最后使用质量较高的样本来计算该地点的古强度值 Bexp (例如, 对筛选后的所有样本取平均值)。Paterson[80]提出通过 Arai 图中曲线的曲率相关的量 k 进行筛选。Cych 等[75]指出, 比起保留或删除这种“非此即彼”的筛选方法, 应该充分利用每一个样本的信息。他们基于 Tauxe 等[81]提出的经验公式(Arai 图存在对直线的偏离时, Bm与 km存在相关性), 将实际古强度 Bexp 视为所有样本对应的 Bm 与 km 的回归参数:
,
其中, 
服从正态分布 
。Cych 等[75]在贝叶斯方法框架下, 计算每个样本对应的参数 Bm 与 km的后验分布, 然后利用贝叶斯回归模型计算 Bexp 的后验分布。这种方法的本质在于, 在筛选样本时, 抛弃非此即彼的保留或放弃策略, 而是将筛选标准(如这里的 km)也作为信息, 纳入真值 Bexp 的计算中。这种策略充分利用每一个样本的信息, 在数据稀疏的情形下, 大大提升了物理量估计的准确性。值得指出的是, 这种思路可以推广到任何数据筛选的问题上, 其核心想法基于以下 3 个方面: 1)每个数据, 哪怕质量较低, 都存在可以利用的信息; 2)质量低的数据信息“价值”低, 反过来, 质量高的数据信息“价值”高, 这可以通过贝叶斯方法来描述; 3)误差与真值之间同样可以建模, 如本例中 Bm 与 km 的线性关系。
总的来说, 尽管依托于电磁学理论, 古地磁学仍然是一个强烈依赖观测的学科。基于统计理论的贝叶斯方法非常适合这样一个样本数据驱动的学科, 探索贝叶斯反演方法在古地磁领域的应用是一个非常具有前景的方向。
地球化学异常识别(geochemical anomalies de-tection)的目的是通过对于地球化学数据的汇编分析, 发现数据在空间分布上的结构与模式[82]。地球化学异常识别应用广泛, 例如地质矿床勘探[83–84]及构造环境机制重建[85–86]等。一方面, 该领域的数据具有空间分布结构(如地质图和地球化学位图等), 另一方面, 数据的空间分析对数据量及准确性的高要求, 使得不确定的量化同样重要。Huang 等[87] 2022 年利用贝叶斯卷积神经网络(Bayesian Convo-lutional Neural Network, BCNN)进行金矿床勘探的研究。卷积神经网络是一种深度学习模型[88], 作为神经网络的一个分支模型族, 卷积神经网络的特点是可以保留数据的空间结构和特征, 因此广泛应用于图像类型的数据处理领域, 如计算机视觉和图像特征提取。Huang 等[87]将卷积神经网络视为正演的预测模型, 将模型的输入(河流沉积物样本中 10 种元素的空间分布位图)映射为对金矿床的空间分布预测:
, (23)
其中, f为卷积神经网络正演模型; x 为模型输入, 即地球化学元素的空间分布位图; y 为模型预测结果, 即金矿床的空间分布; θ 为卷积神经网络的卷积层参数。
在地质预测的反问题框架中, Huang 等[87]利用已知金矿床分布的训练集{xi, yi}, i=1, 2…, n, 求解卷积层参数 θ, 其中 n 为训练集个数。在传统的卷积神经网络中, 可以通过优化方法得到对 θ 的确定性估计值。Huang 等[87]引入贝叶斯框架, 使得求解卷积层参数 θ 的确定性估计值变为求解 θ 的后验分布
。基于训练得到的卷积神经网络参数后验分布, 对于仅了解地球化学元素空间分布 x 的区域金矿空间分布的预测值
, 可以通过对后验分布求数学期望值(即积分)得到:
。 (24)
前面提到, 在地球化学异常识别研究领域, 不确定性的定量分析是重中之重。将贝叶斯方法引入模型中, 将模型参数处理为概率分布, 拓展了量化不确定的可能性。例如, 文献[87]中的随机不确定性(aleatoric uncertainty, 来自数据)与认知不确定性(epistemic uncertainty, 来自模型)可以直接在贝叶斯方法框架下进行计算和比较。这使得勘探工作得以在地球化学异常值显著且不确定性低的区域优先进行。Huang 等[87]对这两种不确定性在空间上的相关性进行了分析, 发现两者之间存在一定的正相关关系。
贝叶斯方法在机器学习与深度学习的兴起中不仅扮演了理论基础的角色, 还通过模型设计和推断策略的演进直接推动这两个领域的发展。传统机器学习的许多经典模型源于贝叶斯理论。例如, 朴素贝叶斯分类器(Naïve Bayesian Classifier)[89]以其简洁性和强大的概率推理能力, 成为分类任务中的重要工具。在深度学习领域, 贝叶斯方法同样被广泛应用于构建复杂的生成模型和处理不确定性问题。例如, 前面提到的贝叶斯神经网络(Bayesian Neural Networks, BNNs)[90–91]通过将权重视为随机变量, 提供了对模型预测中不确定性的量化。而变分自编码器(Variational Auto-Encoder, VAE)[92]作为生成模型的典范, 通过结合贝叶斯理论中的变分推断(varia-tional inference, VI), 在许多领域取得卓越成果。另一方面, 传统机器学习与深度学习的兴起也促进了贝叶斯方法的发展。深度学习带来强大的计算能力和丰富的数据资源, 使得以前在计算上不可行的贝叶斯推断技术成为可能, 例如通过蒙特卡罗方法进行大规模采样。
从大量例子中可以发现, 正演预测函数 d=f(θ)随着模型与问题的不同大相径庭, 但对概率性误差部分的描述没有什么变化, 即大部分情况下是一个简单的高斯分布误差:
, (25)
其中, μ和 Σ 分别是高斯分布的均值和协方差矩阵。基于高斯误差的普遍性(根据中心极限定理), 并出于计算过程的简化, 这样的处理是合理的。然而, 概率性的“预测误差”也存在物理信息, 因此同样具备基于具体问题建模的可能性。下面用一个简化的例子进行说明。
若我们希望通过某一时期俯冲带岛弧内某种造岩矿物(如锆石)的数量来构建反演模型, 假设在特定时段内, 俯冲带岛弧中锆石的模型预测产生量 d由俯冲带岩浆活动强度 M 和岛弧岩浆的成分 Φ 这两个参数共同决定, 则确定性的正演模型部分为
。 (26)
若锆石的实际数量为
, 即使对于最基础的数据类型(锆石数量应该是一个自然数, 但高斯误差可以取任意实数), 以高斯误差作为预测值与锆石真值的偏差
(27)
也不合理。事实上, 考虑到锆石产生的物理意义, 泊松分布(而非高斯分布)可能是一个更加合适的建模方法[24]:
(28)
其中, λ为 Poisson 分布的强度参数。
把以上建模方法推广, 可以将贝叶斯模型的似然函数建模分为两个部分: 1) 确定性的正演模型预测部分, d=f(θ); 2) 概率性的误差生成模型, Y～Ω(d= f(θ)), 其中 Ω 是某种分布, 以预测值 d 为参数。例如, 对于上文中的锆石产率泊松模型, 就有
; (29)
对于现有的高斯误差模型, 则有
。 (30)
注意, 相对于之前的描述, 这里使用“误差生成模型”而非“预测误差”, 希望传达的信息是, 误差本身也存在可描述的生成过程, 其中可能包含物理过程, 甚至可能很重要。我们希望指出, 这可能是将贝叶斯方法作为建模方法时, 让模型能够更真实地反映物理世界的可能途径。
贝叶斯方法的一大特点是可以利用先验信息, 这会直接影响后验结果。从式(3)可以发现, 在极端的情形下, 若某处参数取值的先验概率为 0, 则无论数据和似然函数如何, 该处的后验概率也一定是 0。同样, 先验信息越精确, 后验结果也会越理想(在先验信息与数据不存在矛盾的情形下)。当然, 我们希望尽可能提供更精准的先验信息。从这个角度上讲, 贝叶斯反演框架为我们提供了一条利用定量先验信息的途径。例如, Holbach 等[6]将地震层析成像提供的观测信息作为先验信息, 加入地幔动力学反演模型中, 改进了流变学参数的后验分布反演效果。
同时, 我们应该意识到, 越是精确的先验信息, 一旦存在误解, 对模型推测结果的误导会很大。事实上, 在广泛应用贝叶斯统计方法的其他领域存在一种思路(主要是在数据不那么缺乏的情形下), 即尽可能地减少先验信息的指定, 从而达到最客观的模型推测结果。为此, 存在一系列构建所谓“无信息先验”方法, 诸如 Jefferys[93]的无信息先验分布。
总而言之, 可靠性和精确性是此消彼长的。精确的先验信息可以带来结果的精确性, 但一旦先验信息有偏差, 后验分布也会有非常大的偏差。宽泛的先验信息甚至是无信息先验, 在数据量不足的情形下, 则可能带来信息量不够、结果不够唯一的后验分布。先验分布需要在选择时把握好平衡, 也是在进行贝叶斯反演时需要认真考虑的问题。
对某个具体问题, 先验概率究竟在何种程度下影响后验概率同样值得探讨。很多情形下, 先验概率对后验概率的影响并不大, 后验概率往往主要由数据决定。例如, Pease 等[57]的讨论就为大陆地壳地球化学成分的不同先验信息提供了定量评估。
对于先验分布更加详尽的理论探讨, 可参阅文献[20]。
5.3.1 数值计算方法
前面介绍了一些计算贝叶斯后验的方式, 包括半解析的后验分布高斯近似[4–5]、随机模拟方法(全误差筛选[68–70,94]和拒绝–接受采样[57])以及马尔可夫链蒙特卡洛方法(Metropolis-Hasting 算法[44,48,57,59–62,66,75]和 RJMCMC[15,18–19,35,39,41])。事实上, 针对贝叶斯后验的统计计算是一个活跃的领域, 计算数学和统计学界针对该问题的算法也在持续发展中[21–22]。其中, 有许多针对特定情形的表现优异的先进算法已经应用到地球科学领域(集中在地球物理领域), 例如哈密顿蒙特卡洛模拟等基于梯度的采样算法[95–100]以及基于变分优化逼近后验分布的变分推断(varia-tional inference, VI)[101]等, 有些工作还将代码实现与 GPU 等硬件加速技术相结合[102]。尽管算法在持续进步, 在贝叶斯后验的统计计算领域很难得到一个一劳永逸的最优算法, 关键在于基于具体的科学问题寻找和设计算法, 针对变维问题的 RJMCMC算法[28,33–34]就是一个典型的范例。
5.3.2 数值计算问题
主流 MCMC 的收敛性能是个很大的问题。尽管可以从理论上证明, 足够长的时间后, 马氏链一定会收敛, 但没有严格的数学检验方法证明何时达到收敛[14]。目前的检验方法(如自相关函数图、马氏链轨迹图检验法、分位数检验、Heidelberger检验、Gelman-Rubin 检验以及谱分析[14])都只是收敛的必要非充分条件。换句话说, 他们都只是“诊断方法(diagnositic method)”, 而非“证明方法”。另一方面, 哪怕可以做到收敛, 收敛性能和计算效率也是需要关注的问题。我们注意到, 目前地球科学领域的贝叶斯反演推断应用工作中, 研究者不大注重模型的收敛诊断, 本文列举的大量贝叶斯方法相关文献中, 只有少数提及 MCMC 相关收敛诊断[6,75]。
本文从先验概率和似然函数的描述以及后验概率的描述与计算方法等方面介绍了贝叶斯模型反演在地质与地球化学领域的应用。基于概率建模和贝叶斯统计理论, 贝叶斯反演结合先验知识与观测约束对模型参数进行定量推断, 获得模型参数的后验概率分布以及相互之间的关系, 与传统的优化反演方法相比, 可以获得更多关于参数的不确定性描述 信息。
如今, 贝叶斯模型反演逐渐被地质与地球化学领域的研究者了解, 其应用也有望变得更加广泛。然而, 贝叶斯模型反演的进一步发展面临两大瓶颈, 其一是对所研究的科学问题进行贝叶斯概率建模, 从先验信息的选取、正演模型和似然函数的定义等角度进行贝叶斯反演的数学描述; 其二是贝叶斯后验计算中, 贝叶斯反演面临高维采样的低效问题, 在优化硬件资源的同时, 应用或发展更具针对性和更高效的求解算法同样是一大难题。
无论如何, 理论与技术的发展是地球科学发展进步的核心推动力。或许, 加强统计学与地球科学的跨学科合作, 会是地球科学发展的可靠途径。
致谢 感谢北京大学地球与空间科学学院博士研究生宋文周在本文撰写过程中提供的帮助。
参考文献
[1] Gallagher K, Charvin K, Nielsen S, et al. Markov chain Monte Carlo (MCMC) sampling methods to determine optimal models, model resolution and model choice for Earth Science problems. Marine and Petroleum Geolo-gy, 2009, 26(4): 525–535
[2] Turcotte D L, Schubert G. Geodynamics. Third edition. Cambridge: Cambridge University Press, 2014
[3] Karato S. Deformation of earth materials: an introduc-tion to the rheology of solid earth. Cambridge: Cambri-dge University Press, 2008
[4] Hu J, Rudi J, Gurnis M, et al. Constraining Earth’s no-nlinear mantle viscosity using plate-boundary resol-ving global inversions. Proceedings of the National Academy of Sciences, 2024, 121(28): e2318706121
[5] Rudi J, Gurnis M, Stadler G. Simultaneous inference of plate boundary stresses and mantle rheology using adjoints: large-scale 2-D models. Geophysical Journal International, 2022, 231(1): 597–614
[6] Holbach L, Gurnis M, Stadler G. A Bayesian level set method for identifying subsurface geometries and rheo-logical properties in Stokes flow. Geophysical Journal International, 2023, 235(1): 260–272
[7] Ratnaswamy V, Stadler G, Gurnis M. Adjoint-based es-timation of plate coupling in a non-linear mantle flow model: theory and examples. Geophysical Journal In-ternational, 2015, 202(2): 768–786
[8] Tarantola A. Inverse problem theory and methods for model parameter estimation [EB/OL]. (2005) [2024–10–29]. http://epubs.siam.org/doi/book/10.1137/1.978 0898717921
[9] Goldblatt C, Lenton T M, Watson A J. Bistability of atmospheric oxygen and the Great Oxidation. Nature, 2006, 443: 683–686
[10] Gregory B S, Claire M W, Rugheimer S. Photochemi-cal modelling of atmospheric oxygen levels confirms two stable states. Earth and Planetary Science Letters, 2021, 561: 116818
[11] Dhuime B, Hawkesworth C J, Cawood P A, et al. A change in the geodynamics of continental growth 3 billion years ago. Science, 2012, 335: 1334–1336
[12] Rino S, Komiya T, Windley B F, et al. Major episo- dic increases of continental crustal growth determined from zircon ages of river sands: implications for man-tle overturns in the Early Precambrian. Physics of the Earth and Planetary Interiors, 2004, 146(1/2): 369–394
[13] Belousova E A, Kostitsyn Y A, Griffin W L, et al. The growth of the continental crust: constraints from zircon Hf-isotope data. Lithos, 2010, 119(3/4): 457–466
[14] Gelman A, Carlin J B, Stern H S, et al. Bayesian data analysis third edition [EB/OL]. (2021–02–15) [2024–05–21]. https://sites.stat.columbia.edu/gelman/book/B DA3.pdf
[15] Gallagher K, Bodin T, Sambridge M, et al. Inference of abrupt changes in noisy geochemical records using transdimensional changepoint models. Earth and Pla-netary Science Letters, 2011, 311(1/2): 182–194
[16] Spencer C J, Davies N S, Gernon T M, et al. Compo-sition of continental crust altered by the emergence of land plants. Nature Geoscience, 2022, 15(9): 735–740
[17] Gallagher K. Transdimensional inverse thermal history modeling for quantitative thermochronology. Journal of Geophysical Research: Solid Earth, 2012, 117(B2): 2011JB008825
[18] Licciardi A, Gallagher K, Clark S A. A Bayesian app-roach for thermal history reconstruction in basin mo-deling. Journal of Geophysical Research: Solid Earth, 2020, 125(7): e2020JB019384
[19] Gallagher K, Parra M. A new approach to thermal history modelling with detrital low temperature ther-mochronological data. Earth and Planetary Science Letters, 2020, 529: 115872
[20] van de Schoot R, Depaoli S, King R, et al. Bayesian statistics and modelling [J/OL]. Nature Reviews Me-thods Primers. (2021) [2024–05–21]. https://www.na ture.com/articles/s43586-020-00001-2
[21] Robert C, Casella G. A short history of Markov Chain Monte Carlo: subjective recollections from incomplete data. Statistical Science, 2011, 26(1): 102–115
[22] Robert C P, Elvira V, Tawn N, et al. Accelerating MCMC algorithms. WIREs Computational Statistics, 2018, 10(5): e1435
[23] Gelfand A E, Smith A F M. Sampling-based approaches to calculating marginal densities. Journal of the Ameri-can Statistical Association, 1990, 85: 398–409
[24] Ross S M. Introduction to probability theory // Introdu-ction to probability models (twelfth edition) [M/OL]. (2019) [2022–04–30]. https://linkinghub.elsevier.com/ retrieve/pii/B9780128143469000068
[25] Metropolis N, Rosenbluth A W, Rosenbluth M N, et al. Equation of State Calculations by Fast Computing Ma-chines. The Journal of Chemical Physics, 1953, 21(6): 1087–1092
[26] Hastings W K. Monte Carlo sampling methods using Markov Chains and their applications. Biometrika, 1970, 57(1): 97–109
[27] Geman S, Geman D. Stochastic relaxation, Gibbs di-stributions, and the Bayesian restoration of images. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1984, PAMI-6: 721–741
[28] Green P J. Reversible jump Markov chain Monte Carlo computation and Bayesian model determination. Bio-metrika, 1995, 82(4): 711–732
[29] Hoffman M D, Gelman A. The No-U-turn sampler: adaptively setting path lengths in hamiltonian Monte Carlo. Journal of Machine Learning Research, 2014, 15: 1593–1623
[30] Marinari E, Parisi G. Simulated tempering: a new Monte Carlo scheme. Europhysics Letters, 1992, 19 (6): 451–458
[31] Gelfand A E, Smith A F M. Sampling-based approaches to calculating marginal densities. Journal of the Ame-rican Statistical Association, 1990, 85: 398–409
[32] Sambridge M, Bodin T, Gallagher K, et al. Trans-dimensional inference in the geosciences. Philosophi-cal Transactions of the Royal Society A: Mathemati-cal, Physical and Engineering Sciences, 2013, 371: 20110547
[33] Green P. A Primer on Markov Chain Monte Carlo // Barndorff-Nielsen O, Cox D, Klüppelberg C. Complex stochastic systems, vol 87 [M/OL]. (2000) [2024–10–28]. http://www.crcnetbase.com/doi/abs/10.1201/9781 420035988.ch1
[34] Green P J. Trans-dimensional Markov Chain Monte Carlo // Highly Structured Stochastic Systems: Vol 27. Oxford: Oxford University Press, 2003: 179–198
[35] Malinverno A, Briggs V A. Expanded uncertainty quantification in inverse problems: hierarchical Bayes and empirical Bayes. Geophysics, 2004, 69(4): 1005–1016
[36] Livermore P W, Fournier A, Gallet Y, et al. Transdi-mensional inference of archeomagnetic intensity chan-ge. Geophysical Journal International, 2018, 215(3): 2008–2034
[37] Ruggieri E. A Bayesian approach to detecting change points in climatic records: the Bayesian change point algorithm. International Journal of Climatology, 2013, 33(2): 520–528
[38] Sambridge M. Reconstructing time series and their uncertainty from observations with universal noise. Journal of Geophysical Research: Solid Earth, 2016, 121(7): 4990–5012
[39] Jasra A, Stephens D A, Gallagher K, et al. Bayesian mixture modelling in geochronology via Markov Chain Monte Carlo. Mathematical Geology, 2006, 38(3): 269–300
[40] Sambridge M, Gallagher K, Jackson A, et al. Trans-dimensional inverse problems, model comparison and the evidence. Geophysical Journal International, 2006, 167(2): 528–542
[41] Hopcroft P O, Gallagher K, Pain C C. Inference of past climate from borehole temperature data using Bayesian Reversible Jump Markov chain Monte Carlo: Bayesian inversion of borehole data for past climate. Geophy-sical Journal International, 2007, 171(3): 1430–1439
[42] Allègre C J, Rousseau D. The growth of the continent through geological time studied by Nd isotope analysis of shales. Earth and Planetary Science Letters, 1984, 67(1): 19–34
[43] Zindler A, Hart S. Chemical geodynamics. Annual Re-view of Earth and Planetary Sciences, 1986, 14(1): 493–571
[44] Liu B, Liang Y. An introduction of Markov chain Monte Carlo method to geochemical inverse problems: reading melting parameters from REE abundances in abyssal peridotites. Geochimica et Cosmochimica Acta, 2017, 203: 216–234
[45] Philpotts A R, Ague J J. Principles of igneous and me-tamorphic petrology. Third edition. Cambridge: Cam-bridge University Press, 2022
[46] Shaw D M. Trace elements in magmas: a theoretical treatment. Cambridge: Cambridge University Press, 2006
[47] Liang Y, Liu B. Simple models for disequilibrium frac-tional melting and batch melting with application to REE fractionation in abyssal peridotites. Geochimica et Cosmochimica Acta, 2016, 173: 181–197
[48] Brown E L, Petersen K D, Lesher C E. Markov chain Monte Carlo inversion of mantle temperature and sour-ce composition, with application to Reykjanes Penin-sula, Iceland. Earth and Planetary Science Letters, 2020, 532: 116007
[49] Torrance K E, Turcotte D L. Thermal convection with large viscosity variations. Journal of Fluid Mechanics, 1971, 47(1): 113–125
[50] Mckenzie D P, Roberts J M, Weiss N O. Convection in the earth’s mantle: towards a numerical simulation. Journal of Fluid Mechanics, 1974, 62(3): 464–538
[51] Zhong S, McNamara A, Tan E, et al. A benchmark study on mantle convection in a 3-D spherical shell using CitcomS: benchmarks of 3-D spherical convec-tion models [J/OL]. Geochemistry, Geophysics, Geo-systems. (2008) [2022–03–28]. http://doi.wiley.com/ 10.1029/2008GC002048
[52] Rudnick R L, Gao S. Composition of the continental crust // Holland H D, Turekian K K. Treatise on Geo-chemistry (Second Edition). Oxford: Newnes (Elsevi-er), 2014: 1–51
[53] Cawood P A, Hawkesworth C J, Dhuime B. The contin-ental record and the generation of continental crust. Geological Society of America Bulletin, 2013, 125 (1/2): 14–32
[54] Hawkesworth C J, Cawood P A, Dhuime B, et al. Earth’s continental lithosphere through time. Annual Review of Earth and Planetary Sciences, 2017, 45(1): 169–198
[55] Korenaga J. Crustal evolution and mantle dynamics through Earth history. Philosophical Transactions of the Royal Society A: Mathematical, Physical and En-gineering Sciences, 2018, 376: 20170408
[56] Hawkesworth C, Cawood P A, Dhuime B, et al. Tec-tonic processes and the evolution of the continental crust. Journal of the Geological Society, 2024, 181(4): jgs2024-027
[57] Pease G, Gelb A, Lee Y, et al. A Bayesian formulation for estimating the composition of Earth’s crust. Journal of Geophysical Research: Solid Earth, 2023, 128(7): e2023JB026353
[58] Connolly J A D. Computation of phase equilibria by linear programming: a tool for geodynamic modeling and its application to subduction zone decarbonation. Earth and Planetary Science Letters, 2005, 236(1/2): 524–541
[59] Zhang T, Keller C B, Hoggard M J, et al. A Bayesian framework for subsidence modeling in sedimentary basins: a case study of the Tonian Akademikerbreen group of Svalbard, Norway. Earth and Planetary Science Letters, 2023, 620: 118317
[60] McDannell K T, Keller C B. Cryogenian glacial erosion of the central Canadian Shield: the “late” great uncon-formity on thin ice. Geology, 2022, 50(12): 1336–1340
[61] Cox A A, Keller C B. A Bayesian inversion for emissions and export productivity across the end-Cretaceous boundary. Science, 2023, 381: 1446–1451
[62] Edwards G H, Keller C B, Newton E R, et al. An early giant planet instability recorded in asteroidal meteo-rites [J/OL]. Nature Astronomy. (2024) [2024–10–03]. https://www.nature.com/articles/s41550-024-02340-6
[63] Greber N D, Dauphas N, Bekker A, et al. Titanium isotopic evidence for felsic crust and plate tectonics 3. 5 billion years ago. Science, 2017, 357: 1271–1274
[64] Greber N D, Dauphas N. The chemistry of fine-grained terrigenous sediments reveals a chemically evolved Pa-leoarchean emerged crust. Geochimica et Cosmochi-mica Acta, 2019, 255: 247–264
[65] Tang M, Chen K, Rudnick R L. Archean upper crust transition from mafic to felsic marks the onset of plate tectonics. Science, 2016, 351: 372–375
[66] Ptáček M P, Dauphas N, Greber N D. Chemical evo-lution of the continental crust from a data-driven in-version of terrigenous sediment compositions. Earth and Planetary Science Letters, 2020, 539: 116090
[67] Chen K, Rudnick R L, Wang Z, et al. How mafic was the Archean upper continental crust? Insights from Cu and Ag in ancient glacial diamictites. Geochimica et Cosmochimica Acta, 2020, 278: 16–29
[68] Rosas J C, Korenaga J. Rapid crustal growth and efficient crustal recycling in the early Earth: impli-cations for Hadean and Archean geodynamics. Earth and Planetary Science Letters, 2018, 494: 42–49
[69] Guo M, Korenaga J. Argon constraints on the early growth of felsic continental crust. Science Advances, 2020, 6(21): eaaz6234
[70] Guo M, Korenaga J. The combined Hf and Nd isotope evolution of the depleted mantle requires Hadean con-tinental formation. Science Advances, 2023, 9(12): eade2711
[71] Butler R F. Paleomagnetism: magnetic domains to geologic terranes [EB/OL]. (2004) [2022–04–30]. https://www.geo.arizona.edu/Paleomag/tocpref.pdf
[72] Tauxe L. Essentials of paleomagnetism. Berkeley: Uni-versity of California Press, 2010
[73] Vine F J, Matthews D H. Magnetic anomalies over oceanic ridges. Nature, 1963, 199: 947–949
[74] Heslop D, Roberts A P, Oda H, et al. An automatic model selection‐based machine learning framework to estimate FORC distributions. Journal of Geophysical Research: Solid Earth, 2020, 125(10): e2020JB020418
[75] Cych B, Morzfeld M, Tauxe L. Bias corrected estima-tion of paleointensity (BiCEP): an improved methodo-logy for obtaining paleointensity estimates. Geochemi-stry, Geophysics, Geosystems, 2021, 22(8): e2021GC 009755
[76] Schnepp E, Obenaus M, Lanos P. Posterior archaeoma-gnetic dating: an example from the Early Medieval site Thunau am Kamp, Austria. Journal of Archaeological Science: Reports, 2015, 2: 688–698
[77] Hellio G, Gillet N, Bouligand C, et al. Stochastic mo-delling of regional archaeomagnetic series. Geophy-sical Journal International, 2014, 199(2): 931–943
[78] Lanos P. Bayesian Inference of calibration curves: application to archaeomagnetism // Buck C E, Millard A R. Tools for Constructing Chronologies: vol 177 [M/OL]. (2004) [2024–10–29]. http://link.springer. com/10.1007/978-1-4471-0231-1_3
[79] Nagata T, Arai Y, Momose K. Secular variation of the geomagnetic total force during the last 5000 years. Journal of Geophysical Research, 1963, 68(18): 5277–5281
[80] Paterson G A. A simple test for the presence of mul-tidomain behavior during paleointensity experiments. Journal of Geophysical Research, 2011, 116(B10): B10104
[81] Tauxe L, Santos C N, Cych B, et al. Understanding nonideal paleointensity recording in igneous rocks: insights from aging experiments on lava samples and the causes and consequences of “fragile” curvature in arai plots. Geochemistry, Geophysics, Geosystems, 2021, 22(1): e2020GC009423
[82] Wang J, Zuo R. Uncertainty quantification in geoche-mical mapping: a review and recommendations. Geo-chemistry, Geophysics, Geosystems, 2024, 25(3): e2023GC011301
[83] Arthur G. International geochemical mapping: a new global project. Journal of Geochemical Exploration, 1990, 39(1/2): 1–13
[84] Reimann C, Melezhik V, Niskavaara H. Low-density regional geochemical mapping of gold and palladium highlighting the exploration potential of northernmost Europe. Economic Geology, 2007, 102(2): 327–334
[85] Batanova V G, Sobolev A V. Compositional heteroge-neity in subduction-related mantle peridotites, troodos massif, cyprus. Geology, 2000, 28(1): 55–58
[86] Pearce J A, Stern R J, Bloomer S H, et al. Geochemical mapping of the Mariana arc-basin system: implications for the nature and distribution of subduction compo-nents. Geochemistry, Geophysics, Geosystems, 2005, 6(7): Q07006
[87] Huang D, Zuo R, Wang J. Geochemical anomaly identi-fication and uncertainty quantification using a Baye-sian convolutional neural network model. Applied Geochemistry, 2022, 146: 105450
[88] Younesi A, Ansari M, Fazli M, et al. A comprehensive survey of convolutions in deep learning: applications, challenges, and future trends. IEEE Access, 2024, 12: 41180–41218
[89] Kaur G, Oberai Er N. A review article on Naive Bayes classifier with various smoothing techniques. Interna-tional Journal of Computer Science and Mobile Com-puting, 2014, 3(10): 864–868
[90] Bharadiya J. A Review of Bayesian machine learning principles, methods, and applications. International Journal of Innovative Research in Science Engineering and Technology, 2023, 8: 2033–2038
[91] Jospin L V, Laga H, Boussaid F, et al. Hands-on Baye-sian neural networks — a tutorial for deep learning users. IEEE Computational Intelligence Magazine, 2022, 17(2): 29–48
[92] Girin L, Leglaive S, Bie X, et al. Dynamical variational autoencoders: a comprehensive review. Foundations and Trends in Machine Learning, 2021, 15(1/2): 1–175
[93] Jeffreys H. The theory of probability. Oxford: Oxford University Press, 1998
[94] Korenaga J. Estimating the formation age distribution of continental crust by unmixing zircon ages. Earth and Planetary Science Letters, 2018, 482: 388–395
[95] Gebraad L, Boehm C, Fichtner A. Bayesian elastic full-waveform inversion using Hamiltonian Monte Carlo. Journal of Geophysical Research: Solid Earth, 2020, 125(3): e2019JB018428
[96] Fichtner A, Zunino A, Gebraad L. Hamiltonian Monte Carlo solution of tomographic inverse problems. Geo-physical Journal International, 2019, 216(2): 1344–1363
[97] Zunino A, Gebraad L, Ghirotto A, et al. HMCLab: a framework for solving diverse geophysical inver- se problems using the Hamiltonian Monte Carlo me-thod. Geophysical Journal International, 2023, 235(3): 2979–2991
[98] Biswas R, Sen M K. Transdimensional 2D Full-wave-form inversion and uncertainty estimation [EB/OL]. (2022) [2024–10–29]. http://arxiv.org/abs/2201.09334
[99] Sen M K, Biswas R. Transdimensional seismic inver-sion using the reversible jump Hamiltonian Monte Carlo algorithm. Geophysics, 2017, 82(3): R119–R134
[100] Peng R, Han B, Hu X, et al. 2-D probabilistic inversion of MT data and uncertainty quantification using the Hamiltonian Monte Carlo method. Geophysical Jour-nal International, 2024, 237(3): 1826–1841
[101] Zhao X, Curtis A. Bayesian inversion, uncertainty analysis and interrogation using boosting variational inference. Journal of Geophysical Research: Solid Earth, 2024, 129(1): e2023JB027789
[102] Mandolesi E, Ogaya X, Campanyà J, et al. A rever-sible-jump Markov Chain Monte Carlo algorithm for 1D inversion of magnetotelluric data. Computers & Geosciences, 2018, 113: 94–105
Application of Bayesian Method in Geology and Geochemistry
Abstract The introduction to Bayesian methods is presented from the perspective of inversion modeling, outlining the fundamental principles, computational techniques, development history, and current research status. Various methods and their applications across different fields of geology and geochemistry are analysed in detail. Finally, the future development directions and challenges are discussed.
Key words Bayesian method; model inversion; stochastic simulation; prior probability; posterior probability; likelihood
国家自然科学基金(92155204)资助
收稿日期: 2024–11–11;
修回日期: 2025–01–15