1 引言

已有研究表明，地球介质从浅部到深部具有广泛的地震各向异性.复杂山地勘探中，起伏地表是影响勘探效果的一项重要因素.因此，发展考虑起伏地表和地震各向异性的地震波射线追踪方法，对地震波正演模拟、走时层析成像、偏移成像等正反演研究具有重要的理论意义和实际应用价值.

近年来，众多学者不断研究已经将许多射线追踪算法推广到了各向异性介质.比如打靶法(Farra，2005)、最短路径射线追踪法(Zhou and Greenhalgh， 2005； 赵爱华等，2006)、龙格-库塔法(孙福利等，2009)、试射法(李建国等，2010)、波前构建法(白海军等，2011)、线性走时插值方法(Cardarelli and Cerreto， 2002；马德堂等，2011)等等.其中，最短路径射线追踪方法起源于网络理论，最早由Moser(1991)发展成熟，该方法适合任意复杂程度的模型，并且是无条件稳定的.

目前在各向同性地震波走时计算中，起伏地表情况下射线追踪研究开展得较多.如利用三次样条插值函数拟合起伏地表，实现了初值法射线追踪(岳玉波等，2007)；采用块状结构建立复杂地质体模型实现地震波的快速追踪(徐涛等，2004；吴军和贾雨，2008)；采用非规则网格(梯形网格和三角网格)拟合地表，结合分区多步计算实现复杂层状起伏介质的走时计算(赵瑞和白超英，2010；黄国娇和白超英，2010)等.各向异性地震波射线追踪的研究仍然以简单的层状模型居多(Qin and Schuster， 1993；Zhou and Greenhalgh， 2006；Jiang and Zhou， 2011)，也有少数研究(李晓玲等，2010)考虑到了起伏地表的情况.

本文将最短路径追踪方法由传统规则网格发展到非规则网格，实现了非规则网格和规则网格混合最短路径射线追踪计算，在带起伏地表复杂的各向异性介质模型射线追踪实算中取得了较好效果.

由于考虑了起伏地表，本文研究需要选取合适的网格处理方式.贴体网格技术是各向同性介质领域计算精度比较高的处理方法之一，但是在各向异性介质中，随着贴体网格的各向异性增强，用坐标变换法求解地表起伏区域的走时计算误差增大且计算效率降低(兰海强等，2012).因此，本文算法中采用的是正方形网格与不规则网格构成的混合网格.

首先将整体模型用网格离散化成一系列正方形小单元，每个单元具有独自的Thomsen参数值(v_p0、δ、ε)，再根据边界切割正方形网格的情况，将计算单元分为两种：正常网格(未被边界线穿过的网格)和边界网格(被边界线穿过的网格).

对于正常网格，在小单元每个边界等间距设置K个(不含端点)节点，同时在小单元四个角点也设置节点(上述两类节点统称正常节点)，并约定子波可在任意两个能够直接到达的节点之间传播.这种设置节点的方法既方便编程计算，又可使子波源(二次源)射线传播方向全方位覆盖且空间分布较均匀.为了保证计算精度，K的取值越大越好，但同时也会不可避免地降低计算效率.模型实验的结果(图 1)表明，K≥7时平均相对误差已经降至1‰以下，可以提供较高的计算精度.因此，本文的算例统一取K=9.

图 1

Fig. 1

图 1 K的取值与计算精度和效率的关系示意图 Fig. 1 A sketch to illustrate the relationship between the value of K and accuracy as well as efficiency

本文算法根据节点所处的位置，计算中将子波震源划分成行源、列源、角源三类，其全方位子波传播路径如图 2所示(张美根等，2006；Zhang et al，2009).

图 2

Fig. 2

图 2 正常网格行源(a)、列源(b)、角源(c)传播路径示意图 Fig. 2 A sketch to illustrate the travel path from row source(a)，column source(b) and corner source(c)

对于边界网格，可视为起伏边界切割正方形网格构成的不规则网格.因此，不规则网格设置节点时，只需在正方形网格原有节点的基础上增加设置边界线上的节点(简称边界节点)即可.理想情况下，边界网格内的边界线可以是一条任意的曲线段.但是，考虑到实际勘探中道间距较小，可将边界网格内的地形边界视为线性变化，用一条直线段来近似表 示(图 3).为了描述方便，我们称这条直线段为局部拟合线，局部拟合线依次连接就形成了拟合地形.实际资料处理中，取定基准面，将采集的高程值先简单换算成深度值(Z值)，再进行插值计算，按照网格间距取出水平距离等间距变化的数据点，相邻两点间的连线即为局部拟合线.

图 3

Fig. 3

图 3 边界网格内的节点设置与部分子源传播路径示意图 Fig. 3 A sketch to illustrate the node settings and travel path from some sub-sources in boundary cells

在局部拟合线上设置节点时，要保证子波源射线具有足够的出射方向且出射方向空间分布较均匀，实际计算中以正常网格(正方形网格)上子波源的正常射线方向为参考标准.以图 3局部拟合线P₁P₂为例，当子源是角源S₁时，参照正常网格的角源全方位子波传播路径(图 2)可以算出局部拟合线上应该设置的节点数为2*K+1(不含端点).这种情况下，局部拟合线上的节点密度已经大于正常网格边界上的节点密度，走时计算精度可以得到保证.虽然存在局部拟合线上设置节点数更多的情况，比 如子源位于S₂局部拟合线对应的节点数为3*K+1(不含端点)，但是实际勘探中道间距较小，计算中划分的网格间距也不大，设置的接收点通常也是位于网格的边界，所以局部拟合线上设置过多的节点对于提高接收点初至到时的计算精度作用不大.

另一方面，当局部拟合线切割正方形小单元一角且只有角点上的节点位于边界外侧时，如图 3中局部拟合线P₂P₃，局部拟合线的长度非常短，设置1个节点(不含端点)已经使得局部拟合线上的节点密度大于正常网格边界上的节点密度，再设置2*K+1个节点是没有必要的.所以，本文算法对于这种情况只在局部拟合线上设置3个节点(含端点)，剩余情况设置2*K+3个节点(含端点).节点设置的具体实现过程如下：①判断局部拟合线是否切割正方形网格一角，若是执行②，否则设置2*K+3个节点；②判断步骤①中形成的三角形直角边长度是否均小于λ，λ=小单元边长/(K+1)，若是设置3个节点，否则设置2*K+3个节点.

走时计算中若边界节点与正常节点重合，重合节点只计算一次；若局部拟合线与正方形单元边界的一部分重合，则重合部分不再设置正常节点.同样约定子波可在边界网格内任意两个能够直接到达的节点之间传播，部分边界节点的全方位子波传播路径见图 3.

由于最短路径算法中每段射线路径都是由两节点间的连线近似表示的，而正常网格的节点分布是确定的，因此计算中可以提前将正常网格内二次源所有可能存在射线段的长度计算出来存于数组中，以节省走时计算的耗时.但是边界网格的形状不会完全相同，所以需要提前计算所有边界网格内的所有可能存在射线段传播的局部走时.当同一模型设置的炮点很多时，这种做法可节省大量运算时间.

从震源点开始，波首先到达与震源点相邻的计算节点.从这些节点中选出走时最小的一点作为子波源，波再传向与子波源相邻的计算节点，再选出下一个子波源，重复以上过程就实现了各个节点最小走时的计算.追踪射线路径时，根据互易原理，从接收点开始依次确定前一级子波源直到震源点结束，将各点顺序依次连接就得到射线路径.这就是最短路径射线追踪算法的基本原理.

在计算最小走时时，需要将所有计算节点(包括正常节点和边界节点)T分成四部分：①已经确定最小走时的节点，即已经做过子震源的节点集合S；②已计算过至少一次局部最小走时，但未确定最小走时的节点集合P；③尚未计算走时的节点，即所有节点中除去S、P后剩余节点的集合R；④与震源(或子震源)直接相邻的节点，即每次迭代中需要计算局部最小走时的节点集合N.实现过程如下：

(1)初始化及初次计算走时.置S、P为空集，R=T.将震源走时取为0，计算震源点到相邻节点的走时，震源点放入S中并从R中剔除，算过走时的节点放入P中从R中剔除.

局部最小走时的计算公式如下：

(2)循环计算.从集合P中选取走时最小的点作为子源i，根据公式计算得到子源点i到相邻节点j的局部走时.若j∈R，t(j)=t(i)+ l_ij/V_p，将节点j 放入P中并从 R 中剔除；若j∈P，t(j)=MIN{t(j)，t(i)+ l_ij/V_p} .

(3)循环判断.如果S=T或者 R 为空集时循环结束，否则重复进行步骤(2).

在上述循环计算过程中，需要不断选取走时最 小点作为下一个子源，这是算法最耗时的一步.对于无序序列，快速排序算法是平均性能最好的算法；而 对于正序序列，插入排序算法的耗时可以达到O(n)量级(常旭和刘伊克，2001；张建中等，2003).所以，本算法在初次计算走时时利用快速排序法将波前点排序，之后增加波前点时采用二分-插入排序法进行排序，以提高波前点排序的效率.

在各向异性介质中，相速度与群速度往往是不同的，并且推导准确解析解相当复杂，不适合进行数值计算.1986年，Thomsen提出用一些各向异性参数描述TI介质，并给出了速度的计算公式.在计算走时过程中直接使用到的是群速度，众多学者不断深入研究群速度的近似公式(Byun et al，1989；Sena，1991；赵爱华和丁志峰，2005；吴国忱等，2009)，使得数值计算变得更加方便.在弱各向异性假设下，P波相速度可近似表示为(Thomsen，1986)：

群速度可由相速度表示(Thomsen，1986)：

在弱各向异性线性近似下，群角与相角满足如下关系(Thomsen，1986)：

在具有倾斜对称轴的横向各向同性介质(TTI介质)中，介质对称轴与z轴夹角一般不为0，记作φ₀.定义φ₀的取值范围是[-90°，90°]，且约定从z轴逆时针旋转到对称轴时φ₀取正，反之取负.射线传播方向与z轴的夹角记作η，且约定从z轴逆时针旋转到射线传播方向时η取正，反之取负.则TTI介质中的群速度表示为(Jiang and Zhou， 2011)：

图 4

Fig. 4

图 4 TTI介质群角示意图 η表示射线传播方向与z轴的夹角，φ0表示 对称轴与z轴的夹角，η-φ0即为群角大小. Fig. 4 A sketch to illustrate the group angle in TTI media η is the angle between ray path and z axis，φ0 is the angle between symmetry axis and z axis and η-φ0 is the size of group angle.

为了检验本文算法的准确性，我们以一个规则的均匀各项异性模型(v_p0=1000 m·s^-1、δ=0.1、 ε=0.2、φ₀=0°)计算理论走时与追踪走时的相对误差.模型规模1000 m×1000 m，网格间距5 m，K=9，炮点位于(500 m，500 m)处，接收点均匀分布于模型外边界，计算结果如图 5所示.其中，最大相对 误差0.24631845%，最小相对误差仅为0.00000842%，表现出了很高的计算精度.

图 5

Fig. 5

图 5 本文算法计算误差 Fig. 5 Calculation error of the algorithm in this paper

图 6是传统方法(阶梯状地表)与本文方法计算的初至到时之差.计算模型为带有起伏地表的均匀各向异性模型(v_p0=2000 m·s^-1、δ=0.15、ε=0.2、 φ₀=30°)，炮点位于地下点(0 m，50 m)处，在1000 m 长的测线上等道间距设置201个接收点，起伏地表形态见图 7.除了传统方法中阶梯状边界引起的计算接收点位置与真实接收点位置偏差外，计算中其他各项参数保持完全相同.

图 6

Fig. 6

图 6 传统方法与本文方法计算的初至到时之差 Fig. 6 The difference between the first arrivals calculated by traditional method and the method in this paper

图 7

Fig. 7

图 7 复杂各向异性介质模型I Fig. 7 Complex anisotropic medium model I

图中整体趋势可以看出，本文方法相比传统方法追踪的初至到时更小.当然，图中也存在传统方法计算初至到时更小的个别点.这种情况是传统方法中计算接收点的位置低于真实边界引起的，并不代表走时计算更加准确.传统方法中，起伏边界用阶梯 状的边界来近似，这就不可避免地导致部分或全部计算接收点与真实接收点位置不重合，引入计算误差.本文方法的近似边界是局部拟合线依次连接形成的，局部拟合线两端处的深度值可由实际测量的高程数据转换而来，而算例中的所有接收点均位于局部拟合线端点处，因此，本文方法的计算接收点与真实接收点位置完全重合，从而避免了边界拟合引入的接收点位置偏差，使得初至到时的计算精度得到提高.

进一步，我们采用两个带有起伏地表的复杂各向异性介质模型(图 7，图 8)来验证本文的可靠性和对起伏界面的处理能力.模型I规模1000 m×1000 m，网格间距5m，网格小单元每个边界上设置K=9个节点.该模型由五部分介质组成，各部分介质的具体参数见表 1.模型II是带有起伏地表的各向异性Marmousi模型，模型规模3840 m×1220 m，网格间距10 m，网格小单元每个边界上设置K=9个节点.该模型的v_p0变化范围是1500~5500 m·s^-1(深颜色代表低速)，其它参数见表 2.

图 8

Fig. 8

图 8 复杂各向异性介质模型II Fig. 8 Complex anisotropic medium model II

表 1(Table 1)

表 1 各向异性介质模型I参数 Table 1 The parameters of anisotropic medium model I vp0/(m·s-1) δ ε φ0/(°) 介质1 3000 0.2 0.15 -30 介质2 3400 0.15 0.3 30 介质3 2000 0.25 0.05 -45 介质4 4000 0.1 0.1 45 介质5 4600 0.1 0.25 15

表 1 各向异性介质模型I参数 Table 1 The parameters of anisotropic medium model I

表 2(Table 2)

表 2 各向异性介质模型II参数 Table 2 The parameters of anisotropic medium model II vp0/(m·s-1)范围 δ ε φ0/(°) [1500，2500) 0.1 0.1 45 ［2500，3500) 0.25 0.05 15 ［3500，4500) 0.15 0.3 30 ［4500，5500] 0.2 0.15 -30

表 2 各向异性介质模型II参数 Table 2 The parameters of anisotropic medium model II

图 9和图 10是按照本文算法追踪到的炮点到边界点的初至波射线路径和初至波等时线图，模型I中炮点位于地表点(440 m，46 m)，模型II中炮点位于地表点(1260 m，61 m).从两图结果可以看出，本文方法在模型峰谷位置都获得了准确合理的走时路径和圆滑的走时等值线.在模型内部射线表现出了“趋向高速介质，避离低速介质”的特征，符合射线传播的基本规律，证实本文方法是可靠的.模型II的追踪结果也表明，本文方法对于横向和纵向均变化剧烈的复杂各向异性模型同样具有很好的适用性.

图 9

Fig. 9

图 9 模型I边界点初至射线路径(a)和等时线图(b) Fig. 9 Raypaths to boundary nodes(a) and isochron map of the first arrivals(b)in model I

图 10

Fig. 10

图 10 模型II边界点初至射线路径(a)和等时线图(b) Fig. 10 Raypaths to boundary nodes(a) and isochron map of the first arrivals(b)in model II

本文将最短路径追踪方法由传统规则网格发展到非规则网格，实现了混合网格各向异性地震波最 短路径射线追踪算法.在正方形网格离散化的基础上，通过合理的地形近似拟合实现了各个节点初至走时的快速计算，节省了计算时间.最短路径射线追踪算法的模型实算表明，该算法适合复杂的二维各向异性介质，射线路径追踪结果符合射线传播的基本规律，显示了很好的应用前景，有望推广到三维复杂各向异性介质走时计算中.