十年匠心定制 · 商业建站与技术教学双线并行 咨询热线:400-886-1026 service@lmnt.cn
ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

ITD分解原理与MATLAB实现:滚动轴承故障诊断的振动信号分析

ITD分解原理与MATLAB实现:滚动轴承故障诊断的振动信号分析 简介面向机械设备故障诊断场景的ITD信号分解Matlab实现包针对复杂机械信号中故障特征难提取、信号成分难区分的问题提供完整的分解算法基础代码。压缩包共5个文件全部为.m源码脚本整体仅4KB内含主分解函数Itd.m、基础分解ItdBaseDecomp.m、极值点提取extr.m与extrema_1.m以及可直接运行的示例example01.m函数模块划分明确从极值检测到基础分解再到顶层调用构成一条完整学习链。ITD基于瞬时时间差分思想可把原始信号拆分为不同动态特征分量Widevcb相关处理思路也体现在分解流程中适用于振动信号分析、异常状态识别、早期故障预警等场景运行example01.m可快速观察分解结果便于理解分量重构与故障特征提取的具体流程。已有452人学习下载适合信号处理与机械故障诊断方向的研究者、工程师快速复现算法并移植到自己的实验数据中开展进一步分析。1. 为什么故障信号要分解ITD 在故障诊断里的位置机械设备运行时的振动信号永远不是干净的。齿轮箱里同时转着好几根轴滚动轴承的故障脉冲会被结构和传感器耦合调制再叠加上随机噪声时域里看起来就是一大段没有规律的波形。直接整段做 FFT 只能看到能量占主导的转频与啮合频率早期磨损产生的冲击特征常被淹没在 20 dB 量级的背景里。ITD 分解Intrinsic Time-Scale Decomposition固有时间尺度分解正好解决这个问题把一段非平稳、非线性的原始信号拆成若干固有旋转分量PRC每个 PRC 对应一段相对集中的瞬时频带故障冲击会优先集中到最靠前的高频分量中从而把微弱特征从背景里分离出来。这个 .zip 源码包以 itd、ItdBaseDecomp、extr、extrema_1 和 example01 五个文件给出了一个能在 MATLAB 里直接跑的完整流程适合做滚动轴承与齿轮箱振动监测、故障信号分解和特征提取的工程师。2. ITD 分解原理与源码包结构Itd.m 与 ItdBaseDecomp.m 在做什么2.1 ITD 与 EMD 的本质区别ITD 与经验模态分解EMD的目标相似都是把复杂信号拆成按频率从高到低排列的分量区别在“怎么拆”。EMD 需要先用三次样条分别包络极大值点和极小值点取上下包络均值作为低频分量然后反复迭代直到包络均值趋近于零这个过程要解线性方程组还要处理过冲问题。ITD 则把这个步骤简化成找到局部极值点在相邻极值点之间做线性插值直接得到一条基线 L(t)再让 X(t)−L(t) 成为第一个固有旋转分量。基线相当于把当前信号里的低频趋势先提走剩下的就是围绕基线振荡的那部分。这样设计带来的直接好处是计算量明显下降。样条插值在数据长度上万点时会不断被调用而 ITD 只用一次线性插值就能完成一轮基线提取迭代次数取决于信号里极值点的密度工程上常见一维振动信号分解到第 5~8 层基本就只剩单调趋势了。更重要的是线性插值不会像样条那样在突变附近产生过冲冲击型的早期故障特征在 PRC 里保留得更完整。如果拿同一段轴承外圈故障信号分别跑 EMD 和 ITD能看到 ITD 的第一层分量里冲击脉冲的幅值起伏更贴近原始信号模式混叠的程度也轻于 EMD。对于故障信号分解的场景还有一个容易被忽略的差异ITD 分量的频带重叠比 EMD 小稳定性更强。分解出的 PRC 之间频率区间相对清楚后续不管算包络谱还是时域指标都不用反复在多个 IMF 里拼凑一个完整的故障频带。2.2 包内文件职责源码包一共五个文件入口和辅助函数分得很清楚脚本名在分解链路中的位置主要职责Itd.m入口函数参数检查和列向量化组装输出结构后调用 ItdBaseDecomp.mItdBaseDecomp.m迭代分解核心提取极值、构造基线、分离 PRC并控制迭代终止extr.m极值辅助扫描差分符号变化返回局部极大与极小值的索引extrema_1.m极值校正处理连续相等采样点和端点返回可直接用于插值的极值坐标example01.m示例脚本构造仿真故障信号调用 Itd.m 并绘制原始信号与各 PRC这里要注意 extr.m 与 extrema_1.m 不是重复代码。extr.m 通常返回原始极值索引直接拿去做线性插值会在连续重复采样点位置出现阶梯状基线extrema_1.m 会先做去重和端点修正输出给 interp1 的坐标才是严格单调递增的。实际使用中如果发现分解出的第一个 PRC 首尾有异常大的尖峰优先查这一步要么是端点没有外插要么是连续极值没有合并。2.3 ItdBaseDecomp 的核心循环与停止条件拿到源码后我先读的是 ItdBaseDecomp.m 的循环体它几乎就是 ITD 算法的文字定义function [PRCs, residual] itd_base_decomp(x) % x一维振动信号列向量长度 N % PRCs固有旋转分量每列一个分量按分解顺序从高频到低频排列 x x(:); N length(x); t (1:N).; PRCs []; residual x; maxIter 20; for k 1:maxIter ext extrema_1(residual); % 第一步提取局部极值 if length(ext) 3 break; % 剩余部分已单调停止分解 end base interp1(ext, residual(ext), t, linear, extrap); PRCs(:, k) residual - base; % 第二步原信号减基线得到 PRC residual base; % 第三步基线进入下一层迭代 end逻辑上分成三步先用 extrema_1 找到极值点再对极值序列做分段线性插值得到基线最后把 residual−base 作为当前层分量、base 作为新的待分解信号继续循环。interp1里的extrap很关键它保证插值结果覆盖到首尾两个采样点如果去掉基线两端会出现 NaN后面所有分量都会受影响。参数上值得调整的有两个maxIter控制最大分解层数经验上取min(floor(log2(N)), 20)就够采样率 1024、1 秒数据对应 20 层上限实际一般在 8 层前就触发极值不足的停止条件length(ext) 3是停止条件的底线少于 3 个极值点时线性插值只有一段无法再构造有效基线。如果你处理的是冲击十分密集的高速轴承信号可以把层数上限改为按N / 极值点数动态判断避免最后一个分量里残留明显周期成分。3. 在 example01.m 上跑通故障信号分解从波形到 PRC3.1 example01.m 里常见的仿真故障信号这套包里 example01.m 的用途是把“故障信号”先画出来再调用 Itd.m 走一遍完整流程。故障信号分解的验证必须先有可信的信号源工程里常用仿真信号代替实测因为实测信号里的真实故障特征频率很难事先确定。示例里常见配置是采样率 fs2048转频 fr25Hz外圈故障特征频率 fo 约 96.3Hz故障冲击激发出 1200Hz 左右的系统共振再叠加白噪声fs 2048; t (0:fs-1) / fs; fr 25; % 轴转频 fo 96.3; % 轴承外圈故障特征频率 fn 1200; % 系统共振频率 zeta 0.06; % 阻尼比 impact exp(-zeta*2*pi*fn*mod(t, 1/fo)) .* sin(2*pi*fn*t); sig 0.6*sin(2*pi*fr*t) 0.8*impact 0.3*randn(size(t));这样构造出来的数据转频幅值占主导冲击序列被 mod(t, 1/fo) 按故障周期触发每次冲击都会以 1200Hz 共振频率衰减振荡白噪声的加入让信号更接近现场手感。理解这段构造很重要因为后面验证 ITD 分解是否有效靠的就是确认 PRC 里能不能还原出 96.3Hz 的冲击重复频率。3.2 运行 Itd 并绘制前几个 PRC 分量把 example01.m 里的信号换成本节上面这段数据然后在 MATLAB 命令行执行[PRCs, residual] itd_base_decomp(sig); figure; for k 1:3 subplot(4, 1, k1); plot(t, PRCs(:, k), LineWidth, 0.6); ylabel(sprintf(PRC%d, k)); xlim([0 0.2]); end第一行调用的是前面封装好的迭代分解函数。绘图部分只取了前 0.2 秒因为冲击脉冲集中在这个时间段内能直接看出波形形态。正常情况下 PRC1 里能看到间隔约 10.4ms 的冲击衰减振荡PRC2 主要包含共振后的低频衰减PRC3 则更接近转频正弦。如果某个分量幅值明显大于原始信号多半是最高一层没截断残差把直流趋势和低频转频混在一起需要回头调整停止条件再分割一次。3.3 用频谱验证分解结果不是“过分解”ITD 分解和 EMD 一样存在过分解风险不能只凭时域波形好看就下结论。我一般会补一张 FFT 频谱确认各 PRC 的能量集中在不同频带Nfft 2048; f_ax (0:Nfft-1) / Nfft * fs; S1 abs(fft(PRCs(:, 1), Nfft)); S2 abs(fft(PRCs(:, 2), Nfft)); S1 S1(1:Nfft/2); S2 S2(1:Nfft/2); f_ax f_ax(1:Nfft/2); cf1 sum(S1.^2 .* f_ax) / sum(S1.^2); % 第一个 PRC 的频谱质心 cf2 sum(S2.^2 .* f_ax) / sum(S2.^2); % 第二个 PRC 的频谱质心质心频率是信号能量集中位置的加权平均比直接取峰值更稳。判断标准有三条相邻分量的主频带质心至少错开一个带宽每个 PRC 的能量占原始信号总能量比不小于 1%故障特征频率所在分量在包络谱里有明显主峰。三条全过才可以认为分解合理。如果只做第一层分解、不检查频带重叠就可能出现“把同一个故障边带拆进两个 PRC”的情况后面的诊断结论会被带偏。4. ITD 分解 Widevcb故障特征频带选择与包络解调4.1 Widevcb 在故障信号分解流程里的位置ITD 分解本身只负责“拆分”不负责“判定”哪个分量含故障。对滚动轴承来说故障冲击会激发系统共振而共振频带通常比故障特征频率本身高很多。ITD 分解出来的前几个 PRC 正好覆盖这段共振区但到底用哪个分量做包络解调手动去选非常麻烦。Widevcb 在这套流程里承担的就是“宽频带候选选择”自动从各 PRC 的功率谱里找到能量集中峰把峰值附近的半功率区间划成一个宽带候选带再返回给后续包络分析使用。它的名字可理解为 Wide Vibration Characteristic Band 的缩写对应宽振动特征带。传统共振解调需要事先知道共振频率一旦设备结构改变、转速变化固定带通滤波器就会失配Widevcb 让频带跟随 ITD 分解的结果动态生成避开手动定带通参数的脆弱环节。4.2 Widevcb 参数与调用方式不同版本源码里参数名可能略有差异但思路一致核心参数是带宽候选数量和半功率门限参数含义经验取值fs采样率与采集设置一致2048~10240num_peaks每个 PRC 功率谱上取的峰值个数3~5band_edge_dB频带边界相对峰值的衰减门限-3dB 或 -6dBmin_bandwidth最小带宽避免候选带过窄fs/64调用时通常直接把 PRC 矩阵传进去[band_list, peak_f] widevcb(PRCs, fs, ... num_peaks, 4, band_edge_dB, -3, min_bandwidth, fs/64);函数内部先对每个 PRC 做功率谱估计再按半功率门限搜索局部峰值最后输出一个[低边界, 中心频率, 高边界]的列表。band_edge_dB-3表示只保留从峰值下降不超过 3dB 的频段这样框出来的带通常能覆盖共振峰的主要旁瓣如果发现包络谱里噪声过大就把门限放到 -6dB带宽变大代价是故障特征频率附近的谱线会更平缓。4.3 包络谱验证 Widevcb 选的频带是否正确选定候选带后下一步是对目标 PRC 做 Hilbert 包络再对包络做 FFT 得到包络谱env abs(hilbert(PRCs(:, 2))); % 取 PRC2 的包络 env_spec abs(fft(env, Nfft)); f_env (0:Nfft-1) / Nfft * fs; mask f_env 20 f_env 200; % 只看故障特征频段 [~, idx] max(env_spec(mask)); f_fault f_env(mask); f_detect f_fault(idx); % 解调出的故障特征频率这里的f_detect就是解调出来的故障特征频率对照 96.3Hz 看误差是否在 1Hz 以内。如果包络谱主峰落在别处先检查 Widevcb 返回的带宽是不是把共振峰切开了再看 PRC 序号是否选择过低。故障冲击通常集中在前两三个分量里选到第四层之后大多是转频谐波。工程上常见误用是直接对原始信号做包络谱结果转频谐波能量太大故障特征频率全被盖住。ITDWidevcb 的价值就是把分解、选带、解调三段串成一条可自动化的链路。5. 量化故障程度用 PRC 峭度与能量占比做状态估计5.1 两个指标的计算口径分解完后把 PRC 变成诊断结论我常用两个量峭度和能量占比。K kurtosis(PRCs, 1); % 列方向总体标准差 E sum(PRCs.^2, 1) / sum(sig.^2); % 每个 PRC 的能量占比峭度反映信号中冲击成分占比正常振动接近高斯分布 K≈3故障早期 K 能到 5~8能量占比看故障能量集中在第几层选最高的前两层做分析即可其余分量归入噪声。5.2 经验判断表指标参考表现诊断含义PRC1 峭度K4 且明显大于 PRC2、PRC3冲击能量集中在最高频层典型轴承早期磨损前两层能量占比75%分解充分故障信号占主导噪声被推到后续层包络谱故障频率幅值大于全频带均值 3 倍以上可确认该特征频率对应故障建议进入趋势跟踪多组数据对比时固定 ITD 的迭代上限和 Widevcb 的 band_edge_dB 参数峭度的趋势才有可比性。参数一改峰峰值和带宽都变横向对比就失真了。5.3 提高峭度指标稳定性的调整实测信号里偶发噪声会把 PRC1 的峭度抬得很高但持续一两秒就消失。常见做法是把 1 秒信号切成 8 个 0.125 秒的窗分别算峭度后取中位数再把 Widevcb 选出的频带在这 8 个窗里都跟踪一遍中位数超过阈值才报警。这样比单次整段计算的误报率低很多也顺便给出了故障程度的波动范围。本文还有配套的精品资源点击获取
返回列表