企业官网建设流程全解析

搞设备状态监测和故障诊断这些年，齿轮箱这个“硬骨头”我啃的次数最多。齿轮箱一响，黄金万两的代价不是谁都付得起，但真正麻烦的不是坏了，而是坏了之后你找不准它在哪里坏、振动是怎么传出来的。大多数现场工程师会遇到一个很实际的困惑：传感器就装在轴承座上，传感器测到的振动信号里，混着齿轮啮合激励、轴承故障冲击、轴系不平衡、箱体共振等一大堆东西，你想定位的故障源可能离测点还有两三个传递环节。这时候如果只盯着传感器信号做频谱分析，往往看到的是多条传递路径叠加后的“混合体”，很难直接判断故障到底出在内圈还是外圈、是齿轮齿面磨损还是轴弯。这个项目要解决的，就是用传递路径分析（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条路径的频率响应函数（FRF），Fi是第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 Hz，4 s的数据基本够用。如果数据长度不够，别硬加窗长，改用ZFFT（Zoom-FFT）对啮合频率附近做频谱细化，也能在不增加总时长的前提下把边频带分开。

4.4 噪声干扰严重，低频贡献被掩盖

现场传感器噪声、工频干扰、其他设备振动都会污染TPA结果。响应用的是功率谱或互谱，不是直接对时域做除法——先用H1估计或者H2估计把响应谱做平滑，再做载荷识别，抗噪能力会好很多。此外，如果关心的是啮合频率附近的贡献量，可以用窄带带通滤波器先把信号滤一遍，再做TPA，能有效压低低频大幅值分量对矩阵求逆的“支配效应”。

4.5 常见问题速查表

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给了你“路径贡献”这个维度，但别把它当成唯一的金标准——和频谱、时频、包络谱配合使用，齿轮箱的故障才能定位得更准、更快。