CN101859146B - Satellite fault prediction method based on predictive filtering and empirical mode decomposition - Google Patents

Abstract

The invention relates to a satellite fault prediction method based on predictive filtering and empirical mode decomposition, which belongs to the method for safe operation and fault prediction of the space satellite, and solves the problems of serious noise influence and inaccurate prediction of the fault trend in the traditional satellite fault prediction and diagnosis method. The method comprises the following specific steps: 1, estimating errors of the satellite control system through predictive filtering according to the dynamic relationship of the nonlinear attitude of the satellite, thereby obtaining a system model error; 2, carrying out empirical mode decomposition for the system model error obtained in Step 1 to obtain the IMF component and the residual component of the n-order intrinsic mode function; and 3, establishing a fault trend model related to the residual component obtained in Step 2 through time series analysis, thereby completing the prediction and detection of minor and slow-variation faults. The invention can be applied in the field of fault diagnosis of the satellite attitude control system.

Description

A kind of satellite failure Forecasting Methodology based on predictive filtering and empirical modal decomposition

Technical field

The present invention relates to the safe operation and the failure prediction method of Aerospace Satellite, be specifically related to a kind of satellite attitude control system failure prediction method based on predictive filtering and empirical modal decomposition.

Background technology

Attitude of satellite control is to obtain and keep satellite in aspect-stabilized method and process.Present high precision three-axis attitude stabilization satellite during operate as normal, generally adopts the main topworks of momenttum wheel as attitude control system in orbit, realizes control to the attitude of satellite by the momentum-exchange between momenttum wheel and the satellite.The reliability of momenttum wheel will directly influence the feasibility and the reliability of the attitude control of whole satellite.The current method that mainly is based on hardware redundancy for the fault diagnosis and the reconstruct of satellite attitude control system, this satellite posture control system for complexity is far from being enough.

Safe operation and failure prediction to satellite adopts the predictive filtering method at present, the predictive filtering method is a kind of method of estimation that is applicable to the nonlinear system with unknown input or model error, the non-linear Proctor Central that the viewpoint that its thought source is controlled in Lu from system proposes, on this basis, Crassidis and Markley have proposed a kind of new Real-Time Filtering algorithm according to the least model error criterion---and predictive filtering (Predictive Filtering, PF).Predictive filtering is exported the model error of real-time estimating system by comparing and measuring output and prediction, thereby the correction wave filter state is realized the estimation to time of day.Because predictive filter has the ability of estimation model sum of errors system state simultaneously, domestic Li Ji, a big vast battle-axe used in ancient China have been incorporated into fault diagnosis field by fault being considered as a kind of special model error with the predictive filtering method.The noise that exists in the estimated result of predictive filtering can influence diagnosis performance, can adopt the low-pass filter method to suppress high frequency noise, but because characteristic the unknown of noise, the method that adopts the low-pass filter method to suppress noise often lost efficacy, cause and to make accurate prediction to fault, be unfavorable for further research ensuing fault detection and diagnosis.If can accurately predict, before taking place, fault, can improve the feasibility and the reliability of the attitude control of whole satellite to the estimation of time of day to satellite failure.

Summary of the invention

There is serious, the problem that can't accurately predict fault trend affected by noise in the present invention in order to solve the conventional satellite failure prediction method, and a kind of satellite failure Forecasting Methodology of decomposing based on predictive filtering and empirical modal is provided.

Detailed process of the present invention is as follows:

Step 1: utilize satellite nonlinear attitude kinetics relation, adopt the method for predictive filtering that the satellite control system error is estimated, obtain system model error term;

Step 2: the system model error term that step 1 obtains is carried out the empirical modal decomposition, E rank eigenmode state function IMF component and residual component before obtaining;

Step 3: utilize Time series analysis method to set up the model of the fault trend of the residual component that obtains about step 2, finish the forecast and the detection of small gentle change fault.

The present invention is according to the nonlinear attitude kinetics relation of satellite, fault amount and model uncertainty sum are considered as model error, compare with low-pass filter disposal route in the past, can eliminate The noise effectively, utilize predictive filtering method estimation model error, the model error estimated value that predictive filtering obtains is carried out the empirical modal decomposition, obtain its several in solid model attitude component and trend component, improved precision of prediction, utilize Time series analysis method that the trend of fault is predicted, can be to gradual early stage, small fault forecasts accurately and detects that method is concisely effective.The fault diagnosis field that is used for satellite attitude control system.

Description of drawings

Fig. 1 is the process flow diagram based on the satellite failure Forecasting Methodology of predictive filtering and empirical modal decomposition; Fig. 2 is embodiment three process flow diagrams; Predictive filtering estimated result when Fig. 3 for pitch axis gradual fault takes place; Fig. 4 is the empirical modal decomposition result of yaw axis (non-fault); Fig. 5 is the empirical modal decomposition result of pitch axis (gradual fault); Fig. 6 is the empirical modal decomposition result of wobble shaft (non-fault); Predicting the outcome of gradual fault takes place for pitch axis in Fig. 7; Fig. 8 predicts the outcome for yaw axis generation micromutation fault.

Embodiment

Embodiment one: a kind of satellite failure Forecasting Methodology based on predictive filtering and empirical modal decomposition, detailed process is as follows:

Step 1: utilize satellite nonlinear attitude kinetics relation, adopt the method for predictive filtering that the satellite control system error is estimated, obtain system model error term;

Step 2: the system model error term that step 1 obtains is carried out the empirical modal decomposition, E rank eigenmode state function IMF component and residual component before obtaining;

Step 3: utilize Time series analysis method to set up the model of the fault trend of the residual component that obtains about step 2, finish the forecast and the detection of small gentle change fault.

Empirical mode decomposition method (Empirical Mode Decomposition, EMD) be a kind of new signal decomposition method that proposed by scholar Huang E in 1998, can utilize the variation of signal internal time yardstick to do the parsing of energy and frequency, with signal be launched into solid model state function in several (Intrinsic Mode Function, IMF).Be different from and use the classic method of solid form window for the boundary basis function, the basis function of EMD extracts from signal and obtains, and promptly uses IMF to do substrate.And IMF must satisfy following condition:

1) in whole function, the number of extreme point equates with the number that passes through zero point or differs 1;

2) be zero by the defined envelope local mean value of local extremum envelope at any time.

Wherein, in first condition and the traditional gaussian stationary process narrow frequency range require similar.Second condition is a new idea: globality is required to change into the locality requirement, make instantaneous frequency can not cause unnecessary rocking because of the existence of asymmetric waveform.The EMD and the HHT that rely on these two conditions to make up are considered to handle forcefully adaptive approach non-linear, non-stationary signal, be in recent years to based on the linearity of Fourier transform and the important breakthrough of stable state analysis of spectrum, and obtained using widely.Utilize the empirical mode decomposition method can be adaptive with the composition of signal decomposition for different instantaneous frequencys, thus the noise contribution adaptively in the erasure signal.

The forecast model of setting up fault trend is the main contents of fault forecast.Autoregressive model (Autoregressive model, AR model) is a kind of common model in the time series analysis.Advantages such as it is simple that the AR model has modeling, and calculated amount is little.The above AR model of single order is applicable to stochastic process stably, and single order AR model has m＞different characteristic of 1 model wherein with AR (m), and as the special case of AR model, single order AR model can be predicted nonstationary random process.

Basic thought of the present invention is major part or a kind of important component that fault is considered as the system model error, utilize the model error item of satellite nonlinear attitude kinetics relation design predictive filter estimating system, utilize empirical mode decomposition method that estimated result is handled then, so that can diagnose the small gentle change initial failure of satellite attitude control system.

The objective of the invention is to be achieved through the following technical solutions: the model error item that utilizes satellite nonlinear attitude kinetics relation design predictive filter estimating system, estimated result is carried out empirical modal to be decomposed, obtain some rank IMF and trend term, utilize Time series analysis method to set up the model of trend term, carry out the forecast and the detection of small gentle change fault.

The present invention compared with prior art has following advantage:

1) failure prediction method proposed by the invention utilizes empirical mode decomposition method, compares with the low-pass filter disposal route, can more effectively eliminate The noise, extracts the tendency information of fault amount.

2) fault failure prediction method proposed by the invention is utilized empirical mode decomposition method, and signal is divided into the residual component of several IMF and non-stationary, sets up the AR model on this basis, can improve precision of prediction.

3) the present invention forecasts the fault amount and detects, and directly compares based on the diagnostic method of threshold value, can improve the responsive ability for gradual early stage, small fault.

Embodiment two, present embodiment are further specifying embodiment one, utilize satellite nonlinear attitude kinetics relation in the step 1, adopt the method for predictive filtering the satellite control system error to be estimated the process that obtains system model error term is:

Set the model error of satellite control system and be made up of satellite executing mechanism fault and model uncertainty, the measurement equation of predictive filtering system and real system is respectively:

State equation:

$\stackrel{\·}{x}\left(t\right)=f\left[x\left(t\right)\right]+g\left[x\left(t\right)\right]d\left(t\right)$

${\tilde{y}}_k=c\left[x\left(t_k\right)\right]+v_k$

The predictive filtering equation is:

$\hat{x}\left(t\right)=f\left[\hat{x}\left(t\right)\right]+g\left[\hat{x}\left(t\right)\right]\hat{d}\left(t\right)$

$\hat{y}\left(t\right)=c\left[\hat{x}\left(t\right)\right]$

X (t) ∈ R wherein ^sBe state vector,

Be the estimated value of state vector, f ∈ R ^sBe function of state that can be little,

Be known model error distribution matrix,

Be the measurement functions vector,

Be the estimation of unknown error, the measurement of real system is output as discrete form,

Be shown in t _kMeasured value constantly, v _k∈ R ^mBe the measurement noise, and set v _kBe that average is zero, covariance matrix is Q ∈ R ^{M * m}White Gaussian noise;

Constantly will measure function at t+ Δ t and carry out Taylor expansion, obtain:

$\hat{y}(t+\mathrm{\Δt})=\hat{y}\left(t\right)+z[\hat{x}\left(t\right),\mathrm{\Δt}]+\mathrm{\Λ}\left(\mathrm{\Δt}\right)S\left[\hat{x}\left(t\right)\right]d\left(t\right)$

Δ t=t wherein _K+1-t _kIt is the sampling period; This sampling period is normal value, y (t+ Δ t)=y _K+1, matrix

I element be:

$z_i[\hat{x}\left(t\right),\mathrm{\Δt}]=\underset{k=1}{\overset{p_i}{\mathrm{\Σ}}}\frac{{\mathrm{\Δt}}^k}{k!}{L^k}_f\left(c_i\right)$

P wherein _iDifferential exponent number when occurring d (t) in the Taylor expansion first,

Be c _iK rank Lie derivative;

Matrix Λ (Δ t) is a diagonal matrix, and its diagonal element is:

${\mathrm{\λ}}_{\mathrm{ii}}=\frac{{\mathrm{\Δt}}^{p_i}}{p_i!},$ i＝1，2，…，m

Matrix $S\left[\hat{x}\left(t\right)\right]\∈R^{m\×q},$ The element that its i is capable is:

$s_i=\left\{L_{g_1}\right[L_f^{p_i-1}\left(c_i\right)],L,L_{g_q}[L_f^{p_i-1}\left(c_i\right)\left]\right\},$ i＝1，2，…，m

Get performance index function:

$J\left[d\left(t\right)\right]=\frac{1}{2}{\{\tilde{y}(t+\mathrm{\Δt})-\tilde{y}(t+\mathrm{\Δt})\}}^T{Q}^{-1}\{\tilde{y}(t+\mathrm{\Δt})-\tilde{y}(t+\mathrm{\Δt})\}+\frac{1}{2}{d}^T\left(t\right)\mathrm{Wd}\left(t\right)$

W ∈ R wherein ^{Q * q}Be positive semidefinite weighting matrix, adopt the gradient optimizing algorithm that performance index are optimized, the estimated value that obtains the model error item is:

$\hat{d}\left(t\right)=-{\{{\left[\mathrm{\Λ}\left(\mathrm{\Δt}\right)S\left(\hat{x}\right)\right]}^T{Q}^{-1}\mathrm{\Λ}\left(\mathrm{\Δt}\right)S\left(\hat{x}\right)+W\}}^{-1}\left[\mathrm{\Λ}\left(\mathrm{\Δt}\right)S\left(\hat{x}\right)\right]Q^{-1}[z(\hat{x},\mathrm{\Δt})-\tilde{y}(t+\mathrm{\Δt})+\hat{y}\left(t\right)].$

Because the estimated value of model error item

Contain bigger noise component, can not be directly used in fault diagnosis, be necessary to carry out next step and handle.

Embodiment three, present embodiment are to the further specifying of concrete enforcement one, and the process of E rank eigenmode state function IMF component and residual component was before step 2 obtained:

The estimated value of initialization system model error item is Time t=1,2 ..., N,

Step a, IMF decomposable process initialization: n=1, and satisfy relational expression

Set up, wherein r _N-1(t) be remaining residual error function after (n-1) inferior decomposition;

Step b, screening process initialization: k=1, and satisfy relational expression h _{N (k-1)}(t)=r _N-1(t) set up, wherein h _{N (k-1)}(t) be through the survival function after the k-1 time screening during the n time IMF decomposes;

Step c, obtain the estimated value of system model error term according to screening sequence Survival function h after screening through the k time in the remaining residual error function of decomposing through the n time intrinsic mode function _Nk(t);

Steps d, the survival function h that adopts the judgement of standard deviation criterion to obtain _Nk(t) whether satisfy the condition of eigenmode state function, promptly

Whether less than threshold value T, 0.2≤T≤0.3;

Judged result is for being, execution in step e, and judged result is not for, and then k=k+1 returns execution in step c,

Step e, the n time intrinsic mode function IMF component c of acquisition _n(t)=h _Nk(t);

Step f, obtain the estimated value of system model error term

Remaining residual error function through the n time intrinsic mode function decomposition $r_n\left(t\right)=d\left(t\right)-{\mathrm{\Σ}}_{\mathrm{\α}=1}^n{c}_{\mathrm{\α}}\left(t\right);$

Step g, make n=n+1, return execution in step b, E rank eigenmode state function IMF component and residual component before satisfying n=E and obtaining.

Embodiment four, present embodiment are that step c obtains the estimated value of system model error term according to screening sequence to the further specifying of concrete enforcement three

Survival function h after screening through the k time in the remaining residual error function of decomposing through the n time intrinsic mode function _Nk(t) process is:

Step c1, the survival function h after utilizing cubic spline function to obtain in the residue trend function that system model error term d (t) decomposes through the n time intrinsic mode function to screen through the k-1 time _{N (k-1)}(t) upper and lower enveloping curve;

Step c2, the described survival function h of calculating _{N (k-1)}(t) upper and lower enveloping curve is at time t=1, and 2 ..., the average in the N

Step c3, obtain the estimated value of system model error term Survival function after screening through the k time in the residue trend function of decomposing through the n time intrinsic mode function $h_{\mathrm{nk}}\left(t\right)=h_{n(k-1)}\left(t\right)-{\stackrel{\‾}{m}}_{n(k-1)}\left(t\right).$

Embodiment five, present embodiment are to the further specifying of concrete enforcement three, and among the step f system model error term x (t) are carried out 4 intrinsic mode functions and decompose, and obtain 4 eigenmode state function component IMF1, IMF2, IMF3 and IMF4:{c ₁(t), c ₂(t), c ₃(t), c ₄(t) }; And acquisition is through the remaining residual error function r of 4 intrinsic mode function decomposition ₄(t).

Embodiment six, present embodiment are further specifying concrete enforcement one, utilize Time series analysis method to set up the model of the fault trend of the residual component that obtains about step 2 in the step 3, finish the forecast of small gentle change fault and the process of detection and be:

Residual component as the trend signal, is set up its autoregressive model, and expression formula is:

r(t)＝a ₁r(t-1)+a ₂r(t-2)+…+a _pr(t-p)+ε(t)

Wherein, r (t) is the time series of residual error function, and p is a model order, a _iBe model parameter, ε (t) is a white noise sequence;

Observed reading is r (0), r (1) ..., r (N), the order of forecast model are p, get following system of equations by the autoregressive model expression formula:

$\left\{\begin{array}{c}r\left(p\right)=a_{1}r(p-1)+a_{2}r(p-2)+...+a_{p}r\left(0\right)+\mathrm{\ϵ}\left(p\right)\\ r(p+1)=a_{1}r\left(p\right)+a_{2}r(p-1)+...+a_{p}r\left(1\right)+\mathrm{\ϵ}(p+1)\\ \·\\ \·\\ \·\\ r\left(N\right)=a_{1}r(N-1)+a_{2}r(N-2)+...+a_{p}r(N-p)+\mathrm{\ϵ}\left(N\right)\end{array}\right.$

Order:

$R=\left[\begin{array}{c}r\left(p\right)\\ r(p+1)\\ \·\\ \·\\ \·\\ r\left(N\right)\end{array}\right],$ $A=\left[\begin{array}{c}{a}_1\\ a_2\\ \·\\ \·\\ \·\\ a_p\end{array}\right],$

$B=\left[\begin{array}{cccc}r(p-1)& r(p-2)& \·\·\·& r\left(0\right)\\ r\left(p\right)& r(p-1)& \·\·\·& r\left(1\right)\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ r(N-1)& r(N-2)& \·\·\·& r(N-p)\end{array}\right],$ $\mathrm{\Δ}=\left[\begin{array}{c}\mathrm{\ϵ}\left(p\right)\\ \mathrm{\ϵ}(p+1)\\ \·\\ \·\\ \·\\ \mathrm{\ϵ}\left(N\right)\end{array}\right]$

Then the matrix form of following formula is:

R＝BA+Δ

The least square solution of model parameter is:

A＝(B ^TB) ^-1B ^TR

Obtain the model parameter A of autoregressive model, utilize model parameter A to carry out failure prediction and detection.

The above AR model of single order is applicable to stochastic process stably, and single order AR model has and the different characteristic of AR (m) (m＞1) model, and as the special case of AR model, single order AR model can be predicted nonstationary random process.

Set forth the specific embodiment of the present invention below by satellite executing mechanism Fault Diagnosis Simulation example:

Execution in step one: estimate topworks's fault amount with the predictive filtering method.

Dynamic (dynamical) state equation form of the attitude of satellite and discrete measurement equation are as follows:

$\left\{\begin{array}{c}{\stackrel{\·}{x}}_1\left(t\right)=\left(\frac{I_y-I_z}{I_x}\right)x_2{x}_3+\left(\frac{1}{I_x}\right)T_x+\left(\frac{1}{T_x}\right)f_x\left(t\right)\\ {\stackrel{\·}{x}}_2\left(t\right)=\left(\frac{I_x-I_z}{I_y}\right)x_1{x}_3+\left(\frac{I}{I_y}\right)T_y+\left(\frac{1}{T_y}\right)f_y\left(t\right)\\ {\stackrel{\·}{x}}_3\left(t\right)=\left(\frac{I_x-I_y}{I_z}\right)x_1{x}_2+\left(\frac{I}{I_z}\right)T_z+\left(\frac{1}{I_z}\right)f_z\left(t\right)\end{array}\right.$

$\tilde{y}\left(t_k\right)=x\left(t_k\right)+v\left(t_k\right)$

In the formula, state variable x (t) is the attitude angular velocity of satellite, I _x, I _y, I _zBe the main shaft inertia of satellite, f _x(t), f _y(t), f _z(t) be actuator momenttum wheel fault.According to the predictive filtering theory, obtain the Fault Estimation amount and be:

$f\left(t\right)=-{\{{\left[\mathrm{\Λ}\left(\mathrm{\Δt}\right)S\left(\hat{x}\right)\right]}^T{Q}^{-1}\mathrm{\Λ}\left(\mathrm{\Δt}\right)S\left(\hat{x}\right)+W\}}^{-1}\·$

$\left[\mathrm{\Λ}\left(\mathrm{\Δt}\right)S\left(\hat{x}\right)\right]Q^{-1}[z(\hat{x},\mathrm{\Δt})-\tilde{y}(t+\mathrm{\Δt})+\hat{y}\left(t\right)]$

Wherein

$z(\hat{x},\mathrm{\Δt})=\left[\begin{array}{c}\left(\frac{I_y-I_z}{I_x}\right){\hat{x}}_2{\hat{x}}_3+\left(\frac{1}{I_x}\right)T_x\\ \left(\frac{I_z-I_x}{I_y}\right){\hat{x}}_1{\hat{x}}_3+\left(\frac{I}{I_y}\right)T_y\\ \left(\frac{I_x-I_y}{I_z}\right){\hat{x}}_1{\hat{x}}_2+\left(\frac{I}{I_z}\right)T_z\end{array}\right]\mathrm{\Δt},$ $\mathrm{\Λ}\left(\mathrm{\Δt}\right)=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]\mathrm{\Δt},$

$S\left(\hat{x}\right)=\left[\begin{array}{ccc}\frac{1}{I_x}& 0& 0\\ 0& \frac{1}{{\mathrm{<figure-callout\; id="0"\; label="I\; y"\; filenames="100715164241.png,100715164245.png"\; state="\{\{state\}\}">I</figure-callout>}}_y}& 0\\ 0& 0& \frac{1}{I_z}\end{array}\right]$

With pitch axis generation slope when the t=19.6s is the gradual fault f of actuator of 0.004Nm/s _yBe example, the Fault Estimation result of each as shown in Figure 3, wherein, three curves of Fig. 3 respectively are x, y, three change in coordinate axis direction of z (driftage, pitching, wobble shaft) Fault Estimation value, horizontal ordinate express time, unit is second, and ordinate is represented the Fault Estimation value, and unit is a Newton meter.Owing to there is bigger noise contribution in the estimated result, can't be directly used in fault diagnosis, must handle to extract fault signature signal.

Execution in step two: the predictive filtering result is carried out empirical modal decompose, obtain several IMF and residual component.

For avoiding the fault signature of mass loss rates gyro, promptly carry out an empirical modal through 128 samplings and decompose.Empirical modal decomposes to the 4th layer can be finished.

When the t=19.6s pitch axis breaks down, the empirical modal decomposition result of yaw axis (non-fault) Fault Estimation amount is Fig. 4, Fig. 5 is the empirical modal decomposition result of the Fault Estimation amount of pitch axis, and Fig. 6 is the empirical modal decomposition result of wobble shaft (non-fault) Fault Estimation amount.Fig. 4, Fig. 5, Fig. 6 are respectively driftage, pitching and lift-over three axis signals of gathering are carried out the result that empirical modal decomposes, each figure is respectively the residual component after 1,2,3,4,5 rank IMF (interior solid model state function) component and the decomposition from top to bottom, the horizontal ordinate express time, unit is second, ordinate is represented IMF component or residual component, and unit is a Newton meter.

Execution in step three: set up the AR model of residual component, carry out failure prediction and detection.

With the forecast model of single order AR model, utilize aforementioned the least square estimation method to obtain the model parameter value of AR model as residual component.Predicting the outcome of driftage, pitching and wobble shaft as Fig. 7.Fig. 7 is for to the described example of Fig. 3 (pitch axis breaks down when the t=19.6s), the x that adopts the inventive method to obtain, y, three change in coordinate axis direction of z (driftage, pitching, wobble shaft) Fault Estimation value, the horizontal ordinate express time, unit is second, and ordinate is represented the Fault Estimation value, and unit is a Newton meter.Adopt suitable threshold (0.02Nm), can realize detecting in advance small, initial failure.

In addition, also can verify the validity of failure prediction method proposed by the invention by small mutation failure.When t=20s, the yaw axis fault of undergoing mutation, amplitude is 0.01Nm, the failure prediction result of then driftage, pitching, wobble shaft such as Fig. 8.Fig. 8 is when 20s, when yaw axis generation amplitude is the mutation failure of 0.01N.m, and the x that adopts the inventive method to obtain, y, three change in coordinate axis direction of z (driftage, pitching, wobble shaft) Fault Estimation value, horizontal ordinate express time, unit is second, and ordinate is represented the Fault Estimation value, and unit is a Newton meter.

Above-mentioned fault diagnosis result can be verified the validity of failure prediction method proposed by the invention.

Claims (6)

1. satellite failure Forecasting Methodology of decomposing based on predictive filtering and empirical modal is characterized in that detailed process is as follows:

Step 1: utilize satellite nonlinear attitude kinetics relation, adopt the method for predictive filtering that the satellite control system error is estimated, obtain system model error term;

Step 2: the system model error term that step 1 obtains is carried out the empirical modal decomposition, E rank eigenmode state function IMF component and residual component before obtaining;

Step 3: utilize Time series analysis method to set up the model of the fault trend of the residual component that obtains about step 2, finish the forecast and the detection of small gentle change fault.

2. a kind of satellite failure Forecasting Methodology according to claim 1 based on predictive filtering and empirical modal decomposition, it is characterized in that utilizing in the step 1 satellite nonlinear attitude kinetics relation, adopt the method for predictive filtering the satellite control system error to be estimated the process that obtains system model error term is:

Set the model error of satellite control system and be made up of satellite executing mechanism fault and model uncertainty, its nonlinear state equation and measurement equation are respectively:

State equation:

$\stackrel{\&CenterDot;}{x}\left(t\right)=f\left[x\left(t\right)\right]+g\left[x\left(t\right)\right]d\left(t\right)$

${\tilde{y}}_k=c\left[x\left(t_k\right)\right]+v_k$

The predictive filtering equation is:

$\hat{x}\left(t\right)=f\left[\hat{x}\left(t\right)\right]+g\left[\hat{x}\left(t\right)\right]\hat{d}\left(t\right)$

$\hat{y}\left(t\right)=c\left[\hat{x}\left(t\right)\right]$

X (t) ∈ R wherein ^sBe state vector,

Be the estimated value of state vector, f ∈ R ^sBe function of state that can be little,

Be known model error distribution matrix,

Be the measurement functions vector, d (t) ∈ R ^qBe unknown error,

Be its estimated value,

Be illustrated in t _kThe measurement output of real system constantly is discrete form, v _k∈ R ^mBe the measurement noise, and set v _kBe that average is zero, covariance matrix is Q ∈ R ^{M * m}White Gaussian noise;

Constantly will measure function at t+ Δ t and carry out Taylor expansion, obtain:

$\hat{y}(t+\mathrm{\&Delta;t})=\hat{y}\left(t\right)+z[\hat{x}\left(t\right),\mathrm{\&Delta;t}]+\mathrm{\&Lambda;}\left(\mathrm{\&Delta;t}\right)S\left[\hat{x}\left(t\right)\right]d\left(t\right)$

Δ t=t wherein _K+1-t _kIt is the sampling period; This sampling period is normal value, y (t+ Δ t)=y _K+1, matrix

I element be:

$z_i[\hat{x}\left(t\right),\mathrm{\&Delta;t}]=\underset{k=1}{\overset{p_i}{\mathrm{\&Sigma;}}}\frac{{\mathrm{\&Delta;t}}^k}{k!}{L^k}_f\left(c_i\right)$

P wherein _iDifferential exponent number when occurring d (t) in the Taylor expansion first,

Be c _iK rank Lie derivative;

Matrix A (Δ t) is a diagonal matrix, and its diagonal element is:

${\mathrm{\&lambda;}}_{\mathrm{ii}}=\frac{{\mathrm{\&Delta;t}}^{p_i}}{p_i!},$ i＝1，2，…，m

Matrix $S\left[\hat{x}\left(t\right)\right]\&Element;R^{m\&times;q},$ The element that its i is capable is:

$s_i=\left\{L_{g_1}\right[L_f^{p_i-1}\left(c_i\right)],L,L_{g_q}[L_f^{p_i-1}\left(c_i\right)\left]\right\},$ i＝1，2，…，m

Get performance index function:

$J\left[d\left(t\right)\right]=\frac{1}{2}{\{\tilde{y}(t+\mathrm{\&Delta;t})-\tilde{y}(t+\mathrm{\&Delta;t})\}}^T{Q}^{-1}\{\tilde{y}(t+\mathrm{\&Delta;t})-\tilde{y}(t+\mathrm{\&Delta;t})\}+\frac{1}{2}{d}^T\left(t\right)\mathrm{Wd}\left(t\right)$

W ∈ R wherein ^{Q * q}Be positive semidefinite weighting matrix, adopt the gradient optimizing algorithm that performance index are optimized, obtain being estimated as of model error item:

$\hat{d}\left(t\right)=-{\{{\left[\mathrm{\&Lambda;}\left(\mathrm{\&Delta;t}\right)S\left(\hat{x}\right)\right]}^T{Q}^{-1}\mathrm{\&Lambda;}\left(\mathrm{\&Delta;t}\right)S\left(\hat{x}\right)+W\}}^{-1}\&CenterDot;$

$\left[\mathrm{\&Lambda;}\left(\mathrm{\&Delta;t}\right)S\left(\hat{x}\right)\right]Q^{-1}[z(\hat{x},\mathrm{\&Delta;t})-\tilde{y}(t+\mathrm{\&Delta;t})+\hat{y}\left(t\right)].$

3. a kind of satellite failure Forecasting Methodology of decomposing according to claim 2 based on predictive filtering and empirical modal, it is characterized in that step 2 obtain before the process of E rank eigenmode state function IMF component and residual component be:

Initialization system model error item estimated value is

Time t=1,2 ..., N,

Step a, IMF decomposable process initialization: n=1, and satisfy relational expression

Set up, wherein r _N-1(t) be remaining residual error function after (n-1) inferior decomposition;

Step b, screening process initialization: k=1, and satisfy relational expression h _{N (k-1)}(t)=r _N-1(t) set up, wherein h _{N (k-1)}(t) be through the survival function after the k-1 time screening during the n time IMF decomposes;

Step c, obtain system model error term estimated value according to screening sequence

Survival function h after screening through the k time in the remaining residual error function of decomposing through the n time intrinsic mode function _Nk(t);

Steps d, the survival function h that adopts the judgement of standard deviation criterion to obtain _Nk(t) whether satisfy the condition of eigenmode state function, promptly

Whether less than threshold value T, 0.2≤T≤0.3;

Judged result is for being, execution in step e, and judged result is not for, and then k=k+1 returns execution in step c,

Step e, the n time intrinsic mode function IMF component c of acquisition _n(t)=h _Nk(t);

Step f, obtain system model error term estimated value Remaining residual error function through the n time intrinsic mode function decomposition $r_n\left(t\right)=\hat{d}\left(t\right)-{\mathrm{\&Sigma;}}_{\mathrm{\&alpha;}=1}^n{c}_{\mathrm{\&alpha;}}\left(t\right);$

Step g, make n=n+1, return execution in step b, E rank eigenmode state function IMF component and residual component before satisfying n=E and obtaining.

4. a kind of satellite failure Forecasting Methodology based on predictive filtering and empirical modal decomposition according to claim 3 is characterized in that step c obtains system model error term estimated value according to screening sequence

Survival function h after screening through the k time in the remaining residual error function of decomposing through the n time intrinsic mode function _Nk(t) process is:

Step c1, utilize cubic spline function to obtain system model error term estimated value Survival function h after screening through the k-1 time in the residue trend function of decomposing through the n time intrinsic mode function _{N (k-1)}(t) upper and lower enveloping curve;

Step c2, the described survival function h of calculating _{N (k-1)}(t) upper and lower enveloping curve is at time t=1, and 2 ..., the average in the N

Step c3, obtain system model error term estimated value Survival function after screening through the k time in the residue trend function of decomposing through the n time intrinsic mode function

5. a kind of satellite failure Forecasting Methodology of decomposing based on predictive filtering and empirical modal according to claim 3 is characterized in that among the step f estimated value to system model error term

Carry out 4 intrinsic mode functions and decompose, obtain 4 eigenmode state function component IMF1, IMF2, IMF3 and IMF4:{c ₁(t), c ₂(t), c ₃(t), c ₄(t) }; And acquisition is through the remaining residual error function r of 4 intrinsic mode function decomposition ₄(t).

6. a kind of satellite failure Forecasting Methodology according to claim 3 based on predictive filtering and empirical modal decomposition, it is characterized in that utilizing in the step 3 Time series analysis method to set up the model of the fault trend of the residual component that obtains about step 2, finish the forecast of small gentle change fault and the process of detection and be:

Residual component as the trend signal, is set up its autoregressive model, and expression formula is:

r(t)＝a ₁r(t-1)+a ₂r(t-2)+…+a _pr(t-p)+ε(t)

Wherein, r (t) is the time series of residual error function, and p is a model order, and ai is a model parameter, and ε (t) is a white noise sequence;

Observed reading is r (0), r (1) ..., r (N), the order of forecast model are p, get following system of equations by the autoregressive model expression formula:

$\left\{\begin{array}{c}r\left(p\right)=a_{1}r(p-1)+a_{2}r(p-2)+...+a_{p}r\left(0\right)+\mathrm{\&epsiv;}\left(p\right)\\ r(p+1)=a_{1}r\left(p\right)+a_{2}r(p-1)+...+a_{p}r\left(1\right)+\mathrm{\&epsiv;}(p+1)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ r\left(N\right)=a_{1}r(N-1)+a_{2}r(N-2)+...+a_{p}r(N-p)+\mathrm{\&epsiv;}\left(N\right)\end{array}\right.$

Order:

$R=\left[\begin{array}{c}r\left(p\right)\\ r(p+1)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ r\left(N\right)\end{array}\right],$ $A=\left[\begin{array}{c}{a}_1\\ a_2\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ a_p\end{array}\right],$

$B=\left[\begin{array}{cccc}r(p-1)& r(p-2)& \&CenterDot;\&CenterDot;\&CenterDot;& r\left(0\right)\\ r\left(p\right)& r(p-1)& \&CenterDot;\&CenterDot;\&CenterDot;& r\left(1\right)\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ r(N-1)& r(N-2)& \&CenterDot;\&CenterDot;\&CenterDot;& r(N-p)\end{array}\right],$ $\mathrm{\&Delta;}=\left[\begin{array}{c}\mathrm{\&epsiv;}\left(p\right)\\ \mathrm{\&epsiv;}(p+1)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ \mathrm{\&epsiv;}\left(N\right)\end{array}\right]$

Then the matrix form of following formula is:

R＝BA+Δ

The least square solution of model parameter is:

A=(B ^TB) ^-1B ^TR obtains the model parameter A of autoregressive model, utilizes model parameter A to carry out failure prediction and detection.

Images