地球物理学进展  2014, Vol. 29 Issue (5): 2019-2028   PDF    
引用本文
李兴旺, 白超英, 于子超. 菲涅耳体有限频射线多震相走时同时反演成像[J]. 地球物理学进展, 2014, 29(5): 2019-2028, doi: 10.6038/pg20140505  
LI Xing-wang, BAI Chao-ying, YU Zi-chao. Simultaneous inversion for velocity and reflector geometry based on the multi-phase fresnel volume rays. Progress in Geophysics, 2014, 29(5): 2019-2028, doi: 10.6038/pg20140505   
菲涅耳体有限频射线多震相走时同时反演成像
李兴旺1, 白超英1,2 , 于子超1    
1. 中国西安长安大学地质工程与测绘学院地球物理系, 710054;
2. 中国西安长安大学计算地球物理研究所, 710054
收稿日期: 2013-11-28 修回日期: 2014-6-16 .
基金项目：教育部高校博士学科点基金项目(20110205110010)资助.
作者简介：李兴旺, 男, 1989年9月生于四川安岳, 长安大学地球探测与信息技术专业博士研究生, 主要从事地震正、反演方法及应用研究. (E-mail: lixingwang19890913@163.com)
通讯作者：白超英, 男, 长安大学地测学院地球物理系教授, 目前主要从事地震学相关领域的研究. (E-mail: baicy@chd.edu.cn)
摘要：高频假设下的地震射线理论、以及相应的地震成像理论, 使得传统射线成像方法不可能具有较高的成像分辨率.相反, 基于有限频的射线理论则更符合实际地震的传播规律, 即地震波的走时不仅与中心射线(传统的几何射线)上的速度分布有关, 而且与中心射线一定范围(称其为第一菲涅尔体)内的速度异常分布有关.鉴于此, 本文提出了计算多震相菲涅耳体有限频地震射线的方法, 同时给出了利用多震相菲涅耳体有限频射线进行速度和反射界面同时反演成像的迭代公式.利用多震相走时资料, 分别使用传统射线走时层析成像方法和菲涅尔体有限频射线走时层析成像方法进行速度和界面的同时反演成像, 结果表明: 当射线密度较稀疏时, 无论是对速度模型的重建还是对反射界面几何形状的更新,菲涅尔体有限频射线层析成像方法均优于传统的射线走时层析成像方法, 而变频菲涅尔体有限频射线层析成像方法又优于单一频率的菲涅尔体有限频射线层析成像方法.
关键词： 第一菲涅尔体     菲涅尔体有限频射线     多震相走时     同时反演     地震射线层析成像    
Simultaneous inversion for velocity and reflector geometry based on the multi-phase fresnel volume rays
LI Xing-wang1, BAI Chao-ying1,2 , YU Zi-chao1    
1. Department of Geophysics, College of Geology Engineering and Geomatics, Chang'an University, Xi'an 710054, China;
2. Institute of Computing Geophysics, Chang'an University, Xi'an 710054, China
Abstract: Traditional ray tomography method based on high frequency assumption is unable to obtain a high resolution tomographic picture. In contrast, the finite-frequency ray theory is more suitable for real seismic propagation law, that is to say that the travel time is dependent not only on the velocity distribution along a central ray (or traditional geometric ray), but also on the velocity anomaly within a region (referred as first Fresnel Volume), which embrace the central ray. In this paper we first put forward an algorithm to calculate the multi-phase Fresnel Volume finite-frequency ray, and then give a inversion method to simultaneous invert both velocity and reflector geometry by using these multi-phase arrival time information. In synthetic example, both the traditional ray tomographic and finite-frequency ray tomographic methods are used to simultaneously update both velocity field and reflector geometry with multi-phase arrival times, the results show that the finite-frequency ray tomographic method is advantageous over the traditional ray tomographic method in terms of velocity reconstruction and reflector geometry updating when the ray density is relatively low. The finite-frequency ray tomographic method with varying frequency is again better than the finite-frequency ray tomographic method with single frequency.
Key words: first Fresnel Volume     Fresnel Volume finite frequency ray     multi-phase travel times     simultaneous inversion     seismic ray tomography    
0 引 言

层析成像始于医学领域的CT诊断技术，Aki和Lee(1976)首次将其应用于地震学的研究.经过三十多年的不断发展和完善，地震层析成像已成为目前研究地球内部结构行之有效的主要技术途径之一，出现了多种成像方法技术(成谷等，2002)，例如：体波成像、面波成像、全波形反演成像、自由振荡、有限频成像和噪声成像等.20世纪80年代提出的波动方程层析成像方法，虽然理论上具有更高的分辨率，但由于波动方程层析成像计算效率较低，反演目标函数和速度扰动之间强烈的非线性关系，且反演解对初始模型的依赖程度很高，加之实际地震波传播难以准确描述等诸多现实问题，波动方程层析成像在实际地震成像中的应用并不普遍.目前应用较广的依然是体波走时成像.除了高频假设的制约外，体波走时成像分辨率主要受射线密度或射线交错概率的影响较大.为了在现有地震台网和地震分布的情况下提高走时成像的分辨率，多震相(如：反射波、转换波)走时的引入则至关重要.

不同于射线理论的高频近似，实际地震波的频率是有限的.如：由断层滑动产生的地震波，断层破裂面越大就越容易生成较低频率的地震波，加之地球本身并非完全弹性体，地震波在传播过程中高频部分吸收、衰减严重，使得接收到的地震信号的频率局限于一定的范围.所以，仅用无限高频的观点来看待地震波具有一定的局限性.普林斯顿大学研究团队发展了可有效运用于层析成像研究的有限频理论，在有限频理论下，地震波对速度敏感的范围将不再仅仅局限于中心射线路径上，中心射线附近一定范围内速度异常的影响也将反映在走时数据上，这种影响程度可用走时敏感核来表示.Woodward(1992)分别利用Born和Rytov两种近似给出了单频地震波的走时与空间速度场分布的关系.结果发现，中心射线上的速度结构的敏感度为零，中心射线附近一定范围内的速度异常对走时的敏感度较大，而大于一定距离的速度结构的敏感度趋近于零.因此，影响地震波走时的速度结构主要集中在传统射线一定的范围内，称为第一菲涅尔体.敏感核的形状像一个空心的香蕉，而在垂直于射线路径的截面上，核的形状如同甜甜圈，所以有限频理论又称作香蕉-甜甜圈理论.鲁来玉等(2011)对浅层面波有限频灵敏度核进行了分析，刘玉柱等(2011)探讨了菲涅耳体地震层析成像的分辨率问题.关于有限频理论与传统射线理论的优劣之争，已经有了大量的研究成果(Montelli et al., 2004; Hung et al., 2004; Zhao et al., 2005; Tromp et al., 2005; 刘玉柱等，2009; 杨峰等，2010)，有关有限频层析成像的综述可参见文献(徐小明等，2009；江燕和陈晓非，2011).

敏感核的计算是一个复杂而庞大的过程，将其运用于地震层析成像会极大地增加计算成本.由于影响地震走时的速度结构主要集中在第一菲涅尔体范围内，故可通过求取第一菲涅尔体近似代替计算敏感核分布.Husen和Kissling(2001)讨论了初至波有限频射线(称其为：Fat Ray)成像、以及在地方震定位中的应用，Ke等(2007)也进行了初至波Fat Ray的讨论.本文将这一概念进一步推广，讨论利用多震相菲涅尔体有限频射线进行速度和反射界面的同时反演问题.鉴于此，本文首先讨论多震相菲涅尔体有限频射线的追踪计算问题，然后讨论利用共轭梯度法求解带约束的阻尼最小二乘问题的反演算法，给出利用多震相有限频射线进行同时反演的迭代公式，进而利用数值模拟实例验证反演算法的有效性和正确性.结果表明：无论是对速度结构还是对界面形态，有限频射线方法的结果均优于传统射线方法的结果.

1 多震相菲涅尔体有限频地震射线追踪

白超英等在改进型最短路径算法的基础上，采用分区多步计算技术，实现了2D/3D复杂层状介质中多次波射线追踪正演计算(Bai et al., 2009)，随后又将该算法推广，实现了不规则最短路径算法下多震相地震射线的追踪计算(Bai et al., 2010; 赵瑞和白超英，2010).本文在上述多震相地震射线追踪正演算法的基础上，采用一种近似方法求取多震相地震波第一菲涅尔体分布范围，给出了计算多震相有限频地震射线的方法技术.

1.1 第一菲涅尔体的定义

有限频理论下，影响地震波走时的空间速度结构主要集中在中心(或传统)射线的邻近区域，即第一菲涅尔体内，因此求取第一菲涅尔体是关键.Cerveny和Soares(1992)对第一菲涅尔体定义为

其中x为空间任意一点，s为震源，r为检波器，T为地震波周期.tsr为震源与检波器之间的最快走时，tsx为震源与点x之间的最快走时，trx为检波器与点x之间的最快走时.若点x满足(1)式，则点x属于震源s和检波器r对应的第一菲涅尔体Ω(见图 1所示).

图 1 震源s到检波器r的第一菲涅尔体示意图 Fig. 1 Diagram showing the first Fresnel Volume between source and receiver

从以上定义可以明显地看出，菲涅尔体的计算是可逆的，也就是说交换震源和检波器的位置，菲涅尔体保持不变.第一菲涅尔体的边界由以下关系确定：

1.2 菲涅尔体有限频射线追踪算法

1.2.1 传统菲涅尔体有限频射线追踪算法

追踪具有一定宽度的菲涅尔体射线，传统的方法是从第一菲涅尔体的定义出发，所求的第一菲涅尔体即为有限频射线，现对其算法做简单介绍.对于特定的震源s和检波器r，直达波菲涅尔体可通过以下方法求得：(1)以震源s为源点，计算全区所有节点的初至波走时Ts(图 2a所示)，同时也获得了检波器点的走时tsr；(2)以检波器点r为源点，反向计算全区所有节点的初至波走时Tr(图 2b所示)；(3)计算各节点上的走时差：dt=Ts+Tr-tsr；(4)满足dt≤T/2的点均位于炮(s)检(r)相应的菲涅尔体内，而其中满足dt=T/2的点的集合构成了菲涅尔体射线的边界(图 2c).图 2c给出了五种不同频率(0.5 Hz、1 Hz、2 Hz、5 Hz、20 Hz)下的菲涅尔体有限频直达波射线的边界.

图 2 传统菲涅尔体射线追踪原理，模型速度v=5.0 km/s
a)由震源点计算所得走时场;(b)由检波器点计算所得走时场;(c)不同频率地震波的菲涅尔体边界，从外到内频率分别为0.5 Hz、1 Hz、2 Hz、5 Hz、20 Hz，虚线代表传统射线路径. Fig. 2 Principle of traditional finite frequency Fresnel
Volume ray tracing，velocity model is v=5.0 km/s(a)Travel time field calculated from the source point;(b)Travel time field calculated from the receiver point and (c)The boundaries of first Fresnel Volume with varying frequencies，from inside to

outside the frequency is 20 Hz，5 Hz，2 Hz，1 Hz and 0.5 Hz，respectively. The dashed line is the traditional ray path.

上面介绍的方法原理简单、计算方便，且精度较高，但是需要计算和存储M+N(M为震源总数，N为检波器总数)次全区走时场信息，所需计算量和内存较大，加之算法自身的局限性，不易计算后续波菲涅尔体射线.下面介绍一种近似求取第一菲涅尔体半径的方法来确定菲涅尔体有限频射线.由于本文仅讨论了二维情形下的菲涅尔体，故称其为菲涅尔带更加合适.

1.2.2 第一菲涅尔带半径

由式(1)可知，对于二维情形，均匀各向同性介质中的菲涅尔带为标准椭圆，其标准椭圆方程为

根据(2)式可得

其中L为从震源到检波器的射线长度，λ为地震波波长．设l∈[0，L]为射线路经上的自然坐标，整理(3)式，并使震源位于原点(l=0)处，则射线上任意点l处第一菲涅尔带半径为

在λL的假设下，a≈L/2，b≈ Lλ 4，则第一菲涅尔带半径可近似表达为

其中v(l)为点l处的速度，f为地震波频率.

利用多震相地震射线追踪正演算法，得到射线路径和走时后，根据(7)式，求得中心射线上各点处第一菲涅尔带的半径，据此得到对应射线的近似第一菲涅尔带，即准轴菲涅尔带(Paraxial Fresnel Zone).关于射线追踪算法的详细介绍，可参阅文献(唐小平和白超英，2009).图 3给出了用精确计算法(1.2.1节)和近似计算法(7式)得到的均匀速度场(左图)和线性速度场(右图)中的直达波有限频射线、以及相应的第一菲涅尔带的边界.

图 3 使用传统方法得到的第一菲涅尔带与近似求取第一菲涅尔带半径所得到的准轴菲涅尔带对比
(a)均匀速度场(v=5.0 km/s)的结果;(b)线性速度场(v=0.06z+5.0 km/s)的结果). 地震波频率为2 Hz，图中闭合实线为精确法得到的菲涅尔带边界，短实线为近似法得到的菲涅尔带半径，虚线为中心射线. Fig. 3 The comparison of first Fresnel Zone calculated by the traditional and the approximated ways diagram
(a)uniform velocity field: v=5.0 km/s; diagram(b)linear velocity field: v=0.06z+5.0 km/s). Seismic frequency is 2Hz，in figure the closed line is the result for traditional way and the short line is the results for the approximated way. The dashed line is the central ray path.

从图 3中可以看出精确法和近似法的结果基本相同.在均匀速度模型中(图 3a)，两种算法所得到的第一菲涅尔带几乎完全一致；在线性速度模型中(图 3b)，二者差异也很小.这一结果表明使用近似法求取第一菲涅尔带半径是可行的.另外，从图 2c和图 3a-b中还可以得出这样的结论：菲涅尔带有限频射线的宽度除了与所用地震波的频率有关外(即：频率越高，菲涅尔带有限频射线的宽度越小，反之亦然；当频率趋于无穷高时则菲涅尔带射线褪变成传统的中心或几何射线)，还与速度场的分布有关(即：速度越大，菲涅尔带有限频射线的宽度越大，反之亦然).这一结论可从(7)式中得出，当地震波频率(f)一定时，则速度越大，第一菲涅尔带的半径(R)也就越大.

1.2.3 多震相菲涅尔带有限频射线追踪实例

首先讨论相对简单的一次反射或反射转换波菲涅尔带有限频射线的追踪，选用图 4所示的三层速度模型，自上而下速度分别为：5.0 km/s，6.5 km/s和8.0 km/s.炮点位置(0.0 km，0.0 km)，为了直观清晰起见，地面上只安放了两个检波器.分别计算了三种不同频率的来自上、下界面的反射波.这里震相符号中上标代表上行波，下标代表下行波，数字1、2、3表示分区号.图 5给出了相同三层速度模型中三种不同频率的多次反射或反射转换波的菲涅尔带射线路径.同样的方法可以得到更为复杂的多次波菲涅尔带有限频射线，这里就不在累述.

图 4 三层速度模型中一次反射或反射转换波有限频射线路径(左图: 反射波P1P1和P1P2P2P1; 右图: 反射转换波P1S1和P1P2S2S1). 图中自上而下地震波频率分别为2Hz，5Hz，10Hz. 中心虚线代表传统射线，短实线代表菲涅尔带半径. Fig. 4 The primary finite-frequency reflected or reflected and converted raypaths in three-layer

velocity model(Left panels: P1P1 and P1P2P2P1 primary reflections and Right panels: P1S1 and P1P2S2S1 primary reflected and converted waves). In figure from the top to the bottom the frequency is 2Hz，5Hz and 10Hz，respectively. The central dashed line is the traditional raypath and the short line is the radius of first Fresnel Zone.
图 5 三层速度模型中多次反射或反射转换波有限频射线路径(左图: 多次反射波P1P2P2P2P2P1; 右图: 多次反射转换波P1S2P2S2P2S1). 图中自上而下地震波频率分别为2 Hz，5 Hz，10 Hz. 中心虚线代表传统射线，短实线代表菲涅尔带半径 Fig. 5 The finite-frequency multiple reflections or reflections and conversions in three-layer velocity model(Left panels: P1P2P2P2P2P1 phase and Right panels: P1S2P2S2P2S1 phase). In figure from the top to the bottom the frequency is 2 Hz，5 Hz and 10 Hz，respectively. The central dashed line is the traditional raypath and the short line is the radius of first Fresnel Zone
2 多震相走时联合同时反演成像方法

2.1 方法原理

一般情况下，地下介质的速度分布和界面形态、位置均是未知的，先验信息只能给出粗略的结果，因此进行速度和界面同时反演是必需的.速度和界面同时反演是一项比较困难的工作，因为速度扰动将影响界面更新，而界面的细微变化或将引起射线路经发生较大的改变，导致获得不真实的反演结果.为了克服上述耦合关系，目前主要采用模型参数分离或模型参数归一化的方法(白超英等，2011).

我们使用另一种方法来进行速度和界面的同时反演.通常，非线性多震相走时同时反演问题，可归结为带约束的阻尼最小二乘最优化问题，其反演问题的目标函数为

式中：μ是阻尼因子； A 是雅克比矩阵；Δ m = λ1(δv1，δv2，…，δvN1);λ2(δz1，δz2，…，δzN2)T是模型参数的改变量，其中(δv1，δv2，…，δvN1)为速度模型的改变量，(δz1，δz2，…，δzN2)为离散反射界面深度的改变量，N1为未知速度节点数，N2为离散界面点数，(λ1，λ2)为两种不同量纲模型参数的归一化因子； d =[p1(δt11，δt12，…，δt1M1);p2(δt21，δt22，…，δt2M2);…;pn(δtn1，δtn2，…，δtnMn)]T是n种震相观测数据和理论数据之间的残差向量，其中：p1，p2，…，pn为反演时各震相所占权重，M1，M2，…，Mn是各震相的射线条数； a和b是反演速度模型可容许的上下边界值；e和f 为离散反射界面深度可容许的上下边界值.式(8)解的约束为

上述方程可采用共轭梯度算法求解(Bai and Greenhalgh， 2005).

2.2 Jacobi矩阵

同时反演速度和界面时，Jacobi矩阵包含两部分内容：一是走时对速度的偏导数；二是走时对反射点深度的偏导数：

式中ti是第i条射线的走时；vk为模型空间中第k个单元的速度，N为速度单元个数；zj为第j个离散反射点深度，M为离散反射点数.

走时对速度的偏导数，二维情形下，依据射线理论，某一条射线L的走时ti是对沿射线路经的速度v(x，z)倒数的线积分：

其中dl为沿射线路经的线积分元.而依据有限频理论，菲涅尔带有限频射线具有一定宽度，第一菲涅尔带内的所有空间速度结构均影响走时，相应的有限频走时方程为(Vasco et al.， 1995)

其中R为垂直于射线路经L的菲涅尔带的半径；dr为R上的面积分元；AR(l)为菲涅尔带的面积.当地震波频率趋于无穷大时，菲涅尔带退化成射线，公式(12)则退化成(11)式.

我们将模型离散为K个网格单元，第i条射线的走时为ti，对应的第一菲涅尔带的面积为Ai；第k个单元格的速度分布为vk，与第i条射线对应菲涅尔带重合部分的面积为aki，根据式(12)，可得

其中Δvk为第k个单元的速度增量，Δti为第i条射线的走时残差.

若把划分单元的网格节点作为未知量，考虑各个节点对走时的贡献大小，则有

式中vj和Δvj分别为第j个网格节点的速度及速度增量，pij为第j个节点对第i条射线走时贡献的权系数：

其中Q为包含节点j的单元格的集合，二维情形下Q中元素个数只能为1、2或4；Mki为同时位于单元格k和第i条射线对应的第一菲涅尔带的未知网格节点数.

根据式(14-15)式可得到Jacobi矩阵元素中走时关于速度的偏导数为

走时关于反射点深度的偏导数，可参阅文献(黄国娇和白超英，2010)，这里不再重复.Jacobi矩阵既包含了走时对速度的偏导数，又包含了走时对离散反射点深度的偏导数，二者量纲不同，数值大小差异较大，且后者远比前者稀疏.这种情况下，反演时主要更新速度模型，对界面的更新非常有限.为使速度和界面都具有相当程度的更新，本文采用对两种参数各自归一的方法进行归一化处理，详情可参见文献(黄国娇和白超英，2010).

3 数值模拟实例

本例中选用了一个较为复杂的层状起伏速度模型(图 6a所示)，速度模型随深度线性递增，并含有两个反射界面，上界面类似余弦函数与速度起伏重合，中心位于20 km深度，下界面则是位于36~38 km深度之间的倾斜台阶状界面.模型网格单元间距2 km×2 km，次级节点间距为0.2 km.模型顶部等间距排列11个炮点，21个检波器.反演选用水平层状速度模型作为初始模型，即：模型速度随深度线性递增(v=4.0+0.1z km/s)，初始上界面置于水平，下界面为倾斜.所用震相资料为初至P波以及分别来自上、下界面的一次反射PP波走时.为了进行对比研究，同时反演中采用了两种不同的算法，一种是传统的射线方法，另一种则是本文提出的菲涅尔带有限频射线方法.反演中两种算法采用相同的阻尼因子(0.05)和反演次数(30次)，有限频算法中采用了三种频率(4 Hz、13 Hz、25 Hz).

图 6 同时反演速度成像结果
(a)真实模型和反射界面;(b)传统射线走时反演结果;(c)-(e)分别为4 Hz、13 Hz、25 Hz有限频射线反演结果;(f)变频有限频射线反演结果). 图中白线代表界面. Fig. 6 Simultaneous inversion results for velocity field
(a)Real velocity model and reflected interfaces;(b)Results for traditional ray method;(c)-(e)Results for finite-frequency ray method with 4 Hz，13 Hz and 25 Hz frequency，respectively;(f)Finite-frequency ray method with varying frequency). In figure the white lines are the reflected interfaces.

图 6给出了相应的速度模型反演结果，尽管炮检排列较为稀疏(11个炮点，21个检波器)，但由于使用了三种震相的地震走时数据，我们依然得到了较好的成像效果.三种不同频率的有限频成像结果(图 6c-6e)均好于传统的射线类成像结果(图 6b)，特别是对于上界面附近和上、下界面之间的速度分布.就三种不同频率有限频反演结果来看，13 Hz频率的结果(图 6d)要好于其它两种频率的成像结果.另一方面，上界面以上区域内，射线密度较大，所以传统射线法和有限频射线法对这一区域内的速度重建大体相同，而在上、下界面之间的区域，由于主要是来自下界面的反射波射线，因此，射线密度相对较小，则有限频射线法的优势就比较明显.

为了进一步讨论有限频射线层析成像的特点，我们同时还进行了变频有限频射线层析成像的实验，仿照传统走时成像的自适应反演(Thurber and Eberhart-Phillips， 1999)，即先从粗网格单元开始，反演过程中逐步向细网格单元过渡，也就是先反演大尺度结构，然后反演小尺度(细微结构).在有限频射线反演中，具体来说就是先从低频(对应大尺度结构)开始，反演迭代过程中不断增大频率(对应小尺度或精细结构).本实验中频率从4 Hz开始，然后每次增加4 Hz，即：4 Hz、8 Hz、12 Hz、16 Hz、20 Hz，每个频率迭代6次(共计30次)，其结果见图 6f，总的来说，这种变频反演结果要优于单一固定频率的反演结果.这一点从图 6d和图 6f的对比就可以看出，图 6d是三种固定频率中反演结果最好的，但变频反演的结果(图 6f)则好于图 6d.从这一结果来看，采用变频反演是可行的.

为了直观评价反演效果，即速度模型的更新程度，图 7给出了反演速度与真实速度的百分比误差，即：百分比误差越小，速度场重建恢复率越高.从中可以得出这样的结论：(1)有限频射线法要优于传统的射线法；(2)但就单一有限频来讲，13 Hz的结果(图 7d)是三种频率中最好的；(3)变频有限频射线法要优于单一频率有限频射线法.

图 7 反演速度与真实速度的百分比误差
(a)真实速度模型与初始模型的百分比误差，也即扰动异常;(b)传统射线结果;(c)-(e)分别是4 Hz、13 Hz、25 Hz有限频射线结果;(f)变频有限频射线结果. Fig. 7 The percentage residual between the inverted and the real velocity model
(a)Percentage residual between the real and the initial velocity model;(b)Result for traditional ray method;(c)-(e)Results for finite-frequency ray method with respective 4 Hz，13 Hz and 25 Hz;(f)Result for finite-frequency ray method with varying frequency.

图 8给出了两种不同方法对上、下反射界面更新的细节，从图中对界面的更新情况可以看出，无论是那种单一频率有限频射线法均好于传统射线法.相比之下，13 Hz有限频射线法的结果(图 8c)是三种不同频率结果中最好的，但变频有限频射线法的界面更新结果(图 8e)要优于13 Hz单频有限频射线法的结果，尤其是对类似余弦函数的上界面的更新几乎与原始界面一样.

为了进一步了解反演解的收敛情况，表 1罗列了反演结果的误差值.从表中也可以明显的看出，无论是变频还是单一频率有限频射线法的收敛均好于传统射线法，且13 Hz有限频射线法的收敛在三种频率中是最好的，但变频有限频射线法的收敛又优于任意一个单频有限频射线法.这一结果与上述的讨论是一致的.因此，可以得出这样的结论：不管是对速度结构还是反射界面的更新，单频有限频射线的成像效果均优于传统射线的成像效果，而变频有限频射线的成像效果又好于单频有限频射线的成像效果.

图 8 界面更新结果
(a)传统射线法结果;(b)-(d)分别为使用4 Hz、13 Hz、25 Hz有限频射线法结果;(e)变频有限频射线法结果). 图中白线为真实界面，黑线为更新后的界面. Fig. 8 The updated reflected interface
(a)Result for traditional ray method and (b)-(d)Results for 4 Hz，13 Hz and 25 Hz finite-frequency ray method，respectively;(e)Result for finite-frequency ray method with varying frequency). In figure the white lines are the real interfaces and dark lines are the updated interfaces.
表 1 速度和界面同时反演结果误差 Table 1 The inverted residuals for simultaneous inversion
4 结果与讨论

地震射线密度稀疏、射线相互交叉概率较小的情况下，传统射线层析成像方法不可能具有较高的分辨率.针对这一问题，我们在多震相地震射线追踪正演算法的基础上，提出了计算具有一定宽度(此宽度与频率和速度有关)的多震相菲涅尔带有限频射线的方法，同时给出了利用多震相菲涅尔带有限频射线进行速度模型和反射界面形状同时反演的算法，并与传统的射线法进行了对比分析.当射线密度较大时，利用传统射线法和有限频射线法进行同时反演成像的结果相当，而当射线密度稀疏时，无论是对速度场的重建还是对反射界面几何形状的更新，有限频射线法要优于传统的射线法.另一方面采用变频有限频射线层析成像的物理意义是明显的，即首先进行大尺度(粗网格单元)结构重建，然后有机地过渡到小尺度(精细结构)的重建.数值模拟反演实验也验证了变频有限频射线层析成像的可行性.下一步工作的重点是将其推广至三维情形，在利用多震相走时数据的同时，考虑引入多个频率(变频)成分的地震波，从而在利用天然地震资料对地球内部精细结构进行成像时，能够有效提高成像分辨率.

致 谢 感谢审稿专家和编辑部的支持.
参考文献
[1] Aki K, Lee K W K. 1976. Determination of three-dimensional velocity anomalies under a seismic array using first P arrival times from local earthquakes: 1. A homogeneous initial model[J]. Journal of Geophysical Research, 81(23): 4381-4399.
[2] Bai C Y, Huang G J, Li Z S. 2011. Simultaneous inversion combining multiple-phase travel times within 3D complex layered media [J]. Chinese J. Geophys. (in Chinese), 54(1): 182-192.
[3] Bai C Y, Huang G J, Zhao R. 2010. 2-D/3-D irregular shortest-path ray tracing for multiple arrivals and its applications [J]. Geophys. J. Int., 183(3): 1596-1612.
[4] Bai C Y, Tang X P, Zhao R. 2009. 2-D/3-D multiply transmitted, converted and reflected arrivals in complex layered media with the modified shortest path method[J]. Geophys. J. Int., 179(1): 201-214.
[5] Bai C Y, Greenhalgh S. 2005. 3-D non-linear travel-time tomography: Imaging high contrast velocity anomalies[J]. Pure Appl. Geophys., 162(11): 2029-2049.
[6] Červeny V, Soares J P. 1992. Fresnel volume ray tracing[J]. Geophysics, 57(7): 902-915.
[7] Cheng G, Ma Z T, Geng J H, et al. 2002. A review on the growth of seismic tomography[J]. Progress in Exploration Geophysics (in Chinese), 25(3): 6-12.
[8] Huang G J, Bai C Y. 2010. Simultaneous inversion with multiple traveltimes within 2-D complex layered media[J]. Chinese J. Geophys. (in Chinese), 53(12): 2972-2981.
[9] Hung H S, Shen Y, Chiao Y L. 2004. Imaging seismic velocity structure beneath the Iceland hot spot: A finite frequency approach[J]. Journal of Geophysical Research, 109(B08305): DOI: 10.1029/2003JB002889.
[10] Husen S, Kissling E. 2001. Local earthquake tomography between rays and waves: fat ray tomography[J]. Phys. Earth Planet. Inter., 123(2-4): 127-147.
[11] Jiang Y, Chen X F. 2011. Review on the comparative study between finite-frequency tomography and ray-theoretical tomography[J]. Progress in Geophysics (in Chinese), 26(5): 1566-1575.
[12] Ke B X, Zhang J Z, Chen B F, et al. 2007. Fat-ray first arrival seismic tomography and its application[C]//77th Ann Internat Mtg Soc Expl Geophys, Expanded A. San Antonio: 2822-2826.
[13] Liu Y Z, Dong L G, Li P M, et al. 2009. Fresnel volume tomography based on the first arrival of the seismic wave[J]. Chinese J. Geophys. (in Chinese), 52(9): 2310-2320.
[14] Liu Y Z, Dong L G, Zhu J P. 2011. Research on the seismic resolution of Fresnel volume tomography[J]. Chinese J. Geophys. (in Chinese), 54(9): 2338-2347.
[15] Lu L Y, Ding Z F, He Z Q. 2011. Analysis of 3D sensitivity kernels of the finite frequency surface wave tomography in shallow subsurface[J]. Chinese J. Geophys. (in Chinese), 54(1): 55-66.
[16] Montelli R, Nolet G, Dahlen F A, et al. 2004. Finite-frequency tomography reveals a variety of plumes in the mantle[J]. Science, 303(5656): 338-343.
[17] Tang X P, Bai C Y. 2009. Multiple ray tracing in 2D layered media with the shortest path method[J]. Progress in Geophys. (in Chinese), 24(6): 2087-2096.
[18] Thurber C, Eberhart-Phillips D. 1999. Local earthquake tomography with flexible gridding[J]. Computers & Geosciences, 25(7): 809-818.
[19] Tromp J, Tape C, Liu Q Y. 2005. Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels[J]. Geophys. J. Int., 160(1): 195-216.
[20] Vasco D W, Peterson J J, Majer L E. 1995. Beyond ray tomography: Wavepaths and Fresnel volumes[J]. Geophysics, 60(6): 1790-1804.
[21] Woodward M J. 1992. Wave-equation tomography[J]. Geophysics, 57(1): 15-26.
[22] Xu X M, Shi D N, Li X F. 2009. Research progress of the finite-frequency tomography method[J]. Progress in Geophysics (in Chinese), 24(2): 432-438.
[23] Yang F, Huang J L, Yang T. 2010. Upper mantle structure beneath the Chinese capital region from teleseismic finite-frequency tomography[J]. Chinese J. Geophys. (in Chinese), 53(8): 1806-1816.
[24] Zhao L, Thomas H J, Olsen K B, et al. 2005. Frechet kernels for imaging regional earth structure based on three-dimensional reference models[J]. Bulletin of the Seismological Society of America, 95(6): 2066-2080.
[25] Zhao R, Bai C Y. 2010. Fast multiple ray tracing within complex layered media: The shortest path method based on irregular grid cells[J]. Acta Seismologica Sinica (in Chinese), 32(4): 433-444.
[26] 白超英, 黄国娇, 李忠生. 2011. 三维复杂层状介质中多震相走时联合反演成像[J]. 地球物理学报, 54(1): 182-192.
[27] 成谷, 马在田, 耿建华,等. 2002. 地震层析成像发展回顾[J]. 勘探地球物理进展, 25(3): 6-12.
[28] 黄国娇, 白超英. 2010. 二维复杂层状介质中地震多波走时联合反演成像[J]. 地球物理学报, 53(12): 2972-2981.
[29] 江燕, 陈晓非. 2011. 有限频与射线层析成像比较研究综述[J]. 地球物理学进展, 26(5): 1566-1575.
[30] 刘玉柱, 董良国, 李培明,等. 2009. 初至波菲涅尔体地震层析成像[J]. 地球物理学报, 52(9): 2310-2320.
[31] 刘玉柱, 董良国, 朱金平. 2011. 菲涅耳体地震层析成像分辨率研究[J]. 地球物理学报, 54(9): 2338-2347.
[32] 鲁来玉, 丁志峰, 何正勤. 2011. 浅层有限频率面波成像中的3D灵敏度核分析[J]. 地球物理学报, 54(1): 55-66.
[33] 唐小平, 白超英. 2009. 最短路径算法下二维层状介质中多次波追踪[J]. 地球物理学进展, 24(6): 2087-2096.
[34] 徐小明, 史大年, 李信富. 2009. 有限频层析成像方法研究进展[J]. 地球物理学进展, 24(2): 432-438.
[35] 杨峰, 黄金莉, 杨挺. 2010. 应用远震有限频率层析成像反演首都圈上地幔速度结构[J]. 地球物理学报, 53(8): 1806-1816.
[36] 赵瑞, 白超英. 2010. 复杂层状模型中多次波快速追踪: 一种基于非规则网格的最短路径算法[J]. 地震学报, 32(4): 433-444.