2D_3D起伏地表多震相地震波走时的因式分解程函方程算法.pdf

2023 年 8 月第 58 卷 第 4 期2D/3D起伏地表多震相地震波走时的因式分解程函方程算法张云1，李夕海1，白超英2，3，牛超1，王艺婷1，曾小牛1（1.火箭军工程大学，陕西西安 710025；2.长安大学地质工程与测绘学院地球物理系，陕西西安 710054；3.长安大学计算地球物理研究所，陕西西安 710054）摘要：起伏地表条件的地震波走时计算方法是研究该类地表区地下结构的基础工具。快速行进法和快速扫描法均是基于有限差分求解程函方程而发展起来的地震波走时计算方法，由于震源附近波前曲率较大，这两种算法均存在震源奇异性问题。研究成果表明，对于复杂模型快速行进法的计算效率高于快速扫描法。为此，借鉴快速扫描法解决震源奇异性的思路，采用快速行进法求解因式分解程函方程，从而规避了震源奇异性问题。具体而言，将地震波走时分解为一个距离函数 T0与一个走时扰动值 T1乘积的形式，通过快速行进法求解 T1，并与T0相乘，得到地震波走时；同时，为弥补规则网格迎风差分格式不适用于地表/界面起伏的缺陷，构建了适用于不规则网格的不等距迎风差分格式，进而结合分区多步技术形成了全局多震相地震波走时计算方法。2D 和 3D 数值模拟实例表明，所提算法从根本上解决了快速行进法的震源奇异性问题，显著提高了原算法的计算精度与效率，可精确计算多震相地震波走时。关键词：因式分解程函方程，走时扰动因子，不等距迎风差分格式，多震相地震波走时计算中图分类号：P631 文献标志码：A doi：10.13810/ki.issn.1000-7210.2023.04.015Multiphase seismic traveltime computation in 2D/3D undulated surface model using factored eikonal equationZHANG Yun1，LI Xihai1，BAI Chaoying2，3，NIU Chao1，WANG Yiting1，ZENG Xiaoniu1（1.Rocket Force University of Engineering，Xi an，Shaanxi 710025，China；2.Department of Geophysics，College of Geological Engineering and Geomatics，Chang an University，Xi an，Shaanxi 710054，China；3.Institute of Computational Geophysics，Chang an University，Xi an，Shaanxi 710054，China）Abstract：The seismic traveltime computation scheme in undulating surface conditions is a basic tool to study the underground structures of such surface areas.The fast marching method（FMM）and the fast sweeping method（FSM）are both developed based on solving the eikonal equation with finite difference.They have the problem of source singularity due to the high curvature of the wavefront around the source.Previous studies show that the computational efficiency of FMM is higher than FSM for complex models.Thus，this paper employs the FMM to solve the factorization equation and avoid source singularity.Specifically，the original eikonal equation can be transformed into the factored eikonal equation，in which the seismic traveltime can be regarded as the product of a distance function T0 and a correction factor T1 of traveltime.The correction factor T1 of traveltime can be solved by the FMM algorithm and then the distance function T0 is multiplied to obtain the traveltime（T）.To address the problem that the upwind finite difference formula with even grid spacing is not applicable to surface/interface undulation，this paper constructs an upwind finite difference formula with uneven grid spacing.Finally，the multistage computational technique is adopted to propose a computation method for global multiphase seismic traveltime.The simulation tests indicate that the new algorithm solves the source singularity of FMM，significantly improves the computational accu 地震模拟 文章编号：1000-7210（2023）04-0857-15*陕西省西安市灞桥区同心路 2号火箭军工程大学，710025。Email：xihai_本文于 2022年 8月 30日收到，最终修改稿于 2023年 5月 15日收到。石 油 地 球 物 理 勘 探2023 年racy and efficiency of the original algorithm，and can accurately calculate the multiphase seismic traveltime.Keywords：factored eikonal equation，correction factor of traveltime，upwind finite difference formula with uneven grid spacing，multiphase seismic traveltime computation张云，李夕海，白超英，等.2D/3D起伏地表多震相地震波走时的因式分解程函方程算法 J.石油地球物理勘探，2023，58（4）：857-871.ZHANG Yun，LI Xihai，BAI Chaoying，et al.Multiphase seismic traveltime computation in 2D/3D undulated surface model using factored eikonal equation J.Oil Geophysical Prospecting，2023，58（4）：857-871.0引言基于几何地震学的地震波走时计算技术具有算法便捷、计算效率高等优点而被广泛应用于地震层析成像14、地震定位57 等领域。地震波走时计算方法主要可分为三种：一是基于运动学方程的射线追踪类方法89，但当遇到复杂介质时，这类方法容易陷入局域解，存在阴影区问题；二是基于费马原理的最短路径算法（Shorest Path Method，SPM）1011，这类方法通过采用大量的网格线段近似逼近地震波射线，为保证结果精度，需用比其他方法更庞大数量的网络节点，这也导致了更高的计算空间要求；三是基于有限差分求解程函方程类方法，它有效避免了运动学方程类方法的不足，且在计算精度与计算效率方面具有显著优势，自提出以来即发展迅速，作为正演算法多次解决地震波走时反演问题1214，近年有不少学者将其推广至各向异性介质1520。快速行进法（Fast Marching Method，FMM）21 与快速扫描法（Fast Sweeping Method，FSM）22 均是通过有限差分求解程函方程，从而得到全局地震波走时的计算方法。FMM 采用窄带技术近似波前扩展，即隐含了这样一个假设：群速度给出的能量传播方向与相速度给出的波前传播方向相同。然而，该假设仅对各向同性介质是成立的，对于各向异性介质则不成立。虽然近年来也有不少学者通过加入其他约束性条件将其推广至各向异性介质2325，但过多的约束条件也使得算法变得复杂且低效。因此，对于各向异性介质而言，人们更愿意采用 FSM。FSM由 Zhao22 率先提出，通过沿交替方向扫描计算域求解程函方程，在有限次数的迭代中收敛得到数值解，比 FMM更灵活且稳健性更高。有关FMM与FSM的详细比较，可参考文献 26。尽管FSM在各向异性介质及稳健性上明显优于FMM，但研究表明2728，对于复杂模型或网 格 节 点 较 多 的 模 型 而 言，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 FiniteDifference Formula）计算 3D 起伏地表条件的直达波走时，随后又引入群行进法（Group Marching Method）34 以 提 高 FMM 计 算 效率35。然而，震源点奇异性问题一直伴随着 FMM的发展。为解决此问题，Rawlinson等31 采用源点网格加密技术，分析表明，经过对震源附近网格细化后的FMM 可明显降低源点附近的计算误差，但相对于整个计算区域，源点附近的计算误差仍然较大，因此震源区网格细化并不能根本解决震源点的奇异性问题。此外，震源区网格细化的同时也意味着额外增加了计算量。Alkhalifah 等36 则在球坐标系下求解程函方程，但这种算法复杂，实际应用中难以推广。Fomel等37 提出将地震波走时因式分解，用 FSM 求解因式分解程函方程以解决震源奇异性问题，数值模拟表明，求解因式分解程函方程数值解的准确性高于直接计算原始程函方程，且震源附近改善效果更明显。周小乐等38 推导建立了曲线坐标系因式分解程函方程，结合FSM，在解决震源奇异性问题的同时，可得到起伏地表条件的地震波走时。858第 58 卷 第 4 期张云，等：2D/3D起伏地表下多震相地震波走时的因式分解程函方程算法鉴于 FMM 对复杂模型计算效率的优势及求解因式分解程函方程可解决震源奇异性问题，本文提出了一种起伏地表条件下基于FMM求解因式分解程函方程的数值算法，同样可从根源上解决震源奇异性问题。与传统 FMM 相同，本文所提算法可在保证 O（NlgN，N 为网络节点数）计算成本内使用一阶或二阶有限差分公式求解因式分解程函方程。同时，考虑到加入后续震相约束可显著提高反演分辨率，本文又结合分区多步计算技术实现了多震相地震波走时的计算。文中首先阐释如何建立因式分解程函方程；然后介绍利用FMM求解因式分解程函方程；最后通过2D/3D均匀模型、速度线性增加模型和起伏界面模型验证了所提算法的有效性与鲁棒性。为表述方便，将本文所提算法称为基于因式分解程函方程 FMM，简称FMMF。1因式分解程函方程的导出限于篇幅，在推导因式分解程函方程过程中，仅展示2D直角坐标系情形。对于3D，只需进行类似的扩展即可。首先引入2D直角坐标系程函方程|T(x，z)|2=S2(x，z)（1）式中 T（x，z）和 S（x，z）均为空间位置坐标（x，z）的函数，分别表示震源点 O 到空间点（x，z）的地震波走时与该点的波慢度。将T（x，z）分解为一个关于点（x，z）到震源点O的距离函数T0与一个地震波走时扰动量T1的乘积形式，即T=T0T1（2）这样，根据链式求导法则，标准程函方程式（1）便可写成|T0T1+T1T0|2=S2(x，z)（3）式中T0可通过解析求得，T1则需通过数值求解得到，再通过式（2）便可精确求得地震波走时，从而规避震源奇异性问题。本文选择T0为震源(x0，z0)与待求点（x，z）之间的距离函数，即T0(x，z)=(x-x0)2+(z-z0)2（4）则对于震源点(x0，z0)有T0(x0，z0)=0（5）对于任意一点，有(T0)2=1（6）则对于震源点(x0，z0)，由因式分解程函方程式（3），可得T1(x0，z0)=S(x0，z0)（7）上述T0(x0，z0)与T1(x0，z0)即为本算法的初始化条件。2数值求解T1传统FMM采用有限差分近似走时梯度，利用窄带模拟波前扩展，采用堆排序算法寻找波前中走时最小点作为次级震源。本文算法同样基于以上技术，这些技术已在许多 FMM 算法相关文献21，3132 中有详细阐述，此处不再赘述。本文算法的关键在于如何利用FMM 求解因式分解程函方程，下面分别阐述两种方法求解T1。2.1规则网格求解T1首先，对因式分解程函方程进行离散，则有（为表述方便，用点（x，z）的编号（i，j）表示其位置）maxD-xnT()i，j，D+xnT()i，j，02+maxD-znT()i，j，D+znT()i，j，02=S2(x，z)式中：max（A，B，C）表示求取 A、B、C 中最大值；D-lnT(i，j)和D+lnT(i，j)分别表示空间点（i，j）处向l负方向和l正方向的n阶差分算子。数值计算时分别采用一阶或二阶差分进行计算。严格地讲，二阶差分属于混合差分，因为当二阶差分条件不满足时，即采用一阶差分。以下以x方向为例介绍一阶和二阶差分求解因式分解程函方程过程。规则网格（图 1）中，假定此时（i，j）点的走时为待求点走时，引入符号Tknow表示T值已知、Tunknow表示T图 1规则网格走时扰动因子 T1求解示意图黑圆点为待求点，灰圆点为一阶差分，小方块为二阶差分。（8）859石 油 地 球 物 理 勘 探2023 年值未知。对于一阶差分，可得到式（9）；根据链式求导法，则有式（10）max(x)=D-x1T(i，j)当 Tknow(i-1，j)，Tunknow(i+1，j)或 Tknow(i-1，j)，Tknow(i+1，j)且 Tknow(i-1，j)Tknow(i+1，j)D+x1T(i，j)当 Tknow(i+1，j)，Tunknow(i-1，j)或 Tknow(i-1，j)，Tknow(i+1，j)且 Tknow(i+1，j)Tknow(i-1，j)0 当 Tunknow(i-1，j)，Tunknow(i+1，j)D-x1T(i，j)=T0(i，j)T1(i，j)-T1(i-1，j)h+P(i，j)T1(i，j)D+x1T(i，j)=T0(i，j)T1(i+1，j)-T1(i，j)h+P(i，j)T1(i，j)（10）式中：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为网格剖分间距；P(i，j)=T0(i，j)x。这样，式（10）中就只有一个待求未知量T1(i，j)。对于二阶差分，实际为混合阶数差分（当不满足二阶差分条件时，即运用一阶差分）。与一阶差分相似，仍需先判断待求点(i，j)邻近点状态。同样，以 x方向为例，具体如下max(x)=D-x1T(i，j)当 Tknow(i-1，j)，Tunknow(i+1，j)，Tunknow(i-2，j)或 Tknow(i-1，j)，Tknow(i+1，j)，Tunknow(i-2，j)，且Tknow(i-1，j)Tknow(i+1，j)D+x1T(i，j)当 Tknow(i+1，j)，Tunknow(i-1，j)，Tunknow(i+2，j)或 Tknow(i+1，j)，Tknow(i-1，j)，Tunknow(i+2，j)，且Tknow(i+1，j)Tknow(i-2，j)或 Tknow(i-1，j)，Tknow(i+1，j)，Tknow(i-2，j)，且Tknow(i+1，j)Tknow(i-1，j)Tknow(i-2，j)D+x2T(i，j)当 Tknow(i+1，j)Tknow(i+2，j)，且Tknow(i+1，j)Tknow(i+2，j)或 Tknow(i+1，j)，Tknow(i-1，j)，Tknow(i+2，j)，且Tknow(i-1，j)Tknow(i+1，j)Tknow(i+2，j)0 当 Tunknow(i+1，j)，Tunknow(i-1，j)式中：D-x1T(i，j)、D+x1T(i，j)表达式与一阶（式（10）相同；D-x2T(i，j)、D+x2T(i，j)也由链式求导法得到 D-x2T(i，j)=T0(i，j)3 T1(i，j)-4 T1(i-1，j)+T1(i-2，j)h+P(i，j)T1(i，j)D+x2T(i，j)=T0(i，j)-3 T1(i，j)+4 T1(i+1，j)-T1(i+2，j)h+P(i，j)T1(i，j)（12）（9）（11）860第 58 卷 第 4 期张云，等：2D/3D起伏地表下多震相地震波走时的因式分解程函方程算法式中各变量物理意义与一阶差分（式（9）类似。z方向差分算子的推导，与x方向类似。至此，再结合窄带技术21，便可实现直达波走时计算。2.2不规则网格求解T1为提高走时反演成像的分辨率，采用多震相走时联合反演是一个不错的选择，这就要求正演算法具有计算多震相地震波走时的能力。分区多步计算技术自提出以来，就与多种网格类地震波走时计算方法相结合，用于模拟多震相地震波传播。本算法同样作为网格单元类算法，也可结合分区多步计算技术实现多震相地震波走时的计算。然而实际地表（或地下界面）多呈不规则性，若用规则网格剖分，则地下界面（或地表）附近会产生不规则网格（图2b中M1M2界面）。这样，上述规则网格的迎风差分公式显然已不适用，因此需要针对界面附近不规则网格，推导建立新的局部走时计算公式。图 2不规则网格中的点可划分为两类：一类是界面附近点（C点）；另一类是界面点（M点）。下面分别针对这两类点建立局部走时计算公式。（1）C点为当前待计算点此时以C点为中心的z方向出现CM h，即在C点的x轴方向与z轴上半轴方向仍可采用常规迎风差分格式，而 z轴下半轴常规迎风差分格式显然已不能适用，为此建立C点z轴向下方向的不等距一阶、二阶迎风差分公式（记|CM|=h1）D-z1T(C)=T0(C)T1(C)-T1(M)h1+P(C)T1(C)D-z2T(C)=T0(C)h21 T1(D)-h2 T1(M)+(h2-h21)T1(C)h2 h1-h21 h+P(C)T1(C)（13）而二阶差分算子的推导过程如下，利用泰勒函数将T1(M)与T1(D)在C点附近展开T1(M)=T1(C)-h1T 1(C)+h21T1(C)+O(C)T1(D)=T1(C)-hT 1(C)+h2T1(C)+O(C)（14）式中：T 1(C)、T1(C)分别为T1(C)的一阶与二阶导数；O（C）为高阶无穷小量。消掉T(C)项，则有T 1(C)=h21 T1(D)-h2 T1(M)+(h2-h21)T1(C)h2 h1-h21 h（15）将式（15）代入式（14），可得不等距二阶迎风差分公式。此外，由于|CM|BC，导致T(M)与T(B)的大小已无可比性，则规则网格的迎风差分条件在此已不适用，因此建立如下不等距迎风差分格式 T(C)=minT0(C)Tup1(C)，T0(C)Tdown1(C)当 Tknow(B)，Tknow(M)T0(C)Tup1(C)当 Tknow(B)，Tunknow(M)T0(C)Tdown1(C)当 Tknow(M)，Tunknow(B)（16）式中Tup1和Tdown1分别为将 z 轴向上或向下的差分算子与x方向差分算子相结合，计算得到的走时扰动因子。（2）M点为当前待计算点此时M点周围的四个网格均为非正交，导致程函方程不成立。Rawlinsnon等31 采用三角网格缝合界面，而后在三角网格中求解程函方程，但实现较为困难。本文综合不等距迎风差分公式、费马原理与惠更斯原理，求解界面走时。对于M点的z轴下方向，同样有不等距一阶、二阶差分公式图 2模型界面及网格剖分示意图（a）及其局部放大图（b）红色实线为界面，蓝色矩形块对应放大区域。861石 油 地 球 物 理 勘 探2023 年 D-z1T(M)=T0(M)T1(M)-T1()Dh2+P(M)T1()MD-z2T(M)=T0(M)h22-()h2+h2T1(M)-()h2+h2T1(D)+h22T1(E)h22()h2+h-()h2+h2h2+P(M)T1()M（17）式中h2为M、D之间的距离。对于M点z轴上方向，表达式类似。此外对于M点x方向，由于界面起伏，x方向常规差分公式已不再适用。本文采用惠更斯原理，即将M1、M2视为子波源，就有T(M)=T(Mi)+lis MiM (i=1，2)（18）式中：li为Mi与M之间的距离；s MiM为Mi与M之间的波慢度。最后，再利用费马原理，得到 T(M)=minT(Mi)+Tup1(M)，T0(M)Tdown1(M)（19）式中Tup1(M)、Tdown1分别为M点z轴上方向或下方向采用一阶或二阶差分算子计算得到的走时扰动因子。综上所述，由于起伏界面的存在，可将界面上的点分为三类：第一类为规则网格中节点，采用常规迎风差分公式（式（8）式（10）计算修正因子T1，然后根据式（2）计算该点地震波走时；第二类为界面附近节点（C），对于该点远离界面方向和x方向，走时扰动因子的计算方法与规则网格相同，对于该点的界面方向，采用式（13）计算该方向的走时扰动因子，最后根据费马原理，由式（16）计算该点走时；第三类为界面上点（M），对于该点z轴方向，建立不等距迎风差分公式（式（17）计算该点在 z方向的走时扰动因子，最后结合费马原理、惠更斯原理，由式（19）计算 M 点地震波走时。3多震相地震波走时计算方法基于前面构建的起伏地表模型下规则网格与不规则网格的迎风差分格式，再结合分区多步技术与窄带技术，便可实现多震相地震波走时计算。所谓分区多步计算技术即指按照地下起伏界面的具体情况，将模型分成独立的层状（或起伏层状）区域，相邻区域由地下界面相连接18（图 2a）。窄带技术根据模型中网格节点所处状态将所有节点分为三种：“完成点”、“窄带点”及“远离点”。其中“完成点”表示走时计算已结束的网格点；“窄带点”表示当前波前上的网格节点，其走时大小还可更新；“远离点”表示走时计算还未实施的网格节点。表1给出了本文所提FMMF算法计算直达波走时的实现框图，当要计算其他震相的地震波走时时，只需将界面上走时最小的点作为新的子震源，而后与计算直达波相同，计算当前所需震相的地震波走时即可。现以图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）步时，自界面上走时最小的点向透射区域进行波前扩展即可；若要追踪转换波，则自离散界面上走时最小的点开始向转换波所在区域进行波前扩展扫描时时，采用不同的速度模型即可。简而言之，分区多步计算的原理就是根据所要计算地震波的类型，设计要调用的独862第 58 卷 第 4 期张云，等：2D/3D起伏地表下多震相地震波走时的因式分解程函方程算法立区域的先后顺序及对应的速度模型，达到运用简单组合完成复杂模型中多震相地震波走时计算之目的。表 1 FMMF算法实施过程初始化T(x)=x ，为计算区域内所有网格节点T0(x0)=0，T1(x0)=S(x0)完成点：，波前点：x0，远离点：，x x0当波前点 ,do（1）寻找波前中走时最小的点：xmin=min(T(x)，x front)（2）将xmin点从波前中移除，并将其状态更新为完成点：（3）将与xmin相邻的远离点加入波前点；（4）判断xmin点所在网格是否为规则网格：1）若xmin点所在网格为规则网格，则将式（9）、式（10）（一阶差分）或式（11）、式（12）（二阶差分）代入式（8），求解xmin周边网格点走时；2）若xmin点所在网格为非规则网格，则判断当前计算点是否为界面上的点：（a）若当前计算点不是界面上的点，则采用式（16）计算当前点走时；（b）若当前计算点为界面上的点，则采用式（19）计算当前点走时；结束4 数值模拟及结果分析采用几组模型进行数值模型计算，通过本算法与Multistage FMM 算法31 详尽（主要是计算精度与计算效率）对比，验证本文算法的优越性。当选用模型存在解析解时，将以解析解作为参考解；当选用模型不存在解析解时，将以发展较成熟的不规则最短路径算法（ISMP）37 作为参考解，且通过加密ISMP算法网格间距，确保该算法计算结果的精度。计算所用处理器配置为：Intel（R）Core（TM）i710750H CPU 2.60 GHz。另外，在分析算法所用CPU时间时，均是设置程序循环执行 100 次后取均值得到每次所用CPU时间，以避免偶然性。4.1模型 1（均匀模型）2D均匀模型大小为 4 km4 km，地震波传播速度为 1 km/s，震源位于（2 km，2 km）；3D 均匀模型大小 为 10 km10 km10 km，地 震 波 传 播 速 度 为2 km/s，震源位于（5 km，5 km，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。本文选用平均绝对误 差 与 平 均 相 对 误 差 两 个 指 标 评 估 计 算 精 度，平 均 绝 对误差计算公式为：()i=1N|Texact-TcalN；同 时，平 均 相 对 误 差 的 计 算 公 式 为：()i=1N|Texact-Tcal Texact 100%N。其中Texact为解析解，Tcal为通过 FMM 或 FMMF所计算的结果，N为模型中网格节点总数。表2（2D）和表3（3D）分别给出均匀模型不同网格间距剖分下FMM与FMMF的计算结果对比。从该表可见，对于均匀模型，本文所提FMMF计算结果误差为零（注意：误差分析结果已消除因计算机自身小数保留位数而导致的计算误差）；同时，可看到随着网表 22D均匀模型中 FMM 与 FMMF计算地震波走时误差及 CPU时间对比网格间距km1/101/201/501/1001/200节点总数41418181201201401401801801阶数1212121212FMM平均绝对误差/ms26.9415.41818.9142.91217.7222.3901011.7411.1641017.6931015.970102平均相对误差/%2.4460.5471.7710.0750.7510.0440.2490.0230.1300.011CPU 时间/s2.911022.911020.0940.0950.3780.3851.4891.5005.9546.143FMMF平均绝对误差/ms0000000000平均相对误差/%0000000000CPU 时间/s6.141048.191042.911033.141031.201021.521029.671020.1950.2750.309863石 油 地 球 物 理 勘 探2023 年格间距的减小，两种算法的计算效率都降低了，但在相同的网格剖分间距下，FMMF计算效率明显高于FMM，考虑是由于 FMM 采用震源网格加密方式保证精度，从而使得FMM计算效率降低。另外，当采用二阶差分时，两种算法计算效率都降低了，但降低并不明显。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 Texact=1arccosh 1+12v02v(x)x-x0（20）式中：为速度增加梯度；v0为震源处地震波传播速度；v(x)为x点的速度；x0为震源位置坐标。表 4和表 5分别给出不同网格间距剖分计算 2D和3D线性模型时FMM与FMMF的计算结果对比。同样采用震源网格加密方式保证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和表5可见，随着网格间距的不断减小，不论是 FMM 还是 FMMF，其平均绝对误差与平均相对误差都不断减小，计算精度不断提高。但同时，网格间距的减小，导致网格节点数变多，计算效率也随之下降。另外还发现，当采用二阶差分时，两种算法的计算精度都有大幅度提高，且与均匀模型相同，计算效率降低并不明显，因此建议实际应用中采用二阶差分格式。重要的是在两种算法对比中发现，网格间距相同时，2D 或 3D 情形下，无论是计算精度还是CPU计算时间，FMMF优于 FMM。这是因为计算精度方面，FMMF成功地避免了震源点的奇异性问题；而计算效率方面则是传统FMM为了克服震源点表 33D均匀模型中 FMM 与 FMMF计算地震波走时误差及 CPU时间对比网格间距km1/51/101/201/401/50节点总数515151101101101201201201401401401501501501阶数1212121212FMM平均绝对误差/ms)72.3636.32536.2563.02530.9131.6329.4630.4372.9270.304平均相对误差/%3.1090.3461.6370.1651.2956.5271020.2622.9211020.1431.532102CPU 时间/s1.2501.32817.82819.50084.89087.187687.312715.9531389.2031452.703FMMF平均绝对误差/ms0000000000平均相对误差/%0000000000CPU 时间/s9.3751029.3751021.0311.35912.14014.359177.515188.765393.656433.062图 3速度线性增加的 2D（a）和 3D（b）模型864第 58 卷 第 4 期张云，等：2D/3D起伏地表下多震相地震波走时的因式分解程函方程算法的奇异性需对震源附近网格细化，导致计算量增加。为进一步分析误差分布情况，又给出当 h=0.05 m 时两种算法分别采用一阶（图 4（2D）、图 6（3D）、二阶（图5（2D）、图7（3D）差分所得的绝对误差和相对误差分布图。可发现传统FMM计算所得地震波走时在震源对角线附近区域误差较大，当采用一阶差分格式时，其最大绝对误差分别可达 4.80 ms（2D）和 20 ms（3D），最大相对误差则分别高达 7%（2D）和 4%（3D）；而 FMMF 最大绝对误差仅分别为 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）；而 FMMF最大绝对误差则分别进一步减至0.04 ms（2D）和0.1 ms（3D），最大相对误差则分别减至 0.009%（2D）和 0.02%（3D）。可见，本文的FMMF不论采用一阶还是二阶差分格式，其计算精度高于传统 FMM；同时，通过误差分布图可发现FMMF完全解决了震源点的奇异性问题。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表 5图 3b模型中 FMM 与 FMMF计算地震波走时误差及 CPU时间对比网格间距km1/51/101/201/251/50节点总数515151101101101201201201251251251501501501阶数1212121212FMM平均绝对误差/ms65.2585.38235.3552.81729.2371.52115.431.3582.0920.261平均相对误差/%2.6730.2221.4590.1221.2626.4191020.64