改进的FCT方法在粘滞各向异性介质地震波模拟中的研究

李勤武， 周怀来

(成都理工大学 a.地球物理学院 b.油气藏地质及开发工程国家重点实验室

c.地球探测与信息技术教育部重点实验室，成都　610059)



改进的FCT方法在粘滞各向异性介质地震波模拟中的研究

李勤武a,b,c， 周怀来a,b,c

(成都理工大学 a.地球物理学院 b.油气藏地质及开发工程国家重点实验室

c.地球探测与信息技术教育部重点实验室，成都610059)

摘要：在地震勘探中，粘弹性各向异性介质相比弹性介质更能实际表征地下介质的性质，这里主要针对Kelvin模型从本构方程、几何方程和运动微分方程，推导其一阶应力速度方程，运用交错网格有限差分法对其进行数值模拟，研究该介质中地震波的波场特征与传播规律。同时分析了品质因子Q对地震波在振幅以及频率上的衰减、吸收作用的影响。在数值模拟中，使用交错网格算法必然会造成数值频散或假频现象，从而干扰数值模拟结果的正确性。为了提高数值模拟效果，对常规FCT(通量校正传输法)方法进行优化来压制频散，模拟实例证明，优化后的FCT方法要比常规FCT方法在消除频散方面更加有效，能有效地改善数值模拟精度。

关键词：Kelvin模型； 一阶应力速度方程； 数值模拟； 品质因子Q； 优化后FCT方法

0引言

1845年Stocks第一次对粘弹性介质做了深入的研究，从那以后关于粘弹性介质理论和应用不断发展。上世纪美国地球物理学家 N.H Rick[1]在Stokes 方程中引入地震结构子波理论，用统一完整的观点将地震波在地层中衰减包含在理论中。目前出现了许多有关描述粘弹性介质的数学模型，主要包括 Kelvin-Voigt模型、Maxwell模型、标准线性模型等[2]。在这些模型中Kelvin-Voigt模型的性质更接近于固相介质，称之为“固体型粘弹性介质”；而Maxwell模型其性质与流体相类似，也被称为“流体型粘弹性介质”。在国内,张剑锋等[3]在1994年对水平层状介质中粘弹性波进行了数值模拟；毕玉英等[4]实现了二维粘弹性介质中的波场正演模拟；孟凡顺等[5]通过有限差分算法对复杂地质体中的粘弹性波做了研究；杜启振等[6]研究了裂缝性地层粘弹性介质的波动方程；牛滨华等[2]在出版的《半空间介质与地震波传播》一书中详细研究了弹性介质、粘弹性介质以及双向介质中地震波的传播规律；唐启军[7]实现了对Von Karman型随机各向同性粘弹性单斜各向异性波动方程的数值模拟；贺同江等[8]将Forsyte广义正交多项式微分算子与褶积算子相结合，实现了粘弹性介质的数值模拟。

为了压制数值模拟中造成的数值频散， Boris等[9]首次把流体动力学中通量校正传输法引入到地球物理中；杨顶辉等[10]将通量校正传输法应用到弹性波动方程正演模拟中；董良国等[11]对高阶差分算法中各种因素对频散的影响进行过理论探讨；陈可洋[12]对通量校正传输法进行了改进。作者在Kelvin模型基础上，运用优化后的FCT方法来压制数值模拟中产生的频散效应，应用效果证明了该方法有效性和正确性。

1方法原理

在粘弹性各向异性介质研究中，开尔芬(Kelvin)介质是人们在描述地震波衰减时普遍应用的一种介质模型，其衰减形式更接近于固相衰减。这里针对开尔芬(Kelvin)介质模型进行数值模拟，并运用优化后的FCT方法进行频散压制。

1.1一阶应力速度方程

基于弹性力学理论，结合几何方程、本构方程和运动微分方程可以推导粘弹性介质的一阶应力速度方程。具体方程为式(1)。

(1)

其中：vx表示x方向速度；vy表示y方向速度；vz表示z方向速度；σxx表示x方向正应力；σyy表示y方向正应力；σzz表示z方向正应力；σxy表示xy方向切应力；σxz表示xz方向切应力；σyz表示yz方向切应力；ρ 表示介质密度。

1.2一阶应力速度方程离散化

在有限差分数值模拟中，相同的逼近阶数情况下，交错网格差分格式的精度大约是常规差分格式精度的4倍，通过提高逼近阶数会相应地提高差分精度，但同时也会降低运算速度。当空间差分逼近阶数提高到一定阶数后再提高逼近阶数，对精度的影响就不大了,因此在数值模拟时，需要兼顾计算速度和逼近精度。

对具有2N+1阶导数的连续函数f(x)离散化，其2N阶精度的交错网格差分格式为[14-17]式(2)。

(2)

其中：x是变量；Δx是步长；ci是差分权系数；i=0,1,2…,N。式(3)是式(1)离散化后时间2 阶空间2N阶方程

(3)

1.3频散压制方法

在求解一阶应力速度方程时，以差分形式代替微分形式对方程进行网格剖分，原本连续的介质和波场函数就被离散化，因而不可避免地引入了数值频散也称为伪波动，这种伪波动对研究地震波的波场特征和传播规律干扰非常大。因此对数值频散的有效压制方法研究，则尤为重要。

为了消除这种频散，常用的手段是引入FCT(通量校正传输)方法[8]进行压制，虽然常规的FCT可以比较好地压制频散，但是η1、η2(η1扩散通量校正参数，η2反扩散通量校正参数)两个参数不好选择，而且经过FCT后难免会留有一些未被压制的伪波动(尤其是η1、η2选择不当的时候，这种情况更加普遍)。这里对常规的FCT方法加以改进，进而提高频散压制效果。主要包括七个步骤：

1)通过离散化后的一阶应力速度方程，理论上可以计算出任意时刻的速度和应力值，以对速度分量vx做FCT频散压制为例，求出得到n和n+1时刻速度分量vx的值。

2)计算n时刻的弥散通量Qx、Qz：

(4)

其中，η1∈[0,1]可以通过考查地震道的振幅变化情况得到。

(5)

与η1一样η2的值也是在0～1之间，但是η2一般要比η1稍大。

4)利用弥散通量Qx、Qz对速度分量vx进行修正，这主要是对速度分量vx的波形行平滑处理，同时也会造成振幅损失。

(6)

(7)

6)通量反弥散处理，得到校正后的速度vx。

(8)

其中：

(9)

7)滤波处理。得到校正后的速度vx，选择一个矩形窗口作为vx i,k的实心邻域。得到实心邻域内所有的速度值后，计算所得它们的平均值用于代替vx i,k的值。

(10)

2数值模拟结果

本次研究建立了粘弹性各向异性介质模型，然后通过交错网格算法来实现介质中地震波的数值模拟，并运用改进的FCT方法来压制数值频散效应，以便于研究其波场特征和传播规律。

2.1单层模型

理论研究表明与高速层比较而言，对于低速层所能引起的频散就更为严重。所以，为了检验改进后的FCT方法频散压制效果，设计了一个单层低速模型。模型大小为2 000m×2 000m，空间采样步长Δx=Δz=10m，时间采样率为1ms，使用雷克子波作为震源子波，主频30HZ，震源坐标为(1 000m，1 000m)。其余参数设定见表1。

表1　单层模型介质参数

图1为粘弹性介质中对于不同参数(η1，η2)条件下单层低速模型在500 ms时的时间切片,为了便于对比且不失一般性，在此只对比水平方向速度分量的频散压制效果。

从图1中可以看出，当η1=0.0、η2=0.0此时相当于没做品散压制处理，时间切片图1(a)、图1(b)两者内部数值频散比较严重，信噪比很低，真实波场模糊不清；当η1=0.014、η2=0.017数值频散得到一定的压制，但是由于压制程度不足，图1(c)中的频散仍然很明显。不过与图1(c)相比，图1(d)使用优化后的FCT方法时间切片整体变得更加干净、清晰；随着参数的取值逐渐增大，频散压制效果得到进一步改善；当η1=0.035、η2=0.041时，图1(e)、图1(f)波场更加圆滑，信噪比得到提高，内部的伪波动得到较好地压制，也可以看出这时候的图1(f)效果要比图1(e)明显更好，对图1(f)而言，参数η1、η2已经快接近频散压制效果的最优参数。当η1=0.051、η2=0.060时，地震波的真实波场得到了很好地恢复，能量也得到了补偿，图1(g)、图1(h)的内部已经基本看不到频散、此时可以将η1、η2认为是最优参数。通过图1对比说明，η1、η2的选择对频散压制效果有直接影响，只有合适的η1、η2才能有效地对数值频散进行压制。模拟结果表明：在压制效果上优化后的FCT方法要比常规的FCT方法更好；在最优参数的选择上改进后FCT方法也比常规的FCT更具有适用性。

为了进一步对比优化后的FCT方法和常规的FCT方法，对图1中相应的时间切片抽取其第75道记录进行比较，图2(a)由于振动曲线跳动剧烈，从曲线上已经很难将真实波场所引的振动和频散造成的伪振动区分开来，这正好与时间切片图1(a)、图1(b)中严重的数值频散相吻合；图2(b)、图2(c)随着参数η1、η2取值增大，频散压制后的曲线变得平滑，大约在横坐标50～100处，此时可以看到，使用优化后的FCT方法得到的结果，要比常规FCT方法得到的振动曲线质点振动更平缓，说明优化后的FCT方法能更好地消除频散；图2(d)中两条曲线接近重合，可以看出完整的雷克子波波形、此时已经没有了频散带来的干扰振动。我们可以认为在相同的η1、η2取值时，优化后的FCT方法可以取得更好的频散压制效果，该方法比常规的FCT方法能更好地保持真实波场特性。

2.2三层模型

设计一个三层模型，进一步分析研究粘弹性介质中地震波的传播规律、振幅和频率衰减情况。当品质因子取值很大时，表明介质不再具有粘弹性而转化成了弹性介质。模型大小为2 000 m×2 000 m，空间采样步长Δx=Δz=10 m，时间采样率为1 ms，使用雷克子波作为震源子波，主频30 Hz。震源坐标为(1 000 m,100 m)，其余参数设定见表2。

表2　三层水平层状模型介质参数

图1　常规FCT与优化后FCT频散压制效果对比Fig．1　Comparison of the suppress dispersion effect of conventional 　　　FCT and optimized FCT(a)常规FCT ，(b)优化后FCT；(c)常规FCT ，(d)优化后FCT；(e)常规FCT ，(f)优化后FCT；(g)常规FCT ，(h)优化后FCT

图2　图1中相应时间切片抽取其第75道记录振幅对比Fig．2　The seventy-fifth record amplitude comparison of corresponding time slice in figure 1(a)η1=0.0,η2=0.0;(b) η1=0.014,η2=0.017;(c) η1=0.035,η2=0.041;(d) η1=0.051,η2=0.060

图3是三层模型共炮点地震记录，可以看到直达波和来地层反射界面的反射波。对比水平分量地震记录图3(a)、图3(b)，虽然两者的反射同向轴都可以比较清晰地看到，但是图3(b)中经过吸收衰减后同向轴能量比图3(a)的要弱、比较暗淡，大约在1.1 s时刻图3(a)、图3(b)它们的同向轴差异最明显。由于水平方向上炮点左右两边接收到的地震波起振方向不一致，因而出现了极性反转；而此时在垂直方向上炮点左右两边接收到的地震波起震方向都是一致的(与深度方向相反)，所以图3(c)、图3(d)中没有出现极性质反转现象。在垂直方向上的地震记录图3(c)、图3(d)，大约1 s后图3(c)中同向轴可以较为清晰地观察到，但是此时图3(d)中由于粘弹性介质对波场的吸收作用，使得其反射同向轴已经基本上不能识别。

图4是在前面图3中三层模型共炮点地震记录的基础上，分别抽取了弹性介质和粘弹性介质第75道地震记录切除掉直达波后做振幅能量对比(图3中黑线即为抽取的地震道位置)。图4从另一个角度揭示了粘弹性介质对地震波的吸收衰减作用，可以看出，水平方向图4(a)的振幅比垂直方向图4 (b)的稍强。图4中1.1 s后粘弹性记录振幅基本衰减零，这与图3中在1 s左右观察到的粘弹性介质中反射波能量衰减结果相吻合，同时也看出粘弹性介质不仅对地震波的振幅起到了衰减作用，也使得地震波的波形发生了畸变。

为了便于定量分析衰减所引起的能量变化，图5是对弹性与粘弹性介质中的模拟结果中第75道记录进行频谱分析。从图5中观察到，粘弹性记录的主频发生了变化，向小于主频30 Hz的方向移动。同时粘弹性记录的频谱低频部分有所增加而高频部分相应减少。

图3　三层模型共炮点地震记录Fig．3　The common shot seismic record of three-layer model(a)弹性水平分量；(b)粘弹性水平分量；(c)弹性垂直分量；(d)粘弹性垂直

图4　图3中第75道弹性记录与粘弹性记录振幅对比Fig．4　The elastic and viscoelastic seismogram amplitude contrast of the 　　　seventy-fifth record in figure 3(a)水平分量；(b)垂直分量

图5　图3中第75道弹性记录与粘弹性记录频谱对比Fig．5　The elastic and viscoelastic seismogram spectral contrast of the 　　　seventy-fifth record in figure 3(a)水平分量；(b)垂直分量

3结论

通过对粘弹性介质中地震波的数值模拟结果对比分析发现，在对一阶应力速度方程离散时不可必可避地带入了数值频散。作者提出优化后的FCT(通量校正传输)方法能够有效地解决数值频散问题，在频散压制参数η1、η2的选择上更具适用性，对于同样的η1、η2取值优化后的FCT方法在频散压制效果上也要比常规的FCT方法要好。

粘弹性介质对地震波的波场有吸收衰减作用，这种吸收衰减与品质因子有着一定联系。粘弹性介质除了会使得地震波的振幅能量发生衰减外，也会导致地震波的相位发生变化，并且主频也会随之衰减；在水平方向上地震记录有极性反转现象。

地震波数值模拟对解决和认识实际波场特征问题有重要的指导意义，同时也因为地下介质真实情况较为复杂，要彻底认识复杂粘弹性介质中地震波的传播规律还需要进一步研究。

参考文献：

[1]N．H．瑞克．粘弹性介质中的地震波[M]．许云，译．北京：地质出版社，1981．

N H RICK. Seismic wave in viscoelastic media[M]. XU Y，translation. Beijing: Geological Publishing House,1981.(In Chinese)

[2]牛滨华，孙春岩． 半无限空间各向同性粘弹性介质与地震波传播[M]． 北京: 地质出版社，2007．

NIU B H, SUN C Y. Semi-infinite space isotropic viscoelastic medium and seismic wave propagation[M].Beijing: Geological Publishing House, 2007.(In Chinese)

[3]张剑锋，李幼铭．水平层状介质中粘弹性波的计算[J]．计算结构力学及其应用，1994，11(3)：314-324．

ZHANG J F, LI Y M. Horizontally layered viscoelastic wave of computing[J]. Computational Structural Mechanics and its Applications, 1994，11(3)：314-324．(In Chinese)

[4]毕玉英，杨宝俊．二维粘弹性介质中完全波场正演模拟[J]．石油地球物理勘探，1995，30(3)：351-362．

BI Y Y, YANG B J. 2-D viscoelastic medium completely wave forward modeling [J]. oil geophysical prospecting, 1995，30(3)：351-362．(In Chinese)

[5]孟凡顺，郭海燕．复杂地质体粘滞弹性波正演模拟的有限差分法[J]．青岛海洋大学学报，2000，30(2)：315-320．

MENG F S, GUO H Y. Finite difference method for forward modeling of complex geologic body with viscous elastic wave[J]. Journal of Qingdao ocean university, 2000，30(2)：315-320．(In Chinese)

[6]杜启振，杨慧珠.裂缝性地层粘弹性地震多波波动方程[J].地球物理学报，2001(08)：2801-2806.

DU Q Z, YANG H Z. The fracture formation of Multiwave viscoelastic seismic wave equation[J]. Chinese Journal of Geophysics，2001(08)：2801-2806. (In Chinese)

[7]唐启军,韩立国,王恩利.基于随机各向同性背景的粘弹性单斜介质二维三分量正演模拟[J].西北地震学报,2009,31(1):35-39.

TANG Q J, HAN L G, WANG E L. Forward modeling of three component of viscoelastic medium based on random isotropic background[J].Journal of Northwest China,2009,31(1):35-39. (In Chinese)

[8]贺同江，刘红艳，李小凡.粘弹性介质地震波传播的褶积微分算子法数值模拟研究[J].西北地震学报,2010,32(4):318-324.

HE T J, LIU H Y, LI X F. Numerical simulation of the convolution differential operator method for seismic wave propagation in viscoelastic media[J].Journal of Northwest seismology,2010,32(4):318-324. (In Chinese)

[9]BOOK D L，BORIS J P，HAIN K．Flux-corrected trans-port,II: generalization of the method[J]．Computational physics，1975，18: 248-283．

[10]杨顶辉，藤吉文.各向异性介质中三分量地震记录的FCT 有限差分模拟[J]．石油地球物理勘探,1997,32( 2):181-190．

YANG D H,TENG J W. FCT finite difference simulation of three component seismic records in anisotropic media[J].oil geophysical prospecting,1997,32( 2):181-190．(In Chinese)

[11]董良国,李培明.地震波传播数值模拟中的频散问题[J].天然气工业, 2004, 24(6): 53-56.

DONG L G, LI P M. Frequency dispersion in numerical simulation of seismic wave propagation[J]. Gas Industry, 2004, 24 (6): 53 -56. (In Chinese)

[12]陈可洋. 地震波数值模拟中优化的通量校正传输方法[J].内陆地震,2012, 26(2): 169-179.

CHEN K Y, Optimal flux corrected transport method in numerical simulation of seismic wave[J].Inland earthquake,2012, 26(2): 169-179. (In Chinese)

[13]李录明,罗省贤.多波多分量地震勘探原理及数据处理方法[M].成都：成都科技大学出版社,1997.

LI L M, LUO S X. Multi component seismic exploration principle and data processing method[M].Chengdu:Chengdu University of Science and Technology Publishing House.1997. (In Chinese)

[14]韩翀.地震波数值模拟及波场特征分析[D].成都：成都理工大学,2007.

HAN C. Seismic wave numerical simulation and wave field characteristic analysis[D]. Chengdu:Chengdu University of Technology, 2007. (In Chinese)

[15]KINDELAN M，KAMEL A，SGUAZZERO P．On the construction and ef ciency of staggered numerical differentiators for the wave equation[J]．Geophysics，1990，55(1) : 107 -110．

[16]ROBERTSSON J O A，BLANCH J O，SYMES W W．Viscoelastic finite-difference modeling[J]． Geophysics，1994，59(9) : 1444 -1456．

[17]VIRIEUX J．P-SV wave propagation in heterogeneous Velocity-stress finite-difference method[J] ．Geophysics，1986，51: 889-901．

Study on seismic wave simulation of viscous anisotropic media using the improved FCT method

LI Qin-wua,b,c, ZHOU Huai-laia,b,c

(Chengdu University of Technolog a.College of Geophysics，b.State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation,c.Key Lab of Earth Exploration & Information Techniques of Ministry of Education，Chengdu610059，China)

Abstract:Compared to elastic media, the viscoelastic anisotropic media is more practical subsurface characterization in seismic exploration. In this paper mainly aimed at the Kelvin model, Derived the first order stress velocity equation from constitutive equations, geometric equations and differential equations. Then, the numerical simulation of the staggered grid finite difference is used to discrete it and research wave field characteristics and propagation of seismic waves. At the same time, the analysis of the influence of quality factor Q on the attenuation and absorption of seismic wave in amplitude and frequency. Use staggered grid algorithm would inevitably result in numerical dispersion or aliasing phenomenon in numerical simulation, and thus the correctness of the results of the numerical simulation is disturbed. In order to improve the results of numerical simulation, optimized the conventional FCT (flux corrected transport) method to suppress dispersion. Simulation examples demonstrate that the optimized FCT method is more effective than the conventional FCT method in eliminating frequency dispersion.

Key words:Kelvin model; first-order rate equation stress; numerical simulation; quality factor Q; the optimized FCT

中图分类号：P 631.4

文献标志码：A

DOI：10.3969/j.issn.1001-1749.2016.01.11

文章编号：1001-1749(2016)01-0074-09

作者简介：李勤武(1991-)，男，硕士，从事地震数据处理及波场正演模拟研究，E-mail：liqinwu04@163.com。

基金项目:国家自然科学基金(41204091)

收稿日期：2015-08-07改回日期：2015-09-08