井间地震粘弹性波场特征的数值模拟研究

井间地震粘弹性波场特征的数值模拟研究（精选2篇）

篇1：井间地震粘弹性波场特征的数值模拟研究

井间地震粘弹性波场特征的数值模拟研究

以数值模拟方法为基础,研究了井间地震粘弹性波场特征及传播规律,并重点研究了井间地震波的粘滞衰减和散射衰减.从本构方程和运动方程出发,给出了二维粘弹性波方程交错网格任意偶数阶有限差分数值模拟方法.对均匀模型、含洞均匀背景模型和含洞层状模型进行了井间地震粘弹性波场数值模拟.结果表明,该方法模拟精度高,边界吸收效果好.介质的.粘滞性使得波振幅明显衰减、波形和相位畸变严重,主频向低频偏移显著,且有效频带变窄.洞的散射波同相轴为双曲型,但其同相轴曲率大于直达波的曲率.洞的散射引起直达波振幅衰减明显,而对直达波主频的影响较小.当洞的直径小于一个波长时,洞散射波的主频高于直达波的主频,且与洞的大小关系密切.

作 者：宋常瑜 裴正林 Song Changyu Pei Zhenglin  作者单位：中国石油大学CNPC物探重点实验室,北京,100083 刊 名：石油物探  ISTIC PKU英文刊名：GEOPHYSICAL PROSPECTING FOR PETROLEUM 年，卷(期)： 45(5) 分类号：P631.4 关键词：井间地震   粘弹性方程   粘滞衰减   散射衰减   交错网格   有限差分  

篇2：井间地震粘弹性波场特征的数值模拟研究

上世纪80年代以来, 制作VSP模型日益受到勘探学者们的重视, 很多学者陆续发表了一些相关文献[1]。使用一维波动方程的分解来合成一维VSP炮记录是K.D.Wyatt (1981) [2]首次提出的SVSP方法;D.C.Ganley (1981) [3]在合成一维粘弹性VSP炮记录时对地层的吸收和品质因子的频散效应进行了考虑;K.D.Wyaff和S.B.Wyatt (1981) [4]又提出了二维的射线追踪模型, 并用于VSP-CDP迭加;随后更多的学者开始VSP合成记录方法的研究。

1 射线追踪正演模拟VSP炮记录

射线追踪方法可以快速有效的进行波场的近似模拟[5]。射线的概念来源于光学, 射线的几何扩散和传播时间的概念早在300年以前就已经被惠更斯和费马提出了[6]。将射线方法应用于地震学比该方法本身晚。早期的射线法在地震学和地震勘探中只有运动学方面的应用, 直到1958年, 巴比奇与阿列克谢耶夫, 以及1959年Karal和Keller分别独立地将其在电磁学中的研究成果引入到弹性波领域。此后, 很多学者又对射线级数法的理论及其应用作了进一步的研究[7]。

1.1 射线追踪方法的原理

根据射线追踪的理论, 每个射线追踪方法都是围绕斯奈尔定理、费马原理、惠更斯原理展开的[8]。下面用改进的试射法正演模拟出粘弹性VSP炮记录。

(1) 非零井源距观测系统

将非零井源距的震源固定在距离观测井具有某一距离的位置上, 检波器位于垂直方向, 然后根据目的层的埋藏深度来设定井源距的大小。非零井源距观测系统如图一所示, 从图中可以看出井源距与可能观测到的范围大小之间的关系。

本文采用改进的试射法作两点间射线追踪[9]。改进的试射法的基本思路是:用组射线追踪的方法代替单条射线追踪的方法。组射线追踪就是一次发出几百条甚至上万条的一组射线 (角度差很小) , 利用这组射线进行射线追踪, 然后利用射线追踪计算每条射线的横向移动距离。在众多数值中寻找最接近偏移距的值, 并将该值对应的出射角记录下来;计算横向移动距离和偏移距之间的差值, 如果差值在误差范围内, 则将该初射角当作纵波的入射角, 以此追踪出波的射线路径和旅行时间。如果不在误差范围内, 则改变出射角的角度差使其更小, 再次计算, 直到差值在误差范围内为止。在具体的计算过程中, 分以下步骤来实现。

(2) 射线入射角

在2D xoz空间, 当介质是横向各向同性介质时, 地震波在同一层中的射线路径是直线如图二所示, 假设炮检距为SN=x, 炮点坐标为S (0, 0) , 在界面上, 射线产生反射和投射, 并服从Snell定律[8]。

式 (1) 中, β0、βi-2、βi-1分别是入射角和透射角, V0、Vi-2、V1分别是界面两侧的纵波 (或横波) 速度, 第j个检波器的深度为Hj。反射波满足方程:

联立式 (1) 和式 (2) , 并且有, 因为VSP模型的入射角在 (0°, 90°) , 因此可以得到:

求解 (3) 可以解出出射到第j个检波器的入射波的入射角β0j。按照该思路可以求出不同检波器一次反射波的入射角, 同时也可以求出直达波的入射角。分段求出射线与界面的交点, 图二中 (xi, zi) 是传播射线与第i层界面的交点坐标。

(3) VSP响应的计算

VSP响应的计算包括两个部分:旅行时的计算与地震波振幅的计算[9]。

①旅行时的计算

旅行时的计算如下:

式 (4) 中, N是射线总段数, 第i段射线路径的始点与终点分别是 (xi-1, zi-1) 和 (xi, zi) , Vi表示第i段射线上的地震波的速度。

②振幅的计算

振幅响应使用粘弹性Zoeppritz方程在频率域内进行计算, 本文只计算透射波和一次反射波。首先通过粘弹性Zoeppritz方程计算反射系数和透射系数, 然后合成频率域地震记录, 最后反变换回时间域得到时间域的合成地震记录。

(4) 反射系数和透射系数的计算

本文以P波入射到粘弹性界面的Zoeppritz方程为基础, 计算得到反射系数和透射系数。

地震波的传播受到地层的吸收作用, 因此, 岩石不但有弹性还具有粘滞性。当平面P波入射到粘弹性分界面R时, 产生反射P波、透射P波、反射SV波和透射SV波, 如图三所示。图三中的参数vp1, vp2, vs1, ρ1, ρ2, θp1, θs1, θp2, θs2, Qp1, Qs1, Qp2, Qs2, γ分别表示界面的上下两侧介质的纵横波速度、地层密度、入射角、反射角和透射角、品质因子以及衰减角[8]。

正演模拟时, 速度模型、各层密度、震源位置、每个检波器的坐标都是已知的, 射线与界面的交点是待求的参数。通过改进的试射法计算得到射线与速度密度界面的交点和入射角。这里采用的粘弹性介质Zoeppritz方程为[10]:

式 (5) 中, k表示波数, γ, μ是拉梅系数, Rpp, Rps, Tpp, Tps分别是地震波入射到粘弹性界面的纵波和横波的反射系数和透射系数。其中, 复波数k和γ, μ的表达式为:

对于复波数k来说, P为传播矢量, A是衰减矢量, γ是随波的传播发生变化的衰减角, 但是当纵波的入射角小于临界角时, 反射系数受衰减角的影响非常小;只有入射角位于临界角的附近位置, 反射系数受到的影响才比较大。本节中假定介质1和介质2中的衰减角不变, 均为γ。

(5) 合成地震记录

弹性介质中平面波U (x, ω) 传播一段距离△x后可以表示为[11]:

在粘弹性介质中, 平面波的传播过程中涉及到地震波的衰减, 因此假设震源点波场为U0 (0, ω) , 传播一段距离x后表示为:

式 (7) 中, R (ω) 为由粘弹性精确Zoeppritz方程计算得到的反射系数或透射系数, 在波数k (ω) 引入Q值 (品质因子) 来表达地层的Q衰减效应[12]:

将 (8) 代入 (7) 中, 利用时间τ代替旅行距离x, 可得到波场表达式为:

把不同频率代入式中可得到τ时刻的波场值, 然后将其变换到时间域可得到时间域地震记录。

(6) 射线追踪流程

首先, 将震源处的地震波变换到频率域;其次, 进行射线追踪, 得到地震波传播的入射角;然后, 计算精确的Zoeppritz方程, 得到纵波的反射系数和透射系数;再次, 计算地震波传播的射线长度, 在频率域合成粘弹性介质的地震记录;最后, 将频率域变换到时间域, 得到时间域的合成地震记录。

利用上面的方法原理, 在具体的实现过程中, 射线追踪的实现流程如图四所示。

1.2 VSP地震波场模拟与特征分析

对图五中的三层水平层状介质进行VSP正演模拟。该模型在xoz坐标系中, 其中, 该模型为505m×1855m, 横向采样间隔为5m, 纵向检波点间距为5m, 模型的网格点数为101×371。第一层厚度为752m, 品质因子Q为200, 第二层厚度为1107m, 品质因子Q为30, 第三层厚度为1907m, 品质因子Q为150。各层纵横波速度及密度参数如表一所示。对图五所示的速度模型进行正演模拟, 采用左边放炮的激发方式, 炮检距为50m, 第一个检波点位于距离地表5m的地方, 共371个检波器。在保证计算稳定的前提下选取1ms的时间采样间隔, 选用主频为25Hz的Ricker子波。正演模拟得到的VSP炮记录如图五所示。

图六所示的是弹性介质VSP炮记录和粘弹性介质VSP炮记录, 其中 (a) 中没有增益, (b) 中增益为1。图七和图八分别是对应的波形图。由于该模型只包含两个反射界面, 所以在合成的记录中有两层上行反射波与下行波同相轴在界面的深度位置有交点。因为是水平层状的速度模型, 所以不同界面的反射波同相轴是相互平行的。通过对比图六 (a) 和 (b) 以及对应的波形图可以发现, 粘弹性介质中传播的地震波受到大地滤波作用能量明显减少。为了更直观的观察地震波振幅和相位的变化, 在这里抽取几道地震数据分别在时间域和频率域作比较, 如图九至图十一所示。从图九至图十一中可以看出, 随地震波传播时间 (深度) 的增加, 地震波振幅明显减小, 有效频带范围变窄, 从而导致地震波能量衰减的增大。

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

有限差分方法 (Finite Difference Method, 简称FDM) 是建立在差分原理基础之上的一种数值计算方法。使用有限差分方法对波动方程进行正演模拟的基本思想是:首先用适于数值计算的差分代替波动方程中的微分, 得到相应的差分方程, 然后由给定的初始条件对差分方程进行求解, 就能根据空间中过去时刻波场值的分布及其变化推导出当前时刻的波场值分布。

使用有限差分方法对波动方程进行正演模拟时, 首先需要对模拟区域进行网格划分, 常见的网格划分方法有矩形网格、三角形网格、六边形网格、交错网格以及平面型广义网格等等。为了提高数值模拟精度, 降低数值频散效应, 改善模拟效果, 我们采用一种交错网格高阶差分方法来模拟对称轴沿Z轴的横向各向同性 (VTI) 介质中弹性波的传播。首先介绍了这种数值方法的基本原理与算法, 将这种方法应用于二维弹性波数值模拟中, 研究了这种交错网格有限差分法的稳定性条件, 分析了这种模拟方法的稳定性和频散效应。寻找到了适用于二维速度—应力弹性波方程的吸收边界条件, 并从理论和实例上说明了它们的吸收效果。

2.1 方法原理与算法

在VTI介质中, 如果只考虑x方向和z方向上的分量而不考虑y方向上的分量, 就得到了二维粘弹性波方程, 该方程中包含应力分量和速度分量对时间的导数, 其速度—应力方程可表示如下:

式 (10) ~ (14) 中的参数含义为:vx和vz分别表示水平方向和垂直方向的速度分量;σxx和σzz分别为水平方向和垂直方向上的正应力;σxz表示剪切应力;ρ表示介质的密度;γ和μ表示拉梅常数, γ'和μ'表示粘滞系数。

在VTI介质中, 使用交错网格有限差分方法对弹性波方程进行正演模拟的过程一般包括以下四个步骤:

第一步:计算定义在主体网格上的应力分量和应变分量;

第二步:将速度和应力交错地分布在模型中的各个网格点上, 速度和应力在模型网格点上的具体位置如图十二所示;

第三步:用二阶或者高阶差分来代替速度—应力粘弹性波方程中所有的一阶偏导数, 求得相应的差分方程;

第四步:对第三步得到的差分方程按照时间一步一步地向前推进, 计算各个网格点上的速度或者应力。

下面分别介绍时间和空间上的高阶差分格式的推导过程, 首先介绍时间上的2M阶精度的差分方程格式的推导。

交错网格的高阶有限差分方法的一阶速度应力粘弹性波方程同时在时间和空间上交错, 在t+△t/2时刻和t时刻分别计算速度和应力, 其中, △t表示时间步长。我们用Taylor公式将在t处展开[11], 分别见式 (15) 和式 (16) :

用式 (15) 减去式 (16) 就可以得到精度为2M阶的时间差分格式[11], 见式 (17) :

当2M=2, 即M=1时, 将M=1代入公式 (17) 就可以得到传统的二阶精度时间差分格式。本文用交错网格有限差分方法对粘弹性波波动方程进行正演模拟时就是采用传统的二阶精度时间差分[13], 其中速度vx、vz和应力σzz、σzz、σxz具体的二阶精度时间差分格式见式 (18) ~ (22) 。

如果对公式 (17) 直接进行计算, 在计算时就会牵扯到好几个时间层, 需要较大的运算量。因此可以转嫁速度 (应力) 对时间的奇数阶的高阶的导数为应力 (速度) 对空间的导数, 以减小计算量。这样在计算时就不会涉及到过多的时间层, 从而节省了内存。

下面介绍空间上的2N阶精度差分格式的推导。

我们利用公式 (23) 来计算方程 (10) ~ (14) 中的一阶偏导数:

式 (23) 中, f表示vx、vz、σxx、σzz或者σxz, △x表示空间步长, Cn (N) 表示差分权系数。为了保证一阶空间导数的差分精度达到2N阶, 需要准确地求取差分权系数。我们用Taylor公式将在x处展开, 对方程组 (24) 进行求解就能得到不同阶数的差分权系数。

下面给出几种常用的不同阶数的空间差分权系数:

当2N=4时, C1 (2) =1.125, C2 (2) =-0.04166667;当2N=6时, C1 (3) =1.171875, C2 (3) =-0.06510416, C3 (3) =-0.0046875;当2N=8时, C1 (4) =1.196289, C2 (4) =-0.0797526, C3 (4) =0.009570313, C4 (4) =-0.0006975447;当2N=10时, C1 (5) =1.211243, C2 (5) =-0.08972168, C3 (5) =0.001384277, C4 (5) =-0.00176566, C5 (5) =0.0001186795。

对模型区域进行离散之后, 采用上文给出的交错网格时间和空间上的差分格式, 令x=i△x, z=i△z, t=k△t, 其中, i和j表示空间网格点, k表示时间网格点。我们用分别表示速度vx、vz与应力σxx、σzz、σxz的离散值。为了方便起见, 假定△x=△z, 则一阶速度—应力粘弹性波方程的时间二阶、空间2N阶精度的差分格式[14]见式 (25) ~ (32) :

式中, i和j是指空间离散指数, n是时间离散指数, △t是时间离散步长;Rxx, Rzz和Rxz为记忆变量, 引入的目的是消除松弛方程中的褶积以绕开耗时的褶积运算[15]。在推导方程 (30) ~ (32) 时, 需要考虑下面的关系[14]:

本文在进行正演模拟时选用的时间差分精度为二阶, 空间差分精度为十阶, 所以这种平均作法并不影响数值模拟的精度[16]。同时式 (29) 中的最后一项运算需要考虑 (35) 。方程 (27) ~ (32) 中的系数表达式[14]分别如下所示:

式 (36) 中, τpε和τsε分别是指纵波和横波的应变松弛时间, τσ是应力松弛时间, 下面是松弛时间和品质因子之间的关系表达式[17]:

QP和QS分别表示纵波和横波的品质因子, f0表示子波的中心频率。方程 (25) ~ (32) 中和方向的离散化微分算子分别是Dx和Dz, 其表达式为:

式 (38) 中, am是有限差分加权系数, △x和△z分别指x方向的空间步长和z方向的空间步长, N是倒数算子长度, 对于空间十阶精度的情况, N取5。具体的差分权系数已在上文给出。

2.2 稳定性与频散分析

粘弹性介质的稳定性条件和弹性介质的类似[16]。使用有限差分方法进行正演模拟时, 如果某些参数的选取不合理, 有可能会产生按指数形式增长的没有实际物理意义的数值, 造成模拟结果网格频散严重, 甚至还有可能导致模拟结果的数据溢出进而无法正常计算。所以, 为了确保模拟算法的稳定性, 应该计算出参数的合理取值范围。在进行正演模拟时, 差分方程的稳定性分析是必须要考虑的问题, 数值稳定是一个好的差分格式的必要条件。

本文在进行正演模拟时, 选用董良国 (2000) 提出的一阶弹性波方程交错网格高阶差分稳定性条件[18]。在TI介质中, 二维情况下的稳定性条件见式 (39) 和式 (40) :

式中, 2M表示时间差分的阶数, 2N表示空间差分的阶数, △t表示时间步长, △x和△z分别表示水平方向和垂直方向上的空间步长。

令, 则公式 (39) 和公式 (40) 可以表示为:

上式即为TI介质中一阶弹性波方程交错网格有限差分时间2M阶、空间2N阶精度的稳定性条件。本文进行正演模拟时, 选用的时间差分的精度为二阶, 空间差分的精度为十阶, 即M=1、N=5, 代入公式 (41) 和公式 (42) 可以计算出稳定性条件为:

又因为, 其中, Vp表示纵波速度, Vs表示横波速度, 将其代入上面两式就可以得到由纵波速度、横波速度、时间步长和空间步长表示的稳定性条件, 见式 (45) 和式 (46) :

使用交错网格有限差分方法进行正演模拟时, 最大的纵波速度和最大的横波速度是已知的, 所以我们可以由上述稳定性条件选择合理的时间步长和空间步长来保证正演模拟的稳定性。

2.3 边界条件设置

实际中的地震波是在半无限空间介质中传播的, 但是由于受到计算机存储容量以及计算速度等运行环境的限制, 在无限的区域内进行正演模拟是不切实际的, 我们只能限定一个有限的模拟区域, 这样模拟区域的边界即成为人工边界。如果在正演模拟时对人工边界不作任何处理, 那么地震波传播到人工边界时就会产生很强的边界反射, 对正演模拟的结果产生很大的干扰。因此, 消除边界反射就成为正演模拟中一个比较关键的问题。比较经典的边界条件主要包括吸收边界条件和衰减边界条件两大类, 常用的边界条件有波动方程分解法、旁轴近似法、阻尼衰减法和完全匹配层法。本文在正演模拟时选用阻尼衰减边界条件。

由于阻尼衰减边界条件的形式比较简单, 易于编程实现, 成为应用最为广泛的边界条件之一。Cerjan等人 (1985) 在模拟各向同性介质中的地震波时首先使用了衰减边界条件[19]。它的基本原理是在计算区域的周围设置一层衰减带, 如图十三所示, 在正演模拟的过程中, 对衰减区域每个网格点上的数值解都乘以一个衰减因子, 见式 (47) :

式 (47) 中, u (x, z, t) 表示引入阻尼项之前波动方程的解, U (x, z, t) 表示引入阻尼项之后波动方程的解, G (x, z) 表示衰减因子。这样, 地震波传播到边界时就会在衰减区域迅速发生衰减, 并不会在边界处产生较强的边界反射。

阻尼衰减边界条件的衰减效果主要取决于衰减函数的选取, 下面给出了几种比较常用的衰减函数:

(1) Cerjan给出了指数型衰减函数G (i) =e-α2 (L-i) 2。式中, L表示衰减带的宽度, i表示计算点距离边界的网格点数, α表示衰减系数。通过大量的试验, Cerjan给出了一个经验值:i0=20, α=0.015[20]。

(2) 周竹生 (2007) 给出了衰减函数。式中, L表示衰减带的宽度, i表示计算点距离边界的网格点数[21]。

(3) 杜启振 (2009) 给出了余弦衰减函数以及抛物线衰减函数。式中, L表示衰减带的宽度, i表示计算点距离边界的网格点数, β和γ表示衰减系数[22]。

衰减函数的选取准则为:衰减函数的取值由内向外逐渐减小, 其一阶导数连续, 当总衰减量和衰减带的宽度确定时, 主要衰减区 (即衰减函数的取值为1到0.4) 覆盖的网格数越多越好[23]。图十四是上述几种衰减函数的对比图, 我们可以看到当总衰减量较小时, 三种衰减函数曲线的走势大致相同, 而当总衰减量较大时, 抛物线衰减函数最符合衰减准则。因此, 本文在正演模拟时选用抛物线衰减函数。

2.4 VSP地震波场模拟与特征分析

本文使用交错网格高阶有限差分法对一阶速度—应力粘弹性波方程进行正演模拟, 差分精度选用传统的时间二阶、空间十阶。模拟过程中, 边界条件选用阻尼衰减边界条件, 衰减函数选用抛物线衰减函数, 震源选用Ricker子波, 其表达式见下式:

式 (48) 中, f表示震源的主频。

对图五中的三层水平层状介质模型采用交错网格有限差分法进行正演模拟, 其中, 模型的大小为505米×2255米, 横向采样间隔为5米, 纵向采样间隔为5米, 即模型的网格点数为101×451。速度模型中各层的厚度、纵波速度、横波速度、品质因子和密度如表一所示。

对图五所示的速度模型进行正演模拟, 采用左侧放炮的激发方式, 震源放置在地面, 其位于模型网格点中坐标为 (0, 0) 的位置上。利用上文给出的稳定性条件, 在保证计算稳定的前提下选取1毫秒时间采样间隔, 总记录长度为2.0秒。选用主频为25Hz的Ricker子波进行正演模拟, 检波器垂直分布, 第一个检波器的深度为5米, 这样利用交错网格有限差分法模拟得到弹性和粘弹性介质VSP单炮地震记录, 如图十五和图十六所示。

图十五和图十六分别是弹性和粘弹性VSP记录的x分量和z分量。图十七 (a) 和图十七 (b) 对应的分别是弹性VSP炮记录z分量和粘弹性VSP炮记录z分量的波形图。可以清晰的辨认出图十五和图十六中包含的直达纵波、转换横波、反射纵波等多种波场现象。可以注意到VSP炮记录上直达波能量很强, 而反射波能量与直达波能量相比要小很多。从图十五和图十六中清楚地看到, x分量存在的干扰因素较多, 因此采用VSP炮记录的z分量。为了能够更直观的对比有限差分法得到的弹性VSP炮记录和粘弹性VSP炮记录, 现从炮记录中抽取单道进行时间域和频率域比较, 如图十八至图二十所示。

图十八至图二十分别为从Z分量炮记录中抽取的第50道、150道、300道地震数据在时间域和频率域的响应。从以上三个图中可以发现, 随着检波器深度的增加, 地震波传播时间逐渐增大, 粘弹性介质中地震波的振幅衰减也越来越大, 从频率域显示结果发现, 有效频带范围随着地震波传播时间的增长而减小。

3 结束语

对于粘弹性VSP炮记录的正演模拟, 我们分别采用了射线法和交错网格有限差分法, 在射线法中我们只考虑了直达波和一次反射波, 不考虑转换波的反射, 而在交错网格有限差分法中则包含了转换波反射, 在有限差分数值模拟中, 存在转换波的反射, 这对直达波和纵波反射波也存在一定的影响, 而在射线追踪方法中只考虑了纵波透射和反射。射线追踪方法的缺点是对物性参数变化较大的速度模型不能很好的进行地震波场模拟。

摘要：本文采用射线法和交错网格有限差分法对粘弹性VSP炮记录进行正演模拟, 分别观察地震波振幅、相位、能量和有效频带范围的变化情况, 从而得出粘弹性介质中传播的地震波场的特征情况。