CN104778376B - A kind of hypersonic gliding bullet Skipping trajectory Forecasting Methodology of near space - Google Patents

A kind of hypersonic gliding bullet Skipping trajectory Forecasting Methodology of near space Download PDF

Info

Publication number
CN104778376B
CN104778376B CN201510220893.5A CN201510220893A CN104778376B CN 104778376 B CN104778376 B CN 104778376B CN 201510220893 A CN201510220893 A CN 201510220893A CN 104778376 B CN104778376 B CN 104778376B
Authority
CN
China
Prior art keywords
mrow
msub
mtd
mtr
mtable
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.)
Active
Application number
CN201510220893.5A
Other languages
Chinese (zh)
Other versions
CN104778376A (en
Inventor
周荻
秦雷
李君龙
邹昕光
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Harbin Institute of Technology
Original Assignee
Harbin Institute of Technology
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Harbin Institute of Technology filed Critical Harbin Institute of Technology
Priority to CN201510220893.5A priority Critical patent/CN104778376B/en
Publication of CN104778376A publication Critical patent/CN104778376A/en
Application granted granted Critical
Publication of CN104778376B publication Critical patent/CN104778376B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • 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

  • Pharmaceuticals Containing Other Organic And Inorganic Compounds (AREA)
  • Navigation (AREA)

Abstract

A kind of hypersonic gliding bullet Skipping trajectory Forecasting Methodology of near space, the present invention relates to a kind of Flight Vehicle Trajectory Forecasting Methodology.The present invention is to solve existing method it is low to maneuvering target ballistic prediction precision the problem of, and provide a kind of near space it is hypersonic gliding bullet Skipping trajectory Forecasting Methodology.It is realized in the steps below：First, the trajectory of hypersonic gliding bullet is established；2nd, the Kalman filter of the hypersonic gliding bullet movement locus of real-time tracking is designed；3rd, according to position, speed and the acceleration of the hypersonic gliding bullet of tracking finish time, the angle of attack and roll angle during gliding bullet flight hypersonic with reference to trajectory estimation, the hypersonic gliding bullet of near space carries out waiting the angle of attack and waits roll angle flight within the subsequent flight time, recurrence calculation is circulated to subsequent time using trajectory, obtains the ballistic prediction value of the hypersonic gliding bullet after certain time.Belong to target following technical field.

Description

A kind of hypersonic gliding bullet Skipping trajectory Forecasting Methodology of near space
Technical field
It is particularly a kind of for the hypersonic gliding bullet jump of near space the present invention relates to a kind of ballistic prediction method The Forecasting Methodology of trajectory.
Background technology
The hypersonic gliding bullet of near space is that one kind can be thousands of to more than 10,000 public affairs in the unpowered gliding of near space In, and with stronger lateral maneuvering capability and anti-performance of dashing forward, can carry nuclear warheads or conventional warhead is implemented quickly to beat at a distance Hit, this kind of aircraft with higher lift-drag ratio, and in endoatmosphere flying for long time, its movement locus due to often showing " jump " feature.
From the point of view of the domestic and international research conditions of trajectory prediction at present, the external research to ballistic prediction method is only limitted to trajectory In terms of guided missile, the not disclosed report of research to maneuvering target Forecasting Methodology.And the domestic research to ballistic prediction method is also led In terms of concentrating on ballistic missile, only result of study to maneuvering target ballistic prediction method is current using target maneuver Statistical model, simply calculated forward using Kalman's forecast principle, do not account for the complex situations such as the aerodynamic force of target.Always It, there is presently no the method for carrying out high-precision ballistic prediction to hypersonic gliding bullet.
The content of the invention
The present invention is to solve existing method it is low to maneuvering target ballistic prediction precision the problem of, and one kind is provided and closes on sky Between it is hypersonic gliding bullet Skipping trajectory Forecasting Methodology.
A kind of hypersonic gliding bullet Skipping trajectory Forecasting Methodology of near space, it is realized according to the following steps：
First, the foundation of the trajectory of hypersonic gliding bullet；
Establish the trajectory of the hypersonic gliding bullet of near space, including aerodynamics evaluation equation, acceleration calculation Equation, speed and position accounting equation, speed and trajectory angle accounting equation, and attitude angle accounting equation；
COMPREHENSIVE CALCULATING equation, obtain hypersonic gliding bullet centroid velocity and position equation such as formula (1)~(3) institute Show：
Wherein, xI、yIAnd zIPosition of the hypersonic gliding bullet in geocentric inertial coordinate system is represented,Represent it is hypersonic gliding bullet position to earth center distance,WithThen represent speed of the hypersonic gliding bullet in geocentric inertial coordinate system, aaxI、aayIAnd aazIRepresent aerodynamic force production Component of the raw hypersonic gliding bullet Maneuver Acceleration in geocentric inertial coordinate system, μ is gravitational parameter；
2nd, the tracking filter design of hypersonic gliding bullet；
The tracing of the movement Kalman filter of hypersonic gliding bullet is designed, tracks in real time and estimates high ultrasound Position, speed and the acceleration of ski-running Xiang bullet；
Wherein, the tracing of the movement Kalman filter of the hypersonic gliding bullet is
In above formula,The state estimation of wave filter i kth step is represented,Represent wave filter i from kth step to The state forecast value of the step of kth+1, Pi(k) the state estimation error co-variance matrix of wave filter i kth step, P are representedi(k+1/k) generation Table wave filter i walks the state forecast error co-variance matrix to the step of kth+1, Q from kthiFor wave filter i model prediction errors association side Poor matrix；
The measurement update equation of the tracing of the movement Kalman filter is
Wherein, subscriptTRepresent matrix transposition computing, subscript-1Represent matrix inversion operation, wave filter i calculation matrix
Hi=[1 0 0], i=x, y, z (7)
RiFor wave filter i measurement error covariance matrixes, I represents unit square, Ki(k+1) step of wave filter i kth+1 is represented Filtering gain matrix, ηi(k+1) metrical information of the step of wave filter i kth+1, P are representedi(k+1) shape of the step of wave filter i kth+1 is represented State evaluated error covariance matrix；
3rd, the ballistic prediction of hypersonic gliding bullet；
(3 1) terminate the hypersonic cunning at kT moment according to tracing of the movement Kalman filter tracking in step 2 The acceleration of Xiang bullet, wherein, k is the integer more than 0, and T is calculating cycle, typically desirable T=0.01s, in applying step one Near space it is hypersonic gliding bullet acceleration calculation equation obtain tracking terminate kT moment near spaces it is hypersonic Glide acceleration caused by the aerodynamic force of bullet, and terminates the speed of kT moment hypersonic gliding bullet, application according to tracking The trajectory angle accounting equation of the hypersonic gliding bullet of near space in step 1 calculates its trajectory angle；
(three or two) accelerated according to caused by the aerodynamic force of the hypersonic gliding bullet of the kT moment near spaces calculated Degree and trajectory angle, calculated with reference to the aerodynamics evaluation equation and attitude angle of the hypersonic gliding bullet of the near space in step 1 Equation estimation tracking terminate the kT moment it is hypersonic gliding bullet flight when the angle of attack and attitude angle；
(three or three) are located at subsequent kT to (k+1) T, and in (k+2) T ..., (k+N) T time, N is just whole much larger than 2 Number, the hypersonic gliding bullet of near space are carried out waiting the angle of attack and wait roll angle flight, then jump up and down is showed in fore-and-aft plane The trajectory of jump formula；Meanwhile hypersonic gliding bullet uses BTT control modes, yaw angle β is equal to zero, is rolled by controlling Corner in its system, laterally turned by one Maneuver Acceleration of output, i.e. roll angle instruction γcKeep constant value；
(three or four) calculate and added caused by kT moment aerodynamic force according to waiting the angle of attack, zero yaw angle and waiting roll angle flying condition Speed, the position x at kT momentI、yIAnd zI, the kT moment speed vxI、vyIAnd vzISubstitution formula (1)~(3), accumulated by a step Numerical Partite transport obtains hypersonic gliding bullet in speed of (k+1) the T moment in geocentric inertial coordinate system and position after calculating；
The aerodynamic force that hypersonic gliding bullet is subject at (k+1) T moment is calculated, calculates hypersonic gliding bullet in (k + 1) the aerodynamic force acceleration that the T moment is subject to, then aerodynamic force acceleration (k+1) T moment and the position at (k+1) T moment and Speed substitutes into (1)~(3), calculates position and the speed at (k+2) T moment, by that analogy, calculates until (k+N) T moment Position and speed.
Invention effect：
Under these conditions, the hypersonic gliding bullet final position for emulating to obtain by 100 Monte-Carlo is pre- The statistical result of report is as shown in figure 13, and it is 95% that probability of the terminal location prediction error less than 50km, which is calculated,.
For the hypersonic gliding bullet of near space, the present invention considers target by complex situations such as aerodynamic force, proposes A kind of new jump ballistic prediction method, ballistic prediction precision is improved compared to conventional method, reduces ballistic prediction mistake Difference, the blank near space hypersonic aircraft trajectory prediction field is filled up.
Brief description of the drawings
Fig. 1 is geocentric inertial coordinate system and sensing point inertial coodinate system figure in embodiment two；
Fig. 2 is hypersonic gliding bullet trajectory figure in emulation experiment, and upper figure is fore-and-aft plane trajectory, figure below For lateral plane trajectory；
Fig. 3 is hypersonic gliding bullet X-axis acceleration estimation figure in emulation experiment；
Fig. 4 is hypersonic gliding bullet Y-axis acceleration estimation figure in emulation experiment；
Fig. 5 is hypersonic gliding bullet Z axis acceleration estimation figure in emulation experiment；
Fig. 6 is hypersonic gliding bullet X-axis velocity estimation figure in emulation experiment；
Fig. 7 is hypersonic gliding bullet Y-axis velocity estimation figure in emulation experiment；
Fig. 8 is hypersonic gliding bullet Z axis velocity estimation figure in emulation experiment；
Fig. 9 be emulation experiment in it is hypersonic gliding bullet X-axis location estimation figure, upper figure be X-axis position actual value with Estimate, figure below are the actual value of X-axis position and the difference of estimate；
Figure 10 be emulation experiment in it is hypersonic gliding bullet Y-axis location estimation figure, upper figure be Y-axis position actual value with Estimate, figure below are the actual value of Y-axis position and the difference of estimate；
Figure 11 be emulation experiment in it is hypersonic gliding bullet Z axis location estimation figure, upper figure be Z axis position actual value with Estimate, figure below are the actual value of Z axis position and the difference of estimate；
Figure 12 is hypersonic gliding bullet position estimation error figure in emulation experiment；
Figure 13 is hypersonic gliding bullet terminal location prediction error statistical result in emulation experiment.
Embodiment
Embodiment one：A kind of hypersonic gliding bullet Skipping trajectory prediction side of near space of present embodiment Method, it is realized according to the following steps：
First, the foundation of the trajectory of hypersonic gliding bullet；
Establish the trajectory of the hypersonic gliding bullet of near space, including aerodynamics evaluation equation, acceleration calculation Equation, speed and position accounting equation, speed and trajectory angle accounting equation, and attitude angle accounting equation；
COMPREHENSIVE CALCULATING equation, obtain hypersonic gliding bullet centroid velocity and position equation such as formula (1)~(3).
2nd, the tracking filter design of hypersonic gliding bullet；
The tracing of the movement Kalman filter of hypersonic gliding bullet is designed, tracks in real time and estimates high ultrasound Position, speed and the acceleration of ski-running Xiang bullet；
Wherein, the state a step of forecasting equation of the tracing of the movement Kalman filter of the hypersonic gliding bullet For (4), the measurement update equation of the tracing of the movement Kalman filter is (6)；
3rd, the ballistic prediction of hypersonic gliding bullet；
(3 1) terminate the hypersonic cunning at kT moment according to tracing of the movement Kalman filter tracking in step 2 The acceleration of Xiang bullet, wherein, k is the integer more than 0, and T is calculating cycle, typically desirable T=0.01s, in applying step one Near space it is hypersonic gliding bullet acceleration calculation equation obtain tracking terminate kT moment near spaces it is hypersonic Glide acceleration caused by the aerodynamic force of bullet, and terminates the speed of kT moment hypersonic gliding bullet, application according to tracking The trajectory angle accounting equation of the hypersonic gliding bullet of near space in step 1 calculates its trajectory angle；
(three or two) accelerated according to caused by the aerodynamic force of the hypersonic gliding bullet of the kT moment near spaces calculated Degree and trajectory angle, calculated with reference to the aerodynamics evaluation equation and attitude angle of the hypersonic gliding bullet of the near space in step 1 Equation estimation tracking terminate the kT moment it is hypersonic gliding bullet flight when the angle of attack and attitude angle；
(three or three) are located at subsequent kT to (k+1) T, and in (k+2) T ..., (k+N) T time, N is just whole much larger than 2 Number, the hypersonic gliding bullet of near space are carried out waiting the angle of attack and wait roll angle flight, then jump up and down is showed in fore-and-aft plane The trajectory of jump formula；Meanwhile hypersonic gliding bullet uses BTT control modes, yaw angle β is equal to zero, is rolled by controlling Corner in its system, laterally turned by one Maneuver Acceleration of output, i.e. roll angle instruction γcKeep constant value；
(three or four) calculate and added caused by kT moment aerodynamic force according to waiting the angle of attack, zero yaw angle and waiting roll angle flying condition Speed, the position x at kT momentI、yIAnd zI, the kT moment speed vxI、vyIAnd vzISubstitution formula (1)~(3), accumulated by a step Numerical Partite transport obtains hypersonic gliding bullet in speed of (k+1) the T moment in geocentric inertial coordinate system and position after calculating；
The aerodynamic force that hypersonic gliding bullet is subject at (k+1) T moment is calculated, calculates hypersonic gliding bullet in (k + 1) the aerodynamic force acceleration that the T moment is subject to, then aerodynamic force acceleration (k+1) T moment and the position at (k+1) T moment and Speed substitutes into formula (1)~(3), calculates position and the speed at (k+2) T moment, by that analogy, calculates until (k+N) T moment Position and speed.
Embodiment two：Present embodiment is unlike embodiment one：
The definition of first step coordinate system
In order to draw the Forecasting Methodology of the hypersonic gliding bullet Skipping trajectory of near space, it is necessary to establish a prescription journey The motion of the hypersonic gliding bullet of near space is described, therefore, firstly the need of the several relevant coordinate systems of definition.
(1) geocentric inertial coordinate system (oIxIyIzI)
As shown in figure 1, origin oIPositioned at ground in the heart, oIyIPoint to ground radar sensing point, oIxIAxle position is attacked flat in target In face, with oIyIAxle is vertical, and it is just o to point to targetIzIDetermined by the right-hand rule.
(2) sensing point inertial coodinate system (oFxyz)
Origin is ground radar point oF, oFY-axis plummet is upward, oFX-axis is attacked in plane in target, perpendicular to oFy Axle, it is just o to point to targetFZ-axis is determined by the right-hand rule.This coordinate system is fixed on inertial space in MISSILE LAUNCHING moment.It with The relation of geocentric inertial coordinate system is as shown in figure 1, i.e. oFx、oFY and oFZ-axis is respectively parallel to oIxI、oIyIAnd oIzIAxle.
(3) target inertial coordinate system (otFxtytzt)
Origin is objective emission point otF, by MISSILE LAUNCHING point inertial coodinate system around oFY-axis rotates 180 °, that is, obtains this seat Mark system.This coordinate system is used for the trajectory angle for determining target.
(4) objective body coordinate system (ox1y1z1)
Origin o is located on target centroid.ox1Axle overlaps with its longitudinal axis, and it is just oy to point to head1Axle is longitudinally asymmetric flat at its In face, vertical ox1Axle, point up as just, oz1Direction of principal axis is determined by the right-hand rule.
(5) transformational relation between several coordinate systems
(6) transformation matrix by sensing point inertial coodinate system to target inertial coordinate system
(7) transformation matrix by target inertial coordinate system to objective body coordinate system
Use the coordinate conversion matrix that 231 turns of sequences are derived for
Wherein,, ψ and γ be respectively target the angle of pitch, yaw angle and roll angle.
Step 1 is specially：
(1) hypersonic gliding bullet aerodynamic force is calculated
In glide phase, projection of the aerodynamic force that hypersonic gliding bullet is subject in system calculates according to the following formula：
Wherein, α is the angle of attack,
cyFor lift coefficient, cxFor resistance coefficient, their velocity correlations with hypersonic gliding bullet, it is respectively necessary for looking into Tables 1 and 2 is tried to achieve.Ma in Tables 1 and 2 represents Mach number, i.e., the speed v divided by velocity of sound of hypersonic gliding bullet are superb Speed of the velocity of sound gliding bullet when near space glides is about 15Ma；S is characterized area of reference, if hypersonic gliding bullet The feature reference area of head is S=0.4837m2；
Q=ρ v2/2 (9)
Wherein, q is dynamic head, and ρ is atmospheric density, and it changes with flying height h, and the atmospheric density table that can look into standard obtains Arrive；Hypersonic gliding bullet flying height h is calculated with following formula:
Wherein, xI、yIAnd zIRepresent position of the hypersonic gliding bullet in geocentric inertial coordinate system, Re=6378km For earth mean radius；
The lift coefficient c of the hypersonic gliding bullet of table 1y
The resistance coefficient c of the hypersonic gliding bullet of table 2x
It is located at in-flight yaw angle β=0, therefore aerodynamic force side force Rz1=0；
(2) speed of hypersonic gliding bullet is calculated
It is hypersonic gliding bullet speed be
Wherein, vxI、vyIAnd vzIRepresent speed of the hypersonic gliding bullet in geocentric inertial coordinate system；
The velocity of hypersonic gliding bullet is projected into target inertial coordinate system vxt0、vyt0With vzt0, i.e.,
Wherein,
(3) the trajectory elevation angle theta and trajectory deflection angle ψ of hypersonic gliding bullet are calculatedv：
The attitude angle of hypersonic gliding bullet, ψ and γ calculate according to the following formula：
Wherein,For the angle of pitch, ψ is yaw angle, and γ is roll angle, γcInstructed for roll angle；
(4) projection of the hypersonic gliding bullet Maneuver Acceleration in geocentric inertial coordinate system is calculated
Wherein,
M is hypersonic gliding warhead mass, it is assumed that m=907.2kg；
(5) hypersonic gliding bullet centroid velocity and position equation such as formula (1)~(3) are shown, wherein containing The acceleration of gravity that hypersonic gliding bullet is subject to, i.e.,
Wherein, gxI、gyIWith gzIRepresent projection of the acceleration of gravity in Earth central inertial system.
When Maneuver Acceleration caused by the aerodynamic force of the hypersonic gliding bullet of target has fully showed, tracking is filtered Ripple device comes into stable state, after the Maneuver Acceleration for fully estimating target, is transferred to and is carried out using trajectory prediction algorithm Prediction.
Other steps and parameter are identical with embodiment one.
Embodiment three：Present embodiment is unlike embodiment one or two：Step 2 is specially：
(1) on three axles of sensing point inertial coodinate system, hypersonic gliding is described respectively with Singer models The change of the component of acceleration of bullet, i.e.,
Wherein, λx、λyAnd λzThe o under sensing point inertial coodinate system is represented respectivelyFX-axis, oFY-axis and oFTarget on z-axis direction The inverse of time kept in reserve constant, wx、wyAnd wzO under sensing point inertial coodinate system is represented respectivelyFX-axis, oFY-axis and oFOn z-axis direction Zero mean Gaussian white noise；Acceleration ax、ayAnd azIn both include acceleration a caused by aerodynamic forceaxI、aayIAnd aazI, wrap again Containing acceleration caused by gravity；The o of sensing point inertial coodinate systemFX-axis, oFY-axis and oFZ-axis and three axles of Earth central inertial system oIxI、oIyIAnd oIzIParallel, projection a of the acceleration caused by aerodynamic force in Earth central inertial system respectivelyaxI、aayIAnd aazIIt is exactly Projection of this acceleration in sensing point inertial system；
(2) following three groups of equations of motion are set up under sensing point inertial coodinate system
Wherein, x, y and z represent position of the hypersonic gliding bullet in sensing point inertial coodinate system, WithThen represent projection of the hypersonic gliding bullet speed in sensing point inertial coodinate system.Because sensing point inertia is sat Mark the o of systemFX-axis, oFY-axis and oFZ-axis and three axle o of Earth central inertial systemIxI、oIyIAnd oIzIIt is parallel respectively, so
And because sensing point inertial coodinate system only with respect to Earth central inertial system sensing point inertial coordinate ties up to oIyIDirection of principal axis On translate, so having
For system (18), x-axis subsystem state variable X is definedx=[x vx ax]T, then following state equation is obtained
Wherein,
After carrying out discretization to formula (23) with certain period Δ t, the state equation obtained from kth step to the step of kth+1 is
Xx(k+1)=ΦxXx(k)+ωx(k) (25)
Wherein,
ωx(k) x-axis subsystem state noise vector is represented, is zero mean Gaussian white noise vector.
The position of hypersonic gliding bullet is obtained by Ground Target Tracking measurement device, then measures equation and be
ηx=x+ υx (27)
Wherein, υxX-axis subsystem measurement noise vector is represented, is zero mean Gaussian white noise vector.
Similarly, for system (19), y-axis subsystem state variable X is definedy=[y vy ay]T, then following state side is obtained Journey
Wherein,
After carrying out discretization to formula (28) with certain period Δ t, the state equation obtained from kth step to the step of kth+1 is
Xy(k+1)=ΦyXy(k)+ωy(k) (30)
Wherein,
ωy(k) y-axis subsystem state noise vector is represented, is zero mean Gaussian white noise vector.
Measuring equation is
ηy=y+ υy (32)
Wherein, υyY-axis subsystem measurement noise vector is represented, is zero mean Gaussian white noise vector.
For system (20), z-axis subsystem state variable X is definedz=[z vz az]T, then following state equation is obtained
Wherein,
After carrying out discretization to formula (28) with certain period Δ t, the state equation obtained from kth step to the step of kth+1 is
Xz(k+1)=ΦzXz(k)+ωz(k) (35)
Wherein,
ωz(k) z-axis subsystem state noise vector is represented, is zero mean Gaussian white noise vector.
Measuring equation is
ηz=z+ υz (37)
Wherein, υzZ-axis subsystem measurement noise vector is represented, is zero mean Gaussian white noise vector.
The subsystem on three axles of sensing point inertial coodinate system, i.e., the x-axis subsystem that formula (23) and (27) are formed, formula (28) and (32) form y-axis subsystem, and formula (33) and (37) composition z-axis subsystem, separately design Kalman filter Device, its prognostic equation are formula (4), i.e.,
In above formula,The state estimation of wave filter i kth step is represented,Represent wave filter i from kth step to The state forecast value of the step of kth+1, Pi(k) the state estimation error co-variance matrix of wave filter i kth step, P are representedi(k+1/k) generation Table wave filter i walks the state forecast error co-variance matrix to the step of kth+1, Q from kthiFor wave filter i model prediction errors association side Poor matrix,
The measurement update equation of wave filter is formula (6), i.e.,
Wherein, subscriptTRepresent matrix transposition computing, subscript-1Matrix inversion operation is represented, wave filter i calculation matrix is Formula (7), i.e.,
Hi=[1 0 0], i=x, y, z
RiFor wave filter i measurement error covariance matrixes, I represents unit square, Ki(k+1) step of wave filter i kth+1 is represented Filtering gain matrix, Pi(k+1) the state estimation error co-variance matrix of the step of wave filter i kth+1, η are representedi(k+1) filtering is represented The metrical information of the step of device i kth+1, i.e. ηx(k+1)=x (k+1)+υx(k+1), ηy(k+1)=y (k+1)+υy(k+1),ηz(k+1) =z (k+1)+υz(k+1)。
When no Ground Target Tracking measurement device information, prognostic equation (4) is only run, there is Ground Target Tracking device During metrical information, then prognostic equation (4) and measurement update equation (6) are run simultaneously.
Other steps and parameter are identical with embodiment one or two.
Embodiment four：Unlike one of present embodiment and embodiment one to three：When hypersonic Maneuver Acceleration caused by gliding bullet aerodynamic force fully shows, and tracking filter comes into stable state, fully After the Maneuver Acceleration for estimating target, trajectory prediction algorithm is transferred to.
As shown in Figure 1, step (3 1) is specific as follows：
First, obtained tracking the hypersonic gliding bullet of finish time kT time space according to the estimated result of three subfilters Position [x y z] of the head in sensing point inertial systemT, the position in sensing point inertial coodinate system and the position in Earth central inertial system Relation is
[xI yI zI]T=[x y z]T+[0 Re 0]T (38)
Wherein, Re=6378km is earth mean radius；
According to the acceleration estimation value [a of the hypersonic gliding bullet of tracking finish timex ay az]T, due to sensing point Three of inertial coodinate system are parallel with three axles of Earth central inertial system, can obtain
[axI ayI azI]T=[ax ay az]T (39)
Wherein, axI、ayIAnd azIRepresent projection of the acceleration of hypersonic gliding bullet in Earth central inertial system
The aerodynamic force that currently hypersonic gliding bullet is subject to is calculated with following formula：
When 15Mach or so, lift-drag ratio is about 3:1, due to Rz1=0, it can obtain
According to the velocity estimation value [v of the hypersonic gliding bullet of tracking finish timex vy vz]T, because sensing point is used to Three of property coordinate system are parallel with three axles of Earth central inertial system, can obtain
[vxI vyI vzI]T=[vx vy vz]T (42)
According to the speed [v of hypersonic gliding bullet at the end of trackingxI vyI vzI]T, it is calculated with formula (12) and (13) Trajectory angle θ and ψv。
Other steps and parameter are identical with one of embodiment one to three.
Embodiment five：Unlike one of present embodiment and embodiment one to four：Step (three or two) It is specific as follows：
According to the speed [v of hypersonic gliding bullet at the end of trackingxI vyI vzI]T, its speed is calculated with formula (11) V, formula (9) is recycled to try to achieve current dynamic head q；
Calculating current lift coefficient with following formula is
When the angle of attack is 10 °, lift coefficient 0.37, so estimating the current angle of attack with following formula：
Hypersonic gliding bullet uses banked turn control mode, i.e. yaw angle β=0, then by before in formula (14) Two formulas can calculate the angle of pitchAnd yaw angle ψ；
Estimate current roll angle γ：Projection of the aerodynamic force in inertial system has been calculated by formula (40), it has been projected To objective body coordinate system, i.e.,
Above formula can further be write
Order
It can then be obtained according to formula (45)
Wherein, Ry1Estimated by formula (41), Rz1=0, sin γ and cos γ are tried to achieve with formula (46), then with following formula only One ground determines current roll angle
γ=tan2-1(sinγ,cosγ) (47)
Other steps and parameter are identical with one of embodiment one to four.
Emulation experiment：
For the hypersonic gliding bullet of near space, the present invention considers target by complex situations such as aerodynamic force, proposes A kind of new jump ballistic prediction method, ballistic prediction precision is improved compared to conventional method, reduces ballistic prediction mistake Difference, the blank near space hypersonic aircraft trajectory prediction field is filled up.The present invention is verified below by numerical simulation Effect.
If emulate initial time t0=82s, initial position of the hypersonic gliding bullet in sensing point inertial coodinate system For xI(0)=1161km, yI(0)=- 76km, zI(0)=- 1.1km, component of the initial velocity in sensing point inertial coodinate system For vxI(0)=- 3024m/s, vyI(0)=558m/s, vzI(0)=0m/s.It is hypersonic gliding bullet keep etc. the angle of attack and wait roll Corner flies, and the angle of attack is α=10 °, and roll angle is γ=45 °.Emulation finish time is tf=400s.
The initial estimate of hypersonic gliding bullet tracking Kalman filter is taken as
Xx(0)=[xI(0)+300 vxI(0)+20 0]T, Xy(0)=[yI(0)+300 vyI(0)+20 0]T,
Xz(0)=[zI(0)+300 vzI(0)+20 0]T
The initial value of state estimation variance matrix is taken as
Ground range error is taken as 300m, and angle error is taken as pitch orientation 0.5mrad, yaw direction 0.3mrad.Filtering The device a step of forecasting cycle is taken as 0.01s, and the filter measurement amendment cycle is that the data updating rate of ground survey information is 0.2s. Preceding 118s is modified using metrical information.200s carries out trajectory prediction afterwards.
By taking a simulation result as an example, the trajectory of hypersonic gliding bullet is as shown in Fig. 2 acceleration estimation result As shown in Fig. 3-Fig. 5, velocity estimation result is as shown in Fig. 6-Fig. 8, and location estimation is as shown in Fig. 9-Figure 11, total site error As shown in figure 12.
Under these conditions, the hypersonic gliding bullet final position for emulating to obtain by 100 Monte-Carlo is pre- The statistical result of report is as shown in figure 13, and it is 95% that probability of the terminal location prediction error less than 50km, which is calculated,.

Claims (1)

1. a kind of hypersonic gliding bullet Skipping trajectory Forecasting Methodology of near space, it is characterised in that it is real according to the following steps It is existing：
First, the foundation of the trajectory of hypersonic gliding bullet；
Establish near space it is hypersonic gliding bullet trajectory, including aerodynamics evaluation equation, acceleration calculation equation, Speed and position accounting equation, speed and trajectory angle accounting equation, and attitude angle accounting equation；
COMPREHENSIVE CALCULATING equation, obtain shown in hypersonic gliding bullet centroid velocity and position equation such as formula (1)~(3)：
<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mover> <mi>x</mi> <mo>&amp;CenterDot;</mo> </mover> <mi>I</mi> </msub> <mo>=</mo> <msub> <mi>v</mi> <mrow> <mi>x</mi> <mi>I</mi> </mrow> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mover> <mi>v</mi> <mo>&amp;CenterDot;</mo> </mover> <mrow> <mi>x</mi> <mi>I</mi> </mrow> </msub> <mo>=</mo> <mo>-</mo> <mfrac> <mrow> <msub> <mi>&amp;mu;x</mi> <mi>I</mi> </msub> </mrow> <msubsup> <mi>r</mi> <mi>I</mi> <mn>3</mn> </msubsup> </mfrac> <mo>+</mo> <msub> <mi>a</mi> <mrow> <mi>a</mi> <mi>x</mi> <mi>I</mi> </mrow> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>
<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mover> <mi>y</mi> <mo>&amp;CenterDot;</mo> </mover> <mi>I</mi> </msub> <mo>=</mo> <msub> <mi>v</mi> <mrow> <mi>y</mi> <mi>I</mi> </mrow> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mover> <mi>v</mi> <mo>&amp;CenterDot;</mo> </mover> <mrow> <mi>y</mi> <mi>I</mi> </mrow> </msub> <mo>=</mo> <mo>-</mo> <mfrac> <mrow> <msub> <mi>&amp;mu;y</mi> <mi>I</mi> </msub> </mrow> <msubsup> <mi>r</mi> <mi>I</mi> <mn>3</mn> </msubsup> </mfrac> <mo>+</mo> <msub> <mi>a</mi> <mrow> <mi>a</mi> <mi>y</mi> <mi>I</mi> </mrow> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow>
<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mover> <mi>z</mi> <mo>&amp;CenterDot;</mo> </mover> <mi>I</mi> </msub> <mo>=</mo> <msub> <mover> <mi>v</mi> <mo>&amp;CenterDot;</mo> </mover> <mrow> <mi>z</mi> <mi>I</mi> </mrow> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mover> <mi>v</mi> <mo>&amp;CenterDot;</mo> </mover> <mrow> <mi>z</mi> <mi>I</mi> </mrow> </msub> <mo>=</mo> <mo>-</mo> <mfrac> <mrow> <msub> <mi>&amp;mu;z</mi> <mi>I</mi> </msub> </mrow> <msubsup> <mi>r</mi> <mi>I</mi> <mn>3</mn> </msubsup> </mfrac> <mo>+</mo> <msub> <mi>a</mi> <mrow> <mi>a</mi> <mi>z</mi> <mi>I</mi> </mrow> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow>
Wherein, xI、yIAnd zIPosition of the hypersonic gliding bullet in geocentric inertial coordinate system is represented, Represent it is hypersonic gliding bullet position to earth center distance,WithThen represent high ultrasound Speed of the ski-running Xiang bullet in geocentric inertial coordinate system, aaxI、aayIAnd aazIRepresent hypersonic gliding bullet caused by aerodynamic force Component of the head Maneuver Acceleration in geocentric inertial coordinate system, μ is gravitational parameter；
2nd, the tracking filter design of hypersonic gliding bullet；
The tracing of the movement Kalman filter of hypersonic gliding bullet is designed, tracks in real time and estimates hypersonic cunning Position, speed and the acceleration of Xiang bullet；
Wherein, the state a step of forecasting equation of the tracing of the movement Kalman filter of the hypersonic gliding bullet is
<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mover> <mi>X</mi> <mo>&amp;OverBar;</mo> </mover> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>&amp;Phi;</mi> <mi>i</mi> </msub> <msub> <mover> <mi>X</mi> <mo>^</mo> </mover> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>P</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>/</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>&amp;Phi;</mi> <mi>i</mi> </msub> <msub> <mi>P</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <msubsup> <mi>&amp;Phi;</mi> <mi>i</mi> <mi>T</mi> </msubsup> <mo>+</mo> <msub> <mi>Q</mi> <mi>i</mi> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>,</mo> <mi>i</mi> <mo>=</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow>
In above formula,The state estimation of wave filter i kth step is represented,Wave filter i is represented from kth step to kth+1 The state forecast value of step, Pi(k) the state estimation error co-variance matrix of wave filter i kth step, P are representedi(k+1/k) filtering is represented Device i walks the state forecast error co-variance matrix to the step of kth+1, Q from kthiFor wave filter i model prediction error covariance squares Battle array,
<mrow> <msub> <mi>&amp;Phi;</mi> <mi>i</mi> </msub> <mo>=</mo> <mi>exp</mi> <mrow> <mo>(</mo> <msub> <mi>F</mi> <mi>i</mi> </msub> <mi>&amp;Delta;</mi> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mn>1</mn> </mtd> <mtd> <mrow> <mi>&amp;Delta;</mi> <mi>t</mi> </mrow> </mtd> <mtd> <mrow> <mn>1</mn> <mo>/</mo> <msubsup> <mi>&amp;lambda;</mi> <mi>i</mi> <mn>2</mn> </msubsup> <mrow> <mo>(</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msub> <mi>&amp;lambda;</mi> <mi>i</mi> </msub> <mi>t</mi> </mrow> </msup> <mo>+</mo> <msub> <mi>&amp;lambda;</mi> <mi>i</mi> </msub> <mi>&amp;Delta;</mi> <mi>t</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>1</mn> </mtd> <mtd> <mrow> <mn>1</mn> <mo>/</mo> <msub> <mi>&amp;lambda;</mi> <mi>i</mi> </msub> <mn>1</mn> <mo>-</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msub> <mi>&amp;lambda;</mi> <mi>i</mi> </msub> <mi>&amp;Delta;</mi> <mi>t</mi> </mrow> </msup> <mo>)</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msub> <mi>&amp;lambda;</mi> <mi>i</mi> </msub> <mi>&amp;Delta;</mi> <mi>t</mi> </mrow> </msup> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </mrow>
The measurement update equation of the tracing of the movement Kalman filter is
<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>K</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>P</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>/</mo> <mi>k</mi> <mo>)</mo> </mrow> <msubsup> <mi>H</mi> <mi>i</mi> <mi>T</mi> </msubsup> <msup> <mrow> <mo>&amp;lsqb;</mo> <msub> <mi>H</mi> <mi>i</mi> </msub> <msub> <mi>P</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>/</mo> <mi>k</mi> <mo>)</mo> </mrow> <msubsup> <mi>H</mi> <mi>i</mi> <mi>T</mi> </msubsup> <mo>+</mo> <msub> <mi>R</mi> <mi>i</mi> </msub> <mo>&amp;rsqb;</mo> </mrow> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mover> <mi>X</mi> <mo>^</mo> </mover> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mover> <mi>X</mi> <mo>&amp;OverBar;</mo> </mover> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>K</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>&amp;lsqb;</mo> <msub> <mi>&amp;eta;</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>H</mi> <mi>i</mi> </msub> <msub> <mover> <mi>X</mi> <mo>&amp;OverBar;</mo> </mover> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>P</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>=</mo> <mo>&amp;lsqb;</mo> <mi>I</mi> <mo>-</mo> <msub> <mi>K</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <msub> <mi>H</mi> <mi>i</mi> </msub> <mo>&amp;rsqb;</mo> <msub> <mi>P</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>/</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>,</mo> <mi>i</mi> <mo>=</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow>
Wherein, subscriptTRepresent matrix transposition computing, subscript-1Represent matrix inversion operation, wave filter i calculation matrix
Hi=[1 0 0], i=x, y, z (7)
RiFor wave filter i measurement error covariance matrixes, I represents unit matrix, Ki(k+1) filter of the step of wave filter i kth+1 is represented Ripple gain matrix, ηi(k+1) metrical information of the step of wave filter i kth+1, P are representedi(k+1) state of the step of wave filter i kth+1 is represented Evaluated error covariance matrix；
3rd, the ballistic prediction of hypersonic gliding bullet；
(3 1) terminate the hypersonic gliding bullet at kT moment according to tracing of the movement Kalman filter tracking in step 2 The acceleration of head, wherein, k is the integer more than 0, and T is calculating cycle, typically desirable T=0.01s, facing in applying step one The acceleration calculation equation of the hypersonic gliding bullet of near space obtains tracking and terminates the hypersonic gliding of kT moment near spaces Acceleration caused by the aerodynamic force of bullet, and according to the speed of tracking end kT moment hypersonic gliding bullet, applying step The trajectory angle accounting equation of the hypersonic gliding bullet of near space in one calculates its trajectory angle；
(three or two) according to caused by the aerodynamic force of the hypersonic gliding bullet of the kT moment near spaces that calculate acceleration and Trajectory angle, with reference to the aerodynamics evaluation equation and attitude angle accounting equation of the hypersonic gliding bullet of the near space in step 1 Estimation tracking terminate the kT moment it is hypersonic gliding bullet flight when the angle of attack and attitude angle；
(three or three) are located at subsequent kT to (k+1) T, and in (k+2) T ..., (k+N) T time, N is the positive integer more than 2, is closed on The hypersonic gliding bullet in space is carried out waiting the angle of attack and waits roll angle flight, then flying for vertical bounce formula is showed in fore-and-aft plane Row trajectory；Meanwhile hypersonic gliding bullet uses BTT control modes, yaw angle β is equal to zero, by controlling roll angle at it System laterally turned by one Maneuver Acceleration of output, i.e. roll angle instruction γcKeep constant value；
(three or four) according to waiting the angle of attack, zero yaw angle and waiting roll angle flying condition, calculate acceleration caused by kT moment aerodynamic force, The position x at kT momentI、yIAnd zI, the kT moment speed vxI、vyIAnd vzISubstitution formula (1)~(3), integrate and transport by a step Numerical Hypersonic gliding bullet is obtained after calculation in speed of (k+1) the T moment in geocentric inertial coordinate system and position；
The aerodynamic force that hypersonic gliding bullet is subject at (k+1) T moment is calculated, calculates hypersonic gliding bullet in (k+1) T The aerodynamic force acceleration that moment is subject to, then aerodynamic force acceleration and the position at (k+1) T moment and speed (k+1) T moment Substitution formula (1)~(3), position and the speed at (k+2) T moment are calculated, by that analogy, is calculated until the position at (k+N) T moment Put and speed；
Step 1 is specially：
(1) hypersonic gliding bullet aerodynamic force is calculated
In glide phase, projection of the aerodynamic force that hypersonic gliding bullet is subject in system calculates according to the following formula：
<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>R</mi> <mrow> <mi>x</mi> <mn>1</mn> </mrow> </msub> <mo>=</mo> <mo>-</mo> <msub> <mi>c</mi> <mi>x</mi> </msub> <mi>q</mi> <mi>S</mi> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>R</mi> <mrow> <mi>y</mi> <mn>1</mn> </mrow> </msub> <mo>=</mo> <msub> <mi>c</mi> <mi>y</mi> </msub> <mi>s</mi> <mi>i</mi> <mi>g</mi> <mi>n</mi> <mrow> <mo>(</mo> <mi>&amp;alpha;</mi> <mo>)</mo> </mrow> <mi>q</mi> <mi>S</mi> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>R</mi> <mrow> <mi>z</mi> <mn>1</mn> </mrow> </msub> <mo>=</mo> <mn>0</mn> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>8</mn> <mo>)</mo> </mrow> </mrow>
Wherein, α is the angle of attack,
<mrow> <mi>s</mi> <mi>i</mi> <mi>g</mi> <mi>n</mi> <mrow> <mo>(</mo> <mi>&amp;alpha;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mn>1</mn> <mo>,</mo> <mi>&amp;alpha;</mi> <mo>&gt;</mo> <mn>0</mn> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>-</mo> <mn>1</mn> <mo>,</mo> <mi>&amp;alpha;</mi> <mo>&lt;</mo> <mn>0</mn> </mrow> </mtd> </mtr> </mtable> </mfenced> </mrow>
cyFor lift coefficient, cxFor resistance coefficient, S is characterized area of reference；
Q=ρ v2/2 (9)
Wherein, q is dynamic head, and v is the speed of hypersonic gliding bullet, and ρ is atmospheric density, and it changes with flying height h, The atmospheric density table that standard can be looked into obtains；Hypersonic gliding bullet flying height h is calculated with following formula:
<mrow> <mi>h</mi> <mo>=</mo> <msqrt> <mrow> <msubsup> <mi>x</mi> <mi>I</mi> <mn>2</mn> </msubsup> <mo>+</mo> <msubsup> <mi>y</mi> <mi>I</mi> <mn>2</mn> </msubsup> <mo>+</mo> <msubsup> <mi>z</mi> <mi>I</mi> <mn>2</mn> </msubsup> </mrow> </msqrt> <mo>-</mo> <msub> <mi>R</mi> <mi>e</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>10</mn> <mo>)</mo> </mrow> </mrow>
Wherein, xI、yIAnd zIRepresent position of the hypersonic gliding bullet in geocentric inertial coordinate system, Re=6378km is the earth Mean radius；
It is located at in-flight yaw angle β=0, therefore aerodynamic force side force Rz1=0；
(2) speed of hypersonic gliding bullet is calculated
It is hypersonic gliding bullet speed be
<mrow> <mi>v</mi> <mo>=</mo> <msqrt> <mrow> <msubsup> <mi>v</mi> <mrow> <mi>x</mi> <mi>I</mi> </mrow> <mn>2</mn> </msubsup> <mo>+</mo> <msubsup> <mi>v</mi> <mrow> <mi>y</mi> <mi>I</mi> </mrow> <mn>2</mn> </msubsup> <mo>+</mo> <msubsup> <mi>v</mi> <mrow> <mi>z</mi> <mi>I</mi> </mrow> <mn>2</mn> </msubsup> </mrow> </msqrt> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>11</mn> <mo>)</mo> </mrow> </mrow>
Wherein, vxI、vyIAnd vzIRepresent speed of the hypersonic gliding bullet in geocentric inertial coordinate system；
The velocity of hypersonic gliding bullet is projected into target inertial coordinate system vxt0、vyt0With vzt0, i.e.,
<mrow> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mi>v</mi> <mrow> <mi>x</mi> <mi>t</mi> <mn>0</mn> </mrow> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>v</mi> <mrow> <mi>y</mi> <mi>t</mi> <mn>0</mn> </mrow> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>v</mi> <mrow> <mi>z</mi> <mi>t</mi> <mn>0</mn> </mrow> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>=</mo> <msub> <mi>C</mi> <mrow> <mn>0</mn> <mo>&amp;RightArrow;</mo> <mi>t</mi> <mn>0</mn> </mrow> </msub> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mi>v</mi> <mrow> <mi>x</mi> <mi>I</mi> </mrow> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>v</mi> <mrow> <mi>y</mi> <mi>I</mi> </mrow> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>v</mi> <mrow> <mi>z</mi> <mi>I</mi> </mrow> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>12</mn> <mo>)</mo> </mrow> </mrow>
Wherein,
<mrow> <msub> <mi>C</mi> <mrow> <mn>0</mn> <mo>&amp;RightArrow;</mo> <mi>t</mi> <mn>0</mn> </mrow> </msub> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>1</mn> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </mtd> </mtr> </mtable> </mfenced> </mrow>
(3) the trajectory elevation angle theta and trajectory deflection angle ψ of hypersonic gliding bullet are calculatedv：
<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mi>&amp;theta;</mi> <mo>=</mo> <msup> <mi>sin</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mrow> <mo>(</mo> <msub> <mi>v</mi> <mrow> <mi>y</mi> <mi>t</mi> <mn>0</mn> </mrow> </msub> <mo>/</mo> <mi>v</mi> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>&amp;psi;</mi> <mi>v</mi> </msub> <mo>=</mo> <mo>-</mo> <msup> <mi>tan</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mrow> <mo>(</mo> <msub> <mi>v</mi> <mrow> <mi>z</mi> <mi>t</mi> <mn>0</mn> </mrow> </msub> <mo>/</mo> <msub> <mi>v</mi> <mrow> <mi>x</mi> <mi>t</mi> <mn>0</mn> </mrow> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>13</mn> <mo>)</mo> </mrow> </mrow>
The attitude angle of hypersonic gliding bulletψ and γ are calculated according to the following formula：
Wherein,For the angle of pitch, ψ is yaw angle, and γ is roll angle, γcInstructed for roll angle；
(4) projection of the hypersonic gliding bullet Maneuver Acceleration in geocentric inertial coordinate system is calculated
<mrow> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mi>a</mi> <mrow> <mi>a</mi> <mi>x</mi> <mi>I</mi> </mrow> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>a</mi> <mrow> <mi>a</mi> <mi>y</mi> <mi>I</mi> </mrow> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>a</mi> <mrow> <mi>a</mi> <mi>z</mi> <mi>I</mi> </mrow> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>=</mo> <msubsup> <mi>C</mi> <mrow> <mn>0</mn> <mo>&amp;RightArrow;</mo> <mi>t</mi> <mn>0</mn> </mrow> <mi>T</mi> </msubsup> <msubsup> <mi>C</mi> <mrow> <mi>t</mi> <mn>0</mn> <mo>&amp;RightArrow;</mo> <mn>1</mn> </mrow> <mi>T</mi> </msubsup> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <msub> <mi>R</mi> <mrow> <mi>x</mi> <mn>1</mn> </mrow> </msub> <mo>/</mo> <mi>m</mi> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>R</mi> <mrow> <mi>y</mi> <mn>1</mn> </mrow> </msub> <mo>/</mo> <mi>m</mi> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>R</mi> <mrow> <mi>z</mi> <mn>1</mn> </mrow> </msub> <mo>/</mo> <mi>m</mi> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>15</mn> <mo>)</mo> </mrow> </mrow>
Wherein,
M is hypersonic gliding warhead mass；
(5) speed and the position of hypersonic gliding bullet are calculated
Shown in hypersonic gliding bullet centroid velocity and position equation such as formula (1)~(3), wherein containing hypersonic The acceleration of gravity that gliding bullet is subject to, i.e.,
<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>g</mi> <mrow> <mi>x</mi> <mi>I</mi> </mrow> </msub> <mo>=</mo> <mo>-</mo> <mfrac> <mrow> <msub> <mi>&amp;mu;x</mi> <mi>I</mi> </msub> </mrow> <msubsup> <mi>r</mi> <mi>I</mi> <mn>3</mn> </msubsup> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>g</mi> <mrow> <mi>y</mi> <mi>I</mi> </mrow> </msub> <mo>=</mo> <mo>-</mo> <mfrac> <mrow> <msub> <mi>&amp;mu;y</mi> <mi>I</mi> </msub> </mrow> <msubsup> <mi>r</mi> <mi>I</mi> <mn>3</mn> </msubsup> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>g</mi> <mrow> <mi>z</mi> <mi>I</mi> </mrow> </msub> <mo>=</mo> <mo>-</mo> <mfrac> <mrow> <msub> <mi>&amp;mu;z</mi> <mi>I</mi> </msub> </mrow> <msubsup> <mi>r</mi> <mi>I</mi> <mn>3</mn> </msubsup> </mfrac> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>16</mn> <mo>)</mo> </mrow> </mrow>
Wherein, gxI、gyIWith gzIRepresent projection of the acceleration of gravity in Earth central inertial system；
Step 2 is specially：
(1) on three axles of sensing point inertial coodinate system, hypersonic gliding bullet is described respectively with Singer models Component of acceleration change, i.e.,
<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mover> <mi>a</mi> <mo>&amp;CenterDot;</mo> </mover> <mi>x</mi> </msub> <mo>=</mo> <mo>-</mo> <msub> <mi>&amp;lambda;</mi> <mi>x</mi> </msub> <msub> <mi>a</mi> <mi>x</mi> </msub> <mo>+</mo> <msub> <mi>w</mi> <mi>x</mi> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mover> <mi>a</mi> <mo>&amp;CenterDot;</mo> </mover> <mi>y</mi> </msub> <mo>=</mo> <mo>-</mo> <msub> <mi>&amp;lambda;</mi> <mi>y</mi> </msub> <msub> <mi>a</mi> <mi>y</mi> </msub> <mo>+</mo> <msub> <mi>w</mi> <mi>y</mi> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mover> <mi>a</mi> <mo>&amp;CenterDot;</mo> </mover> <mi>z</mi> </msub> <mo>=</mo> <mo>-</mo> <msub> <mi>&amp;lambda;</mi> <mi>z</mi> </msub> <msub> <mi>a</mi> <mi>z</mi> </msub> <mo>+</mo> <msub> <mi>w</mi> <mi>z</mi> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>17</mn> <mo>)</mo> </mrow> </mrow>
Wherein, λx、λyAnd λzThe o under sensing point inertial coodinate system is represented respectivelyFX-axis, oFY-axis and oFTarget maneuver on z-axis direction The inverse of time constant, wx、wyAnd wzO under sensing point inertial coodinate system is represented respectivelyFX-axis, oFY-axis and oFZero on z-axis direction Average white Gaussian noise；Acceleration ax、ayAnd azIn both include acceleration a caused by aerodynamic forceaxI、aayIAnd aazI, and include weight Acceleration g caused by powerxI、gyIAnd gzI；The o of sensing point inertial coodinate systemFX-axis, oFY-axis and oFZ-axis and the three of Earth central inertial system Individual axle oIxI、oIyIAnd oIzIParallel, projection a of the acceleration caused by aerodynamic force in Earth central inertial system respectivelyaxI、aayIAnd aazI It is exactly projection of this acceleration in sensing point inertial system；
(2) following three groups of equations of motion are set up under sensing point inertial coodinate system
<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mover> <mi>x</mi> <mo>&amp;CenterDot;</mo> </mover> <mo>=</mo> <msub> <mi>v</mi> <mi>x</mi> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mover> <mi>v</mi> <mo>&amp;CenterDot;</mo> </mover> <mi>x</mi> </msub> <mo>=</mo> <msub> <mi>a</mi> <mi>x</mi> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mover> <mi>a</mi> <mo>&amp;CenterDot;</mo> </mover> <mi>x</mi> </msub> <mo>=</mo> <mo>-</mo> <msub> <mi>&amp;lambda;</mi> <mi>x</mi> </msub> <msub> <mi>a</mi> <mi>x</mi> </msub> <mo>+</mo> <msub> <mi>w</mi> <mi>x</mi> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>18</mn> <mo>)</mo> </mrow> </mrow>
<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mover> <mi>y</mi> <mo>&amp;CenterDot;</mo> </mover> <mo>=</mo> <msub> <mi>v</mi> <mi>y</mi> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mover> <mi>v</mi> <mo>&amp;CenterDot;</mo> </mover> <mi>y</mi> </msub> <mo>=</mo> <msub> <mi>a</mi> <mi>y</mi> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mover> <mi>a</mi> <mo>&amp;CenterDot;</mo> </mover> <mi>y</mi> </msub> <mo>=</mo> <mo>-</mo> <msub> <mi>&amp;lambda;</mi> <mi>y</mi> </msub> <msub> <mi>a</mi> <mi>y</mi> </msub> <mo>+</mo> <msub> <mi>w</mi> <mi>y</mi> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>19</mn> <mo>)</mo> </mrow> </mrow>
<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mover> <mi>z</mi> <mo>&amp;CenterDot;</mo> </mover> <mo>=</mo> <msub> <mi>v</mi> <mi>z</mi> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mover> <mi>v</mi> <mo>&amp;CenterDot;</mo> </mover> <mi>z</mi> </msub> <mo>=</mo> <msub> <mi>a</mi> <mi>z</mi> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mover> <mi>a</mi> <mo>&amp;CenterDot;</mo> </mover> <mi>z</mi> </msub> <mo>=</mo> <mo>-</mo> <msub> <mi>&amp;lambda;</mi> <mi>z</mi> </msub> <msub> <mi>a</mi> <mi>z</mi> </msub> <mo>+</mo> <msub> <mi>w</mi> <mi>z</mi> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>20</mn> <mo>)</mo> </mrow> </mrow>
Wherein, x, y and z represent position of the hypersonic gliding bullet in sensing point inertial coodinate system,WithThen represent projection of the hypersonic gliding bullet speed in sensing point inertial coodinate system；Because sensing point inertial coordinate The o of systemFX-axis, oFY-axis and oFZ-axis and three axle o of Earth central inertial systemIxI、oIyIAnd oIzIIt is parallel respectively, so
<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>v</mi> <mi>x</mi> </msub> <mo>=</mo> <msub> <mi>v</mi> <mrow> <mi>x</mi> <mi>I</mi> </mrow> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>v</mi> <mi>y</mi> </msub> <mo>=</mo> <msub> <mi>v</mi> <mrow> <mi>y</mi> <mi>I</mi> </mrow> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>v</mi> <mi>z</mi> </msub> <mo>=</mo> <msub> <mi>v</mi> <mrow> <mi>z</mi> <mi>I</mi> </mrow> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>21</mn> <mo>)</mo> </mrow> </mrow>
And because sensing point inertial coodinate system only with respect to Earth central inertial system sensing point inertial coordinate ties up to oIyIMake on direction of principal axis Translation, so have
<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mi>x</mi> <mo>=</mo> <msub> <mi>x</mi> <mi>I</mi> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>y</mi> <mo>=</mo> <msub> <mi>y</mi> <mi>I</mi> </msub> <mo>-</mo> <msub> <mi>R</mi> <mi>e</mi> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>z</mi> <mo>=</mo> <msub> <mi>z</mi> <mi>I</mi> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>22</mn> <mo>)</mo> </mrow> </mrow>
For system (18), x-axis subsystem state variable X is definedx=[x vx ax]T, then following state equation is obtained
<mrow> <msub> <mover> <mi>X</mi> <mo>&amp;CenterDot;</mo> </mover> <mi>x</mi> </msub> <mo>=</mo> <msub> <mi>F</mi> <mi>x</mi> </msub> <msub> <mi>X</mi> <mi>x</mi> </msub> <mo>+</mo> <msub> <mi>G</mi> <mi>x</mi> </msub> <msub> <mi>w</mi> <mi>x</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>23</mn> <mo>)</mo> </mrow> </mrow>
Wherein,
<mrow> <msub> <mi>F</mi> <mi>x</mi> </msub> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>1</mn> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>1</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mrow> <mo>-</mo> <msub> <mi>&amp;lambda;</mi> <mi>x</mi> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>,</mo> <msub> <mi>G</mi> <mi>x</mi> </msub> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>1</mn> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>24</mn> <mo>)</mo> </mrow> </mrow>
After carrying out discretization to formula (23) with certain period Δ t, the state equation obtained from kth step to the step of kth+1 is
Xx(k+1)=ΦxXx(k)+ωx(k) (25)
Wherein,
<mrow> <msub> <mi>&amp;Phi;</mi> <mi>x</mi> </msub> <mo>=</mo> <mi>exp</mi> <mrow> <mo>(</mo> <msub> <mi>F</mi> <mi>x</mi> </msub> <mi>&amp;Delta;</mi> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mn>1</mn> </mtd> <mtd> <mrow> <mi>&amp;Delta;</mi> <mi>t</mi> </mrow> </mtd> <mtd> <mrow> <mn>1</mn> <mo>/</mo> <msubsup> <mi>&amp;lambda;</mi> <mi>x</mi> <mn>2</mn> </msubsup> <mrow> <mo>(</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msub> <mi>&amp;lambda;</mi> <mi>x</mi> </msub> <mi>&amp;Delta;</mi> <mi>t</mi> </mrow> </msup> <mo>+</mo> <msub> <mi>&amp;lambda;</mi> <mi>x</mi> </msub> <mi>&amp;Delta;</mi> <mi>t</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>1</mn> </mtd> <mtd> <mrow> <mn>1</mn> <mo>/</mo> <msub> <mi>&amp;lambda;</mi> <mi>x</mi> </msub> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msub> <mi>&amp;lambda;</mi> <mi>x</mi> </msub> <mi>&amp;Delta;</mi> <mi>t</mi> </mrow> </msup> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msub> <mi>&amp;lambda;</mi> <mi>x</mi> </msub> <mi>&amp;Delta;</mi> <mi>t</mi> </mrow> </msup> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>26</mn> <mo>)</mo> </mrow> </mrow>
ωx(k) x-axis subsystem state noise vector is represented, is zero mean Gaussian white noise vector；
The position of hypersonic gliding bullet is obtained by Ground Target Tracking measurement device, then measures equation and be
ηx=x+ υx (27)
Wherein, υxX-axis subsystem measurement noise vector is represented, is zero mean Gaussian white noise vector；
Similarly, for system (19), y-axis subsystem state variable X is definedy=[y vy ay]T, then following state equation is obtained
<mrow> <msub> <mover> <mi>X</mi> <mo>&amp;CenterDot;</mo> </mover> <mi>y</mi> </msub> <mo>=</mo> <msub> <mi>F</mi> <mi>y</mi> </msub> <msub> <mi>X</mi> <mi>y</mi> </msub> <mo>+</mo> <msub> <mi>G</mi> <mi>y</mi> </msub> <msub> <mi>w</mi> <mi>y</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>28</mn> <mo>)</mo> </mrow> </mrow>
Wherein,
<mrow> <msub> <mi>F</mi> <mi>y</mi> </msub> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>1</mn> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>1</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mrow> <mo>-</mo> <msub> <mi>&amp;lambda;</mi> <mi>y</mi> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>,</mo> <msub> <mi>G</mi> <mi>y</mi> </msub> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>1</mn> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>29</mn> <mo>)</mo> </mrow> </mrow>
After carrying out discretization to formula (28) with certain period Δ t, the state equation obtained from kth step to the step of kth+1 is
Xy(k+1)=ΦyXy(k)+ωy(k) (30)
Wherein,
<mrow> <msub> <mi>&amp;Phi;</mi> <mi>y</mi> </msub> <mo>=</mo> <mi>exp</mi> <mrow> <mo>(</mo> <msub> <mi>F</mi> <mi>y</mi> </msub> <mi>&amp;Delta;</mi> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mn>1</mn> </mtd> <mtd> <mrow> <mi>&amp;Delta;</mi> <mi>t</mi> </mrow> </mtd> <mtd> <mrow> <mn>1</mn> <mo>/</mo> <msubsup> <mi>&amp;lambda;</mi> <mi>y</mi> <mn>2</mn> </msubsup> <mrow> <mo>(</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msub> <mi>&amp;lambda;</mi> <mi>y</mi> </msub> <mi>&amp;Delta;</mi> <mi>t</mi> </mrow> </msup> <mo>+</mo> <msub> <mi>&amp;lambda;</mi> <mi>y</mi> </msub> <mi>&amp;Delta;</mi> <mi>t</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>1</mn> </mtd> <mtd> <mrow> <mn>1</mn> <mo>/</mo> <msub> <mi>&amp;lambda;</mi> <mi>y</mi> </msub> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msub> <mi>&amp;lambda;</mi> <mi>y</mi> </msub> <mi>&amp;Delta;</mi> <mi>t</mi> </mrow> </msup> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msub> <mi>&amp;lambda;</mi> <mi>y</mi> </msub> <mi>&amp;Delta;</mi> <mi>t</mi> </mrow> </msup> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>31</mn> <mo>)</mo> </mrow> </mrow>
ωy(k) y-axis subsystem state noise vector is represented, is zero mean Gaussian white noise vector；
Measuring equation is
ηy=y+ υy (32)
Wherein, υyY-axis subsystem measurement noise vector is represented, is zero mean Gaussian white noise vector；
For system (20), z-axis subsystem state variable X is definedz=[z vz az]T, then following state equation is obtained
<mrow> <msub> <mover> <mi>X</mi> <mo>&amp;CenterDot;</mo> </mover> <mi>z</mi> </msub> <mo>=</mo> <msub> <mi>F</mi> <mi>z</mi> </msub> <msub> <mi>X</mi> <mi>z</mi> </msub> <mo>+</mo> <msub> <mi>G</mi> <mi>z</mi> </msub> <msub> <mi>w</mi> <mi>z</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>33</mn> <mo>)</mo> </mrow> </mrow>
Wherein,
<mrow> <msub> <mi>F</mi> <mi>z</mi> </msub> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>1</mn> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>1</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mrow> <mo>-</mo> <msub> <mi>&amp;lambda;</mi> <mi>z</mi> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>,</mo> <msub> <mi>G</mi> <mi>z</mi> </msub> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>1</mn> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>34</mn> <mo>)</mo> </mrow> </mrow>
After carrying out discretization to formula (28) with certain period Δ t, the state equation obtained from kth step to the step of kth+1 is
Xz(k+1)=ΦzXz(k)+ωz(k) (35)
Wherein,
<mrow> <msub> <mi>&amp;Phi;</mi> <mi>z</mi> </msub> <mo>=</mo> <mi>exp</mi> <mrow> <mo>(</mo> <msub> <mi>F</mi> <mi>z</mi> </msub> <mi>&amp;Delta;</mi> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mn>1</mn> </mtd> <mtd> <mrow> <mi>&amp;Delta;</mi> <mi>t</mi> </mrow> </mtd> <mtd> <mrow> <mn>1</mn> <mo>/</mo> <msubsup> <mi>&amp;lambda;</mi> <mi>z</mi> <mn>2</mn> </msubsup> <mrow> <mo>(</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msub> <mi>&amp;lambda;</mi> <mi>z</mi> </msub> <mi>&amp;Delta;</mi> <mi>t</mi> </mrow> </msup> <mo>+</mo> <msub> <mi>&amp;lambda;</mi> <mi>z</mi> </msub> <mi>&amp;Delta;</mi> <mi>t</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>1</mn> </mtd> <mtd> <mrow> <mn>1</mn> <mo>/</mo> <msub> <mi>&amp;lambda;</mi> <mi>z</mi> </msub> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msub> <mi>&amp;lambda;</mi> <mi>z</mi> </msub> <mi>&amp;Delta;</mi> <mi>t</mi> </mrow> </msup> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msub> <mi>&amp;lambda;</mi> <mi>z</mi> </msub> <mi>&amp;Delta;</mi> <mi>t</mi> </mrow> </msup> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>36</mn> <mo>)</mo> </mrow> </mrow>
ωz(k) z-axis subsystem state noise vector is represented, is zero mean Gaussian white noise vector；
Measuring equation is
ηz=z+ υz (37)
Wherein, υzZ-axis subsystem measurement noise vector is represented, is zero mean Gaussian white noise vector；
The subsystem on three axles of sensing point inertial coodinate system, i.e., formula (23) and (27) are formed x-axis subsystem, formula (28) and (32) the y-axis subsystem formed, and the z-axis subsystem that formula (33) and (37) are formed, separately design Kalman filter, its is pre- Report equation is formula (4), i.e.,
<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mover> <mi>X</mi> <mo>&amp;OverBar;</mo> </mover> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>&amp;Phi;</mi> <mi>i</mi> </msub> <msub> <mover> <mi>X</mi> <mo>^</mo> </mover> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>P</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>/</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>&amp;Phi;</mi> <mi>i</mi> </msub> <msub> <mi>P</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <msubsup> <mi>&amp;Phi;</mi> <mi>i</mi> <mi>T</mi> </msubsup> <mo>+</mo> <msub> <mi>Q</mi> <mi>i</mi> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>,</mo> <mi>i</mi> <mo>=</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> </mrow>
In above formula,The state estimation of wave filter i kth step is represented,Wave filter i is represented from kth step to kth+1 The state forecast value of step, Pi(k) the state estimation error co-variance matrix of wave filter i kth step, P are representedi(k+1/k) filtering is represented Device i walks the state forecast error co-variance matrix to the step of kth+1, Q from kthiFor wave filter i model prediction error covariance squares Battle array,
The measurement update equation of wave filter is formula (6), i.e.,
<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>K</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>P</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>/</mo> <mi>k</mi> <mo>)</mo> </mrow> <msubsup> <mi>H</mi> <mi>i</mi> <mi>T</mi> </msubsup> <msup> <mrow> <mo>&amp;lsqb;</mo> <msub> <mi>H</mi> <mi>i</mi> </msub> <msub> <mi>P</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>/</mo> <mi>k</mi> <mo>)</mo> </mrow> <msubsup> <mi>H</mi> <mi>i</mi> <mi>T</mi> </msubsup> <mo>+</mo> <msub> <mi>R</mi> <mi>i</mi> </msub> <mo>&amp;rsqb;</mo> </mrow> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mover> <mi>X</mi> <mo>^</mo> </mover> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mover> <mi>X</mi> <mo>&amp;OverBar;</mo> </mover> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>K</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>&amp;lsqb;</mo> <msub> <mi>&amp;eta;</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>H</mi> <mi>i</mi> </msub> <msub> <mover> <mi>X</mi> <mo>&amp;OverBar;</mo> </mover> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>P</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>=</mo> <mo>&amp;lsqb;</mo> <mi>I</mi> <mo>-</mo> <msub> <mi>K</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <msub> <mi>H</mi> <mi>i</mi> </msub> <mo>&amp;rsqb;</mo> <msub> <mi>P</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>/</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>,</mo> <mi>i</mi> <mo>=</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> </mrow>
Wherein, subscriptTRepresent matrix transposition computing, subscript-1Matrix inversion operation is represented, wave filter i calculation matrix is formula (7), i.e.,
Hi=[1 0 0], i=x, y, z
RiFor wave filter i measurement error covariance matrixes, I represents unit square, Ki(k+1) filtering of the step of wave filter i kth+1 is represented Gain matrix, Pi(k+1) the state estimation error co-variance matrix of the step of wave filter i kth+1, η are representedi(k+1) wave filter i is represented The metrical information of the step of kth+1, i.e. ηx(k+1)=x (k+1)+υx(k+1), ηy(k+1)=y (k+1)+υy(k+1),ηz(k+1)=z (k+1)+υz(k+1)；
When no Ground Target Tracking measurement device information, prognostic equation (4) is only run, there is Ground Target Tracking measurement device During information, then prognostic equation (4) and measurement update equation (6) are run simultaneously；
Step (3 1) in step 3 is specific as follows：
First, the hypersonic gliding bullet of tracking finish time kT time space is obtained according to the estimated result of three subfilters to exist Position [x y z] in sensing point inertial systemT, the position in sensing point inertial coodinate system and the position relationship in Earth central inertial system For
[xI yI zI]T=[x y z]T+[0 Re 0]T (38)
Wherein, Re=6378km is earth mean radius；
According to the acceleration estimation value [a of the hypersonic gliding bullet of tracking finish timex ay az]T, due to sensing point inertia Three of coordinate system are parallel with three axles of Earth central inertial system, can obtain
[axI ayI azI]T=[ax ay az]T (39)
Wherein, axI、ayIAnd azIRepresent projection of the acceleration of hypersonic gliding bullet in Earth central inertial system
The aerodynamic force that currently hypersonic gliding bullet is subject to is calculated with following formula：
<mrow> <mi>R</mi> <mo>=</mo> <mi>m</mi> <mo>(</mo> <mrow> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mi>a</mi> <mrow> <mi>x</mi> <mi>I</mi> </mrow> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>a</mi> <mrow> <mi>y</mi> <mi>I</mi> </mrow> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>a</mi> <mrow> <mi>z</mi> <mi>I</mi> </mrow> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mi>g</mi> <mrow> <mi>x</mi> <mi>I</mi> </mrow> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>g</mi> <mrow> <mi>y</mi> <mi>I</mi> </mrow> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>g</mi> <mrow> <mi>z</mi> <mi>I</mi> </mrow> </msub> </mtd> </mtr> </mtable> </mfenced> </mrow> <mo>)</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>40</mn> <mo>)</mo> </mrow> </mrow>
When 15Mach or so, lift-drag ratio is about 3:1, due to Rz1=0, it can obtain
<mrow> <msub> <mi>R</mi> <mrow> <mi>y</mi> <mn>1</mn> </mrow> </msub> <mo>=</mo> <mfrac> <mn>3</mn> <msqrt> <mn>10</mn> </msqrt> </mfrac> <mo>|</mo> <mo>|</mo> <mi>R</mi> <mo>|</mo> <mo>|</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>41</mn> <mo>)</mo> </mrow> </mrow>
According to the velocity estimation value [v of the hypersonic gliding bullet of tracking finish timex vy vz]T, because sensing point inertia is sat Three of mark system are parallel with three axles of Earth central inertial system, can obtain
[vxI vyI vzI]T=[vx vy vz]T (42)
According to the speed [v of hypersonic gliding bullet at the end of trackingxI vyI vzI]T, its trajectory is calculated with formula (12) and (13) Angle θ and ψv；
Step (three or two) in step 3 is specific as follows：
According to the speed [v of hypersonic gliding bullet at the end of trackingxI vyI vzI]T, its speed v is calculated with formula (11), then Current dynamic head q is tried to achieve using formula (9)；
Calculating current lift coefficient with following formula is
<mrow> <msub> <mi>c</mi> <mi>y</mi> </msub> <mo>=</mo> <mfrac> <msub> <mi>R</mi> <mrow> <mi>y</mi> <mn>1</mn> </mrow> </msub> <mrow> <mi>q</mi> <mi>S</mi> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>43</mn> <mo>)</mo> </mrow> </mrow>
When the angle of attack is 10 °, lift coefficient 0.37, so estimating the current angle of attack with following formula：
<mrow> <mi>&amp;alpha;</mi> <mo>=</mo> <mfrac> <msub> <mi>c</mi> <mi>y</mi> </msub> <mn>0.37</mn> </mfrac> <mfrac> <mn>10</mn> <mn>57.3</mn> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>44</mn> <mo>)</mo> </mrow> </mrow>
Hypersonic gliding bullet uses banked turn control mode, i.e. yaw angle β=0, then by preceding two formula in formula (14) The angle of pitch can be calculatedAnd yaw angle ψ；
Estimate current roll angle γ：Projection of the aerodynamic force in inertial system has been calculated by formula (40), has projected this at mesh Standard type coordinate system, i.e.,
<mrow> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mi>R</mi> <mrow> <mi>x</mi> <mn>1</mn> </mrow> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>R</mi> <mrow> <mi>y</mi> <mn>1</mn> </mrow> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>R</mi> <mrow> <mi>z</mi> <mn>1</mn> </mrow> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>=</mo> <msub> <mi>C</mi> <mrow> <mi>t</mi> <mn>0</mn> <mo>&amp;RightArrow;</mo> <mn>1</mn> </mrow> </msub> <msub> <mi>C</mi> <mrow> <mn>0</mn> <mo>&amp;RightArrow;</mo> <mi>t</mi> <mn>0</mn> </mrow> </msub> <mi>R</mi> </mrow>
Above formula can further be write
Order
It can then be obtained according to formula (45)
<mrow> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mi>b</mi> <mn>3</mn> </msub> </mtd> <mtd> <msub> <mi>b</mi> <mn>2</mn> </msub> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>-</mo> <msub> <mi>b</mi> <mn>2</mn> </msub> </mrow> </mtd> <mtd> <msub> <mi>b</mi> <mn>3</mn> </msub> </mtd> </mtr> </mtable> </mfenced> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <mi>s</mi> <mi>i</mi> <mi>n</mi> <mi>&amp;gamma;</mi> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>cos</mi> <mi>&amp;gamma;</mi> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mi>R</mi> <mrow> <mi>y</mi> <mn>1</mn> </mrow> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>R</mi> <mrow> <mi>z</mi> <mn>1</mn> </mrow> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>46</mn> <mo>)</mo> </mrow> </mrow>
Wherein, Ry1Estimated by formula (41), Rz1=0, sin γ and cos γ are tried to achieve with formula (46), current rolling is determined with following formula Corner
γ=tan2-1(sinγ,cosγ) (47)。
CN201510220893.5A 2015-05-04 2015-05-04 A kind of hypersonic gliding bullet Skipping trajectory Forecasting Methodology of near space Active CN104778376B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201510220893.5A CN104778376B (en) 2015-05-04 2015-05-04 A kind of hypersonic gliding bullet Skipping trajectory Forecasting Methodology of near space

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201510220893.5A CN104778376B (en) 2015-05-04 2015-05-04 A kind of hypersonic gliding bullet Skipping trajectory Forecasting Methodology of near space

Publications (2)

Publication Number Publication Date
CN104778376A CN104778376A (en) 2015-07-15
CN104778376B true CN104778376B (en) 2018-03-16

Family

ID=53619836

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201510220893.5A Active CN104778376B (en) 2015-05-04 2015-05-04 A kind of hypersonic gliding bullet Skipping trajectory Forecasting Methodology of near space

Country Status (1)

Country Link
CN (1) CN104778376B (en)

Families Citing this family (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105676638B (en) * 2016-01-11 2018-06-29 北京航空航天大学 Steady gliding/quasi- natural frequency jump gliding combined maneuver is dashed forward ballistic planing method
CN105760573B (en) * 2016-01-21 2019-07-02 中国工程物理研究院总体工程研究所 Along the disturbance gravitation extension approximation method of a wide range of Maneuver Ballistic Trajectory of near space
CN109948304B (en) * 2019-04-17 2022-07-22 哈尔滨工业大学 Method for predicting near space hypersonic aerocraft AHW trajectory
CN110065649B (en) * 2019-05-10 2022-06-07 哈尔滨工业大学 Method for designing near space hypersonic speed aircraft trajectory by adopting virtual aiming point
CN110320510B (en) * 2019-06-14 2022-06-24 南京理工大学 Ballistic missile structure parameter estimation method based on centroid height parameter elimination
CN110990765B (en) * 2019-12-03 2022-07-22 南京理工大学 Method and system for calculating miss distance based on trajectory equation
CN111239722B (en) * 2020-02-12 2023-05-05 哈尔滨工业大学 Tracking algorithm for maneuvering mutation of near space high-speed maneuvering target
CN111428343A (en) * 2020-02-26 2020-07-17 北京电子工程总体研究所 Method for estimating state of glider in near space, storage medium and computing device
CN111783358B (en) * 2020-07-02 2022-10-04 哈尔滨工业大学 Bayesian estimation-based long-term trajectory prediction method for hypersonic aircraft
JP7394802B2 (en) 2021-02-19 2023-12-08 三菱電機株式会社 Gliding flying object identification method, flying object tracking system, flying object countermeasure system, and ground system
CN116909303B (en) * 2023-07-14 2024-02-02 中国人民解放军国防科技大学 Process noise self-adaptive adjusting method for near space target tracking

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103528587A (en) * 2013-10-15 2014-01-22 西北工业大学 Autonomous integrated navigation system
CN103884237A (en) * 2014-04-08 2014-06-25 哈尔滨工业大学 Several-for-one collaborative guidance method based on target probability distribution information

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6536350B2 (en) * 2001-03-07 2003-03-25 The United States Of America As Represented By The United States Department Of Energy Stagnation pressure activated fuel release mechanism for hypersonic projectiles
US8436283B1 (en) * 2008-07-11 2013-05-07 Davidson Technologies Inc. System and method for guiding and controlling a missile using high order sliding mode control

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103528587A (en) * 2013-10-15 2014-01-22 西北工业大学 Autonomous integrated navigation system
CN103884237A (en) * 2014-04-08 2014-06-25 哈尔滨工业大学 Several-for-one collaborative guidance method based on target probability distribution information

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
《临近空间非弹道式目标HTV-2跟踪滤波与预报问题》;秦雷 等;《航天控制》;20150430;第33卷(第2期);第2-5节 *

Also Published As

Publication number Publication date
CN104778376A (en) 2015-07-15

Similar Documents

Publication Publication Date Title
CN104778376B (en) A kind of hypersonic gliding bullet Skipping trajectory Forecasting Methodology of near space
CN107544067B (en) Hypersonic reentry vehicle tracking method based on Gaussian mixture approximation
CN105717942B (en) A kind of unmanned vehicle space barrier-avoiding method and the online planing method of introductory path
CN104035335B (en) Steady glide reentry guidance method based on the longitudinal and transverse journey analytical Prediction method of high accuracy
CN109144084B (en) A kind of VTOL Reusable Launch Vehicles Attitude tracking control method based on set time Convergence monitoring device
CN102706217B (en) Method for controlling attack angle and attack time of multiple missiles
CN106586033A (en) Adaptive segmentation multistage linear spectrum generalized standard control missdistance reentry guidance method
CN108153323B (en) A kind of high-altitude unmanned vehicle high-precision reentry guidance method
CN106681348A (en) Guidance and control integrated design method considering all-strapdown seeker view field constraint
CN105182985A (en) Hypersonic flight vehicle dive segment full amount integration guidance control method
CN104777844B (en) Method for tracking trajectories of hypersonic velocity near space aircraft
CN105486307B (en) For the line-of-sight rate by line method of estimation of maneuvering target
CN105180728B (en) Front data based rapid air alignment method of rotary guided projectiles
CN109446582A (en) A kind of high-precision depression of order considering earth rotation steadily glides dynamic modeling method
CN108153330A (en) Unmanned aerial vehicle three-dimensional track self-adaptive tracking method based on feasible region constraint
CN105486308B (en) Estimation plays the design method of the rapid convergence Kalman filter of line of sight angular speed
Bao et al. Integrated guidance and control for hypersonic morphing missile based on variable span auxiliary control
CN110187713A (en) A kind of longitudinally controlled method of hypersonic aircraft based on aerodynamic parameter on-line identification
CN105115508A (en) Post data-based rotary guided projectile quick air alignment method
CN107368085A (en) Model prediction-based method for controlling height of stratospheric airship in wind field
CN113847913A (en) Missile-borne integrated navigation method based on ballistic model constraint
CN106802570A (en) A kind of method and apparatus of depopulated helicopter position tracking
CN116663263A (en) Aircraft tail end falling speed and falling angle constraint guidance method and system based on sliding die surface
CN110162818A (en) Parachute-bomb ballistic calculation
CN109445283B (en) Control method for fixed-point tracking of under-actuated aerostat on plane

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
EXSB Decision made by sipo to initiate substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant