简介：基于 RA-AF 特征的高斯混合聚类裂纹模式识别资料，面向材料力学、声发射检测与工程诊断人员，用于在无需人工划分分界线的前提下，自动区分拉伸裂纹与剪切裂纹。资源以 PDF 文档承载，内含完整 Matlab 实现代码、设计思路与逐段注释；文件共 1 个，大小 45KB，轻量易读，可配合本地 Matlab 环境直接运行或改造。已有 63 人学习浏览，适合快速搭建 GMM 聚类识别流程的读者参考。 代码围绕 EM 迭代算法展开：先构造示例 RA-AF 数据，再初始化权重、均值和协方差矩阵，通过 E 步计算各样本对每个高斯分量的责任度，M 步更新模型参数，并依据收敛阈值自动停止迭代。输出部分包含拉伸/剪切裂纹数量统计、样本类别标签、隶属概率值，以及以 RA 为横轴、AF 为纵轴的分布图。读者可借此理解 GMM 与 EM 在裂纹模式识别中的实际组织方式，也可将数据生成部分替换为自身采集的声发射参数，用于科研复现或工程验证。

1. 不画分界线的裂纹聚类：RA-AF加上高斯混合模型意味着什么

做声发射裂纹检测的工程师应该都有过这种经历：拿到一批RA-AF散点，按老传统在图上拉一条斜线，斜线上方算一类裂纹，下方算另一类。换了材料、换了传感器，斜线斜率就得重新拍脑袋，经验成分很大；数据分布一旦有重叠，那条线的说服力就很弱。基于RA-AF的高斯混合聚类走的是另一条路——把RA-AF散点当成若干个高斯子群混合生成的分布，用EM迭代求解每个点归属于哪个子群。整个过程不需要你设定分界线，决策边界是数据自己跑出来的。这篇笔记面向做结构健康监测、混凝土/岩石/金属声发射实验、以及想用无监督方法替换人工判据的从业者，从原理到Matlab代码，从参数设定到排错，按一条线讲完。

2. RA-AF特征与GMM：要聚类的到底是什么数据

2.1 RA与AF的定义和物理对应关系

RA是上升角的缩写，计算方式是上升时间除以幅值；AF是平均频率，计算方式是振铃计数除以持续时间。这两个参数都来自声发射波形的简化参数，不需要重建完整波形，只要在每个检出窗口里读出四个量：上升时间、峰值幅值、振铃数、持续时间。工程上常用代码块1这样的函数，把一段波形变成RA、AF两个标量，最后拼成N行2列的特征矩阵。

物理上，不同裂纹模式产生的声发射波形不同。拉伸裂纹通常以纵波为主，波前陡，上升时间短，但持续时间长、振铃多；剪切裂纹往往伴随滑移，波形更拖沓，上升时间偏长、频率成分偏低。已有判据正是利用这个差异，在RA-AF平面上用一条过原点的斜线划分两类模式。问题在于这条斜线的斜率在不同材料之间差得很多：在某种混凝土试块上成立的斜率，挪到金属试件上就不成立；同一种试块在加载速率变化后，散点云团的形状也会跟着变。

这一点恰好是GMM能接住的地方。把一批AE事件的RA-AF做成散点图，一般情况下能看到两到三个云团，云团的形状大致是椭圆，符合用二维高斯分布描述每一类的统计假设。GMM要做的，就是把整体分布建模成若干个高斯分布的加权相加，每个样本按照后验概率被归到各个分量上。

2.2 EM迭代求解高斯混合模型：一个E步一个M步的概率建模

这里不做完整公式推导，讲直觉。GMM假设观测数据X由K个高斯分量生成，每个分量有自己的权重πk、均值μk和协方差Σk。难点在于：我们看不到每个样本到底来自哪个分量，这是一个隐变量；而要求模型参数，又得先知道样本归属。这形成了一个先有鸡还是先有蛋的问题，EM迭代就是解开这个死结的常规套路。

初始化一组参数后，E步用当前参数算每个样本属于分量k的后验概率，文献里也叫责任；M步把这些责任当作软权重，重新估计每个分量的均值、协方差和权重。新参数会让数据的对数似然变大，然后进入下一轮E步，如此循环，直到对数似然变化量小于设定阈值。Matlab的fitgmdist内部实现的就是这个过程，使用者只需要输入特征矩阵和聚类数K。

迭代的性质在实操中要注意两点。第一，EM只能保证收敛到局部最优，不保证全局最优，要多组随机初值并行搜索，选对数似然最大的那组解。第二，协方差估计完全由数据分布形态决定，某个分量的有效样本太少，协方差就会退化到奇异矩阵，这就是后面要加正则化参数的原因。理解这两点，后续调参就不会觉得是玄学。

2.3 无需划分裂纹分界线：软分类对硬边界的替代

传统RA-AF判据的物理直觉本身没错，问题落在实现细节上，那条斜线是人为给定的模型，不是数据给的。边界一旦定死，落在边界附近又被强行归到某一侧的样本，误差就产生了。GMM给出的边界是等后验概率线，也就是说，真正落在边界的样本其实处于两个分布的重叠区域，应该接受“不确定”这个事实，而不是被一刀切。

用工程语言表达：GMM的决策面由数据的协方差矩阵结构决定，椭圆重叠程度直接反映两类裂纹信号的相似程度。决策边界可以是曲线，可以是倾斜椭圆相交出来的弧线，这比一根固定斜率的直线能描述的物理场景更丰富。这也是这个方法经常被当作卖点的原因——从聚类角度来说，它做到让数据说话，而不是让手动画线说话。

2.4 样本量与离群点：这个模型的适用边界在哪里

GMM不是万能的，它有两个明显的适用边界。一个是样本量，每个类别的有效事件数少于几十个时，协方差矩阵估计会非常不稳，椭圆画出来都是歪的；如果总体事件只有三五十个，强行聚类出来的结果没有统计意义，此时建议缩小分析窗口或者换用更保守的人工判据。另一个是离群点，电噪声、撞击噪声产生的异常事件会拉偏某个分量的均值，甚至单独形成一个人工分量，聚类结果看起来漂亮，实际是在给噪声分类。

实操里我会在特征提取后先做一遍物理滤波，比如把持续时间过短或过长的记录剔除，把幅值明显低于检出门槛的事件剔除，再用过滤后的矩阵进入GMM。这一步不属于算法部分，却是聚类结果能否反映真实裂纹模式的前提。常看到有人跳过预处理，直接聚类出四个分量，然后花大量时间解释每个分量是什么，实际上其中一两个分量只是噪声的统计影子。

3. Matlab实现路径：从声发射信号到GMM聚类的最小可运行代码

3.1 从单通道AE波形里提取RA-AF特征

第一步是把波形片段转换成特征向量。下面这个函数接收一个声发射事件的波形列向量和采样率，返回RA、AF两个值。注意上升时间的起点用相对阈值定位，而不是绝对阈值，这样换传感器、换增益时特征值不会系统性漂移。

function feat = extract_ra_af(hitWave, fs) % hitWave: 单个AE事件的波形片段（列向量） % fs: 采样率, Hz % 输出: feat = [RA, AF] ampPeak = max(abs(hitWave)); % 峰值幅值 thr = 0.1 * ampPeak; % 相对阈值 % 上升时间: 从阈值穿越点到峰值的时间 idxPeak = find(abs(hitWave) == ampPeak, 1, 'first'); startIdx = find(abs(hitWave(1:idxPeak)) <= thr, 1, 'last'); if isempty(startIdx) startIdx = 1; end riseTime = (idxPeak - startIdx) / fs; % 秒 % 持续时间: 相对阈值在波形两侧的穿越范围 endIdx = idxPeak - 1 + find(abs(hitWave(idxPeak:end)) <= thr, 1, 'first'); if isempty(endIdx) endIdx = length(hitWave); end duration = (endIdx - startIdx) / fs; % 秒 % 振铃计数: 越过阈值的极值点个数 cross = diff(sign(diff(abs(hitWave) > thr))); ringCount = sum(cross > 0); % RA 单位 ms/V, AF 单位 kHz RA = (riseTime * 1e3) / max(ampPeak, 1e-12); AF = ringCount / max(duration, 1e-12) / 1e3; feat = [RA, AF]; end

逻辑上要说明三点。第一，上升时间用的是相对阈值，10%峰值幅值作为起点比固定电压门槛稳健得多，尤其适合幅值动态范围大的声发射数据。第二，振铃计数这里用“差分的差分统计正穿越”近似，精度够聚类用，但如果要严格按计数器硬件定义，可以换成自相关或固定阈值过零统计。第三，RA和AF的量纲差得很远，RA通常是零点几到几（ms/V），AF往往是几十到几百（kHz），这为后面的标准化埋下伏笔。

批量处理时，对每个事件的波形调用这个函数，把返回的[RA, AF]逐行拼进矩阵，就得到N×2的特征矩阵。

3.2 fitgmdist完成GMM拟合：参数怎么设

Matlab自带的fitgmdist是Statistics and Machine Learning Toolbox里的函数，内部实现了EM迭代。拟合前我会做三件事：过滤异常事件、标准化、固定随机种子。固定随机种子是血泪经验，不固定的情况下同一份数据不同次运行结果可能完全不一样，论文和数据报告里尤其要把这一步写上。

% featAll: N×2 特征矩阵, 来自extract_ra_af批量拼接 % 第1步: 按物理范围过滤异常事件 valid = featAll(:,2) > 1e-3 & featAll(:,2) < 100; % AF物理合理区间 feat = featAll(valid, :); % 第2步: 标准化, 避免AF数值范围淹没RA mu = mean(feat); sd = std(feat); featStd = (feat - mu) ./ sd; % 第3步: 固定随机种子, 保证结果可复现 rng(2024); K = 2; % 先按拉伸/剪切两类处理 gm = fitgmdist(featStd, K, ... 'CovarianceType', 'full', ... 'Replicates', 10, ... 'RegularizationValue', 1e-5, ... 'Options', statset('Display', 'final', 'MaxIter', 500)); % 硬分类与软概率 idx = cluster(gm, featStd); post = posterior(gm, featStd);

参数说明按实际经验展开。CovarianceType默认full，允许每个分量的椭圆有任意朝向和大小，这对RA-AF这种两个特征相关但不完全线性的情况最合适；如果样本量少，可以改成diagonal，椭圆轴与坐标轴对齐，参数少一半，更不容易过拟合。Replicates是随机初始化搜索次数，至少设5到10，取对数似然最大的解。RegularizationValue是给协方差矩阵对角元加一个小的正偏置，防止奇异，我通常在1e-6到1e-3之间选，设太大椭圆会被人为撑圆。MaxIter 500一般够用，EM迭代到这里要么已经收敛要么陷入振荡，更早暴露问题反而是好事。

cluster输出的是硬分类标签，也就是把每个点归给后验概率最大的分量；posterior输出的是N×K的后验概率矩阵。后面绘制椭圆和做在线分类时要用posterior，cluster只适合快速查看归类数量。

3.3 输出分布图与统计数据：椭圆置信区间是验证好工具

聚类完成后最直观的输出是带置信椭圆的分布图。椭圆要画在原始RA-AF坐标里才方便和物理判据对比，所以需要把标准化后的均值和协方差反变换回去：

figure; gscatter(feat(:,1), feat(:,2), idx, 'rb', '.', 12); hold on; colors = {'r', 'b'}; for k = 1:K % 均值与协方差反标准化回原始尺度 muOrig = gm.mu(k,:) .* sd + mu; SigmaOrig = gm.Sigma(:,:,k) .* (sd' * sd); [V, D] = eig(SigmaOrig); % 95%置信椭圆: 卡方缩放系数 t = linspace(0, 2*pi, 100)'; circle = [cos(t), sin(t)] * sqrt(chi2inv(0.95, 2)); ell = circle * sqrt(D) * V'; plot(muOrig(1) + ell(:,1), muOrig(2) + ell(:,2), ... 'Color', colors{k}, 'LineWidth', 1.8); end xlabel('RA (ms/V)'); ylabel('AF (kHz)'); legend('类别1', '类别2', 'Location', 'best');

反标准化这步特别容易写错。均值反变换是gm.mu乘以sd再加mu，协方差反变换是两边各乘sd，写成gm.Sigma .* (sd' * sd)。我见过有人只反变换均值、椭圆还是标准化的，画出来位置对但形状完全是错的。chi2inv(0.95,2)是二维高斯95%置信椭圆的半径缩放系数，约等于2.45；如果想画3σ椭圆，直接换成3.0。

统计数据另存一份表，通常需要输出每类的样本占比、类中心原始坐标和95%椭圆的长短轴长度。用fprintf直接打印或者写到CSV都行，后面做报告时这些数字比图更有说服力。

% 输出核心统计: 权重、均值、协方差 for k = 1:K fprintf('分量%d: 占比=%.2f%%, ' , k, gm.ComponentProportion(k)*100); fprintf('RA均值=%.3f ms/V, AF均值=%.1f kHz\n', ... gm.mu(k,1)*sd(1)+mu(1), gm.mu(k,2)*sd(2)+mu(2)); end

输出参数表如下，图中椭圆和这张表配合，基本够写一段方法描述了：

4. 让聚类结果靠谱：K值、标准化与正则化的三个关键设定

4.1 聚类数K：用BIC曲线而不是直觉定类别数

K取2对应拉伸/剪切两类，这是常见物理预设；但真实数据经常出现第三类，比如混合模式裂纹或者噪声团。K不能只靠经验拍，要用信息准则对比。BIC在似然基础上加了参数量的惩罚，K增大带来的复杂度如果没有相应提升拟合度，BIC会上升，取最小值附近的K通常是稳的：

rng(42); Ks = 1:5; bicVals = nan(size(Ks)); aicVals = nan(size(Ks)); for i = 1:numel(Ks) gmT = fitgmdist(featStd, Ks(i), ... 'Replicates', 5, ... 'RegularizationValue', 1e-5); bicVals(i) = gmT.BIC; aicVals(i) = gmT.AIC; end figure; plot(Ks, bicVals, '-o', 'LineWidth', 1.5); hold on; plot(Ks, aicVals, '-s', 'LineWidth', 1.5); xlabel('聚类数 K'); ylabel('信息准则值'); legend('BIC', 'AIC', 'Location', 'best');

需要提醒的是，BIC最小只是统计意义上的最优，不代表物理上可分。K=3的BIC低于K=2，但第三个分量可能只占5%的样本量、椭圆巨大，这种分量说明数据里有一小撮离群点，并不是一个新裂纹模式。我一般把BIC当作筛选范围而不是最终答案，K在BIC最低点前后各取一个值，画图对比，优先看分量的样本占比和椭圆重叠程度再定。

4.2 标准化与对数变换：RA长尾的两种处理方式

RA特征几乎总是右偏的，少数高RA事件会把均值拉得很远，不做变换时聚类出来的椭圆会被拽成长条。常规做法是zscore标准化，让两个特征同权重，但标准化只解决量纲问题，不解决长尾问题。我给两组方案：如果数据RA范围在0到5之间，直接标准化就够；如果出现超过10的极端值，先做log1p变换再标准化，效果更好。

% 方案A: 直接标准化 featStd = (feat - mean(feat)) ./ std(feat); % 方案B: 对付长尾, 先log1p再标准化 featLog = [log1p(feat(:,1)), feat(:,2)]; featStd = (featLog - mean(featLog)) ./ std(featLog);

注意方案B里只对RA做log1p，AF一般不需要。变换后的模型参数反还原尺度时，均值要按exp(muLog) - 1还原，椭圆边界则直接按变换后尺度画，图里坐标轴标注成原值，但要说明坐标经过对数刻度处理。

4.3 Replicates和RegularizationValue：两个最容易救命的参数

Replicates的默认值是1，也就是只从一个随机初值开始EM，这对GMM来说几乎必坑。EM从不同初值出发可能收敛到不同局部最优，Replicates设成10的意思是跑10个不同初值、保留似然最高的那个解。代价是运行时间涨10倍，但RA-AF这种二维小样本数据集，10个replicates也就是几秒钟的事，不值得在这上面省时间。

RegularizationValue的默认值是0，样本充足时没问题；一旦某个分量的有效样本太少，或者两个特征高度相关，协方差矩阵会奇异，fitgmdist直接报错。这个参数本质是给协方差对角线加一个小偏置，保证可逆。设1e-5通常是安全的，设到1e-2以上椭圆会被明显撑圆，边界失真。矩阵奇异报错时先看数据：是不是某类事件数少于特征维数，或者过滤后样本太少，这些是退化根因，正则化只是兜底。

5. 避坑排查：GMM裂纹识别最常见的五个实际翻车点

5.1 聚类出两团，但不是拉伸/剪切，而是把噪声单独分了出来

现象：结果图里有一类占比很小、椭圆位置在RA或AF的极端角落，物理上解释不通。

原因：噪声事件在特征提取阶段混了进来。电噪声往往有极高的AF、极短的持续时间；撞击噪声往往RA巨大、AF低。它们数量不多，但因为距离中心较远，恰好被GMM当作一个独立高斯分量吸收了。

解决：聚类前先按物理范围过滤。常见做法是按持续时间和能量加门限，持续时间低于10微秒或高于几十毫秒的事件先用valid索引剔除。另外，聚完类之后检查每个分量的样本量，占比小于5%的分量多考虑是噪声团，不要强行解释成新裂纹模式。

提示：特征提取后画一张未经聚类的高密度散点图，先肉眼确认有没有明显的噪声边缘云团，比聚类后发现异常再回头省事得多。

5.2 同一份数据每次跑聚类结果都不一样，连类别编号都变

现象：连续跑三次fitgmdist，类中心略有差异，严重时两个分量的椭圆朝向都变了；有时同一组数据里“类别1”和“类别2”的编号互换。

原因：EM是随机初始化，默认单次搜索收敛到局部最优；另外聚类标签本身没有语义顺序，编号是任意的。

解决：固定随机种子，代码里在fitgmdist前加rng(2024)；同时把Replicates调到10以上，减少局部最优的影响。编号互换不影响后验概率和类中心，但做报告时要避免写“类别1是拉伸裂纹”这种绝对表述，要说“RA较高、AF较低的分量为拉伸裂纹”。

5.3 报错“协方差矩阵奇异”或“无法计算pseudo-inverse”

现象：fitgmdist运行到一半中断，提示大意为协方差矩阵不是正定矩阵，无法求逆。

原因：某个分量有效样本太少，或者标准化后的两个特征存在近似线性关系。RA-AF散点如果所有点都紧贴一条直线，二维协方差矩阵的第二特征值接近零，数值上就奇异了。

解决：先加RegularizationValue=1e-5试试，多数情况下能过；还不行就把CovarianceType改成diagonal，直接放弃两特征的相关性建模。根因是数据，如果某类事件只有十几个点，减K或者增加数据量比调参更有效。

5.4 椭圆重叠面积巨大，聚类出的类别没有工程意义

现象：两个95%置信椭圆有一半以上重叠，后验概率矩阵里大部分样本的值都在0.5附近。

原因：数据本身可能没有双峰结构，或者K取多了。强行把一坨均匀分布劈成两半，GMM也能收敛，但分出来的边界完全由初始化和正则化决定，没有任何物理意义。

解决：先看标准化的特征直方图，RA和AF的联合分布如果只有一个中心峰，也就是单模态，就别用GMM分两类。BIC曲线也会给提示，K=1时BIC最低就说明数据不支持多类划分。此时应该回去检查实验设计，是不是加载速率太低、裂纹模式没有分化，而不是继续调聚类参数。

5.5 椭圆画出来位置对不上散点图，或形状离谱

现象：图上椭圆中心好像偏了，或者椭圆横跨整个图幅。

原因：大概率是反标准化写错了。常见错误是只对均值做了反变换、协方差矩阵还是标准化的；或者把sd用成列向量但忘了转置，导致矩阵维度对不上。

解决：反标准化协方差时统一写成SigmaOrig = gm.Sigma(:,:,k) .* (sd' * sd)，sd是1×2行向量，sd'*sd是2×2外积。用plot把椭圆画出来之后，顺手验证一件事：椭圆的主轴方向应该与散点云团的主方向一致，长短轴比例应该与云团的扁度相当。如果这个对不上，数据没变的情况下先怀疑坐标变换，不要怀疑算法。

6. 进阶技巧：把离线GMM结果导出成在线实时分类器

离线聚类只是第一步，真正想把这个方法变成检测工具，要把GMM模型从训练数据里带出来。Matlab里直接用gm对象即可，它存储了每个分量的权重、均值和协方差。新来一个AE事件，提取RA-AF后做同样变换，然后调用posterior(gm, newFeatStd)，得到的就是这个新事件属于每个裂纹模式的概率。

在线使用时要给判定加一个阈值，我习惯的做法是大于0.9算明确分类，0.5到0.9算倾向某类但需要人工关注，低于0.5直接标为“模式不明确”。这样比硬分类更能反映真实状态，现场人员看到“不明确”标识反而会去复查波形，不容易漏检。

聚类结果本身也需要验证。一个可复现的做法是：把事件序列按时间对半分，两半分别跑同样的GMM流程，比较两组类中心的马氏距离。若两组类中心互相落在对方95%置信椭圆内，说明模式稳定；若差得很远，说明这个数据集的状态不平稳，离线聚类结果只对这个时间段有效。另一个更快的方法是把每个分量的平均波形拿出来看，拉伸和剪切模式平均波形在上升沿陡峭度和频率成分上应该有可辨识差异，如果没有，聚类结果仍需校准。

我自己的习惯是每次实验固定随机种子，把RA与AF的范围、过滤门限全部记录在文件头里。这让我有一次省下了大麻烦：换了一批传感器后聚类出的两个分量与旧模型完全对应不上，翻记录才发现标准物理范围变了，根源是传感器频率响应差异，不是算法坏了。

这些步骤做完，基于RA-AF的高斯混合聚类裂纹模式识别才算真正落地成了一套可以用、可解释、能交代别人复现的方案。希望帮到你。

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