地震走时层析成像是近地表速度建模的重要方法, 而快速、准确且能适应复杂模型的射线追踪方法是其中的关键技术之一^[1-2]。一方面, 射线追踪被用来合成理论走时, 另一方面, 追踪得到的射线路径可进一步用于构建Fréchet层析核函数。传统的三维射线追踪有试射(Shooting)法和弯曲(Bending)法^[3-4], 但这两种算法在复杂模型中常遇到阴影区与焦散面等问题, 并且由于每次只能处理一对炮检点, 计算精度与效率均比较低。NAKANISHI等^[5]和MOSER^[6]将图论和DIJKSTRA^[7]算法引入走时计算而提出了最短路径法^[8], KLIMES等^[9]将该算法扩展到求解三维模型。该类算法的局部走时计算公式简洁, 算法稳定性高, 然而在复杂构造模型中容易陷入局部极值, 并且为保证计算精度需要较大的计算开销, 因此不适用于大规模三维复杂模型。

射线追踪的另一种实现方式是先求程函方程, 然后从接收点出发沿走时梯度方向反向追踪到炮点得到射线路径^[10]。VIDALE^[11]首先提出了采用有限差分法求解程函方程, 但其方法采用盒式扩展来追踪波前, 不满足波前传播的物理规律, 并且在强速度界面处可能由于负数开平方而出现不稳定现象^[12]。为了解决这些问题, 许多学者对传统的有限差分算法进行了大量的研究和改进^[13-14]。SETHIAN^[15]和POPOVICI等^[16]首次将快速推进算法(Fast marching method, FMM)引入地震走时计算, 并得到了广泛的应用^[17-18]。该方法利用满足波前演化“熵条件”的迎风差分格式求解程函方程, 采用窄带技术模拟波前扩展, 利用堆排序算法保存窄带中节点的走时, 被公认为目前精度最高、效率最高、且对任意复杂模型均无条件稳定的射线追踪算法之一^[19-20]。

但是, 由于标准FMM算法只将坐标轴方向上相邻网格点引入局部差分计算, 因而在对角线方向上存在较大的计算误差。即使引入高阶迎风差分, 其误差分布仍然存在明显的对角线方向特征, 且高阶差分本身也需要消耗更多的计算成本。针对这一问题, 多位学者对FMM算法在计算精度上进行了改进^[21-23]。DANIELSSON等^[24]提出了以计算网格点全部邻点(包含对角线方向邻点)建立混合差分模板的方法, 有效改善了FMM算法的精度, 但是该方法难以推广到三维。HASSOUNA等^[25]提出了多模板快速推进法(Multi-stencils fast marching method, MSFM), 其原理是通过坐标变换引入2个差分模板分别计算坐标轴方向以及对角线方向上的网格点, 并选取其中满足迎风条件的解, 该方法显著提高了标准FMM的计算精度且容易扩展到高维。

在完成走时计算的基础上, 射线路径可通过沿走时梯度逐步反向追踪获得^[26-28]。王飞等^[29]采用角度定义走时梯度方向, 并分别讨论了不同入射方向和入射角度下射线追踪的过程。该算法计算效率较高, 但计算精度受到网格点处走时梯度精度的制约, 并且根据入射方向和角度来判断射线在网格穿出位置的算法仅适用于二维, 在高维情况下将变得异常繁琐。

本文提出了一种新的适用于复杂近地表走时层析反演的射线追踪算法。首先在二维MSFM算法的基础上将其扩展到三维, 并给出了起伏地表条件下的波前扩展策略。随后从接收点出发, 采用三线性插值和Runge-Kutta方程求取走时梯度ΔT, 并沿梯度方向逐步反向追踪获得射线路径。文中分别以均匀模型和复杂近地表模型作为测试模型进行了射线追踪试算。

1 方法原理和实现 1.1 三维FMM算法

MSFM是FMM的改进算法, 首先对FMM算法进行简单介绍。

在各向同性介质中, 程函方程(1)描述了波前传播时间T与慢度S的关系:

$\nabla T\left( {x,y,z} \right) \cdot \nabla T\left( {x,y,z} \right) = {S^2}\left( {x,y,z} \right)$

(1)

走时场T并非在空间所有点都可微, 在复杂介质中常会出现波前自交现象, 即初至波前和后续波前扭结成一个不连续点。解决此问题的方法为引入“熵满足”条件来平滑不连续性, 如迎风差分公式:

$\begin{array}{l} \max {\left( {D_a^{ - x}T_{i,j}^k, - D_b^{ + x}T_{i,j}^k,0} \right)^2} + \\ \max {\left( {D_c^{ - y}T_{i,j}^k, - D_d^{ + y}T_{i,j}^k,0} \right)^2} + \\ \max {\left( {D_e^{ - z}T_{i,j}^k, - D_f^{ + z}T_{i,j}^k,0} \right)^2} = {\left( {S_{i,j}^k} \right)^2} \end{array}$

(2)

式中:T_{i, j}^k表示坐标(i, j, k)处的走时; D^-和D⁺分别为前向和后向差分算子; a, b, c, d, e, f代表迎风差分算子的阶数。RAWLINSON等^[28]详细给出了不同阶数下的差分格式, 以二阶差分为例, 有:

$\begin{array}{l} D_2^{ - x}T_{i,j}^k = \frac{{3T_{i,j}^k - 4T_{i - 1,j}^k + T_{i - 2,j}^k}}{{2\delta x}}\\ D_2^{ - y}T_{i,j}^k = \frac{{3T_{i,j}^k - 4T_{i,j - 1}^k + T_{i,j - 2}^k}}{{2\delta y}}\\ D_2^{ - z}T_{i,j}^k = \frac{{3T_{i,j}^k - 4T_{i,j}^{k - 1} + T_{i,j}^{k - 2}}}{{2\delta z}} \end{array}$

(3)

式中:δx, δy, δz分别为x, y, z方向上的网格间距。

将(3)式代入(2)式, 整理可得:

$\begin{array}{*{20}{l}} {\max {{\left[ {\frac{3}{{2\delta x}}\left( {T_{i,j}^k - {T_1}} \right),0} \right]}^2} + }\\ {\max {{\left[ {\frac{3}{{2\delta y}}\left( {T_{i,j}^k - {T_2}} \right),0} \right]}^2} + }\\ {\max {{\left[ {\frac{3}{{2\delta z}}\left( {T_{i,j}^k - {T_3}} \right),0} \right]}^2} = {{\left( {S_{i,j}^k} \right)}^2}} \end{array}$

(4)

式中:T₁, T₂, T₃分别表示计算网格点三个坐标轴方向邻点走时中的较小值。FMM采用类似于DIJKSTRA^[7]算法的窄带技术模拟波前扩展过程。

1.2 三维MSFM算法

三维情况下, 标准FMM算法仅利用了计算网格点周围坐标轴方向上的8个相邻网格点, 因而在对角线方向上存在较大计算误差。MSFM的基本思想是利用坐标变换产生多个差分模板, 将对角线方向邻点引入计算, 从而提高迎风差分法的计算精度^[25]。

如图 1所示, 模板S_w与网格点相交于l₁, l₂, p₁, p₂, q₁和q₂, 其中r₁=[r₁₁ r₁₂ r₁₃]^T, r₂=[r₂₁ r₂₂ r₂₃]^T, r₃=[r₃₁ r₃₂ r₃₃]^T分别为沿l₁l₂, p₁p₂, q₁q₂方向的单位向量, 令U₁, U₂, U₃分别为沿r₁, r₂, r₃的方向导数, 则有如下关系:

$\begin{array}{l} {U_i} = {\mathit{\boldsymbol{r}}_i} \cdot \nabla T\left( {x,y,z} \right) = {r_{i1}}{T_x} + {r_{i2}}{T_y} + {r_{i3}}{T_z}\\ \;\;\;\;\;\;i = 1,2,3 \end{array}$

(5)

式中:T_x, T_y, T_z分别表示T对x, y, z的偏导数。令U=[U₁ U₂ U₃], (5)式可写成:

$\mathit{\boldsymbol{U}} = \mathit{\boldsymbol{R}}\nabla T\left( {x,y,z} \right)$

(6)

将(6)式代入程函方程(1), 整理可得:

$\begin{array}{l} {\left| {\nabla \mathit{\boldsymbol{T}}\left( {x,y,z} \right)} \right|^2} = \nabla T{\left( {x,y,z} \right)^T}\nabla \mathit{\boldsymbol{T}}\left( {x,y,z} \right)\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; = {\mathit{\boldsymbol{U}}^T}{\left( {\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{R}}^T}} \right)^{ - 1}}\mathit{\boldsymbol{U}} = {S^2}\left( {x,y,z} \right) \end{array}$

(7)

令α为r₁和r₂间的夹角, β为r₂和r₃间的夹角, γ为r₁和r₃间的夹角, 公式(7)可简化为:

$\begin{array}{l} U_1^2 + U_2^2 + U_3^2 + 2{U_1}{U_2}\cos \alpha + \\ 2{U_1}{U_3}\cos \gamma + 2{U_2}{U_3}\cos \beta = {S^2}\left( {x,y,z} \right) \end{array}$

(8)

图 2为三维MSFM算法采用的6个迎风差分模板S_w(w∈[1, 6]), 其中S₁包含了坐标轴方向上计算网格点的邻点, 对角线方向邻点则包含在其它5个模板中。以模板S₂为例, 采用二阶差分近似U₁, U₂, U₃, 则有:

$\begin{array}{l} {U_\upsilon } = \max \left[ {\frac{3}{{2\left\| {\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{x}}_\upsilon }} \right\|}}\left( {T_{i,j}^k - {T_\upsilon }} \right),0} \right]\\ \upsilon = 1,2,3 \end{array}$

(9)

式中:x_v为对角线方向邻点中较小走时的网格点位置; T_v为该位置对应的走时。对于模板S₂, 有:

$\begin{array}{l} {T_1} = \min \left( {\frac{{4T_{i - 1,j}^k - T_{i - 2,j}^k}}{3},\frac{{4T_{i + 1,j}^k - T_{i + 2,j}^k}}{3}} \right)\\ {T_2} = \min \left( {\frac{{4T_{i,j + 1}^{z - 1} - T_{i,j + 2}^{z - 2}}}{3},\frac{{4T_{i,j + 1}^{z - 1} - T_{i,j - 2}^{z + 2}}}{3}} \right)\\ {T_3} = \min \left( {\frac{{4T_{i,j + 1}^{z + 1} - T_{i,j + 2}^{z + 2}}}{3},\frac{{4T_{i,j - 1}^{z - 1} - T_{i,j - 2}^{z - 2}}}{3}} \right) \end{array}$

(10)

将(9)式代入(8)式求取T_{i, j}^k值, 可得2个解。由于波前传播总是从较小走时网格点到较大走时网格点, 因而需要检查计算的T_{i, j}^k是否大于邻近的T_v。只有T_{i, j}^k < T_v的解才是可接受的, 否则应采用如下公式:

$\begin{array}{l} \min \left[ {{T_\upsilon } + \frac{{\left\| {\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{x}}_\upsilon }} \right\|}}{{\upsilon \left( \mathit{\boldsymbol{x}} \right)}}} \right]\\ \upsilon = 1,2,3 \end{array}$

(11)

式中:v(x)为计算网格点处的速度值。通过联合6个模板求解最小走时T_{i, j}^k, MSFM算法将计算网格点周围所有26个邻点引入了差分计算。此外, MSFM保持了FMM的算法结构, 因而可以在FMM算法的框架下实现。

1.3 起伏地表MSFM实现流程

上述分析建立了MSFM算法的局部走时计算公式, 实际上若想完成整个计算区域波前传播时间的计算, 还需要建立算法的整体实现策略, 即波前扩展流程。此外考虑到层析成像针对于解决起伏地表条件下的复杂速度模型建模, 本文重点研究了起伏地表下MSFM算法的实现流程。

首先将所有网格点分成4个子集。第一类位于起伏地表之上, 不参与波前扩展计算; 第二类位于波前之后, 走时已经确定, 称为完成点; 第三类位于波前面上, 走时已计算出来但还需要更新, 称为窄带点; 第四类位于波前面之前, 走时还未计算, 称为远离点。整个波前扩展的实现流程如下:

1)首先将激发点处走时初始化为0, 并将其标记为完成点; 根据MSFM差分模板计算激发点周围起伏地表以下所有邻点的走时, 并将其标记为窄带点。

2)在窄带子集中选取走时最小的网格点, 并将该网格点标记为完成点。

3)判断步骤2)中选取出的最小走时点的邻点, 若邻点为完成点或位于起伏地表以上, 则保持该点原有属性不变; 若邻点为窄带点, 则采用MSFM差分模板计算走时，根据走时大小选择保持或更新原走时值; 若邻点为远离点, 则将其标记为窄带点。

4)循环步骤2)和步骤3)直到所有网格点被标记为完成点。

1.4 三维射线追踪的实现

在采用三维MSFM算法完成整个计算区域的地震走时计算后, 射线追踪可以被描述为一个初值问题, 即从接收点沿走时梯度方向出射射线, 然后不断地根据梯度调整反向追踪方向直至射线到达激发点^[11]。令r_n表示当前追踪点的位置, h表示追踪步长, 则下一追踪点的位置r_n+1可以表示为:

${\mathit{\boldsymbol{r}}_{n + 1}} = {\mathit{\boldsymbol{r}}_n} + h \cdot \nabla T$

(12)

因此, 走时梯度的准确求取成为射线追踪的关键, 下面介绍如何根据走时场求取走时梯度。首先, 网格点处走时梯度向量在三个坐标轴方向的分量可以用迎风差分公式(2)来近似, 差分阶数的选取应与MSFM算法保持一致。而对于计算区域中的非网格点位置, 走时梯度可根据网格点处的梯度值采用三线性插算法近似。以图 3为例, 模型中任意点(x, y, z)处的梯度ΔT(x, y, z)可以近似为:

$\begin{array}{l} \nabla T\left( {x,y,z} \right) = \nabla {T_1}\left( {\Delta x - x} \right)\left( {\Delta y - y} \right)\left( {\Delta z - z} \right) + \\ \nabla {T_2}\left( {\Delta x - x} \right)\left( {\Delta y - y} \right)z + \nabla {T_3}\left( {\Delta x - x} \right)y\left( {\Delta z} \right. - \\ \left. z \right) + \nabla {T_4}\left( {\Delta x - x} \right)yz + \nabla {T_5}x\left( {\Delta y - y} \right)\left( {\Delta z - z} \right) + \\ \nabla {T_6}x\left( {\Delta y - y} \right)z + \nabla {T_7}xy\left( {\Delta z - z} \right) + \nabla {T_8}xyz \end{array}$

(13)

式中:ΔT₁~ΔT₈分别为计算点周围8个相邻网格点处的梯度值; Δx, Δy, Δz分别为3个坐标轴方向的网格间隔。为了提高走时梯度的计算精度, 采用四阶Runge-Kutta方法, 将多个预估点处走时梯度的加权平均作为r_n处的走时梯度, 则公式(12)可改进为:

${\mathit{\boldsymbol{r}}_{\mathit{\boldsymbol{n}} + 1}} = {\mathit{\boldsymbol{r}}_n} + \frac{h}{6}\left( {{\mathit{\boldsymbol{k}}_1} + 2{\mathit{\boldsymbol{k}}_2} + 2{\mathit{\boldsymbol{k}}_3} + {\mathit{\boldsymbol{k}}_4}} \right)$

(14)

其中:

$\begin{array}{l} {\mathit{\boldsymbol{k}}_1} = \nabla T\left( {{\mathit{\boldsymbol{r}}_n}} \right)\;\;\;{\mathit{\boldsymbol{k}}_2} = \nabla T\left( {{\mathit{\boldsymbol{r}}_n} + \frac{h}{2}{\mathit{\boldsymbol{k}}_1}} \right)\\ {\mathit{\boldsymbol{k}}_3} = \nabla T\left( {{r_n} + \frac{h}{2}{\mathit{\boldsymbol{k}}_2}} \right)\;\;\;\;\;{\mathit{\boldsymbol{k}}_4} = \nabla T\left( {{\mathit{\boldsymbol{r}}_n} + h{\mathit{\boldsymbol{k}}_3}} \right) \end{array}$

(15)

当r_n与激发点的距离小于一定误差范围时射线追踪结束。相比较于王飞等^[28]采用的最速下降法, 本文方法避免了三维情况下繁琐的出射位置确定, 此外通过三线性插值和Runge-Kutta方法近似走时梯度, 具有更高的计算效率和稳健性。该方法还可以通过调节步长h来平衡计算效率与精度。

2 数值试算

为了从计算精度和效率两方面验证本文算法, 我们设计了2个具有代表性的模型进行验证。

2.1 模型一

模型一为v=500m/s的三维均匀模型, 大小为1000m×1000m×1000m, 激发点位于(0, 0, 0)处。选用均匀模型的原因是其解析解求取方便, 便于对算法的精度进行分析测试。采用计算走时与解析解误差的L₂范数和L_∞范数作为精度指标:

$\begin{array}{l} {L_2} = \frac{1}{n}\sum\limits_n {{{\left| {T - {T_a}} \right|}^2}} \\ {L_\infty } = \max \left( {\left| {T - {T_a}} \right|} \right) \end{array}$

(16)

式中:T为计算值; T_a为解析解; n为总的网格节点数。

本文对比了网格间距分别为20, 10, 5m时FMM算法和MSFM算法的计算精度和计算效率。表 1分别给出了模型一在3种网格间距下2种算法对应的L₂范数和L_∞范数以及CPU运行时间。分析表 1中L₂范数和L_∞范数的变化规律可知, 2种算法的计算精度均受网格间距影响, 间距越小计算精度越高; MSFM算法的平均误差和最大误差均小于FMM算法, 即使在20m网格间距下, MSFM算法的精度也要比FMM算法的精度高出约一个数量级。考虑到MSFM主要是在对角线方向上对FMM算法进行了优化, 如图 4所示, 我们提取三维模型(20m网格间距)对角线方向的剖面(图 4红线所示)进行了误差分析。图 5a和图 5b分别给出了FMM与MSFM算法的误差对比, 分析表明, MSFM算法明显减小了对角线方向的误差, 此外由于走时误差的累计效应, 对角线方向的改善也显著提升了MSFM算法的整体计算精度。图 6分别为FMM算法和MSFM算法计算的波前面与射线路径, 激发点位于模型左上角, 7个接收点均匀布设在右边界。可以明显看出, 根据MSFM算法追踪得到的射线路径几乎与真实路径完全重合, 具有更高的计算精度。

表 1同时列出了不同网格间距下2种算法的CPU运行时间, 分析表中数据可以发现20m网格间距下MSFM算法要比FMM算法增加了约2倍的计算时间, 不过随着模型网格数增加两者的差异逐渐减小。从理论角度分析, MSFM算法和FMM算法具有相似的算法结构, 其计算消耗包括两部分:一部分来自计算迎风差分格式; 另一部分来自窄带技术中的排序操作。在差分格式计算上, MSFM算法消耗的时间应大于FMM算法, 因为MSFM算法需要计算6个差分模板, 而常规的FMM算法则只需要计算1个。不过无论是MSFM算法还是FMM算法在波前扩展上都采用了DIJSTRA^[7]算法, 通过堆排序来实现窄带的添加和删除操作, 算法时间复杂度为O(NlgN)。因而模型规模越大, 用于维持窄带队列结构的删除和插入操作所消耗的时间占总计算时间的比重越高, MSFM算法和FMM算法的计算效率差异越小。三维模型通常为大规模的复杂模型, 因而采用本文提出的三维MSFM算法可以在较少增加计算成本的基础上大大提高射线追踪的精度。

2.2 模型二

为了验证本文算法在起伏地表复杂模型下的计算效果, 我们引入一个由二维复杂近地表模型生成的三维模型(图 7a), 该模型近地表速度结构非常复杂且地表起伏, 能够较好地验证算法的计算精度。模型的尺寸23.8km×20.0km×4.0km, 激发点位于地表(18km, 10km)处, 接收点均匀分布于地表。虽然复杂速度模型中无法求得完全精确的解析解, 但考虑到有限差分类算法的计算精度随网格的加密而增加, 因此可以通过对比2种方法在不同网格间距下走时等值线以及射线路径的差异来验证本文算法的计算精度。

分别采用50m×50m×20m(大网格距)与5m×5m×2m(小网格距)对模型进行离散化, 并用2种算法对模型进行测试。图 7b为小网格距下采用MSFM算法计算的走时场, 从图中可以看出地震波前面的传播过程。图 8给出了不同网格间距下激发点所在的Y=10km剖面上FMM和MSFM算法计算的走时等值线, 其中红线和蓝线分别为小网格和大网格下的计算结果。总体而言, 2种算法都清晰反映了波前面在低速层“内凸”、在高速层“外凸”的传播规律。对比图 8a和图 8b可见, 不同网格距下MSFM算法计算所得的走时差小于FMM算法, 说明MSFM算法在大网格间距下能获得更高的计算精度。

图 9为不网格间距下FMM和MSFM算法的射线追踪结果, 绿线和蓝线分别为小网格和大网格下的计算结果。为了更好的对比两种算法计算结果的差异, 我们将视点放置在了2个不同位置对三维数据体进行显示。可以明显看到, 射线追踪能够较好地获取回转波、折射波等不同类型的初至波路径, 对复杂介质有良好的适应性, 证实了本文提出的射线追踪算法的有效性。对比图 9a和图 9b以及图 9c和图 9d可见, 不同网格距下由MSFM算法得到的两组射线路径差别更小, 特别对于模型中的剥蚀面与高速顶界面(图 9c和图 9d中虚线所示), 速度的剧烈变化往往会导致明显的波前畸变, 造成射线路径的“扭转”, 而MSFM算法由于考虑到了对角方向的走时信息, 在大网格间距下具有比FMM算法更高的精度, 更适用于复杂近地表的射线追踪。

3 结论

本文对传统基于三维FMM的射线追踪算法进行了两点改进。一是在走时正演中结合坐标变换将MSFM算法拓展到三维起伏地表情况, 通过同时考虑坐标轴方向和对角线方向上的信息, 有效改善了FMM算法的计算精度和计算效率; 二是结合三线性插值和Runge-Kutta方程优化了走时梯度的求取, 显著提高了射线追踪的计算精度与稳定性。模型试算表明本文方法能够为层析反演提供高效高精度的射线追踪, 有很好的应用价值。