CN101859146B

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

Info

Publication number: CN101859146B
Application number: CN2010102287440A
Authority: CN
Inventors: 沈毅; 张迎春; 王振华
Original assignee: Harbin Institute of Technology Shenzhen
Current assignee: Harbin Institute of Technology Shenzhen
Priority date: 2010-07-16
Filing date: 2010-07-16
Publication date: 2011-11-30
Anticipated expiration: 2030-07-16
Also published as: CN101859146A

Abstract

A satellite fault prediction method based on predictive filtering and empirical mode decomposition, which involves the safe operation of space satellites and fault prediction methods, and solves the problem that traditional satellite fault prediction and diagnosis methods are seriously affected by noise and cannot accurately predict fault trends The specific process is as follows: 1. Using the nonlinear attitude dynamics relationship of the satellite, the error of the satellite control system is estimated by the method of predictive filtering, and the error item of the system model is obtained; 2. The error item of the system model obtained in step 1 is empirically modeled State decomposition to obtain the first n-order intrinsic mode function IMF components and residual components; third, use the time series analysis method to establish a fault trend model for the residual components obtained in step 2, and complete the prediction and detection of small and slowly changing faults . It is used in the field of fault diagnosis of 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:

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:

\overset{\cdot}{x} (t) = f [x (t)] + g [x (t)] d (t)

{\tilde{y}}_{k} = c [x (t_{k})] + v_{k}

The predictive filtering equation is:

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

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

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 + Δt) = \hat{y} (t) + z [\hat{x} (t), Δt] + Λ (Δt) S [\hat{x} (t)] d (t)

Δ 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} (t), Δt] = Σ_{k = 1}^{p_{i}} \frac{{Δt}^{k}}{k!} {L^{k}}_{f} (c_{i})

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:

λ_{ii} = \frac{{Δt}^{p_{i}}}{p_{i}!},

i＝1，2，…，m

Matrix

S [\hat{x} (t)] &Element; R^{m \times q},

The element that its i is capable is:

s_{i} = {L_{g_{1}} [L_{f}^{p_{i} - 1} (c_{i})], L, L_{g_{q}} [L_{f}^{p_{i} - 1} (c_{i})]},

i＝1，2，…，m

Get performance index function:

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

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} (t) = - {{[Λ (Δt) S (\hat{x})]}^{T} Q^{- 1} Λ (Δt) S (\hat{x}) + W}^{- 1} [Λ (Δt) S (\hat{x})] Q^{- 1} [z (\hat{x}, Δt) - \tilde{y} (t + Δt) + \hat{y} (t)] .

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} (t) = d (t) - Σ_{α = 1}^{n} c_{α} (t);

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_{nk} (t) = h_{n (k - 1)} (t) - {\overset{&OverBar;}{m}}_{n (k - 1)} (t) .

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:

\{\begin{matrix} r (p) = a_{1} r (p - 1) + a_{2} r (p - 2) + . . . + a_{p} r (0) + ϵ (p) \\ r (p + 1) = a_{1} r (p) + a_{2} r (p - 1) + . . . + a_{p} r (1) + ϵ (p + 1) \\ \cdot \\ \cdot \\ \cdot \\ r (N) = a_{1} r (N - 1) + a_{2} r (N - 2) + . . . + a_{p} r (N - p) + ϵ (N) \end{matrix}

Order:

R = [\begin{matrix} r (p) \\ r (p + 1) \\ \cdot \\ \cdot \\ \cdot \\ r (N) \end{matrix}],

A = [\begin{matrix} a_{1} \\ a_{2} \\ \cdot \\ \cdot \\ \cdot \\ a_{p} \end{matrix}],

B = [\begin{matrix} r (p - 1) & r (p - 2) & \cdot \cdot \cdot & r (0) \\ r (p) & r (p - 1) & \cdot \cdot \cdot & r (1) \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ r (N - 1) & r (N - 2) & \cdot \cdot \cdot & r (N - p) \end{matrix}],

Δ = [\begin{matrix} ϵ (p) \\ ϵ (p + 1) \\ \cdot \\ \cdot \\ \cdot \\ ϵ (N) \end{matrix}]

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:

\{\begin{matrix} {\overset{\cdot}{x}}_{1} (t) = (\frac{I_{y} - I_{z}}{I_{x}}) x_{2} x_{3} + (\frac{1}{I_{x}}) T_{x} + (\frac{1}{T_{x}}) f_{x} (t) \\ {\overset{\cdot}{x}}_{2} (t) = (\frac{I_{x} - I_{z}}{I_{y}}) x_{1} x_{3} + (\frac{I}{I_{y}}) T_{y} + (\frac{1}{T_{y}}) f_{y} (t) \\ {\overset{\cdot}{x}}_{3} (t) = (\frac{I_{x} - I_{y}}{I_{z}}) x_{1} x_{2} + (\frac{I}{I_{z}}) T_{z} + (\frac{1}{I_{z}}) f_{z} (t) \end{matrix}

\tilde{y} (t_{k}) = x (t_{k}) + v (t_{k})

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 (t) = - {{[Λ (Δt) S (\hat{x})]}^{T} Q^{- 1} Λ (Δt) S (\hat{x}) + W}^{- 1} \cdot

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

Wherein

z (\hat{x}, Δt) = [\begin{matrix} (\frac{I_{y} - I_{z}}{I_{x}}) {\hat{x}}_{2} {\hat{x}}_{3} + (\frac{1}{I_{x}}) T_{x} \\ (\frac{I_{z} - I_{x}}{I_{y}}) {\hat{x}}_{1} {\hat{x}}_{3} + (\frac{I}{I_{y}}) T_{y} \\ (\frac{I_{x} - I_{y}}{I_{z}}) {\hat{x}}_{1} {\hat{x}}_{2} + (\frac{I}{I_{z}}) T_{z} \end{matrix}] Δt,

Λ (Δt) = [\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix}] Δt,

S (\hat{x}) = [\begin{matrix} \frac{1}{I_{x}} & 0 & 0 \\ 0 & \frac{1}{I_{y}} & 0 \\ 0 & 0 & \frac{1}{I_{z}} \end{matrix}]

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

1. a satellite failure prediction method based on predictive filtering and empirical mode decomposition, is characterized in that concrete process is as follows:

Step 1: Estimate the error of the satellite control system by using the nonlinear attitude dynamic relationship of the satellite, and use the method of predictive filtering to obtain the error term of the system model;

Step 2: Perform empirical mode decomposition on the error item of the system model obtained in step 1 to obtain the IMF component and residual component of the former E-order eigenmode function;

Step 3: Use the time series analysis method to establish a model of the fault trend of the residual component obtained in Step 2, and complete the prediction and detection of small and slowly changing faults.

2. a kind of satellite fault prediction method based on predictive filtering and empirical mode decomposition according to claim 1, it is characterized in that utilize satellite nonlinear attitude dynamics relation in the step 1, adopt the method for predictive filtering to satellite control system error The process of estimating and obtaining the error term of the system model is:

It is assumed that the model error of the satellite control system is composed of satellite actuator failure and model uncertainty, and its nonlinear state equation and measurement equation are respectively:

Equation of state:

\overset{\cdot &Center Dot;}{x x} ((t t)) = = f f [[x x ((t t))]] + + g g [[x x ((t t))]] d d ((t t))

{\overset{~ ~}{y the y}}_{k k} = = c c [[x x (({t t}_{k k}))]] + + {v v}_{k k}

The prediction filter equation is:

\overset{^^}{x x} ((t t)) = = f f [[\overset{^^}{x x} ((t t))]] + + g g [[\overset{^^}{x x} ((t t))]] \overset{^^}{d d} ((t t))

\overset{^^}{y the y} ((t t)) = = c c [[\overset{^^}{x x} ((t t))]]

where x(t)∈R ^s is the state vector,

is the estimated value of the state vector, f∈R ^s is a differentiable state function,

is the known model error distribution matrix,

is the measurement function vector, d(t)∈R ^q is the unknown model error,

for its estimated value,

Indicates the measurement output of the actual system at time t _k in discrete form, v _k ∈ R ^m is the measurement noise, and v _k is set to be Gaussian white noise with zero mean and covariance matrix Q ∈ R ^m×m ;

Taylor expand the measurement function at time t+Δt to get:

\overset{^^}{y the y} ((t t + + Δt Δt)) = = \overset{^^}{y the y} ((t t)) + + z z [[\overset{^^}{x x} ((t t)),, Δt Δt]] + + Λ Λ ((Δt Δt)) S S [[\overset{^^}{x x} ((t t))]] d d ((t t))

Among them, Δt=t _k+1 -t _k is the sampling period; the sampling period is a constant value, y(t+Δt)=y _k+1 , the matrix

The i-th element of is:

{z z}_{i i} [[\overset{^^}{x x} ((t t)),, Δt Δt]] = = {Σ Σ}_{k k = = 11}^{{p p}_{i i}} \frac{{Δt Δt}^{k k}}{k k!!} {L L}^{k k}_{f f} (({c c}_{i i}))

Where p _i is the differential order when d(t) appears for the first time in the Taylor expansion,

is the k-order Lie derivative of _ci ;

Matrix A(Δt) is a diagonal matrix whose diagonal elements are:

λ_{i} = \frac{{Δt}^{p_{i}}}{p_{i}!},

i=1,2,...,m

matrix

S [\hat{x} (t)] &Element; R^{m \times q},

The elements of the i-th row are:

{the s}_{i} = {L_{g_{1}} [L_{f}^{p_{i} - 1} (c_{i})], L, L_{g_{q}} [L_{f}^{p_{i} - 1} (c_{i})]},

i=1,2,...,m

Take the performance index function:

J J [[d d ((t t))]] = = \frac{11}{22} {{\overset{~ ~}{y the y} ((t t + + Δt Δt)) - - \overset{~ ~}{y the y} ((t t + + Δt Δt))}}^{T T} {Q Q}^{- - 11} {{\overset{~ ~}{y the y} ((t t + + Δt Δt)) - - \overset{~ ~}{y the y} ((t t + + Δt Δt))}} + + \frac{11}{22} {d d}^{T T} ((t t)) Wd Wd ((t t))

Where W∈R ^q×q is a positive semidefinite weighting matrix, and the gradient optimization algorithm is used to optimize the performance index, and the estimation of the model error term is obtained as:

\overset{^^}{d d} ((t t)) = = - - {{{[[Λ Λ ((Δt Δt)) S S ((\overset{^^}{x x}))]]}^{T T} {Q Q}^{- - 11} Λ Λ ((Δt Δt)) S S ((\overset{^^}{x x})) + + W W}}^{- - 11} \cdot &Center Dot;

[[Λ Λ ((Δt Δt)) S S ((\overset{^^}{x x}))]] {Q Q}^{- - 11} [[z z ((\overset{^^}{x x},, Δt Δt)) - - \overset{~ ~}{y the y} ((t t + + Δt Δt)) + + \overset{^^}{y the y} ((t t))]] . .

3. a kind of satellite fault prediction method based on predictive filtering and empirical mode decomposition according to claim 2, it is characterized in that step 2 obtains the process of preceding E order eigenmode function IMF component and residual component as:

Set the estimated error term of the system model to be

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

Step a, IMF decomposition process initialization: n=1, and satisfy the relation

Established, where r _n-1 (t) is the residual function remaining after the (n-1)th decomposition;

Step b, initialization of the screening process: k=1, and the relation h _n(k-1) (t)=r _n-1 (t) is satisfied, where h _n(k-1) (t) is the nth time The remaining function after the k-1th screening in the IMF decomposition;

Step c. Obtain the estimated value of the error term of the system model according to the screening procedure

The residual function h _nk (t) after the kth screening of the remaining residual functions after the nth eigenmode function decomposition;

Step d, using the standard deviation criterion to judge whether the obtained residual function h _nk (t) satisfies the condition of the intrinsic mode function, namely

Whether it is less than the threshold T, 0.2≤T≤0.3;

If the judgment result is yes, execute step e, if the judgment result is no, then k=k+1, return to execute step c,

Step e, obtaining the nth intrinsic mode function IMF component c _n (t)=h _nk (t);

Step f. Obtain the estimated value of the error term of the system model The remaining residual function after the nth eigenmode function decomposition

r_{no} (t) = \hat{d} (t) - Σ_{α = 1}^{no} c_{α} (t);

Step g, set n=n+1, return to step b, until n=E is satisfied to obtain the IMF component and residual component of the first E order intrinsic mode function.

4. a kind of satellite fault prediction method based on predictive filtering and empirical mode decomposition according to claim 3, is characterized in that step c obtains system model error term estimated value according to screening procedure

The process of the residual function h _nk (t) after the k-th screening of the remaining residual functions after the n-th eigenmode function decomposition is:

Step c1, using the cubic spline function to obtain the estimated value of the error term of the system model The upper and lower envelope curves of the residual function h _n(k-1) (t) after the k-1 screening in the residual trend function of the n-th eigenmode function decomposition;

Step c2, calculating the mean value of the upper and lower envelope curves of the residual function h _n(k-1) (t) at time t=1, 2, ..., N

Step c3. Obtain the estimated value of the error term of the system model The residual function after the kth screening in the residual trend function after the nth eigenmode function decomposition

5. a kind of satellite fault prediction method based on predictive filtering and empirical mode decomposition according to claim 3 is characterized in that in the step f, to the estimated value of system model error term

Decompose the intrinsic mode function four times to obtain four intrinsic mode function components IMF1, IMF2, IMF3 and IMF4: {c ₁ (t), c ₂ (t), c ₃ (t), c ₄ (t )}; and obtain the remaining residual function r ₄ (t) after four eigenmode function decompositions.

6. a kind of satellite fault prediction method based on predictive filtering and empirical mode decomposition according to claim 3 is characterized in that in step 3, utilizes time series analysis method to set up the model about the fault tendency of the residual error component that step 2 obtains , the process of completing the prediction and detection of small and slowly changing faults is:

The residual component is used as a trend signal, and its autoregressive model is established, and the expression is:

r(t)＝a ₁ r(t-1)+a ₂ r(t-2)+...+a _p r(tp)+ε(t)

Among them, r(t) is the time series of the residual function, p is the model order, ai is the model parameter, ε(t) is the white noise sequence;

The observed values are r(0), r(1), ..., r(N), and the order of the prediction model is p. The following equations are obtained from the expression of the autoregressive model:

\{\begin{matrix} r r ((p p)) = = {a a}_{11} r r ((p p - - 11)) + + {a a}_{22} r r ((p p - - 22)) + + . . . . . . + + {a a}_{p p} r r ((00)) + + ϵ ϵ ((p p)) \\ r r ((p p + + 11)) = = {a a}_{11} r r ((p p)) + + {a a}_{22} r r ((p p - - 11)) + + . . . . . . + + {a a}_{p p} r r ((11)) + + ϵ ϵ ((p p + + 11)) \\ \cdot &Center Dot; \\ \cdot &Center Dot; \\ \cdot &Center Dot; \\ r r ((N N)) = = {a a}_{11} r r ((N N - - 11)) + + {a a}_{22} r r ((N N - - 22)) + + . . . . . . + + {a a}_{p p} r r ((N N - - p p)) + + ϵ ϵ ((N N)) \end{matrix}

make:

R R = = [\begin{matrix} r r ((p p)) \\ r r ((p p + + 11)) \\ \cdot &Center Dot; \\ \cdot &Center Dot; \\ \cdot &Center Dot; \\ r r ((N N)) \end{matrix}],,

A A = = [\begin{matrix} {a a}_{11} \\ {a a}_{22} \\ \cdot &Center Dot; \\ \cdot &Center Dot; \\ \cdot &Center Dot; \\ {a a}_{p p} \end{matrix}],,

B B = = [\begin{matrix} r r ((p p - - 11)) & r r ((p p - - 22)) & \cdot &Center Dot; \cdot &Center Dot; \cdot \cdot & r r ((00)) \\ r r ((p p)) & r r ((p p - - 11)) & \cdot \cdot \cdot \cdot \cdot \cdot & r r ((11)) \\ \cdot &Center Dot; & \cdot \cdot & \cdot &Center Dot; & \cdot \cdot \\ \cdot &Center Dot; & \cdot &Center Dot; & \cdot \cdot & \cdot &Center Dot; \\ \cdot \cdot & \cdot &Center Dot; & \cdot \cdot & \cdot &Center Dot; \\ r r ((N N - - 11)) & r r ((N N - - 22)) & \cdot \cdot \cdot \cdot \cdot \cdot & r r ((N N - - p p)) \end{matrix}],,

Δ Δ = = [\begin{matrix} ϵ ϵ ((p p)) \\ ϵ ϵ ((p p + + 11)) \\ \cdot \cdot \\ \cdot &Center Dot; \\ \cdot &Center Dot; \\ ϵ ϵ ((N N)) \end{matrix}]

Then the matrix form of the above formula is:

R=BA+Δ

The least squares solution for the model parameters is:

A=(B ^T B) ^-1 B ^T R Obtain the model parameter A of the autoregressive model, and use the model parameter A to carry out fault prediction and detection.