一种小样本装备故障在线预测方法

1.本发明涉及装备故障预测技术领域，具体为一种小样本装备故障在线预测方法。


背景技术：

2.故障预测技术能够立足于故障机理分析，假定故障发生模式，利用历史退化数据，深度挖掘潜在故障信息，构建基于物理或数据驱动的预测模型，实现对故障程度的预测，是装备健康状态判断，剩余使用寿命分析，以及基于状态维修的基础，装备退化呈现较大差异，前期故障难以提取特征，后期故障难以反应当前装备健康状态，较近时刻故障数据的量又难以满足拟合精度的需求，增加了故障预测的难度，小样本、快收敛、高精度的故障在线预测技术已成为故障预测领域研究的热点；
3.但是当前的故障预测技术的适应能力差，对故障复杂性高、样本数据量少、预测时效性强的装备的故障预测效果差。


技术实现要素：

4.本发明提供一种小样本装备故障在线预测方法，可以有效解决上述背景技术中提出当前的故障预测技术的适应能力差，对故障复杂性高、样本数据量少、预测时效性强的装备的故障预测效果差的问题。
5.为实现上述目的，本发明提供如下技术方案：一种小样本装备故障在线预测方法，包括如下预测步骤：
6.s1、数据处理：通过wtd降噪算法中的阈值函数优化以及引入bic评估分解层数对wtd复杂度的影响，提出一种改进wtd算法用于在线滤除故障信号中的噪声；
7.s2、故障程度识别：通过mest构建系统或装备的非参数模型，通过对观测向量与历史记忆矩阵的最优重构估计获得估计向量，利用估计向量和观测向量之间的差异性来反应故障程度，并引入csfi进行平滑处理；
8.s3、故障在线预测：通过tsfm非统计性分析方法，消除数据的偶然变动，并引入梯度下降在线更新平滑因子，引入自适应滑动时间窗动态截取时间序列数据，提升tsfm的拟合能力；
9.s4、实验分析：验证小样本条件下装备故障预测模型的有效性和可行性；
10.s5、数据处理分析：通过预设bic对不同分解层数的改进wtd算法复杂性进行评估；
11.s6、故障程度识别分析：以改进wtd得到各小波分解系数作为改进mest的观测变量，设置采样频率、健康状态和退化状态进行故障程度的识别分析，采用双重csfi处理数据，消除曲线中导数不存在的“尖点”，得平滑后故障程度点；
12.s7、故障在线测试分析：预设自适应滑动时间窗、自适应平滑因子、学习因子、最大训练迭代次数、最小允许误差，将平滑后的轴承故障程度值输入故障在线预测模型中，得自适应平滑因子变化趋势、自适应滑动时间窗长度变化趋势、预测误差变化趋势。
13.s8、预测结果总结：对在线预测模型故障预测的效果进行总结。
14.2.根据权利要求1所述的一种小样本装备故障在线预测方法，其特征在于，所述s1中，通过改进wtd算法用于在线滤除故障信号中的噪声，其原理如下：
[0015][0016]
式中，λ为阈值；
[0017]
ω
j,k
为故障信号小波系数；
[0018]
为预估小波系数；
[0019]
j为分解尺度，且1≤j≤j，j为最大尺度；
[0020]
sgn()为符号函数；
[0021]
此阈值函数在小波域内具有连续性，当ω
j,k
→
λ-时，当ω
j,k
→
λ
+
时，
[0022]
阈值λ的选择应满足：
[0023][0024]
式中，n表示信号长度；
[0025]
σj表示第j层高斯白噪声标准差，其表达式为：
[0026][0027]
式中，cd
j,k
为第j层小波分解的高频部分；
[0028]
p为该尺度下小波系数的个数；
[0029]
wtd认为故障信号存在于低频部分ad
j,k
中，噪声存在于高频部分cd
j,k
中；
[0030]
由于噪声的振幅服从高斯分布，以最大分解层数高频部分cd
j,k
为评估模型复杂性的数据，引入bic对模型的复杂性进行评估：
[0031]
bic＝qln(n)-2ln(l)
[0032]
式中，q为模型参数个数；
[0033]
n为样本数量；
[0034]
l为服从高斯分布的最大似然函数，即：
[0035][0036]
根据上述技术方案，所述s2中，估计向量和观测向量具体计算如下：
[0037]
设某一时刻t观测到装备中有n个相互关联的变量，将其记为观测变量x
t
，即
[0038]
x
t
＝[x
t,1
,x
t,2
…
x
t,n
]
t
[0039]
式中，x
t,n
为t时刻状态变量的观测值；
[0040]
构建具有m个历史时刻，n个关联状态变量的历史记忆矩阵d，即
[0041][0042]
由历史记忆矩阵d中m个观测向量x
obs
的线性加权，可获得估计向量x
est
，即
[0043]
x
est
＝dw＝w1x1+w2x2…
wmxm[0044]
式中，w＝[w1,w2...wm]
t
为一个m维的权值向量，代表输入观测向量x
obs
与历史记忆矩阵d的相似度，即
[0045][0046]
式中，为非线性运算符，用以代替普通矩阵中的乘积运算；
[0047]
将d
t
与x
obs
之间的马氏距离(md)作为mest中的非线性运算符，即
[0048][0049]
式中，∑-1
为多维随机变量协方差矩阵的逆矩阵；
[0050]
当两个状态矩阵越相似时，其md越小；
[0051]
当两个状态矩阵差异性越大时，其非线性运算结果越大；
[0052]
将(9)式带入(8)式中，可得mest模型估计向量的最终表达式为：
[0053][0054]
通过对比观测向量x
obs
与估计向量x
est
之间的差值，可得反应装备故障程度的残差值ε，即：
[0055]
ε＝x
est-x
obs
[0056]
通过对比各类故障指标的使用范围，选用均方根反映故障程度；
[0057]
通过求n个维度x
est
与x
obs
残差ε的均方根rmsv，即可得到反应装备故障程度指标dr：
[0058][0059]
运用mest获得的装备故障程度指标dr是由多个离散点构成，dr组成的曲线包含多个导数不存在的“尖点”，引入csfi进行平滑处理。
[0060]
根据上述技术方案，所述s3中，tsfm表达式为：
[0061][0062]
式中，dr
t
为原始序列数据；
[0063]
为第t+t次预测的第i次训练自适应平滑因子；
[0064]
为一次平滑值；
[0065]
为二次平滑值；
[0066]
若t表示预测时间，表示t+t时刻第i次训练的预测值，则预测公式为：
[0067][0068]
其中：
[0069][0070]
通过对比各类梯度下降算法的适用范围，选用随机梯度下降(sgd)自适应更新平滑因子则第t+t次预测的第i次训练自适应平滑因子的表达式为：
[0071][0072]
式中，为第t+t次预测的第i-1次训练自适应平滑因子；
[0073]
β为学习因子；
[0074]
为第t+t次预测的第i次训练预测值；
[0075]
dr
t+t
为第t+t个实际值；
[0076]
为第t+t次预测的第i次训练损失函数对α的偏导：
[0077][0078]
式中，为第t+t次预测的第i次训练时间序列数据长度，即自适应滑动时间窗；
[0079]
绘制滑动时间窗示意图，滑动时间窗示意图中data_mode_i为第i次训练的训练数据，data_test_i为第i次训练的测试数据，
[0080][0081][0082]
式中，μ为调整因子；
[0083]
当时间序列数据不能完全反映故障信息时，即损失函数的偏导一直为同方向，则应持续增长或减小，最大按照1+μ倍增长或1-μ倍缩短；
[0084]
当时间序列数据能部分反映故障信息时，即损失函数的偏导不均为同方向，则长度取决于的方向，为正时，长度增长，为负时，长度缩短；
[0085]
选用方差值varv代表最终预测误差
[0086][0087]
式中len(t)为第t+t次预测故障程度指标数据的个数。
[0088]
根据上述技术方案，所述s3中，故障在线预测的流程如下：
[0089]
step1：第t+t次预测的第i次训练，预测模型先将第i-1次预设滑动时间窗划分为第i次的训练数据data_mode_i和测试数据data_test_i，将第i-1次训练的和带入式中预测第i次训练的值；
[0090]
step2：将带入式计算第i次训练的损失函数若满足最大循环次数maxtrans或小于最小误差minerror，则停止训练，若不满足上述条件，则根据式和式更新自适应平滑因子和自适应滑动时间窗重复step1，开始第t+t次预测的i+1次训练。
[0091]
step3：将step2得到的带入式，求得最终值，并根据式计算预测误差e。
[0092]
根据上述技术方案，所述s4中，以轴承性能退化数据为验证对象，设置采样频率为25.6khz/min，radial force为12kn，转速为2100rpm，运行157.44s，选用水平方向上的振动信号反映被测轴承的故障程度。
[0093]
根据上述技术方案，所述s5中，不同分阶层数wtd效果对比如下：
[0094]
当小波分解层数j＝7时，bic值最小；
[0095]
当小波分解层数j＜7时，随着j的增大，改进wtd降噪的效果也随之增强；
[0096]
当小波分解层数j＞7时，随着j的增大，改进wtd降噪的效果增强不明显；
[0097]
当j＝7时，改进wtd算法既能保证降噪效果满足要求，又能有效防止因精度过高造成的模型复杂度过高的问题；
[0098]
以噪声振幅服从高斯分布为评估模型复杂度的依据，因此将最大分解层数高频部分cd
j,k
输出，绘制高斯白噪声振幅分布图；
[0099]
当j＜7时，改进wtd最大分解层数高频部分cd
j,k
服从高斯分布，此时改进wtd中的阈值还未破坏高斯白噪声的分布特性，即小波分解不充分；
[0100]
当j≥7时，改进wtd最大分解层数高频部分cd
j,k
不服从高斯分布，此时改进wtd中的阈值破坏了高斯白噪声的分布特性，即小波分解充分，说明了设置bic评估分解层数的可行性和预设最高小波分解层数的高频部分cd
j,k
为评估数据的科学性；
[0101]
设置小波分解层数j＝7，以改进wtd对轴承性能退化数据进行降噪处理。
[0102]
根据上述技术方案，所述s6中，以改进wtd得到各小波分解系数作为改进mest的观测变量x
t
，即当j＝7时，x
t
＝[ca
t,7
,cd
t,7
,cd
t,6
,cd
t,5
,cd
t,4
,cd
t,3
,cd
t,2
,cd
t,1
]，设置采样频率fs＝25600hz，即1秒，将轴承性能退化数据前2秒设置为健康状态，即历史时刻m＝2，历史记忆矩阵d＝[x1,x2]，其他时刻每2秒设置为性能退化状态，即观测向量x
obs
，将轴承性能退化数据带入，共计评估77个故障程度dr；
[0103]
同时，为了消除dr曲线中多个导数不存在的“尖点”，采用双重csfi处理数据，插值数为原数据量的10倍，共得7700个平滑后故障程度点。
[0104]
根据上述技术方案，所述s7中，预设自适应滑动时间窗length初始值为100，自适应平滑因子α初始值为0.05，学习因子β为0.5，最大训练迭代次数maxtrans为1000次，最小允许误差minerror为10-8
，将平滑后的7700个轴承故障程度dr值输入故障在线预测模型中，比较自适应平滑因子α变化趋势，自适应滑动时间窗length长度变化趋势，以及预测误差e变化趋势。
[0105]
根据上述技术方案，所述s8中，预测结果总结如下：
[0106]
bic能够准确找出wtd算法的最佳分解层数；
[0107]
改进wtd算法降噪效果优异；
[0108]
md改进的mest算法和双重csfi算法能有效将故障信号转化为故障程度指标；
[0109]
提出的在线预测模型，能够实现对装备故障的在线预测。
[0110]
与现有技术相比，本发明的有益效果：
[0111]
通过提出小样本装备故障在线预测模型，并通过滚动轴承生命周期振动数据验证了本模型的有效性和可靠性，bic能够准确找出wtd算法的最佳分解层数，为改进wtd模型参数设置提供依据，改进wtd算法降噪效果优异，保证了故障数据的可靠性，md改进的mest算法和双重csfi算法能有效将故障信号转化为故障程度指标，为后续故障预测提供了高质量的数据保障，提出的具备滑动时间窗自适应截取故障数据长度和参数自适应更新的故障在线预测模型，能够不断挖掘数据中潜在的故障程度信息，实现对装备故障的在线预测，满足对故障复杂性高、样本数据量少、预测时效性强的额装备的故障预测。
附图说明
[0112]
附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例一起用于解释本发明，并不构成对本发明的限制。在附图中：
[0113]
图1是本发明装备故障在线预测模型步骤图；
[0114]
图2是本发明wtd原理图；
[0115]
图3是本发明滑动时间窗示意图；
[0116]
图4是本发明故障在线预测流程图；
[0117]
图5是本发明轴承疲劳试验台结构图；
[0118]
图6是本发明实验轴承振幅信号图；
[0119]
图7是本发明bic值变化图；
[0120]
图8是本发明不同分解层数wtd效果图；
[0121]
图9是本发明高斯白噪声振幅分布图；
[0122]
图10是本发明各种wtd效果对比图；
[0123]
图11是本发明轴承全寿命dr值变化趋势图；
[0124]
图12是本发明自适应α取值变化趋势图；
[0125]
图13是本发明自适应滑动时间窗长度变化趋势图；
[0126]
图14是本发明预测误差变化趋势图；
[0127]
图15是本发明预测拟合效果图；
[0128]
图16是本发明故障测试的步骤图。
具体实施方式
[0129]
以下结合附图对本发明的优选实施例进行说明，应当理解，此处所描述的优选实施例仅用于说明和解释本发明，并不用于限定本发明。
[0130]
实施例：如图1和16所示，本发明提供一种技术方案，一种小样本装备故障在线预测方法，包括如下预测步骤：
[0131]
s1、数据处理：
[0132]
常用的降噪算法有经验模态分解(emd)和wtd两种；
[0133]
emd基本原理是将故障信号分解成各阶的本征模态函数(imf)，再从中提取故障信号的特征，但emd算法在含有大量噪声的背景下，存在模态混叠和末端效应等问题；
[0134]
wtd主要包括硬阈值降噪法和软阈值降噪法两种，wtd虽然存在小波分解层数依靠
经验确定，硬阈值函数不连续，易造成降噪后产生振荡，软阈值降噪小波系数存在偏差，易导致重构信号误差大的问题，但由于其计算量小，降噪效率高，其和其改进算法广泛应用于故障在线预测中，如通过改进小波阈值函数，有效地滤除了信号中的噪声，提高了输出信号的信噪比，但其存在降噪效果有待提高和小波分解层数如何确定的问题；
[0135]
因此通过对wtd中的阈值函数优化以及引入bic评估分解层数对wtd复杂度的影响，提出一种改进wtd算法用于在线滤除故障信号中的噪声，其原理如下：
[0136][0137]
式中，λ为阈值；
[0138]
ω
j,k
为故障信号小波系数；
[0139]
为预估小波系数；
[0140]
j为分解尺度，且1≤j≤j，j为最大尺度；
[0141]
sgn()为符号函数；
[0142]
此阈值函数在小波域内具有连续性，当ω
j,k
→
λ-时，当ω
j,k
→
λ
+
时，
[0143]
阈值λ的选择应满足：
[0144][0145]
式中，n表示信号长度；
[0146]
σj表示第j层高斯白噪声标准差，其表达式为：
[0147][0148]
式中，cd
j,k
为第j层小波分解的高频部分；
[0149]
p为该尺度下小波系数的个数；
[0150]
如图2所示，wtd认为故障信号存在于低频部分ad
j,k
中，噪声存在于高频部分cd
j,k
中；
[0151]
由于噪声的振幅服从高斯分布，以最大分解层数高频部分cd
j,k
为评估模型复杂性的数据，引入bic对模型的复杂性进行评估，确保在满足降噪精度的条件下，有效降低模型的复杂度
[0152]
bic＝q ln(n)-2 ln(l)
[0153]
式中，q为模型参数个数；
[0154]
n为样本数量；
[0155]
l为服从高斯分布的最大似然函数，即
[0156][0157]
s2、故障程度识别：
[0158]
mest的核心思想是构建系统或装备的非参数模型，通过对观测向量与历史记忆矩阵的最优重构估计获得估计向量，利用估计向量和观测向量之间的差异性来反应系统或装备的故障程度，与神经网络相比，具有计算量小、模型物理意义明确、模型结构简单的优点；
[0159]
设某一时刻t观测到装备中有n个相互关联的变量，将其记为观测变量x
t
，即
[0160]
x
t
＝[x
t,1
,x
t,2
…
x
t,n
]
t
[0161]
式中，x
t,n
为t时刻状态变量的观测值；
[0162]
构建具有m个历史时刻，n个关联状态变量的历史记忆矩阵d，即
[0163][0164]
由历史记忆矩阵d中m个观测向量x
obs
的线性加权，可获得估计向量x
est
，即
[0165]
x
est
＝dw＝w1x1+w2x2…
wmxm[0166]
式中，w＝[w1,w2…
wm]
t
为一个m维的权值向量，代表输入观测向量x
obs
与历史记忆矩阵d的相似度，即
[0167][0168]
式中，为非线性运算符，用以代替普通矩阵中的乘积运算，避免所产生的不可逆现象，扩大式的适应范围；
[0169]
为提升mest算法的多维数据处理能力，将d
t
与x
obs
之间的马氏距离(md)作为mest中的非线性运算符，即
[0170][0171]
式中，∑-1
为多维随机变量协方差矩阵的逆矩阵，直观可以看出，当两个状态矩阵越相似时，其md越小；
[0172]
当两个状态矩阵差异性越大时，其非线性运算结果越大；
[0173]
将(9)式带入(8)式中，可得mest模型估计向量的最终表达式为
[0174][0175]
通过对比观测向量x
obs
与估计向量x
est
之间的差值，可直观得到反应装备故障程度的残差值ε，即
[0176]
ε＝x
est-x
obs
[0177]
通过对比各类故障指标的使用范围，选用均方根反映故障程度，如表1所示：
[0178]
表1误差指标表
[0179][0180]
其中，varv为方差值；
[0181]
rmsv为均方根；
[0182]
sf为形状因子；
[0183]
mf为边际因子；
[0184]
e为能量；
[0185]
se为香农熵；
[0186]
re为仁义熵；
[0187]
te为察里斯熵；
[0188]
通过求n个维度x
est
与x
obs
残差ε的rmsv，即可得到反应装备故障程度指标dr：
[0189][0190]
由于运用mest获得的装备故障程度指标dr是由多个离散点构成，不可避免造成dr组成的曲线包含多个导数不存在的“尖点”，大幅增加了故障预测的难度，csfi能够直接利用各节点连接处的光滑性与连接条件构建关系式，确定待定系数，构造出插值多项式，完成使存在“尖点”的曲线变光滑的操作，具备构造简单，计算量小，可靠性高的优点，因此引入csfi进行平滑处理。
[0191]
s3、故障在线预测：
[0192]
tsfm是一种非统计性分析方法，能够消除数据的偶然变动，提高近期数据在预测中的重要程度，符合数据在故障预测中的客观规律，并且由于其数据需求量小，计算过程简单，模型构造方便，特别适用于短期或中长期的故障在线预测，但由于其缺乏对数据转折点的鉴别能力，以及平滑因子的选择过于主观性，严重限制了其故障预测的效果，因此引入梯度下降在线更新平滑因子，引入自适应滑动时间窗动态截取时间序列数据，提升tsfm的拟合能力，tsfm表达式为：
[0193][0194]
式中，dr
t
为原始序列数据；
[0195]
为第t+t次预测的第i次训练自适应平滑因子；
[0196]
为一次平滑值；
[0197]
为二次平滑值；
[0198]
若t表示预测时间，表示t+t时刻第i次训练的预测值，则预测公式为：
[0199][0200]
其中：
[0201][0202]
通过对比各类梯度下降算法的适用范围，选用随机梯度下降(sgd)自适应更新平滑因子如表2所示：
[0203]
表2梯度下降算法对比表
[0204][0205]
其中，bgd为批量梯度下降(bgd)，momentum-sgd为带动量的随机梯度下降(sgdm)；
[0206]
则第t+t次预测的第i次训练自适应平滑因子的表达式为：
[0207][0208]
式中，为第t+t次预测的第i-1次训练自适应平滑因子；
[0209]
β为学习因子；
[0210]
为第t+t次预测的第i次训练预测值；
[0211]
dr
t+t
为第t+t个实际值；
[0212]
为第t+t次预测的第i次训练损失函数对α的偏导：
[0213][0214]
如图3所示，式中，为第t+t次预测的第i次训练时间序列数据长度，即自适应滑动时间窗；
[0215]
图6中data_mode_i是指第i次训练的训练数据，data_test_i是指第i次训练的测试数据，
[0216][0217]
式中，μ为调整因子；
[0218]
由和(18)式易得，当时间序列数据不能完全反映故障信息时，即损失函数的偏导一直为同方向，则应持续增长或减小，以满足提取故障信息的需求，最大按照1+μ倍增长或1-μ倍缩短；
[0219]
当时间序列数据能部分反映故障信息时，即损失函数的偏导不均为同方向，则长度取决于的方向，为正时，长度增长，为负时，长度缩短。
[0220]
选用中varv代表最终预测误差
[0221][0222]
式中len(t)为第t+t次预测故障程度指标数据的个数。
[0223]
如图4所示，step1：第t+t次预测的第i次训练，预测模型先将第i-1次预设滑动时间窗划分为第i次的训练数据data_mode_i和测试数据data_test_i，将第i-1次训练的和带入式中预测第i次训练的值；
[0224]
step2：将带入式计算第i次训练的损失函数若满足最大循环次数maxtrans或小于最小误差minerror，则停止训练，若不满足上述条件，则根据式和式更新自适应平滑因子和自适应滑动时间窗重复step1，开始第t+t次预测的i+1次训练；
[0225]
step3：将step2得到的带入式，求得最终值，并根据式计算预测误差e。
[0226]
s4、实验分析：
[0227]
为进行一步验证小样本条件下装备故障预测模型的有效性和可行性，利用滚动轴承性能退化数据进行分析与验证；
[0228]
如图5所示，为了获得型号为ldk uer204的轴承故障数据，将两块pcb 352c33加速度计放置在被测轴承的外壳上，彼此之间的夹角为90
°
，即一个放置在横轴，另一个放置在纵轴，设置采样频率为25.6khz/min，radial force为12kn，转速为2100rpm，运行157.44s，由于载荷是施加在水平方向上的，此方向上的加速度计更能准确反映被测轴承的退化信息，因此，选用水平方向上的振动信号反映被测轴承的故障程度；
[0229]
如图6所示，在轴承性能退化实验前期，振幅变化不大，轴承处于正常运行阶段，即a点前，但此时包含有大量噪声信号；
[0230]
在性能退化实验后期，轴承振幅随工作时间的增加而增大，轴承处于性能退化阶段，即a点和c点之间，被检测信号中包含有大量的轴承性能退化信号和噪声信号；
[0231]
在性能退化实验末期，轴承振幅陡然增大，轴承处于故障阶段，即c点之后；
[0232]
因此，是否能剔除噪声信号影响，准确判断轴承性能退化时刻a和发生故障时刻c，真实反映轴承退化性能的趋势，保证较高的预测精度和较短的预测时间是验证所提出模型的重点。
[0233]
s5、数据处理分析：
[0234]
如图7和8所示，通过预设bic对不同分解层数的改进wtd算法复杂性进行评估；
[0235]
当小波分解层数j＝7时，bic值最小；
[0236]
当小波分解层数j＜7时，随着j的增大，改进wtd降噪的效果也随之增强；
[0237]
当小波分解层数j＞7时，随着j的增大，改进wtd降噪的效果增强不明显；
[0238]
因此，当j＝7时，改进wtd算法既能保证降噪效果满足要求，又能有效防止因精度过高造成的模型复杂度过高的问题；
[0239]
如图9所示，以噪声振幅服从高斯分布为评估模型复杂度的依据，因此将最大分解层数高频部分cd
j,k
输出；
[0240]
当j＜7时，改进wtd最大分解层数高频部分cd
j,k
服从高斯分布，此时改进wtd中的阈值还未破坏高斯白噪声的分布特性，即小波分解不充分；
[0241]
当j≥7时，改进wtd最大分解层数高频部分cd
j,k
不服从高斯分布，此时改进wtd中的阈值破坏了高斯白噪声的分布特性，即小波分解充分；
[0242]
因此，说明了设置bic评估分解层数的可行性和预设最高小波分解层数的高频部分cd
j,k
为评估数据的科学性；
[0243]
如图10所示，设置小波分解层数j＝7，以提出的改进wtd和改进wtd(简称引用改进wtd)对轴承性能退化数据进行降噪处理；
[0244]
图10中(1)、(4)、(2)、(5)、(3)、(6)分别为轴承性能退化原始时域图、原始频谱图、引用改进wtd降噪后时域图、引用改进wtd降噪后频谱图、改进wtd降噪后时域图、改进wtd降噪后频域图；
[0245]
提出的改进wtd算法降噪效果优于引用改进wtd算法，即能大幅抑制正常运行阶段噪声的影响，又能保留性能退化阶段和故障阶段的轴承故障程度信息，有利于准确判断性能退化阶段和故障发生阶段的开始时刻。
[0246]
s6、故障程度识别分析：
[0247]
如图11所示，以改进wtd得到各小波分解系数作为改进mest的观测变量x
t
，即当j＝7时，x
t
＝[ca
t,7
,cd
t,7
,cd
t,6
,cd
t,5
,cd
t,4
,cd
t,3
,cd
t,2
,cd
t,1
]，设置采样频率fs＝25600hz，即1秒，将轴承性能退化数据前2秒设置为健康状态，即历史时刻m＝2，历史记忆矩阵d＝[x1,x2]，其他时刻每2秒设置为性能退化状态，即观测向量x
obs
，将轴承性能退化数据带入，共计评估77个故障程度dr；
[0248]
同时，为了消除dr曲线中多个导数不存在的“尖点”，采用双重csfi处理数据，插值数为原数据量的10倍，共得7700个平滑后故障程度点；
[0249]
图11中(1)、(2)、(3)分别为轴承性能退化曲线图，故障程度dr曲线图，以及平滑后故障程度dr曲线图；
[0250]
由图11中(1)和(2)可知当轴承在a点性能发生退化，故障程度dr同时增大，改进mest模型能够准确判断性能退化阶段的开始时刻，并且能于故障发生点c前，即将发生故障点b时增大故障程度dr值，给出预警信号；
[0251]
由图11中(2)和(3)可知，采用的双重csfi算法能够不改变原故障程度dr曲线的趋势，消除曲线中导数不存在的尖点，增强数据的可预测性。
[0252]
s7、故障在线预测分析：
[0253]
实验硬件cpu为i7-10875h，运行内存为32gb，显卡为rtx2060，预设自适应滑动时间窗length初始值为100，自适应平滑因子α初始值为0.05，学习因子β为0.5，最大训练迭代次数maxtrans为1000次，最小允许误差minerror为10-8
，将平滑后的7700个轴承故障程度dr值输入故障在线预测模型中，则自适应平滑因子α变化趋势，自适应滑动时间窗length长度变化趋势，以及预测误差e变化趋势如图所示。
[0254]
如图12所示，在正常运行阶段，即a点前，轴承故障程度dr值变化幅度较小，此阶段自适应平滑因子α除前期短暂调整外，迅速稳定在0.5左右运行；
[0255]
在性能退化阶段，即a点后，随着故障程度dr增大，此阶段自适应平滑因子α变化较为明显，并在检测到轴承将发生故障的时刻b点，自适应平滑因子α调整较为剧烈，表明所提出的自适应平滑因子α能随着故障程度dr实时变化，准确判断性能退化时刻a点和即将完全故障时刻b点，寻优能力强，满足故障在线预测的需求。
[0256]
如图13所示，在正常运行阶段，即a点前，轴承故障程度dr值变化幅度较小，此阶段自适应滑动时间窗length稳定在初始值100左右；
[0257]
在性能退化阶段，即a点后，随着故障程度dr增大，正常运行阶段的性能退化数据已不能表征此时轴承的故障程度，甚至产生负效应，因此，自适应滑动时间窗length长度迅速缩短，并在检测到轴承即将发生故障的时刻b点开始产生小范围波动，表明所提出的自适应滑动时间窗length能随着故障程度dr实时变化，在平稳阶段，保持滑动时间窗长度稳定，在变化剧烈阶段，调整滑动时间窗长度迅速，能够准确判断性能退化时刻a点和将发生故障时刻b点，满足故障在线预测的需求；
[0258]
如图14和15所示，预测误差e的变化趋势与轴承故障程度dr的变化趋势相近，即误差随着故障程度dr趋势的复杂程度逐渐增大，到达即将发生故障时刻b点误差增大较为明显，符合自适应平滑因子α和自适应滑动时间窗length的变化趋势，预测值与实际值几乎完全重合，最终预测误差e＝0.068％，拟合精度高，满足故障在线预测的需求。
[0259]
此实验条件下，单次平均预测时间约为0.0277s，单次平均预测时间约为两次故障程度dr间隔时间的1.385％，即模型能在较短时间内给出装备故障程度的预测结果，预测时间较短，满足故障预测的需求。
[0260]
s8、预测结果总结：
[0261]
针对故障复杂性高、样本数据量少、预测时效性强的装备故障预测问题，提出一种小样本装备故障在线预测模型，并通过滚动轴承生命周期振动数据验证了本模型的有效性和可靠性，预测结果总结如下：
[0262]
(1)bic能够准确找出wtd算法的最佳分解层数，为改进wtd模型参数设置提供依据；
[0263]
(2)改进wtd算法降噪效果优异，保证了故障数据的可靠性；
[0264]
(3)md改进的mest算法和双重csfi算法能有效将故障信号转化为故障程度指标，为后续故障预测提供了高质量的数据保障；
[0265]
(4)提出的具备滑动时间窗自适应截取故障数据长度和参数自适应更新的故障在线预测模型，能够不断挖掘数据中潜在的故障程度信息，实现对装备故障的在线预测。
[0266]
最后应说明的是：以上所述仅为本发明的优选实例而已，并不用于限制本发明，尽管参照前述实施例对本发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

技术特征：
1.一种小样本装备故障在线预测方法，其特征在于：包括如下预测步骤：s1、数据处理：通过wtd降噪算法中的阈值函数优化以及引入bic评估分解层数对wtd复杂度的影响，提出一种改进wtd算法用于在线滤除故障信号中的噪声；s2、故障程度识别：通过mest构建系统或装备的非参数模型，通过对观测向量与历史记忆矩阵的最优重构估计获得估计向量，利用估计向量和观测向量之间的差异性来反应故障程度，并引入csfi进行平滑处理；s3、故障在线预测：通过tsfm非统计性分析方法，消除数据的偶然变动，并引入梯度下降在线更新平滑因子，引入自适应滑动时间窗动态截取时间序列数据，提升tsfm的拟合能力；s4、实验分析：验证小样本条件下装备故障预测模型的有效性和可行性；s5、数据处理分析：通过预设bic对不同分解层数的改进wtd算法复杂性进行评估；s6、故障程度识别分析：以改进wtd得到各小波分解系数作为改进mest的观测变量，设置采样频率、健康状态和退化状态进行故障程度的识别分析，采用双重csfi处理数据，消除曲线中导数不存在的“尖点”，得平滑后故障程度点；s7、故障在线测试分析：预设自适应滑动时间窗、自适应平滑因子、学习因子、最大训练迭代次数、最小允许误差，将平滑后的轴承故障程度值输入故障在线预测模型中，得自适应平滑因子变化趋势、自适应滑动时间窗长度变化趋势、预测误差变化趋势。s8、预测结果总结：对在线预测模型故障预测的效果进行总结。2.根据权利要求1所述的一种小样本装备故障在线预测方法，其特征在于，所述s1中，通过改进wtd算法用于在线滤除故障信号中的噪声，其原理如下：式中，λ为阈值；ω
j,k
为故障信号小波系数；为预估小波系数；j为分解尺度，且1≤j≤j，j为最大尺度；sgn()为符号函数；此阈值函数在小波域内具有连续性，当ω
j,k
→
λ-时，当ω
j,k
→
λ
+
时，阈值λ的选择应满足：式中，n表示信号长度；σ
j
表示第j层高斯白噪声标准差，其表达式为：
式中，cd
j,k
为第j层小波分解的高频部分；p为该尺度下小波系数的个数；wtd认为故障信号存在于低频部分ad
j,k
中，噪声存在于高频部分cd
j,k
中；由于噪声的振幅服从高斯分布，以最大分解层数高频部分cd
j,k
为评估模型复杂性的数据，引入bic对模型的复杂性进行评估：bic＝qln(n)-2ln(l)式中，q为模型参数个数；n为样本数量；l为服从高斯分布的最大似然函数，即：3.根据权利要求1所述的一种小样本装备故障在线预测方法，其特征在于，所述s2中，估计向量和观测向量具体计算如下：设某一时刻t观测到装备中有n个相互关联的变量，将其记为观测变量x
t
，即x
t
＝[x
t,1
,x
t,2
…
x
t,n
]
t
式中，x
t,n
为t时刻状态变量的观测值；构建具有m个历史时刻，n个关联状态变量的历史记忆矩阵d，即由历史记忆矩阵d中m个观测向量x
obs
的线性加权，可获得估计向量x
est
，即x
est
＝dw＝w1x1+w2x2…
w
m
x
m
式中，w＝[w1,w2…
w
m
]
t
为一个m维的权值向量，代表输入观测向量x
obs
与历史记忆矩阵d的相似度，即式中，为非线性运算符，用以代替普通矩阵中的乘积运算；将d
t
与x
obs
之间的马氏距离(md)作为mest中的非线性运算符，即式中，∑-1
为多维随机变量协方差矩阵的逆矩阵；当两个状态矩阵越相似时，其md越小；当两个状态矩阵差异性越大时，其非线性运算结果越大；将(9)式带入(8)式中，可得mest模型估计向量的最终表达式为：通过对比观测向量x
obs
与估计向量x
est
之间的差值，可得反应装备故障程度的残差值ε，即：ε＝x
est-x
obs
通过对比各类故障指标的使用范围，选用均方根反映故障程度；
通过求n个维度x
est
与x
obs
残差ε的均方根rmsv，即可得到反应装备故障程度指标dr：运用mest获得的装备故障程度指标dr是由多个离散点构成，dr组成的曲线包含多个导数不存在的“尖点”，引入csfi进行平滑处理。4.根据权利要求1所述的一种小样本装备故障在线预测方法，其特征在于，所述s3中，tsfm表达式为：式中，dr
t
为原始序列数据；为第t+t次预测的第i次训练自适应平滑因子；为一次平滑值；为二次平滑值；若t表示预测时间，表示t+t时刻第i次训练的预测值，则预测公式为：其中：通过对比各类梯度下降算法的适用范围，选用随机梯度下降(sgd)自适应更新平滑因子则第t+t次预测的第i次训练自适应平滑因子的表达式为：式中，为第t+t次预测的第i-1次训练自适应平滑因子；β为学习因子；为第t+t次预测的第i次训练预测值；dr
t+t
为第t+t个实际值；为第t+t次预测的第i次训练损失函数对α的偏导：式中，为第t+t次预测的第i次训练时间序列数据长度，即自适应滑动时间窗；绘制滑动时间窗示意图，滑动时间窗示意图中data_mode_i为第i次训练的训练数据，data_test_i为第i次训练的测试数据，data_test_i为第i次训练的测试数据，
式中，μ为调整因子；当时间序列数据不能完全反映故障信息时，即损失函数的偏导一直为同方向，则应持续增长或减小，最大按照1+μ倍增长或1-μ倍缩短；当时间序列数据能部分反映故障信息时，即损失函数的偏导不均为同方向，则长度取决于的方向，为正时，长度增长，为负时，长度缩短；选用方差值varv代表最终预测误差式中len(t)为第t+t次预测故障程度指标数据的个数。5.根据权利要求4所述的一种小样本装备故障在线预测方法，其特征在于，所述s3中，故障在线预测的流程如下：step1：第t+t次预测的第i次训练，预测模型先将第i-1次预设滑动时间窗划分为第i次的训练数据data_mode_i和测试数据data_test_i，将第i-1次训练的和带入式中预测第i次训练的值；step2：将带入式计算第i次训练的损失函数若满足最大循环次数maxtrans或小于最小误差minerror，则停止训练，若不满足上述条件，则根据式和式更新自适应平滑因子和自适应滑动时间窗重复step1，开始第t+t次预测的i+1次训练。step3：将step2得到的带入式，求得最终值，并根据式计算预测误差e。6.根据权利要求1所述的一种小样本装备故障在线预测方法，其特征在于，所述s4中，以轴承性能退化数据为验证对象，设置采样频率为25.6khz/min，radial force为12kn，转速为2100rpm，运行157.44s，选用水平方向上的振动信号反映被测轴承的故障程度。7.根据权利要求1所述的一种小样本装备故障在线预测方法，其特征在于，所述s5中，不同分阶层数wtd效果对比如下：当小波分解层数j＝7时，bic值最小；当小波分解层数j＜7时，随着j的增大，改进wtd降噪的效果也随之增强；当小波分解层数j＞7时，随着j的增大，改进wtd降噪的效果增强不明显；当j＝7时，改进wtd算法既能保证降噪效果满足要求，又能有效防止因精度过高造成的模型复杂度过高的问题；以噪声振幅服从高斯分布为评估模型复杂度的依据，因此将最大分解层数高频部分cd
j,k
输出，绘制高斯白噪声振幅分布图；当j＜7时，改进wtd最大分解层数高频部分cd
j,k
服从高斯分布，此时改进wtd中的阈值还未破坏高斯白噪声的分布特性，即小波分解不充分；当j≥7时，改进wtd最大分解层数高频部分cd
j,k
不服从高斯分布，此时改进wtd中的阈值破坏了高斯白噪声的分布特性，即小波分解充分，说明了设置bic评估分解层数的可行性和预设最高小波分解层数的高频部分cd
j,k
为评估数据的科学性；
设置小波分解层数j＝7，以改进wtd对轴承性能退化数据进行降噪处理。8.根据权利要求1所述的一种小样本装备故障在线预测方法，其特征在于，所述s6中，以改进wtd得到各小波分解系数作为改进mest的观测变量x
t
，即当j＝7时，x
t
＝[ca
t,7
,cd
t,7
,cd
t,6
,cd
t,5
,cd
t,4
,cd
t,3
,cd
t,2
,cd
t,1
]，设置采样频率fs＝25600hz，即1秒，将轴承性能退化数据前2秒设置为健康状态，即历史时刻m＝2，历史记忆矩阵d＝[x1,x2]，其他时刻每2秒设置为性能退化状态，即观测向量x
obs
，将轴承性能退化数据带入，共计评估77个故障程度dr；同时，为了消除dr曲线中多个导数不存在的“尖点”，采用双重csfi处理数据，插值数为原数据量的10倍，共得7700个平滑后故障程度点。9.根据权利要求1所述的一种小样本装备故障在线预测方法，其特征在于，所述s7中，预设自适应滑动时间窗length初始值为100，自适应平滑因子α初始值为0.05，学习因子β为0.5，最大训练迭代次数maxtrans为1000次，最小允许误差minerror为10-8
，将平滑后的7700个轴承故障程度dr值输入故障在线预测模型中，比较自适应平滑因子α变化趋势，自适应滑动时间窗length长度变化趋势，以及预测误差e变化趋势。10.根据权利要求1所述的一种小样本装备故障在线预测方法，其特征在于，所述s8中，预测结果总结如下：bic能够准确找出wtd算法的最佳分解层数；改进wtd算法降噪效果优异；md改进的mest算法和双重csfi算法能有效将故障信号转化为故障程度指标；提出的在线预测模型，能够实现对装备故障的在线预测。

技术总结
本发明公开了一种小样本装备故障在线预测方法，包括如下预测步骤：S1、数据处理：通过WTD降噪算法中的阈值函数优化以及引入BIC评估分解层数对WTD复杂度的影响，提出一种改进WTD算法用于在线滤除故障信号中的噪声，本发明通过提出小样本装备故障在线预测模型，并通过滚动轴承生命周期振动数据验证了本模型的有效性和可靠性，BIC能够准确找出WTD算法的最佳分解层数，为改进WTD模型参数设置提供依据，改进WTD算法降噪效果优异，保证了故障数据的可靠性，MD改进的MEST算法和双重CSFI算法能有效将故障信号转化为故障程度指标，为后续故障预测提供了高质量的数据保障。预测提供了高质量的数据保障。预测提供了高质量的数据保障。


技术研发人员：张保山 郭基联 周峰 周章文 张明亮 李波 魏圣军 顾金玲
受保护的技术使用者：中国人民解放军空军工程大学
技术研发日：2023.04.06
技术公布日：2023/7/7