杨顶辉等：高分辨率全波形地震成像研究——进展、挑战与展望

高分辨率全波形地震成像研究——进展、挑战与展望

杨顶辉¹，董兴朋²，黄建东¹，方志龙³，黄雪源⁴，刘少林⁵，刘梦雪¹，蒙伟娟¹

1 清华大学数学科学系

2 中国地震局地质研究所

微信号：ziranziyuan1
更多资料，添加微信
复制微信号

3 电子科技大学资源与环境学院

4 北京工商大学数学与统计学院

5 应急管理部国家自然灾害防治研究院

第一作者/通讯作者：杨顶辉，清华大学数学科学系，教授。1982.9－1986.7云南大学数学系计算数学专业，获学士学位1990.9－1993.7中国科学院计算中心计算数学专业，获理学硕士学位1993.9－1996.6中国科学院地球物理所固体地球物理专业，获理学博士学位。（说明：作者介绍来源于百度百科。本文发表于《中国科学 地球科学》，作者都是著名专家学者，论文中没有作者介绍）

导读（摘要）：

本文系统性综述了高分辨率非线性FWI成像方法的基本理论、研究进展、多学科应用、存在问题及解决思路，展望了未来的发展方向和研究重点。

基金项目：国家自然科学基金重点项目(42330801)资助.

说明：参考文献以原文为准，本推文未作详细标注。

------内容提纲------

1 引言

2 FWI理论和方法研究进展

2.1 FWI基本理论

2.2 地震波场计算

2.3 初始模型构建

2.4 目标函数

2.5 优化算法

2.6 正则化与凸约束优化

2.7 多尺度和多数据反演

2.8 不确定度分析

3 多学科FWI成像实际应用

3.1 油气勘探

3.2 深部结构成像

3.3 地震孕育和发生

3.4 工程地球物理勘探

3.5 医学成像

4 面临的挑战问题与重点研究内容

4.1 解的非唯一性问题

4.2 计算和存储需求的问题

4.3 多参数反演耦合问题和联合成像方法

4.4 震相提取和匹配困难

4.5 人工智能技术的广泛应用

5 FWI成像方法未来应用领域

5.1 4D地震成像

5.2 地球动力学模拟

5.3 行星成像

5.4 医学和工程领域成像

6 总结与展望

6.1 总结

6.2 展望

-----------

1  引言

地球科学既是人类探索地球之谜的古老学科，同时也是问题驱动的前沿科学。由于地球内部的不可入性，使得深入认识地球内部结构、物质组成和演化机制面临巨大挑战。澳大利亚于1999年率先提出“玻璃地球”概念，目标是使澳大利亚大陆地表以下1000m深度以内地质情况变得透明(刘树臣,2003)。2001年，美国启动了“地球透镜”计划，旨在探索北美大陆的构造与演化，并揭示地震和火山喷发的形成机制（https：//www。earthscope-program-2003-2018.org）。2014年，在中国召开了以“中国玻璃地球建设核心技术及发展战略”为主题的香山科学会议，与会者一致认为，“玻璃地球”建设可为探索新一代巨型矿床、揭示地质灾害形成机制提供科学支持，提升国家在资源探测和灾害防御方面的能力。2017年，中国地震局提出了“透明地壳”计划（http：//neobji.ac.cn/article/1256），计划在10年内完成重点地区地下结构、构造和地球物理场变化的观测；随后又提出“解剖地震”计划，以探索地震孕育机理为目标，逐步深化对地震孕育发生规律的认识。类似于美国南加州地震中心(SCEC)，当前正在建设中的中国地震科学实验场，其科学目标之一是揭示大陆型强震孕震体介质结构和空间分布特征，探索孕震过程与应力应变的关系以及大陆型强震孕育演化过程。这一系列计划的总体目标都是为了深入了解和认识地球，其中地震反演成像理论和技术占有重要地位，它是实现上述各种地球探测计划和发展战略目标的关键技术，是认识地球内部结构和物质组成、揭示地震孕育发生和地球动力学演化过程的重要手段。

自地震层析成像方法被用于揭示地球内部的非均匀结构以来，射线类成像方法一直是地震反演成像的主要工具（Aki和Lee，1976；Aki等，1977）。过去四十多年中，许多求解地震波传播路径和地震波走时的射线类方法相继提出，并在实践中占据了主导地位，其中打靶法和弯曲法是两种常见的射线追踪方法（Julian和Gubbins，1977；Langan等，1985）。Zhao等（1992）发展了三维动态射线追踪方法，能够快速计算存在速度不连续面的三维非均匀介质中地震波的传播路径和时间，并被广泛应用（Huang和Zhao，2006；Tian等，2009）。然而，传统射线理论假设地震波为无限高频，这与地震波是有限频率的实际情况不符，导致其理论上难以分辨尺度小于地震波主波长的异常体（Rawlinson等，2010）。针对这一问题，考虑了地震波有限频率特征的有限频率层析成像方法得到了快速发展（Dahlen等，2000；Hung等，2000，2004；Zhao等，2000；任和Shen，2008；Liang等，2011；徐小兵等，2018）。尽管如此，Dahlen等(2000)发展的有限频成像方法是基于一维层状地球模型计算射线路径以及二维或三维敏感核，当模型存在很强的非均匀性时(例如地壳浅部)，该方法的敏感核精度可能不足，影响了其对复杂地下结构的解析能力。此外，该方法在处理转换震相方面也存在困难，因此相比于传统的射线成像方法，其成像分辨率的提升仍然有限。

为了进一步深化对地球内部结构和动力学机制的认识，推进矿产、石油、天然气、地热等资源的勘探开发，以及揭示地震、火山等自然灾害的发生机制，科学有效地进行防灾减灾，亟需发展具有更高分辨率的地震成像理论和方法。在这些实际问题的驱动下，Tarantola于20世纪80年代率先提出了基于广义最小二乘的时间域全波形反演（FWI）成像方法及其基础理论（Tarantola，1984）。相比传统射线追踪方法，FWI能充分利用地震波形记录中所包含的振幅、相位等信息，具备高分辨地震成像的潜力。随着计算机技术的发展，FWI技术研究和应用进入了快速发展期。在理论研究方面，Tromp等(2005)将波形层析成像和伴随方法置于统一框架之下，即在每一次迭代中，只需要进行基于已知信息进行的正演计算和利用波形残差、走时残差等与目标函数相关的伴随震源逆时反传所需要的伴随波场计算。在这一框架下，Tromp等（2005）给出了使用不同信息如波形、走时、振幅等作为目标函数的梯度计算方法，以及反演地震波各向异性、衰减、震源机制和震源参数以及地下不连续面的方法，极大地推动了基于三维时间域弹性和黏弹性波动方程的FWI的发展。研究人员逐步将基于时间域二维声波方程的FWI技术，推广并发展至基于三维时间域和频率域弹性波方程的FWI技术。通过与应用数学、计算机科学等学科的交叉融合，持续优化和改进FWI技术，逐步解决全波形反演所面临的计算量和储存量大、解不唯一、易陷入局部最优、多参数耦合等关键技术问题。在实际应用方面，FWI技术已被广泛应用于深部地球结构成像、油气资源勘探、二氧化碳地质封存等多个领域，并取得了显著进展，显示出相较于传统射线追踪方法具有更高的分辨率和广泛的应用前景。

本文系统性综述了高分辨率非线性FWI成像方法的基本理论、研究进展、多学科应用、存在问题及解决思路，展望了未来的发展方向和研究重点，以期进一步推动高分辨率FWI成像技术在理论方法研究和实际应用方面的新发展。

2  FWI理论和方法研究进展

FWI成像方法是当前研究地球内部结构和岩石动力学的重要手段之一，在数学和物理两方面都具有显著优势。从数学模型角度看，FWI依赖于波动方程(如声波方程、弹性波方程和黏弹性波动方程等)，与传统的射线类成像方法相比，波动方程理论能更准确地刻画地震波在地下介质中的传播，其模拟的地震波传播过程也更贴近实际情况。从物理模型角度看，FWI成像不仅充分利用了地震波的走时信息(运动学信息)，还包含了各种震相(如直达波、反射波、折射波、散射波和面波等)的丰富动力学信息(如振幅和相位)以及震源的频谱信息等。这些地震波形携带了地球内部介质的各类物理参数信息，从而在理论上为优化反演提供了更多的数据和信息约束。

尽管20世纪Tarantola（1984，1986）就已经提出了利用地震波形数据的FWI成像基本理论和方法，但由于当时计算条件和台站数据的限制，难以开展大规模的全波形反演工作，相关研究未能取得实质性进展。进入21世纪，随着大规模高性能并行计算设备的出现，以及一系列非均匀地球介质中地震波传播模拟的高精度快速正演计算方法的不断提出，加之随着地震仪精度和稳定性的提高，高质量宽频带地震数据日益丰富，这些进步使得FWI成像方法在地球科学中的应用成为可能。另一方面，传统的地震成像技术难以满足高分辨率成像的地球物理勘探需求，尤其是难以满足复杂地质条件下油气勘探的实际需要。进入21世纪以来，具有高分辨率成像优势的FWI成像方法得到了快速发展，成为了学术和工业界广泛关注的焦点。图1展示了射线层析成像与全波形反演成像结果的对比（Pratt等，2002）。结果显示，FWI成像结果的分辨率显著高于射线层析成像的结果，能够更好地重建地下模型的细节，而射线层析成像的结果则相对模糊。

图1 射线层析成像与全波形反演结果对比

从上至下，每一行依次代表真实模型(上)、射线层析成像结果(中)和FWI成像结果(下)。其中第一行左、中、右三个真实检测板模型分别代表具有不同非均匀块体大小的三种情况。模型顶部和底部分别均匀布设了101个震源和101个接收器，距离边界均为150m，因此每个模型共计生成101×101=10201条地震记录。三种检测板棋盘格模型的单元尺寸分别为10。5、21和30m。地震波的传播模拟基于有限差分声波方程。在射线层析成像中，使用的地震波主频为300Hz；在全波形反演中，地震波的频带范围为50~500Hz。修改自Pratt等（2002）中的图1和2。

2.1 FWI基本理论

FWI是一个在偏微分波动方程约束条件下的非线性优化问题，其核心思想是通过逐步迭代模型参数m，最小化观测数据d与正演合成数据F(m)之间的残差。在给定观测数据dd的情况下，FWI通过求解以下最优化 

其中，n_src来源代表震源数量；ϕ(⋅,⋅)φ(⋅,⋅)代表用以度量观测数据和合成数据差异的误差函数，是影响FWI成像效果的重要因素。Tarantola（1984）基于广义最小二乘理论，采用L₂范数度量误差的目标函数ϕ(⋅,⋅)φ(⋅,⋅)，最先提出了如下形式的FWI最优化问题：

其中，正演算子F(m)需要在时间域或者频率域求解对应的波动方程。采用矩阵符号的表达方式（Marfurt，1984；Carcione等，2002），在时间域和频率域中，FWI分别需要求解如下方程：

时间域：

其中，s(x，t)、s(x，ω)、u(x，t)和u(x，ω)u分别代表了时间域和频率域的震源项和地震波的波场，t、x和ω分别代表了时间、空间坐标和角频率。M(x)、A(x)和B(x)分别代表了质量矩阵、刚度矩阵和阻抗矩阵，它们是模型参数m的函数，以声波方程为例，模型参数mm为慢度的平方，相应的时间域和频率域波动方程为

时间域：    

针对式（2）中的目标函数ϕ₂(m)，Tarantola（1984）通过将检波器点的波场残差作为伴随震源并计算其对应的伴随波场，提出了用于计算梯度向量g和海森矩阵∇²ϕ₂(m)的伴随方法，从而奠定了求解FWI优化问题的理论基础。具体来说，以频率域FWI为例，考虑单个震源和单个频率，其梯度向量g有如下的数学表达式：

其中，v为以波形数据残差为震源的伴随波场，P为从全波场中提取检波器位置信息的投影矩阵，R(⋅)为对复数取实部的算子。从公式（7）和（8）中可以发现，针对一个震源产生的数据，计算其相应的目标函数和对应的梯度向量，需要求解一个正演波场u和一个伴随波场v，这构成了FWI最主要的计算量。时间域FWI的梯度计算公式也可由相应的矩阵符号表示推导获得（Tromp等，2005），其具体表达形式如下：

其中，u(x，t)和v(x，t)分别代表时间域的正演波场和伴随波场，L(x)代表时间域波动方程的微分算子，L(x)T代表其伴随算子。

求解最优化问题(2)，可以选择全局搜索算法，如遗传算法、模拟退火算法、马尔科夫-蒙特卡罗算法等(魏超等，2008；Martin等，2012），也可以采用局部搜索算法，如最速下降法、共轭梯度法（Akcelik，2002）、牛顿法、拟牛顿法（Nocedal和Wright，2006）和高斯-牛顿法（Pratt等，1998；Epanomeritakis等，2008）等。虽然全局搜索算法具备获得全局最优解的能力，但是其往往需要局部搜索算法几十倍甚至上百倍的迭代次数（Bui-Thanh等，2013）。鉴于FWI每步反演迭代都需要花费很大的计算量来数值求解波动方程，因此全局搜素算法难以在大规模实际应用中推广。相比之下，局部搜索算法因其在计算效率上的优势，目前被广泛应用于各类实际问题的处理。在每一次迭代计算中，局部搜索算法利用目标函数的一阶梯度信息和二阶海森矩阵信息对模型参数的迭代更新量进行计算，具体可用如下通用格式表示：

度向量，矩阵C代表通用的更新量预处理矩阵。特别地，当C=αI、∇²ϕ₂(mk)或有限记忆Broyden-Fletcher-Goldfarb-Shanno（L-BFGS）Hessian矩阵时，公式（12）分别对应于使用最速下降法、牛顿法和L-BFGS方法来迭代更新模型参数。

综上所述，FWI的关键步骤包括正演波场和伴随波场的计算、优化算法的选择、目标函数的选取以及初始模型的构建。图2展示了FWI成像的主要处理流程。在下面小节中，我们将就这几方面的研究进展进行详细介绍。

图2 FWI流程图

2.2 地震波场计算

正演波场和伴随波场的计算是FWI成像过程中至关重要的两步。伴随波场的计算可以采用与正演波场模拟相同的数值计算方法。我们通常使用数值方法来计算这些波场，这构成了FWI最主要的计算量来源。目前，求解波动方程的正演方法主要包括：有限差分方法（如：Virieux，1986；Graves，1996年）、有限元方法（如：Marfurt，1984）、谱元法（如：Komatitsch和Tromp，1999，2002）、边界元法（如：Chen等，2012；Pourahmadian和Guzina，2015）、有限体积法（如：Morton，1998；Schneider等，1999），以及其他类型的数值方法，包括频率-空间域的混合方法（如：Shin和Sohn，1998；Hustedt等，2004；Fan等，2018）、近似解析离散化及其优化方法（如：Yang等，2003，2006；Tong等，2013）、保辛几何结构算法（如：马等，2011，2018），近年来迅速发展的间断有限元（DG）方法（如：Käser和Dumbser，2006；Käser等，2007；He等，2015，2019，2021；Huang等，2023a，2023b）和基于深度学习的人工智能方法（Song等，2023；Zheng和Wang，2023；Chang等，2024）等。这些方法各有优缺点。例如，频率域差分方法能够有效处理多震源问题和刻画介质的衰减性，但其并行性差，难以实现大规模三维地震成像（Song和Williamson，1995）。谱元法作为一种高阶的有限元方法，易于处理起伏地表和间断面，但其在处理三维问题时，相比有限差分方法，同样面临着巨大的计算量和存储量挑战，这在一定程度上限制了其应用(例如，2010年Fichtner等对澳大利亚大陆上地幔的研究，仅使用了57个地震的2137条波形数据)。在相同条件下，传统的有限差分方法具有存储量小，计算速度快，易于并行等多种优点而广受欢迎。然而，传统的有限差分方法存在严重的数值频散（Dablain，1986；Yang等，2002），且网格剖分不够灵活，这影响了正演计算的准确性并可能导致成像结果失真。为了压制数值频散，有限差分方法需要减小空间和时间步长，但这将大幅增加存储需求并降低计算效率，使得大规模三维成像变得更加困难。

实现大尺度三维成像的一个关键挑战是如何在粗网格情况下有效抑制数值频散并提高计算效率。近二十余年来，一系列能有效压制数值频散，且具备存储量小、计算效率高、易并行化计算（因采用局部离散算子）以及能同时计算位移和（或）粒子速度及其梯度场信息等优点的数值计算正演方法被提出。此外，通过优化差分系数来压制数值频散的方法也在不断发展。针对起伏地表问题，Zhang等（2012）发展了贴体网格有限差分方法，使得波场计算更为准确可靠。这些正演波场模拟方法有望在FWI地震成像及其应用中提高计算效率和成像精度。值得注意的是，间断有限元方法作为一种集有限元和有限体积法优势于一身的算法，近年来显现出巨大的潜力。间断有限元方法不仅具有较高的数值精度，还能在时间域对大规模线性矩阵进行有效解耦，具有高并行性，特别适合在GPU环境下进行高效并行计算。近年来，该方法在地球物理领域取得了显著进展，有望在未来的FWI成像中得到应用。然而，间断有限元方法面临的主要问题是其严格的稳定性条件，未来需要进一步研究以提高正演波场模拟的效率。

2.3 初始模型构建

FWI通过迭代优化的方法寻找并构建与实际观测数据最匹配的地下介质模型，其反演成像的关键之一是选取与真实地下情况接近的初始模型。由于波形数据与模型参数之间存在着极强的非线性关系，只有当初始模型足够接近真实模型时，FWI才可能正确地收敛到全局最优解。理论上，当初始模型生成的波形数据与观测波形数据之间存在超过半个波长的相位差时，基于L₂范数的FWI会遇到所谓的周期跳跃（Cycle-skipping）问题，导致反演收敛到局部最优解。然而，由于地下介质的复杂性，如介质的非均匀性、各向异性和衰减等因素的影响，在实际应用中，获取一个相对准确的初始模型往往比较困难。以中国华北克拉通东缘为例，超过10km厚的沉积层S波速度和鄂尔多斯块体下地壳S波速度，相比于PREM参考地球模型，速度异常远超过10%（Zheng等，2008）。在这些构造复杂且具有强烈非均匀区域进行FWI成像时，将面临反演结果较难收敛到真实模型的问题。

在实际应用中，通常利用走时层析成像（Hauser等，2008；Liu和Gu，2012）方法为FWI构建一个相对准确的初始速度模型。此外，在反演初期利用观测数据中的低频信息进行FWI，为后续高频数据的FWI提供初始模型也是一种广泛使用的技巧。相比于高频数据，低频数据的波长更长，因此更不容易发生周期跳跃。然而，在很多实际场景中，足够低频的数据往往比较难以获得。特别地，在地震勘探中3Hz以下的低频数据的信噪比往往非常低，几乎无法提取可用于FWI成像的有效信号。为了解决低频信号缺失的问题，Li和Demanet（2016）提出采用相位追踪和事件分离结合的方法，对低频数据进行恢复。随后，Sun和Demanet（2020）进一步采用深度学习的方法实现了低频信号恢复，缓解了周期跳跃的问题。

2.4 目标函数

FWI作为一个非线性优化问题，其目标函数的选取决定了反演结果的收敛性和可靠性。一方面，不同物理量（如到时、振幅、相位等）对同一个模型参数（如P波速度、S波速度、密度、弹性模量等）的敏感度不同，而且同一个物理量针对不同模型参数的敏感度也存在差异（Tong等，2014b）；另一方面，不同目标函数与模型参数之间的非线性关系，以及它们对反演结果非唯一性的约束能力也各不相同。因此，针对特定的模型参数，选取合适的物理量构建恰当的目标函数进行反演，对于模型参数的准确重构至关重要。

传统的FWI方法采用基于全波形数据残差的L₂范数作为目标函数，容易陷入局部最优解，进而导致反演结果的多解性。针对这一问题，研究人员从不同角度对目标函数进行了改进：

(1)通过选取不同于波形的其他物理量，包括振幅、相位、到时等来构建目标函数。Tromp等（2005）分别针对波形、振幅、衰减和到时信息提出了相应的目标函数，并采用伴随方法给出了相应模型参数梯度向量的计算方法。Wu等（2014）提出采用波形包络数据的反演方法，可以为基于波形数据的反演提供更为准确的初始模型。相比于原始波形数据，包络数据包含更低频的信息，可有效缓解周期跳跃问题。该目标函数被广泛应用于FWI成像研究中（敖瑞德等，2015；Xiong等，2023）。此外，利用波形数据提取走时信息进行反演也被广泛应用（Tong等，2014a）。Luo和Schuster（1991）首先提出基于波动方程的互相关走时层析成像反演方法。由于该方法只利用了波形中的走时信息，可认为其为一种拟波形成像(黄雪源等，2021)。随后Tape等（2009，2010）和vanLeeuwen等（2014）进一步发展了加权互相关走时层析成像方法，提高了反演精度，加强了对深部小尺度构造的刻画能力。为了减小震源参数不确定性引起的系统误差，Yuan等（2016）提出了双差波形误差函数，从而提高了台阵区域的成像分辨率。此外，在使用噪声面波数据时，该方法还能减少由于噪声源分布不均匀导致的对面波提取和成像质量的影响（Liu等，2023；Yan等，2024）。

（2）通过采用不同于L₂范数的其他度量方式来构建目标函数。Crase等（1990）针对非高斯噪声，提出了基于L₁范数的FWI目标函数。Aravkin等（2011）针对服从t分布的数据噪声问题，提出了基于t分布的目标函数，可有效处理观测数据中的异常噪声。Warner和Guasch（2016）提出了自适应波形反演方法。该方法构建了基于维纳滤波（Wienerfilter）测量两个波形相似度的目标函数。相关数值模拟实验显示该方法能够为基于波形差的FWI提供较为准确的初始模型（Zhu和Fomel，2016）。Engquist和Froese（2014）首次将最优输运理论中的Wasserstein度量（Villani，2003）应用于地震信号处理和反演。之后，Engquist等（2016）将二次Wasserstein（W₂）度量与FWI方法结合，获得了凸性更好的目标函数（图3），有效缓解了FWI方法存在的周期跳跃现象，避免了优化过程陷入局部极小值的问题，并通过数值实验验证了该方法的有效性（Yang和Engquist，2018；Yang等，2018）。鉴于最优输运理论在缓解FWI局部最优问题中的潜力，近年来开展了一系列相关的理论和实践研究。其中，Métivier等（2015，2016）将Kantorovic-Rubinstein（KR）度量（W₁度量的推广形式）与FWI相结合，获得了优于传统L₂度量的反演结果(图4)。Dong等(2019)提出了基于W₂度量的被动源FWI成像方法，并应用于川滇地区的地震成像。虽然基于最优输运理论的FWI相较于传统基于L₂范数的FWI成像具有优势，但在实际应用中仍面临诸多挑战，包括：(ⅰ)观测与合成信号的非负性要求；(ⅱ)观测波形与合成波形需保持质量守恒；(ⅲ)高维最优输运映射没有解析解，需依赖于数值求解复杂的Monge-Ampère方程，导致计算量大增。

图3 一维地震信号f和g在不同时移情况下的L₂和W₂距离

（a）观测信号f和合成信号g；（b）f和g信号在L度量和W度量下的距离，基于L度量(蓝色实线)的目标函数存在多个局部极小值点，而基于W度量(红色实线)的目标函数具有很好的凸性。修改自Dong等（2019）中的网络版附图S1和S3。

图4 BP模型

从左向右从上到下依次为：真实模型、初始模型、L₂度量反演结果和KR度量反演结果。修改自Métivier等（2015）中图13和16。

(3)通过在正问题方程中增加参数的额外自由度来扩展FWI的搜索空间，使合成数据匹配所需的精度，构建凸性更好的目标函数。基于这一思想，研究人员通过在模型域扩展参数，提出了扩展波恩建模法(Extended-Bornmodeling)(Symes，2008；Biondi和Almomin，2014；Barnier等，2023a，2023b)。同时，通过在震源域扩展参数，发展了包括对比源法(Contrast-sourcemethod)(vandenBerg和Kleinman，1997；Abubakar等，2011)、波形重构反演法(Wavefieldreconstructioninversion)(vanLeeuwen和Herrmann，2013b，2016；Fang等，2018a，2018b)，以及扩展源FWI(Extended-source FWI)(Wang等，2016；Huang等，2018)在内的多种算法。Fang和Demanet(2020)同时在模型域和震源域对参数进行扩展，提出了抬升松弛波形反演法(Lift and relaxwave for minversion)，进一步缓减了传统FWI易于陷入局部最优解的问题，减轻了其对于高质量低频数据和精确初始模型的依赖。

2.5 优化算法

FWI作为一种迭代优化的反演成像方法，需要选取适当的优化算法，以提升整体优化进程的收敛速度。常用的优化方法包括最速下降法、共轭梯度法、牛顿法、高斯牛顿法和L-BFGS拟牛顿法等。其中，最速下降法(Tong，2021)通过一阶泰勒展开式对函数近似，然后将变量沿着负梯度方向(下降最快的方向)移动一定的步长(通过直接求解获得或一维搜索算法获得)，使目标函数值不断下降直至接近极小值，但该方法存在收敛速度慢，对深部模型参数更新慢等问题。共轭梯度法(Hestenes和Stiefel，1952)是求解大型线性/非线性方程组最有效的算法之一，仅利用一阶导数信息，既克服了最速下降法收敛慢的缺点，又避免了牛顿法需要存储和计算Hessian矩阵及其求逆的缺点，但对于非二次问题，可能存在无法收敛的情况。L-BFGS方法(Liu和Nocedal，1989)是最有效的拟牛顿法之一，使用不包含二阶导数的矩阵近似牛顿法中Hessian矩阵的逆矩阵，在获得接近二阶收敛的同时，只需保存最近的k次迭代信息，从而提升了优化求解的收敛速度并减小了存储量。牛顿和高斯-牛顿算法作为二阶算法，能够获得更高的收敛率，但是其每一步迭代的计算量却是一阶算法和L-BFGS算法的十几到几十倍(Virieux和Operto，2009)。考虑到存储Hessian矩阵和敏感核矩阵的巨大存储量，实际计算中一般不显式地构建这两个矩阵，而是在共轭梯度算法的迭代中执行Hessian矩阵或者高斯-牛顿Hessian矩阵作用于一个向量的乘法运算(Virieux和Operto，2009)，该运算需要对每个炮点进行两次正演问题的计算，分别对应入射波场和伴随波场(Akcelik，2002)。因此，二阶优化算法的计算量与每个模型迭代中求解线性子问题所使用的共轭梯度迭代次数呈正比。

近年来，随着随机优化算法的快速发展，人们逐步将其应用于FWI中，极大地降低了FWI的计算量。由于声波、弹性波方程对于震源具有线性可加性，Haber等（2012）结合样本平均估计方法（Sample average approximation）和同时激发震源方法，提出了求解FWI的快速随机优化算法，受到广泛关注和应用。同时激发震源法需要所有震源共享同样的接收器，因此不适用于海洋勘探的拖缆采集。Van Leeuwen和Herrmann（2013a）提出了更具有普适性的随机L-BFGS算法，并应用于FWI。Fang等（2014）进一步将该方法应用于FWI的分辨率测试与参数不确定性的定量分析中。基于随机优化的算法能够将FWI提速至少一个量级（van Leeuwen等，2014），非常适合三维大规模问题的求解。

局部优化算法需要反演参数相对于目标函数的微分信息，即使采用了伴随方法，针对新参数和新目标函数的微分信息的建立仍然需要复杂的理论推导和繁重的计算。近年来，随着深度学习相关软件的快速发展，其内置的自动微分（auto-differentiation）工具（Paszke等，2017）为FWI的微分信息计算提供了新的实现手段。通过利用Tensorflow（Abadi等，2016）和Pytorch（Paszke等，2019）等深度学习软件实现FWI，可以利用其内置的自动微分工具快速获得任意目标函数对任意反演参数的微分计算，从而避免了复杂的目标函数微分信息理论数学推导和对应的编程实现，极大地方便了相关研究工作的开展。

2.6 正则化与凸约束优化

FWI面临多解性和局部最优的问题，模型参数的正则化是处理和解决这类问题的主要方法之一。此外，正则化方法还可以有效压制由数据中的噪声和数据采集不完备引入的假象，提高反演结果的信噪比。目前，被广泛使用的方法是Tikhonov正则化方法（Tikhonov和Arsenin，1977）。该方法通过在目标函数中引入诸如模型一阶梯度、二阶散度、全变差以及特定变换域的稀疏性等正则项，从而在拟合实际观测数据的同时，保证每一步模型更新都满足先验信息（Anagaw和Sacchi，2011；Asnaashari等，2013）。然而，Tikhonov正则化方法需要精细化地选择正则项的加权系数，这为其在实际场景中的应用带来了困难。

随着L₁谱投影梯度法(Spectral Projected Gradient for L₁ minimization，SPGL₁)(van Den Berg和Friedlander，2007)、投影拟牛顿算法(Projectedquasi-Newton，PQN)(Schmidt等，2009)、交替方向乘子法(Alternating Direction Method of Multipliers，ADMM)(Boyd等，2010)等凸约束优化算法的快速发展，人们开展了基于模型参数凸约束的FWI研究。相比于Tikhonov正则化方法，凸约束优化中的各模型约束条件更具物理解释意义，因此相关约束条件的参数也更容易选择。基于这种思想，Li等(2012)采用SPGL1优化方法提出了基于Curvelet域特征向量稀疏约束的FWI；Guitton等(2012)发展了基于Laplace变换域参数约束的FWI方法；Peters等(2019)采用PQN和ADMM方法提出了一整套适用于求解包含多项凸约束的FWI方法；Esser等(2016)提出了基于模型参数全变差约束的FWI方法，并首次对盐丘模型实现了高精度自动化重构，受到油气勘探界的广泛关注。

近年来，随着人工智能的迅速发展，尤其是生成式网络的快速发展，采用生成式网络总结地层结构先验信息并将其应用于FWI成像的方法逐渐兴起。Fang等(2020)和Mosser等(2020)分别提出了利用生成式对抗神经网络(GAN)(Goodfellow等，2014)对FWI进行正则化约束的方法。通过利用GAN的高分辨生成图像的先验信息，该类方法的反演结果往往拥有比传统理论值更高的图像分辨率。此外，该类方法还具有自动迭代更新先验知识库的优点，有助于对地层结构知识的长期积累、优化和应用(Marano等，2023)。

2.7 多尺度和多数据反演    

在FWI成像中，针对不同的研究对象选择合适的反演尺度是一个难题。为此，Fichtner等(2013)发展了多尺度网格的FWI方法，不同尺度的模型对应不同的网格空间。Dantas等(2019)提出了一种隐式的多尺度FWI方法，在时间域设计对不同波长敏感的目标函数，并使用相同的数值算法和数据集进行分步骤反演，获得了比传统FWI方法更好的成像结果。鉴于FWI方法可以使用全部震相进行反演，无需单独反演直达波、反射波或者面波等各种震相，因此相较于射线成像等传统成像方法，FWI方法可进一步提高成像结果的分辨率。在地震数据覆盖不足的区域，综合利用多种震相对提升成像质量尤为重要。Huang等(2016)基于FWI并利用多震相(P和S以及PmP和SmS)获得了1992年美国南加州Landers地震(M_w7.3)区的高分辨率地壳结构。这些研究表明，完善和发展多数据、多尺度的FWI成像方法对获得高质量成像结果非常重要。

2.8 不确定度分析

模型参数的不确定度分析对于解释成像结果的稳定性至关重要。广泛存在于数据、观测系统、震源和正演模拟中的不确定因素会通过反演传递到最终的反演结果中，因此需要对反演结果的可信度进行全面分析和评价，以便给出更为准确、科学的地质和物理解释。然而，对FWI的结果进行不确定分析仍然是一个极具挑战的问题。传统的棋盘测试计算成本太高，而“点扩散函数”测试方法只能提供特定感兴趣位置的不确定性分析（Fichtner和Trampert，2011；Zhu等，2016）。目前，FWI的不确定性分析技术主要基于贝叶斯推理方法构建模型参数的完整后验概率密度分布函数，通过研究该概率密度函数对各参数进行不确定性分析。理论上，通过马尔科夫链-蒙特卡罗算法（Martin等，2012；Fang等，2014）对后验概率密度函数进行精细采样，可以得到反演参数的完整概率分布。然而，该方法随着问题规模的增大会遇到“维度噩梦”的问题，计算量急剧增加，无法应用于三维大规模成像问题。为了应对这一挑战，研究人员尝试在全局最小值点附近进行采样，然而考虑到实际问题的复杂性，全局最小值点很难确定，因此通常在反演所得的近似解附近对后验模型分布进行采样，近似地对反演参数进行不确定性分析（Bui-Thanh等，2013）。目前，已有包括随机探测（Fichtner和vanLeeuwen，2015；Fang等，2018b），卡尔曼滤波（Thurin等，2019；Eikrem等，2019）和平方根变量度量方法（Liu等，2019；Liu和Peter，2019）等方法被相继提出，并在合成数据实验和实际数据应用中取得了成功（Hoffmann等，2024）。近期，结合生成式神经网络的不确定性分析方法为解决维度噩梦提供了新的思路。通过对现有地层模型进行训练学习，有望在高维模型空间中获得一个表征地层模型的低维特征流形，从而极大地降低所需采样的样本空间维度，减轻模型参数采样的计算压力（Fang和Demanet，2020；Mosser等，2020年）。尽管对地层模型进行有效训练学习需要额外的计算资源，但一旦成功训练出一个可靠的地层模型生成器，它可以显著简化后续反演成像过程中不确定性计算的复杂度。

3  多学科FWI成像实际应用

3.1 油气勘探

FWI成像理论的发展受益于油气勘探领域的迫切需求，并最先被应用于油气储层反演成像。由于浅层气云的覆盖(衰减屏蔽作用)，挪威北海油田的建模和成像一直是勘探领域的难点，甚至被称为勘探的盲区。Sirgue等（2010）在挪威北海油田利用OBC数据率先实现了三维FWI成像，相比反射波走时层析成像，FWI成像结果具有更高的分辨率。Operto等（2015）利用横向各向同性介质中的黏声波方程和频率域FWI方法反演获得了挪威北海Valhall油田的储层速度结构，并采用分频带策略进行反演，获得了分辨率更高的反演结果（图5），不仅精细刻画了气云形态，还准确展示了周围的断裂构造。周等（2015）将反射波FWI应用于挪威北海Valhall油田的实际数据也获得了高分辨率的成像结果（图6）。FWI反演成像在勘探领域的成功，为油气储层精确定位和解释提供了重要参考，极大地推动了FWI研究的热潮。

图5 北海Valhall油田基于不同频率的FWI成像对比

（a）~（d）分别是基于3。5、5、7和10Hz的反演结果。修改自Operto等（2015）中的图14.

图6 不同波场数据的梯度和反演结果对比

(a)~(d)依次是全波场梯度、反射波场梯度、反射波+反射率波场梯度，以及联合全波场梯度；（e）~（f）分别是反射波和联合全波场方程的阻抗模型；（g）~（i）依次是常规FWI、反射波FWI和联合全波场结果。修改自周等（2015）中的图3和9.

3.2 深部结构成像

FWI成像除了在油气勘探方面得到了广泛应用之外，在天然地震领域也备受关注，并分别在美洲、澳洲、欧洲、亚洲等区域取得成功的应用范例。Chen等(2007)基于散射积分的波形层析成像方法，获得了洛杉矶地区三维精细地壳结构模型。Tape等(2009，2010)利用谱元法高精度模拟地震波正传波场和伴随波场在复杂地质构造条件下的传播，并借助于高性能计算平台，开展了南加州区域高分辨率地壳波形成像研究，获得了该区域精细的3D地壳结构模型(图7a)。成像结果表明，和南加州地震中心提供的初始模型相比，局部非均匀性扰动可达30%。此外，一些尺度更小的结构在反演结果中也被清晰地揭示出来，结合南加州地区的地质背景和动力学分析，推测该地区地壳深层的异常区域可能是法拉隆板块俯冲导致的海洋地壳碎片。Lee等(2014)应用基于散射积分方法和伴随波场方法相结合的全三维层析成像改进了南加州地壳模型，揭示了强烈的地壳非均匀性，其结果可提高地震灾害评估的准确性。此外，Gao和Shen(2014)在Cascadia地区通过全波形背景噪声面波层析成像方法，揭示了火山弧后方地幔楔中的低速异常体，表明板块俯冲引起的小尺度对流影响了火山活动的分布。Wang等(2020)在Tape等(2009)的地壳模型基础上，基于背景噪声层析成像方法获得了径向各向异性剪切波速度模型(图7b)。反演结果显示，在圣安德烈斯断层的西面存在强的负各向异性，而前人研究结果并未发现这一特征，他们推测这是由于下地壳的斜长石与黑云母和橄榄石的各向异性特征不同而导致。Huang等(2016)关于1992年美国南加州Landers地震(M_w7。3)区FWI成像结果显示岩浆扩张中心在索尔顿海槽附近呈北东方向扩张，且存在扩张边界；发现研究区中下地壳中的高P波-低S波速度异常，这可能与地幔中残存的古法拉隆板块脱水和下沉驱动有关，揭示了古板块消亡过程对上地幔和地壳速度结构和Moho面形态的影响(图8)。

图7南加州地区不同深度的剪切波速度模型（a）和径向各向异性模型（b)

     剪切波速度模型修改自Tape等（2009）中的图4，径向各向异性模型修改自Wang等（2020）中的图3.

图8 使用直达波和Moho面反射波的多震相层析成像结果

（a）、（b）和（c）分别是南加州Landers地震区反演得到的P波、S波和泊松比纵向剖面图，其中BBE和LE分别代表BigBear地震和Landers地震，SJF和SAF代表圣哈辛托断层和圣安德烈斯断层。修改自Huang等（2016）中的图4.

Fichtner等（2009，2010）开展了澳大利亚地区上地幔速度和各向异性结构研究，获得了高分辨的地震成像结果。Zhu等（2012，2013）通过FWI成像获得了欧洲地区的高分辨率地壳地幔模型。他们的结果揭示了许多较小尺度的结构(包括俯冲板块、上升流等)，特别是存留在地幔过渡带中的滞留板块，为加深该区域深部结构和动力学演化过程的认识提供了重要约束(图9)。

图9 欧洲板块层析成像结果

（a）~（g）表示不同测线的垂直剖面：（a）亚得里亚-迪纳里德斯板块和弗朗西亚板块，（b）希腊板块，（c）阿尔卑斯板块俯冲，（d）纳维亚地区岩石圈拆沉，（e）埃菲尔热点及其伴生的低波速储层，（f）卡拉布里亚板块拆离，（g）北非板块拆离。修改自Zhu等（2012）中的图3

Chen等(2014)利用背景噪声伴随层析成像方法获得了青藏高原东缘精细的三维地壳模型，给出了该地区低速带的空间展布，揭示了该区域地壳流特征远比地壳通道流模型所预测的复杂。Chen等（2015）基于1869个台站记录到的227个地震事件数据，利用FWI成像方法对东亚壳幔结构开展反演成像研究，在经过800万个机时的计算后，获得了东亚地区径向各向异性模型（EARA2014）。该模型显示，浅部区域与地表地质构造相吻合，而深部结构异常揭示了太平洋俯冲板块在地幔过渡带中的形态和动力学过程。

在中国，Zhang等（2018）将背景噪声波形成像方法应用到华北克拉通地区的密集线性台阵，获得了比基于大圆路径和1D分层模型面波成像更高分辨率的结果，显示出华北地壳S波速度存在横向不均匀性的特征。其中东部的渤海湾盆地地壳薄、速度低，存在连续的正径向各向异性，反映了东部克拉通在中生代晚期经历的大规模伸展变形和岩浆作用；西部的鄂尔多斯盆地地壳厚、速度较高，径向各向异性较弱，表明中生代以来的构造事件对西部克拉通没有显著的影响。Xiao等（2020）基于ChinArray台阵的1335个密集流动台站和246个固定台站收集到的地震波形数据，采用伴随层析成像方法获得了青藏高原高分辨率三维径向各向异性速度结构，揭示了青藏高原西部印度岩石圈以低角度俯冲至帕米尔高原，而东部则以近垂直高角度进入地幔转换带，同时拉萨地块中部岩石圈的显著减薄暗示该区域可能发生了地幔对流剥离作用。Dong和Yang（2020）将基于最优传输理论（Villani，2003）的FWI方法应用于实际地震数据，获得华北地区岩石圈径向各向异性新模型。与前人结果相比，新模型揭示渤海盆地结晶地壳表现出更连续的高速特征。Tang等（2023）基于改进的FWI方法给出了川滇地区3D地壳-上地幔的多尺度剪切波衰减模型。结果表明，大地震通常发生于具有强衰减或衰减过渡带区域。Dong等（2019）的速度模型与Tang等（2023）的衰减模型所揭示的该区结构特征一致，从不同角度（速度和衰减结构）为青藏高原东南缘地区地壳通道流假说提供了新的证据。

FWI除了被广泛应用于全球各个局部区域的地球深部构造的成像中，也被应用于全球尺度的地震成像。French等（2013）利用FWI成像方法研究全球地幔结构，发现在上地幔200~350km深度存在与绝对板块运动方向近乎平行的低S波速带，延伸数千公里。Bozdağ等（2016）利用分布于全球的253个5级以上地震以及分布于全球各地的台站所记录到的380万条三分量地震数据，应用FWI伴随成像方法，构建了全球尺度的速度模型（图10）。该模型通过利用不同频率的体波和面波数据，在迭代的过程中不断提升所用数据的频率，从而不断提升模型的分辨率。经过15次迭代，获得了全球尺度高分辨率层析成像模型（GLAD-M15）（图10（右））。从该模型可以清晰地分辨出已知的俯冲板块、地幔热柱和热点等，为全球板块构造和地幔对流的认识提供了新的见解。Lei等（2020）使用震级在5.5~7.2之间的1480个事件进行FWI成像，在GLAD-M15模型的基础上经过10次拟牛顿算法迭代，获得了更高精度的全球层析成像模型（GLAD-M25），发现了更多新的地幔柱和俯冲带模式。

图10 两种模型在250km深度下SV波速度结构对比

左侧为已有S362ANI模型，右侧为使用FWI成像获得的全球尺度层析成像模型GLAD-M15.与S362ANI模型相比，FWI的成像结果进一步揭示了一些诸如俯冲板块、地幔热柱、热点等的主要构造体。修改自Bozdağ等（2016）中的图8和9

3.3 地震孕育和发生

FWI还被应用于强震区和典型断裂带的精细结构探测，以深化对地震孕育和发生机制的认识。Liu等(2023)利用差分伴随噪声层析成像揭示了洛杉矶盆地含水结构，并发现了与地震活动相关的高孔隙度区域；Yan等(2024)采用线性密集台阵的双差伴随背景噪声面波波形成像，获得了郯庐断裂带潍坊段上地壳的精细二维速度结构模型。相较于传统的背景噪声面波走时成像，该研究更好地揭示了断裂破碎带的低速结构异常，为理解断裂带的孕震构造和地震动评估提供了重要的参考模型。Dong和Yang(2022)基于青藏高原东缘密集的地震事件和观测台站数据，利用多尺度FWI成像方法获得了该区域高分辨率的S波速度、径向各向异性和泊松比结构。结果表明，该区域内强震发震位置多位于速度和径向各向异性高梯度区或高、低梯度值过渡带区，表明应力优先在这些构造异常区积累，使得这些区域易于孕育和发生地震(图11)。此外，研究发现2022年泸定M_s6.8地震的震源区下方存在显著的低S波速、负径向各向异性(V_SHSV)和高泊松比异常区(图12)，其可能与深部地壳和上地幔流体(包括部分熔融)活动有关。这些熔/流体沿倾斜的Moho面和/或壳内断裂向上迁移，致使中、上地壳孕震层的弱化和应力集中，进而导致泸定地震的发生，表明深部熔/流体迁移在泸定地震孕育和破裂过程中发挥着重要作用(董兴朋等，2023)。

图11 青藏高原东缘强震位置（五角星）与波速等参数之间的关系

（a）S波速度；（b）径向各向异性；（c）V_P/V_S值；（d）~（f）分别为S波速度、径向各项异性和V_P/V_S的梯度。修改自Dong和Yang（2022）中的图2.

图12 沿北纬29.6°的地表高程（a）、S波速度扰动（b）、径向各向异性（c）和泊松比（d）垂向切片图

横坐标为经度，白色五角星代表泸定地震；灰线代表Moho面埋深。修改自董兴朋等（2023）中的图3.

3.4 工程地球物理勘探

在前述的实际案例中，反演成像均聚焦于“千米”尺度的地震成像结果，但FWI反演方法同样可以应用于以“米”，甚至更小尺度的工程地球物理勘探领域，包括近地表结构成像（Amrouche和Yamanaka，2015；Dokter等，2017；Wang等，2017）、桥梁桩基检测（Tran等，2019）、铁路隧道检测（Yu等，2021）、板状和管道厚度检测（Rao等，2016）等。Jalinoos等（2017）将弹性波FWI方法应用于等间距的桥墩桩基反演成像，他们使用频带宽度为200~1200Hz的声波信号进行反演，获得的高分辨结果清晰地显示出桩基的形态和位置（如图13所示），说明FWI成像技术在工程地球物理勘探领域有着巨大的应用潜力。

图13 P波和S波合成模型

（a）~（d）依次是真实模型、初始模型、500Hz反演结果和1000Hz反演结果。修改自Jalinoos等（2017）中图7.

3.5 医学成像

在医学领域，FWI成像技术同样备受关注。Guasch等（2020）将FWI方法应用于大脑成像（图14），通过使用超声信号获得了“亚毫米”尺度的成像结果。基于FWI方法的超声成像技术可为中风快速诊断、头部创伤，以及神经系统疾病提供常规监测，未来有望替代传统的X光检查。这一超声成像技术可以提供更多与诊断相关的定量组织特性。与CT成像和X射线相比，传统的X射线和CT成像都具有一定的电离辐射，对人体会造成较大的危害，而超声FWI成像技术利用超声波震源和FWI方法在不损害人体组织的情况下即可实现大脑的高分辨率成像。因此，这种超声FWI成像更加安全可靠，具有很强的应用价值。此外，频率域超声全波形反演和深度学习驱动的超声全波形反演等方法也已逐步应用于医学成像研究中，提高了计算效率和成像结果的分辨率（Robins等，2021；Wu等，2023；周C等，2023）。

图14 头骨外脑部数据反演

（a）脑部真实模型；（b）初始模型；（c）超声波层析成像结果；（d）超声波FWI成像结果。红色椭圆表示换能器位置。可以看出，FWI反演结果既精确又有很高的空间分辨率。修改自Guasch等（2020）中的图2.

4  面临的挑战问题与重点研究内容

尽管FWI在过去三十年内获得了巨大发展，但不论是在FWI的理论方面还是在实际应用方面仍面临着许多重大挑战。这些极具挑战性的问题主要包括：(1)解的非唯一性问题；(2)计算和存储需求巨大的问题；(3)多参数反演耦合问题；(4)震相提取和匹配困难的问题；(5)数据处理智能化不高的问题。针对这些理论和实际的挑战，需要我们针对性地开展重点研究和攻关。

4.1 解的非唯一性问题

理论上，FWI是一个波动方程约束下的非线性优化问题。由于地震波波形数据与模型参数之间的非线性、震源信号的有限频性质和数据采集系统的不完备性，这一优化问题通常是欠定的，解是非唯一的。同时，基于L₂范数目标函数的FWI还存在由“周期跳跃”现象所导致的局部最优问题。这两种因素一并降低了FWI成像结果的可靠性。因此，解决FWI成像结果的非唯一性和优化求解过程中遇到的局部最优难题，扩大反演方法的有效收敛范围，提升反演结果的可靠性仍然是我们面临的最大挑战。

针对该问题，除了加强高密度实际地震观测之外，还需要从数据选取和反演理论两方面开展针对性的研究。一方面，需结合信号处理、图像处理以及人工智能的最新技术，实现多种震相数据的准确提取、低频有效信号的准确恢复、环境噪声的有效压制等目标，从而在数据层面更为有效地约束反演优化问题。另一方面，需开展更为深入的反演优化理论的研究和创新。通过引入凸性更强的目标函数、结合凸约束优化引入更为准确的先验模型信息（如带方向的Total Variation（TV）约束和基于Generative Adversarial Networks（GAN）的模型约束）、结合多物理机制进行联合反演（如地震电磁联合反演）等手段，以求解决解的非唯一性和局部最优的问题（Esser等，2018；Yang和马，2023年）。

近年来为解决FWI的多解性问题，人们将Wasserstein度量引入到目标函数的构建中，尽管基于该度量的FWI已被应用于油气勘探和被动源地震成像，但在Wasserstein度量理论及其在地球科学中的应用仍有待深入研究：

（1）现有的基于质量平衡的W₂度量、Kantorovich-Rubinstein（KR）度量以及Wasserstein-Fisher-Rao（WFR）度量等方法要求对正信号进行处理（Yang等，2018；Messud等，2021；周D等，2023）。然而，地震波信号具有实值、复值乃至向量值的复杂性。因此未来发展能够直接应用于真实地震波信号处理的新型度量方式，是一个极具挑战性的理论问题，亟需深入研究。

（2）与局部、逐点比较信号的L₂度量相比，因W₂度量是在全局层面比较两个信号，故基于W₂度量或WFR度量的计算过程涉及更大的计算量。这一特性要求针对W₂度量或WFR度量方法，发展新的快速计算方法，以提高计算效率。

(3)目前，已将具备解析解的一维最优传输映射应用到实际地震数据的成像研究中。然而，高维情况下的最优传输映射尚不能给出解析解，要获得这种情况下的映射，必须借助数值方法求解Monge-Ampère方程。因此，未来需要发展高精度快速数值求解Monge-Ampère方程的计算方法，以实现基于高维Wasserstein度量的FWI成像。

4.2 计算和存储需求的问题

FWI在优化求解过程中，需要反复求解波动方程以计算理论波形与目标函数梯度，这导致FWI的计算量巨大，成为FWI在大规模3D问题应用中的主要技术瓶颈，因此发展快速、稳定、高精度的正演波场模拟方法是克服FWI成像技术瓶颈的关键之一。

尽管频率域FWI成像具有处理多震源和刻画介质衰减性等的优点，但需要求解大型病态的线性代数方程组(Song和Williamson，1995)。求解这样的方程组，不仅计算非常耗时，而且存在数值求解方法的稳定性、收敛性等数学理论问题。因此，发展快速、稳定和精确地求解大型稀疏线性代数方程组的新方法是频率域FWI成像的主要挑战。

FWI在计算梯度与Hessian矩阵时，需要同时存储正演波场和伴随波场，这给计算机的存储性能提出了更高的要求。特别是针对三维FWI，存储三维时间域波场切片所需要的海量存储量也成为限制FWI大规模应用的瓶颈之一。

针对FWI面临的计算和存储需求巨大的问题，我们需要从软硬件两方面开展针对性研究。一方面，应继续开展基于包括有限元法、有限差分法、谱方法、间断有限元(DG)方法等在内的正演算法研究，以期实现在不规则粗网格下的地震波高精度快速数值模拟。对于频率域FWI问题，由于每一次正演模拟都需要求解一个大型稀疏线性代数方程组，导致频率域FWI的计算量巨大。为了实现包括Helmholtz方程在内的频率域波动方程组求解，需要发展具有良好数学性质(如系数矩阵具有比较小的带宽、对称、对角占优、良态或稀疏等)的空间离散方法和大型线性代数方程组的快速求解方法。此外，还应开展适用于求解FWI的随机优化算法研究。通过减少模型更新过程中波动方程的求解次数，可以显著加快FWI的计算速度。进一步地，在每次模型迭代时，施加适当的模型约束条件，可以在加速FWI的同时，保持反演精度。同时，为了缓解由伴随算法计算梯度向量所引起的存储量激增的问题，需要开展一系列基于检查点(Checkpoint)思想的动态波场存储技术，以时间换空间，通过额外的波场计算降低所需存储的波场切片数量，从而减轻存储需求。另一方面，高性能计算机技术的快速发展，为计算量巨大的FWI地震成像提供了快速计算的可能，也为FWI的快速实现开辟了新的途径。借助这些先进的高性能计算平台，并结合计算机系统和地震数据的特点，充分利用多核CPU和GPU并行计算能力，研究适用于FWI正反演的高性能并行计算方法，以提高FWI地震成像的计算效率。

4.3 多参数反演耦合问题和联合成像方法

随着对地球内部介质信息需求的不断提高，FWI逐渐从对单一的纵波速度反演发展到对多个弹性参数以及黏弹性参数的反演。多参数反演面临着多参数耦合的问题，即在反演过程中，各个参数之间会发生相互串扰，影响最终的反演结果。针对多参数反演耦合问题，一方面需要深入研究波形数据中各物理量与不同模型参数间的依赖关系，在设计反演问题时，尽量选择模型相互串扰弱的波形物理量-模型参数对；另一方面，二阶Hessian矩阵表征了各模型参数对各波形物理量的影响，通过构建Hessian矩阵并用其逆计算模型迭代量可有效减少各个参数之间的相互干扰。然而，二阶Hessian矩阵的构建相当复杂且计算量很大，因此研究构造近似的Hessian逆矩阵是解决多参数反演耦合问题的一个重要方向。

由于破坏性地震活动主要集中于地壳内，在特定研究区域内可能出现近震数据不足或地震活动性弱的情况，可以结合远震成像方法与区域台阵背景噪声面波成像方法(Yao等，2008；Gao和Shen，2014），实现多尺度、多震相、多数据的区域地震、远震和背景噪声面波数据联合反演，为地壳和岩石圈反演成像提供更多约束，进而提高反演结果的分辨率和可靠性。因此，结合各类传统反演方法，开展能有效提高FWI计算效率、成像精度和分辨率的联合反演成像算法的研究是未来研究的重要内容。

4.4 震相提取和匹配困难

如4.1所述，准确提取观测数据中各类有效震相数据，有助于对FWI优化问题进行合理的约束，缓解FWI解的非唯一性问题。然而，地球内部结构的非均匀性、岩石物性的复杂性以及地震观测系统的非一致性等因素导致地震信号呈现出高度复杂性，观测数据中不仅包含各种来源的噪声，而且各种不同震相高度混叠在一起，使得对不同震相的准确分离和拾取变得极其困难。另一方面，现有的地震波传播模型与实际地球介质中的波传播模型之间存在误差。地震波传播数学-物理模型(PDE)的建模方式基本基于守恒律、物理原理或者说唯像理论。其建模过程均是建立在各种简化假设和近似基础之上，很难构建准确刻画地震波在地球内部复杂介质中传播规律的定量模型。地震波传播的数学-物理模型的不准确和地球内部真实结构的复杂性导致地震波模型计算的理论波形无法和实际地震观测数据的多震相完美匹配。因此，震相的匹配问题成为了FWI成像的又一个关键问题和主要挑战。

针对震相提取和匹配困难的问题，一方面我们应该充分利用正演计算工具，对各类震相在不同地球介质内的传播规律开展深入研究，从而指导震相提取算法的研发。同时，我们可以结合现有的高质量实际地震数据和合成数据，开展基于深度学习和人工智能的智能震相提取研究。另一方面，针对震相匹配困难的问题，应进一步开展地震波传播的理论研究，建立更接近实际地质情况的地震波传播理论模型，为高精度地震学的发展奠定理论基础，进而获得更为接近真实情况的模拟数据。

4.5 人工智能技术的广泛应用    

近十年来，深度学习和人工智能技术在全球范围内被高度重视和广泛研究，并在许多行业获得了成功应用。未来，人工智能技术也将助力FWI技术的理论研究和实际应用。人工智能技术可以在如下几个方面助力FWI。

(1)正演计算：通过将深度神经网络与波动方程求解方法相结合，可以研发出非均匀网格下的高精度快速正演算法；

(2)模型先验知识：通过利用GAN等生成式模型，结合已知的一系列地球内部结构模型，可以建立全球地球内部结构先验知识库和所对应的生成式模型，将其应用于FWI可以更有效地对优化问题进行约束，降低成像结果对初始模型的强依赖性，缓解多解性和局部最优的问题；

（3）多震相的智能提取：通过目前已获得的各震相海量数据，结合卷积神经网络（Convolutional Neural Network，CNN）在内的各类深度学习模型，构造自动化分离提取地震信号中各震相的智能算法；

(4)低频信号的智能恢复：通过对实际数据和合成数据的学习，结合CNN在内的各类深度学习模型，构造自动恢复低频地震信号的智能算法；

(5)自动微分计算：在FWI中，各参数相对于目标函数的微分信息需要复杂的理论推导和计算。借助现今深度学习各类软件中的自动微分工具，可以极大地简化FWI中相应微分信息的推导和计算，加速研究的推进和应用的推广。

人工智能在FWI理论研究和实际应用中的潜力远不止上述方面，还需要进一步深入地探讨和研究。

5  FWI成像方法未来应用领域

5.1 4D地震成像

4D地震成像(又称为时移地震成像)指的是对同一地区的不同时间点进行三维成像来追踪地下结构的演变过程。该方法在资源勘探与开发、地下储气库监测、二氧化碳地质封存监测、大型水库安全监测以及自然环境监测与灾害预警等方面有着广泛的应用。现今地震数据采集技术和处理计算速度的显著提高，进一步推动了4D地震成像技术的迅猛发展。4D地震成像技术的一个主要应用领域是油藏生产过程的量化监测，通过差异分析和可视化技术描述油藏内部各储层参数(孔隙度、渗透率、饱和度、压力、温度等)的变化，提高油田的采收率，实现经济效益最大化。二氧化碳地下封存后的动态监测以及对周围环境的影响监测也离不开4D地震成像技术。此外，该项技术还可用于研究地震和火山触发的过程，通过对孕震区和火山深处地壳内由地下应力场、流变场变化造成的地震波速度的时空变化进行计算，我们可以实现对地震和岩浆活动进行监测。地震是由于地下应力的积累超过断层的强度所致，因此地壳发震区内随时间演变的应力是控制地震孕育和触发的最重要物理参数。通过实时反演地球内部物理参数变化情况，将有望成为预测地震发生、岩浆活动、火山触发的重要手段，揭示地球动力学机制和地球演化过程。因此，4D地震成像是未来FWI的重要应用领域。

5.2 地球动力学模拟

大尺度动力学数值模拟理论和方法在超级计算机发明之后得到了广泛应用和快速发展，成为研究区域乃至全球板块构造演化及其深部动力学过程不可或缺的重要手段之一。数值模拟的准确性在很大程度上依赖于地学观测数据，而地震成像结果是最具约束力的数据之一。与传统的地震学方法相比，FWI成像技术能够揭示更为精细的地球内部结构。一方面，在数值模型构建之初，上述精细的壳幔结构成像，结合地质学及板块重建的结果，能够为初始模型构建时岩石圈的形态、壳幔分层及流变特征、板块运动特征等的选择提供重要参考；同时又能有效地节省数值计算所需的经济和时间成本。另一方面，FWI的成像结果还能够为数值模型结果提供有效约束。动力学数值模型通常采用正交试验的方法测试不同控制参数的具体影响，因而，如何在一系列数值模型中选取最优者，为研究区域提供更为准确的动力学解释，亟需现今高精度的地下结构作为指导和验证。因此，将FWI应用于动力学数值模拟的研究，对于更深入理解地球内部时空演化过程具有重要意义。

5.3 行星成像

行星探测对于探索行星的起源和形成机制、地球生命起源、恒星演化和行星形成关系具有重要的科学意义。每个行星都展现出独特的物理、化学和矿物质特性，利用FWI技术对行星进行星震成像，获取行星内部的构造特征和物质组成及分布信息，对于揭示太阳系起源及演化等重大科学问题具有关键意义。行星内部蕴含着丰富的稀有金属和矿产资源，是潜在的“地外矿藏”，可供人类开发利用。随着相关技术的发展和成熟，行星将成为人类开发矿产资源等的下一个目的地。因此，将FWI技术应用于行星成像，拓展人类生存空间，是下一步需要关注和研究的重要科学问题。

5.4 医学和工程领域成像

FWI成像技术在医学领域和桥梁、轨道等工程领域应用也在逐步开展。相对于区域和全球尺度的地震成像，医学和工程领域的FWI成像尺度非常小，量级通常为毫米或亚毫米。使用主动源波形层析成像方法和高频的震源子波，如超声波等，对反演结果的分辨率也有更高的要求和限制。超声FWI成像具有低成本、便携、安全等特点，能够满足快速、精确、无损的医学诊断和轨道路基等工程检测。然而，在实际应用中FWI成像分辨率还达不到理论水准，与医学和工程领域实际要求的精度和分辨率还存在着差距，这需要地球物理学和数学专家与医学和工程领域专家进一步加强交流合作，促进学科交叉融合，共同推动FWI成像技术的基础理论发展和场景应用。

6  总结与展望

6.1 总结

以FWI为代表的高分辨地震成像方法和技术，是近20年来地球物理领域的研究热点。它不仅能对地球内部结构进行高精度成像，同时还可被应用于医学、轨道检测和桥梁工程等领域。目前，FWI成像技术已被越来越多地应用于局部、区域和全球尺度的地球内部构造成像研究中，并获得了一系列重要的发现，为深入研究地球内部结构、物质组成和分布、板块俯冲过程和地球动力学机制、地球演化过程、油气资源勘探等提供了新的技术手段。

本文对高分辨率成像，特别是FWI成像基本原理和研究进展进行了系统性的总结，并针对所存在的问题，提出了新的见解和解决思路。主要包括以下几个方面：

(1)非线性FWI理论所面临的挑战性问题主要是解的非唯一性、收敛性，以及反演目标函数的高度非凸性等。通过构建基于最优输运理论中的Wasserstein度量或其他类型目标函数，可以提升目标函数凸性，扩大反问题的收敛半径，进而提高反演成像结果的精度和可靠性。

(2)正演计算是FWI成像的基础，当前亟需发展时间域和频率域波动方程高效、高精度的数值求解方法，进而提高反演成像的计算效率，这是实现三高(高效、高精度、高分辨率)FWI成像的重要基础。

(3)采用多尺度、多震相和多数据融合的FWI反演策略是实现高分辨率成像的有效途径。尤其是多震相和多数据融合可为成像结果提供更多约束，确保反演结果的可靠性和精度。

(4)多方法交叉，充分结合各方法的优势。如：FWI和全局收敛方法的结合以获取可靠的成像结果，将FWI与射线类成像方法相结合以提升深部成像分辨率等，实现高精度地震成像。

(5)将FWI技术用于解决重大的科学问题，包括开展地球深浅部高分辨率结构和衰减与强震之间的关系，以及行星探测等研究，不仅能够加深我们对地球和其他星球内部结构与动力学过程的理解，还可为国家防震减灾和公共安全提供科技支撑。    

(6)以深度神经网络模型(深网模型)为代表的人工智能技术展现了巨大的应用潜力，并已逐渐应用到各学科中，未来将基于地震波传播模型的FWI成像与人工智能相结合，是实现高分辨率且快速地震成像的重要途径。

6.2 展望

在过去近40年里，地震层析成像技术作为揭示地球内部非均匀结构的重要工具，为人类理解地球内部结构和动力学过程、板块俯冲机制以及地震和火山活动机理等提供了宝贵的信息。随着对地球及其他星球探索的深入，地震层析成像预计将继续在未来的地球科学研究中发挥其不可替代的作用，特别是近20年来快速发展的高分辨率FWI成像技术，被视为实现“透明地球”目标的关键。然而，当前FWI在计算速度、反演问题的高度非线性和不适定性等方面仍然面临挑战。未来，围绕解决这些挑战并结合诸如具有“天然实验室”之称的青藏高原、中国地震危险区和地震多发区(如南北地震带、川滇地区等)重点研究区域，以及矿产、油气等地球资源问题，开展地壳上地幔结构成像和衰减结构的高精度高分辨率成像研究。

未来工作重点围绕以下目标展开：(1)构建中国大陆壳幔三维精细结构模型，为“玻璃地球”计划服务；(2)为地震预测、预警提供精确地球模型，认识强震孕育和发生过程；(3)为深部能源和矿产资源的勘测提供重要的高精度高分辨率“地图”；(4)探索其他星球(如月球、火星等)的内部结构，拓展人类生存空间；(5)为其他领域(如医学、工程地球物理、材料科学等)应用提供理论和技术支撑。

随着地球物理学、计算机科学、计算数学和信息科学等学科的进一步发展与交叉融合，FWI成像理论和方法不仅能够在更精细的时空尺度上发挥作用，而且在自然灾害预测、地球资源探测以及探索其他星球的内部结构和物质组成等方面具有重要的科学和实际意义，并推动地球物理学及相关学科的交叉研究，促进新理论和新学科方向的形成，进而对推动社会科技发展具有积极意义。

致谢 感谢责任编委和两位匿名审稿专家的宝贵意见与建议，这些意见极大地提升了本文质量。

-------END------    

原文来源：杨顶辉，董兴朋，黄建东，等. 高分辨率全波形地震成像研究——进展，挑战与展望[J]. 中国科学：地球科学，2025(2). 

封面标题、导读评论和排版整理：《覆盖区找矿》公众号.

推荐读者下载、阅读和引用原文！

------往期精彩回顾------

1. 姚华建等：面波走时三维结构成像方法综述、应用案例与发展方向

2. 6种面波层析成像方法及其在青藏高原的应用与对比

3. 地球科学：高分辨率地震层析成像显示青藏高原软流圈物质运动

---关注“覆盖区找矿”，拥有更多新方法---

宣讲成果，助力转化，激励创新！