CN110174261B - 多退化量监测的齿轮实时剩余寿命预测方法 - Google Patents
多退化量监测的齿轮实时剩余寿命预测方法 Download PDFInfo
- Publication number
- CN110174261B CN110174261B CN201910432767.4A CN201910432767A CN110174261B CN 110174261 B CN110174261 B CN 110174261B CN 201910432767 A CN201910432767 A CN 201910432767A CN 110174261 B CN110174261 B CN 110174261B
- Authority
- CN
- China
- Prior art keywords
- function
- copula
- degradation
- gear
- gearbox
- 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.)
- Active
Links
- 238000006731 degradation reaction Methods 0.000 title claims abstract description 74
- 238000000034 method Methods 0.000 title claims abstract description 38
- 238000012544 monitoring process Methods 0.000 title claims abstract description 21
- 238000012360 testing method Methods 0.000 claims abstract description 93
- 241000039077 Copula Species 0.000 claims abstract description 64
- 230000015556 catabolic process Effects 0.000 claims abstract description 53
- 230000001133 acceleration Effects 0.000 claims abstract description 47
- 238000005315 distribution function Methods 0.000 claims abstract description 47
- 238000011156 evaluation Methods 0.000 claims abstract description 7
- 238000001914 filtration Methods 0.000 claims abstract description 7
- 238000005070 sampling Methods 0.000 claims description 33
- 238000009499 grossing Methods 0.000 claims description 9
- 230000008878 coupling Effects 0.000 claims description 6
- 238000010168 coupling process Methods 0.000 claims description 6
- 238000005859 coupling reaction Methods 0.000 claims description 6
- 230000000694 effects Effects 0.000 claims description 6
- 230000007850 degeneration Effects 0.000 claims description 4
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 3
- 238000007476 Maximum Likelihood Methods 0.000 claims description 3
- 238000004364 calculation method Methods 0.000 claims description 3
- 238000010219 correlation analysis Methods 0.000 claims description 3
- 238000012216 screening Methods 0.000 claims description 3
- 238000005299 abrasion Methods 0.000 claims description 2
- 238000004458 analytical method Methods 0.000 claims 1
- 238000012545 processing Methods 0.000 claims 1
- 238000012423 maintenance Methods 0.000 abstract description 6
- 238000000605 extraction Methods 0.000 abstract 1
- 230000003449 preventive effect Effects 0.000 abstract 1
- 230000008569 process Effects 0.000 description 3
- 230000008859 change Effects 0.000 description 2
- 238000004026 adhesive bonding Methods 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 230000036541 health Effects 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 230000008092 positive effect Effects 0.000 description 1
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/02—Gearings; Transmission mechanisms
- G01M13/021—Gearings
-
- 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/02—Gearings; Transmission mechanisms
- G01M13/025—Test-benches with rotational drive means and loading means; Load or drive simulation
-
- 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/02—Gearings; Transmission mechanisms
- G01M13/028—Acoustic or vibration analysis
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Acoustics & Sound (AREA)
- Testing Of Devices, Machine Parts, Or Other Structures Thereof (AREA)
Abstract
一种多退化量监测的齿轮实时剩余寿命预测方法，属于机械可靠性技术领域。实施步骤如下：1、利用加速度传感器和噪声传感器对主试齿轮箱内齿轮退化实时监测；2、对齿轮退化状态进行特征提取和衰退评估；3、采用核估计和随机滤波理论的方法分别对齿轮箱的振动加速度和噪声进行建模，获得齿轮箱剩余寿命概率密度函数，得到单退化量剩余寿命边缘分布函数；4、利用Copula函数表示齿轮箱的振动加速度和噪声之间的随机相关性，求得齿轮箱剩余寿命的联合分布函数；5、根据齿轮箱剩余寿命联合分布函数求得其剩余寿命联合概率密度函数，最后得到齿轮箱剩余寿命预测值；优点是有效地预测齿轮退化状态及实时剩余寿命，为齿轮预防性维修提供依据。
Description
技术领域
本发明属于机械可靠性领域，具体涉及一种多退化量监测的齿轮实时剩余寿命预测方法。
背景技术
齿轮是机械工业中应用广泛的传动系统中的关键部件；当齿轮发生断齿、齿面疲劳、胶合等故障时，常常会引起整个机械设备灾难性的破坏，以风力发电机组为例，齿轮故障率是整个风电机组中最高的，约占60％，而且其维护费用也较高，约占40％，因此，对齿轮提出合理有效的维修方案已成为风电行业急需解决的问题，而在整个维修方案制定过程中，齿轮的剩余寿命预测是重点和难点，随着信息传感设备的发展，可对齿轮的运行状态进行实时监测，利用接收到的大量实时监测信息更为准确预测系统的退化状态及其剩余寿命，能够提供有关健康状态的关键信息，进而识别和管理故障的发生、规划维修活动，为更合理地制定基于状态的维护维修策略提供依据。
目前，齿轮剩余寿命的预测方法分为四类：基于物理模型的预测方法、基于统计经验的预测方法、基于知识的预测方法和基于数据驱动的预测方法；现有预测方法存在下列问题：首先，现有预测方法需要进行状态退化模型结构假设，需要假设其作为判断依据的样本符合某种特定的模型结构，这些模型结构的假设与实际的物理模型之间常常存在较大的差距；其次，预测模型中所涉及到的参数估计问题，大多不能保证全局收敛；最后，由于齿轮处于变化的环境中，它的状态退化模型会发生改变，单一的预测模型不能适应环境的变化，需要多种预测模型相结合。
发明内容
本发明目的是提供一种多退化量监测的齿轮实时剩余寿命预测方法，可有效地克服了一些现有技术中存在的缺点，所述多退化量是指2～4种退化量。
本发明是这样实现的，其特征是包括以下步骤：
步骤1、通过试验获取表征主试齿轮箱1内齿轮状态的实时监测数据：
采用如图2所示的试验台架，试验台架主、陪试齿轮箱1、2的中心距a＝150mm；试验采用机械杠杆6加载，扭矩采用转矩传感器13#进行测量；主陪试齿轮箱1、2中均为正反交错搭接的一对齿轮，齿轮的断齿状态等效为齿轮的失效；
试验共布置八个加速度传感器1#～8#,两个噪声传感器9#和10#，一个温度传感器11#,一个转速传感器12#，一个转矩传感器13#；有四个加速度传感器1#～4#分别布置在主试齿轮箱1轴承座的径向位置，有两个加速度传感器7#和8#分别布置在主试齿轮箱1的轴向位置，有两个加速度传感器5#和6#分别布置在陪试齿轮箱2轴承座的径向位置；有两个噪声传感器9#和10#分别悬挂在主试齿轮箱1和陪试齿轮箱2正上方40cm处；一个温度传感器11#布置在主试齿轮箱1内部；一个转速传感器12#布置在驱动电机8与陪试齿轮箱2联接轴的中部；一个转矩传感器13#布置在主试齿轮箱1与陪试齿轮箱2联接轴的中部；试验过程中加载了八级载荷；八级载荷的大小在330～850扭矩之间，每级载荷的加载时间为9～12个小时，在第八级载荷时发生断齿，通过用加速度传感器4#记录齿轮的加速度数据，剩余寿命预测选取是从第八级加载开始到断齿的加速度测试点整个时域信号进行分析；采样信息如下：采样频率为20～30kHz，每次采样持续50～70秒，每隔8～10分钟记录一次采样文本；
步骤2、对主试齿轮箱1内齿轮的退化状态进行特征提取，利用均方幅值对齿轮磨损退化性能进行衰退评估，对于每次采样时间长度内，采样信号的均方幅值特征值表示为：
式中：∑为求和号，n为每个采样周期的采样点数,yi为ti时刻齿轮的状态信息均方幅值，yj为齿轮每个采样周期的数据；
步骤3、根据主试齿轮箱1内齿轮的初始数据，利用核密度估计方法以及随机滤波理论，建立单退化量的退化轨迹，得到主试齿轮箱1内齿轮的剩余寿命概率密度函数；
对于主试齿轮箱1的剩余寿命概率密度函数fi,k(xi,k|Yi,k)，通过贝叶斯定理获得fi,k(xi,k|Yi,k)的递推形式；
由贝叶斯定理得：
式中：xi,k表示在单退化量i下tk时刻主试齿轮箱1的剩余寿命；Yi,k表示第i个退化量tk时刻的历史监测数据，i＝1表示单退化量为加速度，i＝2表示单退化量为噪声；
因为：
fi,k(yi,k|xi,k,Yi,k-1)＝fi,k(yi,k|xi,k) (3)
所以：
在tk时刻的剩余寿命等于tk-1时刻的剩余寿命减去时刻tk和时刻tk-1之间的间隔，即：
根据式(5)和式(6)，得到随机滤波方程为：
在t1时刻fi,0(xi,0-t1+t0|Yi,0)＝fi,0(xi,0-t1+t0)，可求得：
在t2时刻，求得：
估计出fi,0(xi,0)和fi,k(yi,k|xi,k)，通过递推公式(7)计算得到fi,k(xi,k|Yi,k)，
用核密度估计fi,0(xi,0)和fi,k(yi,k|xi,k)，即：
设X1，X2，…，Xn为独立同分布的表征退化的随机变量，其概率密度函数为f(x)，则：
式中：h为平滑参数,K(u)为核函数，n为随机变量X的样本数；
由式(10)得：
式中：
(
为样本方差)为平滑参数，
为核函数，ni为随机变量Xi的样本数；
式中：αi、βi为平滑参数,K1(a,b)、K2(u)为核函数；
根据核密度估计结果以及式(7)，得到主试齿轮箱1在时间tk处基于加速度的剩余寿命概率密度函数为f1,k(x1,k|Y1,k)和基于噪声的剩余寿命概率密度函数为f2,k(x2,k|Y2,k)；
根据概率论中概率密度函数与边缘分布函数的关系，基于加速度的剩余寿命的边缘分布函数为
和基于噪声的剩余寿命的边缘分布函数为
步骤4、齿轮的单退化量之间的关系是不知道的，Copula函数是一种相关性分析方法，将单退化量的边缘分布函数与其联合分布函数联系起来，即：
F(t1,t2)＝C(F1(t1),F2(t2)；γ)＝C(u1,u2；γ) (13)
式中：F(t1,t2)为变量t1和t2的联合分布函数；u1＝F1(t1)和u2＝F2(t2)分别为变量t1和t2的边缘分布函数；γ为Copula函数的参数。
变量t1和t2的边缘分布函数F1(t1)和F2(t2)已知，Copula函数的未知参数γ采用极大似然估计法估计：
式中：c为C的密度函数，∑为求和号，n为每个采样周期的采样点数；
在工程实践中，Copula常用的3种函数形式，如表1所示。
表1常用的Copula函数形式
在同一工程实践中选择不同形式的Copula函数，出现不同的结果。经过初步筛选，选择符合尾部相关性的Copula函数。常用的3种Copula函数中不同参数γ下尾部的非对称性不明显，通过模型评价的方法选择最优的Copula函数，AIC准则是最常用的模型评价准则，具有良好的适用性和实用性，即：
式中：L为Copula函数的似然值；n为样本个数；K为被估计参数个数。
常用Copula函数的似然函数如下：
Gumbel函数的似然函数：
Clayton函数的似然函数：
Frank函数的似然函数：
AIC表示Copula函数对数据的拟合效果，AIC值越小，对数据的拟合效果越好。
在求得主试齿轮箱1的加速度和噪声作为退化量下剩余寿命的边缘分布后，分别通过Gumbel Copula、ClaytonCopula、Frank Copula这三种常用的Copula函数对基于振动加速度和噪声的剩余寿命边缘分布函数进行融合，得到主试齿轮箱1的剩余寿命联合分布，根据AIC准则选择最优的Copula函数；
利用所选的最优的Copula函数的性质，融合基于加速度和噪声的边缘分布函数，得到主试齿轮箱1的剩余寿命联合分布函数F(t1,t2)＝C(F1(t1),F2(t2)；γ)，根据概率论中概率密度函数与边缘分布函数的关系，得到主试齿轮箱1基于多退化量下的剩余寿命概率密度函数
步骤5、通过基于Copula函数的多退化量的寿命预测模型预测齿轮剩余寿命；最后得出齿轮的平均剩余寿命为：
本发明优点及积极效果是：
本发明应用到齿轮的剩余寿命预测中，建立预测模型，能够对齿轮剩余寿命做出有效预测；
附图说明
图1为本发明实施例中齿轮剩余寿命预测方法流程图；
图2为本发明实施例中试验台架示意图；
图3为本发明实施例中主试齿轮箱1加速度监测数据；
图4为本发明实施例中主试齿轮箱1噪声监测数据；
图5为本发明实施例中第40小时剩余寿命概率密度；
图6为本发明实施例中第70小时剩余寿命概率密度；
图7为本发明实施例中第40小时剩余寿命的边缘分布函数；
图8为本发明实施例中第70小时剩余寿命的边缘分布函数；
图9为本发明实施例中第40小时剩余寿命的联合分布函数；
图10为本发明实施例中第70小时剩余寿命的联合分布函数；
图11为本发明实施例中第40小时剩余寿命概率密度函数；
图12为本发明实施例中第70小时剩余寿命概率密度函数；
图13为本发明实施例中剩余寿命估计；
图中：1—主试齿轮箱，2—陪试齿轮箱，3—第一联轴器，4—第二联轴器，5—第三联轴器，6—机械杠杆，7—控制系统，8—驱动电机，a—主试齿轮箱1与陪试齿轮箱2的中心距；
1#～8#—加速度传感器，9#、10#—噪声传感器，11#—温度，12#—转速传感器，13#—转矩传感器；
具体实施方式
下面结合附图对本发明实施例做进一步说明：
本发明实施例中，基于多退化量监测的齿轮实时剩余寿命预测方法，方法流程图如图1所示，包括以下步骤：
步骤1、通过试验台架获取表征主试齿轮箱1内齿轮状态的实时监测数据：
采用如图2所示的试验台架，试验台架的中心距a＝150mm；试验采用机械杠杆6加载；主陪试齿轮箱1、2中均为正反交错搭接的一对齿轮，齿轮的断齿状态等效为齿轮的失效；
试验共布置13个传感器，如图2所示，其中1#～8#为加速度传感器，9#和10#为噪声传感器，11#为温度传感器,12#为转速传感器，13#为转矩传感器。1#～4#加速度传感器布置在主试齿轮箱1轴承座的径向位置，7#和8#加速度传感器布置在主试齿轮箱1的轴向位置，5#和6#加速度传感器布置在陪试齿轮箱2轴承座的径向位置；9#和10#噪声传感器分别悬挂在主试齿轮箱1和陪试齿轮箱2正上方40cm处位置；11#温度传感器布置在主试齿轮箱1内部；12#转速传感器布置在驱动电机8与陪试齿轮箱2联接轴的中部；13#转矩传感器布置在主试齿轮箱1与陪试齿轮箱2联接轴的中部；试验过程中对齿轮箱加载了八级载荷，八级载荷的大小分别为：349.5扭矩，430.7扭矩，492.2扭矩，555.6扭矩，612.9扭矩，693.4扭矩，734扭矩，822.7扭矩，每级载荷的时间为10个小时，在第八级载荷时发生断齿，用加速度传感器4#记录齿轮的加速度数据，剩余寿命预测选取是从第八级加载开始到断齿的加速度测试点整个时域信号进行分析；采样信息如下：采样频率为25.6kHz，每次采样持续60秒，每隔9分钟记录一次采样文本；
步骤2、对齿轮的退化状态进行特征提取，利用均方幅值对齿轮磨损退化性能进行衰退评估，对于每次采样时间长度内，采样信号的均方幅值特征值表示为：
式中：∑为求和号，n为每个采样周期的采样点数,yi为ti时刻齿轮的状态信息，yj为齿轮每个采样周期的数据；
通过均方幅值法对加速度和噪声数据的处理，得到如图3和图4所示的主试齿轮箱1退化量随监测时间变化曲线图；由图3和图4可知主试齿轮箱1在发生故障时加速度的阈值为y1＝76.325mm/s2，噪声的阈值为y2＝6.4661mA；
步骤3、基于多退化量监测的齿轮实时剩余寿命预测方法在单退化量下预测主试齿轮箱1的剩余寿命，得到如图5和图6所示的主试齿轮箱1内齿轮的剩余寿命概率密度函数，得到如图7和图8所示的边缘分布函数；
对于主试齿轮箱1的剩余寿命概率密度函数fi,k(xi,k|Yi,k)，通过贝叶斯定理获得fi,k(xi,k|Yi,k)的递推形式；
由贝叶斯定理得：
式中：xi,k表示在单退化量i下tk时刻主试齿轮箱1的剩余寿命；Yi,k表示第i个退化量tk时刻的历史监测数据，i＝1表示单退化量为加速度，i＝2表示单退化量为噪声；
因为：
fi,k(yi,k|xi,k,Yi,k-1)＝fi,k(yi,k|xi,k) (3)
所以：
在tk时刻的剩余寿命等于tk-1时刻的剩余寿命减去时刻tk和时刻tk-1之间的间隔，即：
根据式(5)和式(6)，得到随机滤波方程为：
在t1时刻，可知fi,0(xi,0-t1+t0|Yi,0)＝fi,0(xi,0-t1+t0)，可求得：
在t2时刻，可求得：
估计出fi,0(xi,0)和fi,k(yi,k|xi,k)，通过递推公式(7)计算得到fi,k(xi,k|Yi,k)；
用核密度估计fi,0(xi,0)和fi,k(yi,k|xi,k)，即：
设X1，X2，…，Xn为独立同分布的表征退化的随机变量，其概率密度函数为f(x)，则：
式中：h为平滑参数,K(u)为核函数，n为随机变量X的样本数；
由式(10)得：
式中：
(
为样本方差)为平滑参数，
为核函数，ni为随机变量Xi的样本数；
式中：αi、βi为平滑参数,K1(a,b)、K2(u)为核函数；
根据核密度估计结果以及式(7)，得到主试齿轮箱1在时间tk处基于加速度的剩余寿命概率密度函数为f1,k(x1,k|Y1,k)和基于噪声的剩余寿命概率密度函数为f2,k(x2,k|Y2,k)；
根据概率论中概率密度函数与边缘分布函数的关系，基于加速度的剩余寿命的边缘分布函数为
和基于噪声的剩余寿命的边缘分布函数为
步骤4、齿轮的单退化量之间的关系是不知道的，Copula函数是一种相关性分析方法，将单退化量的边缘分布函数与其联合分布函数联系起来，即：
F(t1,t2)＝C(F1(t1),F2(t2)；γ)＝C(u1,u2；γ) (13)
式中：F(t1,t2)为变量t1和t2的联合分布函数；u1＝F1(t1)和u2＝F2(t2)分别为变量t1和t2的边缘分布函数；γ为Copula函数的参数。
变量t1和t2的边缘分布函数F1(t1)和F2(t2)已知，Copula函数的未知参数γ采用极大似然估计法估计：
式中：c为C的密度函数，∑为求和号，n为每个采样周期的采样点数；
在工程实践中，Copula常用3种函数形式，如表1所示。
表1常用的Copula函数形式
在同一工程实践中选择不同形式的Copula函数，出现不同的结果。经过初步筛选，选择符合尾部相关性的Copula函数。常用的3个Copula函数中不同参数γ下尾部的非对称性不明显，通过模型评价的方法选择合适的Copula函数。AIC准则是最常用的模型评价准则，具有良好的适用性和实用性，即：
式中：L为Copula函数的似然值；n为样本个数；K为被估计参数个数。
常用Copula函数的似然函数如下：
Gumbel函数的似然函数：
Clayton函数的似然函数：
Frank函数的似然函数：
AIC表示Copula函数对数据的拟合效果。即：AIC值越小，对数据的拟合效果越好。
在求得主试齿轮箱1的加速度和噪声作为退化量下剩余寿命的边缘分布后，分别通过Gumbel Copula、Clayton Copula、Frank Copula这三种常用的Copula函数对基于振动加速度和噪声的剩余寿命边缘分布函数进行融合，得到齿轮箱的剩余寿命联合分布；由图9和图10可知，3种常用的Copula函数得到的齿轮箱剩余寿命联合分布函数比较接近，难以判断哪种Copula函数所得结果更接近齿轮箱剩余寿命的真实分布。根据AIC准则选择最优的Copula函数。
由式(15)得到如表2所示的3种常用Copula函数的AIC值。
表2常用Copula函数的AIC值
从表2中可以看出在40小时和70小时时刻，Clayton Copula函数具有最小的AIC值。因此，Clayton Copula函数是最优的Copula函数，可以最好的描述齿轮箱的加速度与噪声之间的相关性。
利用Clayton Copula函数的性质，融合基于加速度和噪声的边缘分布函数，得到主试齿轮箱1的剩余寿命联合分布函数F(t1,t2)＝C(F1(t1),F2(t2)；γ)，根据概率论中概率密度函数与边缘分布函数的关系，得到主试齿轮箱1基于多退化量下的剩余寿命概率密度函数
并与单退化量下主试齿轮箱1的剩余寿命概率密度函数进行比较，如图11和图12所示；
由图11和图12可知，基于Copula函数的多退化量下主试齿轮箱1的剩余寿命概率密度函数的方差比基于单退化量下主试齿轮箱1的剩余寿命概率密度函数的方差要小，其剩余寿命的预测值更接近剩余寿命的真实值，说明基于Copula函数的多退化量下主试齿轮箱1的剩余寿命预测值更准确；
步骤5、通过基于Copula函数的多退化量的寿命预测模型预测齿轮剩余寿命；最后得出齿轮的平均剩余寿命为：
图13和表1给出了基于核估计和随机滤波的单退化量下主试齿轮箱1的剩余寿命预测值与基于Copula函数的剩余寿命预测值以及剩余寿命真实值的比较；
在参数相同条件下，基于核估计和随机滤波的单退化量模型、基于Copula的多退化量模型的齿轮剩余寿命预测值与实际值对比误差如表3所示：
表3基于核估计和随机滤波的单退化量模型预测值、基于Copula的多退化量模型预测值与真实值误差比较
对比表3中的数据能够得出，随着系统运行时间的增长，状态监测信息的增多，剩余寿命预测值与实际值的绝对误差逐渐减小，同时，基于Copula函数的多退化量下主试齿轮箱1的剩余寿命概率密度函数的方差比基于核估计和随机滤波集成的单退化量下主试齿轮箱1的剩余寿命概率密度函数的方差要小，其剩余寿命的预测值更接近剩余寿命的真实值，说明基于Copula函数的多退化量下主试齿轮箱1的剩余寿命预测值更准确；说明本发明所提方法能够很好的进行实时剩余寿命预测；
综上所述，本发明提出基于多退化量监测的齿轮实时剩余寿命预测方法，该方法首先利用传感器得到不同退化量的数据，核密度估计方法利用这些数据对齿轮连续退化状态的概率密度函数进行非参数估计，得到齿轮基于不同退化量下的实时退化状态概率密度函数；利用实时状态监测数据来更新不同退化量下的随机滤波递推模型参数，预测齿轮在监测不同退化量下的剩余寿命，最后利用所选的最优的Copula函数，将单退化量的边缘分布函数与其联合分布函数联系起来，得到齿轮在多退化量下的剩余寿命分布。
Claims (1)
1.基于多退化量监测的齿轮实时剩余寿命预测方法，实施步骤如下：
步骤1、通过试验获取表征主试齿轮箱(1)内齿轮状态的实时监测数据：
采用试验台架的中心距为a＝150mm，试验采用机械杠杆(6)加载，扭矩采用转矩传感器(13#)进行测量，主陪试齿轮箱(1)(2)中均为正、反面交错搭接啮合的一对齿轮，齿轮的断齿状态等效为齿轮的失效，
试验共布置八个加速度传感器(1#～8#),两个噪声传感器(9#、10#)，一个温度传感器(11#),一个转速传感器(12#)，一个转矩传感器(13#)；有四个加速度传感器(1#～4#)分别布置在主试齿轮箱(1)轴承座的径向位置，有两个加速度传感器(7#、8#)分别布置在主试齿轮箱(1)的轴向位置，有两个加速度传感器(5#、6#)分别布置在陪试齿轮箱(2)轴承座的径向位置；有两个噪声传感器(9#、10#)分别悬挂在主试齿轮箱(1)和陪试齿轮箱(2)正上方40cm处；一个温度传感器(11#)布置在主试齿轮箱(1)内部；一个转速传感器(12#)布置在驱动电机(8)与陪试齿轮箱(2)联接轴的中部；一个转矩传感器(13#)布置在主试齿轮箱(1)与陪试齿轮箱(2)联接轴的中部；试验中依次加载八级载荷，八级载荷的大小在330～850扭矩之间，每级载荷的加载时间为9～12个小时，在第八级载荷发生断齿，通过用加速度传感器(4#)记录齿轮的加速度数据，剩余寿命预测选取从第八级加载开始到断齿的加速度测试点整个时域信号进行分析，采样信息如下：采样频率为20～30kHz，每次采样持续50～70秒，每隔8～10分钟记录一次采样文本；
步骤2、对主试齿轮箱(1)内齿轮的退化状态进行特征提取，利用均方幅值对齿轮磨损退化性能进行衰退评估，对于每次采样时间长度内，采样信号的均方幅值特征值表示为：
式中：∑为求和号，n为每个采样周期的采样点数,yi为ti时刻齿轮的状态信息，yj为齿轮每个采样周期的数据；通过均方幅值对加速度和噪声数据处理；
步骤3、根据主试齿轮箱(1)内齿轮的初始故障数据，利用核密度估计方法以及随机滤波理论，建立单退化量的退化轨迹，得到主试齿轮箱(1)内齿轮的剩余寿命概率密度函数；
对于主试齿轮箱(1)的剩余寿命概率密度函数fi,k(xi,k|Yi,k)，通过贝叶斯定理获得fi,k(xi,k|Yi,k)的递推形式；
由贝叶斯定理得：
式中：xi,k表示在单退化量i下tk时刻主试齿轮箱(1)的剩余寿命；Yi,k表示第i个退化量tk时刻的历史监测数据，i＝1表示单退化量为加速度，i＝2表示单退化量为噪声；
因为：
fi,k(yi,k|xi,k,Yi,k-1)＝fi,k(yi,k|xi,k) (3)
所以：
在tk时刻的剩余寿命等于tk-1时刻的剩余寿命减去时刻tk和时刻tk-1之间的间隔，即：
根据式(5)和式(6)，得到随机滤波方程为：
在t1时刻，可知fi,0(xi,0-t1+t0|Yi,0)＝fX,0(xi,0-t1+t0)，可求得：
在t2时刻，可求得：
估计出fi,0(xi,0)和fi,k(yi,k|xi,k)，通过递推公式(7)计算得到fi,k(xi,k|Yi,k)；
用核密度估计fi,0(xi,0)和fi,k(yi,k|xi,k)，即：
设X1，X2，…，Xn为独立同分布的表征退化的随机变量，其概率密度函数为f(x)，则：
式中：h为平滑参数,K(u)为核函数，n为随机变量X的样本数；
由式(10)得：
式中：
(
为样本方差)为平滑参数，
为核函数，ni为随机变量Xi的样本数；
式中：αi、βi为平滑参数,K1(a,b)、K2(u)为核函数；
根据核密度估计结果以及式(7)，得到齿轮箱在时间tk处基于加速度的剩余寿命概率密度函数为f1,k(x1,k|Y1,k)和基于噪声的剩余寿命概率密度函数为f2,k(x2,k|Y2,k)；
根据概率论中概率密度函数与边缘分布函数的关系，基于加速度的剩余寿命的边缘分布函数为
和基于噪声的剩余寿命的边缘分布函数为
步骤4、齿轮的单退化量之间的关系是不知道的，Copula函数是一种相关性分析方法，将单退化量的边缘分布函数与其联合分布函数联系起来，即：
F(t1,t2)＝C(F1(t1),F2(t2)；γ)＝C(u1,u2；γ) (13)
式中：F(t1,t2)为变量t1和t2的联合分布函数；u1＝F1(t1)和u2＝F2(t2)分别为变量t1和t2的边缘分布函数；γ为Copula函数的参数；
变量t1和t2的边缘分布函数F1(t1)和F2(t2)已知，Copula函数的未知参数γ采用极大似然估计法估计：
式中：c为C的密度函数，∑为求和号，n为每个采样周期的采样点数；
在工程实践中，Copula常用的3种函数形式，如表1所示：
表1常用的Copula函数形式
在同一工程实践中选择不同形式的Copula函数，出现不同的结果；经过初步筛选，选择符合尾部相关性的Copula函数；常用的3种Copula函数中不同参数γ下尾部的非对称性不明显，通过模型评价的方法选择最优的Copula函数，AIC准则是最常用的模型评价准则，具有良好的适用性和实用性，即：
式中：L为Copula函数的似然值；n为样本个数；K为被估计参数个数；
常用Copula函数的似然函数如下：
Gumbel函数的似然函数：
Clayton函数的似然函数：
Frank函数的似然函数：
AIC表示Copula函数对数据的拟合效果，AIC值越小，对数据的拟合效果越好，在求得主试齿轮箱(1)的加速度和噪声作为退化量下剩余寿命的边缘分布后，分别通过GumbelCopula、Clayton Copula、Frank Copula这三种常用的Copula函数对基于振动加速度和噪声的剩余寿命边缘分布函数进行融合，得到齿轮箱的剩余寿命联合分布，根据AIC准则选择最优的Copula函数；
利用所选的最优的Copula函数的性质，融合基于加速度和噪声的边缘分布函数，得到主试齿轮箱(1)的剩余寿命联合分布函数F(t1,t2)＝C(F1(t1),F2(t2)；γ)，根据概率论中概率密度函数与边缘分布函数的关系，得到主试齿轮箱(1)基于多退化量下的剩余寿命概率密度函数
步骤5、通过基于Copula函数的多退化量的寿命预测模型预测齿轮剩余寿命；最后得出齿轮的平均剩余寿命为：
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910432767.4A CN110174261B (zh) | 2019-05-23 | 2019-05-23 | 多退化量监测的齿轮实时剩余寿命预测方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910432767.4A CN110174261B (zh) | 2019-05-23 | 2019-05-23 | 多退化量监测的齿轮实时剩余寿命预测方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110174261A CN110174261A (zh) | 2019-08-27 |
| CN110174261B true CN110174261B (zh) | 2020-09-08 |
Family
ID=67691932
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910432767.4A Active CN110174261B (zh) | 2019-05-23 | 2019-05-23 | 多退化量监测的齿轮实时剩余寿命预测方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110174261B (zh) |
Families Citing this family (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110990788A (zh) * | 2019-12-02 | 2020-04-10 | 宁海县浙工大科学技术研究院 | 一种基于三元维纳过程的轴承剩余寿命预测方法 |
| CN111896255B (zh) | 2020-08-12 | 2021-06-04 | 上海理工大学 | 基于多轴随机道路载荷下轮毂轴承服役寿命快速评估方法 |
| CN113468801B (zh) * | 2021-06-07 | 2024-03-26 | 太原科技大学 | 一种齿轮核密度估计剩余寿命预测方法 |
| CN113468721B (zh) * | 2021-06-07 | 2024-03-29 | 太原科技大学 | 一种对齿轮减速箱中齿轮和轴承剩余寿命的预测方法 |
| CN115758094B (zh) * | 2023-01-09 | 2023-05-09 | 河北工业大学 | 一种减速器退化指标提取及寿命预测方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR20010000116A (ko) * | 2000-05-06 | 2001-01-05 | 손원열 | 가스 분석, 농도추정 및 측정, 측정 데이터 보정방법과그의 표시방법 |
| CN105740625A (zh) * | 2016-01-31 | 2016-07-06 | 太原科技大学 | 一种齿轮的实时剩余寿命预测方法 |
| CN108645615A (zh) * | 2018-04-08 | 2018-10-12 | 太原科技大学 | 一种自适应模糊神经网络齿轮剩余寿命预测方法 |
| CN109085814A (zh) * | 2018-07-23 | 2018-12-25 | 西安热工研究院有限公司 | 一种火电汽轮机组整体设备系统延寿评估方法 |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JPH07181784A (ja) * | 1993-12-22 | 1995-07-21 | Minolta Co Ltd | 露光量制御方法およびその装置 |
| CN102279928B (zh) * | 2011-07-20 | 2013-04-03 | 北京航空航天大学 | 基于支持向量机和模糊信息粒化的产品性能退化区间预测方法 |
| CN102542155B (zh) * | 2011-12-05 | 2014-09-17 | 北京航空航天大学 | 基于加速退化数据的粒子滤波剩余寿命预测方法 |
| CN103246803B (zh) * | 2013-04-07 | 2016-05-04 | 河南科技大学 | 一种滚动轴承性能变异过程的显著性检验方法 |
-
2019
- 2019-05-23 CN CN201910432767.4A patent/CN110174261B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR20010000116A (ko) * | 2000-05-06 | 2001-01-05 | 손원열 | 가스 분석, 농도추정 및 측정, 측정 데이터 보정방법과그의 표시방법 |
| CN105740625A (zh) * | 2016-01-31 | 2016-07-06 | 太原科技大学 | 一种齿轮的实时剩余寿命预测方法 |
| CN108645615A (zh) * | 2018-04-08 | 2018-10-12 | 太原科技大学 | 一种自适应模糊神经网络齿轮剩余寿命预测方法 |
| CN109085814A (zh) * | 2018-07-23 | 2018-12-25 | 西安热工研究院有限公司 | 一种火电汽轮机组整体设备系统延寿评估方法 |
Non-Patent Citations (2)
| Title |
|---|
| Remaining Useful Life Prediction of Gearbox;Guoyu Lin等;《2013 International Conference on Quality, Reliability, Risk, Maintenance, and Safety Engineering (QR2MSE)》;20130718;全文 * |
| 基于核密度估计和随机滤波理论的齿轮箱;石慧等;《计算机集成制造系统》;20190121;全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110174261A (zh) | 2019-08-27 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110174261B (zh) | 多退化量监测的齿轮实时剩余寿命预测方法 | |
| CN109883691B (zh) | 核估计和随机滤波集成的齿轮剩余寿命预测方法 | |
| CN109376401B (zh) | 一种自适应多源信息优选与融合的机械剩余寿命预测方法 | |
| Shi et al. | Rolling bearing initial fault detection using long short-term memory recurrent network | |
| Ghasemi et al. | Evaluating the reliability function and the mean residual life for equipment with unobservable states | |
| Saidi et al. | An integrated wind turbine failures prognostic approach implementing Kalman smoother with confidence bounds | |
| CN110414154B (zh) | 一种带有双测点的风机部件温度异常检测和报警方法 | |
| CN108304348B (zh) | 一种基于二元维纳过程的轴承剩余寿命预测方法 | |
| CN107797537A (zh) | 一种应用于自动化生产线的故障预测与健康管理方法 | |
| KR102119661B1 (ko) | 회전기기 고장 예지를 위한 건전성 지표 추이 및 잔존수명 예측 기법 | |
| Zhao et al. | An integrated prognostics method for failure time prediction of gears subject to the surface wear failure mode | |
| Wang | An intelligent system for machinery condition monitoring | |
| WO2015052408A1 (fr) | Surveillance d'un moteur d'aéronef pour anticiper les opérations de maintenance | |
| CN113468721B (zh) | 一种对齿轮减速箱中齿轮和轴承剩余寿命的预测方法 | |
| Ghasemi et al. | Estimating mean residual life for a case study of rail wagon bearings | |
| Salunkhe et al. | Prediction of Remaining Useful Life of mechanical components-a Review | |
| Anis | Towards remaining useful life prediction in rotating machine fault prognosis: an exponential degradation model | |
| Orth et al. | Accuracy and robustness of decision making techniques in condition based maintenance | |
| Lybeck et al. | Validating prognostic algorithms: a case study using comprehensive bearing fault data | |
| CN109580218B (zh) | 一种基于似然学习机的风机齿轮箱状态识别方法 | |
| CN110889186B (zh) | 灰关联分析的加速贮存与自然贮存退化数据一致性检验法 | |
| CN113468801B (zh) | 一种齿轮核密度估计剩余寿命预测方法 | |
| CN110895624A (zh) | 基于最大熵谱估计的加速贮存与自然贮存退化数据一致性检验法 | |
| CN116070368A (zh) | 一种海上风电机组高速轴承剩余寿命预测方法 | |
| CN110889077A (zh) | 基于Kendall相关系数的加速贮存与自然贮存退化数据一致性检验方法 |
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 | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| TR01 | Transfer of patent right | ||
| TR01 | Transfer of patent right |
Effective date of registration: 20231018 Address after: No. 9 Xinhe Road, Yangqu Industrial Park, Shanxi Transformation Comprehensive Reform Demonstration Zone, Taiyuan City, Shanxi Province, 030100 Pb108, South Building, Xinhe Science and Technology Park Patentee after: Shanxi Zhida Intelligent Equipment Co.,Ltd. Address before: 030024 No. 66 tile Road, Wan Berlin District, Shanxi, Taiyuan Patentee before: TAIYUAN University OF SCIENCE AND TECHNOLOGY |