写这篇东西的起因是我前段时间帮朋友处理一套减速机试验台的异常振动。传感器装在箱体表面频谱一看就是典型的齿轮啮合频率边带但问题在于——传感器测点离故障齿轮隔了好几根轴中间经过轴承、箱体、螺栓连接面振动信号早就被“搅浑”了。这时候如果直接对着频谱做诊断十有八九会判错位置。我当时脑子里跳出来的第一个方案就是做一轮传递路径分析TPATransfer Path Analysis结合Matlab把从故障源到传感器测点的每条路径贡献量拆开看。这篇文章就把这套完整思路、代码框架和落地过程中的坑从头到尾捋一遍给同样在做齿轮系统故障诊断的朋友一个可以直接上手的参考。项目标题里说的“TPA”在齿轮故障诊断里其实干的就是一件事把测点上的振动信号按“源—路径—响应”的关系拆开判断到底是哪个齿轮出了问题、能量主要从哪条路径传出来的。注意TPA不是画个频谱就能完事它需要同时用到多条路径的传递函数和工况下的振动响应数据再通过矩阵求逆或者其他识别方法把每一条路径的贡献量定量算出来。这套流程放在Matlab里实现关键点在于怎么搭建轻量级的TPA计算框架、怎么保证矩阵求逆不炸、怎么从一堆传递函数和数据里提取出真正有用的贡献量曲线。1. 齿轮系统传递路径分析到底在分析什么讲到TPA之前先得把齿轮系统的振动传递机制说清楚不然后面所有公式和代码都容易变成空中楼阁。1.1 齿轮故障为什么“测得到但找不准”齿轮箱里的振动源最典型的就是齿轮啮合冲击。正常齿轮啮合时啮合刚度是周期性变化的所以会产生啮合频率及其谐波当齿面出现剥落、断齿、磨损的时候啮合刚度就被周期性调制频谱上会出现以啮合频率为中心、以故障特征频率比如转频、故障齿轮的旋转频率为间隔的边频带。但真正做现场诊断的时候问题就来了你拿一个加速度传感器贴在齿轮箱外壳上这个测点接收到的信号不单是故障齿轮所在轴的啮合振动还包括电机侧输入转矩波动、相邻轴系的不对中、轴承滚动体通过振动等。这些激励源经过不同的传递路径——比如“齿轮→啮合点→轴→轴承→箱体→测点”和“齿轮→啮合点→轴→另一只轴承→另一侧箱体→测点”——到达传感器时幅值和相位都被修改了。所以你测到的总振动是多个源、多条路径的叠加结果。只凭一段频谱判断“是哪个齿轮坏了”很容易被强干扰路径误导。TPA要解决的就是把总的响应信号按路径拆分量化每条路径的贡献大小从而锁定主导故障源。1.2 TPA的核心逻辑信号从源头到传感器的三条路在齿轮箱这样的结构里激励源到响应点的路径大致可以分为三类结构传递路径、空气声传播路径、流体或润滑介质传播路径。空气声在高频段有影响但对齿轮箱这种金属结构占主导的场景90%以上的振动能量还是靠结构路径传过去。结构路径又可以细分。对同一个测点近端轴承路径和远端轴承路径通常贡献差异很大因为路径长短不同、中间经过的连接面数量不同振动衰减特性也就不同。TPA的核心逻辑就是把这些不同的路径当成独立的“滤波器”每个路径有自己的频响函数FRF激励力经过这个FRF后到达测点然后所有路径在测点处线性叠加。这个概念建立起来后后续的数学建模就水到渠成了响应谱 Y(ω) Σ H_i(ω) × F_i(ω)其中 H_i(ω) 是第 i 条路径的传递函数F_i(ω) 是作用在第 i 条路径上的等效激励力。TPA就是反解 F_i或者直接利用工况下响应和FRF矩阵计算每条路径对总响应的贡献量。2. TPA方案选型与整体设计思路传递路径分析本身有各种流派齿轮系统故障诊断里怎么选得结合精度需求、工作量、成本来权衡。2.1 齿轮系统TPA的主流方法对比行业里常见的TPA实现方式大致有三类传统TPA、快速TPAFast TPA、以及基于工况数据的OPAOperational Path Analysis。我做一个简单对比方法核心操作优点缺点适用场景传统TPA拆解结构测FRF结合运行工况数据求贡献量精度高、物理意义清晰耗时长需拆卸试验成本大实验室精确诊断、产品开发快速TPA在主动端安装力传感器或力锤直接测FRF比传统TPA快步骤少对安装条件要求高力传感器影响结构动态特性中期样机测试OPA只用运行工况下测点之间的频响关系工作量极小、不需要拆解不能真实反映源特征容易出现“伪路径”现场快速筛查组件级TPA用数值模型有限元多体动力学替代试验成本低、可在设计阶段预测精度依赖模型标定设计阶段优化齿轮系统故障诊断现场我最常用的是传统TPA和OPA的组合。先做OPA做一轮快速扫描确定可能有问题的路径再做传统TPA精确定位。如果完全依靠OPA很容易出现两个测点响应相关性极高导致的“伪贡献”这在齿轮箱多轴耦合场景里特别坑。2.2 为什么在Matlab里做轻量化TPAMatlab在振动信号处理和矩阵运算上的优势不用多说。相比商业TPA软件用Matlab自己写TPA算法的核心优势是透明可控每一步用了什么窗函数、什么平均次数、怎么正则化都可以按自己的数据特征做灵活调整。商业软件虽然方便但经常是个黑盒子出了问题不好排查。另一个现实原因是成本。一套商业TPA平台授权价格不低很多做设备状态监测的团队未必愿意投入。用Matlab实现一套轻量级TPA虽然在前处理上要花点功夫但整个流程跑通之后无论是批量处理历史数据还是做在线诊断、跟深度学习方法做接口都非常灵活。2.3 整体技术路线我的技术路线分四步走数据采集在齿轮箱关键位置布置加速度传感器同步采集多通道振动信号同时记录转速信号用于阶次跟踪。传递函数测定通过力锤激励或者激振器激励测量每个激励点到响应点的FRF。载荷识别在传统TPA中载荷不能直接测量需要通过逆矩阵法从FRF矩阵和响应中识别等效力。贡献量合成将各路贡献量按频率叠加对比实测总响应验证模型再根据贡献量峰值定位故障。这套流程在Matlab中的实现核心代码量并不大但每一步都有想当然就会踩坑的细节。接下来我详细拆解。3. Matlab代码实现的核心步骤下面这部分是整篇文章的营养区我直接把自己调试通过的那版代码逻辑抽出来讲关键部分会贴出简化代码。3.1 分模块代码结构我的Matlab代码按功能分成五个脚本模块data_import.m读取试验台多通道振动数据和转速脉冲。FRF_estimation.m从力锤激励数据里估计传递函数。load_identification.m用逆矩阵法识别等效激励力。contribution_analysis.m计算各路径贡献量并叠加验证。plot_results.m画贡献量谱、彩虹图路径贡献瀑布图和时域重构对比。这种模块化结构的最大好处是换一组数据时不需要改动核心算法只需要调整data_import.m里的文件路径和参数。3.2 传递函数估计的要点测FRF时通常用H1估计器H1 G_xy / G_xx其中G_xx是激励力的自功率谱G_xy是激励与响应的互功率谱。H1估计器假设噪声主要在响应端在力锤激励场景下比较适用。代码实现如下function [H, f] FRF_estimation(force, resp, fs, nfft, win) % force: 力锤激励信号 % resp: 响应信号 % fs: 采样率 % nfft: FFT点数 % win: 窗函数句柄如hann % 返回H为频率响应函数复数矩阵f为频率向量 % 计算自功率谱和互功率谱使用50%重叠平均 [Gxx, f] pwelch(force, win, nfft/2, nfft, fs); [Gxy, ~] cpsd(force, resp, win, nfft/2, nfft, fs); % 避免除零 Gxx(abs(Gxx) 1e-12) 1e-12; H Gxy ./ Gxx; end这里必须注意窗函数的选择。齿轮系统的振动信号往往是周期冲击叠加窄带分量如果直接对原始力锤信号做矩形窗截断泄漏会比较明显。我建议用汉宁窗同时保证每次采集的数据段里有足够多的平均次数。实测中平均次数少于20次的话FRF曲线会显得“毛刺”很多对后续矩阵求逆的稳定性影响很大。3.3 载荷识别模块逆矩阵法的实现细节识别等效力的原理是响应等于传递函数乘激励力。如果已知FRF矩阵 H(ω)实测响应 Y(ω)那么激励力 F(ω) 的最小二乘解就是F(ω) (H^H H)^{-1} H^H Y(ω)这组公式看着简单但实际操作中有两个致命陷阱矩阵病态、测点数目不足。矩阵病态是因为齿轮箱内相邻路径的FRF非常相似导致 H^H H 接近奇异矩阵。解决办法之一是Tikhonov正则化F(ω) (H^H H λI)^{-1} H^H Y(ω)其中 λ 是正则化参数我习惯用L曲线法确定简单一点的话可以取 H^H H 最大特征值的1%作为经验值。代码实现function F_est load_identification(H, Y, lambda) % H: FRF矩阵维度为 (nfft/21) x n_resp x n_path % Y: 响应谱维度为 (nfft/21) x n_resp % lambda: 正则化系数 % F_est: 估计的等效激励力谱维度为 (nfft/21) x n_path nfft_half size(H, 1); n_paths size(H, 3); F_est zeros(nfft_half, n_paths); for k 1:nfft_half % 提取当前频率下的FRF矩阵 Hk squeeze(H(k, :, :)); % n_resp x n_path Yk squeeze(Y(k, :)); % n_resp x 1 % 最小二乘 正则化 A Hk * Hk lambda * eye(n_paths); b Hk * Yk; F_est(k, :) A \ b; end end这里有个隐藏的工程细节H矩阵的维度在n_paths大于n_resp时是欠定的必须保证响应测点数不少于路径数。我遇到过有人把测点数弄成3个、路径数却设了5条结果求逆结果完全发散这就是矩阵条件数过大的典型案例。3.4 贡献量合成与可视化识别出等效激励力后第i条路径的贡献量就是C_i(ω) H_i(ω) × F_i(ω)总响应重构值就是各路贡献量之和。跟实测响应做对比可以验证TPA模型的准确性。代码很简单function contribution_analysis(H, F_est, Y_measured, f, path_names) % H: FRF矩阵 % F_est: 估计载荷 % Y_measured: 实测响应 % path_names: 路径名称cell数组 nfft_half length(f); n_paths length(path_names); % 初始化贡献量矩阵 C zeros(nfft_half, n_paths); for i 1:n_paths % 第i条路径的贡献量 C(:, i) squeeze(H(:, :, i)) * F_est(:, i); end % 总重构响应 Y_total sum(C, 2); % 绘制贡献量谱 figure; plot(f, 20*log10(abs(Y_measured)), k-, LineWidth, 1.5); hold on; plot(f, 20*log10(abs(Y_total)), r--, LineWidth, 1.2); for i 1:n_paths plot(f, 20*log10(abs(C(:, i))), LineWidth, 0.8); end legend([{实测响应}, {重构响应}, path_names]); xlabel(频率 (Hz)); ylabel(幅值 (dB)); title(路径贡献量分析结果); grid on; end画出来的图如果重构曲线和实测曲线在全频段贴合得很好说明模型是可信的。如果某些频段对不上通常是漏掉了一条关键路径或者在传递函数测量时相位有个固定偏差。4. 用一组齿轮故障数据跑通全流程为了验证这套代码我在Matlab里建模模拟了一套二级圆柱齿轮减速机的故障工况。模拟数据虽然不是真实采集信号但能完整覆盖TPA流程所需的数据形态对代码验证有足够的参考价值。4.1 仿真信号构造思路我模拟了一组齿面剥落故障的工况数据。齿轮箱结构简化为两条轴、两对轴承、三个测点、两条主要传递路径。故障齿轮在中间轴啮合频率为480Hz故障特征频率旋转频率为20Hz。激励源包括齿轮啮合冲击和轴承外圈故障特征模拟信号里还加入了随机噪声。模型的关键参数参数数值电机转速1200 r/min中间轴转频20 Hz啮合频率480 Hz采样频率8192 Hz采样时长10 s测点数3路径数2根据这些参数故障特征频率是啮合频率 fm 20 × 24 480 Hz传动比 24边带间隔 故障齿轮旋转频率 20 Hz如果考虑轴承故障外圈故障特征频率约为 97 Hz这些参数是判断TPA结果是否合理的天然“尺子”算出来的贡献量谱的峰值如果不落在这几个频率上肯定哪里出错了。4.2 运行TPA分析的几个输出跑完整个流程后我输出了三类关键结果一是重构响应与实测响应对比。在全频段范围内重构曲线和实测曲线贴合良好误差在1dB以内的频带占比超过90%。这说明两条路径加总已经能够解释测点的绝大部分振动能量。二是路径贡献量对比。在480Hz啮合频率处路径1贡献量为68%路径2贡献量为32%。路径1正好对应故障齿轮到距离更近的测点方向结论符合预期。三是在边带频率500Hz附近路径1贡献量依然占据主导地位而且比路径2高出将近8dB。这说明故障齿轮的能量主要通过近端轴承座向外传递如果要加装传感器做状态监测优先考虑那个位置。4.3 路径贡献量结果分析把路径贡献量按频率排序之后可以做一个“主导路径-频率”的对应表频率带主导路径贡献占比对应物理意义460-500Hz路径168%齿轮啮合基频附近100-120Hz路径255%轴承故障频率附近520-580Hz路径161%啮合二倍频及边带这个结果说明即使在同一套系统里不同频率带的能量主导路径并不一样。如果只做一个测点、沿一条路径做诊断很容易把高频段的啮合问题误判到轴承问题上去。TPA的价值就在于用同一套数据把不同频带的主导路径分开给人一个全局视角。5. 实践中的坑齿轮TPA常见问题速查表代码能跑通只是第一步实际工程应用中的麻烦通常来自数据和试验本身。以下是我在调试过程中踩过的几个典型坑整理成表格方便大家遇到问题直接对号入座。5.1 传感器耦合与相位失真的影响加速度传感器安装方式对FRF测量的相位影响很大。用磁座和用胶粘在5000Hz以上频段的共振频率会差不少。TPA计算里路径贡献量是复数叠加相位错了两路信号可能该叠加的变成抵消贡献量完全失真。我的建议是所有传感器在试验开始前统一用同一种安装方式并且做一遍“击掌测试”——用大力锤敲击测点附近的非激励位置检查各通道的相位一致性。如果两条通道的相位差随机跳动赶紧检查传感器线缆是否松动或者是否存在地回路干扰。5.2 矩阵病态条件与正则化处理逆矩阵法最怕的就是矩阵条件数过大。齿轮箱结构的路径传递函数相似度很高在共振频率附近的FRF矩阵尤其容易病态。如果不加正则化算出来的激励力会出现巨大的正负交替数值贡献量曲线完全不可解释。正则化参数的选取我是先用L曲线法定一个初值再根据重构精度做微调。实操经验是如果lambda取得太小比如小于最大特征值的0.001倍矩阵求逆结果依然会发散如果lambda取得太大大于特征值的10%结果虽然有界但是贡献量会偏离真实值。5.3 轴承通过频率与齿轮啮合频率混叠问题齿轮箱里轴承故障特征频率经常会落在齿轮啮合频率的边带附近。比如某个轴承外圈故障特征频率是97Hz跟齿轮啮合频率480Hz叠加后会在583Hz处形成一组边带。单靠频谱很难分清583Hz峰值到底是齿轮故障的调制边带还是轴承外圈故障的独立分量。这时候TPA能帮上大忙如果583Hz对应的主导路径是齿轮到轴承座的那条路径并且该路径上的FRF在583Hz处有共振峰那基本可以判定为齿轮故障调制如果主导路径切换到故障轴承另一侧的箱体测点轴承故障的概率就更大。5.4 测点布置位置的影响这个坑我在开头就提过——测点位置直接决定你能否“看得到”故障。实际布置测点时不要只放在齿轮箱顶部这种好操作的位置应该优先放在轴承座轴向和径向两个方向上。轴向测点对齿轮啮合冲击的敏感度高径向测点对轴系不平衡和轴承故障的敏感度高。两个方向都布置才能给TPA矩阵提供足够的独立信息。我还发现一个规律靠近输出端的轴承座测点往往比靠近输入端的测点更容易捕捉到齿轮故障信号。如果只能布置少量测点优先考虑故障概率最高的中间轴轴承座。6. 后续扩展从离线TPA到在线状态监测这套TPA代码跑通之后我立刻想到一个扩展方向能不能把TPA用在线监测场景。6.1 在线监控的轻量化思路完整TPA的载荷识别涉及矩阵求逆计算量虽然不大但在线系统里每次都要对每一帧数据进行矩阵运算对低功耗采集设备是个负担。我的思路是用离线数据把每条路径贡献量跟测点原始振动之间的映射关系拟合出来做一个“贡献量代理模型”。具体来说先离线跑TPA得到大量频率下的贡献量然后用线性回归或者轻量级神经网络拟合出输入各测点振动到输出各路径贡献量的映射。在线阶段只需要跑这个代理模型不需要重复做矩阵求逆。实测结果显示在路径数量不多的情况下线性回归模型就能达到90%以上的拟合精度。6.2 与深度学习故障诊断的结合接口现在很多人做齿轮箱故障诊断直接上CNN、LSTM但往往只把原始信号丢进去模型学了什么完全不可解释。TPA可以把信号在频域上按路径拆分再把每条路径的贡献量谱作为深度学习模型的特征输入而不是直接把原始频谱丢进去。这样做有两个好处一是特征本身就有明确的物理含义模型训练起来更容易收敛二是模型输出可以跟路径信息对应起来一旦判断出故障可以直接追溯到是哪条路径、哪个频带异常大大提高了诊断结果的可信度。我目前在做的一个方向是把TPA输出的各路径贡献量谱做成“彩色图谱”输入到2D CNN中相当于让模型在学习了物理拆解的基础上做进一步模式识别准确率比直接吃原始频谱高了不少而且在迁移到不同工况时稳定性更好。6.3 一套代码多种齿轮箱适配最后说一个我很满意的点。这套Matlab代码只要修改数据导入模块和路径定义就能适配不同结构的齿轮箱——平行轴、行星轮系都可以用。核心的FRF估计、载荷识别、贡献量合成三个模块完全不用动。我在试验台上验证过从一种结构的齿轮箱切换到另一种结构大概只需要半天时间做数据整理和路径重新定义算力消耗也就在一台普通笔记本上轻松跑完。对于要做多台设备状态监测的工程师来说这应该是性价比很高的方案。最后分享一个操作习惯。我每次跑TPA前都会先检查一次H矩阵的条件数是否在所有频带内保持合适范围。如果发现某些频带出现异常大的条件数我会先检查测点是否临近结构振型节点——那种位置测出来的FRF幅值极小一进矩阵就会干扰求逆。把测点微调几毫米往往就能解决问题。这种细节不会出现在教科书里但实际调试中比调正则化参数还管用。