基于Matlab的CNC刀具RUL预测：特征提取与神经网络实时回归实践

简介：面向CNC机床状态监测与剩余寿命（RUL）预测需求，这套基于Matlab的代码资源实现了刀具状态的实时监测与寿命预估功能。系统采用参数化编程，参数可灵活修改，代码注释明细，附赠可直接运行的案例数据，特别适合机械电子、计算机、数学等专业学生用于课程设计、期末大作业和毕业设计。

压缩包共47个文件，以33个.m主程序为核心，辅以5个.py脚本、2个.mat案例数据、2个.csv数据表及3张运行效果图，整体大小仅2.08MB，结构清晰、部署便捷。运行程序即可复现从信号采集、特征提取到模型预测的完整流程。

已有58人学习使用这套资源。对于希望借实际工程案例提升Matlab编程与状态监测建模能力的读者，这套代码提供了直观的参考实现，也方便在此基础上进一步改进和扩展。

1. 为什么刀具RUL预测在CNC机床上是典型的实时回归问题

数控车间里，刀具崩刃往往是毫无先兆的：上一刀表面质量还合格，下一刀切削力突变，工件直接报废，严重时连主轴和刀塔一起损伤。刀具剩余使用寿命（RUL）预测的核心就是把这个"毫无先兆"变成"提前半小时知道"。问题本质上是一个回归任务：把振动、声音、温度等传感信号映射成连续时间值（还能加工多少分钟）。难处不在用什么网络，而在三点：信号本身强噪声、机床工况频繁切换、离线训练好的模型放到在线滑窗上很容易失真。

这套基于Matlab实现的代码资源把特征提取和神经网络预测串成了一条完整链路，同时给出案例数据，不用自己采集就能把整个流程跑起来。它针对Matlab 2014/2019a/2024a做了适配，参数采用集中定义的方式，适合做课程设计、毕业设计，也适合设备维护方向的工程预研。下面按"特征提取→模型训练→实时系统→验证技巧"拆开讲，每一段都能独立复用。

2. 从振动信号到特征向量：Matlab特征提取模块拆解

2.1 为什么只取时域和频域特征而不直接灌原始波形

很多人拿到振动信号第一反应是塞给LSTM或一维CNN。对离线竞赛数据这样可以，但在CNC实时监测里要谨慎。原始波形采样率如果到20kHz，做一次2048点的窗口，每秒就要处理接近10帧数据，直接输入原始波形的模型参数量大，而且现场工控机的Matlab环境不一定有GPU。把信号压缩成低维特征后，喂一个只有十来个神经元的前馈网络就能跑出可用的RUL预测，这对实时性和部署都是友好的。

另一个理由是特征可解释。现场工程师不会只信一个黑盒输出，他们需要知道预测值到底是"振动幅值涨了"还是"频带能量变了"导致的。时域特征对应切削力波动，频域特征对应主轴转频和刀齿通过频率附近的能量分布，这些和刀具磨损机理是对应的。特征提取阶段保留这些物理含义，模型输出异常时才能回溯。

2.2 可直接改动的特征提取函数

下面这个函数包含去直流、时域特征和频域能量比，全部用Matlab基础函数实现，不依赖额外工具箱，2014a也能直接运行：

function feat = extractFeatures(sig, fs, bank) % sig: 一帧振动信号，列向量 % fs: 采样频率(Hz) % bank: 结构体，含 f_low 和 f_high 两个频带边界 seg = sig(:) - mean(sig); % 去直流成分，避免FFT泄漏 N = length(seg); % 时域特征 rms_val = sqrt(mean(seg.^2)); % 均方根，反映振动能量 peak_val = max(abs(seg)); % 峰值，捕捉冲击 std_val = std(seg); % 标准差 kurt_val = mean((seg - mean(seg)).^4) / (std_val^4 + eps); % 峭度 crest_f = peak_val / (rms_val + eps); % 峰值因子 % 频域特征：目标频带能量占总能量比例 Y = fft(seg); P2 = abs(Y(1:floor(N/2)+1)).^2; % 单边功率谱 f = fs * (0:floor(N/2)) / N; idxBand = (f >= bank.f_low) & (f <= bank.f_high); bandRatio = sum(P2(idxBand)) / (sum(P2) + eps); feat = [rms_val, peak_val, kurt_val, crest_f, bandRatio]; end

这段代码的逻辑是先把振动信号去直流，否则FFT会在零频位置出现一个很大的分量，把其他频段能量都压下去。时域特征中，RMS比峰值的鲁棒性好，峰值的幅值容易受偶然电磁干扰影响，所以后面用峰值因子把峰值归一化到RMS上，既保留冲击信息又削弱量纲差异。峭度计算时用std_val^4 + eps防止静止状态下标准差为零造成除零错误。

使用这个函数时，bank结构体需要手动指定。比如主轴转频30Hz、三齿刀铣削，则边带集中在大约90Hz附近，bank.f_low=85、bank.f_high=95。如果不知道具体转频，可以先采一段正常信号画频谱图，找到峰值位置后再定频带。特征输出的维度固定为5，后接神经网络输入层节点数时必须对应。

2.3 特征随磨损的变化趋势

不同特征对磨损的敏感区间不一样，理解这一点对后续调模型阈值很重要。下表汇总了常见特征在刀具从新刀到报废过程中的典型趋势：

RMS是所有特征里最稳定的，但也最容易受到切削参数改变的影响。比如进给量加大时RMS会立刻变大，这种变化和磨损无关。峭度和峰值因子在磨损初期上升很快，因为刀具微崩刃会产生冲击脉冲；但当刀具进入剧烈磨损期，连续大振幅信号占主导，冲击相对淡化，峭度反而下降。所以不能用单一特征做阈值报警，必须靠神经网络来组合使用。

2.4 特征工程中的常见误用

一个常见误用是直接把原始信号的一段统计值当成特征，却不考虑时间对齐。离线训练时数据是完整生命周期，每个特征帧对应一个RUL标签；在线预测时特征帧是边采边算的，如果直接调用extractFeatures而没有维护滑窗，标签和信号会错位。建议把函数改成接收固定长度sig的调用方式，由上层滑窗负责切帧，函数内不要自己等待数据。

另一个误用是在时域特征里加入温度信号却不做时间常数匹配。温度传感器响应慢，振动信号变化快，同一时刻的温度其实反映的是几分钟前的热量积累。如果非要用温度，需要在提取特征前对温度做滞后补偿，或者干脆先用振动特征建模，温度只作为后台参考。

3. 神经网络训练与RUL预测模型：从特征到剩余寿命

3.1 训练标签怎么构造

监督学习需要每个样本对应一个RUL标签。假设一把刀具从安装到报废总共能加工T分钟，在第t分钟取了一段信号，那么这段信号对应的标签就是 T - t。实际项目里更常见的做法是用加工件数：一把刀能加工N个零件，当前加工完第k件，RUL标签就是 N-k。注意RUL标签必须是"剩余量"，而不是"已磨损量"，两个量相差一个常数，训练出来的网络虽然预测结果可以换算，但误差分布会不一样，直接预测剩余量更直观。

如果案例数据没有给完整的生命周期，只给了几段不同磨损状态的信号，就别硬套回归。可以把阶段标签变成伪RUL：比如早期磨损、中期磨损、晚期磨损分别标为90、50、10，这样网络学到的是相对排序，不是绝对时间。等到有完整生命周期数据后再换成真实分钟数。这套Matlab资源附带的是完整案例数据，所以可以直接按真实RUL标签训练。

3.2 网络设计与训练代码

选型上，输入是5个特征，输出是1个RUL值，数据规模通常几千到几万帧，用两层隐藏层的前馈网络已经足够。代码用Matlab神经网络工具箱的fitnet，它默认做回归，不需要手动改激活函数。下面是可直接替换的训练主脚本：

% featMat: nFeatures x nSamples 特征矩阵 % rulLabels: 1 x nSamples RUL标签向量 rng(0); % 固定随机种子，让结果可复现 net = fitnet([12 6], 'trainlm'); % 12和6分别是两个隐藏层神经元数 net.divideFcn = 'dividerand'; % 随机划分训练/验证/测试 net.divideParam.trainRatio = 0.7; net.divideParam.valRatio = 0.15; net.divideParam.testRatio = 0.15; net.trainParam.epochs = 500; % 最大训练轮数 net.trainParam.showWindow = false; % 关闭训练窗口，适合批处理 net.trainParam.min_grad = 1e-6; % 梯度低于阈值则停止 [net, tr] = train(net, featMat, rulLabels); % 用测试集评估 pred = net(featMat(:, tr.testInd)); trueVal = rulLabels(tr.testInd); rmse = sqrt(mean((pred - trueVal).^2)); r2 = 1 - sum((pred - trueVal).^2) / sum((trueVal - mean(trueVal)).^2); fprintf('RMSE=%.2f min, R^2=%.3f\n', rmse, r2);

dividerand会把所有样本随机打乱后按比例分成三份，这样能保证训练集和测试集分布一致。trainlm是Levenberg-Marquardt算法，收敛快，适合样本量中等、网络层数浅的场景。min_grad设成1e-6比默认值更严格，可以延长一点训练时间，但能避免在误差曲面鞍点上提前退出。

代码里featMat(:, tr.testInd)的索引方式要注意：fitnet内部已经用divideParam完成了划分，tr.testInd保存的是全局索引，直接用这个索引取原始特征矩阵即可，不需要手动切出测试集。

3.3 训练参数调整要点

如果训练集只有几百个样本，把隐藏层从[12 6]改成[8 4]，并改用trainbr贝叶斯正则化，能有效抑制过拟合。trainbr训练速度比trainlm慢，但小样本下RUL预测的泛化能力明显更好。反过来，如果样本量有几万帧，12个节点可能不够，可以把第一层加宽到20个节点，但第二层不要超过8个，否则参数量增大反而引入噪声。

特征标准化也是一个容易被忽略的坑。fitnet内部会对输入做标准化，但它保存的是训练集的均值和方差。在线预测时如果直接对原始特征调用net，即使Matlab自动处理了，也要确认新数据的特征量纲和训练集一致。我习惯在训练前自己做一遍标准化，并把均值和标准差保存到scaler结构体里，在线预测时先手动标准化再输入网络，这样换版本、换电脑都不会出问题。

4. 实时监测系统：滑窗参数与版本兼容

4.1 在线预测的代码骨架

离线训练完成后，实时系统要做的事就是循环：取一段信号、算特征、标准化、预测RUL、判断告警。下面的代码模拟了从原始信号流里滑窗读取的场景，example_stream.mat是一个模拟的连续振动数据文件，实际项目里可替换成数据采集卡读取函数：

% 初始化参数 P.fs = 20000; % 采样率(Hz) P.windowLen = 2048; % 滑窗长度(采样点数) P.hop = 512; % 滑窗步长，越小实时性越高 P.band.low = 85; % 频带下界 P.band.high = 95; % 频带上界 P.rulWarn = 30; % RUL低于30分钟时报警 load('trained_net.mat', 'net', 'scaler'); % 加载训练好的模型和标准化参数 stream = load('example_stream.mat', 'sig'); sigStream = stream.sig(:); idx = 1; while idx + P.windowLen - 1 <= length(sigStream) frame = sigStream(idx : idx + P.windowLen - 1); feat = extractFeatures(frame, P.fs, P.band); featScaled = (feat - scaler.mu) ./ scaler.sigma; % 手动标准化 rulNow = net(featScaled); rulNow = max(rulNow, 0); % 预测值不为负 if rulNow < P.rulWarn warning('当前RUL约为 %.1f 分钟，请准备换刀', rulNow); end idx = idx + P.hop; pause(P.hop / P.fs); % 模拟实时采集间隔 end

这段代码的核心是滑窗步长P.hop。windowLen决定了单次特征提取需要累积多少数据点，hop决定了窗口每次后移多少点。窗口重叠越多，RUL预测曲线越平滑，但计算量也越大。feat - scaler.mu里的mu和sigma是训练时按每个特征维度计算出来的，必须和训练特征完全一致，否则预测值会系统性偏移。

网络输出的rulNow是一个数值，可能在极端情况下出现负值或远大于新刀寿命的值。加一句max(rulNow,0)只解决了负值问题，上限最好也做饱和处理，比如取min(rulNow, T_max)，T_max可以在参数区定义。

4.2 关键参数速查表

采样率和窗口长度不是独立参数。频率分辨率大约是fs / windowLen，比如20kHz除以2048约等于9.8Hz，这意味着频带边界精度不到10Hz。如果主轴转频是30Hz，三齿刀通过频率是90Hz，边带范围85~95Hz刚好能覆盖，但换到5齿刀时通过频率变成150Hz，原来的频带就完全失效了。

4.3 Matlab不同版本适配

这套代码资源标注兼容Matlab 2014、2019a和2024a。跨版本最容易踩的坑是两类函数：一是normalize这类后来引入的函数，在2014a里不存在，所以标准化要用(x - mu)./sigma手工计算；二是深度网络工具箱的函数如trainNetwork在2014a里没有，所以代码用经典fitnet，它的接口从2013b开始就没变过。只要在代码里避开tall、normalize、reshape的新增调用形式，老版本也能跑。

另一个兼容性问题出现在中文路径和文件编码上。Matlab 2014a在Windows下对UTF-8编码的.m文件支持不友好，中文注释可能乱码。如果遇到这种情况，把.m文件用Matlab编辑器重新保存为GBK编码，或者把代码里的中文注释改成英文注释。参数名用拼音缩写也可以，关键是保证代码逻辑不变。

4.4 直接运行案例数据的流程

资源里Codes_Feature_extraction和Matlab_Feature_extraction & NN分别对应特征提取和网络训练两大块。建议先打开Codes_Feature_extraction目录，运行里面的特征提取脚本，让它生成特征矩阵和标签并保存到.mat文件；再打开Matlab_Feature_extraction & NN目录，运行网络训练脚本，得到模型文件。最后把模型文件和上节的实时预测脚本放在同一目录，就可以直接跑通实时监测流程。

要注意目录名里有空格，Matlab对路径空格的处理有时会出现引号错误。遇到undefined function错误时，优先检查当前目录是否在路径里，而不是怀疑代码本身。案例数据文件不要改名，脚本里如果用了相对路径，改名会导致load失败。

5. 验证RUL模型时容易忽视的三个细节

5.1 不要随机打乱时间序列

前文训练用的dividerand是随机划分，这在特征独立同分布时没问题。但刀具磨损数据本质是时间序列，相邻窗口的特征几乎一致，如果训练集和测试集互相穿插，模型相当于见过测试样本前后的信息，评估指标会虚高。工程化验证必须改成时间顺序划分：前70%训练，中间15%验证，最后15%测试，模拟"用过去预测未来"。修改方法很简单：

net.divideFcn = 'divideblock'; net.divideParam.trainRatio = 0.7; net.divideParam.valRatio = 0.15; net.divideParam.testRatio = 0.15;

换成divideblock后，训练集是整段数据的前70%，测试集是最后15%。这时候测试集里包含刀具磨损最剧烈的阶段，RUL预测误差通常会比随机划分大不少，但这才是真实在线场景的水平。

5.2 预测结果要做指数平滑

实时滑窗每步都会输出一个RUL，单帧特征受切屑、冷却液冲击影响，预测值会在真实值附近剧烈抖动。直接用这个值触发报警会发生误报。常见做法是加一阶指数平滑：

alpha = 0.1; rulSmooth = alpha * rulNow + (1 - alpha) * rulSmoothPrev;

alpha越小越平滑，但滞后越严重。刀具寿命通常还有几十分钟时才报警，滞后几十秒无碍，所以alpha取0.05到0.15之间比较合适。如果想让系统在刀具快报废时反应更快，可以用动态alpha，当rulNow < 2 * P.rulWarn时把alpha提高到0.3。

5.3 换工况前先校准频带

同一把刀加工不同材料、不同转速时，目标频带完全不同。如果机床有主轴转速输出，可以在线计算刀齿通过频率，并自动更新P.band。最省事的校准方法是从当前信号频谱里找峰值作为主频：

P2seg = abs(fft(frame)).^2; fAxis = (0:floor(P.windowLen/2)) * P.fs / P.windowLen; searchIdx = fAxis > 50 & fAxis < 1000; [~, peakIdx] = max(P2seg(searchIdx)); fSpindle = fAxis(searchIdx); fPass = fSpindle(peakIdx); % 如果是多齿刀，此值可能是刀齿通过频率 P.band.low = fPass * 0.85; P.band.high = fPass * 1.15;

这段代码在每次更换工况后运行一次，把频带中心自动对到当前频谱最强峰上。要注意的是，频谱最强峰也可能来自主轴轴承故障，所以校准后最好人工看一眼峰值位置是否符合当前转速设定。把校准逻辑写成一个独立函数，在线运行时每隔几分钟或收到换刀信号时调用一次，RUL预测准确率会比固定频带高出一截。

本文还有配套的精品资源，点击获取