目前油气开采中重点研究的页岩气主要聚集吸附在致密页岩层或泥岩层中.研究表明页岩的各向异性可达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倾角(θ)为零，然后再推广到任意的空间转角，因此本文交错网格差分方程的适应性比较好.

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

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

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

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

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

式中，ε是纵波各向异性强度的度量，ε越大，纵波在介质中的各向异性程度越强，ε=0表示纵波没有各向异性；γ是横波各向异性强度的参数，γ越大，横波在介质中的各向异性程度也越强，γ=0表示横波没有各向异性，γ同时也是横波分裂强度的度量；v_Pz与v_Sz分别是纵波与横波在垂直介质的各向同性面上的波场相速度；δ是描述纵波层内NMO速度相对强度的参数，表示纵波层内的NMO速度.横向各向同性的弹性矩阵共有6 个参数，其中有一个是非独立的，所以只需要Thomsen 引入的5个参数就可以完全表示横向各向同性介质的弹性关系.公式(4)可推导弹性参数如下(郝重涛等，2006)：

其中是表示横纵波速比的参数.这样就得到了TI介质中的弹性参数表达.从公式(3)不难看出，这种TI介质默认的是垂直方向对称轴的，也称为是VTI介质.由于TTI介质与VTI介质相比仅仅是对称轴的方向不同，TTI介质的对称轴会相对于空间垂直方向上存在转角.因此只要在方程中进行一次坐标变换即可，并不用改变刚度矩阵的参数.

弱各向异性下的伪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波的相速度公式如下：

其中，δ，ε是之前所述的Thomsen各向异性参数，v₀²是P波沿着对称轴的相速度，α表示P波同对称轴所成的相角.将(6)式的P-SV色散关系变换到频率波数域则写为

(7)式是4阶的色散公式，其中，.为了简化计算可以引入一个辅助波场，将(7)式降阶成两个二阶的公式.选择不同的辅助波场，就能表现出不同的二元二阶方程组(Fletcher et al.，2008; Fowler et al.，2010).

记H_xx，H_yy，H_zz为二阶的空间偏导数算子，为横纵波速平方比系数，分别表示P波的压力场，以及辅助波场，利用不同的原理可以构造如下三种不同的伪P波方程：

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

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

(2)依据弹性波投影：如果将弹性波投影到波数域中得到解耦合的P-SV波，可得到伪P波在P-SV平面上的两个方向分量方程，此时辅助场和压力场 同样具有实际的物理意义，二者表示伪P波的两个正交方向分量(Kang and Cheng，2011; Fowler et al.，2010)

(3)依据相速度色散方程：作为简化对色散方程进行降阶，需要在压力场的基础上引入辅助场，不同的辅助场设计会产生不同的分解方程(Zhou et al.，2006)，本文取

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

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

公式(9)、(10)、(12)中的二阶偏导数H_xx，H_yy，H_zz都是针对于基本坐标系而言的，也就是适用于VTI介质.由于TTI介质相对于VTI仅仅是多了空间上的转角，因此将这里的H_xx，H_yy，H_zz对应一个空间旋转映射就可以得到TTI介质中的伪P波方程.令TTI对称轴在空间中的倾角为θ，方位角为，则坐标旋转后的坐标为

所以二阶偏导数H_xx，H_yy，H_zz为

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

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

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

代入公式(15)的二阶方程，并结合(16)可以得到应力场 ， 为

其中一阶速度场空间差分为

公式(17)和(18)就是可以用于交错网格的一阶速度应力方程.这里TTI空间旋转的计算是加载在速度场变量的空间微分上(如公式(19)所示)，因此也不需要对各向异性参数进行空间旋转插值.另外此时的速度场微分已经是一阶微分，所以不存在交叉导数项的计算.具体的差分离散格式如下：设变量的上标为时刻，下标为空间网格点位置，例如 _i，kⁿ表示变量场在n时刻，空间网格(i，k)点上的波场值；Δt表示时间差分步长；Δx，Δz表示不同方向上空间差分步长；又记f=v_Sz²/v_Pz²表示横纵波速度平方的比值，v²=v_Pz²表示纵波速度的平方，N表示差分精度为2(N+1)阶.

首先将速度场变量vpx，vpz，vrx，vrz定义在半网格点上，应力场，定义在整网格点上.根据速度场公式(17)得到速度场差分公式为

速度场vpx，vpz空间导数计算为

速度场vrx，vrz空间导数计算同理为

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

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

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

和

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

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

表 1(Table 1)

表 1 二维模型下不同正演方式的内存需求和计算效率 Table 1 Memory requirement and computational efficiency of different methods in 2D model 正演方式 内存需求/M 临时内存/M 平均计算时间/s 伪P波规则网格 15.36 5.12 311.64 伪P波交错网格 15.36 0 270.74 声波规则网格 7.68 0 56.19 声波交错网格 7.68 0 74.95

表 1 二维模型下不同正演方式的内存需求和计算效率 Table 1 Memory requirement and computational efficiency of different methods in 2D model

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

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

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

图 1

Fig. 1

图 1 三种正演公式模拟的波场快照 (A)公式(10)： (B)公式(9)； (C)公式(12)； (a1—c1)p成分; (a2—c2)r成分; (a3—c3)p+r成分Fig. 1 The snapshots generated by these 3 forward equations (A) Equation (10); (B) Equation (9); (C) Equation (12); (a1—c1) p component; (a2—c2) r component; (a3—c3) p+r component.

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

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

图 2

Fig. 2

图 2 不同公式合成波场的道记录(A)及不同公式合成波场的频散(B) (a1) 第101道;(a2) 第52道; (b1) 公式(9);(b2) 公式(10); (b3) 公式(12).Fig. 2 Record of combined wavefield with different equations (A) and Dispersion of combined wavefield with different equations (B) (a1) Trace No.101; (a2) Trace No.52 ; (b1) Equation 9; (b2) Equation 10; (b3) Equation 12.

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

图 3

Fig. 3

图 3 二维TTI正演波场中的伪横波 (a) 横波速度vPz为0; (b) 横波速度vPz不为0; (c) 横波速度vPz不为0并使用滤波.Fig. 3 Pseudo S wave in 2D TTIforwarding wavefield (a) Velocity vPz is 0; (b) Velocity vPz isn′t 0; (c) Velocity vPz isn′t 0 with filtering.

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

图 4

Fig. 4

图 4 三维空间中对称轴不同转角下的垂向波场(a) 转角:θ=0°, Φ=0°; (b) 转角:θ=45°, Φ=0°; (c) 转角:θ=45°,Φ=30°.Fig. 4 The vertical wavefield in 3D space with different symmetry axis angles(a) Angle: θ=0°, Φ=0°; (b) Angle: θ=45°, Φ=0°; (c) Angle: θ=45°, Φ=30°.

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

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

图 5

Fig. 5

图 5 TTI各向异性参数模型 (a) P波速度vPz； (b) δ； (c) ε； (d) θ.Fig. 5 TTI anisotropic parameter model (a) v<sub>Pz； (b) δ； (c) ε； (d) θ.

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

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

图 6

Fig. 6

图 6 (A) 单一震源的合成波场(p+r)和(B)辅助波场(r)在不同时刻的波场 从左至右时间依次是：0.2 s，0.3 s，0.4 s，0.6 s.Fig. 6 (A) Single source combine wavefield (p+r) at different time, (B) single source auxiliarywavefield (r) at different time From left to right time is 0.2 s, 0.3 s, 0.4 s, 0.6 s, respectively.

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

图 7

Fig. 7

图 7 单震源炮记录Fig. 7 The shot gather of a single source

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

图 8

Fig. 8

图 8 (a) 三维P波vPz速度模型； (b) TTI各向异性模型δ； (c) TTI各向异性模型εFig. 8 (a) 3D P wave velocity model vPz； (b) TTI anisotropic model δ； (c) TTI anisotropic model ε

图 9

Fig. 9

图 9 三维TTI模型的地表观测系统(a)及对应的地震事件的炮记录(b)Fig. 9 Surface geometry of 3D TTI model (a) and shot gather of the seismic events (b)

假定震源位置为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所示.

图 10

Fig. 10

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

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

图 11

Fig. 11

图 11 使用本文正演方式的逆时偏移(RTM)结果Fig. 11 RTM result by our forwarding method

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