TI介质一种稳定的标量波方程及其逆时偏移
杨国权1, 李河昭1, 徐文才1, 王红振1, 刘泽民2
1. 中国石油大学(华东) 地震波传播与成像实验室,山东 青岛 266580
2. 中国石油新疆油田公司采气一厂,新疆 克拉玛依 834000

作者简介: 杨国权(1962-),男,副教授,主要从事地震资料处理与解释的研究工作。

摘要

地下地层存在普遍的各向异性,随着集群计算性能和野外地震采集的迅猛发展,使得逆时偏移技术在各向异性介质中的应用成为可能。基于声学近似的标量波逆时偏移技术既能克服弹性波逆时偏移建模复杂、计算量巨大的限制,又能对地下构造准确成像,是工业界的热点问题之一。笔者从TI介质的弹性波理论出发,基于声学近似的思想,提出了一套稳健的TI介质标量波方程及逆时偏移策略。通过对模型的标量波数值模拟表明,文中提出的标量波方程能较好地压制横波干扰;Hess VTI模型的逆时偏移结果验证了将拟声波方程应用到逆时偏移的可行性。

关键词: 各向异性; 声学近似; 标量波场; 逆时偏移
中图分类号:P631.4 文献标志码:A 文章编号:1000-8918(2016)06-1237-07 doi: 10.11720/wtyht.2016.6.29
A stable scalar wave equation for reverse-time migration in TI media
YANG Guo-Quan1, LI He-Zhao1, XU Wen-Cai1, WANG Hong-Zhen1, LIU Ze-Min2
1.SWPI,China University of Petroleum(East China),Qingdao 266580,China
2. No.1 Gas Recovery Factory,Xinjiang Oilfield Company,Karamay 834000,China
Abstract

Anisotropy is widespread in subsurface.The rapid development of PC-Cluster and field seismic acquisition makes it possible to take anisotropic medium into accounted.Scalar wave equation and its RTM,based on acoustic approximation,has becoame a hot issue of the industry.It can overcome the complexity of elastic modeling and computational limitations.In this paper,the authors derived the scalar acoustic equation on the basis of elastic wave theory in anisotropic medium.According to this equation,the authors put forward a reverse time migration strategy.The numerical simulation shows that the scalar wave equations in this paper can better suppress shear wave interference.RTM result of Hess model demonstrates the feasibility when it is used to RTM.

Keyword: anisotropic medium; acoustics approximation; scalar wave field; reverse time migration

随着地震勘探程度的深入, 以及MPI、GPU技术对计算效率的极大提升[1], 叠前逆时偏移成像技术在生产中得到了广泛应用, 显著提升了横向速度剧变区域、高陡倾角构造的成像效果。地下介质存在广泛的各向异性[2], 故将逆时偏移应用到各向异性介质, 实现对各向异性介质的逆时偏移成像是勘探领域的一大热点[3, 4, 5]

声波方程和弹性波方程是实现地震波场正反传的两种传播算子, 严格上讲, 弹性波方程控制下的地震波场更加接近于野外接收波场。一些学者研究了直接运用弹性波控制方程对弹性波场做逆时外推进行RTM成像。1986年, Robert Sun和GA McMechan[6]首先提出了弹性波逆时偏移; 1994年, WF Chang和GA McMechan[7]运用激发时间成像条件将其扩展到三维; 为避免多波串扰的产生, Yan和Sava[8, 9]提出了把多分量数据作输入, 运用矢量波方程进行波场构建的方法。弹性波RTM方法存在一定的问题, 首先对耦合的方程组进行外推需要巨大的计算量, 其次是多参数的建模难度较大。此外, 野外三分量数据难以获取也限制了此方法在野外勘探中的应用。

与各向同性介质类似, 我们需要一个能够代表各向异性介质P波运动学特征的控制方程进行P波波场外推与成像, 这个方程称标量波方程(也称拟声波方程)。Alkhalifah[10]首先提出了声学近似的思想, 他利用Thomsen[11]参数ε δ 、TI介质同相轴方向的速度VP0VS0表征介质的各向异性, 在耦合的频散方程中令VS0=0简化频散关系, 得到TI介质标量波方程; Alkhalifah[12]、H Zhou[13]、Du[14]分别利用不同辅助变量对该方程作了降阶简化; Grechka[15]对各向异性中的横波干扰作了全面叙述, 认为简单令VS0=0并不能使SV波相速度处处为零; Y Zhang[16]、H Guan[17]分别对H Zhou的方程作了校正处理; Reynam、Ge zhan等通过求解耦合频散关系得到纯P波和纯SV波波动方程, 但由于隐格式求解的效率问题, 应用到逆时偏移中不合适[18, 19]

基于各向异性介质的弹性波基本理论, 笔者提出了一套稳定的TI介质标量波方程, 引入新的横波项消除了横波残留的不稳定现象。通过对HESS模型的逆时偏移表明, 该标量波方程能准确描述TI介质中P波运动学特征, 逆时偏移成像结果验证了该方程的可行性与稳定性。

1 基本原理
1.1 VTI介质弹性波速度— 应力方程

根据弹性波动力学的本构方程、运动微分方程和几何方程, 非均匀各向异性介质的弹性波波动方程可表示为:

ρ2t2U=L(CLTU)+ρF, (1)

对于VTI介质:

C=c11c11-2c66c13000c11-2c66c11c13000c13c13c33000000c44000000c44000000c662

其中:ρ 为密度, t为时间变量, U=(ux, uy, uz)T为位移矢量, C为弹性矩阵, L为偏微分算子矩阵, F=(fx, fy, fz)T为体力向量。对于完全弹性的各向异性介质, 式(1)描述了介质各质点的位移情况和传播规律。

将式(1)表示为速度— 应力形式, VTI介质的弹性波方程(假设体力项为零)改写为:

ρv1t=σ11x1+σ12x2+σ13x3, ρv2t=σ21x1+σ22x2+σ23x3, ρv3t=σ31x1+σ32x2+σ33x3, σ11t=C11v1x1+(C11-2C66)v2x2+C13v3x3, σ22t=(C11-2C66)v1x1+C11v2x2+C13v3x3, σ33t=C13v1x1+C13v2x2+C33v3x3, σ23t=C44v2x3+v3x2, σ13t=C44v1x3+v3x1, σ12t=C66v1x2+v1x13

其中:v1v2v3分别对应x1x2x3方向的速度分量; σ 11σ 22σ 33分别对应x1x2x3方向的正应力; σ 12σ 23σ 13为剪切应力。

1.2 VTI介质标量波速度— 应力方程

为方便描述, Thomsen提出了一套用于表征各向异性的参数, 表示形式为:

VP0=c33ρ, VS0=c55ρ, ε=c11-c332c33, γ=c66-c442c44, δ=(c13+c44)2-(c33-c44)22c33(c33-c44) (4)

对VTI介质的弹性波速度— 应力方程作声学近似, 即令VS0=0, 可得到弹性矩阵元素的表示形式:

c11=ρVP02(1+2ε), c13=ρVP021+2δ, c33=ρVP02, c44=0, c66=05

代入弹性波方程, 得到声学近似下的VTI介质标量波的速度— 应力方程为:

ρv1t=σ11x1, ρv2t=σ22x2, ρv3t=σ33x3, σ11t=ρvP02(1+2ε)v1x1+v2x2+1+2δv3x3, σ22t=ρvP02(1+2ε)v1x1+v2x2+1+2δv3x3, σ33t=ρvP021+2δv1x1+v2x2+v3x36

较传统的标量波方程, 该方程组为一阶形式, 易于差分格式求解; 此外以应力表示的方程组更具有物理意义。

1.3 TTI介质标量波速度— 应力方程

从VTI介质到TTI介质, 虽没有引入新的波现象, 却会使控制方程变得更为复杂, 最为显著的就是, TTI介质会引入更多的交叉导数项, 这会大大的增加计算量, 增强数值模拟的不稳定性[20]。三维观测系统下, 运用坐标变换的方法可以将VTI介质的速度— 应力方程推广到TTI介质中, 转换矩阵A表示为:

A=cosθ0sinθ-sinθsinφcosφcosθsinφ-sinθcosφ-sinφcosφcosθ, (7)

其中:θ φ 分别为TTI介质对称轴在观测坐标系下的极化角与方位角。经坐标变换后的TTI介质的速度— 应力方程表示为:

ρv1/t=d1(σ11), ρv2/t=d2(σ22), ρv3/t=d3(σ33), σ11/t=ρVP02{(1+2ε)[d1(v1)+d2(v2)]+1+2δd3(v3)}, σ22/t=ρVP02{(1+2ε)[d1(v1)+d2(v2)]+1+2δd3(v3)}, σ33/t=ρVP02[1+2δd1(v1)+d2(v2)+d3(v3)]8d1(·)=x1=cos(θ)x1'-sin(θ)sin(φ)x2'-sin(θ)cos(φ)x3', d2(·)=x2=cos(φ)x2'-sin(φ)x3', d3(·)=x3=sin(θ)x1'+cos(θ)sin(φ)x2'+cos(φ)cos(θ)x3'9

其中:d1(· )、d2(· )、d3(· )分别为x1x2x3轴方向的偏微分算子, (x1, x2, x3)表示VTI介质的坐标, ( x1', x2', x3')表示TTI介质的坐标。

1.4 TI介质标量波波场校正

利用VS0=0进行声学近似时, 各向异性参数需要满足ε δ 的条件, 否则残留的横波将会造成数值求解的不稳定。利用椭圆各向异性介质中胀缩震源不产生转化SV波的特性, 本文引入横波校正项以减小横波分量的影响。在应力项中引入横波项, 改写式(6)中的应力项为:

σ11t=ρVP02(1+2ε)v1x1+v2x2+α1+2δv3x3+ρVS02v1x1+v2x2-αv3x3, σ22t=ρVP02(1+2ε)v1x1+v2x2+α1+2δv3x3+ρvS02v1x1+v2x2-αv3x3, σ33t=ρVP021α1+2δv1x1+v2x2+v3x3-ρVS021αv1x1+v2x2-v3x310

其中α = VP0VS0(ε -δ ), 参数α 反映了介质中SV波的动力学特征。引入的横波校正项可以消除SV波波前三角结形态, 压制横波残留, 保证数值模拟波传播的稳定性。经过同1.3中的坐标变换可将此校正方程组推广到TTI介质中。

1.5 TI介质标量波逆时偏移

逆时偏移通过求解双程波动方程进行波场延拓, 利用标量波方程从震源点做正向延拓得到顺时波场S(x, y, z, t); 然后在接收点将炮记录作为边值条件做反向延拓得到逆时波场R(x, y, z, t), 对两个波场施加互相关成像条件[21, 22, 23]得到成像结果:

其中, I(x, y, z)为成像结果, S(x, y, z, t)、R(x, y, z, t)分别代表顺时波场和逆时波场。

2 模型试算

首先运用式(6)、式(8)对均匀TI介质模型进行数值模拟, 模型纵、横向采样点数为501× 501, 采样间隔为10 m, 模型参数如表1所示, 震源选择主频为30 Hz的雷克子波, 时间采样间隔1 ms。

表1 均匀TI介质参数

图1为TI介质标量波数值模拟波场拍照, 通过本文方程得到的波场拍照(图1a、1b)与传统标量波方程得到的波场拍照(图1c、1d)的对比可知:两种方程在波场模拟中均能保持稳定; 两种方程均能较准确描述P波的运动学特征; 本文TI介质标量波速度— 应力方程在横波干扰的压制上具有明显优势。

图1 TI介质标量波数值模拟波场拍照
a— VTI标量波速度— 应力方程; b— TTI标量波速度— 应力方程; c— VTI传统标量波方程; d— TTI传统标量波方程

接下来通过对均匀TI介质(ε < δ )的数值模拟来测试式(9)的横波校正效果, 模型纵、横向采样点数为501× 501, 采样间隔为10 m, 模型参数如表2所示, 震源选择主频为30 Hz的雷克子波, 时间采样间隔1 ms。由图2 TI介质(ε < δ )标量波数值模拟波场拍照可知:当各向异性参数不满足ε δ 时, 数值模拟会产生不稳定的现象(图2a、2b); 而加入横波校正后的标量波方程(图2c、2d)能够保持数值模拟的稳定性, 且能准确描述P波的运动学特征。

图2 TI介质(ε < δ )标量波数值模拟波场拍照
a— VTI传统标量波方程波场拍照; b— TTI传统标量波方程波场拍照; c— VTI波场校正方程波场拍照; d— TTI波场校正方程波场拍照

表2 均匀TI介质(ε < δ )参数

下面选用层状模型对波场校正方程进行验证, 层状模型及其角度参数如图3所示, 模型大小为6 000 m× 6 000 m, 采样间隔为10 m× 10 m, 各层的模型参数如表3所示。在地表3 000 m、深度1 000 m处施加一个主频为30 Hz的雷克子波, 得到的波场记录如图4所示。图4a为利用各向异性完全弹性波方程所得的波场拍照, 其中箭头指示了波场中的横波项; 图4b为利用TTI波场校正方程所得的波场拍照。由图4a、4b箭头处对比可知, TTI波场校正方程能够适用于非均匀介质, 且能稳定地压制横波干扰。

表3 层状TTI模型参数

图3 层状模型(左)及其角度参数(右)

图4 层状模型数值模拟波场拍照
a— 各向异性完全弹性波方程波场拍照; b— TTI波场校正方程波场拍照

最后通过Hess VTI模型验证本文方程应用到逆时偏移中的可行性, Hess VTI模型的P波速度场、各向异性参数ε δ 如图5所示。其中P波速度范围为1 500~4 500 m/s, ε δ 的范围分别为0~0.28 和0~0.16, 300炮激发, 炮间距30 m, 每炮排列数为600, 检波器间隔为10 m, 震源采用主频 25 Hz的雷克子波, 采样间隔4 ms。逆时偏移中, 采用互相关成像条件以及laplacian滤波处理。

图5 Hess VTI模型
a— P波速度场; b— 各向异性参数ε ; c— 各向异性参数δ

分别利用传统标量波方程与本文方程作为控制方程对Hess VTI模型进行逆时偏移成像, 所得成像剖面如图6所示。从成像结果的对比看出:传统标量波方程产生的横波残余影响了RTM的成像结果; 而本文新推导的标量波速度— 应力方程横波残留小、同相轴清晰、地下盐丘轮廓明显, 且对右侧高陡构造能够更准确成像。

图6 Hess VTI模型逆时偏移成像结果
a— 传统标量波方程成像结果; b— 标量波速度— 应力方程成像结果

3 结论

基于弹性波的基本理论, 本文推导出了声学近似下的标量波速度— 应力方程, 数值模拟结果表明, 该方程能够准确描述P波的传播规律。该方程具有更直接的物理意义且易于求解, 与传统标量波方程对比表明该方程在数值模拟中具有更高的稳定性; 引入的横波校正项在压制横波干扰上具有较好的效果; 对Hess VTI模型的逆时偏移试算取得了理想的结果, 显示了本文方程的优越性和实用性。

The authors have declared that no competing interests exist.

参考文献
[1] Kim Y, Cho Y, Jang U, et al. Acceleration of stable TTI P-wave reverse-time migration with GPUs[J]. Computers & Geosciences, 2013, 52: 204-217. [本文引用:1]
[2] 吴国忱. 各向异性介质地震波传播与成像[M]. 北京: 石油大学出版社, 2006. [本文引用:1]
[3] Baysal E, Kosloff D D, Sherwood J W C. Reverse time migration[J]. Geophysics, 1983, 48(11): 1514-1524. [本文引用:1]
[4] Dellinger J A. Anisotropic seismic wave propagation[D]. Palo Alto: Stanford University, 1991. [本文引用:1]
[5] Farmer P A, Jones I F, Zhou H, et al. Application of reverse time migration to complex imaging problems[J]. First Break, 2006, 24(9): 1-6. [本文引用:1]
[6] Sun R, McMechan G A. Line sources for seismic modeling by finite differences in inhomogeneous media[J]. Geoexploration, 1987, 24(3): 183-196. [本文引用:1]
[7] Chang W F, McMechan G A. 3-D elastic prestack, reverse-time depth migration[J]. Geophysics, 1994, 59(4): 597-609. [本文引用:1]
[8] Yan J, Sava P. Elastic wavefield imaging with scalar and vector potentials[C]//77th Annual International Meeting, SEG, Expand ed Abstracts, 2007. [本文引用:1]
[9] Yan J, Sava P. Elastic wave-mode separation for VTI media[J]. Geophysics, 2009, 74: WB19-WB32. [本文引用:1]
[10] Alkhalifah T. Acoustic approximations for processing in transversely isotropic media[J]. Geophysics, 1998, 63(2): 623-631. [本文引用:1]
[11] Thomsen L. Weak elastic anisotropy[J]. Geophysics, 1986, 51(10): 1954-1966. [本文引用:1]
[12] Alkhalifah T. An acoustic wave equation for anisotropic media[J]. Geophysics, 2000, 65(4): 1239-1250. [本文引用:1]
[13] Zhou H, Zhang G, Bloor R. An anisotropic acoustic wave equation for VTI media[C]//68th EAGE Conference and Exhibition incorporating SPE EUROPEC, 2006. [本文引用:1]
[14] Du X, Fletcher R P, Fowler P J. A new pseudo-acoustic wave equation for VTI media[C]//70th EAGE Conference and Exhibition incorporating SPE EUROPEC, 2008. [本文引用:1]
[15] Grechka V, Zhang L, Rector III J W. Shear waves in acoustic anisotropic media[J]. Geophysics, 2004, 69(2): 576-582. [本文引用:1]
[16] Zhang Y, Zhang G. One-step extrapolation method for reverse time migration[J]. Geophysics, 2009, 74(4): A29-A33. [本文引用:1]
[17] Guan H, Dussaud E, Denel B, et al. Techniques for an efficient implementation of RTM in TTI media[C]//81st Annual International Meeting, SEG, Expand ed Abstracts, 2011. [本文引用:1]
[18] Pestana R C, Ursin B, Stoffa P L. Separate P-and SV-wave equations for VTI media[C]//12nd International Congress of the Brazilian Geophysical Society, 2011. [本文引用:1]
[19] Zhan G, Pestana R C, Stoffa P L. An acoustic wave equation for pure P wave in 2D TTI media[C]//12nd International Congress of the Brazilian Geophysical Society, 2011. [本文引用:1]
[20] 张岩, 吴国忱. TTI 介质叠前逆时偏移成像研究综述[J]. 地球物理学进展, 2013, 28(1): 409-420. [本文引用:1]
[21] Claerbout J F. Toward a unified theory of reflector mapping[J]. Geophysics, 1971, 36(3): 467-481. [本文引用:1]
[22] Kosloff D D, Baysal E. Migration with the full acoustic wave equation[J]. Geophysics, 1983, 48(6): 677-687. [本文引用:1]
[23] Schleicher J, Costa J C, Novais A. A comparison of imaging conditions for wave-equation shot-profile migration[J]. Geophysics, 2008, 73(6): S219-S227. [本文引用:1]