Method for Remote Determination of Object Coordinates in Space Based on Exact Analytical Solution of Hyperbolic Equations

Accurate remote determination of the object coordinates in 3D space is one of the main questions in many applications. In one of the most popular methods, such determination of the location of an object uses the measurement by receiving an electromagnetic signal transmitted by several spatially distributed base stations (BS). The main problem is that it is necessary to reduce errors and computation time. To overcome these difficulties, an analytical method for determining the position of an object based on the analysis of time difference of arrival (TDoA) of signals from the transmitter of the object to the receivers of the BS is proposed. One of the main advantages of this method is that it is possible to eliminate the ambiguity in determining the coordinates of the object in space and to increase the accuracy of determining the coordinates when the TDoA measurement between base stations fluctuates. Applications for autonomous automotive vehicles and space-based positioning systems are analyzed. The results obtained show that the proposed algorithm has an accuracy of determining coordinates several times higher than the method of linearization of hyperbolic equations and is less sensitive to TDoA fluctuations at base stations.


Introduction
Many modern applications need accurate 2D/3D object coordinates, which are determined by remote sensing methods. These applications include, for example, driving unmanned vehicles or security systems in "smart city" projects [1,2]. In one of the most popular methods, such determination of the location of an object uses measurements by receiving an electromagnetic signal transmitted by several spatially distributed base stations, taking into account shapes and other properties of objects [3]. Commonly used measurements are time delay for active scenarios [4] and time differences for passive environments [5]. The main aim of positioning tasks is to determine the coordinates (position) of the object/target, which can be carried out using various algorithms [6,7] and the most popular of them are discussed as follows.
For example, in [8], a control algorithm for an unmanned aerial vehicle (UAV) to circumnavigate an unknown target at a fixed radius when the UAV is unable to determine its location based on measurement errors and modeling the system by stochastic differential equations. Among other algorithms, it can be pointed out that with RSSI (Received Strength Signal Indication), the distance to the object is estimated by the signal power, but this algorithm usually has successful application for investigated in [58]. For the proposed new multi-target detection method for FMCW radar, the effect of the technical parameters of the vehicle radars on the required signal-to-noise ratio (SNR) of the receiver is estimated to ensure the probability of true determination of target parameters at 98%.
An exact analytical solution of equations for the task of localization of an object using a number of sensors is often challenged by outlier observations when the number of TDoA measurements is equal to the number of unknown transmitter coordinates [59][60][61]. The measurement equations are nonlinear, and an exact algebraic solution exists only in the specific case of a predetermined scenario in which the number of measurements is equal to the number of unknowns. The solution in [60] cannot use TDoA measurements with additional base stations improving positioning accuracy, and it has a high computational requirement for solving nonlinear equations. Compared to the solutions in [60], the solution proposed in [61] is more computationally attractive, requires only the roots of a quadratic equation, is more general and robust, and can work with arbitrary sensor arrangements.
A method proposed in [62] is based on Spherical-Intersection using only the root of a quadratic equation, but this method is not robust and fails to produce a reasonable solution for some sensor arrangements, due to the need to invert a matrix with the coordinates of sensors. Exact solutions can also be found in [63][64][65], and authors sometimes formulated an appropriated optimization task in form of a polynomial system, the solution of which was found by numerical algebraic methods or polynomial continuation techniques.
A good description of analytical methods is presented in [66][67][68]. Transformation of the coordinates of a hyperbola during a shift and rotation on a plane can be found in [66]. According to this transformation in [67], the effectiveness of two methods of TDOA analytical solution was compared based on the task of finding intersection points of hyperboloids (possible positions of a target). The first method analyzed was based on coordinate transformation from the initial system to a new system to simplify equations solving, and the second one was based on matrix. In [68], the results of an experimental verification of an analytical method based on coordinate transformation are presented. Finally, in [69], a high-precision analytical TDoA algorithm for determining the coordinates of an object on a plane (2D model) is proposed for eliminating the ambiguity of coordinates determination. Of great interest are recent works especially for indoor and AoA localization [70][71][72][73].
In this paper, an analytical method for determining the position of a target based on the analysis of TDoA of microwave radar signals from the transmitter to Base Stations (BS) receivers is proposed. The main feature of the method is that it can eliminate the ambiguity in determining of 3D coordinates of a target and improves the accuracy of determining coordinates when the TDoA measurements fluctuate on BS.
The materials of the article are presented in the following form: In Section 2.1, the fundamentals of the linearization method are presented without going into the specifics of solving linear equations by matrix methods. Section 2.2 sets out the proposed analytical method for determining the coordinates of an object in a positioning system. Section 3.1 describes the spatial ambiguity problem inherent in simple analytical methods. A comparison of the accuracy of determining coordinates by the linearization method and the proposed method at various levels of TDoA fluctuations for local and global positioning systems are carried out in Section 3.2. The results are discussed in Section 4. The paper ends with conclusions (Section 5) from the work done.

Linearization Method of Hyperbolic Equations
The method of linearizing hyperbolic equations was first proposed in [47]. The widespread use of this method at the present time is explained by the fact that it allows one to avoid solving nonlinear equations, which means that significant computing power of the positioning system processor is not required. The method consists in transforming nonlinear equations into a set of linear equations Ax = b in matrix form and further solving the system of equations by matrix methods [25,46,47,50,61,65]. We will compare our proposed algorithm (Section 2.2) with the linearization method; therefore, we will briefly outline the basics of the linearization method without delving into the specifics of solving linear equations by matrix methods.
In the absence of errors in measuring TDoA and receiving line-of-sight (LOS) signals, the real values of the difference in the arrival times of the TDoA signal between BS i and BS 1 are determined by the expression: where M-amount of BS, c l -speed of light, r i = (x i − x) 2 + (y i − y) 2 + (z i − z) 2 -distance between object and BS i , r 1 -distance between object and BS 1 , [x i , y i , z i ]-BS i coordinates, [x, y, z]-object coordinates. From r i1 = r i − r 1 follows r i1 + r 1 = r i . Substitution of the coordinates of the Cartesian system in Equation (1), squaring the right and left sides of the equation, expansion in a Taylor series, and linearization leads to the equation: which in matrix form has the form Aθ = b, where The superscript [ ] T denotes a matrix transposition operation. The solution of the equation Aθ = b can be carried out in various ways of solving the system of linear equations. After calculating the values of the vector θ, the coordinates of the target are determined x = θ 1 + x 1 ; y = θ 2 + y 1 ; z = θ 3 + z 1 ; where x 1 , y 1 , z 1 -known coordinates of BS 1 .

The Proposed Analytical Method for Determining the Coordinates of an Object in a Positioning System with 5 Base Stations
The variant of building a positioning system in space with 4 BSs is the most economically profitable; however, it has a disadvantage that the intersection of two hyperboloids can occur along two lines in space. In this case, it becomes impossible to identify the true position of the object. To eliminate the ambiguity in determining the position of the object, it is necessary to include the fifth BS in the positioning system. With the addition of the fifth BS, it becomes possible to solve four systems of equations, which will have one common root corresponding to the true position of the target. An analytical solution to the problem of determining the location of an object in 3D for a positioning system of 4 BSs was obtained in [67], but in [67], not all values of the coefficients are given, which complicates the direct use of the algorithm. In addition, in [67], the problem of eliminating ambiguity in determining the coordinates was not solved.
In what follows, the development of an analytical algorithm that eliminates the spatial zones of ambiguity in determining the coordinates of the target in space with high accuracy coordinates determination of the object is presented.
Let us consider a spatial positioning system that includes five microwave signal receivers BC C , BS L , BS R , BS U , BS D (Figure 1).

of 18
In what follows, the development of an analytical algorithm that eliminates the spatial zones of ambiguity in determining the coordinates of the target in space with high accuracy coordinates determination of the object is presented.
Let us consider a spatial positioning system that includes five microwave signal receivers BС С , BS L , BS R , BS U , BS D (Figure 1). The coordinates of the transmitter position are solutions to the system of hyperbolic equations: where ∆τ LC , ∆τ RC , ∆τ UC , ∆τ DC -TDoA of the microwave signal between BS L,R,U,D and reference station BS С .
After substitution of X = x − xC, Y = y − yС, Z = z − zС the system of equations can be rewritten as: The coordinates of the transmitter position are solutions to the system of hyperbolic equations: where ∆τ LC , ∆τ RC , ∆τ UC , ∆τ DC -TDoA of the microwave signal between BS L,R,U,D and reference station BS C . After substitution of X = x − x C , Y = y − y C , Z = z − z C the system of equations can be rewritten as: where K > 0: Sensors 2020, 20, 5472 6 of 18 The system of Equations (5) can be rewritten as: Squaring and reducing the general terms leads to the form: where In matrix form, the four systems of equations are: It should be noted that when using 4 BSs in the positioning system, there will be only one system of equations.
Solution of the first system of equations (i = 1) is: Sensors 2020, 20, 5472 7 of 18 The general determinant of the first system of equations is the following: The solutions of the second (i = 2), third (i = 3), and fourth (i = 4) systems of equations are obtained from the solution of the first system by successive replacement of the notation (Table 1): Table 1. Replacement of the notation.
Substitution of Equations (10)- (12) into Equation (6) defines a quadratic equation with respect to the variable K: a i K 2 + b i K + c i = 0, where i corresponds to the choice of a pair of hyperboloids, i = 1, 2, 3, 4.
The roots of the quadratic equation are the following (b i 2 − 4a i c i ≥ 0): In the case that one of the two roots K i1 and K i2 is negative for the same value of i, it can be immediately excluded from the solution, and then another root of the equation remains in the algorithm. However, a situation is possible when both roots K i1 and K i2 are positive. In this case, it becomes Sensors 2020, 20, 5472 8 of 18 impossible to determine the coordinates of the object. It is for such a fairly common case that the fifth BS has to be used in the 3D positioning system based on the analytical method.
Substitution of roots K i1 and K i2 in Equations (10)-(12) defines eight possible sets of x; y; z coordinates of an object, defining eight points in space: In the absence of TDoA measurement errors from eight calculated sets of x; y; z -coordinates, the coordinates of four points will be the same. It is this decision that is the true decision. To identify it, we carry out the following steps of the Algorithm 1: Calculation of 48 = (8 points × 6 combinations) distances D between points-candidates for the true value, determined by Equation (16): where both indices i and j take the values 1, 2, 3, 4, and i j.

2.
Sorting the elements of a string consisting of six distances related to one candidate point in ascending order.

3.
Summing the first three elements of the row. 4.
For all eight candidate points, writing the sum of the first three elements to a line. The meaning of each sum is the sum of the distances between the three nearest points.

5.
In parallel with step 4 of the algorithm, forming three lines by x; y; z from eight elements in each line. Each of the eight line elements corresponds to the coordinates of the candidate points. 6.
Making a cycle on 5 configurations of the spatial arrangement of BS. That is, each BS once becomes a reference. For each configuration, its own line of eight "sums of minimum distances" is formed. The strings are combined into a matrix. 7.
In the resulting matrix of "sums of minimum distances", seeking the element with the minimum value and its indices in the matrix. 8.
In parallel with item 7, forming three coordinate matrices from the rows of item 5. Each row of the matrix corresponds to the configuration of the location of the BS. 9.
Selecting, from the matrix of item 8 by the indices of item 7, the coordinates of the estimation of the position of the object x; y; z. The selected coordinates are the closest to the true coordinates of the target. The rest of the coordinates are either false or determine the coordinates of the target with a greater error.

Spatial Ambiguity Problem
In [67], an example of calculating target coordinates is given, in which the values of K i1 and K i2 are both positive. Two variants of target coordinates were obtained. The authors chose one of the options based on the disposition of the system of receivers to determine final target coordinates; however, they did not specify what kind of disposition it is. Most likely, it meant the negative z coordinate of the target, which means height. It should be noted that K i1 and K i2 are both positive and often occur with a positive z coordinate. In this case, the algorithm [56] is generally unable to give the correct values of the target coordinates. The method of linearizing hyperbolic equations also leads to a gross error in determining the coordinates of the object. In this article, a method is proposed that eliminates the ambiguity of determining the coordinates. Figure 2 shows the areas in which K i1 and K i2 are both positive, and therefore, there is no way to determine the true value of the target coordinates. The calculations were carried out at BS coordinates The method proposed by us is free from the indicated problem of ambiguity in determining the coordinates of an object. The method allows you to determine the coordinates of the object at all points in Figure 2.

Influence of the TDoA Fluctuations on the Accuracy of Coordinate Estimation
TDoA measurement errors occur mainly due to time desynchronization at base stations, as well as the conditions of propagation and reception of radio signals, and are the main reason for inaccurate determination of the object's position. Using a computer experiment in the LabVIEW environment, we investigated the dependence of the root mean square (RMS) deviation of the object coordinates on the standard deviation of the TDoA of microwave signals in the positioning radar. For this, Gaussian white noise with the same standard deviation value for all BSs was added to the TDoA value at the input of each receiver of the positioning system. To determine the standard deviation of the object position estimate, n = 500 computational experiments were carried out for each value of the standard deviation of the Gaussian white noise. The standard deviation of the object coordinates estimate was determined in accordance with the expression: where index k-experiment number, [ k x ; k y ; k z ] -estimation of object coordinates in k experiment, [x; y; z]-object's true coordinates. Figure 3 (time is in picoseconds) shows the dependences of the RMS deviation of the object coordinates estimate on the standard deviation of the Gaussian white noise of the difference in the The method proposed by us is free from the indicated problem of ambiguity in determining the coordinates of an object. The method allows you to determine the coordinates of the object at all points in Figure 2.

Influence of the TDoA Fluctuations on the Accuracy of Coordinate Estimation
TDoA measurement errors occur mainly due to time desynchronization at base stations, as well as the conditions of propagation and reception of radio signals, and are the main reason for inaccurate determination of the object's position. Using a computer experiment in the LabVIEW environment, we investigated the dependence of the root mean square (RMS) deviation of the object coordinates on the standard deviation of the TDoA of microwave signals in the positioning radar. For this, Gaussian white noise with the same standard deviation value for all BSs was added to the TDoA value at the input of each receiver of the positioning system. To determine the standard deviation of the object position estimate, n = 500 computational experiments were carried out for each value of the standard deviation of the Gaussian white noise. The standard deviation of the object coordinates estimate was determined in accordance with the expression: where index k-experiment number, [x k ;y k ;z k ] -estimation of object coordinates in k experiment, [x; y; z]-object's true coordinates.   From the simulation results, it follows that the proposed algorithm provides the accuracy of determining the coordinates of the object 10 times higher than the method of linearization of hyperbolic equations. The proposed method works with high accuracy (5 cm) with BS standard deviation of TDoA up to 200 ps, while the linearization method with such BS standard deviation of TDoA provides an accuracy of 60 cm. Gross failures in the linearization algorithm start at 500 ps standard deviation of TDoA, while the proposed analytical algorithm works without gross failures at 1000 ps standard deviation of TDoA.
The following RMS dependences of the coordinate estimates were obtained for the configuration of five BSs located on a circle as in [74]  The following RMS dependences of the coordinate estimates were obtained for the configuration of five BSs located on a circle as in [74] with a radius of 50 m, which corresponds to the Local Positioning System (LPS). Coordinates were as follows: BS C [x C ; y C ; z C ] = [29.4; −40.45; 20]  It should be noted that the algorithm with one equation [56], as well as the linearization method, produces a gross error at the points of the 3D space at which the determinant of the system of equations becomes zero. An example of such a situation is shown in Figure 5.  It should be noted that the algorithm with one equation [56], as well as the linearization method, produces a gross error at the points of the 3D space at which the determinant of the system of equations becomes zero. An example of such a situation is shown in Figure 5. An effective way to eliminate such catastrophic errors is to use a cycle of five configurations of the spatial arrangement of BS, that is, to replace the reference BS. If a configuration with some reference BS has the effect shown in Figure 5, then when changing the reference BS, this effect will no longer be present. This method is implemented in clause 6 of the proposed algorithm. All graphs below are calculated using a cycle for five configurations of the spatial arrangement of the BS ( Figures  6 and 7). An effective way to eliminate such catastrophic errors is to use a cycle of five configurations of the spatial arrangement of BS, that is, to replace the reference BS. If a configuration with some reference BS has the effect shown in Figure 5, then when changing the reference BS, this effect will no longer be present. This method is implemented in clause 6 of the proposed algorithm. All graphs below are calculated using a cycle for five configurations of the spatial arrangement of the BS (Figures 6 and 7).
An effective way to eliminate such catastrophic errors is to use a cycle of five configurations of the spatial arrangement of BS, that is, to replace the reference BS. If a configuration with some reference BS has the effect shown in Figure 5, then when changing the reference BS, this effect will no longer be present. This method is implemented in clause 6 of the proposed algorithm. All graphs below are calculated using a cycle for five configurations of the spatial arrangement of the BS ( Figures  6 and 7).  The RMS area distributions when Gaussian noise is added to the TDoA values with a standard deviation of more than 10 ps have a similar form to Figure 7, and the differences lie in the values of RMSmin and RMSmax, which are summarized in Table 2. RMSmin for the linearization method is approximately three times higher than for the analytical method; however, there are gross errors associated with zeroing the determinant of the system at some points in space. Therefore, RMSmin and RMSmax for the linearization method are not included in Table 1. As can be seen from the graphs, the accuracy of determining the coordinates strongly depends on the choice of the BS position.
Configurations of global satellite positioning systems are of practical interest. Let us apply the proposed algorithm to simulate the accuracy of determining coordinates in the global positioning system (GPS). Orbiting satellites of the system should be perceived as BS, and the object, the coordinates of which are determined, is on the Earth. We will accept the coordinates of the BS  The RMS area distributions when Gaussian noise is added to the TDoA values with a standard deviation of more than 10 ps have a similar form to Figure 7, and the differences lie in the values of RMS min and RMS max , which are summarized in Table 2. RMS min for the linearization method is approximately three times higher than for the analytical method; however, there are gross errors associated with zeroing the determinant of the system at some points in space. Therefore, RMS min and RMS max for the linearization method are not included in Table 1. As can be seen from the graphs, the accuracy of determining the coordinates strongly depends on the choice of the BS position.
Configurations of global satellite positioning systems are of practical interest. Let us apply the proposed algorithm to simulate the accuracy of determining coordinates in the global positioning system (GPS). Orbiting satellites of the system should be perceived as BS, and the object, the coordinates of which are determined, is on the Earth. We will accept the coordinates of the BS (satellites) as follows:  Figure 8. The target height is z = 0 m. The number of computational experiments n = 100. The satellites are equipped with atomic clocks with a daily instability of no worse than 10 −13 ; therefore, in the calculations, we limited ourselves to the maximum value of the standard deviation of the Gaussian noise TDoA 10 ps (Figures 9 and 10). Such a model can be viewed as a simplified model of the location of satellites and targets on Earth in a GPS positioning system.        RMS min and RMS max for the analytical method and linearization method for the last problem (GPS) are summarized in Table 3. It can be seen that the proposed algorithm is significantly superior to the widely used positioning algorithm based on the linearization of hyperbolic equations.

Discussion
The new method, based on the analytical solution of hyperbolic equations, in contrast to many existing methods, makes it possible to obtain an accurate solution for determining the coordinates of an object based on the TDoA measurement. Considerable efforts of researchers have been aimed at finding ways to transform expressions with radicals into linear relations and to ensure fast convergence of solutions. In part, this direction of research is explained by the limited capabilities of the processor technology and the complexity of the implementation of computations of expressions with radicals, which is necessary in the analytical exact method that we have proposed. Of course, the limited mathematical capabilities of inexpensive processors is the main limitation to the widespread implementation of the proposed method. However, the computational element base, both based on processors and FPGAs, is rapidly developing, and we should expect in the near future the appearance of inexpensive processors that allow calculating mathematical functions with radicals in real-time. In addition, in some positioning systems, the function of mathematical processing can be transferred to a system controller (server), in the software of which the proposed method is implemented. In these cases, the restriction on the possibilities of mathematical data processing is removed.
To accurately determine the moment of arrival of a radio signal from an object to base stations, various forms of radio signals are used. The wider the spectrum of the radio signal, the more accurately it is localized in time. Therefore, ultra-wideband (UWB) signals with a wide spectrum have advantages over signals with a shorter spectrum and are used in positioning systems [75]. In our article, the question of the form of the used radio signal is not a subject of research. Radio signal propagation conditions impose a significant influence on the functioning of positioning systems. Multipath propagation, broadband interference, and atmospheric phenomena in global positioning systems lead to TDoA fluctuations. In our work, we did not consider the issues of propagation and reception of radio signals. In our model, the standard deviation of the TDoA is used as a measure of fluctuations, the value of which is determined by the factors listed above.
Within the framework of the adopted model, we have established that the proposed method has the ability to exclude spatial zones of ambiguity in determining the coordinates, which are characteristic of approximate analytical algorithms with one system of equations. The algorithm also eliminates the gross error at the points of the 3D space, at which the determinant of the system of equations becomes zero.
In future work, we propose to calculate the Cramer-Rao Lower Bound (CRLB) for the proposed method, carry out practice implementations to compare the performance with the simulations results, investigate the resistance of the method to non-line-of-sight (NLOS) radio signal propagation, and also consider the applicability of the method for determining the coordinates of a set of objects.

Conclusions
Modern systems of local and global positioning put forward high requirements for the accuracy of determining the coordinates of objects. This article proposes an analytical method for the exact solution of a system of hyperbolic equations, which, with a high accuracy of less than a few cm, allows one to determine the coordinates of objects and at the same time exclude spatial zones in which approximate analytical algorithms with one system of equations are unable to unambiguously determine the coordinates of an object or give a gross error at the points of the 3D space, at which point the determinant of the system of equations becomes zero. The absence of any omitted terms in the solution provided the developed algorithm with significant resistance to fluctuations in TDoA caused by the conditions of propagation and reception of radio signals and time desynchronization of the BS while maintaining a high accuracy in determining the coordinates. It was found that the proposed method works with high accuracy (5 cm) with BS standard deviation of TDoA in time up to 200 ps, while the linearization method with such BS standard deviation of TDoA provides an accuracy of 60 cm for LPS. The performance of the algorithm has been confirmed by computational experiments for both LPS and GPS.