风力机高速轴承预测

打开实时脚本

这个例子展示了如何建立一个指数退化模型来实时预测风力涡轮机轴承的剩余使用寿命(RUL)。指数退化模型基于参数先验和最新测量(运行到故障的历史数据可以帮助估计模型参数先验，但它们不是必需的)来预测RUL。该模型能够实时检测显著的退化趋势，并在新的观测数据可用时更新其参数先验。该实例遵循典型的预测工作流程:数据导入与探索、特征提取与后处理、特征重要性排序与融合、模型拟合与预测、性能分析。

数据集

数据集采集自20齿小齿轮[1]驱动的2MW风力发电机高速轴。连续50天每天采集一个6秒的振动信号(3月17日有2次测量，在本例中作为两天处理)。在50天的时间内，一个内圈故障发展并导致轴承故障。

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

拷贝文件(.．.fullfile (matlabroot“工具箱”，“predmaint”，.．.“predmaintdemos”，“windTurbineHighSpeedBearingPrognosis”)，.．.“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”) timeUnit =“天”；结束

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

数据导入

创建一个fileEnsembleDatastore风力涡轮机的数据。数据包含振动信号和转速表信号。的fileEnsembleDatastore将解析文件名并提取日期信息为IndependentVariables．有关详细信息，请参阅与此示例关联的支持文件中的helper函数。金宝app

hsbearing =文件集合数据存储(.．.fullfile (“。”，“WindTurbineHighSpeedBearingPrognosis-Data-master”)，.．.“.mat”）;hsbearing。DataVariables = [“振动”，“一环”];hsbearing。IndependentVariables =“日期”；hsbearing。选择edVariables = [“日期”，“振动”，“一环”];hsbearing。ReadFcn = @helperReadData;hsbearing。WriteToMemberFcn = @helperWriteToHSBearing;高(hsbearing)

ans = M×3高表日期振动性心动过速  ____________________ _________________ _______________ 07 - 3月- 2013 01:57:46(585936×1双)(2446×1双)08 - mar - 2013 02:34:21(585936×1双)(2411×1双)09 - 3月- 2013 02:33:43(585936×1双)(2465×1双)10 - 3月- 2013年03:01:02(585936×1双)(2461×1双)11 - 3月- 2013 03:00:24(585936×1双)(2506×1双)12 - 3月- 2013年06:17:10(585936×1双)(2447×1双)13 - 3月- 2013 06:34:04(585936×1双)(2438×1双)14 - 3月- 201306:50:41 [585936×1 double] [2390×1 double]::::::

振动信号的采样率为97656 Hz。

Fs = 97656;%赫兹

数据探索

本节探讨了时域和频域的数据，并寻求灵感，以提取哪些特征用于预后。

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

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

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

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

hsbearing。DataVariables = [“振动”，“一环”，“SpectralKurtosis”];颜色= parula(50);图保存在复位(hsbearing)天= 1;而Hasdata (hsbearing) data = read(hsbearing);Data2add = table;获得振动信号及测量数据V = data.vibration{1};计算窗口大小为128的谱峰度Wc = 128;[SK, F] = pkurtosis(v, fs, wc);data2add。spectral峭度={表(F, SK)};绘制谱峰度plot3(F, day*ones(size(F))， SK，“颜色”，颜色(天，:))写入谱峰度值writeToLastMemberRead (hsbearing data2add);增加天数天=天+ 1;结束持有从包含(的频率(赫兹)) ylabel (的时间(天)) zlabel (“谱峰态”网格)在查看(- 45,30)cbar = colorbar;ylabel (cbar“故障级别(0 -正常，1 -故障)”）

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

特征提取

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

首先，在将DataVariables写入fileEnsembleDatastore之前，在DataVariables中预先分配特性名称。

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

计算每个集成成员的特性值。

hsbearing。选择edVariables = [“振动”，“SpectralKurtosis”];重置(hsbearing)而Hasdata (hsbearing) data = read(hsbearing);V = data.vibration{1};SK = data.SpectralKurtosis{1}.SK;Features = table;%时域特性特性。均值=均值(v);特性。Std = Std (v);特性。偏度=偏度(v);特性。峰度=峰度(v);特性。peak = 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 = std(SK);特性。SKSkewness =偏度(SK);特性。峰度=峰度(SK);%将派生的特性写入相应的文件writeToLastMemberRead (hsbearing、特点);结束

选择自变量日期并将所有提取出来的特征构造成特征表。

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

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

featureTable =收集(高(hsbearing));

使用并行池“本地”评估tall表达式:-通过1 / 1:在1秒内完成评估

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

featureTable = table2schedule (featureTable)

featureTable =50×15时间表日期的意思是性病偏态峰度Peak2Peak RMS CrestFactor ShapeFactor ImpulseFactor MarginFactor能源SKMean SKStd SKSkewness SKKurtosis  ____________________ _______ ______ ___________ ________ _________ ______ ___________ ___________ _____________ ____________ __________ __________ ________ __________ __________ 07 - 3月- 2013 01:57:46 0.34605 2.2705 0.0038699 2.9956 21.621 2.2967 4.9147 1.2535 6.1607 3.3625 3.0907 0.001253 e + 06 - mar - 2013 02:34:21 08年0.025674 -0.22882 3.362 0.24409 2.06210.0030103 3.0195 19.31 2.0765 4.9129 1.2531 2.5266e+06 0.0046823 0.020888 0.057651 3.3508 09-Mar-2013 02:33:43 0.21873 2.1036 -0.0010289 3.0224 21.474 2.1149 5.2143 1.2539 6.5384 3.50004 4.9953 10-Mar-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-Mar-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.0135590.032193 0.11355 3.848 12-Mar-2013 06:17:10 2.0791 -0.008428 3.018 20.05 2.0997 4.7966 1.2537 3.5907 3.5833 e+06 0.0017539 0.028326 -0.1288 3.8072 - mar -2013 06:34:04 0.21293 1.972 -0.0014294 3.0174 18.837 1.9834 4.8496 1.2539 3.8441 2.3051e+06 0.0039353 0.029292 -1.4734 8.1242 14-Mar-2013 06:50:41 0.24401 1.8114 0.0022161 3.0057 17.862 1.8278 4.8638 1.9575 4.1821 1.9575e+06 0.0013107 0.022468 0.27438 2.8597 15-Mar-2013 06:50:03 0.20984 1.9973 0.001559 3.0711 21.122.0083 5.4323 1.2568 6.8272 4.2724 2.3633e+06 0.0023769 0.047767 -2.5832 20.171 16-Mar-2013 06:56:43 0.2342 -0.0019594 3.0072 18.832 1.9979 5.0483 0.030418 0.52322 4.0082 17-Mar-2013 06:56:04 0.21657 2.113 - 0.0013727 0.030418 0.52322 4.0082 17-Mar-2013 06:56:04 0.21657 2.1247 21.858 2.1241 5.4857 1.2587 6.9048 4.0916 2.6437e+06 0.0084907 -0.012744 3.0934 21.589 2.1423 4.7574 1.2575 5.9825 2.6891e+06 0.019624 0.055537 3.1986 17.79618-Mar-2013 18:47:15 0.21919 2.1284 -0.0002039 3.1647 24.051 2.1396 5.7883 1.2595 4.2914 2.6824e+06 0.016315 0.064516 2.8735 11.632 20-Mar-2013 00:33:54 0.35865 2.2536 -0.002308 3.0817 22.633 2.2819 5.2771 1.2571 6.6339 3.0511e+06 0.0011097 0.051539 -0.056774 7.0712 21-Mar-2013 00:33:14 0.1908 2.1782 -0.00019286 3.1548 25.515 2.1865 6.0521 1.26 7.6258 4.3945 2.8013e+06 0.0040392 0.066254 -0.39587 12.111 22-Mar-2013 00:39:50 0.20569 2.1861 0.0020375 3.2691 26.439 2.1958 6.17531.2633 7.8011 4.4882 2.825e+06 0.020676 0.077728 2.6038 11.088

特征后处理

提取的特征通常与噪声有关。相反趋势的噪声有时会对RUL预测产生不利影响。此外，接下来要介绍的特征性能指标之一，单调性，对噪声的鲁棒性不强。因此，对提取的特征应用滞后窗口为5步的因果移动均值滤波器，其中“因果”意味着在移动均值滤波中没有使用未来值。

variableNames = featutable . properties . variableNames;featureTable = varfun(@(x) movmean(x， [5 0])， featureTable);featutablessmooth . properties . variableNames = variableNames;

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

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

移动均值平滑会给信号带来时间延迟，但在RUL预测中选择适当的阈值可以缓解延迟效应。

训练数据

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

Breaktime = datetime(2013, 3,27);breakpoint = find(featutable平滑。日期<休息时间，1，“最后一次”）;trainData = featuttable平滑(1:断点，:);

特征重要性排序

在本例中，使用[3]提出的单调性来量化特征的优点，用于预测。

单调性的 $我$ th特性 $x_{我}$ 计算为

$单调性（ x_{我} ）＝ \frac{1}{米} \sum_{j ＝ 1}^{米} \frac{| 数量的积极的 diff （ x_{我}^{j} ） - 数量的负 diff （ x_{我}^{j} ） |}{n - 1}$

在哪里 $n$ 在这种情况下，测量点的数量是多少 $n ＝ 50$ ． $米$ 在这种情况下，是否监视了机器的数量 $米＝ 1$ ． $x_{我}^{j}$ 是 $我$ 上测量的特征 $j$ 机。 $diff （ x_{我}^{j} ）＝ x_{我}^{j} （ t ） - x_{我}^{j} （ t - 1 ）$ ，即信号之差 $x_{我}^{j}$ ．

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

信号的峰度是基于单调性的顶部特征。

下一节选择特征重要度评分大于0.3的特征进行特征融合。

trainDataSelected = trainData(:， featureImportance{:，:}>0.3);featureSelected = featuttable平滑(:，featureImportance{:，:}>0.3)

featureSelected =50×5时间表Date Mean Kurtosis ShapeFactor MarginFactor SKStd ____________________ _______________ ___________ ____________________ 07-Mar-2013 01:57:46 0.34605 2.9956 1.3625 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.7722 0.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.2333.0272 1.2545 3.8282 0.027977 14-Mar-2013 06:50:41 0.23299 3.0249 1.2544 3.9047 0.02824 15-Mar-2013 06:50:03 0.2315 3.033 1.2548 3.9706 0.032417 16-Mar-2013 06:56:43 0.23475 3.0273 1.2546 3.9451 0.031744 17-Mar-2013 06:56:04 0.23498 3.0407 1.2551 3.9924 0.032691 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.047942 21-Mar-2013 00:33:14 0.23537 3.1044 1.2578 3.9862 0.051023 22-Mar-2013 00:39:50 0.23079 3.1481 1.2593 4.072 0.058908 ⋮

降维与特征融合

在这个例子中，主成分分析(PCA)用于降维和特征融合。在执行PCA之前，将特征归一化为相同的尺度是一个很好的做法。请注意，归一化中使用的PCA系数以及平均值和标准差是从训练数据中获得的，并应用于整个数据集。

mean(trainDataSelected{:，:});sdTrain = std(trainDataSelected{:，:});trainDataNormalized = (trainDataSelected{:，:} - meanTrain)./sdTrain;coef = pca(trainDataNormalized);

平均数，标准偏差和PCA系数被用来处理整个数据集。

PCA1 = (featureSelected{:，:} - meanTrain) ./ sdTrain * coef(:， 1);PCA2 = (featureSelected{:，:} - meanTrain) ./ sdTrain * coef(:， 2);

在前两个主成分的空间中可视化数据。

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

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

healthIndicator = PCA1;

可视化运行状况指示器。

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

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

指数退化模型定义为

$h （ t ）＝ ϕ + θ 经验值（ β t + ϵ - \frac{σ^{2}}{2} ）$

在哪里 $h （ t ）$ 是作为时间函数的运行状况指标。 $ϕ$ 截距项是一个常数。 $θ$ 而且 $β$ 是随机参数决定模型的斜率，在哪里 $θ$ 是对数正态分布的 $β$ 是系统。在每个时间步 $t$ 的分布。 $θ$ 而且 $β$ 是根据最新的观察更新到后验的 $h （ t ）$ ． $ϵ$ 高斯白噪声屈服于 $N （ 0 ， σ^{2} ）$ ．的 $- \frac{σ^{2}}{2}$ 指数中的一项是 $h （ t ）$ 满足

$E ［ h （ t ） | θ ， β ］＝ ϕ + θ 经验值（ β t ）$ ．

在这里，指数退化模型适合于上一节中提取的健康指标，并在下一节中评估性能。

首先移动运行状况指示器，使其从0开始。

healthIndicator = healthIndicator - healthIndicator(1);

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

threshold = healthIndicator(end);

如果有历史数据，请使用适合方法由exponentialDegradationModel估计先验和截距。然而，历史数据是不可用的这个风力涡轮机轴承数据集。斜率参数的先验是任意选择的，且方差较大( $E （ θ ）＝ 1 ， Var （ θ ）＝ 10^{6} ， E （ β ）＝ 1 ， Var （ β ）＝ 10^{6}$ )，因此该模型主要依赖于观测数据。基于 $E ［ h （ 0 ）］＝ ϕ + E （ θ ）$ ,拦截 $ϕ$ 设置为 $- 1$ 所以模型也会从0开始。

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

$Δ h （ t ） \approx （ h （ t ） - ϕ ） Δ ϵ （ t ）$

这里假设噪声的标准差在接近阈值时引起健康指标变化的10%。因此，噪声的标准差可以表示为 $\frac{10 ％ \cdot 阈值}{阈值 - ϕ}$ ．

指数退化模型还提供了评估斜率的重要性的功能。一旦检测到健康指标的显著斜率，模型将忘记之前的观察结果，并基于原始先验重新进行估计。检测算法的灵敏度可通过指定进行调整SlopeDetectionLevel．如果p值小于SlopeDetectionLevel，则表示斜率被检测到。在这里SlopeDetectionLevel设置为0.05。

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

mdl =指数退化模型(.．.“θ”, 1.．.“ThetaVariance”1 e6,.．.“β”, 1.．.“BetaVariance”1 e6,.．.“φ”, 1.．.“NoiseVariance”， (0.1*threshold/(threshold + 1))^2，.．.“SlopeDetectionLevel”, 0.05);

使用predictRUL而且更新方法预测RUL并实时更新参数分布。

在每次迭代中保持记录totalDay = length(healthIndicator) - 1;estRULs = 0 (totalDay, 1);trueRULs = 0 (totalDay, 1);CIRULs = 0 (totalDay, 2);pdfRULs = cell(totalDay, 1);为绘图更新创建图形和轴图ax1 = subplot(2,1,1);Ax2 = subplot(2,1,2);为currentDay = 1:totalDay更新模型参数后验分布update(mdl， [currentDay healthIndicator(currentDay)])预测剩余使用寿命[estRUL, CIRUL, pdfRUL] = predictRUL(mdl，.．.[currentDay healthIndicator (currentDay)),.．.阈值);trueRUL = totalDay - currentDay + 1;%正在更新RUL分布图helperplotrend (ax1, currentDay, healthIndicator, mdl, threshold, timeUnit);helpplotrul (ax2, trueRUL, estRUL, CIRUL, pdfRUL, timeUnit)%保留预测结果estRULs(currentDay) = estrull;trueRUL (currentDay) = trueRUL;CIRULs(currentDay，:) = CIRUL;pdfRUL {currentDay} = pdfRUL;暂停0.1秒使动画可见暂停(0.1)结束

这是实时RUL估计的动画。

性能分析

$α$ - $λ$ Plot用于预后性能分析[5]，其中 $α$ Bound设置为20%。估计的RUL在 $α$ true RUL的边界作为模型的性能度量来计算:

$公关（ r^{＊} （ t ） - α r^{＊} （ t ） < r （ t ） < r^{＊} （ t ） + α r^{＊} （ t ） | Θ （ t ））$

在哪里 $r （ t ）$ 当时的估计RUL是多少 $t$ ， $r^{＊} （ t ）$ 真正的规则是什么时候 $t$ ， $Θ （ t ）$ 估计的模型参数是否及时 $t$ ．

Alpha = 0.2;detectTime = mdl.SlopeDetectionInstant;prob = helperAlphaLambdaPlot(alpha, trueRULs, estRULs, CIRULs，.．.pdfRULs, detectTime, breakpoint, timeUnit);标题(“α- \ \λ阴谋”）

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

可视化预测RUL的概率 $α$ 绑定。

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