CN105740625A - Real time residual life prediction method of gear - Google Patents

Abstract

The invention discloses a real time residual life prediction method of a gear, and belongs to the technical field of mechanical reliability. The method comprises following implementation steps of 1, using a sensor to monitor degeneration of the power gear in real time; 2, carrying out feature extraction to the fatigue state of the gear; carrying out degeneration evaluation to the wear degeneration performance of the gear; 3, building a parameter variable gear degeneration state space model; 4, detecting a sudden change state point in the degeneration process of the gear; 5, correcting a prediction model according to the detected sudden change point information; carrying out state evaluation again; and 6, predicting the residual life of the gear according to the gear state evaluation and a gear fault threshold value. The method has the advantages that the accuracy of predicting the degeneration state and the real time residual life of the gear can be effectively improved, thus providing the basis for preventative maintenance of the gear.

Description

A kind of real-time method for predicting residual useful life of gear

Technical field

The invention belongs to Mechanical Reliability field, be specifically related to the real-time method for predicting residual useful life of a kind of gear.

Background technology

Gear is the critical component in mechanical industry in widely used rotating machinery drive system.Gear engaging stress situation is complex, runs when variable working condition, varying load, it is easier to break down.When the fault such as gear generation broken teeth, the flank of tooth be tired, glued, the catastrophic destruction of whole plant equipment usually can be caused.Simultaneously along with the development of information sensing equipment, the running status of gear can be monitored in real time, utilize degenerate state and the residual life thereof of a large amount of monitoring information more Accurate Prediction systems in real time received, can be provided with the key message closing health status, and then identify and the management generation of fault, Maintenance Planning activity, provide foundation for more reasonably formulating the maintenance and repair strategy based on state；

The wear degradation process of gear is a dynamic stochastic process.Existing gear wear condition prediction major part is all based on Monitoring Data and sets up a random degradation model to analyze its degenerative process.Also have some to propose multistage degenerate state forecast model, the degeneration of gear is divided into normal deterioration and two stages of incipient fault, it is believed that increase at incipient fault stage catagen speed, only seek its residual life in this stage forecast degenerate state.And the degenerative process of gear changes along with the change of apparatus of load and operating mode, beginning gear engagement stage and later stage abrasion boost phase are multistage degenerative processes complicated and changeable, it does not have the invariable stage, it is difficult to describe with a degradation model.During simultaneity factor is run, its degenerate state is difficult to directly measure, and the data that can monitor are often some indirect operation data relevant with state, for instance the signal such as vibration signal, temperature and pressure.And observe data and often because being subject to the interference of random noise, observation error occurs.Carry out estimating can solving this problem well with predicting to the state of system according to observation data configuration state-space model；

In gear degenerative process, often due to its degeneration wear process is sudden change, qualitative change by gradual change, quantitative development, causes that the material of system, structure etc. change, retrogressive mutation point occurs.Deterioration velocity after catastrophe point increases and then makes residual life generation large change.Catastrophe point is the starting point that degenerate state changes, and a large amount of residual life information that it comprises can provide useful information for the formulation of predicting residual useful life and preventive maintenance maintenance policy.

Summary of the invention

It is an object of the present invention to provide a kind of real-time method for predicting residual useful life of gear, prediction based on degenerate state Singularity detection and residual life is carried out joint study, catastrophe point information can be utilized faster the dynamically change of system to be tracked, improve gear degenerate state precision of prediction, can efficiently against the shortcoming existed in prior art.

The present invention is achieved in that and it is characterized in that including following enforcement step:

Step 1, gear-box install sensor, obtain characterize gear degenerate Real-time Monitoring Data；

Acceleration transducer is arranged on the bearing block position of main examination case, mounting temperature sensor in gear-box, installs noise transducer in the surface of main examination case,

Step 2, fatigue state to gear carry out feature extraction, and gear wear degraded performance is carried out slump evaluations；

Utilizing mean square amplitude that gear wear degraded performance is carried out slump evaluations, in each sampling time Δ t/length, the mean square amplitude of time series of discrete random signal is represented by:

$x_{rms}\left(\mathrm{\Δ}t\right)={\left(\frac{1}{n}\underset{i=1}{\overset{n}{\Σ}}{x}_i^2\right(t\left)\right)}^{1/2}---\left(1\right)$

In formula, Δ t is the employing time；N=F_s× Δ t, F_sFor sample frequency, n is sampling number, and ∑ represents summation, i ∈ (1,2,3....n), x_iT () is sampled value；

Step 3, set up the gear degenerate state spatial model of variable element:

Carrying out state estimation by non-linear state space modeling, the general type of model is:

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

Y_t=h (X_t)+v_t(3)

In formula, X_tFor the system mode vector of t, t >=1, f (X_t-1) for systematic state transfer model；Y_tFor the systematic observation vector of t, h (X_t) for systematic observation model, system noise w_tWith observation noise v_tObeying the Gauss distribution of zero-mean, its covariance matrix is: Var (w_t)=Q_t, Var (v_t)=R_t；

Adopting linearisation to simplify, the state-space model after linearisation is:

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

Y_t=H_tX_t+v_t(5)

(4) formula represents with a parameter vector θ with the unknown parameter in (5) formula, i.e. θ={ F_t|t-1,H_t,Q_t,R_t, wherein, F_t|t-1For systematic state transfer matrix, H_tFor systematic observation matrix, the state vector of system is carried out i observation, thus obtaining observation sequence Y_1:i={ Y₁,Y₂,…,Y_i, when system is in running, along with the change of external environment, model parameter changes, therefore can monitoring in real time when, utilize Monitoring Data to carry out the estimation of state-space model parameter and the prediction of degenerate state, try to achieve at t system mode vector X_tBest estimator

Initially with expectation-maximization algorithm (Expectation Maximization, EM) from initial time real-time estimated state spatial model parameter in an iterative manner, this algorithm includes E step and M step two step, specifically comprises the following steps that

E step is to estimate in unknown parameter vector θ process, by the state initial value X to be estimated in Kalman filtering₀As unknown " missing data ", its average and variance respectively a₀,The unknown, filtering estimation difference initial value is P₀, the initial value of parameter to be estimated is:

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

For making degenerate state Initial value choice more reasonable, the method that the filtering of Kalman's forward direction and fixed interval inverse filtering combine is adopted to choose the degenerate state initial value of optimum while parameter estimation,

During the filtering of Kalman's forward direction, at t₁,t₂,...,t_i,t_i+1,...,t_i+nIn the moment, its corresponding degenerate state estimated value is: ${\hat{X}}_t={\hat{X}}_1,{\hat{X}}_2,\mathrm{...},{\hat{X}}_i,{\hat{X}}_{i+1},\mathrm{...},{\hat{X}}_{i+n},$ One-step prediction valueFor:

${\hat{X}}_{t|t-1}=F_{t|t-1}{\hat{X}}_{t-1}---\left(7\right)$

The state estimation that correction predicted estimate obtains is:

${\hat{X}}_t={\hat{X}}_{t|t-1}+K_t\[Y_t-H_t{\hat{X}}_{t|t-1}\]---\left(8\right)$

$K_t=P_{t|t-1}{H}_t^T{\[H_t{P}_{t|t-1}{H_t}^T+R_t\]}^{-1}---\left(9\right)$

K in formula_tFor filtering gain, P_t|t-1For the error covariance matrix of prediction, P_tFor filtering variance of estimaion error battle array:

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

$P_t=\[I-K_t{H}_t\]P_{t|t-1}{\[I-K_t{H}_t\]}^T+K_t{R}_t{K}_t^T---\left(11\right)$

Known to t_i+nIn the moment, receiving the monitoring information of i+n monitoring point, trying to achieve its logarithm maximum likelihood function is:

$\ln L\left(\mathrm{\θ}\right)=\ln \[\underset{t=t_0}{\overset{t_{i+n}}{\Π}}p\left(Y_t\right|X_t,\mathrm{\θ})\]=\underset{t=t_1}{\overset{t_{i+n}}{\Σ}}\ln p\left(Y_t\right|\mathrm{\θ})-\ln p\left(X_0\right|\mathrm{\θ})\underset{t=t_1}{\overset{t_{i+n}}{\Σ}}\ln p\left(X_t\right|X_{t-1},\mathrm{\θ})---\left(12\right)$

In formula, ∏ represents quadrature, a₀For X₀Average,For X₀Variance, the expression formula that can try to achieve likelihood function is:

$\begin{array}{c}\ln L\left(\mathrm{\θ}\right)=-\left\{\ln \right|{\mathrm{\σ}}_0^2|+\frac{{(x_0-a_0)}^2}{{\mathrm{\σ}}_0^2}+(i+n)\ln |Q_{t_{i+n}}|\\ +\underset{t=t_1}{\overset{t_{i+n}}{\Σ}}\frac{{(x_t-F_{t_{i+n}|t_{i+n-1}}{x}_{t-1})}^2}{Q_{t_{i+n}}}+(i+n)\ln \left|R_{t_{i+n}}\right|\\ +\underset{t=t_1}{\overset{t_{i+n}}{\Σ}}\frac{{(Y_t-H_{t_{i+n}}{X}_t)}^2}{R_{t_{i+n}}}\}\end{array}---\left(13\right)$

Adopt expectation-maximization algorithm to estimate parameter, the expected value of calculating formula (13), be designated as

$\begin{array}{c}{E}_{X|Y,{\mathrm{\θ}}_k}\[\ln L\left(\mathrm{\θ}\right)\]=E\{-\[\ln |{\mathrm{\σ}}_0^2|+\frac{{(x_0-a_0)}^2}{{\mathrm{\σ}}_0^2}+(i+n)\ln |Q_{t_{i+n}}|\\ +\underset{t=t_1}{\overset{t_{i+n}}{\Σ}}\frac{{(x_t-F_{t_{i+n}|t_{i+n-1}}{x}_{t-1})}^2}{Q_{t_{i+n}}}+(i+n)\ln \left|R_{t_{i+n}}\right|\\ +\underset{t=t_1}{\overset{t_{i+n}}{\Σ}}\frac{{(Y_t-H_{t_{i+n}}{X}_t)}^2}{R_{t_{i+n}}}\]\}\end{array}---\left(14\right)$

The smooth algorithm of Kalman is the process seeking expected value, first tries to achieve:

$E_{X|Y,{\mathrm{\θ}}_k}\left(X_t\right)={\hat{X}}_{t|T}---\left(15\right)$

$S_{t|T}=E_{X|Y,{\mathrm{\θ}}_k}\left({X_t}^2\right)={\hat{X}}_{t|T}^2+P_{t|T}---\left(16\right)$

$S_{t,t-1|T}=E_{X|Y,{\mathrm{\θ}}_k}\left(X_t{X}_{t-1}\right)={\hat{X}}_{t|T}{\hat{X}}_{t-1|T}+J_{t-1}{P}_{t|T}---\left(17\right)$

Wherein T=t_i+n,t_i+n-1,...,t_i,t_i-1,...,t₂,t₁, Y_TRepresent all monitoring informations,With P_t|TCan be tried to achieve by fixed interval Kalman's inverse filtering:

$J_t=P_t{F}_{t|t-1}^T{\left(P_{t|t-1}\right)}^{-1}---\left(18\right)$

${\hat{X}}_{t|T}={\hat{X}}_t+J_t\[{\hat{X}}_{t+1|T}-{\hat{X}}_{t+1|t}\]---\left(19\right)$

$P_{t|T}=P_t-J_t(P_{t+1|T}-P_{t+1|t}){J_t}^T---\left(20\right)$

(15) formula-(20) formula is substituted into (14) formula, can try to achieve:

$\begin{array}{c}{E}_{X|Y,{\mathrm{\θ}}_k}\[\ln L\left(\mathrm{\θ}\right)\]=-\left\{\ln \right|{\mathrm{\σ}}_0^2|+\frac{S_{t_0|T}-2a_0{\hat{X}}_{t|T}+a_0^2}{{\mathrm{\σ}}_0^2}+(i+n)\ln |Q_{t_{i+n}}|\\ +\underset{t=t_1}{\overset{t_{i+n}}{\Σ}}\frac{S_{t|T}+F_{t_{i+n}|t_{i+n-1}}^2{S}_{t-1|T}-2F_{t_{i+n}|t_{i+n-1}}{S}_{t,t-1|T}}{Q_{t_{i+n}}}+(i+n)\ln \left|R_{t_{i+n}}\right|\\ +\underset{t=t_1}{\overset{t_{i+n}}{\Σ}}\frac{Y_t^2-2H_{t_{i+n}}{Y}_t{\hat{X}}_{t|T}+H_{t_{i+n}}^2{S}_{t|T}}{R_{t_{i+n}}}\}\end{array}---\left(21\right)$

M step is the estimated value trying to achieve each parameter by solving maximum-likelihood estimation:

${\widehat{\mathrm{\θ}}}_{k+1}=\arg \underset{\mathrm{\θ}}{\max }\{E_{X|Y,{\mathrm{\θ}}_k}\[\ln L\left(\mathrm{\θ}\right)\]\}---\left(22\right)$

Expectation-maximization algorithm constantly carries out the iteration of E step and M step, until likelihood function value stabilization no longer changes, namely solves t_i+nThe estimates of parameters in moment

Step 4, the mutation status point in gear degenerative process is detected:

In gear degenerative process, owing to retrogressive mutation point occurs in its degeneration wear process, cumlative chart method is utilized to carry out Singularity detection, to the accumulative summation of the result repeatedly observed to judge whether monitoring process is in state in cont rol, detect mutation status point, then recurrence formula is:

$S_n^{+}=max\underset{i=1}{\overset{n}{\Σ}}\[y_i-{\mathrm{\μ}}_0-k\]=S_{n-1}^{+}+max(y_n-{\mathrm{\μ}}_0-k)---\left(23\right)$

Wherein Represent the detection accumulation that offsets up of average and, y_iFor the observation in i moment, y_nFor the observation in n moment, h is mutation status detection threshold value,Δ μ is minimum early warning amount of degradation, and sample average is μ₀, variance isWhenTime, system normal operation is described, it does not have catastrophe point；WhenTime mean that the abnormal catastrophe point of system,

Abrupt climatic change threshold value can adopt average run-length ARL approximate formula to determine:

$ARL=\frac{\exp (-c)+c-1}{2{\left(\mathrm{\Δ}\mathrm{\μ}\right)}^2}---\left(24\right)$

In formula, parameter_cFor:

C=-(Δ μ) (| h |/σ₀+ 1.166)=-2k (| h |/σ₀+1.166)(25)

The catastrophe point Information revision forecast model that step 5, basis detect, re-starts state estimation；

As in figure 2 it is shown, degenerate state is not undergone mutation, after A point, system degradation state is I curve.But after state mutation point A, deterioration velocity increases, and actual degenerate state is III curve.System receives a large amount of monitoring information, utilizes accumulation summation algorithm state mutation point A to be detected, and although as the passage of time, increasing of Real-time Monitoring Data, the precision that filtering is estimated is more and more higher.Trend with curve II is stablized and converges on actual value by Kalman filtering value gradually.But in the process of convergence, the degenerate state after catastrophe point be cannot be carried out quick tracking, it may appear that state estimation forecast error Δ x, affect the prediction of real-time residual life and the formulation of system maintenance maintenance decision；

Gear degenerate state forecast model is carried out following correction:

Utilize the information that catastrophe point provides to carry out after catastrophe point while system model parameters revision, reselect filtering initial value so that status predication value compares with non-correction model faster close to estimative time of day；

At catastrophe point correction filtering gain K_t, by big filtering gain K_tState estimation can be modified；

Increase the monitoring information after catastrophe point, update the system model parameter during state estimation, carry out adaptive state estimation；

When at t_i+nMoment detects t_iMoment is degenerate state catastrophe point, and gear degenerate state forecast model is carried out following correction:

After catastrophe point, first the correction of state-space model is that catastrophe point place degenerate state initial value is modified.Catastrophe point and t_iMoment filtering estimation difference initial value is:Utilize from initial time to t_i+nAll monitor values received by moment adopt (6) formula-(22) formula expectation-maximization algorithm parameter estimation algorithm to try to achieve estimates of parameters and be designated as:

${\mathrm{\θ}}_{t_{i+n}}=\{F_{t_{i+n}|t_{i+n-1}},H_{t_{i+n}},Q_{t_{i+n}},R_{t_{i+n}}\}---\left(26\right)$

When at t_i+nMoment detects t_iThere is catastrophe point A in the moment, need to revise the Kalman filter model initial value after catastrophe point A and parameter,

Catastrophe point t_iAfterwards, time t is designated as again:Its corresponding degenerate state estimated value isThe degenerate state value monitored after utilizing catastrophe point reappraises parameter and is designated as θⁱ:

Utilize t_iMoment and all monitoring informations being previously received adopt the estimates of parameters that EM parameter estimation algorithm is tried to achieve to be designated as:

${\mathrm{\θ}}_{t_i}=\{F_{t_i|t_{i-1}},H_{t_i},Q_{t_i},R_{t_i}\}---\left(27\right)$

t_iAfter moment correction, initial parameter values to be estimated is designated asTherefore the initial value correction of Kalman filtering parameter to be estimated is designated as:

${\mathrm{\θ}}_0^i=\{F_0^i,H_0^i,Q_0^i,R_0^i\}={\mathrm{\β\θ}}_{t_{i+n}}+{\mathrm{\γ\θ}}_{t_i}---\left(28\right)$

Wherein β and γ is the parameters revision factor, and β+γ=1.Select γ < β, strengthened the weight of data after catastrophe point by the parameters revision factor, the weight of the corresponding data reduced before catastrophe point,

After catastrophe point being detected, increase the gain matrix K after catastrophe point by adaptive gain function_tModulus value so that state estimation is than the time of day converging on component degradation faster.Assume t to be detected_iPlace's degenerate state is undergone mutation, and constructs an adaptive gain function g_i(t_i) it is:

$g_i\left(t_i\right)=\frac{1}{K_{t_i}\{1+\exp \&lsqb;-{\mathrm{at}}_i/b\&rsqb;\}}---\left(29\right)$

Wherein a and b is the parameter of adaptive gain function, and the method increasing filtering gain is with the optimality sacrificing filtering for cost, and therefore the selection of parameter a and the b of adaptive gain function should consider the problem of its convergence rate and filtering optimality simultaneously,

Then revised Kalman gainIt is denoted as:

$K_{t_i}^{\&prime;}=K_{t_i}{g}_i\left(t_i\right)---\left(30\right)$

t_iAfter the initial parameter values correction to be estimated in moment, utilize t_iN monitoring information after moment, is denoted asRepeat above-mentioned Kalman filtering and parameter estimation procedure, by the method for adaptive-filtering, constantly to t_iSystem model parameter unknown after moment and noise statistics θⁱEstimate,

${\mathrm{\&theta;}}^i=\{F_{t_i|t_{i-1}}^i,H_{t_i}^i,Q_{t_i}^i,R_{t_i}^i\}---\left(31\right)$

The correction of state-space model parameter is considered as the error impact that state estimation after catastrophe point is produced by the Monitoring Data before catastrophe point, should manage the impact strengthening the Monitoring Data received after catastrophe point to forecast model simultaneously.Make state estimation than the time of day converging on component degradation faster, and carry out degenerate state prediction accurately.t_iMoment state initial value to be estimatedIts average and variance are respectivelyFiltering estimation difference initial value is

After the catastrophe point Information revision forecast model detected, re-start state estimation, one-step prediction valueFor:

${\hat{X}}_{t|t-1}^i=F_{t|t-1}^i{\hat{X}}_{t-1}^i---\left(32\right)$

${\hat{X}}_t^i={\hat{X}}_{t|t-1}^i+K_t^i\&lsqb;Y_t^i-H_t^i{\hat{X}}_{t|t-1}^i\&rsqb;---\left(33\right)$

$K_t^i=P_{t|t-1}^i{H}_t^{iT}{\&lsqb;H_t^i{P}_{t|t-1}^i{H}_t^{iT}+R_t^i\&rsqb;}^{-1}---\left(34\right)$

In formulaFor the error covariance matrix of the prediction after catastrophe point A,Filtering variance of estimaion error battle array for after catastrophe point A:

$P_{t|t-1}^i=F_{t|t-1}^i{P}_{t-1}^i{F_{t|t-1}^i}^T+Q_{t-1}^i---\left(35\right)$

$P_t^i=\&lsqb;I-K_t^i{H}_t^i\&rsqb;P_{t|t-1}^i{\&lsqb;I-K_t^i{H}_t^i\&rsqb;}^T+K_t^i{R}_t^i{K}_t^{iT}---\left(36\right)$

Wherein t_iMoment adopts revised Kalman filtering gain, and revised Kalman filtering gain is designated asFrom (30) formula:

${K_{t_0}^i}^{\&prime;}=K_{t_0}^i{g}_i\left(t_i\right)---\left(37\right)$

Step 6, according to gear condition estimate and known gear distress threshold value carry out gear predicting residual useful life:

Utilize revised state space equation can the degenerate state after the monitoring point of system be predicted, can solve, by the degenerate state value predicted and known degenerate state fault threshold, the time arriving fault threshold first.And then try to achieve at t with the analysis method of Monte Carlo_pThe average remaining lifetime u of the system that the moment is predicted_N(t_p) it is:

$u_N\left(t_p\right)=\&lsqb;\underset{m=1}{\overset{N}{\&Sigma;}}(T_{fm}-t_p)\&rsqb;/N---\left(38\right)$

T in formula_fmBeing that the m time prediction degenerate state arrives the time utilizing gear fatigue test gained system failure threshold value, N is Monte Carlo simulation number of times, and by the value of LMS control N, namely N meets:

${\left(\frac{u_N\left(t_p\right)-u_{N-1}\left(t_p\right)}{u_N\left(t_p\right)}\right)}^2<\mathrm{\&epsiv;}---\left(39\right)$

Wherein the value of ε is a minimum；

Utilize formula (38) that the average remaining lifetime u of the system that gear is predicted can be tried to achieve after n times Monte Carlo simulation_N(t_p) and residual life distribution function.

The invention have the advantages that

The present invention proposes a kind of real-time method for predicting residual useful life of gear, the present invention is directed to gear wear degenerative process and set up state space forecast model, utilize the gear real-time monitoring vibration information real-time update model parameter received, mutation status point in degenerative process is detected simultaneously, and the life information provided according to catastrophe point adopts the filtering of Kalman's forward direction and smoothing algorithm constantly state-space model parameter to be modified while filtering in conjunction with expectation maximization parameter estimation algorithm, change the filter effect after retrogressive mutation, carry out real-time status estimation and biometry.The dynamically change of system can be tracked after utilizing catastrophe point information that forecast model is modified by the present invention faster, improves prediction gear degenerate state and the accuracy of real-time residual life.

Accompanying drawing explanation

Fig. 1 is the real-time method for predicting residual useful life flow chart of embodiment of the present invention middle gear；

Fig. 2 is an embodiment of the present invention degenerate state forecast model correction map；

Fig. 3 is the mean square amplitude curve figure of characteristics extraction of No. 4 sensor reception data in the embodiment of the present invention；

Accumulation and statistic variation diagram when Fig. 4 is report to the police in the embodiment of the present invention；

Fig. 5 is that the degenerate state do not predicted in the same time before and after forecast model correction in the embodiment of the present invention compares, Fig. 5 (a) is prediction degenerate state comparison diagram before and after 69 hours place's Modifying model, Fig. 5 (b) is prediction degenerate state comparison diagram before and after 70 hours place's Modifying model, and Fig. 5 (c) is 73 hours degenerate state comparison diagrams do not predicted in the same time before and after Modifying model；Fig. 5 (d) is 76 hours degenerate state comparison diagrams do not predicted in the same time before and after Modifying model；

Fig. 6 is by not predicted before and after Modifying model in the embodiment of the present invention that residual life probability density function compares in the same time；Fig. 6 (a) is that before within 70 hours, revising, residual life probability density function figure, Fig. 6 (b) are that after within 70 hours, revising, residual life probability density function figure, Fig. 6 (c) are residual life probability density function figure before revising for 73 hours；Fig. 6 (d) is that after within 73 hours, revising, residual life probability density function figure, Fig. 6 (e) are that before within 76 hours, revising, residual life probability density function figure, Fig. 6 (f) are residual life probability density function figure after revising for 76 hours.

Detailed description of the invention

Below in conjunction with accompanying drawing, an embodiment of the present invention is described further:

In the embodiment of the present invention, the real-time method for predicting residual useful life of gear, method flow diagram as shown in Figure 1: comprise the following steps:

Step 1, gear is carried out fatigue test, obtains and characterize the Real-time Monitoring Data that gear is degenerated:

Gear fatigue life test have employed power stream blocking test stand.The centre-to-centre spacing of testing stand is 150mm, and motor speed is 1200r/min.Casing vibration, oil temperature and noise etc. are monitored by process of the test.Adopting material in test is steel alloy, and tooth face hardness is the hardened face gear of 58-61HRC, and surface treatment is carburizing and quenching.Adopt positive and negative staggered overlap joint engagement system.Main examination case module m=3, the number of teeth is z1=z2=50, pressure angle α=20 °, facewidth 29mm, the real work facewidth 13～14mm；Accompanying examination case number of gear teeth is z3=z4=24, pressure angle α=20 °, and work facewidth 20mm.Gear mesh frequency has 2, and main examination case gear is 1000Hz, and accompanying examination case gear is 480Hz.Test lubricating oil adopts L-CKC320 Industrial Closed gear oil；

11 sensors are arranged in test altogether, and acceleration transducer is arranged on the bearing block position of main examination case, mounting temperature sensor in gear-box, installs noise transducer in the surface of main examination case.No. 1～No. 8 sensors are acceleration transducer, and No. 1～No. 4 transducer arrangements are in the radial direction of main examination axle box bearing seat, and axial at main examination case of No. 7 and No. 8 transducer arrangements, No. 5 and No. 6 transducer arrangements are in the radial direction accompanying examination axle box bearing seat；No. 9 and No. 10 sensors are sonic transducer, are arranged in main examination case and accompany directly over examination case about 40cm place；No. 11 sensors are the temperature sensors of test lubricating oil temperature, are arranged in main examination casing, test lubricating oil temperature in test.Sample frequency 25.6kHz, sampling time 60s, sampling interval 9min, adopt the mode of conventional method and dead load in groups to carry out in this test, torque is 822.7N.M, judges this gear failure when testing gear generation broken teeth；

Step 2, fatigue state to gear carry out feature extraction, and gear wear degraded performance is carried out slump evaluations；

Select when carrying out biometry detachment tooth position in gear fatigue test recently and carry out feature extraction at 463 groups of vibration signals of No. 4 sensors outputs of bearing block location arrangements.Sample frequency 25.6kHz in test, sampling time 60s, sampling interval 9min, be the monitoring time by the conversion of sampled point number of times, after removing singular point, formula (1) just can obtain the eigenvalue of shown in Fig. 3 No. 4 sensor output with monitoring time changing curve figure.

In Fig. 3, No. 4 mean square amplitude curve of sensor overall variation trend after removing indivedual singular points can reflect the relation that test gear each monitoring point abrasion condition is corresponding with vibrational energy.This curve contains from starting running in stage to completing the fatigue test mean square amplitude of vibration signal at 77.2h generation broken teeth, mean square amplitude y during by gear fatigue test known broken teeth^*(T^*)=77.375, are designated as degenerationization threshold value, T^*For arriving the time of fault threshold, T first^*=77.3h；

Step 3, set up the gear degenerate state spatial model of variable element；Setting up in gear fatigue test state-space model between observation and the system mode of No. 4 sensor outputs is (4) formula and (5) formula,

Wherein X_tFor the gear-like state value of t, Y_tObservation for No. 4 sensor outputs of t.The expectation maximization iterative algorithm utilizing (6) formula-(22) formula carries out model parameter estimation, and the observation received is different, selects initial value F_1|0=[1], H₁=[1], Q₀=[1], R₀=[1], a₀=0,Initial error P is estimated in filtering₀=100.In iterative process, when the difference of adjacent twice maximum likelihood function is less than 10^-5Time, iteration stopping；

Mutation status point in gear degenerative process is detected by step 4；

Gear brings into operation to enter through engagement stage and carries out Singularity detection after the gear normal deterioration stage, selects false alarm rate η=0.25%, and average run length ARL=400, taking minimum early warning amount of degradation is Δ μ=20, degeneration mean μ₀=51.9520, variances sigma₀ ²=3.0440, sudden change threshold value of warning can be solved abrupt climatic change threshold value of warning h=59.7023 by (23) formula-(25) formula, and when obtaining reporting to the police, accumulation and statistic are as shown in Figure 4.

In Fig. 4, A point is catastrophe point, and B point is early-warning point, detects that the catastrophe point A point of degenerate state occurs at 66.3h place, and B point exceedes threshold value of warning h report to the police in the accumulation of 69h place and statistic, owing to gear wear degenerate state changes,

Step 5, according to the catastrophe point Information revision forecast model detected, re-starts state estimation；

Utilize from initial time to detecting that all monitor values received for catastrophe point 69h adopt EM parameter estimation algorithm to try to achieve estimates of parameters and are designated asCatastrophe point 66.3h place and all monitoring informations of being previously received are utilized to adopt the estimates of parameters that EM parameter estimation algorithm is tried to achieve to be designated as

When at t_i+nMoment detects t_iMoment is degenerate state catastrophe point, and gear degenerate state forecast model is carried out following correction:

After catastrophe point, first the correction of state-space model is that catastrophe point place degenerate state initial value is modified.Catastrophe point and t_iMoment filtering estimation difference initial value is:Utilize from initial time to t_i+nAll monitor values received by moment adopt (6) formula-(22) formula expectation maximization parameter estimation algorithm to try to achieve estimates of parameters and be designated as:

${\mathrm{\&theta;}}_{t_{i+n}}=\{F_{t_{i+n}|t_{i+n-1}},H_{t_{i+n}},Q_{t_{i+n}},R_{t_{i+n}}\}---\left(40\right)$

When at t_i+nMoment detects t_iThere is catastrophe point A in the moment, need to revise the Kalman filter model initial value after catastrophe point A and parameter,

Catastrophe point t_iAfterwards, time t is designated as again:Its corresponding degenerate state estimated value isThe degenerate state value monitored after utilizing catastrophe point reappraises parameter and is designated as θⁱ:

Utilize t_iMoment and all monitoring informations being previously received adopt the estimates of parameters that expectation maximization parameter estimation algorithm is tried to achieve to be designated as:

${\mathrm{\&theta;}}_{t_i}=\{F_{t_i|t_{i-1}},H_{t_i},Q_{t_i},R_{t_i}\}---\left(41\right)$

t_iAfter moment correction, initial parameter values to be estimated is designated asTherefore the initial value correction of Kalman filtering parameter to be estimated is designated as:

${\mathrm{\&theta;}}_0^i=\{F_0^i,H_0^i,Q_0^i,R_0^i\}={\mathrm{\&beta;\&theta;}}_{t_{i+n}}+{\mathrm{\&gamma;\&theta;}}_{t_i}---\left(42\right)$

Wherein β and γ is the parameters revision factor, and β+γ=1.Select γ < β, strengthened the weight of data after catastrophe point by the parameters revision factor, the weight of the corresponding data reduced before catastrophe point,

After catastrophe point being detected, increase the gain matrix K after catastrophe point by adaptive gain function_tModulus value so that state estimation is than the time of day converging on component degradation faster.Assume t to be detected_iPlace's degenerate state is undergone mutation, and constructs an adaptive gain function g_i(t_i) it is:

$g_i\left(t_i\right)=\frac{1}{K_{t_i}\{1+\exp \&lsqb;-{\mathrm{at}}_i/b\&rsqb;\}}---\left(43\right)$

Wherein a and b is the parameter of adaptive gain function, the method increasing filtering gain is with the optimality sacrificing filtering for cost, therefore the selection of parameter a and the b of adaptive gain function should consider the problem of its convergence rate and filtering optimality simultaneously, and embodiment of the present invention filtering gain parameter value is a=8 × 10^-2, b=2.

Then revised Kalman gainIt is denoted as:

$K_{t_i}^{\&prime;}=K_{t_i}{g}_i\left(t_i\right)---\left(44\right)$

t_iAfter the initial parameter values correction to be estimated in moment, utilize t_iN monitoring information after moment, is denoted asRepeat above-mentioned Kalman filtering and parameter estimation procedure, by the method for adaptive-filtering, constantly to t_iSystem model parameter unknown after moment and noise statistics θⁱEstimate,

${\mathrm{\&theta;}}^i=\{F_{t_i|t_{i-1}}^i,H_{t_i}^i,Q_{t_i}^i,R_{t_i}^i\}---\left(45\right)$

The correction of state-space model parameter is considered as the error impact that state estimation after catastrophe point is produced by the Monitoring Data before catastrophe point, should manage the impact strengthening the Monitoring Data received after catastrophe point to forecast model simultaneously.Make state estimation than the time of day converging on component degradation faster, and carry out degenerate state prediction accurately.t_iMoment state initial value to be estimatedIts average and variance are respectivelyFiltering estimation difference initial value is

After the catastrophe point Information revision forecast model detected, re-start state estimation, one-step prediction valueFor:

${\hat{X}}_{t|t-1}^i=F_{t|t-1}^i{\hat{X}}_{t-1}^i---\left(46\right)$

${\hat{X}}_t^i={\hat{X}}_{t|t-1}^i+K_t^i\&lsqb;Y_t^i-H_t^i{\hat{X}}_{t|t-1}^i\&rsqb;---\left(47\right)$

$K_t^i=P_{t|t-1}^i{H}_t^{iT}{\&lsqb;H_t^i{P}_{t|t-1}^i{H}_t^{iT}+R_t^i\&rsqb;}^{-1}---\left(48\right)$

In formulaFor the error covariance matrix of the prediction after catastrophe point A,Filtering variance of estimaion error battle array for after catastrophe point A:

$P_{t|t-1}^i=F_{t|t-1}^i{P}_{t-1}^i{F_{t|t-1}^i}^T+Q_{t-1}^i---\left(49\right)$

$P_t^i=\&lsqb;I-K_t^i{H}_t^i\&rsqb;P_{t|t-1}^i{\&lsqb;I-K_t^i{H}_t^i\&rsqb;}^T+K_t^i{R}_t^i{K}_t^{iT}---\left(50\right)$

Wherein t_iMoment adopts revised Kalman filtering gain, and revised Kalman filtering gain is designated asFrom (30) formula:

${K_{t_0}^i}^{\&prime;}=K_{t_0}^i{g}_i\left(t_i\right)---\left(51\right)$

Step 6, according to gear condition estimate and known gear distress threshold value carry out gear predicting residual useful life；

Utilize revised state space equation can the degenerate state after the monitoring point of system be predicted, can solve, by the degenerate state value predicted and known degenerate state fault threshold, the time arriving fault threshold first.And then utilize formula (38) to try to achieve at t with the analysis method of Monte Carlo_pThe average remaining lifetime u of the system that the moment is predicted_N(t_p), when determining Monte Carlo simulation times N in its Chinese style (39), the value of ε is ε=e^-10。

From above-described embodiment:

The degenerate state curve that known state-space model before and after not revising in the same time is predicted as shown in Figure 5 is all different, and within 69 hours, place's predicted state value deviation actual value is more, and error is bigger.Along with monitoring information increases, after filtering time is fully grown, the degenerate state that the state-space model before and after revising is predicted is all close to actual value, and error reduces.In identical prediction time, state-space model before and after relatively revising predicts the outcome, the gear degenerate state that revised state-space model is predicted closer to actual value, can well prognoses system degenerate state undergo mutation after the change of its degenerate state and time of failure.The feature that uncorrected state-space model predicts can not suddenly change well after degenerate state change, it was predicted that residual life error is bigger.By Fig. 5 (c) it can be seen that when mutation status point just being detected, owing to initial value just revised by revised state-space model, the state estimation that the recurrence calculation of beginning is obtained is relatively coarse, estimation difference is bigger, but state estimation ratio is faster near estimative time of day after a few step recursive operations, it is possible to follow the tracks of it well and dynamically change.

The residual life probability density function do not predicted in the same time before and after state-space model correction is tried to achieve as shown in Figure 6 by Monte Carlo method.Known receive the increase of monitoring information along with the passage of predicted time, at the residual life probability density function that 76h utilization state spatial model predicts, no matter be that before and after Modifying model, residual life probability density function peak value all moves closer to actual value.The average remaining lifetime predicted after same prediction time, state-space model correction is closer to real residual life, and its probability density function prediction variance compared with the residual life probability density function that the state-space model before revising is predicted reduces.

Forecasting Methodology for comparing the real-time residual life of a kind of gear of the present invention further predicts the accuracy of residual life, introduces the concept of prediction accuracy, it was predicted that accuracy PA is:

$PA=(1-\frac{|u_N\left(t_p\right)-T_a|}{T_a})\&times;100\%---\left(40\right)$

Wherein | u_N(t_p)-T_a| represent absolute value, u_N(t_p) represent t_pMoment predicts the average remaining lifetime tried to achieve, T_aFor gear real surplus life-span, T^*=77.2h is the gear physical fault time, then T_a=T^*-t_p, thus can calculate the comparison of its prediction accuracy before and after state-space model correction as shown in table 1:

Before and after table 1 Modifying model, prediction accuracy compares

After patent of the present invention is undergone mutation in gear running as shown in Table 1, utilize the information that catastrophe point provides that model is modified improving the prediction accuracy of real-time residual life in system operation.

Along with increasing of monitoring information, the predicted predicting residual useful life error of state-space model before and after revising all is gradually reduced, it was predicted that accuracy gradually steps up.But comparing at different future position 69h, 70h and 73h and 75h place, the residual life accuracy predicted after state-space model correction is above utilizing the state-space model before revising to be predicted the prediction accuracy tried to achieve.The residual life predicted after forecast model residual life moves closer to Modifying model is not revised at 75h place, it can thus be appreciated that patent of the present invention after system degradation performance is undergone mutation can Accurate Prediction gear residual life, and for adjusting the foundation that preventive maintenance maintenance policy provides necessary in time.

In sum, the present invention sets up system state space forecast model according to the real-time monitoring information of the gear received.Consider the relatedness of outlier detection and predicting residual useful life simultaneously, after mutation status point in gear wear degenerative process is detected, utilize the life information that degenerate state catastrophe point provides to adopt the filtering of Kalman's forward direction and smoothing algorithm to revise state-space model while filtering in conjunction with expectation maximization parameter estimation algorithm, carry out state estimation and biometry.Faster the dynamically change of system can be tracked by its revised forecast model of gear fatigue life verification experimental verification, improve real-time predicting residual useful life accuracy.

Claims (1)

1. the real-time method for predicting residual useful life of a gear, it is characterised in that implementing step is:

Step 1, gear-box install sensor, obtain characterize gear degenerate Real-time Monitoring Data；

Acceleration transducer is arranged on the bearing block position of main examination case, mounting temperature sensor in gear-box, installs noise transducer in the surface of main examination case；

Step 2, fatigue state to gear carry out feature extraction, and gear wear degraded performance is carried out slump evaluations:

Utilizing mean square amplitude that gear wear degraded performance is carried out slump evaluations, in each sampling time Δ t/length, the mean square amplitude of time series of discrete random signal is represented by:

$x_{rms}\left(\mathrm{\&Delta;}t\right)={\left(\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i^2\right(t\left)\right)}^{1/2}---\left(1\right)$

In formula, Δ t is the employing time；N=F_s× Δ t, F_sFor sample frequency, n is sampling number, and ∑ represents summation, i ∈ (1,2,3....n), x_iT () is sampled value；

Step 3, set up the gear degenerate state spatial model of variable element:

Carrying out state estimation by non-linear state space modeling, the general type of model is:

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

Y_t=h (X_t)+v_t(3)

In formula, X_tFor the system mode vector of t, t >=1, f (X_t-1) for systematic state transfer model；Y_tFor the systematic observation vector of t, h (X_t) for systematic observation model, system noise w_tWith observation noise v_tObeying the Gauss distribution of zero-mean, its covariance matrix is: Var (w_t)=Q_t, Var (v_t)=R_t

Adopting linearisation to simplify, the state-space model after linearisation is:

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

Y_t=H_tX_t+v_t(5)

(4) formula represents with a parameter vector θ with the unknown parameter in (5) formula, i.e. θ={ F_t|t-1,H_t,Q_t,R_t, wherein, F_t|t-1For systematic state transfer matrix, H_tFor systematic observation matrix, the state vector of system is carried out i observation, thus obtaining observation sequence Y_1:i={ Y₁,Y₂,…,Y_i, when system is in running, along with the change of external environment, model parameter changes, therefore can monitoring in real time when, utilize Monitoring Data to carry out the estimation of state-space model parameter and the prediction of degenerate state, try to achieve at t system mode vector X_tBest estimator

Initially with expectation-maximization algorithm from initial time real-time estimated state spatial model parameter in an iterative manner, this algorithm includes E step and M step two step, specifically comprises the following steps that

E step is to estimate in unknown parameter vector θ process, by the state initial value X to be estimated in Kalman filtering₀As unknown " missing data ", its average and variance are respectivelyThe unknown, filtering estimation difference initial value is P₀, the initial value of parameter to be estimated is:

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

For making degenerate state Initial value choice more reasonable, the method that the filtering of Kalman's forward direction and fixed interval inverse filtering combine is adopted to choose the degenerate state initial value of optimum while parameter estimation,

During the filtering of Kalman's forward direction, at t₁,t₂,...,t_i,t_i+1,...,t_i+nIn the moment, its corresponding degenerate state estimated value is:One-step prediction valueFor:

${\hat{X}}_{t|t-1}=F_{t|t-1}{\hat{X}}_{t-1}---\left(7\right)$

The state estimation that correction predicted estimate obtains is:

${\hat{X}}_t={\hat{X}}_{t|t-1}+K_t\&lsqb;Y_t-H_t{\hat{X}}_{t|t-1}\&rsqb;---\left(8\right)$

$K_t=P_{t|t-1}{H}_t^T{\&lsqb;H_t{P}_{t|t-1}{H_t}^T+R_t\&rsqb;}^{-1}---\left(9\right)$

K in formula_tFor filtering gain, P_t|t-1For the error covariance matrix of prediction, P_tFor filtering variance of estimaion error battle array:

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

$P_t=[I-K_t{H}_t]P_{t|t-1}{[I-K_t{H}_t]}^T+K_t{R}_t{K}_t^T---\left(11\right)$

Known to t_i+nIn the moment, receiving the monitoring information of i+n monitoring point, trying to achieve its logarithm maximum likelihood function is:

$\ln L\left(\mathrm{\&theta;}\right)=\ln \&lsqb;\underset{t=t_0}{\overset{t_{i+n}}{\mathrm{\&Pi;}}}p\left(Y_t\right|X_t,\mathrm{\&theta;})\&rsqb;=\underset{t=t_1}{\overset{t_{i+n}}{\mathrm{\&Pi;}}}\ln p\left(Y_t\right|\mathrm{\&theta;})-\ln p\left(X_0\right|\mathrm{\&theta;})\underset{t=t_1}{\overset{t_{i+n}}{\mathrm{\&Sigma;}}}\ln p\left(X_t\right|X_{t-1},\mathrm{\&theta;})---\left(12\right)$

In formula, ∏ represents quadrature, a₀For X₀Average,For X₀Variance, the expression formula that can try to achieve likelihood function is:

$\begin{array}{c}\ln L\left(\mathrm{\&theta;}\right)=-\left\{\ln \right|{\mathrm{\&sigma;}}_0^2|+\frac{{(x_0-a_0)}^2}{{\mathrm{\&sigma;}}_0^2}+(i+n)\ln |Q_{t_{i+n}}|\\ +\underset{t=t_1}{\overset{t_{i+n}}{\mathrm{\&Sigma;}}}\frac{{(x_t-F_{t_{i+n}|t_{i+n-1}}{x}_{t-1})}^2}{Q_{t_{i+n}}}+(i+n)\ln \left|R_{t_{i+n}}\right|\\ +\underset{t=t_1}{\overset{t_{i+n}}{\mathrm{\&Sigma;}}}\frac{{(Y_t-H_{t_{i+n}}{X}_t)}^2}{R_{t_{i+n}}}\}\end{array}---\left(13\right)$

Adopt expectation-maximization algorithm to estimate parameter, the expected value of calculating formula (13), be designated as

$\begin{array}{c}{E}_{X|Y,{\mathrm{\&theta;}}_k}\&lsqb;\ln L\left(\mathrm{\&theta;}\right)\&rsqb;=E\{-\&lsqb;\ln |{\mathrm{\&sigma;}}_0^2|+\frac{{(x_0-a_0)}^2}{{\mathrm{\&sigma;}}_0^2}+(i+n)\ln |Q_{t_{i+n}}|\\ +\underset{t=t_1}{\overset{t_{i+n}}{\mathrm{\&Sigma;}}}\frac{{(x_t-F_{t_{i+n}|t_{i+n-1}}{x}_{t-1})}^2}{Q_{t_{i+n}}}+(i+n)\ln \left|R_{t_{i+n}}\right|\\ +\underset{t=t_1}{\overset{t_{i+n}}{\mathrm{\&Sigma;}}}\frac{{(Y_t-H_{t_{i+n}}{X}_t)}^2}{R_{t_{i+n}}}\&rsqb;\}\end{array}---\left(14\right)$

The smooth algorithm of Kalman is the process seeking expected value, first tries to achieve:

$E_{X|Y,{\mathrm{\&theta;}}_k}\left(X_t\right)={\hat{X}}_{t|T}---\left(15\right)$

$S_{t|T}=E_{X|Y,{\mathrm{\&theta;}}_k}\left({X_t}^2\right)={\hat{X}}_{t|T}^2+P_{t|T}---\left(16\right)$

$S_{t,t-1|T}=E_{X|Y,{\mathrm{\&theta;}}_k}\left(X_t{X}_{t-1}\right)={\hat{X}}_{t|T}{\hat{X}}_{t-1|T}+J_{t-1}{P}_{t|T}---\left(17\right)$

Wherein T=t_i+n,t_i+n-1,...,t_i,t_i-1,...,t₂,t₁, Y_TRepresent all monitoring informations,With P_t|TCan be tried to achieve by fixed interval Kalman's inverse filtering:

$J_t={\mathrm{PF}}_{t|t-1}^T{\left(P_{t|t-1}\right)}^{-1}---\left(18\right)$

${\hat{X}}_{t|T}={\hat{X}}_t+J_t\&lsqb;{\hat{X}}_{t+1|T}-{\hat{X}}_{t+1|t}\&rsqb;---\left(19\right)$

$P_{t|T}=P_t-J_t(P_{t+1|T}-P_{t+1|t}){J_t}^T---\left(20\right)$

(15) formula-(20) formula is substituted into (14) formula, can try to achieve:

$\begin{array}{c}{E}_{X|Y,{\mathrm{\&theta;}}_k}\&lsqb;\ln L\left(\mathrm{\&theta;}\right)\&rsqb;=-\left\{\ln \right|{\mathrm{\&sigma;}}_0^2|+\frac{S_{t_0|T}-2a_0{\hat{X}}_{t|T}+a_0^2}{{\mathrm{\&sigma;}}_0^2}+(i+n)\ln |Q_{t_{i+n}}|\\ +\underset{t=t_1}{\overset{t_{i+n}}{\mathrm{\&Sigma;}}}\frac{S_{t|T}+F_{t_{i+n}|t_{i+n-1}}^2{S}_{t-1|T}-2F_{t_{i+n}|t_{i+n-1}}{S}_{t,t-1|T}}{Q_{t_{i+n}}}+(i+n)\ln \left|R_{t_{i+n}}\right|\\ +\underset{t=t_1}{\overset{t_{i+n}}{\mathrm{\&Sigma;}}}\frac{Y_t^2-H_{t_{i+n}}{Y}_t{\hat{X}}_{t|T}+H_{t_{i+n}}^2{S}_{t|T}}{R_{t_{i+n}}}\}\end{array}---\left(21\right)$

M step is the estimated value trying to achieve each parameter by solving maximum-likelihood estimation:

${\widehat{\mathrm{\&theta;}}}_{k+1}=\arg \underset{\mathrm{\&theta;}}{max}\{E_{X|Y,{\mathrm{\&theta;}}_k}\&lsqb;\ln L\left(\mathrm{\&theta;}\right)\&rsqb;\}---\left(22\right)$

Expectation-maximization algorithm constantly carries out the iteration of E step and M step, until likelihood function value stabilization no longer changes, namely solves t_i+nThe estimates of parameters in moment

Step 4, the mutation status point in gear degenerative process is detected:

In gear degenerative process, owing to retrogressive mutation point occurs in its degeneration wear process, cumlative chart method is utilized to carry out Singularity detection, to the accumulative summation of the result repeatedly observed to judge whether monitoring process is in state in cont rol, detect mutation status point, then recurrence formula is:

$S_n^{+}=max\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\&lsqb;y_i-{\mathrm{\&mu;}}_0-k\&rsqb;=S_{n-1}^{+}+max(y_n-{\mathrm{\&mu;}}_0-k)---\left(23\right)$

Wherein Represent the detection accumulation that offsets up of average and, y_iFor the observation in i moment, y_nFor the observation in n moment, h is mutation status detection threshold value,Δ μ is minimum early warning amount of degradation, and sample average is μ₀, variance isWhenTime, system normal operation is described, it does not have catastrophe point；WhenTime mean that the abnormal catastrophe point of system,

Abrupt climatic change threshold value can adopt average run-length ARL approximate formula to determine:

$ARL=\frac{\exp (-c)+c-1}{2{\left(\mathrm{\&Delta;}\mathrm{\&mu;}\right)}^2}---\left(24\right)$

In formula, parameter c is:

C=-(Δ μ) (| h |/σ₀+ 1.166)=-2k (| h |/σ₀+1.166)(25)

The catastrophe point Information revision forecast model that step 5, basis detect, re-starts state estimation:

When at t_i+nMoment detects t_iMoment is degenerate state catastrophe point, and gear degenerate state forecast model is carried out following correction:

After catastrophe point, first the correction of state-space model is that catastrophe point place degenerate state initial value is modified.Catastrophe point and t_iMoment filtering estimation difference initial value is:Utilize from initial time to t_i+nAll monitor values received by moment adopt (6) formula-(22) formula expectation-maximization algorithm parameter estimation algorithm to try to achieve estimates of parameters and be designated as:

${\mathrm{\&theta;}}_{t_{i+n}}=\{F_{t_{i+n}|t_{i+n-1}},H_{t_{i+n}},Q_{t_{i+n}},R_{t_{i+n}}\}---\left(26\right)$

When at t_i+nMoment detects t_iThere is catastrophe point A in the moment, need to revise the Kalman filter model initial value after catastrophe point A and parameter,

Catastrophe point t_iAfterwards, time t is designated as again:Its corresponding degenerate state estimated value isThe degenerate state value monitored after utilizing catastrophe point reappraises parameter and is designated as θⁱ:

Utilize t_iMoment and all monitoring informations being previously received adopt the estimates of parameters that expectation-maximization algorithm parameter estimation algorithm is tried to achieve to be designated as:

${\mathrm{\&theta;}}_{t_i}=\{F_{t_i|t_{i-1}},H_{t_i},Q_{t_i},R_{t_i}\}---\left(27\right)$

t_iAfter moment correction, initial parameter values to be estimated is designated asTherefore the initial value correction of Kalman filtering parameter to be estimated is designated as:

${\mathrm{\&theta;}}_0^i=\{F_0^i,H_0^i,Q_0^i,R_0^i\}={\mathrm{\&beta;\&theta;}}_{t_{i+n}}+{\mathrm{\&gamma;\&theta;}}_{t_i}---\left(28\right)$

Wherein β and γ is the parameters revision factor, and β+γ=1.Select γ < β, strengthened the weight of data after catastrophe point by the parameters revision factor, the weight of the corresponding data reduced before catastrophe point,

After catastrophe point being detected, increase the gain matrix K after catastrophe point by adaptive gain function_tModulus value so that state estimation is than the time of day converging on component degradation faster.Assume t to be detected_iPlace's degenerate state is undergone mutation, and constructs an adaptive gain function g_i(t_i) it is:

$g_i\left(t_i\right)=\frac{1}{K_{t_i}\{1+\exp \&lsqb;-{\mathrm{at}}_i/b\&rsqb;\}}---\left(29\right)$

Wherein a and b is the parameter of adaptive gain function, and the method increasing filtering gain is with the optimality sacrificing filtering for cost, and therefore the selection of parameter a and the b of adaptive gain function should consider the problem of its convergence rate and filtering optimality simultaneously,

Then revised Kalman gainIt is denoted as:

$K_{t_i}^{\&prime;}=K_{t_i}{g}_i\left(t_i\right)---\left(30\right)$

t_iAfter the initial parameter values correction to be estimated in moment, utilize t_iN monitoring information after moment, is denoted asRepeat above-mentioned Kalman filtering and parameter estimation procedure, by the method for adaptive-filtering, constantly to t_iSystem model parameter unknown after moment and noise statistics θⁱEstimate,

${\mathrm{\&theta;}}^i=\{F_{t_i|t_{i-1}}^i,H_{t_i}^i,Q_{t_i}^i,R_{t_i}^i\}---\left(31\right)$

The correction of state-space model parameter is considered as the error impact that state estimation after catastrophe point is produced by the Monitoring Data before catastrophe point, should manage the impact strengthening the Monitoring Data received after catastrophe point to forecast model simultaneously.Make state estimation than the time of day converging on component degradation faster, and carry out degenerate state prediction accurately.t_iMoment state initial value to be estimatedIts average and variance are respectivelyFiltering estimation difference initial value is

After the catastrophe point Information revision forecast model detected, re-start state estimation, one-step prediction valueFor:

${\hat{X}}_{t|t-1}^i=F_{t|t-1}^i{\hat{X}}_{t-1}^i---\left(32\right)$

${\hat{X}}_t^i={\hat{X}}_{t|t-1}^i+K_t^i\&lsqb;Y_t^i-H_t^i{\hat{X}}_{t|t-1}^i\&rsqb;---\left(33\right)$

$K_t^i=P_{t|t-1}^i{H}_t^{iT}{\&lsqb;H_t^i{P}_{t|t-1}^i{H}_t^{iT}+R_t^i\&rsqb;}^{-1}---\left(34\right)$

In formulaFor the error covariance matrix of the prediction after catastrophe point A,Filtering variance of estimaion error battle array for after catastrophe point A:

$P_{t|t-1}^i=F_{t|t-1}^i{P}_{t-1}^i{F_{t|t-1}^i}^T+Q_{t-1}^i---\left(35\right)$

$P_t^i=\&lsqb;I-K_t^i{H}_t^i\&rsqb;P_{t|t-1}^i{\&lsqb;I-K_t^i{H}_t^i\&rsqb;}^T+K_t^i{R}_t^i{K}_t^{iT}---\left(36\right)$

Wherein t_iMoment adopts revised Kalman filtering gain, and revised Kalman filtering gain is designated asFrom (30) formula:

${K_{t_0}^i}^{\&prime;}=K_{t_0}^i{g}_i\left(t_i\right)---\left(37\right)$

Step 6, according to gear condition estimate and known gear distress threshold value carry out gear predicting residual useful life:

Utilize revised state space equation can the degenerate state after the monitoring point of system be predicted, can solve, by the degenerate state value predicted and known degenerate state fault threshold, the time arriving fault threshold first.And then try to achieve at t with the analysis method of Monte Carlo_pThe average remaining lifetime u of the system that the moment is predicted_N(t_p) it is:

$u_N\left(t_p\right)=\&lsqb;\underset{m=1}{\overset{N}{\mathrm{\&Sigma;}}}(T_{fm}-t_p)\&rsqb;/N---\left(38\right)$

T in formula_fmBeing that the m time prediction degenerate state arrives the time utilizing gear fatigue test gained system failure threshold value, N is Monte Carlo simulation number of times, and by the value of LMS control N, namely N meets:

${\left(\frac{u_N\left(t_p\right)-u_{N-1}\left(t_p\right)}{u_N\left(t_p\right)}\right)}^2<\mathrm{\&epsiv;}---\left(39\right)$

Wherein the value of ε is a minimum,

Utilize formula (38) that the average remaining lifetime u of the system that gear is predicted can be tried to achieve after n times Monte Carlo simulation_N(t_p) and residual life distribution function.