三维粘弹性介质地震波场数值模拟技术研究

陈玉玺，王 璐，章 阳

（1.中国海洋大学，山东青岛266100；2.江西地震局，江西南昌330039）

三维粘弹性介质地震波场数值模拟技术研究

陈玉玺*1，王 璐1，章 阳2

（1.中国海洋大学，山东青岛266100；2.江西地震局，江西南昌330039）

为了适应油气田勘探的发展需求，提高勘探精度，地震勘探由二维向三维发展。传统的地震波数值模拟一般假设地层介质是均匀、完全弹性的，但在实际情况中，由于介质是非完全弹性、各向异性、物性参数可连续变化的，因此地震波的能量在传播过程中是逐渐衰减的。在非均匀介质情况下进行地震波的数值模拟，可以更真实地反映地震波的传播和能量变化规律。推导出三维交错网格高阶有限差分算法，求取优化差分算子，实现三维粘弹性波数值模拟，获得地震记录，并分析其波场特征和能量变化情况。

非完全弹性；交错网格高阶有限差分；差分算子；数值模拟

地震波数值模拟技术是研究复杂油气田地球物理勘探最基本也是最重要的一种方法，是反演的基础。数值模拟的精度也决定了反演的精度，对于我们认识地质构造，油气田勘探，识别裂隙型油藏具有重要的指导作用。三维复杂粘弹性介质的数值模型是地震波数值模拟的主要工作，也是当今地震勘探的重点。

1 三维粘弹性波方程的建立

目前使用最广泛的粘弹性介质模型为佛各特（Viogt）模型[1]，它将粘弹性介质的应力分为2部分：完全弹性应力和粘性应力。

假设λ和μ为完全弹性拉梅常数，λ′和μ′为粘性拉梅常数，则粘弹性波动方程中的拉梅常数等效为将完全弹性波方程中的拉梅常数等效替换成粘弹性波方程中的拉梅常数，就得到了粘弹性波方程。本文主要研究一阶速度-应力形式的三维粘弹性波方程，在均匀各向同性介质中（假设体力为0），其表达式为式（1）。

其中，λ和μ是弹性拉梅常数，λ′和μ′为是粘性拉梅常数，我们知道根据前人对粘弹性介质理论的研究有：由这几个式子可以得出式（1）中的拉梅常数：

其中，ρ是介质的密度，vp和vs是介质的纵波和横波速度，Qp和Qs是纵、横波品质因子，w是圆频率。

2 交错网格高阶有限差分法

交错网格法在常规网格的基础上引入半网格点，可以有效地处理速度分量和应力分量之间的耦合关系，提高数值模拟的精度和稳定性。有限差分法中用离散化的差商近似代替连续的微商，这样就会产生无法避免的数值频散，而高阶有限差分法可以有效地降低由于网格离散化造成的数值频散，提高数值模拟的精度。本文用交错网格高阶有限差分法对三维粘弹性波方程进行数值模拟，其网格剖分如图1所示。

图1 三维粘弹性介质交错网格示意图

其中将正应力分量σxx、σyy、σzz定义在点( )i,j,k上；剪应力分量σxy定义在点上，σxz定义在点上，σyz定义在点上；速度分量vx和速度变化率分量定义在点上，vy和定义在点上，vz和定义在点上。

在地震波场数值模拟中，有几个问题是需要着重考虑的，首先是稳定性条件，本文采用裴正林等提出的稳定性条件，其表达式为[3]：

式中：Δx、Δy、Δz——空间x、y、z三个方向的网格间距；

Δt——时间网格间距；

vmax——模型最大速度；

am——高阶差分系数，本文选取优化差分算子（以时间2阶，空间10阶为例）：a1=1.250，a2=-0.1200，a3=0.0314，a4=-0.0092，a5=0.0018。

其次是吸收边界条件，实际地震勘探中，地震波往往是在无限空间的地下地层中传播。而在用计算机数值模拟时，我们所选择的理论物理模型通常是在有限区域范围内，这样便产生了人工边界反射问题[4]。本文采用PML吸收边界条件实现对三维粘弹性波方程的边界处理[5]。PML的核心思想是在计算区域外面构造有限厚度的吸收层，吸收或衰减向外传播的波，同时在计算模型区域与吸收层的链接处是透明的，使得产生最小可能的虚假反射来解决模拟时的人工边界问题。

三维模型需要在6个面、12条棱和8个角点处分别添加PML吸收边界。在x、y、z三个方向均加载衰减项d(x)、d(y)、d(z)，其计算公式为：

其中，R为理论反射系数，δ为PML厚度[2]。在8个角点区域，d(x)≠0，d(y)≠0，d(z)≠0；在12条棱区域，xoy平面和 yoz平面相交的棱域，d(x)≠0，d(y)=0，d(z)≠0，xoy平面与xoz平面相交的棱域，d(x)=0，d(y)≠0，d(z)≠0，xoz平面与yoz平面相交的棱域，d(x)≠0，d(y)≠0，d(z)=0；在6个面区域，垂直x轴的2个平面，d(x)≠0，d(y)=0，d(z)=0，垂直y轴的2个平面，d(x)=0，d(y)≠0,d(z)=0,垂直z轴的2个平面，d(x)=0,d(y)=0,d(z)≠0。这样整个模型区域都引入PML吸收边界，地震波传播经过模型边界到达PML吸收层时，能量逐渐减弱，因此不会产生任何边界反射。

以vx和σxx为例，根据PML的分裂思想，则三维粘弹性波方程的精度为Ο(Δ t2+Δx2N)的引入PML吸收边界条件的交错网格高阶差分格式为[6]：

其中：

其中：

3 三维复杂粘弹性介质模型

在三维French速度模型的基础上[7]，加入纵横波品质因子Qp、Qs从而构建三维French粘弹性介质速度模型。模型包含1个水平层面、1个倾斜界面和2个隆起构造，整个模型σ=0.25，ρ=2000kg/m3。模型大小为1000m × 1000m × 1000m，空 间 网 格 步 长Δx=Δy=Δz=5m，时间步长Δt=0.5ms，震源子波选择主频为30Hz的雷克子波，计算精度为O(Δt2+Δx10)，图2为相应各个平面的二维垂直切片的纵波速度模型。

图2 各平面上的二维垂直剖面的纵波速度

图3、图4、图5分别为t=500ms时，不同平面上的波场快照，分析各个平面上的粘弹性波波场快照切片，由于三维波场之间的相互干涉和相互影响，形成的波场十分复杂。xoz平面上，隆起构造的侧面反射波和透射波明显，能量强，yoz平面上同样能观察到能量很强的侧面反射波和透射波，xoy平面上能观察到2个隆起构造产生的垂直反射波信息和多次侧面反射波信息，还产生能量较强的散射波。随着传播深度的不断增大，由于存在粘滞吸收衰减作用，能量产生一定的衰减，而且记录到的散射波的能量明显比反射波弱。

图3 t=500ms时刻y=500m处xoz平面上的粘弹性波波场快照

4 结论

以French速度模型为基础，加入纵横波品质因子Qp、Qs构建三维粘弹性介质速度模型，以粘弹性介质为例，根据其应力-应变关系，结合完全弹性波方程和表征介质吸收衰减特性的品质因子Q，推导出三维粘弹性波方程，并进一步给出其一阶速度-应力关系式。实验发现，地震波在粘弹性介质中传播时存在明显的能量衰减现象，地震波的振幅明显衰减，反射波同相轴能量明显减弱，同时随着地震波传播深度的不断增加，能量衰减越发严重。反射波波形发生畸变，而且吸收衰减作用越明显，地震波波形畸变越严重。并且这种粘滞吸收衰减作用对横波影响比对纵波强。同时本文研究的三维粘弹性交错网格高阶有限差分法模拟技术不仅能够更详细全面地描述地震波的传播及其全波场特征，而且该算法稳定，无明显的频散现象，计算精度高，对以后实际的三维地震勘探技术会有一定的帮助。

图4 t=500ms时刻x=500m处yoz平面上的粘弹性波波场快照

图5 t=500ms时刻z=500m处xoy平面上的粘弹性波波场快照

[1]Stanke F E,Bumidge R.Comparison of Spatial and Ensemble Averaging Methods Applied to Wave Propagation in Finely Lay⁃ered Media[C]//60th Ann.Internat.Mtg.，Soc.expl.Geophys.Ex⁃ panded Abstracts，1990:1062-1065.

[2]董敏煜.地震勘探[M].东营：中国石油大学出版社，2006.

[3]牟永光，裴正林.三维复杂介质地震数值模拟[M].北京：石油工业出版社，2005.

[4]陆基孟.地震勘探原理[M].东营：中国石油大学出版社，2009.

[5]Collino F,Tsogka C.Application of the Perfectly Matched Ab⁃sorbing layer Model to the Linear Elastodynamic Problem in Anisotropic Heterogeneous Media[J].Geophysics，2001，66(1): 294-307.

[6]董良国，马在田，曹景忠.一阶弹性波方程交错网格高阶差分解法稳定性研究[J].地球物理学报，2000，43(6):856-864.

[7]Saenger E H，Gold N，Shapiro S A.Modeling the Propagation of Elastic Waves Using a Modified Finite-difference Grid[J]. Wave Motion，2000(31):77-92.

P631

A

1004-5716(2016)12-0026-04

2016-03-04

2016-03-07

陈玉玺（1990-），男（汉族），山东临沂人，中国海洋大学在读硕士研究生，研究方向：油气田与煤田地球物理勘探方法与技术。