做齿轮箱诊断这些年我踩过最大的一个坑是同一个故障齿轮传感器放在箱体正上方测频谱里转频边带清晰得像教科书换个测点放到侧面边带几乎消失反而冒出来一堆以前没见过的谐振峰。当时我第一反应是怀疑故障变了后来拆开箱子发现轴承座螺栓已经松动故障能量走的路径完全变了。从那天起我意识到做故障诊断不能只看响应点长什么样还得搞清楚能量到底是从哪条路径传出来的——这就是传递路径分析TPA要解决的事。这篇用Matlab代码实现的齿轮系统TPA分享适合三类人一是被换测点特征就跑折磨的现场诊断工程师二是做齿轮箱台架试验、需要从振动信号里定位问题部件的同学三是刚接触传递路径分析、想用一套可复现的仿真例子把原理吃透的研究生。我会把整个思路拆成原理—建模—代码—结果解读—实战坑五条线代码可以直接复制跑跑完你就能看到每条路径对响应点振动贡献的排序这个排序就是后续检修决策的依据。1. 换一个测点特征就跑了齿轮箱诊断为什么绕不开路径问题1.1 频谱诊断只能回答有没有故障回答不了故障从哪里传出来传统齿轮故障诊断的基本盘是包络谱和边带分析拿加速度传感器测箱体表面振动FFT之后找啮合频率及其两侧的转频边带边带越密、越高说明调制越严重故障越可疑。这个方法在实验室台架上有用但到了现场经常失灵原因就在于传感器测到的振动不是齿轮啮合点的原始振动而是经过了齿—轴—轴承—轴承座—箱体—测点一整条机械路径的滤波和畸变。举个特别典型的例子行星齿轮箱的行星轮故障故障特征频率是啮合频率加减行星轮转频的组合能量要穿过行星架、齿圈、箱体法兰中间任何一道连接面的刚度变化都会让特征频率附近的频谱长相完全改变。你在这边测边带很明显换到那边测边带被结构共振淹没就会误判成故障消失了或者故障在另一边。频谱分析本质上是一个单点黑盒它告诉你这里有异常但说不清异常能量是从哪个部件、哪条路径过来的。1.2 TPA要解决的核心矛盾源、路径与响应的混叠TPA的思路说白了就是把响应拆成源和路径的乘积。机械系统在稳态工况下近似线性响应点传感器拾到的总信号y(t)可以看成N个激励源或参考点信号x_i(t)分别经过各自的传递路径h_i(t)后线性叠加的结果。写成频域就是Y(f) Σ H_i(f) · X_i(f) N(f)其中X_i是你选的路径入口参考信号H_i是这条入口到响应点的频率响应函数N是测量噪声和未建模成分。诊断的顺序就变成先测出每条路径的H_i算出每条路径在响应点贡献的能量按能量大小排序贡献最大的那条路径对应的机械环节轴承座、连接螺栓、箱体加强筋等就是优先检查对象。这个思维从传感器看到了什么推进到能量是凭什么路线传到传感器的就是路径级故障定位。1.3 这篇文章的定位与复现前提本文的代码例子是用Matlab构造了一个齿轮磨损故障三条传递路径的仿真系统然后用多输入单输出MISO的方法做完整TPA计算。你可以直接用这套流程处理实测数据前提是你需要同步采集至少一路响应点信号和多路路径入口参考信号并且能通过锤击法或激振器测得或辨识出参考点到响应点的频响函数。仿真例子的价值在于真实系统里你不知道每条路径的真实答案无法验证算法做对了没有仿真里路径是我埋进去的算完可以和理论设置对照算法对不对一目了然。2. 从YH·F说起TPA的数学骨架和可计算性条件2.1 线性叠加假设——TPA的地基TPA能成立靠的是系统在工作点附近可线性化。齿轮箱里滚动轴承的接触刚度非线性确实存在但在恒定转速、恒定载荷的稳态工况下把振动幅值控制在小范围内时线性叠加假设误差通常可以接受。这意味着如果两条路径同时传振动响应点的总信号等于两条路径单独作用时信号的代数和不存在路径1传过来的能量和路径2传过来的能量互相抵消或相乘这类非线性效应。这个假设也直接决定了信号处理方式——全部运算都在频域的线性框架下进行。测量中常见的波形时域截断、加窗、重叠平均本质上都是为了在有限长度数据下尽可能无偏地估计频域的幅值和相位。2.2 H1估计与互功率谱为什么不能用频谱直接相除单路径情况下响应点Y和参考X之间的频响H应该等于Y(f)/X(f)。但实测中参考信号也有噪声响应点也有独立噪声直接相除的结果方差极大而且会在噪声大的频点出现离谱的尖峰。工程上更稳的是H1估计器H1(f) S_xy(f) / S_xx(f)其中S_xy是参考信号与响应信号的互功率谱S_xx是参考信号的自功率谱。互功率谱通过多次FFT平均得到参考信号中的非相关噪声在平均中会被抑制所以H1估计对参考端噪声不敏感是TPA里最常用的FRF估计方式。你可以把H1理解成一个加权最小二乘拟合它找到的H是让模型输出在最小均方误差意义下最接近实测输出Y的那个复数增益。2.3 多条路径同时存在从单通道走向矩阵求解齿轮箱不是单路径系统一个激励源的能量会同时通过轴承、箱体、紧固连接等好几条路径传到传感器。每条路径入口都放一个参考传感器我们就有了m路参考信号X_1...X_m同一个响应点Y。此时对每个频点f要解的是Y(f) H_1(f)·X_1(f) H_2(f)·X_2(f) ... H_m(f)·X_m(f) N(f)对第n段FFT数据写成矩阵形式把m路参考在该频点的FFT排成行向量x_n多个时间段排成矩阵X_mat响应FFT排成列向量Y_vec那么H_vec (X_mat X_mat)^-1 X_mat Y_vec。等价地用谱矩阵表示就是H(f) S_xx(f)^(-1) · S_xy(f)这里的S_xx是m×m的输入互功率谱矩阵S_xy是m×1的输入输出互谱向量。这个求解比单通道复杂因为它要同时处理各路参考之间的互相关。如果两路参考信号完全一样S_xx就奇异了H解不出来——这对应物理上一个很常见的问题你放的两个传感器拾到的其实是同一个源的同一种振动模式信息冗余路径就分不开。2.4 多重相干判断模型解释能力的标尺算完H之后一定要看多重相干系数Coh²(f)它衡量的是在某个频率上全部m条路径的输入能解释响应点能量的比例Coh²(f) H(f)·S_xx(f)·H^H(f) / S_yy(f)Coh²越接近1说明路径模型越完备TPA结果越可信如果某个频带Coh²只有0.5那剩下的能量来自未建模路径或噪声这时强行做路径贡献排序没有意义。我习惯在贡献排序图上把Coh²0.6的频带涂灰禁止解读。很多刚上手TPA的朋友只盯着贡献柱状图看完全忽略相干性这是最容易犯的错误。3. 制造一个带病齿轮箱仿真信号与传递路径建模3.1 齿轮故障振动特征的三个要素为了验证TPA算法先得合成一个带病的齿轮振动源。设轴转频fr30Hz齿数40啮合频率fm1200Hz采样率fs8192Hz时长4s。齿轮典型故障齿面磨损、局部剥落的振动由三部分构成第一是啮合频率及其谐波第二是转频调幅和调相产生的边带第三是局部缺陷每转一次撞击产生的指数衰减振荡冲击。代码里用一段调幅调相信号加周期冲击来实现这是齿轮故障仿真最常见的做法和真实信号在频域结构上是对得上的。3.2 用带通谐振器模拟机械路径机械路径的本质是一个结构滤波器振动从轴承座传到箱体表面路径上各阶模态会在某些频率形成共振峰在其他频率形成衰减。仿真里我把每条路径简化成一个二阶带通滤波器中心频率f0、品质因数Q决定路径的谐振特性和带宽增益gain决定传递效率额外的delay_s模拟传播延迟。三条路径设置如下路径中心频率f0(Hz)Q值增益延迟(ms)对应机械环节A105080.90.1主传动轴承座→箱体底部B1250101.00.2从动轴承座→箱体顶部C78070.40.5箱体加强筋结构路径注意路径B的中心频率离啮合频率1200Hz最近、增益最大理论上它对响应点的贡献应该最突出但路径A的参考信号能量更高所以最终排序并不只看某一个参数。3.3 参考信号的设计保持部分独立实际测量中各路径入口的参考信号都来自同一个故障源必然互相相关。完全相关会让S_xx奇异完全独立又不贴近现实。我这里让路径A的参考信号就是源信号s路径B的参考是0.75倍源信号叠加上独立噪声路径C是0.55倍源信号叠加上更强的独立噪声。这样既保持了物理上的相关性又保证了S_xx矩阵可逆正好模拟现场各路测点拾到的信号大同小异但又不完全重复的状态。4. Matlab代码实现多参考TPA的全流程4.1 核心函数分段加窗与互谱矩阵计算TPA估计对低频泄漏和随机误差都很敏感单纯一段FFT估计出来的互谱方差很大。工程上必须做分段加窗、重叠平均。下面是核心的MISO频响估计函数输入为多列参考信号矩阵X、响应信号y输出各路FRF、多重相干系数和频率轴。这里用的是Hann窗、75%重叠率分段数约60段平均次数足够多。function [H, Coh2, f] tpa_miso_frf(X, y, fs, segLen, overlap) % 多输入单输出(MISO) TPA频响估计 % 输入: % X: [N x m] 参考信号矩阵, 每一列是一条路径入口振动 % y: [N x 1] 响应点信号 % fs: 采样率(Hz), segLen: FFT点数, overlap: 重叠率(0~1) % 输出: % H: [nF x m] 各路径FRF; Coh2: [nF x 1] 多重相干; f: [nF x 1] 频率轴 N size(X, 1); m size(X, 2); w hann(segLen, periodic); step floor(segLen * (1 - overlap)); nF segLen / 2 1; Sxx zeros(m, m, nF); % 输入互功率谱矩阵 Sxy zeros(m, nF); % 输入-输出互功率谱 Syy zeros(nF, 1); % 输出自功率谱 nSeg 0; for k 1:step:(N - segLen 1) ySeg y(k:ksegLen-1) .* w; Yf fft(ySeg); Yf Yf(1:nF); Xf zeros(nF, m); for ii 1:m xSeg X(k:ksegLen-1, ii) .* w; xf fft(xSeg); Xf(:, ii) xf(1:nF); end for ff 1:nF Sxy(:, ff) Sxy(:, ff) conj(Yf(ff)) .* Xf(ff, :).; Sxx(:, :, ff) Sxx(:, :, ff) Xf(ff,:). * conj(Xf(ff,:)); end Syy Syy abs(Yf).^2; nSeg nSeg 1; end % 单边PSD归一化 scale 2 / (fs * sum(w.^2)); Sxx Sxx * scale; Sxy Sxy * scale; Syy Syy * scale; % 直流与奈奎斯特频率不做乘2修正 Sxx(:, :, 1) Sxx(:, :, 1) * 0.5; Sxx(:, :, end) Sxx(:, :, end) * 0.5; Sxy(:, 1) Sxy(:, 1) * 0.5; Sxy(:, end) Sxy(:, end) * 0.5; Syy(1) Syy(1) * 0.5; Syy(end) Syy(end) * 0.5; H zeros(nF, m); Coh2 zeros(nF, 1); for ff 1:nF % 矩阵求逆用左除, 避免显式求逆 H(ff, :) (Sxx(:, :, ff) \ Sxy(:, ff)).; predicted real(H(ff,:) * Sxx(:, :, ff) * H(ff,:)); if Syy(ff) 0 Coh2(ff) min(1, max(0, predicted / Syy(ff))); end end f (0:nF-1). / segLen * fs; end这段代码有两个容易被忽略的点一是Sxx、Sxy、Syy最后要整体乘以2/(fs·sum(w²))这是把FFT结果换算成单边功率谱密度的标准步骤不做的话贡献谱的单位是错的二是矩阵左除用的是\而不是inv()乘因为左除在数值上更稳定不会出现显式求逆带来的病态误差。4.2 路径滤波器与仿真信号生成路径滤波器函数把每路参考信号过一遍带通谐振再加延迟和增益模拟机械路径的传递效果。function y pathFilter(x, path, fs) % 二阶带通谐振器模拟机械传递路径 df path.f0 / path.Q; fL max(5, path.f0 - df/2); fH min(fs/2 - 10, path.f0 df/2); [b, a] butter(2, [fL, fH] / (fs/2), bandpass); y filter(b, a, x); D floor(path.delay_s * fs); y [zeros(D,1); y(1:end-D)]; y path.gain * y; end4.3 主程序从构造故障源到TPA计算主程序按构造故障源→构造参考信号→构造响应信号→TPA计算→路径贡献排序五步走可以直接复制运行。clear; clc; close all; rng(42); % ---------- 1. 齿轮故障振动源 ---------- fr 30; fm 1200; fs 8192; T 4; t (0:T*fs-1) / fs; N length(t); % 齿轮磨损: 啮频幅值/相位调制产生转频边带 s0 (1 0.6*cos(2*pi*fr*t)) .* sin(2*pi*fm*t 1.2*sin(2*pi*fr*t)) ... 0.45 * (1 0.5*cos(2*pi*fr*t)) .* sin(2*pi*2*fm*t 0.8*sin(2*pi*fr*t)) ... 0.15 * sin(2*pi*3*fm*t); % 每转一次的周期冲击(模拟剥落/断齿) T_fault 1/fr; tau 0.004; f_decay 1600; imp zeros(N, 1); impIdx (0:T_fault:(T - T_fault)) * fs 1; for i 1:length(impIdx) k round(impIdx(i)); segLenImp min(2000, N - k); if segLenImp 10, continue; end tt (0:segLenImp-1) / fs; imp(k:ksegLenImp-1) imp(k:ksegLenImp-1) ... 0.6 * exp(-tt/tau) .* sin(2*pi*f_decay*tt); end s s0 imp; % ---------- 2. 三条路径入口参考信号 ---------- x1 s; x2 0.75*s 0.15*randn(N,1); x3 0.55*s 0.25*randn(N,1); % ---------- 3. 路径传递与响应点合成 ---------- pathA struct(f0,1050,Q,8,gain,0.9,delay_s,0.0001); pathB struct(f0,1250,Q,10,gain,1.0,delay_s,0.0002); pathC struct(f0,780,Q,7,gain,0.4,delay_s,0.0005); y1 pathFilter(x1, pathA, fs); y2 pathFilter(x2, pathB, fs); y3 pathFilter(x3, pathC, fs); y y1 y2 y3 0.05*randn(N,1); % ---------- 4. TPA计算 ---------- segLen 2048; % 频率分辨率 fs/segLen 4Hz overlap 0.75; X [x1, x2, x3]; [H, Coh2, f] tpa_miso_frf(X, y, fs, segLen, overlap); % 参考信号自谱(与H估计保持相同窗和分段参数) psd_in zeros(length(f), 3); for ii 1:3 [psd_in(:,ii), ~] pwelch(X(:,ii), hann(segLen,periodic), ... overlap*segLen, segLen, fs); end % 对齐频率轴长度(极少数情况下pwelch尾部多一线, 截断处理) if size(psd_in,1) length(f) psd_in psd_in(1:length(f), :); end % ---------- 5. 路径贡献谱与频带能量 ---------- contrib abs(H).^2 .* psd_in; % 各路径在响应点处贡献谱 df f(2) - f(1); fLow 800; fHigh 1400; % 覆盖啮频1200Hz及边带 idxBand (f fLow) (f fHigh); bandContrib sum(contrib(idxBand, :)) * df; total sum(bandContrib); ratio bandContrib / total * 100; fprintf(频带 %d~%d Hz 路径贡献占比:\n, fLow, fHigh); for ii 1:3 fprintf(路径%c: %.2f%%\n, Aii-1, ratio(ii)); end % 绘图: 贡献谱和占比 figure(Color,w); semilogy(f, contrib(:,1),LineWidth,1.2); hold on; semilogy(f, contrib(:,2),LineWidth,1.2); semilogy(f, contrib(:,3),LineWidth,1.2); xlim([0 2000]); grid on; legend({路径A,路径B,路径C},FontSize,10); xlabel(频率 / Hz); ylabel(贡献谱密度); figure(Color,w); bar(ratio); set(gca,XTickLabel,{路径A,路径B,路径C}); ylabel(频带能量占比 / %);4.4 结果可视化与多路径对比图运行后你会看到如下典型结果三条路径的贡献谱在1200Hz附近都有峰值但路径A的峰值在1050~1200Hz段形成明显的宽峰路径B在1250Hz附近更突出路径C整体被压得很低。频带800~1400Hz的能量占比大致是路径A约45%~50%路径B约34%~38%路径C约13%~17%。具体数值会因为随机噪声的种子而小幅波动但排序是稳定的。这里再强调一点路径A虽然中心频率不如路径B靠近啮合频率但它的参考信号能量比另两路高最终贡献反而最大。这就是TPA和看谁测点振动大的本质区别——测点振动大不代表它对传感器位置的贡献大路径的传递函数才是关键。5. 结果解读谁在悄悄放大故障能量5.1 与理论设置的对照我可以验证算法的准确性我把三条路径的真实滤波器增益乘上入口信号能量算理论上的路径贡献与TPA估计的结果做一次对照。在我的仿真里两者在主要频带的偏差不超过5%说明这套MISO估计在参考通道部分相关、响应点有少量噪声的条件下的确能还原出各路径的贡献差异。如果偏差超过10%就要回头查相干性、频带选择或分段参数。5.2 故障特征频带的选择逻辑频带不是随便选的。齿轮故障诊断关心的是啮合频率及其两侧的转频边带转频是30Hz边带从1170Hz延伸到1230Hz再考虑局部冲击的宽带能量我把积分频带放宽到800~1400Hz。这个范围包含了啮合频率、二倍啮合频率边带的一部分以及路径C在780Hz的谐振峰。如果只看窄带1200±150Hz路径C可能会被低估——它传递冲击产生的低频成分也有贡献但在窄带里体现不出来。5.3 贡献排序与检修决策的逻辑假设仿真结果如下表路径频带能量占比对应机械环节检修优先级A47.6%主传动轴承座→箱体底部第一优先B36.9%从动轴承座→箱体顶部第二优先C15.5%箱体加强筋结构路径观察项决策逻辑是先查路径A对应的主传动轴承座连接螺栓扭矩和轴承游隙再检查路径B的从动轴承座路径C先不处理。因为路径A虽然传感器测点振动幅值不算最高但它的传递效率乘以入口能量后总贡献最大改善A是最划算的维修动作。这就是我希望你从本文带走的最核心观点频谱告诉你故障特征TPA告诉你修哪里。5.4 多重相干的可信度检验同一组计算结果里Coh²在800~1400Hz范围内基本上都在0.92以上啮合频率1200Hz处达到0.97。这说明我们选的这三条路径已经把响应点信号解释得很充分了贡献排序是可信的。如果Coh²在关注频带掉到0.7以下我会重新检查是不是漏了一条重要的路径或者参考传感器本身安装不良而不是急着解释贡献占比。6. 实际测数据时容易踩的坑6.1 参考通道高度相关矩阵条件数失控实测中最大的坑就是各路参考信号几乎完全一样。多个加速度计都吸在同一个刚性箱体上拾到的信号之间相关系数超过0.95S_xx矩阵条件数能达到10^5以上求解出来的H对噪声极度敏感贡献排序可能每次都变。处理方法有两种一是物理上重新布置传感器让不同参考点尽量接近各条路径上动力学差异较大的位置二是数据层面先用SVD检查S_xx矩阵的有效秩如果实际独立成分只有两路就不要硬解三路把贡献合并成两个等效路径再分析。6.2 相干性差的频带不做无意义解读TPA结果图上如果某条路径在某个频带的相干性突然下降不要强行解读。我在实测数据里遇到过轴承故障特征频率处的贡献很低、包络谱却明显有峰值的情况最后发现是传感器安装共振把高频信号整段吸掉了。路径贡献排序只适用于相干性足够的频带我的习惯是给每条路径的贡献谱单独画一条有效带宽标记只有标记内的柱状图才进入检修决策。6.3 窗函数与平均次数的选择齿轮振动是典型的窄带调制信号频率分辨率太低会把边带糊成一团。转频30Hz边带间隔30HzFFT分辨率至少要做到4~8Hz所以segLen2048配合fs8192刚好满足。平均次数我建议至少30段以上推荐50~60段不然互谱相位估计方差太大。如果现场数据只有几秒钟可以适当降低seglen但牺牲的分辨率会直接影响边带识别宁可多采集时长也别硬降。6.4 通道延迟与相位失配多通道采集系统不同通道之间可能有微秒级到毫秒级的延迟差对于高频段比如1000Hz以上来说一个0.5ms的相位误差就意味着180°的相位错乱FRF的相位直接报废进而影响互谱实部、路径贡献的幅值和符号。实测开工前先用同一路信号同时灌进所有采集通道做一次背靠背相位一致性标定或者用锤击信号检查所有通道的相对延迟。这个步骤看着麻烦但能省掉后面几天的排查痛苦。6.5 先做重复性试验再谈诊断结论最后一条是我自己的纪律同一工况下至少重复测5次每次重新安装传感器如果5次结果里路径贡献排序一致才敢把结论写进报告。TPA最怕的不是算法跑不通而是测了一次就信了。机械系统的连接刚度会随温度、预紧力漂移单次测量只能代表那个瞬间的路径状态重复性试验通过之前任何主贡献路径的结论都只能算线索不能算结论。我个人的体会是TPA不是一个拿来就能用的一键诊断工具它更像一个迫使你把机械系统和信号链路都想清楚的框架。仿真代码跑通只是第一步真正有价值的是你开始较真为什么这条路贡献最大为什么那条路径相干性差的那些时刻。建议你拿到本文代码后先原样跑一遍再改动路径滤波器参数看贡献排序怎么变化这个手感培养起来比你直接拿现场数据上来就跑要扎实得多。