Projectile Impact Point Prediction Based on Self-Propelled Artillery Dynamics and Doppler Radar Measurements

Abstract

Any trajectory calculation method has three primary sources of errors, which are model error, parameter error, and initial state error. In this paper, based on initial projectile flight trajectory data measured using Doppler radar system; a new iterative method is developed to estimate the projectile attitude and the corresponding impact point to improve the second shot hit probability. In order to estimate the projectile initial state, the launch dynamics model of practical 155 mm self-propelled artillery is defined, and hence, the vibration characteristics of the self-propelled artillery is obtained using the transfer matrix method of linear multibody system MSTMM. A discrete time transfer matrix DTTM-4DOF is developed using the modified point mass equations of motion to compute the projectile trajectory and set a direct algebraic relation between any two successive radar data. During iterations, adjustments to the repose angle are made until an agreement with acceptable tolerance occurs between the Doppler radar measurements and the estimated values. Simulated Doppler radar measurements are generated using the nonlinear six-degree-of-freedom trajectory model using the resulted initial disturbance. Results demonstrate that the data estimated using the proposed algorithm agrees well with the simulated Doppler radar data obtained numerically using the nonlinear six-degree-of-freedom model.

1. Introduction

There are different computation methods to estimate the projectile ballistic trajectory and subsequent impact point. These methods have three primary sources of errors, which are model error, parameter error, and initial state error. The model error depends on how much the projectile equations of motion are simplified. Increasing model complexity means decreasing model error, but the parameter and initial state errors will increase [1]. In 1964, the projectile six-degree-of-freedom (6-DOF) equations of motion in terms of direction cosines [2] was developed by the ballistic research laboratory, which is the most accurate model that can simulate the rigid body flight dynamics under the condition that all the aerodynamic forces and moments, and the initial state, are known to a high degree of accuracy. Using this system of equations, a discrete time transfer matrix method DTTM-6DOF [3] has been developed to accurately compute the projectile ballistic trajectory and minimize the corresponding computational time.

During testing new weapons, an accurate instrumentation is needed to measure the in-flight projectile parameters. Projectiles do not have enough space to install on-board sensors to acquire these parameters. Further, considering the severe environment faced by the projectile during flight such as velocity and spinning rate, a low cost and accurate measurement method is required. The most promising method is using a Doppler radar system [4] due to its robustness, portability, accuracy, and providing nearly continuous velocity-time information on the projectile during the flight down range. A number of authors presented different algorithms based on Doppler radar measurements. In 1958, Shapiro [5] proposed three different estimation methods to compute the six parameters that define the orbital parameters describing the ballistic missile trajectory using single radar measurements which include range, range-rate, and the azimuth and elevation angles. These methods are the maximum likelihood, the weighted least squares, and a deterministic method. In 1985, a new iterative method had been developed [6] and applied to a modified version of 6-DOF model to extract actual thrust made during the flight of rocket-assisted projectiles using Doppler radar measurements. In 2008, an iterative algorithm had been developed [7] using maximum likelihood method to predict the drag characteristics of projectiles in motion by processing data acquired by Doppler radar system, where the point mass model was used due to the projectile small angle of attack along the trajectory.

The movements of artillery and projectile are very complex in launch process because of the complex structure of self-propelled artillery and the severe environments, such as, high temperature, high pressure, high speed, instantaneous state, multibody, and mutation, in launch process [8]. Therefore, gun launch dynamics has a primary concern in designing embedded navigation and guidance systems in case of guided projectiles and evaluating the projectile initial disturbance in case of unguided projectiles. In 2000, the launch dynamics of self-propelled artillery [9] based on transfer matrix method of multibody system MSTMM was developed for studying the vibration characteristics and dynamics of self-propelled artillery and the projectile's initial disturbance. In 2001, a finite element analysis [10] had been done for the breech closure for the 155 mm Cannon M199, which is normally mounted on the Towed Howitzer M198. The analysis is for a 9 body problem with 13 contact surfaces and was solved for both static and dynamic load cases. In 2004, a simplified model [11] had been done to predict joint forces and accelerations along the length of a simplified 155 mm projectile using modal superposition. These predicted loads used to locate all guidance equipment, sensors, joints, and computers that must operate reliably after exiting the gun. In 2005, a sophisticated 3D-finite-element model [12] had been developed to investigate the survivability of embedded electrical systems by simulating the launch dynamics of a surrogate Excalibur projectile. This study determined free mounting locations trouble for sensitive components and parametric investigation in identifying sensitive factors affecting the muzzle-exit motion of projectile substructures. In 2008, a new method [13] named as transfer matrix method of linear multibody system MSTMM for linear hybrid multibody system dynamics was developed. This method had been applied to a Multiple Rocket System as a linear multi-rigid-flexible-body system. In 2011, the dynamics problem of a shipboard gun system had been solved [14] using the discrete time transfer matrix method of multibody system MSDTTMM [15].

In this paper, an iterative method is developed to determine the projectile kinematics including attitude and angular motion using Doppler radar measurements. Doppler radar measurements are only available during the first portion of projectile trajectory including the range, the range-rate, and the azimuth and elevation angles (position and velocity vectors), in order not to be recognized by the enemy counterattack systems. Simulated Doppler radar measurements are generated during flight time using 6-DOF trajectory model including the projectile initial disturbance problem. The projectile initial disturbance is computed by solving the self-propelled artillery and projectile dynamics as a multi-rigid-flexible system using Transfer Matrix Method of Multibody System MSTMM [8], where the overall transfer equation of self-propelled artillery was obtained to compute the corresponding vibration characteristics. In order to iteratively compute the projectile attitude, a discrete time transfer matrix DTTM-4DOF is developed based on linearization of the nonlinear modified point mass equations of motion using second-order Taylor series [16]. The system state vector, which defines the projectile trajectory parameters during Doppler radar observations, is defined by 13 state variables, which are the projectile position, velocity, spin rate, repose angle, and wind velocity. During iterations, adjustments to the repose angle are made until agreement with a certain tolerance occurs between the Doppler radar velocity measurements and the value predicted using DTTM-4DOF. Finally, projectile impact point can be predicted precisely depending on the length of radar data available, where the data estimated at the end of radar observations is used as an initial state to the 6-DOF trajectory model.

2. Launch Dynamics Model for Self-Propelled Artillery System

To accurately estimate the projectile initial state, launch process has to be decided precisely by studying the artillery and projectile dynamics from the moment of firing to the state of muzzle point. Four interactions have to be considered, the interaction between projectile and artillery, the interaction between the artillery transverse and longitudinal motion, the interaction between rigid and flexible body, and, finally, the interaction between space coordinate and time coordinate for describing motion parameters.

The launch dynamics model of self-propelled artillery is a multi-rigid-flexible system which has been studied in [8]. This system is composed of 51 elements (23 bodies and 28 joints) as shown in Figure 1. It has 3 boundary ends including ground, gun breech, and muzzle. The ground (0) is regarded as an infinity rigid body and the first input point, the gun breech (44) is regarded as rigid body and the second input point, and the muzzle (52) is the output point. The road wheels (13–24) are regarded as lumped masses. The hull (37), revolving part (39), elevating part (41), and muzzle brake (51) are regarded as rigid bodies; the barrel is divided into 6 segments (45–50) and each is regarded as a beam with equal sectional area. The interactions (1–12) between ground and road wheel and the connections (25–36) between road wheel and the hull are, respectively, modeled with springs and dampers in 3 directions. Connections (38, 40, 42, and 43) are modeled with springs and rotary springs accompanying dampers to represent relative linear motion and relative angular motion in 3 directions at the same time. All connections between barrel and gun breech, barrel and muzzle brake, and each segment of barrel are regarded as fixed. The masses of traversing mechanism, elevating mechanism, and equilibrator fall into elevating part, revolving parts, and hull, respectively. According to the launch dynamics model of self-propelled artillery and the sequence number of each element, there are 49 connection points which are P_{i, i – 12}, P_{i, i + 12} (i = 13–24), P_{37, i} (i = 25–36), P_37,38, P_39,38, P_39,40, P_41,40, P_41,42, P_41,43, P_44,45, P_45,46, P_46,47, P_47,48, P_48,49, P_49,50, and P_51,50 and 14 boundary points which are P_{0, i} (i = 1–12), P_44,0, P_51,52, respectively. The overall transfer equation of the practical 155 mm self-propelled artillery shown in Figure 1 was obtained using the transfer matrix method of linear multibody system MSTMM [8].

Figure 1:

Launch dynamic model of self-propelled artillery [8].

During launch process, all forces and moments acting between barrel and elevating part, elevating part and revolving part, revolving part and the hull, the hull and road wheels, and road wheels and the ground are regarded as internal forces. The forces that acted on the barrel by projectile [8, 17] and the forces that acted on the barrel and muzzle brake by propellant gas [8] are regarded as external forces of the artillery. The recoil resistance between the barrel and elevating part produced by recoil and counter-recoil mechanism is decomposed into external force and internal force, and the internal force is proportional to the relative displacement between barrel and elevating part.

3. Projectile Flight Trajectory Model

Due to the complexity of projectile initial disturbance determination, a modified point mass 4-DOF model for an unguided rigid projectile was stated in [18, 19]. These equations assume that the epicyclic pitching and yawing motion of the projectile are small everywhere along the trajectory; therefore, the yaw and pitch moments considered in 6-DOF trajectory model [19] are neglected.

The projectile modified point mass equations of motion with respect to earth fixed coordinate system (X₁, X₂, X₃) as shown in Figure 2 are given by [19]

$\dot{\mathit{X}}=\mathit{V}={\left[\begin{array}{ccc}{V}_1& V_2& V_3\end{array}\right]}^T,$

(1)

$\dot{\mathit{V}}=-\frac{\rho v{S}_{\text{ref}}{C}_D}{\mathrm{2}m}{\mathit{V}}_{\mathit{t}}+\frac{\rho S_{\text{ref}}{C}_L^{\alpha }}{\mathrm{2}m}{v}^{\mathrm{2}}{\mathit{\alpha }}_R\hspace{1em}+\frac{\rho S_{\text{ref}}d{C}_N^{p\alpha }}{\mathrm{2}m}p\left[{\mathit{V}}_{\mathit{t}}\times {\mathit{\alpha }}_R\right]+\mathbf{g}+\mathrm{\Lambda }\mathbf{,}$

(2)

$\dot{p}=\frac{\rho S_{\text{ref}}{d}^{2}v}{2I_x}p{C}_l^p,$

(3)

where C _D , C _L ^α, C _N ^pα, and C _l ^p are the aerodynamic drag force, lift force, Magnus force, and spin damping moment coefficients, respectively; V is the projectile velocity vector in earth fixed coordinate; V _t is the projectile total velocity vector with respect to air; v is the projectile total aerodynamic speed; ρ is the air density; m, I _x , I _y , d, and S_ref are the projectile mass, axial and transverse moments of inertia, reference diameter, and reference area, respectively; and g and Λ are the earth's gravitational and Coriolis accelerations vectors.

Figure 2:

Coordinate system and directions of ballistic target.

The projectile repose angle α _R is defined by [19]

${\mathit{\alpha }}_R=\mathit{x}-\left(\text{cos}\left|{\alpha }_t\right|\right)\mathit{i},$

(4)

where x is the projectile's axis of symmetry unit vector and i is the unit vector in the direction of V _t .

And the total aerodynamic angle of attack is [3, 18, 19]

${\alpha }_t={\text{cos}}^{-\mathrm{1}}\frac{{\mathit{V}}_{\mathit{t}}·\mathit{x}}{v}.$

(5)

The simplified form of the repose angle is given by [19]

${\mathit{\alpha }}_R=\frac{\mathrm{2}{I}_{x}p}{\rho S_{\text{ref}}d{v}^{\mathrm{4}}{C}_m^{\alpha }}\left[\dot{\mathit{V}}\times {\mathit{V}}_{\mathit{t}}\right],$

(6)

where C _m ^α is the aerodynamic pitching moment coefficient.

For discrete systems, a good approximation was done for the projectile repose angle [20] as

${\mathit{\alpha }}_R=\frac{\mathrm{2}{I}_{x}p}{\rho S_{\text{ref}}d{v}^{\mathrm{4}}{C}_m^{\alpha }}\left[\mathit{g}\times {\mathit{V}}_{\mathit{t}}\right].$

(7)

The Coriolis acceleration due to earth's rotation and gravitational acceleration are given by [3, 19]

$\begin{array}{c}\mathrm{\Lambda }=\mathrm{2}\Omega \left[\begin{array}{ccc}0& -\text{cos}\mu \text{sin}\lambda & -\text{sin}\mu \\ \text{cos}\mu \text{sin}\lambda & 0& \text{cos}\mu \text{cos}\lambda \\ \text{sin}\mu & -\text{cos}\mu \text{cos}\lambda & 0\end{array}\right]\left[\begin{array}{c}{V}_{\mathrm{1}}\\ V_{\mathrm{2}}\\ V_{\mathrm{3}}\end{array}\right],\\ \mathbf{g}=-g·{\left[\begin{array}{ccc}\frac{X_{\mathrm{1}}}{R_e}& 1-\frac{\mathrm{2}{X}_{\mathrm{2}}}{R_e}& 0\end{array}\right]}^T,\end{array}$

(8)

where g = g _o · [1 – 0.0026 cos(2μ)], g _o = 9.80665 m/s²; Ω is the earth's angular velocity; μ, λ are the corresponding latitude and longitude of the firing sight; and R _e is the average radius of the earth (−6356766 m).

Due to spherical earth approximation, the instantaneous projectile altitude H is defined as [19]

$H=X_{\mathrm{2}}+\frac{X_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{2}·R_e}.$

(9)

4. Discrete Time Transfer Matrix for Modified Point Mass Equations of Motion

To develop a discrete time transfer matrix DTTM-4DOF using the modified point mass equations of motion mentioned above, a system state vector S is defined by 13 state variables to define the projectile trajectory during flight time. These state variables are the projectile position, velocity, spin rate, repose angle, and wind velocity, and are as follows:

${\mathit{S}}_{\mathrm{13}\times \mathrm{1}}={\left[\begin{array}{ccccccccccccc}{X}_{\mathrm{1}}& X_{\mathrm{2}}& X_{\mathrm{3}}& V_{\mathrm{1}}& V_{\mathrm{2}}& V_{\mathrm{3}}& p& {\alpha }_{\mathrm{1}}& {\alpha }_{\mathrm{2}}& {\alpha }_{\mathrm{3}}& W_{\mathrm{1}}& W_{\mathrm{2}}& W_{\mathrm{3}}\end{array}\right]}^T.$

(10)

DTTM-4DOF is the relation between any two successive state vectors, where

${\mathbf{S}}_{\mathrm{10}\times \mathrm{1}}^{i\mathbf{+}\mathrm{1}}=\text{DTTM-}\mathrm{4}\text{DO}{\text{F}}_{\mathrm{10}\times \mathrm{13}}·{\mathbf{S}}_{\mathrm{13}\times \mathrm{1}}^i.$

(11)

The 2nd-order Taylor series is given by [16]

${\mathit{S}}_{\mathrm{10}\times \mathrm{1}}^{i+\mathrm{1}}={\mathit{S}}_{\mathrm{10}\times \mathrm{1}}^i+{\dot{\mathit{S}}}_{\mathrm{10}\times \mathrm{1}}^i·\mathrm{\Delta }t+{\stackrel{‥}{\mathit{S}}}_{\mathrm{10}\times \mathrm{1}}^i·\frac{\mathrm{\Delta }{t}^{\mathrm{2}}}{\mathrm{2}},$

(12)

where

${\dot{\mathbf{S}}}_{\mathrm{10}\times \mathrm{1}}·\mathrm{\Delta }t={\mathbf{I}}_{\mathrm{10}\times \mathrm{13}}{\mathbf{A}}_{\mathrm{13}\times \mathrm{13}}·{\mathbf{S}}_{\mathrm{13}\times \mathrm{1}},$

(13a)

${\stackrel{‥}{\mathbf{S}}}_{\mathrm{10}\times \mathrm{1}}^i·\frac{\mathrm{\Delta }{t}^{\mathrm{2}}}{\mathrm{2}}=\left[\frac{\partial {\dot{\mathbf{S}}}^i}{\partial t}·\frac{\mathrm{\Delta }{t}^{\mathrm{2}}}{\mathrm{2}}\right]=\left[{\left(\frac{\partial \dot{\mathbf{S}}}{\partial \mathbf{S}}\right)}_{\mathrm{10}\times \mathrm{13}}^i·{\dot{\mathbf{S}}}_{\mathrm{13}\times \mathrm{1}}^i·\frac{\mathrm{\Delta }{t}^{\mathrm{2}}}{\mathrm{2}}\right]=\left[{\mathbf{T}}_{\mathrm{10}\times \mathrm{10}}{\mathbf{J}}_{\mathrm{10}\times \mathrm{13}}{\mathbf{A}}_{\mathrm{13}\times \mathrm{13}}{\mathbf{S}}_{\mathrm{13}\times \mathrm{1}}\right].$

(13b)

The equations of motion matrix A, summarizing (1)–(9) multiplied by the time step Δt, are given by

$\mathit{A}=\left[\begin{array}{ccccc}{\mathit{G}}_{3\times 3}\frac{{\mathrm{\Delta }t}^2}{2}& {\mathit{I}}_3·\mathrm{\Delta }t& {\mathit{O}}_{3\times 1}& {\mathit{O}}_{3\times 3}& {\mathit{O}}_{3\times 3}\\ {\mathit{G}}_{3\times 3}·\mathrm{\Delta }t& {\dot{\mathit{V}}}_{\mathit{V}}·\mathrm{\Delta }t& {\dot{\mathit{V}}}_p·\mathrm{\Delta }t& {\dot{\mathit{V}}}_{\mathit{\alpha }}·\mathrm{\Delta }t& {\dot{\mathit{V}}}_{\mathit{W}}·\mathrm{\Delta }t\\ {\mathit{O}}_{1\times 3}& {\mathit{O}}_{1\times 3}& {\dot{p}}_p·\mathrm{\Delta }t& {\mathit{O}}_{1\times 3}& {\mathit{O}}_{1\times 3}\\ {\mathit{\alpha }}_{\dot{\mathit{V}}}{\mathit{G}}_{3\times 3}& {\mathit{\alpha }}_{\dot{\mathit{V}}}{\dot{\mathit{V}}}_{\mathit{V}}& {\mathit{\alpha }}_{\dot{\mathit{V}}}{\dot{\mathit{V}}}_p& {\mathit{\alpha }}_{\dot{\mathit{V}}}{\dot{\mathit{V}}}_{\mathit{\alpha }}& {\mathit{\alpha }}_{\dot{\mathit{V}}}{\dot{\mathit{V}}}_{\mathit{W}}\\ {\mathit{O}}_{3\times 3}& {\dot{\mathit{W}}}_{\mathit{V}}·\mathrm{\Delta }t& {\mathit{O}}_{3\times 3}& {\mathit{O}}_{3\times 3}& {\mathit{O}}_{3\times 3}\end{array}\right],$

(14)

where I₃ is the 3 × 3 identity matrix.

All elements of matrix A in (14) are stated as follows:

${\mathit{G}}_{\mathrm{3}\times \mathrm{3}}=\left[\begin{array}{ccc}-\stackrel{-}{g}& 0& 0\\ 0& \mathrm{2}\stackrel{-}{g}-g_{X_{\mathrm{2}}}& 0\\ 0& 0& 0\end{array}\right],$

(15)

where g _{X 2} = g/X₂, $\stackrel{-}{g}=g/R_e$ ,

${\dot{\mathit{V}}}_{\mathit{V}}={\mathrm{\Lambda }}_{\mathit{V}}-\underset{{\dot{\mathit{V}}}_{\mathit{W}}}{\underbrace{\rho a_1{C}_{D}v{\mathit{I}}_3}},$

(16)

where ${\mathrm{\Lambda }}_{\mathit{V}}=\mathrm{2}·\Omega ·\left[\begin{array}{ccc}0& -\cos \mu ·\sin \lambda & -\sin \mu \\ \cos \mu ·\sin \lambda & 0& \cos \mu ·\cos \lambda \\ \sin \mu & -\cos \mu ·\cos \lambda & 0\end{array}\right]$ ,

$\begin{array}{c}{\dot{\mathit{V}}}_p=\rho a_2{C}_N^{p\alpha }{\mathit{B}}_1,\\ {\dot{\mathit{V}}}_{\mathit{\alpha }}=\rho a_1{C}_L^{\alpha }{v}^2{\mathit{I}}_3,\\ {\dot{p}}_p=\rho b_{1}v{C}_l^p,\\ {\mathit{\alpha }}_{\dot{\mathit{V}}}=\frac{-b_{2}p}{\rho v^4{C}_m^{\alpha }}{\left[{\mathit{V}}_{\mathit{t}}\right]}_{\times },\end{array}$

(17)

where a₁ = (S_ref/2m), a₂ = ((S_refd)/2m), b₁ = ((S_ref · d²)/(2 · I _x )), b₂ = (d/b₁), and anjan65.

The change in wind velocity vector with the projectile position ${\dot{\mathit{W}}}_{\mathit{V}}$ is given by

${\dot{\mathit{W}}}_{\mathit{V}}=\frac{\partial \mathit{W}/\partial t}{\partial \mathit{V}}=\frac{\partial \mathit{W}/\partial t}{\partial \mathit{X}/\partial t}=\frac{\partial \mathit{W}}{\partial \mathit{X}}.$

(18)

The Jacobian matrix J of the system, mentioned in (13b), which represents discretization errors in matrix A, is given by

$\mathbf{J}=\left[\begin{array}{ccccc}{\mathbf{O}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{I}}_{\mathrm{3}}& {\mathbf{O}}_{\mathrm{3}\times \mathrm{1}}& {\mathbf{O}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{O}}_{\mathrm{3}\times \mathrm{3}}\\ {\mathbf{O}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{j}}_{\mathrm{22}}& {\dot{\mathit{V}}}_p& {\mathbf{O}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{j}}_{\mathrm{25}}\\ {\mathbf{O}}_{\mathrm{1}\times \mathrm{3}}& {\mathbf{j}}_{\mathrm{32}}& {\dot{p}}_p& {\mathbf{O}}_{\mathrm{1}\times \mathrm{3}}& -{\mathbf{j}}_{\mathrm{32}}\\ {\mathbf{O}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{j}}_{\mathrm{42}}& {\mathbf{j}}_{\mathrm{43}}& {\mathbf{O}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{j}}_{\mathrm{45}}\end{array}\right].$

(19)

The partial derivatives of the projectile acceleration in (2) with respect to wind velocity and projectile velocity are given by

${\mathit{j}}_{25}=\rho a_1{C}_{D}v{\mathit{I}}_3+[\phantom{\frac{\rho a_{2}p{C}_{N_M}^{p\alpha }M}{v^2}}\rho a_1(\left(C_D+C_{D_M}M\right)\frac{{\mathit{V}}_{\mathit{t}}}{v}-\phantom{\frac{{\mathit{V}}_{\mathit{t}}}{v}}\left(2C_L^{\alpha }+C_{L_M}^{\alpha }M\right){\mathit{\alpha }}_R\phantom{\frac{V_t}{v}})-\frac{\rho a_{2}p{C}_{N_M}^{p\alpha }M}{v^2}{\mathit{B}}_1]\hspace{0.17em}{\mathit{V}}_{\mathit{t}}^T+\rho a_{2}p{C}_N^{p\alpha }{\left[{\mathit{\alpha }}_R\right]}_{\times },$

(20)

where ${C_i}_i=\partial C_i/\partial M$ is the rate of change of any aerodynamic coefficient with respect to Mach number

${\mathit{j}}_{22}={\mathrm{\Lambda }}_{\mathit{V}}-{\mathit{j}}_{25}.$

(21)

Similarly, using (3), the rate of change of projectile spin rate with respect to its velocity is given by

${\mathit{j}}_{32}=\frac{\rho b_{1}p}{v}\left(C_l^p+C_{l_M}^{p}M\right){\mathit{V}}_{\mathit{t}}.$

(22)

The partial derivatives of the projectile yaw of repose in (7) with respect to projectile velocity and spin rate are as follows:

$\begin{array}{c}{\mathit{j}}_{42}=\frac{b_{2}p}{\rho v^4{C}_m^{\alpha }}{\left[\mathit{g}\right]}_{\times }+\frac{b_{2}p\left(4C_m^{\alpha }+C_{m_M}^{\alpha }M\right)}{\rho v^6{\left(C_m^{\alpha }\right)}^2}{\mathit{B}}_2{\mathit{V}}_{\mathit{t}}^T,\\ {\mathit{j}}_{43}=\frac{b_2}{\rho v^4{C}_m^{\alpha }}{\mathit{B}}_2,\end{array}$

(23)

where B₂ = g × V _t .

The time matrix T mentioned in (13b) is given by

${\mathit{T}}_{10\times 10}=\left[\begin{array}{cccc}{\mathit{I}}_3\frac{\mathrm{\Delta }t}{2}.& {\mathit{O}}_{3\times 3}& {\mathit{O}}_{3\times 1}& {\mathit{O}}_{3\times 3}\\ {\mathit{O}}_{3\times 3}& {\mathit{I}}_3.\frac{\mathrm{\Delta }t}{2}& {\mathit{O}}_{3\times 1}& {\mathit{O}}_{3\times 3}\\ {\mathit{O}}_{1\times 3}& {\mathit{O}}_{1\times 3}& \frac{\mathrm{\Delta }t}{2}& {\mathit{O}}_{1\times 3}\\ {\mathit{O}}_{3\times 3}& {\mathit{O}}_{3\times 3}& {\mathit{O}}_{3\times 1}& {\mathit{I}}_3\end{array}\right].$

(24)

By Applying (14), (19), and (25) into (12), the overall discrete time transfer matrix is given by

${\mathit{S}}_{\mathrm{10}\times \mathrm{1}}^{i+\mathrm{1}}=\left[{\mathbf{M}}_{\mathrm{10}\times \mathrm{13}}+\left[{\mathbf{I}}_{\mathrm{10}\times \mathrm{13}}+{\mathbf{T}}_{\mathrm{10}\times \mathrm{10}}·{\mathbf{J}}_{\mathrm{10}\times \mathrm{13}}\right]·{\mathbf{A}}_{\mathrm{13}\times \mathrm{13}}\right]·{\mathit{S}}_{\mathrm{13}\times \mathrm{1}}^i,$

(25)

where, ${\mathbf{M}}_{\mathrm{10}\times \mathrm{13}}=\left[\begin{array}{cc}{\mathbf{I}}_{\mathrm{7}\times \mathrm{7}}& {\mathbf{O}}_{\mathrm{7}\times \mathrm{6}}\\ {\mathbf{O}}_{\mathrm{3}\times \mathrm{7}}& {\mathbf{O}}_{\mathrm{3}\times \mathrm{6}}\end{array}\right]$ ,

$\text{DTTM-}\mathrm{4}\text{DOF}=\left[\begin{array}{ccccc}{\mathbf{m}}_{\mathrm{11}}& {\mathbf{m}}_{\mathrm{12}}& {\mathbf{m}}_{\mathrm{13}}& {\mathbf{m}}_{\mathrm{14}}& {\mathbf{m}}_{\mathrm{15}}\\ {\mathbf{m}}_{\mathrm{21}}& {\mathbf{m}}_{\mathrm{22}}& {\mathbf{m}}_{\mathrm{23}}& {\mathbf{m}}_{\mathrm{24}}& {\mathbf{m}}_{\mathrm{25}}\\ {\mathbf{m}}_{\mathrm{31}}& {\mathbf{m}}_{\mathrm{32}}& {\mathbf{m}}_{\mathrm{33}}& {\mathbf{m}}_{\mathrm{34}}& {\mathbf{m}}_{\mathrm{35}}\\ {\mathbf{m}}_{\mathrm{41}}& {\mathbf{m}}_{\mathrm{42}}& {\mathbf{m}}_{\mathrm{43}}& {\mathbf{m}}_{\mathrm{44}}& {\mathbf{m}}_{\mathrm{45}}\end{array}\right],$

(26)

${\mathit{S}}_{\mathrm{10}\times \mathrm{1}}^{i+\mathrm{1}}=\left[\text{DTTM-}\mathrm{4}\text{DOF}\right]·{\mathit{S}}_{\mathrm{13}\times \mathrm{1}}^i.$

(27)

All elements m _ij of the matrix DTTM-4DOF are listed in the appendix.

5. Projectile Attitude Determination Algorithm

The algorithm presented here is an iterative method to compute (1) the projectile muzzle velocity and line of sight by using the first projectile displacement observed by radar, and (2) the projectile flight attitude and the corresponding angular momentum by using the Doppler velocity observed by radar. At each time step, the difference between the data measured by radar and the estimated parameters by using DTTM-4DOF are linearly related as follows:

${\left[\begin{array}{c}\mathrm{\Delta }\mathit{X}\\ \mathrm{\Delta }\mathit{V}\end{array}\right]}_{i+1}=\left[\begin{array}{cccc}{\mathit{m}}_{11}& {\mathit{m}}_{12}& {\mathit{m}}_{13}& {\mathit{m}}_{14}\\ {\mathit{m}}_{21}& {\mathit{m}}_{22}& {\mathit{m}}_{23}& {\mathit{m}}_{24}\end{array}\right].{\left[\begin{array}{c}\mathrm{\Delta }\mathit{X}\\ \mathrm{\Delta }\mathit{V}\\ \mathrm{\Delta }p\\ \mathrm{\Delta }{\mathit{\alpha }}_R\end{array}\right]}_i.$

(28)

At the beginning of flight time, the initial guess of the projectile velocity V _o is iteratively corrected by adding ΔV _o using the first radar data available (i = 1) until the position difference ΔX₁ between the estimated and observed value is small enough (less than a given tolerance) as follows:

$\mathrm{\Delta }{\mathit{V}}_o={\left[{\mathit{m}}_{12}\right]}^{-1}\mathrm{\Delta }{\mathit{X}}_{\left(i=1\right)}.$

(29a)

The projectile muzzle spin rate is corrected by adding Δp _o using the velocity difference ΔV _o as

$\mathrm{\Delta }{p}_o=\frac{2\pi ·\left|{\mathit{V}}_o\right|-\left|{\mathit{V}}_o+\mathrm{\Delta }{\mathit{V}}_o\right|}{n·d},$

(29b)

where n is the gun twist rate.

At the beginning of each iteration cycle, the DTTM-4DOF matrix is computed and hence the projectile states ${\mathit{S}}_{\mathrm{10}\times \mathrm{1}}^{i+\mathrm{1}}$ in (27) to obtain the projectile velocity difference ΔV_{i + 1} which will be reduced until reaching a predefined tolerance by adding $\mathrm{\Delta }{{\mathit{\alpha }}_R}_i$ to the initial state as

$\mathrm{\Delta }{{\mathit{\alpha }}_R}_i={\left[{\mathit{m}}_{\mathit{24}}\right]}^{-1}\mathrm{\Delta }{\mathit{V}}_{i+1},$

(30a)

${{\mathit{\alpha }}_R}_i={{\mathit{\alpha }}_R}_i+\mathrm{\Delta }{{\mathit{\alpha }}_R}_i.$

(30b)

Then, the estimated parameters are used as initials for the next time step.

As shown in Figure 3, during each time step, the projectile attitude and the corresponding angular momentum can be estimated using the following procedure.

1. The projectile total angle of attack can be computed by

   ${\alpha }_t=\left|{\mathit{\alpha }}_R\right|.$

   (31)

2. The projectile's axis of symmetry unit vector x is computed by [19]

   $\mathit{x}={\mathit{\alpha }}_R+\left(\text{cos}{\alpha }_t\right)\mathit{i},$

   (32)

   where i = V _t /|V _t |.

3. The rate of change of the projectile total velocity unit vector is given by [19]

   $\dot{\mathit{i}}=\frac{\rho S_{\text{ref}}v{C}_L^{\alpha }}{2m}{\mathit{\alpha }}_R+\frac{1}{v}\left[\mathit{i}\times \left(\mathit{g}\times \mathit{i}\right)\right].$

   (33)

4. Due to the small variation in projectile repose angle during flight [19], its rate of change is approximated to be linear during each time step as

   ${\dot{\mathit{\alpha }}}_R\approx \frac{{\mathit{\alpha }}_{R_{i+\mathrm{1}}}-{\mathit{\alpha }}_{R_i}}{\mathrm{\Delta }t}.$

   (34)

5. The rate of change of projectile's axis of symmetry unit vector $\dot{\mathit{x}}$ is computed by [19]

   $\dot{\mathit{x}}={\dot{\mathit{\alpha }}}_R+\left(\text{cos}{\alpha }_t\right)\dot{\mathit{i}}-\left(\text{sin}{\alpha }_t\right)\mathit{i}.$

   (35)

6. The projectile angular momentum divided by I _y in earth fixed coordinate is given by [19]

   $\mathit{h}=\left(\frac{I_x·p}{I_y}\right)\mathit{x}+\left[\mathit{x}\times \dot{\mathit{x}}\right].$

   (36)

Figure 3:

Projectile attitude and impact point determination flow chart.

6. Simulated Trajectory Generation

In order to validate the proposed algorithm, a simulated noiseless ballistic trajectory for 155 mm high explosive HE spin stabilized projectile is generated to simulate the Doppler radar velocity measurements data. The projectile mass, length, center of mass measured from the nose, and the axial and transverse mass moments of inertia are 46.5 kg, 0.902 m, 0.593 m, 0.1585 kg/m², and 1.8816 kg/m², respectively. This trajectory is generated using the nonlinear 6-DOF trajectory model [2, 19] including projectile initial disturbance problem with 52° firing elevating angle.

6.1. Natural Vibration Characteristics of the Self-Propelled Artillery

To compute the vibration characteristics for 155 mm self-propelled artillery, the input points corresponding to sequence numbers of connection points and elements shown in Figure 1 are identified in Table 1.

Table 1:

Corresponding relation between the input point and the sequence number.

Input point	Sequence no.
1	Ground
2–13	P13,1, P14,2, …, P24,12
14	P 37,25
15	P 37,38
16	P 39,38
17	P 39,40
18	P 41,40
19	P 41,42
20	P 44,0
21	P 44,45
22–29	45 (beam 1)
30	P 45,46
31–68	46 (beam 2)
69	P 46,47
70–104	47 (beam 3)
105–126	48 (beam 4)
127–148	49 (beam 5)
149–229	50 (beam 6)
230	P 51,52

All simulation and test results of natural frequencies for the first sixteen ranks are listed in [8], where the simulation results have a good agreement with that obtained by experiments. Figures 4 and 5 show part of the resulted system mode shapes, where the 1st mode mainly represents the recoil motion of barrel, the 2nd mode represents the horizontal revolving motions of turret, cradle, and barrel, the 3rd mode mainly represents the elevating motions of cradle and barrel, and the 16th mode mainly represents the transverse motion of barrel and complex motions of other elements. The low rank mode shapes do not represent the transverse elastic deformation of the barrel that is only represented in high rank [8].

Figure 4:

The 1st and 2nd mode shapes in X- and Z-direction, respectively.

Figure 5:

The 3rd and 16th mode shape in Y-direction.

6.2. Computation of Projectile Initial Disturbance

The resulted projectile launch dynamics parameters with respect to gun tube coordinate system are shown in Figures 6–10. The gun tube coordinate system O₃x _O y _O z _O is defined as follows: the origin O₃ is the intersection point between axis of gun tube and the vertical plane of projectile's axis of symmetry, the x _O -axis is along gun tube axis of symmetry and points to muzzle, the y _O -axis (plumb direction) is vertical to x _O -axis and points up, the direction of z _O -axis (swing direction) is determined by the right-hand role. The projectile longitudinal displacement, velocity, and acceleration with respect to gun tube axis are shown in Figures 6 and 7, where L _b is the total barrel length. Figure 8 shows the projectile center of gravity plumb and sidewise displacement with respect to gun tube axis as function of nondimensionalized longitudinal displacement.

Figure 6:

Projectile longitudinal displacement.

Figure 7:

Projectile longitudinal velocity and acceleration.

Figure 8:

Projectile lateral displacement.

Figure 9:

Projectile angular displacement during launch.

Figure 10:

Projectile angular velocity during launch.

The projectile axes orientation with respect to line of sight LOS, which includes the projectile roll, plumb and sidewise swing components, angular displacement, and velocity, are illustrated in Figures 9 and 10, respectively. The corresponding projectile angular displacement and velocity components at muzzle are shown in Table 2, where the deflection angle is the angle between the projectile velocity vector and LOS while projectile is just flying off from muzzle, and the swing angle is the angle between projectile axes and LOS at the same moment.

Table 2:

Simulation results of initial disturbance of projectile.

Initial disturbance	Plumb	Sidewise
Deflection angle/(10−3 rad)	0.71131	1.4600
Swing angle/(10−3 rad)	2.11265	−0.06758
Swing angle velocity/(rad/s)	−2.92724	−3.31593

The projectile initial flight conditions corresponding to the self-propelled artillery launch dynamics, and hence, the projectile initial disturbance calculations proposed before are listed in Table 3, which include the projectile initial velocity V₁, V₂, V₃, initial roll p, pitch q and yaw r rates, initial angle of attack α, and side-slip angle β.

Table 3:

Simulation results of the projectile parameters at muzzle.

V1 (m/s)	V2 (m/s)	V3 (m/s)	α (10−3 rad)	β (10−3 rad)	p (rad/s)	q (rad/s)	r (rad/s)
429.3847	551.9871	1.021021	1.40134	−1.52758	1279.72	−2.92724	−3.15983

7. Estimation Results

A computer program has been developed by using the nonlinear 6-DOF dynamic equations stated in [19]. All simulations were computed with fixed time step size of 0.001 s and integrated using Runge-Kutta-Gill method.

Firstly, a simulated noiseless Doppler radar velocity measurements data 6DOF-R are generated, including the projectile's position and velocity vectors, using the initial flight conditions listed in Table 3. The Doppler radar output data frequency is assumed to be 100 Hz. The 6DOF-R data is assumed to be the real projectile trajectory data. The corresponding impact point range, drift, maximum altitude, velocity, and flight time are 19689.2 m, 1027.85 m, 8180.9 m, 338.76 m/s, and 81.932 s, respectively.

Secondly, a nondisturbed trajectory simulation has been performed as nominal trajectory 6DOF-N, where the initial flight conditions used are 700 m/s, 52°, and 1279.7 rad/s for muzzle velocity, elevating angle, and spin rate, respectively. The corresponding impact point range, drift, maximum altitude, velocity, and flight time are 19729 m, 995 m, 8162 m, 338.66 m/s, and 81.829 s, respectively.

Finally, the iterative algorithm shown in Figure 3 has been applied to estimate the projectile attitude, using the generated Doppler radar data 6DOF-R, where the initial flight conditions are as the data used to simulate the projectile nominal trajectory 6DOF-N. Figures 11, 12, and 13 show that the estimated projectile velocity, AOA, and attitude, respectively, are matching very well with the simulated Doppler radar trajectory 6DOF-R, although the nominal (non-disturbed) trajectory results 6DOF-N are mismatching the 6DOF-R. Figures 14 and 15 show the rate of change of the projectile AOA and attitude, respectively, during flight time, where the estimated values are matched very well with 6DOF-R except at some peaks/dips data points due to the large difference between the data frequency obtained from radar and the simulation speed of the real trajectory data (1: 10). Therefore, the final estimated value shown in Figure 16, which represents the projectile angular momentum per-unit I _y has mismatched with 6DOF-R data during these peaks/dips data points.

Figure 11:

Projectile velocity vector versus flight time.

Figure 12:

Projectile total angle of attack vector versus flight time.

Figure 13:

The projectile's axis of symmetry unit vector.

Figure 14:

Rate of change of the projectile AOA.

Figure 15:

Rate of change of the projectile's axis of symmetry unit vector.

Figure 16:

Projectile angular momentum per unit I _y .

For nominal trajectory data 6DOF-N, the impact range ∊ _R and drift ∊ _D errors relative to 6DOF-R are 0.202% and −3.2%, respectively. All impact point parameters inaccuracies for the proposed algorithm as function of simulated Doppler radar data availability are shown in Figure 17, where these errors, ∊ _R and ∊ _D , can be reduced by approximately 8 and 7 times, respectively, as the radar data is available for the first 10% of total flight time and approximately by 10 and 70 times to be −0.0205% and −0.0448%, respectively, as the radar data availability is 15% of total flight time.

Figure 17:

Impact range and drift error versus Doppler Radar data availability.

8. Conclusion

The main objective of this paper is determining the effect of self-propelled artillery dynamics on the projectile impact point prediction accuracy, and hence, an iterative projectile attitude determination algorithm was developed to accurately compute the projectile ballistic trajectory and the corresponding impact point parameters by using a Doppler radar data. The projectile initial disturbance was computed by solving the interaction between the self-propelled artillery and the projectile dynamics during launch process. A discrete time transfer matrix DTTM-4DOF was developed by linearizing the projectile nonlinear modified point mass equations using second-order Taylor series. The main advantages of using DTTM-4DOF are (1) linearly relating any two successive radar data, and (2) increasing the computation speed due to solving a set of algebraic equations instead of using traditional numerical methods. An iterative projectile attitude determination algorithm was developed to accurately compute the projectile ballistic trajectory and predict the impact point parameters by using a Doppler radar data. In order to validate the proposed algorithm, a simulated noiseless trajectory data was generated as Doppler radar velocity measurements data 6DOF-R for 155 mm HE spin stabilized projectile using the nonlinear six-degree-of-freedom trajectory model including projectile initial disturbance problem. A nominal ballistic trajectory 6DOF-N was generated with 52° firing elevating angle and without including the effect of launch dynamics to show the importance of the projectile initial disturbance. Based on the results obtained, the proposed algorithm estimated the projectile velocity and attitude very well compared to the 6DOF-R trajectory data. The impact range and drift inaccuracies were reduced from 0.202% and −3.2% in case of nominal trajectory 6DOF-N to be −0.0205% and −0.0448% by applying the proposed algorithm as the radar data availability is 15% of total flight time. As the projectile dynamic stability increased, the radar data availability can be decreased. Finally, the proposed algorithm can estimate the projectile impact point accurately without computing the projectile initial disturbance, where different projectiles have different initial disturbances based on many factors, such as, gun barrel wear and erosion, projectile initial position and attitude inside barrel, and the propellant charge amount, composition, and homogeneity.