基于VETMRRN的空间滚动轴承剩余寿命预测方法

基于vetmrrn的空间滚动轴承剩余寿命预测方法
技术领域
1.本发明涉及一种基于vetmrrn的空间滚动轴承剩余寿命预测方法。


背景技术：

2.空间滚动轴承作为空间运动机构的关键零部件，其寿命和可靠性在很大程度上影响了空间站、卫星等空间飞行器的机械部件能否正常运转、实现预定功能和达到预期服役寿命。因此，预测空间滚动轴承剩余寿命对于适应空间飞行器的高可靠、长寿命，避免或减少飞行器产生故障具有重要意义。然而，由于空间在轨环境中难以获取轴承运行状态信息，同时空间环境与地面常规环境之间存在较大区别，因此国内外航天工业部门和相关研究机构通常采用地面模拟空间环境开展空间滚动轴承寿命鉴定试验来预测空间滚动轴承剩余寿命。在(模拟)空间环境下，空间滚动轴承的失效一般是由于固体润滑膜的磨损造成的精度失效(不是如裂纹等破坏性故障)，相比于早期或典型破坏性故障状态，其有效寿命阶段内的振动特征更为微弱。其次，模拟空间环境的设备(如真空泵)空间狭小，且在该狭小空间内多组轴承同时试验，导致空间轴承振动信号受到强烈的环境噪声干扰，这些干扰成分与空间轴承真实振动信号的低频段分量混叠、耦合。另一方面，空间滚动轴承一般工作在变工况下：模拟空间环境下滚动轴承往往采取径向无加载、轴向加载方式开展寿命鉴定试验，轴承振动受轴向载荷大小影响明显，随着轴承磨损，轴承间隙改变，轴向载荷将发生持续性的变化；同时，由于真空室容量限制，难以实施电机转速的闭环控制，轴承转速波动。以上地面模拟空间环境下空间滚动轴承寿命试验的特殊性使空间滚动轴承的运行、退化和失效过程与地面常规环境下的传统滚动轴承有明显差异，导致地面模拟空间环境下进行空间滚动轴承剩余寿命预测具有较大的挑战性。
3.将人工智能引入地面模拟空间环境下空间滚动轴承寿命预测已经成为全新的研究热点，由于起步不久，目前的相关研究案例为数不多。张焱等采用多尺度变异粒子群优化多核最小二乘支持向量机(multi-kernel least-square support vector machine,mk-lssvm)对地面模拟空间状态下固体润滑滚动轴承寿命进行预测；陈仁祥等采用加噪样本扩展深度稀疏自编码神经网络(dsae-nn)识别真空泵中固体润滑滚动轴承的寿命阶段；董绍江等结合特征噪声能量比、自相关函数和深度残差收缩网络进行模拟空间环境下空间滚动轴承剩余寿命预测；高宁等采用自适应kriging模型作为无偏估计模型，同时采用改进自我中心网络(iego)作为学习算法加速构建模型来预测航空齿轮泵滑动轴承的剩余寿命和可靠性。然而，以上方法均存在一些缺陷。对于拥有较多参数的多尺度变异粒子群优化mk-lssvm，它采取多尺度变异粒子群优化这种遍历优化法来选取lssvm参数，这种做法时间复杂度偏高、运行效率低。dsae-nn等深度学习方法擅长特征缩放，但推理能力较弱。iego难以处理较长时间跨度的动态图模式的空间滚动轴承性能退化数据，换言之，它不能可靠地进行较长时间跨度(即较多时间步)的性能退化趋势和剩余寿命预测。
4.长短时记忆神经网络(long short-term memory,lstm)等经典时间循环神经网络(sequence recurrent neural networks,srnns)具有长时记忆等优点，在解决剩余寿命预
测方面是相对更好的选择，但srnns在时间序列的预测中存在长时依赖的缺陷，因此泛化能力较差；同时，srnns通过遍历整个训练数据集以实现有监督式学习的做法增加了时间复杂度。


技术实现要素：

5.本发明的目的在于克服现有技术的不足，提供一种可以实现较高的预测精度、较好的泛化性能和较高的计算效率的基于vetmrrn的空间滚动轴承剩余寿命预测方法。
6.本发明的目的是通过以下技术方案来实现的：基于vetmrrn的空间滚动轴承剩余寿命预测方法，包括以下步骤：
7.s1、从空间滚动轴承的原始振动加速度数据中提取时、频域特征，进行shapely值特征融合，得到shapely值融合特征作为空间滚动轴承性能退化特征；
8.s2、将得到的空间滚动轴承性能退化特征输入到vetmrrn中用于训练vetmrrn的超参数和网络参数；
9.s3、利用训练好的vetmrrn来多步预测空间滚动轴承以后的性能退化特征趋势；
10.s4、结合训练用的空间滚动轴承性能退化特征和预测得到的性能退化特征趋势，建立威布尔分布可靠度模型，预测出空间滚动轴承的精度失效阈值时间点和剩余寿命。
11.进一步地，所述步骤s1中提取的时、频域特征包括平均幅值、均方幅值、峰值、裕度、偏斜度、峭度、能量谱、倒频谱和功率谱。
12.进一步地，步骤s2中所述vetmrrn为变分资格迹元强化循环网络，vetmrrn包括l层网络结构，每层网络结构中包括n-k+1个时序运算单元；
13.假设vetmrrn的初始隐层状态为h0；vetmrrn的输入时间序列为x＝[x1,x2,
…
,xn]，将该时间序列重构为矩阵形式：
[0014][0015]
式中，k为回顾序列长度，表示vetmrrn完成时间序列预测所需回顾的历史时间步长，也就是输入数据维度；将每个行向量视作当前时刻的输入时间序列，即x
t
表示一次预测过程中t时刻的输入时间序列，1≤t≤n-k+1；
[0016]
x1,x2,
…
,x
n-k
,x
n-k+1
分别输入第一层网络结构的n-k+1个时序运算单元；初始化隐层状态h0＝[0,0,
…
,0]，并分别输入每层网络结构的第一个时序运算单元；假设一次预测过程中t时刻的输入为x
t
，将x
t
输入vetmrrn的时序运算单元中得到该时刻的隐层状态h
t(1)
＝[h
1(1)
,h
2(1)
,h
3(1)
,
…
,h
l(1)
]，l表示当前隐层状态的维度；然后将该隐层状态h
t(1)
和下一网络层上一时刻的隐层状态h
t-1(2)
作为时序运算单元的输入，来计算下一隐层状态由于网络输入维度总是大于隐层状态维度，因此每一层隐层维度都将小于上一层隐层维度，所以最终隐层维度总能减少至1，此时该维度为1的最终隐层状态就作为vetmrrn的输出y
n-k+1
，也就是当前输入时间序列x的预测值；n-k+1时刻第l
层，即最终隐层的状态由n-k时刻该层的隐层状态和n-k+1时刻第l-1层的隐层状态共同决定，即：
[0017][0018]
式中，u
h(l)
表示第l层网络的输入权重矩阵，w
h(l)
表示第l层网络的隐层权重矩阵，b
h(l)
表示第l层网络的偏置，f为时序运算单元输入输出映射函数；
[0019]
n-k+1时刻第l层的隐藏层状态以计算为：
[0020][0021]
⊙
表示求hadamard积，tanh()为激活函数；
[0022]
将和继续向前一网络层和前一时刻展开，将展开为关于和的表达式，展开为关于和的表达式；以此类推，得到最终隐层状态关于输入时间序列x＝[x1,x2,x3,
…
x
n-k+1
]和初始隐层状态h0的表示式。
[0023]
进一步地，所述步骤s2中，训练vetmrrn的超参数和网络参数的具体方法为：设m时刻输入的空间滚动轴承性能退化特征序列为x＝[xm,x
m+1
,
…
,xn]，将该时间序列重构为训练样本集x和训练目标集y，表示如下：
[0024][0025]
将训练样本集x的每一行作为一个k-m+1维训练样本，该训练样本的训练目标为训练目标集y中对应的元素，故共有n-k个训练样本；
[0026]
vetmrrn的训练过程包括以下步骤：
[0027]
s21、初始化参数集ψ，即随机初始化vetmrrn每层网络的权重，并将所有偏置都初始化为0；同时随机初始化vetmrrn的超参数集
[0028]
s22、采用基于神谕门机制的元学习超参数自初始化网络对超参数集进行循环更新，得到vetmrrn的最优超参数；
[0029]
s23、将训练样本集x输入vetmrrn得到性能退化特征趋势预测值y'；
[0030]
s24、使用变分自编码元策略梯度学习算法对vetmrrn的参数集ψ进行更新，将更新后的参数集ψ作为下一次训练的初始值；
[0031]
s25、重复步骤s23和步骤s24的更新过程，直至关于训练目标集y和性能退化特征趋势预测值y'的均方差损失函数l(ψ)小于设定阈值，就完成了vetmrrn的训练。
[0032]
进一步地，所述步骤s22具体实现方法为：将待训练的vetmrrn的超参数组成超参数集其中l表示vetmrrn网络层数，k表示回顾序列长度，表示1到l层隐层维度；为寻找vetmrrn的最优超参数，构建元学习超参数自初始化网络来进行超参数更新，并设计神谕门机制来加快超参数自初始化过程；寻找最优超参数过程通过多
次训练任务来完成，一次训练任务过程具体如下：
[0033]
首先用本次任务中的超参数集的值来初始化待优化元学习超参数自初始化网络中的参数θ
t-1
；然后从原始输入时间序列x＝[x1,x2,
…
,xn]中随机截取一组连续时序子集x(
i)
＝[xm,x
m+1
,x
m+2
,
…
xn]作为超参数自初始化网络的输入y
t-1
并利用随机梯度下降法更新vetmrrn的权重和偏置；接下来将θ
t-1
、y
t-1
与超参数自初始化网络的隐层状态h
t
'-1
一同输入该超参数自初始化网络进行运算得到更新后的参数θ
t
和输出隐层状态h
t
'；随后，将参数θ
t
输入神谕门进行评估并获取损失函数，具体如下：
[0034]
定义一个酉算子o，其计算方法为：
[0035][0036]
式中，|θ
t
》为一个指标寄存器用来寄存参数θ
t
的量子位信息；表示模2加法，q为一个量子比特，h(q)为控制函数，当搜索到最优更新方向时令h(q)＝1翻转量子比特，否则令h(q)＝0不进行翻转；由于量子比特q始终保持不变，这里略去上式对量子比特的依赖，简写为：
[0037][0038]
然后施加神谕门变换，得到：
[0039][0040]
式中，表示haramard变换，表示条件相移，即对|0》以外的基本状态施加-1相位移动，|θ
t
》表示《θ
t
|的对偶向量，i表示单位矩阵；
[0041]
θ
t
输入神谕门得到的输出为g《θ
t
》＝2|θ
t
》《θ
t
|-i，于是将本次训练任务的损失函数构造为：
[0042][0043]
在当前时间步上的参数更新损失函数由当前的神谕门评估值g《θ
t
》和历史最佳评估值min
j＜t
[g《θj》]的差值来决定；如果当前任务的神谕门评估值与上一次任务比较没有得到改善，就认为当前更新方式对损失函数没有贡献，不更新超参数集否则对损失函数进行微分，求出对参数θ
t
的梯度，按该梯度方向对超参数集进行更新，即：
[0044][0045]
式中，表示学习率，表示梯度算子；
[0046]
至此，完成一次训练任务的vetmrrn超参数的更新过程，多次执行训练任务至神谕门评估值不再得到改善，认为vetmrrn的超参数已经达到最优值。
[0047]
进一步地，所述步骤s24具体实现方法为：针对vetmrrn的参数，即权重和偏置，设计变分自编码元策略梯度学习算法来对它们进行训练；
[0048]
将vetmrrn所有网络层的全部权重和偏置组成参数集并且定义vetmrrn的损失函数为预测值和
真实值的均方损失其中为t时刻的预测值，y
t
为t时刻的真实值；为更新该参数集ψ，将vetmrrn作为可与环境交互的智能体，将参数集ψ作为环境的状态值，于是定义参数集ψ的更新方式为：
[0049][0050]
式中，j(ψ)表示性能度量指标，用来度量参数集ψ的优劣；是一个随机估计，其期望近似于j(ψ)对于参数ψ
t
的梯度；α为学习率；
[0051]
将该随机估计实例化为对损失函数相反数的梯度，于是参数集ψ更新过程表示为：
[0052][0053]
将随机估计表示为时序差分误差的形式，于是更新方式转换成如下策略梯度形式：
[0054][0055]
式中，表示s
t
状态下执行动作a
t
得到的策略函数值，为时序差分误差，其中r
t+1
表示执行动作a
t
的收益，v(s
t+1
,w
t
)表示动作价值函数，w
t
表示动作价值函数的权重向量，γ∈[0,1]表示折扣系数；
[0056]
为使两种随机估计下得到的更新过程(2)和(3)一致，令r
t+1
＝0，于是式(2)进一步推导为：
[0057][0058]
然后，为加速参数更新，提高计算效率，在式(4)中构造资格迹算子来推广参数集ψ的更新方式；具体如下：
[0059]
将状态价值函数的权重w作为长时记忆向量，同时给定两个初始值为0的短时记忆向量即资格迹向量和并定义和的更新方式为：
[0060][0061][0062]
式中，λw和表示迹的衰减率，用来控制迹向量的权重，0＜λw＜1，ε为历史梯度贴现因子，初始化为1，并依据下式更新：ε
t+1
＝ε
t
·
γ；
[0063]
于是结合式(4)和式(5)，参数集ψ的更新方式进一步转换为：
[0064][0065]
类似地，根据式(6)和式(7)，权重w的更新方式定义为：
[0066][0067]
式中，α
ψ
和αw表示学习率；
[0068]
构造变分自编码器来编码n个元学习训练任务m＝(m1,m2,
…
,mn)用于训练参数集ψ；在每个训练任务mi中根据变分自编码器原理生成k个与时间序列x分布相似的虚拟时间序列将每个虚拟时间序列用于参数集ψ的更新，于是式(7)、式(8)分别表示为：
[0069][0070][0071]
式中，表示在任务mi中，将虚拟时间序列输入vetmrrn得到的t时刻的预测值；
[0072]
通过式(9)和(10)得到的参数集ψ
t+1j
来预测下一时刻的损失函数l(ψ
t+1j
)，并选择使该损失函数最小时的更新策略作为最优策略来更新参数集ψ，此时的虚拟时间序列为将该虚拟时间序列代入式(9)和(10)，得到：
[0073][0074][0075]
式(9)-式(12)完成了一个元学习训练任务对参数集ψ的更新过程，重复该更新过程直至参数集ψ不再更新，此时认为vetmrrn的参数集ψ已经达到全局最优解，也即变分自编码元策略梯度学习算法被视为训练至收敛。
[0076]
进一步地，所述步骤s3具体实现方法为：首先将作为测试样本输入vetmrrn，得到n+1时刻的预测值
[0077]
接着，利用该预测值构成新的测试样本并将该测试样本输入vetmrrn中以得到n+2时刻的预测值
[0078]
然后，用n+2时刻的预测值构成测试样本将该新测试样本输入vetmrrn中以得到n+3时刻的预测值
[0079]
以此类推，用n+n-1时刻的预测值构成测试样本并输入vetmrrn以得到n+n时刻的预测值最终得到n个性能退化特征预测值。
[0080]
进一步地，所述步骤s4具体实现方法为：结合训练用的性能退化特征序列x和预测得到的性能退化特征趋势值建立空间滚动轴承剩余寿命预测可靠度模型；
[0081]
首先，空间滚动轴承的基本故障率函数h(t,s)视为服从三参数威布尔分布，故该故障率函数计算如下：
[0082][0083]
式中，表示性能退化特征集合，χ为协变量回归参数，β表示形状参数，η表示尺度参数，t表示时间变量；
[0084]
然后，由于威布尔分布可靠度函数r(t,s)与故障率函数h(t,s)满足
则t时刻的可靠度表示为：
[0085][0086]
至此，完成威布尔分布可靠度模型的建立；假设起始预测时刻为ts，预测得到的空间滚动轴承精度失效阈值对应的时刻为t
fp
，则所预测的剩余寿命tr计算为:
[0087]
tr＝t
fp-ts+1
×
δt
[0088]
＝r-1
(r
fp
)-ns×
δt+1
×
δt
[0089]
＝(n
fp-ns+1)
×
δt
[0090]
式中，δt表示性能退化特征曲线中相邻两点之间的时间间隔，r
fp
表示所预测的精度失效阈值时间点对应的可靠度值，ns表示起始预测时间点，n
fp
表示所预测的精度失效阈值时间点。
[0091]
本发明的有益效果是：本发明在vetmrrn理论模型中，构建了新型时间循环网络结构以增大历史信息记忆量，从而提高了vetmrrn的长时记忆能力。另一方面，本发明设计了基于神谕门机制的元学习超参数自初始化网络来加速搜索vetmrrn的最优超参数以快速寻找vetmrrn的最优回顾序列长度，因此vetmrrn具有较好的非线性逼近能力、泛化性能和计算效率。此外，为提高对vetmrrn参数的训练速度和全局优化能力，设计了变分自编码元策略梯度学习算法，该算法构造了资格迹算子来加速参数训练速度(即进一步提高vetmrrn的计算效率)，并利用变分自编码器构建虚拟任务用于元策略梯度训练来寻找参数的全局最优解(即进一步提高vetmrrn的泛化性能)。由于vetmrrn在非线性逼近能力、泛化性能和计算效率等方面的优势，基于vetmrrn的空间滚动轴承剩余寿命预测方法可以实现较高的预测精度、较好的泛化性能和较高的计算效率。
附图说明
[0092]
图1为本发明的基于vetmrrn的空间滚动轴承剩余寿命预测方法的流程图；
[0093]
图2是vetmrrn的网络结构图；
[0094]
图3是元学习超参数自初始化网络结构图；
[0095]
图4是地面模拟空间环境下空间滚动轴承振动监测平台示意图；
[0096]
图5是本实施例中1号空间滚动轴承的shapely值融合特征；
[0097]
图6是vetmrrn预测的shapely值融合特征；
[0098]
图7是由威布尔分布可靠度模型得到的可靠度；
[0099]
图8是基于vetmrrn的方法预测的剩余寿命；
[0100]
图9是五种机器学习方法的学习曲线；
[0101]
图10是基于lstm的方法预测的剩余寿命，(a)是lstm预测的shapely值融合特征，(b)是基于lstm的方法预测的剩余寿命；
[0102]
图11是基于gru的方法预测的剩余寿命，(a)是gru预测的shapely值融合特征，(b)是基于gru的方法预测的剩余寿命；
[0103]
图12是基于dsae-nn的方法预测的剩余寿命，(a)是dsae-nn预测的shapely值融合特征，(b)是基于dsae-nn的方法预测的剩余寿命；
[0104]
图13是基于mk-lssvm的方法预测的剩余寿命，(a)是mk-lssvm预测的shapely值融合特征，(b)是基于mk-lssvm的方法预测的剩余寿命；
[0105]
图14是五种剩余寿命预测方法的计算时间对比。
具体实施方式
[0106]
元学习能寻找网络最优回顾序列长度，增强网络长时记忆能力，避免长时依赖缺陷从而有助于网络泛化性能的提高。强化学习能快速逼近最优学习策略，具有加速学习算法收敛的作用。结合时间循环网络、元学习和强化学习的各自优势，本发明设计了一种新型时间循环神经网络——变分资格迹元强化循环网络(variational eligibility trace meta-reinforce recurrent network，vetmrrn)。vetmrrn具有良好的非线性逼近能力、泛化性能及收敛速度。鉴于vetmrrn的以上优势，本发明专利进而发明了基于vetmrrn的空间滚动轴承剩余寿命预测方法：首先采用shapely值特征融合方法从空间滚动轴承的振动加速度数据中提取时、频域融合特征作为性能退化特征；然后将性能退化特征输入vetmrrn中预测空间滚动轴承的性能退化特征趋势；最后由性能退化特征趋势值建立威布尔分布可靠度模型来预测空间滚动轴承的精度失效阈值点和剩余寿命。下面结合附图进一步说明本发明的技术方案。
[0107]
如图1所示，本发明的基于vetmrrn的空间滚动轴承剩余寿命预测方法，包括以下步骤：
[0108]
s1、从空间滚动轴承的原始振动加速度数据中提取时、频域特征，进行shapely值特征融合，得到(1维)shapely值融合特征作为空间滚动轴承性能退化特征；提取的时、频域特征包括平均幅值、均方幅值、峰值、裕度、偏斜度、峭度、能量谱、倒频谱和功率谱。该shapely值融合特征反映了振动时、频域特征对时间的累积变化：空间滚动轴承正常状态下该融合特征保持平稳上升；而空间滚动轴承退化阶段的时、频域特征开始变得不再稳定，即与时间的相关性消失，于是该融合特征将保持不变。因此shapely值融合特征能定量描述空间滚动轴承振动时、频域特征的变化趋势，即能反映空间滚动轴承的性能退化趋势。
[0109]
s2、将得到的空间滚动轴承性能退化特征输入到vetmrrn中用于训练vetmrrn的超参数和网络参数；
[0110]
所述vetmrrn为变分资格迹元强化循环网络，在vetmrrn中，首先构建了新型时间循环网络结构以增大历史信息记忆量，并设计了基于神谕门机制的元学习超参数自初始化网络来加速搜索vetmrrn的最优超参数以快速寻找最优回顾序列长度，因此vetmrrn具有较好的非线性逼近能力、泛化性能和计算效率。另一方面，为提高对vetmrrn参数的训练速度和全局优化能力，设计了变分自编码元策略梯度学习算法，该算法构造了资格迹算子来加速参数训练速度，同时利用变分自编码器构建虚拟任务用于元策略梯度训练来寻找参数的全局最优解。vetmrrn包括l层网络结构，每层网络结构中包括n-k+1个时序运算单元，如图2所示；在vetmrrn中，首先引入更新单元和重置单元作为vetmrrn的基本时间序列运算单元；然后通过元学习的方式训练vetmrrn的待优化超参数；最后设计神谕门机制以加速vetmrrn的待优化超参数的自初始化过程。
[0111]
假设vetmrrn的初始隐层状态为h0；vetmrrn的输入时间序列为x＝[x1,x2,
…
,xn]，将该时间序列重构为矩阵形式：
[0112][0113]
式中，k为回顾序列长度，表示vetmrrn完成时间序列预测所需回顾的历史时间步长，也就是输入数据维度；将每个行向量视作当前时刻的输入时间序列，即x
t
表示一次预测过程中t时刻的输入时间序列，1≤t≤n-k+1；
[0114]
则更新单元t时刻的计算式如下：
[0115]zt
＝σ(x
t
uz+h
t-1
wz+bz)
ꢀꢀꢀꢀꢀ
(2)
[0116]
式中，z
t
为更新单元输出值，uz表示输入时间序列的权值矩阵，wz表示决定上一个时间步有多少遗留信息将要参与当前时间步运算的权值矩阵，bz表示偏置向量，σ表示sigmoid函数。
[0117]
重置单元t时刻的计算式为：
[0118]rt
＝σ(x
t
ur+h
t-1
wr+br)
ꢀꢀꢀꢀꢀ
(3)
[0119]
式中，r
t
为重置单元输出值。
[0120]
更新单元和重置单元共同组成vetmrrn的时序运算单元，于是可对vetmrrn的隐层变量更新如下：
[0121]ht
＝z
t
⊙ht-1
+(1-z
t
)
⊙
tanh(whx
t
+r
t
⊙
uhh
t-1
)
ꢀꢀꢀ
(4)
[0122]
式中，
⊙
表示hadamard积，tanh()为激活函数。由式(4)可知，当前隐层状态与上一时刻隐层状态相关，因此vetmrrn具有长时记忆的特点。
[0123]
x1,x2,
…
,x
n-k
,x
n-k+1
分别输入第一层网络结构的n-k+1个时序运算单元；初始化隐层状态h0＝[0,0,
…
,0]，并分别输入每层网络结构的第一个时序运算单元；假设一次预测过程中t时刻的输入为x
t
，将x
t
输入vetmrrn的时序运算单元中得到该时刻的隐层状态h
t(1)
＝[h
1(1)
,h
2(1)
,h
3(1)
,
…
,h
l(1)
]，l表示当前隐层状态的维度；然后将该隐层状态h
t(1)
和下一网络层上一时刻的隐层状态h
t-1(2)
作为时序运算单元的输入，来计算下一隐层状态由于网络输入维度总是大于隐层状态维度，因此每一层隐层维度都将小于上一层隐层维度，所以最终隐层维度总能减少至1，此时该维度为1的最终隐层状态就作为vetmrrn的输出y
n-k+1
，也就是当前输入时间序列x的预测值；n-k+1时刻第l层，即最终隐层的状态由n-k时刻该层的隐层状态和n-k+1时刻第l-1层的隐层状态共同决定，即：
[0124][0125]
式中，u
h(l)
表示第l层网络的输入权重矩阵，w
h(l)
表示第l层网络的隐层权重矩阵，b
h(l)
表示第l层网络的偏置，f为时序运算单元输入输出映射函数；
[0126]
n-k+1时刻第l层的隐藏层状态以计算为：
[0127][0128]
将将式(6)中的和继续向前一网络层和前一时刻展开，将展开为关于和的表达式，展开为关于和的表达式；以此类推，得到最终隐层状态关于输入时间序列x＝[x1,x2,x3,
…
x
n-k+1
]和初始隐层状态h0的表示式。
[0129]
由式(6)可知，隐层状态由前一个时间步的隐层状态和上一层网络的隐层状态共同决定，因而vetmrrn的网络结构能最大程度地记忆历史信息，并且随着网络层数的增加，历史信息记忆量随之增加。
[0130]
训练vetmrrn的超参数和网络参数的具体方法为：设m时刻输入的空间滚动轴承性能退化特征序列为x＝[xm,x
m+1
,
…
,xn]，将该时间序列重构为训练样本集x和训练目标集y，表示如下：
[0131][0132]
将训练样本集x的每一行作为一个k-m+1维训练样本，该训练样本的训练目标为训练目标集y中对应的元素，故共有n-k个训练样本；
[0133]
vetmrrn的训练过程包括以下步骤：
[0134]
s21、初始化参数集ψ，即随机初始化vetmrrn每层网络的权重，并将所有偏置都初始化为0；同时随机初始化vetmrrn的超参数集
[0135]
s22、采用基于神谕门机制的元学习超参数自初始化网络对超参数集进行循环更新，得到vetmrrn的最优超参数；具体实现方法为：将待训练的vetmrrn的超参数组成超参数集其中l表示vetmrrn网络层数，k表示回顾序列长度，表示1到l层隐层维度；为寻找vetmrrn的最优超参数，构建元学习超参数自初始化网络来进行超参数更新，并设计神谕门机制来加快超参数自初始化过程，该基于神谕门机制的元学习超参数自初始化网络结构如图3所示。寻找最优超参数过程通过多次训练任务来完成，一次训练任务过程具体如下：
[0136]
首先用本次任务中的超参数集的值来初始化待优化元学习超参数自初始化网络中的参数θ
t-1
；然后从原始输入时间序列x＝[x1,x2,
…
,xn]中随机截取一组连续时序子集x(i)＝[xm,x
m+1
,x
m+2
,
…
xn]作为超参数自初始化网络的输入y
t-1
并利用随机梯度下降法更新vetmrrn的权重和偏置；接下来将θ
t-1
、y
t-1
与超参数自初始化网络的隐层状态h’t-1
一同输入该超参数自初始化网络进行运算得到更新后的参数θ
t
和输出隐层状态h’t
；随后，将参数θ
t
输入神谕门进行评估并获取损失函数，具体如下：
[0137]
定义一个酉算子o，其计算方法为：
[0138]
[0139]
式中，|θ
t
》为一个指标寄存器用来寄存参数θ
t
的量子位信息；表示模2加法，q为一个量子比特，h(q)为控制函数，当搜索到最优更新方向时令h(q)＝1翻转量子比特，否则令h(q)＝0不进行翻转；由于量子比特q始终保持不变，这里略去式(7)对量子比特的依赖，简写为：
[0140][0141]
然后对式(8)的结果施加神谕门变换，得到：
[0142][0143]
式中，表示haramard变换，表示条件相移，即对|0》以外的基本状态施加-1相位移动，|θ
t
》表示《θ
t
|的对偶向量，i表示单位矩阵；
[0144]
由式(9)可知，θ
t
输入神谕门得到的输出为g《θ
t
》＝2|θ
t
》《θ
t
|-i，于是将本次训练任务的损失函数构造为：
[0145][0146]
由(10)可知，在当前时间步上的参数更新损失函数由当前的神谕门评估值g《θ
t
》和历史最佳评估值min
j＜t
[g《θj》]的差值来决定；如果当前任务的神谕门评估值与上一次任务比较没有得到改善，就认为当前更新方式对损失函数没有贡献，不更新超参数集否则对损失函数进行微分，求出对参数θ
t
的梯度，按该梯度方向对超参数集进行更新，即：
[0147][0148]
式中，表示学习率，表示梯度算子；
[0149]
至此，完成一次训练任务的vetmrrn超参数的更新过程，多次执行训练任务至神谕门评估值不再得到改善，认为vetmrrn的超参数已经达到最优值。由于基于神谕门机制的元学习超参数自初始化网络能让vetmrrn快速寻找到最优超参数(包括最优回顾序列长度)，因此vetmrrn能够更好地适应当前输入数据时序长度，从而避免了传统时间序列预测模型长时依赖的缺陷。
[0150]
s23、将训练样本集x输入vetmrrn得到性能退化特征趋势预测值y'；
[0151]
s24、使用变分自编码元策略梯度学习算法对vetmrrn的参数集ψ进行更新，将更新后的参数集ψ作为下一次训练的初始值；具体实现方法为：针对vetmrrn的参数，即权重和偏置，设计变分自编码元策略梯度学习算法来对它们进行训练；
[0152]
将vetmrrn所有网络层的全部权重和偏置组成参数集并且定义vetmrrn的损失函数为预测值和真实值的均方损失其中为t时刻的预测值，y
t
为t时刻的真实值；为更新该参数集ψ，将vetmrrn作为可与环境交互的智能体，将参数集ψ作为环境的状态值，于是定义参数集ψ的更新方式为：
[0153]
[0154]
式中，j(ψ)表示性能度量指标，用来度量参数集ψ的优劣；是一个随机估计，其期望近似于j(ψ)对于参数ψ
t
的梯度；α为学习率；
[0155]
将该随机估计实例化为对损失函数相反数的梯度，于是参数集ψ更新过程表示为：
[0156][0157]
将随机估计表示为时序差分误差的形式，于是式(12)的更新方式转换成如下策略梯度形式：
[0158][0159]
式中，表示s
t
状态下执行动作a
t
得到的策略函数值，
[0160]
为时序差分误差，其中r
t+1
表示执行动作a
t
的收益，v(s
t+1
,w
t
)表示动作价值函数，w
t
表示动作价值函数的权重向量，γ∈[0,1]表示折扣系数；
[0161]
比较式(13)和式(14)，为使两种随机估计下得到的更新过程(13)和(14)一致，令r
t+1
＝0，于是式(13)进一步推导为：
[0162][0163]
然后，为加速参数更新，提高计算效率，在式(15)中构造资格迹算子来推广参数集ψ的更新方式；具体如下：
[0164]
将状态价值函数的权重w作为长时记忆向量，同时给定两个初始值为0的短时记忆向量即资格迹向量和并定义和的更新方式为：
[0165][0166][0167]
式中，λw和表示迹的衰减率，用来控制迹向量的权重，0＜λw＜1，ε为历史梯度贴现因子，初始化为1，并依据下式更新：ε
t+1
＝ε
t
·
γ；
[0168]
于是结合式(15)和式(16)，参数集ψ的更新方式进一步转换为：
[0169][0170]
类似地，根据式(17)和式(18)，权重w的更新方式可定义为：
[0171][0172]
式中，α
ψ
和αw表示学习率；
[0173]
至此，完成资格迹算子的构造。由式(15)-(19)可知，资格迹向量追踪了参数集ψ的历史更新信息，因此能为ψ
t
的更新过程寻找到更优更新方向，从而进一步提高了变分自编码元策略梯度学习算法用于vetmrrn参数训练时的收敛速度。
[0174]
为提高对vetmrrn参数集ψ的全局优化能力，构造变分自编码器来编码n个元学习
训练任务m＝(m1,m2,
…
,mn)用于训练参数集ψ；在每个训练任务mi中根据变分自编码器原理生成k个与时间序列x分布相似的虚拟时间序列将每个虚拟时间序列用于参数集ψ的更新，于是式(18)、式(19)分别表示为：
[0175][0176][0177]
式中，表示在任务mi中，将虚拟时间序列输入vetmrrn得到的t时刻的预测值；
[0178]
通过式(20)和(21)得到的参数集ψ
t+1j
来预测下一时刻的损失函数l(ψ
t+1j
)，并选择使该损失函数最小时的更新策略作为最优策略来更新参数集ψ，此时的虚拟时间序列为将该虚拟时间序列代入式(20)和(21)，得到：
[0179][0180][0181]
式(20)-式(23)完成了一个元学习训练任务对参数集ψ的更新过程，重复该更新过程直至参数集ψ不再更新，此时认为vetmrrn的参数集ψ已经达到全局最优解，也即变分自编码元策略梯度学习算法被视为训练至收敛。显然，由于变分自编码器能够编码多个元学习训练任务用于参数集ψ的更新，因而参数集ψ将有更多可供选择的更新策略来指导该参数集收敛至全局最优解，因此，变分自编码元策略梯度学习算法具有对vetmrrn参数集良好的全局优化能力。
[0182]
s25、重复步骤s23和步骤s24的更新过程，直至关于训练目标集y和性能退化特征趋势预测值y'的均方差损失函数l(ψ)小于设定阈值，就完成了vetmrrn的训练。
[0183]
s3、利用训练好的vetmrrn来多步预测空间滚动轴承以后的性能退化特征趋势；
[0184]
具体实现方法为：首先将作为测试样本输入vetmrrn，得到n+1时刻的预测值
[0185]
接着，利用该预测值构成新的测试样本并将该测试样本输入vetmrrn中以得到n+2时刻的预测值
[0186]
然后，用n+2时刻的预测值构成测试样本将该新测试样本输入vetmrrn中以得到n+3时刻的预测值
[0187]
以此类推，用n+n-1时刻的预测值构成测试样本并输入vetmrrn以得到n+n时刻的预测值最终得到n个性能退化特征预测值。
[0188]
s4、结合训练用的空间滚动轴承性能退化特征和预测得到的性能退化特征趋势，建立威布尔分布可靠度模型，预测出空间滚动轴承的精度失效阈值时间点和剩余寿命；
[0189]
具体实现方法为：结合训练用的性能退化特征序列x和预测得到的性能退化特征趋势值建立空间滚动轴承剩余寿命预测可靠度模型；
[0190]
首先，空间滚动轴承的基本故障率函数h(t,s)视为服从三参数威布尔分布，故该故障率函数计算如下：
[0191][0192]
式中，表示性能退化特征集合，χ为协变量回归参数，β表示形状参数，η表示尺度参数，t表示时间变量；
[0193]
然后，由于威布尔分布可靠度函数r(t,s)与故障率函数h(t,s)满足则t时刻的可靠度表示为：
[0194][0195]
至此，完成威布尔分布可靠度模型的建立；假设起始预测时刻为ts，预测得到的空间滚动轴承精度失效阈值对应的时刻(即预测的可靠度为90％所对应的时刻)为t
fp
，则所预测的剩余寿命tr计算为:
[0196][0197]
式中，δt表示性能退化特征曲线中相邻两点之间的时间间隔，r
fp
表示所预测的精度失效阈值时间点对应的可靠度值，ns表示起始预测时间点，n
fp
表示所预测的精度失效阈值时间点；式中，r是r(t,s)的简写；由于r(t,s)中的表示性能退化特征集合，是个给定值，因此r(t,s)是关于自变量时间t的函数；r-1
(r
fp
)表示预测的精度失效阈值时间点对应的可靠度函数对应的反函数，该反函数的输出值是时间参数，所以r-1
(r
fp
)的输出值就是时间n
fp
×
δt。
[0198]
下面通过具体实验进一步验证本发明的技术效果。
[0199]
实验平台：利用自行搭建的地面模拟空间状态下空间滚动轴承振动监测(即寿命试验)平台采集的振动加速度数据，来验证所发明的基于vetmrrn的空间滚动轴承剩余寿命预测方法的有效性。该实验平台利用真空泵来模拟空间环境，轴承振动监测试验台和传感器放置如图4所示，将型号为c36018的空间滚动轴承安装在轴承振动监测试验台上。该试验台转轴转速为1000r/min，并通过弹簧机构对轴承和转轴施加7kn的轴向载荷。使用高灵敏度的icp加速度计对振动加速度进行采集，采样频率为25.6khz，采样长度为20480个点，采样间隔为2h。
[0200]
在转速为1000r/min条件下，1号空间滚动轴承连续运行1488h后完全止动失效，利用这期间采集的全寿命样本提取性能退化特征(即shapely值融合特征)，如图5所示。从起始点1到111点处，shapely值融合特征不断上升，表明该空间轴承处于正常运行阶段；从第112点(即利用2.4节的威布尔分布可靠度模型求出的可靠度为90％所对应的实际精度失效
阈值时间点nf＝112)开始，shapely值融合特征基本停止上升并保持在较高的数值，这表明该空间滚动轴承处于初始退化阶段，即早期故障阶段，此时可以认为该空间滚动轴承已经发生精度失效。
[0201]
对于空间滚动轴承而言，从正常阶段到初始退化阶段前期是剩余寿命预测中最关键的时间区间。将该区间内前84个shapely值融合特征作为训练样本，将第ns＝85个及其后续65个shapely值融合特征(即从正常阶段后期到初始退化阶段前期的特征数据)作为测试样本。对vetmrrn中的参数设置如下：在vetmrrn的网络结构中，由1.1.2节的元学习超参数自初始化网络来训练获得vetmrrn的超参数；在变分自编码元策略梯度学习算法中，折现因子γ＝0.08，学习率α＝0.05，迹衰减率λw＝0.8和均方误差阈值mse
min
＝5
×
10-4
。完成vetmrrn的参数设置后，就可依据2.2节和2.3节的vetmrrn训练和预测方法得到性能退化特征趋势预测结果，如图6所示。显然，在时间点区间[85,150]内，vetmrrn预测的性能退化特征值与实际值非常接近，二者变化趋势较一致。
[0202]
然后，由威布尔分布可靠度模型，利用作为训练样本的shapely值融合特征序列和vetmrrn预测的shapely值融合特征序列可得到1号空间滚动轴承的可靠度预测曲线，如图7所示。其中，可靠度为90％所对应的时间点，即精度失效阈值时间点预测值为n
fp
＝110.1。
[0203]
因此，所要预测的1号空间滚动轴承的剩余寿命就是预测的精度失效阈值点n
fp
＝110.1和起始预测点ns＝85之间的时间差值，即tr＝(n
fp-ns+1)
×
δt＝52.20h；同时，实际的剩余寿命为t
ra
＝(n
f-ns+1)
×
δt＝56.00h，其中δt＝2h为图8所示性能退化特征曲线中相邻两点的时间间隔。
[0204]
(1)为验证所发明的基于vetmrrn的空间滚动轴承剩余寿命预测方法的优势，首先，将vetmrrn与lstm、门控循环单元(gru)、深度稀疏自编码神经网络(dsae-nn)、多核最小二乘支持向量机(mk-lssvm)的非线性逼近能力及收敛速度进行比较。为了在相同条件下比较这五种机器学习方法的性能，后四种方法的性能退化特征也用shapely值融合特征表示，同时将后四种方法的最大训练步长和mse阈值分别设置为与vetmrrn相同的值(在下文的对比分析中保持相同的做法)。同时，将lstm和gru的学习率都设为η＝0.05；dsae-nn的分层神经元个数为n＝5，稀疏性惩罚因子的权重β＝0.3和稀疏性参数ρ＝0.02；mk-lssvm的惯性权重最大值和最小值分别设为ω
max
＝0.9和ω
min
＝0.4，多项式核阶数设为d＝3，其余参数通过粒子群优化算法(pso)优化选取。以上所有参数设置好后就保持不变。于是得到五种机器学习方法的学习曲线如图9所示。可以看出，在第10步之前，vetmrrn的mse相比其他方法能更快地下降到0.25左右，在后续的训练过程中，vetmrrn的mse相较其他四种方法也明显较小，这表明vetmrrn具有更好的非线性逼近能力和收敛速度。
[0205]
(2)比较所发明的基于vetmrrn的剩余寿命预测方法与基于lstm、gru、dsae-nn和mk-lssvm的预测方法的预测精度。为保证条件相同，后四种方法采用与vetmrrn相同的可靠度模型(即2.4节的威布尔分布可靠度模型)。后四种方法分别预测得到的1号空间滚动轴承的shapely值融合特征和剩余寿命分别如图10-13所示。图10是基于lstm的方法预测的剩余寿命，其中，(a)是lstm预测的shapely值融合特征，(b)是基于lstm的方法预测的剩余寿命。图11是基于gru的方法预测的剩余寿命，其中，(a)是gru预测的shapely值融合特征，(b)是基于gru的方法预测的剩余寿命。图12是基于dsae-nn的方法预测的剩余寿命，其中，(a)是dsae-nn预测的shapely值融合特征，(b)是基于dsae-nn的方法预测的剩余寿命。图13是基
于mk-lssvm的方法预测的剩余寿命，其中，(a)是mk-lssvm预测的shapely值融合特征，(b)是基于mk-lssvm的方法预测的剩余寿命。
[0206]
图6、图10(a)、图11(a)、图12(a)和图13(a)表明vetmrrn预测的shapely值融合特征与lstm、gru、dsae-nn和mk-lssvm相比最接近实际的shapely值融合特征，故依据vetmrrn预测的shapely值融合特征所建立的可靠度模型也更为准确、可靠。由图8、图10(b)、图11(b)、图12(b)和图13(b)也可看到，基于vetmrrn的剩余寿命预测方法和后四种方法预测的剩余寿命分别为52.20h、40.68h、46.18h、49.82h和60.34h，显然，基于vetmrrn的方法预测的剩余寿命与实际值(即56.00h)的误差最小，这表明vetmrrn具有更高的剩余寿命预测精度。
[0207]
(3)为定量评价所发明的vetmrrn、lstm、gru、dsae-nn和mk-lssvm这五种预测方法的泛化性能，以均方误差mse、平均绝对误差mae、平均绝对百分误差mape和拟合度r2作为泛化性能评价指标，即：
[0208][0209][0210][0211][0212]
式中，n为被预测样本的数量，yi′
和yi分别是被预测的shapely值融合特征及其对应的实际值，和var分别是n个实际shapely值融合特征的均值与方差。r2∈(0,1)越接近1，说明拟合程度越高，即泛化性能越好。为减少随机偶然性误差，对上述五种预测方法的泛化性能评价指标分别进行了统计分析，得到在50次重复预测中的泛化性能评价指标平均值，即和如表1所示。
[0213]
表1五种机器学习方法的泛化性能对比
[0214][0215]
由表1可知，相比于其他四种方法，vetmrrn的值最小，同时vetmrrn的值最接近1，这表明vetmrrn在多次重复实验中预测精度更高、泛化性能更好。
[0216]
(4)比较所发明的基于vetmrrn的剩余寿命预测方法的计算时间(即shapely值融合特征计算时间、vetmrrn训练时间、shapely值融合特征预测时间和剩余寿命计算时间之和)与使用lstm、gru、dsae-nn和mk-lssvm进行剩余寿命预测的计算时间。采用pycharm 2020.1.5在以下硬件配置环境下测量计算时间：16gram和3.2ghz intel cpu。记录每种预测方法执行50次的平均计算时间，如图14所示，图中从左到右依次为本发明的方法、基于lstm的方法、基于gru的方法、基于mk-lssvm和基于dsae-nn的方法的计算时间。其中，基于vetmrrn的预测方法的计算时间仅为7.02s，而基于lstm的方法、基于gru的方法、基于dsae-nn的方法和基于mk-lssvm的方法的计算时间分别为12.95s、11.07s、8.31s和12.06s。显然，所发明的基于vetmrrn的空间滚动轴承剩余寿命预测方法的计算时间要短于后四种方法，因而它具有更高的计算效率。
[0217]
本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

技术特征：
1.基于vetmrrn的空间滚动轴承剩余寿命预测方法，其特征在于，包括以下步骤：s1、从空间滚动轴承的原始振动加速度数据中提取时、频域特征，进行shapely值特征融合，得到shapely值融合特征作为空间滚动轴承性能退化特征；s2、将得到的空间滚动轴承性能退化特征输入到vetmrrn中用于训练vetmrrn的超参数和网络参数；s3、利用训练好的vetmrrn来多步预测空间滚动轴承以后的性能退化特征趋势；s4、结合训练用的空间滚动轴承性能退化特征和预测得到的性能退化特征趋势，建立威布尔分布可靠度模型，预测出空间滚动轴承的精度失效阈值时间点和剩余寿命。2.根据权利要求1所述的基于vetmrrn的空间滚动轴承剩余寿命预测方法，其特征在于，所述步骤s1中提取的时、频域特征包括平均幅值、均方幅值、峰值、裕度、偏斜度、峭度、能量谱、倒频谱和功率谱。3.根据权利要求1所述的基于vetmrrn的空间滚动轴承剩余寿命预测方法，其特征在于，步骤s2中所述vetmrrn为变分资格迹元强化循环网络，vetmrrn包括l层网络结构，每层网络结构中包括n-k+1个时序运算单元；假设vetmrrn的初始隐层状态为h0；vetmrrn的输入时间序列为x＝[x1,x2,
…
,x
n
]，将该时间序列重构为矩阵形式：式中，k为回顾序列长度，表示vetmrrn完成时间序列预测所需回顾的历史时间步长，也就是输入数据维度；将每个行向量视作当前时刻的输入时间序列，即x
t
表示一次预测过程中t时刻的输入时间序列，1≤t≤n-k+1；x1,x2,
…
,x
n-k
,x
n-k+1
分别输入第一层网络结构的n-k+1个时序运算单元；初始化隐层状态h0＝[0,0,
…
,0]，并分别输入每层网络结构的第一个时序运算单元；假设一次预测过程中t时刻的输入为x
t
，将x
t
输入vetmrrn的时序运算单元中得到该时刻的隐层状态h
t(1)
＝[h
1(1)
,h
2(1)
,h
3(1)
,
…
,h
l(1)
]，l表示当前隐层状态的维度；然后将该隐层状态h
t(1)
和下一网络层上一时刻的隐层状态h
t-1(2)
作为时序运算单元的输入，来计算下一隐层状态由于网络输入维度总是大于隐层状态维度，因此每一层隐层维度都将小于上一层隐层维度，所以最终隐层维度总能减少至1，此时该维度为1的最终隐层状态就作为vetmrrn的输出y
n-k+1
，也就是当前输入时间序列x的预测值；n-k+1时刻第l层，即最终隐层的状态由n-k时刻该层的隐层状态和n-k+1时刻第l-1层的隐层状态共同决定，即：式中，u
h(l)
表示第l层网络的输入权重矩阵，w
h(l)
表示第l层网络的隐层权重矩阵，b
h(l)
表示第l层网络的偏置，f为时序运算单元输入输出映射函数；
n-k+1时刻第l层的隐藏层状态以计算为：
⊙
表示求hadamard积，tanh()为激活函数；将和继续向前一网络层和前一时刻展开，将展开为关于和的表达式，展开为关于和的表达式；以此类推，得到最终隐层状态关于输入时间序列x＝[x1,x2,x3,
…
x
n-k+1
]和初始隐层状态h0的表示式。4.根据权利要求3所述的基于vetmrrn的空间滚动轴承剩余寿命预测方法，其特征在于，所述步骤s2中，训练vetmrrn的超参数和网络参数的具体方法为：设m时刻输入的空间滚动轴承性能退化特征序列为x＝[x
m
,x
m+1
,
…
,x
n
]，将该时间序列重构为训练样本集x和训练目标集y，表示如下：将训练样本集x的每一行作为一个k-m+1维训练样本，该训练样本的训练目标为训练目标集y中对应的元素，故共有n-k个训练样本；vetmrrn的训练过程包括以下步骤：s21、初始化参数集ψ，即随机初始化vetmrrn每层网络的权重，并将所有偏置都初始化为0；同时随机初始化vetmrrn的超参数集s22、采用基于神谕门机制的元学习超参数自初始化网络对超参数集进行循环更新，得到vetmrrn的最优超参数；s23、将训练样本集x输入vetmrrn得到性能退化特征趋势预测值y'；s24、使用变分自编码元策略梯度学习算法对vetmrrn的参数集ψ进行更新，将更新后的参数集ψ作为下一次训练的初始值；s25、重复步骤s23和步骤s24的更新过程，直至关于训练目标集y和性能退化特征趋势预测值y'的均方差损失函数l(ψ)小于设定阈值，就完成了vetmrrn的训练。5.根据权利要求4所述的基于vetmrrn的空间滚动轴承剩余寿命预测方法，其特征在于，所述步骤s22具体实现方法为：将待训练的vetmrrn的超参数组成超参数集其中l表示vetmrrn网络层数，k表示回顾序列长度，表示1到l层隐层维度；为寻找vetmrrn的最优超参数，构建元学习超参数自初始化网络来进行超参数更新，并设计神谕门机制来加快超参数自初始化过程；寻找最优超参数过程通过多次训练任务来完成，一次训练任务过程具体如下：首先用本次任务中的超参数集的值来初始化待优化元学习超参数自初始化网络中的参数θ
t-1
；然后从原始输入时间序列x＝[x1,x2,
…
,x
n
]中随机截取一组连续时序子集x
(i)
＝[x
m
,x
m+1
,x
m+2
,
…
x
n
]作为超参数自初始化网络的输入y
t-1
并利用随机梯度下降法更新
vetmrrn的权重和偏置；接下来将θ
t-1
、y
t-1
与超参数自初始化网络的隐层状态h'
t-1
一同输入该超参数自初始化网络进行运算得到更新后的参数θ
t
和输出隐层状态h
t
'；随后，将参数θ
t
输入神谕门进行评估并获取损失函数，具体如下：定义一个酉算子o，其计算方法为：式中，|θ
t
>为一个指标寄存器用来寄存参数θ
t
的量子位信息；表示模2加法，q为一个量子比特，h(q)为控制函数，当搜索到最优更新方向时令h(q)＝1翻转量子比特，否则令h(q)＝0不进行翻转；由于量子比特q始终保持不变，这里略去上式对量子比特的依赖，简写为：然后施加神谕门变换，得到：式中，表示haramard变换，表示条件相移，即对|0>以外的基本状态施加-1相位移动，|θ
t
>表示<θ
t
|的对偶向量，i表示单位矩阵；θ
t
输入神谕门得到的输出为g<θ
t
>＝2|θ
t
><θ
t
|-i，于是将本次训练任务的损失函数构造为：在当前时间步上的参数更新损失函数由当前的神谕门评估值g<θ
t
>和历史最佳评估值min
j＜t
[g<θ
j
>]的差值来决定；如果当前任务的神谕门评估值与上一次任务比较没有得到改善，就认为当前更新方式对损失函数没有贡献，不更新超参数集否则对损失函数进行微分，求出对参数θ
t
的梯度，按该梯度方向对超参数集进行更新，即：式中，表示学习率，表示梯度算子;至此，完成一次训练任务的vetmrrn超参数的更新过程，多次执行训练任务至神谕门评估值不再得到改善，认为vetmrrn的超参数已经达到最优值。6.根据权利要求4所述的基于vetmrrn的空间滚动轴承剩余寿命预测方法，其特征在于，所述步骤s24具体实现方法为：针对vetmrrn的参数，即权重和偏置，设计变分自编码元策略梯度学习算法来对它们进行训练；将vetmrrn所有网络层的全部权重和偏置组成参数集并且定义vetmrrn的损失函数为预测值和真实值的均方损失其中为t时刻的预测值，y
t
为t时刻的真实值；为更新该参数集ψ，将vetmrrn作为可与环境交互的智能体，将参数集ψ作为环境的状态
值，于是定义参数集ψ的更新方式为：式中，j(ψ)表示性能度量指标，用来度量参数集ψ的优劣；是一个随机估计，其期望近似于j(ψ)对于参数ψ
t
的梯度；α为学习率；将该随机估计实例化为对损失函数相反数的梯度，于是参数集ψ更新过程表示为：将随机估计表示为时序差分误差的形式，于是更新方式转换成如下策略梯度形式：式中，π
ψt
(a
t
|s
t
)表示s
t
状态下执行动作a
t
得到的策略函数值，为时序差分误差，其中r
t+1
表示执行动作a
t
的收益，v(s
t+1
,w
t
)表示动作价值函数，w
t
表示动作价值函数的权重向量，γ∈[0,1]表示折扣系数；为使两种随机估计下得到的更新过程(2)和(3)一致，令r
t+1
＝0，于是式(2)进一步推导为：然后，为加速参数更新，提高计算效率，在式(4)中构造资格迹算子来推广参数集ψ的更新方式；具体如下：将状态价值函数的权重w作为长时记忆向量，同时给定两个初始值为0的短时记忆向量即资格迹向量和并定义和的更新方式为：的更新方式为：式中，λ
w
和表示迹的衰减率，用来控制迹向量的权重，0＜λ
w
＜1，ε为历史梯度贴现因子，初始化为1，并依据下式更新：ε
t+1
＝ε
t
·
γ；于是结合式(4)和式(5)，参数集ψ的更新方式进一步转换为：类似地，根据式(6)和式(7)，权重w的更新方式定义为：式中，α
ψ
和α
w
表示学习率；构造变分自编码器来编码n个元学习训练任务m＝(m1,m2,
…
,m
n
)用于训练参数集ψ；在每个训练任务m
i
中根据变分自编码器原理生成k个与时间序列x分布相似的虚拟时间序列将每个虚拟时间序列用于参数集ψ的更新，于是式(7)、式(8)分别表示为：
式中，表示在任务m
i
中，将虚拟时间序列输入vetmrrn得到的t时刻的预测值；通过式(9)和(10)得到的参数集来预测下一时刻的损失函数并选择使该损失函数最小时的更新策略作为最优策略来更新参数集ψ，此时的虚拟时间序列为将该虚拟时间序列代入式(9)和(10)，得到：将该虚拟时间序列代入式(9)和(10)，得到：式(9)-式(12)完成了一个元学习训练任务对参数集ψ的更新过程，重复该更新过程直至参数集ψ不再更新，此时认为vetmrrn的参数集ψ已经达到全局最优解，也即变分自编码元策略梯度学习算法被视为训练至收敛。7.根据权利要求1所述的基于vetmrrn的空间滚动轴承剩余寿命预测方法，其特征在于，所述步骤s3具体实现方法为：首先将作为测试样本输入vetmrrn，得到n+1时刻的预测值接着，利用该预测值构成新的测试样本并将该测试样本输入vetmrrn中以得到n+2时刻的预测值然后，用n+2时刻的预测值构成测试样本将该新测试样本输入vetmrrn中以得到n+3时刻的预测值以此类推，用n+n-1时刻的预测值构成测试样本并输入vetmrrn以得到n+n时刻的预测值最终得到n个性能退化特征预测值。8.根据权利要求1所述的基于vetmrrn的空间滚动轴承剩余寿命预测方法，其特征在于，所述步骤s4具体实现方法为：结合训练用的性能退化特征序列x和预测得到的性能退化特征趋势值建立空间滚动轴承剩余寿命预测可靠度模型；首先，空间滚动轴承的基本故障率函数h(t,s)视为服从三参数威布尔分布，故该故障率函数计算如下：式中，表示性能退化特征集合，χ为协变量回归参数，β表示形状参数，η表示尺度参数，t表示时间变量；然后，由于威布尔分布可靠度函数r(t,s)与故障率函数h(t,s)满足
则t时刻的可靠度表示为：至此，完成威布尔分布可靠度模型的建立；假设起始预测时刻为t
s
，预测得到的空间滚动轴承精度失效阈值对应的时刻为t
fp
，则所预测的剩余寿命t
r
计算为：式中，δt表示性能退化特征曲线中相邻两点之间的时间间隔，r
fp
表示所预测的精度失效阈值时间点对应的可靠度值，n
s
表示起始预测时间点，n
fp
表示所预测的精度失效阈值时间点。

技术总结
本发明公开了一种基于VETMRRN的空间滚动轴承剩余寿命预测方法，包括以下步骤：S1、从空间滚动轴承的原始振动加速度数据中提取时、频域特征，进行shapely值特征融合，作为空间滚动轴承性能退化特征；S2、将空间滚动轴承性能退化特征输入VETMRRN中训练VETMRRN的超参数和网络参数；S3、利用VETMRRN来多步预测空间滚动轴承以后的性能退化特征趋势；S4、建立威布尔分布可靠度模型，预测出空间滚动轴承的精度失效阈值时间点和剩余寿命。本发明构建了VETMRRN，其具有较好的非线性逼近能力、泛化性能和计算效率，使得基于VETMRRN的空间滚动轴承剩余寿命预测方法具有较高的预测精度、较好的泛化性能和较高的计算效率。的泛化性能和较高的计算效率。的泛化性能和较高的计算效率。


技术研发人员：李锋 姜沛轩 汪永超
受保护的技术使用者：四川大学
技术研发日：2021.12.13
技术公布日：2022/3/8