2013, Vol. 56
2. 长安大学计算地球物理研究所, 西安 710054
2. Institute of Computing Geophysics, Chang'an University, Xi'an 710054, China
地震层析成像经历了30多年的发展, 形成了众多的方法技术, 已成为人们研究地球内部结构及各向异性最为行之有效的方法技术之一.其成像尺度含盖范围很广, 小到岩石组分, 大至全球构造.地震走时成像技术的分辨率除了受制于高频假设外, 主要是射线覆盖率较低.因此, 在现有地震和地震台网分布情况下, 要想提高分辨率, 采用多震相走时资料是唯一的有效途径.多震相走时成像较为典型的一个例子是赵大鹏等利用美国Landers地震余震进行的多震相(S、SmS、sSmS)走时同时对S波速度和Moho面成像[1], 其中仅用了相距200km的两个观测台站和其间发生的180个余震, 结果表明使用多震相走时成像其空间分辩率可提高至台站间距的1/8~1/6, 本例中为25~35km.
传统的走时同时反演成像方法大多基于双参数的同时更新, 即速度模型和震源位置的同时反演[2-5]或是速度模型和反射界面的同时更新[6-10].然而, 双参数反演是一项比较困难的工作, 因为界面(或震源位置)的微小变化会引起射线(初至波或反射、折射波)路径的较大变化, 加上速度和界面(或震源位置)的强耦合关系, 可能会导致反演失败.即在一个方程组中同时反演两种(或多种)模型参数, 解会向系数大的方向偏离.解决这一问题, 目前主要有4种不同的方法:一是参数分离的方法, 例如:刘福田提出的正交投影算子分解法, 对速度和震源位置(或界面)两种不同的模型参数进行解耦[11-14]; 二是Kennett等提出的子空间(Subspace)反演算法[15], 能同时有效地反演多种模型参数; 三是对速度和界面(或震源位置)两种不同的模型参数进行加权[10]; 四是反演前对两种不同的模型参数进行归一化处理[8, 16-19].
目前为止走时三参数(速度、反射点深度和震源位置)同时反演研究并不多见, Sambridge (1990)或许是首位讨论走时三参数同时反演的研究学者, 但其侧重点主要是在反演中如何引入正则化约束来寻求三维光滑模型, 而不是三参数同时反演算法本身[20].其中正演采用边界值两点射线追踪算法[21], 反演则采用子空间反演方法.随后, Preston (2003)在他的博士论文中也进行了走时三参数同时反演的研究工作[22], 其中正演分为两部分:采用有限差分解程函方程的算法进行初至波射线追踪[23], 采用两步法进行反射波的追踪计算[24], 反演算法则采用共轭梯度求解最小二乘算法(CGLS)[25], 同时反演中引入了正则化并对三种不同的模型参数进行了加权处理.其结果发表在2003年的Science期刊上[26], 这也从侧面说明三参数同时反演是可行的.
然而, 上述走时三参数同时反演算法, 要么是正演采用传统的射线追踪算法[21], 要么是反演中采用受人为因素影响的加权处理[22], 都有可能影响最终的成像结果.另外, 地震资料主要是利用初至波和反射波走时数据, 并没有考虑多震相走时资料.我们知道, 走时成像的分辨率除了受高频假设制约外, 主要受制于射线覆盖率, 确切地说, 受制于不同震相射线的交错概率.因此, 走时三参数同时反演的成败与否, 其一是尽可能采用多震相走时资料; 其二是采用合适的反演算法.鉴于此, 本文将目前流行的基于网格/单元波前扩展的多震相射线追踪算法之一(分区多步不规则最短路径算法)和地震多参数反演算法(子空间反演法)相结合, 提出了一种可利用地震多震相走时资料进行三参数同时反演的方法技术.与走时双参数同时反演成像结果对比分析, 所提出的三参数同时反演算法是一种行之有效的地震走时三参数同时反演成像方法技术, 同时该算法对走时数据中的随机噪声并不敏感, 说明具有较好的实用价值.
2 地震多震相射线追踪算法简介目前基于网格/单元波前扩展的多震相射线追踪算法主要有两种:一是基于有限差分解程函方程的分区多步快速行进法(Multistage FMM)[27-28]; 其二是基于分区多步不规则最短路径算法(Multistage ISPM)[29-30].我们曾就上述两种算法进行了对比研究, 结果表明, 无论是计算精度还是CPU时间, 分区多步不规则最短路径算法均优于分区多步快速行进法[29].
这里所谓分区多步计算中的分区是指在模型参数化时, 可根据速度界面的具体分布情况将模型分成若干个独立的计算区域(相邻的区域由界面连接); 而所谓多步计算则是根据所要计算地震震相的种类从炮点所在的区域逐区(步)进行计算.具体来说, 自炮点所在的单元开始进行扫描, 按照一定的步长进行波前扩展.等当前区域所有网格单元内节点都计算完毕后, 波前停留在该区的速度离散界面上.若追踪下行透射(或透射转换)波, 则从速度界面中走时最小的点开始(惠更斯原理, 视为次级震源), 继续新区内的波前扩展和射线追踪; 若追踪上行反射(或反射转换)波, 则从速度界面中走时最小的点开始, 继续原区内的波前扩展和射线追踪.若存在转换波, 则调用相应的速度模型.按照上述步骤重复则可实现多震相地震射线的追踪计算[29-30].作为一个例子, 图 1给出了多震相地震射线路径(所用速度模型和界面位置分布参见随后反演实例), 其中图 1a显示了来自两个不同界面的反射波的射线路径; 图 1b则为自炮点激发的P波上行至地表, 经反射后分别下行至第一、二个反射界面, 然后反射且上行至检波器的多次反射波的射线路径.
|
图 1 (a)透射和一次反射波射线路径.图中黑线:反射波震相P1P1; 灰线:反射波震相P1P2P2P1; (b)多次透射、反射波的射线路径.图中黑线:反射震相P1P1P1; 灰线:多次透射、反射波震相P1P1P2P2P1 Fig. 1 (a) Raypaths of the transmitted and primary reflections. In the figure the black lines are the raypaths of phase P1P1 and grey lines are the raypaths of P1P2P2P1 phase; (b) Raypaths of the multiple transmissions, conversions and reflections.In the figure the black lines are the raypaths of P1P1P1 phase and the grey lines are the raypaths of P1P1P2P2P1 phase |
Skilling和Bryan[31]在最大熵图像重建中提出子空间(subspace)反演算法, 随后Kennett等[15]将其概念推广, 并讨论了地震多参数的反演问题.其应用研究主要集中在反射或折射波成像, 即速度和反射界面的同时反演[6, 27, 32-35].子空间反演算法的优点就在于其计算效率和实用功能, 这些都是通过正确的选择跨越模型子空间的基向量达到的.同样, 子空间反演算法也可进行正则化和平滑.
3.1 Subspace反演算法简介广义子空间反演算法的核心是将多种参数的模型扰动量(δm=m-mk∈RmN)用给定的基向量{bi}(i=1, 2, …q, 生成的q维子空间Rq⊂RNm)来表示, 即[36]
|
(1) |
其中
|
(2) |
对(2)式求有关αk的导数, 并取极小值(∂Φ/∂αk=0), 则有
|
(3) |
将二阶导数(3)式带入(2)式, 并引入正则化, 则有
|
(4) |
其中
|
(5) |
|
(6) |
|
(7) |
|
(8) |
Wd和Wm分别为数据加权矩阵和模型加权矩阵, d0和d(m)分别是观测数据和理论计算数据, λ为阻尼因子.由(1)和(4)式, 可得子空间反演算法的迭代形式解为
|
(9) |
显然, 当B=I(单位矩阵)时, 子空间近似就褪变成全空间的近似迭代算法.
从方程(5)-(7)式可以看出, 当模型未知参数(Nm)很大时, 若基向量{bi}(i=1, 2, …, q)中q较小, 且选择合适, 则(5)式反演迭代解中广义逆矩阵的维数(q×q)要比原来的维数(Nm×Nm)小得多.
但是, 若基向量是线性相关的, 则矩阵可能是奇异的.因此, 首先需要将这些向量正交化.选择好基向量后, 迭代近似解(5)式中的转换矩阵AkT WdAk, Q'k和W'm矩阵的维数就小的多, 当给定阻尼因子λ后, 可用奇异值分解法求逆矩阵的方法求解, 或者用共轭梯度法求解也是很有效的.具体计算步骤如下:
第一步:选择基向量B=(b1, b2, …, bq), 开始迭代ml*=m0*, l=0, 1, 2, …, lmax;
第二步:搜索最佳阻尼因子{λi}, λi∈[λmin,
λmax], 并求子空间近似值: 
第三步:获得当前解:
第四步:检验
从上述实现步骤来看, 基向量的选择是关键. Oldenburg [37]给出了一种选择基向量的方法, 将基向量同数据梯度、模型残差函数联系起来:
|
(10) |
|
(11) |
三维情形下根据射线理论, 由第i个震源经反射点到达第j个检波器的走时可表述为下列积分:
|
(12) |
其中ds为沿路径Rij上的线积分元, 而Vc(x, y, z)则为沿射线路径上的速度分布函数.当模型网格单元化后, 沿该射线路径Rij上的走时可写成求和形式:
|
(13) |
其中k为射线穿过的单元总数, Rijk、Vck(x, y, z)分别为射线穿过第k个单元的长度和速度分布.
在速度、界面和震源位置同时反演时, Jacob偏导矩阵包含了三部分:一是走时关于速度变化的导数, 二是走时关于反射点深度变化的导数, 三是走时关于震源参数的导数.
|
(14) |
其中vk和zk分别是第k个单元内的速度分布和第k个反射点的深度, Δvk和Δzk则分别是相应速度和反射点深度的扰动量.Δxk(k=1, 2, 3)是第k个震源分别在X, Y和Z方向的扰动量.
(14)式中第一项, 即有关走时关于速度的一阶偏导数计算如下:为简便起见, 这里略去震源编号i, 则tj是第j条射线Rj的走时; vk是长度为N的模型空间中的第k个单元的速度分布; zk是离散界面点中第k个点的深度值.
当射线穿过某个单元时, 走时对速度的一阶导数可用射线两端点导数的平均值代替[38]:
|
(15) |
其中
|
式中ki, kj是穿过第k单元射线段的两个端点; BV(k)则是该单元上主节点的速度值. (14)式中的第二项, 即走时关于反射点深度变化的偏导数则可由以下步骤计算.若考虑一般情况(反射、折射、以及反射转换), 则走时关于反射点深度变化的导数由下式确定[20]:
|
(16) |
其中zn和zint分别是界面反射点和初始反射点的深度(有关wj, wj+1, wn和wz, 见图 2所示), vj和vj+1则分别是界面两侧的速度值.当仅考虑反射时, wj+1· wj=-wj·wn和vj=vj+1, 则(16)式简化为下式:
|
(17) |
|
图 2 走时关于反射点深度偏导数计算示意图 Fig. 2 Diagrammatic showing the calculation of traveltime derivative with respec tto reflector depth change |
(14)式的第三项, 即到时关于震源参数(位置和发震时刻)变化的偏导数则可由以下步骤计算, 首先到时(t=ta-to)关于发震时刻的导数为
|
(18) |
而到时关于震源位置变化的偏导数的计算可参见图 3所示.若震源的深度扰动值为Δzh, 则到时关于震源位置三方向(x, y, z)变化的偏导数可由下式计算[11, 20]:
|
(19) |
其中θ如图 3所示, φ和ψ则分别是射线水平面投影后与X和Y轴的夹角.
|
图 3 走时关于震源位置变化偏导数计算示意图 Fig. 3 Diagrammatic showing the calculation of traveltime derivative with respect to source location change |
为了验证三参数同时反演的有效性, 我们选取了一个较为复杂的三维速度模型(见图 4所示, 模型尺度:200km×200km×50km).在Z=20km和48km深度存在速度界面, 前者为起伏波浪式界面(幅度:±5.0km), 后者为水平界面(见图 8a所示).模型参数化采用4km×4km×2.5km单元(起伏界面附近则采用尺度相当的不规则单元), 70个震源随机分布在5~12km深度范围内, 121个检波器(20km×20km间距)均匀地分布于地面(见图 5所示), 其中最小震中距为2.8km, 最大震中距有262.6km.为了接近实际情况, 初始震源进行了随机扰动(见图 11a所示), 而初始反射界面均置于水平(见图 9b所示, 其中上界面置于20km深度, 下界面置于45km深度, 初始下界面距真实界面垂直距离为3km).模拟观测走时中引入了随机误差(强度:初至P波为±0.2s, 反射波为±0.3s).三种不同反演算法(见下节说明)采用了三种不同的震相走时资料, 即初至P波和来自上、下界面的反射PP波走时.但是考虑实际资料处理中, 并不是所有的台站都有可识别的记录, 所以本文中选用震中距大于30 km的反射波走时数据, 则数值模拟中采用初至P波走时数据8470个, 上界面和下界面的反射PP波走时数据均为6250个.图 6给出了Z=25.0km水平切片(图 6a)和Y=104.0km的垂直切片(图 6b)上各个单元内穿过本文所采用的所有三种地震射线条数的统计分布图.从图中可以看出, 在上界面的峰值处和下界面中间段射线密度较大(图 6b), 参照真实模型图 4a, 我们发现射线在高速区域密集而在低速区域稀疏(图 6a).
|
图 4 三维速度模型 (a)Z=25km水平速度场分布; (b)Y=104km垂直剖面速度场分布; (c)相同水平速度场和一维层状背景速度场的百分比差值, 或称之为速度异常场; (d)相同垂直剖面速度场和一维层状背景速度场的百分比差值, 或称之为速度异常场. Fig. 4 3D velocity model (a) Horizontal slice at depth of 25km; (b) Cross section view at Y=104km; (c) Percentage perturbation between the real and 1-D layered velocity field at the depth of 25km; (d) Percentage perturbation between the real and 1-D layered velocity field at one cross section view of Y=104km. |
|
图 5 震源(十字)和检波器(三角形)分布 Fig. 5 Source location (cross) and receiver (triangle) distribution |
|
图 6 水平面(a)(Z=25.0km)和垂直剖面(b)(Y=104.0km)切片上各单元内穿过地震射线条数的统计分布 Fig. 6 The numbers of account rays through each cell crossing horizontal slice (a)(Z=25.0km) and vertical slice (b)(Y=104.0km) |
|
图 8 三种不同反演算法在Y=104km垂直剖面的成像结果(a, b, c)和相同剖面三种反演解与真实速度的百分比误差(d, e, f)(a, d) V+R反演; (b, e) V+S反演; (c, f) V+R+S反演. Fig. 8 The recovered velocity fields at Y=104km cross section for three different inversion methods (a, b, c) and the percentage differences between the real and the inverted velocity fields at the same cross section (d, e, f) (a, d) V+R result; (b, e) V+S result; (c, f) V+R+S result. |
|
图 9 两种不同反演算法更新的反射界面分布 (a)真实反射界面; (b)初始反射界面; (c) V+R反演结果; (d) V+R+S反演结果. Fig. 9 The updated reflected interfaces by the two different inversion algorithms (a) Thereal subsurface interfaces; (b) The initial subsurface interfaces; (c) Results for the V+R inversion; (d) Result for the V+R+S inversion. |
|
图 11 两种不同反演算法得到的震源位置相对于真实位置的收敛情况 (a)初始扰动的震源位置分布; (b) V+S反演后相对于真实震源位置的收敛情况; (c) V+R+S反演后相对于真实震源位置的收敛情况. Fig. 11 The convergence to the real source locations after inversion by two different algorithms (a) The distance distribution between the real and the initial source locations; (b) The distance convergence to the real source locationsa fter V+S inversion; (c) The distance convergence to the real source locations after V+R+S inversion. |
为了进行系统的对比研究, 模拟实验中进行了三种不同的反演, 其一是速度和界面的同时反演(简记为:V+R反演), 其二是速度和震源位置的同时反演(简记为:V+S反演), 其三是速度、界面和震源位置的同时反演(简记为:V+R+S反演).不失一般性, 三种不同的反演迭代均为25次, 相应的走时残差分别由初始时的0.467s, 0.465s, 0.496s下降至0.057s (V+S反演), 0.055s (V+R反演), 0.061s (V+R+S反演).
图 7给出了相应水平面(Z=25km)三种不同反演算法的成像结果, 从中可以看出, 三种不同的反演算法重建的速度场与真实速度场有相似的异常结构(参见图 4a和图 7(a-c)), 但只是部分恢复(参见图 4c和图 7(d-f)).这一点从垂直剖面的速度场重建中也可进一步证实, 图 8给出了相应垂直剖面(Y=104km)三种不同反演算法得到的速度场分布.无论是重建的水平还是垂直剖面速度场, 双参数(V+R或V+S)反演略好于三参数(V+R+S)反演结果, 但三种不同的反演算法所得到的图像大体类似, 这也从侧面说明三参数同时反演成像是可行的.
|
图 7 三种不同反演算法在Z=25km深度的成像结果(a, b, c)和相同深度三种反演解与真实速度的百分比误差(d, e, f)(a, d) V+R反演; (b, e) V+S反演; (c, f) V+R+S反演. Fig. 7 The recovered velocity fields at the depth of 25km for three different inversion methods (a, b, c) and the percentage differences between the real and the inverted velocity fields at the same depth (d, e, f) (a, d) V+R result; (b, e) V+S result; (c, f) V+R+S result.` |
图 9中给出了两种反演算法更新后的反射界面分布(图 9c:V+R反演结果; 图 9d:V+R+S反演结果).从更新的界面来看, 下界面(水平界面)的反演结果两种算法基本相同, 而上界面(起伏波浪式界面), 则双参数(V+R算法)反演结果略好于三参数(V+R+S算法)的反演结果.这一点也是在预料之中, 毕竟三参数同时反演的难度要远大于双参数反演的难度.然而, 无论怎样, 波浪式起伏界面能够更新到这种程度(图 9d所示)也再次证明三参数同时反演的可行性和有效性.为更进一步评价界面更新的结果, 我们给出了两种反演更新的界面与真实界面间的差异图(图 10).从图 10中可以看出两种反演对界面的更新程度还是比较理想的, 只是下界面(图 10b和图 10d)在界面的四个角上更新相对较差, 这主要是因为穿过这里的射线较少.
|
图 10 两种不同算法更新界面与真实界面间残差 (a) V+R反演上界面; (b) V+R反演下界面; (c) V+R+S反演上界面; (d) V+R+S反演下界面. Fig. 10 The differences between the real and the updated reflected interfaces by the two different inversion algorithms (a) The upper interface by V+R method; (b) The lower interface by V+R method; (c) The upper interface by V+R+S method; (d) The lower interface by V+R+S method. |
三参数反演中不仅是速度模型和反射界面要同时更新, 同时震源位置也在更新.图 11a中给出了70个初始震源距离真实震源空间距离的频度分布情况(最大距离3.5km), 图 11(b-c)分别是双参数(V+S算法)、三参数(V+R+S算法)同时反演后震源位置收敛于真实震源位置的分布情况, 从中可以看出双参数的收敛程度要略好于三参数的收敛情况.然而, 无论怎样, 三参数同时反演中震源位置的收敛还是较为理想的, 最大空间距离由原来的3.5km收敛至1.25km, 绝大多数震源位置与真实震源位置的距离在0.5km内.
就一般的双参数(V+S算法)同时反演算法而言, 震源位置的重新定位只有相对意义, 而无绝对意义, 加之多采用初至P或P-S到时, 当台网分布不理想时, 水平向的定位精度要远好于垂直向的定位精度.然而, 若采用反射波或深震资料, 则震源深度的定位精度将会有较大的改善.为了说明这一点, 图 12中给出了两种不同反演算法反演后震源深度收敛于真实震源深度的分布情况, 图 12b为V+S反演后的收敛情况, 而图 12c则是V+R+S反演后的收敛情况.从中可以看出, 尽管双参数(V+S算法)反演后震源深度的收敛情况要略好于三参数(V+R+S算法)反演后震源深度的收敛情况.但是, 三参数同时反演过程中使得起初最大震源深度误差2.0km降至1.0km, 且绝大多数震源深度收敛于0.5km内, 这说明无论是双参数还是三参数同时反演算法中震源深度的定位误差与三维空间震源定位误差属于同一精度范围, 即不存在震源深度定位较差的问题.
|
图 12 两种不同反演算法得到的震源深度相对于真实震源深度的收敛情况 (a)初始扰动的震源深度分布; (b) V+S反演后相对于真实震源深度的收敛情况; (c) V+R+S反演后相对于真实震源深度的收敛情况. Fig. 12 The convergence to the real focus depths after inversion by two differental gorithms. (a) The distance distribution between the real and the initial focus depths; (b) The distance convergences to the real focus depths after V+S inversion; (c) The distance convergences to the real focus depths after V+R+S inversion. |
现有的地震台网和地震分布情况下, 为了克服走时成像分辨率较低的问题, 则必须使用多震相走时资料以此来增加地震射线的覆盖率, 或者确切地说是不同震相射线的交错概率.然而, 实际走时成像中不但是速度和反射界面具体细节未知, 就是震源位置也存在不确定性(区域成像中更是如此).但由于多参数反演问题的不确定性及难度, 人们大多将其简化成双参数反演问题, 即速度和震源位置或者是速度和反射界面的同时更新问题.除了上述所说的复杂性外, 还在于缺少实用高效的正、反演算法.随着计算技术的日益更新发展, 已出现了成熟的多震相射线追踪算法和多参数反演算法.鉴于此, 本文将流行的多震相射线追踪正演算法和多参数反演算法相结合, 提出了一种三参数(速度、反射界面和震源位置)同时反演成像的方法技术, 数值模拟实例及与双参数(V+R或V+S)反演结果对比分析表明, 本文提出的三参数同时反演方法技术不失为一种高效、实用的反演算法.另外, 无论是双参数或是三参数反演算法均对到时数据中的随机误差不太敏感, 这点在实际应用中也是至关重要的.
| [1] | Zhao D P, Todo S, Lei J S. Local earthquake reflection tomography of the Landers aftershock area. Earth Planet. Sci. Lett. , 2005, 235(3-4): 623-631. DOI:10.1016/j.epsl.2005.04.039 |
| [2] | Pavlis G L, Booker J R. The mixed discrete-continuous inverse problem: application to the simultaneous determination of earthquake hypocenters and velocity structure. J. Geophys. Res. , 1980, 85(B9): 4801-4810. DOI:10.1029/JB085iB09p04801 |
| [3] | Thurber C H. Earthquake locations and three-dimensional crustal structure in the Coyote Lake area, central California. J. Geophys. Res. , 1983, 88(B10): 8226-8236. DOI:10.1029/JB088iB10p08226 |
| [4] | Protti M, Schwartz S Y, Zandt G. Simultaneous inversion for earthquake location and velocity structure beneath central Costa Rica. Bull. Seism. Soc. Am. , 1996, 86. |
| [5] | Asad A M, Pullammanappallil S K, Anooshehpoor A, et al. Inversion of travel time data for earthquake locations and three-dimensional velocity structure in the Eureka Valley area, eastern California. Bull. Seism. Soc. Am. , 1999, 89: 796-810. |
| [6] | Williamson P R. Tomographic inversion in reflection seismology. Geophys. J. Int. , 1990, 100(2): 255-274. DOI:10.1111/j.1365-246X.1990.tb02484.x |
| [7] | 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(16): 2577-2580. DOI:10.1029/1999GL900545 |
| [8] | McCaughey M, Singh S C. Simultaneous velocity and interface tomography of normal-incidence and wide-aperture seismic traveltime data. Geophys. J. Int. , 1997, 131(1): 87-99. DOI:10.1111/gji.1997.131.issue-1 |
| [9] | Rawlinson N, Houseman G A, Collins C D N. Inversion of seismic reflaction and wide-angle reflection traveltimes for three-dimensional layered crustal structure. Geophys. J. Int. , 2001, 145(2): 381-401. DOI:10.1046/j.1365-246x.2001.01383.x |
| [10] | 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(3): 609-626. DOI:10.1046/j.1365-246X.2003.01919.x |
| [11] | Thurber C H. Hypocenter-velocity structure coupling in local earthquake tomography. Phys. Earth Plan. Int. , 1992, 75(1-3): 55-62. DOI:10.1016/0031-9201(92)90117-E |
| [12] | 刘福田. 震源位置和速度结构的联合反演(I)-理论和方法. 地球物理学报 , 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. |
| [13] | 华标龙, 刘福田. 界面和速度反射联合成像-理论与方法. 地球物理学报 , 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. |
| [14] | 周龙泉, 刘福田, 陈晓非. 三维介质中速度结构和界面的联合成像. 地球物理学报 , 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. |
| [15] | Kennett B L N, Sambridge M S, Williamson P R. Subspace methods for large inverse problems with multiple parameter classes. Geophys. J. , 1988, 94(2): 237-247. DOI:10.1111/j.1365-246X.1988.tb05898.x |
| [16] | Hobro J W D, Singh S C, Minshull T A. Three-dimensional tomographic inversion of combined reflection and refraction seismic traveltime data. Geophys. J. Int. , 2003, 152(1): 79-93. DOI:10.1046/j.1365-246X.2003.01822.x |
| [17] | Huang G J, Bai C Y, Zhu D L, et al. 2D/3D seismic simultaneous inversion for the velocity and interface geometry using multiple classes of arrivals. Bull. Seism. Soc. Am. , 2012, 102(2): 790-801. DOI:10.1785/0120110155 |
| [18] | 黄国娇, 白超英. 二维复杂层状介质中地震多波走时联合反演成像. 地球物理学报 , 2010, 53(12): 2972–2981. Huang G J, Bai C Y. Simultaneous inversion with multiple traveltimes within 2-D complex layered media. Chinese J. Geophys. (in Chinese) , 2010, 53(12): 2972-2981. |
| [19] | 白超英, 黄国娇, 李忠生. 三维复杂层状介质中多震相走时联合反演成像. 地球物理学报 , 2011, 54(1): 182–192. Bai C Y, Huang G J, Li Z S. Simultaneous inversion combining multiple-phase travel times within 3D complex layered media. Chinese J. Geophys. (in Chinese) , 2011, 54(1): 182-192. |
| [20] | Sambridge M S. Non-linear arrival time inversion: Constraining velocity anomalies by seeking smooth models in 3-D. Geophys. J. Int. , 1990, 102(3): 653-677. DOI:10.1111/gji.1990.102.issue-3 |
| [21] | Sambridge M S, Kennett B L N. Boundary value ray tracing in a heterogeneous medium: A simple and versatile algorithm. Geophys. J. Int. , 1989, 101(1): 157-168. |
| [22] | Preston L A. Simultaneous inversion of 3D velocity structure, hypocenter locations, and reflector geometry in Cascadia. University of Washington , 2003. |
| [23] | Hole J. Nonlinear high-resolution three-dimensional seismic travel time tomography. Geophys. J. Int. , 1992, 97(B5): 6553-6562. DOI:10.1029/92JB00235 |
| [24] | Hole J, Zelt B. 3-D finite-difference reflection travel times. Geophys. J. Int. , 1995, 121(2): 427-434. DOI:10.1111/gji.1995.121.issue-2 |
| [25] | Paige C, Saunders M. LSQR: An algorithm for sparse linear equations and sparse least squares. ACM Trans. on Math. Soft. , 1982, 8(1): 43-71. DOI:10.1145/355984.355989 |
| [26] | Preston L A, Creager K C, Crosson R S, et al. Intraslab earthquakes: Dehydration of the Cascadia slab. Science , 2003, 302(5648): 1197-1200. DOI:10.1126/science.1090751 |
| [27] | Rawlinson N, Sambridge M. Multiple reflection and transmission phases in complex layered media using a multistage fast marching method. Geophysics , 2004, 69(5): 1338-1350. DOI:10.1190/1.1801950 |
| [28] | 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(1): 253-270. DOI:10.1111/gji.2006.167.issue-1 |
| [29] | Bai C Y, Huang G J, Zhao R. 2-D/3-D irregular shortest-path ray tracing for multiple arrivals and its applications. Geophys. J. Int. , 2010, 183(3): 1596-1612. DOI:10.1111/j.1365-246X.2010.04817.x |
| [30] | Bai C Y, Li X L, Tang X P. Seismic wavefront evolution of multiply reflected, transmitted, and converted phases in 2D/3D triangular cell model. Journal of Seismology , 2011, 15(4): 637-652. DOI:10.1007/s10950-011-9242-y |
| [31] | Skilling J, Bryan R K. Maximum entropy image reconstruction, general algorithm Mon. Not. R. Astron. Soc. , 1984, 211: 111-24. DOI:10.1093/mnras/211.1.111 |
| [32] | Blundell C A. Resolution analysis of seismic P-wave velocity estimates using reflection tomographic inversion. Monash University , 1993. |
| [33] | Wang Y H, Houseman G A. Inversion of reflection seismic amplitude data for interface geometry. Geophys. J. Int. , 1994, 117(1): 92-110. DOI:10.1111/gji.1994.117.issue-1 |
| [34] | Wang Y, Houseman G A. Tomographic inversion of reflection seismic amplitude data for velocity variation. Geophys. J. Int. , 1995, 123(2): 355-372. DOI:10.1111/gji.1995.123.issue-2 |
| [35] | Wang B, Braile L W. Simultaneous inversion of reflection and refraction seismic data and application to field data from the northern Rio Grande rift. Geophys. J. Int. , 1996, 125(2): 443-458. DOI:10.1111/gji.1996.125.issue-2 |
| [36] | Greenhalgh S, Zhou B. Solutions, algorithms and inter-relations for local minimization search geophysical inversion. J. Geophys. Eng. , 2006, 3(2): 101-113. DOI:10.1088/1742-2132/3/2/001 |
| [37] | Oldenburg D W, McGillicray P R, Ellis R G. Generalised subspace methods for large-scale inverse problems. Geophys. J. Int. , 1993, 114: 12-20. DOI:10.1111/gji.1993.114.issue-1 |
| [38] | Bai C Y, Greenhalgh S. 3-D non-linear travel-time tomography: Imaging high contrast velocity anomalies. Pure Appl. Geophys. , 2005, 162(11): 2029-2049. DOI:10.1007/s00024-005-2703-x |