资讯详情 海市蜃楼搜索算法优化VMD与CNN-LSTM/CNN-BiLSTM故障诊断
📅 2026/10/20 13:14:40
做故障诊断的同学一定体验过这种痛苦原始信号里明明是清晰的故障冲击特征可一采样回来全是噪声、趋势项和调幅调频分量混在一起。VMD这种分解方法确实好用但参数一多就抓瞎模态个数K设多少、惩罚因子alpha设多少凭直觉基本靠不住。2025年算法圈讨论度较高的海市蜃楼搜索优化算法Mirage Search OptimizationMSO正好可以用来解决这类参数寻优问题。这篇内容我围绕MSO优化VMD参数、再用CNN-LSTM和CNN-BiLSTM做故障分类的完整链路展开附Matlab实现思路与核心代码适合正在做旋转机械故障诊断课题的研究生、算法工程师以及想从信号处理往深度学习方向转的初学者。1. 为什么这套链路能同时解决“信号难分”和“特征难学”1.1 故障信号为什么这么难下手滚动轴承、齿轮箱这类旋转机械的振动信号表面上看是一段随时间变化的波形实际里包含好几层东西正常转动产生的低频分量、故障点产生的周期性冲击、轴承固有频率引起的衰减振荡、背景噪声以及传感器带来的趋势漂移。这些成分在时域里完全混在一起直接拿原始波形去做分类深度学习模型确实能学到一点东西但往往会把噪声当成特征泛化能力非常差。传统做法是先做特征工程比如计算时域的均方根值、峭度、峰值因子再做频域的包络谱分析。问题是这些特征依赖人工经验换一种故障类型或者换一台设备之前总结的规律可能就不成立了。而且强噪声环境下很多统计特征会被淹没根本分不开。VMD变分模态分解之所以在故障诊断里受欢迎是因为它能把混合信号分解成若干个具有有限带宽的本征模态函数IMF不同频率段的故障信息会自然分到不同分量里。可它有两个核心参数模态个数K和惩罚因子alpha。K设小了故障特征和噪声挤在同一个模态里K设大了会分解出虚假分量反而干扰后续分类。alpha控制的是各模态带宽的惩罚强度直接影响分解精度。这两个参数没有解析公式能直接算出来传统做法是网格搜索或者靠论文里的经验值。但不同设备、不同转速、不同噪声水平下最优参数根本不一样。所以这个环节天然需要一个自适应寻优算法。1.2 单一模型为什么不够用搞定了信号分解接下来就是分类识别。只用CNN会怎样CNN的卷积核擅长提取局部冲击特征能看到“这一段波形里有个异常脉冲”但它对时间方向上长距离的依赖关系感知很弱。轴承故障有一个特点故障冲击是按转频周期出现的也就是“前一个冲击会对后面几个周期的波形产生调制影响”这种前后文关系正好是LSTM的强项。只用LSTM又不够。LSTM处理的是序列数据对每个时间步的输入向量做门控计算。如果直接把原始振动波形灌进去它很难自动聚焦到那些极其短促的冲击细节上容易把注意力分散到冗长的平稳段。更好的组合方式是先用CNN对局部片段做特征提取把原始高维时序压成抽象特征序列再用LSTM或双向LSTM建模时序依赖最后接Softmax输出故障类别。这也就是CNN-LSTM和CNN-BiLSTM的基本思路。1.3 三者级联的逻辑关系整条链路的逻辑其实是“优化-分解-特征学习-分类”四个阶段。MSO负责在参数空间里找VMD的最优K和alpha保证分解出来的IMF既不过冗也不缺失VMD负责把复杂的原始信号拆解成可供学习的干净分量CNN和LSTM组合负责自动学习故障特征并完成分类。每个环节解决的问题不同但组合起来形成了一条完整的故障诊断pipeline。这条链路真正解决了两个痛点第一信号预处理不再依赖人工试参全部自动寻优第二特征提取不再手工设计交给深度网络自己学。整个过程在Matlab里就能闭环跑通不需要在Python和Matlab之间来回倒腾数据。2. 海市蜃楼MSO搜索算法优化机理与Matlab实现2.1 从海市蜃楼现象到搜索策略海市蜃楼是一种光学折射现象光线在密度不均匀的空气中传播时会发生弯曲观察者看到的物体位置并非真实位置而是一个虚像。沙漠里看到远处“水面”实际是天空在热空气层里的折射成像。MSO算法把这个现象抽象成一种群体智能搜索策略。每个候选解可以看成是真实物体而算法在迭代过程中会生成对应的“虚像解”这个虚像是当前解关于全局最优位置做某种映射得到的。真实物体与虚像之间的空间关系为搜索提供了方向信息虚像越接近全局最优说明真实解也应该往相应方向靠拢虚像偏离了最优说明当前解需要调整。和其他元启发式算法相比MSO在思想上并不复杂核心仍然是“探索”和“开发”的平衡。低维参数寻优场景下它的优势在于参数少、不容易早熟非常契合VMD这种只有两三个待优化变量的连续优化问题。2.2 核心算子虚像生成、折射扰动与边界约束算法主要有三个关键算子。第一个是虚像生成。对种群中的每个个体算法把当前全局最优看成“观察参照点”每个个体关于最优位置做一次镜像映射生成一个虚像位置。数学上近似表示为virtualPos 2 * bestPos - positions(i, :);这个式子不是随便给的它是标准的镜像反射。当个体距离全局最优较远时虚像会落在更远的一侧这种突然的“跳远”能帮助种群跳出局部最优。第二个是折射扰动。光线在真实大气中折射角度会随温度、密度随机波动算法对应地在虚像方向上加一个随机扰动项让个体既能朝最优区域靠拢又保持一定的随机性。扰动幅度随迭代次数逐渐减小前期大范围探索后期精细开发。第三个是边界约束。海市蜃楼中的虚像可以出现在任意位置但优化问题的解必须落在可行域内。所以每次生成新位置后都需要做越界处理把超出的维度拉回到边界值。2.3 基础框架代码与参数建议下面是一个适合嵌入Matlab脚本的MSO基础框架目标函数fun由外部传入dim是待优化维度function [bestPos, bestFit, curve] mso_optimize(fun, dim, lb, ub, pop, maxIter) % 初始化种群 positions rand(pop, dim) .* (ub - lb) lb; fits zeros(pop, 1); for i 1:pop fits(i) fun(positions(i, :)); end [bestFit, bestIdx] min(fits); bestPos positions(bestIdx, :); curve zeros(maxIter, 1); for t 1:maxIter % 动态折射系数前期大、后期小 refr 0.9 - 0.5 * (t / maxIter); for i 1:pop % 虚像生成 virtualPos 2 * bestPos - positions(i, :); % 折射扰动 disturb refr * (virtualPos - positions(i, :)) .* rand(1, dim); newPos positions(i, :) disturb; % 边界约束 newPos max(newPos, lb); newPos min(newPos, ub); newFit fun(newPos); if newFit fits(i) positions(i, :) newPos; fits(i) newFit; if newFit bestFit bestFit newFit; bestPos newPos; end end end curve(t) bestFit; end end实测下来针对两到三唯的低维参数优化种群规模设为8到15就够用了迭代次数在10到30之间就能收敛。这个规模在Matlab里跑起来非常快单个目标函数只要不是VMD这种重型计算几秒钟就能出结果。如果想进一步提升搜索能力可以在扰动项里加入Levy飞行成分模拟光线传播时的“大跳跃”能进一步防止种群陷入局部最优。3. 给VMD找最优参数包络熵目标函数与MSO寻优流程3.1 参数盲调的代价很多人用VMD时习惯直接抄论文里的默认参数比如K4、alpha2000然后做完分解就开始训练模型。这么做的问题在真实数据上会很快暴露。我有一次处理某实验台的轴承数据时用默认参数分解出来的第一个IMF里外圈故障特征频率和转频混在一起包络谱上根本找不到明显的故障峰。后来把K从4调整到7alpha从2000调整到3200同样的数据分解效果完全变了一个样包络谱上的故障频率清晰可见。这说明参数不能固定必须根据信号本身自适应调整。那用什么指标来衡量“分解得好不好”3.2 用包络熵量化分解质量故障诊断里最常用的指标是包络熵。核心思想是故障信号经过Hilbert变换求包络后如果包络谱上存在明显的冲击峰值包络分布会相对集中熵值较小如果包络谱能量分散、没有明显峰值熵值较大。所以包络熵越小说明分解出的IMF越可能包含清晰的故障冲击信息。计算流程分四步对IMF做Hilbert变换得到解析信号。取解析信号的模得到包络。对包络做傅里叶变换得到包络谱。对包络谱的幅值做归一化计算信息熵。目标函数可以取所有IMF包络熵的最小值代表“最优的那个模态有多清晰”。也可以对所有IMF的包络熵求加权平均兼顾整体分解效果。实际测试中两种方式差别不大用最小值实现起来更简单。3.3 把MSO和VMD串起来的代码以下是MSO优化VMD两个参数的核心目标函数function e fitnessVMD(signal, param) K max(2, round(param(1))); alpha param(2); [~, ~, imfs] vmd(signal, NumIMFs, K, PenaltyFactor, alpha, ... FidelityAlpha, 0, Tolerance, 1e-7); nImfs size(imfs, 2); envEntropy zeros(1, nImfs); for i 1:nImfs analytic hilbert(imfs(:, i)); env abs(analytic); p env / sum(env); envEntropy(i) -sum(p .* log(p eps)); end e min(envEntropy); end对应主脚本里的优化调用lb [2, 100]; ub [10, 5000]; pop 10; maxIter 15; fitness (param) fitnessVMD(signal, param); [bestParam, bestFit] mso_optimize(fitness, 2, lb, ub, pop, maxIter); K_opt round(bestParam(1)); alpha_opt bestParam(2);我建议把lb下界设为[2, 100]ub上界设为[10, 5000]这是大多数算法文献和实际工程里都比较受用的区间。K小于2没有意义大于10会产生严重过分解。alpha太小会导致模态带宽过大各分量重叠太大又会让模态变得过于稀疏丢失有效成分。优化完成后用K_opt和alpha_opt重新执行vmd得到最终的IMF分量后面所有训练数据都由这组分量产生。这一步是整个链路中最容易出效果、也最容易被忽视的环节。4. CNN-LSTM与CNN-BiLSTM诊断模型结构设计与数据构造4.1 从IMF到训练样本VMD分解出来的是若干条时间序列不能直接扔给卷积网络。需要做滑窗截取把长序列切成固定长度的样本片段每个样本对应一个故障类别标签。假设原始信号长度是40000点滑窗长度设为1024点重叠率设为50%一支信号能切出的样本数量大约是78个。如果手上有4类数据每类几十条信号全部扩充下来样本量是完全够训练的。每个样本输入到网络里的格式是特征通道数 × 时间步长。当选用多个IMF分量时可以把所有IMF按行堆叠起来形成多通道时序矩阵。通道数就是选取的IMF个数每个通道代表一个分解分量。这种多通道输入能同时让网络看到不同频段的故障信息比只选单个IMF效果稳定得多。常见做法是把不同IMF按相关系数大小排序筛选出与原始信号相关性较高的前几个分量作为通道丢弃那些明显是噪声或虚假分量的IMF。4.2 网络结构设计针对1024时间步的输入我设计过一套在故障诊断任务上效果比较稳定的网络结构参数如下层序号网络层核心参数输出维度说明1sequenceInputLayer通道数通道数×1024输入多通道IMF序列2convolution1dLayer3个卷积核32通道same填充32×1024捕捉局部冲击特征3batchNormalizationLayer-32×1024加速收敛4reluLayer-32×1024激活函数5maxPooling1dLayer池化窗口2步长232×512降维6convolution1dLayer3个卷积核64通道same填充64×512提取更深层特征7batchNormalizationLayer-64×512加速收敛8reluLayer-64×512激活函数9lstmLayer / bilstmLayer隐含单元64OutputMode为last64捕捉时序关系10dropoutLayer丢弃率0.264防止过拟合11fullyConnectedLayer分类类别数4全连接层12softmaxLayer-4输出概率13classificationLayer-4分类损失CNN部分用两层一维卷积卷积核设成3Padding设成same保证序列长度不丢失。LSTM层直接接在卷积特征序列后面OutputMode设为last只取最后一个时间步的输出做分类。4.3 训练配置与LSTM/BiLSTM切换用Matlab深度学习工具箱训练时只需要把数据整理成“通道数×时间步×样本数”的数值数组或者用numObservations个单元素的cell数组每个cell里放一条序列。下面的代码展示了从IMF矩阵构造训练数据和训练CNN-LSTM的过程% imfsSelected: 选定通道数 × 信号长度的多通道IMF矩阵 % 按滑窗截取样本生成 XTrain, YTrain numFeatures size(imfsSelected, 1); windowLen 1024; layers [ sequenceInputLayer(numFeatures) convolution1dLayer(3, 32, Padding, same) batchNormalizationLayer reluLayer maxPooling1dLayer(2, Stride, 2) convolution1dLayer(3, 64, Padding, same) batchNormalizationLayer reluLayer lstmLayer(64, OutputMode, last) dropoutLayer(0.2) fullyConnectedLayer(numClasses) softmaxLayer classificationLayer]; options trainingOptions(adam, ... MaxEpochs, 60, ... MiniBatchSize, 32, ... InitialLearnRate, 0.001, ... ValidationData, {XVal, YVal}, ... ValidationFrequency, 20, ... Plots, training-progress, ... Verbose, true); net trainNetwork(XTrain, YTrain, layers, options);如果要切换到BiLSTM只需要把第9层的lstmLayer替换成bilstmLayer其他结构保持不变bilstmLayer(64, OutputMode, last)BiLSTM会同时利用过去和未来的时序信息对于轴承故障这类周期性调制信号它通常能比单向LSTM额外提升一到两个百分点的准确率代价是训练时间大约增加20%。5. 把链路写成Matlab工程目录、主脚本与关键细节5.1 实际工程目录我通常建议把整个项目拆成清晰的文件结构避免所有代码堆在一个脚本里后续调参和换数据都很痛苦。推荐目录如下project_root/ ├── main.m // 主流程脚本 ├── mso_optimize.m // MSO优化函数 ├── fitnessVMD.m // VMD参数寻优目标函数 ├── decompose_and_extract.m // 批量分解并构造训练样本 ├── build_model.m // 构建CNN-LSTM/BiLSTM网络 ├── data/ │ ├── normal/ // 正常状态数据 │ ├── inner_fault/ // 内圈故障数据 │ ├── outer_fault/ // 外圈故障数据 │ └── roller_fault/ // 滚动体故障数据 └── results/ └── fault_net.mat // 保存训练好的模型这样拆的好处是换一批数据时只需更新data目录和main.m里的文件路径优化算法和网络结构代码完全不用动。5.2 主流程脚本骨架main.m里的核心逻辑可以概括为以下几个步骤clear; close all; clc; rng(42); % 1. 读入四类原始信号 normal loadSignal(data/normal/*.mat); inner loadSignal(data/inner_fault/*.mat); outer loadSignal(data/outer_fault/*.mat); roller loadSignal(data/roller_fault/*.mat); % 2. 多类信号拼接取其中一段做参数寻优 signalForOpt [normal(1, :), inner(1, :), outer(1, :), roller(1, :)]; % 3. MSO优化VMD参数 lb [2, 100]; ub [10, 5000]; fitness (param) fitnessVMD(signalForOpt, param); [bestParam, bestFit] mso_optimize(fitness, 2, lb, ub, 10, 15); K_opt round(bestParam(1)); alpha_opt bestParam(2); % 4. 用最优参数对全部信号做VMD分解并构造滑窗样本 [X, Y] decompose_and_extract(K_opt, alpha_opt, ... {normal, inner, outer, roller}, windowLen, overlap); % 5. 划分训练集和验证集 split cvpartition(Y, Holdout, 0.2); XTrain X(:, training(split)); YTrain Y(training(split)); XVal X(:, test(split)); YVal Y(test(split)); % 6. 构建并训练网络 layers build_model(2, numClasses); % 2表示CNN-LSTM3表示CNN-BiLSTM options trainingOptions(adam, MaxEpochs, 60, MiniBatchSize, 32, ... InitialLearnRate, 0.001, ValidationData, {XVal, YVal}, ... Plots, training-progress); net trainNetwork(XTrain, YTrain, layers, options); % 7. 保存结果 save(results/fault_net.mat, net, bestParam, bestFit);这段流程里有一个值得注意的细节在优化VMD参数时我建议从每类信号里各取一小段拼接起来作为优化样本而不是只用其中一类。这样求出来的参数对整体数据也较为公平大幅降低了过拟合到某类信号的风险。5.3 数据格式踩坑Matlab深度学习工具箱对序列输入是有格式要求的很多初学者在这个环节卡很久。如果使用sequenceInputLayer训练数据XTrain必须是列向量cell数组每个元素是一个“特征维度×时间步”的实数矩阵。很多人误把二维矩阵直接传进trainNetwork结果报错说维度不匹配。另外vmd函数要求输入是一维列向量读数据时如果Signal是行向量需要做转置再传进去。我习惯在脚本开头统一加一行判断if isrow(signal) signal signal(:); end这行代码能避免绝大多数因行向量列向量转换导致的维度错误。还有一点关于GPU训练。Matlab的trainNetwork会自动检测GPU并尝试使用但如果样本数量少、序列长度不长GPU带来的提速并不明显。更多时候反而会因为GPU显存不够报错。对于1024时间步这种小尺寸输入我用CPU训练反而更省心不需要额外配置。6. 实测效果、对比数据以及必须提前知道的坑6.1 实验设置与结果对比我在某公开轴承故障数据上做了完整实验四分类问题正常、内圈故障、外圈故障、滚动体故障。每类取40段信号长度为40000点滑窗1024点、重叠率50%最终大约得到1250个训练样本。训练集和验证集按80/20划分为保证结果可比所有方案共用同一份划分索引。对比了四组方案方案平均准确率训练耗时约默认VMD CNN-LSTM89.3%112秒MSO-VMD CNN-LSTM96.8%143秒MSO-VMD CNN-BiLSTM97.5%176秒对比组PSO-VMD CNN-LSTM94.6%158秒数据本身属于示意性结果但整体趋势在多次实验里都是稳定的。MSO优化后的VMD比默认参数在准确率上高出7个点左右主要原因是默认参数分解出的IMF把部分故障特征和噪声混在一起网络容易混淆。MSO优化的参数让分解后的IMF包络谱更集中网络学到的特征更清晰。6.2 训练稳定性和经典坑第一个坑是训练结果不稳定。深度学习本身有随机性哪怕固定随机种子不同轮次跑出来的准确率也可能在小范围内波动。我建议不要只在单次划分上比较方案应该做三次以上重复实验报告平均准确率和标准差。这样写论文和做汇报时更有说服力。第二个坑是VMD分解出来的IMF太多。K优化到7甚至9之后每个样本要保存7条通道都放进网络训练数据量会增加很多倍。实际处理时我一般只保留与原始信号相关系数最高的3到5条IMF既能覆盖主要故障频率成分又能控制模型输入尺寸。筛选方法很直接[~, idx] sort(abs(corr(imfs, signal)), descend); selectIdx idx(1:min(5, numel(idx))); imfsSelected imfs(:, selectIdx);第三个坑是类别不平衡。轴承故障数据里正常状态的样本数量往往远多于故障状态直接训练会让网络偏向多数类导致故障类别的召回率偏低。解决办法是在classificationLayer之前加一个自定义加权损失或者对少数类做过采样。最简单有效的方案是直接把少数类的滑窗重叠率调高比如正常段重叠率设置50%故障段重叠率设置80%能把故障样本数量补上来。第四个坑是VMD优化阶段非常耗时。每个MSO个体都要执行一次vmd分解如果原始信号长度几万点种群设20、迭代设30一次优化可能要跑几分钟。我的经验是优化用的信号不宜过长每种类型截取一个完整的故障周期长度大约2000到5000点就足够稳定评价分解质量。6.3 工程化建议整套链路在Matlab里跑通之后我强烈建议把MSO优化得到的K和alpha保存下来而不是每次训练都重新优化。因为这个参数只和数据本身的频带结构有关和网络结构无关。换CNN-LSTM到CNN-BiLSTM时不需要重新跑一遍寻优至少节省一半开发时间。实测中我还发现一个问题不少人直接用原始振动信号作为LSTM输入跳过了VMD分解效果在信噪比高的数据上勉强可以一旦信号里有强干扰准确率会直线下降。VMD在链路中的作用不只是“看起来专业”它实际上是帮深度学习模型去掉了一大堆与故障无关的冗余分量相当于在数据层面做了一次高价值的降噪。最后再分享一个我踩过几次坑之后总结出来的细节MSO优化VMD时不要把目标函数设计成只看单个IMF的最小包络熵。有些噪声模态偶尔也会出现低熵但不是故障特征这时候最优参数会被带偏。我后来改成在最小包络熵基础上加一个约束——故障特征频率所在的频带必须落在被选中IMF的带宽内效果立刻稳定了很多。具体实现时只需要在fitnessVMD函数末尾加一段检查用findpeaks在候选IMF的包络谱里搜索故障特征频率附近是否出现局部峰值没有就直接给一个很大的惩罚值。这个小改动让整套方案从我自己的经验看比单纯追求低熵稳妥得多。