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

本示例演示了如何从滚珠轴承的振动信号中提取特征、进行健康监测和执行预测。本示例使用了信号处理工具箱™和系统识别工具箱™的功能，不需要预测性维护工具箱™。gydF4y2Ba

数据描述gydF4y2Ba

负载振动数据存储在gydF4y2BapdmBearingConditionMonitoringData.matgydF4y2Ba（这是一个从MathWorks支持文件站点下载的大型数据集〜88MB）。金宝app数据存储在单元阵列中，该单元阵列使用滚珠轴承信号模拟器在轴承外圈上具有单点缺陷的单点缺陷。它包含用于在不同健康状况（缺陷深度从3mm上升到3mm的缺陷深度增加）的轴承的多个振动信号。每个段存储以20kHz的采样率收集1秒的信号。在pdmbearingcondition monitoring.gydF4y2BaData.matgydF4y2Ba,gydF4y2BadefectDepthVecgydF4y2Ba存储缺陷深度随时间的变化情况gydF4y2BaexpTimegydF4y2Ba以分钟为单位存储相应的时间。gydF4y2Ba

url =gydF4y2Ba'//www.tatmou.com/金宝appsupportfiles/predmaint/condition-onitoring-andostics-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 Hz左右。现在计算故障轴承信号的平均峰值频率。平均峰值频率移至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％初始化向量以存储提取的平均峰值频率。gydF4y2BaMakePeakfreq =零（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“平均峰值频率(Hz)”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的gydF4y2BaExpansionConditionMonitoringStatistics.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

根据此图，设置平均峰值频率的阈值为2000Hz，以区分正常轴承和故障轴承，并最大限度地使用轴承。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'ts'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。使用“idssdata”、“getpvec”、“getcov”等参数及其不确定性。状态:使用SSEST对时域数据“tsFeature”进行估计。拟合估计数据:0.2763%(预测焦点)FPE: 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“平均峰值频率(Hz)”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] = forecast(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“平均峰值频率(Hz)”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（任何（均值划线（Tcur-Batchsize + 1：Tcur）>阈值））Disp（gydF4y2Ba“受监控变量超过故障阈值”gydF4y2Ba）;gydF4y2Ba打破gydF4y2Ba；gydF4y2BaelseifgydF4y2Ba(任何(yF.y >阈值)gydF4y2Ba%估计系统到达故障阈值的时间。gydF4y2BaTalarm = Tforecast（finc（yf.y>阈值，1））;DISP（[gydF4y2Ba“估计会达到失败的门槛”gydF4y2Banum2str（Talarm-thistory（END））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）a = x1 x2 x1 01 x2 0.2624 0.746 C = X1 X2 Y1 1 0 K = Y1 X1 0.3902 X2 0.3002采样时间：600秒参数化：指数规范形式：2.干扰组件：估计自由系数的数量：4使用“idsdata”，“getpvec”，“getcov”参数及其不确定性。状态:使用SSEST对时域数据“tsFeature”进行估计。适合估算数据：92.53％（预测焦点）FPE：499.3，MSE：442.7gydF4y2Ba

拟合优度提高到90%以上，正确地捕捉了趋势。gydF4y2Ba

结论gydF4y2Ba

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

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

数据描述gydF4y2Ba

特征提取gydF4y2Ba

状态监测和预后gydF4y2Ba

结论gydF4y2Ba

相关的话题gydF4y2Ba

预测性维护工具箱文档gydF4y2Ba

金宝app

用Matlab预测维护gydF4y2Ba