Matlab复合故障仿真:轴承与齿轮耦合振动信号建模与诊断

Matlab复合故障仿真:轴承与齿轮耦合振动信号建模与诊断 简介本资源是一套面向机械故障诊断研究者与工程实践者的MATLAB复合故障仿真程序聚焦滚动轴承与齿轮两类典型部件同时发生故障的建模与信号生成问题适用于故障特征提取、诊断算法验证及教学实验等场景。压缩包共7个文件包含3个核心MATLAB脚本如Compound_fault_simulation_signal.m、Envelope.m等用于信号合成、包络解调与频谱分析、2张关键结果图时域波形与包络谱、1张数学公式说明图及1篇中文核心期刊CAJ文献整体大小为3.71MB。已有2262人学习下载可直接运行脚本复现复合故障振动信号获取含噪声的时频域特征数据并结合附带的包络谱图与理论公式理解故障调制机理为后续小波分析、辛几何分解等高级诊断方法提供标准化仿真基准。1. 项目缘起为什么需要复合故障仿真信号在设备状态监测与故障诊断领域滚动轴承和齿轮是旋转机械中最核心、也最易损的部件。我们经常能接触到针对单一故障比如轴承内圈剥落、齿轮断齿的仿真或诊断方法相关的Matlab代码和论文也很多。但在真实的工业现场情况往往复杂得多。一台减速机或传动箱运行多年后轴承的磨损和齿轮的点蚀、裂纹很可能会同时出现甚至相互影响、相互“掩盖”。这种“复合故障”的信号特征绝不是两个单一故障信号的简单叠加。我最初意识到这个问题的重要性是在处理一台大型风机齿轮箱的振动数据时。频谱图上既有轴承故障特征频率的边带又叠加了齿轮的啮合频率及其谐波还有大量的调制边带整个频谱图看起来一团乱麻。当时手头的诊断工具和算法无论是针对轴承的包络分析还是针对齿轮的阶次分析单独应用时效果都不理想经常出现误判或漏判。从那时起我就开始琢磨能不能在Matlab里构建一个更贴近现实的“复合故障”仿真模型这不仅能帮助我们深入理解复合故障信号的产生机理更能为后续开发更鲁棒、更智能的诊断算法提供一个可靠的“试验场”。这个“复合故障仿真信号Matlab程序”项目就是为了解决这个问题而生。它旨在模拟滚动轴承如内圈、外圈、滚动体故障和齿轮如局部损伤故障同时存在时设备振动信号会是什么样子。通过这个程序我们可以灵活地调整故障的类型、严重程度、相对相位、转速波动等参数生成海量的、标签明确的仿真数据。这对于研究复合故障信号的耦合机制、验证和对比不同诊断算法的性能、乃至训练基于深度学习的智能诊断模型都有着不可替代的价值。无论你是从事故障诊断研究的学者、工程师还是相关专业的学生掌握这样一套仿真工具都能让你对复杂故障的理解更深一层。2. 仿真模型构建从物理机理到数学公式要仿真复合故障信号不能拍脑袋随便把两个信号加在一起。我们必须回到振动的物理源头理解轴承和齿轮故障各自是如何激励起系统振动的然后再考虑它们如何在同一个系统中耦合。2.1 滚动轴承故障振动模型滚动轴承的故障点蚀、剥落可以看作一个周期性的冲击激励。当滚动体经过故障点时会产生一个短暂的冲击力。这个冲击力会激发轴承-转子系统的固有振动表现为一系列衰减振荡。其数学模型通常用一个周期性的衰减正弦波序列来近似x_bearing(t) ∑_i A_i * exp(-β(t - τ_i)) * sin(2πf_n (t - τ_i) φ) * u(t - τ_i)其中τ_i是第i个冲击发生的时间。对于恒定转速τ_i i / F_bF_b是轴承的故障特征频率内圈故障频率BPFI、外圈故障频率BPFO、滚动体故障频率BSF等。A_i是第i次冲击的幅值。在理想情况下它是常数但为了更真实我们可以引入小幅随机波动模拟载荷分布不均。β是衰减系数与系统的阻尼有关。f_n是系统被激发起的某阶固有频率通常为高频几千Hz。φ是初始相位。u(t)是单位阶跃函数保证冲击在τ_i之后才发生。关键点在于轴承故障特征频率F_b与轴的转频F_r有确定的比例关系这个比例由轴承的几何参数滚子数、接触角等决定。例如对于内圈故障BPFI (N/2) * F_r * (1 (d/D)*cosα)其中N是滚子数d是滚子直径D是节圆直径α是接触角。2.2 齿轮局部故障振动模型齿轮的局部故障如断齿、点蚀会破坏啮合过程的平稳性。在故障齿参与啮合的瞬间齿轮的啮合刚度会发生突变从而产生一个幅值调制和相位调制统称调幅-调相AM-PM效应。一个常用的简化模型是将齿轮的振动信号表示为x_gear(t) [1 A_m(t)] * cos(2πf_m t φ_m(t) φ_0)其中f_m是齿轮的啮合频率f_m Z * F_rZ是齿轮齿数F_r是齿轮所在轴的转频。A_m(t)是调幅函数。当故障齿啮合时A_m(t)会出现一个峰值因为振动加剧。A_m(t)的周期是故障齿轮的旋转周期1/F_r。φ_m(t)是调相函数。故障导致啮合瞬时错位也会引起相位波动。其变化周期同样为1/F_r。φ_0是初始相位。更精细的模型会考虑啮合刚度的时变性即使无故障由于啮合齿对数的变化刚度也是周期性变化的这构成了信号中的正常啮合成分。故障则是在此基础上叠加了一个额外的、与转频同步的调制。2.3 复合故障信号的耦合方式轴承和齿轮的故障振动最终都是通过轴承座或箱体被加速度传感器测量到。因此最简单的复合模型就是线性叠加并考虑一定的背景噪声x_total(t) x_bearing(t) x_gear(t) n(t)其中n(t)是高斯白噪声用于模拟测量噪声和其他未被建模的随机振动。然而真实的耦合可能更复杂非线性叠加故障冲击可能激发起齿轮箱体的不同模态与齿轮啮合振动产生非线性相互作用。传递路径影响轴承故障产生的冲击波传播到传感器与齿轮故障振动的传播路径不同可能会被结构滤波导致幅值和相位变化。在仿真中我们可以用不同的滤波器来粗略模拟这种效应。转速波动的影响在实际中转速并非绝对恒定。微小的转速波动会使故障特征频率和啮合频率都发生“晃动”频率调制这在频谱上会表现为谱线“涂抹”或边带增多。在仿真时我们可以让转频F_r成为一个随时间缓慢变化的函数F_r(t)然后重新计算所有依赖于F_r的频率BPFI,BPFO,f_m等。注意对于初步研究和算法验证线性叠加模型加上转速波动和噪声已经能够生成极具挑战性且足够真实的仿真信号。过于复杂的物理模型会引入大量难以确定的参数反而降低了仿真的实用性和可重复性。3. Matlab程序实现核心模块详解下面我将分模块拆解这个复合故障仿真程序的关键代码实现。我们将构建一个灵活的函数可以通过输入参数控制故障类型、严重程度、信号长度等。3.1 参数设置与基础信号生成首先我们需要定义所有的物理和仿真参数。function [t, x_composite] generate_composite_fault_signal() % 1. 基本仿真参数 fs 20000; % 采样频率 (Hz) T 2; % 信号时长 (秒) t 0:1/fs:T-1/fs; % 时间向量 N length(t); % 采样点数 % 2. 机械系统参数 % 轴转频 (假设有微小波动更真实) F_r_mean 30; % 平均转频 (Hz) F_r_fluct 0.02; % 转频波动幅度 (相对值) % 生成带有慢变波动的转频曲线 F_r F_r_mean * (1 F_r_fluct * sin(2*pi*0.5*t)); % 0.5Hz的波动 % 轴承参数 (以SKF 6205为例) N_balls 9; % 滚子数 d 7.94e-3; % 滚子直径 (m) D 39e-3; % 节圆直径 (m) alpha 0; % 接触角 (度) % 计算轴承故障特征频率 (相对于平均转频) BPFI (N_balls/2) * F_r_mean * (1 (d/D)*cosd(alpha)); % 内圈故障频率 BPFO (N_balls/2) * F_r_mean * (1 - (d/D)*cosd(alpha)); % 外圈故障频率 BSF (D/(2*d)) * F_r_mean * (1 - (d/D)^2 * cosd(alpha)^2); % 滚动体故障频率 FTF 0.5 * F_r_mean * (1 - (d/D)*cosd(alpha)); % 保持架故障频率 % 齿轮参数 Z 25; % 齿轮齿数 f_m_mean Z * F_r_mean; % 平均啮合频率 % 3. 故障参数 fault_type_bearing BPFI; % 轴承故障类型BPFI, BPFO, BSF, none fault_severity_bearing 0.8; % 轴承故障冲击强度 (0~1) fault_type_gear local; % 齿轮故障类型local(局部), none fault_severity_gear 0.6; % 齿轮故障调制深度 (0~1) % 4. 系统响应参数 f_n 3500; % 系统固有频率 (Hz) beta 800; % 衰减系数 SNR_dB 15; % 信噪比 (dB)这部分代码定义了模型的骨架。注意转频F_r被设计为随时间变化的这能模拟更真实的工况。轴承和齿轮的故障类型和严重程度被参数化便于后续灵活生成不同场景的数据。3.2 轴承故障信号生成模块接下来根据选择的故障类型生成轴承冲击序列。% 生成轴承故障振动信号 x_bearing zeros(1, N); if ~strcmp(fault_type_bearing, none) % 确定故障特征频率 switch fault_type_bearing case BPFI F_fault BPFI; case BPFO F_fault BPFO; case BSF F_fault BSF; otherwise F_fault FTF; end % 计算考虑转速波动后的冲击发生时刻 % 冲击间隔不是常数而是瞬时故障频率的倒数 % 采用累积相位的方法来定位冲击时刻 phase 2 * pi * cumsum(F_fault / fs * ones(1,N)); % 这是一个近似更好的做法是积分瞬时频率 % 更精确的方法计算瞬时故障频率 F_fault_inst(t) (几何比例) * F_r(t) F_fault_inst F_fault / F_r_mean * F_r; % 假设故障频率与转频成正比波动 inst_phase 2 * pi * cumtrapz(t, F_fault_inst); % 积分得到瞬时相位 % 找到相位穿越2π整数倍的点即为冲击发生时刻的索引近似 impact_indices find(diff(mod(inst_phase/(2*pi), 1)) 0) 1; impact_indices impact_indices(impact_indices N*0.95); % 避免末尾索引溢出 % 生成每个冲击响应 A0 fault_severity_bearing; % 基础冲击幅值 for idx 1:length(impact_indices) t_impact t(impact_indices(idx)); t_segment t(impact_indices(idx):end) - t_impact; % 只取衰减振荡的前面一部分避免计算全部效率 valid_len min(length(t_segment), round(fs / (beta/pi)) * 5); % 约5个时间常数 t_segment t_segment(1:valid_len); % 冲击幅值加入随机波动 (±20%) A_i A0 * (0.8 0.4*rand()); % 衰减正弦波 impact_wave A_i * exp(-beta * t_segment) .* sin(2*pi*f_n * t_segment); % 叠加到总信号上 signal_indices impact_indices(idx) : (impact_indices(idx)valid_len-1); if max(signal_indices) N x_bearing(signal_indices) x_bearing(signal_indices) impact_wave; end end fprintf(轴承故障(%s)模拟完成生成了%d次冲击。\n, fault_type_bearing, length(impact_indices)); end这个模块有几个关键实现细节瞬时频率积分使用cumtrapz对瞬时故障频率进行积分来得到瞬时相位这比用平均频率乘以时间更准确尤其当转速波动时。冲击时刻检测通过检测瞬时相位对2π取模的跳变点从接近1跳变到0来定位冲击发生的近似索引。这是一种在离散信号中检测周期性事件的高效方法。衰减振荡生成每次冲击都生成一个以系统固有频率f_n振荡、以β为系数衰减的正弦波。只生成有限长度约5个时间常数以提高计算效率。幅值随机性给每次冲击的幅值A_i加入随机波动模拟实际中因载荷变化导致的冲击强度不一。3.3 齿轮故障信号生成模块然后生成齿轮的调幅-调相信号。% 生成齿轮故障振动信号 x_gear zeros(1, N); if ~strcmp(fault_type_gear, none) % 计算齿轮轴的瞬时相位考虑转速波动 inst_phase_gear_shaft 2 * pi * cumtrapz(t, F_r); % 积分得到轴旋转的瞬时相位 % 1. 调幅函数 AM在故障齿啮合期间振动增大 % 假设故障导致每转一次出现一个幅值凸起 am_frequency F_r_mean; % 调幅频率等于转频 am_phase 2 * pi * cumtrapz(t, am_frequency / F_r_mean * F_r); % 同样考虑波动 % 生成一个周期性的脉冲形状作为AM函数基础 am_pulse_width 0.05; % 脉冲宽度占周期的比例 am_base 0.5 * (1 square(am_phase, am_pulse_width*100)); % 方波但值在0和1之间 % 对基础方波进行低通滤波使其更平滑更像真实的调制 [b_lp, a_lp] butter(2, 2*am_frequency/fs); % 低通滤波器截止频率略高于AM频率 am_smooth filtfilt(b_lp, a_lp, am_base); % 归一化并应用故障严重程度 am_smooth am_smooth / max(am_smooth); A_m fault_severity_gear * am_smooth; % 2. 调相函数 PM故障导致啮合瞬时相位突变 % 同样每转一次产生一个相位脉冲 pm_magnitude fault_severity_gear * 0.5; % 相位调制幅度 (弧度) % 生成一个更窄的脉冲作为PM函数 pm_pulse pm_magnitude * (0.5*(1 square(am_phase pi/2, am_pulse_width*50))); % 相位偏移的方波 [b_lp_pm, a_lp_pm] butter(2, 4*am_frequency/fs); phi_m filtfilt(b_lp_pm, a_lp_pm, pm_pulse); % 3. 生成齿轮啮合振动 % 计算瞬时啮合频率和相位 f_m_inst Z * F_r; % 瞬时啮合频率 inst_phase_mesh 2 * pi * cumtrapz(t, f_m_inst); % 瞬时啮合相位 % 生成信号 x_gear (1 A_m) .* cos(inst_phase_mesh phi_m); % 4. 可选添加无故障时的啮合刚度变化引起的振动成分 % 这会使信号更丰富但也会更复杂。此处为简化暂不添加。 fprintf(齿轮局部故障模拟完成。\n); else % 若无齿轮故障可以只生成一个纯净的啮合振动幅值恒定无调制 f_m_inst Z * F_r; inst_phase_mesh 2 * pi * cumtrapz(t, f_m_inst); x_gear 0.3 * cos(inst_phase_mesh); % 一个较小的恒定幅值啮合振动作为背景 end这个模块模拟了齿轮故障的核心——调幅(AM)和调相(PM)调幅(AM)模拟使用一个与转频同步的、平滑化的脉冲序列作为A_m(t)。square函数生成方波filtfilt进行零相位低通滤波使其平滑模拟故障齿啮合时振动逐渐增大又减小的过程。fault_severity_gear参数控制调制的深度。调相(PM)模拟原理类似用一个脉冲序列直接加到啮合相位φ_m(t)上。相位调制在频谱上会产生对称的边带是齿轮故障的另一个重要特征。瞬时频率处理和轴承部分一样啮合频率f_m_inst也是随时间变化的通过积分得到瞬时相位inst_phase_mesh这是生成非平稳信号的关键。无故障情况即使齿轮无故障通常也存在一个较弱的、恒幅的啮合振动成分。程序中用else分支来生成这个背景成分。3.4 信号合成、加噪与可视化最后将两部分信号合成并添加噪声。% 合成复合故障信号 x_composite x_bearing x_gear; % 添加高斯白噪声达到指定的信噪比(SNR) signal_power rms(x_composite)^2; if signal_power 0 noise_power signal_power / (10^(SNR_dB/10)); noise sqrt(noise_power) * randn(1, N); x_composite x_composite noise; else warning(合成信号功率为零未添加噪声。); end % 信号可视化 (时域波形) figure(Position, [100, 100, 1200, 800]); subplot(3,1,1); plot(t, x_bearing); title(轴承故障振动分量); xlabel(时间 (s)); ylabel(幅值); grid on; xlim([0, min(0.5, T)]); % 只看前0.5秒 subplot(3,1,2); plot(t, x_gear); title(齿轮故障振动分量); xlabel(时间 (s)); ylabel(幅值); grid on; xlim([0, min(0.5, T)]); subplot(3,1,3); plot(t, x_composite); title([复合故障信号 (轴承:, fault_type_bearing, 齿轮:, fault_type_gear, )]); xlabel(时间 (s)); ylabel(幅值); grid on; xlim([0, min(0.5, T)]); % 频谱分析 (以观察频率成分) figure(Position, [100, 100, 1200, 600]); NFFT 2^nextpow2(N); f fs/2 * linspace(0,1,NFFT/21); X fft(x_composite, NFFT) / N; plot(f, 2*abs(X(1:NFFT/21))); title(复合故障信号频谱); xlabel(频率 (Hz)); ylabel(幅值); grid on; xlim([0, 1000]); % 聚焦在低频段观察故障特征频率和啮合频率 % 标记关键频率 hold on; yl ylim; plot([BPFI, BPFI], yl, r--, LineWidth, 1.5, DisplayName, [BPFI, num2str(round(BPFI,1)), Hz]); plot([f_m_mean, f_m_mean], yl, g--, LineWidth, 1.5, DisplayName, [啮合频率, num2str(round(f_m_mean,1)), Hz]); legend(Location, best); hold off; end合成与后处理步骤包括线性合成将轴承冲击信号和齿轮振动信号直接相加。添加噪声根据设定的信噪比SNR_dB计算需要添加的白噪声功率然后加入高斯白噪声。rms函数用于计算信号的有效值功率。可视化生成两个图。第一个图并列显示两个故障分量和合成信号的时域波形前0.5秒便于观察冲击与调制的形态。第二个图显示合成信号的频谱并标记出轴承故障特征频率如BPFI和齿轮平均啮合频率的位置便于在频域观察故障特征。4. 仿真结果分析与故障特征解读运行上述程序例如设置轴承为内圈故障BPFI齿轮为局部故障我们可以得到仿真信号。对结果的分析是理解复合故障的关键。4.1 时域波形特征在时域波形图上前0.5秒我们可以观察到轴承分量呈现出一系列稀疏的、衰减振荡的“针状”脉冲每个脉冲之间的间隔大致相等对应BPFI的周期1/BPFI ≈ 0.0037s但由于转速波动间隔并不严格均匀。齿轮分量呈现为一个幅值被周期性调制的、频率更高的正弦波。调制周期等于轴的旋转周期1/F_r_mean ≈ 0.033s。在调制峰值处故障齿啮合幅值最大。复合信号是上述两者的叠加。轴承的冲击脉冲会“骑”在齿轮的调制波形上。在某些时刻如果冲击恰好发生在齿轮振动的波峰或波谷会导致幅值异常大如果发生在过零点则可能被掩盖。这种“掩蔽效应”是复合故障诊断难度的直观体现。4.2 频域频谱特征在频谱图上0-1000Hz我们可以看到啮合频率及其谐波在f_m_mean 750 Hz(Z25, F_r_mean30Hz) 处有一个明显的谱峰这是齿轮啮合振动的主成分。由于其幅值通常较大在频谱中占主导地位。轴承故障频率在BPFI ≈ 270 Hz附近理论上应该有一个谱峰。但在复合信号频谱中这个峰可能非常不明显甚至被淹没在噪声和啮合频率的旁瓣中。这是因为单个轴承冲击的能量分散在很宽的频带上激发了固有频率f_n3500Hz的衰减振荡其频谱是以f_n为中心的宽峰在低频段的能量很小。边带族这是复合故障尤其是齿轮故障的典型特征。在啮合频率f_m的两侧会出现以转频F_r为间隔的边带。在我们的仿真中由于同时存在调幅和调相这些边带会对称地出现在f_m ± n*F_r处。这是诊断齿轮局部故障的核心依据。潜在的调制边带轴承故障特征频率BPFI也可能作为调制频率出现在系统固有频率f_n的周围形成边带。但我们的频谱观察范围在1000Hz以下而f_n在3500Hz所以在这个低频谱上看不到。要观察这个特征需要做高频段的解调分析如包络谱。实操心得直接观察低频段频谱来诊断轴承故障在复合故障场景下往往失败。必须借助包络分析Envelope Analysis或谱峭度Spectral Kurtosis等方法先将高频共振带的轴承冲击信息解调到低频再观察解调后频谱中的故障特征频率。这是处理复合故障信号的标准流程。4.3 参数变化对信号的影响通过调整程序中的参数我们可以研究不同因素对诊断难度的影响故障严重程度增大fault_severity_bearing轴承冲击幅值变大在时域上更容易识别增大fault_severity_gear调幅调相深度增加边带更明显。信噪比(SNR)降低SNR_dB噪声会掩盖故障特征尤其是能量较弱的轴承故障特征。转速波动增大F_r_fluct频谱上的谱线会变宽“涂抹”效应边带结构变得模糊增加频率识别的难度。故障频率接近如果轴承故障频率如BPFO与齿轮转频F_r或其倍频接近它们的边带可能会交织在一起难以区分。5. 高级扩展与算法验证应用基础的仿真程序已经搭建完成但它的价值在于可扩展性和作为算法“试金石”的用途。5.1 模型扩展更复杂的故障与工况多故障点仿真可以同时模拟轴承内圈和外圈故障或者齿轮多个齿存在故障。只需在对应模块中叠加多个冲击序列或多个调制函数即可。时变故障故障的严重程度可以是时变的。例如让fault_severity_bearing随时间线性增加模拟故障从萌生到发展的过程。非高斯噪声与冲击噪声工业现场噪声不一定是高斯白噪声。可以添加周期性干扰如电机工频、随机冲击噪声模拟其他偶然碰撞来增加仿真难度。传递路径建模用不同的FIR或IIR滤波器分别处理x_bearing(t)和x_gear(t)模拟振动从源到传感器的不同传递特性再合成。5.2 作为诊断算法的测试基准生成的仿真数据可以保存下来用于系统性地测试和比较各种诊断算法传统方法测试包络谱分析对原始信号进行带通滤波围绕系统固有频率f_n然后希尔伯特变换求包络再对包络做FFT。成功的算法应该在包络谱中清晰地显示出轴承故障频率BPFI及其谐波。阶次分析如果仿真中包含了转速波动阶次分析比频谱分析更有效。可以验证算法能否在阶次谱中稳定地识别出齿轮故障对应的阶次Z阶及其边带±1阶。倒频谱分析对于齿轮故障倒频谱能突出调制的周期性即转频1/F_r有助于在复杂频谱中提取故障信息。现代智能诊断算法验证特征提取与分类生成大量不同故障类型、不同严重程度、不同噪声水平的仿真样本。从中提取时域、频域、时频域特征如均方根、峭度、小波包能量、谱峰因子等构建特征数据集。用这个数据集来训练和测试SVM、随机森林等机器学习模型评估其分类精度和泛化能力。深度学习端到端诊断直接将仿真信号的一维时域数据或二维时频图如连续小波变换尺度图作为输入训练CNN、LSTM或Transformer模型。仿真数据的优势在于标签绝对准确、样本量可无限生成、可控制变量是研究网络结构、优化超参数的理想工具。算法鲁棒性评估可以固定一种故障组合逐步增加噪声水平降低SNR观察某种诊断算法如某个特征指标的性能曲线如检测概率 vs. SNR从而定量评估其抗噪能力。可以引入转速波动、载荷变化等非平稳因素测试算法对工况变化的适应性。这个复合故障仿真程序就像一个功能完备的“数字孪生”试验台。它剥离了真实数据中诸多不可控的干扰让我们能够在一个受控的环境下深入探究故障的机理、验证方法的有效性、并暴露出算法的弱点。当你为一个新算法在此仿真数据上取得好效果而欣喜时别忘了这仅仅是迈向真实复杂工业应用的第一步。但这一步因其清晰和可控至关重要。本文还有配套的精品资源点击获取