A radar post-processing trajectory reconstruction method with outliers robustness

Abstract

The present application relates to a kind of radar post-processing trajectory reconstruction method with outlier robustness.The method includes step 1, trajectory missile kinematics equation and radar measurement equation are integrated, and mathematical model for trajectory reconstruction is established;Step 2, according to the mathematical model established in step 1, the forward filtering process of nonlinear smoothing algorithm is designed based on maximum correlation entropy criterion and iterative solution method;Step 3, based on the forward filtering process, the backward estimation process of nonlinear smoothing algorithm is designed in combination with quadratic index, and high-precision smoothing estimation of missile position and velocity in radar post-processing stage is realized.The present application can be applied to missile target tracking and after recovery, trajectory characteristic analysis, trajectory model verification and other fields, enhance the robustness of trajectory estimation for radar measurement outlier, improve the precision of trajectory reconstruction.

Description

Technical Field

[0001] The present 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 reconstructing the missile's flight trajectory using radar data measured during flight and combined with a trajectory model. Trajectory reconstruction can be divided into two types: online and ex post. Trajectory parameters obtained from online reconstruction can be used for trajectory extrapolation and missile impact prediction, while ex post reconstruction can be used to analyze trajectory characteristics and verify trajectory models. For ex post trajectory reconstruction, since radar measurement data from the entire detection phase is available, a smoothing algorithm can be used to estimate the missile's flight state. The classic Rauch-Tung-Striebel (RTS) smoother has been applied to trajectory reconstruction. However, this smoother assumes that both system process noise and measurement noise follow a Gaussian distribution, which can lead to degraded estimation performance in the presence 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-satellite multi-stage ballistic missile based on the straight line reconstruction method" (application number CN202310621656.4) is also limited to the calculation of mean and variance, and lacks analysis of high-order statistics.

[0003] Radar measurement data of missile targets is typically given as radial range, azimuth, and elevation angle in polar coordinates. The measurement noise is typically assumed to be zero-mean Gaussian white noise. However, due to factors such as missile attitude fluctuations, the radar echo position may fluctuate randomly, resulting in the radar measurement data sometimes exhibiting outliers (i.e., outliers). This noise probability distribution exhibits a thick-tailed, non-Gaussian property. While the probability of a Gaussian random variable acquiring extreme values ​​is extremely low, a random variable with a thick-tail distribution has a certain probability of acquiring extreme values. Therefore, when radar measurement data exhibit outliers, the Gaussian noise assumption is no longer valid, necessitating the development of data processing algorithms tailored to thick-tailed, non-Gaussian noise distributions. In recent years, an information measure called correlation entropy has been widely used in the estimation of non-Gaussian random systems. The correlation entropy optimization criterion can be viewed as a generalization of the minimum mean square error criterion to 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 motion state and the analysis of 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. First, a trajectory reconstruction mathematical model is established based on the missile motion model and the radar measurement equation. On this basis, a forward filtering process of a nonlinear smoothing algorithm is designed based on the 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 object, the technical solution adopted by the present invention is as follows:

[0006] A radar post-processing trajectory reconstruction method with outlier robustness includes the following steps:

[0007] Step 1: Integrate the ballistic missile kinematic equations and radar measurement equations to establish a mathematical model for trajectory reconstruction;

[0008] Step 2: Based on the mathematical model constructed 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;

[0009] 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.

[0010] The beneficial effects of the present invention are:

[0011] Aiming at the problem that radar measurements of missile targets contain outliers, which affect the accuracy of trajectory estimation, a new trajectory reconstruction method is proposed. Based on the establishment of a mathematical model for 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, respectively. This method suppresses the influence of outliers in radar measurements on the estimation of missile position and velocity, and effectively improves the accuracy of trajectory reconstruction in the radar data post-processing stage. BRIEF DESCRIPTION OF THE DRAWINGS

[0012] Figure 1 This is a flow chart of a radar post-processing trajectory reconstruction method with outlier robustness according to the present invention;

[0013] Figure 2 A simulated trajectory diagram of a radar post-processing trajectory reconstruction method with outlier robustness according to the present invention;

[0014] 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;

[0015] Figure 4This is a comparison chart of velocity estimation error results of a radar post-processing trajectory reconstruction method with outlier robustness according to the present invention. DETAILED DESCRIPTION

[0016] The present invention will be further described below with reference to the accompanying drawings and examples.

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

[0018] The specific implementation steps of the present invention are as follows:

[0019] 1. Establishment of mathematical model for trajectory reconstruction

[0020] (1) During the unpowered phase (including the free phase and the descent phase), a ballistic missile is mainly affected by the earth's gravity, atmospheric drag, centrifugal force, and Coriolis force. Its motion model can be established as follows:

[0021] ,

[0022] 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, , , 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 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 height h is expressed as , is the radius of the Earth; The ballistic coefficient is the mass-to-drag ratio, which indicates the missile's ability to overcome air resistance and maintain flight. It is approximately considered a constant. When estimating the missile's position and velocity, the ballistic coefficient also needs to be estimated. In order to ensure the positivity of the estimated ballistic coefficient, the parameter to be estimated is further introduced. , and its relationship with the ballistic coefficient is , is the initial value of the ballistic coefficient.

[0023] Defining state variables exist Time is , the superscript T represents transposition, and the kinematic equations of the ballistic missile in the unpowered phase and the parameters to be estimated are The first time derivative of After discretization, the state equation can be obtained ,in is zero-mean Gaussian white noise, and the covariance matrix is , submatrix Defined as:

[0024] ,

[0025] 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.

[0026] (2) Next, we build a radar measurement model for missiles. The longitude, latitude, and altitude of the radar station 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 radar center as the origin is expressed as , then:

[0027] ,

[0028] The radar measurement value is generally expressed in a polar coordinate system with the radar center as the origin, in the form of radial distance, azimuth and elevation angle, which can be expressed 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:

[0029] ,

[0030] Where, Represents the nonlinear function between measurement and state quantity. Adding the variable time k to the above formula, the radar measurement equation is expressed as . Measurement noise The probability distribution of has a fat-tail characteristic, which means that the measured data has a certain probability of being an outlier.

[0031] 2. Forward filtering of trajectory reconstruction algorithm

[0032] (1) The initialization of the missile position and velocity can be given by the radar's first two measurements and noise variance, 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:

[0033] ,

[0034] ,

[0035] in, and is the filtered value and corresponding variance of the state at time k-1, For function Jacobian matrix. Considering the thick-tailed non-Gaussian characteristics of the measurement noise, the following performance index function is constructed based on the correlation entropy:

[0036] ,

[0037] 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 , Represents a vector The i-th component of the intermediate matrix It can be obtained by Cholesky decomposition of the matrix, that is .

[0038] (2) The constructed indicator function By minimizing, the filtered value of the state at time k can be obtained. right The gradient of is as follows:

[0039] ,

[0040] Among them, the intermediate process matrix Defined as a diagonal matrix ;

[0041] 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:

[0042] ,

[0043] Setting the above equation to zero gives

[0044] ,

[0045] where the gain matrix is defined as ;

[0046] Note that is related to , while is related to , so the above equation is actually a fixed-point equation about the variable , which can be solved by using the fixed-point iteration method. Let denote the value of the variable at the th iteration, and similarly define , , and . The initial condition of the iteration is , and the termination condition is , , which is a given threshold parameter. Let denote the time index of the last iteration, and let denote the gain matrix after the iteration, and let denote the filter value, i.e., we have

[0047] ,

[0048] The covariance matrix of the filter error is approximately given by

[0049] ,

[0050] 3. Backward estimation of the trajectory reconstruction algorithm

[0051] Since the process noise is Gaussian distributed and the influence of the non-Gaussian measurement noise has been handled in the forward filtering process, the probability distribution , after the forward filtering can be considered as Gaussian distributed, so the backward estimation process of the trajectory reconstruction algorithm can use the following quadratic index:

[0052] ,

[0053] where the nonlinear function can be approximated as

[0054] ,

[0055] The performance index function The smoothed value of the state at time k can be obtained, and the result is as follows:

[0056] ,

[0057] The covariance matrix of the smoothed error is expressed as:

[0058] ,

[0059] At this point, the radar measurement value The smoothed estimation of the position and velocity of the missile is realized, thereby achieving the goal of high-precision trajectory reconstruction.

[0060] Figures 2 to 4 The simulation result diagram of the above-mentioned embodiment is shown in the figure. Figure 2 The trajectory generated by the simulation, the starting time and the ending time of the radar tracking the missile are the 100th second and the 600th second after the missile is launched, respectively. The radar measurement noise obeys a mixed Gaussian distribution , wherein , and the data rate is 1 Hz. The method (referred to as MCERTS) is compared with the classic extended RTS smoothing algorithm (referred to as ERTS), 100 independent Monte Carlo simulation experiments are run, and the comparison results of the position estimation error and the velocity estimation error are shown in Figure 3 and Figure 4 It can be seen that the method has smaller smoothed estimation error of the position and velocity of the missile, and effectively improves the precision of the trajectory reconstruction.

[0061] The above-mentioned specific embodiments further specifically describe the purpose, technical scheme and beneficial effects of the present application, and it should be understood that the above-mentioned is only a specific embodiment of the present application, and is not used to limit the present application, and any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present application should be included in the protection scope of the present application.

[0062] The contents not described in detail in the specification of the present application belong to the prior art known to those skilled in the art.

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: Based on the mathematical model constructed 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; Wherein, 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, They represent the velocities 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, , , 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; Record state variables exist Time is , the superscript T represents transposition, and the kinematic equations of the ballistic missile in the unpowered phase and the parameters to be estimated are The first time derivative of Discretize and get the state equation ,in, Represents the state transition process, is zero-mean Gaussian white noise, covariance matrix , submatrix Defined as: ， in, represents the sampling period, and Respectively represent the noise variance intensity 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 radar center 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 , Denote the measurement noise of radial distance, azimuth angle and elevation angle respectively, then: ， Where, Represents the nonlinear function between measurement and state quantity. Adding the variable time k to the above equation, the radar measurement equation is expressed as ; The step 2 includes: 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 is 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 index function is constructed: ， in, is the Gaussian kernel function, is the kernel bandwidth parameter; the intermediate parameter , Represents a vector The i-th component of ; the intermediate matrix By matrix The Cholesky decomposition of is obtained, and the superscript -1 indicates the inverse operation; Step 2.2: Change 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 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, Representing 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, and the gain matrix after the iteration is completed , filtered value , then: ， The covariance matrix of the filtering error is approximately expressed as: ， Where I represents the identity matrix; 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 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