纵、横波数值模拟方法及系统

本发明涉及一种纵、横波数值模拟方法及系统。

背景技术：

地震波在地下传播时，纵、横波通常会耦合在一起，所以利用完全弹性波方程进行数值模拟时，得到的是纵、横波的混合波场。利用纵、横波解耦的弹性波方程，可以单独得到纯纵波和纯横波波场，这对研究纵、横波在地下的传播特性以及数据处理有着重要的意义。例如，在利用多分量的地震数据进行逆时偏移时，纵、横波场的分离可以压制两者之间的串扰噪音，提高成像质量。

纵、横波解耦的弹性波方程有多种形式。马德堂和朱光明(2003)给出了均匀各向同性介质中的二阶弹性波方程，该方程可以得到完全分离的纯纵波或纯横波的位移波场。zhang等(2007)根据该方程推导出一阶速度-应力形式的分离方程，并利用交错网格有限差分法进行了数值模拟。xiao和leaney(2010)给出了另一种更为简便的分离方程。李振春等(2007)通过对应力进行分解，获得了非均匀介质中的一阶速度-应力形式的分离方程。tang等(2018)研究了界面处剪切模量对纵横波分离效果的影响，并将与剪切模量空间偏导数相关的“余项”在纵横波分量上重新分配，得到了与李振春等(2007)类似的分离方程。两者不同之处是，李振春等(2007)是将“余项”分配在了横波上，而tang等(2018)是将“余项”平均分配在纵波和横波上。

然而以上的分离方程，获得的纵、横波都是矢量场。在进行逆时偏移时，有效的矢量场成像条件就变成了一个亟待解决的问题。一种简单的做法是，将震源波场和检波点波场中纵横波的各个分量分别做互相关，如二维情况下，震源波场纵波和检波点波场横波的两个分量分别做互相关，得到4个偏移剖面。在三维情况下，将得到9个偏移剖面，然而如此多的偏移剖面，给地质解释造成了麻烦。yong等(2016)提出的内积成像条件，能够获得唯一的偏移剖面，但是其振幅是界面处入射波和反射波夹角的函数。

在二维情况下，利用散度和旋度公式可以得到纵、横波的总场，使得互相关成像变得简单。但是散度和旋度运算使得分离出来的纵、横波的振幅和相位都发生了改变。这些改变使得分离出来的纵波和横波的物理意义变得不明确。专利cn103412328a利用傅里叶变换，将波场变换到波数域，利用归一化波数进行保幅波场分离。专利zl201510559614.8修改了散度和旋度算子使得分离后的纵波和横波的物理意义变得明确，但是该方法得到的不是真实的纵、横波场，而是其时间一阶偏导数。

因此，现有技术中至少存在如下问题：已有的纵横波方程模拟出来的是纵、横波场的矢量场，不利于逆时偏移成像处理；利用散度和旋度算子分离的纵、横波，振幅和相位都发生了改变；波数域进行纵、横波分离，傅里叶正反变换会带来误差和干扰；获得的波场不是真实的纵、横波场，而是场其一阶时间偏导数。

技术实现要素：

有鉴于此，有必要提供一种纵、横波数值模拟方法及系统。

本发明提供一种纵、横波数值模拟方法，该方法包括如下步骤：a.根据各向同性介质中的弹性波方程，计算得到纵、横波分离的弹性波方程；b.从空间上对地质模型进行交错网格剖分，从时间上对总波场和纵、横波场进行交错网格剖分；c.根据纵、横波分离的弹性波方程和剖分出的网格，构建空间和时间差分格式；d.根据构建的空间和时间差分格式，计算得到纵波和横波波场。

其中，所述的步骤a具体包括：

二维情况下，令vx和vz分别表示v的水平分量和垂直分量，则得到纵、横波分离的弹性波方程(13)为：

其中，vp、vs分别是真实的纵、横波波场。

所述的步骤c具体包括：

对应的2n阶空间差分格式表示为：

其中，i和j分别表示的是网格点在z方向和x方向上的坐标；δx和δz分别表示x方向和z方向上的空间步长，n和n是正整数，表示的是精度为2n阶的第n个交错网格有限差分系数。

所述的步骤c具体包括如下步骤：

对应的2n阶空间差分格式表示为：

其中，i和j分别表示的是网格点在z方向和x方向上的坐标；δx和δz分别表示x方向和z方向上的空间步长，n和n是正整数，表示的是精度为2n阶的第n个交错网格有限差分系数。

所述的步骤c具体包括：

所对应的2n阶空间差分格式表示为：

其中，i和j分别表示的是网格点在z方向和x方向上的坐标；δx和δz分别表示x方向和z方向上的空间步长，n和n是正整数，表示的是精度为2n阶的第n个交错网格有限差分系数。

所述的步骤c具体包括：

对应的2n阶空间差分格式表示为：

其中，i和j分别表示的是网格点在z方向和x方向上的坐标；δx和δz分别表示x方向和z方向上的空间步长，n和n是正整数，表示的是精度为2n阶的第n个交错网格有限差分系数。

所述的步骤c具体包括：

令q表示vx或vz，r表示v^p或v^s，则2m阶时间差分格式表示为：

其中，t表示的是时间，δt表示的是时间步长，h表示的是时间步数。

所述的步骤c具体包括：

令q表示vx或vz，r表示v^p或v^s，则2m阶时间差分格式表示为：

其中，t表示的是时间，δt表示的是时间步长，h表示的是时间步数。

本发明提供一种纵、横波数值模拟系统，该系统包括弹性波方程计算模块、网格剖分模块、构建模块、纵、横波场计算模块，其中：所述弹性波方程计算模块用于根据各向同性介质中的弹性波方程，计算得到纵、横波分离的弹性波方程；所述网格剖分模块用于从空间上对地质模型进行交错网格剖分，从时间上对总波场和纵、横波场进行交错网格剖分；所述构建模块用于根据纵、横波分离的弹性波方程和剖分出的网格，构建空间和时间差分格式；所述纵、横波场计算模块用于根据构建的空间和时间差分格式，计算得到纵波和横波波场。

本发明纵、横波数值模拟方法及系统，其模拟出来的纵、横波是单个的标量场而不是矢量形式的两个分量，有利于逆时偏移成像；能够在时间-空间域获得真实的纵、横波场，且其振幅和相位保持不变；并且可以利用交错网格有限差分法进行求解，提高了计算精度。

此外，本发明还可以应用到弹性波的逆时偏移中，用来计算震源波场和检波点波场；还可以用到弹性波全波形反演中，用来反演地下的速度结构。

附图说明

图1为本发明纵、横波数值模拟方法的流程图；

图2(a)～(d)为本发明实施例提供的地质模型空间上的网格剖分4种形式示意图；

图3(a)、(b)为本发明实施例提供的地质模型时间上的网格剖分2种形式示意图；

图4为本发明纵、横波数值模拟系统的硬件架构图；

图5(a)～(d)为本发明实施例模拟出的0.35s时的地震波场示意图；

图6(a)、(b)为本发明实施例对角线上模拟出的纵、横波场和真实的纵、横波波场对比图；

图7(a)、(b)为非均匀介质模型的纵波速度模型、密度模型示意图；

图8为(a)～(d)为利用本发明实施例非均匀介质中模拟出0.5s时的地震波场示意图；

图9(a)、(b)为图8(a)～(d)在地表1.32km处模拟出的纵、横波场和真实的纵、横波波场对比图。

具体实施方式

下面结合附图及具体实施例对本发明作进一步详细的说明。

参阅图1所示，是本发明纵、横波数值模拟方法较佳实施例的作业流程图。

步骤s1，根据各向同性介质中的弹性波方程，计算得到纵、横波分离的弹性波方程。具体而言：

各向同性介质中的弹性波方程可以写为：

其中，λ和μ是拉梅常数，ρ是介质的密度，u是地震波位移波场总量，u^p和u^s分别是纵波位移波场和横波位移波场，表示哈密尔顿算子，t表示时间。

将方程(1)中的第二和第三项代入第一项可得：

利用关系式：

和

其中，α和β分别表示纵波和横波在地下介质中的传播速度。方程(2)变为：

设v为地震波速度波场总量，则其与u的关系式分别为：

设v^p为纵波速度波场，v^s为横波速度波场，其与u的关系式为：

将式(6)、式(7)和式(8)代入方程(5)，可得：

分别对式(7)和式(8)两侧求取一阶时间偏导数，并将式(6)代入可得：

将式(9)、式(10)和式(11)写成方程组的形式，得：

方程(12)便是纵横波分离的弹性波方程。二维情况下，令vx和vz分别表示v的水平分量和垂直分量，则方程(12)可以写为：

其中vs表示二维情况下的横波波场。

在频率-波数域可以证明vp和vs的振幅和相位未发生改变，就是真实的纵横波波场。将方程(13)中的最后两项变换到频率-波数域，可得：

其中，和分别表示频率-波数域中地震波速度波场的水平分量和垂直分量，和分别表示频率-波数域中的纵波和横波，ω表示圆频率，和分别表示纵波的水平波数和垂直波数，和分别表示横波波的水平波数和垂直波数。

根据频散关系，式(14)和式(15)两，可写为：

其中和分别表示纵波的单位水平波数和单位垂直波数，和分别表示横波的单位水平波数和单位垂直波数。由于表示的是纵波的单位偏振方向，表示的是横波的单位偏振方向，所以式(16)和式(17)是将和分别投影到了纵波和横波的单位偏振方向上，其振幅和相位未发生改变。从而证明方程(13)中的vp和vs就是真实的纵、横波波场。

步骤s2，从空间上对地质模型进行交错网格剖分，从时间上对总波场和纵、横波场进行交错网格剖分。具体而言：

在本实施例中，从空间上对地质模型进行交错网格剖分有4种形式，如图2(a)～(d)所示。从时间上对总波场和纵、横波场进行交错网格剖分有2种形式，如图3(a)、(b)所示。

在本实施例中，交错网格剖分的过程就是将地质模型进行离散化的过程，需要根据地质模型以及有限差分的精度来定。比方说一个1km×1km的地质模型，可以用5m的间隔进行剖分，也可以用10m的间隔进行剖分，这取决于计算精度，一般而言，网格间距越小，剖分的越密集，计算精度越高。

步骤s3，根据纵、横波分离的弹性波方程和剖分出的网格，构建空间和时间差分格式。具体而言：

步骤s31，根据纵、横波分离的弹性波方程和剖分出的网格构建空间差分格式。

为了表述方便，令vp＝ραv^p，vs＝ρβv^s则方程(13)可以写为：

图2(a)所对应的2n阶空间差分格式可以写为：

图2(b)所对应的2n阶空间差分格式可以写为：

图2(c)所对应的2n阶空间差分格式可以写为：

图2(d)所对应的2n阶空间差分格式可以写为：

式(19)～(34)中，i和j分别表示的是网格点在z方向和x方向上的坐标；δx和δz分别表示x方向和z方向上的空间步长，n和n是正整数，表示的是精度为2n阶的第n个交错网格有限差分系数。

步骤s32，根据纵、横波分离的弹性波方程和剖分出的网格构建时间差分格式。

所述的方程(18)中左侧的时间差分格式，有两种形式(如图3(a)、(b)所示)。具体为：令q表示vx或vz，r表示v^p或v^s，则图3(a)所示2m阶时间差分格式可以表示为：

当m＝1时，有

当m＝2时，有

当m≥3时，以此类推。

图3(b)所示2m阶时间差分格式可以表示为：

当m＝1时，有

当m＝2时，有

当m≥3时，以此类推。

式(35)～(46)中，t表示的是时间，δt表示的是时间步长，h表示的是时间步数。

步骤s4，根据构建的空间和时间差分格式，计算得到纵波和横波波场。具体而言：

从两种时间差分中选取一种替代纵横波分离的弹性波方程里的时间偏导数项，从四种空间差分格式中选取一种替代方程中的空间偏导数项，计算得到纵波和横波波场。

例如，如果选取图2(a)的空间差分格式，和图3(a)的时间差分格式，则的差分方程可以写为：

参阅图4所示，是本发明纵、横波数值模拟系统10的硬件架构图。该系统包括：弹性波方程计算模块101、网格剖分模块102、构建模块103、纵、横波场计算模块104。

所述弹性波方程计算101用于根据各向同性介质中的弹性波方程，计算得到纵、横波分离的弹性波方程。具体而言：

各向同性介质中的弹性波方程可以写为：

其中，λ和μ是拉梅常数，ρ是介质的密度，u是地震波位移波场总量，u^p和u^s分别是纵波位移波场和横波位移波场，表示哈密尔顿算子，t表示时间。

将方程(1)中的第二和第三项代入第一项可得：

利用关系式：

和

其中，α和β分别表示纵波和横波在地下介质中的传播速度。方程(2)变为：

设v为地震波速度波场总量，则其与u的关系式分别为：

设v^p为纵波速度波场，v^s为横波速度波场，其与u的关系式为：

将式(6)、式(7)和式(8)代入方程(5)，可得：

分别对式(7)和式(8)两侧求取一阶时间偏导数，并将式(6)代入可得：

将式(9)、式(10)和式(11)写成方程组的形式，得：

方程(12)便是纵横波分离的弹性波方程。二维情况下，令vx和vz分别表示v的水平分量和垂直分量，则方程(12)可以写为：

其中vs表示二维情况下的横波波场。

在频率-波数域可以证明vp和vs的振幅和相位未发生改变，就是真实的纵横波波场。将方程(13)中的最后两项变换到频率-波数域，可得：

其中，和分别表示频率-波数域中地震波速度波场的水平分量和垂直分量，和分别表示频率-波数域中的纵波和横波，ω表示圆频率，和分别表示纵波的水平波数和垂直波数，和分别表示横波波的水平波数和垂直波数。

根据频散关系，式(14)和式(15)两，可写为：

其中和分别表示纵波的单位水平波数和单位垂直波数，和分别表示横波的单位水平波数和单位垂直波数。由于表示的是纵波的单位偏振方向，表示的是横波的单位偏振方向，所以式(16)和式(17)是将和分别投影到了纵波和横波的单位偏振方向上，其振幅和相位未发生改变。从而证明方程(13)中的vp和vs就是真实的纵、横波波场。

所述网格剖分模块102用于从空间上对地质模型进行交错网格剖分，从时间上对总波场和纵、横波场进行交错网格剖分。具体而言：

在本实施例中，从空间上对地质模型进行交错网格剖分有4种形式，如图2(a)～(d)所示。从时间上对总波场和纵、横波场进行交错网格剖分有2种形式，如图3(a)、(b)所示。

在本实施例中，交错网格剖分的过程就是将地质模型进行离散化的过程，需要根据地质模型以及有限差分的精度来定。比方说一个1km×1km的地质模型，可以用5m的间隔进行剖分，也可以用10m的间隔进行剖分，这取决于计算精度，一般而言，网格间距越小，剖分的越密集，计算精度越高。

所述构建模块103用于根据纵、横波分离的弹性波方程和剖分出的网格，构建空间和时间差分格式。具体而言：

所述构建模块103根据纵、横波分离的弹性波方程和剖分出的网格构建空间差分格式。

为了表述方便，令vp＝ραv^p，vs＝ρβv^s则方程(13)可以写为：

图2(a)所对应的2n阶空间差分格式可以写为：

图2(b)所对应的2n阶空间差分格式可以写为：

图2(c)所对应的2n阶空间差分格式可以写为：

图2(d)所对应的2n阶空间差分格式可以写为：

式(19)～(34)中，i和j分别表示的是网格点在z方向和x方向上的坐标；δx和δz分别表示x方向和z方向上的空间步长，n和n是正整数，表示的是精度为2n阶的第n个交错网格有限差分系数。

所述构建模块103根据纵、横波分离的弹性波方程和剖分出的网格构建时间差分格式。

所述的方程(18)中左侧的时间差分格式，有两种形式(如图3(a)、(b)所示)。具体为：令q表示vx或vz，r表示v^p或v^s，则图3(a)所示2m阶时间差分格式可以表示为：

当m＝1时，有

当m＝2时，有

当m≥3时，以此类推。

图3(b)所示2m阶时间差分格式可以表示为：

当m＝1时，有

当m＝2时，有

当m≥3时，以此类推。

式(35)～(46)中，t表示的是时间，δt表示的是时间步长，h表示的是时间步数。

所述纵、横波场计算模块104用于根据构建的空间和时间差分格式，计算得到纵波和横波波场。

所述纵、横波场计算模块104从两种时间差分中选取一种替代纵横波分离的弹性波方程里的时间偏导数项，从四种空间差分格式中选取一种替代方程中的空间偏导数项，计算得到纵波和横波波场。

例如，如果选取图2(a)的空间差分格式，和图3(a)的时间差分格式，则的差分方程可以写为：

仿真表明本发明实用有效：

均匀介质中的模型参数为：纵波速度为3000m/s、横波速度为1700m/s、密度为1500kg/cm³。图5(a)～(d)为利用本发明实施例模拟出的0.35s时的地震波场，其中，图5(a)是利用本实施例模拟出的地震总波场的x分量(vx)；图5(b)是利用本实施例模拟出的地震总波场的z分量(vz)；图5(c)是本实施例模拟出的纵波波场(v^p)；图5(d)是本实施例分别是模拟出的横波波场(v^s)。

图6(a)、(b)为对角线(左上角到右下角)上模拟出的纵、横波场和真实的纵、横波波场对比图。图6(a)为纵波波场对比图，其中的实线为利用本发明计算出的纵波波场v^p，点线为利用vx和vz在角度域分离出来的真实的纵波波场。图6(b)为横波波场对比图，其中的实线为利用本发明计算出的纵波波场v^s，点线为利用vx和vz在角度域分离出来的真实的横波波场。

图7(a)、(b)为非均匀介质模型的纵波速度模型、密度模型，其中：图7(a)为纵波速度模型，横波速度为纵波速度的图7(b)为密度模型。

图8(a)～(d)为利用本发明实施例非均匀介质中模拟出0.5s时的地震波场。其中，图8(a)和图8(b)分别为总波场的x分量(vx)和z分量(vz)；图8(c)和图8(d)分别是模拟出的纵波波场(v^p)和横波波场(v^s)。

图9(a)、(b)为图8(a)～(d)在地表1.32km处模拟出的纵、横波场和真实的纵、横波波场对比图。图9(a)为纵波波场对比图，其中的实线为利用本发明计算出的纵波波场v^p，点线为利用vx和vz在角度域分离出来的真实的纵波波场。图9(b)为横波波场对比图，其中的实线为利用本发明计算出的纵波波场v^s，点线为利用vx和vz在角度域分离出来的真实的横波波场。

从对比图上可以看出，模拟出的纵、横波与真实的纵、横波有着非常好的吻合度，充分说明本发明是实用有效的。

虽然本发明参照当前的较佳实施方式进行了描述，但本领域的技术人员应能理解，上述较佳实施方式仅用来说明本发明，并非用来限定本发明的保护范围，任何在本发明的精神和原则范围之内，所做的任何修饰、等效替换、改进等，均应包含在本发明的权利保护范围之内。