简介面向机械故障诊断与智能运维领域的MATLAB项目实例完整实现了基于双向门控循环单元BiGRU的轴承剩余寿命RUL预测流程。资源涵盖数据预处理、退化特征构造、健康指数与RUL标签生成、模型搭建与训练、预测评估以及可交互的GUI界面设计适合具备一定MATLAB和信号处理基础的科研人员、高校研究生及工业界从业者参考。压缩包内含1个docx文档包体约112KB文档以目录化结构组织包含项目背景、模型架构、代码示例及应用领域等章节便于按需查阅。目前已有75人学习适用于风力发电、轨道交通、航空发动机等关键设备的寿命预测与健康管理场景。通过该文档读者可获得从原始振动信号到可部署模型的全套设计思路、关键代码片段、超参数调整方向与部署建议可作为理解BiGRU在时间序列回归任务中应用的工程化教学案例。1. 为什么轴承剩余寿命预测绕不开 BiGRU轴承是旋转机械里最容易失效的部件它的剩余使用寿命RUL预测直接决定设备是“定时换”还是“按状态换”。传统方法依赖振动信号的时域统计量、频域包络或时频域峭度再用回归模型拟合退化趋势但轴承退化过程是非线性的而且早期微弱冲击容易被工况噪声淹没。基于双向门控循环单元BiGRU的做法是把同一段振动信号按正时序和反时序同时建模让网络同时看到“过去怎么退化”和“未来趋势指向哪里”。相比单向 GRUBiGRU 在捕捉退化拐点上的误差更低相比 LSTM它的参数更少在 MATLAB 里跑 CPU 也能在可接受时间内完成一轮训练。这篇文章适合三类人想用深度学习做设备健康管理PHM的研究生、需要在 MATLAB 里交付一个带界面演示的剩余寿命估计原型的工程师、以及看过不少 BiGRU 论文却始终调不出好效果的自学者。我会把从原始振动数据到 GUI 显示预测曲线的完整链路拆开讲重点放在样本构造、网络搭建、训练选项和界面集成上。整个过程不需要 Python也不需要额外的深度学习框架MATLAB 2022b 以后的版本就能直接跑通。2. BiGRU 原理与轴承退化数据怎么组织成训练样本2.1 为什么双向门控能比单向 GRU 多学到半拍趋势GRU 的核心是重置门和更新门。重置门决定过去状态有多少被遗忘更新门决定新输入和旧状态各占多少比例。单向 GRU 在 t 时刻的状态只依赖 t 之前的输入这符合时序预测的基本假设但轴承退化的故障特征往往不是单调变化的。比如外圈剥落初期振动幅值可能先小幅回落再持续攀升单向 GRU 在拐点附近容易产生滞后。BiGRU 的解决办法很简单一层正向 GRU 从头到尾读序列一层反向 GRU 从尾到头读序列然后把两个方向的隐藏状态在每一步拼接成一个特征向量。这样在 t 时刻模型既有 t 之前的历史信息也有 t 之后的“未来”信息——训练时这个未来是已知的预测时我们也是用一整段窗口内的信号来估计当前 RUL所以并不违反因果关系。MATLAB 里实现双向循环网络不用手动写两层。深度学习工具箱提供了bilstmLayer它是 LSTM 的双向版本但如果你要的是 BiGRU常见做法是用lstmLayer配合SequenceInputLayer之后再拼接其实不准确。让我明确一个坑MATLAB 深度学习工具箱直到 2023b 之前bilstmLayer只支持 LSTM没有直接的 BiGRU 层。要在 MATLAB 里实现 BiGRU需要自己写一个自定义层或者用两个gruLayer分别设置Direction参数实际上gruLayer也接受Direction你可以在gruLayer(32,Direction,bidirectional)直接创建双向 GRU 层。这一点很多人不知道因为文档里Direction参数明确支持bidirectional。所以办法很简单用gruLayer加Direction设置就得到 BiGRU 层输出大小是隐藏单元数的两倍。2.2 原始振动信号不能直接喂给 BiGRU先做窗口化轴承剩余寿命预测的数据集常用 PHM2012 或 XJTU-SY。这些数据集里的原始信号是几十万点的连续振动波形如果直接整段输入序列长度太长训练慢而且容易过拟合。我一般先把振动信号按固定长度切窗再用滑动窗口截取最近一段时间的历史数据作为一个样本。假设信号采样率是 25.6kHz每 10 分钟一段记录那么一段就有 15360 个点。把这 1.5 万个点全部作为输入特征网络参数量会大得离谱而且高维输入里噪声占比高。所以常见的做法是先对每个时间片段提取统计特征把高维信号压缩成特征向量序列。每个片段比如每 10 分钟提取 8 个特征均方根RMS、峭度、峰值因数、波形因数、脉冲因数、裕度因数、偏度、标准差。这样一条完整生命周期数据就变成一个长度等于片段数、每个片段 8 维特征的时间序列。BiGRU 接收的输入形状是特征维度 × 时间步数在 MATLAB 里需要格式化为numFeatures × numTimeSteps的 cell 数组每个元素。下面是一段提取特征并构造序列的示例代码function features extractFeatures(signal, fs) % 输入: 一段振动信号 signal, 采样频率 fs % 输出: 9维特征向量 N length(signal); t (0:N-1)/fs; rms_val sqrt(mean(signal.^2)); peak max(abs(signal)); % 峭度: 四阶中心矩除以标准差四次方 kurt kurtosis(signal); % 波形因数: RMS 与整流平均值的比值 mean_abs mean(abs(signal)); shape_factor rms_val / (mean_abs eps); % 峰值因数 crest_factor peak / (rms_val eps); % 脉冲因数 impulse_factor peak / (mean_abs eps); % 裕度因数: 峰值与平方根振幅之比 margin_factor peak / (mean(sqrt(abs(signal)))^2 eps); % 偏度 skew skewness(signal); % 标准差 std_val std(signal); features [rms_val, kurt, shape_factor, crest_factor, ... impulse_factor, margin_factor, skew, std_val, peak]; end这段代码在 MATLAB 里直接保存为函数文件即可。提取的特征中RMS 和峭度是轴承退化最敏感的两个指标RMS 随退化整体上升峭度在早期故障出现时会有明显突增。注意最后我加了峰值和标准差这两个冗余特征它们和前面的因数有相关性但 BiGRU 对这种冗余并不敏感保留它们反而能增强对冲击型故障的辨识。所有特征都要在训练集上做归一化否则 RMS 的量纲是 m/s²峭度是无量纲数数值范围差三个数量级模型会偏向大数值特征。2.3 剩余寿命标签怎么定不要用线性下降剩余寿命标签如果直接用“当前时间到失效时间”的线性函数模型会很难学因为早期阶段退化特征不明显而 RUL 却从几千分钟开始线性下降网络被迫在小特征变化时输出大数值差。工程上常用分段线性 RUL 函数设定一个最大寿命阈值L_max当真实剩余寿命大于 L_max 时标签统一为 L_max小于 L_max 时标签等于真实剩余寿命。比如 L_max 取 150 分钟那么轴承刚开始运行时标签恒为 150当剩余寿命真的大于 150 时标签保持 150只有最后 150 分钟才线性下降到 0。这样做的好处是降低早期样本的回归难度让模型把注意力集中在退化加速的后期。L_max 的取值要看工况转速越高、载荷越大L_max 应该越小一般在 100 到 300 之间。% rul: 真实剩余寿命向量, 单位分钟 % L_max: 分段阈值 function label piecewiseRul(rul, L_max) label min(rul, L_max); % 将标签归一化到 [0,1], 方便网络输出使用 sigmoid label label / L_max; end这里把标签归一化到 0 到 1 之间配合输出层用sigmoid层模型输出就落在这个范围内。推理时再乘以 L_max 得到分钟数。注意训练时损失函数默认回归用均方误差如果标签不归一化大数值标签会导致梯度爆炸归一化后 MSE 的数值范围也小得多。3. MATLAB 中搭建 BiGRU 网络的最小可运行代码3.1 网络架构与层参数选择在 MATLAB 里搭建 BiGRU 网络不需要自定义训练循环直接用trainNetwork配合dlnetwork或layers都可以。我比较推荐用dlnetworktrainnet的新工作流因为trainNetwork对序列输入的支持在旧版本里有尺寸限制。不过为了代码简洁以下用trainNetworklayerGraph的经典方式它仍然能跑通。一个常见的 BiGRU 回归网络结构如下inputSize 9; % 特征维度 numHiddenUnits 64; % BiGRU 隐藏单元数 numClass 1; % 输出 RUL 值 layers [ sequenceInputLayer(inputSize, Normalization, none, Name, input) gruLayer(numHiddenUnits, Direction, bidirectional, Name, bigru1) dropoutLayer(0.2, Name, dropout) fullyConnectedLayer(32, Name, fc1) reluLayer(Name, relu1) fullyConnectedLayer(numClass, Name, fc_out) regressionLayer(Name, reg) ]; options trainingOptions(adam, ... MaxEpochs, 100, ... MiniBatchSize, 32, ... InitialLearnRate, 0.001, ... GradientThreshold, 1, ... Shuffle, every-epoch, ... Verbose, true, ... Plots, training-progress);逐层说明sequenceInputLayer这里接受 9 维特征序列。gruLayer(64,Direction,bidirectional)实际会产生两个方向的 GRU每方向 64 个单元输出维度是 128。dropoutLayer放在循环层后面可以抑制过拟合训练时随机置零 20% 的神经元预测时自动关闭。fullyConnectedLayer(32)将 128 维映射到 32 维再用 ReLU 增加非线性。最后全连接输出 1 维regressionLayer计算均方误差。需要注意MiniBatchSize要根据样本量调整。样本量在几千到几万之间时32 是一个安全值如果每个序列的长度很大超过 500 个时间步显存占用会快速上升此时可以把MiniBatchSize降到 16 或者 8。如果训练时出现 loss 为 NaN先检查是不是学习率太大InitialLearnRate大于 0.01 往往直接梯度爆炸。3.2 训练数据的格式cell 数组是最大的坑MATLAB 的trainNetwork在训练序列网络时要求输入数据是 cell 数组每个 cell 是一个numFeatures × numTimeSteps的矩阵而且每个 cell 的列数可以不同。标签则是1 × numObservations的行向量或矩阵。% X: 总样本, 假设有 N 个训练序列 % 每个序列是一个 9 x T_i 的矩阵 numObservations length(signalFeatures); % signalFeatures 是 cell 数组 X signalFeatures; % 每个元素为 9 x T_i Y labels; % 1 x numObservations 每个值为 [0,1] % 检查维度 assert(size(X{1}, 1) 9, 特征维度必须为9); assert(size(Y, 2) numObservations, 标签数量与样本数不一致); % 划分训练/验证集 idx randperm(numObservations); numTrain floor(0.8 * numObservations); XTrain X(idx(1:numTrain)); YTrain Y(idx(1:numTrain)); XVal X(idx(numTrain1:end)); YVal Y(idx(numTrain1:end));这里assert用来提前检查维度错误。常见报错是Invalid training data. Predictors must be a cell array of sequences说明你把数值矩阵直接传进去了。另一个常见问题是每个序列长度 T_i 不同这在 BiGRU 中是允许的因为深度神经网络工具箱会自动填充或打包但要注意你的显卡驱动如果太旧可能触发 cuDNN 的不兼容错误这时候在trainingOptions里设置ExecutionEnvironment,cpu先确认模型逻辑没问题再换回 GPU。3.3 预测阶段从振动信号到 RUL 输出训练完成后对新的测试轴承数据做预测时需要走完全相同的预处理流程。先切窗、提特征、拼接成 9×T 矩阵然后调用predict% newSignal: 新轴承的一段连续振动信号 % fs: 采样率 % 假设我们已经按固定时间间隔切成了多个片段, 每个片段提取特征 features []; for i 1:length(fragments) features [features, extractFeatures(fragments{i}, fs)]; end % 归一化: 使用训练时的均值 std % 这里需要在训练时保存 meanVal, stdVal featuresNorm (features - meanVal) ./ (stdVal eps); % 预测 XTest {featuresNorm}; % 1x1 cell, 每个元素 9 x T predNorm predict(net, XTest); rul predNorm * L_max; % 反归一化得到分钟数这段代码最后一行就是剩余寿命的估计值。注意predict返回的是 cell 数组对于序列回归predict返回的可能是数值矩阵或 cell。如果net是SeriesNetwork或DAGNetworkpredict对单个序列输入返回一个数值矩阵每个时间步对应一个输出。但在我们的网络结构中regressionLayer默认输出是整个序列的最后一个时间步的预测值吗这里有一个关键细节fullyConnectedLayer作用在每个时间步上所以输出层会对序列每一个时间步都输出一个 RUL 预测值。训练时regressionLayer会比较每个时间步的误差但我们实际只用最后一个时间步的值作为最终 RUL。为了得到最后一个时间步的值需要这样写predNorm predict(net, XTest); predLast predNorm(end); % 取最后一个时间步如果predNorm是矩阵size(predNorm)为1 x Tend取到最后一步。若predict返回的是 cell则predNorm{1}(end)。这是新手最容易忽略的坑整个网络是 sequence-to-sequence 输出不是 sequence-to-one。如果只想要序列一个值也可以用sequenceFoldingLayer或自定义输出层但多数情况下取末尾时间步的预测值就够用而且能观察到每个时间步的退化估计变化。4. 用 MATLAB App Designer 把 BiGRU 包成 GUI 演示程序4.1 界面布局模型加载、数据导入、结果展示三块项目标题要求带 GUI 设计和代码详解实际交付时一般用 MATLAB App Designer 做界面。打开 MATLAB在命令行输入appdesigner新建空白应用从左侧拖入组件。我常用的布局是左侧面板“加载振动数据”、“提取特征”、“加载模型”、“开始预测” 四个按钮下方一个Edit Field显示 L_max 值。右侧面板两个坐标轴UIAxes上面的显示原始振动信号下面的显示 RUL 预测曲线。底部一个Label组件显示最终预测结果比如“当前轴承剩余寿命312.5 分钟”。组件名称要设置以app.开头的回调函数比如app.LoadDataButtonPushed。% 按钮回调示例 function LoadDataButtonPushed(app, event) [file, path] uigetfile(*.mat, 选择振动信号文件); if isequal(file, 0) return; end data load(fullfile(path, file)); app.signalData data.signal; % 存到 app 的公共字段 app.Fs data.fs; % 在第一个坐标轴画原始信号 N length(app.signalData); t (0:N-1) / app.Fs; plot(app.UIAxes1, t, app.signalData); title(app.UIAxes1, 原始振动信号); end4.2 把训练好的模型嵌入 GUI存储为 mat 文件训练好的 BiGRU 模型要保存下来供 GUI 调用。在训练脚本中执行save(bigru_rul_model.mat, net, meanVal, stdVal, L_max, fs);注意net里有训练信息文件可能会比较大但通常不超过 20MB因为 GRU 参数量有限。在 GUI 的“加载模型”按钮回调中直接load这个文件并把net等变量存到app对象上function LoadModelButtonPushed(app, event) [file, path] uigetfile(*.mat, 选择训练好的模型文件); if isequal(file, 0) return; end S load(fullfile(path, file)); app.net S.net; app.meanVal S.meanVal; app.stdVal S.stdVal; app.L_max S.L_max; app.Fs S.fs; app.StatusLabel.Text 模型加载完成; end4.3 在 GUI 里执行预测并刷新曲线“开始预测”按钮要做的事先对app.signalData做切窗、提特征再利用app.net预测绘制 RUL 随时间变化的曲线。注意这里的“随时间变化”指的是测试数据本身是有序时间轴我们把每个时间步的预测值画出来。function PredictButtonPushed(app, event) if isempty(app.net) || isempty(app.signalData) app.StatusLabel.Text 请先加载数据和模型; return; end fs app.Fs; sig app.signalData(:); % 窗口长度 10 分钟 winLen fs * 60 * 10; numWin floor(length(sig) / winLen); features zeros(9, numWin); for i 1:numWin segment sig((i-1)*winLen 1 : i*winLen); features(:, i) extractFeatures(segment, fs); end % 归一化 featuresNorm (features - app.meanVal) ./ (app.stdVal eps); % 预测: 输入 cell predSeq predict(app.net, {featuresNorm}); if iscell(predSeq) predNorm predSeq{1}; else predNorm predSeq; end % 取每个时间步的最后一维度? 实际上每个时间步都有输出 rulTime predNorm(end, :) * app.L_max; timeAxis (1:numWin) * 10; % 分钟 plot(app.UIAxes2, timeAxis, rulTime, LineWidth, 2); xlabel(app.UIAxes2, 运行时间 (分钟)); ylabel(app.UIAxes2, 剩余寿命 (分钟)); title(app.UIAxes2, BiGRU 剩余寿命预测结果); currentRul rulTime(end); app.ResultLabel.Text sprintf(当前剩余寿命: %.1f 分钟, currentRul); end这里predNorm(end, :)取的是输出矩阵的最后一行。因为我们的regressionLayer输出是 1 个特征所以每行对应一个时间步。如果输出是1 x T的行向量那么predNorm(end)就是最后时间步。我在代码里统一用predNorm(end, :)应对矩阵情况行向量时也成立。注意timeAxis每个片段 10 分钟如果 numWin 是 50最后预测点就是 500 分钟。4.4 App Designer 回调里常见的报错报错 “未定义变量 app 或函数 app.myData”说明你没有把变量存到app的公共属性。需要在 App Designer 编辑器里点击“属性”按钮定义properties (Access public)比如net []然后才能在回调里使用app.net。报错 “输入参数太多”回调函数签名应该是function PredictButtonPushed(app, event)注意 MATLAB 会自动传两个参数。如果你在代码中手动调用该函数也要传两个参数。predict出现“使用 Parallel pool”卡住在预测前关闭并行池或者在trainingOptions里设置UseParallel, false。5. 收敛不理想时的三个关键调参点与验证技巧5.1 看训练曲线而不是只看 RMSE用trainingOptions里的Plots回调可以在训练过程中实时看到损失曲线。如果Training loss降到 0.1 附近后不再下降而Validation loss开始上升说明过拟合。此时优先增加dropoutLayer的比例到 0.3或增大L2Regularization到0.0001。如果训练 loss 一直不降先检查特征提取是否有问题把 RMS 或峭度时序画出来看有没有明显的退化趋势。如果特征曲线都在波动模型不可能学到单调递减的 RUL 标签。参数调整顺序建议现象优先调整参数调整方向loss 不降InitialLearnRate从 0.001 降到 0.0003过拟合dropoutLayer比例从 0.2 升到 0.4序列太长训练慢MiniBatchSize从 32 降到 8或缩短时间步预测值偏高L_max减小到 100 附近预测值波动大numHiddenUnits从 64 升到 1285.2 滑动窗口步长与多步预测验证最后再给一个进阶验证技巧。不要把整段测试数据一次性预测而是用滑动步长为 1 个片段的方式从测试数据中间某一点开始连续预测未来每 10 分钟的 RUL然后把预测结果和真实分段标签画在同一张图上。这个验证方法更接近在线监测场景每采集完 10 分钟数据就更新一次 RUL 估计。% 在线预测模拟: 每次只使用到当前时刻为止的数据 startIdx 20; % 从第20个片段开始 numSteps 50; predHist zeros(1, numSteps); trueHist zeros(1, numSteps); for i 1:numSteps currentFeatures featuresNorm(:, 1:startIdx i - 1); predSeq predict(app.net, {currentFeatures}); predHist(i) predSeq(end) * app.L_max; % 真实RUL: 假设已知总生命周期片段数 totalWin trueHist(i) (totalWin - (startIdx i - 1)) * 10; end figure; plot((1:numSteps)*10, predHist, r-o); hold on; plot((1:numSteps)*10, trueHist, b--); xlabel(推进时间 (分钟)); ylabel(剩余寿命 (分钟)); legend(BiGRU预测, 真实RUL, Location, best);这段代码模拟了在线预测随着新数据不断到来输入序列长度每次增加一列。预测结果会随着更接近失效点而逐渐逼近真实值但通常会在早期高估、后期低估这是正常现象。如果你的预测曲线在后期出现突然跳变说明特征提取中某个时刻的峭度或峰值因数出现异常尖峰可以考虑对特征序列做中值滤波featuresNormSmoothed movmedian(featuresNorm, 3, 2);注意这里的移动中值滤波要在归一化之后做而且窗口不宜太大否则会抹掉退化拐点。如果你有多个测试轴承的数据可以计算平均绝对误差MAE和 PHM 竞赛常用的打分函数Score sum(exp(-d/13)-1)当 d0sum(exp(d/10)-1)当 d0这个分数对早期预测d0的惩罚更重比 RMSE 更符合工程需求。在论文或项目汇报里同时给出 RMSE 和 Score 两个指标比单看一条拟合曲线更有说服力。本文还有配套的精品资源点击获取