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

张丽娟 万永革 王福昌 靳志同 崔华伟

1)防灾科技学院，三河 065201 2)河北省地震动力学重点实验室，三河 065201 3)山东省地震局，济南 250102

0 引言

活动断层的几何形状在地球动力学以及地质构造研究中具有非常重要的作用(王鸣等，1992； 万永革等，2008)。针对发震断层面的研究发现断层面的空间几何分布特征往往非常复杂，且介质的性质在大范围内具有非均匀性和地质构造多变性的特点，一般而言并非一个确定的简单平面。对于出露在地表之上的断层，可以测量出露断层面的几何形态以对地震深部的几何形态进行推断(Xuetal.，2002，2009)。然而，来自地表的外营力，如地面的剥蚀、 雨水的冲刷作用、 地面生物生长等会给断层地质调查带来很大的不确定性。对于隐伏断层，可通过开挖探槽等地质调查的方法厘定其空间位置(冉勇康等，1997，2018； 张培震等，2003)。然而，地质调查仅能得到浅部的断层形态，而断层深部和浅部的形态有时差别较大，如张先康等(2002)通过深部地震反射发现1679年的三河-平谷大地震的深部与浅部构造存在巨大差别。人工地震测深是一种确定深部活动断层几何形态的方法。然而，该手段需要使用大量人力物力对深部地震剖面进行探测，且地球物理反演的多解性也会给资料解释带来诸多不确定性。考虑到地质调查只能看到地表的破裂，对于深部破裂信息可通过余震分布来确定(万永革等，2008)。

通常地震常发生在预先存在或新激活的断层面上(Waldhauseretal.，2002； Waldhauser，2009； Rossetal.，2020)，因而最直接通常也最为有效的方法是通过对震源位置进行精准定位从而刻画断层面的精细结构。王鸣等(1992)对1989年大同地震中部分余震目录的空间分布特征和断层面的空间结构进行了研究和讨论。 其采用最小二乘法基于余震分布求解断层面，并用高斯-牛顿法进行了迭代，然而这种方法对初始参数的选择比较敏感，且研究结果也并未给出断层面的空间位置等参数。王福昌等(2010)为避免高斯-牛顿迭代算法过程中的缺陷并提高运算速度，采用主成分分析的数值分析方法对断层面的倾角和走向进行了分析和讨论，但依然存在未解决的问题，如地震断层面的各参数估计值可能会受离群震源点的严重影响，从而导致估计结果失真。有学者通过建立模型并考虑定位过程中的误差，基于应力场及反演理论等地震学有关知识，采用模拟退火算法对1976年唐山地震进行了讨论(万永革等，2008)，给出了断层面的滑动角、 倾角及走向等参数，并进行了误差分析。以上方法一般只能基于存在一个主断层面的前提下进行拟合，而不能对含有多个子断层面的复杂地震破裂结构进行处理。有学者对华北地区2001—2009年的地震进行精确定位(张广伟等，2011)，给出了小地震分布的一些分布规律，并采用小地震空间分布估算得出断层面的各个参数。该方法未明确给出划分小地震群的数学模型和具体实现方法，从侧面对万永革等(2008)所给出方法的实用性及有效性进行了验证。Ouillon等(2008)虽然给出了一种比较繁琐的划分小地震的方法，并将其应用于Landers地震的研究中，但其计算速度仍然较慢。王福昌等(2013)在主成分分析的基础上利用几何模糊聚类算法给出了同时推测多个子断层面的统计聚类算法，极大提高了计算速度和地震断层面精度。

目前关于地震断层面走向等空间精细的定量描述大多给出的都是主要断层的走向，不能同时估计多个断层面的空间分布，对于地震断层面空间位置的刻画依然有很大的改进空间。随着数字地震台站在中国大陆的广泛布设和双差定位方法(Waldhauseretal.，2002； Waldhauser，2009)及主事件定位方法(周仕勇等，1999)的提出，测定的震源数据精确度越来越高，能够测定的地震震级下限也越来越低。在地震数量急剧增加的实际背景下，地震震源目录资料中蕴含的大量关于地震所处断层的几何特征信息，给精准的空间定位的数学算法提供了很好的使用机会。万永革(2022)基于密度的空间聚类方法对2021年云南漾濞地震序列的震源机制数据进行了处理，得到2簇节面中心与地震序列的空间密度线性分布较为一致的分支断层，提供了独立于地震波、 大地测量、 地质资料等之外的确定断裂带形状分布的方法，但并没有对聚类结果来自于哪个地震进行统计分析，也未给出断层面的空间定位。本文采用王福昌等(2013)提出的模糊聚类算法，对云南漾濞地震的震源数据进行聚类分析，并对断层面区域进行了最优化处理，在符合假设条件的前提下，给出了一系列最优解，并通过一系列研究对地震序列的几何形态进行精准推测和评估。

1 资料与算法

1.1 资料与数据来源

中国地震台网正式测定，2021年5月21日21时48分云南大理白族自治州漾濞彝族自治县发生了M6.4地震，震中位置位于(25.67°N，99.87°E)(图 1)。震后，云南省地震局立即成立指挥部并启动三级响应，分配工作队奔赴现场，开展震情监视与灾害损失评估等工作。灾区附近的电力、 交通、 水利等基础设施均受到了一定程度的影响。经初步统计，地震有感范围约 1.2 万平方千米，大理市地区及漾濞周边各县普遍有震感。截至2021年5月27日18时，漾濞地震序列已致全县共15786人受灾，共造成经济损失约2.021亿元。

图 1 云南漾濞地震序列的精定位数据分布图Fig. 1 The relocation results of the Yunnan Yangbi earthquake sequence.云南漾濞的精地震序列定位用圆圈表示； 5级地震用红色五角星表示； 主震用黄圆圈表示。图中用不同颜色的圆圈表示地震深度。 黑色实线F1为维西-乔后断裂。虚线多边形为研究区所在的主要区域，椭圆形区域内为第2类数据区域

地震发生后，许多国内外研究者对此次地震的震源机制进行了研究(龙锋等，2021； 余海琳等，2021)，对余震进行了精准定位(Yangetal.，2021)，并开展了GNSS与InSAR破裂滑动分布特点及同震变形场特征(于书媛等，2021； 张克亮等，2021)、 震源区电性结构特征(叶涛等，2021)、 震源区地壳速度结构(杜广宝等，2021； 李大虎等，2021)等研究，为后续的防震减灾及地震研究等领域提供了基础数据。位于维西-乔后及巍山附近的一条隐伏的断裂带(雷兴林等，2021)是本次地震的发震断层。此次云南漾濞地震主要在四川与云南块体偏西的边缘区域孕震，为该区有记载以来发生的最大地震。本次地震以走滑运动为主，主震之后余震活动频繁，其地震序列发生在地震监测台网较为密集的区域，通过密集的地震观测台站可以较好地确定震级下限较低的地震目录。梁姗姗等(2021)采用双差定位方法对此次地震开展了研究，结果显示地震的前震和余震序列存在显著的时空差异性。其中，前震呈条带状，主要分布在NW-SE向，且前震震中的空间位置还存在迁移和往返的特点。而余震呈现出不对称和共轭形态，主要分布在NW和NE向(图 1)。地震带南端NW向还有多组余震沿不同走向分布，揭示了震源区的介质特性和断层几何结构方面的复杂性。这种复杂特征不能仅模拟为单独一条断层，这为本文采用地震精定位结果同时求解多条断层的方法应用提供了绝好的机会。

本文基于中国地震台网中心给出的(99°～101°E，24.5°～26.5°N)空间范围内的观测报告结果，使用Waldhauser等(2002)提出的双差定位法，对2021年5月18日—7月5日云南漾濞共7092个地震事件的101806条Pg和Sg到时数据进行了处理，对这些地震事件进行重定位。定位时，设置Pg的权重为1.0，Sg的权重为0.8。用10km作为相邻地震事件的最大搜索阈值，相邻地震事件的数量最多为20个，严格限制地震事件，并选择至少6个震相记录实现定位，最终得到了5827个地震事件的精确定位结果(图 1)。其中，90.23%的事件被重新定位，其结果见表1。

表 1 双差重定位使用的一维速度模型Table1 Double-difference relocation results by 1-D velocity model

定位结果显示，漾濞地震序列整体呈NW-SE向分布，长约25km，位于维西-乔后断裂西南侧3～10km处，整体走向约为135°；M6.4主震震中位于(25.688°N，99.877°E)，震源深度约9.6km，M>5.0地震3次，分布在余震区的不同部位，余震区中段的最大前震为M5.6，北西端的主震为M6.4，余震区南东段的最大余震为M5.2。综合地震序列深度剖面和震源机制解结果可知，发震断裂应为NW走向，并且可见整个序列地震的密集区域主要限定在主震的SE侧，仅有少部分的地震分布在NW侧，整体上空间展布呈现北西较窄、 南东较宽的特点。因此，可认为单侧破裂特征较为明显，主震方向为SE(雷兴林等，2021)。此外，地震序列在空间结构整体散落的疏密度上也有显著差异性： 南东段相对较为分散，而北西段则更加密集，且在离开序列主干部分的北西段存在2个小地震丛集，似乎有与整体序列的NW-SE走向相偏离的趋势(龙锋等，2021)。这些复杂的地震序列分布为同时求解多断层的参数提供了应用契机。

1.2 模糊聚类算法

本文采用GK模糊聚类算法(王福昌等，2013)，使用自适应距离对模糊算法进行了修改，每个子类对应距离矩阵Ai用于描述数据的局部结构。在FCM算法中将目标函数修改为GK-FCM算法的目标函数：

(1)

(2)

在给定的聚类数目下，满足10标记终止误差值，符号A表示距离矩阵取值，采用随机或指定方式对初始矩阵进行划分U(0)，按照以下步骤进行处理：

第1步，计算类中心：

(3)

第2步，计算模糊协方差矩阵：

(4)

并加入调整单位阵：

Fi=(1-γ)Fi+γdet(F0)1/pI

(5)

Fi=[φi1，φi2，…，φip]diag(λi1，λi2，…，λip)(φi1，φi2，…，φip)-1

(6)

第3步，计算距离：

(7)

第4步，更新划分矩阵：

(8)

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

利用不同的隶属度取值和数据集的集散程度构造有效函数，用有效函数最优的方法确定初始聚类数，再结合主成分分析给出95%置信水平的子断层面长度和宽度、 断层面走向角和倾角等参数的估计值。

2 确定2021年云南漾濞地震序列的断层面

从2021年云南漾濞地震序列的双差定位结果在地表的投影图(图 1)可以看出，虚线多边形区域和实线椭圆形区域地震序列明显属于不同的地震断层，将2部分数据分开进行处理。处理实线区域内的数据时，误差值取ε=10-6。 处理虚线多边形区域内的数据时，取聚类个数为3，将实线椭圆形区域单独处理，聚为一类。

图 2 确定的断层面和地震位置的平面图(a)和立体图(b)Fig. 2 Determined fault planes and hypocenters in map view(a)and in 3-D(b).方框为确定的断层面； 圆圈为聚类所用的属于断层面的数据； 蓝色叉号为聚类过程中判定的离群事件

图 2 中用叉号表示离群的震源点，其余的各子类震源点以及未被指定给任何子类的震源点均采用圆圈标记。图中还给出了漾濞地震序列各震源点的空间分布及断层面在地面上的投影情况。用主成分方法对这些属于某条断层的事件数据进行分析并确定矩形特点的断层面，将各子断层投影至地表上，均为四边形，并用A、 B、 C对其进行标记。沿NW-SE走向断裂的3个断裂面的走向分别为134.22°、 132.65°、 129.45°，断裂面的倾角分别为87.14°、 81.96°、 74.77°，断层面的范围与余震分布能够较好吻合。

本研究采用模糊聚类算法识别出3条断裂面，与Yang等(2021)推断的“漾濞地震破裂的东南端有至少3条断裂”的结果相吻合。不同之处是Yang等(2021)识别的第3条断裂是由西南侧的一簇地震导致的，但那一簇地震并未形成条状地震带，究竟其是不是一条断层，还需要通过其他地球物理探测手段证实，而本文识别出的3条断裂带中的2条位于主要断裂的方向。万永革(2022)采用震源机制节面聚类方法对2021年云南漾濞地震序列中的震源机制节面进行聚类，识别出2条断裂面。断裂面的分布与维西-乔后断裂带的方向几乎平行，且向SE按照地震分布的密集程度逐渐分为2支断裂。西侧分支对应本文的断裂面A、 B，应为主断裂； 东侧分支对应本文的断层面C，与本文的研究结果基本一致。

第2类数据(图 1 中的实线椭圆形区域)位于维西-乔后断裂及巍山断裂附近，与主断裂面斜交，小地震的丛集深度为7～10km。该区域共有276条数据，分支看起来较为明显，但其面状结构不明显，且震源点较少，故对其进行单独处理。采用模糊聚类算法识别为“断层面”。该“断层面”的走向为235.66°，倾角为66.30°，倾向NW(图 2)。从断层面位置的分布来看，这簇地震带位于维西-乔后断裂带附近，需要指出的是仍需通过其他地球物理探测手段确定该分布究竟是否为一条地震断层面，本文称其为伪断层面。由位于维西-巍山-乔后断裂附近的小地震集中活动，推测该地区后期存在持续地震活动的可能。

表2和图3总结了前人得到的断层面形状结果，包括基于InSAR、 GPS及地震波数据等不同资料利用不同方法得到的结果。前人所给出的断层面倾向既有SWW的结论，也有NEE的结果。从本文得到的3个主要断裂面的结果来看，可得到以下结论： 3个断裂面总体几乎垂直于水平面，与已有结果较为吻合，这验证了本研究结论的正确性。

表 2 不同研究组或机构给出的漾濞地震断层面的几何参数Table2 The seismic fault planes with its geometry parameters determined by the previous researches

图 3 其他作者得到的断层几何形状与本研究结果的比较Fig. 3 Comparison of the research results of other researches with that of this paper.其他作者给出的断层形状采用点弧线表示，其中心形状采用粗弧线表示，实弧线为本研究通过聚类分辨的3个断层面的形状。黑点表示其他作者给出的断层面极点位置

为更好地分析研究区所在的主要区域的断层面走向，取数据区域(图 1 中的虚线多边形所围成的区域)的聚类个数为4，得到如下结果(图 4)： A、 B、 C、 D 4个断层面的走向分别为146.01°、 135.55°、 112.09°、 134.994°，倾角分别为84.03°、 88.03°、 82.77°、 83.82°。从结果来看，B、 D断层面为主断层面，与本文前部分主断层面的方向较为接近，从一定意义上验证了漾濞地震的主断层面方向具有一致性，且从另一侧面验证了算法对主断层识别的可靠性。在东侧数据密集的区域识别出的2个断层面的走向与前文虚线多边形所围研究区所在的主要区域的断层面走向略有差异。考虑到大地震破裂一般具有优势断层断裂面，这会导致余震发生的微小裂纹至少有一部分与优势断层的断裂方向具有一致性。本研究中的模糊聚类算法的主要目的之一是精确刻画地震目录所呈现的优势断层面的空间分布，并探讨其他断层面的空间分布。即使某一部分地震数据识别出“断层面”也不一定表明其为真正的地震断裂面。这是由于主震破裂导致邻域内的应力场发生较大改变，致使相同几何形状裂纹的错动方向有较大的变化范围，从而导致地震断裂面的形状出现较大的变化范围。

图 4 按4个断层面处理的主破裂地震数据确定的断层面和地震位置的平面图(a)和立体图(b)Fig. 4 Determined fault planes and hypocenters in map view(a)and in 3-D(b) in main rupture area and 4 fault planes.图例同图3

3 结论与讨论

本研究结合主成分分析，利用GK模糊聚类方法给出了2021年云南漾濞地震活动断层的网络三维空间结构。以簇来标记密度足够高的区域，依据成丛的小地震在断层面附近分布且在附近密集的原则，得到了与该地震序列分布走向大体一致的3个主要分支断层面和1个隐伏断裂面。3个主断层面的走向分别为134.22°、 132.65°、 129.45°，倾角分别为87.14°、 81.96°、 74.77°，所得结果与地震震源数据的一致性较高。由研究结果推测2021年云南漾濞地震主要断层面的分布与维西-乔后断裂近平行，并逐渐分为3个分支，西侧分支的倾向大体为SWW，对应2个主断裂面，并向SE逐渐分为3个断裂分支。在维西-乔后断裂带附近还有一条断层，走向为235.66°，倾角为66.30°，称之为伪断层。

本研究中使用的模糊聚类法适用于分支断裂面比较复杂的断裂带，可同时识别多条断层面，是基于大量地震震源断裂面数据进行处理的一种聚类分析算法，从理论上可能得到统计意义的空间断裂带形状。使用该方法分析地震断层面时，如能配合地震精定位的地震形态展布、 地质等其他手段进行综合分析，将显著增强研究结果的应用价值。本文使用模糊聚类法在漾濞地震序列中得到了4个地震断层面，位于维西-乔后断裂带附近的断层面与主要断裂面的差异较大，考虑可能是地震发生时与地下断层和地质结构存在某种关系。同样，必须指出的是，本文所述方法在云南漾濞地震序列研究中得出的主断层面的空间位置、 2个断层面的倾向分布、 2条断裂的分支断裂之间的连接性及其之间如何连接等相关科学问题仍须借助如地球物理探测技术等其他手段作进一步研究和考证。

另外，本文所得的结果难以评估地震对结果的贡献度，也不能分辨哪个地震将对结果造成主要影响，而是根据每个震源点的数据隶属度按照模糊聚类原则整体确定的断裂面。算法程序在普通计算机上的运行速度已经较快，且可以实现在短时间内进行多次重复计算，如果再考虑随机给出初始划分矩阵，得到的分类结果可能会略有不同，可根据其他探测技术、 地质条件以及地质结构进行分析讨论，以得到更加准确的结果，这一点在王福昌等(2013)关于Landers地震的处理过程中也曾有描述。这种方法可以很好地依据余震的空间分布特点确定大地震子断层空间的定位及其参数，进而可成功解读断层的几何学特征、 地震序列活动的活动范围，并可在无需人工指定的前提下自动完成对地震断层面的中心位置以及整个断层面空间方位的刻画。若根据小地震分布预先通过人工给定一些断层面的中心并构造初始划分矩阵，收敛速度将会更快，并得到更好、 更加精准的定位结果。

文中所提模糊聚类求解断裂带形状的方法，快速计算分析得到与地震分布大体一致的分支断裂带的形状及其空间分布，进一步验证了该方法在地震断层面空间定位的研究领域具有一定应用价值和参考意义，可将这种基于数据分析的算法应用于对地震断裂带的形状分析和预判。在地震发生后进行灾害评估及开展地震灾前预测时，对断层结构的网络性研究结果将对地震附近区域的应力传递检验及相关假设提供帮助(Wanetal.，2010)。

致谢本文的绘图使用GMT及MATLAB软件绘制完成； 审稿专家为本文提出了建设性修改意见，对稿件的质量提升帮助很大。在此一并表示感谢！