CN104506162A - Fault prognosis method for high-order particle filter on basis of LS-SVR (least squares support vector regression) modeling - Google Patents

Abstract

The invention provides a mechanical fault prognosis method for a high-order particle filter on the basis of LS-SVR (least squares support vector regression) modeling and aims to establish a state equation in a data driven mode by extracting signal characteristics of evolving into a fault state from a normal state and implement online real-time prediction. According to the mechanical fault prognosis method, corresponding signal characteristics, such as six statistical characteristics of a time domain, are extracted for a provided mechanical vibration signal, and an LLE (Locally Linear Embedding) method is applied to carry out low-dimension characteristic extraction on sample characteristics; then an m-order HMM (Hidden Markov Model) is constructed by utilizing LS-SVR and is trained by low-dimension characteristics; the m-order HMM is fused into a particle filtering algorithm frame to form an m-order particle filtrer, and finally, a mechanical fault is predicted. Compared with an existing prognosis method, the mechanical fault prognosis method has the characteristics of high prognosis accuracy and wide application range.

Description

Based on the fault predicting method of the order particles filter of LS-SVR modeling

Technical field

The present invention relates to method for diagnosing faults field, specifically use least square method supporting vector machine (LS-SVM) to set up state equation in order particles filter by the mode of study, for the method for equipment fault indication.

Background technology

Maintenance (CBM) based on state instead of traditional Breakdown Maintenance shown according to schedule, become the preferred version that various engineering system implements maintenance, to ensure their reliability, fail safe and availability, CBM utilizes data running time, determining/predict machine state, thus determine current/following malfunction.The generation of machine failure and catastrophic discontinuityfailure is avoided by indication machine state.CBM actualizing technology mainly comprises sensing and monitoring, information processing, failure diagnosis and fault indication algorithm, accurately and timely can detect the useful life of initial failure and all the other assemblies broken down of prediction.Indication ability is important component part wherein, can prediction component or subsystem future and remaining life exactly.

Mechanical breakdown predicting method can be divided into two large classes: the method based on model and the method based on data-driven.Method based on model predicts the evolution trend of fault.The model of a given system, the method based on model can provide predicted estimate accurately.But usually in practical situations both, be difficult to design model accurately, especially when the process of fault propagation is complicated or does not understand completely.The another kind of method based on data-driven adopts the status data collected to set up fault-traverse technique.There is the method for a lot of data-driven, such as: Markov model, recurrence neuralward network, adaptive neuro-fuzzy inference system, support vector regression and Least square support vector regression (LS-SVR) etc. are all fallout predictors popular in mechanical breakdown indication field.

In the method for data-driven, fuzzy neural network and support vector regression achieve successful application in machine performance degradation prediction.But when machine performance dynamic realtime change in time, if do not consider the dynamic change of state at forecasting process, so fuzzy neural network and support vector regression method can not be made and predicting accurately.

In recent decades, in statistics and various engineering field, many experts and scholars are devoted to the research and analysis of the real-time estimation problem of dynamical system, the bayesian theory rising in the 17th century British vicar T.R.Bayes is that the state estimation problem of dynamical system provides strict theoretical frame, it utilizes all Given informations to carry out the posterior probability density function of tectonic system state variable, namely use the posterior probability of system model predictions state, recycle up-to-date measured value and revise.The various statistical values of state obtain as average, variance etc. all can calculate from posterior probability density function.Early 1990s, along with the rapid lifting of computer computation ability and memory space, a kind of real-time on-line simulation algorithm---particle filter combined based on recursion Bayesian Estimation and monte carlo method, is subject to people's attention gradually.Particle filter is a kind of method based on emulation, in its utilization state space, the random sample particle of one group of Weighted Coefficients approaches the probability density function of dbjective state variable, a possible state of each sample representation target, comprehensively all particle states can obtain the minimum variance estimate of dbjective state.This algorithm is not linear by model, the constraint of Gauss's hypothesis, is applicable to any nonlinear and non-Gaussian dynamical system.

In current disclosed document both at home and abroad, Zio E, Peloni G.Particle filtering prognostic estimationof the remaining useful life of nonlinear components.Reliability Engineering and SystemSafety 2011; The mechanical part method for estimating state based on particle filter is proposed in 96 (3): 403-409..The method application, based on the state equation of model and observational equation, in conjunction with the priori probability density function of each state of Monte Carlo simulation technology degradation estimation parts, avoids the hypothesis of Kalman Filter Technology to the simple linear and Gaussian noise of system.Provide the status predication framework that a robustness is stronger.The method have been applied in the estimation of the remaining life of fracture defect, drawn satisfied result.Zhang L, Li X S, Yu J S, et al.A fault prognostic algorithmbased on Gaussian mixture model particle filter.Acta Aeronautica et Astronautica Sinica2009; Propose a kind of prediction algorithm in 30 (2): 319-324, the state estimation stage of this algorithm, adopt Combined estimator and particle filter to estimate the Posterior distrbutionp of objective system failure evolution model state and unknown parameter simultaneously.In the status predication stage of algorithm, have employed two kinds of different computational methods: a kind of method carries out iteration sampling to the Posterior distrbutionp of state variable current time, thus obtain the prior distribution of the state variable of future time instance; Another kind method is the measurement information adopting the method for data-driven to predict objective system in following a period of time, thus the forecasting problem of the prior distribution of future time instance state variable is converted into the estimation problem that solves Posterior distrbutionp.Adopt gauss hybrid models approximate random variable distribution density, thus the result of calculation of two kinds of methods is carried out effectively mutual under a unified prediction framework, further increase accuracy and the reliability of prediction.In the decision phase of algorithm, in the failure evolution model state variable distributed basis obtained, go out the distribution of objective system residual life in conjunction with certain failure criterion approximate calculation.Failure predication the simulation experiment result demonstrates the validity of put forward algorithm.

But above-mentioned two kinds of Forecasting Methodologies have some not enough:

(1) method based on model needs to carry out mathematical modeling to research to picture, when the process of fault propagation is complicated or does not understand completely, is difficult to draw correct failure state model.

(2) because single order Markov form make use of a small amount of historical information, namely the state value of subsequent time only and current time state value have relation, the growth course describing fault can only be used for.But shape State evolution is not only relevant with the state of eve, and relevant with multiple states before.So application single order Markov model can not be real description Condition forecasting model.

Summary of the invention

In order to overcome the deficiencies in the prior art, the invention provides a kind of order particles filtering mechanical breakdown based on LS-SVR modeling indication new method, effectively to improve fault predicted precision.

Technical scheme of the present invention is:

The fault predicting method of described a kind of order particles filter based on LS-SVR modeling, is characterized in that: comprise the following steps:

Step 1: to one group of signal developed to malfunction from normal condition by time slip-window intercept signal, each segment signal intercepted is as a sample, calculate 6 kinds of Time-domain Statistics features of each sample, Time-domain Statistics feature comprises peak value, average amplitude, root-mean-square amplitude, degree of skewness index and kurtosis index; Then set up sample characteristics space, and with local linear embedding grammar, low-dimensional feature extraction is carried out to the sample characteristics space built, obtain low-dimensional feature set;

Step 2: using the low-dimensional feature set that obtains in step 1 as sample training and test fallout predictor, odd point wherein in low-dimensional feature set is as training set, even number point is as test set, and training set is used to obtain forecast model, and test set is used to the effect of the prediction testing forecast model;

Step 3: adopt training set training least square method supporting vector machine, obtain least square method supporting vector machine forecast model, least square method supporting vector machine forecast model is added that plant noise is combined into high-order HMM:

$x_k={\hat{x}}_k+v_{k+1}$

${\hat{x}}_{k}f(x_{k-1},x_{k-2},......,x_{k-m})$

Wherein f () is least square method supporting vector machine forecast model, and m is exponent number, x _kfor the state value in the k moment, for the state estimation in the k moment, x _k-1, x _k-2..., x _k-mfor 1 moment before the k moment is to the state value in front m moment, v _k-1for the plant noise in k-1 moment;

Step 4: produce particle collection according to the high-order HMM that step 3 obtains the particle number of particle set is N _s, the i-th particle represent the state set in 0 to k moment, the corresponding weights of each particle of particle set are weights by obtain, wherein by k moment dbjective state posterior probability density obtain, z _krepresent the measuring value in k moment; And by weights be normalized to $\underset{i}{\mathrm{\Σ}}{\mathrm{\ω}}_k^i=1;$

Step 5: judge sample degeneracy degree according to the weights that step 4 obtains, impose a condition if sample degeneracy arrives, then carry out particle resampling, and particle weights after calculating resampling;

Step 6: obtain state estimation by particle weights and particle value:

${\hat{x}}_k=\underset{i=1}{\overset{\mathrm{Ns}}{\mathrm{\Σ}}}{\mathrm{\ω}}_k^i{x}_k^i$

Again state estimation is substituted into forecast model obtain state value;

Step 7: carry out subsequent time computing to step 6 according to step 4, and circulation is carried out until mechanical breakdown has indicated.

Further preferred version, the fault predicting method of described a kind of order particles filter based on LS-SVR modeling, is characterized in that: in step 6, obtains state estimation by particle weights and particle value by following multistep status predication:

${\hat{x}}_k=\underset{i=1}{\overset{\mathrm{Ns}}{\mathrm{\Σ}}}{\mathrm{\ω}}_k^i{x}_k^i$

$x_k^i=f(x_{k-1}^i,x_{k-2}^i,...,x_{k-m}^i).$

Beneficial effect

Present invention achieves the method for a kind of mechanical breakdown indication, high dimensional feature can be carried out dimensionality reduction by the method, m rank HMM (HMM model) physical training condition evolutionary model is applied in low dimensional feature space, and in conjunction with particle filter formation order particles filter, non-gaussian, nonlinear state are indicated, comparing with existing predicting method, it is high to have predicted precision, feature applied widely.

Why the present invention has its reason of above-mentioned beneficial effect is: with local linear embedding grammar (LLE) dimensionality reduction, be not merely pure dimensionality reduction, one feature extraction mode effectively especially.It is not form new characteristic vector by selected part component from high dimensional feature vector, but high dimensional feature vector is in a mapping of lower dimensional space.And by such mappings, the characteristic vector of new low-dimensional can be obtained, it can keep topological structure and the locally relevance of legacy data, and the dimension of new vector can enough must be low, greatly reduce amount of calculation.Utilize least square method supporting vector machine (LS-SVR) to set up state equation, relative to neural net and SVR, there is stronger generalization ability and less computation complexity, be preferablyly applied in long prediction.In fault indication, the state evolution of fault not only has relation with previous moment, and and former moment have relation, application m rank HMM model is used for the evolutionary model of description state and can better reflects the truth of phylogeny.LS-SVR is dissolved in particle filter algorithm framework, if m=3, forms 3 rank particle filters.Compare with existing method, it indicates to have higher predicted precision to non-gaussian, nonlinear state, and the scope of application is wider.

Accompanying drawing explanation

Fig. 1 is method flow diagram of the present invention.

Embodiment

Below in conjunction with drawings and Examples, the present invention is further described.

Figure 1 shows the present invention realize the main-process stream of the mechanical breakdown indication of the order particles filter of LS-SVR modeling, this general flow chart contains each key step realized needed for final indication.The object of the invention is to develop to the signal characteristic of malfunction by extracting from normal condition, setting up state equation by the mode of data-driven, and prediction real-time online can be realized.Fallout predictor is realized by program, for the mechanical oscillation signal provided, extract the feature of corresponding signal, such as: 6 kinds of statistical properties of time domain are also applied LLE method and carried out low-dimensional feature extraction to it, then LS-SVR is utilized to build m rank HMM model, and with low dimensional feature, it is trained, the most enough being dissolved in particle filter algorithm framework forms m rank particle filter, finally indicates mechanical breakdown.

Here is concrete performing step

Step 1: to one group of signal developed to malfunction from normal condition by time slip-window intercept signal, each segment signal intercepted is as a sample, calculate 6 kinds of Time-domain Statistics features of each sample, Time-domain Statistics feature comprises peak value, average amplitude, root-mean-square amplitude, degree of skewness index and kurtosis index; Then set up sample characteristics space, and with local linear embedding grammar, low-dimensional feature extraction is carried out to the sample characteristics space built, obtain low-dimensional feature set.

Local linear embedding grammar (LLE) is the known method of this area, provides the process of concrete LLE feature extraction below:

Original feature space is made up of 6 kinds of statistical nature vectors, has peak value, average amplitude, root-mean-square amplitude, degree of skewness index and kurtosis index respectively.The d dimensional feature embedded in 6 dimensional feature space is extracted with LLE algorithm.

LLE algorithm can be summed up as three steps: (1) finds k Neighbor Points of each sample; (2) the partial reconstruction weight matrix of this sample point is calculated by the Neighbor Points of each sample point; (3) output valve of this sample is calculated by the partial reconstruction weight matrix of this sample point and its Neighbor Points.

The first step of algorithm is k the Neighbor Points calculating each sample point.K the sample point nearest relative to required sample point is defined as k Neighbor Points of required sample point.K is a value given in advance.

The second step of LLE algorithm is the partial reconstruction weight matrix (W) calculating sample point.Here an error function is defined, as follows:

$\min \mathrm{\ϵ}\left(W\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{|x_i-\underset{j=1}{\overset{k}{\mathrm{\Σ}}}{w}_j^i{x}_{\mathrm{ij}}|}^2---\left(1\right)$

Wherein N represents number of samples, x _ij(j=1,2 ..., k) be x _ik Neighbor Points, x _iwith x _ijbetween weights.And will satisfy condition: here ask for W matrix, need structure one local covariance matrix Q ⁱ.

$Q_{\mathrm{jm}}^i=(x_i-x_{\mathrm{ij}})(x_i-x_{\mathrm{im}})---\left(2\right)$

By above formula with combine, and adopt method of Lagrange multipliers, suboptimization can be obtained and rebuild weight matrix:

$w_j^i=\frac{\underset{m=1}{\overset{k}{\mathrm{\Σ}}}{\left(Q^i\right)}_{\mathrm{jm}}^{-1}}{\underset{p=1}{\overset{k}{\mathrm{\Σ}}}\underset{q=1}{\overset{k}{\mathrm{\Σ}}}{\left(Q^i\right)}_{\mathrm{pq}}^{-1}}---\left(3\right)$

In actual operation, Q ⁱmay be a singular matrix, now must regularization Q ⁱ, as follows:

Q ⁱ=Q ⁱ+ rI (4) wherein r is regularization parameter, and I is the unit matrix of a k × k.

The final step of LLE algorithm is mapped in lower dimensional space by all sample points.Mapping condition meets as follows:

$\min \mathrm{\ϵ}\left(Y\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{|y_i-\underset{j=1}{\overset{k}{\mathrm{\Σ}}}{w}_j^i{y}_{\mathrm{ij}}|}^2---\left(5\right)$

Wherein, ε (Y) is loss function value, y _ix _ioutput vector, y _ij(j=1,2 ..., k) be y _ik Neighbor Points, and two conditions to be met, that is:

$\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{y}_i=0,\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{y}_i{y}_i^T=I---\left(6\right)$

Wherein I is the unit matrix of m × m.Here can be stored in the sparse matrix W of N × N, work as x _jx _ineighbor Points time, otherwise, W _i,j=0.Then loss function can be rewritten as:

$\min \mathrm{\ϵ}\left(Y\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{M}_{i,j}{y}_i^T{y}_j---\left(7\right)$

Wherein M is the symmetrical matrix of a N × N, and its expression formula is:

M=(I-W) ^t(I-W) (8) will make loss function value reach minimum, then get the characteristic vector corresponding to minimum m nonzero eigenvalue that Y is M.In processing procedure, the characteristic value of M arranged from small to large, first characteristic value, almost close to zero, so casts out first characteristic value.Usually characteristic vector corresponding to the characteristic value between 2nd ~ m+1 is got as Output rusults.

In the present invention, being embedded as dimension is that d=1, k value is chosen and rule of thumb obtained as k=13.

Step 2: using the low-dimensional feature set that obtains in step 1 as sample training and test fallout predictor, odd point wherein in low-dimensional feature set is as training set, even number point is as test set, and training set is used to obtain forecast model, and test set is used to the effect of the prediction testing forecast model.

Step 3: adopt training set training least square method supporting vector machine, obtain least square method supporting vector machine forecast model, least square method supporting vector machine forecast model is added that plant noise is combined into high-order HMM:

$x_k={\hat{x}}_k+v_{k+1}$

${\hat{x}}_{k}f(x_{k-1},x_{k-2},......,x_{k-m})$

Wherein f () is least square method supporting vector machine forecast model, and m is exponent number, x _kfor the state value in the k moment, for the state estimation in the k moment, x _k-1, x _k-2..., x _k-mfor 1 moment before the k moment is to the state value in front m moment, v _k-1for the plant noise in k-1 moment.

The training study providing concrete LS-SVR below sets up state equation process:

Classical SVMs is finally solution convex quadratic programming problem, and its computational methods are relatively complicated, have higher amount of calculation.And LS-SVR used here adds least square method on the basis of classical SVR, quadratic programming problem is converted to the problem of a system of linear equations, simplifying the process that it solves, can matrix be upgraded simultaneously, making it be applicable to dynamic prediction when ensureing precision of prediction.

A young waiter in a wineshop or an inn takes advantage of SVMs (LS-SVR) to be that the people such as Belgian mathematician Suykens propose to improve the one of classical SVM algorithm.He changes the inequality constraints condition in the SVR algorithm of standard into equality constraint, with the loss function of the quadratic term of error as training, a convex quadratic programming problem is converted into solution quadratic linear equation group.As long as get for linear we

LSSVM regression problem can be expressed as:

$\min J(\mathrm{\ω},\mathrm{\ξ})=\frac{1}{2}{\left|\right|\mathrm{\ω}\left|\right|}^2+\frac{1}{2}\mathrm{\γ}\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{\mathrm{\ξ}}_i^2---\left(10\right)$

Wherein ω is weight vector, for slack variable, γ is the regularization parameter for balancing error of fitting and model complexity, and b is biased.

For optimization problem, introduce Lagrangian, finally obtain following equation

A in formula _i(i=1 ..., l) be Lagrange multiplier.By optimal condition, respectively with regard to ω, b, ξ, a, ask partial derivative, and be set to 0, can obtain:

$\frac{\∂L}{\∂b}=0\⇒\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{\mathrm{\α}}_i=0---\left(12\right)$

$\frac{\∂L}{\∂e_k}=0\⇒{\mathrm{\α}}_k={\mathrm{Ce}}_k,k=1,...,l$

In order to ask optimum α and b, can be obtained by KKT condition

$\left[\begin{array}{cc}0& I^T\\ I& \mathrm{\Ω}+{\mathrm{\γ}}^{-1}I\end{array}\right]\left[\begin{array}{c}b\\ \mathrm{\α}\end{array}\right]=\left[\begin{array}{c}0\\ y\end{array}\right]---\left(13\right)$

Wherein, y=[y ₁..., y _l]; α=[α ₁..., α _l] ^l; I=[1 ..., 1] _l; Ω is nuclear matrix, and the element that the i-th row j arranges is finally by system of linear equations (5) solution α and b out, least square regression function can be shown out:

$f\left(x\right)=\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{\mathrm{\α}}_{i}K(x,x_i)+b---\left(14\right)$

Wherein K (x, y) is kernel function,

Algorithm for Training process is as follows:

(1) input N=n training sample point, n is input amendment number, to N number of sample training least square method supporting vector machine;

(2) right | α _i| (i=1 ... .., N) sort, beta pruning M | α _i| minimum training sample (in general, M=N × 5%);

(3) establish N=N-M, remaining N number of training sample is formed new training sample;

(4) sample set after cutting down with least square method supporting vector machine training;

(5) then, get back to second step, iterate, till operation result is to the last deteriorated.

Set up LS-SVR forecast model

1, definite kernel function

Polynomial kernel function: K (x, y)=((xy)+c) ^d, wherein c>=0, d is any positive integer.

Gaussian radial basis function kernel function: $K(x_i,y)=\exp \left\{\frac{-{\left|\right|y-x_i\left|\right|}^2}{2{\mathrm{\σ}}^2}\right\}i=\mathrm{1,2},...,l.$

Fourier kernel function: wherein, x, y ∈ R, 0 < q < 1.

As long as any function meets Mercer condition, the kernel function of SVMs can be used as, adopt different functions to be kernel function, the Learning machine of non-linear decision surface that can be dissimilar in the constitution realization input space.The most frequently used kernel function is RBF gaussian radial basis function kernel function.

2, the variance of Modling model, solves parameter alpha and b, determines support vector, and beta pruning.

3, matched curve: according to the support vector obtained, simulate trend curve, and draw last conclusion by error analysis, namely solve support vector by least square method, then the support vector that " beta pruning " part is unnecessary, goes to fit to curve with remaining support vector and goes when meeting certain error to give a forecast.

High-order Hidden Markov (HMM) model:

Describe a mechanical vibrating system with certain suitable model, to analysis, to study this system extremely important.In practical engineering application, often employing dynamic space model describes problem wherein.Dynamic space model is a very important statistical and analytical tool, the Gauss-Markov linear model adopted as Kalman filter is exactly a good example, it describes with system equation the process that state develops in time, and measures equation to describe the noise variance relevant with state.Same, as long as write Gauss-Markov linear model as general mathe-matical map, just with two equations, more generally dynamical system can be described.

Mechanical breakdown indication problem becomes problem when being, analyzes common dynamic state-space model and is divided into state transition model p (x _k| x _k-1) and measurement model p (z _k| x _k), wherein, representative system in the state variable (implicit variable or parameter) of time k, for system is at the measured value of time k.To non-linear, nongausian process, its model can be expressed as:

$\begin{array}{c}{x}_k=f(x_{k-1},v_{k-1})\\ z_k=h(x_k,n_k)\end{array}---\left(15\right)$

In formula, with be respectively process noise and measurement noise, and be that separate, covariance is respectively Q _kand R _kzero-mean additive noise sequence.

In fault indication, the state evolution of fault not only has relation with previous moment, and and former moment have relation, application m rank HMM model is used for the evolutionary model of description state and can better reflects the truth of phylogeny.Then need m rank Markov model.M rank Markov model is as follows:

$\begin{array}{c}{x}_k={\hat{x}}_k+v_{k-1}\\ {\hat{x}}_k=f(x_{k-1},x_{k-2},......,x_{k-m})\end{array}---\left(16\right)$

X _kfor the state value in the k moment, for the state estimation in the k moment, x _k-1, x _k-2..., x _k-mfor 1 moment before the k moment is to the state value in front m moment.

Step 4: produce particle collection according to the high-order HMM that step 3 obtains the particle number of particle set is N _s, the i-th particle represent the state set in 0 to k moment, the corresponding weights of each particle of particle set are weights by obtain, wherein by k moment dbjective state posterior probability density obtain, z _krepresent the measuring value in k moment; And by weights be normalized to $\underset{i}{\mathrm{\&Sigma;}}{\mathrm{\&omega;}}_k^i=1.$

Step 5: judge sample degeneracy degree according to the weights that step 4 obtains, impose a condition if sample degeneracy arrives, then carry out particle resampling, and particle weights after calculating resampling.

Step 6: obtain state estimation by particle weights and particle value:

${\hat{x}}_k=\underset{i=1}{\overset{\mathrm{Ns}}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_k^i{x}_k^i$

Again state estimation is substituted into forecast model obtain state value.

Step 7: carry out subsequent time computing to step 6 according to step 4, and circulation is carried out until mechanical breakdown has indicated.

Provide the concrete derivation based on the order particles filtering indication algorithm of LS-SVR modeling in the present embodiment below:

The basic thought of particle filter is: first the empirical condition of foundation system mode vector is distributed in state space and produces one group of random sample set, these samples are claimed to be particle, then according to measuring and the weight of continuous adjustment particle and position, by the empirical condition distribution that the particle information correction after adjustment is initial.Its essence is and utilize the Discrete Stochastic of particle and weight composition thereof to estimate approximate relevant probability distribution, and estimate according to algorithm recursion renewal Discrete Stochastic.When sample size is very large, this Monte Carlo describes and is just similar to the real posterior probability density function of state variable.This technology is applicable to the nonlinear stochastic system of any non-gaussian background can stated with state-space model, and precision can approach optimal estimation, is a kind of effectively nonlinear filtering technique.

The state priori conditions Probability p (x of known dynamical system ₀), utilize k moment dbjective state x is described _kposterior probability distribution being corresponding weights is particle collection, δ () is sampling function, the i-th particle represent the state set in 0 to k moment.The Posterior probability distribution of k moment dbjective state can be weighted to discretely:

$p\left(x_{0:k}\right|z_{1:k})\&ap;\underset{i=1}{\overset{N_s}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_k^i\mathrm{\&delta;}(x_{0:k}-x_{0:k}^i)---\left(17\right)$

Wherein, weights are determined by importance sampling:

${\mathrm{\&omega;}}_k^i\&Proportional;\frac{p\left(x_{0:k}^i\right|z_{1:k})}{q\left(x_{0:k}^i\right|z_{1:k})}---\left(18\right)$

Corresponding q (x _0:k| z _1:k) be important density function, importance density function decomposition is:

$\begin{array}{c}q\left(x_{0:k}\right|z_{1:k})=q\left(x_k\right|x_{0:k-1},z_{1:k})q\left(x_{0:k-1}\right|z_{1:k-1})\\ =q\left(x_k\right|x_{k-m:k-1},z_{1:k})q\left(x_{k-m:k-1}\right|z_{1:k-1})\end{array}---\left(19\right)$

Then by obtaining particle by (19) with by q (x _0:k-1| z _1:k-1) the particle collection that obtains new particle collection can be obtained

Because posterior probability density function can be expressed as:

$\begin{array}{c}p\left(x_{0:k}\right|z_{1:k})=\frac{p\left(z_k\right|x_{0:k},z_{1:k-1})p\left(x_{0:k}\right|z_{1:k-1})}{p\left(z_k\right|z_{1:k-1})}\\ =\frac{p\left(z_k\right|x_k)p\left(x_k\right|x_{k-m:k-1})}{p\left(z_k\right|z_{1:k-1})}p\left(x_{k-m:k-1}\right|z_{1:k-1})\\ \&Proportional;p\left(z_k\right|x_k)p\left(x_k\right|x_{k-m:k-1})p\left(x_{k-m:k-1}\right|z_{1:k-1})\end{array}---\left(20\right)$

Then can obtain importance weight more new formula be:

${\mathrm{\&omega;}}_k^i\&Proportional;\frac{p\left(z_k\right|x_k^i)p\left(x_k^i\right|x_{k-m:k-1}^i)p\left(x_{k-m:k-1}^i\right|z_{1:k-1})}{q\left(x_k^i\right|x_{k-m:k-1}^i,z_{1:k})q\left(x_{k-m:k-1}^i\right|z_{1:k-1})}={\mathrm{\&omega;}}_{k-1}^i\frac{p\left(z_k\right|x_k^i)p\left(x_k^i\right|x_{k-m:k-1}^i)}{q\left(x_k^i\right|x_{k-m:k-1}^i,z_{1:k})}---\left(21\right)$

In order particles filtering algorithm, select the priori probability density being easy to most realize as importance density function, namely

$q\left(x_k^i\right|x_{k-m:k-1}^i,z_{1:k})=p\left(x_k^i\right|x_{k-m:k-1}^i)---\left(22\right)$

Formula (22) is substituted into (21), importance weight can be reduced to:

${\mathrm{\&omega;}}_k^i\&Proportional;{\mathrm{\&omega;}}_{k-1}^{i}p\left(z_k\right|x_k^i)---\left(23\right)$

By weights normalization, namely

${\mathrm{\&omega;}}_k^i={\mathrm{\&omega;}}_k^i/\underset{i=1}{\overset{N_s}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_k^i---\left(24\right)$

The maximum problem of particle filter is particle degeneracy phenomenon, namely along with the increase of time, importance weight likely focuses on minority particle, only can not go out posterior probability density function by effective expression, in order to avoid this degradation phenomena GORDON proposes method for resampling by these examples.In resampling, first calculate then N is judged _eff< N _th, N _thfor setting threshold, if N _eff< N _th, then particle will by resampling, the particle that weights are large is replicated, little being rejected of weights.Particle is selected the number of times copied to be directly proportional to its weights size.Particle after resampling has equal weights 1/N _s.The particle obtained after this step represents the Posterior distrbutionp of k moment dbjective state.

And posteriority probability density p (x _k| z _1:k) can be expressed as:

$p\left(x_k\right|z_{1:k})\&ap;\underset{i=1}{\overset{N_s}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_k^i\mathrm{\&delta;}(x_k-x_k^i)---\left(25\right)$

Work as N _sduring → ∞, there is law of great number can ensure the true posterior probability p (x of above formula programmable single-chip system _k| z _1:k).

The present embodiment is incorporated into LS-SVR forecast model in order particles filter frame, and when the exponent number getting particle filter is m=3, the 3 rank HMM describing fault growth tendency can be expressed as:

$\begin{array}{c}{x}_k={\hat{x}}_k+v_{k-1}\\ {\hat{x}}_k=f(x_{k-1},x_{k-2},x_{k-3})\end{array}---\left(26\right)$

Wherein, f (x _k-1, x _k-2, x _k-3) be the forecast model set up by LS-SVR, use current time state value, and the state value in the first two moment removes prediction subsequent time state value.If prediction p walks state, then recursive operation is carried out to formula (26), then can obtain the status predication value of p step.In order particles filtering, the population of each state is 100, and the initial weight of each particle corresponding is 0.01.

Derive according to above-mentioned theory, the present embodiment can be summarized as follows:

1, then one group of signal developed to malfunction from normal condition is set up to the feature space of sample by time slip-window intercept signal, and carry out low-dimensional feature extraction with the sample characteristics space of local linear embedding grammar (LLE) to the higher-dimension built, obtain low-dimensional feature set.

2, mechanical breakdown evolutionary model is set up with existing status data training LS-SVR.

3, particle collection is produced with the failure evolution modular form (26) adding process noise trained and 3 rank particle filters.The weights of each particle are calculated according to formula (23) and current measured value, then judge that sample degeneracy degree is the need of resampling, finally can calculate probability density function and state estimation respectively according to formula (25) and (27).According to particle value and current weights, a step status predication can calculate with following formula:

${\hat{x}}_k=\underset{i=1}{\overset{\mathrm{Ns}}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_k^i{x}_k^i---\left(27\right)$

Multistep status predication can be obtained by the continuous recurrence calculation of following formula:

${\hat{x}}_k=\underset{i=1}{\overset{\mathrm{Ns}}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_k^i{x}_k^i$

$x_k^i=f(x_{k-1}^i,x_{k-2}^i,...,x_{k-m}^i).$

4,3 are repeated until mechanical breakdown has indicated.

Evaluation index:

In order to evaluate estimated performance of the present invention and the result comparing emulation experiment, we use following index to illustrate.

Root-mean-square error (RMSE):

The total quantity at M index strong point, y _iwith represent i actual value and predicted value respectively, the value of RMSE is less, shows that precision of prediction is higher.

Simulated conditions: the rolling bearing data in emulation experiment are from the database of NASA (NASA).The data of this rolling bearing be four bearings by the test value of normal condition to running in continuous 35 days of malfunction, wherein within last 3 days, be considered to malfunction.Invariablenes turning speed is 2000 turns, and axle and bearing addition of the radial load of 6000 pounds.The vertical direction of each bearing and horizontal direction have respectively installed a vibration acceleration sensor, and sample frequency is 20kHz.In this experiment, the checking of fault indication algorithm is made of the inner ring fault data of the 3rd bearing in these data.

Simulation result: as shown in table 1, table 2.P in this table refers to that forward prediction p walks, ' RNN ' is the root-mean-square error that recurrent neural network p walks forward prediction, ' LS-SVR ' is that the p of least square method supporting vector machine walks the root-mean-square error of forward prediction, and ' LS-SVR 1 rank particle filter ' is the root-mean-square error that p that least square method supporting vector machine is combined with single order particle filter walks forward prediction.' the present invention ' is the root-mean-square error of the p step forward prediction of method of the present invention.

Early stage (first 32 days) the Condition forecasting experimental data table of table 1. bearing fault

Table 2. bearing fault state (latter 3 days) indication experimental data table

The root-mean-square error of the p step forward prediction that recurrent neural network, Least square support vector regression and the present invention obtain is listed, the root-mean-square error overstriking of same experiment condition Xia higher primary school in table 1 ~ table 2.From experimental data, the present invention adopt the order particles filtering method based on LS-SVR modeling to carry out mechanical breakdown indication than using traditional predicting method there is higher precision.Particularly when 5 step forward direction indication of bearing fault state, root-mean-square error reduces 0.102.

Claims (2)

1. based on a fault predicting method for the order particles filter of LS-SVR modeling, it is characterized in that: comprise following

Step:

Step 1: to one group of signal developed to malfunction from normal condition by time slip-window intercept signal, each segment signal intercepted is as a sample, calculate 6 kinds of Time-domain Statistics features of each sample, Time-domain Statistics feature comprises peak value, average amplitude, root-mean-square amplitude, degree of skewness index and kurtosis index; Then set up sample characteristics space, and with local linear embedding grammar, low-dimensional feature extraction is carried out to the sample characteristics space built, obtain low-dimensional feature set;

Step 2: using the low-dimensional feature set that obtains in step 1 as sample training and test fallout predictor, odd point wherein in low-dimensional feature set is as training set, even number point is as test set, and training set is used to obtain forecast model, and test set is used to the effect of the prediction testing forecast model;

Step 3: adopt training set training least square method supporting vector machine, obtain least square method supporting vector machine forecast model, least square method supporting vector machine forecast model is added that plant noise is combined into high-order HMM:

$x_k={\hat{x}}_k+v_{k-1}$

${\hat{x}}_k=f(x_{k-1},x_{k-2},......,x_{k-m})$

Wherein f () is least square method supporting vector machine forecast model, and m is exponent number, x _kfor the state value in the k moment, for the state estimation in the k moment, x _k-1, x _k-2..., x _k-mfor 1 moment before the k moment is to the state value in front m moment, v _k-1for the plant noise in k-1 moment;

Step 4: produce particle collection according to the high-order HMM that step 3 obtains the particle number of particle set is N _s, the i-th particle represent the state set in 0 to k moment, the corresponding weights of each particle of particle set are weights by obtain, wherein by k moment dbjective state posterior probability density obtain, z _krepresent the measuring value in k moment; And by weights be normalized to $\underset{i}{\mathrm{\&Sigma;}}{\mathrm{\&omega;}}_k^i=1;$

Step 5: judge sample degeneracy degree according to the weights that step 4 obtains, impose a condition if sample degeneracy arrives, then carry out particle resampling, and particle weights after calculating resampling;

Step 6: obtain state estimation by particle weights and particle value:

${\hat{x}}_k=\underset{i=1}{\overset{\mathrm{Ns}}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_k^i{x}_k^i$

Again state estimation is substituted into forecast model obtain state value;

Step 7: carry out subsequent time computing to step 6 according to step 4, and circulation is carried out until mechanical breakdown has indicated.

2. the fault predicting method of a kind of order particles filter based on LS-SVR modeling according to claim 1, is characterized in that: in step 6, obtains state estimation by particle weights and particle value by following multistep status predication:

${\hat{x}}_k=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_k^i{x}_k^i$

$x_k^i=f(x_{k-1}^i,x_{k-2}^i,...,x_{k-m}^i).$