Ballistic Target Tracking with Use of Cinetheodolites

The paper addresses a problem of ballistic object tracking with the use of the cinetheodolite electro-optical tracking system. Electrooptical systems are applied for acquiring the trajectory data of missiles, satellites, and rockets used for delivery of satellites to their prevised orbits. Despite the importance of such systems and their applications, in the open literature there are no publications describing tracking algorithms processing data from cinetheodolites. The paper describes a model-based algorithm of estimation of position and parameters of target motion for such a system, developed by the authors. The model of the system, nonlinear both in its description of the target dynamics and the measurement equations, is presented in detail. The proposed algorithm of estimation is also described, and chosen simulation results are included in the paper. Furthermore, a comparison of the proposed estimation algorithm with other possible, but simpler algorithms is presented.


Introduction
The main motivation for the work presented in this paper was a necessity to obtain possibly accurate and continuous data on position and other parameters of motion of ballistic objects in practical tests of weaponry conducted by the Military Institute of Armament Technology (MIAT), Zielonka, Poland.Among the areas of expertise of the Institute, tests of anti-aircraft missiles are of great importance.One of the methods to assess the capabilities of the missile is the live firing with trajectory recording.It allows estimating the position and velocity of the missile and the target and eventually the mutual geometrical relationship between them, particularly the miss distance.
To estimate the position of a missile or a target, either an accurate on-board positioning system or an external tracking system is necessary.To assess the miss distance, an aerial target imitator fulfilling demanding requirements, between the others in terms of velocity, must be used.At MIAT, a special aerial target, ICP-89, based on the S-5 unguided missile, was developed.Its feasibility of deploying, high initial velocity, and market availability make it adequate for rocket missile tests.ICP-89 has proved its usefulness in numerous experiments with weapons so far.
In practice, an accurate on-board positioning system with an appropriate downlink is unlikely to be installed on a combat missile (particularly a small one) or an aerial target such as ICP-89.Therefore, to assess the performance of the missile during live tests, electro-optical tracking systems are widely used.A cinetheodolite-based electro-optical tracking system (EOTS), acquiring precise angular data to assess flight paths of disparate objects, is currently used at MIAT.It is a core system to determine the performance of the tested missiles.
In the open literature, there are no publications describing tracking algorithms for processing EOTS data; therefore, custom solutions were elaborated at MIAT.So far, straightforward geometrical computations have been used for position calculations of the missile and the aerial target.Such a method, however, has several significant drawbacks.Incomplete measurement data, occasionally occurring during some experiments, require supplementary techniques to attain only approximate and reduced-accuracy results.Moreover, the applied algorithm is not model-based and does not involve any data filtering; thus, its accuracy is not optimal.In this paper, a new algorithm of estimation of the position and parameters of ICP-89 target motion, based on the extended Kalman filter (EKF) [1][2][3][4], is proposed.It improves the accuracy of the estimated position, provides additional parameters of the target motion, such as its velocity and angular orientation, which are not directly available from the previously applied algorithm, and ensures availability of highly accurate, continuous data even in case of temporary outages of the measurements.
The layout of the further part of this paper is as follows.The EOTS measurement system is presented in Section 2. The proposed algorithms of estimation of a ballistic target position and parameters of motion together with a complete state-space model of the tracking system are given in Section 3. The paper includes also chosen results of simulations, demonstrating accuracy of the proposed algorithms (in Section 4) as well as discussion of the obtained results and conclusions given in Section 5.

EOTS Measurement System
The EOTS measurement system used for fire tests consists of at least two cinetheodolite stations, located at positions measured before the experiments, and providing for azimuth α and elevation λ angles of the tracked object.A photography of such a station is shown in Figure 1.
Positions of the stations and tracked object as well as all the parameters of motion are expressed in a local horizontal frame of reference which must be defined at the beginning of the experiments.The geometrical relationships in the EOTS system with two cinetheodolite stations S 1 and S 2 are explained in Figure 2 [5].
The frame of reference is established as follows.Firstly, the positions of S 1 and S 2 stations are measured with the use of an accurate geodetic-class GNSS receiver.These coordinates are obtained in a geodetic frame of reference [4].At the firing stand, about one second before the ICP-89 is launched, its tracer is ignited, and it burns for the whole period of flight of the imitator.Therefore, the ICP-89 is clearly visible from the cinetheodolite stations and the EOTS can establish the angles of its observation.Based on the measured azimuth and elevation angles, its position expressed in the geodetic frame is also calculated.The local frame of reference is assumed to begin at point O, corresponding to the calculated ICP-89 position, and its local coordinates are set to zero.The orientation of the X-axis is assumed to correspond to the direction of shooting, which is approximated by a projection of the longitudinal axis of the launcher on the horizontal plane.The Z-axis is the local vertical, and the Y-axis completes a right-handed Cartesian frame of reference.
After such a setting of the local frame of reference, the local coordinates of S 1 and S 2 stations are calculated from the already known geodetic coordinates of O, S 1 , and S 2 .Local coordinates of these points are as follows: O = 0, 0, 0 , S 1 = x 1 , y 1 , z 1 , and S 2 = x 2 , y 2 , z 2 .For simplicity of further calculations, S 1 and S 2 are assumed to be located on the same horizontal plane as O; thus, their coordinates are S 1 = x 1 , y 1 , 0 and S 2 = x 2 , y 2 , 0 .This assumption is not very restrictive, and the equations presented throughout the paper can be easily extended to a more general case of nonzero altitudes of S 1 and S 2 .
The cinetheodolite stations are synchronized by a common registration-triggering signal and perform acquisition of the images of the observed object as well as the angles of its observation from separate locations.The angles provided by two stations include azimuths α 1 and α 2 along with the elevations λ 1 and λ 2 , which are synchronized and timestamped.The angle measurements are calculated from readouts obtained from encoders installed in the stations.The encoders register only changes of respective angles from an arbitrary point of reset.To obtain meaningful angular data, a known reference point must be recorded.
A pair of measurements of azimuth α and elevation λ of the tracked object, obtained in two cinetheodolite stations of known locations, is sufficient for determination of the object positions and reconstruction of its trajectory.From the geometry of the system shown in Figure 2, the following formulae can be determined: Figure 2: Geometrical relationships in the EOTS measurement system [5].

2
International Journal of Aerospace Engineering where d is the distance between both EOTS stations (base of the system) and d 1 and d 2 are horizontal distances from S 1 and S 2 stations to the tracked object.Based on the above equations, the following set of four equations with three unknowns x, y, z can be formulated: This is an overdetermined set of equations, which can be approximately solved for x, y, z , if the sum of azimuth angles α 1 + α 2 ≠ 0. The requirement can be easily fulfilled by appropriately locating the cinetheodolite stations behind the ICP-89 launcher, so that the tracked object would not pass over the baseline S 1 S 2 of the system.Equation ( 7) can be simplified, if we notice that in the assumed local-level frame of reference, the target moves only in the OXZ plane.Thus, to represent its position we need only two coordinates x, z and the third coordinate y is always 0. Such an assumption was used for elaboration of the models and algorithms presented in the paper, and it resulted in the following simplified version of Equation (7): In practice, disturbing factors and an inaccurate setting of the reference frame may result in a nonzero value of y, which can be considered a cross-track error.However, it does not significantly affect the estimation results for the along-track position x and altitude z.Thus, the results and conclusions presented in the paper would remain valid also for a more general case of modelling motion in all three axes of OXYZ.

Trajectory Estimation of Target Imitator
Estimation of the ICP-89 position can be realized with the use of various algorithms, which can be generally divided into two groups, i.e., non-model-based algorithms and modelbased ones.The simplest are single-point solution algorithms, such as ordinary-least-squares (OLS) and weighted-leastsquares (WLS) methods, which iteratively estimate the target position [6][7][8][9], individually for each new set of measurements.The model-based estimation methods make use of the previous estimates and project them ahead via model equations up to the time when the new measurements have been acquired.Then, the projected estimates are fused with the new measurements to provide for a more accurate solution.An example of such an algorithm is the EKF used by the authors for the ICP-89 trajectory estimation and described in this paper.
3.1.OLS Algorithm.The primary measurements realized in the EOTS measurement system are azimuth and elevation angles, which are obtained simultaneously, with a period of about 33 msec.For every time instance, they can be grouped into the primary measurement vector z ′ as follows: As the relationship between the original measurement vector z ′ and the unknown target position x, z in OXZ is nontrivial, we propose here to construct a new measurement vector z as follows: where t * is a transformation function of the primary (original) measurement vector z′ to the new form represented by the measurement vector z.The transformed vector z and the unknown state vector of the target position x = x z T are related by a nonlinear, but relatively simple formula z = h x + v, which can be deduced from Equation ( 8): The above equation, apart from the geometrical relationships, includes also a vector of measurement errors v, modelling limited accuracy of the measurement vector z.Equation ( 11) can be approximately solved with respect to the unknown position vector x using the OLS algorithm [6][7][8][9] as follows: 3 International Journal of Aerospace Engineering (1) Assume an initial estimated target position vector x = x ẑ T .This can be almost any reasonable guess, but the nearer it is to the true target position, the fewer iterations of OLS are necessary to achieve its required accuracy.The first guess can be based on knowledge of a typical target position, where the first EOTS images and angular data become available during firing tests.The subsequent guesses can be based on the previous target positions, estimated for the previous measurement vector z (2) Calculate the expected measurement vector ẑ = h x (3) Calculate the Jacobian matrix H = ∂h x /∂x x=x of the nonlinear transformation z = h x , at the latest estimated target position x.This yields the following linearized measurement equation: ) Solve Equation ( 12) and calculate a vector of estimation errors Δx.For a full-rank matrix H, Δx is calculated with the use of the Moore-Penrose matrix pseudoinverse: (5) Improve the previous estimate of x by subtracting the calculated correcting term Δx: x ← x − Δx 14 (6) Repeat steps 2 to 5 until Δx falls below a predefined threshold Δx ≤ Δx tr The above OLS algorithm is run for every newly acquired measurement vector z.The Jacobian matrix H in this algorithm is obtained as a matrix of partial derivatives of elements of the measurement vector z with respect to the elements of x, and they can be calculated based on the relationships given in Equation (11) as follows: , , International Journal of Aerospace Engineering Finally, the Jacobian matrix H is as follows: 3.2.WLS Algorithm.The WLS algorithm is another singlepoint solution method [6][7][8][9].It contains the same steps as the OLS, but it additionally weights the measurements depending on their expected accuracy, which improves the accuracy of estimation if the accuracies of measurements are different.If the standard deviations of angle measurement errors are known and given as σ α 1 , σ α 2 , σ λ 1 , σ λ 2 , and there is no cross-correlation of these errors, the following measurement errors covariance matrix R ′ can be written down for the primary measurement vector z ′ : As the measurement vector z is obtained through the transformation z = t z ′ of the primary measurement vector z ′ , the measurement errors of azimuth and elevation angles v ′ also undergo a transformation, which can be obtained through a linearization of the function t * and realized as a multiplication of the v ′ vector through the Jacobian matrix T: Therefore, the covariance matrix R of the artificially constructed measurement vector z is also calculated through a transformation of the covariance matrix of the primary measurement errors R′ as follows: The Jacobian matrix T is obtained as a matrix of partial derivatives of elements of the measurement vector z with respect to the elements of the primary measurement vector z ′ , which are calculated as follows: , 27 , 31 All the other partial derivatives in the Jacobian matrix T are zeroes; thus, finally this matrix is and it should be recalculated for every new set of the measured azimuth and elevation angles.The WLS algorithm can be summarized as follows: (1) Assume an initial estimated target position vector x = x ẑ T in a similar manner as in the OLS algorithm Δẑ in accordance with the assumed accuracy of measurements given by the matrix R: Improve the previous estimate of x by subtracting the calculated correcting term Δx: x ← x − Δx 37 (8) Repeat steps 2 to 7 until Δx falls below a predefined threshold Δx ≤ Δx tr 3.3.EKF Algorithm.The OLS and WLS algorithms are simple, but they do not use information from the previous steps of data processing nor do they comprise any filtering.In contrast, the EKF presented in this section realizes model-based data filtering [1,3].Apart from the measurement model, based on the same geometrical relationships as in OLS or WLS, the EKF requires also modelling of the dynamics of ICP-89 and uses this model for prediction of its positions and other parameters of motion.This way, it preserves past information gained in processing of the older measurement data, and thanks to the use of more information, it is expected to be more accurate than any single-point algorithm.Moreover, in case of outages or incomplete measurement data, the EKF still provides predicted or filtered estimates of reduced accuracy, this way filling the gaps in the missing measurement history, whereas the OLS and WLS do not output any solution in such a case.
3.3.1.Dynamics Model.The use of EKF requires formulation of the dynamics model of the ICP-89 motion in the following form [1]: where f is a nonlinear function, describing changes of the target position and its parameters of flight, x is the state vector, containing variables describing the flight of the ICP-89, G is the matrix of process disturbances, and u is the vector of process disturbances.
The elaboration of the dynamics model for the aerial target imitator ICP-89 requires consideration of its physical properties and its way of operation as well as properties of the atmosphere.The photography of the imitator at the initial phase of flight during firing tests is shown in Figure 3.
During tests, the imitator is launched from the firing stand due to the thrust force generated by its solid propellant.The propellant burns for approximately 0.7 seconds, and for the rest of the flight, the target moves without propulsion or control and follows a free-flight motion pattern.Therefore, the dynamics model of the considered system should be based on equations describing motion of a ballistic object (rocket) [10,11].A full description of such movement consists of six degrees of freedom (6-DOF) equations of motion (EOM), comprising linear and rotational components in the OXYZ frame of reference [11].From the previous considerations, however, we may conclude that the target motion can be restricted to three degrees of freedom (3-DOF) movement in the vertical plane OXZ.ICP-89 in 2D space is illustrated in Figure 4.
The 3-DOF EOM are derived from Newton's second law of dynamics and assume that the motion of the target of a mass m is caused by external forces acting on it, which include gravity (weight) F g = mg and aerodynamical forces: thrust F T , drag F d , and lift F l .The resultant model is nonlinear and can be very sophisticated; therefore, in practice it is often significantly simplified based on engineering experience.A frequently applied simplified set of 3-DOF EOM is as follows [5,10,11]: where x is the horizontal coordinate in OXZ, z is the vertical coordinate in OXZ, v is the target velocity, γ is the pitch angle, α is the angle of attack, F g , F T , F d , F l are external forces (weight, thrust, drag, and lift), C f is the fuel consumption coefficient, and g is the gravitational acceleration.
The drag and lift forces depend on the squared target velocity with respect to the air v, air density ρ, equivalent surface of the target S, and dimensionless drag C d and lift C l coefficients: If we consider that, the first valid measurements of azimuth and elevation angles from the EOTS system become available when the target already moves without propulsion, i.e., F T = 0; the above model given by Equations ( 39)-( 45) can be further simplified.As ICP-89 in its free-flight phase of motion does not consume fuel and does not change its weight, Equation (43) can be eliminated, and the following dynamics model can be formulated: Equations ( 46)-(49) describe an undisturbed target motion in the OXZ plane.In practice, this ideal motion would be disturbed by unknown factors which can be modelled as stochastic processes.The dominant random effects are connected with the moving atmosphere.
If we augment the undisturbed model of the target motion, taking into account the influence of wind, the changes of the target position are given as follows [10]: where w x is the horizontal velocity of the wind and w z is the vertical velocity of the wind.The changeable character of wind influences not only kinematic equations of position (Equations (50)-(51)) but also the equations for velocity and pitch angle as follows [7]: There exist various sophisticated stochastic models of wind, which generally assume existence of strong spatial and temporal correlation of wind velocity [12].However, the usefulness of such models is clearly visible only in estimation of long-term effects of the atmospheric mass movements, e.g., in prediction of power generated by wind farms [13].Their application in the proposed system would significantly complicate the equations, without guarantee of significant improvements of the EKF accuracy.For this reason, authors of some publications neglect the influence of atmospheric mass movement or assume a constant wind model for the whole period of flight [14].
Instead, we propose using a simple stochastic model of wind, which takes into account the random temporal and spatial variability of wind velocity.Our model, however, is only an approximation as it does not consider a correlation which normally exists between the horizontal and vertical velocity of wind fluctuations.The proposed model assumes that the horizontal and vertical components of wind are two independent Wiener processes of unknown initial values, modelled as follows [1]: where u wx and u wz are the independent Gaussian white noises with power spectral densities S wx and S wz .Based on the above assumptions and considerations, we finally obtain the following continuous dynamic model of the target motion, which has the standard form required by the theory and suitable for immediate use in the EKF algorithm:

International Journal of Aerospace Engineering
It is worth noticing that in contrast to the OLS and WLS, the variables in the state vector x include not only target coordinates x, z in the OXZ plane but also other variables, such as velocity v, orientation (pitch angle) γ, and wind velocity components, which influence future states of the object, and therefore, their knowledge is necessary for appropriate prediction of the future target positions and parameters of motion.Thus, the state vector estimated by the EKF contains 6 variables and it is as follows: The above continuous model given by Equation ( 56) can be used directly for propagation ahead of the state vector x.Such an operation requires use of one of the methods of numerical integration [15,16] and is usually realized as a part of data processing in the EKF algorithm [1].The EKF, however, calculates also error covariance matrices, and these operations require linearization of the nonlinear dynamics model.Moreover, the Kalman filter is a discrete-time algorithm and its equations are based on a discrete version of the dynamics model.For the above reasons, it is necessary to perform linearization and discretization (sampling) of the dynamics model.
The linearization is realized through calculation of the Jacobian matrix F = ∂f x /∂x x=x of the nonlinear function f, at the latest estimated state vector x [1].This yields the following, so-called, fundamental matrix of the system: The model sampling consists in calculation of two matrices: the state transition matrix Φ k + 1, k and the covariance matrix of discrete process disturbances Q k .The transition matrix describes the state propagation for one sampling period Δt ahead, from a time step k to the next step of calculation k + 1.There exist several methods of its calculation based on the known fundamental matrix F [1].In the presented algorithm, it was calculated as follows: The covariance matrix of discrete process disturbances Q k represents the uncertainty introduced to the previous knowledge of the state vector at a time step k by cumulative influence of the vector of continuous process disturbances u in the period Δt, between the time steps k and k + 1. Various methods of calculation of this matrix, based on the covariance matrix of continuous process disturbances Q c = E uu T , exist [1].In the presented algorithm, it was calculated as follows: The measurement model used in the EKF is given by Equation ( 11) and is the same as the one used in OLS and WLS.Some EKF equations require additional knowledge of the Jacobian matrix H = ∂h x /∂x x=x of the nonlinear function h * .This matrix can be calculated similarly to the Jacobian matrix used in the OLS, given by Equations ( 15)-( 23).However, it is necessary to take into account that the state vector in the EKF consists of six states, whereas in the OLS of only two states, thus the number of derivatives and the size of the H matrix will be different: Finally, it is clear from Equation ( 11) that all the derivatives with respect to v, γ, w x , w z in Equation ( 63) are zeroes, and on the basis of already calculated derivatives over x and z given by Equations ( 16)-( 22), the H matrix in the EKF is stated as follows: 8 International Journal of Aerospace Engineering The covariance matrix of measurement errors R in the EKF is calculated the same way as in the WLS algorithm, with the use of Equations ( 26) and (34).
3.3.3.Algorithm of Filtration.Estimation of the state vector x, changing with time in accordance with the above formulated dynamics model (Equation ( 56)), based on the measurement vector z, related to the state vector via the measurement model (Equation ( 11)), can be realized by a nonlinear suboptimal filter [3,5].An example of such a filter is the EKF, implemented for that purpose by the authors of this paper.The EKF algorithm can be summarized as follows: (1) Assume an initial estimated state vector x 0 | 0 .
The target position can be initialized in a similar manner as in the OLS and WLS.The rest of the states can be initially assumed to be zeroes (2) Assume an initial value of the covariance matrix of filtration errors P 0 | 0 .Typically, it is initialized as a diagonal matrix with large diagonal values, modelling uncertainty of the initial knowledge of the state vector (3) Calculate the predicted state vector x k + 1 | k at a time step t = k + 1 Δt, through the numerical integration of the deterministic part of the dynamics model x t = f x t , subject to the initial condition (6) Calculate the covariance matrix of discrete process disturbances Q k (Equations ( 60)-( 62)).The Q k matrix must be recalculated at each step due to changing values of the Jacobian matrix F (7) Calculate the covariance matrix of prediction errors: Calculate the Jacobian matrix T k + 1 for the last acquired measurement vector z k + 1 (from Equation (34)).

67
(12) Calculate the corrected state vector x k + 1 | k + 1 at a time step t = k + 1 Δt, by fusing the predicted state vector x k + 1 | k and the current measurement vector z k + 1 : Calculate the covariance matrix of filtration errors: The above set of steps realized during filtration can be interpreted as the initialization (steps 1st and 2nd), the time update or prediction (steps 3rd to 7th), and the measurement update or filtration (steps 8th to 13th).The initialization is executed only once at the beginning of the filter operation, whereas the steps of prediction and filtration are executed recursively for the whole required period of the EKF operation.

Simulation Methodology and Results
After formulating the dynamics and the observation model of the EOTS tracking system, it is possible to run the simulation and analyze the operation of designed estimation algorithms.For the purpose of comparison, two of the presented algorithms, i.e., OLS and EKF, were implemented and tested in the Matlab® scientific environment.The WLS algorithm was not implemented, as the standard deviations of all angle measurements realized in the EOTS are the same; thus, the accuracy offered by WLS should be the same as for OLS.
The appropriate choice of simulation parameters is crucial if we want to obtain valuable and meaningful results.These parameters include the dynamics and the observation model coefficients, as well as initialization parameters for the OLS and EKF.
The dynamics model presented in this paper contains several parameters characteristic for the tracked object and the environment.The geometrical dimensions of the object were chosen to be identical with the real ICP-89 projectile; thus, the calculated surface of the object S = 0 011 m 2 and its mass m = 7 5 kg.The drag coefficient C d was set to 0 12, a value quantitatively consistent with tabular data for conical surfaces.The lift coefficient C l was assumed to be zero, since the surface of the object is axially symmetric with respect to the velocity vector.For simplicity, constant values for the gravitational acceleration g = 9 80665 m/s 2 and air density ρ = 1 2255 kg/m 3 were chosen.The assumed values are representative for typical conditions met during firing tests.
The elements of the covariance matrix of continuous process disturbances Q c , comprising power spectral densities of the white noises in the wind model, were calculated as follows: where σ w is the assumed standard deviation of the Wiener processes of horizontal and vertical wind components at the terminal stage of flight, after a period of Δt.In the simulated scenario σ w = 5m/s and Δt = 24s, which are typical values, representative to the ones met during tests.The primary measurement error covariance matrix R ′ in the EKF was determined based on the exhaustive technical documentation of the EOTS measurement system provided by its developer.The standard deviations of all angle measurements realized in the EOTS are equal σ α 1 = σ α 2 = σ λ 1 = σ λ 2 = 0 68 ′′ in degrees measure (approx.1.05•10 -6 rad).The diagonal elements of R ′ were obtained as squares of the above values, previously converted from arc seconds into radians.
The methodology of simulations included several steps.Firstly, a reference flight trajectory was generated, using a numerical 4th order Runge-Kutta solution to the dynamic equations of motion (Equations ( 50)-( 55)).For such a trajectory, ideal (noiseless) primary measurement vectors z ′ were calculated, with a period of Δt = 33 ms, representing the time separation between two successive measurements in the EOTS system.Next, pseudorandom zero-mean Gaussian sequences of standard deviations 0 68 ′′ were added to the noiseless measurements of observation angles.Such simulated sets of noisy measurements were further processed by the implemented OLS and EKF algorithms.The results of estimation by both algorithms were subsequently analyzed and compared.
The following figures represent a selection of the obtained results.In Figure 5, the reference trajectory, generated in the simulation process, and the ICP-89 flight trajectories estimated by the OLS and EKF are presented.It is noticeable that all the trajectories are close to each other; however, the scale of the figure makes a direct comparison of the estimation accuracy impossible.Such a comparison will be presented further on in this section of the paper.
Figures 6 and 7 contain horizontal and vertical coordinates of the ICP-89 versus time, generated and estimated by the OLS and EKF.
The estimation of velocity, pitch angle, and wind components is not possible in the OLS, but only in the EKF; therefore, the estimation results for those quantities, shown in Figures 8-11, comprises only EKF estimates, together with their reference values, generated during simulations.
To facilitate a comparison of accuracy of the OLS and EKF algorithms, errors of horizontal and vertical coordinates of the ICP-89 target imitator are presented in Figures 12 and  13.These figures give a better insight into the differences between the results from EKF and OLS algorithms than those previously presented in Figures 5-7.As can be seen, the EKF significantly outperforms the OLS solution.One can also observe that the estimation accuracy is not constant for the whole period of target tracking.
This effect is due to changeable conditions of observation, i.e., changing geometrical dilution of the precision factor, which depends on the angle of intersection of the lines between the EOTS stations and the target.For the assumed simulation scenario, this angle decreases in the terminal phase of flight, which increases the geometrical dilution of precision and the target localization errors.
It is not possible to directly compare the two implemented algorithms with respect to the accuracy of estimation of the other state variables (velocity, pitch angle, and wind velocity components), because the OLS does not provide their estimates.Therefore, in the following Figures 14-17 only EKF errors are presented.They are calculated as differences between respective estimates and the actual state values generated in the process of simulation.
The quantitative comparison between the OLS and EKF can be based on the RMS estimation errors of both algorithms.The results of such a comparison for the estimates

Conclusions
The paper presented a problem of estimation of a ballistic object trajectory and its solutions based on the OLS and EKF algorithms.Both algorithms designed, implemented, and tested by the authors of the paper proved to be useful tools for estimation of an aerial target position and its parameters of motion.The usefulness of the presented approach was demonstrated for a practical problem of estimation of a trajectory and parameters of motion of the aerial target imitator ICP-89, based on the measurements from the electrooptical tracking system with two cinetheodolite stations.The simulations of such a system have shown that the target trajectory can be estimated with high accuracy.As expected, tracking errors depend on the geometry of the EOTS system and increase when the ICP-89 is far from the cinetheodolite stations.This conclusion should be considered            when the shooting tests are planned and used for optimal location of the observation stations in the firing range.The presented comparisons of the OLS and EKF accuracy show that the EKF can estimate the target position with 2-3 times better accuracy than the OLS.Moreover, it can directly estimate the target velocity, pitch angle, and components of wind, which is not possible in the OLS.Therefore, the EKF algorithm should be a preferred solution in realization of a practical tracking system.
The authors are currently continuing their works on the presented system aimed at implementation of the presented algorithms, especially the EKF, for processing real measurements from the EOTS system, recorded during live firing tests.The works are aimed also at improving the assumed ballistic model of ICP-89 as well as the wind model.It should be emphasized that the results obtained so far are very promising and it is expected that the EKF algorithm, elaborated by the authors, can be implemented in the software used by MIAT and replace algorithm currently used for fire tests.

Figure 4 :
Figure 4: ICP-89 imitator outline in 2D space with essential variables describing its movement.

4 )( 5 )
Calculate the Jacobian matrix F (Equation (58)), for x = x k | k Calculate the state transition matrix Φ k + 1, k of the discrete version of the dynamics model (Equation (59)).

Figure 10 :
Figure 10: Simulated and estimated horizontal component of wind.

Figure 11 :
Figure 11: Simulated and estimated vertical component of wind.

Figure 12 :
Figure 12: Error of estimation of horizontal position of ICP-89.

Figure 13 :
Figure 13: Error of estimation of vertical position of ICP-89.

Figure 16 :
Figure 16: Error of estimation of horizontal component of wind.

Figure 17 :
Figure 17: Error of estimation of vertical component of wind.

Table 1
. It should be noticed, however, that the velocity and pitch angle are not directly available from OLS and they have been calculated based on three successive position estimates.Thus, the figures shown in Table1, representing RMS velocity and pitch angle errors of OLS, are approximate and may only serve for a crude assessment of the OLS accuracy.

Table 1 :
Comparison of OLS and EKF target localization errors.