引言

50万年前印度板块与欧亚板块相互碰撞，随后印度板块向北持续推挤，而后遭受到阿拉善和鄂尔多斯等坚硬地块的阻挡，导致青藏高原的地壳缩短增厚以及强烈的地表形变，青藏高原地区因此成为板块碰撞、物质运移和地球动力学研究的理想区域。

青藏高原东北缘是青藏高原、鄂尔多斯和阿拉善等3个地块的结合部位，也是青藏高原板块向大陆内部扩展的前缘部位，该区域构造活动强烈，横向不均匀性很强，存在多个活动断裂。余大新等（2014）根据地震面波层析成像结果推断亚洲板块俯冲的前缘位置已抵达柴达木盆地。该区域还经历过强烈的地壳缩短增厚运动，地壳厚度变化明显，从西部的60 km逐渐减薄至东部的40 km （赵金仁等，2005；刘启民等，2014；Li et al，2015 ）。关于该区域的地壳增厚机制一直以来都受到学术界的关注，但尚处于争议阶段，目前主要有两种模型：青藏高原多阶段隆升模型（Tapponnier et al，2001 ）、中下地壳流模型（Clark，Royden，2000）。已有研究显示青藏高原东北缘存在中下地壳低速体（陈九辉等，2005；李永华等，2006；嘉世旭，张先康，2008；Li et al，2014 ）；Jiang等（2014）认为中下地壳流是青藏高原东北缘的主要生长机制，但是Pan和Niu （2011）指出青藏高原东北缘的物质组成与青藏高原中心存在差异，中下地壳流不能很好地解释其机制；刘启民等（2014）也认为中下地壳流在青藏高原东北缘可能并不存在；Li等（2014）的研究显示祁连造山带的地壳低速层与华北克拉通和青藏高原之间的地壳缩短有关；Li等（2015）推测青藏高原西部的地壳增厚是由于欧亚板块与印度板块的重叠作用引起，而东部地壳的减薄由热物质上涌导致；Zheng等（2016）的研究表明青藏高原东北缘的生长机制不符合中下地壳流模型。因此对于青藏高原东北缘地区来说，开展精细的区域速度结构研究，厘清物质运移和地壳增厚机制具有重大意义。

本文将基于研究区域内固定地震台站的波形数据，使用噪声层析成像方法给出青藏高原东北缘不同周期的群速度分布图像和三维剪切波速度结构，以期为地壳增厚机制和区域地球动力学研究提供参考约束条件。

1.   数据与方法

本文收集了青藏高原东北缘的青海、甘肃、宁夏等3个省级地震台网2009年9月至2011年9月期间87个宽频带地震台站所记录到的垂直向波形数据，所使用地震台站的分布见图1。

本文按照Bensen等（2007）的方法进行数据处理，主要流程为：① 原始垂直向地震波形记录采样率为100 Hz，使用低通滤波器滤波后重采样为1 Hz，以北京时间为准，单个台站一天的数据生成一个文件，并按照规定的命名格式重新命名使得后续数据的处理实现完全自动化。这一步需要考虑数据间断的影响，若当天数据间断超过10%则被剔除；② 使用SAC软件对波形数据去平均值，再扣除仪器响应后得到速度振幅；③ 使用绝对振幅滑动平均法剔除地震和其它强振幅所带来的影响，滑动窗长度为101个采样点；④ 进行5—50 s周期的带通滤波和频率域谱白化，获得“白噪声”；⑤ 按天计算台站对的互相关函数并进行线性叠加。图2给出了玉树地震台（YUS）与其它地震台站之间的经验格林函数，并按照台站对距离排列，可以看出整个包络的平均速度大约为3 km/s，但有些波形未显示出完全对称，这是由于噪声源分布随时间变化所致，因此本文对正负互相关分支进行叠加以提高信噪比；同时使用两年的台站对波形进行互相关叠加，以减小噪声源空间分布随季节变化的影响。在数据质量控制方面，采用经验格林函数的信噪比。图3给出了8，10，15，20，30和40 s的信噪比统计图，可以看出周期为10，15和20 s的信号信噪比最好，且随着周期增加信噪比衰减。本文仅使用信噪比高于10的互相关波形。

图  1  青藏高原东北缘地震台站分布及区域构造简图

F₁：青铜峡—固原断裂；F₂：六盘山断裂；F₃：海原断裂；F₄：门源断裂；F₅：西秦岭北缘断裂；F₆：东昆仑断裂；F₇：主峰断裂；F₈：龙门山断裂带

Figure  1.  Distribution of seismic stations and regional tectonic settings in northeastern Tibetan Plateau

F₁：Qingtongxia-Guyuan fault；F₂：Liupanshan fault；F₃：Haiyuan fault；F₄：Menyuan fault；F₅：North edge of west Qinling fault；F₆： Eastern Kunlun fault；F₇：Zhufeng fault；F₈：Longmenshan fault zone

下载: 全尺寸图片 幻灯片

图  2  玉树台与其它台站噪声互相关结果展示图

Figure  2.  Demonstration of cross correlations of long time waveforms between the station YUS and other ones

下载: 全尺寸图片 幻灯片

图  3  不同周期T上台站对的互相关信噪比统计图

Figure  3.  Histogram of cross correlation signal-noise-ratio between station pairs at different periods T

下载: 全尺寸图片 幻灯片

获得台站对经验格林函数后测量频散曲线，关于频散曲线Bensen等（2007）已经进行了较详细的描述，这里大致介绍方法。频散曲线测量使用传统的时频分析法，鉴于台站对的数量较多，实际处理使用计算机程序自动完成。通过希尔伯特（Hilbert）变换构建解析函数，然后使用傅里叶变换并且应用多个高斯滤波器得到二维的包络函数以及对应的相位函数，再分别求取群速度和相速度。自动测量频散曲线需要克服两方面的影响：一方面来自于频谱缺陷、多重反射波和散射波的影响，由噪声互相关得到的经验格林函数受其影响少于天然地震；另一方面来自于前驱波的影响，可通过相位匹配滤波器减少其影响。测定频散曲线前需进行质量控制，满足信噪比大于10且台站间隔大于2倍波长的条件。图4给出了青海省德令哈台与花土沟台的频散曲线测量实例。

图  4  青海德令哈台与青海花土沟台经验格林函数（上）及其时频分析图像（下）

Figure  4.  Empirical Green’s function （top） between the stations DLH and HTG and their frequency-time analysis diagram （bottom）

下载: 全尺寸图片 幻灯片

图  5  8—40 s周期群速度成像所使用的射线路径分布

Figure  5.  Distribution of ray paths used in the group velocity tomography for periods from 8 s to 40 s

下载: 全尺寸图片 幻灯片

图  6  周期为8，10，15，20，30 和40 s的瑞雷面波群成像检测板结果

Figure  6.  Checkerboard test of Rayleigh wave group velocity tomography at periods 8，10，15，20，30 and 40 s

下载: 全尺寸图片 幻灯片

2.   群速度层析成像

本文利用FMST程序（Rawlinson，Sambridge，2004a，b，2005）来开展各周期面波群速度层析成像。研究区域被划分成0.5°×0.5°的空间网格，FMST程序使用快速行进算法（fast marching method，简写为FMM）计算走时，计算效率较高。FMM方法不需要假设射线路径沿大圆路径传播，而且考虑了射线路径在复杂介质中传播的弯曲情况。Fang等（2015）使用FMM正演，提出了基于小波变换的稀疏约束层析成像方法，该方法适用于介质不均匀性强的区域。本研究区域的地壳横向不均匀性较强，上地壳沉积层厚度变化较大，中下地壳存在局部低速层，地壳深度分布横向变化显著，这些特点会使波动传播偏离大圆路径，因此宜使用FMM方法。

图5给出了周期为8，10，15，20，30和40 s的面波群速度二维反演所用的射线路径，基于这些射线路径本文开展了相应周期的检测板试验。速度模型的初始速度设为3 km/s，添加正负交替分布且异常幅度为10%的速度异常体，其大小为2°×2°（图6左上）。测试结果（图6）显示，各周期面波反演得到的速度结构在台站分布密集处均能得到较好的反演结果，其中10，15和20 s周期面波得到的速度结构更可靠，这是由于10—20 s周期范围内可用的射线更多。对于研究区域边缘射线路径覆盖不到的区域（尤其是研究区域西部），反演结果的可靠度偏低。

图7给出了周期为8，10，15，20，30和40 s的群速度二维分布。由图7a可以看出，当周期为8 s时，群速度较低区域对应于较厚的沉积层，高速体、低速体的分布呈明显的条带状特征。研究区域内祁连山以南的柴达木盆地、北部河西走廊和东部的鄂尔多斯地块西缘均表现为低速特征，而高速体范围较小，祁连山山脉、西秦岭造山带和青藏高原结合处表现为高速体。在10 s周期群速度图像（图7b）上，高速体向北部河西走廊以及南部的柴达木盆地以东扩展，范围显著扩大；同时高速异常向鄂尔多斯地块扩展，但鄂尔多斯西缘仍然呈现为低速异常。整体而言，低速体与沉积盆地对应，高速体与造山带对应，与Yang等（2007）给出的8 s和10 s周期的群速度与地壳沉积层有密切关系的结论相一致。

由图7c可以看到：高速体在周期为15 s的群速度分布图上继续扩展，祁连山以北的河西走廊西部已表现为高速异常，且继续向东靠近鄂尔多斯地块；河西走廊东部持续表现为低速。周期为20 s的群速度分布（图7d）表现出其独特性：群速度整体呈低速异常，这与该区域中地壳普遍存在低速层的研究结果（Bao et al，2013 ；Li et al，2014 ）相对应；高速异常分布面积大大减小，但祁连山造山带可见不连续分布的片状高速带，阿拉善地块、青藏地块和鄂尔多斯地块结合处均呈高速异常。

周期为30 s的群速度与莫霍面深度有关；鉴于该区域莫霍面较深且空间差异大（图8，根据杨志高和张雪梅（2016）绘制），40 s周期面波群速度对于研究莫霍面埋深有参考价值。由图7e可见，周期为30 s的群速度分布呈明显特征：青藏高原东北缘广大地区表现为低速异常，向北延伸至河西走廊大部分地区，向东延伸至青藏东缘；鄂尔多斯地块西缘及其以西区域表现为高速异常，向北延伸至阿拉善地块；阿拉善地块南缘与河西走廊结合处也表现为高速异常。这说明青藏地块的莫霍面较深，且其深度向北、向东逐渐减小，尤其是青藏高原以东莫霍面深度减小的现象更明显。因为该区域莫霍面较深，40 s周期的面波群速度能更好地反映莫霍面深度信息。在40 s周期面波群速度图像（图7f）上，高低速异常的分界面更加明显，较30 s周期群速度图像所示的高速异常进一步发展，表现出更高的速度值。与30 s周期群速度特征不同的是，40 s周期群速度在祁连山东部向南出现了明显的高速体异常，大体上位于秦岭造山带与东昆仑造山带之间的过渡区。40 s周期群速度相比30 s周期群速度，其高异常值范围和程度继续加深，显示出莫霍面厚度自青藏高原向北和向东呈现减薄特征。

3.   剪切波速度结构

由8—40 s周期二维速度结构可以得到每个0.5°×0.5°网格点下方8—12 s周期的频散曲线，利用这些频散曲线可以反演每个网格点下方的一维剪切波速度（Herrmann，Ammon，2002）。初始速度模型选取Crust1.0 （Laske et al，2015 ）速度模型对应的每个网格点下方的层状速度模型。将每个网格点反演的一维速度结构合并后，得到研究区域以0.5°×0.5°网格划分的三维速度结构。Crust1.0模型包括水层、冰层、沉积层、上地壳、中地壳和下地壳层，由于本文的研究区域处于内陆地区，排除了水层和冰层。初始速度模型从地表到60 km深度范围采用厚度为5 km均匀分层结构，层速度值等于Crust1.0模型中对应深度范围的S波速平均值，地壳厚度采用接收函数的计算结果（杨志高，张雪梅，2016）。

图9给出了5—55 km深度范围内的剪切波三维速度分布，整体上来看研究区域的速度结构表现出较强的横向不均匀性。为确保该结果的可靠性，图9也给出了射线路径覆盖程度好的网格点（每个网格点内至少有50条射线穿过），而且所给出的网格点分布与检测板试验分辨率高的区域相对应。可以看出：处于5—10 km深度的低速体主要集中在沉积层，例如柴达木盆地和祁连山地块北部的河西走廊表现为低速特征；处于10—15 km深度的柴达木盆地继续保持低速特征，这与柴达木盆地较厚的沉积层相对应；处于20—25 km深度的速度结构上，甘孜地块大体表现为低速异常，柴达木盆地的低速体规模减小，柴达木地块东部呈高速特征；处于30—35 km深度的甘孜地块和祁连山地块均表现为低速异常，柴达木盆地出现高速异常体；处于40—45 km深度的祁连地块北东部保持低速异常，南西部呈高速异常，柴达木盆地呈高速异常；50—55 km深度上，研究区域的东部以高速体为主，西南部以低速体为主。

4.   结果分析

以往的研究显示青藏高原东北缘存在低速体异常，低速体分布在羌塘地块、松潘地块和祁连山造山带（Bao et al，2013 ；Li et al，2013 ；Li et al，2014 ；Jiang et al，2014 ；Zheng et al，2016 ）。Li等（2014）认为祁连造山带的低速体速度值较低，是由于华北克拉通和青藏地块的地壳缩短增厚所致；但Bao等（2013）认为祁连地块的低速体有相对高的速度值，推测祁连造山带处于高原生长的初级阶段。本文的层析成像结果（图9）显示甘孜地块和祁连地块的低速体在25—45 km深度范围内分布广泛，甘孜地块的低速体向北止于东昆仑断裂带以北，祁连地块北西低速体的规模大于其南东部分，这与前人的研究结果（Bao et al，2013 ; Li et al，2013 ）基本一致。柴达木地块10—35 km深度大体表现出高速特征，并将祁连地块的低速体“孤立”，但其40—45 km深度所表现出的低速异常将祁连地块与甘孜地块的低速体相连接，使得祁连造山带的低速体不再孤立，这与Bao等（2013）的结果有一定差别。Li等（2014）认为祁连造山带的地壳低速层是由于华北克拉通和青藏高原地块地壳缩短产生的热导致部分地壳熔融。为了更好地探查祁连山地块深部的三维速度结构，图10给出了速度结构的3个深度剖面。图10a和b展示出祁连山地块下方存在一个低速体，但从图10c上可以看出祁连造山带呈高速特征。王琼等（2013）使用远震波形数据获得青藏高原东北缘台站下方的各向异性，其结果显示祁连山和河西构造区域的上地幔各向异性特征与地壳存在一定差别，表明二者的壳幔变形机制可能有所差别；而隆中盆地及其周缘的上地幔各向异性与地壳各向异性特征比较一致，说明其壳幔变形机制可能相同。郭桂红等（2015）关于地壳各向异性的研究显示祁连山—河西活动构造区和甘东南活动构造区的变形不同，展示出祁连山北西部到鄂尔多斯西部地壳和地幔各向异性的差异，这与本文结果祁连山北西部下方低速体分布相对应。结合图10和各向异性结果推测，祁连山下方的低速体可能是中朝板块与青藏地块相互作用而产生的局部下地壳下沉所致。

Zhu和Helmberger （1998）通过远震体波研究了柴达木盆地和青藏高原的地壳厚度，其结果显示柴达木盆地与青藏高原交界处5 km条带范围内存在15—20 km的地壳厚度陡变，这可能意味着青藏高原受到南部印度板块的侵入，被北部坚硬的柴达木地块所阻挡而使地壳缩短增厚。Zheng等（2016）认为甘孜地块和柴达木地块的波速比显著高于祁连地块，意味着该区域地下岩石以铁镁质矿物为主或存在部分熔融现象。李永华等（2006）的研究显示，祁连地块的波速比随着地壳厚度的增加而降低，这意味着该区域的地壳增厚主要发生在上地壳，且岩石矿物以长英质为主。从本文所得速度结构来看，祁连地块和甘孜地块的低速体由柴达木地块东侧的深部物质相连，并非完全“孤立”，因此柴达木地块东部似乎可能成为青藏高原物质北移的一个通道。祁连地块的物质组成与柴达木地块和甘孜地块有明显不同，这可能是因为祁连地块尚处于青藏高原生长的早期，青藏高原物质的北向逃逸程度较弱。另一方面，该区域的地壳增厚主要是由于区域应力作用下发生在上地壳的横向缩短变形所致，而非青藏高原物质的北向逃逸所引起。柴达木地块和祁连地块的横向不均匀性高于甘孜地块，甘孜地块相对较弱的横向不均匀性可能意味着该区域处在高原生长的“稳定”阶段，柴达木地块和祁连地块更复杂的横向不均匀性则代表生长前缘更复杂的动力学过程。Bao等（2013）给出的速度模型显示祁连地块的低速体速度高于甘孜地块，并且祁连山的高程低于青藏高原，因此认为祁连地块处于青藏高原生长的初级阶段，是该高原生长的前缘。从GPS结果（Zhang et al，2004）来看，本文研究区域的速度场方向在95°E—100°E之间由大体北东向快速转变成东向，意味着柴达木地块存在物质北向运移的可能性。柴达木地块东部深处的低速体在青藏高原的生长过程中所起的作用是否为青藏高原物质北向入侵提供通道，尚需更高分辨率的层析成像结果参与。

从本文的群速度图像上来看，在研究区域东北部的阿拉善地块和鄂尔多斯地块西缘并未发现明显的低速体分布，但该区域的横向不均匀性较强。阿拉善地块和鄂尔多斯地块被认为是古老的稳定地块，从速度结构上来看该区域的波速显著高于研究区域的西南部，验证了该区域构造的稳定性。但是，从40—55 km深度的速度结构图上可以看出明显的横向不均匀性，根据Pan和Niu（2011）的接收函数结果，阿拉善地块的平均地壳厚度为42 km，而本文结果意味着该区域上地幔速度结构的不均匀性较强，因此有必要开展更深入的研究。

图  7  周期为8 s （a），10 s （b），15 s （c），20 s （d），30 s （e）和40 s （f）的瑞雷面波群速度图

Figure  7.  Rayleigh wave group velocity map at periods 8 s （a），10 s （b），15 s （c），20 s （d），30 s （e） and 40 s （f）

下载: 全尺寸图片 幻灯片

图  8  青藏东北缘地区的地壳厚度分布图

黑点表示计算接收函数所用台站

Figure  8.  Distribution of crustal thickness in the northeastern Tibetan Plateau

Black dots are seismic stations used for calculating receiver functions

下载: 全尺寸图片 幻灯片

图  9  5—55 km深度范围内的剪切波速度结构

Figure  9.  Shear wave velocity in the depth range of 5—55 km

（a） 5—10 km；（b） 10—15 km；（c） 20—25 km；（d） 30—35 km；（e） 40—45 km；（f） 50—55 km

下载: 全尺寸图片 幻灯片

图  10  沿图9f中剖线AA′，BB′和CC′的剪切波速度结构剖面图

Figure  10.  Vertical cross sections of shear wave velocity along the lines AA′，BB′ and CC′ in Fig. 9f

下载: 全尺寸图片 幻灯片

5.   讨论与结论

FMST程序计算的理论走时考虑到射线路径偏离大圆弧的效应，相比传统基于射线传播的层析成像方法，FMST算法更适用于横向不均匀性强的区域。本文的速度结构与先前的研究大体保持一致，浅层低速体与沉积层有较好的对应关系，浅层高速结构与造山带也有相关性。研究区域的中地壳广泛分布有低速体，包括甘孜地块、东昆仑断裂带以北部分区域和祁连造山带。祁连造山带的低速体大体上是孤立的，它与甘孜地块的高速体被柴达木地块分割开来，但柴达木地块东部40—45 km深处存在低速体，该低速层将甘孜地块与祁连地块的低速体相连。祁连地块与甘孜地块的物质组成和地壳增厚过程可能有明显差别：甘孜地块代表青藏高原的相对稳定和成熟地块，受青藏高原物质逃逸的影响较大；而祁连地块处于青藏高原生长的早期，地壳增厚可能主要源于中上地壳的缩短增厚；祁连造山带可能存在下地壳局部下沉现象。研究区域的西南部台站密度较低，下一步将引入流动观测数据，提高该区域的成像分辨能力。地下介质的变形可能在介质的各向异性上有所体现，后续研究将引入勒夫波噪声层析成像结果，联合反演地壳介质的径向各向异性，以进一步讨论区域地壳介质形变特征。

审稿专家提出了建设性的修改意见，M. H. Ritzwoller教授及其研究小组成员在测量面波频散曲线上给予了帮助，Nick Rawlinson教授提供了面波层析成像程序FMST，本研究成果图件使用GMT软件（Wessel et al，2013 ）绘制，作者在此表示衷心的感谢！