轴承故障诊断和趋势预测是工业设备健康管理的核心内容，频域特征提取在这方面发挥着至关重要的作用。
1. 频域分析的基本原理
轴承振动信号的频域分析基于傅里叶变换，将时域信号转换为频域表示，从而揭示信号的频率组成特征。轴承故障会产生特定的频率成分，这些成分可以作为故障诊断的依据。
1.1 轴承故障特征频率计算
不同类型的轴承故障会产生不同的特征频率：
- 内圈故障频率(BPFI):
BPFI = (n/2) * f_r * (1 + (d/D) * cosφ) - 外圈故障频率(BPFO):
BPFO = (n/2) * f_r * (1 - (d/D) * cosφ) - 滚动体故障频率(BSF):
BSF = (D/(2*d)) * f_r * (1 - (d/D)^2 * cos^2φ) - 保持架故障频率(FTF):
FTF = (1/2) * f_r * (1 - (d/D) * cosφ)
其中：
n= 滚动体数量f_r= 轴旋转频率 (Hz)d= 滚动体直径D= 轴承节径φ= 接触角
2. MATLAB实现频域特征提取
2.1 基本频域特征提取函数
functionfeatures=extractFrequencyFeatures(signal,fs)% 提取信号的频域特征% 输入:% signal - 输入信号% fs - 采样频率% 输出:% features - 包含各种频域特征的结构体% 计算FFTN=length(signal);fft_data=fft(signal);fft_magnitude=abs(fft_data(1:floor(N/2)+1);fft_magnitude=fft_magnitude/max(fft_magnitude);% 归一化freq_axis=(0:floor(N/2))*fs/N;% 计算功率谱密度 (PSD)[psd,freq_psd]=pwelch(signal,hanning(256),128,256,fs);% 初始化特征结构体features=struct();% 1. 重心频率 (Spectral Centroid)features.centroid=sum(freq_axis.*fft_magnitude)/sum(fft_magnitude);% 2. 均方频率 (Mean Square Frequency)features.mean_square_freq=sum(freq_axis.^2.*fft_magnitude)/sum(fft_magnitude);% 3. 频率方差 (Frequency Variance)features.freq_variance=sum((freq_axis-features.centroid).^2.*fft_magnitude)/sum(fft_magnitude);% 4. 频率标准差 (Frequency Standard Deviation)features.freq_std=sqrt(features.freq_variance);% 5. 频谱偏度 (Spectral Skewness)features.skewness=sum((freq_axis-features.centroid).^3.*fft_magnitude)/...(sum(fft_magnitude)*features.freq_std^3);% 6. 频谱峭度 (Spectral Kurtosis)features.kurtosis=sum((freq_axis-features.centroid).^4.*fft_magnitude)/...(sum(fft_magnitude)*features.freq_variance^2);% 7. 频带能量特征bands=[0,0.1*fs/2;0.1*fs/2,0.2*fs/2;0.2*fs/2,0.3*fs/2;0.3*fs/2,0.4*fs/2;0.4*fs/2,0.5*fs/2];fori=1:size(bands,1)band_mask=freq_axis>=bands(i,1)&freq_axis<=bands(i,2);features.(sprintf('band_energy_%d',i))=sum(fft_magnitude(band_mask));end% 8. 峰值频率及其幅度[pks,locs]=findpeaks(fft_magnitude,'SortStr','descend','NPeaks',5);iflength(pks)>=5fori=1:5features.(sprintf('peak_freq_%d',i))=freq_axis(locs(i));features.(sprintf('peak_amp_%d',i))=pks(i);endend% 9. 频谱熵 (Spectral Entropy)pdf=fft_magnitude/sum(fft_magnitude);% 概率密度函数features.spectral_entropy=-sum(pdf.*log2(pdf+eps));% 10. 包络谱特征 (Envelope Spectrum Features)[env_features,env_freq]=extractEnvelopeFeatures(signal,fs);features.envelope_peak_ratio=env_features.peak_ratio;features.envelope_centroid=env_features.centroid;endfunction[env_features,freq_axis]=extractEnvelopeFeatures(signal,fs)% 提取包络谱特征% 希尔伯特变换获取包络analytic_signal=hilbert(signal);envelope=abs(analytic_signal);% 计算包络谱N=length(envelope);env_spectrum=abs(fft(envelope-mean(envelope)));env_spectrum=env_spectrum(1:floor(N/2)+1);freq_axis=(0:floor(N/2))*fs/N;% 归一化env_spectrum=env_spectrum/max(env_spectrum);% 提取包络谱特征env_features.centroid=sum(freq_axis.*env_spectrum)/sum(env_spectrum);% 计算峰值比例 (故障特征)[pks,locs]=findpeaks(env_spectrum,'MinPeakHeight',0.1);iflength(pks)>1sorted_pks=sort(pks,'descend');env_features.peak_ratio=sorted_pks(2)/sorted_pks(1);elseenv_features.peak_ratio=0;endend2.2 轴承故障诊断系统
functionbearingFaultDiagnosis()% 轴承故障诊断主函数% 加载数据（这里使用模拟数据，实际应用中应替换为真实数据）[healthy_data,faulty_data,fs]=loadBearingData();% 提取健康轴承特征healthy_features=[];fori=1:size(healthy_data,2)features=extractFrequencyFeatures(healthy_data(:,i),fs);healthy_features=[healthy_features;struct2array(features)];end% 提取故障轴承特征faulty_features=[];fori=1:size(faulty_data,2)features=extractFrequencyFeatures(faulty_data(:,i),fs);faulty_features=[faulty_features;struct2array(features)];end% 创建标签labels=[zeros(size(healthy_features,1),1);ones(size(faulty_features,1),1)];all_features=[healthy_features;faulty_features];% 特征标准化all_features=zscore(all_features);% 训练分类器（使用SVM）rng(1);% 设置随机种子以确保可重复性cv=cvpartition(labels,'HoldOut',0.3);idx_train=training(cv);idx_test=test(cv);% 使用PCA进行特征降维[coeff,score,~,~,explained]=pca(all_features(idx_train,:));num_components=find(cumsum(explained)>=95,1);% 保留95%方差train_features=score(:,1:num_components);test_features=(all_features(idx_test,:)-mean(all_features(idx_train,:)))*coeff(:,1:num_components);% 训练SVM分类器svm_model=fitcsvm(train_features,labels(idx_train),...'KernelFunction','rbf','Standardize',true,...'KernelScale','auto','BoxConstraint',1);% 预测predictions=predict(svm_model,test_features);% 评估性能accuracy=sum(predictions==labels(idx_test))/length(labels(idx_test));confusion_mat=confusionmat(labels(idx_test),predictions);fprintf('诊断准确率: %.2f%%\n',accuracy*100);disp('混淆矩阵:');disp(confusion_mat);% 可视化特征空间visualizeFeatures(train_features,labels(idx_train),svm_model);endfunction[healthy_data,faulty_data,fs]=loadBearingData()% 加载轴承数据（这里使用模拟数据）% 在实际应用中，应替换为真实数据加载代码fs=12000;% 采样频率 12kHzt=0:1/fs:1-1/fs;% 生成健康轴承数据（主要是随机噪声）healthy_data=zeros(length(t),10);fori=1:10healthy_data(:,i)=0.1*randn(size(t))+0.05*sin(2*pi*50*t);end% 生成故障轴承数据（包含故障特征频率）faulty_data=zeros(length(t),10);fault_freq=100;% 假设故障特征频率为100Hzfori=1:10% 基础噪声base_signal=0.1*randn(size(t))+0.05*sin(2*pi*50*t);% 添加故障冲击（周期性冲击）impulse_train=zeros(size(t));impulse_period=round(fs/fault_freq);impulse_train(1:impulse_period:end)=1;% 冲击响应（衰减正弦波）impulse_response=exp(-500*t).*sin(2*pi*2000*t);fault_signal=conv(impulse_train,impulse_response,'same');faulty_data(:,i)=base_signal+0.3*fault_signal/max(abs(fault_signal));endendfunctionvisualizeFeatures(features,labels,svm_model)% 可视化特征空间和分类边界figure;% 如果特征维度大于2，使用前两个主成分ifsize(features,2)>2[~,score]=pca(features);features_2d=score(:,1:2);elsefeatures_2d=features;end% 创建网格用于绘制决策边界h=0.02;x1_min=min(features_2d(:,1))-1;x1_max=max(features_2d(:,1))+1;x2_min=min(features_2d(:,2))-1;x2_max=max(features_2d(:,2))+1;[xx1,xx2]=meshgrid(x1_min:h:x1_max,x2_min:h:x2_max);% 预测网格点的类别Z=predict(svm_model,[xx1(:),xx2(:)]);Z=reshape(Z,size(xx1));% 绘制决策区域contourf(xx1,xx2,Z,'AlphaData',0.3);colormap([0.80.90.9;0.90.80.9]);hold on;% 绘制数据点gscatter(features_2d(:,1),features_2d(:,2),labels,'br','xo',10);xlabel('主成分 1');ylabel('主成分 2');title('轴承状态分类结果');legend('健康','故障','Location','best');grid on;hold off;end2.3 趋势预测与健康指标构建
functionhealthIndicatorTrendPrediction()% 轴承健康指标构建与趋势预测% 生成模拟时间序列数据（实际应用中应使用真实数据）[time_series,time]=generateTimeSeriesData();% 提取每个时间点的频域特征features_over_time=[];fori=1:size(time_series,2)feature_vector=extractFrequencyFeatures(time_series(:,i),12000);features_over_time=[features_over_time;struct2array(feature_vector)];end% 构建健康指标（使用多个特征的加权组合）% 这里使用第一主成分作为健康指标[coeff,score,~,~,explained]=pca(features_over_time);health_indicator=score(:,1);% 健康指标归一化（0-1范围，1表示健康，0表示故障）health_indicator=1-(health_indicator-min(health_indicator))/...(max(health_indicator)-min(health_indicator));% 趋势预测（使用ARIMA模型）num_train=floor(0.7*length(health_indicator));train_data=health_indicator(1:num_train);test_data=health_indicator(num_train+1:end);% 训练ARIMA模型arima_model=arima(5,1,2);% ARIMA(5,1,2)estimated_model=estimate(arima_model,train_data);% 预测[yF,yMSE]=forecast(estimated_model,length(test_data),'Y0',train_data);lower=yF-1.96*sqrt(yMSE);upper=yF+1.96*sqrt(yMSE);% 计算预测误差mse=mean((test_data-yF).^2);rmse=sqrt(mse);mae=mean(abs(test_data-yF));fprintf('趋势预测性能:\n');fprintf('RMSE: %.4f\n',rmse);fprintf('MAE: %.4f\n',mae);% 可视化结果plotTrendResults(time,health_indicator,train_data,test_data,yF,lower,upper);% 预测剩余使用寿命 (RUL)predictRUL(health_indicator,time);endfunction[time_series,time]=generateTimeSeriesData()% 生成模拟时间序列数据（轴承退化过程）fs=12000;measurement_interval=3600;% 每小时的测量间隔（秒）total_time=30*24*3600;% 30天的总时间（秒）time=0:measurement_interval:total_time;time_series=zeros(floor(fs*0.5),length(time));% 每次测量0.5秒数据% 生成退化过程数据fori=1:length(time)% 基础信号（随着时间衰减）t_signal=0:1/fs:0.5-1/fs;base_signal=0.1*randn(size(t_signal))+0.05*sin(2*pi*50*t_signal);% 故障发展（随时间增强）degradation_level=1-exp(-0.0001*time(i));% 退化水平fault_freq=100+50*degradation_level;% 故障频率随时间变化% 故障冲击impulse_train=zeros(size(t_signal));impulse_period=round(fs/fault_freq);impulse_train(1:impulse_period:end)=1;impulse_response=exp(-500*t_signal).*sin(2*pi*2000*t_signal);fault_signal=conv(impulse_train,impulse_response,'same');% 组合信号time_series(:,i)=base_signal+degradation_level*0.5*fault_signal/max(abs(fault_signal));endendfunctionplotTrendResults(time,health_indicator,train_data,test_data,yF,lower,upper)% 绘制趋势预测结果figure;% 转换时间为天time_days=time/(24*3600);num_train=length(train_data);% 绘制实际健康指标plot(time_days,health_indicator,'b-','LineWidth',1.5,'DisplayName','实际健康指标');hold on;% 绘制预测结果plot(time_days(num_train+1:end),yF,'r--','LineWidth',1.5,'DisplayName','预测值');fill([time_days(num_train+1:end),fliplr(time_days(num_train+1:end))],...[lower',fliplr(upper')],'r','FaceAlpha',0.2,'EdgeColor','none',...'DisplayName','95%置信区间');% 添加分界线plot([time_days(num_train),time_days(num_train)],[0,1],'k--','LineWidth',1,'DisplayName','训练/测试分界');xlabel('时间 (天)');ylabel('健康指标');title('轴承健康趋势预测');legend('show');grid on;hold off;endfunctionpredictRUL(health_indicator,time)% 预测剩余使用寿命 (RUL)% 设置故障阈值failure_threshold=0.2;% 找到健康指标首次低于阈值的时间点failure_idx=find(health_indicator<failure_threshold,1);ifisempty(failure_idx)fprintf('在当前监测期内未检测到故障发展至阈值\n');return;end% 使用线性回归预测RULlast_n_points=20;% 使用最后20个点进行预测iflength(health_indicator)<last_n_points last_n_points=length(health_indicator);end% 选择用于拟合的数据点fit_indices=(length(health_indicator)-last_n_points+1):length(health_indicator);x=time(fit_indices)';y=health_indicator(fit_indices)';% 线性拟合p=polyfit(x,y,1);y_fit=polyval(p,x);% 求解达到故障阈值的时间% p(1)*t + p(2) = failure_thresholdfailure_time=(failure_threshold-p(2))/p(1);% 计算RULcurrent_time=time(end);rul=failure_time-current_time;fprintf('剩余使用寿命预测:\n');fprintf('当前时间: %.2f 天\n',current_time/(24*3600));fprintf('预测故障时间: %.2f 天\n',failure_time/(24*3600));fprintf('剩余使用寿命: %.2f 天\n',rul/(24*3600));% 绘制RUL预测图figure;plot(time/(24*3600),health_indicator,'b-','LineWidth',1.5,'DisplayName','健康指标');hold on;plot(x/(24*3600),y_fit,'r--','LineWidth',1.5,'DisplayName','趋势拟合');plot([min(x),failure_time]/(24*3600),[failure_threshold,failure_threshold],'k:','LineWidth',1,'DisplayName','故障阈值');plot([failure_time,failure_time]/(24*3600),[0,failure_threshold],'k:','LineWidth',1);xlabel('时间 (天)');ylabel('健康指标');title(sprintf('剩余使用寿命预测: %.2f 天',rul/(24*3600)));legend('show');grid on;hold off;end参考代码 一维信号频域特征提取可用于轴承故障诊断和趋势预测www.3dddown.com/csa/50931.html
3. 实际应用建议
3.1 数据预处理注意事项
降噪处理：在实际工业环境中，振动信号通常包含大量噪声，建议使用小波降噪或自适应滤波等方法进行预处理。
转速变化处理：如果设备转速变化较大，应考虑使用阶比分析(order analysis)代替传统的FFT分析。
数据标准化：不同传感器、不同时间采集的数据可能存在尺度差异，应进行适当的标准化处理。
3.2 特征选择与优化
相关性分析：计算各个特征与故障状态的相关性，选择相关性高的特征。
递归特征消除：使用递归特征消除(RFE)方法选择最优特征子集。
领域知识引导：结合轴承故障机理知识，优先选择与故障物理过程直接相关的特征。
3.3 系统部署考虑
实时性要求：根据实际应用的实时性要求，优化算法复杂度。
自适应阈值：设置自适应故障阈值，避免因工况变化导致的误报。
模型更新机制：建立模型定期更新机制，适应设备老化等长期变化。
4. 总结
频域特征提取在轴承故障诊断和趋势预测中具有重要作用。通过提取频谱重心、频带能量、包络谱特征等多种频域特征，并结合机器学习方法，可以构建高效的故障诊断和预测系统。
