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

齿轮箱一旦振动超标工程师最头疼的事情不是“振动大”而是说不清振动到底从哪个齿轮啮合点出来、经过哪条结构路径传到测点。同一个测点上的信号包含了电机转速波动、各级齿轮啮合激励、轴承故障冲击、箱体共振等多个源头再经过轴、轴承、箱体等不同路径的放大或衰减最后混在一起。只用频谱找峰值、做包络解调经常会出现“幅值最大的频段并不是故障源”的误导结果。传递路径分析TPA就是用来解决这个问题的把一个测点的总响应拆分成“哪些激励、走哪些路、贡献多少”的定量结果。这篇文章结合齿轮系统的特点用Matlab代码走一遍TPA的完整流程从原理到仿真验证再到实测迁移的注意点希望对做故障诊断的工程师和研究生有帮助。1. 齿轮箱故障诊断为什么要做传递路径分析1.1 齿轮箱振动信号是一个“多源叠加混音现场”齿轮箱可以说是典型的“多激励、多路径、单响应”系统。电机输入轴上有第一级小齿轮啮合到大齿轮后再通过轴承传递到下一级轴二级齿轮箱通常就有两对啮合副、四到六个轴承支撑点。每一个啮合副都是宽带激励源每经过一对轴承振动信号就被调制一次最终在箱体表面测点的加速度信号是所有源经过所有路径叠加后的结果。这里的“叠加”不是简单的幅值相加而是带有相位关系的矢量叠加。假设某个测点在某频率下收到了啮合频率的振动这条路径上经过了轴、轴承、箱体每次经过一个结构界面都会发生幅值变化和相位移动。如果另一条路径上的同频信号相位相反两条路径叠加后反而会抵消。这就是为什么现场经常出现“某个齿轮坏得很明显但箱体表面测点上的特征频率幅值反而不高”的怪象。只用窄带频谱或者包络谱去看单一测点本质上是在看“混音现场”的最终录音不同声部的贡献被混在一起根本无法还原哪个乐器在哪个时间弹了什么音。1.2 TPA解决的三个问题源、路径、贡献传递路径分析本质上是一个“系统辨识贡献量分解”的过程它的任务被拆成三个问题源是谁齿轮啮合激励、轴承故障冲击、轴不对中引起的周期力这些都是候选源。路径是什么源产生的力通过哪些物理传力通道到达接收点常见路径包括齿轮→主动轴→主动轴轴承→箱体→测点或者齿轮→从动轴→从动轴轴承→箱体→测点还有通过空气传播的声学路径。贡献有多大每条路径在接收点上产生的响应占总体响应的百分比或者说按照幅值和相位把总响应重新合成为各路径分量。TPA的输出不是一张频谱图而是一张“贡献量排行表”。故障诊断人员拿到这个表以后就可以直接锁定“这个测点的过量振动主要来自第二条路径”再结合齿轮箱结构图把第二条路径上的中间环节拆开来检查。这种方法比单纯对照故障特征频率库去猜要可靠得多尤其在多级齿轮箱这种多源强耦合场景下优势非常明显。1.3 齿轮系统中传递路径的划分思路齿轮箱的路径划分没有唯一答案需要根据分析目的调整。一般来说路径划分应当满足两个原则物理上可解释数学上可叠加。物理上可解释意味着每一条路径都能在齿轮箱结构中找到对应传力链数学上可叠加意味着把一个源沿一条完整路径拆成多段时各段的传递函数相乘等于整条路径的传递函数。以一个二级平行轴齿轮箱为例我通常这样划分路径路径编号路径描述主要传力环节P1高速级啮合激励至箱体测点A小齿轮→高速轴→高速轴轴承→箱体P2高速级啮合激励至箱体测点A另一侧大齿轮→中间轴→中间轴轴承→箱体P3低速级啮合激励至箱体测点A中间轴齿轮→低速轴→低速轴轴承→箱体P4输入轴不平衡力至箱体测点A输入端联轴器→高速轴→高速轴轴承→箱体P5空气声路径齿轮啮合点→内部空气→箱体表面这里有一个容易被忽略的细节同一个源走同一条轴但是左右两侧各有一个轴承那么应该算两条路径还是一条路径从纯数学角度可以合并成一条等效传递路径但工程上建议保留为两条。原因是左右轴承的支撑刚度和故障状态往往不同实测频响函数差异很大。如果把两个轴承合并成一条路径反演出来的路径力会变成某种“平均力”反而丢失了哪一侧轴承劣化的信息。关于空气声路径在齿轮箱 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开发时调通一套算法比较复杂建议先用仿真数据验证算法正确性再换实测数据。仿真思路是先定义“真实的”激励力和路径传递函数合成接收点响应然后用合成响应反推路径力看反推力能否还原真实激励最后做贡献量分解。这个流程能有效隔离算法问题和数据质量问题。模拟齿轮啮合激励时我会构造啮合频率及其谐波并用相位调制模拟齿面故障产生的边带。齿数取 Z120、Z260主动轴转速为 30 r/s啮合频率计算为 30×20600 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/21); 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_pathlambda 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 流程和数据库才是在实际项目中真正能依赖的工具。