2011, Vol. 54
2. 长安大学计算地球物理研究所,西安 710054
2. Institute of Computing Geophysics, Chang'an University, Xi'an 710054, China
自Aki和Lee将医学层析成像用于地震学研究以来[1]，经历了30多年的发展，已形成众多的地震反演方法与技术.地震层析成像技术已成为人们研究地球内部结构及各向异性行之有效的方法技术之一.其成像尺度涵盖范围很广，小到岩石组分，大至全球构造.地震层析成像按照所采用的资料类型可大体分为:体波成像、面波成像、全波型反演成像、地球自由振荡等.目前应用较多的依然是体波成像，尽管全球成像工作中使用了大量多震相走时资料[2~4]和部分区域成像使用了多震相资料[5, 6]，但多数研究者的区域成像依然是单一震相的.
随着油气矿产资源的不断开发、以及地壳及上地幔结构的深入研究，对地震成像分辨率也不断提出更高的要求.然而，在全波型反演成像技术还无法广泛应用的今天，如何提高走时成像的分辨率则成为首要问题.一般而言，单一震相走时成像的分辨率与观测台网的间距和地震分布有关，提高分辨率的途径主要有两种:其一是加密台站；其二是利用多震相走时资料.加密台站受经济及运行成本的制约短期内不太可行，因此，利用多震相走时资料就显得尤为重要.
多震相走时成像较为典型的一个例子是赵大鹏等利用美国Landers 地震余震进行的多震相(S、 SmS、sSmS)走时S波的成像，其中仅用了相距200km的两个观测台站和其间发生的180 个余震，结果表明使用多震相走时成像其空间分辨率可提高至台站间距的1/8~1/6，本例中为25~35km[7].
采用多震相走时数据进行联合反演时，若包含反射(或折射)波走时，则必须知道反射界面的位置才能反演出地下介质的速度分布；然而，地下反射界面的位置和形态往往反演前是未知的，这就要求反演算法应具备同时更新速度分布和反射界面的能力，即:进行速度和界面的同时反演.多震相中的反射波、折射波、或者透射波不仅包含了地下介质速度分布的信息，还包含了速度界面的信息，因此，利用多震相走时进行速度和界面的同时反演是必要的、也是可行的.
然而速度和界面的同时反演是一项比较困难的工作，因为界面的微小变化会引起反射(或折射)波射线路径的较大变化，加上速度和界面的强耦合关系，可能会导致反演失败.即在一个方程组中同时反演两种参数，解会向系数大的方向偏离.解决这一问题，目前主要有三种不同的方法:一是参数分离的方法，如:刘福田等人提出的正交投影算子分解法，对速度和界面参数进行解耦[8~10]；李松林等人[11]则运用Pavlis等[12]提出的参数分离法进行同时反演成像.二是Kennett等[13]提出的Subspace算法，能同时有效地反演多种模型参数.随后Williamson运用 subspace算法来确定速度和界面形状[14].三是对速度和界面参数进行加权归一，如:Zelt等人[15~22]对两种参数进行加权处理同时反演了速度和界面的几何结构.
本文将新近开发的多震相地震射线追踪正演算法与反演算法相结合，采用直达P 波、来自不同界面的一次反射波、及地表与界面间的二次反射波的走时资料联合进行同时反演成像，与单一震相走时成像相比，可有效地提高成像的空间分辨率.反演中对两种不同模型参数进行了归一化处理，同时引入了不同震相走时数据的加权系数，分别讨论了反射界面的精确定位、速度模型重建、以及速度和反射界面的同时反演问题.
2 正演算法简介Rawlinson和Sambridge等提出了分区多步计算技术，并采用有限差分解程函方程的快速行进法实现了层状介质中多次波的追踪计算[23, 24].白超英、唐小平等则采用改进型最短路径算法结合分区多步计算技术实现了2D/3D 复杂层状介质中多次波的追踪计算[25~27].随后又将该算法推广，实现了不规则最短路径算法下的多次波的追踪[28, 29]，其结果表明:无论是计算精度，还是CPU 时间，改进型最短路径算法均优于解程函方程的快速行进算法(相同模型和计算精度下，快速行进法的CPU 用时一般是改进型最短路径算法的5~6倍).鉴于此，有必要将改进型最短路径算法与反演算法相结合开发一套实用性强的多震相走时同时反演成像的方法技术.有关分区多步改进型最短路径算法的基本原理简述如下.
图 1所示为一个地表起伏、含有两个起伏反射界面的三维速度模型示意图.模型参数化时，起伏界面附近采用不规则单元(图 1b)，而其他区域则采用规则单元(图 1c).速度采样在主节点上进行，次级节点(含炮点及检波器)上的速度值则由插值函数得到[28].分区多步计算技术的基本原理中，模型“分区"是指按照模型的地表起伏和速度界面的具体情况将模型分成相对独立的层状区域，相邻区域由速度界面相连接(图 1a).“多步"计算是指按照所要追踪射线的类型，从炮点所在的区域开始，沿着目标射线所要通过的区域，逐区进行波前扫描.具体来说，自炮点所在的单元开始进行扫描，按照一定的步长进行波前扩展.等当前区域所有网格单元内节点都计算完毕后，波前停留在该区的速度离散界面上.若追踪下行透射(或透射转换)波，则从速度界面中走时最小的点开始，继续新区内的波前扩展和射线追踪；若追踪上行反射(或反射转换)波，则从速度界面中走时最小的点开始，继续原区内的波前扩展和射线追踪.若存在转换波，则调用相应的速度模型.按照上述步骤重复则可实现多次波的追踪计算[25~29].
|
图 1 模型中主节点和速度界面分布(a)，模型中不规则网格单元(b)和规则网格单元(c)(b)和(c)中黑色球表示主节点，白色球表示所加入的次级节点. Fig. 1 Distribution of primary nodes and interfaces in 3D velocity model (a)，and irregular cell (b) and regular cell (c) In diagrams (b, c) the black and white circles are primary and secondary nodes, respectively. |
图 2给出了上述模型(图 1所示)中多次波射线路径，其中图 2a显示了来自两个不同界面的反射波(P1P1 及P1P2P2P1)的射线路径(约定:上标代表上行波，下标代表下行波，数字1、2 表示区号)；图 2b则为自炮点激发的P 波上行至地表，经反射后分别下行至第一、二个反射界面，然后反射且转换成S波上行至检波器的多次反射转换波(P1P1S1 及P1P1S2P2S1)的射线路径[28].
|
图 2 (a)透射和一次反射波射线路径，(b)多次透射、反射及转换波的射线路径(a)中黑线:反射震相P1P1；灰线:反射波震相P1P2P2P1.(b)中黑线:多次反射震相P1P1S1；灰线:多次透射、转换及反射波震相P1P1S2P2S1. Fig. 2 Raypaths of (a) the transmitted and primary reflections and (b) the multiple transmissions, conversions and reflections In the figure (a) the black lines are the raypaths of P1P1 phase and grey lines are the raypaths of P1P2P2P1 phase.In the figure (b) the black lines are the raypaths of P1P1S1 phase and the grey lines are the raypaths of P1P1S2P2S1 phase. |
一般而言，非线性多震相(含反射或折射波)走时联合同时反演，可归结为下列带约束的阻尼最小二乘最优化问题，其反演问题的二范数目标函数为
|
(1) |
式中:d= [p1(δt11，δt21，…，δtM11)；p2(δt12，δt22，…，δtM22)；…；pn(δt1n，δt2n，…，δtMnn)]是n种震相实际观测数据和理论数据的走时残差向量，其中:p1，p2，…，pn为反演中考虑各震相所占的权重(或相应震相走时拾取时误差的精度类别)，M1，M2，…，Mn则为各震相的射线条数；m= [λ1(δv1，δv2，…，δvN1)；λ2(Δz1，Δz2，…，ΔzN2)]是模型空间参数的相对改变量，其中:(δv1，δv2，…，δvN1)为速度模型的相对改变量，(Δz1，Δz2，…，ΔzN2)则为离散反射界面位置的相对改变量，λ1，λ2 则为考虑两种不同量纲模型参数的归一化因子.A是Jacobi矩阵，其维数为(M1+M2+…+Mn)×(N1+N2).μ 是阻尼因子；(a，b)和(c，d)分别是反演解中速度模型、离散反射界面位置可容许的上下边界值.
Menke得到公式(1)解的一般形式为[30]
|
(2) |
其中Cm, Cd 分别为模型和数据空间的协方差矩阵.(2)式是非线性问题局部线性化的基本反演公式，其解具有局域解的特征，为了得到具有实际物理意义的解，可采用Cm、Cd 及Zm 的先验信息.因此，选择不同组合的(Cm, Cd, Zm)可得到不同形式的反演解.
若μ=0，Cd=I；则(2)式变为经典的最小二乘问题；
若Cd=Cm=I，Zm=W，A=GW，这里W是表征射线宽度的带宽矩阵，则(2)式转化成Meyerholtz和Szpakowski提出的卷积压制算法[31]；
若Cd=I，Cm=DTD，这里D是一阶或二阶差分算子，则(2)式转化成由Shaw 和Orcutt提出的一阶或二阶正则化方法[32].
本文采用Tarantola和Valette提出的广义带约束的阻尼最小二乘反演解，即:Cm, Cd 分别取模型和数据空间的协方差矩阵的逆阵，则解为[33]
|
(3) |
其解的约束为
|
上述矩阵方程的解可用共轭梯度法采用循环迭代算法求解[34, 35].目前关键的问题是如何求解具有偏导性质的Jacobi矩阵中的元素.
3.2 偏导Jacobi矩阵元计算三维情形下根据射线理论，沿某一射线路径Ri的走时tj可表述为下列积分:
|
(4) |
其中ds为沿路径Ri上的线积分元，而Vc(x，z)则为沿射线路径上的速度分布函数.当模型网格单元化后，沿该射线路径Ri上的走时可写成求和形式:
|
(5) |
其中:k为射线穿过的单元总数，Rj，k、Vc, k(x，z)分别为射线穿过第k个单元的长度和速度分布.
在界面与速度同时反演时，Jacobi偏导矩阵包含了两部分:一是走时关于速度变化的导数；二是走时关于反射点深度变化的导数.
|
(6) |
式中tj是第j条射线Rj的走时；vk是长度为N的模型空间中的第k个单元的速度分布；zk是离散界面点中第k个点的深度值.
当射线穿过某个单元时，走时对速度的一阶导数可用射线两端点导数的平均值代替[35]:
|
(7) |
其中:
|
式中:ki，kj是穿过第k单元射线段的两个端点；BV(k)则是该单元上主节点的速度值.
在反射及透射情形，三维介质中走时对界面位置p点的一阶偏导数可按周龙泉等[10]给出的计算公式.作为一种近似，当只考虑z方向情形，可得Zelt等[16]推导的反射波走时对界面深度的偏导公式:
|
(8) |
当只考虑反射时，(8)式可进一步简化为Bishop等[36]推导的反射波走时对界面反射点的偏导关系式:
|
(9) |
本文中走时对速度的一阶偏导数可以用公式(7)计算，走时关于反射点深度变化((9)式)可由射线追踪正演后利用有限差分法或三次线性插值函数求得，这里采用更为便捷的方法，仿照地震精定位中的做法[37]，即正演时在初始反射界面上下引入两条平行的辅助离散界面，这样(9)式中走时关于反射点深度变化的导数在正演过程中即可求出.即
|
(10) |
其中Δt1，Δt2，…，Δtf为上述两个辅助离散界面相对应点的走时差，而Δz1，Δz2，…，Δzf则为相应z方向的距离，f为离散反射界面上反射点的总数.
3.3 反演参数的归一化处理反演中Jacobi偏导矩阵中既包括走时关于速度的偏导数，又包括走时关于离散反射点深度的偏导数.然而两者的量纲是不同的，数值大小也不在同一量级.加之一条射线要穿过许多个单元，但只到达几个界面.因此，Jacobi偏导矩阵中走时关于速度偏导的矩阵元总数要大于走时关于界面深度偏导的矩阵元总数.这种情况下，走时残差函数就主要作用在速度参数上，而界面深度参数的改变量就会很小，或者几乎不变.为了避免这个问题，McCaughey 和 Singh[18]、James等[19]提出了几种解决办法.总的思想就是要将两种参数的一阶导数等量化.
一种做法是每次进行求解时，速度参数改变量的总和要与界面深度改变量的总和在数值上相当.采用的办法就是对每个界面深度改变量乘上因子w.w的值可通过下面的式子计算[18]:
|
(11) |
另一种做法就是将两种参数的一阶导数值归一到同一个数量级[19].其中又分成两种不同的归一算法.一种是对两种参数各自归一；另一种是不分参数类型进行归一.本文对上述各种归一算法进行了试算，最终采用反演效果较好的两种参数各自归一的方法.
3.4 共轭梯度法求解带约束的阻尼最小二乘问题对于方程(3)，选择不同算子、阻尼因子就会得到不同的解.有一些计算广义逆矩阵的算法求解非稀疏的病态矩阵效果很好，例如奇异值分解法、最小二乘法等，但是它们计算相对费时、且需占用大量计算机内存.而那些代数重建法、联合迭代重建法以及迭代的最小二乘法，由于其对噪声的敏感性和解的非惟一性导致这些方法在解非线性问题时，常常得不到比较接近真实模型的解.Scales成功应用共轭梯度法求解最小二乘问题，并证明了这种方法是简单、快速、有效，适合于稀疏矩阵的求解[38].事实上，共轭梯度法对带约束的阻尼最小二乘问题也是有效的，因为物理解的存在性，系数矩阵的条件数可通过先验信息得到改善，因此，本文采用共轭梯度求解带约束的阻尼最小二乘的迭代算法[34, 35].
4 数值模拟实例为了验证多震相走时联合(或同时)反演算法的有效性，我们选用了一个三维速度模型(图 3)，30km以上为起伏速度分布，而其下则为线性增加.图 3a为某一水平面(z= -24km)的速度分布，图 3b 则为某一垂直剖面(y=40km)的速度分布.为了进行两个界面的同时反演，模型中引入两个反射界面；上界面为深度z= -30km 的水平界面(模拟康氏面)，下界面则为倾斜台阶状界面(模拟莫氏面，深度位于z=-44~-48km 之间，图 4 所示).模型网格化时选用4.0km×4.0km×4.0km 正方体单元进行剖分，界面附近则采用尺度相同的不规则单元，次级节点的间距为0.5km.为了模拟天然地震，在模型外围周边z=-4.0km, -8.0km 和-12.0km处等间距放置了48个炮点(每层16 个)，其余6 个置于模型内部，共计54 个炮点.121 个检波器均匀地置于地表(间距为8.0km).这里主要讨论初至P波(射线条数为6534)、康氏面及莫氏界面的一次反射波(P1P及PmP，射线条数为13068)、以及地表与康氏面和莫氏界面的二次反射波(pP1P 及pPmP，射线条数同样为13068).在以下的数值模拟实验中，反演中各震相走时的权系数相同.但在实际问题中可选用不同的权系数.
|
图 3 三维速度模型 (a)z=-24km 水平面速度分布；(b)y=40km 垂直剖面速度分布. Fig. 3 3D velocity model (a) Horizontal velocity at z= -24 km； (b) Vertical velocity at v=40 km |
|
图 4 炮点(十字)、检波器(三角形)及两个反射界面(深灰色)分布图 Fig. 4 Source (cross) and receiver (triangle) layout and distribution of the two reflective interfaces (dark grey) |
为了分析讨论多震相走时联合反演成像的空间分辨率及其优越性，此模拟实验对上述五种波的不同组合进行了成像(反演过程中反射界面不更新).这样做的理由是检验反演算法的能力，为后面速度和反射界面的同时更新提供参考.选用一维层状线性速度模型作为初始模型，反演结果以与图 3 相同的形式给出(图 5:水平面结果；图 6:垂直剖面结果).其中，图 5b和6b是仅用直达P 波的结果；图 5c和6c是仅用一次反射波(P1P及PmP)的结果；图 5d和6d是仅用二次反射波(pP1P及pPmP)的结果；图 5e和6e则是结合直达P 波和一次反射波(P1P 及PmP)的结果；图 5f和6f则是结合直达P 波、一次反射波(P1P及PmP)和二次反射波(pP1P及pPmP)的结果.
|
图 5 不同震相走时联合反演成像(z=-24km)的结果 (a)真实模型；(b)直达P波；(c)反射P1P及PmP 波；(d)二次反射pP1P及pPmP波；(e)直达P波+反射P1P和PmP波；(f)直达P波+反射P1P及PmP波+二次反射pP1P及pPmP波. Fig. 5 Seismic tomography (horizontal slice at z=-24 km) combining different sets of traveltimes (a) Real model； (b) Direct P arrivals； (c) Reflected P1P and PmP arrivals； (d) Two-time reflected pP1P and pPmP arrivals； (e) Direct P arrivals + relectedP1P andPmP arrivals； (f) DirectP arrivals+relectedP1P andPmP arrivals + two-ime refected pP1P and pPmP arrivals. |
从图 5和图 6中可看出:仅用直达P 波走时资料进行成像，由于受炮检距的制约，其成像深度较浅，但最为主要的是成像的空间分辨率较低(图 5b和6b所示).仅用一次反射波或仅用二次反射波走时资料时，其成像结果大体相同，均好于直达P 波的结果.而采用多震相资料的联合反演结果则均好于单一震相的结果.成像分辨率的提高，不仅需要增加射线条数，更重要的是要增加不同射线的交叉概率，这样有利于速度异常体空间位置的锁定.另一方面，由于反射波资料的使用，成像的深度显然也就增加了，当然这取决于反射界面的深度.但从反演解与实际模型的差别来看，图 5f和6f的结果较之图 5e和6e的结果更接近真实速度模型(图略).从这个模拟实例可以得到这样的结论:采用不同震相的走时资料进行联合反演成像的优势是明显的，即(1)由于增加了射线密度，故而可显著提高成像的空间分辨率；(2)由于增加了交叉射线的概率，因而可减少重建速度场中的畸变.
4.2 层速度已知时反射界面的精确定位
地震勘探的主要目的就是确定不同波阻抗(反射)界面的位置及起伏细节，用地震射线进行偏移成像时，首先要确定(或估计)层间波速，然后叠加，进而偏移成像.这里假定用地震勘探的方法已初步得到了层间波速分布，则可用上述反演算法进行多个反射界面的同时精确定位.速度模型及炮检排列同上，利用一次反射P1P及PmP 波、或二次反射pP1P及pPmP波走时、或两者结合，采用上述反演算法(只更新界面)可对多个反射界面进行同时精确定位(图 7给出了双界面的定位结果).图 7b是仅用一次反射波的结果；图 7c是仅用二次反射波的结果；图 7d则是利用一次、二次反射波联合定位时的结果.从图 7中可看出，即使初始界面远离真实界面，一般经过几次迭代循环后均能够较好地确定反射界面的位置与起伏细节.采用一次反射波、二次反射波走时单独定位时，其平均误差分别为0.25km 和0.23km；而采用一次反射波和二次反射波联合定位时，其平均误差下降至0.12km.
|
图 7 层速度已知条件下对双反射界面的精确定位 (a)真实界面(黑色)和初始界面(灰色)；(b)一次反射波定位结果；(c)二次反射波定位结果；(d)一次和二次反射波联合定位结果. Fig. 7 The accurate location of reflective interface with known internal velocity (a) Real (dark) and initial (grey) interface； (b) Results for primary reflections； (c) Results for two-time reflections；(d)Results for combination of the primary and two-time reflections. |
在同时反演中，采用了两种参数归一化处理(见3.3节).初始模型的选用、炮检排列与4.1节相同，初始反射上界面水平置于真实界面下2km处，而初始反射下界面置于倾斜(参见图 10a所示).图 8b和9b为仅用反射P1P及PmP波走时进行同时反演的结果，图 8c和9c则为直达P波和反射P1P及PmP波走时联合同时反演的结果.从图 8 和图 9 中可看出，采用直达波和反射波走时联合的成像结果要好于单一的反射波走时的成像结果.与图 5、6相比较，即使在同时反演时，本文提出的反演算法也能得到满意的成像结果.与subspace多参数反演算法[13]相比，本文提出的反演算法要优于subspace多参数反演算法.有关细节将另文详述.
|
图 8 不同震相走时联合同时反演成像(水平面:z=-24km)结果 (a)真实模型；(b)反射P1P及PmP波；(c)直达P波+反射P1P及PmP波. Fig. 8 Simultaneous inversion (Horizontal slice at z=-24 km) with different sets of traveltimes (a) Real model； (b) Reflected P1P and PmP arrivals； (c) Direct P arrivals+reflected P1P and PmParrivals. |
图 10则为速度和界面同时反演时反射界面更新的结果(图 10b为仅用反射P1P及PmP 波走时的更新结果，图 10c则为用直达P 波和反射P1P 及 PmP波走时联合时的更新结果).从图中可看出，采用两种走时的联合要好于单一反射波走时的更新情况.采用单一反射波时，其界面定位的平均误差为0.92km, 而采用直达P波和反射波联合时，其定位平均误差下降至0.58km.从这个模拟实验中可初步得到这样的结论:(1)采用不同震相走时联合进行同时反演成像要好于采用单一反射波走时的成像结果；(2)当反射界面不太复杂时，同时反演成像的速度分布图像与联合成像的速度分布图像相当.值得注意的是:反演过程中初始反射界面的位置不能远离真实界面，否则，将造成反演解的不收敛或是重建的速度场及反射界面失真.当然，初始界面的选择可通过先验的信息进行粗略的估计，这一点在反演前是不难做到的.
|
图 10 利用不同震相走时组合同时反演成像中反射界面的更新 (a)真实界面(黑色)和初始界面(灰色)；(b)仅用反射P1P 及 PmP波走时；(c)直达P波+ 反射P1P及PmP波走时联合.图(b, c)中，黑色:真实界面；浅灰色:更新界面. Fig. 10 The subsurface interface updated in simultaneous inversion using different setsof traveltimes (a) Real (black) and initial (grey) interlace； (b) Relected P1P and PmP arrivals only； (c) Direct P arrivals ^reflected P1P and PmP arrivals.In diagrams (b, c)，dark: real subsurface and grey: updated subsurface rnterface. |
本文利用新近提出的三维复杂层状介质中多次波射线追踪正演技术，结合共轭梯度求解带约束阻尼最小二乘问题的反演算法，分析讨论了三维速度模型(含双反射界面)中多震相走时的联合(或同时)反演成像技术.为了兼顾各种震相走时的拾取误差不同，反演算法中引入了不同震相种类数据的权系数；为了平衡两种不同模型参数对走时的贡献不同，反演中引入了不同参数的归一化因子.数值模拟结果表明:多波走时联合成像是一种提高走时成像空间分辨率，进而降低重建模型失真度行之有效的方法技术.加之，宽频带数字地震仪的广泛使用，该算法具有广泛的实际应用价值.另一方面，多波走时联合同时成像技术也是深地震测深(或广角折反射地震剖面)技术首选的算法.它不仅能够重建速度场分布，又能够对反射(如:康氏和莫氏)界面进行精确定位.二维情形下，我们进行了反演算法对随机噪声的敏感性试验，结果表明:该反演算法在走时中随机误差强度不大(≤0.2s)，或模型中随机相对扰动强度不大(≤5.0%)时有效[39].目前工作的重点是如何在反演算法中引入地震射线密度及射线交叉概率，以满足利用天然地震资料进行细结构成像的实际需要.另外，对于天然地震资料，震源定位也存在不确定性，因此，如何在反演算法中同时更新三种模型参数，即:速度分布、震源位置及反射界面则是需要进一步深入研究的课题.
| [1] | Aki K, Lee W H K. Determination of three dimensional velocity anomalies under a seismic array using first P arrival times from local earthquakes. J Geophys Res , 1976, 81: 4381-4399. DOI:10.1029/JB081i023p04381 |
| [2] | Zhao D. Global tomographic images of mantle plumes and subducting slabs: insight into deep Earth dynamics. Phys Earth Planet Inter , 2004, 146: 3-34. DOI:10.1016/j.pepi.2003.07.032 |
| [3] | Lei J, Zhao D. A new insight into the Hawaiian plume. Earth Planet Sci Lett , 2006, 241: 438-453. DOI:10.1016/j.epsl.2005.11.038 |
| [4] | Lei J, Zhao D. Global P-wave tomography: on the effect of various mantle and core phases. Phys Earth Planet Inter , 2006, 154: 44-69. DOI:10.1016/j.pepi.2005.09.001 |
| [5] | Zhao D, Hasegawa A, Horiuchi S. Tomographic imaging of P and S wave velocity structure beneath northeastern Japan. J Geophys Res , 1992, 97: 19909-19928. DOI:10.1029/92JB00603 |
| [6] | Lei J, Xie F, Lan C, et al. Seismic images under the Beijing region inferred from P and PmP data. Phys Earth Planet Inter , 2008, 168: 134-146. DOI:10.1016/j.pepi.2008.06.005 |
| [7] | Zhao D, Todo S, Lei J. Local earthquake reflection tomography of the Landers aftershock area. Earth Planet Sci Lett , 2005, 235: 623-631. DOI:10.1016/j.epsl.2005.04.039 |
| [8] | 刘福田. 震源位置和速度结构的联合反演(1)----理论和方法. 地球物理学报 , 1984, 27(2): 167–175. Liu F T. Simultaneous inversion of earthquake hypocenters and velocity structure (I)—Theory and method. Chinese J Geophys (in Chinese) , 1984, 27(2): 167-175. |
| [9] | 华标龙, 刘福田. 界面和速度反射联合成像----理论与方法. 地球物理学报 , 1995, 38(6): 750–756. Hua B L, Liu F T. Joint reflection imaging of velocity and interface—Theory and method. Chinese J Geophys (in Chinese) , 1995, 38(6): 750-756. |
| [10] | 周龙泉, 刘福田, 陈晓非. 三维介质中速度结构和界面的联合成像. 地球物理学报 , 2006, 49(4): 1062–1067. Zhou L Q, Liu F T, Chen X F. Simultaneous tomography of 3-D velocity structure and interface. Chinese J Geophys (in Chinese) , 2006, 49(4): 1062-1067. |
| [11] | 李松林, 吴宁远, 宋占隆, 等. 速度分布和界面位置的联合反演. 地震学报 , 1997, 19(4): 383–392. Li S L, Wu N Y, Song Z L, et al. Simultaneous inversion of velocity distribution and interface positions. Acta Seismologica Sinica (in Chinese) , 1997, 19(4): 383-392. |
| [12] | Pavlis G L, Booker J R. The mixed discrete-continuous inversion problem: application to the simultaneous determination of earthquake hypocenter and velocity structure. J Geophys Res , 1980, 85: 4801-4810. DOI:10.1029/JB085iB09p04801 |
| [13] | Kennett B L N, Sambridge M S, Williamson P R. Subspace methods for large inverse problem with multiple parameter classes. Geophysical Journal , 1988, 94: 237-247. DOI:10.1111/j.1365-246X.1988.tb05898.x |
| [14] | Williamson P R. Tomographic inversion in reflection seismology. Geophys J Int , 1990, 100: 255-274. DOI:10.1111/j.1365-246X.1990.tb02484.x |
| [15] | Zelt C A, Hojka A M, Flueh E R, et al. 3D simultaneous seismic refraction and reflection tomography of wide-angle data from the central Chilean margin. Geophys Res Lett , 1999, 26: 2577-2580. DOI:10.1029/1999GL900545 |
| [16] | Zelt C A, Sain K, Naumenko J V, et al. Assessment of crustal velocity models using seismic refraction and reflection tomography. Geophys J Int , 2003, 153: 609-626. DOI:10.1046/j.1365-246X.2003.01919.x |
| [17] | Zelt C A, Ellis R M, Zelt B C. Three-dimensional structure across the Tintina strike-slip fault, northern Canadian Cordillera, from seismic refraction and reflection tomography. Geophys J Int , 2006, 167: 1292-1308. DOI:10.1111/gji.2006.167.issue-3 |
| [18] | McCaughey M, Singh S C. Simultaneous velocity and interface tomography of normal-incidence and wide-aperture seismic traveltime data. Geophys J Int , 1997, 131: 87-99. DOI:10.1111/gji.1997.131.issue-1 |
| [19] | James W, Hobro D, Satish C, et al. Three-dimensional tomographic inversion of combined reflection and refraction seismic traveltime data. Geophys J Int , 2002, 152: 79-93. |
| [20] | Zhang J, ten Brink U, Toksoz M. Nonlinear refraction and reflection travel time tomography. J Geophys Res , 1998, 103(B12): 29743-29757. DOI:10.1029/98JB01981 |
| [21] | Korenaga J, Holbrook W, Kent G, et al. Crustal structure of the southeast Greenland margin from joint refraction and reflection seismic tomography. J Geophys Res , 2000, 105(B9): 21591-21614. DOI:10.1029/2000JB900188 |
| [22] | 毛伟健, StuaG W. 多震相的层析成像与在介质中的多波型射线跟踪. CT理论与应用研究 , 1996, 5(4): 42–48. Mao W J, Stua G W. On multi-phase tomography and wavetype ray tracing in earth-media. CT Theory and Applicaiion (in Chinese) , 1996, 5(4): 42-48. |
| [23] | Rawlinson N, Sambridge M. Multiple reflection and transmission phases in complex layered media using a multistage fast marching method. Geophysics , 2004, 69: 1338-1350. DOI:10.1190/1.1801950 |
| [24] | De Kool M, Rawlinson N, Sambridge M. A practical grid- based method for tracking multiple refraction and reflection phases in three-dimensional heterogeneous media. Geophys J Int , 2006, 167: 253-270. DOI:10.1111/gji.2006.167.issue-1 |
| [25] | Bai C Y, Tang X P, Zhao R. 2D/3D multiply transmitted, converted and reflected arrivals in complex layered media with the modified shortest path method. Geophys J Int , 2009, 179: 201-214. DOI:10.1111/gji.2009.179.issue-1 |
| [26] | 唐小平, 白超英. 最短路径算法下三维层状介质中多次波追踪. 地球物理学报 , 2009, 52(10): 2635–2643. Tang X P, Bai C Y. Multiple ray tracing within 3-D layered media with the shorted path method. Chinese J Geophys (in Chinese) , 2009, 52(10): 2635-2643. |
| [27] | 唐小平, 白超英. 利用改进后的最短路径法追踪二维复杂层状介质中的多次波. 地球物理学进展 , 2009, 24(6): 2087–2096. Tang X P, Bai C Y. Multiple ray tracing within 2D layered media with the shorted path method. Progress in Geophys (in Chinese) , 2009, 24(6): 2087-2096. |
| [28] | 赵瑞, 白超英. 复杂层状模型中多次波快速追踪----一种基于非规则网格的最短路径算法. 地震学报 , 2010, 32(4): 433–444. ZhaoR, Bai C Y. Fast multiple ray tracing within complex layered media： the shortestpath methodbasedon irregular grid cells. Acta Seismologica Sinica (in Chinese) , 2010, 32(4): 433-444. |
| [29] | Bai C Y, Huang G J, Zhao R. 2D/3D irregular shortest-path raytracing for multiple arrivals and its applications. Geophys J Int , 2010, 183: 1596-1612. DOI:10.1111/j.1365-246X.2010.04817.x |
| [30] | Menke W. Geophysical Data Analysis: Discrete Inverse Theory. San Diego. California: Academic Press Inc., 1984 : 40 -126. |
| [31] | Meyerholtz K A, Szpakowski S A. Convolutional quelling in seismic tomography. Geophysics , 1989, 54: 570-580. DOI:10.1190/1.1442684 |
| [32] | Shaw P R, Orcutt J A. Waveform inversion of seismic refraction data and applications to young Pacific crust. Geophys J R Astr Soc , 1985, 82: 375-414. DOI:10.1111/j.1365-246X.1985.tb05143.x |
| [33] | Tarantola A, Valette B. Generalized non-linear inverse problems solved using the least squares criterion. Rev Geophys Space Phys , 1982, 20: 219-232. DOI:10.1029/RG020i002p00219 |
| [34] | Zhou B, Sinadinovski C, Greenhalgh S A. Iterative algorithm for the damped minimum norm, least-squares and constrained problem in seismic tomography. Explor Geophys , 1992, 23: 497-505. DOI:10.1071/EG992497 |
| [35] | Bai C Y, Greenhalgh S A. 3-D non-linear travel time tomography: imaging high contrast velocity anomalies. Pure Appl Geophys , 2005, 162: 2029-2049. DOI:10.1007/s00024-005-2703-x |
| [36] | Bishop T N, Bube K P, Cutler R T, et al. Tomographic determination of velocity and depth in laterally varying media. Geophysics , 1985, 50: 903-923. DOI:10.1190/1.1441970 |
| [37] | Bai C Y, Greenhalgh S. 3-D local earthquake hypocenter determination with an 'irregular' shortest-path method. Bull Seism Soc Am , 2006, 96: 2257-2268. DOI:10.1785/0120040178 |
| [38] | Scales J. Tomographic inversion via the conjugate gradient method. Geophysics , 1987, 52: 179-185. DOI:10.1190/1.1442293 |
| [39] | 黄国娇, 白超英. 二维复杂层状介质中地震多波走时联合反演成像. 地球物理学报 , 2010, 53(12): 2972–2981. Huang G J, Bti C Y. Simultaneous inversion with multiple traveltimes within 2-D complex layered media. Chinese J Gephys (in Chinese) , 2010, 53(12): 2972-2981. |