CN104750954A - 一种在复杂各向异性介质中模拟地震波的方法及装置 - Google Patents

Abstract

本发明提供了一种在复杂各向异性介质中模拟地震波的方法及装置，方法包括：给数值化的三维地质模型的每个网格点的弹性常数和密度赋值；获取三维地质模型中网格点沿三个坐标轴的速度分量的一阶空间导数；将速度分量的一阶空间导数代入到弹性介质中弹性波方程组中的第一方程式，并进行求解获取应力分量；获取三维地质模型中网格点沿三个坐标轴的应力分量的一阶空间导数；将应力分量的一阶空间导数代入到弹性介质中弹性波方程组中的第二方程式，并进行求解获取速度分量；根据通过第二方程式获取的速度分量得到当前时间第一时刻的速度分量的值，作为当前时刻的地震波的记录；获取下一时刻的地震波的记录，直至获取时间长度t的地震波的波场值。

Description

一种在复杂各向异性介质中模拟地震波的方法及装置

技术领域

本发明涉及有限差分菱形交错网格优化领域，特别涉及一种在复杂各向异性介质中模拟地震波的方法及装置。

背景技术

关于数值求解波动方程问题主要是体现在如何求解波动方程中空间导数。以下从弹性波方程求解问题的角度出发，探讨目前常用的求解算法。在弹性介质中弹性波方程的表现形式可用一阶速度应力方程来描述，如下：

$\left\{\begin{array}{c}{\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{xx}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{yy}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{zz}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{yz}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{xz}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{xy}}\end{array}\right\}=\left\{\begin{array}{cccccc}{c}_{11}& c_{12}& c_{13}& c_{14}& c_{15}& c_{16}\\ c_{12}& c_{22}& c_{23}& c_{24}& c_{25}& c_{26}\\ c_{13}& c_{23}& c_{33}& c_{34}& c_{35}& c_{36}\\ c_{14}& c_{24}& c_{34}& c_{44}& c_{45}& c_{46}\\ c_{15}& c_{25}& c_{35}& c_{45}& c_{55}& c_{56}\\ c_{16}& c_{26}& c_{36}& c_{46}& c_{56}& c_{66}\end{array}\right\}\left\{\begin{array}{c}{v}_{x,x}\\ v_{y,y}\\ v_{z,z}\\ v_{y,z}+v_{z,y}\\ v_{z,x}+v_{x,z}\\ v_{y,x}+v_{x,y}\end{array}\right\}---\left(1\right)$

$\left\{\begin{array}{c}\mathrm{\ρ}{\stackrel{\·}{v}}_x={\mathrm{\σ}}_{\mathrm{xx},x}+{\mathrm{\σ}}_{\mathrm{xy},y}+{\mathrm{\σ}}_{\mathrm{xz},z}\\ \mathrm{\ρ}{\stackrel{\·}{v}}_y={\mathrm{\σ}}_{\mathrm{xy},x}+{\mathrm{\σ}}_{\mathrm{yy},y}+{\mathrm{\σ}}_{\mathrm{yz},z}\\ \mathrm{\ρ}{\stackrel{\·}{v}}_z={\mathrm{\σ}}_{\mathrm{xz},x}+{\mathrm{\σ}}_{\mathrm{zy},y}+{\mathrm{\σ}}_{\mathrm{zz},z}\end{array}\right.---\left(2\right)$

速度应力公式含有两组方程，即式（1）和式（2）。在式（1）和式（2）中，c和ρ表示介质的弹性常数和密度；σ表示波场的应力分量，其双下标表示应力的6个分量，其第三下表表示对x、y或z的一阶空间导数。v表示波场的质点速度分量，其单下标表示应力的3个分量，其第二下表表示对x,y,z的一阶空间导数，其上面一点表示一阶时间导数。求解式（1）和式（2）主要体现在如何求解方程右边的18个空间导数。

常用求解上述波动方程的数值方法有虚谱法，传统的有限差分算法，和有限差分交错网格方法。

虚谱法使用离散FFT求空间导数，这样，虚谱法精度高，但是计算量大，尤其在实现大型三维地质模型的模拟，大型三维FFT变换计算量太大，实现计算不现实。

传统的有限差分算法的特点为数值模型中只有一种数值网格，使用中心差分实现简单，方便快捷，但是精度低。

有限差分交错网格方法的特点为数值模型中含多种数值网格，相互交错，使用中心差分，这样克服了传统算法的缺点，即快捷，精度高，尤其是提高差分算子阶数之后，能有效提高计算精度。按照波场分量不同的交错方式，有限差分交错网格方法又分为三种方法，即标准交错网格（Standard Staggered Grid，简称SSG），旋转交错网格（Rotated Staggered Grid，简称RSG），和菱形交错网格（Diamond Staggered Grid，简称DSG）。

a、标准交错网格（SSG）

Virieux首先提出了二维空间中弹性波方程标准交错网格。如图1所示，为标准交错网格三维示意图。在三维标准交错网格中，有7种不同的网格（每个数值代表一个网格）相互交错在一起，这样交错的目的是方便使用中心差分方法求解在各个应力分量位置的速度的一阶空间导数（即方程（1）的解），以及在各个速度分量位置的应力分量的一阶空间导数（即方程（2）的解）。例如：求图1中方格2位置（应力分量位置）对x方向的一阶导数，即可使用左右两个“7”网格位置的速度分量值通过中心差分求解。标准交错网格容易实现，计算快捷，并且占用内存小，目前在实际应用中多数声波和弹性波数值模拟使用标准交错网格。标准交错网格适合于常规声波方程，各项同性介质中弹性波方程的模拟，也适合介质对称性高于正交各向异性介质（包含正交各向异性介质）中弹性波的模拟，但是在解决介质对称性低于正交各向异性（如单斜或三斜介质）时，在求解某些波场分量的空间导数时标准交错网格需要对其插值。为了说明这一问题，这里举例说明，以方程（1）中第一个方程为例，

${\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{xx}}=c_{11}{v}_{x,x}+c_{12}{v}_{y,y}+c_{13}{v}_{z,z}+c_{14}(v_{y,z}+v_{z,y})+c_{15}(v_{x,z}+v_{z,x})+c_{16}(v_{x,y}+v_{y,x})---\left(3\right)$

在图1中，为了计算应力位置某些速度分量的一阶空间导数（方程（3）右边6项中后面3项对应的空间导数），需要求解这些分量在其他位置的空间导数，之后再插值到应力位置。这样引入的问题就是导入了插值计算误差，而且这种误差会传播到下一时间片波场的计算之中，从而降低了计算的精度，另外标准交错网格之中涉及到多达7种网格，网格越多，在复杂各向异性介质中有限差分数值模拟地震波时复杂度就越大。

b、旋转交错网格

Saenger等首次提出了有限差分旋转交错网格的方法。旋转交错网格的主要特点是：在求解各波场分量的一阶空间导数时，通过组合沿图2中四个对角线各个顶点波场在四个对角线方向的导数，并且将其映射到某个坐标方向达到计算一阶空间导数的目的。图2中旋转交错网格中只涉及到两种网格，与标准交错网格相中多达七种网格的计算量相比较，旋转交错网格方法实现简单，而且这种网格交错的方式能保证在求解任意各向异性介质（包括单斜和三斜介质）不用进行波场插值，因而精度高。虽然旋转交错网格有居多优点，但是，从图2中不难看出，因为通过组合求解方式求导无疑增加了求取每个空间导数的计算量，从而增加总的计算量；另外，在相当的频散要求的条件下，沿对角线求取导数方式要求适当缩小计算网格的大小，这样一来，在模拟计算时，网格的数量会增多，内存消耗增大。因此，旋转交错网格的特点是：方格种类只有2种，在复杂各向异性介质中有限差分数值模拟地震波时精度高，适合任意各向异性中波动方程求解，但是计算量大（10倍于标准交错网格），占用内存多（5.2倍于标准交错网格）。

c、菱形交错网格

上面讨论的两种网格一个共同的特点是每种网格都是矩形。但是，在菱形交错网格中，包括两个网格，每个网格都是菱形，如图3所示。Lisitsa和Vishnevskiy首次将菱形交错网格用于各向异性介质的地震波场的模拟。菱形交错网格的好处是不需要旋转交错网格中的每种组合即可求出一阶空间导数，而且只有两种网格。可处理任意各向异性介质中的波动方程传播问题，无需插值，总体上由于菱形网格的固有特点，计算量比旋转交错网格少（4倍于标准交错网格），需要的计算内存比旋转交错网格少（4倍于标准交错网格）。另一方面，菱形网格本身引出的问题就是网格是菱形的，计算机无法处理菱形网格，只能想法将其转化成矩形网格。

Lisitsa和Vishnevskiy（2010）给出了菱形交错网格的实现结果，但是并没有明确提出这种网格是如何计算机实现的。故，如何将菱形网格转化为矩形网格，且在复杂各向异性介质中有限差分数值模拟地震波时不失菱形交错网格的优点这一难题继续得到解决。

发明内容

为解决上述问题，本发明提出一种在复杂各向异性介质中模拟地震波的方法及装置，将菱形网格转化为矩形网格，又不失菱形交错网格的优点，便于算法理解和计算机实现。

为实现上述目的，本发明提供了一种在复杂各向异性介质中模拟地震波的方法，所述方法包括：

构造三维地质模型，给模型的每个网格点的弹性常数c和密度ρ赋值；

通过获取三维地质模型中在网格点坐标为（i，j，k）的沿三个坐标轴的速度分量的一阶空间导数；其中，当网格点的坐标i+j+k=偶数，T_i,j,k=v；当网格点的坐标i+j+k=奇数，T_i,j,k的速度分量为网格点（i，j，k）的上、下、左、右、前、后的六个坐标值为偶数的网格点的速度分量的平均值；

将获取的所述速度分量v的一阶空间导数代入到弹性介质中弹性波方程组中的第一方程式，并对第一方程式进行求解获取应力分量σ；其中，应力分量σ为速度分量v的表达式；

利用所述应力分量σ，根据 $T_{i,j,k}=\left\{\begin{array}{cc}\mathrm{\σor\; c\; or\ρ},& (i+j+k)\%2=1\\ v,& (i+j+k)\%2=0\end{array}\right.,$ 通过获取三维地质模型中在网格点坐标为（i，j，k）的沿三个坐标轴的应力分量σ的一阶空间导数；

将获取的所述应力分量σ的一阶空间导数代入到弹性介质中弹性波方程组中的第二方程式，并对第二方程式进行求解获取速度分量v；

根据通过第二方程式获取的速度分量得到当前时刻的速度分量的值，作为当前时刻的地震波的记录；

获取下一时刻的地震波的记录，直至获取时间长度t的地震波的波场值。

可选的，在本发明一实施例中，所述弹性介质中弹性波方程组中的第一方程式为：

$\left\{\begin{array}{c}{\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{xx}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{yy}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{zz}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{yz}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{xz}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{xy}}\end{array}\right\}=\left\{\begin{array}{cccccc}{c}_{11}& c_{12}& c_{13}& c_{14}& c_{15}& c_{16}\\ c_{12}& c_{22}& c_{23}& c_{24}& c_{25}& c_{26}\\ c_{13}& c_{23}& c_{33}& c_{34}& c_{35}& c_{36}\\ c_{14}& c_{24}& c_{34}& c_{44}& c_{45}& c_{46}\\ c_{15}& c_{25}& c_{35}& c_{45}& c_{55}& c_{56}\\ c_{16}& c_{26}& c_{36}& c_{46}& c_{56}& c_{66}\end{array}\right\}\left\{\begin{array}{c}{v}_{x,x}\\ v_{y,y}\\ v_{z,z}\\ v_{y,z}+v_{z,y}\\ v_{z,x}+v_{x,z}\\ v_{y,x}+v_{x,y}\end{array}\right\}$

上式中，c为弹性常数，ρ为密度，σ表示波场的质点应力分量，其双下标表示应力的6个分量，其上面一点表示一阶时间导数；v表示波场的质点速度分量，其单下标表示应力的3个分量，其第二下表表示对x,y,z的一阶空间导数。

可选的，在本发明一实施例中，所述弹性介质中弹性波方程组中的第二方程式为：

$\left\{\begin{array}{c}\mathrm{\ρ}{\stackrel{\·}{v}}_x={\mathrm{\σ}}_{\mathrm{xx},x}+{\mathrm{\σ}}_{\mathrm{xy},y}+{\mathrm{\σ}}_{\mathrm{xz},z}\\ \mathrm{\ρ}{\stackrel{\·}{v}}_y={\mathrm{\σ}}_{\mathrm{xy},x}+{\mathrm{\σ}}_{\mathrm{yy},y}+{\mathrm{\σ}}_{\mathrm{yz},z}\\ \mathrm{\ρ}{\stackrel{\·}{v}}_z={\mathrm{\σ}}_{\mathrm{xz},x}+{\mathrm{\σ}}_{\mathrm{zy},y}+{\mathrm{\σ}}_{\mathrm{zz},z}\end{array}\right.$

上式中，ρ为密度，σ表示波场的质点应力分量，其双下标表示应力的6个分量，其第三下表表示对x,y,z的一阶空间导数；v表示波场的质点速度分量，其单下标x,y,z表示应力的3个分量，其上面一点表示一阶时间导数。

可选的，在本发明一实施例中，所述对第一方程式进行求解获取应力分量σ的方法为常规中心差分方法。

可选的，在本发明一实施例中，所述对第二方程式进行求解获取速度分量v的方法为常规中心差分方法。

为实现上述目的，本发明还提供了一种在复杂各向异性介质中模拟地震波的装置，所述装置包括：

三维地质模型获取单元，用于构造三维地质模型，给模型的每个网格点的弹性常数c和密度ρ赋值；

速度分量的一阶空间导数获取单元，用于通过获取三维地质模型中在网格点坐标为（i，j，k）的沿三个坐标轴的速度分量v的一阶空间导数；其中，当网格点的坐标i+j+k=偶数，T_i,j,k=v；当网格点的坐标i+j+k=奇数，T_i,j,k的速度分量为网格点（i，j，k）的上、下、左、右、前、后的六个坐标值为偶数的网格点的速度分量的平均值；

应力分量获取单元，用于将获取的所述速度分量v的一阶空间导数代入到弹性介质中弹性波方程组中的第一方程式，并对第一方程式进行求解获取应力分量σ；其中，应力分量σ为速度分量v的表达式；

应力分量的一阶空间导数获取单元，用于利用所述应力分量σ，根据 $T_{i,j,k}=\left\{\begin{array}{cc}\mathrm{\σor\; c\; or\ρ},& (i+j+k)\%2=1\\ v,& (i+j+k)\%2=0\end{array}\right.,$ 通过获取三维地质模型中在网格点坐标为（i，j，k）的沿三个坐标轴的应力分量σ的一阶空间导数；

速度分量获取单元，用于将获取的所述应力分量σ的一阶空间导数代入到弹性介质中弹性波方程组中的第二方程式，并对第二方程式进行求解获取速度分量v；

当前时刻地震波记录获取单元，用于根据通过第二方程式获取的速度分量得到当前时间第一时刻的速度分量的值，作为当前时刻的地震波的记录；

地震波的波场值获取单元，用于获取下一时刻的地震波的记录，直至获取时间长度t的地震波的波场值。

可选的，在本发明一实施例中，所述应力分量获取单元使用的弹性介质中弹性波方程组中的第一方程式为：

$\left\{\begin{array}{c}{\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{xx}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{yy}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{zz}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{yz}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{xz}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{xy}}\end{array}\right\}=\left\{\begin{array}{cccccc}{c}_{11}& c_{12}& c_{13}& c_{14}& c_{15}& c_{16}\\ c_{12}& c_{22}& c_{23}& c_{24}& c_{25}& c_{26}\\ c_{13}& c_{23}& c_{33}& c_{34}& c_{35}& c_{36}\\ c_{14}& c_{24}& c_{34}& c_{44}& c_{45}& c_{46}\\ c_{15}& c_{25}& c_{35}& c_{45}& c_{55}& c_{56}\\ c_{16}& c_{26}& c_{36}& c_{46}& c_{56}& c_{66}\end{array}\right\}\left\{\begin{array}{c}{v}_{x,x}\\ v_{y,y}\\ v_{z,z}\\ v_{y,z}+v_{z,y}\\ v_{z,x}+v_{x,z}\\ v_{y,x}+v_{x,y}\end{array}\right\}$

上式中，c为弹性常数，ρ为密度，σ表示波场的质点应力分量，其双下标表示应力的6个分量，其上面一点表示一阶时间导数；v表示波场的质点速度分量，其单下标表示应力的3个分量，其第二下表表示对x,y,z的一阶空间导数。

可选的，在本发明一实施例中，所述速度分量获取单元使用的弹性介质中弹性波方程组中的第二方程式为：

$\left\{\begin{array}{c}\mathrm{\ρ}{\stackrel{\·}{v}}_x={\mathrm{\σ}}_{\mathrm{xx},x}+{\mathrm{\σ}}_{\mathrm{xy},y}+{\mathrm{\σ}}_{\mathrm{xz},z}\\ \mathrm{\ρ}{\stackrel{\·}{v}}_y={\mathrm{\σ}}_{\mathrm{xy},x}+{\mathrm{\σ}}_{\mathrm{yy},y}+{\mathrm{\σ}}_{\mathrm{yz},z}\\ \mathrm{\ρ}{\stackrel{\·}{v}}_z={\mathrm{\σ}}_{\mathrm{xz},x}+{\mathrm{\σ}}_{\mathrm{zy},y}+{\mathrm{\σ}}_{\mathrm{zz},z}\end{array}\right.$

上式中，ρ为密度，σ表示波场的质点应力分量，其双下标表示应力的6个分量，其第三下表表示对x,y,z的一阶空间导数；v表示波场的质点速度分量，其单下标x,y,z表示应力的3个分量，其上面一点表示一阶时间导数。

可选的，在本发明一实施例中，所述应力分量获取单元对第一方程式进行求解获取应力分量σ的方法为常规中心差分方法。

可选的，在本发明一实施例中，所述速度分量获取单元对第二方程式进行求解获取速度分量v的方法为常规中心差分方法。

为实现上述目的，本发明提供了一种有限差分菱形交错网格优化方法，所述方法包括：

定义T(i,j,k)为在网格点（i,j,k）处的任何一个分量；并定义三维空间的坐标系z向下，x向右方，y向纸面外的方向；

根据 $T_{i,j,k}=\left\{\begin{array}{cc}\mathrm{\σor\; c\; or\ρ},& (i+j+k)\%2=1\\ v,& (i+j+k)\%2=0\end{array}\right.,$ 将菱形交错网格中两种菱形网格合并为一种网格，合并后的网格是规则的矩形网格，网格点（i,j,k）沿三个坐标轴的一阶空间导数为：

$\frac{T_{i+1,j,k}-T_{i-1,j,k}}{\mathrm{\Δx}},\frac{T_{i,j+1,k}-T_{i,j-1,k}}{\mathrm{\Δy}},\frac{T_{i,j,k+1}-T_{i,j,k-1}}{\mathrm{\Δz}}.$

上述技术方案具有如下有益效果：改进后的有限差分菱形交错网格将2种网格简化成了1种，将非常规菱形的网格简化成常规的矩形网格，提供了简化的条件，又不失菱形网格的基本特性，求解波场分量一阶空间导数极为方便，极大的简化了在复杂各向异性介质中有限差分数值模拟地震波的实现过程，便于理解和计算机实现，提高了程序设计效率和计算效率。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为标准交错网格三维示意图；

图2为旋转交错网格三维示意图；

图3为菱形交错网格三维示意图；

图4为改进型菱形交错网格示意图；

图5为本发明提供了一种在复杂各向异性介质中模拟地震波的方法流程图；

图6为本发明提供了一种在复杂各向异性介质中模拟地震波的装置框图；

图7a为本发明提供的模拟结果三维空间中速度x分量立体显示波场图；

图7b为本发明提供的模拟结果三维空间中速度y分量立体显示波场图；

图7c为本发明提供的模拟结果三维空间中速度z分量立体显示波场图；

图8为本发明提供的模拟结果二维平面显示不同分量的波场图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本申请的技术方案的工作原理：针对有限差分中菱形交错网格本身的局限，即网格是菱形的（不适合计算机处理），提出合并菱形交错网格中两种非常规菱形网格成一种矩形网格的改进方法。如图4所示，为改进型菱形交错网格示意图。在图4中，定义T(i,j,k)为在网格点（i,j,k）处的任何一个分量。定义三维空间的坐标系z向下，x向右方，y向纸面外的方向。假设i+j+k之和为奇数时，网格点的分量对应于应力分量、介质常量或介质密度；网格点坐标之和为偶数时，网格点的分量均为速度分量。依据这一规律，我们改进方法图3的菱形网格，将两种网格和为一种网格，如图4的网格示意图。

但是这种合并时有条件的，其条件是：

$T_{i,j,k}=\left\{\begin{array}{cc}\mathrm{\σor\; C\; or\ρ},& (i+j+k)\%2=1\\ v,& (i+j+k)\%2=0\end{array}\right.---\left(3\right)$

这样一来两种菱形的网格合并成为了一种网格，合并后的网格是规则的矩形网格，但是注意一点，这样做的目的是为了简化网格，便于计算和计算机实现，本质上还是菱形网格，并没有改变菱形网格的特性。好处是：在求取空间导数时，无需参考图3或图4，因为在求取任意一点的空间导数时，使用网格点位置的上下左右前后的六个点的值即可，从而使用的下标在合并后极为方便。例如：在空间坐标点为（i，j，k）的沿三个坐标轴的一阶空间导数为：

$\frac{T_{i+1,j,k}-T_{i-1,j,k}}{\mathrm{\Δx}},\frac{T_{i,j+1,k}-T_{i,j-1,k}}{\mathrm{\Δy}},\frac{T_{i,j,k+1}-T_{i,j,k-1}}{\mathrm{\Δz}}$

这里，T表示波场应力分量或速度分量，Δx，Δy，Δz表示x，y，z方向的网格间距，上式使用的是中心差分表示的方法。在实现标准交错网格和旋转交错网格时，一般必须仔细参考图1和图2，以便正确的找出对应的空间导数，而在使用改进后的菱形网格时，只需知道其网格下标之和的奇偶特点，就能顺利求出三个坐标方向的空间导数。

这里有两点需要注意的是：在模拟时，记录每个时间片的波场速度分量值的时候，在分量下标之和为偶数位置是速度分量（v），但是在最后我们输出数据时需要输出规则矩形网格的数据，所以，根据上式（3）可知，网格点坐标之和为奇数的位置没有速度分量的值，此时，可通过任何一个坐标点坐标之和为奇数的位置取上、下、左、右、前、后共六个坐标点坐标之和为偶数的位置的速度分量的平均值，这样最后输出的波场速度也是矩形形状，方便处理，这种插值和标准交错网格插值不同。

如图5所示，为本发明提供了一种在复杂各向异性介质中模拟地震波的方法流程图。所述方法包括：

步骤501）：构造三维地质模型，给模型的每个网格点的弹性常数c和密度ρ赋值

步骤502）：通过获取三维地质模型中在网格点坐标为（i，j，k）的沿三个坐标轴的速度分量的一阶空间导数；其中，当网格点的坐标i+j+k=偶数，T_i,j,k=v；当网格点的坐标i+j+k=奇数，T_i,j,k的速度分量为网格点（i，j，k）的上、下、左、右、前、后的六个坐标值为偶数的网格点的速度分量的平均值；

步骤503）：将获取的所述速度分量v的一阶空间导数代入到弹性介质中弹性波方程组中的第一方程式，并对第一方程式进行求解获取应力分量σ；其中，应力分量σ为速度分量v的表达式；

步骤504）：利用所述应力分量σ，根据 $T_{i,j,k}=\left\{\begin{array}{cc}\mathrm{\σor\; c\; or\ρ},& (i+j+k)\%2=1\\ v,& (i+j+k)\%2=0\end{array}\right.,$ 通过获取三维地质模型中在网格点坐标为（i，j，k）的沿三个坐标轴的应力分量σ的一阶空间导数；

步骤505）：将获取的所述应力分量σ的一阶空间导数代入到弹性介质中弹性波方程组中的第二方程式，并对第二方程式进行求解获取速度分量v；

步骤506）：根据通过第二方程式获取的速度分量得到当前时刻的速度分量的值，作为当前时刻的地震波的记录；

步骤507）：获取下一时刻的地震波的记录，直至获取时间长度t的地震波的波场值。

可选的，在本发明一实施例中，所述弹性介质中弹性波方程组中的第一方程式为：

$\left\{\begin{array}{c}{\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{xx}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{yy}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{zz}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{yz}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{xz}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{xy}}\end{array}\right\}=\left\{\begin{array}{cccccc}{c}_{11}& c_{12}& c_{13}& c_{14}& c_{15}& c_{16}\\ c_{12}& c_{22}& c_{23}& c_{24}& c_{25}& c_{26}\\ c_{13}& c_{23}& c_{33}& c_{34}& c_{35}& c_{36}\\ c_{14}& c_{24}& c_{34}& c_{44}& c_{45}& c_{46}\\ c_{15}& c_{25}& c_{35}& c_{45}& c_{55}& c_{56}\\ c_{16}& c_{26}& c_{36}& c_{46}& c_{56}& c_{66}\end{array}\right\}\left\{\begin{array}{c}{v}_{x,x}\\ v_{y,y}\\ v_{z,z}\\ v_{y,z}+v_{z,y}\\ v_{z,x}+v_{x,z}\\ v_{y,x}+v_{x,y}\end{array}\right\}$

上式中，c为弹性常数，ρ为密度，σ表示波场的质点应力分量，其双下标表示应力的6个分量，其上面一点表示一阶时间导数；v表示波场的质点速度分量，其单下标表示应力的3个分量，其第二下表表示对x,y,z的一阶空间导数。

可选的，在本发明一实施例中，所述弹性介质中弹性波方程组中的第二方程式为：

$\left\{\begin{array}{c}\mathrm{\ρ}{\stackrel{\·}{v}}_x={\mathrm{\σ}}_{\mathrm{xx},x}+{\mathrm{\σ}}_{\mathrm{xy},y}+{\mathrm{\σ}}_{\mathrm{xz},z}\\ \mathrm{\ρ}{\stackrel{\·}{v}}_y={\mathrm{\σ}}_{\mathrm{xy},x}+{\mathrm{\σ}}_{\mathrm{yy},y}+{\mathrm{\σ}}_{\mathrm{yz},z}\\ \mathrm{\ρ}{\stackrel{\·}{v}}_z={\mathrm{\σ}}_{\mathrm{xz},x}+{\mathrm{\σ}}_{\mathrm{zy},y}+{\mathrm{\σ}}_{\mathrm{zz},z}\end{array}\right.$

上式中，ρ为密度，σ表示波场的质点应力分量，其双下标表示应力的6个分量，其第三下表表示对x,y,z的一阶空间导数；v表示波场的质点速度分量，其单下标x,y,z表示应力的3个分量，其上面一点表示一阶时间导数。

可选的，在本发明一实施例中，所述对第一方程式进行求解获取应力分量σ的方法为常规中心差分方法。

可选的，在本发明一实施例中，所述对第二方程式进行求解获取速度分量v的方法为常规中心差分方法。

如图6所示，为本发明提供了一种在复杂各向异性介质中模拟地震波的装置框图。所述装置包括：

三维地质模型获取单元601，用于构造三维地质模型，给模型的每个网格点的弹性常数c和密度ρ赋值；

速度分量的一阶空间导数获取单元602，用于通过获取三维地质模型中在网格点坐标为（i，j，k）的沿三个坐标轴的速度分量v的一阶空间导数；其中，当网格点的坐标i+j+k=偶数，T_i,j,k=v；当网格点的坐标i+j+k=奇数，T_i,j,k的速度分量为网格点（i，j，k）的上、下、左、右、前、后的六个坐标值为偶数的网格点的速度分量的平均值；

应力分量获取单元603，用于将获取的所述速度分量v的一阶空间导数代入到弹性介质中弹性波方程组中的第一方程式，并对第一方程式进行求解获取应力分量σ；其中，应力分量σ为速度分量v的表达式；

应力分量的一阶空间导数获取单元604，用于利用所述应力分量σ，根据 $T_{i,j,k}=\left\{\begin{array}{cc}\mathrm{\σor\; c\; or\ρ},& (i+j+k)\%2=1\\ v,& (i+j+k)\%2=0\end{array}\right.,$ 通过获取三维地质模型中在网格点坐标为（i，j，k）的沿三个坐标轴的应力分量σ的一阶空间导数；

速度分量获取单元605，用于将获取的所述应力分量σ的一阶空间导数代入到弹性介质中弹性波方程组中的第二方程式，并对第二方程式进行求解获取速度分量v；

当前时刻地震波记录获取单元606，用于根据通过第二方程式获取的速度分量得到当前时间第一时刻的速度分量的值，作为当前时刻的地震波的记录；

地震波的波场值获取单元607，用于获取下一时刻的地震波的记录，直至获取时间长度t的地震波的波场值。

可选的，在本发明一实施例中，所述应力分量获取单元603使用的弹性介质中弹性波方程组中的第一方程式为：

$\left\{\begin{array}{c}{\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{xx}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{yy}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{zz}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{yz}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{xz}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{xy}}\end{array}\right\}=\left\{\begin{array}{cccccc}{c}_{11}& c_{12}& c_{13}& c_{14}& c_{15}& c_{16}\\ c_{12}& c_{22}& c_{23}& c_{24}& c_{25}& c_{26}\\ c_{13}& c_{23}& c_{33}& c_{34}& c_{35}& c_{36}\\ c_{14}& c_{24}& c_{34}& c_{44}& c_{45}& c_{46}\\ c_{15}& c_{25}& c_{35}& c_{45}& c_{55}& c_{56}\\ c_{16}& c_{26}& c_{36}& c_{46}& c_{56}& c_{66}\end{array}\right\}\left\{\begin{array}{c}{v}_{x,x}\\ v_{y,y}\\ v_{z,z}\\ v_{y,z}+v_{z,y}\\ v_{z,x}+v_{x,z}\\ v_{y,x}+v_{x,y}\end{array}\right\}$

上式中，c为弹性常数，ρ为密度，σ表示波场的质点应力分量，其双下标表示应力的6个分量，其上面一点表示一阶时间导数；v表示波场的质点速度分量，其单下标表示应力的3个分量，其第二下表表示对x,y,z的一阶空间导数。

可选的，在本发明一实施例中，所述速度分量获取单元605使用的弹性介质中弹性波方程组中的第二方程式为：

$\left\{\begin{array}{c}\mathrm{\ρ}{\stackrel{\·}{v}}_x={\mathrm{\σ}}_{\mathrm{xx},x}+{\mathrm{\σ}}_{\mathrm{xy},y}+{\mathrm{\σ}}_{\mathrm{xz},z}\\ \mathrm{\ρ}{\stackrel{\·}{v}}_y={\mathrm{\σ}}_{\mathrm{xy},x}+{\mathrm{\σ}}_{\mathrm{yy},y}+{\mathrm{\σ}}_{\mathrm{yz},z}\\ \mathrm{\ρ}{\stackrel{\·}{v}}_z={\mathrm{\σ}}_{\mathrm{xz},x}+{\mathrm{\σ}}_{\mathrm{zy},y}+{\mathrm{\σ}}_{\mathrm{zz},z}\end{array}\right.$

上式中，ρ为密度，σ表示波场的质点应力分量，其双下标表示应力的6个分量，其第三下表表示对x,y,z的一阶空间导数；v表示波场的质点速度分量，其单下标x,y,z表示应力的3个分量，其上面一点表示一阶时间导数。

可选的，在本发明一实施例中，所述应力分量获取单元603对第一方程式进行求解获取应力分量σ的方法为常规中心差分方法。

可选的，在本发明一实施例中，所述速度分量获取单元605对第二方程式进行求解获取速度分量v的方法为常规中心差分方法。

使用上面改进技术，实现在三维各向异性介质中地震波的模拟，模型大小是300x300x300,网个间距是4米，模拟的时间间隔是0.001秒，波源使用的是20Hz的雷克子波，震源在模型中心位置，介质的密度是2克每立方厘米，介质的36个弹性常数c为：

如图7a所示，为本发明提供的模拟结果三维空间中速度x分量立体显示波场图；如图7b所示，为本发明提供的模拟结果三维空间中速度y分量立体显示波场图；如图7c所示，为本发明提供的模拟结果三维空间中速度z分量立体显示波场图；如图8所示，为三维介质空间中三个速度分量在三个坐标平面中平面图。在上述四幅图中，记录时间为0.35秒，显示出三维介质空间的三个速度分量在各个平面的切片，可以看出改进后的有限差分菱形交错网格方法，能快速准确的模拟地震波在复杂各向异性介质中传播现象。

最后应说明的是：上述仅用以说明本发明而并非限制本发明所描述的技术方案；尽管本说明书对本发明已进行了详细的说明，但是，本领域的技术人员仍然可以对本发明进行修改或等同替换，一切不脱离本发明的精神和范围的技术方案及其改进，其均应涵盖在本发明的权利要求范围中。

Claims (11)

1.一种在复杂各向异性介质中模拟地震波的方法，其特征在于，所述方法包括：

给数值化的地质模型的每个网格点的弹性常数c和密度ρ赋值；

通过获取三维地质模型中在网格点坐标为（i，j，k）的沿三个坐标轴的速度分量的一阶空间导数；其中，当网格点的坐标i+j+k=偶数，T_i,j,k=v；当网格点的坐标i+j+k=奇数，T_i,j,k的速度分量为网格点（i，j，k）的上、下、左、右、前、后的六个坐标值为偶数的网格点的速度分量的平均值，我们称这种网格为菱形交错网格；

将获取的所述速度分量v的一阶空间导数代入到弹性介质中弹性波方程组中的第一方程式，并对第一方程式进行求解获取应力分量σ；其中，应力分量σ为速度分量v的表达式；

利用所述应力分量σ，根据 $T_{i,j,k}=\left\{\begin{array}{cc}\mathrm{\σor\; c\; or\ρ},& (i+j+k)\%2=1\\ v,& (i+j+k)\%2=0\end{array}\right.,$ 通过获取三维地质模型中在网格点坐标为（i，j，k）的沿三个坐标轴的应力分量σ的一阶空间导数；

将获取的所述应力分量σ的一阶空间导数代入到弹性介质中弹性波方程组中的第二方程式，并对第二方程式进行求解获取速度分量v；

根据通过第二方程式获取的速度分量得到当前时刻的速度分量的值，作为当前时刻的地震波的记录；

获取下一时刻的地震波的记录，直至获取时间长度t的地震波的波场值。

2.如权利要求1所述的方法，其特征在于，所述弹性介质中弹性波方程组中的第一方程式为：

$\left\{\begin{array}{c}{\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{xx}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{yy}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{zz}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{yz}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{xz}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{xy}}\end{array}\right\}=\left\{\begin{array}{cccccc}{c}_{11}& c_{12}& c_{13}& c_{14}& c_{15}& c_{16}\\ c_{12}& c_{22}& c_{23}& c_{24}& c_{25}& c_{26}\\ c_{13}& c_{23}& c_{33}& c_{34}& c_{35}& c_{36}\\ c_{14}& c_{24}& c_{34}& c_{44}& c_{45}& c_{46}\\ c_{15}& c_{25}& c_{35}& c_{45}& c_{55}& c_{56}\\ c_{16}& c_{26}& c_{36}& c_{46}& c_{56}& c_{66}\end{array}\right\}\left\{\begin{array}{c}{v}_{x,x}\\ v_{y,y}\\ v_{z,z}\\ v_{y,z}+v_{z,y}\\ v_{z,x}+v_{x,z}\\ v_{y,x}+v_{x,y}\end{array}\right\}---\left(1\right)$

上式中，c为弹性常数，ρ为密度，σ表示波场的质点应力分量，其双下标表示应力的6个分量，其上面一点表示一阶时间导数；v表示波场的质点速度分量，其单下标表示应力的3个分量，其第二下表表示对x,y,z的一阶空间导数。

3.如权利要求1所述的方法，其特征在于，所述弹性介质中弹性波方程组中的第二方程式为：

$\left\{\begin{array}{c}\mathrm{\ρ}{\stackrel{\·}{v}}_x={\mathrm{\σ}}_{\mathrm{xx},x}+{\mathrm{\σ}}_{\mathrm{xy},y}+{\mathrm{\σ}}_{\mathrm{xz},z}\\ \mathrm{\ρ}{\stackrel{\·}{v}}_y={\mathrm{\σ}}_{\mathrm{xy},x}+{\mathrm{\σ}}_{\mathrm{yy},y}+{\mathrm{\σ}}_{\mathrm{yz},z}\\ \mathrm{\ρ}{\stackrel{\·}{v}}_z={\mathrm{\σ}}_{\mathrm{xz},x}+{\mathrm{\σ}}_{\mathrm{zy},y}+{\mathrm{\σ}}_{\mathrm{zz},z}\end{array}\right.---\left(2\right)$

上式中，ρ为密度，σ表示波场的质点应力分量，其双下标表示应力的6个分量，其第三下表表示对x,y,z的一阶空间导数；v表示波场的质点速度分量，其单下标x,y,z表示应力的3个分量，其上面一点表示一阶时间导数。

4.如权利要求1所述的方法，其特征在于，所述对第一方程式进行求解获取应力分量σ的方法为常规中心差分方法。

5.如权利要求1所述的方法，其特征在于，所述对第二方程式进行求解获取速度分量v的方法为常规中心差分方法。

6.一种在复杂各向异性介质中模拟地震波的装置，其特征在于，所述装置包括：

三维地质模型获取单元，用于构造三维地质模型，给模型的每个网格点的弹性常数c和密度ρ赋值；

速度分量的一阶空间导数获取单元，用于通过获取三维地质模型中在网格点坐标为（i，j，k）的沿三个坐标轴的速度分量v的一阶空间导数；其中，当网格点的坐标i+j+k=偶数，T_i,j,k=v；当网格点的坐标i+j+k=奇数，T_i,j,k的速度分量为网格点（i，j，k）的上、下、左、右、前、后的六个坐标值为偶数的网格点的速度分量的平均值；

应力分量获取单元，用于将获取的所述速度分量v的一阶空间导数代入到弹性介质中弹性波方程组中的第一方程式，并对第一方程式进行求解获取应力分量σ；其中，应力分量σ为速度分量v的表达式；

应力分量的一阶空间导数获取单元，用于利用所述应力分量σ，根据 $T_{i,j,k}=\left\{\begin{array}{cc}\mathrm{\σor\; c\; or\ρ},& (i+j+k)\%2=1\\ v,& (i+j+k)\%2=0\end{array}\right.,$ 通过获取三维地质模型中在网格点坐标为（i，j，k）的沿三个坐标轴的应力分量σ的一阶空间导数；

速度分量获取单元，用于将获取的所述应力分量σ的一阶空间导数代入到弹性介质中弹性波方程组中的第二方程式，并对第二方程式进行求解获取速度分量v；

当前时刻地震波记录获取单元，用于根据通过第二方程式获取的速度分量得到当前时间第一时刻的速度分量的值，作为当前时刻的地震波的记录；

地震波的波场值获取单元，用于获取下一时刻的地震波的记录，直至获取时间长度t的地震波的波场值。

7.如权利要求6所述的装置，其特征在于，所述应力分量获取单元使用的弹性介质中弹性波方程组中的第一方程式为：

$\left\{\begin{array}{c}{\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{xx}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{yy}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{zz}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{yz}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{xz}}\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{xy}}\end{array}\right\}=\left\{\begin{array}{cccccc}{c}_{11}& c_{12}& c_{13}& c_{14}& c_{15}& c_{16}\\ c_{12}& c_{22}& c_{23}& c_{24}& c_{25}& c_{26}\\ c_{13}& c_{23}& c_{33}& c_{34}& c_{35}& c_{36}\\ c_{14}& c_{24}& c_{34}& c_{44}& c_{45}& c_{46}\\ c_{15}& c_{25}& c_{35}& c_{45}& c_{55}& c_{56}\\ c_{16}& c_{26}& c_{36}& c_{46}& c_{56}& c_{66}\end{array}\right\}\left\{\begin{array}{c}{v}_{x,x}\\ v_{y,y}\\ v_{z,z}\\ v_{y,z}+v_{z,y}\\ v_{z,x}+v_{x,z}\\ v_{y,x}+v_{x,y}\end{array}\right\}---\left(3\right)$

上式中，c为弹性常数，ρ为密度，σ表示波场的质点应力分量，其双下标表示应力的6个分量，其上面一点表示一阶时间导数；v表示波场的质点速度分量，其单下标表示应力的3个分量，其第二下表表示对x,y,z的一阶空间导数。

8.如权利要求6所述的装置，其特征在于，所述速度分量获取单元使用的弹性介质中弹性波方程组中的第二方程式为：

$\left\{\begin{array}{c}\mathrm{\ρ}{\stackrel{\·}{v}}_x={\mathrm{\σ}}_{\mathrm{xx},x}+{\mathrm{\σ}}_{\mathrm{xy},y}+{\mathrm{\σ}}_{\mathrm{xz},z}\\ \mathrm{\ρ}{\stackrel{\·}{v}}_y={\mathrm{\σ}}_{\mathrm{xy},x}+{\mathrm{\σ}}_{\mathrm{yy},y}+{\mathrm{\σ}}_{\mathrm{yz},z}\\ \mathrm{\ρ}{\stackrel{\·}{v}}_z={\mathrm{\σ}}_{\mathrm{xz},x}+{\mathrm{\σ}}_{\mathrm{zy},y}+{\mathrm{\σ}}_{\mathrm{zz},z}\end{array}\right.---\left(4\right)$

上式中，ρ为密度，σ表示波场的质点应力分量，其双下标表示应力的6个分量，其第三下表表示对x,y,z的一阶空间导数；v表示波场的质点速度分量，其单下标x,y,z表示应力的3个分量，其上面一点表示一阶时间导数。

9.如权利要求6所述的装置，其特征在于，所述应力分量获取单元对第一方程式进行求解获取应力分量σ的方法为常规中心差分方法。

10.如权利要求6所述的装置，其特征在于，所述速度分量获取单元对第二方程式进行求解获取速度分量v的方法为常规中心差分方法。

11.一种有限差分菱形交错网格优化方法，其特征在于，所述方法包括：

定义T(i,j,k)为在网格点（i,j,k）处的任何一个分量；并定义三维空间的坐标系z向下，x向右方，y向纸面外的方向；

根据 $T_{i,j,k}=\left\{\begin{array}{cc}\mathrm{\σor\; c\; or\ρ},& (i+j+k)\%2=1\\ v,& (i+j+k)\%2=0\end{array}\right.,$ 将菱形交错网格中两种菱形网格合并为一种网格，合并后的网格是规则的矩形网格，网格点（i,j,k）沿三个坐标轴的一阶空间导数为：

$\frac{T_{i+1,j,k}-T_{i-1,j,k}}{\mathrm{\Δx}},\frac{T_{i,j+1,k}-T_{i,j-1,k}}{\mathrm{\Δy}},\frac{T_{i,j,k+1}-T_{i,j,k-1}}{\mathrm{\Δz}}.$