基于背景噪声研究青藏高原东缘三维横波速度结构

谭友恒， 于湘伟， 宋倩， 章文波

中国科学院大学地球与行星科学学院， 北京 100049

0 引言

青藏高原东缘是青藏高原向东扩展的前缘部位，也是青藏高原现今地表抬升和地壳增厚最强烈的地区之一.该区域地质构造复杂，发育有规模宏大的造山带和一系列大规模的深大断裂，如六盘山断裂(F1)、西秦岭断裂带(F2)、东昆仑断裂带(F3)、塔藏断裂(F4)、龙门山断裂带(F5)和鲜水河断裂带(F6)等(图1a).在复杂构造变形的作用下，青藏高原东缘强震频发，自2008年5月12日汶川MS8.0地震以来，又先后发生两次七级地震，即2013年4月20日芦山MS7.0地震和2017年8月8日九寨沟MS7.0地震.青藏高原东缘主要由祁连地块、柴达木地块、巴颜喀拉地块、川滇地块北部、鄂尔多斯地块以及部分华南地块等构造单元组成(图1).因此，青藏高原东缘是探索青藏高原隆升过程、变形机制和强震孕震背景的关键地区，其深部速度结构对于进一步了解青藏高原的构造演化具有重要的科学意义.

近年来，基于背景噪声的面波成像法在地球物理研究领域得到了广泛应用.理论研究表明，在噪声源分布均匀的扩散场中，通过对两个地震台站长时间的噪声记录进行互相关计算，可以提取台站对之间的经验格林函数，进而获取短周期面波频散信息并应用于反演地球内部三维横波速度结构(Campillo and Paul, 2003; Roux et al., 2005; Weaver, 2005; Shapiro et al., 2005; Yao et al., 2008).与地震面波成像法相比，背景噪声面波成像法具有两点独特的优势.第一，不依赖于天然地震或人工震源激发信号，是一种无源成像；第二，该方法可以获得短周期(小于30 s)的面波频散信息，从而提高了地壳浅层的分辨能力，进一步拓宽了数据的频带范围.

目前，围绕青藏高原东缘速度结构利用不同研究方法和数据获得了一系列研究结果.如：背景噪声成像法(Yao et al., 2008; Li H Y et al., 2009; Zheng et al., 2010; Yang et al., 2012; Xie et al., 2013; Li X F et al., 2014; Jiang et al., 2014;范文渊等，2015；杨志高等，2019；Bao et al., 2020)、聚束分析成像法(Wang et al., 2020)、地震体波走时成像法(Riaz et al., 2019；夏思茹等，2021)以及接收函数与瑞利面波频散联合反演法(Liu et al., 2014; Bao et al., 2015; Zheng et al., 2016)等.这些研究结果均表明，青藏高原东缘中、下地壳内部存在明显的低速体和正径向各向异性的特征，与下地壳流动模型(Royden et al., 1997; Clark and Royden, 2000; Royden et al., 2008)更为符合.此外，人工深地震探测方法(Wang et al., 2007)的研究结果表明青藏高原下地壳有较高的衰减值，可能是下地壳低速流动的结果.尽管这些研究结果在青藏高原东缘存在壳内低速体方面达成共识，但对低速体空间分布特征和成因方面有不同看法.联合面波频散曲线和接收函数反演方法的结果(Zheng et al., 2016)显示，青藏高原中、下地壳低速体向东北方向扩展但并未跨越东昆仑断裂带(F3)，表明下地壳流动模型并非青藏高原东北部发育的主要机制，范文渊等(2015)的背景噪声研究结果也得到了类似的结论.但也有部分研究(Zheng et al., 2010; Yang et al., 2012; Xie et al., 2013; Jiang et al., 2014; Deng et al., 2015)表明，青藏高原中、下地壳低速体在向北部扩展过程中明显跨过了东昆仑断裂带(F3)，其宽度向高原东缘逐渐缩小，但未到达西秦岭断裂带(F2).此外，Bao等(2020)背景噪声成像法得到了与以往研究不同的结果，认为巴颜喀拉地块地下并非存在大规模的低速体，而是沿多个孤立带分布，与青藏高原东缘大尺度下地壳流动模型不一致.这些均表明青藏高原东缘低速体的成因、空间分布等特征仍存在较大争议，还需要更深入地研究和分析.因此，青藏高原东缘，特别是巴颜喀拉地块附近的地壳和上地幔精细的三维结构以及变形机制仍是一个需要继续深入研究的问题.

本文利用四川、重庆、甘肃、青海、陕西地震台网的105个宽频带地震台站在2017年8月—2018年8月期间记录的连续波形数据，采用背景噪声成像法反演了青藏高原东缘地壳和上地幔顶部高分辨率的三维横波速度结构，探究青藏高原东缘中、下地壳低速体的成因和空间的分布特征，为深入了解该区域地下介质信息和中、下地壳低速体分布提供重要参考依据.

1 数据与方法

本文收集了研究区域(99°E—109°E、28°N—38°N)内105个宽频带固定地震台(图1a)在2017年8月—2018年8月期间记录的连续波形数据，对其垂直分量的背景噪声进行互相关处理并提取基阶瑞利面波相速度频散曲线.基于Bensen等(2007)提出的关于背景噪声成像法的处理流程，对提取的背景噪声数据进行单台数据预处理、互相关数据的提取和叠加、瑞利面波频散曲线的测量及其严格控制和筛选，从而获得基阶瑞利面波相速度频散曲线.

图1 研究区地块、主要活动断裂和地震台站分布图(a)图为(b)图中红色实线框区域的放大图，图中断裂数据来自邓起东等(2002)，白色实心圆为2000年以来发生的MS≥6.0地震事件，震源参数来源于中国地震台网中心正式地震目录(https:∥data.earthquake.cn/gcywfl/index.html).主要活动断裂：F1：六盘山断裂，F2：西秦岭断裂带，F3：东昆仑断裂带，F4：塔藏断裂，F5：龙门山断裂带，F6：鲜水河断裂带.(b)图中地块划分：TB：塔里木盆地，ALB：阿拉善地块，ODB：鄂尔多斯地块，QB：柴达木地块，QLB：祁连地块，BYB：巴颜喀拉地块，WQO：西秦岭造山带，QTB：羌塘地块，CDB：川滇地块，SCB：四川盆地，LSB：拉萨地块，SChB：华南地块.Fig.1 The distribution of blocks, the major active faults and seismic stations in the study area (a) is an enlargement of the red solid boxed area in (b), and the faults in the figure are from Deng et al. (2002). The white solid circles are the seismic events with MS≥6.0 since 2000, and the focal parameters are from the official earthquake catalogue of China Earthquakes Networks Center (https:∥data.earthquake.cn/gcywfl/index.html). The major active faults: F1: Liupanshan fault, F2: West Qinling fault zone, F3: East Kunlun fault zone, F4: Tazang fault, F5: Longmenshan fault zone, F6: Xianshuihe fault zone. The block in (b): TB: Tarim Basin, ALB: Alxa Block, ODB: Ordos Block, QB: Qaidam Block, QLB: Qilian Block, BYB: Bayan Har Block, WQO: West Qinling Orogen, QTB: Qiangtang Block, CDB: Sichuan-Yunnan Block, SCB: Sichuan Basin, LSB: Lhasa Block, SChB: South China Block.

关于背景噪声成像的理论和方法很多文献已做过详细介绍(Weaver and Lobkis, 2002; Shapiro and Campillo, 2004; Weaver and Lobkis, 2004; Roux et al., 2005; Shapiro et al., 2005; Yao et al., 2006; Bensen et al., 2007; Yang et al., 2008; Li H Y et al., 2009; Fang et al., 2010; Bao et al., 2013).传统的面波成像法主要分为两步：(1)利用路径频散数据反演获得面波相速度或群速度二维分布；(2)利用每个网格点下方的频散数据反演获得三维横波速度.本文采用面波频散数据直接反演法(Fang et al., 2015)，可以直接从不同周期不同路径的面波频散走时数据反演青藏高原东缘地壳和上地幔顶部三维横波速度结构.与传统方法相比，该方法具有以下优点：(1)不需要构造二维相速度或群速度图；(2)对结构的横向变化有较好的平滑约束；(3)与全波形反演或伴随层析成像相比，其计算效率更高.

1.1 单台数据预处理

本文将收集的105个固定台以天为单位、频率为100 Hz的原始连续波形数据进行降采样(5 Hz)处理后并依次进行去均值、去线性趋势、去除仪器响应、带通滤波(0.018～0.625 Hz)、谱白化等数据处理，采用滑动绝对值平均归一化的方法对数据在时间域内进行归一化，目的是减少地震信息、仪器不统一和台站附近非平稳噪声源的影响(Bensen et al., 2007).图2a为四川盆地内台站HYS(图3a)记录到的波形原始频谱图，低频区域存在的几个波峰可能是由于微震信号造成的，且零频处存在明显的异常值，这些信号对后续频散曲线的提取均会产生较大干扰.经过上述预处理后，零频处的异常值和低频部分的波峰被有效去除(图2b).因此，通过上述处理方式可以有效地降低微震信号的影响，去除异常信息，拓宽背景噪声信号的频带.

图2 原始波形(a)与谱白化(b)后的频谱图Fig.2 The frequency spectrum diagram of the original waveform (a) and the spectrum whitening (b)

1.2 互相关函数的提取

背景噪声成像法中的一个关键步骤是对预处理后的波形数据做互相关计算，即通过对两个台站记录到的长时间噪声资料进行互相关处理以提取台站间介质的格林函数(Campillo and Paul, 2003).对上述预处理后的任一台站对的垂直噪声数据进行以天为单位的互相关计算，通过13个月的线性叠加从而获得最终的互相关函数.图3为台站HYS与其他台站之间路径分布和13个月叠加后的互相关函数，不同方位角的互相关函数中呈现清晰的信号(图3b中虚线对应的范围).理论上，若噪声源分布均匀，互相关函数正负半轴信号应具有对称性，但实际测量中，由于单条互相关函数受到噪声源分布不均匀的影响，其正负两支往往会出现不对称的现象(Stehly et al., 2006)(图3b).因此在实际数据处理中，通常把正负两支反序叠加后再形成“对称分量”互相关函数(Yang et al., 2007).最后，将“对称分量”互相关函数做希尔伯特转换，以此获得面波经验格林函数(Yao et al., 2011)，作为接下来提取频散曲线的数据基础.

1.3 瑞利面波相速度测量

本文采用基于图像分析的面波相速度频散快速提取的方法(姚华建等, 2004; Yao et al., 2005)，并根据面波格林函数的远场表示定理来自动提取整条频散曲线(Yao et al., 2006, 2008).具体分为如下四个步骤：

(1)用多重滤波的方法计算基阶瑞利面波群走时，并以此设计余弦盾的窗函数，即每个周期窗函数的中心时刻对应该周期的群走时，由此获得加窗后的经验格林函数，这一步可有效地消除高阶面波和其他干扰信号.

(2)以周期带宽为0.4 s的凯泽窗(Kaiser window)窄带滤波器对加窗后的经验格林函数进行滤波，获得周期5～50 s内近似单周期的瑞利面波频散信号.

(3)按周期大小的顺序将单周期信号从小到大排列，获得相走时-周期的振幅矩阵，并根据基阶面波格林函数在远场的近似展开，得到相走时-相速度的转换公式(1)，从而计算相速度-周期的振幅矩阵，

(1)

其中T表示测量周期，cAB(T)代表周期T对应的相速度，t为相走时，Δ是台间距.使用三次样条插值的方法得到了振幅矩阵的图像(图4)，振幅矩阵图单列中不同颜色对应做窄带滤波后单周期的经验格林函数的归一化振幅大小.

(4)在经验速度范围窗内搜索振幅矩阵图像中振幅的最大值，获得基阶瑞利面波相速度频散曲线(图4，白色实线).

图4 台站对CHK-REG的相速度频散曲线测量示意图白色实线表示初次提取的相速度频散曲线.蓝色实心圆表示经过“相邻周期无明显跳变”的原则筛选后的频散点，黑色实心圆表示按照信噪比的原则筛选后的频散点，红色实心圆表示经过面波格林函数的远场表示定理计算出来的频散点.台站对CHK-REG的路径分布见图5a中黄色实线.Fig.4 Phase velocity dispersion curve measurement of the station pair CHK-REG The white solid line represents the phase velocity dispersion curve extracted for the first time. The blue solid circles represent the dispersion points sifted by the principle of “no significant jump in adjacent periods”, the black solid circles represent the dispersion points sifted by the principle of signal-to-noise ratio, and the red solid circles represent the dispersion points calculated by the far-field representation theorem of the surface-wave Green′s function. The path distribution of station pair CHK-REG is shown as the yellow solid line in Fig.5a.

为了增加频散数据的可靠性，本文在筛选频散曲线时采用如下三个标准：①单频信噪比≥5；②台间距>2倍波长；③相邻周期相速度值无明显跳变，最终获得5416条频散曲线.图5b展示了研究区域内台站对筛选后的频散曲线，值得注意的是川滇地块附近分布的台站对BTA-MBI、巴颜喀拉地块附近分布的台站对DAW-GZI，以及巴颜喀拉地块和柴达木地块附近分布的台站对LJS-RTA所对应的频散曲线在周期约为7～15 s范围内出现负斜率，一定程度上说明这些地块下方可能存在低速体.然而，四川盆地和鄂尔多斯地块附近分布的台站对HZHG-YGD、HCH-QCH和HXT-WXT所对应的频散曲线随着周期的增加呈现单调递增的特征，说明不同地块下方速度结构存在显著差异.本文将频散数据点按各周期的路径数进行统计(图6)，显示频散数据的射线路径数在5～40 s周期范围内数量最多，且在周期为12～13 s达到峰值，当周期超过45 s后，平均相速度值出现了明显的上下跳动，这可能是由于长周期数据的不足导致.本文最终截取周期范围为5～40 s的相速度频散值作为三维横波速度反演的数据.不同周期的射线分布情况(图7)显示研究区域内部的路径方位覆盖较好，截取的各周期均得到了充足的频散曲线.

图5 (a)台站对的路径分布和(b)筛选后瑞利面波相速度频散曲线Fig.5 (a) The path distribution of the station pairs and (b) phase velocity dispersion curve of the Rayleigh surface wave after sifting

图6 基阶瑞利面波相速度频散数据和射线路径数随周期的分布图圆点为所有台站所有周期的相速度频散数据，实线表示所有台站对的平均相速度与周期的分布图，虚线表示基阶瑞利面波频散数据在不同周期射线路径数量.Fig.6 The distribution of dispersion data of fundamental-order Rayleigh surface-wave phase velocity and the number of ray paths with period The dots represent the phase velocity dispersion data of all stations for all periods. The solid line represents the distribution diagram of the average phase velocity versus period of all station pairs. The dashed line represents the number of ray paths of the fundamental-order Rayleigh surface-wave dispersion data at different periods.

图7 不同周期瑞利面波射线路径分布图Fig.7 The distribution of Rayleigh surface-wave ray paths with different periods

1.4 初始一维横波速度模型的选定

由于本文采用面波频散数据直接反演三维横波速度的方法(Fang et al., 2015)，其三维反演结果对初始一维横波速度模型具有较强的依赖性，这是由于三维反演结果是在一维速度模型基础上做微扰动.因此，选择合适的初始一维模型，对于提高反演结果的精度至关重要.由于不同区域地质构造差异较大，全球或大尺度一维速度模型对小尺度研究区域并不是最佳的.为了尽量减小反演结果对初始一维速度模型的依赖，本文首先对适合本研究区域的一维速度模型进行了严格筛选.经反复测算和综合分析，最终根据研究区域内的相关研究结果(赵珠和张润生，1987；Wang et al., 2007; 黄晓萍等，2012；王同军等，2017；朱介寿等，2017；李敏娟等，2018；Yao et al., 2019; 宋倩等，2020)，并利用介质物理参数之间的经验转换公式(Brocher, 2005)，构建了8种不同的一维横波速度模型，如图8a所示.将上述模型分别作为初次迭代的一维横波速度模型，利用传统面波相速度计算方法(Herrmann, 2013)和观测平均相速度频散曲线(图6中的实线)反演得到相应的8个输出一维横波速度模型.根据正演获得的相速度值与平均观测相速度值的均方根(RMS)残差分布(图8b)，选择相速度值均方根残差最小(数值为1.9 m·s-1)的模型3反演后的一维横波速度模型作为三维面波相速度直接反演法的初始一维横波速度模型(图9a).

图8 研究区域内一维横波速度分析结果(a) 不同颜色实线代表不同模型，其中模型2、4、5可以直接在对应的文献中获得一维横波速度模型，而其他一维横波速度模型是通过经验公式转换而获得； (b) 不同一维速度模型的相速度均方根(RMS)残差分布.Fig.8 The analyzed results of one-dimensional shear wave velocity in the study area(a) Different colored solid lines represent different models, where models 2, 4, and 5 can be obtained directly from corresponding literatures for one-dimensional shear wave velocity models, while the other one-dimensional shear wave velocity models are obtained by conversion of empirical formula; (b) The distribution of root mean square (RMS) residual for phase velocity in the different one-dimensional velocity models.

图9 (a)初始一维横波速度模型，(b)不同周期基阶瑞利面波相速度对横波速度的深度敏感度曲线Fig.9 (a) The initial one-dimensional shear wave velocity model and (b) the depth sensitivity curves of fundamental-order Rayleigh surface-wave phase velocity at different periods to shear wave velocity

1.5 反演方法

本文采用面波频散数据直接反演法(Fang et al., 2015)开展研究.与传统面波反演法相比，该方法无需构建二维相速度分布图，而是直接基于所有面波频散走时数据反演获得三维横波速度结构.值得一提的是，该方法使用快速行进法(Rawlinson and Sambridge, 2005)计算某一周期的瑞利面波相走时和射线路径，从而考虑了复杂的地下速度结构引起面波射线弯曲的效应.

一般情况下，基阶瑞利面波相速度对其波长约1/3左右深度的横波速度最为敏感，因此，周期5～20 s的瑞利面波主要对0～30 km左右的横波速度结构敏感(图9b)，近似于研究区域的上、中地壳分布范围，周期30～40 s的瑞利面波的敏感深度在40～80 km左右，接近于研究区域下地壳和上地幔顶部横波速度结构.本文的频散数据周期范围为5～40 s，能够很好地反映深度60 km以内(地壳和上地幔顶部)的横波速度结构.

为了提高反演结果的稳定性，在反演中引入光滑系数来权衡目标函数中数据残差大小和模型光滑度，因此，选择合适的光滑系数对于三维反演结果非常重要.本文通过吉洪诺夫正则化法(Tikhonov regularization，Aster et al., 2013)测试了多个光滑系数，获得了数据残差归一化范数和模型方差归一化范数相关关系的L曲线(图10)，最终选择光滑系数为20用于反演计算.

图10 选择最佳平滑系数的L曲线Fig.10 Selecting the L curve with the best smoothing coefficient

2 分辨率测试

为了测试数据的空间分辨能力，本文进行了不同网格间隔的检测板测试(Humphreys and Clayton, 1988; Zhao et al., 1992).首先，将初始的一维横波速度模型(图9a)以等间隔加±6%的速度扰动构成三维棋盘格的模型；然后，根据观测数据所对应的路径计算台站对之间相速度对应的理论走时数据，并将随机高斯噪声添加到理论走时数据中，得到合成理论数据集；最后，采用与实际反演中相同的射线路径分布（图7）和相同的反演参数进行反演．

本文通过大量的检测板测试，分别设置横向分辨率为0.25°×0.25°、0.3°×0.3°、0.4°×0.4°、0.5°×0.5°、1°×1°，最终确定本研究区域横向分辨率为0.5°×0.5°(图11).检测板测试的结果表明本文用于反演的频散数据对研究区域60 km以上深度有很好的空间分辨力，大部分节点上的速度扰动量绝对值均恢复到6%，能够分辨水平方向0.5°×0.5°的速度异常.考虑到边缘区域的恢复效果有所下降，因此在之后的三维反演结果中将截掉未出现正负间隔扰动的边缘区域.

图11 检测板测试横剖面结果Fig.11 The results of cross section of the checkerboard tests

3 三维横波速度结构

图12显示了青藏高原东缘0～60 km深度范围的横波速度结构.其结果表明：青藏高原东缘地下横波速度结构具有显著的横向不均匀性特征，且地下浅部的横波速度结构与地表地质单元体有非常好的对应关系.在0～10 km深度层中(图12a)，鲜水河断裂带(F6)南端与龙门山断裂带(F5)交汇处呈现局部高速特征，可能与峨眉山大型火成岩省有关(图12a中红色框区域).龙门山断裂带(F5)两侧速度变化强烈，西侧青藏高原为高速特征，东侧四川盆地为低速特征；研究区域东北部的鄂尔多斯地块也呈现出非常明显的低速特征.这是由于四川盆地和鄂尔多斯地块普遍存在较厚的沉积层，相比于高原地区，沉积层岩石通常是由碳酸盐类(沉积岩)组成，而高原浅层物质主要是由地幔中侵入地壳的铁镁质基性岩(岩浆岩)组成.通常情况下，沉积岩中的地震波速度明显低于岩浆岩中的地震波速度，因此，研究区域西部浅层横波速度相较于东部更高.四川盆地和鄂尔多斯地块局部区域均呈现低速特征，但低速异常特征存在明显差异，四川盆地横波速度是研究区域的最低值，约为2.7 km·s-1，而鄂尔多斯地块的横波速度最低值高于四川盆地，约为3.0 km·s-1，这可能是由于鄂尔多斯地块是非常稳定的克拉通地块，沉积层的压实作用比较强.

图12 不同深度层横波速度分布图黑色实线表示主要活动断裂，空心圆表示2000—2021年之间MS≥4.0地震在不同深度层的分布情况.图(a)中白色虚线分别表示峨眉山大型火成岩省(ELIP)中间和外部带，红色框表示高速异常，其他字母和符号表示同图1中.Fig.12 The distribution of shear wave velocity at different depths layer The black solid lines represent major active faults, the hollow circles represent the distribution of earthquakes with MS≥4.0 at different depths layer between 2000 and 2021. The white dashed lines represent the middle and outer zones of Emeishan Large Igneous Province (ELIP), and the red box represents high-speed anomalies in (a). Other letters and symbols are the same as Fig.1.

在10～20 km深度层中(图12b)，龙门山断裂带(F5)和六盘山断裂(F1)表现为高低速特征的分界线，且地震多分布在龙门山断裂(F5)带附近.巴颜喀拉地块开始出现整体的低速特征，而四川盆地低速区增大，其分布特征呈倾斜的倒“U”型，鄂尔多斯地块的低速特征也进一步向四周扩展.

在20～30 km深度层中(图12c)，速度结构特征与10～20 km深度层相比出现了明显变化，整体上呈现西低东高的特征，基本反映了研究区中、下地壳深度的横波速度结构.巴颜喀拉地块的低速特征进一步增大，速度值也明显降低，尤其是九寨沟地区横波速度值非常低，约为3.2 km·s-1.而四川盆地则呈现了与10～20 km深度层截然相反的速度分布，由浅地表的低速异常转变为明显的高速异常，盆地北部高速异常达到3.7 km·s-1.鄂尔多斯地块的低速区渐渐减小，其南部呈现高速异常.

与上一层相比，在30～40 km深度层(图12d)，研究区域横波速度呈现西北低、东南高的特征.巴颜喀拉地块的低速特征区域基本覆盖了整个地块，向北扩展跨过东昆仑断裂带(F3)，其中北东方向的低速体接近西秦岭断裂带(F2)，而北西方向的低速体分布范围在柴达木地块下方迅速变窄，同时速度值增大，与祁连地块下方低速体具有一定的连通性.四川盆地的高速区域已扩展至整个盆地.

在40～60 km深度层上(图12e，f)，研究区域的横波速度分布的整体特征与30～40 km很相似，横波速度结构依旧是东南部高、西部低分布特征，这主要受研究区域地壳厚度控制.东南部的四川盆地已进入上地幔深度范围，而西部的青藏高原地区还处于地壳范围，形成了研究区这种东南高、西低的速度特征.与龙门山断裂带(F5)呈现为明显的高低速转换带不同，鲜水河断裂带(F6)并未截断巴颜喀拉地块与其南部的川滇地块的低速异常.鄂尔多斯地块呈现的高速异常区域逐渐扩展至整个地块，并在50～60 km向南与四川盆地的高速异常区连为一体.

为了进一步分析研究区域横波速度结构随深度的纵向变化特征，本文给出了沿不同方向纵剖面的横波速度成像结果(图13)，并将剖面两侧各30 km范围内的地震震源投影在该剖面上.在研究区域北端东西向的剖面AA′(图13a，g)中，西侧的柴达木地块和祁连地块0～15 km呈现低速特征，这与柴达木盆地东侧沉积层有关，研究区域东侧的低速特征下边界约束了鄂尔多斯地块沉积层厚度.位于研究区域中部东西向的剖面BB′(图13b，g)，西侧的巴颜喀拉地块是青藏高原的一个组成部分，在30～45 km深度范围内呈现非常显著的水平低速层的特征，并向东延伸.值得注意的是，该低速层受坚硬的华南地块北部西秦岭造山带的阻挡作用，截止于塔藏断裂(F4)附近.塔藏断裂(F4)东西两侧地壳厚度也存在明显的差异，西侧约50 km，而东侧约为45 km.郑晨等(2018)联合面波频散与接收函数反演也获得了同样的地壳厚度.剧烈的构造变形通常伴随强震的发生，2017年MS7.0九寨沟地震及其余震即发生在该高低速分界线附近，Riaz等(2019)对九寨沟地震震源区速度结构的研究也显示了同样的震源分布特征.在研究区域南端东西向的剖面CC′(图13c，g)中，巴颜喀拉地块与四川盆地的边界线——龙门山断裂带(F5)为非常明显的高低速分界线，显示出两个地块间地壳厚度的显著差异，这与塔藏断裂(F4)附近的结构非常相似.此外，2008年5月12日MS8.0汶川大地震及其余震也发震于龙门山断裂带，该断裂带东侧的四川盆地浅部约15 km厚度的低速层特征显示了沉积层厚度.根据已有的研究结果，青藏高原东缘壳幔边界的横波速度约为4.0 km·s-1(Yang et al., 2012)，将这一壳幔边界横波速度值与本文的研究结果结合，表明四川盆地莫霍界面的深度约为40 km，相关研究显示四川盆地的地壳厚度分布范围为35～45 km(Li H Y et al., 2009)，涵盖了本文的研究结果.值得注意的是，龙门山断裂带西侧的巴颜喀拉地块在深度为30～45 km处显示出较为明显的约为15 km厚的水平低速层，其地壳厚度明显高于40 km.

图13 研究区域内不同方向剖面的横波速度纵剖面图图(a)，(b)，(c)，(d)，(e)，(f)的剖面位置分别对应图(g)中红色剖面线，图(g)中黑色实线表示主要活动断裂，剖面图中黑色空心圆表示2000—2021年之间MS≥4.0地震分布；图(c)中的白色五角星表示2008年5月12日MS8.0汶川地震，其余图中的白色五角星均表示2017年8月8日MS7.0九寨沟地震；剖面中深度与水平向的比例尺为2∶7，而地形图高程与水平向的比例尺为1∶10；其余图中字母符号意义均同图1.Fig.13 The profiles of shear wave velocity for different directions in the study area The locations of the profiles in figures (a), (b), (c), (d), (e) and (f) correspond to the red profile lines in figure (g). The black solid lines in figure (g) represent the major active faults. The black hollow circles in the profile figure indicate the distribution of earthquakes with MS≥4.0 between 2000 and 2021. The white stars in figure (c) represent the MS8.0 Wenchuan earthquake on May 12, 2008, and the white stars in the other figures represent the MS7.0 Jiuzhaigou earthquake on August 8, 2017. The scale of depth to horizontal in the cross section is 2∶7, while the scale of elevation to horizontal in the topographic map is 1∶10. The other symbols and letters in the figure are the same as those in Fig.1.

由BB′和CC′这两条东西向的剖面(图13b，c)可以发现，巴颜喀拉地块下方30～45 km深度内存在显著的水平层状低速异常，但向东延伸的范围从北向南有明显缩短，由北侧水平向约为400 km(图13b)逐渐缩短到南端300 km的尺度(图13c).这可能是青藏高原东缘地下软弱物质向东运移过程中受到四川盆地的阻挡转而向南，导致了这两条剖面低速层分布范围发生变化，这也符合GPS测量相对于华南地块运动的速度场反映的运移情况(张培震，2008).

DD′剖面是横跨祁连地块、柴达木地块和巴颜喀拉地块等构造单元的南北向剖面(图13d，g)，整个剖面下方30～45 km深度范围内存在从南至北贯通的低速层，这也与BB′和CC′两个东西向剖面的结果互相印证.不同地块下方30～45 km深度范围内的低速层厚度不一致，柴达木地块的低速体厚度逐渐向北减薄，在祁连地块下方低速体厚度达到最薄，两个地块之间的低速层存在一定的连通性.图13e的南北向纵剖面EE′结果显示在32.3°N—34.5°N范围内深度约为30～45 km也存在低速层，厚度表现为北薄南厚.由这两个南北向的剖面也可以看出巴颜喀拉地块下方30～45 km深度范围内的低速层由西向东范围逐渐减小的特征.

GG′剖面为SW-NE向(图13f，g)，穿过研究区域的中心点，横跨川滇地块、巴颜喀拉地块、柴达木地块、鄂尔多斯地块等多个构造单元.该剖面西南部的川滇地块和巴颜喀拉地块地下30～45 km深度呈现大范围低速层，在30.5°N—32.2°N范围内低速特征尤为显著.

4 讨论

4.1 青藏高原东缘地壳低速异常体的空间分布特征

青藏高原东缘的深部结构与动力学过程是了解青藏高原构造演化的关键，许多学者用不同研究方法做了大量研究(Zheng et al., 2010, 2016; Li H Y et al., 2012; Yang et al., 2012; Bao et al., 2013, 2020; Xie et al., 2013; Jiang et al., 2014; Deng et al., 2015; Tan et al., 2015; 范文渊等，2015；Li Y H et al., 2017; Wang et al., 2020; 夏思茹等，2021).这些研究结果均表明青藏高原东缘中、下地壳存在明显的低速异常体，但对低速异常体的空间分布特征有不同结论.

本文的背景噪声成像结果表明在川滇地块北部、巴颜喀拉地块和柴达木地块下方广泛存在横波速度低速异常体(图13a—f)，速度值为3.1～3.4 km·s-1，深度为30～45 km，祁连地块下方虽然也存在低速异常体(图13d)，但是其分布范围小，速度值更高.其中川滇地块北部中、下地壳显著的横波低速异常体特征(图13d)在其他学者的研究结果中(Yao et al., 2008; Li et al., 2009)也有出现.巴颜喀拉地块中、下地壳广泛存在的低速体东侧大致以龙门山断裂带(F5)为边界(图13c)，在北侧跨越东昆仑断裂带(F3，图13d)，且巴颜喀拉地块和川滇地块北部低速体并未受到鲜水河断裂带(F6)的影响而被分隔开(Yang et al., 2012).从低速体的分布深度来看，本文的研究结果表明青藏高原东缘的横波低速异常体深度范围为30～45 km(图13b—d，f)，整体深度深于Li等(2014)的研究结果(20～35 km)，与大地电磁测深法(Sun et al., 2020)发现的高导层深度范围较为一致.此外，青藏高原东缘中、下地壳低速体在北部的边界分布极为复杂，一直存有争议.有些学者认为青藏高原东缘中、下地壳低速体未跨越东昆仑断裂带(F3)，且在柴达木地块下方未发现明显的低速体特征(Li et al., 2012, 2014; Zheng et al., 2016)，正径向各向异性的特征(Tan et al., 2015)与该低速体分布特征能够很好的吻合.还有些学者则认为该低速体跨越了东昆仑断裂带(F3)，且继续向北渗透到柴达木地块，但这些研究结果有个共同的特点——低速体紧邻西秦岭断裂带(F2)南侧(Zheng et al., 2010; Jiang et al., 2014; Li Y H et al., 2017; Wang et al., 2020).由于柴达木地块位于本研究区域分辨率很高的区域，因此清晰地刻画了低速体在北部复杂的分布形态.北部的边界分为两部分(图12d)，北偏东的一部分边界接近西秦岭断裂带(F2)(图13e)，北偏西的一部分低速体的宽度在柴达木地块下方迅速变窄，且速度值变大(图13d)，与祁连地块下方低速体具有一定的连通性.祁连地块下方低速体速度值大约为3.5 km·s-1，略高于巴颜喀拉地块低速体的速度值(约3.4 km·s-1)，柴达木地块下方低密度异常(Li H L et al., 2017)也给出了相同的分布范围，因此这可能代表了祁连地块低速体相当于巴颜喀拉地块低速体的早期发育阶段(Tapponnier et al., 2001; Bao et al., 2013).此外，祁连地块中、下地壳的横波低速异常体特征主要分布在西北端，其东南端并无明显的低速异常体特征(图12d)，这一现象也得到了相关研究的佐证(Li X F et al., 2014; Bao et al., 2020).

4.2 青藏高原东缘中、下地壳低速流

关于在青藏高原东缘中、下地壳存在的低速体的成因一直是研究的重点.青藏高原东缘的高地面热流的特征(Hu et al., 2000; Wang, 2001)在一定程度上表明该地区地壳具有部分熔融的条件.此外，实验地球物理表明，地壳中部分熔融体能吸收弹性波能量并造成速度降低(Müller and Raab, 1997; Wang et al., 2020).其他学者研究中也能证实青藏高原地壳中存在部分熔融的现象，比如青藏高原南部中地壳部分熔融(Nelson et al., 1996; Unsworth et al., 2005)、青藏高原中部下地壳和上地幔部分熔融(Hacker et al., 2000, 2014; Meissner et al., 2004)、青藏高原东缘中地壳部分熔融(Yang et al., 2012; Xie et al., 2013)和青藏高原北部中、下地壳部分熔融(Owens and Zandt, 1997; Wei et al., 2001).正是由于青藏高原东缘中、下地壳温度升高以至于发生部分熔融，从而导致巴颜喀拉地块中、下地壳横波速度低于周围的速度.在本文研究区域西南处，Bai等(2010)用大地电磁成像法得到的剖面结果表明，巴颜喀拉地块中、下地壳存在一个高导(低阻)层，是由于部分熔融物质导致，并且这些地壳熔融体在地质时间尺度上可以发生流动.本文沿着该剖面观测到的低速体分布特征(图13f)与高导层的分布特征相似，且巴颜喀拉地块表现出较高的泊松比特征(王兴臣等, 2017)，这些均表明低速体特征可能与地壳物质部分熔融有关，也反映了这部分低速体对应的物质比较软弱，有利于流动.

针对青藏高原东缘地表变形和地壳增厚的现象，有学者提出了一种地球动力学模型——下地壳流动模型(Royden et al.,1997,2008; Clark and Royden, 2000)，通过数值模拟实验，认为如果青藏高原下地壳中软弱的低速物质在发生流动的过程中被坚硬的四川盆地阻挡，反映在地表会形成狭窄陡峭的地形变化，如龙门山断裂带(F5).本文高分辨率三维横波速度结构很好地反映了青藏高原东缘低速体的空间分布.其中，巴颜喀拉地块中、下地壳低速体在东部被四川盆地下地壳高速体阻挡，截止于龙门山断裂带(F5)附近(图13c，e)，且在龙门山断裂带(F5)两侧横向距离变化不到100 km的范围内，地形高程迅速从高原的5 km下降到盆地的几百米(图1).而在巴颜喀拉地块北部，其中、下地壳低速体沿着四川盆地和柴达木盆地之间宽阔的区域继续运移，并未受到坚硬高速体的阻挡，地表地形高度变化非常缓慢，这一特征与Clark和Royden(2000)的研究结果一致，也很好地解释了青藏高原东缘上地壳无显著的缩短而地壳增厚的现象(Burchfiel et al., 1995; Wang et al., 1998).Shapiro等(2004a)从地震径向各向异性分析中，发现青藏高原中部中、下地壳可能存在通道流，这一观点同样也得到在该研究区域内获得的高电导率的论证(Wei et al., 2001).

因此，为了更加清楚地了解青藏高原东缘中、下地壳低速体的运移方向，本文做了横波速度三维示意图(图14)，为图示清楚和重点分析中、下地壳速度结构，图中未绘制0～20 km深度的横波速度.从图14中可明显发现，在研究区域的北部，来自于青藏高原的中、下地壳低速体向东运移过程中逐渐转变为向东北方向运移，且跨越东昆仑断裂带(F3)，到达柴达木地块，形成目前低速体复杂的形态.而在研究区域中南部，青藏高原中、下地壳低速体向东运移的过程中，在靠近龙门山断裂带(F5)的下方受到四川盆地刚性地壳的阻挡作用，出现分流现象：一支沿着鲜水河断裂带(F6)向南运移，进入青藏高原东南缘；而另一支则沿着龙门山断裂带(F5)向北方向运移，并截止于西秦岭断裂带(F2).这样的运移方式与GPS观测到的地表速度场分布特征(张培震，2008)具有很好的一致性.

图14 青藏高原东缘中、下地壳横波速度三维示意图Fig.14 The diagram of three-dimensional shear wave velocity in the middle-lower crust of the eastern margin of the Tibetan Plateau

值得注意的是，祁连地块西北端中、下地壳，存在着较为明显低速特征，其他研究结果(Yang et al., 2012; Bao et al., 2013; Jiang et al., 2014)也给出了同样的特征，但是对于该低速体成因的研究还存在争议.鉴于这是本研究区域的边缘地区，没有足够的分辨率，还需要利用更多的台站波形数据和结合其他研究资料并向北扩大研究区域，才能对低速体的成因做进一步深入的探讨.

4.3 中、下地壳低速体对塔藏断裂附近地下速度结构的影响

青藏高原东缘中、下地壳低速体对塔藏断裂(F4)附近地下速度结构造成了显著的影响，以往的相关研究由于分辨率和研究区域等因素对塔藏断裂(F4)附近地下的速度结构和成因未能给予足够的重视.从图13b、e、f中可以明显发现，塔藏断裂(F4)附近地壳相比于周围呈现低速异常的特征，且在30～45 km深度范围内与巴颜喀拉地块下方的低速体相连接.因此，我们认为这种异常现象可能是由于沿着龙门山断裂带(F5)向北运移的低速物质和来自于青藏高原内部的向东北运移的低速物质受到华南地块北部的西秦岭造山带较硬中、下地壳的阻挡作用，在塔藏断裂(F4)附近沿着裂缝垂直运移到上地壳，导致该断裂地壳内速度相比于周围区域明显偏低(图13f).在东昆仑断裂带(F3)东段的三维电阻率的研究结果(Sun et al., 2020)表明，中、下地壳高电导率层在塔藏断裂(F4)附近出现局部上涌的现象，这从侧面印证了本文研究成果中塔藏断裂(F4)附近地壳巨厚低速体的存在(图13e)，同时也佐证了塔藏断裂(F4)附近上地壳低速特征可能是由于中、下地壳的低速物质上涌导致.正是由于这些低速物质充填进断裂附近地下裂隙后，降低了岩石有效正应力而导致岩石脆性破坏，从而可能诱发了九寨沟地震的发生(Chen et al., 2021).这种中、下地壳低速物质在断裂带附近局部上涌的异常现象在2008年汶川MS8.0地震震源区速度结构的研究中存在(Wang et al., 2014).

此外，塔藏断裂(F4)附近莫霍界面出现明显的下沉(图13e)，应是低速体在塔藏断裂(F4)处上涌和物质的积累导致局部地形的上升和地壳的增厚，由于重量的增加，造成物质向深部有下沉作用，使得塔藏断裂(F4)附近的莫霍界面出现了下凹的现象.

5 结论

本次研究根据台站分布和研究目的选取99°E—109°E、28°N—38°N为研究区域，利用四川、重庆、甘肃、青海、陕西地震台网的105个宽频带固定地震台站在2017年8月—2018年8月记录的连续垂直分量背景噪声数据，通过互相关计算和时频分析，并对频散曲线的质量进行严格控制和筛选，共获得5416条周期为5～40 s的基阶瑞利面波相速度频散曲线.利用面波频散数据直接反演法(Fang et al., 2015)获得了青藏高原东缘地壳和上地幔顶部高精度的三维横波速度结构，主要结论如下：

(1)青藏高原东缘地下浅部的速度结构与地表地质单元体有非常好的对应关系.四川盆地和鄂尔多斯地块的低速特征主要与沉积岩有关；而青藏高原东部则表现为与铁镁质基性岩相关的高速特征.

(2)巴颜喀拉地块中、下地壳存在大范围的低速体，可能由于地壳内物质发生部分熔融，吸收弹性波能量并导致速度降低；而柴达木地块下地壳局部也存在低速体，应为巴颜喀拉地块低速体的延伸部分.

(3)青藏高原东缘中、下地壳横波低速体分布特征概括为：在北部，来自于高原中、下地壳低速体向东运移过程中逐渐转变为向东北方向运移，在跨越东昆仑断裂带后低速体范围迅速变窄，而速度值逐渐增大，与祁连地块的低速体具有一定的连通性，且该低速体速度值略高于巴颜喀拉地块低速体的速度值，这可能意味着祁连地块低速体相当于巴颜喀拉地块低速体的早期发育阶段；在中南部，中、下地壳低速体受到四川盆地刚性地壳的阻挡作用，在龙门山断裂带附近出现不同方向的运移现象，一支沿着鲜水河断裂带向南方向运移，而另一支则沿着龙门山断裂带向北方向运移，并截止于西秦岭断裂带.

(4)塔藏断裂附近地壳内的低速体分布特征可能是由于青藏高原东缘中、下地壳低速体运移到该断裂处，受到华南地块北部的西秦岭造山带较硬的中、下地壳的阻挡作用沿着裂隙垂直运移到上地壳，导致速度低于周围区域，同时由于重量的增加使得莫霍界面下沉.

致谢中国地震台网中心为本研究提供了所用的数据资料，中国科学技术大学姚华建教授团队提供了背景噪声成像法的程序包(Yao_trainingpackage)和面波直接反演法程序(DSurfTomo)，本研究使用了圣路易斯大学地震中心开发并维护的地震学软件包CPS330(Herrmann, 2013)，两位审稿专家为本文提出了宝贵的建议，文中图1、图3、图5—图13使用GMT软件(Wessel et al., 2013)制作，图2和图4使用matlab软件制作，图14使用ParaView软件(Ayachit, 2015)制作，在此一并表示衷心感谢.