室内研究与野外观测表明大多数岩石和地层存在明显的各向异性.引起地层各向异性的主要因素有：岩石的结构各向异性，地层方向应力引起的各向异性，岩性各向异性，岩石晶体定向排列引起的各向异性，以及岩石裂隙裂缝引起的各向异性(Crampin，1984; Helbig，1981; 李芳等，2012).地震波在各向同性介质与各向异性介质中的传播特征不同，在各向异性介质有三种体波：qP波，qSV波和qSH波，这三种波各自以自身的传播速度与相互垂直的极化方向传播(Zhou and Greenhalgh，2006).因此，如果地层具有各向异性，那么地震走时层析成像方法必须适用于这种介质模型，否则，地震走时反演结果可能会产生错误的地质解释.

经过40多年的发展历程，地震各向同性介质走时层析成像方法已经比较成熟，很多学者利用各种射线追踪算法、数学反演算法和走时扰动方程得到了不同的地震走时层析成像方法.这里列举一些在勘探地球物理学方面有代表性的地震走时层析成像方法.Zelt等(2006)的有限差分初至波走时层析成像方法得到了广泛的应用，利用这种初至波走时层析成像方法在地下水污染的圈定和调查方面取得了良好的应用效果.Tryggvason等(2009)结合初至波走时层析和静校正方法，将其运用于核废料处置的勘探选址研究，所得地震剖面的可靠性与分辨率均得到了较大提高.Zhou(2006)提出了地震多尺度可变形地层走时层析成像方法，数值测试结果表明该方法比基于晶格或者网格节点的走时层析成像方法具有较高的成像精度.刘玉柱和杨积忠(2014)利用先验信息在初至波走时层析成像过程中加入正则化 约束条件，反演结果使原来隐藏的特征变得更加清晰.

关于地震各向异性介质走时层析成像的研究也在不断升温，周辉和何樵登(1995)利用层状TI介质的P波和SH波走时，结合双曲线拟合法和遗传算法反演得到了五个弹性模量参数值.张文生和何樵登(1999)利用跨孔qP波和qSV波进行了统计法走时层析成像，提高了计算结果的稳定性和精度.卢明辉等(2005)利用正交各向异性介质P波走时反演得到了单层介质的Thomsen参数.刘玉柱等(2014)利用多参数联合走时层析成像方法对VTI介质进行走时反演取得了较好的数值测试结果.国外对地震各向异性介质走时层析成像的研究相对较多，Cerveny(2001)是较早从事各向异性走时层析成像的学者，他运用各向异性介质“程函方程”的哈密尔顿形式，推导出了各向异性介质走时扰动方程.后续学者利用该走时扰动方程和不同的观测系统对各种各向异性介质模型进行反演取得了不同的效果.Jech(1988)对该走时扰动方程进行了修改，并且实 现了三维各向异性介质走时层析成像方法.Chapman和Pratt(1992)、 Pratt和Chapman(1992)提出了各向异性介质走时反演的线性化方法，利用井间观测方式进行了走时反演.Wang和Tsvankin(2013)提出了网格化的各向异性介质走时层析成像方法，该方法极大地提高了成像的精度.Wang(2014)对各向异性介质弯曲射线法产生的高度非线性问题，提出了修改的牛顿迭代法用于稳定非线性问题的解，并利用井间观测方式反演了各向异性介质模型参数.纵观大多数各向异性介质走时层析成像方法，它们都是采用了弱各向异性介质假设条件.实际上，岩石和地层的各向异性强度是不确定的，将它们假设为弱各向异性介质具有不合理性.Zhou等(2008)结合各向异性介质初至波射线追踪算法、走时扰动方程和非线性共轭梯度，实现了地震各向异性初至波走时层析成像方法.这种地震走时层析成像方法的优点是能够适用于任意强度的各向异性介质模型，它的各向异性介质走时扰动方程避免了qSV波和qSH波的奇异值问题，非线性共轭梯度算法可以同时反演出几个弹性模量参数(或者Thomsen参数).

国内外关于地震各向异性介质走时层析成像的文章不是很多，其技术远远没有达到成熟的阶段，因此仍然有很多的理论和数值模拟研究值得探索.本文以Zhou等(2008)的各向异性初至波走时层析成像方法为基础，利用反射波射线追踪算法实现了TI介质qP反射波非线性走时层析成像.由于地表反射波观测方式相比井间和垂直地震剖面观测方式的成本更低，因此TI介质qP反射波非线性走时层析成像方法具有更加广泛的实用价值.在数值模拟部分，利用qP反射波走时重建层状TI介质模型和块状异常体TI介质模型的图像，用于测试各向异性介质参数重建和成像的可能性.通过分析各种模型的弹性模量参数与Thomsen参数反演结果，验证TI介质qP反射波走时层析成像方法的正确性.

Cerveny推导出了各向异性介质走时扰动方程的表达式.很多学者对该公式做了不同的修改，Zhou推导出了各向异性介质走时扰动方程的三种表达形式，这三种表达式被称为一阶走时扰动方程的波速偏微分形式，其中最为简单的是走时扰动方程的相速度偏微分形式(Zhou et al.，2008)

式中，v_p和v_g 分别代表各向异性介质的相速度和群速度，如果是二维介质情形，可以利用六个物理量 m_v∈ c₁₁，c₁₃，c₃₃，c₄₄，c₆₆，θ₀(五个弹性模量参数和地层对称轴倾角)来表述各向异性介质模型.∂v_p/∂m_v 是相速度对各向异性参数的偏微分表达式，它被称为速度敏感函数，它是各向异性介质走时 扰动方程的重要组成成分.Zhou对速度敏感函数计 算方法做过详细介绍(Zhou and Greenhalgh，2005). 由于相速度和群速度均是波前法向向量 n和模型参数向量 m 的函数，所以(1)式可以写成如下形式:

各向异性模型参数化后，可以得到相速度 v_p(m，n)与群速度 v_g(m，n)模型.

利用(2)式进行各向异性介质走时层析成像时，将研究区域 Ω 划分成 N_Ω 个小单元 ，这样公式(2)可以写成分段求和的形式(Zhou et al.，2008)

式中，R_k 是射线路径 R 在第 k 个单元∂Ω_k 内的射线长度，x_A(k)和x_B(k) 是射线段 R_k 两个端点的坐标值.根据网格四个端点的数值，可以利用拉格朗日插值公式计算 v_p，v_g和∂v_p/∂m_v 在射线段 R_k 两个端点的值(Zhou et al.，2008; Kreyszig，1993).

当震源与检波器的数量比较大时，它们之间存在大量的射线路径，公式(3)可以写成非线性反演的矩阵形式(Zhou et al.，2008)

式中上标-g表示广义逆矩阵，W_d和W_m分别代表数据与模型参数的权重矩阵，λ 是平衡数据拟合度与模型粗糙度之间的变量.(4)式可以利用非线性共轭梯度算法进行求解(Zhou et al.，1992)，其他变量的表达式如下:

式中，τ_j^ob和τ_j^syn(m _l－1)分别为第 j 条射线的观测走时和理论计算走时； Z_m(m _l)是由每个各向异性参数 m_v 的上限 a_vi和下限 b_vi 定义的约束算子； J _l 为雅可比矩阵，它的计算表达式如下：

下面利用各向异性介质qP反射波走时层析成像进行数值模拟.数值模拟包括两部分：第一部分是利用qP反射波走时反演TI介质模型的弹性模量参数；第二部分是利用qP反射波走时反演TI介质模型的Thomsen参数.测试模型包括层状介质模型和块状异常体模型.Zhou对各向异性介质初至波走时层析成像做过相应的研究工作，并且指出当地层对称轴倾角不固定时，算法很难反演得到合理的倾角参数分布.原因是相速度和正演走时均对地层对称轴倾角非常敏感，因此本文假设模型的地层对称轴倾角为已知，只反演模型的弹性模量参数或者Thomsen参数分布.数值模拟实验的观测走时是射线追踪得到的正演走时，并且添加了一定百分比的高斯噪声.在非线性走时反演时，W_d和W_m是根据数据和模型先验信息确定的权重矩阵，用以减小反演的多解性.文中假设先验信息的权重相同，即数据 的权重矩阵和模型的权重矩阵均为单位矩阵(W_d= I_d和W m= I m). 阻尼因子λ是通过多次数值测试，比较反演的收敛情况而选取的数值大小.在反演迭代过程中，利用了阿尔法中值滤波器剔除每次迭代结果的高波数噪声.地震走时层析成像的反演算法是非线性共轭梯度，Zhou等(1992)将其与不同反演算法做过数值和实际数据试验，证明了它是一种有效而快速的反演算法.

下面分别对层状介质模型和块状异常体模型进行TI介质反射波非线性走时反演，求取模型的弹性模量参数.

(1)层状介质模型的弹性模量参数反演.层状介质模型的长度为800 m，深度为600 m，第一层介质的弹性模量参数为c₁₁=8.81 GPa，c₁₃=2.78 GPa，c₃₃=7.34 GPa，c₄₄=2.28 GPa，c₆₆=3.88 GPa，对称轴倾角为θ₀=0°.第二层介质的弹性模量参数为 c₁₁=21.49 GPa，c₁₃=10.85 GPa，c₃₃=19.54 GPa，c₄₄=8.71 GPa，c₆₆=8.45 GPa，对称轴倾角为θ₀=90°. 第三层介质的弹性模量参数为 c₁₁=13.86 GPa，c₁₃=6.23 GPa，c₃₃=12.36 GPa，c₄₄=4.35 GPa，c₆₆=5.09 GPa，对称轴倾角为θ₀=0°.采用地表激发和接收的观测方式，炮点与检波点位置相同，道间距为25 m，炮点与检波点均为33个，模型的网格间距为50 m.将第一层介质参数作为初始模型用以反演第二层和第三层介质的弹性模量参数.为了改善反演成像的效果，我们对模型设置了三个反射界面，充分利用三个反射界面的反射波走时信息.图 1a为层状弹性模量参数模型的反射波射线路径分布图，射线覆盖了整个模型范围，图 1b是利用相同的各向异性介质正演走时结合纵波速度走时层析成像所得速度反演结果，可以看出速度反演结果效果不理想，第一层速度值变化较大，第二层下界面位置不正确.图 2为层状介质模型的弹性模量参数(a—c—e—g)与反演结果(b—d—f—h)的对比图，左侧为真实弹性模量参数，右侧为反演结果，从整体上分析，c₁₁，c₁₃，c₃₃和c₄₄均得到了较好的反演效果.由迭代反演收敛曲线(图 9)可知反演得到了较好的收敛水平，表明反演结果与真实模型相接近.

图1

Fig. 1

图 1 层状介质弹性参数模型的射线分布与速度反演结果 Fig. 1 The raypath distribution and inverted velocity profile of the elastic moduli layered model

图2

Fig. 2

图 2 层状介质模型的弹性模量参数与反演结果 Fig. 2 The true elastic moduli parameters and inverted results for the elastic moduli layered model

(2)块状异常体模型的弹性模量参数反演.块状异常体模型的长度为800 m，深度为500 m，背景介 质的弹性模量参数为c₁₁=12.0 GPa，c₁₃=5.0 GPa，c₃₃=10.0 GPa，c₄₄=4.0 GPa，c₆₆=6.0 GPa，对称 轴倾角为θ₀=0°.块状异常体的弹性模量参数为 c₁₁=20.0 GPa，c₁₃=10.0 GPa，c₃₃=15.0 GPa，c₄₄=8.0 GPa，c₆₆=12.0 GPa，对称轴倾角为θ₀=45°.采用地表激发和接收的观测方式，炮点与检波点位置重合，水平道间隔为16 m，炮点与检波点均为51个，模型的网格间距为50 m.将背景介质参数作为初始模型用以反演块状异常体的弹性模量参数.为了改善反演成像的效果，我们对模型设置了两个反射界面，一个是块状异常体的上边界，另一个是模型的底界面.图 3a为块状弹性模量参数模型的反射波射线路径分布图，射线覆盖了块状异常体区域，图 3b是利用相同的各向异性介质正演走时结合纵波速度走时层析成像所得速度反演结果，反演所得异常体的形状能够得到识别，但是异常体的数值变化非常大，整体效果不太理想.图 4为块状异常体模型的弹性模量参数(a—c—e—g)与反演结果(b—d—f—h)的对比图，左侧为真实弹性模量参数，右侧为反演结果，c₁₃，c₃₃，c₄₄均得到了较好的反演效果，由于c₁₁的敏感性较弱，它的反演效果相对较差.如果某个参数的速度偏导数相对于其他参数的速度偏导数不灵敏，那么利用观测走时数据来重建该参数剖面是不太现实的.国内外文献均是针对敏感性较强的各向异性参数进行反演，而很难把灵敏度弱的各向异性参数反演出来.由迭代反演收敛曲线(图 9)可知反演的走时数据拟合较好，说明反演结果与真实模型相接近.

图3

Fig. 3

图 3 块状异常体弹性参数模型的射线分布与速度反演结果 Fig. 3 The raypath distribution and inverted velocity profile of the elastic moduli blocky abnormal body model

图4

Fig. 4

图 4 块状异常体模型的弹性模量参数与反演结果 Fig. 4 The raypath distribution and inverted velocity profile of the Thomsen layered model

下面分别对层状介质模型和块状异常体模型进行TI介质反射波非线性走时反演，求取模型的Thomsen参数.

(1)层状介质模型的Thomsen参数反演.层状介质模型的长度为800 m，深度为600 m，第一层介质的Thomsen参数为α₀=2.71 km·s^-1，β₀=1.51 km·s^-1，ε=0.10，δ^*=－0.13，γ=0.35，对称轴倾角为θ₀=0°；第二层介质的Thomsen参数为α₀=4.72 km·s^-1，β₀=2.89 km·s^-1，ε=0.26，δ^*=0.17，γ=0.17，对称轴倾角为θ₀=90°；第三层介质的Thomsen参数为α₀=3.37 km·s^-1，β₀=1.83 km·s^-1，ε=0.18，δ^*=0.02，γ=0.26，对称轴倾角为θ₀=0°.采用地表激发和接收的观测方式，炮点与检波点位置相同，水平道间距为25 m，炮点与检波点均为33个，模型的网格间距为50 m.将第一层参数作为初始模型用以反演第二层和第三层的Thomsen参数.为了改善反演成像的效果，我们对模型设置了三个反射界面，充分利用了三个反射界面的反射波走时信息.图 5a为层状Thomsen参数模型的反射波射线路径分布图，射线覆盖了整个模型范围，图 5b是利用相同的各向异性介质正演走时结合纵波速度走时层析成像所得速度反演结果，反演所得结果的中间层速度值变化比较大，并且中间层上下两个边界位置不易识别.图 6为层状介质模型的Thomsen参数(a—c—e—g)与反演结果(b—d—f—h)的对比图，左侧为真实Thomsen参数，右侧为反演结果，α₀，ε和δ^* 均得到了较好的反演效果，因为 β₀ 的敏感性较弱，它的反演效果稍差，其界面位置很难分辨清楚.由迭代反演收敛曲线(图 9)可知观测走时与正演走时之间的残差较小，反演得到的结果与真实模型相接近.

图5

Fig. 5

图 5 层状介质Thomsen参数模型的射线分布与速度反演结果 Fig. 5 The raypath distribution and inverted velocity profile of the Thomsen layered model

图6

Fig. 6

图 6 层状介质模型的Thomsen参数与反演结果 Fig. 6 The true elastic moduli parameters and inverted results for the Thomsen layered model

(2)块状异常体模型的Thomsen参数反演.块状异常体模型的长度为800 m，深度为500 m，背景介质的Thomsen参数为α₀=3.16 km·s^-1，β₀=2.0 km·s^-1，ε=0.25，δ^*=0.3，γ=0.25，对称轴倾角为θ₀=0°；块状异常体的Thomsen参数为α₀=3.87 km·s^-1，β₀=2.83 km·s^-1，ε=0.18，δ^*=1.14，γ=0.25，对称轴倾角为θ₀=45°.采用地表激发和接收的观测方式，炮点与检波点位置重合，水平道间隔为16 m，炮点与检波点均为51个，模型的网格间距为50 m.将背景介质作为初始模型用以反演块状异常体的Thomsen参数.为了改善反演成像的效果，我们对模型设置了两个反射界面，一个是块状异常体的上边界，另一个是模型的底界面.图 7a为块状异常体Thomsen参数模型的反射波射线路径 分布图，射线覆盖了块状异常体区域，图 7b是利用相同的各向异性介质正演走时结合纵波速度走时层析成像所得速度反演结果，可以看出速度反演结果效果非常不理想，完全无法识别异常体的形态和位置.图 8为块状异常体模型的Thomsen参数(a—c—e)与反演结果(b—d—f)的对比图，左侧为真实Thomsen参数，右侧为反演结果，参数 α₀，β₀和δ^* 均得到了较好的反演效果，由于块状异常体的 ε 值比背景介质的 ε 值更低，反演算法比较难反演出 ε 参数.由图 9的收敛曲线可知该模型的走时反演得到了较好的收敛效果，表明反演模型与真实模型的相似性较好.

图7

Fig. 7

图 7 块状异常体Thomsen参数模型的射线分布与速度反演结果 Fig. 7 The raypath distribution and inverted velocity profile of the Thomsen blocky abnormal body model

图8

Fig. 8

图 8 块状异常体模型的Thomsen参数与反演结果 Fig. 8 The true elastic moduli parameters and inverted results for the Thomsen blocky abnormal body model

图9

Fig. 9

图 9 四种模型的迭代反演收敛曲线 Fig. 9 The iterative inversion convergence curve for these four models

各向异性现象在地球介质中普遍存在，地震走时层析成像需要考虑地层的各向异性.本文的各向异性TI介质qP反射波走时层析成像方法适用于任意强度的各向异性介质，该方法相比井间和垂直地震剖面走时层析成像方法具有更加广泛的应用价值.由于qP波相速度对弹性模量参数和Thomsen参数的速度偏微分不同，所以利用qP反射波可以分别反演出弹性模量参数和Thomsen参数.在数值模拟部分，利用qP反射波分别反演层状介质模型和块状异常体模型，均得到了较好的反演结果.这些模型的迭代反演收敛曲线均得到了较好的收敛效果，表明反演结果与真实模型之间相接近.与各向同性介质走时反演结果相比较，各向异性介质走时反演结果具有较好的识别能力.数值模拟结果证明了各向异性TI介质qP反射波走时层析成像算法的正确性.由于受地表观测系统的限制，反射波走时层析成像相比井间直达波走时层析成像的反演效果较差，在固定检波器数量的情况下，应最大限度地增加qP反射波射线覆盖范围以减小走时反演的多解性和非线性.根据反演理论，已知的信息量越多反演的不确定性越弱，如果能够综合利用qP波、qSV波和qSH波进行多波走时联合反演，将有助于改善各向异性参数的反演效果，因此各向异性TI介质多波走时联合反演值得更深入的研究.