周成峰

(福建省建筑设计研究院有限公司 福建福州 350001)

0 引言

目前建立在经典弹性波理论基础之上的反射波法只适合研究固相、液相或气相等单相介质中地震波传播规律，而岩石裂隙发育的富水区及大部分土是由固相和流相(包括气相和液相)组成的双相或多相介质[1]。与单相介质不同，地震波在双相或多相介质中传播时有其特殊的传播规律及波动响应特征。Biot介质是一种典型的双相介质，高阶交错网格波场数值模拟技术可有效地用来研究Biot介质波动响应特征。

Biot介质高阶交错网格波场数值模拟是用有限的离散值来代替连续介质，离散代替连续在频率上总会存在计算误差，频率上的计算误差会导致地震波具有不同的相速度和群速度，即频散现象，其中数值频散包括空间频散和时间频散。数值频散严重干扰地震波场，极大影响数值模拟的精度和分辨率[2]。消除数值频散的方法和原则主要有：

(1)在兼顾计算量的前提下尽可能减少空间网格间距；

(2)当空间网格间距一定时，在保证空间分辨率的情况下尽可能降低子波主频；

(3)在保证稳定性的前提下，减小时间步长；

(4)适当提高时间和空间差分阶数，从而提高模拟精度，压制数值频散；

(5)对差分算子进行校正，补偿离散化引起的误差，但是利用差分算子校正进行压制数值频散会增加一定的计算量，在时间域递推显格式的差分算法有一定难度[3]；

(6)使用通量传输校正(FCT)方法平滑波场，压制数值频散[4]。

在使用FCT方法时会损失一定的能量和频率，并需要更多的内存量，但FCT方法可以用较低阶的差分方程就可达到较高精度的正演模拟，大大提高运算效率。同时，FCT方法可以放大时间和空间步长，从而抵消FCT带来的计算量增加，因此FCT是很好的一种消除数值频散方法[5]。

本文将在不考虑耗散的二维Biot各向同性介质中，运用Matlab编程对交错网格差分格式具有空间四阶精度和时间二阶精度的数值频散问题和FCT校正方法进行分析、探讨。

1 频散计算及FCT校正

1.1 空间与时间频散计算

为了考察网格离散对波动传播速度的影响，引入式(1)的无量纲物理量：

(1)

c为弹性波在真实连续介质中的传播速度；

物理量q可以反映数值速度与真实速度的关系。

为了方便公式推导和分析，引入无量纲参数ζ和H：

(2)

式中λ为波长；

Δt为时间步长；

Δx为x方向的空间步长；

Δz为z方向的空间步长。

(3)

其中：

Ai=cisin(πHicosγi)+c2sin(3πHicosγi)，i=x,z。

则γx、γz为平面波传播方向与坐标轴x、z所成夹角，c1、c2为空间四阶差分权系数。

当x和z方向的空间步长Δx=Δz=Δ时，引入控制数值频散的无量纲参数ζ=VP-solidΔt/Δ和网格尺寸无量纲参数H=Δ/λ，结合不考虑耗散的二维Biot各向同性介质中三类波的速度表达式及式(3)可推得规则网格快纵波数值频散qp-fast,慢纵波数值频散qp-slow和横波数值频散qs：

(4)

(5)

(6)

其中：

Ax=c1sin(πHcosγx)+c2sin(3πHcosγx)

Az=c1sin(πHcosγz)+c2sin(3πHcosγz)

式(4)～式(6)中，vshear-solid为骨架横波速度；

vp-solid为骨架纵波速度；

vf为孔隙流体纵波速度；

ρs、ρf分别为固体基质(固相)密度和孔隙流体(流相)密度；

密度比ξ定义为ρf/ρs；

φ为孔隙度；

弯曲度τ是独立于固体和流体密度的几何度量；

γ是与流体运动中骨架的微观模型有关的因子，当流体中的固体颗粒为球形时，γ=1/2。

1.2 FCT校正方法

FCT方法基本原理是假设所有的极值点都是由数值频散引起的，对所有网格点进行漫射校正处理，达到消除数值频散的效果。再对非局部极值点进行补偿的反漫射校正，使任何不需要漫射的地方得以恢复。实际上，漫射校正是一个线性的平滑过程，反漫射校正由于是有选择的，所以是一个非线性的过程[6]。对于交错网格有限差分速度-应力波动方程，应用FCT方法对差分数值解进行校正的过程中，可以只对速度进行计算，也可以只对应力进行校正，因为在完成对速度或应力的校正后，通过交错网格有限差分速度-应力波动方程对应力或速度进行了校正。显然，两种方法都需对4个参数进行校正，运算量是相同的，本文选择对速度进行校正。

利用FCT方法进行校正，广义地说主要有3步：交错网格差分计算、通量漫射校正(漫射因子η1)和通量反漫射校正(反漫射因子η2)，具体执行步骤详见参考文献[7]。

2 频散分析

在确保数值计算稳定的前提下，数值频散及FCT校正方法分析的基本流程如图1所示。

图1 数值频散及FCT校正方法分析基本流程

2.1 空间频散分析

应用Matlab编程，计算式(4)～式(6)可得，空间差分网格离散引起的快纵波(快P波)、慢纵波(慢P波)和横波的数值频散曲线，如图2～图4所示。其中，横坐标为单位波长内网格点数控制参数H，纵坐标为数值频散qp-fast、qp-slow和qs。在二维平面内考虑到Ax、Az关于角度γx、γz取值的对称性，考察3个不同的传播方向：沿x轴(γx=0°，γz=90°)，沿xz平面三分线(γx=30°，γz=60°)，沿xz平面对角线(γx=45°，γz=45°)。相关参数取值如表1所示。

表1 空间频散计算参数

其中，vp-solid、vshear-solid和vf的单位为m/s。

图2 空间四阶有限差分快P波数值频散曲线

图3 空间四阶有限差分慢P波数值频散曲线

图4 空间四阶有限差分横波数值频散曲线

分析图2～图4可知，无论是快纵波、慢纵波还是横波，因空间网格离散，数值速度随传播方向不同大于或小于真实速度，即发生频散；一般情况下，当00.2时，随着单位波长内网格点数的减少而数值频散现象越来越严重。

2.2 时间频散分析

应用Matlab编程，计算式(4)～式(6)可得，时间差分网格离散引起的快纵波、慢纵波和横波的数值频散曲线，如图5～图8所示，其中横坐标为时间步长Δt，纵坐标为数值频散qp-fast、qp-slow和qs。在二维平面内，考虑到Ax、Az关于角度γx、γz取值的对称性，考察3个不同的传播方向：沿x轴(γx=0°，γz=90° )，沿xz平面三分线(γx=30°，γz=60° )，沿xz平面对角线(γx=45°，γz=45°)。相关参数取值如表2所示。

表2 时间频散计算参数

其中，vp-solid、vshear-solid和的单位为m/s。

图5 时间二阶有限差分快P波数值频散曲线

图6 时间二阶有限差分慢P波数值频散曲线

图7 时间二阶有限差分横波数值频散曲线

分析图5～图7可知，无论是快纵波、慢纵波还是横波，因时间离散，数值速度均大于真实速度，即发生频散，但在不同传播方向上频散情况基本相同；一般情况下，当0<Δt≤0.2×10-3时，基本上不发生频散，当Δt>0.2×10-3时，随着时间步长的增加，数值频散现象越来越严重。基此可知，空间步长和时间步长越小，数值频散现象越弱，但相应地计算量会增加，计算效率降低，因此我们在选取空间步长和时间步长的值时，应在减弱数值频散的同时尽可能地减少计算量。

3 FCT校正分析及探讨

在如图8所示的Biot双相各向同性介质模型中(边界长度单位为m)，基于FCT校正公式详见文献[7]，应用Matlab编程，对比分析有、无使用FCT法压制数值频散的效果图，如图9～图12所示。

图8 FCT压制效果计算模型

其中，η1是漫射因子，它是常量或一个线性函数，其取值取决于有限差分阶段频散误差的大小。在实际应用中，漫射因子的值可以通过简单地质模型的数值模拟运算来确定，当漫射因子取值合理时，数值模拟的结果不随漫射因子的变化而明显变化。一般情况下，0.01≤η1≤0.05比较合适，模拟效果较好[8]。η2是反漫射因子，η2的取值与η1不同，这是因为振幅和分辨率的损失主要有两方面的原因，一是传统有限差分运算引起的频散，二是人为加入的漫射。反漫射运算不仅要补偿人为加入漫射引起的损失，也要补偿传统有限差分运算带来的振幅损失。因此，η2的取值应比η1大10%～15%[7]。

应用FCT校正方法，取η1=0.001,η2=0.0011，得到的频散压制效果图如图9所示。

(a)固相80ms波场快照

(b)固相检波点波形图

(c)流相80ms波场快照

(d)流相检波点波形图图9 频散压制效果图一

取η1=0.01，η2=0.011，得到的频散压制效果图如10所示。

(a)固相80ms波场快照

(b)固相检波点波形图

(c)流相80ms波场快照

(d)流相检波点波形图图10 频散压制效果图二

取η1=0.1，η2=0.11，得到的频散压制效果图如11所示。

(a)固相80ms波场快照

(b)固相检波点波形图

(c)流相80ms波场快照

(d)流相检波点波形图图11 频散压制效果图三

不使用FCT校正，得到的压制效果图如12所示。

分析图9～图12可知：当不使用FCT方法或η1、η2值相对设置较小时，波场模拟中数值频散现象较严重；当η1=0.01，η2=0.011时，波场模拟中数值频散得到很好压制；当η1、η2值相对设置较大时，漫射通量和反漫射通量校正过度，波场模拟中的有效波受到干扰，从而使得波场快照和检波点波形显得模糊。以上说明，可以用FCT法消除数值频散，但η1、η2值应控制在合理范围，取值过小则失去压制数值频散的效果；取值过大则真实波场会被平滑变得模糊不清，分辨率显著下降。本文建议η1取[0.01，0.05]，η2比η1大10%～15%。

(a)固相80ms波场快照

(b)固相检波点波形图

(c)流相80ms波场快照

(d)流相检波点波形图图12 频散压制效果图四

4 结论

在无耗散的二维Biot双相各向同性介质中，对具有空间四阶精度与时间二阶精度的交错网格差分格式的频散问题，及FCT校正方法，经Matlab软件编程计算分析，本文的结论有：

(1)空间和时间的网格离散都会带来数值频散。对于空间网格离散，一般情况下，当00.2时，随着单位波长内网格点数的减少而数值频散现象越来越严重。对于时间网格离散，一般情况下，当0<Δt≤0.2×10-2时，基本上不发生数值频散；当Δt>0.2×10-3时，随着时间步长的增加而数值频散现象越来越严重，但在不同传播方向上频散情况基本相同。

(2)空间步长和时间步长越小，数值频散现象越弱，但相应的计算量会增加，计算效率降低。为了进一步提高计算效率，减弱数值频散，根据能量守恒定律引进通量校正传输法(FCT)，对波场漫射通量和反漫射通量进行校正。经大量试算发现，当FCT方法的系数η1、η2取值过小时，会失去压制数值频散的效果；当η1、η2取值过大时，会把真实波场平滑变得模糊不清，分辨率显著下降；一般情况下，η1取[0.01，0.05]，η2比η1大10%至15%，压制数值频散效果较好。