基于Matlab的齿轮箱传递路径分析与故障诊断贡献量分解

齿轮箱一旦振动超标，工程师最头疼的事情不是“振动大”，而是说不清振动到底从哪个齿轮啮合点出来、经过哪条结构路径传到测点。同一个测点上的信号，包含了电机转速波动、各级齿轮啮合激励、轴承故障冲击、箱体共振等多个源头，再经过轴、轴承、箱体等不同路径的放大或衰减，最后混在一起。只用频谱找峰值、做包络解调，经常会出现“幅值最大的频段并不是故障源”的误导结果。传递路径分析（TPA）就是用来解决这个问题的：把一个测点的总响应拆分成“哪些激励、走哪些路、贡献多少”的定量结果。这篇文章结合齿轮系统的特点，用Matlab代码走一遍TPA的完整流程，从原理到仿真验证再到实测迁移的注意点，希望对做故障诊断的工程师和研究生有帮助。

1. 齿轮箱故障诊断为什么要做传递路径分析

1.1 齿轮箱振动信号是一个“多源叠加混音现场”

齿轮箱可以说是典型的“多激励、多路径、单响应”系统。电机输入轴上有第一级小齿轮，啮合到大齿轮后再通过轴承传递到下一级轴，二级齿轮箱通常就有两对啮合副、四到六个轴承支撑点。每一个啮合副都是宽带激励源，每经过一对轴承，振动信号就被调制一次，最终在箱体表面测点的加速度信号是所有源经过所有路径叠加后的结果。

这里的“叠加”不是简单的幅值相加，而是带有相位关系的矢量叠加。假设某个测点在某频率下收到了啮合频率的振动，这条路径上经过了轴、轴承、箱体，每次经过一个结构界面都会发生幅值变化和相位移动。如果另一条路径上的同频信号相位相反，两条路径叠加后反而会抵消。这就是为什么现场经常出现“某个齿轮坏得很明显，但箱体表面测点上的特征频率幅值反而不高”的怪象。只用窄带频谱或者包络谱去看单一测点，本质上是在看“混音现场”的最终录音，不同声部的贡献被混在一起，根本无法还原哪个乐器在哪个时间弹了什么音。

1.2 TPA解决的三个问题：源、路径、贡献

传递路径分析本质上是一个“系统辨识+贡献量分解”的过程，它的任务被拆成三个问题：

• 源是谁：齿轮啮合激励、轴承故障冲击、轴不对中引起的周期力，这些都是候选源。
• 路径是什么：源产生的力通过哪些物理传力通道到达接收点，常见路径包括齿轮→主动轴→主动轴轴承→箱体→测点，或者齿轮→从动轴→从动轴轴承→箱体→测点，还有通过空气传播的声学路径。
• 贡献有多大：每条路径在接收点上产生的响应占总体响应的百分比，或者说按照幅值和相位，把总响应重新合成为各路径分量。

TPA的输出不是一张频谱图，而是一张“贡献量排行表”。故障诊断人员拿到这个表以后，就可以直接锁定“这个测点的过量振动主要来自第二条路径”，再结合齿轮箱结构图，把第二条路径上的中间环节拆开来检查。这种方法比单纯对照故障特征频率库去猜要可靠得多，尤其在多级齿轮箱这种多源强耦合场景下，优势非常明显。

1.3 齿轮系统中传递路径的划分思路

齿轮箱的路径划分没有唯一答案，需要根据分析目的调整。一般来说，路径划分应当满足两个原则：物理上可解释，数学上可叠加。物理上可解释意味着每一条路径都能在齿轮箱结构中找到对应传力链；数学上可叠加意味着把一个源沿一条完整路径拆成多段时，各段的传递函数相乘等于整条路径的传递函数。

以一个二级平行轴齿轮箱为例，我通常这样划分路径：

这里有一个容易被忽略的细节：同一个源走同一条轴，但是左右两侧各有一个轴承，那么应该算两条路径还是一条路径？从纯数学角度可以合并成一条等效传递路径，但工程上建议保留为两条。原因是左右轴承的支撑刚度和故障状态往往不同，实测频响函数差异很大。如果把两个轴承合并成一条路径，反演出来的路径力会变成某种“平均力”，反而丢失了哪一侧轴承劣化的信息。

关于空气声路径，在齿轮箱 TPA 中通常贡献很小，因为箱体钢板对中高频空气声有隔离作用，结构传声占据主导。但如果是开式齿轮传动或者带大尺寸散热孔的箱体，空气声路径就不能忽略，可以在后期作为误差项分析。处理方式是单独测箱体表面振动到麦克风位置的声学传递函数，然后并到总贡献里做合成验证。

2. 传递路径分析的数学模型与参数选择

2.1 线性叠加模型

传统TPA的理论基础是线性时不变系统的叠加原理。在频域里，接收点响应可以写成：

Y(ω) = Σ Hi(ω) × Fi(ω)

其中 Y 是接收点响应（加速度、速度或力均可），Hi 是第 i 条路径的频响函数，Fi 是第 i 条路径上的激励力，ω 是角频率。这里所有量都是复数，包含幅值和相位。

这个模型隐含三个假设。第一，系统在工作状态下的动态特性与离线测量频响函数时的状态一致。齿轮箱的滚动轴承、齿轮啮合接触都会随转速和负荷发生变化，严格来说运行状态和静态测量状态存在差异，这是误差的主要来源之一。第二，各路径之间互不干涉，即一条路径上的激励不会改变另一条路径的传递特性。第三，源与结构之间的耦合作用可以忽略，源被视为理想力源。齿轮系统中齿轮啮合激励比较接近力源特性，而电机不平衡产生的激励则有一定刚度耦合，处理时要注意。

实际做TPA时不需要把每一个细微环节都精确建模，因为“预测”不是目的，“拆分”才是。只要模型能把总响应拆成若干有物理意义的子响应，并且在合成回测时误差在可接受范围内，就可以用于定位主导路径了。

2.2 频响函数的获取与质量校验

频响函数是 TPA 的核心输入，它定义了力在路径入口处作用时，接收点产生的响应。获取方式有两种：

第一种是锤击法，使用力锤敲击轴承座或者路径入口位置，同时采集力信号和接收点加速度信号，用互谱除以自谱得到频响函数估计。这种方法操作简单、速度快，适合齿轮箱这种体积不大的结构，但激励能量有限，对于大型箱体可能出现低频信噪比不足。

第二种是激振器法，将激振器通过弹性杆连接到结构上，输入随机或扫频信号。激振器法可以施加更大、更可控的激励力，适合体积大、阻尼高的结构，但安装比较复杂，还需要注意激振器顶杆与结构之间的连接刚度，避免引入额外质量效应。

Matlab里可以直接用 tfestimate 函数做频响估计。我整理的常用参数如下：

nfft = 4096; win = hann(nfft, 'periodic'); overlap = round(nfft * 0.75); [H1, f] = tfestimate(force, acc, win, overlap, nfft, fs); [Cxy, f] = mscohere(force, acc, win, overlap, nfft, fs);

频响函数估计好后，必须检查相干系数。相干系数低于0.9的频率点，说明响应中有大量噪声或者出现非线性，这些点上的频响函数不可信，通常需要剔除或者在后续加权中降低权重。如果实测过程中某个频段的相干性一直上不去，我建议先检查激励是否足够（加大锤击力度或激振器输出）、传感器是否松动，而不是急着换更复杂的算法。

2.3 路径力估计的三种方法

频响函数可以测，但路径力在工作状态下很难直接测。路径力本质上是轴承作用在箱体上的动态力，总不能把轴承拆了装力传感器再运行。工程上常用三种间接估计方法：

直接测量法：在试制样品上把路径环节串联力传感器，直接测量力。这个方法最准确，但在齿轮箱上几乎不可行，因为轴承和箱体之间没有安装力传感器的空间，而且加了传感器就改变了接触刚度和预紧状态。一般在台架验证阶段才会用。

动态刚度法：利用路径两侧的相对位移或者相对加速度乘以连接件的动态刚度来估算力。公式是 F = Kd × (U_s - U_r)，其中 Kd 是轴承或支撑结构的动态刚度，U_s 和 U_r 分别是源侧和接收侧的位移。难点在于动态刚度通常不是常数，轴承刚度随转速和负载变化明显，需要查厂家数据或者专门标定。

矩阵反演法：在接收点附近布置多个指示点加速度传感器，建立指示点响应与路径力之间的转送矩阵 T，然后通过矩阵求逆得到路径力估计。这是工程上最常用的方法。基本过程是：

a_ref = T × F

其中 a_ref 是 m 个指示点的加速度向量，F 是 n 条路径的力向量，T 是 m×n 的转送矩阵，矩阵每个元素通过离线测量得到。当 m ≥ n 时，可以求最小二乘解：

F_est = (T'×T)^(-1) × T' × a_ref

实际操作中转送矩阵常常存在数值病态，我会在代码里加一个 Tikhonov 正则化项。这部分后面展开。

2.4 贡献量计算中的相位与矢量叠加

各路径贡献量本质上是复数矢量。假设第 i 条路径的贡献量是 C_i = Hi × Fi，那么接收点总响应的重构值是 Σ C_i。因为 C_i 是矢量，合成时不是简单的幅值加减。

这里有一个典型的错误做法：只取各路径贡献的幅值排行，把幅值最大的路径当作主导路径。但在相位相反的情况下，两条幅值接近的路径可能相互抵消，此时单看幅值会得出完全错误的结论。正确做法是计算各路径贡献矢量的“合成效果”，也就是将每条路径的贡献投影到总响应矢量方向上进行排序。Matlab代码里可以这样做：

total_recon = sum(contrib, 1); proj = real(contrib .* conj(total_recon)); % 投影到总响应方向的实部

这个投影值才真正反映了该路径对最终响应的“推动”作用。

3. Matlab代码实现：构建可验证的齿轮箱TPA流程

3.1 仿真数据结构设计

直接上实测数据进行TPA开发时，调通一套算法比较复杂，建议先用仿真数据验证算法正确性，再换实测数据。仿真思路是：先定义“真实的”激励力和路径传递函数，合成接收点响应，然后用合成响应反推路径力，看反推力能否还原真实激励，最后做贡献量分解。这个流程能有效隔离算法问题和数据质量问题。

模拟齿轮啮合激励时，我会构造啮合频率及其谐波，并用相位调制模拟齿面故障产生的边带。齿数取 Z1=20、Z2=60，主动轴转速为 30 r/s，啮合频率计算为 30×20=600 Hz。边带间隔等于主动轴转频 30 Hz，这是齿轮故障特征的常见表现形式。具体生成代码如下：

fs = 8192; N = fs * 4; t = (0:N-1)' / fs; fr = 30; % 主动轴轴频 30 Hz z1 = 20; fm = fr * z1; % 啮合频率 600 Hz f_source = zeros(N, 1); for k = 1:4 amp = 1 / k^1.5; % 谐波幅值递减 mod_depth = k * 0.8; % 调制指数，模拟齿轮故障 f_source = f_source + amp * sin(2*pi*k*fm*t + ... 2*mod_depth*sin(2*pi*fr*t) + k*pi/6); end

这里不要把调制指数设置过大，否则边带会跨到其他频段，仿真看起来更像严重故障，但不利于算法验证。我建议先设置小一点，保证主要能量集中在啮合频率和二阶谐波附近。

路径传递函数我用二阶系统叠加来模拟结构共振峰。每个轴承路径在 400Hz、1800Hz、4200Hz 附近设置共振峰，阻尼比取 0.03 到 0.06，不同路径的参数错开，以体现不同传力路径的频率选择性。

function H = poleSys(f, f0, zeta, gain) w = f / f0; H = gain ./ (1 - w.^2 + 1j * 2 * zeta * w); end f_axis = (0:fs/2) / fs * fs * 0.5; % 简化示意，实际按 NFFT 生成 H_p1 = poleSys(f_axis, 400, 0.04, 0.8) + poleSys(f_axis, 1800, 0.05, 0.5);

合成接收点响应时按频域相乘再逆变换回时域。最终得到测点信号、参考指示点信号，以及用于反演的真实路径力。

3.2 FRF与路径力反演核心代码

仿真数据准备完之后，第一步是用 tfestimate 估计各条路径的频响函数。这里有一个关键点：每条路径的频响函数必须在“路径入口”到“接收点”之间测量，而不是直接用 CAE 计算值替代。代码层面可以写一个循环，逐路径估计并存储：

nfft = 2048; win = hann(nfft, 'periodic'); noverlap = round(nfft * 0.75); H_path = zeros(n_path, nfft/2+1); for k = 1:n_path [H_path(k,:), f] = tfestimate(... force_path(:,k), acc_measure, win, noverlap, nfft, fs); [Cxy, ~] = mscohere(... force_path(:,k), acc_measure, win, noverlap, nfft, fs); H_path(k, Cxy < 0.85) = NaN; % 相干差的位置置空 end

反演路径力时，需要先构建转送矩阵。转送矩阵的元素表示第 j 条路径入口施加一个单位力时，第 i 个指示点的响应。这个矩阵同样通过实验测量或仿真获得。假设有 n_path 条路径、n_ref 个指示点，矩阵维度是 n_ref × n_path：

lambda = 1e-3; T_reg = T' * T + lambda * eye(n_path); F_est = T_reg \ (T' * a_ref_fft);

这里的 a_ref_fft 是实测指示点响应的频谱向量。正则化参数 lambda 不能盲目取，我一般从 1e-6 开始，以 10 倍步长递增，观察 F_est 的稳定性，选一个既保持幅值不过度收缩又能控制振荡的区间。手写循环自动搜索也可以，但工程上先手动看几个频点更快。

3.3 路径贡献分解和结果可视化

得到路径力估计 F_est 后，各路径贡献量就是路径频响函数与路径力的逐频点乘积：

contrib = zeros(n_path, length(f)); for k = 1:n_path contrib(k,:) = H_path(k,:) .* F_est(k,:); end total_recon = sum(contrib, 1);

注意这里 H_path 和 F_est 是逐频点的行向量，.* 是数组点乘，乘法结果是复数矢量，幅值和相位都保留。

可视化输出建议做三个图。第一张是总响应的实测频谱与重组频谱对比，如果算法正确，两条曲线基本重合，误差在 3dB 以内。第二张是目标频段内各路径贡献量柱状图，显示投影量或幅值。第三张是各路径贡献的相位谱图，用于检查相位抵消情况。绘制投影贡献的代码如下：

[~, idx] = min(abs(f - fm)); % 取啮合频率附近的频点 proj_contrib = real(contrib(:, idx) .* conj(total_recon(idx))); bar(proj_contrib); set(gca, 'XTickLabel', {'P1', 'P2', 'P3', 'P4', 'P5'}); ylabel('Projected contribution');

在做仿真验证时，我会额外画一张“真实路径力 vs 反演路径力”的对比图。这是验证反演算法内在一致性的关键，至少要在目标频段达到单频点误差小于 20%。

3.4 迁移到实测数据的几个关键改造点

仿真验证完之后换实测数据，主要有几个坑。第一个坑是数据采样率不一样，实测可能 51200 Hz，仿真用的 8192 Hz，需要重采样到统一分析频率范围，或者直接按实测数据重新计算 NFFT。第二个坑是转速存在波动，齿轮箱启动和负载变化时轴频会漂移，单纯用 Hz 频率分析会出现峰值模糊。这时需要引入转速计信号，做阶次跟踪重采样。第三个坑是实测的频响函数一般只测到线性平均的传递函数，而运行工况下的轴承载荷和温度会让实际传递特性偏移，尤其在共振峰附近，偏差最明显。

如果实测数据的通道数有限，无法在每条路径入口都布置传感器，我建议至少保证指示点数量不小于路径数量。指示点要布置在刚度较强的箱体筋位或轴承座附近，而不是箱体盖板中央这样柔性大的位置。否则反演矩阵条件数会很差，路径力估计一片发散。

从仿真转向实测时，代码层面其实改动不大，主要是数据对齐和异常值处理。我在实际项目中会封装一个数据加载函数，自动读取 TDMS 或 DAT 格式的时域数据，然后统一做抗混叠滤波和下采样，确保所有通道在同一采样率下工作，以免后续频域分析时频点对不上。

4. 常见问题与排查技巧实录

4.1 矩阵病态导致路径力反演发散

这是 TPA 里最容易踩的坑。转送矩阵反演时，如果两条路径的传递特性过于相似，矩阵就接近奇异，微小噪声会让路径力估计结果产生巨大振荡。现象是 F_est 的值在相邻频点忽高忽低，重构总响应曲线看起来还行，但单独看路径贡献量毫无物理意义。

解决办法有三个方向：一是增加指示点数量，利用冗余信息把条件数压下来；二是精简路径数量，把物理上接近、传递函数相似的路径合并；三是加 Tikhonov 正则化，牺牲一点拟合精度换取稳定性。我在实际项目中经常三种手段并用，正则化系数取 1e-3 左右，然后用贡献量的平滑度来检验是否过正则化。

4.2 FRF测量相干性低

锤击法测量时，如果相干性在某个频段长期低于 0.85，通常是激励能量不足或者结构存在非线性。我的排查顺序是先看力谱是否在目标频段有足够平坦的覆盖，再看传感器固定是否可靠，最后才怀疑结构非线性。力锤的锤头材料不是随便选的，硬锤头激励高频好但低频能量弱，软锤头低频激励好但高频不足。做 TPA 一般需要分别用不同硬度的锤头测两次，然后按频段拼接，我强烈建议不要图省事只测一次。

4.3 相位同步与通道校准问题

多通道数据采集时，各通道之间的相位差会直接影响复数矢量合成。如果采集卡本身具备同步采样功能，这个问题不大。但有些模块化采集系统使用了扫描方式扫描各个通道，通道间会有固定时延，在 600Hz 以上频段，这个时延造成的相位误差就非常明显。解决办法是在实验前做通道校准，用标准信号同时输入所有通道，计算各通道间的传递函数并做相位补偿。Matlab 里可以用一次简单的正弦扫频信号完成这个校准。

4.4 贡献量合不上、重构误差大

当总响应的实测值与 TPA 重构值差异超过 3dB 时，第一反应应该是路径不完整，遗漏了某一条重要传力路径，而不是频响函数不准。有经验的工程师会让重构误差反过来指导路径模型修正：哪些频段误差大，就去找这个频段对应的结构模态，看是否有一个轴承点或者联轴器点没有纳入模型。另一个容易被遗漏的是“接地路径”，比如齿轮箱安装底脚通过螺栓传到基础平台再弹回来，这条路径在低频段往往贡献不可忽视。

4.5 问题速查表

根据我个人这几年的实操体会，TPA 运行马马虎虎，算出一个漂亮的结果并不难，真正考验人的是对结果的解释和验证。每次完成贡献量分解后，我都会挑选一到两条判定为主导的路径，做一次“阻断验证”，比如临时调整对应路径上的预紧力或者更换该路径的衬垫，看测点响应是否按预测的方向发生变化。这种实验与计算互相支撑的做法，比反复调参数更有说服力。代码跑通只是第一步，建立一套适合自己设备的 TPA 流程和数据库，才是在实际项目中真正能依赖的工具。