据中国地震台网测定，北京时间2016年1月21日01时13分在青海省门源县(37.68°N，101.62°E)发生M_S6.4地震，震源深度10 km.震中及周边地区震感强烈，给当地造成了人员受伤和重大财产损失.震源机制解显示此次地震为逆冲型地震，余震序列沿NW—SE方向分布，结合地震现场考察结果，认为其发震断层可能为冷龙岭断裂和民乐—大马营断裂之间的一条次级断裂(胡朝忠等，2016；梁姗姗等，2017).

门源地区(图 1)位于青藏高原东北缘，长期受到印度板块与欧亚板块碰撞挤压作用的影响，新构造运动十分强烈.区域内发育了一系列活动时代、运动性质、空间展布不同的断裂，其中距离此次地震最近的为震中南侧的冷龙岭断裂，该断裂属于北祁连山活动断裂带的一部分，全长127 km，总体走向为NW60°—NW70°，是一条以左旋走滑特征为主，兼具逆冲分量的全新世活动断裂(何文贵等，2000；Lasserre et al.，2002；袁道阳等，2004).震中北侧为民乐—大马营断裂，断裂呈NWW向延伸，倾向SW，距离震中约15 km，此断裂为晚更新世逆冲断裂，全新世以来无明显活动(陈文彬，2003；郭鹏等，2017).此次地震位于祁连山中东段，是我国地震活动最强烈的地区之一.自1900年以来震中附近多次发生6级以上地震，包括破坏性的1927年M_S8.0古浪地震、1928年M_S6.0古浪地震、1986年M_S6.4门源地震等，地震活动频度高、强度大(Chen et al., 2016).

图 1

Fig. 1

图 1 本研究中使用的台站及震中分布图 图中红色三角形为固定台站，蓝色菱形为流动台阵，黑点为地震的震中位置，红色五角星为MS6.4主震位置，紫色五角星为MS≥7.0的历史地震，红线为主要断层. Fig. 1 Distribution of epicenters and seismic stations used in the study Red triangles stand for fixed seismic stations, blue rhombus stand for temporary seismic stations, black dots represent the epicenter of earthquakes, red star is the epicenter of MS6.4 earthquake, purple stars are the epicenters of historical earthquakes with MS≥7.0, and red lines are major faults.

地震精定位和层析成像是研究地下隐伏断裂、获得深部介质结构的重要方法(王长在等，2011).精定位后的震源位置可以更准确地揭示地震的空间分布特征，刻画断层在地壳深部的展布形态，反映深部的地质构造背景(Lei et al., 2012；曲均浩等，2015).而震源区的三维速度结构能够反映地壳介质属性及变化情况，将其与地震活动的空间分布相结合，可为认识地震发生的孕震环境及孕震机理提供重要的依据(李大虎等，2015).2013年启动的“中国地震科学台阵探测——南北地震带北段(ChinArray Ⅱ)”项目在青藏高原东北缘及周边地区布设了大量的宽频带流动地震台站，平均台间距35 km，为青藏高原东北缘的地球物理研究提供了丰富可靠的观测资料，据此产生了一批具有重要意义的研究成果(常利军等，2016；郭慧丽等，2017；王兴臣等，2017).本文利用双差层析成像方法联合反演门源地区(36°N—39°N，100°E—103°E)三维速度结构和震源参数，研究门源M_S6.4地震孕育与发生的深部构造背景，并分析该地区速度结构和地震活动性之间的关系，为未来该区的地震危险性判断和地震趋势研究提供一定的参考.

2008年以来，中国国家数字地震网络系统、甘肃和青海区域数字地震台网相继建成，这些台站在青藏高原东北缘形成了较密的地震观测台网并产出了海量的观测数据；2013年由中国地震局兰州地震研究所协助承担的“中国地震科学台阵探测——南北地震带北段(ChinArray Ⅱ)”项目在甘肃、青海地区布设了密集的流动台站，开展了为期两年半(2014年1月至2016年6月)的野外台阵观测，这些丰富的地震资料为研究该区的精细速度结构提供了有力的数据支撑.

本文收集了2009年1月至2017年6月甘肃、青海数字地震台网及2014年1月至2016年6月中国地震科学探测台阵记录到的发生在门源地区(35°N—40°N，99°E—104°E)的18566个地震的初至P波和初至S波到时资料.图 1为研究中使用的台站及震中分布图.其中，甘肃、青海数字地震台网的固定台站32个，台站仪器为BBVS-60、CMG-3ESPC、CTS-1、BBVS-120等型号，观测频带为50 Hz~60 s或50 Hz~120 s，采样率均为100 Hz；中国地震科学探测台阵的流动台站54个，地震计为CMG-3ESP，频带范围为50 Hz~60 s，采样率为50 Hz.

为了保证地震数据的可靠性，我们对其进行了检查与筛选.首先根据拟合出的地震波走时曲线去除偏差较大的震相数据(如图 2所示)，并要求参与计算的每个地震至少有8个到时记录.研究区内震中距覆盖范围为4~500 km，在地震对匹配时，限定地震对之间的最大距离为20 km，每个地震最多可以和20个地震组成地震对，最终我们挑选出13270个地震的189249条绝对到时数据(P波98704条，S波90545条)及1681996条相对到时数据(P波872780条，S波809216条)用于门源地区震源参数和层析成像的联合反演，其地震对之间的平均距离为4.93 km.挑选后的地震震中分布如图 1所示.考虑到P波和S波到时拾取的精度不同，反演过程中赋予P波震相的权重为1.0，S波震相的权重为0.6.

图 2

Fig. 2

图 2 初至P波、初至S波走时曲线 为保证数据质量，剔除黑线之外的震相数据. Fig. 2 Time-distance curves of initial P and S waves Data outside the black line are removed to ensure data quality.

Zhang和Thurber(2003；2006)将双差定位方法(Waldhauser and Ellsworth, 2000)与地震层析成像相结合提出了双差层析成像方法(tomoDD).该方法通过使用绝对走时数据和相对走时数据来联合反演三维速度结构和震源位置参数.根据地震射线理论，台站k记录到的地震i的地震波到时T_kⁱ可以表示为

(1)

其中τⁱ为地震i的发震时刻，u为地震波的慢度.使用截断的泰勒展开对公式(1)进行线性化，得到地震波从地震i到台站k的观测到时与理论到时之差r_kⁱ，公式为

(2)

其中，Δx_lⁱ(l=1, 2, 3)为震源位置参数在三个坐标轴方向上的扰动，Δτⁱ为地震i发震时刻扰动，δu为地震波的慢度扰动.

因此，地震i与地震j到同一个观测台站k的走时残差之差为

(3)

若两个地震距离很近，则它们到同一个台站的射线路径在震源区之外基本上是重合的，使用相对到时数据能够消除震源区外共同射线路径上速度误差造成的影响(Waldhauser and Ellsworth, 2000).在反演中可以根据公式(2)使用绝对到时数据来确定研究区内大尺度的速度结构和绝对震源位置，根据公式(3)使用相对到时数据来确定震源区附近精细的速度结构与相对震源位置.此外，由于震源位置和速度结构存在耦合效应(Thurber，1992)，震源位置的精度在一定程度上会影响速度结构的反演精度，反之亦然.双差层析成像方法对震源位置和速度结构进行联合反演可以使得二者相互促进，得到更高精度的地震定位结果和三维速度结构(邓文泽等，2014).通过与传统方法对比发现，双差层析成像方法得到的定位结果更加精确，获得的速度结构更加精细(Zhang and Thurber, 2003；Okada et al., 2006；于湘伟等，2010).

双差层析成像方法使用三维规则网格进行模型参数化，模型内任意一点的速度可以通过周围节点的速度线性插值获得.为了获得尽可能精细的三维地壳速度结构，我们根据所使用的地震数据和台站的分布情况，在反演之前尝试了多种不同的网格节点划分并对其进行分辨率测试分析.通过比较不同的网格节点划分所得到的结果，最终将反演区域水平方向节点间距设为0.3°×0.3°，垂直方向上节点分别位于0、5、10、15、20、25、35、45、55 km.水平方向上射线分布及网格节点划分如图 3所示.

图 3

Fig. 3

图 3 射线路径及反演网格划分 黑色实线代表射线路径，白色加号代表反演节点. Fig. 3 Ray paths and grid nodes in the study region Black solid lines represent ray paths, white plus signs represent inversion nodes.

本文采用的初始速度模型参考了王椿镛等(1995)的地震测深剖面结果和周民都等(2006)的体波走时反演结果，表 1给出了所使用的一维P波初始速度模型.利用和达法拟合P波与S波震相走时数据，得到研究区内P波与S波的波速比约为1.70.

表 1 Table 1

表 1 门源地区P波初始一维速度模型 Table 1 Initial 1D velocity model of P waves in the Menyuan region 深度/km 0.00 5.00 10.00 15.00 20.00 25.00 35.00 45.00 55.00 速度/km·s-1 4.50 5.60 5.90 6.00 6.20 6.30 6.40 6.60 8.00

表 1 门源地区P波初始一维速度模型 Table 1 Initial 1D velocity model of P waves in the Menyuan region

双差层析成像在反演中引入了阻尼因子和光滑因子来压制数据误差对反演结果的干扰，以保证结果的稳定性.其中，光滑因子用来约束模型慢度参数，阻尼因子则可以同时约束地震位置参数和模型慢度参数(Thurber et al., 2009).我们使用不同的阻尼因子和光滑因子组合构建权衡曲线(Eberhart-Phillips，1986)来选择最佳的参数值，使数据方差较小的同时保证解的变化量不会严重增加.根据权衡曲线(图 4)，最终选定阻尼因子为500、光滑因子为50进行反演计算.

图 4

Fig. 4

图 4 使用不同的阻尼因子(a)和光滑因子(b)得到的L曲线 箭头所指的数字代表选取的最佳参数值. Fig. 4 Tradeoff curves of different damping coefficients (a) and smoothing coefficient (b) The numbers that arrows point to denote optimal coefficients.

为了衡量反演得到的速度结构的可靠性，我们对其进行了棋盘格测试分析(Humphreys and Clayton, 1988；Zhao et al., 1992).首先在初始模型(表 1)的基础上对相邻网格节点添加±5%的速度扰动得到棋盘格模型，然后基于该棋盘格模型和实际震源位置及台站分布通过正演来计算理论走时，使得在棋盘格测试时，地震台站、地震事件和相应的射线分布与实际反演时用到的观测数据完全一致(徐小明等，2015).最后以未添加扰动的初始模型作为参考模型，根据合成的理论数据使用双差层析成像方法进行反演，分析棋盘格的恢复情况.棋盘格测试结果如图 5、图 6所示，在地表附近，受地震和台站分布等因素的影响，分辨率相对较低.在5~25 km深度处，分辨率比较高.而在深度35 km及以下区域，由于地震分布少、Pn和Sn震相数据有限且主要集中在上地幔顶部传播，因此可以看出P波速度结构在35 km、45 km深度处速度异常的恢复能力有所降低，而在55 km深度处(上地幔顶部)因为可以获得较高的射线密度而具有较高的分辨率.对于S波速度结构，由于Sn震相过少(图 2)，此深度区域已无法恢复棋盘格模型正负相间的速度变化.

图 5

Fig. 5

图 5 不同深度P波棋盘格测试结果 Fig. 5 Resolutions of checkboard tests of P waves at different depths

图 6

Fig. 6

图 6 不同深度S波棋盘格测试结果 Fig. 6 Resolutions of checkboard tests of S waves at different depths

经过14次迭代计算，最终得到了12456个地震的重定位结果.重定位前后地震的走时残差均方根分布如图 7所示，其平均值由重定位前的0.281 s下降到了0.033 s.图 8给出了研究区内重定位后地震的震中位置分布，其中门源M_S6.4主震的震中为37.66°N、101.61°E，震源深度为9.65 km，该结果与黄浩等(2017)利用逆时成像技术得到的定位结果基本一致.重定位后地震在东西、南北及深度三个方向上误差的平均值分别为230.21 m、237.08 m和417.32 m.重定位前后地震震源深度变化明显(图 9)，重定位后地震主要集中在2~13 km深度范围内，在垂直剖面上地震呈条带状分布(图 13).

图 7

Fig. 7

图 7 重定位前(a)后(b)地震走时残差均方根分布 Fig. 7 Histograms of travel time residuals RMS before (a) and after (b) relocation

图 8

Fig. 8

图 8 研究区内重定位后地震震中分布图 黑点为重定位后地震震中位置，红色五角星为MS6.4主震位置，红色沙滩球为主震震源机制解，黄色叉形符号为未成功定位的地震.红线为主要断层：F1冷龙岭断裂；F2托莱山断裂；F3肃南—祁连断裂；F4民乐—大马营断裂；F5皇城—双塔断裂；F6武威—天祝断裂；F7古浪断裂；F8金强河断裂；F9龙首山断裂；F10达坂山断裂；F11日月山断裂；F12拉脊山断裂；F13青海南山北缘断裂. Fig. 8 Distribution of relocated events in the study region Black dots represent epicenters of relocated earthquakes, red star is the epicenter of MS6.4 earthquake, red beach ball represent the focal mechanism solution of MS6.4 earthquake, yellow crosses represent the epicenter of earthquakes that not been relocated, and red lines are major faults: F1 Lenglong Ling fault; F2 Tuolai Shan fault; F3 Sunan-Qilian fault; F4 Minle-Damaying fault; F5 Huangcheng-Shuangta fault; F6 Wuwei-Tianzhu fault; F7 Gulang fault; F8 Jinqiang He fault; F9 Longshou Shan fault; F10 Daban Shan fault; F11 Riyue Shan fault; F12 Laji Shan fault; F13 North-edge fault of Nanshan, Qinghai.

图 9

Fig. 9

图 9 重定位前(a)后(b)地震深度分布 Fig. 9 Histograms of depth distribution before(a) and after(b) relocation

图 13

Fig. 13

图 13 BB′剖面上P波和S波速度结构 黑点为发生在剖面两侧各0.1°范围内的地震在该剖面上的垂直投影，红色五角星为门源MS6.4主震位置. Fig. 13 Velocity structure of P and S waves on BB′ profile Black dots denote the vertical projection of earthquakes within 0.1° either side of profile, red star is the epicenter of MS6.4 earthquake.

由Global CMT提供的门源M_S6.4主震的震源机制解为节面Ⅰ：走向335°，倾角47°，滑动角96°；节面Ⅱ：走向146°，倾角43°，滑动角83°.根据重定位结果，余震序列展布的长轴方向与震源机制解的节面Ⅱ较为吻合，且地震现场考察结果也显示震害分布的椭圆长轴走向为NE120°—NE140°(胡朝忠等，2016)，因此我们推断节面Ⅱ为此次地震的震源断层面.

图 10、图 11给出了反演得到的不同深度处P波和S波速度分布，同时将发生在层面上下各2.5 km范围内的地震垂直投影在该层面上.由反演结果可以看出门源地区地壳速度结构存在显著的横向不均匀性.结合研究区内的地质构造情况(图 8)，对比发现0~5 km的浅层速度结构与地表地质构造特征和断层分布具有一定的对应关系.尤其是在分辨率较高的5 km深度处，P波与S波高速区的走向与山脉及主要活动断裂带的走向基本一致，包括祁连—海原断裂、拉脊山断裂、龙首山断裂及其附近区域均为高速异常区.P波与S波低速区则主要分布在西宁、青海湖区域及张掖—民乐一带等，在地形上均表现为盆地或山谷凹陷区.这种特征与前人在其他地区的研究结果相符合(于湘伟等，2010；陈兆辉等，2014；王小娜等；2015)，在一定程度上也佐证了反演结果具有较高的分辨率和可靠性(王长在等，2013).

图 10

Fig. 10

图 10 不同深度水平方向上P波速度分布 黑点表示发生在层面上下各2.5 km范围内的地震在该层面上的垂直投影，红色五角星代表门源MS6.4主震位置，紫色五角星为MS≥7.0的历史地震. Fig. 10 Velocity variations on horizontal slices of P waves at different depths Black dots denote the vertical projection of earthquakes within 2.5 km above and below horizontal slices, red star is the epicenter of MS6.4 earthquake, purple stars are the epicenter of historical earthquakes with MS≥7.0.

图 11

Fig. 11

图 11 不同深度水平方向上S波速度分布 黑点表示发生在层面上下各2.5 km范围内的地震在该层面上的垂直投影，红色五角星代表门源MS6.4主震位置，紫色五角星为MS≥7.0的历史地震. Fig. 11 Velocity variations on horizontal slices of S wave at different depths Black dots denote the vertical projection of earthquakes within 2.5 km above and below horizontal slices, red star is the epicenter of MS6.4 earthquake, purple stars are the epicenter of historical earthquakes with MS≥7.0.

在10 km深度处，门源M_S6.4震中区左侧、青海共和的西南区域和龙首山断裂北部表现为显著的P波速度高值异常.在15 km深度处，其P波速度结构异常分布特征与10 km处相似.在20~55 km深度处，P波速度结构异常分布特征则发生了明显改变，尤其是在35 km深度处，门源地震震中所对应区域出现了显著的低值异常.对于S波速度结构，在10~25 km深度处各层速度结构的异常分布特征总体差异不大，高值速度异常主要分布在门源—古浪—武威一带、青海共和的西南区域和研究区东北部的龙首山断裂附近.根据前面棋盘格测试结果，S波速度结构在35 km及以下区域可靠性较低，因此这部分区域的S波速度结构不予考虑.

研究区内中小地震主要分布在P波和S波高速体内及高低速异常过渡区.门源M_S6.4主震位于高速体与低速体交界处并偏向高速体一侧，其周围表现为明显的P波和S波高速异常.地震波高速异常表明介质具有相对较高的强度、不易发生破裂(吴建平等，2009)，高速坚硬介质的存在有利于高强度应力的积累和集中释放(刘启元等，2009；易桂喜等，2011)，为强震的孕育和发生创造了条件.因此门源地震的发生可能与高速异常体所代表的障碍体有关，这与Pei和Chen(2012)对玉树地震的研究结果一致.此外，小震活动频繁的古浪地震窗同样处在速度高值区域，尤其是S波速度结构更为明显.该区域为1927年古浪M_S8.0地震的震源区，也是青藏高原东北缘小震活动最为显著的区域(冯建刚等，2012).地壳速度结构同样表现为高速异常的龙首山断裂北部地区和青海共和西南地区近几十年来均发生过强震活动，如1954年甘肃山丹M_S7.3地震、1990年青海共和M_S7.0地震等.

为了进一步研究门源M_S6.4地震孕育发生的深部构造背景以及研究区内速度结构的分布与地震活动性的关系，我们经过门源M_S6.4地震震中位置沿着断层走向和垂直断层走向构建了两条速度剖面AA′和BB′(图 8中绿色虚线即为两条剖面的位置)，并将发生在剖面两侧各0.1°范围内的地震垂直投影在该剖面上.

根据AA′剖面(图 12)，可以看出P波速度结构在肃南—祁连断裂东段(100.8°E—101°E)和门源M_S6.4地震震源区存在两个明显的高速异常体，且在M_S6.4地震震源区下方30~40 km处有低速层存在.1982年由中国地震局组织的门源—平凉—渭南人工地震探测剖面其结果同样显示门源地区地壳40 km处存在低速区(王椿镛等，1995).Huang和Zhao(2004)研究发现大地震通常发生在中上地壳的高速体内，并且在震源区的下地壳到上地幔之间存在与流体有关的低速异常体，认为这种构造条件可能导致了中上地壳中发震层的减弱，从而促进了地震成核.此外，在AA′剖面上，S波速度结构也显示门源M_S6.4地震震源区为高速异常区域.

图 12

Fig. 12

图 12 AA′剖面上P波和S波速度结构 黑点为发生在剖面两侧各0.1°范围内的地震在该剖面上的垂直投影，红色五角星为门源MS6.4主震位置. Fig. 12 Velocity structure of P and S waves on AA′ profile Black dots denote the vertical projection of earthquakes within 0.1° either side of profile, red star is the epicenter of MS6.4 earthquake.

图 13给出了BB′剖面上P波和S波的速度分布.在此剖面上，P波高值异常主要分布在37.3°N—37.7°N区域，此高速异常体范围在深度上约从5 km延伸到25 km.中上地壳的S波速度结构(25 km以下区域分辨率较低不考虑)其高速异常区的分布位置与P波基本一致.门源M_S6.4地震及其余震位于此高速体的边缘.余震序列呈近直立分布，与主震的震源机制解倾角存在较大差异，可能与该区复杂的区域构造作用有关(胡朝忠等，2016；黄浩等，2017).剖面右侧龙首山断裂带上(约38.5°N处)地震分布呈南西向倾斜，较好的反映了龙首山断裂带的地下展布特征.

本文使用甘肃、青海数字地震台网及中国地震科学探测台阵记录到的门源地区地震的初至P波和初至S波到时资料，应用双差层析成像方法联合反演得到了该地区的三维精细速度结构和高精度的震源位置参数.结果显示：

(1) 反演之后门源地区地震的走时残差均方根显著降低，其平均值由重定位前的0.281 s降到了0.033 s.地震在东西、南北及深度三个方向上误差的平均值分别为230.21 m、237.08 m和417.32 m.重定位后地震分布在垂直剖面方向上呈现出与断层位置有关的条带状分布.

(2) 门源地区地壳速度结构呈现出明显的横向不均匀性，浅层P波和S波速度结构与地表地形特征及断裂分布密切相关.整体来看，研究区内中上地壳速度结构高值异常主要分布在冷龙岭断裂附近、古浪—武威一带、青海共和的西南区域和龙首山断裂北侧.P波与S波速度结构其高速异常分布总体上具有较好的一致性，说明反演得到的研究区内地壳速度结构具有较高的可信性(谭皓原和王志，2014).

(3) 地震活动性与地壳速度结构具有很强的对应关系.地震主要分布在高速异常体内及高低速异常过渡区.门源M_S6.4地震震中附近的P波、S波速度结构表现出明显的高速异常，深度上从5 km延伸到约25 km，P波速度结构显示在其下方存在一个低速层.M_S6.4地震及其余震位于此高速异常体边缘，高速异常体所代表的坚硬介质容易产生应力积累和集中释放，为强震的孕育和发生创造了条件.

感谢张海江博士为本研究提供TomoDD程序，感谢评审专家的宝贵意见，本文使用了GMT绘图软件，在此深表感谢.感谢中国地震局地球物理研究所“地震科学探测台阵数据中心”为本研究提供地震波形数据.