MSO-VMD-CNN-LSTM/BILSTM故障诊断方案详解

去年我做设备状态识别项目的时候，拿到现场振动数据后第一反应是直接上卷积神经网络，结果效果并不理想。后来换了思路，先把信号做变分模态分解，再用优化算法自动找分解参数，最后把CNN和LSTM/BILSTM串起来做分类，测试集准确率从83%左右提到了96%。这套流程就是标题里写的MSO-VMD-CNN-LSTM/BILSTM故障诊断方案。其中MSO是海市蜃楼群优化算法，用来解决VMD参数难确定的问题；CNN负责从分解后的信号里提取局部特征，LSTM或BILSTM负责继续学习时序依赖；整个实验在Matlab里跑通。这篇文章就是把这套方案的思路、算子设计、代码实现和踩坑过程完整拆开讲一遍，适合正在做机械故障诊断、轴承和齿轮箱状态识别方向的研究生或工程师参考。

1. 为什么故障诊断要倒腾成“优化+分解+深度学习”三步走

1.1 直接把原始信号喂给深度模型为什么不好使

很多人拿到振动信号后，第一反应就是把原始波形切成长度为1024或者2048的样本，然后直接丢给CNN或者LSTM。这样做不是完全不行，但很容易卡在一个尴尬的准确率上，尤其是当信号里混着大量工况噪声、转频干扰和随机冲击的时候。

原始振动信号本质上是一个非平稳、非线性的时间序列。故障冲击成分往往被淹没在背景噪声里，幅度可能只占信号总能量的百分之几。CNN虽然能靠卷积核提取局部特征，但面对这种信噪比很低的输入，第一层卷积核学到的往往是噪声模式而不是故障模式。LSTM也有类似问题，直接处理长序列时，冗余信息太多，门控单元很难在几百个时间点里定位到几个微弱的冲击峰值。

另外还有一个实际问题：故障样本通常少。现场设备大多数时间处于正常状态，真正故障状态的数据可能要等很久才能积累到足够的量。在样本量不足的情况下，直接把原始高维信号喂给深层网络，非常容易过拟合。训练集准确率能到99%，测试集只有80%出头，这种情况我见过太多次了。

所以，问题不在深度学习模型本身，而在输入表示。如果能在进入网络之前，把故障相关的关键成分从混合信号里分离出来，把“找特征”的难度降下来，后面的网络就能把有限的拟合能力用在真正重要的模式上。

1.2 VMD在这里面到底干了什么事

VMD全称是变分模态分解，它可以把一个复杂信号自适应地分解成若干个本征模态函数IMF。每一个IMF对应信号在某个中心频率附近的一个调幅调频成分，这就像想办法把一个混合声音拆成人声、吉他、鼓点等多个独立的音轨。

和EMD（经验模态分解）相比，VMD最大的优势是有一个严格的数学框架。EMD靠极值包络插值递归筛选，对噪声敏感，容易出现模态混叠；VMD则是在频域里通过迭代求解变分问题，把各模态的中心频率和带宽自动找出来。只要模态数K和惩罚因子alpha设置合理，VMD分解出来的各个IMF在频域上是相对分离的，这对接下来的特征提取非常有利。

在故障诊断场景下，VMD的主要作用是“放大故障特征”。比如滚动轴承外圈故障，信号里会出现周期性的冲击成分，这些冲击对应着高频共振频带；齿轮故障则对应着啮合频率及其边频带。经过VMD分解之后，故障冲击成分通常会集中到某一个或某几个IMF里，而噪声被分散到其他分量中。这样，后续网络要学的模式就清晰多了。

1.3 为什么要在VMD前面加一个MSO寻优

VMD有两个参数对分解效果影响极大：一个是模态数K，另一个是惩罚因子alpha。

K设得太小，会欠分解，两个频率接近的成分叠在同一个IMF里，故障特征被掩盖；K设得太大，会过分解，一个完整成分被硬生生拆成好几个虚假分量，引入额外的干扰。alpha控制各模态的带宽约束，alpha越大，带宽约束越强，各模态的中心频率更容易被限制在很窄的范围内；alpha太小，模态带宽会变宽，容易出现频谱重叠。

人工调这两个参数不是不行，但很费劲。因为实际故障信号的频率结构是未知的，不同工况、不同转速、不同损伤程度下，最佳K和alpha都在变。用网格搜索穷举K从2到10、alpha从1000到3000，每个组合都要做一次完整的VMD分解，非常耗时。

这里就该启发式优化算法出场了。把VMD的两个参数当作优化变量，把某种能反映“分解效果好坏”的指标当作适应度函数，用优化算法自动搜索最优参数组合。MSO就是其中一个比较新的选择。它模拟海市蜃楼形成过程中“虚像折射”和“蜃景迁移”两个机制，在搜索能力和收敛速度之间做了一个相对平衡的取舍。

2. MSO海市蜃楼优化算法的实现细节与调参

2.1 海市蜃楼现象与算法算子的对应关系

海市蜃楼大家都不陌生，沙漠里看到远处好像有一片湖，湖面上倒映着山影，其实那是光线穿过不同密度的空气层后发生折射形成的虚像。虚像的位置和真实物体之间并没有直接对应关系，它是由大气折射率梯度决定的。

MSO的全称我把定为Mirage Swarm Optimization，也就是海市蜃楼群优化算法。它把每个候选解看成是“蜃景粒子”，把真正的最优解看成是“实体”。算法要做的就是不断根据虚像的漂移和折射特征，反向逼近实体位置。

在这个映射关系下，算法设计了两个核心算子：

• 折射算子：模拟光线在大气层中弯曲传播的过程。每个粒子以当前最优解为参照，生成一条偏向最优方向的曲线搜索路径，相当于在局部区域做细致的开发；
• 迁移算子：模拟蜃景虚像随观测角度变化在视觉上发生漂移的过程。部分粒子的位置会向群体历史平均位置方向移动，保留一定的随机性，相当于维持群体多样性、避免过早收敛。

这两个算子一个偏向开发，一个偏向探索，配合种群的迭代更新，就能在参数空间里自动寻优。

2.2 种群初始化、折射算子和迁移算子的计算逻辑

先说明，这里给出的是我在实际实现时使用的一种简化版本。不同的人实现MSO可能会在算子公式上有差异，但核心思路是一致的，你可以按照自己的问题做调整。

种群初始化时，假设优化维度是2，分别对应VMD的K和alpha。设置种群规模N，最大迭代次数MaxIt。每个粒子的位置在上下界之间随机生成：

N = 20; % 种群规模 MaxIt = 30; % 最大迭代次数 dim = 2; % 优化维度 [K, alpha] lb = [3, 1000]; % K最小值，alpha最小值 ub = [10, 3000]; % K最大值，alpha最大值 pop = lb + rand(N, dim) .* (ub - lb);

在上面的初始位置里，K通常会被当连续值处理，但在调用VMD之前会做四舍五入，因为模态数必须是整数。

迭代过程中，每个粒子的新位置主要用两个公式更新。第一个是折射更新：

currentBest = bestPos; % 当前全局最优位置 % R为折射系数，随迭代递减 R = Rmax - (Rmax - Rmin) * (t / MaxIt); for i = 1:N beta = 0.5 + rand * 0.5; r = randn(1, dim); newPos(i, :) = currentBest + beta * (currentBest - pop(i, :)) .* r; end

R在这里相当于大气折射强度的变化。迭代初期R比较大，粒子探索范围更广，可以跑到离当前最优较远的地方；迭代后期R变小，粒子主要在最优解附近精细搜索。这种参数衰减思路和惯性权重类似，但把随机方向乘到了差值项上，搜索轨迹不完全是一条直线，确实不容易卡在局部极值。

第二个是迁移更新。每隔一定迭代步数，取一部分粒子向群体平均位置靠拢：

meanPos = mean(pop, 1); idx = randperm(N, round(N * 0.3)); pop(idx, :) = pop(idx, :) + rand(size(pop(idx, :))) .* (meanPos - pop(idx, :));

这一步的作用相当于让部分“虚像”重新向群体中心漂移，避免所有粒子都挤在同一个局部区域里。

2.3 适应度函数选择与Matlab代码片段

优化算法只是搜索工具，真正决定搜索结果好坏的是适应度函数。我常用的适应度是“包络熵值”，它反映信号的稀疏性。故障冲击成分在包络谱上会集中在少数几个频点附近，包络熵相对较低；而纯噪声的包络谱比较平坦，包络熵偏高。所以把包络熵的最小化作为搜索目标，是合理的。

计算包络熵的简化代码可以这么写：

function envEntropy = EnvelopeEntropy(x) analytic = hilbert(x); env = abs(analytic); envEnv = abs(hilbert(env)); % 对包络再做一次Hilbert，得到包络谱 p = envEnv(:); p = p / sum(p); % 计算信息熵 p = p(p > 0); envEntropy = -sum(p .* log(p)); end

注意，如果你直接对原始信号计算包络熵，会受到噪声影响；更好的做法是在每个粒子候选参数下，先对信号做VMD分解，然后选择主要模态的包络熵作为适应度。比如取前K个模态包络熵的平均值，或者取最小的包络熵值。我实际使用时，会在适应度函数里加一个模态混叠惩罚项，当相邻IMF中心频率的差值过小的时候，惩罚适应度。这样优化出来的参数不容易出现过分解。

调用VMD并计算适应度的核心流程如下：

function fitness = MSOCostFunction(signal, pos) K = round(pos(1)); alpha = pos(2); % 调用VMD函数 [u, ~, omega] = VMD(signal, alpha, 0, K, 0, 1, 1e-7); % 计算每个IMF的中心频率差，惩罚模态混叠 omegaDiff = diff(omega); overlapPenalty = sum(max(0, 100 - omegaDiff)); % 计算分解结果的平均包络熵 entropySum = 0; for i = 1:K entropySum = entropySum + EnvelopeEntropy(u(i, :)); end fitness = entropySum / K + overlapPenalty; end

2.4 MSO的完整搜索流程

整个MSO优化VMD参数的流程，我列在下面：

1. 初始化种群，随机生成N组(K, alpha)候选组合；
2. 对每组候选参数，调用VMD对原始信号做分解；
3. 根据IMF包络熵和模态混叠惩罚项计算每组参数的适应度；
4. 记录当前全局最优位置；
5. 进入迭代循环：
   • 更新折射系数R；
   • 用折射算子更新所有粒子的位置；
   • 每隔若干步，用迁移算子更新部分粒子的位置；
   • 重新计算每个粒子的适应度，更新全局最优；
6. 达到最大迭代次数后，返回全局最优位置，得到最优K和alpha。

这个流程在信号长度大概1秒、采样率12000Hz的情况下，30次迭代大概需要几分钟到十几分钟，完全可以接受。一旦参数确定下来，后面每段信号都可以用同一组参数分解，不用重复优化。

3. VMD分解与训练集构造的完整流程

3.1 VMD参数设定与分解效果判断

用MSO获得最优K和alpha之后，再对全部训练信号做一次性分解。这里要注意，不要对每一段信号都单独做一次MSO优化，一是耗时太长，二是每段信号都优化会导致不同样本的IMF数量不一致，模型没法统一输入格式。正确做法是从训练数据里选一段代表性信号做参数寻优，然后把得到的K和alpha固定下来，应用到所有样本的分解上。

调用VMD时的其他参数也需要固定好。我习惯这样设置：

tau = 0; % 噪声容忍度 DC = 0; % 不强制第一个分量为直流分量 init = 1; % 中心频率均匀初始化 tol = 1e-7; % 收敛容差 [IMF, ~, omega] = VMD(signal, alpha, tau, K, DC, init, tol);

分解结果怎么看？我一般会看两个东西。第一是中心频率序列，如果相邻IMF的中心频率差很小，说明K可能偏大。第二是各IMF的时域波形和频谱，如果某个IMF和原信号在频带上高度重叠，或者出现明显的“断裂感”，就需要调整alpha。

这里给一个实操判断标准：当相邻模态中心频率间隔小于100Hz时，基本可以认为K设大了。当然这个阈值和具体采样率、信号频带有关，大原则是中心频率应该从低到高均匀分布，且每个IMF的频谱主要能量集中在窄带内。

3.2 从分解结果到CNN-LSTM输入数据

VMD分解完之后，每个样本得到K个IMF序列。但这里有个关键问题：怎么把这些IMF变成CNN-LSTM能吃的输入？

我在项目里用的是“特征矩阵”思路。具体做法是：先把每个IMF切成若干个等长的小段，比如1秒信号切成8段，每段取有效值、峭度、峰值因子、包络峰峰值等4个特征，那么一个IMF就能得到一个8×4的特征矩阵。把K个IMF的特征矩阵叠起来，就得到了一个按时间排列的特征序列：8个时间步，每个时间步的维度是K×4。

这样做的好处是，输入LSTM的是经过压缩和提炼的特征，而不是原始的几万个采样点，训练速度快很多，也更容易收敛。CNN在这里可以作为一个局部特征提取层，对这个特征矩阵做1D卷积，进一步提取相邻时间步骤之间的关联信息。

如果你更追求端到端效果，也可以把每个样本的K个IMF直接堆成一个K×L的数字矩阵，当作一张“图像”输入2D-CNN。但这样做对LSTM的序列建模不太友好。我更推荐先用统计特征做降维，形成“时间步×特征维度”的结构，再接入CNN-LSTM。

3.3 标签处理与数据集划分

数据集的构造要特别小心“泄漏”问题。我一开始踩过一个大坑：直接把整个连续采集的信号做VMD，然后不分工况地在整段上开滑动窗口，结果训练集和测试集里可能包含来自同一个原始信号段的样本，模型在测试集上的表现虚高，一换新数据就原形毕露。

正确的做法是先把原始数据按工况切分成不同的记录段，例如一段正常数据、一段内圈故障数据、一段外圈故障数据。在每一段内部再做滑动窗口切片，保证窗口之间有重叠或不重叠都行，但训练集和测试集必须来自完全不同的记录段，不能交叉。

假设我们有正常、内圈故障、外圈故障、滚动体故障四类样本，每类80个样本，按7比3划分训练集和测试集。每个样本的标签就是对应的故障类型。把每个样本的IMF特征矩阵整理成固定形状，通常实现时先用一个元胞数组存放不同长度的特征，再填充到统一尺寸。

划分代码比较简单：

idxTrain = randperm(size(features, 1), round(size(features, 1) * 0.7)); idxTest = setdiff(1:size(features, 1), idxTrain); XTrain = features(idxTrain, :, :); YTrain = labels(idxTrain); XTest = features(idxTest, :, :); YTest = labels(idxTest);

这里需要稍微解释一下features的维度。如果样本数量是M，时间步是T，每个时间步的特征维度是D，那么features在Matlab里可以存成M×T×D的三维数组，也可以存成M×T×D的元胞。用LSTM训练时，更常见的是把每个样本转为cell中的张量，我后面会写具体做法。

4. CNN-LSTM/BILSTM网络搭建与Matlab实现

4.1 网络结构的常规接法

CNN-LSTM的经典思路，先用CNN在局部特征维度上做卷积，压缩冗余信息，再用LSTM在时间维度上建模依赖关系。在Matlab里搭建这类网络，用Deep Learning Toolbox的layer数组就能完成。

我在项目里的网络层结构如下：

numFeatures = size(features, 3); % 每个时间步的特征维度 numClasses = 4; layers = [ sequenceInputLayer(numFeatures) convolution1dLayer(3, 32, 'Padding', 'same') batchNormalizationLayer reluLayer maxPooling1dLayer(2, 'Stride', 2) lstmLayer(100, 'OutputMode', 'last') fullyConnectedLayer(numClasses) softmaxLayer classificationLayer];

需要注意，sequenceInputLayer的输入格式通常是“特征维度×时间步”的顺序，也就是每个样本是numFeatures×T的矩阵。而前面构造的features如果是M×T×D，需要用permute转成T×D或D×T。这里是最容易出维度错误的地方，建议先在代码里认真核对size。

如果只是把原始统计特征直接喂给LSTM，不用CNN也可以；但加一层1D卷积的收益是明显的。CNN能在相邻时间步之间提取局部跨特征模式，比如某个故障在有效值和峭度两个特征上同时出现跳变，CNN卷积核能更快捕捉到这种组合关系。

4.2 用LSTM还是用BILSTM

LSTM的优点是结构简单、训练速度快，适合绝大多数故障诊断场景。BILSTM则通过两个方向的LSTM同时处理序列，一个从前向后读，一个从后向前读，最后把两个方向的隐状态拼接起来，因此能同时利用每个时间步前后的上下文信息。

故障信号里的冲击特征并不总是单向的，比如一个冲击发生后，后续波形会有衰减振荡，前面的波形也会因为冲击带来的频率改变而发生变化。BILSTM在理论上比LSTM更能捕捉这种双向依赖，实际测试下来也确实是BILSTM的准确率略高一些。

但BILSTM不能无脑用。它的计算量是LSTM的两倍以上，训练时间明显变长，在小的数据集上还更容易过拟合。我个人的判断标准是：如果LSTM在验证集上的最高准确率低于90%，可以换BILSTM试试；如果LSTM已经能达到95%以上，再换BILSTM的收益很有限，性价比不高。

在Matlab里，把LSTM层换成BILSTM层只需要一行：

bilstmLayer(100, 'OutputMode', 'last')

其余结构完全一样。需要注意的是，BILSTM隐含层节点数一般取LSTM的一半也能达到差不多的效果，因为双向拼接后的输出维度已经是单向的两倍。

4.3 训练细节与超参数设置

训练网络之前，特征要先做归一化。我通常采用z-score归一化，也就是每个特征维度减去均值再除以标准差，而不是简单的min-max归一化。因为统计特征里有值域相差很大的项，比如峭度最高能到几十，而有效值只有零点几，不做z-score归一化的话，网络会把数值大的特征当成主导因素。

训练参数我建议这样设置：优化器adam，初始学习率0.001，mini-batch大小32，最大迭代次数50到80。如果损失曲线在后期震荡明显，学习率可以降到0.0005，继续训练。

另一个重要的训练技巧是加入早停机制。基于验证集的准确率监控，如果连续10轮准确率没有提升，就停止训练并恢复到最佳权重。Matlab训练选项里可以通过输出网络训练信息并手动判断，或者用ValidationFrequency参数配合验证集来做。

在实际工程里，我还习惯在训练完成后检查一下每个类别的召回率，而不只看整体准确率。比如故障诊断里，如果把内圈故障样本误判成外圈故障，虽然整体准确率可能看着还行，但实际使用中会漏掉真正的故障类型。所以分类报告比单一准确率指标更重要。

5. 实验结果对比与分析

5.1 对比方案怎么设计

为了验证MSO-VMD-CNN-LSTM/BILSTM这套方案到底有没有必要，我建议至少做四组对比。这四组我都在模拟故障信号和自己采集的样本上跑过，结果很有代表性。

• 方案A：直接CNN-LSTM，输入原始信号切片，不做任何分解；
• 方案B：固定K=5、alpha=2000的VMD + CNN-LSTM，不做参数优化；
• 方案C：MSO-VMD + CNN-LSTM，也就是用MSO优化VMD参数；
• 方案D：MSO-VMD + CNN-BILSTM。

每组用相同的训练集、测试集划分方式和训练超参数，保证公平。如果样本量比较小，可以采用多次随机划分后取平均的方式，减少偶然性。

5.2 四组实验结果怎么看

我在某设备振动信号上的实验结果表明，方案A的测试集准确率通常在80%到85%之间，而且训练过程很不稳定，个别随机种子下只能到78%。原因就是之前说的是原始信号特征太脏、故障冲击淹没在噪声里。

方案B换上VMD分解之后，准确率直接到了92%左右，说明分解本身带来的收益非常明显。但是固定参数K=5存在隐患，我用enga不同故障类型的样本分别做VMD，发现最佳K并不一样，有的类型最佳K是4，有的是7。固定K=5对部分样本会欠分解或过分解，准确率也就是“勉强能用”的水平。

方案C用MSO优化参数后，准确率进一步提升到95%以上，而且方差变小了，不同随机种子下的结果比较稳定。这说明参数寻优确实能帮VMD找到更适合当前信号的分解方式。

方案D使用BILSTM后，准确率在方案C的基础上又高了1到2个百分点，但训练时间大约是方案C的1.7倍。如果诊断任务是离线分析，建议直接用BILSTM；如果是要部署到实时监测系统里，LSTM的性价比更高。

比较有意思的是，我在换用另一组不同转速的数据时，方案B的准确率明显下降，但方案C的准确率基本没怎么变。这正好说明了用优化算法自动适配信号特征的价值，而不是靠一组固定参数打天下。

6. 常见问题与踩坑排雷

6.1 六大高频问题速查

我把实际运行中遇到最多的问题整理成了表格，方便大家快速定位。

6.2 几个真正值钱的实操经验

第一个经验，不要把适应度函数设成单一的峭度指标。峭度确实能反映冲击特征，但对随机强冲击干扰特别敏感，一个孤立噪声尖峰就能让峭度值很高，导致MSO选出一组并不适合整体数据的参数。我后来改用包络熵为主、模态混叠惩罚为辅的适应度函数，结果稳定很多。

第二个经验，优化VMD参数用的信号，最好选一段有代表性的轻度故障信号，而不是正常信号。正常信号本身包络熵偏低，优化出来的参数可能会倾向于把信号分解得过于稀疏，对后续故障识别不利。选轻度故障代表性样本做寻优，得到的参数在各类样本上泛化效果更好。

第三个经验，Matlab里跑LSTM数据格式转换是最容易烦躁的环节，没有之一。一定要记住sequenceInputLayer接收的是每个样本一个个的cell，每个cell里存的是特征维×时间步的单精度矩阵。很多人直接把三维数组塞进去，报错了半天不知道怎么回事。我的建议是先把样本数量写清楚，用一个cellfun统一完成转换，不要用循环一个个转。

第四个经验，如果信号长度非常长，不要一次性把所有样本先VMD分解再存盘，会占用大量内存。正确做法是流水线处理：先切样本，逐个对样本做VMD分解和特征提取，把提取后的特征矩阵和标签保存到本地，后面加载特征训练。这样即使原始信号有几十万个采样点，内存也能撑得住。

最后分享一个小操作。每次训练结束，不要太着急看准确率，先花几分钟打印出VMD分解后的中心频率和每个IMF的能量占比，在频谱图上一眼扫过去。如果相邻模态中心频率差小于100Hz，就要警惕过分解；如果某个IMF能量特别大且时域波形和原始信号几乎一样，说明欠分解。这个人工检查方法，比单纯调网络参数更能帮你发现问题。