CN110109177B - 基于旋转时空双变网格有限差分法的地震波正演模拟方法 - Google Patents

Abstract

本发明涉及一种基于旋转时空双变网格有限差分法的地震波正演模拟方法，在地震波模拟的过程中采用旋转双变网格算法，使该正演模拟具有更高的模拟精度与计算效率，同时提高网格剖分的灵活性，压制頻散。本发明需要细化的区域通过时间变步长和空间变网格距实现针对不同模拟情况采取不同的间距和步长，提升模拟精度以及计算效率。并利用全新的时间空间网格布局来处理网格的变化，空间网格中参数和变量分布与时空变化交界面避免了交错分布的情况，简化时间变步长和空间变网格这两个过程。该方法极大简化了空间变网格距和时间变步长算法，更易实现，同时不用对参数进行平均或者插值处理，提高模拟精度的同时避免了非均匀性较强的情况下不稳定性的产生。

Description

基于旋转时空双变网格有限差分法的地震波正演模拟方法

技术领域

本发明属于地震勘探技术领域，具体涉及一种时间域有限差分正演方法，特别涉及一种基于旋转时空双变网格有限差分法的地震波正演模拟方法。

背景技术

通过地震波数值模拟能够有效研究地震波的传播规律，同时在地震波反演过程中起着重要作用。有限差分法因为计算效率高和易于实现等优势成为波场正演模拟的常用方法。1968年，Alterman与Karal将有限差分法引入到地震正演模拟中。在1976年时Madariaga提出了一维速度-应力交错网格有限差分法。1984、1986年Virieux应用交错网格的差分格式对二维非均匀介质中的横波(SH波)以及转换横波(P-SV波)进行模拟。有限差分通常采用单一尺度的网格步长，但是对于大部分的模型如复杂的储层与构造和低速带等会导致局部采样不足影响模拟精度并带来很强的頻散或者局部过采样影响计算效率。Moczo首次利用空间可变网格进行数值模拟。Jastram和Behle提出了针对二维声波模型在某一深度采取变网格步长的算法。随后，Jastram和Tessemer将其推广到了二维弹性波的模拟中。Falk(1996,1998)将变网格步长算法用于井间地震中的管波模拟。Arben Pitarka对速度-应力方程采用不规则交错网格的方法，给出了四阶空间精度的变网格差分系数计算公式。当空间网格发生变化时，为了满足稳定性条件，小网格需要采用较小的时间步长，统一的时间步长会导致大网格区间过采样，大幅降低计算效率。Falk和Tessmer采用了局部时间采样变化的算法，但只允许2n倍的时间采样变化，因此并不适用于交错网格中。Tae-seob等采用插值的手段实现局部时间变化，解决了交错网格算法中时间维度上的过采样问题，时空双变网格得到了很好的应用和推广。

但是，因为交错网格中速度和应力交错分布在大小网格的交界处，因此在大网格波场值传递给小网格或者小网格对应的波场值反馈给大网格的时候显得十分麻烦。也是因为如此，交界处过度时间层中的速度和应力都需要插值求取。同时，因为插值需要大网格中相应的前后两个时刻的波场值，因此大网格进行一次迭代所得到的速度值和波场值并不能满足要求，还应该对过度区域中大网格节点再次迭代运算，得到下一时刻的值用于差值计算。这样的计算步骤显然十分麻烦，不仅编程实现比较麻烦，一定程度地增加了计算量，增加了误差。同时，当前所用的交错网格需要对物性参数进行插值或者平均，特别是非均匀性较强的介质时，可能会导致运算不稳定并产生数值噪音，遇到自由界面还需要进行单独的显式处理。

发明内容

本发明的目的就在于针对上述现有技术的不足，提供一种基于旋转时空双变网格有限差分法的地震波正演模拟方法，该正演方法不用对参数平均或插值且双变处理过程更简单。

本发明的目是通过以下技术方案实现的：

本发明核心是在需要细化的区域通过时间变步长和空间变网格距的方法实现针对不同的模拟情况采取不同的间距和步长，以此提升模拟精度以及计算效率。在处理这两个步骤之前先将双变网格旋转45度角，在这个角度按照新的网格距和差分方向进行差分，此时正好能简化时间变步长和空间变网格这两个过程，同时在赋值参数的时候不同进行额外的平均和插值处理。本发明通过旋转的手段极大简化空间变网格距和时间变步长算法，使之更容易实现且避免了非均匀性较强的情况下不稳定性的产生。

一种基于旋转时空双变网格有限差分法的地震波正演模拟方法，包括以下步骤：

A、要对波动方程进行旋转处理，在旋转之后差分方向不再是沿x和z轴方向，而是与之呈45度夹角的新差分方向。以弹性波波动方程为例，旋转过程如下：

应力是τ_xx,τ_zx,τ_zz速度参数是v_x,v_z,，λ,μ是拉梅常数，ρ是密度。

B、在旋转交错网格中沿坐标轴的差分算子求取过程如下：

C、用上述公式进行整理得到：

式中

和

代表对角线方向，Δx和Δz表示沿坐标轴方向的空间步长，Δr是对角线的长度。

D、采用上述公式进行弹性波正演模拟，具体的旋转双变网格算法实现过程如下：

a、已知粗网格t_k时刻的速度初始值，先更新粗网格区域中的应力和速度值V^k+1和T^{k +3/2}。然后判断波场是否传递到细网格，若波场传递到细网格，进行变网格处理，这样可以节省计算空间，提升计算效率，也可以直接进行网格细化处理。

b、将粗网格中的点直接传递到细网格边界处，而细网格中没有的值采用插值的方法求取。插值公式为：

式中：nk表示的是粗细网格间距的倍数，F₁和F₂是粗网格点上并直接赋值给细网格上的点，f_i就是待差值求取的细网格的波场值。由差值公式，通过由T^k+1/2传递给细网格中的τ^k+3/6求得边界处其余的应力值。

c、利用降阶的方式更新计算细网格中v^k+4/6。

d、在常规双变网格中，由于缺少细网格边界处k+4/6时刻的v_x值，因此需要k和k+1时刻的V_x进行线性插值获得k+4/6时刻的v_x。而这里因为旋转之后v_x已全部求得，故直接更新k+5/6时刻的应力值。

e、k+5/6时刻的τ_xx,τ_zz和τ_xz需要插值求取，插值公式：

f、更新k+1时刻的速度值。同时这里需要将该时刻得到的细网格中的速度值反传给对应的粗网格，这里只有一种情况，两个速度分量不需要分开处理。

g、在常规交错双变网格中，在细网格中的速度值传递给粗网格之后，需要利用这些值计算T^k+3/2，再用T^k+3/2计算出过渡区域边界上的V^k+2，而旋转双变网格不需要再次去计算粗网格这些值。因此这里直接计算更新细网格中的τ^k+7/6。

h、插值求取k+7/6时刻边界上的τ_xx,τ_zz和τ_xz值。

i、更新细网格中k+8/6时刻下的速度值。

j、同样，对于速度值这里不再需要插值，直接计算k+3/2时刻的应力值。

k、将τ^k+3/2传递给的T^k+3/2并且不需要对三个应力分量分开处理。

以上就是旋转双变网格算法的一次迭代过程，循环此过程得到更高精度的模拟过程。

与现有技术相比，本发明的有益效果在于：本发明提出的旋转双变网格有限差分法在时间变步长和空间变网格距之前将双变网格旋转45度角，利用全新的网格布局来处理网格的变化，极大简化了整个算法实现过程，使之更容易实现且避免了非均匀性较强的情况下不稳定性的产生，提升模拟精度，降低存储空间，提升计算效率，提高网格剖分的灵活性，压制頻散。具体具有以下优势：

1、有效地避免了在非均匀性较强的介质中由于运算不稳定产生数值噪声。通常为了提升计算精度，压制頻散，采用交错网格的方式进行有限差分。此时需要对物性参数进行插值或者平均，这也是数值噪声产生的原因。而旋转双变网格将同一物理量定义在同一位置，不需要进行平均处理，避免了稳定性问题的出现。

2、能够有效地避免在复杂构造或者低速带处因为单一网格尺度引起的局部采样不足。因为地下介质在各部分的复杂程度不一，采用单一网格距过大则存在部分区域采样不足并在起伏界面产生頻散，网格距过小则会大量增加存储空间以及降低计算效率，因此旋转双变网格采用空间变网格距的方式只在复杂区域进行空间网格加密，提升计算效率的同时保证了模拟精度。

3、避免时间过采样，进一步提升计算效率。当空间网格变化之后，在保证计算稳定的前提下，需要按照小网格来取时间步长，这样就会在大网格区域产生时间过采样，大幅降低计算效率，而旋转双变网格在时间层也采取变步长策略，以此提升计算效率。

4、进一步简化时间变步长和空间变网格步骤，旋转双变网格通过旋转差分方向，将原来分布在四个不同位置的应力和速度分量从新调整到两个网格点上，并且在该差分方向上进行网格距变化处理和时间变步长处理是十分简洁的，处理的步骤流程只有原来的三分之二。

本发明基于旋转时空双变网格有限差分法的地震波正演模拟方法具有较好的灵活性，能够提高模拟精度的同时将计算效率同步提高，而且整个过程得到进一步简化，对于实际地下复杂介质的模拟效果有显著提高。

附图说明

图1空间网格变化示意图；

图2横坐标为xi时变时间步长示意图；

图3横坐标为xi+1/2时变时间步长示意图；

图4交错网格变量分布图；

图5旋转交错网格变量分布图；

图6旋转双变空间变网格示意图；

图7旋转双变时间变步长示意图；

图8-图9速度z分量在0.36s时波场快照；

图10-图11速度x分量在0.36s时波场快照；

图12-图13速度z分量的地震记录；

图14-图15速度z分量第300道地震记录；

图16向斜模型；

图17-图22不同算法下的波场快照以及地震记录。

具体实施方式

下面结合附图和实例对本发明进一步的详细说明：

本发明基于旋转时空双变网格有限差分法的地震波正演模拟方法是在MATLAB2014b平台上实现的。

一种基于旋转时空双变网格有限差分法的地震波正演模拟方法，包括以下步骤：

首先，要对波动方程进行旋转处理，在旋转之后差分方向不再是沿x和z轴方向，而是与之呈45度夹角的新差分方向。方程旋转过程如下：

应力是τ_xx,τ_zx,τ_zz速度参数是v_x,v_z,，λ,μ是拉梅常数，ρ是密度。

在旋转交错网格中沿坐标轴的差分算子求取过程如下：

用上述公式进行整理得到：

式中

和

代表对角线方向，Δx和Δz表示沿坐标轴方向的空间步长，Δr是对角线的长度，具体的参数分布如图4和图5所示，所有的速度分量和应力分量集中分布在两个不同的节点上。采用上述公式进行弹性波正演模拟，具体的旋转双变网格算法实现过程如下：

a、已知粗网格t_k时刻的速度初始值，先更新粗网格区域中的应力和速度值V^k+1和T^{k +3/2}。然后判断波场是否传递到细网格，若波场传递到细网格，进行变网格处理，这样可以节省计算空间，提升计算效率，也可以直接进行网格细化处理。

b、将粗网格中的应力值直接传递到细网格边界处，而细网格边界中没有的值，需要插值。如图1和图6所示，对于交界处无法直接传递的点中的变量需要用左右两侧的值进行插值，双变网格中空间插值存在两种不同的情况，而且在一个整网格中的变量插值需要附近两个网格中的变量，而旋转双变网格中只存在一种情况，且一个网格中的变量插值也只需要该网格中的赋值变量即可。插值公式为：

式中：nk表示的是粗细网格间距的倍数，F₁和F₂是粗网格点上并直接赋值给细网格上的点，f_i就是待差值求取的细网格的波场值。由差值公式，通过由T^k+1/2传递给细网格中的τ^k+3/6求得边界处其余的应力值。

c、利用降阶的方式更新计算细网格中v^k+4/6。

d、图2和3所示，在常规交错双变网格中，由于缺少细网格边界处k+4/6时刻的v_x值，因此需要k和k+1时刻的V_x进行线性插值获得k+4/6时刻的v_x，而且同样是两种情况，需要分开处理。而这里因为旋转之后v_x已全部求得，故直接更新k+5/6时刻的应力值，如图7所示只有一种时间变步长情况。

e、k+5/6时刻的τ_xx,τ_zz和τ_xz需要插值求取，插值公式：

f、更新k+1时刻的速度值。同时这里需要将该时刻得到的细网格中的速度值反传给对应的粗网格，这里只有一种情况，两个速度分量不需要分开处理。

g、在常规交错双变网格中，在细网格中的速度值传递给粗网格之后，需要利用这些值计算T^k+3/2，再用T^k+3/2计算出过渡区域边界上的V^k+2，而旋转双变网格不需要再次去计算粗网格这些值。因此这里直接计算更新细网格中的τ^k+7/6。

h、插值求取k+7/6时刻边界上的τ_xx,τ_zz和τ_xz值。

i、更新细网格中k+8/6时刻下的速度值。

j、同样，对于速度值这里不再需要插值，直接计算k+3/2时刻的应力值。

k、将τ^k+3/2传递给的T^k+3/2并且不需要对三个应力分量分开处理。

以上就是旋转双变网格算法的一次迭代过程，循环此过程得到更高精度的模拟过程。在交界面处，因为只对边界上的参数进行插值，在边界外的点都是差分计算求取的，因此只能采用二阶精度差分。当我们用更高阶精度模拟的时候，需要对细网格内部的值进行插值，也有另外一种方法就是在交界处采用降阶的方式，以八阶精度为例，依次在靠近边界处降为六阶、四阶、二阶，这样减少了存储空间与计算步骤，因此采用降阶法处理交界面，而每一个时间层的空间交界面都需要对细网格进行降阶处理。

实施例1

根据勘探要求，建立常见的层状模型，上层是低速层，因此采用细网格。为了更好地观察双变算法的虚假反射，故该低速层为水平层，纵波速度2200m/s，横波速度1090m/s，密度为1600Kg/m3；下层纵波速度3000m/s，横波速度1400m/s，密度为2000Kg/m3。加载震源是雷克子波，主频是25HZ。选取的全局时间步长为0.9ms，分别采用四种方法进行模拟。粗网格尺度大小x方向和z方向都是6m，细网格为2m，精细时间层时间步长为0.3ms。在阻抗界面位于地下380m处，用细网格将其覆盖，并且为了使虚假反射与反射波和转换横波区分开来，粗细网格交界面位于地下540m处，观察它们的模拟结果。

正演结果如图8-图15所示，图8是速度z分量v_z在0.36s时的波场快照，分别采用常规细网格、常规粗网格以及旋转双变网格算法。图8-图9中的1是反射波，2是转换横波，3是由于空间变网格与时间网格步长变化引起的数值反射。图10-图11是速度x分量v_x在0.36s时的波场快照，图10是v_z的地震记录，因为数值反射十分微弱，对上图进行了放大处理。图12-图13中的1、2、3、4分别对应直达波，反射波，数值反射和转换横波。图14-图15提取了第300道。可以看到双变交错网格具有与网格距为2的常规交错网格一样的模拟效果，同样旋转之后的双变网格也是如此，他们的曲线在直达波与反射波的位置是重合的，几乎看不到由于双变算法带来的虚假反射。我们用振幅值来衡量它们之间的强弱，交错网格与旋转双变网格的直达波振幅均为6.804×10-8而反射波的振幅均为3.681×10-9。旋转双变算法的数值反射是1.52×10-10，旋转双变算法的反射能量占直达波能量的0.2233％，占反射波的4.129％。旋转双变算法的引起的虚假反射十分微弱，可以忽略。与带来的可以忽略的数值反射相比，双变算法所带来的高精度与高效率使之成为正演模拟中的研究热点。以该模型为例，在相同的配置下，计算使波传播0.18s，常规细网格计算时间为550.376s，旋转双变网格耗时159.818。采用空间变网格时间变步长的方法是十分高效的，在这个模型中占时不到单一尺度网格的百分之三十。根据所计算的模型的不同，需要精细时间层与空间网格精度的区域大小不同，有时甚至远远小于十分之一，同时使有限差分算法更为灵活，而且从图14-图15中可以看到，双变算法带来的误差是十分微弱可以不计的，保证了算法的精度。这里也证明了旋转双变网格具有高效率与高精度的特点。

实施例2

在实施例1是水平界面，现在设计起伏地层，如图16所示，设计一个向斜模型。该模型中星形代表震源，三角形代表检波器接收点，红色实线是示意纵向上的网格尺度变化，密集的区域采用的是2m×2m的网格间距，0.3ms的时间步长，而稀疏的区域采用的是6m×6m的网格距，0.9ms的时间步长。图17-图22分别是0.45s时的波场快照以及地震记录。从图17-图22可以看到当介质中地层呈起伏状时，波场模拟出现了大量的頻散，因为有限差分方法采用的是矩形网格进行网格剖分，用矩形网格模拟起伏界面时是呈阶梯状，形成一个个突变点，当波传播到此处产生大量的人为绕射并形成頻散。为了压制頻散，通常采用更小的网格间距，但这样会加大计算量。综合这两幅图，可以看到旋转双变算法和与之对应的细网格算法都很好地压制了頻散，具有较高的精度，然而对于单一尺度的粗网格可以看到頻散十分严重。因此在某些复杂的地质结构为压制頻散或者提高精度，需要采用更细的网格间距，这也是需要发展时空双变网格算法的缘由。

Claims (1)

1.一种基于旋转时空双变网格有限差分法的地震波正演模拟方法，其特征在于，包括以下步骤：

A、对波动方程进行旋转，旋转后差分方向从沿x和z轴方向转换至呈45度夹角，方程旋转过程如下：

式中，τ_xx，τ_zx，τ_zz为应力，v_x，v_z，为速度，λ，μ为拉梅常数，ρ为密度；

B、在旋转交错网格中沿坐标轴求取差分算子：

C、利用公式(3)对公式(1)进行整理，得到公式(4)；

式中，

和

为对角线方向，Δx和Δz为沿坐标轴方向的空间步长，Δr为对角线的长度；

D、通过公式(4)进行弹性波正演模拟，其中，旋转双变网格有限差分过程为：

a、已知粗网格t_k时刻的速度初始值，先更新粗网格区域中的应力和速度值V^k+1和T^k+3/2；再判断波场是否传递到细网格，若波场传递到细网格，进行变网格或网格细化；

b、将粗网格中的应力值直接传递到细网格边界处，而细网格边界中没有的值，采用公式(5)插值；

式中：nk为粗细网格间距的倍数，F₁和F₂为粗网格点上并直接赋值给细网格上的点，f_i为待差值求取的细网格的波场值；

由公式(5)，通过由T^k+1/2传递给细网格中的τ^k+3/6求得边界处其余的应力值；

c、利用降阶的方式更新计算细网格中v^k+4/6；

d、更新k+5/6时刻的应力值；

e、通过插值公式(6)插值求取k+5/6时刻的τ_xx,τ_zz和τ_xz；

f、更新k+1时刻的速度值，并将该时刻得到的细网格中的速度值反传给对应的粗网格；

g、计算更新细网格中的τ^k+7/6；

h、插值求取k+7/6时刻边界上的τ_xx,τ_zz和τ_xz值；

i、更新细网格中k+8/6时刻下的速度值；

j、计算k+3/2时刻的应力值；

k、将τ^k+9/6传递给T^k+3/2；

循环步骤a-步骤k，得到更高精度的模拟过程。

Images