简介这份资源面向材料力学、岩石断裂与工程检测方向的研究者及研究生提供一套基于RA-AF参数的高斯混合聚类裂纹模式识别Matlab代码无需人工划分裂纹分界线直接采用EM迭代算法完成GMM聚类用于区分拉伸裂纹与剪切裂纹。压缩包内仅含1个PDF文件约45KB集中呈现完整代码、逐段注释与实现思路说明便于对照理解算法流程。内容涵盖数据生成、GMM参数初始化、E步责任度计算、M步权重与均值协方差更新、收敛判断以及类别标签确定、拉伸与剪切裂纹数量统计和RA-AF分布图绘制并输出每个数据点所属类别的概率值。已有62人学习适合希望快速复现聚类识别流程、理解EM算法在裂纹模式分类中应用的读者参考。1. 从一张散点图说起为什么我放弃了画分界线去年帮某实验室处理一批声发射数据RA-AF 散点图往屏幕上一投拉伸裂纹和剪切裂纹的簇肉眼可见地叠在一起。按老办法得在图上手动画一条斜线线左上归拉伸、右下归剪切。问题是这条线的斜率、截距全靠经验换一批试件就得重调审稿人一句“分界依据是什么”就能把人问住。后来我改用高斯混合模型做聚类让 EM 算法自己去找两个高斯分量的均值和协方差不再需要人为划界。这份 MATLAB 代码干的就是这件事输入 RA、AF 两列数据输出每个点属于拉伸或剪切裂纹的概率、类别标签、数量统计和分布图。适合做材料力学、岩石破裂、混凝土损伤方向的研究生和工程师尤其是被“分界线怎么定”卡过的人。它不解决数据采集问题也不替你判断该分几类但只要你手头有 RA-AF 数据它能把分类这件事从“画线玄学”变成可复现的概率推断。2. RA-AF 特征与 GMM 的匹配逻辑为什么是这两个东西凑一起2.1 RA 和 AF 分别代表什么为什么它们能区分裂纹模式RA 是上升时间与幅值的比值单位通常是微秒每伏或类似量纲反映的是信号从起点爬到峰值有多“慢”。AF 是平均频率即振铃计数除以持续时间反映信号振荡的快慢。在声发射检测里拉伸型裂纹张开型通常表现为高 RA、低 AF因为裂纹张开过程释放的能量持续时间长、上升沿缓剪切型裂纹滑移型则相反低 RA、高 AF信号更“尖”、频率更高。这两个特征放在二维平面上理论上会形成两个朝向不同的椭圆簇。但实际数据没那么干净——试件几何、传感器耦合状态、背景噪声都会让两个簇部分重叠。传统做法是在 RA-AF 平面上画一条斜线线两侧分属两类。这条线的位置依赖具体实验条件换一个加载速率或试件尺寸就得重新标定可重复性差。2.2 GMM 为什么比硬分界更合适高斯混合模型的核心假设是整批数据由若干个高斯分布按不同权重混合而成。每个高斯分量有自己的均值向量和协方差矩阵均值决定簇的中心位置协方差决定簇的形状和朝向。对于 RA-AF 数据拉伸裂纹和剪切裂纹各自形成一个高斯分量两个分量的协方差矩阵可以不同——这意味着一个簇可以比另一个更“扁”或更“斜”更贴近真实分布。EM 算法交替执行两步E 步计算每个数据点属于每个分量的后验概率也叫责任度M 步用这些概率加权重新估计每个分量的均值、协方差和权重。迭代到参数不再明显变化为止。相比硬分界GMM 的输出是概率——每个点有 70% 概率属于拉伸、30% 属于剪切而不是非此即彼。这对后续做可靠性分析或阈值筛选很有用你可以只保留概率大于 0.9 的点做统计把模糊地带单独拿出来看。2.3 代码整体结构从数据到标签的完整链路这份代码的流程很直准备 RA-AF 数据矩阵初始化 GMM 参数跑 EM 迭代输出概率和标签画图统计。数据矩阵是 N×2第一列 RA第二列 AF。初始化时权重设为均匀均值随机取协方差设为单位阵。迭代中 E 步用mvnpdf算每个点在每个分量下的概率密度乘以权重后归一化得到责任度M 步用责任度加权更新参数。收敛判据是均值矩阵的范数变化小于阈值。最后用max取责任度最大的分量作为标签gscatter按标签着色画散点图。整个脚本不到 80 行但每一步都有明确的数学含义改参数的地方也集中适合在此基础上做扩展。3. 把代码跑起来数据准备、参数设置与逐段拆解3.1 数据导入与格式要求代码里默认用randn生成示例数据实际用时你得替换成自己的 RA-AF 矩阵。常见做法是从 Excel 或 CSV 读入确保两列都是数值型没有 NaN。如果 RA 和 AF 量纲差异大比如 RA 在 0-1000AF 在 0-500建议先做标准化否则协方差矩阵会被大量纲主导。标准化用zscore一行搞定但要注意标准化后均值接近 0协方差接近单位阵初始化时可以相应调整。下面这段是替换示例数据的写法% 从 CSV 读取 RA-AF 数据假设第一列 RA第二列 AF raw readmatrix(ra_af_data.csv); data raw(:, 1:2); % 只取前两列 % 检查缺失值 if any(isnan(data(:))) warning(数据中存在 NaN已删除对应行); data data(~any(isnan(data), 2), :); end % 可选标准化消除量纲影响 % data zscore(data);逻辑说明readmatrix是 MATLAB R2019a 之后推荐的读表函数比xlsread干净。isnan检查整矩阵any(...,2)找出有 NaN 的行并删除。标准化那行默认注释掉因为标准化会改变均值和协方差的物理含义如果你需要保留原始量纲做解释就别开。参数方面raw(:,1:2)假设前两列就是 RA 和 AF如果列顺序不同改索引即可。3.2 GMM 参数初始化分量数、迭代上限与收敛阈值代码里num_components 2对应拉伸和剪切两类。如果你不确定该分几类可以试 2 到 4看 BIC 或 AIC 哪个更小但这份代码没内置模型选择需要自己加。max_iterations 100是保险上限实际通常 20 到 50 步就收敛。tolerance 1e-6是均值变化的范数阈值太小会多跑无用迭代太大可能提前停。初始化权重用均匀分布均值用randn随机取协方差用单位阵。这里有个坑随机初始化可能导致每次跑的结果略有不同因为 EM 对初值敏感。常见做法是跑多次取对数似然最大的那次或者用 k-means 先粗分一下再拿簇中心当初始均值。下面是把初始化改成 k-means 引导的写法% 用 k-means 初始化均值减少 EM 陷入局部最优的概率 rng(42); % 固定随机种子保证可复现 [idx_init, C_init] kmeans(data, num_components, Replicates, 5); means C_init; % k-means 的簇中心作为初始均值 % 权重按各簇样本比例初始化 weights histcounts(idx_init, 1:num_components1) / size(data, 1); % 协方差仍用单位阵也可以按各簇样本协方差初始化 for k 1:num_components covariances(:, :, k) eye(2); end逻辑说明kmeans的Replicates参数让 MATLAB 跑 5 次取最优降低陷入差局部最优的概率。histcounts统计各簇样本数除以总数得到初始权重。协方差这里仍用单位阵保守但稳定如果你想更激进可以用cov(data(idx_initk,:))初始化但要注意小簇可能导致协方差奇异。rng(42)固定种子方便复现。3.3 EM 迭代核心E 步责任度计算与 M 步参数更新E 步和 M 步是 EM 的心脏。E 步用mvnpdf算每个点在每个分量下的概率密度注意mvnpdf的第二个参数是均值行向量第三个参数是协方差矩阵。算完概率密度后乘以权重再按行归一化得到责任度矩阵每行和为 1。M 步用责任度加权更新权重、均值和协方差。均值更新是责任度加权平均协方差更新是责任度加权的散布矩阵除以该分量的有效样本数。代码里Nk sum(responsibilities, 1)就是每个分量的有效样本数。收敛检查用norm(means - old_means) tolerance只比较均值没比较协方差和权重这是简化处理够用但不完美。下面把 E 步和 M 步单独拎出来加注释% E 步计算责任度 responsibilities zeros(num_samples, num_components); for k 1:num_components % mvnpdf 计算每个点在分量 k 下的概率密度 responsibilities(:, k) mvnpdf(data, means(k, :), covariances(:, :, k)); end responsibilities responsibilities .* weights; % 乘以混合权重 responsibilities responsibilities ./ sum(responsibilities, 2); % 按行归一化 % 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 - means(k, :); % 去均值 covariances(:, :, k) (diff * (responsibilities(:, k) .* diff)) / Nk(k); % 加权协方差 end逻辑说明mvnpdf要求协方差矩阵正定如果某个分量塌缩到几个点上协方差可能奇异mvnpdf会报错或返回 NaN。常见做法是给协方差加一个小的正则项比如covariances(:,:,k) 1e-6*eye(2)。sum(responsibilities, 2)是按行求和结果是一个列向量用./广播除法实现归一化。M 步里responsibilities(:, k) .* diff是每个点的责任度乘以其去均值向量再转置相乘得到加权散布矩阵。注意Nk(k)可能很小导致除零或数值不稳可以在更新前检查Nk(k) 1e-3就跳过该分量或重新初始化。3.4 结果输出概率、标签、统计量与分布图迭代结束后responsibilities就是每个点属于每个分量的概率。[~, labels] max(responsibilities, [], 2)取每行最大概率对应的分量索引作为硬标签。统计拉伸和剪切数量就是数标签 1 和 2 的个数。画图用gscatter按标签着色横轴 RA、纵轴 AF。如果你想看概率分布可以再用histogram画每个点最大概率的直方图概率接近 1 的点占多数说明聚类分离度好。下面加一段概率直方图和置信椭圆的代码% 输出概率统计 max_prob max(responsibilities, [], 2); fprintf(平均最大概率: %.3f\n, mean(max_prob)); fprintf(概率大于 0.9 的点占比: %.1f%%\n, sum(max_prob 0.9) / num_samples * 100); % 画概率直方图 figure; histogram(max_prob, 20); xlabel(最大后验概率); ylabel(频数); title(聚类置信度分布); % 在散点图上叠加 95% 置信椭圆 figure; gscatter(data(:,1), data(:,2), labels); hold on; for k 1:num_components [V, D] eig(covariances(:, :, k)); t linspace(0, 2*pi, 100); ellipse 2.447 * V * sqrt(D) * [cos(t); sin(t)] means(k, :); plot(ellipse(1,:), ellipse(2,:), k--, LineWidth, 1.5); end hold off; xlabel(RA); ylabel(AF); title(RA-AF 分布与 95% 置信椭圆);逻辑说明max_prob是每个点的最大责任度均值越接近 1 说明聚类越干脆。sum(max_prob 0.9)统计高置信点占比低于 70% 就得警惕两个簇重叠严重。置信椭圆用协方差矩阵的特征分解画2.447是二维 95% 置信对应的卡方分位数平方根。V * sqrt(D)是协方差矩阵的“平方根”用来把单位圆映射成椭圆。这段代码不改变聚类结果只是可视化辅助帮你判断结果可不可信。4. 避坑与排查那些让我重跑一整天的细节4.1 协方差矩阵奇异导致 mvnpdf 报错现象跑着跑着mvnpdf返回 NaN 或直接报“协方差矩阵必须为正定”。原因某个高斯分量只分到极少数点或者这些点几乎共线导致协方差矩阵接近奇异。解决在 M 步更新协方差后加正则项covariances(:,:,k) covariances(:,:,k) 1e-6 * eye(2)同时在 E 步前检查Nk(k)如果小于 3 就重新初始化该分量的均值和协方差或者直接减少分量数。4.2 随机初始化导致每次结果不一样现象同一批数据跑两次标签翻转了拉伸和剪切的编号对调。原因EM 对初值敏感随机初始均值可能让两个分量互换角色。解决固定随机种子rng(42)并用 k-means 初始化均值。如果标签编号对调不影响物理结论因为你可以根据均值判断哪个是拉伸那不用管但如果要写进报告最好在输出后按均值大小重新排序标签比如 RA 均值大的归为拉伸。4.3 量纲差异让协方差矩阵被一个特征主导现象RA 范围 0-1000AF 范围 0-10聚类结果几乎只按 RA 分AF 没起作用。原因协方差矩阵对大量纲特征敏感EM 在最大化似然时会被大量纲方向的方差主导。解决先对两列做zscore标准化让每个特征均值为 0、方差为 1。标准化后协方差矩阵接近单位阵两个特征权重相当。但要注意标准化后均值不再对应原始物理量解释结果时得换算回去。4.4 分量数选错导致过拟合或欠拟合现象设num_components 3时多出来的那个分量只分到几个点协方差塌缩设 2 时明明图上能看到三个簇却硬分成两个。原因GMM 本身不自动选分量数你得自己判断。解决跑 2 到 5 个分量算 BIC 或 AIC选最小的那个。BIC 公式是-2*logL k*log(n)其中k是参数个数每个分量有 2 个均值、3 个协方差参数、1 个权重共 6 个n是样本数。代码里没内置但加几行就能算。4.5 收敛判据只看均值协方差还在变现象迭代停了但协方差矩阵和权重跟上一轮比还有明显变化。原因代码里收敛检查只比较均值矩阵的范数没管协方差和权重。解决把收敛判据改成比较所有参数的变化比如norm(means - old_means) norm(weights - old_weights) sum(covariances(:) - old_covariances(:)).^2 tolerance。或者简单点把tolerance调小到1e-8多跑几步反正 100 步上限兜着。5. 进阶用法从两类到多类、从硬标签到软概率的实战技巧5.1 扩展到多类裂纹模式实际工程里裂纹模式可能不止拉伸和剪切两类比如还有混合型。把num_components改成 3 或 4代码结构不用动但初始化要更小心。我一般先用 k-means 跑一遍看簇中心分布如果某两个中心离得很近说明可能不该分那么多类。另外多类时建议用 BIC 选分量数别凭感觉。下面这段是算 BIC 的代码% 在 EM 迭代结束后计算 BIC logL 0; for i 1:num_samples pdf_sum 0; for k 1:num_components pdf_sum pdf_sum weights(k) * mvnpdf(data(i,:), means(k,:), covariances(:,:,k)); end logL logL log(pdf_sum); end num_params num_components * 6 - 1; % 每个分量 6 个参数权重和约束减 1 BIC -2 * logL num_params * log(num_samples); fprintf(BIC: %.2f\n, BIC);逻辑说明logL是对数似然逐点算混合概率密度再取对数求和。num_params每个分量有 2 个均值、3 个协方差对称矩阵、1 个权重共 6 个但权重和为 1 的约束减 1。BIC 越小越好跑不同num_components比较即可。5.2 用软概率做后续筛选硬标签只给你一个类别但软概率能告诉你这个点有多“确定”。我通常会把最大概率低于 0.7 的点单独拎出来人工检查或做二次分析。比如在疲劳裂纹监测里低置信点可能对应裂纹模式转换的过渡阶段单独看这些点的 RA-AF 轨迹能发现模式演化的前兆。代码里max_prob就是现成的加一行筛选就行ambiguous_idx max_prob 0.7; fprintf(模糊点数量: %d\n, sum(ambiguous_idx)); % 把模糊点画成空心圈叠加在散点图上 figure; gscatter(data(:,1), data(:,2), labels); hold on; scatter(data(ambiguous_idx,1), data(ambiguous_idx,2), 60, k, o, LineWidth, 1.5); hold off; title(模糊点最大概率 0.7标注);逻辑说明ambiguous_idx是逻辑索引sum统计数量。scatter用空心黑圈标出模糊点叠加在gscatter上。这样一眼就能看出模糊点集中在两个簇的交界处还是散落各处——如果集中在交界说明分界本身模糊如果散落可能是噪声或第三类模式。5.3 固定随机种子与结果复现EM 对初值敏感不固定种子每次结果都可能不同。我习惯在脚本开头加rng(42)并且用 k-means 初始化。如果换了数据k-means 的Replicates设 5 到 10确保初始中心稳定。另外把最终的means、covariances、weights保存成.mat文件下次直接加载省得重跑。下面是一段保存和加载的写法% 保存模型参数 save(gmm_model.mat, means, covariances, weights); % 下次加载 loaded load(gmm_model.mat); means loaded.means; covariances loaded.covariances; weights loaded.weights;逻辑说明save把三个变量存成二进制 mat 文件load读回来。注意变量名要一致否则得手动赋值。这样做的另一个好处是如果你调了数据预处理但想复用之前的聚类参数可以直接加载对比。从那以后我每次跑 GMM 都强制走一遍先rng固定种子再 k-means 初始化跑完看最大概率直方图低于 0.7 的点单独标出来最后存模型参数。这套流程帮我省了至少三次重跑一整天的情况。希望帮到你。本文还有配套的精品资源点击获取