Radar post-processing trajectory reconstruction method with outlier robustness

Abstract

The invention relates to a radar post-processing trajectory reconstruction method with outlier robustness. The method comprises the following steps: step 1, integrating a ballistic missile kinematics equation and a radar measurement equation, and establishing a mathematical model for trajectory reconstruction; 2, designing a forward filtering process of a nonlinear smoothing algorithm based on a maximum correlation entropy criterion and an iterative solution method according to the mathematical model established in the step 1; and step 3, based on the forward filtering process, designing a reverse estimation process of a nonlinear smoothing algorithm in combination with a quadratic index, and realizing high-precision smooth estimation of the position and the speed of the missile in a radar post-processing stage. The method can be applied to the fields of missile target tracking and post-recovery, ballistic characteristic analysis, ballistic model verification and the like, the robustness of ballistic estimation to radar measurement outliers is enhanced, and the ballistic reconstruction precision is improved.

Description

Technical Field

[0001] The invention belongs to the field of trajectory restoration, and in particular relates to a radar post-processing trajectory reconstruction method with outlier robustness. Background Art

[0002] Trajectory reconstruction refers to the process of reproducing the missile's flight trajectory using the data obtained from radar measurements during the flight of a ballistic missile combined with a trajectory model. Trajectory reconstruction can be divided into two types: online trajectory reconstruction and ex post trajectory reconstruction. The trajectory parameters obtained by online trajectory reconstruction can be used for trajectory extrapolation, missile landing point prediction, etc., while ex post trajectory reconstruction can be used to analyze trajectory characteristics, verify trajectory models, etc. For ex post trajectory reconstruction, since the measurement data of the radar during the entire detection phase has been mastered, the missile's flight state can be estimated using a smoothing algorithm. The classic Rauch-Tung-Striebel (RTS) smoother has been applied to the trajectory reconstruction problem. However, it is based on the assumption that both the system process noise and the measurement noise obey a Gaussian distribution, which may lead to a decrease in estimation performance in the case of non-Gaussian noise. In addition, the trajectory smoothing algorithm based on cubature Kalman filtering proposed in the Chinese patent application “A method for reconstructing the trajectory of a single-star multi-stage ballistic missile based on the straight line reconstruction method” (application number CN202310621656.4) is also limited to the calculation of the mean and variance, and lacks analysis of higher-order statistics.

[0003] The radar measurement data of missile targets are generally given in the form of radial distance, azimuth and pitch angle in polar coordinate system, and the measurement noise is usually assumed to be zero-mean Gaussian white noise. However, due to the change of missile attitude and other reasons, the radar echo position may oscillate randomly, and the radar measurement data sometimes appear as outliers (i.e., wild values), and the probability distribution of the noise has thick-tailed non-Gaussian characteristics. Compared with the extremely low probability of Gaussian random variables to obtain extreme values, thick-tailed distribution random variables have a certain probability of obtaining extreme values. Therefore, when the radar measurement data has wild values, the Gaussian noise assumption is no longer reasonable, and it is necessary to study data processing algorithms for thick-tailed non-Gaussian noise distribution. In recent years, an information measure called correlation entropy has been widely used in the study of estimation problems of non-Gaussian random systems. The correlation entropy optimization criterion can be regarded as a generalization of the minimum mean square error criterion in the non-Gaussian case. Therefore, considering that the existing trajectory reconstruction methods are sensitive to radar measurement outliers, which leads to reduced trajectory restoration accuracy, a high-precision trajectory reconstruction method with outlier robustness can be studied using the correlation entropy optimization criterion to improve the estimation accuracy of the missile's motion state and the analysis and understanding of the trajectory characteristics in the post-processing stage. Summary of the invention

[0004] Aiming at the problem that outliers in radar measurement data affect the accuracy of trajectory reconstruction in the post-processing stage, the present invention provides a radar post-processing trajectory reconstruction method with outlier suppression capability. Firstly, a trajectory reconstruction mathematical model is established based on a missile motion model and a radar measurement equation. On this basis, a forward filtering process of a nonlinear smoothing algorithm is designed based on a maximum correlation entropy criterion. Finally, a reverse smoothing estimation of the missile position and velocity is designed based on a quadratic index, thereby achieving the goal of high-precision trajectory reconstruction and solving the problem of high-precision trajectory restoration estimation.

[0005] To achieve the above purpose, the technical solution adopted by the present invention is as follows: A radar post-processing trajectory reconstruction method with outlier robustness includes the following steps: Step 1: Integrate the ballistic missile kinematic equations and radar measurement equations to establish a mathematical model for trajectory reconstruction; Step 2: According to the mathematical model built in step 1, a forward filtering process of a nonlinear smoothing algorithm is designed based on the maximum correlation entropy criterion and an iterative solution method; Step 3: Based on the forward filtering process, a reverse estimation process of a nonlinear smoothing algorithm is designed in combination with a quadratic index to achieve high-precision smooth estimation of the missile position and velocity in the radar post-processing stage.

[0006] The beneficial effects of the present invention are: Aiming at the problem that the radar measurement values ​​of missile targets contain wild values, thus affecting the accuracy of trajectory estimation, a new trajectory reconstruction method is proposed. On the basis of establishing a mathematical model of trajectory reconstruction, the maximum correlation entropy criterion and the quadratic index are used to design the forward filtering process and the reverse estimation process of the nonlinear smoothing algorithm, which suppresses the influence of the radar measurement wild values ​​on the missile position and velocity estimation, and effectively improves the trajectory reconstruction accuracy in the radar data post-processing stage. BRIEF DESCRIPTION OF THE DRAWINGS

[0007] Figure 1 It is a flow chart of a radar post-processing trajectory reconstruction method with outlier robustness of the present invention; Figure 2 A simulated trajectory diagram of a radar post-processing trajectory reconstruction method with outlier robustness according to the present invention; Figure 3 A comparison diagram of position estimation error results of a radar post-processing trajectory reconstruction method with outlier robustness according to the present invention; Figure 4 This is a comparison diagram of velocity estimation error results of a radar post-processing trajectory reconstruction method with outlier robustness according to the present invention. DETAILED DESCRIPTION

[0008] The present invention will be further described below in conjunction with the accompanying drawings and embodiments.

[0009] like Figure 1 As shown, the main implementation steps of the present invention include establishing a trajectory reconstruction mathematical model based on a missile motion model and a radar measurement model, designing a forward filtering process based on a related entropy index, and designing a reverse estimation process based on a quadratic index.

[0010] The specific implementation steps of the present invention are as follows: 1. Establishment of mathematical model for trajectory reconstruction (1) When a ballistic missile is in the unpowered phase (including the free phase and the descent phase), it is mainly affected by the earth's gravity, atmospheric drag, centrifugal force and Coriolis force. Its motion model can be established as follows: , in, and They represent the position and velocity of the missile along the x, y, and z axes in the Earth-centered Earth-fixed coordinate system respectively; the first two terms of the acceleration expression represent the acceleration caused by the Earth's gravity and atmospheric resistance. , , They represent the acceleration of the ballistic missile along the x, y, and z axes in the Earth-centered Earth-fixed coordinate system, respectively. , , They represent the first-order time derivatives of the position of the ballistic missile along the x, y, and z axes in the Earth-centered Earth-fixed coordinate system respectively; and The last two terms of represent the acceleration caused by Coriolis force and centrifugal force respectively; is the Earth's gravitational constant; is the angular velocity of the Earth's rotation; is the atmospheric density, is the atmospheric density at sea level, is the atmospheric scale height, and the ballistic missile altitude h is expressed as , is the radius of the Earth; is the ballistic coefficient, i.e. the mass-to-drag ratio, which indicates the missile's ability to overcome air resistance and maintain flight. It is approximately considered to be a constant. When estimating the missile's position and velocity, it is also necessary to estimate the ballistic coefficient. In order to ensure the positivity of the estimated value of the ballistic coefficient, we further introduce the parameter to be estimated. , which is related to the ballistic coefficient as follows: , is the initial value of the ballistic coefficient.

[0011] Defining state variables exist Time is The superscript T represents transposition, which transforms the kinematic equations of the ballistic missile in the unpowered phase and the parameters to be estimated into The first time derivative of After discretization, we can get the state equation ,in is zero-mean Gaussian white noise, and the covariance matrix is , the submatrix Defined as: , in, represents the sampling period, and They represent the noise variance strength corresponding to the parameter part to be estimated and the motion state part respectively.

[0012] (2) Next, we build a radar measurement model for missiles. The longitude, latitude, and altitude of the radar site are , this site is expressed in the Earth-centered Earth-fixed coordinate system as The position of the missile in the northeast celestial coordinate system with the center of the radar as the origin is expressed as , then: , The radar measurement value is generally expressed in the polar coordinate system with the radar center as the origin, in the form of radial distance, azimuth and elevation angle, denoted as , and the radar measurement noise is expressed as , Denote the measurement noise of radial distance, azimuth angle, and elevation angle respectively, then the following equation holds: , In the formula, Represents the nonlinear function between the measurement and the state quantity. Adding the variable time k to the above equation, the radar measurement equation is expressed as . Measurement noise The probability distribution of has a thick tail characteristic, which means that the measured data has a certain probability of being an outlier.

[0013] 2. Forward filtering of trajectory reconstruction algorithm (1) The initialization of the missile position and velocity can be given by the radar's measurement values ​​and noise variance at the first two moments, while the initial value of the ballistic coefficient can be selected based on experience. At time k, the state prediction value and variance As shown below: , , in, and is the filtered value and corresponding variance of the state at time k-1, For function Jacobian matrix of . Considering the thick-tailed non-Gaussian characteristics of measurement noise, the following performance index function is constructed based on the correlation entropy: , in, is the Gaussian kernel function commonly used in the definition of correlation entropy, that is, , is the kernel bandwidth parameter; the intermediate parameter Defined as , Representation vector The i-th component of the intermediate matrix It can be obtained by Cholesky decomposition of the matrix, that is .

[0014] (2) The constructed indicator function By minimizing, we can get the filtered value of the state at time k. right The gradient of is as follows: , Among them, the intermediate process matrix Defined as a diagonal matrix ; Further define the weighted covariance matrix , and Approximately expressed as , For function The Jacobian matrix of , then the above gradient vector can be expressed as: , Setting the above formula to zero yields: , The gain matrix Defined as ; Notice and about, and And So the above formula is actually about the variable The fixed point equations can be solved using the fixed point iteration method. Representation variables In the The value at the iteration is similar to the definition , , and The initial condition for the iteration is The termination condition is , is a given threshold parameter. Indicates the time stamp of the last iteration, and records the gain matrix after the iteration ends , filter value , that is: , The covariance matrix of the filtering error is approximately expressed as: , 3. Reverse estimation of trajectory reconstruction algorithm Since the process noise follows a Gaussian distribution and the influence of non-Gaussian measurement noise has been processed in the forward filtering process, it can be considered that the probability distribution obtained after forward filtering is , represents a Gaussian distribution, so the reverse estimation process of the trajectory reconstruction algorithm can use the quadratic index as shown below: , Among them, the nonlinear function It can be approximately expressed as: , By minimizing the performance index function , we can get the smoothed value of the state at time k, and the results are as follows: , The covariance matrix of the smoothed error is expressed as: , So far, using radar measurements A smooth estimation of the missile's position and velocity is achieved, thus accomplishing the goal of high-precision trajectory reconstruction.

[0015] Figures 2 to 4 Graph showing simulation results of the above implementation example. Figure 2 For the trajectory generated by the simulation, the start and end time of the radar tracking the missile are 100 seconds and 600 seconds after the missile is launched. The radar measurement noise follows a mixed Gaussian distribution ,in The data rate is 1 Hz. The proposed method (MCERTS for short) is compared with the classic extended RTS smoothing algorithm (ERTS for short). 100 independent Monte Carlo simulation experiments are run to obtain the comparison results of position estimation error and velocity estimation error. Figure 3 and Figure 4 It can be seen that the method proposed in the present invention has a smaller smoothing estimation error for the missile position and speed, and effectively improves the accuracy of trajectory reconstruction.

[0016] The specific embodiments described above further illustrate the objectives, technical solutions and beneficial effects of the present invention in detail. It should be understood that the above description is only a specific embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

[0017] The contents not described in detail in the specification of the present invention belong to the prior art known to the professional and technical personnel in this field.

Claims

1. A radar post-processing trajectory reconstruction method with outlier robustness, characterized in that: The following steps are involved: Step 1: Integrate the ballistic missile kinematic equations and radar measurement equations to establish a mathematical model for trajectory reconstruction; Step 2: According to the mathematical model built in step 1, a forward filtering process of a nonlinear smoothing algorithm is designed based on the maximum correlation entropy criterion and an iterative solution method; Step 3: Based on the forward filtering process, a reverse estimation process of a nonlinear smoothing algorithm is designed in combination with a quadratic index to achieve high-precision smooth estimation of the missile position and velocity in the radar post-processing stage.

2. The radar post-processing trajectory reconstruction method with outlier robustness according to claim 1, characterized in that: The step 1 comprises: Step 1.1: The kinematic equation of a ballistic missile during its unpowered flight is expressed as: ， in, They represent the positions of the ballistic missile along the x, y, and z axes in the Earth-centered Earth-fixed coordinate system, respectively. They represent the speed of the ballistic missile along the x, y, and z axes in the Earth-centered Earth-fixed coordinate system respectively; is the Earth's gravitational constant; is the angular velocity of the Earth's rotation; atmospheric density , is the atmospheric density at sea level, is an exponential function, is the atmospheric height, is the height of the ballistic missile; ballistic coefficient , is the initial value of the ballistic coefficient, is the parameter to be estimated; , , They represent the acceleration of the ballistic missile along the x, y, and z axes in the Earth-centered Earth-fixed coordinate system, respectively. , , They represent the first-order time derivatives of the position of the ballistic missile along the x, y, and z axes in the Earth-centered Earth-fixed coordinate system respectively; Remember state variables exist Time is The superscript T represents transposition, which transforms the kinematic equations of the ballistic missile in the unpowered phase and the parameters to be estimated into The first time derivative of Discretize and get the state equation ,in, Represents the state transition process. is zero-mean Gaussian white noise, the covariance matrix , the submatrix Defined as: ， in, represents the sampling period, and Respectively represent the noise variance strength corresponding to the parameter part to be estimated and the motion state part; Step 1.2: Record the longitude, latitude and altitude of the radar station as follows: , which is expressed in the Earth-centered Earth-fixed coordinate system as , let the position of the ballistic missile in the northeast celestial coordinate system with the radar center as the origin be expressed as , then: ， In the polar coordinate system with the center of the radar as the origin, the radar measurement values ​​include radial distance, azimuth and elevation angle, which are recorded as , and the radar measurement noise is expressed as , Respectively represent the measurement noise of radial distance, azimuth angle, and elevation angle, then: ， In the formula, Represents the nonlinear function between the measurement and the state quantity. Adding the variable time k to the above equation, the radar measurement equation is expressed as .

3. The radar post-processing trajectory reconstruction method with outlier robustness according to claim 1, characterized in that: The step 2 comprises: Step 2.1: The initialization of the position and velocity of the ballistic missile is given by the radar's measurement values ​​and noise variance at the first two moments. The initial value of the ballistic coefficient is selected based on experience. At time k, the state prediction value and variance As shown below: ， ， in, and is the filtered value and corresponding variance of the state at time k-1, function The Jacobian matrix of , based on the relevant entropy, the following performance indicator function is constructed: ， in, is the Gaussian kernel function, is the kernel bandwidth parameter; the intermediate parameter , Representation vector The i-th component of By matrix The Cholesky decomposition of is obtained, and the superscript -1 indicates the inverse operation; Step 2.2: Set the indicator function Minimize and get the filtered value of the state at time k, right The gradient of is solved as follows: ， Among them, the intermediate process matrix ; Define the weighted covariance matrix ,Will Approximately expressed as ,function The Jacobian matrix of ,but right The gradient solution of is also expressed as: ， Setting the above formula to zero yields , where the gain matrix Defined as ; Use the fixed point iteration method to solve the problem. Representation variables In the The value at the iteration, the initial condition of the iteration is The termination condition is , is a given threshold parameter, Indicates the time scale of the last iteration. After the iteration, the gain matrix is ​​recorded. , filter value , then: ， The covariance matrix of the filtering error is approximately expressed as: ， Where I represents the unit matrix.

4. The radar post-processing trajectory reconstruction method with outlier robustness according to claim 1, characterized in that: The step 3 comprises: The process noise follows a Gaussian distribution, and after forward filtering, the probability distribution is obtained , represents a Gaussian distribution, so the reverse estimation process of the trajectory reconstruction method adopts the quadratic index function as shown below: ， in, It can be expressed approximately as: ， By minimizing the quadratic index function , and obtain a smooth estimate of the state at time k: ， The covariance matrix of the smoothed error is expressed as: ， High-precision trajectory reconstruction is achieved through smooth estimation of the ballistic missile's position and velocity.

Images