CN105446352A - Proportion guide law recognition filtering method - Google Patents

Abstract

The invention relates to a proportion guide law recognition filtering method. The invention aims at solving problems that a conventional MMAE filer supposes that a model is not accurate because a PN guide law navigation constant is known and does not give consideration to the saturation condition of a pursuer controller, and the estimation precision is low. The method comprises the steps: firstly building state equations of PN guide law saturation and motion models in pitching and yam planes; secondly calculating a system state transfer matrix and designing a Kalman filter equation; thirdly calculating the posterior probabilities of the PN guide law motion and saturation motion models of a pursuer at the current moment; finally calculating the switching moment of the pursuer motion model according to the posterior probabilities, employing the estimation result of the saturation motion model Kalman filter equation, or else, employing the estimation result of a PN guide law motion Kalman filter equation. The method is suitable for the field of spaceflight.

Description

A kind of proportional navigation law identification filtering method

Technical field

The present invention relates to proportional navigation law identification filtering method.

Background technology

In exoatmosphere ballistic missile penetration application, in order to improve the penetration ability of trajectory, ballistic missile can discharge defending missile accompanying flying synchronized with it in exoatmosphere.Controlled and its collision by Guidance Law when defending missile finds interception guided missile, thus improve the penetraton probability of ballistic missile.The design of the optimum evasive strategy of ballistic missile and the optimal guidance law of defending missile all supposes that the Guidance Law of interception guided missile is known.Although control Guidance Law based on sliding formwork the motor-driven of interception guided missile to be processed as unknown external disturbance, if the acceleration of interception guided missile can be estimated, the performance of sliding mode guidance significantly can be improved.Therefore the identification problem studying interception guided missile Guidance Law has actual meaning.

The people such as Shaferman devise one and restrain based on the aircraft Initiative Defense cooperative guidance of MMAE (Multiple-ModelAdaptiveEstimator) in the article " CooperativeMultiple-ModelAdaptiveGuidanceforanAircraftDe fendingMissile " of " JOURNALOFGUIDANCECONTROLANDDYNAMICS " magazine the 6th phase in 2010.Suppose interception guided missile adoption rate guiding (PN) Guidance Law, strengthen the one in proportional guidance (APN) Guidance Law or optimal guidance law (OGL), adopt multi-model self-adapting filtering (MMAE) to carry out identification to the Guidance Law of interception guided missile.But there is two problems in this work.1) the method supposition interception guided missile have employed known Guidance Law constant N.Which has limited the scope of application of this MMAE.In actual scene, the navigation constant of interception guided missile Guidance Law is unknown.2) the method does not consider the situation of interception guided missile actuator saturation.Generally, at the guidance phases initial stage, interception guided missile can temporarily enter actuator saturation state owing to compensating initial alignment error fast, and after this error is compensated, controller exits state of saturation.Do not consider the situation of actuator saturation, MMAE filter model will be made inaccurate, and estimated accuracy is low.

Summary of the invention

The object of the invention is to solve following problem: 1) existing method supposition interception guided missile have employed known Guidance Law constant N, limits the scope of application of MMAE.And in actual scene, the navigation constant of interception guided missile Guidance Law is unknown; 2) may there is saturated conditions at guidance initial time in interception guided missile actuator.Existing method does not consider above actual conditions, makes MMAE filter model inaccurate, the problem that estimated accuracy is low, for above problem, the present invention proposes a kind of proportional navigation law identification filtering method.

Above-mentioned goal of the invention is achieved through the following technical solutions:

Step one: pursuer has two kinds of motion models: PN Guidance Law motion model when saturated motion model during controller saturation and controller unsaturation; Set up PN Guidance Law motion model at the state equation of pitch plane, PN Guidance Law motion model at the state equation of the plane of yaw, controller saturation motion model at the state equation of pitch plane and the controller saturation motion model state equation at the plane of yaw; Evader is ballistic missile; Pursuer is interception guided missile; PN Guidance Law is proportional navigation law;

Step 2: according to the PN Guidance Law motion model of step one state equation at pitch plane, PN Guidance Law motion model is at the state equation of the plane of yaw, controller saturation motion model is at the state equation of pitch plane and the controller saturation motion model state equation at the plane of yaw, design PN Guidance Law motion model is at the state-transition matrix of pitch plane, PN Guidance Law motion model is at the state-transition matrix of the plane of yaw, controller saturation motion model is at the state-transition matrix of pitch plane, controller saturation motion model is at the state-transition matrix of the plane of yaw,

Step 3: according to the PN Guidance Law motion model of the step 2 state-transition matrix at pitch plane, PN Guidance Law motion model is at the state-transition matrix of the plane of yaw, controller saturation motion model is at the state-transition matrix of pitch plane, controller saturation motion model is at the state-transition matrix of the plane of yaw, design PN Guidance Law motion model is at the Kalman filter equation of pitch plane, PN Guidance Law motion model is at the Kalman filter equation of the plane of yaw, controller saturation motion model is at the Kalman filter equation of pitch plane and the controller saturation motion model Kalman filter equation at the plane of yaw, Kalman filter equation is Kalman filter equation,

Step 4: at the Kalman filter equation of the plane of yaw, calculate the posterior probability of pursuer current time PN Guidance Law motion model and the posterior probability of controller saturation motion model at the Kalman filter equation of pitch plane and controller saturation motion model at the Kalman filter equation of the plane of yaw, controller saturation motion model at the Kalman filter equation of pitch plane, PN Guidance Law motion model according to the PN Guidance Law motion model of step 3;

Step 5: the posterior probability of pursuer current time PN Guidance Law motion model obtained according to step 4 and the posterior probability of saturated motion model, calculates the moment of PN Guidance Law motion model when pursuer is switched to controller unsaturation from saturated motion model during controller saturation; Adopted controller saturation motion model at the Kalman filter equation of pitch plane and the controller saturation motion model estimated result at the Kalman filter equation of the plane of yaw before the moment of the PN Guidance Law motion model when pursuer is switched to controller unsaturation from saturated motion model during controller saturation; Otherwise adopt PN Guidance Law motion model at the Kalman filter equation of pitch plane and the PN Guidance Law motion model estimated result at the Kalman filter equation of the plane of yaw.

Invention effect

The proportional navigation law identification wave filter that the present invention proposes have employed multiple model filtering method.Devise the extended Kalman filter respectively for PN Guidance Law motion model and saturated motion model.Wherein PN Guidance Law motion model is using the navigation constant of PN Guidance Law as state to be estimated, thus solves the actual conditions of PN Guidance Law navigation constant the unknown.Above-mentioned two class extended Kalman filter parallel runnings, calculate the probability of two models by Bayesian inference framework, thus pick out the motion model of interception guided missile employing, and then estimate acceleration and the PN Guidance Law navigation constant of interception guided missile.Which solve the problem that interception guided missile may be saturated at its controller of guidance initial stage.The filtering method that the present invention proposes improves about 40-50% relative to traditional its estimated accuracy of EKF method of PN Guidance Law model that only adopts.

Consider pursuer controller in guidance process and produce saturated possibility in the early stage, employ multiple model filtering method and identification has been carried out to pursuer Guidance Law, estimate to the navigation constant under PN Guidance Law with at pitch plane and plane of yaw acceleration.Identification can be carried out more accurately to Guidance Law.In an embodiment of the present invention, if MMAE wave filter is Var to the square error that pursuer acceleration estimation produces _mMAE; Omnidistance use PN Guidance Law is Var to the square error that pursuer acceleration estimation produces _pN.If α=Var _mMAE/ Var _pN.As 0< α <1, to pursuer acceleration estimation, MMAE wave filter comparatively illustrates that the omnidistance PN Guidance Law that uses is more accurate to pursuer acceleration estimation; α is less, illustrates to use PN Guidance Law concerning pursuer acceleration estimation relative to whole process, and MMAE is more accurate to pursuer acceleration estimation.

When evader is without time motor-driven, at pursuer saturation stage, to the component a of pursuer acceleration in Evader inertial coordinates system y-axis _yestimation time α=0.59; To the component a of pursuer acceleration in Evader inertial coordinates system z-axis _zestimation time α=0.56.

When evader carry out constant value motor-driven time, at pursuer saturation stage, to pursuer acceleration a _yestimation time α=0.55; To pursuer acceleration a _zestimation time α=0.62.

When evader carry out sinusoidal motor-driven time, at pursuer saturation stage, to pursuer acceleration a _yestimation time α=0.56; To pursuer acceleration a _zestimation time α=0.60.

Accompanying drawing explanation

Fig. 1 is pursuer and evader relative motion relation schematic diagram, and E is the coordinate of evader under evader inertial system, and P is the coordinate of pursuer under evader inertial system, q _{p ε}for the sight line angle of pitch, q _{p β}for sight line crab angle, r is the relative distance of evader to pursuer, LOS ₀for the initial sight line of evader to pursuer;

Fig. 2 is to Pursuer pitch plane navigation constant N in the motor-driven situation of Evader constant value _εestimated result figure;

Fig. 3 is to Pursuer plane of yaw navigation constant N in the motor-driven situation of Evader constant value _βestimated result figure;

Fig. 4 uses PN Guidance Law Kalman identification wave filter to the component a of Pursuer acceleration in Evader inertial coordinates system y-axis under the inorganic emotionally condition of Evader _yestimated result figure, g are acceleration of gravity;

Fig. 5 uses PN Guidance Law Kalman identification wave filter to the component a of Pursuer acceleration in Evader inertial coordinates system z-axis under the inorganic emotionally condition of Evader _zestimated result figure, g are acceleration of gravity;

Fig. 6 uses MMAE wave filter to Pursuer acceleration a under the inorganic emotionally condition of Evader _yestimated result figure;

Fig. 7 uses MMAE wave filter to Pursuer acceleration a under the inorganic emotionally condition of Evader _zestimated result figure;

Fig. 8 uses PN Guidance Law Kalman identification wave filter to Pursuer acceleration a in the motor-driven situation of Evader constant value _yestimated result figure;

Fig. 9 uses PN Guidance Law Kalman identification wave filter to Pursuer acceleration a in the motor-driven situation of Evader constant value _zestimated result figure;

Figure 10 uses MMAE wave filter to Pursuer acceleration a in the motor-driven situation of Evader constant value _yestimated result figure;

Figure 11 uses MMAE wave filter to Pursuer acceleration a in the motor-driven situation of Evader constant value _zestimated result figure;

Figure 12 is that Evader sine mechanism emotionally uses PN Guidance Law Kalman identification wave filter to Pursuer acceleration a under condition _yestimated result figure;

Figure 13 is that Evader sine mechanism emotionally uses PN Guidance Law Kalman identification wave filter to Pursuer acceleration a under condition _zestimated result figure;

Figure 14 is that Evader sine mechanism emotionally uses MMAE wave filter to Pursuer acceleration a under condition _yestimated result figure;

Figure 15 is that Evader sine mechanism emotionally uses MMAE wave filter to Pursuer acceleration a under condition _zestimated result figure.

Embodiment

Embodiment one: a kind of proportional navigation law identification filtering method of present embodiment, specifically prepare according to following steps:

Step one: pursuer has two kinds of motion models: PN Guidance Law motion model when saturated motion model during controller saturation and controller unsaturation; Set up PN Guidance Law motion model at the state equation of pitch plane, PN Guidance Law motion model at the state equation of the plane of yaw, controller saturation motion model at the state equation of pitch plane and the controller saturation motion model state equation at the plane of yaw; Evader is ballistic missile; Pursuer is interception guided missile; PN Guidance Law is proportional navigation law;

Step 2: according to the PN Guidance Law motion model of step one state equation at pitch plane, PN Guidance Law motion model is at the state equation of the plane of yaw, controller saturation motion model is at the state equation of pitch plane and the controller saturation motion model state equation at the plane of yaw, calculate the state-transition matrix of PN Guidance Law motion model at pitch plane, PN Guidance Law motion model is at the state-transition matrix of the plane of yaw, controller saturation motion model is at the state-transition matrix of pitch plane, controller saturation motion model is at the state-transition matrix of the plane of yaw,

Step 3: according to the PN Guidance Law motion model of the step 2 state-transition matrix at pitch plane, PN Guidance Law motion model is at the state-transition matrix of the plane of yaw, controller saturation motion model is at the state-transition matrix of pitch plane, controller saturation motion model is at the state-transition matrix of the plane of yaw, design PN Guidance Law motion model is at the Kalman filter equation of pitch plane, PN Guidance Law motion model is at the Kalman filter equation of the plane of yaw, controller saturation motion model is at the Kalman filter equation of pitch plane and the controller saturation motion model Kalman filter equation at the plane of yaw, Kalman filter equation is Kalman filter equation,

Step 4: at the Kalman filter equation of the plane of yaw, calculate the posterior probability of pursuer current time PN Guidance Law motion model and the posterior probability of controller saturation motion model at the Kalman filter equation of pitch plane and controller saturation motion model at the Kalman filter equation of the plane of yaw, controller saturation motion model at the Kalman filter equation of pitch plane, PN Guidance Law motion model according to the PN Guidance Law motion model of step 3;

Step 5: the posterior probability of pursuer current time PN Guidance Law motion model obtained according to step 4 and the posterior probability of saturated motion model, calculates the moment of PN Guidance Law motion model when pursuer is switched to controller unsaturation from saturated motion model during controller saturation; Adopted controller saturation motion model at the Kalman filter equation of pitch plane and the controller saturation motion model estimated result at the Kalman filter equation of the plane of yaw before the moment of the PN Guidance Law motion model when pursuer is switched to controller unsaturation from saturated motion model during controller saturation; Otherwise adopt PN Guidance Law motion model at the Kalman filter equation of pitch plane and the PN Guidance Law motion model estimated result at the Kalman filter equation of the plane of yaw.

Embodiment two: present embodiment and embodiment one unlike: in described step one, pursuer has two kinds of motion models: PN Guidance Law motion model when saturated motion model during controller saturation and controller unsaturation; Set up PN Guidance Law motion model at the state equation of pitch plane, PN Guidance Law motion model at the state equation of the plane of yaw, controller saturation motion model at the state equation of pitch plane and the controller saturation motion model state equation at the plane of yaw; Detailed process is:

In missile breakthrough scene, prominent anti-ballistic missile will escape the interception of the interception guided missile that the other side's missile defense systems is launched, and is escape side, is called evader, and its relevant all physical quantity all contains subscript e; Interception guided missile is the side of chasing, and is called pursuer, and its relevant all physical quantity all contains subscript p; The Guidance and control supposing interception guided missile can decoupling zero be fore-and-aft plane and lateral plane.Fig. 1 depicts the relative motion relation of pursuer and evader at pitch plane.

Evader LOS coordinate system defines; LOS coordinate system o ' x ₄y ₄z ₄initial point o ' is positioned at the target seeker centre of gyration, o ' x ₄axle is consistent with target-guided missile sight line, and pointing to target by the target seeker centre of gyration is just, o ' y ₄axle is positioned at and comprises o ' x ₄in the plummet face of axle, with o ' x ₄axle is vertical, and pointing to top is just, o ' z ₄axle is determined by the right-hand rule; And the terminal guidance inertial coordinates system of evader is defined as the evader LOS coordinate system o of guidance initial time ₀x ₀y ₀z ₀;

Pursuer needed to eliminate alignment error at the guidance initial stage, and this likely causes pursuer controller to enter state of saturation, and after error concealment, pursuer controller can exit state of saturation; Therefore pursuer has two kinds of motion models: PN Guidance Law motion model when saturated motion model during controller saturation and controller unsaturation; Provide the model that pursuer two kinds moves at pitch plane and the plane of yaw below;

The PN Guidance Law equation of motion modeling of pursuer, r _xfor the component of evader and pursuer relative position under evader inertial system x-axis; r _yfor the component of evader and pursuer relative position under evader inertial system y-axis; r _zfor the component of evader and pursuer relative position under evader inertial system z-axis; Evader is ballistic missile, and pursuer is interception guided missile;

r _x＝x _p-x _er _y＝y _p-y _er _z＝z _p-z _e(1)

In formula, [x _p, y _p, z _p] ^t[x _e, y _e, z _e] ^tbe the position of pursuer and evader under evader inertial system respectively, T is transposition, x _pfor the position of pursuer under evader inertial system x-axis, y _pfor the position of pursuer under evader inertial system y-axis, z _pfor the position of pursuer under evader inertial system z-axis, x _efor the position of evader under evader inertial system x-axis, y _efor the position of evader under evader inertial system y-axis, z _efor the position of evader under evader inertial system z-axis;

PN Guidance Law motion model at the state equation of pitch plane is

${\stackrel{\&CenterDot;}{r}}_x=v_x$

${\stackrel{\&CenterDot;}{r}}_y=v_y$

${\stackrel{\&CenterDot;}{v}}_x=a_{px}-a_{ex}$

${\stackrel{\&CenterDot;}{v}}_y=a_{py}-a_{ey}---\left(2\right)$

${\stackrel{\&CenterDot;}{a}}_{px}=-\frac{a_{px}}{\mathrm{\&tau;}}$

${\stackrel{\&CenterDot;}{a}}_{py}=N_{\mathrm{\&epsiv;}}\frac{(r_x{v}_x+r_y{v}_y)(r_x{v}_y-r_y{v}_x)}{{(r_x^2+r_y^2)}^{\frac{3}{2}}\mathrm{\&tau;}}-\frac{a_{py}}{\mathrm{\&tau;}}$

${\stackrel{\&CenterDot;}{N}}_{\mathrm{\&epsiv;}}=0$

PN Guidance Law motion model at the state equation of the plane of yaw is

${\stackrel{\&CenterDot;}{r}}_x=v_x$

${\stackrel{\&CenterDot;}{r}}_z=v_z$

${\stackrel{\&CenterDot;}{v}}_x=a_{px}-a_{ex}$

${\stackrel{\&CenterDot;}{v}}_z=a_{pz}-a_{ez}---\left(3\right)$

${\stackrel{\&CenterDot;}{a}}_{px}=-\frac{a_{px}}{\mathrm{\&tau;}}$

${\stackrel{\&CenterDot;}{a}}_{pz}=-\frac{N_{\mathrm{\&beta;}}}{\mathrm{\&tau;}}\frac{(r_x{v}_x+r_z{v}_z)(r_z{v}_x-r_x{v}_z)}{{(r_x^2+r_z^2)}^{\frac{3}{1}}}-\frac{a_{pz}}{\mathrm{\&tau;}}$

${\stackrel{\&CenterDot;}{N}}_{\mathrm{\&beta;}}=0$

In formula, v _xfor the component of evader and pursuer relative velocity under evader inertial system x-axis, v _yfor the component of evader and pursuer relative velocity under evader inertial system y-axis, v _zfor the component of evader and pursuer relative velocity under evader inertial system z-axis; a _pxfor the component of acceleration under evader inertial system x-axis of pursuer, a _pyfor the component of acceleration under evader inertial system y-axis of pursuer, a _pzfor the component of acceleration under evader inertial system z-axis of pursuer; a _exfor the component of acceleration under evader inertial system x-axis of pursuer, a _eyfor the component of acceleration under evader inertial system y-axis of pursuer, a _ezfor the component of acceleration under evader inertial system z-axis of pursuer; N _εthe navigation constant of pursuer at pitch plane; N _βthe navigation constant of pursuer at the plane of yaw; τ is time constant, for r _xfirst order derivative, for r _yfirst order derivative, for r _zfirst order derivative, for v _xfirst order derivative, for v _yfirst order derivative, for v _zfirst order derivative, for a _pxfirst order derivative, for a _pyfirst order derivative, for a _pzfirst order derivative, for N _εfirst order derivative, for N _βfirst order derivative;

Suppose that evader can obtain the position of pursuer under inertial system, and then calculate the position of relative self (i.e. evader); The measurement equation calculating the position of self is relatively

h＝[r _x,r _y,r _z] ^T(4)

Controller saturation motion model at the state equation of pitch plane is

${\stackrel{\&CenterDot;}{r}}_x=v_x$

${\stackrel{\&CenterDot;}{r}}_y=v_y$

${\stackrel{\&CenterDot;}{v}}_x=a_{px}-a_{ex}$

${\stackrel{\&CenterDot;}{v}}_y=a_{py}-a_{ey}---\left(5\right)$

${\stackrel{\&CenterDot;}{a}}_{px}=-\frac{a_{px}}{\mathrm{\&tau;}}$

${\stackrel{\&CenterDot;}{a}}_{py}=0$

${\stackrel{\&CenterDot;}{d}}_1=0$

Controller saturation motion model at the state equation of the plane of yaw is

${\stackrel{\&CenterDot;}{r}}_x=v_x$

${\stackrel{\&CenterDot;}{r}}_z=v_z$

${\stackrel{\&CenterDot;}{v}}_x=a_{px}-a_{ex}$

${\stackrel{\&CenterDot;}{v}}_z=a_{pz}-a_{ez}---\left(6\right)$

${\stackrel{\&CenterDot;}{a}}_{px}=-\frac{a_{px}}{\mathrm{\&tau;}}$

${\stackrel{\&CenterDot;}{a}}_{pz}=0$

${\stackrel{\&CenterDot;}{d}}_2=0$

In formula, d ₁and d ₂it is the occupied state of controller saturation motion model; Object is to be consistent with the dimension of PN Guidance Law motion model.

Other step and parameter identical with embodiment one.

Embodiment three: present embodiment and embodiment one or two unlike: according to the state equation of the PN Guidance Law motion model of step one at pitch plane in described step 2, PN Guidance Law motion model is at the state equation of the plane of yaw, controller saturation motion model is at the state equation of pitch plane and the controller saturation motion model state equation at the plane of yaw, calculate the state-transition matrix of PN Guidance Law motion model at pitch plane, PN Guidance Law motion model is at the state-transition matrix of the plane of yaw, controller saturation motion model is at the state-transition matrix of pitch plane, controller saturation motion model is at the state-transition matrix of the plane of yaw, detailed process is:

PN Guidance Law motion model at the state-transition matrix of pitch plane is

${\mathrm{\&Phi;}}_{k+1,k}^1=\left[\begin{array}{ccccccc}1& 0& T& 0& 0& 0& 0\\ 0& 1& 0& T& 0& 0& 0\\ 0& 0& 1& 0& T& 0& 0\\ 0& 0& 0& 1& 0& T& 0\\ 0& 0& 0& 0& 1-T/\mathrm{\&tau;}& 0& 0\\ a_1& a_2& a_3& a_4& 0& 1-T/\mathrm{\&tau;}& a_5\\ 0& 0& 0& 0& 0& 0& 1\end{array}\right]---\left(7\right)$

$a_1=-\frac{{\mathrm{TN}}_{\mathrm{\&epsiv;}}(r_x^3{v}_x{v}_y-2r_x^2{r}_y{v}_x^2+2r_x^2{r}_y{v}_y^2-5r_x{r}_y^2{v}_x{v}_y+r_y^3{v}_x^2-r_y^3{v}_y^2)}{{(r_x^2+r_y^2)}^{5/2}\mathrm{\&tau;}}---\left(8\right)$

$a_2=-\frac{{\mathrm{TN}}_{\mathrm{\&epsiv;}}(r_x^3{v}_x^2-r_x^3{v}_y^2+5r_x^2{r}_y{v}_x{v}_y-2r_x{r}_y^2{v}_x^2+2r_x{r}_y^2{v}_y^2-r_y^3{v}_x{v}_y)}{{(r_x^2+r_y^2)}^{5/2}\mathrm{\&tau;}}---\left(9\right)$

$a_3=\frac{{\mathrm{TN}}_{\mathrm{\&epsiv;}}(r_x^2{v}_y-2r_x{r}_y{v}_x-r_y^2{v}_y)}{{(r_x^2+r_y^2)}^{3/2}\mathrm{\&tau;}}---\left(10\right)$

$a_4=\frac{{\mathrm{TN}}_{\mathrm{\&epsiv;}}(r_x^2{v}_x+2r_x{r}_y{v}_y-r_y^2{v}_x)}{{(r_x^2+r_y^2)}^{3/2}\mathrm{\&tau;}}---\left(11\right)$

$a_5=\frac{T(r_x{v}_x+r_y{v}_y)(r_x{v}_y-r_y{v}_x)}{{(r_x^2+r_y^2)}^{3/2}\mathrm{\&tau;}}---\left(12\right)$

In formula, T is measuring period, and unit is second;

PN Guidance Law motion model at the state-transition matrix of the plane of yaw is

${\mathrm{\&Phi;}}_{k+1,k}^2=\left[\begin{array}{ccccccc}1& 0& T& 0& 0& 0& 0\\ 0& 1& 0& T& 0& 0& 0\\ 0& 0& 1& 0& T& 0& 0\\ 0& 0& 0& 1& 0& T& 0\\ 0& 0& 0& 0& 1-T/\mathrm{\&tau;}& 0& 0\\ a_6& a_7& a_8& a_9& 0& 1-T/\mathrm{\&tau;}& a_{10}\\ 0& 0& 0& 0& 0& 0& 1\end{array}\right]---\left(13\right)$

$a_6=-\frac{{\mathrm{TN}}_{\mathrm{\&beta;}}(r_x^3{v}_x{v}_z-2r_x^2{r}_z{v}_x^2+2r_x^2{r}_z{v}_z^2-5r_x{r}_z^2{v}_x{v}_z+r_z^3{v}_x^2-r_z^3{v}_z^2)}{{(r_x^2+r_z^2)}^{5/2}\mathrm{\&tau;}}---\left(14\right)$

$a_7=-\frac{{\mathrm{TN}}_{\mathrm{\&beta;}}(r_x^3{v}_x^2-r_x^3{v}_z^2+5r_x^2{r}_z{v}_x{v}_z-2r_x{r}_z^2{v}_x^2+2r_x{r}_z^2{v}_z^2-r_z^3{v}_x{v}_z)}{{(r_x^2+r_z^2)}^{5/2}\mathrm{\&tau;}}---\left(15\right)$

$a_8=\frac{{\mathrm{TN}}_{\mathrm{\&beta;}}(r_x^2{v}_z-2r_x{r}_z{v}_x-r_z^2{v}_z)}{{(r_x^2+r_z^2)}^{3/2}\mathrm{\&tau;}}---\left(16\right)$

$a_9=\frac{{\mathrm{TN}}_{\mathrm{\&beta;}}(r_x^2{v}_x+2r_x{r}_z{v}_z-r_z^2{v}_x)}{{(r_x^2+r_z^2)}^{3/2}\mathrm{\&tau;}}---\left(17\right)$

$a_{10}=\frac{T(r_x{v}_x+r_z{v}_z)(r_x{v}_z-r_z{v}_x)}{{(r_x^2+r_z^2)}^{3/2}\mathrm{\&tau;}}---\left(18\right)$

Controller saturation motion model at the state-transition matrix of pitch plane is

${\mathrm{\&Phi;}}_{k+1,k}^3=\left[\begin{array}{ccccccc}1& 0& T& 0& 0& 0& 0\\ 0& 1& 0& T& 0& 0& 0\\ 0& 0& 1& 0& T& 0& 0\\ 0& 0& 0& 1& 0& T& 0\\ 0& 0& 0& 0& 1-T/\mathrm{\&tau;}& 0& 0\\ 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 1\end{array}\right]---\left(19\right)$

Controller saturation motion model at the state-transition matrix of the plane of yaw is

${\mathrm{\&Phi;}}_{k+1,k}^4=\left[\begin{array}{ccccccc}1& 0& T& 0& 0& 0& 0\\ 0& 1& 0& T& 0& 0& 0\\ 0& 0& 1& 0& T& 0& 0\\ 0& 0& 0& 1& 0& T& 0\\ 0& 0& 0& 0& 1-T/\mathrm{\&tau;}& 0& 0\\ 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 1\end{array}\right]---\left(20\right).$

Other step and parameter identical with embodiment one or two.

Embodiment four: one of present embodiment and embodiment one to three unlike: according to the state-transition matrix of the PN Guidance Law motion model of step 2 at pitch plane in described step 3, PN Guidance Law motion model is at the state-transition matrix of the plane of yaw, controller saturation motion model is at the state-transition matrix of pitch plane, controller saturation motion model is at the state-transition matrix of the plane of yaw, design PN Guidance Law motion model is at the Kalman filter equation of pitch plane, PN Guidance Law motion model is at the Kalman filter equation of the plane of yaw, controller saturation motion model is at the Kalman filter equation of pitch plane and the controller saturation motion model Kalman filter equation at the plane of yaw, Kalman filter equation is Kalman filter equation, detailed process is:

PN Guidance Law identification wave filter and saturated identification wave filter all adopt EKF, and EKF designs for Kalman filter, according to the state-transition matrix that step 2 calculates,

PN Guidance Law motion model at the Kalman filter equation of pitch plane is

${\hat{X}}_{k+1,k+1}={\hat{X}}_{k+1,k}+K_k^{\left(1\right)}(Z_{k+1}-H{\hat{X}}_{k+1,k})---\left(21\right)$

$K_{k+1}^{\left(1\right)}=P_{k+1,k}^{\left(1\right)}{H}^T{({\mathrm{HP}}_{k+1,k}^{\left(1\right)}{H}^T+R^{\left(1\right)})}^{-1}---\left(22\right)$

$P_{k+1,k}^{\left(1\right)}={\mathrm{\&Phi;}}_{k+1,k}^{\left(1\right)}{P}_{k,k}^{\left(1\right)}{\left({\mathrm{\&Phi;}}_{k+1,k}^{\left(1\right)}\right)}^T+Q^{\left(1\right)}---\left(23\right)$

$P_{k+1,k+1}^{\left(1\right)}=(I-K_{k+1}^{\left(1\right)}H)P_{k+1,k}^{\left(1\right)}---\left(24\right)$

PN Guidance Law motion model at the Kalman filter equation of the plane of yaw is

${\hat{X}}_{k+1,k+1}={\hat{X}}_{k+1,k}+K_k^{\left(2\right)}(Z_{k+1}-H{\hat{X}}_{k+1,k})---\left(25\right)$

$K_{k+1}^{\left(2\right)}=P_{k+1,k}^{\left(2\right)}{H}^T{({\mathrm{HP}}_{k+1,k}^{\left(2\right)}{H}^T+R^{\left(2\right)})}^{-1}---\left(26\right)$

$P_{k+1,k}^{\left(2\right)}={\mathrm{\&Phi;}}_{k+1,k}^{\left(2\right)}{P}_{k,k}^{\left(2\right)}{\left({\mathrm{\&Phi;}}_{k+1,k}^{\left(2\right)}\right)}^T+Q^{\left(2\right)}---\left(27\right)$

$P_{k+1,k+1}^{\left(2\right)}=(I-K_{k+1}^{\left(2\right)}H)P_{k+1,k}^{\left(2\right)}---\left(28\right)$

Controller saturation motion model at the Kalman filter equation of pitch plane is

${\hat{X}}_{k+1,k+1}={\hat{X}}_{k+1,k}+K_k^{\left(3\right)}(Z_{k+1}-H{\hat{X}}_{k+1,k})---\left(29\right)$

$K_{k+1}^{\left(3\right)}=P_{k+1,k}^{\left(3\right)}{H}^T{({\mathrm{HP}}_{k+1,k}^{\left(3\right)}{H}^T+R^{\left(3\right)})}^{-1}---\left(30\right)$

$P_{k+1,k}^{\left(3\right)}={\mathrm{\&Phi;}}_{k+1,k}^{\left(3\right)}{P}_{k,k}^{\left(3\right)}{\left({\mathrm{\&Phi;}}_{k+1,k}^{\left(3\right)}\right)}^T+Q^{\left(3\right)}---\left(31\right)$

$P_{k+1,k+1}^{\left(3\right)}=(I-K_{k+1}^{\left(3\right)}H)P_{k+1,k}^{\left(3\right)}---\left(32\right)$

Controller saturation motion model at the Kalman filter equation of the plane of yaw is

${\hat{X}}_{k+1,k+1}={\hat{X}}_{k+1,k}+K_k^{\left(4\right)}(Z_{k+1}-H{\hat{X}}_{k+1,k})---\left(33\right)$

$K_{k+1}^{\left(4\right)}=P_{k+1,k}^{\left(4\right)}{H}^T{({\mathrm{HP}}_{k+1,k}^{\left(4\right)}{H}^T+R^{\left(4\right)})}^{-1}---\left(34\right)$

$P_{k+1,k}^{\left(4\right)}={\mathrm{\&Phi;}}_{k+1,k}^{\left(4\right)}{P}_{k,k}^{\left(4\right)}{\left({\mathrm{\&Phi;}}_{k+1,k}^{\left(4\right)}\right)}^T+Q^{\left(4\right)}---\left(35\right)$

$P_{k+1,k+1}^{\left(4\right)}=(I-K_{k+1}^{\left(4\right)}H)P_{k+1,k}^{\left(4\right)}---\left(36\right)$

In formula, it is the estimated value to k+1 moment state; it is the predicted value to k+1 moment state;

with the filter gain of PN Guidance Law motion model at pitch plane in k moment and k+1 moment respectively;

with the filter gain of PN Guidance Law motion model at the plane of yaw in k moment and k+1 moment respectively;

with the filter gain of controller saturation motion model at pitch plane in k moment and k+1 moment respectively;

with the filter gain of controller saturation motion model at the plane of yaw in k moment and k+1 moment respectively;

for k+1 moment PN Guidance Law motion model is at the covariance matrix of the state forecast of pitch plane;

for k+1 moment PN Guidance Law motion model is at the covariance matrix of the state forecast of the plane of yaw;

for k+1 moment controller saturation motion model is at the covariance matrix of the state forecast of pitch plane;

for k+1 moment controller saturation motion model is at the covariance matrix of the state forecast of the plane of yaw;

R ⁽¹⁾for PN Guidance Law motion model is at the measurement noises covariance matrix of pitch plane;

R ⁽²⁾for PN Guidance Law motion model is at the measurement noises covariance matrix of the plane of yaw;

R ⁽³⁾for controller saturation motion model is at the measurement noises covariance matrix of pitch plane;

R ⁽⁴⁾for controller saturation motion model is at the measurement noises covariance matrix of the plane of yaw;

for PN Guidance Law motion model is at the state-transition matrix in k moment to the k+1 moment of pitch plane;

for PN Guidance Law motion model is at the state-transition matrix in k moment to the k+1 moment of the plane of yaw;

for controller saturation motion model is at the state-transition matrix in k moment to the k+1 moment of pitch plane;

for controller saturation motion model is at the state-transition matrix in k moment to the k+1 moment of the plane of yaw;

Q ⁽¹⁾for PN Guidance Law motion model is at the process noise covariance matrix of pitch plane;

Q ⁽²⁾for PN Guidance Law motion model is at the process noise covariance matrix of the plane of yaw;

Q ⁽³⁾for controller saturation motion model is at the process noise covariance matrix of pitch plane;

Q ⁽⁴⁾for controller saturation motion model is at the process noise covariance matrix of the plane of yaw;

for PN Guidance Law motion model is at the k moment state estimation covariance matrix of pitch plane;

for PN Guidance Law motion model is at the k moment state estimation covariance matrix of the plane of yaw;

for controller saturation motion model is at the k moment state estimation covariance matrix of pitch plane;

for controller saturation motion model is at the k moment state estimation covariance matrix of the plane of yaw;

for PN Guidance Law motion model is at the k+1 moment state estimation covariance matrix of pitch plane;

for PN Guidance Law motion model is at the k+1 moment state estimation covariance matrix of the plane of yaw;

for controller saturation motion model is at the k+1 moment state estimation covariance matrix of pitch plane;

for controller saturation motion model is at the k+1 moment state estimation covariance matrix of the plane of yaw;

for PN Guidance Law motion model in k moment of pitch plane to the transposition of the state-transition matrix in k+1 moment;

for PN Guidance Law motion model in k moment of the plane of yaw to the transposition of the state-transition matrix in k+1 moment;

for controller saturation motion model in k moment of pitch plane to the transposition of the state-transition matrix in k+1 moment;

for controller saturation motion model in k moment of the plane of yaw to the transposition of the state-transition matrix in k+1 moment;

Z _k+1for the measured value in k+1 moment; I is 7 × 7 unit matrixs; H is calculation matrix, H=[1100000], H ^tfor the transposition of H, k is sampling instant.

Other step and parameter identical with one of embodiment one to three.

Embodiment five: one of present embodiment and embodiment one to four unlike: according to the Kalman filter equation of the PN Guidance Law motion model of step 3 at pitch plane in described step 4, PN Guidance Law motion model is at the Kalman filter equation of the plane of yaw, at the Kalman filter equation of pitch plane and controller saturation motion model, at the Kalman filter equation of the plane of yaw, (2 pitching one group are a MMAE to controller saturation motion model, 2 driftage one group be a MMAE), calculate the posterior probability of pursuer current time PN Guidance Law motion model and the posterior probability of saturating control motion model, detailed process is:

If model set is M={M _j| j=1 ..., r}, in formula, r is the number of model, r=2, M ₁represent PN Guidance Law motion model, M ₂representing controller saturation motion model, is P{M at the prior probability of k moment model j _j| Z ^k-1, be P{M in the posterior probability of k moment model j _j| Z ^k, the measured value in kth moment is z (k), and the measured value sequence definition in a front k moment is

$Z^k\stackrel{\mathrm{\&Delta;}}{=}\left\{z\left(k\right)\right|k=1,\mathrm{...},k\}---\left(37\right)$

According to Bayes' theorem,

$P\left\{M_j\right|Z^k\}=P\{M_j|z\left(k\right),Z^{k-1}\}=\frac{p\left(z\right(k\left)\right|Z^{k-1},M_j)P\left\{M_j\right|Z^{k-1}\}}{\underset{i=1}{\overset{r}{\&Sigma;}}p\left(z\right(k\left)\right|Z^{k-1},M_i)P\left\{M_i\right|Z^{k-1}\}}---\left(38\right)$

In formula, P{M _i| Z ^k-1be the prior probability of model i, and p (z (k) | Z ^k-1, M _i) for model i is at the likelihood function in k moment, P{M _j| z (k), Z ^k-1it is k moment model M _jposterior probability, p (z (k) | Z ^k-1, M _j) for model j is at the likelihood function in k moment; Therefore the posterior probability of model j calculates by the prior probability of all models and likelihood function, and the likelihood function of model i is the distribution function that this model respective filter newly ceases; Assuming that it is the normal distribution of zero that new breath obeys average, if the variance in its k moment is S _i(k), namely

p(z(k)|Z ^k-1,M _i)＝N(0,S _j(k))(39)

In formula, N (0, S _j(k)) to represent average be 0, variance is S _jthe probability density function of (k) normal distribution; S _jk () is the variance that model j newly ceases, namely

$v\left(k\right)=z\left(k\right)-H{\hat{X}}_{k+1,k}---\left(40\right)$

S _j(k)＝E[v(k)v ^T(k)](41)

In formula, v (k) is new breath, v ^tk () is the transposition of v (k), E [] is mathematical expectation function;

If the state estimation of the Kalman filter that model i is corresponding is the then state estimation in MMAE wave filter kth moment for

${\hat{X}}_{k,k}=\underset{i=1}{\overset{r}{\&Sigma;}}P\left\{M_i\right|Z^k\}{\hat{X}}_{k,k}^i---\left(42\right)$

I value is 1 or 2.

Other step and parameter identical with one of embodiment one to four.

Embodiment six: one of present embodiment and embodiment one to five unlike: the posterior probability of pursuer current time PN Guidance Law motion model obtained according to step 4 in described step 5 and the posterior probability of saturated motion model, calculate the moment of PN Guidance Law motion model when pursuer is switched to controller unsaturation from saturated motion model during controller saturation; Adopted controller saturation motion model at the Kalman filter equation of pitch plane and the controller saturation motion model estimated result at the Kalman filter equation of the plane of yaw before the moment of the PN Guidance Law motion model when pursuer is switched to controller unsaturation from saturated motion model during controller saturation; Otherwise adopt PN Guidance Law motion model at the Kalman filter equation of pitch plane and the PN Guidance Law motion model estimated result at the Kalman filter equation of the plane of yaw; Detailed process is:

Consider that the external behavior of pursuer is no longer the motion model of PN Guidance Law due to saturated impact, removes the saturated motion model entering PN Guidance Law afterwards in the first stage of guidance.We adopt two MMAE wave filters to estimate the state of two stage pursuer respectively; These two MMAE wave filters are called MMAE_A and MMAE_B wave filter, and wherein the first stage estimated by MMAE_A wave filter, state when namely pursuer is in controller saturation; And subordinate phase estimated by MMAE_B wave filter, state when namely pursuer is in controller unsaturation; Need the moment of the PN Guidance Law motion model estimated when pursuer is switched to controller unsaturation from saturated motion model during controller saturation in addition, thus determine to adopt the estimated value of which MMAE wave filter to represent the state of current pursuer at certain moment t;

If pursuer does not enter state during controller saturation, then the motion model of pursuer whole process is consistent, and the result that MMAE_A with MMAE_B estimates is consistent;

The method of estimation detailed process in the moment of PN Guidance Law motion model when pursuer is switched to controller unsaturation from saturated motion model during controller saturation is:

All containing two models in each MMAE wave filter: saturated motion model when PN Guidance Law motion model during controller unsaturation and controller saturation; Calculate the posterior probability of present filter two models according to multiple model filtering method, then according to the posterior probability of present filter two models, determine saturation stage when whether pursuer is in controller saturation at present; MMAE_B wave filter adopts the mode of scanning to check saturated exit point, namely at interval of the time that check_period is long, judge the motion model of current pursuer according to formula (48), the posterior probability of saturated motion model when PN Guidance Law motion model when supposing MMAE wave filter middle controller unsaturation and controller saturation is respectively P{M _pNGand P{M _sat, then deterministic process is:

mode＝PNGifP{M _PNG}-P{M _Sat}≥ε(48)

mode＝Saturationelse

In formula, ε is threshold value.Work as P{M _pNG-P{M _satduring }>=ε, then the model that current pursuer adopts is PN Guidance Law model, judges that pursuer has exited state of saturation; Otherwise restart and initialization MMAE_B wave filter.Once judge every the check_period time, until judge that pursuer has exited state of saturation position.Check_period is round of visits, and general check_period is set to 0.5 second.

When ε gets larger value, higher to the judging nicety rate of pursuer "current" model; Otherwise it is when ε gets smaller value, lower to the judging nicety rate of pursuer "current" model.But this does not also mean that the value of ε is the bigger the better, after can analyze for this.

Check_period is relevant to parameter ε ε gets very high values, and the check_period of needs is also comparatively large, and it is inaccurate that this will cause pursuer to exit saturated time point estimation; When ε gets small value, check_period can get less value.This can improve the precision that pursuer exits saturation time point estimation outwardly, but can reduce the accuracy of judgement degree of pursuer "current" model when ε gets small value.Therefore need select suitable ε and check_period value just can take into account pursuer "current" model accuracy of judgement degree and exit the degree of accuracy of estimation of saturated moment.

Here is that MMAE_B exits the false code of saturated moment assessment function to pursuer, and algorithm calculates the backed off after random that pursuer exits the estimated value quit_saturation_time in saturated moment.

Other step and parameter identical with one of embodiment one to five.

Following examples are adopted to verify beneficial effect of the present invention:

Embodiment one:

A kind of proportional navigation law identification of the present embodiment filtering method is specifically prepared according to following steps:

The multi-model PN Guidance Law identification wave filter that this patent proposes has carried out simulating, verifying in a missile breakthrough scene.In order to verify the performance of this Guidance Law identification wave filter, the estimated performance of itself and PN Guidance Law Kalman identification wave filter is contrasted.Controller saturation situation do not considered by this wave filter, thinks that pursuer is omnidistance and moves under the constraint of PN Guidance Law.

In this missile breakthrough scene, pursuer tackles evader.Consider that PN Guidance Law is common Guidance Law, assuming that evader knows that pursuer adopts PN Guidance Law in advance, but do not know its navigation constant in pitching and jaw channel.Evader can obtain the position of pursuer, uses multi-model PN Guidance Law identification wave filter to estimate the PN Guidance Law navigation constant of pursuer and the acceleration on pitch plane and the plane of yaw.These information can pass to the defending missile of protection evader, and defending missile utilizes these information to tackle pursuer, to protect evader.The simulation configurations of pursuer describes in Table 1, and the simulation configurations of evader describes in table 2.

Table 1:pursuer simulation configurations information

Table 2:evader simulation configurations information

Process noise matrix and the measurement noises arranged in matrix of PN Guidance Law Kalman identification wave filter and saturated Kalman identification wave filter are

Q _PNG＝Q _Sat＝diag(10,10,10,20,20,20,2)

R _PNG＝R _Sat＝diag(10,10,10)

Fig. 2-15 is the simulation results under this simulating scenes.Fig. 2 and Fig. 3 is that in the motor-driven situation of evader constant value, wave filter MMAE_A is to the estimated value of the navigation constant of pursuer Guidance Law.MMAE_A and MMAE_B exits in saturated rear very short time at pursuer the estimated value of navigation constant and has converged to true value.It is 2.7s and 4.6s respectively that Pursuer exits the saturated time at pitch channel and jaw channel controller.Can see when pursuer controller is in state of saturation from Fig. 2 and Fig. 3, PN Guidance Law Kalman identification wave filter causes not converging to true value to the estimated value of navigation constant due to large model error.Pursuer exit saturated after the very fast estimated value to navigation constant converged to true value, be pitch channel navigation constant be respectively 4, jaw channel navigation constant is 5.

Above-mentioned simulation process considers pursuer controller in guidance process and produces saturated possibility in the early stage, employ multiple model filtering method and identification has been carried out to pursuer Guidance Law, estimate to the navigation constant under PN Guidance Law with at pitch plane and plane of yaw acceleration.Identification can be carried out more accurately to Guidance Law.In an embodiment of the present invention, if MMAE wave filter is Var to the square error that pursuer acceleration estimation produces _mMAE; Omnidistance use PN Guidance Law is Var to the square error that pursuer acceleration estimation produces _pN.If α=Var _mMAE/ Var _pN.As 0< α <1, to pursuer acceleration estimation, MMAE wave filter comparatively illustrates that the omnidistance PN Guidance Law that uses is more accurate to pursuer acceleration estimation; α is less, illustrates to use PN Guidance Law concerning pursuer acceleration estimation relative to whole process, and MMAE is more accurate to pursuer acceleration estimation.

When evader is without time motor-driven, at pursuer saturation stage, to pursuer acceleration a _yestimation time α=0.59; To pursuer acceleration a _zestimation time α=0.56.

When evader carry out constant value motor-driven time, at pursuer saturation stage, to pursuer acceleration a _yestimation time α=0.55; To pursuer acceleration a _zestimation time α=0.62.

When evader carry out sinusoidal motor-driven time, at pursuer saturation stage, to pursuer acceleration a _yestimation time α=0.56; To pursuer acceleration a _zestimation time α=0.60.

Fig. 4, Fig. 5, Fig. 6, Fig. 7 are under the inorganic emotionally condition of evader, the estimation of transshipping at pitch channel and jaw channel to pursuer when whole guidance process all uses PN Guidance Law Kalman identification wave filter to carry out identification to pursuer.When pursuer controller is in state of saturation, the overload estimated value of PN Guidance Law Kalman identification wave filter in pitching and jaw channel all has larger difference with true value.And once pursuer controller exits saturated, the estimated value of PN Guidance Law Kalman identification wave filter has converged to true value soon.Multiple-model estimator method in this paper considers the saturated impact on acceleration identification.Consider the estimated value of MMAE_A and omnidistance PN Guidance Law identification Kalman filter as shown in Figure 4.The saturated time is exited according to the pursuer calculated, the controller of pursuer exit saturated before use the estimated value of MMAE_A wave filter, and the controller of pursuer exit saturated after use the estimated result of omnidistance PN Guidance Law identification Kalman filter.

Comparison diagram 4, Fig. 5, Fig. 6, Fig. 7 to be seen at saturation stage wave filter in this paper the overload estimated accuracy of pursuer pitch channel and jaw channel all higher than the estimated accuracy of PN Guidance Law Kalman identification wave filter.And pursuer controller exit saturated after, the estimated accuracy of two kinds of wave filters is identical.When evader perform constant value motor-driven or sinusoidal motor-driven time, the not motor-driven result of result and evader is consistent.As shown in Fig. 8, Fig. 9, Figure 10, Figure 11, Figure 12, Figure 13, Figure 14, Figure 15.

The present invention also can have other various embodiments; when not deviating from the present invention's spirit and essence thereof; those skilled in the art are when making various corresponding change and distortion according to the present invention, but these change accordingly and are out of shape the protection domain that all should belong to the claim appended by the present invention.

Claims (6)

1. a proportional navigation law identification filtering method, is characterized in that what a kind of proportional navigation law identification filtering method specifically carried out according to following steps:

Step one: pursuer has two kinds of motion models: PN Guidance Law motion model when saturated motion model during controller saturation and controller unsaturation; Set up PN Guidance Law motion model at the state equation of pitch plane, PN Guidance Law motion model at the state equation of the plane of yaw, controller saturation motion model at the state equation of pitch plane and the controller saturation motion model state equation at the plane of yaw; Evader is ballistic missile; Pursuer is interception guided missile; PN Guidance Law is proportional navigation law;

Step 2: according to the PN Guidance Law motion model of step one state equation at pitch plane, PN Guidance Law motion model is at the state equation of the plane of yaw, controller saturation motion model is at the state equation of pitch plane and the controller saturation motion model state equation at the plane of yaw, calculate the state-transition matrix of PN Guidance Law motion model at pitch plane, PN Guidance Law motion model is at the state-transition matrix of the plane of yaw, controller saturation motion model is at the state-transition matrix of pitch plane, controller saturation motion model is at the state-transition matrix of the plane of yaw,

Step 3: according to the PN Guidance Law motion model of the step 2 state-transition matrix at pitch plane, PN Guidance Law motion model is at the state-transition matrix of the plane of yaw, controller saturation motion model is at the state-transition matrix of pitch plane, controller saturation motion model is at the state-transition matrix of the plane of yaw, design PN Guidance Law motion model is at the Kalman filter equation of pitch plane, PN Guidance Law motion model is at the Kalman filter equation of the plane of yaw, controller saturation motion model is at the Kalman filter equation of pitch plane and the controller saturation motion model Kalman filter equation at the plane of yaw, Kalman filter equation is Kalman filter equation,

Step 4: at the Kalman filter equation of the plane of yaw, calculate the posterior probability of pursuer current time PN Guidance Law motion model and the posterior probability of controller saturation motion model at the Kalman filter equation of pitch plane and controller saturation motion model at the Kalman filter equation of the plane of yaw, controller saturation motion model at the Kalman filter equation of pitch plane, PN Guidance Law motion model according to the PN Guidance Law motion model of step 3;

Step 5: the posterior probability of pursuer current time PN Guidance Law motion model obtained according to step 4 and the posterior probability of saturated motion model, calculates the moment of PN Guidance Law motion model when pursuer is switched to controller unsaturation from saturated motion model during controller saturation; Adopted controller saturation motion model at the Kalman filter equation of pitch plane and the controller saturation motion model estimated result at the Kalman filter equation of the plane of yaw before the moment of the PN Guidance Law motion model when pursuer is switched to controller unsaturation from saturated motion model during controller saturation; Otherwise adopt PN Guidance Law motion model at the Kalman filter equation of pitch plane and the PN Guidance Law motion model estimated result at the Kalman filter equation of the plane of yaw.

2. a kind of proportional navigation law identification filtering method according to claim 1, is characterized in that: in described step one, pursuer has two kinds of motion models: PN Guidance Law motion model when saturated motion model during controller saturation and controller unsaturation; Set up PN Guidance Law motion model at the state equation of pitch plane, PN Guidance Law motion model at the state equation of the plane of yaw, controller saturation motion model at the state equation of pitch plane and the controller saturation motion model state equation at the plane of yaw, pursuer is interception guided missile; Detailed process is:

In missile breakthrough scene, prominent anti-ballistic missile will escape the interception of the interception guided missile that the other side's missile defense systems is launched, it is escape side, be called evader, interception guided missile is the side of chasing, be called pursuer, pursuer has two kinds of motion models: PN Guidance Law motion model when saturated motion model during controller saturation and controller unsaturation;

The PN Guidance Law equation of motion modeling of pursuer, r _xfor the component of evader and pursuer relative position under evader inertial system x-axis; r _yfor the component of evader and pursuer relative position under evader inertial system y-axis; r _zfor the component of evader and pursuer relative position under evader inertial system z-axis; Evader is ballistic missile, and pursuer is interception guided missile;

r _x＝x _p-x _er _y＝y _p-y _er _z＝z _p-z _e(1)

In formula, [x _p, y _p, z _p] ^t[x _e, y _e, z _e] ^tbe the position of pursuer and evader under evader inertial system respectively, T is transposition, x _pfor the position of pursuer under evader inertial system x-axis, y _pfor the position of pursuer under evader inertial system y-axis, z _pfor the position of pursuer under evader inertial system z-axis, x _efor the position of evader under evader inertial system x-axis, y _efor the position of evader under evader inertial system y-axis, z _efor the position of evader under evader inertial system z-axis;

PN Guidance Law motion model at the state equation of pitch plane is

${\stackrel{\&CenterDot;}{r}}_x=v_x$

${\stackrel{\&CenterDot;}{r}}_y=v_y$

${\stackrel{\&CenterDot;}{v}}_x=a_{px}-a_{ex}$

${\stackrel{\&CenterDot;}{v}}_y=a_{py}-a_{ey}---\left(2\right)$

${\stackrel{\&CenterDot;}{a}}_{px}=-\frac{a_{px}}{\mathrm{\&tau;}}$

${\stackrel{\&CenterDot;}{a}}_{py}=N_{\mathrm{\&epsiv;}}\frac{(r_x{v}_x+r_y{v}_y)(r_x{v}_y-r_y{v}_x)}{{(r_x^2+r_y^2)}^{\frac{3}{2}}\mathrm{\&tau;}}-\frac{a_{py}}{\mathrm{\&tau;}}$

${\stackrel{\&CenterDot;}{N}}_{\mathrm{\&epsiv;}}=0$

PN Guidance Law motion model at the state equation of the plane of yaw is

${\stackrel{\&CenterDot;}{r}}_x=v_x$

${\stackrel{\&CenterDot;}{r}}_z=v_z$

${\stackrel{\&CenterDot;}{v}}_x=a_{px}-a_{ex}$

${\stackrel{\&CenterDot;}{v}}_z=a_{pz}-a_{ez}---\left(3\right)$

${\stackrel{\&CenterDot;}{a}}_{px}=-\frac{a_{px}}{\mathrm{\&tau;}}$

${\stackrel{\&CenterDot;}{a}}_{pz}=-\frac{N_{\mathrm{\&beta;}}}{\mathrm{\&tau;}}\frac{(r_x{v}_x+r_z{v}_z)(r_z{v}_x-r_x{v}_z)}{{(r_x^2+r_z^2)}^{\frac{3}{2}}}-\frac{a_{pz}}{\mathrm{\&tau;}}$

${\stackrel{\&CenterDot;}{N}}_{\mathrm{\&beta;}}=0$

In formula, v _xfor the component of evader and pursuer relative velocity under evader inertial system x-axis, v _yfor the component of evader and pursuer relative velocity under evader inertial system y-axis, v _zfor the component of evader and pursuer relative velocity under evader inertial system z-axis; a _pxfor the component of acceleration under evader inertial system x-axis of pursuer, a _pyfor the component of acceleration under evader inertial system y-axis of pursuer, a _pzfor the component of acceleration under evader inertial system z-axis of pursuer; a _ex, a _eyand a _ezthe component of acceleration under evader inertial system of pursuer, N _εthe navigation constant of pursuer at pitch plane; N _βthe navigation constant of pursuer at the plane of yaw; τ is time constant, for r _xfirst order derivative, for r _yfirst order derivative, for r _zfirst order derivative, for v _xfirst order derivative, for v _yfirst order derivative, for v _zfirst order derivative, for a _pxfirst order derivative, for a _pyfirst order derivative, for a _pzfirst order derivative, for N _εfirst order derivative, for N _βfirst order derivative;

Suppose that evader obtains the position of pursuer under inertial system, and then calculate relatively self position, the measurement equation calculating the position of self is relatively

h＝[r _x,r _y,r _z] ^T(4)

Controller saturation motion model at the state equation of pitch plane is

${\stackrel{\&CenterDot;}{r}}_x=v_x$

${\stackrel{\&CenterDot;}{r}}_y=v_y$

${\stackrel{\&CenterDot;}{v}}_x=a_{px}-a_{ex}$

${\stackrel{\&CenterDot;}{v}}_y=a_{py}-a_{ey}---\left(5\right)$

${\stackrel{\&CenterDot;}{a}}_{px}=-\frac{a_{px}}{\mathrm{\&tau;}}$

${\stackrel{\&CenterDot;}{a}}_{py}=0$

${\stackrel{\&CenterDot;}{d}}_1=0$

Controller saturation motion model at the state equation of the plane of yaw is

${\stackrel{\&CenterDot;}{r}}_x=v_x$

${\stackrel{\&CenterDot;}{r}}_z=v_z$

${\stackrel{\&CenterDot;}{v}}_x=a_{px}-a_{ex}$

${\stackrel{\&CenterDot;}{v}}_z=a_{pz}-a_{ez}---\left(6\right)$

${\stackrel{\&CenterDot;}{a}}_{px}=-\frac{a_{px}}{\mathrm{\&tau;}}$

${\stackrel{\&CenterDot;}{a}}_{pz}=0$

${\stackrel{\&CenterDot;}{d}}_2=0$

In formula, d ₁and d ₂it is the occupied state of controller saturation motion model.

3. a kind of proportional navigation law identification filtering method according to claim 2, it is characterized in that: according to the state equation of the PN Guidance Law motion model of step one at pitch plane in described step 2, PN Guidance Law motion model is at the state equation of the plane of yaw, controller saturation motion model is at the state equation of pitch plane and the controller saturation motion model state equation at the plane of yaw, calculate the state-transition matrix of PN Guidance Law motion model at pitch plane, PN Guidance Law motion model is at the state-transition matrix of the plane of yaw, controller saturation motion model is at the state-transition matrix of pitch plane, controller saturation motion model is at the state-transition matrix of the plane of yaw, detailed process is:

PN Guidance Law motion model at the state-transition matrix of pitch plane is

${\mathrm{\&Phi;}}_{k+1,k}^1=\left[\begin{array}{ccccccc}1& 0& T& 0& 0& 0& 0\\ 0& 1& 0& T& 0& 0& 0\\ 0& 0& 1& 0& T& 0& 0\\ 0& 0& 0& 1& 0& T& 0\\ 0& 0& 0& 0& 1-T/\mathrm{\&tau;}& 0& 0\\ a_1& a_2& a_3& a_4& 0& 1-T/\mathrm{\&tau;}& a_5\\ 0& 0& 0& 0& 0& 0& 1\end{array}\right]---\left(7\right)$

$a_1=-\frac{{\mathrm{TN}}_{\mathrm{\&epsiv;}}(r_x^3{v}_x{v}_y-2r_x^2{r}_y{v}_x^2+2r_x^2{r}_y{v}_y^2-5r_x{r}_y^2{v}_x{v}_y+r_y^3{v}_x^2-r_y^3{v}_y^2)}{{(r_x^2+r_y^2)}^{5/2}\mathrm{\&tau;}}---\left(8\right)$

$a_2=-\frac{{\mathrm{TN}}_{\mathrm{\&epsiv;}}(r_x^3{v}_x^2-r_x^3{v}_y^2+5r_x^2{r}_y{v}_x{v}_y-2r_x{r}_y^2{v}_x^2+2r_x{r}_y^2{v}_y^2-r_y^3{v}_x{v}_y)}{{(r_x^2+r_y^2)}^{5/2}\mathrm{\&tau;}}---\left(9\right)$

$a_3=\frac{{\mathrm{TN}}_{\mathrm{\&epsiv;}}(r_x^2{v}_y-2r_x{r}_y{v}_x-r_y^2{v}_y)}{{(r_x^2+r_y^2)}^{3/2}\mathrm{\&tau;}}---\left(10\right)$

$a_4=\frac{{\mathrm{TN}}_{\mathrm{\&epsiv;}}(r_x^2{v}_x+2r_x{r}_y{v}_y-r_y^2{v}_x)}{{(r_x^2+r_y^2)}^{3/2}\mathrm{\&tau;}}---\left(11\right)$

$a_5=\frac{T(r_x{v}_x+r_y{v}_y)(r_x{v}_y-r_y{v}_x)}{{(r_x^2+r_y^2)}^{3/2}\mathrm{\&tau;}}---\left(12\right)$

In formula, T is measuring period, and unit is second;

PN Guidance Law motion model at the state-transition matrix of the plane of yaw is

${\mathrm{\&Phi;}}_{k+1,k}^2=\left[\begin{array}{ccccccc}1& 0& T& 0& 0& 0& 0\\ 0& 1& 0& T& 0& 0& 0\\ 0& 0& 1& 0& T& 0& 0\\ 0& 0& 0& 1& 0& T& 0\\ 0& 0& 0& 0& 1-T/\mathrm{\&tau;}& 0& 0\\ a_6& a_7& a_8& a_9& 0& 1-T/\mathrm{\&tau;}& a_{10}\\ 0& 0& 0& 0& 0& 0& 1\end{array}\right]---\left(13\right)$

$a_6=-\frac{{\mathrm{TN}}_{\mathrm{\&beta;}}(r_x^3{v}_x{v}_z-2r_x^2{r}_z{v}_x^2+2r_x^2{r}_z{v}_z^2-5r_x{r}_z^2{v}_x{v}_z+r_z^3{v}_x^2-r_z^3{v}_z^2)}{{(r_x^2+r_z^2)}^{5/2}\mathrm{\&tau;}}---\left(14\right)$

$a_7=-\frac{{\mathrm{TN}}_{\mathrm{\&beta;}}(r_x^3{v}_x^2-r_x^3{v}_z^2+5r_x^2{r}_z{v}_x{v}_z-2r_x{r}_z^2{v}_x^2+2r_x{r}_z^2{v}_z^2-r_z^3{v}_x{v}_z)}{{(r_x^2+r_z^2)}^{5/2}\mathrm{\&tau;}}---\left(15\right)$

$a_8=\frac{{\mathrm{TN}}_{\mathrm{\&beta;}}(r_x^2{v}_z-2r_x{r}_z{v}_x-r_z^2{v}_z)}{{(r_x^2+r_z^2)}^{3/2}\mathrm{\&tau;}}---\left(16\right)$

$a_9=\frac{{\mathrm{TN}}_{\mathrm{\&beta;}}(r_x^2{v}_x+2r_x{r}_z{v}_z-r_z^2{v}_x)}{{(r_x^2+r_z^2)}^{3/2}\mathrm{\&tau;}}---\left(17\right)$

$a_{10}=\frac{T(r_x{v}_x+r_z{v}_z)(r_x{v}_z-r_z{v}_x)}{{(r_x^2+r_z^2)}^{3/2}\mathrm{\&tau;}}---\left(18\right)$

Controller saturation motion model at the state-transition matrix of pitch plane is

${\mathrm{\&Phi;}}_{k+1,k}^3=\left[\begin{array}{ccccccc}1& 0& T& 0& 0& 0& 0\\ 0& 1& 0& T& 0& 0& 0\\ 0& 0& 1& 0& T& 0& 0\\ 0& 0& 0& 1& 0& T& 0\\ 0& 0& 0& 0& 1-T/\mathrm{\&tau;}& 0& 0\\ 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 1\end{array}\right]---\left(19\right)$

Controller saturation motion model at the state-transition matrix of the plane of yaw is

${\mathrm{\&Phi;}}_{k+1,k}^4=\left[\begin{array}{ccccccc}1& 0& T& 0& 0& 0& 0\\ 0& 1& 0& T& 0& 0& 0\\ 0& 0& 1& 0& T& 0& 0\\ 0& 0& 0& 1& 0& T& 0\\ 0& 0& 0& 0& 1-T/\mathrm{\&tau;}& 0& 0\\ 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 1\end{array}\right]---\left(20\right).$

4. a kind of proportional navigation law identification filtering method according to claim 3, it is characterized in that: according to the state-transition matrix of the PN Guidance Law motion model of step 2 at pitch plane in described step 3, PN Guidance Law motion model is at the state-transition matrix of the plane of yaw, controller saturation motion model is at the state-transition matrix of pitch plane, controller saturation motion model is at the state-transition matrix of the plane of yaw, design PN Guidance Law motion model is at the Kalman filter equation of pitch plane, PN Guidance Law motion model is at the Kalman filter equation of the plane of yaw, controller saturation motion model is at the Kalman filter equation of pitch plane and the controller saturation motion model Kalman filter equation at the plane of yaw, Kalman filter equation is Kalman filter equation, detailed process is:

PN Guidance Law identification wave filter and saturated identification wave filter all adopt EKF, and EKF designs for Kalman filter, according to the state-transition matrix that step 2 calculates,

PN Guidance Law motion model at the Kalman filter equation of pitch plane is

${\hat{X}}_{k+1,k+1}={\hat{X}}_{k+1,k}+K_k^{\left(1\right)}(Z_{k+1}-H{\hat{X}}_{k+1,k})---\left(21\right)$

$K_{k+1}^{\left(1\right)}=P_{k+1,k}^{\left(1\right)}{H}^T{({\mathrm{HP}}_{k+1,k}^{\left(1\right)}{H}^T+R^{\left(1\right)})}^{-1}---\left(22\right)$

$P_{k+1,k}^{\left(1\right)}={\mathrm{\&Phi;}}_{k+1,k}^{\left(1\right)}{P}_{k,k}^{\left(1\right)}{\left({\mathrm{\&Phi;}}_{k+1,k}^{\left(1\right)}\right)}^T+Q^{\left(1\right)}---\left(23\right)$

$P_{k+1,k+1}^{\left(1\right)}=(I-K_{k+1}^{\left(1\right)}H)P_{k+1,k}^{\left(1\right)}---\left(24\right)$

PN Guidance Law motion model at the Kalman filter equation of the plane of yaw is

${\hat{X}}_{k+1,k+1}={\hat{X}}_{k+1,k}+K_k^{\left(2\right)}(Z_{k+1}-H{\hat{X}}_{k+1,k})---\left(25\right)$

$K_{k+1}^{\left(2\right)}=P_{k+1,k}^{\left(2\right)}{H}^T{({\mathrm{HP}}_{k+1,k}^{\left(2\right)}{H}^T+R^{\left(2\right)})}^{-1}---\left(26\right)$

$P_{k+1,k}^{\left(2\right)}={\mathrm{\&Phi;}}_{k+1,k}^{\left(2\right)}{P}_{k,k}^{\left(2\right)}{\left({\mathrm{\&Phi;}}_{k+1,k}^{\left(2\right)}\right)}^T+Q^{\left(2\right)}---\left(27\right)$

$P_{k+1,k+1}^{\left(2\right)}=(I-K_{k+1}^{\left(2\right)}H)P_{k+1,k}^{\left(2\right)}---\left(28\right)$

Controller saturation motion model at the Kalman filter equation of pitch plane is

${\hat{X}}_{k+1,k+1}={\hat{X}}_{k+1,k}+K_k^{\left(3\right)}(Z_{k+1}-H{\hat{X}}_{k+1,k})---\left(29\right)$

$K_{k+1}^{\left(3\right)}=P_{k+1,k}^{\left(3\right)}{H}^T{({\mathrm{HP}}_{k+1,k}^{\left(3\right)}{H}^T+R^{\left(3\right)})}^{-1}---\left(30\right)$

$P_{k+1,k}^{\left(3\right)}={\mathrm{\&Phi;}}_{k+1,k}^{\left(3\right)}{P}_{k,k}^{\left(3\right)}{\left({\mathrm{\&Phi;}}_{k+1,k}^{\left(3\right)}\right)}^T+Q^{\left(3\right)}---\left(31\right)$

$P_{k+1,k+1}^{\left(3\right)}=(I-K_{k+1}^{\left(3\right)}H)P_{k+1,k}^{\left(3\right)}---\left(32\right)$

Controller saturation motion model at the Kalman filter equation of the plane of yaw is

${\hat{X}}_{k+1,k+1}={\hat{X}}_{k+1,k}+K_k^{\left(4\right)}(Z_{k+1}-H{\hat{X}}_{k+1,k})---\left(33\right)$

$K_{k+1}^{\left(4\right)}=P_{k+1,k}^{\left(4\right)}{H}^T{({\mathrm{HP}}_{k+1,k}^{\left(4\right)}{H}^T+R^{\left(4\right)})}^{-1}---\left(34\right)$

$P_{k+1,k}^{\left(4\right)}={\mathrm{\&Phi;}}_{k+1,k}^{\left(4\right)}{P}_{k,k}^{\left(4\right)}{\left({\mathrm{\&Phi;}}_{k+1,k}^{\left(4\right)}\right)}^T+Q^{\left(4\right)}---\left(35\right)$

$P_{k+1,k+1}^{\left(4\right)}=(I-K_{k+1}^{\left(4\right)}H)P_{k+1,k}^{\left(4\right)}---\left(36\right)$

In formula, it is the estimated value to k+1 moment state; it is the predicted value to k+1 moment state;

with the filter gain of PN Guidance Law motion model at pitch plane in k moment and k+1 moment respectively;

with the filter gain of PN Guidance Law motion model at the plane of yaw in k moment and k+1 moment respectively;

with the filter gain of controller saturation motion model at pitch plane in k moment and k+1 moment respectively;

with the filter gain of controller saturation motion model at the plane of yaw in k moment and k+1 moment respectively;

for k+1 moment PN Guidance Law motion model is at the covariance matrix of the state forecast of pitch plane;

for k+1 moment PN Guidance Law motion model is at the covariance matrix of the state forecast of the plane of yaw;

for k+1 moment controller saturation motion model is at the covariance matrix of the state forecast of pitch plane;

for k+1 moment controller saturation motion model is at the covariance matrix of the state forecast of the plane of yaw;

R ⁽¹⁾for PN Guidance Law motion model is at the measurement noises matrix of pitch plane;

R ⁽²⁾for PN Guidance Law motion model is at the measurement noises matrix of the plane of yaw;

R ⁽³⁾for controller saturation motion model is at the measurement noises matrix of pitch plane;

R ⁽⁴⁾for controller saturation motion model is at the measurement noises matrix of the plane of yaw;

for PN Guidance Law motion model is at the state-transition matrix in k moment to the k+1 moment of pitch plane;

for PN Guidance Law motion model is at the state-transition matrix in k moment to the k+1 moment of the plane of yaw;

for controller saturation motion model is at the state-transition matrix in k moment to the k+1 moment of pitch plane;

for controller saturation motion model is at the state-transition matrix in k moment to the k+1 moment of the plane of yaw;

Q ⁽¹⁾for PN Guidance Law motion model is at the process noise matrix of pitch plane;

Q ⁽²⁾for PN Guidance Law motion model is at the process noise matrix of the plane of yaw;

Q ⁽³⁾for controller saturation motion model is at the process noise matrix of pitch plane;

Q ⁽⁴⁾for controller saturation motion model is at the process noise matrix of the plane of yaw;

for PN Guidance Law motion model is at the k moment state estimation covariance matrix of pitch plane;

for PN Guidance Law motion model is at the k moment state estimation covariance matrix of the plane of yaw;

for controller saturation motion model is at the k moment state estimation covariance matrix of pitch plane;

for controller saturation motion model is at the k moment state estimation covariance matrix of the plane of yaw;

for PN Guidance Law motion model is at the k+1 moment state estimation covariance matrix of pitch plane;

for PN Guidance Law motion model is at the k+1 moment state estimation covariance matrix of the plane of yaw;

for controller saturation motion model is at the k+1 moment state estimation covariance matrix of pitch plane;

for controller saturation motion model is at the k+1 moment state estimation covariance matrix of the plane of yaw;

for PN Guidance Law motion model in k moment of pitch plane to the transposition of the state-transition matrix in k+1 moment;

for PN Guidance Law motion model in k moment of the plane of yaw to the transposition of the state-transition matrix in k+1 moment;

for controller saturation motion model in k moment of pitch plane to the transposition of the state-transition matrix in k+1 moment;

for controller saturation motion model in k moment of the plane of yaw to the transposition of the state-transition matrix in k+1 moment;

Z _k+1for the measured value in k+1 moment; I is 7 × 7 unit matrixs; H is calculation matrix, H=[1100000], H ^tfor the transposition of H, k is sampling instant.

5. a kind of proportional navigation law identification filtering method according to claim 4, it is characterized in that: according to the Kalman filter equation of the PN Guidance Law motion model of step 3 at pitch plane in described step 4, PN Guidance Law motion model is at the Kalman filter equation of the plane of yaw, controller saturation motion model is at the Kalman filter equation of pitch plane and the controller saturation motion model Kalman filter equation at the plane of yaw, calculate the posterior probability of pursuer current time PN Guidance Law motion model and the posterior probability of saturating control motion model, detailed process is:

If model set is M={M _j| j=1 ..., r}, in formula, r is the number of model, r=2, M ₁represent PN Guidance Law motion model, M ₂representing controller saturation motion model, is P{M at the prior probability of k moment model j _j| Z ^k-1, be P{M in the posterior probability of k moment model j _j| Z ^k, the measured value in kth moment is z (k), and the measured value sequence definition in a front k moment is

$Z^k\stackrel{\mathrm{\&Delta;}}{=}\left\{z\left(k\right)\right|k=1,...,k\}---\left(37\right)$

According to Bayes' theorem,

$P\left\{M_j\right|Z^k\}=P\{M_j|z\left(k\right),Z^{k-1}\}=\frac{p\left(z\right(k),Z^{k-1},M_j)P\left\{M_j\right|Z^{k-1}\}}{\underset{i=1}{\overset{r}{\mathrm{\&Sigma;}}}p\left(z\right(k),Z^{k-1},M_j)P\left\{M_j\right|Z^{k-1}\}}---\left(38\right)$

In formula, P{M _i| Z ^k-1be the prior probability of model i, and p (z (k) | Z ^k-1, M _i) for model i is at the likelihood function in k moment, P{M _j| z (k), Z ^k-1it is k moment model M _jposterior probability, p (z (k) | Z ^k-1, M _j) for model j is at the likelihood function in k moment; Therefore the posterior probability of model j calculates by the prior probability of all models and likelihood function, and the likelihood function of model i is the distribution function that this model respective filter newly ceases; Assuming that it is the normal distribution of zero that new breath obeys average, if the variance in its k moment is S _i(k), namely

p(z(k)|Z ^k-1,M _i)＝N(0,S _j(k))(39)

In formula, N (0, S _j(k)) to represent average be 0, variance is S _jthe probability density function of (k) normal distribution; S _jk () is the variance that model j newly ceases, namely

$v\left(k\right)=z\left(k\right)-H{\hat{X}}_{k+1,k}---\left(40\right)$

S _j(k)＝E[v(k)v ^T(k)](41)

In formula, v (k) is new breath, v ^tk () is the transposition of v (k), E [] is mathematical expectation function.

6. a kind of proportional navigation law identification filtering method according to claim 5, it is characterized in that: the posterior probability of pursuer current time PN Guidance Law motion model obtained according to step 4 in described step 5 and the posterior probability of saturated motion model, calculate the moment of PN Guidance Law motion model when pursuer is switched to controller unsaturation from saturated motion model during controller saturation; Adopted controller saturation motion model at the Kalman filter equation of pitch plane and the controller saturation motion model estimated result at the Kalman filter equation of the plane of yaw before the moment of the PN Guidance Law motion model when pursuer is switched to controller unsaturation from saturated motion model during controller saturation; Otherwise adopt PN Guidance Law motion model at the Kalman filter equation of pitch plane and the PN Guidance Law motion model estimated result at the Kalman filter equation of the plane of yaw.Detailed process is:

Two MMAE wave filters are adopted to estimate the state of two stage pursuer respectively; These two MMAE wave filters are called MMAE_A and MMAE_B wave filter, and wherein the first stage estimated by MMAE_A wave filter, state when namely pursuer is in controller saturation; And subordinate phase estimated by MMAE_B wave filter, state when namely pursuer is in controller unsaturation; Need the moment of the PN Guidance Law motion model estimated when pursuer is switched to controller unsaturation from saturated motion model during controller saturation in addition, thus determine to adopt the estimated value of which MMAE wave filter to represent the state of current pursuer at certain moment t;

If pursuer does not enter state during controller saturation, then the motion model of pursuer whole process is consistent, and the result that MMAE_A with MMAE_B estimates is consistent;

The method of estimation detailed process in the moment of PN Guidance Law motion model when pursuer is switched to controller unsaturation from saturated motion model during controller saturation is:

All containing two models in each MMAE wave filter: saturated motion model when PN Guidance Law motion model during controller unsaturation and controller saturation; Calculate the posterior probability of present filter two models according to multiple model filtering method, then according to the posterior probability of present filter two models, determine saturation stage when whether pursuer is in controller saturation at present; MMAE_B wave filter adopts the mode of scanning to check saturated exit point, namely at interval of the time that check_period is long, check_period is round of visits, judge the motion model of current pursuer according to formula (48), the posterior probability of saturated motion model when PN Guidance Law motion model when supposing MMAE wave filter middle controller unsaturation and controller saturation is respectively P{M _pNGand P{M _sat, then deterministic process is:

$\begin{array}{cc}\mathrm{mod}e=PNG& ifP\left\{M_{PNG}\right\}-P\left\{M_{Sat}\right\}\&GreaterEqual;\mathrm{\&epsiv;}\\ \mathrm{mod}e=Saturation& else\end{array}---\left(48\right)$

In formula, ε is threshold value; Work as P{M _pNG-P{M _satduring }>=ε, then the model that current pursuer adopts is PN Guidance Law model, judges that pursuer has exited state of saturation; Otherwise restart and initialization MMAE_B wave filter, check_period is set to 0.5 second.