CN116611160A - Online real-time characteristic parameter identification and trajectory prediction method for uncontrolled aircraft based on measured trajectory parameters - Google Patents
Online real-time characteristic parameter identification and trajectory prediction method for uncontrolled aircraft based on measured trajectory parameters Download PDFInfo
- Publication number
- CN116611160A CN116611160A CN202310442699.6A CN202310442699A CN116611160A CN 116611160 A CN116611160 A CN 116611160A CN 202310442699 A CN202310442699 A CN 202310442699A CN 116611160 A CN116611160 A CN 116611160A
- Authority
- CN
- China
- Prior art keywords
- trajectory
- ballistic
- aircraft
- coefficient
- uncontrolled
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 51
- 238000001914 filtration Methods 0.000 claims abstract description 38
- 238000004364 calculation method Methods 0.000 claims abstract description 36
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 6
- 238000005259 measurement Methods 0.000 claims description 40
- 239000011159 matrix material Substances 0.000 claims description 10
- 238000012545 processing Methods 0.000 claims description 10
- 230000008569 process Effects 0.000 claims description 9
- 230000009466 transformation Effects 0.000 claims description 6
- 238000004590 computer program Methods 0.000 claims description 4
- 238000013016 damping Methods 0.000 claims description 4
- 230000001902 propagating effect Effects 0.000 claims description 3
- 230000001133 acceleration Effects 0.000 claims description 2
- 230000006870 function Effects 0.000 claims description 2
- 230000014509 gene expression Effects 0.000 claims description 2
- 230000005484 gravity Effects 0.000 claims description 2
- 230000010354 integration Effects 0.000 claims description 2
- 238000005096 rolling process Methods 0.000 claims description 2
- 230000003068 static effect Effects 0.000 claims description 2
- 230000001131 transforming effect Effects 0.000 claims description 2
- 101100272279 Beauveria bassiana Beas gene Proteins 0.000 claims 1
- 230000008859 change Effects 0.000 description 6
- 238000005516 engineering process Methods 0.000 description 5
- 230000000694 effects Effects 0.000 description 4
- 230000007547 defect Effects 0.000 description 3
- 238000009499 grossing Methods 0.000 description 3
- RZVHIXYEVGDQDX-UHFFFAOYSA-N 9,10-anthraquinone Chemical compound C1=CC=C2C(=O)C3=CC=CC=C3C(=O)C2=C1 RZVHIXYEVGDQDX-UHFFFAOYSA-N 0.000 description 2
- 238000007476 Maximum Likelihood Methods 0.000 description 2
- 238000012937 correction Methods 0.000 description 2
- 230000004069 differentiation Effects 0.000 description 2
- 239000006185 dispersion Substances 0.000 description 2
- 230000010355 oscillation Effects 0.000 description 2
- 239000002245 particle Substances 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 238000013459 approach Methods 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 238000010304 firing Methods 0.000 description 1
- 230000000644 propagated effect Effects 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 230000036962 time dependent Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/10—Geometric CAD
- G06F30/15—Vehicle, aircraft or watercraft design
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
- G06F17/12—Simultaneous equations, e.g. systems of linear equations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2113/00—Details relating to the application field
- G06F2113/08—Fluids
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Computational Mathematics (AREA)
- Data Mining & Analysis (AREA)
- General Engineering & Computer Science (AREA)
- Algebra (AREA)
- Geometry (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- Computer Hardware Design (AREA)
- Computing Systems (AREA)
- Evolutionary Computation (AREA)
- Operations Research (AREA)
- Fluid Mechanics (AREA)
- Automation & Control Theory (AREA)
- Aviation & Aerospace Engineering (AREA)
- Life Sciences & Earth Sciences (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Evolutionary Biology (AREA)
- Probability & Statistics with Applications (AREA)
- Feedback Control In General (AREA)
Abstract
The invention discloses a method for identifying aerodynamic parameters and forecasting trajectory of an uncontrolled aircraft based on a section of measured trajectory parameters, which comprises the steps of firstly establishing a simplified motion state equation suitable for online rapid calculation of the aircraft according to the properties of measured trajectory data and the characteristics of the subsequent flight trajectory of the aircraft, and establishing a trajectory filtering aerodynamic parameter identification model by applying an unscented Kalman filtering algorithm; then, inputting a section of measured ballistic data into the model to finish ballistic filtering, and identifying the aircraft resistance coincidence coefficient and the lift coincidence coefficient; and finally, taking the filter state value of the last point and the identified coincidence coefficient as initial values of subsequent flight trajectory calculation to carry out trajectory calculation, so as to realize the forecast of the subsequent flight trajectory. The method can carry out ballistic filtering and pneumatic parameter identification on a section of actually measured ballistic parameters on line in real time, rapidly and accurately forecast the subsequent flight trajectory of the aircraft, and has the characteristics of high ballistic calculation precision, small calculation amount and short time.
Description
Technical Field
The invention belongs to the technology of calculation of external trajectory of an uncontrolled aircraft, in particular to an online real-time characteristic parameter identification and trajectory prediction method of the uncontrolled aircraft based on a section of measured trajectory parameters.
Background
For various uncontrolled aircrafts (such as shells, rocket shells and the like), after a ballistic model is established, if a set of parameters (such as structural parameters, aerodynamic coefficients, launching condition parameters and the like) are given, a determined trajectory can be calculated, which is called as a theoretical trajectory, and the process is ballistic calculation. In practical engineering applications, theoretical trajectory is obtained before the aircraft shoots. Taking the projectile as an example, after the projectile is launched, parameters such as actual initial velocity, aerodynamic coefficient, initial disturbance and the like of the projectile not only are different from a given theoretical value, but also are different from each projectile, so that the actual trajectory of each projectile has a random difference, and in this sense, the actual trajectory of each projectile can be called as a random trajectory. The actual trajectory of each bullet can be measured by a certain sensor or an external measuring device (such as a missile-borne satellite signal receiving device or a trajectory tracking radar, etc.), namely, the related information of the random trajectory is measured. For a segment of actually measured ballistic data, one or two ballistic characteristic parameters can be selected to represent the random information of the trajectory, the ballistic characteristic parameters are identified from the actually measured ballistic data by utilizing a ballistic filtering technology, the random ballistic characteristic of the projectile is obtained in an equivalent way, and the subsequent flight trajectory state of the projectile can be predicted by combining the ballistic characteristic parameters with ballistic calculation, and the process is called 'ballistic prediction'. Clearly, the ballistic prediction is distinguished from the ballistic calculation, and the core of the ballistic prediction is the identification (or extraction) of the ballistic characteristic parameters, so that the randomness of the actual ballistic trajectory of the rocket is represented as the difference of the ballistic characteristic parameters of each projectile.
Therefore, the trajectory calculation is a process of performing data calculation on the flight trajectory of the rocket by the external trajectory model under known or set condition state parameters, wherein the condition state parameters include initial velocity, firing angle, rocket aerodynamic coefficient, rocket structure parameter, meteorological parameter and the like. The ballistic prediction is based on ballistic calculation, and according to a section of measured rocket motion information, the ballistic characteristic parameters are identified by utilizing a ballistic filtering technology, and the subsequent ballistic prediction process is performed, so that the real-time performance is emphasized or online processing is required. For example, for a long-range grenade, after shooting, a piece of ballistic data on an arc-up section can be obtained by using a ballistic tracking radar, and for the requirement of the flight control of the grenade, we are required to accurately forecast the actual landing point of the grenade by using the piece of ballistic data before the grenade lands (even before entering the arc-down section). Under the condition, the resistance coincidence coefficient and the lift coincidence coefficient in the filtering model can be selected as ballistic characteristic parameters, namely, all random disturbance equivalents on the trajectory are the changes of the resistance coefficient. When the Kalman filtering is constructed, the coincidence coefficients are used as a state variable, the optimal estimation of the coincidence coefficients can be obtained by utilizing measured ballistic data, so that the original resistance coefficient and lift coefficient (theoretical value) are corrected, and then ballistic calculation is carried out, so that the actual drop point of the bullet can be obtained rapidly. Since the actual landing point is known before the shot landing point, time is available for realizing flight control, and the flight trajectory is changed to approach the target point. The ballistic prediction is widely applied to modern bullets and arrows, particularly to controlled bullets and arrows, intelligent bullets and arrows and the like, and can effectively improve the precision, reduce the scattering and increase the range. From the above, the basis of ballistic prediction is pneumatic parameter identification.
Currently, the pneumatic parameter identification method mainly adopts a parameter differentiation method (also called a C-K method) and a maximum likelihood method. The method divides the trajectory into a plurality of small sections, the Mach number and the pneumatic parameters are regarded as constants on each small section, generally, tens to twenty measurement data are arranged on each small section, pneumatic parameter identification is carried out by establishing a conjugate equation set and solving the conjugate equation set, the pneumatic parameters corresponding to the Mach number on the small section can be obtained, and then a smooth technology is adopted to obtain a smooth curve of the pneumatic parameters changing along with the Mach number. The more variables of the original equation set of the method are, the more undetermined parameters are, the number of equations in the conjugate equation set is expanded sharply; and the measuring trajectory is required to be divided into a plurality of small sections, and after the pneumatic parameters on each small section are identified, the smoothness of the data is realized through a smoothing technology. The method has the advantages of large calculated data volume, long time consumption, suitability for data processing after shooting test and difficulty in meeting the requirement of on-line real-time rapid pneumatic parameter identification of a closed loop projectile calibrating system.
At present, for the ballistic calculation of conventional shells, because of dispersion, and the factor causing the dispersion is random, the ballistic calculation of each shell is difficult to accurately calculate in real time; however, since conventional shells belong to surface-killed weapons, the average ballistic effect is usually concerned, and therefore, the theoretical reference trajectory is mostly calculated offline under known shooting conditions according to determined aerodynamic ballistic parameters in the calculation of the trajectory. The method for calculating the theoretical reference trajectory offline is difficult to accurately forecast the subsequent actual flight trajectory of the uncontrolled aircraft according to the actual measured trajectory parameters of the uncontrolled aircraft in air flight.
Disclosure of Invention
The invention aims to provide an on-line real-time characteristic parameter identification and trajectory prediction method for an uncontrolled aircraft based on a section of measured trajectory parameters.
The technical solution for realizing the purpose of the invention is as follows: in a first aspect, the invention provides a method for identifying and forecasting the trajectory of an uncontrolled aircraft on line in real time based on a section of measured trajectory parameters, comprising the following steps:
step 1, selecting aerodynamic parameters of an uncontrolled aircraft as characteristic parameters, and establishing a ballistic filtering aerodynamic parameter identification model; firstly, according to the nature of measured trajectory data and the characteristics of the subsequent flight trajectory of an uncontrolled aircraft, a simplified five-degree-of-freedom rigid motion state equation is established, then a unscented Kalman filtering algorithm is applied to conduct linearization and discretization processing on the established simplified rigid motion state equation, and finally a Kalman ballistic filtering pneumatic parameter identification model is established;
step 2, filtering the measured ballistic parameters by using the ballistic filtering aerodynamic parameter identification model established in the step 1, and identifying the resistance coincidence coefficient and the lift coincidence coefficient of the uncontrolled aircraft to obtain initial conditions, resistance and lift coincidence coefficient for calculating the subsequent flight trajectory;
and 3, forecasting the subsequent flight trajectory by using the initial condition and the coincidence coefficient of the trajectory calculation obtained in the step 2, and obtaining the subsequent actual flight trajectory and trajectory drop point information of the uncontrolled aircraft.
In a second aspect, the present invention provides a computer device comprising a memory, a processor and a computer program stored on the memory and executable on the processor, the processor implementing the steps of the method of the first aspect when the program is executed.
In a third aspect, the present invention provides a computer readable storage medium having stored thereon a computer program which when executed by a processor performs the steps of the method of the first aspect.
Compared with the prior art, the invention has the remarkable advantages that: (1) According to the method, the unscented Kalman filtering algorithm is applied to perform ballistic filtering on the measured ballistic data, so that the noise of the measured data is reduced; according to the method, the theoretical dynamic model of the physical system is fully considered, the theoretical model is used for prediction, and then the measurement result is used for correction, so that the state value precision is improved; the method overcomes the defects that the traditional data smoothing method, the least square method and the like do not consider the characteristics of a physical system and data testing errors, but only carry out pure digital processing. (2) The method disclosed by the invention applies the established unscented Kalman ballistic filtering pneumatic parameter identification model, can identify and process the actual resistance coincidence coefficient and lift coincidence coefficient of the air flight uncontrolled aircraft which are influenced by various factors comprehensively according to a section of actually measured flight trajectory parameters, and further can accurately forecast the subsequent actual flight trajectory of the uncontrolled aircraft; the algorithm is simple in processing and high in calculation speed, overcomes the defects that the number of the equations of the conjugate equation set established by the parameter differentiation method and the maximum likelihood method is large, the sectional processing is needed, the identification parameter smoothing processing is needed, the processing process is complex, the number of the solved equations is large, and the like, and is not suitable for the requirement of online real-time rapid processing. (3) The uncontrolled aircraft established by the method simplifies a five-degree-of-freedom motion state equation, fully considers the influence of gyroscopic effect on lateral motion of the uncontrolled aircraft, performs trajectory calculation by using the identified actual resistance coincidence coefficient and lift coincidence coefficient, and has high accuracy of calculated trajectory falling points; in addition, the ballistic calculation can take a larger step length (such as 50 milliseconds), the ballistic calculation is very short (can be completed within tens of milliseconds in the current computer hardware environment), and the method is very suitable for online real-time rapid and accurate ballistic calculation; the method overcomes the defects that a six-degree-of-freedom trajectory model and a particle trajectory model in the existing outside ballistics are suitable for offline calculation of theoretical reference trajectory under known shooting conditions according to determined aerodynamic trajectory parameters, and accurate trajectory calculation is difficult to carry out on an actual uncontrolled aircraft flying in any air; the calculation step length required by the six-degree-of-freedom trajectory model is small (about several milliseconds), the calculation amount is large, the required calculation time is too long, and the method is not suitable for the requirement of online real-time rapid trajectory calculation; the influence of gyroscopic effect on lateral motion of the uncontrolled aircraft is not considered by the simple particle trajectory model, and the solution accuracy is poor although the trajectory solution is fast.
The invention is described in further detail below with reference to the accompanying drawings.
Drawings
FIG. 1 is a flow chart of the pneumatic parameter identification and ballistic prediction method based on a segment of measured ballistic parameters according to the present invention.
Fig. 2 is a graph of the change over time of the input horizontal distance actual measurement value and its filtered value.
FIG. 3 is a graph of the variation of the ballistic high measured value of the input versus time with its filtered value.
FIG. 4 is a graph of the change in the measured value of the input yaw and its filtered value over time.
Fig. 5 is a graph showing the change of the measured horizontal velocity value and the filtered value thereof with time.
Fig. 6 is a graph showing the change of the measured vertical velocity value and the filtered value with time.
FIG. 7 is a graph of the change over time of the input lateral velocity actual measurement and its filtered value.
FIG. 8 is a graph of the identified drag compliance coefficient and lift compliance coefficient versus time.
Fig. 9 is a graph of predicted subsequent flight trajectory changes for an uncontrolled aircraft.
Fig. 10 is a graph of predicted subsequent flight trajectory bias current variation for an uncontrolled aircraft.
Fig. 11 is a graph of predicted subsequent flight trajectory speed change for an uncontrolled aircraft.
Detailed Description
Referring to fig. 1, the method for on-line real-time pneumatic parameter identification and ballistic prediction based on a section of measured ballistic parameters comprises the following steps:
first, according to the measured ballistic data And the characteristics of the subsequent flight trajectory, and based on a state equation describing the flight of the uncontrolled aircraft in the outer ballistics, establishing a simplified motion state equation adapting to the online rapid calculation of the uncontrolled aircraft. And establishing an unscented Kalman ballistic filter pneumatic parameter identification model by using an unscented Kalman filtering algorithm.
(1) The influence of the gyroscopic effect on the lateral movement of the uncontrolled aircraft is fully considered, the correction of the actual resistance coincidence coefficient and the lift coincidence coefficient on the trajectory is considered, and a simplified five-degree-of-freedom movement state equation of the uncontrolled aircraft under the nonstandard condition is established, as shown in the formulas (1) to (12).
The auxiliary equations and the expressions of force and moment are as follows:
δ r ＝arccos(v r ξ/v r )
w x ＝-wcos(α W -α N )
w z ＝-wsin(α W -α N )
g 0 ＝9.78034×(1+5.28001×10 -3 ·sin 2 Λ 1 )
c x ＝c x0 (1+k·δ r 2 )＝c x0 +c x2 ·δ r 2 ，c y ≈c′ y ·δ r ，m z ≈m′ z ·δ r
wherein: v is the flight speed of the uncontrolled aircraft; θ a Speed high-low angle (also called ballistic dip); psi phi type 2 Is the velocity direction angle (also known as ballistic deflection angle); omega ξ 、ω η 、ω ζ Is the rotational angular velocity of the uncontrolled aircraft;a vertical axis high-low angle (also called a high-low swing angle); />Is the longitudinal axis direction angle (also called direction swing angle); gamma is the roll angle of the uncontrolled aircraft; x, y and z are the spatial coordinates (horizontal distance, high trajectory and lateral deflection) of the uncontrolled aircraft; v (V) x 、V y 、V z Three components of the flight speed of an uncontrolled aircraft (horizontal speed, vertical speed, lateral speed); t is the time of flight; />Three components in a ballistic coordinate system for all the resultant forces acting on the uncontrolled aircraft; m is M ξ 、M η 、M ζ Three components on the bullet axis coordinate system for all resultant moments acting on the uncontrolled aircraft; m is the mass of the uncontrolled aircraft; a is equatorial moment of inertia; c is polar moment of inertia; delta 2 Is the angle of attack of the direction; delta 1 The attack angle is high and low; b cx Is the resistance coincidence coefficient; b cy Is the lift conforming coefficient; d is the maximum cross-sectional diameter of the uncontrolled aircraft, l is the length of the uncontrolled aircraft, eta is the rifling entanglement, and S is the maximum cross-sectional area of the uncontrolled aircraft; g is gravity acceleration; ρ is the air density, w x Is the longitudinal wind, w z Is crosswind; c x ,c y ,c z The drag coefficient, the lift coefficient and the lateral force coefficient of the uncontrolled aircraft respectively; m is m z ,m′ zz ,m′ y ,m′ xz Respectively a static moment coefficient, a pitching damping moment coefficient derivative, a magnus moment coefficient derivative and a rolling damping moment coefficient derivative of the uncontrolled aircraft; omega shape E Is the rotation angular velocity of the earth; Λ is the local latitude; r is R E Is the average radius of the earth; h is a 0 Is high in altitude;
(2) Establishing a filtering state equation, wherein the state variable quantity of the uncontrolled aircraft is as follows:
the simplified rigid body trajectory model is rewritten into a state space form as
Wherein:
f 5 ＝v cos 2 sinθ a ，f 6 ＝v sinψ 2 ，
f 11 ＝0，f 12 ＝0
(3) Establishing a filter measurement equation
For an uncontrolled aircraft, the ballistic parameters that can be measured are
Z＝[v x v y v z x y zωξ] T (15)
Since the directly measured data are the position and velocity three components of the uncontrolled aircraft, and the ballistic tilt θ in the filtered model a And ballistic deflection angle ψ 2 . Thus, the measurement equation of the ballistic filter system can be expressed as
Wherein: n is the measurement noise.
Since the state variable X is a 12X 1 dimensional matrix, Z is a 7X 1 dimensional matrix and n is a 7X 1 dimensional matrix.
The variance matrix of the measurement noise n is
Wherein: the symbol sigma represents the standard deviation.
(4) Establishing unscented Kalman filtering discrete-time nonlinear system
Given a discrete-time nonlinear system with n state variables, we note that
X k+1 ＝f(X k ,u k ,t k )+w k (18) Wherein: w (w) k ～(0,Q k ) Representing the variance of the process noise, u k Representing the control vector.
The system measurement equation can be written as
y k ＝h(X k ,t k )+v k (19)
Wherein: v k ～(0,R k ) Representing the variance of the measurement noise.
(5) Filter system initialization
Initial value of optimal estimate for given state variable
And gives the initial value of the state variable covariance optimal estimation
(6) Calculation of filter system state variables and covariance thereof
The estimated value and covariance of the state variable can be propagated from one measurement point (point k-1) to the next (point k) using the following steps.
(a) To propagate from step k-1 to step k, the optimal estimate of step k-1 is utilizedSum of covariance->Constructing sigma point as
Wherein: n is the number of state variables.
(b) Transforming the sigma point by using a known nonlinear function f (, i.e. reduced rigid body ballistic system) to obtain
That is to say, willAs an initial value, the reduced rigid body trajectory equation set is integrated, and the integration step length is exactly deltat=t k -t k-1 Obtaining sigma point +.>Variation value of +.>
(c) Obtaining the predicted value of the optimal estimation of the state variable at the kth moment by using unscented transformation
(d) Similarly, the unscented transformation is used to obtain the predicted value of the covariance at the kth moment
(7) Updating of measurement equations for a filtering system
(a) Constructing sigma points using the most recent state variable optimal estimate predictions, with
(b) Nonlinear transformation is carried out on the sigma point by utilizing a measurement equation, thus obtaining
(c) Obtaining approximate mean value of measurement quantity
(d) Solving approximate covariance of measurement quantity
(e) Cross covariance estimation between state variables and measurement quantities
(f) Kalman gain and basic recursion relation
And secondly, inputting a section of actually measured ballistic data into the unscented Kalman ballistic filter pneumatic parameter identification model established in the first step, and identifying the ballistic filter, the resistance coincidence coefficient and the lift coincidence coefficient.
Step 2-1, establishing a nonlinear filtering system by applying the formula (18) and the formula (19);
step 2-2, initializing state variables and covariance matrix elements by using the formulas (20) and (21);
step 2-3, applying equations (22) to (25), propagating the estimated value and covariance of the state variable from one measurement point (kth-1) to the next measurement point (kth);
step 2-4, applying the formulas (26) to (31), updating the measurement equation to obtain the optimal estimated value of each state variable
I.e. the resistance is identified as conforming to the coefficientAnd lift coefficient +.>
Third, the ballistic filter value of the last point is obtainedAnd the identified resistance coincidence coefficient +.>Coefficient of lift>As initial conditions for calculation of subsequent flight trajectory, the established simplified motion state equations (1) - (12) of the uncontrolled aircraft are applied to predict the subsequent flight trajectory, and any point t on the subsequent flight trajectory of the uncontrolled aircraft is obtained j Ballistic information of the place-> Landing point coordinates (X) of uncontrolled aircraft C 、Z C )。
The invention is described in further detail below with reference to examples:
examples
The on-line real-time aerodynamic parameter identification and ballistic prediction method based on a section of measured ballistic parameters comprises the following contents:
1. inputting a section of measured ballistic parameters
Taking 20140114-01 of a section of ballistic data of a certain uncontrolled aircraft measured by a GPS satellite positioning device of 1 st and 14 th 2014 as an example, taking parameters (parameters) of 6 seconds in 12-18 secondsThe measurement interval is 100 milliseconds), specifically including the horizontal distance x k High y of trajectory k Lateral deviation z k Horizontal speedVertical speed +.>And lateral speed->Over time, as shown by the solid black lines in fig. 2-7, respectively. The measured trajectory parameters comprise measurement error noise of the GPS satellite positioning device, and as can be seen in fig. 4 and 7, the measurement error noise data have obvious oscillation, and the other figures are less obvious in oscillation due to larger values and relatively smaller measurement error noise.
2. Actually measured ballistic parameter filtering and resistance coincidence coefficient and lift coincidence coefficient identification calculation
According to the second step of the specific implementation mode, the input 6 seconds of measurement data are utilized to filter the actually measured ballistic parameters, and the parameters are identified, so that the resistance coincidence coefficient and the lift coincidence coefficient of various influences on the comprehensive trajectory are obtained.
The dashed lines in FIGS. 2-7 are filtered values of the input parameters, respectively horizontal distancesBallistic height->Lateral deviation->Horizontal speed->Vertical speed +.>And sideways directionSpeed->Time-dependent changes. The broken line in fig. 2-7 shows that the method performs ballistic filtering on the measured ballistic data, reduces the noise of the measured data, fully considers the theoretical dynamic model of the physical system, predicts the filtered value by using the theoretical model, corrects the filtered value by using the measured result, and improves the accuracy of the state value.
FIG. 8 is a graph of the identified drag compliance coefficient and lift compliance coefficient versus time. The dynamic adjustment between the measured ballistic parameter value and the theoretical predicted value after a period of time can be seen, and the recognized resistance coincidence coefficient and lift coincidence coefficient gradually converge and tend to a certain constant value.
3. Obtaining initial value of trajectory calculated by subsequent flight trajectory
And taking the ballistic filtering value, the resistance coincidence coefficient and the lift coincidence coefficient of the last point in the ballistic filtering and aerodynamic parameter identification process as initial values of subsequent flight trajectory calculation.
The ballistic filter values, drag coincidence coefficients and lift coincidence coefficients for 18 seconds in fig. 2-8 are respectively: horizontal distance filter valueIs 3722.9m, ballistic high filtering value->Is 2004.1m and lateral deviation filtering value ∈>Is-30.0 m, horizontal velocity filtering value +.>190.3m, vertical velocity filter value +.>16.2m, lateral velocity filter value +.>Is-0.1 m, the resistance accords with the coefficient +.>For 0.9389 and lift compliance coefficient +.>1.0011.
4. Subsequent flight trajectory and landing point coordinate prediction of uncontrolled aircraft
And simplifying a motion state equation and calculating an initial value by using the uncontrolled aircraft, calculating a numerical value by using a fourth-order Dragon-Gregorian tower method, forecasting the subsequent flight trajectory, and obtaining landing point coordinate information of the uncontrolled aircraft.
The dashed lines in fig. 9-11 are the trajectory, bias current and velocity profiles that forecast the trajectory of the subsequent flight of the uncontrolled aircraft. Calculating to obtain landing point coordinates (X C 、Z C ) The measured landing point coordinates of the uncontrolled aircraft were (7670.3 m, 9.1 m) for (7658.6 m, 10.7 m). The method has higher resolving precision when the follow-up actual flight trajectory is forecast according to a section of actually measured flight trajectory parameters of the uncontrolled aircraft, and can meet the engineering use requirements.
Claims (7)
1. An on-line real-time characteristic parameter identification and trajectory prediction method for an uncontrolled aircraft based on a section of measured trajectory parameters is characterized by comprising the following steps:
step 1, selecting aerodynamic parameters of an uncontrolled aircraft as characteristic parameters, and establishing a ballistic filtering aerodynamic parameter identification model; firstly, according to the nature of measured trajectory data and the characteristics of the subsequent flight trajectory of an uncontrolled aircraft, a simplified five-degree-of-freedom rigid motion state equation is established, then a unscented Kalman filtering algorithm is applied to conduct linearization and discretization processing on the established simplified rigid motion state equation, and finally a Kalman ballistic filtering pneumatic parameter identification model is established;
step 2, filtering the measured ballistic parameters by using the ballistic filtering aerodynamic parameter identification model established in the step 1, and identifying the resistance coincidence coefficient and the lift coincidence coefficient of the uncontrolled aircraft to obtain initial conditions, resistance and lift coincidence coefficient for calculating the subsequent flight trajectory;
and 3, forecasting the subsequent flight trajectory by using the initial condition and the coincidence coefficient of the trajectory calculation obtained in the step 2, and obtaining the subsequent actual flight trajectory and trajectory drop point information of the uncontrolled aircraft.
2. The method for on-line real-time feature parameter identification and ballistic prediction of an uncontrolled aircraft based on a measured ballistic parameter of claim 1, wherein the measured ballistic data of step 1 includes a corresponding time t i Space coordinates (x) i ，y i ，z i ) And three-dimensional speedWherein i=1, 2, …, n and n are the number of measuring points.
3. The method for on-line real-time feature parameter identification and trajectory prediction of an uncontrolled aircraft based on a measured ballistic parameter according to claim 2, wherein the building of the ballistic filter pneumatic parameter identification model is specifically:
step 1-1, taking an uncontrolled aircraft as an object, establishing a simplified five-degree-of-freedom rigid motion state equation, specifically:
the auxiliary equations and the expressions of force and moment are as follows:
δ r ＝arccos(v rξ /v r )
w x ＝-wcos(α W -α N )
w z ＝-wsin(α W -α N )
g 0 ＝9.78034×(1+5.28001×10 -3 .sin 2 Λ 1 )
c x ＝c x0 (1+k·δ r 2 )＝c x0 +c x2 ·δ r 2 ，c y ≈c′ y ·δ r m z ≈m′ z ·δ r
wherein: v is the flight speed of the uncontrolled aircraft; θ a Is the high and low angle of the speed; psi phi type 2 Is the velocity direction angle; omega ξ 、ω η 、ω ζ Is the rotational angular velocity of the uncontrolled aircraft;a vertical axis is a high-low angle; />Is the angle of the longitudinal axis; gamma is the roll angle of the uncontrolled aircraft; x, y and z are the space coordinates of the uncontrolled aircraft; v (V) x 、V y 、V z Is three components of the flight speed of the uncontrolled aircraft; t is the time of flight;three components in a ballistic coordinate system for all the resultant forces acting on the uncontrolled aircraft; m is M ξ 、M η 、M ζ Three components on the bullet axis coordinate system for all resultant moments acting on the uncontrolled aircraft; m is the mass of the uncontrolled aircraft; a is equatorial moment of inertia; c is polar moment of inertia; delta 2 Is the angle of attack of the direction; delta 1 The attack angle is high and low; b cx Is the resistance coincidence coefficient; b cy Is the lift conforming coefficient; d is the maximum cross-sectional diameter of the uncontrolled aircraft, l is the length of the uncontrolled aircraft, eta is the rifling entanglement, and S is the maximum cross-sectional area of the uncontrolled aircraft; g is gravity acceleration; ρ is the air density, w x Is the longitudinal wind, w z Is crosswind; c x ,c y ,c z The drag coefficient, the lift coefficient and the lateral force coefficient of the uncontrolled aircraft respectively; m is m z ,m′ zz ,m′ y ,m′ xz Respectively a static moment coefficient, a pitching damping moment coefficient derivative, a magnus moment coefficient derivative and a rolling damping moment coefficient derivative of the uncontrolled aircraft; omega shape E Is the rotation angular velocity of the earth; Λ is the local latitude; r is R E Is the average radius of the earth; h is a 0 Is high in altitude;
step 1-2, a filtering state equation is established, and the state variable quantity of the uncontrolled aircraft is as follows:
the simplified rigid body trajectory model is rewritten into a state space form as
Wherein:
f 5 ＝vcosψ 2 sinθ a ，f 6 ＝vsinψ 2 ，
f 11 ＝0，f 12 ＝0；
step 1-3, establishing a filtering measurement equation
For an uncontrolled aircraft, the measurable ballistic parameter is z= [ v x v y v z x y z ω ξ ] T (15)
Since the directly measured data are the position and velocity three components of the uncontrolled aircraft, and the ballistic tilt θ in the filtered model a And ballistic deflection angle ψ 2 Therefore, the measurement equation of the ballistic filter system is expressed as
Wherein: n is measurement noise;
since the state variable X is a 12X 1 dimensional matrix, Z is a 7X 1 dimensional matrix, and n is a 7X 1 dimensional matrix;
the variance matrix of the measurement noise n is
Wherein: symbol sigma represents standard deviation;
step 1-3, establishing a unscented Kalman filtering discrete time nonlinear system
Given a discrete-time nonlinear system with n state variables, we note that
X k+1 ＝f(X k ,u k ,t k )+w k (18)
Wherein: w (w) k ～(0,Q k ) Representing the variance of the process noise, u k Representing a control vector;
the system measurement equation can be written as
y k ＝h(X k ,t k )+v k (19)
Wherein: v k ～(0,R k ) Representing the variance of the measurement noise;
step 1-3, filter System initialization
Initial value of optimal estimate for given state variable
And gives the initial value of the state variable covariance optimal estimation
Step 1-4, calculation of Filter System State variable and covariance thereof
Propagating the estimated value and covariance of the state variable from one measurement point to the next measurement point using the following steps;
(1) To propagate from step k-1 to step k, the optimal estimate of step k-1 is utilizedSum of covariance->Constructing sigma point as
Wherein: n is the number of state variables;
(2) Transforming the sigma point by using a known nonlinear function f (,) to obtain
Will beAs an initial value, the reduced rigid body trajectory equation set is integrated, and the integration step length is exactly deltat=t k -t k-1 Obtaining sigma point +.>Variation value of +.>
(3) Obtaining the predicted value of the optimal estimation of the state variable at the kth moment by using unscented transformation
(4) Similarly, the unscented transformation is used to obtain the predicted value of the covariance at the kth moment
Step 1-5, updating measurement equation of filtering system
(1) Constructing sigma points using the most recent state variable optimal estimate predictions, with
(2) Nonlinear transformation is carried out on the sigma point by utilizing a measurement equation, thus obtaining
(3) Obtaining approximate mean value of measurement quantity
(4) Solving approximate covariance of measurement quantity
(5) Cross covariance estimation between state variables and measurement quantities
(6) Kalman gain and basic recursion relation
4. The method for identifying and predicting the trajectory of the uncontrolled aircraft on line and in real time based on a segment of measured trajectory parameters according to claim 3, wherein in the step 2, the measured trajectory parameters are filtered, and the resistance coincidence coefficient and the lift coincidence coefficient are identified specifically as follows:
step 2-1, establishing a nonlinear filtering system by applying the formula (18) and the formula (19);
step 2-2, initializing state variables and covariance matrix elements by using the formulas (20) and (21);
step 2-3, applying equations (22) to (25), propagating the estimated value and covariance of the state variable from one measurement point to the next measurement point;
step 2-4, applying the formulas (26) to (31), updating the measurement equation to obtain the optimal estimated value of each state variable
I.e. the resistance is identified as conforming to the coefficientAnd lift coefficient +.>
5. The method for identifying and predicting the trajectory of the uncontrolled aircraft on line in real time based on a measured trajectory parameter according to claim 4, wherein the predicting the trajectory of the subsequent flight in step 3 specifically comprises:
ballistic filter value of last pointAnd the identified resistance coincidence coefficient +.>Lift force conforming to coefficient +.>As initial conditions for calculation of subsequent flight trajectory, the established simplified rigid motion state equation (1) of the uncontrolled aircraft is applied to predict subsequent flight trajectory, and any point t on the subsequent flight trajectory of the uncontrolled aircraft is obtained j Ballistic information of the place->And uncontrolled aircraft landing point coordinates (X C 、Z C )。
6. An electronic device comprising a memory, a processor and a computer program stored on the memory and executable on the processor, characterized in that the processor implements the steps of the method according to any of claims 1-5 when the program is executed.
7. A computer readable storage medium, on which a computer program is stored, characterized in that the program, when being executed by a processor, implements the steps of the method according to any of claims 1-5.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310442699.6A CN116611160A (en) | 2023-04-23 | 2023-04-23 | Online real-time characteristic parameter identification and trajectory prediction method for uncontrolled aircraft based on measured trajectory parameters |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310442699.6A CN116611160A (en) | 2023-04-23 | 2023-04-23 | Online real-time characteristic parameter identification and trajectory prediction method for uncontrolled aircraft based on measured trajectory parameters |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN116611160A true CN116611160A (en) | 2023-08-18 |
Family
ID=87680806
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310442699.6A Pending CN116611160A (en) | 2023-04-23 | 2023-04-23 | Online real-time characteristic parameter identification and trajectory prediction method for uncontrolled aircraft based on measured trajectory parameters |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116611160A (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117951922A (en) * | 2024-03-26 | 2024-04-30 | 西安现代控制技术研究所 | Remote guidance rocket online aerodynamic coefficient identification method |
-
2023
- 2023-04-23 CN CN202310442699.6A patent/CN116611160A/en active Pending
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117951922A (en) * | 2024-03-26 | 2024-04-30 | 西安现代控制技术研究所 | Remote guidance rocket online aerodynamic coefficient identification method |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109033493B (en) | Unscented Kalman filtering-based method for identifying pneumatic parameters of high-speed rotation bullet | |
| CN109933847B (en) | Improved active segment trajectory estimation algorithm | |
| CN107367941B (en) | Method for observing attack angle of hypersonic aircraft | |
| CN107423556B (en) | Remote rocket gun emission data calculation method based on radial basis function proxy model | |
| CN108073742B (en) | Method for estimating flight state of intercepted missile tail section based on improved particle filter algorithm | |
| CN116611160A (en) | Online real-time characteristic parameter identification and trajectory prediction method for uncontrolled aircraft based on measured trajectory parameters | |
| CN112711816B (en) | Flight projectile trajectory correction method based on meteorological grid | |
| CN112710298B (en) | Rotating missile geomagnetic satellite combined navigation method based on assistance of dynamic model | |
| CN112013726A (en) | Three-order model-based full strapdown guidance control integrated design method | |
| CN114462293B (en) | Hypersonic speed target medium-long term track prediction method | |
| CN114020019A (en) | Guidance method and device for aircraft | |
| Guan et al. | Modeling of dual-spinning projectile with canard and trajectory filtering | |
| CN117610466A (en) | Shell segment aerodynamic parameter identification method based on model prediction static programming algorithm | |
| CN114001602A (en) | Rocket gun disturbance detection method based on quaternion Kalman filtering denoising fusion | |
| CN109376364B (en) | High-speed rotation bullet pneumatic parameter identification method based on extended Kalman filtering | |
| Gite et al. | Estimation of yaw angle from flight data using extended Kalman filter | |
| CN115542746A (en) | Energy control reentry guidance method and device for hypersonic aircraft | |
| CN115342815A (en) | Estimation method for visual angle rate of maneuvering target in anti-atmosphere or near space | |
| CN114935277A (en) | Online planning method for ideal trajectory of gliding extended-range guided projectile | |
| CN112949150A (en) | Variable structure-based adaptive multi-model box particle filter ballistic target tracking method | |
| CN115046433B (en) | Aircraft time collaborative guidance method based on deep reinforcement learning | |
| Liu et al. | Aerodynamic Parameters Identification of EKF for High-Spinning Projectile Based on Geomagnetic Measurement Data | |
| Fiot et al. | A velocity observer for exterior ballistics using an embedded frequency detection of pitch and yaw aerodynamics | |
| CN113687096B (en) | Crosswind estimation method based on embedded atmospheric data system | |
| CN117519257B (en) | Supersonic speed cruising altitude control method based on back-stepping method |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |