Fast and Fault-Tolerant Passive Hyperbolic Localization Using Sensor Consensus

TL;DR

This paper introduces a fast and reliable method for accurate location estimation in challenging environments using sensor consensus.

Contribution

A novel branch-and-bound computational method with proven global convergence for optimizing consensus functions in localization.

Findings

The proposed algorithm achieves global maximum of consensus functions efficiently.

Simulation and real measurements confirm the method's robustness and accuracy in NLOS conditions.

Abstract

The accuracy of passive hyperbolic localization applications using Time Difference of Arrival (TDOA) measurements can be severely compromised in non-line-of-sight (NLOS) situations. Consensus functions have been successfully used to provide robust and accurate location estimates in such challenging situations. In this paper, a fast branch-and-bound computational method for finding the global maximum of consensus functions is proposed and the global convergence property of the algorithm is mathematically proven. The performance of the method is illustrated by simulation experiments and real measurements.

Figures

Full text

1. Introduction

Determining the location of an object is a key service in a wide range of applications, including navigation [1,2], robotics [3,4], tracking [5,6], safety and security [7,8], healthcare [9,10], and defense [11,12], just to name a few. Global Navigation Satellite Systems (GNSSs) can be used in various outdoor applications [1], while for indoor and other GNSS-deprived areas, several other technologies have been proposed, based on, e.g., acoustic signals [8,11,12,13], radio signals [5,9], visible light [4,14], or images [15,16].

The localization systems utilize a number of spatially distributed anchors, and most often measure quantities directly related to distances, distance differences, angles, or angle differences between the anchors and the target. The location estimate is then calculated using Time of Arrival (TOA), Time Difference of Arrival (TDOA), angle of arrival (AOA), or angle difference of arrival (ADOA) methods, respectively [17].

The TDOA localization, also known as hyperbolic positioning, is a fairly common method in geolocation [1,18,19], sensor networks [8,11,20,21], indoor localization [22,23], or the Internet of Things [24,25]. In the TDOA scenario, differences of distances between the target and the anchors are used, most often calculated from differences of measured times of arrivals of signals (sound or electromagnetic waves) between the anchors and the target, using the propagation speed of the signal. The one-stage or Direct Position Determination (DPD) methods determine the position estimate without explicitly calculating the Time Difference of Arrival (TDOA) values. In contrast, the two-stage methods estimate the TDOA values in stage 1 and then estimate the location from them in stage 2. The location estimation process is also known as position fixing. Although two-stage methods are not optimal, they are asymptotically equivalent to DPD methods and are popular in systems where the transmission of large amounts of raw data between sensors and the central processing unit is prohibitive [26]. This paper focuses on the position fixing of two-stage methods.

The process of position fixing most often involves closed-form solutions [18,27,28,29,30,31,32], iterative techniques [33,34], or various consensus-based methods [35,36,37,38,39,40]. These methods yield satisfactory results when the time difference measurements contain only small errors (i.e., measurement noise). However, in practical scenarios, measurements can contain significant errors (outliers), especially when there is no direct line of sight between the anchors and the target, or when the environment produces signal reflections, or when signal detection, and thus time measurements, is not reliable. If such outlier measurements are present, the position estimates of closed-form and iterative techniques will be biased unless additional measures are taken to remove the outliers prior to the estimation [41,42]. Consensus-based methods can provide accurate estimates even when outlier measurements are present, and were successfully utilized in applications where non-line-of-sight conditions and reflections cause unreliable measurements [38,39].

In this paper, a consensus-based method is studied, where the position estimate is determined by identifying the location with the highest value of the consensus function [38]. A fast computation method was recently proposed that was observed to reliably find the global optimum of the consensus function [40]. This paper reviews the consensus function and introduces an updated search method. The contributions of this paper are the following:

• A new, fast evaluation method is introduced.
• The main contribution is the theoretical proof that ensures that the proposed algorithm always finds the global maximum of the consensus function over a finite grid.
• Finally, a comprehensive performance analysis is provided using simulations and real measurements.

The symbols used in this paper are enumerated in the Section: List of Symbols.

2. Related Work

In Section 2.1, the TDOA problem is formulated, with special emphasis on potential large measurement errors, e.g., due to non-line of sight (NLOS). The consensus-function-based solution is reviewed in Section 2.2, and various fast evaluation methods are discussed in Section 2.3.

2.1. The TDOA Localization Problem

The emitter $[eqn]$ is located at an unknown position $[eqn]$ and it emits a signal at an unknown time $[eqn]$ . Sensors $[eqn]$ are deployed at known positions $[eqn]$ . The emitted signal is detected by sensor $[eqn]$ at time instant $[eqn]$ . If the propagation speed of the signal is $[eqn]$ and there is line of sight (LOS) between $[eqn]$ and $[eqn]$ , then the detection time $[eqn]$ can be expressed as follows:

$$[eqn]$$

A measurement scenario is illustrated in Figure 1, where sensors $[eqn]$ are in LOS of the source, while $[eqn]$ and $[eqn]$ are in NLOS positions and provide measurements through a reflective surface. In the example, measured times $[eqn]$ satisfy (1), but for $[eqn]$ and $[eqn]$ , the measured time is larger than the ideal time calculated from (1).

The objective is to determine the position $[eqn]$ using the measurements $[eqn]$ and positions $[eqn]$ . The task becomes more challenging when measurements contain outliers, such as those caused by NLOS conditions.

Several solutions have been proposed to solve the TDOA problem. Closed-form solutions are available for a limited number of sensors [27,28,29,30,31,32], while iterative techniques are used for general cases [33,34]. Particle filters [35,43], Hough transform-based solutions [35,36,37], and consensus functions [38,39,40] have also been proposed.

In most TDOA localization problems, the propagation speed of the signal is assumed to be known. In [44], a joint estimation method is proposed, which is capable of estimating both the source location and the propagation speed. The placement of anchors has a significant effect on the localization accuracy. The optimal geometry, considering the communication constraints that are important in practice, has been analyzed in [45].

Most solutions assume that measurements are correct and do not include outliers. In practical cases, however, outliers are common, due to NLOS situations, reflections, sensor faults, etc. Therefore, in real-world scenarios, the presence of outliers can severely bias location estimates unless they are removed from the estimation [41,42].

Several works address the NLOS problem. In [46], an iterative solution, based on the maximum correntropy criterion with a variable center, was proposed. Various methods using convex optimization and semidefinite programming have also been proposed [47,48,49]. However, these methods are computationally intensive. In [50], a neurodynamic optimization approach is proposed, which provides a trade-off between accuracy and computational complexity.

Other approaches for providing accurate location estimates in the presence of outliers involve sensor consensus. RANSAC-based [51,52] methods randomly select a core group including a small number of sensors (typically 4–5) and use it to perform localization. The measurements from the other sensors are then checked against the core estimate, and the location estimate is refined. In an iterative process, several core groups are selected, and the best result is chosen. By conducting a high number of trials, the results are statistically correct with a high probability [53].

Hough transform-based solutions utilize a voting system in which each pair of sensors votes for a set of possible location estimates. The position with the highest number of votes becomes the location estimate. Fast calculation methods exist when the unknown location is in a plane [37].

Consensus functions also use a voting system, but unlike Hough transform and RANSAC, which use groups of a small number of sensors in the voting process, here, all of the sensors cast a single cooperative vote. The method was successfully utilized in three dimensional acoustic systems with a high number of unreliable measurements [38,39]. The detailed operation is reviewed in Section 2.2.

2.2. Consensus-Function-Based Localization

In this section, consensus functions are reviewed based on [38,39]. The distance $[eqn]$ between sensor $[eqn]$ and an arbitrary point $[eqn]$ is

$$[eqn]$$

If the source position is $[eqn]$ , then the emission time $[eqn]$ can be estimated from measurement $[eqn]$ using (1):

$$[eqn]$$

Thus, $[eqn]$ is the estimate of the emission time by a single sensor $[eqn]$ , provided the source is at location $[eqn]$ . Using all sensor measurements, the set of $[eqn]$   $[eqn]$ estimates are available for an arbitrary point $[eqn]$ . From the set of $[eqn]$ , the consensus group is formed, according to Definition 1.

Definition **1.**A consensus group at position  $[eqn]$  with window length  $[eqn]$  contains an S subset of  $[eqn]$  such that for any two elements  $[eqn]$  it is true that  $[eqn]$

In the ideal case when there are no measurement errors present, for the true location $[eqn]$ , all $[eqn]$ estimates are equal to the true emission time $[eqn]$ . In real cases when measurement noise is present, the $[eqn]$ estimates will be in the close vicinity of $[eqn]$ , forming a large consensus group. The extent of the ‘vicinity’ is dependent on the level of noise and characterized by the window length $[eqn]$ . Sensors that produce outlier measurements may have an estimated value $[eqn]$ that is significantly different from $[eqn]$ , and these estimates will not be members of the consensus group.

An illustration is shown in Figure 2a, where estimates $[eqn]$ form a five-element consensus group around time $[eqn]$ . Estimates $[eqn]$ and $[eqn]$ , corresponding to outlier measurements, are not part of the consensus group (they form one-element consensus groups).

Definition **2.***The value of the consensus function at point $[eqn]$ * is defined as the cardinality of the largest consensus group for location $[eqn]$ .

In the example of Figure 2a, the value of the consensus function is 5. For locations $[eqn]$ other than $[eqn]$ , estimates $[eqn]$ are likely to differ significantly, resulting in the formation of small consensus groups rather than one large group. Figure 2b shows an example where the consensus function is evaluated at a point other than the true source location. The estimates $[eqn]$ are scattered and the largest consensus groups contain only two sensors; thus, the value of the consensus function is 2.

Note that sensors corresponding to a consensus group calculated at point $[eqn]$ agree on the hypothesis that the source is located at position $[eqn]$ . The more sensors in consensus, i.e., the higher the consensus function, the higher the likelihood that the source is actually at position $[eqn]$ . Therefore, the estimated source location is the one with the highest consensus function.

To create consensus groups, the maximum difference $[eqn]$ allowed between group members is defined. In Figure 2, a sliding window with width $[eqn]$ is shown. The window function $[eqn]$ of width $[eqn]$ is defined as follows:

$$[eqn]$$

Using the window function, the consensus function can be expressed as follows:

$$[eqn]$$

The set of points where the consensus function $[eqn]$ takes its maximum is

$$[eqn]$$

The location estimate $[eqn]$ is calculated as the mean of the positions in $[eqn]$ :

$$[eqn]$$

The window width $[eqn]$ is essential in calculating the consensus function as it determines the required proximity of estimated emission times $[eqn]$ within a consensus group. At the true location, the estimates $[eqn]$ should be treated as a single group, despite the small perturbations in the values caused by measurement noise, as shown in Figure 2a. This requirement leads to a lower bound on $[eqn]$ , as follows. The accuracy of $[eqn]$ depends on the accuracy of measurements $[eqn]$ and the precision of the sensor locations $[eqn]$ . If the maximum sensor location error is $[eqn]$ and the maximum time measurement error is $[eqn]$ , then, from (3), the maximum error of $[eqn]$ is

$$[eqn]$$

Since the maximum difference between any two values of $[eqn]$ is $[eqn]$ , this is the lower bound on $[eqn]$ . To keep $[eqn]$ as small as possible, the window width is set to

$$[eqn]$$

In practical scenarios, a rough a priori estimate of the measurement errors $[eqn]$ and $[eqn]$ is required. From these estimates, the design parameter $[eqn]$ can be determined using (9).

2.3. Calculation of the Consensus-Function-Based Location Estimate

The source location is determined by finding the maximum of the consensus function. However, this can be challenging because the consensus function is not smooth, being an integer-valued function, and it may have several local maxima. Traditional gradient-based search methods do not work for such functions. An exhaustive search on a finite grid provides the global maximum; however, for large target areas, especially in 3D, it is not practical. To provide a faster solution, a fast heuristic calculation was proposed in [38,39]. This method is based on integer arithmetic and the Generalized Bisection method and searches for the source location and emission time in a four-dimensional space (x, y, z, t). The performance of the search method has been experimentally validated [38], but the correctness of the algorithm has not been proven.

Another fast solution was proposed in [53] using a RANSAC-based approach. Random sets of five measurements were utilized to calculate the core estimate, which was then compared to the other measurements using the consensus function. This method is able to provide the correct estimates with high probability if the number of random trials is high enough [53].

A fast Branch and Bound-type calculation method has been proposed recently [40]. A global convergence property of the method has been reported, supported by an experimental validation. The following section presents an improved version of [40] and provides a mathematical proof of its global convergence property.

3. Fast Consensus-Based Localization

The fast search method proposed in [40] performs the search on a grid for the maximum of the consensus function. The grid is iteratively refined until the required grid size is reached. During the refinement process, areas that are unlikely to contain the global maximum are excluded from a further search; thus, the search method converges quickly. The method’s correctness was validated through simulation experiments. Now, an enhanced and more detailed version of the algorithm [40] is proposed in Section 3.1. The correctness of the algorithm is proven in Section 3.2 and Section 3.3.

3.1. Fast Search on a Finite Grid

The proposed fast search method is described in Algorithm 1, followed by an illustration of its operation.Algorithm 1 Branch and Bound algorithm to find the maximum of the consensus function Input:

• 1.measurements: $[eqn]$
• 2.sensor positions: $[eqn]$
• 3.required (fine) grid size: $[eqn]$ Initialization:
• 4.Define a grid $[eqn]$ over the search area with a grid cell size equal to $[eqn]$ . This is the fine grid where the maximum of the consensus function is searched for.
• 5.Cover the search area with a coarse grid $[eqn]$ using grid cells of size $[eqn]$ where $[eqn]$ is an integer greater than 1, so that the grid points of $[eqn]$ overlap with some of those of $[eqn]$ . The search starts from this grid.
• 6.Mark each cell of $[eqn]$ as active.
• 7.Calculate the upper bound $[eqn]$ of $[eqn]$ for each cell Iteration: repeat Branch and Bound while active cells exist****Branch:
• 8.Select the active cell with the highest $[eqn]$ value and refer to it as $[eqn]$ . If multiple cells have the same highest value, select the one with the smallest size.
• 9.Replace $[eqn]$ with smaller cells $[eqn]$ that have a linear size half of the original ( $[eqn]$ in 2D and $[eqn]$ in 3D).
• 10.If $[eqn]$ is completely outside of the search area, mark it passive.
• 11.If the size of $[eqn]$ is larger than $[eqn]$ then mark $[eqn]$ as active and recalculate the upper bound $[eqn]$ for $[eqn]$
• 12.If the size of $[eqn]$ is equal to $[eqn]$ then mark $[eqn]$ as final and calculate $[eqn]$ in the center of the cell Bound:
• 13.Calculate the maximum $[eqn]$ of the $[eqn]$ values in the final cells.
• 14.Mark active cells with $[eqn]$ as passive Output:
• 15.Collect the final cells with $[eqn]$ to form $[eqn]$ $[eqn]$

The algorithm finds the maximum of the consensus function using a finite grid with a grid size of $[eqn]$ (step 3), i.e., in 2D, the grid cells are squares of size $[eqn]$ and in 3D, the cells are cubes of size $[eqn]$ . The consensus function is evaluated in the center of the cells. Note that the grid size determines the resolution of the search.

Figure 3 illustrates the operation, using a 2D example. Figure 3a shows the fine grid over the search space (step 4). Then, a coarse grid (in red) is placed over the fine grid (step 5). The size of the coarse grid must be a $[eqn]$ multiple of $[eqn]$ ; in the example $[eqn]$ , the coarse cells are thus of size $[eqn]$ , as shown in Figure 3b. Note that in practical cases, $[eqn]$ is typically much larger than the small value used in this example.

The further operation is illustrated in Figure 4. In Figure 4a, the three starting cells are active, indicated by the white background, and the upper bounds $[eqn]$ are displayed within each cell (according to steps 6–7). Figure 4b shows the first branching step. The highest upper bound is seven, so the corresponding cell is selected and divided into four smaller cells (steps 8–9). One cell is outside of the search region and is marked as passive (step 10), shown by an orange background. The new cell size is still larger than $[eqn]$ , so the upper bounds are calculated for the new cells (step 11). As there are no exact consensus values calculated yet, the bound step is idle. Figure 4c illustrates the subsequent branching step, where the cell with an upper bound of 6 is divided into four smaller cells. The new cell size is now $[eqn]$ , so the new cells are marked as final (shown by green background) and the exact consensus values are calculated in these cells (step 12). The next bound step is shown in Figure 4d. The value $[eqn]$ is 5 (step 13), and the maximum values are shown by a yellow color in the corresponding cells. Several cells with $[eqn]$ are marked as passive (step 14), shown by an orange background color. Figure 4e shows the next branching step of the active cell with an upper bound of 5. The subsequent bound step is shown in Figure 4f: $[eqn]$ is 5, and thus the cells with upper bound values below 5 are marked as passive. In Figure 4g, the next branch step is shown, where the cell with $[eqn]$ is divided into four final cells. Since there are no more active cells, the bound step is idle and the iteration is finished. Figure 4h shows the cells with $[eqn]$ , forming $[eqn]$ (step 15).

The design parameter $[eqn]$ defines the required resolution and thus accuracy, since the search method finds the maximum of the consensus function of the defined grid. Small $[eqn]$ results in a finer resolution but potentially requires more computation. It is therefore inadvisable to select values for $[eqn]$ that are excessively small.

In practical situations, the size of $[eqn]$ can be determined along with parameter $[eqn]$ , as follows: The window width $[eqn]$ is determined from the time and distance measurement errors, according to (9). The window width $[eqn]$ reflects the measurement errors in terms of timing error, while the quantity $[eqn]$ quantifies the measurement errors in terms of distance. Although the overall positioning accuracy that can be achieved is dependent on the placement of the sensors [44], the order of magnitude can be estimated by $[eqn]$ . Thus, a practically good design choice for grid size $[eqn]$ lies in the range of $[eqn]$ .

3.2. Upper Bound of the Consensus Function

In the operation of Algorithm 1, the upper bound $[eqn]$ of the consensus function in a cell is essential. In [40], a quantity was proposed for $[eqn]$ and it was validated through simulations. However, to date, no proof has been found to support that the proposed quantity is indeed an upper bound. In the following, a mathematical proof will be provided.

The following theorem, Theorem 1, provides an upper bound for the consensus function in a spherical region.

Theorem **1.***Let point $[eqn]$ * *be the center of a sphere $[eqn]$ * *with radius $[eqn]$ * *and let $[eqn]$ * *be a point inside $[eqn]$ * *. If $[eqn]$ *, then the following inequality holds:

$$[eqn]$$

Proof of Theorem **1.**Points $[eqn]$ , $[eqn]$ , and sensor $[eqn]$  are shown in Figure 5, along with the sphere $[eqn]$  around $[eqn]$ . The distance between points $[eqn]$  and $[eqn]$  is denoted as $[eqn]$ . Let us define the quantity $[eqn]$  as follows:

$$[eqn]$$

Since $[eqn]$ is inside $[eqn]$ , trivially

$$[eqn]$$

By applying the triangle inequality to the triangle $[eqn]$ of Figure 5, it also follows that

$$[eqn]$$

By combining (11)–(13), it follows that

$$[eqn]$$

Let the value of the consensus function $[eqn]$ at point $[eqn]$ be denoted by $[eqn]$ . Thus, according to (5) and (4), there exists a time instant $[eqn]$ and a set $[eqn]$ of $[eqn]$ sensors for which the following holds:

$$[eqn]$$

For the sake of simplicity, but without loss of generality, let us assume that the sensors in $[eqn]$ are $[eqn]$ . In the proof, index $[eqn]$ refers to these sensors, specifically $[eqn]$ . Using (3), the inequality (15) can be written as

$$[eqn]$$

Using (11), the following holds for point $[eqn]$ :

$$[eqn]$$

By utilizing the basic inequality $[eqn]$ from (17), it follows that

$$[eqn]$$

Applying (14) and (16), inequality (18) leads to

$$[eqn]$$

and thus using (3), it follows that

$$[eqn]$$

Comparing (4) and (20) leads to

$$[eqn]$$

Note that (21) holds for all $[eqn]$ sensors $[eqn]$ . Therefore, using (5) and (21), it follows that

$$[eqn]$$

Thus, for any point $[eqn]$ inside $[eqn]$ , $[eqn]$ , which leads directly to (10). □

Using the results of Theorem 1, Theorem 2 provides an upper bound for the consensus function in a cell.

Theorem **2.***Let us denote the dimension of the search space with $[eqn]$ **Let point $[eqn]$ * *be the center of a cell $[eqn]$ * *of size $[eqn]$ **, where $[eqn]$ * *is a cube for $[eqn]$ * *or a square for $[eqn]$ **. Then, for all points $[eqn]$ * *inside $[eqn]$ *,

$$[eqn]$$

with

$$[eqn]$$

Proof of Theorem **2.**For a cube of size $[eqn]$ , the maximum distance between its center $[eqn]$  and any point $[eqn]$  inside the cube is $[eqn]$ . For a square of size $[eqn]$ , the maximum distance between $[eqn]$  and any point $[eqn]$  inside the square is $[eqn]$ . Thus, for a cell of dimension $[eqn]$  ( $[eqn]$ ), the cell is inside of a sphere with center $[eqn]$  and radius $[eqn]$ . With the choice of

$$[eqn]$$

Theorem 1 yields (23). □

3.3. Global Convergence

Finally, the global convergence property of Algorithm 1 is proven.

Theorem **3.***If the upper bound $[eqn]$ * *for a cell with center $[eqn]$ * *and size $[eqn]$ * *is computed as $[eqn]$ * *with $[eqn]$ **, then Algorithm 1 finds the global maximum of $[eqn]$ * on the finite grid $[eqn]$ .

Proof of Theorem **3.**Let $[eqn]$  be the maximum of the consensus function inside the search area, calculated over $[eqn]$ , and let $[eqn]$  denote a particular cell of $[eqn]$  inside the search area, where the consensus function is $[eqn]$ . In Algorithm 1, the initial cells are iteratively divided into smaller cells in the branching steps until the required size of $[eqn]$  is reached, the cells marked as final being a subset of the cells of $[eqn]$ . It will be shown that at the end of Algorithm 1, $[eqn]$  is in the set of final cells. In Algorithm 1, the consensus values are calculated in the final cells (step 12), which are the cells of the fine grid $[eqn]$ . Therefore, for the value $[eqn]$  calculated in step 13, it is true that $[eqn]$ . Note also that for any cell of size $[eqn]$  that contains $[eqn]$ , according to Theorem 2, $[eqn]$   $[eqn]$ . Therefore, in the bounding phase (steps 13–14) for a cell containing $[eqn]$ , $[eqn]$ , so this cell cannot be marked passive in step (14). The same cell cannot be made passive in step 10 either, because $[eqn]$  is located inside the search area. At the beginning of the algorithm, $[eqn]$  is part of a larger active cell (step 6). Active cells are eventually (a) divided into smaller active cells (steps 8–9) or (b) smaller final cells (step 12) or (c) marked as passive (steps 10 and 14). In case (a), $[eqn]$  is preserved in one of the active cells, while in case (b), $[eqn]$  becomes one of the final cells. Case (c) does not apply to a cell containing $[eqn]$ , as was shown above. Thus, a cell containing $[eqn]$  is either active or final. Since the algorithm always divides the active cells into smaller cells, at the end of the iteration, the cell containing $[eqn]$  necessarily becomes one of the final cells. It has been shown that every cell of $[eqn]$  that gives the maximum consensus function is part of the set of final cells. Since the maximum is searched in this set (step 15), Algorithm 1 finds all of the global maxima over the grid $[eqn]$ . □

4. Performance Evaluation

The performance of the proposed consensus-function-based method (CF) is evaluated by simulation examples and real measurements. For comparison, we use the iterative LS algorithm and the three-dimensional extension of the COM-W algorithm [32], which provides a closed-form solution.

4.1. Simulations

The simulations utilize a sensor setup that was previously used in a real shooter localization experiment [38]. The bird’s eye views of the simulation setups are shown in Figure 6, using 15, 25, and 35 sensors, denoted by grey circles. The elevations of the sensors were in the range of −0.4 m to 10.9 m. The six simulated target positions are shown by red crosses and are listed in Table 1. As Table 1 shows, there were 2D and 3D experiments. In the 2D experiment, the target elevation was known to be at 0 m. This situation is common when the target, such as an autonomous vehicle, is moving on a plane. In the 3D experiments, the target elevation was also to be determined, which is a typical scenario when trying to locate a flying object. It is important to note that the sensor positions were three-dimensional in both the 2D and 3D experiments.

The simulations were carried out as follows:

• For each target position, we calculated the exact distances $[eqn]$ between the target position and the sensors.
• We calculated the exact times of arrivals $[eqn]$ as $[eqn]$ .
• We added measurement noise with normal distribution $[eqn]$ to $[eqn]$ and added additional measurement noise with normal distribution $[eqn]$ to $[eqn]$ to emulate the faulty (outlier) sensors. The number of outliers was $[eqn]$ For each target position, 100 independent measurements were created.
• Using measurements $[eqn]$ and the sensor positions $[eqn]$ , the estimated target positions were calculated by LS, COM-W, and CF.
• The tests were conducted using Matlab version R2021b on a computer with i5-8265 CPU with clock frequency of 1.6 GHz, and 24 GB of RAM.
• The LS algorithm was started from a random position within 1 m of the true position.
• The CF method was implemented in Matlab according to Algorithm 1. Apart from Matlab’s built-in vector operations, no acceleration methods (e.g., multithreading) were used.
• The final resolution of the CF method was $[eqn]$ .

In the first experiment, there were no outliers present and $[eqn]$ was $[eqn]$ (equivalent to a distant measurement error of $[eqn]$ ). The estimated positions for the 3D case are shown in Figure 6 with markers ‘x’, close to the true positions. Table 2 and Table 3 summarize the results for the 2D and 3D experiment, respectively: for each target position, the square root of the Cramer–Rao lower bound (CRLB) is shown, followed by the root-mean-square error (RMSE) of the tested algorithms. In this test, all three methods performed very well with an RMSE close to the theoretical optimum ( $[eqn]$ ). According to the results, there is no practical difference in accuracy among the three methods in this outlier-free experiment.

In the second experiment, we created outlier measurements as well. We selected the number $[eqn]$ of outliers randomly among $[eqn]$ , independently in each experiment. In each experiment, we generated $[eqn]$ correct measurements with a noise level of $[eqn]$ and $[eqn]$ outliers with a noise level of $[eqn]$ . The results for the 3D case are shown in Figure 7. Comparing Figure 6 and Figure 7, it is apparent that the LS and COM-W methods have significantly increased variance, with errors occasionally reaching tens of meters. The CF method, however, has a small estimation error. The detailed results for the 2D and 3D cases are listed in Table 4 and Table 5, respectively. Clearly, neither the LS nor the COM-W tolerated well the presence of outliers; both algorithms had several meters of errors (the LS solver also diverged several times, but these cases are omitted from the analysis). In 2D, the mean error of LS and COMW methods increased to 2–3 m, while in 3D, the error of these methods was between 3 and 8 m. The CF method, however, provided low error levels, in 2D, around 0.2 m and in 3D, between 0.3 m and 0.5 m, which was close to the theoretical optimum.

The distribution of the error is somewhat visible in Figure 7, but for a more comprehensive visualization, Figure 8a presents the cumulative distribution function (cdf) of the localization error for the LS, COM-W, and CF methods. To create the figure, we conducted 1000 simulation experiments for position #1 in 3D, with $[eqn]$ . The vertical dashed lines show the mean error values, corresponding to the values shown in Table 5.

The comparison of the results of the CF method in the first and second experiments shows that the level of RMSE increases when outliers are present. Note that the theoretical limit (CRLB) corresponds to the case where every sensor provides good measurements, which is not the case here. The reduction in the number of functional sensors led to an increase in the error level.

We measured the mean execution times of the algorithms as a function of the number of sensors and they are presented in Table 6. In the experiments, the LS method was the fastest, with an execution time of approximately 4 ms in the 2D case and 5 ms in the 3D case. The CF method required approximately 20 ms in 2D and 90 ms in 3D. The slowest method was COM-W with an average execution time of more than 200 ms in 2D and 1.3 s in 3D. It is also apparent that the execution time of the COM-W method strongly depends on $[eqn]$ , whereas the LS and CF methods did not show significant dependence on it. As an example, the cumulative distribution functions of the execution times for the 3D experiment at position #1 are shown in Figure 8b. The vertical dashed lines show the mean execution times for this particular experiment. The sharp transition of COM-W indicates that the execution time of COM-W is quite deterministic, with an average of approximately 950 ms. The iterative LS has a run-time that varies between 3 and 20 ms, with a mean value of 6 ms. The execution time of CF is also dependent on the specific scenario, with the majority of values falling between 10 and 400 ms, with a mean value of 165 ms.

The experiment in Figure 9 illustrates the fault tolerance of CF. In the test, we used target position #5 with 35 sensors, and varied the number of outliers from 0 to 25, with the outlier sensors chosen randomly. For each outlier number, we conducted 100 independent experiments, and measured the RMSE for the methods CF, LS, and COM-W. The figure shows the mean estimation errors with solid lines, and the upper and lower edges of the shaded areas correspond to the 90th percentile and the 50th percentile (median) of the error, respectively.

As Figure 9 shows, the error levels for all three methods are close to the Cramer–Rao lower bound, when there are no outliers present. When the first outliers appear, the error of LS and COM-W increases to several meters and the error continues to increase for a higher number of outliers. It is noteworthy that the COM-W method tolerates outliers better than the LS method. The error of the CF method increased slightly with the presence of outliers but remained close to the theoretical value until the number of bad measurements exceeded 25 (i.e., 71% of the total number of 35). Above this outlier number, the performance of the CF method decreased significantly. The ability of the CF method to tolerate more than 70% of sensors being outliers clearly demonstrates its robustness.

4.2. Measurements

The performance of the CF method was evaluated using two real-world measurements. The first measurement was obtained from the public database UTIL [54], where measurements were performed using ultrawideband (UWB) radios. The data were collected in an indoor flight arena of size 7 m × 8 m × 3.5 m. The setup comprised 8 fixed UWB anchors, and the target was a UWB unit mounted on a flying quadrotor platform. The reference trajectory of the target was measured by a millimeter-accuracy optical system, using 10 Vicon Vantage+ cameras. We used measurement record const4-trial6-tdoa2-traj3, which contained obstacles (wooden and metal boxes) to create challenging NLOS measurements [54]. The measurement setup is shown in Figure 10, where the grey circles indicate the beacons and the red line represents the target trajectory.

Note that in the tests, we estimated the target position in every measurement position separately, and no filters (e.g., Kalman filters) were utilized to enhance the estimation quality. The estimated target positions are shown as colored dots in Figure 11, for the LS, COM-W, and CF methods. According to the results, the vertical inaccuracy is significantly higher than the inaccuracy in either the x or y direction. This is due to the particular setup of the beacons. The horizontal-only (RMSE-xy) and three-dimensional (RMSE-xyz) values are shown in Table 7. The 2D errors are around 0.2 m, while the 3D error is around 0.5 m. In this experiment, the LS method had the highest RMSE of 0.62 m, while the COM-W handled the NLOS situations better with an RMSE of 0.47 m. The CF provided the most accurate estimates with an error of 0.41 m. Note that the accuracy of the reference algorithm in [54] is 0.45 m, which was also outperformed by CF.

The execution times are also shown in Table 7. In this setup, both LS and COM-W required approximately 4 ms to provide an estimate, while the execution time of CF was around 11 ms.

The data of the second measurement are from a counter-sniper experiment [38]. The experiment was conducted in a village that had been utilized for military training. The village comprised several streets and numerous buildings. The 57 acoustic sensors were deployed on the street level and the window sills. The shooter fired a sniper weapon from various reference positions and the muzzle blast of the weapon was detected by the sensors. The ten reference shooter positions are listed in Table 8. The sensors are represented by grey circles in the bird’s eye view of Figure 12, where the elevations were between −0.4 m and 10.9 m. Note that the setup is a challenging NLOS situation: in each experiment, multiple sensors lacked a direct line of sight. These sensors were either unable to detect the muzzle blast, or they detected signals that were reflected by one or more of the walls of the buildings. The reflected signals arrived significantly later than the signals that should have arrived in a line-of-sight scenario (see Figure 1), so these measurements must be considered outliers.

In Figure 12, the estimated positions are shown by magenta, cyan, and blue x markers for the LS, COM-W, and CF, respectively. The enlarged image illustrates that the CF method produced estimations that closely match the true positions, whereas the LS and COM-W methods had significantly larger errors. The corresponding estimation errors are listed in Table 9, along with the number $[eqn]$ of sensors providing measurement results and the consensus function value $[eqn]$ . In every case, there is a difference between $[eqn]$ and $[eqn]$ , which means measurements with large errors are present. The significant errors are due to the lack of line of sight for a number of the sensors. The outliers cause high estimation errors for LS and COM-W. The mean error was 6.6 m for COM-W and 6.8 m for LS. The maximum error was 12.5 m for COM-W and 14.2 m for LS. The CF method tolerated outliers much better: the mean error was below 0.9 m, while the maximum error was 1.9 m.

The execution times of the algorithms are also shown in Table 9. The LS method was the fastest, with a mean run-time of 18 ms. The CF required 153 ms on average, while for COM-W, the mean execution time was 949 ms. Note that the execution time of COM-W depends on the number N of sensors, while for LS and CF methods, the speed of convergence depends on the actual shape of the error surface, which is hard to predict.

In the case of the indoor measurement, the CF method demonstrated comparable speed with the LS and COM-W methods. In the context of the shooter problem, however, the CF method was significantly slower than the LS method, yet proved to be considerably faster than the COM-W method. The indoor measurement is a typical small-scale problem, whereas the shooter application can be considered a large-scale problem. According to experiments, the CF scales with the number of sensors reasonably well, and the method can be applied for large-scale problems.

The results of the run-time measurements indicate that the execution time of the LS method is only slightly affected by the actual measurement setup, and that it is not directly dependent on the number of sensors. The execution time of the COM-W is influenced by the number of sensors only. The run-time of the CF is in part dependent on the number of sensors and the actual shape of the consensus function. Consequently, the speed of the CF method is the least predictable of the three methods, especially in dynamically changing environments.

The shooter test database allows for a comparison of the proposed CF and previous consensus-based methods. In [38], the Generalized Bisection method was applied to accelerate the computation, while in [53], a mixed RANSAC–consensus approach was utilized. The brute-force computation of the consensus function in the $[eqn]$ search space with a grid size of $[eqn]$ = 0.1 m would require a $[eqn]$ evaluation of the consensus function. The Generalized Bisection method was reported to use $[eqn]$ consensus function calls [38], while the RANSAC–consensus solution required $[eqn]$ function calls to solve the shooter localization problem. The proposed CF method required an average of $[eqn]$ consensus function calls in the shooter test database, being the most efficient from among the consensus-based approaches. Note that the CF method guarantees that the global maximum is found, while the global convergence of other methods is not guaranteed.

5. Conclusions

A fast and reliable calculation method has been proposed to find the maximum of the consensus function on a finite grid. The algorithm has been theoretically proven to find the global maximum of the consensus function on the grid. The proposed CF method is capable of solving the hyperbolic localization problem even in the presence of a large number of outlier measurements (e.g., from NLOS measurements). The method is extremely robust: it provides estimates with an error level close to the theoretical optimum (Cramer–Rao lower bound) when a large percentage (even as much as 70%) of the measurements are outliers.

The application of the CF method is beneficial in situations where measurements include outliers (e.g., NLOS measurements or bad sensor detections), even a large number of them, and the necessary level of accuracy is still near the theoretical minimum. One disadvantage of the CF method is that its computational cost is higher than that of the simple LS method. Also, its speed depends on the actual error surface; thus, it cannot be predicted and guaranteed. Therefore, the application of the CF method is suggested in soft real-time systems.

The consensus-based approach is not limited to TDOA measurements only. The derivation of hybrid consensus functions including TDOA and angle of arrival (AOA) or angle difference of arrival (ADOA) measurements requires further research. Another topic for future work is the more efficient implementation of the proposed method through the use of parallel implementation.