CN113238218A - Near space hypersonic target tracking method based on PHD filtering - Google Patents

Near space hypersonic target tracking method based on PHD filtering Download PDF

Info

Publication number
CN113238218A
CN113238218A CN202110438668.4A CN202110438668A CN113238218A CN 113238218 A CN113238218 A CN 113238218A CN 202110438668 A CN202110438668 A CN 202110438668A CN 113238218 A CN113238218 A CN 113238218A
Authority
CN
China
Prior art keywords
target
measurement
initial
tracking method
near space
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202110438668.4A
Other languages
Chinese (zh)
Other versions
CN113238218B (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.)
Sun Yat Sen University
Original Assignee
Sun Yat Sen University
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 Sun Yat Sen University filed Critical Sun Yat Sen University
Priority to CN202110438668.4A priority Critical patent/CN113238218B/en
Publication of CN113238218A publication Critical patent/CN113238218A/en
Application granted granted Critical
Publication of CN113238218B publication Critical patent/CN113238218B/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
    • 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
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
    • G01S7/023—Interference mitigation, e.g. reducing or avoiding non-intentional interference with other HF-transmitters, base station transmitters for mobile communication or other radar systems, e.g. using electro-magnetic interference [EMI] reduction techniques
    • 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
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
    • G01S7/41—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section
    • G01S7/418—Theoretical aspects

Landscapes

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

Abstract

The invention discloses a near space hypersonic target tracking method based on PHD filtering, which comprises the following steps: initializing a filter by using measurement of an initial moment, transferring radar measurement values of the first four moments from a radar station spherical coordinate system to a radar station ENU rectangular coordinate system, and obtaining an initial target intensity function by adopting a two-point difference method; inputting the initial target into the initialized filter, and calculating according to the process of the GM-PHD filter to obtain a prediction target set; dividing the measurement into live target measurement and clutter measurement according to the Mahalanobis distance between the current time measurement and the target position predicted by the last time; obtaining the current updated value of the filter by using the measurement value at the current moment, and completing the update of the survival target; pruning and merging the Gaussian terms in the update formula of the step S4; and calculating the intensity function after pruning and merging, extracting the target state and finishing the estimation of the target state. The method realizes the accurate tracking of the unknown number of the adjacent space aircrafts under the strong clutter.

Description

Near space hypersonic target tracking method based on PHD filtering
Technical Field
The invention relates to the technical field of signal processing, in particular to a near space hypersonic target tracking method based on PHD filtering.
Background
In recent years, with the development of aerospace technology, hypersonic aircrafts based on the concept of boost-glide trajectory gradually become hot spots of domestic and foreign research. The hypersonic flight vehicle can glide and fly in the atmosphere, can rapidly maneuver in a large range, and has the advantages of long range, high speed, strong maneuverability and the like. For single-Target high-speed supersonic aircraft Tracking, a nonlinear filtering algorithm represented by Extended Kalman Filter (EKF) plays an important role in achieving good Tracking effect on high maneuvering targets, and high-precision Tracking of ballistic reentry targets is realized, such as documents [ Jawahar A, Koteswara S.Modified Polar Extended Kalman Filter (MP-EKF) for bearing-on Target Tracking [ J ]. Indian Journal of Science and Technology,2016,9(26):1-5], [ Belik B V, Belov S G.Using of Extended Kalman Filter for Mobile Target Tracking [ J ]. derived Target Tracking System radial [ J ]. program, 2017,103 ].
In a strong clutter environment, clutter measurement and target measurement are difficult to distinguish, the initial target number is unknown, and meanwhile, hypersonic multi-target tracking is necessary to be discussed aiming at the problem that a hypersonic aircraft carries a plurality of gliding warheads and the prevention capability of the hypersonic aircraft is enhanced by utilizing a bait bomb. The traditional multi-target Tracking takes Data Association as a core, a representative algorithm comprises a multi-level Hypothesis testing Method (MHT) and Joint Probability Data Association (JPDA), and when the number of targets and the number of clutters are increased, the complexity of the algorithm is exponentially increased, namely 'combined explosion' occurs. Multi-target tracking algorithms based on the Random Finite Set (RFS) theory avoid data correlation, with the most important work being the Probabilistic Hypothetical Density (PHD) filter [ Mahler R.P.S.. Multi target Bayesian filtering first-order quantities [ J ]. IEEE Transactions on Aerospace and Electronic Systems,2003,39(4): 1152-. The main implementation methods include Gaussian mixture-PHD (GM-PHD) filter and Sequential Monte Carlo-PHD (SMC-PHD) filter.
At present, the following two problems exist in the method for performing hypersonic multi-target tracking by using a PHD filter:
(1) conventional PHD filters typically assume that the initial position information of the initial object and the new object intensity function are known and, in practical applications, are not directly available.
(2) Clutter interference near target measurement often causes filter performance reduction, wrong tracking and over-estimation of target number.
Disclosure of Invention
The invention provides a near space hypersonic target tracking method based on PHD filtering, aiming at solving the problems that the traditional GMPDH algorithm needs to know the initial position information and the new target intensity function of an initial target and the performance of a filter is reduced due to clutter interference, and realizing accurate tracking on an unknown number of near space aircrafts under strong clutter.
In order to achieve the purpose of the invention, the technical scheme is as follows: a near space hypersonic target tracking method based on PHD filtering comprises the following steps:
s1: initializing a filter by using measurement of an initial moment, transferring radar measurement values of the first four moments from a radar station spherical coordinate system to a radar station ENU rectangular coordinate system, and obtaining an initial target intensity function by adopting a two-point difference method;
s2: inputting an initial target into the initialized filter, and calculating according to the process of the GM-PHD filter in the predicting step to obtain a predicted target set;
s3: dividing the measurement into live target measurement and clutter measurement according to the Mahalanobis distance between the current time measurement and the target position predicted by the last time;
s4: obtaining the current updated value of the filter by using the measurement value at the current moment, and completing the update of the survival target;
s5: setting a pruning threshold and a merging threshold to prune and merge the Gaussian items in the updating formula of the step S4;
s6: and calculating the intensity function after pruning and merging, extracting the target state and finishing the estimation of the target state.
Preferably, in step S1, the expression of the initial target intensity function is as follows:
Figure BDA0003034135100000021
wherein ,NP(x-m) a Gaussian probability density function with an argument of x, a mean of m, and a covariance of P,
Figure BDA0003034135100000022
is the initial weight of the jth gaussian term,
Figure BDA0003034135100000023
is the initial covariance of the jth gaussian term.
Figure BDA0003034135100000031
Representing an initial target state, the expression for which is as follows:
Figure BDA0003034135100000032
in the formula ,x1j(k)，x2j(k)，x3j(k) Measuring coordinates in x, y and z directions under an ENU rectangular coordinate system of the radar station for the jth measurement at the moment k;
Figure BDA0003034135100000033
Figure BDA0003034135100000034
is an inverse function of the measurement equation in the x, y, z dimensions,
Figure BDA0003034135100000035
j-th measurement at time k;
Figure BDA0003034135100000036
Figure BDA0003034135100000037
Figure BDA0003034135100000038
set all initial Gaussian item weights
Figure BDA0003034135100000039
Is 1.
Further, the initial target is input into the initialized filter for prediction, wherein the equation of the state transition of the initial target is expressed as follows:
X(k+1)＝F(k)X(k)+W(k)
wherein F (k) is a state transition matrix having
Figure BDA00030341351000000310
wherein ,FijIs a zero matrix, and i ≠ j,
Figure BDA00030341351000000311
in the formula ：
p1＝(2-2αT+α2T2-2e-αT)/(2α3)
q1＝(αT-1+e-αT)/α2
r1＝(1-e-αT)/α
s1＝e-αT
w (k) represents a Gaussian white noise sequence with a mean of 0 and a covariance of Q (k);
the specific calculation formula of the prediction step is as follows:
vk|k-1(xk|Z1:k-1)＝vs,k|k-1(xk|Z1:k-1)+yk(xk)
wherein
Figure BDA0003034135100000041
Figure BDA0003034135100000042
Figure BDA0003034135100000043
wherein ,yk(xk) As a function of the nascent object intensity.
Still further, in step S3, the dividing measure is specifically as follows:
starting at time 2, the measurement at time k is divided into:
Figure BDA0003034135100000044
wherein ,skIs the number of Gauss terms at time k, ZkCollecting all the measurements;
Figure BDA0003034135100000045
Figure BDA0003034135100000046
wherein ,RkTo measure the noise covariance;
the measurement sets falling within the threshold are taken as the surviving target measurements,
Figure BDA0003034135100000047
the measurement sets falling outside the threshold are used as clutter measurements,
Figure BDA0003034135100000048
where T represents a threshold value.
Still further, the threshold value is determined by the following equation if PGTo correctly measure the probability of falling within the validation region, there are
T＝-2ln(1-PG)。
Still further, in step S4, a measurement set is used
Figure BDA0003034135100000049
And updating the survival target according to the following formula:
Figure BDA00030341351000000410
wherein ,
Figure BDA0003034135100000051
Figure BDA0003034135100000052
Figure BDA0003034135100000053
Figure BDA0003034135100000054
Figure BDA0003034135100000055
Figure BDA0003034135100000056
wherein ,HkA Jacobian matrix at the time k for the measurement equation h (·); kk(zk) Representing clutter intensity, pDIndicates the target detection probability, RkMeasuring a noise covariance matrix;
Figure BDA0003034135100000057
is a Gaussian function; hk TIs HkFurther, the clutter intensity K needs to be recalculatedk(zk)，
Let V be the whole observation region area, λ be the clutter average, VkFor the new observation region, V is set without considering the overlapping conditionkIs taken as the sum of the areas of the threshold regions corresponding to all measurements, i.e.
Figure BDA0003034135100000058
Figure BDA0003034135100000059
Intensity of the clutter
Figure BDA00030341351000000510
Still further, in step S5, the pruning and merging specifically include the following steps:
let l be 0 and make l be 0,
Figure BDA00030341351000000511
t _ prun is a pruning threshold;
when set I is not empty, the following process is repeated:
l＝l+1，
Figure BDA00030341351000000512
Figure BDA00030341351000000513
Figure BDA00030341351000000514
Figure BDA00030341351000000515
Figure BDA00030341351000000516
i ═ I \ L, up to
Figure BDA0003034135100000061
Wherein U _ merg denotes a combining threshold.
Still further, in step S5, the combined clipping strength function is:
Figure BDA0003034135100000062
still further, in step S5, the target state is extracted according to the following formula:
Figure BDA0003034135100000063
the invention has the following beneficial effects:
the method can realize the tracking of a plurality of high maneuvering targets, particularly has good performance under the condition of high nonlinearity of target flight tracks, can overcome the defect that the traditional GMPDH algorithm needs initial target information, and can track the unknown number of near space vehicles in a strong clutter environment.
Drawings
Fig. 1 is a flowchart illustrating steps of a method for tracking a hypersonic target in a near space according to the embodiment.
Fig. 2 is a schematic view of a boost-glide trajectory according to the present embodiment.
Fig. 3 is a diagram illustrating a tracking result by the method according to the present embodiment.
Fig. 4 is a diagram showing the effect of the estimation of the target number of the first 100 seconds by the method according to the present embodiment, in which a represents the conventional EK-GMPHD method and B represents the target tracking method according to the present embodiment.
Fig. 5 is a diagram of the effect of evaluating the target state estimation using the Wasserstein distance.
Detailed Description
The invention is described in detail below with reference to the drawings and the detailed description.
Example 1
As shown in fig. 1, a method for tracking a hypersonic target in a close space based on PHD filtering includes the following steps:
step S1: initializing a filter by using measurement of an initial moment, transferring radar measurement values of the first four moments from a radar station spherical coordinate system to a radar station ENU rectangular coordinate system, and obtaining an initial target intensity function by adopting a two-point difference method;
the expression of the initial target intensity function is as follows:
Figure BDA0003034135100000064
wherein ,NP(x-m) a gaussian probability density function with an argument of x, a mean of m, and a covariance of P;
Figure BDA0003034135100000071
the initial weight of the jth Gaussian item;
Figure BDA0003034135100000072
is the initial covariance of the jth gaussian term.
Figure BDA0003034135100000073
Representing an initial target state, the expression for which is as follows:
Figure BDA0003034135100000074
in the formula ,x1j(k)，x2j(k)，x3j(k) Measuring coordinates in x, y and z directions under an ENU rectangular coordinate system of the radar station for the jth measurement at the moment k;
Figure BDA0003034135100000075
Figure BDA0003034135100000076
is an inverse function of the measurement equation in the x, y, z dimensions,
Figure BDA0003034135100000077
j-th measurement at time k;
Figure BDA0003034135100000078
Figure BDA0003034135100000079
Figure BDA00030341351000000710
set all initial Gaussian item weights
Figure BDA00030341351000000711
Is 1.
Step S2: inputting an initial target into the initialized filter, and calculating according to the process of the GM-PHD filter in the predicting step to obtain a predicted target set;
wherein the equation for the state transition of the initial target is expressed as follows:
X(k+1)＝F(k)X(k)+W(k)
wherein F (k) is a state transition matrix having
Figure BDA00030341351000000712
wherein ,FijIs a zero matrix, and i ≠ j,
Figure BDA00030341351000000713
in the formula ：
p1＝(2-2αT+α2T2-2e-αT)/(2α3)
q1＝(αT-1+e-αT)/α2
r1＝(1-e-αT)/α
s1＝e-αT
w (k) represents a Gaussian white noise sequence with a mean of 0 and a covariance of Q (k).
The specific calculation formula of the prediction step is as follows:
vk|k-1(xk|Z1:k-1)＝vs,k|k-1(xk|Z1:k-1)+yk(xk)
wherein
Figure BDA0003034135100000081
Figure BDA0003034135100000082
Figure BDA0003034135100000083
wherein ,yk(xk) As a function of the nascent object intensity.
Step S3: dividing the measurement into live target measurement and clutter measurement according to the Mahalanobis distance between the current time measurement and the target position predicted by the last time;
starting at time 2, the measurement at time k is divided into:
Figure BDA0003034135100000084
wherein ,skIs the number of Gauss terms at time k, ZkCollecting all the measurements;
Figure BDA0003034135100000085
Figure BDA0003034135100000086
wherein ,RkTo measure the noise covariance;
the measurement sets falling within the threshold are taken as the surviving target measurements,
Figure BDA0003034135100000087
the measurement sets falling outside the threshold are used as clutter measurements,
Figure BDA0003034135100000088
where T represents a threshold value.
The threshold value is determined by the following formula if PGTo correctly measure the probability of falling within the validation region, there are
T＝-2ln(1-PG)。
S4: obtaining the current updated value of the filter by using the measurement value at the current moment, and completing the update of the survival target;
usage measurement collections
Figure BDA0003034135100000091
And updating the survival target according to the following formula:
Figure BDA0003034135100000092
wherein ,
Figure BDA0003034135100000093
Figure BDA0003034135100000094
Figure BDA0003034135100000095
Figure BDA0003034135100000096
Figure BDA0003034135100000097
Figure BDA0003034135100000098
wherein ,HkA Jacobian matrix at the time k for the measurement equation h (·); kk(zk) Representing clutter intensity, pDIndicates the target detection probability, RkTo measure the noise covariance matrix, Hk TIs HkThe transposing of (1).
Clutter intensity K in the above equation setk(zk) Recalculation is required, as follows:
let V be the whole observation region area, λ be the clutter average, VkFor the new observation region, V is set without considering the overlapping conditionkIs taken as the sum of the areas of the threshold regions corresponding to all measurements, i.e.
Figure BDA0003034135100000099
Figure BDA00030341351000000910
Intensity of the clutter
Figure BDA00030341351000000911
Step S5: setting a pruning threshold and a merging threshold to prune and merge the Gaussian items in the updating formula of the step S4;
let l be 0 and make l be 0,
Figure BDA00030341351000000912
t _ prun is a pruning threshold;
when set I is not empty, the following process is repeated:
l＝l+1，
Figure BDA00030341351000000913
Figure BDA00030341351000000914
Figure BDA0003034135100000101
Figure BDA0003034135100000102
Figure BDA0003034135100000103
i ═ I \ L, up to
Figure BDA0003034135100000104
Wherein U _ merg denotes a combining threshold.
S6: and calculating the intensity function after pruning and merging, extracting the target state and finishing the estimation of the target state.
The intensity function after pruning and merging is:
Figure BDA0003034135100000105
extracting the target state according to:
Figure BDA0003034135100000106
in order to further improve the technical effect of the method for tracking a hypersonic target in a near space according to this embodiment, the following specific examples are given:
tracking two high maneuvering targets under an ENU three-dimensional rectangular coordinate system of a radar station, wherein a scanning period T is 1s, and the initial states of the targets are as shown in a table 1:
TABLE 1 target initial State
Figure BDA0003034135100000107
The aerodynamic model of the near space hypersonic aerocraft is a set of nonlinear ordinary differential equations,
Figure BDA0003034135100000108
Figure BDA0003034135100000109
Figure BDA00030341351000001010
wherein u is the earth gravity constant, and is 3.986005 × 1014m3/s2，ReWhich is the radius of the earth, is,
Figure BDA0003034135100000111
adas a resistance parameter, atAs a bending force parameter, acIs a climbing force parameter. ρ is the air density, an exponential function of the target height h,
ρ(h)＝ρ0e-κh
wherein ρ0＝1.225 0kg/m3Sea level atmospheric density, k 1.186 × 10-4。
The numerical solution of the differential equation set is solved by using a fourth-order Rungku tower method to obtain a simulated trajectory, the initial value of the differential equation set is shown in Table 1, and the boosting-gliding trajectory is shown in figure 2.
Setting the maneuvering frequency alpha to be 0.5, and setting the jerk variance of the target in the x, y and z directions
Figure BDA0003034135100000112
Figure BDA0003034135100000112
Figure BDA0003034135100000112
7, 2 and 9 respectively, and the standard deviation of the measured noise is respectively
Figure BDA0003034135100000113
Probability of target survival ps0.99, target detection probability pDProbability P of measurement falling into the confirmation area of 0.98G0.99, the gaussian pruning threshold T _ prun is 10-5The combining threshold U _ merg is 4, and the clutter rate is 2. The specific simulation experiment steps are as follows:
(1) initializing variables according to the simulation parameter settings described above, and initializing a filter according to step S1 of the method described in this embodiment;
(2) obtaining a prediction target set according to step S2 of the method of this embodiment;
(3) dividing the obtained measurements according to step S3 of the method of this embodiment;
(4) step S4 of the method according to this embodiment is to obtain the current updated value of the filter by using the measured value at the current time, where the measurement equation h (-) is
Figure BDA0003034135100000114
wherein ,ZK＝[R A E]TAnd (4) measuring the value at the moment k, wherein R is the distance from the target to the radar station, A is an azimuth angle, and E is a pitch angle. Measurement noise vk R，vk A，vk EIndependent of each other, obey a zero-mean gaussian distribution. Jacobian matrix H of measurement equation H (-)kIs composed of
Figure BDA0003034135100000115
wherein ,
Figure BDA0003034135100000121
(5) pruning and merging the gaussian terms according to step S5 of the method described in this embodiment;
(6) extracting a target state according to step S6 of the method of the present embodiment;
(7) steps S2-S6 loop until finished.
The final tracking result is shown in fig. 3, and the final target number estimation of the first 100 seconds is shown in fig. 4, wherein a represents the conventional EK-GMPHD method, and B represents the target tracking method according to the present embodiment.
It can be seen that the conventional EK-gmpld method may cause the target number to be overestimated, and the time average of the absolute error of the target estimation is used as the evaluation index of the target number estimation, that is:
Figure BDA0003034135100000122
wherein ,
Figure BDA0003034135100000123
is an estimate of the k-th second target number, XkThe target number estimation error of the conventional EK-gmpld method is 0.1875 for the k-th second real target number, and the target number estimation error of the target tracking method described in this embodiment is 0.0643, which is only 34.27% of that of the conventional method. Therefore, the target tracking method has the advantages that the effect is obviously better than that of the traditional method, and the target number estimation is more accurate.
The Wasserstein distance is an index for measuring the similarity degree of two sets of point sets, the smaller the value of the Wasserstein distance is, the higher the target estimation precision is, and the experimental result is shown in the attached figure 5, wherein A represents the conventional EK-GMPDH method, and B represents the target tracking method described in this embodiment.
It can be seen that the target state estimation accuracy of the target tracking method described in this embodiment is better than that of the conventional method. The time average of the Wasserstein distance of the conventional EK-GMPLD method is 37704, and the time average of the Wasserstein distance of the object tracking method described in this embodiment is 19340, which is only 51.29% of that of the conventional EK-GMPLD method.
It should be understood that the above-described embodiments of the present invention are merely examples for clearly illustrating the present invention, and are not intended to limit the embodiments of the present invention. Any modification, equivalent replacement, and improvement made within the spirit and principle of the present invention should be included in the protection scope of the claims of the present invention.

Claims (10)

1. A near space hypersonic target tracking method based on PHD filtering is characterized in that: the method comprises the following steps:
s1: initializing a filter by using measurement of an initial moment, transferring radar measurement values of the first four moments from a radar station spherical coordinate system to a radar station ENU rectangular coordinate system, and obtaining an initial target intensity function by adopting a two-point difference method;
s2: inputting an initial target into the initialized filter, and calculating according to the process of the GM-PHD filter in the predicting step to obtain a predicted target set;
s3: dividing the measurement into live target measurement and clutter measurement according to the Mahalanobis distance between the current time measurement and the target position predicted by the last time;
s4: obtaining the current updated value of the filter by using the measurement value at the current moment, and completing the update of the survival target;
s5: setting a pruning threshold and a merging threshold to prune and merge the Gaussian items in the updating formula of the step S4;
s6: and calculating the intensity function after pruning and merging, extracting the target state and finishing the estimation of the target state.
2. The PHD filtering-based near space hypersonic target tracking method according to claim 1, characterized in that: in step S1, the expression of the initial target intensity function is as follows:
Figure FDA0003034135090000011
wherein ,NP(x-m) a Gaussian probability density function with an argument of x, a mean of m, and a covariance of P,
Figure FDA0003034135090000012
the initial weight of the jth Gaussian item;
Figure FDA0003034135090000013
representing an initial target state, the expression for which is as follows:
Figure FDA0003034135090000021
in the formula ,x1j(k)，x2j(k)，x3j(k) Measuring coordinates in x, y and z directions under an ENU rectangular coordinate system of the radar station for the jth measurement at the moment k;
Figure FDA0003034135090000022
Figure FDA0003034135090000023
is an inverse function of the measurement equation in the x, y, z dimensions,
Figure FDA0003034135090000024
j-th measurement at time k;
Figure FDA0003034135090000025
Figure FDA0003034135090000026
Figure FDA0003034135090000027
set all initial Gaussian item weights
Figure FDA0003034135090000028
Is 1.
3. The PHD filtering-based near space hypersonic target tracking method according to claim 2, characterized in that: inputting an initial target into the initialized filter for prediction, wherein the equation of the state transition of the initial target is expressed as follows:
X(k+1)＝F(k)X(k)+W(k)
wherein F (k) is a state transition matrix having
Figure FDA0003034135090000029
wherein ,FijIs a zero matrix, and i ≠ j,
Figure FDA00030341350900000210
in the formula ：
p1＝(2-2αT+α2T2-2e-αT)/(2α3)
q1＝(αT-1+e-αT)/α2
r1＝(1-e-αT)/α
s1＝e-αT
w (k) represents a Gaussian white noise sequence with a mean of 0 and a covariance of Q (k);
the specific calculation formula is as follows:
vk|k-1(xk|Z1:k-1)＝vs,k|k-1(xk|Z1:k-1)+yk(xk)
wherein
Figure FDA0003034135090000031
Figure FDA0003034135090000032
Figure FDA0003034135090000033
wherein ,yk(xk) As a function of the nascent object intensity.
4. The PHD filtering-based near space hypersonic target tracking method according to claim 3, characterized in that: in step S3, the division measurement is as follows:
starting at time 2, the measurement at time k is divided into:
Figure FDA0003034135090000034
wherein ,skIs the number of Gauss terms at time k, ZkCollecting all the measurements;
Figure FDA0003034135090000035
Figure FDA0003034135090000036
wherein ,RkTo measure the noise covariance;
the measurement sets falling within the threshold are taken as the surviving target measurements,
Figure FDA0003034135090000037
the measurement sets falling outside the threshold are used as clutter measurements,
Figure FDA0003034135090000038
where T represents a threshold value.
5. The PHD filtering-based near space hypersonic target tracking method according to claim 4, characterized in that: the threshold value is determined by the following formula if PGTo correctly measure the probability of falling within the validation region, there are
T＝-2ln(1-PG)。
6. The PHD filtering-based near space hypersonic target tracking method according to claim 4, characterized in that: step S4, using the measurement set
Figure FDA0003034135090000039
And updating the survival target according to the following formula:
Figure FDA0003034135090000041
wherein ,
Figure FDA0003034135090000042
Figure FDA0003034135090000043
Figure FDA0003034135090000044
Figure FDA0003034135090000045
Figure FDA0003034135090000046
Figure FDA0003034135090000047
wherein ,HkA Jacobian matrix at the time k for the measurement equation h (·); kk(zk) Representing clutter intensity; p is a radical ofDIndicates the target detection probability, RkMeasuring a noise covariance matrix;
Figure FDA0003034135090000048
is a Gaussian function; hk TIs HkThe transposing of (1).
7. The PHD filtering-based near space hypersonic target tracking method according to claim 6, characterized in that: need to recalculate clutter intensity Kk(zk)，
Let V be the whole observation region area, λ be the clutter average, VkFor the new observation region, V is set without considering the overlapping conditionkIs taken as the sum of the areas of the threshold regions corresponding to all measurements, i.e.
Figure FDA0003034135090000049
Figure FDA00030341350900000410
Intensity of the clutter
Figure FDA00030341350900000411
8. The PHD filtering-based near space hypersonic target tracking method according to claim 7, characterized in that: step S5, the pruning and merging are specifically as follows:
let l be 0 and make l be 0,
Figure FDA00030341350900000412
t _ prun is a pruning threshold;
when set I is not empty, the following process is repeated:
l＝l+1，
Figure FDA00030341350900000413
Figure FDA00030341350900000414
Figure FDA0003034135090000051
Figure FDA0003034135090000052
Figure FDA0003034135090000053
i ═ I \ L, up to
Figure FDA0003034135090000054
Wherein U _ merg denotes a combining threshold.
9. The PHD filtering-based near space hypersonic target tracking method according to claim 8, characterized in that: in step S5, the combined pruning strength function is:
Figure FDA0003034135090000055
10. the PHD filtering-based near space hypersonic target tracking method according to claim 8, characterized in that: step S5, extracting a target state according to the following formula:
Figure FDA0003034135090000056
CN202110438668.4A 2021-04-22 2021-04-22 Near space hypersonic speed target tracking method based on PHD filtering Active CN113238218B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110438668.4A CN113238218B (en) 2021-04-22 2021-04-22 Near space hypersonic speed target tracking method based on PHD filtering

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110438668.4A CN113238218B (en) 2021-04-22 2021-04-22 Near space hypersonic speed target tracking method based on PHD filtering

Publications (2)

Publication Number Publication Date
CN113238218A true CN113238218A (en) 2021-08-10
CN113238218B CN113238218B (en) 2023-10-20

Family

ID=77128964

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110438668.4A Active CN113238218B (en) 2021-04-22 2021-04-22 Near space hypersonic speed target tracking method based on PHD filtering

Country Status (1)

Country Link
CN (1) CN113238218B (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116909303A (en) * 2023-07-14 2023-10-20 中国人民解放军国防科技大学 Process noise self-adaptive adjusting method for near space target tracking

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2014169865A (en) * 2013-03-01 2014-09-18 Hitachi Ltd Target tracking device, target tracking program and target tracking method
US20140372073A1 (en) * 2011-05-04 2014-12-18 Jacques Georgy Two-stage filtering based method for multiple target tracking
CN106407677A (en) * 2016-09-09 2017-02-15 南京理工大学 Multi-target tracking method in case of loss of measurement data
CN107831490A (en) * 2017-12-01 2018-03-23 南京理工大学 A kind of improved more extension method for tracking target
CN110320512A (en) * 2019-07-09 2019-10-11 大连海事大学 A kind of GM-PHD smothing filtering multi-object tracking method based on tape label
CN111722214A (en) * 2020-06-03 2020-09-29 昆明理工大学 Radar multi-target tracking PHD implementation method
CN111856442A (en) * 2020-07-03 2020-10-30 哈尔滨工程大学 Multi-target tracking method for self-adaptively estimating strength of newborn target based on measured value driving
CN112328959A (en) * 2020-10-14 2021-02-05 哈尔滨工程大学 Multi-target tracking method based on adaptive extended Kalman probability hypothesis density filter

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20140372073A1 (en) * 2011-05-04 2014-12-18 Jacques Georgy Two-stage filtering based method for multiple target tracking
JP2014169865A (en) * 2013-03-01 2014-09-18 Hitachi Ltd Target tracking device, target tracking program and target tracking method
CN106407677A (en) * 2016-09-09 2017-02-15 南京理工大学 Multi-target tracking method in case of loss of measurement data
CN107831490A (en) * 2017-12-01 2018-03-23 南京理工大学 A kind of improved more extension method for tracking target
CN110320512A (en) * 2019-07-09 2019-10-11 大连海事大学 A kind of GM-PHD smothing filtering multi-object tracking method based on tape label
CN111722214A (en) * 2020-06-03 2020-09-29 昆明理工大学 Radar multi-target tracking PHD implementation method
CN111856442A (en) * 2020-07-03 2020-10-30 哈尔滨工程大学 Multi-target tracking method for self-adaptively estimating strength of newborn target based on measured value driving
CN112328959A (en) * 2020-10-14 2021-02-05 哈尔滨工程大学 Multi-target tracking method based on adaptive extended Kalman probability hypothesis density filter

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
QIAN ZHANG等: "Improved Bearings-Only Multi-Target Tracking with GM-PHD Filtering", 《SENSORS》, vol. 16, no. 9, pages 1 - 18 *
王鲁平等: "Kernel density estimation and marginalized-particle based probability hypothesis density filter for multi-target tracking", 《JOURNAL OF CENTRAL SOUTH UNIVERSITY》, vol. 22, no. 3, pages 956 - 965, XP035469277, DOI: 10.1007/s11771-015-2606-7 *
袁常顺等: "基于PHD滤波的相控阵雷达多目标跟踪算法", 《系统工程与电子技术》, vol. 38, no. 3, pages 539 - 544 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116909303A (en) * 2023-07-14 2023-10-20 中国人民解放军国防科技大学 Process noise self-adaptive adjusting method for near space target tracking
CN116909303B (en) * 2023-07-14 2024-02-02 中国人民解放军国防科技大学 Process noise self-adaptive adjusting method for near space target tracking

Also Published As

Publication number Publication date
CN113238218B (en) 2023-10-20

Similar Documents

Publication Publication Date Title
EP1610152B1 (en) Tracking of a moving object for a self-defence system
CN112613532B (en) Moving target tracking method based on radar and cyclic neural network complement infrared fusion
CN109791200B (en) Wire and tower classification based on trajectory tracking
CN114859339A (en) Multi-target tracking method based on millimeter wave radar
CN106990447B (en) A kind of multiple mobile object body monitoring method based on gravitational vectors and its gradient tensor
CN104199022A (en) Target modal estimation based near-space hypersonic velocity target tracking method
Agate et al. Road-constrained target tracking and identification using a particle filter
CN111830501B (en) HRRP history feature assisted signal fuzzy data association method and system
CN115204212A (en) Multi-target tracking method based on STM-PMBM filtering algorithm
CN113702940B (en) Spatial cluster target resolution method based on multi-element characteristic information hierarchical fusion and application
CN113030940B (en) Multi-star convex type extended target tracking method under turning maneuver
CN113238218A (en) Near space hypersonic target tracking method based on PHD filtering
CN111931287B (en) Near space hypersonic target trajectory prediction method
CN116047495B (en) State transformation fusion filtering tracking method for three-coordinate radar
CN105373805A (en) A multi-sensor maneuvering target tracking method based on the principle of maximum entropy
CN116224320B (en) Radar target tracking method for processing Doppler measurement under polar coordinate system
CN112986978A (en) Method for obtaining trust degree of radar target tracking filtering
CN111896946A (en) Continuous time target tracking method based on track fitting
CN116500602A (en) Multi-target tracking track management method based on passive distributed radar system
Zhang et al. Improved interacting multiple model-new nearest neighbor data association algorithm
CN114565020A (en) Aircraft sensor signal fusion method based on deep belief network and extended Kalman filtering
Lu et al. A new performance index for measuring the effect of single target tracking with Kalman particle filter
RU2726189C1 (en) Device for recognition of targets, which are not objects of reconnaissance
Zhang et al. EMD-based gray association algorithm for group ballistic target
Hardiman et al. Nonlinear estimation techniques for impact point prediction of ballistic targets

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