印度和欧亚板块在50 Ma发生汇聚碰撞(Molnar & Tapponnier，1975; Rowley，1996)，造就了地球上最宏伟的高原——“青藏高原”，为我们研究陆陆碰撞过程和高原隆升机制提供了最佳实验场.目前，青藏高原西部主要的地球物理探测剖面如图 1a所示，包括中法合作探测剖面(Wittlinger et al.，2004)、“Hi-CLIMB”计划(Nábělek et al.，2009)、“ANTILOPE”计划(Zhao et al.，2010，2014)及南侧印度平原两条剖面(Caldwell et al.，2013; Rai et al.，2006).

图 1

Fig. 1

图 1 研究区域和宽频带流动台站位置 (a) 红色三角形为TW-80试验宽频带地震台站，不同颜色菱形代表其他研究者布置的宽频带台站; (b) 紫色圆圈代表1900—2015年发生的MS>4.0 级以上地震，右上角红色圆圈是本文所用地震事件. HB：喜马拉雅块体；LB：拉萨块体；QB：羌塘块体；TSH：甜水海块体；SG：松潘—甘孜块体；Qaidm：柴达木盆地；Tarim：塔里木盆地；MBT：主边界逆冲断裂；MCT：主中央逆冲断裂；IYS：雅鲁藏布江缝 合带；BNS：班公湖—怒江缝合带；JS：金沙江缝合带；AKMS：阿尼玛卿缝合带；LMF：龙木错断裂；DWT：多玛—乌江逆冲断裂；MT：曼 冬错北逆冲断裂；SF：狮泉河逆冲断裂；KF：喀喇昆仑断裂. Zhada：札达县；Shiquanhe：狮泉河地区； Rutog：日土县；Domar：多玛乡. Fig. 1 Study area and passive-source broad-band seismic stations (a) Red triangles denote broad-band seismic stations used in this study, the other coloured diamonds denote previous broad-band seismic stations; (b) The purple circles represent the seismic events MS>4.0 between 1900 and 2015 year. The red circles in the right-corner denote the seismic events used to compute the receiver functions. HB: Himalaya block; LB: Lhasa block; QB: Qiangtang block; TSH: Tianshuihai block; SG: Songpan-Garzê block; Qaidm: Qaidm Basin; Tarim: Tarim Basin; MBT: Main Boundary Thrust of the Himalayan system; MCT: Main Central Thrust of the Himalayan system; IYS: Yarlung Zangbo River suture; BNS: Bangong-Nujiang suture; JS：Jinsha River suture；AKMS: A′nyêmaqên-Kunlun-Mustagh suture; LMF: Longmu Co fault; DWT: Domar-Wujiang thrust; MT: Mandong-Cuobei thrust; SF: Shiquanhe fault; KF: Karakorum fault. Zhada: Zanda county; Shiquanhe: Shiquanhe city; Rutog: Rutog county; Domar: Domar country.

印度板块在青藏高原下的俯冲位置历来是地球科学研究的热点问题.远震体波研究显示在青藏高原下方，印度板块俯冲距离存在东西向差异(Zhao et al.，2010; Li et al.，2008).横波分裂研究可能为此提供了一种解释：印度岩石圈板片向北俯冲过程中，俯冲角度东陡西缓，导致俯冲距离东近西远，板片发生撕裂和断离(Chen et al.，2015).关于高原内部各块体边界处的地表断裂在地下空间如何展布也尚未达成共识，它们与地壳变形密切相关.为此，张忠杰等在青藏高原西部开展了人工源和天然源联合地震探测工作(TW-80试验，Zhang et al.，2014; Wu et al.，2015).综合已有的观测剖面，共组成4条地球物理探测走廊，如图 1a所示.TW-80试验接收函数研究揭示该剖面块体边界下方莫霍面存在错断，中法合作剖面接收函数和层析成像结果显示莫霍面在班公湖—怒江缝合带错断~10 km，羌塘块体下方加深至90 km，北至阿尔金断裂带突然抬升为50~60 km(Wittlinger et al.，2004).而INDEPTH深部地震探测结果显示雅鲁藏布江缝合带和班公湖—怒江缝合带下方没有莫霍层错动(赵文津等，2008；Kind et al.，2002)，表明青藏高原东西部地壳变形方式存在明显差异.为研究西部地壳变形方式，我们在接收函数研究的基础上进一步反演S波速度，分析壳内速度结构特征.

TW-80试验在青藏高原西缘80°E附近构建了一条南北向剖面，实施人工源宽角地震和天然源宽频带地震联合探测(Zhang et al.，2014).宽频带台站位置如图 1所示，共22台地震仪(Reftek-72A 数据采集器，周期范围50Hz-30/60s CMG3-ESP 地震仪)，台间距15~20 km，剖面全长400多公里，记录时间从2011年11月至2013年11月，采样频率40 Hz.图 1b中显示剖面横跨数个构造单元：甜水海块体、羌塘块体、拉萨块体及喜马拉雅块体.块体之间以缝合带或断裂带为界，构造线走向均呈北西向，都可以和高原东部的主要构造带相接.此地印度和欧亚板块碰撞十分强烈，发育大型走滑喀拉昆仑 断裂和龙木错断裂，数条逆冲断裂，及新生代拉分盆地.由北向南出露三条蛇绿岩带：班公湖—怒江蛇绿岩带、狮泉河蛇绿岩带和雅鲁藏布江蛇绿岩带，均为特提斯洋演化的产物(Lacassin et al.，2004).

甜水海块体和羌塘块体以龙木错断裂隔开，该断裂带SW走向，分布火山岩(李海兵等，2006).羌塘块体和拉萨块体以班公湖—怒江缝合带为界，研究区内是该缝合带西段.很多地球物理研究认为印 度板块岩石圈地幔已俯冲到该边界附近(Zhao et al.，2010; Kind and Yuan，2010;余大新等，2014).南端雅鲁藏布江缝合带是拉萨块体和喜马拉雅块体的分界线，为印度和欧亚两大板块于晚侏罗—白垩纪汇聚后，残留的新特提斯洋洋壳和岩石圈(Tapponnier et al.，1981).其北部是右旋走滑喀拉昆仑断裂带，图 1b中可看到印度河(印度境内)—狮泉河(中国境内)被右行错开约120 km(李海兵等，2008).两侧发育新生代盆地：札达盆地和噶尔盆地，其中札达盆地新生代地层厚约750 m，磁性地层学研究表明该盆地受喀拉昆仑断裂活动控制(王世锋等，2008).

利用接收函数研究地壳上地幔及地幔过渡带已成为一种重要的地球物理探测手段(Kind et al.，2002; Chen and Ai，2009;Tian et al.，2005a，b; Ai et al.，2007; Shi et al.，2009; 吴庆举和曾融生，1998；刘启元等，1996；李永华等，2006).本文数据处理包含两步：计算单个台站叠加接收函数，反演台站下方S波速度.我们采用时域迭代反褶积算法计算接收函数(Ligorría and Ammon，1999)，远震事件的震中距为30°~90°，回折点位于下地幔，可有效避开上地幔过渡带三重相或地核影区的影响(Çakır et al.，2000).为了提高反演结果的稳定性，计算接收函数时采用两种高斯滤波器——滤波系数分别为1.0和2.5，对应截止频率0.5 Hz、1.2 Hz.由于事件分布方位不均匀，为了提高信噪比，压制横向不均匀性，以台站为中心将方位角划分为不同小区域(每10°划为一个网格)，小区域内不同地震事件的接收函数进行叠加，剔除明显异常的结果，把剩下的不同区域叠加结果再相加，用最终的高信噪比叠加波形进行反演.图 2a是AL04台挑选出的接收函数，对应高斯系数1.0、2.5，最上面为叠加波形，可以看到信噪比很高，Moho面产生的转换震相PmS非常清晰.每个台站做类似处理，用两个叠加接收函数反演S波速度.图 2b为单个台站对应的两种叠加波形，其中AL08位于狮泉河断裂带上(图 1)，接收函数信噪比低，我们用邻近台站接收函数的平均值来代替.

图 2

Fig. 2

图 2 单个台站的接收函数 (a) AL04台挑选出的两种接收函数，最上面为叠加波形，左侧为各接收函数相应的射线参数； (b) 单个台站两种叠加接收函数，分别对应高斯系数2.5和1.0. Fig. 2 The selected receiver functions for each station (a) The receiver functions for the station AL04 filtered by the Gaussian filter 1.0 and 2.5, respectively. The top panels are the stacked results. The left panel gives the corresponding ray parameter for every single receiver function; (b) Two stacked results for each station, corresponding with the Gaussian filter 1.0 and 2.5, respectively.

接收函数反演是很强的非线性反演问题，目前主要有两种处理方法：一种是线性化反演方程，转化为线性反演问题；另一种是直接进行非线性反演.对于接收函数线性化反演，Ammon等(1990)的研究最具代表性，由于采用R and all算法计算微分地震图和跳跃反演技术，效率极高而被地球物理学家广泛采用(徐鸣洁等，2005；彭恒初等，2012；刘启民等，2014；Rai et al.，2006).Ammon等研究表明，线性反演强烈依赖初始模型.为克服此缺陷，很多全局优化算法被引入该领域，如模拟退火(Sen and Stoffa，1991)、遗传算法(Sambridge and Drijkoningen，1992)等. 此外还有一些新方法，刘启元等(1996)根据Tarantola(1987)的非线性反演理论，提出接收函数径向与垂向分量复谱比的非线性反演.Julià等(2000)、胡家富等(2005)进一步发展了接收函数和面波频散联合反演技术.

本文采用相邻算法(Neighbourhood algorithm，简称NA 算法)进行接收函数反演.NA算法是Sambridge(1999a，1999b)基于非线性反演理论提出的一种全局算法，具有较强的自适应搜索能力，用一组模型集代替搜索全局最优解，可避免反演解陷入目标函数的局部最小域，在地球物理反演中得到 广泛应用(Sherrington et al.，2004; Bannister et al.，2004; Hetényi et al.，2006; Snoke and Sambridge，2002; 贺传松等，2004).算法引入几何学概念——“维诺图”.在d维模型空间中给定一组随机采样点n_s，这些离散点将模型空间分割成n_s个区域，称为维诺单元(Voronoi cell)，每个维诺单元内的点到该采样点的距离都是最近的，用L₂准则测算.对每个维诺单元构建失配函数，从失配值最低的n_r个区域重采样，自适应实现高失配、低失配区域的稀少采样和密集采样.

模型空间内的地层包含如下参数：地层厚度、顶面S波速度、底面S波速度及层内P-S波速度比，层内速度为线性变化.图 3是两个模型实例.图 3a、3b对应模型1，模型参数及反演搜索空间见表 1.图 3b 中黑实线为合成的含噪声接收函数(采样频率25 Hz，记录时长30 s，共876个数据点)，相当于观测值，红实线为最佳反演结果合成的接收函数.可见，NA算法可以很好重建接收区的地壳速度结构.本研究区构造复杂，实际计算得到的接收函数除了Moho面的转换震相PmS外，还含有多个壳内间断面转换震相，致使PmS震相模糊或消失.因此，我们引入含沉 积层、结晶基底的模型2，与图 3c、3d相对应，模型参数及反演搜索空间见表 2.结果显示，此方法可以很好拟合壳内转换波震相.

图 3

Fig. 3

图 3 接收函数反演模型实例 (a)(b)对应模型1，(c)(d)对应模型2.(b)(d)中黑实线为合成的含噪声接收函数，红实线为最佳反演结果合成的理论值；(a)(c)中两条虚线指示模型搜索空间，黑实线为设定的真实地球模型，灰色区域为拟合度高的1000个反演模型集， 红实线为最佳结果，白实线为1000个模型集的平均值. Fig. 3 Two model examples of the NA algorithm to obtain the S-wave velocities (b)(d)The black line denotes the input receiver function, the red line denotes the synthetics based on the best data-fitting inversion result；(a)(c)The black dashed lines denote the parameter space, the gray lines denote the 1000 acceptable inversion results, the red line denotes the best data-fitting model, the white line denotes the average of the ensemble of 1000 models.

表 1(Table 1)

表 1 模型1(上)及模型空间参数(下) Table 1 Model 1 (Upper) and associated parameter space (Lower) 层 厚度(km) 顶面S波速度(km·s -1) 底面S波速度(km·s -1) V P/ V S 地壳 35 2.8 3.4 1.73 地幔 35 4.5 4.5 1.80 上地壳 1~10 1.7~3.2 1.7~3.2 1.6~1.8 中地壳 5~15 2.2~3.2 2.5~3.4 1.6~1.8 下地壳 5~25 2.6~3.5 3.4~3.6 1.6~1.9 地幔 20~40 3.8~4.5 4.0~5.0 1.6~1.9

表 1 模型1(上)及模型空间参数(下) Table 1 Model 1 (Upper) and associated parameter space (Lower)

表 2(Table 2)

表 2 模型2(上)及模型空间参数(下) Table 2 Model 2 (Upper) and associated parameter space (Lower) 层 厚度(km) 顶面S波速度(km·s -1) 底面S波速度(km·s -1) V P/ V S 沉积层 1 2.75 2.8 2.8 结晶基底 1 2.8 2.95 1.9 上地壳 2.5 2.95 3.7 1.7 中地壳 16 3.7 3.5 1.7 下地壳 17 3.6 4.2 1.7 地幔 20 4.41 4.85 1.7 沉积层 0~2 1.75~3.0 1.75~3.0 2.0~3.0 结晶基底 0~3 1.5~3.5 1.5~3.5 1.65~2.0 上地壳 1~15 2.6~3.6 2.8~4.0 1.65~1.8 中地壳 5~20 3.2~4.5 3.2~4.5 1.65~1.8 下地壳 5~20 3.2~4.5 3.2~4.5 1.65~1.8 地幔 5~30 4.0~5.0 4.0~5.0 1.7~1.9

表 2 模型2(上)及模型空间参数(下) Table 2 Model 2 (Upper) and associated parameter space (Lower)

利用上述NA算法反演单个台站下方S波速度.为增强反演结果可对比性，设定相同的模型搜索空间：沉积层、结晶基底、上地壳、中地壳和下地壳，模型参数见表 3.针对反演过程的参数设置为n_s=100，n_r=20，迭代次数为500，共获得50100个反演模型，取失配值最小的1000个模型计算平均值作为最终解.图 4为台站AL00、AL1011和AL21的反演结果.高斯系数2.5计算的接收函数含更多高频成分，除了来自莫霍面的转换震相，还包含多个壳内震相，而1.0计算的接收函数PmS信噪比很高，却损失了壳内探测精度，两者相互约束可以提供更丰富的速度信息.三个台站横跨喜马拉雅块体、拉萨块 体和羌塘块体，PmS-P波的到时差从~7.2 s(AL00)增至~9.5 s(AL21)，说明由南向北莫霍面逐渐加深，反演结果与之一致(图 7).图 5是所有台站接收函数与波形拟合结果.

表 3(Table 3)

表 3 反演设定的模型空间 Table 3 Parameters space used in the inversion 层 厚度(km) 顶面S波速度(km·s -1) 底面S波速度(km·s -1) V P/ V S 沉积层 0~5 1.5~3.0 1.5~3.0 2.0~3.0 结晶基底 5~20 1.7~3.2 1.7~3.2 1.65~2.0 上地壳 10~20 2.0~3.5 2.4~4.0 1.65~1.8 中地壳 10~20 2.5~4.0 2.8~4.5 1.65~1.8 下地壳 10~20 3.0~4.5 3.0~4.5 1.65~1.8 地幔 5~30 3.2~5.0 4.0~5.0 1.7~1.9

表 3 反演设定的模型空间 Table 3 Parameters space used in the inversion

图 4

Fig. 4

图 4 台站AL00、AL1011、AL21的接收函数反演结果 左图均为模型搜索空间及拟合度高的1000个模型集，线条含义参考图3.右图黑实线为实际观测值，红实线为拟合最好的反演模型合成的理论接收函数，上面波形对应高斯系数2.5，下面波形对应高斯系数1.0. Fig. 4 The inversion results of the AL00，AL1011，AL21 stations The left panel denotes the parameter space and the ensemble of good data-fitting 1000 models, the meaning of the colorful lines is similar with the Fig.3. The right panel denotes the observed receiver functions (the black line) and the synthetics (the red line) obtained from the best data-fitting inversion result.

图 5

Fig. 5

图 5 单个台站接收函数反演与波形拟合结果 左图黑线为单台采用高斯系数1.0计算的叠加接收函数，红线为拟合最佳的反演模型合成的理论值； 右图类似，对应高斯系数2.5. Fig. 5 The fitness between the synthetics obtained from the inversion result and the observed receiver functions for all stations The left panel denotes the stacked receiver functions with Gaussian filter 1.0 for each station, the black lines represent stacked result and the red lines represent the synthetics obtained from the best data-fitting inversion result. The right panel is similar, corresponding to the Gaussian filter 2.5.

图 6为本文反演的单个台站S波速度.图 7为剖面下方二维S波速度结构.TW-80试验用接收函数方法系统研究了该剖面下方的壳内结构和莫霍面形态，为便于结果对比，我们把接收函数叠加成像的结果投影在该速度图上.据观察，莫霍面与速度值4.0的等值线趋势大致吻合，研究区西侧Rai等(2006)反演结果类似.剖面下方莫霍面平均深~70 km，拉萨块体下方50 km 附近存在S波速度为4.0~4.1 km·s^-1的局部高速体.从地震活动(图 1b中紫色圆圈)的分布来看，地震主要集中于该剖面两端.南端喀拉昆仑断裂带下方，震源深度达53 km，北端龙木错断裂带下方，震源深度达60 km，表明这两条走滑断裂至少向下延伸到50~60 km.断裂活动在其周围形成的剪切破碎区域在S波速度图上表现为低速特征，对应剖面两端中上地壳的低速异常，而两条断裂带之间的中上地壳缺乏连续低速异常分布.

图 6

Fig. 6

图 6 单台S波速度反演结果 两条虚线指示模型搜索空间，灰色区域为拟合度高的1000个反演模型集，红实线为最佳反演结果，而白实线为1000个模型集的平均值. Fig. 6 The shear-wave velocity structure for each station The black dashed lines denote the parameter space, the gray lines denote the 1000 acceptable inversion results, the red line denotes the best data-fitting model, the white line denotes the average of the ensemble of 1000 models.

图 7

Fig. 7

图 7 测线下方二维S波速度结构 白色圆圈为区域79°E—81°E，31°N—35°N范围内自1900年至2015年发生的4级以上地震(地震目录来自USGS：http://earthquake.usgs.gov/earthquakes/search/)，壳内速度间断面和莫霍面结构(黑线)是该剖面接收函数结果(Zhang et al., 2014). Fig. 7 The shear-wave velocity map along the survey line The white circles denote the MS>4 earthquakes in the region 79°E—81°E，31°N—35°N, the crustal velocity discontinuity and the Moho structure plotted by black lines are derived from the previous receiver function study (Zhang et al., 2014).

图 8为研究区莫霍面结构及其横向变化，测线 位置见图 1a.西侧○ a 线综合 Rai等(2006)和Wittlinger等(2004)北段剖面的接收函数结果，○ b线由本文剖面、Caldwell等(2013)观测位置及Wittlinger等(2004)观测剖面的中段组成，○ c 、○ d线分别为Zhao等(2010)、Nábělek等(2009)接收函数结果.可见，莫霍面南北向变化趋势相似，从南端印度板块~40 km 向北逐渐加深，羌塘块体达~80 km，北至阿尔金断裂带下方又抬升至~60 km.整体来看，青藏高原西部主要块体边界切穿地壳，莫霍面发生错断(如喀拉昆仑断裂、班公湖—怒江缝合带和阿尔金断裂，图 8 中黑色虚线).图 7速度结构显示拉萨块体下方~50 km 处存在不连续高速间断面，该剖面接收函数结果在此深度范围有相应的转换波震相(Zhang et al.，2014)，综合南侧Caldwell等(2013)接收函数结果，推测此为MHT滑脱面向北延伸.Hi-CLIMB探测剖面(~85°E)观测到拉萨块体南部也存在一组壳内 转换波震相(Nábělek et al.，2009)，往东接收函数研究显示青藏高原南部下地壳存在速度间断面(Schulte-Pelkum et al.，2005)，结合本研究结果，推测该间断面下方为高速异常.综合分析认为，印度板块在青藏高原下的俯冲可能已到达班公湖—怒江缝合带附近，向北俯冲过程中印度下地壳发生榴辉岩化，根据藏南地区大地电磁结果推测由于不同部位含水量差异致使榴辉岩相变程度有所不同(金胜等，2007；魏文博等，2009).

图 8

Fig. 8

图 8 研究区Moho结构及其横向变化, 剖面位置参看图1a Fig. 8 The Moho structures based on the passive-source seismic survey lines shown in Fig.1a and the lateral variations in the study area

本文利用青藏高原西部TW-80试验宽频带流动台阵，在接收函数基础上反演S波速度结构.综合现有的青藏高原西部宽频带剖面探测结果认为：

(1)喜马拉雅块体和拉萨块体下方，S波速度结构显示下地壳存在10~20 km 厚的高速层，可能是印度下地壳向北俯冲过程中发生榴辉岩化所致；

(2)地震活动集中在剖面两端：喀拉昆仑断裂带和龙木错断裂带附近.南端喀拉昆仑断裂带下方震源深度达53 km，北端龙木错断裂带下方震源深度达60 km，表明两者均为深断裂，断裂活动造成的剪切破碎域致使剖面两端中上地壳表现为低速异常；

(3)与青藏高原东北缘上地壳增厚模式(Tian and Zhang，2013; Tian et al.，2014；Xu et al.，2014; Zhang et al.，2011)、东南缘下地壳流增厚模式(Royden et al.，1997; Clark et al.，2005)不同，喀拉昆仑断裂带和龙木错断裂带之间的中上地壳没有观测到连续的低速异常，说明可能缺乏解耦层，支持青藏高原西部地壳以整体缩短增厚模式为主，正在开展的人工源宽角地震资料和背景噪声成像研究有望进一步提供约束.