CN105572730B - 三维复杂结构声波正演方法 - Google Patents

三维复杂结构声波正演方法 Download PDF

Info

Publication number
CN105572730B
CN105572730B CN201510931296.3A CN201510931296A CN105572730B CN 105572730 B CN105572730 B CN 105572730B CN 201510931296 A CN201510931296 A CN 201510931296A CN 105572730 B CN105572730 B CN 105572730B
Authority
CN
China
Prior art keywords
mrow
mfrac
msubsup
grid
difference equation
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.)
Expired - Fee Related
Application number
CN201510931296.3A
Other languages
English (en)
Other versions
CN105572730A (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.)
Institute of Geology and Geophysics of CAS
Original Assignee
Institute of Geology and Geophysics of CAS
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 Institute of Geology and Geophysics of CAS filed Critical Institute of Geology and Geophysics of CAS
Priority to CN201510931296.3A priority Critical patent/CN105572730B/zh
Publication of CN105572730A publication Critical patent/CN105572730A/zh
Application granted granted Critical
Publication of CN105572730B publication Critical patent/CN105572730B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Medicines Containing Antibodies Or Antigens For Use As Internal Diagnostic Agents (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明实施例提供一种三维复杂结构声波正演方法，该方法包括：获取改进后的GPU‑CPU异构系统，该改进后的GPU‑CPU异构系统包括以下改进内容：将波动方程离散为6个差分方程，包括计算X轴方向的应力的第一差分方程、计算Y轴方向的应力的第二差分方程、计算Z轴方向的应力的第三差分方程、计算X轴方向的质点运动速度的第四差分方程、计算Y轴方向的质点运动速度的第五差分方程以及计算Z轴方向的质点运动速度的第六差分方程；在GPU内设计4个计算程序分别用于承载第一差分方程、第二差分方程、第三差分方程、第四差分方程、第五差分方程以及第六差分方程的计算；采用改进后的GPU‑CPU异构系统对待处理的地震数据进行正演。

Description

三维复杂结构声波正演方法
技术领域
本发明涉及地震勘探技术领域，特别涉及一种三维复杂结构声波正演方法。
背景技术
石油和天然气的产量关乎国计民生，对地下构造的精细认识是指导油气勘探的关键环节，现今主要的方法是地震勘探。通过观测和分析人工地震产生的地震波对地下构造高精度成像，可有效得到地层、构造的分布，进而估算油气层厚度及分布范围。
地震勘探包含数据采集、资料处理和资料解释三个环节，数据处理是将地震数据变成地震剖面或者构造图，确定地下构造，指出有利的含油气地区。地震数据处理需要进行大量数据运算，对计算效率要求极高。
地震数据处理中包含去噪、反褶积、动校正、静校正、速度分析、叠加、偏移成像。地震偏移成像是地震数据处理的核心，也是计算量最大的环节，通常也是实际生产的瓶颈。
近年来，偏移成像技术不断发展，由人工偏移转变到计算机数字化偏移，由射线偏移延伸到波动方程偏移，由叠后偏移拓展到叠前偏移，由二维扩展到三维等。
逆时偏移能够对任意类型的波进行模拟成像，不存在倾角限制，并适用于任意复杂介质，成为目前最精确的偏移成像技术之一。其原理是对双程波动方程进行精确求解，相对于射线方法，不做过分近似。逆时偏移主要的计算量在于炮点的正传和检波点的正传，相当于两次正演过程，因此三维复杂结构正演的效率决定了逆时偏移计算的效率。目前通过大规模CPU集群以及OpenMP(用于共享内存并行系统的多处理器程序设计的一套指导性的编译处理方案)、MPI(Message Passing Interface，消息传递接口)等并行实现节点并行，可以很大程度上提高运算效率，但在巨大计算需求的逆时偏移计算面前仍出现较多难题。因此通过GPU(Graphics Processing Unit，图像处理器)-CPU异构系统进行三维复杂结构声波正演加速显得尤为重要。
近些年，全波形反演作为一种高精度速度模型构建的方法得到人们的普遍关注，三维全波形反演技术中正演模拟部分计算量是限制其大规模应用的瓶颈之一，因此加速优化三维复杂结构声波正演对三维全波形反演技术也至关重要。
目前主要通过CPU集群以及OpenMP、MPI并行实现节点并行，也有使用GPU-CPU异构系统加速，但未见使用较好的加速策略来使得GPU-CPU异构系统在三维复杂结构声波正演中取得较高加速比。
发明内容
本发明实施例提供了一种三维复杂结构声波正演方法，以解决现有技术中使用GPU-CPU异构系统加速但加速比低的技术问题。该方法包括：获取改进后的图形处理器GPU-CPU异构系统，其中，该改进后的GPU-CPU异构系统包括以下改进内容：将波动方程离散为6个差分方程，该6个差分方程包括一个用于计算X轴方向的应力的第一差分方程、一个用于计算Y轴方向的应力的第二差分方程、一个用于计算Z轴方向的应力的第三差分方程、一个用于计算X轴方向的质点运动速度的第四差分方程、一个用于计算Y轴方向的质点运动速度的第五差分方程以及一个用于计算Z轴方向的质点运动速度的第六差分方程；在GPU内设计4个计算程序，其中，一个计算程序用于承载第一差分方程的计算，一个计算程序用于承载第二差分方程的计算，一个计算程序用于承载第三差分方程的计算，一个计算程序用于承载第四差分方程、第五差分方程以及第六差分方程的计算；采用所述改进后的GPU-CPU异构系统对待处理的地震数据进行正演。
在一个实施例中，所述第一差分方程为以下形式：
所述第二差分方程为以下形式：
所述第三差分方程为以下形式：
所述第四差分方程为以下形式：
所述第五差分方程为以下形式：
所述第六差分方程为以下形式：
其中，vx为X轴方向的质点运动速度；vy为Y轴方向的质点运动速度；vz为Z轴方向的质点运动速度；τx为X轴方向的应力；τy为Y轴方向的应力；τz为Z轴方向的应力；ρ代表介质密度，取值为1；am为常量；Δt为模拟过程中的时间间隔；Δx为x轴方向网格间隔；Δy为y轴方向网格间隔；Δz轴为z方向网格间隔；i,j,k代表网格点的坐标位置；x,y,z分别代表X轴，Y轴和Z轴方向；m为计算一个网格点所需要的基准网格点个数；L为计算一个网格点所需要的基准网格点的最大个数；t为时刻，t+1为t+1时刻，表示向后推一个时刻，t-1代表t-1时刻，表示向前推一个时刻，代表时刻，表示后推半个时刻，代表时刻，表示向前推半个时刻；d(x)是X轴方向的衰减系数，d(x)＝d0(i/pML)n，n为指数阶数，n为自然数，d0＝(ln(1/R))(n+1)v/(2pML)，i代表从区域边界到三维速度模型最内层的网格点数，pML为匹配层的网格厚度；d(y)是Y轴方向的衰减系数，d(y)＝d0(i/pML)n；d(z)是Z轴方向的衰减系数，d(z)＝d0(i/pML)n；v为网格点的声波速度值；τ是网格点的应力。
在一个实施例中，设计用于承载第一差分方程计算的计算程序，包括：将输入的待处理的地震数据按照X轴方向连续的方式进行存储，数据块的大小为X轴方向上待处理的网格数；采用连续编号的GPU线程计算X轴方向连续的网格的应力，网格与线程一一对应，通过与当前计算网格在X轴方向上相邻的网格的速度分量来计算当前计算网格的应力。
在一个实施例中，设计用于承载第二差分方程计算的计算程序，包括：将输入的待处理的地震数据按照X轴方向连续的方式进行存储，将输入的待处理的地震数据在X轴和Z轴组成的平面上划分为二维数据块；采用连续编号的GPU线程计算Y轴方向的网格的应力，网格与线程一一对应，通过与当前计算网格在Y轴方向上相邻的网格的速度分量来计算当前计算网格的应力。
在一个实施例中，设计用于承载第四差分方程、第五差分方程以及第六差分方程计算的计算程序，包括：将输入的待处理的地震数据按照X轴方向连续的方式进行存储，将输入的待处理的地震数据在X轴和Y轴组成的平面上划分为二维数据块；采用连续编号的GPU线程分别计算X轴、Y轴以及Z轴方向的网格的质点运动速度，网格与线程一一对应，通过与当前计算网格在X轴、Y轴以及Z轴方向上相邻的网格的应力来计算当前计算网格的质点运动速度。
在一个实施例中，在采用所述改进后的GPU-CPU异构系统对待处理的地震数据进行正演之前，还包括：对所述待处理的地震数据，在X轴方向数据的首地址处填充字节，填充后的X轴方向数据的字节数是32的倍数；在所述待处理的地震数据的头部填充数据，填充后的待处理的地震数据的首地址按128字节对齐。
在一个实施例中，通过以下公式计算在X轴方向数据的首地址处填充字节的长度：pad_x＝32-(nx+2pML+2r)％32；其中，pad_x是在X轴方向数据的首地址处填充字节的长度，nx是X轴方向的网格数，pML是匹配层的网格厚度，r是计算当前计算网格时需要的相邻网格层数，r＝d/2，d是空间导数离散精度。
在一个实施例中，通过以下公式计算在所述待处理的地震数据的头部填充数据的大小：lead＝32-pad_x-r；其中，lead是在所述待处理的地震数据的头部填充的数据大小，pad_x是在X轴方向数据的首地址处填充字节的长度，r是计算当前计算网格时需要的相邻网格层数，r＝d/2，d是空间导数离散精度。
在一个实施例中，在采用所述改进后的GPU-CPU异构系统对待处理的地震数据进行正演之后，还包括：实时分析4个所述计算程序对计算资源和内存带宽的使用情况，根据所述使用情况调整4个所述计算程序的计算方式。
在本发明实施例中，通过对GPU-CPU异构系统进行改进，将波动方程离散为6个差分方程，并设计4个计算程序来承载上述6个差分方程的计算，由于计算程序的数目小，降低了启动计算程序带来的开销，同时，采用一个计算程序来承载计算速度的三个差分方程的计算，实现可以重用应力变量，减少对全局内存的访问次数，从而提高了GPU-CPU异构系统的计算效率，进而提高了三维复杂结构声波正演的计算效率，对石油天然气勘探的实际生产意义重大，有助于实现三维全波形反演技术的实现。
附图说明
此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，并不构成对本发明的限定。在附图中：
图1是本发明实施例提供的一种三维复杂结构声波正演方法的流程图；
图2是本发明实施例提供的一种数据块划分示意图；
图3是本发明实施例提供的一种数据填充示意图。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚明白，下面结合实施方式和附图，对本发明做进一步详细说明。在此，本发明的示意性实施方式及其说明用于解释本发明，但并不作为对本发明的限定。
在本发明实施例中，提供了一种三维复杂结构声波正演方法，如图1所示，该方法包括：
步骤101：获取改进后的图形处理器GPU-CPU异构系统，其中，该改进后的GPU-CPU异构系统包括以下改进内容：
将波动方程离散为6个差分方程，该6个差分方程包括一个用于计算X轴方向的应力的第一差分方程、一个用于计算Y轴方向的应力的第二差分方程、一个用于计算Z轴方向的应力的第三差分方程、一个用于计算X轴方向的质点运动速度的第四差分方程、一个用于计算Y轴方向的质点运动速度的第五差分方程以及一个用于计算Z轴方向的质点运动速度的第六差分方程；
在GPU内设计4个计算程序，其中，一个计算程序用于承载第一差分方程的计算，一个计算程序用于承载第二差分方程的计算，一个计算程序用于承载第三差分方程的计算，一个计算程序用于承载第四差分方程、第五差分方程以及第六差分方程的计算；
步骤102：采用所述改进后的GPU-CPU异构系统对待处理的地震数据进行正演。
由图1所示的流程可知，在本发明实施例中，通过对GPU-CPU异构系统进行改进，将波动方程离散为6个差分方程，并设计4个计算程序来承载上述6个差分方程的计算，由于计算程序的数目小，降低了启动计算程序带来的开销，同时，采用一个计算程序来承载计算速度的三个差分方程的计算，实现可以重用应力变量，减少对全局内存的访问次数，从而提高了GPU-CPU异构系统的计算效率，进而提高了三维复杂结构声波正演的计算效率，对石油天然气勘探的实际生产意义重大，有助于实现三维全波形反演技术的实现。
具体实施时，所述第一差分方程为以下形式：
所述第二差分方程为以下形式：
所述第三差分方程为以下形式：
所述第四差分方程为以下形式：
所述第五差分方程为以下形式：
所述第六差分方程为以下形式：
其中，vx为X轴方向的质点运动速度；vy为Y轴方向的质点运动速度；vz为Z轴方向的质点运动速度；τx为X轴方向的应力；τy为Y轴方向的应力；τz为Z轴方向的应力；ρ代表介质密度，取值为1；am为常量；Δt为模拟过程中的时间间隔；Δx为x轴方向网格间隔；Δy为y轴方向网格间隔；Δz轴为z方向网格间隔；i,j,k代表网格点的坐标位置；x,y,z分别代表X轴，Y轴和Z轴方向；m为计算一个网格点所需要的基准网格点个数；L为计算一个网格点所需要的基准网格点的最大个数；t为时刻，t+1为t+1时刻，表示向后推一个时刻，t-1代表t-1时刻，表示向前推一个时刻，代表时刻，表示后推半个时刻，代表时刻，表示向前推半个时刻；d(x)是X轴方向的衰减系数，d(x)＝d0(i/pML)n，n为指数阶数，n为自然数，d0＝(ln(1/R))(n+1)v/(2pML)，i代表从区域边界到三维速度模型最内层的网格点数，pML为匹配层的网格厚度；d(y)是Y轴方向的衰减系数，d(y)＝d0(i/pML)n；d(z)是Z轴方向的衰减系数，d(z)＝d0(i/pML)n；v为网格点的声波速度值；τ是网格点的应力。
具体实施时，考虑到GPU-CPU异构系统的特点，GPU-CPU异构系统具有更强大的并行能力。其编程语言为类似C语法的CUDA语言，具体的可参考英伟达公司发布的CUDA手册。在GPU上可执行的程序称为kernel，GPU卡上有不同种类的内存，如下表1所示，例如，Constant Memory(常数存储器)存在于显存，只读，可被grid中所有thread读取；SharedMemory(共享内存)是最重要的内存，类似于CPU中的cache，可读写，非常快，3.X硬件每个SM16KB，只能被thread block内的线程读取；Register(寄存器)3.0硬件有63个寄存器，thread私有；Global Memory(全局内存)访问最慢(400-600周期)。因此在利用内存时，要尽可能一次读入数据做多次运算，重复读写对GPU延迟很严重。
表1
为了提高GPU-CPU异构系统的计算效率，关键就是要隐藏用于存储的时间延迟，可以通过大量任务同时并发以及减少对全局内存的访问来实现。因此，在本实施例中，提出设计上述4个计算程序(即kernel)的实现方法，以实现4个计算程序采取不同的GPU内存使用策略，实现对全局内存的访存满足连续访问，例如，设计用于承载第一差分方程计算的计算程序，包括：
将输入的待处理的地震数据按照X轴方向连续的方式进行存储，数据块的大小为X轴方向上待处理的网格数；即为了尽量减小边界数据的比例，采用一维的Block(数据块)，Block大小等于X轴方向待处理的网格数，自然地，该kernel的Grid大小为(ny+2pML,nz+2pML)；
采用连续编号的GPU线程计算X轴方向连续的网格的应力，网格与线程一一对应，通过与当前计算网格在X轴方向上相邻的网格的速度分量来计算当前计算网格的应力。即通过使用连续编号的GPU线程处理X轴方向连续的网格的应力，以满足对全局内存的连续访问，具体实现时，每个Block把对应的数据从全局内存加载到共享内存，后面计算需要用到的数据都从共享内存中读取。
具体实施时，计算Y轴方向的应力分量时，由于数据在Y轴方向的存储不连续，此时，如果仍然采用与计算第一差分方程的计算程序一样的实现策略，不能满足对全局内存的连续访问，会极大地降低GPU上的性能。因此，在本实施例中，设计用于承载第二差分方程计算的计算程序，包括：
将输入的待处理的地震数据按照X轴方向连续的方式进行存储，将输入的待处理的地震数据在X轴和Z轴组成的平面上划分为二维数据块；即把输入的待处理的地震数据在X轴和Z轴组成的平面上划分为二维数据块，如图2所示(左边是Grid的示意图，右边是Grid中的一个Block的放大示意图)，Grid和Block都设计成二维。可先确定Block的大小，假定为(Bx,Bz)，则Grid的大小为((nx+2pML)/Bx,(nz+2pML)/Bz)，此时，连续的线程访问的也是连续的内存空间；
采用连续编号的GPU线程计算Y轴方向的网格的应力，网格与线程一一对应，通过与当前计算网格在Y轴方向上相邻的网格的速度分量来计算当前计算网格的应力。
设计用于承载第三差分方程计算的计算程序可以采用设计用于承载第二差分方程计算的计算程序的方法，例如，将输入的待处理的地震数据按照X轴方向连续的方式进行存储，将输入的待处理的地震数据在X轴和Y轴组成的平面上划分为二维数据块；采用连续编号的GPU线程计算Z轴方向的网格的应力，网格与线程一一对应，通过与当前计算网格在Z轴方向上相邻的网格的速度分量来计算当前计算网格的应力。
具体实施时，通过以下方法设计用于承载第四差分方程、第五差分方程以及第六差分方程计算的计算程序：将输入的待处理的地震数据按照X轴方向连续的方式进行存储，将输入的待处理的地震数据在X轴和Y轴组成的平面上划分为二维数据块；采用连续编号的GPU线程分别计算X轴、Y轴以及Z轴方向的网格的质点运动速度，网格与线程一一对应，通过与当前计算网格在X轴、Y轴以及Z轴方向上相邻的网格的应力来计算当前计算网格的质点运动速度。
通过上述方法设计上述4个计算程序后，GPU-CPU异构系统相对于同价位CPU集群可以达到450倍的加速。
上述4个kernel(计算程序)设计时，保证了对全局内存的访存满足连续访问，因此，理论上对全局内存的写操作效率能达到100％。对于全局内存的读操作，由于要加载边界数据，理论上读数据效率尽管不能达到100％，但也会非常接近。应用NVIDIA的性能分析工具nvvp分析该kernel，写数据实际达到的效率只有80％，没有达到理论值。
假定计算规模nx＝ny＝nz＝256，匹配层的网格厚度pML＝32，空间导数离散精度为d阶，则需要的邻居网格层数为r＝d/2。实际计算时，还需要在每个方向增加r层网格，如图3(a)所示。此时，X轴方向实际的网格数为nx+2pML+2r，转化为用字节表示为4(nx+2pML+2r)，它通常不是128的倍数。在GPU上，每连续的32个线程组织成一个线程束(warp)同时执行。假如每个线程需要写一个单精度浮点数，即4个字节，整个warp写的数据为128字节。此时，要达到理想的合并访问，即一个读指令，就完成对128字节的写操作，除了要求这128个字节在内存中是连续的以外(这在前面设计kernel时已经得到满足)，还要求目的数据的首地址按128字节对齐。三维数据体在X轴方向连续存储，要求X轴方向数据所占字节为128的整数倍，以保证在Y或者Z平面跳转时，跳转字节长度为128字节的倍数。因此，首先，在采用所述改进后的GPU-CPU异构系统对待处理的地震数据进行正演之前，还包括：对待处理的地震数据，在X轴方向数据的首地址处填充字节，转化为数据单位，填充后的X轴方向数据的字节数是32的倍数，以确保跳转长度为128字节的倍数。具体的，通过以下公式计算在X轴方向数据的首地址处填充字节的长度：
pad_x＝32-(nx+2pML+2r)％32；
其中，pad_x是在X轴方向数据的首地址处填充字节的长度，nx是X轴方向的网格数，pML是匹配层的网格厚度，r是计算当前计算网格时需要的相邻网格层数，r＝d/2，d是空间导数离散精度。
在X轴方向数据的首地址处填充字节后，如图3(b)所示，以浮点数据为单位，保证跳转长度为128字节的倍数。刚才前面已经指出，合并访问的条件之一是首地址按128字节对齐，显然该首地址不符合这一要求。因此，需要在待处理的地震数据的头部填充数据，填充后的待处理的地震数据的首地址按128字节对齐，填充后的待处理的地震数据如图3(c)所示，具体的，通过以下公式计算在所述待处理的地震数据的头部填充数据的大小：
lead＝32-pad_x-r；
其中，lead是在所述待处理的地震数据的头部填充的数据大小，pad_x是在X轴方向数据的首地址处填充字节的长度，r是计算当前计算网格时需要的相邻网格层数，r＝d/2，d是空间导数离散精度。
通过数据填充优化全局内存的访存后，程序性能提高到1.46Gcell/s，相比上个设计计算程序的优化策略GPU实现达到1.32倍加速，相比同价位CPU实现达到608倍加速。全局内存写数据效率提高到100％，读数据效率达到98.6％。
具体实施时，在采用所述改进后的GPU-CPU异构系统对待处理的地震数据进行正演之后，还包括：实时分析4个所述计算程序对计算资源和内存带宽的使用情况，根据所述使用情况调整4个所述计算程序的计算方式。可以使用nvvp分析上述计算程序对计算资源和内存带宽的使用情况，以分析出此时的瓶颈，例如，此时计算资源和内存带宽使用都小于60％，瓶颈在于延迟，目前该kernel的占有率已经达到50％，由SM的寄存器总量决定。此时调整采取的计算策略，由每个线程一次处理一个网格的计算改为同时处理X轴方向相邻的两个网格点的计算。这一优化后，程序性能达到1.70Gcell/s。相应地，内存带宽的利用率由原来的50％提高到75％。
在本发明实施例中，通过对GPU-CPU异构系统进行改进，将波动方程离散为6个差分方程，并设计4个计算程序来承载上述6个差分方程的计算，由于计算程序的数目小，降低了启动计算程序带来的开销，同时，采用一个计算程序来承载计算速度的三个差分方程的计算，实现可以重用应力变量，减少对全局内存的访问次数，从而提高了GPU-CPU异构系统的计算效率，进而提高了三维复杂结构声波正演的计算效率，对石油天然气勘探的实际生产意义重大，有助于实现三维全波形反演技术的实现。
显然，本领域的技术人员应该明白，上述的本发明实施例的各模块或各步骤可以用通用的计算装置来实现，它们可以集中在单个的计算装置上，或者分布在多个计算装置所组成的网络上，可选地，它们可以用计算装置可执行的程序代码来实现，从而，可以将它们存储在存储装置中由计算装置来执行，并且在某些情况下，可以以不同于此处的顺序执行所示出或描述的步骤，或者将它们分别制作成各个集成电路模块，或者将它们中的多个模块或步骤制作成单个集成电路模块来实现。这样，本发明实施例不限制于任何特定的硬件和软件结合。
以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明实施例可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (9)

1.一种三维复杂结构声波正演方法，其特征在于，包括：
获取改进后的图形处理器GPU-CPU异构系统，其中，该改进后的GPU-CPU异构系统包括以下改进内容：
将波动方程离散为6个差分方程，该6个差分方程包括一个用于计算X轴方向的应力的第一差分方程、一个用于计算Y轴方向的应力的第二差分方程、一个用于计算Z轴方向的应力的第三差分方程、一个用于计算X轴方向的质点运动速度的第四差分方程、一个用于计算Y轴方向的质点运动速度的第五差分方程以及一个用于计算Z轴方向的质点运动速度的第六差分方程；
在GPU内设计4个计算程序，其中，一个计算程序用于承载第一差分方程的计算，一个计算程序用于承载第二差分方程的计算，一个计算程序用于承载第三差分方程的计算，一个计算程序用于承载第四差分方程、第五差分方程以及第六差分方程的计算；
采用所述改进后的GPU-CPU异构系统对待处理的地震数据进行正演。
2.如权利要求1所述的三维复杂结构声波正演方法，其特征在于，所述第一差分方程为以下形式：
<mrow> <mfrac> <mrow> <msubsup> <mi>&amp;tau;</mi> <mrow> <mi>x</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mi>t</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> </mrow> </msubsup> <mo>-</mo> <msubsup> <mi>&amp;tau;</mi> <mrow> <mi>x</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mi>t</mi> <mo>-</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> </mrow> </msubsup> </mrow> <mrow> <mi>&amp;Delta;</mi> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <mi>d</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo> </mrow> <mfrac> <mrow> <msubsup> <mi>&amp;tau;</mi> <mrow> <mi>x</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mi>t</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> </mrow> </msubsup> <mo>-</mo> <msubsup> <mi>&amp;tau;</mi> <mrow> <mi>x</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mi>t</mi> <mo>-</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> </mrow> </msubsup> </mrow> <mn>2</mn> </mfrac> <mo>=</mo> <mo>-</mo> <msup> <mi>&amp;rho;v</mi> <mn>2</mn> </msup> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>L</mi> </munderover> <msub> <mi>a</mi> <mi>m</mi> </msub> <mfrac> <mrow> <msubsup> <mi>v</mi> <mrow> <mi>x</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mo>+</mo> <mi>m</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> <mi>t</mi> </msubsup> <mo>-</mo> <msubsup> <mi>v</mi> <mrow> <mi>x</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>-</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mo>-</mo> <mi>m</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> <mi>t</mi> </msubsup> </mrow> <mrow> <mi>&amp;Delta;</mi> <mi>x</mi> </mrow> </mfrac> <mo>;</mo> </mrow>
所述第二差分方程为以下形式：
<mrow> <mfrac> <mrow> <msubsup> <mi>&amp;tau;</mi> <mrow> <mi>y</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mi>t</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> </mrow> </msubsup> <mo>-</mo> <msubsup> <mi>&amp;tau;</mi> <mrow> <mi>y</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mi>t</mi> <mo>-</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> </mrow> </msubsup> </mrow> <mrow> <mi>&amp;Delta;</mi> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <mi>d</mi> <mrow> <mo>(</mo> <mi>y</mi> <mo>)</mo> </mrow> <mfrac> <mrow> <msubsup> <mi>&amp;tau;</mi> <mrow> <mi>y</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mi>t</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> </mrow> </msubsup> <mo>-</mo> <msubsup> <mi>&amp;tau;</mi> <mrow> <mi>y</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mi>t</mi> <mo>-</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> </mrow> </msubsup> </mrow> <mn>2</mn> </mfrac> <mo>=</mo> <mo>-</mo> <msup> <mi>&amp;rho;v</mi> <mn>2</mn> </msup> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>L</mi> </munderover> <msub> <mi>a</mi> <mi>m</mi> </msub> <mfrac> <mrow> <msubsup> <mi>v</mi> <mrow> <mi>y</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mo>+</mo> <mi>m</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> <mi>t</mi> </msubsup> <mo>-</mo> <msubsup> <mi>v</mi> <mrow> <mi>y</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>-</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mo>-</mo> <mi>m</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> <mi>t</mi> </msubsup> </mrow> <mrow> <mi>&amp;Delta;</mi> <mi>y</mi> </mrow> </mfrac> <mo>;</mo> </mrow>
所述第三差分方程为以下形式：
<mrow> <mfrac> <mrow> <msubsup> <mi>&amp;tau;</mi> <mrow> <mi>z</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mi>t</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> </mrow> </msubsup> <mo>-</mo> <msubsup> <mi>&amp;tau;</mi> <mrow> <mi>z</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mi>t</mi> <mo>-</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> </mrow> </msubsup> </mrow> <mrow> <mi>&amp;Delta;</mi> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <mi>d</mi> <mrow> <mo>(</mo> <mi>z</mi> <mo>)</mo> </mrow> <mfrac> <mrow> <msubsup> <mi>&amp;tau;</mi> <mrow> <mi>y</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mi>t</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> </mrow> </msubsup> <mo>-</mo> <msubsup> <mi>&amp;tau;</mi> <mrow> <mi>y</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mi>t</mi> <mo>-</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> </mrow> </msubsup> </mrow> <mn>2</mn> </mfrac> <mo>=</mo> <mo>-</mo> <msup> <mi>&amp;rho;v</mi> <mn>2</mn> </msup> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>L</mi> </munderover> <msub> <mi>a</mi> <mi>m</mi> </msub> <mfrac> <mrow> <msubsup> <mi>v</mi> <mrow> <mi>z</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mo>+</mo> <mi>m</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> <mi>t</mi> </msubsup> <mo>-</mo> <msubsup> <mi>v</mi> <mrow> <mi>z</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>-</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mo>-</mo> <mi>m</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> <mi>t</mi> </msubsup> </mrow> <mrow> <mi>&amp;Delta;</mi> <mi>z</mi> </mrow> </mfrac> <mo>;</mo> </mrow>
所述第四差分方程为以下形式：
<mrow> <mfrac> <mrow> <msubsup> <mi>v</mi> <mrow> <mi>x</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <mo>-</mo> <msubsup> <mi>v</mi> <mrow> <mi>x</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> <mi>t</mi> </msubsup> </mrow> <mrow> <mi>&amp;Delta;</mi> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <mi>d</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo> </mrow> <mfrac> <mrow> <msubsup> <mi>v</mi> <mrow> <mi>x</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <mo>-</mo> <msubsup> <mi>v</mi> <mrow> <mi>x</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> <mi>t</mi> </msubsup> </mrow> <mn>2</mn> </mfrac> <mo>=</mo> <mo>-</mo> <mfrac> <mn>1</mn> <mi>&amp;rho;</mi> </mfrac> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>L</mi> </munderover> <msub> <mi>a</mi> <mi>m</mi> </msub> <mfrac> <mrow> <msubsup> <mi>&amp;tau;</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>+</mo> <mn>1</mn> <mo>+</mo> <mi>m</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> <mrow> <mi>t</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> </mrow> </msubsup> <mo>-</mo> <msubsup> <mi>&amp;tau;</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>-</mo> <mi>m</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> <mrow> <mi>t</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> </mrow> </msubsup> </mrow> <mrow> <mi>&amp;Delta;</mi> <mi>x</mi> </mrow> </mfrac> <mo>;</mo> </mrow>
所述第五差分方程为以下形式：
<mrow> <mfrac> <mrow> <msubsup> <mi>v</mi> <mrow> <mi>x</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <mo>-</mo> <msubsup> <mi>v</mi> <mrow> <mi>x</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> <mi>t</mi> </msubsup> </mrow> <mrow> <mi>&amp;Delta;</mi> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <mi>d</mi> <mrow> <mo>(</mo> <mi>y</mi> <mo>)</mo> </mrow> <mfrac> <mrow> <msubsup> <mi>v</mi> <mrow> <mi>x</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <mo>-</mo> <msubsup> <mi>v</mi> <mrow> <mi>x</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> <mi>t</mi> </msubsup> </mrow> <mn>2</mn> </mfrac> <mo>=</mo> <mo>-</mo> <mfrac> <mn>1</mn> <mi>&amp;rho;</mi> </mfrac> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>L</mi> </munderover> <msub> <mi>a</mi> <mi>m</mi> </msub> <mfrac> <mrow> <msubsup> <mi>&amp;tau;</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>+</mo> <mn>1</mn> <mo>+</mo> <mi>m</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> <mrow> <mi>t</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> </mrow> </msubsup> <mo>-</mo> <msubsup> <mi>&amp;tau;</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>-</mo> <mi>m</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> <mrow> <mi>t</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> </mrow> </msubsup> </mrow> <mrow> <mi>&amp;Delta;</mi> <mi>y</mi> </mrow> </mfrac> <mo>;</mo> </mrow>
所述第六差分方程为以下形式：
<mrow> <mfrac> <mrow> <msubsup> <mi>v</mi> <mrow> <mi>x</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> </mrow> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <mo>-</mo> <msubsup> <mi>v</mi> <mrow> <mi>x</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> </mrow> <mi>t</mi> </msubsup> </mrow> <mrow> <mi>&amp;Delta;</mi> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <mi>d</mi> <mrow> <mo>(</mo> <mi>z</mi> <mo>)</mo> </mrow> <mfrac> <mrow> <msubsup> <mi>v</mi> <mrow> <mi>x</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> </mrow> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <mo>-</mo> <msubsup> <mi>v</mi> <mrow> <mi>x</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> </mrow> <mi>t</mi> </msubsup> </mrow> <mn>2</mn> </mfrac> <mo>=</mo> <mo>-</mo> <mfrac> <mn>1</mn> <mi>&amp;rho;</mi> </mfrac> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>L</mi> </munderover> <msub> <mi>a</mi> <mi>m</mi> </msub> <mfrac> <mrow> <msubsup> <mi>&amp;tau;</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>+</mo> <mi>m</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> <mrow> <mi>t</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> </mrow> </msubsup> <mo>-</mo> <msubsup> <mi>&amp;tau;</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>-</mo> <mi>m</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mrow> <mi>t</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> </mrow> </msubsup> </mrow> <mrow> <mi>&amp;Delta;</mi> <mi>z</mi> </mrow> </mfrac> <mo>;</mo> </mrow>
其中，vx为X轴方向的质点运动速度；vy为Y轴方向的质点运动速度；vz为Z轴方向的质点运动速度；τx为X轴方向的应力；τy为Y轴方向的应力；τz为Z轴方向的应力；ρ代表介质密度，取值为1；am为常量；Δt为模拟过程中的时间间隔；Δx为x轴方向网格间隔；Δy为y轴方向网格间隔；Δz轴为z方向网格间隔；i,j,k代表网格点的坐标位置；x,y,z分别代表X轴，Y轴和Z轴方向；m为计算一个网格点所需要的基准网格点个数；L为计算一个网格点所需要的基准网格点的最大个数；t为时刻，t+1为t+1时刻，表示向后推一个时刻，t-1代表t-1时刻，表示向前推一个时刻，代表时刻，表示后推半个时刻，代表时刻，表示向前推半个时刻；d(x)是X轴方向的衰减系数，d(x)＝d0(i/pML)n，n为指数阶数，n为自然数，d0＝(ln(1/R))(n+1)v/(2pML)，i代表从区域边界到三维速度模型最内层的网格点数，pML为匹配层的网格厚度；d(y)是Y轴方向的衰减系数，d(y)＝d0(i/pML)n；d(z)是Z轴方向的衰减系数，d(z)＝d0(i/pML)n；v为网格点的声波速度值；τ是网格点的应力。
3.如权利要求1所述的三维复杂结构声波正演方法，其特征在于，设计用于承载第一差分方程计算的计算程序，包括：
将输入的待处理的地震数据按照X轴方向连续的方式进行存储，数据块的大小为X轴方向上待处理的网格数；
采用连续编号的GPU线程计算X轴方向连续的网格的应力，网格与线程一一对应，通过与当前计算网格在X轴方向上相邻的网格的速度分量来计算当前计算网格的应力。
4.如权利要求1所述的三维复杂结构声波正演方法，其特征在于，设计用于承载第二差分方程计算的计算程序，包括：
将输入的待处理的地震数据按照X轴方向连续的方式进行存储，将输入的待处理的地震数据在X轴和Z轴组成的平面上划分为二维数据块；
采用连续编号的GPU线程计算Y轴方向的网格的应力，网格与线程一一对应，通过与当前计算网格在Y轴方向上相邻的网格的速度分量来计算当前计算网格的应力。
5.如权利要求1所述的三维复杂结构声波正演方法，其特征在于，设计用于承载第四差分方程、第五差分方程以及第六差分方程计算的计算程序，包括：
将输入的待处理的地震数据按照X轴方向连续的方式进行存储，将输入的待处理的地震数据在X轴和Y轴组成的平面上划分为二维数据块；
采用连续编号的GPU线程分别计算X轴、Y轴以及Z轴方向的网格的质点运动速度，网格与线程一一对应，通过与当前计算网格在X轴、Y轴以及Z轴方向上相邻的网格的应力来计算当前计算网格的质点运动速度。
6.如权利要求1至5中任一项所述的三维复杂结构声波正演方法，其特征在于，在采用所述改进后的GPU-CPU异构系统对待处理的地震数据进行正演之前，还包括：
对所述待处理的地震数据，在X轴方向数据的首地址处填充字节，填充后的X轴方向数据的字节数是32的倍数；
在所述待处理的地震数据的头部填充数据，填充后的待处理的地震数据的首地址按128字节对齐。
7.如权利要求6所述的三维复杂结构声波正演方法，其特征在于，通过以下公式计算在X轴方向数据的首地址处填充字节的长度：
pad_x＝32-(nx+2pML+2r)％32
其中，pad_x是在X轴方向数据的首地址处填充字节的长度，nx是X轴方向的网格数，pML是匹配层的网格厚度，r是计算当前计算网格时需要的相邻网格层数，r＝d/2，d是空间导数离散精度。
8.如权利要求6所述的三维复杂结构声波正演方法，其特征在于，通过以下公式计算在所述待处理的地震数据的头部填充数据的大小：
lead＝32-pad_x-r
其中，lead是在所述待处理的地震数据的头部填充的数据大小，pad_x是在X轴方向数据的首地址处填充字节的长度，r是计算当前计算网格时需要的相邻网格层数，r＝d/2，d是空间导数离散精度。
9.如权利要求1至5中任一项所述的三维复杂结构声波正演方法，其特征在于，在采用所述改进后的GPU-CPU异构系统对待处理的地震数据进行正演之后，还包括：
实时分析4个所述计算程序对计算资源和内存带宽的使用情况，根据所述使用情况调整4个所述计算程序的计算方式。
CN201510931296.3A 2015-12-15 2015-12-15 三维复杂结构声波正演方法 Expired - Fee Related CN105572730B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201510931296.3A CN105572730B (zh) 2015-12-15 2015-12-15 三维复杂结构声波正演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201510931296.3A CN105572730B (zh) 2015-12-15 2015-12-15 三维复杂结构声波正演方法

Publications (2)

Publication Number Publication Date
CN105572730A CN105572730A (zh) 2016-05-11
CN105572730B true CN105572730B (zh) 2017-11-14

Family

ID=55883071

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201510931296.3A Expired - Fee Related CN105572730B (zh) 2015-12-15 2015-12-15 三维复杂结构声波正演方法

Country Status (1)

Country Link
CN (1) CN105572730B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107561585A (zh) * 2017-09-19 2018-01-09 北京大学 一种多核多节点并行三维地震波场生成方法和系统
CN112099936A (zh) * 2019-06-17 2020-12-18 中国石油天然气集团有限公司 三维声波npml算法的异构并行计算实现方法及装置

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102207987A (zh) * 2011-05-31 2011-10-05 中国航天标准化研究所 基于OpenCL的GPU加速三维时域有限差分电磁场仿真的方法
WO2011141440A1 (en) * 2010-05-12 2011-11-17 Shell Internationale Research Maatschappij B.V. Seismic p-wave modelling in an inhomogeneous transversely isotropic medium with a tilted symmetry axis
CN102353988A (zh) * 2011-07-08 2012-02-15 中国科学院地质与地球物理研究所 基于图形处理器计算起伏地表直接叠前逆时偏移的方法
CN102636808A (zh) * 2012-02-07 2012-08-15 浪潮(北京)电子信息产业有限公司 一种地震叠前时间偏移走时处理方法及装置
CN103675908A (zh) * 2012-09-21 2014-03-26 中国石油化工股份有限公司 一种海量数据图形处理器的波动方程逆时偏移成像方法
CN105005072A (zh) * 2015-06-02 2015-10-28 中国科学院地质与地球物理研究所 利用cuda的pml边界三维地震波传播模拟方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2011141440A1 (en) * 2010-05-12 2011-11-17 Shell Internationale Research Maatschappij B.V. Seismic p-wave modelling in an inhomogeneous transversely isotropic medium with a tilted symmetry axis
CN102207987A (zh) * 2011-05-31 2011-10-05 中国航天标准化研究所 基于OpenCL的GPU加速三维时域有限差分电磁场仿真的方法
CN102353988A (zh) * 2011-07-08 2012-02-15 中国科学院地质与地球物理研究所 基于图形处理器计算起伏地表直接叠前逆时偏移的方法
CN102636808A (zh) * 2012-02-07 2012-08-15 浪潮(北京)电子信息产业有限公司 一种地震叠前时间偏移走时处理方法及装置
CN103675908A (zh) * 2012-09-21 2014-03-26 中国石油化工股份有限公司 一种海量数据图形处理器的波动方程逆时偏移成像方法
CN105005072A (zh) * 2015-06-02 2015-10-28 中国科学院地质与地球物理研究所 利用cuda的pml边界三维地震波传播模拟方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
地震叠前逆时偏移算法的CPU/GPU实施对策;李博 等;《地球物理学报》;20101231;第53卷(第12期);第2938-2943页 *
用GPU提速地震资料单程波有限差分叠前深度偏移;刘国峰 等;《APPLIED GEOPHYSICS》;20120331;第9卷(第1期);第41-48页 *

Also Published As

Publication number Publication date
CN105572730A (zh) 2016-05-11

Similar Documents

Publication Publication Date Title
Abdelkhalek et al. Fast seismic modeling and reverse time migration on a GPU cluster
CN103970960B (zh) 基于gpu并行加速的无网格伽辽金法结构拓扑优化方法
Bergmann et al. Algorithmic choices in WARP–A framework for continuous energy Monte Carlo neutron transport in general 3D geometries on GPUs
CN102565854B (zh) 一种海量数据gpu波动方程逆时偏移成像方法
CN105005072B (zh) 利用cuda的pml边界三维地震波传播模拟方法
CN103135132A (zh) Cpu/gpu协同并行计算的混合域全波形反演方法
CN108460195B (zh) 海啸数值计算模型基于gpu并行的快速执行方法
CN103675908A (zh) 一种海量数据图形处理器的波动方程逆时偏移成像方法
CN111638551A (zh) 地震初至波走时层析方法及装置
CN105572730B (zh) 三维复杂结构声波正演方法
Zhang et al. An fast simulation tool for fluid animation in VR application based on GPUs
Shankar et al. GRaM-X: A new GPU-accelerated dynamical spacetime GRMHD code for Exascale computing with the Einstein Toolkit
WO2013033651A1 (en) Full elastic wave equation for 3d data processing on gpgpu
CN109490948A (zh) 地震声学波动方程矢量并行计算方法
Reinhardt et al. Fully asynchronous SPH simulation
CN101963916A (zh) 编译处理方法及装置
CN108072895B (zh) 一种基于gpu的各向异性叠前逆时偏移优化方法
CN109001804A (zh) 一种基于三维地震数据确定有效应力的方法、装置及系统
Araújo et al. Boosting memory access locality of the Spectral Element Method with Hilbert space-filling curves
CN106548512B (zh) 网格模型数据的生成方法
Wei et al. Acceleration of a 2D unsteady Euler solver with GPU on nested Cartesian grid
CN107590353A (zh) 大气紊流场模拟方法和采用knl处理器的服务器的集群
Yalçın et al. GPU algorithms for diamond-based multiresolution terrain processing
Yao et al. A parallel volume rendering method for massive data
Fang et al. Optimizing complex spatially-variant coefficient stencils for seismic modeling on GPU

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20171114

Termination date: 20191215

CF01 Termination of patent right due to non-payment of annual fee