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

2021-09-24李天扬李桐林孟相禹

世界地质 2021年3期

李天扬，李桐林，孟相禹

吉林大学地球探测科学与技术学院，长春 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绘制成图。