据中国地震台网中心正式测定，2016年1月21日01时13分13.0秒(北京时间)在青海海北州门源县(37.68°N，101.62°E)发生M_S6.4地震，震源深度10 km (图 1).据青海地震台网快报目录，截至2016年12月31日，共记录到M_L1.8以上余震总数189个，其中最大余震为5.0级.地震发生后，青海省门源县、甘肃省金昌市永昌县、金川区、武威市凉州区、张掖市民乐县和山丹县均震感强烈，但未造成人员伤亡(http://www.cea.gov.cn/publish/dizhenj/468/553/102161/102163/20160121053537836927385/index.html).本次门源M_S6.4地震位于青藏高原东北缘，处于鄂尔多斯地块、秦祁地块、四川盆地与青藏高原各地块的拼接处(图 1).由于印度板块的碰撞挤压(如Molnar and Tapponnier, 1975, 1981; England and Houseman, 1985；Craig et al., 2012)，青藏高原的深部结构极其复杂(Wang et al., 2010; Lei et al., 2014a; Lei and Zhao, 2016a, 2016b; Zhou and Lei, 2016)，历史上曾经发生多次强震，尤其2008年汶川M_S8.0地震以来，在青藏高原东北缘发生了一系列强震，如2013年芦山M_S7.0地震、2014年康定M_S6.3地震以及2013年岷县—漳县M_S6.6地震(图 1).最近研究者开展了青藏高原东缘强震发生与深部结构之间关系的研究，揭示出地壳强震发生可能与地壳结构的横向不均匀性密切相关(Lei and Zhao, 2016b).由图 2a可以看出，距离此次门源地震震中最近的断裂带为北西向的冷龙岭断裂.野外地质考察表明，此断裂走向NW60°~70°、倾向为NE向，倾角为50°~60°，且该断裂由一组平行的次级断裂组成，其断裂带的宽度约1~3 km(何文贵等，2000).

图 1

Fig. 1

图 1 研究区域的历史地震和区域地质构造图 Fig. 1 Map showing historic earthquakes and tectonics in the study area

本次地震发生后，不同的研究机构迅速开展相关研究工作，分别获得了有意义的研究成果.中国地震局地球物理研究所给出的震源机制解结果显示，此次地震为逆冲型地震，矩心深度为7 km (http://www.cea-igp.ac.cn/tpxw/273560.shtml).中国地震局地质研究所给出的地震发震构造背景图显示出此次地震震中附近有2条断裂，分别是冷龙岭断裂和民乐—大马营断裂(http://www.eq-igl.ac.cn/upload/images/2016/1/2191153838.jpg).冷龙岭断裂是全新世的左旋走滑兼逆冲断裂，断裂走向近西北向.民乐—大马营断裂是晚更新世左旋走滑兼逆冲断裂，断裂走向亦为近西北向.然而，中国地震局地球物理研究所的震源机制解显示本次地震主震属于逆冲型地震，而其邻近的两个断层却是以走滑分量为主，其逆冲分量为辅，这种震源机制与野外地质调查结果显示的断层类型的不一致性给理解此次地震的发震机理带来一定困难，因此有必要重新确定本次地震的序列特征和震源机制，以深入探讨发震机理.

本文将利用青海地震台网、甘肃地震台网的固定地震台站及部分流动地震台站的波形数据，通过重新拾取震相到时，联合使用HYPOINVERSE 2000方法(Klein, 2014)和HypoDD方法(Waldhauser and Ellsworth, 2000)对主震及余震序列M_L≥1.8地震的189个事件进行重新定位，采用gCAP方法(Zhu and Ben-Zion, 2013)求解主震震源机制解，深入探讨本次地震的发震机理与构造意义.

本文收集了2016年1月21日至2016年12月31日期间(最后一个M_L≥1.8地震发生在11月8日)在研究区域内(101.3°E—101.9°E，37.4°N—37.9°N)且被青海、甘肃固定地震台站和部分地震流动台站共计37个台站(图 2a)记录到的波形数据.由图 2a可知，由于联合了流动地震台站，使得此次门源地震的台站空间分布具有较大改善，为本研究获得有意义的余震序列定位结果与主震震源机制解提供了有利条件.流动台站数据来自中国地震局地球物理研究所2013—2016年期间沿南北地震带开展的流动地震观测台网记录到的余震资料，该台网的观测设备主要由型号为Reftek-130的数据采集器和型号为CMG-3T的地震计组成.

图 2

Fig. 2

图 2 (a)本研究所使用的初始地震震中和地震台站分布图.标注有代号的台站为拾取波形时所增加的甘肃台站；(b)时距曲线图；(c)本研究定位采用的速度模型 Fig. 2 (a) Initial earthquake hypocenters and seismic stations used in this study. Stations with codes are added for picking up waveforms from Gansu seismic network; (b) Time-distance curves; (c) Velocity models used for location in this study

本研究拾取震相所用波形数据是以青海台网快报目录为依据.由于较小震级的地震波形信噪比较低，再加上波形上各种反射、折射震相相互叠加，导致到时难以精确提取，因此本研究由波形上提取到时所选取的震级下限为M_L1.8，地震事件总共189个，其中4.0级以上地震有3个、3.0~3.9级之间的地震事件有18个、2.0~2.9级之间的地震事件有133个、1.8~1.9级之间的地震事件有35个.为了保证拾取初至震相是直达波，截取波形时按照青海台网快报目录位置只选取了震中距为200 km以内的地震台站(图 2a).本研究参与定位的地震台站共31个(图 2a).震源机制反演也采用重新定位时所用的宽频带波形资料.在反演主震震源机制时，选择震中距在260 km范围内的台站.根据波形质量和台站方位角分布，筛选出14个台站(图 2a)的波形资料用于开展震源机制研究.需要说明的是，尽管青海省地震局已经给出比较完整的震相到时资料，但是由于此次地震发生在青海和甘肃交界地区，而青海省局给出的震相报告中所含甘肃台站较少，为了充分利用足够多资料进行定位，因此我们截取波形重新拾取初至波到时资料(图 3).

图 3

Fig. 3

图 3 青海门源MS6.4地震主震到时拾取的波形 蓝色台站为增加的甘肃台网台站. Fig. 3 Seismic waveforms for picking the arrival times of the MS6.4 Menyuan, Qinghai mainshock Blue stations are added seismic stations of the Gansu network.

为尽量保证到时数据的精度，在挑选Pg和Sg时采用以下标准：(1) 每个地震事件最少有3条地震台站记录；(2) 每个台站最少记录5个地震事件；(3) 在震相选取时只选择波形上清晰的到时(图 3).以上所述这些标准作为条件，经拾取震相后，我们最终获得2259条Pg震相到时资料和1156条Sg震相到时资料(图 2b).为证实本研究拾取的到时资料的可靠性，我们将本研究震相拾取结果与青海台网和甘肃台网观测报告中的震相到时分别进行了对比，结果显示我们拾取的P波到时与青海台网和甘肃台网报告到时资料差在0.1 s内的资料数可达83.8%和83.3%, 我们拾取的S波到时与青海台网和甘肃台网报告到时资料差在0.1 s内的资料数可达63.9%和60.0%.这些结果说明，观测报告中的到时资料与我们提取的到时资料非常接近，处于误差范围内，由此表明本研究所拾取震相较为可靠.

高精度的定位能够更加准确揭示地震序列的空间分布特征，探讨与断裂构造活动的关系和构造的成因(Lei et al., 2012; 张广伟等，2014；王未来等, 2014a, 2014b；房立华等，2013).本文首先采用HYPOINVERSE 2000定位方法对地震序列进行了绝对定位，该方法的基本原理是基于传统的Geiger法基本思路，采用线性化反演方法，利用最小二乘法求解震源坐标和发震时刻的校正量，在求解过程中采用多种数据类型加权.本研究在进行绝对定位时，初始发震时刻设为首个台站记录到时的前2 s，震中位于第一个到时台站的下方，初始震源深度设为10 km；通过绝对定位获取较为精确的震源位置，尤其在改善震源深度方面.然后，使用HypoDD定位方法对余震序列进行相对定位，该方法则是使用相对走时残差来修订地震位置.若两个地震之间的距离远小于地震事件到台站的距离，那么这两个地震传播到台站的路径近似相同，因此这两个地震的到时差可以归因为地震位置的不同.在进行相对定位时，需要对P波和S波确定先验权重.由于S波震相不容易识别，使得其读取精度要小于P波震相读取精度.因此，在本次研究中将P波走时赋予先验权重为1.0，而S波赋予的先验权重为0.5.为保证定位精确性，取地震对间距最大距离为4 km，最小连接数为8.采用共轭梯度最小二乘法求解，共迭代16次，在前8次迭代中采用标准偏差的10倍值作为截断值，后8次迭代中采用标准偏差的6倍值作为截断值.通过反复迭代，除去残差大于截断值的震相数据.每次迭代利用上次迭代的结果更新震源位置、残差和矩阵，将每次迭代的残差大小作为下一次迭代的加权函数，经过多次迭代，直到得到稳定解.通过相对定位方法获取精确的地震间相对位置，为能够较好地描述发震断层几何结构特征提供可靠的地震学证据.

震源机制解的求解采用gCAP法.该方法是在CAP方法(Zhu and Helmberger, 1996；Xie et al., 2013；Zheng et al., 2009)基础上，将波形反演中包含非双力偶分量，其基本原理是将宽频带数字地震波形记录分为体波Pnl和S波(或面波)两部分并分别进行带通滤波；通过对Pnl和S波(或面波)两部分赋予不同权重，计算理论和实际波形的拟合误差函数；最后采用网格搜索得到最小误差的最优解(谢祖军等，2012; 张广伟和雷建设，2015).为更好地探讨震源破裂机制，本研究分别求取主震的双力偶机制解和全矩张量解.反演中，体波Pnl和S波截取波形窗长分别设置为30 s与60 s，Pnl和S波滤波范围分别为0.05~0.2 Hz和0.02~0.1 Hz，走向、倾角和滑动角的搜索间隔均为5°，深度搜索间隔为1 km.格林函数采用频率波数(FK)法来计算(Zhu and Rivera, 2002)，其中采样间隔设为0.1 s，采样点个数为1024.

本文在地震序列重定位和主震震源机制反演过程中均采用如图 2c所示的一维层状速度模型.该初始一维速度模型主要参考了王椿镛等(1995)的人工勘探和姚志祥等(2014)接收函数的研究结果，并将泊松比设为1.73.

采用HYPOINVERSE 2000定位方法对本研究拾取的震相对门源主震及余震序列共189个地震事件进行了绝对定位，结果显示定位后的主震位置为37.68°N、101.60°E，震源深度为13.12 km.主震水平方向和深度方向的定位误差分别为0.9 km和1.9 km，余震序列的水平方向和深度方向的定位误差分别为1.1 km和4.6 km.然后，将绝对定位结果作为门源地震序列的初始位置，再用HypoDD相对定位法进行重新定位，期望获得更为精细的定位结果，最终获得163个地震事件的重新定位结果.本次门源地震主震震中位置修订为37.67°E、101.61°N，震源深度修订为11.98 km.重新定位后，沿纬度方向、经度方向和深度方向的定位平均误差分别为0.42 km、0.42 km和0.86 km.经过双差定位后，误差比绝对定位明显减小，在深度方向上尤其如此.

由于门源地震所在区域的波速横向不均匀性非常明显，因此为检验速度模型的不确定性对双差定位结果的影响程度，本研究对图 2c所示的速度模型分别施加了-10%和10%的速度扰动.通过双差定位后，两次试验均获得163个定位事件结果，与我们重新定位后获得的重定位个数一致.由图 4和表 1可以清楚看出，模型速度减小10%后重定位结果在水平方向上基本向地震簇中心移动，且最大移动量约1.0 km的事件基本处于地震簇边缘，深度方向最大移动量约1.67 km；而模型增加10%后重定位结果在水平方向上基本向地震簇外移动，且最大移动量约2.0 km的事件也处于事件簇边缘，深度方向最大移动量可达5.5 km.尽管如此，整体位置变化量较小(图 4). 表 1给出了这2次试验沿纬度方向、经度方向和深度方向的定位结果变化量，其平均值分别为0.42 km、0.42 km和0.18 km及0.42 km、0.43km和0.38 km，这一空间位置的变化量完全处于前述误差范围内，说明双差定位方法不会由于速度模型的不确定性对定位结果产生有意义的影响(赵博等，2013).

图 4

Fig. 4

图 4 在图 2c速度模型上分别增加(a—b)和减少(c—d)10%速度变化情况下的定位结果(箭头终点)和本研究重新定位结果(圆圈)的变化图 Fig. 4 Variations of location results (end points of arrows) inferred from the models as shown in Fig. 2c with the velocities increased (a—b) and decreased (c—d) by 10% relative to those (circles) inferred from unchanged models

表 1(Table 1)

表 1 不同速度扰动情况下的定位结果与本研究重新定位结果的变化分析 Table 1 Variations of location results using different velocity models

表 1 不同速度扰动情况下的定位结果与本研究重新定位结果的变化分析 Table 1 Variations of location results using different velocity models

为进一步清晰了解重定位的效果，我们将相对定位前后结果进行对比.重定位前(图 5a)，地震震中在平面图上相对零散分布，在深度方向上主要集中在10 km；重定位后(图 5b)，地震震中在平面图上沿NW—SE向的优势展布更为明显，在深度方向上显示了明显的4~14 km发震层(图 5b).从重定位深度上看，门源M_S6.4主震震源深度由初始的中国地震台网中心正式观测报告10 km改进13.12 km，再经相对定位后变为11.98 km.

图 5

Fig. 5

图 5 重定位前(a)和重定位后(b)地震平面及经纬向侧面图 大星号代表门源主震，小星号代表最大ML 5.0余震，圆圈代表余震，圆圈大小代表余震震级，F表示频度(占总事件的比例). Fig. 5 Earthquake distributions in planar view and along the longitude and latitude directions before (a) and after (b) relocations The big star denotes the Menyuan mainshock. The small star denotes the largest aftershock ML 5.0. Circles are aftershocks, the size of which represents magnitude. F stands for earthquake frequency(the proportion of the total number of events).

采用gCAP方法通过不同深度最佳双力偶机制解和地震全矩张量反演，得到主震的最佳质心深度分别为9 km和8 km(图 6)，说明此次地震最佳质心深度在8~9 km之间，与双差定位结果11.98 km有一定差别.这样的差别可能由于其反映的震源深度概念不同，双差定位采用初至震相到时获得的震源深度反映了震源初始破裂点，而震源机制采用资料包含有面波因而获得的深度反映的是破裂面的质心深度. 图 7展示出9 km深度上主震最佳双力偶机制解(下半球投影)以及理论和实际波形拟合图，可以看出14个台中所用波形段总共70个，其中相关系数在70%以上有63个，占总数的93%，反演方差为6.216. 图 8展示出本次地震全矩张量解(下半球投影)以及理论和实际波形拟合图，也可以看出14个台中所用波形段总共70个，其中相关系数在70%以上有63个，占总数的93%，反演方差却降至5.965.两种反演结果整体上波形拟合效果较好、所获得的震源参数均较为可信，且全矩张量解的残差较小，说明全矩张量解更能客观地反映震源机制.在地震全矩张量反演时得到各向同性(Isotropic, ISO)、补偿线性偶极子(Compensated Linear Vector Dipole, CLVD)和双力偶(Double Couple, DC)三部分所占的百分比分别为5.15%、0.723%和94.61%，以及ζ(ISO部分在全矩张量M_ij中所占的大小)和χ (CLVD部分在偏张量D_ij中所占的大小)具体值分别为0.23和-0.09.这些结果说明，门源M_S6.4地震非双力偶成分较小，主要是双力偶破裂机制.地震全矩张量解中包含有一定的非双力偶成分，表明震源破裂过程较为复杂，这种复杂性可能与本次地震所处构造环境有一定关系，比如震源区介质的各向异性或断层面的非平面结构(Vavryčuk, 2005).另外，通过与表 2中不同研究机构的震源机制结果对比，发现本研究两次结果与其他机构的结果基本一致，均显示出此次门源主震为逆冲型地震.

图 6

Fig. 6

图 6 2016年1月21日MS6.4主震不同深度上双力偶机制解及其反演方差(a)和矩张量解及其反演方差(b) Fig. 6 Double couple mechanism (a) and moment tensor solutions (b) at different depths versus inversion variance of the 21 January 2016 Menyuan MS6.4 earthquake

图 7

Fig. 7

图 7 2016年1月21日MS 6.4主震双力偶机制解及理论(红色)和实际(黑色)波形对比 波形下方数字表示理论相对实际波形的移动时间和二者的相关系数，左侧大写字母表示台站名，台站下方数字为震中距和相对偏移时间. Fig. 7 Double-coupled source mechanism solution and comparison between the theoretical (red) and observed (black) waveforms of the 21 January 2016 MS6.4 mainshock The numbers below the waveforms denote the time movement of theoretical relative to observed seismograms and their correlation coefficients. The capital letters on the left denote seismic stations. The numbers below the stations are epicentral distances and relative time-movement.

图 8

Fig. 8

图 8 2016年1月21日MS 6.4主震矩张量解及理论(红色)和实际(黑色)波形对比 波形下方数字表示理论相对实际波形的移动时间和二者的相关系数，左侧大写字母表示台站名，台站下方数字为震中距和相对偏移时间. Fig. 8 The moment tensor solution and comparison between the theoretical (red) and observed (black) waveforms of the 21 January 2016 MS6.4 mainshock The numbers below the waveforms denote the time movement of theoretical relative to observed seismograms and their correlation coefficients. The capital letters on the left denote seismic stations. The numbers below the stations are epicentral distances and relative time-movement.

表 2(Table 2)

表 2 不同研究机构震源机制结果 Table 2 Focal mechanism solutions from different organizations

表 2 不同研究机构震源机制结果 Table 2 Focal mechanism solutions from different organizations

一般地讲，余震分布大致能够直观地展示断层面的破裂延展范围(万永革等，2008；张广伟和雷建设，2013)，因此本研究利用相对定位结果分别展示出震后24小时和震后292天(截至2016年11月8日)时间范围内分别沿余震展布方向(纵剖面AA′)及垂直方向剖面(纵剖面BB′)(图 9).由震后24小时剖面图(图 9a)可以看出，纵剖面AA′显示此次地震余震在时间尺度上主要沿着东南向展布，说明在此时间段此次地震破裂主要由西北向东南向扩展，约12 km；由纵剖面BB′可以看出倾角约45°或更陡的发震断层面.结合胡朝忠等(2016)现场考察结果，认为发震断层节面Ⅱ，其走向约140°左右(表 2).

图 9

Fig. 9

图 9 主震后24小时(a)和震后292天(b)重定位后地震事件平面及纵剖面图 大星号代表门源主震，小星号代表最大ML 5.0余震，圆圈代表余震，圆圈大小代表余震震级，颜色代表主震后的时间，黑色实线代表可能的发震断层. AA′为沿断裂带方向纵剖面，BB′为横穿主震且垂直于断层的纵剖面.剖面纵横向比例1:1.震源机制解结果取自本研究gCAP双力偶机制解. Fig. 9 Relocated earthquakes in planar view and along the cross-sections in 24 hours (a) and 292 days (b) of the mainshock The big star denotes the Menyuan mainshock. The small star denotes the largest aftershock of ML 5.0. The circles are aftershocks, the size and color of which represent magnitude and time after the mainshock. The black line denotes possible seismogenic fault. AA′ is the cross-section along the fault. BB′ is perpendicular to the fault going through the mainshock. The scale of sections is 1:1. The source mechanism solution is from the gCAP method.

8月13日在门源主震NW方向发生M_L5.0较大余震，为更好地探讨门源主震的破裂范围及余震序列优势发震层，我们给出震后292天余震序列分布的剖面图(图 9b)，结果显示沿着断裂方向AA′的地震延展范围扩展到16 km，与震后24小时时间范围内的余震空间扩展范围有一定差别.尽管沿着断裂AA′东南端余震延展范围变化较小，但随着M_L5.0余震的发生，余震更多的向主震西北向扩展，这与四川芦山地震序列较长时间单侧破裂(Lei et al., 2014b)不同. 1986年8月26日门源M_S6.4地震破裂尺度为20 km(汪进等，1992)，此次地震在震后24小时内破裂尺度仅为10 km，表明此次地震在震后24小时时间段内地震能量并没有完全释放，随着较大余震的发生，地震能量才逐渐释放.余震序列优势发震层在4~14 km深度范围内，浅部地震少，且区域人口分布稀疏，这也许是此次地震震害轻且无人员伤亡的重要原因.

本文将HYPOINVERSE 2000绝对定位方法与HypoDD相对定位方法相结合，获得2016年1月21日青海门源M_S6.4地震以及震后292天内M_L≥1.8共189个地震序列的重定位结果，采用gCAP方法获得了主震双力偶机制解和全矩张量解.定位结果显示，余震优势展布方向为NW—SE，长度约16 km，震源优势分布深度在4~14 km之间，且余震没有浅层地震.深度剖面展示出主震的发震破裂面较陡，结合gCAP震源机制解结果，认为此次地震主震为逆冲型地震.定位结果和震源机制结果显示出与震中区域最近的冷龙岭断裂活动性质不太一致，因此推断本次门源M_S6.4地震序列的发震断层极有可能是位于北西向的冷龙岭断裂和民乐—大马营断裂东南段之间的一条盲断层.

青海门源M_S6.4地震发生在祁连山断裂带东段，属于青藏高原东北缘，长期受到印度板块东北向挤压作用(Lei et al., 2014a; Lei and Zhao, 2016a, 2016b)，使应力不断积累.尤其是2008年汶川M_S8.0地震发生后，引发了南北地震带不同断裂带的应力调整，使得青藏高原东北缘地区地震活动性增强，如巴颜喀拉地块南边界的2010年4月14日青海玉树M_S7.1地震(刘巧霞等，2012)、青藏高原东边界的2013年4月20日芦山M_S7.0地震(张广伟和雷建设，2013；Lei et al., 2014b)和2014年11月22日康定M_S7.3地震(易桂喜等，2015)以及2013年岷县—漳县M_S6.6地震(王未来等，2014b).这些地震印证了青藏高原东北缘是地震活动强烈的地区(周民都等，2006；陈九辉等，2005；Lei and Zhao, 2016b)，均为受青藏高原东北向挤压，造成断层应力失稳而形成的地震.此次地震与1986年门源M_S6.5地震位置(图 1)非常接近，可能也位于冷龙岭断裂和民乐—大马营断裂之间的这条盲断层上，表明该断裂已表现出新的活动迹象.

感谢青海省地震局马建新高级工程师在区域台网资料方面提供的支持；感谢中国地震局地壳应力研究所张广伟在本项研究工作过程中给予的帮助；同时，中国地震台网中心地震台网部的全体同事在工作中予以帮助，在此一并表示诚挚谢意.感谢匿名审稿人对本文提出宝贵的修改意见.