CN116992553A - Whole-course trajectory estimation method of boosting gliding aircraft - Google Patents
Whole-course trajectory estimation method of boosting gliding aircraft Download PDFInfo
- Publication number
- CN116992553A CN116992553A CN202310601814.XA CN202310601814A CN116992553A CN 116992553 A CN116992553 A CN 116992553A CN 202310601814 A CN202310601814 A CN 202310601814A CN 116992553 A CN116992553 A CN 116992553A
- Authority
- CN
- China
- Prior art keywords
- model
- vector
- bgv
- state
- dynamic pressure
- 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.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 43
- 238000005259 measurement Methods 0.000 claims abstract description 32
- 238000001914 filtration Methods 0.000 claims abstract description 13
- 238000013461 design Methods 0.000 claims abstract description 9
- 239000013598 vector Substances 0.000 claims description 76
- 230000001133 acceleration Effects 0.000 claims description 48
- 239000011159 matrix material Substances 0.000 claims description 20
- 230000008569 process Effects 0.000 claims description 13
- 230000008859 change Effects 0.000 claims description 10
- 238000010276 construction Methods 0.000 claims description 10
- 238000004364 calculation method Methods 0.000 claims description 6
- 230000007423 decrease Effects 0.000 claims description 6
- 230000009466 transformation Effects 0.000 claims description 5
- 230000000630 rising effect Effects 0.000 claims description 4
- 238000010924 continuous production Methods 0.000 claims description 2
- 230000007704 transition Effects 0.000 claims description 2
- 238000006243 chemical reaction Methods 0.000 claims 1
- 238000004422 calculation algorithm Methods 0.000 abstract description 13
- 230000003190 augmentative effect Effects 0.000 abstract description 2
- 238000001514 detection method Methods 0.000 description 5
- 238000010586 diagram Methods 0.000 description 4
- 238000012986 modification Methods 0.000 description 3
- 230000004048 modification Effects 0.000 description 3
- 230000005653 Brownian motion process Effects 0.000 description 2
- 230000001174 ascending effect Effects 0.000 description 2
- 230000003416 augmentation Effects 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- 230000008878 coupling Effects 0.000 description 2
- 238000010168 coupling process Methods 0.000 description 2
- 238000005859 coupling reaction Methods 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 238000000342 Monte Carlo simulation Methods 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 238000005311 autocorrelation function Methods 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 230000003111 delayed effect Effects 0.000 description 1
- 230000008034 disappearance Effects 0.000 description 1
- 238000009826 distribution Methods 0.000 description 1
- 230000005484 gravity Effects 0.000 description 1
- 238000003825 pressing Methods 0.000 description 1
- 230000011218 segmentation Effects 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000001052 transient effect Effects 0.000 description 1
- 238000009827 uniform distribution Methods 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/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
- 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
Abstract
The invention discloses a whole course trajectory estimation method of a boosting gliding aircraft, which comprises the following steps: s1, constructing a whole-course trajectory estimation model set comprising a differential polynomial model, an accurate dynamics model and parameterized dynamics based on the whole-course trajectory characteristics of a BGV of a booster gliding aircraft; s2, constructing a variable structure measurement equation based on a low-orbit satellite constellation; and S3, when the BGV is subjected to ballistic estimation based on a variable structure measurement equation, the whole-course ballistic estimation of the BGV is completed based on switching decision and state handover rules of all models in a numerical design model set of real-time dynamic pressure. The invention provides a full-range trajectory estimation method of a booster gliding aircraft, which constructs a model set comprising a differential polynomial model, an accurate dynamics model and parameterized dynamics based on BGV full-range trajectory characteristics, designs a switching decision and a state handover of the model, and further provides a self-adaptive single-model augmented Kalman filtering algorithm for full-range trajectory estimation and parameter identification.
Description
Technical Field
The invention relates to the field of detection of a boosting gliding aircraft. More particularly, the invention relates to a method for estimating the whole course trajectory of a boosted gliding aircraft.
Background
With the increasing maturity of the construction and use of low orbit satellite constellation, it becomes possible to detect the whole course of the assisted gliding aircraft (Boost Glide Vehicle, BGV), but at present, the trajectory estimation algorithm for BGV is mostly based on ground-based radar detection data to estimate a certain segment of trajectory, and no suitable whole course trajectory estimation method is available.
Disclosure of Invention
It is an object of the present invention to address at least the above problems and/or disadvantages and to provide at least the advantages described below.
To achieve these objects and other advantages and in accordance with the purpose of the invention, a method for estimating the global trajectory of a booster gliding aircraft is provided, comprising:
s1, constructing a whole-course trajectory estimation model set comprising a differential polynomial model, an accurate dynamics model and parameterized dynamics based on the whole-course trajectory characteristics of a BGV of a booster gliding aircraft;
s2, constructing a variable structure measurement equation based on a low-orbit satellite constellation;
and S3, when the BGV is subjected to ballistic estimation based on a variable structure measurement equation, the whole-course ballistic estimation of the BGV is completed based on switching decision and state handover rules of all models in a numerical design model set of real-time dynamic pressure.
Preferably, in S1, the differential polynomial model is a uniform acceleration model and is applied to a rising section of BGV, and a construction procedure of the differential polynomial model is configured to include:
in each observation interval, the motion of BGV is set as a uniform acceleration process, and the state quantity X in 9 dimensions of the rising section is as follows:
in the above, the state quantity X includes the position vector [ X Y Z ] of the ground system] T Denoted by r, velocity vectorUse->Indicating acceleration vector +.>Use->A representation;
the 9-dimensional state equation can be expressed as:
wherein O is 3×3 Is a third-order zero square matrix, I 3×3 Is a three-order unit array, tau m For the maneuver time constant, T represents the measurement interval, w (T) is the continuous process noise vector, when T/τ m When < 1, the discretized form is:
in E { w } k }＝O 3×1 ，，Phi represents a state transition matrix, X k Representing the moment kState vector, G 1 Represented by a process noise matrix, O 3×1 Represented is a zero matrix of 3 rows and 1 column, W k Represented is a discretized driving white noise sequence for identifying parameters,/>Is to W k Transpose of vector matrix is performed, E { } is the expectation of random variables in brackets,/>The square of the maximum acceleration is shown.
Preferably, the precise kinetic model is applied to the free segment of the BGV, and the construction flow is configured to include:
constructing and taking into account the non-spherical attraction J of the earth when BGV flies to a height where the dynamic pressure is smaller than a predetermined value 2 The following accurate kinetic model of perturbation:
wherein mu is the gravitational constant, a e Is the radius of the equator of the earth,radial of the center of the earthLatitude of the earth's center->ω e Is the rotation angular velocity of the earth;
defining a 6-dimensional state vector:
discretizing the accurate kinetic model into:
a in the above k (X k ) Substituting the state vector at the moment k into the constant acceleration in the observation interval obtained by calculation of the accurate dynamics model (9); e { W k }＝O 3×1 G 2 The process noise matrix of the accurate dynamics model is represented, C is a constant, and the value 7,u is typically a determined acceleration increment.
Preferably, the parametric dynamics model is applied to the gliding segment of the BGV, and the construction process is configured to include:
after the BGV height is reduced and the dynamic pressure is larger than a preset value, combining the aerodynamic coefficient, the windward cross-sectional area and the mass into aerodynamic parameters to be jointly estimated together with the motion state of the BGV, and constructing the following parameterized dynamic model:
wherein alpha represents an aerodynamic drag parameter, andm represents mass, S represents cross-sectional area, C d The aerodynamic drag coefficient, K the lift-drag ratio, gamma the roll angle, ρ the atmospheric density, ν the velocity, +.>The coordinate transformation matrix from the speed coordinate system to the geocentric fixedly connected coordinate system is shown; the method comprises the steps of carrying out a first treatment on the surface of the
Defining a corresponding 9-dimensional state vector:
discretizing the parameterized kinetic model into a sectional uniform acceleration model:
in which a is k (X k ) Substituting the state vector at the moment k into the constant acceleration in the observation interval obtained by calculation of the parameterized dynamics model (12);
wherein alpha is max Represents the maximum value of aerodynamic resistance parameter, K max The lift-drag ratio maximum value is represented, and Δγ represents the roll angle variation range.
Preferably, the construction flow of the variable structure measurement equation is configured to include:
because of the measurement vector dimension of the measurement equation and the constant change of the orbit parameters of the visible satellite, the measurement equation is represented by n>Observed quantity of 1 satellite (A) j ，E j ) And target state X k+1 The relationship between them can be expressed as:
Z k+1 ((A j ，E j ) j＝1，…，n )h(X k+1 )
let the number of orbits of the j-th visible satellite be (a, e, i, omega, f) j The position vectors of the target in the satellite orbit system are respectively:
in the middle ofRepresented is a coordinate transformation matrix of the ground system to the rail system,/for>Representing the position vector of the satellite on the ground system;
calculating two angular measurements of the satellite infrared detector from the position vector: defining azimuth a as the position vector at x s oy s In-plane projection and oy s The angle of the axis, elevation E, is the position vector and x s oy s The relationship between the observation vector and the position vector of the missile in the satellite orbit system is as follows:
preferably, the switching decision and state handover rule of each model in the model set is as follows:
the use sequence and the switching of the model are determined based on the numerical value of the real-time dynamic pressure, and when the ground-fixed system position vector and the speed vector at the k moment are known, the dynamic pressure calculation method comprises the following formula:
in the above formula, q represents dynamic pressure, L represents geodetic longitude, B represents geodetic latitude, and H represents geodetic altitude;
by combining aerodynamic acceleration with earth J 2 Dynamic pressure with equivalent magnitude of acceleration change caused by perturbation is used as threshold value q of model switching T ；
Starting from rocket engine ignition, the dynamic pressure is increased and then reduced, when the dynamic pressure is gradually reduced to q T Previously, a differential polynomial model was used;
when the dynamic pressure gradually decreases to q T When the model is switched to the accurate dynamics model;
when movingThe pressure decreases and increases again, when it increases to q T At that time, the parameterized kinetic model is switched to and used all the way thereafter.
Preferably, when the differential polynomial model is applied, the satellite measurement data of the first 3 moments need to be filtered and initialized, and then the filtering operation is performed;
when the differential polynomial model performs state handover to the accurate dynamics model, the following formula is adopted to complete handover, and filtering operation is performed after handover:
when the accurate dynamics model performs state handover to the parameterized dynamics model, the following formula is adopted to expand and initialize the state to complete handover, and filtering operation is performed after handover:
the invention at least comprises the following beneficial effects: based on low orbit satellite constellation detection data, a global trajectory estimation and identification algorithm for BGV is provided: based on the global trajectory characteristics of the BGV, a model set comprising a differential polynomial model, an accurate dynamics model and parameterized dynamics is constructed, and a switching decision and a state handover of the model are designed, so that a self-adaptive single-model augmentation Kalman filtering algorithm is provided and used for global trajectory estimation and parameter identification.
Additional advantages, objects, and features of the invention will be set forth in part in the description which follows and in part will become apparent to those having ordinary skill in the art upon examination of the following or may be learned from practice of the invention.
Drawings
FIG. 1 is a flow chart of a model use and switch structure of the present invention;
FIG. 2 is a schematic diagram of a BGV ballistic curve;
FIG. 3 is a diagram of the overall ballistic visible satellite number distribution;
FIG. 4 is a graph showing the comparison of the root mean square error curve of trajectory estimation obtained by the method of the present invention with the position estimation based on the result of the even acceleration single model algorithm;
FIG. 5 is a graph showing the comparison of the root mean square error curve of trajectory estimation obtained by the method of the present invention with the speed estimation based on the results of the homogeneous acceleration single model algorithm; the method comprises the steps of carrying out a first treatment on the surface of the
FIG. 6 is a schematic diagram of the result of estimating the acceleration of the boost segment obtained by the present invention;
fig. 7 is a schematic diagram of a result of estimating pneumatic parameters of a boost segment obtained by adopting the method.
Detailed Description
The present invention is described in further detail below with reference to the drawings to enable those skilled in the art to practice the invention by referring to the description.
Model set design for 1 whole course trajectory estimation
The trajectory estimation state vector is selected to be described in the ground-based system, i.e. the position vector of the ground-based system [ X Y Z ] taking into account the geometrical relationship of the detector and the object] T And velocity vectorThe differential polynomial model is added with the earth-fixed acceleration vector, and the parameterized dynamic model is added with the aerodynamic parameter vector.
Whether it is acceleration in a differential polynomial model or aerodynamic parameters in a parametric dynamics model, it is essentially an augmented kalman filter that recognizes and estimates model parameters while state estimation.
1.1 design of parameter model to be identified
The conventional state augmentation method generally uses Wiener process to describe parameters to be estimated, namely, supposing that increment between parameter sequences to be estimated is uncorrelated, actually, for BGV, whether boost section acceleration or aerodynamic parameters are used as control instructions, the control instructions generally have non-randomness and piecewise correlation, the parameters to be estimated are described to be more consistent with the reality by using first-order Markov process with autocorrelation function in exponential decay form, and the corresponding parameter model to be estimated can be expressed as follows
τ in m W (t) is the drive white noise, which is the maneuver time constant. Discretizing the formula (1) (1) to obtain
(θ i ) k+1 ＝(θ i ) k exp(-T/τ m )+w k (2)
W in k For discretizing the driving white noise sequence for identifying parameters, the variance is
In the middle ofIs the variance of the parameter to be estimated. When T/tau m In the case of < 1, the formulae (3) (3) can be simplified to
I.e. noise variance and parameter varianceProportional to the measurement interval T, to the maneuvering time constant τ m Inversely proportional. />Reflecting the variation range of the parameter, representing the maneuverability of the target, τ m Then reflecting the real-time maneuver strategy of the target, different τ m Corresponding to different real-time motion patterns.
From equations (3) (3) and (4) (4), it is known that the first order Markov process model of the parameter to be estimated will drive the variance Q of the white noise as compared with the Wiener process model θ Expressed as a target maneuver time constant τ m And parameter variance to be estimatedThe function of (2) has clearer physical meaning, and is beneficial to the construction and self-adaptive design of noise variance.
For BGV, the acceleration and aerodynamic parameters of the boosting section are usually limited in value, so that the variation range is determined, and the value is assumed to be within a certain limited interval theta epsilon theta min ,θ max ]And satisfies the uniform distribution in the interval, the variance is
After the parameter model to be identified is provided, a state equation for trajectory estimation of each flight stage can be constructed.
1.2 differential polynomial model
The BGV is acted by thrust, aerodynamic force and earth attraction force when flying in the ascending section, the thrust characteristics (thrust curve, second consumption and the like), flying program and aerodynamic force parameters and the like of the defending side are unknown, and the assumption of gravity turning of the boosting section is not established any more, so the pitch angle is formedMass m, aerodynamic parameters-> And->Is a time-varying unknown parameter, and the thrust P and cross-sectional area S are approximately constant but unknown. The coupling of multiple unknown parameters cannot be accurately estimated one by one based on satellite detection data only, but only the combination of the coupling parameters (namely acceleration along three directions of coordinate axes) can be estimated, and the dynamics model is converted into a differential polynomial model. The invention adopts a uniform acceleration model, namely, in each observation interval, the motion of BGV is assumed to be a uniform acceleration process. The state quantity is 9 dimensions of ground fixation system position, speed and acceleration
The 9-dimensional state equation can be expressed as
In which r= [ X Y Z] T Is the position vector in the earth-fixed system, O 3×3 Is a third-order zero square matrix, I 3×3 Is a three-order unit array. When T tau m When < 1, the discretized form is
In E { w } k }＝O 3×1 ，
1.3 accurate kinetic model
When BGV flies to a height where dynamic pressure is smaller than a certain value, aerodynamic force cannot be used for control at the moment, and the method is similar to a free section of a ballistic missile, and can construct and consider J 2 Accurate dynamic model of perturbation
Discretization is also required for formulas (9) (9). At this time, a state vector is defined
Discretized state model is
In which a is k (X k ) To get the state vector at k timeSubstituting the amount into the constant acceleration in the observation interval calculated by the formulas (9) and (9); e { w k }＝O 3×1 ，C is a constant, typically a value of 7, and the unknown input is no greater than J 2 Perturbation acceleration, set u=0.01 m/s 2 。
1.4 parameterized kinetic model
When the dynamic pressure is larger than a given value, the stage of only utilizing aerodynamic force to carry out maneuvering control is carried out, and the disappearance of thrust force reduces a large number of parameters to be estimated, so that the aerodynamic coefficient, the windward cross section area and the mass can be combined into a parameter (called aerodynamic parameter) to be jointly estimated together with the motion state of the BGV, and a parameterized dynamic model is formed.
In the middle ofThe pneumatic resistance parameter is that K is lift-drag ratio and gamma is roll angle.
Equation (12) (12) is a nonlinear continuous system, and it is necessary to discretize it. Under the condition of small observation interval, the method can adopt a segmentation uniform acceleration model to discretize, so that the solution of the Jacobi matrix is avoided. Defining a state vector
Discretized state model is
In which a is k (X k ) For substituting the state vector at the k moment into the constant acceleration (the discretized state equation still has strong nonlinearity due to the nonlinearity relation between the acceleration and the state vector) in the observation interval calculated by the formulas (121) (12), E { w } k }＝O 3×1 ，
The differential polynomial model for the ascending segment, the precise kinetic model for the free segment and the parameterized kinetic model for the gliding segment together form a global trajectory estimation model set.
2 variable structure measurement equation based on low orbit satellite constellation
With reference to orbital parameters of the satellites, the low-orbit constellation designed for detection is 225/15/0: i.e. 15 orbital planes, 15 satellites in each orbital plane, the orbital tilt angle adopts 97.6deg near the polar orbit, the satellites adopt circular orbits, and the orbit height is 560km.
The variable structure here means that the number and the number of the visible satellites change frequently in one ballistic cycle due to the severe change of the relative position relationship between the BGV and the low-orbit satellites, so that the dimension of the measurement vector of the measurement equation and the orbit parameters of the visible satellites change continuously, i.e. the measurement equation is variable in structure. The observation equation is n (n>1) Observed quantity of visible satellite (A j ,E j ) j＝1,…,n And target state X k+1 Relationship between them.
Z k+1 ((A j ，E j ) j＝1，…，n )h(X k+1 ) (15)
For simplicity, the following formula omits subscript k+1, assuming the number of orbits of the j-th visible satellite is (a, e, i, Ω, ω, f) j Firstly, converting the satellite position vector into a ground inertial system position and velocity vector, and then converting the satellite position vector into a ground fixed system position and velocity vector, wherein the satellite position vector in the ground fixed system isThe position vector of the target in the ground system is [ X Y Z ]] T The position vectors of the targets in the satellite orbit system are respectively
In the middle ofIs a coordinate transformation matrix of a ground system to a track system.
From the position vector, two angular measurements of the satellite infrared detector can be calculated: defining azimuth a as the position vector at x s oy s In-plane projection and oy s The angle of the axis, elevation E, is the position vector and x s oy s Included angle of plane. The relationship between the observation vector and the position vector of the missile in the satellite orbit system is that
3 model switching decision and state handoff design
After the measurement equations of the state model set and the variable structure are constructed, the use sequence and the switching decision of the models are also required to be determined, and the state of the different models is switched (or called the reinitialization of the new model).
3.1 model switching decisions
Determining the use sequence and switching of the model according to the numerical value of the real-time dynamic pressure, and when the state estimation (ground-system position vector and speed vector) at the moment k is known, the dynamic pressure calculation method at the moment is as follows
Determining a dynamic pressure q threshold value: general aerodynamic parametersMagnitude of about 10 -5 m 2 Kg, taking into account the earth J 2 The change in acceleration due to perturbation is about 0.01m/s 2 Definition of the pneumatic acceleration and the Earth J 2 Dynamic pressure with equivalent magnitude of acceleration change caused by perturbation is the threshold value q of model switching T 。
Starting from rocket engine ignition, the dynamic pressure is increased and then reduced, when the dynamic pressure is gradually reduced to q T Before, a uniform acceleration model is adopted; when the dynamic pressure gradually decreases to q T When the model is switched to the accurate dynamics model; when the dynamic pressure decreases and increases, when the dynamic pressure increases to q T At this time, the parameterized kinetic model is switched to and used all the time thereafter. The model usage and switching architecture is shown in fig. 1.
3.2 Filter initialization and State handoff design
The differential polynomial model adopts a uniform acceleration model, and the state vector is 9 dimensions of position, speed and acceleration, so that satellite measurement data at the first 3 moments are required to be filtered and initialized. The initialization process is divided into two steps: (a) And converting the double-star observed data into equivalent position observed data. (b) And obtaining an initial state vector estimation and a covariance matrix by using the weighted least square estimation. Specific method reference publications are not repeated.
The state vector of the differential polynomial model includes 9 dimensions of position, velocity and acceleration, while the state vector of the exact kinetic model is 6 dimensions of position and velocity, so that the following state handoffs can be made
The state vector of the accurate dynamics model is 6-dimensional of the ground-based position vector and the speed vector, and the state vector of the parameterized dynamics model is the ground-based position vector, the speed vector and aerodynamic parameters, so that state expansion and initialization are required during state handover, and alpha is caused 0 ＝10 -5 m 2 /kg，K 0 ＝3，γ 0 ＝0，α max ＝10 -4 m 2 /kg，K max ＝3，Δγ＝10°。
The propagation of the state estimation mean and covariance adopts unscented transformation, a filter adopts unscented Kalman filtering, and a specific algorithm refers to relevant documents.
4 simulation result analysis
A BGV trajectory is designed for estimating the estimation precision of the full trajectory estimation algorithm, and a trajectory curve is shown in figure 2.
The trajectory is detected by using a satellite constellation constructed in a variable structure measurement equation based on a low orbit satellite constellation, and the number of visible satellites for obtaining the whole trajectory is shown in figure 3. It can be seen that the constructed satellite constellation can meet the double coverage of the full trajectory of BGV, with a minimum number of visible satellites of 2 and a maximum of 6. In the flight process of approximately 900sec in BGV, the number of satellite handover times is more than 18, and the average observed arc length of each satellite is about 50 seconds, which puts high demands on the target tracking capability of the satellite and the data transmission processing capability of the target handover.
Assume that the infrared detector has an angular error of 10 -4 rad, generating simulated measurement data of the trajectory, performing full trajectory estimation of BGV by adopting the method, obtaining a root mean square error curve of the trajectory estimation by Monte Carlo simulation, and comparing the root mean square error curve with a result based on a uniform acceleration single model algorithm as shown in figures 4-5.
From the results, it can be seen that:
(1) Since the number of the visible satellites is equal to or greater than 2, the three-dimensional positioning can be formed on the target, so that the filtering is convergent, the stable tracking can be formed on the target, and the tracking estimation accuracy of the target depends on the maneuverability of the target and the number of the visible satellites. The stronger the target mobility, the fewer the number of satellites in view, the greater the estimation error and the lower the accuracy.
(2) In the transient boosting section, the estimation accuracy of the two algorithms is equivalent, but in the free section and the pneumatic control section, the algorithm is significantly smaller than the uniform acceleration model method, when the filtering reaches a steady state, the position estimation error is reduced from 50m to 20m, and the speed estimation error is reduced from 20m/s to about 5m/s, because more correct target dynamics information is introduced in the two stages, and the model is more finely matched.
(3) The algorithm herein is convergence stable even when the model is switched, its estimation error is not significantly increased and is smaller than the uniformly accelerated model approach, which illustrates that the model switching decisions presented herein are accurate and efficient.
(4) In the boosting section, the pulling section and the pressing section, the filtering error is larger due to the severe acceleration change, but the method can still ensure the stable filtering convergence, and if the estimation accuracy is required to be improved, the measurement data rate can be increased or the measurement data accuracy can be improved.
The method not only can obtain the position and speed estimation result with higher precision, but also can further obtain the acceleration information of the boosting section and the pneumatic parameter information of the pneumatic control section, and the comparison with the true value result is shown in fig. 6 and 7.
From the results, it can be seen that: the algorithm can accurately and effectively identify the model parameters, and in the boosting section, even if the acceleration changes severely, the estimation result can still effectively track the true value (slightly delayed); in the pneumatic control section, pneumatic resistance parameters, lift-drag ratio and roll angle can be effectively identified, and the pneumatic control section can be used for pneumatic characteristic analysis of targets and subsequent rapid prediction of reachable areas and tracks.
The above is merely illustrative of a preferred embodiment, but is not limited thereto. In practicing the present invention, appropriate substitutions and/or modifications may be made according to the needs of the user.
The number of equipment and the scale of processing described herein are intended to simplify the description of the present invention. Applications, modifications and variations of the present invention will be readily apparent to those skilled in the art.
Although embodiments of the invention have been disclosed above, they are not limited to the use listed in the specification and embodiments. It can be applied to various fields suitable for the present invention. Additional modifications will readily occur to those skilled in the art. Therefore, the invention is not to be limited to the specific details and illustrations shown and described herein, without departing from the general concepts defined in the claims and their equivalents.
Claims (7)
1. A method for estimating the overall trajectory of a boosted gliding aircraft, comprising:
s1, constructing a whole-course trajectory estimation model set comprising a differential polynomial model, an accurate dynamics model and parameterized dynamics based on the whole-course trajectory characteristics of a BGV of a booster gliding aircraft;
s2, constructing a variable structure measurement equation based on a low-orbit satellite constellation;
and S3, when the BGV is subjected to ballistic estimation based on a variable structure measurement equation, the whole-course ballistic estimation of the BGV is completed based on switching decision and state handover rules of all models in a numerical design model set of real-time dynamic pressure.
2. The method for estimating global trajectory of a booster gliding aircraft according to claim 1, wherein in S1, the differential polynomial model is a uniform acceleration model and is applied to a rising segment of BGV, and a construction process of the differential polynomial model is configured to include:
in each observation interval, the motion of BGV is set as a uniform acceleration process, and the state quantity X in 9 dimensions of the rising section is as follows:
in the above, the state quantity X includes the position vector [ X Y Z ] of the ground system] T Denoted by r, velocity vectorUse->Indicating acceleration vector +.>Use->A representation;
the 9-dimensional state equation can be expressed as:
wherein O is 3×3 Is a third-order zero square matrix, I 3×3 Is a three-order unit array, tau m For the maneuver time constant, T represents the measurement interval, w (T) is the continuous process noise vector, when T/τ m When < 1, the discretized form is:
in the middle ofPhi represents a state transition matrix, X k Representing the state vector at time k, G 1 Represented by a process noise matrix, O 3×1 Represented is a zero matrix of 3 rows and 1 column, W k Represented is a discretized driving white noise sequence for identifying parameters,/>Is to W k Transpose of vector matrix is performed, E { } is the expectation of random variables in brackets,/>The square of the maximum acceleration is shown.
3. The method of whole course trajectory estimation of a boosted gliding aircraft according to claim 1, wherein the accurate dynamics model is applied to the free segment of BGV, the construction procedure being configured to include:
constructing and taking into account the non-spherical attraction J of the earth when BGV flies to a height where the dynamic pressure is smaller than a predetermined value 2 The following accurate kinetic model of term perturbation:
wherein mu is the gravitational constant, a e Is the radius of the equator of the earth,radial of the center of the earthLatitude of the earth's center->ω e Is the rotation angular velocity of the earth;
defining a 6-dimensional state vector:
discretizing the accurate kinetic model into:
a in the above k (X k ) Constant in the observation interval calculated for substituting the state vector at time k into the accurate kinetic model (9)Acceleration of (2);G 2 the process noise matrix of the accurate dynamics model is represented, C is a constant, and the value 7,u is typically a determined acceleration increment.
4. The method of whole course trajectory estimation of a booster gliding aircraft according to claim 1, wherein the parameterized kinetic model is applied to a gliding segment of a BGV, and the construction flow is configured to include:
after the BGV height is reduced and the dynamic pressure is larger than a preset value, combining the aerodynamic coefficient, the windward cross-sectional area and the mass into aerodynamic parameters to be jointly estimated together with the motion state of the BGV, and constructing the following parameterized dynamic model:
wherein alpha represents an aerodynamic drag parameter, andm represents mass, S represents cross-sectional area, C d The aerodynamic drag coefficient, K the lift-drag ratio, gamma the roll angle, ρ the atmospheric density, ν the velocity, +.>Representing the speed coordinate system to the earth's centerA coordinate conversion matrix of a fixed coordinate system;
defining a corresponding 9-dimensional state vector:
discretizing the parameterized kinetic model into a sectional uniform acceleration model:
in which a is k (X k ) Substituting the state vector at the moment k into the constant acceleration in the observation interval obtained by calculation of the parameterized dynamics model (12);
wherein alpha is max Represents the maximum value of aerodynamic resistance parameter, K max The lift-drag ratio maximum value is represented, and Δγ represents the roll angle variation range.
5. The method of whole course trajectory estimation of a boosted gliding aircraft according to claim 1, wherein the construction flow of the variable structure measurement equation is configured to include:
because of the measurement vector dimension of the measurement equation and the constant change of the orbit parameters of the visible satellite, the measurement equation is represented by n>Observed quantity of 1 satellite (A) j ，E j ) And target state X k+1 The relationship between them can be expressed as:
let the number of orbits of the j-th visible satellite be (a, e, i, omega, f) j Position vector of target in satellite orbit systemThe method comprises the following steps of:
in the middle ofRepresented is a coordinate transformation matrix of the ground system to the rail system,/for>Representing the position vector of the satellite on the ground system;
calculating two angular measurements of the satellite infrared detector from the position vector: defining azimuth a as the position vector at x s oy s In-plane projection and oy s The angle of the axis, elevation E, is the position vector and x s oy s The relationship between the observation vector and the position vector of the missile in the satellite orbit system is as follows:
6. the method of estimating the global trajectory of a boosted gliding aircraft according to claim 1, wherein the switching decision and state handover rules of each model in the set of models are:
the use sequence and the switching of the model are determined based on the numerical value of the real-time dynamic pressure, and when the ground-fixed system position vector and the speed vector at the k moment are known, the dynamic pressure calculation method comprises the following formula:
in the above formula, q represents dynamic pressure, L represents geodetic longitude, B represents geodetic latitude, and H represents geodetic altitude;
by combining aerodynamic acceleration with earth J 2 Dynamic pressure with equivalent magnitude of acceleration change caused by perturbation is used as threshold value q of model switching T ；
Starting from rocket engine ignition, the dynamic pressure is increased and then reduced, when the dynamic pressure is gradually reduced to q T Previously, a differential polynomial model was used;
when the dynamic pressure gradually decreases to q T When the model is switched to the accurate dynamics model;
when the dynamic pressure decreases and increases, when the dynamic pressure increases to q T At that time, the parameterized kinetic model is switched to and used all the way thereafter.
7. The method for estimating global trajectory of a booster gliding aircraft according to claim 6, wherein when the differential polynomial model is applied, the satellite measurement data of the first 3 moments is required to be filtered and initialized, and then the filtering operation is performed;
when the differential polynomial model performs state handover to the accurate dynamics model, the following formula is adopted to complete handover, and filtering operation is performed after handover:
when the accurate dynamics model performs state handover to the parameterized dynamics model, the following formula is adopted to expand and initialize the state to complete handover, and filtering operation is performed after handover:
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310601814.XA CN116992553B (en) | 2023-05-25 | 2023-05-25 | Whole-course trajectory estimation method of boosting gliding aircraft |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310601814.XA CN116992553B (en) | 2023-05-25 | 2023-05-25 | Whole-course trajectory estimation method of boosting gliding aircraft |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116992553A true CN116992553A (en) | 2023-11-03 |
| CN116992553B CN116992553B (en) | 2024-02-06 |
Family
ID=88525533
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310601814.XA Active CN116992553B (en) | 2023-05-25 | 2023-05-25 | Whole-course trajectory estimation method of boosting gliding aircraft |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116992553B (en) |
Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20120316819A1 (en) * | 2008-04-22 | 2012-12-13 | United States Government, As Represented By The Secretary Of The Navy | Process for estimation of ballistic missile boost state |
| CN109781374A (en) * | 2018-12-05 | 2019-05-21 | 湖南涉外经济学院 | A kind of method that real-time online quickly estimates aircraft thrust |
| CN109933847A (en) * | 2019-01-30 | 2019-06-25 | 中国人民解放军战略支援部队信息工程大学 | A kind of improved boost phase trajectory algorithm for estimating |
| US10776450B1 (en) * | 2016-09-26 | 2020-09-15 | United States Of America, As Represented By The Secretary Of The Navy | Closed form estimator for ballistic missile flight |
| CN112861250A (en) * | 2020-12-21 | 2021-05-28 | 中国人民解放军96901部队23分队 | Gliding trajectory degradation solution along with energy change based on attack angle and resistance acceleration |
| CN113759956A (en) * | 2020-12-14 | 2021-12-07 | 北京天兵科技有限公司 | Flight trajectory design method for sub-orbital vehicle |
| US20220107160A1 (en) * | 2020-10-02 | 2022-04-07 | United States Of America, As Represented By The Secretary Of The Navy | Glide Trajectory Optimization for Aerospace Vehicles |
| CN114707317A (en) * | 2022-03-25 | 2022-07-05 | 惠州学院 | Method and system for measuring flight parameters of spinning projectile based on trajectory prior knowledge |
-
2023
- 2023-05-25 CN CN202310601814.XA patent/CN116992553B/en active Active
Patent Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20120316819A1 (en) * | 2008-04-22 | 2012-12-13 | United States Government, As Represented By The Secretary Of The Navy | Process for estimation of ballistic missile boost state |
| US10776450B1 (en) * | 2016-09-26 | 2020-09-15 | United States Of America, As Represented By The Secretary Of The Navy | Closed form estimator for ballistic missile flight |
| CN109781374A (en) * | 2018-12-05 | 2019-05-21 | 湖南涉外经济学院 | A kind of method that real-time online quickly estimates aircraft thrust |
| CN109933847A (en) * | 2019-01-30 | 2019-06-25 | 中国人民解放军战略支援部队信息工程大学 | A kind of improved boost phase trajectory algorithm for estimating |
| US20220107160A1 (en) * | 2020-10-02 | 2022-04-07 | United States Of America, As Represented By The Secretary Of The Navy | Glide Trajectory Optimization for Aerospace Vehicles |
| CN113759956A (en) * | 2020-12-14 | 2021-12-07 | 北京天兵科技有限公司 | Flight trajectory design method for sub-orbital vehicle |
| CN112861250A (en) * | 2020-12-21 | 2021-05-28 | 中国人民解放军96901部队23分队 | Gliding trajectory degradation solution along with energy change based on attack angle and resistance acceleration |
| CN114707317A (en) * | 2022-03-25 | 2022-07-05 | 惠州学院 | Method and system for measuring flight parameters of spinning projectile based on trajectory prior knowledge |
Non-Patent Citations (2)
| Title |
|---|
| 吴楠, 王锋, 孟凡坤: "无再入观测弹道导弹气动参数和落点联合预报", 《弹道学报》 * |
| 吴楠: "助推滑翔飞行器预警探测滤波方法与误差链研究", 《中国博士学位论文全文数据库（电子期刊）》 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116992553B (en) | 2024-02-06 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US8639476B2 (en) | Process for estimation of ballistic missile boost state | |
| CN107544067B (en) | Hypersonic reentry vehicle tracking method based on Gaussian mixture approximation | |
| CN109933847B (en) | Improved active segment trajectory estimation algorithm | |
| CN111474960B (en) | Critical altitude supersonic speed target tracking method based on control quantity characteristic assistance | |
| CN104567880A (en) | Mars ultimate approach segment autonomous navigation method based on multi-source information fusion | |
| Chen et al. | Hypersonic boost–glide vehicle strapdown inertial navigation system/global positioning system algorithm in a launch-centered earth-fixed frame | |
| CN110412868B (en) | Non-cooperative spacecraft orbit determination method using inter-satellite optical images | |
| CN109740209A (en) | Hypersonic aircraft on-line parameter identification method and the mechanical model for using it | |
| CN112710298B (en) | Rotating missile geomagnetic satellite combined navigation method based on assistance of dynamic model | |
| CN111121770B (en) | Interactive multi-missile multi-model flight path fusion method | |
| CN114510076A (en) | Target collaborative detection and guidance integrated method and system based on unscented transformation | |
| CN110146092B (en) | Double-body asteroid detection track optimization method based on navigation information evaluation | |
| CN114462293B (en) | Hypersonic speed target medium-long term track prediction method | |
| CN104729510A (en) | Method for determining relative adjoint orbit of space target | |
| CN111598232A (en) | Method for estimating complex micro-motion space cone target parameters by using deep learning convolutional neural network | |
| CN109781374A (en) | A kind of method that real-time online quickly estimates aircraft thrust | |
| CN111409865B (en) | Deep space probe approach segment guidance method based on intersection probability | |
| CN109186614A (en) | Short distance autonomous relative navigation method between a kind of spacecraft | |
| CN110728026B (en) | Terminal trajectory target passive tracking method based on angular velocity measurement | |
| CN116992553B (en) | Whole-course trajectory estimation method of boosting gliding aircraft | |
| Wang et al. | Novel in-flight coarse alignment of low-cost strapdown inertial navigation system for unmanned aerial vehicle applications | |
| Li et al. | Tracklet-to-object Matching for Climbing Starlink Satellites through Recursive Orbit Determination and Prediction | |
| CN114662285A (en) | Intelligent resolving method for fire control model of high-speed aircraft | |
| CN114660587A (en) | Jump and glide trajectory target tracking method and system based on Jerk model | |
| Kim et al. | Ballistic object trajectory and launch point estimation from radar measurements using long-short term memory networks |
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 | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |