CN106662665B - 用于更快速的交错网格处理的重新排序的插值和卷积 - Google Patents
用于更快速的交错网格处理的重新排序的插值和卷积 Download PDFInfo
- Publication number
- CN106662665B CN106662665B CN201580040506.8A CN201580040506A CN106662665B CN 106662665 B CN106662665 B CN 106662665B CN 201580040506 A CN201580040506 A CN 201580040506A CN 106662665 B CN106662665 B CN 106662665B
- Authority
- CN
- China
- Prior art keywords
- stress
- space
- node
- model unit
- model
- 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
Links
- 238000012545 processing Methods 0.000 title claims description 20
- 230000008707 rearrangement Effects 0.000 title description 5
- 238000000034 method Methods 0.000 claims abstract description 37
- 238000009795 derivation Methods 0.000 claims abstract description 31
- 238000004088 simulation Methods 0.000 claims abstract description 5
- 238000003384 imaging method Methods 0.000 claims description 11
- 238000007906 compression Methods 0.000 claims description 7
- 230000006835 compression Effects 0.000 claims description 7
- 230000000638 stimulation Effects 0.000 claims description 2
- 230000008569 process Effects 0.000 description 14
- 230000006854 communication Effects 0.000 description 7
- 238000004891 communication Methods 0.000 description 6
- 230000008859 change Effects 0.000 description 5
- 230000006870 function Effects 0.000 description 4
- 230000004048 modification Effects 0.000 description 4
- 238000012986 modification Methods 0.000 description 4
- 239000002245 particle Substances 0.000 description 4
- 230000000644 propagated effect Effects 0.000 description 4
- 230000008901 benefit Effects 0.000 description 3
- 238000004422 calculation algorithm Methods 0.000 description 3
- 238000010586 diagram Methods 0.000 description 3
- 239000012530 fluid Substances 0.000 description 3
- 230000005012 migration Effects 0.000 description 3
- 238000013508 migration Methods 0.000 description 3
- 239000004215 Carbon black (E152) Substances 0.000 description 2
- 238000004364 calculation method Methods 0.000 description 2
- 230000005284 excitation Effects 0.000 description 2
- 238000001914 filtration Methods 0.000 description 2
- 229930195733 hydrocarbon Natural products 0.000 description 2
- 150000002430 hydrocarbons Chemical class 0.000 description 2
- 230000001902 propagating effect Effects 0.000 description 2
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 2
- 230000001133 acceleration Effects 0.000 description 1
- 238000003491 array Methods 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 244000309464 bull Species 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 238000005094 computer simulation Methods 0.000 description 1
- 238000013500 data storage Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000004069 differentiation Effects 0.000 description 1
- 238000006073 displacement reaction Methods 0.000 description 1
- 230000005672 electromagnetic field Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000009472 formulation Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000011221 initial treatment Methods 0.000 description 1
- 229910052500 inorganic mineral Inorganic materials 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 239000011707 mineral Substances 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 239000003129 oil well Substances 0.000 description 1
- 239000013307 optical fiber Substances 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 230000002085 persistent effect Effects 0.000 description 1
- 238000003908 quality control method Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 239000011435 rock Substances 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000005428 wave function Effects 0.000 description 1
Classifications
-
- 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
-
- 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
-
- 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/673—Finite-element; Finite-difference
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Remote Sensing (AREA)
- General Physics & Mathematics (AREA)
- Geology (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- General Life Sciences & Earth Sciences (AREA)
- Geophysics (AREA)
- Theoretical Computer Science (AREA)
- General Engineering & Computer Science (AREA)
- Geometry (AREA)
- Evolutionary Computation (AREA)
- Computer Hardware Design (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
一种改进的有限差分交错网格波传播系统和方法。一个说明性的基于计算机的波场模拟方法包括：将至少一个信号应用到形成模型空间的模型单元的网格，每个模型单元具有与应力节点相关联的应力值和与和应力节点交错的速度节点相关联的速度值；以及通过交替地更新所述应力值和所述速度值，将所述至少一个信号作为波传播到所述模型空间中，以获得与所述至少一个信号相关联的时间相关波场。应力值更新包括，对于每个模型单元：确定模型单元的速度值的空间导数；将所述空间导数插值到所述模型单元内的多个应力节点；以及对于所述模型单元内的每个应力节点，组合与所述应力节点相关联的空间导数以更新与所述应力节点相关联的至少一个应力值。
Description
相关申请的交叉引用
不适用。
关于联邦赞助的研究或开发的声明
不适用。
背景
技术领域
本发明一般涉及烃的地球物理勘探领域。更具体地，本发明涉及一种地震建模中更快速和更高效的网格处理的方法。
背景技术
科学家和工程师通常对勘探和工程项目采用地球物理勘测。地球物理勘测可以提供有关地下结构的信息，包括地层边界、岩石类型以及是否存在流体储层。这样的信息极大地有助于搜索水、地热储层和诸如烃和矿石的矿床。石油公司特别经常投资于广泛的地震和电磁勘测，以选择勘探油井的地点。
地球物理勘测可以在陆地上或在水中执行。如图1的示例勘测中所示，感兴趣区域104附近的能量源102生成传播到感兴趣区域中并从诸如地层边界的内部特征反射的波106。反射波108能量中的一些最终到达表面112上的接收器阵列110。记录系统114捕获所接收的信号以用于存储和处理。使用许多不同的源位置并且可选地使用不同的接收器位置来重复该过程。虽然存在用于将接收到的波信号转换为地下结构的图像的各种方法，但是最流行的这种技术采用有限差分波场建模，其是使用离散时间步长和波函数导数的快速近似在时间上向前或向后传播波的过程。利用足够小的时间步长和足够精确的近似，有限差分波场建模过程提供了高效、高保真的解决方案。然而，有限差分建模过程增加了相当大的计算需求。
已经进行了许多努力以满足有限差分波场建模过程的计算需求。在最大化运行时效率的尝试中，已经实现了许多组合的硬件-软件解决方案。除了利用由不同的硬件系统提供的不同的并行化策略之外，软件通常被定制为使用例如不同的离散化、波等式近似和传播策略来采用不同的波传播内核。在满足存储器、总线带宽和其它系统限制的同时以可接受的保真度实现更快速的处理(即，更短的运行时间)的尝试，已经导致在文献中公开的无数的内核变化。然而，计算负担仍然保持不合需要地高。
发明内容
因此，本文公开了有限差分交错网格波传播系统和方法，其具有通过在传播内核内的操作的重新排序而实现的减少的计算负担。一种说明性的计算机实现的波场模拟方法包括：将至少一个信号应用到包括多个模型单元的网格，所述网格形成模型空间，每个模型单元具有与多个应力节点相关联的多个应力值，以及与和所述应力节点交错的多个速度节点相关联的多个速度值；以及通过交替地更新所述应力值和所述速度值，将所述至少一个信号作为波传播到所述模型空间中，以获得与所述至少一个信号相关联的时间相关波场。应力值更新包括，对于每个模型单元：确定模型单元的速度值的空间导数；将所述空间导数插值到所述模型单元内的多个应力节点；以及对于所述模型单元内的每个应力节点，组合与所述应力节点相关联的空间导数以更新与所述应力节点相关联的至少一个应力值。
计算节约不需要限于波传播。一种说明性的基于计算机的方法提供用于将第一网格中的值应用于与第一网格交错的第二网格中的节点的组合卷积插值操作，所述方法包括：对于第二网格中的每一个单元：使用卷积滤波器的系数获得在所述第一网格中的对应单元周围的窗口中值的加权和；以及将所述加权和插值到所述第二网格的所述单元中的多个节点中的每一个。
所公开的非暂时性信息存储介质还提供用于地震成像方法的软件，所述地震成像方法包括：将至少一个地震源或接收器信号应用于包括多个模型单元的网格，所述网格形成模型空间，每个模型单元具有与多个应力节点相关联的多个应力值和与和应力节点交错的多个速度节点相关联的多个值；以及通过交替地更新所述应力值和所述速度值，将所述至少一个地震源或接收器信号作为波传播到所述模型空间中，以获得与所述至少一个地震源或接收器信号相关联的时间相关波场；重复所述应用和传播操作以获得多个时间相关波场；分析对应的源并接收一个或多个波场以获得匹配信息；堆叠所述匹配信息以对由所述模型空间表示的地下区域进行成像；以及显示所述图像。应力值更新包括，对于每个模型单元：确定模型单元的速度值的空间导数；将所述空间导数插值到所述模型单元内的多个应力节点；以及对于所述模型单元内的每个应力节点，组合与所述应力节点相关联的空间导数以更新与所述应力节点相关联的至少一个应力值。
附图说明
在附图中：
图1是地震勘测的示意图。
图2是具有有限差分波场模拟器的地震勘测系统的框图。
图3是说明性数字化接收信号迹线的曲线图。
图4是采用说明性建模系统的地质建模器的视图。
图5是适用于地球物理建模的说明性硬件平台的框图。
图6A是表示感兴趣的地下区域的数据空间。
图6B是具有交错网格值的数据空间的单元。
图6C是包括插值速度值的数据空间的单元。
图7是说明性地下成像方法的流程图。
图8是第一种类型的说明性偏移内核的流程图。
图9是第二种类型的说明性偏移内核的流程图。
然而，应当理解，在附图中给出的具体实施例及其具体实施方式不限制本公开。相反，它们为本领域技术人员提供了在随附权利要求的范围内辨别与给定实施例中的一个或多个一起包括的替代形式、等同物和修改的基础。
具体实施方式
图1的地震勘测可以使用图2的系统来执行。地震接收器阵列110将地震波转换为被放大和数字化的电信号。(说明性信号波形在图3中示出)。记录系统114经由总线202或其它通信路径收集数字化信号，并将数字化信号存储在非暂时性信息存储介质上以供以后处理。通常，每个数字化信号与参数204(诸如接收器位置和激发位置以及系统设计者认为有价值的其它信息)相关联。记录系统114可以执行一些初始处理以过滤和/或压缩数据，并且在至少一些情况下，执行质量控制。
记录系统114经由因特网或一些其它通信机制206将地震勘测数据提供给具有足够的计算资源的数据处理中心208，以用于成像过程。数据处理中心包括一个或多个计算机，其可以使用有限差分波场建模来执行偏移，从而将记录的地震信号转换成地下结构的三维地图或“图像”，其可以显示在监视器210上并存储在永久储存器中以供以后使用。
在数据处理中心208中，地震勘测数据的处理由诸如图4所示的地质建模器来协调。地质建模器可以采用工作站404的用户界面402来启动地震数据的处理并且查看和分析得到的地震图像。工作站404是诸如图5所示的地下建模系统的硬件平台的一部分。(其它硬件平台也是可用且合适的，包括超级计算机、大规模并行处理器网络、以及包括GPU或其它专用协处理器的异构架构)。
说明性硬件平台经由局域网(LAN)504将工作站404耦合到一个或多个多处理器计算机506。所述一个或多个多处理器计算机506又经由存储区域网络(SAN)508耦接到一个或多个共享存储单元510。使用个人工作站404，地质建模器能够将地震勘测数据加载到系统中，以配置和监测地震勘测数据的处理，并且从系统，任选地以体积图像的形式检索结果。
个人工作站404可以采取具有示出输入和结果数据的图形表示的显示器，以及具有使用户能够移动文件并执行处理软件的键盘的台式计算机的形式。LAN 504提供多处理器计算机506之间以及与个人工作站404的高速通信。LAN 504可以采取以太网网络的形式。
多处理器计算机(一个或多个)506提供并行处理能力以使得能够适当地促进输入数据的处理以导出结果数据。每个计算机506包括多个处理器512、分布式存储器514、内部总线516、SAN接口518和LAN接口520。每个处理器512对分配的任务进行操作，以解决整个问题的一部分，并且对整个结果的至少一部分做贡献。与每个处理器512相关联的是分布式存储器模块514，其存储应用软件和工作数据集，以供处理器使用。内部总线516经由相应的接口518、520提供处理器间通信和到SAN或LAN网络的通信。不同的计算机506中的处理器之间的通信可以由LAN 504提供。
SAN 508提供对共享存储设备510的高速访问。SAN 508可以采取例如光纤通道(Fibrechannel)或无限带宽(Infiniband)网络的形式。共享存储单元510可以是大型的、独立的信息存储单元，其采用磁盘介质用于非易失性数据存储。为了提高数据访问速度和可靠性，共享存储单元510可以被配置为冗余磁盘阵列(“RAID”)。
无论硬件平台采取所示形式还是另一合适形式，所有这些平台共有的问题是与处理器可以执行计算的速率相对的相对低的数据检索率。也就是说，设计用于处理大量数据的大多数计算机具有限制处理可以发生的速率的存储器瓶颈。因此，尽管通过增加可用存储器可以减少与许多可并行化处理相关联的计算负担，但是这样的解决方案可能无法任何更快速地完成处理。通常，更好的方法是减少所需的存储器访问量，即使这样的策略增加了由处理器承载的计算负担。
为了使硬件平台能够模拟波场并构建图像，软件将感兴趣区域建模为被划分为单元602的网格的二维或三维空间。图6A示出了3D空间的示例。每个单元具有通常与单元内的单个点相关联的代表性的一组值。然而，如文献中所公开的(参见例如R.W.Graves，“Simulating Seismic Wave Propagation in 3D Elastic Media Using Staggered-GridFinite Diffferences”，Bull.Seismological Soc.Am.，v86(4)pp1091-1106，1996年8月)，使用交错网格产生某些效率优点，至少在各向同性和正交各向异性介质中。在深入研究交错网格的细节之前，对成像过程的概述是有帮助的。
图7是可由数据处理中心中的计算机执行的说明性地震成像方法的流程图。在框702中，计算机获得勘测数据，包括数字化信号和相关联的源和接收器位置信息。还可以在该块中获得地下的初始模型，例如，以作为位置的函数指定估计密度和弹性。通常采用均匀密度模型作为起点并在成像过程期间逐步细化。
在框704中，计算机模拟通过将源信号(或合适的替代)从源位置偏移到模型空间中而由源生成的波场的演变(两个相关波参数(诸如力、压力或应力；以及位移、速度或应变)的空间分布)。也就是说，源生成根据波动方程从源位置向外传播的应力和粒子速度场。有限差分模型模拟该传播以作为时间的函数确定源波场(即，模型空间中每个点处的应力和粒子速度)。下面进一步描述传播内核。
在框706中，将相似的操作应用于接收信号。认识到记录的信号表示已经到达接收器位置的波，则相应的粒子速度或应力根据波动方程在时间上向后传播到模型空间中。可以使用与之前相同的有限差分传播内核，即使具有反转的时间索引，以作为时间的函数确定接收波场。(随着时间的反转，接收波形像源波形一样被处理。)
因为地球中的反射体在勘测过程期间将源波场转换为接收波场，所以这两个波场将在反射点处匹配。因此，在框708中，作为时间和位置的函数来分析源和接收波场，以识别匹配，并因此指示反射点的那些区域。对于每次地震激发重复该过程，并且将匹配信息相加(“堆叠”)以对整个感兴趣区域进行成像。
在框710中，分析图像以识别估计的地层参数(密度和/或弹性)中的不准确性。由于这种不准确性在图像数据中产生某些模式，因此它们使得初始估计可以被改进。可以迭代地重复框704-710以逐步精化和改进模型参数(以及得到的图像)。一旦模型参数收敛，则计算机就可以在框712中显示图像。然后可以基于所显示的图像采取动作，包括转向决定、着陆决定和完成决定。
特别要注意的是，偏移过程可能需要在达到收敛之前重复多次，并且偏移过程的每次迭代都需要针对每次激发，针对每个时间步长，传播等式在模型空间中的每个点处的解。因此，有限差分传播内核可以重复数千亿次。虽然该过程包括可以支持在并行处理器上的实现的显著程度的并行性，但是计算负担仍然很高。即使内核计算效率的小幅增益也可以大大节省执行时间。在偏移过程中的这种效率增益的益处不限于地震成像。例如，也可以采用偏移来模拟针对宽范围的物理(声波、电磁波、流体动力学)和应用环境的波场。
如前所述，通过使用交错网格实现了显著的效率增益。在图6B中示出了说明性的交错网格，其中来自图6A的3D空间的单个单元602具有与角节点相关联的正(压缩)应力值σxx、σyy、σzz，与相应面节点相关联的剪切应力值σxy、σxz、σyz以及与边缘节点相关联的速度值vx、vy、vz。边缘节点沿着一个轴从角节点偏移(“交错”)半个单元间距，并且面节点沿着两个轴从角节点偏移半个单元间距。在各向同性或甚至正交各向异性介质中，将速度和应力值相互关连的有限差分算子自然地以它们用于修改的节点值为中心，因此不需要额外的插值。
然而，如在J.W.Rector，G.M.Hoversten，K.T.Nihei，“Seismic Modeling EnginesConsortium”，Lawrence Berkley National Laboratory，2003年12月3日中所述，额外的相互依赖性存在于完全广义的各向异性介质中。因此，必须修改传播内核以做两件事之一。内核可以采用全模板(诸如在旋转交错网格模型中所使用的)，其需要四倍的存储器，或者内核可以根据需要使用插值来提供对中心。具体地，如图6C中加下划线的值所示，速度值可以被插值到单元602的主体中心节点和其它边缘节点，以支持应力值的更新。虽然这种插值增加了额外的计算负担，但是全模板方法是优选的，因为它避免了加剧存储器带宽瓶颈并且它避免了某些网漂移不稳定性。
使用张量符号，粒子速度vi和应力σij的动量守恒和本构等式为：
其中ρ是密度，i、j、k、l是在空间轴X、Y、Z范围内的索引，xi是沿给定轴的距离，fi是沿给定轴应用的力(例如，作为源信号的一部分)，以及cijkl是将应力与应变相关的弹性常数。依照约定，隐式求和发生在仅出现在等式右侧的那些索引上，即，等式(1a)的右侧具有从1到3(轴X到Z)范围内的遍及j的隐式求和，以及等式(1b)的右侧具有遍及k的隐式求和以及遍及l的隐式求和，每个范围从1到3。
等式(1b)当减小到有限差分形式时，需要插值以将差分场集中于交错网格中的应力值节点上。以σxy的等式的有限差分形式为例，我们有：
σ′xy＝σxy+cxyxxDx(vx)+cxyyyDy(vy)+cxyzzDz(vz)+cxyxy(Dx(vy)+Dy(vx))+
cxyxz(Dx(vz)+Dz(vx))+cxyyz(Dy(vz)+Dz(vy))， (1c)
其中Di()是沿给定轴的导数的有限差分近似。为了有限差分适当地集中于针对σxy的面节点上，速度场分量中的每一个都必须被插值。
虽然插值可以是计算上廉价的操作，但是三个速度场分量中的每一个被插值以完成四个不同节点处的信息，需要在单元内进行九次插值。然而，每个有限差分算子对来自多个单元的值进行操作。单元的数量(以及因此所需的插值的数量)随着所选择的有限差分算子的阶而变化。N阶的有限差分算子对来自N+1个单元的值进行运算。因此，用于更新单元的四个应力节点中的每一个的九个有限差分需要在3N+1个单元中执行9次插值，以进行总共27N+9次插值。对于3阶有限差分，需要90次插值来更新一个单元。更高的阶可以用于增加的保真度。
虽然存在不同程度的用于重复使用被插值的值的机会，但需要额外的存储器来存储这些值，直到它们可被重新使用。(在极端情况下，可以在该折衷通常不被认为可接受之前跨整个体积执行插值。因此，仍继续执行90次插值以更新每个单元。
考虑到这一点，图8和9是示出可用于实现偏移块704、706(图7)的对传播内核的不同方法的流程图。偏移块704、706通过将问题的各块(例如，模型空间的不同部分)分配给不同的处理器来并行化偏移操作。在块802开始，每个处理器接收其任务(一个或多个)，其可以表示模型空间的一部分，并且在本地存储器中对数据进行分级以用于到处理器高速缓存和来自处理器高速缓存的加速分页。四个索引可以用于跟踪通过传播过程的进展，这些索引表示三个空间维度和时间维度。在块804中初始化这些索引。
在块806中将第一数据块分页到高速缓存中。该示例假设每个数据块表示对应于给定y索引值的数据体的x-z切片。当传播内核遍历(iterate through)给定的x-z切片时，执行块808-824的集合以更新应力值。块806和826-830使传播内核遍历数据体的切片中的每一个，以完成针对给定时间步长的应力值更新。块831-834使传播内核遍历随后的时间步骤。
单元集合840表示对给定单元执行以更新与该单元相关联的应力值的操作。在块808中，处理器对沿着当前空间导数轴的N+1个单元中的每一个中的交错速度节点(见图6C)执行速度场分量的双线性插值。在块810中，针对当前单元中的应力值节点中的每一个处的速度场分量中的每一个确定沿着该轴的有限差分。在块812中，处理器检查是否存在更多的空间导数轴，并且如果存在，则在块814中选择下一个轴，并且处理器返回到块808。轴全部在当前单元中截断，因此，不需要重复当前单元的值的插值-仅需要插值沿所选择轴的附加的N个单元的值。
一旦已经确定当前单元中的每个应力节点处的空间导数的完整集合，则在块816中，处理器将它们与相应的弹性常数组合以更新当前单元的应力值。见等式(1b)、(1c)。
在块818中，处理器确定x索引是否已经达到其范围的结束，并且如果不是，则在块820中增加该值，并且控制返回到块808。否则，在块822中，处理器确定z索引是否已经到达其范围的结束，并且如果不是，则在块824中增加该值，重置x索引，并且控制返回到块808。否则，当前x-z切片已经完成。在块826中，存储更新的应力值。在块828中，处理器确定y索引是否已达到其范围的结束，并且如果不是，则在块830中增加该值，重置x和z索引，并且控制返回到块806以加载下一数据块。否则，已经针对完整的数据体更新了应力值，并且在块831中，处理器遍历数据体以根据等式(1a)更新速度场。一旦已经完成，数据体的完整的时间步骤就已经完成。在块832中，处理器确定时间索引是否已经到达其范围的结束，并且如果不是，则在块834中增加该值，重置x、y、z索引，并且控制返回到块806。
与图8的传播内核相反，图9的改进传播内核重新排序插值和导数运算。因此，在针对每个单元执行的操作集合(框850)内，传播内核通过在块852中沿给定轴找到速度分量的空间导数(有限差分)开始。对于x轴，块852产生导数Dx(vx)、Dx(vy)、Dx(vz)。由于在该步骤中不存在插值，因此这些导数的中心将变化。导数Dx(vx)、Dy(vy)、Dz(vz)将以角节点为中心，而三个面节点将各自具有两个相关联的导数。例如，导数Dx(vy)和Dy(vx)将以与σxy相关联的面节点为中心。由于交错节点的布置，对于该步骤不需要插值。然后，在块854中，处理器仅在当前单元内执行双线性插值，以确定所有应力节点处的导数值。由于存在4个应力节点并且在块852中计算的每个导数已经以它们中的一个为中心，因此每个导数需要被插值到三个其它节点。在考虑了所有空间轴之后，处理器仅需执行总共9次双线性插值，而不考虑有限差分阶N。因此，此重新排序导致节省27N次插值。
我们注意到，有限差分计算可以被认为是特定类型的滤波或其它卷积运算。现在上面已经强调的益处对于对所插值的值执行的其它滤波或卷积运算是类似可实现的。也就是说，通过互换插值和卷积运算的顺序来实现显著的计算节省。
一旦完全理解了上述公开内容，许多其它变化和修改对于本领域技术人员将变得显而易见。例如，硬件和软件实现的细节经受宽范围的变化、替代实施例和不同形式的优化或近似。此外，尽管前述公开内容是在地震成像系统的上下文中提供的，但是可以在包括声场、电磁场和流体流动力学的其它上下文中模拟波传播。旨在将所附权利要求解释为包括所有这样的变化和修改。
Claims (23)
1.一种计算机实现的波场模拟方法，包括：
将至少一个信号应用到包括多个模型单元的网格，所述网格形成模型空间，每个模型单元具有与多个应力节点相关联的多个应力值，以及与和所述应力节点交错的多个速度节点相关联的多个速度值，其中所述至少一个信号表示由地震传感器阵列接收的地震激发或地震道；
通过交替地更新所述应力值和所述速度值，将所述至少一个信号作为波传播到所述模型空间中，以获得与所述信号相关联的时间相关波场，其中所述更新应力值包括，对于每个模型单元：
确定所述模型单元的速度值的空间导数；
将所述空间导数插值到所述模型单元内的多个应力节点；
对于所述模型单元内的每个应力节点，组合与所述应力节点相关联的空间导数以更新与所述应力节点相关联的至少一个应力值；
重复所述应用操作和所述传播操作以获得多个时间相关波场；
分析对应的源并接收一个或多个波场以获得匹配信息；
堆叠所述匹配信息以生成由所述模型空间表示的地下区域的图像；以及
使用生成的图像采取动作。
2.根据权利要求1所述的方法，其中所述组合包括使用一个或多个弹性常数作为系数确定所述空间导数的加权和。
3.根据权利要求1所述的方法，还包括将所述时间相关波场存储在非暂时性信息存储介质上。
4.根据权利要求1所述的方法，其中每个模型单元的速度值包括在三个空间方向中的每一个方向上的一个或多个速度分量。
5.根据权利要求4所述的方法，其中每个模型单元的空间导数包括所述单元的每个速度分量在所述三个空间方向的每一个方向上的一个或多个空间导数。
6.根据权利要求5所述的方法，其中所述空间导数各自使用N阶的有限差分算子来确定，其中N大于或等于3。
7.根据权利要求5所述的方法，其中所述应力值包括三个空间方向中的每一个方向上的压缩应力值和三个剪切应力值，所述剪切应力值与从压缩应力节点偏移的相应剪切节点相关联，并且其中所述空间导数中的每一个被确定或插值到所述压缩应力节点和所述相应剪切节点中的每一个。
8.根据权利要求7所述的方法，其中所述插值包括双线性插值。
9.根据权利要求1所述的方法，其中所述插值包括线性插值。
10.根据权利要求6所述的方法，其中所述有限差分算子包括卷积滤波器。
11.根据权利要求1所述的方法，其中所述插值是有限差分波传播内核的一部分。
12.一种非暂时性信息存储介质，当以与处理单元可操作关系放置时，配置所述处理单元实现地震成像方法，所述地震成像方法包括：
将至少一个地震源或接收器信号应用于包括多个模型单元的网格，所述网格形成模型空间，每个模型单元具有与多个应力节点相关联的多个应力值，以及与和所述应力节点交错的多个速度节点相关联的多个速度值；以及
通过交替地更新所述应力值和所述速度值，将所述至少一个地震源或接收器信号作为波传播到所述模型空间中，以获得与所述至少一个地震源或接收器信号相关联的时间相关波场，其中所述更新所述应力值包括，对于每个模型单元：
确定所述模型单元的速度值的空间导数；
将所述空间导数插值到所述模型单元内的多个应力节点；
对于所述模型单元内的每个应力节点，组合与所述应力节点相关联的空间导数以更新与所述应力节点相关联的至少一个应力值；
重复所述应用操作和所述传播操作以获得多个时间相关波场；
分析对应的源并接收一个或多个波场以获得匹配信息；
堆叠所述匹配信息以生成由所述模型空间表示的地下区域的图像；以及
使用生成的图像采取动作。
13.根据权利要求12所述的介质，其中所述组合包括使用一个或多个弹性常数作为系数确定所述空间导数的加权和。
14.根据权利要求12所述的介质，其中所述插值包括双线性插值。
15.根据权利要求12所述的介质，其中每个模型单元的速度值包括在三个空间方向中的每一个方向上的一个或多个速度分量。
16.一种波场模拟系统，包括：
处理器；以及
可操作地连接到所述处理器的存储器，所述存储器存储指令，所述指令当被所述处理器执行时，使得所述系统：
将至少一个地震源或接收器信号应用于包括多个模型单元的网格，所述网格形成模型空间，每个模型单元具有与多个应力节点相关联的多个应力值，以及与和所述应力节点交错的多个速度节点相关联的多个速度值；以及
通过交替地更新所述应力值和所述速度值，将所述至少一个地震源或接收器信号作为波传播到所述模型空间中，以获得与所述至少一个地震源或接收器信号相关联的时间相关波场，其中所述更新所述应力值包括，对于每个模型单元：
确定所述模型单元的速度值的空间导数；
将所述空间导数插值到所述模型单元内的多个应力节点；
对于所述模型单元内的每个应力节点，组合与所述应力节点相关联的空间导数以更新与所述应力节点相关联的至少一个应力值；
重复所述应用操作和所述传播操作以获得多个时间相关波场；
分析对应的源并接收一个或多个波场以获得匹配信息；
堆叠所述匹配信息以生成由所述模型空间表示的地下区域的图像；以及
使用生成的图像采取动作。
17.根据权利要求16所述的系统，其中所述组合包括使用一个或多个弹性常数作为系数确定所述空间导数的加权和。
18.根据权利要求16所述的系统，其中所述插值包括双线性插值。
19.根据权利要求16所述的系统，其中每个模型单元的速度值包括在三个空间方向中的每一个方向上的一个或多个速度分量。
20.根据权利要求19所述的系统，其中对于每个模型单元的每个速度分量，该模型单元的空间导数包括在三个空间方向中的每一个方向上的一个或多个空间导数。
21.根据权利要求20所述的系统，其中所述空间导数各自使用N阶的有限差分算子来确定，其中N大于或等于3。
22.根据权利要求21所述的系统，其中所述应力值包括三个空间方向中的每一个方向上的压缩应力值和三个剪切应力值，所述剪切应力值与从压缩应力节点偏移的相应剪切节点相关联，并且其中所述空间导数中的每一个被确定或插值到所述压缩应力节点和所述相应剪切节点中的每一个。
23.根据权利要求16所述的系统，其中时间相关波场被存储在非暂时性信息存储介质上。
Applications Claiming Priority (3)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| US14/447,353 US9928315B2 (en) | 2014-07-30 | 2014-07-30 | Re-ordered interpolation and convolution for faster staggered-grid processing |
| US14/447,353 | 2014-07-30 | ||
| PCT/US2015/018512 WO2016018464A1 (en) | 2014-07-30 | 2015-03-03 | Re-ordered interpolation and convolution for faster staggered-grid processing |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN106662665A CN106662665A (zh) | 2017-05-10 |
| CN106662665B true CN106662665B (zh) | 2019-05-28 |
Family
ID=52808112
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201580040506.8A Expired - Fee Related CN106662665B (zh) | 2014-07-30 | 2015-03-03 | 用于更快速的交错网格处理的重新排序的插值和卷积 |
Country Status (6)
| Country | Link |
|---|---|
| US (1) | US9928315B2 (zh) |
| EP (1) | EP3175263A1 (zh) |
| CN (1) | CN106662665B (zh) |
| AU (1) | AU2015297012A1 (zh) |
| CA (1) | CA2952935A1 (zh) |
| WO (1) | WO2016018464A1 (zh) |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US10901103B2 (en) | 2018-03-20 | 2021-01-26 | Chevron U.S.A. Inc. | Determining anisotropy for a build section of a wellbore |
| CN110187379A (zh) * | 2019-05-22 | 2019-08-30 | 中铁二院工程集团有限责任公司 | 一种基于tsp法隧道超前地质预报效果的测试方法 |
| US11281825B2 (en) * | 2020-06-30 | 2022-03-22 | China Petroleum & Chemical Corporation | Computer-implemented method for high speed multi-source loading and retrieval of wavefields employing finite difference models |
| CN112329311A (zh) * | 2020-11-10 | 2021-02-05 | 中国石油大学(华东) | 地震波传播有限差分模拟方法、装置及计算机存储介质 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US4611312A (en) * | 1983-02-09 | 1986-09-09 | Chevron Research Company | Method of seismic collection utilizing multicomponent receivers |
| CN103412328A (zh) * | 2013-08-01 | 2013-11-27 | 中国石油天然气集团公司 | 基于交错网格有限差分算法的波数域保幅波场分离方法 |
| CN103777238A (zh) * | 2012-10-17 | 2014-05-07 | 中国石油化工股份有限公司 | 一种纯纵波各向异性波场模拟方法 |
Family Cites Families (13)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6295505B1 (en) | 1995-01-10 | 2001-09-25 | Schlumberger Technology Corporation | Method of filter generation for seismic migration using Remez algorithm |
| WO2001029581A1 (en) | 1999-10-21 | 2001-04-26 | Core Laboratories, Inc. | Transfer function method of seismic signal processing and exploration |
| US6687617B2 (en) | 2001-06-28 | 2004-02-03 | Pgs America, Inc. | Method and system for migration of seismic data |
| US6778909B1 (en) | 2002-10-05 | 2004-08-17 | 3Dgeo Development, Inc. | Seismic data processing systems and methods |
| US20040122596A1 (en) | 2002-12-19 | 2004-06-24 | Core Laboratories, Inc. | Method for high frequency restoration of seismic data |
| US7039525B2 (en) | 2003-09-23 | 2006-05-02 | Pgs Geophysical As | Method for seismic migration using explicit depth extrapolation operators with dynamically variable operator length |
| US20050203982A1 (en) | 2004-02-13 | 2005-09-15 | Joseph Kolibal | Method and apparatus for approximating, deconvolving and interpolating data using berstein functions |
| RU2339056C2 (ru) | 2004-04-07 | 2008-11-20 | Вестернджеко Сайзмик Холдингз Лимитед | Обобщенное трехмерное прогнозирование кратных волн от поверхности |
| US7035737B2 (en) | 2004-06-04 | 2006-04-25 | Pgs Americas, Inc. | Method for seismic wavefield extrapolation |
| US20070271040A1 (en) | 2006-05-22 | 2007-11-22 | Jiaxiang Ren | Method for seismic migration in anisotropic media with constrained explicit operators |
| GB0905261D0 (en) | 2009-03-27 | 2009-05-13 | Geco Technology Bv | Processing seismic data |
| GB2478574B (en) | 2010-03-11 | 2012-10-03 | Geco Technology Bv | Processing geophysical data |
| CN103119471A (zh) | 2010-09-20 | 2013-05-22 | 雪佛龙美国公司 | 用于生成地下构造的图像的系统和方法 |
-
2014
- 2014-07-30 US US14/447,353 patent/US9928315B2/en active Active
-
2015
- 2015-03-03 CN CN201580040506.8A patent/CN106662665B/zh not_active Expired - Fee Related
- 2015-03-03 CA CA2952935A patent/CA2952935A1/en not_active Abandoned
- 2015-03-03 EP EP15714054.2A patent/EP3175263A1/en not_active Withdrawn
- 2015-03-03 WO PCT/US2015/018512 patent/WO2016018464A1/en active Application Filing
- 2015-03-03 AU AU2015297012A patent/AU2015297012A1/en not_active Abandoned
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US4611312A (en) * | 1983-02-09 | 1986-09-09 | Chevron Research Company | Method of seismic collection utilizing multicomponent receivers |
| CN103777238A (zh) * | 2012-10-17 | 2014-05-07 | 中国石油化工股份有限公司 | 一种纯纵波各向异性波场模拟方法 |
| CN103412328A (zh) * | 2013-08-01 | 2013-11-27 | 中国石油天然气集团公司 | 基于交错网格有限差分算法的波数域保幅波场分离方法 |
Non-Patent Citations (2)
| Title |
|---|
| Accuracy of Staggered 3-D Finite-Difference Grids for Anisotropic Wave Propagation;Heiner igel 等;《SEG technical program expanded abstracts 1992》;19921231;第1244-1246页 |
| Anisotropic finite-difference algorithm for modeling elastic wave propagation in fractured coalbeds;Zhenglin Pei 等;《Geophysics》;20120229;第77卷(第1期);第C13–C26页 |
Also Published As
| Publication number | Publication date |
|---|---|
| WO2016018464A1 (en) | 2016-02-04 |
| AU2015297012A1 (en) | 2017-01-12 |
| US9928315B2 (en) | 2018-03-27 |
| CA2952935A1 (en) | 2016-02-04 |
| EP3175263A1 (en) | 2017-06-07 |
| US20160034612A1 (en) | 2016-02-04 |
| CN106662665A (zh) | 2017-05-10 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Martin et al. | Gravity inversion using wavelet-based compression on parallel hybrid CPU/GPU systems: application to southwest Ghana | |
| CN105467443A (zh) | 一种三维各向异性弹性波数值模拟方法及系统 | |
| EP1877834A1 (en) | 3d fast fault restoration | |
| Antonietti et al. | Numerical modeling of seismic waves by discontinuous spectral element methods | |
| CN106662665B (zh) | 用于更快速的交错网格处理的重新排序的插值和卷积 | |
| Minisini et al. | Local time stepping with the discontinuous Galerkin method for wave propagation in 3D heterogeneous media | |
| CN102692644A (zh) | 生成深度域成像道集的方法 | |
| Lee et al. | An optimized parallel LSQR algorithm for seismic tomography | |
| Homma et al. | A physics-based Monte Carlo earthquake disaster simulation accounting for uncertainty in building structure parameters | |
| Karavaev et al. | A technology of 3D elastic wave propagation simulation using hybrid supercomputers | |
| WO2016001697A1 (en) | Systems and methods for geologic surface reconstruction using implicit functions | |
| Alkhimenkov et al. | Resolving wave propagation in anisotropic poroelastic media using graphical processing units (GPUs) | |
| Muratov et al. | Grid-characteristic method on unstructured tetrahedral meshes | |
| Moradi et al. | Quantum computing in geophysics: Algorithms, computational costs, and future applications | |
| Seriani et al. | Numerical modeling of mechanical wave propagation | |
| Cruz-Atienza et al. | 3D finite-difference dynamic-rupture modeling along nonplanar faults | |
| CN105531602B (zh) | 利用多个加速处理部件(apc)实现时域有限差分模型的系统和方法 | |
| Cao et al. | Accelerating 2D and 3D frequency-domain seismic wave modeling through interpolating frequency-domain wavefields by deep learning | |
| US10454713B2 (en) | Domain decomposition using a multi-dimensional spacepartitioning tree | |
| CN107807392A (zh) | 一种自适应抗频散的分块时空双变逆时偏移方法 | |
| Zhang et al. | Efficient simulation of wave propagation with implicit finite difference schemes | |
| Suárez-Carreño et al. | Simulation of Wave Propagation Using Finite Differences in Oil Exploration | |
| Smith et al. | Elastic wave propagation in variable media using a discontinuous Galerkin method | |
| Antonietti et al. | Simulation of 3D elasto-acoustic wave propagation based on a Discontinuous Galerkin Spectral Element method | |
| Wu et al. | Parallel earthquake simulations on large-scale multicore supercomputers |
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 | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20190528 |
|
| CF01 | Termination of patent right due to non-payment of annual fee |