1 引言

各向同性逆时偏移理论已趋于成熟，利用GPU加速技术可大幅提高计算效率，采用随机边界技术可避免波场存储的难题^{[1, 2]}.但在考虑各向异性逆时偏移算法时仍存在加速比不足和算子稳定性等困难^[3].其关键需要解决的科学问题是各向异性介质中动力学方程表达形式和算法实现.

本文研究弱各向异性介质中的地震波传播规律，可以有助于从中提取物理属性，在分析地层构造成因及岩性等工作中起到重要作用^[4].沉积地层长期的演化过程是由不同时期的薄层沉积来构成的，一种直观的各向异性假设是在薄层介质中沿层传播的地震波速与垂直于地层传播的波速不同而造成的各向异性称为TTI介质^[5].在研究TTI介质的运动学方程前，首先假设地层倾角为零时的情况，然后再将其拓展到地层的任意倾角情况.在水平层状介质假设的前提下的横向各向同性介质通常称为VTI介质^[6].

鉴于纵波勘探的广泛应用，纵波资料采集也比较成熟，研究纵波逆时偏移的各向异性方程是一个折中的选择，也是现阶段各向异性逆时偏移研究的主要方向^[4].很多地球物理学家已经在此方面做出大量的贡献:Alkhalifah(1998)在假设剪切横波的速度为零的情况下提出一种四阶的纵波各向异性方程，基本满足纵波在各向异性介质中的运动学规律，该方程的频散关系与弹性波方程的纵波解吻合^[7].实际的应用中发现，尽管设定了横波速度为零，但是仍然能够造成一种伪横波的出现，这是由于两个Thomson参数不相等而造成的后果，是方程近似导致的假象(Tsvankin，2001)^[8].在Alkhalifah的频散关系的基础上，周红波(2006)把四阶方程分解为两个二阶方程，有利于有限差分实现算法^[9].另外，还有一些学者得到不尽相同的纵波各向异性方程^[10].Duveneck(2008)从Hooke定律和运动学方程出发，可得到一个新的纵波各向异性的方程，但其数值计算非常的复杂^[11].ZhangH Z(2008)以TTI介质中的纵波相速度表达式为基础推导出一个对称的方程组，可利用高阶有限差分方法计算，但是该方法在VTI介质中可行，若应用于TTI介质中对称轴的角度发生剧烈变化时，仍存在稳定性问题^[12].随后，Fletcher(2009)提出一种加入横波分量的方法解决角度剧变时的稳定性问题，但是加入的横波会影响成像效果^[13].张宇(2009)从VTI介质中弹性波方程出发，通过自共轭旋转方式推导出一个TTI波动方程组^[14]，经证实该方程组在对称轴剧烈变化的时候仍然是稳定的解，但是在实际计算中会产生高频噪音.Duveneck(2011)提出直接从胡克定律中应力和应变的旋转，推导出TTI波动方程组，但实际计算过于复杂，在角度不变的情况下可退化为张宇2009年得到的方程^[15].可见，在TTI介质中，针对纵波研究一种稳定的不会带来横波噪音的传播算法是一个难题，也正是本文的主要研究内容.

2 弹性波方程的Thomsen表征

线性弹性动力学运动方程的一般表达为

(1)

式中:ρ 表示介质的密度､u表示质点的位移场，下标i､j､k､l分别可以表示笛卡尔坐标系中的x､y､z任意一个维度，cijkl表示四阶的刚度张量是Hooker定律中的应力与应变的关系张量.F表示单位质量单位体积介质所受的外力.

利用弹性系数刚度矩阵表示介质各向异性非常不直观.1986 年Thomsen 引入五个各向异性参数来描述弹性介质中的弱各向异性，每个参数都是弹性参数的函数，并且具有明确的物理意义^[16]:

(2)

其中V_pz 和V_sz 分别为P 波和S 波垂直TI介质各向同性面的相速度；ε 是P 波各向异性强度参数的度量，ε=0，纵波无各向异性，ε 越大，介质的纵波各向异性越强；δ 是纵横波速度比的相关参数；γ 是度量横波各向异性和横波分裂强度的参数，γ = 0 时横波无各向异性，γ 越大介质的横波各向异性越大.利用ε､γ 和δ 同样可以表征弱VTI介质中弹性矩阵元素^[17]:

(3)

其中

(4)

3 VTI介质的相速度

相速度是地震波在各向异性介质中传播的一个非常重要的物理量，它描述了波场的相位传播特性.可以通过求解Christoffel方程得到VTI介质中波前相速度公式(Tsvankin，2001；Thomsen，1986)^{[8, 16]}:

(5)

(6)

(7)

(8)

V_p(θ)是纵波的相速度，V_SV(θ)是SV 波的相速度，V_SH(θ)是SH 波的相速度，这些相速度都是传播角θ 的函数.

利用Thomsen 参数将VTI介质平面波P波相速度改写为^[17]

(9)

4 TTI介质中的P-SV 波动方程

相速度表达了实际的地震波速度在TTI介质中的变化规律，描述能量的传播过程.通过对相速度研究可以获得很多关于地震波传播相关的信息，下面介绍如何利用相速度来获得表达TTI介质中的波动方程表达式.

地震波传播角度θ､相速度Vp､频率ω 与波数k之间的有如下关系^[9]:

(10)

(11)

将(10､11)式带入相速度的计算公式(9)可以得到P-SV 波的频散关系为

(12)

将(12)式两边同时加入波函数p，再利用反傅立叶变换可以得到相应的时空域的波动方程表达式.但是由于上式是一个四阶的多项式，对应时空域是一个四阶的偏微分方程.引入辅助变量

(13)

可将(12)式转化为两个二阶的多项式，表达为

(14)

令η=1-f，其物理意义为横波速度与纵波速度之比.(14)式可化为

(15)

利用时空域与频率波数域的对应关系式

(16)

可以将(15)式变换，得到时空域偏微分方程组的表达式:

(17)

(17)式为VTI介质中的波动方程表达式，该表达式满足P-SV 波在VTI介质中的相速度传播规律.利用坐标旋转可以得到相应的TTI介质中的波动方程表达式为

(18)

其中

(19)

其中β 为TTI介质的对称轴与垂直方向的夹角，即VTI介质的对称轴旋转角度.

5 TTI介质的波场正演差分模拟

利用二阶中心差分格式近似二阶时间偏导数，从(18)式可以得到差分方程为

(20)

空间微分利用高阶空间差分近似可得

(21)

其中

(22)

利用N阶中心差分近似^[18-19]，有

(23)

式中分别用Δx､Δz表示差分网格的间距，al､bl分别为二阶偏导数和交叉导数的差分系数.

上述差分方程的稳定性条件，可以先假设方程中ε=δ，此时两个波场相等即p=q，那么求取其中一个方程的稳定性条件即可^[20]:

(24)

利用差分系数表示，并将延拓步长表示为

(25)

(26)

其中V_max 是P 波速度的最大值，h表示最大网格间距.

另一个物理稳定条件为P-SV 波相速度的平方不能小于零^[21]，即

(27)

当假设横波速度为零时，(27)式可化简为

(28)

6 实际测试 6.1 简单模型测试

设计256×256的网格区域，网格间距为10m，使用主频为25Hz的雷克子波，震源位于区域的中心，模型的纵波速度为3000 m/s，各向同性介质中的波场快照如图 1a所示，在ε 为0.3，δ 为0.1 的VTI介质中的波场快照示于图 1b:

可见，VTI介质的波场快照中间区域有菱形的波前面产生，这个是伪SV 波的波前面，其产生原因是由于方程本身是P-SV 波的运动学方程，简单假设V_sz=0并不能使SV 波的相速度为零，其相速度的表达式可表示为^[22]

(29)

P 波逆时偏移中需要将SV 波消除.当ε =δ时，VSV(θ)≡0此时为各向同性或者椭圆VTI介质，在此种介质中该菱形波前面并不存在，如图 1d所示.另一方面，方程中并未考虑震源的各向异性，把震源激发的位置设置为椭圆VTI介质或各向同性介质即可消除该菱形SV 波(如图 1c所示)令震源激发点处有ε=δ，其他区域都保持ε=0.3，δ=0.1.在海洋勘探中，震源是在水中激发，是各向同性介质，所以可以使用此方法消除伪SV 波.利用此方法计算的P 波点源响应显示于图 2.若Vsz ≠0时，将会出现SV 波，实际为弹性波的正演模拟(如图 3所示).其结果与图 2中纯P 波正演所得结果对比，可验证P波各向异性介质的正演模拟的正确性.但是，在陆地勘探中此方法适用条件不能得到满足，所以仍需寻找新的解决办法.后文以复杂模型的计算结果来说明本文解决此问题的方法.

6.2 BP2007TTI模型测试

BP2007TTI模型是典型的以各向异性介质为基础的模型，该模型的数据是利用弹性波正演获得，其P波速度模型和各向异性参数示于图 4中，模型的网格大小为6.25 m×6.25 m，垂直方向包含1801个采样点最深达到11.25km，水平方向包括12536个采样点，横向测线长度为78.35km.利用BP模型中的速度和各向异性参数做单炮的正演模拟，并取得200ms时刻的波场快照(如图 5所示).

如图 5a所示的波场快照，当假设V_sz=0 时在延拓过程中会产生不稳定性问题，去掉Vsz=0的假设可以使方程稳定，但是在波场中会加入SV 波的噪音(如图 5c所示)，控制SV 波的速度可以使SV波的各向异性和反射系数最小化^[23]，令

(30)

其中C为小于1 的常数，本文实验了当C=0.3 和C=0.8时的波场快照(如图 5b､d所示).可见当C=0.8时SV波已经完全消失，与图 5a中的P波的相位和运动学特征相同，此时的波场是一个稳定的纯P波场.

BP2007模型左端盐丘的各向异性很小，用各向同性的逆时偏移也可得到较好的成像效果，在中部的盐丘(图 7a)和右测的断层结构(图 6a)是成像的难点，利用各向同性的偏移算法很难准确成像(图 7b､图 6b).利用本文的VTI介质的逆时偏移方法和TTI介质的逆时偏移方法对这两个部分进行偏移成像，其结果分别示于图 7c､图 7d､图 6c和图 6d.从图中可见TTI逆时偏移的成像效果要优于VTI逆时偏移和各向同性逆时偏移的结果，归位更准确，聚焦能量更强.

7 结论

本文针对以往TTI介质中的纵波模拟的数值稳定性问题，指出在角度变化剧烈时的不稳定性问题产生的原因是:方程的推导与计算过程中人为假设横波速度为零.本文指出引入横波分量可以使数值计算稳定，并以TTI介质下的P-SV 相速度规律和波动方程为基础，利用中心差分格式实现方程的有限差分模拟，并分析其数值稳定性条件.实际测算表明，P-SV 波的方程中包含横波分量，若假设SV的速度为零，对称轴角度剧烈变化时会导致方程的差分格式不稳定；若加入SV 波，选择合适的SV 波速度可以使SV 波的全区各向异性和反射系数达到极小，并可有效的抑制SV 波对纵波勘探的影响.因此，本文的方法是一种稳定的TTI介质中的纵波逆时偏移方法.模型数据计算表明，TTI逆时偏移成像结果归位准确，成像信噪比高.