地球物理学进展  2014, Vol. 29 Issue (3): 1080-1083   PDF    
引用本文
张风雪, 余大新, 张广成. 2014. 地震射线走时层析成像的简单数值例子[J]. 地球物理学进展, 29(3): 1080-1083, doi: 10.6038/pg20140310  
ZHANG Feng-xue, YU Da-xin, ZHANG Guang-cheng. 2014. A simple numerical example for seismic traveltime tomography. Progress in Geophysics, 29(3): 1080-1083, doi: 10.6038/pg20140310   
地震射线走时层析成像的简单数值例子
张风雪1,2, 余大新2, 张广成3    
1. 中国地震局地球物理研究所地震观测与地球物理成像重点实验室, 北京 100081;
2. 中国地震局地球物理研究所, 北京 100081;
3. 山东省地震工程研究院, 济南 250021
收稿日期: 2013-8-26 修回日期: 2014-1-15 .
基金项目：国家自然科学基金资助项目（41304066）和中国地震局地球物理研究所中央级公益性科研院所基本科研业务专项资助项目（DQJB13B13）联合资助.
作者简介：张风雪，男，1984年生，助理研究员，主要从事地震走时层析成像研究.（E-mail：zhangfengxue336@163.com）
摘要：地震射线走时层析成像是研究地下深部速度结构的重要地球物理方法之一，本文通过一个地震射线走时层析成像的简单数值计算实例，结合地震射线走时层析成像的基本理论和步骤，使用数值计算，具体讲述了绝对走时层析成像和相对走时层析成像的基本原理和实施流程.与抽象的概念和公式相比，一步一步的实例计算，使读者对地震射线走时层析成像的认识更具体、更直观、更易于理解，作者也希望这个数值实例对初学者能有一定的帮助和指导作用.
关键词： 层析成像     数值计算    
A simple numerical example for seismic traveltime tomography
ZHANG Feng-xue1,2, YU Da-xin2, ZHANG Guang-cheng3    
1. Key Laboratory of Seismic Observation and Geophysical Imaging, Institute of Geophysics, China Earthquake Administration, Beijing 100081, China;
2. Institute of Geophysics, China Earthquake Administration, Beijing 100081, China;
3. Shandong Institute Of Earthquake Engineering, Jinan 250021, China
Abstract: The seismic ray traveltime tomography is one of the most important geophysics methods exploring the deep underground velocity structure. In this paper, we give a talk about the basic theories and procedures of the seismic ray traveltime tomography by a simple numerical example. By the way of numerical example we explain the basic theories and procedures of the absolute and relative traveltime tomography step by step. This can give readers more concrete, intuitional and comprehensible seismic ray traveltime tomography, comparing to the abstract concepts and pure formula. And the authors hope this can also give some helps and advising to those who are beginner of the seismic traveltime tomography.
Key words: tomography     numerical    
0 引 言

地震射线走时层析成像是研究地下深部速度结构的重要地球物理方法之一(Aki et al., 1976a，1976b，1977；刘福田等，1989；Zhao et al., 1992，1994；Huang et al., 2006)，在具体的实施过程中，该方法涉及以下步骤：剖分研究区域进行网格参数化、拾取地震资料中射线的走时、追踪地震射线、构建反演方程、求解矩阵方程、评价结果的分辨效果及可靠性等.这些基本的流程步骤在相关的教材和文献中都有详细的讲解和表述(Aki et al., 1976a，1976b，1977；Zhao et al., 1992，1994；Rawlinson et al., 2003；田有等，2009；于更新，2006；王兆国，2010；张风雪，2011)，各位读者对地震射线走时层析成像的原理和具体实施过程并不陌生，本文旨在以地震射线走时层析成像理论为基础，通过一个简单的数值计算例子(张风雪，2011)，具体演示地震走时层析成像的整个流程，与抽象的公式相比，这更能使读者对地震射线走时层析成像有直观的理解，也期望这个例子对初学者能有一定的帮助作用.

1 射线走时层析成像的基本理论

在地震射线走时层析成像中，已知的观测数据是射线走时，未知参数是介质的速度.在一个速度连续变化的介质中，地震射线走时与介质速度的关系式为

式中，L是地震波的传播路径，v(x)是与空间位置有关的介质速度，t是地震射线走时.公式(1)的积分路径是L，而L又是一个依赖于速度v(x)的传播路径，所以公式(1)是一个非线性的方程，非线性意味着层析成像反演将会变得较为困难.不过我们可以采用线性化的方式将公式(1)进行线性化，在线性化前，需要有介质的参考速度，我们不妨设定参考速度为v0(x)，公式(1)可以写为

t0是地震波在参考速度模型v0(x)中沿路径L0传播所用的时间.如果在参考速度v0(x)上加一个速度扰动值δv，即v(x)=v0(x)+δv，那么地震波的传播路径和射线走时都会有一个扰动值，即L=L0+δL，t=t0+δt，代入公式(2)得

将公式(3)中的 1 v0(x)+δv 用几何级数展开

将公式(4)代入公式(3)中并略去二阶以上的小量

假定速度扰动量δv很小，参考速度模型v0(x)和扰动速度模型v(x)中的射线路径相同，则公式(5)可写为

结合公式(2)不难得出：

慢度s(x)的定义是

若对公式(8)两边求导，在参考速度模型v0(x)中可写为

如果把公式(1)和(7)中的速度变量换成慢度变量，可得到：

公式(10)和(11)分别是采用绝对走时和走时残差进行层析成像的基本公式.其中t和δt分别是绝对走时和走时残差；s(x)和δs(x)分别是慢度和慢度扰动量.如果将模型介质离散成M个块体，公式(10)和(11)可统一写为

其中m=1，2，…，M；d是观测的地震射线走时(绝对走时或走时残差)，gm是地震射线在第m个块体中的传播路径长度，xm是第m个块体的慢度参数.以上是只有一条射线的情况，如果有N条地震射线，公式(12)变为

其中n=1，2，…，N；dn是第n条地震射线的观测走时(绝对走时或走时残差)，gnm是第n条地震射线在第m个块体中的传播路径长度；xm是第m个块体的慢度参数.公式(13)的矩阵形式为

其简写形式为

其中，D 是N维的观测走时向量，元素为dn； G 是一个N×M维的矩阵，元素为gnm； X 是M维的慢度参数向量，元素为xm.

2 层析成像的计算举例

首先设定一个二维均匀模型，模型中介质的波速为1.0 km/s，此模型不妨称为模型Ⅰ.模型的尺度为2 km×2 km，见图 1，模型被图中的两条虚线均匀分为四个块体，分别是X1，X2，X3，X4，每个块体的尺度均为1 km×1 km.假定有三个激发源和三个接收点，激发源为S1，S2，S3，接收点为A，B，C，激发源和接收点的位置在图 1中分别以五角星和三角形表示.同时我们也设定了五条射线，由于是均匀介质，所以射线路径应该是连接激发源和接收点的直线段，在图中以带箭头的黑色实线表示，为了便于区别不同的射线路径，在此我们定义：S1B表示从激发源S1出发到达接收点B的射线，S1B是所设定的五条射线中的一条，其它四条射线分别为S1C，S2A，S2B，S3B.

图 1 激发源、接收点、射线路径在均匀模型中的分布 Fig. 1 The distribution of shooting points，receivers and raypaths in the homogeneity media

在均匀介质中，射线的走时可以比较容易地以解析的形式给出，即射线路径长度除以波速，五条射线的走时分别为(单位：s)：

如果将均匀模型中的四个块体给以速度扰动，加了扰动速度后四个块体的速度分别为X1:0.9 km/s；X2:1.1 km/s；X3:1.1 km/s；X4:0.9 km/s，此模型不妨称为模型Ⅱ，为了加以区别，我们将速度相对较低的X1和X4块体涂成灰色，见图 2.在扰动模型Ⅱ中，上述五条射线的路径用最小走时射线追踪方法(Moser，1991)求出的分布情况如图 2所示，这些射线的走时分别为(单位：s)：

图 2 扰动模型及扰动射线 Fig. 2 The perturbation model and perturbation raypaths

如果以均匀模型Ⅰ作为参考模型，以扰动模型Ⅱ中的走时为观测数据，可以采用绝对走时或走时残差成像反演的方法把模型Ⅱ反演出来.

2.1 采用绝对走时成像反演模型Ⅱ

以模型Ⅱ中的绝对走时为观测量进行反演，按公式(10)构建方程

上式中的x1，x2，x3，x4分别代表X1，X2，X3，X4块体中的慢度，将上式写成矩阵的形式

采用matlab求解得出

x1，x2，x3，x4是慢度值，按公式(8)转为速度则得到各个块体的速度值：X1：0.9112 km/s，X2：1.1000 km/s，X3：1.1206 km/s，X4：0.9253 km/s.

2.2 采用走时残差成像反演模型Ⅱ

以模型Ⅱ相对模型Ⅰ的走时残差为观测量进行反演.将模型Ⅱ中的射线走时减去模型Ⅰ中的射线走时得到各条射线的走时残差

按公式(11)构建方程

上式中的x1，x2，x3，x4分别代表X1，X2，X3，X4块体中的慢度扰动量，将上式写成矩阵的形式

采用matlab求解得出

x1，x2，x3，x4是慢度扰动值，按公式(9)可转为速度扰动值，将各个扰动值加到参考模型Ⅰ中则得到各个块体的速度值：X1：0.9026 km/s，X2：1.0909 km/s，X3：1.1076 km/s，X4：0.9193 km/s.

将使用绝对走时和走时残差反演的结果与模型Ⅱ的实际速度值做一个对比，如表 1，从表中我们可看出，反演求解出的各个块体的速度误差都在3%以内，块体X2中的射线覆盖交叉情况较好，无论在绝对走时反演中还是走时残差反演中，块体X2的速度都相对较好地被解出.其它块体由于射线覆盖不足或交叉不好，恢复效果相对较差.以上的例子中都只是单次反演，但在实际的层析成像中，一般都是采用多次迭代反演的模式，也就是说将前次反演所得的速度模型作为下次反演的参考模型，也就是不断更新本例中的模型Ⅰ，从理论上讲，多次迭代反演对速度的恢复效果比单次迭代会好些.

表 1 反演结果 Table 1 The inverse results
3 结 论

本文依据地震射线走时层析成像的理论，通过一个数值计算的例子具体介绍了地震射线走时层析成像的基本流程.通过这个例子，作者希望读者对地震射线走时层析成像能有一个直观的理解，也希望对层析成像的初学入门者能有一定的帮助作用.

致 谢 本研究受到国家自然科学基金资助项目(41304066)和中国地震局地球物理研究所中央级公益性科研院所基本科研业务专项资助项目(DQJB13B13)的联合资助.感谢审稿专家提出的宝贵修改意见和建议.

参考文献
[1] Aki K, Christoffersson A, Husebye E S. 1976a. Three-Dimensional seismic structure of the lithosphere under Montana Lasa[J]. Bulletin of the Seismological Society of America, 66(2): 501-524.
[2] Aki K, Lee W H K. 1976b. 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.
[3] Aki K, Christoffersson A, Husebye E S. 1977. Determination of the three-dimensional seismic structure of the lithosphere[J]. Journal of Geophysical Research, 82(2): 277-296.
[4] Huang J L, Zhao D P. 2006. High-resolution mantle tomography of China and surrounding regions[J]. Journal of Geophysical Research, 111: B09305, doi: 10.1029/2005JB004066.
[5] Liu F T, Li Q, Wu H, et al. 1989. On the tomographic inverse method used in velocity image reconstruction[J]. Chinese J. Geophys. (in Chinese), 32(1): 46-61.
[6] Moser T J. 1991. Shortest path calculation of seismic rays[J]. Geophysics, 56(1): 59-67.
[7] Rawlinson N, Sambridge M. 2003. Seismic traveltime tomography of the crust and lithosphere[J]. Advances in Geophysics, 46: 81-198, doi: 10.1016/S0065-2687(03)46002-0.
[8] Tian Y, Zhao D P, Liu C, et al. 2009. A review of body-wave tomography and its applications to studying the crust and mantle structure in China[J]. Earth Science Frontiers (in Chinese), 16(2): 347-360.
[9] Wang Z G. 2010. The tomographyic imaging of the seismological summarized information and the inversion of focal mechanism in Northeast China[Ph. D. thesis]. (in Chinese), Jilin: Jilin University.
[10] Yu G X. 2006. Finite difference simulation of wave field seismic tomography method and its applications[Ms. D. thesis]. (in Chinese), Beijing: China University of Geosciences.
[11] Zhang F X. 2011. Finite-frequency traveltime tomography by body wave and its application in the North China[Ph. D. thesis]. (in Chinese), Beijing: Institute of Geophysics, China Earthquake Administration.
[12] Zhao D P, Hasegawa A, Horiuchi S. 1992. Tomographic imaging of P and S wave velocity structure beneath Northeastern Japan[J]. Journal of Geophysical Research, 97(B13): 19909-19928.
[13] Zhao D P, Hasegawa A, Kanamori H. 1994. Deep structure of Japan subduction zone as derived from local, regional, and teleseismic events[J]. Journal of Geophysical Research, 99(B11): 22313-22329.
[14] 刘福田, 李强, 吴华,等. 1989. 用于速度图象重建的层析成象法[J]. 地球物理学报, 32(1): 46-61.
[15] 田有, 赵大鹏, 刘财,等. 2009. 体波走时层析成像方法及其在中国壳幔结构研究中的应用评述[J]. 地学前缘, 16(2): 347-360.
[16] 王兆国. 2010. 东北地区地震综合信息联合成像与震源机制反演[D]. 吉林: 吉林大学.
[17] 于更新. 2006. 有限差分模拟走时波场的地震层析成像方法与应用研究[D]. 北京: 中国地质大学(北京).
[18] 张风雪. 2011. 有限频体波走时层析成像及其在华北地区的应用[D]. 北京: 中国地震局地球物理研究所.