Home Journals Progress in Geophysics
Progress in Geophysics

Abbreviation (ISO4): Prog Geophy      Editor in chief:

About  /  Aim & scope  /  Editorial board  /  Indexed  /  Contact  / 

Bayesian inversion: an important method and challenge for seismic source parameter inversion research

  • Zhe WANG ,
  • YunHua LIU
Expand

Received date: 2023-10-17

  Online published: 2024-12-19

Copyright

Copyright ©2024 Progress in Geophysics. All rights reserved.

Abstract

The role of geophysical inversion in seismic research and prediction is of paramount importance. This paper seeks to comprehensively outline the constraints associated with conventional inversion techniques, focusing on the introduction of a Bayesian-based uncertainty inversion method. Bayesian inversion involves computing posterior distributions utilizing diverse prior distributions and likelihood functions, with special emphasis on established techniques like the Markov Chain Monte Carlo (MCMC) and variational inference methods, thereby augmenting the reliability of inversion outcomes.The manuscript furnishes an elaborate exposition on pivotal techniques within Bayesian inversion, notably delving into regularization methods (such as Laplace and von Karman regularization) that confine the parameter space in seismic inversion, validated through rigorous case studies. Moreover, it expounds on sampling methodologies (including the Metropolis-Hastings algorithm and Gibbs sampling) that facilitate parameter space sampling and approximate posterior distributions. The application of the Metropolis-Hastings algorithm in seismic inversion is meticulously elucidated.The discussion accentuates the criticality of model parameter selection, notably the influence of uncertainty associated with fault geometric shape selection on inversion results. Additionally, it probes into the challenges encountered in constructing finite fault source models and presents a Bayesian-based case study evaluating the credibility of different slip model clusters.In conclusion, the paper summarizes the limitations inherent in the Bayesian approach and delineates potential avenues for future research directions. In the realm of geophysical inversion, the application of Bayesian methods presents novel prospects for overcoming the constraints of traditional methodologies.

Cite this article

Zhe WANG , YunHua LIU . Bayesian inversion: an important method and challenge for seismic source parameter inversion research[J]. Progress in Geophysics, 2024 , 39(5) : 1771 -1787 . DOI: 10.6038/pg2024HH0340

0 引言

频繁发生的地震对人类的生命安全和财产造成了巨大的威胁,特别是特大地震造成的伤亡人数居高不下.其中几次破坏力极强的特大地震包括震级MW9.5的智利大地震(Barrientos and Ward, 1990)、震级MW9.1的苏门答腊岛地震(Ammon et al., 2005),以及我国发生在2008年震级MW8.0的汶川地震(谢世亮和孙玉军,2021).这些特大地震导致了数以万计的死亡人数和无法估量的财产损失.地球物理反演参数在地震灾害研究中被广泛应用,例如用于地震震源机制研究(郑勇等,2009)、地震波传播模拟(裴正林和牟永光,2004)和地震灾害损失评估(王晓青等,2009)等方面.通过地球物理反演震源参数,可以更准确地确定地震的震源位置、深度和震级等信息.通过研究地震的震源位置、深度和震级等参数,可以更好地了解地震的强度和影响范围,进而评估可能受到影响的地区的损失情况.因此,地球物理反演震源参数对地震灾害研究具有重要意义.
地壳结构的复杂性使我们难以直接通过观测来获取断层面的错动情况.目前已经有多种大地测量技术应用于地震后地表形变的测量,但我们更希望通过测量得到的形变数据通过地球物理反演进一步研究发震机制.反演技术通过搜索正演模型,然后计算正演模型,找到一组模拟形变值与真实形变值最匹配的一组参数作为反演结果.正反演问题在地球物理学中密切相关,反演问题建立在正演问题的基础上(刘博瀚,2021).位错模型理论的发展为利用大地测量数据反演断层参数提供了基础.许多学者使用基于均匀弹性半空间的位错模型进行反演.例如,árnadóttir等(2001)利用大地测量的同震形变数据,基于均匀弹性半空间的位错模型采用非线性拟牛顿法反演确定断层的几何参数和滑动分布.王文萍和王庆良(1999)利用垂直形变数据,基于均匀弹性半空间的位错模型采用遗传算法反演断层滑动的几何参数,然后利用拟牛顿阻尼最小二乘法求解非均匀滑动量.Monelli等(2009)通过拟合强运动和GPS数据运用贝叶斯反演方法对2000年鸟取西部地震的破裂过程进行成像.Fukuda和Johnson(2008)对中国台湾集集地震采用了贝叶斯框架下的全概率反演方法,相较于传统的反演方法,贝叶斯算法可以同时估计滑动分布和平滑参数,且反演得到的结果比传统ABIC方法更加光滑,结果也很贴近传统的交叉验证法.
地球物理反演具有复杂性和不确定性,地球的介质、形状和断层几何复杂性,包括形变观测数据的误差都会使建立的函数模型存在一定的误差(系统误差、粗差等)(韩月萍等,2015).因此针对不同地震,我们需要建立更加符合实际的模型,采用更加可信的反演方法,使研究结果更加可靠,效率更高.下文将首先介绍传统的反演方法及其局限性,进而介绍贝叶斯算法的基本理论及其优缺点、应用实例和未来的研究方向,以期能够为地震形变反演工作提供参考.

1 地球物理反演过程中的传统方法

地球物理反演问题是指通过对地球物理场进行测量,推断出地下介质的性质和结构的过程.由于反演问题通常是非线性且不可逆的,因此需要运用各种数学方法来解决.地球物理反演利用地表观测结果和其他地学观测数据,结合现有的物理模型,计算相关物理参数,修正或提出更精确的地球运动学模型(王乐洋,2012).研究地震断层形变可以通过地表观测结果来探索地壳运动与地震的关系.大地测量反演已经从单一反演逐渐发展为多类观测数据的联合反演,并且反演模型也在不断创新,从连续模型到离散模型,从线性反演到非线性反演.大地测量反演震源参数主要使用解析法和数值法,解析法适用于线性和简单非线性问题,而数值模拟法适用于大部分震源参数模拟问题(张冬丽等,2004).反演断层的几何参数通常采用非线性反演方法,如模拟退火法、遗传算法、神经网络法和粒子群算法.而断层的运动参数通常使用附加约束的非负约束最小二乘法和总体最小二乘法进行求解.InSAR震源参数反演的精度受到多个因素的影响,包括断层几何、断层面离散化方法、平滑因子的选择、位错模型的选择、约束数据的多样性、数据权重的确定方法和降采样方法等.

1.1 粒子群算法

粒子群反演算法是地球物理学中一种重要的数值模拟方法,利用粒子群优化算法对目标函数进行优化,以获得最优解(罗腾等,2019).粒子群反演算法(Particle Swarm Optimization,PSO)是一种基于群体智能优化算法的反演方法.PSO算法最早可以追溯到20世纪90年代,Kennedy和Eberhart (1995)基于对鸟群、鱼群的模拟提出粒子群优化算法.最初,PSO是为了解决连续优化问题,Kennedy和Eberhart (1995)又提出了PSO的离散二进制版用来解决工程实际中的组合优化问题.近年来,粒子群算法在地球物理参数反演研究中得到广泛应用,其收敛效率高,所需的控制参数较少,而且所求得的最优解相对稳定可靠.易远元和王家映(2009)将粒子群算法引入到地球物理学非线性反演的问题中.该算法在模型空间中随机进行全局搜索,对初始值不敏感,因此不易陷入局部最优解,且该算法还具备快速收敛的特点,但可能出现收敛过早的问题.针对基本的PSO算法可能存在的收敛过早现象,特别是在处理多峰值问题时更容易陷入局部最优解,冯万鹏和李振洪(2010)对局部PSO算法进行了进一步改进,将粒子在空间范围内的密度函数作为判别最优值的指标,结合了单纯型(Simplex)算法的高效收敛优势,提出了MPSO(Multipeaks Particle Swarm Optimization)算法.然而,目前粒子群算法在地球物理反演中仍存在一些局限性,包括对非线性问题的反演效果不佳(靳锡波,2022)、对复杂介质结构的反演困难(阳佳慧和刘庆华,2020)、对参数的依赖性较强以及随机性对反演结果的影响较大(戴前伟等,2019).针对这些局限性,需要根据具体问题的特点选择适当的反演方法,或者结合其他算法进行改进,以提高反演结果的精度和可靠性.

1.2 模拟退火法

模拟退火法是一种基于概率的全局优化方法,最早可以追溯到20世纪50年代,Metropolis基于固体的退火过程提出模拟退火法的思想.直到1983年,Kirkpatrick等人才成功将Metropolis准则引入到优化过程中,并得到一种对Metropolis算法进行迭代的优化算法,称之为“模拟退火法”(Kirkpatrick et al., 1983).Rothman(1986)首次将首次将模拟退火法引入地球物理方法中,并解决了地震勘探中的剩余静矫校正问题.此后,在解决局部极小问题上模拟退火法表现良好,引起众多学者关注,从而被广泛应用于地球物理学反演.其中经典的模拟退火法被称为Metropolis模拟退火法(胡桂武等,2005).为了提高算法的效率,人们还提出了其他改进的模拟退火算法,其中较为典型的是热浴(Heat bath)模拟退火法(简称HBSA算法)(Lu et al., 2016)和极快速模拟退火算法(VFSA算法)(Chunduru et al., 1997).同时,国内外研究人员不断尝试从不同的阶段入手对算法进行改进,并将模拟退火法同其他计算智能方法相结合应用到各种复杂的系统建模和优化问题中(张繁昌等,1997; 杨立强等,2005).与传统的反演方法相比,模拟退火法具有更强的全局搜索能力,从而能寻找全局最优解.
然而国内对模拟退火算法的研究历史较短,相比蚁群算法等其他算法,研究相对较少.主要集中于模拟退火在工程等实际应用方面的效果和效率,而对其理论性研究不够深入.尽管模拟退火算法在全局极小值搜索方面具有优势,但也存在收敛速度较慢的缺点(陈晓等,2016).

1.3 最小二乘法

最小二乘法最初用于计算谷神星轨道(Paris, 2012),随后在多个领域得到了广泛应用,特别是在测量领域.作为一种基本的平差方法,最小二乘法通过最小化观测值与模拟值之间的残差平方和来确定模型参数(张彦栋,2018).Legendre(1806)首次简明地阐述了高斯和勒让德对最小二乘法的贡献,高斯-马尔科夫(Gauss-Markov)模型假设函数模型是已知的非随机的,只在观测值向量中含有随机误差.在地球物理学领域,何宗海(1994)在地震学中采用了阻尼最小二乘法算法来计算地震的震源加速度.李志才等(2009)基于GPS观测数据反演了汶川地震断层参数,采用弹性半空间均匀位错理论,并采用有界变量最小二乘法反演了汶川地震断层滑动分布.温扬茂等(2012)利用自动剖分技术获得了断层面的非均匀子块,并通过最小二乘法反演了断层的精细滑动分布.目前,最小二乘法也被广泛应用于大地测量反演中的线性反演问题.然而,现有的滑动分布反演研究多是基于最小二乘法或约束最小二乘准则.如果将总体最小二乘方法应用于滑动分布反演,将会得到更合理的滑动参数,更好地拟合地表形变位移.
尽管最小二乘法仅考虑了观测向量的误差,忽略了系数矩阵的误差,但总体最小二乘法在反演滑动分布时能够同时考虑观测向量和系数矩阵(格林函数矩阵)所包含的误差.目前,国内外许多研究人员对病态总体最小二乘法进行了研究,主要方法包括Tikhonov正则化法(梁勇等,2019)、广义正则化法(杨冠雨等,2020)、岭估计法(Dorugade, 2014)、虚拟观测法(Michalos et al., 2018)等.
尽管最小二乘法被广泛使用,但它也存在一些局限性.最小二乘法的反演结果在很大程度上受到先验信息的影响,需要提前获得一定的先验信息才能进行反演.此外,对数据质量有较高的要求,低质量数据或异常点可能导致不可靠的反演结果.对反演模型参数的数量和质量也有一定要求,需要事先选择适当的参数,并对参数进行约束和平滑处理,以获得可靠的反演结果.为了克服最小二乘法的局限性,学者们提出了许多改进方法,例如在最小二乘法基础上加入正则化项,通过调节正则化系数来平衡模型拟合精度和复杂度之间的关系(周绍民等,2017).正则化最小二乘法可以有效降低过拟合风险,并提高反演结果的可靠性.将贝叶斯方法与最小二乘法结合,可以通过考虑先验信息来减小反演结果的不确定性,提高反演结果的可靠性(王芳和陈生昌,1999).贝叶斯最小二乘法可以有效处理反演模型的不确定性问题.

1.4 遗传算法

遗传算法是一种基于生物进化机制衍生而来的搜索优化算法,通过模拟生物进化过程来寻找最优解.作为非线性反演方法家族的核心成员之一,它在非线性反演的研究和发展中扮演着重要角色.其基本思路是利用自然选择和基因遗传机制,不断演化种群以逐步接近最优解(马随波,2020).在地球物理学反演中,遗传算法被广泛应用于解决非线性和非凸优化问题,具有较高的可靠性和灵活性.遗传算法最早由John Holland于1975年提出,与MPSO算法和模拟退火法类似,它也是一种非线性的启发式全局最优算法.目前,遗传算法已经在多个领域的研究中得到应用.例如,王文萍和王庆良(1999)将遗传算法与普通最小二乘法结合,成功反演了1990年青海共和MS7.0级地震的震源参数;许才军和温扬茂(2008)利用InSAR数据获得了1997年西藏玛尼地震的同震形变场,然后运用遗传算法反演得到断层的均匀滑动参数.师学明和王家映(2008)通过在地球物理资料反演中应用遗传算法的实例,系统总结和归纳了遗传算法的特点及其局限性.周江(2018)应用遗传算法研究了震源机制反演问题,并反演得到了微震的震源机制解和震源参数.
与模拟退火反演算法相比,遗传算法在地球物理反演中具有一些优点:它不依赖于初始模型的选择,能够找到全局最小点而不陷入局部极小,并且在反演过程中无需计算雅克比偏导数矩阵等.此外,遗传算法具有隐性并行性,是一种群智能演化算法.随着并行计算和集群式计算机技术的发展,这种算法将会得到越来越广泛的研究和应用.
然而,遗传算法在地球物理反演中存在一些挑战.它具有全局搜索能力但容易陷入局部最优解,需采取多种策略避免此问题.处理大规模反演需大量计算资源和时间.算法性能受参数选择和调整影响,如交叉率、变异率.遗传算法得出的最优解通常是一组参数,缺乏明确物理解释,限制了可解释性和可视化性.因此,在实际应用中需根据问题选择合适算法、优化参数,并引入多种搜索策略和启发式方法,以提高搜索效果和解释性(王涵,2022).
综上所述可看出,传统的地球物理学反演方法存在多种局限性,例如容易受初始种群影响而出现较强的随机性、容易陷入局部最优解、无法处理非线性问题等.参数反演问题通常涉及多参数、高维和非线性特征,这些特点对传统方法的反演效率和准确性构成挑战.此外,这些方法通常依赖于初始值或先验信息进行反演,因此反演结果容易受初始值或先验信息的影响.针对这些问题,贝叶斯反演方法被提出并逐渐成为地球物理学反演中的重要方法.贝叶斯反演方法能够利用已知的先验信息,并通过与实验数据的比较不断更新先验信息,从而得出更可靠的反演结果.与传统的反演方法相比,贝叶斯反演方法具有考虑先验信息的影响、通过不断更新提高反演结果准确性以及提供反演结果的不确定性信息等优点.

2 贝叶斯反演方法

2.1 贝叶斯方法基本理论

Bayes和Price(1997)首次提出了贝叶斯方法,随后经过Laplace(2012)的进一步推导,发展成一种基于概率统计推断的方法.贝叶斯反演利用先验信息和观测数据构成的似然函数作为约束条件,结合正演计算和适当的采样方法,推测模型参数的后验概率分布,并从中提取统计参数用于模型的预测和评价.
贝叶斯公式可表示为(Tarantola,2005):
式中,m表示模型参数,d表示观测数据.p(m)是模型m的先验概率分布,由先验信息确定.p(d|m)是似然函数,表示给定模型m的情况下获得观测数据d的条件概率大小.在贝叶斯反演中,似然函数相似于确定性反演方法中的目标函数.p(d)被视为归一化常数,使得概率值在整个模型区间内积分大小的总和为1,表示观测数据d的出现概率.p(md)为后验概率,表示给定观测数据d条件下,模型m的条件概率.它表征了在先验信息和观测数据的共同约束下,模型参数m的概率大小.

2.2 贝叶斯方法基本流程

贝叶斯方法作为一种常用的反演方法,可以通过考虑不确定性和先验信息来提高反演的准确性.在贝叶斯反演的第一步中,需要建立一个反演模型,包括待反演的参数、观测数据以及误差模型等.这些模型可以使用数学公式或计算机程序进行表示.第二步是确定先验信息,在进行反演之前,需要确定待反演参数的先验概率分布.这可以通过利用以往的实验数据、已知的地质信息和经验知识等进行获取.先验信息可以以概率密度函数的形式表示.在获得先验信息后,可以使用贝叶斯公式计算后验概率分布,即在给定观测数据后待反演参数的概率分布.具体而言,后验概率分布可以表示为:
其中,p(θ|D)表示给定观测数据D后参数θ的后验概率分布;p(θ|D)表示在参数θ下观测数据D的概率分布,即似然函数;p(θ)表示参数θ的先验概率分布.然后,通过数值方法可以对后验概率分布进行采样,从而得到待反演参数的后验分布.最后,可以对反演结果进行分析和解释,比较不同模型的优劣,以及评估参数的不确定性.具体流程如图 1所示.
图1 贝叶斯方法基本流程图

Figure 1 Basic flowchart of Bayesian method

2.2.1 选择先验分布

先验信息的分布直接影响模型的后验分布,先验信息越丰富,则反演的模型参数就能更精确,从而提高反演的稳定性.针对不同的先验信息,可以采用不同的分布形式.在先验信息较少的情况下,通常采用均匀分布(Bodin et al., 2009);而对于具有一定信息性先验的情况,先验分布更集中于模型参数期望值附近的概率更高,这时候高斯分布更适用.与均匀分布相比,高斯分布能够提供待反演参数的均值和方差,并具有平滑约束的特性.高斯分布常用于地震层位速度的反演(Malinverno and Briggs, 2004)、重力密度测量(冯娟等,2014)和岩石物性反演(王志伟等,2021)等.然而,对于具有离散特征的模型参数,柯西分布比高斯分布更适用.柯西分布的值更多分布在零值附近,而非零值附近的偏离程度比高斯分布更自由,因此柯西分布能够获得含有更多零值的稀疏反演结果(Alemie and Sacchi, 2011).柯西分布常用于获得地层界面信息(林恬,2017)和标记异常物体的研究(Yin and Zhang, 2014).除了高斯分布和柯西分布,还有其他分布函数在不同的研究中得到应用.Theune等(2010)发现拉普拉斯分布也可用于反演层位界面,Visser等(2019)发现Gamma分布可用于反演各层速度值,Dirichlet分布可用于表征层位厚度,Poisson分布可用于表征层位个数.有时也会使用多个分布函数来构建先验分布,这样可以使已知模型的先验信息更加精确,如多维τ分布和Huber分布.多维τ分布是多维高斯分布和多维柯西分布的广义形式,并可通过自由度的转换变为高斯分布和柯西分布.τ分布在表征弹性参数的先验分布上表现较好.Huber分布是高斯分布和拉普拉斯分布的混合分布,在存在噪声的情况下仍能保持反演的稳定性,并在小标准差处具有平滑特性(Guitton and Symes, 2003).当模型参数的性质只与相邻参数相关,与其他参数无关时,可以运用马尔科夫随机场构建先验信息来约束模型参数.马尔科夫随机场能够保证复杂空间模型参数的连续性,常用于储层成像(Avseth et al., 2001)、岩相反演(de Figueiredo et al., 2019)等领域的研究.
这些不同的分布函数和方法的应用能够使先验分布更加符合实际情况,并提供更准确的反演结果.因此,在选择先验分布时需要根据具体问题的特点和已有的先验知识进行合理的选择.

2.2.2 构建似然函数

贝叶斯方法是一种利用概率论和统计学原理进行推断的方法,它的基本思想是在已知先验概率的基础上,通过观测数据来更新我们对未知参数的估计.在贝叶斯方法中,似然函数起着关键作用,它反映了观测数据对参数的约束程度.
在构建似然函数时,通常假设观测数据来自于一个已知的概率分布,而未知参数的取值则是我们要推断的对象.在这种情况下,需要使用已知概率分布和观测数据,通过条件概率的计算来构建似然函数.具体来说,假设有一个参数为θ的概率分布f(θ),同时观测到了一组数据D={x1, x2, …, xn}.那么,给定参数θ的条件下,数据D的概率可以表示为(姚铭等,2020):
其中,f(xi|θ)表示参数为θ时,xi的概率密度函数.这个式子就是似然函数.可以看到,似然函数的值取决于参数θ和观测数据D,因此它可以作为贝叶斯方法中后验概率的核心.需要注意的是,似然函数与后验概率密度函数是不同的概念.似然函数描述的是给定参数的情况下,数据出现的可能性;而后验概率密度函数则是给定数据的情况下,参数的可能取值.通过似然函数和先验概率来计算后验概率密度函数,可以实现对未知参数的推断.
似然函数通常用于描述在给定模型参数的情况下,观测数据出现的概率.在地球物理学反演问题中,似然函数的形式可以采用高斯似然函数、Poisson似然函数、残差平方和似然函数等.例如,在一些情况下,地球物理学反演问题可以看作线性回归问题,因此适合使用高斯似然函数来描述.具体而言,假设观测数据为y,模型预测为,观测误差服从正态分布,那么高斯似然函数可以表示为(Bagnardi and Hooper, 2018):
其中θ表示模型参数,σ表示观测误差的标准差.
在某些地球物理学反演问题中,观测数据可以被视为计数数据,例如地震计记录的地震事件数量.在这种情况下,适用Poisson似然函数来描述.具体而言,假设观测数据为y,模型预测为,观测误差服从Poisson分布,那么Poisson似然函数可以表示为(赵志文等,2008):
其中,θ表示模型参数.
在某些情况下,地球物理学反演问题可以看作是非线性回归问题,此时可以采用残差平方和似然函数来描述.具体而言,假设观测数据为y,模型预测为,那么残差平方和似然函数可以表示为(彭军还,1994):
其中,θ表示模型参数.
需要注意的是,似然函数的具体形式取决于观测数据的类型、误差分布等因素.在选择似然函数时,需要根据具体问题的特点和已有的先验知识进行合理选择.

2.2.3 后验分布的计算

后验分布(posterior distribution)是贝叶斯统计学中的一个概念,用于描述在给定观测数据的条件下,参数的不确定性.具体来说,后验分布是在贝叶斯定理下,由先验分布和似然函数计算得出的参数分布(Gouveia and Scales, 1998).在贝叶斯统计中,通常关心的是参数θ在给定观测数据y的条件下的分布,即后验分布p(θ|y).根据贝叶斯定理,可以将后验分布表示为(Downton and Lines, 2001):
其中, p(θ|y)是似然函数,p(θ)是先验分布,p(y)是边缘似然函数,它可以表示为:
计算后验分布的方法主要有两种:解析计算和数值计算.在简单的模型中,可以通过解析方法得到后验分布.例如,在高斯模型中,如果观测数据和先验分布都服从正态分布,那么后验分布也是正态分布,可以直接计算均值和方差.另一种方法是使用共轭先验,这种先验分布具有良好的数学性质,可以确保后验分布与先验分布属于同一分布族,便于计算.例如,如果先验分布和似然函数构成的模型具有共轭性,那么后验分布也可以表示为同一分布族中的分布.
然而,在大多数实际问题中,后验分布往往无法通过解析方法计算得到.这时就需要使用数值方法来近似计算后验分布.常见的数值方法包括马尔可夫链蒙特卡罗方法(Markov Chain Monte Carlo, MCMC)(Gallagher et al., 2009)和变分推断方法(Blei et al., 2017).马尔可夫链蒙特卡罗方法通过构造一个马尔可夫链,在参数空间中进行随机游走,通过接受率准则来决定是否接受状态转移,从而得到参数样本集合,以近似后验分布.MCMC方法可以处理任意维度的参数空间,但计算效率较低,需要大量的状态转移和样本生成,并需要精细调整接受率准则.尽管如此,MCMC方法生成的后验模型样本有助于定量评估反演模型参数的不确定性.然而,在高维地球物理反演问题中,计算精度和计算效率总是相互制约的.
另一种近年来新兴的计算后验分布的方法是变分推断(variational inference)方法.其思想是通过寻找一个近似分布来近似后验分布.变分推断方法将后验分布表示为参数化的分布族,并通过最小化近似分布与真实后验分布的KL散度来求得最优近似分布的参数(Blei et al., 2017).相比MCMC方法,变分推断方法通常更高效,并且适用于高维参数空间和大规模数据集.然而,变分推断方法也存在一些局限性,例如可能存在偏差问题,且无法处理存在多个不同后验模态的情况.此外,还有一些其他的计算后验分布的方法,例如优化方法和采样重要性抽样方法等,每种方法都有其优缺点,需要根据具体问题选择适当的方法.

2.3 贝叶斯反演方法关键性技术

2.3.1 正则化方法

贝叶斯方法在地球物理学中的应用主要解决反演问题中的不确定性和稳定性问题,通过其正则化能力来实现.地球物理反演问题通常涉及大量参数,参数空间复杂且容易过度拟合,这会导致模型过于贴合观测数据而忽略真实模型的信息(陈晓和于鹏,2014).贝叶斯方法的正则化能力通过引入先验信息,减少参数空间的搜索范围,从而避免过度拟合的问题.此外,地球物理反演数据通常受到噪声、不完整或缺失的影响,这些因素可能导致反演结果不可靠或具有误导性.反演模型通常存在多个可能的解,这些解之间存在不确定性.贝叶斯方法的正则化能力可以通过引入先验信息来填补数据缺失或不完整的部分,从而提高反演结果的可靠性.因此,贝叶斯方法的正则化能力在地球物理学中具有重要意义,可以提高反演结果的可靠性和精度,同时减少不确定性,使得反演结果更加可靠和稳定.
拉普拉斯正则化和冯卡门正则化是贝叶斯反演中常见的正则化方法.拉普拉斯正则化是通过在目标函数中添加参数向量的二阶导数(拉普拉斯算子)的负数来实现的,以克服过度拟合的问题(Williams,1995).在反演问题中,拉普拉斯正则化通过在反演目标函数中添加拉普拉斯正则化项来控制解的平滑度和稳定性,得到更可靠和稳定的反演结果(Amey et al., 2018).冯卡门正则化是通过在目标函数中引入参数向量的先验分布来实现的,先验分布通常为高斯分布或其他分布(Hooper and Bekaert, 2014).这种正则化方法主要用于缩小参数空间,引入先验信息,从而控制反演结果的不确定性.在反演问题中,冯卡门正则化通过在反演目标函数中添加先验分布的负对数来约束和调节解.这两种方法都是常见的正则化方法,但在实现方式和应用场景上略有不同.拉普拉斯正则化主要用于控制解的平滑度和稳定性,而冯卡门正则化主要用于引入先验信息缩小参数空间,控制反演结果的不确定性.具体应用时,应根据反演问题的特点和需求选择适当的正则化方法.
滑动分布在揭示地震的多个属性方面起到重要作用,包括断层几何形态、应力变化、摩擦性质和未来潜在危险.滑动反演方法将滑动量表示为深度的函数.然而,通常这种方法建立在一些假设的基础上,例如以最小化滑动量的二阶空间导数(即拉普拉斯算子)为目标.这些假设可能缺乏明确的物理基础,且可能导致产生偏差.越来越多的证据表明,断层滑动具有分形特性.因此,Amey等(2018)提出将分形滑动作为正则化的一种形式纳入滑动反演中,而不采用拉普拉斯平滑,并开发了一种有效解决滑动问题的贝叶斯方法,其中包括冯卡门正则化(Amey et al., 2018).通过合成测试,该方法能更好地恢复分形滑动,即使在输入滑动不具有分形特性的情况下,Amey等(2018)的方法也表现更佳.他们利用干涉合成孔径雷达(InSAR)和全球定位系统(GPS)数据,并应用这种方法来研究2014年发生在纳帕谷地区的6.0级地震上的断层滑动情况.结果显示,冯卡门和拉普拉斯反演给出了类似的滑动幅度,但滑动位置不同,且冯卡门正则化解决方案的滑动置信区间远比拉普拉斯正则化解决方案更紧凑.
为了进行比较,研究人员还使用了具有拉普拉斯先验的贝叶斯方法来解决断层滑移问题(Minson et al., 2013).在这种方法中,他们同时解决了模型拟合数据、模型拟合拉普拉斯分布以及拉普拉斯平滑超参数的问题.图 2显示了模态解的结果,两种方法都表明滑动分布更多出现在断层西南段而非断层东北段上,峰值滑动超过1 m.然而,拉普拉斯解的峰值滑动略高,并且在西南断层上延伸更远.西南断层上的滑动标准偏差显著大于冯卡门正则化反演的标准偏差,表明拉普拉斯解在该区域存在更大的不确定性.这个结果并不令人意外,因为冯卡门先验提供了比拉普拉斯先验更紧密的约束.需要注意的是,这种更紧密的约束是基于对断层滑动特征的先验认识.冯卡门解在东北断层上略微减小了滑动幅度,并更好地拟合了该区域的数据.
图2 (a) 采用冯卡门先验和(b)拉普拉斯先验的贝叶斯反演的断层滑动分析(众数);(c,d)是对应的断层滑动分布标准偏差(修改自Amey et al., 2018)

两种方法反演的滑动分布大致相似,但最大滑动量位置有所不同,因为拉普拉斯在沿着走向方向上延伸了断层的滑动(如右图所示),并且在断层东北段上施加了更大的滑动.从标准偏差可以看出,拉普拉斯先验正则化的不确定性更明显.

Fig 2 The Bayesian inversion of fault slip using (a) von Mises and (b) Laplace priors (modes); (c, d) Corresponding standard deviation of the fault slip distribution (adapted from Amey et al., 2018)

The slip distributions retrieved by both methods are approximately similar, yet the location of maximum slip differs as Laplace extends slip along the strike direction of the fault (as shown in the right panel) and applies greater slip on the northeast fault. From the standard deviation, it is evident that the Laplace prior regularization exhibits more pronounced uncertainty.

拉普拉斯平滑和冯卡门正则化在反演滑动问题上表现出不同的特点.拉普拉斯平滑将滑动放置在不同的位置,而冯卡门正则化则在纳帕谷地震中提供了更紧密的置信度边界.这种更紧密的置信度边界是通过学者们施加的自相似先验获得的,该先验基于观测证据.因此,使用冯卡门正则化能够捕捉地震滑动的自相似特性.
冯卡门正则化是一种能够引入先验信息、更好地约束反演结果的方法.通过在目标函数中引入参数的先验分布,可以将反演结果限制在先验范围内,从而有效控制反演结果的不确定性,避免得到不现实的解.这对于解决高度非线性和高度欠定的反演问题尤为重要.相对于拉普拉斯正则化,冯卡门正则化的先验分布可以更灵活地约束参数的范围和分布.由于冯卡门正则化能够引入先验信息,因此能够更好地处理数据不完备或噪声较大的情况,从而得到更可靠和稳健的反演结果.Amey等(2018)通过对比分析拉普拉斯平滑和冯卡门正则化在反演断层滑动问题上的表现,发现冯卡门正则化提供了更紧密的置信度边界和更好的拟合效果.冯卡门正则化有助于提高反演结果的可靠性和精度,特别是在具有分形特性的滑动问题中.因此,在地球物理学中的反演研究中,选择合适的正则化方法,如冯卡门正则化,能够为解决问题提供更好的约束和稳定性.

2.3.2 采样方法

在贝叶斯反演方法中,采样方法是一种常用的技术,用于从参数的后验概率分布中抽取样本,以获取参数的概率分布信息(蒋星达等,2022).贝叶斯反演的目标是获得参数的后验概率分布,即在观测数据和先验信息下,参数可能取值的概率分布.采样方法能够生成大量样本,从后验概率分布中获取参数的概率分布信息,包括参数的平均值、方差、高阶统计量等,为地球内部结构和物性的推断提供更全面和准确的结果(余小东等,2020).此外,采样方法能够在高维参数空间中进行抽样,有效地探索参数空间中的不同取值,并在复杂的概率分布中获取高维参数的联合概率分布信息,从而提高反演的准确性和可靠性.
地球物理反演中的观测数据和先验信息通常导致参数的后验概率分布具有复杂的结构,例如非高斯性、多模态性等.采样方法通过随机抽样的方式对复杂的概率分布进行探索,并生成符合后验概率分布的样本,更好地捕捉参数的概率分布信息.此外,采样方法的收敛性和效率对贝叶斯反演的结果质量和计算速度具有重要影响.选择合适的采样方法对地球物理反演的成功至关重要,它直接影响了反演结果的质量、可靠性和计算效率.目前,在贝叶斯反演方法中常用的采样方法为马尔科夫链蒙特卡洛法.MCMC是一种基于马尔科夫链的采样方法,通过在参数空间中随机漫游来生成参数的样本(Doucet et al., 2000),它们可以在高维参数空间中进行采样,并在概率分布存在复杂结构时表现较好.吉布斯采样(Gibbs Sampling)是一种特殊的MCMC方法,用于多维参数空间中的条件分布的抽样.它通过依次从各个参数的条件概率分布中抽取样本,逐步在多维参数空间中收敛到后验概率分布(Bagnardi and Hooper, 2018).
Bagnardi和Hooper(2018)提出了一种贝叶斯反演方法,用于从多个测地数据集中快速获得对震源模型参数的后验概率密度函数(PDFs).这一方法采用了马尔可夫链蒙特卡洛方法(MCMC)来高效地对后验概率密度函数进行采样,其中包括Metropolis-Hastings算法(Hastings, 1970).这种方法还具备自动步长选择功能,可应用于多种地震反演问题(图 3).作为案例研究,该方法已成功应用于2015年中国皮山地震(MW6.4)的震源参数反演,它通过InSAR地表形变数据约束,获取最优震源参数以及与之相关的不确定性(Ainscoe et al., 2017).Metropolis-Hastings算法可以有效控制采样,在足够大的迭代次数后,样本的分布近似于后验分布.它首先从先验分布p(m)中选择一个初始模型参数集 mi=0,可以任意选择,也可以使用直接搜索方法(如模拟退火)预先估计得到.然后估计关联的似然函数p(d|mi).通过在p(m)内随机生成新的模型参数集合,生成一个新的模型试验.当每个模型参数的先验分布是均匀且独立的时候,可以通过将mi中的每个参数按照anΔmj的大小扰动,其中an是从范围[-1, 1]内均匀分布生成的随机值,Δmj是每个参数mj的最大随机步长.如果新模型试验中的模型参数超出了均匀先验概率的范围,可以通过以下方法进行调整(Bagnardi and Hooper, 2018):
图3 (a) 收敛的马尔可夫链和(b)未收敛的马尔可夫链的示例追踪图(修改自Bagnardi and Hooper, 2018)

前20000次迭代也被丢弃作为调整最大步长的燃烧期的一部分,因为变化的步长可能导致对后验概率密度函数的非典型抽样.

Fig 3 (a) Example trace plot of a converged Markov chain and (b) example trace plot of a non-converged Markov chain (modified from Bagnardi and Hooper, 2018)

The first 20000 iterations were also discarded as part of the burn-in period for adjusting the maximum step length, as variable step lengths can lead to non-typical sampling of the posterior probability density function.

其中,UBjLBj分别是给定模型参数mj的均匀概率的上界和下界.然后计算新一组模型参数p(d|mi+1) 的似然函数,如果其值大于先前的似然函数,即p(d|mi+1)>p(d|mi),则接受这一步骤,并保留试验模型值.如果新的似然函数值小于先前的似然函数值,即p(d|mi+1) < p(d|mi),则按照新似然函数和先前似然函数的比值,以概率b的方式接受这一步骤,其中b是从范围为[0, 1]的均匀分布中随机选择的值.否则,保留先前的模型参数集,将mi+1设置为mi.然后,i递增,并从选择新的试验模型参数开始重复上述步骤.这种方法确保了能够保留改善当前模型概率的模型参数集,同时也能够保留概率较低的模型参数集,从而允许算法逃离局部极小值.该过程重复进行,直到获得对后验分布的代表性抽样(例如,进行105到107次迭代).
Metropolis-Hastings算法的效率受到随机步长Δmj的选择影响,这些步长必须考虑到参数的先验分布.步长过小可能导致收敛速度过慢,步长过大则可能导致接受率降低.为了最大化算法的效率和确保适当的收敛性,通常在前几万次迭代中会自动调整每个模型参数的步长大小,以进行敏感性测试.此外,检查马尔科夫链的收敛性非常重要,通常使用迹线图等工具进行可视化检查.燃烧期的确定是关键一步,它表示从初始模型参数到具有合理似然函数值的模型参数的过渡过程,通常涉及数千次迭代.Hooper和Bekaert(2014)展示了迹线图的示例,显示了经过燃烧期后已经收敛的链和未收敛的链.最后,可以依次或并行运行多个独立的MCMC反演,并将结果组合成最终的后验概率分布,这是组合独立马尔科夫链的结果(Anderson and Poland, 2016).这一步还可以用来检查所有链的后验概率分布是否收敛到类似的分布,而不依赖于选择的初始模型参数集.
在InSAR数据反演领域,一种常见的替代方法用于估计模型参数的不确定性,这一方法曾被用于对2015年喀什地震的研究(Ainscoe et al., 2017).Wright等(2003)使用非线性反演技术寻找在给定数据添加模拟噪声后的情况下的最佳拟合解的加权最小二乘法.这个过程重复进行多次(例如102次),使用不同的数据噪声模拟,得到的解的分布表示模型参数的后验概率分布.虽然这种方法和Marco等人的方法在所有先验都是均匀的情况下应该提供类似的结果,但Hooper和Wright (2009)表明,与提出的MCMC方法类似,贝叶斯反演方法在需要运行的前向模型数量方面更高效.因此,在快速响应或操作性环境中(例如火山观测站)以及需要快速和稳健地估计源参数的情况下,贝叶斯反演方法可能更可取.

2.4 贝叶斯反演的参数及模型选取

2.4.1 模型参数

在贝叶斯反演方法中,模型参数通常包括待估计的地球物理或地质参数,这些参数用于描述地球系统的特性.具体而言,模型参数可以包括但不限于以下几类:震源参数,如地震的震源深度、震源位置、震源时刻等参数,用于描述地震的发生和演化过程(颜丙囤等,2023);介质参数,如地下介质的速度、密度、衰减等参数,用于描述地球内部的结构和性质(薛清峰,2021);界面参数,如地层的位置、倾角、分界面形状等参数,用于描述地下不同地层之间的界面形态(周龙泉等,2006);地形参数,如地表地形的高程、坡度、曲率等参数,用于描述地表地形的形态特征(刘甲美等,2015);地球化学参数,如地下岩石的化学成分、岩石圈的热流、地下水的含水量等参数,用于描述地球物质的化学性质和热力学特性(赵振燊,2012);反演模型参数,如反演模型中的网格节点位置、节点属性值等参数,用于描述反演模型的细节和复杂性(伍吉仓和陈永奇,1997).在贝叶斯反演中,这些模型参数需要通过观测数据和先验信息进行估计,以得到最优的参数估计结果和参数不确定性.模型参数之间通常存在耦合关系,即一个参数的取值可能会对其他参数的估计结果产生影响.贝叶斯反演方法可以考虑参数之间的相互关系,将耦合关系纳入模型中,以更好地估计参数的值和不确定性(雷文文,2012).在反演过程中通常存在多种不同的模型可以用于解释观测数据.贝叶斯反演方法可以通过引入模型先验概率来进行模型选择,从而选择最合适的模型.在地球物理反演中对于模型参数的估计和不确定性的处理非常重要,可以提高反演结果的准确性、可靠性和实时性.
在面对非线性问题时,Mosegaard和Tarantola (1995)通过MCMC方法有效地实现了非线性后验分布的采样,既保持了反演问题的非线性特性,又通过全局搜索提高了解的准确性.然而,由于模型解空间通常很大,MCMC方法需要高计算成本,这可能影响到马尔科夫链的收敛速度.如果可以使用统计参数来描述模型参数,就可以显著减小解空间的范围,从而提高反演效率.通过地质统计学的研究,可以有效地描述空间变化现象.Deutsch和Journel (1992)利用地质统计学方法对未知模型参数进行估计,并使用概率分布来表示不确定性.Hansen等(2006)将地质统计学与贝叶斯反演方法相结合,通过地质统计方法进行先验描述,降低了反演的复杂性.张繁昌等(2014)提出了一种地震数据约束的贝叶斯随机反演方法,结合贝叶斯理论和地质统计学.该方法首先建立了先验解空间,然后通过马尔科夫链扰动模拟算法,采样后验概率分布,找到最佳匹配解,同时生成大量样本.这允许对地层反演参数进行评估和不确定性分析.Ragon等(2018)研究了震源反演中与断层几何形态不确定性相关的影响.他们关注断层结构参数,如位置、走向和倾角,对反演结果的影响.研究中提出了一个实用框架,将断层几何形态的不确定性考虑进反演问题,并结合了预测误差和观测误差的失配协方差矩阵.通过贝叶斯和优化的静态滑动反演方法验证了这一方法.H. Vasyura-Bathke等(2021)指出地震学常采用简化模型估计地震参数,但这可能导致不确定性.提出综合的贝叶斯方法,同时估计多个参数,强调处理不确定性和噪声对于准确估计至关重要.文章讨论了地震数据受到的噪声干扰,以及如何估计理论误差.最近的研究集中于已知源几何和位置的线性反演,对于包括非线性参数的改进尚不明确.然而,BEAT软件(Vasyura-Bathke et al., 20192020)采用了通用的贝叶斯框架进行不确定性量化,并提供了多种噪声模型的选择,并已成功应用于具有任意频率内容、质心位置、站点距离和源参数化的数据.用户可以自由更改BEAT或根据需要应用于本地、区域和全球数据,使用DC、DV或完整的CMTs、矩形断层和多段有限断层.
贝叶斯反演方法在地球物理领域扮演重要角色,用于估计多种模型参数,包括地震源、介质、地层、地形、地球化学、反演模型等.这些参数相互关联,因此贝叶斯方法提供了精确的估计和不确定性信息,提高反演结果的准确性和可靠性.近期研究集中于解决非线性问题,并引入MCMC方法来有效采样非线性后验分布,但由于计算成本较高,需要提高计算效率(Bagnardi and Hooper, 2018).地质统计学方法有效描述空间变化,结合贝叶斯反演,可减少复杂性.贝叶斯方法在地震学中也用于选择最适模型.部分研究关注了断层几何不确定性对震源反演的影响,提出实用框架(Vasyura-Bathke et al., 2019).还有部分学者强调了噪声和理论误差的处理,以提高地震参数估计的准确性(Vasyura-Bathke et al., 2021).然而,对于包括非线性参数的更进一步改进仍需深入研究.

2.4.2 模型的选取

在地球物理学中,构建有限断层模型已成为强震后进行断层滑动反演的常规方法.尽管地壳形变和地震观测数据的质量和数量不断增加,对于给定事件推断地下断层滑动仍存在显著的变异性.从地表形变来估计断层滑动的分布,本质上是一个非限定的反问题,因为可能存在多个模型可以同时解释数据.因此,同一地震的不同有限断层模型通常会展示出显著的差异性.在过去10年中,地球物理学界一直在努力研究并表征这种模型的变异性(Mai et al., 2016).
随着完整数据集的可获得性,对1992年MW7.3的兰德斯地震已经发布了多个有限断层模型(Murray et al., 1993; Freymueller et al., 1994; Hudnut et al., 1994; Xu et al., 2016),这些模型的推断显示了相似的滑动分布特征,即主要的滑动发生在断裂中央——Homestead Valley断层.然而,这些模型也存在显著的差异.Gombert等(2018)使用贝叶斯方法研究了1992年兰德斯地震,并提出了一个新的震后模型.他们通过联合反演多个观测数据来评估地下断层滑移分布的不确定性,尤其关注浅层区域的滑移.这些不确定性问题主要来自滑动反演的不适定性和使用平滑或正则化约束的物理限制.为减小模型不确定性,他们设计了一个三维断层几何模型,考虑了弹性结构的不确定性,并生成了合理模型的集合来适应观测数据和不确定性.研究结果表明,浅层滑移不足与未考虑断层附近破损带的剪切模量减小有关.
Gombert等(2018)生成了50万个代表滑动分布后验信息的模型,以及一个后验平均模型,在图(5)中,展示了后验平均模型(即所有采样模型的平均值),以及95%置信度椭圆.后验平均模型是一种常见选择,因为贝叶斯方法考虑集合解决方案,而不是单个模型.
图 4观察到整个断层系统主要是走滑运动,大部分滑动集中在破裂的中央和北部,峰值振幅约为11 m.这些特征与之前的结果相当一致,尽管已发布的模型具有较低的最大滑动幅度(Fialko,2004; Xu et al., 2016)这种差异可能是由于以前的研究中施加的平滑导致最大滑动幅度降低.这一研究基于丰富的地震测量数据、准确的不确定度估计和逼真的断层几何形态,为兰德斯地震提供了一个可信的随机有限断层模型.该研究采用了贝叶斯反演方法,具有两个主要优点:(1)该方法不受任何平滑方法偏差的影响;(2)可以获得后验参数不确定性,从而提供有关模型有效性的有价值信息.该解决方案的预测与各种观测结果表现良好.
图4 后续平均余震滑动模型(修改自Gombert et al., 2018)

每个子断层补丁的颜色表示滑动幅度.箭头及其相关的95%置信度椭圆表示滑动方向和不确定性.左下角的插图显示了按深度分的规范化破裂能.整个断层系统和各个断层段的浅滑动缺陷(SSD)的概率密度函数(PDF)被呈现.同一图上的垂直线表示两个已发表模型(Cotton and Campillo, 1995; Fialko,2004)的SSD.

Fig 4 Shows the time-averaged postseismic sliding model (modified from Gombert et al., 2018)

The color of each subfault patch represents the magnitude of sliding. Arrows and their associated 95% confidence ellipses represent the direction of sliding and its uncertainty. The inset in the bottom left corner displays the normalized rupture energy distribution as a function of depth. The Probability Density Functions (PDFs) of Shallow Slip Deficit (SSD) for the entire fault system and individual fault segments are presented. Vertical lines on the same figure represent the SSDs of two published models (Cotton and Campillo, 1995; Fialko, 2004).

此外,在地震学中,使用真实的地球模型通常需要计算数值或合成的格林函数(GFs),因为解析解很少可用.合成GFs的计算是具有挑战性的,因此在优化的每次迭代中动态重新计算这些GFs是耗时且不切实际的.为了应对这个问题,Heimann等(2019)提出了Pyrocko-GF——一个基于Python的框架和工具包,可以预先计算和高效存储GFs.该框架与计算方法和GF类型无关,允许用户在优化过程中快速获取存储的GF,进行比较、可视化和质量检查.它支持各种地震学应用,包括合成地震波形、静态位移等,并且可以扩展到其他地球物理现象的模拟.Pyrocko-GF提供了Python API,方便用户进行地球物理过程的正演建模.在贝叶斯反演中,模型的选择对于准确估计地震学中的格林函数(GFs)非常关键.Pyrocko-GF框架提供了预先计算和高效存储GFs的功能,使用户能够比较、验证和质量检查各种模型.这对于贝叶斯反演中确定最合适的地震模型非常有意义,因为用户可以使用这些预先计算的GFs进行模型对比和评估,从而选择最优模型并提高反演的准确性和可靠性.

2.5 贝叶斯方法的局限性和未来发展方向

贝叶斯反演方法是一种基于概率论和统计学理论的反演方法,它具有许多优点,但同时也存在一些局限性和问题.首先,贝叶斯反演方法需要确定先验分布,即反演参数的概率分布,而这个先验分布往往需要依赖领域专家或者已知数据的先验知识.如果先验分布的选取不合适或者不准确,将会对反演结果产生较大的影响,导致反演结果的不准确性.其次,贝叶斯反演方法需要进行高维积分计算,这是计算成本高、效率低的一个主要问题.在高维问题中,随着参数维度的增加,积分计算的复杂度呈指数增长,因此在实际应用中,通常需要对参数空间进行采样和离散化,从而降低计算复杂度,但这也可能会影响反演结果的准确性.此外,贝叶斯反演方法还需要处理噪声数据和不确定性信息,这也是一个比较困难的问题.对于噪声数据的处理,需要对数据进行预处理、滤波等操作;而对于不确定性信息的处理,则需要确定合适的先验分布和模型,同时考虑模型误差、测量误差和数值误差等因素.最后,贝叶斯反演方法的应用范围和适用性也存在限制.由于贝叶斯反演方法需要处理大量的数据和参数信息,因此对于数据量小、参数空间简单的问题,贝叶斯反演方法可能会出现过拟合等问题.所以,贝叶斯反演方法虽然具有很多优点,但同时也存在一些局限性和问题,需要根据实际情况和问题类型选择合适的反演方法和技术.
在固体地球物理学中,贝叶斯反演方法已经被广泛应用于参数反演问题.未来的发展方向和研究重点可以聚焦于以下几个方面:首先,随着测量技术的不断发展,实际观测数据的质量和数量都将不断提高,如何在这种情况下更好地利用数据信息并结合先验知识来推断地下结构参数是一个重要的研究方向.其次,针对反演模型中的非线性问题,如何采用高效的计算方法和技术来实现反演参数空间的快速探索和优化也是一个关键的研究领域.第三,如何选择或设计适当的先验模型,以及如何根据实际观测数据来更新先验模型,都是贝叶斯反演方法中的重要问题.未来的研究可以探索新的先验模型,以及将贝叶斯反演方法与其他反演方法结合起来,以更好地实现参数反演.最后,贝叶斯反演方法在固体地球物理学中的应用还有许多挑战和机遇,如海洋地球物理学中的参数反演、非线性反演问题等,这些领域的研究将有助于推动贝叶斯反演方法的发展和应用.

3 结论与展望

贝叶斯反演算法在各领域,包括地球物理学、信号处理、计算机视觉和自然语言处理等,都得到了广泛应用和不断发展.它以其独特的特点在参数估计和模型反演中发挥着关键作用.首先,贝叶斯反演算法在参数不确定性的量化方面具有更强大的能力.它不仅提供点估计,还能够给出完整的后验概率分布,包括参数的概率密度函数.这使得我们能够更全面地了解参数估计的不确定性,而不仅仅是一个确定的数值.这对于决策制定和风险评估至关重要.其次,贝叶斯反演算法允许更好地利用已知的先验信息.在很多情况下,我们对问题的先验知识是有限的,而贝叶斯方法可以有效地将这些知识与观测数据相结合.这有助于提高估计的准确性,尤其是在观测数据有限或噪声较大的情况下.此外,贝叶斯反演算法能够适应不同的噪声模型.它可以处理各种噪声类型,包括高斯噪声、泊松噪声、混合噪声等.这种灵活性使其在各种实际应用场景中都能够有效地工作.贝叶斯方法还能够处理非线性反演问题,这在现实世界中非常常见.与传统的线性方法相比,它具有更广泛的适用性.
然而,贝叶斯反演算法的一个挑战是计算效率问题,尤其在高维反演和大规模数据情况下.解决这个问题需要开发更高效的多参数和高维反演算法,以确保在合理的时间内得到结果.联合反演一直是地球物理反演领域的研究热点.它的目标是将多种类型的数据相互融合,以提供更准确的反演结果.这种联合反演本质上是引入更多的先验信息,可以帮助缩小解的非唯一性,提高反演结果的可靠性.因此,如何建立地球物理联合反演模型,实现多种地球物理数据的有效融合,并通过地质构造信息的量化获取更多的先验信息,是未来需要关注的问题.
总的来说,贝叶斯反演算法在各个领域都具有广泛的应用前景,但也需要克服一些挑战,包括计算效率和联合反演的问题.通过不断的研究和方法改进,贝叶斯方法将在更多应用中发挥关键作用,提高参数估计的准确性和稳定性.

感谢审稿专家提出的修改意见和编辑部的大力支持!

Ainscoe E A , Elliott J R , Copley A . Blind thrusting, surface folding, and the development of geological structure in the MW6.3 2015 Pishan (China) earthquake. Journal of Geophysical Research: Solid Earth, 2017, 122(11 9359 9382.

DOI

Alemie W , Sacchi M D . High-resolution three-term AVO inversion by means of a Trivariate Cauchy probability distribution. Geophysics, 2011, 76(3): R43-R55.

DOI

Amey R M J , Hooper A , Walters R J . A Bayesian method for incorporating self-similarity into earthquake slip inversions. Journal of Geophysical Research: Solid Earth, 2018, 123(7): 6052-6071.

DOI

Ammon C J , Ji C , Thio H K . Rupture process of the 2004 Sumatra-Andaman earthquake. Science, 2005, 308(5725): 1133-1139.

DOI

Anderson K R , Poland M P . Bayesian estimation of magma supply, storage, and eruption rates using a multiphysical volcano model: Kīlauea Volcano, 2000-2012. Earth and Planetary Science Letters, 2016, 447: 161-171.

DOI

árnadóttir T , Hreinsdóttir S , Gudmundsson G . Crustal deformation measured by GPS in the South Iceland Seismic Zone due to two large earthquakes in June 2000. Geophysical Research Letters, 2001, 28(21): 4031-4033

Avseth P , Mukerji T , Jørstad A . Seismic reservoir mapping from 3-D AVO in a North Sea turbidite system. Geophysics, 2001, 66(4): 1157-1176.

DOI

Bagnardi M , Hooper A . Inversion of surface deformation data for rapid estimates of source parameters and uncertainties: A Bayesian approach. Geochemistry, Geophysics, Geosystems, 2018, 19(7): 2194-2211.

DOI

Barrientos S E , Ward S N . The 1960 Chile earthquake: inversion for slip distribution from surface deformation. Geophysical Journal International, 1990, 103(3): 589-598.

DOI

Blei D M , Kucukelbir A , Mcauliffe J D . Variational inference: A review for statisticians. Journal of the American Statistical Association, 2017, 112(518): 859-877.

DOI

Bodin T , Sambridge M , Gallagher K . A self-parametrizing partition model approach to tomographic inverse problems. Inverse Problems, 2009, 25(5): 055009.

DOI

Chen X , Zhou L , Yu P . Bayesian joint inversion of magnetotelluric and seismic data based on simulated annealing method. Journal of East China University of Technology (Natural Science), 2016, 39(1): 59-66

Chunduru R K , Sen M K , Stoffa P L . Hybrid optimization methods for geophysical inversion. Geophysics, 1997, 62(4): 1196-1207.

DOI

Cotton F , Campillo M . Frequency domain inversion of strong motions: Application to the 1992 Landers earthquake. Journal of Geophysical Research: Solid Earth, 1995, 100(B3): 3961-3975.

DOI

Dai Q W , Chen W , Zhang B . Improved particle swarm optimization and its application to full-waveform inversion of GPR. Geophysical and Geochemical Exploration, 2019, 43(1): 90-99

de Figueiredo L P , Grana D , Roisenberg M . Multimodal Markov chain Monte Carlo method for nonlinear petrophysical seismic inversion. Geophysics, 2019, 84(5): M1-M13.

DOI

Deutsch C V , Journel A G . Geostatistical software library and user's guide. New York, 1992, 119(147): 578

Dorugade A V . New ridge parameters for ridge regression. Journal of the Association of Arab Universities for Basic and Applied Sciences, 2014, 15: 94-99

Doucet A , Godsill S , Andrieu C . On sequential Monte Carlo sampling methods for Bayesian filtering. Statistics and Computing, 2000, 10(3): 197-208.

DOI

Feng J , Meng X H , Chen Z X . Review of the gravity density interface inversion. Progress in Geophysics, 2014, 29(1): 223-228.

DOI

Feng W P , Li Z H . A novel hybrid PSO/simplex algorithm for determining earthquake source parameters using InSAR data. Progress in Geophysics, 2010, 25(4): 1189-1196.

DOI

Fialko Y . Evidence of fluid-filled upper crust from observations of postseismic deformation due to the 1992 MW7.3 Landers earthquake. Journal of Geophysical Research: Solid Earth, 2004, 109(B8): B08401.

DOI

Freymueller J , King N E , Segall P . The co-seismic slip distribution of the Landers earthquake. Bulletin of the Seismological Society of America, 1994, 84(3): 646-659.

DOI

Fukuda J , Johnson K M . A fully Bayesian inversion for spatial distribution of fault slip with objective smoothing. Bulletin of the Seismological Society of America, 2008, 98(3): 1128-1146.

DOI

Gallagher K , Charvin K , Nielsen S . Markov chain Monte Carlo (MCMC) sampling methods to determine optimal models, model resolution and model choice for Earth Science problems. Marine and Petroleum Geology, 2009, 26(4): 525-535.

DOI

Gombert B , Duputel Z , Jolivet R . Revisiting the 1992 Landers earthquake: a Bayesian exploration of co-seismic slip and off-fault damage. Geophysical Journal International, 2018, 212(2): 839-852.

DOI

Gouveia W P , Scales J A . Bayesian seismic waveform inversion: Parameter estimation and uncertainty analysis. Journal of Geophysical Research: Solid Earth, 1998, 103(B2): 2759-2779.

DOI

Guitton A , Symes W W . Robust inversion of seismic data using the Huber norm. Geophysics, 2003, 68(4): 1310-1319.

DOI

Han Y P , Luo S M , Chen F C . Effect of model error and result of regional vertical crustal movement. Journal of Geodesy and Geodynamics, 2015, 35(1): 21-25

Hansen T M , Journel A G , Tarantola A . Linear inverse Gaussian theory and geostatistics. Geophysics, 2006, 71(6): R101-R111.

DOI

Hastings W K . Monte Carlo sampling methods using Markov chains and their applications. Biometrika, 1970, 57(1): 97-109.

DOI

Heimann S , Vasyura-Bathke H , Sudhaus H . A Python framework for efficient use of pre-computed Green's functions in seismological and other physical forward and inverse source problems. Solid Earth, 2019, 10(6): 1921-1935.

DOI

Hu G W , Niu Y , Yang S L . New genetic algorithm for MSA solution problem. Journal of Lanzhou University of Technology, 2005, 31(3): 91-93.

DOI

Hudnut K W , Bock Y , Cline M . Co-seismic displacements of the 1992 Landers earthquake sequence. Bulletin of the Seismological society of America, 1994, 84(3): 625-645.

DOI

Jiang X D , Zhang W , Yang H . The research on Bayesian inference for geophysical inversion. Reviews of Geophysics and Planetary Physics, 2022, 53(2): 159-171

Kirkpatrick S , Gelatt C D Jr , Vecchi M P . Optimization by simulated annealing. Science, 1983, 220(4598): 671-680.

DOI

Li Z C , Zhang P , Jin S G . Wenchuan earthquake deformation fault inversion and analysis based on GPS observations. Acta Geodaetica et Cartographica Sinica, 2009, 38(2): 108-113 108-113, 119

DOI

Liang Y , Zhu G B , Liu H T . Regularization of the earth gravity field determination using satellite gravity gradiometry data. Progress in Geophysics, 2019, 34(4): 1323-1327.

DOI

Liu J M , Gao M T , Chen K . On the correlation of ground motion parameters with slope stability incorporating topographic effects. Acta Seismologica Sinica, 2015, 37(5): 865-874

Lu Y X , Peng S P , Du W F . Rayleigh wave inversion using heat-bath simulated annealing algorithm. Journal of Applied Geophysics, 2016, 134: 267-280.

DOI

Luo T , Feng X , Guo Z Q . Seismic inversion of anisotropy parameters of fractured reservoirs by simulated annealing and particle swarm optimization. Journal of Jilin University (Earth Science Edition), 2019, 49(5): 1466-1476

Mai P M , Schorlemmer D , Page M . The earthquake-source inversion validation (SIV) project. Seismological Research Letters, 2016, 87(3): 690-708.

DOI

Malinverno A , Briggs V A . Expanded uncertainty quantification in inverse problems: Hierarchical Bayes and empirical Bayes. Geophysics, 2004, 69(4): 1005-1016.

DOI

Michalos G , Karvouniari A , Dimitropoulos N . Workplace analysis and design using virtual reality techniques. CIRP Annals, 2018, 67(1): 141-144.

DOI

Minson S E , Simons M , Beck J L . Bayesian inversion for finite fault earthquake source models Ⅰ—theory and algorithm. Geophysical Journal International, 2013, 194(3): 1701-1726.

DOI

Monelli D , Mai P M , Jónsson S . Bayesian imaging of the 2000 Western Tottori (Japan) earthquake through fitting of strong motion and GPS data. Geophysical Journal International, 2009, 176(1): 135-150.

DOI

Mosegaard K , Tarantola A . Monte Carlo sampling of solutions to inverse problems. Journal of Geophysical Research: Solid Earth, 1995, 100(B7): 12431-12447.

DOI

Murray M H , Savage J C , Lisowski M . Coseismic displacements: 1992 Landers, California, earthquake. Geophysical Research Letters, 1993, 20(7): 623-626.

DOI

Pei Z L , Mou Y G . Numerical simulation of seismic wave propagation. Progress in Geophysics, 2004, 19(4): 933-941.

DOI

Peng J H . Sum maximum likelihood estimate: a new estimate rule. Journal of Guilin College of Geology, 1994, 14(4): 436-444

Ragon T , Sladen A , Simons M . Accounting for uncertain fault geometry in earthquake source inversions-Ⅰ: theory and simplified application. Geophysical Journal International, 2018, 214(2): 1174-1190.

DOI

Rothman D H . Automatic estimation of large residual statics corrections. Geophysics, 1986, 51(2): 332-346.

DOI

Shi X M , Wang J Y . Lecture on non-linear inverse methods in geophysics (4) - Genetic algorithm method. Chinese Journal of Engineering Geophysics, 2008, 5(2): 129-140.

DOI

Theune U , Jensås I Ø , Eidsvik J . Analysis of prior models for a blocky inversion of seismic AVA data. Geophysics, 2010, 75(3): C25-C35.

DOI

Vasyura-Bathke H , Dettmer J , Steinberg A . The Bayesian earthquake analysis tool. Seismological Research Letters, 2020, 91(2A): 1003-1018.

DOI

Vasyura-Bathke H , Dettmer J , Dutta R . Accounting for theory errors with empirical Bayesian noise models in nonlinear centroid moment tensor estimation. Geophysical Journal International, 2021, 225(2): 1412-1431.

DOI

Visser G , Guo P , Saygin E . Bayesian transdimensional seismic full-waveform inversion with a dipping layer parameterization. Geophysics, 2019, 84(6): R845-R858.

DOI

Wang F , Chen S C . A new strategy of solving non-linear inverse problem. Geophysical Prospecting for Petroleum, 1999, 38(3): 76-85

Wang L Y . Research on theory and application of total least squares in geodetic inversion. Acta Geodaetica et Cartographica Sinica, 2012, 41(4): 629

Wang W P , Wang Q L . Dislocation parameters of Gonghe earthquake jointly inferred by genetic algorithm and least square method. Acta Seismologica Sinica, 1999, 21(3): 285-290.

DOI

Wang X Q , Ding X , Wang L . A study on fast earthquake loss assessment and its application to 2008 Wenchuan M8 earthquake. Acta Seismologica Acta, 2009, 31(2): 205-211

Wang Z W , Fu L Y , Han T C . Review of thermoelasticity theory in rocks and its applications in geophysics. Reviews of Geophysics and Planetary Physics, 2021, 52(6): 623-633

Wen Y M , Xu C J , Liu Y . Source parameters of 2008 Qinghai Dachaidan MW6.3 earthquake from InSAR inversion and automated fault discretization method. Geomatics and Information Science of Wuhan University, 2012, 37(4): 458-462

Williams P M . Bayesian regularization and pruning using a Laplace prior. Neural Computation, 1995, 7(1): 117-143.

DOI

Wright T J , Lu Z , Wicks C . Source model for the MW6.7, 23 October 2002, Nenana mountain earthquake (Alaska) from InSAR. Geophysical Research Letters, 2003, 30(18): 1974.

DOI

Wu J C , Chen Y Q . Inversion of model parameters by considering prior information: Determination of the dislocation model parameters of Lancang-Gengma earthquakes. Crustal Deformation and Earthquake, 1997, 17(2): 27-32

Xie S L , Sun Y J . The co-seismic and post-seismic deformation and stress of Wenchuan earthquake in 2008 and their influence on the Songpinggou Geohazards. Acta Geoscientica Sinica, 2021, 42(5): 691-700

Xu C J , Wen Y M . Nonhomogeneity of the crust from MS7.9 Manyi(Tibet) earthquake with InSAR observation. Geomatics and Information Science of Wuhan University, 2008, 33(8): 846-849

Xu X H , Tong X P , Sandwell D T . Refining the shallow slip deficit. Geophysical Journal International, 2016, 204(3): 1867-1886

Xue Q F . Wave equation based joint inversion of microseismic source location, source excitation time and anisotropic parameters of the VTI media. Earthquake Research in China, 2021, 37(2): 285-299

Yan B T , Yin H T , Feng B . Source parameters and slip distribution of the Luding MS6.8 earthquake in 2022 constrained by InSAR data. Journal of Geodesy and Geodynamics, 2023, 43(3): 221-225

Yang G Y , Luan X W , Meng F S . Seismic data reconstruction based on Shearlet transform and total generalized variation regularization. Chinese Journal of Geophysics, 2020, 63(9): 3465-3477.

DOI

Yang J H , Liu Q H . An inversion method for stratified media based on hybrid optimization algorithm. Modern Radar, 2020, 42(5): 46-51

Yang L Q , Song H B , Hao T Y . Application of impedance inversion based on BP neural network. Progress in Geophysics, 2005, 20(1): 34-37

Yao M , Gao G , Hu R Q . Improved Bayesian inversion algorithm. Progress in Geophysics, 2020, 35(5): 1911-1918.

DOI

Yi Y Y , Wang J Y . Lectures on non-linear inversion methods of geophysical data (10) — Particle swarm optimization inversion method. Chinese Journal of Engineering Geophysics, 2009, 6(4): 385-389

Yin X Y , Zhang S X . Bayesian inversion for effective pore-fluid bulk modulus based on fluid-matrix decoupled amplitude variation with offset approximation. Geophysics, 2014, 79(5): R221-R232

Yu X D , Lu C D , Wang X B . Adaptive trans-dimensional Bayesian inversion of airborne time-domain electromagnetic data. Progress in Geophysics, 2020, 35(5): 2023-2032.

DOI

Zhang D L , Tao X X , Zhou Z H . Comparison of green's function of near-field ground motion due to analytical method and numerical method. Northwestern Seismological Journal, 2004, 26(3): 199-205

Zhang F C , Yin X Y , Wu G C . Impedance inversion by using annealing neural network. Journal of the University of Petroleum, China, 1997, 21(6): 16-18 16-18, 23

Zhang F C , Xiao Z B , Yin X Y . Bayesian stochastic inversion constrained by seismic data. Oil Geophysical Prospecting, 2014, 49(1): 176-182

Zhao Z W , Wang D H , Li H . Approximate Bayesian estimation for parameters in Poisson regression under quadratic loss function. Journal of Jilin University (Science Edition), 2008, 46(5): 836-840

Zheng Y , Ma H S , J . Source mechanism of strong aftershocks (MS≥5.6) of the 2008/05/12 Wenchuan earthquake and the implication for seismotectonics. Science in China Series D: Earth Sciences, 2009, 52(6): 739-753

Zhou L Q , Liu F T , Chen X F . Simultaneous tomography of 3-D velocity structure and interface. Chinese Journal of Geophysics, 2006, 49(4): 1062-1067

Zhou S M , Liu J X , Sun H L . Study of one-dimensional magnetotelluric regularized inversion based on non-negative least squares method. Chinese Journal of Engineering Geophysics, 2017, 14(3): 253-261

, , . 基于模拟退火算法的大地电磁和地震数据贝叶斯联合反演. 东华理工大学学报(自然科学版), 2016, 39(1): 59-66

前伟 , , . 改进型粒子群算法及其在GPR全波形反演中的应用. 物探与化探, 2019, 43(1): 90-99

, 小红 , 召曦 . 重力密度界面反演方法研究进展. 地球物理学进展, 2014, 29(1): 223-228.

DOI

万鹏 , 振洪 . InSAR资料约束下震源参数的PSO混合算法反演策略. 地球物理学进展, 2010, 25(4): 1189-1196.

DOI

月萍 , 三明 , 阜超 . 区域地壳垂直运动研究中的模型误差及其影响. 大地测量与地球动力学, 2015, 35(1): 21-25

宗海 . 中国南北地震带中北部最小完整性震级的探讨. 西北大学学报(自然科学版), 1994, 24(5): 411-416

桂武 , , 胜良 . 求解MSA问题的新型遗传算法. 兰州理工大学学报, 2005, 31(3): 91-93

星达 , , . 地球物理反演问题中的贝叶斯方法研究. 地球与行星物理论评, 2022, 53(2): 159-171

志才 , , 双根 . 基于GPS观测数据的汶川地震断层形变反演分析. 测绘学报, 2009, 38(2): 108-113 108-113, 119

, 广彬 , 洪涛 . 卫星重力梯度数据确定地球重力场的正则化方法研究. 地球物理学进展, 2019, 34(4): 1323-1327.

DOI

甲美 , 孟潭 , . 地形效应影响下地震动参数与斜坡稳定性的相关性研究. 地震学报, 2015, 37(5): 865-874

, , 智奇 . 基于模拟退火粒子群优化算法的裂缝型储层各向异性参数地震反演. 吉林大学学报(地球科学版), 2019, 49(5): 1466-1476

正林 , 永光 . 地震波传播数值模拟. 地球物理学进展, 2004, 19(4): 933-941

军还 . 和极大似然估计——一种新的估计准则. 桂林冶金地质学院学报, 1994, 14(4): 436-444

学明 , 家映 . 地球物理资料非线性反演方法讲座(四) 遗传算法. 工程地球物理学报, 2008, 5(2): 129-140

, 生昌 . 解非线性反演问题的新策略. 石油物探, 1999, 38(3): 76-85

乐洋 . 基于总体最小二乘的大地测量反演理论及应用研究. 测绘学报, 2012, 41(4): 629

文萍 , 庆良 . 利用遗传算法和最小二乘联合反演共和地震位错参数. 地震学报, 1999, 21(3): 285-290

晓青 , , . 四川汶川8级大地震灾害损失快速评估研究. 地震学报, 2009, 31(2): 205-211

志伟 , 力耘 , 同城 . 岩石热弹性理论及其在地球物理中的应用. 地球与行星物理论评, 2021, 52(6): 623-633

扬茂 , 才军 , . 利用断层自动剖分技术的2008年青海大柴旦MW6.3级地震InSAR反演研究. 武汉大学学报(信息科学版), 2012, 37(4): 458-462

吉仓 , 永奇 . 顾及先验信息的模型参数反演——澜沧-耿马地震位错模型参数的求定. 地壳形变与地震, 1997, 17(2): 27-32

世亮 , 玉军 . 2008年汶川大地震同震、震后形变和应力及对松坪沟地质灾害的影响. 地球学报, 2021, 42(5): 691-700

才军 , 扬茂 . 基于InSAR数据的西藏玛尼MS7.9级地震的地壳不均匀性研究. 武汉大学学报(信息科学版), 2008, 33(8): 846-849

清峰 . 基于波动方程的微地震震源位置、震源时间及VTI介质各向异性参数联合反演. 中国地震, 2021, 37(2): 285-299

丙囤 , 海涛 , . InSAR数据约束下的2022年泸定MS6.8地震震源参数及滑动分布. 大地测量与地球动力学, 2023, 43(3): 221-225

冠雨 , 锡武 , 凡顺 . 基于Shearlet变换和广义全变分正则化的地震数据重建. 地球物理学报, 2020, 63(9): 3465-3477.

DOI

佳慧 , 庆华 . 基于混合优化算法的层状介质反演方法. 现代雷达, 2020, 42(5): 46-51

立强 , 海斌 , 天珧 . 基于BP神经网络的波阻抗反演及应用. 地球物理学进展, 2005, 20(1): 34-37

, , 瑞卿 . 一种改进的贝叶斯反演算法. 地球物理学进展, 2020, 35(5): 1911-1918.

DOI

远元 , 家映 . 地球物理资料非线性反演方法讲座(十)——粒子群反演方法. 工程地球物理学报, 2009, 6(4): 385-389

小东 , 从德 , 绪本 . 时间域航空电磁数据的自适应变维贝叶斯反演研究. 地球物理学进展, 2020, 35(5): 2023-2032.

DOI

冬丽 , 夏新 , 正华 . 近场地震动格林函数的解析法与数值法对比研究. 西北地震学报, 2004, 26(3): 199-205

繁昌 , 兴耀 , 国忱 . 用模拟退火神经网络技术进行波阻抗反演. 石油大学学报(自然科学版), 1997, 21(6): 16-18 16-18, 23

繁昌 , 张波 , 兴耀 . 地震数据约束下的贝叶斯随机反演. 石油地球物理勘探, 2014, 49(1): 176-182

志文 , 德辉 , . 二次损失函数下Poisson回归模型中参数的近似Bayes估计. 吉林大学学报(理学版), 2008, 46(5): 836-840

, 宏生 , . 汶川地震强余震(MS≥5.6)的震源机制解及其与发震构造的关系. 中国科学(D辑: 地球科学), 2009, 39(4): 413-426

龙泉 , 福田 , 晓非 . 三维介质中速度结构和界面的联合成像. 地球物理学报, 2006, 49(4): 1062-1067

绍民 , 建新 , 欢乐 . 基于非负最小二乘法的一维MT正则化反演研究. 工程地球物理学报, 2017, 14(3): 253-261

Outlines

/