地球物理学进展  2016, Vol. 31 Issue (2): 836-844   PDF    
基于拟声波一阶应力-速度方程的TI介质叠前逆时偏移
徐文才, 李振春, 王姣, 李河昭    
中国石油大学(华东)地球科学与技术学院, 青岛 266580
摘要:随着计算机硬件技术发展和对野外数据高分辨率勘探需求的增加,叠前逆时偏移技术逐渐的运用到实际生产中;又鉴于实际地层中各向异性的普遍存在性,故实现各向异性介质的叠前逆时偏移成像很有必要.文中从各向异性介质中弹性波的基本理论出发,导出各向异性介质中的拟声波一阶速度-应力方程,并基于此方程提出一套稳健的叠前逆时偏移策略.对逆冲模型的数值模拟表明在旋转角变化的TTI介质(具有倾斜对称轴的横向各向同性介质)中该拟声波方程稳定性要好于基于频散关系推导出的拟声波方程;VTI_HESS模型的逆时偏移结果表明较常规声波方程,新导出的拟声波方程成像更清晰,深处保幅效果更加明显.
关键词各向异性介质     声学近似     速度-应力方程     拟声波方程     逆时偏移    
TI medium reverse time migration based on quasi acoustic first-order velocity- stress equation
XU Wen-cai, LI Zhen-chun, WANG Jiao, LI He-zhao    
School of Geoscience, China University of Petroleum (East china), Qingdao 266580, China
Abstract: With the development of computer computing rate and hardware level, especially the GPU parallel computing technology is introduced, prestack reverse time migration gradually used up to the field of seismic data. In view of the universal existence of anisotropy in the actual formation, considering anisotropy in the reverse time migration imaging is very necessary. In this paper, we will derive quasi acoustic first-order velocity-stress equation based on the basic theory of elastic wave in anisotropic medium, and According to this equation, put forward a set of robust prestack reverse time migration strategy. Numerical simulation of the thrust model shows that this new quasi acoustic wave equation in wave field extrapolation is more stability than the quasi acoustic wave equation which derived from dispersion relation in the TTI medium with the varying rotation. The HESS model of reverse time migration imaging showed that compared with the conventional acoustic wave equation, quasi acoustic wave equation imaging more clearly, the effect of amplitude preserving in depths is more obvious.
Key words: anisotropic medium     acoustics approximation     velocity stress equation     quasi acoustic wave equation     reverse time migration    
0 引 言

随着地震勘探对深层地层的深入,对复杂地区的保幅成像越来越受到处理人员的重视了.鉴于地层各向异性的广泛性,如何对各向异性介质进行保幅成像也逐渐成为地震勘探的热门问题.由于地震波在各向异性介质中一般会发生横波分裂、不同方向传播的地震波速度不同、纵横波耦合强烈等现象;灰岩页岩相间的薄层交互层、有固定排列方向的垂直裂隙以及沉积岩中的旋回性砂泥岩薄互层都会引发地震波的各向异性.针对地球介质的实际情况进行地震勘探和高精度成像已经越来越重要和艰难.叠前逆时偏移技术利用全波长成像原理,对高陡构造、横向变速以及棱柱波等成像有较强的优势,具有更高的成像精确和更丰富的有效信息,故本文将其引入到各向异性介质的成像中优势是十分明显.

各向异性介质的逆时偏移成像一般有两种方式:一是直接利用弹性波的控制方程作为地震波场正向、反向传播的算子,进而进行成像.Dellinger和Etgen(1990),Yan和Sava等(2010)通过求解各向异性介质弹性波的Christoffel方程得到P波和S波的偏振方向,从而求得空间域波场分离算子.逆时偏移中利用分离出纯P波、纯S波进行偏移成像.从某种意义上说,利用这种弹性波方程更加接近实际地质体,但实现的过程中存在许多问题:一是多参数初始模型的建立,这个对野外数据采集要求太高,各向异性介质参数将更多,很难实现;二是弹性波控制方程描述的地震波场是耦合的,这对计算机的计算效率和内存要求加大;三是在成像过程中,每个空间采样点的地震波场正、反传的每个成像时间步长都必须使用空间域波场分离算子实现波场分离,这太消耗时间和计算机的存储空间.考虑到以上困难,在实际生产中,主要还是使用标量波方程进行偏移成像.

另一中方式就是各向异性下的拟声波近似方程,即采用声学近似的方式导出控制方程.Alkhalifah(1998)最先提出qP波的控制方程:在频散方程中令对称轴上的SV波速度为零,简化频散关系,得到qP波方程;还有Zhou等(2006)Du等(2008)等都对频散关系作了相应的优化简化,分别导出各自的qP波方程,但这样的简化也就决定其SV波并不能完全的消除,同时方程必须在介质满足特定的条件下才能保持稳定;Zhang等(2011)等人对Du的拟声波方程进行改进,通过保持能量守恒的方式证明其方程的稳定性;Pestana等(2011)Zhan等(2011)等人则是通过求解出解耦的P波和SV波的频散关系从而得到纯P波和纯SV波波动方程,由于方程需要用隐格式求解,而隐格式求解效率较低,用在计算量较大的逆时偏移中并不是很合适.

本文根据笔者导出的拟声波一阶应力-速度方程,提出一套稳定的各向异性介质叠前逆时偏移技术.为了避免由各向异性参数引起的地震波传播的不稳定现象,文中根据Fletcher等(2009)理论,引入横波项,消除由残留的横波的不稳定现象.由本文导出的拟声波方程为一阶偏微分方程组,能很好的适用交错网格差分形式.通过对VTI_HESS模型和逆冲模型的数值模拟和逆时偏移,表明该拟声波方程能准确的描述TI介质中地震波运动学特征,并能很好的实现各向异性介质逆时偏移.

1 原理与方法1.1 TI介质拟声波一阶速度-应力方程

各向异性介质的拟声波方程是研究该介质下地震波正演、偏移成像、反演等的重要依据,传统的拟声波方程都是由频散关系通过不同方式推导而来,其实它们之间可以相互转化,只是由于引入的辅助函数不同,才使得求解出来的波动方程有所差别.本文将基于弹性波的基本理论,即本构方程、运动微分方程还有几何方程,推导TI介质中的拟声波一阶速度-应力方程.

我们知道TTI介质可由VTI介质通过坐标变换而来,所以VTI介质可认为是TTI介质的一种特例.下面以VTI介质为例,在观测坐标系下,其微分运动方程可写为

本构方程可写为

其中σij(ij=1、2、3)表示应力,εij(ij=1、2、3)为应变,vi(ij=1、2、3)为质点的速度,Cij(ij=1、2、3、4、5、6)为弹性矩阵系数.本构关系(2)式两边对时间求偏导:

由几何方程我们可得到应变和速度分量的关系为

将(4)式带入本构方程的微分式(3)中,再联合运动微分方程可得

(5)式是VTI介质中弹性波一阶速度-应力方程,下面对其进行声学近似,刚度系数与Thomsen参数关系为

其中vpvsεδγ为5个Thomsen参数,Cij(ij=1、2、3、4、5、6)为弹性矩阵系数.令y=0可得

带入弹性波方程,得到声学近似下的的一阶方程组为

(8)式就是所要求的拟声波一阶速度-应力方程(为了更清楚表示物理量,将下标由1、2、3改为xyz).较其传统的拟声波方程,该方程组为一阶形式,更容易求解;另一方面该方程组以应力表示更具有物理意义.

1.2 从VTI到TTI

在三维观测系统内,设TI介质的对称轴和OXZ平面内的观测系统Z轴的夹角为θ,这个角度称为极化角;在OXY平面内与X轴的夹角为φ,称为方位角.这样只要对(8)进行相应的坐标变换,就能将VTI介质中的一阶方程转换到TTI介质中.设转换矩阵为C(详细的推导过程参考附件B)

变量(xyz)到(x′y′z′)的变化可以表示为

其中d1(·)、d2(·)、d3(·)表示在xyz轴方向上的偏微分算子,(xyz)表示VTI介质中的坐标,(x′y′z′)表示为转换到TTI介质中的坐标.

将(9)式带入到方程(8)式中就可得到TTI介质中拟声波一阶速度-应力方程为

这样(10)就是所要求解的一阶速度-应力方程,也是本文正演和偏移成像的基本方程.

对二维均匀各向异性介质试验,各向异性参数:vp0=3000 m/s,ε=0.24,δ=0.1,其中TTI介质的XOZ平面内的旋转角度:θ=45°.本文新推导的qP波数值模拟效果如图 1所示,为了能够对比本文推导的一阶方程与基于频散关系推导的二阶方程差异,笔者选取近年来较为常用的Zhou等(2006)所推导的二阶方程做对比(其方程是由Alkhalifah的四阶方程降阶转变而来),图 2Zhou等(2006)qP波方程波场快照.

图 1 新qP波方程波场快照,VTI(a)、TTI(b) Fig. 1 The wave field snapshot of new qP wave equation, VTI (a), TTI (b)

图 2 Zhou(2006)qP波方程波场快照,VTI(a)、TTI(b) Fig. 2 The wave field snapshot of Zhou’s qP wave equation, VTI (a), TTI (b)

由模型一的波场快照图 12可知,当各向异性参数εδ时,Zhou的方程和新推导的方程的qP波波前面呈椭圆状,符合各向异性介质中qP波的运动学特征;而且都没有出现频散等不稳定现象,只是横波干扰不一样,Zhou的方程的横波干扰明显强于新推导一阶方程,主要是由于二阶方程放大了横波干扰项.

1.3 从能量守恒的角度进行稳定性分析

各向异性介质中的拟声波方程主要由两种方式推导:一是由P-SV波的频散关系简化变换得到,另一就是像本文由弹性波的基本理论推导而来.针对以上两种方式的拟声波方程,笔者将从能量守恒的角度来分析其控制方程的稳定性.

在物理学中我们常把动能表示为,势能表示为Ep(t)=∫Fdh.其中势能即保守力做的功,在弹性波动力学中,其内力就是应力,也就是说其势能为质点在应力作用下产生位移所做的功,可以表示为Ep(t)=∫σdu(其中u是位移矢量,σ是应力矢量),再利用位移和应变关系:,将该关系式带入势能表达式中得到势能密度为

为了方程的简洁,设密度为常数1,这样我们可以将动能密度表示为

由上述(11)、(12)式我们就能定义弹性波的能流密度为(单位质量的质点所携带的能量)

其中σε分别为应力矢量[σ11σ22σ33σ23σ13σ12]和应变矢量[ε11ε22ε33ε23ε13ε12],根据应力应变关系可以将上式能流密度表示为

要证明声学近似的qP一阶方程组在传播过程中能量是守恒的,即证明:

将本文所推导的拟声波一阶速度-应力方程组带入动能和势能表达式中,求取能流密度时间偏导(详细的推导过程请参考附件A,正文中笔者将直接写出计算结果),经过简化计算可以得到动能和势能对时间导数的表达式为(VTI介质)

其中σVσH表示垂向的应力和水平方向的应力,我们设边界为无穷远,即在边界上vxvyvzσVσH的值都趋向于零值,所以(14)中两式相加就得到我们要证明的结论:.然而我们利用相同的方式,求取Du等(2008)拟声波方程所对应的动能和势能时间导数(详细的推导过程请参考附件A),其表达式为

从各个偏导数的系数可以看出:.

2 模型试算2.1 逆冲模型数值模拟和逆时偏移

为了对比新推导的qP方程和传统的qP波方程数值模拟的稳定性及各向异性参数的适应性,对较为复杂的逆冲模型进行试算,将模型改为二维TTI介质模型,如图 3分别是逆冲模型的Vp0εδθ模型,其中θ是指ZOX平面内的旋转角度.

图 3 逆冲模型Vp0(a)、ε(b)、δ(c)和θ(d)模型 Fig. 3 Vp0(a)、ε(b)、δ(c) and θ(d)of Thrust model

图 4图 5分别为新推导的qP方程和Zhou等(2006)推导的qP方程正演结果,由上图可知新推导的方程对逆冲模型的数值稳定性要好于Zhou等(2006)方程,特别是随着波传播的时间增大,差异更明显.这是由于Zhou、X.Du等人推导的TTI介质中的qP波方程为二阶方程,这使得横波残留波场也成二次方增加,特别是在TTI介质的各向异性参数变换较大的区域更为明显.

图 4 新qP方程的波场快照t=1000 ms(a)、t=2000 ms(b) Fig. 4 The wave field snapshot of new qP wave equation, t=1000 ms (a), t=2000 ms (b)

图 5 Zhou(2008)方程的波场快照 t=1000 ms(a)、t=2000 ms(b) Fig. 5 The wave field snapshot of Zhou’s equation, t=1000 ms (a), t=2000 ms (b)

图 6分别为Zhou方程与新qP方程RTM成像结果,对比可见,新方程时的成像结果都比Zhou方程成像结果更加清晰,对各向异性介质的刻画更加符合真实的模型,在高陡构造成像保幅效果更好,由此,体现了新方程的在TI介质逆时偏移成像方面的优越性,同时验证了该方程数值求解的稳定性.

图 6 逆冲模型偏移结果,Zhou方程RTM成像结果(a),新qP方程RTM成像结果(b) Fig. 6 RTM of Thrust model, Zhou RTM imaging results (a), New pseudo acoustic RTM imaging results (b)
2.2 HESS模型数值模拟和逆时偏移

为了验证逆时偏移中考虑各向异性可行性和必要性,利用本文推导一阶拟声波方程对HESS模型进行试验,如图 7所示,分别是HESS模型的纵波速度场,ε模型,δ模型,其中P波速度范围为1500 m/s到4500 m/s,εδ的范围分别为0到0.28和0到0.16,震源为主频25赫兹的雷克子波,总共300炮,炮间距:30 m,每炮排列数:600,检波器间隔:10 m,时间采样长度:4.1 s,时间采样间隔:0.4 ms,空间采样间隔:10 m.采用震源归一化互相关成像条件以及Laplace滤波.

图 7 HESS模型Vp0(a)、ε(b)、δ(c)模型 Fig. 7 The Vp0(a)、ε(b) and δ(c) of Hess model

图 8为Zhou推导的qP方程和新导的qP方程正演结果(其中震源在模型中心位置),由图可知在HESS模型中新导的方程和Zhou方程的qP波波前都能很好的保持运动学特征,但在各向异性参数εδ变化加大区域新导的方程数值稳定性要好于Zhou方程.

图 8 t=1500 ms波场快照 Zhou方程(a)、新qP方程(b) Fig. 8 The wave field snapshot of t=1500 ms, Zhou's equation (a), New pseudo acoustic equation(b)

图 9分别为常规声波逆时偏移、拟声波逆时偏移成像结果.由图可知常规逆时偏移构造形态发生错乱,深部能量弱,成像效果差.VTI逆时偏移对该模型断层构造进行了很好的刻画,高陡倾角断层面处同相轴清晰,地下盐丘轮廓明显,尤其是大倾角地区,保幅成像效果更加明显.

图 9 HESS模型偏移结果,常规声波RTM成像结果(a),新qP方程RTM成像结果(b) Fig. 9 RTM of HESS model, acoustic RTM imaging results (a), New pseudo acoustic RTM imaging results (b)
3 结论与认识

本文基于弹性波的基本理论,推导出各向异性介质拟声波一阶速度-应力方程.文中从能量守恒的角度出发,证明了新推导的方程比基于频散关系推导出的拟声波方程更加稳定且能量保持守恒.通过对均匀二维模型和逆冲模型进行正演模拟,表明新推导的方程不仅能很好的维持qP波的运动学特征,而且能够很好的适应于任意旋转角度的TTI模型,而其他控制方程在都会出现不同程度的数值不稳定现象.本文将该各向异性介质拟声波一阶速度-应力方程引用到叠前逆时偏移技术中,引入归一化互相关成像条件.模型试算表明,各向异性介质拟声波一阶速度-应力方程逆时偏移能够校正介质的各向异性影响,成像剖面更加清晰,深部能量加强,整体能量更均衡,成像效果更加符合真实地层的构造特征.

致 谢 作者感谢两位审稿专家的宝贵意见,感谢国家973项目(2011CB202402),国家自然科学基金(41104069,41274124),山东省自然科学基金(ZR2011DQ016).

附录A1 从能量角度分析本文一阶速度-应力方程稳定性

由文中的定义,我们定义弹性波的能流密度为(单位质量的质点所携带的能量)

其中σε分别为应力矢量[σ11σ22σ33σ23σ13σ12]和应变矢量[ε11ε22ε33ε23ε13ε12],.

1.1 对势能项求时间导数

将(8)式应力应变关系改写为

将(A-2)带入到势能项中有:

从而我们可以将势能表达式写为

对(A-4)求时间导数:

其中(A-5)中应力更新方程可以从文中(8)得到:

将(A-6)带入(A-5)中:

对上式求分部积分:

我们设边界为无穷远,即在边界上vxvyvzσVσH的值都趋向于零值,将(A-7)改写为

1.2 对动能项求时间导数

将(8)中的速度更新方程带入(A-8)可到到:

由(A-10)和(A-8)可知,现在我们证明了本文推导的一阶速度-应力方程组在传播过程中能量是守恒.

2 从能量角度分析基于频散关系导出的拟声波方程的稳定性

基于频散关系推导拟声波方程主要有:Alkhalifah(1998)Zhou等(2006)Du等(2008)和Fowler(2009)的方程.下面于较为常用的Du等(2008)方程为例.

2.1 对势能项求时间导数

同理可得:

对上式求时间导数:

参考文献
[1] Alkhalifah T. 1998. Acoustic approximations for processing in transversely isotropic media[J]. Geophysics, 63(2): 623-631.
[2] Alkhalifah T. 2000. An acoustic wave equation for anisotropic media[J]. Geophysics, 65(4): 1239-1250.
[3] Cheng J B, Kang W. 2014. Simulating propagation of separated wave modes in general anisotropic media, part I: qP-wave propagators[J]. Geophysics, 79(1): C1-C18.
[4] Cheng J B, Kang W, Wang T F. 2013. Description of qP-wave propagation in anisotropic media, Part I: Pseudo-pure-mode wave equations[J]. Chinese Journal of Geophysics (in Chinese), 56(10): 3474-3486, doi: 10.6038/cjg20131022.
[5] Dellinger J, Etgen J. 1990. Wave-field separation in two-dimensional anisotropic media[J]. Geophysics, 55(7): 914-919.
[6] Du Q Z, Qin T. 2009. Multicomponent prestack reverse-time migration of elastic waves in transverse isotropic medium[J]. Chinese Journal of Geophysics (in Chinese), 52(3): 801-807.
[7] Du X, Fletcher R P, Fowler P J. 2008. A new pseudo-acoustic wave equation for VTI media[C].//70th EAGE Conference and Exhibition incorporating, EAGE. Extended Abstracts, H033.
[8] Fletcher R P, Du Xiang, Fowler P J. 2009. Reverse time migration in tilted transversely isotropic media[J]. Geophysics, 74(6): WCA179-WCA187.
[9] Huang C, Dong L G. 2009. Staggered-grid high-order finite-difference method in elastic wave simulation with variable grids and local time-steps[J]. Chinese Journal of Geophysics (in Chinese), 52(11): 2870-2878, doi: 10.3969/j.issn.0001-5733.2009.11.022.
[10] Kenneth P B, Nemeth T, Stefani J P, et al. 2012. On the instability in second-order systems for acoustic VTI and TTI media[J]. Geophysics, 77(5): T171-T186.
[11] Li B, Li M, Li H W, et al. 2012. Stability of reverse time migration in TTI media[J]. Chinese Journal of Geophysics (in Chinese), 55(4): 1366-1375, doi: 10.6038/j.issn.0001-5733.2012.04.032.
[12] Li Z C, Li N, Huang J P, et al. 2013. The quantity study of the factors that influence the Shear-wave splitting time in Fracture Media[J]. Chinese Journal of Geophysics (in Chinese), 28(1): 240-249, doi: 10.6038/pg20130125.
[13] Pestana R C, Ursin B, Stoffa P L.. 2011. Separate P-and SV-wave equations for VTI media[C].//SEG Technical Program Expanded Abstracts. Extended Abstracts, 163-167.
[14] Shen M C, Li M L, Zhou H L. 2014. Numerical simulation of decoupled VTI media acoustic wave equation[J]. Progress in Geophysics (in Chinese), 29(3): 1218-1223, doi: 10.6038/pg20140330.
[15] Wang J, Li Z C. 2014. Vector wave imaging and wave separation in the TI media[J]. Progress in Geophysics (in Chinese), 29(6): 2754-2760, doi: 10.6038/pg20140642.
[16] Yan J, Sava P. 2009. Elastic wave-mode separation for VTI media[J]. Geophysics, 74(5): WB19-WB32.
[17] Yan J, Sava P. 2011. Improving the efficiency of elastic wave-mode separation for heterogeneous tilted transverse isotropic media[J]. Geophysics, 76(4): T65-T78.
[18] Zang H, Li Z C. 2011. Seismic wave simulation method based on dual-variable grid[J]. Chinese Journal of Geophysics (in Chinese), 54(1): 77-86, doi: 10.3969/j.issn.0001-5733.2011.01.009.
[19] Zhan G, Pestana R C, Stoffa P L. 2011. An acoustic wave equation for pure P wave in 2D TTI media[C].//SEG Technical Program Expanded Abstracts. Extended Abstracts, 168-173.
[20] Zhan G, Pestana R C, Stoffa P L. 2012. Decoupled equations for reverse time migration in tilted transversely isotropic media[J]. Geophysics, 77(2): T37-T45.
[21] Zhang Y, Zhang H Z, Zhang G Q. 2011. A stable TTI reverse time migration and its implementation[J]. Geophysics, 76(3): WA3-WA11.
[22] Zhou H B, Zhang G Q, Bloor R. 2006. An anisotropic acoustic wave equation for modeling and migration in 2D TTI media[C].//SEG Technical Program Expanded Abstracts. Expanded Abstracts, 194-198.
[23] Zhou W, Guo Q S, Liu X Y, et al. 2014. The influence of Thomsen parameters on Kirchhoff prestack depth migration in VTI media[J]. Progress in Geophysics (in Chinese), 29(6): 2866-2873, doi: 10.6038/pg20140657.
[24] 程玖兵, 康玮, 王腾飞. 2013. 各向异性介质qP波传播描述I: 伪纯模式波动方程[J]. 地球物理学报, 56(10): 3474-3486, doi: 10.6038/cjg20131022.
[25] 楚泽涵, 张颖. 1998. 横波在各向异性复合介质中的传播[J]. 中国石油大学学报(自然科学版), 22(1): 21-24. (请补充英文文献信息)
[26] 杜启振, 秦童. 2009. 横向各向同性介质弹性波多分量叠前逆时偏移[J]. 地球物理学报, 52(3): 801-807.
[27] 黄超, 董良国. 2009. 可变网格与局部时间步长的交错网格高阶差分弹性波模拟[J]. 地球物理学报, 52(11): 2870-2878, doi: 10.3969/j.issn.0001-5733.2009.11.022.
[28] 李博, 李敏, 李红伟,等. 2012. TTI介质有限差分逆时偏移的稳定性探讨[J]. 地球物理学报, 55(4): 1366-1375, doi: 10.6038/j.issn.0001-5733.2012.04.032.
[29] 李振春, 李娜, 黄建平,等. 2013. 裂缝介质横波分裂时差影响因素定量研究[J]. 地球物理学进展, 28(1): 240-249, doi: 10.6038/pg20130125.
[30] 沈铭成, 李录明, 周怀来. 2014. VTI介质解耦合声波方程数值模拟[J]. 地球物理学进展, 29(3): 1218-1223, doi: 10.6038/pg20140330.
[31] 王娟, 李振春. 2014. TTI介质矢量波成像及波型分解[J]. 地球物理学进展, 29(6): 2754-2760, doi: 10.6038/pg20140642.
[32] 张慧, 李振春. 2011. 基于双变网格算法的地震波正演模拟[J]. 地球物理学报, 54(1): 77-86, doi: 10.3969/j.issn.0001-5733.2011.01.009.
[33] 周巍, 郭全仕, 刘旭跃,等. 2014. Thomsen参数对VTI介质克希霍夫叠前深度偏移的影响[J]. 地球物理学进展, 29(6): 2866-2873, doi: 10.6038/pg20140657.