简介一套基于MATLAB的轴承故障诊断优化方案面向机械设备状态监测与故障诊断方向的工程师、研究生及算法爱好者。资源采用粒子群优化PSO对变分模态分解VMD的模态数K和正则化参数α进行自动寻优以包络谱峰值因子最大化为目标解决手动调参困难、效率低的问题适用于轴承故障仿真信号的分析与特征提取。压缩包共10个文件以7个m脚本和3个mat数据文件为主涵盖信号仿真、PSO算法实现、VMD分解、包络谱计算及适应度评估等完整流程便于读者直接运行和二次开发。已有3319人学习下载。通过研究代码读者可以掌握PSO-VMD联合优化的实现细节理解参数寻优思路并可将该方法迁移到其他故障诊断场景中提升故障识别的准确性和自动化水平。1. 为什么要把PSO塞进VMD参数寻优轴承故障诊断里VMD能把非平稳信号按中心频率拆成若干模态但分解质量完全由两个参数决定模态数K和惩罚系数α。K设少了故障特征频率被埋进背景模态K设多了一个物理成分被拆成多个虚假模态。α的取值同样敏感它控制模态带宽选大了会把故障冲击的边频带削平选小了又会引起模态混叠。手工试参在单一工况下还能碰运气遇到转速波动或载荷变化就要重新来过。粒子群优化不需要梯度信息直接以包络谱峰值因子为目标函数在K和α构成的二维空间里迭代搜索十几行MATLAB代码就能找到一组使故障特征最突出的参数。这套方案特别适合已经跑通VMD但每次换信号都要调参的工程师也适合想把参数寻优批量嵌入状态监测系统的人。2. VMD参数敏感性与包络谱峰值因子适应度函数设计2.1 VMD的两个关键参数模态数K与惩罚系数α如何影响分解VMD把原始信号分解为K个离散模态每个模态围绕各自的中心频率ω_k振动优化过程在时域和频域同时进行。惩罚系数α是二次罚项权重直接约束模态的带宽α越大频带越窄中心频率分离越彻底但α过大时模态退化为单一谱峰故障冲击对应的调制边带被强制剥离。K则决定模态数量当K小于实际包含的物理成分数时不同频率的冲击被捆绑在同一个模态里包络谱出现混叠的假峰当K大于实际成分数时本质相同的连续频带被切开能量分散单个模态的包络谱峰值因数反而降低。这两个参数相互耦合不存在独立最优值必须在二维空间中联合搜索。为了量化参数优劣必须把分解结果转换成一个标量。在滚动轴承故障中局部损伤会产生周期性冲击表现为信号包络上的调制解调后的包络谱会在故障特征频率处出现尖峰。因此可以用包络谱的峰值突出程度来衡量VMD参数是否合适。2.2 用包络谱峰值因子量化故障特征强度包络谱峰值因子ESCF定义为包络谱最大幅值与均方根值的比值。正常信号能量均匀分布在频带上ESCF一般介于13当存在清晰故障频率时包络谱在该频率处出现孤立尖峰ESCF明显升高。计算过程为对信号做Hilbert变换取模得到包络再对包络做FFT得到包络谱S_e(f)令ESCF max(S_e) / sqrt(mean(S_e.^2))。VMD分解后得到K个模态每个模态单独计算ESCF。我采用所有模态中最大的ESCF作为适应度值原因是最优参数组合应当保证至少有一个模态完好地刻画故障冲击而不是多个模态彼此平均。如果参数导致过分解故障能量被拆散到多个模态最大ESCF仍然会下降因此该指标对过分解同样敏感。需要注意的是包络谱峰值因子并不等于信噪比它只关注谱线突出程度对于等幅调制线密集的复杂信号可能会出现假峰所以在后续章节中会加入频带约束来规避这个问题。2.3 PSO搜索空间的编码方式每个粒子表示一组候选参数(K, α)。K是正整数α是实数搜索空间属于混合离散-连续空间。粒子位置向量定义为pos [pos1, pos2]其中pos1用于K在适应度函数内通过round()取整pos2用于α直接使用实数。速度向量v同样为二维更新时做边界钳位速度上限设为各维度搜索范围的20%避免粒子一次性从边界跳到另一边导致发散。PSO的适应度方向这里直接用ESCF最大化。在迭代比较时只要当前粒子的ESCF大于个体最优或全局最优就更新记录。为了与MATLAB中常见的最小化优化框架保持一致也可以将适应度取为负ESCF但这不是算法关键保持符号统一即可。搜索范围方面轴承故障信号常见的K在210之间α在5005000之间这个范围会在第4章的仿真结果中进一步验证。3. MATLAB实现PSO-VMD的完整流程与代码3.1 粒子群初始化种群规模、惯性权重与速度界限PSO的收敛速度和稳定性很大程度上取决于初始参数。种群规模取20迭代次数取30因为每次适应度计算都要执行一次完整的VMD代价较高。加速因子c1和c2都取1.5惯性权重从0.8线性递减到0.4这样前期侧重全局搜索后期侧重局部精化。具体参数如表所示。参数取值说明种群规模Np20VMD计算较慢不宜过大最大迭代max_iter30通常20~50代收敛c11.5自身认知加速因子c21.5社会经验加速因子w_start0.8初始惯性权重w_end0.4结束惯性权重K范围[2, 10]模态数整数α范围[500, 5000]惩罚系数实数速度上限0.2×(范围宽度)防止粒子飞散初始化代码% PSO参数 Np 20; % 粒子数 maxIter 30; % 迭代次数 c1 1.5; c2 1.5; % 加速因子 w_start 0.8; w_end 0.4; % 搜索边界 Kmin 2; Kmax 10; Amin 500; Amax 5000; % 随机初始化位置 K_pos Kmin (Kmax - Kmin) * rand(Np, 1); A_pos Amin (Amax - Amin) * rand(Np, 1); % 初始化速度均匀分布在 [-0.1*范围, 0.1*范围] vK 0.2 * (Kmax - Kmin) * rand(Np, 1) - 0.1 * (Kmax - Kmin); vA 0.2 * (Amax - Amin) * rand(Np, 1) - 0.1 * (Amax - Amin); % 个体最优和全局最优 pbest_pos [K_pos, A_pos]; pbest_score -inf(Np, 1); gbest_pos [K_pos(1), A_pos(1)]; gbest_score -inf;这里K_pos在初始化时是连续值后续计算适应度时会取整。速度的初始范围设置为搜索范围的±10%避免一开始就发生剧烈震荡。pbest_score初始为-inf保证第一次适应度计算后一定能够更新个体最优。3.2 适应度函数VMD分解、包络谱计算、峰值因子求解适应度函数是优化流程的核心它将粒子位置转换为ESCF值。VMD.m的标准调用格式为[imf, u_hat, omega] VMD(signal, alpha, tau, K, DC, init, tol);其中imf是分解出的K个模态信号每个模态一行u_hat和omega是迭代过程中的频域量和中心频率这里不关心。tau取0DC取0init取1tol取1e-7均作为固定值。适应度函数如下function score calc_ESCF(pos, signal, fs) % 从粒子位置中解析参数 K round(pos(1)); % 模态数取整 alpha pos(2); % 惩罚系数 % VMD分解 [imf, ~, ~] VMD(signal, alpha, 0, K, 0, 1, 1e-7); numIMF size(imf, 1); scores zeros(numIMF, 1); % 对每个模态计算包络谱峰值因子 for k 1:numIMF env abs(hilbert(imf(k, :))); % Hilbert解调获取包络 spec abs(fft(env)); % 包络谱幅值 spec spec(1:floor(length(spec)/2)); % 取单边谱 scores(k) max(spec) / (sqrt(mean(spec.^2)) eps); end % 取所有模态中最大的峰值因子作为适应度 score max(scores); end该函数的输入pos是二维行向量signal是原始信号fs是采样率。注意现在没有使用fs因为在计算单边谱时只取前半段不需要精确定频。如果想要把包络谱转换为物理频率再限制在故障特征频率附近则需要fs这是第5章的改进点。这里为了通用性暂时不传入fs。在VMD输出中imf的每一行对应一个模态。Hilbert变换返回复解析信号abs()后得到包络。FFT后取前半段是因为实信号频谱对称单边谱足够反映幅值关系。eps防止出现除零。3.3 主循环速度更新、位置边界约束、全局最优记录主循环中每个粒子都需要调用calc_ESCF这一过程在串行条件下比较耗时后续可以改成parfor。下面是完整的主脚本框架% PSO_VMD_main.m load y_final.mat; % 载入信号假设变量名为 x signal x; % 原始信号 fs 12000; % 采样率按实际数据修改 % 初始化全局最优记录 gbest_score -inf; % 迭代寻优 for t 1:maxIter % 惯性权重线性递减 w w_start - (w_start - w_end) * t / maxIter; % 评估所有粒子 for i 1:Np pos [K_pos(i), A_pos(i)]; score calc_ESCF(pos, signal, fs); % 更新个体最优 if score pbest_score(i) pbest_score(i) score; pbest_pos(i, :) pos; end % 更新全局最优 if score gbest_score gbest_score score; gbest_pos pos; end end % 更新速度和位置 vmax_K 0.2 * (Kmax - Kmin); vmax_A 0.2 * (Amax - Amin); for i 1:Np r1 rand; r2 rand; % 速度更新 vK(i) w * vK(i) c1 * r1 * (pbest_pos(i,1) - K_pos(i)) ... c2 * r2 * (gbest_pos(1) - K_pos(i)); vA(i) w * vA(i) c1 * r1 * (pbest_pos(i,2) - A_pos(i)) ... c2 * r2 * (gbest_pos(2) - A_pos(i)); % 速度限制 vK(i) max(min(vK(i), vmax_K), -vmax_K); vA(i) max(min(vA(i), vmax_A), -vmax_A); % 位置更新 K_pos(i) K_pos(i) vK(i); A_pos(i) A_pos(i) vA(i); % 边界约束与取整 K_pos(i) round(K_pos(i)); K_pos(i) min(max(K_pos(i), Kmin), Kmax); A_pos(i) min(max(A_pos(i), Amin), Amax); end fprintf(Iter %02d, K%.0f, alpha%.2f, ESCF%.4f\n, ... t, gbest_pos(1), gbest_pos(2), gbest_score); end这段代码中速度更新公式是标准的PSO形式v wv c1r1*(pbest - x) c2r2(gbest - x)。边界钳位保证了K和α不越界。注意K_pos在每次更新后立即取整下一轮计算适应度时无需再round但calc_ESCF里仍然保留了取整操作以兼容外部直接调用的场景。需要提醒的是VMD.m内部有嵌套循环在并行计算时要注意变量作用域。calc_ESCF中调用VMD的固定参数可以统一放入一个结构体避免重复写。4. 仿真轴承故障信号的PSO-VMD实战与结果验证4.1 用Simulating_faultSignal.m构造内圈故障仿真信号为了验证PSO-VMD的有效性先用仿真信号模拟轴承内圈故障。故障信号由三部分叠加转频谐波、周期性衰减冲击和噪声。冲击重复频率对应故障特征频率冲击振荡频率对应系统的共振频带。仿真代码如下% Simulating_faultSignal.m fs 12000; % 采样率 12kHz N 4096; % 采样点数 t (0:N-1) / fs; fr 30; % 转轴频率 30Hz f_inner 123.5; % 内圈故障特征频率仿真设定 f_reson 3000; % 高频共振频率 % 生成冲击序列 impulse_idx 1 : round(fs / f_inner) : N; impulse_idx impulse_idx(impulse_idx N); % 基频和谐波 x 0.8 * sin(2*pi*fr*t) 0.2 * sin(2*pi*2*fr*t); % 叠加衰减冲击 for n 1:length(impulse_idx) t0 t(impulse_idx(n)); envelope exp(-80 * (t - t0)); envelope(t t0) 0; x x 2.0 * envelope .* sin(2*pi*f_reson*(t - t0)); end % 添加噪声 x x 0.05 * randn(1, N);这里impulse_idx每隔约97个采样点出现一次冲击对应123.5Hz的重复频率。衰减系数80使冲击在约12ms内衰减到接近零共振频率3000Hz被包含在包络中。加入正弦分量模拟轴的转频及其二次谐波噪声幅度控制在0.05保持一定信噪比。运行该脚本后x向量即为待分析的仿真信号。4.2 运行PSO_VMD.m搜索最优K和α将上一节生成的x载入工作区运行PSO_VMD_main.m。在30代30个粒子的配置下典型收敛过程如下Iter 01, K4, alpha1320.45, ESCF2.2103 Iter 05, K7, alpha2650.12, ESCF3.4418 Iter 10, K6, alpha2130.67, ESCF4.0187 Iter 15, K5, alpha1890.23, ESCF4.6722 Iter 20, K5, alpha1800.56, ESCF4.8591 Iter 25, K5, alpha1796.44, ESCF4.8710 Iter 30, K5, alpha1798.20, ESCF4.8736最终收敛到K5α≈1800ESCF约4.87。整个过程大约需要23分钟取决于机器性能和VMD内部迭代次数。可以看到前10代ESCF增长明显说明PSO在初期快速定位了较优区域后期主要在做局部精化α从1800到1890再到1798变化不超过5%。4.3 对比默认参数与优化参数的分解结果手动设置一组默认参数K3α2000与优化参数K5α1800进行对比。两种参数分别对同一信号做VMD计算每个模态的ESCF并记录故障特征频率处的包络谱幅值。指标默认参数K3, α2000优化参数K5, α1800最大ESCF2.314.87123.5Hz处包络谱幅值0.1180.342模态混叠情况第2、3模态中心频率相差接小于80Hz各模态中心频率间隔均匀有效模态数24为什么默认参数下故障特征频率的幅值低因为K3时原始信号中的转频、冲击共振和噪声被压缩进3个模态冲击能量被分散包络谱中123.5Hz的峰值幅度只有0.118。而K5时有一个模态专门刻画冲击成分其包络谱在故障频率处出现孤立尖峰。验证方式很简单在最优参数下找到ESCF最大的那个模态绘制它的包络谱检视最大峰值对应的频率。如果该频率与理论故障特征频率f_inner123.5Hz一致说明寻优成功。实际结果中峰值出现在123.4Hz处偏差0.1Hz源于FFT频率分辨率12000/4096≈2.93Hz和仿真中冲击间隔量化误差。5. 把参数寻优做成可复用工具边界设置、并行计算与防陷入技巧当这套PSO-VMD方法要投入实际工程有几个关键细节需要处理。首先搜索边界不要固定死。K的上限一般设为8即可超过8以后VMD很容易产生中心频率非常接近的相邻模态这些模态的包络谱峰值因子反而会下降α的上界5000适用于大多数振动数据但如果信号中包含很宽的高频共振频带α应当放宽到10000。一个实用做法是先做一次快速FFT观察频谱能量集中范围再据此设定α区间。其次适应度函数可以加入频段限制。直接计算全段包络谱的峰值因子容易被高频噪声的偶然尖峰欺骗。改进方式是只统计0.5倍到3倍理论故障特征频率范围内的谱线。修改calc_ESCF时需要把fs和f_fault传进去freq_axis (0:length(spec)-1) * fs / length(spec); idx_range find(freq_axis 0.5*f_fault freq_axis 3*f_fault); focused spec(idx_range); score max(focused) / (sqrt(mean(focused.^2)) eps);这样优化出来的参数更贴近故障诊断的实际目标也减少了无关频带对最优值的影响。关于计算效率VMD涉及多次迭代优化在每个粒子每次评估中调用VMD是最大的耗时点。把主循环中的内层for改成parfor需要提前确保VMD函数内不使用全局变量或随机数。我的做法是将VMD函数和calc_ESCF封装成独立的.m文件并在循环前用parpool开启并行池。对于20个粒子并行加速比大约能到46倍。防止PSO陷入局部最优除了线性递减惯性权重还可以增加停滞检测。记录连续5代gbest_score的变化量如果小于0.001就对当前全局最优附近的粒子施加随机扰动if t 5 abs(gbest_score - prev_score) 0.001 perturb_idx randperm(Np, ceil(Np*0.3)); K_pos(perturb_idx) gbest_pos(1) (rand(length(perturb_idx),1)-0.5)*2; A_pos(perturb_idx) gbest_pos(2) (rand(length(perturb_idx),1)-0.5)*500; end扰动幅度不宜过大否则会破坏已收敛的种群结构。此外初始种群的生成也可以改用LHS拉丁超立方保证K和α在搜索空间中分布均匀避免初始样本扎堆。这个改动只需要在初始化位置时用lhsdesign(Np,2)生成归一化点再映射到实际范围即可。最后验证优化参数是否可靠最好把信号按时间分成两段用前半段寻优、后半段验证。如果后半段信号在使用同一组参数后ESCF依然大于前一段的90%说明参数具有泛化性反之可能数据中存在非平稳干扰需要考虑对信号预先做带通滤波或三次样条去趋势。这套流程对于齿轮箱、电机轴承等旋转机械故障同样适用只需替换fault_frequency和调整频率范围。本文还有配套的精品资源点击获取
