基于切比雪夫伪谱法的各向异性衰减面波模拟方法

摘要

本发明公开基于切比雪夫伪谱法的各向异性衰减面波模拟方法，包括：获取浅层地球物理参数数据、定义观测系统和震源函数；根据自由界面垂向和切向应力连续条件，由波动方程获得自由界面处面波波场更新方程；在垂直和水平方向采用切比雪夫和傅里叶伪谱法计算自由界面处面波波场更新方程的垂向和水平导数；四阶龙格?库塔离散自由界面处面波波场更新方程的各向异性衰减介质面波方程的时间微分；计算模型水平方向两个边界的波场吸收；将波动方程分解为入射波方程和出射波方程，并得到底边界的波场更新方程，且据以计算当前时刻的面波，以产生面波炮记录和波场快照。通过本发明，以解决现有技术存在的浅层各向异性衰减面波模拟问题。

说明书

技术领域

本发明涉及地震勘探的技术领域，尤其涉及一种基于切比雪夫伪谱法的各向异性 衰减面波模拟方法。

背景技术

高精度三维地震技术在各个油田得以成功的推进和创新，我国东部老油田，例如， 胜利油田、大庆油田等均处于平原地区，而其他主要产区，例如塔河、鄂尔多斯盆地多为沙 漠、黄土塬等复杂地表区，探区内资源潜力大，巨厚、疏松的黄土层和浅层沙漠地层不利于 有效地震波的产生，且浅地表对地震波的吸收衰减作用强烈，严重影响地震数据处理。

但浅层面波发育，面波波场包含浅地表横波速度、介质对地震波的吸收强弱、各向 异性等信息，通过弹性介质的面波已经可实现反演浅层横波速度信息，并可将反演结果应 用到静校正流程，但对面波的各向异性衰减研究有限，浅层各向异性衰减的面波模拟研究， 可用于建立浅层介质衰减和速度各向异性关系，是我国黄土塬、沙漠地区复杂地质三维地 震数据处理技术的理论基础。

发明内容

本发明的主要目的在于提供一种，以解决现有技术存在的浅层各向异性衰减面波 模拟的问题。

为解决上述问题，本发明实施例提供一种基于切比雪夫伪谱法的各向异性衰减面 波模拟方法，包括：获取浅层地球物理参数数据、定义观测系统和震源函数；根据自由界面 垂向和切向应力连续条件，由波动方程获得自由界面处面波波场更新方程；在垂直方向采 用切比雪夫伪谱法计算所述自由界面处面波波场更新方程的垂向导数；在水平方向采用傅 里叶伪谱法计算所述自由界面处面波波场更新方程的水平导数；所述自由界面处面波波场 更新方程的各向异性衰减介质面波方程的时间四阶龙格-库塔离散计算；所述自由界面处 面波波场更新方程的模型水平方向的两个边界波场吸收计算；将所述波动方程分解为入射 波方程和出射波方程，并将入射波在底边界吸收掉，从而得到底边界的波场更新方程；根据 所述底边界的波场更新方程，计算当前时刻的面波，以产生面波炮记录和波场快照；输出所 述面波炮记录和波场快照。

根据本发明的技术方案，通过获取浅层地球物理参数数据、定义观测系统和震源 函数，再根据自由界面垂向和切向应力连续条件，获得自由界面处面波波场的更新方程，接 着，垂直方向采用切比雪夫伪谱法计算导数，水平方向采用傅里叶伪谱法计算导数，各向异 性衰减介质面波的时间离散计算，模型水平方向的左右边界吸收计算，而模型的底边界，采 用类似于自由边界的方法，将反射到模型内部的地震波场吸收掉，从而得到底界面的波场 更新方程，且根据底界面的波场更新方程计算面波波场，并输出面波炮记录和波场快照。如 此一来，有效地解决了浅层各向异性衰减面波模拟问题，尤其适用于各向异性衰减介质面 波的正演分析工作，还能够快速实现各向异性衰减介质面波模拟，准确记录面波正演结果， 并提高地质解释的精度。

附图说明

此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，本发 明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：

图1是根据本发明实施例的基于切比雪夫伪谱法的各向异性衰减面波模拟方法的 流程图；

图2是根据本发明实施例的切比雪夫伪谱法网格剖分的示意图；

图3是根据本发明实施例的一维切比雪夫不均匀网格剖分的示意图；

图4a和图4b分别是根据本发明实施例的勘探区采集的微测井数据的示意图；

图5a和图5b分别根据本发明实施例的面波炮记录的示意图；

具体实施方式

本发明的主要思想在于，基于通过获取浅层地球物理参数数据、定义观测系统和 震源函数，再根据自由界面垂向和切向应力连续条件，获得自由界面处面波波场的更新方 程，接着，垂直方向采用切比雪夫伪谱法计算导数，水平方向采用傅里叶伪谱法计算导数， 各向异性衰减介质面波的时间离散计算，模型水平方向的左右边界吸收计算，而模型的底 边界，采用类似于自由边界的方法，将反射到模型内部的地震波场吸收掉，从而得到底界面 的波场更新方程，且根据底界面的波场更新方程计算面波波场，并输出面波炮记录和波场 快照。如此一来，有效地解决了浅层各向异性衰减面波模拟问题，尤其适用于各向异性衰减 介质面波的正演分析工作，还能够快速实现各向异性衰减介质面波模拟，准确记录面波正 演结果，并提高地质解释的精度。

为使本发明的目的、技术方案和优点更加清楚，以下结合附图及具体实施例，对本 发明做进一步地详细说明。

根据本发明的实施例，提供了一种基于切比雪夫伪谱法的各向异性衰减面波模拟 方法。

图1是根据本发明实施例的基于切比雪夫伪谱法的各向异性衰减面波模拟方法的 流程图。

在步骤S102中，获取浅层地球物理参数数据、定义观测系统和震源函数。其中，所 述浅层地球物理参数数据包括纵波速度、横波速度、Thomson各向异性参数、纵波品质因子 和横波品质因子。所述观测系统包括震源位置、检波器起始位置、检波器个数、检波器间距， 定义地震模拟的总时间长度和时间采样率，记录波场快照的时间间隔。所述震源函数包括 震源类型和震源主频，并且震源类形还包括爆炸震源和点震源。

进一步来说，震源类型可如式(1)、式(2)所示：

$\overrightarrow{F}=Z\overrightarrow{j},---\left(1\right)$

$\overrightarrow{F}=\frac{\partial \phi }{\partial x}\overrightarrow{i}+\frac{\partial \phi }{\partial y}\overrightarrow{j},---\left(2\right)$

其中，表示体力项，Z表示体力的垂向分量，分别为x和y方向的单位向量， φ表示一标量函数，是关于时间的函数，x、y分别表示相互垂直的二维坐标轴。

在步骤S104中，根据自由界面垂向和切向应力连续条件，由波动方程获得自由界 面处面波波场更新方程。

首先，根据各向异性衰减理论，采用标准线性模型表征介质粘弹性，通过粘弹性波 动方程中引入的记忆变量表征介质衰减效应，二维直角坐标系下，各向异性衰减介质qP- qSV波动方程，如式(3)、(4)、(5)和(6)所示：

$\begin{array}{c}{\stackrel{·}{v}}_x=\frac{1}{\rho }(\frac{\partial {\sigma }_{xx}}{\partial x}+\frac{\partial {\sigma }_{xz}}{\partial z})+f_x\\ {\stackrel{·}{v}}_z=\frac{1}{\rho }(\frac{\partial {\sigma }_{xz}}{\partial x}+\frac{\partial {\sigma }_{zz}}{\partial z})+f_z\end{array},---\left(3\right)$

其中，v_x、v_z为质点在x、z方向的速度分量，σ_xx、σ_zz、σ_xz为应力分量，ρ为介质的密度， f_x、f_z表示外力；

$\begin{array}{c}{\stackrel{·}{\sigma }}_{xx}={\hat{c}}_{11}{\stackrel{·}{ϵ}}_{xx}+{\hat{c}}_{13}{\stackrel{·}{ϵ}}_{zz}+2c_{15}{\stackrel{·}{ϵ}}_{xz}+Y\underset{l=1}{\overset{L_1}{\Sigma }}{\stackrel{·}{e}}_{1l}+c_{55}\underset{l=1}{\overset{L_2}{\Sigma }}{\stackrel{·}{e}}_{2l}\\ {\stackrel{·}{\sigma }}_{zz}={\hat{c}}_{13}{\widehat{ϵ}}_{xx}+{\hat{c}}_{33}{\stackrel{·}{ϵ}}_{zz}+2c_{35}{\stackrel{·}{ϵ}}_{xz}+Y\underset{l=1}{\overset{L_1}{\Sigma }}{\stackrel{·}{e}}_{1l}-c_{55}\underset{l=1}{\overset{L_2}{\Sigma }}{\stackrel{·}{e}}_{2l}\\ {\stackrel{·}{\sigma }}_{xz}=c_{15}{\stackrel{·}{ϵ}}_{xx}+c_{35}{\stackrel{·}{ϵ}}_{zz}+2{\hat{c}}_{55}{\stackrel{·}{ϵ}}_{xz}+c_{55}\underset{l=1}{\overset{L_2}{\Sigma }}{\stackrel{·}{e}}_{3l}\end{array},---\left(4\right)$

其中，

$\begin{array}{c}{\hat{c}}_{11}=c_{11}-X+{\mathrm{YM}}_{u1}+c_{55}{M}_{u2}\\ {\hat{c}}_{13}=c_{13}+2c_{55}-X+{\mathrm{YM}}_{u1}-c_{55}{M}_{u2}\\ {\hat{c}}_{33}=c_{33}-X+{\mathrm{YM}}_{u1}+c_{55}{M}_{u2}\\ {\hat{c}}_{55}=c_{55}{M}_{u2}\\ X=\frac{c_{11}+c_{33}}{2},Y=X-c_{55}\end{array},---\left(5\right)$

其中，表示未松弛的介质弹性系数，c_ij表示松弛的介质弹性系数，为应变的时间导数，e_1l、e_2l、e_3l为记忆变量，M_uv，v＝1,2为广义标准线性模型的松 弛函数，v＝1对应P波松弛函数，v＝2对应SV波松弛函数，X、Y为常数；

$\begin{array}{cc}{\stackrel{··}{e}}_{1l}={\phi }_{1l}(\frac{\partial v_x}{\partial x}+\frac{\partial v_z}{\partial z})-\frac{{\stackrel{·}{e}}_{1l}}{{\tau }_{\sigma l}^{\left(1\right)}},& l=1,\mathrm{...}{L}_1\\ {\stackrel{··}{e}}_{2l}={\phi }_{2l}(\frac{\partial v_x}{\partial x}-\frac{\partial v_z}{\partial z})-\frac{{\stackrel{·}{e}}_{2l}}{{\tau }_{\sigma l}^{\left(2\right)}},& l=1,\mathrm{...}{L}_2\\ {\stackrel{··}{e}}_{3l}={\phi }_{2l}(\frac{\partial v_x}{\partial z}+\frac{\partial v_z}{\partial x})-\frac{{\stackrel{·}{e}}_{3l}}{{\tau }_{\sigma l}^{\left(2\right)}},& l=1,\mathrm{...}{L}_2\end{array},---\left(6\right)$

其中，为介质的应变与应力松弛时间，为记 忆变量的一阶、二阶导数。

然后，根据各向异性衰减波动方程对应的一阶双曲方程，将波动方程在自由界面 的波场分解为入射分量和出射分量，从已知波场计算入射波，并依据自由界面σ_zz＝σ_xz＝0， 相应的入射波对应的特征值为λ₂＝-c_p,λ₄＝-c_s，则可获得自由界面面波波场更新方程，如 式(7)所示：

其中，上标(new)表示重新计算的波场，上标(old)表示最初由波动方程计算的波 场。

在步骤S106中，在垂直方向采用切比雪夫伪谱法计算所述自由界面处面波波场更 新方程的垂向导数。通过切比雪夫伪谱法计算，式(3)、(4)、(6)可转化为式(8)、(9)、(10)， 如下所示：

其中，式(8)、式(9)和式(10)实现了垂向导数的切比雪夫网格剖分计算，这种网格 在模型边界的网格间距小，模型中间网格间距最大，利于提高自由界面面波模拟精度。

进一步的，切比雪夫不规格网格剖分图如图2所示，在垂直方向的两个边界网格密 集，可以提高面波模拟的精度，而且式(7)所示的面波更新方程尤其适用于这种边界。切比 雪夫不规格网格剖分的区间是[-1,+1]之间，如图3所示，在区间内由余弦函数定义式(11)， 如下所示：

$\begin{array}{cc}{x}_i=cos\left(\frac{\pi i}{M}\right),& i=0,\mathrm{...},M,\end{array}---\left(11\right)$

通过切比雪夫多项式离散函数f(x)后，可得式(12)，如下所示：

$\begin{array}{cc}f\left(x_i\right)=\underset{j=0}{\overset{M}{\Sigma }}{b}_j{T}_j\left(x_i\right),& i=0,...M,\end{array}---\left(12\right)$

并且，函数f(x_i)关于x的空间导数如式(13)所示：

$\frac{df\left(x_i\right)}{dx}=\underset{k=0}{\overset{M}{\Sigma }}{d}_k{T}_k\left(x_i\right),---\left(13\right)$

其中，b_j、d_k为其系数，T_j(x_i)为切比雪夫函数，T_k(x_i)为余弦函数。切比雪夫网格剖 分区间为引入扩展函数，并将区间[-1,1]映射到实际物理区间[0,z_max]，则在实 际物理区间内，函数的空间如式(14)所示：

其中，为区间扩展函数。

在步骤S108中，在水平方向采用傅里叶伪谱法计算所述自由界面处面波波场更新 方程的水平导数。在本实施例中，假设某函数为f(x)，根据傅里叶变换性质，水平导数计算 过程如式(15)所示：

其中，FFT、FFT^-1分别为正反傅里叶变换，F(w)为f(x)的傅里叶变换，k_x为波数，通 过式(15)可计算函数f(x)关于x的空间导数。

在步骤S110中，所述自由界面处面波波场更新方程的各向异性衰减介质面波方程 的时间四阶龙格-库塔离散计算。在本实施例中，各向异性衰减介质面波方程的时间离散计 算是采用四阶龙格-库塔方法。

首先，将各向异性衰减介质波动方程改写成式(16)，如下所示：

$-\frac{\partial V}{\partial t}+A\frac{\partial V}{\partial x}+B\frac{\partial V}{\partial z}+D=0,---\left(16\right)$

其中，矢量V为{v_x,v_z,σ_xx,σ_zz,σ_xz,e_1l,e_2l,e_3l}^T，分别为关于时间、方 向x和z的偏导数，A、B表示弹性系数和函数φ_1l、φ_2l的矩阵，D表示波动方程(9)中的常量 矩阵。

然后，各向异性衰减介质波动方程的时间微分表达式变为式(17)，如下所示：

$V^{n+1}=V^n+\frac{1}{6}dt(H_1+2H_2+2H_3+H_4),---\left(17\right)$

并且，式(17)中的参数如式(18)所示：

$\begin{array}{c}{H}_1={\mathrm{MV}}^n+D^n\\ H_2=M(V^n+\frac{dt}{2}{H}_1)+D^{n+1/2}\\ H_3=M(V^n+\frac{dt}{2}{H}_2)+D^{n+1/2}\\ H_4=M(V^n+{\mathrm{dtH}}_3)+D^{n+1}\end{array},---\left(18\right)$

其中，H₁表示Vⁿ点处的斜率；H₂表示利用Δ₁求得的(Vⁿ+dt/2)点处的斜率；H₃表示利 用H₂求得的(Vⁿ+dt/2)点处的斜率；H₄表示利用Δ₃求得的(Vⁿ+dt)点处的斜率；M表示速度在 空间三个方向的导数算子，Dⁿ表示各向异性衰减方程中常量在开始时刻(第n时刻)的数值， D^n+1/2表示常量在时间段中点(第n+1/2时刻)的数值，Dⁿ⁺¹表示常量在时间段终点(第n+1时 刻)的数值。

在步骤S112中，所述自由界面处面波波场更新方程的模型水平方向的两个边界波 场吸收计算。在本实施例中，模型水平方向的两个边界波场吸收计算是采用卷积完全匹配 层(CPML)吸收边界法，如式(19)所示：

$\frac{\partial u}{\partial x}=\frac{\partial u}{\partial x}+{\xi }_x\left(t\right)*\frac{\partial u}{\partial x}=(1+b)\frac{\partial u}{\partial x}+a\left({\xi }_x\right(t)*\frac{\partial u}{\partial x}),---\left(19\right)$

其中，参数a和b分别为边界吸收系数，u为匹配层处的地震波场，ξ_x(t)为匹配层内 线性变换的衰减函数。

在步骤S114中，将所述波动方程分解为入射波方程和出射波方程，并将入射波在 底边界吸收掉，从而得到底边界的波场更新方程。其中，所述底边界的波场更新方程如式 (20)所示：

${\stackrel{·}{v}}_x^{\left(new\right)}=\frac{1}{2}({\stackrel{·}{v}}_x^{\left(old\right)}-\frac{1}{{\mathrm{\rho c}}_s}{\stackrel{·}{\sigma }}_{xz}^{\left(old\right)})$

${\stackrel{·}{v}}_z^{\left(new\right)}=\frac{1}{2}({\stackrel{·}{v}}_z^{\left(old\right)}-\frac{1}{{\mathrm{\rho c}}_p}{\stackrel{·}{\sigma }}_{zz}^{\left(old\right)})$

${\stackrel{·}{\sigma }}_{xx}^{\left(new\right)}={\stackrel{·}{\sigma }}_{xx}^{\left(old\right)}-\frac{{\hat{c}}_{13}}{2{\hat{c}}_{33}}({\stackrel{·}{\sigma }}_{zz}^{\left(old\right)}+{\mathrm{\rho c}}_p{\stackrel{·}{v}}_z^{\left(old\right)})$

${\stackrel{·}{\sigma }}_{zz}^{\left(new\right)}=\frac{1}{2}({\stackrel{·}{\sigma }}_{zz}^{\left(old\right)}-{\mathrm{\rho c}}_p{\stackrel{·}{v}}_z^{\left(old\right)})$

${\stackrel{·}{\sigma }}_{xz}^{\left(new\right)}=\frac{1}{2}({\stackrel{·}{\sigma }}_{xz}^{\left(old\right)}-{\mathrm{\rho c}}_s{\stackrel{·}{v}}_x^{\left(old\right)})$

${\stackrel{··}{e}}_{1l}^{\left(new\right)}={\stackrel{··}{e}}_{1l}^{\left(old\right)}-\frac{{\phi }_{1l}}{2{\hat{c}}_{33}}({\stackrel{·}{\sigma }}_{zz}^{\left(old\right)}+{\mathrm{\rho c}}_p{\stackrel{·}{v}}_z^{\left(old\right)})$

${\stackrel{··}{e}}_{2l}^{\left(new\right)}={\stackrel{··}{e}}_{1l}^{\left(old\right)}+\frac{{\phi }_{2l}}{2{\hat{c}}_{33}}({\stackrel{·}{\sigma }}_{zz}^{\left(old\right)}+{\mathrm{\rho c}}_p{\stackrel{·}{v}}_z^{\left(old\right)})$

${\stackrel{··}{e}}_{3l}^{\left(new\right)}={\stackrel{··}{e}}_{3l}^{\left(old\right)}-\frac{{\phi }_{2l}}{2{\hat{c}}_{55}}({\stackrel{·}{\sigma }}_{xz}^{\left(old\right)}+{\mathrm{\rho c}}_s{\stackrel{·}{v}}_x^{\left(old\right)}),---\left(20\right)$

其中，c_p、c_s分别为弹性情况下纵波和横波速度，上标(new)表示重新计算的波场， 上标(old)表示最初由波动方程计算的波场。

在步骤S116中，根据所述底边界的波场更新方程，计算当前时刻的面波，以产生面 波炮记录和波场快照。在步骤S118中，输出所述面波炮记录和波场快照。

上述已说明了基于切比雪夫伪谱法的衰减各向异性面波模拟方法，以下将提供一 些实例来验证上述方法的适用性及准确性。

图4a和图4b是根据本发明实施例的勘探区采集的微测井数据的示意图。利用微测 井单道记录，拾取初至，并转化为垂直的t0时间，再与距离绘制时深曲线，然后，利用时深曲 线计算、划分速度界面，从而得到近地表速度模型，模型定义为v_P＝1826m/s，v_S＝1000m/s，ρ ＝2300g/cm³，Q_P＝100，Q_S＝60。Thomson各向异性参数可以选择为：ε＝0.195，δ＝-0.22；模 型计算尺寸为1000×1000m，水平方向网格间距dx＝1m，垂直方向最大网格间距dz＝2m，震 源位于N_x＝500m和N_z＝500m，震源类型采用爆炸震源，子波主频为45Hz。通过本实施例的上 述步骤进行模拟后，可获得面波炮记录，如图5a所示。在图5a和图5b中，R表示各向异性衰减 介质的面波波场，ux为水平分量，uz为垂向分量。

综上所述，根据本发明的技术方案，通过获取浅层地球物理参数数据、定义观测系 统和震源函数，再根据自由界面垂向和切向应力连续条件，获得自由界面处面波波场更新 方程，将当前时刻的地震波场更新为各向异性衰减面波的传播方程，接着，垂直方向采用具 有不规格网格剖分的切比雪夫伪谱法计算导数，水平方向采用规格网格剖分的傅里叶伪谱 法计算导数，各向异性衰减介质面波的时间微分采用四阶龙格-库塔法离散，针对二维模型 水平方向的左右边界，采用非分裂CPML吸收较强的边界反射，而模型的底边界，采用类似于 自由边界的方法，将反射到模型内部的地震波场吸收掉，从而得到底界面的波场更新方程， 以上步骤完成后，叠代计算面波波场，并输出面波炮记录和波场快照。如此一来，有效地解 决了浅层各向异性衰减面波模拟问题，尤其适用于各向异性衰减介质面波的正演分析工 作，还能够快速实现各向异性衰减介质面波模拟，准确记录面波正演结果，并提高地质解释 的精度。

显然，本领域的技术人员应该明白，上述的本发明的基于切比雪夫伪谱法的衰减 各向异性面波模拟方法和各步骤可以用通用的计算装置来实现，它们可以集中在单个的计 算装置上，或者分布在多个计算装置所组成的网络上，可选地，它们可以用计算装置可执行 的程序代码来实现，从而，可以将它们存储在存储装置中由计算装置来执行。这样，本发明 不限制于任何特定的硬件和软件结合。存储装置为非易失性存储器，如：ROM/RAM、闪存、磁 盘、光盘等。

以上所述仅为本发明的实施例而已，并不用于限制本发明，对于本领域的技术人 员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、 等同替换、改进等，均应包含在本发明的权利要求范围之内。