简介这份资源面向材料力学、岩石断裂与工程检测方向的研究人员及研究生提供一套基于RA-AF参数与高斯混合模型GMM的裂纹模式识别Matlab代码无需人工划分拉伸与剪切裂纹的分界线直接通过EM迭代算法完成聚类。压缩包为单个PDF文件约45KB内含完整代码、注释说明与运行结果便于对照理解算法流程。内容涵盖数据生成、GMM参数初始化、E步责任度计算、M步权重与均值协方差更新、收敛判断以及拉伸裂纹与剪切裂纹的数量统计、分布图绘制和每个数据点所属类别的概率输出。读者可据此掌握从RA-AF特征到裂纹类型自动判别的实现思路并直接迁移到自己的实验数据中快速获得分类结果与可视化图表。目前已有62人学习关注适合需要开展裂纹识别建模或学习EM-GMM聚类实战的读者参考。1. 裂纹模式识别为什么可以绕开“分界线”做结构健康监测的同行大多踩过同一个坑手里攒了几百条声发射AE波形想按裂纹类型分成几类第一反应是画一条 RA-AF 分界线把点分成“拉伸型”和“剪切型”。这条线怎么画有人翻论文抄一个斜率有人拿经验值硬套结果换一批试件、换一个加载速率分类结果就全乱。更麻烦的是真实裂纹往往不是非此即彼而是两种模式混在一起硬切一刀等于人为制造误差。这篇要讲的方案核心思路就是不画那条线。它把 RA上升时间/幅值和 AF平均频率当成二维特征直接丢给高斯混合模型GMM用 EM 迭代算法去拟合数据的真实分布让算法自己“看出”有几个簇、每个簇长什么样。输出的是概率归属、分布图和统计表而不是一条拍脑袋的直线。适合谁手里有 AE 参数数据、想做无监督裂纹模式分类、又不想被分界线玄学折磨的工程师和研究生。Matlab 环境下这套代码不需要额外工具箱核心就是几段矩阵运算和迭代循环。2. RA-AF 特征与 GMM 的匹配逻辑为什么是它俩2.1 RA 和 AF 到底在描述什么RA 是 Rise Time 除以 Amplitude单位通常是 μs/V 或 ms/mV反映信号从起点爬到峰值有多“慢”。AF 是 Average Frequency一般用振铃计数除以持续时间单位 kHz反映信号在时域上的“密集程度”。这两个量放在一起物理意义很直接拉伸型裂纹比如张开型通常表现为高 RA、低 AF因为信号上升慢、持续时间长剪切型裂纹则相反低 RA、高 AF信号陡而密。但这里有个容易被忽略的点RA 和 AF 并不是严格独立的它们都受传播路径、传感器耦合、阈值设置影响。你如果直接拿原始值做聚类量纲差异会让 RA可能几百完全压过 AF可能几十距离计算就废了。所以代码里第一步一定是标准化常见做法是 z-score让两个特征均值为 0、方差为 1站在同一起跑线上。2.2 GMM 相比 K-means 和硬阈值强在哪K-means 是硬聚类一个点要么属于 A 要么属于 B边界是垂直平分线本质还是“画线”。GMM 是软聚类每个点给出属于每个高斯分量的概率比如“70% 拉伸、30% 剪切”。这个概率信息在裂纹模式识别里非常值钱因为很多信号本身就是过渡态硬分只会丢掉这种模糊性。更关键的是GMM 允许每个簇有自己的协方差矩阵也就是说簇可以是椭圆、可以倾斜而 K-means 只能假设球形簇。RA-AF 散点图里拉伸簇和剪切簇往往不是正圆而是沿某个方向拉长的椭圆GMM 能自然拟合这种形状。EM 算法就是用来迭代求解这些高斯参数的E 步算每个点属于各簇的概率M 步用这些概率加权更新均值、协方差和权重反复直到对数似然收敛。2.3 用 BIC 定簇数别再用眼睛数了GMM 需要你事先指定簇数 K。很多人凭感觉设 K2但实际数据可能有三类甚至四类。代码里我一般会跑一个 K 从 1 到 6 的循环每次算 BIC贝叶斯信息准则选 BIC 最小的那个 K。BIC 会惩罚模型复杂度避免 K 太大导致过拟合。这一步是自动的不需要人工干预也是这套方案比“画线”更客观的地方。% 假设 data 是 n×2 矩阵第一列 RA第二列 AF data [RA(:), AF(:)]; % z-score 标准化 data_norm (data - mean(data)) ./ std(data); % BIC 选簇数 maxK 6; bic zeros(maxK,1); for k 1:maxK try gm fitgmdist(data_norm, k, Replicates, 5, RegularizationValue, 1e-6); bic(k) gm.BIC; catch bic(k) inf; % 拟合失败给无穷大 end end [~, bestK] min(bic); fprintf(BIC 选出的最优簇数: %d\n, bestK);这段代码里fitgmdist是 Matlab 统计与机器学习工具箱的函数Replicates设 5 表示用 5 个不同初始值各跑一遍取对数似然最大的那次避免陷入局部最优。RegularizationValue加一个很小的正数到协方差对角线上防止某个簇只有一两个点时协方差矩阵奇异。如果你没有这个工具箱后面我会给纯手写 EM 的版本。3. 手写 EM 迭代不依赖工具箱的 GMM 实现3.1 EM 算法的两步循环拆解EM 的核心就两个公式。E 步对每个点 i 和每个簇 j算责任度 r_ij正比于簇 j 的权重 π_j 乘以高斯密度 N(x_i | μ_j, Σ_j)。M 步用 r_ij 更新 π_j、μ_j、Σ_j。听起来简单但协方差矩阵更新时容易出数值问题比如某个簇只剩一个点Σ 变成零矩阵下次迭代直接 NaN。所以实际代码里必须加正则项。function [mu, Sigma, pi, gamma] myGMM(X, K, maxIter, tol) [n, d] size(X); % 初始化随机选 K 个点作为均值 idx randperm(n, K); mu X(idx, :); Sigma repmat(eye(d), [1, 1, K]); % 每个簇初始协方差为单位阵 pi ones(1, K) / K; gamma zeros(n, K); logLik_old -inf; for iter 1:maxIter % E 步算责任度 for j 1:K diff X - mu(j, :); % 马氏距离平方 invS inv(Sigma(:,:,j) 1e-6 * eye(d)); maha sum((diff * invS) .* diff, 2); logpdf -0.5 * (d * log(2*pi) log(det(Sigma(:,:,j))) maha); gamma(:, j) log(pi(j)) logpdf; end % 对数域归一化防止下溢 maxLog max(gamma, [], 2); gamma exp(gamma - maxLog); gamma gamma ./ sum(gamma, 2); % M 步更新参数 Nk sum(gamma, 1); for j 1:K mu(j, :) (gamma(:, j) * X) / Nk(j); diff X - mu(j, :); Sigma(:,:,j) (diff * (diff .* gamma(:, j))) / Nk(j) 1e-6 * eye(d); pi(j) Nk(j) / n; end % 算对数似然判断收敛 logLik sum(maxLog log(sum(exp(gamma - maxLog), 2))); if abs(logLik - logLik_old) tol fprintf(EM 在第 %d 次迭代收敛\n, iter); break; end logLik_old logLik; end endE 步里我用对数域计算再归一化这是血泪经验直接算 exp 当数据点离均值远时容易下溢成 0整个责任度矩阵全零后面就没法更新了。M 步更新协方差时加1e-6 * eye(d)是后悔药防止奇异。tol一般设 1e-6 到 1e-4太小迭代次数多太大提前停。3.2 初始化的坑随机选点可能翻车上面用randperm随机选 K 个点当初始均值简单但风险大。如果碰巧选到两个非常近的点两个簇会抢同一批数据最后收敛到一个很差的局部最优。我一般会跑 10 次不同初始化选对数似然最大的那次。或者用 K-means 先跑一遍把 K-means 的簇中心当 GMM 初始均值稳得多。% 用 K-means 初始化 GMM 均值 [kmIdx, kmMu] kmeans(data_norm, bestK, Replicates, 5); % 把 kmMu 作为 myGMM 的初始 mu 传进去K-means 虽然也是硬聚类但它的结果拿来当 GMM 起点很合适因为 K-means 收敛快簇中心已经大致在数据密集区。这一步不增加多少计算量但能明显降低 EM 跑飞的概率。3.3 协方差类型选择全协方差还是对角GMM 的协方差矩阵可以设成全协方差每个簇有自己的完整 Σ、对角协方差只保留对角线假设特征独立、或者共享协方差所有簇共用一个 Σ。RA-AF 数据里RA 和 AF 通常负相关全协方差能捕捉这种相关性分类边界更贴合实际。但如果你的数据量很少比如少于 100 个点全协方差参数太多容易过拟合这时候对角协方差更稳。代码里我默认全协方差你可以在fitgmdist里用CovarianceType,diagonal切换。4. 从聚类结果到分布图与统计表输出怎么读4.1 散点图叠加置信椭圆聚类跑完最直观的验证就是画图。把 RA 作横轴、AF 作纵轴每个点按最大责任度染色再叠加每个簇的 1σ 和 2σ 置信椭圆。如果椭圆重叠严重说明簇数可能多了或者数据本身就不适合分那么细。figure; gscatter(data_norm(:,1), data_norm(:,2), clusterIdx); hold on; theta linspace(0, 2*pi, 100); for j 1:bestK [V, D] eig(Sigma(:,:,j)); for scale [1, 2] ellipse scale * V * sqrt(D) * [cos(theta); sin(theta)]; plot(mu(j,1) ellipse(1,:), mu(j,2) ellipse(2,:), k--, LineWidth, 1); end end xlabel(RA (标准化)); ylabel(AF (标准化)); title(sprintf(GMM 聚类结果 (K%d), bestK));gscatter按类别染色eig分解协方差矩阵得到椭圆轴方向和长度。1σ 椭圆大概覆盖 68% 的点2σ 覆盖 95%。如果两个簇的 2σ 椭圆几乎重合那说明它们统计上不可分强行分开没意义。4.2 统计表均值、协方差、归属概率图看完还得有数。每个簇的均值向量告诉你这个簇在 RA-AF 空间的位置协方差矩阵告诉你它的形状和方向。另外我一般会输出每个点的最大归属概率概率低于 0.6 的点标记为“模糊点”这些点往往是过渡态单独拎出来分析比硬塞进某一类更有价值。簇编号RA 均值AF 均值点数占比平均归属概率1-0.820.7538%0.9120.95-0.6845%0.8830.120.0517%0.63这张表里簇 3 的平均归属概率只有 0.63说明它可能不是一个独立模式而是簇 1 和簇 2 的过渡区。你可以试着把 K 降到 2 再跑一次看 BIC 是不是反而更小。4.3 把标准化结果还原到物理量前面所有计算都在标准化空间做但写报告时得用原始单位。还原很简单均值乘回标准差加均值协方差乘标准差的外积。注意协方差还原是Sigma_orig diag(std) * Sigma_norm * diag(std)不是简单乘一个数。mu_orig mu .* std(data) mean(data); Sigma_orig zeros(2,2,bestK); for j 1:bestK Sigma_orig(:,:,j) diag(std(data)) * Sigma(:,:,j) * diag(std(data)); end还原后你就能说“拉伸型裂纹的 RA 集中在 320 μs/VAF 集中在 45 kHz”这种具体结论而不是一堆标准化后的无量纲数。5. 避坑与排查这五个地方最容易翻车5.1 现象EM 迭代几次后责任度全变 NaN原因某个簇的点数趋近于零协方差矩阵行列式趋近于零log(det) 变成 -inf后面全乱。解决在 M 步更新后检查Nk(j)如果小于某个阈值比如 1e-3把这个簇的协方差重置为单位阵或者直接把这个簇去掉重新初始化。更稳妥的做法是加正则项前面代码里的1e-6 * eye(d)就是干这个的。5.2 现象BIC 选出的 K 等于数据点总数原因数据量太少或者特征维度太高BIC 公式里的惩罚项压不住似然增益。解决限制 K 的最大值不超过 sqrt(n/2)n 是点数。另外检查数据里有没有重复点或离群点一个极端离群点就能撑起一个簇。5.3 现象两个簇的均值几乎一样但协方差差很多原因GMM 允许簇有不同形状有时候两个簇中心重合但一个胖一个瘦这在实际物理意义上很难解释。解决改用共享协方差SharedCovariance,true强迫所有簇形状一致只区分位置。如果这样 BIC 变差很多说明数据确实需要不同形状的簇那就保留但要在报告里说明。5.4 现象聚类结果每次跑都不一样原因EM 对初始值敏感随机初始化导致每次收敛到不同局部最优。解决固定随机种子rng(42)或者用 K-means 初始化再跑多次取对数似然最大的。我一般跑 20 次初始化选最好的那次虽然慢一点但结果可复现。5.5 现象RA 和 AF 量纲差异导致某个特征完全主导原因没做标准化RA 数值范围可能是 AF 的十倍欧氏距离里 RA 贡献了 99%。解决一定要 z-score 标准化。如果你想让某个特征权重更大可以在标准化后再乘一个权重系数但要有物理依据别乱调。6. 进阶把 GMM 概率输出接到后续流程里这套代码跑通之后最有价值的不是那张散点图而是每个点的归属概率矩阵。你可以拿这个概率做很多事比如设一个阈值只把高置信度的点拿去训练有监督分类器低置信度的点人工标注或者把概率作为特征输入到裂纹扩展预测模型里比硬标签包含更多信息。另一个方向是处理时序数据。AE 信号是连续采集的相邻信号的 RA-AF 往往有相关性。你可以把 GMM 的概率输出做滑动平均得到一条“拉伸模式概率随时间变化”的曲线这条曲线在裂纹从拉伸主导转向剪切主导时会有明显爬升比单点分类更能反映损伤演化趋势。% 假设 gamma 是 n×K 归属概率矩阵按时间顺序排列 windowSize 20; probSmooth movmean(gamma(:,1), windowSize); % 对第一类做滑动平均 plot(1:n, probSmooth, LineWidth, 1.5); xlabel(信号序号); ylabel(拉伸模式概率 (平滑后)); title(裂纹模式随时间的演化趋势);movmean是 Matlab 自带的滑动平均函数窗口大小根据采样率调一般取 10 到 50。窗口太小曲线毛刺多太大又会把真实转变点抹平。我一般会试几个窗口选那个既能看出趋势又不失细节的。最后说个我自己的习惯每次跑完 GMM我都会把 BIC 曲线、散点图、概率演化曲线三张图放在一起看。如果 BIC 在某个 K 处明显下凹散点图上对应 K 个簇分得开概率曲线在物理预期的时间点发生转变三者一致我才认为这次聚类是可信的。三者对不上宁可回去检查数据质量也不硬解释。希望帮到你。本文还有配套的精品资源点击获取