CN107561585A - 一种多核多节点并行三维地震波场生成方法和系统 - Google Patents
一种多核多节点并行三维地震波场生成方法和系统 Download PDFInfo
- Publication number
- CN107561585A CN107561585A CN201710845834.6A CN201710845834A CN107561585A CN 107561585 A CN107561585 A CN 107561585A CN 201710845834 A CN201710845834 A CN 201710845834A CN 107561585 A CN107561585 A CN 107561585A
- Authority
- CN
- China
- Prior art keywords
- data
- wave field
- parallel
- wave
- mpi
- 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.)
- Pending
Links
Landscapes
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明公布了一种多核多节点并行三维地震波场生成方法和系统，基于高阶有限差分和多核多节点并行加速，实现高精度高效率地求解三维地震波动方程；包括：获取数据信息；确定地震波传播三维正演模拟模型的初始条件、生成边界条件和算法稳定性条件，对地震波传播方程进行高阶有限差分和数值模拟；结合MPI对多节点并行及OpenMP对多核计算的处理能力，设计基于MPI/OpenMP的混合多核多节点并行的正演模拟算法结构；由MPI控制多炮数据循环并行计算，并输出多炮波场切片结果和地面地震记录数据，模拟生成勘探地震波场。本发明在大大提高计算效率的同时，有效地提高计算精度，满足三维地震波场生成对计算效率和精度的要求。
Description
技术领域
本发明属于数值模拟与高性能计算技术领域，涉及地震波场正演模拟方法，具体涉及一种基于多核多节点并行加速的三维地震波场的正演模拟方法及其系统。
背景技术
地震波场正演模拟是有效反映地下介质构造变化的一种数值模拟方法。随着油气勘探开发技术的发展，以及计算机运算能力的提升，地震波场正演模拟的方法已得到了广泛的应用。
地震波场数值正演模拟是在地下介质结构和参数已知的情况下，利用理论计算的方法研究地震波在地下介质中的传播规律，并获得人工合成地震记录的一种技术。在现有的方法中，有限差分方法具有很好的灵活性而被广泛应用于数值计算之中。
其中，有限差分方法的研究始于1968年，Alterman和Karal在文献[1]中，首先将弹性波有限差分模拟的方法应用到了拟合地层的层状介质中。之后，有限差分技术快速发展，文献[2]记载了Alford等比较了高阶差分与低阶差分保持同样精度时对网格间距的要求，深入研究有限差分法算法精度。随后，Kelly等文献[3]中研究了用有限差分法人工合成了地震记录，将方法与实际数据紧紧连系在了一起。在文献[4]中，Virieux提出了稳定有限差分交错网格的差分形式，适用于任何泊松比的介质，同时，交错网格提高了局部计算精度和算法的稳定性。之后，文献[5]记载了Levander将Virieux的结果进行了差分阶数的提升，得到更高精度的结果。在文献[6]中，Crase则将精度提高到任意阶数，进行高阶差分，因而对CPU、内存以及储存空间的要求也大大提升。在文献[7]中，Graves又在三维空间下完成了这一方法的实践应用，使该方法更好的结合到我们实际的勘探工作中。董国良等在文献[8]～[9]中对一阶弹性波交错网格高阶差分的方法以及其稳定性问题进行了深入研究。文献[10]中记载了张剑锋和刘铁林将网格的研究进一步发展，提出了矩形网格和三角网格相结合的方法，提高有限差分法在构造较复杂地区结果的精度。在文献[11]中，董国良和李培明又在前人的基础上针对频散这一问题进行了更加详细的讨论。文献[12]记载了李胜军等对地震波数值模拟中的频散压制方法的总结分析。之后，有限差分方法越来越多的应用到了勘探工业界，在实际工作中起着重要作用。目前，有限差分法可以说是最常用的数值模拟方法，已经较为成熟，它正向着高精度和高效率的波场模拟方向发展。
正演模拟的方法还包括限元法、伪谱法、克希霍夫积分法等等。但从方法上而言，在有限差分法中，我们用相应的空间、时间的差分代替波场函数的空间、时间导数，对于各种条件下的介质都有着广泛的应用，能够真实并准确地模拟地震波在地下介质中的运动方式。同时，其具有运算快、数据量小的优势，这是其他方法都难以比拟的。因此，该方法是勘探应用最重要的方法之一。
上述的地震波场数值模拟技术是在地下介质结构和参数已知的情况下，利用理论计算的方法研究地震波在地下介质中的传播规律，并获得人工合成地震记录的一种技术。
近年来，多核资料处理技术得到了迅速的发展，以往的处理方式多是串行执行的，而现今，并行模式已经得到广泛的利用，MPI(Message Passing Interface，信息传递接口)和OpenMP(Open Multi-Processing，共享存储并行编程)技术目前占据很主导的地位。
随着版本不断更新，MPI和OpenMP的计算内核也在不断改进，目前的计算语言已大大简化，同时，操作的方式也更加人性化，特别是OpenMP并行语言已经集成到了多种语言上，尤其是C++语言。对于这两种语言，MPI更加注重于节点之间的并行实现，OpenMP则是针对多核技术的并行语言。结合两者的优势，取长补短，能更好地提升计算的效率。
目前，在这一技术领域的应用中，三维弹性波场正演模拟结果能更加接近野外采集的实际数据，具有很好的参照性和对比性，同时，也能为之后的三维反演提供模型数据基础。然而，三维数据也带来了庞大的计算量，保证计算精度的同时，计算效率的改善也是必须考虑的问题。
参考文献：
[1]Alterman,Z.,and Karal,F.C..Propagation of elastic wave in layeredmedia by finite difference method[J].Bulletin of the Seismological Society ofAmerica,1968,58(1):367-398.
[2]Alford,R.M.,Kelly,K.R.,Boore,D.M..Accuracy of finite differencemodeling of the acoustic wave equation[J].Geophysics.1974,39(6):834-842.
[3]Kelly,K.R.,et al.Synthetic seismograms；a finite differenceapproach[J].Geophysics.1976,41(1):2-27.
[4]Virieux,J..P-SV wave propagation in heterogeneous media:Velocity-stress finite-difference method(Shear waves)[J].Geophysics,1986,51(4):889-901.
[5]Levander,A.R..Fourth-order finite-difference P-SV seismograms[J].Geophysics.1988.53(11):1425-1436.
[6]Crase,E..High-order(space and time)finite-difference modeling ofelastic wave equation[J].Expanded Abstracts of 60th SEG AnnualMeeting.1990.987-991.
[7]Graves,R.W..Simulating Seismic Wave Propagation in 3D ElasticMedia Using Staggered-Grid Finite-Differences[J].Bulletin of seismologicalSociety of America.1996,86(4):1091-1106.
[8]董国良,马在田,曹景忠等.一阶弹性波方程交错网格高阶差分解法[J].地球物理学报.2000a,43(3):411-419.
[9]董国良,马在田,曹景忠等.一阶弹性波方程交错网格高阶差分解法稳定性研究[J].地球物理学报.2000b,43(6):856-864.
[10]张剑锋,刘铁林.各向异性介质中弹性波的数值模拟[J].固体力学学报.2000,21(3):234-242.
[11]董国良,李培明.地震波传播数值模拟中的频散问题[J].天然气工业.2004,24(6):53-56.
[12]李胜军,孙成禹,高建虎等.地震波数值模拟中的频散压制方法分析[J].石油物探.2008,47(5):444-449.
发明内容
为了解决计算效率的问题，本发明提供一种多核多节点并行三维地震波场生成方法和系统，基于高阶有限差分和多核多节点并行加速，实现高精度、高效率地求解三维地震波动方程，大大提高计算效率的同时，有效地提高计算精度，满足对计算效率和精度的要求。
本发明的核心是：本发明将MPI/OpenMP多核多节点混合并行模式用于多炮集三维的勘探地震波正演模拟生成。基于多核多节点并行生成勘探地震波，首先获取地震子波及观测系统数据、获取地质参数信息；根据地震子波数据和地质参数信息，再根据初始条件、边界条件以及算法稳定性条件，对地震波传播方程进行三维高阶有限差分，通过数值模拟得到高精度的勘探地震波传播数据；基于MPI/OpenMP多核多节点混合并行模式改进算法结构，完成三维并行加速处理，实现高精度结果的高效输出；最后输出波场时间切片结果和地面地震记录结果。因此，本发明从高阶差分和多核多节点并行加速入手，是三维地震正演模拟和MPI/OpenMP多核多节点混合并行加速技术的有效结合，能够数十倍加快传统三维正演模拟的计算速度，大大提高计算效率的同时，保证有效的计算精度。
本发明提供的技术方案如下：
多核多节点并行三维地震波场生成方法，对地震波传播方程进行高阶有限差分，通过数值模拟计算得到高精度的勘探地震波传播数据。在计算过程中，通过多核多节点并行加速处理，高精度、高效地生成地面地震记录数据和波场切片结果：包括如下步骤：
1)获取地震子波及观测系统数据，获取地质参数信息；
2)确定地震波传播三维正演模拟模型的初始条件、生成边界条件和算法稳定性条件，对地震波传播方程进行高阶有限差分，通过数值模拟得到高精度的勘探地震波传播数据；
3)设计基于MPI/OpenMP的多核多节点并行的正演模拟算法结构，通过多节点设置、数据一维化并行、数据计算顺序分配三种并行加速处理方法，实现高精度结果的高效输出；
4)MPI控制多炮数据循环并行计算，依次输出多炮波场切片结果和地面地震记录数据，模拟生成勘探地震波场。
本发明还提供一种多核多节点并行三维地震波场生成系统，包括震源子波获取模块、观测系统模块、地质参数获取模块、非地质参数调整模块、边界条件生成模块、多核多节点正演模块和地震波场生成模块；具体地：
A.震源子波获取模块，用于获取地震子波数据；
B.观测系统模块，用于建立地面数据收集系统；
C.地质参数获取模块，用于获取模型所对应对的地质参数信息；
C01.勘探数据采集子模块，采集地震数据；
C02.处理地震剖面子模块，利用地震工业处理软件处理所述地震数据得出地震剖面和建模数据；
C03.分析合成模型子模块，对所述地震剖面和建模数据进行分析得出对应的地质参数信息，合成模型数据。
D.非地质参数调整模块，用于调整合适的参数，满足稳定性条件；
E.边界条件生成模块，用于生成模型所对应的完全匹配层边界条件；
F.多核多节点正演模块，用于完成弹性波波动方程高阶差分离散和MPI/OpenMP多核多节点混合并行加速处理，达到高精度、高效率；
F01.多核高阶离散子模块，利用高阶差分和OpenMP对波动方程进行一维化数值离散分流并行计算；
F02.MPI多炮并行子模块，利用MPI分配多节点并行处理多炮数据。
G.地震波场生成模块，用于生成地面地震记录和波场切片结果。
与现有技术相比，本发明的有益效果是：
本发明提供一种多核多节点并行三维地震波场生成方法和系统，基于高阶差分和MPI/OpenMP多核多节点混合并行加速，实现高精度、高效率地求解三维地震波动方程，大大提高计算效率的同时，保证最有效的精度要求。本发明在具体实施时，应用了MPI/OpenMP多核多节点混合并行加速方法，在保证稳定性的前提下，大大提高计算的效率。
本发明中，高阶有限差分方法具有较高的精度和较好的数值稳定性，特别是在三维计算中，微小数值变化所带来的影响更为严重，需要高阶高精度的计算，才能得到更加准确，逼近真实数据的结果。同时，针对高阶有限差分的大数据量，基于MPI/OpenMP多核多节点的并行加速计算能大大提高其计算效率。另外，本发明实施例基于MPI/OpenMP编程语言为基础，优化了模拟过程中的并行构架，提高计算效率并保证计算稳定，高效且高精度地完成地震波场正演模拟。再而，完全匹配边界条件可以用于波场边界上进行振幅衰减。本发明针对高阶差分对边界层进行调整，保证了边界吸收的效果。
由此，本发明提供的多核多节点并行三维地震波场生成方法能够为三维复杂区地震波传播规律分析研究提供高精度、高效率、稳定性好的数值模拟计算。
附图说明
图1是本发明提供的多核多节点并行三维地震波场生成方法的流程框图。
图2是本发明提供的多核多节点并行三维地震波场生成系统的结构框图。
图3是本发明提供的多核多节点并行加速正演模拟方法的流程框图。
图4是本发明提供的完全匹配层边界条件的示意图；
其中，密集横竖线标示的阴影部分为完全匹配层边界区域。
图5是本发明实施例采用多核多节点并行三维地震波场模拟的实验中模型图、震源及观测系统分布图：
其中，(a)为实验模型纵波速度参数分布；(b)为实验模型横波速度参数分布；(c)为震源及观测系统分布，星号表示震源，延X方向序号为1至10，虚线表示三维地震记录接收线，在地震记录标准格式SEGY(SEG：Society of Exploration Geophysicists，勘探地球物理学家协会；Y是一种储存格式的标识；SEGY是一种由勘探地球物理学家协会提出的标准磁带数据格式之一，是石油勘探行业地震数据的最为普遍的格式之一)记录中延Y、X方向依次排列。
图6是本发明实施例采用多核多节点并行三维地震波场模拟的实验中计算结果地震记录展示图：
其中，(a)为第一炮的地震记录，以标准地震记录格式SEGY显示，纵轴为时间，以秒单位，横轴是记录排列号；(b)为第六炮的地震记录，以标准地震记录格式SEGY显示，纵轴为时间，以秒单位，横轴是记录排列号；(c)为第十炮的地震记录，以标准地震记录格式SEGY显示，纵轴为时间，以秒单位，横轴是记录排列号。
图7是本发明实施例采用现有的单核单节点与采用本发明的多核多节点并行数值模拟三维地震波场计算时间的结果对比。
具体实施方式
下面结合附图，通过实施例进一步描述本发明，但不以任何方式限制本发明的范围。
本发明提供一种基于多核多节点并行三维地震波场生成方法和系统，基于高阶有限差分和MPI/OpenMP多核多节点混合并行加速，对地震波传播方程进行高阶有限差分和多核多节点并行加速数值模拟，得到高精度的勘探地震波传播数据，从而高精度、高效地生成地面地震记录数据和波场切片结果，由此实现高精度、高效率地求解三维地震波动方程，大大提高计算效率的同时，保证最有效的精度要求的目的。
本发明中，高阶有限差分方法具有较高的精度和较好的数值稳定性，特别是在三维计算中，微小数值变化所带来的影响更为严重，需要高阶高精度的计算，才能得到更加准确，逼近真实数据的结果。同时，针对高阶有限差分的大数据量，基于MPI/OpenMP多核多节点的并行加速计算能大大提高其计算效率。另外，本发明实施例基于MPI/OpenMP编程语言为基础，优化了模拟过程中的并行构架，提高计算效率并保证计算稳定，高效且高精度地完成地震波场正演模拟。再而，完全匹配边界条件可以用于波场边界上进行振幅衰减。本发明针对高阶差分对边界层进行调整，保证了边界吸收的效果。如图4所示，对三维数据而言，在其坐标方向的切片上均满足图示的吸收边界层设计。
本发明采用的三维地震正演模拟的计算方程为地震弹性波应力-速度方程，如公式(1)：
式中，ρ为密度；Vi为速度分量(i：x、y、z)；σij为应力分量(i、j：x、y、z)；λ、μ为描述弹性介质的拉梅常数；x、y、z为笛卡尔坐标系的三分量；t为时间。
时间和空间差分方式定义为2阶时间、2L阶空间差分；i、j、r分别代表计算网格的坐标位置；n表示时间迭代中时间网格的位置。差分格式如下：
式中，Δx、Δy、Δz为笛卡尔坐标系中三个方向上的差分间距；Δt为时间差分间距。差分定义为2阶时间、2L阶空间差分。为差分系数。
本发明在具体实施时，应用了MPI/OpenMP多核多节点并行加速方法，在保证稳定性的前提下，大大提高计算的效率。图1所示为本发明提供的多核多节点并行三维地震波场生成方法的流程，包括以下步骤：
1)获取地震子波及观测系统数据，获取地质参数信息；
2)根据所述的地质参数信息，确定正演模拟模型的初始条件、生成边界条件以及算法稳定性条件，对地震波传播方程进行高阶有限差分，数值模拟得到高精度的勘探地震波传播数据；
三维地震正演模拟的计算方程为地震弹性波应力-速度方程(公式(1))。利用高阶有限差分对波动方程进行数值离散计算，具体通过步骤3)中MPI/OpenMP多核多节点并行加速处理数值模拟得到高精度的勘探地震波传播数据，生成地面地震记录数据和波场切片结果；通过公式(2)-(11)进行方程计算和差分，得到差分计算的结果；
其中，生成边界条件是生成模型所对应的完全匹配层边界条件，通过边界波场衰减达到吸收边界波场的效果，能很大程度上减少人工反射的干扰，边界条件的生成主要包括如下步骤：首先设置衰减系数得到衰减因子，再分解运算方程并加入衰减项，最后代入原方程重新进行差分计算。完全匹配层是基于衰减因子在边界处的衰减作用实现的，本发明采用衰减因子如下：
其中，Vmax的取值为速度模型中最大纵波速度，δ为匹配层宽度，x、y、z为三个方向上的网格点位置，Δx、Δy、Δz为笛卡尔坐标系中三个方向上的差分间距，R为理想边界层反射系数，即衰减系数(一般取值介于10-4-10-10)；为调节系数，可以自行调节边界吸收的强度；ddx(x)、ddy(y)、ddz(z)为衰减因子，当ddx(x)、ddy(y)、ddz(z)不为零时是衰减状态，ddx(x)、ddy(y)、ddz(z)均为零时是不衰减状态。加入完全匹配层边界条件的具体过程为：首先设置衰减系数和调节系数，进行测试和对比，得到合适的衰减因子；再将方程(1)分解在x、y、z三个方向上，只要有相关方向的差分计算，均需要在该方向做分解；而后在x、y、z三个方向上分别加入变量衰减项，衰减项中包含衰减因子；最后带入差分方程在x、y、z三个方向上分别进行数值差分计算。
稳定性是弹性波数值模拟中需关注的重要标准之一。数值算法一般可分为无条件稳定算法和有条件稳定算法。本发明采用的递推算法属于有条件稳定算法，对于本发明计算公式(2)-(10)所用的交错网格而言，其差分格式要满足以下式15的稳定性条件：
式中，Δx、Δy、Δz为笛卡尔坐标系中三个方向上的差分间距；Δt为时间差分间距。差分定义为2阶时间、2L阶空间差分。此处研究的是各项同性弹性介质，地震波传播速度与传播方向无关。Vmax的取值为速度模型中最大纵波速度。ε是稳定性因子，由差分的方式决定，同时也可以自行调节稳定性的强度。模型需要完成了稳定性的检测，才能保证计算结果的稳定。
3)基于MPI/OpenMP多核多节点并行加速模式，修改正演模拟算法结构，完成并行加速处理，实现高精度结果的高效输出；
根据三维地震正演模拟的弹性波应力-速度方程，本发明在OpenMP分流并行计算中，将变量数据计算顺序分为速度分量计算和应力分量计算两部分，先后进行速度分量和应力分量的计算，在时间循环中按顺序进行迭代，使得高阶有限差分算法适用于OpenMP并行加速；
图3是本发明具体实施中多核多节点并行加速正演模拟的流程框图，如图3所示，本发明中，地震波动方程MPI/OpenMP混合并行计算生成包括如下步骤：
31)多节点设置：根据炮域信息对总炮数进行平均分配，并按照MPI节点划分方式，即通过MPI_Comm_size和MPI_Comm_rank函数将数据对数据炮集数进行划分，将不同炮域的信息划分到不同的节点上进行计算；
32)数据一维化并行：为了提高并行计算的效率，在计算中需要将三维数据体一维化，即将三个维度的数据以一维数组的形式顺序排列，并将三维数组对应的差分计算方程转换为相应一维数组对应的差分计算方程，再应用OpenMP进行数据分流并行计算；
33)数据计算顺序分配：对整个波场空间的模拟而言，每个网格点均需要大量的迭代运算，为了避免计算数据互相干扰而产生紊乱的现象，需要通过并行断点进行计算顺序分配，在OpenMP分流并行中，将数据计算分为速度分量计算和应力分量计算两部分，先后进行速度分量和应力分量的计算，即在循环中以先后顺序进行速度分量和应力分量迭代；
34)计算结果输出及数据返回：一维数据计算完成后，对数据进行重构，按照32)中三维数据转一维数据的方法，反向进行逆排列，重新将数据的维度以三维SEGY数据的形式进行输出，三维SEGY数据是三维实际地震数据的标准格式。
4)MPI控制多炮数据循环并行计算，依次输出多炮波场切片结果和地面地震记录数据，模拟生成勘探地震波场。
本发明中，高阶有限差分方法具有较高的精度和较好的数值稳定性。在三维计算中，微小数值变化所带来的影响非常严重，需要高阶高精度的计算，才能得到更加准确，逼近真实数据的结果。同时，针对高阶有限差分的大数据量，基于多核多节点的并行加速计算能大大提高其计算效率。另外，本发明以MPI/OpenMP混合并行计算为基础，优化了模拟过程中的并行构架，高效且高精准地完成弹性波场正演模拟计算。再而，完全匹配边界条件可以用于波场边界上进行振幅衰减。本发明针对高阶差分对边界层进行调整，保证了边界吸收的效果。图4是本发明完全匹配层边界条件的示意图。如图4所示，对三维数据而言，在其坐标方向的切片上均满足图示的完全匹配吸收边界层设计。
为了让本发明的目的、技术和优势更加清晰，下面将结合实施方式和附图，对本发明做进一步详细说明。在此，本发明的示意性实施方式及其说明用于解释本发明，但并不作为对本发明的限定。
图2是本发明提供的地震波场模拟系统的结构框图，包括震源子波获取模块、观测系统模块、地质参数获取模块、非地质参数调整模块、边界条件生成模块、多核多节点正演模块和地震波场生成模块；具体地：
A.震源子波获取模块，用于获取地震子波数据，地震子波通常选取雷克子波，以震源的形式导入波动方程正演模拟计算中；
B.观测系统模块，用于建立地面数据收集系统，其包括地震炮点信息、地震检波点信息以及地震炮点和检波点的关系信息；
C.地质参数获取模块，用于获取模型所对应对的地质参数信息，包括速度、密度等重要的模型参数；
C01.勘探数据采集子模块，采集地震数据；
C02.处理地震剖面子模块，利用地震工业处理软件处理所述地震数据得出地震剖面和建模数据；
C03.分析合成模型子模块，对所述地震剖面和建模数据进行分析得出对应的地质参数信息，合成模型数据。
D.非地质参数调整模块，用于调整合适的参数，满足稳定性条件，参数包括地震子波的主频、时间空间差分网格的大小；
E.边界条件生成模块，用于生成模型所对应的完全匹配层边界条件，通过边界波场衰减达到吸收边界波场的效果，能很大程度上减少人工反射的干扰，边界条件的生成主要包括步骤如下：首先设置衰减系数，再分解运算方程并加入衰减项，最后代入原方程重新进行差分计算；
F.多核多节点正演模块，用于完成弹性波波动方程高阶差分离散和MPI/OpenMP多核多节点混合并行加速处理，达到高精度、高效率；
F01.多核高阶离散子模块，利用高阶差分和OpenMP对波动方程进行一维化数值离散分流并行计算，首先通过高阶差分离散方程，再通过一维化方法转换数据和方程，最后通过OpenMP进行数据分流并行计算；
F02.MPI多炮并行子模块，利用MPI分配多节点并行处理多炮数，对炮数进行划分，通过MPI给不同炮分配节点，循环计算，直到所有炮计算完成。
G.地震波场生成模块，用于生成地面地震记录和波场切片结果，该模块的两个结果计算如下：波场切片的结果是在每一个时间节点时进行保存，即保存地震波场传播过程中的状态；地面地震记录的结果是通过计算每一时间节点地面的波场情况，在最后时间单节点循环计算结束后进行保存。
下面通过实例对本发明的效果做进一步说明。
具体实施采用三维盐丘速度模型，利用本发明提供方法生成10炮三维地震波。图5是本发明实施例采用多核多节点并行三维地震波场模拟的实验中模型图、震源及观测系统分布图：其中，(a)为实验模型纵波速度参数分布；(b)为实验模型横波速度参数分布；(c)为震源及观测系统分布，星号表示震源，延X方向序号为1至10，虚线表示三维地震记录接收线，在地震记录标准格式SEGY记录中延Y、X方向依次排列。图6是本发明实施例采用多核多节点并行三维地震波场模拟的实验中计算结果地震记录展示图，以其中三炮计算结果为例进行展示：其中，(a)为第一炮的地震记录，以标准地震记录格式SEGY显示，纵轴为时间，以秒单位，横轴是记录排列号；(b)为第六炮的地震记录，以标准地震记录格式SEGY显示，纵轴为时间，以秒单位，横轴是记录排列号；(c)为第十炮的地震记录，以标准地震记录格式SEGY显示，纵轴为时间，以秒单位，横轴是记录排列号。
实施例针对三维盐丘模型，采用多核多节点并行三维地震波场生成方法进行正演模拟，具体包括如下步骤：
1)首先将实验分为四组，输入同样的模型、地震子波和观测系统文件，并对模型、地震子波和观测系统数据进行同样的处理；
2)对四组实验，设置相同的初始条件、边界条件，并验证算法的稳定性；
3)第一、二组实验采用传统的单节点进行计算，炮点间采用顺序计算；第三、四组实验应用本发明中的MPI炮点并行方法，将10炮的数据分配到10个节点上进行计算；
4)第一、三实验采用传统的单核进行计算，直接针对三维数据、三维方程进行循环计算；第二、四组实验采用本发明中的OpenMP多核并行计算技术，将数据一维化后，在4个核上进行分流并行计算；
5)输出地面地震记录数据，将所有时间节点保存的观测系统接收数据(波场时间切片的地面部分)在数值计算结束后输出，并记录四组实验的计算时间。
实验结果记录了0.4s波场单核单节点、4核单节点、单核10节点、4核10节点计算时间的对比，如图7所示。结果表明，采用本发明提出的多核多节点(此处为4核10节点)并行三维地震波场生成方法，相比于本实施实验所用的原始单核单节点，能够大幅加快正演模拟的速度，加速可达约24.2倍。在实际应用中，如果使用更高级别的计算设备，可取得更佳的加速效果。采用本发明提出的多核多节点并行三维地震波场生成方法有很好的应用前景。
需要注意的是，公布实施例的目的在于帮助进一步理解本发明，但是本领域的技术人员可以理解：在不脱离本发明及所附权利要求的精神和范围内，各种替换和修改都是可能的。因此，本发明不应局限于实施例所公开的内容，本发明要求保护的范围以权利要求书界定的范围为准。
Claims (8)
1.一种多核多节点并行三维地震波场生成方法，对地震波传播方程进行高阶有限差分，通过数值模拟计算得到高精度的勘探地震波传播数据；在计算过程中，通过多核多节点并行加速处理，高精度、高效地生成地面地震记录数据和波场切片结果；包括如下步骤：
1)获取地震子波及观测系统数据，获取地质参数信息；
2)确定地震波传播三维正演模拟模型的初始条件、生成边界条件和算法稳定性条件，对地震波传播方程进行高阶有限差分，通过数值模拟得到高精度的勘探地震波传播数据；包括：
2A)采用地震弹性波应力-速度方程作为三维地震正演模拟的计算方程；
2B)利用高阶有限差分对波动方程进行数值离散计算；其中，时间和空间差分方式采用2阶时间、2L阶空间差分；
2C)生成边界条件：生成地震波传播三维正演模拟模型所对应的完全匹配层边界条件，通过边界波场衰减达到吸收边界波场的效果；生成边界条件首先设置衰减系数得到衰减因子，分解运算方程并加入衰减项，再代入原方程重新进行差分计算；基于衰减因子在边界处的衰减作用实现完全匹配层；
2D)设置数值模拟中差分方式的稳定性条件；
3)进行地震波动方程MPI/OpenMP混合并行计算：设计基于MPI/OpenMP多核多节点混合并行的正演模拟算法结构，采用并行加速处理方法；包括：
3A)多节点设置：根据炮域信息对总炮数进行平均分配，并按照MPI节点划分方式，将不同炮域的信息划分到不同的节点上进行计算；
3B)数据一维化并行：在计算中将三维数据体一维化，即将三个维度的数据以一维数组的形式顺序排列，并将三维数组对应的差分计算方程转换为相应一维数组对应的差分计算方程，再应用OpenMP进行数据分流并行计算；
3C)数据计算顺序分配：通过并行断点进行计算顺序分配，避免计算数据互相干扰而产生紊乱；具体在OpenMP分流并行计算中，将变量分为速度分量和应力分量，变量数据计算顺序为采用先后进行速度分量和应力分量的计算，在时间循环中按顺序进行迭代，使得高阶有限差分算法适用于OpenMP并行加速；
3D)计算结果输出及数据返回：一维数据计算完成后，对数据进行重构，按照3B)中三维数据转一维数据的方法，反向进行逆排列，重新将数据的维度以三维SEGY数据的形式进行输出；
4)通过MPI控制多炮数据循环并行计算，依次输出多炮波场切片结果和地面地震记录数据，模拟生成勘探地震波场。
2.如权利要求1所述多核多节点并行三维地震波场生成方法，其特征是，步骤2C)中，衰减因子表示如式12～14：
其中，Vmax的取值为速度模型中最大纵波速度，δ为匹配层宽度，x、y、z为三个方向上的网格点位置，Δx、Δy、Δz为笛卡尔坐标系中三个方向上的差分间距，R为理想边界层反射系数，即衰减系数(一般取值介于10-4-10-10)；为调节系数，可自行调节边界吸收的强度；ddx(x)、ddy(y)、ddz(z)为衰减因子，当ddx(x)、ddy(y)、ddz(z)不为零时是衰减状态，ddx(x)、ddy(y)、ddz(z)均为零时是不衰减状态；
完全匹配层边界条件的实现过程具体为：
首先设置衰减系数和调节系数，进行测试和对比，得到合适的衰减因子；再将地震弹性波应力-速度方程分解在差分计算的各方向上；而后在各方向上分别加入变量衰减项，衰减项中包含衰减因子；最后带入差分方程在各方向上分别进行数值差分计算。
3.如权利要求1所述多核多节点并行三维地震波场生成方法，其特征是，2D)中设置稳定性条件表示为式15：
<mrow>
<msub>
<mi>&Delta;tV</mi>
<mi>max</mi>
</msub>
<msqrt>
<mrow>
<mfrac>
<mn>1</mn>
<mrow>
<msup>
<mi>&Delta;x</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mo>+</mo>
<mfrac>
<mn>1</mn>
<mrow>
<msup>
<mi>&Delta;y</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mo>+</mo>
<mfrac>
<mn>1</mn>
<mrow>
<msup>
<mi>&Delta;z</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
</mrow>
</msqrt>
<mo>&le;</mo>
<mi>&epsiv;</mi>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>15</mn>
<mo>)</mo>
</mrow>
</mrow>
式中，Δx、Δy、Δz为笛卡尔坐标系中三个方向上的差分间距；Δt为时间差分间距；Vmax的取值为速度模型中最大纵波速度；ε是稳定性因子，由差分的方式决定，同时也可自行调节稳定性的强度。
4.如权利要求1所述多核多节点并行三维地震波场生成方法，其特征是，3B)中MPI节点划分根据炮域信息对总炮数进行平均分配，具体通过MPI_Comm_size函数和MPI_Comm_rank函数将数据对数据炮集数进行划分。
5.如权利要求1所述多核多节点并行三维地震波场生成方法，其特征是，步骤4)中通过MPI控制多炮数据循环并行计算，具体利用MPI分配多节点并行处理多炮数据，对炮数进行划分，通过MPI为不同炮分配节点，循环计算，直到所有炮计算完成。
6.如权利要求1所述多核多节点并行三维地震波场生成方法，其特征是，步骤4)中输出多炮波场切片结果和地面地震记录数据，具体是：在每一个时间节点保存波场切片的结果；计算每一时间节点地面的波场情况，在最后时间单节点循环计算结束后，保存地面地震记录的结果。
7.利用权利要求1～6所述多核多节点并行三维地震波场生成方法实现的多核多节点并行三维地震波场生成系统，包括震源子波获取模块、观测系统模块、地质参数获取模块、非地质参数调整模块、边界条件生成模块、多核多节点正演模块和地震波场生成模块；具体地：
A.震源子波获取模块，用于获取地震子波数据；
B.观测系统模块，用于建立地面数据收集系统；
C.地质参数获取模块，用于获取模型所对应对的地质参数信息；
D.非地质参数调整模块，用于调整合适的参数，满足稳定性条件；
E.边界条件生成模块，用于生成模型所对应的完全匹配层边界条件；
F.多核多节点正演模块，用于完成弹性波波动方程高阶差分离散和MPI/OpenMP多核多节点混合并行加速处理，达到高精度、高效率生成三维地震波场；包括多核高阶离散子模块和MPI多炮并行子模块；
多核高阶离散子模块利用高阶差分和OpenMP对波动方程进行一维化数值离散分流并行计算；MPI多炮并行子模块利用MPI分配多节点并行处理多炮数据；
G.地震波场生成模块，用于生成地面地震记录和波场切片结果。
8.如权利要求7所述多核多节点并行三维地震波场生成系统，其特征是，所述地质参数获取模块包括：勘探数据采集子模块、处理地震剖面子模块和分析合成模型子模块；勘探数据采集子模块用于采集地震数据；处理地震剖面子模块通过利用地震工业处理软件处理所述地震数据得出地震剖面和建模数据；分析合成模型子模块用于对所述地震剖面和建模数据进行分析得出对应的地质参数信息，合成模型数据。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710845834.6A CN107561585A (zh) | 2017-09-19 | 2017-09-19 | 一种多核多节点并行三维地震波场生成方法和系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710845834.6A CN107561585A (zh) | 2017-09-19 | 2017-09-19 | 一种多核多节点并行三维地震波场生成方法和系统 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN107561585A true CN107561585A (zh) | 2018-01-09 |
Family
ID=60981242
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710845834.6A Pending CN107561585A (zh) | 2017-09-19 | 2017-09-19 | 一种多核多节点并行三维地震波场生成方法和系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107561585A (zh) |
Cited By (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108445537A (zh) * | 2018-02-07 | 2018-08-24 | 中国地质大学（武汉） | 基于Spark的叠前地震数据AVO弹性参数反演方法及系统 |
| CN109342575A (zh) * | 2018-09-06 | 2019-02-15 | 北京大学 | 弹性波偏振系统 |
| CN109490948A (zh) * | 2018-11-16 | 2019-03-19 | 西南石油大学 | 地震声学波动方程矢量并行计算方法 |
| CN110927784A (zh) * | 2019-11-29 | 2020-03-27 | 中国地质大学（武汉） | 一种基于OpenMP的频率域高斯束格林函数计算方法 |
| CN111766623A (zh) * | 2020-05-25 | 2020-10-13 | 国家超级计算天津中心 | 一种基于mpi进程拓扑的地震波场三维正演模拟方法 |
| CN112578450A (zh) * | 2020-10-16 | 2021-03-30 | 中国石油大学(华东) | 基于频散介质标量波方程有限差分的正演模拟方法及装置 |
| CN112859162A (zh) * | 2019-11-28 | 2021-05-28 | 中国石油天然气集团有限公司 | 基于伪谱法的三维声波方程空间域正演方法及装置 |
| CN113536638A (zh) * | 2021-07-22 | 2021-10-22 | 北京大学 | 基于间断有限元和交错网格的高精度地震波场模拟方法 |
| CN113960663A (zh) * | 2021-10-25 | 2022-01-21 | 中国地质大学（北京） | 基于并行计算的三维各向异性衰减正演模拟方法及系统 |
| WO2024087002A1 (en) * | 2022-10-25 | 2024-05-02 | Saudi Arabian Oil Company | Methods and systems for determining attenuated traveltime using parallel processing |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2013055840A3 (en) * | 2011-10-11 | 2013-10-24 | Saudi Arabian Oil Company | High performance and grid computing with quality of service control |
| CN105572730A (zh) * | 2015-12-15 | 2016-05-11 | 中国科学院地质与地球物理研究所 | 三维复杂结构声波正演方法 |
| WO2017034988A1 (en) * | 2015-08-25 | 2017-03-02 | Saudi Arabian Oil Company | High performance and grid computing of hydrocarbon exploration and production data with partitioning quality of service control |
| CN106646593A (zh) * | 2016-08-30 | 2017-05-10 | 国家超级计算天津中心 | 一种跨节点并行的三维起伏地表声波正演模拟方法 |
| CN106842320A (zh) * | 2017-01-19 | 2017-06-13 | 北京大学 | Gpu并行三维地震波场生成方法和系统 |
-
2017
- 2017-09-19 CN CN201710845834.6A patent/CN107561585A/zh active Pending
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2013055840A3 (en) * | 2011-10-11 | 2013-10-24 | Saudi Arabian Oil Company | High performance and grid computing with quality of service control |
| WO2017034988A1 (en) * | 2015-08-25 | 2017-03-02 | Saudi Arabian Oil Company | High performance and grid computing of hydrocarbon exploration and production data with partitioning quality of service control |
| CN105572730A (zh) * | 2015-12-15 | 2016-05-11 | 中国科学院地质与地球物理研究所 | 三维复杂结构声波正演方法 |
| CN106646593A (zh) * | 2016-08-30 | 2017-05-10 | 国家超级计算天津中心 | 一种跨节点并行的三维起伏地表声波正演模拟方法 |
| CN106842320A (zh) * | 2017-01-19 | 2017-06-13 | 北京大学 | Gpu并行三维地震波场生成方法和系统 |
Non-Patent Citations (1)
| Title |
|---|
| 宋鹏 等: "基于MPI+OpenMP的三维声波方程正演模拟", 《中国海洋大学学报》 * |
Cited By (17)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108445537A (zh) * | 2018-02-07 | 2018-08-24 | 中国地质大学（武汉） | 基于Spark的叠前地震数据AVO弹性参数反演方法及系统 |
| CN108445537B (zh) * | 2018-02-07 | 2019-05-31 | 中国地质大学（武汉） | 基于Spark的叠前地震数据AVO弹性参数反演方法及系统 |
| CN109342575A (zh) * | 2018-09-06 | 2019-02-15 | 北京大学 | 弹性波偏振系统 |
| CN109342575B (zh) * | 2018-09-06 | 2024-03-22 | 北京大学 | 弹性波偏振系统 |
| CN109490948A (zh) * | 2018-11-16 | 2019-03-19 | 西南石油大学 | 地震声学波动方程矢量并行计算方法 |
| CN109490948B (zh) * | 2018-11-16 | 2020-04-28 | 西南石油大学 | 地震声学波动方程矢量并行计算方法 |
| CN112859162B (zh) * | 2019-11-28 | 2024-04-30 | 中国石油天然气集团有限公司 | 基于伪谱法的三维声波方程空间域正演方法及装置 |
| CN112859162A (zh) * | 2019-11-28 | 2021-05-28 | 中国石油天然气集团有限公司 | 基于伪谱法的三维声波方程空间域正演方法及装置 |
| CN110927784A (zh) * | 2019-11-29 | 2020-03-27 | 中国地质大学（武汉） | 一种基于OpenMP的频率域高斯束格林函数计算方法 |
| CN110927784B (zh) * | 2019-11-29 | 2020-09-29 | 中国地质大学（武汉） | 一种基于OpenMP的频率域高斯束格林函数计算方法 |
| CN111766623B (zh) * | 2020-05-25 | 2023-01-31 | 国家超级计算天津中心 | 一种基于mpi进程拓扑的地震波场三维正演模拟方法 |
| CN111766623A (zh) * | 2020-05-25 | 2020-10-13 | 国家超级计算天津中心 | 一种基于mpi进程拓扑的地震波场三维正演模拟方法 |
| CN112578450A (zh) * | 2020-10-16 | 2021-03-30 | 中国石油大学(华东) | 基于频散介质标量波方程有限差分的正演模拟方法及装置 |
| CN113536638A (zh) * | 2021-07-22 | 2021-10-22 | 北京大学 | 基于间断有限元和交错网格的高精度地震波场模拟方法 |
| CN113536638B (zh) * | 2021-07-22 | 2023-09-22 | 北京大学 | 基于间断有限元和交错网格的高精度地震波场模拟方法 |
| CN113960663A (zh) * | 2021-10-25 | 2022-01-21 | 中国地质大学（北京） | 基于并行计算的三维各向异性衰减正演模拟方法及系统 |
| WO2024087002A1 (en) * | 2022-10-25 | 2024-05-02 | Saudi Arabian Oil Company | Methods and systems for determining attenuated traveltime using parallel processing |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106842320B (zh) | Gpu并行三维地震波场生成方法和系统 | |
| CN107561585A (zh) | 一种多核多节点并行三维地震波场生成方法和系统 | |
| Abdelkhalek et al. | Fast seismic modeling and reverse time migration on a GPU cluster | |
| Crempien et al. | UCSB method for simulation of broadband ground motion from kinematic earthquake sources | |
| Shi et al. | Reverse time migration of 3D vertical seismic profile data | |
| CN104570082B (zh) | 一种基于格林函数表征的全波形反演梯度算子的提取方法 | |
| CN105137486A (zh) | 各向异性介质中弹性波逆时偏移成像方法及其装置 | |
| CN105319581A (zh) | 一种高效的时间域全波形反演方法 | |
| Xu et al. | Implementation of prestack reverse time migration using frequency-domain extrapolation | |
| Raknes et al. | Efficient 3D elastic full-waveform inversion using wavefield reconstruction methods | |
| Xue et al. | An efficient GPU implementation for locating micro-seismic sources using 3D elastic wave time-reversal imaging | |
| CN105911584A (zh) | 一种隐式交错网格有限差分弹性波数值模拟方法及装置 | |
| CN107340540A (zh) | 弹性波场的方向波分解方法、装置以及计算机存储介质 | |
| Ba et al. | Near-fault broadband seismograms synthesis in a stratified transversely isotropic half-space using a semi-analytical frequency-wavenumber method | |
| He et al. | Modeling 3-D elastic wave propagation in TI media using discontinuous Galerkin method on tetrahedral meshes | |
| Chen | Staggered‐grid pseudospectral viscoacoustic wave field simulation in two‐dimensional media | |
| Wenk et al. | Regional wave propagation using the discontinuous Galerkin method | |
| CN114460640B (zh) | 有限差分模拟弹性波全波形反演方法和装置 | |
| Etienne et al. | Three-dimensional acoustic full waveform inversion-algorithm and application to Valhall | |
| Zheng et al. | Finite difference method for first-order velocity-stress equation in body-fitted coordinate system | |
| Song et al. | Multi-scale seismic full waveform inversion in the frequency-domain with a multi-grid method | |
| Guan et al. | Preconditioned Conjugate Gradient Algorithm‐Based 2D Waveform Inversion for Shallow‐Surface Site Characterization | |
| Sarajaervi et al. | Ray-based modeling and imaging in viscoelastic media using graphics processing units | |
| Wu et al. | Least-squares reverse time migration of simultaneous sources with deep-learning-based denoising | |
| Liu et al. | Seismic wavefield simulation with irregular surface topography and graphic processing unit (GPU) implementation |
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 | ||
| WD01 | Invention patent application deemed withdrawn after publication | ||
| WD01 | Invention patent application deemed withdrawn after publication |
Application publication date: 20180109 |