0 引言

2020年7月12日6时38分，在河北省唐山市古冶区发生M_S 5.1地震，震源深度10 km，震中位于（39.78°N，118.44°E）。此次地震发生于NE向唐山断裂带中唐山—古冶断裂北侧，巍山—长山断裂延伸处。唐山—古冶断裂是全新世走滑活动断层，巍山—长山断裂走向NE、倾向NW，是一条逆断型断裂（朱琳等，2021）。2021年4月16日，在此次古冶M_S 5.1地震震中东侧约30 km处发生滦县M_S 4.3地震，震中位于（39.75°N，118.71°E），震源深度9 km。2次地震均位于唐山M_S 7.8地震附近，震中距离较近且发生间隔时间不到1年，因此古冶M_S 5.1地震对滦县M_S 4.3地震是否存在触发作用值得研究。

由于此次古冶地震发生在1976年唐山M_S 7.8地震震源区，唐山区域地震活动性得到诸多学者的密切关注，相关研究有：徐志国等（2021）根据台网观测报告，利用双差定位的方法进行重定位，并结合波形资料采用近震全波形方法得到此次地震的震源机制解，分析得到此次地震为走滑型，发震断层面为SW向节面，结合1976年唐山地震发震构造，认为此次地震应为唐山主震的一次余震；朱琳等（2021）据1976年唐山强震群震源位错模型，采用Burgers流变模型，计算唐山强震群在2020年古冶M_S 5.1地震震源区产生的库仑应力变化，结果显示，震后15年前后，库仑应力演化状态差异显著；张苏祥等（2022）基于改进的DBSCAN算法自动识别断层的方法对唐山地区的断裂进行了识别；唐杰等（2023）综合构造背景、多年小震时空图像背景、多次地震的震源机制解结果等资料分析表明：此次地震发震断裂为与唐山断裂带共轭的抹轴峪断裂，未来该区域仍具有一定发震风险；余海琳（2023）基于密集台站观测的数据，使用双差定位的方法，给出古冶地震序列的精定位结果，密集台站观测的数据为万永革科研团队在古冶M_S 5.1地震发生后，基于其团队在古冶地震震源区布设的短周期地震台、其他研究团队布设的长周期地震台，以及固定台波形数据研究所得；严少鹏等（2024）利用震后布设的密集台阵观测资料和区域固定台网观测资料，识别出20 505个P波初动符号，得到96个震源机制解和该地区0.25°×0.25°的精细地壳应力场，结果表明，唐山古冶地震震源区小震震源机制以走滑型为主，其次是正断型；P轴走向为NNE向，多数倾伏角偏大，T轴走向为NNW向，倾伏角偏小。

地震发震构造及应力触发是地震学研究的重要课题之一。在Lee等（2024）的最新研究中，强调了断层的几何形状和粗糙度对断层滑移的重要性，并且表明断层网络中的大规模复杂性在断层破裂过程中具有重要作用。

本文使用余海琳（2023）给出的唐山古冶M_S 5.1地震及其余震的精定位结果。基于改进的模糊聚类方法（万永革等，2023），给出古冶M_S 5.1地震丛集的断层几何形状及其误差范围；并基于此次地震震源区的震源机制数据（王晓山等，2020），使用网格搜索法（Wan et al，2016）反演构造应力场，进一步计算发震断层滑动角。最后，基于本文研究成果（古冶地震发震断层的走向、倾角和滑动角）以及余海琳（2023）给出的滦县M_S 4.3地震的走向、倾角和滑动角等信息，计算古冶M_S 5.1地震对滦县M_S 4.3地震以及其他余震的应力触发作用。这些研究结果对了解此次古冶地震的发震断层构造、该震源区构造应力场以及主震对余震的应力触发提供了基础资料，具有一定科学意义。

1 研究方法

本文采用的方法是基于GK改进的模糊聚类算法（Gustafson et al，1979），并采用自适应距离对模糊算法进行了优化（万永革等，2023）。该方法在门源地震序列发震断层、九寨沟地震序列断层面几何形状和芦山地震序列断层面几何形状和滑动性质的求解中应用效果良好（万永革等，2023；Jin et al，2023；靳志同等，2024）。

该方法首先给出聚类个数、猜测的初始类中心和每个震源点属于每个类的概率（大于50%）以及事件之间的距离，该距离为使用空间两点间的距离公式进行计算使用欧式距离计算；其次，采用Gustafson等（1979）提出的改进模糊聚类分析方法，当属于某个类的所有震源点的距离之和最小时，即可确定某类的中心位置和属于每一类的概率，初步得到断层面上地震的分类情况；然后，计算震源点协方差矩阵的特征值，若最大特征值和最小特征值之差较大，则该类中存在一个子断层面，并求出该断层面的走向和倾角，以及断层面的4个顶点。若最小特征值和其余2个特征值差别不大，则该类震源点内部无明显断层面；最后，以上一步给出的断层面走向和倾角作为初始值，采用万永革等（2008）方法中的高斯—牛顿法，对断层面走向和倾角进行精细求解，最终确定断层面相关参数的误差和4个顶点的位置。采用该方法能够直接清晰呈现该区域多个断层面的三维分布以及断层面参数信息。

GK模糊聚类的目标函数为

$\min J_m(\boldsymbol{X}, \boldsymbol{U}, \boldsymbol{V}, \boldsymbol{A})=\sum\limits_{i=1}^c \sum\limits_{k=1}^n\left(\mu_{i k}\right)^m D_{i k d_i}^2$

(1)

式中，$\boldsymbol{X}=\left\{x_1, x_2, \cdots, x_n\right\}$为由n个数据点$X_k=\left\{x_{k 1}, x_{k 2}, \cdots, x_{k p}\right\} \in R^p(k=1, 2, \cdots, n)$构成的数据集，i为类的下标$(i=1, 2, \cdots, c)$，k为数据样本值下标$(k=1, 2, \cdots, c)$；$\boldsymbol{U}=\left(\mu_{i k}\right)_{c \times n}$为隶属度矩阵，μik表示第k个元素与第i个类的隶属度；$\boldsymbol{V}=\left\{v_1, v_2, \cdots, v_c\right\}$，$v_i \in R^p$为第i个类的中心$(i=1, 2, \cdots, n)$；$m \in[1, \infty)$为模糊指数，常令1.5≤m≤2.5；$x_k-v_i$的内积为$d_{i k}^2=\left\|x_k-v_i\right\|^2=\left(x_k-v_i, x_k-v_i\right)$。$\left\|\boldsymbol{A}_i\right\|=\rho_i, \rho>0$，对于每个类，ρ_i为常数。

由拉格朗日乘数法可得$\boldsymbol{A}=\left[\rho_i \operatorname{det}\left(\boldsymbol{F}_i\right)\right] \frac{1}{n} \boldsymbol{F}_i^{-1}$，其中F_i为第i类的模糊协方差矩阵：

$\boldsymbol{F}_i=\frac{\sum\limits_{k=1}^n\left(\mu_{i k}\right)^m\left(x_k-v_i\right)\left(x_k-v_i\right)^{\mathrm{T}}}{\sum\limits_{i=1}^N\left(\mu_{i k}\right)^m}$

(2)

在给定聚类数目下，满足i＜c＜n。文中采用X表示数据集，取参数$m \in[1.5, 2.5]$定义加权指数。用ε＞0标记终止误差值，符号A表示距离矩阵取值，采用随机或指定方式对初始矩阵U ⁽⁰⁾进行划分，按照以下步骤进行处理：

第1步，计算类中心：

$v_i^{(l)}=\frac{\sum\limits_{i=1}^n\left(\mu_{i k}^{(l-1)}\right)^m x_k}{\sum\limits_{i=1}^n\left(\mu_{i k}^{(l-1)}\right)^m} \quad(1 \leqslant i \leqslant c)$

(3)

第2步，计算模糊协方差矩阵并加入调整单位阵：

$\boldsymbol{F}_i^{(l)}=\frac{\sum\limits_{k=1}^n\left(\mu_{i k}^{(l-1)}\right)^m\left(x_k-v_i^{(l)}\right)}{\sum\limits_{i=1}^N\left(\mu_{i k}^{(l-1)}\right)^m} \quad(1 \leqslant i \leqslant c)$

(4)

$\boldsymbol{F}_i=(1-\gamma) F_i+\lambda \operatorname{det}\left(F_0\right)^{1 /p} I$

(5)

然后，计算F_i的特征值λ_ij和特征向量φ_ij。记$\lambda_{i \max }=\max _{1 \leqslant j \leqslant p} \lambda_{i j}$，当$\lambda_{i \max } /\lambda_{i j}>\beta$时，令$\lambda_{i j}=\lambda_{i \max } /\beta$，β可取10¹⁵，再次调整F_i的取值：

$\boldsymbol{F}_i\left[\phi_{i 1}, \phi_{i 2}, \cdots, \phi_{i p}\right] \operatorname{diag}\left(\lambda_{i 1}, \lambda_{i 2}, \cdots, \lambda_{i p}\right)\left(\phi_{i 1}, \phi_{i 2}, \cdots, \phi_{i p}\right)^{-1}$

(6)

第3步，计算距离并更新划分矩阵：

$D_{i k L_i}^2\left(x_k, v_i\right)=\left(x_k-v_i^{(l)}\right)^{\mathrm{T}}\left[\rho_i \operatorname{det}\left(F_i\right)^{1 /p} F_i^{-1}\right]\left(x_k-v_i^{(l)}\right)$

(7)

$\mu_{i k}^{(l)}=\frac{1}{\sum\nolimits_{j=1}^c\left(D_{i k A_i}\left(x_k, v_i\right)\right) /\left(D_{j k A_i}\left(x_k, v_i\right)\right)^{2 /(m-1)}} \quad(1 \leqslant i \leqslant c, 1 \leqslant k \leqslant n)$

(8)

直至||U ^(l) −U ^(l−1)||＜ε使目标函数最优，从而得到每个类的中心位置以及每个样本点对于每个类的隶属程度，把所有满足隶属度大于阈值T的震源点归并到同一类中。

利用改进的模糊聚类方法（万永革等，2023）可以根据小震多发区的小震定位结果拟合地下的可能断层面，特别是对于无地表出露的地下隐伏断层，该方法具有独特优势。

2 断层几何形状和滑动性质 2.1 断层几何形状

基于余海琳（2023）提供的古冶地震序列精定位结果（图 1），使用万永革等（2023）改进的模糊聚类方法，计算此次古冶M_S 5.1地震所在断层的几何形状（图 2），得到其发震断层A的中心坐标及其走向、倾角及断层面的4个顶点坐标（表 1），其中断层A的中心坐标为（118.44°E，39.75°N），深度为12.48 km。

2.2 古冶地震序列震源区应力场及发震断层面滑动性质

为研究聚类得到的此次古冶地震发震断层（断层A）的滑动性质，首先需要确定研究区应力场，通过应力场投影到断层上的剪应力方向即可确定在此研究区应力场作用下的断层滑动性质。

搜集此次古冶地震震中附近具有震源机制的地震（王晓山等，2020），绘制在图 3(a)中。使用震源机制分类方法（万永革，2024），将表 2中所列地震在球形三角形中分类表示[图 3(b)]。结果显示，其中有61个走滑型地震、31个正断型地震、12个正走滑型地震、3个逆走滑型地震和1个逆断型地震（表 2）。

基于表 2所示地震相关参数，采用根据震源机制求解应力场的网格搜索法（Wan et al，2016），得到该地区构造应力场[图 3(c)、(d)]。古冶M_S 5.1地震震源区构造应力场的最优压轴走向和倾伏角，中间轴走向和倾伏角，以及张轴走向和倾伏角见表 3。该应力场结果与高原等（1995）、严少鹏等（2024）的研究结果较一致，且与华北地区总体构造应力特征基本一致。

将所求的此次古冶地震震源区构造应力场结果投影到表 1中得到的断层面上，得到发震断层滑动角为-177.2°，该结果与诸多学者或机构所得结果（表 4）基本一致。

3 主震对余震的应力触发 3.1 破裂模型

地震本质上是破裂以及位错在断层面上开始、传播和终止的过程。解析地震的破裂模型是研究地震成因、破裂动力学机制和探究活动断层结构等方向的基础，在防震减灾、应急响应和危险性评估等实际问题中发挥着重要作用（岳汉等，2020）。鉴于此次古冶地震震级相对较小，本研究采用的破裂模型是基于Wells等（1994）提出的经验公式计算所得。

据上文计算结果可知，此次古冶地震滑动角为-177.2°，表明其为走滑型地震。据Wells等（1994）给出的走滑地震震源参数经验关系式，计算得到断层破裂的长度、破裂宽度和平均滑动量分别为3.9 km、4.1 km和1.9 cm。

3.2 地壳速度结构

在计算格林函数时需使用地壳速度结构模型，本研究使用的是由CRUST1.0提供的地壳分层速度结构模型，具体参数见表 5。

3.3 主震对余震的应力触发

由图 2可知，此次地震序列的余震集中分布在6 km、12 km及15 km深度附近。已有研究（盛书中等，2015；靳志同等，2019）表明，不同深度上的库仑应力分布存在显著差异。采用Wang等（2006）给出的PSGRN/PSCMP软件包，计算以上3个不同深度上的库仑应力，研究主震对余震的应力触发作用。

PSGRN/PSCMP软件包利用正交归一法计算地震应力场Green函数，理论模拟周围各点运动状态和受力状态。同时，该软件包考虑了变形和地球重力场的耦合效应，弥补了早期出版物中常见的公式错误等不足。

采用余海琳（2023）给出的120个余震精定位结果，将上文计算所得发震断层的走向、倾角和滑动角（230.2°，86.2°，-177.2°）作为接收断层参数，计算主震在不同深度上对余震的应力触发，得到此次古冶M_S 5.1地震在6 km、12 km和15 km深度上的库仑应力分布（图 4），其中余震深度范围分别为：3—7 km、9—14 km和14—18 km（表 6）。

由图 4、表 6可知，在深度6 km附近，余震触发状况不理想。在选择的20个余震中，仅有一个位于库仑应力变化为正区域。可能是因为该深度上余震所在断层与本文设置的发震断层参数不一致。也就是说，这些余震可能发生在其他断层上。可见，接收断层设置参数不同，对库仑应力分布的影响较大。而且，在7—9 km深度无余震发生，说明该区域地质状况较为复杂。由上述区域震源机制分类结果可知，该区域走滑型地震和正断型地震占比均较大，揭示了该区域断层类型分布的多样性和地震活动的复杂性。

据上文所述，此次古冶地震序列所在断层几何中心深度为12.48 km（古冶地震震源深度为10 km），因此选择计算并分析12 km深度上的库仑应力分布，并统计此次古冶地震对[9，14) km范围内余震的应力触发。经统计，在选择的46个余震中，有29个处于库仑应力变化为正区域，17个处于库仑应力变化为负区域，主震对余震的触发比例达63%。同时，计算15 km深度上的库仑应力分布，统计主震对[14，18) km深度范围内余震的应力触发。在选择的54个余震中，有11个处于库仑应力为正区域，43个处于库仑应力为负区域，主震对余震的触发比例约为20%。

综合评估可知，古冶M_S 5.1地震序列余震数目较少，且对120个余震的触发比例仅约34%，触发效果不理想，但对于12 km深度上余震的触发效果较好。经分析，可能原因如下：唐山—古冶断裂属于唐山断裂带的子断裂，古冶地震发生在该子断裂带上。该地区的地质状况较为复杂，李枭等（2024）在古冶—滦县地区上地壳各向异性研究中发现：唐山断裂带附近台站快剪切波偏振方向与断裂走向一致，该区域各向异性程度较强。此次地震对余震的触发状况不好也可能受其影响。

3.4 主震对M_S 4.3地震的应力触发

滦县M_S 4.3地震在M_S 5.1地震后发生，且距其震中较近，分析此次地震发生是否受到古冶M_S 5.1地震的影响具有一定意义。因此，本节主要讨论M_S 5.1地震对M_S 4.3地震的触发作用。

将滦县M_S 4.3地震震源深度（9 km）以及余海琳（2023）给出的该地震的走向、倾角和滑动角（121.23°、85.18°和-0.9°）作为接收断层参数，使用PSGRN/PSCMP软件包，计算古冶M_S 5.1地震在9 km深度上的库仑应力分布，结果见图 5。

由图 5可见，滦县M_S 4.3地震发生在古冶M_S 5.1地震触发库仑应力变化为正区域。经计算，该处库仑应力为0.000 2 MPa，虽未达到0.01 MPa的触发阈值，但库仑应力变化为正，说明古冶M_S 5.1地震对滦县M_S 4.3地震的发生起到促进作用。而且，0.01 MPa只是前人根据经验给出的临界触发值。据Ziv等（2000）对地震累积静态应力变化的分析，当库伦应力变化小于0.01 MPa时，累积应力变化有可能增加，进而影响区域地震活动性。因此表明，古冶M_S 5.1地震对滦县M_S 4.3地震的发生具有触发作用。

4 讨论和结论

在应力触发研究中，库仑应力计算存在一定不确定性，比如：地震的破裂参数、分层地壳介质模型、摩擦系数、接收断层参数等影响因素，需要对库仑应力计算的可靠性和不确定度进行量化。另外，本研究未考虑1976年唐山大地震对此次M_S 5.1地震库仑破裂应力的影响。李枭等（2024）通过对唐山断裂北段区域台站快剪切波偏振优势方向的分析，指出：1976年唐山M_S 7.8地震序列引发的尖端效应在唐山断裂带北段应力集中区的应力场并未发生改变；王想等（2021）结合1976年唐山地震的时间序列，从其余震空间分布以及发震构造角度分析，认为：此次古冶M_S 5.1地震是1976年唐山M_S 7.8地震后的又一次显著地震活动。关于此方面的研究有待进一步讨论。

本文主要就唐山古冶M_S 5.1地震发震构造和M_S 5.1地震对余震的触发作用进行分析，得到以下结论：

（1）基于古冶地震震源区余震精定位，使用改进的模糊聚类方法，得到发震断层的几何参数：走向和倾角分别约为230.2°、86.2°。断层面的4个顶点坐标分别为（118.479°E，39.771°N，4.41 km）、（118.473°E，39.777°N，16.59 km）、（118.403°E，39.732°N，16.59 km）和（118.409°E，39.726°N，4.41 km）。基于震源区构造应力场，得到发震断层滑动角约为-177.2°。

（2）基于经验公式计算得到发震断层破裂模型：断层长约3.9 m，宽约4.1 m，平均滑动量约0.02 m；基于破裂模型和地壳速度模型，将计算所得走向、倾角和滑动角作为接收断层参数，分析认为：M_S 5.1主震对M_S 4.3余震有触发作用，且对分布在9—14 km深度上的余震触发效果较好。

文中图形采用GMT软件（Wessel et al，1998）和MATLAB软件绘制而成，特此致谢。本研究使用的PAGRN/PSCMP软件由汪荣江教授提供，在此表示感谢。对审稿专家提出的建设性意见表示感谢。