PyFWI：地震全波形反演和油藏监测的Python软件包

软件X 22（2023）101384原始软件出版物PyFWI：用于全波形反演和油藏监测的Python软件包Amir Mardana，Bernard Girouxa，Gabriel Fabien-OuelletbaINRS-ETE，490 Rue de la Couronne，Québec，QC，G1K 9A9，Canadab蒙特利尔理工学院，2900 Boulevard Edouard-Montpetit，蒙特利尔，魁北克，H3 T 1 J 4，加拿大ar t i cl e i nf o文章历史记录：2022年11月30日收到收到修订版，2023年2月15日接受，2023年保留字：全波形反演储层监测时移地震PythonGPU编程a b st ra ct地震数据的全波形反演（FWI）是一种可用于对地下成像以及监测地下中的时移变化（TL-FWI）的技术。PyFWI是一个软件包，旨在在GPU上高效地执行FWI和TL-FWI，以用于研究目的。在PyFWI中实现了几种延时策略，如并行，双差，级联，中心差分，交叉更新，同时和加权平均。TL-FWI的一个重要挑战是不同年份的参数之间的为了缓解这个问题，PyFWI允许使用不同的参数化。PyFWI是用Python编写的，依赖于OpenCL来实现GPU上的计算使用OpenCL使PyFWI可以跨使用来自不同制造商的GPU构建的系统进行版权所有2023作者。由爱思唯尔公司出版这是CC BY许可下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）中找到。代码元数据当前代码版本0.1.9用于此代码版本的代码/存储库的永久链接https://github.com/ElsevierSoftwareX/SOFTX-D-22-00399Reproducible Capsule的永久链接法律代码许可证GPL-3.0使用git的代码版本控制系统使用python的软件代码语言、工具和服务、OpenCL编译要求、操作环境依赖性如果可用，链接到开发人员文档/手册https://pyfwi.readthedocs.io/en问题支持电子邮件mardan.amir. gmail.com1. 动机和意义地震方法被用来成像和表征地下。他们依靠机械波探测地球[1]。这些波基于地球的（粘）弹性特性传播。全波形反演（FWI）已显示出其直接从记录的波场中提取地下（粘）弹性特性的高分辨率图像的能力[2，3]。时间推移FWI（TL-FWI）也显示了监测地下性质变化的有前景的结果，这对CO2封存和石油天然气应用很重要[4*通讯作者。电子邮件地址：martan.amir. gmail.com（Amir Mardan），bernard. inrs.ca（Bernard Giroux），gabriel. polymtl.ca（Gabriel Fabien-Ouellet）。https://doi.org/10.1016/j.softx.2023.101384FWI允许通过最小化从初始模型估计的数据与接收器位置处的观测数据之间的残差来获得地下模型用局部梯度下降算法进行最小化。为了估计失配函数的梯度，使用伴随状态方法[8]，其需要多一次建模运行。由于FWI在计算上是昂贵的，已经做了大量的工作来使这种方法适应GPU编程[9，10]。PyFWI在其核心中依赖于OpenCL [11]，通过利用GPU编程来加速计算由于PyFWI是用Python编写的，因此这种集成是通过PyOpenCL进行的[12]。与CUDA相反，OpenCL与大多数现有处理器兼容[10]，这使得PyFWI成为一个可移植的选择。除了计算成本之外，FWI也是一个不适定问题[13]。另一个挑战是参数之间的串扰[3，14]。为了使问题更好地定态，本文介绍了四种常用的正则化方法，如Tikhonov，全变分，2352-7110/©2023作者。 由Elsevier B.V.出版。这是一篇开放获取的文章，使用CC BY许可证（http://creativecommons.org/licenses/by/4.0/）。可在ScienceDirect上获得目录列表SoftwareX期刊主页：www.elsevier.com/locate/softxAmir Mardan、Bernard Giroux和Gabriel Fabien-Ouellet软件X 22（2023）1013842阿：←阿克阿克拉阿拉拉阿克阿拉τzz但是，u布勒姆阿阿克拉阿阿克拉阿阿克拉阿拉阿克我IJ1= −，τ、⎧⎪⎪⎨ 阿勒特 =z 公司简介 +，拉斯阿阿克τxx先验信息和参数关系[4，15-解决串扰问题的最有效方法是使用最佳参数化[14，18]和使用Hessian优化[19，20]进行反演。为了模拟各向同性弹性地球，需要三个参数[18]。根据手头的问题，可以使用不同的参数化（都与弹性属性相关）来执行计算。在当前版本中，PyFWI可以用于在四个不同的参数中执行TL-FWI，对于边界条件，使用完全匹配层[27]（PML）来减少模型边缘的反射能量。为了估计偏微分方程约束的代价函数的梯度，最优化问题，如FWI，可以使用伴随状态方法[8]。波场的伴随，u，写为，u=[vx，vz，τxx，τzz，τxz]T，（ 4）是方程的伴随解（2），表示为诸如DV（密度、P波速度和S波速度）、PCS（孔隙度、粘土含量和含水饱和度）、LMD（孔隙度、粘土含量和(Lamε模量和密度）和KMD（体积模量、剪切模量和密度）。v −<$λτzz−<$µτxz−rvx，v<$z <$λτ<$xx<$（λ+2µ）τ<$zz<$µτ<$xzrvzlus和密度）。应该注意的是，用户可以轻松地增加此列表与自定义的参数化。除了除了上述技术，梯度缩放能力-⎨<$τ<$xx=−<$ρ−1v<$x−r、=− −阿拉阿克τxz（五）[21]《易经》中的卦卦，是卦的卦⎩⎪∂τ˜xz =−ρ−1vx −ρ−1vz−r、[22]这是为了更好地解决这个问题。crosstalk的问题PyFWI是一个多功能和灵活的软件包，可以执行TL-其中T 是转置运算符，r表示残差在地震数据的组成部分成本函数的梯度可以通过下式获得：FWI 已经提出了不同的方法来提供准确的使用TL-FWI的地下时间推移估计。PyFWI允许生物学家使用七个⟨∂F(m)⟩方法，如并行，双差[23]，级联[24]，交叉更新[25]，同步[25]，中心差分[5]，和加权平均反演[6]。这些方法在其是前向波场与伴随波场之间的零滞后互相关乘以散射矩阵，F（m）。当量（2）用于执行前向传播，更详细[26]。布勒姆它在时间（t）：0→T），其中T是记录2. 软件描述PyFWI是一个Python包，于2022年 1月在GNU通用公共许可证（GPL v3）下的Python包索引（PyPi）上首次发布源代码托管在GitHub上，并对贡献开放。2.1. 软件功能PyFWI的主要功能是提供所需的工具来执行FWI和TL-FWI与各种技术，考虑弹性介质中的地震波传播FWI是接收器位置处观测和估计波场之间残差的局部最小化考虑到正则化项（χREG），FWI问题的代价函数可以写为：时间相反，Eq. （5）执行反向传播，时间的倒数（t0T）。因此，前向波场u应该在求解方程时可用（5）使用等式（1）计算成本函数的梯度（六）、在PyFWI中，波动方程在时间域中求解，并且u应该在正演模拟的每个时间步存储，以便在获得相应时间步的伴随波场之后使用虽然这种技术很简单，但这种数据量的存储和输入/输出要求很高，需要更有效的策略。在[28]中讨论了不同的策略。在PyFWI中，采用检查点策略，其中在某些时间步长存储前向PyFWI允许用户指定应存储的前向波场的百分比由于其易于实现和与GPU编程的极大兼容性，有限差分法（FDM）被用于求解方程组。（2）和（5）。这些方程在空间上离散化χ（m）WRF（m）−2如图所示，使用交错网格[29]。凌晨1随着时间的推移，衍生品2016 -05 -2200：00：00（d）X2+XREG，（1）使用蛙跳法估计，其中速度和其中R在接收器的位置处对地震波场进行采样向量m和d分别表示模型参数和观测数据，以及作为数据失配的加权算子的Wd波场u是偏微分算子F（m）的解。在时域中，对于2D弹性情况，该算子可以写为ρvx=τxx+τxz+svx，ρvz <$τ zz <$τxzsvz2018年10月20日 =（λ+2µ）vx+λvz+sτx，（二它们的伴随以整数时间步长更新，应力及其伴随以半时间步长更新。图图1b示出了用于波传播的时间离散化方案2.2. 软件构架算法1描述了计算成本函数的梯度的过程（等式1）。（1））关于模型参数使用伴随状态方法。该算法是作为PyFWI的核心在一个名为wave_propagation的模块中实现的（图1）。2）的情况。此模块允许用户执行波传播（对于-∂τzz =λvx+（λ+2µ）vz +sτ，我爱你xz =µ（μvx +πvz），向前和向后）使用类Waveguerator。挥 手-传播器用于执行前向传播（步骤3-其中s是源函数，并且其中x和z方向上的粒子速度（vx和vz）加上法向应力（τxx和τzz）和剪切应力（τxz）构建波场矢量u。问题的初始条件由下式给出：vt=0=τt=0（三）共轭梯度法，提出了最小化mχ=、包括两种主要方法，正演模拟⎪⎧Amir Mardan、Bernard Giroux和Gabriel Fabien-Ouellet软件X 22（2023）1013843用于执行反向传播并计算梯度（步骤15为了使用PyFWI估计梯度，首先需要输入采集几何结构并定义参数，例如源和接收器的位置以及源函数。该步骤在模块采集的帮助下执行。Amir Mardan、Bernard Giroux和Gabriel Fabien-Ouellet软件X 22（2023）1013844Fig. 1. 用于离散化方程的网格的几何形状。（2）和（5）。(a)空间离散采用交错网格。用于（b）前向和（c）伴随建模的时间离散化，其中k是时间步长。图二、（a）PyFWI的结构。PyFWI的核心是wave_propagation，用于执行波传播和计算梯度。 该模块在fwi和tl_fwi中用于执行FWI和TL-FWI。（b）每个模块的应用程序，其中使用采集和处理来定义采集并在当前版本中执行增益。除了观测数据和初始模型之外，这些参数是用于计算梯度的输入。PyFWI为前向和后向波传播准备缓冲区。在下一步中，该软件调用OpenCL内核，该内核应用有限差分格式将计算卸载到GPU。前向波传播在GPU上执行在下一步骤中，在接收器的位置（d est）处对波形进行采样为了计算代价函数的梯度，调用OpenCL内核来执行反向传播，并使用伴随状态方法计算梯度。最后，将梯度传输到CPU，并在Python中执行最小化步骤。在步骤28估计梯度之后在算法1中，模块grad_switcher可用于提供根据不同参数化的成本函数的梯度。这种操作的计算成本可以忽略不计。因此，更改参数化是在CPU上执行的。在wave_propagation之上，PyFWI由模块fwi和tl_fwi组成，分别执行全波形反演和时移全波形反演。前向和后向传播是执行这些反演的本质。模块tl_fwi中的类TimeLapse是从类WaveLaparator派生的，因此受益于后者来模拟波的传播。 类TimeLapse也在其计算中使用来自模块fwi的类FWI。成本函数的最小化主要发生在FWI类中。优化过程可以使用模块fwi_tools中定义的正则化进行正则化。3. 说明性实例在本节中，我们将展示一个简单的合成示例，展示PyFWI的主要功能。在这个例子中，我们模拟了两个常见的场景，并使用MarmousiAmir Mardan、Bernard Giroux和Gabriel Fabien-Ouellet软件X 22（2023）1013845≤≤←←2˜˜图三. 本文以Marmousi模型为例进行研究。真实（a-c）基线，（d-f）监控模型，和（g-i）时间推移变化。从左到右，分别显示VP、VS和ρ算法1用于使用伴随状态方法计算梯度的伪代码，其中Ns、Nt、dest和检查点是源的数量、时间样本的数量、估计的数据和前向波场应当被存储的时间步长的列表。1：输入：观测数据（d obs）、初始模型（m0）和采集参数2：初始化OpenCL3：s14：whiles Nsdo5：t16：whilet Ntdo7：为vi调用kernel_updatev8：打电话 τij的内核更新9：如果t在检查点中，则10：存储vi和τij11：增量t12：增量s第13章：一夜情14：计算残差dest−dobs15：当1≤s≤Ns时，16：t←Nt17：当1≤t≤Nt时，储层中的VP和VS对应于它们各自的背景速度。这意味着在一个储层中速度随时间减小（模拟储层中的CO2注入），而在另一个储层中速度随时间增加（模拟石油生产）。图3g-i显示了真实的时间推移变化介绍了基线模型和监测模型模拟的地震资料及其差异如图4所示。延时数据（图4c）被放大了5、更好的视觉效果该图显示了在部署7个各向同性爆炸源的采集中，中间源的炮集。对于这些源，Ricker小波与采用30 Hz的峰值频率。我们还在模型周围使用PML，以避免来自模型边缘的反射为了执行TL-FWI，我们在多尺度策略中使用了截止频率为[10，20，25，35，40，55] Hz的低通巴特沃兹滤波器，每个频率迭代10次（列表1）。多尺度策略被提供以减轻周期跳跃的问题，并且包括执行从较低频率向较高频率的反转。本研究使用同步时间推移算法，通过考虑基线模型和监测模型之间差异的惩罚项，最大限度地减少不同年份的观测数据和估计数据之间的差异，χ（mb，mm）=λb<$RF（mb，ub）−db<$2+λm<$RF（mm，um）−dm<$218：如果t在检查点中，则19：读作vi和τij20：其他第21章：打电话 τij的内核更新2+λmm−mb<$2，2（七）22：调用kernel_updatev forvi第23章：打电话 τij的内核更新24：调用kernel_updatev forvi第25章：一夜情 梯度的kernel_updateg第26章：一夜情第27章：一夜情第28章：一夜情第29章：一夜情模型[30]，其是用于分析FWI算法的效率的众所周知的模型[4，6，25]。图图3a-c示出了真实基线模型的性质，3d-f真正的监视器模型Marmousi模型包含两个砂岩储层。为了创建基线模型，我们改变了其中λ和下标b和m表示每个项、基线和监控模型的正则化权重。1234567891011上市 一曰：通过创建对象TL并通过提供初始模型、从PyFWI开始。tl_fwiimportTimeLapseinv_freqs=[10，20，25，35，40，55]iterations=[10freqinv_freqs]#CreatingTL-FWIobjecttl=TimeLapse（db_obs，dm_obs，in pa，src，rec_loc，model_shape）#通过在初始化模块中添加数据来实现目标，从而实现时间间隔的调整dm，rms=tl（m0，iterations，inv_freqs，tl_method =“sim“，n_params =1，k_0=3，k_end =4）Amir Mardan、Bernard Giroux和Gabriel Fabien-Ouellet软件X 22（2023）1013846见图4。来自（a）基线、（b）监测模型和（c）时移差的地震数据。延时数据的幅度被放大了5倍，以便更好地可视化。图五. 使用（a-c）DV参数化和（d-f）PCS参数化估计的延时图像。从左到右，列显示了VP，VS，和ρ分别。PCS参数化导致在这项研究中的准确性提高了11.7%每个频率的迭代，以及用于执行多尺度反演的频率列表TL-FWI结果如图所示。5使用两种不同的参数化。图图5a-c是使用DV参数化获得的，而图5a-c是使用DV参数化获得的。使用PCS参数化估计5比较估计的时间推移变化（图。 5）真实的变化（图）。图3g-i虽然估计仍然有工件，PyFWI为用户提供了选择适合手头项目的参数化和策略的可能性4. 影响PyFWI可用于进行全波形反演和时移全波形反演。使用OpenCL允许此包通过利用大多数供应商的GPU的计算能力虽然作者致力于进一步开发该软件包，但当前版本只能用于2D问题。该软件包于2022年1月首次发布，根据PePy网站的统计，它已经被下载了13000次。作者在各种研究中使用了该软件包，并发表了不同的论文[6，31-这应该增加了作者研究小组之外的这个软件包的可见度。5. 结论本文介绍了PyFWI，它是一个用于进行全波形反演和时移全波形反演的Python包。为了减轻一些固有的挑战，对于这些问题，PyFWI提供了不同的选项，例如不同类型的正则化、多尺度策略、梯度缩放和不同的参数化。该软件包最重要的方面是能够使用不同的参数化进行反演，这对结果有显着影响。使用OpenCL使该软件包与大多数GPU制造商兼容。源代码托管在GitHub上，并开放供贡献和进一步开发。PyFWI也在Python包索引（PyPi）网站上提供，使用Python包安装程序（pip）可以轻松安装。竞合利益作者声明，他们没有已知的竞争性财务利益或个人关系，可能会影响本文报告的工作数据可用性数据将根据要求致谢该 项 目 得 到 了 NSERC 的 支 持 ， 加 拿 大 发 现 基 金 B 。 Giroux（RGPIN-2017-06215）。A.马尔丹还得到了加拿大SEG/Landmark奖学金的部分支持。引用[1] 吉 鲁 湾 ttcrpy ： 用 于 走 时 计 算 和 光 线 跟 踪 的 Python 包 。 SoftwareX2021;16：100834. http://dx.doi.org/10.1016/j.softx.2021的网站。100834[2]Maharramov M ，Biondi BL ， Meadows MA.时间 推移逆 理论及 其应 用。Geophysics 2016;81 （ 6 ） ： R485-501. http://dx.doi.org/10.1190/geo2016-0131.1.Amir Mardan、Bernard Giroux和Gabriel Fabien-Ouellet软件X 22（2023）1013847[3] 放大图片作者：J.时域粘弹性全波形反演Geophys J Int 2017;209（3）：1718-34. http://dx.doi.org/10.1093/gji/ggx110.[4] Asnaashari A，Brossier R，Garambois S，Audebert F，Thore P，Virieux J.使用先验模型的正则化全波形反演的时移地震成像：哪种策略？GeophysProspect 2015;63（1）：78网址：//dx.doi.org/10.1111/1365-2478.12176网站。[5] 周伟，卢明.中心差分时移四维地震全波形反演地球物理2021;86（2）：R161-72。http://dx.doi.org/10的网站。1190/geo2019-0834.1。[6] 放大图片作者：Mardan A，Giroux B，Fabien-Ouellet G.改进级联法时移地震全波形反演第二届EAGE地震反演会议。欧洲地球科学家工程师协会; 2022年，第1-5页。http://dx.doi.org/10.3997/2214-4609.202229003网站。[7] 塔兰托拉河声波近似中地震反射数据的反演。地球物理1984;49（8）：1259-66. http://dx.doi.org/10.1190/1的网站。一四四一七五四。[8] 普莱西斯RE. 的审查伴随态方法为计算具有地球物理应用的泛函的梯度。GeophysJInt2006;167（2）：495-503.http://dx.doi.org/10.1111/j.1365-246X.2006的网站。02978.x。[9] Shin J，Ha W，Jun H，Min D-J，Shin C.使用单个GPU卡进行3D拉普拉斯域全波形反演Comput Geosci 2014;67：1网址：//dx.doi.org/10.1016/j.cageo.2014.02.006网站。[10]放大图片作者：J.时域地震模拟利用OpenCL在非均匀计算平台上ComputGeosci 2017;100：142网址：//dx.doi.org/10.1016/j.cageo.2016.12.004网站。[11]石杰，五原D，石G。 OpenCL：异构计算系统的并行编程标准。Comput SciEng 2010;12（3）：66-73. http://dx.doi.org/10.1109/MCSE.2010.69网站。[12][10] Klöckner A，Pinto N，Lee Y，Catanzaro B，Ivanov P，Fasih A. PyCUDA和PyOpenCL：基于脚本的GPU运行时代码生成方法。并行计算2012;38（3）：157-74. http://dx.doi.org/10.1016/j.parco.2011的网站。九点零一分[13] 作者：J.勘探地球物理全波形反演综述Geophysics 2009;74（6）：WCC1-26. http://dx.doi.org/10.1190/1的网站。3238367。[14]Operto S，Gholami Y，Prieux V，Ribodetti A，Brossier R，Metivier L等人，多分量数据多参数全波形反演指南：从理论到实践。Lead Edge 2013;32（9）：1040-54. http://dx.doi.org/10.1190/tle32091040.1网站。[15]Anagaw AY，Sacchi MD.用全变分方法进行边缘保持地震成像。J GeophysEng2012;9（2）：138-46.[16]Esser E，Guasch L，van Leeuwen T，Aravkin AY，Herrmann FJ.全波形反演中的全SIAMJImagingSci2018;11（1）：376-406.http://dx.doi.org/10.1137/17M111328X网站。[17]胡Q，Innanen K.具有岩石物理约束的弹性全波形反演。第一届应用地球科学能 源 国 际 会 议 。 勘 探 地 球 物 理 学 家 协 会 ; 2021 年 ， 第 662-6 页 。http://dx.doi.org/10的网站。1190/segam2021-3594471.1.[18] 塔兰托拉河地震反射数据的非线性弹性反演策略地球物理1986;51（10）：1893-903. http://dx.doi.org/10.1190/1的网站。1442046[19]Métivier L，Brossier R，Virieux J，Operto S.全波形反演和截断牛顿法。SIAMJSciComput2013;35（2）：B401-37.http://dx.doi.org/10.1137/120877854网站。[20]杨P，Brossier R，Métivier L，Virieux J，Zhou W.粘声多参数全波形反演的时域预处理截断牛顿法。SIAM J Sci Comput 2018;40（4）：B1101网址：//dx.doi.org/10.1137/17M1126126网站。[21]Nocedal J，Wright SJ.数值优化。New York，NY：Springer; 2006，p. 664.http://dx.doi.org/10.1007/978-0-387-40065-5_1网站。[22]Lavoué F，Brossier R，Métivier L，Garambois S，Virieux J.通过多偏移距GPR数据的全波形反演进行二维介电常数和电导率成像：频域拟牛顿方法。Geophys J Int 2014;197（1）：248-68. http://dx.doi.org/10.1093/gji/ggt528网站。[23] 黄文忠，陈晓，陈文忠.时移井间地震数据的差分波形层析成像及其在天然气水合物生产监测中的应用。2004年：SEG技术计划扩展摘要。勘探地球物理学家协会; 2005年，第2323-6页。http://dx.doi.org/10.1190/1.1845221网站。[24] [10]杨文，杨文，杨文.用于改进成像和储层表征的时移数据的全波场反演在：SEG技术计划扩展摘要2012。勘探地球物理学家协会; 2012年，第1-6页。http://dx.doi.org/10.1190/segam2012-1043.1网站。[25] Maharramov M，Biondi B.基于全变分正则化的时移地震数据稳健联合全波形反演2014年，http：//dx.doi.org/10.48550/arXiv.1408.0645网站，arXiv：地球物理学。[26]放大图片作者：Mardan A，Giroux B，Fabien-Ouellet G.加权平均时移地震全 波 形 反 演 Geophysics 2023;88 （ 1 ） ： R25-38. http://dx.doi 的 网 站 。org/10.1190/geo2022-0090.1。[27]Berenger J-P.用于电磁波吸收的完美匹配层。J Comput Phys 1994;114（2）：185-200. http://dx.doi.org//jcph.1994.1159.[28]Nguyen BD，McMechan GA.二维弹性叠前逆时偏移中避免存储震源波场快照的五种方法Geophysics2015;80（1）：S1-S18.http://dx.doi.org/10.1190/geo2014-0014.1网站。[29]非均匀介质中的P-SV波传播;速度-应力有限差分法。地球物理1986;51（4）：889-901. http://dx.doi的网站。org/10.1190/1.1442147。[30]张 文 辉 ， 张 文 辉 .Marmousi ， 模 型 和 数 据 。 1990 ，http://dx.doi.org/10.3997/2214-4609.201411190。[31]放大图片作者：Mardan A，Giroux B，Fabien-Ouellet G.用于流体饱和度监测的时移全波形反演第二届EAGE地震反演会议。欧洲地球科学家工程师协会;2022年，第1-5页。http://dx.doi.org/10.3997/2214-4609.202210635网站。[32]放大图片作者：Mardan A，Giroux B，Fabien-Ouellet G.非重复性对时移全波形反演的影响。第二届EAGE地震反演会议。欧洲地球科学家工程师协会;2022年，第1-5页。http://dx.doi.org/10.3997/2214-4609.202211009网站。[33] 杨文伟，王文伟，王文伟.利用时移全波形反演直接监测第二届国际应用地球 科 学 会 议 能 源 扩 展 摘 要 。勘 探 地 球 物 理 学 家 协 会 ; 2022 年 ，http://dx.doi.org/10.1190/image2022-3746685.1。