青藏高原东北缘地壳及上地幔顶部速度结构研究

夏思茹, 石磊, 李永华, 郭良辉

1 中国地震局地球物理研究所, 北京 100081 2 中国地质大学(北京), 地球物理与信息技术学院 , 北京 100083 3 中国地震局震源物理重点实验室, 北京 100081

0 引言

青藏高原东北缘是青藏高原向北东方向扩展的前缘部位，晚新生代构造变形强烈、地震活动频繁(Meyer et al., 1998; Tapponnier et al., 2001)(图1a)，不仅是研究大陆动力学演化的热点区域，也是研究大陆强震孕育相关过程的理想场所.前人在青藏高原东北缘及邻区开展了大量的地球物理探测与研究工作(Zhang et al., 2010; 孟小红等, 2012; Xiao et al., 2016; Feng et al., 2020)，相关成果不仅推进了对青藏高原隆升过程与岩石圈变形机制的认识，也为认识强震孕育环境及其发生机理提供了重要的深部资料.

但由于不同研究给出的结果不尽一致，部分深部构造问题并未取得共识.如，深地震测深剖面探测(张先康等, 2008; 王帅军等, 2019)及接收函数(李永华等, 2006; Wang et al., 2016)等研究表明，青藏高原东北缘地壳厚度变化剧烈，整体地壳厚度介于42～65 km之间，其中松潘—甘孜地块地壳厚度约为47～55 km，西秦岭造山带地壳厚度介于42～60 km之间，祁连构造带地壳厚度较厚，约为44～65 km.地壳平均波速比略低于全球平均波速比，地壳物质组成以长英质为主，不存在大规模熔融.地震体波走时成像显示青藏高原东北缘存在明显壳内低速层(肖卓和高原, 2017; 李敏娟等, 2018; 莘海亮等, 2020; 瞿辰等, 2020)，但是由于体波成像垂向分辨率较低，关于低速层的深度分布仍存在争议(Guo et al., 2019; Xin et al., 2019).此外，体波成像研究主要基于线性反演展开，对初始模型依赖明显，不同研究给出了相对速度异常，但不同研究结果的绝对速度差异很难比较(Ye et al., 2016; Sun et al., 2019).地震面波和噪声成像研究显示，松潘—甘孜地块和祁连构造带下方地壳呈明显的S波低速异常，但关于壳内低速是否连续分布、及其成因仍存在争议(Bao et al., 2013; Li et al., 2017; 付媛媛和肖卓, 2020).

前人利用体波走时与面波联合反演，很好地描绘了壳幔的速度结构(Fang et al., 2016; Guo et al., 2018).本文基于面波和噪声成像结果(Xie et al., 2013; Shen et al., 2016)构建三维初始模型，采用双差层析成像方法(Zhang and Thurber, 2003, 2006)反演青藏高原东北缘地区固定和流动地震台站的P、S波走时资料，获得了研究区可以解释体波走时与面波观测资料的较高分辨率的三维P、S波速度结构和地震重定位结果，并对速度模型的可能地质含义及研究区地震的孕育背景进行了分析与探讨.

1 数据与方法

本文使用的地震走时数据分为两部分，第一部分是由中国地震台网中心提供的自2008年1月至2019年12月71个固定台站记录的P、S波震相报告；第二部分是利用“中国地震科学探测台阵-南北地震带北段”418个流动台站(图1b)于2013年9月至2016年3月记录的近震波形资料，根据中国地震台网中心提供的地震目录截取地震波形，手动拾取1.5级以上地震事件的P、S波到时.本研究采用双差走时成像算法(TomoDD)开展体波走时反演(Zhang and Thurber, 2003, 2006)，此方法将双差定位(Waldhauser and Ellsworth, 2000)和绝对走时层析成像方法(Thurber, 1983)相结合，同时反演单一地震事件的绝对走时数据和两个相邻事件的相对走时数据，以实现对三维速度结构与震源位置的联合反演.假设两个相邻地震事件i与j之间的距离远小于这两个事件到台站k的距离，可以默认两个地震事件的传播路径相似，将两个事件的走时残差相减即为双差.由于地震位置与到时的关系是非线性的，因此需要用Taylor级数展开的方法使其线性化，一般将地下速度结构划分为三维网格节点(x1,x2,x3)，因此表示为如下线性关系式：

图1 (a) 研究区地形与构造背景， (b) 研究区台站位置分布图 黑色实线为构造边界(Taylor and Yinn, 2009)，F1:昆仑断裂带，F2:西秦岭北缘断裂带，F3:海原断裂带，F4:祁连北缘断裂带，Songpan-Garzê：松潘—甘孜地块，West Qinling：西秦岭造山带，Qilian：祁连造山带，Alxa：阿拉善地块；红色五角星为甘肃岷彰6.7级地震，绿色五角星为九寨沟7.0级地震，蓝色五角星为玉树7.3级地震.图(a)中蓝色圆点为地震事件，图(b)中蓝色方形表示固定台站，黑色三角星表示 流动台站，白色虚线表示图8中垂直切面位置.Fig.1 Topography and tectonic background (a) and distribution of seismic stations (b) of the research area Black solid lines are tectonic boundary (Taylor and Yinn, 2009), F1: Kunlun Faults, F2: West Qinling Faults, F3: Haiyuan Faults, F4: Qilian northern marginal Faults, Songpan-Garzê: Songpan-Garzê Block, West Qinling: West Qinling Orogen, Qilian: Qilian Orogenic, Alxa: Alxa Block, red star is the epicenter of Minzhang MS6.7, green star is the epicenter of Jiuzhaigou MS7.0, blue star is the epicenter of Yushu MS7.3. Blue dots in Fig(a) are seismic events, blue squares and black triangles in Fig(b) represent the fixed and mobile stations, respectively. White dashed lines represent the position of vertical sections in Fig.8.

(1)

为了确保反演结果的可靠性，反演过程中要求每个事件被至少5个台站记录到，共获取地震事件25454个.本文根据地震波走时时距曲线挑选初至震相(图2)，包括震中距≤200 km的Pg、Sg波震相，和震中距>200 km的Pn、Sn波震相.用于反演的Pg和Pn波到时共计179584条，Sg和Sn波到时共计145131条.反演过程中将Pg和Pn波权重设置为1，Sg和Sn波的权重设置为0.5，以此进行联合反演.地震事件对与台站之间的最大距离为500 km，事件对之间的最小距离为0.1 km，最大距离为50 km，其中每一个事件对所需要的震相数的最小值设为8，最大值设为120，且将每个地震事件的最大邻居数设为50，剔除数据中存在的明显异常，最终用于反演的地震事件19424个.

图2 地震波走时时距曲线 黑色圆点为P波震相，灰色圆点为S波震相.Fig.2 Travel-time curves for P and S waves The black dots are P wave phases, the gray dots are S wave phases.

2 初始模型构建与参数选取

双差地震成像反演基于线性反演开展，对初始模型有较强的依赖性，初始模型与实际结果偏差较大时，有时经过多次迭代反演结果也不收敛.研究中我们将背景噪声和面波成像得到的S波速度结果(Xie et al., 2013; Shen et al., 2016)，利用Brocher(2005)给出的VP-VS转换关系式(公式2)，得到初始P波速度模型.初始模型水平网格间距为0.5°×0.5°，深度方向上0～60 km以5 km为间隔共计13层.

(2)

在双差层析成像反演计算中，光滑权重与阻尼因子的取值对反演结果的稳定性有直接的影响.本文采用折中曲线来确定这两个正则化参数(图3a,b).本文测试表明，当光滑权重为100，阻尼值为400时，既可以使得速度模型不过于复杂，又可以拟合本研究的大部分观测数据.此外，由于双差反演计算采用的是快速稳定的最小平方QR分解(LSQR)迭代算法，该算法没有明确的收敛条件，因此需要对残差的变化趋势判断是否继续迭代，经过多次试验，当迭代次数为10次时，残差变化趋势曲线趋于平缓，因此本文最终选取的迭代次数为10次.

图3 (a) 光滑权重均衡曲线； (b) 阻尼因子均衡曲线； (c) 不同迭代次数残差变化趋势Fig.3 (a) The optimum smoothing parameter selected by trade-off curve, (b) The optimum damping factor parameter selected by trade-off curve, (c) Retrieval of residual variation trend by travel time

图4给出反演前的残差与反演后的残差分布直方图，可以看出，反演前P波残差主要分布±2 s范围内，均方根为0.81 s；反演后在±1 s范围内，均方根降为0.29 s；反演前S波残差主要分布在±3 s范围内，均方根为1.17 s；反演后在±1.6 s范围内，均方根降为0.51 s.

图4 (a) P波反演前、后残差分布； (b) S波反演前、后残差分布 灰色实框为反演前残差，黑色方框为反演后残差.Fig.4 (a) P waves residual distribution before and after inversion, (b) S-wave residual distribution before and after inversion, the gray solid boxes are residuals before inversion, and the black hollow boxes are residuals after inversion.

3 分辨率测试

本文采用棋盘测试的方法，在初始模型基础上加入±5%相间的速度扰动，再使用与实际反演相同的射线分布和反演参数进行反演，根据模型的恢复效果来评价分辨能力.图5为部分深度棋盘测试结果.本研究中P波走时数据总量多于S波，因此P波成像分辨效果优于S波成像结果.整个反演深度范围(0～60 km)，除30～40 km深度外(多数地震分布在25 km以浅)，P波速度成像横向分辨率可达0.5°×0.5°.0～20 km深度S波成像横向分辨率与P波成像分辨率相当，约为0.5°×0.5°，由于Sn走时数据较少，40 km以深横向分辨率降低为1.0°×1.0°.本研究反演所用的Pn和Sn震相均为折射波震相，所形成的地震射线经Moho面传至地表，在60 km深度地震射线密度较高，因此60 km的恢复效果优于40 km.

图5 P、S波成像分辨率测试 P波模型水平网格间距为0.5°×0.5°(a—d)，10 km深度和20 km深度的S波模型水平网格间距为0.5°×0.5°(e—f)，40 km和60 km 深度的S波模型水平网格间距为1.0°×1.0°(g—h). P、S波成像检测测试垂向间距为5～10 km (i—l).Fig.5 The checkerboard resolution test at different depths for P and S wave tomographic imaging The P wave velocity model has been parameterized with a grid spacing of 0.5°×0.5° in the horizontal direction (a—d), the S wave velocity model has been parameterized with a grid spacing of 0.5°×0.5° in the 10 km、20 km depth (e—f) and 1.0°×1.0° for 40～60 km depth in the horizontal direction (g—h). P- and S-wave checkerboard anomaly test results in vertical sections, with a grid of 5～10 km (i—l).

4 结果

4.1 地震事件重定位

本文对研究区2008—2019年19424个1.5级以上地震进行了重定位，重定位后位于阿拉善地块与祁连构造带的地震事件逐渐收敛到海原断裂带，松潘—甘孜地块与西秦岭造山带内的地震分布更加集中(图6a).图6b为中国台网统一地震目录的事件深度，受区域台网分布间距较大且速度模型简单的影响，主要集中在5～15 km范围内，占事件总数的55%.重定位后，整体的震源分布有变深的趋势，0～15 km深度事件较为均匀分布(图6c).前人研究(高见等, 2013; 李敏娟等, 2018; 黄焱羚等, 2020)也显示重定位后地震事件呈现出收敛到断裂带周缘的特征，研究区地震震源深度主要发生在20 km以浅.

图6 (a) 重定位前、后震中分布； (b) 重定位前震源深度； (c) 重定位后震源深度Fig.6 (a) The distribution of epicenters before(blue bots) and after(orange bots) relocation； (b) The epicenters depth before relocation； (c) The epicenters depth after relocation

4.2 P和S波速度分布

本文通过联合反演获得了青藏高原东北缘及邻区地下0～60 km深度范围内的P、S波速度结果，展示了5 km、10 km、15 km、20 km、30 km、50 km深度的P、S波速度分布(图7).在上地壳深度范围，沉积层较厚的地区(如阿拉善、柴达木和松潘—甘孜地块)表现为P、S波低速异常，而沉积层缺失的山区则呈现为高速异常，这一点在以往的地震成像研究中也有发现(Li et al., 2017).其中5～10 km深度范围，柴达木和松潘—甘孜地块都表现为低速异常，而阿拉善地块在5 km深度表现为低速异常，10 km深度则呈现为高速异常，这暗示阿拉善地块的沉积层厚度较柴达木盆地和松潘—甘孜地块要薄.深地震测深剖面研究结果显示，阿拉善地块沉积层厚度约4～5 km(郭文斌等, 2016)，与本研究结果一致.

图7 不同深度水平切片 黑色实线为构造边界，红色五角星形为岷漳6.7级地震，绿色五角星形为九寨沟7.0级地震.Fig.7 Velocity images of different layers The black solid lines are structural boundary, the red star is the epicenter of Minzhang MS6.7 earthquake, the green star is the epicenter of Jiuzhaigou MS7.0 earthquake.

15～30 km深度范围，阿拉善和柴达木等地块表现为P、S波高速异常，而青藏高原东北缘其它地块则呈现低速异常，这与前人采用地震体波和面波成像得到的结果相同(Li et al., 2014; Jiang et al., 2014; Wei et al., 2017).40 km深度及以下，阿拉善地块表现为P、S波高速异常，而青藏高原东北缘各地块则整体表现为低速异常，这主要是因为在该深度范围，青藏东北缘地区属于地壳范围，而阿拉善地块则已经进入上地幔顶部.这与接收函数研究得到的阿拉善地块地壳较薄(～40 km)，青藏东北缘地壳厚度较厚(42～65 km)的观测结果相符(王兴臣等, 2017).

从垂向剖面来看，祁连构造带和松潘—甘孜地块下方中下地壳范围存在明显的P、S波低速带，但这两个壳内低速带并不相互连通.这一观测结果与前人的面波成像结果(Li et al., 2017; Shen et al., 2016)相一致.此外，本文研究还显示，松潘—甘孜地块的壳内低速带P、S波速度分别为5.5 km·s-1和3.2 km·s-1，较祁连构造带分别要高0.1 km·s-1和0.2 km·s-1.

图8 典型剖面垂直切片(剖面位置标示于图1b) 红色五角星形为岷漳6.7级地震，绿色五角星形为九寨沟7.0级地震.Fig.8 Velocity images of typical vertical sections The red star is the epicenter of Minzhang MS6.7 earthquake, the green star is the epicenter of Jiuzhaigou MS7.0 earthquake.

岷漳6.7级地震与九寨沟7.0级地震震源区地处不同块体的边界，两侧地壳P、S速度特征明显不同.这种地震分布在高、低速过渡带的现象在其他地区也有发现(Hauksson and Haase, 1997; Zhao et al., 2000; Huang et al., 2002; Lei and Zhao, 2009; 王长在等, 2011).研究区的壳内低速层可能处于部分熔融或易于蠕变的状态，脆性的上地壳部分更容易积累应变能，从而导致地震的发生.

5 讨论和分析

5.1 与已有模型比较

为验证结果的可靠性，利用Herrmann(2013)的程序分别基于Crust1.0模型、Litho1.0模型和本研究获得的速度模型合成了对应的Rayleigh波相速度频散，并与Shen 等(2016)利用噪声和面波成像方法得到的Rayleigh波相速度频散结果进行对比.我们任意选取了研究区中2个点的速度模型正演对比结果进行示例(图9).可以看出，基于Crust1.0模型、Litho1.0模型合成的Rayleigh波相速度与Shen 等(2016)利用噪声和面波成像方法得到的相速度存在较大差异，这表明Crust1.0模型和Litho1.0模型无法反映研究区真实的地下结构.基于本研究速度模型获得的合成频散与Shen 等(2016)得到的Rayleigh波频散较为一致，均在±0.1 km·s-1误差范围内.本研究反演得到的VP和VS模型不仅能很好地解释体波走时，也可以很好的拟合面波观测数据，这也表明本研究得到速度模型较Crust1.0和Litho1.0模型更为可靠.

图9 观测与理论面波相速度频散比较 带误差棒的绿色点为相速度观测值(Shen等，2016)，红色、蓝色和黄色实线分别代表 基于本研究模型、Crust1.0模型和Litho1.0模型的正演频散.Fig.9 Comparison of observed and predicted Rayleigh-wave phase velocity dispersion curves at three selected grids The green dots represent the observed dispersions (Shen et al., 2016). The red, blue and yellow lines represent synthetic dispersions based on our model, Crust1.0 model, and Litho1.0 model, respectively.

5.2 壳内低速层分布及成因

本文研究结果(图7)显示松潘—甘孜地块及祁连构造带中下地壳表现为显著的P、S波低速异常.这与前人采用背景噪声和面波成像(Yang et al., 2012; Li et al., 2014，2017; Bao et al., 2013; Zheng et al., 2016;王琼和高原, 2018)、体波层析成像(肖卓等, 2017; Sun et al., 2019; Guo et al., 2019)及深地震测深剖面探测(Ye et al., 2015; Zhang et al., 2013)得到的研究结果相一致.祁连构造带低速层主要分布在西北部，P波速度低速层深度范围为10～40 km，S波速度低速层主要分布在20～40 km深度.由B-B′剖面(图8c)可进一步看出，松潘—甘孜地块在5～20 km和20～40 km存在双层P波低速层，大致以100°E为界，东西两侧分别存在低速层，这两个低速层并不联通.

目前关于松潘—甘孜地块壳内低速层的成因及地壳流是否存在是讨论的热点.背景噪声成像反演显示松潘—甘孜地块存在低速层，认为低速层形成可能与地壳流有关(范文渊等, 2015; Tan et al., 2015; Wei et al., 2017).但接收函数(李永华等, 2006; Pan and Niu, 2011; 王兴臣等, 2017)与宽角反射、折射地震剖面(Gao et al., 2013；Zhang et al., 2013)研究表明，青藏高原东北缘泊松比及地震波速度均较低，不支持下地壳流的存在.我们的结果显示松潘甘孜地块壳内低速层的S波速度小于3.2 km·s-1，这与前人采用背景噪声层析成像方法所得到松潘—甘孜地块整体S波速度较低的结论相一致，他们据此推测这种显著的壳内低速异常可能与地壳部分熔融有关(Yang et al., 2012; Xie et al., 2013; Li et al., 2014; 潘佳铁等, 2017; Wei et al., 2017).大地电磁研究资料表明，研究区电导率较高，存在部分熔融可能(赵凌强等, 2015).地热研究结果显示，松潘—甘孜地块东部中地壳地下温度大致在700～800 ℃，下地壳大致在1000 ℃左右(Jiménez-Munt et al., 2008)，这样的温度足以使得壳内物质发生熔融.

与松潘—甘孜地块相对，祁连构造带的壳内低速层范围较小，其P波速度值低于5.6 km·s-1，S波低速层速度值低于3.4 km·s-1，且低速层主要分布在祁连构造带西北部.肖卓和高原(2017)研究认为祁连构造带上中地壳(10～30 km)存在低速层，本文研究结果显示P波速度低速层深度范围为10～40 km.祁连构造带西北部大地热流值与全球平均大地热流值大致相同(Wang, 2001)，且泊松比较低(王兴臣等, 2017)，因此推断不存在部分熔融(Li et al., 2014, 2017).对于该区低速层产生原因，背景噪声和面波成像研究认为，祁连构造带下的低速体可能是由于地幔岩浆活动底侵作用造成(Li et al., 2017).祁连构造带西北部地形较高，接收函数研究结果显示该区地壳厚度较厚(Yue et al., 2012; Wang et al., 2017)，低速带可能与地壳增厚和地表隆升有关，被认为是壳内响应(Tseng et al., 2009).

5.3 中强地震与速度结构的关系

青藏高原东北缘由于长期受到东北方向的挤压作用，应变能不断积累.自2008年以来，研究区发生4级以上地震138次，5级以上地震29次，6.5级以上地震3次，分别是2010年4月14日青海玉树7.3级地震，2013年7月22日甘肃岷漳6.7级地震和2017年8月8日四川九寨沟7.0级地震.本研究将岷漳6.7级地震与九寨沟7.0级地震重定位后震源位置投影到水平与垂直速度切片上，进一步讨论研究区中强地震与速度结构的关系.

2013年7月22日在甘肃省定西市岷县、漳县交界处发生6.7级地震，中国地震台网中心快速测定的震源位置为东经104.21°、北纬34.54°、震源深度15 km.重定位后震源位置在水平方向没有太大变化，但深度约为8.5 km，与前人研究(孙蒙等, 2015)提供的震源深度基本一致，都位于脆性的上地壳范围.岷漳6.7级地震分布在西秦岭北缘断裂带以南西秦岭地块内，该地震震源以西的西秦岭地块下方有明显的壳内低速带，但其东侧下方则不存在明显的壳内低速异常带(图8a,b).

重定位后的九寨沟7.0级地震震源深度为12 km，与德国地球科学研究中心(GFZ)发布的深度及后续震源参数研究结果(张旭等, 2017; 宋秀青, 2017; 郑国栋等, 2019)类似.赵博等(2018a,b)利用地震波形反演得到的破裂面的矩心深度约为8～9 km，较运动学方法得到的深度要浅，这与易桂喜等(2017)九寨沟地震震源机制研究结果一致.该地震同样处于脆性的上地壳范围，其下伏为明显的壳内低速度带(图8c,d).九寨沟7.0级地震地处昆仑断裂带东段附近，该地震震源位置以南的松潘—甘孜地块存在明显的壳内低速异常，其北侧的西秦岭地块下方则不存在明显的壳内低速异常带.

6 结论

本文利用青藏高原东北缘71个固定台站和手动拾取的“中国地震科学探测台阵-南北地震带北段”418个流动台站所记录到的25454个天然地震事件资料，采用双差层析成像方法对近震走时数据进行反演，获得了研究区0～60 km高分辨率(0.5°×0.5°)的三维P、S波速度结构和重定位结果.

研究结果表明，基于Crust1.0和Litho1.0模型正演所得的频散曲线与实际Rayleigh波相速度频散曲线存在较大误差.本文给出的P、S波速度模型既可以很好的解释体波走时又可以拟合面波相速度观测资料，较已有全球模型更为可靠.

松潘—甘孜地块和祁连构造带下方20～40 km深度范围表现为显著的P、S波低速异常.松潘—甘孜地块壳内低速层可能与地壳部分熔融有关.祁连构造带的壳内低速层主要分布在祁连构造带西北部，可能与地壳增厚有关.

精定位后的岷漳6.7级地震和九寨沟7.0级地震震源深度都位于脆性的上地壳.两个地震震源区地处不同块体的边界，两侧地壳P、S速度特征也都明显不同，一侧存在明显的壳内低速异常，另一侧则不存在明显的壳内低速异常带.震源区的壳内低速层可能处于部分熔融或易于蠕变的状态，脆性的上地壳部分更容易积累应变能，从而导致地震的发生.

致谢感谢中国地震局地球物理研究所“中国地震科学探测台阵数据中心”为本研究提供地震波形数据，感谢中国科学技术大学张海江教授提供TomoDD程序.本文部分图件使用GMT软件绘制(Wessel and Smith, 1998).