☰
ITD信号分解MATLAB工具包:旋转机械故障诊断实战指南
2026/10/31 10:27:38 网站建设 项目流程

简介:面向机械设备故障诊断与信号处理研究者的ITD分解算法实现,适用于旋转机械、滚动轴承等场景的振动信号分析。代码基于瞬时时间差分思路,将故障信号拆分为不同动态特性的分量,以分离正常运行基频、磨损引起的高频冲击及噪声成分,帮助定位轴承、齿轮等部件的早期异常。压缩包共5个文件,均为Matlab源文件,整体仅4KB,包含主分解脚本、基础分解模块、极值点提取函数及可直接运行的示例程序,分别承担核心迭代、底层差分运算、局部特征定位与流程演示,结构清晰、接口精简,便于对照修改和二次开发。已有452人学习,适合信号分解方向初学者快速上手,也适合工程师用于故障特征提取验证。通过运行示例可观察各分量输出,结合Widevcb宽频带处理思路,能更深入地理解非平稳信号的分离机制,为实际设备状态监测与故障诊断提供可复用的代码基础。

1. ITD 信号分解:一份能直接跑通的故障信号分解 MATLAB 工具包

做设备振动诊断的人,大概率翻车过 EMD:分解慢、模态混叠、端点飞翼,样本一长就卡。ITD(Intrinsic Time-Scale Decomposition,固有时间尺度分解)是另一种拆解思路——它不迭代求包络,而是通过极值点构造基线信号,一次分解就能得到旋转分量和基线,速度比 EMD 快一个量级。这份 ITD 分解工具包就是干这个事的:里面有 Itd.m 主函数、基线分解函数、极值提取函数和一个能直接跑的 example01.m 示例脚本。适合正在做轴承、齿轮箱等旋转机械故障诊断的工程师,或者刚接触信号分解想找一条能尽快复现路径的研究生。核心价值一句话:把原始振动信号拆成可解释的旋转分量,故障特征不用靠玄学。

2. ITD 分解原理与包内代码职责:先从「基线 + 旋转分量」看起

2.1 ITD 与 EMD 的分家之处:为什么用极值线性基线

先说清楚 ITD 到底在做什么。给定一段离散振动信号 x(t),ITD 把它拆成两部分:

x(t) = L(t) + H(t)

其中 L(t) 是基线信号,代表信号的局部均值趋势;H(t) 是旋转分量,代表围绕基线的振荡部分。和 EMD 拆出若干个 IMF 不同的是,ITD 的拆分是严格的一步到位——一次运算就得到一组基线和一个旋转分量。如果对残余基线继续做同样的操作,就能像剥洋葱一样得到多层旋转分量,这就是多层 ITD 分解。

基线怎么构造?按 ITD 原始论文的做法,先找信号的所有极值点(极大值和极小值),记极值时刻为 τ_k、极值大小为 X(τ_k)。然后在相邻两个极值之间定义一个控制点:

L(τ_{k+1}) = α · X(τ_k) + (1 − α) · X(τ_{k+1})

α 默认取 0.5,意思是基线控制点落在两个相邻极值的中间位置,既不超前也不滞后。得到这些控制点之后,再在控制点之间做分段线性插值,就得到了整条基线 L(t)。最后 H(t) = x(t) − L(t)。

这个设计和 EMD 最大的区别在于:EMD 每一层都要用三次样条包络迭代多次,直到满足停止条件;ITD 只做一次极值查找和一次线性插值,计算量和信号长度严格线性相关。我实际对比过,同样一段 10 万点的轴承振动数据,EMD 跑到第 4 层就要十几秒,ITD 跑满 60 层也就一两秒。所以只要不是对分解精度有极端要求,实时监测场景里 ITD 的性价比明显更高。

那这个包里的文件是怎么对应的?从文件名能猜出大概:

文件名职责关键点
Itd.m主入口函数通常负责调度整段分解流程,返回多层分量或直接出图
ItdBaseDecomp.m基线分解函数对应上述 α 控制点与插值逻辑,输入一段信号,输出基线与旋转分量
extr.m极值查找返回极值位置和极值大小,供基线控制点使用
extrema_1.m异名极值辅助函数处理端点外延和相邻极值去重,避免基线在端点处发散
example01.m示例脚本生成或读取一组信号,调用主函数跑一遍完整分解并绘图

如果你打开 ItdBaseDecomp.m,大概率看到的流程就是:x = x(:) 转成列向量 → 调用 extrema_1.m 找极值 → 按 α 公式算基线控制点 → 线性插值得到 L → H = x − L。这一步是整份资源的命根子,后面所有特征提取都建立在它产出的旋转分量上。

2.2 Itd.m 与 ItdBaseDecomp.m 的分工:包内文件职责对照

很多第一次拿到压缩包的人会困惑:到底该调 Itd.m 还是 ItdBaseDecomp.m?我的习惯是把它俩当成入口和内核的关系。

ItdBaseDecomp.m 是单层内核,只负责拆一次。你给它一个列向量,它返回 [L, H]。要拆多层,你得自己写循环:把上一轮的分量再喂进去。而 Itd.m 大概率是封装好的完整入口,内部循环迭代、控制停止条件、甚至直接出图。example01.m 里多半就是调 Itd.m。

拆开压缩包后我建议按这个顺序做三件事:先打开 example01.m 看它调了哪个函数、传了什么参数;再打开 Itd.m 看它的循环体里是用 residual 作为下一层输入,还是用 baseline 作为下一层输入;最后打开 ItdBaseDecomp.m 确认 α 和插值方式写死在代码里还是可配置。

这里有个值得留意的细节:ITD 的层间输入有两种流派。一种把每层的旋转分量 H 继续往下拆,得到的是平行的细节分量;另一种把每层的基线 L 继续往下拆,得到的是一组从高频到低频逐渐过渡的旋转分量。故障诊断里我更推荐后一种,因为低频段的旋转分量往往对应转频及边带,高频段对应冲击特征,物理意义更直观。example01.m 如果画出了从上到下频率依次降低的多层波形,那它用的应该就是基线递归。

2.3 用 example01.m 跑通第一次分解:代码骨架与参数解释

拿到包以后不要急着改代码,先把示例跑通。我一般会在 MATLAB 里把 example01.m 拆成下面这个骨架来理解:

% example01.m 的典型骨架:构造一段含故障特征的信号并执行 ITD 分解 clear; clc; close all; fs = 1000; % 采样率,单位 Hz t = (0:999) / fs; % 1 秒信号,共 1000 点 f_base = 50; % 基频分量,模拟设备转频 f_fault = 137; % 故障特征频率,比如某个轴承部件的通过频率 x = sin(2*pi*f_base*t) + 0.35*sin(2*pi*f_fault*t) + 0.08*randn(size(t)); % 调用 ITD 基线分解,第二参数是最大分解层数 [baseline, rotation] = ItdBaseDecomp(x, 60); figure; subplot(2, 1, 1); plot(t, x); title('原始信号'); xlabel('时间/s'); subplot(2, 1, 2); plot(t, rotation); title('第一层旋转分量'); xlabel('时间/s');

这段代码的逻辑很简单:先合成一个仿真信号,基频 50 Hz、故障特征频率 137 Hz 加少量噪声,然后调 ItdBaseDecomp.m 得到基线和旋转分量,最后画图对比。

参数层面的注意点有两个。第一个是 fs 必须和真实采集设备一致,示例里 fs=1000 只是为了跑通流程,换成实际数据后如果不改采样率,后面所有频谱分析都会错位。第二个是 ItdBaseDecomp 第二个参数 60 的语义——它到底是「最多迭代 60 层」还是「只拆 60 个点」?我建议你打开函数体确认一下,按文件命名习惯推断大概率是分解层数。如果不确定,可以先用小参数 5 跑一遍,看返回数组的行数和 x 是否一致。跑通之后,下一步就是把这个分解过程接到真实故障诊断流程里。

3. 搭一个故障诊断流程:从振动波形到可解释的故障特征

3.1 预处理:去直流、去漂移,让分解吃到干净信号

ITD 对输入信号的平稳性比 EMD 宽容,但宽容不等于可以乱喂。实测中最常遇到的是两类情况:传感器输出的信号带直流偏置,或者低频的趋势项把基线整体抬走。无论是哪种,分解出来的第一层旋转分量都会被污染。

我处理真实采集数据时一般按这个顺序来:

% 假设 x 是读入的原始振动数据,fs 是采样率 x = x(:); % 强制转为列向量 fs = 25600; % 按采集设备实际标定值修改 % 第一步:去直流 x = x - mean(x); % 第二步:去除低频趋势项,用二次多项式拟合 t = (0:numel(x)-1)' / fs; p = polyfit(t, x, 2); x_detrend = x - polyval(p, t); % 第三步:可选,带通滤波,只保留关注频段 % [b, a] = butter(4, [500, 5000] / (fs / 2), 'bandpass'); % x_filtered = filtfilt(b, a, x_detrend);

第一步去直流是把信号的整体均值拉回零。第二步去趋势针对的是轴心位置缓慢漂移、温度变化造成的传感器零位移动。polyfit 的阶数不要太高,我一般用 2 阶,阶数太高会把低频故障特征一起拟掉。

如果你读到的信号本身已经是带通滤波后的输出,第三步可以直接跳过。但注意:标题里提到的 widevcb 场景,也就是需要处理很宽频带振动信号的情况,预处理更需要谨慎。宽频带意味着从 0.5 Hz 到 5 kHz 都有能量,如果不先用带通把无关低频和高频噪声剔掉,ITD 会把大量精力浪费在分解环境噪声上,故障分量反而排到很后面。

3.2 多层分解与分量筛选:哪些旋转分量真正值得看

预处理干净之后进入正式分解。前面说过,单层 ITD 只产生一个旋转分量,要做故障诊断需要多层分解。我常用的循环写法如下:

% 多层 ITD 分解:每一轮对残余基线继续拆分 layers = 60; % 最大分解层数 residual = x_detrend; % 初始输入 comps = zeros(numel(x_detrend), layers); % 存储每层旋转分量 for k = 1:layers [baseline, comps(:, k)] = ItdBaseDecomp(residual, 1); residual = baseline; % 残差更新为基线,下一轮继续拆 if std(residual) < 1e-6 * std(x_detrend) break; % 基线能量降到很低时提前停止 end end comps(:, k+1:end) = []; % 截掉未分解的空列

这段代码每一轮把基线 residual 继续传入 ItdBaseDecomp,分解出的旋转分量存到 comps 矩阵里,直到基线的标准差降到原始信号标准差的百万分之一就停止。我把 ItdBaseDecomp 的第二参数传 1,是因为循环外已经控制了层数,函数内部不需要再额外迭代。

分解完之后不能把所有分量都拿去看,否则会被噪声分量带偏。我一般先用相关系数筛一遍:

% 计算每个分量与原始信号的相关性,保留强相关分量 corr_all = zeros(size(comps, 2), 1); for k = 1:size(comps, 2) tmp = corrcoef(x_detrend, comps(:, k)); corr_all(k) = abs(tmp(1, 2)); end keep_idx = find(corr_all > 0.5); keep_comps = comps(:, keep_idx);

阈值 0.5 是我个人的经验值。你可能会发现 60 层里只有前 5 到 8 层的相关系数超过 0.5,剩下的和原始信号相关度很低,基本可以视为分解产生的伪分量或噪声。筛选这一步能显著降低后续误判率,尤其当噪声比较大时,不加筛选直接求特征容易造成虚警。

3.3 特征指标与频谱对照:峭度、RMS、包络谱怎么配合

筛选出有效分量之后,接下来就是从这些分量里提取故障特征。最常用的三个指标:RMS 反映振动能量水平,峭度反映冲击特性,包络谱反映故障特征频率成分。

RMS 计算公式很简单:sqrt(mean(component.^2))。健康状态下 RMS 相对稳定,一旦出现磨损或裂纹,冲击能量会让 RMS 明显上升。峭度是四阶统计量,健康信号接近 3,出现冲击故障后可能跳到 5 甚至 10 以上。包络谱则是把旋转分量做 Hilbert 变换得到包络信号,再对包络做 FFT,就能看出故障特征频率及其倍频。

这三者的配合逻辑是这样:先用 RMS 和峭度做快速筛查——RMS 异常升高说明能量变大,峭度异常说明信号里有冲击成分;再用包络谱定位具体频率,判断故障源在哪个部件。比如轴承外圈故障的包络谱会在 BPFO 及其谐波处出现谱峰,内圈故障则会在 BPFO 基础上叠加转频调制。

实际使用时建议做一个表格记录每个分量的指标值,与健康状态的基线值对比。对比时针式表可以这样设计:每列是一个特征,每行是一个分量。第一次拿到正常数据时先建立基线,之后每次测量都往表里追加一行。一旦某个分量的峭度超过基线值两倍以上,基本可以列入重点关注对象。

4. 工程化封装:把 ITD 从脚本变成可复用的分析函数

4.1 端点外延方式:镜像延拓与线性延拓的取舍

ITD 依赖极值点构造基线,所以极值点在信号端点附近的处理方式直接决定分解结果会不会在两端「飞起来」。这是新手最容易忽略的参数,但恰恰是故障诊断里最影响判断力的细节。

默认的极值查找函数 extrema_1.m 如果只在原始信号内部找极值,那么信号开头和结尾处各只有一个方向的极值点,基线插值到端点附近会缺少参考点,导致旋转分量在两端剧烈摆动。我一般会在调用极值查找前先把信号做镜像延拓:把开头一段信号以第一个点为轴翻转接到信号前面,把结尾一段以最后一个点为轴翻转接到信号后面。延拓长度取信号长度的 2% 到 5%,至少覆盖一个主周期。

如果在封装函数里做,可以在预处理阶段内置这个逻辑。注意:镜像延拓引入的伪信号不要参与最后的特征计算,分解完成后要把结果截回原来的长度。另外,线性延拓虽然实现简单,但会让基线在端点处呈现明显折线效果,对低频分量影响很大,我基本不用。

4.2 迭代停止条件:不要只数层数,要看残余能量

很多人喜欢把层数固定成 60、80 甚至 100,认为分得越细越好。实际层数过高时,后续层分量会被分解成形状规则的伪振荡,频率上几乎没有物理意义,反而干扰筛选。

更可靠的做法是把停止条件建立在残余能量上:

% 停止条件示例:残余能量低于原信号 1% 就停 energy_ratio = sum(residual.^2) / sum(x_detrend.^2); if energy_ratio < 0.01 break; end

这个 1% 的阈值可以根据信号特性调整。信噪比高、故障特征明显的信号可以放宽到 2%;背景噪声大、故障特征微弱的信号建议收紧到 0.5%。注意 I 的书写:层数上限还是要设一个安全值,防止极端情况下循环陷入死循环,我一般上限 60 层。

4.3 封装函数 itd_analyze:输入输出与边界处理

把脚本变成可复用函数,是从「能跑」到「能反复用」的关键一步。我把整套流程封装成一个函数,输入原始信号和参数,输出分量矩阵和筛选结果:

function [comps, keep_idx, baseline] = itd_analyze(x, fs, maxLayer, alpha) % ITD 分析封装:预处理 -> 多层分解 -> 相关性筛选 % 输入: % x 原始振动信号,列向量 % fs 采样率,仅用于趋势项移除时的时间轴构造 % maxLayer 最大分解层数,默认 60 % alpha 基线控制点系数,默认 0.5 % 输出: % comps 筛选后的旋转分量矩阵 % keep_idx 保留的层索引 % baseline 最后残余基线 if nargin < 3, maxLayer = 60; end if nargin < 4, alpha = 0.5; end x = x(:); t = (0:numel(x)-1)' / fs; x = x - mean(x); p = polyfit(t, x, 2); x = x - polyval(p, t); residual = x; comps_full = zeros(numel(x), maxLayer); for k = 1:maxLayer [baseline, comps_full(:, k)] = ItdBaseDecomp(residual, 1); residual = baseline; if sum(residual.^2) / sum(x.^2) < 0.01 break; end end comps_full(:, k+1:end) = []; % 相关性筛选 corr_all = zeros(size(comps_full, 2), 1); for k = 1:size(comps_full, 2) tmp = corrcoef(x, comps_full(:, k)); corr_all(k) = abs(tmp(1, 2)); end keep_idx = find(corr_all > 0.5); comps = comps_full(:, keep_idx); end

这个函数把前面三步整合到了一起。调用方式很简单:[comps, idx] = itd_analyze(x, 25600),拿到的 comps 就是筛选后的有效分量。第 2 章提到的 α 参数在这里也留了接口,想改变基线控制点位置时可以直接传值。封装完成后,我处理任何一批新数据都只用一行代码,不用再翻脚本改参数了。

5. 常见问题与避坑记录:五个让 ITD 翻车的细节

5.1 高频踩坑记录:现象、原因、解决方案

第一坑:分解结果在信号两端剧烈发散,波形像喇叭口。

现象:前 50 个点和最后 50 个点的旋转分量幅值远超中间区域,甚至淹没真实故障特征。

原因:极值查找没有处理端点,基线插值在端点处缺参考点。常见做法是端点不延拓直接插值。

解决:先镜像延拓再分解,分解后截回原长。延拓长度取信号长度的 2%~5%,至少覆盖一个主周期。

第二坑:层数设到 80 层,后 20 层的分量形态规则但完全没有物理意义。

现象:所有层都被筛成高相关分量,包络谱上出现大量幽灵频率,看不出故障在哪。

原因:过度分解后,基线被打成周期性伪振荡,能量占比虽小但相关度不低。

解决:用残余能量占比作为停止条件,阈值设置为原信号能量的 1%;同时把层数上限锁在 60 层以内。

第三坑:100 万点长信号分解到一半 MATLAB 直接内存耗尽。

现象:报错「Out of memory」或者电脑风扇狂转,进程卡死。

原因:comps 矩阵预分配为 1000000 × 60 的 double 数组,占内存约 480 MB,加上中间变量轻松超过 1 GB。

解决:先把信号分段分解,每段 20 万点;或者把 comps 改成单精度;也可以边分解边保存到磁盘。分解长信号时分段是最省事的路子。

第四坑:原始信号带明显转子转速波动,转频随时间缓慢变化。

现象:固定转频对照失效,包络谱峰值被展宽,故障特征频率找不准。

原因:ITD 假设局部振荡是平稳的,转速变化让故障特征频率不断漂移。

解决:先对信号做阶比跟踪或转速归一化,再用 ITD 分解。如果手头没有转速信号,至少把信号切成长度较短的时段分别分析。

第五坑:采样率低,冲击波形被削平,分解后的旋转分量里看不到高频冲击。

现象:峭度偏低,包络谱只有转频,找不到故障特征频率。

原因:采样率低于故障特征频率的 5 倍以上,冲击信号被采样点抹平。

解决:检查采集系统采样率是否足够。轴承诊断一般要保证故障特征频率的 10 倍以上,冲击类故障建议采样率不低于 20 kHz。

5.2 拿到分解结果后必须做的三项验证

分解结果不是直接拿来就算数的。我以前吃过亏:拿一个分量算完峭度直接报警,最后发现是分解伪影。现在每批数据分解完,固定做三项验证。

第一项是重建校验。把所有保留分量的矩阵按列求和,再加上残余基线,得到的信号必须和预处理后的 x 高度一致。误差的标准是逐点绝对误差最大值不超过原始信号幅值的 1e-6。如果对不上,说明某层分解有 bug 或者 comps 矩阵被截断过。

第二项是正交性检查。理想情况下不同旋转分量之间近似正交,即 comps 的协方差矩阵非对角元应该远小于对角元。如果第 2 层和第 3 层高度相关,说明分解层间信息没完全剥开,两层的物理意义会重复,后续特征提取会重复计算同一个故障模式。

第三项是与包络谱对照。取第一有效分量做包络谱,谱峰频率必须能跟已知的转频或故障特征频率对得上。如果第一有效分量的包络谱完全找不到任何已知频率的整数倍,那多半是分解参数没设对,不要继续往深处分析了。

这三项验证做下来不超过五分钟,但能把大部分翻车情况拦在前面。

6. 用包络谱锁定故障特征频率:一个可复用的验证技巧

整个 ITD 分解流程走完之后,最后一步是把故障特征频率从分量里挖出来。我常用的手法是把筛选后的前几个有效旋转分量做能量加权叠加,再对叠加信号做 Hilbert 包络谱分析。这样做比单看一个分量更稳,因为故障冲击往往分布在相邻多个分量里。

% 取前两个有效分量叠加,再做包络谱分析 sig = comps(:, 1) + comps(:, 2); env = abs(hilbert(sig)); % 包络信号 N = numel(env); f = (0:floor(N/2)-1)' * fs / N; % 频率轴 EnvSpec = abs(fft(env)); halfSpec = EnvSpec(1:floor(N/2)); figure; plot(f, halfSpec); xlabel('频率/Hz'); ylabel('包络幅值'); title('ITD 旋转分量叠加包络谱');

这里的叠加系数我一般直接取 1:如果两个分量都通过了相关性筛选,能量占比基本在同一量级,等权叠加不会引入额外偏差。如果你发现早期分量幅值明显大于后期分量,可以按能量的平方根加权,但实测差异不大。

拿到包络谱后,第一件事是对照转频及其倍频——转频处有峰说明存在转子不平衡或轴弯曲;第二件事是对照轴承特征频率 BPFO、BPFI、BSF、FTF,看是否在故障特征频率处出现谱峰。完整诊断时我还会把原始信号的直接频谱、ITD 分量的包络谱、以及健康状态的基线包络谱叠在一张图上,故障频率的增量会非常清楚。

从那以后,我每次拿到新的 ITD 工具包,都强制走一遍完整流程:先跑 example01.m 验证包内函数对得上一篇原始论文的分解逻辑,再用仿真信号算一遍分量筛选和包络谱,最后才碰真实数据。这套验证顺序帮我避开了不少「看起来像故障、其实是参数没校准」的尴尬。希望你拿到这份 ITD 分解 zip 后,也能用同样的思路快速上手,少走我走过的弯路,希望帮到你。

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

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

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

立即咨询