简介这是一套K均值聚类分析的可执行源码与演示数据包面向机器学习初学者和故障诊断研究人员。针对故障类型未知的问题利用K均值自动划分样本并以CH指数筛选最佳聚类数避免人工指定带来的主观性。压缩包内共2个文件一个M脚本负责聚类实现、指数计算与可视化结果输出一个电子表格存放异常加热片段等原始数据包体约57KB数据经简单预处理后可直接作为算法输入便于快速复现。已有4508人学习/浏览适合通过实例掌握K均值聚类、聚类数评估以及散点图和统计图表展示流程无需额外配置环境可用于课程设计、毕业设计或入门实践也适合在故障样本类别数量不详时快速建立聚类分析流程。代码结构清晰可以直接运行能输出不同类别的分布图可作为故障分类和模式识别项目的参考实现。1. KMeans聚类故障类型分析里那些没人手写的分类规则在设备故障诊断场景里最头疼的事情不是数据不够而是你根本不知道眼前的故障到底该分成几类。轴承磨损、电气老化、冷却异常、负载突变这些故障模式在波形上互相交叠靠人眼去分类完全不可行。KMeans聚类作为无监督学习中应用最广的算法不需要任何标签直接对样本的分布结构下手把相似的数据点自动归堆。这份资源包里的Copy_of_Cluster_Flow.m脚本配合ML_Cluster.zip中的数据文件就是一个完整的故障数据聚类分析流程——从读入Excel文件、标准化处理、自动评估最佳聚类数到最终把分类结果以不同颜色的散点图呈现出来。适合那些拿到了设备监测数据但不知道从哪儿开始做分类的工程师也适合正在学聚类算法、需要一个能跑通的完整例子的学生。2. KMeans核心原理与Calinski-Harabasz指数聚类个数到底怎么定2.1 KMeans迭代过程的工程视角KMeans的数学定义谁都背得出来但在工程实现里真正影响结果的是迭代过程中的细节。算法的流程是这样的先从样本中随机挑K个点作为初始聚类中心然后计算每个样本到所有中心的欧氏距离把它们划归到最近的中心接着重新计算每个簇内所有样本的均值把均值点作为新的中心重复这两步直到聚类中心的变化量小于某个阈值或达到最大迭代次数。从工程角度看这里有几个直接影响结果质量的参数。首先是距离度量方式KMeans默认用的是欧氏距离这意味着所有特征都必须在同一个量纲尺度下参与计算否则数值范围大的特征会完全主导距离计算。其次是初始中心的选取MATLAB的kmeans函数默认使用kmeans算法来初始化中心它会让初始中心彼此尽可能分散避免收敛到局部最优解。最后是迭代终止条件MaxIter和Replicates这两个参数分别控制最大迭代轮数和重复运行的次数。% 典型KMeans调用方式 % X: m×n矩阵m为样本数n为特征数 % k: 聚类数 % idx: 每个样本所属簇的编号m×1向量 % C: 最终聚类中心坐标k×n矩阵 [idx, C] kmeans(X, k, Distance, sqeuclidean, ... MaxIter, 500, Replicates, 5, Start, plus);这段代码里Replicates设为5意思是算法会从不同初始中心开始独立跑5次只保留误差平方和最小的一次结果。MaxIter设为500给足迭代次数以免提前收敛。这两项是KMeans工程落地最实用的后悔药——前者对抗随机性后者对抗不收敛。2.2 Calinski-Harabasz指数不用拍脑袋选KKMeans最大的痛点是K值必须提前给定。故障分析场景下你根本不知道设备到底有几种典型状态拍脑袋设一个K3可能把两种故障混成一个簇也可能把一个连续退化过程强行拆碎。Calinski-Harabasz指数CH指数就是为了解决这个问题。它的计算逻辑是类内平方和Within-cluster sum of squares, WSS衡量样本到所属簇中心的距离总和这个值越小说明簇内越紧凑类间平方和Between-cluster sum of squares, BSS衡量各簇中心到全局中心的加权距离这个值越大说明不同簇之间分得越开。CH指数就是这两者的比值经过自由度修正CH(k) [BSS / (k-1)] / [WSS / (n-k)]n是样本总数k是当前聚类数。计算方式不需要自己写MATLAB的evalclusters函数直接内置了CalinskiHarabasz评价对象。实际操作中把K从2扫到10对每个K跑一遍KMeans计算出对应的CH指数取最大值的K就是数据结构的最佳分割数。% 自动评估最优聚类数 % X: 标准化后的样本矩阵 eva evalclusters(X, kmeans, CalinskiHarabasz, ... KList, 2:10); % KList: 要评估的聚类数范围 % eva.OptimalK: 自动选取的最佳聚类个数 % eva.CriterionValues: 每个K对应的CH数值 optimal_k eva.OptimalK;这段代码跑完之后optimal_k就是脚本自动选出的最佳聚类数。CriterionValues里存着每个K对应的CH指数可以用来画折线图观察CH值随K的变化趋势。如果曲线在某个K值处出现明显的尖峰说明数据在这个分割数下结构最清晰如果曲线单调上升没有尖峰说明数据本身不存在天然的簇结构这时候强行聚类没有意义。2.3 为什么CH指数比肘部法则靠谱肘部法则Elbow Method看的是WSS随K变化的拐点但拐点是一个主观判断的结果——不同的人看同一张图可能选出不同的K。CH指数给的是一个数值化的标准K从2增加时如果CH值持续增长说明多分出来的簇确实让数据结构更清晰当CH值开始回落说明多余的分割只是在把本来就集中的数据强行拆开这就该停了。在故障分析里这个客观标准尤为重要因为故障类型划分直接影响后续的维护决策主观性必须被压缩到最低。3. 数据预处理与导入流程把Excel里的故障样本喂给KMeans3.1 xlsx数据文件的读取与清洗资源包中的abnormal_heating_fragments11111.xlsx是一个典型的故障数据工作簿存储设备异常加热过程的片段信息。这类数据的常见形态是行表示一个时间片段或一次事件列表示不同的监测特征——温度、温升速率、电流、振动幅值、运行时长等。直接用readtable读取后第一件事不是送进聚类算法而是做数据体检。故障数据最常见的坑有三个缺失值、离群点和重复行。缺失值会让距离计算失效离群点会因为欧氏距离对极值的敏感性而扭曲聚类中心重复行则会在总量中产生权重偏差。% 读取Excel数据 clear; clc; T readtable(abnormal_heating_fragments11111.xlsx); % 检查缺失值 fprintf(缺失值统计:\n); disp(ismissing(T)); % 删除包含缺失值的行 T_clean rmmissing(T); % 删除完全重复的行 T_unique unique(T_clean); % 提取数值矩阵 X_raw table2array(T_unique); fprintf(清洗后数据维度: %d 行 × %d 列\n, size(X_raw, 1), size(X_raw, 2));rmmissing直接删掉含NaN的行unique在这里能把完全一致的行去重。这两步做完数据维度通常会缩水一些这是正常的——合格的样本集比体面的样本量更重要。3.2 标准化KMeans里最容易被跳过的关键步骤如果直接拿原始数据跑KMeans结果大概率是灾难性的。比如温度特征的数值在几十到几百之间而电流百分比的数值在0到1之间欧氏距离几乎完全被温度主导电流信息形同虚设。这就是为什么标准化这一步在聚类流程中不可跳过。工程上最常用的标准化方法是Z-score标准化用每列均值减除后除以标准差让所有特征都落在近似均值为0、标准差为1的分布上。MATLAB里的zscore函数一行搞定。有一种常见误用是把zscore写成了normalize且不指定方法这会导致默认范围映射而不是Z-score结果完全不同。% Z-score标准化 % 核心目的消除量纲差异避免大数值特征主导距离计算 X_std zscore(X_raw); % 标准化后的数据均值为0标准差为1 assert(abs(mean(X_std(:, 1))) 1e-10, 标准化失败均值不为0); assert(abs(std(X_std(:, 1)) - 1) 1e-10, 标准化失败标准差不为1);标准化之后的X_std才能作为KMeans的输入。如果原始数据存在明显的量纲差异但跳过标准化聚类结果往往表现为一个簇模板是大数值特征的堆叠其他簇全是小数值特征的簇拥簇间区分度极差CH指数也高不了。3.3 数据维度检查与特征筛选拿到X_raw之后还应该看一下数据的数值范围这决定了是否需要额外的处理。用summary或min/max扫一遍各列的分布范围如果某个特征的方差几乎为零说明它在这个数据集里没有区分能力留着只会增加计算量。直接剔除即可。% 检查各列方差剔除零方差特征 variance_vals var(X_std); zero_var_cols find(variance_vals 1e-8); if ~isempty(zero_var_cols) fprintf(剔除零方差特征列: %d\n, zero_var_cols); X_std(:, zero_var_cols) []; end这一步不算必须但做一次能避免一种隐蔽的问题某些传感器在故障期间记录的数据恒定为某个值这种列对聚类没有贡献却会在距离计算中引入噪声。特征选择的目标是让每个参与计算的特征都携带区分信息。4. 聚类执行与可视化从散点图到故障模式判读4.1 运行KMeans并获得样本标签在标准化矩阵X_std和自动选出的最佳聚类数optimal_k都已就绪后正式执行KMeans聚类。这里要注意Replicates建议至少设5因为KMeans的初始中心选择有随机性不同随机种子可能收敛到不同的局部最优解多跑几次保留最优结果能显著提升稳定性。% 执行KMeans聚类 % 注意此时X_std是标准化后的数据optimal_k由CH指数自动确定 rng(42); % 固定随机种子保证结果可复现 final_idx kmeans(X_std, optimal_k, ... Distance, sqeuclidean, ... MaxIter, 1000, ... Replicates, 10, ... Start, plus); % 统计每个簇的样本数 cluster_counts histcounts(final_idx, 1:optimal_k1); fprintf(最佳聚类数: %d\n, optimal_k); for k 1:optimal_k fprintf(簇 %d: %d 个样本\n, k, cluster_counts(k)); endrng(42)这一行很关键。故障分析项目往往需要向别人展示或复现结果如果不固定随机种子每次运行得到的标签顺序可能不一样——虽然聚类结构大体相同但样本点编号会变给后续排查工作带来困扰。4.2 二维散点图与三维投影的可视化实现聚类结果的可视化是让故障模式显形的关键一步。原始数据可能有五六个甚至十几个特征直接画高维图不现实通常做法是取最重要的两个主成分做二维散点图或者取三个主成分做三维图。主成分分析PCA在这里不是算法步骤而是可视化工具——它能把高维数据投影到低维空间并尽量保留方差结构。MATLAB内置的pca函数可以直接使用。投影完成后不同簇用不同颜色的散点标记出来横纵轴标注为 PC1/PC2。% 降维到二维用于可视化 % score投影后的主成分得分第一列PC1第二列PC2 % latent各主成分的解释方差占比 [coeff, score, latent] pca(X_std); % 二维散点图每个簇一个颜色 figure(Position, [100, 100, 900, 600]); hold on; % 定义颜色映射最多支持到10个簇 colorMap lines(optimal_k); for k 1:optimal_k mask (final_idx k); scatter(score(mask, 1), score(mask, 2), 40, ... colorMap(k, :), filled, MarkerFaceAlpha, 0.7); end grid on; xlabel(PC1); ylabel(PC2); title(sprintf(KMeans聚类结果可视化 (K%d), optimal_k)); legend(cellstr(num2str((1:optimal_k), 簇 %d))); hold off;lines(optimal_k)自动生成一组区分度较好的颜色序列40是散点尺寸MarkerFaceAlpha设为0.7让重叠点呈半透明状避免大片覆盖。图例用簇 1、簇 2这种编号比直接标故障类型A更稳妥——因为聚类结果本身并不知道每个簇对应的故障名称标记编号能防止误导。4.3 各簇特征的统计分析散点图能直观看出簇的分割情况但要真正理解每个簇代表什么状态还需要回到原始特征空间去看各簇的统计特征。这个时候箱线图最合适对每个聚类簇绘制其关键特征的分布。暖温、电流、振动这类的分布差异直接决定了这个簇是正常工况还是异常加热。% 箱线图对比各簇在某个原始特征上的分布 % 以X_raw第2列假设为温度特征为例 figure(Position, [150, 150, 800, 500]); feature_idx 2; % 换成实际温度对应列 boxplot(X_raw(:, feature_idx), final_idx, ... Labels, cellstr(num2str((1:optimal_k), 簇 %d))); xlabel(聚类簇); ylabel(温度特征值); title(各簇温度分布对比); grid on;把每个特征列都跑一遍这样的箱线图就能逐一确认簇1温度中位数明显偏高温升速率也快判定为异常加热模式簇2温度平缓、电流稳定可能是正常运行工况簇3如果特征混杂说明这一簇内部还有细分空间可以考虑在局部数据上再做一次聚类。这种先分簇、再翻译的流程是把数学结果转成设备诊断结论的桥梁。5. KMeans聚类避坑指南数据标准化、空簇与初始中心的三类翻车5.1 翻车现场一聚类结果全是大数特征的天下现象聚类完成之后画散点图发现样本被按照某个单一特征从低到高切成了几条带状分布簇的形状又窄又长跟预想的圆形簇完全不沾边。原因没做标准化。欧氏距离对量纲敏感数值范围大的特征在距离计算中权重过高其他特征的信息被完全淹没。这种情况在故障数据里尤其常见——温度动辄几十上百电流可能只有零点几距离几乎完全被温度主导。解决跑聚类之前强制执行X_std zscore(X_raw)并且用assert验证标准化后的均值和标准差确保执行正确。我在自己的项目里已经把这个检查写进固定流程宁可多写两行断言也不省这一步。5.2 翻车现场二某次聚类运行出现空簇现象histcounts统计结果显示某个簇的样本数为零聚类数明明是5画图的时候只剩4个可见的簇。原因初始中心选得不好两个中心初始位置过于接近迭代过程中其中一个中心被挤到了另外一个簇导致一个簇被吃掉。KMeans对初始中心非常敏感每次运行的结果都可能不一样。解决Replicates参数从默认的1提到5以上。MATLAB的kmeans会独立运行多次每次都使用不同的初始中心集合最后选择误差平方和最小的结果。同时把Start指定为plus来启用 kmeans 初始化让初始中心彼此尽可能分散。这两种手段搭配使用空簇概率能压到非常低。5.3 翻车现场三每次运行聚类结果都不一样现象同一个脚本同一个数据集今天跑出来簇1有200个样本明天跑出来簇1有153个样本簇内样本成员也变了导致向领导汇报的结果无法复现。原因KMeans的初始中心是随机生成的即使Replicates设为10不同随机种子下跑出来的最优结果也可能是不同的局部最优解。解决在脚本开头固定随机种子写上rng(42)这一行。42只是一个任意值任何固定数字都可以关键是让随机数生成器的起点固定下来。这样一来每次运行脚本都会经历相同的随机序列获得完全一致的结果。同一份数据同一份脚本任何人运行都能得到相同输出——这在故障诊断项目交付时是硬性要求。5.4 翻车现场四CH指数选出的K不符合业务预期现象evalclusters自动选出的最佳K是8但业务上故障类型只分了4种多出来的簇不知道对应什么物理意义。原因CH指数是在纯数学层面寻找最优分割它不知道你的业务约束。有些故障模式在特征空间里是渐进变化的同一个故障可以被数据自然切成多个子簇。解决CH指数是最佳聚类数的参考但不是唯一标准。看CriterionValues的完整曲线如果K4和K8的CH值差距不到10%选4是明智的业务决策。另一种做法是做完K8的聚类后检查每个簇的中心特征差异差异极小的簇可以考虑合并。聚类是为业务服务的不是为指标服务的。5.5 翻车现场五Excel数据读入后全是NaN现象readtable读入后ismissing显示几乎所有单元格都是缺失值数据矩阵没法用。原因最常见的是文件的单元格区域不规范比如表头在第二行、表头有合并单元格、或者工作簿里有多个sheet没选对。readtable默认读第一个sheet如果你的数据在第二个sheet上就会出错。解决用readtable时显式指定Sheet参数并且用VariableNamingRule处理非标准列名。读取后先用head(T)看前几行数据是否正常再往下走。我的习惯是任何从Excel读数据的脚本读取之后第一步永远是打印行列数和前几行内容确认数据形状正确之后才继续。6. 自动化K值选择脚本把CH指数曲线和最优聚类串成一条命令前面的流程已经把每一步拆开了但如果每次换数据都要手动改参数、重新跑一遍效率会很低。我通常会把这些步骤封装成一个内部函数输入标准化数据矩阵输出最优聚类数和聚类标签同时自动画出CH指数曲线和聚类散点图。function [best_k, idx, criteria] auto_kmeans_cluster(X_std, k_range) % 输入: % X_std - 标准化后的样本矩阵m×n % k_range - 要搜索的聚类数范围如2:10 % 输出: % best_k - CH指数确定的最佳聚类数 % idx - 最佳聚类数下的样本标签 % criteria - 每个K对应的CH值用于绘图 % 第一步用CH指数评估所有候选K eva evalclusters(X_std, kmeans, CalinskiHarabasz, ... KList, k_range); best_k eva.OptimalK; criteria eva.CriterionValues; % 第二步在最佳K下执行最终聚类 rng(42); % 固定种子保证可复现 idx kmeans(X_std, best_k, ... Distance, sqeuclidean, ... MaxIter, 1000, ... Replicates, 10, ... Start, plus); % 第三步绘制CH指数曲线 figure(Position, [50, 50, 800, 400]); plot(k_range, criteria, bo-, LineWidth, 1.5); hold on; plot(best_k, criteria(k_range best_k), r*, ... MarkerSize, 12); xlabel(聚类数 K); ylabel(Calinski-Harabasz 指数); title(最优K值选择); grid on; fprintf(Optimal K %d\n, best_k); end这个函数把评估、聚类、绘图三步绑定在一起日常使用时只需要一行调用[best_k, idx] auto_kmeans_cluster(X_std, 2:10);整个流程跑完最优聚类数和聚类标签就直接到手了。不过要注意这个函数有一个隐含假设最佳K在k_range范围内。如果你的数据特征复杂可能最优K会落在k_range的边界上——比如K10的时候CH值还在涨那说明你给的上界设小了需要扩大范围重新跑。另外聚类数是整数evalclusters给出的结果从数学上是最高CH对应的K但如果两条曲线走势接近还是回到第5章的教训结合实际业务做决断别被数字牵着走。这套流程我用了很多次最深的感触是跑通一个聚类分析轻而易举但让聚类结果在业务上站得住脚九十以上的功夫花在数据清洗、标准化的边界确认和聚类数的合理性验证上。从那以后我每次拿到一批新的故障数据第一件事必定是先做标准化和缺失值扫描然后带着全程复现的随机种子跑CH评估。这样出来的结论给你能复现给别人也能复现才算是真正落在实处的分析。希望帮到你也欢迎带着你的数据在评论区验证这套KMeans流程的边界。本文还有配套的精品资源点击获取 KMeans聚类实战：从标准化到CH指数自动确定故障分类
简介这是一套K均值聚类分析的可执行源码与演示数据包面向机器学习初学者和故障诊断研究人员。针对故障类型未知的问题利用K均值自动划分样本并以CH指数筛选最佳聚类数避免人工指定带来的主观性。压缩包内共2个文件一个M脚本负责聚类实现、指数计算与可视化结果输出一个电子表格存放异常加热片段等原始数据包体约57KB数据经简单预处理后可直接作为算法输入便于快速复现。已有4508人学习/浏览适合通过实例掌握K均值聚类、聚类数评估以及散点图和统计图表展示流程无需额外配置环境可用于课程设计、毕业设计或入门实践也适合在故障样本类别数量不详时快速建立聚类分析流程。代码结构清晰可以直接运行能输出不同类别的分布图可作为故障分类和模式识别项目的参考实现。1. KMeans聚类故障类型分析里那些没人手写的分类规则在设备故障诊断场景里最头疼的事情不是数据不够而是你根本不知道眼前的故障到底该分成几类。轴承磨损、电气老化、冷却异常、负载突变这些故障模式在波形上互相交叠靠人眼去分类完全不可行。KMeans聚类作为无监督学习中应用最广的算法不需要任何标签直接对样本的分布结构下手把相似的数据点自动归堆。这份资源包里的Copy_of_Cluster_Flow.m脚本配合ML_Cluster.zip中的数据文件就是一个完整的故障数据聚类分析流程——从读入Excel文件、标准化处理、自动评估最佳聚类数到最终把分类结果以不同颜色的散点图呈现出来。适合那些拿到了设备监测数据但不知道从哪儿开始做分类的工程师也适合正在学聚类算法、需要一个能跑通的完整例子的学生。2. KMeans核心原理与Calinski-Harabasz指数聚类个数到底怎么定2.1 KMeans迭代过程的工程视角KMeans的数学定义谁都背得出来但在工程实现里真正影响结果的是迭代过程中的细节。算法的流程是这样的先从样本中随机挑K个点作为初始聚类中心然后计算每个样本到所有中心的欧氏距离把它们划归到最近的中心接着重新计算每个簇内所有样本的均值把均值点作为新的中心重复这两步直到聚类中心的变化量小于某个阈值或达到最大迭代次数。从工程角度看这里有几个直接影响结果质量的参数。首先是距离度量方式KMeans默认用的是欧氏距离这意味着所有特征都必须在同一个量纲尺度下参与计算否则数值范围大的特征会完全主导距离计算。其次是初始中心的选取MATLAB的kmeans函数默认使用kmeans算法来初始化中心它会让初始中心彼此尽可能分散避免收敛到局部最优解。最后是迭代终止条件MaxIter和Replicates这两个参数分别控制最大迭代轮数和重复运行的次数。% 典型KMeans调用方式 % X: m×n矩阵m为样本数n为特征数 % k: 聚类数 % idx: 每个样本所属簇的编号m×1向量 % C: 最终聚类中心坐标k×n矩阵 [idx, C] kmeans(X, k, Distance, sqeuclidean, ... MaxIter, 500, Replicates, 5, Start, plus);这段代码里Replicates设为5意思是算法会从不同初始中心开始独立跑5次只保留误差平方和最小的一次结果。MaxIter设为500给足迭代次数以免提前收敛。这两项是KMeans工程落地最实用的后悔药——前者对抗随机性后者对抗不收敛。2.2 Calinski-Harabasz指数不用拍脑袋选KKMeans最大的痛点是K值必须提前给定。故障分析场景下你根本不知道设备到底有几种典型状态拍脑袋设一个K3可能把两种故障混成一个簇也可能把一个连续退化过程强行拆碎。Calinski-Harabasz指数CH指数就是为了解决这个问题。它的计算逻辑是类内平方和Within-cluster sum of squares, WSS衡量样本到所属簇中心的距离总和这个值越小说明簇内越紧凑类间平方和Between-cluster sum of squares, BSS衡量各簇中心到全局中心的加权距离这个值越大说明不同簇之间分得越开。CH指数就是这两者的比值经过自由度修正CH(k) [BSS / (k-1)] / [WSS / (n-k)]n是样本总数k是当前聚类数。计算方式不需要自己写MATLAB的evalclusters函数直接内置了CalinskiHarabasz评价对象。实际操作中把K从2扫到10对每个K跑一遍KMeans计算出对应的CH指数取最大值的K就是数据结构的最佳分割数。% 自动评估最优聚类数 % X: 标准化后的样本矩阵 eva evalclusters(X, kmeans, CalinskiHarabasz, ... KList, 2:10); % KList: 要评估的聚类数范围 % eva.OptimalK: 自动选取的最佳聚类个数 % eva.CriterionValues: 每个K对应的CH数值 optimal_k eva.OptimalK;这段代码跑完之后optimal_k就是脚本自动选出的最佳聚类数。CriterionValues里存着每个K对应的CH指数可以用来画折线图观察CH值随K的变化趋势。如果曲线在某个K值处出现明显的尖峰说明数据在这个分割数下结构最清晰如果曲线单调上升没有尖峰说明数据本身不存在天然的簇结构这时候强行聚类没有意义。2.3 为什么CH指数比肘部法则靠谱肘部法则Elbow Method看的是WSS随K变化的拐点但拐点是一个主观判断的结果——不同的人看同一张图可能选出不同的K。CH指数给的是一个数值化的标准K从2增加时如果CH值持续增长说明多分出来的簇确实让数据结构更清晰当CH值开始回落说明多余的分割只是在把本来就集中的数据强行拆开这就该停了。在故障分析里这个客观标准尤为重要因为故障类型划分直接影响后续的维护决策主观性必须被压缩到最低。3. 数据预处理与导入流程把Excel里的故障样本喂给KMeans3.1 xlsx数据文件的读取与清洗资源包中的abnormal_heating_fragments11111.xlsx是一个典型的故障数据工作簿存储设备异常加热过程的片段信息。这类数据的常见形态是行表示一个时间片段或一次事件列表示不同的监测特征——温度、温升速率、电流、振动幅值、运行时长等。直接用readtable读取后第一件事不是送进聚类算法而是做数据体检。故障数据最常见的坑有三个缺失值、离群点和重复行。缺失值会让距离计算失效离群点会因为欧氏距离对极值的敏感性而扭曲聚类中心重复行则会在总量中产生权重偏差。% 读取Excel数据 clear; clc; T readtable(abnormal_heating_fragments11111.xlsx); % 检查缺失值 fprintf(缺失值统计:\n); disp(ismissing(T)); % 删除包含缺失值的行 T_clean rmmissing(T); % 删除完全重复的行 T_unique unique(T_clean); % 提取数值矩阵 X_raw table2array(T_unique); fprintf(清洗后数据维度: %d 行 × %d 列\n, size(X_raw, 1), size(X_raw, 2));rmmissing直接删掉含NaN的行unique在这里能把完全一致的行去重。这两步做完数据维度通常会缩水一些这是正常的——合格的样本集比体面的样本量更重要。3.2 标准化KMeans里最容易被跳过的关键步骤如果直接拿原始数据跑KMeans结果大概率是灾难性的。比如温度特征的数值在几十到几百之间而电流百分比的数值在0到1之间欧氏距离几乎完全被温度主导电流信息形同虚设。这就是为什么标准化这一步在聚类流程中不可跳过。工程上最常用的标准化方法是Z-score标准化用每列均值减除后除以标准差让所有特征都落在近似均值为0、标准差为1的分布上。MATLAB里的zscore函数一行搞定。有一种常见误用是把zscore写成了normalize且不指定方法这会导致默认范围映射而不是Z-score结果完全不同。% Z-score标准化 % 核心目的消除量纲差异避免大数值特征主导距离计算 X_std zscore(X_raw); % 标准化后的数据均值为0标准差为1 assert(abs(mean(X_std(:, 1))) 1e-10, 标准化失败均值不为0); assert(abs(std(X_std(:, 1)) - 1) 1e-10, 标准化失败标准差不为1);标准化之后的X_std才能作为KMeans的输入。如果原始数据存在明显的量纲差异但跳过标准化聚类结果往往表现为一个簇模板是大数值特征的堆叠其他簇全是小数值特征的簇拥簇间区分度极差CH指数也高不了。3.3 数据维度检查与特征筛选拿到X_raw之后还应该看一下数据的数值范围这决定了是否需要额外的处理。用summary或min/max扫一遍各列的分布范围如果某个特征的方差几乎为零说明它在这个数据集里没有区分能力留着只会增加计算量。直接剔除即可。% 检查各列方差剔除零方差特征 variance_vals var(X_std); zero_var_cols find(variance_vals 1e-8); if ~isempty(zero_var_cols) fprintf(剔除零方差特征列: %d\n, zero_var_cols); X_std(:, zero_var_cols) []; end这一步不算必须但做一次能避免一种隐蔽的问题某些传感器在故障期间记录的数据恒定为某个值这种列对聚类没有贡献却会在距离计算中引入噪声。特征选择的目标是让每个参与计算的特征都携带区分信息。4. 聚类执行与可视化从散点图到故障模式判读4.1 运行KMeans并获得样本标签在标准化矩阵X_std和自动选出的最佳聚类数optimal_k都已就绪后正式执行KMeans聚类。这里要注意Replicates建议至少设5因为KMeans的初始中心选择有随机性不同随机种子可能收敛到不同的局部最优解多跑几次保留最优结果能显著提升稳定性。% 执行KMeans聚类 % 注意此时X_std是标准化后的数据optimal_k由CH指数自动确定 rng(42); % 固定随机种子保证结果可复现 final_idx kmeans(X_std, optimal_k, ... Distance, sqeuclidean, ... MaxIter, 1000, ... Replicates, 10, ... Start, plus); % 统计每个簇的样本数 cluster_counts histcounts(final_idx, 1:optimal_k1); fprintf(最佳聚类数: %d\n, optimal_k); for k 1:optimal_k fprintf(簇 %d: %d 个样本\n, k, cluster_counts(k)); endrng(42)这一行很关键。故障分析项目往往需要向别人展示或复现结果如果不固定随机种子每次运行得到的标签顺序可能不一样——虽然聚类结构大体相同但样本点编号会变给后续排查工作带来困扰。4.2 二维散点图与三维投影的可视化实现聚类结果的可视化是让故障模式显形的关键一步。原始数据可能有五六个甚至十几个特征直接画高维图不现实通常做法是取最重要的两个主成分做二维散点图或者取三个主成分做三维图。主成分分析PCA在这里不是算法步骤而是可视化工具——它能把高维数据投影到低维空间并尽量保留方差结构。MATLAB内置的pca函数可以直接使用。投影完成后不同簇用不同颜色的散点标记出来横纵轴标注为 PC1/PC2。% 降维到二维用于可视化 % score投影后的主成分得分第一列PC1第二列PC2 % latent各主成分的解释方差占比 [coeff, score, latent] pca(X_std); % 二维散点图每个簇一个颜色 figure(Position, [100, 100, 900, 600]); hold on; % 定义颜色映射最多支持到10个簇 colorMap lines(optimal_k); for k 1:optimal_k mask (final_idx k); scatter(score(mask, 1), score(mask, 2), 40, ... colorMap(k, :), filled, MarkerFaceAlpha, 0.7); end grid on; xlabel(PC1); ylabel(PC2); title(sprintf(KMeans聚类结果可视化 (K%d), optimal_k)); legend(cellstr(num2str((1:optimal_k), 簇 %d))); hold off;lines(optimal_k)自动生成一组区分度较好的颜色序列40是散点尺寸MarkerFaceAlpha设为0.7让重叠点呈半透明状避免大片覆盖。图例用簇 1、簇 2这种编号比直接标故障类型A更稳妥——因为聚类结果本身并不知道每个簇对应的故障名称标记编号能防止误导。4.3 各簇特征的统计分析散点图能直观看出簇的分割情况但要真正理解每个簇代表什么状态还需要回到原始特征空间去看各簇的统计特征。这个时候箱线图最合适对每个聚类簇绘制其关键特征的分布。暖温、电流、振动这类的分布差异直接决定了这个簇是正常工况还是异常加热。% 箱线图对比各簇在某个原始特征上的分布 % 以X_raw第2列假设为温度特征为例 figure(Position, [150, 150, 800, 500]); feature_idx 2; % 换成实际温度对应列 boxplot(X_raw(:, feature_idx), final_idx, ... Labels, cellstr(num2str((1:optimal_k), 簇 %d))); xlabel(聚类簇); ylabel(温度特征值); title(各簇温度分布对比); grid on;把每个特征列都跑一遍这样的箱线图就能逐一确认簇1温度中位数明显偏高温升速率也快判定为异常加热模式簇2温度平缓、电流稳定可能是正常运行工况簇3如果特征混杂说明这一簇内部还有细分空间可以考虑在局部数据上再做一次聚类。这种先分簇、再翻译的流程是把数学结果转成设备诊断结论的桥梁。5. KMeans聚类避坑指南数据标准化、空簇与初始中心的三类翻车5.1 翻车现场一聚类结果全是大数特征的天下现象聚类完成之后画散点图发现样本被按照某个单一特征从低到高切成了几条带状分布簇的形状又窄又长跟预想的圆形簇完全不沾边。原因没做标准化。欧氏距离对量纲敏感数值范围大的特征在距离计算中权重过高其他特征的信息被完全淹没。这种情况在故障数据里尤其常见——温度动辄几十上百电流可能只有零点几距离几乎完全被温度主导。解决跑聚类之前强制执行X_std zscore(X_raw)并且用assert验证标准化后的均值和标准差确保执行正确。我在自己的项目里已经把这个检查写进固定流程宁可多写两行断言也不省这一步。5.2 翻车现场二某次聚类运行出现空簇现象histcounts统计结果显示某个簇的样本数为零聚类数明明是5画图的时候只剩4个可见的簇。原因初始中心选得不好两个中心初始位置过于接近迭代过程中其中一个中心被挤到了另外一个簇导致一个簇被吃掉。KMeans对初始中心非常敏感每次运行的结果都可能不一样。解决Replicates参数从默认的1提到5以上。MATLAB的kmeans会独立运行多次每次都使用不同的初始中心集合最后选择误差平方和最小的结果。同时把Start指定为plus来启用 kmeans 初始化让初始中心彼此尽可能分散。这两种手段搭配使用空簇概率能压到非常低。5.3 翻车现场三每次运行聚类结果都不一样现象同一个脚本同一个数据集今天跑出来簇1有200个样本明天跑出来簇1有153个样本簇内样本成员也变了导致向领导汇报的结果无法复现。原因KMeans的初始中心是随机生成的即使Replicates设为10不同随机种子下跑出来的最优结果也可能是不同的局部最优解。解决在脚本开头固定随机种子写上rng(42)这一行。42只是一个任意值任何固定数字都可以关键是让随机数生成器的起点固定下来。这样一来每次运行脚本都会经历相同的随机序列获得完全一致的结果。同一份数据同一份脚本任何人运行都能得到相同输出——这在故障诊断项目交付时是硬性要求。5.4 翻车现场四CH指数选出的K不符合业务预期现象evalclusters自动选出的最佳K是8但业务上故障类型只分了4种多出来的簇不知道对应什么物理意义。原因CH指数是在纯数学层面寻找最优分割它不知道你的业务约束。有些故障模式在特征空间里是渐进变化的同一个故障可以被数据自然切成多个子簇。解决CH指数是最佳聚类数的参考但不是唯一标准。看CriterionValues的完整曲线如果K4和K8的CH值差距不到10%选4是明智的业务决策。另一种做法是做完K8的聚类后检查每个簇的中心特征差异差异极小的簇可以考虑合并。聚类是为业务服务的不是为指标服务的。5.5 翻车现场五Excel数据读入后全是NaN现象readtable读入后ismissing显示几乎所有单元格都是缺失值数据矩阵没法用。原因最常见的是文件的单元格区域不规范比如表头在第二行、表头有合并单元格、或者工作簿里有多个sheet没选对。readtable默认读第一个sheet如果你的数据在第二个sheet上就会出错。解决用readtable时显式指定Sheet参数并且用VariableNamingRule处理非标准列名。读取后先用head(T)看前几行数据是否正常再往下走。我的习惯是任何从Excel读数据的脚本读取之后第一步永远是打印行列数和前几行内容确认数据形状正确之后才继续。6. 自动化K值选择脚本把CH指数曲线和最优聚类串成一条命令前面的流程已经把每一步拆开了但如果每次换数据都要手动改参数、重新跑一遍效率会很低。我通常会把这些步骤封装成一个内部函数输入标准化数据矩阵输出最优聚类数和聚类标签同时自动画出CH指数曲线和聚类散点图。function [best_k, idx, criteria] auto_kmeans_cluster(X_std, k_range) % 输入: % X_std - 标准化后的样本矩阵m×n % k_range - 要搜索的聚类数范围如2:10 % 输出: % best_k - CH指数确定的最佳聚类数 % idx - 最佳聚类数下的样本标签 % criteria - 每个K对应的CH值用于绘图 % 第一步用CH指数评估所有候选K eva evalclusters(X_std, kmeans, CalinskiHarabasz, ... KList, k_range); best_k eva.OptimalK; criteria eva.CriterionValues; % 第二步在最佳K下执行最终聚类 rng(42); % 固定种子保证可复现 idx kmeans(X_std, best_k, ... Distance, sqeuclidean, ... MaxIter, 1000, ... Replicates, 10, ... Start, plus); % 第三步绘制CH指数曲线 figure(Position, [50, 50, 800, 400]); plot(k_range, criteria, bo-, LineWidth, 1.5); hold on; plot(best_k, criteria(k_range best_k), r*, ... MarkerSize, 12); xlabel(聚类数 K); ylabel(Calinski-Harabasz 指数); title(最优K值选择); grid on; fprintf(Optimal K %d\n, best_k); end这个函数把评估、聚类、绘图三步绑定在一起日常使用时只需要一行调用[best_k, idx] auto_kmeans_cluster(X_std, 2:10);整个流程跑完最优聚类数和聚类标签就直接到手了。不过要注意这个函数有一个隐含假设最佳K在k_range范围内。如果你的数据特征复杂可能最优K会落在k_range的边界上——比如K10的时候CH值还在涨那说明你给的上界设小了需要扩大范围重新跑。另外聚类数是整数evalclusters给出的结果从数学上是最高CH对应的K但如果两条曲线走势接近还是回到第5章的教训结合实际业务做决断别被数字牵着走。这套流程我用了很多次最深的感触是跑通一个聚类分析轻而易举但让聚类结果在业务上站得住脚九十以上的功夫花在数据清洗、标准化的边界确认和聚类数的合理性验证上。从那以后我每次拿到一批新的故障数据第一件事必定是先做标准化和缺失值扫描然后带着全程复现的随机种子跑CH评估。这样出来的结论给你能复现给别人也能复现才算是真正落在实处的分析。希望帮到你也欢迎带着你的数据在评论区验证这套KMeans流程的边界。本文还有配套的精品资源点击获取