企业官网建设流程全解析

简介：面向材料力学与工程领域研究者的Matlab实现，基于RA-AF特征的高斯混合模型（GMM）用于识别拉伸裂纹与剪切裂纹，重点在于无需预先划分裂纹分界线，借助EM迭代算法完成聚类，直接输出两类裂纹分布图、统计数据及各样本所属类别的概率值。资源包仅含1个PDF文件，体积约45KB，内容为完整可运行的Matlab代码、示例数据生成逻辑及逐段注释，涵盖GMM参数初始化、E步骤/M步骤推导、收敛判断与结果可视化，适合正在学习聚类算法或开展声发射裂纹模式识别的读者参考。其中包含演示用的RA-AF数据生成逻辑，替换为实测数据即可运行；目前已有62人学习浏览此资源。代码结构清晰，可帮助理解EM迭代的数值实现细节，并为后续改进RA-AF判据或扩展多类型裂纹识别提供基础模板，教学演示或科研预实验均可快速上手。

1. 裂纹模式识别为什么绕开分界线：先看GMM聚类能做什么

在声发射检测里，RA-AF图是识别拉伸裂纹和剪切裂纹最常用的二维判据。传统做法是人工画一条分界线，落在哪边就归哪类。只要两类数据重叠得厉害，这条线就变成玄学：线往左偏，误判一批；线往右偏，又误判另一批。基于RA-AF的高斯混合聚类裂纹模式识别则完全绕开分界线——它用GMM拟合两类裂纹的数据分布，再用EM迭代算法从无标签数据里反推出均值和协方差，最终给出的是每个样本属于拉伸或剪切裂纹的概率值，而不是一个生硬的是非判断。这套Matlab代码从数据生成、EM迭代到分布图与统计输出全部包含，适合正在做声发射信号处理的研究生，也适合做材料无损检测的工程师拿去做模式识别的起点模板。

2. RA-AF特征与GMM模型：从“画线分类”到“概率归属”的切换

2.1 RA-AF特征：声发射信号里最关键的两个参数是什么

RA 是上升时间与幅值的比值：上升时间指声发射波形从起始到峰值经历的时间，幅值是波峰振幅，放在一起表示波形上升的陡峭程度。AF 是平均频率：振铃计数除以持续时间，表示波形在时间轴上的疏密程度。拉伸裂纹产生的声发射波通常比较“尖锐”，上升快、频率高，对应的 RA 值偏低、AF 值偏高；剪切裂纹的波形相对“平缓”，上升慢、频率低，RA 值偏高、AF 值偏低。把这两个参数放到二维平面，两类裂纹会各自聚成一片区域。问题是区域之间有重叠，尤其在实际构件中，裂纹混合出现时叠加噪声，重叠部分比想象中大得多。

RA 和 AF 的原始数据来源一般是声发射采集系统导出的通道数据，每个事件一条记录。做聚类之前，先要把原始波形转成这两个标量参数。常见做法是在采集软件里直接设置 RA、AF 的导出公式，也可以在 Matlab 里自己按公式算：

% 计算RA和AF的常见写法 % rt：上升时间(us)，amp：幅值(dB)，counts：振铃计数，duration：持续时间(us) RA = rt ./ (amp * 1e-3); % 上升时间/幅值，单位us/mV AF = counts ./ (duration * 1e-6); % 振铃计数/持续时间，单位kHz

这里的单位换算按采集系统常用量纲走：幅值从 dB 换算成 mV，持续时间从微秒换算成秒。如果你的采集系统已经直接给出 RA、AF，那这两行就不用重复算，但要注意单位——不少系统给的是 us/V 而不是 us/mV，差 1000 倍聚类结果完全不同。换算好之后，把两列数据合并成一个 N×2 矩阵，这就是后续 GMM 的输入。

2.2 GMM为什么适合裂纹模式识别：混合分布的数学描述

GMM假设数据由 K 个高斯分布混合而成。放在裂纹识别场景里，K 通常取 2，对应拉伸裂纹和剪切裂纹。每个高斯分量有自己的均值向量、协方差矩阵和权重，权重表示该裂纹模式在总体中占的比例。这样建模的好处是：它不要求两类数据有明确边界，重叠区域由两个高斯分布的密度共同覆盖，归属关系由概率而非硬阈值决定。

从数学上看，GMM 的概率密度是多个高斯分量的加权叠加：

p(x) = Σ_k π_k · N(x | μ_k, Σ_k)

这里 π_k 是第 k 分量的混合权重，满足所有 π_k 之和等于 1；μ_k 是均值向量，决定该分量的中心位置；Σ_k 是协方差矩阵，决定该分量在二维平面上的椭圆形状和方向。对裂纹识别来说，拉伸裂纹和剪切裂纹各自的数据分布各有倾向：拉伸裂纹的椭圆可能偏向 RA 小、AF 大的方向，剪切裂纹则反过来。两个椭圆重叠的部分就是分类模糊的区域，这也是人工分界线最头疼的地方——分界线只能线性切开平面，而两个椭圆的重叠不一定是线性可分的。GMM 用两个二维高斯函数直接描述这种非线性的分布形态，重叠区域由两个分量共同分配概率，这是它在物理直观上比线性分界线更贴合的一层原因。

另一个重要点在于，GMM 拟合的并不是“哪个样本属于哪一类”，而是“这些样本整体服从什么样的混合分布”。一旦分布参数确定，任何一个新样本都可以代入公式计算它属于每个分量的后验概率。这就给后续的在线监测留了接口——可以先离线训练好 GMM，再对实时获取的声发射事件算概率，不需要每次重新聚类。

2.3 EM算法：在没有标签的情况下如何反推混合参数

GMM 的似然函数没有解析解，因为每个样本不知道来自哪个高斯分量，这就是所谓的“隐变量”问题。EM 算法的核心思路是：先猜一组参数，根据这组参数给每个样本分配属于各分量的概率（E步），再基于这些概率重新估计参数（M步），反复迭代直到参数变化小于阈值。E步算出的责任度就是每个样本属于每个分量的后验概率，M步相当于做一次带权重的最大似然估计。

展开来看，E步做的事情是：固定当前均值、协方差、权重，对每个样本分别计算它在每个高斯分量下的概率密度，乘以对应权重后再归一化。归一化后的值就是责任度，取值范围在 0 到 1 之间，表示这个样本有多大概率属于该分量。M步做的事情是：把责任度当作软标签，重新计算每个分量的加权均值、加权协方差和权重。加权均值就是所有样本按责任度加权后的平均位置，责任度高的样本贡献更大。

EM 迭代的关键属性是：每次迭代都会让对数似然函数单调不减。也就是说，理论上它保证收敛，但不保证收敛到全局最优。初始参数给得不合适，EM 很可能收敛到局部最优，所以初始化策略在工程里非常关键。常见做法是用 kmeans 先做一次粗糙的硬聚类，把聚类中心和类内散度作为 EM 的初值，这样比纯随机初值稳定得多。

2.4 和人工分界线方法的对比：三个关键差距

人工分界线是最早出现的方法：在 RA-AF 散点图上画一条直线或折线，线一侧算拉伸裂纹，另一侧算剪切裂纹。优点是直观，缺点是重叠区域的误判率几乎无法控制。有人提出用 SVM 代替人工画线，但 SVM 需要带标签的训练数据，在工程现场给每个声发射事件标注类型成本不低，而且标注本身的可靠性也是个问题。

GMM + EM 的最大不同在于，它不关心分界线在哪里，只关心数据分布拟合得有多好。下表整理了三种方法的差异：

一个容易被忽视的细节是，GMM 并不追求“分得越干净越好”，它追求的是“数据分布拟合得越准越好”。分布拟合得准，即使两类数据重合程度很高，每个样本的概率输出仍然有参考价值；而人工画线一旦画在密度重叠区域，结果基本没法稳定复现。同一批数据，不同工程师画出来的线可能差很远，GMM 在给定初始化策略和分量数后，结果是可重复的。

3. Matlab代码逐段拆解：数据生成、参数初始化与EM迭代核心

3.1 数据生成：构造二维RA-AF样本集

在实际项目中，RA-AF 数据来自声发射采集系统。为了演示聚类流程，原代码用高斯随机数构造了两类裂纹的二维样本：

num_samples = 200; % 每类样本数量 % 拉伸裂纹：RA偏低，AF偏高（代码里用均值为正的第一维模拟） tensile_data = [randn(num_samples, 1)*0.5 + 1; randn(num_samples, 1)*0.5 + 2]; % 剪切裂纹：RA偏高，AF偏低 shear_data = [randn(num_samples, 1)*0.5 - 1; randn(num_samples, 1)*0.5 - 2]; % 合并为二维特征矩阵：每行是一个样本，列分别对应RA和AF data = [tensile_data, shear_data];

这段代码构造的是一个 400 行 2 列的数据矩阵。第一列模拟 RA 方向的特征值，第二列模拟 AF 方向的特征值，每行代表一个声发射事件。0.5是标准差，控制类内离散程度；均值偏移控制两类之间的间距。实际使用中，data应该换成从声发射系统导出的真实 RA-AF 矩阵，格式保持每一行一个事件、两列分别对应 RA 和 AF 即可，GMM 不需要额外信息。

需要注意这里的写法容易造成误解：tensile_data先纵向拼接了两批各 200 行的高斯样本，之后用方括号横向合成 400 行 2 列。真正决定样本类别的是两列的组合关系，而不是单列数值。如果你刚开始接触这段代码，可以先记住“行是事件，列是特征”这个约定，后面避坑章还会展开讲。

3.2 GMM参数初始化：权重、均值、协方差怎么给初值

num_components = 2; % 高斯分量数，对应拉伸和剪切两种裂纹 max_iterations = 100; % 最大迭代次数 tolerance = 1e-6; % 收敛阈值 weights = ones(num_components, 1) / num_components; % 初始权重，各0.5 means = randn(num_components, 2); % 初始均值：2行2列，每行是一个分量的中心 covariances = repmat(eye(2), [1, 1, num_components]); % 初始协方差：2x2x2的三维数组

分量数设为 2，对应拉伸和剪切两类裂纹。如果实际损伤模式不止两类，比如出现混合裂纹，可以改成 3 或更多，但输出统计的代码要同步调整。初始权重取均匀分布是最保守的选择，它不引入任何先验偏见。初始均值用randn随机生成，这是原代码里风险最高的一行——随机初始值最容易让 EM 收敛到局部最优，后面避坑章会给出替代方案。初始协方差用单位矩阵，对二维特征来说，每个分量的初始分布是一个半径为 1 的正圆，这个选择在没有任何先验信息时是合理的。

repmat(eye(2), [1, 1, num_components])生成的是一个 2×2×2 的三维数组，第三维索引对应分量编号。用covariances(:, :, k)访问第 k 个分量的协方差矩阵，这种存储结构在代码里很常见，不要把它当成普通的二维矩阵拼接。

3.3 E步：责任度计算与归一化

for iter = 1:max_iterations num_samples = size(data, 1); responsibilities = zeros(num_samples, num_components); for k = 1:num_components % 每个样本在第k个高斯分量下的概率密度 responsibilities(:, k) = mvnpdf(data, means(k, :), covariances(:, :, k)); end % 乘以混合权重 responsibilities = responsibilities .* weights'; % 按行归一化，使每个样本对所有分量的责任度之和为1 responsibilities = responsibilities ./ repmat(sum(responsibilities, 2), [1, num_components]);

E 步的关键是mvnpdf函数，它计算每个样本在指定均值、协方差下的多元高斯概率密度。mvnpdf(data, means(k, :), covariances(:, :, k))返回一个 N×1 的列向量，对应 N 个样本在第 k 个分量下的密度值。三个分量的密度值并列存放后，乘上对应权重，再按行归一化，得到的就是责任度矩阵，维度 N×K，每一行之和为 1。

weights'这里把 K×1 的权重向量转置成 1×K，然后利用 Matlab R2016b 以后引入的隐式扩展机制和 responsibilities 做逐元素乘法，效果是每一列乘上对应分量的权重。在旧版本 Matlab 上，这类操作必须显式写成repmat(weights', [num_samples, 1])，否则会报尺寸不匹配。如果读者用的是老版本，建议统一写成responsibilities .* repmat(weights', [num_samples, 1])更稳妥。

3.4 M步：均值、协方差和权重的更新

% M步骤 Nk = sum(responsibilities, 1); % 各类别的有效样本数（软计数） weights = Nk / num_samples; % 更新混合权重 for k = 1:num_components % 更新均值：按责任度加权平均 means(k, :) = sum(responsibilities(:, k) * data, 1) / Nk(k); % 更新协方差 diff = data - repmat(means(k, :), [num_samples, 1]); covariances(:, :, k) = (diff' * (responsibilities(:, k) * diff)) / Nk(k); end

M 步的核心是软计数：Nk不是整数值，而是每个类别责任度之和，表示“等效样本数”。比如拉伸分量的 400 个责任度加起来等于 250，那就相当于有 250 个等效样本属于拉伸裂纹，剩下 150 个等效样本被分配给剪切分量。权重更新就是拿这个等效样本数除以总样本数。

均值更新的写法sum(responsibilities(:, k) * data, 1) / Nk(k)，在 Matlab 中responsibilities(:, k)是 N×1 列向量，data是 N×2 矩阵，两者直接相乘会报维度错误。正确的做法是用逐元素乘法：

% 正确写法：责任度加权求和 means(k, :) = sum(responsibilities(:, k) .* data, 1) / Nk(k);

协方差更新那一行同样存在问题，responsibilities(:, k) * diff在矩阵乘法规则下不成立，应该改成responsibilities(:, k) .* diff。这是原代码里最容易直接报错的位置，很多同学把项目下载下来第一跑就卡在这里。避坑章会再展开一次。

3.5 收敛判据：只看均值变化够不够

% 检查收敛 if iter > 1 && norm(means - old_means) < tolerance break; end old_means = means; end

收敛判断只比较了新旧均值的 Frobenius 范数变化。实际工程中，协方差和权重在迭代后期也会有小幅抖动，尤其在数据量少的时候。更稳妥的做法是把均值、协方差、权重拼接成一个向量，统一计算变化量：

% 更全面的收敛判断 param = [means(:); weights(:)]; for k = 1:num_components param = [param; covariances(:, :, k)(:)]; end if iter > 1 && norm(param - old_param) < tolerance break; end old_param = param;

不过从工程角度说，均值是最敏感的指标，均值一旦稳定，协方差和权重也基本到稳态。如果你看到迭代次数跑到max_iterations才有输出，先不要急着加迭代次数，重点检查协方差矩阵是否退化、数据是否标准化，这两点原因占比最高。

4. 结果输出与可视化：分布图、概率值与统计报告的整理

4.1 类别归属：从责任度到离散标签的转换

% 输出每个数据所属类别的概率值 probabilities = responsibilities; % 确定每个数据点所属的类别 [~, labels] = max(responsibilities, [], 2);

max(responsibilities, [], 2)沿第二维取最大值，labels是 N×1 的列向量，每个值是 1 或 2，表示该样本在哪个分量下责任度最高。probabilities直接保存了责任度矩阵，这是后续做置信度分析的基础数据。

需要特别提醒的是：聚类标签本身是任意的。第一次跑可能分量 1 对应拉伸裂纹，第二次跑因为随机初始化，分量 1 变成了剪切裂纹。两个分量没有先验语义，想让labels和物理裂纹类型对应，必须在后面加一个映射步骤。常见做法是聚类结束后比较两个分量的均值：RA 值较低、AF 值较高的分量定义为拉伸裂纹；RA 值较高、AF 值较低的定义为剪切裂纹。这一步不做的话，统计输出就是错误的。

4.2 分布图绘制：gscatter的参数细节

figure; gscatter(data(:, 1), data(:, 2), labels); title('拉伸裂纹和剪切裂纹分布图'); xlabel('RA'); ylabel('AF');

gscatter是 Matlab 自带的分组散点图函数，第三个参数labels是分组依据，每一组自动分配不同颜色和符号。这里data第 1 列是 RA，第 2 列是 AF，x 轴对应 RA，y 轴对应 AF。如果真实数据的 RA 和 AF 数量级相差很大，建议先标准化再画图，或者把坐标轴改成对数刻度。

gscatter默认会自动生成图例，对两类样本来说图例没问题，但如果分量数增加到 5 类以上，图例会占掉很大一块画布。可以用'legend', 'off'参数关闭默认图例，然后在图上手动标注类别含义。另外一个经常被忽略的点是：gscatter的符号形状由类别索引决定，不是由你的物理类型决定。假如分量 1 实际上是剪切裂纹，图上圆形点就是剪切裂纹，图例文字要按实际对应关系写，不能想当然写成拉伸。

4.3 统计输出：拉伸裂纹、剪切裂纹数量和比例的解读

tensile_count = sum(labels == 1); shear_count = sum(labels == 2); fprintf('拉伸裂纹数量: %d\n', tensile_count); fprintf('剪切裂纹数量: %d\n', shear_count);

这段统计直接依赖前面标签与裂纹类型的对应关系。如果分量 1 是剪切裂纹、分量 2 是拉伸裂纹，打印结果就会颠倒。所以代码看起来逻辑清晰，实际上藏着一条隐藏前提——必须先确认分量顺序。

一个更稳妥的做法是先打印两个分量的均值，人工确认后再输出统计：

disp('各分量均值（第一列RA，第二列AF）：'); disp(means);

看到均值之后，再决定labels == 1对应的是拉伸还是剪切。对自动化流程来说，可以把判定规则写进代码：比较means(1,1)和means(2,1)的 RA 大小，RA 小的那个分量定义成拉伸裂纹，然后重新映射标签。这样整个流程不需要人工介入，也能保证统计口径一致。

4.4 概率值到哪里去了：后续分析能做什么

概率值是 GMM 比传统画线法多出来的信息。每个样本都有一对概率值，比如 [0.8, 0.2] 表示 80% 属于拉伸、20% 属于剪切。对处于重叠区域的样本，这份信息比单纯的标签更有分析价值。你可以在聚类完成后统计置信度低于某个阈值的样本比例，比如把所有最大概率小于 0.7 的样本标记为“低置信度”，这些样本往往是裂纹模式转变期或者噪声干扰较大的事件。

另一个常见用途是时序分析。把声发射事件按时间窗口分组，每组内统计每类裂纹的概率均值或数量占比，然后绘制损伤演化曲线。典型趋势是加载初期拉伸裂纹占比高，中后期剪切裂纹逐渐增多。GMM 给的概率值可以天然地做加权统计，比硬标签的 0/1 计数平滑得多，噪声事件对曲线的扰动也更小。

5. EM聚类实战避坑：五条容易翻车的细节记录

5.1 现象：M步报“矩阵维度必须一致”

运行原代码时，在协方差更新那一行经常会直接抛错，报错信息类似“矩阵维度必须一致”。原因很明确：responsibilities(:, k) * diff里，responsibilities(:, k)是 N×1 列向量，diff是 N×2 矩阵，矩阵乘法要求左边矩阵的列数等于右边矩阵的行数，这里 N 不等于 2，直接报错。原代码把矩阵乘法和逐元素乘法的概念混用了。

解决方法是把内部的乘号改成点乘：

covariances(:, :, k) = (diff' * (responsibilities(:, k) .* diff)) / Nk(k);

这里.*表示逐元素乘法，责任度向量的每个元素对应乘到diff的每一行上，相当于给每个样本的偏差向量加权。如果用的是老版本 Matlab，建议用bsxfun(@times, responsibilities(:, k), diff)代替；新版本里点乘加隐式扩展已经够用。

5.2 现象：迭代到上限不收敛，协方差矩阵出现奇异

EM 迭代次数跑满 100 次仍然没有触发收敛条件，或者某次迭代时协方差矩阵变成奇异矩阵，mvnpdf计算出 NaN。原因有几个方向：一是数据量太少，某一分量的等效样本数 Nk 趋近于 0，协方差更新时分母为 0；二是数据本身存在大量重复或近似重复的样本，导致协方差矩阵不可逆；三是特征维度里存在一个方差远大于另一个的维度，协方差矩阵条件数过大。

解决办法分两步。第一步，给协方差矩阵加正则项：

covariances(:, :, k) = covariances(:, :, k) + 1e-6 * eye(2);

这个操作相当于给协方差矩阵对角线加上一个很小的正数，防止矩阵奇异，代价是引入一点偏差，但对聚类结果影响可忽略。第二步，检查数据质量。我一般会先跑一次sum(isnan(data))和sum(isinf(data))，把异常值清理干净再重新聚类。数据量太少时不要强行用 GMM，减少分量数或者增加样本量才是根本解决办法。

5.3 现象：两个高斯分量最终合并成一个

聚类结果里两个分量的均值非常接近，协方差椭圆几乎重叠，输出的分布图看起来像只有一簇数据。这种现象的本质是 EM 收敛到了局部最优。初始均值是随机生成的，如果两个初始点落在同一条数据流形的同一侧，EM 会沿着这条路径收敛到“一个分量吃掉两类数据”的局部解。

解决方法是改用 kmeans 初始化。先对数据做一次硬聚类，把聚类中心作为means的初值，把每个簇的样本协方差作为covariances的初值：

% 更稳健的初始化：先用kmeans预聚类 [idx, c] = kmeans(data, num_components); means = c; for k = 1:num_components cluster_data = data(idx == k, :); if size(cluster_data, 1) > 1 covariances(:, :, k) = cov(data); else covariances(:, :, k) = eye(2); end end

我还习惯做多组初始化对比：把 EM 跑 5 到 10 次，每次用不同的随机种子，取对数似然最大的那一次。这个做法不复杂，但对结果稳定性的提升非常明显。

5.4 现象：聚类结果把拉伸和剪切裂纹标签对调了

第一次运行得到的结果里拉伸裂纹对应分量 1，第二次运行变成分量 2，统计输出的数量完全反了。这不是代码 bug，而是 GMM 的固有特性：分量索引是模型内部约定，没有物理意义绑定。任何一次初始化不同，分量编号都可能交换。

解决办法是固定映射规则。聚类结束后比较两个分量的均值：

% 固定映射：认为RA较小、AF较大的分量为拉伸裂纹 if means(1, 1) < means(2, 1) % 分量1的RA更小 tensile_label = 1; shear_label = 2; else tensile_label = 2; shear_label = 1; end labels_physical = zeros(size(labels)); labels_physical(labels == tensile_label) = 1; % 1表示拉伸 labels_physical(labels == shear_label) = 2; % 2表示剪切

如果项目里还有少量已知类型的样本，更直接的做法是用这些样本校验分量对应关系，然后固定映射。从那以后我做自动批量处理时，都会在聚类代码后强制加上这段映射逻辑，否则统计表里的数量随时可能反过来。

5.5 现象：RA和AF量纲差太大，聚类结果被单维主导

真实数据里 RA 可能是几十到几百的数量级，AF 可能只有几到几十，或者反过来。GMM 的协方差矩阵对不同维度的方差一视同仁，方差大的维度会在距离计算中占据绝对主导，另一个维度几乎不起作用，聚类结果等效于只按单变量切分。

解决办法是聚类前做 z-score 标准化：

mu = mean(data); sigma = std(data); data_norm = (data - mu) ./ sigma;

聚类结束后，把标签映射回原始数据坐标画图。注意标准化不会改变数据的相对分布结构，只是把两个维度的尺度统一了。如果在标准化后聚类结果有明显变化，说明原来某个维度确实被压制了，这个现象本身也值得记录到实验报告里。

6. 进阶验证：真实数据跑通、分量数选择与迁移建议

6.1 从合成数据迁移到真实声发射数据的处理流程

把代码里合成数据部分替换成真实数据，先要做一次清洗，顺序不能乱：第一步删除 RA、AF 为 NaN 或 Inf 的记录；第二步排除幅值接近 0 的异常事件，这类事件算出来的 RA 会大到离谱；第三步做 z-score 标准化；第四步喂给 GMM；最后用原始坐标画分布图。清洗这一步不能省，真实数据里一条幅值近 0 的记录可以把一个分量的协方差矩阵拉爆，EM 完全不收敛。

6.2 分量个数选择：从“知道两类”扩展到“不知道几类”

如果面对的是未知损伤模式，分量数不是拍脑袋定 2，而是要做选择。常见做法是遍历 K=1 到 K=6，分别计算 BIC，选 BIC 最小的 K。BIC 对模型复杂度有惩罚项，在样本量不确定时比 AIC 更保守，不容易出现过拟合。

% 用BIC选择GMM分量数，代码示意 for k = 1:6 gmm_model = fitgmdist(data_norm, k, 'RegularizationValue', 1e-5); bic(k) = gmm_model.BIC; aic(k) = gmm_model.AIC; end [~, best_k] = min(bic); fprintf('BIC选择的分量数: %d\n', best_k);

自己实现 EM 时，每跑一个 K 要记录最终的对数似然，再套 BIC 公式。注意 BIC 里的参数个数不是 K 本身，而是所有自由参数的总数。二维 GMM 下每个分量包含均值向量 2 个、协方差矩阵的上三角元素 3 个、权重 1 个，共 6 个参数，再加上权重和为 1 的约束，实际自由参数是 6K-1 个。很多初学同学在这里算错，导致 BIC 曲线完全失真。

6.3 最终验证：聚类结果与裂纹实际形貌对应

聚类完成后，建议把自动划分的结果和样本的原始波形时间戳对齐，查看拉伸裂纹和剪切裂纹的出现顺序。在一段典型的加载过程中，早期往往拉伸裂纹占主导，中后期剪切裂纹比例上升。如果 GMM 聚类结果能复现这种趋势，说明分出来的两个模式有物理意义。验证时可以用时间窗口统计：

% 按时间窗口统计裂纹类型比例 window_size = 100; num_windows = floor(size(data, 1) / window_size); ratio_shear = zeros(num_windows, 1); for w = 1:num_windows seg = labels_physical((w-1)*window_size + 1 : w*window_size); ratio_shear(w) = sum(seg == 2) / window_size; end plot(ratio_shear);

这条思路能确认聚类不是纯统计产物，而是对应了裂纹演化过程。我踩过最深的坑就是只看聚类图、没看时间序列，结果把两段完全独立的传感器数据混在一起聚了类，得到漂亮的分布图但没有物理意义。从那以后我每次跑 GMM 聚类都会强制把时间顺序引入验证环节，至少画一张按时间分窗的裂纹类型比例曲线，确认趋势合理再下结论。希望这个流程能帮到你。

本文还有配套的精品资源，点击获取