2D/3D 起伏地表多震相地震波走时的因式分解程函方程算法

张云，李夕海＊，白超英，牛超，王艺婷，曾小牛

（1.火箭军工程大学，陕西西安 710025； 2.长安大学地质工程与测绘学院地球物理系，陕西西安 710054； 3.长安大学计算地球物理研究所，陕西西安 710054）

0 引言

基于几何地震学的地震波走时计算技术具有算法便捷、计算效率高等优点而被广泛应用于地震层析成像［1-4］、地震定位［5--7］等领域。地震波走时计算方法主要可分为三种：一是基于运动学方程的射线追踪类方法［8-9］，但当遇到复杂介质时，这类方法容易陷入局域解，存在阴影区问题； 二是基于费马原理的最短路径算法（Shorest Path Method，SPM）［10-11］，这类方法通过采用大量的网格线段近似逼近地震波射线，为保证结果精度，需用比其他方法更庞大数量的网络节点，这也导致了更高的计算空间要求； 三是基于有限差分求解程函方程类方法，它有效避免了运动学方程类方法的不足，且在计算精度与计算效率方面具有显著优势，自提出以来即发展迅速，作为正演算法多次解决地震波走时反演问题［12-14］，近年有不少学者将其推广至各向异性介质［15-20］。

快速行进法（Fast Marching Method，FMM）［21］与快速扫描法（Fast Sweeping Method，FSM）［22］均是通过有限差分求解程函方程，从而得到全局地震波走时的计算方法。FMM 采用窄带技术近似波前扩展，即隐含了这样一个假设：群速度给出的能量传播方向与相速度给出的波前传播方向相同。然而，该假设仅对各向同性介质是成立的，对于各向异性介质则不成立。虽然近年来也有不少学者通过加入其他约束性条件将其推广至各向异性介质［23-25］，但过多的约束条件也使得算法变得复杂且低效。因此，对于各向异性介质而言，人们更愿意采用FSM。FSM 由Zhao［22］率先提出，通过沿交替方向扫描计算域求解程函方程，在有限次数的迭代中收敛得到数值解，比FMM 更灵活且稳健性更高。有关FMM 与FSM 的详细比较，可参考文献［26］。尽管FSM 在各向异性介质及稳健性上明显优于FMM，但研究表明［27-28］，对于复杂模型或网格节点较多的模型而言，FSM 计算效率低于FMM。主要原因是对于较复杂模型FSM 需更多的迭代次数才能收敛，对于具有大量网格节点的复杂模型，此问题更严重。如Li等［29］为保证反演效率，选择FMM 作为正演算法以解决反演算法的效率问题；Benaichouche 等［30］则选用FSM 作为正演地震波走时计算方法，然而在反演计算雅可比矩阵时，选用FMM求解原始程函方程，这样虽然提高了效率，但在一定程度上牺牲了计算精度。

FMM 是通过有限差分求解程函方程，采用窄带近似波前而扩展起来的一种方法，Rawlinson 等［31］将FMM 与分区多步计算技术（Multistage Computational Technique）相结合，实现了2D 复杂地质结构下的多震相地震波走时计算（Multistage FMM）。随后，de Kool 等［32］考虑地球曲率影响，又将该算法推广至球坐标系下进行多震相地震波走时计算。孙章庆等［33］基于FMM，提出不等距迎风差分格式（Uneven Grid Upwind Finite-Difference Formula）计算3D 起伏地表条件的直达波走时，随后又引入群行进法（Group Marching Method）［34］以提高FMM 计算效率［35］。然而，震源点奇异性问题一直伴随着FMM 的发展。为解决此问题，Rawlinson 等［31］采用源点网格加密技术，分析表明，经过对震源附近网格细化后的FMM 可明显降低源点附近的计算误差，但相对于整个计算区域，源点附近的计算误差仍然较大，因此震源区网格细化并不能根本解决震源点的奇异性问题。此外，震源区网格细化的同时也意味着额外增加了计算量。Alkhalifah 等［36］则在球坐标系下求解程函方程，但这种算法复杂，实际应用中难以推广。Fomel等［37］提出将地震波走时因式分解，用FSM 求解因式分解程函方程以解决震源奇异性问题，数值模拟表明，求解因式分解程函方程数值解的准确性高于直接计算原始程函方程，且震源附近改善效果更明显。周小乐等［38］推导建立了曲线坐标系因式分解程函方程，结合FSM，在解决震源奇异性问题的同时，可得到起伏地表条件的地震波走时。

鉴于FMM 对复杂模型计算效率的优势及求解因式分解程函方程可解决震源奇异性问题，本文提出了一种起伏地表条件下基于FMM 求解因式分解程函方程的数值算法，同样可从根源上解决震源奇异性问题。与传统FMM 相同，本文所提算法可在保证O（NlgN，N为网络节点数） 计算成本内使用一阶或二阶有限差分公式求解因式分解程函方程。同时，考虑到加入后续震相约束可显著提高反演分辨率，本文又结合分区多步计算技术实现了多震相地震波走时的计算。

文中首先阐释如何建立因式分解程函方程； 然后介绍利用FMM 求解因式分解程函方程； 最后通过2D/3D 均匀模型、速度线性增加模型和起伏界面模型验证了所提算法的有效性与鲁棒性。为表述方便，将本文所提算法称为基于因式分解程函方程FMM，简称FMM-F。

1 因式分解程函方程的导出

限于篇幅，在推导因式分解程函方程过程中，仅展示2D 直角坐标系情形。对于3D，只需进行类似的扩展即可。首先引入2D 直角坐标系程函方程

式中T（x，z）和S（x，z）均为空间位置坐标（x，z）的函数，分别表示震源点O到空间点（x，z）的地震波走时与该点的波慢度。将T（x，z）分解为一个关于点（x，z）到震源点O的距离函数T0与一个地震波走时扰动量T1的乘积形式，即

这样，根据链式求导法则，标准程函方程式（1）便可写成

式中T0可通过解析求得，T1则需通过数值求解得到，再通过式（2）便可精确求得地震波走时，从而规避震源奇异性问题。本文选择T0为震源(x0，z0)与待求点（x，z）之间的距离函数，即

则对于震源点(x0，z0)有

对于任意一点，有

则对于震源点(x0，z0)，由因式分解程函方程式（3），可得

上述T0(x0，z0)与T1(x0，z0)即为本算法的初始化条件。

2 数值求解T1

传统FMM 采用有限差分近似走时梯度，利用窄带模拟波前扩展，采用堆排序算法寻找波前中走时最小点作为次级震源。本文算法同样基于以上技术，这些技术已在许多FMM 算法相关文献［21，31-32］中有详细阐述，此处不再赘述。本文算法的关键在于如何利用FMM 求解因式分解程函方程，下面分别阐述两种方法求解T1。

2.1 规则网格求解T1

首先，对因式分解程函方程进行离散，则有（为表述方便，用点（x，z）的编号（i，j）表示其位置）

式中：max（A，B，C）表示求取A、B、C中最大值；和分别表示空间点（i，j）处向l负方向和l正方向的n阶差分算子。数值计算时分别采用一阶或二阶差分进行计算。严格地讲，二阶差分属于混合差分，因为当二阶差分条件不满足时，即采用一阶差分。以下以x方向为例介绍一阶和二阶差分求解因式分解程函方程过程。

规则网格（图1）中，假定此时（i，j）点的走时为待求点走时，引入符号Tknow表示T值已知、Tunknow表示T值未知。对于一阶差分，可得到式（9）； 根据链式求导法，则有式（10）

图1 规则网格走时扰动因子T1求解示意图

式中：T1(i-1，j)和T1(i+1，j)分别为(i-1，j)点和(i+1，j)点的走时扰动因子，可结合T(i-1，j)和T(i+1，j)由式（2）求出；T0(i，j)可由式（4）解析求出；h为网格剖分间距；。这样，式（10）中就只有一个待求未知量T1(i，j)。

对于二阶差分，实际为混合阶数差分（当不满足二阶差分条件时，即运用一阶差分）。与一阶差分相似，仍需先判断待求点(i，j)邻近点状态。同样，以x方向为例，具体如下

式中各变量物理意义与一阶差分（式（9））类似。

z方向差分算子的推导，与x方向类似。至此，再结合窄带技术［21］，便可实现直达波走时计算。

2.2 不规则网格求解T1

为提高走时反演成像的分辨率，采用多震相走时联合反演是一个不错的选择，这就要求正演算法具有计算多震相地震波走时的能力。分区多步计算技术自提出以来，就与多种网格类地震波走时计算方法相结合，用于模拟多震相地震波传播。本算法同样作为网格单元类算法，也可结合分区多步计算技术实现多震相地震波走时的计算。

然而实际地表（或地下界面）多呈不规则性，若用规则网格剖分，则地下界面（或地表）附近会产生不规则网格（图2b 中M1M2界面）。这样，上述规则网格的迎风差分公式显然已不适用，因此需要针对界面附近不规则网格，推导建立新的局部走时计算公式。图2不规则网格中的点可划分为两类：一类是界面附近点（C点）； 另一类是界面点（M点）。下面分别针对这两类点建立局部走时计算公式。

图2 模型界面及网格剖分示意图（a）及其局部放大图（b）

（1）C点为当前待计算点

此时以C点为中心的z方向出现CM≠h，即在C点的x轴方向与z轴上半轴方向仍可采用常规迎风差分格式，而z轴下半轴常规迎风差分格式显然已不能适用，为此建立C点z轴向下方向的不等距一阶、二阶迎风差分公式（记|CM|=h1）

而二阶差分算子的推导过程如下，利用泰勒函数将T1(M)与T1(D)在C点附近展开

将式（15）代入式（14），可得不等距二阶迎风差分公式。

此外，由于|CM|≠BC，导致T(M)与T(B)的大小已无可比性，则规则网格的迎风差分条件在此已不适用，因此建立如下不等距迎风差分格式

（2）M点为当前待计算点

此时M点周围的四个网格均为非正交，导致程函方程不成立。Rawlinsnon 等［31］采用三角网格缝合界面，而后在三角网格中求解程函方程，但实现较为困难。本文综合不等距迎风差分公式、费马原理与惠更斯原理，求解界面走时。对于M点的z轴下方向，同样有不等距一阶、二阶差分公式

式中h2为M、D之间的距离。对于M点z轴上方向，表达式类似。

此外对于M点x方向，由于界面起伏，x方向常规差分公式已不再适用。本文采用惠更斯原理，即将M1、M2视为子波源，就有

式中：li为Mi与M之间的距离；MiM为Mi与M之间的波慢度。

最后，再利用费马原理，得到

综上所述，由于起伏界面的存在，可将界面上的点分为三类：第一类为规则网格中节点，采用常规迎风差分公式（式（8）～式（10））计算修正因子T1，然后根据式（2）计算该点地震波走时； 第二类为界面附近节点（C），对于该点远离界面方向和x方向，走时扰动因子的计算方法与规则网格相同，对于该点的界面方向，采用式（13）计算该方向的走时扰动因子，最后根据费马原理，由式（16）计算该点走时；第三类为界面上点（M），对于该点z轴方向，建立不等距迎风差分公式（式（17））计算该点在z方向的走时扰动因子，最后结合费马原理、惠更斯原理，由式（19）计算M点地震波走时。

3 多震相地震波走时计算方法

基于前面构建的起伏地表模型下规则网格与不规则网格的迎风差分格式，再结合分区多步技术与窄带技术，便可实现多震相地震波走时计算。所谓分区多步计算技术即指按照地下起伏界面的具体情况，将模型分成独立的层状（或起伏层状）区域，相邻区域由地下界面相连接［18］（图2a）。窄带技术根据模型中网格节点所处状态将所有节点分为三种：“完成点”、“窄带点”及“远离点”。其中“完成点”表示走时计算已结束的网格点；“窄带点”表示当前波前上的网格节点，其走时大小还可更新；“远离点”表示走时计算还未实施的网格节点。

表1给出了本文所提FMM-F 算法计算直达波走时的实现框图，当要计算其他震相的地震波走时时，只需将界面上走时最小的点作为新的子震源，而后与计算直达波相同，计算当前所需震相的地震波走时即可。

表1 FMM-F 算法实施过程

现以图2a所示模型计算上界面反射波为例，说明反射波走时计算原理（震源所在位置为第一分区），其他震相的走时计算也可类推。

（1）首先设定算法初始化条件，给定T0(x0，z0)=0，T1(x0，z0)=S(x0，z0)，并根据窄带技术将震源点记为“完成点”，即该点地震波走时不再更新；

（2）判断当前震源点所在网格是否包含界面或起伏地表：若是，则采用非规则网格迎风差分格式（式（13）～式（19））计算震源周围点走时扰动因子T1；若否，则采用规则网格迎风差分格式（式（8）～式（12））计算T1，通过式（2）和式（4）得到震源周围点地震波走时，这些点构成的面即为当前波前所在位置，其状态记为“窄带点”；

（3）在“窄带点”中搜寻走时最小点，将此点作为次级震源并将其状态更新为“完成点”继续循环第（2）步计算，直到得到震源所在分区的全局走时；

（4）读取界面离散点地震波走时并将其状态均赋为“窄带点”，返回第（3）步，要注意的是由于是计算反射波，所以在第（2）步时只需考虑界面以上的网格节点即可。

若要追踪透射波，则在第（4）步返回第（3）步时，自界面上走时最小的点向透射区域进行波前扩展即可； 若要追踪转换波，则自离散界面上走时最小的点开始向转换波所在区域进行波前扩展扫描时时，采用不同的速度模型即可。简而言之，分区多步计算的原理就是根据所要计算地震波的类型，设计要调用的独立区域的先后顺序及对应的速度模型，达到运用简单组合完成复杂模型中多震相地震波走时计算之目的。

4 数值模拟及结果分析

采用几组模型进行数值模型计算，通过本算法与Multistage FMM 算法［31］详尽（主要是计算精度与计算效率）对比，验证本文算法的优越性。当选用模型存在解析解时，将以解析解作为参考解； 当选用模型不存在解析解时，将以发展较成熟的不规则最短路径算法（ISMP）［37］作为参考解，且通过加密ISMP算法网格间距，确保该算法计算结果的精度。计算所用处理器配置为： Intel（R） Core（TM） i7-10750H CPU @2.60 GHz。另外，在分析算法所用CPU 时间时，均是设置程序循环执行100 次后取均值得到每次所用CPU 时间，以避免偶然性。

4.1 模型1（均匀模型）

2D 均匀模型大小为4 km×4 km，地震波传播速度为1 km/s，震源位于（2 km，2 km）；3D 均匀模型大小为10 km×10 km×10 km，地震波传播速度为2 km/s，震源位于（5 km，5km，5 km）。对于FMM，采用震源网格加密方式保证精度，加密区域范围为[x0-1 km，x0+1 km]、[z0-1 km，z0+1 km]（2D）或[x0-1 km，x0+1 km]、[y0-1 km，y0+1 km]、[z0-1 km，z0+1 km（]3D），(x0，y0，z0）（为震源位置坐标，加密网格间距为0.005 km。本文选用平均绝对误差与平均相对误差两个指标评估计算精度，平均绝对误差计算公式为：；同时，平均相对误差的计算公式为：。其中Texact为解析解，Tcal为通过FMM 或FMM-F 所计算的结果，N为模型中网格节点总数。

表2（2D）和表3（3D）分别给出均匀模型不同网格间距剖分下FMM 与FMM-F 的计算结果对比。从该表可见，对于均匀模型，本文所提FMM-F计算结果误差为零（注意：误差分析结果已消除因计算机自身小数保留位数而导致的计算误差）； 同时，可看到随着网格间距的减小，两种算法的计算效率都降低了，但在相同的网格剖分间距下，FMM-F 计算效率明显高于FMM，考虑是由于FMM 采用震源网格加密方式保证精度，从而使得FMM 计算效率降低。另外，当采用二阶差分时，两种算法计算效率都降低了，但降低并不明显。

表2 2D 均匀模型中FMM 与FMM-F 计算地震波走时误差及CPU 时间对比

表3 3D 均匀模型中FMM 与FMM-F 计算地震波走时误差及CPU 时间对比

4.2 模型2（速度线性增加模型）

再采用速度线性增加模型（图3）进行试算。2D模型（图3a）大小为8 km×4 km，速度自地表4 km/s线性增至底界面的6 km/s，震源位于(4 km，0)； 3D模型（图3b）大小为10 km×10 km×10 km，速度自地表2 km/s 线性增至底界面的4 km/s，震源位于（5 km，5 km，0）。该模型的地震波走时解析解由下式求出［37］

图3 速度线性增加的2D（a）和3D（b）模型

式中：α为速度增加梯度；v0为震源处地震波传播速度；v(x)为x点的速度；x0为震源位置坐标。

表4 和表5 分别给出不同网格间距剖分计算2D和3D 线性模型时FMM 与FMM-F 的计算结果对比。同样采用震源网格加密方式保证FMM 计算精度。水平轴加密区域范围为[x0-1 km，x0+1 km（]2D）或[x0-1 km，x0+1 km]、[y0-1 km，y0+1 km（]3D），由于此模型震源置于地表，所以深度方向加密区域范围为[z0，z0+1 km]，（x0，y0，z0）为震源位置坐标，加密网格间距为5 m。同样，选用平均绝对误差与平均相对误差两个指标评估计算精度。

表4 图3a 模型中FMM 与FMM-F 计算地震波走时误差及CPU 时间对比

表5 图3b 模型中FMM 与FMM-F 计算地震波走时误差及CPU 时间对比

从表4和表5可见，随着网格间距的不断减小，不论是FMM 还是FMM-F，其平均绝对误差与平均相对误差都不断减小，计算精度不断提高。但同时，网格间距的减小，导致网格节点数变多，计算效率也随之下降。另外还发现，当采用二阶差分时，两种算法的计算精度都有大幅度提高，且与均匀模型相同，计算效率降低并不明显，因此建议实际应用中采用二阶差分格式。重要的是在两种算法对比中发现，网格间距相同时，2D 或3D 情形下，无论是计算精度还是CPU 计算时间，FMM-F 优于FMM。这是因为计算精度方面，FMM-F 成功地避免了震源点的奇异性问题； 而计算效率方面则是传统FMM 为了克服震源点的奇异性需对震源附近网格细化，导致计算量增加。

为进一步分析误差分布情况，又给出当h=0.05 m 时两种算法分别采用一阶（图4（2D）、图6（3D））、二阶（图5 （2D）、图7（3D））差分所得的绝对误差和相对误差分布图。可发现传统FMM 计算所得地震波走时在震源对角线附近区域误差较大，当采用一阶差分格式时，其最大绝对误差分别可达4.80 ms（2D）和20 ms （3D），最大相对误差则分别高达7%（2D）和4% （3D）； 而FMM-F 最大绝对误差仅分别为0.16 ms （2D）和1.5 ms （3D），而最大相对误差也分别仅为0.016%（2D）和0.04%（3D）。当采用二阶差分格式时，两种算法的计算精度都得到了很大提高，FMM 最大绝对误差分别降至0.4 ms（2D）和1.5 ms（3D），同样最大相对误差也分别减小到2.7%（2D 和2%（3D）； 而FMM-F 最大绝对误差则分别进一步减至0.04 ms（2D）和0.1 ms（3D），最大相对误差则分别减至0.009%（2D）和0.02%（3D）。可见，本文的FMM-F 不论采用一阶还是二阶差分格式，其计算精度高于传统FMM； 同时，通过误差分布图可发现FMM-F完全解决了震源点的奇异性问题。

图4 2D 模型FMM（a、b）、FMM-F（c、d）算法一阶差分计算的绝对误差（a、c）和相对误差（b、d）

图5 对应图4 的二阶差分计算结果

图6 3D 模型中FMM（上）和FMM-F（下）算法采用一阶差分格式计算误差分布

图7 对应图6 采用二阶差分格式计算结果

4.3 模型3（起伏界面模型）

在上文分析了直达波的走时计算精度的基础上，通过一个数值例子验证本文算法的多震相走时计算精度。因没有解析解，故将Multistage ISPM 算法［39］的计算结果作为参考值。为了保证Multistage ISPM算法较高的计算精度，网格主节点间距设置为0.25 km，且在主节点间加入19 个次级节点，这样通过Multistage ISPM 算法求取结果的精度较高，可作为参考值。所用模型如图8a （2D）、图9a （3D）所示，其大小分别为100 km×50 km （2D） 和100 km×100 km×50 km （3D）。

图8 2D 起伏模型计算精度分析

图9 3D起伏模型及不同剖面的速度分布

2D 模型地表的起伏函数设置为：z=2+2×sin (2πx/50 )，且在25～35 km 深度范围内含有一个倾斜断层，断层大约位于x=50km（图8a 中黑线），震源置于点（50 km，15 km）处，并在地面以2 km 间距均匀布放51个检波器。模型速度由地表的2 km/s线性增至断层界面的4 km/s，而后又线性增至底界面的6 km/s（图8a）。

由于常规Multistage FMM［31］采用三角网格离散界面，为满足迎风差分条件，需对钝角三角网格进行剖分，导致常规Multistage FMM 对界面的离散与Multistage ISPM 算法和本文算法不一致，且无可比性，故此部分不再将Multistage FMM 作为对比的算法，而仅将Multistage ISPM 算法作为参照，以验证本文算法的有效性。

图8c、图8d 分别为采用图8a 的炮检排列与速度结构模型时，自震源激发后的直达波和自底界面反射上行波的走时等值线图，可见Multistage ISPM 算法计算结果对应的黑色虚线、本文算法采用一阶差分计算结果对应的红色虚线及采用二阶差分对应的绿色虚线三种等值线的重合度较高，表明本文算法在计算直达波与反射波时精度较高。为进一步突显对比效果，还对等值线图做了局部放大显示（图8c、图8d右下角）及与地面检波器观测反射波到时的相对误差（图8b，计算公式为，通过局部放大图清晰可见采用二阶差分时的等值线与Multistage ISPM 算法的重合度更高，表明采用二阶差分的计算精度更高，通过与地面检波器的相对误差对比，可得出同样结论。

针对3D 模型，模拟计算了8 种震相地震波的走时，并以Multistage ISPM 算法作为标准计算了相对误差（误差计算公式与2D 相同，图10）。它们分别为经上界面反射后的纯反射波P1P （图10a）、S1S（图10b）与反射转换波P1S （图10c）、S1P （图10d），经底界面反射后的纯反射波P2P（图10e）与反射转换波S2P （图10f），以及经上界面与下界面多次透射、反射的转换波P2S1P2S、S2S1P2P。对比图10上图与下图可见，经上界面反射后的地震波走时的相对误差明显大于经底界面反射的地震波走时，这是缘于上界面存在一个较大孔洞，而底界面仅为一平面； 同时，8 种地震波震相的走时相对误差均在0.4%以内，可见本文算法对复杂起伏地表（含地下孔穴等剧烈起伏地质异常体）模型中多震相地震波走时的计算精度也较高，可作为一种高效、计算精度高的多震相地震波走时场计算的有效方法。

图10 3D 起伏模型地表检波器接收的多震相走时的相对误差分布

5 结论

由于原始程函方程的FMM 采用平面波近似波前，其波前曲率较大，导致震源奇异性问题，且利用三角形（2D）或四面体（3D）缝合界面，实现较为困难。为此，本文发展了一种利用快速行进法求解因式分解程函方程的多震相地震波走时计算算法。该方法采用因式分解求解程函方程解决震源奇异性问题，以不等距差分格式精确定位地表位置，通过迎风混合网格差分公式计算地震波走时，并用窄带技术作为算法的整体波前扩展技术。数值模拟实例表明，以解析解、Multistage ISPM 算法［31］作为参考解，与Multistage FMM 算法［18］对比分析，得到以下认识：

（1）所提算法计算精度较高，不论是平均绝对误差还是平均相对误差都显著降低；

（2）通过弃用震源网格细化以保证计算精度，且该算法计算效率较高；

（3）利用该算法可灵活、稳定地计算复杂地表条件下的走时。