做振动信号分析这些年遇到的最棘手的问题往往不是仪器精度而是怎么从一团乱麻一样的非平稳信号里把特征抠出来。轴承点蚀、齿轮裂纹、转子碰摩这些故障的早期特征都藏在一段段瞬态冲击里单纯看时域波形容易被噪声淹没看频谱又很难定位冲击发生的时间。ITDIntrinsic Time-Scale Decomposition固有时间尺度分解就是我在这种场景下用得越来越顺手的一种时频分析方法。它不需要预先设定基函数直接从信号自身的时间尺度出发把非平稳信号逐层分解成若干个固有旋转分量Proper Rotation ComponentPRC再配合多图像绘制技术把原始信号、各层分量、频谱特征摆在同一张画布里信号的非平稳特征就能被一层层剥开、看得明明白白。这篇文章适合正在做机械故障诊断、振动监测、生物医学信号处理或者刚接触非平稳信号分解算法的读者。我会把ITD的原理、算法流程、Python实现以及“分解完之后怎么用多图把结果讲清楚”这一整套经验都写出来。代码是可以在本地直接跑的案例是仿真出来的所有参数和判断依据我都会说明白希望你看完能少踩几个坑。1. ITD是什么从非平稳信号处理的痛点说起1.1 为什么我放弃纯时域/频域分析先说一个我自己的观察很多做故障诊断的人一开始接触信号分析就是FFT快速傅里叶变换画个频谱找峰值。FFT确实好用但它有一个天然短板——它假设信号是平稳的或者至少在分析窗口内是平稳的。可现实里的振动信号几乎都不平稳转速有波动、负载有变化、故障冲击是瞬态的、背景噪声是非高斯的。你拿FFT去处理一个包含变转速成分的信号谱线会被展宽特征频率被淹没在边频带里有时候连基频都看不清楚。后来我转向了小波变换和短时傅里叶变换STFT这两个属于时频分析方法比纯频域强不少但问题也很明显小波分析的基函数需要先选选错了结果能差出好几个量级STFT的窗长一旦定死就陷入了时间分辨率与频率分辨率不可兼得的困局。说白了这些方法都要求你先假设“信号长什么样”再做匹配分析。但实际工程里我们往往不知道信号到底是什么形态这就很尴尬。1.2 ITD的核心思想把信号“一层层剥开”ITD是Frei和Osorio在2007年提出的一种自适应时频分析方法。它的核心思想非常朴素任何一个复杂信号都可以看成是由不同时间尺度的分量叠加而成的。时间尺度大的是低频趋势时间尺度小的是高频振荡而ITD就像一个“信号分拣员”根据信号本身的极值点分布把一个整体信号逐层剥离成多个反映不同局部时间尺度的分量外加一个单调趋势项。整个过程很像剥洋葱第一层剥出来的是局部波动最剧烈的成分剩下的“洋葱芯”继续被当作新的输入信号再次剥皮如此反复直到最后一层只剩一个单调序列基本没法再拆为止。每一层剥出来的分量就叫固有旋转分量PRC。它跟经验模态分解EMD里的本征模态函数IMF类似但生成机制完全不同。EMD靠的是“包络均值”做筛分需要反复迭代ITD靠的是“基线提取”一次计算就出分量速度快得多。1.3 ITD与EMD的对比快在哪里、稳在何处我自己在工作中经常拿ITD和EMD做对比。最直观的差异就是速度同一段一万点的振动信号EMD可能要做几十次包络拟合、筛分迭代CPU忙活半天ITD因为只用极值点之间的线性插值构造基线一次分解算下来速度通常是前者的好几倍。对需要实时或准实时监测的场景这种计算效率优势非常实际。稳定性上ITD也有一个隐形的优势——它不会像EMD那样容易因为包络拟合过冲而在分量里产生伪振荡。当然这不是说ITD没有缺点它对极值点选取非常敏感端点处容易出现扭曲如果信号噪声太大极值点会成倍增加导致分解出大量无意义的分量。这些我在后面的“常见问题”里会逐一展开。总体下来我的经验是在工程信号分析里ITD配合后续的时域统计特征或谱分析往往比直接上EMD要省心。2. 算法流程深入拆解基线提取与固有旋转分量2.1 基线提取算子三个极值点定一段基线ITD最关键的步骤是“基线提取”。很多人第一次看ITD的公式容易被吓退但说白了它干的事就一句话在任意两个相邻极值点之间的区间上用一条折线或更平滑的曲线把信号的整体趋势“接”出来作为当前尺度上的基线。这个基线怎么定标准做法是先找到信号的局部极大值和局部极小值点然后对于相邻的极值点序列用三个极值点 X_{k-1}、X_k、X_{k1} 的幅值关系计算第 k 个基线控制点的数值。控制点大致可以理解成“这一段基线在这个位置应该落在哪”。控制点算出来后在极值点之间做线性插值就能得到整条基线。原始信号减去这条基线得到的第一个分量就是最高频的那个PRC。我画过很多次这条基线的示意图它对信号形态的跟随效果非常直观当信号向上冲时基线也会跟着略微上抬当信号向下俯冲时基线会稍微下压。这样减去基线后剩下来的PRC就自然地保留住了信号在该尺度上的局部振荡特征。第一次看这个原理时我觉得它其实不像一个高深的数学工具更像是“用手顺着波形画了一条中线”只是这条中线不是拍脑袋画的而是有严格的极值点约束。2.2 PRC的生成与迭代分解PRC的生成流程是这样设原始信号为 x(t)先提取它的极值点再按上面的规则算出基线信号 L(t)然后令第一个分量c_1(t) x(t) - L(t)这个 c_1(t) 就是第一个固有旋转分量。它代表了信号中最局部、最快速变化的成分。接下来把基线信号 L(t) 当成新的“原始信号”重新做极值点提取和基线计算继续分离出第二个分量 c_2(t)。一直重复直到剩下的残余信号 R(t) 是单调序列或者极值点数量太少无法再继续分解为止。整个过程用代码伪码表示就是residual x while not is_monotonic(residual): baseline compute_baseline(residual) prc residual - baseline prcs.append(prc) residual baseline要注意的是每层分量的频率是逐步降低的。第一层PRC往往捕捉冲击、尖端噪声、高频共振越往后的PRC越接近调幅调频慢变趋势最后的残余项基本对应信号的直流偏置或全局趋势。整套分解就像把一个复杂信号里的“快变量”和“慢变量”剥离开来这正是一般FFT做不到的事。2.3 终止条件、端点效应与参数设定ITD的分层不会无休止地进行下去实际使用中需要设定终止条件。我常用的规则有两条一是如果剩余残差的极值点数量少于3个就停止迭代二是人为限定最大分解层数比如默认分解5层或8层就够了超过之后的分量往往幅值极小物理意义也不清晰。端点效应是ITD绕不开的坑。因为极值点都是内部的信号最左端和最右端那一段基线控制点没法用“三个极值点”的正常规则算出来。很多实现里会直接把端点信号值当作控制点或者做镜像延拓处理。实测下来端点效应会导致分解出的PRC在首尾位置出现扭曲甚至小幅扩散但影响范围通常只在信号长度的2%到5%左右。如果信号足够长直接忽略两端如果信号很短最好在分解前做一段镜像延长分解完再截断。参数方面最值得注意的是极值点的检测方式。Python里可以用scipy.signal.argrelextrema也可以用find_peaks。对噪声敏感的信号我先做一轮简单的滑动平均或低通滤波再提取极值否则毛刺会产生大量伪极值点把分解结果搞得支离破碎。3. 多图像绘制实操用Python把ITD结果画到位3.1 环境准备与数据仿真所谓多图像绘制并不是简单地把几幅图拼在一起而是要把ITD分解结果中的有效信息通过合理的排版和标注呈现成一张“看了就能讲清楚问题”的图。我在实际项目里最常用的组合是顶层画原始信号中间层逐个画PRC底层画残差项如果还要分析频域特征就在每个时域子图右侧对齐放一个对应的频谱图或者单独用一行画包络谱。先我把环境准备一下Python 3.8以上装上NumPy、Matplotlib、SciPy三个库就够了。整个流程不依赖特别冷门的包这意味着你可以快速在自己电脑上复现。然后是数据仿真。为了演示非平稳信号的特征提取我构造一段带周期性冲击的滚动轴承外圈故障振动信号import numpy as np import matplotlib.pyplot as plt from scipy.signal import argrelextrema, hilbert fs 20000 # 采样率 20kHz duration 1.0 # 1秒 t np.arange(0, duration, 1/fs) n len(t) # 转频与故障特征频率 fr 25.0 # 轴转频 25Hz对应 1500rpm f_outer 94.0 # 外圈故障特征频率约 3.76×fr # 平稳的转频成分与二倍频 sig 0.8*np.sin(2*np.pi*fr*t) 0.4*np.sin(2*np.pi*2*fr*t) # 构造周期性冲击序列 impact_interval int(fs / f_outer) impulse_train np.zeros(n) impulse_train[::impact_interval] 1 # 系统冲击响应指数衰减正弦振荡 impulse_response np.exp(-1800*t[:500]) * np.sin(2*np.pi*1200*t[:500]) # 卷积得到故障冲击成分 impact_component np.convolve(impulse_train, impulse_response, modesame) # 叠加噪声形成仿真信号 signal sig impact_component 0.15*np.random.randn(n) signal signal / np.max(np.abs(signal)) * 2.0这段代码生成了1秒、20kHz采样率的信号包含25Hz转频及其二倍频、94Hz外圈故障冲击、以及18kHz带宽内的衰减振荡响应和随机噪声。从波形上看它就是一条“看起来有点规律的乱线”但里面其实埋着典型的非平稳故障特征。接下来把它丢给ITD分解。3.2 实现ITD分解核心代码ITD分解的Python实现重点就在“极值点检测”和“基线控制点计算”上。我把核心代码写在这里def extract_extrema(x): 提取局部极大值和极小值对应的索引与幅值 max_idx argrelextrema(x, np.greater)[0] min_idx argrelextrema(x, np.less)[0] idx np.sort(np.concatenate([max_idx, min_idx])) return idx, x[idx] def compute_baseline(x, idx, ext_val): 根据极值点计算基线控制点并插值得到整条基线 k len(idx) ctrl np.zeros(k) for i in range(1, k - 1): x_left ext_val[i - 1] x_mid ext_val[i] x_right ext_val[i 1] denom x_right - x_left if abs(denom) 1e-12: alpha 0.5 else: alpha (x_mid - x_left) / denom ctrl[i] alpha * x_left (1 - alpha) * x_right # 端点控制点直接取极值点本身 ctrl[0] ext_val[0] ctrl[-1] ext_val[-1] # 线性插值生成完整基线 baseline np.interp(np.arange(len(x)), idx, ctrl) return baseline def itd_decompose(x, max_imf8): 固有时间尺度分解返回PRC列表和残差项 prcs [] residual x.copy() for _ in range(max_imf): idx, ext_val extract_extrema(residual) if len(idx) 4: break baseline compute_baseline(residual, idx, ext_val) prc residual - baseline prcs.append(prc) residual baseline return prcs, residual这段代码的关键在于alpha的计算它等于中间极值点在左右极值之间的相对位置决定了基线控制点偏向哪一侧。如果中间点更靠近左极值控制点就更偏向左极值如果更靠近右极值控制点就偏向右侧这样就能保证基线随信号走势自适应地调整。跟我前面说的“顺着波形画中线”完全对应。分解时我加了一个很朴素的简化端点控制点直接取极值点值。这样写代码最短、最容易理解适合学习和二次开发。你在工程中如果觉得端点影响大可以改成镜像延拓先把信号首尾各延长一段极值再分解最后裁剪掉延展部分。3.3 多子图绘制与特征标注分解完成之后最重要的一步就是多图像绘制。我通常会画一张“多层时域堆叠图”第一行放原始信号后面依次放PRC1到PRC几最后一格放残差。用共享x轴的方式保证时间对齐每种成分标注清楚这样一眼就能看到不同时间尺度下的分量形态prcs, residual itd_decompose(signal, max_imf5) num_rows len(prcs) 2 fig, axes plt.subplots(num_rows, 1, figsize(12, 2.0*num_rows), sharexTrue) axes[0].plot(t, signal, lw0.6, colork) axes[0].set_ylabel(原始信号) axes[0].set_xlim([0, 1.0]) for i, prc in enumerate(prcs): axes[i1].plot(t, prc, lw0.8) axes[i1].set_ylabel(fPRC{i1}) axes[-1].plot(t, residual, lw0.8, colortab:red) axes[-1].set_ylabel(残差) axes[-1].set_xlabel(时间 (s)) plt.tight_layout() plt.show()画这种图有个规律越靠前的PRC波形看起来越毛糙频率越高越靠后的PRC波形越平滑。假如信号里有周期冲击你通常能看到高频PRC里出现一串间隔均匀的尖峰间隔的倒数就是冲击频率也就是94Hz。为了进一步确认我通常还会把FFT频谱叠加在PRC旁边或者单独画一个频谱子图。这个习惯在故障诊断里特别重要因为“时域图看得出有没有冲击”和“频谱图能读出周期频率”是两个互补的维度。3.4 绘图细节信息密度、时间轴、频谱联动多图绘制最考验排版和信息管理我详细说说几个要点。第一子图数量不要贪多。分解出8个分量不等于要把8个图全画出来。通常只画包含主要能量或与目标特征相关的4到5个子图其余用文字描述或放进汇总表。画太多只会让图变成一条条细线审阅的人根本看不清。第二时间轴必须对齐且共享x轴。手动设置每个子图自己的x范围容易造成视觉错位错位之后“同一时刻冲击对应哪个分量”就完全没法看了。直接用sharexTrue是省力又不出错的做法。第三频谱联动时要注意纵轴尺度的统一。PRC1和PRC4的幅值可能差几十倍如果不分别做归一化低频分量的频谱会高到把高频分量的细节压成一条平线。我通常在每个分量频谱里做“幅值归一化到该层最大值”并把这个说明写在图注里避免误导。第四标注特征频率。在频谱图上用垂直虚线标出转频、倍频、外圈故障特征频率这些位置再配上文字说明比让读者自己找峰值要友好得多。这也是“多图像绘制”在解析非平稳特征时的价值——你不仅画出了分解结果还把判断逻辑一起画进去了。4. 案例解析滚动轴承故障信号的非平稳特征提取4.1 仿真信号构造与故障特征频率我用3.1节生成的仿真信号做案例。它的构成可以归纳成一个表成分频率/参数物理含义正弦分量125 Hz轴转频旋转体不平衡特征正弦分量250 Hz二倍转频可能对应不对中特征周期冲击94 Hz外圈故障特征频率冲击间隔约10.6ms瞬态振荡1200 Hz 衰减振荡轴承系统固有频率处的调制响应高斯噪声幅值0.15模拟传感器噪声和环境干扰94Hz这个故障特征频率是经验设定值。真实轴承的故障特征频率一般根据滚动体个数、节径、滚子直径和接触角计算公式大致是 BPFO (n_ball / 2) × f_r × (1 - (B_d / P_d) × cos α)。我这里为了演示直接取了近似值算出物理背景合理即可。工程上你必须用轴承手册里的参数精确计算不能拍脑袋。4.2 ITD分解结果分析把这段信号丢进上面的itd_decompose实测分解出5个PRC和一个残差项。各个分量的特征我归纳如下PRC1频率最高波形最杂乱能量集中在约1200Hz附近的瞬态振荡。它忠实捕捉了每个冲击到来时的短时高频响应是判断“何时发生冲击”的关键分量。PRC2幅值相对PRC1小一些主要包含与冲击间隔相关的周期性成分波形上能隐约看到一次冲击之间的包络调制。PRC3频率进一步降低能分辨出25Hz转频及其倍数是旋转部件的基频信息层。PRC4非常平滑主要呈现信号整体慢变趋势。PRC5及残差几乎是一条直线反映直流偏置或极低频漂移对本案例没有诊断价值。这个结果符合我对ITD的预期高频故障冲击被快速剥离到第一层低频旋转信息保留在靠后的分量里。配合多图你可以明显看到“PRC1的冲击间隔是均匀的”数出10个左右尖峰对应1秒时长周期性一目了然。要是信号里同时存在多个故障源比如外圈故障叠加内圈故障不同冲击频率会在不同PRC层被拆开这就是ITD自适应的价值所在。4.3 结合频谱和包络谱的联合诊断单纯看堆叠图只能判断“有冲击”但要确认冲击频率是94Hz还需要频谱验证。对信号做FFT时由于噪声和调制的影响94Hz附近不一定出现干净的谱峰这时候有一个很有用的技巧——对PRC1做包络解调也就是先希尔伯特变换求瞬时幅值再对瞬时幅值做FFT得到包络谱。故障冲击的周期性会在包络谱上形成一个清晰谱峰。我在实际项目里就是这么做的from scipy.signal import hilbert analytic hilbert(prcs[0]) envelope np.abs(analytic) env_spectrum np.abs(np.fft.rfft(envelope)) freqs np.fft.rfftfreq(n, 1/fs)画包络谱时94Hz位置会看到一个明显的峰值旁边可能还有25Hz的边带。这正好说明冲击的重复频率是94Hz同时冲击幅度受到轴转频的调制。这个“PRC1时域上看到均匀尖峰包络谱上看到94Hz主峰”的组合就是一次完整的非平稳特征解析闭环。我建议大家在做这类分析时把“原始信号波形 ITD分层堆叠图 PRC分量包络谱”三张图拼在一起形成一套证据链。只给时域图别人会质疑你凭眼睛判断只给频谱图又会丢失瞬态定位信息。多图联动的意义就在于此。5. 常见问题与避坑实录5.1 分解结果不理想怎么办问题1分解出来的PRC看起来全是锯齿没有明显物理意义。这是我最早用ITD时最常遇到的问题原因基本都是噪声太强。极值点对噪声极敏感一个含噪信号的局部极大值可能有一半是噪声毛刺这会让基线控制点毫无规律分量里自然全是碎屑。解决办法很直接在ITD之前先做轻度平滑处理比如滑动平均、Savitzky-Golay滤波或者只提取幅值超过某个阈值的极值点。注意平滑窗口不要太大否则会抹掉真实的冲击尖峰。问题2端点产生剧烈震荡边缘分量幅值异常大。这就是端点效应在作祟。我的建议是分解前对信号首尾做一定比例的镜像延拓比如各延长一个波长分解完再截掉。如果懒得延拓也可以在画图时把时间轴范围向内缩一点忽略两端失真区域。判断失真区间的经验标准是观察PRC最两端的幅值是否明显大于中间段且呈喇叭状如果是就说明端点效应影响到了这里。问题3某个PRC和原始信号的相关系数极低看不出它到底代表什么。这种情况往往是分解层数过多导致的“过分解”。后面几层PRC幅值极小能量占比连1%都不到纯属把数值噪声也拆开了。实用的做法是每分解一层就算一下该PRC与原始信号的相关系数或者能量占比低于某个阈值比如5%就停止分解。这样能有效避免最后画出一堆无意义的细线。5.2 可视化阶段最容易踩的坑多图绘制里我踩过最大的一个坑是不同PRC的幅值尺度差异。ITD前面几层分量的幅值通常远大于后面几层如果所有子图用同一个y轴范围后面的分量会“一条直线”贴在上面完全看不出波形。这个问题的解法在3.4节已经提过每层可以单独归一化但一定要在图里明确标出每层的幅值范围否则看图的人会误以为所有分量的真实幅值大小接近。另一个坑是子图之间的视觉干扰。当子图数量超过6个每层线宽又差不多时整张图看上去就像一捆意大利面。我的个人习惯是高频分量用细线低频分量和残差用相对粗的线关键特征所在的分量用深色突出辅助分量用浅灰色降低存在感。还有一个许多人会忽略的细节——时间轴刻度上的冲击对齐。比如你要说明“PRC1里的每个尖峰对应原始信号里的冲击”最好在原始信号图上用浅色竖线标注冲击发生的时刻并让这些竖线延伸到下面的PRC子图。这样比口述“你看这个尖峰”有说服力得多。5.3 进阶思路ITD不是终点我接触过很多项目用ITD分解完就把结果撂在一边这是浪费了它的价值。真正有意义的做法是把分解出来的分量当作“特征输入”继续往下走。比如对每个PRC提取时域统计特征均方根值、峰值因子、峭度、脉冲指标。轴承早期故障时高频PRC的峭度会明显上升这个信号比原始波形峭度更灵敏。对PRC1这类高频分量做包络谱分析提取故障特征频率及其边带分布。将多个PRC的能量占比、频谱质心、瞬时带宽整理成特征向量喂给SVM、随机森林或轻量神经网络做自动故障识别。在监测系统里把ITD分解结合归一化脉冲能量做成“健康指数”趋势曲线比单看原始信号均方根值更能捕捉早期微弱变化。我个人的体会是ITD的价值不在于它“比EMD高级”而在于它在计算效率和时频分辨率之间找到了一个很实用的平衡点特别适合工程场景里“数据量大、故障微弱、趋势不明”的窘境。配合多图像绘制你既能给领导讲“这个信号确实有问题”又能给算法工程师提供可量化的特征输入算是一套从分析到表达都比较完整的方法链。最后分享一个小习惯不管用哪种时频分解方法拿到结果先不要急着堆代码而是在纸上大概画出信号要分成哪几层。一两个小时后当你发现ITD拆出来的结果和你预判基本一致时你对这个算法的信任感自然就建立起来了。这个信任感是你后续在复杂故障信号里敢不敢依赖ITD的关键。