摘 要：本文利用二维介质中的射线追踪完成在地球内部地震波的三维的射线追踪，根据走时残差，建立地震波旅行时层析成像系统方程组。通过川滇地区宽频带地震仪所记录的远震P波走时残差，获得该区域的三维P波速度结构，同时结合当地地质构造和已有的层析成像结果，对本论文的成像结果进行分析，并对青藏高原东缘地区深部速度结构和壳慢动力学特点进行讨论与研究。

关键词：三维地震层析成像；射线追踪；地壳结构

1 前言

地球内部蕴含着丰富的资源的同时其内部的构造活动又给人类带来巨大的灾难，而利用天然地震所产生的地震波信息，已经获得了许多对地球科学产生深远影响的发现，比如地球内部成层结构的建立、深源地震的发现等[1-5]。由于地震波的速度取决于介质的密度、物质组成、弹性模量、温度等因素，因而地震波场所携带信息的差异反应了地球内部介质特性的变化，地震层析成像技术就是利用地震波的到时、波形或频散特征等，建立地球内部的三维速度结构影像。但是由于研究对象的特殊性，如震源定位不是很准确，地震波射线的透射角度有限，在非均匀介质中地震波传播的复杂性等种种原因，致使地震层析成像技术很难达到医学CT的效果。

2 IASP91地球模型中的射线追踪

如果我们把椭圆的地球近似看做球形，且内部成层状构造，可以认为地球是由具有对称性的无限多个均匀的同心圈层构成，地震波在地球介质内部传播，因而，可以利用二维介质中的射线追踪完成地震波在地球内部的射线追踪，包括地震波在圈层介质分界面上的反射波等震相的追踪。

根据IASP91地球模型的速度参数表可知，其内部速度是地球半径的函数，并且每个地层的速度梯度不等，为保证地层速度变化的连续性，网格结点之间的速度值利用三次样条插值获得，进而对介质的纵向非均匀特征进行描述。

针对上述地球模型，我们可以将地震波在地球内部传播的纵波及横波速度分别表示为

（1）

在二维介质的情况下，地震波传播的射线轨迹满足程函方程

（2）

可以推导出二维情形下射线轨迹满足的一阶微分方程组。

在数值计算中，利用数值微分可以写为

（3）

其中（x0，z0）为射线起始点E的坐标（图1），为射线的离源角，及为速度场在x、z两个坐标轴方向上的速度梯度，为给定的时间步长。射线经过时间后，射线将到达A点，其空间坐标为（x，z），离源角为。当A点位于某一界面时，需要根据斯奈尔定律来计算经过折射后的离源角。

图1 射线路径示意图

（1）当射线由震源E以离源角α0出发，按规定的步长dx逐步前进，到达A点，计算A点的坐标，以及走时TA，并判断点A是否到达或者超过第一层下界面。

（2）若A点没有达到下一层界面，则根据点所在地层的层序数，点到地心的距离R1，获取射线当前所在位置的速度值以及震相。

（3）根据层序数判断是射线是否到达地层的最底层，并对地层或块体内部进行追踪，确定射线传播的路径及旅行时间，并返回射线终点的坐标及离源角α等参数。

（4）从A点出发，继续以步长dx前进，到达B点，计算B点的坐标以及走时TB，同样进行判断B点所处位置，若依然没有到达或者超过界面，则重复以上计算，直到射线到达和超过界面为止；若射线超过界面，则采用“二分法”求射线与界面的交点，以及交点处界面倾角，并判断是否进行廻折波射线追踪。

（5）如果到达的不是目的界面，则要计算出入射线与法线的夹角、折射线与横轴的夹角等相关参数；射线到达目的界面后，则进行反射波上行射线追踪，直到追踪到地面为止。

（6）改变震源处的离源角，重复上述步骤就得到了一条条经过地下介质到达地面的射线。

3 层析成像的基本原理

地震波旅行时层析成像是典型的地球物理反演问题，依据地震波穿过介质内部所需的时间来重建介质的速度结构。

通常情况下，射线的走时t 可以写成：

（4）

其中，i是射线条数，l是射线元，v、s为分别介质的速度、慢度函数。

如果假设介质的慢度分布函数为，其中为介质的慢度初始模型，为介质的慢度扰动量，那么，根据（4）式有

（5）

从而有

（6）

所以，由上式可知，如果知道地震波的走时扰动量（走时残差），就可以研究介质中慢度（速度）的扰动分布特征。将积分离散化则层析成像系统的方程组具有如下形式

（7）

写成矩阵，有

D=CU （8）

其中为走时残差，由实际观测值得到；为介质分块中的慢度分布，为待求向量；C为不同单元格内所有线段的长度所形成的一个大型稀疏系数矩阵，由射线的传播路径确定。

（9）

4 地震层析成像的应用

青藏高原的形成原因一直备受国内外地学家的关注，尤其是青藏高原岩石圈的深部结构。2003-2004年期间，美国麻省理工大学与国土资源部成都地质矿产研究所合作，在三江地区进行了宽频带地震观测，本文利用远震P波走时残差，获得该区域的三维P波速度结构。

我们研究的区域范围为24.0°N～31.0°N，99.5°E～104.0°E，本文使用了25台宽频带地震仪所记录的地震数据，从中挑选出了信噪比较高的714个地震事件，震级大于4.5级，震中距为25.0°～180.0°，震源深度为0～641km，并且每一个远震事件至少有10个台站具有清晰的地震记录，每个记录的走时残差绝对值小于5s。本文按照上文阐述的射线追踪以及层析成像方法对研究区域宽频带地震数据记录进行地震波旅行时层析成像。

根据众多学者对三江地区的研究资料 [5-10] 可知，该区域地壳比较厚，莫霍面的深度达到了70km左右，并且该地壳有上、中、下三层结构，其厚度分别为10km、25km、35km，所以本文在IASP91地球模型基础上稍微做了修改，修改后的IASP91地球模型参数我们称为初始模型。网格在经度方向和纬度方向上均以1.0°间距作水平方向的划分，在深度方向上以10、20、30、40、50、65、80、100、120、150、200、250km进行划分。

我们对三江地区观测数据进行层析成像，则由0～125 km深度处的速度扰动水平剖面图中可以看出，川西高原地壳和上地幔的速度存在明显的横向不均匀性，川中块体的速度相对较高，而川青块体以及川滇菱形块体的速度相对较低，并且川滇菱形块体相对较低的速度区域呈南北向分布；大型断裂带两侧存在明显的速度差异，如龙门山断裂带、安宁河断裂带、金沙江断裂带等为研究区域中几大块体的速度分界线。

纵观图2的4个深度剖面，我们发现，在金沙江断裂带和安宁河断裂带之间的研究区域，随着深度的增加，特别是深度150～250km区间，低速度区域的面积逐渐增大。在川滇菱形块体东侧的低速区域，已穿过其速度分界线龙门山断裂、安宁河断裂带延伸至四川盆地内部。根据曾融生、Blackman、Sliver 等[11-14]提出青藏高原物质有横向流动的可能性，，推断出现这一现象的原因可能为：印度板块向青藏高原的下部俯冲，使青藏高原的地壳抬升，在巨大的挤压下，增多的地壳物质向青藏高原东部地区的上地幔流动。其中东流的物质在龙门山断裂带附近遇到了四川盆地的阻挡，部分改向东南方向流动，但是仍有部分向下侵入到四川盆地的西南区域，致使区域 P 波速度较低。并且在深部断裂带呈现负的速度异常，更有助于地壳块体沿断裂的侧向挤出。

图3 P 波震相速度扰动水平剖面

从102°E经度P波震相的速度扰动垂直剖面（图4）可以看出，该剖面自南往北依次穿过滇西块体、滇中块体、川西块体、川东南块体、川青块体。在0～150 km深度上，23°N～31°N范围内的地壳和上地幔速度存在明显的不均匀性，低速和高速相间。产生这一现象的原因可能是印度板块在与欧亚板块碰撞、俯冲、挤压过程中，由于地壳应力场的非均匀性和鲜水河断裂带的存在，在高速的四川盆地边缘发生了斜上方的逆冲运动造成的。在深度150km以下，该剖面的速度值逐渐变小，且低速度区域的面积逐渐增大。这一结果也在一定程度上验证了川西高原在下地壳至上地幔的范围内存在较软物质的可能性。

图4 102°E经度P 波震相速度扰动垂直剖面

图5 29°N纬度P 波震相速度扰动垂直剖面

从29°N纬度P波震相的速度扰动垂直剖面（图5）可以看出，在0～100 km深度范围内，地壳和上地幔的速度存在横向的不均匀性。在川滇菱形块体内101°E～102°E范围的80km深度以上，该处存在一个高速区域。100km以下，川滇菱形块体低速度区域的面积逐渐增大。并且低速体在西部以金沙江断裂为界，东部以小江断裂为界，据此我们推测，在青藏高原受到印度板块NNE方向的挤压后，川滇地区受到连带影响，由于断裂带的岩石张力的破碎，使四川盆地沿着小江断裂逆向俯冲（图5箭头所示）。

综上所述，我们推测：印度板块以NNE方向运动与欧亚板块碰撞，并向青藏高原下部俯冲，致使青藏高原地壳增厚，并向东强烈挤压，由于地壳应力场的非均匀性和断裂带的存在，在这种挤压下，川滇地区的地壳发生变形，区内的次级块体及边界产生逆冲或者扭转，对川滇地区的地质构造造成影响。由研究区域大面积的低速体的存在，推测该现象可能是由地壳和上地幔存在高导层、高热流值造成的。

纵观整个研究区域，并结合前人的研究成果，我们得出川滇地区地壳结构的总体特征是：在上地壳速度异常分布中，存在明显的横向不均匀性，但是总体上川滇地区地壳、上地幔的速度均速度较低；龙门山断裂带、安宁河断裂带，金沙江断裂带在地壳和上地幔一定深度内的速度异常中显示出构造分界特征，并且深部断裂带呈现负速度异常；地壳厚度变化剧烈，地壳和上地幔存在高导层、高热流值。

5 存在问题

地震波旅行时层析成像包括数据处理，正演模拟和数据反演等环节，而每个环节都有可能导致最终结果的不确定。该方法在应用中主要存在问题具体如下：

（1）反演的数据量不足。由于台站的数量有限，参与层析成像反演的数据量不足，从而导致成像结果的分辨率不高。

（2）震源定位不准确。由于台网的稀疏，使得地震的定位存在很大的误差，尤其是震源的深度，那么对地震层析成像反演结果有很大影响的的旅行时就会变的不准确

（3）初始参考模型引起的误差。在进行射线追踪时，由于不知道真实的射线路径，我们首先引入了参考模型，然后通过逐步迭代对参数进行修正。如果成像反演时利用了折射或者反射震相，则初始参考模型中速度界面信息的准确性尤为重要，因为速度间断面对这些震相的射线传播路径影响很大，所以可靠的速度间断面信息对地震层析成像反演是至关重要的。

参考文献

[1] Wang Y X， Jiang M， Nábelek J，et al. The new insight of the velocity structure beneath the Indian-Eurasian continental collision zone. 2006， AGU West. Pac. Geophys. Meet， Suppl.， Abstract S45A～04.

[2] Zhao Dapeng， Akira Hasegawa. P wave tomographic imaging of the crust and upper mantle

beneath the Japan islands [J]. J. Geophys. Res.， 1993， 98（B3）： 4333～4353

[3] 贺世杰，郭峰. 地幔柱的识别和演化综合述评 [J]. 地球科学进展，2003，18（3）：433～439

[4] 刘建华，刘福田，吴华，等. 中国南北带地壳和上地幔的三维速度图像 [J]. 地球物理学报，1989，32（2）：143～151

[5] 王有学，姜枚，熊盛青，等. 西昆仑岩石圈的拆沉作用及其深部构造含义—地震层

析成像及航磁异常证据 [J]. 中国地质，2006，33（2）：299～308

[6] 姚振兴，李白基，梁尚鸿，等. 青藏高原地区瑞利波群速度和地壳构造[J]. 地球物理学报，1981，24（3）：287～29

[7] 丁志峰，何正勤，吴建平，等. 青藏高原地震波三维速度结构的研究[J]. 中国地震， 2001，17（2）：202～209

[8] 钱辉，姜枚， Chen Wangping，等. 青藏高原吉隆—鲁谷（Hi-Climb）层析成像与引藏碰撞的消减作用[J]. 地球物理学报，2007，50（5）：1427～1436

[9] 吴庆举，曹融生，赵文津. 喜马拉雅-青藏高原的上地幔倾斜构造与陆-陆碰撞过程[J]. 中国科学D辑，2004，34（10）：919～925

[10] 贺日政，赵大鹏，高锐，等. 西昆仑造山带下岩石圈地幔速度结构[J]. 地球物理学报，2006，49（3）：778～787

[11] 曾融生，朱介寿，周兵，等. 青藏高原及其东部邻区的三维地震波速度结构与大陆碰撞模型 [J]. 地震学报，1992，11（S）：520～533

[12] 曾融生，丁志峰， 吴庆举. 青藏高原岩石圈构造及动力学过程研究[J]. 地球物理学报， 1994，37（S）：99～116

[13] Blackman D K. Seicmic anisotrophy in the mantel beneath an oceanic spreading center[J]. Nature，1993， 366： 675～677

[14] Sliver P G， Chan W W. Implication for continental structure and evolution from seismic anisotropy [J]. Nature， 1998， 335（1）： 34～39