西藏西部地区的远震P波层析成像研究＊

樊 杰 杨文采

（中国杭州 310027 浙江大学地球科学学院）

引言

青藏高原是世界上面积最大、海拔最高、形成时代最晚的高原，它形成于印度板块与亚欧板块的碰撞过程中，该碰撞始于距今50 Ma前.持续的碰撞使得青藏高原现今的地壳仍在增厚，并且在其边缘发生了很多强震（Yin，Harrison，2000）.青藏高原独特的地理位置和地质背景使其成为板块碰撞及俯冲相关的地球动力学研究的天然实验室.

前人已运用体波成像、面波成像、背景噪声、接收函数等多种方法在青藏高原开展了大量的深部结构探测研究，地球动力学、高原内部壳幔结构地区差异及其两侧板块俯冲状态等方面取得了丰富的成果.

关于青藏高原地球动力学方面的研究，Molnar和Tapponnier （1975）认为印度板块的俯冲可能贯穿整个西藏地区地壳，但还未扩展至300 km或400 km深度以上.Razi等（2014）采用体波层析成像得到了西藏西部地区上地幔的速度结构，其结果显示该区上地幔存在狭长的高速异常，这一高速异常可能是印度板块俯冲到青藏高原下方之后地壳底部物质榴辉岩化所致.Bao等（2015）关于印度板块撕裂的研究表明印度板块向北俯冲的过程中，青藏高原南部东西两边的俯冲角度较小，中间位置的俯冲角度较大，因而推断其在俯冲过程中发生了撕裂，这与喜马拉雅地区地面裂谷的形成具有一定关联.Chen等（2015）提出板块撕裂或停滞模型，该模型指出印度板块存在不同的俯冲角度，可以很好地解释拉萨地块东部所布设台阵相对西侧具有明显的SKS到时滞后现象.Peng等（2016）关于喜马拉雅造山带东构造结地区的远震P波层析成像结果表明：向北俯冲的印度板块在西侧比较平缓，东侧比较陡峭，板块撕裂发生在平缓与陡峭区域的交界处，从而导致软流层物质上涌，上覆大陆受到底部侵蚀；向北俯冲的印度板块呈小角度的近水平俯冲，向东北俯冲的印度板块以大角度俯冲进入缅甸板块，这种俯冲会加速板块的水平扩张，并导致最终的撕裂.Wang等（2019）提出印度板块先向北俯冲，然后在东西两侧受一定水平推力的影响发生拆离，掉入上地幔，从而形成喜马拉雅弧形造山带.综上，板块撕裂按照延伸方向分为两种方式：一种是水平撕裂，其延伸方向垂直于碰撞方向，通常在青藏高原东部由于重力失稳而导致板块拆沉或断裂（Rosenbaumet al，2008）；另一种是垂直撕裂，撕裂方向近乎平行于延伸方向，在青藏高原的中部和西部通常会由于不同的俯冲速率造成剪切撕裂从而导致板块拆离或板块破裂（Huanget al，2010）.

目前对于青藏高原地球动力学的解释主要分为以下模型：高原隆升模型，该模型认为下地壳和上地幔流导致了高原隆升（Deweyet al，1988；Molnaret al，1993）；变形增厚模型，该模型认为该地区岩石圈本身在新生代发生塑性缩短增厚，而且增厚的岩石圈地幔在重力和热对流作用下发生剥离（England，Houseman，1989）；陆内俯冲模型，该模型认为青藏高原下方发生了亚洲岩石圈的陆内俯冲（Meyeret al，1998）；印度板块俯冲模型（Kosarevet al，1999；Tapponnieret al，2001；周华伟等，2002；Baoet al，2015；Li，Song，2018）等.

关于青藏高原壳幔结构的地区差异，从多年的地球物理研究（Chen，Molnar，1981；Barazangi，Ni，1982；Sandvolet al，1997；Chen，Yang，2004；Hearnet al，2004）中得到的一个最主要的特征是尽管青藏高原的地形开阔平坦，但其不同地区的深部壳幔结构存在差异，这导致不同地区深部的物理性质有一定差异，南北差异主要表现在：北部上地幔P波和S波的速度较低，而南部速度较高（Chen，Molnar，1981）；Sn在北部传播比较慢，而在南部传播较快（Barazangi，Ni，1982）；横波分裂在北部地区较明显，而在南部地区相对较弱（Sandvolet al，1997）；Pn波速度在北部地区较低，而在南部地区较高（Hearnet al，2004）.南部与北部的物理参数差异说明了壳幔结构的南北差异. 东西差异主要表现在：拉萨地块东西向展现出明显的SKS横波分裂横向变化，在拉萨地块东侧SKS到时滞后很多，可达1.0—2.0 s （Chenet al，2015），这表明东侧下覆上地幔为低速；在西藏地区的东南部，印度板块呈近水平俯冲状态，同时伴有水平方向的撕裂，至少被撕裂成两部分；而在印度板块俯冲的最西侧，俯冲板块发生断裂并拆沉至软流圈内部（Ceylanet al，2012）.印度板块东西部不同的俯冲状态导致东西壳幔结构的差异.

关于青藏高原南北两侧板块俯冲状态的讨论中，Kind等（2002）的接收函数研究结果显示青藏高原北缘存在转换界面（33°N—37°N），此转换界面被认为是亚欧板块向南俯冲的证据.然而，Liang等（2012）利用布设在西藏地区内部的流动台站得到的地震数据进行地震体波成像，其研究结果显示藏北地区上地幔由非常均匀的低速层组成，不存在亚欧板块向南俯冲的痕迹.Riaz等（2019）的研究结果也显示西藏东北部的主要构造活动是由于印度板块对亚欧大陆的挤压造成南北方向缩短和东西方向伸长所致，并且由持续挤压造成的该区域物质流出所形成的转换断层可能与2017年8月九寨沟地震的发生有关.

本研究所关注的西藏西部地区主要包括拉萨地块和羌塘地块（图1）.拉萨地块处于青藏高原最南端，南界为雅江缝合带，该地块东西长约2 000 km，南北宽约300 km，北部隔班公—怒江缝合带与羌塘地块相接，羌塘地块北界为金沙缝合带，该地块在青藏高原中部宽500—600 km，而在东西两侧宽150 km （Yin，Harrison，2000）.拉萨地块和羌塘地块在晚古生代之前曾是冈瓦纳大陆的一部分（Kappet al，2003，2007），早中生代从冈瓦纳大陆分离，开始向北漂移；晚侏罗纪拉萨地块与羌塘地块发生碰撞（Yin，Harrison，2000；Kappet al，2003）；白垩纪之后，由于两地块之间的碰撞，拉萨地块的地壳部分缩短增厚，表面隆起.

图1 西藏西部地区构造背景及本文所用台阵分布Fig.1 Tectonic settings of western Tibet and the seismic arrays used in this study

本文的研究区域聚焦西藏西部地区，在该区域有两组研究人员的研究结果值得关注.Razi等（2014）利用区域体波层析成像对西藏西部地区120 km深度以浅的壳幔速度结构进行了研究，其结果表明：在中地壳深度，青藏高原内部呈现为明显的低速异常，而在喜马拉雅山地区呈现为明显的高速异常，雅江缝合带则为低速异常南部的分界线.他们的剖面结果显示：一高速体以40°向北俯冲，最终俯冲至34°N，高低速异常边界比较明显；另一高速体以40°向东俯冲，最终俯冲至82°E，根据该高速体所处的位置推测其可能为撕裂的印度板块.Li和Song （2018）通过Pn波和S波层析成像对西藏西部地区的地下结构进行研究，其研究结果显示：印度板块在上地幔被分成四部分，且撕裂的印度板块与青藏高原接触地区的地震活动比较频繁；在上地幔130 km的深度上，S波速度分布显示印度板块撕裂更加明显，在撕裂的印度板块与青藏高原接触的地方应力比较集中.

确定西藏西部地区的深部结构对于研究印度板块的俯冲状态（包括俯冲角度和向北俯冲范围）和青藏高原的地质构造演化均具有极其重要的意义.本文研究区域为西藏西部地区，尽管前人对该地区的浅部结构已经有所研究，但是若要更清楚地分析印度板块的俯冲状态以及西藏西部地区的一些深部构造特征，仅有浅部结果远不够.此外，西藏西部的地理位置特殊且自然环境恶劣，导致该地区的地震数据一直缺失，对前期研究结果造成一定的影响.另外，尽管前人对于青藏高原整体大背景进行了较为深入的研究，但是对于青藏高原西部深部结构的研究结果不够深入且目前仍存在较大争议.为此，本研究拟利用大量西藏西部台阵数据，采用快速行进法（fast marching teleseismic tomography，缩写为FMTT）正演计算理论走时，并采用子空间迭代反演算法得到西藏西部地区更高分辨率的深部速度相对分布，以期更加清楚地认识青藏高原地区的地球动力学背景，了解印度板块向北俯冲的范围及形态，推断青藏高原内部是否存在地壳管道流，并以此来推测青藏高原未来的演化趋势.

1 数据与方法

本文所采用的数据分别包括来自于地震学联合研究机构（Incorporated Research Institutions for Seismology，缩写为IRIS）的YT （2001年9月到2001年10月），Y2 （2007年7月到2010年12月）和XF （2002年9月到2005年8月）等三个台阵，以及国际地震中心（International Seismological Centre，缩写为ISC）网站在印度、尼泊尔地区的一些地震数据（2007年1月到2009年6月）.本文所选用的这四个台阵基本上覆盖了整个西藏西部地区（图1），确保能够获得分辨率更高的速度结构.数据预处理时，按照以下条件来筛选地震事件：震中距为30°—90°，采用0.02—1 Hz的频段进行滤波，信噪比高于8，同一个地震事件能够被不少于5个台站所记录.本研究共下载了1万5 103个地震事件，经过预处理最终筛选出2 191个地震事件，342个台站，4万4 453个P波震相.

本研究所采用的地震事件分布如图2，可以看到以研究区域为中心，地震事件在各个方位覆盖完整，这样可以综合反映各个方向的地震事件对于速度相对分布的影响.

图2 反演所用2191个地震事件的分布Fig.2 Distribution of 2 191 seismic events used in this study

图3给出了经过预处理之后地震事件的波形，可见P波震相起跳明显，易于识别，说明预处理所采用的参数合适.

图3 2007年7月18日18点01分25秒地震事件经预处理之后的波形展示Fig.3 Seismic waveforms of an event occurred at 18:01:25 on July 18，2007 after pre-processing

得到适于拾取P波初至的地震波形之后我们进行到时残差拾取，首先选用AK135模型利用快速行进法计算理论到时.快速行进法是一种射线追踪方法，可用来求解程函方程，进而以此来计算地震波走时场和射线路径（李天扬等，2021）.直角坐标系下的二维程函方程可以表示为

式中，t为地震波走时，s为地震波在介质中传播的慢度.解程函方程就是求解上式中的t，采用迎风差分算法将上式简化为

通过求解式（2）即可得到走时t.

由于所采用的地震事件均为远震，对于同一个地震事件而言，每一个台站记录到的波形基本一致，因此可以采用互相关方法来确定到时残差，此处采用自适应叠加方法（Rawlinson，Sambridge，2004）.为了平衡研究区域之外地层的横向不均匀性、震源定位不准确等因素对于层析成像结果的影响，本文采用相对到时残差，即用每个事件所对应每个台站的到时残差减去接收到该事件的所有台站的残差平均值，以期通过此方法来去除上述影响.图4给出了采用自适应叠加方法得到的叠加前后的波形对比，可见10次迭代后波形起跳对齐较好，有利于提高到时拾取效率，节约到时残差的计算时间.

图4 自适应叠加方法效果示意图（a）基于AK135模型的初始对齐数据；（b）10次迭代之后的对齐数据Fig.4 Adaptive stacking example（a）Alignment of initial traces given by AK135 model；（b）Alignment of final traces achieved after ten iterations

2 反演参数选取与分辨率测试

反演前需要选取最佳的阻尼因子和光滑因子参数，选取的参数既要确保结果模型与初始模型相差不大，又要确保结果模型适用于实际地震数据.图5为选取两个参数时所使用的折中曲线（Rawlinsonet al，2006）.首先按照经验先固定光滑因子保持20不变，在0—500之间变换阻尼因子，绘制数据方差关于模型方差的折中曲线，从曲线中选取最佳阻尼因子为20；保持最佳阻尼因子20不变，在0—500之间变换光滑因子，绘制数据方差关于模型粗糙度的折中曲线，从曲线中选取最佳光滑因子为20，与最初预估的光滑因子基本一致，因此最终选择的最佳阻尼因子和最佳光滑因子均为20.

图5 最优阻尼因子ε和光滑因子η的估计Fig.5 Trade-off curves of optimum smoothing weight factor η and damping parameter ε

由于我们研究的是上地幔速度结构，地壳部分选用不同的速度模型，可能会对远震层析成像结果产生一定的影响，因此本文分别选用CRUST1.0和CRUST2.0两种地壳模型来进行层析成像，结果如图6所示.可见，采用两种模型除了对速度相对值有一定的影响之外，并未对研究结果造成较大不同，因此可以任意选择其中一种地壳模型.鉴于CRUST1.0精度更高，本研究选用CRUST1.0模型建立地壳部分的速度结构模型.

图6 基于CRUST1.0 （a）和CRUST2.0模型（b）获得的西藏西部地区在100 km和250 km深度处的地壳层析成像结果（深度标于子图的左下角，下同）Fig.6 Tomography results of the crust in western Tibet at the depth of 100 km and 250 km based on CRUST1.0 （a） and CRUST2.0 （b） models. The depth is labeled at the lower-left corner of each subfigure，the same below

由反演前后得到的到时残差分布对比（图7）可以看出，采用四个台阵进行六次迭代之后残差分布呈现向中心收敛的趋势，并且计算结果显示，残差方差从1.103下降至0.551，整体下降了50%左右，表明本文的反演方法收敛效果较好，反演方法有效.

图7 初始模型（a）和最终模型（b）的相对到时残差分布图Fig.7 Distribution of relative arrival time residuals for the initial（a）and final（b）models

为了验证本研究反演结果的分辨率，进行了西藏西部层析成像的棋盘格测试（图8），所采用的棋盘网格大小为1° × 1°，相邻两个网格间隔为1°.由图8可以看到：对于100 km以上深度即中上地壳深度，棋盘格测试仅可恢复出台站下方的异常，150—410 km深度范围内研究区域内的棋盘格测试均恢复较好.这表明本文给出的西藏西部地区台站下方下地壳和整个研究区域上地幔的远震层析成像结果较为可靠.

图8 西藏西部层析成像棋盘格测试恢复结果（左上角为对比模型）Fig.8 Checkerboard tests in this study （The compared model is on the upper-left）

3 层析成像结果

图9为研究区域的P波层析成像结果.50 km深度的层析成像结果（图9a）显示青藏高原内部出现大范围低速区，表明青藏高原内部大部分区域的下地壳以低速为主，这可能意味着下地壳流的存在（薛光琦等，2006）.Copley和McKenzie （2007）指出西藏地壳的流变状态，特别是青藏高原中下地壳存在的管道流导致了青藏高原的地壳物质向外流动.Leech （2008）基于大地电磁资料的研究结果并结合花岗岩火山年代认为，喀喇昆仑断层可能是一道阻碍青藏高原内部管道流流出的屏障，而本文结果显示50—150 km深度的喀喇昆仑断层以西低速异常很弱，这恰好验证了该观点.在青藏高原西南部边缘50—150 km深度范围内存在一些高低速异常的相间分布，可能与印度板块的撕裂相关，高速区（图9b红色箭头所示）对应于撕裂的印度板块.且本文在150 km深度处层析成像结果（图9b）所显示的高速体被低速体所截断部位与Li和Song （2018）所得印度板块被撕裂的位置基本吻合，印证了此猜测.将低速区的位置与发育在该区的新生代裂谷系相比，可以看到83°E和85°E两条裂谷系均与本研究150 km深度处的低速异常区相吻合.杨文采等（2019）关于区域重力场的三维密度扰动成像和地震层析成像研究结果表明西藏地区的新生代裂谷系与上中地壳低密度带、下地壳低波速的物质蠕动流有关.基于此，本文认为该关联性可以延伸至上地幔深度.有趣的是，周华伟等（2002）也通过高分辨率全球层析成像模型发现雅江缝合带以北的85°E—93°E之间存在大范围的低速异常区，该低速异常被认为是地幔物质经俯冲板块脆弱部位上隆的证据.

图9 西藏西部不同深度（a—d）处的层析成像结果红色箭头代表俯冲过程中撕裂的印度板块Fig.9 Tomography results of western Tibet at different depths （a−d）The red arrows represent the position of the torn Indian Plate

上地幔150—400 km深度范围内，青藏高原内部出现一些局部高速异常，这些高速异常大多分布于研究区域的西南部，可能与印度板块俯冲过程中发生的拆沉有关，因为学界普遍认为印度板块密度较大，在相对速度分布中呈现高速异常.但也有可能是印度板块在俯冲过程中发生撕裂之后与青藏高原底部接触的部位发生榴辉岩化，导致接触部位呈现高速异常（Náběleket al，2009；Raziet al，2014）.武振波等（2020）基于成像结果和岩石学研究成果证明藏南地块下方自西向东均存在俯冲印度板块的榴辉岩化现象，且印度板块在80°E这一经度下，俯冲板块前缘到达班公湖—怒江缝合带，向东递减至拉萨地块中部.从本文150—400 km深度的层析成像结果（图9b，c，d）均可看出，印度板块向北的俯冲范围从西到东逐渐减小，这一点与Kumar等（2006）关于西藏西部S−P转换震相结果以及Liang和Song （2006）的体波层析成像结果相一致.Li等（2008）的研究也表明印度岩石圈向北俯冲在西部贯穿了喜马拉雅山脉甚至整个青藏高原，而在东部的俯冲甚至未到达雅江缝合带.这些研究均表明印度板块在青藏高原西部俯冲年代较早或俯冲速率相对较快.

关于印度板块与亚欧板块碰撞后印度板块究竟俯冲至整个青藏高原下方还是俯冲到一半发生拆沉进入地幔这个问题长期以来争论不休.Kosarev等 （1999）的研究表明，印度板块确实俯冲到了整个青藏高原下方，然而喜马拉雅—青藏大陆岩石层地震探测（Himalayan-Tibetan continental lithosphere during mountain building，缩写为Hi-Climb）项目采用远震和重力数据联合反演的研究结果表明，印度岩石圈在青藏高原中部近乎垂直地延伸至雅江缝合带和班公—怒江缝合带下方400 km处（Basuyauet al，2013）.本文的成像结果只聚焦于西藏西部，50 km深度以下本文结果所呈现的青藏高原西部持续的高速异常说明印度板块应该不是垂直俯冲，而是近水平俯冲，且高速异常体向北至多延伸至34°N （图9b，c，d），即班公—怒江缝合带以南，这说明印度板块向北俯冲最远到达班公—怒江缝合带，并未贯穿整个青藏高原内部，该俯冲范围与Razi等（2014）利用近震得到的西藏西部地区的层析成像结果一致.Kind等（2002）利用P波接收函数方法对垂直于喜马拉雅山脉走向剖面上所获取的地震数据进行研究，其结果显示在青藏高原南部下地壳和莫霍面存在两个Ps转换震相，表明印度板块俯冲跨越了雅江缝合带并接近班公—怒江缝合带（赵志丹等，2003），本文通过上地幔高速异常的范围最终确定西藏西部地区印度板块的俯冲最北侧到达班公—怒江缝合带.与此不同的是，郑洪伟等（2007）的层析成像研究表明沿88°E剖面位置，印度板块的俯冲前缘已经到达羌塘地块的中部，之后进入上地幔深处，这种差异有待更深入的研究.吕庆田等（1998）根据层析成像结果推测印度板块仅俯冲到雅江缝合带以南，而高原腹地的地壳缩短增厚是通过陆内变形完成，本文认为得出这种结论的原因在于该研究所采用的台站基本都位于青藏高原中部，而本文研究中靠近88°E的150 km深度也未观测到由于印度板块俯冲所呈现的高速异常（图9b）.

基于以上结果及讨论，本研究提出印度板块在西藏西部地区的俯冲机制：印度板块在青藏高原内部近水平向北俯冲，西部向北俯冲范围相对东部较大，但最远仅俯冲至34°N附近.且俯冲过程中发生了板块撕裂现象，被撕裂的印度板块拆沉进入上地幔，而撕裂间隙由于应力释放影响了西藏西部地区新生代裂谷的形成.

4 讨论与结论

本文首次将快速行进法应用于西藏西部地区，得到了分辨率较高且较为可靠的远震层析成像结果.根据远震P波层析成像结果，青藏高原内部下地壳广泛的低速异常暗示该地区存在下地壳管道流，青藏高原内部上地幔顶部表现出一定范围的低速异常，靠近喜马拉雅山地区表现为大范围的高速异常.青藏高原内部上地幔明显的局部高速异常可能是由于印度板块俯冲过程中拆沉或榴辉岩化所致.在上地幔深度，印度板块向青藏高原俯冲过程中出现了板块撕裂现象，在速度相对分布上表现为多个高速条带，撕裂后各部分的俯冲范围各不相同，自西向东俯冲至青藏高原内部的范围逐渐减小.本文仅对西藏西部地区大尺度的地球动力学背景进行了一定的分析，后续若要进行局部细致的地球动力学研究，需增加该区域地震台站的密度.