地下介质普遍存在各向异性特征, 在各向异性较弱的地区, 采用各向同性假设能够得到较好的成像结果, 但是如果各向异性较强, 各向同性成像结果在深度和横向上都会产生较大误差。因此, 研究各向异性参数建模和偏移成像是非常必要的, 而各向异性射线追踪作为建模和成像的基础, 决定了建模精度和成像效果, 因此, 提高各向异性射线追踪的精度和稳定性是研究的重点。

各向异性介质非常复杂, Thomsen^[1]在1986年提出了弱各向异性（VTI）介质模型, 建立了Thomsen参数。在这个模型基础上, Alkhalifah^[2]提出了近似的qP波相速度方程和频散关系式。Č ervený ^[3]解决了各向异性介质中射线追踪问题, 极大的推动了射线层析反演的发展。Zheng和Psencik^[4]获得了更为普遍的关于qP波和qS波的解析表达式。Psencik^[5]等利用扰动理论给出了在光滑非均匀各向异性介质中qP波的近似射线方程。Farra^[6]利用扰动理论导出了qP波的一阶射线追踪方程和傍轴射线追踪方程。上述射线追踪方法的声波方程都是在Alkhalifah^[2]提出的令横波等于零的假设下得到的, 该假设条件在一些特殊条件下精度较低、稳定性较差, 因此不能适用于所有地质情况, 在实用性上存在一定的局限性。

本文在前人研究的基础上, 借鉴Alkhalifah^[2]的思路, 从VTI介质准确的相速度方程出发, 将根号项用Taylor级数展开, 得到P-SV波分离后的相速度, 以及VTI介质的频散关系。基于这个频散关系式, 采用傅里叶反变换得到qP波方程, 带入平面波方程后通过高频近似得到相应的程函方程和射线方程。最后在二维和三维不同介质中对Alkhalifah^[2]射线追踪和本文P-SV波分离射线追踪两种方法的旅行时等值线进行了绘制, 并且对其稳定性和精度进行了对比分析。

Alkhalifah^[7~9]认为, 即使是在很强的各向异性介质中, P波速度和旅行时是独立于横波速度V_SO的, 因此P波旅行时仅是纵波速度V_PO与各向异性参数ε 、δ 的函数。从弹性波方程的准确频散关系出发, 将V_SO=0代入其中, 得到VTI介质的频散关系：

式中：NMO是正常时差校正, V_NMO是NMO速度, V_PO为纵波速度, , P是射线参数, , τ 是时间参数, η 是非椭圆率参数, 。

根据频散关系进一步导出相应的程函方程：

根据Aki^[10]等提出的运动学方程解满足的关系可以得到射线方程：

式中：d表示微分求导, T是旅行时, ds是空间步长, , a_x= , 表示参数对x, y, z的偏导数。

上述方法由于直接令横波等于零（V_SO）近似得到VTI介质的频散关系, 当δ 为负值, 旅行时计算会出现误差; 另外, 当η < 0时, 射线追踪会出现不稳定现象^[6]。针对这些缺点, 本文基于Pestana等^[11]给出的频散关系, 进而推导得到程函方程和射线方程。

Pestana等推导出新的频散关系如下^[11]：

式中：ω 是频率, V_h表示纵波横向速度, =(1+2ε ) , V_n表示纵波垂向速度, , k是波数, 。

上式两边同乘以傅里叶域波场φ (k_x, k_y, k_z, ω ), 分别对k_x, k_y, k_z, ω 做反傅里叶变换可得三维VTI介质的qP波方程：

将平面波解带入上式qP波方程, 基于高频近似, 取实部中ω ⁴项前的系数对应相等, 得到相应的程函方程^[12~14]：

根据Aki 和Richards（1980）运动学方程的解满足如下关系^[8]：

其中, φ 是程函方程, 利用上式对程函方程进行推导即可得到新的射线方程：

（1）稳定性对比

Alkhalifah^[2]认为令横波等于零得到的qP波方程解的振幅项是正负指数形式, 即当η < 0时, 四组解中至少有一个解的振幅会随时间呈指数快速增大, 这会严重影响数值计算的稳定性, 因此, 当η < 0时, 该方法不能用于波场模拟。而新方法的射线方程则弥补了这一缺陷, 在包括η < 0的任何情况下都不会出现极度不稳定现象, 能够更好的用于波场模拟。

为了验证这两种射线方程的稳定性, 针对η < 0, 选取各向异性参数ε =0.1, δ =-0.51进行数值实验。图1a为常速模型中利用新射线方程求解得到的旅行时等值线, 震源在模型的中心点位置; 图1b为在常梯度速度模型中利用新射线方程求解得到的旅行时等值线, 震源在模型底部的中心位置。

通过数值实验可知, 当η < 0时, Alkhalifah^[2]的射线方程数值计算不稳定, 得到的射线点坐标和相应旅行时为无限大（这种情况无法用绘图表示）, 与Alkhalifah给出的结论一致。而从图1可以看出, 新的射线方程在η < 0时仍可以进行准确的旅行时计算, 且该方程不存在稳定性条件的限制, 对任意η 均适用。可见, 新的射线追踪算法稳定性更好、适用性更强。

（2）精度对比

程函方程和射线方程的精度直接受频散关系的控制, 而频散关系又由相速度公式得到, 因此可以用相速度检验两种方法的精度。以下将两种近似方法得到的相速度跟真实相速度作对比来验证方法的精度。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 新射线方程旅行时等值线图Fig.1 Traveltime contours using new ray equation

设计3种模型, 具有相同的纵、横波速度, V_PO=2000m/s, V_SO=1000m/s, 相角范围为0° ~ 90° , 各向异性参数有所区别, 分别为ε =δ =0、ε =0.065, δ =0.059和ε =0.11, δ =-0.035。最终计算得到的相速度曲线如图2~4所示。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 各向同性介质相速度曲线(ε =δ =0)Fig.2 The curves of phase velocity for isotropic medium

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 ε =0.065, δ =0.059时的VTI介质相速度曲线Fig.3 The curve of phase velocity for Mesa clay Shale VTI medium

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 ε =0.11, δ =-0.035时的VTI介质相速度曲线Fig.4 The curve of phase velocity for Taylor Sand VTI medium

由图2可知, 当介质为各向同性时, Alkhalifah方法得到的相速度和新方法得到的相速度与准确相速度曲线完全相同。由图3可知, 当各向异性较弱时, Alkhalifah方法得到的相速度和新方法得到的相速度几乎没有差别, 二者与准确相速度的误差随相角的增大而变大。由图4可知, 当各向异性逐渐增大, 且δ 变为负值时, Alkhalifah方法得到的相速度和新方法得到的相速度从20° 左右开始出现了偏差, 虽然偏差较小, 但是新方法的相速度更加接近于真实相速度, 所以当δ < 0时, 新方法的精度更高。

综上所述, 新方法采用P-SV波分离手段近似得到的qP波相速度比直接令V_S0=0得到的qP波相速度更稳定、更精确, 能够更好的模拟波场传播, 有利于提高VTI介质参数建模的精度和偏移成像的效果。

对常梯度速度模型（2000~4000m/s, 梯度为1m/s）进行射线追踪, 源点坐标为(1000m, 0m), 射线扫描范围为-90° ~90° , 角度间隔为0.1° , 各向异性参数分别为ε =0.065, δ =-0.059和ε =0.11, δ =-0.035。射线追踪结果分别如图5和图6所示。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 VTI介质ε =0.065, δ =-0.059常梯度速度模型等值线图Fig.5 Analysis on VTI medium velocity isoline in constant gradient model and composite

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 VTI介质ε =0.11, δ =-0.035常梯度速度模型等值线图Fig.6 Analysis on VTI medium velocity isoline in constant gradient model and composite

从图5可以看出, 当δ > 0时, 2种方法的旅行时波前面重合, 说明2种方法得到的旅行时一致。从图6可以看出, 当δ < 0时, 2种方法的旅行时波前面不再重合, 略有偏差, 说明2种方法计算得到的旅行时不一致, 从上文分析可以知道, 本文方法得到的旅行时更加准确。

各向异性介质射线追踪是各向异性参数建模和偏移成像的基础, 它的精度和稳定性直接影响建模和成像的效果, 是各向异性研究中的关键问题。本文基于前人提出的一种新的VTI介质频散关系, 导出了相应的程函方程和射线方程。经过分析, 该方程能够弥补传统方法在η < 0时出现的不稳定现象, 适用于各种参数的VTI介质, 并且对于δ < 0时的旅行时计算更加精确。模型试算表明了该方法的优越性。