CN107462925A - 一种三维孔隙弹性介质中快速波场模拟方法 - Google Patents

一种三维孔隙弹性介质中快速波场模拟方法 Download PDF

Info

Publication number
CN107462925A
CN107462925A CN201710643565.5A CN201710643565A CN107462925A CN 107462925 A CN107462925 A CN 107462925A CN 201710643565 A CN201710643565 A CN 201710643565A CN 107462925 A CN107462925 A CN 107462925A
Authority
CN
China
Prior art keywords
mrow
msub
difference
munderover
msup
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.)
Granted
Application number
CN201710643565.5A
Other languages
English (en)
Other versions
CN107462925B (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.)
Xian Jiaotong University
Original Assignee
Xian Jiaotong University
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 Xian Jiaotong University filed Critical Xian Jiaotong University
Priority to CN201710643565.5A priority Critical patent/CN107462925B/zh
Publication of CN107462925A publication Critical patent/CN107462925A/zh
Application granted granted Critical
Publication of CN107462925B publication Critical patent/CN107462925B/zh
Active 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
    • G01V1/30—Analysis
    • G01V1/301—Analysis for determining seismic cross-sections or geostructures
    • G01V1/302—Analysis for determining seismic cross-sections or geostructures in 3D data cubes
    • 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
    • G01V1/30—Analysis
    • G01V1/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
    • 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
    • G01V1/34—Displaying seismic recordings or visualisation of seismic data or attributes
    • G01V1/345—Visualisation of seismic data or attributes, e.g. in 3D cubes
    • 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/62—Physical property of subsurface
    • G01V2210/624—Reservoir parameters
    • G01V2210/6242—Elastic parameters, e.g. Young, Lamé or Poisson
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/70—Other details related to processing
    • G01V2210/74—Visualisation of seismic data

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)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明公开了一种三维孔隙弹性介质中快速波场模拟方法，该方法针对的地质模型为三维孔隙弹性介质，基于频散参数，通过一定的约束条件使其小于误差阈值，从而为不同的速度选取相应的差分阶数，实现自适应差分阶数有限差分方法，有效减少运行时间，提高方法的有效性。首先，针对三维孔隙弹性介质波动方程，基于平面波理论和交错网格有限差分方法，得到该波动方程对应的数值频散关系；其次基于数值频散关系得到频散参数，通过使其满足一定的约束条件，实现不同速度采用不同的差分阶数，即变差分阶数有限差分方法；再次，引入基于MPI和OpenMP的并行算法，进一步缩短运行时间。最后通过模型算例验证本发明既可以缩短运行时间，又可以保证计算精度。

Description

一种三维孔隙弹性介质中快速波场模拟方法
技术领域
本发明属于地震勘探技术领域，涉及一种模拟方法，尤其是一种三维孔隙弹性介质中快速波场模拟方法，其针对三维孔隙弹性介质中高效的正演模拟技术。
背景技术
我国的常规油气藏开采已接近尾声，开发的重点将转向非常规致密油气藏。地震波探测是进行油气探测最有效的方法之一。波场数值模拟技术是地震波探测的重要基础，也是认识储层地震波响应的关键工具，更是地震资料解释中参数反演和成像的前提条件。目前的数值模拟方法主要针对声学、弹性、粘弹性等简单的单相介质的波动方程，但是这些介质理论不能描述真实的地质特征，特别是不能准确描述非常规致密油气藏的储层特征。孔隙弹性介质理论是现有的可准确描述真实油气藏物性特征的理论之一。因此开展基于孔隙弹性介质波动方程的波场数值模拟方法研究，具有重要的理论意义和应用价值。
此外，由于油气储藏存在于三维空间，并且其横向和纵向都有非均匀性。因此必须开展三维数值模拟方法的研究。但是，三维实际地球物理模型尺寸庞大，并且相对于声学、弹性介质，孔隙弹性介质的波动方程更加复杂。因此三维孔隙弹性介质的数值模拟将耗费大量的运行时间，使得方法的运行效率大大降低。如何提高其运算效率是将三维孔隙弹性介质波场模拟应用于实际中首要解决的问题之一。
有限差分方法实现简单，计算效率高，是最常用的数值模拟方法之一。有限差分方法是将波动方程中时间和空间导数用差商的形式表示，从而实现了时间和空间的离散。有限差分方法有很多种，包括显式和隐式有限差分方法，时间域和频率域有限差分方法，交错网格和旋转交错网格有限差分方法等等。
为减少有限差分方法的运行时间，即提高有限差分方法的有效性，有很多出色的研究工作。Jensen采用变网格有限差分方法减少运行时间，同时可更好地刻画复杂边界。Tessmer在模型的不同区域，使用不同时间步长同样很可观地减少了计算时间。Wang在博士论文中，研究了变网格大小和变时间步长的有限差分方法，并将其应用到逆时偏移中，取得了可喜成果。Liu和Sen提出自适应的变差分阶数有限差分方法，根据模型中不同区域的速度来确定差分阶数。结果表明，该方法不仅能够减少运行时间，并且不会降低结果的准确性。
然而，在提高三维孔隙弹性介质波场模拟的有效性方面，还未有相关方法与策略。
发明内容
本发明的目的在于克服上述现有技术的缺点，提供一种三维孔隙弹性介质中快速波场模拟方法，该方法可减少正演模拟的运行时间，从而可应用于实际地球物理模型的正演模拟。
本发明的目的是通过以下技术方案来实现的：
这种三维孔隙弹性介质中快速波场模拟方法：首先，针对三维孔隙弹性介质波动方程，基于平面波理论和交错网格有限差分方法，得到三维孔隙弹性介质波动方程对应的数值频散关系；其次基于数值频散关系得到频散参数，通过使频散参数满足约束条件，实现不同速度采用不同的差分阶数，即变差分阶数有限差分方法；再次，引入基于MPI和OpenMP的并行算法，缩短运行时间。
进一步，以上三维孔隙弹性介质中快速波场模拟方法具体包括以下步骤：
1)三维孔隙弹性介质波动方程的数值频散关系推导
基于平面波理论及交错网格有限差分方法，得到三维孔隙弹性介质波动方程的数值频散关系：
其中，下标i＝1,2或3；当i＝1时表示慢P波，i＝2时表示S波，i＝3时表示快P波；k为波数，ω为角频率；θ是波传播方向与z轴夹角，φ是xoy平面与z轴夹角；2M为有限差分方法的阶数，τ为时间步长，h是空间网格大小，am(m＝1,2,...,M)是交错网格有限差分方法的差分系数；
2)频散参数的定义
数值频散由频散参数来衡量，表示为经有限差分方法离散后的速度与真实速度的比值：
其中：
当δ等于1时，没有数值频散；当δ大于或小于1时，会产生数值频散；
频散参数定义为经有限差分方法离散后的波传播一个离散网格的时间与实际传播时间之差：
当ε等于0时，没有数值频散；当ε大于或小于0时，会产生数值频散；
3)制定约束条件，实现自适应差分阶数有限差分方法
定义快P波、S波及慢P波频散参数的参数ξ：
在频率fmax内，通过约束参数ξ使其小于误差阈值η；
通过公式(6)，自适应地为不同速度选取相应的差分阶数，从而实现自适应差分阶数有限差分方法，减少运行时间且不降低精度；
4)在正演模拟过程中，引入并行算法
针对三维模型的有限差分方法，采用消息传递接口(MPI)及OpenMP，对算法进行并行化，全面缩短运行时间。
本发明具有以下有益效果：
本发明的三维孔隙弹性介质中快速波场模拟方法针对的地质模型为三维孔隙弹性介质，基于频散参数，通过一定的约束条件使其小于误差阈值，从而为不同的速度选取相应的差分阶数，实现自适应差分阶数有限差分方法，有效减少运行时间，提高方法的有效性。
进一步，本发明方便灵活，可提前为模型中不同速度计算相应的差分阶数，在正演模拟过程中不需要在边界处进行额外处理。
进一步，本发明可有效减少正演模拟的运行时间，并且不会降低方法的精度。
进一步，本发明引入并行算法，可进一步缩短运行时间。
附图说明
图1是频散参数log10|ε|随频率变化曲线，其中τ＝0.001s，h＝25m；
图2是不同速度、误差阈值的差分阶数，其中τ＝0.001s，h＝25m；
图3是并行处理时数据划分示意图；
图4是进程间通信示意图；
图5是固体速度z分量在0.6s的波场快照，由(a)固定差分阶数(M＝11)(b)变差分阶数(M＝5,11)有限差分方法计算的结果，及(c)二者的差剖面；
图6是三维推覆模型示意图；
图7是基于推覆模型的固体速度z分量的地震记录，由(a)固定差分阶数(M＝12)，(b)变差分阶数(M＝4，4，4，5，6，6，12)及(c)二者差剖面，其中接收器位于平面y＝10km和平面z＝400m的交叉线上；
图8是基于推覆模型的固体速度z分量的地震记录，由(a)固定差分阶数(M＝12)，(b)变差分阶数(M＝4，4，4，5，6，6，12)及(c)二者差剖面，其中接收器位于平面x＝10km和平面z＝400m的交叉线上；
图9是基于推覆模型的固体速度z分量的地震记录，由(a)固定差分阶数(M＝12)，(b)变差分阶数(M＝4，4，4，5，6，6，12)及(c)二者差剖面，其中接收器位于平面x＝10km和平面y＝10km的交叉线上。
具体实施方式
本发明针对的地质模型为三维孔隙弹性介质，基于该波动方程对应的数值频散关系，通过约束三种波的频散参数使其小于误差阈值，进而实现对不同速度采用不同的差分阶数，有效减少正演模拟的运行时间。首先基于平面波理论和交错网格有限差分方法，得到三维孔隙弹性介质波动方程的数值频散关系；其次基于数值频散关系定义频散参数，通过约束三种波的频散参数使其小于误差阈值，从而实现自适应差分阶数有限差分方法；最后，引入基于MPI和OpenMP的并行算法，进一步缩短三维孔隙弹性介质波场模拟的运行时间。
下面结合附图和实例对本发明进行详细的描述。
1)三维孔隙弹性介质波动方程的数值频散关系的推导
首先，三维孔隙弹性介质的波动方程可写为
其中，
Q＝[vx，vy，vz，wx，wy，wz，τxx，τyy，τzz，τxy，τxz，τyz，s]T (14)
(vx,vy,vz)是固体速度分量，(wx,wy,wz)是流体速度分量，(τxx,τxz,τxy,τyy,τyz,τzz)是固体应力，s与流体压力有关。
将上述时间空间域波动方程变换到频率波数域，可得该波动方程的频散关系：
这里，vS，vsP和vfP表示S波、慢P波及快P波的速度。
基于平面波理论和交错网格有限差分方法及公式(9)，得到数值频散关系：
其中，下标i＝1,2,3。当i＝1时表示慢P波，i＝2时表示S波，i＝3时表示快P波；k为波数，ω为角频率；θ是波传播方向与z轴夹角，φ是xoy平面与z轴夹角；2M为有限差分方法的阶数，τ为时间步长，h是空间网格大小，am(m＝1,2,...,M)是交错网格有限差分方法的差分系数。
2)频散参数的定义
数值频散可用频散参数衡量，基于公式(10)，频散参数可表示为有限差分方法离散后的速度与真实速度的比值
其中，
当δ等于1时，没有数值频散；当δ大于或小于1时，会产生数值频散。
频散参数也可定义为经有限差分方法离散后的波传播一个离散网格的时间与实际传播时间之差
当ε等于0时，没有数值频散；当ε大于或小于0时，会产生数值频散。
3)制定约束条件，实现自适应差分阶数有限差分方法
为了衡量三维孔隙弹性介质的数值频散，定义包含了快P波、S波及慢P波频散参数的参数ξ
在一定频率fmax内，通过约束参数ξ使其小于误差阈值η
通过公式(15)，可以自适应地为不同速度选取相应的差分阶数，从而实现自适应差分阶数有限差分方法，有效减少运行时间但不会降低精度。图1给出固定差分阶数及自适应差分阶数有限差分方法计算的频散曲线。图1(a)中不同速度采用相同的差分阶数；而(b)中根据公式(15)为不同速度采用不同的差分阶数。对比两图可看出，为高速度采用低阶有限差分方法，并不会降低方法的精度，说明了本发明的可行性。
图2是不同速度、不同误差阈值对应的差分阶数，可看出，速度和误差阈值越小，差分阶数越大。
4)在正演模拟过程中，引入并行算法
在并行计算中，计算区域被分成更小的子域，每个小子域由相应的进程进行处理；进程之间的通信由MPI完成。有限差分方法的实现过程中，仅仅需要区域边界周围的几个点，因此有限差分方法很容易进行并行化处理。
一般而言，对于三维并行有限差分方法，计算区域应当在三个坐标轴上同时划分。并且应均匀划分，以平衡进程之间的工作量及通信量，如图3所示。
图4为进程间通信的示意图，每个进程在自己的子区域处理问题，并在每个时刻点，相邻的进程进行通信并交换信息。为此，需在每个进程的边界处开辟一定大小的空间，以存储从相邻进程接收的信息。该空间的厚度等于差分阶数的一半，即M。
此外，在程序实现时，空间循环部分可采用OpenMP，对线程进行优化，进一步提高运算效率。
数值结果
层状介质
该模型大小为[0,5000m]×[0,5000m]×[0,5000m]，分界面位于深度1250米处，模型参数如表1所示。数值模拟时，时间步长为0.001s，网格大小为25m。震源为主频10Hz，时延0.1s的Ricker子波，位于模型中心。
表1层状介质的参数
基于公式(21)为该模型中不同速度计算对应的差分阶数，层1的差分阶数为10，层2的差分阶数为22。22阶空间精度的有限差分方法模拟结果作为参考结果。固体速度z分量的波场快照如图5所示，由于模型是横向均匀的，因此只显示过震源的yoz平面的切片。从图中可看出，变差分阶数和固定差分阶数有限差分方法的数值结果无明显差别，二者的误差仅在0.7％。
方法的有效性由运行时间来测试。表2列出了不同层厚情况下，两种有限差分方法的运行时间。可看出，本发明的运行效率与模型特征非常相关。高速区越大、低速区越小，有效性越高，即越省时。
表2不同层厚的计算时间
简化的推覆模型
下面通过模拟一个稍复杂模型来测试本发明的有效性。该模型为包含断层和推覆层的SEG推覆模型，如图6所示。模型大小为[0,20.025km]×[0,20.025km]×[0,4.675km]，时间步长为0.001s，空间网格大小为25m。震源为主频10Hz，时延0.1s的Ricker子波，位于(10km,10km,400m)。根据公式(21)为本模型七组数据计算相应的差分阶数，如表3所示，其中η＝1e-3。可以看出，对于该模型，最大的差分阶数为24，最小的差分阶数为8。24作为固定差分阶数有限差分方法的差分阶数。
表3推覆模型中不同速度对应的差分阶数
两种有限差分方法计算的固体速度z分量的地震记录及其差剖面如图7-9所示，其中接收器分别位于(…,10km,400m)，(10km,…,400m)，(10km,10km,…)。从这三幅图可以看出，两种有限差分方法计算的结果并无明显差异，并且误差分别仅为0.34％，0.34％，0.28％。因此，变差分阶数有限差分方法并不会降低方法的精度。
表4位两种有限差分方法的运行时间，可以看出变差分阶数有限差分方法可以节省33.54％的时间，有效性得以大大提高。
表4两种有限差分方法的运行时间

Claims (2)

1.一种三维孔隙弹性介质中快速波场模拟方法，其特征在于，首先，针对三维孔隙弹性介质波动方程，基于平面波理论和交错网格有限差分方法，得到三维孔隙弹性介质波动方程对应的数值频散关系；其次基于数值频散关系得到频散参数，使频散参数满足约束条件，实现不同速度采用不同的差分阶数，即变差分阶数有限差分方法；再次，引入基于MPI和OpenMP的并行算法，缩短运行时间。
2.根据权利要求1所述的三维孔隙弹性介质中快速波场模拟方法，其特征在于，具体包括以下步骤：
1)三维孔隙弹性介质波动方程的数值频散关系推导
基于平面波理论及交错网格有限差分方法，得到三维孔隙弹性介质波动方程的数值频散关系：
<mrow> <mtable> <mtr> <mtd> <mrow> <msup> <mrow> <mo>(</mo> <munderover> <mi>&amp;Sigma;</mi> <mrow> <mi>m</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </munderover> <msub> <mi>a</mi> <mi>m</mi> </msub> <mi>sin</mi> <mo>(</mo> <mrow> <mrow> <mo>(</mo> <mrow> <mi>m</mi> <mo>-</mo> <mn>0.5</mn> </mrow> <mo>)</mo> </mrow> <msub> <mi>k</mi> <mi>i</mi> </msub> <mi>h</mi> <mi> </mi> <mi>sin</mi> <mi>&amp;theta;</mi> <mi>cos</mi> <mi>&amp;phi;</mi> </mrow> <mo>)</mo> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <munderover> <mi>&amp;Sigma;</mi> <mrow> <mi>m</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </munderover> <msub> <mi>a</mi> <mi>m</mi> </msub> <mi>s</mi> <mi>i</mi> <mi>n</mi> <mo>(</mo> <mrow> <mrow> <mo>(</mo> <mrow> <mi>m</mi> <mo>-</mo> <mn>0.5</mn> </mrow> <mo>)</mo> </mrow> <msub> <mi>k</mi> <mi>i</mi> </msub> <mi>h</mi> <mi> </mi> <mi>s</mi> <mi>i</mi> <mi>n</mi> <mi>&amp;theta;</mi> <mi>s</mi> <mi>i</mi> <mi>n</mi> <mi>&amp;phi;</mi> </mrow> <mo>)</mo> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>+</mo> <msup> <mrow> <mo>(</mo> <munderover> <mi>&amp;Sigma;</mi> <mrow> <mi>m</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </munderover> <msub> <mi>a</mi> <mi>m</mi> </msub> <mi>s</mi> <mi>i</mi> <mi>n</mi> <mo>(</mo> <mrow> <mrow> <mo>(</mo> <mrow> <mi>m</mi> <mo>-</mo> <mn>0.5</mn> </mrow> <mo>)</mo> </mrow> <msub> <mi>k</mi> <mi>i</mi> </msub> <mi>h</mi> <mi> </mi> <mi>c</mi> <mi>o</mi> <mi>s</mi> <mi>&amp;theta;</mi> </mrow> <mo>)</mo> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>=</mo> <msup> <mrow> <mo>(</mo> <msubsup> <mi>r</mi> <mi>i</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msubsup> <mi>s</mi> <mi>i</mi> <mi>n</mi> <mo>(</mo> <mrow> <mn>0.5</mn> <mi>&amp;omega;</mi> <mi>&amp;tau;</mi> </mrow> <mo>)</mo> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </mtd> </mtr> </mtable> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>
其中，下标i＝1,2或3；当i＝1时表示慢P波，i＝2时表示S波，i＝3时表示快P波；k为波数，ω为角频率；θ是波传播方向与z轴夹角，φ是xoy平面与z轴夹角；2M为有限差分方法的阶数，τ为时间步长，h是空间网格大小，am(m＝1,2,…,M)是交错网格有限差分方法的差分系数；
2)频散参数的定义
数值频散由频散参数来衡量，表示为经有限差分方法离散后的速度与真实速度的比值：
<mrow> <msub> <mi>&amp;delta;</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>&amp;theta;</mi> <mo>,</mo> <mi>&amp;phi;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <msub> <mi>v</mi> <mrow> <msub> <mi>FD</mi> <mi>i</mi> </msub> </mrow> </msub> <msub> <mi>v</mi> <mi>i</mi> </msub> </mfrac> <mo>=</mo> <mfrac> <mn>2</mn> <mrow> <msub> <mi>r</mi> <mi>i</mi> </msub> <msub> <mi>k</mi> <mi>i</mi> </msub> <mi>h</mi> </mrow> </mfrac> <mi>arc</mi> <mi> </mi> <mi>s</mi> <mi>i</mi> <mi>n</mi> <mrow> <mo>(</mo> <msub> <mi>r</mi> <mi>i</mi> </msub> <msqrt> <mi>q</mi> </msqrt> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow>
其中：
<mrow> <mtable> <mtr> <mtd> <mrow> <mi>q</mi> <mo>=</mo> <msup> <mrow> <mo>(</mo> <munderover> <mi>&amp;Sigma;</mi> <mrow> <mi>m</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </munderover> <msub> <mi>a</mi> <mi>m</mi> </msub> <mi>sin</mi> <mo>(</mo> <mrow> <mrow> <mo>(</mo> <mrow> <mi>m</mi> <mo>-</mo> <mn>0.5</mn> </mrow> <mo>)</mo> </mrow> <msub> <mi>k</mi> <mi>i</mi> </msub> <mi>h</mi> <mi> </mi> <mi>sin</mi> <mi>&amp;theta;</mi> <mi>cos</mi> <mi>&amp;phi;</mi> </mrow> <mo>)</mo> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <munderover> <mi>&amp;Sigma;</mi> <mrow> <mi>m</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </munderover> <msub> <mi>a</mi> <mi>m</mi> </msub> <mi>s</mi> <mi>i</mi> <mi>n</mi> <mo>(</mo> <mrow> <mrow> <mo>(</mo> <mrow> <mi>m</mi> <mo>-</mo> <mn>0.5</mn> </mrow> <mo>)</mo> </mrow> <msub> <mi>k</mi> <mi>i</mi> </msub> <mi>h</mi> <mi> </mi> <mi>s</mi> <mi>i</mi> <mi>n</mi> <mi>&amp;theta;</mi> <mi>s</mi> <mi>i</mi> <mi>n</mi> <mi>&amp;phi;</mi> </mrow> <mo>)</mo> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>+</mo> <msup> <mrow> <mo>(</mo> <munderover> <mi>&amp;Sigma;</mi> <mrow> <mi>m</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </munderover> <msub> <mi>a</mi> <mi>m</mi> </msub> <mi>s</mi> <mi>i</mi> <mi>n</mi> <mo>(</mo> <mrow> <mrow> <mo>(</mo> <mrow> <mi>m</mi> <mo>-</mo> <mn>0.5</mn> </mrow> <mo>)</mo> </mrow> <msub> <mi>k</mi> <mi>i</mi> </msub> <mi>h</mi> <mi> </mi> <mi>c</mi> <mi>o</mi> <mi>s</mi> <mi>&amp;theta;</mi> </mrow> <mo>)</mo> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </mtd> </mtr> </mtable> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow>
当δ等于1时，没有数值频散；当δ大于或小于1时，会产生数值频散；
频散参数定义为经有限差分方法离散后的波传播一个离散网格的时间与实际传播时间之差：
<mrow> <msub> <mi>&amp;epsiv;</mi> <mi>i</mi> </msub> <mo>=</mo> <mfrac> <mi>h</mi> <msub> <mi>v</mi> <mrow> <msub> <mi>FD</mi> <mi>i</mi> </msub> </mrow> </msub> </mfrac> <mo>-</mo> <mfrac> <mi>h</mi> <msub> <mi>v</mi> <mi>i</mi> </msub> </mfrac> <mo>=</mo> <mfrac> <mi>h</mi> <msub> <mi>v</mi> <mi>i</mi> </msub> </mfrac> <mrow> <mo>(</mo> <msubsup> <mi>&amp;delta;</mi> <mi>i</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msubsup> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow>
当ε等于0时，没有数值频散；当ε大于或小于0时，会产生数值频散；
3)制定约束条件，实现自适应差分阶数有限差分方法
定义快P波、S波及慢P波频散参数的参数ξ：
<mrow> <mi>&amp;xi;</mi> <mo>=</mo> <mrow> <mo>(</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mn>3</mn> </munderover> <mo>|</mo> <mrow> <msub> <mi>&amp;epsiv;</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mrow> <msub> <mi>v</mi> <mi>i</mi> </msub> <mo>,</mo> <mi>M</mi> <mo>,</mo> <mi>f</mi> </mrow> <mo>)</mo> </mrow> </mrow> <mo>|</mo> <mo>)</mo> </mrow> <mo>/</mo> <mn>3</mn> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </mrow>
在频率fmax内，通过约束参数ξ使其小于误差阈值η；
<mrow> <msub> <mi>&amp;xi;</mi> <mrow> <mi>m</mi> <mi>a</mi> <mi>x</mi> </mrow> </msub> <mrow> <mo>(</mo> <msub> <mi>v</mi> <mrow> <mi>s</mi> <mi>P</mi> </mrow> </msub> <mo>,</mo> <msub> <mi>v</mi> <mi>S</mi> </msub> <mo>,</mo> <msub> <mi>v</mi> <mrow> <mi>f</mi> <mi>P</mi> </mrow> </msub> <mo>,</mo> <mi>M</mi> <mo>,</mo> <msub> <mi>f</mi> <mrow> <mi>m</mi> <mi>a</mi> <mi>x</mi> </mrow> </msub> <mo>)</mo> </mrow> <mo>=</mo> <munder> <mrow> <mi>m</mi> <mi>a</mi> <mi>x</mi> </mrow> <mrow> <mi>f</mi> <mo>&amp;Element;</mo> <mo>&amp;lsqb;</mo> <mn>0</mn> <mo>,</mo> <msub> <mi>f</mi> <mrow> <mi>m</mi> <mi>a</mi> <mi>x</mi> </mrow> </msub> <mo>&amp;rsqb;</mo> <mo>,</mo> <mi>&amp;theta;</mi> <mo>&amp;Element;</mo> <mo>&amp;lsqb;</mo> <mn>0</mn> <mo>,</mo> <mi>&amp;pi;</mi> <mo>&amp;rsqb;</mo> <mo>,</mo> <mi>&amp;phi;</mi> <mo>&amp;Element;</mo> <mo>&amp;lsqb;</mo> <mn>0</mn> <mo>,</mo> <mn>2</mn> <mi>&amp;pi;</mi> <mo>&amp;rsqb;</mo> </mrow> </munder> <mrow> <mo>(</mo> <mi>&amp;xi;</mi> <mo>)</mo> </mrow> <mo>&lt;</mo> <mi>&amp;eta;</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow>
通过公式(6)，自适应地为不同速度选取相应的差分阶数，从而实现自适应差分阶数有限差分方法，减少运行时间且不降低精度；
4)在正演模拟过程中，引入并行算法
针对三维模型的有限差分方法，采用消息传递接口(MPI)及OpenMP，对算法进行并行化，全面缩短运行时间。
CN201710643565.5A 2017-07-31 2017-07-31 一种三维孔隙弹性介质中快速波场模拟方法 Active CN107462925B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710643565.5A CN107462925B (zh) 2017-07-31 2017-07-31 一种三维孔隙弹性介质中快速波场模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710643565.5A CN107462925B (zh) 2017-07-31 2017-07-31 一种三维孔隙弹性介质中快速波场模拟方法

Publications (2)

Publication Number Publication Date
CN107462925A true CN107462925A (zh) 2017-12-12
CN107462925B CN107462925B (zh) 2018-10-30

Family

ID=60546969

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710643565.5A Active CN107462925B (zh) 2017-07-31 2017-07-31 一种三维孔隙弹性介质中快速波场模拟方法

Country Status (1)

Country Link
CN (1) CN107462925B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107908913A (zh) * 2017-12-22 2018-04-13 中国海洋大学 基于并行计算机的地球动力数模算法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103630933A (zh) * 2013-12-09 2014-03-12 中国石油天然气集团公司 基于非线性优化的时空域交错网格有限差分方法和装置
CN105911584A (zh) * 2015-09-25 2016-08-31 中国科学院地质与地球物理研究所 一种隐式交错网格有限差分弹性波数值模拟方法及装置
US20160334530A1 (en) * 2013-12-30 2016-11-17 Denis Evgenievich SYRESIN Method and system for processing acoustic waveforms
CN106353797A (zh) * 2015-07-17 2017-01-25 中国石油化工股份有限公司 一种高精度地震正演模拟方法
CN106646597A (zh) * 2016-12-14 2017-05-10 中国石油大学(北京) 基于弹簧网络模型的正演模拟方法及装置

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103630933A (zh) * 2013-12-09 2014-03-12 中国石油天然气集团公司 基于非线性优化的时空域交错网格有限差分方法和装置
US20160334530A1 (en) * 2013-12-30 2016-11-17 Denis Evgenievich SYRESIN Method and system for processing acoustic waveforms
CN106353797A (zh) * 2015-07-17 2017-01-25 中国石油化工股份有限公司 一种高精度地震正演模拟方法
CN105911584A (zh) * 2015-09-25 2016-08-31 中国科学院地质与地球物理研究所 一种隐式交错网格有限差分弹性波数值模拟方法及装置
CN106646597A (zh) * 2016-12-14 2017-05-10 中国石油大学(北京) 基于弹簧网络模型的正演模拟方法及装置

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
周建科: "波动方程有限元法数值模拟及井震标定研究", 《中国石油大学硕士学位论文》 *
周艳杰: "求解波动方程的低数值频散方法及数值模拟", 《清华大学博士学位论文》 *
张邦: "三维复杂介质弹性波场有限差分数值模拟方法研究", 《长安大学硕士学位论文》 *
蔡其新等: "有限差分数值模拟的最小频散算法及其应用", 《石油地球物理勘探》 *
高晗: "黏弹孔隙介质波场数值模拟方法研究", 《吉林大学硕士学位论文》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107908913A (zh) * 2017-12-22 2018-04-13 中国海洋大学 基于并行计算机的地球动力数模算法
CN107908913B (zh) * 2017-12-22 2020-12-11 中国海洋大学 基于并行计算机的地球动力数模方法

Also Published As

Publication number Publication date
CN107462925B (zh) 2018-10-30

Similar Documents

Publication Publication Date Title
Seriani et al. Spectral element method for acoustic wave simulation in heterogeneous media
CN103149585B (zh) 一种弹性偏移地震波场构建方法及装置
CN105044771B (zh) 基于有限差分法的三维tti双相介质地震波场数值模拟方法
Hægland et al. Improved streamlines and time-of-flight for streamline simulation on irregular grids
CN105549068A (zh) 一种三维各向异性微地震干涉逆时定位方法及系统
CN112327358B (zh) 一种粘滞性介质中声波地震数据正演模拟方法
Chen et al. A normalized wavefield separation cross-correlation imaging condition for reverse time migration based on Poynting vector
CN107561585A (zh) 一种多核多节点并行三维地震波场生成方法和系统
CN107678062A (zh) 双曲Radon域综合预测反褶积和反馈循环方法压制多次波模型构建方法
Ganis et al. Modeling fluid injection in fractures with a reservoir simulator coupled to a boundary element method
CN105911584B (zh) 一种隐式交错网格有限差分弹性波数值模拟方法及装置
Vamaraju et al. Enriched Galerkin finite element approximation for elastic wave propagation in fractured media
CN103926619A (zh) 一种三维vsp数据的逆时偏移方法
Sun et al. 2-D poroelastic wave modelling with a topographic free surface by the curvilinear grid finite-difference method
CN109946742A (zh) 一种TTI介质中纯qP波地震数据模拟方法
CN104965222A (zh) 三维纵波阻抗全波形反演方法及装置
CN108279437A (zh) 变密度声波方程时间高阶精度交错网格有限差分方法
CN105676280A (zh) 基于旋转交错网格的双相介质地质数据获取方法和装置
Wang et al. A high-efficiency spectral element method based on CFS-PML for GPR numerical simulation and reverse time migration
CN107340537A (zh) 一种p-sv转换波叠前逆时深度偏移的方法
CN107462925B (zh) 一种三维孔隙弹性介质中快速波场模拟方法
Wang et al. Acoustic reverse time migration and perfectly matched layer in boundary-conforming grids by elliptic method
Wang et al. A fourth order accuracy summation-by-parts finite difference scheme for acoustic reverse time migration in boundary-conforming grids
CN105842745A (zh) 重力异常的小波域优化位变滤波分离方法
CN105319594A (zh) 一种基于最小二乘参数反演的傅里叶域地震数据重构方法

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