风力发电机高速轴承预测

打开生活的脚本

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

数据集

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

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

复制文件(...完整文件（matlabroot，“工具箱”，“predmaint”，...“predmaintdemos”，“windTurbineHighSpeedBearingPrognosis”），...“WindTurbineHighSpeedBearingPrognosis-Data-master”)fileattrib（完整文件(“WindTurbineHighSpeedBearingPrognosis-Data-master”，‘* .mat‘），“+w”）

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

时间单位=“\次5天”；

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

如果存在(“WindTurbineHighSpeedBearingPrognosis-Data-master.zip”，“文件”)解压(“WindTurbineHighSpeedBearingPrognosis-Data-master.zip”) timeUnit =“一天”；结束

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

数据导入

创建一个文件集成数据存储风力发电机的数据。数据包括一个振动信号和一个转速表信号。的文件集成数据存储将解析文件名并提取日期信息，如下所示IndependentVariables。有关详细信息，请参阅与此示例关联的支持文件中的帮助器函数。金宝app

hsbearing = fileEnsembleDatastore (...fullfile (“。”，“WindTurbineHighSpeedBearingPrognosis-Data-master”），...“.mat”); hsbearing.DataVariables=[“振动”，“转速表”]; hs.独立变量=“日期”; hsbearing.SelectedVariables=[“日期”，“振动”，“转速表”];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 double] 12-Mar-2013 06:17:10 [585936×1 double] [2447×1 double] 13-Mar-2013 06:34:04 [585936×1 double] [2438×1 double] 14-Mar-2013 06:50:41 [585936×1 double] [2390×1 double]::::::

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

fs=97656；%赫兹

数据探索

本节探讨时域和频域数据，并寻求为预测目的提取哪些特征的灵感。

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

复位（hs轴承）t开始=0；数字保持在…上而Hasdata (hsbearing) data = read(hsbearing);v = data.vibration {1};T = tstart + (1:length(v))/fs;向下采样信号以减少内存使用绘图（t（1:10:end），v（1:10:end））t开始=t（end）；结束持有从xlabel('时间，每天6秒，总共50天')伊拉贝尔(‘加速度（g）’）

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

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

hsbearing。DataVariables = [“振动”，“转速表”，“SpectralKurtosis”];颜色= parula (50);图保存在…上重置(hsbearing) day = 1;而Hasdata (hsbearing) data = read(hsbearing);data2add =表;获得振动信号和测量数据v=数据振动{1}；%计算谱峰度，窗宽= 128wc = 128;[SK, F] = pkurtosis(v, fs, wc);data2add。谱峭度= {table(F, SK)};%绘制光谱峰度图地块3（F，日*个（尺寸（F）），SK，“颜色”，颜色（天，：）%写入光谱峰度值writeToLastMemberRead (hsbearing data2add);%增加天数Day = Day + 1;结束持有从xlabel(的频率(赫兹))伊拉贝尔(的时间(天)) zlabel (“光谱峰度”)网格在…上视图（-45,30）cbar=彩色条；ylabel（cbar，“故障级别(0 -正常，1 -故障)”）

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

特征提取

基于上一节中的分析，将从时域信号和频谱峰度中提取一系列统计特征。有关这些特征的更多数学细节，请参见[2-3]。

首先，在将功能名称写入fileEnsembleDatastore之前，在DataVariables中预先指定这些功能名称。

hsbearing.DataVariables=[hsbearing.DataVariables；...“平均值”；“性病”；“偏斜”；“峰度”；“Peak2Peak”；...“RMS”；“CrestFactor”；“ShapeFactor”；“脉冲发生器”；“MarginFactor”；“能源”；...“我的意思是”；“SKStd”；“倾斜度”；“峰度”];

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

hsbearing。选择edVariables = [“振动”，“SpectralKurtosis”]；重置（hs轴承）而hasdata（hsbearing）data=read（hsbearing）；v=data.vibration{1}；SK=data.SpectralKurtosis{1}.SK；features=table；%时域特征features.Mean=Mean（v）；features.Std=Std（v）；features.skutosis=skutosis（v）；features.Peak2Peak=Peak2Peak（v）；features.RMS=RMS（v）；features.CrestFactor=max（v）/features.RMS.ShapeFactor=features.RMS/Mean（abs（v））；features.pulsefactor=max（v）/Mean（abs（v））^2；特征。能量=总和（v.^2）；%光谱峰度相关特征特性。SKMean =意味着(SK);特性。SKStd =性病(SK);特性。SKSkewness =偏态(SK);特性。SKKurtosis =峰度(SK);%将派生的特性写入相应的文件WriteLastMemberRead（hs轴承、功能）；结束

选择自变量日期并将提取的所有特征构建特征表。

hsbearing。选择edVariables = [“日期”，“平均值”，“性病”，“偏斜”，“峰度”，“Peak2Peak”，...“RMS”，“CrestFactor”，“ShapeFactor”，“脉冲发生器”，“MarginFactor”，“能源”，...“我的意思是”，“SKStd”，“倾斜度”，“峰度”];

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

featureTable =收集(高(hsbearing));

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

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

featureTable = table2timetable (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从xlabel(“时间”)伊拉贝尔(“特征值”)传奇(在平滑的，“平滑后”)头衔(“SKMean”）

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

训练数据

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

breaktime=datetime（2013,3,27）；breakpoint=find（featureTableSmooth.Date“最后一次”）;trainData = featuretables平滑(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}$ ．

%由于已完成移动窗口平滑，请将“窗口大小”设置为0以%关闭函数中的平滑特性重要性=单调性（列车数据，“窗口大小”, 0);helperSortedBarPlot (featureImportance的单调性）;

信号的峰度是信号单调性的最主要特征。

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

trainDataSelected=trainData（：，Feature重要性{：，：}>0.3）；featureSelected=featureTableSmooth（：，Feature重要性{：，：}>0.3）

精选特色=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);

平均值、标准偏差和PCA系数用于处理整个数据集。

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

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

图numData = size(featureTable, 1);PCA1, PCA2， []， 1:numData，“填充”)xlabel(“PCA 1”)伊拉贝尔(“PCA 2”)cbar=颜色条；ylabel（cbar[“时间”时间单位')'])

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

healthIndicator=PCA1；

可视化健康指标。

图表（特征选择日期、健康指标、，“o”)xlabel(“时间”)头衔(“健康指标”）

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

指数退化模型定义为

$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）；

阈值的选择通常基于机器的历史记录或某些特定领域的知识。由于此数据集中没有可用的历史数据，因此选择健康指标的最后一个值作为阈值。建议根据平滑（历史）数据选择阈值，以便部分缓解平滑的延迟效应。

阈值=健康指标（结束）；

如果有历史数据，请使用适合提供的方法exponentialDegradationModel估计先验和截距。然而，历史数据不能用于该风力涡轮机轴承数据集。斜率参数的先验是任意选取的，方差很大( $E （ θ ）＝ 1 ，变量（ θ ）＝ 10^{6} ， E （ β ）＝ 1 ，变量（ β ）＝ 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 = exponentialDegradationModel (...“θ”1....“先锋”，1e6，...“贝塔”1....“BetaVariance”，1e6，...“Phi”, -1,...“噪音回避”，(0.1*阈值/(阈值+ 1))^2，...“SlopeDetectionLevel”, 0.05);

使用预测和使现代化方法预测RUL并实时更新参数分布。

在每次迭代中保持记录totalDay = length(healthIndicator) - 1;estRULs = 0 (totalDay, 1);trueRULs = 0 (totalDay, 1);CIRULs = zeros(totalDay, 2);pdfRULs = cell(totalDay, 1); / /累计天数%创建用于打印更新的图形和轴图ax1 = subplot(2,1,1);Ax2 = subplot(2, 1, 2);为currentDay = 1: totalDay%更新模型参数后验分布更新(mdl [currentDay healthIndicator (currentDay)))%预测剩余使用寿命[estRUL，CIRUL，pdfRUL]=预测值（mdl，...[currentDay healthIndicator (currentDay)),...阈值）；trueRUL=totalDay-currentDay+1；%更新RUL分布图HelperPlotRend（ax1、currentDay、healthIndicator、mdl、阈值、时间单位）；helperPlotRUL（ax2、trueRUL、estRUL、CIRUL、pdfRUL、时间单位）%保留预测结果estRULs（currentDay）=estRUL；特鲁尔（当前日）=特鲁尔；CIRULs（当前日期：）=CIRUL；pdfRULs{currentDay}=pdfRUL；%暂停0.1秒以使动画可见暂停（0.1）结束

这是实时RUL估计的动画。

性能分析

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

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

哪里 $r （ t ）$ 是当时的估计RUL $t$ ， $r^{＊} （ t ）$ 这是当时真正的规则吗 $t$ ， $Θ （ t ）$ 估计的模型参数是否在时刻 $t$ ．

alpha=0.2；检测时间=mdl.slopedetectionstant；prob=helperAlphaLambdaPlot（alpha，trueRULs，estRULs，CIRULs，...pdfRULs, detectTime, breakpoint, timeUnit);标题(“α- \ \λ阴谋”）

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

将预测RUL的概率可视化到 $α$ 跳跃

图t=1:totalDay；保持在…上Plot (t, prob) Plot ([breakpoint breakpoint]， [0 1]，“k -”。)持有从xlabel([“时间”时间单位')']) ylabel (“概率”)传奇('预测RUL在\alpha界限内的概率'，“Train-Test断点”)标题('概率在\alpha范围内，\alpha='1'num2str（α*100）'%'])