基于程函方程快速行进法的地震走时层析成像法

李天扬，李桐林，孟相禹

吉林大学 地球探测科学与技术学院，长春 130026

0 引言

地震走时层析成像法是一种利用不同地震台站连续接收到的地震事件，利用地震走时断层扫描反演地球内部速度分布的方法。由1977年发展至今已有40多年的历史，已成为如今对地下速度结构探测的重要研究手段，在火山活动、板块运动学和地幔流等领域的研究中扮演着重要的角色[1]。

地震层析成像过程主要包括正演和反演两个环节。早期的层析成像法正演采用的是传统的射线追踪方法，例如试射法、弯曲法，但由于其追踪过程中存在多路径问题,致使反演不稳定且多解，故针对该问题相继提出了基于程函方程数值解的地震波走时计算方法，主要包括有限差分法[2]、快速推进法[3--4]和快速扫描法[5--7]。近几十年来，最短路径法[8--10]等方法发展起来。反演环节早期的地震层析成像法中应用的阻尼最小二乘法、奇异值分解法，随后发展的模拟退火法[11]、共轭梯度法、迭代重建法及子空间反演法等诸多方法都可以针对不同模型应用于反演环节中。

地震层析成像中不同的波也具有不同的用处，其中最常见的有初至波、反射波。Rawlinson et al.[12]将FMM正演算法与Subspace反演算法相结合利用远震初至波资料研究了澳大利亚Murray盆地的上地幔速度结构；2015年，付翠[13]正演利用旅行时线性插值射线追踪法以及反演的最小二乘因式分解法将初至波与反射波联合反演应用到井间地震中；2017年,俞岱等[14]利用旅行时线性插值法以及初至波旅行时反演实现了初至波层析成像并行运算；查树贵[15]通过对反射波走时信息加约束条件提高成像的分辨率；Wang et al.[16]利用反射波振幅信息对模型速度的空间分布进行反演计算；张云等[17]在2019年基于多步快速行进法和子空间反演法实现了对二维起伏层状介质中的P波和反射波的反演。国内目前对于FMM算法以及利用其进行射线追踪不多见，孙章庆等[18--20]利用快速推进法迎风双线性插值法实现了三维地震波走时计算以及在复杂山地下和复杂海底中地震波的射线追踪和走时计算；王晗等[21]反演通过双曲Radon变换实现了对地震数据的重建；肖汉等[22]利用快速匹配法计算得到了VTI介质的走时。笔者在前人的工作基础上，正演采用基于程函方程快速行进法(FMM)进行射线追踪，实现了对入射波以及反射波的射线追踪，反演采用了子空间反演法，并且设计了理论模型进行反演，并对各反演结果进行了分析。

1 正演理论

1.1 快速行进法解程函方程

快速行进法是射线追踪方法的一种，主要用来求解非线性的程函方程来计算地震波走时场以及射线路径。程函方程是指导波前传播的物理现象表达式，其二维直角坐标系下的程函方程为：

(1)

式中：t为地震波的旅行时间；s为地震波在介质中传播的慢度。程函方程的解即为对t的求解，采用迎风差分方法对其离散并简化为：

(2)

(3)

Rawlinson and Sambridge将一、二、三阶的算子进行了精度上和计算效率上的比较，得到了一阶算子精度达不到标准，二阶、三阶的精度差异很小，但是二阶算子的计算效率远远高于三阶算子的结论，故通常应用时选择二阶快速行进法。二阶迎风差分算子为：

(4)

走时t可通过对式(2)的离散程函方程解得，其流程是选取需要计算的节点以及该节点上下左右4个方向结点并向波前拓展的方向建立差分格式，解出该方程的黏性解即为该节点的走时。

1.2 快速行进法解程函方程求取波前的方法

通过快速行进法得到的程函方程的解即为波前。由于快速行进法是一种基于网格的射线追踪方法，因此要描述其在网格中如何将波前推进。

如图1a所示，网格中所有的点划分为3个区域，分别为上风区域、窄带及下风区域。其中上风区域中的点是在波前的后面，称为alive点，其走时是已知；窄带中的点是在波前上，称为close点，走时虽然已经算出，但是还没有进行更新；下风区域中的点在波前的前面，称为far点，其走时还未进行计算。

图1b所示为窄带原理在源点进行波前演变的步骤。alive点的走时是已经计算出来的，使用离散化公式计算出close点走时，然后选择其中值最小的变成alive点，far点没有计算值，窄带是由计算出的close点不停地向前演变。在这个过程中，新的alive点不断的被标记，并持续更新邻近的close点和far点，窄带的传播形状近似为初至波前的形状，快速行进法使窄带沿着网格传播，计算结束的标志是所有的点变成alive点。

1.3 射线求取以及反射走时雅克比矩阵的计算

从接收点开始逐点寻找走时场梯度的最大下降方向，该方向点的连线即为接收点到源点之间的射线路径。首先要找到接收点位置并记下，然后设置每一步射线增量路径长度进行梯度的计算，计算完每一步之后得到的点记下并继续上一步的计算直到设置的界面位置(此过程采用的是已经计算过的反射波场)，记下界面上的点并继续以上操作直到源点为止(此过程采用的是已经计算过的入射波场)，连接此过程的点即为一条入射波和反射波的射线路径。

在二维直角坐标系中，模型任一点P可表示为(x,z)，根据射线理论，由第i个源点经反射点到第j个接收点的走时ti,j可表示为：

(5)

式中：ds为第i个震源到第j个检波器的射线路径Rij上的线积分元；v(x,z)为沿着该射线路径上的速度函数，当模型单元网格化后，沿着射线路径Rij上的走时可以写成求和的形式：

(6)

式中：N为射线穿过的单元总数；Rij,k、v(x,z)分别为射线穿过第k个单元的长度和速度分布。在对速度进行反演时，Frechet偏导数矩阵包含了走时关于速度变化的导数，其公式为：

(7)

当射线穿过某一个网格单元时，走时对速度的一阶偏导数可用射线两端点导数的平均值代替：

(8)

式中的：

(9)

式中：c=ki,kj，ki,kj是穿过第k个单元射线段的两个端点；BV(k)即为网格单元速度；fc(x,z)即为所求Frechet导数矩阵也就是雅克比矩阵。

2 反演理论以及算法流程

2.1 反演理论

子空间算法是Skilling et al.[23]在最大熵图像重建中提出的，随后被诸多学者将其概念推广和改进，并用其讨论了地震多参数的同时反演问题，主要用于对速度以及界面的反演中。

走时层析成像方法反演可归结为如下目标函数求极小值的问题，即：

(10)

子空间反演法是在多个模型参数的子空间内同时沿着多个方向去寻找目标函数S(m)的最小值，算法的核心是将多个参数的模型扰动量用p维子空间的一维基向量{aj},j=1,2,…，p来表示，即：

(11)

由公式(9)可推出：

(12)

(13)

(14)

②如果p≧2，则构造基向量公式如下：

(15)

③采用奇异值分解法将①、②中构建的基向量标准正交化。

相比于传统的其他梯度类算法，其优点是只需要对一个维数等于子空间维数的矩阵进行反演，并正确选择跨越模型的子空间的基向量，这个优点也使其实用性高于其他算法。

2.2 算法流程

依据快速行进法射线追踪理论以及子空间反演法为基础，在Windows环境下运用Fortran语言编程实现。实现了对二维模型快速行进法的射线追踪以及对模型的反演。其中射线追踪过程为：①给定初始模型并设置好网格范围以及源点接收点；②在源点附近进行源的细化计算，提高计算准确度；③设置网格所有点为far点并开始计算；④计算到网格底面使所有点计算成为alive点并备份计算所得网格时间场；⑤将计算结束的网格区域再次设置为far点，从保留了时间场的网格底面开始向接收点进行计算；⑥到接收点计算结束后完成对入射波以及反射波波场的计算；⑦如1.3所述完成对射线的追踪。反演流程如图2所示。正演反演结束后结果数据通过matlab和surfer绘制成图。

图2 反演流程图Fig.2 Inversion flow chart

3 模型试算

采用理论模型进行试算，网格剖分情况为26 km×30 km。异常体情况如图3、图4所示，低速异常体和高速异常体大小形状位置均相同，其中低速异常体背景速度为5 km/s，低速模型为速度3 km/s；高速异常体背景速度为5 km/s，高速模型速度为8 km/s。

图3 低速异常体理论模型Fig.3 Theoretical model of low velocity abnormal body

图4 高速异常体理论模型Fig.4 Theoretical model of high velocity abnormal body

图5、6、7为分别为低速异常体和入射波场、反射波场及射线路径图，由正演计算所得，在图中可以看出对应异常体的位置；图8、9、10分别为高速异常体的入射波场、反射波场及射线路径图，图中可以看出对应位置异常体的位置。

图5 低速异常体入射波场示意图Fig.5 Schematic diagram of incident wave field of low velocity abnormal body

图6 低速异常体反射波场示意图Fig.6 Schematic diagram of reflection wave field of low velocity abnormal body

图7 低速异常体射线路径示意图Fig.7 Schematic diagram of low velocity abnormal body ray path

图8 高速异常体入射波场示意图Fig.8 Schematic diagram of incident wave field of high velocity abnormal body

图9 高速异常体反射波场示意图Fig.9 Schematic diagram of reflection wave field of high velocity abnormal body

图11和图12为低速异常体和高速异常体的反演结果，由反演结果可以看出无论是低速或者是高速异常体的空间位置与给定模型基本吻合，但是异常体的形状和大小有些许差异，这是由于射线的密度和射线角度的覆盖率所导致。此外异常体的速度与背景速度的差异也会对反演结果产生影响。

图10 高速异常体射线路径示意图Fig.10 Schematic diagram of high velocity abnormal body ray path

图11 低速异常体反演结果Fig.11 Inversion results of low velocity abnormal body

图12 高速异常体反演结果Fig.12 Inversion results of high velocity abnormal body

4 结论

(1)快速行进法(FMM)不仅能追踪初至波，而且能追踪反射波。

(2)通过对初始模型以及异常体模型的反射波射线追踪，反射波能够反映速度结构的变化。

(3)子空间算法反演得到的反射波成像结果，能够准确的恢复出异常体的空间位置，但是由于射线密度和射线角度的原因导致异常体的形状与大小出现了些许差异。

(4)本文所设计固定界面进行射线追踪以及异常体反演，在实际地质条件下界面往往更为复杂，后续还需对复杂界面下的速度进行射线追踪和反演以及对界面进行反演。