2017, Vol. 60
2. 国土资源部应用地球物理综合解释理论开放实验室—波动理论与成像技术实验室, 长春 130026;
3. 重庆市地质矿产勘查开发局208水文地质工程地质队, 重庆 408300
2. Laboratory for Integrated Geophysical Interpretation Theory of the Ministry for Land and Resources of China-Laboratory for Wave Theory and Imaging Technology, Changchun 130026, China;
3. Chongqing Bureau of Geology and Minerals Exploration 208 Hydrogeological & Engineering Geology Brigade, Chongqing 408300, China
三维复杂山地地质条件下地震波的传播规律非常复杂，通常有很多种类型的地震波同时存在.对三维复杂山地条件下的多波型走时进行计算，有利于分析三维复杂山地条件下地震波的运动学特征.同时，地震波走时计算技术还可以为地震层析成像 (井西利等，2007)、偏移成像 (Sun，2000)、走时反演 (瞿辰等，2007；白超英等，2011；黄国娇和白超英, 2013;黄国娇等，2015)等数据处理技术提供有利的正演模拟工具.因此，本文将围绕三维复杂山地条件下的多波型走时计算问题展开.
目前，地震走时计算方法主要有两点射线追踪法 (Julian and Gubbins, 1977)、最短路径射线追踪法 (Nakanishi and Yamaguchi, 1986；王辉和常旭，2000；Bai et al., 2007)、波前构建法 (Vinje et al., 1993)、走时插值法 (Asakawa and Kawanaka, 1993；Wang and Ma, 1999)、有限差分法 (Vidale, 1988, 1990；Sethian and Popovici, 1999)等.其中，两点射线追踪法的优点是可处理多值问题，而且计算精度通常较高，但其计算效率很有限，不易处理三维空间计算量大的问题，且还存在盲区问题；最短路径射线追踪法将图形理论引入射线追踪中，方法新颖且具有很高的走时计算精度，但因为需要设置很多网格线上的网格节点，所以其在处理三维空间时需要花费很大的内存空间；波前构建法能同时计算走时和射线路径，且能处理多值走时问题，但由于需要做各种插值问题和网格定位问题，因此其在三维空间非常复杂.综合考虑，走时插值法近些年在三维空间内被广泛应用 (张东等, 2009, 2013；梅胜全等，2010；刘锋等，2012；李培明等，2013；孙章庆等，2015).但是却很少有基于这些算法研究三维复杂山地条件下的多波型走时计算问题的.实际上，三维复杂山地条件下的多波型走时计算是一个非常复杂的问题，其待解决的问题主要包括：①三维复杂地表模型的建立与网格剖分；②三维复杂地表条件下地震波走时计算公式建立；③三维复杂地表条件下地震波走时计算整体实现步骤；④三维复杂地表条件下各种不同波型走时计算实现策略.关于问题①，实际上是一个采用怎样的网格建立复杂山地数学模型的问题，该网格不仅需要能准确刻画复杂地表形态，同时还能满足后续问题②③④的计算需要；关于问题②，实际上涉及的是算法的复杂程度和计算精度的问题；关于问题③④，实际上涉及的是算法的灵活性、适应性及整体实现步骤的问题.
针对上述问题①—③，白超英等提出了三维层状介质中多次波追踪的最短路径算法 (白超英等，2011；唐小平和白超英，2009)、兰海强等提出了基于曲网格的复杂地表条件下的地震波走时计算方法 (Lan and Zhang, 2013a, 2013b; Lan et al., 2014)；孙章庆等(2012) 也曾采用不等距网格剖分三维复杂地表模型，采用不等距迎风差分法简洁地建立三维复杂地表条件下的地震波走时计算公式，采用快速推进法作为算法的整体实现步骤.该算法具有计算简洁和灵活的优势，但是相对走时插值法其计算精度相对有限，且无法处理多波型走时计算的问题.孙章庆等(2015) 也曾研究过三维水平地表条件下的双线性插值法，该算法具有很好的精度.但是，由于受限于其公式推导的基本原理，该算法相对有限差分法则不易简洁地处理复杂地表问题.鉴于此，为了简洁地处理三维复杂地表问题，并保证算法的计算精度，在此提出一种迎风混合法.该方法在地表及附近采用简洁的不等距迎风差分法，而在除地表以外的其他区域采用计算精度更高的迎风双线性插值法.综合分析，该迎风混合法既能用不等距网格准确定位地表高程位置，也能采用不等距迎风差分法简洁地处理三维复杂地表附近的走时计算公式的建立问题，同时还能保证算法整体的精度和简洁性 (因为算法整体上的绝大部分计算工作由计算精度较高的迎风双线性插值法在均匀的正方体网格中实施).针对上述问题④的多波型走时计算问题，唐小平和白超英 (2009) 采用分区多步计算技术解决.此外，Rawlinson和Sambridge (2004) 采用多级次的快速推进法处理三维水平地表条件下的多波型问题，这些算法具有非常灵活的特点，但是该算法每迭代一次，仅选取波前上的一个最小走时点作为扩展点.与之相比，Kim和Folie (2000) 提出的群快速推进法，在原理上则更接近于真实模拟地震波实际波前的传播过程，其在每次波前扩展计算时不是选择一个全局最小走时点作为扩展点，而是选择一组满足波前扩展条件的扩展点.但是，群快速推进法不具备计算多波型走时计算的能力.鉴于此，本文将快速推进法中的多级次策略 (Rawlinson and Sambridge, 2004) 与群推进法相结合，提出多级次群推进法来处理多波型走时计算问题.
综上所述，为了解决三维复杂山地条件下多波型走时计算问题，本文综合多种算法的优势，提出一种多级次群推进迎风混合法.下面将首先阐述三维复杂山地模型的网格剖分问题；其次，阐述三维复杂山地条件下的迎风混合法；然后，阐述三维复杂山地条件下的多级次群推进法，包括三维复杂山地条件下的群推进法、多波型走时计算问题的描述和针对各种波型的多级次算法等问题；最后，给出算法分析与计算实例，用于检验算法的计算精度和其在三维复杂山地介质中的适应性和稳定性.
2 三维复杂山地模型的网格剖分如图 1所示，为了简洁地实现三维复杂山地条件下的地震波走时计算，采用不等距网格剖分三维复杂山地模型.该网格剖分方法的基本步骤为：①读入三维地表高程数据，并根据该高程数据确定整体的计算空间，即以地表最高点为起点向地下延伸到研究所关心的目标深度作为计算空间；②根据计算精度要求，采用合适网格间距的正方体网格均匀地剖分整个计算空间；③根据地表高程数据，去除地表以上部分的网格，同时保留地表上采样点在网格上的精确高程位置.
|
图 1 三维复杂山地模型的建立 Fig. 1 Establishment of a 3D complex mountainous model |
上述网格剖分步骤生成的网格如图 1所示，该网格的大部分区域为网格间距 (如上步骤②可知:该网格间距的大小是根据计算精度的要求来确定的) 相等的均匀的正方体网格，而在地表处为不等距网格，这些不等距网格的网格间距是通过计算地表高程点与最近的规则网格点的距离来确定的.与其他的一些复杂网格和阶梯状近似的网格相比，该网格具有如下优势：①准确定位地表上采样点的高程位置；②在地表附近无需采用其他复杂的多面体网格 (如四面体)；③地表附近的简洁处理，不会给地表附近的局部算法带来过多的复杂性；④整体算法的大部分计算工作在常规均匀的正方体网格中进行，进而能保证整体算法的简洁性.
3 三维复杂山地条件下的迎风混合法目前在三维走时计算方面，有限差分方法是最简洁的一种算法，同时也很易推广到解决三维复杂山地问题，但是该算法的计算精度相对有限.要想获得更高的计算精度只能通过采用更高阶精度的差分格式或提高网格密度来实现，但这两种途径均需要大幅度提高计算量，尤其是在处理三维空间问题时.针对该问题，孙章庆等(2015) 曾提出了迎风双线性插值的三维地震波走时计算方法，该方法能获得不错的计算精度.但是，该方法却不易处理三维复杂地形问题，因为在复杂的网格中该算法的公式建立非常复杂.综合上述因素，在此提出一种迎风混合法.该算法，在地表附近采用简洁的不等距迎风差分法，而在地表以外的其他计算空间，采用迎风双线性插值法.
如图 2a所示，采用不等距网格处理三维复杂山地模型，在进行走时计算公式建立时，主要有三种类型的网格节点：①地表上的点 (如图 2a所示的点S)；②与地表直接相邻的规则正方体网格节点 (如图 2a所示的点Sz)；③与地表不相邻的规则正方体网格节点 (如图 2a所示的点D).由于不同类型的点所处的位置不同，且周围网格节点的已知条件也不一样，因此下面分别阐述不同类型网格节点走时计算公式的建立方法.
|
图 2 迎风混合法 (a) 三维复杂地表附近的局部网格； (b) 不等距迎风差分法； (c) 迎风双线性插值法. Fig. 2 Upwind hybrid method (a) Local grid nearby the 3D complex earth′s surface; (b) Upwind finite difference method with un-equal grid spacing; (c) Upwind bilinear interpolation method. |
在公式建立之前，首先需要阐述一个走时计算过程中的关键问题：在走时计算的每一步扩展计算时，由于波的传播扩展方向不同 (尤其是在复杂介质中)，不同被算点周围走时值已知的网格节点的个数和位置通常是不同的；在这种情况下，选取周围的哪些网格节点进行被算点的走时计算就显得尤为重要而又棘手，因为其事关算法对复杂介质的适应性.迎风算法就是解决该问题的一种非常有效的算法，对复杂介质有无条件稳定性 (Sethian and Popovici, 1999).在物理上迎风算法满足地震波前传播的熵守恒定律，其最早用于处理传播曲面在复杂介质中遇到的拓扑结构变化、拐角、尖角等问题.在描述地震波波前面传播问题方面，其实际上对应着Huygens原理和Fermat原理，即地震波总是从走时值更小的区域和方向传入，而该方向即可视为局部计算的子震源方向，该方向上走时值最小的那些邻近点即可作为被算点的已知条件，也即迎着地震波波前面传来的方向去建立局部走时计算公式.下面结合迎风算法的思想，来建立上述三种不同位置网格节点的走时值的计算公式：
(1) 当计算地表上点的走时值时 (如图 2a所示的点S)
此种情况，因为被算点位于地表上，所以地震波只能从地表以下的方向传入或从邻近地表点的方向传入.当地震波由地表以下的方向传入时，如图 2a所示此时点Sz的走时值为已知，直接在点Sz处采用不等距迎风差分策略离散Eikonal方程，可知：
|
(1) |
其中，tSx=min (tSx1, tSx2)，tSy=min (tSy1, tSy2)，h为正方体网格的网格间距，d为线段SSz的长，s为当前计算网格单元的平均慢度 (速度的倒数).对方程 (1) 进行整理可知地表上点S的走时值的计算公式为：
|
(2) |
当地震波由邻近地表点传入时则直接利用Fermat原理和Huygens原理，视其他的地表点为子震源点，直接计算点S的走时值：
|
(3) |
其中，Si为被算点周围的邻近的地表点，s为当前网格单元的平均慢度.通常情况下，对于地表点，地震波通常从一个斜下方传入，因此最后利用Fermat原理将公式 (2)、(3) 统一在一起，即可获得地表上S点的走时计算公式：
|
(4) |
(2) 与地表不相邻的规则正方体网格节点 (如图 2a所示的点D)
因为与地表直接相邻的规则正方体网格节点 (如图 2a所示的点Sz) 是同时包含不等距网格和规则网格的一种较为复杂的情况，所以在此先阐述与地表不相邻的规则正方体网格节点的情况.此时在当前正方体网格处建立局部坐标系 (如图 2c所示)，然后直接采用水平地表条件下的迎风双线性插值法 (孙章庆等，2015) 即可获得被计算点D的走时值tD的计算公式：
|
(5) |
其中，k1=(tA－tO)/h，k2=(tB－tO)/h，Δt1=tB－tO，Δt2=tA－tO，tO=min (tD1i)i=1, 2, …, 6，(D1i)i=1, 2, …, 6分别表示与点D直接相邻且相距一个网格间距的6个规则正方体网格节点，tA=min (tOx+, tOx－)，tB=min (tOy+, tOy－)，tOx+, tOx－分别表示O点在x方向上相邻的两个网格节点的走时值，tOy+, tOy－分别表示O点在y方向上相邻的两个网格节点的走时值，tA、tB、tO分别表示点A、B、O的走时值且为已知条件.公式 (5) 看似简洁，但由该公式后面的参数说明可知，其隐含着Fermat原理的利用，体现了迎风思想，即地震波由走时值更小的区域传入，所以选择走时值分布更小的网格单元作为已知条件去计算当前被算点的走时值.此外，公式 (5) 并没有包括k1 > 0或k2>0的情况，因为：综合利用双线性假设 (平面波假设)和Fermat原理，此时需要计算点D的走时值，说明波是从其他方向传向点D的.此时正方形单元OACB为已知条件，即暗含波由该网格单元传向点D，又基于tA=min (tOx+, tOx－)，tB=min (tOy+, tOy－)，所以必有波传经该网格单元时先经过A、B再经过O点或同时经过O点 (此时波的传播方向为垂直向上)，所以必有tA≤tO, tB≤tO，也即仅存在k1≤0, k2≤0的情况.
(3) 与地表直接相邻的规则正方体网格节点 (如图 2a所示的点Sz)
此时，因为被算点Sz上、下方网格的形态不一样，所以采用的计算策略不一样.首先，考虑假设波从上方传向点Sz，此时同样在点Sz处直接采用不等距迎风差分策略离散Eikonal方程，即方程 (1).不同的是方程的未知变量变为了tSz，求解该方程并同时综合考虑当地震波由邻近地表点传入 (则直接利用Fermat原理和Huygens原理) 点Sz的情况，可获得由上向下方向，tSz的计算公式：
|
(6) |
其中，Si、tSi、tS、tSx、tSy代表的意思和公式 (1)—(4) 是一致的，θ=h2/d2.其次，考虑波从下方传向点Sz.此时，tSz的计算公式为类似于公式 (5) 的由迎风双线性插值获得的公式，设其为tSz－.最后，同样利用Fermat原理，综合考虑波从各个方向传向点Sz，可获得计算点Sz的走时值的计算公式：
|
(7) |
其中，tSz+为公式 (6)，tSz－为类似于公式 (5) 的计算公式.
综上所述，为了计算三维复杂山地条件下的地震波走时，在此提出了一种迎风混合法.当计算地表上的走时值时，采用简洁的不等距迎风差分法，即公式 (4)；当计算地表以外的规则网格节点处的走时值时，则采用计算精度更高的迎风双线性插值法，即公式 (5)；当计算地表与规则网格间的过渡区域的走时值时，则同时采用不等距迎风差分法和迎风双线性插值法并通过Fermat原理将它们统一在一起，即公式 (7).混合法具有如下优点：①能采用相对简洁的公式和方案，处理复杂地表问题，并确保不对地表高程位置做近似处理；②能保证整个算法的绝大部分计算工作在规则均匀的正方体网格中进行；③能保证整个算法的绝大部分计算工作，由计算精度较高的算法来执行；④由于引入了迎风算法作为控制条件，所以该算法在面对任意复杂介质时具有无条件稳定性 (Sethian and Popovici, 1999)，因此其并不会出现复杂介质中数值求解Eikonal方程时，所遇到的因奇异值而引起的计算发散的问题.除了如上的局部走时扩展计算外，完成三维复杂山地条件下的地震波走时计算，还需要另外一个重要的算法，即整体波前扩展算法，下面将阐述该问题.
4 三维复杂山地条件下的多级次群推进法 4.1 三维复杂地表条件下的群推进法群推进法实际上是快速推进法的一种更高级的波前扩展方式，与快速推进法类似，其也是通过一个波前曲面的扩展演化来模拟地震波的传播过程 (Kim and Folie, 2000).不同之处在于，其在选择波前曲面内的当前扩展点时，不是一次选择一个全局最小走时点，而是在一定限定条件下一次同时选择很多点作为扩展点.这种处理方式实际上与地震波的传播规律是更相吻合的，因为地震波前在传播过程中也是整体向前推进的.
如图 3所示，群推进法在实施时，将计算区域内的全部网格节点的属性分为三种，并用idTT的数值来表示：其中，idTT=2表示计算已经完成的点，图中为黑色充填的圆点；idTT=1表示当前波前上的点，图中为灰色充填的圆点；idTT=0表示计算还未扩展到达的区域的网格节点，图中为白色充填的圆点.其中包括所有波前点的曲面近似地表示当前波前，图中为灰色填充的窄带面.群推进法在实现时分为初始化和扩展推进两个步骤.其中，初始化主要确定震源点的位置及属性、初始波前面的组成和走时值，而扩展推进阶段主要需要处理各网格节点属性的转化、走时值的计算以及整个计算是否完成的判断等问题.下面将分别详细阐述这两个步骤：
|
图 3 复杂山地条件下的群推进法 Fig. 3 Group marching method in complex mountainous conditions |
·初始化：
① 设震源点的坐标为 (si, sj, sk)，其走时值设为T (si, sj, sk)=0.0，其属性设为idTT=2；
② 设置与震源点相连的且在计算区域范围内的网格节点 (最多为6个，分别为：(si-1, sj, sk)、(si+1, sj, sk)、(si, sj-1, sk)、(si, sj+1, sk)、(si, sj, sk-1)、(si, sj, sk+1))的属性为idTT=1，其走时值通过局部解析公式计算而得，它们共同构成了初始波前面，并取它们中的最小走时作为扩展推进时的起始值TM；
③ 设置其他所有网格节点的属性为idTT=0，其走时值设置为一个相对非常大的数 (该数值远大于最终可以算出来的可能的最大的走时值，之所以如此设置是为了实施以上阐述的迎风性质，因为其值为无穷大，所以其永远不会被选为最小走时点作为走时扩展计算的已知条件)；
④ 为了判断波前点是否能够被纳入群组，而设定
·扩展推进：
① 设定：TM=TM+ΔT；
② 将所有波前面内的走时值T (i, j, k)小于TM的网格节点(i, j, k)纳入到扩展群组中，并计算与它们直接相连的属性idTT≤1的节点的走时值，若这些节点的属性为idTT=0，则将它们的属性改为idTT=1，即将它们纳入波前曲面之内；
③ 将步骤②中被纳入扩展群组的所有网格节点 (i, j, k) 的属性由idTT=1改为idTT=2，即将它们从波前曲面中移除，进而这些点的走时计算已完成；
④ 判断波前曲面上的点是否为空，若是则计算完成，否则跳回步骤①继续扩展推进.
上述初始化的过程实际上是模拟地震波的激发状态和设置相应的计算参数，而扩展推进过程实际上是模拟地震波的传播过程，其是通过网格节点属性的不断更新变化来实现的.这个过程的效果是地震波传播经过的区域网格节点的属性最终都变为idTT=2，而当所有网格节点的属性均变为idTT=2时则计算结束.此外，为了让上述算法适应三维复杂地表条件，将地表以上的网格节点的属性也设为idTT=2.根据以上群推进法的实现过程可知，这样处理的效果是地表以上的网格节点永远不参与计算，这与地震波不会穿越地空界面传入空气中的物理事实是相吻合的.
上述复杂山地条件下改进的群推进算法，如果不加入任何控制条件，仅能计算出初至波的走时值，即最小走时值.这显然是不能满足对各种类型的地震波的走时计算.因此，下面首先阐述多波型走时计算问题.
4.2 多波型走时计算问题的描述通常情况下，常规走时算法一般计算出来的是初至波走时，即最小走时.如图 4a，在此不失一般性地以层状介质为例，若下覆地层中地震波传播速度低，则在下覆地层中仅存在透射波.但是，通常实际地层中下覆地层的地震波传播速度是递增的，即如图 4a所示的三层介质的速度递增v3 > v2 > v1，则此时会在下覆地层的界面处产生走时值更小的首波 (折射波)，如图 4a中A (界面1的折射波)、B区 (界面2的折射波) 的走时信息即为折射波的走时信息.如此一来，在地表上接收的仅为初至波信息.但是实际复杂介质中，还存在一些更有用的波型.如图 4b所示，在三维复杂介质中，地震波由震源点S激发后在地下介质的传播过程中，由于受地下复杂介质构造的影响，会产生各种类型的地震波，如图 4b所示的折射波 (S，Z1，Z2，Z3)、透射波 (S，T1，T2，T3)、反射波 (S，R1，R2，R3，R4)、多次反射波 (S，M1，M2，M3，M4，M5，M6)、绕射波 (S，D0→Dp→(D1，D2，…，Di，…，Dn))等的走时信息.
|
图 4 多波型走时计算问题的描述 (a) 初至波； (b) 多波至. Fig. 4 Description of traveltime computational problem of various seismic waves (a) First-arrival; (b) Multiple-arrival. |
常规群推进法是在Fermat原理控制下求取最小走时的.显然，将该算法直接用于计算其他类型的地震波的走时值是不可行的.与其他走时计算方法相比，针对各种类型地震波的走时算法需要解决如下核心问题：如何将各种不同类型的地震波的运动学特征融入到各种地震波的走时算法中，以保证算法的正确性.为了很好地处理多波型走时计算问题，下面将在以上阐述的三维复杂山地条件下的群推进法中引入多级次算法.
4.3 针对各种波型的多级次算法如图 5所示，为了实现不同波型的走时计算，特在上述复杂地表条件下的群推进法的基础上，提出如下的多级次算法，该算法的核心思想是以不同类型地震波产生的机理和传播特征为基础，通过分不同级次和不同的模型重新初始化过程，来实现对不同波型的走时的计算.下面将结合图 4a、4b和图 5，分别阐述常见的不同地震波型对应的多级次算法：
|
图 5 多级次算法 Fig. 5 Multistage algorithm |
折射波 (S，Z1，Z2，Z3)
实际上，常规初至波走时计算过程中就隐含着计算出了折射波走时信息，在此还将该问题单独提出来是因为：当需要单独计算某一地层的折射波走时时，必须屏蔽掉来自更深高速层的折射波，如图 4a所示地表上接收的折射波信息同时包含着界面1和界面2的折射波信息，而非仅为界面1产生的折射波信息.为了解决该问题，仅需要在初至波走时计算时仅在第一、二层介质中进行即可，具体操作为：在模型初始化时，将界面2及以下网格节点的属性均设置为idTT=2.根据群推进法的基本原理，如上操作的效果是界面2及以下网格节点在群推进法的扩展推进计算过程中均不参与计算，如此一来自然就屏蔽掉了来自第三高速层的折射波信息.
透射波 (S，T1，T2，T3)
同样如图 4a所示，实际上常规的群推进法能够计算出部分第二层介质中的透射波信息，但是由于下覆高速层折射波的干扰，第二层大部分区域的透射波走时无法被计算，为解决该问题采用逐层计算的多级次策略.如图 4b所示，具体步骤如下：(1) 仅在第一层介质中计算直达波的走时值 (同上，将界面1及以下所有网格节点的属性设为idTT=2，然后采用常规群推进法即可实现)；(2) 提取界面1上所有网格节点的直达波走时为初始条件，重新初始化整个模型 (界面1以上和界面2以下的所有网格节点的属性均设置为idTT=2)，仅在第二层中采用常规群推进法进行计算，即可获得第二层介质中和界面2上的透射波走时值；(3) 向下重复步骤 (2) 的操作即可获得第三层中的透射波走时.依次类推，逐层计算，结合逐层重新初始化，即可全部计算出下覆地层和界面上的透射波走时值.
反射波 (S，R1，R2，R3，R4)
根据反射波的产生机理和传播规律，反射波走时计算时首先需要计算入射波的走时值.实际上，反射界面上的入射波走时计算的方法即为上述界面上的透射波的走时值计算方法.所以如图 4b所示，反射波走时计算的步骤如下：(1) 按照计算透射波走时的步骤计算到达界面2的入射波走时 (S，R1，R2)；(2) 重新初始化模型，将界面2以下的所有网格节点的属性设为idTT=2，提取界面2上的入射波走时值作为初始条件；(3) 同样按照计算透射波的逐层计算策略，不同的是此时计算方向为向上，直至地表.其中步骤 (2) 重新初始化模型的依据是入射波到达界面的时刻即为反射波开始的时刻，所以界面上入射波的走时信息可以作为反射波走时计算的初始条件.
多次反射波 (S，M1，M2，M3，M4，M5，M6)
多次反射波的走时计算的步骤，与反射波走时计算的步骤类似，不同之处仅在于其需要更多的计算步骤和更多的模型初始化的步骤，同时其计算方向也是根据多次反射波的具体传播方向而定的，关于该问题在此不做具体重复阐述了.
绕射波 (S，D0→Dp→(D1，D2，…，Di，…，Dn))
与反射波走时计算的多级次策略略有不同，如图 4b所示以点绕射为例，绕射波走时计算的具体步骤如下：(1) 根据计算入射波的策略，计算到达绕射点的入射波的走时值；(2) 重新初始化模型，仅保留绕射点的入射波走时值为初始条件，将其他网格节点的属性均设为idTT=0；(3) 以绕射点为新的震源，计算由其产生的其他空间区域的绕射波走时值.实际上，上述绕射波走时计算的多级次策略是根据Huygens原理，将绕射点作为次级震源点进行计算的.
上述各种类型的地震波的走时计算，是以不同类型波的传播规律为控制条件，以群推进法为基本算法来实现的.它们的被实现使得本文的走时算法具备了计算多波型走时的能力.下面将给出算法分析与计算实例来验证算法的有效性.
5 算法分析与计算实例 5.1 计算精度与效率对比分析为了对比分析本文提出的新算法的计算精度与效率，采用如图 6a所示的三维山峰模型.模型的大小为8.0 km×8.0 km×6.0 km，地震波在介质中的传播速度为1.0 km·s-1.计算时采用的网格间距为20.0 m，震源点置于 (4.0 km，4.0 km，0.2 km) 处，并且震源点可以通过直线直接连接到地表上的任意位置 (可以获得精确的地表上的解析走时值).
|
图 6 计算精度对比分析 (a) 三维复杂地表模型； (b，c) 采用不等距迎风有限差分法计算得出的地表上和地表以下Y=4.0 km剖面上走时计算结果的相对误差分布； (d，e) 采用本文提出的迎风混合法计算得出的地表上和地表以下Y=4.0 km剖面上走时计算结果的相对误差分布. Fig. 6 Contrast analysis of computational accuracy (a) 3D complex surface model; (b，c) Relative error distribute on the earth′s surface and in the section of Y=4.0 km obtained by the upwind finite-difference method using un-equal grid spacing; (d，e) Relative error distribute on the earth′s surface and in the section of Y=4.0 km obtained by upwind hybrid method presented in this paper. |
图 6b，6c和图 6d，6e分别给出了单独采用不等距迎风差分法 (孙章庆等，2012)和采用本文算法时，地表上走时计算结果的相对误差和剖面Y=4.0 km上的走时计算结果的相对误差.其中，在地表上单独采用不等距有限差分法的最大误差为8.49%，平均相对误差为1.57%；而本文算法的最大相对误差为6.40%，平均相对误差为1.17%.对比分析两种算法的计算精度，可以发现本文算法的计算精度有了明显的提高，其无论是在复杂地表及附近，还是在地下介质空间中均能保证不错的计算精度.以上计算精度的对比分析表明：本文新提出的迎风混合算法在解决复杂地表走时计算问题时，采用简洁的不等距有限差分法处理地表及附近的走时计算问题，采用计算精度更高的迎风双线性插值法 (孙章庆等，2015) 处理地表以下均匀网格中的走时计算问题，是可行的，同时算法还能保证很好的计算精度.
如表 1所示，其给出了不同方法：快速推进迎风有限差分法 (Fast marching method with upwind finite-difference，FMM-UFD)、快速推进迎风混合差分法 (Fast marching method with upwind hybrid method，FMM-UHM)、群推进迎风有限差分法 (Group marching method with upwind finite-difference，GMM-UFD)、本文的群推进迎风混合法 (Group marching method with upwind hybrid method，GMM-UHM)，在分别采用不同网格间距 (计算量不同) 时，图 6a模型走时计算的CPU耗时对比 (计算的硬件设备参数为：Intel Core i5-3230 2.60GH CPU，4.00 G内存).分析该表格可以发现：①将计算耗时的第二列与第一列、第四列与第三列作对比，迎风混合法因为在绝大多数计算区域采用迎风双线性插值法，所以其保留延续了迎风双线性插值法快速准确收敛的特性 (孙章庆等，2015)，计算效率更高；②将计算耗时的第三列与第一列、第四列与第二列作对比，可以发现群推进法比快速推进法具有更高的计算效率.综上所述，本文提出的群推进迎风混合法具有相对更高的计算效率.
|
表 1 计算效率对比 (CPU耗时，s) Table 1 Comparison of calculation efficiency (CPU running time，s) |
为了检验本文提出的三维复杂山地多级次群推进迎风混合法在计算多波型走时时的有效性，采用如图 7所示的一个由复杂地表和起伏界面构造的三层层状介质模型，大小为8.0 km×8.0 km×6.0 km，从上至下三层介质中地震波的传播速度分别为1.0 km·s-1、1.5 km·s-1、2.0 km·s-1；计算时采用的网格间距为20.0 m，震源点位于地表上其坐标为 (0.0 km，0.0 km，1.0 km).
|
图 7 三维复杂地表起伏界面模型 Fig. 7 3D complex surface topographic model with rugged interfaces |
如图 8a，8b所示，首先采用常规的非多级次群推进法计算整个模型中的初至波走时.从初至波走时中，可以发现明显的来自于下覆高速层的折射波走时信息.如果需要计算折射波的走时信息，则采用常规计算初至波走时的计算策略即可.但是如果需要计算透射波和对应于界面反射波的入射波的走时信息，该计算策略则会失效.因为如图 8a，8b所示，在第一层和第二层介质中，由于下覆高速层的存在，直达波、透射波的走时信息被走时值更小的折射波 (首波) 所取代，所以无法计算到界面1和界面2上完整的直达波 (入射波和透射波) 的走时信息.因此，需要屏蔽掉来自下覆地层的折射波信息，在此采用逐层向下计算透射波走时的多级次策略.图 8c，8d分别给出了X-Y=0和Y=0剖面上直达波的走时分布，其中包括了分别通过界面1和界面2的透射波的走时信息.对比图 8c，8d和图 8a，8b可以发现，逐层向下计算的多级次策略成功地屏蔽掉了来自下覆高速层的折射波信息，能够获得界面1和界面2上完整的入射波走时信息，其是计算反射波走时的基础.
|
图 8 初至波、折射波、直达波、透射波的走时分布 (a) X-Y=0剖面上初至波走时分布，其中包括折射波信息； (b) Y=0剖面上初至波走时分布，其中包括折射波信息； (c) X-Y=0剖面上直达波走时分布，其中包括经过界面进入下覆地层的透射波信息； (d) Y=0剖面上直达波走时分布，其中包括经过界面进入 下覆地层的透射波信息. Fig. 8 Traveltime distributions of first-arrival, refracted, direct-arrival, and transmitted waves (a) Traveltime contours of first-arrival waves superimposed on the section X-Y=0, and refracted wave included; (b) Traveltime contours of first-arrival waves superimposed on the section Y=0, and refracted wave included; (c) Traveltime contours of direct-arrival waves superimposed on the section X-Y=0, and transmitted waves included; (d) Traveltime contours of direct-arrival waves superimposed on the section Y=0 and transmitted waves included. |
由以上的阐述可知，为了计算某界面上反射波的走时值，首先需要计算界面上准确的入射波的走时值，采用逐层向下的多级次策略.图 9a和图 10a分别给出了界面1和界面2上的直达波的走时值，及入射波的走时值，而以此为基础即可计算反射波的走时值.图 9b给出了界面1的反射波的走时分布.同理，可以采用如上的多级次策略，计算多次波的走时值，只不过是多次反射波走时计算时采用的级次和步骤更多.图 10b—10d分别给出了第二层介质中多次反射波的走时分布.最后采用绕射波多级次策略即可获得绕射波走时值.图 9c给出了绕射波的走时分布.
|
图 9 X-Y=0剖面上反射波与绕射波的走时分布 (a) 到达界面1的入射波的走时分布； (b) 界面1上的反射波的走时分布； (c) 界面1某一绕射点的绕射波的走时分布. Fig. 9 Traveltime distributions of reflected and diffracted waves in section X-Y=0 (a) Traveltime contours of incident waves at interface 1; (b) Traveltime contours of reflective waves at interface 1; (c) Traveltime contours of diffracted waves at interface 1. |
|
图 10 X-Y=0剖面上多次反射波的走时分布 (a) 到达界面2的入射波的走时分布； (b) 界面2到界面1的反射波走时分布；(c) 界面1到界面2的反射波走时分布； (d) 界面2到地表上的反射波走时分布. Fig. 10 Traveltime distributions of multiple waves in section X-Y=0 (a) Traveltime contours of incident waves at interface 2; (b) Traveltime contours of reflected waves from interface 2 to interface 1; (c) Traveltime contours of reflected wave from interface 1 to interface 2; (d) Traveltime contours of reflected wave from interface 2 to the surface. |
上述计算实例充分说明了本文算法能在很好适应复杂地表条件和起伏界面的前提下，灵活地计算各种地震波型的走时值.
6 结论与讨论为了计算三维复杂地表条件下的多波型走时，提出了一种多级次群推进迎风混合法.该方法采用不等距网格精确定位地表位置，采用迎风混合法计算三维复杂地表条件下的地震波走时，采用群推进法作为算法的整体波前扩展方式，采用多级次策略计算各种地震波型的走时值.基于以上策略，新提出的算法具备如下特点：①不对三维复杂地表位置做任何近似处理，而是采用不等距网格对地表位置进行简洁且准确地定位；②在地表附近采用简洁的不等距迎风有限差分法处理走时计算问题，而在其他更多的均匀网格区域采用计算精度更高的迎风双线性插值法，所以总体上算法简洁且精确；③通过在三维复杂地表条件下的群推进法中引入多级次策略，进而有针对性地提出了三维复杂地表条件下灵活计算各种地震波型的走时算法.针对算法的精度对比分析和计算实例分析，可以得出如下结论：①新算法具有很好的计算精度，其无论是最大误差，还是平均相对误差都低于单独采用不等距迎风有限差分法，同时其无论是处理地表及附近区域还是其他地下区域均具有很好的计算精度；②新算法能稳定、灵活地适应三维复杂地表问题；③因为综合采用了收敛速度更快的迎风双线性插值的局部走时算法和迭代推进效率更高的群推进法，所以新算法具有相对更高的计算效率；④新算法能灵活稳定地计算各种波型的走时值.通过三维复杂地表条件下多波型走时算法的研究和实例分析，发现三维复杂地表条件下地震波的传播特性非常复杂，同时存在很多波型，因此采用本文算法对这些波型的传播特征进行分析，有利于更加深入地分析地震波在三维复杂地表条件下的传播规律，进而为三维复杂地表区域地震数据采集观测系统的设计和三维复杂地表条件下地震数据处理提供重要的参考信息和正演模拟工具.
| Asakawa E, Kawanaka T. 1993. Seismic ray tracing using linear traveltime interpolation. Geophysical Prospecting, 41(1): 99-111. DOI:10.1111/gpr.1993.41.issue-1 | |
| Bai C Y, Greenhalgh S, Zhou B. 2007. 3D ray tracing using a modified shortest-path method. Geophysics, 72(4): T27-T36. DOI:10.1190/1.2732549 | |
| Bai C Y, Huang G J, Li Z S. 2011. Simultaneous inversion combining multiple-phase travel times within 3D complex layered media. Chinese J. Geophys. (in Chinese), 54(1): 182-192. DOI:10.3969/j.issn.0001-5733.2011.01.019 | |
| Huang G J, Bai C Y. 2013. Simultaneous inversion of three model parameters with multiple classes of arrival times. Chinese J. Geophys. (in Chinese), 56(12): 4215-4225. DOI:10.6038/cjg20131224 | |
| Huang G J, Bai C Y, Qian W. 2015. Simultaneous inversion of three model parameters with multiple phases of arrival times in spherical coordinates. Chinese J. Geophys. (in Chinese), 58(10): 3627-3638. DOI:10.6038/cjg20151016 | |
| Jing X L, Yang C C, Wang S Q. 2007. A improved seismic reflection tomographic method. Chinese J. Geophys. (in Chinese), 50(6): 1831-1836. | |
| Julian B R, Gubbins D. 1977. Three-dimensional seismic ray tracing. J. Geophys., 43: 95-113. | |
| Kim S, Folie D. 2000. The group marching method: An O (N) level set eikonal solver.//70th Annual International Meeting, SEG, Expanded Abstracts, 2297-2300. | |
| Lan H Q, Zhang Z J. 2013a. A high-order fast-sweeping scheme for calculating first-arrival travel times with an irregular surface. Bulletin of the Seismological Society of America, 103(3): 2070-2082. DOI:10.1785/0120120199 | |
| Lan H Q, Zhang Z J. 2013b. Topography-dependent eikonal equation and its solver for calculating first-arrival traveltimes with an irregular surface. Geophys. J. Int., 193(2): 1010-1026. DOI:10.1093/gji/ggt036 | |
| Lan H Q, Chen J Y, Zhang Z J. 2014. A fast sweeping scheme for calculating P wave first-arrival travel times in transversely isotropic media with an irregular surface. Pure Appl. Geophys., 171(9): 2199-2208. DOI:10.1007/s00024-014-0836-5 | |
| Li P M, Mei S Q, Ma Q B. 2013. An improved bilinear interpolation travel-time ray-tracing method. Oil Geophysical Prospecting (in Chinese), 48(4): 553-558. | |
| Liu F, Zhang D, Yang Y, et al. 2012. A fast numerical solution to minimum equation in 3-D seismic LTI ray tracing. Journal of Wuhan University (Natural Science Edition) (in Chinese), 58(5): 395-400. | |
| Mei S Q, Deng F, Zhong B S, et al. 2010. The 3D ray tracing method base on the improved bilinear traveltime interpolation. Computing Techniques for Geophysical and Geochemical Exploration (in Chinese), 32(2): 152-157. | |
| Nakanishi I, Yamaguchi K. 1986. A numerical experiment on nonlinear image reconstruction from first-arrival times for two-dimensional island arc structure. Journal of Physics of the Earth, 34(2): 195-201. DOI:10.4294/jpe1952.34.195 | |
| Qu C, Zhou H L, Zhao D P. 2007. Deep structure beneath the west margin of Philippine Sea Plate and South China Sea from P and S wave travel time tomography. Chinese J. Geophys. (in Chinese), 50(6): 1757-1768. | |
| Rawlinson N, Sambridge M. 2004. Multiple reflection and transmission phases in complex layered media using a multistage fast marching method. Geophysics, 69(5): 1338-1350. DOI:10.1190/1.1801950 | |
| Sethian J A, Popovici A M. 1999. 3-D traveltime computation using the fast marching method. Geophysics, 64(2): 516-523. DOI:10.1190/1.1444558 | |
| Sun J G. 2000. Limited-aperture migration. Geophysics, 65(2): 584-595. DOI:10.1190/1.1444754 | |
| Sun Z Q, Sun J G, Han F X. 2012. Traveltime computation using the upwind finite difference method with nonuniform grid spacing in a 3D undulating surface condition. Chinese J. Geophys. (in Chinese), 55(7): 2441-2449. DOI:10.6038/j.issn.0001-5733.2012.07.028 | |
| Sun Z Q, Sun J G, Yue Y B, et al. 2015. 3D traveltime computation using fast marching upwind bilinear interpolation method. Chinese J. Geophys. (in Chinese), 58(6): 2011-2023. DOI:10.6038/cjg20150616 | |
| Tang X P, Bai C Y. 2009. Multiple ray tracing within 3-D layered media with the shortest path method. Chinese J. Geophys. (in Chinese), 52(10): 2635-2643. DOI:10.3969/j.issn.0001-5733.2009.10.024 | |
| Vidale J. 1988. Finite-difference calculation of travel times. Bulletin of the Seismological Society of America, 78(6): 2062-2076. | |
| Vidale J E. 1990. Finite-difference calculation of traveltimes in three dimensions. Geophysics, 55(5): 521-526. DOI:10.1190/1.1442863 | |
| Vinje V, Iversen E, Gjøystdal H. 1993. Traveltime and amplitude estimation using wavefront construction. Geophysics, 58(8): 1157-1166. DOI:10.1190/1.1443499 | |
| Wang H, Chang X. 2000. 3-D ray tracing method based on graphic structure. Chinese J. Geophys. (in Chinese), 43(4): 534-541. | |
| Wang H Z, Ma Z T. 1999. Traveltime calculation in 3D media with arbitrary velocity distribution.//69th Annual International Meeting, SEG, Expanded Abstracts, 1263-1266. | |
| Zhang D, Fu X R, Yang Y, et al. 2009. 3-D seismic ray tracing algorithm based on LTI and partition of grid interface. Chinese J. Geophys. (in Chinese), 52(9): 2370-2376. DOI:10.3969/j.issn.0001-5733.2009.09.023 | |
| Zhang D, Zhang T T, Qiao Y F, et al. 2013. A 3-D ray tracing method based on B-spline traveltime interpolation. Oil Geophysical Prospecting (in Chinese), 48(4): 559-566. | |
| 白超英, 黄国娇, 李忠生. 2011. 三维复杂层状介质中多震相走时联合反演成像. 地球物理学报, 54(1): 182–192. DOI:10.3969/j.issn.0001-5733.2011.01.019 | |
| 黄国娇, 白超英. 2013. 多震相走时联合三参数同时反演成像. 地球物理学报, 56(12): 4215–4225. DOI:10.6038/cjg20131224 | |
| 黄国娇, 白超英, 钱卫. 2015. 球坐标系下多震相走时三参数同时反演成像. 地球物理学报, 58(10): 3627–3638. DOI:10.6038/cjg20151016 | |
| 井西利, 杨长春, 王世清. 2007. 一种改进的地震反射层析成像方法. 地球物理学报, 50(6): 1831–1836. | |
| 李培明, 梅胜全, 马青坡. 2013. 一种改进的双线性插值射线追踪方法. 石油地球物理勘探, 48(4): 553–558. | |
| 刘锋, 张东, 杨艳, 等. 2012. 三维LTI射线追踪极小值方程的快速数值解法. 武汉大学学报 (理学版), 58(5): 395–400. | |
| 梅胜全, 邓飞, 钟本善, 等. 2010. 基于改进的双线性旅行时插值的三维射线追踪. 物探化探计算技术, 32(2): 152–157. | |
| 瞿辰, 周蕙兰, 赵大鹏. 2007. 使用纵波和横波走时层析成像研究菲律宾海板块西边缘带和南海地区的深部结构. 地球物理学报, 50(6): 1757–1768. | |
| 孙章庆, 孙建国, 韩复兴. 2012. 三维起伏地表条件下地震波走时计算的不等距迎风差分法. 地球物理学报, 55(7): 2441–2449. DOI:10.6038/j.issn.0001-5733.2012.07.028 | |
| 孙章庆, 孙建国, 岳玉波, 等. 2015. 基于快速推进迎风双线性插值法的三维地震波走时计算. 地球物理学报, 58(6): 2011–2023. DOI:10.6038/cjg20150616 | |
| 唐小平, 白超英. 2009. 最短路径算法下三维层状介质中多次波追踪. 地球物理学报, 52(10): 2635–2643. DOI:10.3969/j.issn.0001-5733.2009.10.024 | |
| 王辉, 常旭. 2000. 基于图形结构的三维射线追踪方法. 地球物理学报, 43(4): 534–541. | |
| 张东, 傅相如, 杨艳, 等. 2009. 基于LTI和网格界面剖分的三维地震射线追踪算法. 地球物理学报, 52(9): 2370–2376. DOI:10.3969/j.issn.0001-5733.2009.09.023 | |
| 张东, 张婷婷, 乔友锋, 等. 2013. 三维旅行时场B样条插值射线追踪方法. 石油地球物理勘探, 48(4): 559–566. | |