自1976年Aki和Lee将医学层析成像技术应用于地震层析成像研究以来，地震层析成像技术经过近40年的不断发展和完善，已成为研究地球内部结构、以及介质各向异性行之有效的方法技术之一.其成像尺度涵盖范围很广，小到岩石组分，大至全球构造.对于区域或全球层析成像而言，地球曲率不可忽略，地球表面不能简单视作平面来处理(Snoke and Lahr，2001)，而必须采用球坐标系才能保证所需的计算精度.同时，地球内部的结构是复杂变化的，需用三维球坐标系下的速度模型来描述.在速度异常(-2%)的南太平洋地区，Zhao和Lei(2004)分别采用3D速度模型和1D速度模型计算走时和射线路径，结果显示PcP震相的射线路径改变了77 km，其他震相的射线路径偏差也达到了几十千米，走时也存在数秒的差异.因此，如果所研究区域达到一定的尺度时(例如区域、特别是全球走时成像)，不考虑地球曲率的影响和介质的横向非均匀性，计算的射线路径和走时就会包含较大的误差.

目前应用较广的走时层析成像，高频假设等因素的制约导致其成像的空间分辨率较低.因此，在现有地震和地震台网分布情况下，采用多震相走时成像是提高分辨率非常有效的途径.多震相走时成像较为典型的一个例子是Zhao等(2005)利用美国L and ers地震余震进行多震相(S，SmS，sSmS)走时同时反演S波速度和Moho面，其中仅采用相距200 km的两个观测台站和其间发生的180个余震，其结果表明使用多震相走时成像其空间分辨率可提高至台站间距的1/8～1/6(即25~35 km).

在多震相走时成像中，由于反射或反射转换波震相的走时由震源参数、3D速度分布和反射界面共同决定，且这三类参数反演前往往是未知的，并存在强耦合关系(如：震源位置较小的变化将引起射线路径、以及射线入射到界面上反射点位置的较大变化，反过来说，反射点位置较小的变化也将同样引起射线路径的变化等)，因此要求反演算法应具有同时反演三参数(速度、震源参数、反射界面)的能力.

对于区域或者全球走时层析成像，为保证计算精度和提高成像分辨率，需要进行三维球坐标系下多震相走时三参数同时反演成像.由于全球走时成像中需要成千上万的走时数据(大多采用不仅仅是单一震相的走时)，且缺少高效和高精度的三维正演射线追踪计算方法及相应软件，故目前采用较多的还是基于1D球坐标系下的两点射线追踪算法，此算法反演中假设射线路径不变，然后采用扰动的方法更新速度模型( Van der Hilst et al.，1997; Su et al.，1994; Grand et al.，1997).近些年来人们已意识到采用三维球坐标系的重要性，相继出现了几种不同的射线追踪方法，例如：Bijwaard和Spakman(1999)研发的快速运动学射线追踪算法，在网格算法中引入了弯曲法，能够追踪3D球坐标系下的(临界)折射、多次反射、以及一些绕射(Pdiff)波，且计算的各主要震相绝对走时误差约为0.1 s，但实质上该算法依然是两点间的射线追踪.Zhao和Lei(2004)将球坐标系下的伪弯曲射线追踪算法(Zhao et al.，1992)进行了推广，可进行区域或全球多震相走时反演成像.Tian等(2007)开发了一套可用于全球有限频成像的动力学射线追踪程序，所计算的主要震相走时绝对误差约为0.1 s，但其算法依然是基于两点间的射线追踪问题(局域解问题)，计算速度仍然存在问题(每个震源和台站都需要重新进行射线追踪).

网格/单元波前扩展算法由于其具有全局解、算法稳健、计算效率高等优势越来越受到人们的关注，其中主要包括有限差分解程函方程算法(Vidale，1988，1990)和最短路径算法(Moser，1991; Klimeš and KvasniČka，1994). Rawlinson和Sambridge(2004)采用有限差分解程函方程的快速行进法与分区多步计算技术相结合，实现了直角坐标系下层状介质中多次波的追踪计算.后来，De Kool等人(2006)将其推广到3D球坐标系下，能够快速、准确、稳健地计算复杂区域模型中的多种震相，包括后续反射波，但该算法目前为止未见用于全球走时成像的研究报道.白超英等则采用改进型最短路径算法结合分区多步计算技术实现了直角坐标系下规则/不规则最短路径算法的多次波的追踪(唐小平和白超英，2009a，2009b; Bai et al.，2010; 赵瑞和白超英，2010)，研究结果表明，无论是计算精度，还是CPU 时间，改进型最短路径算法均优于解程函方程的快速行进算法(赵瑞和白超英，2010).随后又将该算法应用于直角坐标系下2D/3D地震多震相走时联合双参数(黄国娇和白超英，2010；白超英等，2011)同时反演成像.

目前为止，走时三参数(速度、反射界面和震源参数)同时反演研究并不多见，Sambridge(1990)虽然是首位讨论走时三参数同时反演的研究学者，反演采用子空间反演方法(Kennett et al.，1988)，但其正演采用了传统的边界值两点射线追踪算法(Sambridge and Kennett，1990)；Preston(2003)在他的博士论文中也进行了走时三参数同时反演的研究工作，其研究成果已在Science上发表(Preston et al.，2003)，其中利用正则化加权处理方式对三种参数进行解耦，容易引入人为因素.在直角坐标系下，黄国娇和白超英(2013)采用分区多步最短路径算法(Multiple ISPM)进行多震相的射线追踪，反演采用子空间算法，实现了3D多震相走时三参数同时反演成像.

鉴于此，针对区域或全球层析成像，本文将上述分区多步最短路径算法推广至球坐标系中，进行多震相走时和射线路径的追踪计算.然后将其与地震多参数反演算法(子空间反演算法)相结合，提出了一种3D球坐标系下可利用地震多震相走时资料进行三参数同时反演的方法技术.为验证该三参数同时反演算法的有效性和正确性，我们选用一个区域模型进行数值模拟实验，在合成理论走时数据中加入随机噪声并对其进行三参数同时反演，将反演成像结果与双参数同时反演进行对比分析.结果表明，三参数和双参数同时反演方法技术具有相近的反演能力，都能较好地反演出震源位置、速度分布、界面形态和位置.

利用分区多步不规则最短路径算法进行射线追踪之前，需要将模型参数化.三维球坐标系下模型参数化如图 1所示，首先根据速度间断面的起伏形态将模型分为四个区(图 1a)，然后对每个区域用规则单元或者不规则单元将其参数化.地心点所在单元形成如图 1c所示的五面体，图 1d为两极极点所在单元，而不包含地心、极点且远离起伏地表和界面的区域则可采用如图 1b所示的规则六面体单元参数化.单元的角点为主节点，主节点间等距离的插入次级节点(这点很重要，因为越往地心，相同扇形的单元弧长越短)，次级节点只分布在网格单元的边上，单元内无节点，但震源和地震台站可以位于其中.

图 1

Fig. 1

图 1 球坐标系下模型参数化示意图(图a, 为清晰起见, 次级节点未标出)和所用参数化单元(图b—d), 图b—d中大球表示主节点, 小球表示所加入的次级节点 Fig. 1 Model parametrization in 3D spherical coordinates (diagram a, for clear representation, the secondary nodes are not depicted) and three kinds of cells used to divide the spherical model (diagram b, c and d, in the figure the larger circles are primary nodes and the smaller circles are the secondary nodes, respectively)

速度采样在主节点上进行，次级节点的速度可由主节点上的速度值通过三次线性插值函数(3D)得到.当节点(震源或者地震台站)位于单元内部时，需分规则单元和非规则单元两种情况对待：

(1)当节点位于规则单元中，如图 1b所示，节点速度则由所在单元八个主节点上的速度值通过三线性拉格朗日插值函数给出，插值函数为：

x_k^e和v(x_k^e)(k=1，2，…，8)分别为立方体单元八个角点上的矢量坐标及速度.

(2)若节点位于不规则单元内，如图 1c—1d所示，则需先将单元剖分为几个四面体(四面体的四个角点必须是主节点，编号1—4)，然后再由节点所在四面体的四个主节点上定义的四面体体积插值函数给出(Huang et al.，2013)：

式中u_j为点在四面体的体积坐标，v(x_j^e)为主节点上的采样速度值.值得注意的是，本文采用两个不同的速度值对速度不连续面(反射界面)进行定义，即界面处存在两个速度，界面上速度和界面下速度.就像在AK135模型一样，当计算界面的反射波时，就采用界面上速度，而当计算透射或者折射时，则需用界面下速度.

分区多步计算首先是模型分区，然后按照给定的震相，自震源所在区域开始，沿着目标射线通过的区域分段分区进行射线追踪计算.举例来说，如图 1a中所示，根据反射界面将模型分为了四个区，速度不连续(间断)面连接相邻区域.若计算图 1a中界面一(深度660 km)的透射或反射震相时，首先自震源所在的单元开始波前扩展计算，然后按照一定的步长进行第一区的波前扩展计算.等第一区所有节点上的走时计算完成之后，波前停留在第一区的速度离散界面上，并保存了震源到界面离散点上的最小走时及相应的射线路径.此时如果计算透射波P₁P₂(上、下标分别代表上行波和下行波，1和2表示区号)，则调用第二区的P波速度(如果追踪P₁S₂震相则调用第二区的S波速度)，自界面一上走时最小的点开始进行第二区向下的波前扩展扫描计算；若是计算反射波P₁P¹(或反射转换波P₁S¹)，同样自界面一上走时最小的点开始，并调用第一区相应的速度模型向上进行波前扩展计算.同理，依据多次波是由简单的入射、透射、反射及转换波按照一定的规律组合而成的原理，即可实现任意多震相地震射线的追踪计算.图 2是3D球坐标系下采用上述分区多步不规则最短路径算法追踪计算的多震相射线路径图(所用速度模型和反射界面位置详见第4节反演实例)，更复杂的多次波均可采用该算法进行追踪计算，但是需要在输入文件中给出多次波的传播途径，即射线穿过哪些区域，与哪些界面相交，遇到界面是反射、折射还是转换.

图 2

Fig. 2

图 2 球坐标系下多震相射线路径图 浅灰色线：初至P波；灰线：深度410 km处界面的反射波P410P；黑线：深度660 km处界面的反射波P660P. Fig. 2 Multiple ray paths in a spherical coordinate The directed arrivals (light grey lines), the reflection from the 410 km subsurface (gray lines) and the reflections from 660 km subsurface (black lines) are shown.

Huang等(2013)将该算法计算的49种不同的全球主要震相走时与AK135全球走时表(Kennet et al.，1995)进行对比分析，讨论该算法的计算精度问题.结果表明绝大多数震相的最大绝对走时误差小于0.1 s.在全球或者区域走时层析成像中，还要求正演算法具有较高的计算效率.因此，我们分析了图 2所示模型中三种不同震相(P，P₄₁₀P和P₆₆₀P)在不同尺度单元参数化和次级节点情况下的计算效率，结果如图 3所示.其中震源和台站的个数及位置分布均与第4节的反演实例中相同，单元大小分别为4°×4°×200 km(图 3a)、2°×2°×100 km(图 3b)和1°×1°×50 km(图 3c)，计算机配置为Intel i5-2300CPU处理器、2.8 GHz的主频和4 G内存.

图 3

Fig. 3

图 3 计算多震相(P，P410P和P660P)的CPU时间图 (a) 4°×4°×200 km；(b) 2°×2°×100 km；(c) 1°×1°×50 km.图中括号内数字为当前模型参数化下节点总数. Fig. 3Results of computational efficiency tests on tracing the multiple phases ( P, P410P and P660P) (a) Cell scale of 4°×4°×200 km; (b) Cell scale of 2°×2°×100 km; (c) Cell scale of 1°×1°×50 km. In the figures, the numbers in the parentheses indicate the total number of nodes for the corresponding parameterized model.

从图 3c中可以看出，当单元大小为1°×1°×50 km，加入11个次级节点(地表次级节点间隔约为8 km)时，计算这三种震相几乎要花费近8个小时.但是，通常加入5~7个次级节点(次级节点间距为17~22.5 km)就能满足精度需求，此时CPU时间就减少至1.0~1.5个小时.值得注意的是，反演成像中每次迭代基本上有95%的时间花费在计算这三种震相的走时和射线路径上.因此，我们可大致估计采用上述三种震相同时反演三参数的时间.

三参数的同时反演是一项很困难的工作，因为三者具有不同的物理量纲，并且具有强耦合关系.若在同一个方程中求解这三种不同的模型参数，解会偏向系数大的一方，很可能导致反演的失败.因此，如何解耦成为三参数同时反演中急需解决的问题.直角坐标系中，我们将共轭梯度法求解带约束的阻 尼最小二乘算法(DMNCL-CG)和子空间法(Subspace)进行了对比，结果表明Subspace更适合于多参数的同时反演.因此，本文采用Subspace算法对球坐标系中三参数进行同时反演.

多震相走时同时反演可归结为求下列目标函数的极小值问题(Rawlinson and Sambridge，2003)：

其中g(m)和d_obs分别代表理论走时和观测走时； m和m₀分别是更新的模型和初始模型； C_d^－1和C_m^－1分别为数据协方差和模型协方差矩阵的逆.ε为阻尼因子，η是光滑因子，D 为二阶导数差分算子.阻尼因子ε和光滑因子η用来平衡模型估计 m 计算的理论走时与观测走时间的拟合度、模型估计 m 与初始模型 m ₀的接近程度以及 m 的光滑程度.

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

式中为步长因子，和分别代表梯度向量矩阵(一阶偏导数)和Hessian矩阵向量(二阶偏导数).由公式(3)可以推导出一、二阶导数的表达式：

其中Frechet偏导数矩阵.将表达式(5a)(5b)代入公式(4)中就可得到反演每次迭代时模型的更新量：

公式(6)可直接用于每次迭代中计算模型的更新量，但是每次都需要重新计算子空间基向量A、梯度向量和Frechet偏导数矩阵 G .其中最关键的问题是计算Frechet偏导数矩阵和构建基向量，Frechet偏导数的计算将在接下来的部分详细介绍，而大多数子空间算法构建的基向量{a^j}均是利用模型空间中的最速上升向量 γ = C_m、以及它的变化梯度(Rawlinson and Sambridge，2003)，本文也采用此方法.对于给定的子空间维数p(p必须大于等于待反演的参数种类)，构造子空间基向量{a^j}(j=1，2，…，p)的具体步骤如下：

第一步：构建基向量 a¹=C_m

第二步：若p≥1，则采用下列公式构建其余的基向量：

第三步：采用奇异值分解法(SVD)将第一步和第二步中构建的基向量{a^j} 标准正交化.

三维球坐标系中，地球内任一点P可表示为(r，θ，φ)，实际应用中通常用径向深度、纬度和经度来表示，记为(z，Lat，Long)，并且有r=R-z(R=6371.5 km为地球半径)；θ=90°-Lat和φ=Long.根据射线理论，由第i个震源经反射点到达第j个检波器的走时可表述为下列积分：

其中ds为沿路径R_ij上的线积分元，而V_c(r，θ，φ)则为沿射线路径上的速度分布函数.当模型网格单元化后，沿该射线路径R_ij上的走时可写成求和形式：

其中：n为射线穿过的单元总数，R_ij，k、V_c，k(r，θ，φ)分别为射线穿过第k个单元的长度和速度分布.

在速度、界面和震源位置同时反演时，Frechet偏导矩阵包含了三部分：一是走时关于速度变化的导数，二是走时关于反射点深度变化的导数，三是走时关于震源参数的导数.

其中v_k和r_k分别是第k个单元内的速度分布和第k个反射点的深度，Δv_k和Δr_k则分别是相应速度和反射点深度的扰动量，Δx_k(k=1，2，3)是第k个震源分别在r，θ和φ方向的扰动量.

(9)式中第一项，即有关走时关于速度的一阶偏导数计算如下(为简便起见，这里略去震源编号i，则t_j是第j条射线R_j的走时)：

当射线穿过某个单元时，走时对速度的一阶导 数可用射线两端点导数的平均值代替(Bai and Greenhalgh，2005):

其中：，k_i，k_j是穿过第k单元射线段的两个端点；BV(k)则是该单元上主节点的速度值.

(9)式中的第二项，即：走时关于反射点深度变化的偏导数可由以下步骤计算：若考虑一般情况(反射、折射、以及反射转换)，则走时关于反射点深度变 化的导数由(11)式确定(Rawlinson and Sambridge，2003)：

其中r_n是离散界面节点的深度坐标，r_int是反射点的深度坐标(有关w_j，w_j+1，w _n和w _r，见图 4所示)，球坐标系中 w _r=(1，0，0)，根据界面节点上的深度插值函数计算.v_j和v_j+1分别是界面两侧的速度值.当仅考虑反射时，w _j+1· w _n=－ w _j· w _n和v_j=v_j+1，则(16)式简化为：

图 4

Fig. 4

图 4 走时关于反射点深度偏导数计算示意图 Fig. 4 Diagrammatic showing the calculation of traveltime derivative with respect to reflector depth change

(9)式的第三项，即走时关于震源参数(位置和发震时刻)变化的偏导数可由以下步骤计算，首先走时(t=t_a－t_o)关于发震时刻的导数为：

而走时关于震源深度变化的偏导数的计算可参见图 5所示.假设震源的深度扰动值为Δr_h，则走时关于震源深度变化的偏导数可表示为：

图 5

Fig. 5

图 5 走时关于震源位置变化偏导数计算示意图 Fig. 5Diagrammatic showing the calculation of traveltime derivative with respect to source location change

同理可得走时关于震源位置其余两个方向(Lat，Long)变化的偏导数表达式：

其中α_h如图 5所示，代表射线与r轴的夹角；β_h和γ_h则是射线在θφ 平面投影后分别与θ和φ轴的夹角.

为检验该三参数同时反演算法的有效性及正确性，我们选择一个尺度为20°×20°×700 km的区域模型，如图 6所示.图 6a是经度-深度剖面的起伏速度图，速度沿着纬度方向不变化.模型中引入两个起伏反射界面模拟地幔中两个速度不连续面(如图 6b所示)，深度分别位于410 km和660 km处，起伏幅度均为±20 km.模型参数化选用1°×1°×50 km的单元，次级节点间距为10 km(相当于在单元面的不同方向加了9个次级节点).30个震源随机分布在深度20 km至100 km之间(图 7a是震源在经纬度面上的投影；图 7b则是震源在经度-深度剖面上的投影)，441个检波器均匀(间距为1°×1°)分布于地表(模型顶面).

图 6

Fig. 6

图 6 (a)某一纬度时, 经度-深度剖面速度图; (b)两个起伏反射界面图 Fig. 6 (a) Velocity model showing one cross section and (b) two undulating subsurface interfaces

图 7

Fig. 7

图 7 震源在经纬度剖面(a)和经度-深度剖面(b)投影 Fig. 7 The sources projected on the horizontal plane (a) and the longitude-depth cross section (b)

本文的同时反演中采用了3种不同的震相：直达P波、反射界面410 km和660 km的一次反射波(P₄₁₀P、P₆₆₀P).为检验反演算法抗噪声能力，三种震相走时数据中加入了随机噪声(直达P波、P₄₁₀P和P₆₆₀P分别加入±0.2 s、±0.3 s和±0.3 s的随机误差).选用背景层状速度作为初始速度模型，两个 初始界面均为水平(深度分别为410 km和660 km). 为了模拟实际地震定位中存在误差这一现实，初始震源位置的选择则是分别对实际震源位置在三个方向(经度、纬度、深度)进行扰动后得到，其相应扰动最大幅度分别为±2°、±2°和±10 km.

为了进行系统的对比研究，模拟实验中进行了三种不同的反演，其一是速度和界面的同时反演(简记为：V+R反演)，其二是速度和震源位置的同时反演(简记为：V+S反演)，其三是速度、界面和震源位置的同时反演(简记为：V+R+S反演).不失一般性，三种不同的反演选用相同的反演参数，即阻尼因子和子空间维数分别为0.001和6(对2~8之间逐一尝试，选择较合适的子空间维数为6).迭代均为15次，相应的走时残差分别由初始时的 1.834 s，9.15 s，9.269 s下降至0.087 s(V+R反演)，0.25 s(V+S反演)，0.289 s(V+R+S反演).

图 8给出了三种不同反演算法在Lat=10°垂直剖面上速度场的反演结果.从中可以看出，三种不同的反演算法重建的速度场分布与真实速度(图 8a)具有相似的异常形态，但是异常起伏的幅度似乎未完全恢复，且模型下部的更新相对较差，其主要原因是模型下部只有P₆₆₀P震相的射线穿过，射线覆盖率和交叉概率有限所致.对比图 8b—8d，会发现相对于双参数的反演结果(图 8b和8c)，三参数同时反演的结果(图 8d)要略差些，但是差异不明显，仍然能够反演出异常的形态和结构.

图 8

Fig. 8

图 8 三种反演算法在Lat=10°垂直剖面的成像结果 (a) 真实模型; (b) V+R; (c) V+S; (d) V+R+S. Fig. 8 The reconstructed velocity fields in Lat=10 ° cross section with three kinds of inversion methods (a) Real velocity model; (b) V+R; (c) V+S; (d) V+R+S.

为更加直观显示反演更新的速度场与真实速度场之间的差异，图 9给出了与图 8对应的垂直剖面上反演的速度场与真实速度场之间的百分比差异.从图 9中可以看出，虽然三种反演方法都能重建出速度异常的形态和结构，但是速度异常幅度并未完全恢复.并且，可以明显地看出双参数同时反演结果要略好于三参数同时反演结果，特别是模型下部速度异常的重建结果.

图 9

Fig. 9

图 9 三种不同反演算法的反演结果在Lat=10°垂直剖面上与真实速度的百分比误差 (a)初始速度模型; (b) V+R反演结果; (c) V+S反演结果; (d) V+R+S反演结果. Fig. 9 The difference (in percentage) between the initial and true velocity model in Lat=10° cross section, between the true velocity and the inverted velocity with three different inversion methods in the same cross section (a) Initial velocity model; (b) Results for the V+R inversion; (c) Results for the V+S inversion; (d) Results for the V+R+S inversion.

图 10为两种反演算法对反射界面的更新结果.可以看出，两种反演算法都能较好地更新反射界面的形态和位置，与真实反射界面(图 6b)很相似，但仅从图 10中无法判断哪种算法的更新结果更好.因此，为进一步评价反射界面的更新结果，我们给出了两种反演算法更新的反射界面与真实反射界面间的平均差异(表 1).由于三参数(V+R+S)同时反演的难度远大于双参数(V+R)同时反演的难度，故从表 1中可知，双参数(V+R)同时更新反射界面的结果要略好于三参数(V+R+S)同时更新的结果，但通过图 10的对比可得出三参数也能较好地更新反射界面几何形态和位置.

图 10

Fig. 10

图 10 两种不同反演算法更新的反射界面 (a) V+R反演结果; (b) V+R+S反演结果. Fig. 10 The updated reflected interfaces by the two different inversion algorithms (a) Results for the V+R inversion; (b) Results for the V+R+S inversion.

表 1(Table 1)

表 1 界面深度平均绝对误差统计表 Table 1 Mean absolute error statistics of the inverted depth for the two interfaces 平均绝对误差 (km) 初始值 V+R反演 V+R+S反演 410 km界面 12.76 1.55 4.03 660 km界面 12.76 1.42 3.39

表 1 界面深度平均绝对误差统计表 Table 1 Mean absolute error statistics of the inverted depth for the two interfaces

图 11是两种不同反演算法对震源位置的更新结果.这里我们用三种收敛情况来评价震源位置的更新结果，更新的震源位置相对于真实震源位置距离的收敛情况(图 11左)、更新震源的水平位置相对于真实震源水平位置距离的收敛情况(图 11中)和更新的震源深度相对于真实震源深度的收敛情况(图 11右).从图 11中可以看出，双参数的反演结果(图 11d，11e，11f)要略好于三参数的反演结果(图 11g，11h，11i).无论是双参数(V+S)的同时反演结果(图 11d—11f)，还是三参数(V+R+S)的同时反演结果(图 11g—11i)，其中对震源水平位置(图 11e和11h)的更新结果均明显好于震源深度(图 11f和11i)的结果.可知在该模型下存在震源深度的定位误差大于水平面内的定位误差的问题，其主要原因是对于区域或者全球成像而言，虽然利用了反射波走时资料，但模型尺度相对较大，走时对震源深度较小的变化不敏感，从而导致震源深度的更新相对较差.然而，在初始震源偏离真实震源较远(此例中最大距离约为320 km)的情况下，三参数的同时反演中，对震源的更新结果还是较为理想的.

图 11

Fig. 11

图 11 两种不同反演算法得到震源的位置(a, d, g)、水平位置(b, e, h)和深度(c, f, i)分别相对于真实震源位置的收敛情况 (a, b, c)初始值; (d, e, f) V+S反演结果; (g, h, i) V+R+S反演结果.DZ表示在深度方向上距离真实震源的距离. Fig. 11 The convergence to the real source locations (a, d, g), the real horizontal position (b, e, h) and the real focus depths (c, f, i) by two different algorithms (a, b, c) Initial; (d, e, f) Results for the V+S inversion; (g, h, i) Results for the V+R+S inversion.

对于区域或全球层析成像，为保证计算精度、提高成像分辨率，本文将直角坐标系下的多震相射线追踪正演(Multistage ISPM)算法推广至球坐标系中，并与比较流行的多参数反演算法(Subspace)相结合，提出一种3D球坐标系中多震相走时三参数同时反演的方法技术.数值模拟实例表明，三参数同时反演方法技术和双参数反演方法具有相近的反演能力，都能较好地反演出震源位置、速度分布、界面形态和位置.另外，无论是双参数或是三参数同时反演算法均对走时数据中的随机误差不太敏感，这点在实际应用中也是至关重要的.因此，本文提出的球坐标系下的三参数同时反演算法是一种利用多震相走时资料同时反演速度、震源位置以及反射界面行之有效且实用性较强的方法技术.