数据集

数据集采集自一台2MW风力机高速轴，由一个20齿小齿轮[1]驱动。连续50天每天采集6秒的振动信号(3月17日有两次测量，在本例中被视为两天)。在50天的时间内，出现了内圈故障，导致轴承失效。

工具箱中提供了数据集的精简版本。若要使用精简数据集，请将数据集复制到当前文件夹并启用其写权限。

拷贝文件(．..fullfile (matlabroot“工具箱”,“维护前”,．..“predmaintdemos”,“风机高速轴承预测”），．..“WindTurbineHighSpeedBearingPrognosis-Data-master”) fileattrib (fullfile (“WindTurbineHighSpeedBearingPrognosis-Data-master”,‘* .mat‘），' + w ')

压缩数据集的测量时间步长为5天。

timeUnit =“\ * 5天”;

有关完整数据集，请转到此链接https://github.com/mathworks/WindTurbineHighSpeedBearingPrognosis-Data，将整个存储库下载为zip文件，并将其保存在与此活动脚本相同的目录中。使用这个命令解压文件。完整数据集的测量时间步长为1天。

如果存在(“WindTurbineHighSpeedBearingPrognosis-Data-master.zip”,“文件”)解压缩(“WindTurbineHighSpeedBearingPrognosis-Data-master.zip”)时间单位=“天”;结束

此示例中的结果是从完整数据集生成的。强烈建议下载完整数据集以运行此示例。从压缩数据集生成的结果可能没有意义。

数据导入

创建一个fileEnsembleDatastore风力涡轮机数据的完整性。数据包含振动信号和转速表信号。这个fileEnsembleDatastore将解析文件名并提取日期信息为独立变量．有关更多细节，请参阅与本示例相关的支持文件中的helper函数。金宝app

hsbearing=文件集成数据存储(．..完整文件('.',“WindTurbineHighSpeedBearingPrognosis-Data-master”），．..“.mat”）;hsbearing。DataVariables = [“振动”,“一环”];hsbearing.independentvariables =.“日期”;hsbearing.selectedvariables = [“日期”,“振动”,“一环”];hsbearing。ReadFcn = @helperReadData;hsbearing。WriteToMemberFcn = @helperWriteToHSBearing;高(hsbearing)

这是一个M X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X 2013年3月10日03:01:02[585936×1双倍][2461×1双倍]2013年3月11日03:00:24[585936×1双倍][2506×1双倍]2013年3月12日06:17:10[585936×1双倍][2447×1双倍]2013年3月13日06:34:04[585936×1双倍][2438×1双倍]2013年3月14日06:50:41[585936×1双倍][2390×1双倍]：

振动信号采样率为97656 Hz。

fs = 97656;%赫兹

数据探索

本节探讨时域和频域的数据，并寻找用于预测目的的特征提取的灵感。

首先在时域可视化振动信号。在这个数据集中，连续50天有50个6秒的振动信号。现在把50个振动信号一个接一个地画出来。

重置(hsbearing) tstart = 0;图保存在而Hasdata (hsbearing) data = read(hsbearing);v = data.vibration {1};T = tstart + (1:length(v))/fs;%降低信号采样以减少内存使用绘图（T（1：10：结束），v（1：10：结束））tstart = t（结束）;结束持有从包含(“时间(秒)，每天6秒，共50天”) ylabel (“加速度(g)”)

振动信号在时域上呈现出信号冲击性增大的趋势。量化信号脉冲性的指标，如峰度、峰峰值、峰值因子等。，是该风力涡轮机轴承数据集[2]的潜在预测特征。

另一方面，频谱峰度被认为是在频域预测风力涡轮机的有力工具[3]。为了可视化频谱峰度随时间的变化，绘制频谱峰度值作为频率和测量日的函数。

hsbearing。DataVariables = [“振动”,“一环”,“幽灵荨麻疹”];颜色= parula (50);图保存在重置(hsbearing) day = 1;而Hasdata (hsbearing) data = read(hsbearing);data2add =表;获得振动信号和测量数据v = data.vibration {1};%计算窗口大小为128的光谱峰度wc = 128;[SK, F] = pkurtosis(v, fs, wc);data2add。谱峭度= {table(F, SK)};绘制光谱峰度plot3 (F,天*的(大小(F)), SK,“颜色”颜色(天,:))%写谱峰度值writeToLastMemberRead (hsbearing data2add);%增加天数天=天+1；结束持有从包含(的频率(赫兹)) ylabel (‘时间（天）’) zlabel (“谱峰态”网格)在视图(- 45,30)cbar = colorbar;ylabel (cbar'故障严重性（0-正常，1-故障）')

颜色条中显示的故障严重程度是归一化为0到1刻度的测量日期。观察到，随着机器条件的退化，光谱峰度值在10 kHz左右逐渐增加。谱峰度的统计特征，如平均值、标准差等。，将是轴承退化[3]的潜在指标。

特征提取

基于上一节的分析，我们将从时域信号和谱峰度中提取一组统计特征。[2-3]提供了更多关于这些特征的数学细节。

首先，在将特性名称写入文件ensembledatastore之前，在datavvariables中预先分配它们。

hsbearing。DataVariables = [hsbearing.DataVariables;．..“的意思是”;“性病”;“偏斜”;“峰度”;“Peak2Peak”;．..“RMS”;“CrestFactor”;“ShapeFactor”;“ImpulseFactor”;“MarginFactor”;“能量”;．..“SKMean”;“SKStd”;“SKSkewness”;“SKKurtosis”];

计算每个集合成员的特征值。

hsbearing.selectedvariables = [“振动”,“幽灵荨麻疹”];重置(hsbearing)而Hasdata (hsbearing) data = read(hsbearing);v = data.vibration {1};SK = data.SpectralKurtosis {1} .SK;特点=表;％时域功能特性。的意思= (v);特性。Std =性病(v);特性。偏态=偏斜度(v);特性。峰度=峰度(v);特性。Peak2Peak = Peak2Peak (v); features.RMS = rms(v); features.CrestFactor = max(v)/features.RMS; features.ShapeFactor = features.RMS/mean(abs(v)); features.ImpulseFactor = max(v)/mean(abs(v)); features.MarginFactor = max(v)/mean(abs(v))^2; features.Energy = sum(v.^2);光谱峰度相关特征特性。SKMean =意味着(SK);特性。SKStd =性病(SK);特性。SKSkewness =偏态(SK);特性。SKKurtosis =峰度(SK);%将衍生要素写入相应的文件writeToLastMemberRead (hsbearing、特点);结束

选择自变量日期以及构建要素表的所有提取功能。

hsbearing.selectedvariables = [“日期”,“的意思是”,“性病”,“偏斜”,“峰度”,“Peak2Peak”,．..“RMS”,“CrestFactor”,“ShapeFactor”,“ImpulseFactor”,“MarginFactor”,“能量”,．..“SKMean”,“SKStd”,“SKSkewness”,“SKKurtosis”];

由于特征表足够小，足以容纳内存(50 * 15)，所以在处理之前收集该表。对于大数据，建议以高格式执行操作，直到您确信输出足够小，可以装入内存。

featureTable=聚集（高（hs方向））；

通过1 of 1: Completed in 1 sec在1秒内完成

将表转换为时间表，以便时间信息始终与特性值相关联。

featureTable = table2timetable (featureTable)

featureTable =50×15时间表日期的意思是性病偏态峰度Peak2Peak RMS CrestFactor ShapeFactor ImpulseFactor MarginFactor能源SKMean SKStd SKSkewness SKKurtosis  ____________________ _______ ______ ___________ ________ _________ ______ ___________ ___________ _____________ ____________ __________ __________ ________ __________ __________ 07 - 3月- 2013 01:57:46 0.346052.2705 0.0038699 2.9956 21.621 2.2967 4.2535 6.1607 3.3625 3.0907e+06 0.001253 0.025674 -0.22882 3.362 08- 3 -2013 02:34:21 0.24409 2.0621 0.0030103 3.0195 19.31 2.0765 4.9129 1.2545 6.163 3.7231 2.5266e+06 0.0046823 0.020888 0.057651 3.3508 09- 3 -2013 02:33:43 0.21873 2.1036 -0.0010289 21.474 2.1149 5.2143 6.5384 3.87662.6208e+06 -0.0011084 0.022705 -0.50004 4.9953 10- march -2013 03:01:02 0.21372 2.0081 0.001477 3.0415 19.52 2.0194 5.286 1.2556 6.637 4.1266 2.3894e+06 0.0087035 0.034456 1.4705 8.1235 11- march -2013 03:00:24 0.21518 2.0606 0.0010116 3.0445 21.217 2.0718 5.0058 1.2554 6.2841 3.8078 2.515e+06 0.013559 0.032193 0.11355 3.848 12- march -2013 06:17:10 0.293352.0791 -0.008428 3.018 20.05 2.0997 4.7966 1.2537 6.0136 3.5907 2.5833e+06 0.0017539 0.028326 -0.1288 3.8072 13.3 -2013 06:34:04 0.21293 1.972 -0.00142941.9575e+06 0.0013107 0.022468 0.27438 2.8597 15- march -2013 06:50:03 0.20984 1.9973 0.001559 3.0711 21.12 2.0083 5.4323 1.2568 6.8272 4.2724 2.3633e+06 0.0023769 0.047767 -2.5832 20.171 16- march -2013 06:56:43 0.23318 1.9842 -0.0019594 3.0072 18.832 1.9979 5.0483 1.254 6.3304 3.9734 2.3387e+06 0.0076327 0.030418 0.52322 4.00820.21657 2.113 -0.0013711 3.1247 21.858 2.1241 5.4857 1.2587 6.9048 4.0916 2.6437e+06 0.0084907 0.037876 2.3753 11.491 17- 3 -2013 18:47:56 0.19381 2.1335 -0.012744 3.0934 21.589 2.1423 4.7574 1.2575 5.9825 3.5117 2.6891e+06 0.019624 0.055537 3.1986 17.796 18- 3 -2013 18:47:15 0.21919 2.1284 -0.0002039 3.1647 24.051 2.1396 1.2595 7.29024.2914 2.6824e+06 0.016315 0.064516 2.8735 11.632 20- march -2013 00:33:54 0.35865 2.2536 -0.002308 3.0817 22.633 2.2819 5.2771 1.2571 6.6339 3.6546 3.0511e+06 0.0011097 0.051539 -0.056774 7.0712 21- march -2013 00:33:14 0.1908 2.1782 -0.00019286 3.1548 25.515 2.1865 1.26 7.6258 4.3945 2.8013e+06 0.0040392 0.066254 -0.39587 12.11100:39:50 0.20569 2.1861 0.0020375 3.2691 26.439 2.1958 6.1753 1.2633 7.8011 4.4882 2.825e+06 0.020676 0.077728 2.6038 11.088⋮

特征后处理

提取的特征通常与噪声有关。相反趋势的噪声有时会对RUL预测造成危害。此外，接下来要介绍的性能指标之一单调性对噪声不具有鲁棒性。因此，对提取的特征采用滞后窗为5步的因果移动平均滤波器，其中“因果”表示在移动平均滤波中不使用未来值。

variableNames=featureTable.Properties.variableNames；featureTableSmooth=varfun（@（x）movmean（x，[5 0]），featureTable）；featureTableSmooth.Properties.variableNames=variableNames；

下面是一个显示平滑前后特征的例子。

图保存在情节(featureTable。日期,featureTable.SKMean) plot(featureTableSmooth.Date, featureTableSmooth.SKMean) hold从包含(“时间”) ylabel (“特征值”)传说(“平滑前”,平滑后的)标题(“SKMean”)

移动平均平滑引入了信号的时延，但在RUL预测中选择合适的阈值可以缓解时延效应。

训练数据

在实践中，在开发预测算法时，无法获得整个生命周期的数据，但可以合理地假设已经收集了生命周期早期的一些数据。因此，前20天(生命周期的40%)收集的数据被视为训练数据。下面的特征重要性排序和融合仅基于训练数据。

Breaktime = datetime(2013, 3, 27);断点=找到(featureTableSmooth。日期< breaktime, 1，“最后”）;trainData = featuretables平滑(1:断点，:);

特征重要性排序

在这个例子中，[3]提出的单调性被用来量化特征对预后的价值。

单调性 $\mathit{我}$th特性 ${\mathit{x}}_{\mathit{我}}$被计算为

在哪里 $\mathit{N}$在这种情况下，是测量点的数量吗 $\mathit{N}=50$． $\mathit{M}$在这种情况下，是否监控机器的数量 $\mathit{M}=1.$． ${\mathit{x}}_{\mathit{我}}^{\mathit{J}}$是 $\mathit{我}$测量的特征 $\mathit{J}$机。 $\mathrm{差异}\left({\mathit{x}}_{\mathit{我}}^{\mathit{J}}\right)={\mathit{x}}_{\mathit{我}}^{\mathit{J}}\left(\mathit{T}\right)-{\mathit{x}}_{\mathit{我}}^{\mathit{J}}\left(\mathit{T}-1.\right)$，即信号的差值 ${\mathit{x}}_{\mathit{我}}^{\mathit{J}}$．

由于移动窗口平滑已经完成，设置'WindowSize'为0%关闭函数内的平滑featureImportance =单调性(trainData,“WindowSize”, 0);helperSortedBarPlot (featureImportance“单调性”）;

信号的峰度是基于单调性的最高特征。

在下一节中，特征重要性评分大于0.3的特征被选择进行特征融合。

trainDataSelected = trainData(:， featureImportance{:，:}>0.3);featureSelected = featuretablessmooth (:， featureImportance{:，:}>0.3)

featureSelected =50×5时间表Date Mean Kurtosis ShapeFactor MarginFactor SKStd ____________________ _______________ ___________ ____________________ 07-Mar-2013 01:57:46 0.34605 2.9956 1.2535 3.3625 0.025674 08-Mar-2013 02:34:21 0.29507 3.0075 1.254 3.5428 0.023281 09-Mar-2013 02:33:43 0.26962 3.0125 1.254 3.6541 0.023089 10-Mar-2013 03:01:02 0.25565 3.0197 1.2544 3.77220.025931 11-Mar-2013 03:00:24 0.24756 3.0247 1.2546 3.7793 0.027183 12-Mar-2013 06:17:10 0.25519 3.0236 1.2544 3.7479 0.027374 13-Mar-2013 06:34:04 0.233 3.0272 1.2545 3.82820.031744 17-Mar-2013 18:47:56 0.21839 3.0533 1.2557 3.9792 0.037226 18-Mar-2013 18:47:15 0.21943 3.0778 1.2567 4.0538 0.043097 20-Mar-2013 00:33:54 0.23854 3.0905 1.2573 3.9658 0.0479424.072 - 0.058908⋮

降维与特征融合

本例中采用主成分分析(PCA)进行降维和特征融合。在进行PCA之前，将特征归一化到相同的尺度是一个很好的做法。注意，归一化使用的PCA系数和均值、标准差都是从训练数据中得到的，并应用于整个数据集。

meanTrain =意味着(trainDataSelected {:,:});sdTrain =性病(trainDataSelected {:,:});trainDataNormalized = (trainDataSelected{:，:} - meanTrain)./sdTrain;系数= pca (trainDataNormalized);

采用均值、标准差和主成分分析系数对整个数据集进行处理。

PCA1=（特征选择{：，：}-meanTrain）。/sdTrain*coef（：，1）；PCA2=（特征选择{：，：}-meanTrain）。/sdTrain*coef（：，2）；

在前两个主组件的空间中可视化数据。

图numData = size(featureTable, 1);PCA1, PCA2， []， 1:numData，“填充”)包含(“PCA 1”) ylabel (《PCA 2》) cbar = colorbar;ylabel (cbar [“时间”timeUnit“)”]）

该图表明，随着机器接近故障，第一个主成分正在增加。因此，第一主成分是一个有前途的融合健康指标。

healthIndicator = PCA1;

可视化运行状况指示器。

图绘制(featureSelected。日期、healthIndicator“o”)包含(“时间”)标题(“健康指示器”)

剩余使用寿命(RUL)估计拟合指数退化模型

指数退化模型定义为

在哪里 $\mathit{H}\left(\mathit{T}\right)$为时间函数的运行状况指示器。 $\mathrm{\phi .}$为作为常数的截距项。 $\theta$和 $\beta$是确定模型的斜率的随机参数，其中 $\theta$lognormal-distributed和 $\beta$是系统。在每个时间步 $\mathit{T}$，以及 $\theta$和 $\beta$根据最新观察到的 $\mathit{H}\left(\mathit{T}\right)$． $\mathrm{\epsilon .}$高斯白噪声是否屈服于 $\mathit{N}\left(0,{\mathrm{\sigma .}}^{2.}\right)$．这 $-\frac{{\mathrm{\sigma .}}^{2.}}{2.}$指数项中的项是 $\mathit{H}\left(\mathit{T}\right)$满足

$\mathit{E}\left[\mathit{H}\left(\mathit{T}\right)|\theta ,\beta \right]=\mathrm{\phi .}+\theta \mathrm{经验值}\left(\beta \mathit{T}\right)$．

这里，指数衰减模型适合于上一节中提取的运行状况指示器，在下一节中评估性能。

首先转移健康指标，使其从0开始。

healthIndicator = healthIndicator - healthIndicator(1);

阈值的选择通常基于机器的历史记录或一些特定领域的知识。由于此数据集中没有可用的历史数据，因此将选择运行状况指示器的最后一个值作为阈值。建议根据平滑(历史)数据选择阈值，以部分缓解平滑的延迟效应。

阈值= healthIndicator(结束);

如果历史数据可用，请使用适合提供的方法指数退化模型估计先验值和截距。但是，此风力涡轮机轴承数据集没有历史数据。斜率参数的先验值是任意选择的，方差很大( $\mathit{E}\left(\theta \right)=1.,\mathrm{Var}\left(\theta \right)={10}^{6.},\mathit{E}\left(\beta \right)=1.,\mathrm{Var}\left(\beta \right)={10}^{6.}$)，因此该模型主要依赖于观测数据。基于 $\mathit{E}\left[\mathit{H}\left(0\right)\right]=\mathrm{\phi .}+\mathit{E}\left(\theta \right)$,拦截 $\mathrm{\phi .}$被设置为 $-1.$这样模型也会从0开始。

健康指标的变化与噪声变化之间的关系可以推导为

这里假设噪声的标准偏差在接近阈值时导致健康指标变化的10%。因此噪声的标准差可以表示为 $\frac{10\%\cdot \mathrm{阈值}}{\mathrm{阈值}-\mathrm{\phi .}}$．

指数退化模型还提供了评估斜率显著性的功能。一旦检测到健康指标的显著斜率，模型将忘记以前的观察结果，并根据原始先验重新开始估计。可通过指定SlopeDetectionLevel．如果p值小于SlopeDetectionLevel，则宣布探测到坡度。在这里SlopeDetectionLevel为0.05。

现在用上面讨论的参数创建一个指数退化模型。

mdl=指数退化模型(．..“西塔”, 1．..“ThetaVariance”1 e6,．..“β”, 1．..“BetaVariance”1 e6,．..'phi', 1．..“NoiseVariance”，(0.1*阈值/(阈值+ 1))^2，．..“SlopeDetectionLevel”, 0.05);

使用refiCtrul.和更新提出了一种实时预测RUL和更新参数分布的方法。

在每次迭代中保持记录totalDay=长度（healthIndicator）-1；estRULs=零（总天数，1）；trueRULs=零（totalDay，1）；CIRULs=零（totalDay，2）；pdfRULs=单元格（总天数，1）；%创建图形和坐标轴用于绘图更新图ax1 = subplot(2,1,1);Ax2 = subplot(2, 1, 2);为currentDay=1:totalDay%更新模型参数后验分布更新(mdl [currentDay healthIndicator (currentDay)))%预测剩余的使用寿命[estRUL, CIRUL, pdfRUL] = prestrul (mdl, pdl)．..[currentDay healthIndicator (currentDay)),．..阈值);true = totalDay - currentDay + 1;%更新RUL分布图helperplotrend (ax1, currentDay, healthIndicator, mdl, threshold, timeUnit);helperPlotRUL(ax2, trueRUL, estRUL, CIRUL, pdfRUL, timeUnit)保持预测结果发感（截止日期）= estrul;Trueruls（截止日）= Truerul;线圈（截止日期:) = cirul;pdfruls {curstday} = pdfrul;暂停0.1秒使动画可见暂停(0.1)结束

这是实时RUL估计的动画。

性能分析

$\alpha$- $\lambda$Plot用于预测性能分析[5]，其中 $\alpha$Bound设置为20%。估计的RUL在 $\alpha$真实RUL的界限计算为模型的性能指标：

在哪里 $\mathit{R}\left(\mathit{T}\right)$是估计的RUL时间 $\mathit{T}$, ${\mathit{R}}^{*}\left(\mathit{T}\right)$是真正的规则吗 $\mathit{T}$, $\Theta \left(\mathit{T}\right)$估计的模型参数是否在时刻 $\mathit{T}$．

α= 0.2;detectTime = mdl.SlopeDetectionInstant;prob = helperAlphaLambdaPlot(alpha, trueRULs, estRULs, CIRULs，．..pdfRULs, detectTime, breakpoint, timeUnit);标题(“\alpha-\lambda绘图”)

由于预先设定的先验不能反映真实的先验，模型通常需要几个时间步长来调整到一个合适的参数分布。随着数据点的增多，预测也会变得更加准确。

将预测的RUL的概率可视化 $\alpha$边界。

图t = 1:totalDay;持有在Plot (t, prob) Plot ([breakpoint breakpoint]， [0 1]，“k -”。)举行从包含([“时间”timeUnit“)”]) ylabel (“概率”)传说(“在\alpha范围内预测RUL的概率”,“列车测试断点”)标题('概率在\ α范围内，\ α = 'num2str(α* 100)“%”]）

参考

[1] Bechhoefer, Eric, Brandon Van Hecke, David He。“改进光谱分析的处理。”预测和健康管理协会年会，新奥尔良，洛杉矶，10月．2013.

[2] Ali, Jaouher Ben等。基于无监督机器学习的风电轴承在真实实验条件下逐步退化的在线自动诊断应用声学132（2018）：167-181。

[3] 通过频谱峰度衍生指数和SVR预测风力涡轮机高速轴轴承健康状况应用声学120(2017): 1 - 8。

我是杰米·巴利斯。合并数据源来预测剩余的有用寿命——一种识别预后参数的自动化方法。”（2010）.

[5] Saxena，Abhinav等，“离线评估预后表现的指标。”国际预测与健康管理杂志1.1(2010): 4-23。