CN113075652B - Three-dimensional tracking method for hypersonic aircraft - Google Patents

Three-dimensional tracking method for hypersonic aircraft Download PDF

Info

Publication number
CN113075652B
CN113075652B CN202110212919.7A CN202110212919A CN113075652B CN 113075652 B CN113075652 B CN 113075652B CN 202110212919 A CN202110212919 A CN 202110212919A CN 113075652 B CN113075652 B CN 113075652B
Authority
CN
China
Prior art keywords
measurement
trajectory
value
model
aircraft
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
CN202110212919.7A
Other languages
Chinese (zh)
Other versions
CN113075652A (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.)
Air Force Early Warning Academy
Original Assignee
Air Force Early Warning Academy
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 Air Force Early Warning Academy filed Critical Air Force Early Warning Academy
Priority to CN202110212919.7A priority Critical patent/CN113075652B/en
Publication of CN113075652A publication Critical patent/CN113075652A/en
Application granted granted Critical
Publication of CN113075652B publication Critical patent/CN113075652B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/66—Radar-tracking systems; Analogous systems
    • G01S13/72—Radar-tracking systems; Analogous systems for two-dimensional tracking, e.g. combination of angle and range tracking, track-while-scan radar
    • G01S13/723—Radar-tracking systems; Analogous systems for two-dimensional tracking, e.g. combination of angle and range tracking, track-while-scan radar by using numerical data
    • G01S13/726—Multiple target tracking
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/93—Radar or analogous systems specially adapted for specific applications for anti-collision purposes
    • G01S13/933—Radar or analogous systems specially adapted for specific applications for anti-collision purposes of aircraft or spacecraft

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Physics & Mathematics (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Physics & Mathematics (AREA)
  • Electromagnetism (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

The invention discloses a hypersonic aircraft three-dimensional tracking method, which is applied to the field of aircraft trajectory tracking and comprises the following specific steps: establishing a motion equation of the hypersonic aircraft in the glide stage under a half-speed coordinate system; dividing a maneuvering trajectory of the hypersonic aircraft into a longitudinal trajectory and a transverse trajectory, modeling acceleration on the longitudinal trajectory into an attenuation oscillation model for tracking, and tracking on the transverse trajectory by adopting an interactive multi-model formed by a Singer model and a uniform acceleration model; measuring information of the radar under a spherical coordinate system is obtained, and the measuring information is converted into a Cartesian coordinate system by unbiased measurement conversion; and finally, designing a self-adaptive Kalman filtering method of multiple fading factors on the basis of the conventional Kalman filtering algorithm, and enhancing the self-adaptive tracking capability of the model to the strong maneuver of the aircraft. Compared with the existing method, the method effectively improves the tracking precision and enhances the adaptability of the model to strong maneuverability when the hypersonic aircraft is tracked in three dimensions.

Description

Three-dimensional tracking method for hypersonic aircraft
Technical Field
The invention relates to the technical field of aircraft trajectory tracking, in particular to a three-dimensional tracking method for a hypersonic aircraft.
Background
The hypersonic aircraft usually flies in a 20-100 km near airspace, has the flying speed of more than 5Ma, is a novel weapon with strategic deterrence and tactical attack capabilities, and has the characteristics of high flying height, high flying speed, strong defense capability, large attack threat and the like. The development of the hypersonic aircraft technology leads the aerospace safety of all countries to face new challenges, and the early warning, stable tracking and defense interception of the hypersonic aircraft become urgent requirements. Due to the high maneuvering characteristic of the hypersonic aircraft, the phenomena of accuracy reduction and even filtering divergence easily occur during tracking by adopting a conventional tracking model and Kalman filtering. Therefore, the establishment of a maneuvering model according with the motion characteristics of the hypersonic aircraft and the design of a more robust filtering algorithm are feasible paths for improving the tracking precision of the hypersonic aircraft and realizing stable tracking.
The hypersonic aircraft tracking mainly comprises two parts, namely a maneuvering model and a filtering algorithm. The maneuvering models of the existing hypersonic flight vehicle are mainly divided into a single model and a multi-model. The single model mainly simplifies the maneuvering of the hypersonic aircraft into uniform motion, uniform accelerated motion, circular motion, sine wave motion and the like by performing different maneuvering statistical modeling on the acceleration, but the single model cannot depict maneuvering patterns of the hypersonic aircraft in three dimensions, and the tracking error is large; the idea of multiple models is mainly to combine multiple different single models into a model set, and to cover all possible maneuvering patterns of the hypersonic aircraft as much as possible by using the model set so as to realize the tracking effect superior to that of the single model. Meanwhile, in a filtering algorithm, the traditional Kalman filtering or nonlinear filtering cannot correct system prediction according to the maneuvering condition of the hypersonic aircraft, and has no self-adaptive adjustment capability, so that when the aircraft maneuvers, the filtering error is fluctuated greatly, and even the filtering divergence is caused. Most of the existing solutions of adaptive filtering algorithms introduce fading factors, but because hypersonic aircrafts are three-dimensional maneuvering targets, single fading factors cannot be accurately adjusted in three dimensions respectively, and the effect is poor.
Therefore, the three-dimensional motion trajectory of the hypersonic aircraft is decomposed into a longitudinal trajectory and a transverse trajectory by fully combining the maneuvering characteristics of the hypersonic aircraft, and tracking models suitable for the motion characteristics of the two trajectories are respectively established. Secondly, a correction threshold is set based on the measurement noise covariance and the residual error, multiple fading factors are introduced to adjust the weight of each channel in real time, and an exponential weighting innovation covariance estimation method for limiting the memory length is designed to carry out limited memory on old measurement information, so that the influence of recent measurement data is improved.
Disclosure of Invention
The invention provides a hypersonic aircraft three-dimensional tracking method, aiming at solving the problems that the maneuvering mode of a hypersonic aircraft in the near space is complex, and a single kinematics model and a conventional filtering algorithm can not realize stable and high-precision three-dimensional tracking in the prior art, effectively improving the tracking precision and enhancing the adaptability of the model to strong maneuvering.
In order to achieve the above object, the present invention provides a method for three-dimensional tracking of a hypersonic flight vehicle, which is characterized in that the method comprises:
the method comprises the following steps: establishing a motion equation of the hypersonic aircraft in the glide section under a half-speed coordinate system;
step two: the maneuvering trajectory of the hypersonic aircraft is divided into a longitudinal maneuvering trajectory and a transverse maneuvering trajectory:
modeling a motion state vector of the hypersonic aircraft in the longitudinal direction, wherein a discrete form longitudinal ballistic model of the hypersonic aircraft is represented as:
X(k+1)＝F DO X(k)+W(k)
wherein X (k +1) and X (k) are state variables at time k and time k +1, respectively, and F DO A vertical state transition matrix, W (k) process noise,
Figure BDA0002952076000000021
wherein each variable respectively represents the position, the speed, the acceleration and the jerk of the aircraft on the longitudinal trajectory at the moment k;
the modeling of the transverse trajectory of the hypersonic aerocraft by adopting a model set consisting of a uniform acceleration model and a Singer model in the transverse direction is represented as follows:
X(k+1)＝F CA X(k)+W(k)
wherein X (k +1) and X (k) are state variables at time k and time k +1, respectively, and F CA W (k) is process noise,
Figure BDA0002952076000000031
wherein each variable respectively represents the position, the speed and the acceleration of the aircraft in the x and y directions on the transverse trajectory at the moment k;
step three: obtaining radar measurement, and converting measurement data in a spherical coordinate system into a Cartesian coordinate system by adopting correction unbiased measurement conversion to obtain a converted measurement value;
step four: and (3) performing tracking filtering on the model predicted value obtained in the step two and the measurement value obtained in the step three by adopting a self-adaptive Kalman filtering method of multiple fading factors, wherein the method specifically comprises the following steps: establishing a filtering equation, setting a correction threshold value, selecting multiple fading factors, estimating covariance of an innovation sequence and tracking filtering.
Preferably, in the first step, the equation of motion of the hypersonic aircraft in the glide section under the half-speed coordinate system is as follows:
Figure BDA0002952076000000032
wherein r, lambda, phi, V, theta and sigma respectively represent geocentric distance, longitude, latitude, speed inclination angle and azimuth angle; upsilon is c Represents the roll angle; a is D 、a L Respectively, the resistance acceleration and the lift acceleration.
Preferably, the longitudinal state transition matrix F in the second step DO Is composed of
Figure BDA0002952076000000041
In the formula,
Figure BDA0002952076000000042
Figure BDA0002952076000000043
alpha is the maximum correlation attenuation and has a value range of
Figure BDA0002952076000000044
Beta is motor-driven oscillationThe frequency, the span of values is between (0,1), and T is the sampling interval.
Preferably, the transverse state transition matrix F in the second step CA Is composed of
Figure BDA0002952076000000045
Where T is the sampling interval.
Preferably, the measured value in the third step is
Figure BDA0002952076000000051
Wherein A is an azimuth; e is a pitch angle; r is m Is the target distance;
Figure BDA0002952076000000052
Figure BDA0002952076000000053
in order to square the error in the pitch angle,
Figure BDA0002952076000000054
Figure BDA0002952076000000055
is the square of the azimuthal error.
Preferably, the covariance of the measurement noise in the measurement values is
Figure BDA0002952076000000056
Wherein R is a symmetric matrix, and each element in R after correction unbiased measurement conversion is
Figure BDA0002952076000000057
In the formula,
Figure BDA0002952076000000058
Figure BDA0002952076000000059
is the square of the range error.
Preferably, the filter equation established in step four is
X(k)＝F(k-1)X(k-1)+Γ(k-1)W(k-1)
Z(k)＝H(k)X(k)+V(k)
Wherein, X (k) and Z (k) respectively represent a state vector and a measurement vector at the target k moment; f (k) represents a state transition matrix at time k; h (k) represents the measurement matrix at time k; Γ (k-1) is the system noise matrix at time k-1; w (k) and v (k) represent the state noise and the measurement noise at time k, respectively.
Preferably, the modified threshold value in step four is set to
Figure BDA0002952076000000061
Wherein chi is a threshold coefficient and has a value range of 0-1, R i,i (k) Representing the measurement noise at time k for the measurement vector at row i and column i.
Preferably, the multiple fading factors ρ (k) selected in step four adjust each channel in real time according to the actual filtering condition of each channel:
ρ(k)＝diag{ρ 1 ,ρ 2 ,…,ρ i ,…,ρ n }
wherein n is the state vector dimension; ρ (i) represents an extinction factor corresponding to the ith-dimension state vector; the multiple fading factors are calculated in the following way
Figure BDA0002952076000000062
In the formula, H i (k) Representing the ith element in the measurement matrix at the k moment; when rho (k) is I, the multiple evanescent filtering degenerates to conventional Kalman filtering; rho i (k) Value of (a) is evaluated by an innovation covariance estimate
Figure BDA0002952076000000063
The influence of the impact is such that,
Figure BDA0002952076000000064
the estimate directly affects p i (k) Has a regulating effect of
Figure BDA0002952076000000065
Is influenced by the measured value, and is thus given by rho i (k) And setting the upper value limit as a real number larger than 1 to avoid filtering oscillation.
Preferably, the recurrence formula of the innovation covariance estimate is:
Figure BDA0002952076000000066
the size of the innovation covariance estimation value depends on the values of the memory length N and the forgetting factor b, wherein the selection of the values of N and b depends on the intensity of the change of the system state, and for a system with more gradual change, the value of the memory length N is larger, and the value of the forgetting factor b is also larger, so that the innovation covariance estimation value is obtained.
Aiming at the maneuvering characteristics of the hypersonic aerocraft in the glide section, the invention decomposes the three-dimensional motion trajectory into a longitudinal trajectory and a transverse trajectory, tracking models suitable for the motion characteristics of the two trajectories are respectively established, so that the problem that the real motion state of the target cannot be reflected due to the mismatch of the tracking models when the maneuvering mode of the hypersonic aerocraft is changed is solved, and aiming at the phenomenon that the accuracy of the conventional Kalman filtering is easy to decline when the tracking error is larger, the adaptive Kalman filtering with multiple fading factors is constructed, the fading factor can be adjusted in real time according to the filtering divergence condition, the predicted value is corrected by using the latest measurement value, therefore, the tracking effect is superior to that of the conventional Kalman filtering, stable filtering performance can still be kept after filtering convergence, the convergence speed and the tracking precision are improved, and finally, simulation experiments show that the method provided by the invention has higher precision and stability when the hypersonic aircraft is tracked.
Drawings
FIG. 1 is a flow chart of a hypersonic aircraft three-dimensional tracking algorithm based on adaptive filtering according to the invention;
FIG. 2 is a relationship between a spherical observation coordinate system and a Cartesian coordinate system of a radar;
FIG. 3 is a hypersonic aircraft maneuver trajectory;
FIG. 4 is a tracking error diagram of the longitudinal trajectory of the glide section of the hypersonic aerocraft for 0-300 s, which is obtained by the method in a simulation experiment;
FIG. 5 is a tracking error diagram of the transverse trajectory of the glide section of the hypersonic aerocraft for 0-300 s, which is obtained by the method in a simulation experiment;
FIG. 6 is a graph of total tracking error of 0-300 s in the glide section of a hypersonic flight vehicle obtained by the method in a simulation experiment;
FIG. 7 is a velocity estimation diagram of 0-300 s of a glide section of a hypersonic flight vehicle obtained by the method in a simulation experiment;
FIG. 8 is a graph of the total velocity error of the hypersonic flight vehicle in 0-300 s of the glide section obtained by the method in a simulation experiment;
FIG. 9 shows a roll angle control mode of the glide section of the hypersonic aircraft for 0-300 s in a simulation experiment;
FIG. 10 shows the attack angle control mode of the hypersonic flight vehicle in the glide section for 0-300 s in the simulation experiment.
Detailed Description
The invention is described in further detail below with reference to the figures and specific embodiments.
The invention provides a hypersonic aircraft three-dimensional tracking method, which is realized by adopting the following technical scheme:
the method comprises the following steps: establishing a motion equation of the hypersonic aircraft in a glide section under a half-speed coordinate system;
step two: the maneuvering trajectory of the hypersonic aircraft is divided into a longitudinal maneuvering trajectory and a transverse maneuvering trajectory. On a longitudinal trajectory, according to the acceleration attenuation oscillation characteristic, modeling the acceleration attenuation oscillation characteristic into a zero-mean random process of second-order time autocorrelation; on a transverse trajectory, the acceleration change is relatively simple, and a Singer model and a uniform acceleration model are adopted to construct a model set;
step three: acquiring measurement data of the radar in a spherical coordinate system, and converting the measurement data in the spherical coordinate system into a Cartesian coordinate system by correction unbiased measurement conversion to obtain an observed value of the hypersonic aircraft;
step four: a self-adaptive Kalman filtering method based on multiple fading factors is designed, the process noise covariance weight of each channel is self-adaptively adjusted, a divergence degree judgment threshold is set, the influence of measurement errors is reduced, the self-adaptive tracking capability of a filtering algorithm is enhanced, and self-adaptive tracking filtering is performed based on the tracking model established in the second step and the observed value obtained in the third step.
In one embodiment of the invention, the following embodiments are employed:
step one, establishing a motion equation of the hypersonic aircraft in a glide section
The influence of earth rotation and oblateness is not considered, and the motion equation of the hypersonic aerocraft is established under a half-speed coordinate system as follows:
Figure BDA0002952076000000081
wherein r, λ, φ, V, θ, σ represent the geocentric distance, longitude, latitude, velocity inclination angle and azimuth angle, respectively. Alpha is alpha c Representing angle of attack, u c Representing the roll angle, wherein the attack angle and the roll angle are aircraft control variables; g represents the gravitational acceleration. a is D 、a L Respectively, the resistance acceleration and the lift acceleration. The expressions of lift acceleration and drag acceleration are
Figure BDA0002952076000000091
Wherein, C L (α)、C D And (alpha) is a lift coefficient and a drag coefficient respectively, the two coefficients are functions related to an attack angle, S is a target equivalent sectional area, m is a target mass, and rho is an atmospheric density.
Step two, respectively modeling the longitudinal maneuvering trajectory and the transverse maneuvering trajectory of the hypersonic aircraft
2.1 longitudinal ballistic model
According to the hypersonic aircraft motion equation in the step one, analysis shows that the maneuvering trajectory of the aircraft mainly shows oscillation characteristics in the longitudinal direction, and the oscillation amplitude of the maneuvering trajectory gradually decreases along with the time extension, namely, the maneuvering trajectory shows an attenuating oscillation form. Accordingly, the motion state vector of the aircraft in the longitudinal direction (i.e., z-axis) is modeled, and its discrete form longitudinal ballistic model is represented as:
X(k+1)＝F DO X(k)+W(k) (3)
in the formula,
Figure BDA0002952076000000092
each variable represents the position, velocity, acceleration, and jerk of the aircraft in the longitudinal direction (i.e., z-axis) at time k.
Its state transition matrix F DO Is composed of
Figure BDA0002952076000000093
In the formula,
Figure BDA0002952076000000094
Figure BDA0002952076000000101
alpha is the maximum correlation attenuation and is usually in the range of
Figure BDA0002952076000000102
Beta is the maneuvering oscillation frequency, the value range is usually between (0,1), and T is the sampling interval.
The process noise is
Figure BDA0002952076000000103
τ is the time difference, which is related to the aircraft cruise cycle.
2.2 transverse ballistic model
Analysis shows that the maneuvering trajectory of the hypersonic aerocraft mainly shows linear, C-shaped or S-shaped maneuvering in the transverse direction (namely the x axis and the y axis) and the acceleration change is relatively simple, so that the transverse trajectory of the aerocraft is modeled by a model set consisting of a uniform acceleration model and a Singer model, and the state variable of the model set is expressed as
Figure BDA0002952076000000104
Wherein each variable is represented as position, velocity, acceleration of the aircraft at time k on the lateral trajectory (i.e., x-axis and y-axis), respectively.
Modeling the transverse trajectory as a uniform acceleration model, the discretization form being expressed as:
X(k+1)＝F CA X(k)+W(k) (8)
in the formula, a state transition matrix F CA Is composed of
Figure BDA0002952076000000111
W (k) is the process noise, its covariance matrix Q CT Is composed of
Figure BDA0002952076000000112
S w Is the power spectral density of white noise.
Modeling the lateral trajectory as a Singer model, the discretized form is represented as:
X(k+1)＝F Singer X(k)+W(k) (11)
wherein,
Figure BDA0002952076000000113
the inverse of the maneuvering time constant, i.e., the maneuvering frequency, of β here is similar to the meaning of β in the longitudinal ballistic motion model. W (k) is the process noise, its covariance matrix Q Singer Is composed of
Figure BDA0002952076000000121
In the formula
Figure BDA0002952076000000122
Figure BDA0002952076000000123
Is the acceleration variance. e is approximately equal to 2.7183 and is a natural constant, and q 1-q 6 are intermediate coefficients.
Step three, obtaining radar measurement
During tracking, radar measurement data is typically established in a spherical coordinate system, while the kalman filter operates in a cartesian coordinate system. To solve the problem, the invention adopts correction unbiased measurement conversion to convert the measurement data in the spherical coordinate system into the Cartesian coordinate system, and obtains the converted measurement value of
Figure BDA0002952076000000124
Wherein A is an azimuth; e is a pitch angle; r is m Is the target distance;
Figure BDA0002952076000000125
Figure BDA0002952076000000126
in order to square the error in the pitch angle,
Figure BDA0002952076000000127
Figure BDA0002952076000000128
is the square of the azimuthal error. Measured noise covariance R of
Figure BDA0002952076000000129
Wherein R is a symmetric matrix, and each element in R after correction unbiased measurement conversion is
Figure BDA0002952076000000131
In the formula,
Figure BDA0002952076000000132
Figure BDA0002952076000000133
is the square of the range error.
Step four, designing a self-adaptive Kalman filtering method of multiple fading factors to carry out tracking filtering
4.1 Filter equation
Assuming a sampling period of T, the time state equation and the measurement equation of the discrete system can be expressed as
X(k)＝F(k-1)X(k-1)+Γ(k-1)W(k-1) (17)
Z(k)＝H(k)X(k)+V(k) (18)
Wherein, X (k) and Z (k) respectively represent a state vector and a measurement vector at the target k moment; f (k) represents a state transition matrix at time k; h (k) represents the measurement matrix at time k; Γ (k-1) is the system noise matrix at time k-1; w (k) and v (k) represent the state noise and the measurement noise at time k, respectively.
A typical discrete-time Kalman filtering process may be represented as
Figure BDA0002952076000000134
P(k+1/k)＝F(k)P(k/k)F T (k)+Γ(k)Q(k)Γ T (k) (20)
Figure BDA0002952076000000135
K(k+1)＝P(k+1/k)H T (k+1)[H(k+1)P(k+1/k)H T (k+1)+R(k+1)] -1 (22)
Figure BDA0002952076000000136
P(k+1/k+1)＝[I-K(k)H(k)]P(k+1/k) (24)
In the formula, K (K +1) is a kalman filter gain matrix, P (K +1/K) is a state prediction error covariance matrix, P (K/K) is a state estimation error covariance matrix, R (K +1) is a filter residual error, q (K) is a process noise covariance matrix, and R (K +1) is a measurement noise variance.
After the initial value of the target state vector and the initial matrix of the estimation error covariance of the state vector are given, filtering can be carried out continuously and recurrently according to the latest measurement information.
4.2 setting correction threshold
Because the measured value has errors, in order to reduce the influence of the measurement errors and ensure the excellent performance of the filter, a correction threshold is set before the predicted value is corrected. And starting the multiple fading self-adaptive filtering only when the error between the predicted value and the measured value is larger than the correction threshold, otherwise, continuing to execute the conventional Kalman filtering. Here, the error between the predicted value and the measured value is measured by r (k), the maximum error of the measured value is represented by the mean square value of the covariance of the measured error, and the threshold is set
Figure BDA0002952076000000144
Wherein χ is threshold coefficient, the value range is generally 0-1, R i,i (k) Representing the measurement noise at time k for the measurement vector at row i and column i. In the present invention, the threshold coefficient χ is set to 0.5.
4.3 selection of multiple Elimination factors
If a scalar fading factor p is introduced in the formula (20), it can be expressed as
P(k+1/k)＝ρ(k)F(k)P(k/k)F T (k)+Γ(k)Q(k)Γ T (k) (25)
In the formula, rho (k) is more than or equal to 1, the standard Kalman filtering becomes the single fading factor Kalman filtering.
The method for determining the scalar fading factor rho comprises the following steps
Figure BDA0002952076000000141
In the formula,
Figure BDA0002952076000000142
for the estimated value of covariance of innovation, tr (-) represents trace-finding operator, N (k), M (k) are intermediate parameters, and the estimated value obtained by windowing is
Figure BDA0002952076000000143
Wherein N is the sliding window length.
As can be seen from equations (25) to (27), when dealing with multi-channel filtering, the single fading factor ρ (k) has the same adjustment capability for each channel filtering. In order to avoid the defects, the invention can adjust each channel in real time according to the actual filtering condition of each channel by designing multiple fading factors rho (k).
ρ(k)＝diag{ρ 1 ,ρ 2 ,…,ρ i ,…,ρ n } (28)
Wherein n is the state vector dimension; ρ (i) represents an extinction factor corresponding to the i-th dimension state vector.
P(k+1/k)＝ρ(k)F(k)P(k/k)F T (k)+Γ(k)Q(k)Γ T (k) (29)
Define matrix j (k) ═ Γ (k-1) P (k-1) Γ T (k-1) then M (k) ═ H (k) J (k) H T (k) .1. the The measurement matrix h (k) usually has two values of 0 or 1, so that the multiple fading factors can be obtained from ρ (k) m (k) ═ n (k)The calculation method is
Figure BDA0002952076000000151
In the formula, H i (k) Representing the ith element in the measurement matrix at time k.
When ρ (k) ═ I, the multiple evanescent filtering degrades to conventional kalman filtering. From equation (30), p i (k) Value of (a) is evaluated by an innovation covariance estimate
Figure BDA0002952076000000152
The influence is large, therefore
Figure BDA0002952076000000153
The estimate directly affects p i (k) The regulating effect of (1). While
Figure BDA0002952076000000154
And is influenced by the measured value, so that the rho can be given according to the actual situation of the application scene in order to eliminate the instability brought by the measurement error to the filtering as much as possible i (k) And setting an upper value limit, which is usually a real number greater than 1, to avoid too much oscillation of filtering.
4.4 innovation sequence covariance estimation
As can be seen from equation (27), the innovation covariance estimate calculated by windowing
Figure BDA0002952076000000156
Is an arithmetic average of historical data, and cannot increase the weight value of recent measurement information. In order to increase the weight of the latest measured data and make the proportion of the old data which is elongated along with the time lower and lower, an exponential weighting innovation covariance estimation method for limiting the memory length is designed.
Setting the memory length of the covariance sequence of the defined innovation as N (N is less than or equal to k), the weighting coefficient as x (i), and satisfying the condition at the time k
Figure BDA0002952076000000155
In the formula, b is a forgetting factor, and the value range is usually 0.7-0.95. After the forgetting factor b is determined, according to the nature of the array, there are
Figure BDA0002952076000000161
Therefore, the temperature of the molten metal is controlled,
Figure BDA0002952076000000162
the recurrence formula for innovation covariance is derived below.
Innovation covariance estimation at time k
Figure BDA0002952076000000163
Is composed of
Figure BDA0002952076000000164
Innovation covariance estimation at time k-1
Figure BDA0002952076000000165
Is composed of
Figure BDA0002952076000000166
As can be seen from the formulas (33) and (34), the term I can be used
Figure BDA0002952076000000167
To obtain a recursion formula of the innovation covariance estimation value
Figure BDA0002952076000000168
The formula (35) can show that the size of the innovation covariance estimation value depends on the values of the memory length N and the forgetting factor b, wherein the selection of the values of N and b mainly depends on the intensity of the change of the system state, for a system with more gradual change, the value of the memory length N should be larger, and the value of the forgetting factor b should also be larger, so as to obtain a more accurate innovation covariance estimation value.
4.5 tracking filtering
And tracking and filtering the model predicted value obtained in the step two and the observed value obtained in the step three according to the self-adaptive filtering method.
Simulation experiment
The hypersonic aircraft target is boosted to 50km height by a transmitter and released for gliding, the reentry speed of the hypersonic aircraft is 15Ma, and the control parameter setting of the hypersonic aircraft in the gliding stage is shown in figures 8 and 9, wherein the roll angle is a constant mode, and the attack angle is a function mode related to the speed. The geographical coordinates of the observation radar are assumed to be [12 degrees, 1.5 degrees and 1km ], the sampling interval is 0.1s, the distance error is 100m, the pitch angle and azimuth angle errors are 0.1 degrees, the radar is observed in the gliding stage of the aircraft, the total observation time is 300s, no ground shielding exists in the observation process, and the distance between the target and the radar is about 300 km. The method comprises the following steps of (1) setting parameters of a motion model: the maneuvering oscillation frequency beta is 0.06rad/s, the maximum relevant attenuation quantity alpha is 1/300, the prior probability of the model on the transverse trajectory is 0.5, and the state is updated through the adaptive Kalman filtering of the multiple extinction factors. The monte carlo number of times is 100, and the corresponding position Root Mean Square Error (RMSE) is obtained as shown in fig. 4 to 6, and the speed estimation and the speed RMSE are respectively shown in fig. 7 and 8. The method for tracking the hypersonic aircraft in three dimensions can keep higher precision and stability during tracking.
Finally, it should be noted that the above detailed description is only for illustrating the technical solution of the patent and not for limiting, although the patent is described in detail with reference to the preferred embodiments, it should be understood by those skilled in the art that the technical solution of the patent can be modified or replaced by equivalents without departing from the spirit and scope of the technical solution of the patent, which should be covered by the claims of the patent.

Claims (4)

1. A three-dimensional tracking method for a hypersonic aircraft is characterized by comprising the following steps: the method comprises the following steps:
the method comprises the following steps: establishing a motion equation of the hypersonic aircraft in a glide section under a half-speed coordinate system;
step two: the maneuvering trajectory of the hypersonic aircraft is divided into a longitudinal maneuvering trajectory and a transverse maneuvering trajectory:
modeling a motion state vector of the hypersonic aircraft in the longitudinal direction, wherein a discrete form longitudinal ballistic model of the hypersonic aircraft is represented as:
X(k+1)＝F DO X(k)+W(k)
wherein X (k +1) and X (k) are state variables at time k and time k +1, respectively, and F DO A vertical state transition matrix, W (k) process noise,
Figure FDA0003653880920000011
wherein each variable respectively represents the position, the speed, the acceleration and the jerk of the aircraft on the longitudinal trajectory at the moment k;
modeling the transverse trajectory of the hypersonic aerocraft by adopting a model set consisting of a uniform acceleration model and a Singer model in the transverse direction is represented as follows:
X(k+1)＝F CA X(k)+W(k)
wherein X (k +1) and X (k) are state variables at time k and time k +1, respectively, and F CA W (k) is process noise,
Figure FDA0003653880920000012
wherein each variable respectively represents the position, the speed and the acceleration of the aircraft in the x and y directions on the transverse trajectory at the moment k;
the longitudinal state transition matrix F DO Is composed of
Figure FDA0003653880920000013
In the formula,
Figure FDA0003653880920000021
Figure FDA0003653880920000022
a is the maximum correlated attenuation and the value range is
Figure FDA0003653880920000023
Beta is the maneuvering oscillation frequency, the value range is between (0,1), and T is the sampling interval;
the transverse state transition matrix F CA Is composed of
Figure FDA0003653880920000024
Wherein T is a sampling interval;
step three: obtaining radar measurement, and converting measurement data in a spherical coordinate system into a Cartesian coordinate system by adopting correction unbiased measurement conversion to obtain a converted measurement value;
step four: and (3) performing tracking filtering on the model predicted value obtained in the step two and the measurement value obtained in the step three by adopting a self-adaptive Kalman filtering method of multiple fading factors, wherein the method specifically comprises the following steps: establishing a filtering equation, setting a correction threshold value, selecting multiple fading factors, estimating covariance of an innovation sequence and tracking filtering;
the filter equation is
X(k)＝F(k-1)X(k-1)+Γ(k-1)W(k-1)
Z(k)＝H(k)X(k)+V(k)
Wherein, X (k) and Z (k) respectively represent a state vector and a measurement vector at the target k moment; f (k) represents a state transition matrix at time k; h (k) represents the measurement matrix at time k; Γ (k-1) is the system noise matrix at time k-1; w (k) and V (k) respectively represent state noise and measurement noise at the time k;
the correction threshold value in the fourth step is set to
Figure FDA0003653880920000031
Wherein χ is a threshold coefficient, the value range is 0-1, R i,i (k) Representing the measurement noise of the ith row and the ith column of the measurement vector at the k moment;
and adjusting the multiple fading factors rho (k) selected in the fourth step in real time according to the actual filtering condition of each channel:
ρ(k)＝diag{ρ 1 ,ρ 2 ,…,ρ i ,…,ρ n }
wherein n is the state vector dimension; ρ (i) represents an extinction factor corresponding to the ith-dimension state vector; the multiple fading factors are calculated in the following way
Figure FDA0003653880920000032
In the formula, H i (k) Representing the ith element in the measurement matrix at the k moment; when rho (k) is I, the multiple evanescent filtering degenerates to conventional Kalman filtering; rho i (k) Value of (a) is evaluated by an innovation covariance estimate
Figure FDA0003653880920000033
The influence of the impact is such that,
Figure FDA0003653880920000034
the estimate directly affects p i (k) Has a regulating effect of
Figure FDA0003653880920000035
Is influenced by the measured value, and is thus given by rho i (k) Setting an upper value limit to be a real number larger than 1, and avoiding filtering oscillation;
the recurrence formula of the innovation covariance estimate is:
Figure FDA0003653880920000036
the size of the innovation covariance estimation value depends on the values of the memory length N and the forgetting factor b, wherein the selection of the values of N and b depends on the intensity of the change of the system state, and for a system with more gradual change, the value of the memory length N is larger, and the value of the forgetting factor b is also larger, so that the innovation covariance estimation value is obtained.
2. The hypersonic aircraft three-dimensional tracking method according to claim 1, characterized in that: in the first step, the equation of motion of the hypersonic aerocraft in the glide section under the half-speed coordinate system is as follows:
Figure FDA0003653880920000041
wherein r, lambda, phi, V, theta and sigma respectively represent geocentric distance, longitude, latitude, speed inclination angle and azimuth angle; a is c Representing angle of attack, g gravity acceleration, upsilon c Representing the roll angle; a is D 、a L Respectively a resistance acceleration and a lift acceleration.
3. The hypersonic aircraft three-dimensional tracking method according to claim 1, characterized in that: the measured value in the third step is
Figure FDA0003653880920000042
Wherein A is an azimuth; e is a pitch angle; r is m Is the target distance;
Figure FDA0003653880920000043
Figure FDA0003653880920000044
in order to square the error in the pitch angle,
Figure FDA0003653880920000045
Figure FDA0003653880920000046
is the square of the azimuthal error.
4. The hypersonic aircraft three-dimensional tracking method according to claim 3, characterized in that: the covariance R of the measurement noise in the measured values is
Figure FDA0003653880920000047
Wherein R is a symmetric matrix, and each element in R after correction unbiased measurement conversion is
Figure FDA0003653880920000051
In the formula,
Figure FDA0003653880920000052
Figure FDA0003653880920000053
is the square of the range error.
CN202110212919.7A 2021-02-25 2021-02-25 Three-dimensional tracking method for hypersonic aircraft Active CN113075652B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110212919.7A CN113075652B (en) 2021-02-25 2021-02-25 Three-dimensional tracking method for hypersonic aircraft

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110212919.7A CN113075652B (en) 2021-02-25 2021-02-25 Three-dimensional tracking method for hypersonic aircraft

Publications (2)

Publication Number Publication Date
CN113075652A CN113075652A (en) 2021-07-06
CN113075652B true CN113075652B (en) 2022-08-05

Family

ID=76610010

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110212919.7A Active CN113075652B (en) 2021-02-25 2021-02-25 Three-dimensional tracking method for hypersonic aircraft

Country Status (1)

Country Link
CN (1) CN113075652B (en)

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107561503B (en) * 2017-08-28 2020-08-14 哈尔滨工业大学 Adaptive target tracking filtering method based on multiple fading factors
CN109163720A (en) * 2018-08-27 2019-01-08 广西科技大学 Kalman filter tracking method based on fading memory exponent
CN110061716B (en) * 2019-01-15 2021-05-04 河海大学 Improved kalman filtering method based on least square and multiple fading factors
CN110567490B (en) * 2019-08-29 2022-02-18 桂林电子科技大学 SINS initial alignment method under large misalignment angle
CN111474960B (en) * 2020-03-23 2023-03-31 中国人民解放军空军工程大学 Critical altitude supersonic speed target tracking method based on control quantity characteristic assistance
CN111783307B (en) * 2020-07-07 2024-03-26 哈尔滨工业大学 Hypersonic aircraft state estimation method

Also Published As

Publication number Publication date
CN113075652A (en) 2021-07-06

Similar Documents

Publication Publication Date Title
CN107544067B (en) Hypersonic reentry vehicle tracking method based on Gaussian mixture approximation
JP6328789B2 (en) Method and apparatus for determining angle of arrival (AOA) in a radar warning receiver
CN105785359B (en) A kind of multiple constraint maneuvering target tracking method
CN111797478B (en) Strong maneuvering target tracking method based on variable structure multi-model
CN107192995A (en) A kind of Pure orientation underwater target tracking algorithm of multi-level information fusion
CN111274740B (en) Multi-aircraft cooperative penetration trajectory optimization design method
CN107193009A (en) A kind of many UUV cooperative systems underwater target tracking algorithms of many interaction models of fuzzy self-adaption
CN111474960B (en) Critical altitude supersonic speed target tracking method based on control quantity characteristic assistance
CN110749891B (en) Self-adaptive underwater single beacon positioning method capable of estimating unknown effective sound velocity
CN112947541B (en) Unmanned aerial vehicle intention track prediction method based on deep reinforcement learning
CN110794409A (en) Underwater single beacon positioning method capable of estimating unknown effective sound velocity
CN114740497B (en) UKF multisource fusion detection-based unmanned aerial vehicle deception method
CN107621632A (en) Adaptive filter method and system for NSHV tracking filters
CN114489101A (en) Terminal guidance control method and system for unmanned aerial vehicle
CN113075652B (en) Three-dimensional tracking method for hypersonic aircraft
CN113761662B (en) Generation method of trajectory prediction pipeline of gliding target
CN114485676B (en) Track planning method of distributed flying radar platform
CN109084772B (en) Unscented Kalman based sight line conversion rate extraction method and system
CN114565020B (en) Aircraft sensor signal fusion method based on deep belief network and extended Kalman filtering
CN114580301A (en) Hypersonic aircraft trajectory online prediction method based on parameter extrapolation
CN115686059A (en) Hypersonic aircraft flight-forbidden region avoidance guidance method based on pseudo-spectrum method
Jeong et al. A study on the performance comparison of three optimal Alpha-Beta-Gamma filters and Alpha-Beta-Gamma-Eta filter for a high dynamic target
CN112698323B (en) Full-automatic landing radar guiding noise suppression method based on alpha-beta-gamma filter
CN114966667A (en) Low-altitude maneuvering target tracking method based on interactive multiple models
CN113790727A (en) Pulse maneuver detection method based on auxiliary state parameters

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