CN104732093B

CN104732093B - 一种基于弥散黏滞性波动方程的fct‑fdm正演模拟方法

Info

Publication number: CN104732093B
Application number: CN201510145190.0A
Authority: CN
Inventors: 高静怀; 赵海霞
Original assignee: Xian Jiaotong University
Current assignee: Xian Jiaotong University
Priority date: 2015-03-30
Filing date: 2015-03-30
Publication date: 2018-04-17
Anticipated expiration: 2035-03-30
Also published as: CN104732093A

Abstract

本发明公开了一种基于弥散黏滞性波动方程的FCT‑FDM正演模拟方法，从数值模拟的角度研究了弥散黏滞性波的传播特性。首先阐述了FCT‑FDM的计算步骤，设计一种流体饱和质模型利用该方法模拟弥散黏滞性波的传播，并将数值结果与声波模拟结果做对比分析。结果表明：弥散黏滞性波相对于声波在振幅上有明显的衰减，相位也会发生显著的变化。弥散黏滞性波的衰减程度主要取决于弥散衰减系数和黏滞性衰减系数的取值。本发明为地震波的数值模拟提供了一种有效的工具，尤其对于含有流体的介质，它能用于刻画地震波在流体介质中的衰减和相移特性。此外，该技术方案易于实现，可操作性强。

Description

一种基于弥散黏滞性波动方程的FCT-FDM正演模拟方法

技术领域

本发明属于地球物理勘探领域，涉及一种地震波正演模拟方法，特别涉及一种基于弥散黏滞性波动方程的FCT-FDM正演模拟方法。

背景技术

地震数值模拟是以地震波在地下介质中的传播理论为基础，不论是在勘探地震还是天然地震中，地震波的传播理论都得到了广泛的应用。地震波模拟技术是地球物理数据处理的基础，它是一种假设已知地下地质结构和相应物理参数的情况下，模拟地下地质结构中波传播规律并得到各观测点的地震记录的一种地震模拟技术。地震波模拟技术在地球物理勘探中有着重要的意义，是一种地震资料解释的强有力的辅助工具，为地震数据反演以及地震观测系统的设计和评价提供了理论基础和依据。随着计算机技术和地震波理论的迅速发展，出现了许多种地震数值模拟技术和方法。在各类地震模拟技术中有限差分法是出现最早的一种方法，它也是波动方程求解最常用的一种数值格式，由于其原理简单、容易实现和计算效率高的优点而被各个领域广泛应用。Alterman et al.首次利用有限差分法求解层状介质中弹性波的传播问题。随后基于Alterman和Karal的方法又发展了很多有限差分技术，如：Alford et al.对于声波方程的有限差分模拟进行了研究，并对其数值频散特性和差分格式的精度进行了分析；Kelly et al.利用均匀和非均匀有限差分格式模拟了弹性介质中波的传播，并给出了一种合成地震记录的方法；Virieux et al.将有限差分推广到交错网格中提出了一阶速度-应力方程交错网格有限差分法，交错网格差分算法不仅在弹性介质的波场模拟中得到了广泛的应用，还进一步被用于黏弹性介质和各向异性介质的波场正演中。这种方法消除了部分假象从而提高了地震模拟的精度，但当地震波在具有强反差界面的介质中传播时，比如裂缝、裂隙和孔隙等，交错网格有限差分法会产生不稳定现象。有很多学者研究了当介质参数存在强反差的界面(比如空气和岩石间的界面)时波的传播问题。为了处理由介质参数强反差界面导致的不稳定现象，Saenger et al.提出了一种旋转交错网格有限差分法，并将其应用于模拟各向同性介质中弹性波的传播。旋转交错网格可以处理介质中含有参数具有强反差的不连续界面，该方法不用对界面附近网格点上的弹性模量求平均，可以将其用于研究介质更加复杂情形中波的传播问题，比如含裂缝的介质。Saenger et al.还利用旋转交错网格有限差分法模拟了黏弹性介质和各向异性介质中波的传播。Krüger et al.利用旋转交错网格模拟裂缝处的散射波场，并将结果和解析解进行比较，两者相吻合，证明了旋转交错网格正演方法的正确性。何洋洋和高静怀等利用有限差分法对比研究了弹性介质和黏弹性介质中Rayleigh面波的频散特性。

然而，在常规有限差分法中，当单位波长内的网格点数太少时会导致严重的数值频散问题。Boris et al.在求解流体力学方程时提出了一种通量校正传输有限差分法(Flux Corrected Transport-Finite Difference Method—FCT-FDM)，并将它用于求解声波方程。FCT-FDM能有效地降低由大梯度的变化、间断等因素引起数值解的不稳定性，还可以消除粗网格条件下的数值频散。杨顶辉和滕吉文等将FCT技术与各向异性介质中求解波动方程的有限差分法结合，得到一种适合于声波和弹性波方程求解的通量校正传输有限差分法。

目前已经存在很多种描述地震数据所观察到的不同物理现象的理论模型，如弹性和黏弹性介质理论、Biot理论、喷射流理论等。在烃类储层中观察到的依赖频率变化的反射现象的物理机理尚未确定，Korneev et al.提出了一种弥散黏滞性模型，即在标量波动方程中增加了弥散衰减项和黏滞性阻尼项，并将实验观察的结果与弥散黏滞性模型的结果进行对比，发现低Q值(小于5)能解释实验和实际数据中所观察到的现象。另一方面，常规的扩散方程都没有考虑惯性流效应的影响，常规的波动方程也没有考虑扩散流效应的影响。由达西流描述的准静态扩散过程主要发生在低频频段，而由Biot-Gassmann公式描述的波动过程主要在高频频段起作用，这使得在其过渡频段会引起严重问题，因为在过渡频段存在扩散流和惯性流的耦合作用。因此，在描述多孔多相介质中波的传播机理时，应该考虑扩散流与惯性流的耦合作用，而弥散黏滞性波动方程考虑了扩散流和惯性流的耦合作用。

发明内容

本发明的目的在于提供一种基于弥散黏滞性波动方程的FCT-FDM正演模拟方法，将FCT技术与有限差分法结合得到一种数值模拟方法用于求解弥散黏滞性波动方程，进而从数值模拟的角度研究弥散黏滞性波的传播特性。

为了实现上述目的，本发明采用的技术方案包括以下步骤：

1)根据实际地质背景条件、岩石物理测试及其测井资料建立介质模型的参数；

2)对所设计的介质模型进行网格离散；

3)震源函数在空间上采用高斯函数，时间上采用Ricker子波，形式为：

s(x,z,t)＝g(x,z)·f(t)

式中：

f(t)＝(1-2(πf₀t)²)exp(-(πf₀t)²)

式中：f₀表示Ricker子波的中心频率，模型计算中f₀＝30Hz；β为常数；(x₀,z₀)表示震源的空间位置；

4)将二维弥散黏滞性波动方程中的微分用差分近似替代，得到其相应的有限差分格式；其中空间采样步长和时间采样步长必须满足以下数值格式的稳定性条件：

其中，h为空间网格步长；

5)对二维弥散黏滞性波动方程进行FCT-FDM计算，得到每一时刻的波场值，其具体方法如下：

首先，给出弥散黏滞性波动方程的数学描述，二维弥散黏滞性波动方程其形式为：

式中：u为波场函数；γ，η分别为弥散和黏滞性衰减系数，它们与岩石的孔隙度、渗透率以及流体的密度、粘度等有关；υ为非频散介质中波的传播速度；x，t分别为空间和时间变量；

式(1)中左端第一项表示惯性项，第二项为弥散耗损力即扩散项，第三项表示黏滞性阻尼，第四项为波动方程的弹性部分；

弥散黏滞性波动方程FCT-FDM方法的具体的计算步骤如下：

5-1)定义式(1)的数值解为给定初值为震源函数；其中，表示第n时间步在网格点(x_j,z_k)处的波场值；n为时间采样点，j为空间x方向的采样点，k为空间z方向的采样点；

5-2)利用弥散黏滞性波动方程有限差分格式计算

式中：h为空间网格步长；Δt为时间步长；

5-3)计算漫射通量：

a)计算第n时间步的漫射通量

其中，λ₁为参数；

b)利用漫射通量平滑差分方程的解

其中，为平滑以后的解；

5-4)计算反漫射通量：

a)计算第n+1时间步的漫射通量

式中：λ₂与λ₁类似，为参数；

b)利用和计算反漫射通量

c)利用反漫射通量修正平滑以后的解

式中：为修正以后的解，且

(11)式中“sign”为符号函数；(8)式中max(·)，min(·)分别表示物理量的最大值和最小值；

5-5)重复步骤5-2)～5-4)，直到最大的计算时间步，得到整个计算区域内不同时间的波场值；

6)在计算区域的边界处使用吸收边界条件并对其进行处理，模拟实际地下介质中波的传播。

所述步骤2)中，采用矩形网格对模型进行离散。

与现有技术相比，本发明具有以下有益效果：

本发明针对弥散黏滞性波动方程提出了一种地震波数值模拟方法，也即将FCT(通量校正传输)技术与FDM(有限差分方法)结合得到一种FCT-FDM方法。由于弥散黏滞性波动方程包含了弥散衰减项和黏滞性阻尼项，也即考虑了扩散项和惯性项的耦合作用，利用本发明所得方法求解该方程可用于研究含流体饱和介质中地震波的传播特性，刻画地震波在流体介质中的衰减和相移特性；该方法中引入了FCT技术使得在粗网格下也能得到较为精确的数值结果，也即引入FCT技术用于降低数值频散效应，进而提高数值模拟的精度。本发明所得方法的原理简单，算法易于实现且易于并行处理。

附图说明

图1本发明中流体饱和介质模型示意图；

图2弥散黏滞性波和声波在不同时刻的波场快照；其中，(a)、(c)、(e)、(g)、(i)和(k)表示弥散黏滞性波，(b)、(d)、(f)、(h)、(j)和(l)表示声波；

图3弥散黏滞性波和声波沿着含水饱和介质层垂向布置检波器接收到的地震记录对比图；其中，(a)、(b)、(c)和(d)分别为在空间位置(450m，396m)、(450m，420m)、(450m，456m)和(450m，510m)处接受到的地震记录；

图4弥散黏滞性波和声波沿着含油饱和介质层垂向布置检波器接收到的地震记录对比图；其中，(a)、(b)、(c)和(d)分别为在空间位置(150m，396m)、(150m，420m)，(150m，456m)和(150m，510m)处接受到的地震记录；

图5流体饱和介质模型的地面记录；其中，(a)为弥散黏滞性波，(b)为声波；

图6流体饱和介质模型的VSP记录；其中，(a)为弥散黏滞性波，(b)为声波。

具体实施方式

下面结合附图对本发明做进一步详细的说明

本发明的基础是地震波的传播及其数值模拟方法，本发明基于弥散黏滞性波动方程给出了一种FCT-FDM正演模拟方法，具体实施步骤分别为：

1)介质模型的建立

所建立的介质模型包含地质模型和数学物理模型，其中地质模型是地质人员依据实际地质背景条件建立，本说明的介质模型如图1所示；而数学物理模型本发明涉及弥散黏滞性波动方程，并结合实际地质背景、岩石物理实验、测井等资料给出所建立介质模型的参数，模型参数如表1所示。

表1流体饱和介质模型参数

2)模型离散化

对所设计的介质模型进行网格离散，本发明采用矩形网格对模型进行离散。

3)震源函数

本发明中震源函数在空间上采用高斯函数，时间上采用Ricker子波，形式为：

s(x,z,t)＝g(x,z)·f(t)

式中：

f(t)＝(1-2(πf₀t)²)exp(-(πf₀t)²)

式中：f₀表示Ricker子波的中心频率f₀＝30Hz，模型计算中；β为常数；(x₀,z₀)表示震源的空间位置。

4)弥散黏滞性波动方程的数值格式

将弥散黏滞性波动方程中的微分用差分近似替代，得到其相应的有限差分格式，如(2)式所示。其中空间采样步长和时间采样步长必须满足该数值格式的稳定性条件：

5)FCT-FDM计算

利用上面所描述的FCT-FDM计算步骤计算每一时刻的波场值。

弥散黏滞性波动方程FCT-FDM方法主要包含三个步骤：有限差分方程的迭代、差分方程解的平滑和反漫射通量校正。其中，有限差分方程的迭代使得地震波向外传播；平滑差分解的过程主要用于降低数值频散，进而提高数值解的精度；反漫射通量的校正主要是修正平滑以后方程的解。

式中：u为波场函数；γ，η分别为弥散和黏滞性衰减系数，它们与岩石的孔隙度、渗透率以及流体的密度、粘度等有关；υ为非频散介质中波的传播速度；x，t分别为空间和时间变量。

(1)式中左端第一项表示惯性项，第二项为弥散耗损力即扩散项，第三项表示黏滞性阻尼，第四项为波动方程的弹性部分。且方程(1)式中既包含了惯性项也包含了扩散项，也就是耦合了惯性流与扩散流，或者说耦合了低频扩散与高频波现象。

弥散黏滞性波动方程FCT-FDM方法的具体的计算步骤如下：

1.定义弥散黏滞性波动方程(1)的数值解为给定初值也即震源函数；

2.利用弥散黏滞性波动方程有限差分格式(2)计算

式中：h为空间网格步长；Δt为时间步长。

3.计算漫射通量：

a)计算第n时间步的漫射通量P，Q

其中，λ₁为参数，在数值模拟实验中，我们发现它的取值为0.008≤λ₁≤0.05可以达到好的结果。

b)利用漫反射通量P，Q平滑差分方程的解

其中，为平滑以后的解。

4.计算反漫射通量：

a)计算第n+1时间步的漫射通量

式中：λ₂与λ₁类似，两者的取值可以相同也可以不同。在下节的数值模拟中取λ₁＝λ₂＝0.01。

b)利用和计算反漫射通量X，Y

c)利用反漫射通量X，Y修正平滑以后的解

式中：为修正以后的解，且

(11)式中“sign”为符号函数；(8)式中max(·)，min(·)分别表示物理量的最大值和最小值。

5.返回第2步，重复计算步骤2—4，直到最大的计算时间步。这样就可以得到整个计算区域内不同时间的波场值。

6)吸收边界处理

实际的地下介质无限大，而数值模拟的计算区域是有限的，这样使得地震波传播到计算区域边界会产生人工边界反射波，因此，需要在计算区域的边界处使用合理的吸收边界条件并对其进行处理，才能模拟实际地下介质中波的传播。本发明利用赵海霞、高静怀等提出的弥散黏滞性波动方程吸收边界条件进行处理。数值模拟结果分析

本发明利用FCT-FDM对二维流体饱和介质模型分别用弥散黏滞性波动方程和声波方程进行了数值模拟，从波场快照和地震记录进行对比分析，阐述了弥散黏滞性波的衰减特性。二维流体饱和介质模型如图1所示，模型大小为[0，900m]×[0,900m]，两个均匀介质层中间为60m厚的含油、含水饱和介质层，层位于深度390m，其中含水饱和介质层位于薄层中间。介质参数如表1所示。震源位于(450m，9m)处，Ricker子波中心频率为f₀＝30Hz，空间采样间隔为dx＝dz＝3m，时间采样间隔为dt＝0.001s，时间延迟50ms，这些参数的取值满足弥散黏滞性波动方程数值格式的稳定性条件(14式)。

图2为不同时刻的波场快照，左边对应于弥散黏滞性波(如图2(a)，(c)，(e)，(g)，(i)，(k))，右边对应于声波(如图2(b)，(d)，(f)，(h)，(j)，(l))。图3为沿着含水饱和介质层垂向布置四个检波器接受的地震记录，检波器的位置分别为(450m，396m)、(450m，420m)，(450m，456m)和(450m，510m)，即检波器位于流体饱和介质层的中间和下方。图4为沿着含油饱和介质层垂向布置四个检波器接受的地震记录，检波器的位置分别为(150m，396m)、(150m，420m)，(150m，456m)和(150m，510m)。图3和图4中实线为弥散黏滞性波，虚线为声波。图5和图6分别为地面地震记录和VSP记录(检波器如图1所示)。

从图2明显地看出当弥散黏滞性波和声波未传到流体饱和介质层时，两者的波场是相同的，因为两者都在均匀介质中传播，如图2(a)和(b)所示。当波在流体饱和介质层中和它的下方传播时，弥散黏滞性波相对于声波其振幅有明显的衰减，如图2(a)-(l)所示。

从图3和图4可知，当弥散黏滞性波在流体饱和介质层中和穿过该层时振幅会产生明显的衰减，相位也有一定的变化，因为含水饱和介质的衰减参数比含油饱和介质的衰减参数大，所以沿含水饱和介质接受的地震记录的振幅衰减更强。在弥散黏滞性波情况下由于弥散和黏滞性衰减参数的影响，层间多次反射波被衰减的几乎不可见，而在声波情形下，层间多次反射波可以明显地看见。

通过将图5和图6所示的弥散黏滞性波和声波的地面记录、VSP记录进行对比，发现当波在流体饱和介质层中和穿过该层时，弥散黏滞性波动方程情况下的直达波和反射波都有很大程度的衰减，而声学介质中这些波都没有发生衰减。

由此可见，弥散衰减系数和黏滞性衰减系数对数值结果会产生很大的影响，对于较大的衰减参数会产生更强的衰减。而这两个衰减参数与介质的特性有关，它们与介质的孔隙度、渗透率、介质密度以及其它流体饱和岩石参数有关，实际应用时需要针对特定的地质模型来确定这两个参数。

以上内容仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明权利要求书的保护范围之内。

Claims

1.一种基于弥散黏滞性波动方程的FCT-FDM正演模拟方法，其特征在于，包括以下步骤：

2)对所设计的介质模型进行网格离散；

s(x,z,t)＝g(x,z)·f(t)

式中：

f(t)＝(1-2(πf₀t)²)exp(-(πf₀t)²)

<mrow> <mn>0</mn> <mo>&le;</mo> <mi>&Delta;</mi> <mi>t</mi> <mo>&le;</mo> <mi>m</mi> <mi>i</mi> <mi>n</mi> <mrow> <mo>(</mo> <mfrac> <msqrt> <mn>6</mn> </msqrt> <mn>4</mn> </mfrac> <mfrac> <mi>h</mi> <mi>&upsi;</mi> </mfrac> <mo>,</mo> <mfrac> <msup> <mi>h</mi> <mn>2</mn> </msup> <mrow> <msup> <mi>&gamma;h</mi> <mn>2</mn> </msup> <mo>+</mo> <mn>8</mn> <mi>&eta;</mi> </mrow> </mfrac> <mo>)</mo> </mrow> </mrow>

其中，h为空间网格步长；

<mrow> <mfrac> <mrow> <msup> <mo>&part;</mo> <mn>2</mn> </msup> <mi>u</mi> </mrow> <mrow> <mo>&part;</mo> <msup> <mi>t</mi> <mn>2</mn> </msup> </mrow> </mfrac> <mo>+</mo> <mi>&gamma;</mi> <mfrac> <mrow> <mo>&part;</mo> <mi>u</mi> </mrow> <mrow> <mo>&part;</mo> <mi>t</mi> </mrow> </mfrac> <mo>-</mo> <mi>&eta;</mi> <mrow> <mo>(</mo> <mfrac> <mrow> <msup> <mo>&part;</mo> <mn>3</mn> </msup> <mi>u</mi> </mrow> <mrow> <mo>&part;</mo> <msup> <mi>x</mi> <mn>2</mn> </msup> <mo>&part;</mo> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <mfrac> <mrow> <msup> <mo>&part;</mo> <mn>3</mn> </msup> <mi>u</mi> </mrow> <mrow> <mo>&part;</mo> <msup> <mi>z</mi> <mn>2</mn> </msup> <mo>&part;</mo> <mi>t</mi> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>-</mo> <msup> <mi>&upsi;</mi> <mn>2</mn> </msup> <mrow> <mo>(</mo> <mfrac> <mrow> <msup> <mo>&part;</mo> <mn>2</mn> </msup> <mi>u</mi> </mrow> <mrow> <mo>&part;</mo> <msup> <mi>x</mi> <mn>2</mn> </msup> </mrow> </mfrac> <mo>+</mo> <mfrac> <mrow> <msup> <mo>&part;</mo> <mn>2</mn> </msup> <mi>u</mi> </mrow> <mrow> <mo>&part;</mo> <msup> <mi>z</mi> <mn>2</mn> </msup> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>=</mo> <mn>0</mn> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>

弥散黏滞性波动方程FCT-FDM方法的具体的计算步骤如下：

5-2)利用弥散黏滞性波动方程有限差分格式计算

<mrow> <mtable> <mtr> <mtd> <mrow> <msubsup> <mi>u</mi> <mrow> <mi>j</mi> <mo>,</mo> <mi>k</mi> </mrow> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <mo>=</mo> <mo>&lsqb;</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>-</mo> <mn>4</mn> <mi>a</mi> <mo>-</mo> <mi>&gamma;</mi> <mo>(</mo> <mi>&Delta;</mi> <mi>t</mi> <mo>)</mo> </mrow> <mo>-</mo> <mn>4</mn> <mi>b</mi> <mo>&rsqb;</mo> <msubsup> <mi>u</mi> <mrow> <mi>j</mi> <mo>,</mo> <mi>k</mi> </mrow> <mi>n</mi> </msubsup> <mo>-</mo> <mo>&lsqb;</mo> <mn>1</mn> <mo>-</mo> <mi>&gamma;</mi> <mrow> <mo>(</mo> <mi>&Delta;</mi> <mi>t</mi> <mo>)</mo> </mrow> <mo>-</mo> <mn>4</mn> <mi>a</mi> <mo>&rsqb;</mo> <msubsup> <mi>u</mi> <mrow> <mi>j</mi> <mo>,</mo> <mi>k</mi> </mrow> <mrow> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </msubsup> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>-</mo> <mi>a</mi> <mrow> <mo>(</mo> <msubsup> <mi>u</mi> <mrow> <mi>j</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>k</mi> </mrow> <mrow> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </msubsup> <mo>+</mo> <msubsup> <mi>u</mi> <mrow> <mi>j</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <mi>k</mi> </mrow> <mrow> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </msubsup> <mo>+</mo> <msubsup> <mi>u</mi> <mrow> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> <mrow> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </msubsup> <mo>+</mo> <msubsup> <mi>u</mi> <mrow> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> <mrow> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>+</mo> <mrow> <mo>(</mo> <mi>a</mi> <mo>+</mo> <mi>b</mi> <mo>)</mo> </mrow> <mo>&lsqb;</mo> <mrow> <mo>(</mo> <msubsup> <mi>u</mi> <mrow> <mi>j</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>k</mi> </mrow> <mi>n</mi> </msubsup> <mo>+</mo> <msubsup> <mi>u</mi> <mrow> <mi>j</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <mi>k</mi> </mrow> <mi>n</mi> </msubsup> <mo>+</mo> <msubsup> <mi>u</mi> <mrow> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>n</mi> </msubsup> <mo>+</mo> <msubsup> <mi>u</mi> <mrow> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> <mi>n</mi> </msubsup> <mo>)</mo> </mrow> <mo>&rsqb;</mo> </mrow> </mtd> </mtr> </mtable> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow>

式中：h为空间网格步长；Δt为时间步长；

5-3)计算漫射通量：

a)计算第n时间步的漫射通量

<mrow> <mtable> <mtr> <mtd> <mrow> <msubsup> <mi>P</mi> <mrow> <mi>j</mi> <mo>,</mo> <mi>k</mi> </mrow> <mi>n</mi> </msubsup> <mo>=</mo> <msub> <mi>&lambda;</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <msubsup> <mi>u</mi> <mrow> <mi>j</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>k</mi> </mrow> <mi>n</mi> </msubsup> <mo>-</mo> <msubsup> <mi>u</mi> <mrow> <mi>j</mi> <mo>,</mo> <mi>k</mi> </mrow> <mi>n</mi> </msubsup> <mo>-</mo> <msubsup> <mi>u</mi> <mrow> <mi>j</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>k</mi> </mrow> <mrow> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </msubsup> <mo>+</mo> <msubsup> <mi>u</mi> <mrow> <mi>j</mi> <mo>,</mo> <mi>k</mi> </mrow> <mrow> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>,</mo> </mrow> </mtd> <mtd> <mrow> <mn>0</mn> <mo>&le;</mo> <msub> <mi>&lambda;</mi> <mn>1</mn> </msub> <mo>&le;</mo> <mn>1</mn> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msubsup> <mi>Q</mi> <mrow> <mi>j</mi> <mo>,</mo> <mi>k</mi> </mrow> <mi>n</mi> </msubsup> <mo>=</mo> <msub> <mi>&lambda;</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <msubsup> <mi>u</mi> <mrow> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>n</mi> </msubsup> <mo>-</mo> <msubsup> <mi>u</mi> <mrow> <mi>j</mi> <mo>,</mo> <mi>k</mi> </mrow> <mi>n</mi> </msubsup> <mo>-</mo> <msubsup> <mi>u</mi> <mrow> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> <mrow> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </msubsup> <mo>+</mo> <msubsup> <mi>u</mi> <mrow> <mi>j</mi> <mo>,</mo> <mi>k</mi> </mrow> <mrow> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>,</mo> </mrow> </mtd> <mtd> <mrow> <mn>0</mn> <mo>&le;</mo> <msub> <mi>&lambda;</mi> <mn>1</mn> </msub> <mo>&le;</mo> <mn>1</mn> </mrow> </mtd> </mtr> </mtable> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow>

其中，λ₁为参数；

b)利用漫射通量平滑差分方程的解

其中，为平滑以后的解；

5-4)计算反漫射通量：

a)计算第n+1时间步的漫射通量

<mrow> <mtable> <mtr> <mtd> <mrow> <msubsup> <mover> <mi>P</mi> <mo>~</mo> </mover> <mrow> <mi>j</mi> <mo>,</mo> <mi>k</mi> </mrow> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <mo>=</mo> <msub> <mi>&lambda;</mi> <mn>2</mn> </msub> <mrow> <mo>(</mo> <msubsup> <mi>u</mi> <mrow> <mi>j</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>k</mi> </mrow> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <mo>-</mo> <msubsup> <mi>u</mi> <mrow> <mi>j</mi> <mo>,</mo> <mi>k</mi> </mrow> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <mo>-</mo> <msubsup> <mi>u</mi> <mrow> <mi>j</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>k</mi> </mrow> <mi>n</mi> </msubsup> <mo>+</mo> <msubsup> <mi>u</mi> <mrow> <mi>j</mi> <mo>,</mo> <mi>k</mi> </mrow> <mi>n</mi> </msubsup> <mo>)</mo> </mrow> <mo>,</mo> </mrow> </mtd> <mtd> <mrow> <mn>0</mn> <mo>&le;</mo> <msub> <mi>&lambda;</mi> <mn>1</mn> </msub> <mo>&le;</mo> <mn>1</mn> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msubsup> <mover> <mi>Q</mi> <mo>~</mo> </mover> <mrow> <mi>j</mi> <mo>,</mo> <mi>k</mi> </mrow> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <mo>=</mo> <msub> <mi>&lambda;</mi> <mn>2</mn> </msub> <mrow> <mo>(</mo> <msubsup> <mi>u</mi> <mrow> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <mo>-</mo> <msubsup> <mi>u</mi> <mrow> <mi>j</mi> <mo>,</mo> <mi>k</mi> </mrow> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <mo>-</mo> <msubsup> <mi>u</mi> <mrow> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>n</mi> </msubsup> <mo>+</mo> <msubsup> <mi>u</mi> <mrow> <mi>j</mi> <mo>,</mo> <mi>k</mi> </mrow> <mi>n</mi> </msubsup> <mo>)</mo> </mrow> <mo>,</mo> </mrow> </mtd> <mtd> <mrow> <mn>0</mn> <mo>&le;</mo> <msub> <mi>&lambda;</mi> <mn>2</mn> </msub> <mo>&le;</mo> <mn>1</mn> </mrow> </mtd> </mtr> </mtable> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </mrow>

式中：λ₂与λ₁类似，为参数；

b)利用和计算反漫射通量

c)利用反漫射通量修正平滑以后的解

式中：为修正以后的解，且

<mrow> <mtable> <mtr> <mtd> <mrow> <msubsup> <mover> <mi>X</mi> <mo>~</mo> </mover> <mrow> <mi>j</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mo>,</mo> <mi>k</mi> </mrow> <mi>n</mi> </msubsup> <mo>=</mo> <mrow> <mo>(</mo> <msub> <mi>s</mi> <mi>x</mi> </msub> <mo>&CenterDot;</mo> <mi>max</mi> <mo>{</mo> <mn>0</mn> <mo>,</mo> <mi>min</mi> <mo>&lsqb;</mo> <msub> <mi>s</mi> <mi>x</mi> </msub> <mo>&CenterDot;</mo> <msubsup> <mi>X</mi> <mrow> <mi>j</mi> <mo>-</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mo>,</mo> <mi>k</mi> </mrow> <mi>n</mi> </msubsup> <mo>,</mo> <mo>|</mo> <msubsup> <mover> <mi>P</mi> <mo>~</mo> </mover> <mrow> <mi>j</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mo>,</mo> <mi>k</mi> </mrow> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <mo>|</mo> <mo>,</mo> <msub> <mi>s</mi> <mi>x</mi> </msub> <mo>&CenterDot;</mo> <msubsup> <mi>X</mi> <mrow> <mi>j</mi> <mo>+</mo> <mfrac> <mn>3</mn> <mn>2</mn> </mfrac> <mo>,</mo> <mi>k</mi> </mrow> <mi>n</mi> </msubsup> <mo>&rsqb;</mo> <mo>}</mo> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msubsup> <mover> <mi>Y</mi> <mo>~</mo> </mover> <mrow> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> </mrow> <mi>n</mi> </msubsup> <mo>=</mo> <mrow> <mo>(</mo> <msub> <mi>s</mi> <mi>y</mi> </msub> <mo>&CenterDot;</mo> <mi>max</mi> <mo>{</mo> <mn>0</mn> <mo>,</mo> <mi>min</mi> <mo>&lsqb;</mo> <msub> <mi>s</mi> <mi>y</mi> </msub> <mo>&CenterDot;</mo> <msubsup> <mi>Y</mi> <mrow> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>-</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> </mrow> <mi>n</mi> </msubsup> <mo>,</mo> <mo>|</mo> <msubsup> <mover> <mi>Q</mi> <mo>~</mo> </mover> <mrow> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> </mrow> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <mo>|</mo> <mo>,</mo> <msub> <mi>s</mi> <mi>y</mi> </msub> <mo>&CenterDot;</mo> <msubsup> <mi>Y</mi> <mrow> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>+</mo> <mfrac> <mn>3</mn> <mn>2</mn> </mfrac> </mrow> <mi>n</mi> </msubsup> <mo>&rsqb;</mo> <mo>}</mo> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> </mtable> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>8</mn> <mo>)</mo> </mrow> </mrow>

2.根据权利要求1所述的基于弥散黏滞性波动方程的FCT-FDM正演模拟方法，其特征在于：所述步骤2)中，采用矩形网格对模型进行离散。