HO-VMD-TCN西储大学轴承故障诊断
HO-VMD-TCN西储大学轴承故障诊断
方法概述
HO-VMD-TCN方法是一种结合高阶变分模态分解（High-Order Variational Mode Decomposition, HO-VMD）与时序卷积网络（Temporal Convolutional Network, TCN）的轴承故障诊断方法。该方法特别适用于处理非线性和非平稳振动信号，通过频率分量的精确提取与时序模式的深度学习分析，提升故障诊断的准确性。
核心原理
-
HO-VMD分解信号：
- HO-VMD是一种改进的信号分解技术，可将复杂振动信号分解为若干本征模态函数（IMFs），以实现频谱的精准分离。
- 高阶变分正则化进一步抑制分解过程中的模态混叠，提高对噪声信号的鲁棒性。
-
TCN提取时序特征：
- TCN是一种基于卷积的时序建模方法，使用因果卷积和扩张卷积（dilated convolution）处理序列数据。
- 它具备高效的长时间依赖建模能力，能够从分解后的模态信号中提取关键的时序模式，完成故障特征的深度学习。
诊断流程
-
信号采集：
- 从西储大学轴承数据集中提取振动信号，包含不同运行工况下的正常与故障数据。
-
预处理：
- 对信号进行降噪与归一化处理，以提高后续分析的精度。
-
信号分解：
- 使用HO-VMD将原始振动信号分解为若干本征模态，提取各频率分量的关键特征。
-
特征提取与分类：
- 将分解后的IMFs输入TCN网络，通过卷积操作学习时序特征。
- 使用全连接层与分类器进行健康状态与故障类型的识别。
-
结果输出：
- 输出诊断准确率、混淆矩阵、以及不同故障类型的分类结果。
方法优势
- 高效信号分解：HO-VMD能够准确分离频率分量，适应非平稳信号的复杂特性。
- 鲁棒时序建模：TCN在处理长序列数据时表现出色，避免了循环神经网络（RNN）中的梯度消失问题。
- 适应复杂工况：该方法对轴承信号中的噪声和多种运行工况具有良好的适应性。
应用场景
HO-VMD-TCN可用于工业设备的实时健康监测与故障诊断，特别是在轴承、齿轮等旋转机械系统中，帮助预测潜在故障，减少意外停机损失。
%% 以最小包络熵、最小样本熵、最小信息熵、最小排列熵，排列熵/互信息熵，为目标函数（任选其一），采用HO算法优化VMD，求取VMD最佳的两个参数
clear
clc
close all
addpath(genpath(pwd))
%% 选择适应度函数
xz = 1;
if xz == 1
fobj=@EnvelopeEntropyCost; % 最小包络熵
elseif xz == 2
fobj=@SampleEntropyCost; % 最小样本熵
elseif xz == 3
fobj=@infoEntropyCost; % 最小信息熵
elseif xz == 4
fobj=@PermutationEntropyCost; % 最小排列熵
elseif xz == 5
fobj=@compositeEntropyCost; % 复合指标：排列熵/互信息熵
end
%% 读取数据（一共三种转速数据，已整理完，实际使用时选择一种即可）
res = xlsread('西储大学驱动端振动数据(转速1797).xlsx');
%% 设置参数
D = 2; % 优化变量数目
lb = [100 3]; % 下限值，分别是a,k
ub = [2500 10]; % 上限值
T = 8; % 最大迭代数目
N = 6; % 种群规模
vmddata = [];
samplenum = size(res,1)/10;
for i = 1:10 % 十种故障状态
disp(['正在对第',num2str(i),'个故障类型的数据进行VMD优化！'])
% 一种状态是120/10个样本，每次选120/10个样本进行VMD优化和特征提取
every_data = res(1+samplenum*(i-1):samplenum*i,:);
% 从当前状态的数据中任选一组数据进行VMD优化
da = every_data(1,:);
% 调用HO算法
[Best_score,Best_pos,Bestidx,curve] = HO(N,T,lb,ub,D,fobj,da');
% 输出最佳位置
display(['第',num2str(i),'个故障类型的最佳VMD参数是：', num2str(fix(Best_pos)),'，最佳IMF分量是：IMF',num2str(Bestidx)]);
% 最佳位置取整
bbh = fix(Best_pos);
% 进行9种时域指标特征提取，将优化得到的两个参数和最小适应度的索引值带回VMD中，提取得到当前状态的特征向量
new_data = featureExtraction(bbh(1),bbh(2),Bestidx,every_data);
% 将每个状态提取得到的特征向量都放在一起
vmddata = [vmddata;new_data];
end
data = vmddata; % data变量中为处理好的数据
bv = 120; % 每种状态数据有120组
%% 加标签值
hhh = size(data,2);
for i = 1:size(data,1)/bv
data(1+bv*(i-1):bv*i,hhh+1) = i;
end
%% 将处理好的处理输出成Excel
xlswrite('HO-VMD特征提取后的数据集.xlsx', data);
%% 分析数据
num_class = length(unique(data(:, end))); % 类别数（Excel最后一列放类别）
num_dim = size(data, 2) - 1; % 特征维度
num_res = size(data, 1); % 样本数（每一行，是一个样本）
num_size = 0.7; % 训练集占数据集的比例
data = data(randperm(num_res), :); % 打乱数据集（不打乱数据时，注释该行）
flag_conusion = 1; % 标志位为1，打开混淆矩阵（要求2018版本及以上）
%% 设置变量存储数据
P_train = []; P_test = [];
T_train = []; T_test = [];
%% 划分数据集
for i = 1 : num_class
mid_res = data((data(:, end) == i), :); % 循环取出不同类别的样本
mid_size = size(mid_res, 1); % 得到不同类别样本个数
mid_tiran = round(num_size * mid_size); % 得到该类别的训练样本个数
P_train = [P_train; mid_res(1: mid_tiran, 1: end - 1)]; % 训练集输入
T_train = [T_train; mid_res(1: mid_tiran, end)]; % 训练集输出
P_test = [P_test; mid_res(mid_tiran + 1: end, 1: end - 1)]; % 测试集输入
T_test = [T_test; mid_res(mid_tiran + 1: end, end)]; % 测试集输出
end
%% 数据转置
P_train = P_train'; P_test = P_test';
T_train = T_train'; T_test = T_test';
%% 得到训练集和测试样本个数
M = size(P_train, 2);
N = size(P_test , 2);
%% 数据归一化
[p_train, ps_input] = mapminmax(P_train, 0, 1);
p_test = mapminmax('apply', P_test, ps_input);
t_train = categorical(T_train);
t_test = categorical(T_test );
%% 数据格式转换
for i = 1: M
pc_train{1, i} = p_train(:, i);
tc_train{1, i} = t_train(:, i);
end
for i = 1: N
pc_test {1, i} = p_test(:, i);
tc_test {1, i} = t_test(:, i);
end
%% 设置网络参数
numFilters = 16; % 卷积核个数
filterSize = 3; % 卷积核大小
dropoutFactor = 0.05; % 空间丢失因子
numBlocks = 2; % 残差块个数
numFeatures = num_dim; % 特征个数
numClasses = num_class; % 类别个数
%% 输入层结构
layer = sequenceInputLayer(numFeatures, Normalization = "rescale-symmetric", Name = "input");
%% 将输入层加入空白网络
lgraph = layerGraph(layer);
outputName = layer.Name;
%% 建立网络结构 -- 残差块
for i = 1 : numBlocks
dilationFactor = 2 ^(i - 1); % 膨胀因子
% 建立网络结构
layers = [
convolution1dLayer(filterSize, numFilters, DilationFactor = dilationFactor, ...
Padding = "causal", Name="conv1_" + i) % 一维卷积层
layerNormalizationLayer % 层归一化
spatialDropoutLayer(dropoutFactor) % 空间丢弃层
convolution1dLayer(filterSize, numFilters, ...
DilationFactor = dilationFactor, Padding = "causal") % 一维卷积层
layerNormalizationLayer % 层归一化
reluLayer % 激活层
spatialDropoutLayer(dropoutFactor) % 空间丢弃层
additionLayer(2, Name = "add_" + i)]; % 合并层
lgraph = addLayers(lgraph, layers); % 将卷积层加入到网络
lgraph = connectLayers(lgraph, outputName, "conv1_" + i); % 将模块的卷积层首层和残差结构连接
% 残差连接 -- 首层
if i == 1
layer = convolution1dLayer(1, numFilters, Name = "convSkip"); % 建立残差卷积层
lgraph = addLayers(lgraph, layer); % 将残差卷积层加入到网络
lgraph = connectLayers(lgraph, outputName, "convSkip"); % 将残差卷积层
lgraph = connectLayers(lgraph, "convSkip", "add_" + i + "/in2"); % 将跳跃残差层连接到 addtion 层 输入口2
else
lgraph = connectLayers(lgraph, outputName, "add_" + i + "/in2"); % 将残差层连接到 addtion 层 输入口2
end
% 更新输出层名字
outputName = "add_" + i;
end
%% 网络输出层
layers = [
fullyConnectedLayer(num_class, Name = "fc")
softmaxLayer
classificationLayer];
lgraph = addLayers(lgraph, layers); % 将输出层加入到网络
lgraph = connectLayers(lgraph, outputName, "fc"); % 连接输出层到网络最后
%% 查看网络结构
figure
plot(lgraph)
title("Temporal Convolutional Network")
set(gcf,'color','w')
%% 设置训练参数
options = trainingOptions("adam", ... % 优化器Adam
InitialLearnRate = 5e-3, ... % 初始学习率为0.005
MaxEpochs = 100, ... % 最大训练次数
LearnRateSchedule = 'piecewise', ... % 学习率下降
LearnRateDropFactor = 0.8, ... % 学习率下降因子 0.8
LearnRateDropPeriod = 50, ... % 经过50次训练后 学习率为 0.005*0.8
Plots = "training-progress", ... % 显示训练曲线
Verbose = 0); % 关闭命令行
%% 训练网络
[net, Loss] = trainNetwork(pc_train, tc_train, lgraph, options);
%% 仿真预测
t_sim1 = predict(net, pc_train);
t_sim2 = predict(net, pc_test );
%% 反归一化
for i = 1: M
T_sim1(i) = cell2mat(vec2ind(t_sim1(i)));
end
for i = 1: N
T_sim2(i) = cell2mat(vec2ind(t_sim2(i)));
end
%% 性能评价
error1 = sum((T_sim1 == T_train)) / M * 100 ;
error2 = sum((T_sim2 == T_test )) / N * 100 ;
%% 绘图
figure
plot(1: M, T_train, 'r-*', 1: M, T_sim1, 'b-o', 'LineWidth', 1)
legend('真实值', 'TCN预测值')
xlabel('预测样本')
ylabel('预测结果')
string = {'训练集预测结果对比'; ['准确率=' num2str(error1) '%']};
title(string)
grid
set(gcf,'color','w')
figure
plot(1: N, T_test, 'r-*', 1: N, T_sim2, 'b-o', 'LineWidth', 1)
legend('真实值', 'TCN预测值')
xlabel('预测样本')
ylabel('预测结果')
string = {'测试集预测结果对比'; ['准确率=' num2str(error2) '%']};
title(string)
grid
set(gcf,'color','w')
%% 混淆矩阵
if flag_conusion == 1
figure
cm = confusionchart(T_train, T_sim1);
cm.Title = 'Confusion Matrix for Train Data';
cm.ColumnSummary = 'column-normalized';
cm.RowSummary = 'row-normalized';
set(gcf,'color','w')
figure
cm = confusionchart(T_test, T_sim2);
cm.Title = 'Confusion Matrix for Test Data';
cm.ColumnSummary = 'column-normalized';
cm.RowSummary = 'row-normalized';
set(gcf,'color','w')
end
%% 损失函数曲线
figure
subplot(2, 1, 1)
plot(1 : length(Loss.TrainingAccuracy), Loss.TrainingAccuracy, 'r-', 'LineWidth', 1)
xlabel('迭代次数')
ylabel('准确率')
legend('训练集正确率')
title ('训练集正确率曲线')
grid
set(gcf,'color','w')
subplot(2, 1, 2)
plot(1 : length(Loss.TrainingLoss), Loss.TrainingLoss, 'b-', 'LineWidth', 1)
xlabel('迭代次数')
ylabel('损失函数')
legend('训练集损失值')
title ('训练集损失函数曲线')
grid
set(gcf,'color','w')