Nonlinear filtering trajectory drop point prediction method based on Runge-Kutta algorithm

Abstract

The invention discloses a nonlinear filtering trajectory drop point prediction method based on a Runge-Kutta algorithm, and the method comprises the steps: firstly, preprocessing the trace point information of a shell, and converting a radar array plane rectangular coordinate system into an east-north-sky coordinate system; then, constructing a projectile mass point trajectory model, and if a reliable track exists, carrying out point track association on the preprocessed point track data and the reliable track based on an extended Kalman filtering algorithm; using a logic method track initiation algorithm to associate the trace point data which are not associated with the temporary track which is not prohibited so as to realize track initiation; carrying out track extrapolation on unupdated reliable tracks; performing flight path extinction on the reliable flight path meeting the extinction condition; and finally, resolving a trajectory by utilizing a Runge-Kutta method and combining a trajectory equation set, and further calculating a trajectory drop point. According to the method, target tracking filtering is accurate, trajectory extrapolation errors are reduced, the number of required measurement parameters is small, the calculation speed is high, rapid drop point prediction can be achieved, and the real-time requirement of actual engineering is met.

Description

Technical Field

[0001] The invention belongs to the field of trajectory landing point prediction, and in particular relates to a nonlinear filtering trajectory landing point prediction method based on the Runge-Kutta algorithm. Background Art

[0002] In range tests, the traditional manual search method is time-consuming, labor-intensive and inefficient because the projectile's landing point has a wide range and is difficult to accurately determine, and the recovery process is complicated. Effective landing point prediction can not only improve the hit rate, but also help optimize the deployment of weapon platforms, correct launch conditions, and reduce casualties.

[0003] There are generally two types of existing artillery point prediction technologies based on observation data. One type is a machine learning method based on prior training samples. Because the number of artillery launch tests is limited, it is difficult to provide enough accurate prior training samples for machine learning, and the artillery trajectory identification is required to have strong real-time performance in the process of ballistic point prediction, so this type of method is not currently suitable for point prediction. The other type is a Kalman filter estimation algorithm based on a theoretical ballistic model. Due to the nonlinear characteristics of the radar measurement equation and the ballistic state equation, although the tracking accuracy can be improved, due to the influence of meteorological factors and system errors, the radar measurement data often contains noise and discontinuity, which makes it difficult for the traditional Kalman filter method to effectively filter out the noise, and may have problems of insufficient accuracy or divergence under the nonlinear ballistic model. The original Kalman filter method is no longer applicable, so it is technically necessary to establish a random nonlinear filter. Summary of the invention

[0004] The purpose of the present invention is to provide a nonlinear filtering trajectory impact point prediction method based on the Runge-Kutta algorithm, based on the numerical analysis method of the fourth-order Runge-Kutta algorithm, combined with radar measurement data, particle trajectory model and meteorological parameters, using the extended Kalman algorithm for tracking filtering, providing accurate state estimation and dynamic adjustment in a highly nonlinear environment, and realizing real-time tracking and impact point prediction of the projectile trajectory.

[0005] The technical solution to achieve the purpose of the present invention is: a nonlinear filtering trajectory drop point prediction method based on the Runge-Kutta algorithm, comprising the following steps:

[0006] Step 1, pre-processing the point information of the shell;

[0007] Step 2, constructing a projectile particle trajectory model;

[0008] Step 3: If there is no reliable track at present, go directly to step 4; if there is a reliable track, perform point-to-point association between the preprocessed point track data and the reliable track based on the extended Kalman filter algorithm;

[0009] Step 4, using the logical method track initiation algorithm, the point track data that has not been associated is associated with the temporary track that has not yet set sail, so as to realize the track initiation;

[0010] Step 5: Extrapolate the track of the unupdated reliable track;

[0011] Step 6, the reliable tracks that meet the extinction conditions are extincted;

[0012] Step 7, repeat steps 1 to 6, and track the projectile in real time until the target reaches the radar blind spot. When the point trace cannot be obtained, the fourth-order Runge-Kutta method is used in combination with the ballistic equations to solve the trajectory and then calculate the trajectory landing point.

[0013] Furthermore, the point information of the shells is preprocessed; the rotation matrix of the radar array rectangular coordinate system is converted to the northeast sky coordinate system:

[0014]

[0015] where ε A , ε E are the array azimuth and array elevation, R Z (ε A ), R X (ε E ) represent the rotation around the Z axis ε A The rotation matrix of angle, rotate ε around the X axis E The rotation matrix of the angle.

[0016] Furthermore, the projectile particle trajectory model is constructed as follows:

[0017]

[0018] in, is the cross-sectional area of ​​the projectile, d is the projectile diameter; m is the mass of the projectile; ρ is the air density; c x is the drag coefficient; w x , w y are longitudinal wind and cross wind respectively, g is the acceleration due to gravity; v r is the relative velocity between the projectile and the air, v x , v y , v z is the component of the projectile's center-of-mass velocity in the northeast celestial coordinate system.

[0019] Furthermore, when the point track data and the reliable track are associated, according to the ellipsoid gate rule, if

[0020]

[0021] The converted measurement value Z c(k+1) is the candidate echo; is the measurement prediction value, S(k+1) is the measurement prediction covariance, and the parameter γ is given by χ 2 The distribution table is obtained.

[0022] Furthermore, when the point-to-point data and the reliable track are associated based on the extended Kalman filter algorithm, the weighted norm of the new information is adjusted according to the nearest neighbor criterion.

[0023]

[0024] Achieve the minimum conversion measurement value Z c (k+1) will be used in the filter to update the target state.

[0025] Furthermore, when the point-to-point data and the reliable track are associated based on the extended Kalman filter algorithm, the reliable track is updated according to the state update equation and the covariance update equation; the state equation is

[0026]

[0027] in is the one-step prediction of the state, K(k+1) is the Kalman gain matrix, is the measured predicted value, Z(k+1) is the measured value;

[0028] The covariance update equation is

[0029]

[0030] Where I is the identity matrix of the same dimension as the covariance, P(k+1|k) is the one-step prediction of the covariance, K′(k+1) is the transpose of the Kalman gain matrix, and h x (k) is the Jacobian matrix of the measurement equation, and R(k+1) is the measurement noise covariance.

[0031] Furthermore, the fourth-order Runge-Kutta method is used in combination with the ballistic equations to solve the trajectory and then deduce the landing point of the trajectory.

[0032] Assume that the state variables are (x1, x2, …, x m ), the system of equations at the previous time step t n The value of each variable is (t n ,x 1,n ,x 2,n ,…,x m,n ), then at the next time step t n+1 The fourth-order Runge-Kutta algorithm formula for calculating each variable is:

[0033]

[0034] in

[0035]

[0036] where f i ′(t n ,x i,n ) is the derivative result of the state equation, and h is the time step.

[0037] A computer device comprises a memory, a processor and a computer program stored in the memory and executable on the processor, wherein the steps of the above method are implemented when the processor executes the program.

[0038] A computer-readable storage medium stores a computer program, which implements the steps of the above method when executed by a processor.

[0039] A computer program product comprises a computer program, which implements the steps of the above method when executed by a processor.

[0040] Compared with the prior art, the present invention has the following significant advantages:

[0041] (1) The present invention adopts the extended Kalman filter algorithm to make the nonlinear approximation more accurate, so as to make the tracking filter of the target more precise.

[0042] (2) The present invention uses the fourth-order Runge-Kutta method to solve the trajectory for landing point prediction, which reduces the trajectory extrapolation error and meets the requirement of improving the trajectory extrapolation accuracy.

[0043] (3) The point-of-fall prediction model adopted by the present invention is a particle trajectory model, which requires fewer measurement parameters and has a faster calculation speed. It can achieve rapid point-of-fall prediction and meet the real-time requirements of actual engineering. BRIEF DESCRIPTION OF THE DRAWINGS

[0044] Figure 1 The figure is a flow chart of the nonlinear filtering trajectory landing point prediction method based on the Runge-Kutta algorithm in the present invention.

[0045] Figure 2 This is the logical principle diagram of the track starting sliding window method.

[0046] Figure 3 It is the iterative calculation process of the fourth-order Runge-Kutta method.

[0047] Figure 4 A three-dimensional view of the ballistic landing point prediction results obtained based on measured data.

[0048] Figure 5 This is a top view of the ballistic landing point prediction results obtained based on the measured data. DETAILED DESCRIPTION

[0049] In order to further explain the technical means and effects adopted by the present invention to achieve the purpose of the invention, the following is a detailed description of a nonlinear filtering trajectory landing point prediction method based on the Runge-Kutta algorithm proposed by the present invention in combination with the accompanying drawings and specific implementation methods.

[0050] See also Figure 1 , Figure 1 This is a flowchart of an extended Kalman filter algorithm provided by an implementation example of the present invention to establish a track, and then use the fourth-order Runge-Kutta method to predict the landing point. The specific steps are as follows:

[0051] Step 1: Generate the point trace information of the artillery shell. Pre-process the artillery shell data measured by the radar system, convert the measured data from the radar polar coordinate system to the radar array rectangular coordinate system (RFC), and then convert the radar array rectangular coordinate system to the northeast celestial (ENU) coordinate system to obtain the point trace information of the artillery shell.

[0052] According to the position coordinates of the target in the radar polar coordinate system (θ A ,θ E ,d), where θ A is the azimuth of the target, θ E is the pitch angle of the target, and d is the distance from the target to the origin O of the radar coordinate system. Then the calculation formula of the target (x′, y′, z′) in the rectangular coordinate system of the radar array is as follows:

[0053]

[0054] The conversion of the radar array rectangular coordinate system to the northeast sky coordinate system requires two rotation processes. The rotation matrix is

[0055]

[0056] where ε A , ε E are the array azimuth and array elevation, R Z (ε A ), R X (ε E ) represent the rotation around the Z axis ε A Angle rotation matrix, rotate ε around the X axis E The rotation matrix of the angle.

[0057] The target position coordinates (x, y, z) in the northeast sky coordinate system are obtained as follows:

[0058]

[0059] Step 2: Construct a projectile particle trajectory model. Use the extended Kalman filter method to estimate the state of the projectile particle trajectory model.

[0060] Considering the projectile structural parameters, aerodynamic parameters, non-standard meteorological conditions and other information, the following particle trajectory model with time t as the independent variable in the northeast sky coordinate system is established:

[0061]

[0062] in, is the cross-sectional area of ​​the projectile, d is the projectile diameter; m is the mass of the projectile; ρ is the air density; c x (Ma) is the drag coefficient; w x , w y are longitudinal wind and cross wind respectively, g is the acceleration due to gravity; v r is the relative velocity between the projectile and the air, v x , v y , v z is the component of the projectile's center-of-mass velocity in the northeast celestial coordinate system.

[0063] The state variable at the last moment k of trajectory identification is

[0064] X(k)=[x k ,vx k ,y k ,vy k ,z k ,vz k ] T (13)

[0065] Then the state equation of the system obtained from the particle ballistic model is:

[0066] X(k+1)=F(k)X(k)+G(k)+V(k) (14)

[0067] Where F(k) is the state transfer matrix, G(k) is the coefficient matrix, V(k) is the process noise sequence, and is a Gaussian white noise sequence with zero mean and covariance Q(k).

[0068] Let T = t k+1 -t k ,have

[0069]

[0070] Then the one-step prediction of the state can be obtained as

[0071]

[0072] The one-step forecast of the covariance is

[0073] P(k+1|k)=[F(k)+G(k)f X (k)]P(k|k)[F(k)+G(k)f X(k)]′+Q(k)(18)

[0074] where f X (k) is the Jacobian matrix of F(k).

[0075]

[0076] The measurement equation is

[0077] Z(k)=h(X(k))+W(k) (20)

[0078] The measurement matrix h(X(k)) is

[0079]

[0080] Here W(k) is the measurement noise sequence, and it is assumed to be white Gaussian noise with zero mean and covariance R(k).

[0081] Then the measured prediction value can be obtained as

[0082]

[0083] The associated measurement prediction covariance is

[0084] S(k+1)=h X (k+1)P(k+1|k)h X ′(k+1)+R(k+1) (23)

[0085] Among them, h X (k) is the Jacobian matrix of h(k).

[0086] The gain is

[0087] K(k+1)=P(k+1|k)h X ′(k+1)S -1 (k+1) (24)

[0088] The state update equation is

[0089]

[0090] The covariance update equation is

[0091]

[0092] Where I is the identity matrix of the same dimension as the covariance.

[0093] Step 3: If there is no reliable track at present, go directly to step 4; if there is a reliable track, the preprocessed point track data and the reliable track are associated based on the extended Kalman filter algorithm.

[0094] First, a track that has not been processed in this inner loop is selected from the reliable tracks as the current track, and the last updated point track (associated or extrapolated) of the current track is used as the state X(k|k) and covariance P(k|k) of the previous moment.

[0095] Next, point-to-point association is performed based on the nearest neighbor algorithm (NNSF). First, a tracking gate is set, and the echo obtained by the preliminary screening of the ellipsoid gate becomes a candidate echo to limit the number of echoes involved in the correlation judgment. For the ellipsoid gate, if the conversion measurement value Z in the target rectangular coordinate system is c (k+1) satisfies

[0096]

[0097] The converted measurement value Z c (k+1) is the candidate echo, and equation (27) is called the ellipsoid gate rule. is the measurement prediction value, S(k+1) is the measurement prediction covariance, and the parameter γ is given by χ 2 The distribution table is obtained.

[0098] If there is only one measurement value that falls into the relevant gate, then the measurement value can be directly used for track update; however, if there is more than one echo that falls into the relevant gate of the tracked target, the candidate echo with the smallest statistical distance should be taken as the target echo, that is, in the nearest neighbor standard filter, the weighted norm of the innovation is

[0099]

[0100] The minimum amount of measurement is used to update the target state in the filter.

[0101] Finally, the first-order extended Kalman filter model updates the state and covariance of the reliable track at the current moment. The inner loop is to completely traverse all the reliable tracks once. The point tracks that are not associated at the current moment are put into step 4;

[0102] Step 4: According to the logical track initiation algorithm, the point track data that has not been associated is associated with the temporary track that has not yet set sail to achieve track initiation.

[0103] The logical method uses multiple assumptions to identify possible tracks through prediction and correlation gates. The steps are as follows:

[0104] Step 4.1: Use the measurements obtained in the first scan to establish a threshold for the track head, use the velocity method to establish an initial correlation gate, and establish possible tracks for the second scan measurements that fall into the initial correlation gate;

[0105] Step 4.2: Extrapolate each possible track, with the extrapolation point as the center. The size of the subsequent correlation gate is determined by the track extrapolation error covariance; the third scan measurement falls into the subsequent correlation gate closest to the extrapolation point and is interconnected;

[0106] Step 4.3: If there is no measurement in the subsequent related wave gate, cancel this possible track, or use the expanded related wave gate with acceleration limit to check whether the third scan measurement falls within it;

[0107] Step 4.4: Continue the above steps until a stable track is formed. The track initiation is considered complete.

[0108] Step 4.5: In all previous scans, the measurements that did not fall into the relevant gates and participate in the data interconnection judgment are taken as new track heads, and go to step 4.1. The schematic diagram is as follows: Figure 2 shown.

[0109] Step 4.6: Convert the temporary track that meets the required length into a reliable track.

[0110] The establishment of the temporary track of the three-coordinate radar realizes the estimation of the six-dimensional state vector through the initialization of the Kalman filter. The state variables of the trajectory identification are:

[0111] X(k)=[x k ,vx k ,y k ,vy k ,z k ,vz k ] T (29)

[0112] The measured value Z(k) in the northeast celestial coordinate system is

[0113]

[0114] At this time, the initial state of the system only needs to be determined by the measured values ​​Z(0) and Z(1) at the first two moments, that is,

[0115]

[0116] And the filter starts working from the moment k=2.

[0117] Step 5: If the reliable track is not updated at the current moment, that is, there is no point track associated with it, the reliable track is updated by track extrapolation, and the extrapolated value is the one-step prediction of the state and the one-step prediction covariance.

[0118] Step 6: Eliminate the reliable tracks that meet the extinction conditions (e.g., extrapolated n times);

[0119] Step 7: Use the Runge-Kutta method in numerical analysis for reliable trajectory, transform the predicted landing point into a discrete iterative process, combine the ballistic equations to solve the trajectory, gradually calculate the position and speed of the trajectory at each time point, and then calculate the trajectory landing point. The steps are as follows:

[0120] Step 7.1: Determine the initial conditions. First, we need to determine the initial conditions for trajectory calculation, including the target initial velocity and initial position. These initial conditions will be used as input parameters for the Runge-Kutta algorithm. Let the differential equations and initial values ​​be:

[0121]

[0122] x i (t0) = x i0 (i=1,2,…,m) (34)

[0123] Step 7.2: Set the time step. The time step is used to control the accuracy of the iteration. Generally speaking, the smaller the time step, the higher the calculation accuracy, but the amount of calculation will also increase accordingly.

[0124] Step 7.3: Iterative calculation. Based on the state of the previous time step, the predicted value is calculated in stages using linear interpolation. According to the predicted value and the differential equation, the fourth-order Runge-Kutta method uses four intermediate slope calculations to estimate the gradient at different time points and explicitly calculate the velocity increment and position increment at each stage, such as Figure 3 According to formula (35), the final new state is calculated by combining the results of all stages using weighted average as the starting point of the next time step; the values ​​of the variables in the equation group at the previous time step are (t n ,x 1,n ,x 2,n ,…,x m,n ), then the fourth-order Runge-Kutta algorithm formula for calculating each variable in the next time step is:

[0125]

[0126] in:

[0127]

[0128] Where h is the time step.

[0129] Step 7.4: Output the result. Continue iterating until the bull's eye height is reached, and then output the ballistic landing point.

[0130] The above algorithm is simulated based on real measurement data. The simulation object is a 460mm training projectile with a mass of 500kg. In the actual flight of the projectile, it will be affected by the wind, and the wind speed and direction are related to the height. In the simulation study, an average constant wind is used to replace the constantly changing wind under actual conditions, so that the influence of the constant wind on the trajectory of the projectile is equivalent to the ever-changing wind. This equivalent average wind is set as: longitudinal wind w x =1.82m / s, crosswind w y =3.44m / s, and the rest of the meteorological conditions are calculated according to the standard meteorological conditions for artillery. The actual target result is 4m, and the predicted landing point result is 6.99m. Figure 4 , Figure 5 As shown in the figure, segment AB is the tracking filter segment, the green circle is the measured data, and the red star is the updated value after filtering. Segment BC is the track extrapolation segment, point C reaches the radar blind area, segment CD is the landing point prediction segment, and the red cross at point D is the landing point prediction result. The blue cross is the actual landing point, and the red five-pointed star is the bull's eye. It can be seen from the figure that the error of the predicted landing point is less than 3 meters, which meets the engineering requirements.

[0131] The above contents are further detailed descriptions of the present invention in combination with specific preferred embodiments, and it cannot be determined that the specific implementation of the present invention is limited to these descriptions. For ordinary technicians in the technical field of the present invention, several simple deductions or substitutions can be made without departing from the concept of the present invention, which should be regarded as falling within the scope of protection of the present invention.

Claims

1. A nonlinear filtering trajectory drop point prediction method based on Runge-Kutta algorithm, characterized in that: The steps include: Step 1, pre-processing the point information of the shell; Step 2, constructing a projectile particle trajectory model; Step 3: If there is no reliable track at present, go directly to step 4; if there is a reliable track, perform point-to-point association between the preprocessed point track data and the reliable track based on the extended Kalman filter algorithm; Step 4, using the logical method track initiation algorithm, the point track data that has not been associated is associated with the temporary track that has not yet set sail, so as to realize the track initiation; Step 5: Extrapolate the track of the unupdated reliable track; Step 6, the reliable tracks that meet the extinction conditions are extincted; Step 7, repeat steps 1 to 6, and track the projectile in real time until the target reaches the radar blind spot. When the point trace cannot be obtained, the fourth-order Runge-Kutta method is used in combination with the ballistic equations to solve the trajectory and then calculate the trajectory landing point.

2. The nonlinear filtering trajectory drop point prediction method based on the Runge-Kutta algorithm according to claim 1 is characterized in that: Preprocess the point information of the shell; the rotation matrix of the radar array rectangular coordinate system converted to the northeast sky coordinate system is where ε A , ε E are the array azimuth and array elevation, R Z (ε A ), R X (ε E ) represent the rotation around the Z axis ε A The rotation matrix of angle, rotate ε around the X axis E The rotation matrix of the angle.

3. The nonlinear filtering trajectory drop point prediction method based on the Runge-Kutta algorithm according to claim 1 is characterized in that: The projectile particle trajectory model is constructed as follows: in, is the cross-sectional area of ​​the projectile, d is the projectile diameter; m is the mass of the projectile; ρ is the air density; c x is the drag coefficient; w x , w y are longitudinal wind and cross wind respectively, g is the acceleration due to gravity; v r is the relative velocity between the projectile and the air, v x , v y , v z is the component of the projectile's center-of-mass velocity in the northeast celestial coordinate system.

4. The nonlinear filtering trajectory drop point prediction method based on the Runge-Kutta algorithm according to claim 1 is characterized in that: When the point track data and the reliable track are associated, according to the ellipsoid gate rule, if the The converted measurement value Z c (k+1) is the candidate echo; is the measurement prediction value, S(k+1) is the measurement prediction covariance, and the parameter γ is given by χ 2 The distribution table is obtained.

5. The nonlinear filtering trajectory drop point prediction method based on the Runge-Kutta algorithm according to claim 4 is characterized in that: When the point-to-point data and reliable tracks are associated based on the extended Kalman filter algorithm, the weighted norm of the new information is used according to the nearest neighbor criterion. Achieve the minimum conversion measurement value Z c (k+1) will be used in the filter to update the target state.

6. The nonlinear filtering trajectory drop point prediction method based on the Runge-Kutta algorithm according to claim 5 is characterized in that: When the point-to-point data and the reliable track are associated based on the extended Kalman filter algorithm, the reliable track is updated according to the state update equation and the covariance update equation; the state equation is: in is the one-step prediction of the state, K(k+1) is the Kalman gain matrix, is the measured predicted value, Z(k+1) is the measured value; The covariance update equation is Where I is the identity matrix of the same dimension as the covariance, P(k+1|k) is the one-step prediction of the covariance, K′(k+1) is the transpose of the Kalman gain matrix, and h x (k) is the Jacobian matrix of the measurement equation, and R(k+1) is the measurement noise covariance.

7. The nonlinear filtering trajectory drop point prediction method based on the Runge-Kutta algorithm according to claim 1 is characterized in that: The fourth-order Runge-Kutta method is combined with the ballistic equations to solve the trajectory and then deduce the trajectory point; Let the state variables be (x1, x2, ..., x m ), the system of equations at the previous time step t n The value of each variable is (t n ,x 1,n ,x 2,n ,...,x m,n ), then at the next time step t n+1 The fourth-order Runge-Kutta algorithm formula for calculating each variable is: in where f i ′(t n ,x i,n ) is the derivative result of the state equation, and h is the time step.

8. A computer device comprising a memory, a processor and a computer program stored in the memory and executable on the processor, characterized in that: When the processor executes the program, the steps of the method according to any one of claims 1 to 7 are implemented.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the program is executed by a processor, the steps of the method described in any one of claims 1 to 7 are implemented.

10. A computer program product, comprising a computer program, characterized in that When the computer program is executed by a processor, the steps of the method according to any one of claims 1 to 7 are implemented.

Images