CN109490948A - 地震声学波动方程矢量并行计算方法 - Google Patents

Abstract

本发明属于地球探测与信息技术领域，尤其涉及一种地震声学波动方程矢量并行计算方法，包括：(一)地震模拟；用计算机和数值方法模拟地震波场在地下介质中的传播过程，包括以下分步骤：(二)计算并行加速；对地震模拟项目做并行化需求分析，然后针对串行代码中用时最多的，且可以用矢量运算代替的部分用AVX指令集并行代码代替。本发明提供的地震声学波动方程矢量并行计算方法，不仅可以模拟出高密度采集情形的地震响应特征和变化规律，而且能大幅节省计算指令，最终实现约3～7倍的二次加速比的能力。

Description

地震声学波动方程矢量并行计算方法

技术领域

本发明属于地球探测与信息技术领域，尤其涉及一种地震声学波动方程矢量并行计算方法。

背景技术

波动方程正演模拟是地球物理研究的基础工作，它服务于地震资料采集、处理、解释的各个阶段，其重要性不言而喻。但在通常情况下，大规模地震正演模拟需要消耗大量的计算资源，如较大的计算时间花费和海量的内存消耗。

高性能计算是推动地震勘探技术进步的重要引擎。近年来，并行计算技术的发展使得大规模地震正演模拟的快速、高精度实现具有了现实可能性。一方面，适用于分布式存储系统的消息传递MPI并行以及适用于共享式存储系统的OpenMP并行发展迅速，已被引入地震波场模拟以实现加速。另一方面，基于CPU模板缓冲区优化技术的SIMD指令集不断发展并逐步趋向成熟，通过将数据矢量化的方式提高CPU计算带宽，从而实现在一个CPU指令周期内对多个数据的同时操作，即实现矢量运算，因此能节约计算指令并有效减轻一级指令缓存压力。时间域地震正演所具有的两个特征，即时间维度迭代的顺序性和空间维度的无序性，采用传统的任务级并行不可避免地会产生时间同步、模拟计算任务切割、周围局部波场信息在节点间的传递等额外花费从而影响加速效果。本专利基于多核CPU平台，利用常规OpenMP拆分任务，然后利用AVX指令完成数据并行，从而实现二次加速以提高正演效率。

发明内容

针对现有技术存在的问题，本发明提供了一种基于AVX指令集的地震声学波动方程矢量并行计算方法，它能隐式给出完整的地震波场。该方法对传播介质的空间变化没有限制，并且在使用足够细的网格时其计算结果精度很高。此外，这种方法还能处理介质的不同流变性质，以及生成地震波场快照，从而对地震数据的解释来说是一种重要的辅助工具。本发明的另一个优点是在计算上比其他方法更加高效。

本发明提供一种地震声学波动方程矢量并行计算方法，所述地震声学波动方程矢量并行计算方法包括：

(一)地震模拟；

(二)计算并行加速。

其中(一)地震模拟，包括以下分步骤；

(1)提供已知的地质构造信息；

可以根据测井资料用平均法获得地震级的速度和密度场，以及通过地震相干体的生成、油藏单元的定义以及研究区域断层数据的融合进行地质模拟。使用联合测井(声波和密度测井)信息可以进一步提高模拟精度。

(2)通过地震正演预测一组检波器在空间特定位置所能记录到的振动。

在第(2)步中，为了保证模拟计算的正常运行，地震正演参数(网格间距、网格样点数、边界处理、时间步长)的选择及模型试验配置应按以下分步执行：

(a)根据最大震源频率和最小速度，确定网格间距：

其中dx表示网格间距，c_min表示最小速度，f_max表示最大频率。实际网格间距的确定还依赖于模拟方法的选择。对于时间、空间域的地震模拟有限差分方法，在每一波长范围内至少需要5个网格点。

(b)根据地质模型大小确定网格点数。

(c)在模型两侧、顶部和底部边界为吸收带(用于衰减人工反射波场)分别增加一定的网格点。海绵边界吸收方法需要宽度达到4个波长，这里的波长λ＝2c_max/f_d，其中f_d是地震信号的主频。

(d)根据以下稳定条件和精度标准，选择时间步长。

其中dt表示时间采样间隔，c_max表示最大地下介质速度，dx_min表示最小横向采样间距。

(e)定义炮点-检波点排列。若是用于地表采集试验，则将炮点、检波点按照等间隔的方式布设在地表。

针对声学介质情形，各方向上的法向应力均等于声压的负值，介质属性参数只需要外加密度，即

τ_xz→0,τ_zz＝τ_xx→-p

μ→0,λ→k

其中τ_xx、τ_zz为法向应力；τ_xz为剪切应力；λ、μ为弹性常数；k为体积模量。对于速度场v和声压场p，运动的向量方程是

其中，r是位置向量；ρ(r)表示介质密度；f(r,t)为外源力；p(r,t)表示声压场。

为了实现地震模拟，求解上述波动方程采用有限差分法的交错网格形式，计算公式如下：

其中，Δt是时间采样间隔；h表示空间采样间距；为时间n、空间(i,k)处的声压场；表示空间位置(i+1/2,k)处的介质浮力。k_i,k为空间(i,k)处的体积模量；为时间n+1/2、空间(i+1/2,k)处的水平速度分量。

其中(二)计算并行加速，包括：首先要对地震模拟项目做并行化需求分析，然后针对串行代码中用时最多的，且可以用矢量运算代替的部分用AVX指令集并行代码代替。AVX代码实现通常按照以下流程：

第一步：包含头AVX文件#include<immintrin.h>；

第二步：定义__m256矢量类型的常量、变量或指向数组的指针，注意数组内存地址需要32位对齐。若起始地址不是32的倍数，调试时可能报错误segfault，并且内存不对齐时不利于CPU访问，造成性能下降。

第三步：用_mm256_load或_mm256_setzero等函数给__m256矢量类型变量赋初值，并将*__m256指针首地址指向32位对齐数组的首地址。

第四步：矢量运算，使用_mm256_add、_mm256_sub等函数对AVX数据类型做矢量运算。

值得注意的是，AVX指令函数操作的对象为256bit矢量类型数据。__m256矢量类型矢量可存储16个整型数据、8个单精度浮点数，或者4个双精度浮点数。

本发明提供的地震声学波动方程矢量并行计算方法，不仅可以模拟出高密度采集情形的地震响应特征和变化规律，而且能大幅节省计算指令，最终实现约3～7倍的二次加速比的能力。

附图说明

通过下面结合附图进行的详细描述，本发明的上述和其它目的、特点和优点将会变得更加清楚，其中：

图1是实施例二维网格声学场以及介质属性的位置；

图2为常规OpenMP并行计算示意图；

图3是实施例的地震声学波动方程矢量并行计算方法的波场传播过程计算流程图；

图4是实施例提供的声学波动方程有限差分法的时间推进计算过程；

图5是实施例提供的AVX矢量化加速比曲线。

图6是实施例提供的三维盐丘标准模型；

图7是实施例提供的波场快照切片，t＝1.6s；

图8是实施例提供的计算效率对比图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

下面结合附图及具体实施例对本发明的应用原理作进一步描述。

本发明的地震声学波动方程矢量并行计算方法包括：用计算机和数值方法模拟地震波场在地下介质中的传播过程。首先提供已知的地质构造信息，然后通过地震正演预测一组检波器在空间特定位置所能记录到的振动。

利用CPU模板缓冲区优化技术实现内存/计算密集的大规模地震声学波动方程正演模拟矢量并行计算方法，网格配置如图1所示。图1显示的是二维网格声学场以及介质属性的位置，其垂向表示时间轴。左栏显示压力场p以及介质压缩模量k的位置。考虑半时间步长时，右栏显示的是x方向(圆形)、z方向(正方形)的速度分量和介质浮力b。

如图2所示，其实现过程仅将工作任务分解到不同CPU上，而在单一CPU中每次指令仅对一个数据进行操作。在OpenMP方案基础上，使用SIMD技术将模板缓冲区数据矢量化，进一步挖掘CPU的运算潜力。如图3所示，其实现过程不仅将工作任务分解到不同CPU上，而且在单一CPU中每次指令能对多个数据同时进行操作。在地震模拟实现过程中，采用256位AVX指令集对波场传播模拟采用8路矢量并行，将增加数倍CPU运算带宽，从而大幅节省计算指令，最终得到约3～7倍的二次加速比。

本发明的实施例的地震声学波动方程矢量并行计算方法首先使用时刻t_-1和t处的波场值，计算时间t₊₁的整个空间范围(x_min到x_ma)。图中，更新后时间t₊₁的点用上数第二行的实线圆圈显示，而有待更新的数据点用上数第四行的虚线圆圈显示；沿空间轴采用串行顺序的算法更新数据。从图4可以看出，一旦时间t₊₁最左边的值可用，将会得到时间t₊₂部分空间域的所有信息。传统算法没有使用这种策略，在更新时间t₊₂任一空间位置的波场以前，先确定时间t₊₁整个空间范围内的波场。

图5是本发明实施例提供的AVX矢量化加速比曲线。当Z方向网格固定；点线为沿X方向矢量化，虚线代表沿Y方向矢量化；模型网格尺寸是：Nx＝Ny＝676，Nz＝201。当Z方向网格固定；蓝色曲线为沿X方向矢量化，黄色代表沿Y方向矢量化。指令带宽的数倍提升会对寄存器数据的加载速度有更高的要求，为了提升寄存器加载数据速度，需要注意三点：一是CPU的访问速度由快到慢依次是暂存器、缓存、内存，存储量由小到大依次也是暂存器、缓存、内存；二是n维数组在内存中只有一个方向的地址连续，而其余n-1个维度方向的地址不连续，而AVX函数加载非连续地址数据到YMM矢量暂存器会比地址连续数据更加耗时；三是当数据量超过缓存大小时，访问大数据中的一个元素通常比访问小数据中的一个元素要慢。在计算x和y方向偏导数时，每一个场分量在z方向是对齐的，只需要将其邻近的场值一次性赋值给若干个__m256类型数组即可实现快速推进。当x方向网格点数大于3倍y方向网格点数时，沿y方向矢量化加速比相对沿x方向大约提高1倍。

图6是本发明实施例提供的三维盐丘标准模型，当Z方向网格固定；点线为沿X方向矢量化，虚线代表沿Y方向矢量化；模型网格尺寸是：Nx＝Ny＝676，Nz＝201。

图7是本发明实施例提供的波场快照切片，t＝1.6s。正演参数：PML的吸收层数是32，时间采样间隔1ms，空间采样间隔△x＝△y＝△z＝10m，震源坐标(338，338，3)，震源主频：30Hz。模拟精度为时间二阶差分、空间8阶差分。

图8是本发明实施例提供的计算效率对比图。使用常规的OpenMP并行和基于OpenMP+AVX的非常规方案进行正演模拟。试验硬件参数：Intel(R)Xeon(R)CPU E5-26650@2.40GHz(8核*2CPU)，RAM 16GB，L1级缓存：8*32KB指令缓存，8*32KB数据缓存；L2级缓存：8*256KB；三级缓存：20MB。经过测试，采用OpenMP最大可达约14倍加速，在OpenMP基础上利用AVX优化可达约3.8倍的二次加速。

根据本发明实施例，模拟精度的选择由本领域技术人员自行设定，如根据实际地质数据采集的频率范围、计算硬件环境灵活设定。在本发明所述方法的计算步骤中，所述依赖模型的速度取值范围、模型尺寸、地震记录长度与计算设备的内存环境有关，而计算效率与参与计算的CPU核数及AVX的优化程度保持一致。

本发明提供一种基于AVX指令集的地震声学波动方程矢量并行计算方法，它能处理介质的不同流变性质，以及生成地震波场快照，从而对地震数据的解释来说是一种重要的辅助工具。

根据本发明的地震声学波动方程矢量并行计算方法，能隐式给出完整的地震波场。该方法对传播介质的空间变化没有限制，并且在使用足够细的网格时其计算结果精度很高。本发明的另一个优点是在计算上比其他方法更加高效。

本发明的地震声学波动方程矢量并行计算方法，具有可靠地模拟地震响应特征和规律的优越性，不需要增加任何硬件配置便能将正演效率提高约3-7倍。此外，计算指令的减少节约了CPU指令缓存资源。

尽管已经参照其示例性实施例具体显示和描述了本发明，但是本领域的技术人员应该理解，在不脱离权利要求所限定的本发明的精神和范围的情况下，可以对其进行形式和细节上的各种改变。

Claims (3)

1.地震声学波动方程矢量并行计算方法，其特征在于，包括以下两个步骤：

(一)地震模拟；

用计算机和数值方法模拟地震波场在地下介质中的传播过程，包括以下分步骤：

(1)提供已知的地质构造信息；

根据测井资料用平均法获得地震级的速度和密度场，以及通过地震相干体的生成、油藏单元的定义以及研究区域断层数据的融合进行地质模拟；

(2)通过地震正演预测一组检波器在空间特定位置所能记录到的振动；

(二)计算并行加速；

对地震模拟项目做并行化需求分析，然后针对串行代码中用时最多的，且可以用矢量运算代替的部分用AVX指令集并行代码代替。

2.根据权利要求1所述的地震声学波动方程矢量并行计算方法，其特征在于，在第(一)步的第(2)分步骤中，为了保证模拟计算的正常运行，地震正演参数的选择及模型试验配置应按以下分步骤执行：

(a)根据最大震源频率和最小速度，确定网格间距：

其中dx表示网格间距，c_min表示最小速度，f_max表示最大频率。实际网格间距的确定还依赖于模拟方法的选择。对于时间、空间域的地震模拟有限差分方法，在每一波长范围内至少需要5个网格点；

(b)根据地质模型大小确定网格点数；

(c)在模型两侧、顶部和底部边界为吸收带分别增加一定的网格点；海绵边界吸收方法需要宽度达到4个波长，这里的波长λ＝2c_max/f_d，其中f_d是地震信号的主频；

(d)根据以下稳定条件和精度标准，选择时间步长：

其中dt表示时间采样间隔，c_max表示最大地下介质速度，dx_min表示最小横向采样间距；

(e)定义炮点、检波点排列；若是用于地表采集试验，则将炮点、检波点按照等间隔的方式布设在地表。

3.根据权利要求1所述的地震声学波动方程矢量并行计算方法，其特征在于，步骤(二)中AVX代码实现通常按照以下流程：

第一步：包含头AVX文件#include<immintrin.h>；

第二步：定义__m256矢量类型的常量、变量或指向数组的指针，注意数组内存地址需要32位对齐。若起始地址不是32的倍数，调试时可能报错误segfault，并且内存不对齐时不利于CPU访问，造成性能下降；

第三步：用_mm256_load或_mm256_setzero等函数给__m256矢量类型变量赋初值，并将*__m256指针首地址指向32位对齐数组的首地址；

第四步：矢量运算，使用_mm256_add、_mm256_sub等函数对AVX数据类型做矢量运算。