CN113945994A - 使用有限差分模型进行高速多源加载和波场检索的方法 - Google Patents

使用有限差分模型进行高速多源加载和波场检索的方法 Download PDF

Info

Publication number
CN113945994A
CN113945994A CN202110734477.2A CN202110734477A CN113945994A CN 113945994 A CN113945994 A CN 113945994A CN 202110734477 A CN202110734477 A CN 202110734477A CN 113945994 A CN113945994 A CN 113945994A
Authority
CN
China
Prior art keywords
processing unit
grid
source
receiver
central processing
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202110734477.2A
Other languages
English (en)
Other versions
CN113945994B (zh
Inventor
张昌华
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China Petroleum and Chemical Corp
Original Assignee
China Petroleum and Chemical Corp
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China Petroleum and Chemical Corp filed Critical China Petroleum and Chemical Corp
Publication of CN113945994A publication Critical patent/CN113945994A/zh
Application granted granted Critical
Publication of CN113945994B publication Critical patent/CN113945994B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V20/00—Geomodelling in general
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V11/00—Prospecting or detecting by methods combining techniques covered by two or more of main groups G01V1/00 - G01V9/00
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/62—Physical property of subsurface
    • G01V2210/624—Reservoir parameters
    • G01V2210/6248—Pore pressure
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/67—Wave propagation modeling
    • G01V2210/679—Reverse-time modeling or coalescence modelling, i.e. starting from receivers

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Geophysics (AREA)
  • General Engineering & Computer Science (AREA)
  • Geometry (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

一种方法，用于在GPU上针对波场模拟使用有限差分模型来有效注入震源和检索接收器波场，其中震源和接收器不在数值网格点上，因而处于任意位置。为此，该方法采用GPU纹理内存来增加内存读取带宽，然后将震源定位在有限差分网格中的任意或模拟位置，并将它们扩展到新生成的网格点数上。

Description

使用有限差分模型进行高速多源加载和波场检索的方法
技术领域
本发明大体上涉及使用高速计算方法的地球物理勘测和地震数据处理，其中在具有带宽和存储限制的图形处理单元中采用了有限差分模型。
背景技术
1.概述
即使使用当今的高性能计算机器，也可能没有足够的计算能力来使用单个统一计算方案对一整片勘探区域进行建模。因此，本领域的技术人员将转而利用混合处理单元来实现高级数值编程解决方案。这些混合方法需要有效且实用的实现方式，这是因为这些单元的几个不同建模内核在地震勘探模型的不同区域(或子域)中运行，甚至可能在不同的有限差分网格分辨率下运行。
两种主要的处理单元是图形处理单元(GPU)和中央处理单元(CPU)。GPU以比同类CPU更低的价格提供高性能的计算能力，并采用各种硬件模块，如台式机或中心。GPU的主要目的是提供改进的数值性能，因此GPU和CPU之间的性能有着显著的不同。例如，一个特定的GPU可以在一秒内进行一万亿次浮点运算，而CPU的性能由每个插槽的内核数量(例如4核)和所述插槽的性能(即每秒时钟周期数)决定。
然而，GPU依靠线程级并行来隐藏较长的片外内存访问延迟，因此明智地利用片上内存资源(包括寄存器文件、共享内存和数据缓存)对应用程序的性能至关重要。相比之下，GPU的一个片上特性是它们可以并行地寻址不同的片上内存，从而增加内存带宽和延迟。因此，管理GPU片上内存资源对于应用程序的开发人员来说是一项非常重要的任务。更重要的是，由于GPU的不同代的片上内存资源不同，性能可移植性成为一项艰巨的挑战。然而，在本领域中众所周知，GPU对通用编程提供的支持非常有限，这是因为在GPU中执行的代码不能自分配内存，并且需要由主机CPU来完成此操作。
使用混合方法的这些复杂性意味着从CPU到GPU的编程很困难，并且需要在相同或不同时间生成高级算法。显然，这影响了本领域的技术人员在复杂介质中计算大规模地震波动方程模型和地震波传播模拟时可以在石油和天然气工业中执行的工作量。
2.现阶段的GPU架构
使用统一计算设备架构(CUDA)编程模型的C/C++语言中使用最广泛的高性能并行应用程序之一是NVIDIA的Tesla架构。NVIDIA架构的一个例子通常是一组流式多处理器，每个处理器都有自己的处理器和共享内存(用户管理的缓存)。这些处理器完全能够执行整数和单精度浮点运算，并带有用于双精度的附加内核。然而，多处理器中的所有流处理器都与未被硬件缓存的全局设备内存共享获取、控制和内存单元。
在这些架构中，CUDA将线程排列成线程块，可以在分配给所述线程块的任何共享内存位置处读写。因此，线程块内的线程可以通过共享内存进行通信或使用共享内存作为用户管理的缓存，这是因为共享内存的延迟比全局内存的延迟低两个数量级。
然而，由于在大多数架构中前向和后向阶段访问和修改相同的数据，因此混合CPU/GPU实现需要在主机和GPU内存之间进行连续的内存传输，以保持数据的一致性。由于频繁和大量的内存传输，当前的GPU性能受到极大的影响，然后在GPU内核中执行步骤，以访问驻留在图形双倍数据速率(GDDRAM)模块中的数据。这提供了一个额外的好处，因为GDDR模块可以在同一时钟周期上请求和接收数据。因此，主机只需要计算执行环境、内核调用，以及在必要时往来于磁盘的I/O传输。
3.声波传播
声波的计算建模一直是一个活跃的研究领域，并且已经从航空航天等领域发展到地震学、地球物理学甚至气象学。事实上，空间阵列中声波的计算建模和模拟是许多科学和工程应用的基础，在石油和天然气壳层和储层勘探中也是一样。然而，迄今为止，在大型或复杂结构中基于声波传播模型来获得良好的结果仍然是一个主要的计算挑战。
该领域中最常用的模拟声波传播的方法通常包括两个阶段的过程：(1)计算代表声学空间的脉冲响应；以及(2)将脉冲响应与无回声记录的或合成生成的源信号卷积。然而，由于这些脉冲响应依赖于对时变空间声场建模的精确计算，因此它们的主要挑战仍然是精确求解方法(例如有限差分)的计算要求。在大多数情况下，用于求解声学传播的方法要求来自空间离散化的每个波长大约3-5个节点(即网格点)。假设最小频率范围为3Hz，最大频率为100Hz，节点间距为5-50m，结果是10立方公里的声学空间充满了约8百万到80亿个节点。因此，声波方程的数值求解器仅限于相当简单、受限或小的空间，这通常被认为是耗时的。为了克服这个问题，大多数当前的声学仿真系统使用几何声学技术，它们仅对较高频率和早期反射来说是准确的，并且可能无法对衍射效应进行建模。可以看到，克服了一个问题会产生另一个问题。
4.全波场反演(FWI)
FWI是一种偏微分方程约束优化方法，它迭代地最小化测量波场和计算波场之间的失配范数。地震FWI涉及多次迭代，单次迭代可能涉及以下计算：(1)求解前向方程；(2)求解伴随方程；以及(3)对这些前向方程和伴随方程的解进行卷积，以产生成本函数的梯度。应当注意，对于二阶优化方法，例如Gauss-Newton方法来说，还需要进行(4)求解扰动前向方程。
经典时域FWI最初由Tarantola于1984年提出，通过最小二乘意义上的预测数据和观测数据之间的差异中的能量最小化来改进速度模型(参见Tarantola,A.,1984,Inversion of seismic reflection data in the acoustic approximation,Geophysics,49,1259-1266，以及Symes,W.W.,2008,Migration velocity analysis andwaveform inversion,Geophysical Prospecting,56,765-790)。Tarantola于1986年对其进行了进一步开发，并应用于弹性情况(参见Pica,A.,J.Diet和A.Tarantola,1990，Nonlinear inversion of Peace in alaterally invariant medium,Geophysics,55,284-292)。此后，频域FWI成为一个活跃的研究领域，并为鲁棒反演提供了层次框架(参见Pratt,G.,C.Shin,et al.,1998,Gauss-newton and full newton methods infrequency-space seismic waveform inversion,Geophysical Journal International,133,341-362)。最近，Shin和Cha引入了拉普拉斯域FWI和拉普拉斯傅立叶域变体(参见Shin,C.和YH Cha,2008,Waveform inversion in the Laplace domain,GeophysicalJournal International,173,922-931,2009；以及Waveform inversion in the Laplace-Fourier domain,Geophysical Journal International,177,1067-1079)。然而，它仍然是一个具有挑战性的问题，并吸引了地球物理学家越来越多的努力(参见Virieux,J.和S.Operto,2009,An overview of full-waveform inversion in explorationgeophysics,Geophysics,74,WCC1-WCC26)。
然而，最近通过利用GPU在计算能力和硬件方面取得的进步使得FWI方法成为一种广泛使用的方法，以减轻地震成像中的计算缺陷(参见Micikevicius,P.,2009,3D finitedifference computation on GPUs using CUDA:Proceedings of 2nd Workshop onGeneral Purpose Processing on Graphics Processing Units,ACM,79-84)和反演(见Boonyasiriwat,C.,G.Zhan,M.Hadwiger,M.Srinivasan和G.Schuster,2010,Multisourcereverse-time migration and full-waveform inversion on a GPGPU,第72届EAGE会议及展览)。但是，GPU计算的一个关键问题仍然存在，它不一定与GPU的计算资源有关，而是与主机和客户端之间传输数据的通信时间以及它们的内存资源有关。
5.逆时偏移(RTM)
逆时偏移是本领域众所周知的地震成像偏移方法，它主要用于利用记录的地震波形来绘制地下反射率。RTM能够处理具有高速度对比的结构倾角与盐盆地和其他具有复杂结构和速度分布的地质盆地中常见的条件的最恶劣的组合。这些地下条件会严重扭曲任何所传播的波前和图像，因此使地震数据的建模和成像复杂化。因此，在这些类型的条件下，需要从多个震源进行炮击，从而创建可由传感器捕获的波场类型的宽广阵列。
RTM方法的另一个优点包括在与形成图像的位置相同的空间位置评估图像的质量。考虑到当今大多数地震偏移应用仍然使用初级反射作为唯一信号，RTM使用所有可计算波类型的能力是独一无二的，有助于减少由于将非初级波误认为初级反射而导致的成像伪影。
RTM的基本思想是一个三步过程：(1)通过适当的速度模型对波场进行正向建模；(b)通过完全相同的模型反向传播测量数据；以及(c)使用成像条件进行超级定位。完整的波动方程用于模拟。由于使用了包括多次反射在内的整个波场，因此成像没有倾角限制。因此，甚至可以对目标的垂直边界或下边缘进行成像。
在大多数情况下，RTM方法根据估计的速度模型和源参数在时间反转波形数据的外推和模拟预测之间寻找地下反射率的初始图像，作为图像空间中的最佳匹配。在数学意义上，RTM方法基本上是在时间上向后外推初始为零的x-y-z体积，将地震数据P(x,y,z＝0,t)作为每个时间步长处z＝0的边界条件，以便计算x-y-z体积在不同时间的快照。在时间t＝0，这个x-y-z体积生成图像I(x,y,z)的偏移结果。
6.有限差分模型
有限差分模型是石油和天然气行业中广泛使用的一种方法，用于复杂介质中的大规模地震波动方程建模和地震波传播模拟，用于合成数据的建模和偏移(如逆时偏移)或反演(例如全波形反演)。在这个模型中，地震波场的时间演化是在离散网格点的体积上采样的。然而，在这些网格点上不需要这种波场的激发波场或记录波场。
常见的商业应用包括参考合成地震数据生成、逆时偏移和全波形反演。商业的有限差分建模器通常在跨模型域应用公共内核的方案中求解声学或各向异性伪声波方程。这些方案通常是高度优化的，包括负载平衡，这是需要考虑吸收某些域中的边界区域的额外成本的。然而，计算硬件的可用性正在导致勘探区域中更大且日益异构的集群环境。这部分是因为多核CPU的变化速度足够快，集群通常包含跨越几代CPU技术的节点，还因为混合CPU和GPU集群以及混合内核的出现。因此，计算节点的这种非一致性使得负载平衡更难以实现。
但是，在使用有限差分模型时，无论是在CPU还是在GPU上，内存访问的成本通常是性能的瓶颈。例如，对于波场f的网格点i处的一阶空间导数
Figure BDA0003139882820000051
的数值计算来说，必须从计算机内存中加载f的N+1个值，导致N+1:1的读取-计算比，其中N是有限差分的阶数。因此，减少内存访问的次数或增加内存读写的带宽对于有限差分法的效率至关重要。
当激励源点处于有限差分网格点上时，注入一个这样的源只需要一次读和一次写的内存访问。但是，当源位置不在网格点上时，必须找到与位于网格点之外的源等效的网格点源组合。此网格点源的有效数量通常为nxnynz，其中nx(y,z)是x(y,z)方向上源位置的网格点数量，通常取等于精度要求的有限差分阶数N的一半的相同值。因此，对于16阶有限差分来说，如果nx(y,z)＝N/2＝8，则组合网格点源的数量最多为512个。更重要的是，这些网格点不在连续内存中，并且内存访问成本至少增加nxnynz倍。在逆时偏移(RTM)或全波形反演(FWI)中，数据成为反向波传播过程中的源。在这种情况下，有效网格点源的数量为nxnynzNxNy，其中Nx和Ny分别是x方向和y方向上的接收器数量。类似地，为了在任意位置检索波场，需要nxnynz个网格点上的波场值的组合。
在GPU上，与缓冲的片上内存访问相比，片外全局内存访问非常广泛。因此，当源或接收器不在网格点上时，如何增加GPU上的源加载或波场检索的带宽对有限差分码的性能有显著影响。
由于上述原因，需要找到一种更好的方法来增加源加载或波场检索过程的带宽，以在GPU上增加可用计算机资源之间的计算负载。
发明内容
因此，本发明通过提供一种使用有限差分模型的高速多源加载和波场检索的计算机实现方法，满足了上述需要并克服了现有技术中的一项或多项缺陷。
通常，根据本发明的一个实施例，空间局部性表明线程可能在空间上读取其他附近线程的内存地址。然而，在有限差分方案中，GPU线程仅被映射以处理网格点。类似地，用于源加载表达式(4)或波场检索表达式(5)的线程将仅访问其附近线程的内存。如果上述示例的内存位于片外全局内存资源上，则此类访问将成为执行成本高昂的操作。因此，通过使用根据本发明的实施例的纹理内存，线程可以更高的带宽从附近的内存地址读取，从而以10的倍量级大大减少了加载和检索时间。
根据本发明的一个实施例，在GPU上在有限差分法中于任意位置处快速多源加载的步骤序列通常包括：(1)使用nx(y,z)来确定每个实际源的组合网格点的有效数量，其中nx(y,z)通常取N/2，N是有限差分阶数；(2)将每个有效网格点源项视为具有权重的独立点源；(3)在片外全局内存上分配独立的网格点源网格位置、权重以及波场变量；(4)将那些全局内存映射到纹理内存；(5)使用CUDA编程的纹理内存获取函数和API来读取网格点源网格位置，以及这些网格点上的字段值和权重；以及(6)在该网格站点的字段上加载具有权重的网格点源。
类似地，本发明的另一个实施例需要在GPU上使用有限差分方法在任意位置处进行快速波场检索的步骤序列。该序列通常包括：(1)使用nx(y,z)来确定每个接收器的组合网格点的有效数量，其中nx(y,z)通常取N/2，N是有限差分阶数；(2)计算贡献权重；(3)在GPU的片外全局内存上分配这些网格点位置、权重以及波场变量；(4)将那些全局内存映射到纹理内存；(5)使用CUDA编程的纹理内存获取函数和过程来读取这些网格点上的字段值和权重；(6)将加权波场相加而得到检索波场。
最后，虽然这里公开的实施例描述了在利用标量声波方程的有限差分近似的波传播中的多源注入和波场检索的方法和系统，但是本领域的普通技术人员将清楚，它们还可以替代地应用于各向同性和各向异性介质中的矢量波动方程和弹性方程，而不会背离如下面的权利要求所定义的本发明的真实范围。
符号列表
Figure BDA0003139882820000071
Figure BDA0003139882820000081
附图说明
通过结合附图并考虑以下描述，可以容易地理解本发明的教导。
图1示意性显示了根据本发明的某些实施例的说明性环境的横截面图的示意图，其包含地震源、地震接收器、井眼位置、井筒、各种传输射线的入射点，以及各种入射角；
图2是显示了根据本发明的某些实施例的勘探区域的俯视图的示意图，其中示出了具有接收器和炮点网格点和线的采集形式，其中接收器和炮点均位于网格上的任意位置处；
图3是根据本发明的某些实施例的在网格上的任意位置处使用有限差分模型来高速加载震源入射位置点的计算机实现方法的流程图；
图4是根据本发明的某些实施例的在网格上的任意位置处使用有限差分模型来高速检索接收器位置的计算机实现方法的流程图；以及
图5示出了根据本发明的某些实施例的计算机实现方法的图形用户界面，显示了在不同的时间步长间隔处的压力波场，以及其各自的震源点和接收器点的位置。
具体实施方式
现在将详细参考本发明的几个实施例，其示例在附图中示出。应注意的是，在适用之处，在附图中使用相同或类似的附图标记来指示相同或类似的功能。附图仅出于说明的目的描绘了本发明的实施例。本领域的技术人员将从以下描述中容易地认识到，在不脱离本文描述的原理的情况下，可以采用其中示出的结构、系统和方法的替代实施例。
本发明的计算机实现方法引入了一种有效的方式，以便在震源或接收器不在数值网格点上时使用有限差分来来注入震源和检索波场，用于在GPU上进行波场模拟。在现有技术中，与计算指令相比，内存访问是一个瓶颈。此外，当震源位于有限差分网格中的任意位置时，它们必须在多个网格点上扩展，从而有效地增加注入到模拟中的震源点的数量。由于这些扩展的网格点往往不在一个连续内存区域中，因此这大大增加了内存访问成本，但当震源的实际数量很少时，所增加的内存访问成本是有限的。然而，在执行逆时偏移(RTM)或全波形反演(FWI)时，在反向波模拟的过程中实际震源的数量等于数据轨迹的数量，从而显著增加了内存访问成本。类似地，在前向波模拟中，当接收器不在网格点上时，在接收器处记录的波场必须从其相邻的网格点中获得。在本发明的当前实施例之一中，该计算机实现方法旨在通过使用GPU纹理内存来克服这些内存访问成本。尽管用于震源注入或波场检索的GPU线程必须读取非连续的内存，然而内存访问模式表现出大量的空间局部性。在计算应用程序中，这粗略地意味着线程可能从附近线程所读取的地址“附近”的地址中进行读取。GPU纹理内存设计用于增加此空间局部性内存读取的有效带宽。通过使用这种纹理内存，可以将震源注入和波场检索的效率提高十倍。
1.部署增加GPU带宽的技术
为了通过本计算机实现方法来计算声波传播方程，以简化形式开发了以下算法：
Figure BDA0003139882820000091
其中，VP是声波速度，
Figure BDA0003139882820000092
是因
Figure BDA0003139882820000093
处的点激励处产生的
Figure BDA0003139882820000095
处的压力波场，而s(t)是震源时间小波。
在有限差分数值解中，方程(1)在空间和时间上离散，X、Y、Z和t维度上的网格大小分别为(Nx,Ny,Nz,Nt)。然而，压力波场只在这些网格点上有值：
Figure BDA0003139882820000094
因此，在针对
Figure BDA0003139882820000101
的该有限差分模型中，从内存资源中读取和写入内存资源的计算始终是个瓶颈。
当震源位置
Figure BDA0003139882820000102
在网格点(is,js,ks)上时，有限差分模型中时间步长值为n下的震源加载可简单计算为：
Figure BDA0003139882820000103
对于每个震源来说，此震源加载的内存成本只是一次读取和一次写入内存资源。因此，与其他操作相比，本领域的大多数技术人员会忽略该成本。
A.当接收器位置处于网格点上时
另一方面，当接收器位置
Figure BDA0003139882820000104
在网格点(ir,jr,kr)上时，该接收器在时间步长值n下记录的压力波场
Figure BDA0003139882820000105
可以简单地检索为：
Figure BDA0003139882820000106
尽管如此，对于每个接收器来说，这一检索过程的内存成本只是在内存中的一次读取和一次写入。因此，与其他操作相比，本领域的大多数技术人员也会忽略该成本。
然而，当源点
Figure BDA0003139882820000107
不在有限差分网格点上时，必须找到一组网格点(i,j,k)，其组合起来等价于位于
Figure BDA0003139882820000108
处的震源。因此，必须实现以下等式：
Figure BDA0003139882820000109
其中，
Figure BDA00031398828200001010
以及
Figure BDA00031398828200001011
Figure BDA00031398828200001012
和
Figure BDA00031398828200001013
是针对任何网格外震源的最近网格点，
Figure BDA00031398828200001014
Figure BDA00031398828200001015
和
Figure BDA00031398828200001016
是权重函数，而Δx(Δy,Δz)是x(y,z)方向的网格步长。
理论上，每个实际点震源的网格点(i,j,k)的集合通常扩展到整个空间上的网格空间。然后，根据表达式(2)来加载表达式(4)中的每个网格点震源，这显著增加了成本。但在实际操作中，出于精度考虑，使用了表达式(4)的窗口网格点震源，使得对于
Figure BDA00031398828200001017
处的每个点震源使用：-nx/2+1≤i≤nx/2,-ny/2+1≤j≤ny/2,-nz/2+1≤k≤nz/2组网格点震源。因此，前面提到的网格点源的有效数量是nxnynz。同样，出于精度考虑，在实践中，2nx(y,z)通常取有限差分算子的阶数N，因此，对于16阶有限差分来说，nx(y,z)＝8，组合网格点震源项的有效数量最终高达512，这显然将加载成本增加了512倍。此外，由于对于每个点震源来说扩展网格点-nx/2+1≤i≤nx/2,-ny/2+1≤j≤ny/2,-nz/2+1≤k≤nz/2不在连续的内存资源中，因此内存访问成本也进一步增加。这种增加的内存访问成本在RTM或FWI的反向波模拟阶段变得更糟，这是因为实际点震源有NxNy个，而扩展网格点震源的数量变为nxnynzNxNy个。
B.当接收器位置不在网格点上时
类似地，当接收器
Figure BDA0003139882820000111
不在有限差分网格点上时，必须找到一组网格点波场项(i,j,k)，其可组合起来以得到位于
Figure BDA0003139882820000112
的波场，可使用以下表达式：
Figure BDA0003139882820000113
在表达式(5)下，接收器处的波场R(xr,yr,zr,n)与
Figure BDA0003139882820000114
的分散在不连续内存网格点-nx/2+1≤i≤nx/2,-ny/2+1≤j≤ny/2,-nz/2+1≤k≤nz/2中的nxnynz个网格点波场值组合起来，因此也导致了昂贵的内存访问成本。
虽然对于在震源或接收器不在网格点处或不在计算机内存中的连续区域上的情况下的震源加载或波场检索来说需要网格点：-nx/2+1≤i≤nx/2,-ny/2+1≤j≤ny/2,-nz/2+1≤k≤nz/2，然而如果从震源或接收器位置中扩展网格点，它们仍然可以在空间上定位。这可以通过GPU上的纹理缓存内存来实现，其最初是为图形应用程序设计的，其中内存访问模式表现出大量的空间局部性，其有效读写带宽比片外全局内存的读取高约10倍。此外，这种纹理内存也可用于通用计算，但到目前为止还没有被利用。
2.加载或注入震源或入射点
一般而言，本发明的实施例之一中提出的用于在GPU上运行的有限差分方法中在任意位置处进行快速多源注入的步骤序列包括：
a)确定各实际震源的组合网格点震源的有效数量。当震源点
Figure BDA0003139882820000115
不在有限差分网格点上时，必须找到一组网格点(i,j,k)震源，它们组合起来等价于位于
Figure BDA0003139882820000116
处的震源，这可根据表达式(4)和
Figure BDA0003139882820000117
其为网格外震源的最接近网格点，Δx(Δy,Δz)是x(y,z)方向的网格步长，
Figure BDA0003139882820000118
Figure BDA0003139882820000119
和
Figure BDA00031398828200001110
是权重函数，在结合于计算机实现方法中时用如下算法表示：
Figure BDA0003139882820000121
其针对具有坐标
Figure BDA0003139882820000122
的网点点i，满足：
Figure BDA0003139882820000123
其中nx是用户输入的窗口大小，或者nx＝N/2，即有限差分阶数的一半，I0是零阶修正的贝塞尔函数，b＝4.31，是形状控制参数。类似地，可以得到y方向上的网格点j和
Figure BDA0003139882820000124
以及z方向上的k和
Figure BDA0003139882820000125
由此，组合点震源的有效数量在三维中为nxnynz，每个有效的网格点震源项都被视为一个独立的点震源，其权重为
Figure BDA0003139882820000126
出于精度考虑，2n(x,y,z)通常取有限差分算子的阶数N。因此，对于16阶有限差分，如果取nx,y,z＝N/2＝8，组合的网格点震源项的数量达到512个。对于总数为Ns个的震源，独立震源的有效数量为nxnynzNs。
b)分配这些nxnynzNs个有效的独立网格点震源网格位置、它们的权重以及片外GPU全局内存上的波场变量，由以下源代码表示：
//allocate the source position on grid point or no on grid points
cudaMalloc((void**)&config->srcpos_shot,config->NSRC_shot*sizeof(int4))；
//Allocate source weight when offset on grid points
cudaMalloc((void**)&config->srcSw8_shot,config->NSRC_shot*sizeof(float))；
c)使用具有以下源代码的CUDA复制API来从CPU内存复制到片外GPU全局内存资源，将点震源位置从CPU复制到GPU
cudaMemcpy(config->srcpos_shot,&config->srcpos_shot_thread[0],config->NSRC_shot*sizeof(int4),cudaMemcpyHostToDevice)；
//copy source weigh from CPU to GPU
cudaMemcpy(config->srcSw8_shot,&config->srcSw8_shot_thread[0],config->NSRC_shot*sizeof(float),cudaMemcpyHostToDevice)；
d)使用CUDA API(应用程序编程接口)将那些分配的片外全局内存映射到纹理内存，如以下源代码所示；
cudaBindTexture(0,d_srcSw8_tex,config->srcSw8,&channelDesc,
nx*ny*nz*config->NREC*sizeof(float))；
e)使用CUDA纹理内存获取API读取独立点震源网格位置，以及波场
Figure BDA0003139882820000131
和这些网格点上的权重
Figure BDA0003139882820000132
f)将具有权重的有效的独立网格点震源加载或注入在波场上的网格点(is+i，js+j，ks+k)上，如下述表达式所示，
Figure BDA0003139882820000133
源代码：
//fetch source location
int nslowsrc＝tex1Dfetch(d_srcpos_tex,isrc)-RADIUS；
int nmiddlesrc＝tex1Dfetch(d_srcpos_tex,1*nsrc+isrc)；
int nfastsrc＝tex1Dfetch(d_srcpos_tex,2*nsrc+isrc)；
int j＝tex1Dfetch(d_srcpos_tex,3*nsrc+isrc)；
size_t globalAddr＝(nslowsrc+RADIUS)*stridemf+nmiddlesrc*stridef+nfastsrc；
//fetch weighted source
float tmp＝tex1Dfetch(src,it*ntsrc+j)；
//fetch wavefield
float pp＝tex1Dfetch(d_p_tex,globalAddr)；
//inject the wavefield
atomicAdd(&p[globalAddr],pp)；
3.检索接收器波场
在GPU上运行的有限差分法中在任意位置处进行快速波场检索的步骤序列通常包括：
a)确定针对每个接收器的组合网格点的有效数量，并计算贡献权重。当接收器
Figure BDA0003139882820000134
不在有限差分网格点上时，必须找到一组网格点波场项(i,j,k)，它们组合起来可以获得位于
Figure BDA0003139882820000135
处的波场，可使用如下表达式：
Figure BDA0003139882820000136
这里，nx(y,z)是本领域普通技术人员通常认为的窗口大小，2nx(y,z)等于有限差分阶数N，表示为：
Figure BDA0003139882820000141
其为网格外接收器的最近网格点，
Figure BDA0003139882820000142
和
Figure BDA0003139882820000143
是权重函数
Figure BDA0003139882820000144
针对
Figure BDA0003139882820000145
类似地，计算
Figure BDA0003139882820000146
和
Figure BDA0003139882820000147
对于总数为Nr个的接收器来说，用于获取接收器站点上的波场的网格点总数为nxnynzNr。
b)在片外全局内存上分配这些nxnynzNr个网格点位置、权重以及这些网格位置上的波场变量；
c)将那些所分配的片外全局内存的映射到纹理内存；
d)使用CUDA纹理内存提取API来读取波场
Figure BDA0003139882820000148
和这些网格点上的权重
Figure BDA0003139882820000149
e)汇总由
Figure BDA00031398828200001410
加权的波场
Figure BDA00031398828200001411
根据上述表达式得到时间步长n下的位置(xr,yr,zr)上的R(xr,yr,zr,n)的检索波场。
f)使用CUDA复制API来从片外GPU全局内存资源复制到CPU内存。
4.附图的详细说明
图1是可以使用本发明的优选实施例的勘探区域101上方的地球的一部分的横截面图。重要的是应当注意，图1的勘探区域101是基于陆地的区域102。本领域的普通技术人员将认识到，地震勘探区域101产生局部地质的详细图像，以确定可能的碳氢化合物(石油和天然气)储层的位置和大小，从而确定潜在的井位103。在这些勘探区域中，声波在爆炸期间从地下岩层在不同的入射点104处反弹，反射回地表的波由地震数据记录传感器105捕获，由数据传输系统305从所述传感器105处无线传输303，然后存储到中央处理单元上的内存资源中以供后续处理，并由中央处理单元和图形处理单元进行处理，它们通常通过计算机系统设备来控制，计算机系统设备具有输出设备，例如显示器、键盘、鼠标、无线和有线网卡、打印机等。
参照图1，勘探区域上方的地球的一部分的横截面图示出了不同类型的地层102、103、104，它们将包括本发明中的地震勘探数据。特别是，本领域的普通技术人员很快就会意识，本示例示出了共中点型道集，其中地震数据轨迹按表面几何形状分类，以模拟地球中的单个反射点。在这个例子中，来自几个炮点和接收器的数据可以组合成一个图像道集，或者根据要执行的分析的类型而单独地使用。尽管本示例可以说明平面反射器和相应的图像道集类别，但是也可以使用本领域中已知的其他类型或类别的图像道集，其选择可以取决于各种地球条件或事件的存在。
此外，如图1所示，来自多个入射点或震源104的地震能量将从不同地层之间的界面反射。然后，这些反射将被多个地震数据记录接收传感器105和井103捕获，每个传感器布置在彼此的不同位置偏移110处。因为所有入射点104和所有地震数据记录传感器105都放置在不同的偏移110处，勘探地震数据或轨迹(在本领域中也称为道集)将以各种入射角108来记录。入射点104在地球中生成向下的传输射线105，其通过记录传感器105由向上传输的反射来捕获。在该示例中，井位103用附接到井筒109上的现有钻井示出，使用本领域已知的技术可获得沿着井筒109的多次测量。该井筒109用于获得震源信息(如小波)以及测井数据，包括P波速度、S波速度、密度，以及其他地震数据。图1未示出的其他传感器布置在勘探区域101内，以便还捕获解释者和本领域的普通技术人员执行各种地球物理分析所需的层位数据信息。在本示例中，将根据现场记录来对道集进行分类，以检查振幅、信噪比、动校正、频率成分、相位和其他地震属性相对于入射角108、偏移测量110、方位角以及对数据处理和成像来说重要的且本领域的普通技术人员已知的其他几何属性的依赖性。
尽管入射点或炮点104在图1中表示为常见的共中点炮点几何形状，且炮点线大多是水平的，然而本领域的普通技术人员很快就会意识到，所述模式可以很容易地以其他方式表示，例如垂直的、对角线的或这三者的组合，这又会产生不同的炮点几何形状。然而，由于操作条件的原因，如图1所示的用于小波、输入参数模型和地震数据的均匀采集的接收传感器105的均匀覆盖通常无法在整个调查区域内实现。
参见图2，其示出了可以使用本发明的优选实施例的地震勘探区域网格201。重要的是应当注意，图2的网格是基于陆地的区域102，包括任意位置处的线条状的炮点104和也处于任意位置的接收器105的一个完整的勘探计划可能因目标、预算、资源和时间等勘探特征而变化。
本领域的普通技术人员将认识到，例如网格201可产生局部地质的详细图像，以确定可能的碳氢化合物(石油和天然气)储层的位置和大小，从而确定潜在的井位103。由图2示出的陆地采集网格几何形状通常通过条带式炮击来进行，其中接收器电缆以平行线(共线方向)布置，而炮点位于垂直方向(横线方向)。在这些网格中，声波在爆炸期间从地下岩层在不同的入射点或炮点104处反弹，反射回地表的波由地震数据记录传感器105捕获，由数据传输系统305从所述传感器105处无线传输，然后存储以供在具有内存资源的中央处理单元上后续处理，然后由具有图形处理单元、片外全局内存资源和纹理内存资源的数字高性能计算系统来进行分析。虽然炮点104在图2中表示为具有几乎水平延伸的炮点线206的交叉排列模式的几何图形，然而本领域的普通技术人员很快就会意识到，所述模式可以很容易地以其他方式表示，例如垂直的、对角线的或这三者的组合。类似地，记录传感器105布置在接收器线207上，其显示为穿过炮点线206，但也可以表示为具有不同的图案(例如，对角线)。线条状炮击方法产生大范围的震源-接收器方位角，这在图形处理单元的分析期间可能是一个问题。震源-接收器方位角是参考线(例如接收器线或倾角线)与通过震源及接收器站的线之间的角度。然而，由于操作条件的原因，如图2所示的均匀覆盖通常无法在整个勘探区域内实现。
图3示出了根据本发明的某些实施例的计算机实现方法的流程图，该方法用于在网格上的任意位置使用有限差分模型来高速加载震源入射点。计算机实现方法301使用中央处理单元和图形处理单元的组合，这两者都具有用于存储某些所生成的数据的内存资源。特别是，计算机实现方法301始于中央处理单元通过外部源来接收需要存储的信息。这样的信息在步骤302处存储到中央处理单元的内存资源中，其通常包括震源时间小波、震源点和建模网格，所有这些都具有位于任意位置的震源网格点303。中央处理单元存储303信息到其内存资源中的方式是通过使用3D坐标中的声波速度模型。之后，中央处理单元在步骤304处从其内存资源中仅检索建模网格，并在步骤305处开始提取具有任意震源网格点位置的压力波场变量。然后，中央处理单元在步骤306处存储压力波场变量，并开始在步骤307处计算针对具有任意网格点位置的所存储的压力波场变量的有效震源网格点位置。
根据有限差分模型的算法来计算具有坐标
Figure BDA0003139882820000161
的点i的有效震源网格点位置，这些参数满足：
Figure BDA0003139882820000171
其中nx＝N/2，即有限差分阶数的一半。类似地，用于网格点j的有限差分模型的算法包括计算y方向上的
Figure BDA0003139882820000172
用于网格点k的有限差分模型的算法包括计算z方向上的
Figure BDA0003139882820000173
因此，组合点震源的有效数量在三维中为nxnynz，每个有效的点震源项都被算法视为一个独立的点震源，其具有权重
Figure BDA0003139882820000174
出于精度考虑，步骤307通常采用2n(x,y,z)作为有限差分算子的阶数N。因此，对于16阶有限差分，中央处理单元在步骤307取nx,y,z＝N/2＝8，这又会导致组合的网格点震源项的数量高达512个。然后，中央处理单元在步骤309处将压力波场变量的有效震源网格点位置存储到其内存资源中。在完成步骤309后，中央处理单元向图形处理单元发送消息钩子，以在步骤310启动将所存储的压力波场变量的有效震源网格点位置从中央处理单元上的内存资源中复制到位于图形处理单元上的片外全局内存资源，这可使用CUDA复制应用程序编程接口。然后，图形处理单元在步骤311处使用以下映射CUDA API来将其片外全局内存资源映射到其纹理内存资源：
cudaBindTexture(0,d_srcSw8_tex,config->srcSw8,&channelDesc,
nx*ny*nz*config->NREC*sizeof(float))；
一旦映射步骤311完成，图形处理单元在步骤312处生成初始时间步长值，其包括从震源时间小波、震源点和建模网格中获得的最小值，所有这些都包含以下三个不同的变量。这些变量是：a)稳定性要求；b)输入数据采样间隔；以及c)成像条件间隔。然后，图形处理单元在步骤313处将初始时间步长值存储到片外全局内存资源，并将图形用户界面发送到连接到图形处理单元的显示监视器，以便用户(通常是本领域的普通技术人员)在步骤314处输入用户定义的最大时间步长值。使用与步骤312相同的变量，用户然后将其最大时间步长值定义为以下各项中的最小值：a)0.1ms到2ms的稳定性要求；b)1ms到4ms的输入数据采样间隔；c)2ms到8ms的成像条件间隔。在输入最大时间步长值后，用户通过图形用户界面确认，发出信号以在步骤315处将最大时间步长值存储到其片外全局内存资源。图形处理单元在步骤316处检索初始时间步长值和用户定义的最大时间步长值，以便在步骤317处开始计算波传播算法，其中在检索到的初始时间步长值和用户定义的最大时间步长值之间的每个输入的时间步长值处使用压力波场变量的有效震源网格点位置。在步骤317执行的波传播计算是一个迭代过程，以解决由表达式(1)提出的算法。迭代步骤317基本上使用本发明的高效的计算机实现方法来注入在网格上生成波场的震源，然后传播它们，直到最后输入的时间步长值等于用户定义的最大时间步长值。图形处理单元将在步骤318验证最后输入的时间步长值是否等于用户定义的最大时间步长值。如果不是，则图形处理单元重复在步骤317处完成的计算，直到用完所有输入的时间步长值。同时，图形处理单元利用其并行处理能力，验证最后在步骤317处执行的计算是否等于用户定义的最大时间步长值，并完成波传播的迭代过程317-318，并在步骤319将新的震源网格点加载或注入到具有有效震源压力波场变量的新生成的网格上。在步骤319完成后，图形处理单元便在步骤320开始使用CUDA复制应用程序编程接口来将具有来自每个计算的波传播算法的有效位置处的压力波场变量的新震源网格点的新网格复制到中央处理单元中的随机存取存储器设备中。一旦图形处理单元完成了步骤321，它就向中央处理单元发送信号，以开始将具有来自每个计算的波传播算法的有效位置处的压力波场变量的新震源网格点的新网格存储到其内存资源(通常是硬盘驱动器)中。然后，在连接到中央处理单元上的监视器中显示图形用户界面，向用户指示该计算机实现方法已经完成了它的第一阶段，包括使用有限差分模型来高速加载震源入射点，最终数据已被存储，并且正在准备开始该计算机实现方法的第二阶段，其包括使用有限差分模型来高速检索接收器位置。这些数据可用于进一步的地球物理和地震成像分析。
来看图4，其示出了根据本发明的某些实施例的计算机实现方法的流程图，用于在网格上的任意位置使用有限差分模型来高速检索接收器位置。计算机实现方法401使用中央处理单元和图形处理单元的组合，这两者都具有用于存储某些所生成的数据的内存资源。特别是，计算机实现方法401始于中央处理单元通过外部源来接收需要存储的信息。这样的信息在步骤402处存储到中央处理单元的内存资源中，其通常包括震源时间小波、接收器点(有时仅是一个点)和建模网格，所有这些都具有位于任意位置的接收器网格点403。中央处理单元存储403信息到其内存资源中的方式是通过使用3D坐标中的声波速度模型。之后，中央处理单元在步骤404处从其内存资源中仅检索建模网格，并在步骤405处开始提取具有任意接收器网格点位置的压力波场变量。然后，中央处理单元在步骤406处存储压力波场变量，并开始在步骤407处计算针对具有任意网格点位置的所存储的压力波场变量的有效接收器网格点位置。
根据有限差分模型的算法来计算具有坐标
Figure BDA0003139882820000191
的网格点i的有效接收器网格点位置，这些参数满足：
Figure BDA0003139882820000192
其中nx＝N/2，即有限差分阶数的一半。类似地，用于网格点j的有限差分模型的算法包括计算y方向上的
Figure BDA0003139882820000193
用于网格点k的有限差分模型的算法包括计算z方向上的
Figure BDA0003139882820000194
因此，组合点震源的有效数量在三维中为nxnynz，每个有效的点震源项都被算法视为一个独立的点震源，其具有权重
Figure BDA0003139882820000195
出于精度考虑，步骤407通常采用2n(x,y,z)作为有限差分算子的阶数N。因此，对于16阶有限差分，中央处理单元在步骤407取nx,y,z＝N/2＝8，这又会导致组合的网格点接收器项的数量高达512个。然后，中央处理单元在步骤409处将压力波场变量的有效接收器网格点位置存储到其内存资源中。在完成步骤409后，中央处理单元向图形处理单元发送消息钩子，以在步骤410启动将所存储的压力波场变量的有效接收器网格点位置从中央处理单元上的内存资源中复制到位于图形处理单元上的片外全局内存资源，这可使用CUDA复制应用程序编程接口。然后，图形处理单元在步骤411处使用以下映射CUDA API来将其片外全局内存资源映射到其纹理内存资源：
cudaBindTexture(0,&d_srcSw8_tex,config->srcSw8,&channelDesc,
nx*ny*nz*config->NREC*sizeof(float))；
一旦映射步骤411完成，图形处理单元在步骤412处生成初始时间步长值，其包括从接收器点和建模网格中获得的最小值，所有这些都包含以下三个不同的变量。这些变量是：a)稳定性要求；b)输入数据采样间隔；以及c)成像条件间隔。然后，图形处理单元在步骤413处将初始时间步长值存储到片外全局内存资源，并将图形用户界面发送到连接到图形处理单元的显示监视器，以便用户(通常是本领域的普通技术人员)在步骤414处输入用户定义的最大时间步长值。使用与步骤412相同的变量，用户然后将其最大时间步长值定义为以下各项中的最小值：a)0.1ms到2ms的稳定性要求；b)1ms到4ms的输入数据采样间隔；c)2ms到8ms的成像条件间隔。在输入最大时间步长值后，用户通过图形用户界面确认，发出信号以在步骤415处将最大时间步长值存储到其片外全局内存资源。图形处理单元在步骤416处检索初始时间步长值和用户定义的最大时间步长值，以便在步骤417处开始计算波传播算法，其中在检索到的初始时间步长值和用户定义的最大时间步长值之间的每个输入的时间步长值处使用压力波场变量的有效接收器网格点位置。在步骤417执行的波传播计算是一个迭代过程，以解决由表达式(1)提出的算法。迭代步骤417基本上使用本发明的高效的计算机实现方法来检索在网格上生成的接收器波场，然后传播它们，直到最后输入的时间步长值等于用户定义的最大时间步长值。图形处理单元将在步骤418验证最后输入的时间步长值是否等于用户定义的最大时间步长值。如果不是，则图形处理单元重复在步骤417处完成的计算，直到用完所有输入的时间步长值。同时，图形处理单元利用其并行处理能力，验证最后在步骤417处执行的计算是否等于用户定义的最大时间步长值，并完成波传播的迭代过程417-418，并在步骤419将新的接收器网格点检索到具有有效接收器压力波场变量的新生成的网格上。在步骤419完成后，图形处理单元便在步骤420开始使用CUDA复制应用程序编程接口来将具有来自每个计算的波传播算法的有效位置处的压力波场变量的新接收器网格点的新网格复制到中央处理单元中的随机存取存储器设备中。一旦图形处理单元完成了步骤421，它就向中央处理单元发送信号，以开始将具有来自每个计算的波传播算法的有效位置处的压力波场变量的新接收器网格点的新网格存储到其内存资源(通常是硬盘驱动器)中。然后，在连接到中央处理单元上的监视器中显示图形用户界面，向用户指示该计算机实现方法已经完成，最终数据已被存储，并且可用于进一步的地球物理和地震成像分析。
参考图5，其示出了根据本发明的某些实施例的计算机实现方法的图形用户界面，显示了在不同时间步长间隔501处的压力波场，以及其各自的震源点和接收器点位置。特别是，图形用户界面显示了建模网格502，其具有震源炮点或入射点503、接收器504，以及由接收器504沿地球地下层所捕获的压力波场505。当计算机实现方法执行从其初始时间步长值506到其最大用户定义时间步长值510的波传播步骤时，图形处理单元将网格外的震源(在图5中由爆炸符号表示)加载或注入到网格中，并且将网格外的接收器网格点(在图5中由倒三角形表示)检索到网格中，得到了清晰和精确的波传播，如附图标记507、509和511所示。
根据本发明的优选实施例，仅作为示例性实施例而详细描述了某些硬件和软件，它们并不限制所公开的实施例的实现结构。例如，尽管已经描述了图3中的接收系统设备的许多内部和外部部件，但是本领域的普通技术人员将会理解，这样的部件及其互连是众所周知的。另外，所公开的发明的某些方面可以体现在使用一个或多个接收系统、计算机系统设备或非暂时性计算机可读存储设备来执行的软件中。可以将技术的程序方面视为通常以某种类型的机器可读介质来承载或体现的可执行代码和/或关联数据的形式的“产品”或“制品”。有形的非暂时性“存储”类型的介质和设备包括用于计算机、进程等的任何或所有的内存或其他存储器，或其相关的模块，例如各种半导体存储器、磁带驱动器、磁盘驱动器、光盘或磁盘，以及可以随时为软件编程提供存储的部件。
如本文所使用的术语“勘探区域”是指地质感兴趣的区域或体积，并且可以在任何测量尺度下与该区域或体积的几何形状、姿态和布置相关联。一个区域可能具有其中已发生的诸如折叠、断层、冷却、卸载和/或断裂的特征。
如本文所使用的术语“计算”包括各种各样的动作，包括计算、确定、处理、推导、勘探、查找(例如，在表、数据库或另一数据结构中查找)、确信等。它还可以包括接收(例如接收信息)、访问(例如访问内存中的数据)等。而且，“计算”可以包括解析、选择、选定、构建等。
如本文所用的术语“地表下”和“地下”是指在任何海拔高度上或在海拔范围内的任何一块陆地的顶面之下(无论是高于、低于还是处于海平面)，和/或低于任何地面的地表(无论是高于、低于还是处于海平面)。
除非另有明确说明，否则诸如“定义”、“创建”、“包括”、“代表”、“预分析”、“预定义”、“选择”、“构建”、“分配”、“创建”、“引入”、“消除”、“重新网格化”、“整合”、“发现”、“执行”、“预测”、“确定”、“输入”、“输出”、“识别”、“分析”、“使用”、“分发”、“干扰”、“增加”、“调整”、“合并”、“模拟”、“减少”、“分发”、“指定”、“提取”、“显示”、“执行”、“实现”和“管理”等术语可以指检索系统或其他电子设备的动作和过程，该系统或其他电子设备将某些电子设备的存储器(例如内存资源或非临时性计算机可读存储器)中的表示为物理量(电子，磁或光)的数据转换为存储器或传输设备或显示设备中的类似地表示为物理量的其他数据。
本文公开的实施例还涉及计算机实现的系统，用作用于执行本文所述的操作的检索系统的一部分。该系统可以被特别地构造用于所需目的，或者它可以包括由存储在内存资源或非暂时性计算机可读存储器中的计算机程序或代码选择性地激活或重新配置的通用计算机。这样，计算机程序或代码可以被存储或编码在计算机可读介质中，或在某种类型的传输介质上实现。计算机可读介质包括用于以诸如计算机之类的机器可读的形式存储或传输信息的任何介质或机构(“机器”和“计算机”在本文中可同义地使用)。作为一个非限制性示例，计算机可读介质可包括计算机可读存储介质(例如只读存储器ROM、随机存取存储器RAM、磁盘存储介质、光学存储介质、闪存设备等)。传输介质可以是双绞线、同轴电缆、光纤或其他一些合适的有线或无线传输介质，用于传输信号，例如电、光、声或其他形式的传播信号(例如载波、红外信号、数字信号等)。
本文所使用的接收系统或传感器105通常至少包括能够执行机器可读指令的硬件，以及用于执行产生期望结果的动作(通常是机器可读指令)的软件。另外，检索系统可以包括硬件和软件的混合体，以及计算机子系统。
硬件通常至少包括具有处理器功能的平台，例如客户端计算机(也称为服务器)和手持处理设备(例如智能电话、个人数字助理PDA或个人计算设备PCD)。此外，硬件可包括可以存储机器可读指令的任何物理设备，例如内存或其他数据存储器。其他形式的硬件包括硬件子系统，例如包括诸如调制解调器、调制解调器卡、端口和端口卡之类的传输设备。
软件包括存储在任何存储介质(例如RAM或ROM)中的任何机器代码，以及存储在其他设备(例如非暂时性计算机可读介质，例如外部硬盘驱动器或闪存)上的机器代码。软件可以包括源代码或目标代码，包含能够在客户端计算机、服务器计算机、远程桌面或终端中执行的任何指令集。
软件和硬件的组合也可以用于为所公开的发明的某些实施例提供增强的功能和性能。一个示例是直接将软件功能制造到硅芯片中。因此，应当理解，硬件和软件的组合也包括在检索系统的定义内，因此本发明将其设想为可能的等效结构和等效方法。
计算机可读介质或内存资源包括被动数据存储器，例如随机存取存储器(RAM)，以及半永久数据存储器，例如外部硬盘驱动器和外部数据库。另外，本发明的实施例可以体现在计算机的RAM中，以将标准计算机转换为新的特定计算机器。
数据结构是可以实现本发明的实施例的限定的数据组织。例如，数据结构可以提供数据的组织，或可执行代码的组织。数据信号可以跨非暂时性传输介质承载，并且跨各种数据结构存储和传输，因此可以用于传输本发明的实施例。
可以将系统计算机设计为可在任何特定的体系结构上工作。例如，该系统可以在高性能计算系统上执行，该高性能计算系统通常包括物理上连接的或通过局域网、客户端-服务器网络、广域网、互联网和其它便携式和无线设备及网络来连接的多个单台计算机的集合。
“输出设备”包括导致产生的直接行为，以及促进产生的任何间接行为。间接行为包括向用户提供软件，维护使用户能够通过其影响显示的网站，超链接至此类网站，或与执行此类直接或间接行为的实体合作。因此，用户可以单独操作或与第三方供应商合作，以便在显示设备上生成参考信号。显示设备可以作为输出设备包括在内，并且应适合于显示所需的信息，例如但不限于CRT监视器、LCD监视器、等离子设备、平板设备或打印机。显示设备可以包括已经通过使用旨在用于评估、校正和/或改善显示结果的任何常规软件进行校准的设备(例如，已经使用监视器校准软件进行了调整的彩色监视器)。作为在显示设备上显示参考图像的附加或替代，根据本发明的方法可以包括向对象提供参考图像。“提供参考图像”可以包括通过物理、电话或电子传递的方式将参考图像创建或分发给对象，通过网络提供对参考图像的访问，或者向被配置为在对象工作站上运行的对象或包含参考图像的计算机创建或分发软件。在一个示例中，提供参考图像可以涉及使对象能够经由打印机获得硬拷贝形式的参考图像。例如，信息、软件和/或指令可以被传输(例如，通过数据存储器或硬拷贝以电子或物理方式)和/或以其他方式可用(例如，通过网络)，以便于主体使用打印机来打印硬拷贝形式的参考图像。在这样的示例中，打印机可以是已经通过使用旨在用于评估、校正和/或改善打印结果的任何常规软件进行了校准的打印机(例如，已经使用颜色校正软件来调节的彩色打印机)。
一个数据库或多个数据库可以包含任何标准或专有数据库软件，例如Oracle，Microsoft Access，SyBase或DBase II。数据库可以具有字段、记录、数据和其他数据库元素，它们可以通过数据库专用软件进行关联。另外，可以映射数据。映射是将一个数据条目与另一个数据条目相关联的过程。例如，可以将包含在字符文件位置中的数据映射到第二个表中的字段。数据库的物理位置没有限制，并且数据库可以是分布式的。例如，数据库可能位于服务器的远程位置，并在单独的平台上运行。此外，可以在局域网和因特网的无线网络上访问数据库。
此外，模块、特征、属性、方法和其他方面可以被实现为软件、硬件、固件或其任何组合。在本发明的部件被实现为软件时，该部件都可以被实现为独立程序，作为较大程序的一部分，多个独立程序，静态或动态链接库，内核可加载模块，设备驱动程序，和/或计算机编程领域的技术人员现在或将来已知的所有其他方式。另外，本发明不限于在任何特定的操作系统或环境中的实现。
下面定义本文所使用的各种术语。对于在下面未对定义的在权利要求中使用的术语来说，应给予相关领域的技术人员所给出的反映在至少一个印刷出版物或授权专利中的最广泛可能的定义。
如本文所使用的置于第一实体和第二实体之间的“和/或”是指(1)第一实体、(2)第二实体以及(3)第一实体和第二实体中的一个。用“和/或”列出的多个元素应以相同的方式解释，即，如此连接的元素中的“一个或多个”。
另外，附图中的流程图和框图示出了根据本发明的各种实施例的系统、方法和计算机程序产品的可能实现的架构、功能和操作。还应注意，在一些替代实施方式中，方框中指出的功能可以不按图中指出的顺序发生。例如，取决于所涉及的功能，实际上可以基本上同时执行连续示出的两个方框，或者有时可以相反的顺序执行这些方框。还应注意，框图和/或流程图中的每个方框以及框图和/或流程图中的方框的组合可以由执行指定的硬件功能或动作的基于专用硬件的系统或者专用硬件和计算机指令的组合来实现。
尽管在前述说明书中已经针对本发明的某些优选实施例描述了本发明，并且出于说明的目的已经阐述了许多细节，但是本发明不应不当地限于已经出于说明性目的而阐述的前述内容。相反，在不背离如以下权利要求书所限定的本发明的真实范围的情况下，多种修改和替代实施例对于本领域技术人员将是显而易见的。另外应当理解，在本文的任何一个实施例中示出或描述的结构特征或方法步骤也可以在其他实施例中使用。

Claims (18)

1.一种计算机实现方法，用于在网格上的任意位置处使用有限差分模型来高速加载震源入射点，所述方法包括：
使用3D坐标中的声波速度模型来将震源时间小波、震源点以及具有任意位置处的震源网格点的建模网格存储到中央处理单元上的内存资源中；
通过所述中央处理单元从所述内存资源中检索所述建模网格；
通过所述中央处理单元从所检索的建模网格中提取具有任意震源网格点位置的压力波场变量；
将所提取的具有任意震源网格点位置的压力波场变量存储到所述中央处理单元上的内存资源中；
通过所述中央处理单元来计算具有任意网格点位置的所存储的压力波场变量的有效震源网格点位置；
将所计算出的压力波场变量的有效震源网格点位置存储到所述中央处理单元上的内存资源中；
使用CUDA复制应用程序编程接口来将所存储的压力波场变量的有效震源网格点位置从所述中央处理单元上的内存资源复制到位于图形处理单元上的片外全局内存资源中；
将所述图形处理单元的片外全局内存资源映射到所述图形处理单元的纹理内存资源中；
通过所述图形处理单元来生成初始时间步长值；
将所生成的初始时间步长值存储到所述图形处理单元的片外全局内存资源中；
将用户定义的最大时间步长值输入到所述图形处理单元中；
将所述输入的用户定义的最大时间步长值存储到所述图形处理单元的片外全局内存资源中；
从所述图形处理单元的片外全局内存资源中检索所述初始时间步长值和所述用户定义的最大时间步长值；
使用在所检索到的初始时间步长值和所述用户定义的最大时间步长值之间的每个输入的时间步长值处的压力波场变量的所复制的有效震源网格点位置，通过所述图形处理单元来计算波传播算法；
重复通过所述图形处理单元来计算波传播算法的步骤，直到最后输入的时间步长值等于所述用户定义的最大时间步长值；
通过所述图形处理单元来加载具有来自每个所计算的波传播算法的有效位置处的压力波场变量的新震源网格点的新网格；
使用CUDA复制应用程序编程接口来将具有来自每个计算的波传播算法的有效位置处的压力波场变量的新震源网格点的新网格从所述图形处理单元复制到所述中央处理单元中；和
将所复制的具有来自每个计算的波传播算法的有效位置处的压力波场变量的新震源网格点的新网格存储到所述中央处理单元的内存资源中。
2.根据权利要求1所述的计算机实现方法，其特征在于，使用3D坐标中的声波速度模型来将震源时间小波、震源点以及具有任意位置处的震源网格点的建模网格存储到中央处理单元上的内存资源中的步骤还包括由表达式s(t)在时域中表示的震源小波。
3.根据权利要求1所述的计算机实现方法，其特征在于，使用3D坐标中的声波速度模型来将震源时间小波、震源点以及具有任意位置处的震源网格点的建模网格存储到中央处理单元上的内存资源中的步骤还包括由下式表示的各个震源点：
Figure FDA0003139882810000021
4.根据权利要求1所述的计算机实现方法，其特征在于，使用3D坐标中的声波速度模型来将震源时间小波、震源点以及具有任意位置处的震源网格点的建模网格存储到中央处理单元上的内存资源中的步骤还包括由表达式VP表示的声波速度模型。
5.根据权利要求1所述的计算机实现方法，其特征在于，通过所述中央处理单元从所检索的建模网格中提取具有任意震源网格点位置的压力波场变量的步骤还包括表达式：
Figure FDA0003139882810000022
6.根据权利要求1所述的计算机实现方法，其特征在于，通过所述中央处理单元来计算具有任意网格点位置的所存储的压力波场变量的有效震源网格点位置的步骤还包括针对每个网格点外震源的多达512个网格点位置。
7.根据权利要求1所述的计算机实现方法，其特征在于，通过所述中央处理单元来计算具有任意网格点位置的所存储的压力波场变量的有效震源网格点位置的步骤还包括为每个网格点位置分配权重，由下式表示：
Figure FDA0003139882810000031
8.根据权利要求1所述的计算机实现方法，其特征在于，将用户定义的最大时间步长值输入到所述图形处理单元中的步骤包括输入以下中的最小值：0.1ms到2ms的稳定性要求；1ms到4ms的输入数据采样间隔；以及2ms到8ms的成像条件间隔。
9.根据权利要求1所述的计算机实现方法，其特征在于，使用在所检索到的初始时间步长值和所述用户定义的最大时间步长值之间的每个输入的时间步长值处的压力波场变量的所复制的有效震源网格点位置，通过所述图形处理单元来计算波传播算法的步骤还包括表达式：
Figure FDA0003139882810000032
10.根据权利要求1所述的计算机实现方法，其特征在于，通过所述图形处理单元来加载具有来自每个所计算的波传播算法的有效位置处的压力波场变量的新震源网格点的新网格的步骤还包括表达式：
Figure FDA0003139882810000033
11.一种计算机实现方法，用于在网格上的任意位置处使用有限差分模型来高速检索接收器波场位置，所述方法包括：
使用3D坐标中的声波速度模型来将震源时间小波、接收器位置以及具有任意位置处的接收器网格点的建模网格存储到中央处理单元上的内存资源中；
通过所述中央处理单元从所述内存资源中检索所述建模网格；
通过所述中央处理单元从所检索的建模网格中提取具有任意接收器网格点位置的接收器压力波场变量；
将所提取的具有任意接收器网格点位置的接收器压力波场变量存储到所述中央处理单元上的内存资源中；
通过所述中央处理单元来计算具有任意网格点位置的所存储的接收器压力波场变量的有效接收器网格点位置；
将所计算出的接收器压力波场变量的有效接收器网格点位置存储到所述中央处理单元上的内存资源中；
使用CUDA复制应用程序编程接口来将所存储的接收器压力波场变量的有效接收器网格点位置从所述中央处理单元上的内存资源复制到位于图形处理单元上的片外全局内存资源中；
将所述图形处理单元的片外全局内存资源映射到所述图形处理单元的纹理内存资源中；
通过所述图形处理单元来生成初始时间步长值；
将所生成的初始时间步长值存储到所述图形处理单元的片外全局内存资源中；
将用户定义的最大时间步长值输入到所述图形处理单元中；
将所述输入的用户定义的最大时间步长值存储到所述图形处理单元的片外全局内存资源中；
从所述图形处理单元的片外全局内存资源中检索所述初始时间步长值和所述用户定义的最大时间步长值；
使用在所检索到的初始时间步长值和所述用户定义的最大时间步长值之间的每个输入的时间步长值处的接收器压力波场变量的所复制的有效接收器网格点位置，通过所述图形处理单元来计算波传播算法；
重复通过所述图形处理单元来计算波传播算法的步骤，直到最后输入的时间步长值等于所述用户定义的最大时间步长值；
通过所述图形处理单元来加载具有来自每个所计算的波传播算法的有效位置处的接收器压力波场变量的新接收器网格点的新网格；
使用CUDA复制应用程序编程接口来将具有来自每个计算的波传播算法的有效位置处的接收器压力波场变量的新接收器网格点的新网格从所述图形处理单元复制到所述中央处理单元中；和
将所复制的具有来自每个计算的波传播算法的有效位置处的接收器压力波场变量的新接收器网格点的新网格存储到所述中央处理单元的内存资源中。
12.根据权利要求11所述的计算机实现方法，其特征在于，使用3D坐标中的声波速度模型来将接收器位置以及具有任意位置处的接收器网格点的建模网格存储到中央处理单元上的内存资源中的步骤还包括由下式表示的各个接收器点：
Figure FDA0003139882810000051
13.根据权利要求11所述的计算机实现方法，其特征在于，使用3D坐标中的声波速度模型来将接收器位置以及具有任意位置处的接收器网格点的建模网格存储到中央处理单元上的内存资源中的步骤还包括由表达式VP表示的声波速度模型。
14.根据权利要求11所述的计算机实现方法，其特征在于，通过所述中央处理单元从所检索的建模网格中提取具有任意接收器网格点位置的接收器压力波场变量的步骤还包括表达式：
Figure FDA0003139882810000052
15.根据权利要求11所述的计算机实现方法，其特征在于，通过所述中央处理单元来计算具有任意网格点位置的所存储的接收器压力波场变量的有效接收器网格点位置的步骤还包括针对每个网格点外接收器的多达512个网格点位置。
16.根据权利要求11所述的计算机实现方法，其特征在于，将用户定义的最大时间步长值输入到所述图形处理单元中的步骤包括输入以下中的最小值：0.1ms到2ms的稳定性要求；1ms到4ms的输入数据采样间隔；以及2ms到8ms的成像条件间隔。
17.根据权利要求11所述的计算机实现方法，其特征在于，使用在所检索到的初始时间步长值和所述用户定义的最大时间步长值之间的每个输入的时间步长值处的接收器压力波场变量的所复制的有效接收器网格点位置，通过所述图形处理单元来计算波传播算法的步骤还包括表达式：
Figure FDA0003139882810000053
18.根据权利要求11所述的计算机实现方法，其特征在于，通过所述图形处理单元来加载具有来自每个所计算的波传播算法的有效位置处的接收器压力波场变量的新接收器网格点的新网格的步骤还包括表达式：
Figure FDA0003139882810000054
CN202110734477.2A 2020-06-30 2021-06-30 使用有限差分模型进行高速多源加载和波场检索的方法 Active CN113945994B (zh)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
US16/917,545 2020-06-30
US16/917,545 US11281825B2 (en) 2020-06-30 2020-06-30 Computer-implemented method for high speed multi-source loading and retrieval of wavefields employing finite difference models

Publications (2)

Publication Number Publication Date
CN113945994A true CN113945994A (zh) 2022-01-18
CN113945994B CN113945994B (zh) 2023-07-04

Family

ID=79030960

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110734477.2A Active CN113945994B (zh) 2020-06-30 2021-06-30 使用有限差分模型进行高速多源加载和波场检索的方法

Country Status (2)

Country Link
US (1) US11281825B2 (zh)
CN (1) CN113945994B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115081267B (zh) * 2022-05-18 2023-08-29 内蒙古农业大学 基于优化的时空变步长有限差分地震波数值模拟方法
CN116578416A (zh) * 2023-04-26 2023-08-11 中国人民解放军92942部队 一种基于gpu虚拟化的信号级仿真加速方法

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2013033651A1 (en) * 2011-08-31 2013-03-07 Spectraseis Ag Full elastic wave equation for 3d data processing on gpgpu
WO2015155597A2 (en) * 2014-04-07 2015-10-15 Cgg Services Sa Attenuating pseudo s-waves in acoustic anisotropic wave propagation
WO2016018464A1 (en) * 2014-07-30 2016-02-04 Chevron U.S.A. Inc. Re-ordered interpolation and convolution for faster staggered-grid processing
CN105467443A (zh) * 2015-12-09 2016-04-06 中国科学院地质与地球物理研究所 一种三维各向异性弹性波数值模拟方法及系统
CN106199692A (zh) * 2015-05-30 2016-12-07 中国石油化工股份有限公司 一种基于gpu的波动方程反偏移方法
US20180203144A1 (en) * 2015-05-20 2018-07-19 Optasense, Inc. Interferometric Microseismic Imaging Methods and Apparatus
US20190033479A1 (en) * 2017-07-25 2019-01-31 Advanced Geophysical Technology Inc. Memory efficient q-rtm computer method and apparatus for imaging seismic data
WO2019030659A1 (en) * 2017-08-10 2019-02-14 Seismic Apparition Gmbh METHOD FOR ACQUIRING AND PROCESSING SEISMIC DATA
CN110988988A (zh) * 2019-11-25 2020-04-10 中国矿业大学(北京) 基于垂直裂缝介质的地震波波场模拟方法及装置

Family Cites Families (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB2381314B (en) 2001-10-26 2005-05-04 Westerngeco Ltd A method of and an apparatus for processing seismic data
US8549500B2 (en) 2007-02-14 2013-10-01 The Mathworks, Inc. Saving and loading graphical processing unit (GPU) arrays providing high computational capabilities in a computing environment
US9995835B2 (en) 2013-07-17 2018-06-12 Chevron U.S.A. Inc. System and method of implementing finite difference time domain models with multiple accelerated processing components (APCS)
CN104570081B (zh) * 2013-10-29 2017-12-26 中国石油化工股份有限公司 一种积分法叠前时间偏移地震资料处理方法及系统
CN104570080A (zh) * 2013-10-29 2015-04-29 中国石油化工股份有限公司 一种海量数据叠前逆时偏移多gpu卡协同快速计算方法
CN105319581B (zh) * 2014-07-31 2018-01-16 中国石油化工股份有限公司 一种高效的时间域全波形反演方法
CN105717539B (zh) * 2016-01-28 2018-01-30 中国地质大学（北京） 一种基于多gpu计算的三维tti介质逆时偏移成像方法
US10678553B2 (en) 2017-10-10 2020-06-09 Apple Inc. Pro-active GPU hardware bootup
US11105942B2 (en) * 2018-03-27 2021-08-31 Schlumberger Technology Corporation Generative adversarial network seismic data processor

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2013033651A1 (en) * 2011-08-31 2013-03-07 Spectraseis Ag Full elastic wave equation for 3d data processing on gpgpu
WO2015155597A2 (en) * 2014-04-07 2015-10-15 Cgg Services Sa Attenuating pseudo s-waves in acoustic anisotropic wave propagation
WO2016018464A1 (en) * 2014-07-30 2016-02-04 Chevron U.S.A. Inc. Re-ordered interpolation and convolution for faster staggered-grid processing
US20180203144A1 (en) * 2015-05-20 2018-07-19 Optasense, Inc. Interferometric Microseismic Imaging Methods and Apparatus
CN106199692A (zh) * 2015-05-30 2016-12-07 中国石油化工股份有限公司 一种基于gpu的波动方程反偏移方法
CN105467443A (zh) * 2015-12-09 2016-04-06 中国科学院地质与地球物理研究所 一种三维各向异性弹性波数值模拟方法及系统
US20190033479A1 (en) * 2017-07-25 2019-01-31 Advanced Geophysical Technology Inc. Memory efficient q-rtm computer method and apparatus for imaging seismic data
WO2019030659A1 (en) * 2017-08-10 2019-02-14 Seismic Apparition Gmbh METHOD FOR ACQUIRING AND PROCESSING SEISMIC DATA
CN110988988A (zh) * 2019-11-25 2020-04-10 中国矿业大学(北京) 基于垂直裂缝介质的地震波波场模拟方法及装置

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
谭嘉言等: "基于GPU计算的地震波场高阶有限差分正演研究", 《CT理论与应用研究》 *
谭嘉言等: "基于GPU计算的地震波场高阶有限差分正演研究", 《CT理论与应用研究》, vol. 25, no. 01, 29 February 2016 (2016-02-29), pages 1 - 12 *

Also Published As

Publication number Publication date
CN113945994B (zh) 2023-07-04
US11281825B2 (en) 2022-03-22
US20210406427A1 (en) 2021-12-30

Similar Documents

Publication Publication Date Title
EP2497043B1 (en) Seismic imaging systems and methods employing a 3d reverse time migration with tilted transverse isotropy
US20210103065A1 (en) Determining properties of a subterranean formation using an acoustic wave equation with a reflectivity parameterization
US11733413B2 (en) Method and system for super resolution least-squares reverse time migration
CN113945994B (zh) 使用有限差分模型进行高速多源加载和波场检索的方法
US20210262329A1 (en) Method for Generating Initial Models For Least Squares Migration Using Deep Neural Networks
Yuan et al. Seismic source tracking with six degree‐of‐freedom ground motion observations
US11474267B2 (en) Computer-implemented method and system employing compress-sensing model for migrating seismic-over-land cross-spreads
WO2022187685A1 (en) Method and system for faster seismic imaging using machine learning
WO2013033651A1 (en) Full elastic wave equation for 3d data processing on gpgpu
US20160047925A1 (en) Method of Determining Seismic Acquisition Aperture
US20210018638A1 (en) Generating enhanced seismic velocity models using geomechanical modeling
US20190204463A1 (en) Variable Aperture Estimation using Bottom-Up Ray Tracing
US11573347B2 (en) Computing program product and method that interpolates wavelets coefficients and estimates spatial varying wavelets using the covariance interpolation method in the data space over a survey region having multiple well locations
US20240142649A1 (en) Method and system for determining migration data using multiblock gathers
US11614555B2 (en) Method and system for connecting elements to sources and receivers during spectrum element method and finite element method seismic wave modeling
WO2024087126A1 (en) Method and system for determining migration data using cross-spread gathers and parallel processing
US20220137248A1 (en) Computing program product and method for prospecting and eliminating surface-related multiples in the beam domain with deghost operator
Ogarko et al. Tomofast-x 2.0: an open-source parallel code for inversion of potential field data, to recover density, susceptibility and magnetisation vector, with topography and wavelet compression
US20240184007A1 (en) Method and system for determining attenuated seismic time
Ning Multicomponent Distributed Acoustic Sensing: Concept, Theory, and Applications
US20240125961A1 (en) Method and system for determining wavefield components using independent component analysis
US12000971B2 (en) Method and system for seismic processing using virtual trace bins based on offset attributes and azimuthal attributes
WO2024087002A1 (en) Methods and systems for determining attenuated traveltime using parallel processing
US20230184972A1 (en) Method and system for seismic processing using virtual trace bins based on offset attributes and azimuthal attributes
Ogarko et al. Tomofast-x 2.0: an open-source parallel code for inversion of potential field data with topography using wavelet compression

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant