王璐琛，常旭，王一博

1 中国科学院页岩气与地质工程重点实验室，中国科学院地质与地球物理研究所，北京　100029　2　中国科学院大学，北京　100049



TTI介质的交错网格伪P波正演方法

王璐琛1,2，常旭1*，王一博1

1 中国科学院页岩气与地质工程重点实验室，中国科学院地质与地球物理研究所，北京1000292中国科学院大学，北京100049

摘要研究了三维弱各向异性近似下，利用伪P波(伪纵波)模拟弹性波场P分量在倾斜对称轴的横向各向同性(TTI)介质中的传播过程，并对比了分别基于弹性Hooke定律、弹性波投影和运动学色散方程所建立的三种二阶差分伪P波方程的正演特点.目前这些伪P波方程数值计算主要采用规则网格差分，但是规则网格在TTI模拟中有低效率、低精度以及不稳定的缺点.为了提高计算的精度，本文构建出相应方程的交错网格有限差分格式.通过对比伪P波方程在三维TTI介质中不同的数值模拟的表达形式，本文认为基于色散方程所建立的伪P波方程在模拟弹性波中P波传播的过程中具有最小的噪声.本文分析不同的各向同性对称轴空间角度的频散特征，并引入适当的横波速度维持计算的稳定.二维模型算例表明，本文提出的交错网格正演算法可以得到稳定光滑的伪P波正演波场.使用本文交错网格算法对二维BP TTI模型的逆时偏移也具有较稳定的偏移结果.

关键词伪P波； 倾斜各向异性； 交错网格

1引言

目前油气开采中重点研究的页岩气主要聚集吸附在致密页岩层或泥岩层中.研究表明页岩的各向异性可达70%(Vernik and Nur,1992; Johnston and Christensen,1995)，因此使用倾斜对称轴的横向各向同性(TTI)介质(Qin et al.,2003a; 2003b; 吴国忱等,2010; 黄翼坚等,2011)，能更好地描述地震波在页岩中的传播规律.页岩气的开采需要利用被动源微地震进行监测.被动源地震的震源不同于主动源地震，其震源存在多种破裂形式.这些破裂形式受到所处地下环境的介质各向异性所影响，尤其是在最近研究比较热门的微地震领域，其震源基本都是处在很明显的TTI沉积层介质中.研究震源处介质的各向异性有助于了解地震区域震源机制的构成.

由于弹性波的传播本身已经比较复杂，介质的各向异性又会增加非零的Hooke弹性参数从而加重这种复杂的状况，会大幅度增加波场正演的计算时间以及计算量.同时考虑到纵波勘探应用广泛，技术也更加成熟，是比较常用的研究方法.因此，只研究类似P波的波场传播(称为伪P波)，无疑是在各向异性研究中一种比较折中的方式(Ekkehart,2010; Yoon et al.,2010).针对伪P波的波场传播，可以通过限定不同的约束条件(Alkhalifah,1998; Zhou et al.,2006; Fletcher et al.,2008)或不同的分离方式从TTI弹性波方程中建立出伪P波的传播方程(Fletcher et al.,2009; Zhan et al.,2011; Kang and Cheng,2011).

这些方程普遍都是复杂的二阶微分方程，所以在数值实现上基本采用最直观的离散方式，用规则网格差分将二阶方程进行数值离散.然而规则网格由于差分精度有限，一方面不利于在有限的存储与计算时间下提高计算精度；另一方面在各向异性参数变化较为剧烈的区域，由于差分网格点间距过大，容易造成计算的不稳定性.同时有些规则网格算法会将各向异性参数的微商引入导数项的计算中，这样容易因各向异性参数微商计算的精度不够而造成数值计算发散.

由于规则网格差分在TTI模拟中存在如上弊端，构造出结构简约的交错网格差分格式就显得尤为必要.首先，交错网格本身的计算精度就要高于同样差分阶数的规则网格，另外由于交错网格是在半网格点上进行差分的，能更好的适应各向异性的复杂变化，所以数值计算的稳定性也能提高.

本文主要研究在弱各向异性(Thomsen,1986)近似下，地震波传播交错网格差分的实现方法，分析了TTI介质中交错网格相比于规则网格所具有的高效、低存储等计算优势，建立了统一的交错网格伪P波差分方式.本文对比了三种主流的伪P波正演方程，并分别建立了交错网格的离散公式.本文的交错网格差分方法将各向异性的旋转特性加载在空间导数上，这样一方面保留了各向异性的物性参数，便于提取物理属性，方便研究地层的岩性(杜启振和秦童,2009)，另一方面避免了交叉导数项的出现，简化了计算，提高了计算稳定性.本文在研究TTI介质的性质时是先假定TTI倾角(θ)为零，然后再推广到任意的空间转角，因此本文交错网格差分方程的适应性比较好.

2三维TTI介质的波场传播理论

地震波在各向异性介质中的传播是一个复杂的弹性波传播过程，各向异性会使得Hooke弹性系数矩阵变得更加复杂，同时也不利于波场模拟.由于弹性各向异性介质波动方程的复杂，所以大多数各向异性正演和偏移方法通常都是在Thomsen弱各向异性的近似下使用伪P波方程来实现.

2.1TTI介质的弱各向异性

弹性波波动方程的一般数学表示为

(1)

其中，u通常代表质点的位移场，ρ和F分别表示介质的密度和所受到的外力.cijkl表示4阶Hooke定律中的刚度矩阵，表示应力与应变的对应关系.由于应力的对称性，所以空间上9个应力分量中只有Pxx,Pyy,Pzz,Pxy,Pyz,Pzx，这6个是独立的，分别表示3个正应力以及3个剪切应力.又由于弹性系数张量的对称性，所以刚度矩阵中一共只有21个是独立的变量.刚度矩阵具体的形式如下：

(2)

如果是横向各向同性的情况，弹性刚度矩阵则化简如公式(3)所示：

(3)

式(3)中的弹性参数，无法直接对应物理量，因此按照Thomsen引入的5个弱各向异性弹性参数(Thomsen,1986)进行如下替换，可赋予各参数可测的物理意义：

(4)

(5)

2.2弱各向异性下的P波运动学色散方程

弱各向异性下的伪P波方程，并不是简单的在弹性波动方程中令横波波速为0.因为各向异性本身就是对弹性的一种描述，不相等的Thomsen弱各向异性参数(例如ε≠δ，此时称为非椭圆各向异性)(Thomsen,1986)会导致方程在不同方向的传播上出现不同的横波假象.所以如设定了横波波速为0，还是会产生伪横波的出现(Tsvankin et al.,2001).同时根据Alkhalifah提出的基于横波速度为0的纵波各向异性方程(Alkhalifah,1998; 2000)，被假定为0的也仅仅是沿着TTI介质中各向同性对称轴方向的横波速度.虽然纵波的传播在运动学上基本满足各向异性传播性质，但是不同的推导方式还是会有不同的振幅以及频散表现，尤其是在不同的TTI对称轴的空间角度上，表现会更加剧烈.

首先研究VTI介质的情况.在各向异性的情况下，波的相速度能够描述波在介质中传播的特点，因此可以从相速度的角度建立波场的运动学传播公式.先求解Christoffel方程可以得到精确的P-SV波和SH波相速度公式(Dellinger and Etgen,1990; Tsvankin.,1996;2001)，其中P-SV波的相速度公式如下：

(6)

(7)

2.3伪P波方程组构造

(1)依据弹性Hooke定律：由公式(3)和(5)所表示的弹性Hooke方程进行变量代换(Zhang et al.,2011)，设

(8)

其中，ux，uy，uz表示三维空间三个方向上的位移.将公式(8)代入TTI弹性Hooke方程就可以得到如下的伪P波时域方程组：

(9)

(10)

(11)

那么就可以得到如下的运动学伪P波时域方程组(Zhang et al.,2011; Fletcher et al.,2009; Zhan et al.,2011; 张岩和吴国忱,2013)：

(12)

这里要强调的是，之所以称是伪P波，因为简单的令沿着TTI对称轴的横波速度为0并不能将介质其他部分的横波相速度消除(Grechka et al.,2004).当然可以通过辅助波场压制伪横波传播，或者通过震源处理(Duveneck et al.,2008)减小伪横波的影响.

(13)

(14)

这样就得到了TTI的伪P波方程.其中包含一个通常意义上的声波压力场和一个计算引入的辅助场，本质上是对TTI弹性波方程的一种计算上的简化.

2.4伪P波方程交错网格差分格式

公式(9)、(10)、(12)均给出的是二阶位移方程，按照公式格式进行数值离散即可得到相应的VTI数值方程，进一步借由公式(14)旋转得到TTI数值方程，这样就可以实现公式的规则网格差分格式.这样处理的一个好处是，从VTI过度到TTI所需的空间旋转计算被加载在了变量的空间二阶微分上，避免了对各向异性参数的空间旋转插值，提高了计算的稳定性以及精确性，缺点就是二阶微分中包含了交叉导数项，不利于高精度的计算.为了提高计算精度，现在考虑如何实现公式(9)、(10)、(12)的应力、速度交错网格差分格式.这里以二维情况的公式(12)为例说明.二维方程与空间导数的旋转形式如公式(15)和(16)所示：

(15)

(16)

对比声波方程交错网格格式，记速度场vpx,vpz,vrx,vrz满足

(17)

(18)

(19)

(20)

(21)

速度场v rx,v rz空间导数计算同理为

(22)

代入公式(22)到公式(19)中，得旋转后的速度场空间导数差分计算为

(23)

所以公式(18)所对应的应力方程差分格式由上式代入可得

(24)

公式(20)和公式(24)合在一起就是公式(12)所对应的方程的一阶速度应力方程交错网格差分格式.我们很容易发现，公式(12)的方程形式仅仅体现在交错网格差分格式的应力方程(24)中，这样的目的在于可以最大限度的简化速度场方程的差分形式，同时也易于扩展到其他形式的伪P波方程.因此，公式(9)和(10)所对应方程的交错网格速度方程同公式(20)完全一样，应力方程同理公式(24)其格式分别如下:

(25)

和

(26)

这样就构建出了三种方法所对应方程(公式(9)(10)(12))的交错网格差分格式.其中对TTI对称轴的空间旋转效应的计算仅加载在了应力方程中的速度场空间导数差分上.

3伪P波在TTI介质中的数值计算

3.1伪P波的正演

先考察TTI介质中不同正演差分策略的计算效率.以二维均匀模型正演计算为例，计算区域网格数为400×400，正演步数为3000步，对伪P波交错网格、伪P波规则网格两种正演策略各进行40次计算，结果显示于表1.

作为对比，对同样大小的速度模型采用声波交错网格和声波规则网格两种正演策略各进行相同的40次正演计算，结果也显示于表1.

表1　二维模型下不同正演方式的内存需求和计算效率

从表1可以看出，对于伪P波的正演模拟，交错网格方法无论是在内存需求上还是计算时间上都要少于规则网格算法，这一点同传统声波方程的有限差分是不同的.对于声波方程而言，交错网格的计算精度比同阶差分的规则网格精度要高，但是计算效率要慢，这是因为一阶速度应力方程比二阶声波方程需要更多的差分计算.由于TTI伪P波的规则网格正演方程需要计算二阶空间微分的交叉导数，因此对每个网格点需要额外的一次循环计算(先计算一个方向的一阶差分，再计算另一个方向上的差分)，并需要引入相应的交叉导数临时变量.这样即增加了额外计算量导致总的计算效率下降，以至于平均计算时间反而超过TTI伪P波交错网格，又增大了内存需求(三维情况下计算的交叉导数项更多，对应需要更多的临时变量).另外由于交错网格差分精度比同阶的规则网格精度要高，因此TTI伪P波交错网格算法相比于规则网格算法有着更优越的应用效果.

再考察TTI不同正演方程的模拟效果.通常情况下降阶得出的二元二阶耦合方程中，辅助波场可以不具有实质上的物理含义，也就是说在此情况下并不会影响到P波运动学的特性.而如公式(10)的方程(Kang and Cheng,2011)，又是从各向异性弹性波运动方程中推导的P波方程，分别表示两个应力方向上的波场传播(二维情况下为沿着倾斜对称轴的方向和垂直于倾斜对称轴的方向)，其波场快照如图1a所示(使用了声波近似条件，即vSz=0).其中P波速度vPz=2500 m·s-1，各向异性参数和TTI对称轴的空间角度分别是ε=0.24，δ=0.1，θ=45°.将其二者相加组合就是完整的P波波场.

可以看到，在图1中各个分量波场中((a1—c1)和(a2—c2)图)内部尖角所表示的伪S波正好呈现振幅相反的特点，相加后的波场中(a3—c3图)，最外层的P波成分变得更加完整清晰，同时内部的伪S波尖角也得到了衰减.所以这里就将耦合公式中的两个分量进行相加来表示模拟的伪P波波场.

将此处理应用到公式(9)，公式(12)所表示的三种伪P波耦合方程，二维情况下的各自的波场快照对比结果如图1B和图1C，其中P波速度vPz=2500 m·s-1，各向异性参数和TTI对称轴的空间角度分别是ε=0.24，δ=0.1，θ=45°，同样使用声波近似条件，即令vSz=0.

对比图1B、1C和图1A，不难看出，不管采用何种公式，辅助波场r都体现着更多的内部伪S波成分，对于外圈的P波成分，或者没有实际意义(图1B)，或者起到补偿消减的作用(图1C).所以说明，将两个波场组合相加，可以提高伪P波波场的精确性.由于两个波场所表示的主体不同，辅助场主要体现的是次生的伪S波场成分，所以对这两个波场就可以进行各自的滤波处理，从而进一步压制伪S波的成分(张岩和吴国忱,2013).

分析图1以及对应的公式(公式9、10、12)，可 以看出，图1A(公式10)是按照波场运动的方向进行波场分解的，两个分量场表示伪P波波场两个方向的运动分量.图1B(公式9)是按照动力学应力场的方向分量进行波场分解的，两个分量场表示伪P波波场的两个方向的应力分量.其波场反映着动力学的振幅信息.图1C(公式12)是按照波场的P、SV成分进行分解的，两个分量场分别表示P波场和SV波场.图2表示的是三个伪P波公式(公式9、10、12)正演合成波场的振幅和频散对比：

图1　三种正演公式模拟的波场快照(A)公式(10)： (B)公式(9)； (C)公式(12)； (a1—c1)p成分; (a2—c2)r成分; (a3—c3)p+r成分

从图2A中可以明显看出公式(10)与公式(12)对应的伪P波波场是比较相似的，而基于应力的公式(9)对应的绿色实线同另外两个相差较大，无论是波形的形状还是到时都存在较明显的差异.它反映了TTI介质中质点的应力场同速度场之间的差异.图2B中可看出，频散主要集中在波场中间的伪横波尖角处.并且公式(9)的频散明显要强于另外两个公式正演的波场.所以本文之后的模型正演部分，采用的是频散较小的公式(12)构造的伪P波方程.同时采用有约束的横波速度(Tsvankin et al.,2001)(不是简单设为0)，以及特征投影滤波的方式(Zhou and Zhang,2009)压制正演中的伪横波噪声，可以进一步降低数值频散，提高计算的稳定性.设vPz=2500 m·s-1，ε=0.24，δ=0.1，图3展示了伪横波数值频散的压制效果.其中外圈椭圆为正演的伪P波波前，内部波场为伪横波噪声.

除了不同降阶公式本身的影响，TTI介质中伪P波的传播还会受到对称轴角度的影响.这种影响主要是源于数值计算上产生的频散现象，角度越倾斜，频散的现象越严重，尤其是在三维模拟的情况下.令θ表示TTI对称轴相对于竖直的倾角，φ表示TTI对称轴的方位角，图4所示的是对称轴在不同空间转角下的波场垂直切片.

可以看到，在空间转角为0°的时候，数值频散现象是非常弱的，但是随着转角的出现，在内部伪S波尖角的地方还是出现了频散现象，右图的空间角度更加复杂，频散的现象也比中间的要更加复杂.

3.2模型正演与偏移结果

首先考察二维情况下本文所述方法的正演效果.这里作为简化，我们只研究固定TTI对称轴倾角的情况.所用模型如图5所示.

模型纵向深2000 m，横向长4000 m，网格大小为4 m，其中TTI各向异性仅仅存在于中间的两层，也就是说最上层和最下层是普通的各向同性介质，中间两层中TTI的对称轴角度设定为固定值 45°.此模型表示一个斜插入的各向异性层，并且该层内还具有构造起伏的特征.

图2　不同公式合成波场的道记录(A)及不同公式合成波场的频散(B)(a1) 第101道;(a2) 第52道; (b1) 公式(9);(b2) 公式(10); (b3) 公式(12).

图3　二维TTI正演波场中的伪横波(a) 横波速度vPz为0; (b) 横波速度vPz不为0; (c) 横波速度vPz不为0并使用滤波.

图4　三维空间中对称轴不同转角下的垂向波场(a) 转角:θ=0°,φ=0°; (b) 转角:θ=45°,φ=0°; (c) 转角:θ=45°,φ=30°.

图5　TTI各向异性参数模型(a) P波速度vPz； (b) δ； (c) ε； (d) θ.

采用伪P波模拟的单一地下震源激发的不同时刻地震波的波场快照如图6所示.其震源位置为横向2000 m，深度1000 m处.图6中从左至右的波场时刻依次是0.2、0.3、0.4、0.6 s.

从图6B中可以发现，在非各向异性介质中，辅助波场直接衰减为0，从另一方面也说明了本文中使用的伪P波方程的辅助场描述的是伪S波的产生与传播.因此也就是说，合成波场在地表的检波器处同伪P波波场是相同的，波场合成并不影响最终的道记录.如果说存在影响的话，影响也仅会发生在TTI介质的井中接收情况中.此地下震源对应的伪P波单炮记录如图7所示，从图中可见模拟的地震记录波形稳定，并且没有计算频散与发散.说明我们的正演效果较好.

对于三维情况，速度模型与各向异性参数模型如图8所示(只显示非恒定值的量)，采用前面图5的二维4层速度模型在Y方向上平推而得，因此三维模型的中间两层TTI的对称轴角度为固定值θ=45°，φ=0°.地面观测系统如图9a所示，用水平地表“米”字型共4条测线进行观测.模型大小为X方向4000 m，Y方向800 m，深度2000 m，网格大小4 m.

假定震源位置为X方向2000 m，Y方向400 m，深度1000 m处.此震源事件的单炮记录如图9b所示.从图9b中可见，三维TTI情况下，本文正演的波场也具有良好的数值稳定性.

最后再用标准TTI模型的逆时偏移(RTM)来进一步测试本文的正演算法.RTM偏移采用炮点沿时间正向延拓震源子波，检波点同时沿时间反向延拓地震记录的方式，得到炮点与检波点波场.接着利用成像条件与这两个波场对地下反射点进行成 像.其中，波形正确、稳定、频散小的正演差分求解，是RTM偏移成像的基础.BP TTI模型是国际上普遍采用的二维各向异性模型，其模型网格大小为6.25 m×6.25 m，横向距离78.7 km，最大深度11.25 km，我们用模型横向0～51 km的部分进行计算.其速度的模型与各向异性参数模型如图10所示，使用发布的数据以及本文所述的正演方式对其进行偏移的结果如图11所示.

图6　(A) 单一震源的合成波场(p+r)和(B)辅助波场(r)在不同时刻的波场从左至右时间依次是：0.2 s，0.3 s，0.4 s，0.6 s.

图7　单震源炮记录

图8　(a) 三维P波vPz速度模型； (b) TTI各向异性模型δ； (c) TTI各向异性模型ε

图9　三维TTI模型的地表观测系统(a)及对应的地震事件的炮记录(b)

图10　BP TTI模型的速度与各向异性参数(a) vPz； (b) ε； (c) δ； (d) θ.

图11　使用本文正演方式的逆时偏移(RTM)结果

图11偏移结果同相轴连续性较好，并且偏移结果的盐丘位置同给定模型中的盐丘位置能相互匹配，说明本文的正演算法能正确的模拟TTI介质中地震波波形.另外，成像结果无发散，说明本文的正演差分算法具有良好的稳定性.

4结论

在现阶段对复杂地质的地震资料的处理和解释过程中，忽略各向异性将会导致很大的偏差，无法得到精度更高的结果.介质的各向异性，尤其是更贴近实际情况的TTI介质，能够显著地改善资料处理的质量.对于TTI介质的研究，首要的任务是估计正确的各向异性参数，然后是使用恰当的方式模拟波场的传播.本文所使用的伪P波方程能够简化TTI波场的复杂性，并能较快捷的应用到三维的成像、偏移中.这其中要注意的是，不同的伪P波方程具有不尽相同的物理意义，同时也有不同的频散现象.本文对比了基于位移场和应力场的伪P波方程各自的特点，推导并实现了正演模拟三维各向异性介质中P波成分的伪P波交错网格有限差分算法.从正演的波场上看，基于色散方程降阶分离出的伪P波波场频散较小，正演模拟效果最好.最后，由于频散在TTI介质中的影响会远远大于各向同性介质，因此如何有效的去除频散，认识频散在不同方程条件下的物理意义，仍将是后续研究的重点.

致谢感谢国家自然科学基金资助(41230317、41274112)，感谢英国石油公司(BP)提供的二维TTI理论模型与其合成数据.

References

Alkhalifah T.1998.Acoustic approximations for processing in transversely isotropic media.Geophysics,63(2): 623-631.

Alkhalifah T.2000.An acoustic wave equation for anisotropic media.Geophysics,65(4): 1239-1250.

Dellinger J,Etgen J.1990.Wave-field separation in two-dimensional anisotropic media.Geophysics,55(7): 914-919.Du Q Z,Qin T.2009.Multicomponent prestack reverse-time migration of elastic waves in transverse isotropic medium.Chinese J.Geophys.(in Chinese),52(3): 801-807.

Duveneck E,Milcik P,Bakker P M,et al.2008.Acoustic VTI wave equations and their application for anisotropic reverse-time migration.∥ 78th Annual International Meeting,SEG,Expanded Abstracts,2186-2190.

Ekkehart T.2010.Reverse-time migration in TTI media.∥ 80th Annual International Meeting,SEG,Expanded Abstracts,3193-3197.

Fletcher R,Du X,Fowler P J.2008.A new pseudo-acoustic wave equation for TI media.∥ 78th Annual International Meeting,SEG,Expanded Abstracts,2082-2086.

Fletcher R P,Du X,Fowler P J.2009.Reverse time migration in tilted transversely isotropic (TTI) media.Geophysics,74(6): WCA179-WCA187.

Fowler P J,Du X,Fletcher R P.2010.Coupled equations for reverse time migration in transversely isotropic media.Geophysics,75(1): S11-S22.Grechka V,Zhang L B,Rector III J W.2004.Shear waves in acoustic anisotropic media.Geophysics,69(2): 576-582.

Hao C T,Yao C,Wang X.2006.The characteristics of velocities with azimuth variation for arbitrary spatial orientation TI media.Progress in Geophysics (in Chinese),21(2): 524-530.

Huang Y J,Zhu G M,Liu C Y.2011.An approximate acoustic wave equation for VTI media.Chinese J.Geophys.(in Chinese),54(8): 2117-2123,doi: 10.3969/j.issn.0001-5733.2011.08.019.

Johnston J E,Christensen N I.1995.Seismic anisotropy of shales.J.Geophys.Res.,100(B4): 5991-6003.

Kang W,Cheng J B.2011.Prestack scalar reverse time migration of elastic seismic data in TI media.∥ 81th Annual International Meeting,SEG,Expanded Abstracts,3367-3371.

Qin Y L,Zhang Z J,Li S L.2003a.CDP mapping in tilted transversely isotropic (TTI) media.Part I: Method and effectiveness.Geophysical Prospecting,51(4): 315-324.

Qin Y L,Zhang Z J,Xu S H.2003b.CDP mapping in tilted transversely isotropic (TTI) media.Part II: Velocity analysis by combining CDP mapping with a genetic algorithm.Geophysical Prospecting,51(4): 325-332.

Thomsen L.1986.Weak elastic anisotropy.Geophysics,51(10): 1954-1966.

Tsvankin I.1996.P-wave signatures and notation for transversely isotropic media: An overview.Geophysics,61(2): 467-483.

Tsvankin I.2001.Seismic Signatures and Analysis of Reflection Data in Anisotropic Media.Amsterdam: Pergamon.

Vernik L,Nur A.1992.Ultrasonic velocity and anisotropy of hydrocarbon source rocks.Geophysics,57(5): 727-735.

Wu G C,Liang K,Yin X Y.2010.The analysis of phase velocity and polarization feature for elastic wave in TTI media.Chinese J.Geophys.(in Chinese),53(8): 1914-1923,doi: 10.3969/j.issn.0001-5733.2010.08.017.

Yoon K,Sang S,Jean J,et al.2010.Stability and speedup issues in TTI RTM implementation.∥ 80th Annual Meeting,SEG,Expanded Abstracts,3221-3225.

Zhan G,Pestana R C,Stoffa P L.2011.An acoustic wave equation for pure P wave in 2D TTI media.∥ 81th Annual International Meeting,SEG,Expanded Abstracts,168-173.

Zhang Y,Wu G C.2013.Methods of removing pseudo SV-wave artifacts in TTI media qP-wave reverse-time migration.Chinese J.Geophys.(in Chinese),56(6): 2065-2076,doi: 10.6038/cjg20130627.

Zhang Y,Zhang H Z,Zhang G Q.2011.A stable TTI reverse time migration and its implementation.Geophysics,76(3): WA3-WA11.

Zhou H B,Bloor R,Zhang G Q.2006.An anisotropic acoustic wave equation for modeling and migration in 2D TTI media.∥ 76th Annual International Meeting,SEG,Expanded Abstracts,194-198.

Zhou H Z,Zhang G Q,Zhang Y.2009.Removing S-wave noise in TTI reverse time migration.∥ 79th Annual International Meeting,SEG,Expanded Abstracts,2849-2853.

附中文参考文献

杜启振,秦童.2009.横向各向同性介质弹性波多分量叠前逆时偏移.地球物理学报,52(3): 801-807.

郝重涛,姚陈,王迅.2006.任意空间取向TI介质中速度随方位变化特征.地球物理学进展,21(2): 524-530.

黄翼坚,朱光明,刘池洋.2011.一个近似的VTI介质声波方程.地球物理学报,54(8): 2117-2123,doi: 10.3969/j.issn.0001-5733.2011.08.019.

吴国忱,梁锴,印兴耀.2010.TTI介质弹性波相速度与偏振特征分析.地球物理学报,53(8): 1914-1923,doi: 10.3969/j.issn.0001-5733.2010.08.017.

张岩,吴国忱.2013.TTI介质qP波逆时偏移中伪横波噪声压制方法.地球物理学报,56(6): 2065-2076,doi: 10.6038/cjg20130627.

(本文编辑汪海英)

基金项目国家自然科学基金(41230317、41274112)资助.

作者简介王璐琛，男，1989年生，在读博士，主要从事地震波成像，微地震的反演和震源机制研究.E-mail：wangluchen@mail.iggcas.ac.cn *通讯作者常旭，研究员，研究方向为地震波传播、地震成像、地球物理监测.E-mail：changxu@mail.iggcas.ac.cn

doi:10.6038/cjg20160325 中图分类号P631

收稿日期2014-12-03，2016-01-15收修定稿

Forward modeling of pseudo P waves in TTI medium using staggered grid

WANG Lu-Chen1,2，CHANG Xu1*，WANG Yi-Bo1

1KeyLaboratoryofShaleGasandGeoengineeringInstituteofGeologyandGeophysics,ChineseAcademyofSciences,Beijing100029,China2UniversityofChineseAcadamyofSciences,Beijing100049,China

AbstractA tilted transversely isotropic (TTI) medium is a good approximation for anisotropic problems.We discuss the property of a 3D tilted transversely isotropic (TTI) medium and use pseudo P waves to simulate the P component of the elastic wave field in the weak elastic anisotropy approximation.We compare forward features of 3 second-order difference equations of pseudo P waves based on Hooke′s law,elastic wave projection and dispersion equation,respectively.Now numerical calculation of these equations uses regular grid finite difference method,which is of low efficiency,low accuracy and is unstable in TTI simulation.In order to improve the calculating accuracy,this paper builds the staggered grid finite difference format for each forwarding simulation.Comparing three different numerical expressions of pseudo P wave equation in a 3D TTI medium,we find that the method based on dispersion equation is more convenient to simulate elastic P wave propagation and has lowest numerical noise.By analyzing the frequency dispersion characteristics based on different angles of isotropic symmetry axis,we introduce appropriate SV wave velocity in calculation to keep calculation stable.An example of a 2D simple synthetic model shows that the staggered grid forwarding method allow us to obtain a smooth and stable pseudo P wave field.With this strategy,2D BP TTI reverse time migration can also achieve basically stable migration results.

KeywordsPseudo P wave; Tilted transversely isotropic; Staggered grid

王璐琛，常旭，王一博.2016.TTI介质的交错网格伪P波正演方法.地球物理学报，59(3):1046-1058,doi:10.6038/cjg20160325.

Wang L C,Chang X,Wang Y B.2016.Forward modeling of pseudo P waves in TTI medium using staggered grid.Chinese J.Geophys.(in Chinese),59(3):1046-1058,doi:10.6038/cjg20160325.