0 引言

震源机制解是描述震源几何特性的重要参数，是研究断层错动方式和分析地壳应力场的基础资料(Gephart et al，1984；许忠淮，1985)。为捕捉强震前的震兆信息，既可以利用不同阶段不同类型的震源机制解来反映强震前后构造应力场的时空变化(许忠淮等，1989；刁桂苓等，2011)，也可以利用震源机制的一致性规律判别强震前的小震群序列类型(陈颙，1978；程万正等，2006)。但无论以上何种方法，均需要对震源机制解进行有效的类型划分。早期对震源机制解的对比分析多依赖专业知识和经验，如今，对于大量震源机制解的产出，直观观察震源机制解的异同或单纯利用人工处理分析显然不够精确。因此，利用机器学习对震源机制解进行准确快速的类型区分具有重要意义。

聚类分析是一种多元统计分析方法，其针对大量数据或者样品，根据数据本身特性将个体或者对象分类，使得同一类对象之间相似性比其他类对象的相似性更强(方开泰等，1982；何晓群，2004)。震源机制解是一组具有多参数、多维度的数据，聚类分析方法可以有效地对大样本下的震源机制解集进行聚类划分，并提取出潜在的地震学特征。刁桂苓等(1992)提出了一种借助模式识别中系统聚类的方法，把两个解的P轴夹角与T轴夹角之和作为距离衡量，对震源机制解进行聚类划分。通过研究海城地震序列震源机制解，验证了该方法的可靠性和实用性，但由于其相似性距离的度量标准过于简单，迭代终止条件选取的较为主观，对聚类结果产生了一定影响。随后，许多学者针对震源机制的聚类做了大量研究工作，例如采用体波谱振幅相关系数聚类分析震源机制的变化(Lund，2002；朱航等，2006；崔子健等，2012)、利用波形相关性和聚类分析确定小震群的震源机制解(Shelly et al，2016)、通过震源机制节面聚类提取断层几何参数(万永革等，2019)、震源机制一致性的显著性检验法研究(郭祥云等，2019)等。上述工作多是针对大震前后或小区域事件集的震源机制相似性程度进行研究，无法体现大空间和大时间跨度下震源机制的差异性。另外，聚类的方法仍采用层次聚类等传统方法，类与类之间的划分取舍也存在一定的主观性。

谱聚类(Spectral Clustering，简称SC)算法是一种基于图论的聚类方法，是目前聚类分析中的一个研究热点，其本质是将聚类问题转化为图的最优划分问题。其优势在于能在任意形状的样本空间上聚类，且收敛于全局最优解，因此该方法适用于三维球空间的震源机制解(高琰等，2007；蔡晓妍等，2008)。同时，由于谱聚类采用了主成分分析(PCA)中降维的思想，这对于处理多维的震源机制数据也十分有利。基于此，本文尝试采用谱聚类方法对震源机制解进行聚类分析，利用最小旋转角为相似度矩阵，以间隔统计量法的最优解确定聚类数，从而实现对海量震源机制解的自动快速类型划分。为检验这种方法的准确性与实用性，本文利用小区域、小数据量以及大区域、大数据量这两组震源机制解实例，分别进行了谱聚类分析，并针对聚类结果，给出该方法的优势以及适用范围。

1 方法 1.1 谱聚类

谱聚类算法来源于谱图划分理论(Fiedler，1973)，其主要思想是将n个样本数据(x₁，x₂，…，x_n)看作空间中的点，用边将这些点之间连接起来，组成一个基于样本相似度的无向加权图G=(V，E)。其中，V(v₁，v₂，…，v_n)表示图G的顶点集合，E表示具有权重值的边的集合，定义权重w_ij为点v_i和点v_j之间边的权重，由于是无向图，所以w_ij=w_ji。假设距离较远的两点间的边权重较小，而距离较近的两点间的边权重较大，寻找一种切图的方式令切图后子图的边权重和尽可能低，而子图间的权重和尽可能高，从而达到聚类的效果。目前，谱聚类的算法已在很多场景中得到应用，其具体原理也在前人研究中经过反复讨论(Ng et al，2002；高琰等，2007；蔡晓妍等，2008；张宪超，2017)。因此，本文首先将在简要介绍该方法原理基础上，研究针对复杂的震源机制解数据所采取的聚类策略和方法。

谱聚类的基础之一就是确定相似矩阵W的生成方式。相似矩阵主要通过样本点之间的距离度量来构建，在三维空间中一般采用欧式距离等函数。但对于震源机制解，普通的空间距离无法真实客观地描述不同断层的空间差异。Kagan(1991)采用了一种三维旋转方法，可以利用4种旋转模式来实现一个震源机制到另一个震源机制的旋转。这种旋转角度客观地反映了两个震源机制在三维空间的差别。因此，本文采用Kagan(1991)的算法，将两震源机制的最小旋转角作为距离度量。

利用距离构建相似矩阵一般有三种方法：$\epsilon $- 邻近法、K邻近法和全连接法(张宪超，2017)。$\epsilon $- 邻近法和K邻近法均给出了一种近似关系，针对大数据量具有计算速度快的优势，但不可避免会遗漏掉部分细节信息。为真实反映每个震源机制解的差异性，本文采用全连接法中的高斯核函数(RBF)将距离值转换为相似矩阵，即

$ w_{i j}=\exp \left(-\frac{\left\|x_i-x_j\right\|_2^2}{2 \sigma^2}\right) $

(1)

其中，x_i-x_j为两个震源机制的最小旋转角；σ为高斯核比例因子，用来控制样本点的邻域宽度。将w_ij组成矩阵形式，就可以得到相似矩阵W。

根据边的权重值还可以计算出度矩阵D。度d的定义是样本点和其相连的所有边的权重之和，公式表示为

$ d_i=\sum\limits_{j=1}^n w_{i j} $

(2)

利用度可以组成一个n×n的对角线矩阵，即度矩阵D为

$ \boldsymbol{D}=\left(\begin{array}{cccc} d_1 & \cdots & \cdots & \cdots \\ \vdots & d_2 & \vdots & \vdots \\ \vdots & \vdots & \ddots & \vdots \\ \cdots & \cdots & \cdots & d_n \end{array}\right) $

(3)

其中，度矩阵的对角线表示每个样本点的度，其他矩阵位置均为0。

谱图理论实现数据聚类问题的关键在于分割子图，即通过划分准则寻找最优解，这是一个NP难问题(Non-deterministic Polynominal-Hard)，通常解决的方法是采用谱方法将该问题松弛化为连续实数值，使其转换为求解拉普拉斯矩阵的谱分解，然后再将松弛化的问题离散化，得到最终的聚类结果，其本质是对图划分准则的逼近(Zhou et al，2003)。目前，划分准则包括最小割集准则(Minimum cut)、规范割集准则(Normalized cut)和比例割集准则(Ratio cut)等(蔡晓妍等，2008；Wu et al，1993；Shi et al，2000；Hagen et al，2002)。本文采用常见的规范割集准则(简称NCut)计算归一化的随机游走拉普拉斯矩阵，NCut切图的函数表达式为

$ \operatorname{NCut}\left(A_1, A_2, \cdots, A_k\right)=\frac{1}{2} \sum\limits_{i=1}^k \frac{W\left(A_i, \overline{A_i}\right)}{{vol}\left(A_i\right)} $

(4)

这里假设聚类数为k，A_i表示第i类样本数据集，A_i为除A_i子集外其他子集的并集，W(A_i，A_i)表示A_i和A_i之间切图的权重之和，vol(A_i)表示第i类所有样本的度之和。W(A_i，A_i)和vol(A_i) 表示为

$ W(A, B)=\sum\limits_{i \in A, j \in B} w_{i j} $

(5)

$ {vol}(A)=\sum\limits_{i \in A} d_i $

(6)

谱聚类最终的优化目标是最小化NCut函数，这里需要用到维度规约的思想去近似求解NP难问题，即计算标准拉普拉斯矩阵最小的前k个特征，然后求出对应的特征向量，并标准化，再对最后的特征向量矩阵进行传统的k-means聚类，即可得到最终的聚类结果。标准化拉普拉斯矩阵为

$ \boldsymbol{L}_{s y m}=\boldsymbol{D}^{-\frac{1}{2}} \boldsymbol{L} \boldsymbol{D}^{-\frac{1}{2}}=\boldsymbol{I}-\boldsymbol{D}^{-\frac{1}{2}} \boldsymbol{W} \boldsymbol{D}^{-\frac{1}{2}} $

(7)

其中，D表示度矩阵，L为拉普拉斯矩阵，数学表达式为$ \boldsymbol{L}=\boldsymbol{D}-\boldsymbol{W}$，L_sym为标准化后的拉普拉斯矩阵。

1.2 间隔统计量法确定聚类数

在聚类分析中，聚类数目的多少直接影响聚类的质量。从谱聚类的原理中可以发现，聚类数k需要用户给出，这对于简单数据可以很直观或有目的地得到。但是对于震源机制这种多维多属性的复杂数据集，显然无法直观获取聚类的最优数目。因此如何自动地确定最佳聚类数目至关重要。通常解决这一关键问题的方法是利用评价指标穷举聚类数，然后对每一个聚类结果进行评估，根据评估结果确定最优的聚类数目。这一类方法虽然计算量大且耗时严重，但给出的结果是最为可靠的。目前评价方法主要有平方误差和法(SSE)、Calinski-Harabasz准则法、轮廓系数法和间隔统计量法(Gap statistic)等(Calinski et al，1974；Rousseeuw，1987)。本文主要采用间隔统计量法(以下简称GS)评估震源机制解聚类数目。

GS聚类评价方法是由Tibshirani在2000年提出的，其思想主要来源于方差分析(Tibshirani et al，2001)，通过待分数据集的离散度与由参考分布生成的数据集的离散度的差异构造一个间隔统计量，用以刻画间隔统计量关于聚类数目的变化情况，当间隔统计量取得最大值时，对应的k为最佳聚类数(窦婷，2017)。

在应用于震源机制解的聚类数时，假设震源机制解数据集已有谱聚类算法聚成k类：C₁，C₂，…，C_k，震源机制解x_i和x_j的距离平方即最小旋转角的平方表示为d_ij，则C_r类(r=1，2，…，k) 的类内紧凑度为

$ D_r=\sum\nolimits_{x_i, x_j \in C_r} d\left(x_i, x_j\right)=\sum\nolimits_{x_i, x_j \in C_r}\left\|x_i-x_j\right\|^2 $

(8)

所有类内紧凑度总和，即平方误差和W_k可以表示为

$ W_k=\sum\limits_{r=1}^k \frac{1}{2 n_r} D_r $

(9)

其中，n_r表示第r类中震源机制解的数量。从式(9)可以看到，随着聚类数的增大，每一类震源机制解越来越相似，距离平方和越来越小，因此W_k值必定随着震源机制的聚类数增大而减小。但是，当W_k减小的较为缓慢时，可以认为增加聚类数并不能进一步改善聚类效果，这个减小缓慢点即为“肘部”拐点，此时的聚类数k可以视为最佳的震源机制解类别数。如何定量判断“肘部”拐点是GS方法中的一项关键工作。

根据震源机制解的维度，利用蒙特卡洛抽样方法建立参考数据集，并计算每个参考数据集平方误差和，建立Gap值的关系式，即

$ \operatorname{Gap}_n(k)=E_n^*\left\{\lg \left(W_k\right)\right\}-\lg \left(W_k\right) $

(10)

$ E^*\left\{\lg \left(W_k\right)\right\}=\frac{1}{B} \sum\limits_{b=1}^B \lg W_{k b}^* $

(11)

其中，E^*{lg(W_k)}表示lg(W_k)的期望，B为蒙特卡洛抽样次数，W_kb^*表示第b个参考数据集内平方误差和。由于参考数据集随机数的非线性特征，需要利用标准误差sd(k)进行修正，即

$ s d(k)=\left[\frac{1}{B} \sum\limits_{b=1}^B\left(\left(\lg \left(W_{k b}^*\right)-\bar{l}\right)^2\right]^{\frac{1}{2}}\right. $

(12)

$ \bar{l}=\frac{1}{B} \sum\limits_{b=1}^B \lg W_{k b}^* $

(13)

定义$s_k=\sqrt{1+\frac{1}{B}} s d(k) $，使式(14)成立的最小的k_min就是通过GS法得到的最佳聚类数。即当最佳聚类数为k_min时，Gap值应为Gap(k) 中的最大值

$ \operatorname{Gap}(k) \geqslant \operatorname{Gap}(k+1)-s_{k+1} $

(14)

1.3 震源机制平均解的求取

在完成震源机制解的聚类后，需要对每一类震源机制的空间状态进行描述，这就需要寻找一个平均解来体现断层面的空间几何特性。本文主要采用刁桂苓等(1992)求平均解的方法，对同类震源机制的P轴或T轴进行矢量合成，计算出同类解的力轴平均位置，从而得到震源机制的平均解。

1.4 数值实验

为研究将谱聚类方法应用于震源机制解的聚类中是否可行以及结果是否可靠，本文根据生成的一个合成震源机制目录进行方法验证。通过改变3个断层面参数：走向={90°，135°，180°}、倾角={30°，60°}和滑动角={-90°，0°，90°}来生成一个三维的格状网(3×2×3)(Martinez-Garzon et al，2014)。这些参数的不同组合生成了多种断层模式，从纯走滑到斜滑再到纯正断或逆断，每个网格点随机生成符合正态分布的100个震源机制。为直观区分每一类的震源机制，将标准差σ设置为5°。利用上文给出的震源机制解谱聚类方法和GS方法，得到结果如图 1所示。

图 1(b)为利用GS方法得到的Gap值随聚类数变化曲线，图中可以明显看到，在聚类数k为18时，Gap值最大，因此可以认为18是此震源机制数据的最佳聚类数，这也符合之前的预设(网格点3×2×3=18)。图 1(a)是当k=18时，合成震源机制的聚类情况，图中不同颜色表示不同的类别，可以看到每一簇聚类都基本被区分出，聚类效果良好。通过数值实验可以看到，本文研究的震源机制解谱聚类方法是可行的，结果稳定且可靠。

2 应用实例与计算 2.1 数据选取

为检验本文震源机制解谱聚类方法的实用性，本文选取了两个区域的地震震源机制解集作为主要研究对象，用以分析震源机制解的聚类效果，分别为1975年海城M_S7.3地震序列和川滇及周边区域地震事件集。这主要基于以下3个原因：①两组不同区域的地震均有明确的震源机制解集，其计算结果已被不同学者或机构验证，可以确保震源机制解的准确性和可靠性；②1975年海城M_S7.3地震序列已有学者做过类似的聚类分析(刁桂苓等，1992)，通过聚类结果可以对比分析震源机制解谱聚类方法相较于其他方法的优劣；③川滇块体和周边区域研究范围较广，且数据量较大，区内断裂特征明显，有利于通过聚类结果分析区域地震活动与构造特征的关系，从而探讨揭示区域内不同的孕震环境。

1975年海城M_S7.3地震序列采用顾浩鼎等(1976)的结果，共有24个震源机制解，其具体分布如图 2(a)所示。川滇及周边区域震源机制数据主要采用罗钧(2013)的结果，其利用CAP方法反演得到了该区域2007—2013年间共75次3.5级以上地震震源机制解。为了扩大聚类时的数据量，另外补充了1976年至2020年12月31日GCMT记录的该区域141个震源机制解，总计216个震源机制结果，具体分布如图 2(b)所示。

2.2 谱聚类计算

震源机制解是一个包含节面和P轴、T轴和B轴的12维的数据：2个节面的走向(strike)、倾角(dip)和滑动角(rake)，以及P轴、T轴和B轴的方位角(trend)和倾伏角(plunge)。尽管12个参数中只有3个是独立的，即可以由节面的3个参数或P轴2个参数加上T轴的1个参数求出其余的参数(中国地震局监测预报司，2017)，但是为了讨论震源机制的异同，需要考虑数据的整体完备性。

利用谱聚类方法对震源机制解进行聚类分析，其具体步骤如下：

(1) 对震源机制数据按地震事件的时间顺序进行排列，并标注序号，用以建立震源机制数据的索引。

(2) 构建相似矩阵和度矩阵。利用Kagan算法求出两两地震对的最小旋转角，组成一个n×n的对称矩阵W(n为区域内震源机制解的数目)。这里的对称矩阵是由无向图的性质所决定。根据全连接方法中高斯核函数的定义，对距离矩阵进行转换，得到同样对称的相似矩阵。这里规定核参数σ=40，可以保证震源机制数据之间的权重值处于0~1区间内。同样，度矩阵也可以通过高斯核函数求得的权重值计算得到。

(3) 利用相似矩阵和度矩阵，计算得到拉普拉斯矩阵，然后对其标准化，得到标准化后的拉普拉斯矩阵L_sym(n×n)。

(4) 计算标准化拉普拉斯矩阵最小的k个特征值以及对应的特征向量f。将各自对应的特征向量f组成矩阵形式，然后按行进行标准化，最终组成n×k维的特征矩阵F。这时，震源机制解转换为了一组具有k维的数据矩阵。

(5) 利用传统的k-means聚类方法，对k维的n个新数据进行聚类，得到k个类别，利用序号索引还原到这k个类别中，可以得到最终的聚类结果。

在聚类过程中，采用GS方法分别对两个区域震源机制解的聚类结果进行评价，获取最佳聚类个数，如图 3所示。可以看到当聚类数k均为7时聚类效果最佳。比较两个区域数据聚类时的Gap值和误差条显示，海城地震数据Gap值较小且误差范围较大，这主要是因为地震样本量过少导致的；川滇及周边地区地震样本量较多，且地震分布范围较广，类间相似程度差异大，聚类结果较为准确且稳定。从以上两个区域地震数据Gap值分布状态分析可知，地震样本数量直接影响Gap值的误差分布，这主要是由震源机制的特性所决定的。将最佳聚类数代入震源机制谱聚类的计算步骤中，可以得到最终的聚类结果。下面将着重对聚类结果进行研究和分析。

3 结果分析 3.1 1975年海城M_S7.3地震序列

表 1给出了1975年海城M_S7.3地震序列的地震目录和震源机制解(顾浩鼎等，1976)。根据上文提到的方法和步骤，对海城区域24个震源机制解聚类分析，计算得到7类震源机制解，聚类结果见表 1。刁桂苓等(1992)曾利用系统聚类中的最长距离法将该区域24个震源机制划分为了3种类型，其中，系统聚类得到的第Ⅰ类解包含于本文的第Ⅰ类解，第Ⅱ类解对应了本文第Ⅰ、Ⅱ和Ⅲ类解，第Ⅲ类解对应于本文第Ⅳ、Ⅴ和Ⅶ类解。观察两个方法得到的聚类结果，主要有以下3点不同：

(1) 序号20的震源机制解其聚类结果间存在明显差异。系统聚类法将其归于前震相似的震源机制类别(第Ⅰ类)，而谱聚类方法将其划分为一个新的类别(第Ⅵ类)。通过分析相似度矩阵可知，谱聚类结果中第Ⅰ类解类内平均最小旋转角为6.1°，而序号20的震源机制解和第Ⅰ类平均解最小旋转角为39°，两者相差甚远。因此将序号20的震源机制解划分为第Ⅰ类解是不合适的。同样情况存在于序号为19的地震，经分析相似度矩阵，其与第Ⅰ类解距离更近。出现这种偏差可能是由于系统聚类采用的样本距离仅考虑了P、T轴，虽然简化了计算，但是对于空间的震源机制解数据未做到更好的约束。

(2) 谱聚类方法将系统聚类中的第Ⅱ类解进一步细分为两类，如表 1和图 4所示。为了突出每一类震源机制的特征，利用P轴和T轴的矢量和计算得到了每一类震源机制的平均解，图 4中用黑色圆弧和大实心圆来表示。从图 4中可以清晰看到，两类解还是有较大的差异，谱聚类方法得到的第Ⅱ类解地震为纯走滑型，P轴和T轴倾伏角较小，接近于水平；而第Ⅲ类解地震具有一定的倾滑分量，P轴和T轴倾伏角则明显变陡。谱聚类得到的结果较于系统聚类要更为精确一些。

(3) 顾浩鼎等(1976)和刁桂苓等(1992)提到的5个“错动过头”的地震在谱聚类中也区分了出来，而且进一步细分了这5个地震的类型(第Ⅳ、Ⅴ和Ⅶ类)，如图 4所示。三类解体现出明显的差异性，其中第Ⅳ和第Ⅴ类解具有明显的正断层倾滑性质，可能为“错动过头”所致，也可能受不同的区域构造背景所控制，如葠窝水库的垂向水重力影响。如果简单地归于一类解，体现不出部分震源机制的差异性，不利于区域内复杂地震活动的力学成因分析和研究。

通过上述分析可以发现，谱聚类方法在处理复杂震源机制解聚类过程中所得结果较好，不仅给出了类间震源机制的差异，而且利用GS法给出了聚类的最优解。但从图 3的Gap值波动范围可知，由于海城研究区仅有24个地震事件，样本量过少，因此结果有一定的误差，但无论如何，震源机制解的谱聚类方法得到的聚类结果仍是较为准确且可靠的。

3.2 川滇及周边区域

川滇地区位于青藏高原的东南缘，是近年来中国大陆地震活动最强烈的地区之一(闻学泽等，2003)。该区域范围较广，其构造活动以大规模剪切变形为主，兼有强烈的隆升运动，发育多条活动大断裂，构造背景非常复杂(罗钧等，2014)。对上文提到的216个震源机制解进行聚类分析，得到聚类结果空间分布(图 5)，图中震源机制解的聚类结果用不同颜色进行区分，每种类型的平均解如图例所示。

根据研究区内主要断裂带，简要分析和讨论聚类结果的空间分布特征：

金沙江断裂带主要分布两种类型的震源机制解：第Ⅳ类正断型和第Ⅶ类走滑型。其北段与理塘断裂带交汇区域以及南段与德钦—中甸—大具断裂带交汇区域均发育为正断型的第Ⅳ类解，而该断裂带东侧的理塘断裂带中段地震也为第Ⅳ类，说明两者的震源特征相同，这与罗钧等(2014)给出的结论一致。第Ⅶ类解多分布于金沙江断裂带以北区域，受构造活动影响走滑特征显著。鲜水河断裂带北段和中段同样分布第Ⅶ类解，但从甘孜往南，震源机制解变为第Ⅵ类。尽管两类解均为走滑型，但是分属两个不同的类别，说明鲜水河断裂带南北段受到不同的构造背景影响，差异特征明显，这在前人研究中得到了验证(罗钧等，2014；吴萍萍等，2016)。

龙门山断裂带主要分布第Ⅲ类逆冲型解，这主要受青藏高原隆起、深部物质向东流动的影响，形成“挤压”型地震。但由于该区域应力环境十分复杂，由图 5可以清晰地看到自SW至NE方向震源机制解类型发生多次变化，分段特征显著(杨宜海等，2021)。荥经—马边—盐津断裂带以逆冲型的第Ⅴ类解为主，和龙门山断裂带西南段地震形成原因相似。但由于两条断裂带断层走向不同，因此两个区域内震源机制解应分属于不同的类型，聚类结果准确且合理。

安宁河—则木河—小江断裂带及周边主要分布走滑型的第Ⅵ类解，其NS向节面的左旋走滑特征与该区域地质活动一致(段梦乔等，2019；苏珊，2019)。在巧家附近可以明显看到两个不同于附近震源机制解类型的地震，主要表现为正断型和正断兼走滑分量的特征，说明该区域两次地震均受到一定拉张作用的影响(罗钧等，2014)。谱聚类方法准确地将这两个地震区分出来。昭通-莲花山断裂带西南区域分布着与小江断裂带一致的第Ⅵ类解，但在东北区域受坚固的华南块体阻挡，地震断层类型变为NE向节面右旋走滑且部分带有逆冲分量，在更靠近四川盆地交接地时，震源机制解变为第Ⅲ类逆冲型，体现出了明显的分段特征(李君等，2019)。

丽江—小金河断裂带中北段南侧靠近安宁河断裂带主要分布走滑型的第Ⅵ类解，而北侧靠近德钦—中甸—大具断裂带主要分布正断型的第Ⅰ类解。两类解以丽江—小金河断裂带为界分隔开，呈现明显的空间分布特征，说明丽江—小金河断裂带中北段两侧受不同构造动力作用的影响。丽江—小金河断裂带南段分布有三种类型震源机制解，其中第Ⅰ类为该断裂带主要分布的震源机制解，第Ⅱ类为西南端大盈江和瑞丽断裂带主要分布的震源机制解，第Ⅵ类为东侧滇中次级块体主要分布的震源机制解，丽江—小金河断裂带南段夹于三个区域交汇处，体现出该段断层活动和发震断层破裂的复杂性(罗钧等，2014；Tian et al，2019)。

总体上看，川滇及周边区域震源机制解的聚类类型分区性明显，其空间分布符合各个断裂带的构造活动和块体间相互作用的控制，聚类结果准确可靠。

4 结论与讨论

本文利用谱聚类方法对震源机制解数据进行了聚类分析，该方法针对震源机制解的数据结构特点，采用震源机制之间的空间最小旋转角作为相似度矩阵，利用GS方法得到的最优解约束聚类的数量，基本解决了震源机制复杂数据解集的聚类问题。在谱聚类方法研究中，主要得到以下认识：

(1) 震源机制解的数据结构复杂多维，仅从单个或几个元素或者直观观察很难分辨出震源机制解之间的异同。谱聚类是近年来较为流行的一种聚类方法，针对多维数据有着良好的适应性，这对于解决多维多属性的震源机制解无疑是有利的。通过模拟数据验证，可以观察到不同类别之间明显的差异性，说明谱聚类方法针对震源机制解的聚类结果是可信的。

(2) 相似矩阵的构建是谱聚类算法的前提，震源机制解表达的是一种空间球体的概念，他们之间的最小空间旋转角可以真实反映两者之间的亲疏远近，因此采用最小空间旋转角表达震源机制解之间的差别最为科学(万永革，2019)。

(3) 本研究中的构图方式采用的是全连接法，其可以真实反映每个震源机制解的差异性，适用于区域性震源机制解的聚类，即样本量不高的情形。但如果样本过大，此方法迭代速度会受到严重制约，此时可考虑$ \epsilon$-邻近法，K邻近法等，尽管可能会造成部分信息缺失，但计算速度会大幅提升。

(4) 切图准则是谱聚类的核心思想，在现阶段应用方法很多，每个方法都有其解决NP难问题的优势范围。本文采用的NCut方法应用较为广泛，切图结果可靠性也较高，其归一化的特征可以平衡类间的相异度，避免对小区域的分割。针对震源机制解，NCut方法有其合理的解决范围，但其他切图方法或许更为适用，这需要在以后的实践中具体验证。

(5) 解决聚类数目问题的GS方法，其本质是统计样本类内和类间的离差平方和，通过统计量刻画样本观察值和参考分布下期望值之间的差异。这种复杂的大量统计方法可靠性高，应用在小样本中可以快速实现，但对于大样本的统计其时间成本过高，还需要高性能硬件的支撑。现阶段确定聚类数目方法很多，效果也因人而异，如何挑选适用于自身数据集的方法仍是当前聚类应用的一个难点工作，这需要更为深入的研究，从而寻找解决问题的最优解。

本研究不仅选择利用小区域小数据量的震源机制解数据进行谱聚类分析，还对大区域大数据量的震源机制解集进行了聚类分析。研究结果表明，相较于系统聚类方法，震源机制解的谱聚类方法在实际地震中更为适用，同时，GS方法也对聚类数做出了合理的判断，聚类结果更为准确且符合实际；在大数据量处理上，震源机制解的谱聚类方法不需要考虑震源机制类型，如Zoback(1992)的震源机制解分类方法，也可以合理准确地划分出震源机制解的不同类型，应用在川滇及周边地区中可以明显看到其聚类效果。由此可见，震源机制解的谱聚类方法不仅适用于小区域小数据量的地震事件，还适用于大区域大数据量的地震事件，且聚类结果基本准确且可靠。但在实际应用过程中，震源机制解的谱聚类方法仍存在一定的局限性，如川滇及周边地区的安宁河断裂带和曲江断裂带，断层实际走向差别巨大，但由于地震的震源机制很相似，谱聚类方法将其聚为同一类，说明分类结果的解释仍然需要结合断层实际走向，否则可能出错。解决该问题是改进谱聚类方法及今后发展的一个方向。