简介基于 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/VAF往往是几十到几百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的后验概率矩阵。后面绘制椭圆和做在线分类时要用posteriorcluster只适合快速查看归类数量。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输出参数表如下图中椭圆和这张表配合基本够写一段方法描述了输出项变量说明混合权重gm.ComponentProportion每个类别的样本占比类中心gm.mu反标准化后对应物理量协方差矩阵gm.Sigma反标准化后椭圆形状来源硬分类标签idxcluster输出软概率postposterior输出N×K4. 让聚类结果靠谱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最小只是统计意义上的最优不代表物理上可分。K3的BIC低于K2但第三个分量可能只占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做log1pAF一般不需要。变换后的模型参数反还原尺度时均值要按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散点如果所有点都紧贴一条直线二维协方差矩阵的第二特征值接近零数值上就奇异了。解决先加RegularizationValue1e-5试试多数情况下能过还不行就把CovarianceType改成diagonal直接放弃两特征的相关性建模。根因是数据如果某类事件只有十几个点减K或者增加数据量比调参更有效。5.4 椭圆重叠面积巨大聚类出的类别没有工程意义现象两个95%置信椭圆有一半以上重叠后验概率矩阵里大部分样本的值都在0.5附近。原因数据本身可能没有双峰结构或者K取多了。强行把一坨均匀分布劈成两半GMM也能收敛但分出来的边界完全由初始化和正则化决定没有任何物理意义。解决先看标准化的特征直方图RA和AF的联合分布如果只有一个中心峰也就是单模态就别用GMM分两类。BIC曲线也会给提示K1时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的高斯混合聚类裂纹模式识别才算真正落地成了一套可以用、可解释、能交代别人复现的方案。希望帮到你。本文还有配套的精品资源点击获取