地下介质普遍存在各向异性特征。在地层接近水平或各向异性程度较弱的地区, 各向同性或VTI介质假设条件能够得到较好的地下介质成像效果。但是, 如果各向异性程度较强, 并且地层倾角较大, 基于各向同性和VTI介质的假设则无法得到地下介质的高精度成像结果。此种情况下通常会出现成像道集无法拉平、构造模糊或失真、成像位置和深度存在误差等现象, 给后续油藏描述和储层预测带来困难, 从而降低钻井成功率^[1-4]。因此, 随着地震勘探的不断发展以及勘探目标越来越复杂, 各向同性和VTI介质的假设条件已经不能够满足勘探精度的需求, 需要研究更加符合地下地质特征的TTI介质假设条件下的建模成像技术以提高复杂构造的勘探精度, 而TTI射线追踪是TTI建模成像的关键技术, 它决定了建模成像的精度和效率, 是复杂构造精确成像研究的重要内容之一。

目前各向异性介质射线追踪方法主要有^[5]:打靶法、波前构建法、插值法、最短路径法和程函方程法等。打靶法基于渐进射线理论^[6], 适用于一般各向异性介质, 但对于TTI这种复杂介质则会出现阴影区问题。波前构建法是根据当前的波前利用射线理论估计新的波前, 白海军等^[7]虽然利用群速度和相速度表达的射线追踪方程避免了复杂的特征值求解问题, 但是实现方案较为复杂。插值法基于费马原理^[8-10], 分为向前确定走时和向后确定射线路径两个过程, 在各向异性介质中往往得不到精确的出射点位置。最短路径法基于惠更斯原理和费马原理^[11-13], 能同时求取震源到模型所有节点的走时和相应的射线路径, 但只适用于各向异性介质初至波射线追踪。VIDALE^[14]基于扩张波前的思想建立了经典的程函方程射线追踪系统, 之后, ALKHALIFAH^[15]提出了各向异性程函方程, 并在此基础上建立了各向异性射线追踪系统, 该方法稳定性强、精度高, 是VTI介质射线追踪最主要的方法, 但是在TTI介质中, 程函方程则非常繁琐, 无论是对于推导射线追踪方程还是对于射线追踪效率, 影响都很大, 因此该方法不适用于TTI介质。

本文在前人研究的基础上, 利用Hamilton相速度公式推导出相速度射线追踪方程, 再将TTI介质群速度和相速度公式代入相速度射线追踪方程得到TTI介质射线追踪方程, 然后利用射线参数与相角的关系将射线追踪方程中射线参数与时间的偏微分公式转换为射线出射角与时间的偏微分公式, 最后采用龙哥库塔法求解射线追踪系统的微分方程组, 从而得到射线追踪每一步延拓的坐标、出射角和时间。这种方法相比较于程函方程法, 公式简单、物理意义明确, 运算效率高, 并且解微分方程得到的是出射角而不是程函方程方法中的射线参数, 可以提供后期层析反演所需要的数据, 更有利于运算。

图 1为TTI介质相角与群角示意图, 其中, TTI介质中出射角(射线方向与z轴的夹角)为θ, 对称轴倾角(对称轴与z轴的夹角)为θ_c, 在本文公式中记为θ′, 射线方位角(射线方向水平投影与0方位的夹角)为φ, 对称轴方位角(对称轴水平投影与0方位的夹角)为φ_c, 在本文公式中记为φ′, 则TTI介质中相角γ可表示为:

图 1(Fig. 1)

图 1 TTI介质相角与群角示意

(1)

通过相角可以得到TTI介质相速度公式:

(2)

其中:

式中:v为地震波相速度; v_P0是沿对称轴方向的纵波速度; ε和δ是Thomsen各向异性参数。

BERRYMAN提出射线速度可以表示为波动矢量和相速度的函数, 并给出了射线速度公式^[16]:

(3)

式中:v_G表示地震波射线速度, 即群速度; v是地震波相速度; 波动矢量k=k_xI+k_yJ+k_zK; I, J, K代表直角坐标系三个方向的单位向量。

吴国忱等^[17]以公式(3) 为基础, 推导出的TTI介质三维群速度公式为:

(4)

基于相速度的射线追踪方法是从Hamilton相速度公式出发推导运动学射线追踪系统, 然后将TTI介质相速度和群速度代入该系统, 得到TTI介质射线追踪方程, 此时的方程是坐标、射线参数和旅行时的微分方程, 解此方程就会得到当前点的坐标、射线参数及传播时间, 可以用作计算下一个传播点的初始值。但是射线参数只是中间变量, 与速度和时间并无联系, 无法直接作为初始值使用, 因此利用射线参数与角度的关系进一步将TTI介质射线追踪方程变为坐标、出射角和旅行时的微分方程, 解方程得到的坐标、出射角及传播时间可以直接作为初始值进行计算。通过该微分方程就可以从出射点一直追踪到检波点, 计算出一条射线曲线及旅行时信息。基于相速度的射线追踪方法具体推导过程如下。

首先从Hamilton相速度公式出发^[18]:

(5)

得到相应的运动学射线追踪方程:

(6)

式中:x_i是坐标; T是旅行时; τ是时间步长; K是Hamilton相速度; p是射线参数; v是TTI介质相速度; v_G是TTI介质群速度。

将TTI介质相速度((2) 式)和群速度((4) 式)公式代入射线追踪系统((6) 式), 可以得到基于相速度的三维TTI介质射线追踪方程:

(7)

将射线参数与出射角的关系式p_x=sinθcosφ/v, p_y=sinθsinφ/v, p_z=cosθ/v代入公式(7), 即可得到三维TTI介质射线追踪方程出射角形式:

(8)

其中:

公式(8) 即为基于相速度的TTI介质三维射线追踪方程, 它是以等时间为步长进行延拓的, 即每一点旅行时一样。先计算每一步延拓点的坐标, 然后采用相速度对坐标的偏导数计算出射角。基于相速度的TTI介质三维射线追踪方程形式简单, 方程中的变量物理意义明确, 并且只计算5个参数, 运算稳定、计算效率高。

最后用龙哥库塔法求解微分方程组(公式(8)), 可以得到射线追踪每一步延拓的坐标、出射角和时间。龙哥库塔公式为:

(9)

式中:Y_n和Y_n+1分别是n次和n+1次延拓得到的值; h是延拓步长; K为龙哥库塔四阶计算值; t_n是n次延拓的时间。

在验证方法的精度和效率之前, 首先验证方法的正确性。将TTI声波正演模拟产生的单炮记录与TTI射线追踪得到的旅行时曲线进行叠合, 如果完全重合, 说明本文方法所得到的旅行时信息正确; 如果不完全重合, 说明本文方法存在问题, 需要修改。

设计一个双层模型, 模型大小为横向4000m, 纵向4000m, 网格大小为横向10m, 纵向10m, 网格数量为横向401个, 纵向401个。模型第一层速度v_P0=3000m/s, 各向异性参数ε=0.1, δ=0.1, 对称轴倾角θ=45°, 第二层速度v_P0=3250m/s, 各向异性参数ε=0.12, δ=0.12, 对称轴倾角θ=25°, 反射界面深度为2000m, 炮点位于地表(2000m, 10m)处, 两组检波器分别放置在深度为1990m和3990m两个层位上，检波器间距均为10m。将TTI声波正演所得单炮记录(黑白同相轴)与相速度法得到的旅行时曲线(绿色曲线)进行叠合, 如图 2所示, 其中图 2a是1990m处的叠合图, 图 2b是3990m处的叠合图。

图 2(Fig. 2)

图 2 旅行时与正演记录叠合显示 a 1990m; b 3990m

由图 2a可以看出, 在第一层均匀模型中, 利用本文提出的基于相速度的射线追踪法得到的旅行时曲线与声波正演模拟产生的单炮记录吻合程度非常高, 说明本文方法在均匀模型中的应用是正确的。由图 2b可以看出, 射线从地表炮点出发, 在两种介质分界面发生了一次透射, 最终到达3990m处的检波点, 利用基于相速度的射线追踪法得到的旅行时曲线与声波正演模拟产生的单炮记录吻合程度同样非常高, 说明该方法在复杂介质中依然准确。

通过二维模型验证方法的精度和效率。程函方程法的精度很高, 但是运算速度较慢, 这里分别采用恒速度模型和变速度模型绘制射线路径和等时线(波前面)图, 并通过等时线与程函方程法进行精度对比, 如果两者等时线重合率高, 说明基于相速度的射线追踪法的精度同样很高; 如果两者等时线出现了较大的偏差, 则需要进一步验证哪种方法精度更高。再采用一个较为复杂的二维模型利用密集射线追踪验证基于相速度的射线追踪法的计算效率, 如果相速度法的效率明显高于程函方程法, 说明相速度法更具有实用性; 如果两种方法的计算效率相当, 则说明相速度法与程函方程法具有相同的实用性。

设计一个二维模型, 模型大小为横向4000m, 纵向4000m, 网格大小为横向10m, 纵向10m, 网格数量为横向401个, 纵向401个。图 3中的模型速度v_P0=2000m/s, 各向异性参数ε=0.3, δ=-0.1, 对称轴倾角θ=45°, 炮点在模型正中间, 坐标为(2000m, 2000m)。图 4中的模型速度v_P0=(2000+0.5z)m/s, 各向异性参数ε=0.3, δ=-0.1, 对称轴倾角θ=45°, 炮点在地表中间位置, 坐标为(2000m, 0)。

图 3(Fig. 3)

图 3 恒速度场(速度2000m/s)射线路径与等时线(a)以及相速度法与程函方程法所得等时线叠合显示(b)

图 4(Fig. 4)

图 4 变速度场射线路径与等时线(a)以及相速度法与程函方程法所得等时线叠合显示(b)

图 3和图 4分别为恒速度场、变速度场射线路径与等时线以及基于相速度的射线追踪法与程函方程法所得等时线叠合显示图。图 3a和图 4a为基于相速度的射线追踪法得到的射线路径和等时线, 背景为速度场。从等时线可以看出, 各向异性情况下波前面不再是圆形, 而是椭圆形, 并且随着对称轴的倾斜发生旋转, 变速情况下射线发生了弯曲, 出现了回转现象。图 3b和图 4b为相速度法和程函方程法得到的等时线叠合图, 其中黑色为程函方程的等时线, 红色为相速度的等时线。因为两种方法的等时线完全重合, 所得到的旅行时信息一致, 所以图中只能看到红色的相速度等时线, 说明相速度法射线追踪的效果和精度跟程函方程法一致, 验证了本文方法的有效性。

采用复杂的二维TTI模型进一步验证本文方法的有效性和运算效率。图 5a为速度场, 图 5b为ε场, 图 5c为δ场, 图 5d为倾角场, 图 5e为射线路径图, 图 5f为等时线图。该模型包含高速盐丘和高陡构造等复杂介质, 射线在浅层出现了回转现象, 在盐丘顶部出现了透射和全反射, 等时线也体现了射线的传播方向和传播速度。

图 5(Fig. 5)

图 5 复杂二维TTI介质的参数场、射线路径和等时线 a速度场; b ε场; c δ场; d倾角场; e射线路径; f等时线

表 1为该TTI模型采用两种方法的计算时间, 射线密度X°/180°表示射线出射角范围是0~180°, 方向向下, 每隔X°出射一条射线, 旅行时单位为s。我们采用3组射线密度(10°/180°, 2°/180°和0.5°/180°)分别测试程函方程法和相速度法射线追踪的用时。从表 1可以看出, 随着射线密度的增大, 两种方法的用时都是成倍数增长, 射线密度越大, 两者用时差距也越大, 当射线密度为0.5°/180°时, 相速度射线追踪法用时比程函方程法少了12.62s, 可见, 相速度法的计算效率明显优于程函方程法, 这种优势在三维模型等大规模数据运算时更为明显。

表 1(Tab. 1)

表 1 3种射线密度下两种方法计算时间对比

表 1 3种射线密度下两种方法计算时间对比

通过模型试算发现, 虽然相速度法射线追踪与程函方程法射线追踪具有相同的精度和效果, 都能够为建模和成像提供准确的射线路径和旅行时信息, 但是相速度法的计算效率远超程函方程法, 能够在很大程度上缩短TTI建模和成像的时间, 比程函方程法更适用于大规模实际数据的处理。

通过三维模型试算验证方法在三维空间中的射线追踪路径和旅行时信息, 以及等时线在空间中的旋转情况。设计一个三维模型, 模型大小为横1000m, 宽1000m, 深度1000m, x, y, z方向网格大小均为10m, x, y, z方向网格数量均为101个, 炮点在模型正中间, 坐标为(500m, 500m, 500m)。

图 6中, 模型速度v_P0=2000m/s, 各向异性参数ε=0.3, δ=-0.1, 对称轴倾角θ=45°, 方位角φ=45°; 图 7中, 模型速度v_P0=(2000+0.5z)m/s, 各向异性参数ε=0.3, δ=-0.1, 对称轴倾角θ=45°, 方位角φ=45°。图 6a和图 7a分别为模型的速度场; 图 6b和图 7b为射线路径图; 图 6c和图 7c为等时线图。由图 6b和图 7b可以看出, 恒速度场中的射线呈直线传播, 变速度场中的射线呈弯曲状, 出现了回转现象; 由图 6c和图 7c可以看出, 波前面不但垂向上旋转了45°(倾角), 水平方向也旋转了45°(方位角)。三维模型测试结果说明相速度射线追踪方法在三维空间中也能取得正确有效的射线路径和旅行时信息。

图 6(Fig. 6)

图 6 三维模型恒速度试算结果 a速度场; b射线路径; c等时线

图 7(Fig. 7)

图 7 三维模型变速度试算结果 a速度场; b射线路径; c等时线

本文在程函方程法TTI介质射线追踪的基础上, 提出了一种简洁、高效的基于相速度的TTI介质射线追踪方法。该方法以Hamilton相速度为基础推导出了射线追踪系统, 并利用群速度公式和出射角公式将其转换为TTI介质出射角格式, 得到基于相速度的TTI介质射线追踪方程, 其公式简单、物理意义明确, 并且所得出参数能够用于层析反演等后续处理。二维、三维模型试算结果表明, 本文提出的相速度法射线追踪的精度较高、效果较好, 能提供准确的射线路径和旅行时信息, 计算效率较高, 非常适用于大规模实际数据处理。

在实际生产中, 射线追踪可用于射线类偏移和建模, 本文所述射线追踪方法只反映透射波、回转波和折射波, 不能够反映反射波, 因此, 在实际使用过程中需要确定反射点位置和反射界面的倾角和方位角, 再遵循反射定律向炮点和检波点进行射线追踪, 从而完成一个炮检对的射线路径模拟。由于无法模拟反射波, 因此也无需考虑边界效应, 但是, 必须考虑全反射现象, 这种现象是我们不想看到的, 为了避免发生这种现象, 在射线追踪之前一般需要将模型进行平滑, 消除界面的突变带来的不稳定性, 再进行射线追踪计算。