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

摘要

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

权利要求

1.一种三维孔隙弹性介质中快速波场模拟方法，其特征在于，首先，针对三维孔隙弹性介质波动方程，基于平面波理论和交错网格有限差分方法，得到三维孔隙弹性介质波动方程对应的数值频散关系；其次基于数值频散关系得到频散参数，使频散参数满足约束条件，实现不同速度采用不同的差分阶数，即变差分阶数有限差分方法；再次，引入基于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)频散参数的定义

数值频散由频散参数来衡量，表示为经有限差分方法离散后的速度与真实速度的比值：其中：vi为波的真实速度， 为有限差分方法离散后波的速度；

当δ等于1时，没有数值频散；当δ大于或小于1时，会产生数值频散；

频散参数定义为经有限差分方法离散后的波传播一个离散网格的时间与实际传播时间之差：当ε等于0时，没有数值频散；当ε大于或小于0时，会产生数值频散；

3)制定约束条件，实现自适应差分阶数有限差分方法定义快P波、S波及慢P波频散参数的参数ξ：

在频率fmax内，通过约束参数ξ使其小于误差阈值η；

其中，vS，vsP，vfP表示S波、慢P波及快P波的速度；

通过公式(6)，自适应地为不同速度选取相应的差分阶数，从而实现自适应差分阶数有限差分方法，减少运行时间且不降低精度；

4)在正演模拟过程中，引入并行算法

针对三维模型的有限差分方法，采用消息传递接口(MPI)及OpenMP，对算法进行并行化，全面缩短运行时间。

说明书

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

技术领域

[0001] 本发明属于地震勘探技术领域，涉及一种模拟方法，尤其是一种三维孔隙弹性介质中快速波场模拟方法，其针对三维孔隙弹性介质中高效的正演模拟技术。

背景技术

[0002] 我国的常规油气藏开采已接近尾声，开发的重点将转向非常规致密油气藏。地震波探测是进行油气探测最有效的方法之一。波场数值模拟技术是地震波探测的重要基础，也是认识储层地震波响应的关键工具，更是地震资料解释中参数反演和成像的前提条件。目前的数值模拟方法主要针对声学、弹性、粘弹性等简单的单相介质的波动方程，但是这些介质理论不能描述真实的地质特征，特别是不能准确描述非常规致密油气藏的储层特征。
孔隙弹性介质理论是现有的可准确描述真实油气藏物性特征的理论之一。因此开展基于孔隙弹性介质波动方程的波场数值模拟方法研究，具有重要的理论意义和应用价值。

[0003] 此外，由于油气储藏存在于三维空间，并且其横向和纵向都有非均匀性。因此必须开展三维数值模拟方法的研究。但是，三维实际地球物理模型尺寸庞大，并且相对于声学、弹性介质，孔隙弹性介质的波动方程更加复杂。因此三维孔隙弹性介质的数值模拟将耗费大量的运行时间，使得方法的运行效率大大降低。如何提高其运算效率是将三维孔隙弹性介质波场模拟应用于实际中首要解决的问题之一。

[0004] 有限差分方法实现简单，计算效率高，是最常用的数值模拟方法之一。有限差分方法是将波动方程中时间和空间导数用差商的形式表示，从而实现了时间和空间的离散。有限差分方法有很多种，包括显式和隐式有限差分方法，时间域和频率域有限差分方法，交错网格和旋转交错网格有限差分方法等等。

[0005] 为减少有限差分方法的运行时间，即提高有限差分方法的有效性，有很多出色的研究工作。Jensen采用变网格有限差分方法减少运行时间，同时可更好地刻画复杂边界。Tessmer在模型的不同区域，使用不同时间步长同样很可观地减少了计算时间。Wang在博士论文中，研究了变网格大小和变时间步长的有限差分方法，并将其应用到逆时偏移中，取得了可喜成果。Liu和Sen提出自适应的变差分阶数有限差分方法，根据模型中不同区域的速度来确定差分阶数。结果表明，该方法不仅能够减少运行时间，并且不会降低结果的准确性。

[0006] 然而，在提高三维孔隙弹性介质波场模拟的有效性方面，还未有相关方法与策略。

发明内容

[0007] 本发明的目的在于克服上述现有技术的缺点，提供一种三维孔隙弹性介质中快速波场模拟方法，该方法可减少正演模拟的运行时间，从而可应用于实际地球物理模型的正演模拟。

[0008] 本发明的目的是通过以下技术方案来实现的：

[0009] 这种三维孔隙弹性介质中快速波场模拟方法：首先，针对三维孔隙弹性介质波动方程，基于平面波理论和交错网格有限差分方法，得到三维孔隙弹性介质波动方程对应的数值频散关系；其次基于数值频散关系得到频散参数，通过使频散参数满足约束条件，实现不同速度采用不同的差分阶数，即变差分阶数有限差分方法；再次，引入基于MPI和OpenMP的并行算法，缩短运行时间。

[0010] 进一步，以上三维孔隙弹性介质中快速波场模拟方法具体包括以下步骤：

[0011] 1)三维孔隙弹性介质波动方程的数值频散关系推导

[0012] 基于平面波理论及交错网格有限差分方法，得到三维孔隙弹性介质波动方程的数值频散关系：

[0013]

[0014] 其中， 下标i＝1,2或3；当i＝1时表示慢P波，i＝2时表示S波，i＝3时表示快P波；k为波数，ω为角频率；θ是波传播方向与z轴夹角，φ是xoy平面与z轴夹角；2M为有限差分方法的阶数，τ为时间步长，h是空间网格大小，am(m＝1,2,...,M)是交错网格有限差分方法的差分系数；

[0015] 2)频散参数的定义

[0016] 数值频散由频散参数来衡量，表示为经有限差分方法离散后的速度与真实速度的比值：

[0017]

[0018] 其中：

[0019]

[0020] 当δ等于1时，没有数值频散；当δ大于或小于1时，会产生数值频散；

[0021] 频散参数定义为经有限差分方法离散后的波传播一个离散网格的时间与实际传播时间之差：

[0022]

[0023] 当ε等于0时，没有数值频散；当ε大于或小于0时，会产生数值频散；

[0024] 3)制定约束条件，实现自适应差分阶数有限差分方法

[0025] 定义快P波、S波及慢P波频散参数的参数ξ：

[0026]

[0027] 在频率fmax内，通过约束参数ξ使其小于误差阈值η；

[0028]

[0029] 通过公式(6)，自适应地为不同速度选取相应的差分阶数，从而实现自适应差分阶数有限差分方法，减少运行时间且不降低精度；

[0030] 4)在正演模拟过程中，引入并行算法

[0031] 针对三维模型的有限差分方法，采用消息传递接口(MPI)及OpenMP，对算法进行并行化，全面缩短运行时间。

[0032] 本发明具有以下有益效果：

[0033] 本发明的三维孔隙弹性介质中快速波场模拟方法针对的地质模型为三维孔隙弹性介质，基于频散参数，通过一定的约束条件使其小于误差阈值，从而为不同的速度选取相应的差分阶数，实现自适应差分阶数有限差分方法，有效减少运行时间，提高方法的有效性。

[0034] 进一步，本发明方便灵活，可提前为模型中不同速度计算相应的差分阶数，在正演模拟过程中不需要在边界处进行额外处理。

[0035] 进一步，本发明可有效减少正演模拟的运行时间，并且不会降低方法的精度。

[0036] 进一步，本发明引入并行算法，可进一步缩短运行时间。

附图说明

[0037] 图1是频散参数log10|ε|随频率变化曲线，其中τ＝0.001s，h＝25m；

[0038] 图2是不同速度、误差阈值的差分阶数，其中τ＝0.001s，h＝25m；

[0039] 图3是并行处理时数据划分示意图；

[0040] 图4是进程间通信示意图；

[0041] 图5是固体速度z分量在0.6s的波场快照，由(a)固定差分阶数(M＝11)(b)变差分阶数(M＝5,11)有限差分方法计算的结果，及(c)二者的差剖面；

[0042] 图6是三维推覆模型示意图；

[0043] 图7是基于推覆模型的固体速度z分量的地震记录，由(a)固定差分阶数(M＝12)，(b)变差分阶数(M＝4，4，4，5，6，6，12)及(c)二者差剖面，其中接收器位于平面y＝10km和平面z＝400m的交叉线上；

[0044] 图8是基于推覆模型的固体速度z分量的地震记录，由(a)固定差分阶数(M＝12)，(b)变差分阶数(M＝4，4，4，5，6，6，12)及(c)二者差剖面，其中接收器位于平面x＝10km和平面z＝400m的交叉线上；

[0045] 图9是基于推覆模型的固体速度z分量的地震记录，由(a)固定差分阶数(M＝12)，(b)变差分阶数(M＝4，4，4，5，6，6，12)及(c)二者差剖面，其中接收器位于平面x＝10km和平面y＝10km的交叉线上。

具体实施方式

[0046] 本发明针对的地质模型为三维孔隙弹性介质，基于该波动方程对应的数值频散关系，通过约束三种波的频散参数使其小于误差阈值，进而实现对不同速度采用不同的差分阶数，有效减少正演模拟的运行时间。首先基于平面波理论和交错网格有限差分方法，得到三维孔隙弹性介质波动方程的数值频散关系；其次基于数值频散关系定义频散参数，通过约束三种波的频散参数使其小于误差阈值，从而实现自适应差分阶数有限差分方法；最后，引入基于MPI和OpenMP的并行算法，进一步缩短三维孔隙弹性介质波场模拟的运行时间。

[0047] 下面结合附图和实例对本发明进行详细的描述。

[0048] 1)三维孔隙弹性介质波动方程的数值频散关系的推导

[0049] 首先，三维孔隙弹性介质的波动方程可写为

[0050]

[0051] 其中，

[0052] Q＝[vx，vy，vz，wx，wy，wz，τxx，τyy，τzz，τxy，τxz，τyz，s]T   (14)[0053] (vx,vy,vz)是固体速度分量，(wx,wy,wz)是流体速度分量，(τxx,τxz,τxy,τyy,τyz,τzz)是固体应力，s与流体压力有关。

[0054] 将上述时间空间域波动方程变换到频率波数域，可得该波动方程的频散关系：

[0055]

[0056] 这里，vS，vsP和vfP表示S波、慢P波及快P波的速度。

[0057] 基于平面波理论和交错网格有限差分方法及公式(9)，得到数值频散关系：

[0058]

[0059] 其中， 下标i＝1,2,3。当i＝1时表示慢P波，i＝2时表示S波，i＝3时表示快P波；k为波数，ω为角频率；θ是波传播方向与z轴夹角，φ是xoy平面与z轴夹角；2M为有限差分方法的阶数，τ为时间步长，h是空间网格大小，am(m＝1,2,...,M)是交错网格有限差分方法的差分系数。

[0060] 2)频散参数的定义

[0061] 数值频散可用频散参数衡量，基于公式(10)，频散参数可表示为有限差分方法离散后的速度与真实速度的比值

[0062]

[0063] 其中，

[0064]

[0065] 当δ等于1时，没有数值频散；当δ大于或小于1时，会产生数值频散。

[0066] 频散参数也可定义为经有限差分方法离散后的波传播一个离散网格的时间与实际传播时间之差

[0067]

[0068] 当ε等于0时，没有数值频散；当ε大于或小于0时，会产生数值频散。

[0069] 3)制定约束条件，实现自适应差分阶数有限差分方法

[0070] 为了衡量三维孔隙弹性介质的数值频散，定义包含了快P波、S波及慢P波频散参数的参数ξ

[0071]

[0072] 在一定频率fmax内，通过约束参数ξ使其小于误差阈值η

[0073]

[0074] 通过公式(15)，可以自适应地为不同速度选取相应的差分阶数，从而实现自适应差分阶数有限差分方法，有效减少运行时间但不会降低精度。图1给出固定差分阶数及自适应差分阶数有限差分方法计算的频散曲线。图1(a)中不同速度采用相同的差分阶数；而(b)中根据公式(15)为不同速度采用不同的差分阶数。对比两图可看出，为高速度采用低阶有限差分方法，并不会降低方法的精度，说明了本发明的可行性。

[0075] 图2是不同速度、不同误差阈值对应的差分阶数，可看出，速度和误差阈值越小，差分阶数越大。

[0076] 4)在正演模拟过程中，引入并行算法

[0077] 在并行计算中，计算区域被分成更小的子域，每个小子域由相应的进程进行处理；进程之间的通信由MPI完成。有限差分方法的实现过程中，仅仅需要区域边界周围的几个点，因此有限差分方法很容易进行并行化处理。

[0078] 一般而言，对于三维并行有限差分方法，计算区域应当在三个坐标轴上同时划分。并且应均匀划分，以平衡进程之间的工作量及通信量，如图3所示。

[0079] 图4为进程间通信的示意图，每个进程在自己的子区域处理问题，并在每个时刻点，相邻的进程进行通信并交换信息。为此，需在每个进程的边界处开辟一定大小的空间，以存储从相邻进程接收的信息。该空间的厚度等于差分阶数的一半，即M。

[0080] 此外，在程序实现时，空间循环部分可采用OpenMP，对线程进行优化，进一步提高运算效率。

[0081] 数值结果

[0082] 层状介质

[0083] 该模型大小为[0,5000m]×[0,5000m]×[0,5000m]，分界面位于深度1250米处，模型参数如表1所示。数值模拟时，时间步长为0.001s，网格大小为25m。震源为主频10Hz，时延0.1s的Ricker子波，位于模型中心。

[0084] 表1层状介质的参数

[0085]

[0086] 基于公式(21)为该模型中不同速度计算对应的差分阶数，层1的差分阶数为10，层2的差分阶数为22。22阶空间精度的有限差分方法模拟结果作为参考结果。固体速度z分量的波场快照如图5所示，由于模型是横向均匀的，因此只显示过震源的yoz平面的切片。从图中可看出，变差分阶数和固定差分阶数有限差分方法的数值结果无明显差别，二者的误差仅在0.7％。

[0087] 方法的有效性由运行时间来测试。表2列出了不同层厚情况下，两种有限差分方法的运行时间。可看出，本发明的运行效率与模型特征非常相关。高速区越大、低速区越小，有效性越高，即越省时。

[0088] 表2不同层厚的计算时间

[0089]

[0090] 简化的推覆模型

[0091] 下面通过模拟一个稍复杂模型来测试本发明的有效性。该模型为包含断层和推覆层的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作为固定差分阶数有限差分方法的差分阶数。

[0092] 表3推覆模型中不同速度对应的差分阶数

[0093]

[0094] 两种有限差分方法计算的固体速度z分量的地震记录及其差剖面如图7-9所示，其中接收器分别位于(…,10km,400m)，(10km,…,400m)，(10km,10km,…)。从这三幅图可以看出，两种有限差分方法计算的结果并无明显差异，并且误差分别仅为0.34％，0.34％，0.28％。因此，变差分阶数有限差分方法并不会降低方法的精度。

[0095] 表4位两种有限差分方法的运行时间，可以看出变差分阶数有限差分方法可以节省33.54％的时间，有效性得以大大提高。

[0096] 表4两种有限差分方法的运行时间

[0097]