多彩编程 多彩编程MZPH · CODE BLOG
ARTICLE DETAIL

文章详情

深耕前端与后端开发技术的一线实战笔记与踩坑复盘。

SSA-VMD参数自动寻优与皮尔逊模态筛选的振动信号降噪方法

SSA-VMD参数自动寻优与皮尔逊模态筛选的振动信号降噪方法 简介面向信号处理方向的学生、科研人员与工程技术人员这份资源提供了一套将麻雀搜索算法SSA与变分模态分解VMD结合的完整MATLAB实现方案用于解决信号分解与降噪问题。其核心思路是先用SSA自动寻优VMD的模态数K与惩罚因子α再对分解得到的IMF分量计算皮尔逊相关系数据此筛选噪声分量并施加小波软/硬阈值降噪最后完成信号重构与效果评估同时配套IMF分量图、频谱图、Hilbert边际谱及3D视图等可视化输出。资源包共14个文件以11个m脚本为主体辅以mat数据、xlsx与csv结果文件整体约190KB结构紧凑、开箱即用。目前已有149人学习适合希望快速复现SSA-VMD联合降噪流程、理解参数优化与分量筛选逻辑的读者参考借鉴。1. SSA-VMD 这套组合拳到底在解决什么信号难题一段从轴承、齿轮箱或者超声探头采回来的振动信号十有八九是糊的。工频干扰、电磁噪声、机械冲击混在一起直接做包络谱或者做故障特征提取谱线要么被噪声淹没要么被模态混叠搅成一团。变分模态分解VMD就是冲着这个场景来的它把信号拆成若干个中心频率互不重叠的带限模态比经验模态分解EMD那种递归筛分更抗模态混叠端点效应也轻得多。但 VMD 有个绕不开的坎——分解层数 K 和惩罚因子 alpha 得人来定K 选大了过分解选小了欠分解alpha 选偏了模态带宽要么糊要么碎。麻雀搜索算法SSA就是来替你把这两个参数自动搜出来的搜完再用皮尔逊系数挑出真正含故障信息的模态最后小波阈值降噪加信号重构把干净的特征波形交到你手上。这套流程适合做旋转机械故障诊断、超声检测、EEG/EMG 预处理的工程师MATLAB 完整代码和数据能让你当天就跑通自己的数据。2. 为什么用 SSA 而不是网格搜索来定 VMD 参数2.1 VMD 的两个参数为什么不能拍脑袋VMD 的本质是求解一个约束变分问题把输入信号 f 分解成 K 个模态分量 u_k每个模态围绕自己的中心频率 omega_k 震荡约束是所有模态加起来能重构回原信号。它的核心迭代在频域完成靠交替方向乘子法ADMM更新 u_k 和 omega_k。这里 K 决定了模态数量alpha 决定了每个模态的带宽惩罚力度。K 的敏感性很好理解一个含内圈故障的轴承信号理论上故障冲击会调制出若干边频带K 取 3 可能只拆出主频和两个边带K 取 8 就会把噪声也拆成独立模态出现中心频率挨得很近的虚假分量。alpha 的敏感性更隐蔽alpha 小模态带宽宽容易把相邻频率成分吸进同一个模态alpha 大带宽窄模态容易碎成尖峰重构误差反而上升。工程上常见的做法是固定 alpha 扫 K或者固定 K 扫 alpha但这两个参数是耦合的单变量扫描经常错过真正的最优组合。网格搜索的问题在于代价。K 取 2 到 10alpha 取 500 到 5000 步长 500就是 9 乘 10 等于 90 次 VMD每次 VMD 在几千点的信号上跑几十毫秒到几百毫秒加上你要算适应度整体时间还能忍。但如果信号长、K 范围宽、alpha 精度要求高网格搜索的维度灾难就出来了。更麻烦的是网格搜索没有记忆它不会根据已经试过的点调整下一步往哪走。2.2 SSA 的寻优逻辑和适应度函数怎么定麻雀搜索算法是 2020 年前后提出来的一种群智能优化算法模仿麻雀种群的觅食和反捕食行为。种群分成发现者和跟随者发现者负责找食物丰富区域跟随者跟着发现者走另外随机选一部分麻雀当警戒者感知到危险就带着种群往安全区跳。它的更新公式里有两个关键随机数一个控制发现者的搜索范围一个控制警戒者的逃逸方向整体收敛速度比粒子群PSO快陷入局部最优的概率也低一些。用 SSA 优化 VMD核心是适应度函数怎么定。常见的有三类包络熵最小、样本熵最小、重构误差最小。做故障诊断我一般用包络熵因为故障冲击越明显包络越稀疏包络熵越小。包络熵的定义是对每个模态做希尔伯特变换取包络归一化后算信息熵K 个模态的包络熵取平均或者取最小。如果你做的是降噪重构用重构误差加模态相关性的组合更合适。% SSA 优化 VMD 参数的主循环骨架 % 输入signal 一维信号dim 优化维度K 和 alphalb/ub 下上界 % 输出best_K, best_alpha, convergence_curve function [best_K, best_alpha, curve] ssa_vmd_optimize(signal, dim, lb, ub, maxIter, popNum) % 初始化种群位置第一维是 K整数第二维是 alpha pop repmat(lb, popNum, 1) rand(popNum, dim) .* repmat(ub - lb, popNum, 1); fitness zeros(popNum, 1); for i 1:popNum fitness(i) vmd_fitness(signal, round(pop(i,1)), pop(i,2)); end [bestFit, idx] min(fitness); bestPos pop(idx, :); curve zeros(maxIter, 1); for t 1:maxIter [~, sortIdx] sort(fitness); pop pop(sortIdx, :); fitness fitness(sortIdx); % 发现者占 20%位置更新偏向当前最优 PD round(popNum * 0.2); for i 1:PD if rand 0.5 pop(i,:) pop(i,:) .* exp(-i / (rand * maxIter eps)); else pop(i,:) pop(i,:) randn(1, dim); end end % 跟随者跟随发现者警戒者随机跳 for i PD1:popNum if i popNum / 2 pop(i,:) randn(1, dim) .* exp((pop(popNum,:) - pop(i,:)) / i^2); else pop(i,:) pop(1,:) abs(pop(i,:) - pop(1,:)) .* randn(1, dim); end end % 边界处理K 取整并限制在 [2, 10] pop(:,1) min(max(round(pop(:,1)), 2), 10); pop(:,2) min(max(pop(:,2), lb(2)), ub(2)); for i 1:popNum fitness(i) vmd_fitness(signal, pop(i,1), pop(i,2)); end [curBest, idx] min(fitness); if curBest bestFit bestFit curBest; bestPos pop(idx, :); end curve(t) bestFit; end best_K round(bestPos(1)); best_alpha bestPos(2); end这段代码里pop的第一列是 K第二列是 alphavmd_fitness是你自己写的适应度函数内部调用 VMD 并返回包络熵。PD控制发现者比例0.2 是常用值种群大可以调到 0.3。maxIter一般取 20 到 50再大收敛曲线就平了。注意 K 必须取整否则 VMD 会报错边界处理那两行不能省。2.3 适应度函数和 VMD 调用的具体写法适应度函数决定了 SSA 往哪个方向搜写错了整个优化就是白跑。下面这个版本用最小包络熵适合故障冲击明显的信号。function fit vmd_fitness(signal, K, alpha) % 调用 VMD返回 K 个模态 % 常见做法是设置 DC0, init1, tol1e-7 try [u, ~, ~] VMD(signal, alpha, 0, K, 0, 1, 1e-7); catch fit 1e10; % VMD 失败给一个极大惩罚值 return; end entropySum 0; for k 1:K env abs(hilbert(u(k, :))); env env / (sum(env) eps); pk env(env 0); entropySum entropySum - sum(pk .* log(pk eps)); end fit entropySum / K; endVMD函数的参数顺序在不同版本里可能不一样常见的是VMD(signal, alpha, tau, K, DC, init, tol)tau 取 0 表示无噪声容忍DC 取 0 表示不强制第一个模态为直流init 取 1 表示均匀初始化中心频率。hilbert做希尔伯特变换取解析信号包络归一化后算熵。如果 VMD 因为 K 太大报错catch 里给 1e10 让 SSA 自动避开这个区域。这个适应度函数越小越好SSA 内部用min找最优逻辑一致。3. 皮尔逊系数筛模态和小波阈值降噪怎么接3.1 用皮尔逊系数挑出含故障信息的模态VMD 分解完你手上是 K 个模态但真正含故障特征的往往只有一两个剩下的要么是工频、要么是噪声。皮尔逊相关系数衡量两个变量的线性相关程度取值在 -1 到 1 之间绝对值越接近 1 越相关。把每个模态和原信号算皮尔逊系数系数大的模态和原信号共享更多成分通常就是主频和故障调制成分所在的模态。但这里有个坑如果原信号本身噪声很大噪声模态和原信号的皮尔逊系数也可能不低。所以常见做法是结合峭度一起看皮尔逊系数高且峭度大的模态优先保留。峭度对冲击成分敏感故障冲击会让峭度明显偏离 3。我一般设一个阈值皮尔逊系数绝对值大于 0.3 且峭度大于 3 的模态进入重构其余丢弃。% 皮尔逊系数筛选模态 corrCoef zeros(1, K); kurtVal zeros(1, K); for k 1:K % corrcoef 返回 2x2 矩阵取非对角元素 cc corrcoef(signal, u(k, :)); corrCoef(k) cc(1, 2); kurtVal(k) kurtosis(u(k, :)); end % 筛选条件相关系数绝对值大于 0.3 且峭度大于 3 selectedIdx find(abs(corrCoef) 0.3 kurtVal 3); if isempty(selectedIdx) [~, selectedIdx] maxk(abs(corrCoef), 2); % 兜底取前两个 end signal_recon sum(u(selectedIdx, :), 1);corrcoef返回的是相关系数矩阵cc(1,2)才是原信号和模态的相关系数。kurtosis默认算的是超额峭度还是原始峭度取决于 MATLAB 版本老版本返回的是原始峭度正态分布为 3新版本默认减 3用之前先help kurtosis确认一下。maxk是取最大的几个兜底逻辑防止筛选条件太严导致一个模态都不剩。3.2 小波阈值降噪的参数怎么设重构信号里还有残余噪声小波阈值降噪是收尾的常规操作。它的逻辑是对信号做多层小波分解高频细节系数里噪声占主导用阈值把小于阈值的系数置零或收缩再逆变换回来。关键参数有三个小波基、分解层数、阈值规则。小波基选db4或sym8是工程上的常见做法db4支撑长度短、计算快sym8对称性好、重构失真小。分解层数一般取floor(log2(N))再减 1 到 2N 是信号长度层数太多会把有用高频也压掉。阈值规则里sqtwolog是固定阈值rigrsure是无偏似然估计heursure是前两者结合minimaxi是极小极大准则。振动信号我一般用rigrsure或heursure阈值函数用软阈值硬阈值容易在重构时产生附加震荡。% 小波阈值降噪 wname db4; level floor(log2(length(signal_recon))) - 1; % 使用启发式阈值和软阈值 [thr, sorh, keepapp] ddencmp(den, wv, signal_recon); signal_denoised wdencmp(gbl, signal_recon, wname, level, thr, sorh, keepapp);ddencmp自动给出阈值、软硬阈值选择和是否保留近似系数den表示降噪wv表示小波。wdencmp的gbl表示全局阈值每个细节层用同一个阈值。如果你想让每层阈值不同把gbl换成lvd并传入每层阈值向量。降噪完的信号就是最终输出可以直接做包络谱或者特征提取。3.3 重构和效果验证的对照方法重构就是把筛选后的模态加起来再经过小波降噪。验证效果不能只看波形好不好看得用指标说话。常用的有信噪比SNR、均方根误差RMSE、平滑度指标。如果你有干净信号做参考SNR 和 RMSE 直接算如果没有就看包络谱里故障特征频率及其倍频是否清晰突出。% 效果对比原始信号 vs 重构降噪后 snr_val 10 * log10(sum(signal.^2) / sum((signal - signal_denoised).^2)); rmse_val sqrt(mean((signal - signal_denoised).^2)); % 包络谱对比 env_orig abs(hilbert(signal)); env_deno abs(hilbert(signal_denoised)); f (0:length(signal)-1) * fs / length(signal); figure; subplot(2,1,1); plot(f, abs(fft(env_orig - mean(env_orig)))); subplot(2,1,2); plot(f, abs(fft(env_deno - mean(env_deno))));fs是采样频率包络谱里故障特征频率fault_freq及其 2 倍频、3 倍频如果比原始信号更突出说明这套流程有效。注意 SNR 计算需要干净参考信号实际工程里往往没有那就退而求其次看包络谱峰值信噪比或者看峭度提升量。4. 避坑与排查这套流程最容易翻车的五个地方4.1 VMD 报错矩阵维度不一致或直接卡死现象SSA 跑到一半 MATLAB 报错或者 VMD 内部迭代不收敛程序卡住。原因K 被 SSA 更新成了非整数或者超出信号允许范围或者 alpha 太小导致模态中心频率初始化冲突。另外VMD函数版本不同参数顺序不一样传错位置也会报维度错误。解决在 SSA 边界处理里强制round取整并限制 K 在 2 到 10 之间alpha 限制在 500 到 5000。调用 VMD 前先which VMD确认函数来源用help VMD核对参数顺序。如果还卡把tol从 1e-7 放宽到 1e-6迭代上限设 500。4.2 皮尔逊系数筛选后模态全被丢掉现象selectedIdx是空数组重构信号全是零或者报错。原因阈值设得太严比如相关系数要求大于 0.5 且峭度大于 5实际信号里没有一个模态同时满足。解决加兜底逻辑筛选为空时按相关系数绝对值排序取前两个。或者把阈值降到 0.2 和 2.5。更稳的做法是先用相关系数排个序画出来看看拐点在哪再定阈值。4.3 小波降噪后故障冲击被削平现象降噪后波形光滑了但包络谱里故障特征频率反而变弱。原因分解层数太多或者用了硬阈值把故障冲击对应的高频细节系数也压掉了。解决层数减 1 到 2阈值函数换成软阈值或者改用rigrsure这种自适应阈值。如果还不行跳过小波降噪直接对重构信号做包络谱很多时候 VMD 加皮尔逊筛选已经够干净了。4.4 SSA 收敛曲线震荡不下降现象curve画出来上下跳最后的最优值还不如手动设的 K5、alpha2000。原因种群太小或者发现者比例太低搜索不充分或者适应度函数里 VMD 失败返回的惩罚值太大把种群多样性压死了。解决种群从 20 加到 30maxIter从 20 加到 30发现者比例从 0.2 调到 0.3。惩罚值从 1e10 降到 1e6避免个别失败点主导排序。另外检查适应度函数里有没有对模态做归一化没归一化的话熵计算会受幅值影响。4.5 换一组数据整套参数就失效现象在轴承数据上跑得好好的 K 和 alpha换到齿轮箱数据上分解结果一塌糊涂。原因SSA 每次都是针对当前信号重新寻优的如果你把上一次的最优参数硬编码进去等于没优化。另外不同信号的采样率、长度、噪声水平不同适应度函数的尺度也不一样。解决每次换数据都重新跑 SSA不要复用参数。如果嫌慢把 SSA 种群和迭代次数降下来先粗搜再精搜。适应度函数里对包络熵做归一化让不同信号之间的适应度值可比。5. 把 SSA-VMD 用到你自己的数据上三个调参习惯第一个习惯是先看信号再定搜索范围。拿到一段新信号先画时域波形和频谱估一下主频带在哪、噪声水平多高。如果主频集中在 1kHz 以下K 的上限设 8 就够了如果信号本身很干净alpha 下限可以提到 1000避免过宽的模态。我一般会先用 K5、alpha2000 手动跑一次 VMD看看分解出来的模态长什么样再决定 SSA 的搜索边界往哪偏。第二个习惯是适应度函数别只用一种指标。包络熵对冲击敏感但对连续振荡信号不敏感样本熵对复杂度敏感但计算慢重构误差直接反映分解质量但容易偏向 K 大的解。我通常把包络熵和重构误差加权组合权重 0.7 和 0.3这样既能突出冲击成分又不会让 K 无限增大。加权公式写在vmd_fitness里改一行就行。第三个习惯是留一组对照。跑完 SSA-VMD 加皮尔逊筛选加小波降噪之后别急着下结论把原始信号、VMD 重构信号、最终降噪信号三条包络谱画在一张图上对比。如果最终结果的故障特征频率峰值比原始信号高出 3dB 以上这套流程就算有效如果只高出 1dB 不到检查一下是不是皮尔逊筛选阈值太松把噪声模态也重构进去了。这个对照习惯帮我省了很多次返工也让我在换数据时心里有底。我自己的教训是一开始总想用一套参数打天下结果在齿轮箱数据上翻车了好几次后来老老实实每次重新寻优反而稳定了。这套 SSA-VMD 加皮尔逊系数加小波阈值降噪的流程代码量不大但参数之间的耦合关系得亲手调过才有感觉。希望帮到你。本文还有配套的精品资源点击获取
返回列表