利用振动信号状态监测和预测gydF4y2Ba

这个例子展示了如何从滚珠轴承从振动信号中提取特征,进行健康监测,并进行预测。这个示例使用功能从信号处理工具箱™和系统辨识工具箱™,并且不需要预测维修工具箱™。gydF4y2Ba

数据描述gydF4y2Ba

存储在负载振动数据gydF4y2BapdmBearingConditionMonitoringData.matgydF4y2Ba(这是一个大型数据集~ 88 mb的下载MathWorks支持文件网站)。金宝app数据存储在一个单元阵列,是使用一个球轴承信号模拟器生成与单点缺陷轴承的外环。它包含多个部分的轴承振动信号模拟在不同的健康状况(从上面3嗯缺陷深度增加3毫米)。每个段存储信号的采样率为1秒收集20 kHz。在pdmBearingConditionMonitoringgydF4y2BaData.matgydF4y2Ba,gydF4y2BadefectDepthVecgydF4y2Ba商店有缺陷深度和时间而变化gydF4y2BaexpTimegydF4y2Ba商店中相应的时间分钟。gydF4y2Ba

url =gydF4y2Ba“//www.tatmou.com/金宝appsupportfiles/predmaint/condition-monitoring-and-prognostics-using-vibration-signals/pdmBearingConditionMonitoringData.mat”gydF4y2Ba;websave (gydF4y2Ba“pdmBearingConditionMonitoringData.mat”gydF4y2Baurl);负载gydF4y2BapdmBearingConditionMonitoringData.matgydF4y2Ba%定义处理数据点的数量。gydF4y2BanumSamples =长度(数据);gydF4y2Ba%定义采样频率。gydF4y2Bafs = 20 e3;gydF4y2Ba%单位:赫兹gydF4y2Ba

情节如何缺陷深度不同部分的数据的变化。gydF4y2Ba

情节(expTime defectDepthVec);包含(gydF4y2Ba的时间(分钟)gydF4y2Ba);ylabel (gydF4y2Ba“缺陷深度(m)”gydF4y2Ba);gydF4y2Ba

情节的健康和错误的数据。gydF4y2Ba

时间= linspace (0, 1, fs)”;gydF4y2Ba%健康轴承信号gydF4y2Ba次要情节(2,1,1);情节(时间、数据{1});包含(gydF4y2Ba“时间(s)”gydF4y2Ba);ylabel (gydF4y2Ba“加速度(m / s ^ 2)”gydF4y2Ba);传奇(gydF4y2Ba“健康轴承信号”gydF4y2Ba);gydF4y2Ba%轴承故障信号gydF4y2Ba次要情节(2,1,2);情节(时间、数据{结束});包含(gydF4y2Ba“时间(s)”gydF4y2Ba);ylabel (gydF4y2Ba“加速度(m / s ^ 2)”gydF4y2Ba);传奇(gydF4y2Ba轴承故障信号的gydF4y2Ba);gydF4y2Ba

特征提取gydF4y2Ba

在本节中,代表特性提取每一段的数据。这些特性将被用于健康监测和预测。轴承诊断和预测包括时域特性的典型特征(均方根、峰值信号峰度等)或频域特性(频率、峰值频率,等等)。gydF4y2Ba

在选择使用哪个功能之前,绘制振动信号谱图。可视化在时域或频域和时频域信号可以帮助发现信号的模式表明退化或失败。gydF4y2Ba

首先计算健康轴承的谱图数据。使用一个窗口大小为500数据点和一个重叠的比率为90%(相当于450数据点)。设置点FFT的数量是512。gydF4y2BafsgydF4y2Ba表示前面定义的采样频率。gydF4y2Ba

[~,fvec tvec P0) =光谱图(500450512年数据{1},fs);gydF4y2Ba

P0gydF4y2Ba光谱图,gydF4y2BafvecgydF4y2Ba频率向量和吗gydF4y2BatvecgydF4y2Ba是时候向量。gydF4y2Ba

健康的轴承信号的谱图。gydF4y2Ba

clf;显示亮度图像(tvec fvec, P0)包含(gydF4y2Ba“时间(s)”gydF4y2Ba);ylabel (gydF4y2Ba的频率(赫兹)gydF4y2Ba);标题(gydF4y2Ba“健康轴承信号谱图”gydF4y2Ba);轴gydF4y2BaxygydF4y2Ba

现在情节的振动信号的谱图有错误的模式发展。你可以看到信号能量集中在更高的频率。gydF4y2Ba

[~,fvec tvec Pfinal] =光谱图(数据{}结束,500450512年,fs);显示亮度图像(tvec fvec Pfinal)包含(gydF4y2Ba“时间(s)”gydF4y2Ba);ylabel (gydF4y2Ba的频率(赫兹)gydF4y2Ba);标题(gydF4y2Ba轴承故障信号谱图的gydF4y2Ba);轴gydF4y2BaxygydF4y2Ba

自谱图数据从健康和有缺陷的轴承是不同的,具有代表性的特征可以从谱图中提取,用于状态监测和预测。在这个例子中,从谱图中提取平均峰值频率作为卫生指标。谱图的表示gydF4y2Ba $PgydF4y2Ba (gydF4y2Ba tgydF4y2Ba,gydF4y2Ba ωgydF4y2Ba)gydF4y2Ba$ 。峰值频率在每一次实例被定义为:gydF4y2Ba

$PgydF4y2Ba egydF4y2Ba 一个gydF4y2Ba kgydF4y2Ba FgydF4y2Ba rgydF4y2Ba egydF4y2Ba 问gydF4y2Ba (gydF4y2Ba tgydF4y2Ba)gydF4y2Ba =gydF4y2Ba 一个gydF4y2Ba rgydF4y2Ba ggydF4y2Ba 米gydF4y2Ba 一个gydF4y2Ba {xgydF4y2Ba}_{ωgydF4y2Ba} PgydF4y2Ba (gydF4y2Ba tgydF4y2Ba,gydF4y2Ba ωgydF4y2Ba)gydF4y2Ba$

平均峰值频率峰值频率上面定义的平均值。gydF4y2Ba

$米gydF4y2Ba egydF4y2Ba 一个gydF4y2Ba ngydF4y2Ba PgydF4y2Ba egydF4y2Ba 一个gydF4y2Ba kgydF4y2Ba FgydF4y2Ba rgydF4y2Ba egydF4y2Ba 问gydF4y2Ba =gydF4y2Ba \frac{1gydF4y2Ba}{TgydF4y2Ba} {\intgydF4y2Ba}_{0gydF4y2Ba}^{TgydF4y2Ba} PgydF4y2Ba egydF4y2Ba 一个gydF4y2Ba kgydF4y2Ba FgydF4y2Ba rgydF4y2Ba egydF4y2Ba 问gydF4y2Ba (gydF4y2Ba tgydF4y2Ba)gydF4y2Ba dgydF4y2Ba tgydF4y2Ba$

健康滚珠轴承计算的平均峰值频率的信号。gydF4y2Ba

[~,I0] = max (P0);gydF4y2Ba%找出峰值频率。gydF4y2BameanPeakFreq0 =意味着(fvec(钱数))gydF4y2Ba%计算平均峰值频率。gydF4y2Ba

meanPeakFreq0 = 666.4602gydF4y2Ba

健康的轴承振动信号平均峰值频率在650赫兹。现在计算故障轴承的平均峰值频率的信号。平均峰值频率2500赫兹以上变化。gydF4y2Ba

[~,Ifinal] = max (Pfinal);meanPeakFreqFinal =意味着(fvec (Ifinal))gydF4y2Ba

meanPeakFreqFinal = 2.8068 e + 03gydF4y2Ba

检查数据在中间阶段,当缺陷深度不是很大但开始影响振动信号。gydF4y2Ba

[~,fvec tvec Pmiddle] =光谱图(数据{/ 2}结束,500450512年,fs);显示亮度图像(tvec fvec Pmiddle)包含(gydF4y2Ba“时间(s)”gydF4y2Ba);ylabel (gydF4y2Ba的频率(赫兹)gydF4y2Ba);标题(gydF4y2Ba轴承信号谱图的gydF4y2Ba);轴gydF4y2BaxygydF4y2Ba

高频噪声分量分布谱图。这样的现象是混合的影响原始振动和振动引起的小缺陷。准确地计算平均峰值频率,滤波器的数据删除这些高频组件。gydF4y2Ba

应用振动信号的中值滤波去除高频噪声分量以及保持高频率的有用的信息。gydF4y2Ba

dataMiddleFilt = medfilt1(数据{/ 2}结束,3);gydF4y2Ba

中值滤波后的阴谋声谱图。高频分量受到抑制。gydF4y2Ba

[~,fvec tvec Pmiddle] =光谱图(dataMiddleFilt 500450512 fs);显示亮度图像(tvec fvec Pmiddle)包含(gydF4y2Ba“时间(s)”gydF4y2Ba);ylabel (gydF4y2Ba的频率(赫兹)gydF4y2Ba);标题(gydF4y2Ba过滤后的轴承信号谱图的gydF4y2Ba);轴gydF4y2BaxygydF4y2Ba

成功以来平均峰值频率区分健康从错误的球轴承滚珠轴承,从每段数据中提取平均峰值频率。gydF4y2Ba

%定义一个进度条。gydF4y2Bah = waitbar (0,gydF4y2Ba“开始提取特征”gydF4y2Ba);gydF4y2Ba%初始化向量存储提取的平均峰值频率。gydF4y2BameanPeakFreq = 0 (numSamples, 1);gydF4y2Ba为gydF4y2Bak = 1: numSamplesgydF4y2Ba%获得最新的数据。gydF4y2BacurData = {k}数据;gydF4y2Ba%应用中值滤波。gydF4y2BacurDataFilt = medfilt1 (curData 3);gydF4y2Ba%计算光谱图。gydF4y2Ba[~,fvec tvec P_k] =光谱图(curDataFilt 500450512 fs);gydF4y2Ba%计算峰值频率在每一次实例。gydF4y2Ba[~,我]= max (P_k);meanPeakFreq (k) =意味着(fvec(我));gydF4y2Ba%显示进度条指示了多少样品处理。gydF4y2Bawaitbar (k / numSamples h,gydF4y2Ba“特征提取”gydF4y2Ba);gydF4y2Ba结束gydF4y2Ba关闭(h);gydF4y2Ba

提取的平均峰值频率与时间。gydF4y2Ba

情节(expTime meanPeakFreq);包含(gydF4y2Ba的时间(分钟)gydF4y2Ba);ylabel (gydF4y2Ba“意思是峰值频率(赫兹)”gydF4y2Ba);网格gydF4y2Ba在gydF4y2Ba;gydF4y2Ba

状态监测和预测gydF4y2Ba

在本节中,状态监测和预测执行使用一个预定义的阈值和动态模型。对于状态监测,创建一个警报,这个警报触发如果平均峰值频率超过预定义的阈值。预测,确定一个动态模型预测平均峰值频率的值在接下来的几个小时。创建一个报警触发,如果预测意味着峰值频率超过了预定义的阈值。gydF4y2Ba

预测可以帮助我们更好地准备一个潜在的故障,甚至停止机器故障前。考虑平均峰值频率作为时间序列。我们可以估计平均峰值频率和时间序列模型使用模型来预测未来值。使用前200的意思是峰值频率值创建一个初始时间序列模型,一旦10新值,使用过去的100值更新时间序列模型。这批处理模式的更新时间序列模型捕捉瞬间的趋势。更新后的时间序列预测模型用于计算10先行一步。gydF4y2Ba

tStart = 200;gydF4y2Ba%开始时间gydF4y2BatimeSeg = 100;gydF4y2Ba%的长度数据构建动态模型gydF4y2BaforecastLen = 10;gydF4y2Ba%定义预测时间范围gydF4y2BabatchSize = 10;gydF4y2Ba%为更新动态模型定义批量大小gydF4y2Ba

对于预测和状态监测,您需要设置一个阈值来决定何时停止机器。在这个例子中,使用健康的统计数据和故障轴承产生仿真来确定阈值。pdm的gydF4y2BaBearingConditionMonitoringStatistics.matgydF4y2Ba商店平均峰值频率的概率分布为健康和有缺陷的轴承。概率分布是由微扰计算的缺陷深度健康和有缺陷的轴承。gydF4y2Ba

url =gydF4y2Ba“//www.tatmou.com/金宝appsupportfiles/predmaint/condition-monitoring-and-prognostics-using-vibration-signals/pdmBearingConditionMonitoringStatistics.mat”gydF4y2Ba;websave (gydF4y2Ba“pdmBearingConditionMonitoringStatistics.mat”gydF4y2Baurl);负载gydF4y2BapdmBearingConditionMonitoringStatistics.matgydF4y2Ba

情节平均峰值频率的概率分布为健康和有缺陷的轴承。gydF4y2Ba

情节(pFreq pNormal,gydF4y2Ba“g——”gydF4y2BapFreq pFaulty,gydF4y2Ba“r”gydF4y2Ba);包含(gydF4y2Ba“意思是峰值频率(赫兹)”gydF4y2Ba);ylabel (gydF4y2Ba的概率分布gydF4y2Ba);传奇(gydF4y2Ba正常轴承的gydF4y2Ba,gydF4y2Ba“有缺陷的轴承”gydF4y2Ba);gydF4y2Ba

基于这一情节,设置阈值的平均峰值频率2000赫兹区分正常轴承和轴承故障以及轴承的使用最大化。gydF4y2Ba

阈值= 2000;gydF4y2Ba

计算采样时间和其单位转换为秒。gydF4y2Ba

samplingTime = 60 * (expTime (2) -expTime (1));gydF4y2Ba%单位:秒gydF4y2BatsFeature = iddata (meanPeakFreq (1: tStart), [], samplingTime);gydF4y2Ba

情节最初的200年平均峰值频率数据。gydF4y2Ba

情节(tsFeature.y)gydF4y2Ba

情节表明初始数据的结合常数和噪音水平。这个预计最初的轴承是健康和平均峰值频率预计不会大幅改变gydF4y2Ba

确定一个二阶状态空间模型使用第一个200数据点。获得模型规范形式和指定采样时间。gydF4y2Ba

past_sys = ss (tsFeature 2gydF4y2Ba“t”gydF4y2BasamplingTime,gydF4y2Ba“形式”gydF4y2Ba,gydF4y2Ba“规范”gydF4y2Ba)gydF4y2Ba

past_sys =离散时间状态空间模型发现:x (t + Ts) = x (t) + K e (t) y (t) = C x (t) + e (t) = (x1, x2) x1 0 1 x2 0.9605 - 0.03942 C = (x1, x2)日元1 0 K = 0.003656日元x1 -0.003899 x2样品时间:600秒参数化:规范形式指数:2。干扰组件:估计很多免费的系数:4使用“idssdata”、“getpvec”、“getcov”参数及其不确定性。对时域数据状态:估计使用党卫军“tsFeature”。适合估算数据:0.2763%(预测聚焦)消防工程:640年,MSE: 602.7gydF4y2Ba

最初的估计动态模型拟合优度低。拟合优度指标的归一化均方根误差(NRMSE)计算gydF4y2Ba

$NgydF4y2Ba RgydF4y2Ba 米gydF4y2Ba 年代gydF4y2Ba EgydF4y2Ba =gydF4y2Ba 1gydF4y2Ba - - - - - -gydF4y2Ba \frac{为gydF4y2Ba {xgydF4y2Ba}_{tgydF4y2Ba rgydF4y2Ba ugydF4y2Ba egydF4y2Ba} - - - - - -gydF4y2Ba {xgydF4y2Ba}_{pgydF4y2Ba rgydF4y2Ba egydF4y2Ba dgydF4y2Ba} 为gydF4y2Ba}{为gydF4y2Ba {xgydF4y2Ba}_{tgydF4y2Ba rgydF4y2Ba ugydF4y2Ba egydF4y2Ba} - - - - - -gydF4y2Ba 米gydF4y2Ba egydF4y2Ba 一个gydF4y2Ba ngydF4y2Ba (gydF4y2Ba {xgydF4y2Ba}_{tgydF4y2Ba rgydF4y2Ba ugydF4y2Ba egydF4y2Ba})gydF4y2Ba 为gydF4y2Ba}$

在哪里gydF4y2Ba ${xgydF4y2Ba}_{tgydF4y2Ba rgydF4y2Ba ugydF4y2Ba egydF4y2Ba}$ 是真正的价值,gydF4y2Ba ${xgydF4y2Ba}_{pgydF4y2Ba rgydF4y2Ba egydF4y2Ba dgydF4y2Ba}$ 预测的价值。gydF4y2Ba

当评估数据的结合常数和噪声水平,gydF4y2Ba ${xgydF4y2Ba}_{pgydF4y2Ba rgydF4y2Ba egydF4y2Ba dgydF4y2Ba} \approxgydF4y2Ba 米gydF4y2Ba egydF4y2Ba 一个gydF4y2Ba ngydF4y2Ba (gydF4y2Ba {xgydF4y2Ba}_{tgydF4y2Ba rgydF4y2Ba ugydF4y2Ba egydF4y2Ba})gydF4y2Ba$ ,给NRMSE接近0。验证模型,绘制残差的自相关。gydF4y2Ba

渣油(tsFeature past_sys)gydF4y2Ba

看到,残差不相关的和生成的模型是有效的。gydF4y2Ba

使用识别模型gydF4y2Bapast_sysgydF4y2Ba预测平均峰值频率值和计算预测的标准偏差值。gydF4y2Ba

[yF, ~, ~, yFSD] =预测(past_sys, tsFeature forecastLen);gydF4y2Ba

情节的预测价值和置信区间。gydF4y2Ba

tHistory = expTime (1: tStart);forecastTimeIdx = (tStart + 1): (tStart + forecastLen);tForecast = expTime (forecastTimeIdx);gydF4y2Ba%绘制历史数据,预测价值和95%置信区间。gydF4y2Ba情节(tHistory meanPeakFreq (1: tStart),gydF4y2Ba“b”gydF4y2Ba,gydF4y2Ba…gydF4y2BatForecast yF.OutputData,gydF4y2Ba“kx”gydF4y2Ba,gydF4y2Ba…gydF4y2Ba[tHistory;tForecast],阈值*的(1,长度(tHistory) + forecastLen),gydF4y2Ba“r——”gydF4y2Ba,gydF4y2Ba…gydF4y2BatForecast yF.OutputData + 1.96 * yFSD,gydF4y2Ba“g——”gydF4y2Ba,gydF4y2Ba…gydF4y2BatForecast yf.outputdata - 1.96 * yFSD,gydF4y2Ba“g——”gydF4y2Ba);ylim([400, 1.1 *阈值]);ylabel (gydF4y2Ba“意思是峰值频率(赫兹)”gydF4y2Ba);包含(gydF4y2Ba的时间(分钟)gydF4y2Ba);传奇({gydF4y2Ba“过去的数据”gydF4y2Ba,gydF4y2Ba“预测”gydF4y2Ba,gydF4y2Ba故障阈值的gydF4y2Ba,gydF4y2Ba“95% C.I”gydF4y2Ba},gydF4y2Ba…gydF4y2Ba“位置”gydF4y2Ba,gydF4y2Ba“northoutside”gydF4y2Ba,gydF4y2Ba“定位”gydF4y2Ba,gydF4y2Ba“水平”gydF4y2Ba);网格gydF4y2Ba在gydF4y2Ba;gydF4y2Ba

情节表明平均峰值频率的预测价值远低于阈值。gydF4y2Ba

现在作为新数据更新模型参数,并重新评估预测的值。还创建一个报警检查或失败的预测价值超过阈值的信号。gydF4y2Ba

为gydF4y2BatCur = tStart: batchSize: numSamplesgydF4y2Ba%最新功能到iddata对象中。gydF4y2BatsFeature = iddata (meanPeakFreq ((tCur-timeSeg + 1): tCur), [], samplingTime);gydF4y2Ba%的新数据时更新系统参数。用以前的模型gydF4y2Ba%参数作为初始猜测。gydF4y2Basys = ss (tsFeature past_sys);past_sys =系统;gydF4y2Ba%更新状态空间模型的预测输出。也计算gydF4y2Ba%的标准差预测输出。gydF4y2Ba[yF, ~, ~, yFSD] =预测(sys, tsFeature forecastLen);gydF4y2Ba%找到时间相应历史数据和预测价值。gydF4y2BatHistory = expTime (1: tCur);forecastTimeIdx = (tCur + 1): (tCur + forecastLen);tForecast = expTime (forecastTimeIdx);gydF4y2Ba%绘制历史数据,预测平均峰值频率值和95%gydF4y2Ba%置信区间。gydF4y2Ba情节(tHistory meanPeakFreq (1: tCur),gydF4y2Ba“b”gydF4y2Ba,gydF4y2Ba…gydF4y2BatForecast yF.OutputData,gydF4y2Ba“kx”gydF4y2Ba,gydF4y2Ba…gydF4y2Ba[tHistory;tForecast],阈值*的(1,长度(tHistory) + forecastLen),gydF4y2Ba“r——”gydF4y2Ba,gydF4y2Ba…gydF4y2BatForecast yF.OutputData + 1.96 * yFSD,gydF4y2Ba“g——”gydF4y2Ba,gydF4y2Ba…gydF4y2BatForecast yf.outputdata - 1.96 * yFSD,gydF4y2Ba“g——”gydF4y2Ba);ylim([400, 1.1 *阈值]);ylabel (gydF4y2Ba“意思是峰值频率(赫兹)”gydF4y2Ba);包含(gydF4y2Ba的时间(分钟)gydF4y2Ba);传奇({gydF4y2Ba“过去的数据”gydF4y2Ba,gydF4y2Ba“预测”gydF4y2Ba,gydF4y2Ba故障阈值的gydF4y2Ba,gydF4y2Ba“95% C.I”gydF4y2Ba},gydF4y2Ba…gydF4y2Ba“位置”gydF4y2Ba,gydF4y2Ba“northoutside”gydF4y2Ba,gydF4y2Ba“定位”gydF4y2Ba,gydF4y2Ba“水平”gydF4y2Ba);网格gydF4y2Ba在gydF4y2Ba;gydF4y2Ba%显示警报当实际监控变量或超过预期值gydF4y2Ba%故障阈值。gydF4y2Ba如果gydF4y2Ba(任何(meanPeakFreq (tCur-batchSize + 1: tCur) >阈值))disp (gydF4y2Ba“监控变量超过失败的阈值”gydF4y2Ba);gydF4y2Ba打破gydF4y2Ba;gydF4y2BaelseifgydF4y2Ba(任何(yF.y >阈值))gydF4y2Ba%估计系统将达到故障阈值的时候。gydF4y2BatAlarm = tForecast(找到(yF.y >阈值,1));disp ([gydF4y2Ba“估计达到故障阈值”gydF4y2Banum2str (tAlarm-tHistory(结束)gydF4y2Ba分钟从现在。gydF4y2Ba]);gydF4y2Ba打破gydF4y2Ba;gydF4y2Ba结束gydF4y2Ba结束gydF4y2Ba

估计在80分钟达到故障阈值。gydF4y2Ba

检查最近的时间序列模型。gydF4y2Ba

sysgydF4y2Ba

sys =离散时间状态空间模型发现:x (t + Ts) = x (t) + K e (t) y (t) = C x (t) + e (t) = (x1, x2) x1 0 1 x2 0.2624 - 0.746 C = (x1, x2)日元1 0 K = 0.3002日元x1 0.3902 x2样品时间:600秒参数化:规范形式指数:2。干扰组件:估计很多免费的系数:4使用“idssdata”、“getpvec”、“getcov”参数及其不确定性。对时域数据状态:估计使用党卫军“tsFeature”。适合估算数据:92.53%(预测聚焦)消防工程:499.3,MSE: 442.7gydF4y2Ba

拟合优度增加到90%以上,这一趋势是正确地捕获。gydF4y2Ba

结论gydF4y2Ba

这个例子展示了如何从测量数据中提取特征进行状态监测和预测。基于提取的特征,动态生成模型,验证和用于预测故障的时间,这样可以实际故障发生之前采取行动。gydF4y2Ba

利用振动信号状态监测和预测gydF4y2Ba

数据描述gydF4y2Ba

特征提取gydF4y2Ba

状态监测和预测gydF4y2Ba

结论gydF4y2Ba

相关的话题gydF4y2Ba