长安大学课程设计

课程设计

题目: 有限频射线层析成像专业：地球物理学

姓名：

学号：

指导教师：白超英教授

完成时间：2012年03月01日

有限频射线层析成像

1.引言

与探索太空一样，长久以来，了解地下世界的面貌也是人类的美好愿景，这也是地震学研究的主要任务之一。上天容易入地难！一直以来，人们主要通过地震资料，特别是天然地震资料来认识地球深部构造。地震层析成像方法广泛应用于探测地壳和地幔结构，其发展相对较早，技术方法也比较成熟。我们所要介绍的有限频射线层析成像方法就是其中的一种。

众所周知，在地震层析成像中，基于射线理论的正演解决方案具有局限性。为正确地表达地震层析成像中的正演问题，需要应用三维地震波理论。但就地震研究的规模而言，目前为止，这种方法仍是不可行的。使用地震波菲涅尔体的概念，产生了一套更加完整、更符合实际情况、更加精确的正演解决方案。一种在地震层析成像中使用有限频射线的方法，称为有限频射线层析成像。有限频理论更加符合实际情况，有限频射线层析成像产生的反演结果一定程度上优于射线层析成像产生的反演结果。

2.现有方法的缺点

自上世纪70年代以来，地震层析成像主要是根据射线理论，假设地震波为无限高频（或无限频宽），按照费马原理，即地震波沿着耗时最少的路径传播，在不同台站上观测到的到时的差异完全取决于射线路径上的速度结构。然而，实际地震波并不是无限频宽，其频率成分是有限的，影响地震波传播的主要构造也不限于中心射线，而是集中于射线邻近区域，射线理论只是无限高频条件下的一种近似。该理论只适用于地球内部速度变化平缓且速度异常体变化的尺度远大于地震波本身波长的情形，这在很大程度上限制了地震层析成像对细微速度结构的解析能力。

20世纪80年代就已经提出了不需要高频近似的波动方程层析成像方法，在理论上比走时层析成像具有更高的反演分辨率。由于波动方程层析成像计算效率非常低，其反演的目标函数和速度摄动之间表现为强烈的非线性关系，对初始模型要求很高，再加上地震子波反演困难、地震信号的信噪比较低、实际地震波传

播难以准确描述等诸多现实问题，在实际地震反演中的应用受到严重制约。

普林斯顿大学研究团队发展了可有效运用于层析成像研究的有限频理论。该理论针对地震波频带本身所具有的有限频宽特性，考虑地震波在非均匀介质中由衍射作用所产生的波前复原效应，以及不同频率的散射波经互相干涉对地震波走时的影响。在有限频理论中，相关走时完全不受射线路径上速度结构的影响，而是对环绕在射线路径周围区域的三维速度结构最为敏感。有限频理论符合地震波在介质中传播的实际情况，能更精确地反映了地下速度结构，提高了对速度异常的分辨能力。基于上述原因，发展出了有限频射线层析成像方法。

3.地震层析成像

地震层析成像主要包含两方面的内容：一是正演，即计算走时与射线追踪；二是反演求解，修正模型。不断重复上述两个步骤，直至模型满足一定的要求。

3.1 正演方法

一般使用的正演方法有两种：射线追踪和波场数值模拟。

3.1.1 射线追踪方法

射线追踪的方法种类较多。传统方法有基于初值问题的打靶法（Shooting method, i.e., Julian and Gubbins，1977）和基于边值问题的弯曲法（Bending method, i.e., Um and Thurber, 1987），但是他们不能处理介质中较强的速度变化，有时无法求出全局最小走时，计算效率较低，阴影区内无射线。随着射线追踪方法的发展，出现了直接从Huygens原理或Fermat原理出发，采用等价波前来描述地震波场特征的方法。Vidale（1988，1990）和Podvin（1990）等人则从程函方程出发，首先求出走时场分布，再用计算走时场的最速下降方向的办法，得到每一条接收点到震源的射线路径。随后，Qin（1992）等人对Vidale的方法作了改进，提出了波前扩展方法。黄联捷、李幼铭、吴如山(1992)基于Huygens 原理提出WFRT 方法，根据计算精度的要求，逐次细化网格，进而求得最小走时。Sava 和Fomel （1998）提出了HWT（huygencs wavefront tracing）法。Moser（1991）提出了根据费马原理及图论的最短路径法。

上世纪80年代后期发展起来的基于网格单元的射线追踪算法（例如：有限差分解程函方程算法，简称FD；最短路径射线追踪算法，简称SPM）由于其诸多优

点，因而倍受人们的青睐。与传统的射线追踪算法相比而言，基于网格单元的算法具有四大优点：①可利用波振面向外扩展传播的原理一次性计算出速度模型中所有网格节点的射线走时及相应的路径，并能正确的追踪检波器位于射线阴影区的衍射波射线；②算法数值计算稳定，所得到的解总是全局最佳射线路径及相应最小走时；③多炮多道接收时表现得更为高效，精度也比传统射线追踪方法要高；

④在连续介质中网格单元算法始终能找到初至波走时，而传统射线方法则只能找出唯一走时，且很难判别这个走时是属于初至波还是属于后续波（De Kool et al，2006）。目前，基于网格单元的射线追踪算法已成功用于地震定位、地震偏移成像和地震层析成像中。

3.1.2 波场模拟方法

波场模拟实质是解波动方程，基于波动方程的层析成像一般有理论地震图法和接收函数法。由于波动方程包含了地震波场的全部信息，比仅利用走时资料射线追踪层析成像更能客观地反映地下结构的信息，因此对于研究复杂条件下的各种波场最为有效，具有广阔的发展前景。但因其计算量大、计算速度慢、且易引进干扰波，目前还存在许多困难问题有待解决。

目前常用的方法有：伪谱法、有限单元法、有限差分等方法。有限差分法具有原理简单、易于编程、计算速度快、精度高等特点，然而缺点是不能保证地表起伏较大或地质结构复杂时的计算精度。有限单元法处理复杂边界时有突出优势，它可采用灵活的三角网格划分方式，较为理想地拟合起伏地表和任意形状的界面，但计算量大。伪谱法是一种整体方法，但由于它的实现过程主要是利用快速离散傅里叶变换，所以，自然地引入了周期性边界条件，使得自由边界条件很难满足，人工边界的反射也不易消除。在实际应用中往往将两种不同的数值模拟方法相结合，利用二者的优点克服缺点，提高数值模拟的精度和对模型的适应性。近年来，关于波动方程波场模拟的研究较多，如王秀明等(2003)，裴正林(2004)的研究。

3.2 反演方法

地震层析成像反演方法可以分为两类：第一类是基于算子的线性或拟线性反演方法，又称为“确定性反演方法”；另一类是基于模型的完全非线性反演方法，又称为“随机反演方法”。 ①线性反演方法，主要包括代数重建法（ART ）、同时迭代重建法（SIRT ）、及其改进型；②拟线性反演方法，主要包括高斯-牛顿法（Gauss-Newton method or Quasi-Newton method ）、阻尼最小二乘法（Damped-least squares ）、最速下降法（Steepest descent method ）、共扼梯度法（Conjugate gradients ）、亚空间法（Subspace method ）、及其改进型；③完全非线性反演方法，主要包括蒙特卡罗类方法（Monto Carlo-type method ），如：遗传算法（Genetic algorithm ）、神经网络法（Neural network method ）、模拟退火法（Simulated annealing ）、及其改进型。上述反演算法又可与空间平滑（smoothing ）、滤波（filtering ）、约束（constrained ）、规则化（regulation ）等组合形成不同的反演方法，有关局域解的反演算法及其相互关系，相关文献有Zhao and Kayal ，2000；Rawlinson and Sambirdge ，2003；Aster et al ，2005；Tarantola ，2005；Greenhalgh et al ，2006。有关全局最佳解估计法的相关文献有Sambridge and Mosegaard ，2002；Rawlinson and Sambridge ，2003。尽管全局最佳估计法理论上是解决非线性反演的最佳途径，但其惊人的计算耗时限制了实际应用，目前使用较多的还是拟线性反演方法。

4.有限频射线的概念

我们已经知道，从介质中的震源传播到接收器的地震波实际是一个体而不是一条射线，这个体称为第一菲涅尔体。地震波的第一菲涅尔体定义为产生地震能量的构造界面的核心空间区域，因此，第一菲涅尔体内每一点的散射都对接收器观察到的信号有相当的贡献，而第一菲涅尔体之外的点对走时的贡献为零。在真实地球观测到的走时反映了地震能量在第一菲涅尔带内的传播。第一菲涅尔体的定义如下：

2/T | - t + t | t sr rx sx (1) T 为地震波卓越周期，sr t 为震源与接收点之间的最短走时；sx t ，rx t 分别为菲涅尔

体内一点x 到震源和接收点之间的走时度。

有限频射线的概念与第一菲涅尔体类似，它也可以用上式定义。震源与接收器之间的走时受到它们之间的第一菲涅尔体上所有的点的影响，第一菲涅尔体之外的点对走时的贡献为零。将有限频理论应用于层析成像即有限频射线层析成像。

图1 有限频射线概念图。走时sx t ，rx t 之和小于2/T t sr +的点集定义为有限频射线(sr t 为震

源到接收点的走时，T 为主波周期)。

5. 有限频射线层析成像

简单来说，有限频射线层析成像就是在传统射线层析成像的基础上加入有限频射线的概念，使接收器记录的到时受到整个有限频射线区域而不是射线路径介质的影响。

射线理论只适用于地球内部速度变化平缓且速度异常体变化的尺度远大于地震波本身波长的情形，很大程度上限制了地震层析成像对细微速度结构的解析能力；波动方程层析成像方法条件苛刻，实际地震反演中的应用受到严重制约；有限频层析成像其敏感核的计算是一个相当复杂的过程，计算成本非常高。因此提出了有限频射线层析成像。

2

/T | - t + t | t sr rx sx =0

=| - t + t | t sr rx sx

图2简单两层速度模型中首播的有限频射线，a)为二维情况,b)为3维情况。图a)中轮廓表

示有限频射线宽度0.05。

实际情况下，接收器的分布常常是不均匀的，因此层析成像的其中一个问题就是层析成像矩阵较为稀疏，导致约束不足。射线层析成像依据沿无限窄射线路径线积分测量走时，相比较而言，有限频射线层析成像的走时是由体积分获得，相关的体由第一菲涅尔带来近似。由于有限频射线对速度模型进行了更多采样，可减少层析成像矩阵的零空间维度，削弱反演的多解性。运用有限频射线理论，有时候也可以仅仅选用部分地震数据或者使用大尺寸的单元格进行层析成像，可大幅度降低计算开销。

在此之前，已有很多学者做过有限频射线层析成像方面的研究，相对于传统射线层析成像，它具有如下优点：

①模型规则化的要求更低；

②有限频射线层析成像使用大网格单元，能够大幅度降低计算开销；

③反演结果的分辨率能够得到一定程度的提高。

图3. 有限频射线层析成像与射线层析成像反演结果对比。上图为给定速度模型，中间为射线层析成像反演结果，下图为有限频射线层析成像反演结果。

6.主要工作内容

随着社会的发展进步，人类进一步详细了解地球内部构造的愿望更加强烈。从某种角度看来，当前广泛使用的射线理论已经逐渐不能满足人们的需要。计算机及相关领域技术的不断发展，为运用新理论提供了很好的条件。相比射线理论，有限频射线理论更精确，更符合实际情况，在一定程度上提高对精细构造的分辨能力。当然，随着科学技术的日益发展，还会有更加合理、更加精确的理论运用于实际当中。

由于时间有限，在有限频射线层析成像方面，我将要尝试以下工作：

(1)运用已有的射线追踪方法，计算速度模型中所有节点的走时及追踪射线

路径，在此基础上求得每一条射线对应的有限频射线范围；

(2)将现有层析成像程序作相应的改进，使之适用于有限频射线理论；

(3)使用模型进行测试，比较其与射线层析成像结果的差异。

7.结束

完成本次课程设计工作之前，我阅读了部分相关的文献，进一步了解了地震层析成像技术的发展现状及发展方向。通过这次课程设计，我弄清了毕业设计将要进行的工作，理清了思路，为毕业设计及做好了准备。

感谢相关学者所做的研究和努力，也感谢我的导师白超英教授的悉心指导。

参考文献

Cerveny V and Soares J E P, 1992. Fresnel volume ray tracing [J]. Geophysics,57, 902-915. Dahlen F A, Hung S H. and Nolet G., 2000. Fréchet kernels for finite-frequency traveltimes–I.

Theory [J]. Geophys J Int,141,157-174.

Hole J A and Zelt B C, 1995. 3-D finte-difference reflection traveltimes [J]. Geophys J Int,121, 427-434.

Hung S H, Dahlen F A and Nolet G., 2000. Fréchet kernels for finite-frequency traveltimes – II.

Examples [J]. Geophys. J Int,141, 175-203.

Klimes L and Kvasnicka M., 1994. 3-D network ray tracing [J]. Geophys J Int,116, 726-738. Vasco D W, Peterson J E and Majer E L, 1995. Beyond ray tomography: wave paths and Fresnel volumes [J]. Geophysics, 60, 1790-1804.

Vidale J E, 1988. Finte-difference travel time calculation [J]. Bull Seism Soc Am,78, 2062-2076.

Vidale J E, 1990. Finite-difference calculations of travel times in three dimensions [J].

Geophysics,55, 521-526.

Virieux J, 1991. Fast and accurate ray tracing by Hamiltonian perturbation.[J]. J Geophys Res, 96, 579-594.

Virieux J and Farra V, 1991. Ray tracing in 3-D complex isotropic media: An analysis of the problem[J]. Geophysics, 56, 2057-2069.