地震学全波形反演进展

祝贺君，刘沁雅，杨继东

1 美国得克萨斯州立大学 达拉斯分校 地球科学系，美国得克萨斯州理查德森 75080

2 加拿大多伦多大学 物理与地球科学系，加拿大安大略省多伦多市 M5H 2N2

3 中国石油大学（华东）地球科学与技术学院，青岛 266580

0 引言

研究地球内部物质结构及其与板块运动的相互作用是地球科学的一个重要研究内容.目前由于钻探技术的限制，无法直接测量大部分地球深部的物理参数，包括地震波传播速度、电导率、温度和岩石成分等.所以我们需要使用间接的物理和化学方法来探测地球甚至于其他行星的内部构造.由于地震波传播的穿透性，地震成像给我们提供了一个很好的间接探测工具（Dziewonski and Romanowicz,2015;Liu and Gu,2012;Nolet,2011;Rawlinson et al.,2010;Ritsema and Lekic,2020;Romanowicz,2003;Tromp,2020）.从1970 年代开始，地球物理学家们就尝试使用地表记录的地震信号来推测地球内部的温度以及物质组成的三维分布，以用于研究地幔对流、岩石圈演化和板块俯冲等重要的地球科学问题（Aki et al.,1977;Dziewonski and Woodhouse,1987;Grand et al.,1997;Masters et al.,1996;Su et al.,1994;Van der Hilst et al.,1997;Woodhouse and Dziewonski,1984;Zhao,2004）.这种使用地表记录的地震波信号来反演地球内部三维结构的方法被称为地震层析成像，该方法与医学和声学成像有很多的共同点.它能够为地球科学家们提供关于地球内部的地震波速度、各向异性、密度乃至于地震波衰减的三维分布信息.通过研究这些三维模型以及结合岩石物理实验和地球动力学模拟结果，我们可以进一步推测地球内部的温度、部分熔融、物质组成以及水和挥发成分的分布.这些结果能够进一步地帮助我们了解地球内部的结构和演化，进而更好地了解地幔对流和板块运动.

在天然地震学中，所使用的地表记录大多来自于天然地震，又被称为被动震源.也可以使用人工震源激发的弹性波或者声波来进行成像，比如能源勘探中所使用的气枪和可控震源车、炸药以及核爆等，这些震源被称之为主动源.与主动震源和医学成像不同，在被动源地震学研究当中，由于地质构造以及研究区域的限制，无法自由地选择天然震源以及接收台站的位置.这就直接导致了在天然地震成像中，很难获得均匀的成像覆盖率.比如，目前全球地震台网以及区域台网主要分布在北半球的陆地上，所以相对而言，在南半球以及海洋区域的成像分辨率还比较低，这对于全球尺度的地幔成像来说是一个很大的问题.目前正在广泛应用和发展中的海底地震仪和海洋流动地震仪（Nolet et al.,2019），以及利用海底光纤电缆提取震动和应变信号（Lindsey et al.,2019;Zhan et al.,2021），在未来会很好地帮助我们解决这一问题.

早期的地震层析成像主要是基于高频射线理论和采用P 波或者S 波的走时信息.通过运动学和动力学射线追踪，可以计算地震体波在假设的一维或三维地球模型中的走时和振幅（Cerveny,2001）.通过测量实际地震信号的走时和振幅，并且与模拟结果进行对比，试图找到一个最优的三维地球模型来解释所观测得到的地震信号.该方法可以帮助我们研究地壳、地幔以及地核的地震学结构.另外，相对于体波到时，拟合面波频散能够帮助我们更好地研究地壳、岩石圈以及上地幔的结构（Barmin et al.,2001;Ekstrom et al.,1997）.另一种被用于地震成像的信号是地球的自由振荡.通过拟合观测和模拟的自由振荡频率，可以推测地幔乃至于地核的地震学参数（Dziewonski,1971;Masters et al.,1982;Ishii and Tromp,1999）.在早期的走时层析成像方法当中，通常使用人工或者自动相位识别方法，比如长时平均/短时平均（STA/LTA）算法（Allen,1978），去测量实际的地震震相到时.这些测量方法主要提取的是高频的到时信息，所以与高频射线理论相匹配.目前，通过人工智能的方法也可以很好地提取高频到时信息（Ross et al.,2018）.

随着宽频带地震仪精度和稳定性的提高，地震学家们开始使用实际观察到的地震波波形和理论计算所得到的地震图之间的互相关来测量特定震相的到时差.该方法所得到的到时差与信号的频率有关，而射线理论是基于无限高频假设，所以出现了理论与实际观测之间的不匹配.从2000 年左右开始，Dahlen 等提出有限频率层析成像（finite-frequency tomography）的概念（Dahlen et al.,2000）.这一理论结合了射线理论和散射波的波恩近似，相关的结论显示对于一维模型，互相关所提取的到时只受到传统射线周围的速度结构的影响，反而与射线路径本身上的速度扰动无关，并且所影响的范围与菲涅尔带（Fresnel zone）的宽度有关（Hung et al.,2000;Zhou et al.,2004）.几乎同时，Zhao 等（2000）使用自由振荡叠加理论，同样得到了有限频率敏感核.相对于前者，使用自由振荡叠加的计算量较大，但是由于使用了波动方程的严格解，所以它适用于一切震相，包括临界反射/折射波，核幔边界衍射以及面波等.这一有限频率敏感核（finite-frequency sensitivity kernel）也被形象地称为香蕉—甜圈敏感核（banana-doughnut kernel）.该有限频率敏感核的概念把地震波散射融入到层析成像当中，并考虑了波前弥合（wavefront healing）的问题（Nolet and Dahlen,2000）.在这一理论框架下，成像结果与所使用的地震波的波长有关.从理论上而言，它能够提高成像的分辨率和精确度（Montelli et al.,2004a），尤其是能够更好地研究低速带分布，比如地幔柱（mantle plume）等.值得注意的是，在早期的有限频成像当中，仍然使用高频射线或者自由振荡叠加理论来计算一维速度模型下的理论地震图以及敏感核.另外，相关的敏感核在勘探地震学中也有所讨论和应用（Luo and Schuster,1991;Woodward,1992），但是直到2000 年之后这两个方向的理论才得以统一.

另一方面，在勘探地震学领域，Lailly（1984）及Tarantola（1984）革命性地提出在声波介质当中，通过直接拟合观测和理论计算所得到的地震图来反演地球内部的波速结构.这一理论随后被拓展到二维弹性/黏弹性介质反演当中（Gauthier et al.,1986;Mora,1987;Tarantola,1986,1988）.在这一全波形反演（full waveform inversion）理论框架下，Tarantola 使用伴随方法（adjoint method）来计算目标函数的梯度，然后通过局部最优化算法，比如最速下降法或者共轭梯度法（Fletcher and Reeves,1964），来估计局部最优的三维速度模型，进而拟合实际所观察到的波形.值得注意的是，在使用伴随方法来数值计算目标函数的梯度时，需要使用逆时间传播来计算伴随波场.这一过程与勘探地震学中所使用的逆时偏移成像（reverse time migration）是相似的（Baysal et al.,1983;Claerbout,1971;McMechan,1983）.所不同的是在逆时偏移成像当中，需要计算正传以及反射波场，其中后者使用的是地表记录所得到的首次反射波（primary reflections）波形信号，然后使用互相关成像条件（cross correlation imaging condition）来得到地壳中的反射面位置和振幅信息.在全波形反演当中，同样需要计算正传波场.然而在计算反传（伴随）波场中，使用的是实际记录和理论计算所得到的波形残差来构建伴随波场的激发震源.可以使用多种震相，包括直达波、反射波和面波等，来构建地球的三维结构模型，而不仅仅像逆时偏移成像一样来寻找地下反射界面.另外，传统的逆时偏移成像不需要使用非线性迭代算法来求解一个局部最优化问题，它可以被看成是全波形反演的首次迭代结果（Tarantola,1984）.

Tromp 等（2005）将以上的有限频成像、伴随方法、逆时偏移和香蕉—甜甜圈（banana-doughnut）敏感核统一到一个伴随层析成像的框架中.在这一框架下，可以自由地选用各种不同的反演目标函数，比如全波形残差、有限频到时或者振幅等.在伴随成像中，可以使用相同的算法去构建这些目标函数的梯度.对于不同的目标函数，只需要修改相应的激发伴随波场的震源时间函数即可.同时，Tromp等（2005）也提出了使用伴随成像的框架来反演地震波衰减、各向异性以及震源参数（包括地震位置以及震源机制解）.在伴随成像的每次迭代当中，对于每一个主动或者被动震源，只需要计算两次地震波场以获得目标函数的梯度.一次是用原有的震源信息来计算正传波场，另一次则是使用与目标函数有关的伴随震源来计算反传的伴随波场（Liu and Tromp,2006;Tape et al.,2007），所以它的计算量与震源的数目成正比，而与所使用的台站和震相的数量无关.这种方法可以帮助我们有效地计算目标函数的梯度.但是仍然无法得到相应的二阶导数，也就是Hessian 矩阵的信息（Plessix,2006）.另一种相似的方法则是由Zhao 等（2005）提出的散射积分方法（scattering-integral）.利用地震波动方程格林函数所具有的源与观测点之间的互易性，通过使用硬盘存储接收点格林函数，不仅可以计算每个源和台站对之间的敏感核，而且能够同时构建目标函数的梯度及其二阶导数.这样能够帮助我们分析反演结果的分辨率和不确定性信息.但是由于需要存储四维地震波场库，而且每一次迭代时都需要重新计算所有接收点格林函数，所以这种方法相对于伴随成像需要更大的存储空间（Chen et al.,2007a）.而在有些特定情况下，其计算量可能与伴随成像接近.另外，由于伴随层析成像和全波形反演具有许多相似性.所以在本文中，我们有时会相互替代使用这两种名称.全波形反演理论上需要使用记录的全部波形信号.但是考虑到实际记录当中的噪声成分以及反演当中的非线性，很少有研究将全部记录所得的波形直接应用于反演当中.

自从全波形反演理论的提出，这种新方法已经被应用到各种尺度的地球内部结构成像研究当中.比如，Tape 等和Chen 等使用伴随/全波形成像来研究南加州地壳的三维速度结构（Chen et al.,2007b;Lee et al.,2014;Tape et al.,2009,2010）.他们的模型能够帮助我们很好地研究这一区域的地质构造，并且已被美国南加州地震中心（Southern California Earthquake Center）作为区域标准的三维公共速度模型（比如最新的CVM-H15.1.0 和CVMS4.26）.同时这些高分辨率高精度的三维速度模型能够帮助我们更好地研究南加州地震的震源信息（Lee et al.,2011;Liu et al.,2004;Wang and Zhan,2020;Zhao et al.,2006）.另一方面，全波形反演方法也被广泛地应用到区域上地幔成像，包括欧洲（Fichtner et al.,2013;Fichtner and Villasenor,2015;Zhu et al.,2012,2015;）、美洲（Krischer et al.,2018;Zhu et al.,2017,2020b）、大西洋（Colli et al.,2013）、亚洲（Chen et al.,2015;Simutė et al.,2016;Tao et al.,2018）、澳大利亚（Fichtner et al.,2009,2010）、新西兰（Chow et al.,2022）、南极洲（Lioyd et al.,2019）等.目前也被应用于全球地幔研究 当中（Bozdağ et al.,2016;French et al.,2013;French and Romanowicz,2015;Lei et al.,2020）.同时，全波形反演也在油气资源开发过程中发挥越来越重要的作用（Operto et al.,2015;杨积忠等,2014）.我们将会在2.8—2.10 节中详细讨论这些具体的例子.迄今为止，已经有一些全波形反演的开源软件可以供研究人员下载和使用，比如Seisflow（https://github.com/rmodrak/seisflows）、LASIF（https://dirkphilip.github.io/LASIF_2.0/）、PySIT（https://pysit.readthedocs.io/）和IFOS3D（https://git.scc.kit.edu/GPIAG-Software/IFOS3D/）等.

1 理论框架

本文只是回顾全波形反演的主要理论框架和计算过程，具体的推导可以参考其他文献和专著（Chen and Lee,2015;Fichtner,2010;蒋梦凡等,2021;Liu and Gu,2012;Tromp et al.,2005;Virieux and Operto,2009;杨午阳等,2013;杨勤勇等,2014）.

1.1 目标函数

全波形反演在数学上是求解一个非线性最优化问题（Liu and Gu,2012;Plessix,2006;Tromp et al.,2005）.由于所涉及的模型参数量以及正演数值计算量十分巨大，目前在实际应用当中，只能通过使用目标函数的导数与局部最优化算法来求解这一问题.常用的方法有共轭梯度（Fletcher and Reeves,1964）或者L-BFGS 算法（Nocedal and Wright,2006）.在目前阶段尚无法使用全局优化算法（Sen and Stoffa,2013;Tarantola,2005），比如模拟退火法或者蒙特卡罗法，来解决大型的三维反演问题.这一局限性致使目前大部分研究只能给出某一个问题的局部最优解，而无法准确地提供关于这一问题的后验概率分布（posteriori probability density distribution）以及不确定性分析（uncertainty quantification），这仍然是当前研究的前沿问题（详见2.5 节）.

由于以上的问题，所以选取一个最佳的目标函数/残差函数来求解这一局部最优化问题至关重要.这一目标函数是我们想要最小化的标量，它也被用来衡量最优化过程是否收敛.地震学家们主要使用的观测数据是地表位移记录，也就是地震图.它本质上是记录地表震动的时间和空间的函数.目前大部分区域尚没有密集的台站，只有少数区域具有覆盖面积大且分布密集的台网，比如USArray 和Hinet 等.所以不失一般性，我们假设所采集的地震图是一维的时间函数.使用目标函数来衡量实际观测与理论模拟的地震记录之间的差别.在伴随层析成像和全波形反演中，最直接也是最早使用的目标函数是波形之间的最小二乘残差，即：

式中，目标函数J是一个标量.在公式（1）中它的单位是 m2.sw和dw分别表示模拟和观测所得到的某个分量w上的地震图.e和r分别表示某一震源和台站.Ns、Nr、Nw分别表示地震、台站和时间窗的数目.由于观测所得到的地震图中难免有噪声，所以我们可以使用特定的时间窗来测量某一个时间段中的波形差（Maggi et al.,2009）.时间窗的选取可以根据波形之间的相似度以及所需要研究的特定震相来决定.其中值得注意的是，时间窗中的信号可以是单一的传统震相，也可以是多个不同路径的震相叠加.这一优势使得全波形反演和伴随成像能够更加充分地利用地震记录信息来约束地球内部结构.由于陆地地震台可以记录三分量地表震动，所以目标函数可以根据三分量位移来计算求得.

全波形目标函数的好处在于它能够最直接地衡量观测与理论计算之间的细微差别.所以它能够帮助我们更好地约束地球模型，以及提供更高分辨率的反演结果（Virieux and Operto,2009）.但是它的缺点包括：（1）如果记录地震图的信噪比较低，则反演过程主要是去拟合不必要的噪声，这些噪声最终会被投射到反演模型当中，导致错误的结果和解释；（2）如果观测和模拟波形之间的到时差大于信号的半周期，则会出现所谓的周期跳跃（cycle skipping）现象（Virieux and Operto,2009）.从而使得全波形目标函数中出现许多局部最小值（董良国等,2013;Engquist and Yang,2019;Métivier et al.,2016;Yang et al.,2018）.因为目前我们大多数时候使用的是目标函数的梯度来求解这一局部最优化问题，所以如果起始模型远离全局最优解，只能获得起始模型周围的局部最优解，而非全局最优解.这一问题在高频全波形反演中尤为严重.

考虑到以上全波形目标函数所存在的诸多问题，另一种被广泛使用的目标函数是基于到时信息.这一点与传统的P 波或S 波射线走时层析成像的概念相同.可以设计如下的目标函数：

式中， ΔTw(e,r;m)是观测以及预测的某个震相的到时差.相对于全波形最小二乘目标函数[公式（1）]，使用到时残差目标函数可以更好地帮助我们在反演过程中避免陷入局部最小值（Luo and Schuster,1991）.

另一类目标函数被称之为双差函数，它使用两个台站接收到的同一震源的到时差来反演地球结构.这与传统的双差层析成像（Zhang and Thurber,2003）和地震定位（Waldhauser and Ellsworth,2000）的概念是相似的，该方法目前也被用于伴随层析成像当中（Yuan et al.,2016）.它可以减少震源对地球结构反演的影响，同时主要反映台站之间的速度结构.

1.2 模型参数化

在设立了目标函数之后，第二个问题就是要确定我们所要求解问题的模型参数.它涉及到以下几个方面的问题：

（1）如何离散化所关心的研究目标.其中最直接的方法就是使用正演模拟当中的数值网格来参数化连续性的地球模型，这也是目前大多数研究中所采用的方法.它的优点是正演和反演使用的是同一套网格点，易于操作.缺点则是在正演模拟当中，为了满足数值计算的精度，稳定性以及频散要求（Komatitsch and Tromp,1999;Virieux et al.,2011），网格点的数目一般十分巨大.比如，在地壳或地幔反演当中，所使用的正演网格点可以轻松到达百万级别（Tape et al.,2010;Zhu et al.,2015）.另一种方法就是用全球层析成像中所广泛使用的基函数来离散化所要研究的目标（Nolet,2011;Rawlinson et al.,2010）.常用的基函数包括水平方向上的球谐函数以及深度方向上的样条函数（Dziewonski and Romanowicz,2015），使用这些基函数可以大大地减少模型参数量.另一个优点则是当地球模型本身具有有限的波数分布时，这样也可以大大地减少模型参数量.这点在全球地幔成像当中尤为突出，比如，地幔的主要结构用6～8 阶的球谐函数就可以很好地描述（Su and Dziewonski,1992）.另外使用基函数还可以起到对反演模型进行平滑或者是低通滤波的效果.

（2）传统地震成像最关心的地球模型参数是各向同性的P 波和S 波速度（Grand et al.,1997;Van der Hilst et al.,1997;Zhao,2004）.通过反演求得的这些波速信息可以帮助我们分析地球内部的温度、物质组成以及部分熔融和挥发成分等.我们可以通过各种方法来构造目标函数与模型参数之间的关系如下：

式中，α、β 和ρ 分别表示P 波和S 波的速度以及密度.Kα、Kβ和Kρ分别是目标函数相对于以上三个模型参数的梯度.如果我们使用到时差[公式（2）]作为目标函数，同时用高频射线或者有限频方法来计算这些导数，则以上的公式与传统的射线走时层析成像是一致的.与全波形反演所不同的是，我们使用三维地球模型以及数值伴随算法来计算上述目标函数的梯度（见1.3 节）.

使用伴随方法的另一个好处是可以考虑多个不同的模型参数，包括各向异性波速度以及地震波衰减.即可以在一个统一的理论框架下来求解地球的多参数模型（Fichtner and Driel,2014;Sieminski et al.,2007a,2007b;Tromp et al.,2005）.这也是目前的一个研究重点，本文将会在2.2、2.3 和2.6 节中具体讨论.

1.3 使用伴随方法计算目标函数的梯度

从公式（3）中可以看到，伴随成像/全波形反演与高频射线和有限频成像的一个重要区别在于如何计算目标函数的梯度.在伴随层析成像和全波形反演当中，使用的是伴随方法来计算这一三维空间函数（Plessix,2006）.根据不同的目标函数，这些梯度的单位也有所不同.如果目标函数是全波形最小二乘残差，梯度的单位是而如果我们选取到时残差作为目标函数，那么它们的单位则是同时值得注意的是，这里的目标函数的梯度与前面所描述的敏感核是不一样的，后者不包括实际观测的数据和残差.使用伴随方法来推导目标函数的梯度大体上可以分为以下两类：一种是基于地震波场的波恩近似展开来推导（Tarantola,1984;Tromp et al.,2005）.另一种则是使用拉格朗日乘数（Lagrange multiplier）的方法来推导（Liu and Tromp,2006;Virieux and Operto,2009）.以上两种方法所得到的结果是相同的.具体的数学推导过程可以参考相关的文献和专著（Chen and Lee,2015;Fichtner,2010）.

基于Tromp 等（2005）的推导，可以使用正传和伴随波场的卷积来构建目标函数的梯度，其中伴随波场可以通过计算反传波场得到（如图1 所示）.用包含测量信息的伴随时间函数在接收台站上作为震源来激发伴随波场.所以我们可以使用正演的数值算法和程序来求解伴随波场.对于二阶波动方程，正传和伴随波场之间的差别只在于波传播的时间方向以及震源时间函数.根据不同的模型参数，可以使用正传和伴随波场的时间和空间导数来计算它们的梯度.如果选取密度ρ、体积模量 κ和剪切模量 µ作为模型参数，它们的梯度可以通过以下的公式求得：

图1 二维均匀速度模型下构建SH 波的敏感核.从左到右分别是正传波场、伴随波场、相互作用波场和剪切波速度的敏感核.五角星和方块分别表示震源和接收台站的位置.时间从下往上分别是8 s、16 s、24 s、32 s 和44 s.地表使用的是自由表面边界条件（修改自Tromp et al.,2005 中的图3）Fig.1 Construction of an SH sensitivity kernel in a 2D homogeneous velocity model.From left to right are forward wavefield,adjoint wavefield,interaction wavefield and shear wave velocity sensitivity kernel.The star and rectangular represent source and receiver.The time steps from bottom to top are 8,16,24,32 and 44 seconds.A traction free boundary condition is applied to the Earth's surface (modified from Figure 3 in Tromp et al.,2005)

式中，s和s+分别是正传和伴随波场，D和D+则是正传和伴随应变偏量（strain deviator），T是记录时间.

如果我们选取波阻抗Zα，以及P 波和S 波速度作为模型参数，则可以通过以下的公式来计算它们的梯度：

值得注意的两点是：（1）在有的文献中可以看到使用的是卷积[如公式（4）]，而在另一些论文中使用的则是互相关（主要在勘探地震学中）.它们本质上是一样的，区别只在于对伴随波场的时间传播方向的定义不同.（2）公式（4）中的密度梯度与勘探地震学中的逆时偏移成像非常相似，这一点在Tarantola（1984）中已经指出.Claerbout（1971）得到的成像条件是基于物理的直观考虑，就是地下反射面存在于下行（正传）波场与反射（伴随）波场相位一致的时候.在这一点上数学推导与物理直观达到了高度统一.只不过Claerbout 的成像条件是目标函数相对于密度的梯度.而后通过数值计算，Zhu 等（2009）发现相对于密度梯度而言，波阻抗梯度能够更好地反映模型的高波数信息（Wu and Aki,1985），也就是反射率.这一点已被工业界用于实际勘探开发当中，被称为逆反射成像（inverse scattering）（Douma et al.,2010）或者是能量模式成像条件（Rocha et al.,2016）.

1.4 最优化迭代

通过以上的伴随方法可以求得相关模型的目标函数梯度.然后通过基于导数的局部最优算法，比如共轭梯度或者是L-BFGS 方法，可以获得初始模型附近的局部最优解（Akcelik et al.,2002,2003）.其中，共轭梯度法使用的是本次迭代的梯度以及上次迭代的方向之间的加权平均来求得本次迭代的方向.相对于最简单的最速下降法，它具有更好的收敛性.另一类方法被称为拟牛顿迭代法，其中一个典型代表是L-BFGS 方法（Nocedal and Wright,2006）.由于使用了前几次迭代的方向信息，LBFGS 方法能够帮助我们进一步构建近似的二阶导数信息.所以相对于共轭梯度和最速下降算法，它具有更好的收敛性（图2），进而加快整个反演的进度（Modrak and Tromp,2016）.这一点在实际应用中至关重要，因为相对于射线和有限频成像，伴随层析成像和全波形反演的计算量要求仍然较高.所以如果能够减少迭代的次数，则能更加有效地反演地下模型，进而更好地研究相关的科学问题.

图2 对比九个不同理论模型下使用三种迭代方法所得到的收敛性，包括最速下降法（黑线），非线性共轭梯度法（绿线）和L-BFGS 法（红线）.不同的理论模型如每个图的标题所示（修改自 Modrak and Tromp,2016 中的图2）Fig.2 Convergency comparison of three iterative methods for nine synthetic velocity models,including the steepest descent method(black lines),nonlinear conjugate gradient method (green lines) and L-BFGS method (red lines).The 2D synthetic velocity models are shown in the titles of each panel (modified from Figure 2 in Modrak and Tromp,2016)

另外两种可以加速迭代收敛的方法是：（1）使用近似方法来构建Hessian 矩阵.因为Hessian 矩阵本身十分巨大，它的大小与模型的网格点数的平方成正比（Pratt et al.,1998）.所以目前仍然很难直接地计算这一矩阵.但是通过各种近似，可以构建它的对角量，使用该信息可以帮助我们更好地平衡模型浅部与深部的梯度大小，并且可以减少在震源以及接收台站上模型梯度数值较大的问题.（2）使用高斯-牛顿方法来迭代地考虑Hessian 矩阵的作用（Epanomeritakis et al.,2008;刘璐等,2013;Métivier et al.,2014），这被称之为截断高斯牛顿（Truncated Gauss-Newton）方法.在这一方法中，在外部共轭梯度或者L-BFGS 循环当中加入一个内循环，通过求解法方程式（Normal equation）来迭代地考虑Hessian 矩阵的作用.这一方法已被用于勘探地震学的数值试验当中，但是相对于一般的共轭梯度或者L-BFGS 方法，它涉及到内外两个迭代循环，所以计算量仍然十分巨大.

通过以上的迭代，如果目标函数的数值不断减小，而且迭代模型趋于稳定，则在一定迭代次数之后，我们可以得到一个最终优化的模型.另外，在迭代过程中所广泛使用的技巧是多尺度反演（Bunks et al.,1995;王毓玮等,2016;张文生等,2015）.也就是在迭代的早期，我们使用长周期的信号来反演大尺度的结构.随着反演的推进，我们逐渐地加入短周期的信号来进一步约束小尺度的结构.使用这一方法能够更好地帮助我们规避目标函数在初始模型周围的局部最小值以及反演当中的强非线性问题.

2 研究进展

2.1 目标函数的设计

如1.1 节所述，设计一个好的目标函数对于全波形反演是否成功具有十分重要的意义.而波形最小二乘残差具有前述的诸多缺点，所以如何设计一个好的目标函数是目前的一个研究热点.本文总结几个大的目标函数类型.

（1）直接基于时间信号的特征.这一类型包括前面的波形最小二乘差（Tarantola,1984）以及后来提出的波形包络差（胡勇等,2017;Wu et al.,2014）.相关的研究显示使用波形包络可以很好地增加信号中的低频信息，进而减小反演问题的非线性以及避免陷入局部最小区域；（2）从地震图中提取到时和频率的信息.这一类包括使用互相关来提取观测与预测的时间差（Leeuwen and Mulder,2010;Luo et al.,2016）.可以直接在时间域来测量，也可以通过小波变换，进而在时间—频率域中测量（Fichtner et al.,2008;Yuan and Simons,2014）.该方法也可以用来测量振幅差；（3）自适应反演.Warner 等提出使用维纳滤波（Wiener filter）的概念来测量两个波形之间的差（Warner and Guasch,2016;Zhu and Fomel,2016），这一方法可以不使用时间窗来分割地震信号.相关的数值试验显示很好的反演效果.其他类似的方法包括动态拟合（dynamic warping）（Hale,2013;Ma and Hale,2013）和地震图匹配（seismogram registration）（Baek et al.,2013;Zhu,2018a）；（4）最优传输方法.这一方法是基于数学中的最优传输路径来测量两个时间序列之间的差别（Engquist and Yang,2019;Métivier et al.,2016;Yang et al.,2018）.相对于最小二乘波形差，这一目标函数具有更宽的收敛带，从而对初始模型的要求更低（图3）.对于一维时间序列，有很好的数学方法来求解这一最优化问题，进而得到相关的伴随震源信号.相关的研究也被应用于实际地壳反演当中（Dong et al.,2022）.但是目前对于二维图像的对比，需要求解Monge-Ampere 方程（Engquist and Froese,2014）.这是一个非线性方程组，它的稳定数值解还有待进一步研究.同时这一方法要求两个信号之间具有相同的能量，这在实际的地震记录中是很难得到满足的.其他在勘探地震学中开发的用于解决周期跳跃的方法还有时间轴扩展（Biondi and Almomin,2014）和波场重构（van Leeuwen and Herrmann,2013）等.

图3 对比不同目标函数的特征.（a）显示的是子波波形.（b）和（c）分别显示的是基于最小二乘波形残差和最优化路径所得到的目标函数随着不同时间移动（s）的特征（修改自Engquist and Froese,2014 中的图1）Fig.3 Comparison of misfit functions based on least-square waveform differences and optimal transport distances.Panel (a) shows the input source wavelet.Panels (b) and(c) are the misfits as functions of time shifts for the leastsquares waveform differences and optimal distances,respectively (modified from Figure 1 in Engquist and Froese,2014)

2.2 各向异性反演

我们使用各向异性来描述地震波沿不同方向所具有的不同传播速度，在地震学研究中包括径向各向异性和方位各向异性（Backus,1965;Hess,1964;Montagner and Nataf,1986;Smith and Dahlen,1973）.这些各向异性来自于形状或是晶格的定向排列[shape-preferred orientation (SPO)/lattice-preferred orientation (LPO)]（Jung and Karato,2001;Karato et al.,2008;Nicolas and Christensen,1987;Zhang and Karato,1995）.地壳和地幔中的各向异性矿物包括橄榄石、辉石以及云母，它们是导致地震波各向异性的主要原因.目前可以使用剪切波分裂（Silver and Chan,1991;Silver,1996）以及面波频散（Lin et al.,2011;Marone and Romanowicz,2007;Moschetti et al.,2010）来研究地球内部的各向异性分布.相对来说，如果有很好的台站分布，横波分裂可以提供很好的横向分辨率.但是由于使用的是垂向入射的远震体波信号，所以对于垂向上的各向异性分布具有很差的分辨率.相反，因为不同频率的面波对不同深度的速度具有不同的敏感性,使用面波频散可以帮助我们提高垂向上的分辨率.但是由于在面波相速度或群速度反演当中，使用最小二乘法来求解一个层析成像问题，通常需要加入正则化约束.同时由于面波的横向传播，相对来说，它的横向分辨率要低于垂向分辨率.通过研究地球内部的各向异性，可以讨论地壳、岩石圈与软流圈中的应力、形变以及流动性特征，同时也可以用来研究地幔中的对流以及地核构造（Long and Silver,2008;Long and Becker,2010;Park and Levin,2002）.

如前所述，伴随层析成像和全波形反演提供了一个统一的框架来反演各种不同模型参数的组合.所以它也可以用于反演各向异性（Rusmanugroho et al.,2017;Sieminski et al.,2007a,2007b）.比如在地壳成像方面，Wang 等（2020）在Tape 等的南加州地壳模型基础之上，加入了径向各向异性.他们发现在圣安德烈斯断层西面的中下地壳中，横向偏振的剪切波速度比纵向偏振要低（负径向各向异性～-6%）.而跨过断层，这一比例正好相反.这与上地壳和上地幔中常见的径向各向异性特征正好相反.为了解释这一现象，他们提出一种可能性是由于下地壳的斜长石具有与黑云母和橄榄石不一样的各向异性特征（Wang et al.,2020），导致它们的快速轴与应变方向垂直而非平行.在地幔成像方面，Zhu 等（2020c）使用伴随方法反演上地幔的方位各向异性，进而研究俯冲板块与地幔对流的关系.例如，他们发现在中南美洲和卡斯凯迪亚俯冲带中，当出现板块断裂时，上地幔物质会流动通过这些板块的间断空间（图4），从而在板块窗以及板块周围产生环流（Hall et al.,2000;Russo and Silver,1994）.这些观测与岩石物理实验（Buttles and Olson,1998;Kincaid and Griffiths,2003）和地球动力学模拟相吻合（Schellart,2004;Schellart et al.,2007;Stegman et al.,2006）.另外，在欧洲上地幔反演结果当中（Zhu and Tromp,2013），发现反演所得到的方位各向异性方向与过去三千万年的大地构造过程有很好的相关性.同时通过分析方位各向异性的大小在垂向上的变化，他们推测出脆性和韧性变形在岩石圈中的分布，这一结果与实验和震源深度分布结果相吻合（Chen and Molnar,1983;Kohlstedt et al.,1995）.另外，通过同时反演Rayleigh 波和Love 波，Fichtner 等（2010）发现澳大利亚大陆和海洋岩石圈在径向各向异性上存在差别.比如在海洋岩石圈中，正径向各向异性（横向偏振剪切波快于纵向偏振剪切波）主要存在于80～150 km 之间.而在大陆岩石圈中，径向各向异性比较弱且主要出现在小于100 km 左右的区域.这与早期的全球面波成像结果相吻合（Nettles and Dziewonski,2008）.

图4 三维俯冲板块导致的地幔对流.（a）和（b）分别显示中美洲和卡斯凯迪亚俯冲带的反演结果.绿色块体显示的是剪切波速度扰动大于1.5%的结果.黄色箭头表示通过伴随方位各向异性成像所得到的地幔对流情况[（a）和（b）分别修改自 Zhu et al.,2020a 中的图7 和 Zhu et al.,2020c 中的图5 ]Fig.4 Three dimensional subducting slabs and induced mantle flows.Panels (a) and (b) are results for the Middle American and Cascadian subduction zones.Green bodies represent regions with shear wave velocity perturbations greater than 1.5%.Yellow arrows represent horizontal mantle flows constrained by azimuthal anisotropy adjoint tomography [Panels (a) and (b) are modified from Figure 7 in Zhu et al.,2020a and Figure 5 in Zhu et al.,2020c,respectively]

2.3 地震波衰减成像

地震波在地球内部传播过程中经受衰减作用，其中动能和弹性势能被转换成热能.相对于过去30多年地震波波速层析成像的显著进步来说（Grand et al.,1997;Meschede and Romanowicz,2015;Schaeffer and Lebedev,2013），地震波衰减成像的进步相对较小（Dalton et al.,2008;Gung and Romanowicz,2004;Romanowicz,1995）.相对于以前的结果，在最近的一些大陆和全球尺度的衰减成像当中（图5），能够看到相对较好的可比性（Adenis et al.,2017a,2017b;Bao et al.,2016a,2016b;Ma et al.,2016;Karaoglu and Romanowicz,2017）.相对于波速成像，地震波衰减成像有以下的几个难点：（1）需要有高质量的地震波振幅信息，这就需要具有较高信噪比的地震记录，同时要对地震仪器响应进行校正；（2）体波走时和面波频散主要受地震波波速的影响，然而地震波振幅受到多重因素的影响，其中包括介质衰减、波速、震源大小、辐射模式、弹性聚焦与散焦（Durek et al.,1993;Romanowicz,1987;Ruan and Zhou,2010,2012）、散射以及台站局部的放大效应等（Eddy and Ekstrom,2014;Lin et al.,2012）；（3）目前尚难以区分内在衰减和由地震波散射所导致的表象衰减.所以在衰减成像当中，需要仔细地考虑这些因素.

目前，全波形反演也被用于研究地球内部地震波衰减的分布（Karaoglu and Romanowicz,2017,2018;Zhu et al.,2013）.Tromp 等揭示了在伴随成像/全波形反演框架下，如果使用标准线性固体模型（standard linear solid model）（Liu et al.,1976）来表示一定频段下与频率无关的品质因子Q，那么Q的目标函数的梯度计算公式与剪切模量的公式是一样的.只需要把计算伴随波场的震源时间函数改成新的计算公式即可（Tromp et al.,2005）.所以，在这一框架下，相对于波速成像，需要使用两个不同的伴随震源时间函数来计算伴随波场，从而同时得到目标函数对于波速和衰减的梯度信息，这也就使得反演的计算量翻倍.这一方法被Zhu 等（2013）应用到欧洲地幔反演当中.他们发现在上地幔浅部，品质因子Q模型与地震波波速分布具有较好的相关性，这一点与全球尺度的地震波衰减成像结果相吻合（图5）.比如，洋中脊地区具有较高的衰减和低速度带，而在大陆克拉通下则是低衰减和高波速带.另外，在北大西洋的地幔转换带当中，他们发现Q值比周围较小，表示这一区域的衰减较大.为了解释这一现象，他们认为这一区域应受含水量的影响较大，这与地幔转换带水过滤的科学假设相吻合（Bercovici and Karato,2003;Huang et al.,2005;Schmandt et al.,2014）.另一方面，Komatitsch 等（2016）改进了伴随波场的存储方法，从而大大地改进了计算衰减敏感核的稳定性和精度问题.

图5 对比三个全球尺度的上地幔地震波衰减成像结果.（a-c）分别来自于QRLW8（Gung and Romanowicz,2004）、SEMUCB-UMQ（Karaoglu and Romanowicz,2018）和QRFSI12（Dalton et al.,2008）.每一个图的下方第一行和第二行分别显示的是各个深度的剪切模量衰减值及其扰动量，其中SEMUCB-UMQ 是全波形衰减反演的结果（修改自Karaoglu and Romanowicz,2018 中的图9 ）Fig.5 Comparison of three global scale upper mantle seismic attenuation tomography models.Panels (a-c) are results from QRLW8(Gung and Romanowicz,2004)，SEMUCB-UMQ (Karaoglu and Romanowicz,2018) and QRFSI12 (Dalton et al.,2008).The first and second lines in each panel represent absolute and relative perturbations of shear modulus attenuation (modified from Figure 9 in Karaoglu and Romanowicz,2018)

标准线性固体模型只是众多表示地震波衰减理论模型当中的一个，它具有以下几个缺点：（1）品质因子Q并不在所推导的地震波公式当中，相反，需要求得相应的应力和应变松弛时间因子（Carcione et al.,1988;Robertsson et al.,1994）.这给反演和正演都带来了一定的困难；（2）如上所述，如果要同时获得波速和衰减的梯度信息，需要计算两次伴随波场，从而加大了反演的计算量.为了解决上述问题，在过去几年当中，一些用以描述地震波在地球内部衰减现象的新的理论模型被提出.比如Fichtner 和Driel（2014）提出一个改进的标准线性固体模型，使得Q直接出现在波场公式中，进而简化了计算目标函数相对于Q的梯度的过程.同时这种方法也可以帮助我们表示与频率相关的Q的情况（Lekic et al.,2009）.另一方面，Zhu 等和Yang 等从不同的角度提出了常量Q的声波方程（Yang and Zhu,2018;Zhu and Harris,2014）.而后这些方程被推广到黏弹性介质（Zhu and Sun,2017），以及品质因子Q的全波形反演过程当中（Yang et al.,2020）.

另外在地震波衰减反演当中，选取好的反演策略和目标函数也十分重要.相关的研究表明使用两个反演阶段具有更好的效果，即第一阶段主要使用到时信息来约束地震波速度，而在第二阶段中加入振幅并同时反演速度和衰减（Kamei and Pratt,2013;Zhu et al.,2013）.另外，相对于波形和到时目标函数，基于振幅的目标函数，比如波形包络差和振幅频谱比，对衰减反演具有更好的效果（Pan and Wang,2020）.

2.4 模型正则化

模型正则化是求解反演问题当中一个必不可少的重要环节，它在射线和有限频成像中都起到平滑模型以及提供先验信息的作用（Aster et al.,2018;Backus and Gilbert,1968,1970;Tarantola,2005）.同样在全波形反演当中，正则化也被考虑进来.因为和所有其他的反演问题一样，我们也要面对问题的非唯一性，数据的不准确性以及采样的不完备性.目前主要有以下几种途径来提供正则化约束：（1）对目标函数相对于模型参数的梯度进行高斯平滑或者低通滤波.这样做的效果与使用Tikhonov正则化是一致的（Operto et al.,2006;Tape et al.,2007）；（2）显性地使用Tikhonov正则化.这样能够帮助我们在拟合实际观测数据的同时，提供一个光滑的或者是具有最小能量/振幅的模型；（3）可以在正则化约束当中提供相关的先验信息，比如Asnaashari 等（2013）将测井所得到的一维速度模型考虑到正则化约束当中，这样能够帮助我们在拟合地表观测地震图的同时也拟合实际得到的测井数据；（4）使用稀疏正则化.相对于Tikhonov 正则化，如果我们的模型在物理空间或者某个特定转换空间下具有系数稀疏性的特征.那么可以在反演过程中构建一个稀疏化约束，从而使用更少的模型数量来拟合实际观测数据（Lin and Huang,2014）.如图6 所示，使用改进的全变差正则化方法所得到的弹性波全波形反演模型能够更好地约束强速度间断面，同时具有提高模型信噪比的作用.另外，可以使用小波变换或seislet 来将模型转换到相应的变换空间，然后用全变差（total variation）来做相关的稀疏性约束（Askan and Bielak,2008;Xue et al.,2017）.

图6 不同正则化约束下二维弹性波全波形反演所得到的P 波速度结果.（a-c）分别表示基于Tikhonov、全变差和改进的全变差正则化的结果（修改自 Lin and Huang,2014 中的图10）Fig.6 P wave velocity models from an elastic full waveform inversion with different regularization schemes.Panels(a-c) are results constrained with Tikhonov,Total variation and modified Total variation regularization,respectively (modified from Figure 10 in Lin and Huang,2014)

2.5 分辨率和不确定性

相对于射线层析成像，全波形反演能够帮助我们通过拟合实际观察和模拟所得到的波形来更好地提取地球内部物质信息.在这一过程当中，我们能够考虑复杂的物理过程，比如波前弥合，有限频和多路径传播等.但是对于伴随层析成像和全波形反演，目前的一大挑战是如何评估所得到的模型的分辨率以及不确定性.由于在全波形反演当中，我们需要很高的计算量以及较大的模型空间，目前我们尚无法得到传统的分辨率矩阵的信息（Backus and Gilbert,1968,1970）.同时因为庞大的计算量，我们也很难像射线层析成像一样进行棋盘（checkerboard）实验（Leveque et al.,1993）.在过去的一段时间里，我们有以下几种方法来分析相关的成像分辨率和不确定性：（1）Fichtner 和Trampert（2011）以及Zhu 等（2016）使用点扩散函数来分析某个区域的分辨率，以及多参数反演下模型参数之间的互相干扰（trade off）.这些点扩散函数是Hessian 矩阵与模型在某个区域的扰动量相乘的结果.比如，图7 显示在一个模拟反演当中，不同的地区由于射线覆盖、背景速度等的影响，所得到的点扩散函数具有比较大的差别.这一分析方法的缺点是计算量比较大，尤其是当我们要分析多个区域的分辨率时；（2）Bui-Thanh 等（2013）在贝叶斯反演框架下重新构建全波形反演.通过使用low-rank 近似，他们能够得到后验协方差矩阵的特征值和特征向量.这样能帮助我们不仅得到某个局部最优解，而且可以得到在这一特定结果下的后验概率分布信息（Zhu et al.,2016）；（3）Fichtner 和van Leeuwen（2015）使用一些随机样本的互相关来得到方向以及位置相关的分辨率信息，以及各个参数之间的互相干扰.他们的结果显示使用多个随机样本就可以提取出相关的分辨率信息，这一方法被应用到欧洲上地幔成像当中；（4）Liu 和Peter（2019）使用square-root variable metric 的方法来取代传统的共轭梯度算法，从而在不需要矩阵参与计算的情况下来构建Hessian矩阵，进而分析后验概率分布的信息；（5）Thurin等（2019）使用卡曼滤波（Karman filtering）来分析全波形反演的不确定性.这些方法都在一定程度上得到全波形反演的分辨率以及不确定信息.

图7 使用点扩散函数分析二维全波形反演所得到的不同位置上的分辨率.（a-i）分别显示点扩散函数在右下角图中不同区域的结果（修改自 Fichtner and van Leeuwen,2015 中的图5）Fig.7 Using point-spread functions to analyze resolution at different locations in an 2D full waveform inversion.Panels (a-i) illustrate results at different locations as shown in the bottom right panel (modified from Figure 5 in Fichtner and van Leeuwen,2015)

2.6 多参数反演

如2.2 和2.3 节所示，能够在全波形反演这一统一的框架下构建不同地震学模型参数的梯度，进而反演它们的三维分布.但是目前大多数的方法只是基于一阶导数的信息，所以就难以分析并且减少各个不同参数之间的互相串扰（孙史磊等，2020）.而要提取这一信息，需要考虑Hessian 矩阵的影响（Pratt et al.,1998）.一种可以帮助我们更好地设计不同参数组合的方法是通过理论分析各个参数组合对于入射地震波的脉冲响应.希望找到一个参数组合，使得各个参数的脉冲响应在方位上的重合尽可能地小（Wu and Aki,1985）.如图8 所示，相对于P 波速度和密度的组合，P 波速度和波阻抗的组合能够更好地区分宽角度和窄角度的辐射响应（Operto et al.,2013）.所以后一组合在多参数反演中具有更好的效果.另外，在多参数反演下，Hessian矩阵的大小由参数的个数以及模型网格点数目所决定，它是一个区块化的对称矩阵.区块之间的数值表示了各个参数对反演的影响，以及它们之间互相干扰的情况.如果能够构建Hessian 矩阵以及计算它的逆矩阵，进而在迭代中考虑这一因素，则可以更好地进行多参数联合反演（刘璐等,2013;Operto et al.,2013）.但是如前2.5 节所述，目前我们仍然很难得到Hessian 矩阵以及它的逆矩阵，所以大部分的研究在于近似Hessian 矩阵，比如：（1）使用理论方法来近似Hessian 矩阵的对角分量，并在共轭梯度迭代当中将他们当作预调节算子（preconditioner）.这样能够在多参数联合反演当中减少各个参数之间的互相干扰，同时能够更好地平衡深部与浅部的反演结果；（2）在迭代过程中使用不同的算法来近似地构建Hessian 矩阵.这些算法包括L-BFGS、truncated Gauss Newton（Métivier et al.,2014）以及前面所提到的square root metric 方法（Liu and Peter,2019）.它们都可以相对有效地减少各个参数之间的互相干扰.

图8 散射波对于不同模型参数组合的辐射样式.（a,b）分别表示散射波对于声波速度和密度组合所得到的辐射样式.（c,d）分别表示声波速度和波阻抗组合的结果.蓝色和白色五角星分别表示激发震源和散射点.绿色曲线表示由射线加上波恩近似所得到的理论波前（修改自 Operto et al.,2013 中的图2）Fig.8 Radiation patterns of scattered waves for different combinations of model parameters.Panels (a,b) show the radiation patterns for P wave velocity and density,respectively.Panels (c,d) are results for the combination of P wave velocity and impedance.Blue and white stars denote exciting source and scattering point.Green curves represent the amplitudes of scattered waves from Ray+Born approximation (modified from Figure 2 in Operto et al.,2013)

2.7 与背景噪声记录的结合

从2000 年开始，背景噪声成像已经被广泛地应用于研究地壳波速（Lin et al.,2008;Shapiro et al.,2005;Yang et al.,2007）、各向异性（Huang et al.,2010;Lin et al.,2011;Moschetti et al.,2010;Yao et al.,2010），以及观测地壳速度随时间变化（Brenguier et al.,2008a,2008b）.在传统的背景噪声成像中，首先计算背景噪声互相关函数（Bensen et al.,2007;Yao et al.,2006），然后提取面波的相速度或者群速度的频散曲线（Barmin et al.,2001）.通过使用一维速度反演（比如马尔科夫链蒙特卡罗MCMC）来拟合这些频散曲线（Shen et al.,2013）.最后通过插值这些一维速度模型来得到地球的三维速度结构.

背景噪声成像数据也被应用于全波形反演当中.最常见的方法是假设通过背景噪声互相关所得到的记录能够很好地近似两个台站之间的格林函数.通过拟合观测的经验格林函数和理论地震图，可以用全波形反演来约束三维速度模型.这一方法已经被应用于研究西藏东南部的地壳速度（Chen et al.,2014）以及卡斯凯迪亚地壳与上地幔结构当中（Gao and Shen,2014）.另外，Lee 等（2014）和Wang 等（2020）结合天然地震记录以及背景噪声互相关结果，分别反演美国南加州地壳速度和各向异性分布.如图9 所示，相对于只使用天然地震记录所得到的模型，加入背景噪声互相关能够更好地约束模型区域的速度扰动.另外，Zhu（2018b）使用Rayleigh 面波互相关来研究美国得克萨斯州北部和俄克拉何马州的地壳速度.相对于天然地震记录，使用背景噪声成像能够研究构造活动比较弱、地震震源分布不均匀区域的地壳波速结构.美国得克萨斯州北部和俄克拉何马州在2008 年以前没有太多的地震，所以相关的研究比较少.新的地壳模型能够更好地研究这一区域的诱发地震的震源特征（Ellsworth,2013）.以上的方法都是基于垂直分量Rayleigh 面波的经验格林函数，最近Wang 等（2019）将这一方法推广到三分量格林函数噪声伴随成像当中，并在南加州各向异性成像中取得了很好的应用效果（Wang et al.,2020）.

图9 通过使用背景噪声信号改进南加州地壳剪切波速度模型.从左到右分别是M16（只使用天然地震记录）、M21（结合天然地震记录和背景噪声互相关函数），以及两者之间的差（Diff）（修改自 Wang et al.,2019 中的图9）Fig.9 Comparison of shear wave velocity models by incorporating ambient noise cross correlation functions.From left to right are model M16 (only constrained by earthquake records),M21(further constrained by ambient noise cross correlation functions)and differences between M16 and M21 (modified from Figure 9 in Wang et al.,2019)

上述研究中的一个重要假设是使用背景噪声互相关所得到的记录能够很好地近似两个台站之间的经验格林函数.这一假设只有在以下条件满足的情况下才成立：（1）激发背景噪声的震源必须是均匀并且同势能分布.（2）所使用的噪声记录时间足够长，进而得到的互相关函数能够稳定地收敛到经验格林函数.然而这些假设在实际情况下是很难得到满足的.为了解决这一问题，Tromp 等（2010）提出了在伴随层析成像的框架下，通过对比观测和计算的集合互相关函数，来同时反演地球速度结构以及背景噪声的震源分布.这样能够避免使用均匀且同势能噪声来源的假设.这一方法目前被Ermert 等（2016）用来研究激发全球背景噪声的震源分布.

2.8 地壳速度反演

相对于传统的射线走时层析成像，全波形反演具有更高的分辨率以及能够更好地反演具体的速度扰动大小.这里列举一些重要的地壳成像方面的应用.（1）盆 地 结 构.Tape 等（2009）和Lee 等（2014）早期将伴随层析成像和全波形反演应用到研究南加州地壳当中.他们发现相对于以前的走时层析成像所得到的三维地壳模型，通过拟合三分量的地震图，新的模型中所成像的盆地结构更加清晰，尤其是盆地以及近地表的低速带相对于初始模型的扰动可以到达20%～30%.这么大的速度扰动是很难通过传统的射线走时层析成像得到的（图10）；（2）断层结构.南加州具有很多大型的断层，比如圣安德烈斯和Garlock 断层等.在Tape 和Lee 等的三维模型中，我们可以清晰地看到跨过断层带时所出现的大的速度变化，这与早期的到时层析成像的结果相吻合（Lin et al.,2010），但是他们所得到的速度变化更大.相对于Tape 等的模型，Lee 等（2014）结合了背景噪声所得到的格林函数和天然地震数据，所以相对来说具有更好的成像覆盖率.同时上述模型与二维主动源反演速度模型拟合较好（图10）.以上的研究都是使用各向同性的速度反演.最近，Wang 等（2020）在Tape模型的基础上加入了背景噪声数据，反演了各向同性及径向各向异性的南加州地壳结构；（3）下地壳结构.通过结合伴随层析成像以及背景噪声所得到的经验格林函数，Chen 等（2014）反演了西藏东南部的地壳结构.他们发现中下地壳（约20 km 厚）具有很低的剪切波速度（-6%～-12%），这可能来自于印度板块和欧亚板块的强烈碰撞所导致的下地壳部分熔融.另外，使用全波形反演所得到的三维地壳模型可以更好地研究地震震源，比如地震定位以及震源机制解（Lee et al.,2011;Liu et al.,2004;Wang and Zhan,2020）.

图10 对比不同南加州地壳模型在二维LARSE-I（a）和LARSE-II（b）剖面上的P 波速度结果.CVM-S4.26、CVM-S4 和CVM-H11.9 分别是三个南加州地震研究中心的标准公共地壳模型.Lutter 等（1999,2004）以及Fuis 等（2003）分别是两个二维主动源地壳速度反演结果（修改自 Lee et al.,2014 中的图8）Fig.10 Comparisons of P wave velocities in different Southern Californian velocity models for 2D LARSE-I and LARSE-II cross sections.Models CVM-S4.26,CVM-S4 and CVM-H11.9 are three 3D reference crustal velocity models from the Southern California Earthquake Center (SCEC).Lutter et al.(1999,2004) and Fuis et al.(2003) are 2D crustal velocity models constrained by active source experiments (modified from Figure 8 in Lee et al.,2014)

上述地壳反演的另一个重要要求是，需要设计一个模拟区域能够把尽可能多的震源和台站包括其中.对于地震活动不频繁、但有密集地震台阵的区域，常常无法只使用本地的震源来进行全波形反演.一种可以在地壳和上地幔反演中同时考虑远震震源的方法是在三维反演区域之外，使用一些快速算法和较简单的背景模型来计算远震波场（比如FK 或者DSM 方法）.然后将计算所得到的波场当作边界条件，而在所要研究的区域三维模型中仍然使用数值算法（比如有限差分或者谱元法）来计算三维波场以及目标函数梯度.这种混合数值模拟波场的计算方法使得我们可以在三维地壳波速反演以及反射面成像当中使用远震记录（Beller et al.,2017;Monteiller et al.,2015,2021;Pienkowska et al.,2021;Tong et al.,2014;Wang et al.,2016），从而大大地改进密集台阵下面的结构分辨率以及成像覆盖率.例如，Wang 等（2016）展示只要用极少数的几个远震，并使用一个近线性台阵上的P 波及其散射和折射波的记录，就可以反演出欧洲西南部比利牛斯山脉（Pyrenees）下方岩石圈的纵波和横波速度.相似的方法也被应用到反演区域在震源和台站中间的情况.从震源用快速算法和简单模型传到模拟区域的地震波经过三维数值模拟后，仍需要通过边界条件用快速算法和简单地球模型传回到接收台站.这类混合数值计算波场的方法被用于全波形反演中，也被称为box tomography，并被应用于上地幔成像当中（Clouzet et al.,2018;Masson and Romanowicz,2017）.

2.9 地幔速度反演

如上所述，伴随层析成像和全波形反演已经被广泛地应用于很多地幔尺度的研究当中.本文列举一些与板块俯冲带和地幔柱相关的典型研究成果.

（1）板块俯冲带成像.研究海洋板块俯冲的几何特征以及其波速结构对于我们了解板块构造、地球动力学以及矿物物理具有十分重要的意义（Stern,2002）.目前在大陆尺度上，全波形反演能够更好地刻画海洋岩石圈在俯冲过程当中所经历的复杂过程（Tao et al.,2018;Zhu et al.,2012）.相关的成像结果与传统射线走时层析成像结果相当接近.但同时通过提高模型分辨率以及更好地反演具体的波速大小，我们也发现了一些新的现象.比如Zhu 等（2012）在欧洲区域的研究中能够清晰地看到板片滑脱（slab detachment），这一现象与以前走时层析成像的结果相似（Wortel and Spakman,2000），但在新的模型中我们能够找到更多的区域具有此滑脱现象.另外，现实中海洋板块在俯冲过程当中会发生断裂，而且这一断裂面会沿着海沟方向不断扩大，进而导致滑脱不断增大.这些断裂会导致软流圈物质上涌，进而与地表的火山作用相关.最近，Zhu 等（2020c）在中南美洲的全波形反演结果中发现Rivera 和Cocos 板块在俯冲过程中也出现类似的滑脱现象.同时在板块断裂区间填充有大量的低速带物质，其中可能发生部分熔融现象.另一方面，通过反演方位各向异性，他们发现板块的快速后退以及板块之间的空隙会产生巨大的压强差从而导致局部的软流圈物质流动.这一流动会使得橄榄石等各向异性矿物的快速轴方向沿着流动/应变方向定向排布，从而产生可观测的地震波各向异性.这对于我们研究俯冲板块与周围的地幔物质相互对流作用具有十分重要的意义.另一方面，在西太平洋和东北亚，早期的射线走时层析成像已经发现太平洋板块在俯冲过程中与上地幔，特别是地幔转换带具有复杂的关系（Fukao and Obayashi,2013;Zhao et al.,1992,1994）.比如在东北亚，太平洋板块在地幔转换带中可以横向地延展上千千米.同时Tang 等（2014）的研究显示在地幔过渡带中，太平洋板块发生断裂，导致地幔物质上涌，进而与地表的长白山五大连池火山相连.Tao 等（2018）通过使用全波形反演更好地描述这一俯冲板块在上地幔，特别是地幔过渡带中的特征.如图11 所示，通过拟合地幔过渡带的SH 波形，相对于以前的成像结果，他们的俯冲板块分辨率更高，成像结果更加清晰.在全球尺度上，相关的伴随层析成像结果也更好地显示俯冲板块在下地幔的具体形态.而这些特征在传统的走时层析成像结果当中常常并不十分清楚，进而导致模型解释上的一系列争论.通过更加清晰地描绘俯冲板块在地幔过渡带和下地幔的几何学特征，能够更好地帮助我们了解地幔对流以及地球内部的动力学过程（Lei et al.,2020）.

图11 使用通过地幔转换带的SH 波改进俯冲板块伴随层析成像结果.（a）显示的是二维剖面和所使用的地震台站的位置.（b）对比初始模型（左）和迭代改进之后的模型（右）对SH 地幔转换波的影响.黑色和红色地震图分别是实际观测和模拟计算所得到的结果.（c）和（d）分别是初始模型和迭代改进之后的模型在（a）中所示的二维剖面上的剪切波速度扰动（修改自 Tao et al.,2018 中的图6）Fig.11 Imaging subducting slabs by using SH triplication waveforms.Panel (a) shows the locations of 2D vertical cross sections and stations.Panel (b) compares observed (black) and synthetic (red) SH waveforms from the starting (left) and updated models(right).Panels (c) and (d) show shear wave velocity perturbations for the starting and updated models on the 2D vertical cross section shown in Panel (a) (modified from Figure 6 Tao et al.,2018)

（2）地幔柱成像.相对于海洋板块的俯冲，地幔柱提供了地幔对流的上升区域.从地幔柱假说提出至今（Morgan,1971），地震学界对地幔柱的成像一直十分关注，而且不同的成像方法常常得到了不同结论（Montelli et al.,2004b,2006）.比如有的研究显示冰岛的地幔柱来自于核幔边界（Bijwaard and Spakman,1999），而有的研究显示这一区域的低速带主要局限于上地幔（Ritsema et al.,1999）.相对于高波速的海洋板块俯冲带来说，低波速的地幔柱对于成像精度和分辨率的要求更大.由于波前愈合以及成像分辨率的局限，传统的走时层析成像很难反演几百千米直径的上升低速地幔柱.相反，全波形反演考虑了上述的物理过程，从而能够更好地帮助我们反演尺度更小的低波速带.比如Richers 等（2013）发现冰岛地幔柱在上地幔中向东部偏移，可能与地幔对流方向有关 .在全球尺度上，新的成像结果发现了更多的地幔柱，它们连接核幔边界和地表的热点.French 和Romanowicz（2015）还发现这些地幔柱与软流圈上垂直于洋中脊的低速带相连.值得注意的是这一研究使用的是谱元方法来正演模拟地震波传播，但是使用的是近似的方法（nonlinear asymptotic coupling theory,NACT）来计算敏感核.相对于伴随层析成像和全波形反演，这样能够大大地减少计算量.如图12 所示，Lei 等（2020）在全球伴随层析成像模型当中发现南非地幔柱在上升到地幔转换带时，分成多个更小的地幔柱，然后上升连接到地表的火山和热点.这些结果大大地提高了对地幔柱的认识，从而为更好地了解地幔对流提供了巨大的帮助.

图12 对比三个全球尺度地幔剪切波成像模型在6 个二维剖面上的结果.从左到右分别是来自于GLAD-M25（Lei et al.,2020）、TX2015（Lu and Grand,2016）和SEMUCB-WM1（French and Romanowicz,2015）.（a-f）分别对应于以下的热点:（a）阿法尔州;（b）百慕大群岛和加那利群岛;（c）佛得角和哈加尔高原;（d）冰岛和艾费尔高原;（e）复活节岛 和加拉帕戈斯群岛；（f）马里昂县和凯尔盖朗群岛（修改自 Lei et al.,2020 中的图15）Fig.12 Comparisons of shear wave velocity perturbations in three global-scale mantle tomography models.From left to right are results from GLAD-M25 (Lei et al.,2020)，TX2015 (Lu and Grand,2016) and SEMUCB-WM1 (French and Romanowicz,2015).Panels (a) to (f) show results for the following hotspots: (a) Afar;(b) Bermuda and Canary;(c) Cape Verde and Hoggar;(d) Iceland and Eifel;(e) Easter and Galapagos and (f) Marion and Kergulen (modified from Figure 15 in Lei et al.,2020)

2.10 油气储层反演

虽然我们对全波形反演和伴随方法早在1980 年代就已经在理论和实践上进行了探讨（Tarantola,1984,1988），但是直到2007 年左右，因为其在二维盲测实验当中所取得的成功（Brenders and Pratt,2007），才又一次得到人们的关注.这来源于过去几十年中计算能力和算法精度的提高.目前全波形反演已经被广泛地应用于能源勘探项目当中，为复杂条件下寻找油气资源提供了巨大的帮助（Operto et al.,2015;Operto and Miniuss,2018）.比如，图13 对比在北海Vahall 油田，通过拟合海底地震台站观测波形所得到的速度模型与传统的反射到时层析成像的结果，其中全波形反演所得到的结果能够更好地寻找油气储层.相对于经典的全波形反演，目前在勘探领域有以下的一些发展.

图13 对比全波形反演在北海Valhall 油田勘探当中的应用.左图和右图分别是使用反射波到时层析成像和全波形反演所得到的结果.（a-c）以及（h-j）分别是在175 m、500 m 和1 000 m 所得到的P 波速度结果.（d-g）以及（k-n）分别是相应虚线位置的二维剖面结果（修改自 Operto et al.,2015 中的图4 和图11）Fig.13 Comparison of velocity models for the ocean bottom node data from the Valhall oilfield.The left panel and right panel are models from reflection travel time tomography and full waveform inversion,respectively.(a-c) and (h-j) in each panel are P wave velocities at 175 m,500 m and 1000 m.(d-g) and (k-n) demonstrate results in corresponding vertical profiles as shown by dashed lines (modified from Figures 4 and 11 in Operto et al.,2015)

（1）时间与频率域的选择.相对于早期的时间域正演和反演，使用频率域、Laplace 域以及Fourier-Laplace 域的反演具有以下几个优点（卞爱飞等,2010;Pratt and Worthington,1990;Pratt,1999;Pratt and Shipp,1999;Shin and Cha,2008,2009）：第一，如果能够在三维空间下快速而准确地求解Helmholtz 方程，那么在频率域以及Laplace 域的计算量将与震源数目无关.这大大地减少了全波形反演的计算量，对大型三维勘探项目尤其重要（Engquist and Ying,2011;Operto et al.,2007;Poulson et al.,2013）.第二，在频率域以及Laplace 域中，可以更好地使用多尺度反演策略.相关的研究表明，只需要使用几个特定的频率或者频率组合就能够很好地覆盖所要研究对象的波数分布（Bunks et al.,1995）.第三，相对于时间域，在频率域下可以更直接地描述地震波的衰减，进而更好地反演品质因子Q的结构（Kamei and Pratt,2013;Prieux et al.,2013）.

（2）反射波全波形反演.早期的二维全波形理论实验表明，如果只使用透射波，比如井间层析成像，则可以很好地约束低波数信息，即光滑的速度结构.相反，如果只使用反射波（地表接收信号），则主要只能约束高波数的速度分布（Gauthier et al.,1986;Mora,1987）.如何使用反射波来约束低波数的速度模型一直以来是一个重要的研究问题.通过构建反射波的敏感核，可以使用反射波的到时以及波形来达到约束低波数的速度模型的效果.相对于直达P 波的敏感核，反射波的敏感核被形象地称为Rabbit ear（Ma and Hale,2013;姚刚和吴迪,2017;Zhou et al.,2015）.它与传统的偏移脉冲响应（migration impulse response）不同，后者主要描述高波数的速度结构（反射率）（图14）.

图14 在二维均匀速度模型下对比直达波和反射波的梯度.（a）和（b）分别是直达波和反射波的梯度，其中反射界面如图（b）中的紫色直线所示.（c）直达波+反射波的梯度.（d）改进的联合全波形梯度结果（修改自 Zhou et al.,2015 中的图3）Fig.14 Gradients for direct and reflected waves in a 2D homogeneous velocity model.Panels (a) and (b) are gradients for direct and reflected P waves.The reflector is shown as the purple line in Panel (b).Panels (c) and (d) are gradients for direct+reflected waves,and joint full waveform inversion gradient (modified from Figure 3 in Zhou et al.,2015)

（3）高频成像与反演.在传统勘探开发当中，我们使用层析成像来反演光滑的低波数速度结构.然后使用偏移方法，比如逆时偏移成像（McMechan,1983），来找到相应的高波数反射界面以及油气储层位置.随着全波形反演的发展，它被逐渐地推广到高频段.目前的反演频率可以达到30～50 Hz.如果仍然能够使这一非线性反演问题得到收敛，而且避免前述的周期跳跃问题，在未来有希望使用全波形反演来替代以前的层析成像+偏移的勘探开发步骤.

（4）与全局反演的结合.如前所述，全波形反演本质上是一种求解局部非线性最优化的问题.它只能找到目标函数的局部最优解，而无法保证得到全局最优，同时也很难提供相关的不确定性信息.目前，有相关的研究结合数值正演以及全局反演方法，比如MCMC 或者Hamiltonian MCMC 方法，来进行速度模型的采样（Gebraad et al.,2020;Stuart et al.,2019;Zhao and Sen,2021）.后者相对于前者，由于使用了目标函数梯度的信息，能够使得反演过程更快地收敛（图15）.这一点对该类方法十分重要.由于这些取样涉及到上千甚至上万的正演模拟，在目前的计算能力下做三维规模化推广还是十分具有挑战性的.所以目前大部分类似的研究还是在二维模拟反演当中.如果能够使用基函数来更好地参数化模型，则有希望大幅地减少模型空间的大小，进而加快采样的收敛速率.

图15 使用基于目标函数梯度的马尔科夫链蒙特卡罗方法反演二维Marmousi 模型.（a-c）分别表示平均模型、标准差模型以及平均与真实模型之间的差.（d）显示6 个采样点的一维边界后验概率分布情况（修改自 Zhao and Sen,2021 中的图9 和10）Fig.15 Results for the 2D Marmousi model from a gradient based Markov chain Monte Carlo sampling.Panels (a) to (c) are results for the mean velocity model,standard deviation and differences between the mean and true velocity models.Panel (d) shows the 1D marginal posteriori probability density functions at six different locations (shown in Panel a) (modified from Figures 9 and 10 in Zhao and Sen,2021)

其他的勘探领域的全波形反演问题涉及到前面所述的设计优化的目标函数，多参数联合反演，与机器学习的结合以及在实际数据中提取低频信号等，这里不再做赘述.

3 结论与展望

全波形反演目前已经被广泛地应用于多尺度的地球内部成像研究当中，它可以更好地寻找油气资源，描述岩石圈的波速结构进而更好地研究地震震源分布和机制，同时还可以用于研究板块俯冲过程以及地幔柱的结构.这些成果已经帮助我们更好地研究地球内部的物质组成以及动力学过程.相关的研究成果对岩石矿物物理和地球动力学具有十分重要的科学意义.目前主要的研究问题还包括选取优化的目标函数，如何评价反演结果的分辨率以及不确定性，多参数联合反演以及反演结果的解释等.相关方面的突破在未来可以更好地了解地球以及其他行星的内部物理及动力学过程.

致谢

衷心感谢中国科学技术大学姚华建教授的邀请，以及三位评审专家提出的宝贵修改意见.同时感谢编辑部老师在稿件修改和发表过程中的大力支持.