KSVD-WMSDL加权多尺度字典学习在轴承故障诊断中的MATLAB实现
2026/9/15 7:35:48 网站建设 项目流程

简介:面向轴承故障检测的加权多尺度字典学习KSVD-WMSDL算法Matlab仿真资源,适合本科、研究生及教研人员开展故障诊断、特征提取与稀疏表示方向的算法验证和参数调整。内容围绕KSVD字典学习、多尺度分解与加权策略三大核心模块展开,同时附带可操作的仿真录像,方便按步骤复现,快速理解算法流程。包内共238个文件,压缩包仅3.31MB;其中95个m脚本承载主流程与实验代码,14个c文件对应omp等核心函数的mex源码,15个mexw64与15个mexmaci64为预编译动态库,可直接调用而免去手动编译,另有PDF文档、PNG结果图和avi操作录像辅助学习。资源目录组织清晰,文件类型分明,便于按需查阅和二次开发。已有276人学习下载,对于需要快速复现算法、研究字典学习机制或用于课程设计与学术对照的实验者,是一份紧凑实用的参考资料。

1. 轴承故障检测为什么绕不开字典学习

滚动轴承故障诊断里最磨人的不是分类器选型,而是特征提不出来。早期微弱故障的信号都被齿轮啮合振动和噪声盖住,直接做FFT,故障特征频率常常淹没在边带里。稀疏表示是另一条路:把信号拆成过完备字典上少数原子的线性组合,冲击成分会被分离到与它形态最接近的原子上。KSVD从训练数据里自适应学字典,比固定小波基更贴合现场信号。这个MATLAB工程把KSVD扩展成加权多尺度版本(KSVD-WMSDL),先做小波包分解再逐频带学字典,最后按故障特征频率带宽加权,用稀疏系数做识别。工程带录屏和完整OMP工具箱,适合刚接触稀疏表示又想快速复现结果的本硕学生。

2. KSVD字典更新与OMP稀疏编码的MATLAB实现

2.1 稀疏表示模型与KSVD的求解框架

设振动信号被切成m段,每段长度n,排成矩阵Y ∈ R^(n×m)。目标是用字典D ∈ R^(n×k)和稀疏系数X ∈ R^(k×m)去逼近:

min ||Y - D X||_F²,s.t. ||x_i||_0 ≤ T₀

其中||x_i||_0是非零元素个数,T₀是稀疏度上限。这个目标函数同时优化D和X,是非凸问题。KSVD的处理方式是交替迭代:固定D用OMP解稀疏系数,固定X逐列更新字典。更新某个原子d_k时,先找出所有用到该原子的样本索引,计算去掉当前原子后的误差矩阵,再对该矩阵做SVD,取最大奇异值对应的分量作为新原子和该行系数。之所以比MOD(最优方向法)稳定,是因为它更新原子时同步修正了对应行的系数,误差下降更直接。单脉冲冲击在时域上的形态恰好就是少数原子的线性叠加,这是稀疏表示能压住噪声、保留冲击特征的根本原因。

2.2 OMP每一步在做什么

OMP(正交匹配追踪)每次迭代做四件事:计算残差与所有原子的内积、挑内积绝对值最大的原子、用最小二乘更新已选原子上的系数、重算残差。重复T₀次后结束。和MP(匹配追踪)的关键区别在于,OMP每次迭代都会把所有已选原子上的系数重新投影到它们张成的子空间上,保证残差始终与已选原子正交,收敛行为更稳定。

% 单样本OMP核心逻辑(对应ompmex.c的内部流程) r = y; % 残差,初始为原始信号段 selected = []; % 已选原子索引 for iter = 1:T0 scores = abs(D' * r); % 残差在所有原子方向上的投影绝对值 [~, idx] = max(scores); selected = [selected, idx]; coeffs = D(:, selected) \ y; % 最小二乘求解已选原子上的系数 r = y - D(:, selected) * coeffs; % 更新残差 end

scores越大说明该原子与当前残差越相关,也就是信号段里最显著的成分被逐步剥离。第6行的反斜杠是MATLAB最小二乘求解,比手写正规方程数值更稳。注意这套逻辑在工具箱里被C语言mex实现并加速了,实际工程不需要用MATLAB循环去复刻,但理解它有助于排查“为什么系数里全是0”这类问题。

2.3 OMP工具箱里那些C文件的实际分工

工程压缩包里有一组C文件和MEX封装,来自KSVD工具箱的底层库,其中ompmex.c被调用最频繁,需要在运行前用mex命令编译成平台相关文件。

文件作用说明
ompmex.c批处理OMP稀疏编码入口多列信号同时求解,返回系数矩阵
omp2mex.c双稀疏OMP变体字典本身稀疏时加速
myblas.c底层BLAS相关运算矩阵乘等线性代数基础操作
collincomb.c / rowlincomb.c列/行线性组合字典更新时的组合运算
im2colstep.c / col2imstep.c信号分块与重组长信号按滑窗切成样本矩阵,供后续稀疏编码使用

ompmex.c的输入不是字典本身,而是D'*YD'*D两个矩阵。OMP每次迭代都需要残差与所有原子的内积,预计算D'*D之后,每次只需从Gram矩阵里查表,省掉重复矩阵乘的开销。im2colstep.c负责把一维振动信号按指定窗长和步长切成矩阵,相当于给字典学习准备训练样本。

2.4 MATLAB里实际调用OMP的写法

主脚本中稀疏编码部分一般是这样的:

% 参数设置 windowLen = 256; % 窗长,对应一个样本的点数 stepLen = 128; % 滑窗步长,控制样本数 sparsity = 8; % 稀疏度T0 % 滑窗切分,得到样本矩阵Y(每列一个样本段) Y = im2colstep(signal, [windowLen], [stepLen]); % 预计算互相关与Gram矩阵 A = D' * Y; % k x m,每列是字典与样本的内积 G = D' * D; % k x k,Gram矩阵,对称正定 % 调用mex版OMP求稀疏系数 X = omp(A, G, sparsity);

sparsity = 8表示每段信号最多用8个原子逼近。噪声占比高时可以加到12~15;训练样本量很大时建议保守一些。Y的列数决定样本数,滑窗步长越小样本越多,但相邻样本相关性也越强,字典学到冗余原子的概率变大。常见做法是窗长取256~512、步长取窗长的一半,跑一轮后看重构误差再微调。

3. 加权多尺度字典学习:从频带拆分到原子加权

3.1 为什么单字典处理轴承信号不够

轴承故障冲击会激起轴承座和传感器的多阶共振,内圈、外圈、滚动体故障的冲击周期不同,频谱重心也不同。单个KSVD字典在时域上学习,原子形态受最强共振分量主导,弱故障分量很容易被丢掉。这就是WMSDL多尺度结构的引入动机:先把信号按频带拆开,每个频带用KSVD单独学一套原子,再按故障特征频率所在频带的重要性加权。这样做的好处有两个。第一,频带内信号相对平稳,字典原子更纯粹,不会出现“一个原子同时拟合两个频带成分”的混叠。第二,加权落在字典合并环节,不改OMP求解器,工程实现成本很低。WMSDL里的“多尺度”落在频带上,“加权”落在字典合并环节。

3.2 小波包分解构造子带信号

多尺度字典的第一步是小波包分解。三层小波包把信号分成8个等带宽子带,每个子带对应一段频率区间。

% 三层小波包分解,db4小波基 wpt = wpdec(signal, 3, 'db4'); % 提取8个子带的重构信号 subSignals = zeros(8, length(signal)); for i = 1:8 subSignals(i, :) = wprcoef(wpt, [3, i-1]); end

小波包与普通小波分解的区别在于,它对高频细节也做二分,所以低频转频特征和高频故障冲击都能保留。wprcoef(wpt, [3, i-1])取第3层、节点i-1的重构信号,节点编号从0到7。db4是Daubechies长度4的小波,对冲击类信号波形保真度好;如果冲击衰减很快,换成sym5效果往往更好。里层数超过3层后,子带间隔变小,每个子带内的能量和样本量都下降,字典学到的原子容易过拟合到噪声上。

3.3 加权策略:把故障特征频率变成权重锚点

加权是KSVD-WMSDL区别于普通多尺度字典学习的核心。计算分三步:先对每个子带重构信号做Hilbert包络谱;再在包络谱中搜索该子带对应的轴承故障特征频率(BPFI、BPFO、BSF)及其倍频;最后把包络谱中落在特征频率邻域内的能量与子带总能量的比值作为该子带的权重。

% 子带权重计算:以内圈故障BPFI为例 fs = 12000; % 采样率 BPFI = 236.4; % 内圈故障特征频率,单位Hz w = zeros(8, 1); for i = 1:8 env = abs(hilbert(subSignals(i, :))); % 包络 spec = abs(fft(env)); % 包络谱 f = (0:length(spec)-1) / length(spec) * fs; band = f > (BPFI-3) & f < (BPFI+3); % 主频±3Hz邻域 w(i) = sum(spec(band)) / sum(spec(f > 20)); % 去除直流后的能量占比 end % 归一化并作为子字典加权系数 w = w / sum(w);

权重取值范围0到1且和为1,保证加权字典不改变整体稀疏表示的数量级。f > 20是为了滤掉包络谱中接近直流的低频干扰。工程上一般把这个权重更新过程放进迭代循环,每轮用新字典重新解码训练信号再更新权重,让字典逐步聚焦到故障相关频带。权重迭代3到5轮后基本收敛,再多跑只会让高权重子带持续膨胀。

3.4 子字典加权拼接与OMP复用

得到权重后,把8个子字典按权重缩放后拼接成一个冗余字典。拼接时不要求子字典互相正交,也不需要对原子做额外处理,OMP会自适应选择有用的原子。

% 假设 DictL{i} 是第i个子带的KSVD字典 wDict = []; for i = 1:8 wDict = [wDict, w(i) .* DictL{i}]; end % 用加权字典做稀疏编码,接口与单字典一致 G = wDict' * wDict; A = wDict' * testSignalMat; X = omp(A, G, sparsity);

把权重乘到字典上,等价于在目标函数里给对应原子加了惩罚项:权重小的原子内积贡献被压缩,OMP在选择原子时天然避开它们。这样做的最大好处是不用改mex文件,直接复用omp(A, G, sparsity)。小波包分解的参数推荐如下表,数据量不足时优先降层数而不是降窗长。

参数推荐值依据
分解层数38个子带覆盖共振频带,数据量小时不选4层
小波基db4 / sym5对冲击波形保真度好
子字典原子数16~32兼顾泛化能力与字典冗余度
权重迭代轮数3~5超过5轮高频子带权重会接近1

4. 轴承故障检测完整仿真流程与分类实验

4.1 从原始信号到诊断结果的整体管线

仿真主脚本结构分为五段:数据读取、小波包分解、子字典学习、加权融合、稀疏表示分类。数据来自公开轴承数据集或实验室采集的振动信号,采样率常见12k或48k。整体流程如下:

% 主流程伪代码 signals = loadBearingData('fault_inner.mat'); % 载入内圈故障信号 trainY = im2colstep(signals, [256], [128]); % 训练样本 DictL = cell(8, 1); % 1) 多尺度分解 wpt = wpdec(signals, 3, 'db4'); for i = 1:8 subSig = wprcoef(wpt, [3, i-1]); subY = im2colstep(subSig, [256], [128]); % 2) 子带KSVD字典学习 DictL{i} = ksvd(subY, atomNum, 'maxIter', 20, 'T', sparsity); end % 3) 加权融合 wDict = buildWeightedDictionary(DictL, w, fs, BPFI); % 4) 稀疏编码与特征提取 X = omp(wDict' * testY, wDict' * wDict, sparsity); feat = extractFeatures(X);

每个环节的尺寸要先对齐:DictL{i}大小为256×atomNum,wDict大小为256×(8·atomNum)。8个子带各学atomNum个原子会造成字典冗余,atomNum一般取16~32,8个字典合计128~256个原子。这里用到matlab统计和机器学习工具箱里的函数,2021a版本直接内置。

4.2 数据预处理与样本切分要点

预处理直接影响字典训练。高频噪声分量如果混进训练样本,字典会把噪声当成有意义结构学进去,所以先做带通滤波保留共振频带,再滑窗切分。窗长要与冲击周期匹配:窗太短,一个窗口装不下完整冲击衰减过程;窗太长,两个相邻冲击落进同一窗口,OMP会混乱。

% 带通滤波示例(巴特沃斯二阶带通) fc1 = 500; % 高通截止 fc2 = 4000; % 低通截止 [b, a] = butter(2, [fc1, fc2]/(fs/2), 'bandpass'); sigF = filtfilt(b, a, rawSignal); % 滑窗切分 step = round(256 * 0.5); % 窗长一半作为步长 trainY = im2colstep(sigF, [256], [step]);

这里用filtfilt而不是filter,因为前者做零相位滤波,不引入相位偏移,冲击位置在时间轴上不漂移。滤波截止频率不要卡得太死,轴承座共振频带通常在几百到几千赫兹,保留500~4000Hz基本覆盖大多数情况。如果字典学出的原子全是正弦状波形,说明训练样本频带太窄,滤波器带宽需要放宽。

4.3 稀疏系数特征与分类判别

稀疏系数矩阵X不能直接丢给分类器,因为它的行数等于原子数。工程里常用的特征有四类:

特征计算方式物理含义
系数l2范数sqrt(sum(X.^2,1))信号在该字典下的表示能量
重构误差sum((Y-D*X).^2,1)字典与信号的匹配程度
原子索引分布sum(X~=0,2)哪些子带的原子被频繁使用
实际稀疏度sum(X~=0,1)非零系数个数
% 特征构造 reconErr = sum((testY - wDict * X).^2, 1); % 每个样本的重构误差 coefL2 = sqrt(sum(X.^2, 1)); % 系数l2范数 usedIdx = sum(X ~= 0, 2); % 每个原子被使用的次数 % 特征表送SVM featTable = [reconErr; coefL2]'; svmModel = fitcsvm(featTable(1:100, :), labels(1:100), 'KernelFunction', 'rbf');

usedIdx放在特征表里往往有奇效:正常轴承信号在加权字典上的原子分布是散开的,故障信号的原子会集中在故障子带对应的索引区间。如果分类准确率上不去,先看这一列特征的分布有没有区分度。RBF核适合这种特征维度不高但类别边界不规则的场景。

4.4 运行录屏里复现实验的正确顺序

工程附带的AVI录屏里,从mex编译到出图全流程都有操作。建议按三遍走:第一遍看完整运行过程,第二遍自己运行主脚本,第三遍改参数对比。MATLAB 2021a直接运行主脚本,第一次跑会编译mex文件,需要装好支持的编译器。若报错找不到ompmex,在命令行执行mex -setup选择gcc或MinGW。

% 验证字典是否工作正常:重构前后波形对比 reconSignal = col2imstep(wDict * X, size(sigF), [256], [step]); rsnr = 10 * log10(sum(sigF.^2) / sum((sigF - reconSignal).^2)); fprintf('Reconstruction SNR = %.2f dB\n', rsnr);

重构信噪比在8~20dB之间都算正常。低于6dB说明字典没学好,要回头检查滑窗步长和迭代次数;高于25dB则要怀疑字典过拟合到噪声上了。工程里还有个常见做法是把正常样本和故障样本的usedIdx画成柱状图,看原子使用分布的重叠程度,比直接看混淆矩阵更直观。

5. 参数匹配关系与字典学坏的预判方法

5.1 字典尺寸、稀疏度与权重的匹配关系

训练样本数m、原子数k、稀疏度T₀三者的关系比很多资料里说的严苛。原子总数超过样本数1/3时,字典很容易过拟合,每个样本几乎都有专属原子。我一般用下表的经验范围:

参数建议区间说明
窗长256~512由冲击衰减时间决定,观察两个脉冲间距
原子数/子带16~32总原子数不超过样本数的1/3
稀疏度T₀6~12随噪声水平增大而增大
权重迭代轮数3~5过多会让少数子带权重接近1,失去多尺度意义
小波包层数3~4数据量小时用3层,4层需要更多训练样本

权重迭代轮数是最容易被忽略的参数。每轮全流程结束,用新字典重新解码训练信号,更新一次权重;跑到第4、5轮权重就稳定了,再跑只会让高权重子带持续膨胀。常见做法是固定跑4轮,取第4轮的权重作为最终结果。这本质上是个不动点迭代,4轮足够收敛到工程可用的精度。

5.2 字典学坏的三种典型信号

第一种:字典里大量原子彼此高度相似。计算原子互相关矩阵G,如果非对角元素均值超过0.8,说明字典冗余,原子数设多了。减原子数比加稀疏度更有效,因为稀疏度增加只会让每个样本用更多原子,不解决原子间冗余。

第二种:重构误差随迭代下降极慢甚至上升。检查训练样本中是否有几个异常大冲击把字典带偏。解决方案是把超过3倍中位数的样本段剔除,用干净样本重训。轴承升降速阶段的数据经常触发这个问题,截取平稳转速段再训练。

第三种:权重迭代后某个子带权重超过0.9。这看起来像聚焦效果很好,实际上其他子带信息被全部丢弃,多尺度退化成单字典。可以把该子带的原子数减少,把其他子带的权重下限固定为0.05,强制保留频带多样性。拿你自己的信号跑一遍,如果RSNR稳定在10dB以上,且正常与故障样本的原子索引分布有明显分界,这套KSVD-WMSDL流程就算在数据上站稳了。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询