An Algorithm for Determination of Projectile Attitude Angles in Projectile Trajectory Prediction

， ， ，

Institute of Launch Dynamics，School of Energy and Power Engineering，Nanjing University of Science and Technology，Nanjing 210094，P．R．China

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 meth－od 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，nonrotating Earth．The trajectory was used as a baseline from which perturbations of the actual trajectory were calculated．The missile trajectory was obsened using apair 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 apriori 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 agreat 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 aprojectile［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

wherex，y，zare the coordinates of the projectile geometry center in the ground coordinate system；vx，vy，vzthe components of the velocity of the projectile geometry center in the ground coordinate system；mis the mass of the projectile，Athe equatorial rotational inertia of the projectile，Cthe polar moment of inertia of the projectile，gthe gravitational acceleration，vrthe relative atmospheric velocity，wxthe velocity of the longitudinal wind，wzthe the velocity of the cross wind，vrthe relative acceleration，Rxthe drag force，Rythe lift force，Rzthe Magnus force，ψrthe relative deflection angle，θrthe relative ballistic inclination angle，δrthe relative angle of attack，δr1the projection of the relative angle of attack on the vertical plane，δr2the projection of the relative angle of attack on the lateral plane，ψrthe relative declination angle velocity，θthe trajectory inclination angle velocity，αrthe relative aiming angle，kxDthe polar damping moment coefficient，kzDthe equator torque damping coefficient，Mzthe Magnus moment，andMzthe stable moment．βDηandβDζare the projections of dynamic unbalance angleβDin axisηandζof theo1ξηζ，respectively，LmηandLmζthe projections of mass eccentricityLmin axisηandζof theo1ξηζ，respectively．

2 Determination Algorithm for Pro－jectile Attitude Angles

When computing the 6DOF trajectory model，one needs to provide 12independent initial values including （x，y，z，vx，vy，vz，γ，φa，φ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（γ，will be studied in this paper．

If the radar data at time instanttiandti＋1are known，one can expand the positions of a projectile using a second－order Taylor series attias

whereXiandXi＋1are the position coordinates of the projectile attiandti＋1，respectively；Viis the velocities of the projectile atti，Δt＝ti＋1－tithe difference betweentiandti＋1，Aiis the accelerations of the projectile atti．

Further，it can be deduced thatAi＝Moreover，Aican also be expressed as，where the（x，y，z，vx，vy，vz）iattiand （x，y，z，vx，vy，vz）i＋1atti＋1can 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 （γ，φa，φ2）iare unknown and independent，therefore，can be considered as a nonlinear function with respect toatti，namelyViΔt）．By solving，φa，φ2）i，the problem is transformed into solving a ternary nonlinear equations as follows

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 functionf（x），its gradient（f），and positive definite Hessian matrixB，the Taylor series is

and the Taylor series of the gradient itself（secant equation）

is used to updateB．The DFP formula finds a solution that is symmetric，positive definite and closest to the current approximate value ofBk

where

andBkis a symmetric and positive definite matrix．The corresponding update to the inverse Hessian approximationHk＝B－1kis given by

Bis assumed to be positive definite，and the vectorssTkandymust satisfy the curvature condition

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

We solve，φa，φ2）ibased on the DFP method．

In a similar way，known the radar data atti＋1andti＋2，one can calculate，φa，φ2）i＋1．Hereby，theattican be approximately evaluated by

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 attiti＋1andti＋2，namely （X，V）i，（X，V）i＋1，（X，V）i＋2；

Step 2 Use（X，V）iand（X，V）i＋1to calculate（γ，φa，φ2）iby DFP method；

Step 3 Use（X，V）i＋1and（X，V）i＋2to calculate（γ，φa，φ2）i＋1by DFP method；

Step 4 Use（φa，φ2）iand（φa，φ2）i＋1to calcu－late，）iby difference quotient．

Up to now，has been ob－tained．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 measure－ment，the initial valuesof exterior ballistic equations are all known．Further，the 6DOF trajectory model can be solved，and the impact point can be computed．The cycle stops whennis 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 155mm spin－stabilized projectile is studied for ballistic simulation．The mass of the projectile is 45．678kg；its equator moment of inertia is 1．785kg·m2，polar moment of inertia 0．164 3kg·m2， twist rate of rib rifling 0．141 260 16，length of rifling 6．41m．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 car－ried out based on three elevation angles including 35°，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 1to 2safter launch，the time step is set to be 0．01s．Thereby，total 100positions and velocity coordinates，as well as orientation angles and their time derivatives will be generated at the 100time 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 at－titude angles and their time derivatives．

Figs．3—11show 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 （6degree 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

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°．

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

Table 2 Four impact points of projectiles in the firing test

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

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 3into 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．12—16．

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

3．4 Test verification on the prediction of projec－

tile′s impact point

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

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 xversus time for comparison in the case of 45°elevation angle

Fig．15 Velocity yversus 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

Table 4shows the error analysis of the pre diction 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 theZdirection．From the meteorological conditions of the test，the high altitude wind speed，which can reach 62—77m／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．

The experiment validates that the algorithmhas good prediction effect using the 6DOF trajectory model．

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

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）．

［1］ RUI X T，LIU Y X，YU H L．Launch dynamics of tank and self－propelled artillery［M］．Beijing：Science Press，2008．（in Chinese）

［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．

［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．（in Chinese）

［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．

［8］ ERIC N，MEIR P，STANTON M．Projectile launch point estimation from radar measurement［C］∥ A－merican 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．

［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．（in Chinese）

［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．

［19］ZHAO S．Application of MATLAB programming and optimization design［M］．Beijing：Electronics Industry Publishing House，2013．（in Chinese）