齿轮箱中考虑轴承退化影响的齿轮剩余寿命预测方法

1.本发明属于机械可靠性技术领域，具体涉及一种齿轮箱中考虑轴承退化影响的齿轮剩余寿命预测方法，


背景技术：

2.在工业生产中齿轮减速箱被广泛的应用在旋转机械设备传动系统中，当齿轮发生故障时，会增加齿轮传动的振动，导致转速不稳定，从而增加对轴承的损害，轴承是齿轮箱的重要组成部分，通常需要承受转动方向改变和受力变化，一旦发生失效将会导致整个传动系统发生崩溃，如果轴承的间隙过大，径向振动也会随之变大，从而对齿轮造成冲击；如果间隙太小，齿轮间的摩擦阻力就会变大，严重影响了齿轮的寿命，现代化设备中系统越来越复杂，系统中有许多重要的子系统和零部件协同工作，部件之间退化过程中存在随机相关性影响，


技术实现要素：

3.本发明目的是提供一种齿轮箱中考虑轴承退化影响的齿轮剩余寿命预测方法，该方法可有效地预测齿轮箱的寿命，提供预防措施，
4.本发明是这样实现的，其特征在于实施步骤如下：
5.步骤1、通过试验获取表征斜齿轮箱内的实时监测数据，采用如图1、图2 所示的斜齿轮箱试验台进行试验，其上安装的主试齿轮箱gb1和陪试齿轮箱gb2 采用背靠背的结构进行安装，分别与驱动电机1和负载电机2相连，在主试齿轮箱gb1和陪试齿轮箱gb2上均安装有振动传感器、声学传感器和油温传感器， 16通道数据采集系统das分别与主试齿轮箱gb1和陪试齿轮箱gb2及计算机pc 相连接，驱动电机1同时与电气测量装置4和16通道数据采集系统das相连接，控制柜3同时与电气测量装置4及负载电机2相连接，16通道数据采集系统das，将采集到的振动信号、声学信号和油温信号转换为数字信号，然后传输到计算机pc进行后期msb频谱分析，计算机pc将16通道数据采集系统das采集到的数据进行整理；控制柜3安装在距试验台约3米处，其内安装有电流和电压传感器，进行对负载电机2远距离信号检测，电流传感器其灵敏度是5a/v，频率范围是0～1000hz，对该供电电源在变频、变负荷的瞬时电流进行有效测量，电气测量装置4与控制柜3和驱动电机1相连接，得到控制柜3的指令后将电流信号输入16通道数据采集系统das中，将电压信号输入驱动电机1让驱动电机运转，驱动电机1上安装有振动传感器用于采集振动信号，采集到的信号通过16通道数据采集系统das传输到计算机pc；驱动电机1额定转速为1500rpm，振动信号的采样频率为96khz；
6.在试验期间，首先控制驱动电机1以正弦曲线变化的转速运转30min，然后分别在0％负载、25％负载、50％负载、75％负载、100％负载5种不同负载条件下以全速的50％运转，最后同样的分别在5种负载下以全速的70％运行；
7.步骤2、对齿轮的振动数据进行msb调制信号双谱变换后对齿轮的退化数据进行特征提取并预测其剩余寿命，常规的双谱b(f
x
,fc)通过傅里叶变换后在频域内可表示为：
8.b(f
x
,fc)＝e[x(fc)x(f
x
)x
*
(fc+f
x
)]
ꢀꢀꢀ
(1)
[0009]
其中x(f)为振动信号x(t)的离散傅立叶变换；e(g)为求均值符号；f
x
为调制频率；fc为载波频率；x
*
为x的复共轭，
[0010]
对于振动信号x(t)其msb变换在频域中表示为：
[0011]bms
(f
x
,fc)＝e[x(fc+f
x
)x(f
c-f
x
)x
*
(fc)x
*
(fc)]
ꢀꢀꢀ
(2)
[0012]
msb对常规双谱改进之后在进行退化特征提取时，同时考虑信号的幅值和相位的变化，能够更好地说明调制信号和载波信号之间的非线性关系，充分表示信号的调制特性，要比常规的双谱分析方法更准确，
[0013]
以幅值和相位来表示，可将式(2)写为：
[0014][0015]
公式(3)中msb的总相位可通过以下方式计算：
[0016][0017]
和耦合时，相位关系可以表示为：
[0018][0019][0020]
步骤3、从多部件系统部件间存在的复杂随机相关性的特征以及对部件连续退化状态的不同影响角度出发对部件间的随机相关性进行分析，可将其分为三类：单一单向随机相关性，单一多向随机相关性和双向随机相关性，
[0021]
1)单一单向随机相关性：是指某个部件的退化只会对系统中的单个部件的退化产生单向影响；
[0022]
2)单一多向随机相关性：是指某个部件的退化会对系统中多个部件的退化产生单向影响或者是某个部件的退化会受到系统中其它多个部件退化的单向影响；
[0023]
3)双向随机相关性：是指部件自身的退化会受到系统中其它部件的影响，同时自身也会对其它部件的退化产生影响，本发明考虑的是更为复杂的部件间具有双向随机相关性的多部件系统，假设一个系统中的部件i和部件j具有双向随机相关性，即一个部件的退化会加剧另一个部件的退化，反之亦然；
[0024]
在考虑更为复杂的部件间具有双向随机相关性的多部件系统进行研究时，对双向随机相关性的模型进行构建，假设部件i在t
k-1
时刻的退化量为则它在tk时刻的退化量为：
[0025][0026]
其中为部件i在[0,t
k-1
]内的k-1个独立同分布的退化增量，假设部件j在时间t
k-1
的退化量为则它在tk时刻的退化量为：
[0027][0028]
[0029]
其中，为部件j在tk时刻的退化量；为部件j自身固有的退化增量；为部件i对部件j随机相关影响的退化增量，
[0030]
将其概率密度函数记为则的核密度估计为：
[0031][0032]
式中，∑为求和符号；k(g)为核函数；为窗宽；k-1为部件i在t
k-1
时刻退化增量的样本数；为包含部件i自身固有退化增量及部件j对部件i随机相关影响退化增量的样本值，
[0033]
采用积分均方误差对核估计进行衡量：
[0034][0035]
其中mise代表求积分均方误差；∫为积分符号，
[0036]
不同的核函数对结果的影响很小，选择常用的gaussian核来建模：
[0037][0038]
通过求解mise的最小值可得到窗宽的最优解：
[0039][0040]
其中，σ
k-1
为k-1个已知样本的标准差；n是样本的数量，
[0041]
将公式(11)代入上式，可以求出：
[0042][0043]
同理，部件i在运行过程中自身固有的退化增量为其核密度估计可表示为：
[0044][0045]
部件j在此系统运行期间的退化增量为其核密度估计可表示为：
[0046][0047]
部件j对部件i随机相关影响的退化增量记为将其服从的概率密度函数记为f
k-1
(δx
ji
)，由式(9)知，为考虑部件i自身固有退化及部件j对部件i 随机相关影响的概率密度函数，则f
k-1
(δx
ji
)的核密度估计为：
[0048][0049]
步骤4、引入窗宽因子λn采用自适应窗宽核密度估计更能反映样本数据的真实情况，当时，核密度估计与真实值之间的误差很小，则tk时刻获得新样本数据后的自适应窗宽可表示为：
[0050][0051]
其中，为tk时刻新的自适应窗宽；λk为窗宽因子；为部件i在窗宽下tk时刻的核密度估计值，
[0052]
部件i在tk时刻基于自适应窗宽的核密度估计为：
[0053][0054]
其中，“*”为卷积符号；为部件i在t
k-1
时刻的核密度估计值；为部件i在tk时刻自身退化影响的核密度估计值；为部件j对部件i随机相关影响的核密度估计值，
[0055]
对退化分布进行计算，得到部件i基于自适应窗宽的核密度估计表达式后，可以通过t
k-1
时刻的核密度估计来求部件i在t
k-1
时刻累积退化的概率密度函数，记为
[0056][0057]
部件i在t
k-1
时刻累积退化的概率密度函数可通过t
k-1
时刻核密度估计值的k-1重卷积得到，
[0058]
tk时刻获得新样本数据后，概率密度函数为：
[0059][0060]
t
k-1+m
时刻获得m个样本数据后，概率密度函数为：
[0061][0062]
同理，在任意时刻获得新的样本数据后，都可获得相应的累积退化概率密度函数；
[0063]
步骤5、建立剩余寿命预测模型，给定部件i的失效阈值为设部件i在t
k-1+t
时刻的退化量为时，系统发生故障，那么部件i在t
k-1
时刻的剩余寿命概率分布函数f
t
(t)为：
[0064][0065]
由公式(21)可知：
[0066][0067]
将公式(23)代入公式(22)，可得：
[0068][0069]
则t
k-1
时刻部件i的剩余寿命概率密度函数为：
[0070][0071]
本发明优点及积极效果是：
[0072]
可自适应地选择出更加准确的窗宽，提高了拟合度，可有效地预测考虑轴承退化影响的齿轮剩余寿命，为后续的维修决策提供依据，
附图说明
[0073]
图1为实施本发明采用的测试台示意图；
[0074]
图2为图1中主试齿轮箱gb1与陪试齿轮箱gb2的结构示意图；
[0075]
图3为本发明实施例中模拟设备在实际运行过程中环境与负载的随机变化状况示意图；
[0076]
图4为本发明实施例中msb调制信号双谱分析图；
[0077]
图5为本发明实施例中部件i的退化状态图；
[0078]
图6为本发明实施例中随机相关性示意图；
[0079]
图7为本发明实施例中核密度估计在不同监测时间下剩余寿命预测结果图；
[0080]
图中：1-驱动电机 2-负载电机 3-控制柜 4-电气测量装置
[0081]
5.1-第一联轴器 5.2-第二联轴器 5.3-第三联轴器
[0082]
das-16通道数据采集系统 pc-计算机 gb1-主试齿轮箱 gb2-陪试齿轮箱
[0083]z1-主试齿轮箱中第一对齿轮中主动齿轮的齿数
[0084]z2-主试齿轮箱中第一对齿轮中从动齿轮的齿数
[0085]z3-主试齿轮箱中第二对齿轮中传动齿轮的齿数
[0086]z4-主试齿轮箱中第二对齿轮中从动齿轮的齿数
[0087]z5-陪试齿轮箱中第一对齿轮中主动齿轮的齿数
[0088]z6-陪试齿轮箱中第一对齿轮中从动齿轮的齿数
[0089]z7-陪试齿轮箱中第二对齿轮中传动齿轮的齿数
[0090]z8-陪试齿轮箱中第二对齿轮中从动齿轮的齿数
具体实施方式
[0091]
下面结合附图对本发明实施例做进一步说明：
[0092]
如图1所示，实施步骤如下：
[0093]
步骤1、通过试验获取表征斜齿轮箱内的实时监测数据，采用如图1、图2 所示的斜齿轮箱试验台进行试验，其中主试齿轮箱gb1和陪试齿轮箱gb2采用背靠背的结构进行安装，分别与驱动电机1和负载电机2相连接，16通道数据采集系统das分别与主试齿轮箱gb1和陪试齿轮箱gb2及计算机pc相连接，驱动电机1同时与电气测量装置4和16通道数据采集系统das相连接，控制柜3 同时与电气测量装置4及负载电机2相连接，16通道数据采集系统das，将从传感器获取的模拟信号转换为数字信号，传输到计算机pc进行后期msb频谱分析及数据整理；控制柜3安装在距试验台约3米处，其内安装有电流和电压传感器，进行对负载
电机2远距离信号检测，电流传感器其灵敏度是5a/v，频率范围是0～1000hz，对该供电电源在变频、变负荷的瞬时电流进行有效测量，电气测量装置4与控制柜3和驱动电机1相连接，得到控制柜3的指令后将电流信号输入16通道数据采集系统das中，将电压信号输入驱动电机1让驱动电机运转，驱动电机1上安装有振动传感器用于采集振动信号，采集到的信号通过16通道数据采集系统das传输到计算机pc；驱动电机1额定转速为1500rpm，振动信号的采样频率为96khz；上述试验台上安装的主试齿轮箱gb1与陪试齿轮箱gb2的设备参数值为表1所示，
[0094]
表1齿轮箱规格
[0095][0096]
在试验期间，首先控制驱动电机1以正弦曲线变化的转速运转30min，然后分别在0％负载、25％负载、50％负载、75％负载、100％负载5种不同负载条件下以全速的50％运转，最后同样的分别在5种负载下以全速的70％运行，表2是轴频和啮合频率的公式计算；
[0097]
表2轴频和啮合频率
[0098][0099]
步骤2、对齿轮的振动数据进行msb调制信号双谱变换后对齿轮的退化数据进行特征提取并预测其剩余寿命，本实施例试验运行了838h，监测到陪试齿轮箱gb2的振动出现明显峰值后停止试验，打开齿轮箱后，观察到陪试齿轮箱gb2 在低速阶段发生磨损，齿根部位磨损较为严重，齿轮在运行初期样本特征波动较大，为齿轮啮合阶段，使用这个阶段的数据来分析部件的磨损退化过程是不恰当的，因此选择运行300h之后的数据来分析，对振动数据进行调制信号双谱变换，进行特征提取并预测其剩余寿命，常规的双谱b(f
x
,fc)通过傅里叶变换后在频域内可表示为：
[0100]
b(f
x
,fc)＝e[x(fc)x(f
x
)x
*
(fc+f
x
)]
ꢀꢀꢀ
(1)
[0101]
其中x(f)为信号x(t)的离散傅立叶变换；e(g)为求均值符号；f
x
为调制频率； fc为载波频率；x
*
为x的复共轭，
[0102]
对于振动信号x(t)其msb变换在频域中表示为：
[0103]bms
(f
x
,fc)＝e[x(fc+f
x
)x(f
c-f
x
)x
*
(fc)x
*
(fc)]
ꢀꢀꢀ
(2)
[0104]
msb对常规双谱改进之后在进行退化特征提取时，同时考虑信号的幅值和相位的变化，能够更好地说明调制信号和载波信号之间的非线性关系，充分表示信号的调制特性，要比常规的双谱分析方法更准确，
[0105]
以幅值和相位来表示，可将式(2)写为：
[0106][0107]
公式(3)中msb的总相位可通过以下方式计算：
[0108][0109]
和耦合时，相位关系可以表示为
[0110][0111][0112]
将公式(5)和公式(6)代入公式(4)，可以得到msb的总相位为零，它的幅值由其四个分量幅值的乘积所确定，因此，如果(fc+f
x
)和(f
c-f
x
)来源于f
x
和 fc的非线性影响，在双频b
ms
(f
x
,fc)处会出现明显的双谱峰，这样对msb的表述将更为准确，与此相反，如果f
x
和fc无非线性作用，或者它们是随机噪声，其msb 的相位随时间变化，最终平均结果为零，由此msb具有很好地抑制噪声作用，图4是采用msb对试验台数据特征提取的调制双谱分析图，通过msb双谱分析先从幅值谱发现具有峰值的信号幅值，同时，该坐标点对应的相位谱中也有较大的相干系数，则说明该幅值可以反映部件的退化特征，若相位谱中对应着非常小的相干系数，说明该幅值不是通过调制作用生成，可能是由噪声产生的，从而将其剔除，其在识别非线性影响的同时考虑幅值和相位的影响，从而有效地抑制随机噪声的干扰，提取能够表征部件退化的振动信号特征数据出来，便可得到部件的退化状态曲线图，图2中gb2中的齿轮z5和gb1中的齿轮z3通过主轴进行连接，二者之间具有较强的相关性，故将齿轮z5和齿轮z3作为部件i和部件j进行研究，图5是部件i的退化状态曲线图；
[0113]
步骤3、从多部件系统部件间存在的复杂随机相关性的特征以及对部件连续退化状态的不同影响角度出发对部件间的随机相关性进行分析，可将其分为三类：单一单向随机相关性，单一多向随机相关性和双向随机相关性，图6是本实施中部件间随机相关性示意图，
[0114]
1)单一单向随机相关性：是指某个部件的退化只会对系统中的单个部件的退化产生单向影响；
[0115]
2)单一多向随机相关性：是指某个部件的退化会对系统中多个部件的退化产生单向影响或者是某个部件的退化会受到系统中其它多个部件退化的单向影响；
[0116]
3)双向随机相关性：是指部件自身的退化会受到系统中其它部件的影响，同时自身也会对其它部件的退化产生影响，本实施例考虑的是更为复杂的部件间具有双向随机相关性的多部件系统，假设一个系统中的部件i和部件j具有双向随机相关性，即一个部件的退化会加剧另一个部件的退化，反之亦然，本发明考虑更为复杂的部件间具有双向随机相关性的多部件系统进行研究，假设一个系统中的部件i和部件j具有双向随机相关性，即一个部件的退化会加剧另一个部件的退化，反之亦然，随着现代传感器技术的发展，可以通过
传感器监测系统中各部件的历史退化状态数据及实时运行数据，利用这些样本数据进行随机相关性建模并建立剩余寿命预测模型，
[0117]
在考虑更为复杂的部件间具有双向随机相关性的多部件系统进行研究时，对双向随机相关性的模型进行构建，假设部件i在t
k-1
时刻的退化量为则它在tk时刻的退化量为：
[0118][0119]
为考虑部件之间的随机相关性影响，其部件i在单位时间t
k-t
k-1
的退化增量为：
[0120][0121]
其中，为部件i自身固有的退化增量；为部件j对部件i随机相关影响的退化增量，
[0122]
同理，假设部件j在时间t
k-1
的退化量为则它在tk时刻的退化量为：
[0123][0124][0125]
其中，为部件j在tk时刻的退化量；为部件j自身固有的退化增量；为部件i对部件j随机相关影响的退化增量，
[0126]
以部件i为例，设为部件i在[0,t
k-1
]内的k-1个独立同分布的退化增量，将其概率密度函数记为则的核密度估计为：
[0127][0128]
其中，∑为求和符号；k(g)为核函数；为窗宽；k-1为部件i在t
k-1
时刻退化增量的样本数；为包含部件i自身固有退化增量及部件j对部件i随机相关影响退化增量的样本值，采用积分均方误差对核估计进行衡量：
[0129][0130]
其中mise代表求积分均方误差；∫为积分符号，
[0131]
核函数k(g)和窗宽是核密度估计的两个重要参数，核函数的种类很多，但是不同的核函数对结果的影响很小，选择常用的gaussian核来建模：
[0132][0133]
窗宽的选择对核密度估计结果的准确性有很大影响，通过求解mise的最小值可得到窗宽的最优解：
[0134]
[0135]
其中，σ
k-1
为k-1个已知样本的标准差；n是样本的数量，
[0136]
将公式(11)代入上式，可以求出：
[0137][0138]
同理，部件i在运行过程中自身固有的退化增量为其核密度估计可表示为：
[0139][0140]
部件j在此系统运行期间的退化增量为其核密度估计可表示为：
[0141][0142]
部件j对部件i随机相关影响的退化增量记为将其服从的概率密度函数记为f
k-1
(δx
ji
)，由式(9)知，为考虑部件i自身固有退化及部件j对部件i随机相关影响的概率密度函数，则f
k-1
(δx
ji
)的核密度估计为：
[0143][0144]
步骤4、引入窗宽因子λn采用自适应窗宽核密度估计更能反映样本数据的真实情况，当时，核密度估计与真实值之间的误差很小，则tk时刻获得新样本数据后的自适应窗宽可表示为：
[0145][0146]
其中，为tk时刻新的自适应窗宽；λk为窗宽因子；为部件i在窗宽下tk时刻的核密度估计值，
[0147]
部件i在tk时刻基于自适应窗宽的核密度估计为：
[0148]
[0149]
其中，“*”为卷积符号；为部件i在t
k-1
时刻的核密度估计值；为部件i在tk时刻自身退化影响的核密度估计值；为部件j对部件i随机相关影响的核密度估计值，采用递推来实现核密度估计的实时更新，下一时刻有新增样本时，通过退化增量便可更新计算，大大减少了计算量，降低了计算复杂度，
[0150]
对退化分布进行计算，得到部件i基于自适应窗宽的核密度估计表达式后，可以通过t
k-1
时刻的核密度估计来求部件i在t
k-1
时刻累积退化的概率密度函数，记为
[0151][0152]
部件i在t
k-1
时刻累积退化的概率密度函数可通过t
k-1
时刻核密度估计值的k-1重卷积得到，
[0153]
tk时刻获得新样本数据后，概率密度函数为：
[0154][0155]
t
k-1+m
时刻获得m个样本数据后，概率密度函数为：
[0156][0157]
同理，在任意时刻获得新的样本数据后，都可获得相应的累积退化概率密度函数；
[0158]
步骤5、建立剩余寿命预测模型，给定部件i的失效阈值为设部件i在t
k-1+t
时刻的退化量为时，系统发生故障，那么部件i在t
k-1
时刻的剩余寿命概率分布函数f
t
(t)
为：
[0159][0160]
由公式(21)可知：
[0161][0162]
将公式(23)代入公式(22)，可得：
[0163][0164]
则t
k-1
时刻部件i的剩余寿命概率密度函数为：
[0165][0166]
将msb特征提取后的退化数据代入模型，可根据式22-式25求得部件i在当前监测时刻的剩余寿命概率密度函数，图7是本次试验在7个监测时间下的剩余寿命预测结果，可以看出，随着监测时间的增加，样本数据增大，本发明提出方法的剩余寿命预测结果越来越接近真实值，说明本发明提出的基于自适应核窗宽的核密度估计剩余寿命预测方法可以很好的对部件i剩余寿命的概率密度函数进行估计，
[0167]
表3给出了参数估计和核密度估计方法对部件剩余寿命预测结果的误差分析，从表中可以看出，随着监测时间的推移，获得的样本数据越多，两种方法对部件剩余寿命的预测就越精确，在同一监测时间下，本发明提出的核密度估计的方法要比参数估计方法更精确，进一步验证了非参数核密度估计方法的有效性和精确性，
[0168]
表3参数估计和核估计剩余寿命预测误差分析
[0169][0170]
综上所述，本发明针对多部件系统部件连续退化过程中存在的双向随机相关性，基于自适应窗宽的核密度估计方法建立了实时剩余寿命预测模型，首先针对部件间存在的双向随机相关性影响，引入条件核密度估计来进行建模，然后采用自适应窗宽的核密度估计方法求得相应的概率密度函数，最后建立实时剩余寿命预测模型，得到部件的剩余寿命。