Tracking of Flying Objects on the Basis of Multiple Information Sources

1. Introduction

In the paper, a new approach to the estimation of the localization of flying objects on the basis of multiple radiolocation sources has been proposed. The basic purpose of the radiolocation is the detection and tracking of objects. The main task of a tracking system is the estimation and prediction of localization and motion of tracked objects. In the classical tracking systems, the estimation of localization and motion of tracked objects is based on the series of measurements performed by a single radar. Measurements obtained from other sources are not taken directly into account. It is possible to assess the estimate of a tracking object on the basis of the different sources of information, applying, for example, least squares method (LSM). However, such solution is seldom applied in practice, because necessary formulas are rather complicated. In this paper, a new approach is proposed. The key idea of the proposed approach is so called matrix of precision. This approach makes possible tracking not only on the basis of radar signals, but also on the basis of bearings. It makes also possible the tracking of objects on the basis of multiple information sources. Simplicity is the main attribute of proposed estimators. Their iterative form corresponds well with the prediction-correction tracking model, commonly applied in radiolocation. In the paper numerical examples presenting advantages of the proposed approach are shown.

The paper consists of seven parts. Introduction is the first one.

In the second part the idea of the matrix of precision is presented and it is demonstrated how it can be used to uniformly describe the dispersion of measurement. Traditionally, the measurement dispersion is described by a matrix of covariance. Formally, the matrix of precision is an inverse of the matrix of covariance. However, these two ways of description are not interchangeable. If an error distribution is described by the singular matrix of precision then the corresponding matrix of covariance does not exist, more precisely, it contains infinite values. It means in practice that some components of a measured vector are not measurable, in other words, an error of measurement can be arbitrarily large. The application of the matrix of precision makes possible an uniform description of measurements taken from various sources of information, even if measurements come from different devices and measure different components. Zero precision corresponds to components which are not measured.

In the third part, the problem of estimation of stationary parameters is formulated. Using the matrix of precision, the simple solution of the problem is presented. It appears that the best estimate is a weighted mean of measurements, where the weights are the matrices of the precision of measurements. It has been proved that the proposed solution is equivalent to the least squares method (LSM). Additionally it is simple and scalable.

In the fourth part of this paper the problem of the estimation of states of dynamic systems, such as a flying aircraft, is formulated. Traditionally, for such an estimation the Kalman filter is applied. In this case the uncertainty of measurement is described by the error matrix of covariance. If the matrix of precision is singular, it is impossible to determine the corresponding matrix of covariance and utilize the classical equation of Kalman filter. In the presented approach this situation is typical. It appears that there is such a transformation of Kalman filter equations, that the estimation based on the measurements with error described by the singular matrix of precision, can be performed. Such a transformation is presented and its correctness is proved.

In the fifth part, numerical examples are presented. They show the usability of the concept of matrix of precision.

In the sixth part, the summary and conclusions are shown, as well as the practical application of presented idea is discussed. The practical problems which are not considered in the paper are also pointed out.

2. The concept of precision matrix

Traditionally in order to describe the degree of the dispersion distribution the covariant matrix is used. There exist another statistics which can be used to characterize dispersion of the distribution. However, the covariance matrix is the most popular one and in principle the precision matrix is not used. Formally the precision matrix is the inverse of covariance matrix.

$$W_x=C_x{}^{-1}$$E1

The consideration may take much simpler form if the precision matrix is used. Note that precision matrix is frequently used indirectly, as inverse of covariance matrix. For example, in the well-known equation for the density of multi-dimensional Gauss distribution:

$$\sqrt{\frac{\left|C^{-1}\right|}{{\left(2\pi \right)}^k}}{e}^{\left(-\frac{{(x-m_x)}^T{C}^{-1}(x-m_x)}{2}\right)}$$E2

The equation (1) can not be used if the covariance matrix is singular. As long as the covariance (precision) matrix is not singular, the discussion which statistics is better to describe degree of the dispersion distribution is as meagre as discussion about superiority of Christmas above Easter. Our interest is in analysis of extreme cases (e.g. singular matrix).

At first, the singular covariance matrix will be considered. Its singularity means that some of the member variables (or their linear combinations) are measured with error equal to zero. Further measurement of this member variable makes no sense. The result can be either the same, or we shall obtain contradiction if the result would be different. The proper action in such a case is modifying the problem in such a manner that there is less degrees of freedom and the components of the vector are linearly independent. The missing components computation is based on the linear dependence of the variables. In practice, the presence of the singular covariance matrix means that the linear constraints are imposed on the components of the measured vectors and that the problem should be modified. The second possibility (infinitely accurate measurements) is impossible in the real world.

In distinction from the case of the singular covariance matrix, when the measurements are described by the singular precision matrix, measurements can be continued. The singularity of the precision matrix means that either the measurement does not provide any information of one of components (infinitive variation) or there exist linear combinations of member variables which are not measured. This can results from the wrong formulation of the problem. In such a case, the system of coordinates should be changed in such a manner that the variables which are not measured should be separated from these which are measured. The separation can be obtained by choosing the coordination system based on the eigenvectors of precision matrix. The variables corresponding to non-zero eigenvalues will then be observable.

There is also the other option. The singular precision matrix can be used to describe the precision of measurement in the system, in which the number of freedom degrees is bigger then the real number of components, which are measured. Then all measurements may be treated in a coherent way and the measurements from various devices may be taken into account. Thus all measurements from the devices which do not measure all state parameters can be included (i.e. direction finder). In this paper we are focusing on second option.

Each measuring device uses its own dedicated coordinate system. To be able to use the results of measurements performed by different devices, it is necessary to present all measured results in the common coordinate system. By changing coordinate system the measured values are changed as well. The change involves not only the digital values of measured results, but also the precision matrix describing the accuracy of particular measurement. We consider the simplest case namely the linear transformation of variables. Let $$y$$ denote the vector of measurements in common coordinate system and $$x$$ denote the vector of measurements in the measuring device dedicated coordinate system. Let $$x$$ and $$y$$ be related as follows:

$$y=Ax+b$$E3

In this case the covariance of measurement changes according to formula:

$$C_y=A{C}_x{A}^T$$E4

The precision matrix as inverse of covariance matrix changes according to formula:

$$W_y=C_y^{-1}=A^T{}^{-1}{W}_x{A}^{-1}$$E5

Consider the case when the transformation $$A$$ is singular. In this case, we can act in two ways. If the precision matrix $$W_x$$ is not singular it can be inverted and the covariance matrix $$C_x$$ can be find. Using the formula (4) and the obtained result $$C_y$$ it is possible to invert back, to obtain the precision matrix$$W_y$$. The second method consists in increasing transformation $$A$$ in such way to be invertible. It is done by writing additional virtual lines to the transformation matrix$$A$$. The new transformation matrix is now as follows:

$$A\text{'}=\left[\begin{array}{c}{A}_{}\\ A_v\end{array}\right]$$E6

Additional rows of this matrix should be chosen in such way that the obtained precision matrix has block-diagonal form:

$$W_y\text{'}=\left[\begin{array}{cc}{W}_A& 0\\ 0& W_v\end{array}\right]$$E7

it means that additional virtual lines of transformation $$A\text{'}$$ should be chosen to fulfill the following relationship:

$$A_v^T{}^{-1}{W}_x{A}_v{}^{-1}=0$$E8

In such a way, the additional virtual lines do not have any influence on the computation result of precision matrix in the new coordinate system.

Consider the specific case. Let the measuring device measures target OX and OY positions independently, but with different precision. Measurement in the direction of the axis OX is done with accuracy $$w_{OX}={0.25}^2$$ and in the direction of axis OY measured precision is$$w_{OY}={0.2}^2$$. Then the precision matrix has a form:

$$W_x=\left[\begin{array}{cc}{0.25}^2& 0\\ 0& {0.2}^2\end{array}\right]$$E9

We are interested in accuracy of the new variable which is a linear combination of these variables. Let $$y=Ax$$ where:

$$A=\left[\begin{array}{cc}-1& 2\end{array}\right]$$

Matrix $$A$$ is not invertible. We add an additional row so that the matrix $$A\text{'}$$ is invertible and formula (8) is fulfilled. General form of the matrix fulfilling this condition is of the form:

$$A\text{'}=\left[\begin{array}{cc}-1& 2\\ 50n& 16n\end{array}\right]$$E10

We chose$$n=1$$. In that case:

$$A{\text{'}}^{-1}=\frac{1}{116}\left[\begin{array}{cc}-16& 2\\ 50& 1\end{array}\right]$$E11

The precision matrix in the new coordinate system has the following form:

$$W_y=A{\text{'}}^{-1}{}^T{W}_{x}A{\text{'}}^{-1}=\frac{1}{116}\left[\begin{array}{cc}1& 0\\ 0& 0.0025\end{array}\right]$$E12

We ignore additional variables which have been added to make $$A\text{'}$$ invertible. Finally the accuracy of the measurement of variable $$y$$ is 1/116.

The same result can be obtained using the first method. Then the covariance matrix is of the form:

$$C_x=\left[\begin{array}{cc}{4}^2& 0\\ 0& 5^2\end{array}\right]$$E13

The value of the covariance of variable $$y$$ according to formula (4) is:

$$C_y=\left[\begin{array}{cc}-1& 2\end{array}\right]\left[\begin{array}{cc}16& 0\\ 0& 25\end{array}\right]\left[\begin{array}{c}-1\\ 2\end{array}\right]=116$$E14

While the second method is more complicated it is nevertheless more general, and permits to find the precision matrix in the new coordinate system even if the precision matrix is singular.

In a case when it is necessary to make nonlinear transformation of variables$$y=f(x)$$, the linearization of transformation is proposed. Instead of matrix $$A$$ the Jacobian matrix of transformation $$f$$ is used:

$$J=\left[\frac{\partial f}{\partial x}\right]$$E15

In order to find the value of partial derivative, it is necessary to know point $$x_0$$ around which linearization is done.

3. Estimation of stationary parameters

In engineering it is often necessary to estimate certain unknown value basing on several measurements. The simplest case is, when all measurements are independent and have the same accuracy. Then, the intuitive approach i.e. using the representation of results as arithmetic average of measurements is correct (it leads to minimal variation unbiased estimator). The problem is more complicated if the measured value is a vector and individual measurements have been performed with different accuracy. Usually, to characterize measurements accuracy the covariance matrix of measurement error distribution is used. Application of precision matrix leads to the formulas which are more general and simple. Further the precision matrix describes the precision of measurements when the number of degrees of freedom is smaller than the dimension of the space in which the measurement has been done. For example, direction finder measures the direction in the three dimensional space. The direction finder can be used as a radar with an unlimited error of measured distance. In this case, the covariance matrix of the measurement does not exist (includes infinite values). But there exists the precision matrix (singular). Information included in the precision matrix permits to build up optimal (minimal variation unbiased) linear estimator.

3.1. Formulation of the problem

Consider the following situation. The series of n measurements $$x_i$$ of a vector $$x$$ value is done under the following assumptions:

1. There is no preliminary information concerning the value of the measured quantity$$x$$.

2. The measured quantity is constant.

3. All measurements are independent.

4. All measurements are unbiased (expected value of the measured error is zero).

5. All measurements can be presented in common coordinate system and the precision matrix P_i which characterizes the measurement error distribution in these coordinates, can be computed.

The problem is to find the optimal linear estimator of unknown quantity $$x$$ and to find its accuracy. By optimal estimator, we understand the unbiased estimator having minimal covariance. By estimation of precision we understand presenting its covariance (or precision) matrix.

3.2. The solution

The first step consists in presenting all measurement results in the common coordinate system. The precision matrices must be recalculated to characterize the measurement error distribution in the common coordinate system. The method of the recalculation the precision (covariance) matrixes in the new coordinated system have been presented before. If all measurements have been done in the same coordinate system, the first step is omitted.

The unbiased minimal variation estimator of unknown quantity $$x$$ is a weighed average of all measurements$$x_i$$, where the weights are the precision matrices $$W_i$$ of the individual measurements. The precision matrix $$W_{\stackrel{⌢}{x}}$$ of such estimator is equal to the sum of the precision matrixes $$W_i$$ of all measurements.

$$\widehat{x}=\frac{{\displaystyle \sum W_i{x}_i}}{{\displaystyle \sum W_i}}$$E16

$$W_{\widehat{x}}={\displaystyle \sum W_i}$$E17

In formula (16) the results of the individual measurements are represented as column vectors. Formulas (16-17) also can be presented in the iterative form (18-21) :

$${\stackrel{⌢}{x}}_1=x_1$$E18

$$W_{{\widehat{x}}_1}=W_1$$E19

$${\widehat{x}}_{i+1}=\frac{W_{{\widehat{x}}_i}{\widehat{x}}_i+W_{i+1}{x}_{i+1}}{W_{{\widehat{x}}_i}+W_{i+1}}$$E20

$$W_{{\widehat{x}}_{i+1}}=W_{{\widehat{x}}_i}+W_{i+1}$$E21

Proof:

The presented formula represents the particular case of the general estimator of state $$x$$ of the linear system, in which information about $$x$$ is accessible via$$y=Ax+e$$, where $$y$$ is the vector of the measured results, $$A$$is observation which is the projection of the system state on the measured vector, $$e$$is random noise having average value and covariance$$Ee{e}^T=C$$, which represents the measurement errors. In such case, the best linear unbiased estimator has the following form:

$$\widehat{x}=C^{-1}A{\left(A^T{C}^{-1}A\right)}^{-1}y$$E22

In our case, the formula is radically simplified. All measurement results are presented in the same coordination system, therefore the observation Matrix has the following form:

$$A=\left[\begin{array}{c}\begin{array}{cccc}1& 0& \mathrm{...}& 0\\ 0& 1& \mathrm{...}& 0\\ \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}\\ 0& 0& \mathrm{...}& 1\end{array}\\ \mathrm{.......}\\ \begin{array}{cccc}1& 0& \mathrm{...}& 0\\ 0& 1& \mathrm{...}& 0\\ \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}\\ 0& 0& \mathrm{...}& 1\end{array}\end{array}\right]$$E23

Since all measurements are independent, the covariance matrix has the block-diagonal form:

$$C=\left[\begin{array}{cccc}{C}_1& 0& \mathrm{...}& 0\\ 0& C_2& \mathrm{...}& 0\\ \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}\\ 0& 0& \mathrm{...}& C_n\end{array}\right]$$E24

$$C^{-1}=\left[\begin{array}{cccc}{C}_1^{-1}& 0& \mathrm{...}& 0\\ 0& C_2^{-1}& \mathrm{...}& 0\\ \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}\\ 0& 0& \mathrm{...}& C_n^{-1}\end{array}\right]=\left[\begin{array}{cccc}{W}_1& 0& \mathrm{...}& 0\\ 0& W_2& \mathrm{...}& 0\\ \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}\\ 0& 0& \mathrm{...}& W_n\end{array}\right]$$E25

After simplification and changing of the multiplication order (symmetrical matrix) we obtain:

$$\widehat{x}=C^{-1}A{\left(A^T{C}^{-1}A\right)}^{-1}y={\displaystyle \sum W_i{\left({\displaystyle \sum W_j}\right)}^{-1}{y}_i={\displaystyle \sum W_i{y}_i{\left({\displaystyle \sum W_j}\right)}^{-1}}}$$E26

As we can see, it is another form of formula (16). Its interpretation is easy to read and easy to remember. The optimal linear estimator represents the weighed average of all measurements, when weighs are the precision matrices of individual measurements.

4. Modified Kalman filter

In many cases it is necessary to estimate the changing state vector of the system. The value of the vector changes between successive measurements. Such situation exists in radiolocation when the position of moving target is tracked. Restricting ourselves to linear systems, the solution of this problem is known as the Kalman filter. In the Kalman filter theory, the model of system is presented as the pair of equations:

$$x_{i+1}=F{x}_i+q_{i+1}$$E27

$$z_{i+1}=H{x}_{i+1}+r_{i+1}$$E28

The first equation describes dynamics of system in the discrete time domain. Vector $$q$$ models the influence of internal noise. It is assumed that expected value $$q$$ is zero and it is independent from the vector of state value$$x$$, as well as of the values of vector $$q$$ at the previous iterations:

$${\rm E}{q}_i=0$$E29

$$\underset{\text{i}\ne \text{j}}{{\displaystyle \text{Ε}}}{q}_i{q}_j=0$$E30

$${\rm E}{q}_i{q}_i^T=Q$$E31

The covariance of internal noise $$q$$ is described by matrix$$Q$$. The second equation describes output $$z$$ of the system which depends on actual state of the vector$$x$$. The vector $$r$$ models the influence of noise on the measuring process. Further it is assumed that the expected value $$r$$ is equal to zero and it is independent from the value of the state vector$$x$$, internal noise vector$$q$$, as well as of the vector $$r$$ from the previous iterations:

$${\rm E}{r}_i=0$$E32

$$\underset{\text{i}\ne \text{j}}{{\displaystyle \text{Ε}}}{r}_i{r}_j=0$$E33

$${\rm E}{r}_i{r}_i^T=R$$E34

It is assumed that we have estimate of the value$${\widehat{x}}_i$$, accuracy of which is described by covariance matrix $$P_i$$ and as well as the new measurement $$z_{i+1}$$ which is accurately described by the covariance matrix$$R_{i+1}$$. How taking into consideration the last measurement can we obtain the estimated value of the state vector at the time$$i+1$$? The solution of this problem leads to the classic Kalman filter equation. Kalman filter can be presented in the iteration form:

new estimate =state prediction + correctioncorrection = gain*(new measurement – measurement prediction)

The equation of the state prediction has the following form:

$${\widehat{x}}_{i+1|i}=F{\widehat{x}}_i$$E35

The covariance matrix of the state prediction has the following form:

$$P_{i+1|i}=F{P}_i{F}^T+Q_{i+1}$$E36

The equation of the measurement prediction has the form:

$${\widehat{z}}_{i+1|i}=H{\widehat{x}}_{i+1|i}$$E37

The Kalman gain is defined as:

$$K_{i+1}=P_{i+1|i}{H}^T{\left[H{P}_{i+1|i}{H}^T+R_{i+1}\right]}^{-1}$$E38

The new estimate of the state has the following form:

$${\widehat{x}}_{i+1}={\widehat{x}}_{i+1|i}+K_{i+1}\left(z_{i+1}-{\widehat{z}}_{i+1|i}\right)$$E39

The covariance matrix of new estimate of the state has the following form:

$$P_{i+1}=\left[1-K_{i+1}H\right]P_{i+1|i}$$E40

Equations (27-40) constitutes the description of the classic Kalman filter. The Kalman filter can be considered as the system in which the input consists of string of measurements $$z_i$$ together with the covariance matrix $$R_i$$ which characterizes the accuracy of the i-th measurement and with the output being new estimate of the state vector $${\widehat{x}}_i$$ together with the covariance matrix $$P_i$$ which describes accuracy of this estimate. In the case when the precision of measurements is described by the singular precision matrix$$W_i$$, it is impossible to obtain corresponding covariance matrix $$R_i$$ and the classic Kalman filter equation can not be used (covariance matrix is required). Such situation happens for instance when the movement of the object based on the bearings is tracked. The bearing can be treated as a degenerated plot in which precision of distance measurement is equal zero. Thus bearing is the plot described by the singular precision matrix. Determination of the covariance matrix is impossible. However, it is not necessary. The only equation in which the covariance matrix $$R$$ is used is equation (38). The covariance matrix in this equation is inverted by indirect way. A part of the equation (38) in which matrix $$R$$ exists has form:

$${\left[HP{H}^T+R\right]}^{-1}$$E41

It is easy to prove that:

$$A\left[R^{-1}+A^{-1}\right]R=\left[A+R\right]$$E42

We invert both sides:

$$R^{-1}{\left[R^{-1}+A^{-1}\right]}^{-1}{A}^{-1}={\left[A+R\right]}^{-1}$$E43

Hence:

$${\left[HP{H}^T+R\right]}^{-1}=R^{-1}{\left[R^{-1}+{\left[HP{H}^T\right]}^{-1}\right]}^{-1}{\left[HP{H}^T\right]}^{-1}$$E44

It can be written as:

$${\left[HP{H}^T+R\right]}^{-1}=W{\left[W+{\left[HP{H}^T\right]}^{-1}\right]}^{-1}{\left[HP{H}^T\right]}^{-1}$$E45

Finally, the equation (38) has form:

$$K_{i+1}=P_{i+1|i}{H}^T{W}_{i+1}{\left[W_{i+1}+{\left[H{P}_{i+1|i}{H}^T\right]}^{-1}\right]}^{-1}{\left[H{P}_{i+1|i}{H}^T\right]}^{-1}$$E46

The only problem which can appear here is the possible singularity of the matrix$$H{P}_{i+1|i}{H}^T$$. This is the covariance matrix of the measurement prediction. If this matrix is singular, it means that we have certainty about measurement result of some parts and act of filtration have no sense as it has been discussed before. During normal filtration process such situation never happens. In such transformed form of Kalman filter the measurements which accuracy is described by the precision matrix can be used as the input, even if this matrix is singular. The price which is to pay is the necessity to twice invert the matrix which dimension corresponds to the dimension of the observation vector.

5. Examples

5.1. Estimation of the static system.

In order to demonstrate the proposed technique, the following example is used. The position of tracking target is determined using three measurements. The first measurement comes from the radar. The second one comes from the onboard GPS device. The third one comes from the direction finder. The radar is a device which measures three coordinates: slant range ($$R$$[m]), azimuth ($$\beta$$[rad]), and elevation angle ($$\varphi$$[rad]). For our simulation we assume that standard deviation of particular errors are: $$\left[\begin{array}{ccc}1500& 0.001& 0.001\end{array}\right]$$([m rad rad]) i.e. the radar works as an altimeter precisely measuring elevation and azimuth while the slant range is not precisely measured. As far as GPS device is concerned we know that this device measures position with accuracy of up to 100 m and this measurement can be presented in any Cartesian coordinate system. The altitude is not measured by GPS. The standard deviation of bearing error for direction finder is [0.001] ([rad]). We know also very precisely the location of the direction finder. For this example the simulation has been done. We have chosen Cartesian coordinate system located at the radar as the global coordinate system. Furthermore, we have assumed that the tracking target is at the point$$x=\left[\begin{array}{ccc}30000& 40000& 12000\end{array}\right]$$. The radar position is then at$$R=\left[\begin{array}{ccc}0& 0& 0\end{array}\right]$$. Bearing finder is in the point$$N=\left[\begin{array}{ccc}70000& 10000& 0\end{array}\right]$$. We obtain following measurements results:

$$z_{Rad}=\left[\begin{array}{ccc}\text{50771}& \text{0,64363}& \text{0,23388}\end{array}\right]$$E47

$$z_{GPS}=\left[\begin{array}{cc}\text{30029}& \text{39885}\end{array}\right]$$E48

$$z_{Nam}=\left[\text{-0,92733}\right]$$E49

Basing on these results we want to obtain the estimator of$$x$$. The first step is to present the results of measurements in the common coordinate system. In our case it is $$XYZ$$ coordinate system connected to the radar. Now, some values should be attached to missing variables. If the transformation is linear, we can assign them any values. If the transformation of variables is nonlinear, we linearize around an unknown point. The values of the missing variables are estimated based on results of other measurements. After assigning the missing variables we obtain:

$$z_{Rad}=\left[\begin{array}{ccc}\text{50771}& \text{0}\text{.64336}& \text{0}\text{.23388}\end{array}\right]$$E50

$$z_{GPS}=\left[\begin{array}{ccc}30029& 39885& 11766\end{array}\right]$$E51

$$z_{Nam}=\left[\begin{array}{ccc}51363& \text{-0}\text{.92733}& 0.23113\end{array}\right]$$E52

Transforming to new coordinate system, we obtain following results:

$$z_{Rad}^{XYZ}=\left[\begin{array}{ccc}\text{29638}& \text{39507}& \text{11766}\end{array}\right]$$E53

$$z_{GPS}^{XYZ}=\left[\begin{array}{ccc}\text{30029}& \text{39885}& 11766\end{array}\right]$$E54

$$z_{Nam}^{XYZ}=\left[\begin{array}{ccc}\text{30001}& \text{39997}& \text{11766}\end{array}\right]$$E55

In order to determine the precision matrix we invert the corresponding covariance matrices:

$$W_x=C_x{}^{-1}$$E56

The accuracy of the GPS device is 100 m. Therefore we consider 100 as the standard deviation of error distribution.

$$W_{GPS}=C_{GPS}^{-1}={\left[\begin{array}{cc}{100}^2& 0\\ 0& {100}^2\end{array}\right]}^{-1}=\left[\begin{array}{cc}0.0001& 0\\ 0& 0.0001\end{array}\right]$$E57

The standard deviation of the direction finder is 0.001 rad. The variance and precision are accordingly:

$$W_{Nam}=C_{Nam}^{-1}={\left[{0.001}^2\right]}^{-1}={10}^6$$E58

After inserting zeros to the rows and columns corresponding to missing variables we obtain the following precision matrices:

$$W_{Rad}=\left[\begin{array}{ccc}0.(4)e-6& 0& 0\\ 0& {10}^6& 0\\ 0& 0& {10}^6\end{array}\right]$$E59

$$W_{GPS}=\left[\begin{array}{ccc}0.0001& 0& 0\\ 0& 0.0001& 0\\ 0& 0& 0\end{array}\right]$$E60

$$W_{Nam}=\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& {10}^6\end{array}\right]$$E61

When applying the linear transformation of variables: $$y=Tx+a$$the covariance matrix changes according to the formula:$$C_y=T{C}_x{T}^T$$. In the case of the nonlinear transformation:$$y=Y(x)$$, we replace the matrix $$T$$ by the Jacobian matrix$$J$$, which is the linearization of this transformation.

$$J={\left[\begin{array}{ccc}\frac{\partial y_1}{\partial x_1}& \mathrm{...}& \frac{\partial y_1}{\partial x_n}\\ \mathrm{...}& \mathrm{...}& \mathrm{...}\\ \frac{\partial y_n}{\partial x_1}& \mathrm{...}& \frac{\partial y_n}{\partial x_n}\end{array}\right]}_{x=\stackrel{⌢}{x}}$$E62

$$C_y=J{C}_x{J}^T$$. Similarly we obtain:$$W_y=J^{-1^T}{W}_x{J}^{-1}$$. In our case, the radar and the direction finder (after inserting the missing variables) operate in the spherical coordinates which are related to global system by the set of formulas:

$$X=X_0+R\cos (\phi )\sin (\beta )$$E63

$$Y=Y_0+R\cos (\phi )\cos (\beta )$$E64

$$Z=Z_0+R\sin (\phi )$$E65

The Jacobian matrix of such transformation has the following form:

$$J=\left[\begin{array}{ccc}\cos (\phi )\sin (\beta )& R\cos (\varphi )\cos (\beta )& -R\sin (\phi )\sin (\beta )\\ \cos (\phi )\cos (\beta )& -R\cos (\phi )\sin (\beta )& -R\sin (\phi )\cos (\beta )\\ \sin (\phi )& 0& R\cos (\phi )\end{array}\right]$$E66

Transforming variables from the local coordinate system of the measuring device to global system, we obtain as follows:

For the radar:

$$W_{Rad}^{XYZ}=J_{z_{Rad}}^{-1}{}^T{W}_{Rad}{J}_{z_{Rad}}^{-1}=\left[\begin{array}{ccc}\text{26}\text{.998}& \text{-18}\text{.659}& \text{-5,2424}\\ \text{-18}\text{.659}& \text{16}\text{.124}& \text{-6,9881}\\ \text{-5,2424}& \text{-6}\text{.9881}& \text{36}\text{.713}\end{array}\right]e-5$$E67

In case of the GPS device we do not have any change of variables. We are only adding the zero values corresponding to additional variable to the precision matrix:

$$W_{GPS}^{XYZ}=\left[\begin{array}{ccc}0.0001& 0& 0\\ 0& 0.0001& 0\\ 0& 0& 0\end{array}\right]$$E68

For the direction finder:

$$W_{Nam}^{XYZ}=J_{z_{Nam}}^{-1}{}^T{W}_{Nam}{J}_{z_{Nam}}^{-1}=\left[\begin{array}{ccc}\text{14}\text{.4}& \text{19}\text{.201}& <1\text{e-15}\\ \text{19}\text{.201}& \text{25}\text{.604}& <1\text{e-16}\\ <\text{1e-15}& <\text{e-16}& <1\text{e-32}\end{array}\right]e-5$$E69

Finally, we obtain the position estimator:

$$x^{*}=\frac{z_{Rad}^{XYZ}{W}_{Rad}^{XYZ}+z_{GPS}^{XYZ}{W}_{GPS}^{XYZ}+z_{Nam}^{XYZ}{W}_{Nam}^{XYZ}}{W_{Rad}^{XYZ}+W_{GPS}^{XYZ}+W_{Nam}^{XYZ}}$$E70

resulting in:

$$x^{*}=\left[\begin{array}{ccc}\text{30008}& \text{39973}& \text{11908}\end{array}\right]$$E71

Note that utilization of the measurements from GPS device and direction finder, which do not measure the height, actually results in some improvement of altitude estimation. This is due to the fact that the variables which provide the radar measurements in XYZ coordinates are not independent, but they are strongly correlated. The linear combination of them is measured with high accuracy. Thus, by improving the estimation of XY variables we also improve the estimation of the variable in Z direction.

5.2. Estimation in the dynamic system

The next example presents the application of the proposed technique to the estimation of the parameters of a dynamic system. It is the most typical situation when the racking target is moving. Our purpose is the actualization of motion parameters of tracking target based on three measurements. We assume that the tracking is two dimensional. The following parameters are tracked: $$x$$- the position on axis OX; $$vx$$- the motion speed on axis OX; $$y$$- the position on axis OY; $$vy$$- the motion speed on axis OY. We consider the flight altitude $$z$$ as fixed and equal 2000. Finally, we have the following data:

At time:

$$T0=100$$E72

Based on previous measurements, we obtain estimation of the initial parameters of tracked target motion:

$$x_0=\left[\begin{array}{c}x\\ vx\\ y\\ vy\end{array}\right]=\left[\begin{array}{c}-259.5\\ 126.7\\ 60075\\ -114.3\end{array}\right]$$E73

The covariance matrix of this estimation has the form:

$$C_{x_0}=\left[\begin{array}{cccc}{600}^2& 0& 0& 0\\ 0& {20}^2& 0& 0\\ 0& 0& {600}^2& 0\\ 0& 0& 0& {20}^2\end{array}\right]$$E74

At time:

$$T1=108$$E75

The tracked target has been detected by the radar R1 which is located at the beginning of our coordinate system

$$R1=\left[\begin{array}{c}0\\ 0\\ 0\end{array}\right]$$E76

The parameters of the detection are:

$$z_{R1}=\left[\begin{array}{c}59025\\ 0.1169\\ 0.2122\end{array}\right]$$E77

We assume that the radar measurement is unbiased. The covariance matrix of the measurement error has the following form:

$$C_{R1}=\left[\begin{array}{ccc}{55}^2& 0& 0\\ 0& {0.08}^2& 0\\ 0& 0& {0.15}^2\end{array}\right]$$E78

The next position measurements of the tracking target took place at the time:

$$T2=113.8$$E79

We obtain the bearing of tracked target. The direction finder was located at the point:

$$N2=\left[\begin{array}{c}10000\\ 40000\\ 100\end{array}\right]$$E80

The bearings value is:

$$z_{N2}=\left[-0.3853\right]$$E81

The covariance measurement error is:

$$C_{N2}=\left[{0.05}^2\right]$$E82

At the near time:

$$T3=114$$E83

We obtain bearing from another direction finder which was located an point:

$$N3=\left[\begin{array}{c}30000\\ 50000\\ 120\end{array}\right]$$E84

The bearing value is:

$$z_{N3}=\left[-1.2433\right]$$E85

The covariance measurement error is:

$$C_{N3}=\left[{0.05}^2\right]$$E86

Traditionally, to solve such a problem, the theory of Kalman filter is used. In this theory two assumptions are usually accepted. First that all measurements are carried out by the same device. Second, that the measurements are done in regular time intervals. None of the above assumptions is satisfied in our example. To address the issue of different measuring device, we assume that the measurements are always presented in the common global coordinate system, and their dispersion is described by the precision matrix which can be singular. Using the modified formulas of Kalman filters (46), the estimation can be done, even if the precision matrix is singular. In order to address the issue of the irregular time intervals of the measurements, the Kalman filter has been modified to consider the changes of time interval between individual measurements. It means that the matrix $$F$$ describing dynamics of the tracked process and matrix $$Q$$ described covariance of the internal noise of tracked process depend on time $$\Delta t$$ which elapses from the last measurement. We assume the linear model of the tracked target movement in which velocity does not have systematic changes. Therefore the equation of motion has the following form:

$$\{\begin{array}{l}x(t)=x_0+v{x}_0*t\hfill \\ vx(t)=v{x}_0\hfill \\ y(t)=y_0+v{y}_0*t\hfill \\ vy(t)=v{y}_0\hfill \end{array}$$E87

We introduce the following state vector of Kalman filter:

$$x=\left[\begin{array}{c}x\\ vx\\ y\\ vy\end{array}\right]$$E88

Therefore, the dynamic matrix $$F(\Delta t)$$ has the following form:

$$F\left(\Delta t\right)=\left[\begin{array}{cccc}1& \Delta t& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& \Delta t\\ 0& 0& 0& 1\end{array}\right]$$E89

We also assume that the internal noise of the process has additive character and it is the source of change of motion velocity. In such a case, the velocity variation is increasing linearly with time:

$$Q_v=Q\left[\begin{array}{cc}\Delta t& 0\\ 0& \Delta t\end{array}\right]$$E90

The internal noise related to velocity affects the position through the dynamic process described by matrix$$F\left(\Delta t\right)$$. It can be proved that in this case the covariance matrix $$Q\left(\Delta t\right)$$ of the process internal noise has the following form:

$$Q\left(\Delta t\right)=Q\left[\begin{array}{cccc}\frac{\Delta t^3}{3}& \frac{\Delta t^2}{2}& 0& 0\\ \frac{\Delta t^2}{2}& \Delta t& 0& 0\\ 0& 0& \frac{\Delta t^3}{3}& \frac{\Delta t^2}{2}\\ 0& 0& \frac{\Delta t^2}{2}& \Delta t\end{array}\right]$$E91

The parameter $$Q$$ provides the information on how much the velocity of the motion changes during one second in the medium square sense. In the simulations we assume that$$Q=1$$.

Additionally we assume that the observation vector has the following form:

$$z=\left[\begin{array}{c}x\\ y\end{array}\right]$$E92

Therefore, the observation matrix $$H$$ has a form:

$$H=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right]$$E93

The direction finder, as well as the radar, make measurements at spherical coordinate system. To present measurements results at the global reference system we use formulas (63-65) described in the first example. The necessary elements have been already discussed. The algorithm of the estimation of the motion parameters can be described on the basis of measurements coming from different sources. Utilizing previous measurements we know the estimation $${\widehat{x}}_0$$ of the tracked target motion parameters in time $$t_0$$ and the covariance matrix$$p_0$$, describing the precision of this estimate. When at the time $$t_1>t_0$$ the next measurement occur, the following steps are taken:

• Calculate the matrices $$F\left(\Delta t\right)$$ and $$Q\left(\Delta t\right)$$ according to formulas (89, 91)

• Calculate the prediction of the state vector $${\widehat{x}}_{t1|t0}$$ using formula (35)

• Calculate the covariance matrix $$P_{t1|to}$$ of the prediction of the state using formula (36)

• Calculate the prediction of detection $${\stackrel{⌢}{z}}_{t1|t0}^{XYZ}$$ at the global coordinate system

• Calculate the prediction of detection $${\stackrel{⌢}{z}}_{t1|t0}^{Rad}$$ at the local radar (direction finder) coordinate system

• Expand the measurements vector$${\tilde{z}}_{t1}^{Rad}$$. The values which has not been measured are replaced by predicted values from the previous measurements.

• Create the precision matrix $$W_{t1}^{Rad}$$ of the measurements at the local radar (direction finder) coordinate system. The values which have not been measured have zero precision.

• Transform the expanded vector of measurements $${\tilde{z}}_{t1}^{Rad}$$ and corresponding precision matrix $$W_{t1}^{Rad}$$ to global coordinate system $${\tilde{z}}_{t1}^{XYX}$$ and$$W_{t1}^{XYZ}$$.

• Using formula (39-40, 46) calculate the new estimate of the motion parameters $${\widehat{x}}_{t1}$$ and their corresponding precision matrix$$P_{t1}$$.

Note that tracking is done in two dimensional space. The altitude Z is not tracked and is fixed as Z=2000. If the altitude is required in computation like for example in formula (97) we take Z=2000m.

For the first measurements we have accordingly:

$$\Delta t=108-100=8$$E94

$${\widehat{x}}_{108|100}=F(8){\widehat{x}}_{100}\approx \left[\begin{array}{c}754\\ 127\\ 59161\\ -114\end{array}\right]$$E95

The covariance matrix of the prediction:

$$P_{108|100}=F{\left(8\right)}^T{P}_{100}F\left(8\right)+Q\left(8\right)$$E96

The prediction of the detection at the global coordinate system XYZ:

$${\stackrel{⌢}{z}}_{108|100}^{XYZ}=\left[\begin{array}{c}754\\ 59161\\ 2000\end{array}\right]$$E97

In the last formula we use information on tracking target altitude 2000m. Now, we consider the system of the radar R1. The prediction of the detection at the radar local coordinate system:

$${\stackrel{⌢}{z}}_{108|100}^{Rad}=\left[\begin{array}{c}59200\\ 0.013\\ 0.034\end{array}\right]$$E98

The radar measures all coordinates. It is not necessary to increase measuring vector and to use the prediction.

$${\tilde{z}}_{108}^{Rad}=\left[\begin{array}{c}59025\\ 0.117\\ 0.212\end{array}\right]$$E99

The dispersion of the radar measurements is described by the precision matrix of the form:

$$W_{108}^{Rad}=\left[\begin{array}{ccc}1/{55}^2& 0& 0\\ 0& 1/{0.8}^2& 0\\ 0& 0& 1/{0.15}^2\end{array}\right]$$E100

After transformation to the global coordinate system XYZ we obtain:

$${\tilde{z}}_{108}^{XYZ}=\left[\begin{array}{c}6732.7\\ 57306\\ 2000\end{array}\right]$$E101

The precision matrix is transformed according to formula:

$$W_{108}^{XYZ}=J_{\widehat{z}}^{-1}{}^T{W}_{108}^{Rad}{J}_{\widehat{z}}^{-1}$$E102

where $$J_{\widehat{z}}^{}$$ is Jacobian matrix of transformation (66) calculated at the point$${\widehat{z}}_{108}^{Rad}$$.

Because we are tracking only the components XY, we ignore the last row in the measurement vector$${\tilde{z}}_{108}^{XYZ}$$, as well as the last row and the last column in the precision matrix$$W_{108}^{XYZ}$$.

Now, we can determine Kalman gain $$K$$ and the state estimate$${\widehat{x}}_{108}$$, as well as the corresponding covariance matrix$$P_{108}$$:

$$K_{108}=P_{108|100}{H}^T{W}_{108}{\left[W_{108}+{\left[H{P}_{108|100}{H}^T\right]}^{-1}\right]}^{-1}{\left[H{P}_{108|100}{H}^T\right]}^{-1}$$E103

$${\widehat{x}}_{108}={\widehat{x}}_{108|100}+K_{108}\left({\tilde{z}}_{108}^{XYZ}-{\widehat{z}}_{108|100}\right)$$E103.1

$$P_{108}=\left[1-K_{108}H\right]P_{108|100}$$E105

It is the end of the first step of the determination of the estimate of tracking target motion parameters.

The next step is more interesting because it demonstrates how one can use the measurements which were done by the direction finder. Accordingly:

For the first measurement we have:

$$\Delta t=113.8-108=5.8$$E106

$${\widehat{x}}_{113.8|108}=F(5.8){\widehat{x}}_{108}$$E107

The covariance matrix of prediction is:

$$P_{113.8|108}=F{\left(5.8\right)}^T{P}_{108}F\left(5.8\right)+Q\left(5.8\right)$$E108

The prediction of the detection at the global coordinate system XYZ

$${\stackrel{⌢}{z}}_{113.8|108}^{XYZ}=\left[\begin{array}{c}1571.7\\ 56647\\ 2000\end{array}\right]$$E109

We are taking into consideration the direction finder N2. The prediction of the detection at the direction finder local coordinate system:

$${\stackrel{⌢}{z}}_{113.8|108}^{N2}=\left[\begin{array}{c}18775\\ -0.4687\\ 0.1015\end{array}\right]$$E110

The direction finder measures only azimuth$$z_{N2}=\left[-0.3853\right]$$. We treat the direction finder as a radar which measures all coordinates. The missing coordinates are obtained utilizing the prediction$${\stackrel{⌢}{z}}_{113.8|108}^{N2}$$. We obtain:

$${\tilde{z}}_{113.8}^{N2}=\left[\begin{array}{c}18775\\ -0.3853\\ 0.1015\end{array}\right]$$E111

The precision matrix variables calculated using the prediction have zero precision:

$$W_{113.8}^{N2}=\left[\begin{array}{ccc}0& 0& 0\\ 0& 1/{0.05}^2& 0\\ 0& 0& 0\end{array}\right]$$E112

After the transformation to global XYZ coordinate system we obtain:

$${\tilde{z}}_{113.8}^{XYZ}=\left[\begin{array}{c}2987.1\\ 57291\\ 2000\end{array}\right]$$E113

As before, we determine the precision matrix $$W_{113.8}^{XYZ}$$ at the global XYZ coordinate system and the Kalman gain $$K_{113.8}$$ as well as new estimate $${\widehat{x}}_{113.8}$$ together with the corresponding covariance matrix$$P_{113.8}$$:

$$W_{113.8}^{XYZ}=J_{\widehat{z}}^{-1}{}^T{W}_{113.8}^{N2}{J}_{\widehat{z}}^{-1}$$E114

$$K_{113.8}=P_{113.8|108}{H}^T{W}_{113.8}{\left[W_{113.8}+{\left[H{P}_{1138|108}{H}^T\right]}^{-1}\right]}^{-1}{\left[H{P}_{113.8|108}{H}^T\right]}^{-1}$$E115

$${\widehat{x}}_{113.8}={\widehat{x}}_{113.8|108}+K_{113.8}\left({\tilde{z}}_{113.8}^{XYZ}-{\widehat{z}}_{113.8|108}\right)$$E116

$$P_{113.8}=\left[1-K_{113.8}H\right]P_{113.8|108}$$E117

It is obvious that in the case of the second bearing we act similarly:

$$\Delta t=114-113.8=0.2$$E118

$${\widehat{x}}_{114|113.8}=F(0.2){\widehat{x}}_{108}$$E119

The covariance matrix of prediction:

$$P_{114|113.8}=F{\left(0.2\right)}^T{P}_{113.8}F\left(0.2\right)+Q\left(0.2\right)$$E120

The prediction of the detection in global XYZ coordinate system:

$${\stackrel{⌢}{z}}_{114|113.8}^{XYZ}=\left[\begin{array}{c}2087\\ 56625\\ 2000\end{array}\right]$$E121

Now, we consider the local coordinate system of the direction finder N3. The prediction of the detection at the direction finder local coordinate system:

$${\stackrel{⌢}{z}}_{114|113.8}^{N3}=\left[\begin{array}{c}28750\\ -1.338\\ 0.065\end{array}\right]$$E122

The direction finder measures only azimuth$$z_{N3}=\left[-1.2433\right]$$.

$${\tilde{z}}_{114}^{N3}=\left[\begin{array}{c}28750\\ -1.2433\\ 0.065\end{array}\right]$$E123

The precision matrix variables defined using the prediction have zero precision:

$$W_{114}^{N3}=\left[\begin{array}{ccc}0& 0& 0\\ 0& 1/{0.05}^2& 0\\ 0& 0& 0\end{array}\right]$$E124

After transformation to global XYZ coordinate system we obtain:

$${\tilde{z}}_{114}^{XYZ}=\left[\begin{array}{c}2560.1\\ 59332\\ 2000\end{array}\right]$$E125

As before, we determine the precision matrix $$W_{114}^{XYZ}$$ at the global XYZ coordinate system and the Kalman gain $$K_{114}$$ as well as a new estimate $${\widehat{x}}_{114}$$ together with the corresponding covariance matrix$$P_{114}$$:

$$W_{114}^{XYZ}=J_{\widehat{z}}^{-1}{}^T{W}_{114}^{N3}{J}_{\widehat{z}}^{-1}$$E126

$$K_{114}=P_{114|113.8}{H}^T{W}_{114}{\left[W_{114}+{\left[H{P}_{114|113.8}{H}^T\right]}^{-1}\right]}^{-1}{\left[H{P}_{114|113.8}{H}^T\right]}^{-1}$$E127

$${\widehat{x}}_{114}={\widehat{x}}_{114|113.8}+K_{114}\left({\tilde{z}}_{114}^{XYZ}-{\widehat{z}}_{114|113.8}\right)$$E128

$$P_{114}=\left[1-K_{114}H\right]P_{114|113.8}$$E129

The presented example demonstrates the power of the described technique. Even a single bearing can improve the estimation of motion parameters of the tracking target. If several bearings are utilized, they need not to be simultaneous.

6. Summary

In the proposed technique the measurements carried out by different devices are treated in a simple and uniform manner. The key idea comprises in expansion the vector of the measured values and insertion of missing variables based on the prediction. The precision matrix is used to describe their dispersion. In order to obtain the estimate of the tracked target state we have to present all measurements and corresponding precision matrixes in a common coordinate system. The formulas (16-17, 18-21) should be used for static systems or the modified Kalman filters equations (27-40, 46) should be used for dynamic systems.

For linear objects, the Kalman filters are the optimal estimators. In the presented example dynamics of the tracked target has been described by linear differential equation The presented methodology can be applied also in the case when the dynamics of tracked target is nonlinear. In this case, instead of formula (27) we use linearization. Let the dynamics of the tracked process be described by a nonlinear function:

$$x\left(t+\Delta t\right)=F\left(x\left(t\right),\Delta t\right)$$E130

In order to determine the covariance of the state prognosis we use linearization $$\tilde{F}\left(x\left(t\right),\Delta t\right)$$ of the function $$F\left(x\left(t\right),\Delta t\right)$$

$$P\left(t+\Delta t\right)=\tilde{F}\left(x\left(t\right),\Delta t\right)P\left(t\right)\tilde{F}{\left(x\left(t\right),\Delta t\right)}^T+Q\left(\Delta t\right)$$E131

where:

$$\tilde{F}\left(t,\Delta t\right)=\left[\frac{\partial f_i}{\partial x_j}\right]$$E132

In order to increase accuracy, we can decrease time interval. Calculations (131) could be done several times. Let:

$$\delta t=\frac{\Delta t}{n}$$E133

$$x\left(t+\delta t\right)=F\left(x\left(t\right),\delta t\right)$$E134

$$x\left(t+\Delta t\right)=F\left(F\left(F\left(F\left(\mathrm{...}F\left(x\left(t\right)\right),\delta t\mathrm{...}\right),\delta t\right),\delta t\right),\delta t\right)$$E135

Analogically:

$$P\left(t+\delta t\right)=\tilde{F}\left(x\left(t\right),\delta t\right)P\left(t\right)\tilde{F}{\left(x\left(t\right),\delta t\right)}^T+Q\left(\delta t\right)$$E136

$$P\left(t+\Delta t\right)=\tilde{F}\left(x,\delta t\right)P\left(\tilde{F}\left(x,\delta t\right)P\left(\mathrm{...}\tilde{F}\left(x,\delta t\right)P\left(x\right)\tilde{F}{\left(x,\delta t\right)}^T+Q\left(\delta t\right)\mathrm{...}\right)\tilde{F}{\left(x,\delta t\right)}^T+Q\left(\delta t\right)\right)\tilde{F}{\left(x,\delta t\right)}^T+Q\left(\delta t\right)$$E137

In the past, this approach was rejected due to computer high power requirements. Presently it creates no problem at all.

An interesting case of proposed methodology application is the use of measurements from Doppler radar. The Doppler radar measures frequency shifting of received signal permitting to detect approaching speed of tracking targets. More precisely, it permits to determine speed radial component $$v_R$$ of the tracking target. The methodology is analogical to the example presented before. The velocity vector is expanded in order to include all components:

$$\tilde{v}=\left[\begin{array}{c}{v}_R\\ v_{\beta }\\ v_{\varphi }\end{array}\right]$$E138

The missing values are calculated utilizing the prediction. The corresponding precision matrix has the following form:

$$W=\left[\begin{array}{ccc}{\delta }_{vR}& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]$$E139

After transformation to the global coordinate system XYZ we can use measurements from the Doppler radar to improve estimation of the tracking target motion parameters.

In all so-far considered examples, we assumed that time error does not exist. In practice, the time is always measured with certain error. If for example the measuring time is presented in seconds, then during one second the tracked target can move several hundreds meters. Fortunately, in our case the error of time measuring can be easily taken into consideration. Time error $$\delta t$$ is transformed into the position measurements error $$\delta x$$ according to formula:

$$\delta x=vx*\delta t$$E140

Assuming that the time measuring error is independent from the position measurement error, the variances of both errors are added:

$$C_x^2\text{'}=C_x^2+{\delta }_x^2=C_x^2+v{x}^2{\delta }_t^2$$E141

The error of the time measurement can be easily taken into consideration by the proper increase of the error covariance value $$C_x^2$$ of the measurement component X. Of course, the same apply to the remaining components Y and Z.

When tracking the state of dynamic system by means of the modified Kalman filter the assumption that the measurements are done in regular time intervals has been abandoned. Nevertheless, it has been assumed that the received measurements are sequential:$$t_1\le t_2\le t_3\le \mathrm{....}\le t_{n-1}\le t_n$$. Sometimes even this assumption is unfulfilled. It may happen that due to a transition delay the measurements carried out at earlier time are received later then the ones obtained in later time. In order to consider the measurements delay the very simple but useful trick can be used. We store in memory the series of last measurements together with their estimated parameters values and the corresponding covariance matrices which describe accuracy of estimate. When the delayed measurement reaches the system we retrieve from the memory an estimation of time $$t_0$$ which occurred earlier then the time $$t_d$$ (attributed to the delayed measurement). All measurements later then $$t_0$$ and delayed measurements are rearranged to be sequential in time. Then, the prediction and estimation procedures are performed consecutively for each measurement. This way, the delayed measurement is taken into consideration. It must be noted that storing only the values of estimated parameters is not enough. In order to make the correct estimation, it is necessary to know also the covariance matrix of these parameter estimates.

The presented method is in 100% correspondence with the well-known least squares method (LSM). The careful analysis of the presented examples with taking into consideration the different weights of particular errors, existing correlations between them as well as the flow of time will provide the same results for the estimated values as the LSM method. Therefore the proposed methodology do note give better estimators then the classical theory. However the presented methodology avoids the very complicated and arduous computations and presents all measurements in a uniform way. Such approach is very useful for the automatic data processing by an automatic tracking system.

Since the presented technique is in 100% correspondence with the LSM, it also it discloses the weak points of LSM. In particular, it is not robust methodology. It is not robust for the values drastically odd (different) and for large errors. Even the single erroneous measurement not fitting the assumed model of errors distribution results in the drastic increase of the final error of the estimator. Because of this, it is recommended that the preliminary preselection of tracking measurements should be done. Drastically different values should be rejected.

The Kalman filter represents an optimal estimator for linear system only if the real model of the process corresponds to the one accepted during Kalman filter (27-34) design. It is particularly important when considering the error values. Consider for example the ladar (laser radar) which measures the position of an object from the distance of a few tens of kilometers with accuracy up to 5 m. Naturally, the error $${\delta }_R=5$$ should be taken into consideration. However, usually other sources of errors exist. For example, the position of ladar is known only with the limited precision, the time of measurement is provided with a limited accuracy. The Earth curvature can also be a source of errors during calculations. Based on the author’s experience some various processes which normally are not taken into consideration may be the source of casual and systematic errors. For this reason, asserting the error value $${\delta }_R=5$$ is a wrong idea. According to author it is better to use larger values of errors then the small ones. If erroneous values are assumed in the filter design and are larger then in reality the resulting filter is not optimal, but it is more flexible and robust with respect to drastically different values. For this reason, moderation and caution should be used in choosing the parameters of the model utilized to design a filter.