RA-AF高斯混合聚类裂纹识别：用GMM替代人工分界线的Matlab实现

简介：基于RA-AF的高斯混合聚类裂纹模式识别Matlab代码包，面向材料力学、无损检测及声发射信号处理方向的研究人员与学生，解决传统裂纹分类依赖人工划分分界线、难以处理重叠分布数据的问题。代码采用EM迭代算法对GMM高斯混合模型进行参数估计，无需预设裂纹分界线，即可自动区分拉伸裂纹与剪切裂纹，并输出两类裂纹的分布图、统计数据以及每个样本属于各类别的概率值。资源为1个PDF文档，内含完整Matlab代码及详细注释，适合中高级用户直接阅读并迁移到自身RA-AF数据中使用。压缩包尺寸仅45KB，轻量易下载。该资源已有62人学习，说明其在相关领域具有一定参考价值。通过该代码可快速掌握GMM聚类与EM算法在裂纹模式识别中的具体实现流程，节省自行编写与调试时间。

1. RA-AF高斯混合聚类裂纹识别：它到底省掉了哪一步

做声发射的小伙伴十有八九都干过这件事：把每个撞击事件的RA（上升时间/振幅）和AF（平均频率）算出来，在平面上点一张散点图，然后凭经验画一条斜线，线的左上算拉伸裂纹，右下算剪切裂纹。图纸一换、数据一多，这条线就成了争议焦点——不同的人画出来的比例能差出十个点。基于RA-AF的高斯混合聚类裂纹模式识别，就是用EM迭代算法对RA-AF二维分布做GMM聚类，自动把样本分成拉伸和剪切两类，直接输出分布椭圆、后验概率和统计表。分界线不再靠人手画，而是由数据自己长出来。这套做法适合常做材料损伤监测、结构健康检测和声发射波形数据处理的工程师与研究者，最大的价值是去掉人眼分界的玄学成分，让裂纹识别结果能反复复现，也拿得出手写进检测报告。

2. 为什么要改用GMM而不是K-means：RA-AF数据重叠与EM软划分原理

2.1 RA和AF这两个参数在表达什么

RA是上升时间除以峰值振幅，反映的是信号到达峰值快不快；AF是振铃计数除以持续时间，反映的是信号在时间轴上的疏密程度。工程经验里，拉伸型裂纹的波形上升陡、振铃密集，通常表现为低RA高AF；剪切型裂纹的波形上升缓、振铃稀疏，通常表现为高RA低AF。于是传统的RA-AF关联分析法把这两个变量放在一张散点图上，用一条斜率为负的经验直线把两类点分开。

问题在于，真实采集到的数据并不是两个边界清晰的点团。混凝土、岩石、金属试件在加载过程中的声发射事件，往往存在大量过渡性质的信号，两个点团之间有一大片重叠区。用手画直线时，这条线稍微转一个角度，分类结果就明显变化；不同人处理同一组数据，裂纹占比差几个百分点是常事。这也是很多检测报告里RA-AF图被质疑主观的原因。

2.2 GMM用“软划分”处理重叠区

高斯混合模型假设整批数据由若干个高斯分布叠加产生。以两类为例，每个样本的概率密度写成：

p(x) = w1 * N(x | μ1, Σ1) + w2 * N(x | μ2, Σ2)

其中w是混合权重，μ是每个成分的中心，Σ是协方差矩阵。对RA-AF数据来说，两个高斯成分各自对应一类裂纹的云团，w大概等于两类事件的占比，Σ则描述这个云团的形状和朝向。拟合GMM就是估计这些参数，使得模型生成这批数据的似然最大。

EM迭代算法是拟合GMM的标准做法。E步用当前参数计算每个样本属于每个成分的后验概率，M步用这些后验概率作为权重重新估计μ和Σ，然后循环。每一次迭代对数似然单调上升，直到变化小于阈值。因为后验概率是个0到1之间的连续值，GMM在重叠区不会硬性把点塞给某一类，而是给一个“这个点更像拉伸还是更像剪切”的置信度，这对后续做波形复核非常有用。

2.3 K-means为什么会在RA-AF数据上翻车

K-means做的是硬划分：计算样本到各个聚类中心的欧氏距离，离谁近就归谁。它假设每个类在空间里的形态接近球形，但RA-AF数据两个云团通常是斜椭圆形态，RA低AF高的分布与RA高AF低的分布各自沿着一条负相关方向伸展。K-means不估计协方差结构，也没有置信度输出，边界附近的点被强制归到最近一侧，结果看上去很整齐，但掩盖了“这个点本身就模糊”的事实。

GMM则通过协方差矩阵保留了云团的形状信息。两个高斯成分的协方差矩阵不同、主轴方向不同，决策边界就不再是一条简单直线，而是两个密度曲面相等处形成的二次曲线，通常是一个椭圆边界。这个边界由数据拟合出来，不依赖操作人员的手感和经验。这也是标题里“无需划分裂纹分界线”的含义——不是没有分界线，而是分界线不再由人指定。

2.4 EM收敛后的输出到底有哪些

拟合完成后，GMM给出的产物包括每个成分的均值、协方差矩阵、混合权重，以及每个样本的后验概率矩阵。均值对应两类裂纹在RA-AF平面上的中心位置，协方差矩阵决定95%置信椭圆的长轴方向和胖瘦，权重直接反映两类裂纹在整个事件样本里的占比。Matlab里用cluster函数可以取出硬标签，用posterior函数可以取出后验概率，两套输出配合使用，先看整体比例，再看单点置信度，这是后面所有分析的基础。

3. 输入数据准备：从声发射采集参数到RA-AF特征矩阵

3.1 输入文件的列名与单位约定

GMM聚类本身只认数字矩阵，但RA-AF特征必须从声发射采集参数里算出来。我一般先让采集软件导出一张CSV表格，每一行是一个声发射撞击事件，至少包含四个字段：上升时间、峰值振幅、振铃计数、持续时间。单位不统一，后面算出来的RA和AF就是一场灾难，尤其是振幅列，必须先确认是线性值还是dB值。

RA的计算是上升时间除以峰值振幅，常规单位是μs/V；如果振幅列是mV，那么RA单位就是μs/mV，数值上会差1000倍，所以绘图时坐标轴单位必须始终写清楚。AF的计算是振铃计数除以持续时间，先换算成Hz再除以1000得到kHz。

3.2 由原始参数计算RA与AF的Matlab代码

% 读取声发射导出表，假设列名为 rose_time_us, amplitude_mV, counts, duration_us ae = readtable('ae_events.csv'); % RA = 上升时间 / 峰值振幅（线性幅值） % 振幅列单位是 mV，则 RA 单位为 us/mV，同批数据保持一致即可 RA = ae.rise_time_us ./ ae.amplitude_mV; % AF = 振铃计数 / 持续时间，换算成 kHz % counts / (duration_us * 1e-6) 得到 Hz，再除 1000 得到 kHz AF = ae.counts ./ (ae.duration_us * 1e-6) / 1e3;

这段代码的核心是先把物理量转成两个标量特征。AF那行也可以直接写成 ae.counts ./ ae.duration_us * 1e3，两种写法结果一样，前者更容易看出单位换算过程。读取表格用的readtable在较老的Matlab版本上可能不支持部分CSV格式，遇到报错时先检查列名是否包含特殊字符，或改用readmatrix加范围读取。

3.3 清洗：把0、NaN和极端值挡在聚类之前

振幅为零、持续时间为零的事件会让RA或AF变成Inf或NaN，这些坏点不处理，标准化之后均值偏移，fitgmdist还会给它们分配极端后验概率，直接污染聚类结果。我习惯在计算完RA和AF后先做一次有效性过滤，再做一个MAD离群值过滤。

% 剔除无效值与极端值 valid = RA > 0 & AF > 0 & isfinite(RA) & isfinite(AF); RA = RA(valid); AF = AF(valid); % MAD去离群：超过“中位数 ± 6倍MAD”的视为异常 med_RA = median(RA); mad_RA = median(abs(RA - med_RA)); med_AF = median(AF); mad_AF = median(abs(AF - med_AF)); keep = abs(RA - med_RA) < 6 * mad_RA & abs(AF - med_AF) < 6 * mad_AF; RA = RA(keep); AF = AF(keep); % 拼成特征矩阵：行是事件，两列分别是RA和AF X_raw = [RA(:), AF(:)]; % 标准化，消除量纲差异 mu0 = mean(X_raw); sd0 = std(X_raw); X = (X_raw - mu0) ./ sd0;

MAD比标准差更抗离群值，6倍这个系数适合几百到几千个事件的常规数据量。如果事件数特别大，超过几千个，可以收紧到5倍；事件数很少，比如不到一百个，则不建议做MAD过滤，直接看散点图手工剔除明显异常点。注意先过滤再标准化，否则离群值会把mu0和sd0带偏，后面画椭圆时的逆变换坐标也会跟着错。

4. 基于EM的GMM聚类Matlab实现：主流程代码、分布图与统计输出

4.1 用fitgmdist完成GMM聚类的主流程

Matlab统计与机器学习工具箱里的fitgmdist封装了完整的EM迭代过程，直接调用比手写EM稳定得多。下面这段是主流程，输入X是上一章标准化后的N×2矩阵，输出包括硬标签、后验概率和排序后的成分参数。

% 聚类数量：拉伸裂纹与剪切裂纹两类 k = 2; % fitgmdist：EM迭代拟合高斯混合模型 gm = fitgmdist(X, k, ... 'CovarianceType', 'full', ... % full允许椭圆斜向分布 'Replicates', 15, ... % 多组随机初始化，取似然最大者 'RegularizationValue', 1e-3); % 协方差正则，防奇异 % 提取硬标签与后验概率 idx = cluster(gm, X); post = posterior(gm, X); % 固定标签语义：RA均值较小的成分为第1类（拉伸裂纹） [~, order] = sort(gm.mu(:, 1)); label = zeros(size(idx)); for j = 1:k label(idx == order(j)) = j; end % 后验概率与成分参数按同样顺序重排，保证后续绘图一致 post = post(:, order); mu = gm.mu(order, :); Sigma = gm.Sigma(:, :, order); w = gm.ComponentProportion(order);

fitgmdist输出的成分顺序每次可能不同，这是GMM的数学性质决定的，成分本身没有固定编号。所以一定要按RA均值排序固定一下语义：RA均值低的是拉伸裂纹，RA均值高的是剪切裂纹。Replicates参数设为15，表示从15组随机初始化里挑似然最大的结果，能明显减少EM落到局部最优解的概率。RegularizationValue给1e-3，协方差矩阵不会因为某个成分样本太少而变成奇异矩阵。

4.2 分布椭圆与散点图：把聚类结果翻译成RA-AF图

聚类之后要输出一张能直接放进报告里的分布图，我习惯画两组散点加两组95%置信椭圆，坐标轴恢复成原始的RA和AF单位。椭圆参数从标准化空间逆变换回原始坐标，核心是记住逆变换公式：原始坐标 = 标准化坐标 * 标准差 + 均值。

% 原始坐标下的散点图 figure('Color', 'w', 'Position', [80 80 560 460]); colors = [0.85 0.33 0.10; 0.00 0.45 0.74]; % 类1拉伸、类2剪切 idx1 = find(label == 1); idx2 = find(label == 2); scatter(X_raw(idx1, 1), X_raw(idx1, 2), 18, colors(1, :), 'filled'); hold on; scatter(X_raw(idx2, 1), X_raw(idx2, 2), 18, colors(2, :), 'filled'); % 画95%置信椭圆 theta = linspace(0, 2 * pi, 120); for j = 1:k [V, D] = eig(Sigma(:, :, j)); r = 2.447; % 二维正态分布95%置信半径 sqrt(chi2inv(0.95, 2)) xy = V * diag(sqrt(diag(D))) * [cos(theta); sin(theta)] * r; % 逆标准化：从标准化空间映射回原始坐标 xy(1, :) = xy(1, :) * sd0(1) + mu0(1); xy(2, :) = xy(2, :) * sd0(2) + mu0(2); plot(xy(1, :), xy(2, :), '-', 'Color', colors(j, :), ... 'LineWidth', 1.6, 'HandleVisibility', 'off'); end xlabel('RA (us/mV)'); ylabel('AF (kHz)'); legend({'拉伸裂纹', '剪切裂纹'}, 'Location', 'best'); grid on; title('RA-AF GMM clustering with 95% confidence ellipse');

椭圆半径2.447是二维卡方分布95%分位数的平方根，不是一维的1.96倍标准差，画错了椭圆会偏小一圈。HandleVisibility设为off是为了让椭圆不进入图例，图例只保留两组散点的名称。逆标准化那两行是关键，很多人画完椭圆发现位置完全对不上，多半是忘了把椭圆顶点从标准化坐标变换回原始坐标。

4.3 输出统计表：权重、均值、协方差怎么读

聚类不能只给一张图，统计量要落成表格。下面代码把成分参数从标准化空间还原成原始RA-AF单位，并写出CSV文件。

% 还原到原始坐标下的成分参数 mu_orig = mu .* sd0 + mu0; % 每类中心的RA和AF值 Sigma_orig = zeros(size(Sigma)); for j = 1:k % 协方差逆变换：diag(sd0) * Sigma_std * diag(sd0) Sigma_orig(:, :, j) = Sigma(:, :, j) .* (sd0' * sd0); end % 样本计数与占比 N1 = sum(label == 1); N2 = sum(label == 2); counts = [N1; N2]; % 输出统计表 T = table((1:k)', w(:), mu_orig(:, 1), mu_orig(:, 2), counts, ... 'VariableNames', {'Cluster', 'Proportion', 'Mean_RA', 'Mean_AF', 'Count'}); writetable(T, 'gmm_crack_stats.csv'); % 命令行快速查看 for j = 1:k fprintf('Cluster %d: weight=%.3f, Mean_RA=%.2f, Mean_AF=%.2f\n', ... j, w(j), mu_orig(j, 1), mu_orig(j, 2)); end

这里最容易出错的是协方差的逆变换。标准化相当于对每个维度做了缩放，所以原坐标协方差等于 diag(sd0) * Sigma_std * diag(sd0)，也就是逐元素乘 sd0' * sd0 的外积矩阵。权重那一列可以近似理解成两类裂纹事件的占比，但严格说它是两个高斯成分的混合比例，当两类云团重叠面积较大时，它和按硬标签统计的样本占比会有几个百分点的差异，报告里建议同时写两者。

4.4 没有工具箱时的手写EM兜底代码

有些环境没有统计与机器学习工具箱，fitgmdist用不了。这时候可以手写一个简化版的双高斯EM，作为兜底和理解算法的工具。

function [label, mu, Sigma, w] = em_gmm_2class(X, maxIter) % 简易二维双高斯EM，供基础MATLAB环境兜底 % 输入X为N×2特征矩阵，输出硬标签和成分参数 [N, d] = size(X); rng(0); init = randperm(N, 2); % 随机抽两个样本作为初始均值 mu = X(init, :); Sigma = repmat(cov(X), 1, 1, 2); % 初始协方差用全局协方差 w = [0.5 0.5]; logL = -inf; for iter = 1:maxIter % E步：计算每个样本属于每个成分的后验概率 R = zeros(N, 2); for j = 1:2 R(:, j) = w(j) * mvnpdf(X, mu(j, :), Sigma(:, :, j)); end s = sum(R, 2); % 归一化因子，也是似然贡献 R = R ./ s; % M步：按后验概率加权更新参数 Nk = sum(R, 1); % 每个成分的有效样本数 w = Nk / N; for j = 1:2 mu(j, :) = (R(:, j)' * X) / Nk(j); Xc = X - mu(j, :); Sigma(:, :, j) = (Xc .* R(:, j))' * Xc / Nk(j) + 1e-4 * eye(d); end % 收敛判断：对数似然增量小于阈值 newlogL = sum(log(s)); if abs(newlogL - logL) < 1e-6 logL = newlogL; break; end logL = newlogL; end [~, label] = max(R, [], 2); end

手写EM没有Replicates机制，很容易受初始化影响，实际用的时候建议外层循环跑20次，每次随机初始化，最后保留对数似然最大的一组结果。mvnpdf也属于统计工具箱，如果连它都没有，就只能自己写二维高斯密度函数，但那种环境已经很少见了。手写版主要用于理解EM的迭代逻辑和排查fitgmdist输出的合理性，正式分析还是优先用fitgmdist。

4.5 fitgmdist关键参数速查

fitgmdist的参数直接影响聚类质量，按数据情况调整的优先级从高到低是Replicates、RegularizationValue、CovarianceType。数据量大、初始化影响小的场景可以降低Replicates到5；数据量小、经常报奇异协方差的场景，优先加RegularizationValue或切diagonal。

5. 裂纹模式识别避坑：EM聚类最常见的5个翻车现场

5.1 标签顺序漂移：跑两次结果对调

现象：同一份数据上午跑和下午跑，类1和类2的颜色、均值、后验列全部对调，统计表里拉伸和剪切的比例正好互换。原因：GMM的成分顺序在数学上没有任何语义约定，fitgmdist每次返回的成分编号只代表求解顺序，不代表物理类别。解决：按RA均值排序固定标签，代码在第4.1节已经出现过，这里强调一下，排序不仅要作用在硬标签上，还要同步重排后验概率矩阵post和成分参数mu、Sigma，否则后验置信度那一列会对不上物理类别。

5.2 协方差奇异与局部最优：EM初始化的影响不是玄学

现象：运行时报“ill-conditioned covariance created”或“the covariance of each component is not positive definite”，聚类结果里某个椭圆缩成一条线，或者两个椭圆完全叠在一起。原因：某个成分分到的样本太少，协方差矩阵用少量样本估计接近奇异；或者初始化太差，EM收敛到了对数似然的局部极值点。解决：加大Replicates到20，同时把RegularizationValue提高到1e-2，必要的时候把CovarianceType改成diagonal，放弃斜椭圆换数值稳定。要注意RegularizationValue不是越大越好，太大会把椭圆撑圆，丢失RA-AF云团的方向信息。

5.3 0、NaN、Inf混进数据：一个坏点带偏全局

现象：聚类结果整体看起来没问题，但有一类只分到极少数点，或者两个椭圆的主轴方向奇怪。翻原始数据发现有个事件的持续时间为0或振幅为0。原因：振幅0让RA变成Inf，持续时间0让AF变成Inf，标准化后mu0和sd0被拉偏，后续所有高斯密度计算都受影响，fitgmdist给这些点分配了极端后验。解决：在算完RA和AF后立即做一次总检查，打印无效行数量，趁早发现问题。过滤条件一定要同时包含正值判断和isfinite判断，两者缺一不可，filter逻辑见第3.3节代码。

5.4 强制k=2聚不出两类：BIC和物理机制对不上

现象：固定k=2时，其中一个成分的权重只有0.02，相当于第二个椭圆里只有几个点；或者两个椭圆的中心距离很近，重叠面积超过一半。原因：事件数太少、两个阶段的数据中心过于接近、或数据里混入了第三种损伤机制。解决：用BIC做模型选择快速判断，把k从1扫到4，看BIC曲线拐点在哪。

kList = 1:4; bic = zeros(size(kList)); for i = 1:numel(kList) gmT = fitgmdist(X, kList(i), 'Replicates', 5, ... 'RegularizationValue', 1e-3); bic(i) = gmT.BIC; end figure('Color', 'w'); plot(kList, bic, '-o'); xlabel('Number of components'); ylabel('BIC'); grid on;

BIC最小值对应统计上最优的分量数，但工程上不要盲从。如果k=2和k=3的BIC差距很小，优先选2，因为两种裂纹的物理语义更容易在报告里解释；如果k=3明显更好，就要回到原始波形看是不是出现了第三种损伤形式，比如混合型裂纹或噪声团。

5.5 振幅用了dB：RA直接失去物理意义

现象：RA数值小得离谱，聚类椭圆形态正常但坐标轴完全看不懂，对照文献里的结果对不上。原因：采集软件导出的振幅列是dB值，RA需要线性幅值作分母。dB是取过对数的值，直接拿它除上升时间，得到的量没有任何物理意义。解决：先确认振幅单位，若是dB则换算成线性值再算RA。

% 振幅列是dB时，先转成线性相对值（以1μV为参考） A_lin = 10 .^ (ae.amplitude_dB / 20); RA = ae.rise_time_us ./ A_lin;

如果采集软件同时导出了mV或μV的线性幅值列，直接用那列，不要自己换算，减少一次出错机会。AF的计算不涉及振幅，不受dB影响，所以这个坑只坑RA不坑AF。

6. 用后验概率做置信度分级：把模糊裂纹单独画出来

GMM比人工画线强的一个地方是它天然给每个样本输出后验概率，这个信息在常规聚类图里完全被浪费了。cluster函数返回的硬标签只保留了最大后验概率对应的类别，却丢掉了置信度本身。我现在的习惯是每次聚类完都加一张置信度分级图：后验概率超过0.75的按类别着色，低于0.75的单独画成黑色空心圆，作为“不确定区”。

% 取出每个样本的最大后验概率 conf = max(post, [], 2); uncertain = conf < 0.75; figure('Color', 'w', 'Position', [80 80 560 460]); scatter(X_raw(idx1, 1), X_raw(idx1, 2), 18, colors(1, :), 'filled'); hold on; scatter(X_raw(idx2, 1), X_raw(idx2, 2), 18, colors(2, :), 'filled'); scatter(X_raw(uncertain, 1), X_raw(uncertain, 2), 30, 'k', 'o'); legend({'拉伸裂纹', '剪切裂纹', '不确定区'}, 'Location', 'best'); xlabel('RA (us/mV)'); ylabel('AF (kHz)'); grid on; title('RA-AF clustering with confidence threshold 0.75');

阈值取0.7到0.8都合理，具体看验收要求。不确定区占比高说明两个高斯成分重叠严重，RA-AF参数本身区分度不够，这时候硬分类不如回看波形。我曾在某组应力腐蚀数据上看到不确定区占比19%，回去逐条看波形，发现这些事件的波形确实是拉伸和剪切特征的混合形态，硬归到任何一侧都会误导裂纹比例结论。于是我做了一个习惯：把不确定区比例当作每次试验的“分离质量”指标记录在结果表里，比例超过15%就提醒自己和团队，这批数据的RA-AF区分能力有限，需要结合波形特征或更多参数交叉验证。

这个习惯帮我挡过不少无效结论，也让GMM聚类不只是“画一个漂亮的椭圆图”，而是真正变成可以量化、可复现的裂纹模式识别流程。希望帮到你。

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