唐山地震的破裂过程及其力学分析

张之立1, 李钦祖2, 谷继成1, 靳雅敏2, 杨懋源3, 刘万琴1

摘要: 由 P 波初动符号资料在 DJS-6机上计算了主震及17个较大余震的断层面解,并按照有限移动源模式测定了主震及三个最大余震的震源参数.主震是发生在一个近似直立的右旋走滑断层上,走向 N30°E,破裂方式为不对称的双侧破裂,以2.7公里/秒的平均速度向北东传播70公里,向南西传播45公里.测定的主震震源参数例如平均位错136厘米,地震矩1.24×10²⁷达因·厘米,应力降12巴等.大多数M_L>5.0的余震是发生在主破裂面附近及主破裂面两端的扩展分支上,该扩展分支位于膨胀符号区并与主破裂偏离80°左右.较大余震的多数亦集中在这两个扩展分支上.本文试图从理论上分析这种断裂扩展的力学特征.对于脆性材料的复合变形情形,破裂不再沿原来平面扩展,而是与原来平面偏离一个角度的另一面内扩展.并提出一个力学模型,计算了断层扩展角,计算结果与观测事实比较吻合.根据以上结果,本文讨论了唐山地震特点及发生的力学条件,认为唐山地震不同于发生在大断层上能用粘滑机制解释的那类地震,它和海城地震类似的是,除水平应力场作用外,还可能有地下物质的变迁,由于这种变迁使局部地壳受到垂直力.它和海城地震不同的是,它发生在一个比较均匀的脆性介质内,因而能够积累能量发生大震而没有前震.

Abstract: Fault-plane solutions of the mainshock and 17 large aftershocks are computed on the DJS-6 computer from the data of the first motions of P waves. Based on the model of finite moving source, the fracture processes and the source parameters of the mainshock and three largest aftershocks are estimated. The Mainshock occurred on a nearly vertical right-lateral strike-slip fault, striking N30°E, with an asymmetric bilateral fracture which propagated 70 km northeastward and 45 km southwestward with an average velocity of 2.7 km/sec. Source parameters of the mainshock such as the average dislocation, seismic moment and stress-drop are estimated to be 136 cm, 1.24×10²⁷ dyne, cm and 12 bars respectively.Most of the aftershocks with M_L>5.0 occurred near the main fault plane and along the fracture branches at both of its ends. These fracture branches are located in the dilatational quarters and make about 80 angles with the mam fault plane. On the two fracture branches, most of the major aftershocks are concentrated. A mechanical analysis of this phenomenon of fracture propagation is attempted in this paper.It is shown that for complex deformations, the fracture propagation no longer takes place along the original plane for brittle materials but deviates from it by a certain angle. The mechanical model used to calculate this deviation angle of fault yields results which agree essentially with the observation.Based on the above results, it can be inferred that the characteristics and mechanical conditions of the Tangshan earthquake as well as those of the Haicheng earthquake are quite different from those of the earthquakes occurring on some great faults which may be explained by stick-slip mechanisms. These two earthquakes are not only caused by the horizontal stress fields but also by local vertical forces induced from the movement of the underground materials. It is also worthy of note that Tangshan earthquake is different from Haicheng earthquake in that it occurs in a relatively homogeneous and brittle medium and therefore no foreshocks took place before the mainshock.

引言

图像信息（pattern informatics，简写为PI）法是一种研究地震活动模式的统计方法。该方法是在希尔伯特（Hilbert）空间创建一个与时间相关的代表地震活动性的系统状态向量，该空间中单位矢量随时间的变化确定了小地震在大地震之前随着时间从一个地点转移到另一个地点的特征模式，在对各个模式学习后，识别一定时间内可能的地震活动模式并给出强震的发生概率，最终输出包含地震发生概率的“热点”图（Rundle et al，2000a，b，2002；Tiampo et al，2002）。该方法经过多次改进（Chen et al，2005；Holliday et al，2005）后，被Holliday等（2006）正式命名为图像信息法，后在多个国家得到应用。 PI方法相对于空间非均匀的泊松模型相对强度（relative intensity，简称 RI）法和均匀泊松模型（随机预测）具有更好的预测效能（蒋长胜，吴忠良，2008），且在中长期的地震预测研究中表现优异（Chen et al，2006；Holliday et al，2006；Nanjo et al，2006a，b；Radan et al，2013；Zhang et al，2013，2017；Mohanty et al，2016；宋程等，2018；尼鲁帕尔·买买吐孙，余怀忠，2020；余怀忠等，2020；Chang et al，2020）。

图像信息法在我国大陆中部地区的强震预测研究中也得到了很好的应用（蒋长胜，吴忠良，2008；Zhao et al，2010；Zhang et al，2013，2017；张小涛等，2014）。蒋长胜等（2009）将强震频发的川滇地区以及安达曼—苏门答腊地区组成统一的强震预测区并对该区域进行回溯性研究，结果显示汶川M_S8.0地震前震中附近出现了中长期尺度的地震活动异常。Zhao等（2010）结合PI和前兆逆向追踪方法（Keilis-Borok et al，2004）研究了川滇地区长期的地震活动性，为中国地震局年度会商工作（郑兆苾等，2000；Shi et al，2001；张国民等，2002；马宏生等，2004；Wu et al，2007）采用PI方法作了一定贡献。Zhang等（2013）对2008年汶川M_S8.0地震和于田M_S7.3地震进行回溯性研究，结果表明对于较长的时间窗以及较大的网格，大地震的预测效果更好。张小涛等（2014）对川滇地区进行的回溯性检验显示2013年芦山M_S7.0地震发震前存在热点信息的演化过程，同时还对该地区进行向前预测，并提出“2022年前川滇菱形地块南部热点地区发生M_L≥6.5地震的概率较高”，实际上2021年的漾濞M_S6.4地震即发生在这一预测区内，显示了该方法的有效性。张盛峰（2015）对PI方法在南北地震带的适用性进行了检验，结果表明PI算法在南北地震带地区的预测效果优于仅考虑丛集性的“相对强度”RI算法和随机预测。Zhang等（2017）对青藏高原及周边地区进行了回溯性研究，其结果表明，PI方法在7—10年的中长期预测上具有更好的预测效能，文中讨论部分提出区域的选择可能会影响预测效能，需后续研究。

2016年1月31日至2022年1月31日在南北地震带发生了7次M_S≥6.0地震（不含余震），为我们研究PI方法的可预测性提供了震例条件。由于PI方法的计算过程涉及所选取区域内所有网格参量的归一化处理，所以不同的选取范围在理论上会产生不同的计算结果。但前人尚未针对选取范围对预测结果的影响进行系统研究，因此，本文将选取不同的地理范围对上述7次M_S≥6.0地震进行回溯性预测研究，并利用R值评分法（许绍燮，1989）和受试者工作特征（receiver operating characteristic，缩写为ROC）检验方法（Swets，1973；Molchan，1997）进行地震预测效能检验，以期获得不同的区域选取对图像信息法可预测性影响的认识。

1. PI方法

PI方法的计算过程具体如下：

1）将研究区划分为N个边长为L的网格，网格中心点X_i代表该网格，对每个网格内的地震事件构建一维向量N_i（t），其中N_i（t）是第i个网格及其最近邻的8个网格$ ［ $（摩尔（Moore）近邻$ ］ $内每单位时间大于截止震级M_C的地震数量。

2）定义五个时刻，t₀，t_b，t₁，t₂，t₃，并基于这五个时刻划分三个时段：背景时段$ ［ { t_{{\rm{b}}}\text{，} {t_{1}}} ］ $，t_b在所使用目录的起始时刻t₀与t₁之间滑动；异常学习时段$ ［ { t_{{{1}}}\text{，} {t_{2}}} ］ $，t₂＞t₁；预测时段$ ［ { t_{{{2}}}\text{，} {t_{3}}} ］ $。

3）定义第i个网格的地震活动强度为I_i（t_b，t），求取网格i和摩尔近邻网格在背景时段$ ［ { t_{{\rm{b}}}\text{，} {t_{1}}} ］ $内震级M＞M_C的地震数量的均值。

$$ {I_i} （ {{t_{\mathrm{b}}}\text{，} t} ）＝ \frac{1}{{t － {t_{\mathrm{b}}}}}\sum\limits_{{t{ \text{′}}} ＝ {t_{\mathrm{b}}}}^t {{N_i} （ {{t{ \text{′}}}} ） } $$

(1)

4）对整个研究区所有网格的地震活动强度${\boldsymbol{I}} $_i（t_b，t）进行标准化，定义为：

$$ \mathop {{{\boldsymbol{I}}_i}}\limits^ \wedge （ {{t_{\mathrm{b}}}\text{，} t} ）＝ \frac{{{{\boldsymbol{I}}_i} （ {{t_{\mathrm{b}}}\text{，} t} ）－ \left\langle {{{\boldsymbol{I}}_i} （ {{t_{\mathrm{b}}}\text{，} t} ） } \right\rangle }}{{\sigma （ {{t_{\mathrm{b}}}\text{，} t} ） }} \text{，} $$

(2)

式中，<I_i（t_b，t）>是所有网格的平均地震活动强度，σ（t_b，t）是所有网格地震活动强度的标准差；

5）定义第i个网格的异常地震活动性为异常学习时段偏离背景时段的差值

$$ \Delta {{\boldsymbol{I}}_i} （ {{t_{\mathrm{b}}}\text{，} {t_1}\text{，} {t_2}} ）＝ \mathop {{{\boldsymbol{I}}_i}}\limits^ \wedge （ {{t_{\mathrm{b}}}\text{，} {t_2}} ）－ \mathop {{{\boldsymbol{I}}_i}}\limits^ \wedge （ {{t_{\mathrm{b}}}\text{，} {t_1}} ） \text{.} $$

(3)

6）计算异常地震活动性的平均变化量，以减少噪声对结果造成的影响。

$$ \overline{\Delta {\boldsymbol{I}}_i （ t_{\mathrm{b}}\text{，} t_1\text{，} t_2 ） }＝\frac{1}{t_1-t_0}\sum\limits_{t_{\text{b}}＝t_0}^{t_1}\Delta {\boldsymbol{I}}_i （ t_{\mathrm{b}}\text{，} t_1\text{，} t_2 ） $$

(4)

7）求取平均异常地震活动性的平方，并将其定义为发震概率函数${\boldsymbol{P}}_{i} （ t_{0}\text{，} t_{1}\text{，} t_{2} ） $：

$$ {\boldsymbol{P}}_i （ t_0\text{，} t_1\text{，} t_2 ）＝\overline{\Delta {\boldsymbol{I}}_i （ t_0\text{，} t_1\text{，} t_2 ） }^2\text{.} $$

(5)

8）计算各网格相对于整个区域的发震概率，将研究区所有网格的平均发震概率作为背景概率：

$$ \Delta {{\boldsymbol{P}}_i} （ {{t_0}\text{，} {t_1}\text{，} {t_2}} ）＝ {{\boldsymbol{P}}_i} （ {{t_0}\text{，} {t_1}\text{，} {t_2}} ）－ \left\langle {{{\boldsymbol{P}}_i} （ {{t_0}\text{，} {t_1}\text{，} {t_2}} ） } \right\rangle \text{，} $$

(6)

热点定义为$\Delta {\boldsymbol{P}}_{i} （ t_{0}\text{，} t_{1}\text{，} t_{2} ） $为正值的区域，即发震概率函数${\boldsymbol{P}}_{i} （ t_{0}\text{，} t_{1}\text{，} t_{2} ） $高于背景概率的区域。

Rundle等（2002）提出了“阈值ω”的概念，ω＝lg（ΔP/ΔP_max），其中ΔP为第i个网格概率与背景概率的差值，ΔP_max为所有网格概率增益的最大值。由式（1）—（6）可知，PI算法以地震目录为输入，对研究区域划分网格后计算出各网格的地震活动强度，后续通过标准化、均值计算以及去除背景概率得到“热点”区域，整个计算过程基于地震活动强度进行，也因此区域地震活动性的不同对PI方法的预测性能会产生影响。

2. 研究区域及计算参数选取

2.1 研究区域

本文所选的研究区域为大陆中部地区（21.0°N—41.5°N，96.0°E—110.0°E），涉及南北地震带及其周边区域（图1）。虽然对于南北地震带展布的范围认识不一（康来迅，1991；杨兴悦等，2012；张盛峰，2015），但都包含于本文所选的研究区域。图1中黄色边框矩形是张盛峰（2015）提到的南北地震带（21°N—41.5°N，97.5°E—107.5°E）范围，包括我国内蒙古、四川、重庆、云南、贵州、宁夏、甘肃、陕西等多个省市。根据发震构造不同，南北地震带分为南段（28°N以南）、北段（34°N以北）以及中段（28°N—34°N）三段（杨兴悦等，2012）。南北地震带历史上共记录到M_S≥7.0地震67次，包括7次M_S≥8.0地震，其中北段4次，中段2次，南段1次，这些大震对我国的经济、民生都造成了巨大的损失。从1970年以来的资料来看，该区域共发生M_S6.0—6.9地震73次，M_S7.0—7.9地震18次，M_S≥8.0地震1次（2008年汶川M_S8.0地震），即平均一年发生近2次M_S≥6.0地震，3年发生1次M_S≥7.0地震，50年1次M_S≥8.0地震。

图 1 本文研究区域和相关目标地震图示

Figure 1. The regions designed and related target earthquakes in this study

下载: 全尺寸图片幻灯片

为了获得因区域选取不同对图像信息法可预测性的影响，本文设计了6个不同的研究区域，内含全部或部分2016年以来发生在南北地震带的M_S≥6.0地震，通过ROC检验方法评估区域的选取对图像信息法地震预测效能的影响。研究区域列于表1，并示于图1。

表 1 本文研究区域

Table 1. Regions designed in this study

研究区域序号	经纬度范围	网格（1°×1°）数
1	（21°N—41.5°N，97.5°E—107.5°E）	205
2	（26°N—41.5°N，97.5°E—107.5°E）	155
3	（21°N—41.5°N，100°E—110°E）	205
4	（26°N—41.5°N，100°E—110°E）	155
5	（21°N—41.5°N，96°E—104°E）	164
6	（30°N—41.5°N，96°E—104°E）	92

下载: 导出CSV

| 显示表格

设计的6个研究区分别编号1—6，其中：研究区1范围为南北地震带区域，自北向南依次包括阿拉善、祁连、鄂尔多斯、柴达木、巴颜喀拉、川滇、华南、滇南等8个二级地块；研究区2相对研究区1去除了滇南和部分华南、川滇地块区域；研究区3自北向南依次包括阿拉善、燕山、祁连、鄂尔多斯、柴达木、巴颜喀拉、川滇、华南、滇南等9个二级地块，相对研究区1去除了阿拉善、巴颜喀拉、川滇地块的部分地区，新加了鄂尔多斯、华南地块的部分区域；研究区4相对研究区3去除了滇南地块、华南、川滇地块部分地区；研究区5自北向南依次包括阿拉善、祁连、柴达木、巴颜喀拉、川滇、华南、滇西、滇南等8个二级地块，相对研究区1和3去除了鄂尔多斯地块和华南地块部分地区，新增了滇西地块，扩大了位于阿拉善、祁连、柴达木、巴颜喀拉、川滇地块的区域；研究区6相对研究区5去除了川滇、华南、滇西、滇南地块以及巴颜喀拉的部分区域。

总的来看，研究区域按照1与2，3与4，5与6分为三组，其中研究区2，4，6与1，3，5在东西向上具有相同的经度范围，但在南北向上不同程度地去除了南部部分地区。这样设计便于对比区域变化对预测效能的影响。

2.2 用于回溯性PI预测检验的目标地震

前人曾利用PI方法对南北地震带及周边区域进行强震回溯性检验预测研究（蒋长胜，吴忠良，2008；蒋长胜等，2009；Zhao et al，2010；Zhang et al，2013，2017；张小涛等，2014；张盛峰，2015），针对2008年汶川M_S8.0、2013年芦山M_S7.0、2014年鲁甸M_S6.5等多次地震都取得了较好的回溯性检验预测效果。2016年以来南北地震带地震频发，因此选取2016年1月31日至2022年1月31日期间的M_S≥6.0地震作为回溯性检验研究的震例。

本文研究区域中，2016年1月31日至2022年1月31日共发生7次M_S≥6.0地震，其中2019年四川长宁M_S6.0地震和2021年四川泸县M_S6.0地震的发震构造被认为与工业开采有关（易桂喜等，2019，2021），但无论是天然构造地震，还是水库诱发地震、注水诱发地震等人工诱发地震，地震的发生归根结底是沿某一构造面或结构面发生破裂所致（常祖峰等，2020）。因此，本文将上述两震例纳入目标地震，并通过PI图像的演化过程观察天然构造地震与人工活动所引起的地震PI图像演化过程是否存在差异。所选目标地震的空间分布见图1，具体信息见表2。

表 2 本文研究的震例

Table 2. Earthquake cases in this study

序号	发震时间（北京时间）年-月-日	震中	北纬/°	东经/°	M_S
1	2 017−08−08	四川九寨沟	33.20	103.82	7.0
2	2 019−06−17	四川长宁	28.34	104.90	6.0
3	2 021−05−21	云南漾濞	25.67	99.87	6.4
4	2 021−05−22	青海玛多	34.59	98.34	7.4
5	2 021−09−16	四川泸县	29.20	105.34	6.0
6	2 021−12−24	老挝	22.33	101.69	6.0
7	2 022−01−08	青海门源	37.77	101.26	6.9

下载: 导出CSV

| 显示表格

2.3 PI计算参数的设定

由PI方法原理可知，其算法运行需要进行截止震级M_C、目标震级M_T、地震事件深度、网格尺度$\Delta $x、滑动步长$\Delta $t以及三个时间段（t_b−t₁， t₁−t_2，t₂−t₃）等参数的设定。

对于截止震级M_C的设定，首先需要考虑目录的完备性。本文使用的目录由中国地震台网中心提供。在1970年1月1日至2022年1月31日期间，研究区域共记录到M_S3.0—3.9地震5 763次，M_S4.0—4.9地震1153次，M_S5.0—5.9地震357次，M_S6.0—6.9地震145次，M_S7.0—7.9地震26次，M_S8.0—8.9地震3次。利用地震目录绘制G-R关系图，见图2。该图显示该地区的完备震级为M_S2.0。Holliday等（2005）的研究表明，PI方法的预测目标震级M_T≥M_C＋2.0，因此本文使用的截止震级M_C 和目标震级M_T分别为4.0和6.0。根据前人经验，热点显示的阈值参数取lg（ΔP/ΔP_max）＝−0.6 （Zhang et al，2013）。

图 2 研究区地震目录完整性分析

Figure 2. Analysis of the seismic catalogue completeness in study area

下载: 全尺寸图片幻灯片

对于网格大小的选取，Zhang等（2017）对青藏高原及周边地区的研究表明，预测时间窗为8—10年时，网格为1°×1° 效果更佳，预测时间窗为7—10年时网格为2°×2°效果更佳。张小涛等（2014）对川滇地区9年预测时间窗的回溯性研究表明，2°×2°具有更好的预测效能。夏彩韵等（2015）对新疆两次于田M_S7.3地震进行回溯性研究发现，1°×1°在该地区的预测效能更好。由于本研究所选择的部分区域较小，网格划分为2°×2°时网格数偏少，会对结果稳定性有一定影响，因此，本文采用1°×1°的网格划分，并系统地研究选取不同区域对图像信息方法可预测性的影响。PI算法适用5—10年的中长期预测，针对M_S≥7.0地震，较长的时间窗会取得更好的预测效果（Zhang et al，2017），本文研究区域主要是M_S≥6.0地震，因此选取5年时间窗、1年滑动步长进行回溯性检验。三个时段的选择详列于表3。

表 3 本研究的时间窗选取

Table 3. Time windows designed in this study

序号	t₀	t₁	t₂	t₃
1	19 70	20 11−01−31	20 16−01−31	20 21−01−31
2	19 70	20 12−01−31	20 17−01−31	20 22−01−31
3	19 70	20 13−01−31	20 18−01−31	20 23−01−31
4	19 70	20 14−01−31	20 19−01−31	20 24−01−31
5	19 70	20 15−01−31	20 20−01−31	20 25−01−31
6	19 70	20 16−01−31	20 21−01−31	20 26−01−31
7	19 70	20 17−01−31	20 22−01−31	20 27−01−31
注：表中t₀，t₁，t₂和t₃为PI算法的时间点，其中背景时段为$ ［ t_{{\mathrm{b}}} \text{，} t_{1} ］ $，异常学习时段为$ ［ t_{{{1}}} \text{，} t_{2} ］ $，预测时段为$ ［ t_{{{2}}} \text{，} t_{3} ］ $，详见第一节。共得到7个不同的时间段用于PI计算。

下载: 导出CSV

| 显示表格

3. 计算结果及分析

3.1 区域选取不同对PI预测效能的影响

3.1.1 计算结果分析

根据设定的PI参数以及六个研究区进行PI计算，得到6 （区域数）×7 （预测时段数）共42幅PI热点图，这里仅对部分图幅进行分析。

由于所计算的7个预测时段中仅时段1和2所覆盖的预测时间已经结束，可以进行统计检验，因此下文基于预测时段2 （2017−01−31—2022−01−31）对6个研究区的结果进行分析。图3为6个研究区域在2017年1月31日至2022年1月31日的PI预测图以及在相同时段内研究区实际发生的M_S≥6.0地震（蓝色圆圈）分布。图3a中只有南北地震带南部的滇西南地区出现热点，2021年12月24日老挝M_S6.0地震所在网格的西北摩尔邻近网格出现热点，整个区域的北部在预测时段发生的3次M_S≥6.9地震却没有热点出现。

图 3 6个研究区在2017年1月31日至2022年1月31日预测时段的PI热点图

蓝色圆圈表示预测时段内发生的地震，灰色倒三角表示异常学习时段内发生的地震（a）研究区1；（b）研究区2；（c）研究区3；（d）研究区4；（e）研究区5；（f）研究区6

Figure 3. PI seismic hotspots maps in the period from 31 January 2017 to 31 January 2022 for six regions

Blue circles indicate earthquakes that occurred during the prediction period，and gray inverted triangles represent earthquakes that occurred during the anomaly learning period（a） Region 1；（b） Region 2；（c） Region 3；（d） Region 4；（e） Region 5；（f） Region 6

下载: 全尺寸图片幻灯片

区域2 （图3b）是在区域1 （图3a）的基础上去除了滇南和部分华南、川滇地块区域。这个区域为南北地震带的中段、北段。在图3a中，南北地震带北段没有热点。在其它计算参数相同的情况下，去除了图3a南部部分区域后的图3b研究区就出现了较多PI热点。2019年6月17日长宁M_S6.0地震、2021年9月16日泸县M_S6.0地震以及2022年1月8日门源M_S6.9地震都位于新增热点附近，但M_S≥7.0地震均未报准。由图3a与图3b的对比结果可以看出，区域选择不同会造成PI热点图的不同。

区域3 （图3c）去除了区域1 （图3a）中的阿拉善、巴颜喀拉、川滇地块的部分地区，新加了鄂尔多斯、华南地块的部分区域。由图3c可见：该研究区域在中部、南部出现了较多的PI热点，且2021年12月24日老挝M_S6.0地震和2021年9月16日泸县M_S6.0地震均落在PI热点内，2019年6月17日长宁M_S6.0地震所在网格的摩尔邻近网格也有热点，算报准地震。

区域4 （图3d）是把区域2 （图3b）的研究区边界整体向东平移了2.5°。由图3d可以看到热点的空间分布图有部分是相同的，但热点颜色（发震概率越高的地方越深）有变化，例如：图3d中发震概率最高的地方在南部的西昌至盐源一带，但图3b中发震概率最高的地方在甘肃祁连到青海门源一带。与图3b相比，图3d的热点减少并相对集中。从地震的对应情况看，图3b可报准3次地震，图3d可报对2次地震。

区域5 （图3e）的研究区域是把区域1 （图3a）研究区域的左边界向西平移1.5°，右边界向西平移3.5°。两个区域的PI热点在滇西南地区大致相同，但图3e中川滇交界和云南中部出现了新的PI热点。2021年5月21日的漾濞M_S6.4地震在图3e的研究区域可被预测，但在图3a中被漏报。

区域6 （图3f）是在区域5 （图3e）基础上去除了川滇、华南、滇西南地块以及巴颜喀拉部分区域。在区域5中，30°N以北地区（图3e）并没有热点，但研究区在图3e的基础上仅保留30°N以北地区后，PI热点就出现了，主要集中在四川松潘、理县一带和青海德令哈至都兰一带，以及青海祁连至甘肃山丹一带。且除了2022年1月8日门源M_S6.9地震报准，其它热点所在网格及其摩尔邻近网格均未发生M_S≥6.0地震，属于虚报热点。

上述结果表明，在其它参数一致但研究区不同的情况下，PI预测结果会有很大的差异。这主要是因为研究区域的变化会导致地震活动性的不同，若将地震活动性高的区域和地震活动性低的区域一起计算，会使得后者的异常很难识别。

表4是图3中各研究区域内发生在预测时间窗内M_S≥6.0地震的预测检验情况。可见，不同区域对同一次地震的预测情况是不同的。例如，在包含泸县M_S6.0地震的4个研究区域内，研究区1和4漏报，而研究区2和3报准。

表 4 2017年1月31日至2022年1月31日预测时段不同区域的预测结果与地震的对应关系

Table 4. The corresponding relationship between the earthquake prediction results and earthquakes in different regions during the prediction period from January 31，2017 to January 31，2022

研究区	不同区域的预测与地震的对应关系						各研究区的地震报准率
研究区	四川九寨沟 M_S7.0	四川长宁 M_S6.0	云南漾濞 M_S6.4	青海玛多 M_S7.4	四川泸县 M_S6.0	老挝 M_S6.0	各研究区的地震报准率
1	0	0	0	0	0	1	0.17
2	0	1	/	0	1	/	0.50
3	0	1	/	/	1	1	0.75
4	0	1	/	/	0	/	0.33
5	0	/	1	0	/	0	0.25
6	1	/	/	0	/	/	0.50
各地震的报准比例	0.16	0.75	0.50	0	0.50	0.67
注：表中“1”表示报准地震，“0”表示漏报地震，“/”表示该地震未落在该区域的预测时段。

下载: 导出CSV

| 显示表格

对于不同的研究区域，预测时间窗内发生的地震报准比例也不同。由表4可见，研究区3能报准75%的地震，研究区2和6能报准50%的地震，研究区4能报准33%的地震，研究区5能报准25%的地震，而研究区1的地震报准率最低，只有17%。

从表4还看到：6次M_S≥6.0地震中，四川九寨沟M_S7.0地震和青海玛多M_S7.4地震在这一时间窗的预测效果较差；对于九寨沟M_S7.0地震，除了研究区6能预测外，其它5个研究区都漏报；对于玛多M_S7.4地震，包含该地震的4个研究区在这一时间窗均漏报；M_S6.0地震在各研究区都有较好的对应效果。这很可能与本研究设置的PI参数有关，这是下一步需要研究的问题。

3.1.2 R值评分和ROC检验

地震预测方法的预测效能需要经过统计检验，本文参考前人做法（Holliday et al，2005），采用R值评分法和ROC检验进行效能评估。

1） R值评分检验。R值评分法是由许绍燮（1989）提出，本文采用Shi等（2001）的算法。其主要步骤是：将研究区域网格化，计算报准率减去预报网格数与总网格数之比。随机预测时R＝0，当R大于0时表示该方法的预测效能高于随机预测。

$$ R＝\frac{有震报准数}{总震数}－\frac{预报网格数}{总网格数} $$

(7)

将图3中的6个研究区域在2017年1月31日至2022年1月31日预测时间窗的R值计算出来，列于表5。可以看出：所有研究区的R值评分都大于0，这表明PI方法对M_S≥6.0地震预测总体上是有效的；研究区3的预测效能最高，R值评分达到0.71；研究区6的预测效能次之，R值评分为0.60。从图1可见，这两个区域的共同特点是去除了地震活动性偏高的云南及附近地区，这样就不会出现区域内地震活动差异导致高活动性地区孕震信息压制低活动性地区孕震信息的问题。

表 5 各区域统计检验结果

Table 5. Statistical test results of each region

研究区	R值评分	有效预测系数E_f
1	0.13	0.377 3
2	0.55	0.463 9
3	0.71	0.403 5
4	0.46	0.420 7
5	0.13	0.331 0
6	0.60	0.287 0

下载: 导出CSV

| 显示表格

2） ROC检验。ROC检验主要是通过计算不同阈值ω下的获得的“命中率”（纵轴）和“虚报率”（横轴）并绘制ROC曲线的评估方法（Swets，1973；Molchan，1997）。Zhang等（2013）将ROC检验图的预测检验线与随机预测检验线所包络的面积定义为有效预测系数E_f，E_f值介于0与0.5之间，E_f值越大，预测效果越好。预测时段为2017年1月31日至2022年1月31日PI算法的ROC曲线如图4所示，E_f值列于表5。

图 4 6个区域的PI算法回溯性研究的ROC检验

（a）研究区1；（b）研究区2；（c）研究区3；（d）研究区4；（e）研究区5；（f）研究区6

Figure 4. ROC test of the PI forecasting for the six selected regions

（a） Region 1；（b） Region 2；（c） Region 3；（d） Region 4；（e） Region 5；（f） Region 6

下载: 全尺寸图片幻灯片

由表5和图4可得，PI方法在所有区域的预测效能都高于均匀泊松模型（随机预测），其中研究区2的E_f值最高，接近最大值，之后依次为研究区4，3，1，5，6。研究区2，3，4的R值检验和ROC检验得分较高，这显示以上研究区域用于向前预测时能得到较好效果。研究区6的R值评分与ROC检验得分差异较大，表明该区域向前预测具有不稳定的特征，可能是由于参数等原因，后续可深入研究。研究区域1和5的预测效果在6个区域中排倒数，不适于向前预测。

3.2 PI热点图的时空演化趋势

图3仅显示了单一预测时段的结果，实际上可以通过多个预测时段的结果来获取热点异常的演化趋势，继而进行地震预测。下面将以三个典型区域为例来分析讨论。第一个为研究区1；第二个为研究区3，从表5给出的结果看该区和研究区2在图3预测窗内的预测效能最优；第三个为研究区6，从表5给出的结果看该区在图3预测窗内的预测效能稳定性最差。

3.2.1 研究区1的PI热点演化与强震活动分析

图5为研究区1在7个预测时间窗的PI预测图像。

图 5 研究区1内不同预测时间窗口的PI地震热点演化

蓝色圆圈表示预测时段内发生的地震，灰色倒三角表示异常学习时段内发生的地震

Figure 5. Evolution of PI seismic hotspots for different time windows in region 1

Blue circles indicate earthquakes that occurred during the prediction period，and gray invertedtriangles represent earthquakes that occurred during the anomaly learning period（a） 2016—2020；（b） 2017—2021；（c） 2018—2022；（d） 2019—2023；（e） 2020—2024；（f） 2021—2025；（g） 2022—2026

下载: 全尺寸图片幻灯片

图5b即图3a。在图3a预测时间窗时，2021年漾濞M_S6.4地震、2021年泸县M_S6.0地震震中附近并未出现“热点”，因此在表5中评估预测效能时，为漏报地震。但在图5d，5e和图5f这三个时间窗内，这两次地震震中附近就出现了PI热点，并且随时间发震概率增高，显示这两次地震可以在这三个时间窗口报准。

但2021年玛多M_S7.4地震震中附近无论在哪一个时间窗口的PI热点预测图上均未出现热点，不能被预测出来。这个研究区的PI时空演化图显示南北地震带南段的M_S≥6.0地震可以在震前作出预测。但北段基本未出现PI热点，发生在北段的2017年九寨沟M_S7.0地震、2021年玛多M_S7.4地震和2022年M_S门源6.9地震都未能作出预测。本文参考前人（Holliday et al，2005；Zhang et al，2013）的研究结果，所选的震级下限为M4.0，所预测的目标震级为M6.0以上，预测时间窗为5年，网格尺寸为1°×1°。从检验结果看，此套参数对M_S6.0地震的预测效果较好，但对M_S7.0地震的预测效果较差，这表明对于M7.0以上地震的预测还需要优化参数。同时况春利等（2019）研究显示，二十世纪以来南北地震带中南段的中强震显著多于北段，认为中段和南段的地震分布较为相似，与北段存在明显差异。因此若同时分析南北地震带中、南段和北段，则北部相对较弱的地震异常信息会被掩盖。本研究认为是以上原因导致了北段的地震预测效果较差。

3.2.2 研究区3的PI热点演化与强震活动分析

图6为研究区3在7个预测时间窗的PI预测图像。

图 6 研究区3不同预测时间窗口的PI地震热点演化图

蓝色圆圈表示预测时段内发生的地震，灰色倒三角表示异常学习时段内发生的地震

Figure 6. Evolution of PI seismic hotspots for different time windows in region 3

Blue circles indicate earthquakes that occurred during the prediction period，and gray inverted triangles represent earthquakes that occurred during the anomaly learning period（a） 2016—2020；（b） 2017—2021；（c） 2018—2022；（d） 2019—2023；（e） 2020—2024；（f）2021—2025；（g）2022—2026

下载: 全尺寸图片幻灯片

图6a即图3c。在图3c时间窗口中，2021年长宁M_S6.0、2021泸县M_S6.0以及2021老挝M_S6.0地震的震中及其摩尔近邻网格都出现了“热点”，但热点对应的发震概率较低。在图6d，6e和图6f三个时间窗中，上述三次地震震中对应的PI 热点值随着时间推移逐渐升高，同时PI热点对应的网格数减少且更加集中于地震震中位置。

2022年门源M_S6.9地震震中在图6b-f这5个时间窗中都位于PI热点内，且PI热点呈现一个“出现→达到峰值→逐渐降低”的演化过程，该地震的预测效果很好。

2017年九寨沟M_S7.0地震只存在于图6a和6b的时间窗，且该地震震中都无PI热点产生，现有结果不能预测该地震。针对本文设计的5年预测时间窗，若深入研究需给出剩下包含该地震的预测时间窗的PI预测图像。这个研究区的PI时空演化图显示该区域的M_S6.0地震都可以在震前作出预测。

3.2.3 研究区6的PI热点演化与强震活动分析

图7为研究区6在7个预测时间窗的PI预测图像。

图 7 研究区6不同预测时间窗口的PI地震热点演化图

蓝色圆圈表示预测时段内发生的地震，灰色倒三角表示异常学习时段内发生的地震

Figure 7. Evolution of PI seismic hotspots for different time windows in region 6

Blue circles indicate earthquakes that occurred during the prediction period，and gray inverted triangles represent earthquakes that occurred during the anomaly learning period （a） 2016—2020；（b） 2017—2021；（c） 2018—2022；（d） 2019—2023；（e） 2020—2024；（f） 2021—2025；（g） 2022—2026

下载: 全尺寸图片幻灯片

图7b即图3f。在图3f时间窗口中，2021年玛多M_S7.4 地震震中附近并未出现热点，因此在表5中评估预测效能时，为漏报地震。且在后续的时间窗口中，该地震周围始终未出现热点，不能被预测出来。但图7b-7f中2022年门源M_S6.9地震却始终位于PI热点及其摩尔邻近网格内，说明该区域对于M_S≥6.0地震可以在震前作出预测。

3.3 研究区综合分析

对计算所得的43幅PI热点图进行分析，得出以下认识：

1）在2017年1月31日至2022年1月31日预测时段，除去研究区6对2017年M_S7.0九寨沟地震有一定的预测效能，其它区域针对M_S≥7.0的地震预测效能较差。可能是因为本文的PI参数适合M_S6.0地震的预测，后续可根据Zhang等（2017）的研究结果探究PI方法针对M_S≥7.0的参数设定。

2）一般情况下，强震发生后震中附近的热点会逐渐消失，震后震中附近如果持续存在热点，后续还可能发生强震。例如：2016年1月21日门源发生M_S6.4地震，但研究区域1—4显示地震发生后，该地区的地震热点始终存在，之后发生了2022年门源M_S6.9地震。

3）对各研究区的热点演化过程的分析可知：① 在研究区1中南部滇西南地区始终存在热点且在2021年漾濞M_S6.4地震、2021年老挝M_S6.0地震发生后，高概率热点依然存在；② 研究区2中海原断裂中东部附近存在逐渐演化的高概率热点；③ 在研究区3和4中可以看到2021年泸县M_S6.0地震发生前有一个明显的热点概率逐渐增加达到峰值再衰减的演化过程，但是地震发生后龙门山断裂南部和小江断裂东北部地区的地震热点对应概率反而达到了一个极高值。此外，研究区5中小江断裂处同样出现了高概率热点； ④ 长宁M_S6.0地震震中附近热点趋势呈现“出现→概率降低或消失→震前再现”的演化趋势，泸县地震震中附近热点演化在研究区1和3却呈现两种不同的趋势，在研究区1中为“逐渐出现→概率升高→逐渐消失”的演化趋势，在研究区3中呈现“逐渐出现→消失→再次出现→概率升高→震后热点持续存在”的演化趋势。两次地震均出现了热点反复“出现→消失”的过程，相比于门源、漾濞等天然构造地震呈现的“震前逐渐出现→达到峰值→逐渐消失”的热点演化过程具有一定的差异性。

根据前文PI热点演化与强震之间的关系，上述地区为M_S≥6.0地震值得关注的发震区域。

4. 讨论与结论

本文以中国大陆中部地区（南北地震带及周边区域）为研究范围，通过划分6个不同的研究区，开展了不同空间范围的选取对地震预测效果的影响研究，并基于ROC和R值评分方法进行了统计检验，得到的结论如下：

1）相比于地震活动强度差异较大的区域，图像信息方法在地震活动强度存在差异的区域具有更好的预测效能。这可能是导致PI方法对南北地震带南部的地震具有较好预测效果，但是对北部地区的地震预测效果则较差的原因。

2）根据对所有研究区的“热点”时空演化趋势分析，2027年以前滇西南地区、海原断裂中东部附近、小江断裂中部地区以及龙门山断裂南部和小江断裂东北部地区存在的“热点”地区发生M_S≥6.0地震的概率较高。

PI方法在实际的地震预报过程中仍然存在许多局限性：

1）进行回溯性研究时，往往会通过参数调整以突出震前热点，具有较强的主观性，在实际的预测过程中需提前分析选定研究区域的最优参数设置以确保PI方法预测的稳定性。根据Zhang等（2013）的结果，本文选定的时间窗尺度为5年，后续可深入研究不同时间窗尺度的PI方法的预测效能，以将该方法用于年度危险区判定和地震大形势预测。

2）在本文的研究结果中，所有的研究区域对于2021年玛多M_S7.4地震都未呈现较好的预测效果，推测原因有以下几点：① 震前没有出现偏离平衡态的异常地震活动是可能的，玛多地震可能属于此类地震，需后续综合各因素深入研究；② 玛多地震的震中位于青藏地块内次级地块巴颜喀拉地块内部，本文主要针对南北地震带整体进行研究，后续可以针对巴颜喀拉地块对玛多地震进行深入研究；③ 本文所采用的截止震级和目标震级为M_C4.0和M_T6.0，针对M_S7.0地震，后续可选定合适的截止震级进行研究；④ 针对大地震，相对大的网格和较长的时间窗会取得更好的效果（Zhang et al，2013）。本文主要采用的是1°×1°的网格以及5年的时间窗，后续可选择更大的网格以及更长的时间窗进行研究。

3）在目标震例中，本文选取了被认为与工业开采有关的两次地震，分别为2019年长宁M_S6.0地震和2021年泸县M_S6.0地震。在PI图像演化的分析过程中，两地震震中附近都出现了热点反复出现→消失的现象，相比于门源、漾濞等天然构造地震的热点演化过程具有一定的差异性。推测认为这两次地震的PI图像演化区别于构造地震，可能主要受人工开采活动的影响，但还需后续对相关地区活动深入调查分析。

本文中还有一个值得探讨的现象：在2016年的门源M_S6.0地震发生以后，该地区的热点并未消失，而是始终处于一个发震概率的高值，之后2022门源M_S6.9地震发生。夏彩韵等（2015）针对两次新疆于田M_S7.3地震进行了检验，其结果表明它们都发生在发震概率的峰值处。后续可深入研究两次门源地震的关联性。

中国地震台网中心提供了地震目录数据，中国地震局地球物理研究所张琰博士在作图方面提供了帮助，作者在此谨表谢意。

[1] L. Knopoff, Analytical calculation of the fault—plane problem, publ. Domin. Obs, 24, 309—315, 1961.

[2] A. , J. Wickens and J. H. Hodgson, Computer re—evaluation of earthquake mechanism, Publ.Domin. Obs. 33, 1, 1—10, 1967.

[3] Hoang—Trong Pho, Extended distances and angles of incidence of P waves, Bull. Seism. Soc. Am., 62, 4, 885—902, 1972.

[4] A. Ben—Menahem, Radiation of seismic surface wave from finite moving sources, Bull. Seis. Soc Am., 51, 401, 1961.

[5] N. A. Haskell, Total energy and energy spectral density of elastic wave radiation from pro—pagation faults, Bull. Seism. Soc. Am., 54, 6, 1811—1841, 1964.

[6] T. Hirasawa and W. Stauder, On the seismic body wave from a finite moving source, Bull Seism. Soc. Am., 55, 237—262, 1965.

[7] J. C. Savage, Relation of corner frequency to fault dimensions, J. Geophys. Res., 77, 20, 3788——3795, 1972.

[8] T. C. Hanks and M. Wyss, The use of body wave spectra in the determination of seismic source parameters, Bull. Seis. Soc. Am., 62, 561—589, 1972.

[9] H. Berekhemer and K. H. Jacob, Investigation of the dymayical process in earthquake foci by analyzing the pulse shape of body wave, AD—738953, 1968

[10] 中国科学院地球物理研究所, 地震仪器概论, 科学出版社, 90—104, 1975.

[11] 虢顺民等, 唐山地震区域构造背景和发震模式的讨论, 地质科学, 4, 305—320, 1977.

[12] J. N. Brune, Toctonic stress and the spectra of seismic shear wave from earthquakes. J. Geophys. Res., 75, 26, 4997—5009, 1970.

[13] M. A. Chinnery, Theoretical fault models, Publ. Domin. Obs., 37, 7, 221—223, 1969.

[14] G. S. Stewart, R. Butler, and H. Kanamori, Surface and body wave analysis for the Feb. 4, 1975 Haicheng and July 27, 1976 Tangshan Chinese earthquakes (abstract) EOS, Trans. Amer. Geophys. Union, 57, 12, 953, 1976.

[15] 陈运泰等, 巧家、石棉的小震震掠参数的侧定及其地震危险性的估计, 地球物理学报, 19, 3, 206—231, 1976.

[16] G. C. Sih and H. Liebowitz, Mathematical theories of brittle fracture, Fracture, 11, 68—188, 1968.

[17] G, C. Sih, Mechanics of fracture 1, Methods of analysis and solutions of crack problems, 1973.

[18] О.Г.Шамина, А.А0Павлв, Ю.ф.Копничев, Исследование Процесса подговки Грешены, Изв.АН.СССР.фчз.Зем., 8, 17—30, 1973.

[19] E. Hoek and Z. T. Bieniawski, Brittle fracture propagation in rock under compression. International J. of Fracture Mech., 3, 137—155, 1965.

[20] 傅承义, 地球十讲, 科学出版社, 1976.

[21] 郭增建等, 震源孕育模式的初步讨论, 地球物理学报, 16, 43—48, 1973.

期刊类型引用(2)

1.	TIAN Weixi，ZHANG Yongxian，FENG Maoning，JU Changhui，ZHANG Shengfeng. CSEP-CN: Parameter Optimization of Pattern Informatics Method. Acta Geologica Sinica(English Edition). 2024(S1): 67-69 . 必应学术
2.	LI Baochun，YE Gaofeng，JIN Sheng，WEI Wenbo，ZHANG Letian. Rheological Evidence of the Lithospheric Destruction of the Eastern Block of the North China Craton. Acta Geologica Sinica(English Edition). 2024(S1): 36-38 . 必应学术

唐山地震的破裂过程及其力学分析

计量

THE FRACTURE PROCESSES OF THE TANGSHAN EARTHQUAKE AND ITS MECHANICAL ANALYSIS