采用模糊聚类算法确定2021年玛多地震序列的断层结构

Geometry of Seismogenic Faults Determination of the 2021 Maduo Earthquake Sequence by Fuzzy Clustering Algorithm

摘要

为更深入理解青海玛多区域的地质构造，准确确定该区域的各个分支断层面形状以及断层面参数是关键.本研究按照成丛小震发生在断层面邻近区域的原则，基于模糊聚类算法聚类数据并确定断层面的方法对2021年青海玛多Ms7.4级地震序列进行断层面形状求解，得到了各个分级下的断层结构.然后将研究区域的应力场投影到获得的断层面上，得到了各个断层面的相对剪应力和正应力.结果表明本次地震主干断层面破裂长度为163.71 km，总体走向为285.81°，倾角为85.62°，断层面主要分布于昆仑山口‒江错断裂上.该地震序列的东西两侧表现有明显的分叉现象，针对这一较为复杂的断层结构，通过区域划分和模糊聚类数目的各种尝试，得到了划分为6个断层面的模糊聚类求解断层面结果.将研究区域的构造应力场结果投影到6个断层面上，发现这些断层面上相对剪应力普遍大于相对正应力.由此推断本次地震主要是由构造应力场作用引起，而位于西端分叉处北支的小断裂是主破裂引发的分支破裂，由于本次破裂具有左旋走滑运动性质，使断层东段东北部出现了尾端拉张，导致东段分叉上形成了两个小断层.本研究中通过模糊聚类方法得出的断层参数与其他机构作者数据一致性较高.该方法基于大量地震序列数据可较为准确地得出断层面信息且可操作性高，在未来地震构造、断层面确定以及分析中具有重要意义.

关键词

模糊聚类 / 断裂带形状 / 断层面解 / 玛多Ms7.4级地震 / 发震断层 / 地震学.

Key words

fuzzy clustering / geometry of seismogenic faults / fault plane solution / Maduo 7.4 earthquake / seismogenic fault / seismology

0 引言

据中国地震台网报道（https://news.ceic.ac.cn/CC20210522020411.html），位于青海省果洛藏族自治州的玛多县地区在北京时间2021年5月22日2时4分发生7.4级地震，震中位置为34.59°N、98.34°E，震源深度17 km.美国地质调查局（USGS）给出的地震震级为7.3级，初始破裂点位于34.613°N、98.246°E，震源深度10 km，震源机制解走向N92°E，倾角62°，滑动角为-40°，表现为一条近EW走向，略带拉张分量的左旋走滑错动事件（https://earthquake.usgs.gov/earthquakes/eventpage/us7000e54r/executive）.地震的发生与断层的存在息息相关.本次玛多地震位于青藏高原中北部的巴颜喀拉块体上，巴颜喀拉块体内部发育了包括昆仑山口‒江错断裂、玛多‒甘德断裂、达日断裂、东昆仑断裂等一系列大型活动断裂，而本次玛多地震则处于东昆仑断裂带和巴颜喀拉‒松潘弧形构造带的交会部位，是继2008年汶川Ms8.0地震之后在中国境内发生的震级最高、地表破裂最长的地震事件（李智敏等，2021）.该块体内部发育着较多大型走滑断裂，比如达日断裂、甘德南缘断裂和五道梁‒长沙贡玛断裂等，此次走滑事件更加说明巴颜喀拉块体内部具有大地震孕育和发生的构造条件（梁明剑等，2022）.

研究地震事件的基础工作就是确定该事件所在区域的地震发震断层面的参数信息，比如断层面的位置、走向角、倾角、滑动角等参数.地震中断层面的形状往往不是一个简单的几何平面，所以对于不同类型的断层所采用的研究方法也有很大差异.其中对于一些深部破裂信息，采用深部地震测深需要耗费大量成本，且地球物理反演存在多解性，使得不确定性提高，而普通的地质调查只能看到地表的破裂，因此万永革等（2008）选择通过余震分布的方式来确定参数.

一般认为浅源地震事件中主震发生后，大量余震集中分布在断层面及附近，因此可利用震源点位置的空间分布约束断层面的结构（王福昌等，2013）.王鸣和王培德（1992）根据这种原理分析研究了1989年大同‒阳高地震部分余震的空间分布和综合断层面解，他们将高斯‒牛顿法应用于利用余震分布求断层面最小二乘解的过程中，但是这一方法在迭代过程中对初始参数选择要求十分严格，而且没有给出位置参数以及走向角和倾角的误差估计等参数.万永革等（2008）结合模拟退火算法、反演理论和应力场等相关地震学知识，在对唐山地震的应用中进一步给出了断层面的走向、倾角和滑动角参数的估计值和误差.张广伟等（2011）利用小震空间分布估算出了唐山地震断层面解，从而验证了万永革等（2008）所提方法的准确性及可操作性.王福昌等（2010）为避免万永革等（2008）方法中因使用高斯‒牛顿法造成迭代过程的繁琐，选择主成分分析方法求解断层面走向角和倾角，从而加快了计算速度.近年来，也有许多学者利用余震序列约束断层面，Wang et al.（2020）通过HypoDD定位和中等地震事件震源机制反演的方法分析了2017年西藏米林县地震余震序列并得到地震断层情况.Evans et al. （2020）将新的统计方法用于Community Fault Models（CFM）中，将南加州地震与其可能的发震断层概率相关联，使其更好地辅助识别由主要断层活动引起的小地震群.但上述方法往往只能对一个主断层面进行拟合，对于拥有多个子断层面的复杂破裂结构不能同时处理，并且部分方法操作复杂、工作量大.

随着双差定位法和波形互相关技术的不断发展及其在地震学领域的应用实践，使利用小地震的丛集性来确定断层结构成为可能.王福昌等（2013）在主成分分析算法基础上结合模糊聚类算法、稳健扩充算法和主成分分析法给出了一种可以对活断层网络三维空间结构进行重构的新方法.该方法每次从震源点较密集的区域开始，利用假设检验不断扩充子断层面，并可以获得多个子断层面的参数.本文采用王福昌等（2013）的模糊聚类算法，对青海玛多地震及其余震进行聚类分析，得出符合实际的一系列断层面区域，从而推测本次青海玛多地震事件中断层的几何形态，并与采用模拟退火算法和高斯‒牛顿算法相结合的反演方式所得到的断层面参数以及其他机构作者得出有关玛多地震参数的信息进行对比，结果显示出较高的一致性，验证了模糊聚类算法用于确定断层面的准确性.

1 模糊聚类确定断层面算法

本文采用的方法是基于GK模糊聚类算法（王福昌等，2013）并采用了自适应距离对模糊算法进行修改.在该方法中通过划分区域，设置适当断层面个数并去除离群事件，选择聚类所使用的地震事件处理.利用这种方法，能够直接清晰呈现出该区域多个断层面的三维分布情况以及断层面参数信息.该方法曾被张丽娟等（2022）用于漾濞地震事件研究中，本文参考了其对该方法的描述，GK模糊聚类的目标函数为：

\mathrm{m} \mathrm{i} \mathrm{n} {\mathit{G}}_m {\left(X;M,V,\left.A\right)=\right.{\sum }_{i=\mathrm{1}}^c{\sum }_{k=\mathrm{1}}^n\left({\mu }_{ik}\right.)}^m D_{ik{A}_i}^{\mathrm{2}}，(1)

其中 \mathit{X} = \{ {\mathit{x}}_{\mathrm{1}} , {\mathit{x}}_{\mathrm{2}} , \cdots , {\mathit{x}}_n \}为由n个数据点 {\mathit{x}}_{\mathit{k}} = ( {\mathit{x}}_{\mathit{k}\mathrm{1}} , {\mathit{x}}_{\mathit{k}\mathrm{2}} , \cdots , {\mathit{x}}_{\mathit{k}\mathit{p}} ) \in {\mathit{R}}^{\mathit{p}} （ k = \mathrm{1,2} , \cdots , n）构成的数据集， i 为类号 （i=1，2，…，c）； \mathit{M} = {\left({\mu }_{ik}\right)}_{c\times n}为隶属度矩阵； \mathit{V} = \left\{v_{\mathrm{1}},v_{\mathrm{2}},\cdots ,v_c\right\}， {\mathit{v}}_i \in {\mathit{R}}^p为第i（i=1，2，…，c）个类的中心， R 表示所设置选定的区域；每一子类对应的距离矩阵为 A_i，用于描述数据每一类的局部结构；由拉格朗日乘数法可得 A_i = {\left[{\rho }_i\mathrm{d}\mathrm{e}\mathrm{t}\left(F_i\right)\right]}^{\frac{\mathrm{1}}{n}} S_i^{-\mathrm{1}}，其中 S_i为第i类的模糊协方差矩阵：

{\mathit{S}}_i = \frac{{\sum }_{k=\mathrm{1}}^n{\left({\mu }_{ik}\right)}^m\left(x_k-v_i\right){\left(x_k-v_i\right)}^T}{{\sum }_{i=\mathrm{1}}^N{\left({\mu }_{ik}\right)}^m}，（2）

其中，数据样本序号用 k 表示 （k=1，2，…，n）， {\mu }_{ik}表示第k个元素隶属于第i个类的隶属度， m ( \mathrm{1} \le m \le \mathrm{\infty } )为加权指数； x_k - v_i的内积为 d_{ik}^{\mathrm{2}} = {‖x_k-v_i‖}^{\mathrm{2}} = \left(x_k-v_i\right. , \left.x_k-v_i\right)，因为 ‖A_i‖ = {\rho }_i， \rho > \mathrm{0}，对每个类 {\rho }_i是常数.

当对于数据集 X，首先给定聚类数目，取 ，聚类数目决定了最后结果中的断层数量，加权指数 m，一般取1.5至2.5，终止误差 \epsilon > \mathrm{0}，距离矩阵 A，随机或指定初始划分矩阵 {\mathit{U}}^{\left(\mathrm{0}\right)}.对l=1，2， \cdots，n，则：

Step 1计算聚类中心：v_i^{\left(l\right)} = \frac{{\sum }_{i=\mathrm{1}}^n{\left({\mu }_{ik}^{\left(l-\mathrm{1}\right)}\right)}^m{x}_k}{{\sum }_{i=\mathrm{1}}^n{\left({\mu }_{ik}^{\left(l-\mathrm{1}\right)}\right)}^m},\mathrm{1} \le i \le c .(3)

Step 2计算模糊协方差矩阵：{\mathit{S}}_i^{\left(l\right)} = \frac{{\sum }_{k=\mathrm{1}}^n{\left({\mu }_{ik}^{\left(l-\mathrm{1}\right)}\right)}^m\left(x_k-v_i^{\left(l\right)}\right){\left(x_k-v_i^{\left(l\right)}\right)}^T}{{\sum }_{k=\mathrm{1}}^n{\left({\mu }_{ik}^{\left(l-\mathrm{1}\right)}\right)}^m}， \mathrm{1} \le i \le c，(4)

加入调整单位阵：

{\mathit{S}}_{\mathit{i}} = \left(\mathrm{1}-\gamma \right) {\mathit{S}}_{\mathit{i}} + \gamma \mathrm{d} \mathrm{e} \mathrm{t} {\left({\mathit{S}}_{\mathrm{0}}\right)}^{\raisebox{1ex}{$\mathrm{1}$}\!\left/ \!\raisebox{-1ex}{$p$}\right.} \mathit{I} , (5)

计算 {\mathit{S}}_{\mathit{i}}的特征向量 {\Phi }_{ij}和特征值 {\lambda }_{ij}，记，当 {\lambda }_{i\mathrm{m}\mathrm{a}\mathrm{x}} / {\lambda }_{ij} > \beta时，令 {\lambda }_{ij} = {\lambda }_{i\mathrm{m}\mathrm{a}\mathrm{x}} / \beta， \beta可取为 {\mathrm{10}}^{\mathrm{15}}进行重新调整如下：

{\mathit{S}}_i = \left[{\varphi }_{i\mathrm{1}},{\varphi }_{i\mathrm{2}},\cdots ,{\mathrm{\Phi }}_{ip}\right] \mathrm{d} \mathrm{i} \mathrm{a} \mathrm{g} \left({\lambda }_{i\mathrm{1}},{\lambda }_{i\mathrm{2}},\cdots ,{\mathrm{\lambda }}_{ip}\right) {\left[{\varphi }_{i\mathrm{1}},{\varphi }_{i\mathrm{2}},\cdots ,{\mathrm{\Phi }}_{ip}\right]}^{-\mathrm{1}} . (6)

Step 3计算距离并更新划分矩阵：

D_{ik{A}_i}^{\mathrm{2}} \left(x_k,v_i\right) = {\left(x_k-v_i^{\left(l\right)}\right)}^T \left[{\rho }_i\mathrm{d}\mathrm{e}\mathrm{t}{\left(S_i\right)}^{\raisebox{1ex}{$\mathrm{1}$}\!\left/ \!\raisebox{-1ex}{$p$}\right.}{S}_i^{-\mathrm{1}}\right] \left(x_k-v_i^{\left(l\right)}\right)，\mathrm{1} \le i \le c，\mathrm{1} \le k \le n,(7)

{\mu }_{ik}^{\left(l\right)} = \frac{\mathrm{1}}{{\displaystyle \sum _{j=\mathrm{1}}^c}{\left[D_{ik{A}_i}\left(x_k,v_i\right)/D_{jk{A}_i}\left(x_k,v_j\right)\right]}^{\mathrm{2}/\left(m-\mathrm{1}\right)}}，\mathrm{1} \le i \le c，\mathrm{1} \le k \le n .(8)

使目标函数最优，即 ，从而得到每个类的聚类中心以及每个样本点对于每个类的隶属度.再判断震源点数据是否对某一个类的隶属度大于阈值T，如果大于阈值，就将该震源点数据添加到该类中进行处理，利用隶属度和数据集构建有效性函数，用有效性函数确定初始聚类数，再结合主成分分析最终得出断层面长度和宽度、断层面走向角和倾角等参数.在实际操作中，根据所应用的地震实例，调整聚类数目和类中心得出更确切的结果.

2 2021年青海玛多地震序列断层面的确定

2021年青海玛多地震序列的1 208个余震事件分布如图1所示，数据来自于王未来等（2021）采用双差定位方法对震后8天的余震序列重新进行的精定位研究.从中不难看出本次地震序列整体呈现NWW走向的狭长条带状分布，本次地震主要沿东昆仑断裂带南侧的NWW向江错断裂分布，因为江错断裂向西延伸与发生2001年东昆仑Ms8.1级地震的昆仑山口断裂相连，因此该断裂被统称为昆仑山口‒江错断裂，从区域上看，2021年的青海玛多Ms7.4级地震应该是在2001年东昆仑Ms8.1级地震的破裂基础上，应变能沿着东昆仑断裂带向东南方向进一步迁移的结果.这两次地震与1937年花石峡地震的同震地表破裂存在断裂间隙，江错断裂带的地震危险性应该得到重视（Ha et al.， 2022）.此外，与2008年汶川地震一样，同属于难以预测的大陆内部地震事件（熊熊，2022）.

从地震的分布中可以看出东西两端明显存在分叉现象，基于这一特点可以初步将其分为西、中、东三段.潘家伟等（2021）开展的同震破裂野外调查结果表明此次玛多地震的地表破裂表现为左旋走滑运动性质，而在青藏高原已发生的走滑型地震中尚未报道过这种两端发育较大规模分叉的“扫帚”状同震地表破裂（潘家伟等，2021），且本次事件并不是单一断层面的破裂，这为可以同时求解多个断层面的模糊聚类方法提供了绝佳的应用机会.

为求解本次玛多地震序列的整体破裂特性，笔者先将整个研究区域中数据处理时取聚类个数为1，即把所有地震序列事件看做一条整体断裂处理，进而可以看到其总体形态，结果如图2所示，得到的整体断层面长度为163.71 km，走向为285.81°，倾角为85.62°，发震断层近垂直，这与地表破裂带调查（李智敏等，2021，潘家伟等，2021）、大地电磁探测（詹艳等，2021）、InSAR和GPS同震形变（华俊等，2021； Xiao et al.， 2022）、地震序列精定位研究（王未来等，2021）等相关研究结果大体一致.从中也可明显看出，小地震分布区域与得出的断层面范围较吻合.

为了进一步分析断层面情况，考虑到玛多地震两端具有分叉特征，按照地震分布的疏密情况在之前西、中、东三段的基础上，将西端分为西1段、西2段；将中段分为中西段和中东段；将东段分为东1段、东2段，即总共分为6段，这样的分段也与潘家伟等（2021）的分段方式类似，方便之后断层参数按段进行对比，分段情况如图3所示.

在采用模糊聚类算法的过程中，将聚类个数取4，即采用模糊聚类算法识别出4条断裂面，如图4所示.红色圆圈表示被选择用于聚类的地震数据，并进行主成分分析以确定断层面，按照字母顺序，断层面走向分别为105.34°、105.72°、113.27°、91.93°，断层面倾角分别为81.24°、84.69°、85.23°、89.57°，断裂面长度分别为31.77 km、31.88 km、50.39 km、24.22 km.因为该事件分布特点呈较集中的狭长条带状，所以离群事件较少.

如图4所示，4个断层面基本与图3中西1段、中西段、中东段、东2段吻合，但是当继续增加断层数量后，尝试将断层面个数设置为目标分段6时，西1段和东1段始终不能出现预期中的小断层面，而是在中西段、中东段余震密集的地方出现平行于B、C断层的重复断层.笔者推测可能是西1段及东1段的余震数量相对于中部以及其他部分太少，使得这些区域的地震对于断层面的拟合影响较小从而不能成功识别出局部的小断层面，反而这些分叉处的小地震参与到了中部更为明显的断层面拟合中.为了验证这一猜想并获得分叉处小断层的数据，笔者尝试通过局部选取范围的处理对东西两端开叉的区域进行断层面拟合.

通过对玛多事件进行西中东三段分区域操作，可以克服之前两端分叉处不能出现小断层的情况，进而得到按照图3的6段式断层面分布及参数.断层面分区域选取以及所得断层面如图5所示，利用这种方法所得出的6段断层面的参数如表1所示.

余震的震源性质可以用于约束震源区域复杂的断层结构，本文所使用的聚类模糊方法侧重于使用余震位置进行约束，为更好约束断层面分布，还应考虑用较大余震的震源机制进行约束.张志朋等（2022）使用时间域矩张量反演方法TDMT对本次玛多主震以及9个4.5级以上的余震求解了震源机制解.张建勇等（2022）采用CAP波形反演与P波初动极性反演方法对玛多震源区域的132个余震求取了震源机制.结果显示主震以及大部分余震事件为走滑型，少部分为逆断型，主要分布在东段区域，且大部分余震震源机制解与主震相似，结合余震震源机制与该地区地质构造背景，表明昆仑山口‒江错断裂为本次事件的发震断层，验证了所求结果的准确性.

3 所求断层面结果与其他方法结果对比

为了验证模糊聚类法的准确性，笔者仍然使用相同的地震数据并仍按照图3中的地震分段，采用模拟退火算法并结合高斯‒牛顿算法反演本次青海玛多地震的断层面参数信息，并与上文采用模糊聚类算法所得到的结果进行对比分析.

按照成丛小震发生在断层面邻近区域的原则，利用一定数量小地震的震源位置信息可以反演获得断层面的相关参数.其基本思路是找到一个平面，使所有小地震的震源位置到这个平面距离的平方和尽可能最小（王鸣和王培德，1992）.还是按照图3的6个小段的分法，基于每个小段中所选矩形区域内的余震事件进行断层面的参数反演，结合模拟退火算法和高斯‒牛顿算法的程序（万永革等，2000，2008；王福昌等，2008；Wan，2010），利用反演求得每个小段地震断层面走向角、倾角、距坐标原点的距离以及断层的4个顶点的经纬度、深度等断层面位置信息，所得结果如表2所示.因为各段中余震与断层面关系大致相同，本文仅展示中东段的情况（图6）用于说明性质.因为小地震常常发生在断层面周围，所以小地震与断层面的距离分布具有对称性，也因为这一性质，从图6d中也可以看出，小地震与反演断层面距离的分布频度也大致呈现对称性.可以发现，利用前人方法所得到的反演断层面参数（表2）与利用模糊聚类算法得出的断层面数据（表1）高度一致.本文主要突显模糊聚类算法的结果情况，表2数据仅用于与模糊聚类结果对比，下文仍将采用基于模糊聚类方法的表1中数据进行分析对比.

本次地震之后，潘家伟等（2021）、李智敏等（2021）曾前往现场进行地表破裂地质调查，调查得出本次地震同震地表破裂带总体沿N105°E走向分布，地表破裂带约151 km，这与本文模糊聚类结果总体一致，除了整体情况，本文通过模糊聚类所得到的6段式小断层的走向、倾角、长度等参数也与调查结果基本一致，进一步表明模糊聚类方法的准确性.同时结果的一致性也说明了本次研究中涉及的昆仑山口‒江错断裂的深部断层结构与地表破裂在几何上表现一致，但这种一致性只能维持在本次研究中所用事件的最大深度，在更深处的区域是否能保持一致以及断层深部与表面的其他性质是否一致还需要进一步讨论研究.

除此之外，中国地震局地球物理研究所（https://www.cea-igp.ac.cn/kydt/278249.html）等国内外机构网站以及一些学者给出了该地震的震源机制结果，笔者使用求解中心震源机制解的程序（万永革，2019）分别以各个震源机制为初始解得出了中心震源机制解（表3）.比较得到的最小解作为最终结果，发现将USGS得到的震源机制作为初始解得到的震源机制标准差最小，最终中心解得出的节面Ⅰ：走向101.93°、倾角83.78°、滑动角-3.48°，节面Ⅱ：走向192.3°、倾角86.54°、滑动角-171.62°作为结果，P轴的走向57.23°、倾伏角6.85°，不确定范围分别为40.67°~73.67°和-4.82°~18.47°；T轴的走向327°、倾伏角1.95°，不确定范围分别为310.44°~343.44°和-9.66°~13.54°；B轴的走向221.22°、倾伏角82.88°，不确定范围分别为127.13°~353.03°和79.24°~98.79°.所得出的中心震源机制解及空间三维辐射花样见图7，从表3和图7中观察数据可以看出，本次地震震源机制解距中心解的空间旋转角最大达38.95°，最小为0.52°，这些信息也体现了前人得到的震源机制解比较集中，并与万永革等在Seismology小组公众号中发布的结果高度一致（https://mp.weixin.qq.com/s/60hFBoaBXa-c_Q9BVCNYHw）.笔者利用把所有地震序列事件看做一条整体断裂处理时得到的结果与表3数据对比，结果比较吻合，尤其与张志朋等（2022）、GCMT（https://www.globalcmt.org/CMTsearch.html）结果更为接近.特别注意的是，在这些学者和机构对本次地震给出的地表破裂调查研究以及反演分析的结果中，也存在着各种差异，尤其是USGS给出的震源机制解与其他数据差异较大，且大部分反演结果中走向普遍仅略大于100°而小于地表调查研究结果，造成这种反演误差的原因可能是青藏高原本身地下复杂的地质结构以及反演所使用的速度模型偏差（张喆和许力生，2021）.而本文结果与野外实地调查结果更吻合，更加说明本文结果的准确性.

在本次事件中，西段破裂带有明显的南北两支分叉，南支与中西段相似，更接近主断层，所发生余震数量也更多；相反北支所涉及的余震数少，偏离主断层.中西段与中东段为本次玛多地震余震发生的主要区域，也是主震的发生区域，发生于这段区域的地震事件高达722个，占比约59.77%，这两段的长度分别是26.79 km和51.73 km，占整个破裂带全长的绝大部分，也可以看到从西到东，走向逐渐变大，倾角逐渐变缓，基本与主断层一致.在中东段最东边分布着一些零星、不集中的小余震，这些零星地震在对断层面拟合处理中会影响断层的选取，使得拟合出的结果中断层面不再延伸至最东边，因此将零星的地方不予考虑（图4以及图5e、5f）.东段同样分叉为南北两支，其中南支地震零星稀少，但南支与主断层接近，而北支明显偏离主断层位置，走向呈90°，且东段地震在北支分布集中，所以在整体拟合中能被较准确识别（图4）.潘家伟等（2021）将东段地表破裂带分为北、中、南三支，但实际上大同小异.除此之外，从观察倾角也可看出，6个断层都大致接近垂直.

4 所求断层面的滑动特性

为研究该发震区域的6个断层成因是由应力场构造引起还是被主震破裂所触发的.本文利用模糊聚类的走向和倾角数据（表1），结合Wan （2010）的应力场数据，选择了最接近的经纬度数据，其中压轴走向和倾伏角分别为57°和7°，张轴走向和倾伏角为320°和45°，应力形因子（万永革，2015）R值取0.65，使用滑动角计算程序（万永革等，2008）得出各断层的滑动角数据，再利用相关程序（万永革，2020）对已有数据计算出各断层面的相对剪应力、相对正应力和剪滑角，结果如表4和图8所示.通过对所得结果中相对剪应力和相对正应力的对比，可以看出分布在昆仑山口‒江错断裂上的断层，相对剪应力都普遍远大于相对正应力，估计是由应力场构造导致的剪切运动引起的，而分叉中西2段的这个小断层可能是由本次主震破裂引起.东段的两个小断层的产生原因可能是由于本次破裂是左旋走滑运动性质，进而使东段出现了尾端拉张的情况，导致了东边的分叉情况.从应力情况的分析中，更加能印证青海玛多Ms7.4级地震应该是在2001年Ms8.1级地震的破裂基础上沿着东昆仑断裂带向东南方向应变能进一步迁移的猜测.

5 结论与讨论

本研究利用模糊聚类方法，结合主成分分析对2021年玛多地震进行反演，获得了区域内断层的三维空间结构及具体参数.把足够高密度的区域划分为簇，利用成丛小震发生在断层面附近的原则，直接得到了本次地震序列区域的整体断裂情况和4个主要断层.整体破裂长度为163.71 km，总体走向为285.81°，倾角为85.62°，发震断层近垂直.4个断层面走向分别为105.34°、105.72°、113.27°、91.93°，断层面倾角分别为81.24°、84.69°、85.23°、89.57°，断裂面长度分别为31.77 km、31.88 km、50.39 km、24.22 km.但为了获得更精确的断层情况，将其局部进行处理，进一步得出6个断层面情况下的数据及三维空间结构（图5及表1）.并采用结合模拟退火全局搜索和高斯‒牛顿局部搜索方法的断层面参数反演拟合，获得了6段式情况下的断层面数据（表2）.本次模糊聚类数据与表2数据以及其他学者机构数据进行分析对比，结果较为一致.随后将本次模糊聚类方法得出的数据进行了一些分析应用，得到以下结论与讨论：

（1）本次玛多地震为两端发育有较大规模分叉的“扫帚”状同震地表破裂，笔者推测它们可能是此次玛多地震发震断层两端末梢效应新产生的两条尾端分支断裂，它们与主干断裂整体组成了“树”型，这也基本确定此次地震的发震构造应为昆仑山口‒江错断裂.

（2）采用这种模糊聚类新方法所得出的断层面参数数据与其他学者以及国内外权威机构（表3）给出的数据结果比较，结果基本一致，充分说明本文所使用的这种方法具有准确性和可操作性.

（3）本研究中的模糊聚类方法可以用于比较复杂的断裂带，不用分区域分部处理，可以整体识别出多条断层面，是一种基于大量地震序列数据得出断层面数据的聚类分析方法，并可获得一种统计学意义的空间断裂带形状，使用简单、结果清楚直接.

（4）结合应力场数据，获得了该地震断层的滑动角数据以及相对剪应力和相对正应力数据，据此可大致推测主断裂上的4个断层主要是由构造引力场引起的断层面，而对于两端分叉中的小断层，猜测西2段可能是由主震破裂所触发的小断裂，而东段两个小断层的产生原因可能是由于本次破裂是左旋走滑运动性质，使东段出现了尾端拉张的情况，进而导致了东边的分叉情况.

本文所用的数据来自王未来等（2021）采用双差定位法对震后8天的余震序列重新开展的精定位研究.双差定位的结果在深度上并没有太高分辨率，但可以较好消除地壳速度结构横向不均匀性导致的误差，已经是处理群体事件情况下相对最优的方法，所以这种误差对聚类结果影响不大.另外在事件选择和进行双差定位法时都会丢失一些地震事件，但丢失的事件相对质量较差或不适合重定位，且重定位后震源相对位置的定位精度更高，这些丢失的事件可能对结果的准确性更有益，但对本文模糊聚类的结果影响不大.

由于所得到的断层面是使用隶属度大于阈值的地震数据得出的结果，所以余震的分布会影响聚类的选择，并进一步影响聚类分析的结果，而该结果可以同时得出多个断层面，所以难以分辨是哪些地震在某个具体断层面的聚类结果中起了贡献.并且由于它是将范围内的地震全部都引入考虑，最终按照提前设置的断层面个数给出断层面结构，所以在某些局部由于地震数量太少，受主体区域中大量地震的影响会导致无法得出位于这些局部的小断层，而是在地震密集的区域沿着主断层的走向重复得出断层面.目前本文所采用的方法是分步分区域处理.

参考文献

[1] Deng, Q. D., Zhang, P. Z., Ran, Y. K., et al., 2002. Basic Characteristics of Active Tectonics of China. Science in China (Series D), 32(12): 1020-1030, 1057 (in Chinese).

[2] Evans, W. S., Plesch, A., Shaw, J. H., et al., 2020. A Statistical Method for Associating Earthquakes with Their Source Faults in Southern California. The Bulletin of the Seismological Society of America, 110(1): 213-225. https://doi.org/10.1785/0120190115

[3] Ha, G. H., Liu, J. R., Ren, Z. K., et al., 2022. The Interpretation of Seismogenic Fault of the Maduo Mw7.3 Earthquake, Qinghai Based on Remote Sensing Images—A Branch of the East Kunlun Fault System. Journal of Earth Science, 33(4): 857-868. https:// doi.org/10.1007/s12583-021-1556-2

[4] Hua, J., Zhao, D. Z., Shan, X. J., et al., 2021. Coseismic Deformation Field, Slip Distribution and Coulomb Stress Disturbance of the 2021 Mw7.3 Maduo Earthquake Using Sentinel-1 InSAR Observations. Seismology and Geology, 43(3): 677-691 (in Chinese with English abstract).

[5] Li, Z. M., Li, W. Q., Li, T., et al., 2021. Seismogenic Fault and Coseismic Surface Deformation of the Maduo Ms7.4 Earthquake in Qinghai, China: A Quick Report. Seismology and Geology, 43(3): 722-737 (in Chinese with English abstract).

[6] Liang, M. J., Huang, F. P., Sun, K., et al., 2022. The Holocene Activity and Its Evidence from Paleoearthquake of the Middle Segment of Wudaoliang-Changshagongma Fault Inside the Bayan Har Block. Earth Science, 47(3): 766-778 (in Chinese with English abstract).

[7] Pan, J. W., Bai, M. K., Li, C., et al., 2021. Coseismic Surface Rupture and Seismogenic Structure of the 2021-05-22 Maduo (Qinghai) Ms7.4 Earthquake. Acta Geologica Sinica, 95(6): 1655-1670 (in Chinese with English abstract).

[8] Wan, Y. G., 2010. Contemporary Tectonic Stress Field in China. Earthquake Science, 23(4): 377-386. https://doi.org/10.1007/s11589-010-0735-5

[9] Wan, Y. G., 2015. A Grid Search Method for Determination of Tectonic Stress Tensor Using Qualitative and Quantitative Data of Activate Faults and Its Application to the Urumqi Area. Chinese Journal of Geophysics., 58(9): 3144-3156 (in Chinese with English abstract).

[10] Wan, Y. G., 2019. Determination of Center of Several Focal Mechanisms of the Same Earthquake. Chinese Journal of Geophysics, 62(12): 4718-4728 (in Chinese with English abstract).

[11] Wan, Y. G., 2020. Simulation on Relationship between Stress Regimes and Focal Mechanisms of Earthquakes. Chinese Journal of Geophysics, 63(6): 2281-2296 (in Chinese with English abstract).

[12] Wan, Y. G., Shen, Z. K., Diao, G. L., et al., 2008. An Algorithm of Fault Parameter Determination Using Distribution of Small Earthquakes and Parameters of Regional Stress Field and Its Application to Tangshan Earthquake Sequence. Chinese Journal of Geophysics, 51(3): 793-804 (in Chinese with English abstract).

[13] Wan, Y. G., Wu, Z. L., Zhou, G. W., et al., 2000. How to Get Rake Angle of the Earthquake Fault from Known Strike and Dip of the Two Nodal Planes. Seismological and Geomagnetic Observation and Research, 21(5): 26-30 (in Chinese with English abstract).

[14] Wang, F. C., Cao, H. R., Wan, Y. G., 2010. Application of Linear Errors-in-Variables Model for Determination Main Earthquake’s Fault Parameters of Wenchuan Earthquake. Journal of Applied Statistics and Management, 29(3): 381-390 (in Chinese with English abstract).

[15] Wang, F. C., Wan, Y. G., Hu, S. T., 2008, Application of Particle Swarm Optimization to the Estimation of Mainshock Fault Plane Parameters. Journal of Seismological Research, 31(2): 149-154, 198 (in Chinese with English abstract).

[16] Wang, F. C., Wan, Y. G., Qian, X. S., et al., 2013. A New Method for Determining Fault Planes Parameters According to Earthquake Clustering. Chinese Journal of Geophysics, 56(2): 522-530 (in Chinese with English abstract).

[17] Wang, M., Wang, P. D., 1992. Focal Mechanism and Seismogenic Structure of Datong Yanggao Earthquake on October 18, 1989. Acta Seismologica Sinica, 14(4): 407-415 (in Chinese with English abstract).

[18] Wang, W. L., Fang, L. H., Wu, J. P., et al., 2021. Study on Precise Location of Maduo Ms7.4 Earthquake Sequence in Qinghai in 2021. Science in China (Series D), 51(7): 1193-1202 (in Chinese with English abstract).

[19] Wang, W. P., Yang, J. S., Wang, Y. B., 2020. Seismic Sequences of the 2017 Mainling Mw 6.5 Earthquake: Imaging the Seismogenic Fault by Aftershock Analysis. Pure and Applied Geophysics, 177(7): 3161-3174. https://doi.org/10.1007/s00024-020-02422-2

[20] Xiao, L. H., Zheng, R., Zou, R., 2022. Coseismic Slip Distribution of the 2021 Mw7.4 Maduo, Qinghai Earthquake Estimated from InSAR and GPS Measurements. Journal of Earth Science, 33(4): 885-891. https://doi.org/10.1007/s12583-022-1637-x

[21] Xiong, X., 2022. What is the Mechanism to Control the Spatial Distribution of Earthquakes within the Continent? Earth Science, 47(10): 3906-3907 (in Chinese).

[22] Zhan, Y., Liang, M. J., Sun, X. Y., et al., 2021. Deep Structure and Seismogenic Pattern of the 2021. 5. 22 Madoi (Qinghai) Ms7.4 Earthquake. Chinese Journal of Geophysics, 64(7): 2232-2252 (in Chinese with English abstract).

[23] Zhang, G. W., Lei, J. S., Xie, F. R., et al., 2011. Precise Relocation of Small Earthquakes Occurred in North China and Its Tectonic Implication. Acta Seismologica Sinica, 33(6): 699-714, 843 (in Chinese with English abstract).

[24] Zhang, J. Y., Wang, X., Chen, L., et al., 2022. Seismotectonics and Fault Geometries of the Qinghai Madoi Ms7.4 Earthquake Sequence: Insight from Aftershock Relocations and Focal Mechanism Solutions. Chinese Journal of Geophysics, 65(2): 552-562 (in Chinese with English abstract).

[25] Zhang, L. J., Wan, Y. G, Wang, F. C., et al., 2022. Geometry of Seismogenic Faults of the 2021 Yangbi Earthquake Sequence Determined by Fuzzy Clustering Algorithm. Seismology and Geology, 44(6): 1634-1647 (in Chinese with English abstract).

[26] Zhang, Z., Xu, L. S., 2021. The Centroid Moment Tensor Solution of the 2021 Mw7.5 Maduo, Qinghai, Earthquake. Acta Seismologica Sinica, 43(3): 387-391 (in Chinese with English abstract).

[27] Zhang, Z. P., Li, J., Feng, B., et al., 2022. Precise Location and Focal Mechanism Solutions of the 2021 Maduo, Qinghai Ms7.4 Earthquake Sequence. China Earthquake Engineering Journal, 44(1): 218-226 (in Chinese with English abstract).

基金资助

国家自然科学基金项目(42174074;41674055)

河北省地震科技星火计划项目(DZ20200827053)

中央高校科研业务费专项(ZY20215117)

河北省地震动力学重点实验室开放基金项目(FZ212105)