CN107526105A

CN107526105A - 一种波场模拟交错网格有限差分方法

Info

Publication number: CN107526105A
Application number: CN201710676672.8A
Authority: CN
Inventors: 张懿洁; 高静怀; 彭济根; 柳清伙; 张金淼
Original assignee: Xiamen University; Xian Jiaotong University
Current assignee: Xiamen University; Xian Jiaotong University
Priority date: 2017-08-09
Filing date: 2017-08-09
Publication date: 2017-12-29

Abstract

本发明公开了一种波场模拟交错网格有限差分方法，该方法适用于任何介质的波动方程的波场模拟。本发明针对高阶交错网格有限差分方法中因差分系数求解产生的不稳定现象，给出基于范德蒙矩阵的算法求解差分系数有效克服了这一问题。该方法可提供精确稳定的波场模拟结果；该技术方案易于实现、适用性广泛。

Description

一种波场模拟交错网格有限差分方法

技术领域

本发明属于地震勘探技术领域，涉及数值模拟方法，尤其是一种精确稳定的波场模拟交错网格有限差分方法。

背景技术

波场正演模拟技术是地震波探测的重要基础，也是认识储层地震波响应的关键工具，更是地震资料解释中参数反演和成像的前提条件。由于有限差分方法实现简便、可灵活处理复杂介质，已被广泛应用于地震波正演模拟。有限差分方法是将波动方程中时间和空间导数用差商的形式表示，从而实现了时间和空间的离散。有限差分方法有很多种，包括显式和隐式有限差分方法，时间域和频率域有限差分方法，交错网格和旋转交错网格有限差分方法等等。

交错网格有限差分方法是指在波动方程的离散过程中，把介质参数和方程变量分别放在不同的网格点上，从而提高方法的准确性。为了进一步提高交错网格有限差分方法的准确性以满足实际需求，可使用高阶有限差分方法，或减小时间步长或空间网格大小。此外，还可以采用优化方法，例如模拟退火方法、最小二乘方法等来计算有限差分方法的差分系数已获得更好的准确性。但是，高阶有限差分方法是最方便的。

然而，当差分阶数足够高时，通过频散分析可发现有限差分方法不再稳定。原因是计算差分系数公式中的系数矩阵接近奇异，但仍采用矩阵求逆的方法来求解。目前，并未发现解决该问题的方法。

发明内容

本发明的目的在于克服上述现有技术的缺点，提供一种精确稳定的波场模拟交错网格有限差分方法，其针对高阶交错网格有限差分方法中因差分系数求解产生的不稳定现象，给出基于范德蒙矩阵的算法求解差分系数有效克服了这一问题。该方法可提供精确稳定的波场模拟结果；该技术方案易于实现、适用性广泛。

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

这种波场模拟交错网格有限差分方法，采用高阶交错网格有限差分方法进行数值模拟时，由矩阵求逆方法求解差分系数产生的不稳定现象，采用基于范德蒙矩阵的算法求解差分系数，得到稳定准确的波场结果。

进一步，以上波场模拟交错网格有限差分方法，具体包括以下步骤：

1)基于交错网格有限差分方法和Taylor展开方法，对一阶空间导数进行2M阶精度展开，通过对比等号两端变量的系数，得到2M阶精度空间导数的差分系数a_m求解公式：

2)将步骤1)中的公式(1)中线性方程组进行换元：

那么，公式(2)中的系数矩阵为范德蒙矩阵，写为

其中，x_m＝(2m-1)²。

进一步的，以上步骤2)中，为求解公式(2)中的差分系数a_m，采用以下算法：

首先，进行双重循环，外循环为k从1递增至M-1，内循环为n从M递减至k+1，操作b(n)＝b(n)-x(k)×b(n-1)；

接着，再进行双重循环，外循环为k从M-1递减至1，第一个内循环为n从k+1递增至M，操作b(n)＝b(n)/(x(n)-x(n-k))，第二个内循环为n从k递增至M-1，操作b(n)＝b(n)-b(n+1)；

算法中符号表示为x_M×1，其中x(m)＝(2m-1)²，b_M×1＝[1 0 … 0]；经过此算法处理，最终的b(n)就是公式(2)中线性方程组的变量。

进一步，以上步骤1)中，在交错网格有限差分方法中，2M阶精度空间导数如下所示：

上式中，h为空间网格大小，a_m(m＝1,2,...,M)是交错网格有限差分方法的差分系数；

对和进行Taylor展开，并对比等式左右两端函数f的n阶导数的系数，得到2M阶空间导数的差分系数a_m求解公式，即步骤1)的公式(1)。

本发明具有以下有益效果：

本发明的波场模拟交错网格有限差分方法适用于任何介质的波动方程的波场模拟，其针对高阶交错网格有限差分方法中因差分系数求解产生的不稳定现象，给出基于范德蒙矩阵的算法求解差分系数有效克服了这一问题。

进一步，本发明方便灵活，适用于任何介质波动方程的正演模拟，包括声学、弹性、粘弹性、孔隙弹性、各向异性介质。

进一步，本发明可有效克服矩阵求逆引起的不稳定现象，可计算准确稳定的差分系数，从而得到更稳定准确的波场信息。

附图说明

图1是不同差分阶数的频散曲线，其中τ＝1ms,h＝10m,v＝2000m/s,β＝kh，差分系数由矩阵求逆的方法求得；

图2是不同差分阶数的频散曲线，其中τ＝1ms,h＝10m,v＝2000m/s,β＝kh，差分系数由基于范德蒙矩阵的方法求得；

图3是由(a)矩阵求逆方法和(b)本发明计算的均匀弹性介质速度z分量在0.2s的波场快照；

图4是实际储层地球物理模型的介质参数示意图；

图5是由(a)矩阵求逆方法和(b)本发明计算的实际储层地球物理模型的速度z分量的地震记录。

具体实施方式

本发明针对采用高阶交错网格有限差分方法进行数值模拟时，由矩阵求逆方法求解差分系数产生的不稳定现象，采用基于范德蒙矩阵的算法求解差分系数解决了这一问题，可得到更为稳定准确的波场结果。

下面结合附图和实例对本发明进行详细的描述。

1)2M阶空间导数的差分系数a_m求解公式

在交错网格有限差分方法中，2M阶精度空间导数如下所示

这里，h为空间网格大小，a_m(m＝1,2,...,M)是交错网格有限差分方法的差分系数。

对和进行Taylor展开，并对比等式左右两端函数f的n阶导数的系数，可得到2M阶空间导数的差分系数a_m求解公式：

基于公式(5)，采用矩阵求逆的方法求取差分系数a_m，图1给出了不同差分阶数时的频散曲线，频散参数如下所示：

其中，

显然，当δ等于1时，没有数值频散，当δ大于或小于1时，有数值频散。从图中可看出，30阶交错网格有限差分方法的精度低于8阶交错网格有限差分方法。显然，这与“差分阶数越高，精度越高”的事实不符。

表1给出不同M时，公式(5)中系数矩阵的条件数：

表1 不同M时，公式(5)中系数矩阵的条件数

显然，该矩阵的条件数随着M的增大急剧增加，当M为10时，条件数足够大并且矩阵已接近奇异。因此，采用矩阵求逆的方法求取高阶交错网格有限差分方法的差分系数时，势必会产生不稳定现象。

2)差分系数求取公式的转换

经过代数处理，公式(5)可变为：

那么，上式的系数矩阵变为范德蒙矩阵：

其中，x_m＝(2m-1)²。

3)基于范德蒙矩阵求取差分系数的算法

为求解方程Vx＝b，可采用表2中的算法来实现。本算法也可表述为

接着，再进行双重循环，外循环为k从M-1递减至1，第一个内循环为n从k+1递增至M，操作b(n)＝b(n)/(x(n)-x(n-k))，第二个内循环为n从k递增至M-1，操作b(n)＝b(n)-b(n+1)。

算法中符号表示为x_M×1，其中x(m)＝(2m-1)²，b_M×1＝[10…0]。经过该算法的迭代，最终得到的b就是所要求的未知量x。

表2 方程Vx＝b的求解，其中V为范德蒙矩阵

图2给出了基于范德蒙矩阵的算法计算的频散曲线，可看出差分阶数越高，数值频散越小，即方法更加稳定准确。

稳定性分析

为了验证本发明的稳定性，通过是否满足稳定性条件说明。

在二维情况下，稳定性条件表达式如下所示

其中，τ为时间步长，v_max为模型的最大速度。对于声波介质，其为声波速度；弹性介质，其为P波速度；孔隙弹性介质，其为快P波速度。

表3给出了本发明是否满足上述稳定性条件，其中

表3 不同差分阶数情况下的稳定性条件，其中τ＝1ms,h＝10m,v＝2000m/s

数值结果

将本发明应用于均匀弹性介质及实际储层地球物理模型的数值模拟，以验证有效性。

均匀弹性介质参数：P波速度4000m/s，S波速度2000m/s，密度2600kg/m³。数值模拟时，时间步长为1ms，空间网格大小为10m，模型大小为[0,2000m]×[0,2000m]。震源为主频30Hz，时延1/30s的Ricker子波，位于模型中心。差分阶数为36，差分系数分别由矩阵求逆的方法及本发明方法求得。速度z分量的波场快照如图3所示。可看出，矩阵求逆方法的结果数值频散非常严重，而本发明计算的结果没有数值频散。

实际储层地球物理模型的P波、S波速度及密度如图4所示。改模型自上而下分别为：疏松砂岩、泥岩、致密砂岩、泥岩、致密砂岩、泥岩。其中第三层介质中包含倒梯形目标区，包含气层、油层和水层。数值模拟时，时间步长为1ms，空间网格大小为10m，模型大小为[0,58km]×[0,12.8km]。震源为主频10Hz，时延0.1s的Ricker子波，位于[29km,10m]。差分阶数为36，差分系数分别由矩阵求逆的方法及本发明方法求得。速度z分量的地震记录如图5所示。可看出，矩阵求逆方法的结果数值频散非常严重，而本发明计算的结果没有数值频散。说明本发明计算的差分系数更加准确稳定。

Claims

1.一种波场模拟交错网格有限差分方法，其特征在于，采用高阶交错网格有限差分方法进行数值模拟时，由矩阵求逆方法求解差分系数产生的不稳定现象，采用基于范德蒙矩阵的算法求解差分系数，得到稳定准确的波场结果。

2.根据权利要求1所述的波场模拟交错网格有限差分方法，其特征在于，

2)将步骤1)中的公式(1)中线性方程组进行换元：

那么，公式(2)中的系数矩阵为范德蒙矩阵，写为

其中，x_m＝(2m-1)²。

3.根据权利要求2所述的波场模拟交错网格有限差分方法，其特征在于，步骤2)中，为求解公式(2)中的差分系数a_m，采用以下算法：

算法中符号表示为x_M×1，其中x(m)＝(2m-1)²，b_M×1＝[10…0]；经过此算法处理，最终的b(n)就是公式(2)中线性方程组的变量。

4.根据权利要求2所述的波场模拟交错网格有限差分方法，其特征在于，步骤1)中，在交错网格有限差分方法中，2M阶精度空间导数如下所示：

<mrow> <mfrac> <mrow> <mo>&part;</mo> <mi>f</mi> </mrow> <mrow> <mo>&part;</mo> <mi>x</mi> </mrow> </mfrac> <mo>=</mo> <mfrac> <mn>1</mn> <mi>h</mi> </mfrac> <munderover> <mo>&Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </munderover> <msub> <mi>a</mi> <mi>m</mi> </msub> <mo>{</mo> <mi>f</mi> <mo>&lsqb;</mo> <mi>x</mi> <mo>+</mo> <mfrac> <mi>h</mi> <mn>2</mn> </mfrac> <mrow> <mo>(</mo> <mn>2</mn> <mi>m</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>&rsqb;</mo> <mo>-</mo> <mi>f</mi> <mo>&lsqb;</mo> <mi>x</mi> <mo>-</mo> <mfrac> <mi>h</mi> <mn>2</mn> </mfrac> <mrow> <mo>(</mo> <mn>2</mn> <mi>m</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>&rsqb;</mo> <mo>}</mo> <mo>+</mo> <mi>O</mi> <mrow> <mo>(</mo> <msup> <mi>h</mi> <mrow> <mn>2</mn> <mi>M</mi> </mrow> </msup> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow>