关鹏虎，李 斌，1b，李自红，宋美琴，梁向军

(1.太原理工大学 a.矿业工程学院，b.地震与地质灾害防治研究所，太原 030024；2.山西省地震局，太原 030021)

地震是地下介质(岩层)在应力作用下发生脆性破裂而引起地面震动的一种自然现象，而震源机制解是推断地震发生时震源断层类型、理解断层破裂力学过程与力学机制的重要依据，也是反映震源区应力状态、地震孕育过程的重要参数。关于如何求解震源机制，国内外地震学者做了大量研究工作并提出了不同的求解方法，如P波初动法[1]、振幅比法[2]、P波初动联合振幅比法[1，3]、CAP全波形反演法[4]与矩张量反演法[5]等等，且研究结果对世界上很多地震都作出过较为合理的解释[6-10]。

但较多的震例研究也显示，对于同一地震，基于不同方法求解的震源机制常常存在一定的差异[11-12]。主要原因在于不同求解方法基于的理论基础、逻辑算法、数据类型不同，以及不同研究者对数据的选用标准与用于计算的数据量等存在差异。此外，不同区域记录台站的疏密不同及记录数据的质量差异，可能会进一步增加震源机制的求解误差，从而影响解的精确性与可靠性。特别对于发震频率高的中小地震，因其震级小、记录台站相对少、台站对震中包裹性不理想、波形信噪比不高等问题[13]，其震源机制的准确求解一直是震源参数研究中的一个关键基础问题，也是当前各地震监测中心经常遇到的一个值得深入分析与探讨的实际问题。

中小地震的震源机制求解方法中，目前普遍认可和广泛应用的方法包括P波初动法、P波初动联合振幅比法、CAP波形反演法与矩张量反演法等。针对中小地震震源机制求解的可靠性问题，本文基于上述多种方法与不同数据类型、数据量，系统计算与对比分析了近年来发生在山西地震带的3次代表性中小地震的震源机制解，展示与验证了参与计算的数据类型、数据量及求解方法对中小地震震源机制解可靠性与稳定性的影响与约束，以期为中小地震震源机制解的可靠性分析及不同震级地震的震源机制求解条件与方法选择提供参考依据。此外，鉴于上述求解方法基本原理清晰且以往文献均有详细描述，本文未再赘述。

1 资料选取与计算思路

1.1 地震资料选取

山西地震带历史上强震频发，中小地震活跃，近年来先后发生2009年3月28日原平4.2级、2010年1月25日河津4.8级、2010年6月7日阳曲4.6级、2016年3月12日盐湖4.8级、2016年4月7日原平4.7级、2016年12月8日清徐4.6级等一系列显著地震[14-16]。区域内的山西测震台网，目前由72个实时记录与传输的数字化宽频带地震台站组网而成(山西省内台站57个和邻省15个)，基本上均匀分布于山西省及周边地区，如图1所示。现有台网条件下，ML2.0及以上的地震，通常能被周围几十个台站同时记录到，震中亦能被台站较好地包围[13]。因此，山西地震带一直以来是开展中小地震研究的理想目标区域之一。

图1 山西区域测震台网及3次代表性地震分布图Fig.1 Distribution of Shanxi seismic network and the three representative events selected for this study

本文选取了近年来发生在山西地震带的3次震级不等的代表性地震为研究对象，分别为2016年1月3日侯马ML2.0地震、2015年6月2日太原ML3.2地震与2010年6月5日阳曲ML4.8地震。选取原因：1)发生于山西地震带中部，记录台站多，台站对震中包裹性好；2)记录波形质量好，信噪比高；3)震级不等，能作为不同震级地震的代表。这3次地震能满足系统分析与探讨不同方法求解中小地震震源机制解结果可靠性的数据要求。

1.2 地壳速度结构

以往研究表明[10，14]，无论哪种震源机制求解方法，地壳速度结构都是影响震源机制解的重要因素之一。为减少地壳速度结构对求解精度与可靠性的影响，本文在利用不同方法求解震源机制解时，均采用了山西地震带最新的1D速度结构模型[14]，地壳分层及数据详见表1。

表1 本文计算采用的山西地震带1D速度结构模型Table 1 1D velocity model of Shanxi rift system applied in this study

1.3 计算思路

本文旨在分析与探讨实际求解过程中，中小地震获得可靠、稳定震源机制解的求解条件与方法选择。因此在求解过程中，总体上采用了分类计算、对比分析的研究思路，即针对不同的震级，选择适合该震级范围的求解方法系统求解，验证不同数据类型、不同数据数量、不同方法对震源机制可靠性与稳定性的影响。本文所有结果均采用下半球投影。

2 多方法求解与对比分析

2.1 侯马ML2.0地震

2016年1月3日山西侯马发生ML2.0地震，山西地震台网共有28个台站记录到该次地震，其中16个台站的记录波形在垂直分量上有清晰的P波初动。对于此次侯马ML2.0地震，采用了计算小震震源机制解最为广泛应用和认可的P波初动法、P波初动联合振幅比法2种求解方法。

2.1.1基于P波初动

HAVSKOV和OTTEMÖLLER[17]基于理论波形的P波初动数据分析震源机制解时表明，若利用P波初动获得可靠性可被接受的震源机制解，至少需要6个对震中包裹较好的台站记录到较为清晰的P波初动。鉴于这一认识，对于此次侯马ML2.0地震，首先选取了震中距最小、方位分布最好的5个台站起算，基于FOCMEC程序[3]求解了震源机制解。求解中设置矛盾数为零、格点搜索步长为2°，并将结果与16个台站全方位覆盖条件下的解进行对比分析。经多次不同台站组合与反复计算，结果如图2所示。其中，图2(a)为两组代表性结果，基于5个台站的记录，无论哪个组合，其结果都与台站全方位覆盖条件下的求解结果相差甚远，都无法求得可靠的震源机制解。

图2(b)是在5个震中距最小、方位分布最好的台站基础上，增加1个台站后获得的震源机制解结果。该台站的增加，有助于改善台站的总体方位角分布。增加不同台站的多次求解结果显示，在此情况下震源机制解至少存在两个最优解，一个是与台站全方位覆盖条件下震源机制解结果相近，另一个则相差甚远。问题在于在实际求解过程中，无法分辨哪个结果是最接近实际的可靠结果。因此，6个分布较好的台站记录到清晰的P波初动能求解可靠的震源机制解，只是理论上的分析，而在实际的求解过程中往往难以操作。

在上述6个台站基础上，再增加1个记录台站，以进一步改善台站的方位分布条件。本文共尝试了5组不同的组合，图2(c)是7个台站条件下求解震源机制的代表性结果。结果显示，即使增加到7个方位较好的台站，其结果仍然存在两种可能性，难以选择可靠结果的问题仍然突出。同样的思路增加到8个台站，台站方位分布进一步得到改善，结果见图2(d)，同样存在两种结果难以选择。图2(e)、2(f)是将记录台站增加到10个、12个时的2组代表性震源机制结果，可以看出其结果完全趋于稳定，且与台站全方位覆盖条件下震源机制解结果非常接近。

2.1.2基于P波初动联合振幅比

为进一步验证振幅比在震源机制求解过程中对结果的约束，本文仍以台站记录较好的侯马ML2.0地震为例，首先以6个震中距最小、方位分布最好的台站记录到的P波初动为基础，逐步增加另外2、3、4、5个震中距较小台站的振幅比数据，并将计算结果与台站全方位覆盖条件下的结果进行对比。计算方法仍然采用基于Snoke方法的FOCMEC程序[3]，数据处理较为简单，即在垂向上读取P波极性，在垂向、径向和切向上读取P波、SV波与SH波振幅，用到的振幅比包括SV/P、SH/P、SV/SH.

(a)5个台;(b)6个台;(c)7个台;(d)8个台;(e)10个台;(f)12个台；STR为走向，DIP为倾角，RAK为滑动角，Source为求解方法，蓝虚线是基于全部16个台P波初动的求解结果，红实线是逐步增加P波初动的求解结果图2 侯马ML2.0地震基于P波初动方法的求解结果Fig.2 Results of focal mechanism of Houma ML2.0 earthquake based on the P-wave first motion method

计算结果如图3所示。图3(a)显示在6个P波初动的基础上，增加2个台站的6个振幅比数据后，通过网格搜索在一定误差范围内可以获得基本稳定的计算结果，且与台站全方位覆盖条件下震源机制解结果具有可比性。表明在P波初动数据不足的情况下增加振幅比，可以对结果起到限定作用。图3(b)-3(d)显示在6个P波初动的基础上，逐步增加振幅比数量，求得震源机制解的变化情况。结果表明，增加振幅比数据会显著增加矛盾比，同时需要不断增加振幅比误差才有解，对结果的约束与可靠度提升并不明显。此外，尝试在8个P波初动的基础上，逐步增加另外2、3、4、5个震中距较小台站的振幅比数据求解震源机制解，同样显示了这一问题。

2.2 太原ML3.2地震

2015年6月2日的太原ML3.2地震发生在山西地震带的中部，山西地震台网几乎所有的地震台站都记录到了该次地震，地震速报震源深度16 km.鉴于此次地震震级略大(3级以上)、记录台站数目较多且近震波形清晰，除采用P波初动法外，也尝试了采用基于波形反演的CAP方法与矩张量反演方法求解其震源机制解，具体计算结果如下。

注：STR为走向，DIP为倾角，RAK为滑动角，Source为求解方法；蓝虚线是基于全部16个台站P波初动的求解结果，红实线是在6个台P波初动的基础上逐步增加振幅比的结果图3 侯马ML2.0地震基于P波初动联合振幅比方法的求解结果Fig.3 Results of Houma ML2.0 earthquake based on the combined P-wave first motion method and amplitude ratios

2.2.1基于P波初动

对于此次地震，山西地震台网共有24个台站垂向上记录到了非常清晰的P波初动，如图4(a)所示。这些台站对震中的包裹性好，最大空隙角为26°，完全满足基于P波初动求解可靠震源机制解的条件。故上述24个台站的P波初动数据全部参与计算，未再加入振幅比数据。为确保基于P波初动结果的可信度，分别采用了Snoke与FPFIT方法[3，18]求解，参数设置包括矛盾数为1，格点搜索步长为2°，FPFIT误差为2°。求解结果如图4(b)所示，显示两种方法获得的震源机制解具有可比性。详细的震源机制解参数见表2。

表2 基于不同方法的3次代表性地震的震源机制解参数Table 2 Parameters of focal mechanisms of the three representative earthquakes by different methods

图4 太原ML3.2地震结果Fig.4 Results of Taiyuan ML3.2 earthquake

2.2.2CAP方法求解

此次地震的震中附近台站较多且波形记录清晰，震中100 km范围内的台站有12个，200 km内的台站有23个，定位精度较高，在一类精度范围内(t<0.5 s)。近震台站数据预处理及滤波后波形清晰，基本满足CAP方法的求解条件。

主要数据处理过程包括：

1)选择波形记录质量较好、参与反演的宽频带台站，本文选取了7台站参与反演。

2)进行波形资料预处理，包括去仪器响应、去均值、去倾斜、水平分量旋转至径向和切向，手动标注P波和面波震相，以及按一定时间窗截取震相。本文截取P波长度35 s、面波长度70 s，并对截取后的震相进行带通滤波，P波滤波范围0.02～0.15 Hz，面波滤波范围0.02～0.1 Hz，并对截取后的波形进行两端尖灭处理。

3)理论地震图计算。

4)震源机制解反演。利用准备好的理论和观测波形数据，通过波形拟合在一定参数空间内搜索，得到相对误差最小的解便是最优解。波形拟合时允许P波有5 s，S波有10 s的时间移动，P波权重1.0，面波权重为0.5.CAP方法的波形拟合情况及反演结果见图5，在震源深度20 km处得到最佳震源机制解[13]，具体参数详见表2.

图5 CAP方法求解太原ML3.2地震震源机制解的波形拟合情况及反演结果[13]Fig.5 Waveform fitting and inversion results of focal mechanism of Taiyuan ML3.2 earthquake by CAP method[13]

2.2.3矩张量反演方法

以往震例中，对于ML3.0地震的震源机制求解，更多采用P波初动或初动联合振幅比等方法，矩张量反演应用较少，主要原因在于小震的波形资料往往不能满足滤波的要求，影响后续较好的波形拟合与结果的可信度[15]。考虑到此次太原ML3.2地震的记录台站多、波形记录较好，且震级接近以往矩张量反演震例中的最小震级，本文也尝试了利用矩张量反演方法求解该地震的震源机制解。

主要数据处理过程包括：

1)选择参与反演的地震波形资料(滤波后波形清晰)；

2)对原始波形数据进行预处理，包括去仪器响应、去均值、去倾斜，三分量波形旋转至垂向、切向和径向分量，积分至速度或位移数据以及滤波；

3)基于区域速度结构等参数计算理论格林函数；

4)进行矩张量反演及反演结果置信度的评估。对于本次地震，仅有5个100 km内的台站满足参与矩张量反演的条件。HAVSKOV和OTTEMÖLLER[17]研究表明，对于质量较好的波形数据，利用速度数据或位移数据进行矩张量反演的结果具有很好的一致性。本文在波形数据预处理过程中,直接采用的是速度数据。此外，不同震级的地震频谱特征不同，对于本次ML3.2地震，选择的滤波频带范围为0.02～0.1 Hz.

波形拟合及反演结果如图6所示。显然对于ML3.2地震，理论波形与实际波形的拟合程度不够理想，虽然能给出震源机制结果(具体参数详见表2)，但约化方差VR(Variance Reduction)=16.8%值较低，结果与其他方法差距较大，可信度较低。

图6 矩张量反演方法求解太原ML3.2地震震源机制解的波形拟合情况及反演结果Fig.6 Waveform fitting and inversion results of focal mechanism of Taiyuan ML3.2 earthquake based on moment tensor inversion method

2.3 阳曲ML4.8地震

2010年6月5日阳曲ML4.8显著地震在阳曲和太原地区震感强烈。山西地震台网72个地震台站全部记录到此次地震。本文采用了上述3种不同的方法求解其震源机制解，以实例探讨该震级的中等地震震源机制解的精度与可靠性对计算方法的依赖程度。

2.3.1基于P波初动

对于此次地震，山西地震台网共有36个台站垂直向记录到清晰的P波初动，也完全满足基于P波初动求解可靠的震源机制解的条件。为了确保结果的可信度，同样分别采用了Snoke与FPFIT方法求解了震源机制解，求解过程中参数设置包括P波极性矛盾数为1，格点搜索步长为2°，FPFIT误差为2°，结果显示两种方法获得的震源机制解具有较好的一致性[14]。求解结果见图7，具体参数详见表2.

图7 阳曲ML4.8地震基于P波初动及其他不同方法的震源机制求解结果Fig.7 Focal mechanisms of Yangqu ML4.8 earthquake based on the P-wave first motion and other different method

2.3.2CAP方法求解

图8为CAP方法的波形拟合情况及反演结果。共有7个台站参与反演，波形资料的预处理过程和计算参数的选择与太原ML3.2地震相同，不再赘述。反演结果显示，对于ML4.8地震，波形拟合程度较高，反演结果可信度高，具体震源机制解参数见表2.

图8 CAP方法求解阳曲ML4.8地震震源机制的波形拟合情况及反演结果Fig.8 Waveform fitting and inversion results of focal mechanism of Yangqu ML4.8 earthquake based on CAP

2.3.3矩张量反演方法

对于此次显著地震，LI et al[14]进行了地震矩张量求解，共有11个台站参与了反演，除滤波频带调整为0.02～0.5 Hz外，其余波形资料预处理方法参数选择与太原ML3.2地震的矩张量反演相同。波形拟合情况及矩张量反演结果见图9，理论波形与实际波形拟合程度较高，且在11个台站参与反演的情况下，VR=46.3%，结果符合信度要求[14-15]。具体参数详见表2，显示不同方法获得的震源机制解具有较好一致性。此外，定位结果显示此次地震震中位于交城断裂与系舟山西麓断裂的延长线上，震源机制的两个节面走向与震中附近区域构造基本吻合。

图9 矩张量反演方法求解阳曲ML4.8地震震源机制的波形拟合情况及反演结果[14]Fig.9 Waveform fitting and inversion results of Taiyuan ML4.8 earthquake by moment tensor inversion method[14]

3 讨论与结论

本文基于山西地震台网数字波形资料，利用P波初动、P波初动联合P/SV/SH振幅比、CAP全波形反演与矩张量反演等多种方法，计算了发生在山西地震带的3次震级不等、具有代表性的中小地震的震源机制解。通过由少到多逐步增加参与计算的P波初动与振幅比数据，以及基于不同方法求解结果的对比分析，全面动态展示了参与计算数据的类型、数量及采用方法对中小地震震源机制解可靠性与稳定性的约束。经反复计算与对比分析，取得了以下认识。

1)对于ML2.0及以下小震的震源机制，基于P波初动和P波初动联合振幅比仍是当前最适用的求解方法。实际求解过程中，基于P波初动获得稳定可靠的震源机制解，理想条件下至少需要8个以上对震中包裹性较好的台站记录到较为清晰的P波初动。当震级较小、记录台站数量不足时，往往需要联合振幅比来丰富参加计算的数据量。无论从理论推导、以往震例或本文的计算结果,都表明联合振幅比能在一定程度上提高对计算结果的约束[1，19]。但过多增加振幅比数据的同时会显著增加振幅的矛盾比，需要不断增大振幅比误差才有解，而过大的矛盾比可能会直接影响到震源机制最终解的可靠性。其主要原因，一是振幅比受介质影响较大，特别是SV波又同时受到自由表面反射的较大影响，导致在实际操作中过多地增加振幅比数据时引入了较大的数据误差[17]；二是准确读取振幅比还需要更为准确的震相识别，这对震级较小、信噪比不高的地震来说有时更为困难。因此，在台站记录数量较多、波形记录清晰的情况下，仅基于P波初动即可获得相对可靠的震源断面解，适当增加振幅比可以对结果起到更好的约束，但过多地引入振幅比有时可能并非是最佳的选择。

2)对于ML3.0这一级别的显著地震，其震源机制求解方法除传统的基于P波初动或P波初动联合振幅比方法外，近年来也有不同学者尝试利用CAP波形反演和矩张量反演方法求解震源机制解[10，12]。CAP方法和矩张量反演方法是否是这一级别地震最佳的震源机制求解方法，CAP方法与矩张量反演方法获得可靠性震源机制解的震级下限是多少仍然是求解震源机制解时需要关注的问题。本文以台站记录较为理想的太原ML3.2地震为例，采用3种不同的方法分别进行了震源机制的求解并进行了结果间的对比分析。结果表明对于该ML3.2地震，基于P波初动和CAP波形反演方法均能给出稳定且信度高的结果，但矩张量反演得到的理论波形与实际波形的拟合程度不够理想，VR值偏低，结果置信度较低。

针对矩张量反演方法的震级适用下限问题，笔者系统统计了近年来不同学者基于矩张量反演方法求解不同地区的中小地震震源机制解的计算实例。如LI et al[15]利用矩张量反演求解山西地震带ML>3.0地震的震源机制解时，其震级下限为ML3.7地震；赵韬等[20]利用矩张量反演方法求解了陕西、四川和甘肃省地震台网记录到的165次地震的震源机制解，其计算震例中的震级下限为ML3.5；刘俊清等[21]利用矩张量反演方法求解了2017年9月23日朝鲜ML3.4地震；毛燕等[22]对云南数字地震台网记录的2000年姚安M6.5地震的余震进行了地震矩张量反演，其震级下限为ML3.0等等。尽管前人研究结果中有ML3.0震级的计算实例，但已有大多数震例的震级下限都在ML3.5或更高。结合前人研究结果与本文实例，认为一般情况下对于ML3.5或以上的地震选择矩张量反演方法可能更为理想。

3)对于ML4.8左右的显著地震，记录台站较多、波形数据的信噪比高，满足多种方法求解震源机制的条件[23]。本文实例也证明这一点，即对于震级较大的中等地震，只要记录台站数量满足基本要求，无论采用上述哪种方法一般都能求得较为可靠的震源机制解。基于P波初动的断层面解反映的是初始破裂情况，而基于波形反演的CAP方法与矩张量反演方法,则更好地诠释了整个地震破裂过程，特别对于台网稀疏地区的地震。地震矩张量可以分解为表示断层面剪切错动的双力偶分量(DC)、震源体膨胀或收缩的各向同性分量(ISO)与震源体优势方向的张裂或挤压变形(CLVD)，加之构造地震通常是由断层两盘的剪切位错触发，因此矩张量能更好地反映断层破裂信息。在条件许可的情况下，矩张量反演方法或许是获得可靠震源机制解的理想选择。

4)本文以记录较为理想的3次代表性中小地震为例，从数据的类型、数据量及采用方法上展示了中小地震获得可靠震源机制解的条件和求解建议。震级较小地震的震源机制的准确求解是当前震源参数研究中一个重要而又复杂的问题，除地震波形记录质量、参与计算台站个数、台站对震中的包围程度等因素外，区域地壳速度结构、发震断层的复杂程度、地震定位误差(特别是震源深度误差)、台站与震中的方位误差等，都会对结果准确性与可靠性产生一定影响[23]，相关问题在今后工作中值得进一步研究。