CN105956236B - 一种双重更新的四因素随机退化模型齿轮寿命预测方法 - Google Patents

Abstract

一种双重更新的四因素随机退化模型齿轮寿命预测方法，充分考虑实际齿轮退化过程中广泛存在的时变性、个体差异性、非线性和测量不确定性四方面随机因素，建立贴合其实际退化过程的四因素随机退化模型，并推导出该模型下齿轮剩余寿命基于观测值的概率密度函数；采用MLE方法对该模型参数进行初始化，避免参数初始化的盲目性；引入双重更新方法对模型状态及参数进行更新，获得更准确的更新结果，进而很大程度上提高了齿轮剩余寿命的预测精度。

Description

一种双重更新的四因素随机退化模型齿轮寿命预测方法

技术领域

本发明属于机械设备剩余寿命预测技术领域，具体涉及一种双重更新的四因素随机退化模型齿轮寿命预测方法。

背景技术

剩余寿命预测在齿轮状态监测数据的基础上为齿轮维护决策提供直接的指导和依据，成为齿轮预知性维护策略不可或缺的重要环节。基于随机退化模型的齿轮寿命预测方法将齿轮的剩余寿命视作随机变量，从监测信号中提取齿轮的健康状态信息，继而从概率密度角度建立其剩余寿命的预测分布。该方法能够充分挖掘齿轮状态监测数据的动态特性和随机特性，契合齿轮的失效机理，对齿轮的退化过程进行准确描述，因而受到众多学者的广泛认同。在实际应用中，齿轮退化的随机性往往表现在四个方面：时变性、个体差异性、非线性和测量不确定性，应用随机退化模型对齿轮的退化过程进行建模描述时需要同时考虑这四种因素，然而现有模型大都仅考虑部分随机性因素，加之模型参数未进行初始化，模型参数和系统状态更新方法表现欠佳，致使获得的齿轮剩余寿命概率分布偏离实际，剩余寿命预测精度不尽人意。

发明内容

为了克服上述现有技术的上述缺点，本发明提供一种双重更新的四因素随机退化模型齿轮寿命预测方法，使得齿轮剩余寿命预测结果精度更高。

为了达到上述目的，本发明采取的技术方案为：

一种双重更新的四因素随机退化模型齿轮寿命预测方法，包括以下步骤：

1)使用监测系统采集齿轮振动信号，提取健康指标形成观测值向量y＝[y₁,...,y_s]^T，s为独立观测值个数；

2)构建四因素随机退化模型：

其中，

x_k—t_k时刻齿轮的真实退化状态；

θ_k—t_k时刻衰退系数，表征个体差异性，为随机变量，服从分布

b—模型阶数，表征非线性；

ω_k-1—过程噪声，服从分布N(0,σ²(t_k-t_k-1))，ω_k-1＝σ(B(t_k)-B(t_k-1))，B(t_k)为标准布朗运动，σB(t_k)表征时变性；

y_k—t_k时刻齿轮健康指标观测值；

ν—观测噪声，服从分布N(0,γ²)；

3)采用极大似然估计(MLE)方法对式(1)四因素随机退化模型的参数进行初始化，具体步骤为：

3.1)t_k时刻齿轮健康指标观测值y_k～N(x_k,γ²)，根据步骤1)所得观测值向量y计算获得对数似然函数：

其中，为四因素随机退化模型的未知参数向量，Ω＝σ²M+γ²Ι_s，M＝[min{t_i,t_j}]_1≤i,j≤s，Ι_s为s阶单位矩阵；

3.2)计算式(2)关于μ_θ和的偏导，极大化似然函数：

其中，求解方程组(3)得到：

3.3)将式(4)计算结果回代至式(2)，简化对数似然函数为：

3.4)通过Matlab提供的寻优命令对式(5)进行寻优获得t和回带至式(4)求解出μ_θ和进而根据关系Ω＝σ²M+γ²Ι_s，求解出未知参数向量的极大似然估计，作为式(1)四因素随机退化模型的初始参数

4)将式(1)四因素随机退化模型表示为式(6)形式，并带入初始参数Θ₀：

其中，z_k＝[x_k,θ_k]^T，

5)引入双重更新方法，对四因素随机退化模型的状态和参数进行实时更新，该方法由两次相互补充、嵌套进行的更新过程组成：

5.1)外层更新开始，令k＝0，对正态分布即进行随机采样，获得由N个初始粒子构成粒子集对应的粒子权值集合为其中，

5.2)里层更新开始，令k＝k+1，将P_k-1|k-1和N个未更新粒子带入式(7)迭代更新，计算获得P_k|k和

其中，

5.3)里层更新结束，由步骤5.2)所得P_k|k和形成正态分布列对每个分布进行一次采样获得新的粒子集各粒子权值更新如下：

其中，

5.4)由步骤5.3)所得粒子及权值计算t_k时刻齿轮更新状态和参数：

外层更新结束；

6)将步骤5.4)的更新结果代入式(10)，获得t_k时刻基于观测值的齿轮剩余寿命概率密度f_Lk|y1:k(l_k|y_1:k)，式(10)为t_k时刻基于观测值的齿轮剩余寿命概率密度理论推导结果，

其中，l_k为t_k时刻齿轮剩余寿命，λ为失效阈值，ψ(l_k)＝bl_k(l_k+t_k)^b-1，

7)将步骤6)所得f_Lk|y1:k(l_k|y_1:k)代入式(11)，计算t_k时刻齿轮预测剩余寿命为：

8)判断齿轮更新状态是否达到失效阈值λ，如未达到则获取其下一时刻的健康指标观测值，令k＝k+1，重复步骤5.2)至步骤7)，直到达到失效阈值停止。

本发明的有益效果：充分考虑实际齿轮退化过程中广泛存在的时变性、个体差异性、非线性和测量不确定性四方面随机因素，建立贴合其实际退化过程的四因素随机退化模型，并推导出该模型下齿轮剩余寿命基于观测值的概率密度函数；采用MLE方法对该模型参数进行初始化，避免参数初始化的盲目性；引入双重更新方法对模型状态及参数进行更新，获得更准确的更新结果，进而很大程度上提高了齿轮剩余寿命的预测精度。

附图说明

图1为本发明流程图。

图2为齿轮加速退化寿命实验台结构图。

图3为测试齿轮全寿振动信号：图3(a)为第一组测试齿轮全寿振动信号；图3(b)为第二组测试齿轮全寿振动信号。

图4为测试齿轮RMS指标和起始预测时刻：图4(a)为第一组测试齿轮的RMS指标和起始预测时刻；图4(b)为第二组测试齿轮的RMS指标和起始预测时刻。

图5为五种方法对测试齿轮健康状态的估计结果：图5(a)为五种方法对第一组测试齿轮健康状态估计结果；图5(b)为五种方法对第二组测试齿轮健康状态估计结果。

图6为五种方法对测试齿轮剩余寿命的预测结果：图6(a)为五种方法对第一组测试齿轮剩余寿命预测结果；图6(b)为五种方法对第二组测试齿轮剩余寿命预测结果。

具体实施方式

下面结合附图和实施例对本发明做进一步详细描述。

如图1所示，一种双重更新的四因素随机退化模型齿轮寿命预测方法，包括以下步骤：

1)使用监测系统采集齿轮振动信号，提取健康指标形成观测值向量y＝[y₁,...,y_s]^T，s为独立观测值个数；

2)构建四因素随机退化模型：

其中，

x_k—t_k时刻齿轮的真实退化状态；

θ_k—t_k时刻衰退系数，表征个体差异性，为随机变量，服从分布

b—模型阶数，表征非线性；

ω_k-1—过程噪声，服从分布N(0,σ²(t_k-t_k-1))，ω_k-1＝σ(B(t_k)-B(t_k-1))，B(t_k)为标准布朗运动，σB(t_k)表征时变性；

y_k—t_k时刻齿轮健康指标观测值；

ν—观测噪声，服从分布N(0,γ²)；

3)采用极大似然估计(MLE)方法对式(1)四因素随机退化模型的参数进行初始化，具体步骤包括：

3.1)t_k时刻齿轮健康指标观测值y_k～N(x_k,γ²)，根据步骤1)所得观测值向量y计算获得对数似然函数：

其中，为四因素随机退化模型的未知参数向量，Ω＝σ²M+γ²Ι_s，M＝[min{t_i,t_j}]_1≤i,j≤s，Ι_s为s阶单位矩阵；

3.2)计算式(2)关于μ_θ和的偏导，极大化似然函数：

其中，求解方程组(3)得到：

3.3)将式(4)计算结果回代至式(2)，简化对数似然函数为：

3.4)通过Matlab提供的寻优命令对式(5)进行寻优获得t和回带至式(4)求解出μ_θ和进而根据关系Ω＝σ²M+γ²Ι_s，求解出未知参数向量的极大似然估计，作为式(1)四因素随机退化模型的初始参数

4)将式(1)四因素随机退化模型表示为式(6)形式，并带入初始参数Θ₀：

其中，z_k＝[x_k,θ_k]^T，

5)引入双重更新方法，对四因素随机退化模型的状态和参数进行实时更新，该方法由两次相互补充、嵌套进行的更新过程组成：

5.1)外层更新开始，令k＝0，对正态分布即进行随机采样，获得由N个初始粒子构成粒子集对应的粒子权值集合为其中，

5.2)里层更新开始，令k＝k+1，将P_k-1|k-1和N个未更新粒子带入式(7)迭代更新，计算获得P_k|k和

其中，

5.3)里层更新结束，由步骤5.2)所得P_k|k和形成正态分布列对每个分布进行一次采样获得新的粒子集各粒子权值更新如下：

其中，

5.4)由步骤5.3)所得粒子及权值计算t_k时刻齿轮更新状态和参数：

外层更新结束；

6)将步骤5.4)的更新结果代入式(10)，获得t_k时刻基于观测值的齿轮剩余寿命概率密度f_Lk|y1:k(l_k|y_1:k)，式(10)为t_k时刻基于观测值的齿轮剩余寿命概率密度理论推导结果，

其中，l_k为t_k时刻齿轮剩余寿命，λ为失效阈值，ψ(l_k)＝bl_k(l_k+t_k)^b-1，

7)将步骤6)所得f_Lk|y1:k(l_k|y_1:k)代入式(11)，计算t_k时刻齿轮预测剩余寿命为：

8)判断齿轮更新状态是否达到失效阈值λ，如未达到则获取其下一时刻的健康指标观测值，令k＝k+1，重复步骤5.2)至步骤7)，直到达到失效阈值停止。

实施例：从齿轮加速退化寿命实验台上采集齿轮加速退化实验数据对本发明进行验证，实验台结构图如图2所示，该实验台采用与FZG标准减速机相似的背靠背式结构，主体包括一个驱动电机和两个齿轮箱，齿轮箱分别为陪侍齿轮箱和测试齿轮箱，前者安装斜齿轮，后者安装直齿圆柱齿轮。测试齿轮安装位置为测试齿轮箱齿轮啮合对中小齿轮，小齿轮轴上安装有联轴器实现扭矩加载。实验过程中，驱动电机转速为1460rpm，负载扭矩约为1300Nm。加速度传感器安装在测试齿轮箱箱体上，采用间歇采样的形式，采样频率为20kHz，采样持续时间为6s，即数据长度为120000。实验共进行两组，两组测试齿轮上均未提前加工故障，实验结束后在第一组测试齿轮上观察到磨损、点蚀和裂纹故障，在第二组测试齿轮上观察到严重裂纹故障。加速度传感器采集得到的两组测试齿轮的全寿命振动信号如图3所示。

分别提取两组振动测试信号的RMS值作为健康指标，并选择起始预测时刻如图4所示，分别为2088s和2098s。为获得当前时刻的剩余寿命预测结果，将自起始预测时刻至当前时刻之间的RMS值作为观测值，输入本发明四因素随机退化模型，使用MLE方法对四因素随机退化模型参数进行初始化，循环进行状态和参数双重更新、计算剩余寿命近似概率分布、计算当前时刻齿轮预测寿命、提取下一个健康指标，直到齿轮的状态估计值达到预先设定的失效阈值。

本发明方法概括为“非线性退化模型+KPF更新”方法(M4)，为验证其有效性，引入其他四种寿命预测方法与之进行比较，分别为：“线性退化模型+KF更新”方法(M1)、“非线性退化模型+KF更新”方法(M2)、“非线性退化模型+PF更新”方法(M3)和“Paris-Erdogan模型+KPF更新”方法(M5)。分别应用五种方法对两组齿轮进行剩余寿命预测，状态估计结果和寿命预测结果如图5和图6所示。从图5中可以看出，M3及本发明方法的状态估计结果曲线更加平滑，说明两方法较其他三种方法具有更加抗噪声干扰能力；从图6中可以看出，随着预测时刻接近齿轮寿命终点，本发明方法的剩余寿命预测结果能够更快速的收敛于真实寿命。

为定量验证本发明方法较于以上四种方法有更优良的剩余寿命预测表现，分别计算五种方法的累计相对精度(CRA)，结果如表1所示，可以看出，本发明方法在两组齿轮实例中均取得了最大值，验证了其更高的寿命预测精确度。

表1

本发明提出一种双重更新的四因素随机退化模型齿轮寿命预测方法，适用对象为齿轮，但齿轮与其他常见旋转零部件退化机理相似，因此本发明方法作出调整后同样适用于与本发明模型描述随机退化特征相似的其他零部件或机械设备。此外，由于加速退化实验通过增大负载加速齿轮退化并未改变其退化机理，本发明方法同样适用于实际应用中缓慢随机退化齿轮的剩余寿命预测，通过更换模型输入指标、修改模型参数或拓展模型适用条件等都未脱离本发明的构思，也应视为本发明的保护范围。

Claims (1)

1.一种双重更新的四因素随机退化模型齿轮寿命预测方法，其特征在于，包括以下步骤：

1)使用监测系统采集齿轮振动信号，提取健康指标形成观测值向量y＝[y₁,...,y_s]^T，s为独立观测值个数；

2)构建四因素随机退化模型：

其中，

x_k—t_k时刻齿轮的真实退化状态；

θ_k—t_k时刻衰退系数，表征个体差异性，为随机变量，服从分布N(μ_θ,σ_θ ²)；

b—模型阶数，表征非线性；

ω_k-1—过程噪声，服从分布N(0,σ²(t_k-t_k-1))，ω_k-1＝σ(B(t_k)-B(t_k-1))，B(t_k)为标准布朗运动，σB(t_k)表征时变性；

y_k—t_k时刻齿轮健康指标观测值；

ν—观测噪声，服从分布N(0,γ²)；

3)采用极大似然估计(MLE)方法对式(1)四因素随机退化模型的参数进行初始化，具体步骤为：

3.1)t_k时刻齿轮健康指标观测值y_k～N(x_k,γ²)，根据步骤1)所得观测值向量y计算获得对数似然函数：

其中，为四因素随机退化模型的未知参数向量，Ω＝σ²M+γ²Ι_s，M＝[min{t_i,t_j}]_1≤i,j≤s，Ι_s为s阶单位矩阵；

3.2)计算式(2)关于μ_θ和的偏导，极大化似然函数：

其中，求解方程组(3)得到：

3.3)将式(4)计算结果回代至式(2)，简化对数似然函数为：

3.4)通过Matlab提供的寻优命令对式(5)进行寻优获得t和回带至式(4)求解出μ_θ和进而根据关系Ω＝σ²M+γ²Ι_s，求解出未知参数向量的极大似然估计，作为式(1)四因素随机退化模型的初始参数

4)将式(1)四因素随机退化模型表示为式(6)形式，并带入初始参数Θ₀：

其中，z_k＝[x_k,θ_k]^T，

5)引入双重更新方法，对四因素随机退化模型的状态和参数进行实时更新，该方法由两次相互补充、嵌套进行的更新过程组成：

5.1)外层更新开始，令k＝0，对正态分布即进行随机采样，获得由N个初始粒子构成粒子集对应的粒子权值集合为其中，

5.2)里层更新开始，令k＝k+1，将P_k-1|k-1和N个未更新粒子带入式(7)迭代更新，计算获得P_k|k和

其中，

5.3)里层更新结束，由步骤5.2)所得P_k|k和形成正态分布列对每个分布进行一次采样获得新的粒子集各粒子权值更新如下：

其中，

5.4)由步骤5.3)所得粒子及权值计算t_k时刻齿轮更新状态和参数：

外层更新结束；

6)将步骤5.4)的更新结果代入式(10)，获得t_k时刻基于观测值的齿轮剩余寿命概率密度式(10)为t_k时刻基于观测值的齿轮剩余寿命概率密度理论推导结果，

其中，l_k为t_k时刻齿轮剩余寿命，λ为失效阈值，

7)将步骤6)所得代入式(11)，计算t_k时刻齿轮预测剩余寿命为：

8)判断齿轮更新状态是否达到失效阈值λ，如未达到则获取其下一时刻的健康指标观测值，令k＝k+1，重复步骤5.2)至步骤7)，直到达到失效阈值停止。