震源机制是研究地震构造机理的理想方法，实际震源破裂过程含有断层面的剪切错动、震源体的膨胀或收缩、优势方向的张裂或挤压变形等，因此中小地震、大地震事件的震源机制可用地震矩张量较为准确地描述(Chapman and Leaney, 2012).近年来，矩张量用于震源机制的研究得到广泛应用，如构造地震震源机制的研究(赵翠萍等，2008；郑勇等，2009；李圣强等，2013；易桂喜等，2015)、火山及地热区域地震震源机制的研究(Julian and Sipkin, 1985；Foulger，1988；Foulger et al., 2004；Minson et al., 2007；张广伟和雷建设，2015)和核试验与采矿区诱发地震震源机制的研究(Dreger and Woods, 2002；Walter et al., 2009；Baig and Urbancic, 2010)等.

2013年吉林前郭震群位于太平洋板块俯冲带弧后构造地区(Lei and Zhao, 2005，2006；Lei et al., 2013)的松辽盆地中部北东向扶余—肇东断裂和北西向查干花—道字井断裂的交汇部位附近(图 1)，同时震源区也位于油田开采区内，并且附近分布有多口采油井.该震群具有强度大、频率高、分布范围集中、爆发性强的活动特点.其发震机制是与太平洋板块俯冲引起弧后扩张有关，还是与油田开采有关，目前尚不十分明确.全矩张量解中双力偶分量和矩心深度是区分这两种成因最重要的参数(Cesca et al., 2014).全球矩心矩张量中心的GCMT与吴微微等(2014)的矩张量反演结果均显示与普通的构造地震有很大的差异，主要表现有两个方面：(1) 矩张量解中非常低的双力偶分量，(2) 很浅的矩心深度.然而，GCMT、吴微微等(2014)只给出了前郭震群中震级较大地震的偏矩张量解，其结果对理解该震群发震机理具有一定局限性.由于地震矩张量的研究主要包含有矩张量反演、矩张量分解、参数稳定性评估及其解的图形化(Chapman and Leaney, 2012)，因此本研究采用矩张量反演方法ISOLA(Zahradník et al., 2005)对前郭震群中M_S≥ 5.0地震开展多种模式的矩张量反演研究，并且对所得的偏矩张量解、全矩张量解和纯双力偶解利用F-test进行显著性检验，使用Jackknifing(刀切法)讨论最优解矩张量各分量的稳定性，最后将矩张量解投影到Hudson震源类型图上(Hudson et al., 1989)，深入探讨前郭震群的发震机理，对于防灾减灾具有重要科学意义.

图 1

Fig. 1

图 1 前郭震群地震信息及研究区历史地震活动性 (a)深色区域为松辽盆地，黑色三角为宽频带地震台站，白色三角为短周期地震台站，蓝色圆斑为历史地震，红色圆斑为2013年前郭震群，震级与圆斑半径成正比，F-5-3是第二松花江断裂；(b)是对(a)图红色框区域进行放大图，F-5-1是扶余—肇东断裂，F-5-2/2′是查干花—道字井断裂；(c) HypoDD方法对震群序列重新定位后的震源深度分布直方图(王亮等，2015)；(d)和(e)地震序列的M-T图，时间段分别为2013年10月31日—2013年12月30日和2014年1月1日—2016年4月30日. Fig. 1 The information of Qianguo seismic swarm and the seismicity of historical earthquakes in the study area (a) The largest dark area represents the Songliao Basin, whereas the black triangles and white triangles show the location of broad-band and short-period seismic stations, respectively.The blue circles denote historical earthquakes, where as the red ones denote the 2013 Qianguo seismic swarm and their radius is proportional to the magnitude.F-5-3: the Second Songhuajiang River Fault; (b) The enlarged area shown in(a).F-5-1, Fuyu-Zhaodong Fault, F-5-2/2′, Chaganhua-Daozijing Fault; (c) Histogram for the source depth distribution inferred from the HypoDD method(Wang et al., 2015); (d—e) The M-T map of the seismic swarm sequence from December 30, 2013 to October 31, 2013, and April 30, 2016 to January 1, 2014.

自2013年10月31日开始，吉林省前郭有M_S5.8震群活动以来，在一个月的时间内发生M_S≥5.0地震5次，最大地震为2013年11月23日的M_S5.8地震.震群位于松辽盆地中部北东向扶余—肇东断裂和北西向查干花—道字井断裂的交汇部位附近(图 1a，b)，该区域100 km范围内，1119年曾发生前郭6³/₄地震(唐雅芝，1990；吴戈，1991)(图 1a).本次震群的同一构造部位曾发生2006年吉林省乾安M_S5.0地震(图 1a).到目前为止，M_S5.8震群的余震活动依然没有结束，截至2016年4月30日，共记录余震1351次，其中M_S5.0~5.9地震5次，M_L 4.0~4.9地震15次，M_L 4.0以上地震呈连发特征并伴有大量M_L 1.0以下余震活动(图 1d，1e)，图中标出5级地震发生时刻为：北京时间2013年10月31日11时03分35秒发生首震M_S5.5、紧接着于10分06秒发生M_S5.0地震、11月22日16时18分51秒发生M_S5.3地震、11月23日06时04分23秒发生M_S5.8地震，随后于32分29秒发生M_S5.0地震.序列中大部分地震的波形资料中含有很强的面波成分，如图 2所示的是震中距为142 km、方位角为132°的长春台(CN2) 记录到的5次地震的垂直向波形，可以看出4个地震的面波均非常发育，可见本次前郭Ms 5.8级地震群的震源较浅.

图 2

Fig. 2

图 2 CN2台站的垂直分向原始波形 Fig. 2 The Z-component original seismograms of five events at CN2 station

前郭震群发生后，吉林省地震局迅速布设了7个流动地震观测点，加上吉林省、黑龙江省、内蒙古和辽宁省台网原有的45个固定台站，共计有52个地震台站参与定位.这些台站对震中区域形成很好的方位包围，且记录到了完整的地震序列.我们选择序列中5个M_S≥ 5.0地震作为研究对象(图 1d).为了避免地壳结构横向不均匀性造成的影响，参与矩张量反演的台站选择范围限制在松辽盆地内部，共计使用10个宽频带台站(图 1a)，震中距范围约为86~255 km.

为了考察地壳速度模型对结果的影响，在研究中使用了3个模型VM-1、VM-2和VM-3(表 1和图 3).VM-1是以CRUST2.0导出模型为基础，参考近年来松辽盆地及东北地区地球物理研究成果，如中国满洲里—绥芬河地学断面项目(杨宝俊等，1996)、中朝地台东北缘地区的地震层析成像结果(卢造勋等，2002)、松辽盆地浅部基底推覆伸展作用的地震学证据(单玄龙等，2009)等，对研究区莫霍面深度及中间层速度不连续面进行了细微调整，获得更符合松辽盆地地壳结构的速度模型.VM-2是吴微微等(2014)用到的速度模型，该模型参考了背景噪声互相关、瑞利波频散曲线以及近远震联合层析成象方法得到的东北地区岩石圈速度结构资料.VM-3是根据接收函数研究成果获得的前郭震中区速度模型.

表 1 Table 1

表 1 地壳速度模型信息 Table 1 The information on crustal velocity models VM Top depth of each layer (km) VP (km·s-1) VS (km·s-1) VM-1 0.0 2.50 1.20 1.0 4.00 2.10 3.5 6.20 3.60 15.5 6.40 3.70 33.0 8.20 4.70 VM-2 0.0 2.5 1.45 1.0 4.00 2.32 3.5 6.20 3.60 15.5 6.40 3.74 30.0 8.00 4.61 VM-3 0.0 4.0 2.10 2.0 6.2 3.60 14.0 6.4 3.60 26.5 6.8 3.80 39.0 8.2 4.70

表 1 地壳速度模型信息 Table 1 The information on crustal velocity models

图 3

Fig. 3

图 3 研究用的速度模型 Fig. 3 Crustal velocity models used in the study

地震矩张量反演需要准确的地震定位、可靠的地壳速度模型、多个方位角覆盖好的台站及其高信噪比的地震波形资料(Henry et al., 2002；Ford et al., 2010；Davi and Vavryčuk 2012).为获得较可靠的矩张量结果，本研究首先对5个M_S≥ 5.0地震进行重新定位，定位速度模型使用模型VM-1.然后，通过信噪比选择滤波频段，利用不同的台站组合制定不同的反演方案，进行偏矩张量反演，利用反演结果评价指标确定最佳的反演方案，分析矩张量分解的双力偶分量，最后开展全矩张量反演及纯双力偶反演.

本研究中地震重新定位采用地壳速度模型VM-1(表 1和图 3)和HYPOINVERSE-2000(Klein，2002)定位方法.为了提高定位精度，选取距前郭震中约500 km范围内的地震台站，且采用远、近台站相结合的策略(张天中等，2007)对5个M_S≥ 5.0地震进行重新定位.定位后，震源的水平方向和垂直方向的定位标准差分别为0.4 km和0.8 km，发震时刻标准差为0.25秒.重定位结果(表 2)显示，震源深度范围为5.2~6.3 km.可以看出，相对于多数构造地震事件，这些地震事件均属于较浅的地震.

表 2 Table 2

表 2 5个MS≥ 5.0地震的重新定位结果 Table 2 Relocation of the five MS≥5.0 events No. MS Date Time Lat.(°N)， Lon.(°E) Depth (km) 1 5.5 2013-10-31 03:03:34 44.66，124.12 6.2 2 5.0 2013-10-31 03:10:07 44.66，124.12 5.2 3 5.3 2013-11-22 08:18:49 44.66，124.14 6.0 4 5.8 2013-11-22 22:04:23 44.64，124.12 5.4 5 5.0 2013-11-22 22:32:29 44.63，124.13 6.3

表 2 5个MS≥ 5.0地震的重新定位结果 Table 2 Relocation of the five MS≥5.0 events

矩张量反演需要对观测波形进行带通滤波，在仪器响应的平展段内，选择2个合适的频率作为滤波频率的上下限进行滤波.为了减小地壳速度模型不准确性对研究结果的影响，滤波的下限频率需要在满足反演信噪比要求的基础上尽可能地低，而上限频率的确定是根据地震的震级与震中距选择相应的经验值.对于5级地震，震中距在100~300 km，滤波的上限频率选择通常为0.05或者0.06 Hz(Zahradník et al., 2005；Benetatos et al., 2013)，因此本研究着重利用谱分析法对滤波下限频率进行研究.本研究通过对每个台站三通道波形初动到达后的P波段和S波段开展快速傅立叶变换获得的信号功率谱来计算波形信噪比，将信噪比SNR=10作为阈值(Zahradník et al., 2005).结果显示，大部分台站的波形频率下限高于0.02 Hz时，其信噪比就超过阈值，而个别台站(如SGT台、LUB台)的频率下限必须高于0.04 Hz.我们最终确定了每个台站的滤波频率上下限(表 3中的方案1).图 4是以M_S5.8地震的CN2台站为例，展示三分向波形的谱分析结果，可以看出图 4d显示出CN2台站的两个水平分向波形的信噪比均较高，但垂直分向(Z-component)信噪比较低，因此在反演时就要剔除了这种信噪比低的波形资料，否则会对结果的可靠性产生一定的影响.

表 3 Table 3

表 3 矩张量反演方案及可靠性分析 Table 3 The moment tensor inversion cases and stability analysis Case Station freqlimits MW Depth (km) DC (%) VR (%) SNR CN FMVR STVAR K-angle (°) 1 SGT:0.04~0.08|CN2:0.02~0.06*a BCT:0.03~0.06|WLT:0.02~0.06*a SPT|AGL|YST| WUT:0.03~0.06*b IDR| LUB:0.04~0.06*b 5.1 4.0 44.6 78.0 163 5.8 23 0.28 27 2 SPT|SGT|CN2|YST|WUT*c 5.1 3.0 43.7 92.0 456 6.5 17 0.18 20 3 IDR |BCT|LUB|AGL |CN2 |WUT*c 5.1 4.0 44.7 92.0 515 3.3 15 0.10 21 4 BCT|AGL|WUT*c 5.1 4.0 41.9 92.0 515 3.2 27 0.20 25 5 IDR|BCT|LUB|AGL|SPT|SGT|WUT*c 5.1 4.0 87.3 92.0 431 5.3 14 0.15 20 注：*a，台站之间用“|”隔开，数值表示滤波频率的上下限；*b，多个站使用共同的滤波频段；*c，所有台站采用方案1的滤波频率.

表 3 矩张量反演方案及可靠性分析 Table 3 The moment tensor inversion cases and stability analysis

图 4

Fig. 4

图 4 宽频带台站波形信噪比分析 (a) CN2台站的三通道波形信号功率谱.红色曲线为P/S波位移功率谱，黑色曲线为噪声功率谱；(b) CN2台站的三通道波形信噪比.红色曲线为位移信噪比，蓝色曲线为平均信噪比.黑色虚线表示滤波下限频率. Fig. 4 The SNR of broadband stations seismic waveform (a) The signal power spectrum of three channel waveforms recorded at seismic station CN2, where the red curves represent the P/S wave displacement power spectrum, and the black curves represent the noise power spectrum; (b) The signal to noise ratios of waveforms recorded at seismic station CN2.The red curves are the displacement signal to noise ratios, while the blue curves are the average signal to noise ratios.The black dashed line indicates the low frequency threshold value.

地震矩张量是对多种震源类型的统一描述，用一个二阶张量表示，有多种分解方法，本研究将矩张量分解为ISO(Volumetric)、CLVD(Compensated linear vector dipole)和DC(Double couple)三个部分(Vavryčuk(2001)).ISO为震源体积变化分量，CLVD为补偿线性矢量偶极分量，DC为纯双力偶分量.通常将只包含CLVD与DC两个分量的解称为偏矩张量解(DEV)，包含ISO、CLVD和DC三个分量的解称为全矩张量解.我们首先对5个M_S≥ 5.0地震进行偏矩张量反演，通过5个评估指标选择最佳的台站组合方案，反演结果与GCMT结果进行比较，同时考察结果中双力偶(DC)分量的大小，然后采用最佳方案进行纯双力偶及全矩张量反演.

由于ISOLA方法(Zahradník et al., 2005)是在时间域将迭代反卷积多点源反演方法(Kikuchi and Kanamori, 1991)发展到全波形矩张量反演并采用网格搜索及最小二乘法获得地震的最佳矩张量解，因此该方法具有纯DC、偏矩张量及全矩张量反演三种方式的优势.另外，该方法既可将震源看作单点源也可将震源区视作多点源来处理，且在4维空间中(平面，深度，时间)搜索矩心位置.本研究中，我们采用VM-1(表 1和图 3)使用频率波数法(Bouchon，1981)计算全波场格林函数，在此过程中将震中固定，矩心深度限制在0~15 km范围内，每隔1 km计算一个格林函数，其矩心时间在发震时刻前后2 s范围以0.2 s作为步长来调整.

因为参与反演的台站数量、方位角覆盖范围、台站间方位角间隔及台站几何分布均会对反演结果造成影响(Davi and Vavryčuk, 2012)，所以我们依据VR、CN、SNR、FMVAR、STVAR 5个参数(Sokos and Zahradník，2013)选取最佳反演方案.根据参与反演的台站数量与方位角间隔，制定了5种台站组合方案(图 5).由于篇幅的原因，本研究仅以M_S5.8地震为例进行介绍(表 3).

图 5

Fig. 5

图 5 5种组合方案所用台站示意图 Fig. 5 The distribution of seismic stations used in the five cases

方案1：使用尽可能多的台站.地震台网在盆地内共计10个宽频带台站，震中距最小86 km(SGT台)，最大255 km(LUB台).

方案2：选择台站分布较多的东南方向.由于台站密度分布不均匀，在方位角80°~150°之间有很密集的台站分布，在这个区间选择5个台站.

方案3：选择最大的方位角覆盖范围.最小的方位角间隔，但台站数量尽可能的少，最终选择6个台站参与反演.

方案4：选择最少的台站数量.最大的方位角覆盖范围，最终在三个方位选择3个站.

方案5：远近台站结合.在以上最佳方案中增加一个近台(SGT台).

VR(Variance reduction)为约化方差，用来表示观测波形和合成波形的拟合程度，用下式计算：

(1)

式中，o、s分别表示观测波形和合成波形.上述5个方案的VR值均超过90%，表示观测波形和合成波形有很好的拟合，同时也说明地壳速度模型VM-1能很好地反映实际地壳速度结构.CN(Condition Number)是一个相对综合的指标，来自格林函数矩阵，是最大特征值和最小特征值的比值，值越小，表示矩张量解的稳定性越高，计算公式如下：

(2)

(3)

在(2)、(3) 式中，d表示观测波形，G为格林函数，m表示模型参数.max_eigenval和min_eigenval表示矩阵G的最大与最小特征值.SNR为参与反演的台站所有波形的平均信噪比.FMVAR(Focal Mechanism Variability Index)表示解的离散程度.在反演过程中，不同的震源深度、不同的平面位置、不同的发震时刻会获得多个结果，FMVAR用各个震源机制解与最佳结果的空间夹角K-angle角(Kagan，1991)的方差表示.K-angle表示任意空间两个震源球之间最小夹角，在20°范围内认为是相同的震源机制.STVAR(Space Time Variability Index)是时空比值，即反演中发震时刻的允许范围与平面、垂直震中位置的允许范围组成的4维时空向量的点积(相当于4维空间的体积)与理论上全部解的时空向量点积的比值.在实际工作中往往与FMVAR联合使用，STVAR越小，表明解的可靠性越高.

为了检查各方案解的稳定性，把GCMT发布的M_S5.8地震矩张量结果作为参考，求各个方案解的K-angle角，由表 3可以看出K-angle在20°~27°范围内，变化范围较小，约在7°以内，表明5个方案中双力偶部分的解比较稳定.

由于方案1的VR与SNR值明显小于其它4个方案，故舍掉方案1；在方案2、3、4、5中，根据CN值可以将结果分成2组，即方案2、5和方案3、4，依据CN值最小的原则，舍去值较大的方案2、5，最终选择FMVR、STVAR值相对较小的方案3作为最优反演方案，它的结果与GCMT解的K-angle之差为21°，表明二者的双力偶部分基本类似.在后续的全矩张量反演与纯双力偶反演中均使用该方案.图 6—8展示出M_S5.8地震矩张量反演过程的波形拟合、矩心深度和时空矫正情况，可知矩心深度在4 km处，矩心时间相对于发震时刻向后移动2个时间步长约0.4 s时，得到矩张量反演的最优解，图中显示观测波形和合成波形高度拟合，VR值达到92%，波形互相关系数接近于1，DC分量非常低，约44.7%.

图 6

Fig. 6

图 6 MS5.8地震偏矩张量解的拟合波形与观测波形对比 Fig. 6 Comparison between observed and synthetic waveforms in the seismic moment tensor solution of the MS5.8 event

图 7

Fig. 7

图 7 MS5.8地震矩心深度与波形互相关系数 Fig. 7 Centroid depth and cross-correlation coefficient of the MS5.8 event

按照方案3的台站组合，对其它4个地震进行偏矩张量反演，并将结果进行分解(表 4).一个显著的现象是所有矩张量解具有较低的DC分量，在43%~71%之间，其中3个地震在50%以下，表明震源机制的断层错动成分很少；另一个显著的特征是较浅的矩心深度(4 km).第六列是GCMT发布的前4次地震的矩张量分解后的DC分量，虽然相对于本研究结果均较高，但所有地震的DC分量均低于70%，明显低于典型构造地震双力偶分量之值，说明前郭地震并非典型的构造地震.

表 4 Table 4

表 4 5个MS≥5.0地震的偏矩张量反演结果 Table 4 The results of the Deviational Moment Tensor inversion of the five MS≥ 5.0 events No. MS MW Depth(km) Mrr/Mtt/Mpp/Mrt/Mrp/Mtp CLVD (%) DC (%) GCMT* DC(%) 1 5.5 5.0 4.0 3.116/1.427/-4.534/-0.506/1.250/0.245 56.2 43.8 62 2 5.0 4.8 4.0 0.591/1.465/-2.056/-0.085/-0.469/0.868 53.6 46.4 65 3 5.3 4.8 4.0 1.106/0.508/-1.614/-0.338/-0.755/0.430 38.6 61.4 69 4 5.8 5.1 4.0 2.814/2.253/-5.067/-1.114/1.418/-0.113 55.3 44.7 65 5 5.0 4.4 4.0 4.925/0.418/-5.343/0.622/-1.017/1.686 28.7 71.3 GCMT*，GCMT只发布了前4个地震的矩张量解.

表 4 5个MS≥5.0地震的偏矩张量反演结果 Table 4 The results of the Deviational Moment Tensor inversion of the five MS≥ 5.0 events

为了考察震源机制的体积变化成分，我们用ISOLA方法对5次M_S≥5.0地震进行全矩张量(ISO+CLVD+DC)反演，采用VM-1、VM-2、VM-3三种不同的地壳速度模型(表 1和图 3)建立格林函数，考察不同地壳模型对反演结果的影响.由4.2节中的讨论，本研究对5次地震的全矩张量反演采用台站组合方案3，3个速度模型的矩张量解分解结果均显示70%以下的较低DC分量.由于偏矩张量反演结果表 4中显示M_S5.3地震的DC分量高于其它地震，本文仅以M_S5.3地震为例，讨论其在3个速度模型中的全矩张量反演结果(表 5)，结果显示VR值在90%以上，表明理论波形与观测波形具有很高的拟合程度.另外，矩震级、矩心深度的反演结果在3个速度模型中均相似，DC分量非常稳定，仅在60%~70%之间小范围内变化.3个速度模型的全矩张量反演结果显示，不同的地壳速度模型对反演结果几乎没有影响.表 6列出了5个M_S≥5.0地震采用VM-1的全矩张量反演结果，DC分量均很低，矩心深度为3~4 km.

表 5 Table 5

表 5 MS5.3地震3个速度模型的全矩张量反演结果 Table 5 The results of the MS5.3 event moment tensor inversion using the three velocity models VM MW VR (%) Depth (km) Mrr/Mtt/Mpp/Mrt/Mrp/Mtp ISO (%) DC (%) CLVD (%) VM-1 4.8 92 4.0 2.181/0.826/-1.300/-0.315/-0.741/0.435 23.4 60.7 15.9 VM-2 4.8 91 4.0 4.523/2.535/-5.919/-1.809/2.573/-0.114 5.8 60.1 34.1 VM-3 4.9 90 4.0 7.713/5.386/-4.982/-3.656/3.445/-0.105 24.8 67.9 7.3

表 5 MS5.3地震3个速度模型的全矩张量反演结果 Table 5 The results of the MS5.3 event moment tensor inversion using the three velocity models

表 6 Table 6

表 6 5个MS≥ 5.0地震的全矩张量反演结果 Table 6 The results of the full moment tensor inversion of the five MS≥5.0 events No. MS MW Depth (km) Mrr/Mtt/Mpp/Mrt/Mrp/Mtp ISO (%) DC (%) CLVD (%) 1 5.5 5.0 4.0 4.043/1.748/-4.235/-0.445/1.530/0.226 11.5 47.0 41.5 2 5.0 4.8 4.0 0.652/1.472/-2.035/0.204/-0.080/0.877 1.3 46.9 51.8 3 5.3 4.8 3.0 2.181/0.826/-1.300/-0.315/-0.741/0.435 23.4 60.7 15.9 4 5.8 5.1 4.0 1.106/1.441/-6.365/-0.807/0.859/-0.112 -19.7 25.5 54.8 5 5.0 4.4 3.0 1.439/-1.278/-7.323/0.496/-0.802/1.925 -29.9 28.8 41.3

表 6 5个MS≥ 5.0地震的全矩张量反演结果 Table 6 The results of the full moment tensor inversion of the five MS≥5.0 events

图 8

Fig. 8

图 8 MS5.8地震矩张量在矩心深度平面内搜索的结果 Fig. 8 Grid search results for the MS5.8 event in the plane of the centroid depth

GCMT和吴微微等(2014)曾给出了前郭震群震级较大地震的偏矩张量解，两者的结果具有很好的一致性，而本研究使用全矩张量反演方法得出相对前人研究具有更低的DC分量(表 6)，大多在60%~70%之间，充分说明这些地震与普通构造地震的矩张量解(DC分量大于80%)有很大的差异.然而，在地震矩张量反演中，地震定位、速度模型、波形信噪比等因素会导致矩张量解的不稳定性(Chapman and Leaney, 2012；Adamová and Šileny，2013).为了进一步考察解的不稳定性，需要使用F-test对各个震源类型的矩张量解的显著性进行检验、以及用Jackknifing方法对模型参数的不确定度进行评估(Dreger et al., 2000；Templeton and Dreger, 2006；Adamova et al., 2008).另外，还可以将矩张量解投影到Hudson震源类型图上(Hudson et al., 1989)，以对震源机制类型做进一步判断.最后，我们利用震源深度对前郭震群震源机制进行讨论.

我们使用F-test比较不同的震源类型的统计显著性，包括纯DC、DC+CLVD(偏矩张量)和ISO+DC+CLVD(全矩张量)三种反演模式获得的矩张量解.波形拟合差e_i使用以下公式：

(4)

式中d_i为反演中用到的观测地震波形，s_i为合成地震波形，i为一个数据点或时段.估计方差σ_X²表示如下：

(5)

式中X表示模型类型，分别为DC、DC+CLVD和ISO+ DC+CLVD，t为参与反演的台站数量，N和M是观测量数和模型的参数个数.构造如下3个统计量F₁，F₂和F₃：

(6)

显然，这些统计量服从自由度为(N₁-M₁ -1，N₂-M₂ -1)(Minke，1989)的F分布.检验结果见表 7，F_*α=0.05一列表示相应统计量置信水平为95%的临界值.统计量F₁，F₂和F₃的值均大于对应的临界值，可见5个M_S≥5.0地震的最佳震源类型为ISO+DC+CLVD，即全矩张量反演结果.

表 7 Table 7

表 7 5个地震的F检验结果及最佳震震源模型 Table 7 F-test and the optimal source model of the five events No. MS MW F1 F1α=0.05 F2 F2α=0.05 F3 F3α=0.05 Best Source Model 1 5.5 5.0 3.02 2.09 1.80 1.79 2.00 1.79 ISO+ DC+CLVD 2 5.0 4.8 5.35 2.04 1.62 1.31 1.88 1.31 ISO+ DC+CLVD 3 5.3 4.8 2.05 1.97 1.89 1.48 1.50 1.48 ISO+ DC+CLVD 4 5.8 5.1 3.29 2.09 2.55 1.79 2.35 1.79 ISO+ DC+CLVD 5 5.0 4.4 8.00 1.91 2.80 1.64 2.01 1.64 ISO+ DC+CLVD

表 7 5个地震的F检验结果及最佳震震源模型 Table 7 F-test and the optimal source model of the five events

在矩张量反演的参数不确定度分析中，通常使用Jackknifing(刀切法)、Bootstrap(自举法)等，本研究使用Jackknifing方法(Sokos and Zahradník，2013；Boyd et al., 2015).我们选择盆地内部10个地震台站，在方案3(表 3)的基础上，按每次增减一个台站、增减一个通道进行多次反演，用95%的置信度来确定参数的置信区间.以M_S5.8地震为例，图 9a、9b、9c是双力偶部分的节面参数(走向、倾向、滑动角)的直方图与置信区间，图 9d—9f是矩张量分解后ISO、DC和CLVD分量的直方图与置信区间，可以看出参数的区间估计值较为集中，表明反演结果很稳定.全矩张量双力偶部分的断层面解的2个节面参数分别为走向325°/221°，倾向37°/155°和滑动角70°/56°，它们与震源模型的三分量(表 5)均位于置信区间内，表明反演结果较为可信.

图 9

Fig. 9

图 9 MS5.8地震全矩张量反演参数不确定度的估计 (a—c)双力偶解的断层面参数(走向、倾向和滑动角)不确定度估计；(d—f)全矩张量各分量(ISO分量、DC分量和CLVD分量)的不确定度估计. Fig. 9 The stability assessment of the full moment tensor parameters of the MS5.8 event (a—c) The uncertainties of fault plane parameters, the strike, dip and rake of the double couple solution;(d—f) The uncertainty of the ISO、DC and CLVD components of the full moment tensor.

本研究还使用Hudson震源类型图进一步确定震源矩张量类型.由于矩张量的分解方法不是唯一的(Julian et al., 1998)，即MT=ISO+DEV，而DEV(偏矩张量)可以进一步分解，因此我们根据所研究的问题选择Vavryčuk(2001)提出的分解方法，即将DEV分解为DC和CLVD，相当于人为的加入了先验约束.这种对震源类型的先验约束，可更为客观地将矩张量解投影到Hudson震源类型图上进行研究.Hudson震源类型图是由T和k两参数表示的二维矩张量投影图，具体计算公式如下：

(7)

(8)

式中，M_ISO是矩张量的ISO分量，M′_x、M′_y和M′_z是偏矩张量的3个分量，T表示震源的偏矩张量部分，取值范围为[-1，+1]，即纯负CLVD到纯正CLVD.k表示震源的ISO部分，在图上用-V和+V表示，取值范围为[-1，+1]，即震源的收缩与膨胀.T与k的交点坐标为(0，0)，表示纯DC解.LVD(Linear vector dipole)表示震源区只有一个方向的拉伸与压缩，是CLVD的一种特殊情况，Crack(Tensile crack)是ISO一种特殊情况，表示震源的不对称膨胀与收缩，在某一方向变化量较大.我们将前郭矩张量结果投影到Hudson震源类型图上(图 10)，图中显示大部分地震的投影均远离纯双力偶圆点位置，靠近-CLVD与-LVD，表明震源的力学特征中包含了很少的断层位错机制.

图 10

Fig. 10

图 10 前郭震群5个MS≥ 5.0地震的Hudson震源类型图 Fig. 10 Hudson source type plot of the five MS≥ 5.0 of Qianguo siesmic swarm

通过F-test显著性检验，Jackknifing对模型参数稳定性检验，Hudson震源类型图的震源类型分析，充分说明反演结果中DC分量较稳定，研究获得的全矩张量结果较为可靠.

本文通过开展前郭M_S5.8震群中较大的5个5级以上地震矩张量反演研究，获得了它们的震源机制解.矩张量解的显著特征是很大的非DC分量和很浅的矩心深度，而高的非DC分量表明地震的破裂机制包含了更多的非断层位错因素.

地震矩张量解中相当大的体积变化分量及很浅的震源深度均为采矿区诱发地震的主要特点(Benetatos et al., 2013；Cesca et al., 2014)，诱发地震大多为塌陷或爆炸型震源机制，震源区体积变化较大，通常利用矩张量结果中的非DC分量来判断(Cesca et al., 2014).现有研究表明，大规模开采石油和天然气(包括注水)会引起地下流体运移，改变地下深处孔隙压力，造成有效应力的变化，从而诱发地震(刘元生等，2004).油田作业可诱发5~6级的中等地震，如1952年美国加州凡迪茨附近油田的5.4级地震和1969年3月31日中东苏伊士湾的6.1级地震，均为采油活动诱发的地震(刘元生等，2004).前郭震群震中区位于油田开采区内，附近分布有多口采油井，油田作业很可能在局部区域诱发地震活动.通过本研究区2006年吉林乾安M_S5.0地震的震源深度研究，认为该地震的发生与油气开采有关(高金哲等，2013).综合2013年前郭5次M_S≥5.0地震的矩张量解、震源深度，本研究认为尽管这些地震处于太平洋板块俯冲带的弧后构造区，但2013年前郭M_S5.8震群并非典型的构造地震，而有可能是诱发地震.

捷克查理大学的Jiri教授非常热心的在矩张量反演过程中提供了指导；吉林省地震局台网中心提供了地震波形资料；匿名审稿专家针对研究方法及文章写作提出了非常宝贵的修改意见和建议，在此一并表示感谢！