2.1

OLS – general case

The statistical approach to trilateration (multilateration) leads to the solution of the nonlinear least squares problem (NLS). One common approach to solve this problem is the linearization of nonlinear function. Given this approach, the linear location methodology attempts to transform nonlinear expressions into a set of the linear equations (1) with zero mean interference^2,4,8.

Such a system of linear equations can be written in matrix form as:

where the sizes of individual matrices and vectors: k – number of explained variables; n – number of observations (measurements), ε – random component. Later, β means the true (unknown) value of the parameter vector, and, a its estimation result.

The y vector is the observed (empirical) value of the explained variable, and ŷ is the theoretical value resulting from estimates:

The vector of residuals e is given by:

If the studied phenomenon can be described with a model (2), then unknown parameters can be obtained using the ordinary least squares method. The idea of this method is to find such values of an unknown parameter vector β that minimize the sum of squares of residuals, i.e. differences between the observed and theoretical values:

The sum of squared residuals (SSE) as a function of the searched parameter vector β is e^Te:

To find the that minimizes the sum of squared residuals, we need to take the derivative of (7) with respect to . This gives us the following equation:

Comparing the first derivative, i.e. (8) to zero, we what are called the normal equations:

Finally transforming (9) we get the expression describing the OLS method:

2.2

OLS – case study

For a specific analyzed case of tracking the location of an air object in 3D space, formula (10) takes the following form:

The individual matrices and vectors for the OLS method have the following forms: dP – position correction vector for individual object coordinates:

and matrix H:

where: k – discrete time index, R_i – distance between i – th modul and the localized object at the moment k – 1 estimated coordinates values of the tracked object, x_mi,y_mi, z_mi – coordinate values of i – th module.

In order to determine the best position estimate, the location correction vector of the tracked object was calculated (11). The elements of the dR vector correspond to subsequent differences between the measured values of the distance between and i – th module and the tracked object, and their estimated values:

The dP correction vector represents the estimation error of the object’s position parameters, therefore it is added to the current estimate and an iterative attempt to minimize its value is made:

The iteration shall be terminated when the following condition is met:

where: epsilon – assumed accuracy of position correction.

The advantage of this algorithm is the low mathematical complexity. On the other hand, it needs the appropriate selection of the starting point for the object position, several iterations in order to obtain an precise estimate of parameters and appropriate geospatial configuration of the modules.

2.3

EKF

Kalman filtration is widely well known so in the further part of the study only the elements necessary to describe and analyse the presented system will be cited.

The operation of each Kalman filtration algorithm is based on the estimation of a new process state in two consecutive steps^1,3,5,7: Prediction – the current process status is estimated based on the information from the previous step without taking into account any measurements and Correction – involves the use of new real measurements (noised) to correct the values determined in the previous step, at the prediction stage.

Regardless of the type of filter, its operation scheme (Fig. 4) can be presented in the form of several operations^1,5: initialization performed once at the beginning of the filter operation (x₀,P₀), and then recursively: prediction of the state vector (x) and the covariance matrix of prediction errors (P), calculation of Kalman’s gain matrix (K_k+1) and correction of prediction results using the last current measurement^3,6,9.

When developing the tracking system using Kalman filtration, it is necessary to define the state (dynamics) (Eq.17) and observation (Eq. 18) equations. Often, it turns out that one or both of them have a non-linear nature and then it is necessary to use a non-linear filtration algorithm.

Figure 4.

The Kalman filtration process⁹ – block diagram

In the system under consideration, it is assumed that both measurement errors (v) and process interferences (w) are Gaussian noises (taken into account at a later stage by the appropriate Q_k and R_k+1 matrices). The discreet version of the system model consists of the dynamics and observation equations^1,7:

Matrices Θ and H (transition and observation) are Jacobi matrices composed of partial derivatives of the nonlinear f and h functions with respect to all elements of the state vector.

In the case under consideration, we achieve the linear equation of system dynamics (state) and a non-linear equation of observation^1,7:

therefore it is only necessary to calculation the H matrix as partial derivatives.

The dynamics model adopted in the research looks as follows:

where: x, y, z – object coordinates in rectangular coordinate system, v_x, v_y, v_z – object velocity in rectangular coordinate system, w_x, w_vx, w_y, w_vy, w_z, w_vz – random disturbances of the motion model, T – discretization time.

The observation model^2,6, on the other hand, is the relationship between the measured distances in the system and the coordinates of the air object location in the state vector:

where: R – measured distance between the i-th module and the located object, X_i, Y_i, Z_i – coordinates of the i-th module, x_j, y_j, z_j – coordinates of the j-th point of the analysed trajectory.

3.1

Research organization

Simulation studies of the air object location system were carried out using UWB modules measuring distance. The configuration of the modules placement corresponds to the possibility of placing them in the MUT training area (Fig. 3) at a later time.

The tests were carried out for different configurations of the deployment and installation height of individual modules (installation in a short time in a simple way for high system mobility) and various trajectories of the air object movement.

The images presented later show: trajectory I (with a shape similar to the number eight) and trajectory II (showing the movement of the object outside the area covered by the system placement).

Table 1.

The position configuration of ultrawideband modules in 3D space.

In further analysis of distance data, the ordinary least squares method – OLS (as a reference algorithm), a tracking algorithm using extended Kalman filtration – EKF and its modification were used for the location of the object. This modification consisted of filter initialization with data from the approximation of the first trajectory point obtained using the OLS (EKF_OLS) algorithm.

3.2

Tests results – graphic imaging

The tests were carried out for two types of trajectories generated simulation in the Matlab® environment. Graphic images of the constant measured distance spheres for a circular system and their projections on the XY plane were presented for them. Then the ideal and estimated (using the OLS, EKF and EKF_OLS algorithms) trajectories were presented, as well as the estimated values of individual coordinates of the tracked air object. All these views are used to analyze the operation of the methods used to track the air object.

Figure 5.

The spheres of equal distance circular system for each modules and trajectory I : 3D view (left) and zoom with the exemplary point of the estimated position (right)

Figure 6.

The spheres of equal distance circular system for each modules and trajectory I: projection view on the XY plane (left) and zoom with the exemplary point of the estimated position (right)

Figure 7.

Trajectory I: ideal (blue), OLS estimated (red) and EKF estimated (green)

Figure 8.

Trajectory I: projection view on the XY plane (left) and subsequent values of the estimated height of the object (right)

Figure 9.

Trajectory I: subsequent estimated values of the X coordinate (left) and subsequent estimated values of the Y coordinate (right)

Figure 10.

The spheres of equal distance circular system for each modules and trajectory II: 3D view (left) and zoom with the exemplary point of the estimated position (right)

Figure 11.

The spheres of equal distance circular system for each modules and trajectory II: projection view on the XY plane (left) and zoom with the exemplary point of the estimated position (right)

Figure 12.

Trajectory II: ideal (blue), OLS estimated (red) and EKF estimated (green)

Figure 13.

Trajectory II: projection view on the XY plane (left) and subsequent values of the estimated height of the object (right)

3.3

Tests results – validation

In order to verify the implementation of each algorithms^8,11, depending on the type of air object trajectory (moving inside or outside the system’s operating area) and spatial configuration of the modules (changes in the position and height of mounting the modules), a quantitative assessment was also made using two types of statistical metrics^3,4: root mean square error (25) and mean absolute error (26) of the estimated spatial location of the tracked air object:

where: l – implementation number; Alg – type of algorithm (OLS, EKF lub EKF_OLS); NoP – Number of Points, j – j-th point trajectory, x_j,y_j,z_j – ideal coordinates, x_Alg,j,y_Alg,j,z_Alg,j – coordinates estimated in the selected algorithm, NoR – Number of Realizations.

The values of the respective tracking algorithms errors were calculated by averaging them for NoR = 100 implementations for different realizations of the measurement process disturbances (Table 2 and 3).

Table 2.

Averaged error values (m) for trajectories I (T) and various module configurations (C)

Table 3.

Averaged error values (m) for trajectories II (T) and various module configurations (C)

The tables show the cases where the EKF initialization data are correctly selected and the cases of incorrectly chosen values (marked as EKF_incorrect) which may cause deterioration of location parameters or the inability to determine the location of the object.