基于相对密度核估计的齿轮剩余寿命预测方法与流程

1.本发明属于大数据和智能制造技术领域，具体涉及一种基于相对密度核估计的齿轮剩余寿命预测方法。


背景技术：

2.近年来，随着工业化设备不断地朝着智能化、高性能和复杂化的方向发展。这类设备在受到系统内部因素和外部环境因素的影响下设备发生故障的概率也随之增加，设备故障会导致其性能和健康状态产生一定程度的退化。如果设备出现大的故障，这不仅造成财政和资源损失重大，也给人民的安全带来危险。为此，对设备进行实时监测和剩余寿命预测有着十分重要的研究价值和意义。
3.现研究阶段可将剩余寿命预测的方法大致分为基于物理失效模型、基于知识表示和数据驱动的方法。面对复杂化的设备，获得物理模型是非常困难的事，知识表示的方法更适合于定性推理而不太适合定量计算，且难以获得完整的知识。因此，数据驱动的方法成为剩余寿命预测研究的主流方向。
4.si等将数据驱动方法的伽马分布、回归模型、维纳过程和随机滤波模型等寿命预测方法进行总结分析；zhai等提出自适应维纳过程的剩余寿命预测模型；刘文溢等提出基于高阶隐半马尔科夫模型的剩余寿命预测模型。上述数据驱动方法大部分需要假设退化模型和参数估计，并且参数估计法对模型选择有局限性，过分依赖概率密度函数形式的先验界定，所以不能确保预测模型的精确性和适用性。
5.随着工业系统复杂化的发展，机器学习作为数据驱动的方法，裴洪等对基于机器学习的预测方法从浅层和深度学习两方面进行了详细分析总结。其中，基于深度学习的预测方法有强大的特征提取能力，不需要事先对未知的退化模型假设，是目前剩余寿命预测应用上的主流方法。zhang等提出用lstm神经网络的方法对设备进行预测；li等利用深度卷积神经网络方法进行寿命预测；而chen等结合了guo等和li等两种方法的优势提出了一种基于端到端可训练卷积递归神经网络的机械健康指标构建方法；张继冬等提出一种基于全卷积变分自编码网络的轴承剩余寿命预测方法。上述机器学习的方法，它的算法在适应制造系统和过程的复杂和非线性特性方面往往受到限制，而且其模型内部结构称为“黑盒子”，不能够清楚地表征系统退化特征的变化，且网络在学习过程中随着输入的不断增加，参数调节不能保证全局最优，从而使系统的剩余寿命预测准确性受到影响。
6.核密度估计的方法是对数据分布不附加任何假设，从数据本身出发研究数据分布特点的非参数估计方法。该方法避免了大部分数据驱动方法需要模型假设和参数估计的问题以及避免了机器学习的不足。因此，核密度估计方法在剩余寿命预测技术上的应用受到学者们的高度关注和重视。现有的核密度估计模型中， hu等用非参数核估计的方法对风速建模评估系统可靠性；sidebe等针对不同运行环境对随机系统退化状态的作用，提出用两个函数建模的方法对系统可靠性进行了计算；杨楠等提出一种基于非参数核密度估计的风功率波动性概率密度建模方法，并针对模型带宽选择问题，构造一种以拟合优度检验为约
束条件的带约束带宽优化模型；李存华等将核密度估计方法应用在聚类算法的构造上，提出基于网格数据重心的分箱核估计近似方法。上述文献中核估计窗宽的选取大部分采用的是固定的窗口宽度，将固定值作为窗宽会造成样本点分散的区域拟合度低，密集的区域拟合度过高。
7.为了解决固定窗宽的不足，赵渊等提出一种非参数多变量核密度估计负荷模型研究的方法，该方法在核估计窗宽选择上实现了自适应选取窗宽；颜伟等针对核估计中最优带宽选取的重要性，提出一种不依赖总体真实分布的最优带宽改进模型；为了提高预测的准确性，张卫贞等提出的实时剩余寿命预测方法，将积分均方误差方法引入核估计窗宽的选择上实现了自适应窗宽选取。上述方法能够自适应的选择窗宽和实现了特征退化分布和剩余寿命预测的实时更新，但是在数据分布不均匀的、变密度的样本处选择窗宽的合理性有待提高。而且现有核密度估计模型在有界支持[0, ∞]上的核估计量在0附近是有偏的，并且从支持的原点引入了一个值为-h的左移，会使得剩余寿命估计的准确度降低。


技术实现要素：

[0008]
针对核密度估计在设备剩余寿命预测时，由于在监测数据分布不均匀区域选择的窗宽不够准确和核估计模型中随机变量的有界性产生的边界偏差，从而导致预测不准确的技术问题，本发明提出了一种相对密度核估计的实时剩余寿命预测方法，预测的剩余寿命非常接近实际值，预测的准确性大大增加。
[0009]
为实现上述目的，本发明所采用的技术方案为：基于相对密度核估计的齿轮剩余寿命预测方法，具体步骤如下：
[0010]
一、建立估计未知变量的核密度估计模型；
[0011]
二、确定自适应的相对密度窗宽；
[0012]
三、计算齿轮退化样本的特征退化分布；
[0013]
四、建立齿轮剩余寿命预测模型；
[0014]
tn时刻系统剩余寿命预测的概率密度函数为：
[0015][0016]
式中，为剩余寿命预测的概率密度，表示t
n t
时刻的核微分同胚变换相对密度核估计的n t次卷积，δx为退化增量，为微分同胚变换；
[0017]
在实时变化的系统下获得新的样本数据后，通过可推算出下一时刻的进而实现齿轮实时剩余寿命预测。
[0018]
在步骤一中，核密度估计模型的具体建模过程如下：
[0019]
1.1、设δx1,δx2,
…
δxi,
…
,δxn为n个独立同分布的随机变量，f(δx)为其概率密度函数，则f(δx)的核密度估计为：
[0020][0021]
其中，h为窗宽，k(
·
)为核函数，n为样本数；
[0022]
在公式(1)中选用高斯核函数：
[0023][0024]
1.2、自适应窗宽采用公式(3)积分均方误差求其最小值得到初始最优窗宽 hn：
[0025][0026]
其中，δx为退化增量，为n个初始退化增量样本估计的概率密度函数， fn(δx)为n个初始退化增量样本实际的概率密度函数；
[0027]
由式(3)求得n个初始样本确定的窗宽hn为：
[0028][0029]
将高斯核函数代入式(4)可求出hn为：
[0030][0031]
式中，f
n”(δx)表示对fn(δx)求二次导函数,σn为n个初始样本特征退化增量的方差；
[0032]
在步骤二中，确定自适应相对密度窗宽的过程为：
[0033]
将k近邻思想计算样本点的相对密度引入自适应窗宽来提高核密度估计的准确性；
[0034]
2.1、建立相对密度模型：
[0035]
假设x1,x2...,xi,...xj,...,xn为n个样本点并用数据集a表示，则相对密度的模型建立过程如下：
[0036]
2.1.1、计算样本点xi与样本点xj的欧氏距离d(xi,xj)：
[0037][0038]
2.1.2、计算样本点xi的k近邻距离：k_dist(xi)＝d(xi,xj)，并且满足：
[0039]
a)、对于任意正整数k，在样本中至少有不包括xi在内的k个点x
′j∈a\{xi}，则d(xi,x
′j)≤d(xi,xj)，其中a\{xi}表示集合a中不包含样本点xi；
[0040]
b)、样本中最多有不包括xi在内的k-1个点x
′j∈a\{xi}，满足 d(xi,x
′j)＜d(xi,xj)；
[0041]
2.1.3、已知k_dist(xi)，样本点xi的k距离邻域可表示为：
[0042]
nk(xi)＝{x
′
∈a\{xi}|d(xi,x
′i)≤k_dist(xi)}
ꢀꢀꢀ
(7)
[0043]
式中，nk(xi)表示样本点xi的k距离邻域包含到xi的距离不大于k_dist(xi) 的所有样本，x
′i被称为xi的k近邻，d(xi,x
′i)是xi和x
′i的欧氏距离；
[0044]
2.1.4、计算样本点xi相对于xj的可达距离reach_dist(xi,xj)：
[0045]
reach_dist(xi,xj)＝max{k_dist(xj),d(xi,xj)}
ꢀꢀꢀ
(8)
[0046]
式中，d(xi,xj)为样本点xi与xj之间的欧氏距离，k_dist(xj)为样本点xj的 k近邻
距离，max{
·
,
·
}表示取极大值；
[0047]
2.1.5、计算样本点xi的局部可达密度，可表示为：
[0048][0049]
式中，lrd(xi)表示样本点xi的k距离邻域内点到xi的平均可达距离的倒数， |nk(xi)|表示nk(xi)的绝对值；
[0050]
2.1.6、样本点xi的相对密度ρ(xi)可表示为：
[0051][0052]
其中，lrd(xj)为样本点xj的局部可达密度，通过上述推导将公式(9)代入公式(10)可推出相对密度ρ(xi)的表达式为：
[0053][0054]
采用随机抽取的方式选取样本点xj；
[0055]
相对密度ρ(xi)表示的是样本点xi的局部可达密度与样本点xi的k距离邻域内的样本点nk(xi)的局部可达密度平均值之比，如果ρ(xi)越接近1，则点xi的邻域点密度相对均匀；若ρ(xi)越小于1，则点xi的密度高于其邻域点密度，xi为密集点；ρ(xi)越大于1，则点xi的密度小于其邻域点密度，xi为稀疏点；
[0056]
2.2、确定自适应相对密度窗宽：
[0057]
将相对密度ρ(xi)作为窗宽引入核密度估计的模型中，从而构建相对密度的核密度估计表达式如下：
[0058][0059]
其中，ρ(xi)为不同样本点处核估计的相对密度窗宽，且计算公式为式(11) 所示；k(
·
)为核函数；
[0060]
2.3、实时更新自适应相对密度窗宽：
[0061]
核密度估计的实时更新用已知的n个样本的核估计推导第n 1个样本的核密度估计，第n个样本点的核密度估计表示为：
[0062][0063]
第n 1个样本点的核密度估计为：
[0064][0065]
式中，ρ(x
n 1
)为样本点在x
n 1
处的相对密度窗宽；
[0066]
通过以上递推可知，当任意t
n j
时刻增加j个样本时，n j个样本数据的核密度估计可递推为：
[0067][0068]
在步骤三中，特征退化分布的计算过程如下：
[0069]
每隔单位时间采集一次退化数据，计算其退化增量δx1,δx2,
…
,δxn为n个独立同分布的退化增量样本，的核密度估计表示为：
[0070][0071]
其中，ρ(δxi)(i＝1,2,3,
…
,n)为不同样本点处的相对密度窗宽；
[0072]
核密度估计可得[0,tn]随机退化增量的概率密度推出[0,tn]上随机退化增量累计特征退化量xn(令xn＝nδx)的可由n重卷积得到：
[0073][0074]
当t
n 1
时刻增加一新样本时，[0,t
n 1
]的退化量的概率密度函数如下：
[0075][0076]
当t
n j
时刻新增j个样本时，[0,t
n j
]的退化量的概率密度函数为：
[0077][0078]
在步骤四中，建立剩余寿命预测模型的具体过程为：
[0079]
4.1、估计核微分同胚：
[0080]
设样本x1,x2,
…
,xn作为随机变量x恒等分布有概率密度函数f(x)，可以形成 fn的一个估计量序列，使只依赖于前n个观测值，设f为所有连续且有界的随机变量概率密度函数的集合，其支持范围限于实区间[a,b]，则核微分同胚估计量可表示为：
[0081][0082]
式中，k(
·
)是核函数；hn为带宽；是随机变量从区间[a,b]到的微分同胚变换；当x趋于下限或上限时，它的一阶导数趋于无穷，其中，和x的关系式如下：
[0083][0084]
4.2、建立核微分同胚估计的剩余寿命预测模型：
[0085]
对有界系统退化增量进行对数核微分同胚变换：
[0086][0087]
设tn为当前时刻，[0,tn]监测时间内采集到的当前样本退化数据，其相应的特征退化随时间退化的特征增量可以作为随机变量δx∈[0,x
th
](x
th
为失效阈值)；
[0088]
自适应相对密度的核微分同胚估计为：
[0089][0090]
k(
·
)是核函数，ρ(δxi)是相对密度窗宽；
[0091]
剩余寿命的概率密度分布函数f
t
(t)为：
[0092][0093]
g(x
n t
)为[0,t
n t
]特征退化量的概率密度；
[0094]
通过对已知的单位时间随机退化特征增量核微分同胚估计，用卷积求其特征退化量的概率密度函数，将对数微分同胚变换后[0,t
n t
]特征退化量的概率密度记为gd(x
n t
)：
[0095][0096]
剩余寿命预测的概率密度为：
[0097][0098]
式中，表示t
n t
时刻的核微分同胚变换相对密度核估计的n t次卷积，在实时监测数据的更新下，t
n t
时刻n t个样本核微分同胚变换的相对密度核估计为：
[0099][0100]
根据公式(25)和(27)可推出：
[0101][0102]
从而能够推出tn时刻系统剩余寿命预测的概率密度函数为：
[0103][0104]
在实时变化的系统下获得新的样本数据后，通过可推算出下一时刻的进而实现实时剩余寿命预测。
[0105]
本发明与现有技术相比，具体有益效果体现在：
[0106]
一、本发明该模型将k近邻思想计算出的样本点的相对密度作为核密度估计的窗宽，不仅解决了固定窗宽由样本数据分布不均导致拟合不足的问题，而且能够对任意形状、密度不均匀的数据集自适应地选择出更加准确的窗宽，提高了拟合度。
[0107]
二、本发明引入核微分同胚变换的方法，有效地消除了核密度估计带来的边界偏差问题。
[0108]
三、本发明还建立了基于相对密度核估计的实时更新模型，有效避免了核密度估计的重复计算，提高了计算效率。
[0109]
四、本发明随着样本数据的增加，剩余寿命的预测越来越接近实际值，方差变得越来越小，提高了预测的准确性和有效性。
附图说明
[0110]
图1为本发明的流程框图。
[0111]
图2为核密度估计示意图。
[0112]
图3为固定窗宽下的核密度估计示意图。
[0113]
图4为样本特征退化随时间变化的曲线图。
[0114]
图5为齿轮实验台架的结构图。
[0115]
图6为主齿轮箱的传感器位置示意图。
[0116]
图7为特征值随监测时间变化曲线示意图。
[0117]
图8为t-40h时刻样本点的相对密度窗宽选取示意图，其中，“o”是特征退化增量样本值，“*”是相对密度窗宽值。
[0118]
图9为t-70h时刻样本点的相对密度窗宽选取示意图，其中，“o”是特征退化增量样本值，“*”是相对密度窗宽值。
[0119]
图10为t-77h时刻样本点的相对密度窗宽选取示意图，其中，“o”是特征退化增量样本值，“*”是相对密度窗宽值。
[0120]
图11为不同监测时刻下，概率密度函数的预测曲线图，其中，“*”为估计值，“o”为实际值。
[0121]
图12为t＝40h的核估计概率密度比较图，其中，“*”为估计值，“o”为实际值。
[0122]
图13为t＝70h的核估计概率密度比较图，其中，“*”为估计值，“o”为实际值。
[0123]
图14为t＝40h的核估计偏差比较图，其中，左侧曲线为传统的核密度估计曲线，右侧曲线为核微分同胚的核密度估计曲线。
[0124]
图15为t＝70h的核估计偏差比较图，其中，左侧曲线为传统的核密度估计曲线，右侧曲线为核微分同胚的核密度估计曲线。
[0125]
图16为lstm模型下，周期数为420(10min)，即70h时预测的剩余寿命结果图。
[0126]
图17为lstm模型下，周期数为450(10min)，即75h时预测的剩余寿命结果图。
[0127]
图18为t＝60h时，本发明模型与gamma过程的比较图，“*”为预测值，“o”为实际值。
[0128]
图19为t＝70h时，本发明模型与gamma过程的比较图，“*”为预测值，“o”为实际值。
[0129]
图20为t＝60h时，本发明模型与wiener过程模型的比较图。
[0130]
图21为t＝70h时，本发明模型与wiener过程模型的比较图。
[0131]
图22为特征值随监测时间变化曲线图。
[0132]
图23为不同监测时刻的概率密度函数的预测比较图。
具体实施方式
[0133]
为了使本发明所要解决的技术问题、技术方案及有益效果更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
[0134]
如图1所示，基于相对密度核估计的齿轮剩余寿命预测方法，具体步骤如下：
[0135]
一、建立核密度估计模型：
[0136]
核密度估计方法对数据分布不附加任何假定，是一种从数据本身出发，用来估计未知变量的概率密度函数的方法。其由已知的n个样本点，通过选择任意核函数(如高斯核函数)及窗宽得到n个核函数，再线性叠加形成核密度的估计函数。例如，对样本集为：
[0137]
x＝{x1＝-2.1,x2＝-1.3,x3＝-0.4,x4＝1.9,x5＝5.1,x6＝6.2}的核密度估计拟合结果如图2所示。
[0138]
设δx1,δx2,
…
δxi,
…
,xn为n个独立同分布的随机变量，f(δx)为其概率密度函数，则f(δx)的核密度估计为：
[0139][0140]
式中，h为窗宽；k(
·
)为核函数；n为样本数；其中，核密度估计取决于k(
·
) 和h的选择。
[0141]
核函数作为影响核密度估计的一个因素，一般情况下任何函数都可作为核函数，常用的有四次核、均匀核、三角核和高斯核。核函数的选择对核密度估计的准确度作用不大，本发明选用广泛应用的高斯核函数。
[0142][0143]
窗宽h作为影响核密度估计平滑性和核函数宽度的主要因素，当h较小时，核密度估计曲线不够光滑且曲折，表露了较多细节；当窗宽较大时，核密度估计曲线比较平滑，但掩盖了较多细节。因此，选择合适的窗宽对核密度估计是非常重要的。现有的核估计窗宽分为固定窗宽和自适应窗宽，下面分别对两种窗宽进行介绍，并将其存在的问题进行分析。
[0144]
1、固定窗宽是在核密度估计时采用固定不变的值；这种选取窗宽的方法一般会导致核估计在低密度区拟合不足，高密度区过度拟合，通过图2可以清楚地分析出固定窗宽的不足之处。
[0145]
图3为固定窗宽h1＝12，在四组样本数据密度不同情况下的核密度估计示意图，由图3可知，在数据不足的区间，核估计欠拟合；在数据集中的区域，核估计过度拟合。
[0146]
2、目前自适应窗宽是窗宽选择的主流方向，其本质是随着样本点的增加能够自适应地选取窗宽，提高估计准确度。在现有研究阶段常用的自适应窗宽方法大多采用的是通过式(3)积分均方误差求其最小值得到初始最优窗宽hn：
[0147][0148]
其中，δx为退化增量，为n个初始样本估计的概率密度函数，fn(δx)为 n个初始样本实际的概率密度函数。
[0149]
由式(3)求得hn为：
[0150][0151]
式中：将高斯核函数代入式(4)可求出hn为：
[0152][0153]
式中，f
n”(δx)表示对fn(δx)求二次导函数,σn为n个初始样本特征退化增量的方差；
[0154]
上述自适应窗宽方法在实际应用中，解决了样本数据在实时变化下选取窗宽以及固定窗宽造成过(欠)拟合问题。如果样本接近于正态分布时，选取该方式是最优选择；但是当真实分布为非对称或者多峰时，该方法可能导致过度平滑，准确性有待提高。
[0155]
二、确定相对密度窗宽的确定：
[0156]
自适应窗宽的窗宽选择法已经成为一个发展趋势，该方法能够在监测数据实时变
化的情况下，自适应地选择合理的窗宽，也能解决固定窗宽带来数据拟合不足和过度的问题。而通过k近邻思想计算样本点的相对密度，其能通过样本点之间的k近邻距离范围判断出各样本点与其周围样本点的稀疏和密集程度，从而选择合理的相对密度。此外，相对密度的思想本质就是直接对数据分析计算来判断密度的大小，与核密度估计思想不谋而合，而且能够对任意形状的数据集进行快速、准确地识别样本的中心。因此，可以将k近邻思想计算样本点的相对密度引入自适应窗宽来提高核密度估计的准确性。
[0157]
2.1、相对密度的模型建立：
[0158]
假设x1,x2...,xi,...xj,...,xn为n个样本点并用数据集a表示，则相对密度的模型建立过程如下：
[0159]
2.1.1、计算样本点xi与样本点xj的欧氏距离d(xi,xj)：
[0160][0161]
2.2.2、计算样本点xi的k近邻距离k_dist(xi)：k_dist(xi)＝d(xi,xj)，并且满足：
[0162]
a)、对于任意正整数k，在样本中至少有不包括xi在内的k个点x
′j∈a\{xi}，则d(xi,x
′j)≤d(xi,xj)；
[0163]
b)、样本中最多有不包括xi在内的k-1个点x
′j∈a\{xi}，满足 d(xi,x'j)＜d(xi,xj)；
[0164]
2.1.3、已知k_dist(xi)，样本点xi的k距离邻域可表示为：
[0165]
nk(xi)＝{x
′i∈a\{xi}|d(xi,x
′i)≤k_dist(xi)}
ꢀꢀꢀ
(7)
[0166]
式中，nk(xi)表示样本点xi的k距离邻域包含到xi的距离不大于k_dist(xi) 的所有样本，x
′i被称为xi的k近邻，d(xi,x
′i)是xi和x
′i的殴氏距离。
[0167]
2.1.4、计算样本点xi相对于xj的可达距离reach_dist(xi,xj)：
[0168]
reach_dist(xi,xj)＝max{k_dist(xj),d(xi,xj)}
ꢀꢀꢀ
(8)
[0169]
式中，reach_dist(xi,xj)表示样本点xi相对于xj的可达距离；d(xi,xj)为样本点 xi与xj之间的欧氏距离，k_dist(xj)为样本点xj的k近邻距离，max{
·
,
·
}表示取极大值。
[0170]
2.1.5、计算样本点xi的局部可达密度，可表示为：
[0171][0172]
式中，lrd(xi)表示样本点xi的k距离邻域内点到xi的平均可达距离的倒数， |nk(xi)|为样本点xi的k距离邻域内所有样本，|nk(xi)|表示nk(xi)的绝对值。
[0173]
2.1.6、样本点xi的相对密度ρ(xi)可表示为：
[0174][0175]
其中，lrd(xj)为样本点xj的局部可达密度，通过上述推导将公式(9)代入公式(10)可推出相对密度ρ(xi)的表达式为：
[0176][0177]
本发明采用随机抽取的方法选取样本点xj。
[0178]
从公式不难看出相对密度ρ(xi)表示的是样本点xi的局部可达密度与样本点 xi的k距离邻域内的样本点|nk(xi)|的局部可达密度平均值之比。如果ρ(xi)越接近1，则点xi的邻域点密度相对均匀；若ρ(xi)越小于1，则点xi的密度高于其邻域点密度，xi为密集点；ρ(xi)越大于1，则点xi的密度小于其邻域点密度，xi为稀疏点。
[0179]
2.2、自适应相对密度窗宽的确定：
[0180]
样本点的相对密度ρ(xi)能够通过计算点之间的距离来计算样本的密度，点之间的距离越大，密度越低；点之间的距离越小，密度越高。而通过ρ(xi)来选择窗宽，能够根据ρ(xi)与数值1之间的大小关系来先判断样本的疏密程度再选取窗宽，从而在低密集区域选择大的窗宽，在高密度区域选择小的窗宽。因此，将相对密度ρ(xi)作为窗宽引入核密度估计的模型中，从而构建相对密度的核密度估计表达式如下：
[0181][0182]
其中，ρ(xi)为不同样本点处核估计的相对密度窗宽，且计算公式为式(11) 所示；k(
·
)为核函数。
[0183]
2.3、自适应相对密度窗宽的实时更新
[0184]
由于在实际应用中样本数据都是实时更新的，如果每增加一个样本都对其从头开始计算，那么计算量会随数据的增多变得复杂化。所以，为使核估计的计算性能得到提升，实现核密度估计的实时更新是不可或缺的。
[0185]
核密度估计的实时更新用已知的n个样本的核估计推导第n 1个样本的核密度估计，推导过程如下：
[0186]
第n个样本点的核密度估计表示为：
[0187][0188]
第n 1个样本点的核密度估计为：
[0189][0190]
式中，ρ(x
n 1
)为样本点在x
n 1
处的相对密度窗宽。
[0191]
通过以上递推可知，当任意t
n j
时刻增加j个样本时，
n j
个样本数据的核密度估计
可递推为：
[0192][0193]
3特征退化分布的计算
[0194]
假设tn为系统退化状态时间，可得到[0,tn]时间内所有退化样本的特征退化随时间变化的趋势可由图4所示。由于系统的退化增量是单调非负并且独立同分布的，每隔单位时间采集一次退化数据，计算其退化增量δx1,δx2,
…
,δxn为n个独立同分布的退化增量样本，的核密度估计表示为：
[0195][0196]
其中，ρ(δxi)(i＝1,2,3,
…
,n)为不同样本点处的相对密度窗宽。
[0197]
为使设备的剩余寿命预测准确，由核密度估计可得[0,tn]随机退化增量的概率密度推出[0,tn]上随机退化增量累计特征退化量xn(令xn＝nδx)的可由n重卷积得到：
[0198][0199]
当t
n 1
时刻增加一新样本时，[0,t
n 1
]的退化量的概率密度函数如下：
[0200][0201]
当t
n j
时刻新增j个样本时，[0,t
n j
]的退化量的概率密度函数为：
[0202][0203]
四、剩余寿命预测模型
[0204]
通过引入相对密度建立核估计的剩余寿命估计模型来保证核估计的平滑性和收
敛性，为了消除传统核密度估计的有界偏差问题，本发明在自适应相对密度的核密度估计模型的基础上引入核微分同胚变换的方法，通过空间映射的方式对剩余寿命预测模型进行变换，从而解决核函数剩余寿命预测模型在边界上估计的偏差性和无界性。
[0205]
4.1、核微分同胚估计
[0206]
设样本x1,x2,
…
,xn作为随机变量x恒等分布有概率密度函数f(x)，可以形成 fn的一个估计量序列，使只依赖于前n个观测值，设f为所有连续且有界的随机变量概率密度函数的集合，其支持范围限于实区间[a,b]，则核微分同胚估计量可表示为：
[0207][0208]
式中，k(
·
)是核函数；hn为带宽；是随机变量从区间[a,b]到的微分同胚变换；当x趋于下限或上限时，它的一阶导数趋于无穷。其中，和x的关系式如下：
[0209][0210]
4.2、核微分同胚估计的剩余寿命预测模型
[0211]
设tn为当前时刻，[0,tn]监测时间内采集到的当前样本退化数据，其相应的特征退化随时间退化的特征增量可以作为随机变量δx∈[0,x
th
](x
th
为失效阈值)。
[0212]
为消除核估计的边界偏移，使其准确性提高，对有界系统退化增量进行对数核微分同胚变换：
[0213][0214]
由上述将有界随机变量转换到实数域上，通过这种空间映射方式，建立微分同胚的核密度估计模型有效减少了边界处的自变量偏移问题。
[0215]
由公式(22)可推出，自适应相对密度的核微分同胚估计为：
[0216][0217]
k(
·
)是核函数，ρ(xi)是相对密度窗宽。
[0218]
设t
n t
时刻，特征退化量达到x
th
(如图4所示)时系统失效。要对当前tn时刻的剩余寿命进行预测，可通过初始时刻到当前tn时刻的退化量x
1:n
(记xn＝x(tn)， x
1:n
＝{x1,x2,...,xn})，推出t
n t
时刻的x
n t
。设t为设备的剩余寿命，则剩余寿命的概率密度分布函数f
t
(t)为：
[0219][0220]
其中，g(x
n t
)为[0,t
n t
]特征退化量的概率密度。
[0221]
通过对已知的单位时间随机退化特征增量核微分同胚估计，用卷积求其特征退化
量的概率密度函数，将对数微分同胚变换后[0,t
n t
]特征退化量的概率密度记为gd(x
n t
)：
[0222][0223]
剩余寿命预测的概率密度为：
[0224][0225]
式中表示t
n t
时刻的核微分同胚变换相对密度核估计的 n t次卷积，在实时监测数据的更新下，t
n t
时刻n t个样本核微分同胚变换的相对密度核估计为：
[0226][0227]
根据公式(25)和(27)可推出：
[0228][0229]
从而能够推出tn时刻系统剩余寿命预测的概率密度函数为：
[0230][0231]
随着实时监测的进行，监测到的样本数据不断增多，样本的核密度估计也随着不断更新，采用非实时的寿命预测模型时，每新增一个样本数据，基于这些样本的核密度估计都要重新计算，这样会造成历史样本不断重复计算，计算量也会越来越大，为避免实时监测系统中样本核密度估计不断重复计算的问题，提出对核密度估计模型实时更新的递推算法，进而实现对特征退化分布和实时剩余寿命的不断实时更新。
[0232]
因此，在实时变化的系统下获得新的样本数据后，通过可推算出下一时刻的进而实现实时剩余寿命预测。
[0233]
五、采用本预测方法的实验结果分析：
[0234]
5.1、齿轮磨损试验：
[0235]
以齿轮箱的齿轮为研究对象对本发明模型进行验证，如图5所示，其中心距为150mm，电机转速为1200r/min，本实验过程主要是针对加速度传感器的振动信号进行监测。
[0236]
如图6所示，1#～8#为加速度传感器(4#装设在轴承座的径向)；在主试箱和陪试箱的正上方40cm处安装了9#和10#为声音传感器；齿轮箱的温度是通过主试箱内的11#温度传感器来测。采用快速测定法，加载了八级载荷，在第八级载荷发生断齿。本实验主要对4#传感器记录齿轮箱在第八级载荷的加速度数据进行分析。采样频率为25.6khz，每次持续60s，
每9min记录一次数据。
[0237]
为了能够更好的展现退化趋势，可通过采用均方幅值方法对数据进行特征提取，减少后期预测中的误差。求得采样信号的均方幅值为：
[0238][0239]
式中，n为采样点数，yj为初始振动信号，yi为均方根幅值。
[0240]
如图7所示，该特征提取方法可以很好的把齿轮箱退化趋势展示出来。
[0241]
由图7知，在t∈[0,10]h时，齿轮处于啮齿阶段；t∈[10,68]h时，特征值逐渐增大，齿轮进入正常磨损；t∈[68,77]h时，齿轮磨损加剧，在77.17h发生断齿，此时齿轮的故障阈值为y＝76.325mm/s2。
[0242]
5.1.1、自适应相对密度窗宽的确定
[0243]
通过自适应相对密度方法对核密度估计的窗宽进行选取。如图8、图9和图 10所示，分别给出了t＝40h，70h，77h三个不同时刻样本点的相对密度窗宽选取。从图8、图9和图10中可以看出，随着样本数据的增多，相对密度窗宽方法能够自适应的选取窗宽，在不同样本密度下能够选择出合适的窗宽。
[0244]
5.1.2、齿轮的实时剩余寿命预测
[0245]
用本发明模型对齿轮实时剩余寿命预测，得到图11所示的自适应相对密度核估计的剩余寿命预测估计值与实际值在实时监测系统下不同监测时间的比较图。其中，剩余寿命预测值通过平均剩余寿命求得：
[0246][0247]
由图11分析可得，在初始阶段由于样本数据少，剩余寿命的预测值与实际值之间的误差较大；随着监测时间的变化，样本的增多，剩余寿命的概率密度越来越高，越来越窄，方差逐渐越小，预测值逐渐贴近实际值，表明预测的剩余寿命愈加准确。
[0248]
为了更加清楚的分析不同时刻剩余寿命的概率密度比较，表1给出了本发明所提方法的剩余寿命预测值和实际值的均方根误差。由表中数据可以分析出，在监测系统的实时更新，随着样本数据的不断增多，预测值与实际值之间的误差逐渐变小，验证了本发明方法的可行性。
[0249]
表1剩余寿命预测值与实际值的比较
[0250][0251]
5.1.3、基于固定、积分均方误差和相对密度三种窗宽准确性比较
[0252]
在实时监测系统下，随着数据的不断变化，要使选择的窗宽能够满足自适应地在密度小的区间选择较大的窗宽，在密度大的区间选择较小的窗宽。
[0253]
图12和图13分别给出了监测时间t＝40h和t＝70h时基于固定窗宽、自适应大拇指法则窗宽与自适应相对密度窗宽三种窗宽的核估计概率密度的比较。
[0254]
由图12和图13可知，基于自适应相对密度窗宽的方法与其他两种窗宽的方法相比，其结果更接近实际值。随着监测时间的变化，样本数据的不断增加，三个窗宽下的剩余寿命估计值与实际值之间的误差会逐渐变小，而本发明相对密度窗宽的方法相较于其他方法预测结果的误差相对更小。从整体分析来说，本发明所提的相对密度窗宽在核估计时可以更好的将数据拟合，并能更加准确的对概率密度函数进行估计，从而使剩余寿命预测的可靠性和准确度进一步提高。
[0255]
5.1.4、核估计的边界偏移问题
[0256]
在剩余寿命预测上，为了消除核估计的边界偏移问题，对有界随机变量进行对数核微分同胚变换将定义域转化到实数域上，再用核密度估计进行概率密度估计，从而提高预测的准确性。图14和图15给出在t＝40h和t＝70h两个不同时刻下本发明提出模型与传统核密度估计模型在核估计的偏差问题的比较。
[0257]
从图14和图15可以看出，两个不同时刻都反映出传统核密度估计存在有界性产生的边界偏差问题，而本发明提出的方法经过核微分同胚变换将定义域转化到实数域上，有效地消除了核估计的边界偏差。因此，本发明提出的模型能够有效解决边界偏移问题，提高了剩余寿命预测可靠性和准确性。
[0258]
5.1.5、不同模型比较
[0259]
目前在机器学习的剩余寿命预测方法中lstm神经网络方法对于故障时间序列有着较好的预测效果；gamma过程和wiener过程是两个最常用随机过程退化建模的方法，在寿命预测建模中应用较多。因此，为了验证本发明模型的精确度，采用上述三种方法与本发明模型进行比较分析。
[0260]
(1)、采用lstm神经网络方法对齿轮进行了剩余寿命预测。如图16所示，在周期数为455(10min)时，达到阈值76.325mm/s2，故其预测的剩余寿命为455-420＝35(10min)＝
5.83h。如图17所示，在周期数为456(10min)时，达到阈值76.325mm/s2，故其预测的剩余寿命为456-450＝6(10min)＝1h，与本发明模型具体的结果比较见表2。
[0261]
(2)、齿轮的退化是一个持续累积退化的过程，gamma过程由于具有非负、增长、独立增量的属性，被广泛用于逐渐累积损伤过程的退化建模中。为进一步测试本发明模型的准确性，在相同条件下，采用基于gamma过程的剩余寿命预测方法进行比较说明。图18和图19分别是监测时间t＝60h和t＝70h时两种模型的剩余寿命概率密度。
[0262]
通过图18和图19可分析出，本发明模型方法相对于基于gamma过程的方法来说，其剩余寿命的概率密度函数的方差愈来愈小，可以更好地反映样本数据，预测值更加贴近真实值。t＝70h时能够明显看出本发明方法的预测值更接近实际值。因此，本发明模型在剩余寿命预测的应用上更加准确有效。
[0263]
(3)、wiener过程模型多用于对具有非单调退化过程的设备进行建模。在相同条件下采用wiener过程模型对齿轮的剩余寿命进行预测。图20和图21分别给出了监测时间t＝60h和t＝70h时本发明模型和wiener过程模型的剩余寿命概率密度比较。
[0264]
由图20和图21可以看出，本发明模型与wiener过程模型相比，剩余寿命的概率密度函数的方差愈来愈小，预测值更加贴近真实值，预测的准确性提高。
[0265]
此外，通过引入相对误差指标进一步量化预测的精确度，给出监测时间在 t＝60h，65h，70h，75h，77h五个时刻下四种方法的比较结果。
[0266]
通过表2可知，在不同时刻的相对误差比较下，四种模型的剩余寿命预测结果都逐渐趋于真实值，而本发明模型的预测结果与lstm模型、gamma过程以及基于wiener过程三种方法相比，其相对误差更小，从而验证了本发明模型的优越性。
[0267]
表2不同监测时刻下相对误差的比较结果
[0268][0269]
5.2、滚动轴承的加速寿命试验
[0270]
为验证本发明方法的有效性，利用ieee phm2012提供的轴承全寿命数据对模型进行验证。本发明以转速为1800r/min，载荷为4000n工况下的bearing 1-1 的全寿命振动数据为例进行分析。该数据来源于femto-st研究中心 pronostia试验台对滚动轴承的加速寿命试验，振动信号的采样频率为25.6 khz，10s采集一次，一次采集0.1s，即一次采集2560个样本点。
[0271]
通过均方幅值法处理bearing 1-1的全寿命振动数据，如图21所示的特征值随监测时间变化曲线。从图22可以看出，均方根随时间基本呈现单调增加的趋势，能较好地反映
其退化趋势，该轴承在t＝2.749
×
104s时，磨损开始加剧，且在t＝2.803
×
104s时失效，均方根的失效阈值为5.607mm/s2。
[0272]
采用本发明模型对轴承进行剩余寿命预测，由图23可知，随着系统运行时间的增加，接收到的监测样本不断增多，基于相对密度核估计的剩余寿命的概率密度不断实时更新，剩余寿命的概率密度变窄变高，方差越来越小，说明预测的准确性不断提高。
[0273]
为进一步对本发明提出方法的预测效果进行评估，表3给出了对不同监测时间，实际剩余寿命、本发明模型预测的平均剩余寿命的均方根误差比较。
[0274]
表3剩余寿命预测值与实际值的比较
[0275][0276]
从表3可以看出，随监测时间的增加，rmse呈现逐渐减小的趋势，且随着监测数据的增多，预测的剩余寿命与真实寿命的误差更小，说明本发明模型预测的剩余寿命更接近实际的寿命值。
[0277]
该模型将k近邻思想计算出的样本点的相对密度作为核密度估计的窗宽，不仅解决了固定窗宽由样本数据分布不均导致拟合不足的问题，而且能够对任意形状，密度不均匀的数据集自适应地选择出更加准确的窗宽，提高了拟合度。在核估计的边界偏差上，引入核微分同胚变换的方法有效的消除了核密度估计带来的边界偏差问题。此外，随着样本的增加，为了避免核密度估计的重复计算，建立了基于相对密度核估计的实时更新模型。实例分析表明，随着样本数据的增加，剩余寿命的预测越来越接近实际值，方差变得越来越小，提高了预测的准确性；并且通过与lstm模型、基于wiener过程以及基于gamma过程的预测方法比较，进一步验证了本发明模型的准确性和有效性。
[0278]
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包在本发明范围内。