简介这份资源是一套基于RA-AF特征与高斯混合模型GMM的裂纹模式识别Matlab实现面向材料力学、无损检测及信号处理方向的研究人员与工程师用于在无需划分裂纹分界线的前提下通过EM迭代算法完成GMM聚类区分拉伸裂纹与剪切裂纹并输出分布图、统计数据及每个样本的类别概率。压缩包为单个PDF文档体积仅45KB内含完整的Matlab代码、运行逻辑说明与注释便于对照学习与二次开发。目前已有62人浏览学习适合具备一定Matlab基础和统计学知识的读者使用。通过该代码可实现完整的GMM聚类流程数据生成与预处理、参数初始化、E步与M步迭代更新、收敛判断、分类标签输出以及可视化能够帮助理解EM算法在材料裂纹识别中的实际应用也可为其他聚类任务提供可复用的代码框架。1. 裂纹分界线画不准不如让GMM自己聚类RA-AF 思路到底在解决什么做裂纹模式识别的人大多经历过同一个尴尬先手动画出每一条裂纹的分界线再按线内灰度、宽度、方向去归类。光照一变、噪声一多、裂纹连成网状时分界线就画得越来越玄学换个人画结果能差一大截。RA-AF 加高斯混合聚类GMM的路线是想把这个前置步骤直接拿掉不划任何分界线把图像切成子区域每个子区域提一组 RA-AF 融合特征空间统计属性加频域属性交给 EM 迭代的 GMM 模型自动聚类最后输出分布图和统计表。适合正在做混凝土、沥青或金属疲劳表面裂纹分类的同学也适合想把手动阈值分割换成无监督聚类的工程人员。2. RA-AF 特征提取不画分界线用特征向量描述每个子区域2.1 先弄懂 RA-AF 是哪两类特征的融合RA-AF 在不同资料里展开并不完全统一但落到操作层面指向很一致RA 指区域属性Region Attribute负责描述“这个小窗口在空间域长什么样”AF 指频域属性或角频率属性负责描述“窗口里的纹理密度、方向性有多强”。把两类特征首尾拼接就得到一个融合特征向量。为什么这样做就能不划裂纹分界线因为聚类的对象从“像素连通域”变成了“特征空间里的点”。传统做法要先判断哪些像素属于裂纹再做连通域分析而 RA-AF 路线把每个滑窗当成一个独立样本不管窗口里有没有裂纹、裂纹在什么位置都提取同样长度的特征向量。GMM 聚类时在特征空间里自动把“密集裂纹区”“背景区”“微裂缝过渡区”的样本分开像素级的边界自然就不用画了。这里的关键在于特征设计决定聚类上限GMM 只是把特征空间里真实存在的聚集结构找出来。举个具体的数量级一块 512×512 的表面灰度图取 win64、step16会得到 29×29841 个窗口每个窗口输出一组特征最终喂给 GMM 的就是 841×d 的特征矩阵。样本量、特征维度和聚类类别数之间的关系后面选参数时会反复提到。2.2 在 Matlab 里提取 RA 区域属性的最小代码我一般先把灰度图转成 double再固定步长滑窗取块。对每个窗口RA 部分取四个常用统计量局部均值、局部方差、梯度幅值均值、局部粗糙度相邻像素差绝对值之和。最小实现如下% img_gray: 单通道灰度图, double类型, 数值范围[0,1] % win: 窗口边长(像素), step: 滑动步长(像素) function featRA extractRA(img_gray, win, step) [H, W] size(img_gray); rows 1:step:H-win1; cols 1:step:W-win1; nWin length(rows) * length(cols); featRA zeros(nWin, 4); % 每行对应一个窗口的RA特征 idx 0; for r rows for c cols patch img_gray(r:rwin-1, c:cwin-1); gx diff(patch, 1, 2); % 水平差分: win x (win-1) gy diff(patch, 1, 1); % 竖直差分: (win-1) x win gx gx(1:end-1, :); % 裁成 (win-1) x (win-1) gy gy(:, 1:end-1); % 裁成 (win-1) x (win-1) gradMag sqrt(gx.^2 gy.^2); idx idx 1; featRA(idx, 1) mean(patch(:)); % 区域均值 featRA(idx, 2) var(patch(:)); % 区域方差 featRA(idx, 3) mean(gradMag(:)); % 梯度均值 featRA(idx, 4) sum(abs(diff(patch(:)))); % 粗糙度代理量 end end end逻辑说明每个窗口输出一个 1×4 向量把所有窗口的行堆起来就是 N×4 的 RA 特征矩阵。均值与方差反映窗口明暗和对比度梯度均值反映边缘密集程度最后的粗糙度代理量对裂纹那种随机高频纹理比较敏感。参数说明win 要结合图像分辨率和裂纹宽度取裂纹在图上占 515 个像素时win3264 比较稳妥step 决定样本量step 越小样本越密但相邻窗口相关性越强我一般取 win 的 1/4。注意 diff 会让梯度图每维少一个像素代码里做了统一裁边图像边长不能被窗口整除时循环条件会自动丢掉最后一截边缘影响不大更稳的做法是先 padarray 补边。2.3 AF 频域特征与拼接顺序能不能直接喂给 GMMRA 只能描述“空间上长什么样”两个模式在空间统计量上可能完全相同但一个呈横向条状、一个呈点状噪点这时需要频域属性出场。AF 部分我对窗口做二维 FFT提取低频能量占比、主峰归一化幅值和主峰方向角三个量function featAF extractAF(img_gray, win, step) [H, W] size(img_gray); rows 1:step:H-win1; cols 1:step:W-win1; nWin length(rows) * length(cols); featAF zeros(nWin, 3); % 低频能量占比, 主峰归一化幅值, 方向角 idx 0; for r rows for c cols patch img_gray(r:rwin-1, c:cwin-1); F fftshift(fft2(patch - mean(patch(:)))); % 去直流再FFT amp abs(F); [m, n] size(amp); cy floor(m/2) 1; cx floor(n/2) 1; % 零频位置 [Y, X] meshgrid(1:m, 1:n); dist sqrt((X - cx).^2 (Y - cy).^2); lowMask dist max(m, n) * 0.15; % 低频半径阈值 idx idx 1; featAF(idx, 1) sum(amp(lowMask)) / (sum(amp(:)) eps); [~, imax] max(amp(:)); featAF(idx, 2) amp(imax) / (sum(amp(:)) eps); featAF(idx, 3) atan2(Y(imax) - cy, X(imax) - cx); end end end逻辑说明先减均值再 FFT去掉直流分量对低频占比的干扰fftshift 把零频移到中心方便用半径阈值圈出低频区。低频能量占比高说明窗口纹理平滑占比低说明高频噪声或微裂纹丰富。主峰方向角对条状裂纹很有区分度单条定向裂纹的主峰方向很集中网状裂纹的主峰方向会比较发散。参数说明低频半径取 0.15 倍最大边长是经验起点图像噪声强就降到 0.1主峰幅值除以全谱能量消除窗口面积差异带来的尺度影响。拼接时把 featRA 和 featAF 横向合并成 N×7 特征矩阵喂给 GMM 前必须做 z-score 标准化 —— 这一步不做的话方差大的特征会主导距离计算协方差估计也容易病态。3. 用 EM 迭代跑 GMM 聚类核心代码与三个必调参数3.1 为什么选 GMM 而不是 K-Means 或 DBSCAN高斯混合聚类GMM的核心假设是数据由若干个高斯分布混合生成每个分量代表一种模式。裂纹模式之间没有硬边界同一张图里密集裂纹区到背景区是渐变的边界处的滑窗很可能同时包含裂纹和背景。K-Means 做硬分配过渡区样本会被生硬塞进某一类而且它隐含假设簇是各向同性的球状结构特征尺度不同时表现很差。DBSCAN 对密度参数 eps 和 minPts 太敏感7 维特征空间里密度分布很难估常常出现“调半天参数换一张图就失效”的情况。GMM 的优势是三件事一是软划分给每个样本“属于第 k 类的概率”过渡区的不确定性被保留下来二是每个高斯分量自带协方差矩阵能描述裂纹特征空间里的椭圆分布三是 EM 迭代过程输出对数似然可以直接判断收敛状态也为后面 BIC 选 k 提供依据。常用资料里把 GMM 和 EM 绑在一起说本质上是同一件事先猜一组参数E 步算后验概率M 步更新均值、协方差和混合系数循环到对数似然增量小于阈值为止。用 fitgmdist 时看不到这些细节但要关注的三个参数正好对应上手时最容易出错的地方k 怎么定、协方差结构怎么选、正则化加不加。3.2 fitgmdist 跑通聚类的完整代码% featAll: N x d 的特征矩阵(已用zscore标准化) % k: 想聚成的模式类别数 rng(42); % 固定随机种子, 保证结果可复现 options statset(MaxIter, 1000, Display, final); mdl fitgmdist(featAll, k, ... Start, plus, ... % 用kmeans结果做初始化 CovarianceType, full, ... % 每个分量独立满协方差 RegularizationValue, 1e-5, ... % 防止协方差奇异 Options, options); label cluster(mdl, featAll); % 硬划分: 取后验概率最大的类逻辑说明fitgmdist 用 EM 迭代估计每个分量的均值向量、协方差矩阵和混合系数。cluster 返回的是对后验概率做 argmax 的硬标签适合做后续统计如果需要软概率做不确定度分析改用 mdl.posterior(featAll)。参数说明Start 取 plus 先用 k-means 算一遍初始中心比随机初始化稳定得多CovarianceType 建议先试 full样本不足或特征强相关导致奇异时降级到 diagonalRegularizationValue 手动给 1e-5 到 1e-3这是避免 Ill-conditioned covariance 警告最直接的手段。3.3 初始化方式与协方差结构对结果的影响先说初始化。EM 对初始值敏感初始中心选得差会收敛到局部最优。k-means 按距离平方加权采样分散初始中心是 7 维特征矩阵上的默认首选。数据量大时还可以设 Replicates, 5让 MATLAB 从 5 组不同起点跑 5 遍保留对数似然最大的模型相当于买一份默认的“后悔药”代价是耗时涨约 5 倍。注意 Replicates 本身带随机性所以 rng 必须在 fitgmdist 之前设置否则“跑多遍取最优”又变成不可复现的随机过程。再说协方差结构。三种可选结构的取舍可以这样理解协方差结构特点适用场景full每个分量独立形状与方向各类形状差异大、样本充足diagonal各特征独立椭圆轴平行坐标轴样本少、特征相关性弱tied所有分量共享同一个协方差特征已标准化、追求训练速度full 适合“背景类接近球形、裂纹条带类是很扁的椭圆”这种异质分布但某个类的有效样本数比特征维度还少时full 协方差必然奇异。特征维度是 7每个类少于 50 个样本我就换成 diagonal。tied 训练最快但对各类形状差异明显的纹理数据欠拟合。这三个选项没有绝对正确要靠验证集上的分离度和物理意义判断。4. 输出分布图与统计数据从聚类标签到一张能放进报告的图4.1 主特征空间的散点分布图怎么画GMM 在高维特征空间聚类但人眼只能看懂二维投影。常见做法是从 7 个特征里挑两个类间区分度最高的画散点图。我在每个项目里都先算各特征在类别标签下的组间方差与组内方差之比取比值最大的两列作为投影轴避免拍脑袋选特征导致“图上挤成一团报告里说不清”。画散点用 gscatter% 例: 用第1列(区域均值)和第4列(粗糙度代理量)做投影 colors lines(k); figure; gscatter(featAll(:,1), featAll(:,4), label, colors, ., 18); xlabel(RA: 区域均值 (标准化)); ylabel(RA: 粗糙度代理量 (标准化)); title(GMM 聚类结果的二维投影分布); grid on;逻辑说明gscatter 的第三个参数是标签向量同标签的点自动同色能直观看出类别在投影平面上的分团情况。参数说明点形统一用 . 时字号设 18 不容易糊成一团如果 k 大于 6配色会自动循环建议改用灰度或直接缩小投影特征数。二维投影只是可视化辅助不代表聚类在 7 维空间里的完整结构所以下面的椭圆叠加才是展示协方差的关键。4.2 高斯椭圆等高线与数据点叠加散点能看出聚集趋势但要体现每个高斯分量在特征空间的覆盖范围需要把均值和协方差画成椭圆等高线。做法是取模型里的 mu 和 Sigma对投影维度做 eig 分解得到主轴方向和半径figure; hold on; colors lines(k); for ki 1:k idxKi (label ki); plot(featAll(idxKi,1), featAll(idxKi,4), ., Color, colors(ki,:)); mu2 mdl.mu(ki, [1,4]); % 投影维度的均值 Sigma2 mdl.Sigma([1,4],[1,4],ki); % 投影维度的协方差 [V, D] eig(Sigma2); theta linspace(0, 2*pi, 100); r 2.448; % 95%置信椭圆半径系数 ellipse V * diag(sqrt(diag(D)) * r) * [cos(theta); sin(theta)]; plot(mu2(1)ellipse(1,:), mu2(2)ellipse(2,:), ... Color, colors(ki,:), LineWidth, 1.5); end xlabel(RA: 区域均值 (标准化)); ylabel(RA: 粗糙度代理量 (标准化)); title(GMM 分量椭圆与数据点叠加); hold off;逻辑说明eig 的特征向量 V 是主轴方向特征值开根号是主轴上的标准差乘以 2.448 对应包含约 95% 概率质量的椭圆。中心是高斯分量的均值形状反映协方差大小与朝向。参数说明这里手写椭圆而不是画概率密度曲线是为了只显示聚类结构、不显示无关的背景等高线如果模型用了 diagonalSigma2 自动退化为对角阵椭圆轴向与坐标轴平行代码不用改。投影轴区分度不足时椭圆会长条重叠那说明该换投影特征而不是调椭圆参数。4.3 导出统计表混合比例、均值向量与协方差矩阵报告和论文需要数字支撑。fitgmdist 的模型对象自带 ComponentProportion、mu 和 Sigma再把每个类别内部的特征极值与均值算出来导出一张表k mdl.NumComponents; summaryT array2table([(1:k), mdl.ComponentProportion(:)], ... VariableNames, {类别编号, 混合比例}); for fi 1:size(featAll,2) colMin zeros(k,1); colMax zeros(k,1); colMean zeros(k,1); for ki 1:k sub featAll(labelki, fi); colMin(ki) min(sub); colMax(ki) max(sub); colMean(ki) mean(sub); end summaryT.(sprintf(特征%d最小, fi)) colMin; summaryT.(sprintf(特征%d最大, fi)) colMax; summaryT.(sprintf(特征%d平均, fi)) colMean; end writetable(summaryT, gmm_cluster_stats.xlsx);逻辑说明每一行对应一个聚类类别列包含混合比例以及每个特征在类内的最小、最大、均值直接是报告里表格素材。参数说明协方差矩阵是 d×d×k 的数组不适合压进 xlsx 表格我会单独存成 .mat 文件保留如果特征做过 z-score均值向量的单位是标准差倍数报告里要注明“标准化后统计”。混合比例对裂纹模式识别的工程意义很大比如“某类模式占比 12%”直接对应这种裂纹形态在整个样本中的权重。5. 避坑指南GMM 裂纹模式识别里最容易翻车的五个点5.1 每次聚类结果都不一样怎么办现象同一张图、同一个 k今天跑和明天跑结果不同有时类别合并有时类别裂开甚至标签颜色整体对不上。原因EM 初始化带随机性cluster 输出的是后验 argmax初始化不同收敛点就不同另外类别编号本身是任意的类别 1 和类别 2 互换不算 bug。解决脚本开头固定 rng(42)并把 Replicates 设成 35对比两次实验时用混淆矩阵对标签做重排固定一次实验的标签为基准下一次的每个类别找出与基准重叠样本最多的类别编号做映射重排后一致才算可复现。我经历过“跑三遍三个结果”的混乱后来固定随机种子成了肌肉记忆这类问题基本绝迹。5.2 迭代不收敛警告信息怎么读现象命令行出现 Maximum iterations reached 或 Ill-conditioned covariance 的红色警告但模型对象照常返回后续画图还能出图。原因协方差矩阵奇异常见于某个类样本太少、特征维度高于样本数或两个特征强线性相关导致行列式接近零。解决按顺序排查 —— 先确认特征矩阵做过 z-score然后给 RegularizationValue 加 1e-5 到 1e-3再把 CovarianceType 降到 diagonal最后看每个类别的样本占比低于 5% 就减少 k 或增加样本。一个容易忽略的点看到警告不等于模型不可用但报告里引用统计量前必须用 BIC 对比没警告的模型否则数字可能没有意义。5.3 分布图重叠严重特征没挑对还是分量数不对现象散点图上三类点混在一起椭圆大面积交叠类别边界看不出物理对应。原因一是所选两个投影轴区分度不够裂纹和背景在投影上的分布本来就接近二是 k 选大了GMM 为了拟合密度会生成多余分量多余分量天然落在同一片区域。解决先用 k2 跑一遍确认能把“裂纹/背景”分开能分开再逐步加类别同时把投影轴换成 AF 的低频能量占比和 RA 的梯度均值频域特征对纹理差异更敏感。具体做法是先算每个特征在类别间的组间/组内方差比把比值最小的两三个特征记下来聚类后画图验证。这两个动作能区分“特征问题”和“分量数问题”比盲目调参高效得多。5.4 窗口尺寸改变后类别全变了现象win32 聚出三类win64 聚出四类且每个类别的特征均值对不上。原因RA 和 AF 都是尺度相关特征窗口小捕捉的是微裂纹纹理窗口大捕捉的是整个纹理场的统计平均特征含义已经变了。这不是 bug是特征本身的尺度敏感性。解决窗口尺寸依据裂纹宽度和图像分辨率定让每个窗口至少覆盖 35 条裂纹纹理周期整个数据集只用同一个 win 和 step特征提取脚本跑完一批数据后不要再改窗口尺寸。我习惯在项目记录里注明“本批结果基于 win64/step16”这句备注能省掉后面无数返工。5.5 用 silhouette 评估 GMM 的坑现象用 silhouette 系数选 k得出 k5但第 5 类只是一条渐变过渡带物理上称不上独立裂纹模式。原因silhouette 假设类别是紧凑凸集GMM 允许软重叠过渡区样本天然被 silhouette 惩罚并且 silhouette 偏向球状紧致的划分而裂纹模式往往就是扁椭圆且部分重叠的。解决选 k 主要看 BIC 曲线拐点加椭圆可视化silhouette 只做参考更实用的验证是检查每个类别的特征均值是否对应明确的裂纹形态比如“高梯度均值 高低频能量占比”对应密集微裂纹“低梯度 高低频占比”对应平滑背景。对不上物理含义的模型即使指标好看也要重跑另一个经验是当 BIC 曲线没有拐点时先别急着换聚类算法回到原始图像上看看是不是裂纹形态本身连续可分。6. 验证进阶BIC 定分量数与一个稳定的复现习惯6.1 用 BIC 曲线确定高斯分量数k 是 GMM 里唯一需要人为指定的核心参数。fitgmdist 返回的模型自带 BIC。我跑项目时会把 k1 到 k8 的 BIC 画成曲线拐点就是复杂度与拟合度的平衡位置ks 1:8; bic zeros(size(ks)); for i 1:length(ks) mdl_tmp fitgmdist(featAll, ks(i), Start,plus, ... RegularizationValue,1e-5, Options, statset(MaxIter,1000)); bic(i) mdl_tmp.BIC; end figure; plot(ks, bic, o-, LineWidth, 1.2); xlabel(高斯分量数 k); ylabel(BIC); grid on;逻辑说明BIC 在似然之外加参数数量惩罚k 增大时曲线先降后升拐点的 k 就是推荐值。参数说明BIC 的绝对数值没有意义只看曲线形状如果单调下降没有拐点说明数据本身不支持明显聚类结构再加分量只是拟合噪声要回到特征设计找原因。6.2 复现性习惯设置随机种子与固定初始化GMM 的随机性集中在初始化。固定随机源后整个流程就是确定性的我现在每个脚本开头都固定 rng(42)聚类完立刻把模型和标签存成 .mat后续画图、写报告都从文件读取不再重新跑一遍换机器、换 MATLAB 版本也不怕结果漂移。rng(42); % 固定全局随机种子 mdl fitgmdist(featAll, k, Start, plus, ... Options, statset(MaxIter, 1000, TolFun, 1e-6)); save(gmm_model.mat, mdl, label, featAll);逻辑说明三行脚本解决的是一类很难察觉的问题 —— 随机初始化导致标签整体置换或收敛点漂移。参数说明TolFun 控制 EM 迭代停止阈值默认值通常够用设小一点能多迭代几轮压榨一点对数似然。6.3 把同一套流程搬到声发射信号上这套“不划边界 多域特征 GMM聚类”不只用于图像。做声发射裂纹监测时我把每个时间窗当一个样本RA 部分换成 RMS、峭度、上升时间AF 部分换成功率谱重心频率和带宽特征矩阵还是 N×7 左右fitgmdist 的代码一行不用改。差别是滑窗从空间变成时间聚类后我会按时间轴统计每类占比变化用来观察裂纹从萌生到扩展的阶段转移。我在最先接触这套方案时也翻过车第五节五个坑几乎每个都踩了一遍。现在的固定习惯是先定窗口尺寸和随机种子再跑 BIC 选 k最后看椭圆图确认物理意义 —— 顺序不乱GMM 结果基本都在预期内。希望帮到你。本文还有配套的精品资源点击获取