[发明专利]一种齿轮的实时剩余寿命预测方法

权利要求书

1.一种齿轮的实时剩余寿命预测方法，其特征在于实施步骤是：

步骤1、在齿轮箱安装传感器，获取表征齿轮退化的实时监测数据；

加速度传感器安装在主试箱的轴承座位置，在齿轮箱内安装温度传感器，在主试箱的正上方安装噪声传感器，

步骤2、对齿轮的疲劳状态进行特征提取，对齿轮磨损退化性能进行衰退评估；

利用均方幅值对齿轮磨损退化性能进行衰退评估，对于每次采样时间Δt长度内，离散随机信号的时间序列均方幅值可表示为：

xrms(Δt)=(1nΣi=1nxi2(t))1/2---(1)]]>

式中Δt为采用时间；n＝F_s×Δt，F_s为采样频率，n为采样点数，∑表示求和，i∈(1,2,3....n)，x_i(t)为采样值；

步骤3、建立可变参数的齿轮退化状态空间模型：

通过非线性状态空间建模进行状态估计，模型的一般形式为：

X_t＝f(X_t-1)+w_t-1(2)

Y_t＝h(X_t)+v_t(3)

式中，X_t为t时刻的系统状态向量，t≥1，f(X_t-1)为系统状态转移模型；Y_t为t时刻的系统观测向量，h(X_t)为系统观测模型，系统噪声w_t和观测噪声v_t服从零均值的高斯分布，其协方差阵为：Var(w_t)＝Q_t，Var(v_t)＝R_t；

采用线性化进行简化，线性化后的状态空间模型为：

X_t＝F_t|t-1X_t-1+w_t-1(4)

Y_t＝H_tX_t+v_t(5)

(4)式和(5)式中的未知参数用一个参数向量θ表示，即θ＝{F_t|t-1,H_t,Q_t,R_t}，其中，F_t|t-1为系统状态转移矩阵，H_t为系统观测矩阵，对系统的状态向量进行i次观测，从而得到观测序列Y_1:i＝{Y₁,Y₂,…,Y_i}，当系统在运行过程中，随着外界环境的变化，模型参数发生变化，因此可以在可实时监测条件下，利用监测数据进行状态空间模型参数的估计及退化状态的预测，求得在t时刻系统状态向量X_t的最优估计量

首先采用期望最大化算法从初始时刻以迭代的方式实时估计状态空间模型参数，该算法包括E步骤和M步骤两个步骤，具体步骤如下：

E步骤是估计未知参数向量θ过程中，将卡尔曼滤波中的待估计状态初值X₀作为未知的“缺失数据”，其均值和方差分别为a₀,未知，滤波估计误差初值为P₀，待估计参数的初值为：

θ₀＝{F₀,H₀,Q₀,R₀}(6)

为使退化状态初值选择更为合理，采用卡尔曼前向滤波及固定区间反向滤波相结合的方法在参数估计的同时来选取最优的退化状态初值，

卡尔曼前向滤波时，在t₁,t₂,...,t_i,t_i+1,...,t_i+n时刻，其相应的退化状态估计值为：一步预测值为：

X^t|t-1=Ft|t-1X^t-1---(7)]]>

校正预测估计得到的状态估计值为：

X^t=X^t|t-1+Kt[Yt-HtX^t|t-1]---(8)]]>

Kt=Pt|t-1HtT[HtPt|t-1HtT+Rt]-1---(9)]]>

式中K_t为滤波增益，P_t|t-1为预测的误差方差阵，P_t为滤波估计误差的方差阵，I是单位矩阵：

P_t|t-1＝F_t|t-1P_t-1F_t|t-1^T+Q_t-1 (10)

Pt=[I-KtHt]Pt|t-1[I-KtHt]T+KtRtKtT---(11)]]>

已知在t_i+n时刻，接收到i+n个监测点的监测信息，求得其对数最大似然函数为：

ln L(θ)=ln[Πt=t0ti+np(Yt|Xt,θ)]=Σt=t1ti+nln p(Yt|θ)-ln p(X0|θ)Σt=t1ti+nln p(Xt|Xt-1,θ)---(12)]]>

式中Π表示求积，a₀为X₀的均值，为X₀的方差，可求得似然函数的表达式为：

ln L(θ)=-{ln|σ02|+(x0-a0)2σ02+(i+n)ln|Qti+n|+Σt=t1ti+n(xt-Fti+n|ti+n-1xt-1)2Qti+n+(i+n)ln|Rti+n|+Σt=t1ti+n(Yt-Hti+nXt)2Rti+n}---(13)]]>

采用期望最大化算法估计参数，计算式(13)的期望值，记为

EX|Y,θk[ln L(θ)]=E{-[ln|σ02|+(x0-a0)2σ02+(i+n)ln|Qti+n|+Σt=t1ti+n(xt-Fti+n|ti+n-1xt-1)2Qti+n+(i+n)ln|Rti+n|+Σt=t1ti+n(Yt-Hti+nXt)2Rti+n]}---(14)]]>

卡尔曼光滑算法即为求期望值的过程，首先求得：

EX|Y,θk(Xt)=X^t|T---(15)]]>

St|T=EX|Y,θk(Xt2)=X^t|T2+Pt|T---(16)]]>

St,t-1|T=EX|Y,θk(XtXt-1)=X^t|TX^t-1|T+Jt-1Pt|T---(17)]]>

其中T＝t_i+n,t_i+n-1,...,t_i,t_i-1,...,t₂,t₁，Y_T表示所有监测信息，与P_t|T可通过固定区间卡尔曼反向滤波求得：

Jt=PtFt|t-1T(Pt|t-1)-1---(18)]]>

X^t|T=X^t+Jt[X^t+1|T-X^t+1|t]---(19)]]>

P_t|T＝P_t-J_t(P_t+1|T-P_t+1|t)J_t^T(20)

将(15)式-(20)式代入(14)式，可求得：

EX|Y,θk[ln L(θ)]=-{ln|σ02|+St0|T-2a0X^t|T+a02σ02+(i+n)ln|Qti+n|+Σt=t1ti+nSt|T+Fti+n|ti+n-12St-1|T-2Fti+n|ti+n-1St,t-1|TQti+n+(i+n)ln|Rti+n|+Σt=t1ti+nYt2-2Hti+nYtX^t|T+Hti+n2St|TRti+n}---(21)]]>

M步骤是通过求解最大似然估计算法求得各参数的估计值：

θ^k+1=argmaxθ{EX|Y,θk[ln L(θ)]}---(22)]]>

期望最大化算法不断进行E步骤和M步骤的迭代，直到似然函数值稳定不再变化，即求解到t_i+n时刻的参数估计值

步骤4、对齿轮退化过程中的突变状态点进行检测：

齿轮退化过程中，由于其退化磨损过程出现退化突变点，利用累积和控制图法进行突变点检测，对多次观测的结果累计求和来判断监测过程是否处于统计控制状态，检测出突变状态点，则递推公式为：

Sn+=maxΣi=1n[yi-μ0-k]=Sn-1++max(yn-μ0-k)---(23)]]>

其中表示检测均值向上偏移的累积和，y_i为i时刻的观测值，y_n为n时刻的观测值，h为突变状态检测阈值，Δμ为最小预警退化量，样本均值为μ₀，方差为当时，说明系统运行正常，没有突变点；当时就意味着系统存在异常突变点，

突变检测阈值可以采用平均运行长度ARL近似计算公式来确定：

ARL=exp(-c)+c-12(Δμ)2---(24)]]>

式中，参数c为：

c＝-(Δμ)(|h|/σ₁+1.166)＝-2k(|h|/σ₁+1.166)(25)

步骤5、根据检测到的突变点信息修正预测模型，重新进行状态估计；

当在t_i+n时刻检测到t_i时刻为退化状态突变点，对齿轮退化状态预测模型进行以下修正：

突变点后状态空间模型的修正首先是对突变点处退化状态初始值进行修正，突变点即t_i时刻滤波估计误差初值为：利用从初始时刻到t_i+n时刻所接收到的所有监测值采用(6)式-(22)式期望最大化算法参数估计算法求得参数估计值记为：

θti+n={Fti+n|ti+n-1,Hti+n,Qti+n,Rti+n}---(26)]]>

当在t_i+n时刻检测到t_i时刻存在突变点A，需修正突变点A之后的卡尔曼滤波模型初始值及参数，

突变点t_i之后，时间t重新记为：其相应的退化状态估计值为利用突变点后监测到的退化状态值重新估计参数记为θⁱ：

利用t_i时刻及之前接收到的所有监测信息采用期望最大化算法参数估计算法求得的参数估计值记为：

θti={Fti|ti-1,Hti,Qti,Rti}---(27)]]>

t_i时刻修正后待估计参数初值记为因此卡尔曼滤波待估计参数的初值修正记为：

θ0i={F0i,H0i,Q0i,R0i}=βθti+n+γθti---(28)]]>

其中β与γ为参数修正因子，且β+γ＝1，选择γ＜β，通过参数修正因子来加大突变点之后的数据的权重，相应减少突变点之前的数据的权重，

检测到突变点后，通过自适应增益函数增加突变点后的增益矩阵K_t的模值，使得状态估计比较快的收敛于部件退化的真实状态，假设检测到t_i处退化状态发生突变，构造一自适应增益函数g_i(t_i)为：

gi(ti)=1Kti{1+exp[-ati/b]}---(29)]]>

其中a与b为自适应增益函数的参数，增大滤波增益的方法是以牺牲滤波的最优性为代价的，因此自适应增益函数的参数a与b的选择应同时考虑其收敛速度与滤波最优性的问题，

则修正后的卡尔曼增益记作：

Kti′=Ktigi(ti)---(30)]]>

t_i时刻的待估计参数初值修正后，利用t_i时刻后的n个监测信息，记作重复上述卡尔曼滤波及参数估计过程，通过自适应滤波的方法，不断对t_i时刻之后未知的系统模型参数和噪声统计特性θⁱ进行估计，

θi={Fti|ti-1i,Htii,Qtii,Rtii}---(31)]]>

状态空间模型参数的修正应考虑突变点之前的监测数据对突变点之后状态估计产生的误差影响，同时应设法加大突变点之后接收到的监测数据对预测模型的影响，使得状态估计比较快的收敛于部件退化的真实状态，并进行准确的退化状态预测，t_i时刻待估计状态初值其均值和方差分别为滤波估计误差初值为

根据检测到的突变点信息修正预测模型后，重新进行状态估计，一步预测值为：

X^t|t-1i=Ft|t-1iX^t-1i---(32)]]>

X^ti=X^t|t-1i+Kti[Yti-HtiX^t|t-1i]---(33)]]>

Kti=Pt|t-1iHtiT[HtiPt|t-1iHtiT+Rti]-1---(34)]]>

式中为突变点A后的预测的误差方差阵，P_tⁱ为突变点A后的滤波估计误差的方差阵：

Pt|t-1i=Ft|t-1iPt-1iFt|t-1iT+Qt-1i---(35)]]>

Pti=[I-KtiHti]Pt|t-1i[I-KtiHti]T+KtiRtiKtiT---(36)]]>

其中t_i时刻采用修正后的卡尔曼滤波增益，修正后的卡尔曼滤波增益记为由(30)式可知：

Kt0i′=Kt0igi(ti)---(37)]]>

步骤6、根据齿轮状态估计及已知的齿轮故障阈值进行齿轮剩余寿命预测：

利用修正后的状态空间方程可对系统的监测点之后的退化状态进行预测，通过预测到的退化状态值及已知的退化状态故障阈值可求解首次到达故障阈值的时间，进而以蒙特卡洛的分析方法求得在t_p时刻所预测的系统的平均剩余寿命u_N(t_p)为：

uN(tp)=[Σm=1N(Tfm-tp)]/N---(38)]]>

式中T_fm为第m次预测退化状态到达利用齿轮疲劳试验所得系统故障阈值的时间，N为蒙特卡洛仿真次数，用最小方差控制N的取值，即N满足：

(uN(tp)-uN-1(tp)uN(tp))2<ϵ---(39)]]>

其中ε的取值为一极小值；

N次蒙特卡洛仿真后利用式(38)可求得齿轮所预测的系统的平均剩余寿命u_N(t_p)及剩余寿命分布函数。