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

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

【技术实现步骤摘要】
使用有限差分模型进行高速多源加载和波场检索的方法


[0001]本专利技术大体上涉及使用高速计算方法的地球物理勘测和地震数据处理，其中在具有带宽和存储限制的图形处理单元中采用了有限差分模型。

技术介绍

[0002]1.概述
[0003]即使使用当今的高性能计算机器，也可能没有足够的计算能力来使用单个统一计算方案对一整片勘探区域进行建模。因此，本领域的技术人员将转而利用混合处理单元来实现高级数值编程解决方案。这些混合方法需要有效且实用的实现方式，这是因为这些单元的几个不同建模内核在地震勘探模型的不同区域(或子域)中运行，甚至可能在不同的有限差分网格分辨率下运行。
[0004]两种主要的处理单元是图形处理单元(GPU)和中央处理单元(CPU)。GPU以比同类CPU更低的价格提供高性能的计算能力，并采用各种硬件模块，如台式机或中心。GPU的主要目的是提供改进的数值性能，因此GPU和CPU之间的性能有着显著的不同。例如，一个特定的GPU可以在一秒内进行一万亿次浮点运算，而CPU的性能由每个插槽的内核数量(例如4核)和所述插槽的性能(即每秒时钟周期数)决定。
[0005]然而，GPU依靠线程级并行来隐藏较长的片外内存访问延迟，因此明智地利用片上内存资源(包括寄存器文件、共享内存和数据缓存)对应用程序的性能至关重要。相比之下，GPU的一个片上特性是它们可以并行地寻址不同的片上内存，从而增加内存带宽和延迟。因此，管理GPU片上内存资源对于应用程序的开发人员来说是一项非常重要的任务。更重要的是，由于GPU的不同代的片上内存资源不同，性能可移植性成为一项艰巨的挑战。然而，在本领域中众所周知，GPU对通用编程提供的支持非常有限，这是因为在GPU中执行的代码不能自分配内存，并且需要由主机CPU来完成此操作。
[0006]使用混合方法的这些复杂性意味着从CPU到GPU的编程很困难，并且需要在相同或不同时间生成高级算法。显然，这影响了本领域的技术人员在复杂介质中计算大规模地震波动方程模型和地震波传播模拟时可以在石油和天然气工业中执行的工作量。
[0007]2.现阶段的GPU架构
[0008]使用统一计算设备架构(CUDA)编程模型的C/C++语言中使用最广泛的高性能并行应用程序之一是NVIDIA的Tesla架构。NVIDIA架构的一个例子通常是一组流式多处理器，每个处理器都有自己的处理器和共享内存(用户管理的缓存)。这些处理器完全能够执行整数和单精度浮点运算，并带有用于双精度的附加内核。然而，多处理器中的所有流处理器都与未被硬件缓存的全局设备内存共享获取、控制和内存单元。
[0009]在这些架构中，CUDA将线程排列成线程块，可以在分配给所述线程块的任何共享内存位置处读写。因此，线程块内的线程可以通过共享内存进行通信或使用共享内存作为用户管理的缓存，这是因为共享内存的延迟比全局内存的延迟低两个数量级。
[0010]然而，由于在大多数架构中前向和后向阶段访问和修改相同的数据，因此混合
CPU/GPU实现需要在主机和GPU内存之间进行连续的内存传输，以保持数据的一致性。由于频繁和大量的内存传输，当前的GPU性能受到极大的影响，然后在GPU内核中执行步骤，以访问驻留在图形双倍数据速率(GDDRAM)模块中的数据。这提供了一个额外的好处，因为GDDR模块可以在同一时钟周期上请求和接收数据。因此，主机只需要计算执行环境、内核调用，以及在必要时往来于磁盘的I/O传输。
[0011]3.声波传播
[0012]声波的计算建模一直是一个活跃的研究领域，并且已经从航空航天等领域发展到地震学、地球物理学甚至气象学。事实上，空间阵列中声波的计算建模和模拟是许多科学和工程应用的基础，在石油和天然气壳层和储层勘探中也是一样。然而，迄今为止，在大型或复杂结构中基于声波传播模型来获得良好的结果仍然是一个主要的计算挑战。
[0013]该领域中最常用的模拟声波传播的方法通常包括两个阶段的过程：(1)计算代表声学空间的脉冲响应；以及(2)将脉冲响应与无回声记录的或合成生成的源信号卷积。然而，由于这些脉冲响应依赖于对时变空间声场建模的精确计算，因此它们的主要挑战仍然是精确求解方法(例如有限差分)的计算要求。在大多数情况下，用于求解声学传播的方法要求来自空间离散化的每个波长大约3
‑
5个节点(即网格点)。假设最小频率范围为3Hz，最大频率为100Hz，节点间距为5
‑
50m，结果是10立方公里的声学空间充满了约8百万到80亿个节点。因此，声波方程的数值求解器仅限于相当简单、受限或小的空间，这通常被认为是耗时的。为了克服这个问题，大多数当前的声学仿真系统使用几何声学技术，它们仅对较高频率和早期反射来说是准确的，并且可能无法对衍射效应进行建模。可以看到，克服了一个问题会产生另一个问题。
[0014]4.全波场反演(FWI)
[0015]FWI是一种偏微分方程约束优化方法，它迭代地最小化测量波场和计算波场之间的失配范数。地震FWI涉及多次迭代，单次迭代可能涉及以下计算：(1)求解前向方程；(2)求解伴随方程；以及(3)对这些前向方程和伴随方程的解进行卷积，以产生成本函数的梯度。应当注意，对于二阶优化方法，例如Gauss
‑
Newton方法来说，还需要进行(4)求解扰动前向方程。
[0016]经典时域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 and waveform 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 in frequency
‑
space seismic waveform inversion,Geophysical Journal International,133,341
‑
362)。最近，Sh本文档来自技高网...

【技术保护点】

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

【专利技术属性】
技术研发人员：张昌华，
申请(专利权)人：中国石油化工股份有限公司，
类型：发明
国别省市：