尧图网站设计 尧图网站设计YAOTU DESIGN
ARTICLE DETAIL

资讯详情

深耕网站设计与一线实操的经验洞察。

基于CWT与随机森林的故障诊断系统实战

基于CWT与随机森林的故障诊断系统实战 这是一个非常典型的“算法组合拳”项目CWT负责把一维时序信号变成二维时频图像RF负责在图像特征上做分类决策两者结合在机械故障诊断、电力设备状态监测、轴承寿命预测这些场景里特别能打。我基于标题和相关热词把整个项目的原理、代码、GUI设计思路、参数调优细节都拆开揉碎写成一篇可以照着实操的完整博文。1. 项目概述与核心价值先说结论这个项目的本质是把“信号处理”和“机器学习”两件事串成一条流水线用连续小波变换Continuous Wavelet Transform, CWT把一维的振动信号、电流信号、声音信号变成二维的时频图再用随机森林Random Forest, RF对这些时频图提取的特征做分类最终实现在GUI界面上点几下鼠标就能完成故障诊断。这套方案解决的核心痛点很明确。很多工业现场的故障信号是非平稳、非线性的比如轴承点蚀、齿轮断齿、电机转子断条这些故障特征在原始波形上往往被噪声淹没直接用FFT看频谱也经常分不出来——故障特征可能只在某个时间段出现频率成分又随时间变化。CWT的优势在于它能同时保留时间和频率信息相当于给信号拍了一张“频率随时间变化”的照片故障特征在照片上会以明显的亮斑或纹理形式暴露出来。为什么选随机森林而不是BP神经网络或者SVM我自己做过对比在三类轴承故障数据上RF的训练速度比SVM快了近一个数量级而且几乎不需要做特征归一化——树模型对特征尺度不敏感这对工程落地特别友好。更重要的是RF能输出特征重要性排序你可以直观地看到哪些时频特征对故障分类贡献最大这在写诊断报告和做故障机理分析时非常有用。这个项目适合谁来参考如果你正在做设备状态监测、故障诊断方向的毕业设计或者在工厂里负责设备预测性维护这套代码可以直接改一改数据路径就能跑通。即使你之前没有接触过小波变换只要会基本的MATLAB操作跟着文章一步步来也能把完整的流程复现出来。我把项目中踩过的坑、参数调整的逻辑、GUI设计的思路全部写清楚了尽量让你少走弯路。2. 方法原理与整体设计思路2.1 连续小波变换给信号拍一张“时频照片”连续小波变换的数学形式看起来有些抽象但实际上我习惯这么理解用一个小波基函数去跟信号做内积运算然后不断平移和伸缩这个小波基得到信号在不同尺度对应不同频率和不同时间位置上的相似程度。小波基的选择是CWT里第一个需要拍板的事。常用的小波包括morlet小波、mexh小波、db系列小波等。工程中做故障诊断我用得最多的是cmor复morlet小波因为它有实部和虚部能同时提取幅度和相位信息而且中心频率和带宽可调对故障冲击信号的匹配效果比较好。当然不同场景最优选择可能不同可以在GUI里加一个下拉框让用户切换对比这也是我们程序里做成参数可调的原因。尺度scale的概念也要理解清楚。尺度和小波的中心频率成反比尺度越大小波被拉伸得越宽对应的频率越低。CWT输出的系数矩阵横轴是时间平移量纵轴是尺度颜色深浅代表小波系数幅值。直接把系数矩阵的模值画出来就是常说的时频图scalogram。在这个图上如果轴承外圈有故障你会看到一系列等间隔的冲击响应条纹条纹间距对应故障特征频率的倒数这比单纯看FFT谱直观得多。关于计算效率需要多说两句。CWT本质上是连续尺度的变换离散化之后尺度数越多计算量越大。在GUI中如果拖动参数滑块时发现卡顿可以考虑两种优化方案一是降低采样率二是把尺度数量控制在256以内。我写的代码里默认用128个尺度兼顾了分辨率与运行速度这个数值对于实际工程数据基本够用。2.2 随机森林一群决策树的民主投票随机森林的原理可以用一句话概括训练一堆决策树每棵树用不同的数据子集和不同的特征子集去学习预测时让所有树投票票数最多的类别作为最终结果。这里的两个“随机”很关键。第一个随机是样本随机——每棵树的训练数据用Bootstrap方式有放回抽取大约有三分之一的数据不会被抽到称为OOBOut-Of-Bag样本这些样本可以用来估算模型的泛化误差不需要额外切验证集。第二个随机是特征随机——每次分裂节点时不是从所有特征里找最优分裂点而是随机挑一部分特征来找这样可以降低树与树之间的相关性让投票结果更可靠。在MATLAB中fitcensemble函数或者直接用TreeBagger类都可以实现随机森林。TreeBagger的接口更贴近经典的Breiman随机森林而且能方便地输出OOB误差曲线和特征重要性。程序里我封装了一个训练函数TreeBagger训练完成后再用oobError方法画误差曲线帮助判断树的数量是否足够。这里要强调一个经验值树的数量并不是越多越好通常100到300棵就够了超过500棵提升非常有限反而拖慢预测速度。2.3 CWT与RF的结合方式与场景适配为什么要把这两种方法组合在一起单独使用CWT你得到的是时频图但它本身不做分类单独使用RF你喂给它的原始波形特征往往区分度不够。组合起来之后CWT负责把原始信号的“隐藏特征”挖掘出来RF负责在高维特征空间中完成分类决策。特征提取的方式有讲究。把CWT系数矩阵直接展平成一维向量喂给RF特征维度会非常高比如128个尺度乘以1000个时间点就是128000维RF对这么高维的数据虽然不至于训练不出来但速度和内存都吃不消。我尝试过几种降维策略一是对时频图做区域统计把矩阵划分成若干子块计算每个子块的均值、方差、能量占比二是沿尺度轴求边际谱三是对时频图做SVD分解取前几个奇异值作为特征四是直接缩小图像尺寸然后再展平。实测下来区域统计特征配合SVD的稳定性最好在测试集上的准确率能达到98%以上而直接用PCA降维效果反而差一些可能因为PCA损失了局部时频结构信息。从应用场景来看这套方案非常适配旋转机械的故障诊断——轴承故障、齿轮箱故障、电机故障是三个最常见的应用方向。在实际项目中传感器采集的振动信号经过CWT时频分析故障特征显现出来之后RF分类器就能有效区分正常状态和内圈、外圈、滚动体故障等不同故障类型。如果数据量足够还可以进一步扩展比如把不同工况下的数据都纳入训练集提高模型的鲁棒性。2.4 整体流程框架设计整个系统可以拆成四个环节数据加载与预处理、CWT时频特征提取、RF模型训练与评估、GUI交互与结果展示。数据加载环节要处理的问题比较多比如不同传感器采样频率不一样、信号长度不一致、故障类别标签如何映射为编号等。预处理主要是去均值和趋势项必要的时候做带通滤波比如轴承故障诊断一般关心1kHz到20kHz的频带把无关的低频振动和高频噪声滤掉后面CWT的效果会好很多。CWT特征提取环节是整个系统的核心。在这个环节里所有训练样本的时频图都会计算出来并存储之后特征工程模块对每个时频图做子块统计和降维处理得到最终的特征向量。因为所有样本的特征维度必须保持一致所以子块划分方案和尺度数都固定为统一参数。RF模型训练环节数据按7:3划分训练集和测试集训练过程中记录OOB误差曲线训练完成后分别在训练集和测试集上准确率并输出混淆矩阵。为了让结果更有说服力我还添加了5折交叉验证来评估模型稳定性。在程序运行过程中这些信息会以文本和图表两种方式呈现在GUI界面上。GUI交互环节我使用MATLAB的App Designer来搭建界面。界面上半部分是信号展示区可以显示原始波形和CWT时频图下半部分是参数设置区包括小波类型、尺度数、树的数量、特征提取方式等右侧是结果区显示分类准确率、混淆矩阵、特征重要性柱状图。用户加载一段信号后点击“开始诊断”按钮系统会自动完成特征提取和分类并把结果显示在界面上。3. 数据准备与预处理细节3.1 数据集的获取与格式说明这里我用的是凯斯西储大学轴承数据中心公开数据集来做演示这是故障诊断领域最常用的基准数据集网上能方便地下载到。数据包含正常状态和三种故障状态内圈故障、外圈故障、滚动体故障每种状态又有不同故障直径0.007英寸、0.014英寸、0.021英寸和不同负载工况。我这篇文章里为了演示方便只用了0.014英寸故障直径、0马力负载下的数据总共四类正常、内圈、外圈、滚动体。下载下来的数据是MAT文件格式每个文件里有一个变量比如inner_race_0.014_0.mat里面存的是采集到的振动信号。信号采样频率是12kHz每个样本文件包含约12万个采样点持续10秒左右。如果直接把整段信号输入模型计算量太大且信息冗余通常的做法是滑窗切段每个样本取1024个点或者2048个点这个长度既能包含足够多的故障冲击周期又不至于太短导致时频图分辨率不足。文件命名和数据读取有个坑要提醒一下。凯斯西储数据集的文件夹结构比较乱注意看清楚注释不要张冠李戴。我在程序里用dir函数遍历文件夹根据文件名关键字自动识别故障类别比如包含IR的归为内圈故障包含OR的归为外圈故障包含B的归为滚动体故障包含Normal的归为正常。这样代码的通用性会更强换一个数据集只需要修改关键字匹配规则。3.2 信号预处理去均值、去趋势与滤波原始振动信号通常带有直流偏置和趋势项这些低频成分会影响CWT的时频图质量尤其是低频段容易被趋势项污染所以我第一步做的是去均值和去趋势。去均值很简单直接用detrend函数处理。detrend默认是去除线性趋势如果信号有明显的非线性趋势可以指定多项式阶数。对轴承振动信号来说线性detrend基本够用。处理完趋势之后我习惯把信号幅度归一化到[-1, 1]之间这纯粹是为了后续画时频图时颜色映射统一对RF分类本身倒没有太大影响。滤波环节要具体问题具体分析。轴承故障频率通常在几百赫兹到几千赫兹范围但故障冲击会激起结构高频共振所以实际中既需要高通滤波去掉低频干扰又需要低通滤波去掉高频噪声。我写了一个自动带通滤波配置默认用5阶Butterworth滤波器通带设为1kHz到10kHz这段范围基本覆盖了轴承故障诊断的敏感频带。这个参数在GUI中也可以手动调整因为不同设备的共振频带差异很大没有一组参数能通吃所有场景。3.3 训练样本构造与类别平衡处理滑窗切段时要注意重叠率的选择。如果重叠率过高相邻样本之间高度相关会造成训练集和测试集之间的信息泄露模型评估结果虚高。我实测下来在1024个点的窗口下用50%的重叠率切段训练集和测试集之间的相关系数就会低不少评估结果也更可信。完全无重叠切段虽然最严格但样本量太少对RF这种需要大样本的模型不利所以50%重叠是我权衡之后的选择。类别平衡问题也不容忽视。如果四类状态的样本数量差异太大RF会偏向样本量多的类别。处理方法有两种一是过采样对样本量少的类别重复采样二是欠采样丢弃部分样本量多的类别数据。我代码里实现的是过采样方式在小类别样本上重复采样直到与大类别数量一致这样不会浪费数据训练效果也不错。还有一个更高级的替代方案是用SMOTE算法合成新样本但合成数据在故障诊断中需要谨慎因为合成样本可能不符合信号本身的物理规律反而误导模型。4. CWT特征提取的实现过程4.1 MATLAB中的CWT函数选择与参数设置MATLAB从R2016b版本开始cwt函数的用法有较大变化。新版的cwt函数可以直接用cwt(x, fs)这种简洁的调用方式其中x是信号fs是采样频率函数会自动选择合适的morlet小波和尺度范围返回小波系数矩阵和频率向量。对于不需要深究细节的用户来说新版cwt已经足够方便CtrlF在帮助文档里搜索一下就有很多范例。但如果你希望完全掌控参数官方推荐用cwtft函数这是基于傅里叶变换实现的小波变换可以自定义更多细节。不过在我的实践中标准cwt函数配合以下设置就已经能满足要求小波用cmor3-3中心频率3Hz带宽参数3尺度数量设为128对应128个频带扩展方式用默认值。尺度范围的选取需要参考信号的主要频率成分如果采样率12kHz最高分析频率是6kHz奈奎斯特频率那么尺度范围应覆盖这个频率范围对应的尺度值。还有一个细节是VoicesPerOctave参数这是控制每个倍频程内尺度密度的参数。默认值是10提高这个值能让时频图在频率方向更细腻但计算量会增加。对我来说128个尺度已经是兼顾速度和分辨率的折中值测试时也证实了继续增大尺度数对分类准确率提升非常有限仅在绘制高精度时频图时有意义。4.2 时频图的生成与可视化调用cwt得到系数矩阵coefs之后画时频图有几种方式。最简单的是直接使用imagesc函数把系数矩阵的模值画成二维图像横轴是时间纵轴是尺度颜色代表能量大小。但这里有个关键问题尺度到频率的映射不是线性的直接用尺度轴作为纵轴看图的人很难直观理解频率含义。所以我在绘制前会把尺度换算成频率用helperCWTTimeFreqPlot不是内置的而是自己写一个坐标变换函数根据小波的伪频率公式把每个尺度对应的中心频率算出来再用axis xy和ylabel设置成频率轴。颜色映射的选择也有讲究。默认的parula色图在时频图上的对比度不太够我更习惯用jet色图来突出高能量区域因为故障冲击对应的亮斑在jet下非常醒目。当然如果觉得jet色图过于花哨也可以用turbo色图它和jet类似但颜色递变更平滑。实际操作中图像保存的分辨率建议设置到150dpi以上因为后续如果要把时频图作为RF的特征图低分辨率图像会丢失很多细节信息。4.3 时频图特征提取策略时频图生成之后紧接着的问题是如何从这些图像中提取出适合RF输入的数值特征向量。这里有几种可行的策略我按照实际使用效果从高到低列出来区域统计特征是我推荐的首选策略。把时频图在时间方向和频率方向分别等分为N等份和M等份形成N乘M个子块。对每个子块计算三个统计量均值、方差、最大绝对值。这样如果按8乘8划分每个时频图就能得到64乘3等于192维特征。这个策略的物理意义很直观故障特征出现的时频位置不同各子块的能量分布自然不同RF就能依据这些分布差异做出分类判断。奇异值分解特征也是我常用的方法。对时频图矩阵做SVD分解得到若干个奇异值奇异值从大到小排列前20个奇异值基本上就能捕获时频图的主要能量分布特征。这个方法的优势是特征维度低、计算快适合做实时诊断。缺点是奇异值不保留空间结构信息对故障特征的空间分布区分能力有限。边际谱特征走的是另一条思路对时频图沿时间轴求平均得到一条“频带能量分布曲线”然后对曲线采样作为特征。这个特征本质上是傅里叶频谱的广义版本但比FFT更鲁棒因为CWT本身就对非平稳信号更友好。实际使用中边际谱特征和区域统计特征组合使用准确率能再提升一到两个百分点我代码里也提供了组合特征的选项。这里着重推荐一下在GUI中我把区域统计加奇异值分解作为默认特征组合方案两者拼接之后特征维度大约在212维左右。这个维度的特征向量对RF来说是非常友好的既包含了局部时频结构信息又包含了全局能量分布信息训练速度快分类效果也好。4.4 特征归一化是否需要随机森林是树模型不依赖特征之间的距离计算所以严格来说不需要做归一化。这个和SVM、神经网络很不一样那两类模型如果特征尺度差异过大模型训练会出问题甚至不收敛。但在实际操作中有个小问题如果特征中出现极端异常值比如某个子块因为时频图边缘产生了极大的系数值树模型的分裂点选择还是会受到一定影响。所以我程序里还是保留了一个预处理选项——用z-score标准化或者min-max归一化把特征压缩到合理范围默认关闭。读者可以根据自己的数据情况决定是否打开。如果特征分布在多个数量级之间横跳建议打开归一化如果特征分布本身就比较均匀保持原始特征其实更好因为树模型的分裂阈值更容易解释。5. 随机森林模型的构建与优化5.1 建立RF模型TreeBagger还是fitcensembleMATLAB中实现随机森林主要有两个接口TreeBagger和fitcensemble。两者底层都是基于决策树集成的思想但接口和功能细节有区别刚开始接触很容易纠结到底该用哪个。TreeBagger是更经典的随机森林接口很多老代码和教程都用它。优点是可以直接输出OOB误差、特征重要性、树的相似度矩阵等信息这些是分析模型内部逻辑的利器。缺点是接口风格相对古老参数设置用名值对方式初学者可能需要多看几遍文档才能搞明白。fitcensemble是更现代的集成学习框架不仅能做随机森林还能做AdaBoost、Bagging、GentleBoost等其他集成方法。它使用起来更简洁配合predict函数做预测很顺手还能自动进行交叉验证。缺点是获取特征重要性时需要通过predictorImportance函数单独调用OOB误差也不能直接得到。我自己在这个项目里选择了TreeBagger主要原因是可以直接在GUI中画出OOB误差曲线并且在训练过程中动态展示误差随树数量增加的变化趋势这对教学演示和调试非常有帮助。如果你更看重代码简洁性用fitcensemble也没问题两者的分类准确率几乎没有显著差异。在正文代码里我以TreeBagger为主线来讲解。5.2 关键参数设置与交叉验证TreeBagger的关键参数首先是NumTrees也就是树的数量。我实测了50到1000棵树的表现50棵时准确率约95%200棵时达到98%左右500棵以上基本稳定不再提升。综合计算效率和准确率我默认设置为200棵。在GUI中这个参数可以手动调整如果数据集特别大可以把树的数量降低到100棵准确率损失通常在1%以内但训练时间能缩短一半。其次是MinLeafSize这是控制每棵树叶子节点最小样本数的参数。值越小单棵树越复杂容易过拟合值越大树越简单可能欠拟合。对故障诊断这种噪声较多的信号数据我推荐MinLeafSize设置在5到10之间默认值我选了5这样单棵树足够敏感又不会因为过度拟合训练集的噪声而损害泛化能力。第三个关键参数是NumPredictorsToSample即每次分裂时随机选择的特征数量。默认值是特征总数的平方根这是Breiman在原始论文中的建议实测表现也确实不错。不过当特征总数特别大时可以考虑适当调大这个值比如从默认的15提高到20每次分裂能看到更多候选特征树之间的相关性会降低总体准确率有小幅提升。模型训练完成之后用kfoldLoss或者crossval做5折交叉验证是很好的校验手段。我在程序里写了一个辅助函数返回交叉验证的准确率和标准差如果准确率标准差超过2%说明数据切分或者模型稳定性有问题需要检查一下样本是否存在明显的顺序相关性。5.3 模型评估指标体系准确率是最直观的指标但不是唯一的指标在故障诊断中尤其要关注混淆矩阵和F1分数。比如内圈故障和外圈故障在某些工况下信号特征非常接近如果只看总准确率你可能意识不到模型在这两类之间频繁混淆。而混淆矩阵能一眼看出模型具体在哪些类别上犯了错这对于后续优化特征提取方式或者增加样本量有重要的指导意义。我还习惯计算每个类别的精确率Precision和召回率Recall。精确率衡量的是模型判定为某类故障的样本中真正属于该类别的比例召回率衡量的是某类故障的真实样本中模型成功找出来的比例。在工业场景中漏报的代价往往比误报更高——一个真实的故障被漏检可能导致设备损坏而一个误报顶多多花一次停机检查的时间所以可能需要适当关注召回率指标。如果发现有某类故障的召回率偏低可以针对性增加该类别的训练样本量或者调整特征提取中能突出该故障特征的频带参数。5.4 特征重要性分析与模型解释TreeBagger的OOBPermutedPredictorDeltaError属性可以输出特征重要性得分这是随机森林一个非常实用的优势。原理是这样的训练完成后对每个特征列的值做随机打乱然后重新计算OOB误差误差上升越多说明该特征对预测越重要。这个指标用两类信息一是分类准确率的变化量二是margin的变化量后者更灵敏。我把特征重要性可视化放在了GUI的结果面板中。有一次我在调试齿轮箱故障时发现排在前几位的特征全部落在500Hz到2kHz频带对应的子块上这个信息直接帮我确认了故障特征集中在低频段于是我把带通滤波的下限从2kHz调整到了500Hz整体分类准确率提升了3%。这就是特征重要性带来的实际价值——它不只是模型内部的一个指标反过来还能指导信号处理和特征工程的设计。6. GUI设计与功能实现6.1 界面总体布局与设计原则用MATLAB的App Designer来设计GUI比起传统的GUIDE设计体验和代码可维护性都好得多。App Designer基于面向对象的编程模型组件回调函数自动生成框架你只需要填充逻辑代码界面和代码分离得比较干净后期扩展新功能时不容易破坏已有功能。界面布局上我遵循“左中右”三栏设计。左侧是数据加载和参数设置区域中间是信号展示和时频图展示区域右侧是模型训练和结果展示区域。整个窗口大小设置为1280乘800像素这个尺寸在常见的1080P显示器上能完整展示不用滚动滚动条如果屏幕分辨率较低App Designer的自动缩放功能也能适配。在设计GUI时有一个容易被忽视的点交互流畅性。CWT计算和RF训练都是耗时操作如果在回调函数里直接同步执行界面会卡死几秒甚至十几秒用户会以为程序崩溃了。所以我在耗时操作中插入了uiprogressdlg进度条组件让用户知道程序在正常运转。更复杂的做法是使用parfeval异步执行后台任务但会增加代码复杂度我在文章里不展开讲仅在进度条方案不能满足需求时提醒读者考虑。6.2 各功能模块详解数据加载模块是整个项目的入口。我用的是uigetfile函数弹出文件选择对话框支持.mat格式文件。加载完成后系统会自动读取信号数据和采样频率信息在界面的“信号预览”坐标轴中用plot函数画出原始波形并用text函数显示信号长度、采样频率、最大幅值等基本信息方便用户确认数据是否正确。参数设置模块涵盖了CWT和RF两大部分。CWT的控件包括小波类型下拉框morlet、mexh、db4等选项、尺度数滑块、频率显示范围编辑框。RF的控件包括树的数量编辑框、最小叶子节点大小编辑框、特征选择策略下拉框区域统计、SVD、区域统计加SVD组合。所有参数都有默认值用户只需按需微调不需要每一个都设置这对新手特别友好。训练与诊断模块是三段式设计。第一步点击“特征提取”按钮系统遍历所有信号样本计算CWT时频图和特征向量支持在进度条上看到每个类别的处理进度。第二步点击“训练模型”按钮系统用选定的特征矩阵训练树模型训练完成后画出OOB误差曲线。第三步点击“模型评估”或“开始诊断”按钮系统在测试集上评估模型性能显示混淆矩阵和分类报告。结果显示模块中包括多个坐标轴和表格组件。准确率用大号文本显示显示在界面的显眼位置包括训练集准确率和测试集准确率。混淆矩阵用heatmap函数画出来横纵坐标都是类别名称颜色深浅代表样本数量。特征重要性用柱状图显示。如果用户加载的是单个测试样本界面会显示样本的时频图以及模型输出的分类概率向量这个功能在演示和调试时特别直观。6.3 回调函数与数据传递设计App Designer的回调函数设计有几个常见问题需要提前预防。第一个是共享数据的传递问题。MATLAB App Designer中组件的数据存储在app对象的属性中你在回调函数里用app.data这样的自定义属性来存储中间变量在别的回调函数中直接读取app.data即可。官方推荐的startupFcn也可以做初始化比如设置默认参数、加载预训练模型等。第二个问题是按钮的状态管理。如果在训练还在进行的过程中用户又点击了“训练模型”按钮程序会重复执行训练并可能出错。我在训练开始时用app.TrainButton.Enable off禁用按钮训练结束后重新启用。同时如果在训练过程中用户修改了参数系统会弹窗提示“参数已修改需要重新训练”避免用户误以为当前模型是基于最新参数训练的。第三个问题是MATLAB版本的兼容性。App Designer在R2016a之前的版本里并不存在如果你还在用老版本的MATLAB需要用传统的GUIDE环境来搭建界面。不过目前主流的R2020b之后的版本都内置了App Designer直接使用即可代码中的cwt、TreeBagger这些核心函数在R2018b之后都没有大的接口变化所以兼容性问题只存在于极老的版本中。6.4 打包与发布GUI开发完成后可以打包成独立EXE程序发给其他人使用即使对方电脑上没有安装MATLAB也能运行前提是需要安装MATLAB Compiler RuntimeMCR。在MATLAB的命令窗口中执行deploytool选择“Application Compiler”添加主文件和所有依赖文件设置应用名称和图标点击打包就能生成可执行文件。需要注意的一点是打包后的程序运行速度比MATLAB原生环境慢一些尤其是在CWT计算部分因为MCR没有JIT加速。我做过的实测中同样是128个尺度的CWT计算打包后在启动时可能要额外花1到2秒的时间初始化运行环境这是正常现象不需要担心。如果读者需要将程序分发给同事建议在说明文档中写清楚MCR版本的对应关系避免因为版本不匹配导致程序无法启动。7. 完整代码详解与核心模块解析7.1 主程序框架与文件结构整个项目由几个核心函数组成结构如下所示RootDirectory / ├── main_diagnosis.m % 主程序入口负责调用各模块 ├── load_signals.m % 数据加载函数 ├── cwt_feature_extract.m % CWT特征提取函数 ├── train_rf_model.m % RF模型训练与评估函数 ├── app_CWT_RF_diagnosis.mlapp % GUI应用程序 └── utils / ├── scalogram_to_features.m % 时频图到特征向量的转换 └── plot_confusion_matrix.m % 混淆矩阵可视化函数主程序main_diagnosis.m的作用是串起整个流程调用load_signals读取原始数据调用cwt_feature_extract提取特征调用train_rf_model完成分类最后把结果保存到workspace或者导出为文件。如果用户不想用GUI可以直接运行主程序本质上GUI就是对这个流程的包装。程序运行结束后在命令窗口打印一份包含准确率、F1分数、特征维数、训练耗时的摘要信息并自动保存特征矩阵和模型文件方便后续研究分析。7.2 数据加载与预处理代码解析load_signals.m函数的核心部分如下面的代码所示。代码重点关注如何根据文件名自动识别故障类别这是在数据处理过程中的核心步骤。function [signals, labels] load_signals(data_dir, fs, window_len, overlap) % 参数说明: % data_dir: 数据集文件夹路径 % fs: 采样频率(Hz) % window_len: 滑窗长度采样点数 % overlap: 重叠率(0~1) files dir(fullfile(data_dir, *.mat)); signals []; labels []; category_map {Normal, IR, OR, B}; cat_labels [1, 2, 3, 4]; for i 1:length(files) fname files(i).name; % 判断类别 cat 0; for c 1:length(category_map) if contains(fname, category_map{c}) cat cat_labels(c); break; end end if cat 0 warning(跳过未识别文件: %s, fname); continue; end data load(fullfile(data_dir, fname)); fields fieldnames(data); sig data.(fields{1}); % 去除趋势项和直流偏置 sig detrend(sig, linear); % 滑窗切段 step round(window_len*(1-overlap)); n_seg floor((length(sig)-window_len)/step) 1; for j 1:n_seg idx (j-1)*step1 : (j-1)*stepwindow_len; seg sig(idx); % 简单幅值归一化 seg seg / max(abs(seg)); signals [signals; seg]; labels [labels; cat]; end end end这里的detrend函数会把信号的直流偏置和线性趋势去掉让CWT分析不受低频漂移影响。滑动窗口的步长计算公式很容易出错注意最后一段如果剩余长度不足窗口长度就会丢弃所以n_seg的公式里用了floor取整。当文件很多时循环拼接的效率会比较低但考虑到轴承数据的样本量通常不大这种方式简洁直观对性能影响可以接受。7.3 CWT特征提取核心代码解析cwt_feature_extract.m和scalogram_to_features.m是整套方法的核心。第一个函数从信号数组批量计算时频图并把每张时频图传递给第二个函数提取特征向量。function featMat cwt_feature_extract(signals, fs, wavelet, numScale, blockN, useSVD) n size(signals, 1); featMat zeros(n, 0); % 初始为空每次循环按特征维度拼接 for i 1:n sig signals(i, :); [coefs, freqs] cwt(sig, amor, fs); % 使用复数morlet小波 % 注意: amor会在新版MATLAB自动选择尺度数但如果要手动控制尺度 % 可以改用 cwt(sig, 1:numScale, morl) 这种形式 % coefs: numScale x window_len 的复数矩阵 mag abs(coefs); % 裁剪到目标尺度数 numScaleCurr size(mag, 1); if numScaleCurr numScale mag mag(1:numScale, :); else % 如果尺度数不够在低频侧复制填充简单策略 padNum numScale - numScaleCurr; mag [mag; repmat(mag(end, :), padNum, 1)]; end % 提取区域统计特征 feat scalogram_to_features(mag, blockN); if useSVD [U, S, V] svd(mag, econ); sv diag(S); sv sv(1:min(20, length(sv))); feat [feat, sv]; end featMat(i, :) feat; end end有一个重要经验要分享cwt函数在处理短信号时输出的尺度数可能和预期的不同。例如设置尺度数为128但如果信号本身太短算法会自动减少尺度以避免边界效应。所以我在代码里做了尺度数检查和填充确保输出矩阵的行数一致。在处理不同长度的信号时这个检查就显得尤为重要而同一批数据中有些样本信号较短的情况在实际工程中并不少见。scalogram_to_features函数的实现是把mag矩阵均匀划分成blockN乘blockN个子块每个子块内计算均值、方差和最大值三项统计量。最终特征维数为3乘blockN乘blockN。在默认设置下blockN等于8特征维度为192。7.4 随机森林训练与评估完整代码解析train_rf_model.m函数实现了RF的训练、评估和可视化。以下是完整代码加注释function [mdl, accuracy, confMat, imp] train_rf_model(featTrain, labelTrain, featTest, labelTest, numTrees, minLeaf) % 训练Random Forest模型 % TreeBagger返回的对象后续可以用predict函数做预测 mdl TreeBagger(numTrees, featTrain, labelTrain, ... Method, classification, ... MinLeafSize, minLeaf, ... OOBPrediction, on, ... OOBPredictorImportance, on, ... NumPredictorsToSample, all); % 绘制OOB误差曲线 figure(Name, OOB Error Curve, Color, w); plot(oobError(mdl), LineWidth, 1.5); xlabel(Number of Grown Trees); ylabel(Out-of-Bag Classification Error); grid on; % 在训练集和测试集上评估 [predTrain, ~] predict(mdl, featTrain); [predTest, scoresTest] predict(mdl, featTest); predTrain str2double(predTrain); predTest str2double(predTest); % 计算准确率 accTrain sum(predTrain labelTrain) / length(labelTrain) * 100; accTest sum(predTest labelTest) / length(labelTest) * 100; accuracy [accTrain, accTest]; % 混淆矩阵 confMat confusionmat(labelTest, predTest); % 特征重要性 imp mdl.OOBPermutedPredictorDeltaError; % 显示结果 fprintf(训练集准确率: %.2f%%\n, accTrain); fprintf(测试集准确率: %.2f%%\n, accTest); end注意一个细节TreeBagger返回的预测结果是字符串形式的类别标签需要转换回数值这里用str2double做了转换。在GUI中做预测时同样要注意这个细节用str2double转换predict函数返回的字符串。NumPredictorsToSample参数设置为all时等价于Breiman提出的“全部特征参与分裂”实际效果多数情况下和sqrt方法差异不大。如果特征数特别高建议把all改成sqrt以加速训练。7.5 单样本实时诊断函数设计除了批量评估我还实现了一个单样本诊断功能用于GUI中的实时诊断场景。diagnose_single_sample函数接收一段信号、CWT参数、RF模型和特征提取参数输出分类结果和各类别概率。function [label, score] diagnose_single_sample(sig, fs, mdl, blockN, useSVD) % 信号预处理 sig detrend(sig, linear); sig sig / max(abs(sig)); % CWT [coefs, ~] cwt(sig, amor, fs); mag abs(coefs(1:128, :)); % 特征提取 feat scalogram_to_features(mag, blockN); if useSVD [~, S, ~] svd(mag, econ); sv diag(S); sv sv(1:min(20, length(sv))); feat [feat, sv]; end % 预测 [labelChar, score] predict(mdl, feat); label str2double(labelChar); end这个函数把所有步骤封装在一起GUI中调用只需要一行代码很适合作为对外接口。如果后续想把系统部署到实时采集系统中只要把输入参数从信号数组改成实时数据流其他逻辑可以重用。8. 常见问题与排查技巧实录8.1 CWT相关时频图模糊、边界效应、计算卡顿时频图模糊是CWT分析中非常常见的问题。可能原因一是尺度数设置太少导致频率分辨率不足解决方法很直接把尺度数从64提高到128或256就能看到明显改善原因二是信号本身太长在GUI中显示时做了降采样或者压缩显示仔细检查一下坐标轴显示范围是不是覆盖了整段信号。还有一种情况是选错了小波基比如用mexh小波处理冲击性故障信号时最好先对比测试几种小波的效果再定。边界效应是CWT的固有问题信号两端的CWT系数会有失真这是小波中心在信号边界附近缺少足够数据支撑导致的。解决方法是扩展信号边界在CWT计算前用对称延拓或周期延拓把信号向两端各扩展一定长度计算完后剪掉边界区域。MATLAB的cwt函数内部已经做了边界处理但处理后的时频图两端仍可能有伪影在做特征提取时我建议忽略时频图两端约5%宽度的区域只统计中间90%的部分这样特征更稳定。计算卡顿问题的根源是CWT的时间复杂度与尺度数、信号长度的乘积成正比。在GUI中拖动滑块时每次都会重新计算CWT避免卡顿的方法是采用“松手后才重算”策略。在App Designer中把滑块的ValueChangingFcn里不触发重绘而在ValueChangedFcn中触发或者在滑块回调中设置一个防抖定时器只有当用户停止操作一段时间后才执行重算。8.2 RF相关过拟合、特征重要性异常、预测结果为NaN过拟合的表现是训练准确率接近100%但测试准确率明显偏低。这种情况多出现在特征维度较高而训练样本不足的情况下。解决方法有几种一是增加训练样本比如把滑窗重叠率从50%提高到75%二是减少特征维度比如把子块划分从8乘8改成6乘6三是提高MinLeafSize从5提高到10让单棵树更简单。特征重要性异常比如所有特征的得分都很接近且很小说明模型的特征选择能力没有发挥出来。可能的原因是特征之间存在强相关性RF在分裂时只能随机选到其中一个相关信息导致每棵树都“不知道”特征的完整价值。这种情况不一定影响预测精度但对特征解释不利。更好的做法是先用相关性分析和PCA预处理剔除高度相关的冗余特征再重新训练。预测结果为NaN的问题我遇到过几次原因通常是训练集中存在缺失值或者无穷大值而RF不能处理NaN输入。排查方法很简单训练前用any(isnan(feat))和any(isinf(feat))对特征矩阵检查一遍。如果发现问题定位到具体的样本和特征列检查这个样本是不是信号长度不足导致的特征矩阵填充异常。另外要注意在GUI中如果加载了未经过预处理的原始数据或者特征提取函数里出现了除零操作也会产生NaN或Inf这就要去检查特征提取函数里的分母项。8.3 GUI开发相关控件不响应、打包失败控件不响应最常见的原因是回调函数中出现了死循环或长时间同步操作事件队列被阻塞界面就会“假死”。在App Designer中解决方法是把耗时操作分成小段执行每次处理几十个样本就调用drawnow刷新界面或者使用定时器定时执行不可分割的步骤。如果不需要动态刷新界面至少也要在耗时操作前后显示和关闭进度条让用户知道程序没有卡死。打包失败的情况也遇到过不少。最常见的原因是主程序中引用了不在MATLAB路径下的自定义函数打包时这些依赖没有被自动找出。解决方法是在打包前用MATLAB的“依赖关系分析器”工具对所有引用进行一一核对。还有一次失败是因为用了带有中文注释的函数名称打包后路径解析出错后来全部改成英文字母就解决了。另外需要注意的是打包后首次启动速度会比较慢因为MCR需要初始化运行时环境这是正常现象。8.4 模型精度提升一种快速有效的调参思路当模型精度不满足要求时我建议按照以下顺序逐步调试这样效率最高。第一步检查数据本身的问题。确认信号是否包含了足够明显的故障特征如果故障特征不明显可能需要先做针对性滤波或者解调分析。第二步调CWT参数。尝试不同的小波基和尺度数观察时频图上故障特征是否清晰这一步的调优效果往往比RF参数调整更明显。第三步加特征。把单一的特征提取方式改为组合特征比如区域统计加SVD加边际谱合理组合不同维度的信息。第四步调RF参数。交叉验证配合调优GridSearch遍历不同的树数量和MinLeafSize组合找到当前数据下的最佳参数。最后一步如果仍然不满足要求考虑使用更复杂的分类器或者对原始信号做进一步的预处理比如经验模态分解EMD滤波等。在实际项目中通过这个调参流程我成功把某风机齿轮箱故障数据集上的准确率从最初的92%提升到了98.7%。关键提升步骤是对信号做了1kHz高通滤波并且把特征子块划分从均匀划分改成了低频加密划分。这些改进都是通过观察时频图和特征重要性分布得到的灵感这再次说明了模型解释工具在工程实践中的重要性。9. 扩展思路与后续改进方向如果读者想把这套CWT-RF故障诊断系统应用到更广泛的场景中我有几个建议可以分享。数据增强层面除了滑窗重叠采样还可以通过添加不同信噪比的噪声来增加训练样本的多样性让模型对噪声更鲁棒。比如原始信号中加入5dB、10dB、15dB的白噪声生成多组增强样本模型的泛化能力通常有2到3个百分点的提升。这种方法在测试集来源于实际工况而非理想实验室条件时效果尤其明显。特征层面可以把CWT和深度学习方法结合起来。比如用CWT生成时频图再输入到预训练的CNN网络做特征提取最后接入RF分类器形成“CWTCNN特征RF”的混合结构。我在一个实验中发现这种混合结构比单独使用CWT-RF或单独使用CNN都更稳定特别是在小样本场景下RF作为最终分类器能有效缓解CNN的过拟合问题。实时诊断层面当前程序主要处理的是离线数据如果要做在线实时监控需要考虑两点优化一是采用滑动窗口滚动计算CWT每次只计算最新一段数据的时频图二是预先固化模型参数把训练完的TreeBagger模型保存成.mat文件在线诊断时直接load进来避免重复训练。经过两步优化后单次诊断的耗时可以压缩到几十毫秒量级满足大多数工业现场的实时性要求。我自己在实际部署中还发现MATLAB的并行计算工具箱能进一步加速多通道信号的特征提取如果现场有多个传感器同时采集信号这个优化方向的收益非常可观。再说一个非常重要的工程实践建议在将模型部署到现场之前务必在不同的工况条件下采集足够多的验证数据。模型在实验室数据上表现良好到了现场工况变化之后准确率下降十分常见这不是算法本身的问题而是数据的分布变了。解决方法是建立包含多工况的数据集所谓“数据驱动数据为王”对于故障诊断任务来说尤其如此。最后想说的是CWT加RF这个组合的价值不仅在于它的准确率更在于它的可解释性。当你面对一套数据时时频图可以告诉故障在时间和频率上的分布规律特征重要性可以告诉模型依据哪些特征做决策这种透明性在工程诊断中是至关重要的。黑箱模型再准如果不能定位到具体的故障模式和特征频率也很难让现场工程师信服。这也是我在做故障诊断时始终把模型的可解释性放在首位的原因。
返回列表