Transactions of Nanjing University of Aeronautics and Astronautics  2018, Vol. 35 Issue (2): 361-370   PDF    
An Algorithm for Determination of Projectile Attitude Angles in Projectile Trajectory Prediction
Zha Qicheng, Rui Xiaoting, Yu Hailong, Zhou Qinbo     
Institute of Launch Dynamics, School of Energy and Power Engineering, Nanjing University of Science and Technology, Nanjing 210094, P. R. China
Abstract: A fast and accurate algorithm is established in this paper to increase the precision of ballistic trajectory prediction. The algorithm is based on the six-degree-of-freedom (6DOF) trajectory equations, to estimate the projectile attitude angles in every measuring time. Hereby, the algorithm utilizes the Davidon-Fletcher-Powell (DFP) method to solve nonlinear equations and Doppler radar trajectory test information containing only position coordinates of the projectile to reconstruct the angular information. The ″position coordinates by the test″ and ″angular displacements by reconstruction″ at the end phase of the radar measurement are used as an initial value for the trajectory computation to extrapolate the trajectory impact point. The numerical simulations validate the proposed method and demonstrate that the estimated impact point agrees very well with the real one. Morover, other artillery trajectory can be predicted by the algorithm, and other trajectory models, such as 4DOF and 5DOF models, can also be incorporated into the proposed algorithm.
Key words: ballistic trajectory    impact point prediction    Davidon-Fletcher-Powell method    six-degree-of-freedom trajectory model    Doppler radar data    
0 Introduction

Artillery weapons have drawn great attention from many countries since they are praised as the ″god of war″ and plays a very important role in land battles[1]. The trajectory prediction, which is tantamount to acquire the position of the origin or impact point of a trajectory using the radar data, is extremely significant for a artillery combat chain[2]. The accurate prediction of the impact point can be beneficial to the power of artillery, and can improve the accuracy of shooting and its amending[3]. To predict the impact point of ballistic projectiles is a very critical and relevant problem which has applications in area defense[4].

One of the earliest studies of trajectory prediction and determination was conducted by Johann Heinrich Lambert (1728—1779). His method was used for spacecraft trajectory prediction in order to solve an orbital boundary value problem. Although this research effort focused on the prediction of ballistic trajectories, Lambert′s theorem is helpful to calculate the elliptical parameters of the trajectory to generate an estimate. Utilizing this method, the transfer trajectory between two points can be calculated given a transfer time. The technique is based upon an iterative solution that converges on the semi-major axis of the transfer ellipse and is relevant to the problem considered in this study since it provides a methodology for estimating a trajectory as a solution to a two-point boundary value problem[5].

Isaacson and Vaughan[6] described a method of estimating and predicting ballistic trajectories using a Kalman filter. A template was developed using missile range and altitude data by modeling the ballistic missile flight over a spherical, non-rotating Earth. The trajectory was used as a baseline from which perturbations of the actual trajectory were calculated. The missile trajectory was obsened using a pair of geosynchronous satellites. The method was used to determine uncertainties in the missile launch point and missile position during flight.

Danis[7] developed a method for estimating ballistic missile launch parameters using a single pair of angle measurements taken from two satellites, or utilizing both measurements taken from a single satellite. Using a priori information about the trajectory and launch time, the launch location and missile heading were calculated from the geometric relationship between the missile and satellite position. If sufficient measurements were taken during the missile flight, a Kalman filter was used to assess measurement errors and refine the launch parameter estimation.

Eric et al.[8] used the three-degree-of-freedom trajectory equations for trajectory prediction. Khalil et al. [9, 10] established the modified point mass model for projectile trajectory based on the discrete time transfer matrix method, and further reconstructed the repose angle at any measuring time instant. Wu et al.[11] used the ADS-B fusion algorithm for trajectory prediction. No matter what the method was, the purpose was to rapidly obtain an accurate impact point. It should be mentioned that different ballistic models have different calculation precisions.

Further, the precision of the initial value has a great influence on the calculation precision of all kinds of external ballistic equations[12]. The common models used for the motion of a projectile include: Point mass trajectory model (the 3DOF model, herein, DOF is the degree-of-freedom), modified point mass trajectory model (the 4DOF model), and rigid body trajectory model (the 6DOF model). These models are different from each other in their preliminaries, physical assumptions, aerodynamic parameters involved, computational speeds and accuracies. Therefore, one should make a reasonable choice according to a specific situation. The 3DOF model is only concerned with three position coordinates of a projectile, and used to rapidly estimate the approximate impact point of a relatively poor precision. The 4DOF model further considers the repose angle of a projectile[13], and is generally used when the firing angle is less than 50°, or when the ascent stage of the anti-aircraft gun and the navy gun are against air targets[14]. The 6DOF model not only considers the projectile′s position using three coordinates, but also its orientation using three angular displacements. The 6DOF model has the best accuracy when the initial values of the integral variables are precisely known. The difficulty of applying 6DOF model lies in the acquisition of these initial values. A possible way is to utilize the trajectory radar which can measure position and velocity of a projectile. However, the information of orientation and angular velocity is not included.

It is well known that the Kalman filter is widely used in parameter estimation for predictive target tracking[15-17]. However, if there are position coordinates by the radar test, instead of angle measurements, the projectile attitude angles can not be estimated by the Kalman filter.

In order to improve the precision of ballistic trajectory prediction, we propose an algorithm to estimate the projectile attitude angles based on the 6DOF trajectory model. Hereby, the algorithm utilizes the Davidon-Fletcher-Powell (DFP) method in solving nonlinear equations and radar trajectory test information containing only position coordinates of the projectile to reconstruct the angular displacements.

The DFP formula finds the solution to the secant equation that is the closest to the current estimate and satisfies the curvature condition. It was the first quasi-Newton method to generalize the secant method to a multi-dimensional problem. This update maintains the symmetry and positive definiteness of the Hessian matrix[18].

Due to the particularity of the artillery, and missile environment of high temperature and high pressure, high overload, high speed rotation, and space and the deficiency of the sensor is expensive, installed the sensor on the projectile is difficult. Fewer people are choosing to install sensors on ordinary projectiles. In this paper, an algorithm based on the six-degree-of-freedom (6DOF) trajectory equations is developed and presented to estimate the projectile attitude angles in every measuring time. The proposed method is proved by experiment and the experiment demonstrates that the algorithm has good prediction effect using the 6DOF trajectory model.

The passive trajectory of other artillery or use other trajectory models such as 4DOF and 5DOF models can also be incorporated into the proposed algorithm.

1 Mathematical Model of Exterior Ballistic

Considering the dynamic unbalance, mass eccentric and gust of wind, the equations of motion of the ordinary projectile can be written as

$ \left\{ \begin{array}{l} \frac{{{\rm{d}}x}}{{{\rm{d}}t}} = {v_x},\frac{{{\rm{d}}y}}{{{\rm{d}}t}} = {v_y},\frac{{{\rm{d}}z}}{{{\rm{d}}t}} = {v_z}\\ \frac{{{\rm{d}}\gamma }}{{{\rm{d}}t}} = \dot \gamma ,\frac{{{\rm{d}}{\varphi _a}}}{{{\rm{d}}t}} = {{\dot \varphi }_a},\frac{{{\rm{d}}{\varphi _2}}}{{{\rm{d}}t}} = {{\dot \varphi }_2} \end{array} \right. $ (1)
$ \begin{array}{*{20}{c}} {\frac{{{\rm{d}}{v_x}}}{{{\rm{d}}t}} = - \frac{{{R_x}}}{m}\frac{{{v_x} - {w_x}}}{{{v_r}}} - \frac{{{R_y}}}{{m\sin {\delta _r}}} \cdot }\\ {\left( {\sin {\delta _{{r_1}}}\cos {\delta _{{r_2}}}\sin {\theta _r} + \sin {\delta _{{r_2}}}\sin {\psi _r}\cos {\theta _r}} \right) + }\\ {\frac{{{R_z}}}{{m\sin {\delta _r}}}\left( {\sin {\psi _r}\cos {\theta _r}\cos {\delta _{{r_2}}}\sin {\delta _{{r_1}}} - \sin {\theta _r}\sin {\delta _{{r_2}}}} \right)} \end{array} $ (2)
$ \begin{array}{*{20}{c}} {\frac{{{\rm{d}}{v_y}}}{{{\rm{d}}t}} = - \frac{{{R_x}}}{m}\frac{{{v_y}}}{{{v_r}}} + \frac{{{R_y}}}{{m\sin {\delta _r}}} \cdot }\\ {\left( {\sin {\delta _{{r_1}}}\cos {\delta _{{r_2}}}\cos {\theta _r} - \sin {\delta _{{r_2}}}\sin {\psi _r}\sin {\theta _r}} \right) + }\\ {\frac{{{R_z}}}{{m\sin {\delta _r}}}\left( {\sin {\psi _r}\sin {\theta _r}\cos {\delta _{{r_2}}}\sin {\delta _{{r_1}}} + \cos {\theta _r}\sin {\delta _{{r_2}}}} \right) - g} \end{array} $ (3)
$ \begin{array}{*{20}{c}} {\frac{{{\rm{d}}{v_z}}}{{{\rm{d}}t}} = - \frac{{{R_x}}}{m}\frac{{{v_z} - {w_z}}}{{{v_r}}} + \frac{{{R_y}}}{{m\sin {\delta _r}}}\sin {\delta _{{r_2}}}\cos {\psi _r}}\\ { - \frac{{{R_z}}}{{m\sin {\delta _r}}}\cos {\psi _r}\cos {\delta _{{r_2}}}\sin {\delta _{{r_1}}}g} \end{array} $ (4)
$ \frac{{{\rm{d}}\dot \gamma }}{{{\rm{d}}t}} = - {{\ddot \varphi }_a}\sin {\varphi _2} - {{\ddot \varphi }_a}{{\dot \varphi }_2}\cos {\varphi _2} - {k_{xD}}\left( {\dot \gamma + {{\dot \varphi }_a}\sin {\varphi _2}} \right){v_r} $ (5)
$ \begin{array}{*{20}{c}} {\frac{{{\rm{d}}{{\dot \varphi }_a}}}{{{\rm{d}}t}} = \frac{{{M_z}\left( {\cos {\delta _{{r_1}}}\sin {\delta _{{r_2}}}\sin {\alpha _r} + \sin {\delta _{{r_1}}}\cos {\alpha _r}} \right) + {M_y}\left( {\cos {\delta _{{r_1}}}\sin {\delta _{{r_2}}}\cos {\alpha _r} - \sin {\delta _{{r_1}}}\cos {\alpha _r}} \right)}}{{A\cos {\varphi _2}\sin {\delta _r}}} + }\\ {\frac{{\left( {2A - C} \right){{\dot \varphi }_a}{{\dot \varphi }_2}\sin {\varphi _2} - C\dot \gamma {{\dot \varphi }_2}}}{{A\cos {\varphi _2}}} - {k_{zD}}{{\dot \varphi }_a}\cos {\varphi _2}{v_r} + \left( {1 - \frac{C}{A}} \right)\left( {\ddot \gamma {\beta _{D\xi }} + {{\dot v}_x}{\beta _{D\eta }}} \right) + \frac{{m{L_{{m_\eta }}}{{\dot v}_x}}}{A}} \end{array} $ (6)
$ \begin{array}{*{20}{c}} {\frac{{{\rm{d}}{{\dot \varphi }_2}}}{{{\rm{d}}t}} = \frac{{{M_z}\left( {\cos {\delta _{{r_1}}}\sin {\delta _{{r_2}}}\cos {\alpha _r} - \sin {\delta _{{r_1}}}\sin {\alpha _r}} \right) - {M_y}\left( {\cos {\delta _{{r_1}}}\sin {\delta _{{r_2}}}\cos {\alpha _r} + \sin {\delta _{{r_1}}}\cos {\alpha _r}} \right)}}{{A\sin {\delta _r}}} - }\\ {\left( {1 - \frac{C}{A}} \right)\dot \varphi _a^2\cos {\varphi _2}\sin {\varphi _2} + \frac{{C\dot \gamma {{\dot \varphi }_a}\cos {\varphi _2}}}{A} + - {k_{zD}}{{\dot \varphi }_2}{v_r} + \left( {1 - \frac{C}{A}} \right)\left( {{{\dot \gamma }^2}{\beta _{{D_\xi }}} - \ddot \gamma {\beta _{{D_\eta }}}} \right) + \frac{{m{L_{{m_\xi }}}{{\dot v}_x}}}{A}} \end{array} $ (7)

where x, y, z are the coordinates of the projectile geometry center in the ground coordinate system; vx, vy, vz the components of the velocity of the projectile geometry center in the ground coordinate system; m is the mass of the projectile, A the equatorial rotational inertia of the projectile, C the polar moment of inertia of the projectile, g the gravitational acceleration, vr the relative atmospheric velocity, wx the velocity of the longitudinal wind, wz the the velocity of the cross wind, $ {{\dot v}_r}$ the relative acceleration, Rx the drag force, Ry the lift force, Rz the Magnus force, ψr the relative deflection angle, θr the relative ballistic inclination angle, δr the relative angle of attack, δr1 the projection of the relative angle of attack on the vertical plane, δr2 the projection of the relative angle of attack on the lateral plane, $ {{\dot \psi }_r}$ the relative declination angle velocity, $ {\dot \theta }$ the trajectory inclination angle velocity, αr the relative aiming angle, kxD the polar damping moment coefficient, kzD the equator torque damping coefficient, Mz the Magnus moment, and Mz the stable moment. βDη and βDζ are the projections of dynamic unbalance angle βD in axis η and ζ of the o1ξηζ, respectively, Lmηand Lmζ the projections of mass eccentricity Lm in axis η and ζ of the o1ξηζ, respectively.

2 Determination Algorithm for Projectile Attitude Angles

When computing the 6DOF trajectory model, one needs to provide 12 independent initial values including (x, y, z, vx, vy, vz, γ, φa, φ2, ${\dot \gamma } $, ${{\dot \varphi }_a} $, ${{\dot \varphi }_2} $). On one hand, however, only the first 6 quantities, namely (x, y, z, vx, vy, vz), can be obtained from radar test data; on the other hand, the initial value of the spin angle γ may not influence the motion of the projectile since it does not explicitly appear in the motion equations, Eqs.(7)—(12). Consequently, how to obtain the angular velocity of spinning, two angles of oscillation, and two angular velocities of oscillation ($ {\dot \gamma }$, φa, φ2, ${{\dot \varphi }_a} $, $ {{\dot \varphi }_2}$), will be studied in this paper.

If the radar data at time instant ti and ti+1 are known, one can expand the positions of a projectile using a second-order Taylor series at ti as

$ {X_{i + 1}} = {X_i} + {V_i}\Delta t + \frac{1}{2}{A_i}\Delta {t^2} $ (8)

where Xi and Xi+1 are the position coordinates of the projectile at ti and ti+1, respectively; Vi is the velocities of the projectile at ti, Δt=ti+1ti the difference between ti and ti+1, Ai is the accelerations of the projectile at ti.

Further, it can be deduced that ${A_i} = \frac{2}{{\Delta {t^2}}}({X_{i + 1}} - {X_i} - {V_i}\Delta t) $. Moreover, Ai can also be expressed as ${\left( {\frac{{{\rm{d}}{v_x}}}{{{\rm{d}}t}}, \frac{{{\rm{d}}{v_y}}}{{{\rm{d}}t}}, \frac{{{\rm{d}}{v_z}}}{{{\rm{d}}t}}} \right)_i} $, where the (x, y, z, vx, vy, vz)i at ti and (x, y, z, vx, vy, vz)i+1 at ti+1 can be obtained from radar measurement. On the right side of Eqs.(7)—(9), only the angular velocity of spinning and the two angles of oscillation (${\dot \gamma } $, φa, φ2)i are unknown and independent, therefore, $\left( {\frac{{{\rm{d}}{v_x}}}{{{\rm{d}}t}}, \frac{{{\rm{d}}{v_y}}}{{{\rm{d}}t}}, \frac{{{\rm{d}}{v_z}}}{{{\rm{d}}t}}} \right) $ can be considered as a nonlinear function with respect to (${\dot \gamma } $, φa, φ2)i at ti, namely $ {F_i}(\dot \gamma, {\varphi _a}, {\varphi _2}) = {A_i} - \frac{2}{{\Delta {t^2}}}({X_{i + 1}} - {X_i} - {V_i}\Delta t)$. By solving (${\dot \gamma } $, φa, φ2)i, the problem is transformed into solving a ternary nonlinear equations as follows

$ \left\{ \begin{array}{l} {F_{xi}}\left( {\dot \gamma ,{\varphi _a},{\varphi _2}} \right) = 0\\ {F_{yi}}\left( {\dot \gamma ,{\varphi _a},{\varphi _2}} \right) = 0\\ {F_{zi}}\left( {\dot \gamma ,{\varphi _a},{\varphi _2}} \right) = 0 \end{array} \right. $ (9)

Solving nonlinear equations have a wide range of applications. Typically, the nonlinear equation is numerically solved by iteration methods. Common methods include the fixed point iteration method, Newton iterative method and Quasi Newton iterative method.

DPF method avoids obtaining inverse matrix in its iteration formula, thus increases the speed of calculation[19].

Given a function f(x), its gradient $\left( {\nabla \mathit{\boldsymbol{f}}} \right) $, and positive definite Hessian matrix B, the Taylor series is

$ \mathit{\boldsymbol{f}}\left( {{\mathit{\boldsymbol{x}}_k} + {\mathit{\boldsymbol{s}}_k}} \right) = \mathit{\boldsymbol{f}}\left( {{\mathit{\boldsymbol{x}}_k}} \right) + \nabla \mathit{\boldsymbol{f}}{\left( {{\mathit{\boldsymbol{x}}_k}} \right)^{\rm{T}}}{\mathit{\boldsymbol{s}}_k} + \frac{1}{2}\mathit{\boldsymbol{s}}_k^{\rm{T}}\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{s}}_k} $ (10)

and the Taylor series of the gradient itself (secant equation)

$ \nabla \mathit{\boldsymbol{f}}\left( {{\mathit{\boldsymbol{x}}_k} + {\mathit{\boldsymbol{s}}_k}} \right) = \nabla \mathit{\boldsymbol{f}}\left( {{x_k}} \right) + \mathit{\boldsymbol{B}}{\mathit{\boldsymbol{s}}_k} $ (11)

is used to update B. The DFP formula finds a solution that is symmetric, positive definite and closest to the current approximate value of Bk

$ {\mathit{\boldsymbol{B}}_{k + 1}} = \left( {\mathit{\boldsymbol{I}} - {\gamma _k}{\mathit{\boldsymbol{y}}_k}\mathit{\boldsymbol{s}}_k^{\rm{T}}} \right){\mathit{\boldsymbol{B}}_k}\left( {\mathit{\boldsymbol{I}} - {\gamma _k}{\mathit{\boldsymbol{s}}_k}\mathit{\boldsymbol{y}}_k^{\rm{T}}} \right) + {\gamma _k}{\mathit{\boldsymbol{y}}_k}\mathit{\boldsymbol{y}}_k^{\rm{T}} $ (12)

where

$ {\mathit{\boldsymbol{y}}_k} = \nabla \mathit{\boldsymbol{f}}\left( {{\mathit{\boldsymbol{x}}_k} + {\mathit{\boldsymbol{s}}_k}} \right) - \nabla \mathit{\boldsymbol{f}}\left( {{\mathit{\boldsymbol{x}}_k}} \right) $ (13)
$ {\gamma _k} = \frac{1}{{\mathit{\boldsymbol{y}}_k^{\rm{T}}{\mathit{\boldsymbol{s}}_k}}} $ (14)

and Bk is a symmetric and positive definite matrix. The corresponding update to the inverse Hessian approximation Hk=Bk-1 is given by

$ {\mathit{\boldsymbol{H}}_{k + 1}} = {\mathit{\boldsymbol{H}}_k} - \frac{{{\mathit{\boldsymbol{H}}_k}{\mathit{\boldsymbol{y}}_k}\mathit{\boldsymbol{y}}_k^{\rm{T}}{\mathit{\boldsymbol{H}}_k}}}{{\mathit{\boldsymbol{y}}_k^{\rm{T}}{\mathit{\boldsymbol{H}}_k}{\mathit{\boldsymbol{y}}_k}}} + \frac{{{\mathit{\boldsymbol{s}}_k}\mathit{\boldsymbol{s}}_k^{\rm{T}}}}{{\mathit{\boldsymbol{y}}_k^{\rm{T}}{\mathit{\boldsymbol{s}}_k}}} $ (15)

B is assumed to be positive definite, and the vectors skT and y must satisfy the curvature condition

$ \mathit{\boldsymbol{s}}_k^{\rm{T}}{\mathit{\boldsymbol{y}}_k} = \mathit{\boldsymbol{s}}_k^{\rm{T}}\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{s}}_k} > 0 $ (16)

We solve ${(\dot \gamma, {\varphi _a}, {\varphi _2})_i} $ based on the DFP method.

In a similar way, known the radar data at ti+1 and ti+2, one can calculate $ {(\dot \gamma, {\varphi _a}, {\varphi _2})_{i + 1}}$. Hereby, the $ {({{\dot \varphi }_a}, {{\dot \varphi }_2})_i}$ at ti can be approximately evaluated by

$ {\left( {{{\dot \varphi }_a},{{\dot \varphi }_2}} \right)_i} \approx \frac{{{{\left( {{\varphi _a},{\varphi _2}} \right)}_{i + 1}} - {{\left( {{\varphi _a},{\varphi _2}} \right)}_i}}}{{{t_{i + 1}} - {t_i}}} $ (17)

The schematic diagram for a cycle is shown in Fig. 1.

Fig. 1 Calculating diagram within a loop

Specific steps are as follows:

Step 1 Acquire the radar data from the measurement on the coordinates and velocities of the projectile at titi+1 and ti+2, namely (X, V)i, (X, V)i+1, (X, V)i+2;

Step 2 Use (X, V)i and (X, V)i+1 to calculate $ {(\dot \gamma, {\varphi _a}, {\varphi _2})_i}$ by DFP method;

Step 3 Use (X, V)i+1 and (X, V)i+2 to calculate ${(\dot \gamma, {\varphi _a}, {\varphi _2})_{i + 1}} $ by DFP method;

Step 4 Use (φa, φ2)i and (φa, φ2)i+1 to calcu-late ${({{\dot \varphi }_a}, {{\dot \varphi }_2})_i} $ by difference quotient.

Up to now, ${(\dot \gamma, {\varphi _a}, {\varphi _2}, {{\dot \varphi }_a}, {{\dot \varphi }_2})_i} $ has been obtained. Its basic idea is to use three consecutive test data of position and velocity of the projectile to solve the angles on the first point. Along with the position and velocity data by radar measurement, the initial values (x, y, z, vx, vy, vz, $ {\dot \gamma }$, φa, φ2, $ {{\dot \varphi }_a}, {{\dot \varphi }_2}$)i of exterior ballistic equations are all known. Further, the 6DOF trajectory model can be solved, and the impact point can be computed. The cycle stops when n is the last point of radar test data.

The flow chart of the algorithm is shown in Fig. 2.

Fig. 2 Flow chart of 6DOF trajectory prediction method

3 Smulation and Experimental Verification 3.1 Simulation verification of the projectile′s angles and angular velocities

Since the radar test data does not contain angles information, a simulation method is adopted to verify the calculation precision of the angles. Hereby, a 155 mm spin-stabilized projectile is studied for ballistic simulation. The mass of the projectile is 45.678 kg; its equator moment of inertia is 1.785 kg·m2, polar moment of inertia 0.164 3 kg·m2, twist rate of rib rifling 0.141 260 16, length of rifling 6.41 m. The atmospheric density, the variation of speed of sound, and aerodynamic coefficient can be referred to Ref.[1].

The ballistic data can be generated by using the 6DOF trajectory model as ″the simulated radar data″.

First of all, the trajectory simulation is carried out based on three elevation angles including35°, 49°, 65°(low, common, high), and direction angle 0°. Basic parameters are shown in the Table 1.

Table 1 Simulation data of the ballistic

Assume the simulated trajectory as a ″simulated radar data″, and take the results of the simulation from 1 to 2 s after launch, the time step is set to be 0.01 s. Thereby, total 100 positions and velocity coordinates, as well as orientation angles and their time derivatives will be generated at the 100 time instant. By substituting these 100 ″simulated radar data″ of position and velocity coordinates into the proposed attitude angles determination algorithm, one can finally compare the calculated results with the ″simulated radar data″ of attitude angles and their time derivatives.

Figs. 311 show the comparison between the ″simulated radar data″ and estimated results by the algorithm in the case of 35° low elevation angle, 49° common elevation angle and 65° high elevation angle, respectively. It can be seen that the estimated angles and their time derivatives using the algorithm 6DOF (6 degree of freedom estimated) are in good agreement with the simulated radar data, which validates the algorithm in the case of the common, high and low elevation angles, respectively.

Fig. 3 The spin rate versus time for comparison in the case of 35° elevation angle

Fig. 4 The pitch angle versus time for comparison in the case of 35° elevation angle

Fig. 5 The yaw angle versus time for comparison in the case of 35° elevation angle

Fig. 6 The spin rate versus time for comparison in the case of 49° elevation angle

Fig. 7 The pitch angle versus time for comparison in the case of 49° elevation angle

Fig. 8 The yaw angle versus time for comparison in the case of 49° elevation angle

Fig. 9 The spin rate versus time for comparison in the case of 65° elevation angle

Fig. 10 The pitch angle versus time for comparison in the case of 65° elevation angle

Fig. 11 The yaw angle versus time for comparison in the case of 65° elevation angle

3.2 Results of ballistic test

In order to further verify the algorithm in this article, some firing tests have been implemented with an elevation angle of 45° and a direction angle of 0°.

Hereby impact points of four projectiles are shown in Table 2.

Table 2 Four impact points of projectiles in the firing test

In the aforementioned table, the first 10th measuring filtered points from Doppler radar test for the first projectile are shown in Table 3.

Table 3 The first 10th measuring filtered points of the first projectile from Doppler radar test

3.3 Test verification on projectile′s position and velocity coordinates

Substituting the Doppler radar test data in Table 3 into the trajectory prediction algorithm and computing the impact point using the 6DOF trajectory model with the estimated initial values, one can compare the actual test data (Real radar data) with the simulation data (6DOF estimated), as shown in Figs. 1216.

Fig. 12 Range versus time for comparison in the case of 45° elevation angle

Fig. 13 Altitude versus time for comparison in the case of 45° elevation angle

Fig. 14 Velocity x versus time for comparison in the case of 45° elevation angle

Fig. 15 Velocity y versus time for comparison in the case of 45° elevation angle

Fig. 16 Velocity z versus time for comparison in the case of 45° elevation angle

3.4 Test verification on the prediction of projectile′s impact point

The error of projectile′s impact point can be defined as

$ \varepsilon = \frac{{{\rm{estimated}} - {\rm{real}}}}{{{\rm{real}}}} \times 100 $ (18)

where ESTIMATED is the impact point of estimated by the proposed algorithm, and REAL the real radar date.

Table 4 shows the error analysis of the prediction of the projectile′s impact point with an elevation angle of 45°. The error between the algorithm results and real data is small, with a slight deviation in the Z direction. From the meteorological conditions of the test, the high altitude wind speed, which can reach 62—77 m/s in the vertex height, is very large.In such case, high altitude wind measurement error will be relatively large. However, the accuracy of the prediction is acceptable in engineering scale.This deviation may be caused due to two reasons:One is measuring error of meteorological conditions; the other is the poor state of the artillery used, resulting in a relatively large initial disturbance.

Table 4 Error analysis of the prediction impact point of 45° elevation angle

The experiment validates that the algorithm has good prediction effect using the 6DOF trajectory model.

4 Conclusions

In this paper, a determination algorithm used DFP method to solve nonlinear equations is applied to estimate projectile′s attitude angles at the end phase of the radar data of the projectile trajectory test. We established the determination algorithm of the projectile′s attitude angles for trajectory prediction, and extrapolated the impact point. In order to verify the proposed method, the trajectory simulation was implemented based on three elevation angles, 35°, 49°, 65° as the ″simulated radar data″. By substituting 100 ″simulated radar data″ of position and velocity coordinates into the proposed attitude angles determination algorithm, one can compare the estimated results with the ″simulated radar data″ of attitude angles and their time derivatives. Good agreement has been achieved in the result. The firing test is carried out and the real radar data are substituted into the algorithm for trajectory prediction. The prediction and the actual curves agree very well, and the error is acceptable within engineering scope.

Not only the projectile but also the other artillery trajectory can be predicted by the algorithm. The proposed algorithm can also be extended to trajectory prediction using other trajectory models such as 4DOF and 5DOF models. This article provides a novel technical foundation for evaluating the dynamic performance and improving the accuracy of artillery trajectory prediction.

Acknowledgements

The research was supported by the Research Fund for the Doctoral Program of Higher Education of China (No.20133219110037); the Natural Science Foundation of China(No.11472135); and the Program for New Century Excellent Talents in University (No.NCET-10-0075).

References
[1]
RUI X T, LIU Y X, YU H L. Launch dynamics of tank and self-propelled artillery[M]. Beijing: Science Press, 2008.
[2]
FRESCONI F, COOPER G, COSTELLO M. Practical assessment of real-time impact point estimators for smart weapons[J]. Journal of Aerospace Engineering, 2011, 24(1): 1-11. DOI:10.1061/(ASCE)AS.1943-5525.0000044
[3]
SHAO X, WANG H, ZHU L. A potential efficiency assessment model for reconnaissance and spotting radar[J]. Fire Control and Command Control, 2008, 7(33): 40-43.
[4]
KRAMER K A, STUBBERUD S C. Impact time and point predicted using a neural extended Kalman filter[C]//IEEE Intelligent Sensors, Sensor Networks and Information Processing Conference. Melbourne, Australia: IEEE, 2008: 199-204.
[5]
PRUSSING J E, CONWAY B A. Orbital mechanics[M]. New York: Oxford University Press, 1993.
[6]
ISAACSON J A, VAUGHAN D R. Estimation and prediction of ballistic missile trajectories[J]. Journal of Clinical Psychology, 1996, 31(4): 770-775.
[7]
DANIS N J. Space based tactical ballistic missile launch parameter estimation[J]. IEEE Transactions on Aerospace and Electronic Systems, 1993, 29(2): 412-424. DOI:10.1109/7.210079
[8]
ERIC N, MEIR P, STANTON M. Projectile launch point estimation from radar measurement[C]//American Control Conference. Portland, USA: [s. n. ], 2005: 1275-1282.
[9]
KHALIL M, RUI X T, HENDY H. Discrete time transfer matrix method for projectile trajectory prediction[J]. Journal of Aerospace Engineering, 2015, 28(2): 04014057. DOI:10.1061/(ASCE)AS.1943-5525.0000381
[10]
KHALIL M, RUI X T, ZHA Q C, et al. Projectile impact point prediction based on self-propelled artillery dynamics and Doppler radar measurements[J]. Advances in Mechanical Engineering, 2013(3): 153913.
[11]
WU Z, WANG M, ZHANG R. A new algorithm to real-time optimal estimation of radar biases based on ADS-B[J]. Journal of Southwest Jiaotong University, 2013, 48(1): 1-6.
[12]
MCCOY R. Modern exterior ballistics:The launch and flightdynamics of symmetric projectiles[M]. USA: Schiffer, 2011: 1-11.
[13]
LIESKE R, REITER M. Equations of motion for a modified point mass trajectory: U. S. Army Ballistic Research Laboratory Rep: R1314(AD 485869)[R]. USA: U. S. Army, Aberdeen Proving Ground, 1966.
[14]
PU F, RUI X T. Exterior ballistics[M]. Beijing: National Defense Industry Press, 1989.
[15]
RAVINDRA V C, BAR-SHALOM Y, WILLETT P. Projectile identification and impact point prediction[J]. IEEE Transactions on Aerospace and Electronic Systems, 2007, 46(4): 2004-2021.
[16]
YANG L Q, XIAO Q G, NIU Y, et al. Design of localization system based on reducing Kalman filter[J]. Journal of Nanjing University of Aeronautics and Astronautics, 2012, 44(1): 134-138.
[17]
CAI J, ZHANG L, DONG P. Remaining useful life prediction for aero-engines combining sate space model and kf algorithm[J]. Transactions of Nanjing University of Aeronautics and Astronautics, 2017, 34(3): 265-271.
[18]
FLETCHER R, POWELL M J D. A rapidly convergent descent method for minimization[J]. The Computer Journal, 1963, 6(2): 163-168. DOI:10.1093/comjnl/6.2.163
[19]
ZHAO S. Application of MATLAB programming and optimization design[M]. Beijing: Electronics Industry Publishing House, 2013.