简介这份Matlab数据聚类代码包面向本科、硕士阶段开展故障诊断与模式识别研究的读者基于FCM模糊C均值聚类并结合GK、GG及PCA等方法对西储大学轴承数据进行状态识别可直接对照示例运行学习。压缩包共51个文件以40个m脚本为主体涵盖FCMclust、GKclust、GGclust、Kmeans、PCA等核心算法实现以及聚类有效性评估、归一化与反归一化、投影与降维等配套工具另有6个txt说明、2张png示意图、1个mat数据、1个pdf工具箱文档和1个wk1数据文件总大小约2.45MB。其中提供Demos、clusteringexamples、normexample、optnumber、comparing等多组示例可帮助读者对比不同聚类方法的识别效果、确定最佳聚类数并附FuzzyClusteringToolbox.pdf便于系统理解聚类流程与参数选择。目前已有127人学习下载适合需要快速上手Matlab聚类算法并进行轴承状态识别实验的教研场景也可作为课程设计或毕业设计的参考代码。1. FCM聚类与轴承状态识别为什么模糊划分比硬聚类更适合故障诊断旋转机械的轴承故障诊断本质上是一个模式识别问题从振动信号里提取特征再把特征划分到正常、内圈故障、外圈故障、滚动体故障这几类状态中去。传统K-Means这类硬聚类算法强行把每个样本归到唯一一个簇里但实际采集到的轴承振动数据在故障早期往往存在过渡状态——同一个样本可能既像正常状态又带轻微故障特征硬聚类很容易在这个边界区域发生误判。FCMFuzzy C-Means模糊C均值通过隶属度矩阵让每个样本以概率形式归属于所有簇恰好能刻画这种不确定性这也是FCM在设备状态识别里一直没被淘汰的原因。我在这类项目里的处理流程通常分四步先用西储大学CWRU公开的轴承数据库做数据源对原始振动信号分段提取时域和频域特征再把特征矩阵分别送入FCM及其变体GKGustafson-Kessel、GGGath-Geva聚类器做状态划分最后用PCA做主成分降维来对比不同特征组合下的聚类效果。这套方案的好处是计算量可控、Matlab实现成熟且不需要打标签就能看出数据的内在结构——对有监督标签的CWRU数据正好能用来验证聚类结果和真实故障类型的吻合度。适合的对象是正在做故障诊断课程设计、毕业设计或者刚接触无监督状态识别的工程师。2. 西储大学轴承数据读取与特征构造从.mat文件到聚类输入矩阵2.1 CWRU数据的目录结构与Matlab加载方式西储大学轴承数据中心提供的是采样好的振动信号常见格式是.mat文件文件名里带着载荷、故障类型和故障尺寸信息。比如105.mat表示驱动端风扇端加速度数据采样频率12kHz载荷0马力内圈故障0.007英寸.mat这类命名则直接标明故障位置和损伤程度。做FCM聚类前要先清楚每个.mat文件内部存的是DE驱动端、FE风扇端、BA基座三个通道的时间序列通常取DE通道就能保留最丰富的故障冲击特征。在Matlab里加载这类数据的标准做法是用一个循环遍历文件夹下所有.mat文件把每次的振动信号切成固定长度的样本段。切段不是为了省内存而是让每个样本点都能携带独立的统计特征——西储大学的数据连续记录了几秒到几十秒若把整段信号当成一个样本整个数据集只有几十个样本聚类毫无统计意义。合理做法是每段取1024或2048个点相邻段之间允许50%重叠这样一盘数据能产生数百个样本后续聚类才谈得上稳定。% 遍历CWRU数据文件夹读取.mat并切段 fileList dir(fullfile(CWRU_data, *.mat)); Fs 12000; % 采样频率 12kHz segmentLen 2048; % 每段样本点数 overlap 0.5; % 重叠率50% step round(segmentLen * (1 - overlap)); allFeatures []; allLabels []; for i 1:length(fileList) data load(fullfile(CWRU_data, fileList(i).name)); fieldName fieldnames(data); signal data.(fieldName{1}); % 取第一个通道通常是DE if size(signal, 1) size(signal, 2) signal signal; % 统一为列向量 end nSeg floor((length(signal) - segmentLen) / step) 1; for j 1:nSeg idx (j - 1) * step 1 : (j - 1) * step segmentLen; seg signal(idx); feats extractFeatures(seg); % 自定义特征函数 allFeatures [allFeatures; feats]; allLabels [allLabels; extractLabel(fileList(i).name)]; end end save(bearing_features.mat, allFeatures, allLabels);这段代码里extractFeatures是需要自己实现的特征提取函数extractLabel则是从文件名解析故障类型的辅助函数。注意load之后用fieldnames获取实际变量名因为不同版本的CWRU .mat文件内部变量命名不完全一致硬编码字段名容易报错。切段时使用重叠能平滑相邻段之间的特征突变但重叠率不宜超过75%否则大量样本高度相关聚类结果会显得过于干净而不真实。2.2 特征集选择时域指标与频域指标的搭配策略特征提取决定了聚类输入矩阵的质量这是整个流程里最影响结果的一步。常见的做法是每个样本段提取约12到20个统计特征包括时域的均值、标准差、均方根、峰值因子、峭度、偏度、波形因子以及频域的重心频率、均方频率、频率方差。峭度对早期故障冲击特别敏感而均方根则反映整体振动能量水平两者互补性很强。在西储大学数据上做FCM聚类特征数量不是越多越好。特征过多会导致维数灾难FCM的隶属度在欧氏距离下会被拉平区分度下降特征过少则可能漏掉故障的关键差异。我一般建议先提取全部候选特征然后用PCA做降维观察累计贡献率再决定最终保留多少维。后文会给出具体的PCA调用方式这里先把特征提取函数写完整。function feats extractFeatures(seg) % 时域特征 rms_val sqrt(mean(seg.^2)); peak_val max(abs(seg)); kurt_val kurtosis(seg); skew_val skewness(seg); crest_factor peak_val / (rms_val eps); shape_factor rms_val / (mean(abs(seg)) eps); % 频域特征 N length(seg); fft_amp abs(fft(seg)) / N; fft_amp fft_amp(1 : floor(N/2)); % 取单边谱 freqs (0 : floor(N/2)-1) * (12000 / N); fc sum(freqs .* fft_amp) / (sum(fft_amp) eps); % 重心频率 msf sum(freqs.^2 .* fft_amp) / (sum(fft_amp) eps); % 均方频率 fv sum((freqs - fc).^2 .* fft_amp) / (sum(fft_amp) eps); % 频率方差 feats [rms_val, peak_val, kurt_val, skew_val, crest_factor, ... shape_factor, fc, msf, fv]; end这里用 eps是为了防止零能量段导致除零错误。均值没有放进特征里因为振动信号的均值接近零携带的判别信息有限而峰值因子、峭度这类无量纲指标抗载荷波动能力更强。实际做聚类前特征矩阵每一列都要做Z-score标准化否则量级大的特征如均方频率会完全主导FCM的欧氏距离计算小量级特征如峭度等于没参与。3. 基于FCM的状态识别核心实现隶属度矩阵、聚类中心与参数设定3.1 FCM目标函数与迭代更新规则FCM聚类的思想是让每个样本对所有簇的隶属度之和为1通过最小化加权类内距离来确定聚类中心。目标函数写作J Σ(i1..c) Σ(k1..n) u_ik^m * d(x_k, v_i)其中c是聚类类别数n是样本数m是模糊指数fuzzinessu_ik是样本k对簇i的隶属度v_i是第i个聚类中心。迭代更新规则并不复杂先随机初始化隶属度矩阵然后交替更新聚类中心和隶属度直到目标函数变化小于阈值或达到最大迭代次数。FCM的收敛性已有大量工作验证实际使用中更重要的是m的取值和初始聚类中心的选择。Matlab自带fcm函数在Fuzzy Logic Toolbox里但很多人没有这个工具箱。更通用的做法是自己实现FCM迭代代码量不大而且能自由嵌入后续要做的GK和GG变体。下面给出基于循环的实现同时包含三种数据集的聚类演示一个二维合成数据用于白盒验证CWRU特征数据作为正式输入。3.2 手写FCM聚类器并用合成数据验证正确性function [U, V, objFunc] myFCM(X, c, m, maxIter, tol) % X: n x d 特征矩阵, c: 聚类数, m: 模糊指数(默认2) % U: n x c 隶属度矩阵, V: c x d 聚类中心 if nargin 4, maxIter 100; end if nargin 5, tol 1e-5; end n size(X, 1); % 随机初始化隶属度并归一化 U rand(n, c); U U ./ sum(U, 2); objFunc zeros(maxIter, 1); for iter 1 : maxIter % 更新聚类中心 Um U .^ m; V (Um * X) ./ sum(Um, 1); % 计算距离矩阵 n x c dist zeros(n, c); for i 1 : c diff X - V(i, :); dist(:, i) sum(diff .^ 2, 2); end % 更新隶属度 dist max(dist, eps); invDist dist .^ (-1 / (m - 1)); U invDist ./ sum(invDist, 2); % 目标函数值 objFunc(iter) sum(sum((U .^ m) .* dist)); if iter 1 abs(objFunc(iter) - objFunc(iter-1)) tol objFunc objFunc(1 : iter); break; end end end这段实现里两个关键细节值得注意距离矩阵计算没有用pdist2是为了兼容老版本Matlab而且循环计算在高维特征下反而更节省峰值内存隶属度更新公式里对距离加了eps防止出现0距离导致除零。模糊指数m取2是默认经验值m越大隶属度越模糊m接近1时退化为硬聚类。迭代终止不看固定迭代次数而是看目标函数是否进入平稳段这是FCM收敛判断的可靠信号。用经典Iris数据验证一下这个函数的输出是否与Matlab自带fcm一致注意它们初始化随机聚类中心的排列顺序可能不同但隶属度矩阵的行和应恒为1。load fisheriris; X meas; X zscore(X); % 标准化 [U, V, J] myFCM(X, 3, 2, 100, 1e-5); [~, hardCluster] max(U, [], 2); figure; gscatter(meas(:,1), meas(:,2), hardCluster, rgb, .); title(FCM对Iris数据的硬划分结果);聚类完成后通过max(U, [], 2)取得每个样本隶属度最大的簇得到硬划分结果用于可视化。这里用Iris做基准不是因为轴承数据不能用而是Iris有已知的ground truth能快速确认算法实现没有逻辑错误。在Iris上FCM的正确率通常能达到89%左右和K-Means接近但FCM的隶属度值能告诉你哪些样本处于类间模糊区——这些样本恰恰是状态识别中最值得关注的异常点。3.3 CWRU特征数据的FCM状态识别与可视化对西储大学数据聚类类别数c需要根据实际包含的故障类型来确定。如果选了正常、内圈故障、外圈故障、滚动体故障四类数据c就设为4如果还区分故障直径0.007、0.014、0.021英寸那类别数会膨胀到10个以上此时FCM仍然能运行但聚类边界的可解释性会下降建议先用PCA降到二维或三维做散点图观察自然簇的数量再定c。load bearing_features.mat; % 标准化特征 X zscore(allFeatures); % PCA降维可视化与贡献率分析 [coeff, score, latent] pca(X); explained cumsum(latent) / sum(latent) * 100; figure; plot(explained, o-); xlabel(主成分编号); ylabel(累计贡献率(%)); grid on; title(PCA累计贡献率曲线); % 用FCM聚类 c 4; [U, V, J] myFCM(X, c, 2, 200, 1e-5); [~, clusterLabel] max(U, [], 2); figure; gscatter(score(:,1), score(:,2), clusterLabel); xlabel(PC1); ylabel(PC2); title(FCM聚类结果(PCA降维显示));explained向量明确告诉你在当前特征集下前两个主成分捕获了多少信息量。若PC1PC2不到60%二维投影图里类间重叠会很严重此时的聚类可视化只能看出大致的分离趋势不能作为评判聚类好坏的唯一依据。FCM聚类完成后用clusterLabel对比文件的真实故障标签allLabels可以用混淆矩阵计算准确率。但要注意无监督聚类输出的簇编号是任意的需先做标签匹配——比如簇1可能对应真实的内圈故障而不一定是正常状态。4. GK与GG聚类对比及PCA降维场景三种聚类器的选择与调参4.1 GK与GG从球形簇到椭球形簇的扩展逻辑FCM用欧氏距离隐含假设每个簇是球形的且大小相近。轴承故障数据里不同故障类型的特征方差差异往往很大——正常状态振动信号离散程度远小于故障状态这就造成各个簇的协方差结构不同。GK聚类Gustafson-Kessel通过引入每个簇自己的协方差矩阵用马氏距离替代欧氏距离从而适配椭球形簇GG聚类Gath-Geva进一步基于高斯混合假设用后验概率来估计模糊协方差矩阵能自动感知每个簇的体积和密度差异。从实现上看GK和FCM的唯一区别是距离计算部分。FCM的dist(:,i)直接算欧氏距离而GK需要先计算每个簇的模糊协方差矩阵再算马氏距离。GG更复杂一些距离项里增加了先验概率和高斯分布的归一化因子。三者在Matlab里的运行时差异并不大数据量在几千样本量级都是秒级完成但选型上我会遵循以下原则特征维度较低且各故障类别方差接近时用FCM故障类型间能量差异明显时用GK数据分布严重不平衡、某一类样本数量很少时用GG的密度自适应。4.2 GK聚类器的Matlab实现与特征协方差约束function [U, V, F] myGK(X, c, m, maxIter, tol) % F: 1 x c cell每个元素是d x d协方差矩阵 [n, d] size(X); if nargin 4, maxIter 100; end if nargin 5, tol 1e-5; end U rand(n, c); U U ./ sum(U, 2); for iter 1 : maxIter Um U .^ m; V (Um * X) ./ sum(Um, 1); F cell(1, c); dist zeros(n, c); for i 1 : c diff X - V(i, :); % 模糊协方差矩阵 Fi (diff * (Um(:, i) .* diff)) / sum(Um(:, i)); Fi (Fi Fi) / 2; Fi Fi 1e-8 * eye(d); % 正则化防止奇异 F{i} Fi; % 马氏距离 dist(:, i) sum((diff / Fi) .* diff, 2); end dist max(dist, eps); invDist dist .^ (-1 / (m - 1)); U invDist ./ sum(invDist, 2); end end协方差矩阵正则化这一行是GK实现里最容易踩的坑。当某个簇内样本数过少或特征强相关时Fi可能变成奇异矩阵diff / Fi就会得到NaN聚类迭代直接崩溃。1e-8 * eye(d)是常见的Tikhonov正则化手段但如果你发现聚类结果对正则化系数很敏感说明特征之间共线性严重此时优先去做PCA降维而不是继续加大正则化系数。GK聚类对初始隶属度比FCM更敏感建议用FCM的聚类中心作为GK的初始中心这能明显缩短迭代时间并减少陷入局部最优的概率。4.3 PCA降维作为聚类预处理与后验证的手段PCA在这里有两个作用一是把高维特征压缩到低维消除特征间的相关性以便FCM/GK使用标准欧氏距离二是用降维后的前几个主成分做二维可视化帮助直观判断故障类型的可分性。但PCA和无监督聚类之间存在一个常见的逻辑误区——PCA是无监督降维它只求保留最大方差方向并不保证这些方向对故障分类最优。也就是说PCA得分图上看起来分离得很好的两个簇在原始特征空间里可能反而重叠。% 对比原始特征空间与PCA降维后的聚类纯度 [U_fcm_raw, ~, ~] myFCM(X, c, 2, 200, 1e-5); [~, label_raw] max(U_fcm_raw, [], 2); purity_raw computePurity(label_raw, allLabels); X_pca score(:, 1:2); % 只保留前两个主成分 [U_fcm_pca, ~, ~] myFCM(X_pca, c, 2, 200, 1e-5); [~, label_pca] max(U_fcm_pca, [], 2); purity_pca computePurity(label_pca, allLabels); fprintf(原始特征空间FCM纯度: %.2f%%\n, purity_raw * 100); fprintf(PCA(2维)后FCM纯度: %.2f%%\n, purity_pca * 100);computePurity是一个辅助函数它对聚类结果和真实标签做最佳匹配后计算正确率——因为聚类输出的簇编号是任意的所以需要遍历所有排列组合取最大匹配或者用匈牙利算法实现。这段实验告诉你一个重要结论如果你的特征提取质量足够好PCA降到2维通常会让纯度下降5到15个百分点因为信息损失直接体现在聚类边界上。反之如果PCA降维后纯度反而大幅提升说明原始特征矩阵里存在大量噪声维那些维度的特征不仅没帮助还在干扰距离计算。这个反向判断法在特征筛选时非常实用。4.4 三种聚类器在西储大学数据上的调参对比参数FCMGKGG模糊指数m2.01.5~2.51.5~2.51.5~2.0协方差正则化不需要1e-6~1e-81e-6~1e-8初始中心随机/FCM预训练FCM中心KMeansFCM中心迭代终止阈值1e-4~1e-51e-4~1e-51e-4~1e-5适用场景各簇球形、方差类似簇椭球形、方差不同样本数不均衡在西储大学数据上CWRU的四种状态正常、内圈、外圈、滚动体振动能量差异较大正常状态的均方根通常远低于故障状态的均方根这导致各簇方差结构显著不同。这种情况下GK比FCM更稳因为马氏距离天然能处理这种情况。而GG的优点是密度高的簇会被优先识别例如滚动体故障早期冲击较弱、样本特征较紧凑GG可以避免FCM把它错误地并入内圈故障簇。实际使用中可以用同样的混合数据跑三种聚类器对比各自与真实标签的匹配准确率选择最优者进入最终识别流程。提示模糊指数m不宜超过2.5。m过大会让所有样本的隶属度都趋向均匀分布聚类结果失去辨识力m接近1.2时结果逼近硬聚类但迭代稳定性变差。经验上轴承数据m2即可满足绝大多数场景。5. 状态识别效果评估与故障标签映射从聚类结果到可读结论5.1 用轮廓系数和模糊划分系数验证聚类质量聚类模型不能只看准确率因为无监督场景下没有标签可对比。轮廓系数Silhouette Coefficient可以衡量簇内紧密度和簇间分离度的相对大小在Matlab里直接用silhouette(X, clusterLabel)就能得到每个样本的轮廓值取均值得到整体指标。对FCM来说还有一个专属指标是划分系数Partition CoefficientPC mean(sum(U .^ 2, 2))PC越接近1表示隶属度越清晰越接近1/c表示模糊得几乎无法区分。% 计算轮廓系数 silhVal silhouette(X, clusterLabel); meanSilh mean(silhVal); fprintf(平均轮廓系数: %.3f\n, meanSilh); % 计算模糊划分系数 PC mean(sum(U .^ 2, 2)); fprintf(划分系数PC: %.3f\n, PC); % PC理论上在[1/c, 1]越接近1越好 % 输出每个簇的样本数与典型隶属度 for i 1 : c idx_i find(clusterLabel i); fprintf(簇%d: %d个样本, 平均最大隶属度%.3f\n, ... i, length(idx_i), mean(max(U(idx_i, :), 2))); end轮廓系数和PC两个指标要结合着看。如果轮廓系数高但PC低接近1/c说明硬划分看着清晰但模糊隶属度其实分布很均匀这通常意味着选择了过多的聚类中心c——多出来的簇在硬划分时强行瓜分了样本样本本身并不真正偏向任何簇。反过来PC高但轮廓系数低说明每个样本确实只归属一个簇但不同簇之间的间距太小远近关系不明确。两组指标互相印证能比单一准确率更全面地判断聚类质量。西储大学数据由于真实故障标签存在还可以画混淆矩阵来看聚类错误的具体模式。常见的错误模式是内圈故障和滚动体故障互相混淆因为两者都在旋转周期中产生冲击成分频域特征高度相似。如果出现这种混淆不要急着调聚类参数先检查特征提取中是否包含了足够的频带能量比特征——例如把频段按能量分成低频、中频、高频三个子带并计算各自占比往往能显著拉开内圈和滚动体的距离。5.2 聚类结果与真实故障标签的自动匹配方法当有ground truth标签但聚类簇编号是乱的无法直接比较。需要一个自动匹配流程把每个簇编号映射到真实故障类型。常见做法是构建一个c×c的计数矩阵元素(i,j)表示簇i中真实标签为j的样本数然后对计数矩阵求解最大权值匹配。简单场景下用贪心算法就够——每次找到计数矩阵里最大的元素把对应的簇和标签绑定删除该行该列并重复。c通常不超过10贪心即可获得全局最优。function matchTab mapClusterToLabel(clusterLabel, allLabels, c) matchTab zeros(c, max(allLabels)); for i 1 : c idxInCluster find(clusterLabel i); if isempty(idxInCluster), continue; end tab tabulate(allLabels(idxInCluster)); matchTab(i, tab(:,1)) tab(:,2); % 计数矩阵 end % 贪心匹配 matchPair []; tmp matchTab; for k 1 : c [maxVal, idx] max(tmp(:)); if maxVal 0, break; end [r, col] ind2sub(size(tmp), idx); matchPair [matchPair; r, col]; tmp(r, :) 0; tmp(:, col) 0; end fprintf(簇编号 - 真实标签映射:\n); for k 1 : size(matchPair, 1) fprintf(簇%d - 标签%d\n, matchPair(k,1), matchPair(k,2)); end end匹配完成后重新给每个样本赋一个对齐后的预测标签再与真实标签计算整体准确率和每一类的召回率。在西储大学数据上如果FCM调参得当整体准确率通常能到85%以上GK则在90%上下。如果准确率低于70%先回到特征工程环节检查峭度、峰值因子这些指标是否因为切段长度太短而失去统计意义——切段长度至少要有10个以上的旋转周期12kHz采样下旋转频率30Hz一个周期是400个采样点2048点约5个周期其实有点偏少建议对低速故障加大到4096点。5.3 实际部署时的故障判定逻辑与状态识别输出聚类模型训练好之后面对新采集的一段振动信号状态识别流程可以固化成一个函数[state, conf] bearingStateDetect(signal, model)。它先对信号切特征、标准化然后调用FCM的前向计算来求新样本的隶属度向量而不是重新跑整个聚类迭代——因为聚类中心已经固定只需要把新样本放进算一次距离和隶属度即可。隶属度向量中最大值对应的故障类型就是识别结果最大值本身可以作为置信度。function [state, conf] bearingStateDetect(signal, model) % model 包含聚类中心V、标准化参数mu/sigma、标签映射表 feats extractFeatures(signal); feats (feats - model.mu) ./ (model.sigma eps); dist zeros(1, size(model.V, 1)); for i 1 : length(dist) diff feats - model.V(i, :); dist(i) sum(diff .^ 2, 2); end u (dist .^ (-1/(model.m-1))) / sum(dist .^ (-1/(model.m-1))); [conf, idx] max(u); state model.labelMap(idx); end部署时值得注意的一点是新样本与训练数据分布差异过大的问题。如果新样本对所有簇的最大隶属度都很低说明它不属于任何一个已知状态可能是新的故障模式或传感器异常此时应输出“未知状态”而不是强行归属到最近的簇。判断阈值可以用训练集所有样本最小最大隶属度的统计分布来定例如取0.75分位数作为告警线。这个逻辑虽然简单但在实际设备监测里比盲目分类更有工程价值——它让聚类模型具备了故障预警能力而不只是分类能力。最后建议把这套流程封装成一份数据字典类型的model结构体包含聚类参数m、c、聚类中心V、标准化均值与标准差、标签映射表、告警阈值。这样识别模块和训练模块解耦后续无论是重新跑FCM/GK/GG对比还是换一批数据重新训练都不需要改动线上识别函数——这也是把聚类实验落地成工程模块时值得养成的好习惯。本文还有配套的精品资源点击获取
简介这份Matlab数据聚类代码包面向本科、硕士阶段开展故障诊断与模式识别研究的读者基于FCM模糊C均值聚类并结合GK、GG及PCA等方法对西储大学轴承数据进行状态识别可直接对照示例运行学习。压缩包共51个文件以40个m脚本为主体涵盖FCMclust、GKclust、GGclust、Kmeans、PCA等核心算法实现以及聚类有效性评估、归一化与反归一化、投影与降维等配套工具另有6个txt说明、2张png示意图、1个mat数据、1个pdf工具箱文档和1个wk1数据文件总大小约2.45MB。其中提供Demos、clusteringexamples、normexample、optnumber、comparing等多组示例可帮助读者对比不同聚类方法的识别效果、确定最佳聚类数并附FuzzyClusteringToolbox.pdf便于系统理解聚类流程与参数选择。目前已有127人学习下载适合需要快速上手Matlab聚类算法并进行轴承状态识别实验的教研场景也可作为课程设计或毕业设计的参考代码。1. FCM聚类与轴承状态识别为什么模糊划分比硬聚类更适合故障诊断旋转机械的轴承故障诊断本质上是一个模式识别问题从振动信号里提取特征再把特征划分到正常、内圈故障、外圈故障、滚动体故障这几类状态中去。传统K-Means这类硬聚类算法强行把每个样本归到唯一一个簇里但实际采集到的轴承振动数据在故障早期往往存在过渡状态——同一个样本可能既像正常状态又带轻微故障特征硬聚类很容易在这个边界区域发生误判。FCMFuzzy C-Means模糊C均值通过隶属度矩阵让每个样本以概率形式归属于所有簇恰好能刻画这种不确定性这也是FCM在设备状态识别里一直没被淘汰的原因。我在这类项目里的处理流程通常分四步先用西储大学CWRU公开的轴承数据库做数据源对原始振动信号分段提取时域和频域特征再把特征矩阵分别送入FCM及其变体GKGustafson-Kessel、GGGath-Geva聚类器做状态划分最后用PCA做主成分降维来对比不同特征组合下的聚类效果。这套方案的好处是计算量可控、Matlab实现成熟且不需要打标签就能看出数据的内在结构——对有监督标签的CWRU数据正好能用来验证聚类结果和真实故障类型的吻合度。适合的对象是正在做故障诊断课程设计、毕业设计或者刚接触无监督状态识别的工程师。2. 西储大学轴承数据读取与特征构造从.mat文件到聚类输入矩阵2.1 CWRU数据的目录结构与Matlab加载方式西储大学轴承数据中心提供的是采样好的振动信号常见格式是.mat文件文件名里带着载荷、故障类型和故障尺寸信息。比如105.mat表示驱动端风扇端加速度数据采样频率12kHz载荷0马力内圈故障0.007英寸.mat这类命名则直接标明故障位置和损伤程度。做FCM聚类前要先清楚每个.mat文件内部存的是DE驱动端、FE风扇端、BA基座三个通道的时间序列通常取DE通道就能保留最丰富的故障冲击特征。在Matlab里加载这类数据的标准做法是用一个循环遍历文件夹下所有.mat文件把每次的振动信号切成固定长度的样本段。切段不是为了省内存而是让每个样本点都能携带独立的统计特征——西储大学的数据连续记录了几秒到几十秒若把整段信号当成一个样本整个数据集只有几十个样本聚类毫无统计意义。合理做法是每段取1024或2048个点相邻段之间允许50%重叠这样一盘数据能产生数百个样本后续聚类才谈得上稳定。% 遍历CWRU数据文件夹读取.mat并切段 fileList dir(fullfile(CWRU_data, *.mat)); Fs 12000; % 采样频率 12kHz segmentLen 2048; % 每段样本点数 overlap 0.5; % 重叠率50% step round(segmentLen * (1 - overlap)); allFeatures []; allLabels []; for i 1:length(fileList) data load(fullfile(CWRU_data, fileList(i).name)); fieldName fieldnames(data); signal data.(fieldName{1}); % 取第一个通道通常是DE if size(signal, 1) size(signal, 2) signal signal; % 统一为列向量 end nSeg floor((length(signal) - segmentLen) / step) 1; for j 1:nSeg idx (j - 1) * step 1 : (j - 1) * step segmentLen; seg signal(idx); feats extractFeatures(seg); % 自定义特征函数 allFeatures [allFeatures; feats]; allLabels [allLabels; extractLabel(fileList(i).name)]; end end save(bearing_features.mat, allFeatures, allLabels);这段代码里extractFeatures是需要自己实现的特征提取函数extractLabel则是从文件名解析故障类型的辅助函数。注意load之后用fieldnames获取实际变量名因为不同版本的CWRU .mat文件内部变量命名不完全一致硬编码字段名容易报错。切段时使用重叠能平滑相邻段之间的特征突变但重叠率不宜超过75%否则大量样本高度相关聚类结果会显得过于干净而不真实。2.2 特征集选择时域指标与频域指标的搭配策略特征提取决定了聚类输入矩阵的质量这是整个流程里最影响结果的一步。常见的做法是每个样本段提取约12到20个统计特征包括时域的均值、标准差、均方根、峰值因子、峭度、偏度、波形因子以及频域的重心频率、均方频率、频率方差。峭度对早期故障冲击特别敏感而均方根则反映整体振动能量水平两者互补性很强。在西储大学数据上做FCM聚类特征数量不是越多越好。特征过多会导致维数灾难FCM的隶属度在欧氏距离下会被拉平区分度下降特征过少则可能漏掉故障的关键差异。我一般建议先提取全部候选特征然后用PCA做降维观察累计贡献率再决定最终保留多少维。后文会给出具体的PCA调用方式这里先把特征提取函数写完整。function feats extractFeatures(seg) % 时域特征 rms_val sqrt(mean(seg.^2)); peak_val max(abs(seg)); kurt_val kurtosis(seg); skew_val skewness(seg); crest_factor peak_val / (rms_val eps); shape_factor rms_val / (mean(abs(seg)) eps); % 频域特征 N length(seg); fft_amp abs(fft(seg)) / N; fft_amp fft_amp(1 : floor(N/2)); % 取单边谱 freqs (0 : floor(N/2)-1) * (12000 / N); fc sum(freqs .* fft_amp) / (sum(fft_amp) eps); % 重心频率 msf sum(freqs.^2 .* fft_amp) / (sum(fft_amp) eps); % 均方频率 fv sum((freqs - fc).^2 .* fft_amp) / (sum(fft_amp) eps); % 频率方差 feats [rms_val, peak_val, kurt_val, skew_val, crest_factor, ... shape_factor, fc, msf, fv]; end这里用 eps是为了防止零能量段导致除零错误。均值没有放进特征里因为振动信号的均值接近零携带的判别信息有限而峰值因子、峭度这类无量纲指标抗载荷波动能力更强。实际做聚类前特征矩阵每一列都要做Z-score标准化否则量级大的特征如均方频率会完全主导FCM的欧氏距离计算小量级特征如峭度等于没参与。3. 基于FCM的状态识别核心实现隶属度矩阵、聚类中心与参数设定3.1 FCM目标函数与迭代更新规则FCM聚类的思想是让每个样本对所有簇的隶属度之和为1通过最小化加权类内距离来确定聚类中心。目标函数写作J Σ(i1..c) Σ(k1..n) u_ik^m * d(x_k, v_i)其中c是聚类类别数n是样本数m是模糊指数fuzzinessu_ik是样本k对簇i的隶属度v_i是第i个聚类中心。迭代更新规则并不复杂先随机初始化隶属度矩阵然后交替更新聚类中心和隶属度直到目标函数变化小于阈值或达到最大迭代次数。FCM的收敛性已有大量工作验证实际使用中更重要的是m的取值和初始聚类中心的选择。Matlab自带fcm函数在Fuzzy Logic Toolbox里但很多人没有这个工具箱。更通用的做法是自己实现FCM迭代代码量不大而且能自由嵌入后续要做的GK和GG变体。下面给出基于循环的实现同时包含三种数据集的聚类演示一个二维合成数据用于白盒验证CWRU特征数据作为正式输入。3.2 手写FCM聚类器并用合成数据验证正确性function [U, V, objFunc] myFCM(X, c, m, maxIter, tol) % X: n x d 特征矩阵, c: 聚类数, m: 模糊指数(默认2) % U: n x c 隶属度矩阵, V: c x d 聚类中心 if nargin 4, maxIter 100; end if nargin 5, tol 1e-5; end n size(X, 1); % 随机初始化隶属度并归一化 U rand(n, c); U U ./ sum(U, 2); objFunc zeros(maxIter, 1); for iter 1 : maxIter % 更新聚类中心 Um U .^ m; V (Um * X) ./ sum(Um, 1); % 计算距离矩阵 n x c dist zeros(n, c); for i 1 : c diff X - V(i, :); dist(:, i) sum(diff .^ 2, 2); end % 更新隶属度 dist max(dist, eps); invDist dist .^ (-1 / (m - 1)); U invDist ./ sum(invDist, 2); % 目标函数值 objFunc(iter) sum(sum((U .^ m) .* dist)); if iter 1 abs(objFunc(iter) - objFunc(iter-1)) tol objFunc objFunc(1 : iter); break; end end end这段实现里两个关键细节值得注意距离矩阵计算没有用pdist2是为了兼容老版本Matlab而且循环计算在高维特征下反而更节省峰值内存隶属度更新公式里对距离加了eps防止出现0距离导致除零。模糊指数m取2是默认经验值m越大隶属度越模糊m接近1时退化为硬聚类。迭代终止不看固定迭代次数而是看目标函数是否进入平稳段这是FCM收敛判断的可靠信号。用经典Iris数据验证一下这个函数的输出是否与Matlab自带fcm一致注意它们初始化随机聚类中心的排列顺序可能不同但隶属度矩阵的行和应恒为1。load fisheriris; X meas; X zscore(X); % 标准化 [U, V, J] myFCM(X, 3, 2, 100, 1e-5); [~, hardCluster] max(U, [], 2); figure; gscatter(meas(:,1), meas(:,2), hardCluster, rgb, .); title(FCM对Iris数据的硬划分结果);聚类完成后通过max(U, [], 2)取得每个样本隶属度最大的簇得到硬划分结果用于可视化。这里用Iris做基准不是因为轴承数据不能用而是Iris有已知的ground truth能快速确认算法实现没有逻辑错误。在Iris上FCM的正确率通常能达到89%左右和K-Means接近但FCM的隶属度值能告诉你哪些样本处于类间模糊区——这些样本恰恰是状态识别中最值得关注的异常点。3.3 CWRU特征数据的FCM状态识别与可视化对西储大学数据聚类类别数c需要根据实际包含的故障类型来确定。如果选了正常、内圈故障、外圈故障、滚动体故障四类数据c就设为4如果还区分故障直径0.007、0.014、0.021英寸那类别数会膨胀到10个以上此时FCM仍然能运行但聚类边界的可解释性会下降建议先用PCA降到二维或三维做散点图观察自然簇的数量再定c。load bearing_features.mat; % 标准化特征 X zscore(allFeatures); % PCA降维可视化与贡献率分析 [coeff, score, latent] pca(X); explained cumsum(latent) / sum(latent) * 100; figure; plot(explained, o-); xlabel(主成分编号); ylabel(累计贡献率(%)); grid on; title(PCA累计贡献率曲线); % 用FCM聚类 c 4; [U, V, J] myFCM(X, c, 2, 200, 1e-5); [~, clusterLabel] max(U, [], 2); figure; gscatter(score(:,1), score(:,2), clusterLabel); xlabel(PC1); ylabel(PC2); title(FCM聚类结果(PCA降维显示));explained向量明确告诉你在当前特征集下前两个主成分捕获了多少信息量。若PC1PC2不到60%二维投影图里类间重叠会很严重此时的聚类可视化只能看出大致的分离趋势不能作为评判聚类好坏的唯一依据。FCM聚类完成后用clusterLabel对比文件的真实故障标签allLabels可以用混淆矩阵计算准确率。但要注意无监督聚类输出的簇编号是任意的需先做标签匹配——比如簇1可能对应真实的内圈故障而不一定是正常状态。4. GK与GG聚类对比及PCA降维场景三种聚类器的选择与调参4.1 GK与GG从球形簇到椭球形簇的扩展逻辑FCM用欧氏距离隐含假设每个簇是球形的且大小相近。轴承故障数据里不同故障类型的特征方差差异往往很大——正常状态振动信号离散程度远小于故障状态这就造成各个簇的协方差结构不同。GK聚类Gustafson-Kessel通过引入每个簇自己的协方差矩阵用马氏距离替代欧氏距离从而适配椭球形簇GG聚类Gath-Geva进一步基于高斯混合假设用后验概率来估计模糊协方差矩阵能自动感知每个簇的体积和密度差异。从实现上看GK和FCM的唯一区别是距离计算部分。FCM的dist(:,i)直接算欧氏距离而GK需要先计算每个簇的模糊协方差矩阵再算马氏距离。GG更复杂一些距离项里增加了先验概率和高斯分布的归一化因子。三者在Matlab里的运行时差异并不大数据量在几千样本量级都是秒级完成但选型上我会遵循以下原则特征维度较低且各故障类别方差接近时用FCM故障类型间能量差异明显时用GK数据分布严重不平衡、某一类样本数量很少时用GG的密度自适应。4.2 GK聚类器的Matlab实现与特征协方差约束function [U, V, F] myGK(X, c, m, maxIter, tol) % F: 1 x c cell每个元素是d x d协方差矩阵 [n, d] size(X); if nargin 4, maxIter 100; end if nargin 5, tol 1e-5; end U rand(n, c); U U ./ sum(U, 2); for iter 1 : maxIter Um U .^ m; V (Um * X) ./ sum(Um, 1); F cell(1, c); dist zeros(n, c); for i 1 : c diff X - V(i, :); % 模糊协方差矩阵 Fi (diff * (Um(:, i) .* diff)) / sum(Um(:, i)); Fi (Fi Fi) / 2; Fi Fi 1e-8 * eye(d); % 正则化防止奇异 F{i} Fi; % 马氏距离 dist(:, i) sum((diff / Fi) .* diff, 2); end dist max(dist, eps); invDist dist .^ (-1 / (m - 1)); U invDist ./ sum(invDist, 2); end end协方差矩阵正则化这一行是GK实现里最容易踩的坑。当某个簇内样本数过少或特征强相关时Fi可能变成奇异矩阵diff / Fi就会得到NaN聚类迭代直接崩溃。1e-8 * eye(d)是常见的Tikhonov正则化手段但如果你发现聚类结果对正则化系数很敏感说明特征之间共线性严重此时优先去做PCA降维而不是继续加大正则化系数。GK聚类对初始隶属度比FCM更敏感建议用FCM的聚类中心作为GK的初始中心这能明显缩短迭代时间并减少陷入局部最优的概率。4.3 PCA降维作为聚类预处理与后验证的手段PCA在这里有两个作用一是把高维特征压缩到低维消除特征间的相关性以便FCM/GK使用标准欧氏距离二是用降维后的前几个主成分做二维可视化帮助直观判断故障类型的可分性。但PCA和无监督聚类之间存在一个常见的逻辑误区——PCA是无监督降维它只求保留最大方差方向并不保证这些方向对故障分类最优。也就是说PCA得分图上看起来分离得很好的两个簇在原始特征空间里可能反而重叠。% 对比原始特征空间与PCA降维后的聚类纯度 [U_fcm_raw, ~, ~] myFCM(X, c, 2, 200, 1e-5); [~, label_raw] max(U_fcm_raw, [], 2); purity_raw computePurity(label_raw, allLabels); X_pca score(:, 1:2); % 只保留前两个主成分 [U_fcm_pca, ~, ~] myFCM(X_pca, c, 2, 200, 1e-5); [~, label_pca] max(U_fcm_pca, [], 2); purity_pca computePurity(label_pca, allLabels); fprintf(原始特征空间FCM纯度: %.2f%%\n, purity_raw * 100); fprintf(PCA(2维)后FCM纯度: %.2f%%\n, purity_pca * 100);computePurity是一个辅助函数它对聚类结果和真实标签做最佳匹配后计算正确率——因为聚类输出的簇编号是任意的所以需要遍历所有排列组合取最大匹配或者用匈牙利算法实现。这段实验告诉你一个重要结论如果你的特征提取质量足够好PCA降到2维通常会让纯度下降5到15个百分点因为信息损失直接体现在聚类边界上。反之如果PCA降维后纯度反而大幅提升说明原始特征矩阵里存在大量噪声维那些维度的特征不仅没帮助还在干扰距离计算。这个反向判断法在特征筛选时非常实用。4.4 三种聚类器在西储大学数据上的调参对比参数FCMGKGG模糊指数m2.01.5~2.51.5~2.51.5~2.0协方差正则化不需要1e-6~1e-81e-6~1e-8初始中心随机/FCM预训练FCM中心KMeansFCM中心迭代终止阈值1e-4~1e-51e-4~1e-51e-4~1e-5适用场景各簇球形、方差类似簇椭球形、方差不同样本数不均衡在西储大学数据上CWRU的四种状态正常、内圈、外圈、滚动体振动能量差异较大正常状态的均方根通常远低于故障状态的均方根这导致各簇方差结构显著不同。这种情况下GK比FCM更稳因为马氏距离天然能处理这种情况。而GG的优点是密度高的簇会被优先识别例如滚动体故障早期冲击较弱、样本特征较紧凑GG可以避免FCM把它错误地并入内圈故障簇。实际使用中可以用同样的混合数据跑三种聚类器对比各自与真实标签的匹配准确率选择最优者进入最终识别流程。提示模糊指数m不宜超过2.5。m过大会让所有样本的隶属度都趋向均匀分布聚类结果失去辨识力m接近1.2时结果逼近硬聚类但迭代稳定性变差。经验上轴承数据m2即可满足绝大多数场景。5. 状态识别效果评估与故障标签映射从聚类结果到可读结论5.1 用轮廓系数和模糊划分系数验证聚类质量聚类模型不能只看准确率因为无监督场景下没有标签可对比。轮廓系数Silhouette Coefficient可以衡量簇内紧密度和簇间分离度的相对大小在Matlab里直接用silhouette(X, clusterLabel)就能得到每个样本的轮廓值取均值得到整体指标。对FCM来说还有一个专属指标是划分系数Partition CoefficientPC mean(sum(U .^ 2, 2))PC越接近1表示隶属度越清晰越接近1/c表示模糊得几乎无法区分。% 计算轮廓系数 silhVal silhouette(X, clusterLabel); meanSilh mean(silhVal); fprintf(平均轮廓系数: %.3f\n, meanSilh); % 计算模糊划分系数 PC mean(sum(U .^ 2, 2)); fprintf(划分系数PC: %.3f\n, PC); % PC理论上在[1/c, 1]越接近1越好 % 输出每个簇的样本数与典型隶属度 for i 1 : c idx_i find(clusterLabel i); fprintf(簇%d: %d个样本, 平均最大隶属度%.3f\n, ... i, length(idx_i), mean(max(U(idx_i, :), 2))); end轮廓系数和PC两个指标要结合着看。如果轮廓系数高但PC低接近1/c说明硬划分看着清晰但模糊隶属度其实分布很均匀这通常意味着选择了过多的聚类中心c——多出来的簇在硬划分时强行瓜分了样本样本本身并不真正偏向任何簇。反过来PC高但轮廓系数低说明每个样本确实只归属一个簇但不同簇之间的间距太小远近关系不明确。两组指标互相印证能比单一准确率更全面地判断聚类质量。西储大学数据由于真实故障标签存在还可以画混淆矩阵来看聚类错误的具体模式。常见的错误模式是内圈故障和滚动体故障互相混淆因为两者都在旋转周期中产生冲击成分频域特征高度相似。如果出现这种混淆不要急着调聚类参数先检查特征提取中是否包含了足够的频带能量比特征——例如把频段按能量分成低频、中频、高频三个子带并计算各自占比往往能显著拉开内圈和滚动体的距离。5.2 聚类结果与真实故障标签的自动匹配方法当有ground truth标签但聚类簇编号是乱的无法直接比较。需要一个自动匹配流程把每个簇编号映射到真实故障类型。常见做法是构建一个c×c的计数矩阵元素(i,j)表示簇i中真实标签为j的样本数然后对计数矩阵求解最大权值匹配。简单场景下用贪心算法就够——每次找到计数矩阵里最大的元素把对应的簇和标签绑定删除该行该列并重复。c通常不超过10贪心即可获得全局最优。function matchTab mapClusterToLabel(clusterLabel, allLabels, c) matchTab zeros(c, max(allLabels)); for i 1 : c idxInCluster find(clusterLabel i); if isempty(idxInCluster), continue; end tab tabulate(allLabels(idxInCluster)); matchTab(i, tab(:,1)) tab(:,2); % 计数矩阵 end % 贪心匹配 matchPair []; tmp matchTab; for k 1 : c [maxVal, idx] max(tmp(:)); if maxVal 0, break; end [r, col] ind2sub(size(tmp), idx); matchPair [matchPair; r, col]; tmp(r, :) 0; tmp(:, col) 0; end fprintf(簇编号 - 真实标签映射:\n); for k 1 : size(matchPair, 1) fprintf(簇%d - 标签%d\n, matchPair(k,1), matchPair(k,2)); end end匹配完成后重新给每个样本赋一个对齐后的预测标签再与真实标签计算整体准确率和每一类的召回率。在西储大学数据上如果FCM调参得当整体准确率通常能到85%以上GK则在90%上下。如果准确率低于70%先回到特征工程环节检查峭度、峰值因子这些指标是否因为切段长度太短而失去统计意义——切段长度至少要有10个以上的旋转周期12kHz采样下旋转频率30Hz一个周期是400个采样点2048点约5个周期其实有点偏少建议对低速故障加大到4096点。5.3 实际部署时的故障判定逻辑与状态识别输出聚类模型训练好之后面对新采集的一段振动信号状态识别流程可以固化成一个函数[state, conf] bearingStateDetect(signal, model)。它先对信号切特征、标准化然后调用FCM的前向计算来求新样本的隶属度向量而不是重新跑整个聚类迭代——因为聚类中心已经固定只需要把新样本放进算一次距离和隶属度即可。隶属度向量中最大值对应的故障类型就是识别结果最大值本身可以作为置信度。function [state, conf] bearingStateDetect(signal, model) % model 包含聚类中心V、标准化参数mu/sigma、标签映射表 feats extractFeatures(signal); feats (feats - model.mu) ./ (model.sigma eps); dist zeros(1, size(model.V, 1)); for i 1 : length(dist) diff feats - model.V(i, :); dist(i) sum(diff .^ 2, 2); end u (dist .^ (-1/(model.m-1))) / sum(dist .^ (-1/(model.m-1))); [conf, idx] max(u); state model.labelMap(idx); end部署时值得注意的一点是新样本与训练数据分布差异过大的问题。如果新样本对所有簇的最大隶属度都很低说明它不属于任何一个已知状态可能是新的故障模式或传感器异常此时应输出“未知状态”而不是强行归属到最近的簇。判断阈值可以用训练集所有样本最小最大隶属度的统计分布来定例如取0.75分位数作为告警线。这个逻辑虽然简单但在实际设备监测里比盲目分类更有工程价值——它让聚类模型具备了故障预警能力而不只是分类能力。最后建议把这套流程封装成一份数据字典类型的model结构体包含聚类参数m、c、聚类中心V、标准化均值与标准差、标签映射表、告警阈值。这样识别模块和训练模块解耦后续无论是重新跑FCM/GK/GG对比还是换一批数据重新训练都不需要改动线上识别函数——这也是把聚类实验落地成工程模块时值得养成的好习惯。本文还有配套的精品资源点击获取
