Ballistic Prediction Method for Hypersonic Vehicles in Near-Space Based on Maneuver Pattern Recognition

Abstract

A method for predicting the trajectory of a near-space hypersonic vehicle based on maneuver mode recognition, which belongs to the technical field of vehicle trajectory prediction. The present invention solves the problem that the existing methods do not consider the influence of different maneuver types on the vehicle trajectory, resulting in large errors in the trajectory prediction results obtained when using the existing methods to predict the vehicle trajectory. Based on the prediction of the position and velocity at the current moment, the present invention distinguishes the observed state variables, identifies its maneuver mode, then determines the data fitting scheme, uses the fitted function, and predicts the position and velocity at the next moment by solving the differential equation until the trajectory prediction is completed. Compared with the traditional method, the method of the present invention improves the trajectory prediction accuracy and reduces the trajectory prediction error by identifying the maneuver mode. The present invention can be applied to the technical field of vehicle trajectory prediction.

Description

Technical Field

[0001] The present invention belongs to the technical field of aircraft trajectory prediction, and in particular relates to a near-space hypersonic aircraft trajectory prediction method based on maneuvering pattern recognition. Background Art

[0002] Near-space hypersonic vehicles refer to a type of space vehicle that cruises or glides at a speed of no less than 5 Ma in near-space. They have fast response and strong penetration capabilities, and are a type of weapon with strategic deterrence capabilities. After the hypersonic vehicle is boosted to a predetermined altitude, the hypersonic glider or hypersonic cruise body separates from the booster and continues to attack the predetermined target. The terminal attack speed may reach Mach 5 or above. It has strong long-range gliding and maneuvering flight capabilities in near-space.

[0003] The existing trajectory prediction method for near-space aircraft is to use the Singer model to estimate the aircraft state through Kalman filtering to obtain information on the aircraft position, velocity, and acceleration, and assume that the near-space hypersonic aircraft maintains a constant attitude angle for subsequent trajectory prediction. However, the existing method does not consider the impact of different maneuver types on the aircraft trajectory, resulting in large errors in the trajectory prediction results when the existing method is used to predict the aircraft trajectory. Summary of the invention

[0004] The purpose of the present invention is to solve the problem that the existing methods do not take into account the influence of different maneuver types on the trajectory of the aircraft, resulting in large errors in the trajectory prediction results obtained when the aircraft trajectory is predicted using the existing methods, and to propose a near-space hypersonic aircraft trajectory prediction method based on maneuver pattern recognition.

[0005] The technical solution adopted by the present invention to solve the above technical problems is:

[0006] A near-space hypersonic aircraft trajectory prediction method based on maneuver pattern recognition, wherein the method uses a Kalman filter to track the state information of the hypersonic aircraft during the tracking phase, and the state estimation value at the end of the tracking phase is converted into Denoted as:

[0007]

[0008] Among them, kT represents the end time of the tracking phase, x p (kT),y p (kT) and z p (kT) is the position of the hypersonic vehicle in the x-axis, y-axis and z-axis of the observation inertial coordinate system at the end of tracking kT, v xp (kT), v yp(kT) and v zp (kT) are the x-axis, y-axis and z-axis velocities of the hypersonic vehicle in the observation inertial coordinate system at the time kT when the tracking ends; and They are respectively the tracking end time kT for Z x (kT), Z y (kT) and Z z The estimated value of (kT), Z x (kT), Z y (kT) and Z z (kT) is the state variable at the end of tracking kT, and x p (kT),y p (kT), z p (kT), v xp (kT), v yp (kT) and v zp (kT) is used as the initial value for trajectory prediction;

[0009] According to the state variable Z in the tracking process x , Z y and Z z Determine the maneuver type of the hypersonic vehicle based on the changing characteristics of the vehicle;

[0010] If the maneuver type is a sinusoidal acceleration maneuver, use the model f(t) = a 1 sin(a 2 t+a 3 ) for Z x , Z y , Z z Fitting is performed to obtain the function Z x (t), Z y (t), Z z (t), a 1 、a 2 and a 3 is the fitting coefficient;

[0011] If the maneuver type is a constant acceleration maneuver, the model f(t) = a is used. 4 t+a 5 Z x , Z y , Z z Fitting function Z x (t), Z y (t), Z z (t), a 4 and a 5 is the fitting coefficient;

[0012] According to the fitted function Z x (t), Zy (t) and Z z (t) and the initial value of trajectory prediction are used to predict the trajectory of near-space hypersonic vehicles.

[0013] Furthermore, the observation inertial coordinate system is defined as:

[0014] The observation inertial coordinate system takes the observation base station as the origin O, and the line from the center of the earth to the observation base station pointing to the sky as the positive direction of the y-axis; in the plane formed by the observation base station, the center of the earth and the target, the direction from the origin to the target is the positive direction of the x-axis; and the z-axis direction is determined according to the right-hand rule.

[0015] Furthermore, the state variable Z in the tracking process x , Z y and Z z When judging the maneuver type of a hypersonic vehicle based on the changing characteristics of

[0016] Furthermore, the function Z obtained by fitting x (t), Z y (t) and Z z (t) and the initial value of trajectory prediction are used to predict the trajectory of near-space hypersonic vehicles. The specific process is as follows:

[0017] According to the fitted function Z x (t), Z y (t) and Z z (t), calculate the acceleration a of the hypersonic vehicle in the x-axis, y-axis and z-axis of the observation inertial coordinate system at time (k+1)T (k+1)T,x 、a (k+1)T,y and a (k+1)T,z :

[0018]

[0019] Where ρ is the air density, v t(k+1)T,x 、v t(k+1)T,y and v t(k+1)T,z are the x-axis, y-axis and z-axis velocities of the hypersonic vehicle in the observation inertial coordinate system at time (k+1)T respectively; is the transformation matrix from the ballistic coordinate system to the observation inertial coordinate system; Z x ((k+1)T) is based on the function Z x (t) The state variable at time (k+1)T, Z y ((k+1)T) is based on the function Z y (t) The state variable at time (k+1)T, Z z ((k+1)T) is based on the function Zz (t) The state variable obtained at time (k+1)T;

[0020] According to the initial value of trajectory prediction and a (k+1)T,x 、a (k+1)T,y and a (k+1)T,z , get the predicted values ​​of the position and velocity of the hypersonic vehicle at time (k+1)T in the observation inertial coordinate system;

[0021] Then, based on the predicted values ​​of the position and velocity of the hypersonic aircraft in the observed inertial coordinate system at time (k+1)T and the acceleration of the hypersonic aircraft in the x-axis, y-axis and z-axis directions of the hypersonic aircraft in the observed inertial coordinate system at time (k+2)T, the predicted values ​​of the position and velocity of the hypersonic aircraft in the observed inertial coordinate system at time (k+2)T are obtained, and so on, until the predicted values ​​of the position and velocity of the hypersonic aircraft at each time are obtained, thus completing the trajectory prediction of the near-space hypersonic aircraft.

[0022] Furthermore, the ballistic coordinate system is defined as:

[0023] The ballistic coordinate system takes the observation base station as the origin O′, the target velocity direction as the positive direction of the x′ axis, the y′ axis is in the plane formed by the observation base station, the center of the earth and the target, the y′ axis is perpendicular to the x′ axis, and the y′ axis is positive when it points to the target direction; the direction of the z′ axis is determined according to the right-hand rule.

[0024] Furthermore, the initial value of the trajectory prediction and a (k+1)T,x 、a (k+1)T,y and a (k+1)T,z , and obtain the predicted values ​​of the position and velocity of the hypersonic vehicle at time (k+1)T in the observation inertial coordinate system; the specific process is:

[0025] The differential equation is established as follows:

[0026]

[0027]

[0028]

[0029]

[0030]

[0031]

[0032] Where, μ is the gravitational coefficient; R e is the radius of the earth; v (k+1)T,x 、v (k+1)T,y and v(k+1)T,z are the x-axis, y-axis and z-axis velocities of the hypersonic vehicle in the observation inertial coordinate system at time (k+1)T; (k+1)T ,y (k+1)T and z (k+1)T are the position coordinates of the hypersonic vehicle in the observation inertial coordinate system at time (k+1)T along the x-axis, y-axis and z-axis respectively; For x (k+1)T The first derivative of v (k+1)T,y The first derivative of For z (k+1)T The first derivative of v (k+1)T,x The first derivative of v (k+1)T,y The first derivative of v (k+1)T,z The first derivative of ;

[0033] According to the initial value of trajectory prediction and a (k+1)T,x 、a (k+1)T,y and a (k+1)T,z Solve the differential equation to obtain the predicted values ​​of the position and velocity of the hypersonic vehicle at time (k+1)T in the observed inertial coordinate system.

[0034] Furthermore, the initial value of the trajectory prediction and a (k+1)T,x 、a (k+1)T,y and a (k+1)T,z To solve the differential equations, the Runge-Kutta method is used.

[0035] The beneficial effects of the present invention are:

[0036] The present invention is aimed at near-space hypersonic aircraft, taking into account the influence of complex conditions such as aerodynamic forces on the target, and in the case of unknown enemy aircraft parameters such as aircraft mass, reference area, thrust size, and aerodynamic parameters, based on the prediction of the position and speed at the current moment, the observed state quantities are distinguished, the maneuvering mode is identified, and then the data fitting scheme is determined. The function obtained by fitting is used to predict the position and speed at the next moment by solving differential equations until the trajectory prediction is completed. Compared with the traditional method, the method of the present invention improves the trajectory prediction accuracy and reduces the trajectory prediction error by identifying the maneuvering mode. BRIEF DESCRIPTION OF THE DRAWINGS

[0037] The simulation of the target maneuvering mode is sinusoidal acceleration maneuvering:

[0038] Figure 1 It is the longitudinal plane flight trajectory diagram of the hypersonic vehicle in the observation inertial coordinate system in the simulation experiment;

[0039] Figure 2 It is the trajectory diagram of the hypersonic vehicle flying in the lateral plane in the observation inertial coordinate system in the simulation experiment;

[0040] Figure 3 This is the acceleration estimation error diagram of the X-axis of the hypersonic vehicle in the simulation experiment;

[0041] In the figure: t represents time, Error in a x Represents the acceleration estimation error of the X-axis;

[0042] Figure 4 This is the acceleration estimation error diagram of the Y-axis of the hypersonic vehicle in the simulation experiment;

[0043] Error in a y Represents the acceleration estimation error of the Y axis;

[0044] Figure 5 This is the acceleration estimation error diagram of the Z-axis of the hypersonic vehicle in the simulation experiment;

[0045] Error in a z Represents the acceleration estimation error of the Z axis;

[0046] Figure 6 This is the velocity estimation error diagram of the X-axis of the hypersonic vehicle in the simulation experiment;

[0047] Error in v x Represents the velocity estimation error of the X-axis;

[0048] Figure 7 This is the velocity estimation error diagram of the Y-axis of the hypersonic vehicle in the simulation experiment;

[0049] Error in v y Represents the velocity estimation error of the Y axis;

[0050] Figure 8 This is the velocity estimation error diagram of the Z-axis of the hypersonic vehicle in the simulation experiment;

[0051] Error in v z Represents the velocity estimation error of the Z axis;

[0052] Fig. 9 This is the X-axis position estimation error diagram of the hypersonic vehicle in the observation coordinate system in the simulation experiment;

[0053] Error in x represents the position estimation error of the X-axis;

[0054] Fig.10This is the Y-axis position estimation error diagram of the hypersonic vehicle in the observation coordinate system in the simulation experiment;

[0055] Error in y represents the position estimation error of the Y axis;

[0056] Fig.11 This is the Z-axis position estimation error diagram of the hypersonic vehicle in the observation coordinate system in the simulation experiment;

[0057] Error in z represents the position estimation error of the Z axis;

[0058] Fig.12 This is a comparison chart of the estimated, predicted and actual values ​​of the X-axis acceleration of the hypersonic vehicle in the observation coordinate system in the simulation experiment;

[0059] Fig.13 This is a comparison chart of the estimated, predicted and actual values ​​of the Y-axis acceleration of the hypersonic vehicle in the observation coordinate system in the simulation experiment;

[0060] Fig.14 This is a comparison chart of the estimated, predicted and actual values ​​of the Z-axis acceleration of the hypersonic vehicle in the observation coordinate system in the simulation experiment;

[0061] Fig.15 This is a comparison chart of the estimated, predicted and actual values ​​of the X-axis velocity of the hypersonic vehicle in the observation coordinate system in the simulation experiment;

[0062] Fig.16 This is a comparison chart of the estimated, predicted and actual values ​​of the Y-axis velocity of the hypersonic vehicle in the observation coordinate system in the simulation experiment;

[0063] Fig.17 This is a comparison chart of the estimated, predicted and actual values ​​of the Z-axis velocity of the hypersonic vehicle in the observation coordinate system in the simulation experiment;

[0064] Fig.18 This is a comparison chart of the estimated, predicted and actual values ​​of the X-axis position of the hypersonic vehicle in the observation coordinate system in the simulation experiment;

[0065] Fig.19 This is a comparison chart of the estimated, predicted and actual values ​​of the Y-axis position of the hypersonic vehicle in the observation coordinate system in the simulation experiment;

[0066] Fig. 20 This is a comparison chart of the estimated, predicted and actual values ​​of the Z-axis position of the hypersonic vehicle in the observation coordinate system in the simulation experiment;

[0067] Fig.21 This is the error curve diagram of the position estimation, prediction and actual value of the hypersonic aircraft in the simulation experiment;

[0068] error in R represents the position error between the estimated, predicted and actual values;

[0069] Fig. 22 This is a statistical result diagram of the terminal position prediction error of the hypersonic vehicle in the simulation experiment;

[0070] P Rerror represents the terminal position prediction error;

[0071] Fig.23 is the state quantity Z of the hypersonic vehicle in the simulation experiment x Fitting curve result graph;

[0072] Fig.24 is the state quantity Z of the hypersonic vehicle in the simulation experiment y Fitting curve result graph;

[0073] Fig.25 is the state quantity Z of the hypersonic vehicle in the simulation experiment z Fitting curve result graph;

[0074] The target maneuvering mode is the simulation under the condition of constant acceleration maneuvering:

[0075] Fig.26 This is the longitudinal plane flight trajectory diagram of the hypersonic vehicle in the simulation experiment in the observation inertial coordinate system;

[0076] Fig. 27 It is the trajectory diagram of the hypersonic vehicle flying in the lateral plane in the observation inertial coordinate system in the simulation experiment;

[0077] Fig.28 This is the acceleration estimation error diagram of the X-axis of the hypersonic vehicle in the simulation experiment;

[0078] In the figure: t represents time, Error in a y Represents the acceleration estimation error of the X-axis;

[0079] Fig.29 This is the acceleration estimation error diagram of the Y-axis of the hypersonic vehicle in the simulation experiment;

[0080] Error in a y Represents the acceleration estimation error of the Y axis;

[0081] Fig.30 This is the acceleration estimation error diagram of the Z-axis of the hypersonic vehicle in the simulation experiment;

[0082] Error in a z Represents the acceleration estimation error of the Z axis;

[0083] Fig.31This is the velocity estimation error diagram of the X-axis of the hypersonic vehicle in the simulation experiment;

[0084] Error in v x Represents the velocity estimation error of the X-axis;

[0085] Fig.32 This is the velocity estimation error diagram of the Y-axis of the hypersonic vehicle in the simulation experiment;

[0086] Error in v y Represents the velocity estimation error of the Y axis;

[0087] Fig.33 This is the velocity estimation error diagram of the Z-axis of the hypersonic vehicle in the simulation experiment;

[0088] Error in v z Represents the velocity estimation error of the Z axis;

[0089] Fig.34 This is the X-axis position estimation error diagram of the hypersonic vehicle in the observation coordinate system in the simulation experiment;

[0090] Error in x represents the position estimation error of the X-axis;

[0091] Fig.35 This is the Y-axis position estimation error diagram of the hypersonic vehicle in the observation coordinate system in the simulation experiment;

[0092] Error in y represents the position estimation error of the Y axis;

[0093] Fig.36 This is the Z-axis position estimation error diagram of the hypersonic vehicle in the observation coordinate system in the simulation experiment;

[0094] Error in z represents the position estimation error of the Z axis;

[0095] Fig.37 This is a comparison chart of the estimated, predicted and actual values ​​of the X-axis acceleration of the hypersonic vehicle in the observation coordinate system in the simulation experiment;

[0096] Fig.38 This is a comparison chart of the estimated, predicted and actual values ​​of the Y-axis acceleration of the hypersonic vehicle in the observation coordinate system in the simulation experiment;

[0097] Fig.39 This is a comparison chart of the estimated, predicted and actual values ​​of the Z-axis acceleration of the hypersonic vehicle in the observation coordinate system in the simulation experiment;

[0098] Fig.40 This is a comparison chart of the estimated, predicted and actual values ​​of the X-axis velocity of the hypersonic vehicle in the observation coordinate system in the simulation experiment;

[0099] Fig.41 This is a comparison chart of the estimated, predicted and actual values ​​of the Y-axis velocity of the hypersonic vehicle in the observation coordinate system in the simulation experiment;

[0100] Fig.42 This is a comparison chart of the estimated, predicted and actual values ​​of the Z-axis velocity of the hypersonic vehicle in the observation coordinate system in the simulation experiment;

[0101] Fig.43 This is a comparison chart of the estimated, predicted and actual values ​​of the X-axis position of the hypersonic vehicle in the observation coordinate system in the simulation experiment;

[0102] Fig.44 This is a comparison chart of the estimated, predicted and actual values ​​of the Y-axis position of the hypersonic vehicle in the observation coordinate system in the simulation experiment;

[0103] Fig.45 This is a comparison chart of the estimated, predicted and actual values ​​of the Z-axis position of the hypersonic vehicle in the observation coordinate system in the simulation experiment;

[0104] Fig.46 This is the error curve diagram of the position estimation, prediction and actual value of the hypersonic aircraft in the simulation experiment;

[0105] error in R represents the position error between the estimated, predicted and actual values;

[0106] Fig.47 This is a statistical result diagram of the terminal position prediction error of the hypersonic vehicle in the simulation experiment;

[0107] P Rerror represents the terminal position prediction error;

[0108] Fig.48 is the state quantity Z of the hypersonic vehicle in the simulation experiment x Fitting curve result graph;

[0109] Fig.49 is the state quantity Z of the hypersonic vehicle in the simulation experiment y Fitting curve result graph;

[0110] Fig.50 is the state quantity Z of the hypersonic vehicle in the simulation experiment z Fitting curve result graph;

[0111] Fig.51 This is a schematic diagram of the possible landing point range and flight corridor of the aircraft. DETAILED DESCRIPTION

[0112] The method of the present invention is further described in detail below in conjunction with the accompanying drawings and prior art:

[0113] A near-space hypersonic vehicle trajectory prediction method based on maneuvering pattern recognition, the specific execution steps of the prediction method are as follows:

[0114] Step 1: Establish the state equation of the near-space hypersonic vehicle; the specific process of step 1 is:

[0115] Introduce three state variables Z x , Z y and Z z , and the state variable Z x , Z y and Z z The expressions are:

[0116]

[0117] Where: C x0 , C y0 and C z0 They represent the aerodynamic coefficients of the hypersonic vehicle along the x0 axis, y0 axis and z0 axis in the velocity coordinate system respectively; C p0 is the thrust coefficient; m T is the mass of the hypersonic vehicle; S is the equivalent reference area of ​​the hypersonic vehicle;

[0118] Select the state variable X as:

[0119] X=[xyzv x v y v z Z x Z y Z z ] T (2)

[0120] Where: x, y and z are the position coordinates of the hypersonic vehicle in the x-axis, y-axis and z-axis in the observation inertial coordinate system respectively; v x 、v y and v z are the x-axis, y-axis and z-axis velocities of the hypersonic vehicle in the observation inertial coordinate system, respectively;

[0121] Then the state equation of near-space hypersonic vehicle is established as:

[0122] The state equation of the near-space hypersonic aircraft includes an acceleration calculation equation, a velocity calculation equation, and a position calculation equation;

[0123]

[0124] in: is the first-order derivative of x, is the first-order derivative of y, is the first-order derivative of z, v x The first derivative of v y The first derivative of v z The first derivative of Z x The first derivative of Z y The first derivative of Z z The first derivative of x 、a y and a z are the accelerations of the hypersonic vehicle in the x-axis, y-axis and z-axis in the observation inertial coordinate system; w 1 、w 2 and w 3 They are all Gaussian white noise, representing a slow time-varying random walk process.

[0125] The intermediate variable a tx 、a ty and a tz The intermediate variable matrix and the intermediate variable g x , g y and g z The intermediate variable matrix The expression is as follows:

[0126]

[0127]

[0128] Where: ρ is the air density; represents the speed of the hypersonic vehicle; the velocity coordinate system is replaced by the ballistic coordinate system, then is the transformation matrix from the ballistic coordinate system to the observation inertial coordinate system ( is the transformation matrix from the velocity coordinate system to the observation inertial coordinate system. For axisymmetric aircraft, lateral maneuvers are generally controlled by the yaw channel through the sideslip angle. In actual flight, the roll angle of the aircraft is generally small, so the velocity coordinate system can be replaced by the ballistic coordinate system. Therefore, represents the transformation matrix from the ballistic coordinate system to the observation inertial coordinate system); μ is the gravitational coefficient; R e is the radius of the Earth;

[0129] Transformation matrix from ballistic coordinate system to observation inertial coordinate system The expression is:

[0130]

[0131] in:

[0132]

[0133] Step 2: Design a Kalman filter for near-space hypersonic aircraft state tracking, and use the designed Kalman filter to track the state information of the hypersonic aircraft in real time; the specific process of step 2 is:

[0134] The prediction equation of the Kalman filter for near-space hypersonic vehicle state tracking is:

[0135]

[0136] in, is the Kalman filter’s estimate of the state variable X, represents the state equation of the hypersonic aircraft system using the estimated values ​​of the state variables; G is the system noise matrix; Q is the hypersonic aircraft system state model error covariance matrix; R n is the measurement error covariance matrix; superscript -1 represents matrix inversion operation; superscript T represents the matrix transpose operation; H is the measurement matrix of the Kalman filter; P is the covariance matrix of the Kalman filter, is the first-order derivative of P with respect to time; F is the Jacobian matrix;

[0137]

[0138]

[0139] The Runge-Kutta method is used to solve formula (8) and the intermediate variable is obtained and the value of P(k+1 / k);

[0140] The measurement correction equation of the Kalman filter used for near-space hypersonic aircraft state tracking is:

[0141]

[0142] Where: I represents the identity matrix, K(k+1) represents the filter gain matrix of the k+1th step of the Kalman filter, η(k+1) represents the measurement information of the k+1th step of the Kalman filter, is the estimated value of the Kalman filter at the k+1th step, and P(k+1) is the value of the covariance matrix of the Kalman filter at the k+1th step;

[0143] η(k+1)=HX+R n(12)

[0144] By combining formula (8) and formula (11), a Kalman filter for near-space hypersonic aircraft state tracking is obtained. The designed Kalman filter is used to track the state information of the hypersonic aircraft in real time.

[0145] Step 3: Obtain the state information of the hypersonic aircraft from the start time 0 to the end time kT of tracking according to the Kalman filter designed in step 2, and obtain the state estimation value of the hypersonic aircraft in the tracking stage according to the state information obtained by tracking;

[0146] Wherein, k is an integer greater than 0, T is the calculation period, which can generally be taken as T = 0.01s; the state estimation value of the hypersonic aircraft at the time kT when the tracking ends is used as the initial value of the trajectory prediction;

[0147] The specific process of step three is:

[0148] The Kalman filter designed in step 2 obtains the state information of the hypersonic aircraft from the start time 0 to the end time kT of tracking, and obtains the state estimation value of the hypersonic aircraft at the end time kT of tracking based on the state information obtained by tracking as follows:

[0149]

[0150] Where: x p (kT),y p (kT) and z p (kT) is the position of the hypersonic vehicle on the x-axis, y-axis and z-axis in the observation inertial coordinate system at the time kT when the tracking ends, v xp (kT), v yp (kT) and v zp (kT) are the x-axis, y-axis and z-axis velocities of the hypersonic vehicle in the observation inertial coordinate system at the time kT after the tracking ends; and They are respectively the tracking end time kT for Z x , Z y and Z z Estimated value;

[0151] x p (kT),y p (kT), z p (kT), v xp (kT), v yp (kT) and v zp (kT) is used as the initial value for trajectory prediction.

[0152] From equations (3) to (5), it can be seen that for near-space hypersonic aircraft, the acceleration at the end of the flight is generated by aerodynamic force and the earth's gravity. Before the hypersonic aircraft enters the final downward pressure, it needs to maneuver. Since the attitude angle of the aircraft usually changes; considering the large maneuverability of the hypersonic aircraft, the aerodynamic coefficient C x , C y , C z The changes are relatively large.

[0153] From formula (1), we can see that under maneuvering conditions, Z x , Z y , Z z The change is large, and the maneuvering mode is unknown. Therefore, according to the estimated value of the tracking phase, the maneuvering type of the hypersonic aircraft is identified by using the cluster analysis method, and then the prediction scheme is determined, the state quantity is fitted, and the fitting result is used as the basis for the prediction. x , Z y , Z z The value of is determined by the fitting results.

[0154] In formula (4), ρ, v, The influence of the earth's gravity on the acceleration change is relatively large. The air density ρ is only related to the height of the hypersonic aircraft, and the height of the hypersonic aircraft and the magnitude of the earth's gravity are only related to the state quantities x, y, and z. The velocity v and the coordinate transformation matrix Only with state v x 、v y 、v z related.

[0155] That is, the predicted values ​​of the position and velocity of the x-axis, y-axis and z-axis in the observation coordinate system at time kT can be used to obtain the predicted value of the acceleration at time (k+1)T, and the predicted values ​​of the velocity and position of the x-axis, y-axis and z-axis in the observation coordinate system at time (k+1)T can be obtained according to formula (3), thereby realizing iterative calculation and making predictions throughout the process. This process is a differential equation with 6 parameters, which can be solved using the Runge-Kutta method.

[0156] Step 4: Based on the estimated value in the tracking phase, the maneuver type is identified by the cluster analysis method, and then the prediction scheme is determined, the state quantity is fitted, and the fitting result is used as the basis for prediction. The initial value of the trajectory prediction and the result of the fitting of the state quantity are used to solve the predicted value of the position and speed of the hypersonic aircraft at the time (k+1)T, and the predicted value of the position and speed of the hypersonic aircraft at the time (k+2)T is continued to be solved by using the predicted value of the position and speed at the time (k+1)T and the result of the fitting of the state quantity. The specific process of step 4 is:

[0157] In the process of modeling the target, the state variable Z x , Z y , Z z It is directly related to the target acceleration and uses the Z that enters the steady state during the estimation process. x , Z y , Z z The estimated value of Z x , Z y , Z z The change characteristics of X-51 can be divided into two types of maneuvers, namely sinusoidal acceleration maneuvers and constant acceleration maneuvers. First, the maneuver type of X-51 is determined by the clustering method. If the maneuver type is sinusoidal acceleration maneuvers, the model f(t) = a is used. 1 sin(a 2 t+a 3 ) for Z x , Z y , Z z Fitting function Z x (t), Z y (t), Z z (t); if the maneuver type is a constant acceleration maneuver, the model f(t) = a is used 1 t+a 2 Z x , Z y , Z z Fitting function Z x (t), Z y (t), Z z (t). x (t), Z y (t), Z z Substituting (t) into formula (4) we can get

[0158] Combining formula (3) to formula (7), we can obtain the following differential equation:

[0159]

[0160] According to x p (kT),y p (kT), z p (kT), v xp (kT), v yp (kT), v zp (kT) and The Runge-Kutta method is used to solve differential equation (14) to obtain the predicted values ​​of the position and velocity of the hypersonic vehicle in the observation coordinate system at time (k+1)T; where: x p ((k+1)T),yp ((k+1)T) and z p ((k+1)T) represents the predicted value of the position of the hypersonic vehicle in the x-axis, y-axis and z-axis in the observation coordinate system at time (k+1)T, v xp ((k+1)T),v yp ((k+1)T) and v zp ((k+1)T) represents the predicted value of the speed of the hypersonic vehicle in the x-axis, y-axis and z-axis directions in the observation coordinate system at time (k+1)T;

[0161] Similarly, according to x p ((k+1)T),y p ((k+1)T),z p ((k+1)T),v xp ((k+1)T),v yp ((k+1)T),v zp ((k+1)T) and The Runge-Kutta method is used to solve the differential equation (14) to obtain the predicted values ​​of the position and velocity of the hypersonic vehicle in the x-axis, y-axis and z-axis directions in the observation coordinate system at time (k+2)T.

[0162] The specific process of the clustering method is:

[0163] Step 1: First, we need to obtain a certain sample

[0164] The maneuvering mode of the target aircraft is set to the sine form, and the acceleration is set to sine signals of different frequencies. Only the tracking process is performed, and the data estimated in the tracking process is stored and recorded as X = {x 1 T ,x 2 T ,......,x n T} T , where n is the number of times the tracking process is repeated; the maneuvering mode of the target aircraft is set to a constant value form, and the acceleration is set to a constant value signal of different frequencies, and only the tracking process is performed. The data estimated by the tracking process is stored and recorded as Y = {y 1 T ,y 2 T ,......,y m T} T , where m is the number of times the tracking process is repeated. The classified m+n data are used as row vectors to obtain a matrix A={X T ,Y T} T , note A={A 1T ,A 2 T ,......,A m+n T} T .

[0165] Step 2: Next, input the sample to be tested and classify it together with the data sample obtained in step 1. The classification adopts the linear discriminant analysis (LDA) method, which can be regarded as a linear classifier. LDA hopes that the distance between different categories of the projected data is as large as possible, and the same category is as compact as possible. Assume that a∈A is represented as a one-dimensional vector b through linear projection transformation, and the transformation relationship is:

[0166] b=wa T (15)

[0167] Then we further set a threshold b 0 , if b≥b 0 When , it is judged to belong to the first category, otherwise it belongs to the second category. Therefore, the problem of LDA can be understood as how to adjust the projection direction w to ensure that b can distinguish the two categories as much as possible in one-dimensional space.

[0168] LDA considers both inter-class and intra-class relationships. On the one hand, we hope that the distance between the projected means of the data of two categories is as large as possible. Assume that the mean of category one is The mean of category 2 is Then the distance between the two types of means after projection is

[0169] On the other hand, we hope that the data of the same category are as concentrated as possible. Assume that the data in category i is represented as i=1,2, then the intra-class covariance of the original data is expressed as

[0170]

[0171] The intra-class variance of the projected data is expressed as:

[0172]

[0173] Combining the maximum class distance and the minimum intra-class variance, the optimization problem solved by LDA is:

[0174]

[0175] Solved

[0176] Step 5. Repeat the process of step 4 until the predicted values ​​of the position and speed of the hypersonic aircraft at each moment are obtained, thus completing the trajectory prediction of the near-space hypersonic aircraft.

[0177] Simulation part

[0178] This simulation lasted for 400 s, and the state of the near-space vehicle was estimated from 82 to 282 s;

[0179] like Figures 1 to 51 As shown, for the attacker, the near-space aircraft designed by the present invention, whether it is a cruise hypersonic aircraft or a gliding hypersonic aircraft, will have requirements for terminal velocity and attack angle when attacking the target point, and the terminal velocity and attack angle are usually required to reach a certain minimum value. The simulation results show that if the terminal velocity is required to be greater than 2Mach and the attack angle is greater than 70°, the speed needs to be around 2800m / s when entering the final downward pressure point. Under the constraints of this actual situation, the maneuverability of near-space hypersonic aircraft in the terminal stage has been constrained. If the speed and direction of the aircraft are known at 1000km from the target point, the possible landing point range of the aircraft is approximately a circular area with a radius of 50km with the target point as the center. In this case, if the final landing point of the aircraft is determined, the width of the flight corridor in the final stage is only about 20km (about 10km in the positive and negative directions of the z-axis). The relationship between the possible landing point range of the aircraft and the width of the flight corridor is as follows Fig.51 As shown. Under the condition that the target maneuvering mode is unknown, the target maneuvering mode can be identified based on the state quantity obtained by filter tracking, and this can be used as the basis for selecting the state fitting model. Select an appropriate fitting model to fit the state quantity, and the result is as follows Fig.23 , Fig.24 , Fig.25 and Fig.48 , Fig.49 , Fig.50 shown.

[0180] In the simulation of the present invention, the final attack point of the aircraft is directly set as the origin of the observation coordinate system, so the maximum lateral maneuvering range of the aircraft is about 20km. If the attack point of the aircraft takes any value within the possible landing point range and does not coincide with the origin of the observation coordinate system, it will not affect the results of state estimation and prediction. In addition, for the trajectory prediction of near-space aircraft, it is necessary to know the approximate target point of the aircraft, and detection devices such as radar are usually deployed near the defense target. At the same time, the purpose of trajectory prediction is to prepare for the subsequent interception of near-space aircraft, and only the prediction of the trajectory before the downward pressure point has the meaning of interception. If the approximate attack range of the enemy aircraft is uncertain, the aircraft may enter the downward pressure section to attack at any time. In this case, it is difficult to predict the interception. By Fig.51It can be seen that after the possible attack range of the aircraft is limited to a circle with a radius of 50km, the lateral value range of the final pressure point is about 120km. The significance of the research of this invention lies in that after 200s of state estimation, the corridor range of 120km width before the pressure point at 400s can be narrowed to a range less than 10km through trajectory prediction.

[0181] Under the above conditions, 100 Monte-Carlo simulations were performed when the target maneuvering mode was sinusoidal acceleration maneuver and constant acceleration maneuver. The statistical results of the final position prediction error of the near-space hypersonic aircraft obtained after 100 Monte-Carlo simulations are as follows: Fig. 22 and Fig.47 As shown, when the target maneuvering mode is a sinusoidal acceleration maneuver, the terminal position prediction error calculated by statistics is less than 10km, and when the target maneuvering mode is a constant acceleration maneuver, the terminal position prediction error calculated by statistics is less than 3km.

[0182] The above calculation examples of the present invention are only used to explain the calculation model and calculation process of the present invention in detail, and are not intended to limit the implementation methods of the present invention. For ordinary technicians in the relevant field, other different forms of changes or modifications can be made based on the above description. It is impossible to list all the implementation methods here. All obvious changes or modifications derived from the technical solution of the present invention are still within the scope of protection of the present invention.

Claims

1. A method for predicting the trajectory of a near-space hypersonic vehicle based on maneuver mode recognition, characterized in that, in the tracking stage of the method, a Kalman filter is used to track the state information of the hypersonic vehicle, and the prediction equation of the Kalman filter is: Among them, is the estimated value of the state variable X by the Kalman filter, represents the hypersonic vehicle system state equation using the estimated value of the state variable; G is the system noise matrix; Q is the hypersonic vehicle system state model error covariance matrix; R n is the measurement error covariance matrix; the superscript -1 represents matrix inversion operation; the superscript T represents matrix transpose operation; H is the measurement matrix of the Kalman filter; P is the covariance matrix of the Kalman filter, is the first derivative of P with respect to time; F is the Jacobian matrix; The Runge-Kutta method is used to solve formula (8) to obtain the intermediate variables and the value of P(k+1 / k); The measurement correction equation of the Kalman filter is: Where: I represents the identity matrix, K(k + 1) represents the filtering gain matrix of the (k + 1)-th step of the Kalman filter, η(k + 1) represents the measurement information of the (k + 1)-th step of the Kalman filter, is the estimated value of the (k + 1)-th step of the Kalman filter, and P(k + 1) is the value of the covariance matrix of the (k + 1)-th step of the Kalman filter; η(k + 1) = HX + R n (12) Combining formula (8) and formula (11), a Kalman filter for tracking the state of a near-space hypersonic vehicle is obtained; The state estimate value at the end of the tracking phase is denoted as: where, kT represents the end moment of the tracking phase, x p (kT), y p (kT), and z p (kT) are the positions of the hypersonic vehicle along the x-axis, y-axis, and z-axis in the observed inertial coordinate system at the end moment kT of tracking, and v xp (kT), v yp (kT), and v zp (kT) are the velocities of the hypersonic vehicle along the x-axis, y-axis, and z-axis in the observed inertial coordinate system at the end moment kT of tracking, respectively; and are the estimated values of Z x (kT), Z y (kT), and Z z (kT) at the end moment kT of tracking, Z x (kT), Z y (kT), and Z z (kT) are the state variables at the end moment kT of tracking. Take x p (kT), y p (kT), z p (kT), v xp (kT), v yp (kT), and v zp (kT) as the initial values for ballistic prediction; Based on the state variables Z x , Z y and Z z during the tracking process, the maneuver type of the hypersonic vehicle is judged; If the maneuver type is a sinusoidal acceleration maneuver, then use the model f(t) = a 1 sin(a 2 t + a 3 ) to fit Z x , Z y , and Z z respectively, to obtain the functions Z x (t), Z y (t), Z z (t), where a 1 , a 2 , and a 3 are fitting coefficients; If the maneuver type is a constant acceleration maneuver, then the model f(t) = a 4 t + a 5 is used to fit Z x , Z y , and Z z respectively to obtain the functions Z x (t), Z y (t), and Z z (t), where a 4 and a 5 are fitting coefficients; According to the functions Z x (t), Z y (t) and Z z (t) obtained by fitting and the initial values of the ballistic prediction, the trajectory prediction of the hypersonic vehicle in the near space is carried out.

2. The method for predicting the trajectory of a near-space hypersonic vehicle based on maneuver mode recognition according to claim 1, characterized in that, the definition of the observation inertial coordinate system is: The observation inertial coordinate system takes the observation base station as the origin O, and the positive direction of the y-axis is the direction from the center of the earth to the observation base station pointing to the sky; in the plane formed by the observation base station, the center of the earth and the target, the direction from the origin to the target is the positive direction of the x-axis; then the direction of the z-axis is determined according to the right-hand rule.

3. The method for predicting the trajectory of a near-space hypersonic vehicle based on maneuver mode recognition according to claim 2, characterized in that, When determining the maneuver type of a hypersonic vehicle based on the variation characteristics of the state variables Z x , Z y and Z z , the clustering method is adopted.

4. The method for predicting the trajectory of a near-space hypersonic vehicle based on maneuver mode recognition according to claim 3, characterized in that, The function Z obtained by fitting x (t), Z y (t) and Z z (t) and the initial values of the ballistic prediction are used to perform the ballistic prediction of the hypersonic vehicle in the near space; the specific process is as follows: According to the functions Z x (t), Z y (t) and Z z (t), calculate the accelerations ax (k+1)T,x , ay (k+1)T,y and az (k+1)T,z of the hypersonic vehicle in the x-axis, y-axis and z-axis directions of the observation inertial coordinate system at the moment of (k + 1)T: where ρ is the air density, v t(k+1)T,x 、v t(k+1)T,y and v t(k+1)T,z are the velocities of the hypersonic vehicle at the (k + 1)T moment in the x-axis, y-axis, and z-axis directions of the observation inertial coordinate system, respectively; is the transformation matrix from the ballistic coordinate system to the observation inertial coordinate system; Z x ((k + 1)T) is the state variable at the (k + 1)T moment obtained according to the function Z x (t), Z y ((k + 1)T) is the state variable at the (k + 1)T moment obtained according to the function Z y (t), Z z ((k + 1)T) is the state variable at the (k + 1)T moment obtained according to the function Z z (t); According to the initial values of the ballistic prediction and a (k+1)T,x 、a (k+1)T,y and a (k+1)T,z , the predicted values of the position and velocity of the hypersonic vehicle at the (k + 1)T moment in the observation inertial coordinate system are obtained; Furthermore, based on the predicted values of the position and velocity of the hypersonic vehicle at time (k + 1)T in the observation inertial coordinate system and the accelerations of the hypersonic vehicle in the x-axis, y-axis, and z-axis directions of the observation inertial coordinate system at time (k + 2)T, the predicted values of the position and velocity of the hypersonic vehicle at time (k + 2)T in the observation inertial coordinate system are obtained, and so on, until the predicted values of the position and velocity of the hypersonic vehicle at each moment are obtained, and the trajectory prediction of the near-space hypersonic vehicle is completed.

5. The method for predicting the trajectory of a near-space hypersonic vehicle based on maneuver mode recognition according to claim 4, characterized in that, the definition of the trajectory coordinate system is: The trajectory coordinate system takes the observation base station as the origin O′, the direction of the target velocity is the positive direction of the x′-axis, the y′-axis is in the plane formed by the observation base station, the center of the earth and the target, the y′-axis is perpendicular to the x′-axis, and the positive direction of the y′-axis is the direction pointing to the target; then the direction of the z′-axis is determined according to the right-hand rule.

6. The method for predicting the trajectory of a near-space hypersonic vehicle based on maneuver mode recognition according to claim 5, characterized in that, The initial value predicted according to the trajectory and a (k+1)T,x 、a (k+1)T,y and a (k+1)T,z to obtain the predicted values of the position and velocity of the hypersonic vehicle at the (k + 1)T moment in the observed inertial coordinate system; the specific process is as follows: The following differential equation is established: where μ is the gravitational coefficient; R e is the radius of the Earth; v (k+1)T,x , v (k+1)T,y and v (k+1)T,z are the velocities of the hypersonic vehicle in the x-axis, y-axis, and z-axis directions in the observation inertial coordinate system at the (k + 1)T moment, respectively; x (k+1)T , y (k+1)T and z (k+1)T are the position coordinates of the hypersonic vehicle in the x-axis, y-axis, and z-axis directions in the observation inertial coordinate system at the (k + 1)T moment, respectively; is the first derivative of x (k+1)T , is the first derivative of v (k+1)T,y , is the first derivative of z (k+1)T , is the first derivative of v (k+1)T,x , is the first derivative of v (k+1)T,y , is the first derivative of v (k+1)T,z ; According to the initial values of the ballistic prediction and a (k+1)T,x 、a (k+1)T,y and a (k+1)T,z Solve the differential equation to obtain the predicted values of the position and velocity of the hypersonic vehicle at the (k + 1)T moment in the observed inertial coordinate system.

7. The method for predicting the trajectory of a near-space hypersonic vehicle based on maneuver mode recognition according to claim 6, characterized in that, The initial value predicted according to the trajectory and a (k+1)T,x 、a (k+1)T,y and a (k+1)T,z Solve the differential equation using the Runge-Kutta method.

Images