断裂带震源机制节面聚类确定断裂带产状方法及在2021年漾濞地震序列中的应用

万永革

1 防灾科技学院，河北三河 065201 2 河北省地震动力学重点实验室，河北三河 065201

0 引言

活动断层几何形状是研究地球动力学和地震危险性的基础资料.对于出露地表的断层，可以通过测量出露断层面的几何形态推测地震深部的几何形态(如Xu et al., 2002, 2009)，然而雨水冲刷、地表剥蚀、地面生物作用等外营力的作用导致断层的地质调查具有一定的不确定性.对于隐伏断层，则可以通过挖探槽等地质调查来测量(如张培震等，2003；冉勇康等，1997, 2018)，但浅部破裂形态和深部断层可能有很大差异，如张先康等(2002)通过深地震反射揭示的1679年三河平谷大地震与浅部构造存在巨大差别.人工地震测深确定活动断层几何形态(如刘保金等，2012；李燕等，2017)，需要大量的人力物力，且地球物理反演多解性使得资料解释存在很多不确定性.

大地震发生后一段时间内，大量余震在断层面上及其附近发生，因此余震震源位置的空间分布可以较为精确地勾画出断层面的形状和位置.万永革等(2008)假定地震发震断层可以用一个平面来模拟，且设大多数余震发生在这个断层面的附近，则可以通过余震震源位置参数来求解发震断层的走向、倾角及位置，并根据研究地区的应力场参数估计断层的滑动角.虽然如此，这种方法依赖于丛集的余震，如果余震未出现丛集也使得方法失效，另外这种方法对地震定位的精度要求较高，而地震深度的精确确定一直是地震学的一项富有挑战性的工作.

当周围布设有密集的大地测量测站的大地震发生后，通常可以通过大地测量资料来求解断层的几何形状和滑动性质(陈运泰等，1979；Shen et al., 2009; Wan et al.，2017).然而大地震本来就很稀少，而且能满足这样密集大地测量资料条件的活动断层更少，因此，活动断层形状的精确测定一直是地球科学家孜孜以求的工作目标.

目前震源机制(其中一个节面为断层面)的测定方法是一种省时省力的方法，这方面发展比较迅速，有P波初动求解断层面解的方法(俞春泉等，2009)、P波初动和P/S振幅比结合的方法(Kisslinger et al., 1981; Snoke et al.,1984; 梁尚鸿等，1984; 吴大铭等，1989; Hardebeck and Shearer, 2003)、近震体波波形方法(倪江川等，1991；Dreger and Helmberger, 1993; 姚振兴等, 1994; Herrmann, 2013； 杨宜海等，2017)、CAP方法(Zhao and Helmberger, 1994; Zhu et al., 1996；易桂喜等，2012)、面波方法(Udias, 1971; Aki and Patton, 1978; Patton, 1980; Kanamori and Given, 1981)、远震长周期地震波的矩心矩张量方法(Dziewonski et al., 1981)、W-Phase确定方法(Duputel et al., 2012)等等.随着数字地震台站在中国大陆的广泛布设，使得测定的震源机制解的震级下限越来越低，目前测定的震源机制数量急剧增加.这些地震震源机制资料也蕴含了大量关于地震所处断层几何特征的信息，但震源机制节面的大量数据也增加了选择的困难，并且震源机制节面之一为辅助面(非地震破裂面)更增加了问题的复杂性.目前还未见自丛集地震震源机制节面数据定量提取活动断层几何参数的相关文献.

本文拟对丛集地震震源机制中提取断层几何信息的方法进行研究，对发生在断层上的中小地震震源机制节面进行聚类分析，得到断层面走向和倾角的估计方法，并给出其置信区间.从而提供一种区别于地质方法、拟合断层面和大地测量确定断层面的另外一种方法.

2021年的漾濞地震序列发生在地震监测台网较为密集的区域，并且该地震是一个典型的前震-主震-余震性序列，通过密集的地震观测台站可以较好地确定地震的震源机制，为上述方法的应用提供了绝好的机会.本文对这个地震序列的大量震源机制的节面进行聚类分析，从而推测该地震序列所发生的断层的几何形态.

1 方法

1.1 聚类分析

由于研究问题是：采用基于密度分布的聚类方法，找出密集分布的震源机制解节面的中心作为可能的活动断层的几何形状，本小节首先介绍基于密度的聚类方法.密度聚类方法的指导思想是，只要一个区域中对象的密度大于某个阈值，就把它加到与之相近的聚类中去.对于簇中每个对象，在给定的半径ε的邻域中至少要包含最小数目k个对象.这类算法能克服基于距离的算法只能发现“类圆形”的聚类的缺点，可发现任意形状的聚类，且对噪声数据不敏感.这类算法中比较有代表性的算法为DBSCAN(Density-Based Spatial Clustering of Application with Noise )，它将簇定义为密度相连点的最大集合，能够把具有足够高密度的区域划分为簇，并可在有“噪声”的空间数据库中发现任意形状的聚类(Ester et al., 1996; Sander et al., 1998; Daszykowski et al., 2001).基于密度的聚类方法需要设定两个参数：其一为密度，采用一个聚类对象的ε-邻域至少包含最小数目k个对象来表征，其二为聚类对象的邻域半径ε.根据前人研究(Daszykowski et al., 2001)，最小数目k经验设定为

k=int(m/25),

(1)

其中，m为需要分类的数据总个数.聚类对象的邻域半径ε表达为

(2)

有了上面两个参数，即可按照DBSCAN算法进行聚类.其具体步骤如下：(1)将需要分类的所有数据对象均标记未被访问，给出数据的总个数m和维数n；按照公式(1)和(2)给出类的最小数目k和邻域半径ε；(2)从所有数据对象中逐个抽取数据点，判断该点是否已被访问，如果是，则跳过该点，否则进入步骤(3).(3)首先标记该数据点为已访问数据点，求解该点周围的邻域半径ε内的数据点数是否超过k，如未超过，则标记该数据点为噪声点(邻域内只有一个数据)或边界点(邻域内有多个数据，但未超过k)，继续检查下一个数据；否则该点为核心点，产生一个新类，对该类进行编号，并将该核心点及其邻域点的序号放到该类的序号集合中，进入步骤(4).(4)对该簇中所有数据点逐个寻找其邻域点，每处理一个数据点，则在该序号集合中去除该数据点的序号，免得重复操作.不妨设一个数据点的邻域点为b.若邻域点b已访问，则不做任何处理.若未被访问，首先标记该数据点b已被访问，然后将该点归为这一类中，并序号放入序号集合中.如此递归操作直至处理完该簇中的所有数据点(即序号集合中为空)就得到了该类的所有数据点.(5)转入步骤(2)的下一个数据，直至全体数据点均被访问，并将未被分类的点标记为噪声点(具体过程见图1).

1.2 震源机制中节面的定量差别及类中心平均值的解法

上一小节是一般的DBSCAN聚类方法.其中统计数据点邻域里的点的个数需要计算两个数据点之间的距离，通常用欧式距离来表达.而对于震源机制节面的聚类，数据点是两个节面.为表达它们之间的差别，我们首先求出节面的单位法向矢量.设震源机制节面的走向和倾角分别为φ和δ，则其单位法向在北东下坐标系下可表示为(万永革，2016)

n=[-sinφsinδ,cosφsinδ,-cosδ],

(3)

若两个节面的法向分别n1和n2，则两个向量之间的夹角为

α=arccos(n1·n2),

(4)

这就是两个节面差别的定量表达.采用该距离作为两个震源机制节面数据之间的距离.由于节面的法向与其相反方向对该节面的描述是一致的，如果两个向量之间的角度大于90°，则取180-α.

另外，在DBSCAN聚类方法中，一个聚类中心为其中数据点的平均值.这里的数据对象为震源机制节面的单位法向矢量.但单位法向矢量的反方向与该单位法向矢量是相同的.为此，仿照万永革(2019)求取同一地震多个震源机制中心解的思路，本研究先求解单位法向矢量各维的平均值组成一个尝试的单位平均法向矢量，然后计算尝试的单位法向矢量与类中所有震源机制节面的法向矢量的夹角，对于夹角大于90°的，取其相反方向为求类中心的震源机制节面的单位法向矢量，再将类中震源机制节面的单位法向矢量各维数据进行平均得到类中心的单位法向矢量.

图1 本研究所用的DBSCAN算法流程图Fig.1 Flow chart of the DBSCAN algorithm used in this study

2 在2021年云南漾濞地震序列中的应用

北京时间2021年5月21日21时48分34秒，云南省漾濞县(北纬25.67°，东经99.87°)发生MS6.4地震(Yang et al., 2021b).在6.4级地震发生之前，已经发生多次前震，震级最大可达5.8级，6.4级地震之后又有多次余震活动，最大余震可达5.2级.余海琳等(2021)基于中国地震台网中心提供的5月18日14时28分至6月8日5时53分的漾濞地震序列观测报告，读取了2543个P波初动数据，选择最少有8个P波初动的地震事件求解震源机制解，得到85个ML≥2.5震源机制解.其中矛盾比为0的震源机制有26个，占总数的30.6%，矛盾比在0～0.10范围内的震源机制占总数的32.9%，矛盾比分布在0.1～0.18范围内的震源机制数占总数的28.2%，矛盾比大于0.18的震源机制仅占总数的8%.由于震源机制有两个节面，地质上定义断层类型通常根据滑动角，因此两个节面上不同的滑动角会导致断层类型判别的困难.本研究根据Zoback(1992)在世界应力图中给出的震源机制分类方法(表1)，所确定的漾濞地震序列的震源机制中有正断型14个，走滑型40个，正走滑1个，逆断型4个，不确定型26个，可见本地震序列总体以走滑型震源机制为主.具体漾濞地震的震源机制参数见余海琳等(2021)的文章，震源机制分布见图2.

采用第二节所述方法对震源机制节面(85个震源机制，170个节面)进行聚类，可以得到三簇聚类中心.第一簇节面数为27，其中心节面法向的走向为227.41°，倾伏角为2.96°，标准差为10.67°，其中心节面的走向为317.41°，置信区间为312.26°～322.55°，倾角为87.04°，置信区间为83.19°～89.55°.第二簇节面数为28，其断层面法向的走向为308.64°，倾伏角为 4.27°，标准差为10.49°，其中心节面的断层面走向为 38.64°，置信区间为 34.02°～ 43.25°，倾角为85.73°，置信区间为80.92°～89.46°.第三簇节面数为12，其断层面法向的走向为 29.42°，倾伏角为 4.42°，标准差为9.79°，其中心节面的断层面走向为断层面走向为119.42°，置信区间为118.79°～120.05°，倾角为85.58°，置信区间为80.16°～89.01°(图3a—c).第一类和第二类的类中心节面的夹角为81.04°，第一类和第三类的类中心节面的夹角为19.43°，第二类和第三类的类中心节面的夹角为80.51°.这表明：得到的第一类和第三类的中心节面与第二类中心节面接近垂直.噪声节面数据个数为103个(图3d)，并且在全空间分布，较为随机.

表1 震源机制解分类依据表Table 1 Classification of focal mechanism solutions

图3 2021年云南漾濞地震序列震源机制聚类结果震源机制节面采用绿色弧线表示，聚类中心节面采用红色弧线表示.黑点表示震源机制节面的极点位置，红点表示聚类中心的极点位置.聚类中心极点周围的蓝绿色椭圆为聚类中心的置信区间. (a)—(c) 分别为聚类的第1—3类的聚类结果； (d) 聚类后的噪声节面数据.Fig.3 The result of clustering of the nodal planes of focal mechanisms in the 2021 Yangbi earthquake sequenceGreen curves show nodal planes of focal mechanisms, and red curve shows the central nodal plane of the cluster. Black dots show the poles of the nodal planes and red one shows the pole of the central nodal plane of the cluster. (a)—(c) Result of the 3 clusters of nodal planes. (d) The noise nodal planes after clustering analysis.

为将本文聚类的节面中心与实际地震序列的位置进行比较，将Yang等(2021a)采用2021年5月18日至5月28日发生的地震事件进行双差重定位的2144个地震事件按照震级大小绘于图2.可以看到，地震分布沿与维西—乔后断裂近乎平行的断裂带分布，并且向东南方向逐渐分为两支断裂，按照地震分布的密集程度推测的两支断裂见图2的黑色虚线位置.西侧分支大体对应于本研究得到的第三类聚类中心，倾向为南西西，应该对应于主断裂；东侧分支大体对应于第一类聚类中心，倾向为北东东，这两支断裂相向分布.而第二类聚类中心大体与本文聚类的第一类和第二类中心近乎垂直，根据地质信息判断第二类聚类中心为本次地震序列破裂的辅助面，不是断层面.

Yang等(2021a)根据其地震双差定位结果推断地震破裂的东南端有至少三条断裂，本研究采用震源机制节面聚类方法识别了两条.Yang等(2021a)识别的第三条断裂是由图2的“第一类聚类结果”的西南侧的一簇地震导致的，从地震震源位置分布来看，那一簇地震还没有形成条状地震带分布.究竟是不是一条断层，还需要其他地球物理探测手段证实.

表2总结了前人采用GPS，InSAR和地震波数据等不同资料和方法得到的漾濞地震断层面的几何形状(如果是两个节面，选择北西-南东走向的节面作为断层面)，本研究将这些结果绘于图4，可以看到前人所给的断层面有倾向南西西的结果，也有北东东的结果.将这些结果按照第二节的方法进行法向平均，得到的断层面法向中心为：走向46.01°，倾伏角3.81°，标准差为9.24°.得到平均断层面走向为136.07°，置信区间为135.07°～137.07°，倾角为86.19°，置信区间为79.76°～89.37°.可见总体来讲断层面接近垂直.

表2 前人得到的漾濞地震的断层面几何参数Table 2 Geometry parameters of the seismic fault planes determined by the previous authors

图4 其他作者得到的断层几何形状与本研究结果的比较其他作者给出的断层形状采用点线弧线表示，其中心形状采用粗弧线表示，实弧线为本研究聚类分辨的两个断层面形状.黑点表示其他作者给出断层面的极点位置，五星表示本研究给出的断层形状极点位置.Fig.4 Comparison of the fault planes provided by other authors and that of this studyDotted curves show fault planes by other authors, thick curve shows the average fault plane of the other authors, solid curves show the fault planes determined in this study. Dots show the poles of the fault planes provided by other authors, star shows the pole of the average fault plane provided by other authors.

3 结论与讨论

本研究根据发生在一个地震带上的大量震源机制节面有可能存在与断裂带几何形状一致节面集中区的假设，把具有足够高密度的区域划分为簇，并可在有“噪声”的空间数据库中发现任意形状的聚类的DBSCAN方法进行改造，给出了独立于地震波、大地测量、地质等资料之外求解断裂带走向和倾角的算法，将该算法应用于2021年云南漾濞地震序列中，得到了与该地震序列的地震分布走向大体一致的两个分支断裂带，验证了该种方法的有效性.

本文基于同一断裂带上的震源机制解其中一个节面与断裂带形状基本一致，对震源机制解的节面数据进行聚类分析估计断裂带形状的.这种假设适合于分支断裂不太复杂的断裂带, 是基于大量地震震源机制解节面的一种去除“噪声”的统计结果.该结果从理论上可能得到一种统计意义的断裂带形状，并且可以得到结果的置信区间.

通常主震断层破裂发生的地震震源机制是丰富多彩的.即使是同一应力场，大地震破裂导致的各种形态微小裂纹上的余震震源机制也会呈现多种类型，同一应力场不同断层表现震源机制的模拟就说明了这一点(万永革，2020)，观测的余震序列，如汶川地震(胡幸平等，2008；易桂喜等,2012)、芦山地震(林向东等，2013；罗艳等，2015)、九寨沟地震(杨宜海等，2017)等，就表现具有丰富类型的震源机制.但由于大地震破裂具有优势断层，至少有一部分导致余震发生的微小裂纹跟优势断层方向趋于一致，本文聚类的目的就是找到这一部分节面所刻画的断层面.但即使这一部分导致余震发生的微小裂纹上的震源机制也不一定是一种类型.这是由于主震破裂导致邻域内应力场有较大改变，致使相同几何形状裂纹的错动方向有较大变化范围，导致了震源机制的辅助面形状有较大的变化范围，这些地震震源机制就表现为不同类型.因此不同类型的震源机制也有可能具有相同几何形状的断层面形状，在进行聚类分析时没有必要划分为不同的震源机制类型进行聚类.

对于本次漾濞地震序列所在区域，黄小龙等(2015)通过地质调查发现本次漾濞地震序列的西北部由三条近乎平行的断裂组成的炼铁盆地东缘主边界断裂组成，本文没有分辨出三条断裂，原因之一是那三条断裂近乎平行，在节面上没有区别；其二是因为本文所给的是一种聚类分析的统计方法，只能对断裂形状进行“大轮廓”的，“写意”的分类.

需要指出的是，一个地震震源机制解有两个节面，本研究将所有震源机制节面放在一起进行聚类分析，得到的聚类结果难以分辨出是哪个地震对聚类结果的贡献.这是因为聚类结果是对所有节面去除“噪声”后的结果，而且还有对非断层面的辅助节面的类的贡献.即同一个地震的震源机制即可能出现在与实际断层面较为一致的簇中，也可能出现在与实际断层面垂直的节面簇中.因此本文没有对聚类结果来自于哪个地震进行统计分析.另外，采用Daszykowski 等(2001)的经验设定聚类簇含有的最小数目k按(1)式计算，此时得到的k为3.如果k取4，聚类得到5簇数据，第一簇含有16个数据，对应的断层面走向为311.67°，置信区间为309.36°～313.98°,倾角为85.12°,置信区间为83.18°～87.05°；第二簇含有24个数据，断层面走向为36.64°，置信区间为33.28°～40.00°,倾角为85.29°,置信区间为80.44°～89.85°；第三簇含有9个数据，对应的断层面走向为119.45°，置信区间为118.74°～120.16°,倾角为89.77°,置信区间为86.79°～87.28°；第四簇含有9个数据，对应的断层面走向为180°，倾角为0°,这是一种不可能的断层面；k=4对应的结果图为图5(由于第四簇与地表平行，未绘出)，与文中结果比较同样聚类得出了倾向西北和东南的两个断层面(第一簇和第三簇).如果取k=5，聚类得到两簇数据，第一簇含有17个数据，对应的断层面走向为311.75°，置信区间为308.35°～315.15°,倾角为85.99°,置信区间为82.10°～87.92°，第二簇含有24个数据，对应断层面走向为36.64°，置信区间为 33.28°～40.00°,倾角为85.29°,置信区间为80.44°～89.85°，对应的结果图为图6.可以看出，采用k=5，仅得到一个东南倾向的断层面，这是由于k的取值使得密度和邻域半径的综合效应使得仅分辨出一簇最为主要的断层，可能反映了两断层的综合结果.由此看出，虽然k的不同取值得到的聚类结果有一定差别，但都得到了类似的结果，但某些结果由于邻域最小数目k和邻域半径的关系可能使得某些结果表现不太理想，需要用户在运用该方法的过程中尝试不同的k值.但究竟如何根据数据的多少、误差等确定k值仍然是以后继续研究的问题.

本文的输入数据是同一地震带上的大量地震震源机制解的节面数据，若地震震源机制解数据不准确，也会导致本方法聚类结果的不稳定，因此保证输入的震源机制的准确性是致关重要的.尽管如此，本文以云南漾濞地震序列的较为传统的P波初动资料(不是地震波形或大地测量的精确数据)求解大量的震源机制进行统计分析，也得到了两个与地震分布大体一致的两条分支断层走向和倾角，这表明本研究的算法有一定稳健性.

图5 除k=4，(a)和(c)分别为聚类的第1—3类，(d)为噪声外，与图3相同Fig.5 Same as Fig.3, except that k=4, (a)—(c) are the 3 clusters of nodal planes, and (d) is the noise nodal planes

图6 除k=5，(a)—(b) 分别为聚类的第1类和第2类，(c)为噪声外，与图3相同Fig.6 Same with Fig.3, except that k=5, (a) and (b) are the 2 clusters of nodal planes, and (c) is the noise nodal planes

图7 采用Fit_Fault软件和精确定位的漾濞地震序列地震精定位数据拟合断层面的结果小震分布在水平面(a)、断层面(b)和垂直于断层面的横断面(c)上的投影，(d)小震距断层面距离的分布，圆圈表示精确小震定位，粗线表示的是断层面边界，AA′为断层上边界端点，DD为倾向，DF为距断层面距离，SD为走向距离，F为分布频次.Fig.7 Result of fitting the fault plane according to precisely located Yangbi earthquake sequence by using Fit_Fault software(a)—(c) are the projections of small earthquakes on horizontal plane, fault plane and cross section perpendicular to the fault plane. (d) The distribution of the distance between small earthquakes and fault plane. The circle represents the precise location of small earthquakes. The thick line represents the boundary of fault plane. AA′ is the endpoint of the boundary above the fault. DD is the dip direction. DF is the distance from the fault plane. SD is the strike direction. F is the frequency of distribution.

我们采用万永革等(2008)开发的采用精确小震位置拟合断层面的程序Fit_Fault，同样利用Yang等(2021a)采用双差定位方法得到的漾濞地震序列的位置进行求解，得到断层面走向为136.0°±0.3°，倾角为88.9°±0.8°，对小震的拟合情况见图7.由于万永革等(2008)所开发的程序只能拟合1个断层面，这个拟合结果是我们采用震源机制聚类分析得到的两个断层面的综合，并且与本研究的震源机制节面聚类识别的断层面相差不大.但采用震源机制节面聚类可以识别出两个断层面，说明本文开发的程序有一定的优势.

本研究聚类分析的结果通常需要配合地质、地球物理探测等手段进行综合分析，如本方法应用于漾濞地震序列，还需要结合地震精定位的地震形态展布.如本研究应用于漾濞地震序列得到三类结果，但第二类结果与地震序列的空间分布不一致而舍去，可以认为是地震破裂的辅助面(非断层面).同样在本文应用于漾濞地震序列得到的两条分支断裂相向的倾向分布以及两个分支断裂是否相连、如何连接等问题也需要其他地球物理探测进行进一步的检验和研究.

尽管如此，本研究提出了采用震源机制节面聚类求解断裂带形状的方法，并将其应用于云南漾濞地震序列，分析得到与地震分布大体一致的两条分支断裂带形状，说明该方法具有一定的预示作用，可以将其应用于大量震源机制的断裂带来约束断裂带的形状.

本研究研发的根据同一断裂带上发生地震震源机制进行聚类分析求解断层面形状的程序可向感兴趣的读者提供.

致谢Yang等(2021a)在网络上公布其漾濞地震序列的双差定位结果供本研究使用，审稿人提供了建设性修改意见，防灾科技学院研究生黄少华同学在绘图方面提供了帮助，本文绘图采用MATLAB软件和GMT软件(Wessel and Smith, 1998)绘制而成，特此致谢.