搞设备状态监测和故障诊断这些年齿轮箱这个“硬骨头”我啃的次数最多。齿轮箱一响黄金万两的代价不是谁都付得起但真正麻烦的不是坏了而是坏了之后你找不准它在哪里坏、振动是怎么传出来的。大多数现场工程师会遇到一个很实际的困惑传感器就装在轴承座上传感器测到的振动信号里混着齿轮啮合激励、轴承故障冲击、轴系不平衡、箱体共振等一大堆东西你想定位的故障源可能离测点还有两三个传递环节。这时候如果只盯着传感器信号做频谱分析往往看到的是多条传递路径叠加后的“混合体”很难直接判断故障到底出在内圈还是外圈、是齿轮齿面磨损还是轴弯。这个项目要解决的就是用传递路径分析Transfer Path Analysis, TPA配合Matlab仿真把齿轮系统的振动传递路径拆开、算清贡献量再反推故障源的位置和性质让故障诊断不再是“猜”。这篇内容适合正在做齿轮箱故障诊断研究的学生、设备状态监测工程师以及想把TPA方法用Matlab落地跑通的入门选手。我会从TPA的核心思路、齿轮系统的路径建模、Matlab代码怎么一步步实现再到实际调试中容易踩的坑完整走一遍。1. 方案拆解齿轮故障诊断为什么要引入传递路径分析1.1 传感器测到的信号比你想的复杂得多先说一个被很多人忽略的事实加速度传感器测到的振动响应并不是故障源本身的振动而是故障激励经过若干结构路径衰减、调制、延迟之后到达测点的结果。齿轮箱里齿轮啮合产生周期性的啮合力这个力通过齿轮→轴→轴承→箱体→测点安装座最后才被传感器记录下来。不同路径的刚度、阻尼、质量分布完全不同同一激振力在不同路径上的传播效果可以差出几个数量级。这就带来一个诊断上的麻烦如果一个测点主要受某条路径支配而故障源又不在那条路径上你得到的频谱可能压根反映不出故障特征。我见过不止一个现场案例频谱图上滚动轴承外圈故障特征频率倍频成分非常清晰但拆开后发现故障在齿轮上——因为齿轮啮合激励通过一条低衰减路径把轴承故障的微弱冲击“盖”掉了。TPA的价值就在这里它从“响应”反推“源”把每条路径单独拎出来算贡献告诉你传感器测到的能量到底是从哪个激励点、哪条路径来的。1.2 TPA在齿轮系统上的工作逻辑传递路径分析的本质是建立“激励源—传递路径—响应点”三者的线性关系。在齿轮箱这个场景里激励源可以理解为齿轮啮合产生的动态力、轴承故障产生的冲击力、不平衡产生的离心力等。这些力分别作用在不同位置经过各自的传递路径结构传播在某个响应点叠加形成测点振动。如果把问题放到频域响应点和激励点之间的关系可以写成Y(f) H1(f)*F1(f) H2(f)*F2(f) ... Hn(f)*Fn(f)其中Y是响应点频谱Hi是第i条路径的频率响应函数FRFFi是第i个激励力。TPA要干的事情有两件第一识别出每条路径上的等效激励力Fi第二计算每条路径对响应点贡献量的占比找出最“活跃”的那条路径。在齿轮故障诊断里一旦确定了主导路径对应的激励源位置故障定位就完成了大半。1.3 为什么用经典TPA而不是盲目上深度学习现在很多人一聊故障诊断就想到深度学习、端到端分类但工程上有一个很现实的点深度学习需要大量带标签的故障样本还需要在同一个设备上重新训练才能迁移。而齿轮箱的类型、转速、载荷变化太复杂很多时候你根本没有足够的历史故障数据可训练。TPA这一类基于物理模型的方法样本需求低可解释性强算出的路径贡献量本身就是物理量可以直接和轴承特征频率、齿轮啮合频率对标。当然TPA也不是没有代价——它需要频响函数测量需要一个相对可靠的激励源-路径模型对测点布置、数据质量的要求都不低。我在项目里通常的做法是“TPA做粗定位频谱细化特征频率做精确认证”两者配合比单用任何一种都稳得多。2. 核心原理与实操要点TPA的数学底子和工程细节2.1 频响函数矩阵和矩阵求逆法做TPA绕不开频响函数FRF。FRF描述的是“某个激励点用力激励在某个响应点会产生多大的振动响应”是系统的固有属性和激励源无关。测量FRF最常用的手段是锤击法或者激振器法在激励点用力锤敲一下同时记录激励力和响应点的加速度信号两者的频谱比值就是FRF。齿轮系统TPA的难点在于激励源往往不止一个。齿轮啮合点是源两个轴承座也是源甚至还可能把电机的电磁激励也算进来。每个源都对应一个激励点每个激励点对应一个FRF。如果系统里有n个激励点、m个响应测点就构成一个m×n的FRF矩阵。在工程上为了能反演出n个激励力通常要求响应测点数量m≥n然后通过矩阵求逆或者最小二乘来识别载荷F(f) pinv(H(f)) * Y(f)这里pinv是伪逆。算出等效载荷后第i条路径在测点处的贡献量就是Hi(f)*Fi(f)。用矩阵求逆法做载荷识别有一个隐藏前提各个激励力之间必须是互不相关的或者至少频域上是线性独立的。否则矩阵会出现病态求逆结果对噪声极其敏感算出来的Fi会“爆炸”。2.2 测点布置和FRF测量的工程经验测点布置这件事直接决定TPA的成败。我在做齿轮箱TPA时响应点一般选在轴承座上方或箱体刚性较强的位置尽量贴近力传播的主路径不要放在箱体大平面的中央——那个位置模态很多测出来的FRF锯齿感极强而且对传感器质量、粘贴方式特别敏感。激励点的设置则和激励源建模有关。对于齿轮啮合振动可以在齿轮啮合点附近比如从箱体外壁对应啮合区域的位置做激励对于轴承故障冲击通常把激励点放在轴承座上。这里需要说明实际结构里齿轮啮合点是不方便直接敲击的工程上常常用“等效激励点”替代在力传递主路径上的某个可接触位置测量FRF然后把这个FRF做距离修正或局部修正近似为驱动点FRF。这种做法本身有误差但做故障贡献排序是够用的。2.3 载荷识别与贡献量计算的常见坑载荷识别看着简单实际调试时很容易踩坑。首先是矩阵逆的数值稳定性。FRF矩阵的条件数如果偏大比如超过100求逆结果基本不可信。解决办法有增加测点数量用最小二乘稳定解、对FRF矩阵做正则化处理比如Tikhonov正则化、或者改用工况传递路径分析OPA作为替代。其次是频率分辨率。齿轮啮合频率和故障边频带的间隔往往就是轴频比如转频25 Hz那么要想分辨啮合频率旁边的25 Hz边带频率分辨率至少要达到1 Hz以下对应FFT窗长至少1秒。如果采样时长不够边频和啮合主峰糊在一起贡献量排序会失真。2.4 源相关性和部分相干问题齿轮系统里有一个经常被忽略的问题多个激励源之间并不是完全独立的。比如齿轮啮合激励和轴系不平衡激励都跟转速相关它们的频率在轴频的整数倍上会重叠。此时用矩阵求逆法算出的“等效载荷”并不能严格代表物理激励力而更像是一个“等效集总源”。这不是说TPA就不成立了——对于故障诊断我们关心的是“哪个路径对响应贡献大”而不是“谁是物理上真正的力源”。等效集总源仍然可以指示故障位置所在的路径。但如果两条路径的FRF在关心的频带内高度相似矩阵求逆会出现明显的数值振荡。这种时候我会给FRF矩阵加入少量频带正则化或者在多个工况下分别求解再取平均能有效压低伪峰值。3. Matlab代码走一遍从仿真信号到贡献量排序3.1 参数设置与故障激励源信号模拟仿真阶段的任务是构造一个有“已知答案”的齿轮系统用来验证TPA算法对不对。我习惯先写一个公共参数模块把所有关键特征频率提前算好后面所有环节都引用这些变量避免在代码里到处塞“魔法数”。%% 参数定义 fs 20000; % 采样率 20 kHz T 4; % 信号时长 4 s t (0:1/fs:T-1/fs).; n length(t); % 齿轮副参数 Z1 42; % 输入轴主动齿数 Z2 65; % 输出轴从动齿数 n1 1500 / 60; % 输入轴转速单位 r/s 25 Hz fm Z1 * n1; % 啮合频率 1050 Hz fr n1; % 输入轴转频 25 Hz % 模拟齿轮局部故障齿面剥落/断齿产生的调制边带 % 典型特征: 啮合频率 fm 两侧出现 fr 间隔的边频带 sideband_amp [0.10, 0.22, 0.35, 0.50, 1.00, 0.45, 0.28, 0.15, 0.08]; sideband_freq fm (-4:4) * fr; F_excite zeros(n, 1); for k 1:numel(sideband_freq) F_excite F_excite sideband_amp(k) * sin(2*pi*sideband_freq(k)*t 2*pi*rand); end % 加入少量随机冲击成分模拟轴承早期故障的周期性冲击 T_bearing 1 / (Z2 * n1); % 输出轴轴承故障间隔示意 impulse_train zeros(n, 1); idx round(1 : T_bearing*fs : n); impulse_train(idx) 1; F_excite F_excite 0.15 * conv(impulse_train, exp(-t(1:round(0.002*fs))*1000), same);这段代码本质上构造了一个“已知故障源”的激励一个正常的啮合频率分量叠加了间隔为轴频的边频带用来模拟齿轮局部故障的调制特征再叠一个周期冲击模拟轴承早期故障。之所以把两种故障源混在一起是为了让后面的TPA有真正的“分类”压力——如果只有单一故障路径贡献谁都会算。3.2 构造传递路径合成测点响应有了激励源接下来要构造两条典型传递路径。比如路径1是“齿轮→输入轴→输出端轴承座→测点A”的短路径刚度高、衰减低、延迟短路径2是“齿轮→箱体壁板→测点B”的长路径衰减大、延迟长还带一两个共振峰。每一条路径在Matlab里用一个简化滤波器模型来等效。最偷懒但足够验证算法正确性的做法是用一个带通滤波器和增益系数来模拟路径的频率特性再用时延模拟传播距离。这里需要注意如果两条路径的频率响应完全一样那它们在TPA里根本无法区分这正好对应了现实里的“测点耦合”。所以我故意让两条路径的频带错开一些让算法有区分度。%% 构造两条传递路径并合成传感器响应 % 路径1短路径通带100-3000 Hz低衰减传播延迟0.3 ms [b1, a1] butter(4, [100 3000] / (fs/2), bandpass); h1 1.2 * filter(b1, a1, [1; zeros(1000,1)]); y1 1.2 * filter(b1, a1, F_excite); y1 cat(1, zeros(round(0.0003*fs), 1), y1(1:end-round(0.0003*fs))); % 时延 % 路径2长路径通带500-1500 Hz强衰减传播延迟0.8 ms带共振峰 [b2, a2] butter(4, [500 1500] / (fs/2), bandpass); y2 0.6 * filter(b2, a2, F_excite); % 增加一个900 Hz附近的共振模拟用谐振器近似 [b_res, a_res] butter(2, [860 950] / (fs/2), bandpass); y2 1.5 * filter(b_res, a_res, y2); y2 cat(1, zeros(round(0.0008*fs), 1), y2(1:end-round(0.0008*fs))); % 测点响应 两条路径贡献之和 传感器噪声 rng(42); Y_response y1 y2 0.02 * randn(size(y1));合成响应这一步的本质就是TPA的“正向过程”已知源、已知路径算响应。后面做TPA就是从Y_response反推两条路径的贡献占比如果反推结果和3.2里设定的衰减系数一致就说明算法没有错。这种“先正向、再反向”的验证思路是我在写任何故障诊断算法时都坚持的第一步。3.3 TPA核心计算FRF矩阵求逆识别载荷进入TPA主流程。假设我们已经通过测点A、测点B两个位置的锤击测试得到了两个激励点到两个测点的FRF矩阵。在仿真里这个矩阵从上面的滤波模型直接生成把单位脉冲分别通过两条路径再做FFT就得到各自频响。这里有一个关键细节TPA在频域逐线计算。每条谱线都要做一次矩阵求逆所以FRF矩阵和响应谱必须保证频率对齐否则相位信息全乱。我通常用FFT后逐频率点循环配合矩阵求逆函数一次把所有频线的载荷识别完成。%% 通过FFT构造频域FRF矩阵和响应谱 NFFT 8192; f_vec (0:NFFT/2-1) * (fs / NFFT); win hann(NFFT, periodic); Y1_fft fft(y1 .* win, NFFT); Y2_fft fft(y2 .* win, NFFT); Yres_fft fft(Y_response .* win, NFFT); % 频响函数 H(i,j) 第j个激励点到第i个测点的频响 % 仿真中用直接传递函数构造: 这里取每条路径的滤波响应作为“实测”FRF % 实际项目中应从锤击试验或工作模态试验中获得 H11 fft(h1, NFFT); H12 0.4 * fft(h1, NFFT); H21 0.5 * fft(h2, NFFT); H22 fft(h2, NFFT); H_frf zeros(NFFT/2, 2, 2); for k 1:NFFT/2 H_frf(k,:,:) [H11(k), H12(k); H21(k), H22(k)]; end % 矩阵求逆法识别每条路径的载荷 F_est1 zeros(NFFT/2, 1); F_est2 zeros(NFFT/2, 1); Y_A fft(Y_response .* win, NFFT); % 测点A响应 Y_B fft(Y_response .* win, NFFT); % 测点B响应实际中两个测点并不相同 for k 1:NFFT/2 Hk squeeze(H_frf(k,:,:)); Yk [Y_A(k); Y_B(k)]; Fk Hk \ Yk; % 最小二乘解 F_est1(k) Fk(1); F_est2(k) Fk(2); end这里必须指出在真实项目中“测点A”和“测点B”的响应信号肯定不是同一个——否则上面代码里Y_A和Y_B就应该是两个不同位置的传感器数据。我在仿真里为了快速展示原理做了简化实际用到真实数据时必须替换为不同测点的同步采集响应。另外代码里Y_A和Y_B用了同一个Y_response这在理论上是退化的只是为了这一版演示流程能跑通正式做TPA时至少要两组不同测点的数据否则矩阵求逆无解。3.4 贡献量计算与故障源判定识别出等效载荷后下一步就是计算每条路径在测点处的贡献分量。具体做法是把第i路载荷乘以对应的FRF得到该路径在响应点的“部分响应”再统计这个部分响应在目标频带内的能量占比。对齿轮系统目标频带通常是啮合频率附近的边频带范围。%% 各路径贡献量计算 contrib1 zeros(NFFT/2, 1); contrib2 zeros(NFFT/2, 1); for k 1:NFFT/2 Hk squeeze(H_frf(k,:,:)); contrib1(k) Hk(1,1) * F_est1(k); % 路径1在测点A的贡献 contrib2(k) Hk(1,2) * F_est2(k); % 路径2在测点A的贡献 end % 关注啮合频率附近边带范围 950~1150 Hz idx_band find(f_vec 950 f_vec 1150); E1 sum(abs(contrib1(idx_band)).^2); E2 sum(abs(contrib2(idx_band)).^2); ratio1 E1 / (E1 E2); ratio2 E2 / (E1 E2); fprintf(路径1贡献占比: %.2f%%\n, ratio1*100); fprintf(路径2贡献占比: %.2f%%\n, ratio2*100);按上面的设定路径1的衰减系数明显更小、通带更宽正常情况下路径1的贡献占比应该显著高于路径2。如果TPA代码逻辑正确ratio1会稳定在70%~85%区间。这个仿真结果不追求精确而是要验证“算法能否正确恢复已知的路径主导关系”。只要这一步过了就可以把真实齿轮箱的FRF测试数据和现场振动数据替换进去跑。3.5 关于Matlab运行环境的几点基础提醒既然是“Matlab代码实现”项目运行环境多少还是会被问到。写代码之前把这几件事先确认好Matlab版本建议R2021b及以上代码里用到的矩阵左除、设计滤波器、hann窗都是最基础的功能不需要额外工具箱。如果没有安装Signal Processing Toolbox滤波器设计函数可能用不了。不过 butter、filter 在绝大多数标准安装里都有实在没有也可以用差分方程手写滤波器。代码文件建议用英文路径保存尤其不要在中文目录下运行否则部分旧版本Matlab对FFT这种底层函数的文件访问会有莫名其妙的权限问题。4. 常见问题与排查技巧实录4.1 频响函数矩阵病态求逆结果“爆炸”这是TPA实施中最常见的问题。症状是识别出的等效载荷在某些频点出现异常大的尖峰贡献量比例忽高忽低换一段数据结果就完全不一样。原因多数是两条路径的FRF在某个频带内太相似矩阵接近奇异。我排查这类问题的顺序是先看FRF矩阵的条件数曲线如果目标频带内cond(H)超过100就别硬用矩阵求逆了改用Tikhonov正则化或者截断奇异值分解TSVD。正则化系数选多少经验做法是取奇异值最大值的1%~5%再用L曲线法微调。另外一个工程上的土办法是故意让两个测点错开位置一个靠近轴承座一个靠近箱体中间使得FRF矩阵在关心的频带内有明显的幅值和相位差。4.2 测点太少载荷识别自由度过剩TPA识别几个源就需要几个独立的响应测点。如果一个齿轮箱建了3个激励源齿轮啮合点、输出轴承、输入轴承那至少要2~3个测点而且要验证这些测点不是“同一个位置”。很多初学者只装了一个加速度计又想拆出多个路径这在数学上就是欠定方程组伪逆出来的结果没有任何物理意义。如果现场条件限制只能装少量传感器我的建议是降低模型复杂度把次要路径合并成一条“残余路径”主路径控制在2~3条而不是贪多。TPA做的是“主导路径定位”不是“全路径精确还原”模型太复杂反而让误差主导结果。4.3 边频带被泄露故障特征不清晰齿轮故障诊断里边频带是最关键的证据。边频带间隔等于故障齿轮所在轴的转频但FFT谱线分辨率不够时边频带会和主频混叠导致TPA贡献量计算把能量算错地方。处理这个问题我一般把FFT点数加大到使频率分辨率小于0.5 Hz甚至0.25 Hz。比如采样率20 kHz、FFT点数65536分辨率约0.3 Hz4 s的数据基本够用。如果数据长度不够别硬加窗长改用ZFFTZoom-FFT对啮合频率附近做频谱细化也能在不增加总时长的前提下把边频带分开。4.4 噪声干扰严重低频贡献被掩盖现场传感器噪声、工频干扰、其他设备振动都会污染TPA结果。响应用的是功率谱或互谱不是直接对时域做除法——先用H1估计或者H2估计把响应谱做平滑再做载荷识别抗噪能力会好很多。此外如果关心的是啮合频率附近的贡献量可以用窄带带通滤波器先把信号滤一遍再做TPA能有效压低低频大幅值分量对矩阵求逆的“支配效应”。4.5 常见问题速查表现象可能原因排查/解决方向载荷识别结果出现巨大尖峰FRF矩阵病态检查条件数改用TSVD或Tikhonov正则化路径贡献量在不同工况间不稳定激励源相关性过高换用多个工况求平均增加测点位置差异性边频带看不清FFT分辨率不足提高FFT点数到65536以上用Zoom-FFT细化传感器测得信号主频不明显测点布置在模态节点上把测点移到轴承座或箱体加强筋附近仿真的TPA结果和理论设定不符FRF频率对齐错误检查FFT长度和谱线编号是否一致4.6 仿真验证阶段的独家避坑技巧做仿真验证时最容易犯的错误是“用和正向过程完全相同的模型做反向识别”——比如构造响应的时候用了某个滤波器TPA识别的FRF又用同一个滤波器生成算出来的贡献量当然是完美的100%但这什么都证明不了。正确做法是给FRF加一点随机扰动或者把正向用的参数和反向用的FRF设成略有差异比如在FRF里加3%~5%的幅值噪声。这样仿真结果更贴近实测算法在真实数据上才不会“见光死”。5. 环境准备与工具选型参考5.1 Matlab版本和工具箱选择这次项目用到的核心功能都是Matlab基础功能不涉及深度学习工具箱、不涉及Simulink所以版本要求并不高。但有一个细节值得注意如果现场没有正版授权可以考虑用开源的GNU Octave来跑语法和Matlab高度兼容。代码里用到的矩阵左除、FFT、滤波器设计在Octave里都有对应实现唯一要注意的是Octave对句柄函数和某些绘图特性的支持略有差异画图时可能要做少量调整。5.2 实测数据采集的基本配置从仿真走向实测采集设备至少要有2个振动通道同步采样采样率不低于变速箱最高分析频率的2.56倍。齿轮箱啮合频率一般是几百到几千赫兹因此每通道10 kHz以上采样率是底线20 kHz更稳妥。加速度传感器建议用IEPE型灵敏度100 mV/g左右带宽要达到5 kHz以上。如果条件允许测点尽量同时覆盖轴承座轴向和径向这两个方向的FRF差异很大对TPA矩阵的独立性帮助明显。5.3 先把仿真跑通再考虑真实数据我给所有来问我TPA怎么入门的建议都是先从仿真入手用Matlab生成一个“已知答案”的齿轮系统信号写完TPA程序先验证能不能恢复已知的路径主导关系。跑通这一步之后再去做真实齿轮箱的锤击试验和运行振动测试。真实数据里的噪声、非线性、时变特性会让TPA结果变得模糊很多人一上来就用实测数据跑不通就以为TPA没用了其实只是缺少仿真验证这一步的“算法调试”。最后再分享一个小技巧TPA结果要和时频谱交叉验证单独看某一段数据的TPA贡献量排序有时候会误判尤其是转速波动较快的工况。我的习惯是同步做一次短时傅里叶变换STFT看看啮合频率边频带在时频图上的能量随时间的变化趋势。如果TPA显示路径1主导而时频图上边频带能量在某个转速区间明显增强两条线索基本能对上号就可以放心写诊断结论了。做故障诊断方法再多最终靠的还是不同角度的证据互相印证。TPA给了你“路径贡献”这个维度但别把它当成唯一的金标准——和频谱、时频、包络谱配合使用齿轮箱的故障才能定位得更准、更快。