地球物理学报  2011, Vol. 54 Issue (2): 329-335   PDF    
引用本文
赵海波, 陈百军, 李奎周, 等. 黏弹性介质VSP记录模拟及在估算Q值研究中应用. 地球物理学报, 2011, 54(2): 329-335, doi: 10.3969/j.issn.0001-5733.2011.02.008.  
ZHAO H B, CHEN B J, LI K Z, et al. VSP record numerical modeling in viscoelastic media and its application in the study of Q-value estimation method. Chinese J. Geophys. (in Chinese), 2011, 54(2): 329-335, doi: 10.3969/j.issn.0001-5733.2011.02.008.  
黏弹性介质VSP记录模拟及在估算Q值研究中应用
赵海波, 陈百军, 李奎周, 成德安     
大庆油田有限责任公司勘探开发研究院,大庆 163712
收稿日期 2010-04-06, 2010-05-28 收修定稿
基金项目: 国家重点基础研究发展计划(973)项目(2009CB219300)资助
作者简介: 赵海波,男,1979年生,博士,主要从事地震波传播理论及正演模拟、地震岩石物理分析、叠前地震储层预测等研究. E-mail:zhaohaibo@petrochina.com.cn
摘要: 采用标准线性固体模型,本文建立了黏弹性介质完全匹配层吸收边界的高阶速度-应力交错网格有限差分算法,并对黏弹性介质中的地震波传播进行了数值模拟.基于黏弹性波动方程正演模拟提供的零偏VSP全波场数据,本文进行了质心频移法计算Q值的反演分析.结果表明,反射波、转换波及短程多次波对频谱的影响较大,对Q值反演造成一定误差.本文的结论为实际零偏VSP资料估算地层Q值提供了有益的指导.
关键词: 黏弹性介质      标准线性固体模型      数值模拟      零偏VSP      Q值     
VSP record numerical modeling in viscoelastic media and its application in the study of Q-value estimation method
ZHAO Hai-Bo, CHEN Bai-Jun, LI Kui-Zhou, CHENG De-An     
Exploration and Development Research Institute of Daqing Oilfield Company Ltd., Daqing 163712, China
Abstract: Using the standard linear solid model, a high-order velocity-stress staggered-grid finite-difference scheme with the perfectly matched layer method was proposed for simulating seismic wave propagation in viscoelastic media. With the full-wave field data of viscoelastic modeling of zero-offset VSP from several numerical experiments, the frequency shift method to calculate Q-value was analyzed. The numerical results showed that reflected waves, transformed waves and short-path multiples have relatively large influences on the frequency spectra of received data, and make further impacts on the Q-value inversion. The conclusions taken from the paper can provide helpful guidance for the estimation of subsurface interval Q-value with oilfield zero-offset VSP data.
Key words: Viscoelastic medium      Standard linear solid model      Numerical modeling      Zero-offset VSP      Q-value     
1 引 言

随着地震勘探需求的逐渐增大和技术的快速发展，人们越来越关心黏弹性介质中地震波的传播特征.通过理解黏弹性介质，一方面可消除或补偿地层吸收衰减作用，实现地震资料保幅处理和提高分辨率，另一方面可利用地震波所携带的黏弹性信息进行地层岩性预测和烃类检测.

通过正演模拟技术，可以对黏弹性介质中地震波场有直观、形象的认识和理解.此外，该技术可为地震资料处理解释中各种方法的应用提供标准的先验数据体，用以检验方法的适用性.目前，黏弹性介质数值模拟具有代表性的方法主要有有限差分法[1~3]、有限元法[4]和伪谱法[5, 6].综合考虑处理复杂介质模型的灵活性和计算效率及精度，交错网格有限差分法在勘探地球物理领域的应用最为普遍.

地震勘探中，通常采用地层品质因子(Q值)衡量地下介质对地震波能量衰减的强弱.已有的研究结果表明，Q值与岩石物性、孔隙流体类型及流体饱和度等因素有十分密切的关系[7, 8].利用VSP 数据可以估算地层Q值，已研究出的主要方法有脉冲宽度法[9]、脉冲振幅衰减法[10]、谱比法[11]、质心频移法[12]等.基于小波变换技术一些学者提出了相应的方法.高静怀等[13]在小波域发展了一种计算瞬时频率的方法，并在此基础上提出了计算VSP 资料Q值的方法.以Morlet小波为基本小波，赵伟等[14]提出了利用VSP资料在时-频域计算Q值的方法.在上述方法中，质心频移法在现场应用最为广泛.该方法主要通过分析信号质心频率的变化来计算Q值.以往讨论质心频移法，利用射线追踪建立正演数据来进行讨论[12].射线理论的优点是计算速度快，能很好地反映地震波的运动学特点.但射线理论在反映地震波的动力学特征方面效果不好，不能体现出转换波和多次波等复杂波场对计算Q值的影响.

基于大庆油田徐家围子地区徐深21-1 的VSP数据建立的速度模型，本文采用完全匹配层吸收边界的高阶交错网格有限差分算法进行非均匀黏弹性介质数值模拟.然后，在黏弹性介质VSP 正演模拟基础上，利用多个测试模型对质心频移法计算Q值进行了相应的分析.

2 标准线性黏弹性介质波动方程

在唯象黏弹性理论中，黏弹性模型的基本假设是应力张量与应变张量的时间历史有关，即存在如下关系[15]:

(1)

式中*代表时间褶积；Tij和ekl分别为应力张量和应变张量；Cijkl为与时间有关的四阶张量，称为松弛方程.对于松弛方程Cijkl的具体数学描述，不同的模型有不同的计算表达式.常用的模型有Maxwell黏弹性模型、Kelvin黏弹性模型、标准线性固体(SLS)黏弹性模型或其他三元黏弹性模型(如Poyting-Thomson体等).这些模型中，SLS黏弹性模型常被用于有限差分数值模拟中，主要是因为利用差分算子较容易实现应力和应变关系方程[5].此外，Toverud和Ursin[16]的研究结果表明，SLS 黏弹性模型符合实际地震波传播规律.

对于SLS模型，松弛方程可由式(2)给出

(2)

式中MR 为松弛模量，H(t)为Heaviside函数.该松弛方程等效于N个标准线性单元并联，τσn和τεn分别为第n个耗散机制下应力和应变的松弛时间.

在各向同性、非均匀黏弹性介质情况下，本构关系方程可简化为

(3)

考虑到式(2)，式(3)中的系数表达形式为

(4)

(5)

(6)

式中λ 和μ 为弹性介质的拉梅系数，τεnp和τεns分别是介质第n个机制的纵波和横波黏弹性应变松弛时间，τσn是介质第n个机制的黏弹性应力松弛时间.

黏弹性介质的动量守恒方程与弹性介质的相同，即

(7)

式中Vi为黏弹性介质的质点速度，ρ 为介质密度，fi为体力.

3 标准线性黏弹性介质交错网格有限

差分方程现考虑x-z二维平面问题，且式(5)和(6)中的N=1(即只考虑一种耗散机制)，利用交错网格有限差分算法将方程(3)和(7)离散化，可得到标准线性黏弹性介质控制系统的差分格式方程:

(8a)

(8b)

(9a)

(9b)

(9c)

(10a)

(10b)

(10c)

式中i和j是空间离散指数；n为时间离散指数；Δt为时间离散步长；Rxx、Rzz和Rxz为记忆变量，引入记忆变量的目的是消除式(2)中的褶积以绕开耗时的褶积运算[17].需要说明的是，在推导方程(10)时，考虑了下面的关系:

(11a)

(11b)

(11c)

这种平均作法并不影响数值模拟的精度，因为本文差分格式的时间精度为二阶[18].同时，在数值模拟中，方程(9)等式右边最后一项有关记忆变量的运算需要考虑方程(11).方程(9)和(10)中，系数c1、c2、c3、c4、d1 和d2 的表达式分别为

式中τεp 和τεs 分别是纵波和横波的应变松弛时间，τσ 是应力松弛时间.松弛时间与品质因子存在如下关系[19]:

(12a)

(12b)

(12c)

式中Qp 和Qs 分别是纵波和横波的品质因子，f0 为声源的中心频率.

方程(8)~(10)中，Dx和Dz分别表示x和z方向的离散化微分算子，

式中am为有限差分加权系数；Δx和Δz分别为x方向空间步长和z方向空间步长；M为导数算子长度，对于空间八阶精度情况，M取为4.am在各阶精度时的数值可参考文献[20].

此外，黏弹性介质的稳定性条件与弹性介质的类似[18]，这里不再赘述.吸收边界条件是数值模拟的重要研究课题之一，其效果的好坏对数值模拟精度的影响相当大.本文采用完全匹配层吸收边界条件进行黏弹性介质的数值模拟.该技术的具体描述和差分格式实现可参考文献[21].

4 数值模拟

为考察本文提出算法的有效性及稳定性，考虑图 1所示的非均匀介质模型.图 1 为基于徐家围子地区徐深21-1的VSP数据建立的层状速度模型.图中小方块代表接收器，水平方向代表地面接收器，垂直方向代表VSP接收器.地层品质因子在图中已标出(Qp=Qs)，地层纵波速度vp 从上到下分别为2326、2462、2673、2735、2977、4040、4983、4836、5200m/s.相应的横波速度为vs =vp/m/s，相应的介质密度为ρ =230× (3.28vp)0.25kg/m3[22].图 2 是不同时刻的波场快照图，图 2a是黏弹性介质情况，图 2b为与之对应的完全弹性介质情况.数值模拟中，震源是主频为42 Hz的Ricker子波，空间网格步长Δx和Δz均为4 m，时间步长Δt为0.35 ms.比较而言，黏弹性介质中波场扩散存在能量衰减，同时波形变宽，表现出很强的频散.这些特点在图 3中显现的更加直观，图中给出了两个接收器在黏弹性介质和完全弹性介质情况下的波形和振幅谱比较，两个接收器在垂直方向深度分别为1200 m 和3200 m(VSP接收器线).从图中看到，黏弹性介质中的高频成分衰减快，频谱的中心频率向低频方向移动.这说明，介质的黏弹性对频谱形态和中心频率有很明显的影响.另外，从波场快照上看，人工边界产生的虚假波能量很弱，吸收边界条件吸收外行波效果理想.

图 1 基于徐深21-1零偏VSP资料建立的层状模型 Fig. 1 Layered model established with zero-offsetVSP data of Xushen 21-1 well
图 2 波场快照图 (a)黏弹性介质；(b)完全弹性介质. Fig. 2 Wave field snapshots (a) Viscoelastic model； (b) Fully elastic model.
图 3 两个垂直深度点检波器波形和频谱显示 (a) 1200 m； (b) 3200 m. Fig. 3 Waveforms and spectra of two receivers in VSP receiver array

图 4为共炮点的零偏VSP 与地面地震联合显示结果.图中对波场进行了分析，标出了各种声波模式以及图 1中的各界面反射位置.零偏VSP 与地面地震有非常好的对应关系，表明零偏VSP 可用于地面地震的层位标定以及识别多次波.

图 4 共炮点的VSP与地面地震联合剖面显示 (a)零偏VSP；(b)地面地震. Fig. 4 Joint section of VSPand surface seismic from the same source (a)Zero-offset VSP；(b)Surface seismic
5 质心频移法计算Q值分析

以往讨论质心频移方法时，都是利用射线追踪和褶积手段建立正演数据来进行相应的讨论.这种手段容易实现，但是波场相对简单，不能真正的说明质心频移法计算Q值的特点.下面将基于黏弹性波动方程正演模拟提供的零偏VSP全波场数据，详细讨论质心频移法.

本节共讨论四种模型情况，即各层速度不同、密度相同模型，速度相同、密度不同模型，速度和密度均相同模型，速度和密度均不同模型.模型共有8层，具体参数见表 1~表 4，其中，Qp = Qs.数值模拟中，震源是主频为225Hz的Ricker子波，空间网格步长Δx和Δz均为5m，时间步长Δt为0.1ms.第一个接收点距离井口的距离是10m，共119个接收点，相邻两个接收点间距为5m.

表 1 速度不同、密度相同模型 Table 1 Model with different velocity and identical density
表 2 速度相同、密度不同模型 Table 2 Model with identical velocity and different density
表 3 速度和密度均相同模型 Table 3 Model with identical velocity and identical density
表 4 速度和密度均不同模型 Table 4 Model with different velocity and different density

图 5给出模型Q值和反演Q值的比较.反演过程中，需要沿着初至开一定宽度的时窗，在这个时窗内进行相应处理.另外，由于有些检波器放在了地层界面处，使得计算的Q发生较大波动，因此需要平滑处理.从这四个数值实验可看到，速度和密度均相同情况下的反演结果最好，这主要是因为波场中反射波和转换波的干扰很小(虽然给定模型速度相同，但是各层的Q不同会引起各层的速度发生微小变化，同样会有反射波和转换波的出现).速度和密度均不同情况下的反演结果最差，主要是波场中反射波和转换波的能量强，使得初至波受到的干扰大.同时也注意到，随深度增加误差逐渐增大，地层界面处Q值会发生一定程度的波动.

图 5 反演Q值与模型Q值比较 (a)速度不同、密度相同模型；(b)速度相同、密度不同模型；(c)速度和密度均相同模型；(d)速度和密度均不同模型. Fig. 5 Comparison of inverted Q-value and model Q-value (a) Different velocity and identical density； (b) Identical velocity and different density；(c) Identical velocity and identical density； (d) Different velocity and different density.

对结果的分析表明，质心频移法可较好地计算零偏VSP 资料的Q值.虽然估算的Q值与真实Q值之间有一定的误差，但反演各层Q值的相对关系可以保证.这些结果对于实际零偏VSP资料估算地层Q值，提供了很好地理论指导和解释上的借鉴.当利用反演得到的连续Q值进行反Q滤波时，要按层位划分平均得到各个层段的Q值，同时Q值可以适当调整，但是相对关系不能变动.

6 结 论

采用SLS模型，本文提出了黏弹性介质高阶交错网格有限差分算法，并使用完全匹配层技术处理人工边界.通过非均匀介质模型测试，本文算法稳定.数值模拟结果表明，介质的黏弹性可导致波速改变和能量衰减.基于黏弹性波动方程正演模拟提供的零偏VSP全波场数据，本文对多个模型进行了质心频移法计算Q值的反演和分析.结果表明，反射波、转换波及短程微屈多次波对Q值反演产生较大影响，但反演得到的Q值相对关系较为可靠.

参考文献
[1] Robertsson J O A, Blance J O, Symes W W. Viscoelastic finite-difference modeling. Geophysics , 1994, 59: 1444-1456. DOI:10.1190/1.1443701
[2] Wang X M, Kevin D. Viscoelastic wave modeling using a staggered high-order finite-difference in inhomogeneous media. Expanded Abstracts, 6th SEGJ International Symposium: Imaging and Technology, Tokyo, 2003
[3] Saenger E, Bohlen T. Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid. Geophysics , 2004, 69: 583-591. DOI:10.1190/1.1707078
[4] Eliane B, Abdelaaziz E, Patrick J. A mixed finite element approach for viscoelastic wave propagation. Computational Geosciences , 2004, 8: 255-299.
[5] Carcione J M. Seismic modeling in viscoelastic media. Geophysics , 1993, 58: 110-120. DOI:10.1190/1.1443340
[6] 单启铜, 乐友善. PML边界条件下二维粘弹性介质波场模拟. 石油物探 , 2007, 46(2): 126–130. Shan Q T, Le Y S. Wavefield simulation of 2-D viscoelastic medium in perfectly matched layer boundary. Geophysical Prospecting for Petroleum (in Chinese) , 2007, 46(2): 126-130.
[7] Winkler K W, Nur A. Attenuation: effects of pores and frictional sliding. Geophysics , 1982, 47(1): 1-15. DOI:10.1190/1.1441276
[8] William F, Murphy Ⅲ. Effects of partial water saturation on attenuation in Massilon sandstone and Vycor porous glass. J. Acoust. Soc. Am. , 1982, 71(6): 1458-1467.
[9] Wright C, Hoy D. A note on pulse broadening and anelastic attenuation in near-surface rocks. Phys. Earth Planetary Int., , 1981, 25(1): 1-8.
[10] Brzostowski M, McMechan G. 3-D tomographic imaging of near-surface seismic velocity and attenuation. Geophysics , 1992, 57(3): 396-403. DOI:10.1190/1.1443254
[11] Hauge P S. Measurements of attenuation from vertical seismic profiles. Geophysics , 1981, 46(11): 1548-1558. DOI:10.1190/1.1441161
[12] Quan Y, Harris J M. Seismic attenuation tomography using the frequency shift method. Geophysics , 1997, 62(3): 895-905. DOI:10.1190/1.1444197
[13] 高静怀, 杨森林. 利用零偏移VSP资料估计介质品质因子方法研究. 地球物理学报 , 2007, 50(4): 1198–1209. Gao J H, Yang S L. On the method of quality factors estimation from zero-offset VSP data. Chinese J. Geophys. (in Chinese) , 2007, 50(4): 1198-1209.
[14] 赵伟, 葛艳. 利用零偏移距VSP资料在小波域计算介质Q值. 地球物理学报 , 2008, 51(4): 1202–1208. Zhao W, Ge Y. Estimation of Q from VSP data with zero offset in wavelet domain. Chinese J. Geophys. (in Chinese) , 2008, 51(4): 1202-1208.
[15] Christensen R M. Theory of Viscoelasticity—An Introduction. New York: Academic Press, 1982 : 14 -16.
[16] Toverud T, Ursin B. Estimation of viscoelastic parameters from zero-offset VSP data. Expanded Abstracts, 61st EAGE Conference, Helsinki, Session: 4.29, 1999
[17] Carcione J M, Kosloff D, Kosloff R. Wave propagation simulation in a linear viscoacoustic medium. Geophys.J. Roy. Astr. Soc. , 1988, 93: 393-407.
[18] Robertsson J O A, Blance J O, Symes W W. Viscoelastic finite-difference modeling. Geophysics , 1994, 59: 1444-1456. DOI:10.1190/1.1443701
[19] Hestholm S O, Ruud B O. 2D finite-difference viscoelastic wave modeling including surface topography. Geophys. Prop. , 2000, 48: 341-373.
[20] 张海澜, 王秀明, 张碧星. 井孔的声场和波. 北京: 科学出版社, 2004 : 140 -150. Zhang H L, Wang X M, Zhang B X. Acoustic Field and Wave in Borehole (in Chinese). Beijing: Science Press, 2004 : 140 -150.
[21] 赵海波, 王秀明, 王东, 等. 完全匹配层吸收边界在孔隙介质弹性波模拟中的应用. 地球物理学报 , 2007, 50(2): 581–591.
[22] Gardner G H F, Gardner L W, Gregory A R. The diagnostic basis for stratigraphic straps. Geophysics , 1974, 39: 770-780. DOI:10.1190/1.1440465