地球物理学报  2011, Vol. 54 Issue (1): 200-207   PDF    
引用本文
李桂花, 冯建国, 朱光明. 黏弹性VTI介质频率空间域准P波正演模拟. 地球物理学报, 2011, 54(1): 200-207, doi: 10.3969/j.issn.0001-5733.2011.01.021.  
LI G H, FENG J G, ZHU G M. Quasi-P wave forward modeling in viscoelastic VTI media in frequency-space domain. Chinese J. Geophys. (in Chinese), 2011, 54(1): 200-207, doi: 10.3969/j.issn.0001-5733.2011.01.021.  
黏弹性VTI介质频率空间域准P波正演模拟
李桂花1, 冯建国1, 朱光明2     
1. 山东省沉积成矿作用与沉积矿产重点实验室,山东科技大学地质科学与工程学院, 青岛 266510;
2. 长安大学地质工程与测绘学院, 西安 710054
收稿日期 2009-01-06, 2011-01-05 收修定稿
基金项目: 国家自然科学基金项目(40904009,50974081)和山东科技大学科学研究'春蕾计划'项目联合资助
作者简介: 李桂花,女,1977年生,博士,现从事地震波传播与地球物理数字信号处理及数值模拟等的教学与研究工作.E-mail:lghxa@163.com
摘要: 有限差分方法是波场数值模拟的一个重要方法,时间域有限差分计算方法因按时间片递推计算,每个时间片的舍入误差会累积到下一片中,当时间片较多,最终会导致累积误差太大.而频率域计算是按频率片对空间网格进行整体求解方程组,其计算误差分配到了每个网格点上,并且各个频率片之间是独立计算的,因此不存在累计误差,而且在频率-空间域更易于模拟黏弹各向异性介质的衰减,同时采用25点优化差分算子克服了常规差分算子的数值频散.本文首先推导了黏弹性VTI介质的频率-空间域准P波波动方程,在此基础上,利用25点优化差分算子构造了差分格式,实现了频率空间域的正演模拟.并以断层模型为例,模拟并对比分析了利用弹性VTI介质波动方程差分格式和利用黏弹性VTI介质波动方程差分格式模拟的井间地震准P波波场记录,得出黏弹性介质地震记录由于衰减振幅相对较小,主频相对较低的结论.
关键词: 数值模拟      频率-空间域      各向异性      黏弹性      准P波     
Quasi-P wave forward modeling in viscoelastic VTI media in frequency-space domain
LI Gui-Hua1, FENG Jian-Guo1, ZHU Guang-Ming2     
1. Shandong Provincial Key Laboratory of Depositional Mineralization and Sedimentary Minerals, College of Geological Sciences and Engineering, Shandong University of Science and Technology, Qingdao 266510, China;
2. College of Geological Engineering and Geomatics of Chang'an University, Xi'an 710054 , China
Abstract: The finite difference is an important method of wavefield numerical simulation. It can calculate recursively by time slice, and the round-off errors of each time slice may accumulate to the next time slice in time domain, eventually, which will lead to too large accumulated error when the time slice is enough. However, solving all equation set in frequency domain is based on the frequency slice on the space grid, whose error is assigned to each grid point, it is without accumulated error and easier to simulate the attenuation of the viscoelastic anisotropic media in frequency-space domain. In order to overcome the numerical dispersion of routine finite-difference operators, 25-point optimized finite-difference operators are applied. This paper derives firstly the formula of quasi-P wave in viscoelastic VTI media in frequency-space domain, and then constructs the difference format using 25-point optimized finite difference operators, so it implements forward modeling. A case study of fault model compares and analyzes the crosswell seismic quasi-P wavefield stimulated with the finite difference method for both elastic VTI media and viscoelastic VTI media. In the end, it comes to the conclusion that the amplitude of seismic wave is relatively smaller and the main frequency is decreased in viscoelastic media.
Key words: Forward modeling      Frequency-space domain      Anisotropic      Viscoelastic      Quasi P-wave     
1 引言

许多学者已经对时间空间域黏弹各向异性介质数值模拟做了很多研究[1~3].除时间空间域正演模拟各向异性介质中的波场外，一些学者还研究了频率空间域正演模拟地震波场的方法.频率-空间域有限元解法最初由Lysmer和Drake[4]提出，其后该方法被Marfurt[5]、Marfurt和Shin等[6]所发展.Pratt和Worthington[7~9]将频率域有限差分模拟应用于井间层析反演和地震成像.Shin等[10~12]将频率域9点和25点有限差分应用到2D 标量波的正演模拟中.吴国忱等[13~16]和殷文等[17]将频率空间域有限差分模拟用于VTI介质弹性介质的模拟.他们的研究成果对频率空间域正演模拟都起到了很大的促进作用，但他们的研究还没有涉及黏弹性介质中波的正演模拟，为了能够模拟出更接近实际介质中的波，本文推导了频率空间域黏弹性VTI介质的波动方程及其25点差分格式，并模拟了黏弹介质模型井间地震的波场记录.

2 黏弹性VTI介质频率空间域准P波波动差分方程

弹性介质中的本构方程可写为

(1)

式中TI，EJ分别为应力、应变分量，CIJ为弹性系数.

对于黏弹性介质，由黏弹性理论，基于Kelvin- voigt模型的黏弹性波的本构方程可以由(1)式变为

(2)

其中，$\tilde{C}$IJ，I，J=1，2，…，6，为复弹性系数，它包含了黏弹性衰减系数ηIJ.

对于黏弹性VTI介质，CIJ和ηIJ构成如下矩阵:

(3)

其中，.

则由(2)、(3)式可推导出黏弹性VTI介质中的波动方程:

(4)

设平面波方程U=Aei(kxx+kyy+kzz-ωt) 代入(4)式，并忽略体力项，可得如下的Kelvin-Christoffel方程

(5)

式中

(6)

若要使(6)式有非零解，必须使该式的系数矩阵的行列式为零，即

(7)

通过求解方程式(7)，可获得黏弹性VTI介质纵横波耦合的频散关系方程.为了消除记录中强横波对纵波的干扰，根据Alkhalifah声波假设思想，假设横波速度为零，即VS0 = 0，则用Thomsen 参数表征的弹性参数值变为

(8)

式中ρ，VP0，ε，δ 由正演模型给出.

在二维情况下，ky=0时，将式(8)代入(6)式，得

(9)

因为衰减系数与品质因子具有如下关系:

(10)

式中，QP 为准P波的品质因子，ω 为圆频率.

将式(9)、(10)代入式(6)并整理得黏弹性VTI介质中准P波的频散方程:

(11)

两边乘以准P 波的傅里叶变换P(kx，kz，ω)，同时对方程两边进行关于kx，kz和ω 的反傅里叶变换，得到黏弹性VTI介质的时间-空间域准P 波波动方程:

(12)

对方程作傅里叶变换，将时间-空间域p(z，x，t)变换到频率-空间域F(z，x，ω)，可得二维黏弹性 VTI介质中准P波的频率-空间域波动方程:

(13)

为了书写方便，引入复数

则上式可简化为

(14)

3 黏弹性VTI介质频率空间域准P波差分格式

针对VTI介质频率域准P 波波动方程式(14)的有限差分解用25点差分算子加权平均，充分利用25点的信息有利于提高正演模拟的精度，抑制数值频散，其中差分系数又称为25 点优化差分算子，用高斯-牛顿优化方法进行求取.

图 1a说明代替∂2F/∂x2 的优化差分算子模型.每一行用5 个网格点来构成两个二阶中心差分算子，其中一个差分算子是用每一行的第1、3、5项(圆圈的项)，其间距为2Δx，另一个是用每一行的第2、3、4项(黑点的项)，其间距为Δx.每一行的这两个差分算子用加权系数c(黑点的)和d(圆圈的)来平均，作为这一行的差分算子，这样5行各有一个差分算子，再用加权系数b1，b2，b3 将其加权平均.假设第3行的加权系数为b1，第2 行和第4 行有相同的加权系数b2，第1行和第5行有相同的加权系数b3，这样就构成了25个网格点关于∂2F/∂x2 近似的优化差分算子.图 1b说明代替∂2F/∂z2 的优化差分算子模型，同理可计算关于∂2F/∂z2 的优化差分算子[12].将25点差分算子引入到∂4F/∂x2∂z2 的差分中，构造了两个差分算子，每个差分算子都是9项构成的四阶混合差分算子，一个间隔是Δx和Δz(黑点)，另一个间隔是2Δx和2Δz(圆圈)，再将这两个差分算子用加权系数e(黑点)和f(圆圈)来加权平均，其加权系数排列如图 1c所示，这样就构造了25点关于∂4F/∂x2∂z2 的近似优化差分算子.对于ω4F项中的波场值，把25个网格点内的波场值用加权系数a1，a2，a3，a4，a5，a6，加权平均后的值作为放置点Fi，j的波场值，并假设与放置点距离相同处网格点的加权系数相同，加权系数排列如图 1d 所示.∂2F/∂x2、∂2F/∂z2、∂4F/∂x2∂z2 及Fi，j的25点优化差分算子如下:

(15)

(16)

(17)

(18)

图 1 用于25点差分计算的网格点分布[12](a)∂2F/∂x2 ；(b)∂2F/∂z2 ；(c)∂2F/∂x2∂z2 ；(d)Fi，j Fig. 1 Computational grids used to approximate differential operators by the 25-point weighted-averaging scheme

将式(15)~(18)代入黏弹性VTI介质准P波频率空间域波动方程(14)，同时将VP0，ε，δ，QP 离散化，则离散化后的准P波的频率-空间域差分方程为

(19)

其中

式(19)是放置点Fi，j处采用优化算子的差分格式，对于其中每一个网格点Fi，j都能建立一个这样的方程，再加上震源项G，可建立如同(20)的方程组.

(20)

式中总网格点数为Nz·Nx(Nz，Nx分别为z，x方向的空间采样点数).式(20)是黏弹性VTI介质频率空间域准P波波动方程有限差分解的差分格式，求解上述方程组，可解得频率为ω 的单频波快照.对每一个频率ω 都求解这样一个方程组，就得到每一个频率下所有网格点的波场值，然后对每一个网格点所有频率下的波场值进行反傅里叶变换，就可得到该网格点上的波场值，而每一个时间点上所有网格点的波场值，就是时间波场快照，从而实现黏弹性 VTI介质的准P波的正演.

4 弹性和黏弹性VTI介质断层模型试算

断层模型如图 2所示，黏弹性介质模型参数见表 1.采用主频为50Hz的雷克子波作为震源，振幅系数为100，震源坐标为(0 m, 500 m).检波器位于右边x=1000m处的垂直线上，检波器间隔10m, 共100个检波器.在模拟中采用了特征值分析法吸收边界条件和完全匹配层法衰减边界条件.模拟所用的参数:空间采样间隔10m, 时间采样间隔2ms, 用频率-空间域25点有限差分方法模拟获得的弹性介质和黏弹性介质的井间地震记录见图 3(a, b)，从地震记录上肉眼看不出有什么差别.但从图 4 弹性和黏弹性记录的所有道的平均频率振幅图的对比图上可以看出，弹性介质主频50 Hz, 黏弹性介质主频约40Hz, 不仅主频向低频移动，而且振幅也从弹性介质的1000降低到黏弹性介质的600.

图 2 断层模型 Fig. 2 Fault model
图 3 弹性和黏弹性介质准P波井间地震记录的对比 (a)弹性介质；(b)黏弹性介质. Fig. 3 Comparision the Quasi-P wave records between the elastic media and the viscoelastic media (a) Elastic media；(b) Viscoelastic media.
图 4 弹性和黏弹性井间地震记录的频率振幅图的对比 (a)弹性介质；(b)黏弹性介质. Fig. 4 Comparison frequency-amplitude diagram of the crosswell data between the elastic media and the viscoelastic media (a) Elastic media；(b) Viscoelastic media.
表 1 黏弹性介质模型参数表 Table 1 Model parameters of viscoelastic media
5 结语

综上所述，频率空间域正演模拟中，衰减系数是频率的函数，可较容易实现黏弹性正演模拟.通过对弹性和黏弹性VTI介质中qP波井间地震正演模拟结果的对比，得出黏弹性介质地震记录由于衰减，振幅相对较小，高频吸收比低频多，主频相对较低，频带相对变窄.

参考文献
[1] 宋常瑜, 裴正林. 井间地震黏弹性波场特征的数值模拟研究. 石油物探 , 2006, 45(5): 508–513. Song C Y, Pei Z L. Numerical simulation of viscoelastic wavefield for crosshole seismic exploration. Geophysical Prospecting for Petroleum (in Chinese) , 2006, 45(5): 508-513.
[2] 郭智奇, 刘财, 杨宝俊, 等. 黏弹各向异性介质中地震波场模拟与特征. 地球物理学进展 , 2007, 22(3): 804–810. Guo Z Q, Liu C, Yang B J, et al. Seismic wave fields modeling and feature in viscoelastic anisotropic media. Progress in Geophysics (in Chinese) , 2007, 22(3): 804-810.
[3] 杜启振, 刘莲莲, 孙晶波. 各向异性黏弹性空隙介质地震波场伪谱法正演模拟. 物理学报 , 2007, 56(10): 6143–6148. Du Q Z, Liu L L, Sun J B. Numerical modeling of seismic wavefield in anisotropic viscoelastic porous medium with the pseudo-spectral method. Acta Physica Sinica (in Chinese) , 2007, 56(10): 6143-6148.
[4] Lysmer, Drake. A finite-element method for seismology. In: Bolt B A Eds. Methods in Computational Physics, Volume 11:Seismology: Surface Waves and Earth Oscillations. Academic Press Inc., 1972. 181~216
[5] Marfurt K J. Accuracy of finite-difference and finite-element modeling of scalar and elastic wave equations. Geophysics , 1984, 49(5): 533-549. DOI:10.1190/1.1441689
[6] Marfurt K J, Shin C S. The future of iterative modeling in geophysical exploration. In Eisner E, Ed. Handbook of Geophysical Exploration: I Seismic Exploration, Volume 21: Supercomputers in Seismic Exploration. Pergamon Press, 1989. 203~228
[7] Pratt R G, Worhington M H. Inverse theory applied to multi-source crosshole tomography, Part I:Acoustic wave-equation method. Geophys.Prosp. , 1990, 38: 287-310. DOI:10.1111/gpr.1990.38.issue-3
[8] Pratt R G. Inverse theory applied to multi-source crosshole tomography, part II:Elastic wave-equation method. Geophys.Prosp. , 1990, 38: 287-310. DOI:10.1111/gpr.1990.38.issue-3
[9] Pratt R G. Frequency-domain elastic wave modeling by finite differences: A tool for crosshole seismic imaging. Geophysics , 1990, 55(5): 626-632. DOI:10.1190/1.1442874
[10] Jo C H, Shin C, Suh J H. An optimal 9-point finite-difference frequency space 2-D scalar wave extrapolator. Geophysics , 1996, 61(2): 529-537. DOI:10.1190/1.1443979
[11] Shin C, Sohn H. A frequency-space 2-D scalar wave extrapolator using extended 25-point finite-difference operators. Geophysics , 1998, 63(1): 289-296. DOI:10.1190/1.1444323
[12] Min D J, Shin C, Kwon B D, et al. Improved frequency domain elastic wave modeling using weighted-averaging difference operators. Geophysics , 2000, 65(3): 884-895. DOI:10.1190/1.1444785
[13] 吴国忱, 梁锴. VTI介质频率-空间域准P波正演模拟. 石油地球物理勘探 , 2005, 40(5): 535–545. Wu G C, Liang K. Quasi P-wave forward modeling in frequency-space domain in VTI media. Oil Geophysical Prospecting (in Chinese) , 2005, 40(5): 535-545.
[14] 吴国忱, 梁锴. VTI介质准P波频率空间域组合边界条件研究. 石油物探 , 2005, 44(4): 301–307. Wu G C, Liang K. Combined boundary conditions of quasi-P wave within frequency-space domain in VTI media. Geophysical Prospecting For Petrole (in Chinese) , 2005, 44(4): 301-307.
[15] 吴国忱. 各向异性介质地震波传播和成像. 山东东营: 中国石油大学出版社, 2006 . Wu G C. Seismic Wave Propagation and Imaging in Anisotropy Media (in Chinese). Dongying: China University of Petroleum Press, 2006 .
[16] 吴国忱, 梁锴. VTI介质qP波方程高精度有限差分算子. 地球物理学进展 , 2007, 22(3): 896–904. Wu G C, Liang K. High precision finite difference operators for qP wave equation in VTI media. Progress in Geophysics (in Chinese) , 2007, 22(3): 896-904.
[17] 殷文, 印兴耀, 吴国忱, 等. 高精度频率域弹性波方程有限差分方法及波场模拟. 地球物理学报 , 2006, 49(2): 561–568. Yin W, Yin X Y, Wu G C, et al. The method of finite difference of high precision elastic wave equations in the frequency domain and wave-field simulation. Chinese J.Geophys. (in Chinese) , 2006, 49(2): 561-568.