0 引言

地下岩石普遍存在各向异性特征，地震波在实际地层中的传播速度不仅与空间位置相关，还受传播方向的影响^[1-2]，在一定程度上制约了地震勘探的分辨率和储层预测能力^[3]。因此揭示地震波在各向异性介质中的传播规律有着重要意义，常用方法之一就是地震射线追踪。

研究各向异性介质中地震波的传播过程常采用横向各向同性(Transversely Isotropic，TI)介质假设。Postma^[4]证明当薄层中的单层厚度远小于地震波波长时，周期性薄层组合结构可以等效为均匀TI介质。TI介质可由5个独立弹性参数表示，Thomsen^[5]采用P波和S波的垂向速度以及三个无量纲的各向异性参数替代5个弹性参数，具有更清晰的物理意义。计算TI介质相速度或群速度的研究较为成熟^[6-7]，对应的射线追踪方法也较为丰富^[8-9]。

沉积盆地大多具有薄互层特征，TI介质可以较好地表达其各向异性属性^[4]。然而，在构造应力或人工压裂作用下，地层中容易发育垂向或近似垂向的高角度裂缝组，从而表现出明显的正交各向异性(ORT)特征^[10-11]。目前，对于ORT介质中qP波(或声波)速度各向异性研究已取得一些成果：Song等^[12]给出了弱正交各向异性介质中声波群速度计算公式；Alkhalifah^[13]在声波近似条件下，采用6个独立参数计算ORT介质中P波慢度，简化了计算过程；Abedi等^[14]讨论了一种新的ORT介质中计算声波相速度和群速度的近似公式，具有较高计算精度；Tsvankin^[15]根据TI介质和ORT介质对称平面上Christoffel方程的相似性，参照Thomsen参数的形式推导了一种适合于地震速度分析和反演的ORT介质参数化方法，同时还给出了弱正交各向异性介质中qP波相速度的近似表达式；Hao等^[16]研究了倾斜正交各向异性(Tilted Orthorhombic，TORT)介质中声波慢度曲面。

地震P波旅行时对于近地表建模、微地震定位、裂缝储层特征研究等有着十分重要的意义。目前，ORT介质中地震射线追踪研究多限于计算初至qP波旅行时和射线路径。Mensch等^[17]利用椭球各向异性介质作为参考介质，结合扰动理论求取ORT介质中qP波旅行时及射线路径；Pšenčík等^[18]在弱正交各向异性介质假设下讨论了qP波射线路径及旅行时的计算问题；Ibanez-jacome等^[19]采用快速行进法计算ORT介质中的声波旅行时场；为高效计算三维多层ORT介质中射线路径和旅行时，Blias^[20]依据费马原理采用牛顿法求解射线与弯曲界面的交点，进而实现两点间P波射线追踪，随后又将该方法应用于TORT介质^[21]；Masmoudi等^[22]将P波旅行时近似表示为背景介质中多个各向异性参数(扰动参数)的函数，进而利用扰动理论计算旅行时场。对于反射波，卢明辉等^[23]计算了单层水平反射界面条件下ORT介质模型中全方位角的反射P波旅行时；Abedi等^[24]利用新的声波近似公式求取ORT介质的速度分布，并计算了简单模型中的反射波射线路径。

由于实际地层的结构较为复杂，要求相应的射线追踪算法有较强的适应性，并兼具较高的稳定性和计算效率。针对这一需求，本文基于分区多步改进型最短路径算法，研究了一种TORT介质中计算多次反射qP波射线路径和旅行时的射线追踪算法。在求取速度分布的过程中，引入球面三角剖分算法^[25]确定近似均匀分布的相速度方向采样点，可在保证采样精度的同时减少采样点数，从而降低计算量。为进一步提高计算效率，采用GPU并行计算技术求取各个采样方向上的群速度。

1 TORT介质中分区多步最短路径算法 1.1 TORT介质的群速度计算方法

高频地震波的传播可由Christoffel方程描述^[1]

$ \operatorname{det}\left(\mathit\Gamma_{j k}-G{\delta_{j k}}\right)=0 $

(1)

式中：δ_jk为克罗内克函数；Γ和G分别为Christoffel矩阵及其特征值。在指定的相速度方向n上，有

$ \mathit\Gamma_{j k}(\boldsymbol{n})=a_{i j k l} n_i n_l $

(2)

$ G(\boldsymbol{n})=a_{i j k l} n_i n_l g_j g_k=c^2 $

(3)

式中：a_ijkl(i，j，k，l=1，2，3)为密度归一化的四阶弹性张量系数；c为相速度；g为归一化特征向量(偏振方向)。对于一般各向异性介质，存在三组特征值和特征向量，相速度最大的一组对应qP波，其余两组分别为qS1和qS2波。相应的群速度可表示为

$ v_i=a_{i j k l} p_l g_j g_k $

(4)

式中p=n/c为慢度矢量。ORT介质可由9个Voigt形式的独立参数{a₁₁，a₁₂，a₁₃，a₂₂，a₂₃，a₃₃，a₄₄，a₅₅，a₆₆}表示，将之代入式(1)~式(4)即可求出水平层状ORT介质中的群速度。

受构造运动影响，实际地层大多并不水平。为计算TORT介质中的群速度，可按照地层倾斜情况将群速度矢量进行一定角度的旋转。若以θ₀、φ₀、α分别表示地层的倾角、倾向以及裂缝的方位角，则旋转矩阵为

$ \begin{aligned} \boldsymbol{M} & =\left(\begin{array}{ccc} \;\;\;\cos \alpha & \sin \alpha & 0 \\ -\sin \alpha & \cos \alpha & 0 \\ 0 & 0 & 1 \end{array}\right)\left(\begin{array}{ccc} \cos \theta_0 & 0 & -\sin \theta_0 \\ 0 & 1 & 0 \\ \sin \theta_0 & 0 & \cos \theta_0 \end{array}\right)\left(\begin{array}{ccc} \;\;\;\cos \varphi_0 & \sin \varphi_0 & 0 \\ -\sin \varphi_0 & \cos \varphi_0 & 0 \\ 0 & 0 & 1 \end{array}\right) \\ & =\left(\begin{array}{c} \;\;\;\cos \alpha \cos \theta_0 \cos \varphi_0-\sin \alpha \sin \varphi_0 & \;\;\;\cos \alpha \cos \theta_0 \sin \varphi_0+\sin \alpha \cos \varphi_0 & -\cos \alpha \sin \theta_0 \\ -\sin \alpha \cos \theta_0 \cos \varphi_0-\cos \alpha \sin \varphi_0 & -\sin \alpha \cos \theta_0 \sin \varphi_0+\cos \alpha \cos \varphi_0 & \;\;\;\sin \alpha \sin \theta_0 \\ \sin \theta_0 \cos \varphi_0 & \sin \theta_0 \sin \varphi_0 & \cos \theta_0 \end{array}\right) \\ & =\left(\begin{array}{lll} M_{11} & M_{12} & M_{13} \\ M_{21} & M_{22} & M_{23} \\ M_{31} & M_{32} & M_{33} \end{array}\right) \end{aligned} $

(5)

变换后的群速度可写为

$ \widetilde{v}_i=M_{i j} v_j $

(6)

式中$\widetilde{v}_i$和v_j分别为TORT和ORT介质中的群速度分量。

1.2 相速度方向近似均匀采样

各向异性介质中，采用最短路径法计算旅行时和射线路径时，需要已知各主节点处多个射线出射方向的群速度。目前直接求取指定射线方向的群速度较为复杂，通常的做法是按照一定间隔对相速度方向进行全方位采样(ORT介质仅需在单个象限内采样)，求出所有采样点对应的射线方向及其群速度，随后插值得到指定射线方向上的群速度。为便于计算，常采用倾角θ和方位角φ指定相速度的方向

$ \boldsymbol{n}=(\sin \theta \cos \varphi, \sin \theta \sin \varphi, \cos \theta)^{\mathrm{T}} $

(7)

一般做法是按照等角度间隔对θ和φ进行采样，采样点较多且分布极不均匀(θ值越小，采样越密集)。本文引入基于正八面体的球面三角剖分算法^[25]确定相速度方向采样点，在保证采样精度的前提下，采样点数量可减少一半左右，且分布更均匀，在很大程度上提高了计算效率。图 1a~图 1c给出了一个象限内的球面剖分示意图，其步骤如下：①将三个顶点(分别为(0°，0°)、(90°，0°)、(90°，90°))所确定的采样区域(单个象限)围成球面三角形(图 1a)；②寻找球面三角形中每条边的中点，并将任意两个中点相连接，从而使原来的球面三角形划分为4个新的球面三角形(图 1b)；③依据步骤②对每一个新的球面三角形进行再次剖分，直至剖分精度满足要求，此时每个采样区域被划分为许多小的球面三角形(图 1c)。

依次将所有采样点对应的坐标(θ，φ)代入式(7)，并结合式(2)~式(6)即可求得ORT或TORT介质的群速度分布。对于特定射线方向，其群速度可在包含该方向的球面三角形内通过插值求得。如图 1d所示，球面三角形ABC内P点对应的群速度，可由其所在球面三角形三个顶点处的值内插得到^[26]

$ v_P=\frac{v_A S_{P B C}+v_B S_{A P C}+v_C S_{A B P}}{S_{A B C}} $

(8)

式中S为对应球面三角形的面积。

求取ORT(或TORT)介质的群速度分布属于计算密集型任务，且各个相速度采样方向上的计算是独立的，为进一步提高计算效率，采用了GPU并行计算。

1.3 分区多步改进型最短路径射线追踪算法

最短路径算法^[27]广泛应用于追踪地震波射线路径及计算旅行时，Zhou等^[8]证明了最短路径算法在一般各向异性介质中的适用性，并在TI介质中追踪多种类型的地震波。唐小平等^[28]通过在主节点之间插入次级节点的方式对最短路径算法进行改进，并结合分区多步计算技术，提出了分区多步改进型最短路径算法(Modified Shortest-path Method with Multistage Scheme，MSPM)用以追踪多次反射或反射转换波。本文亦应用该算法计算TORT介质中多次反射qP波的射线路径及旅行时，主要步骤如下：

(1) 采用规则网格进行模型参数化，计算每个网格点(主节点)上的群速度分布；

(2) 按照界面分布形态，将模型划分为多个独立的计算区域；

(3) 由震源出发进行波前扩展，计算波前上每个网格中的最小旅行时，直至波前到达相应的反射界面；

(4) 仅保留界面上的旅行时，根据地震波类型，选择下一个计算区域进行波前扩展；

(5) 重复步骤(4)直至波前到达检波器；

(6) 提取检波器处旅行时及射线路径信息。

计算每个网格中的最小旅行时(步骤(3))，是最短路径算法的核心。在合适的网格间距下，相邻主节点间的速度差异不大，波前点I附近任意节点J处的旅行时可由以下公式(两点间距离除以平均群速度)近似计算

$ t_J=\min _\limits{I \in \boldsymbol{C}_I}\left[t_I+\frac{2 D\left(\boldsymbol{X}_I, \boldsymbol{X}_J\right)}{v\left(\boldsymbol{X}_I, \boldsymbol{r}_{I J}\right)+v\left(\boldsymbol{X}_J, \boldsymbol{r}_{I J}\right)}\right] $

(9)

式中：J为待求节点序号；I是波前上的次级震源节点序号；C_I为节点J的邻近次级震源(节点)的集合；X_I和X_J分别表示节点I和J的坐标；t_I是自炮点到达节点I的最小旅行时；D(X_I，X_J)是节点I与节点J之间的距离；r_IJ为射线方向；v(X_I，r_IJ)和v(X_J，r_IJ)分别是节点I和节点J处r_IJ方向的群速度。

2 模型测试 2.1 TORT介质速度分布

为更直观表征ORT介质弹性参数的物理特征，Tsvankin^[15]将其转化为两个垂向速度和七个无量纲的参数

$ \left\{\begin{array}{l} V_{\mathrm{P} 0}=\sqrt{\frac{c_{33}}{\rho}} \\ V_{\mathrm{S} 0}=\sqrt{\frac{c_{55}}{\rho}} \end{array}\right. $

(10)

$ \left\{\begin{array}{l} \varepsilon_1=\frac{c_{22}-c_{33}}{2 c_{33}} \\ \delta_1=\frac{\left(c_{23}+c_{44}\right)^2-\left(c_{33}-c_{44}\right)^2}{2 c_{33}\left(c_{33}-c_{44}\right)} \\ \gamma_1=\frac{c_{66}-c_{55}}{2 c_{55}} \end{array}\right. $

(11)

$ \left\{\begin{array}{l} \varepsilon_2=\frac{c_{11}-c_{33}}{2 c_{33}} \\ \delta_2=\frac{\left(c_{13}+c_{55}\right)^2-\left(c_{33}-c_{55}\right)^2}{2 c_{33}\left(c_{33}-c_{55}\right)} \\ \gamma_2=\frac{c_{66}-c_{44}}{2 c_{44}} \end{array}\right. $

(12)

$ \delta_3=\frac{\left(c_{12}+c_{66}\right)^2-\left(c_{11}-c_{66}\right)^2}{2 c_{11}\left(c_{11}-c_{66}\right)} $

(13)

式中：ρ为密度；c_ij=ρ·a_ij为归一化前的弹性参数；V_P0、V_S0分别表示qP波、qS波的垂向速度；ε、δ、γ与在VTI介质中Thomsen^[5]参数定义相同，下标“1、2、3”分别表示yz、xz和xy平面。

为了展示ORT (或TORT)介质的群速度分布特征，表 1给出了三组各向异性介质参数：模型M1为ORT介质^[15]，模型M2是模型M1经过一定角度旋转得到的TORT介质，模型M3则是由Schoenberg等^[29]给出的标准正交各向异性介质。表 1中介质参数对应的密度归一化弹性参数如表 2所示。

图 2a所示为模型M1 (ORT介质)对应的群速度分布，可以看出，该分布关于三个正交平面对称，在第一个象限内，群速度最小值所在方向大致为(θ=30°，φ=0°)，而群速度最大值所在方向大致为(θ=90°，φ=90°)。模型M2 (TORT介质)中各向异性参数与模型M1完全一致，仅对其旋转了一定角度，因此对应的群速度分布(图 2b)等同于将图 2a旋转相应的角度。

图 3给出了xy、xz、yz三个正交平面内模型M1(实线)和模型M2 (虚线)的群速度分布曲线。因为ORT介质的三个对称面与上述三个正交平面重合，平面内性质等同于VTI介质，故模型M1的群速度曲线(图 3中实线所示)具有对称的四象限分布特征；而对于TORT介质，其对称面与上述三个正交平面并不重合，因此模型M2的群速度曲线(图 3中虚线所示)并没有此对称现象。

图 4a和图 4b所示分别为模型M1和模型M3中1 s时刻的波前，波前面到震源的距离R能够清晰反映地震波速度在各个方向的差异。总体看来，由于模型M1的垂向P波速度(V_P0参数)较大，群速度更高，其波前面范围也就更大。

2.2 TORT介质中地震射线追踪

数值模拟的模型尺寸为5.0 km×3.0 km×2.5 km，采用0.1 km×0.1 km×0.1 km的正方体网格单元进行剖分，并在主节点之间插入9个次级节点(次级节点间距0.01 km)，震源置于(0，0，0)处(图 5中红色圆点)。

为检验算法的计算精度，首先采用一个特殊的均匀ORT介质模型(表 2中模型M4)，该模型中a₁₂=－a₆₆，a₁₃=－a₅₅，a₂₃=－a₄₄，其直达qP波旅行时的解析解为^[30]

$ t=R \sqrt{\frac{N_1^2}{a_{11}}+\frac{N_2^2}{a_{66}}+\frac{N_3^2}{a_{55}}} $

(14)

式中N为射线(群速度)方向。图 5a所示为模型M4的直达波旅行时的相对误差，最大约为0.34%。由于模型M4较为特殊，为不失一般性，采用均匀M1模型再次计算直达qP波旅行时场。图 5b为计算结果与精确解(由距离除以射线方向上的群速度给出)的相对误差，可以发现，当射线出射方向接近网格对角线时误差稍大，但最大不超过0.35%。误差分析结果表明本文所采用的射线追踪算法具有较高的计算精度。

第二个例子为计算反射qP波旅行时和射线路径。在上例所采用模型M1的基础上，于模型底部设置一个水平反射界面，17个检波器分别位于两条相互垂直的测线上(图 6中灰色圆点)，所得反射波射线路径如图 6所示。由于采用了均匀模型M1 (ORT介质)，反射波射线路径由两条直线段构成。图 7a为直达波旅行时等值线图，其yz平面内不同方向上旅行时场等值线疏密程度不同，表明介质具有各向异性特征；图 7b为反射后的旅行时场等值线图，由于采用MSPM算法求取反射波旅行时，将反射界面(底界面)视为次级震源，因此可以发现图 7a和图 7b中xy平面内旅行时等值线图完全一致。

分区多步最短路径算法亦可追踪多次反射波。在第二个例子的模型基础上，于1.0~1.5 km深度增加一个倾斜反射界面(此时模型中包含两个反射界面，图 8)，将模型划分为上、下两个区域：上层为ORT介质(模型M3)，V_P0沿深度线性增加，即V_P0=0.2z+2.437；下层为TORT介质(模型M2)，V_P0沿深度也线性增加，即V_P0=0.2z+3.0。图 8为本文方法追踪的多次波射线路径，可以看到，地震波在两个界面之间进行了多次透射和反射，射线路径由6条直线段组成。该方法亦可用于计算更复杂的地震波类型。

2.3 计算效率分析

计算ORT/TORT介质的群速度分布较为耗时，在整个射线追踪过程中占较大的比重。本文采用球面三角剖分算法确定近似均匀的相速度方向采样点，能够降低一半左右的计算量。同时采用GPU并行计算求取各个采样方向上的qP波群速度值，可进一步减少计算时间。以图 8为例，考察CPU串行与GPU并行的运行时间。CPU型号为Intel Core i7-9700 K (3.60 GHz)，GPU型号为NVIDIA Quadro P4000。CPU串行程序的总运行时间为1161 s，其中计算群速度耗时706 s (约为总耗时的61%)，三维射线追踪耗时455 s；采用GPU加速后总的计算耗时缩减至481 s，其中计算群速度耗时23 s (约为总耗时的4.8%)，三维射线追踪耗时458 s，说明采用GPU并行计算ORT介质的群速度时效率约为使用CPU的30.7倍。该例中需要追踪的地震波类型较为复杂(射线由6段组成)，因此射线追踪的耗时较长。

3 结论

本文首先讨论了ORT/TORT介质中qP波群速度的计算方法，结合分区多步改进型最短路径算法实现了多次波射线追踪。为提高计算效率，引入基于正八面体的球面三角剖分算法确定相速度方向采样点，在精度大致相同的情形下，可减少大约一半采样点。同时采用GPU并行计算进一步减少了耗时。数值计算实例表明，该算法可追踪复杂的多次波；精度分析结果显示，直达波旅行时计算误差不超过0.35%，表明该算法具有较高的计算精度；计算效率分析结果表明，计算ORT/TORT介质的群速度是一个较为耗时的过程，采用GPU并行效率可提高近30倍。

本文着重讨论ORT/TORT介质中多次反射qP波的追踪问题，但其思想也适用于一般各向异性介质。由于qS1和qS2波的波前存在奇异点，将二者完全分离不易实现，因此本文仅对qP波进行了讨论。发展同时适用于三种波的多次反射或反射转换波射线追踪算法是今后的研究方向。