使用MATLAB对滚动轴承故障信号进行形态学滤波
形态学滤波非常适合于滚动轴承故障信号处理，因为它能有效地从强烈的背景噪声中提取出由局部损伤（点蚀、裂纹）引起的瞬态冲击成分，同时很好地保留冲击的形态特征。
一、 基本概念：数学形态学滤波
数学形态学滤波源于图像处理，但其理念可以很好地应用于一维信号。其核心是通过一个称为结构元素的探针在信号上“移动”，通过特定的运算来提取信号中的特定形状。
两个基本操作：
-
腐蚀：
- 作用： 收缩信号，平滑峰值，消除小的孤立的脉冲噪声。
- 效果：
输出信号 = 输入信号与结构元素重叠区域的最小值。它会“削平”尖峰。
-
膨胀：
- 作用： 扩张信号，填充谷底，连接短的间断。
- 效果：
输出信号 = 输入信号与结构元素重叠区域的最大值。它会“填平”低谷。
由基本操作构成的滤波器：
-
开运算： 先腐蚀后膨胀。
- 效果： 消除信号中尖锐的峰值噪声，同时整体形态不变。非常适合去除正脉冲噪声。
-
闭运算： 先膨胀后腐蚀。
- 效果： 消除信号中尖锐的谷底噪声，同时整体形态不变。非常适合去除负脉冲噪声。
-
形态学梯度： 膨胀结果减去腐蚀结果。
- 效果： 突出信号的边缘（即变化剧烈的地方），非常适合增强故障冲击的上升沿和下降沿。
-
高帽/顶帽变换： 原信号减去开运算结果。
- 效果： 提取出被开运算消除掉的亮细节（即峰值部分）。这是轴承故障诊断中最常用的工具，因为它能完美地将周期性冲击成分从背景信号中“剥离”出来。
二、 步骤与代码
我们将使用一个模拟的轴承故障信号来演示整个过程。您也可以将自己的数据替换进去。
步骤 1：生成或加载轴承故障信号
首先，我们模拟一个包含噪声的轴承故障冲击信号。
% 参数设置
fs = 12000; % 采样频率 12kHz
t = 0:1/fs:1; % 1秒时间轴
N = length(t);
% 1. 模拟轴承故障冲击（周期性脉冲）
f_bpf = 100; % 轴承故障特征频率 (BPFO, BPFI, BSF等，这里假设100Hz)
T = 1/f_bpf; % 冲击周期
duty_ratio = 0.02; % 脉冲占空比
impact_train = pulstran(t, 0:T:1, 'rectpuls', duty_ratio*T);
% 2. 对冲击进行衰减，模拟真实情况（一个衰减振荡）
damp_factor = 800;
osc_freq = 3000; % 轴承系统的共振频率
osc_wave = exp(-damp_factor * t) .* sin(2 * pi * osc_freq * t);
osc_wave = osc_wave / max(abs(osc_wave)); % 归一化
% 将冲击序列与振荡波形卷积，得到真实的故障信号模型
fault_signal = conv(impact_train, osc_wave, 'same');
% 3. 添加噪声
noise_intensity = 0.6;
noisy_signal = fault_signal + noise_intensity * randn(1, N);
% 4. 可视化原始信号与含噪信号
figure;
subplot(3, 1, 1);
plot(t, fault_signal);
title('模拟的轴承故障冲击信号（无噪声）');
xlabel('时间 (s)');
ylabel('幅值');
grid on;
subplot(3, 1, 2);
plot(t, noisy_signal);
title('添加噪声后的信号');
xlabel('时间 (s)');
ylabel('幅值');
grid on;
步骤 2：设计结构元素
结构元素的选择是形态学滤波成败的关键。对于轴承故障信号，通常选择扁平结构元素，其形状为一段恒定幅值的线段。
- 长度： 这是最重要的参数。它应该与您要提取的冲击的宽度相匹配。
- 太短： 无法有效平滑噪声。
- 太长： 会平滑甚至滤除真实的故障冲击。
- 经验值： 通常为共振周期（1/共振频率）的 0.6 ~ 0.9 倍。可以通过信号的频谱分析来估算共振频率。
% 设计结构元素
resonance_period = 1 / osc_freq; % 共振周期
se_length = round(0.7 * resonance_period * fs); % 计算结构元素长度（点数）
se = strel('line', se_length, 0); % 创建一个‘扁平’的线性结构元素，角度为0度
disp(['结构元素长度: ', num2str(se_length), ' 个点']);
步骤 3：应用形态学滤波
这里我们演示最常用的开运算和高帽变换。
% 1. 开运算滤波 (去除正脉冲噪声)
open_signal = imopen(noisy_signal, se.Neighborhood);
% 2. 高帽变换 (提取故障冲击)
tophat_signal = imtophat(noisy_signal, se.Neighborhood);
% 注意：MATLAB的imopen/imtophat默认用于图像（2D），但也可以处理1D向量。
% 另一种更严谨的方法是直接使用腐蚀和膨胀函数：
% open_signal = imdilate(imerode(noisy_signal, se.Neighborhood), se.Neighborhood);
% tophat_signal = noisy_signal - open_signal;
步骤 4：结果可视化与分析
% 可视化滤波结果
subplot(3, 1, 3);
plot(t, tophat_signal);
title('形态学滤波后（高帽变换结果）');
xlabel('时间 (s)');
ylabel('幅值');
grid on;
% 为了更清晰地对比，单独画一个图
figure;
ax1 = subplot(2, 1, 1);
plot(t, noisy_signal, 'b');
title('原始含噪信号');
ylabel('幅值');
grid on;
ax2 = subplot(2, 1, 2);
plot(t, tophat_signal, 'r');
title('高帽变换提取出的冲击成分');
xlabel('时间 (s)');
ylabel('幅值');
grid on;
linkaxes([ax1, ax2], 'x'); % 链接两个子图的x轴，便于缩放对比
% 频谱分析对比
figure;
f = (0:N/2-1)*fs/N;
Y_original = abs(fft(noisy_signal));
Y_original = Y_original(1:N/2);
Y_filtered = abs(fft(tophat_signal));
Y_filtered = Y_filtered(1:N/2);
plot(f, 20*log10(Y_original/max(Y_original)), 'b');
hold on;
plot(f, 20*log10(Y_filtered/max(Y_filtered)), 'r', 'LineWidth', 1.5);
title('频谱对比');
xlabel('频率 (Hz)');
ylabel('归一化幅值 (dB)');
legend('原始信号', '形态学滤波后');
xlim([0, 1500]); % 聚焦在低频段观察故障特征频率及其倍频
grid on;
参考代码 MATLAB 滚动轴承故障信号进行形态学滤波 www.youwenfan.com/contentcnh/54986.html
三、 关键要点与技巧
-
结构元素选择：
- 形状： 对于振动信号，
‘flat’（扁平）是最常用和最有效的。 - 长度： 这是核心参数。需要通过先验知识（如已知的共振频率）或频谱分析来反复调试。可以尝试不同的长度，观察滤波效果，选择能最清晰提取出周期性冲击的那个长度。
- 形状： 对于振动信号，
-
组合滤波：
- 有时单一操作效果不佳，可以组合使用。例如，先进行闭-开运算，再进行高帽变换，可以同时抑制正负脉冲噪声。
-
与包络分析结合：
- 形态学滤波（尤其是高帽变换）的输出是增强后的冲击信号，但这个信号可能仍然包含一些残余噪声。可以对其进一步进行希尔伯特包络解调，然后对包络谱进行分析，从而非常清晰地看到故障特征频率及其倍频，这是诊断轴承故障的“金标准”。
% 对形态学滤波结果进行包络分析 analytic_signal = hilbert(tophat_signal); envelope_signal = abs(analytic_signal); % 计算包络谱 Y_envelope = abs(fft(envelope_signal)); Y_envelope = Y_envelope(1:N/2); f_env = (0:N/2-1)*fs/N; figure; plot(f_env, 20*log10(Y_envelope/max(Y_envelope))); title('滤波后信号的包络谱'); xlabel('频率 (Hz)'); ylabel('幅值 (dB)'); xlim([0, 500]); % 聚焦在故障特征频率附近 grid on; % 在包络谱图中，你应该能在f_bpf (100Hz) 及其倍频处看到明显的谱线。
浙公网安备 33010602011771号