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

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

Info

Publication number
CN114942139A
CN114942139A CN202210292839.1A CN202210292839A CN114942139A CN 114942139 A CN114942139 A CN 114942139A CN 202210292839 A CN202210292839 A CN 202210292839A CN 114942139 A CN114942139 A CN 114942139A
Authority
CN
China
Prior art keywords
component
degradation
time
gear
gear box
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202210292839.1A
Other languages
English (en)
Inventor
石慧
李哲昊
康辉
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Taiyuan University of Science and Technology
Original Assignee
Taiyuan University of Science and Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Taiyuan University of Science and Technology filed Critical Taiyuan University of Science and Technology
Priority to CN202210292839.1A priority Critical patent/CN114942139A/zh
Publication of CN114942139A publication Critical patent/CN114942139A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M13/00—Testing of machine parts
    • G01M13/04—Bearings
    • G01M13/045—Acoustic or vibration analysis
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01D—MEASURING NOT SPECIALLY ADAPTED FOR A SPECIFIC VARIABLE; ARRANGEMENTS FOR MEASURING TWO OR MORE VARIABLES NOT COVERED IN A SINGLE OTHER SUBCLASS; TARIFF METERING APPARATUS; MEASURING OR TESTING NOT OTHERWISE PROVIDED FOR
    • G01D21/00—Measuring or testing not otherwise provided for
    • G01D21/02—Measuring two or more variables by means not covered by a single other subclass
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N20/00—Machine learning
    • G06N20/10—Machine learning using kernel methods, e.g. support vector machines [SVM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N7/00—Computing arrangements based on specific mathematical models
    • G06N7/01—Probabilistic graphical models, e.g. probabilistic networks
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/02—Preprocessing
    • G06F2218/04—Denoising
    • G06F2218/06—Denoising by applying a scale-space analysis, e.g. using wavelet analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/08—Feature extraction
    • G06F2218/10—Feature extraction by analysing the shape of a waveform, e.g. extracting parameters relating to peaks
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
    • Y02E10/00—Energy generation through renewable energy sources
    • Y02E10/70—Wind energy
    • Y02E10/72—Wind turbines with rotation axis in wind direction

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Software Systems (AREA)
  • Computing Systems (AREA)
  • Artificial Intelligence (AREA)
  • Data Mining & Analysis (AREA)
  • Evolutionary Computation (AREA)
  • Mathematical Physics (AREA)
  • General Engineering & Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Algebra (AREA)
  • Probability & Statistics with Applications (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Medical Informatics (AREA)
  • Acoustics & Sound (AREA)
  • Testing Of Devices, Machine Parts, Or Other Structures Thereof (AREA)

Abstract

一种考虑轴承退化影响的齿轮剩余寿命的预测方法，属于机械可靠性技术领域，其特征是实施步骤如下：1、利用传感器对齿轮箱内齿轮和轴承退化实时监测；2、对齿轮箱内齿轮和轴承的退化状态进行特征提取，利用调制信号双谱变换对齿轮和轴承磨损退化性能进行特征提取及退化评估；3、对系统中齿轮和轴承存在连续退化双向随机相关影响进行相关性建模；4、完全从数据的角度构建基于自适应窗宽的条件核密度估计退化模型得到部件在当前时刻的考虑轴承退化影响的齿轮剩余寿命预测模型；5、最后得齿轮剩余寿命概率密度函数；优点是可自适应地选择出窗宽，有效地预测考虑轴承退化影响的齿轮剩余寿命，为后续的维修决策提供依据。

Description

齿轮箱中考虑轴承退化影响的齿轮剩余寿命预测方法
技术领域
本发明属于机械可靠性技术领域，具体涉及一种齿轮箱中考虑轴承退化影响的齿轮剩余寿命预测方法，
背景技术
在工业生产中齿轮减速箱被广泛的应用在旋转机械设备传动系统中，当齿轮发生故障时，会增加齿轮传动的振动，导致转速不稳定，从而增加对轴承的损害，轴承是齿轮箱的重要组成部分，通常需要承受转动方向改变和受力变化，一旦发生失效将会导致整个传动系统发生崩溃，如果轴承的间隙过大，径向振动也会随之变大，从而对齿轮造成冲击；如果间隙太小，齿轮间的摩擦阻力就会变大，严重影响了齿轮的寿命，现代化设备中系统越来越复杂，系统中有许多重要的子系统和零部件协同工作，部件之间退化过程中存在随机相关性影响，
发明内容
本发明目的是提供一种齿轮箱中考虑轴承退化影响的齿轮剩余寿命预测方法，该方法可有效地预测齿轮箱的寿命，提供预防措施，
本发明是这样实现的，其特征在于实施步骤如下：
步骤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；
在试验期间，首先控制驱动电机1以正弦曲线变化的转速运转30min，然后分别在0％负载、25％负载、50％负载、75％负载、100％负载5种不同负载条件下以全速的50％运转，最后同样的分别在5种负载下以全速的70％运行；
步骤2、对齿轮的振动数据进行MSB调制信号双谱变换后对齿轮的退化数据进行特征提取并预测其剩余寿命，常规的双谱B(fx,fc)通过傅里叶变换后在频域内可表示为：
B(fx,fc)＝E[X(fc)X(fx)X*(fc+fx)] (1)
其中X(f)为振动信号x(t)的离散傅立叶变换；E(g)为求均值符号；fx为调制频率；fc为载波频率；X*为X的复共轭，
对于振动信号x(t)其MSB变换在频域中表示为：
BMS(fx,fc)＝E[X(fc+fx)X(fc-fx)X*(fc)X*(fc)] (2)
MSB对常规双谱改进之后在进行退化特征提取时，同时考虑信号的幅值和相位的变化，能够更好地说明调制信号和载波信号之间的非线性关系，充分表示信号的调制特性，要比常规的双谱分析方法更准确，
以幅值和相位来表示，可将式(2)写为：
Figure RE-GDA0003739368750000021
公式(3)中MSB的总相位可通过以下方式计算：
Figure RE-GDA0003739368750000031
Figure RE-GDA0003739368750000032
和
Figure RE-GDA0003739368750000033
耦合时，相位关系可以表示为：
Figure RE-GDA0003739368750000034
Figure RE-GDA0003739368750000035
步骤3、从多部件系统部件间存在的复杂随机相关性的特征以及对部件连续退化状态的不同影响角度出发对部件间的随机相关性进行分析，可将其分为三类：单一单向随机相关性，单一多向随机相关性和双向随机相关性，
1)单一单向随机相关性：是指某个部件的退化只会对系统中的单个部件的退化产生单向影响；
2)单一多向随机相关性：是指某个部件的退化会对系统中多个部件的退化产生单向影响或者是某个部件的退化会受到系统中其它多个部件退化的单向影响；
3)双向随机相关性：是指部件自身的退化会受到系统中其它部件的影响，同时自身也会对其它部件的退化产生影响，本发明考虑的是更为复杂的部件间具有双向随机相关性的多部件系统，假设一个系统中的部件i和部件j具有双向随机相关性，即一个部件的退化会加剧另一个部件的退化，反之亦然；
在考虑更为复杂的部件间具有双向随机相关性的多部件系统进行研究时，对双向随机相关性的模型进行构建，假设部件i在tk-1时刻的退化量为
Figure RE-GDA0003739368750000036
则它在tk时刻的退化量为：
Figure RE-GDA0003739368750000037
其中
Figure RE-GDA0003739368750000038
为部件i在[0,tk-1]内的k-1个独立同分布的退化增量，假设部件j在时间tk-1的退化量为
Figure RE-GDA00037393687500000314
则它在tk时刻的退化量为：
Figure RE-GDA0003739368750000039
Figure RE-GDA00037393687500000310
其中，
Figure RE-GDA00037393687500000311
为部件j在tk时刻的退化量；
Figure RE-GDA00037393687500000312
为部件j自身固有的退化增量；
Figure RE-GDA00037393687500000313
为部件i对部件j随机相关影响的退化增量，
将其概率密度函数记为
Figure RE-GDA0003739368750000041
则
Figure RE-GDA0003739368750000042
的核密度估计为：
Figure RE-GDA0003739368750000043
式中，∑为求和符号；K(g)为核函数；
Figure RE-GDA0003739368750000044
为窗宽；k-1为部件i在tk-1时刻退化增量的样本数；
Figure RE-GDA0003739368750000045
为包含部件i自身固有退化增量及部件j对部件i随机相关影响退化增量的样本值，
采用积分均方误差对核估计进行衡量：
Figure RE-GDA0003739368750000046
其中MISE代表求积分均方误差；∫为积分符号，
不同的核函数对结果的影响很小，选择常用的Gaussian核来建模：
Figure RE-GDA0003739368750000047
通过求解MISE的最小值可得到窗宽
Figure RE-GDA0003739368750000048
的最优解：
Figure RE-GDA0003739368750000049
其中，σk-1为k-1个已知样本的标准差；n是样本的数量，
将公式(11)代入上式，可以求出：
Figure RE-GDA00037393687500000410
同理，部件i在运行过程中自身固有的退化增量为
Figure RE-GDA00037393687500000411
其核密度估计
Figure RE-GDA00037393687500000412
可表示为：
Figure RE-GDA00037393687500000413
部件j在此系统运行期间的退化增量为
Figure RE-GDA00037393687500000414
其核密度估计
Figure RE-GDA00037393687500000415
可表示为：
Figure RE-GDA00037393687500000416
部件j对部件i随机相关影响的退化增量记为
Figure RE-GDA0003739368750000051
将其服从的概率密度函数记为fk-1(Δxji)，由式(9)知，
Figure RE-GDA0003739368750000052
为考虑部件i自身固有退化及部件j对部件i 随机相关影响的概率密度函数，则fk-1(Δxji)的核密度估计为：
Figure RE-GDA0003739368750000053
步骤4、引入窗宽因子λn采用自适应窗宽核密度估计更能反映样本数据的真实情况，当
Figure RE-GDA0003739368750000054
时，核密度估计与真实值之间的误差很小，则tk时刻获得新样本数据后的自适应窗宽
Figure RE-GDA00037393687500000515
可表示为：
Figure RE-GDA0003739368750000055
其中，
Figure RE-GDA0003739368750000056
为tk时刻新的自适应窗宽；λk为窗宽因子；
Figure RE-GDA0003739368750000057
为部件i在窗宽
Figure RE-GDA0003739368750000058
下tk时刻的核密度估计值，
部件i在tk时刻基于自适应窗宽的核密度估计为：
Figure RE-GDA0003739368750000059
其中，“*”为卷积符号；
Figure RE-GDA00037393687500000510
为部件i在tk-1时刻的核密度估计值；
Figure RE-GDA00037393687500000511
为部件i在tk时刻自身退化影响的核密度估计值；
Figure RE-GDA00037393687500000512
为部件j对部件i随机相关影响的核密度估计值，
对退化分布进行计算，得到部件i基于自适应窗宽的核密度估计表达式后，可以通过tk-1时刻的核密度估计
Figure RE-GDA00037393687500000513
来求部件i在tk-1时刻累积退化的概率密度函数，记为
Figure RE-GDA00037393687500000514
Figure RE-GDA0003739368750000061
部件i在tk-1时刻累积退化的概率密度函数
Figure RE-GDA0003739368750000062
可通过tk-1时刻核密度估计值的k-1重卷积得到，
tk时刻获得新样本数据后，概率密度函数
Figure RE-GDA0003739368750000063
为：
Figure RE-GDA0003739368750000064
tk-1+m时刻获得m个样本数据后，概率密度函数为：
Figure RE-GDA0003739368750000065
同理，在任意时刻获得新的样本数据后，都可获得相应的累积退化概率密度函数；
步骤5、建立剩余寿命预测模型，给定部件i的失效阈值为
Figure RE-GDA0003739368750000066
设部件i在tk-1+t时刻的退化量为
Figure RE-GDA0003739368750000067
时，系统发生故障，那么部件i在tk-1时刻的剩余寿命概率分布函数FT(t)为：
Figure RE-GDA0003739368750000068
由公式(21)可知：
Figure RE-GDA0003739368750000071
将公式(23)代入公式(22)，可得：
Figure RE-GDA0003739368750000072
则tk-1时刻部件i的剩余寿命概率密度函数为：
Figure RE-GDA0003739368750000073
本发明优点及积极效果是：
可自适应地选择出更加准确的窗宽，提高了拟合度，可有效地预测考虑轴承退化影响的齿轮剩余寿命，为后续的维修决策提供依据，
附图说明
图1为实施本发明采用的测试台示意图；
图2为图1中主试齿轮箱GB1与陪试齿轮箱GB2的结构示意图；
图3为本发明实施例中模拟设备在实际运行过程中环境与负载的随机变化状况示意图；
图4为本发明实施例中MSB调制信号双谱分析图；
图5为本发明实施例中部件i的退化状态图；
图6为本发明实施例中随机相关性示意图；
图7为本发明实施例中核密度估计在不同监测时间下剩余寿命预测结果图；
图中：1-驱动电机 2-负载电机 3-控制柜 4-电气测量装置
5.1-第一联轴器 5.2-第二联轴器 5.3-第三联轴器
DAS-16通道数据采集系统 PC-计算机 GB1-主试齿轮箱 GB2-陪试齿轮箱
Z1-主试齿轮箱中第一对齿轮中主动齿轮的齿数
Z2-主试齿轮箱中第一对齿轮中从动齿轮的齿数
Z3-主试齿轮箱中第二对齿轮中传动齿轮的齿数
Z4-主试齿轮箱中第二对齿轮中从动齿轮的齿数
Z5-陪试齿轮箱中第一对齿轮中主动齿轮的齿数
Z6-陪试齿轮箱中第一对齿轮中从动齿轮的齿数
Z7-陪试齿轮箱中第二对齿轮中传动齿轮的齿数
Z8-陪试齿轮箱中第二对齿轮中从动齿轮的齿数
具体实施方式
下面结合附图对本发明实施例做进一步说明：
如图1所示，实施步骤如下：
步骤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所示，
表1齿轮箱规格
Figure RE-GDA0003739368750000091
在试验期间，首先控制驱动电机1以正弦曲线变化的转速运转30min，然后分别在0％负载、25％负载、50％负载、75％负载、100％负载5种不同负载条件下以全速的50％运转，最后同样的分别在5种负载下以全速的70％运行，表2是轴频和啮合频率的公式计算；
表2轴频和啮合频率
Figure RE-GDA0003739368750000092
步骤2、对齿轮的振动数据进行MSB调制信号双谱变换后对齿轮的退化数据进行特征提取并预测其剩余寿命，本实施例试验运行了838h，监测到陪试齿轮箱GB2的振动出现明显峰值后停止试验，打开齿轮箱后，观察到陪试齿轮箱GB2 在低速阶段发生磨损，齿根部位磨损较为严重，齿轮在运行初期样本特征波动较大，为齿轮啮合阶段，使用这个阶段的数据来分析部件的磨损退化过程是不恰当的，因此选择运行300h之后的数据来分析，对振动数据进行调制信号双谱变换，进行特征提取并预测其剩余寿命，常规的双谱B(fx,fc)通过傅里叶变换后在频域内可表示为：
B(fx,fc)＝E[X(fc)X(fx)X*(fc+fx)] (1)
其中X(f)为信号x(t)的离散傅立叶变换；E(g)为求均值符号；fx为调制频率； fc为载波频率；X*为X的复共轭，
对于振动信号x(t)其MSB变换在频域中表示为：
BMS(fx,fc)＝E[X(fc+fx)X(fc-fx)X*(fc)X*(fc)] (2)
MSB对常规双谱改进之后在进行退化特征提取时，同时考虑信号的幅值和相位的变化，能够更好地说明调制信号和载波信号之间的非线性关系，充分表示信号的调制特性，要比常规的双谱分析方法更准确，
以幅值和相位来表示，可将式(2)写为：
Figure RE-GDA0003739368750000101
公式(3)中MSB的总相位可通过以下方式计算：
Figure RE-GDA0003739368750000102
Figure RE-GDA0003739368750000103
和
Figure RE-GDA0003739368750000104
耦合时，相位关系可以表示为
Figure RE-GDA0003739368750000105
Figure RE-GDA0003739368750000106
将公式(5)和公式(6)代入公式(4)，可以得到MSB的总相位为零，它的幅值由其四个分量幅值的乘积所确定，因此，如果(fc+fx)和(fc-fx)来源于fx和 fc的非线性影响，在双频BMS(fx,fc)处会出现明显的双谱峰，这样对MSB的表述将更为准确，与此相反，如果fx和fc无非线性作用，或者它们是随机噪声，其MSB 的相位随时间变化，最终平均结果为零，由此MSB具有很好地抑制噪声作用，图4是采用MSB对试验台数据特征提取的调制双谱分析图，通过MSB双谱分析先从幅值谱发现具有峰值的信号幅值，同时，该坐标点对应的相位谱中也有较大的相干系数，则说明该幅值可以反映部件的退化特征，若相位谱中对应着非常小的相干系数，说明该幅值不是通过调制作用生成，可能是由噪声产生的，从而将其剔除，其在识别非线性影响的同时考虑幅值和相位的影响，从而有效地抑制随机噪声的干扰，提取能够表征部件退化的振动信号特征数据出来，便可得到部件的退化状态曲线图，图2中GB2中的齿轮Z5和GB1中的齿轮Z3通过主轴进行连接，二者之间具有较强的相关性，故将齿轮Z5和齿轮Z3作为部件i和部件j进行研究，图5是部件i的退化状态曲线图；
步骤3、从多部件系统部件间存在的复杂随机相关性的特征以及对部件连续退化状态的不同影响角度出发对部件间的随机相关性进行分析，可将其分为三类：单一单向随机相关性，单一多向随机相关性和双向随机相关性，图6是本实施中部件间随机相关性示意图，
1)单一单向随机相关性：是指某个部件的退化只会对系统中的单个部件的退化产生单向影响；
2)单一多向随机相关性：是指某个部件的退化会对系统中多个部件的退化产生单向影响或者是某个部件的退化会受到系统中其它多个部件退化的单向影响；
3)双向随机相关性：是指部件自身的退化会受到系统中其它部件的影响，同时自身也会对其它部件的退化产生影响，本实施例考虑的是更为复杂的部件间具有双向随机相关性的多部件系统，假设一个系统中的部件i和部件j具有双向随机相关性，即一个部件的退化会加剧另一个部件的退化，反之亦然，本发明考虑更为复杂的部件间具有双向随机相关性的多部件系统进行研究，假设一个系统中的部件i和部件j具有双向随机相关性，即一个部件的退化会加剧另一个部件的退化，反之亦然，随着现代传感器技术的发展，可以通过传感器监测系统中各部件的历史退化状态数据及实时运行数据，利用这些样本数据进行随机相关性建模并建立剩余寿命预测模型，
在考虑更为复杂的部件间具有双向随机相关性的多部件系统进行研究时，对双向随机相关性的模型进行构建，假设部件i在tk-1时刻的退化量为
Figure RE-GDA0003739368750000111
则它在tk时刻的退化量为：
Figure RE-GDA0003739368750000112
为考虑部件之间的随机相关性影响，其部件i在单位时间tk-tk-1的退化增量
Figure RE-GDA0003739368750000121
为：
Figure RE-GDA0003739368750000122
其中，
Figure RE-GDA0003739368750000123
为部件i自身固有的退化增量；
Figure RE-GDA0003739368750000124
为部件j对部件i随机相关影响的退化增量，
同理，假设部件j在时间tk-1的退化量为
Figure RE-GDA0003739368750000125
则它在tk时刻的退化量为：
Figure RE-GDA0003739368750000126
Figure RE-GDA0003739368750000127
其中，
Figure RE-GDA0003739368750000128
为部件j在tk时刻的退化量；
Figure RE-GDA0003739368750000129
为部件j自身固有的退化增量；
Figure RE-GDA00037393687500001210
为部件i对部件j随机相关影响的退化增量，
以部件i为例，设
Figure RE-GDA00037393687500001211
为部件i在[0,tk-1]内的k-1个独立同分布的退化增量，将其概率密度函数记为
Figure RE-GDA00037393687500001212
则
Figure RE-GDA00037393687500001213
的核密度估计为：
Figure RE-GDA00037393687500001214
其中，∑为求和符号；K(g)为核函数；
Figure RE-GDA00037393687500001215
为窗宽；k-1为部件i在tk-1时刻退化增量的样本数；
Figure RE-GDA00037393687500001216
为包含部件i自身固有退化增量及部件j对部件i随机相关影响退化增量的样本值，采用积分均方误差对核估计进行衡量：
Figure RE-GDA00037393687500001217
其中MISE代表求积分均方误差；∫为积分符号，
核函数K(g)和窗宽
Figure RE-GDA00037393687500001218
是核密度估计的两个重要参数，核函数的种类很多，但是不同的核函数对结果的影响很小，选择常用的Gaussian核来建模：
Figure RE-GDA00037393687500001219
窗宽的选择对核密度估计结果的准确性有很大影响，通过求解MISE的最小值可得到窗宽
Figure RE-GDA00037393687500001220
的最优解：
Figure RE-GDA0003739368750000131
其中，σk-1为k-1个已知样本的标准差；n是样本的数量，
将公式(11)代入上式，可以求出：
Figure RE-GDA0003739368750000132
同理，部件i在运行过程中自身固有的退化增量为
Figure RE-GDA0003739368750000133
其核密度估计
Figure RE-GDA0003739368750000134
可表示为：
Figure RE-GDA0003739368750000135
部件j在此系统运行期间的退化增量为
Figure RE-GDA0003739368750000136
其核密度估计
Figure RE-GDA0003739368750000137
可表示为：
Figure RE-GDA0003739368750000138
部件j对部件i随机相关影响的退化增量记为
Figure RE-GDA0003739368750000139
将其服从的概率密度函数记为fk-1(Δxji)，由式(9)知，
Figure RE-GDA00037393687500001310
为考虑部件i自身固有退化及部件j对部件i随机相关影响的概率密度函数，则fk-1(Δxji)的核密度估计为：
Figure RE-GDA00037393687500001311
步骤4、引入窗宽因子λn采用自适应窗宽核密度估计更能反映样本数据的真实情况，当
Figure RE-GDA00037393687500001312
时，核密度估计与真实值之间的误差很小，则tk时刻获得新样本数据后的自适应窗宽
Figure RE-GDA00037393687500001313
可表示为：
Figure RE-GDA00037393687500001314
其中，
Figure RE-GDA00037393687500001315
为tk时刻新的自适应窗宽；λk为窗宽因子；
Figure RE-GDA00037393687500001316
为部件i在窗宽
Figure RE-GDA00037393687500001317
下tk时刻的核密度估计值，
部件i在tk时刻基于自适应窗宽的核密度估计为：
Figure RE-GDA0003739368750000141
其中，“*”为卷积符号；
Figure RE-GDA0003739368750000142
为部件i在tk-1时刻的核密度估计值；
Figure RE-GDA0003739368750000143
为部件i在tk时刻自身退化影响的核密度估计值；
Figure RE-GDA0003739368750000144
为部件j对部件i随机相关影响的核密度估计值，采用递推来实现核密度估计的实时更新，下一时刻有新增样本时，通过退化增量便可更新计算，大大减少了计算量，降低了计算复杂度，
对退化分布进行计算，得到部件i基于自适应窗宽的核密度估计表达式后，可以通过tk-1时刻的核密度估计
Figure RE-GDA0003739368750000145
来求部件i在tk-1时刻累积退化的概率密度函数，记为
Figure RE-GDA0003739368750000146
Figure RE-GDA0003739368750000147
部件i在tk-1时刻累积退化的概率密度函数
Figure RE-GDA0003739368750000148
可通过tk-1时刻核密度估计值的k-1重卷积得到，
tk时刻获得新样本数据后，概率密度函数
Figure RE-GDA0003739368750000149
为：
Figure RE-GDA0003739368750000151
tk-1+m时刻获得m个样本数据后，概率密度函数为：
Figure RE-GDA0003739368750000152
同理，在任意时刻获得新的样本数据后，都可获得相应的累积退化概率密度函数；
步骤5、建立剩余寿命预测模型，给定部件i的失效阈值为
Figure RE-GDA0003739368750000153
设部件i在tk-1+t时刻的退化量为
Figure RE-GDA0003739368750000154
时，系统发生故障，那么部件i在tk-1时刻的剩余寿命概率分布函数FT(t)为：
Figure RE-GDA0003739368750000155
由公式(21)可知：
Figure RE-GDA0003739368750000156
将公式(23)代入公式(22)，可得：
Figure RE-GDA0003739368750000161
则tk-1时刻部件i的剩余寿命概率密度函数为：
Figure RE-GDA0003739368750000162
将MSB特征提取后的退化数据代入模型，可根据式22-式25求得部件i在当前监测时刻的剩余寿命概率密度函数，图7是本次试验在7个监测时间下的剩余寿命预测结果，可以看出，随着监测时间的增加，样本数据增大，本发明提出方法的剩余寿命预测结果越来越接近真实值，说明本发明提出的基于自适应核窗宽的核密度估计剩余寿命预测方法可以很好的对部件i剩余寿命的概率密度函数进行估计，
表3给出了参数估计和核密度估计方法对部件剩余寿命预测结果的误差分析，从表中可以看出，随着监测时间的推移，获得的样本数据越多，两种方法对部件剩余寿命的预测就越精确，在同一监测时间下，本发明提出的核密度估计的方法要比参数估计方法更精确，进一步验证了非参数核密度估计方法的有效性和精确性，
表3参数估计和核估计剩余寿命预测误差分析
Figure RE-GDA0003739368750000163
综上所述，本发明针对多部件系统部件连续退化过程中存在的双向随机相关性，基于自适应窗宽的核密度估计方法建立了实时剩余寿命预测模型，首先针对部件间存在的双向随机相关性影响，引入条件核密度估计来进行建模，然后采用自适应窗宽的核密度估计方法求得相应的概率密度函数，最后建立实时剩余寿命预测模型，得到部件的剩余寿命。

Claims (1)

1.一种齿轮箱中考虑轴承退化影响的齿轮剩余寿命预测方法，其特征在于实施步骤如下：
步骤1、通过试验获取表征斜齿轮箱内的实时监测数据，采用斜齿轮箱试验台进行试验，其上安装的主试齿轮箱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；
在试验期间，首先控制驱动电机(1)以正弦曲线变化的转速运转30min，然后分别在0％负载、25％负载、50％负载、75％负载、100％负载5种不同负载条件下以全速的50％运转，最后同样的分别在5种负载下以全速的70％运行；
步骤2、对齿轮的振动数据进行MSB调制信号双谱变换后对齿轮的退化数据进行特征提取并预测其剩余寿命，常规的双谱B(fx,fc)通过傅里叶变换后在频域内可表示为：
B(fx,fc)＝E[X(fc)X(fx)X*(fc+fx)] (1)
其中X(f)为振动信号x(t)的离散傅立叶变换；E(g)为求均值符号；fx为调制频率；fc为载波频率；X*为X的复共轭；
对于振动信号x(t)其MSB变换在频域中表示为：
BMS(fx,fc)＝E[X(fc+fx)X(fc-fx)X*(fc)X*(fc)] (2)
MSB对常规双谱改进之后在进行退化特征提取时，同时考虑信号的幅值和相位的变化，能够更好地说明调制信号和载波信号之间的非线性关系，充分表示信号的调制特性，要比常规的双谱分析方法更准确；
以幅值和相位来表示，可将式(2)写为：
Figure FDA0003560913300000021
公式(3)中MSB的总相位可通过以下方式计算：
Figure FDA0003560913300000022
Figure FDA0003560913300000023
和
Figure FDA0003560913300000024
耦合时，相位关系可以表示为：
Figure FDA0003560913300000025
Figure FDA0003560913300000026
步骤3、从多部件系统部件间存在的复杂随机相关性的特征以及对部件连续退化状态的不同影响角度出发对部件间的随机相关性进行分析，可将其分为三类：单一单向随机相关性，单一多向随机相关性和双向随机相关性，
1)单一单向随机相关性：是指某个部件的退化只会对系统中的单个部件的退化产生单向影响；
2)单一多向随机相关性：是指某个部件的退化会对系统中多个部件的退化产生单向影响或者是某个部件的退化会受到系统中其它多个部件退化的单向影响；
3)双向随机相关性：是指部件自身的退化会受到系统中其它部件的影响，同时自身也会对其它部件的退化产生影响，本发明考虑的是更为复杂的部件间具有双向随机相关性的多部件系统，假设一个系统中的部件i和部件j具有双向随机相关性，即一个部件的退化会加剧另一个部件的退化，反之亦然；
在考虑更为复杂的部件间具有双向随机相关性的多部件系统进行研究时，对双向随机相关性的模型进行构建，假设部件i在tk-1时刻的退化量为
Figure FDA0003560913300000031
则它在tk时刻的退化量为：
Figure FDA0003560913300000032
其中
Figure FDA0003560913300000033
为部件i在[0,tk-1]内的k-1个独立同分布的退化增量，假设部件j在时间tk-1的退化量为
Figure FDA0003560913300000034
则它在tk时刻的退化量为：
Figure FDA0003560913300000035
Figure FDA0003560913300000036
其中，
Figure FDA0003560913300000037
为部件j在tk时刻的退化量；
Figure FDA0003560913300000038
为部件j自身固有的退化增量；
Figure FDA0003560913300000039
为部件i对部件j随机相关影响的退化增量，
将其概率密度函数记为
Figure FDA00035609133000000310
则
Figure FDA00035609133000000311
的核密度估计为：
Figure FDA00035609133000000312
式中，∑为求和符号；K(g)为核函数；
Figure FDA00035609133000000313
为窗宽；k-1为部件i在tk-1时刻退化增量的样本数；
Figure FDA00035609133000000314
为包含部件i自身固有退化增量及部件j对部件i随机相关影响退化增量的样本值，
采用积分均方误差对核估计进行衡量：
Figure FDA00035609133000000315
其中MISE代表求积分均方误差；∫为积分符号，
不同的核函数对结果的影响很小，选择常用的Gaussian核来建模：
Figure FDA00035609133000000316
通过求解MISE的最小值可得到窗宽
Figure FDA0003560913300000041
的最优解：
Figure FDA0003560913300000042
其中，σk-1为k-1个已知样本的标准差；n是样本的数量，
将公式(11)代入上式，可以求出：
Figure FDA0003560913300000043
同理，部件i在运行过程中自身固有的退化增量为
Figure FDA0003560913300000044
其核密度估计
Figure FDA0003560913300000045
可表示为：
Figure FDA0003560913300000046
部件j在此系统运行期间的退化增量为
Figure FDA0003560913300000047
其核密度估计
Figure FDA0003560913300000048
可表示为：
Figure FDA0003560913300000049
部件j对部件i随机相关影响的退化增量记为
Figure FDA00035609133000000410
将其服从的概率密度函数记为fk-1(Δxji)，由式(9)知，
Figure FDA00035609133000000411
为考虑部件i自身固有退化及部件j对部件i随机相关影响的概率密度函数，则fk-1(Δxji)的核密度估计为：
Figure FDA00035609133000000412
步骤4、引入窗宽因子λn采用自适应窗宽核密度估计更能反映样本数据的真实情况，当
Figure FDA00035609133000000413
时，核密度估计与真实值之间的误差很小，则tk时刻获得新样本数据后的自适应窗宽
Figure FDA00035609133000000414
可表示为：
Figure FDA00035609133000000415
其中，
Figure FDA0003560913300000051
为tk时刻新的自适应窗宽；λk为窗宽因子；
Figure FDA0003560913300000052
为部件i在窗宽
Figure FDA0003560913300000053
下tk时刻的核密度估计值，
部件i在tk时刻基于自适应窗宽的核密度估计为：
Figure FDA0003560913300000054
其中，“*”为卷积符号；
Figure FDA0003560913300000055
为部件i在tk-1时刻的核密度估计值；
Figure FDA0003560913300000056
为部件i在tk时刻自身退化影响的核密度估计值；
Figure FDA0003560913300000057
为部件j对部件i随机相关影响的核密度估计值，
对退化分布进行计算，得到部件i基于自适应窗宽的核密度估计表达式后，可以通过tk-1时刻的核密度估计
Figure FDA0003560913300000058
来求部件i在tk-1时刻累积退化的概率密度函数，记为
Figure FDA0003560913300000059
Figure FDA00035609133000000510
部件i在tk-1时刻累积退化的概率密度函数
Figure FDA00035609133000000511
可通过tk-1时刻核密度估计值的k-1重卷积得到，
tk时刻获得新样本数据后，概率密度函数
Figure FDA00035609133000000512
为：
Figure FDA0003560913300000061
tk-1+m时刻获得m个样本数据后，概率密度函数为：
Figure FDA0003560913300000062
同理，在任意时刻获得新的样本数据后，都可获得相应的累积退化概率密度函数；
步骤5、建立剩余寿命预测模型，给定部件i的失效阈值为
Figure FDA0003560913300000063
设部件i在tk-1+t时刻的退化量为
Figure FDA0003560913300000064
时，系统发生故障，那么部件i在tk-1时刻的剩余寿命概率分布函数FT(t)为：
Figure FDA0003560913300000065
由公式(21)可知：
Figure FDA0003560913300000066
将公式(23)代入公式(22)，可得：
Figure FDA0003560913300000071
则tk-1时刻部件i的剩余寿命概率密度函数为：
Figure FDA0003560913300000072
CN202210292839.1A 2022-03-23 2022-03-23 齿轮箱中考虑轴承退化影响的齿轮剩余寿命预测方法 Pending CN114942139A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210292839.1A CN114942139A (zh) 2022-03-23 2022-03-23 齿轮箱中考虑轴承退化影响的齿轮剩余寿命预测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210292839.1A CN114942139A (zh) 2022-03-23 2022-03-23 齿轮箱中考虑轴承退化影响的齿轮剩余寿命预测方法

Publications (1)

Publication Number Publication Date
CN114942139A true CN114942139A (zh) 2022-08-26

Family

ID=82905782

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210292839.1A Pending CN114942139A (zh) 2022-03-23 2022-03-23 齿轮箱中考虑轴承退化影响的齿轮剩余寿命预测方法

Country Status (1)

Country Link
CN (1) CN114942139A (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117423175A (zh) * 2023-10-25 2024-01-19 深圳泰瑞谷科技有限公司 一种车辆诊断数据展示方法、装置、诊断仪及介质
WO2024063693A1 (en) * 2022-09-20 2024-03-28 Hitachi, Ltd. Method and system for remaining useful life prediction of a multi-component operational system
CN117875191A (zh) * 2024-03-08 2024-04-12 东莞市星火齿轮有限公司 一种基于大数据的齿轮箱寿命评估方法及系统、存储介质
CN117930003A (zh) * 2024-03-21 2024-04-26 南昌三瑞智能科技股份有限公司 无人机电机使用寿命的测试方法、处理装置及测试设备

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2024063693A1 (en) * 2022-09-20 2024-03-28 Hitachi, Ltd. Method and system for remaining useful life prediction of a multi-component operational system
CN117423175A (zh) * 2023-10-25 2024-01-19 深圳泰瑞谷科技有限公司 一种车辆诊断数据展示方法、装置、诊断仪及介质
CN117875191A (zh) * 2024-03-08 2024-04-12 东莞市星火齿轮有限公司 一种基于大数据的齿轮箱寿命评估方法及系统、存储介质
CN117875191B (zh) * 2024-03-08 2024-05-28 东莞市星火齿轮有限公司 一种基于大数据的齿轮箱寿命评估方法及系统、存储介质
CN117930003A (zh) * 2024-03-21 2024-04-26 南昌三瑞智能科技股份有限公司 无人机电机使用寿命的测试方法、处理装置及测试设备
CN117930003B (zh) * 2024-03-21 2024-05-31 南昌三瑞智能科技股份有限公司 无人机电机使用寿命的测试方法、处理装置及测试设备

Similar Documents

Publication Publication Date Title
CN114942139A (zh) 齿轮箱中考虑轴承退化影响的齿轮剩余寿命预测方法
Wang et al. Unknown fault feature extraction of rolling bearings under variable speed conditions based on statistical complexity measures
Zappalá et al. Side‐band algorithm for automatic wind turbine gearbox fault detection and diagnosis
de Azevedo et al. A review of wind turbine bearing condition monitoring: State of the art and challenges
El-Thalji et al. A summary of fault modelling and predictive health monitoring of rolling element bearings
CN111259717B (zh) 一种旋转设备状态异常判断方法及系统
CN105181019B (zh) 旋转类机械早期故障预警分析方法
Feng et al. A novel similarity-based status characterization methodology for gear surface wear propagation monitoring
CN106197996A (zh) 基于多元数据的海上起重机齿轮箱故障诊断装置及方法
Zhe et al. Pitting damage levels estimation for planetary gear sets based on model simulation and grey relational analysis
CN110174281B (zh) 一种机电设备故障诊断方法及系统
CN113947017A (zh) 一种滚动轴承剩余使用寿命预测方法
Cocconcelli et al. An algorithm to diagnose ball bearing faults in servomotors running arbitrary motion profiles
Feng et al. A novel cyclic-correntropy based indicator for gear wear monitoring
CN109883691A (zh) 核估计和随机滤波集成的齿轮剩余寿命预测方法
Chorna et al. Identification of changes in the parameters of induction motors during monitoring by measuring the induction of a magnetic field on the stator surface
CN109596349A (zh) 一种基于vmd和pct的减速器故障诊断方法
Shang et al. An intelligent fault diagnosis system for newly assembled transmission
Cao et al. Deterioration state diagnosis and wear evolution evaluation of planetary gearbox using vibration and wear debris analysis
CN110398362B (zh) 一种机器人rv减速器故障诊断和定位方法
Yao et al. Fault detection of complex planetary gearbox using acoustic signals
Zhang et al. Generalized transmissibility damage indicator with application to wind turbine component condition monitoring
CN113204849A (zh) 一种齿轮箱齿轮剥落故障检测方法
Gildish et al. Vibration-Based Estimation of Gearbox Operating Conditions: Machine Learning Approach
CN109580218B (zh) 一种基于似然学习机的风机齿轮箱状态识别方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination