搞过旋转机械故障诊断的朋友应该都有体会,最难的不是把某个单一算法跑通,而是把“信号进来—特征出来—故障判出来—结果让人信服”这条链路完整串起来。这套项目标题里写的ICEEMDAN-CMBE特征提取、Bayes-TCN-BiLSTM故障诊断、SHAP可解释分析,正好覆盖了这条链路里的全部关键节点,而且全部基于MATLAB实现。我把它拆开复现过两遍,今天把每个模块背后的原理、参数取舍、踩坑记录一次性讲清楚,适合正在做轴承或齿轮箱故障诊断课题、想用深度学习替代传统浅层模型、或者已经跑通基础CNN但想知道怎么加可解释性分析的朋友参考。
这套东西本质上解决的是三个问题:第一,原始振动信号是非线性非平稳的,直接扔进神经网络效果不稳定;第二,模型精度上去了但不知道怎么证明它学到了物理特征而不是噪声;第三,超参数调起来靠感觉,复现性差。ICEEMDAN负责把信号分解成稳定分量,CMBE负责从分量里量化特征,TCN-BiLSTM负责时序建模,贝叶斯优化负责自动找超参,SHAP负责把决策过程摊开给人看。下面按模块逐个拆。
1. 这套技术栈到底在解决什么问题
1.1 传统故障诊断的三大痛点
我最早做故障诊断用的是“手工特征+随机森林”这套老组合,当时最头疼的就是特征不知道该提多少、提哪些。振动信号里包含正常运转的周期分量、轴承故障冲击分量、齿轮啮合分量,还有大量环境噪声,时域和频域特征互相重叠。后来开始用深度模型直接端到端学习,精度确实上去了,但新问题又冒出来:网络本质是个黑盒,出了结果得靠“实验经验”担保,写论文和做汇报时总不能只说“模型认为这个是内圈故障”。这就是第一个痛点——模型不可解释。
第二个痛点是特征质量配不上模型复杂度。很多端到端方案把原始信号直接切成固定长度喂给卷积网络,稍微换一个工况、换一个转速,精度就掉得很厉害。因为原始信号里带载波、带噪声、带转频调制,模型需要花费大量参数去对抗这些干扰。如果先用自适应分解把信号拆成物理含义更明确的IMF分量,再用熵类指标做量化,输入网络的特征维度可以大幅压缩,泛化能力反而更好。
第三个痛点是超参数调优。TCN的膨胀系数、卷积核数量、BiLSTM隐层节点数、学习率、Dropout比例,这些变量互相影响,手动调一遍至少要跑几十组实验。用贝叶斯优化之后,虽然每次评估仍然要训练一次网络,但整体实验次数能压到二三十次以内,并且能找到比手动调参更稳定的组合。
1.2 为什么把五种方法拼在一起
看到这个标题的时候,大概率会有人问:ICEEMDAN都做完分解了,直接对IMF提取时域特征不就行了吗,为什么还要算CMBE?这里有个关键逻辑:轴承故障冲击会改变信号在不同尺度上的复杂度分布,单看IMF的方差、峭度这些指标只能反映整体强度变化,但不同故障类型的复杂度模式差异需要用多尺度熵来刻画。CMBE就是把这个“复杂度指纹”提取出来,输入给分类器。
那有了CMBE特征,为什么还要TCN-BiLSTM?因为特征是按时间窗滑动提取的,每个时间窗得到的熵值序列本身就包含故障演化的时间趋势。TCN擅长捕捉局部时序依赖,BiLSTM擅长捕捉长距离双向依赖,两者串联之后,相当于先做局部特征提炼再做全局时序推理。SHAP放在最后一步则服务于“信任”问题:它告诉工程师每一个输入特征对最终故障类别判断贡献了多少,让模型的决策边界可以被审查。
这五个模块的关系用一句话概括:分解是为了净化输入,熵特征是为了压缩表达,时序网络是为了捕捉演化规律,贝叶斯优化是为了公平调参,SHAP是为了让结论站得住脚。下面按执行顺序逐个展开。
2. ICEEMDAN信号分解:把非线性非平稳信号“拆开”
2.1 ICEEMDAN与CEEMDAN/EMD的关键差异
经典EMD最大的问题是模态混叠,就是一次冲击事件被拆进好几个IMF里。为了解决这个问题出现了EEMD,靠添加白噪声来平均掉模式混叠,但EEMD分解不彻底,且每次都要对全信号做多次EMD,计算开销非常大。CEEMDAN虽然改进了噪声添加方式,但在早期迭代中仍然会在IMF里残留部分噪声。ICEEMDAN的改进思路更干脆:不再直接把高斯白噪声加给信号,而是加“经过EMD分解后的噪声分量”,也就是用E_k(w(i))来构造参考信号,这样可以保证分解得到的残差信号中噪声能量被大幅压低,IMF的完备性也更好。
对于轴承振动这种冲击性强、信噪比波动大的信号,ICEEMDAN的优势非常直观:分解出的第一个IMF通常直接对应最强的故障冲击分量,后续IMF才逐步呈现调制谱和低频转频分量。我在实际对比中发现,同样一段内圈故障信号,CEEMDAN可能需要十几个IMF才能把故障特征和噪声分开,而ICEEMDAN通常前五六个IMF就能把主要物理成分拎出来。
2.2 MATLAB中的参数设置与实现要点
MATLAB代码包里ICEEMDAN函数的核心参数有三个:噪声标准差Nstd、集合次数NE、最大迭代数MaxIter。我把常用范围列在下面。
| 参数 | 建议范围 | 说明 |
|---|---|---|
| Nstd | 0.05~0.2 | 越大抗噪能力越强,但过大可能会破坏信号弱冲击成分 |
| NE | 50~200 | 集合次数越多分解越稳定,但计算时间线性增长 |
| MaxIter | 500~1000 | 筛分迭代次数,轴承信号一般500足够 |
| SNRFlag | 1 | 自适应降噪标志位,通常开启 |
调用格式一般是:
imfs = ICEEMDAN(signal, Nstd, NE, MaxIter, SNRFlag);初次复现时建议先对原始信号做一次零均值标准化,再设置Nstd=0.1、NE=100。标准化这段经常被忽略,但如果不做,能量大的通道会让熵值计算失去相对可比性。分解完之后不要着急把所有IMF都留下,先用相关系数和峭度做一个筛选:相关系数反映IMF与原始信号的整体一致性,峭度反映冲击性,两个指标结合可以淘汰掉纯噪声分量。
2.3 模态混叠与端点效应的实战处理
模态混叠在ICEEMDAN中已经大幅减少,但并不是零。当轴承转速波动大或者存在变载荷工况时,相邻IMF仍然会偶发频率重叠。我的处理习惯是对每个IMF做边际谱分析,画出Hilbert边际谱后检查频率带是否出现明显的交叉重叠。如果出现,就把这两个IMF合并成一个特征通道,不再强制要求物理意义上的完全正交。
端点效应是另一个老问题。信号两端在分解时会出现虚假振荡,最直接的办法是数据延拓。我用过镜像延拓和AR预测延拓两种方式,实测下来镜像延拓对周期性较强的齿轮信号更稳,AR预测延拓对突变型轴承信号更自然。但是要注意,延拓长度一般控制在信号总长的5%以内,延拓太长会让计算量膨胀,并且会引入虚假模式。
这段代码里我用的是镜像延拓函数,在调用ICEEMDAN之前先执行:
extLen = round(length(signal) * 0.05); signal_ext = [flipud(signal(1:extLen)); signal; flipud(signal(end-extLen+1:end))]; imfs_ext = ICEEMDAN(signal_ext, Nstd, NE, MaxIter, SNRFlag); imfs = imfs_ext(extLen+1:end-extLen, :);注意延拓完成后要严格裁剪回原始长度,否则后面滑窗截取特征时会产生索引错位。
3. CMBE复合多尺度熵:量化IMF中的故障特征
3.1 从单尺度熵到复合多尺度熵的演进逻辑
单一尺度上的样本熵或模糊熵只描述信号在一个粗粒化粒度下的复杂度。轴承故障信号的特殊之处在于,不同故障类型对多个时间尺度的影响不同:内圈故障的高频冲击在最小尺度上变化最明显,外圈故障的调制特征在中尺度上更突出,滚动体故障则会在多个尺度上产生间歇性复杂度波动。单尺度熵看到的是“快照”,多尺度熵看到的才是“行为模式”。
传统多尺度熵有一个明显的坑:粗粒化过程中数据长度会被压缩,尺度越大,可用的数据点越少,熵值估计的方差就越大。复合多尺度熵的改进是,在每个尺度上不只用第一个粗粒化序列,而是把所有可能的粗粒化起点都算一遍,然后取平均。这样一来,短数据情况下的估计稳定性大幅提升,对滚动轴承这种一次采集只有几千个点的场景特别友好。
需要说明的是,标题里的CMBE属于复合多尺度熵族,MATLAB实现中我默认使用模糊熵算子,因此有些文献里也直接写CMFE。两者核心粗粒化框架完全一致,区别只在于下层的熵算子。代码里我留了算子切换开关,你如果想用样本熵或散布熵,只要把熵函数替换掉就可以,特征矩阵的拼接逻辑不需要改动。
3.2 CMBE特征矩阵的构建流程
CMBE的计算流程可以拆成四步。第一步,对每个保留的IMF做归一化,避免幅值差异主导熵值。第二步,给定最大尺度因子,比如20,对每个尺度生成对应的复合粗粒化序列:
function seqs = compositeCoarseGrain(x, scale) % x: 输入序列 % scale: 当前尺度因子 len = floor(length(x) / scale); seqs = zeros(scale, len); for start = 1:scale seqs(start, :) = mean(reshape(x(start:start+len*scale-1), scale, len), 1); end end第三步,对每个起点的粗粒化序列都在嵌入维度m和相似容差r下计算一遍模糊熵,然后把所有起点的熵值取平均,得到当前尺度下的复合熵。第四步,把所有尺度下的熵值拼接成一个向量,这个向量就是这个IMF的CMBE特征。
实际操作时,我一般保留前8到10个IMF,每个IMF提取20个尺度的熵值,这样单个时间窗的特征向量维度是160到200维。如果后续网络仍然出现过拟合,优先把尺度因子降到12而不是减少IMF数量,因为IMF降维会丢失物理特征,而降低尺度只是减少冗余表达。
3.3 尺度因子与熵参数如何选
尺度因子、嵌入维度、容差这三个参数的搭配直接影响诊断结果。我给出一个经多次实验验证的起始组合:嵌入维度m=5,模糊幂指数n=2,相似容差r=0.15×std(即信号标准差的0.15倍),最大尺度=20。
m=5是对轴承振动信号的折中取值,太小区分不出不同故障的复杂度差异,太大又会让熵值对噪声过于敏感。r取0.15倍标准差是因为太小容差会让熵值波动剧烈,太大又会把故障和正常状态的差异抹平。这里有一个容易犯的错:如果对每个IMF独立标准化后再计算熵,r就必须按标准化后的std来设置,如果沿用原始幅值下的r,会把不同IMF的量纲混在一起。
尺度因子方面我建议先用20跑通全流程,再测试10和15对比。如果发现精度变化在1%以内,就可以放心用15,计算速度能明显提升,因为尺度因子越大复合粗粒化涉及的重叠计算越密集。
4. Bayes-TCN-BiLSTM:超参自动寻优的时序分类网络
4.1 为什么用贝叶斯优化而不是网格搜索
网格搜索在这个场景下基本不可行。假设要优化6个超参数,每个参数只取5个候选值,组合数就是5的6次方等于15625组实验,每组还要训练一遍完整网络,这对普通工作站来说意味着数天甚至数周的计算量。随机搜索虽然减少了实验次数,但完全忽略历史评估信息,浪费了已跑过实验的反馈价值。
贝叶斯优化的核心是用高斯过程代理模型去拟合“超参数组合→模型性能”的映射关系。每次评估完一组超参后,代理模型会更新,并通过采集函数(比如EI)自动选择既有潜力又有不确定性的点作为下一组候选。简单说,它是在“探索”和“利用”之间做了平衡,实际跑下来的效果是:十几组实验就能把精度推向收敛平台,比手动调参少花一半以上的时间。
MATLAB里直接用bayesopt函数,配合自定义目标函数即可:
vars = [ optimizableVariable('numFilters', [16, 64], 'Type', 'integer') optimizableVariable('numBlocks', [2, 6], 'Type', 'integer') optimizableVariable('lstmUnits', [32, 128], 'Type', 'integer') optimizableVariable('learnRate', [1e-4, 1e-2], 'Transform', 'log') optimizableVariable('miniBatchSize', [16, 64], 'Type', 'integer') optimizableVariable('dropout', [0.1, 0.5], 'Type', 'real') ]; results = bayesopt(@objFcn, vars, ... 'MaxObjectiveEvaluations', 30, ... 'AcquisitionFunctionName', 'expected-improvement-plus', ... 'UseParallel', true);注意目标函数内部必须完成“从超参到网络搭建、训练、验证集评估、输出负准确率”的完整流程,并且每次训练时固定随机种子,否则同一组超参两次评估的结果波动会把优化过程带偏。我踩过这个坑,第一次跑贝叶斯优化时没有固定种子,优化曲线一直在抖动,换了种子后波动小了很多,收敛也快了不少。
4.2 TCN+BiLSTM的结构设计思路
TCN用膨胀因果卷积捕捉局部时序特征。膨胀系数按指数增长(1,2,4,8,...),意味着最后一层卷积的等效感受野能覆盖很长的历史区间,但又不会像普通RNN那样出现梯度消失。在故障诊断场景里,TCN的压缩能力很强,可以把160维的CMBE特征序列先映射成一组局部模式向量。
BiLSTM放在TCN后面,负责对TCN输出的局部模式做双向时序建模。轴承故障特征不是对称的,故障冲击会从出现到加剧再到平稳,前向LSTM能捕捉特征的累积趋势,后向LSTM能捕捉相对于后续状态的相对变化,两者拼接后得到的语义更完整。
一个常用的网络结构是:输入层→TCN残差块×N→BiLSTM层→全连接→Softmax。TCN每个残差块内部包含卷积、批归一化、ReLU和Dropout。需要特别注意的是:BiLSTM的输入维度必须等于TCN输出的通道数,设计时要把这两个数字对上。我在代码里将TCN输出通道设置为64,这样BiLSTM层的输入就是64维。
搭建关键层用MATLAB代码表示大致是:
lgraph = layerGraph(); tempLayers = [ sequenceInputLayer(size(features, 2), 'Name', 'input') convolution1dLayer(3, numFilters, 'Padding', 'causal', 'DilationFactor', 1, 'Name', 'conv1') batchNormalizationLayer('Name', 'bn1') reluLayer('Name', 'relu1') convolution1dLayer(3, numFilters, 'Padding', 'causal', 'DilationFactor', 2, 'Name', 'conv2') batchNormalizationLayer('Name', 'bn2') reluLayer('Name', 'relu2') bilstmLayer(lstmUnits, 'OutputMode', 'last', 'Name', 'bilstm') dropoutLayer(dropout, 'Name', 'drop') fullyConnectedLayer(numClasses, 'Name', 'fc') softmaxLayer('Name', 'softmax') classificationLayer('Name', 'output') ];这里只是个简化示意,完整实现里需要用addLayers和connectLayers把残差连接接好。TCN实操中有个细节:causal padding是必要的,保证卷积输出不会出现未来信息的泄漏,这对时序分类任务很关键。
4.3 训练策略与防过拟合细节
训练过程我用了三个技巧:学习率热启动衰减、早停和类别权重均衡。热启动的意思是刚开始学习率比较大(比如0.01),训练几个epoch后按余弦曲线衰减到初始值的十分之一。这个策略在轴承数据集上比固定学习率收敛更快,也更容易冲出初始较差的局部区域。
早停我用验证集损失做监控,连续12个epoch没有下降就停止训练并恢复最佳权重。MATLAB里通过设置trainingOptions的ValidationPatience参数实现。
类别不均衡在故障诊断里很常见,正常样本远远多于故障样本,特别是启动阶段采集的数据。在classificationLayer前面设置类权重,或者在数据采样阶段做滑动窗口重采样都可以解决欠拟合少数类的问题。我的习惯是窗口重叠率设为0.5做样本扩增,代码上比较干净,不需要改网络结构。
训练完成后一定要对模型做混淆矩阵评估,不要只看整体准确率。轴承故障诊断里最怕的是把外圈故障误判为正常,这种错误在整体准确率里看不出来,但在混淆矩阵里会表现出明显的对角线外分布,需要用焦点损失或用代价敏感学习针对性调整。
5. SHAP可解释性分析:让故障诊断不再“黑盒”
5.1 SHAP值的计算逻辑
SHAP值来自合作博弈论中的Shapley值。它的基本思想是:把每一个输入特征看作博弈参与者,模型输出看作参与者的总收益,Shapley值描述的是每个参与者在所有可能组合中对收益的平均边际贡献。放到故障诊断里,就是计算每个CMBE特征对“模型判定为内圈故障”这个输出的贡献量。
具体到计算上,假设有一个特征向量x,我们要算第j个特征的重要度。做法是枚举所有包含第j个特征的特征子集S,分别计算加入第j个特征前后模型输出概率的变化,然后对所有可能的子集组合按合作博弈的权重做平均。这个计算在特征数量超过十几个时会变得极其庞大,所以实际实现中采用的是基于采样和蒙特卡洛的近似算法,通过多次随机抽样逼近真值。
MATLAB的Statistics and Machine Learning Toolbox提供了shapley函数,可以直接对机器学习模型计算Shapley值。对于深度学习模型,需要先包装一个预测函数,把输入特征序列转换成分类得分。代码思路如下:
% 假设 net 是训练好的网络,X 是特征矩阵 modelsh = shapley(@(x) predict(net, x), X_train, 'QueryPoint', X_test(1,:));如果你的输入特征还是原始时序数据而不仅仅是CMBE特征向量,SHAP计算会面临维度灾难。所以实际项目中我都是在特征维度上做SHAP,也就是对160维的CMBE特征做可解释性,而不是对原始振动点做可解释性。这样既能说明哪些尺度上的复杂度特征更关键,又能控制计算开销。
5.2 MATLAB中的实现与可视化
MATLAB中执行SHAP分析后,主要输出三张图:特征重要性柱状图、单个样本的SHAP力图、每个特征随取值变化的SHAP依赖图。特征重要性柱状图直接按平均绝对SHAP值降序排列,可以直观看到哪些尺度和哪些IMF贡献最大。SHAP力图则展示特定样本中,哪些特征把预测结果推向某个故障类别,哪些特征在反向拉低这个类别的概率。
我分享一个自己的分析过程。有一次模型把外圈故障样本错判成内圈故障,我选出错判样本画SHAP力图,发现第3个IMF在第15个尺度上的熵值贡献异常高,把决策推向了内圈。后来检查发现,该样本的故障尺寸较大,冲击成分在信号里产生了很强的调制边带,恰好和某型内圈故障的调制特征相似。这就说明模型不是随机犯错,而是被真实物理特征误导,可解释性直接帮我定位到了误判原因。
5.3 可解释性对实际产线排障的意义
在实验室里做诊断,准确率高就够了;在产线上做诊断,准确率之外还得回答“为什么”。一旦把模型部署到现场,工程师不可能信任一个只会给概率的神秘盒子。SHAP让维修人员能直接看到:设备当前状态主要是由哪些频带、哪些尺度上的异常复杂度引起的。如果某个尺度特征持续走高,就可以提前安排检查对应部件,而不是等报警阈值超限再停机。
SHAP另一个用途是反向验证特征提取合理性。如果分析结果显示高贡献的特征集中在前三个IMF的高尺度区间,说明分解尺度选择没有问题;如果贡献最高的是最后几个IMF,那通常意味着噪声分量没有被过滤干净,模型正在用噪声模式做决策,这时需要回头缩减IMF保留数量并重新训练。
6. 完整代码框架与踩坑实录
6.1 代码模块划分
整个MATLAB项目我按六个脚本组织,每个脚本只负责一个环节,方便改动和复现。第一个是数据加载脚本,负责读取原始振动信号并完成标准化、滑窗切分和数据集划分;第二个是ICEEMDAN分解脚本,处理所有信号的分解并把IMF保存成mat文件;第三个是CMBE特征提取脚本,读取IMF批量计算特征矩阵;第四个是贝叶斯优化脚本,定义超参搜索空间和目标函数;第五个是模型训练与评估脚本,输出准确率、F1、混淆矩阵;第六个是SHAP分析脚本,加载训练好的模型进行可视化。
把数据划分放在最开始特别重要。有人习惯先分解再划分,结果训练集和验证集之间可能出现同一条信号不同窗口的泄漏,导致验证精度虚高。我强烈建议先按设备或按工况划分,确保同一段原始信号的所有时间窗都在同一集合里,再去做分解和特征提取。
6.2 常见问题速查表
| 现象 | 可能原因 | 解决方法 |
|---|---|---|
| ICEEMDAN分解速度极慢 | NE过大或信号太长 | 将NE设为60~80,信号先降采样到20000Hz以内 |
| 分解后IMF数量不稳定 | 信号长度差异大 | 统一截取相同长度,按能量排序后仅保留前K个IMF |
| CMBE特征计算报维度错误 | 不同IMF长度不一致 | 在粗粒化前统一对IMF做长度对齐或裁剪 |
| 贝叶斯优化训练时间过长 | 目标函数里每次完整跑训练 | 加入早停,并将MaxObjectiveEvaluations控制在25~30 |
| 模型验证精度高但测试精度低 | 数据窗口重叠导致泄漏 | 检查数据划分是否按信号段切分,而不是按窗口随机切分 |
| SHAP计算内存不足 | 特征维度过高且QueryPoint太多 | 只对误分类样本计算,或对特征做PCA降维后再分析 |
| GPU训练OOM | 序列过长或batch过大 | 减小miniBatchSize,或改为CPU训练对比 |
6.3 个人调试心得
第一,不要一开始就追求最复杂的配置。我建议先用少量样本、小网络跑通全流程,确认每个模块的数据形状能正确传递,再加数据量扩样本。全流程联调时最容易出问题的地方反而是CMBE特征矩阵的维度约定:网络的输入节点必须和滑窗提取的特征数严格对齐,差一个数字都会在训练时报维度错误。
第二,记录每一步的中间文件格式。分解结果、特征矩阵、训练日志统一用mat文件保存并标注版本号,否则调整一次参数就要从头重算,大量时间花在重复分解上。
第三,固定随机种子并且把数据划分结果保存下来。复现故障诊断实验时,随机性的影响比其他任务更大,因为数据量相对少,随机划分窗口会导致最后的结果上下浮动几个百分点。保存好数据划分索引之后,后续任何对比实验都在同一套划分上跑,结果之间才有可比性。
第四,SHAP分析不是只在最后跑一次就完事。我习惯在调参过程中每隔几次实验就对错误分类样本做一次SHAP,看模型是否有稳定的不合理依赖。如果某次改动后SHAP最重要的特征从物理相关特征跳到了无关特征,基本可以判断训练出了异常,省掉很多后期排查时间。
我现在回头看,这套ICEEMDAN-CMBE加上Bayes-TCN-BiLSTM和SHAP的组合,最大的价值不是单个模型精度能拉到多高,而是把整个诊断过程从信号分解到模型决策变成了一个可审计、可追溯的流程。对我个人来说,最有成就感的一刻是某次现场数据复测中,模型对早期外圈故障给出了高置信度报警,而SHAP图上贡献最大的特征正好对应故障特征频率所在尺度区间,这种“结果+依据”双重验证的方式,让现场工程师第一次真正接受了这个模型。后续如果再扩展,我打算把SHAP分析与在线监测系统对接,在设备运行过程中实时输出关键特征的贡献度变化曲线,这样可以在故障萌芽阶段就给出可解释的预警信号。