引言

据中国地震台网中心测定，2018年5月28日1时50分52.5秒在吉林省松原市宁江区发生M_S5.7地震，震中位置为（45.27°N，124.71°E），震源深度为13 km。哈尔滨、长春白城等地震感明显，地震造成震中地区部分房屋损毁，烈度达Ⅶ度，无人员伤亡。

松原M_S5.7地震发生于松辽盆地内部，处于NNE-NE向的扶余—肇东断裂和NW向第二松花江断裂交会部位（图1）。该地区曾于1119年发生过M6¾地震（吴戈等，1988；唐雅芝，1988，1990）；近年来，松原地区地震活动强度和频度在东北地区最高，于2006年和2013年扶余—肇东断裂南端分别发生前郭M_S5.0地震和M_S5.8震群（图1）。此次松原M_S5.7地震发生前，从2017年7月23日开始该地区先后发生了4次M_S≥4.0地震及系列小震。此外，松原还是东北地区重要的石油开采区和经济文化聚集地，因此开展本区内地震构造、地震活动性的相关研究具有十分重要的意义。

图  1  研究区地质构造、历史地震及台站分布

Figure  1.  Regional geological structure and distribution of earthquake epicenters and seismic stations

下载: 全尺寸图片 幻灯片

关于2013年前郭M_S5.8震群的研究较多，如吴微微等（2014）、刘双庆等（2015）、刘俊清等（2017a，b）、盛书中等（2017）、李君和王勤彩（2018），这些研究从震源机制、小震精定位和烈度分布等方面取得的结果较为一致，均认为发震断层为NW向的查干泡—道字井断裂。而关于2018年松原M_S5.7地震的研究结果较少，震源机制解显示此次地震为走滑型地震，存在NW和NE向两个节面，与第二松花江断裂和扶余—肇东断裂走向基本一致（杨文等，2018；盘晓东等，2018；李君等，2019）。小震精定位结果显示余震分布无明显方向性，不能确定发震断层（盘晓东等，2018；李君等，2019），而杨文等（2018）的研究结果显示M_S5.7地震序列沿NE向展布。李君等（2019）根据松原M_S5.7地震震源机制解的两个节面和断层性质判断该地震的发震断层为第二松花江断裂。而松原M_S5.7地震等烈度线长轴为NE向（盘晓东等，2018）。上述研究结果对发震构造的认识也不尽一致，因此，确定此次地震的发震断层对深入研究该地区发震构造和应力场特征十分必要。

鉴于此，本文利用黑龙江省地震台网，以及吉林、内蒙古部分台站记录到的宽频带三分量地震波形数据，通过ISOLA软件（Sokos，Zahradník，2008，2013）对2018年5月28日发生在松原的M_S5.7地震事件进行地震三分量全波形矩张量反演，获取此次地震的震源机制解和矩心位置。根据Zahradník等（2008a）的震源-矩心（hypocenter-centroid，缩写为H-C）方法分析了松原M_S5.7地震矩心、震源和节面的相对位置从而推断发震断层。

引言

据中国地震台网中心测定，2018年5月28日1时50分52.5秒在吉林省松原市宁江区发生M_S5.7地震，震中位置为（45.27°N，124.71°E），震源深度为13 km。哈尔滨、长春白城等地震感明显，地震造成震中地区部分房屋损毁，烈度达Ⅶ度，无人员伤亡。

松原M_S5.7地震发生于松辽盆地内部，处于NNE-NE向的扶余—肇东断裂和NW向第二松花江断裂交会部位（图1）。该地区曾于1119年发生过M6¾地震（吴戈等，1988；唐雅芝，1988，1990）；近年来，松原地区地震活动强度和频度在东北地区最高，于2006年和2013年扶余—肇东断裂南端分别发生前郭M_S5.0地震和M_S5.8震群（图1）。此次松原M_S5.7地震发生前，从2017年7月23日开始该地区先后发生了4次M_S≥4.0地震及系列小震。此外，松原还是东北地区重要的石油开采区和经济文化聚集地，因此开展本区内地震构造、地震活动性的相关研究具有十分重要的意义。

图  1  研究区地质构造、历史地震及台站分布

Figure  1.  Regional geological structure and distribution of earthquake epicenters and seismic stations

下载: 全尺寸图片 幻灯片

关于2013年前郭M_S5.8震群的研究较多，如吴微微等（2014）、刘双庆等（2015）、刘俊清等（2017a，b）、盛书中等（2017）、李君和王勤彩（2018），这些研究从震源机制、小震精定位和烈度分布等方面取得的结果较为一致，均认为发震断层为NW向的查干泡—道字井断裂。而关于2018年松原M_S5.7地震的研究结果较少，震源机制解显示此次地震为走滑型地震，存在NW和NE向两个节面，与第二松花江断裂和扶余—肇东断裂走向基本一致（杨文等，2018；盘晓东等，2018；李君等，2019）。小震精定位结果显示余震分布无明显方向性，不能确定发震断层（盘晓东等，2018；李君等，2019），而杨文等（2018）的研究结果显示M_S5.7地震序列沿NE向展布。李君等（2019）根据松原M_S5.7地震震源机制解的两个节面和断层性质判断该地震的发震断层为第二松花江断裂。而松原M_S5.7地震等烈度线长轴为NE向（盘晓东等，2018）。上述研究结果对发震构造的认识也不尽一致，因此，确定此次地震的发震断层对深入研究该地区发震构造和应力场特征十分必要。

鉴于此，本文利用黑龙江省地震台网，以及吉林、内蒙古部分台站记录到的宽频带三分量地震波形数据，通过ISOLA软件（Sokos，Zahradník，2008，2013）对2018年5月28日发生在松原的M_S5.7地震事件进行地震三分量全波形矩张量反演，获取此次地震的震源机制解和矩心位置。根据Zahradník等（2008a）的震源-矩心（hypocenter-centroid，缩写为H-C）方法分析了松原M_S5.7地震矩心、震源和节面的相对位置从而推断发震断层。

1.   地震矩张量反演方法

本文选择Sokos和Zahradník （2008，2013）的ISOLA近震全波形矩张量反演方法计算松原M_S5.7地震的矩张量解。ISOLA是基于单点源或多点源模型利用迭代反褶积方法，通过网格搜索得到地震最佳矩张量解。Kikuchi和Kanamori （1991）基于六元矩张量的反问题公式提出了使用方差减少（variance reduction，缩写为VR）快速评估观测波形与合成波形之间的相关性：VR＝1－E/O，其中E＝∑（O_i－S_i）²，O＝∑O_i²，O和S分别表示观测波形和合成波形。VR为表征观测波形与合成波形拟合度的重要参数，最大值为1，VR越接近1，表明观测波形与合成波形的拟合度越高。反演过程中使用Bouchon （2003）的离散波数法来计算格林函数。通过网格搜索计算质心位置和质心时间。矩张量可以分解为三部分：双力偶（double couple，缩写为DC）、补偿线性矢量偶极子（compensated linear vector dipole，缩写为CLVD）和各向同性或体积分量（isotropic or volumetric component，缩写为ISO）。它们的相对大小以百分比表示，并且ISO＋CLVD＋DC＝1。当不考虑或忽略体积变化即ISO＝0时，得到的解称为偏矩张量解。本文更多考虑双力偶成分的变化，且由偏矩张量计算得到的标量矩、走向、倾角和滑动角结果可信（Zahradník et al，2008b），所以本文进行的矩张量反演均为偏矩张量反演。

波形反演频率的下限和上限分别由噪声（自然噪声或仪器噪声）与给定的震中距和速度模型质量决定（Fojtíková，Zahradník，2014）。信噪比（signal-to-noise ratio，缩写为SNR）根据振幅谱来评估（图2）。首先，由用户选择波形图上的信号时间窗口，然后系统自动选择噪声时间窗口，使其与信号时间窗具有相同长度，并在信号窗开始时结束。因此，噪声窗口位于信号窗口之前。对两个窗口分别计算振幅谱，并对地震动的每个分量分别进行平滑和绘制。信噪比是作为频率函数（信噪比曲线）计算的，用户在反演过程中选择滤波频带时，可以参考信噪比曲线。

图  2  IDR，BNX和GAN台站记录的松原M_S5.7地震波形（a）、傅里叶谱（b）和信噪比（c）

Figure  2.  Seismic swarm of Songyuan M_S5.7 event （a），Fourier amplitude spectra of the signals and noises （b） and the signal-to-noise ratios of three-component and mean values （c） for the stations IDR，BNX and GAN

下载: 全尺寸图片 幻灯片

矩张量的可分辨性由条件数（condition number，缩写为CN）量化，CN是格林函数矩阵的最大奇异值与最小奇异值之比（Zahradník，Custódio，2012），是一个相对度量，该值的大小反映反演结果优劣。该方法细节请参阅Krizova等（2013）。

另外，震源机制变异指数（focal-mechanism variability index，缩写为FMVAR）和时空变异指数（space-time variability index，缩写为STVAR）通常也用来评价地震矩张量解的优劣，具体方法解释详见Sokos和Zahradník （2013）。

1.   地震矩张量反演方法

本文选择Sokos和Zahradník （2008，2013）的ISOLA近震全波形矩张量反演方法计算松原M_S5.7地震的矩张量解。ISOLA是基于单点源或多点源模型利用迭代反褶积方法，通过网格搜索得到地震最佳矩张量解。Kikuchi和Kanamori （1991）基于六元矩张量的反问题公式提出了使用方差减少（variance reduction，缩写为VR）快速评估观测波形与合成波形之间的相关性：VR＝1－E/O，其中E＝∑（O_i－S_i）²，O＝∑O_i²，O和S分别表示观测波形和合成波形。VR为表征观测波形与合成波形拟合度的重要参数，最大值为1，VR越接近1，表明观测波形与合成波形的拟合度越高。反演过程中使用Bouchon （2003）的离散波数法来计算格林函数。通过网格搜索计算质心位置和质心时间。矩张量可以分解为三部分：双力偶（double couple，缩写为DC）、补偿线性矢量偶极子（compensated linear vector dipole，缩写为CLVD）和各向同性或体积分量（isotropic or volumetric component，缩写为ISO）。它们的相对大小以百分比表示，并且ISO＋CLVD＋DC＝1。当不考虑或忽略体积变化即ISO＝0时，得到的解称为偏矩张量解。本文更多考虑双力偶成分的变化，且由偏矩张量计算得到的标量矩、走向、倾角和滑动角结果可信（Zahradník et al，2008b），所以本文进行的矩张量反演均为偏矩张量反演。

波形反演频率的下限和上限分别由噪声（自然噪声或仪器噪声）与给定的震中距和速度模型质量决定（Fojtíková，Zahradník，2014）。信噪比（signal-to-noise ratio，缩写为SNR）根据振幅谱来评估（图2）。首先，由用户选择波形图上的信号时间窗口，然后系统自动选择噪声时间窗口，使其与信号时间窗具有相同长度，并在信号窗开始时结束。因此，噪声窗口位于信号窗口之前。对两个窗口分别计算振幅谱，并对地震动的每个分量分别进行平滑和绘制。信噪比是作为频率函数（信噪比曲线）计算的，用户在反演过程中选择滤波频带时，可以参考信噪比曲线。

图  2  IDR，BNX和GAN台站记录的松原M_S5.7地震波形（a）、傅里叶谱（b）和信噪比（c）

Figure  2.  Seismic swarm of Songyuan M_S5.7 event （a），Fourier amplitude spectra of the signals and noises （b） and the signal-to-noise ratios of three-component and mean values （c） for the stations IDR，BNX and GAN

下载: 全尺寸图片 幻灯片

矩张量的可分辨性由条件数（condition number，缩写为CN）量化，CN是格林函数矩阵的最大奇异值与最小奇异值之比（Zahradník，Custódio，2012），是一个相对度量，该值的大小反映反演结果优劣。该方法细节请参阅Krizova等（2013）。

另外，震源机制变异指数（focal-mechanism variability index，缩写为FMVAR）和时空变异指数（space-time variability index，缩写为STVAR）通常也用来评价地震矩张量解的优劣，具体方法解释详见Sokos和Zahradník （2013）。

2.   数据选择及处理

用于地震矩张量反演的波形选自黑龙江省区域测震台网的7个宽频带台站，以及吉林和内蒙古测震台网共享的4个宽频带台站，分布如图1所示。如表1所示，台站所使用的宽频带地震计主要以BBVS-60，Guralp CMG-3ESPC和CTS-1H/VBB为主，采样率为100 Hz，震中距最大为358 km，最小为170 km。

表  1  反演使用的宽频带台站仪器

Table  1.  Broadband station instruments for inversion

序号	台站名	采样率/Hz	地震计	数采	震中距/km
1	IDR	100	BBVS-60	EDAS-24IP	216
2	QAN	100	BBVS-60	EDAS-24IP	269
3	NEH	100	Guralp CMG-3ESPC	EDAS-24IP	358
4	NZN	100	CIS-1H/VBB	EDAS-24IP	290
5	WLT	100	CIS-1H/VBB	EDAS-24L	226
6	GAN	100	CIS-1H/VBB	EDAS-24IP	313
7	WAK	100	BBVS-60	EDAS-24IP	222
8	BNX	100	CIS-1H/VBB	EDAS-24IP	217
9	WUC	100	Guralp CMG-3ESPC	EDAS-24IP	187
10	YST	100	BBVS-60	EDAS-24GN	170
11	FMT	100	BBVS-60	EDAS-24IP	231

下载: 导出CSV 

| 显示表格

本研究对每个台站三分量波形数据进行频谱分析和信噪比分析。图2为部分台站（BNX，GAN和IDR）记录到的M_S5.7地震事件波形和频谱图，可以看出，当频率大于0.04 Hz时（图2b和c），所有台站波形均有一个较高的信噪比，因此，选择0.04 Hz为滤波频率下限。但是，个别台站的个别分量（如2a中GAN台站的UD分量）信噪比较差，在反演过程中得到的合成波形与观测波形拟合较差，因此，在优化反演过程中信噪比明显偏低的分量将不参与反演。滤波频率上限是根据地震的震中距和地震波剪切波波长确定，对于震中距小于300 km的台站反演波形频率上限一般为0.1 Hz （Zahradník et al，2005；Benetatos et al，2013），由于本文研究台站，其最大震中距为358 km，所以选择滤波频率上限为0.09 Hz。因此，本文选取的滤波频率范围为0.04—0.09 Hz，为了求取最优解，反演过程中还会在此范围内选择最优滤波频率。

地壳速度结构能够提供有关震源及震源介质的重要信息，是影响矩张量解的重要因素。本文对比了CRUST1.0、吴微微等（2014）和刘俊清等（2017a）给出的吉林松原地区一维速度模型（图3），通过试错，使用ISOLA软件在其它参数不变情况下，计算了四种速度模型下松原M_S5.7地震的震源机制解，结果列于表2，可以看出，由不同速度模型得到的结果相差很小，本文选用具有较高VR值和DC值的吴微微等（2014）的一维速度模型作为反演的速度模型。

表  2  不同速度模型计算出的地震矩张量结果对比

Table  2.  Comparison of moment tensor solutions calculated by different velocity models

速度模型	节面	走向/°	倾角/°	滑动角/°	VR	DC	MW	矩心深度/km
CRUST1.0	Ⅰ	217	77	162	64%	60.8%	5.1	6
	Ⅱ	311	73	14
吴微微等（2014）	Ⅰ	216	75	159	74%	73.3%	5.1	6
	Ⅱ	312	70	16
刘俊清等（2017a）VM-1	Ⅰ	216	77	160	70%	70.4%	5.1	6
	Ⅱ	311	71	13
刘俊清等（2017a）VM-3	Ⅰ	219	85	171	56%	53.9%	5.2	6
	Ⅱ	310	81	5

下载: 导出CSV 

| 显示表格

图  3  松辽盆地一维速度模型对比

Figure  3.  Comparison of velocity models for Songliao basin

下载: 全尺寸图片 幻灯片

2.   数据选择及处理

用于地震矩张量反演的波形选自黑龙江省区域测震台网的7个宽频带台站，以及吉林和内蒙古测震台网共享的4个宽频带台站，分布如图1所示。如表1所示，台站所使用的宽频带地震计主要以BBVS-60，Guralp CMG-3ESPC和CTS-1H/VBB为主，采样率为100 Hz，震中距最大为358 km，最小为170 km。

表  1  反演使用的宽频带台站仪器

Table  1.  Broadband station instruments for inversion

序号	台站名	采样率/Hz	地震计	数采	震中距/km
1	IDR	100	BBVS-60	EDAS-24IP	216
2	QAN	100	BBVS-60	EDAS-24IP	269
3	NEH	100	Guralp CMG-3ESPC	EDAS-24IP	358
4	NZN	100	CIS-1H/VBB	EDAS-24IP	290
5	WLT	100	CIS-1H/VBB	EDAS-24L	226
6	GAN	100	CIS-1H/VBB	EDAS-24IP	313
7	WAK	100	BBVS-60	EDAS-24IP	222
8	BNX	100	CIS-1H/VBB	EDAS-24IP	217
9	WUC	100	Guralp CMG-3ESPC	EDAS-24IP	187
10	YST	100	BBVS-60	EDAS-24GN	170
11	FMT	100	BBVS-60	EDAS-24IP	231

下载: 导出CSV 

| 显示表格

本研究对每个台站三分量波形数据进行频谱分析和信噪比分析。图2为部分台站（BNX，GAN和IDR）记录到的M_S5.7地震事件波形和频谱图，可以看出，当频率大于0.04 Hz时（图2b和c），所有台站波形均有一个较高的信噪比，因此，选择0.04 Hz为滤波频率下限。但是，个别台站的个别分量（如2a中GAN台站的UD分量）信噪比较差，在反演过程中得到的合成波形与观测波形拟合较差，因此，在优化反演过程中信噪比明显偏低的分量将不参与反演。滤波频率上限是根据地震的震中距和地震波剪切波波长确定，对于震中距小于300 km的台站反演波形频率上限一般为0.1 Hz （Zahradník et al，2005；Benetatos et al，2013），由于本文研究台站，其最大震中距为358 km，所以选择滤波频率上限为0.09 Hz。因此，本文选取的滤波频率范围为0.04—0.09 Hz，为了求取最优解，反演过程中还会在此范围内选择最优滤波频率。

地壳速度结构能够提供有关震源及震源介质的重要信息，是影响矩张量解的重要因素。本文对比了CRUST1.0、吴微微等（2014）和刘俊清等（2017a）给出的吉林松原地区一维速度模型（图3），通过试错，使用ISOLA软件在其它参数不变情况下，计算了四种速度模型下松原M_S5.7地震的震源机制解，结果列于表2，可以看出，由不同速度模型得到的结果相差很小，本文选用具有较高VR值和DC值的吴微微等（2014）的一维速度模型作为反演的速度模型。

表  2  不同速度模型计算出的地震矩张量结果对比

Table  2.  Comparison of moment tensor solutions calculated by different velocity models

速度模型	节面	走向/°	倾角/°	滑动角/°	VR	DC	MW	矩心深度/km
CRUST1.0	Ⅰ	217	77	162	64%	60.8%	5.1	6
	Ⅱ	311	73	14
吴微微等（2014）	Ⅰ	216	75	159	74%	73.3%	5.1	6
	Ⅱ	312	70	16
刘俊清等（2017a）VM-1	Ⅰ	216	77	160	70%	70.4%	5.1	6
	Ⅱ	311	71	13
刘俊清等（2017a）VM-3	Ⅰ	219	85	171	56%	53.9%	5.2	6
	Ⅱ	310	81	5

下载: 导出CSV 

| 显示表格

图  3  松辽盆地一维速度模型对比

Figure  3.  Comparison of velocity models for Songliao basin

下载: 全尺寸图片 幻灯片

3.   矩张量反演结果

首先使用全部11个台站波形，以M_S5.7地震震中位置为中心，在其下方2—15 km 深度（步长1 km），在0.04—0.09 Hz滤波频带内反演矩心位置，得到初步结果，如图4和表3中第一行测试1所示，可见在3—8 km深度矩张量解参数变化不大，具有较好的稳定性，最优解深度为6 km，该深度解显示VR值为74%和较高的双力偶成分占73.3%。该解在深度和时间上具有较大的稳定性，矩张量在最佳走向、倾角和滑动角附近略有变化。由于所选台站中有些震中距较大，有些台站的方位角重叠，而且使用ISOLA进行波形反演对台站方位角的覆盖和台站数量要求不高。于是，在此基础上进一步优化台站选择和滤波频率范围，通过试错，最终反演结果列于表3，可以看出选择NZN，IDR，WLT，WAK，BNX，YST和FMT等7个台站在0.05—0.08 Hz频带内进行反演得到的解进一步优化，震源机制解较为稳定，VR值（88%）和DC值（93.5%）有明显提升。因此，在固定震中位置为矩心位置进行深度反演时的矩张量最优解为：矩心深度为6 km；节面Ⅰ的走向、倾角和滑动角分别为214°，69°和154°；节面Ⅱ的走向、倾角和滑动角分别为313°，66°和23°；DC＝93.5%，VR＝88%，M_W＝5.2。本文反演的两次结果与USGS，GFZ和李君等（2019）结果基本一致（表3），表明反演结果可信。下文进行的反演均在表3测试2的台站选择和滤波频率下进行。

表  3  不同震源机制计算结果对比

Table  3.  Comparison of focal mechanisms from different agency or researchers

来源	MW	矩心深度/km	节面Ⅰ				节面Ⅱ			DC	VR
			走向/°	倾角/°	滑动角/°		走向/°	倾角/°	滑动角/°
测试1	5.1	6	216	75	159		312	70	16	73.3%	74%
测试2	5.2	6	214	69	154		313	66	23	93.5%	88%
李君等（2019）	−	6	220	79	162		314	72	12	−	−
GFZ （2018）	5.2	14	223	77	180		313	90	13	−	−
USGS （2018）	5.3	12	47	68	−169		313	80	−22	64%	−
注：测试1为利用本文选择的全部11个宽频带台站反演的结果，测试2为利用NZN，IDR，WLT，WAK，BNX，YST和FMT等7个 　　宽频带台站反演结果

下载: 导出CSV 

| 显示表格

图  4  松原M_S5.7地震矩心深度与波形互相关系数

Figure  4.  Centroid depth versus cross-correlation coefficient for the M_S5.7 Songyuan event

下载: 全尺寸图片 幻灯片

矩心代表地震破裂区的重心。矩心可以在矩张量反演过程中确定，通过在假定位置周围进行网格搜索找到一个最优的波形拟合的点。对于一次M_S4.0的地震事件，震源-矩心距离与其定位误差相当，甚至更小。因此，对于小的地震事件矩心可视为震源位置。对于M_S5.7地震，其破裂尺度一般在七至十几千米不等，震源与矩心不重合。

在前述确定最佳矩心深度为6 km的基础上，以震中位置为中心，在0.05—0.08 Hz滤波频带内使用偏矩张量反演在水平面上进行网格搜索矩心位置，分别在5 km，6 km，7 km和8 km深度建立11×9个1 km×1 km的网格进行反演矩心位置。观测波形与合成波形相关性空间分布如图5所示，可以明确看出：相对于原始震中定位，矩心发生了偏移；矩心的最佳位置位于震中南侧5 km，西侧2 km处；在深度6 km和7 km处方差减少量最大，即VR＝93%，而这两个深度处得到的解（走向、倾角和滑动角）非常相似。本文选取7 km深度的矩张量解，因为它具有更小的FMVAR值（FMVAR＝12±7）和更高的双力偶成分（96.4%）。该深度矩心位置为（45.225°N，124.685°E）；节面Ⅰ的走向、倾角和滑动角分别为217°，82°和164°，节面Ⅱ的走向、倾角和滑动角分别为309°，74°和8°，DC＝96.4%，VR＝93%，M_W5.2 （图6）。合成波形与观测波形的拟合如图7所示，其中台站IDR的NS向波形未参与反演。

图  5  不同深度网格平面观测波形与合成波形之间相关性

Figure  5.  Plot of correlation between observed and synthetic records for horizontal grids with different depth

下载: 全尺寸图片 幻灯片

图  6  ISOLA软件反演松原M_S5.7地震矩张量最终结果

Figure  6.  The final centroid position of Songyuan M_S5.7 event from seismic moment inversion by ISOLA

下载: 全尺寸图片 幻灯片

图  7  图6中矩张量解的波形拟合情况

右上角蓝色数字代表单测项的方差减少，台站IDR的南北向波形未在反演中使用

Figure  7.  Waveform fit for the centroid solution shown in Fig.6

Blue numbers mark the variance reduction of the individual components and the NS waveforms of station IDR were not used in the inversion

下载: 全尺寸图片 幻灯片

3.   矩张量反演结果

首先使用全部11个台站波形，以M_S5.7地震震中位置为中心，在其下方2—15 km 深度（步长1 km），在0.04—0.09 Hz滤波频带内反演矩心位置，得到初步结果，如图4和表3中第一行测试1所示，可见在3—8 km深度矩张量解参数变化不大，具有较好的稳定性，最优解深度为6 km，该深度解显示VR值为74%和较高的双力偶成分占73.3%。该解在深度和时间上具有较大的稳定性，矩张量在最佳走向、倾角和滑动角附近略有变化。由于所选台站中有些震中距较大，有些台站的方位角重叠，而且使用ISOLA进行波形反演对台站方位角的覆盖和台站数量要求不高。于是，在此基础上进一步优化台站选择和滤波频率范围，通过试错，最终反演结果列于表3，可以看出选择NZN，IDR，WLT，WAK，BNX，YST和FMT等7个台站在0.05—0.08 Hz频带内进行反演得到的解进一步优化，震源机制解较为稳定，VR值（88%）和DC值（93.5%）有明显提升。因此，在固定震中位置为矩心位置进行深度反演时的矩张量最优解为：矩心深度为6 km；节面Ⅰ的走向、倾角和滑动角分别为214°，69°和154°；节面Ⅱ的走向、倾角和滑动角分别为313°，66°和23°；DC＝93.5%，VR＝88%，M_W＝5.2。本文反演的两次结果与USGS，GFZ和李君等（2019）结果基本一致（表3），表明反演结果可信。下文进行的反演均在表3测试2的台站选择和滤波频率下进行。

表  3  不同震源机制计算结果对比

Table  3.  Comparison of focal mechanisms from different agency or researchers

来源	MW	矩心深度/km	节面Ⅰ				节面Ⅱ			DC	VR
			走向/°	倾角/°	滑动角/°		走向/°	倾角/°	滑动角/°
测试1	5.1	6	216	75	159		312	70	16	73.3%	74%
测试2	5.2	6	214	69	154		313	66	23	93.5%	88%
李君等（2019）	−	6	220	79	162		314	72	12	−	−
GFZ （2018）	5.2	14	223	77	180		313	90	13	−	−
USGS （2018）	5.3	12	47	68	−169		313	80	−22	64%	−
注：测试1为利用本文选择的全部11个宽频带台站反演的结果，测试2为利用NZN，IDR，WLT，WAK，BNX，YST和FMT等7个 　　宽频带台站反演结果

下载: 导出CSV 

| 显示表格

图  4  松原M_S5.7地震矩心深度与波形互相关系数

Figure  4.  Centroid depth versus cross-correlation coefficient for the M_S5.7 Songyuan event

下载: 全尺寸图片 幻灯片

矩心代表地震破裂区的重心。矩心可以在矩张量反演过程中确定，通过在假定位置周围进行网格搜索找到一个最优的波形拟合的点。对于一次M_S4.0的地震事件，震源-矩心距离与其定位误差相当，甚至更小。因此，对于小的地震事件矩心可视为震源位置。对于M_S5.7地震，其破裂尺度一般在七至十几千米不等，震源与矩心不重合。

在前述确定最佳矩心深度为6 km的基础上，以震中位置为中心，在0.05—0.08 Hz滤波频带内使用偏矩张量反演在水平面上进行网格搜索矩心位置，分别在5 km，6 km，7 km和8 km深度建立11×9个1 km×1 km的网格进行反演矩心位置。观测波形与合成波形相关性空间分布如图5所示，可以明确看出：相对于原始震中定位，矩心发生了偏移；矩心的最佳位置位于震中南侧5 km，西侧2 km处；在深度6 km和7 km处方差减少量最大，即VR＝93%，而这两个深度处得到的解（走向、倾角和滑动角）非常相似。本文选取7 km深度的矩张量解，因为它具有更小的FMVAR值（FMVAR＝12±7）和更高的双力偶成分（96.4%）。该深度矩心位置为（45.225°N，124.685°E）；节面Ⅰ的走向、倾角和滑动角分别为217°，82°和164°，节面Ⅱ的走向、倾角和滑动角分别为309°，74°和8°，DC＝96.4%，VR＝93%，M_W5.2 （图6）。合成波形与观测波形的拟合如图7所示，其中台站IDR的NS向波形未参与反演。

图  5  不同深度网格平面观测波形与合成波形之间相关性

Figure  5.  Plot of correlation between observed and synthetic records for horizontal grids with different depth

下载: 全尺寸图片 幻灯片

图  6  ISOLA软件反演松原M_S5.7地震矩张量最终结果

Figure  6.  The final centroid position of Songyuan M_S5.7 event from seismic moment inversion by ISOLA

下载: 全尺寸图片 幻灯片

图  7  图6中矩张量解的波形拟合情况

右上角蓝色数字代表单测项的方差减少，台站IDR的南北向波形未在反演中使用

Figure  7.  Waveform fit for the centroid solution shown in Fig.6

Blue numbers mark the variance reduction of the individual components and the NS waveforms of station IDR were not used in the inversion

下载: 全尺寸图片 幻灯片

4.   地震发震断层判断

Zahradník等（2008a）提出了一种简单、有效判断发震断层的方法，称为震源-矩心方法（hypocenter-centroid method，缩写为H-C），其中：震源位置是地震破裂开始的地方；矩心位置是断层上主要滑动区域的点源近似。该方法尤其适用于走滑发震断层的判断。H-C方法的原理是将震源和矩心（在节面Ⅰ与节面Ⅱ交线上）投影到三维空间，然后计算震源到两个节面的距离。如果震源位于一个节面上，那么这个节面就是真实的发震断层。如果震源不在任何一个节面上，那么距离震源最近的那个节面就是真实发震断层。

由于所使用的定位台站和方法不同，震源定位结果也不尽相同，不同机构定位结果见表4。我们使用黑龙江省地震台网的地震台站和内蒙古、吉林台网部分共享台站的资料，利用HYPOINVERSE 2000 （Klein，1989）对松原M_S5.7地震进行了重定位，震中位置为（45.27°N，124.77°E），深度为8.1 km，水平误差为1.1 km，垂直误差为2.0 km，残差为0.18。^①

表  4  松原M_S5.7地震基本参数

Table  4.  Basic parameters of Songyuan M_S5.7 event

编号	发震时间			震中位置		MS	震源深度/km	震中位置	结果来源
	年-月-日	时： 分： 秒		北纬/°	东经/°
1	2018-05-28	01： 50： 52.5		45.31	124.69	5.7	10	松原市宁江区	中国地震台网（2018）
2	2018-05-28	01： 50： 52.6		45.27	124.71	5.7	10	松原市宁江区	刘俊清等①
3	2018-05-28	01： 50： 52.6		45.24	124.64	5.3	10	扶余	USGS（2018）
4	2018-05-28	01： 50： 53.9		45.33	124.42	5.2 （MW）	10	东北	GFZ（2018）
5	2018-05-28	01： 50： 52.5		45.27	124.77	5.7	8.1	松原市宁江区	本文

下载: 导出CSV 

| 显示表格

下面分别使用中国地震台网中心（2018）、刘俊清等^①和本文所获定位结果的震源定位信息，与本文反演得到的矩心位置和两个节面进行分析。图8给出了松原M_S5.7地震事件震源-矩心方法分析结果，显示本文定位结果与矩心的距离是8.43 km，距离节面Ⅰ（绿色）2.46 km，距离节面Ⅱ（红色）7.52 km。刘俊清等^①的定位结果与矩心、节面Ⅰ和距离节面Ⅱ的距离分别为6.17 km，0.98 km和4.12 km，中国地震台网中心（2018）定位结果与矩心、节面Ⅰ和距离节面Ⅱ的距离分别为9.92 km，4.87 km和6.49 km。三个定位结果虽然稍有偏差，但是都显示节面Ⅰ是发震断层面，该节面与扶余—肇东断裂走向一致，所以据此判断扶余—肇东断裂为此次地震的发震断层。

图  8  松原M_S5.7地震震源-矩心位置图

Figure  8.  The hypocenter-centroid plot of Songyuan M_S5.7 event

下载: 全尺寸图片 幻灯片

4.   地震发震断层判断

Zahradník等（2008a）提出了一种简单、有效判断发震断层的方法，称为震源-矩心方法（hypocenter-centroid method，缩写为H-C），其中：震源位置是地震破裂开始的地方；矩心位置是断层上主要滑动区域的点源近似。该方法尤其适用于走滑发震断层的判断。H-C方法的原理是将震源和矩心（在节面Ⅰ与节面Ⅱ交线上）投影到三维空间，然后计算震源到两个节面的距离。如果震源位于一个节面上，那么这个节面就是真实的发震断层。如果震源不在任何一个节面上，那么距离震源最近的那个节面就是真实发震断层。

由于所使用的定位台站和方法不同，震源定位结果也不尽相同，不同机构定位结果见表4。我们使用黑龙江省地震台网的地震台站和内蒙古、吉林台网部分共享台站的资料，利用HYPOINVERSE 2000 （Klein，1989）对松原M_S5.7地震进行了重定位，震中位置为（45.27°N，124.77°E），深度为8.1 km，水平误差为1.1 km，垂直误差为2.0 km，残差为0.18。^①

表  4  松原M_S5.7地震基本参数

Table  4.  Basic parameters of Songyuan M_S5.7 event

编号	发震时间			震中位置		MS	震源深度/km	震中位置	结果来源
	年-月-日	时： 分： 秒		北纬/°	东经/°
1	2018-05-28	01： 50： 52.5		45.31	124.69	5.7	10	松原市宁江区	中国地震台网（2018）
2	2018-05-28	01： 50： 52.6		45.27	124.71	5.7	10	松原市宁江区	刘俊清等①
3	2018-05-28	01： 50： 52.6		45.24	124.64	5.3	10	扶余	USGS（2018）
4	2018-05-28	01： 50： 53.9		45.33	124.42	5.2 （MW）	10	东北	GFZ（2018）
5	2018-05-28	01： 50： 52.5		45.27	124.77	5.7	8.1	松原市宁江区	本文

下载: 导出CSV 

| 显示表格

下面分别使用中国地震台网中心（2018）、刘俊清等^①和本文所获定位结果的震源定位信息，与本文反演得到的矩心位置和两个节面进行分析。图8给出了松原M_S5.7地震事件震源-矩心方法分析结果，显示本文定位结果与矩心的距离是8.43 km，距离节面Ⅰ（绿色）2.46 km，距离节面Ⅱ（红色）7.52 km。刘俊清等^①的定位结果与矩心、节面Ⅰ和距离节面Ⅱ的距离分别为6.17 km，0.98 km和4.12 km，中国地震台网中心（2018）定位结果与矩心、节面Ⅰ和距离节面Ⅱ的距离分别为9.92 km，4.87 km和6.49 km。三个定位结果虽然稍有偏差，但是都显示节面Ⅰ是发震断层面，该节面与扶余—肇东断裂走向一致，所以据此判断扶余—肇东断裂为此次地震的发震断层。

图  8  松原M_S5.7地震震源-矩心位置图

Figure  8.  The hypocenter-centroid plot of Songyuan M_S5.7 event

下载: 全尺寸图片 幻灯片

5.   讨论

通过H-C方法验证了扶余—肇东断裂为松原M_S5.7地震的发震断层，该结果与反演得到的震源机制解的节面I走向、倾角较为一致，且与地震等烈度线长轴方向（盘晓东等，2018）一致。余震的分布情况有助于我们判断发震断层面，但是李君等（2019）和刘俊清等^①的余震精定位分析结果显示余震没有明显的方向性展布，无法帮助判断发震断层。震源机制解显示此次地震为走滑型地震，两个节面分别与扶余—肇东断裂和第二松花江断裂走向一致，并且相关文献（李传友等，1999；刘权峰等，2017；吕晗，2017；刘俊清，2018）显示，对扶余—肇东断裂的描述为逆冲型，看上去似乎与震源机制显示的走滑型相矛盾，李君等（2019）也据此认为第二松花江断裂为发震断层。但是，本文也注意到刘权峰等（2017a，b）和刘俊清（2018）对扶余—肇东断裂最新活动时代描述为晚更新世，而这并不能代表现今的地壳运动状态，中国地震局第一监测中心（2018）的GPS应变图就能够很好地证明（1999—2007年、2009—2013年和2013—2015年三个阶段GPS应变图表明东北地区主应变具有明显变化）。于吉鹏等（2019）的研究表明2012—2017年松原地区主压应力方向为北东向，方向与北东向扶余—肇东断裂的走向一致，因此，有助于断裂发生走滑破裂。常金龙^②利用2018年1—7月GPS基站数据，分别以内蒙古乌兰浩特（NMWL）、吉林长岭（JLCL）和长春（CHAN）三个台站为基准，分析大庆肇源（DQZY）基站的相对位置变化，其结果显示DQZY相对于三个参考点均有北东向运动的特点，长轴方向与扶余—肇东断裂走向基本一致（图9）。考虑到DQZY距离吉林松原M_S5.7地震震中仅几十千米的距离，所以该点的运动状态基本可以反映地震前后震中位置的断裂活动变化。所以GPS资料表明松原M_S5.7地震前后区域应力状态和GPS基站相对变化分析也都支持扶余—肇东断裂为发震断裂的判定。

图  9  松原M_S5.7地震前后GPS基站位置相对变化

上行三小图分别为DQZY站相对于NMWL，JLCL和CHAN三个基站位置变化示意图

Figure  9.  Relative variation of GPS station around Songyuan M_S5.7 event

The three pictures at the top of the figure are the relative changes of DQZY station relative to NMWL，JLCL and CHAN stations，respectively

下载: 全尺寸图片 幻灯片

综上所述，就目前所掌握的信息，本文认为扶余—肇东断裂为2018年5月27日M_S5.7地震的发震断裂。

5.   讨论

通过H-C方法验证了扶余—肇东断裂为松原M_S5.7地震的发震断层，该结果与反演得到的震源机制解的节面I走向、倾角较为一致，且与地震等烈度线长轴方向（盘晓东等，2018）一致。余震的分布情况有助于我们判断发震断层面，但是李君等（2019）和刘俊清等^①的余震精定位分析结果显示余震没有明显的方向性展布，无法帮助判断发震断层。震源机制解显示此次地震为走滑型地震，两个节面分别与扶余—肇东断裂和第二松花江断裂走向一致，并且相关文献（李传友等，1999；刘权峰等，2017；吕晗，2017；刘俊清，2018）显示，对扶余—肇东断裂的描述为逆冲型，看上去似乎与震源机制显示的走滑型相矛盾，李君等（2019）也据此认为第二松花江断裂为发震断层。但是，本文也注意到刘权峰等（2017a，b）和刘俊清（2018）对扶余—肇东断裂最新活动时代描述为晚更新世，而这并不能代表现今的地壳运动状态，中国地震局第一监测中心（2018）的GPS应变图就能够很好地证明（1999—2007年、2009—2013年和2013—2015年三个阶段GPS应变图表明东北地区主应变具有明显变化）。于吉鹏等（2019）的研究表明2012—2017年松原地区主压应力方向为北东向，方向与北东向扶余—肇东断裂的走向一致，因此，有助于断裂发生走滑破裂。常金龙^②利用2018年1—7月GPS基站数据，分别以内蒙古乌兰浩特（NMWL）、吉林长岭（JLCL）和长春（CHAN）三个台站为基准，分析大庆肇源（DQZY）基站的相对位置变化，其结果显示DQZY相对于三个参考点均有北东向运动的特点，长轴方向与扶余—肇东断裂走向基本一致（图9）。考虑到DQZY距离吉林松原M_S5.7地震震中仅几十千米的距离，所以该点的运动状态基本可以反映地震前后震中位置的断裂活动变化。所以GPS资料表明松原M_S5.7地震前后区域应力状态和GPS基站相对变化分析也都支持扶余—肇东断裂为发震断裂的判定。

图  9  松原M_S5.7地震前后GPS基站位置相对变化

上行三小图分别为DQZY站相对于NMWL，JLCL和CHAN三个基站位置变化示意图

Figure  9.  Relative variation of GPS station around Songyuan M_S5.7 event

The three pictures at the top of the figure are the relative changes of DQZY station relative to NMWL，JLCL and CHAN stations，respectively

下载: 全尺寸图片 幻灯片

综上所述，就目前所掌握的信息，本文认为扶余—肇东断裂为2018年5月27日M_S5.7地震的发震断裂。

6.   结论

本文采用ISOLA软件对2018年5月28日吉林松原M_S5.7地震进行了地震矩张量反演，结果显示在5—8 km深度的震源机制解较为稳定，与其他机构结果较为一致，结果可信。进一步使用网格搜索法在5—8 km深度反演搜索最佳矩心位置，最终确定最佳解为：矩心位置为（45.225°N，124.685°E）；矩心深度为7 km；节面Ⅰ的走向、倾角和滑动角分别为217°，82°和164°；节面Ⅱ的走向、倾角和滑动角分别为309°，74°和8°；DC＝96.4%，VR＝93%，矩震级为M_W5.2。

进一步使用H-C方法分析了震源、矩心与节面之间的关系，推断扶余—肇东断裂为此次地震的发震断层。计算得到的震源机制似乎与地震地质资料相矛盾，但是GPS资料分析显示了现今吉林松原地区的应力状态，结果与之前的判断相一致。

总之，本文的计算与分析表明2018年5月28日松原M_S5.7地震的发震断层为NE向扶余—肇东断裂，为走滑型地震。

两位审稿专家对本文修改提出的宝贵的意见，作者在此表示衷心感谢！

6.   结论

本文采用ISOLA软件对2018年5月28日吉林松原M_S5.7地震进行了地震矩张量反演，结果显示在5—8 km深度的震源机制解较为稳定，与其他机构结果较为一致，结果可信。进一步使用网格搜索法在5—8 km深度反演搜索最佳矩心位置，最终确定最佳解为：矩心位置为（45.225°N，124.685°E）；矩心深度为7 km；节面Ⅰ的走向、倾角和滑动角分别为217°，82°和164°；节面Ⅱ的走向、倾角和滑动角分别为309°，74°和8°；DC＝96.4%，VR＝93%，矩震级为M_W5.2。

进一步使用H-C方法分析了震源、矩心与节面之间的关系，推断扶余—肇东断裂为此次地震的发震断层。计算得到的震源机制似乎与地震地质资料相矛盾，但是GPS资料分析显示了现今吉林松原地区的应力状态，结果与之前的判断相一致。

总之，本文的计算与分析表明2018年5月28日松原M_S5.7地震的发震断层为NE向扶余—肇东断裂，为走滑型地震。

两位审稿专家对本文修改提出的宝贵的意见，作者在此表示衷心感谢！