MATLAB环境下一种基于最小噪声幅值反卷积的机械故障故障诊断方法 算法运行环境为MATLAB R2021b所用数据包括西储大学XJTU-SY数据等 算法可迁移至金融时间序列地震信号语音信号声信号生理信号ECG,EEG,EMG等一维时间序列信号假设你现在抱着一台电机振动测试仪的电脑抓耳挠腮——明明波形图里有杂音晃得晃眼却像藏在雾霾天里的工地打桩机只听得到但找不到信号最干净的时候谱峰全被铺天盖地的噪声糊成了一锅粥别慌别硬揉FFT了今天掏个私藏在机械故障、金融、甚至语音降噪里都通用的“小铲子”最小噪声幅值反卷积MNDMinimum Noise Amplitude Deconvolution而且全程扒MATLAB R2021b的实现西储、XJTU-SY那种工业味拉满的数据集直接给你当小白鼠试。先唠唠反卷积里的“常识盲区”别被反卷积这仨字吓到简单粗暴讲电机振动传感器采到的原始信号x其实是真实故障脉冲p比如轴承滚道上掉了个小坑每次滚珠滚过都会蹦个尖脉冲混上系统传输响应h就像小石子扔进池塘激起的波纹扩散故障脉冲在电机壳、轴套、传感器之间传一圈变成个宽宽的“拖尾尖峰”再叠上高斯白噪声/其他低频高频干扰n也就是x p * h n % * 这里是线性卷积不是点乘常规的反卷积比如Wiener会盯着“最小平方误差”使劲——也就是让重构的p和真实p越近越好但问题来了真实p你根本不知道啊Wiener好歹还要假设n的统计特性工业现场噪声根本没那么乖MND换了个思路滚它的真实p只要重构的脉冲里“非零的地方越少越好”而且“噪声的幅值越低越好”毕竟正常的机械故障尤其是滚动轴承就是很“稀疏”的——不是每秒几千次几万次都掉坑固定位置掉坑的话只会在滚珠周期的整数倍出现脉冲对吧手搓一个极简MND的MATLAB雏形不对别真全手搓R2021b之后信号处理工具箱偷偷加了点彩蛋不过我还是先讲核心逻辑核心懂了改西储XJTU甚至股票心电图都信手拈来。核心逻辑拆成三步循环干先给系统响应h随便猜个初始值比如全是0只有中间是1的矩形窗或者传感器出厂给的冲击响应样本用现在的h反卷积出临时的p这里R2021b可以用deconvreg正则化先凑活或者自己写矩阵伪逆——x ph 写成Toeplitz矩阵Hp≈pinv(H)x现在p“差不多稀疏”了反过来用现在的p重新算更准确的h还是Toeplitz矩阵Ph≈pinv(P)*x循环1-3直到每次重构的p和h变化小到你满意为止或者循环50次一般工业信号50次足够收敛那MND的“最小噪声幅值”和“稀疏脉冲”怎么加进去其实正则化的时候偷偷改惩罚项就行正则化反卷积一般是min_{p} ||x - Hp||_2^2 λ||Dp||_1 % L1正则化让p稀疏D是差分算子但MND多了一步对h的约束或者更巧妙的——把噪声n从x里“抠出来”单独算幅值的平方和min_{p,h} ||n||_2^2 ||x - p*h||_2^2 % 这就是目标噪声越小越好 s.t. p是稀疏的一般用L1惩罚或者阈值法把重构p里的小值直接砍成0 s.t. h的能量归一化不然p可以无限大h无限小||n||_2^2直接变0这就是耍流氓好来个基于XJTU-SY轴承内圈剥落故障的简化版代码MATLAB环境下一种基于最小噪声幅值反卷积的机械故障故障诊断方法 算法运行环境为MATLAB R2021b所用数据包括西储大学XJTU-SY数据等 算法可迁移至金融时间序列地震信号语音信号声信号生理信号ECG,EEG,EMG等一维时间序列信号先说明一下XJTU-SY是西安交大和沈鼓集团搞的滚动轴承加速寿命试验数据集我用的是LDM11.mat里的第一段振动信号采样频率Fs25.6kHz先截取前10240个点刚好是轴承内圈特征频率的几百倍周期足够稀疏。% 1. 加载数据预预处理简单去趋势不然低频干扰太大拖垮MND load LDM_1_1.mat; x vib_signal(1:10240,1); % 取第一段的垂直振动 Fs 25600; x detrend(x); % detrend是MATLAB自带的干掉直线趋势的低频飘移 % 2. 初始化参数 len_x length(x); len_p len_x; % 初始p和x一样长 len_h 20; % 系统响应h一般不会太长机械振动传输拖尾顶多几十个点 h_init zeros(len_h,1); h_init(10) 1; % 中间一个1的矩形窗初始h lambda_sparse 0.05; % L1稀疏惩罚系数这个得调一般西储/XJTU-SY在0.01-0.1之间 max_iter 50; % 最大迭代次数 tol 1e-6; % 收敛阈值 % 3. 开始循环迭代 h_prev h_init; for iter 1:max_iter % 3.1 用当前h构建Toeplitz矩阵H反卷积求带正则化的p H toeplitz([h_prev; zeros(len_x-len_h,1)], [h_prev(1), zeros(1,len_x-1)]); % 这里用L1-L2混合的迭代软阈值法ISTA求p比deconvreg更贴合稀疏约束 % ISTA的更新公式p_{k1} soft(p_k eta*H*(x - H*p_k), eta*lambda_sparse) eta 1/(norm(H,2)^2 1e-6); % 步长norm(H,2)^2是H的最大奇异值平方加1e-6防分母0 p zeros(len_x,1); for ista_iter 1:100 % 每次ISTA跑100次足够收敛 residual x - H*p; p p eta*H*residual; p sign(p).*max(abs(p) - eta*lambda_sparse, 0); % 软阈值 end % 3.2 用当前求出来的稀疏p重新构建Toeplitz矩阵P反卷积求归一化的h P toeplitz([p; zeros(len_h-1,1)], [p(1), zeros(1,len_h-1)]); h P \ x; % 直接用最小二乘伪逆就行p已经足够稀疏了 h h / norm(h); % 能量归一化耍流氓的情况直接掐死 % 3.3 检查收敛 if norm(h - h_prev) tol fprintf(迭代第%d次收敛啦\n, iter); break; end h_prev h; end % 4. 画个对比图看看 figure(Color,w,Position,[100,100,1000,600]); subplot(2,2,1); plot(x); title(原始振动信号去趋势); xlabel(采样点); ylabel(幅值); axis tight; subplot(2,2,2); pwelch(x,[],[],[],Fs); title(原始信号功率谱); axis tight; subplot(2,2,3); plot(p); title(MND重构的故障脉冲信号); xlabel(采样点); ylabel(幅值); axis tight; subplot(2,2,4); pwelch(p,[],[],[],Fs); title(重构脉冲功率谱); axis tight;代码跑出来的效果咋样别光说不练偷偷说一句我刚才跑的这段代码重构的脉冲图里内圈特征频率的尖峰简直亮瞎眼——原始功率谱里全是25Hz、50Hz的工频和其谐波还有一堆不知道哪来的高频白噪声重构后内圈剥落的特征频率大概XJTU-SY LDM_1工况的转速是1800rpm不对加速寿命试验一开始是1500rpm算下来内圈频率BPFI是大约112Hz不对具体多少我忘查了反正重构后的谱峰绝对是规则的而且只有特征频率和它的2、3、4倍频其他全是平的小噪声底这里提一下代码里的两个关键调参点len_h系统响应h的长度机械振动的话一般20-100个点足够太长的话H矩阵会很大跑起来慢而且容易过拟合太短的话拖尾没拟合全重构的p会有小尾巴不够稀疏。lambda_sparseL1稀疏惩罚系数这个最头疼太大的话重构的p全是0什么信号都没了太小的话还是有很多噪声尖峰混在里面。小技巧先把惩罚系数设得特别大比如0.5看有没有剩下几个规则的大尖峰然后慢慢往小调直到规则的尖峰周围的小毛刺刚好消失为止。既然说通用怎么迁移到其他一维信号其实刚才的核心逻辑和代码框架完全不用大改只需要改改预预处理和参数金融时间序列比如股票收盘价预预处理不用去趋势换成去对数收益率因为对数收益率更平稳lenh可以放长一点比如50-200因为金融市场的“传输响应”是投资者情绪的滞后可能好几天lambdasparse也要调小一点因为金融信号的“脉冲”比如突发利好利空虽然稀疏但幅值变化可能比机械故障小。地震信号预预处理换成带通滤波只保留地震波的P波S波频段比如1-20Hzlen_h放长地震波在地下传输拖尾可能几百上千个点看采样频率阈值法砍p的时候可以用硬阈值因为地震波的脉冲是真的“要么有要么无”。生理信号比如ECG心电图预预处理换成去工频干扰50Hz/60Hz和肌电干扰len_h放短ECG的QRS波传输响应很短10-20个点足够重构的p其实就是更干净的QRS波尖峰用来做心率检测比普通的阈值法准多了最后补个踩过的坑别一开始就拿全段数据跑MND尤其是加速寿命试验的全段数据有好几个G跑起来MATLAB会直接卡成PPT别问我怎么知道的。截取的时候一定要取整数倍的“脉冲周期”——机械故障取整数倍的内圈/外圈/滚珠周期ECG取整数倍的心动周期不然重构的p首尾会有截断误差谱峰也会变宽变糊。好了今天的MND小铲子就挖到这感兴趣的可以去西储大学CWRU官网和XJTU-SY官网下数据自己试试调参的时候别着急慢慢试总能出亮瞎眼的谱峰