Open AccessArticle

Passive TDOA Emitter Localization Using Fast Hyperbolic Hough Transform

Gyula Simon

^1,*

and

Ferenc Leitold

Alba Regia Technical Faculty, Óbuda University, 8000 Székesfehérvár, Hungary

Institute of Information Technology, University of Dunaújváros, 2400 Dunaújváros, Hungary

Author to whom correspondence should be addressed.

Appl. Sci. 2023, 13(24), 13301; https://doi.org/10.3390/app132413301

Submission received: 12 September 2023 / Revised: 24 November 2023 / Accepted: 4 December 2023 / Published: 16 December 2023

(This article belongs to the Special Issue Machine Perception and Learning)

Download

Browse Figures

Versions Notes

Abstract

A fast Hough transform (HT)-based hyperbolic emitter localization system is proposed to process time difference of arrival (TDOA) measurements. The position-fixing problem is provided for cases where the source is known to be on a given plane (i.e., the elevation of the source is known), while the sensors can be deployed anywhere in the three-dimensional space. The proposed solution provides fast evaluation and guarantees the determination of the global optimum. Another favorable property of the proposed solution is that it is robust against faulty sensor measurements (outliers). A fast evaluation method involving the hyperbolic Hough transform is proposed, and the global convergence property of the algorithm is proven. The performance of the algorithm is compared to that of the least-squares solution, other HT-based solutions, and the theoretical limit (the Cramér–Rao lower bound), using simulations and real measurement examples.

Keywords:

source localization; emitter localization; Hough transform; robustness

1. Introduction

Localization and positioning in GNSS-denied areas are primary services for many emerging applications, including indoor pedestrian tracking [1], robotics [2], UAV tracking [3], healthcare [4], engineering [5], wildlife protection [6], and security [7], just to name a few. In particular, passive time difference of arrival (TDOA) emitter localization methods are widely used in many application fields, e.g., geolocation [8,9,10], mobile user tracking [11], speaker localization [12], damage localization [13], or acoustic shooter localization [7,14,15,16,17]. In one-stage or direct position determination (DPD) methods, the received signals are used to determine the position estimate, without an explicit calculation of the TDOA values [18,19]. In two-stage methods, the TDOA values are determined from the measured signals in stage 1 (e.g., using various cross-correlation techniques [20,21,22] or taking the differences of the direct time of arrival measurements [14,16]) and then the location estimation (also called position fixing) is performed in stage 2.

DPD methods are reported to have higher accuracy when the signal-to-noise ratio of the measurements is low, while their disadvantages include complex computations and large communication bandwidth required to transfer the raw signal from the sensors to the processing unit. Although not optimal, two-stage methods are asymptotically equivalent to DPD methods [23] and are popular in massively distributed sensor systems, because only a small amount of data (basically only time stamps) is generated at the sensors in stage 1, which requires a small communication bandwidth. Additionally, the localization performed in stage 2 may also be performed in a relatively simple manner [14].

Note that for applications in which the time measurement of incoming events is achieved with the sensor units themselves, only the two-stage framework can be utilized (e.g., shooter localization [14,15,16]). In this paper, the focus is on the position fixing of two-stage methods.

For position fixing (stage 2), the proposed solutions include iterative techniques [24,25], closed-form solutions [26,27,28,29,30], particle filters [31,32], consistency functions [14,33], and Hough transform (HT)-based solutions [31,34]. When outlier measurements are present, the location estimate will be biased, unless the outliers are removed [35,36]. The voting approach of HT, however, has the advantage of implicitly filtering out outliers; this makes the location estimates robust against bad measurements resulting from sensor faults or non-line-of-sight measurements.

In this paper, a novel HT-based TDOA emitter localization method is proposed. The location estimate is found at the maximum of the hyperbolic HT (HHT), specially tailored for hyperbolic localization. The proposed method involves a novel hierarchical search to find the maximum of the HHT on a predefined grid. The novelties of the proposed search methods are the following:

The speed of the search is higher than those of the earlier HHT solutions [31];
The search algorithm is guaranteed to find the global maximum on the predefined grid.

The outline of the paper is as follows: In Section 2, related work is reviewed. In Section 3, the proposed fast hyperbolic Hough transform is introduced. In Section 4, the performance of the proposed method is compared to that of earlier solutions through simulation examples and real measurements. Section 5 concludes the paper.

2. Related Work

The original Hough (line) transform was utilized to detect straight lines on images [37,38]. Since then, several variants and generalizations of the original HT have been proposed and applied in various fields [39]. The different versions of HT have also been successfully applied in various localization schemes. The bearing trace HT was used to process sonar data for ship tracking and monitoring [40]. In the acoustic shot localization system [41], the HT of the time-stamped angle of arrival measurements was used to detect projectile trajectories. HT was also utilized to detect and localize reflector surfaces to prevent suboptimal solutions in indoor localization [42]. An HT-based solution was proposed for hyperbolic localization in [31], and the model was also extended to contain various types of measurements [34]. We will refer to the variant of HT adapted for processing TDOA measurements as hyperbolic HT (HHT).

Various solutions were proposed to decrease the processing time of HT: probabilistic HT [43] and randomized HT [44] reduce the computational cost by randomly decreasing the number of processed pixels, while fast HT [45] uses a hierarchical approach, which increases resolution only where detection is possible.

To increase the processing speed of HHT, in [31], two approaches were proposed. The randomized HHT (R-HHT) uses randomly selected points of the parameter space. This approach, however, still requires a high number of samples for reasonable accuracy (in [31] as much as 50% and 75% of the grid points were utilized); thus, the increase in speed is modest. The hybrid HHT (H-HHT) is performed in two stages: in stage 1, an R-HHT is executed. In stage 2, a grid search is performed around the maximum value found in step 1, using a small grid size. According to the experiments, this approach provides faster evaluation, but the speed gain is still moderate. It is worth noting that neither R-HHT nor H-HHT guarantees that the global optimum is indeed found.

In this paper, a fast HHT (F-HHT) is introduced that significantly decreases the required amount of computation and guarantees convergence to the global maximum of the parameter space on a predefined grid.

3. Fast Hyperbolic Hough Transform

3.1. Problem Formulation

The measurement scenario is shown in Figure 1a. Let the sensors be denoted by

S_{i}, i = 1,2, \dots, N

;

S_{i}

is deployed at location

(x_{i}, y_{i}, z_{i})

. The event, emitted by emitter

E

at an unknown time instant

t_{0}

from an unknown position

(x_{0}, y_{0}, 0)

on a given plane, is detected by sensor

S_{i}

at time instant

t_{i}

. Notice that this situation (often referred to as the 2.5D scenario) is commonly encountered in practice when the emitter is located at ground level, but the sensors are placed at positions of arbitrary height.

The signal propagation speed (e.g., the speed of light or the speed of sound, depending on the signal type and medium) is

c

. The distance between

S_{i}

and

E

d_{i}

d_{i} = c (t_{i} + n_{i} - t_{0}),

(1)

where

n_{i}

is the time measurement error of

S_{i}

and

t_{i} + n_{i}

is the ideal time measurement.

The distance difference

d_{i, j}

is defined as follows:

d_{i, j} = d_{i} - d_{j} = c (t_{i} - t_{j} + n_{i} - n_{j}) = c t_{i, j} + Δ d_{i, j},

(2)

where

t_{i, j} = t_{i} - t_{j}

is the measured TDOA between sensors

S_{i}

and

S_{j}

, and

Δ d_{i, j}

is the distance difference error:

Δ d_{i, j} = c (n_{i} - n_{j}) .

(3)

Let us assume that the time measurement error

n_{i}

can be modeled by Gaussian noise

n_{i} = N (0, σ_{t, i})

, and

n_{i}, n_{j}, i \neq j

are independent. Then,

Δ d_{i, j} = N (0, σ_{i, j}),

(4)

where the variance

σ_{i, j}^{2}

of the distance difference

d_{i, j}

σ_{i, j}^{2} = c^{2} (σ_{t, i}^{2} + σ_{t, j}^{2}) .

(5)

Let the reference sensor be

S_{1} .

Let us calculate the distance differences from the TDOA measurements

t_{i, 1}

{\tilde{d}}_{i, 1} = c t_{i, 1} .

(6)

The objective is to provide position estimate

({\tilde{x}}_{0}, {\tilde{y}}_{0}),

given sensor positions

(x_{i}, y_{i}, z_{i})

, measurements

{\tilde{d}}_{i, 1}

, and variances

σ_{i, 1}^{2} .

The estimation process is illustrated in Figure 1b.

3.2. Hyperbolic Hough Transform

To provide location estimates from TDOA measurements, an HHT was proposed in [31]. Now a simpler, and physically more interpretable, version is introduced.

Let us define a grid where the emitter location is searched for. Using (4), the likelihood of the emitter being at grid position

(x, y)

, given a measurement

{\tilde{d}}_{i, 1}

, is the following:

L (x, y| {\tilde{d}}_{i, 1}) = \frac{1}{\sqrt{2 π} σ_{i, 1}} {e^{- \frac{1}{2} (\frac{{\tilde{d}}_{i, 1} - d_{i, 1} (x, y)}{σ_{i, 1}})}}^{2},

(7)

where

{\tilde{d}}_{i, 1}

is the measured distance difference between sensors

S_{i}

and

S_{1}

, according to (6), and

d_{i, 1} (x, y)

is the exact distance difference at point

(x, y)

, as follows:

d_{i, 1} (x, y) = d_{i} (x, y) - d_{1} (x, y)

(8)

and

d_{i} (x, y) = \sqrt{{(x - x_{i})}^{2} + {(y - y_{i})}^{2} + z_{i}^{2}} .

(9)

The HHT uses a voting system, where a pair of sensors casts a vote to each of the grid points: a high/low vote represents a high/low likelihood of the emitter being at the given grid point. In the voting process, the normalized version

L_{n}

of likelihood function (7) is utilized, where the maximum value of each vote is

\frac{1}{N - 1}

L_{n} (x, y | {\tilde{d}}_{i, 1}) = \frac{1}{N - 1} {e^{- \frac{1}{2} (\frac{{\tilde{d}}_{i, 1} - d_{i, 1} (x, y)}{σ_{i, 1}})}}^{2} .

(10)

Thus, the HHT value at grid position

(x, y)

is the following:

A (x, y) = \sum_{i = 2}^{N} L_{n} (x, y | {\tilde{d}}_{i, 1}) .

(11)

The highest attainable value of

A (x, y)

1

; this occurs when all the sensor pairs vote with their maximum possible value. If the measurement error is zero, i.e.,

n_{i} = 0, i = 1, 2, \dots, N

, then at the true source position

A (x_{0}, y_{0}) = 1

. In practical scenarios, where measurement noise is present,

A (x_{0}, y_{0})

is lower than, but close to 1.

The position estimate is at the maximum of

A (x, y) :

({\tilde{x}}_{0}, {\tilde{y}}_{0}) = \underset{(x, y)}{argmax} A (x, y) .

(12)

The HHT is illustrated in Figure 2, where four sensors are used, with

S_{1}

serving as the reference. The sensors in the example are placed on the plane of the source. The sensor positions are the same, but the source is placed in different positions in Figure 2a,b. In this scenario, 3 pairs of sensors vote: darker points indicate higher votes (white corresponds to zero and black corresponds to one). The result of votes of sensor pair

(S_{i}, S_{1})

is denoted by

H_{i, 1}

, forming a hyperbola in Figure 2. The image corresponds well to the well-known two-dimensional geometric interpretation:

L_{n}

is maximal when

{\tilde{d}}_{i, 1} = d_{i, 1} (x, y)

, where

d_{i, 1} (x, y)

is the difference of distances

d_{i}

and

d_{1}

, as shown in (8). Thus, the set of points where

L_{n}

is maximal forms a hyperbola with focal points at

S_{i}

and

S_{1}

and a major axis of

{\tilde{d}}_{i, 1}

. The sum of all votes is maximal where the hyperbolas intersect, providing the location estimate (the darkest point).

Note that the width of the hyperbola depends on both the accuracy of measurements, i.e.,

σ_{i, 1}

, and the geometric dilution of precision (GDOP) [46]. The higher

σ_{i, 1}

, the wider the likelihood function, and thus the hyperbola, as follows from (10). High GDOP makes the hyperbola even wider, particularly noticeable for

H_{2,1}

in Figure 2b:

H_{2,1}

is narrow in proximity to

S_{2}

and becomes wider at the top of the figure. The GDOP naturally effects the accumulated HHT and, consequently, the position estimate: in Figure 2a, the small dark area around the location estimate indicates good GDOP, resulting in high confidence in the location estimate. In Figure 2b, however, the dark area is significantly larger, forming an ellipsoid-like pattern. As a result, small measurement errors can lead to significant differences in the estimate, resulting in lower confidence in the location estimate.

Notes:

In practice, the exact values of the measurement variances are usually not known. In such cases, $σ_{i, 1}^{2}$ can serve as an estimate of the true value. If the estimated value is smaller than the true variance then the source position may be outside of the skirt of the generated “wide hyperbola” $H_{i, 1}$ . On the other side, if the estimated variance is higher than the true value then the source position is safely included inside $H_{i, 1}$ . Therefore, in practice, the estimate $σ_{i, 1}^{2}$ should be an upper estimate of the true variance.
In the 2.5D case, the shape of $H_{i, 1}$ is not necessarily a hyperbola, but rather a general intersection of a 3-dimensional hyperboloid surface and the plane of the source.

Similarly to other variants of Hough transform, the HHT can be evaluated on a grid. The size of the grid has important implications:

(a): If the grid size is large, the maximum may be missed;
(b): The grid size determines possible accuracy: a smaller grid size can provide higher accuracy;
(c): Smaller grid size implies a higher computational cost. Note that if the grid contains $M \times M$ points, then the number of likelihood function evaluations will be $M^{2}$ , and the search for the maximum will also require $M^{2}$ (albeit simpler) operations.

To improve calculation speed, in [31] the randomized R-HHT was proposed. However, this method faces an obstacle resulting from implication (a): the number of required sampling points to achieve a high probability of finding the global maximum remains high; the improvement, reported in [31], was not better than twofold. The hybrid H-HHT method [31] increases the accuracy and reduces the number of trials, using a fine-grid search in the second phase. In the next subsection, the fast HHT will be introduced, which not only significantly reduces the computational costs but also guarantees that the optimum on the given grid is found.

3.3. Fast HHT

The basic idea behind the F-HHT is similar to that of the fast HT [45]: instead of calculating the transform on a dense grid, a hierarchical approach is utilized. First, a search with a large grid size is performed. In the next round, the search is focused on “promising” areas, now with smaller grid size. The process is iteratively repeated until the required accuracy (i.e., grid size) is achieved. The main novelty is the utilization of two metrics: the exact value of the HHT at a given point and the upper bound of the HHT in the vicinity of this point. This approach allows fast and provable convergence to the global maximum of the HHT.

The search space is divided into rectangular tiles with grid points at the centers of the tiles, as shown in Figure 3. In the first round, the tile size is

D_{1}

; in the second round, “promising” tiles are divided into four smaller tiles of size

D_{2} = D_{1} / 2

. In general, in the n-th round, the tile size is reduced to

D_{n} = \frac{D_{1}}{2^{n - 1}} .

(13)

To ensure that the global optimum is found, the F-HHT utilizes two metrics at each grid point

(x, y)

. The first is the exact value

A (x, y)

, according to (11), while the second is the upper bound

\bar{A} (x, y)

of the HHT in the actual tile surrounding

(x, y)

Consider a tile of size

D_{n}

with center

P = (x, y)

, as shown in Figure 4. The difference of distances

P S_{i} = d_{i}

and

P S_{1} = d_{1}

d_{i, 1}

, according to (8). Let

P^{'}

be an arbitrary point inside the tile, where

P^{'} S_{1} = d_{1}^{'}

P^{'} S_{i} = d_{i}^{'}

, and

d_{i, 1}^{'} = d_{i}^{'} - d_{1}^{'}

\begin{matrix} d_{i, 1}^{'} - d_{i, 1} = \\ (d_{i}^{'} - d_{1}^{'}) - (d_{i} - d_{1}) = (d_{i}^{'} - d_{i}) - (d_{1}^{'} - d_{1}) \end{matrix}

(14)

From the triangle inequality applied to

S_{i} P P^{'}

, it follows that

- P P^{'} \leq d_{i}^{'} - d_{i} \leq P P^{'} .

Since

P P^{'} \leq \frac{D_{n}}{\sqrt{2}}

, the following inequity holds, for all

i

- \frac{D_{n}}{\sqrt{2}} \leq d_{i}^{'} - d_{i} \leq \frac{D_{n}}{\sqrt{2}} .

(15)

From (14) and (15), the following inequality follows:

- \sqrt{2} D_{n} \leq d_{i, 1}^{'} - d_{i, 1} \leq \sqrt{2} D_{n} .

(16)

The likelihood function at any point

P^{'}

in the tile is

L_{n} (P^{'} | {\tilde{d}}_{i, 1}) = \frac{1}{N - 1} {e^{- \frac{1}{2} (\frac{{\tilde{d}}_{i, 1} - d_{i, 1}^{'}}{σ_{i, 1}})}}^{2} .

(17)

Since

|{\tilde{d}}_{i, 1} - d_{i, 1}^{'}| = |({\tilde{d}}_{i, 1} - d_{i, 1}) - (d_{i, 1}^{'} - d_{i, 1})|,

(18)

using the reverse triangle theorem, it follows that

|{\tilde{d}}_{i, 1} - d_{i, 1}^{'}| \geq ||{\tilde{d}}_{i, 1} - d_{i, 1}| - |d_{i, 1}^{'} - d_{i, 1}|| .

(19)

Let us denote

Δ d_{i, 1} = |{\tilde{d}}_{i, 1} - d_{i, 1}| .

(20)

Δ d_{i, 1} > \sqrt{2} D_{n}

then, using (16), the inequity (19) can be rewritten as

|{\tilde{d}}_{i, 1} - d_{i, 1}^{'}| \geq Δ d_{i, 1} - \sqrt{2} D_{n},

(21)

otherwise it is certainly true that

|{\tilde{d}}_{i, 1} - d_{i, 1}^{'}| \geq 0 .

(22)

Notice that

L_{n} (P^{'} | {\tilde{d}}_{i, 1})

in (17) increases if

|{\tilde{d}}_{i, 1} - d_{i, 1}^{'}|

decreases. Thus, using (21) and (22), in the tile with center

P

, an upper bound

{\bar{L}}_{n} (P | {\tilde{d}}_{i, 1})

for

L_{n} (P^{'} | {\tilde{d}}_{i, 1})

can be obtained, as follows:

{\bar{L}}_{n} (P | {\tilde{d}}_{i, 1}) = \{\begin{matrix} \frac{1}{N - 1} {e^{- \frac{1}{2} (\frac{Δ d_{i, 1} - \sqrt{2} D_{n}}{σ_{i, 1}})}}^{2} & i f Δ d_{i, 1} > \sqrt{2} D_{n} \\ \frac{1}{N - 1} & o t h e r w i s e, \end{matrix}

(23)

while the exact value of the likelihood function in center

P

L_{n} (P | {\tilde{d}}_{i, 1}) = \frac{1}{N - 1} {e^{- \frac{1}{2} (\frac{Δ d_{i, 1}}{σ_{i, 1}})}}^{2} .

(24)

Thus, the HHT in a grid point

P

is the following:

A (P) = \sum_{i = 2}^{N} L_{n} (P | {\tilde{d}}_{i, 1}),

(25)

and the upper bound of the HHT in this tile is

\bar{A} (P) = \sum_{i = 2}^{N} {\bar{L}}_{n} (P | {\tilde{d}}_{i, 1}) .

(26)

Using the above results, Algorithm 1 describes the operation of the F-HHT.

Algorithm 1 Operation of the F-HHT

Input:

measurements ${\tilde{d}}_{i, 1}, σ_{i, 1}, i = 2,3, \dots, N$
sensor positions $(x_{i}, y_{i}, z_{i}), i = 1,2, \dots, N$
desired grid size $D_{f i n a l}$

Initialization:

4.: create a grid over the search area $S$ with grid size $D_{1}$
5.: add each tile to the PromisingList
6.: ${G r i d S i z e : = D}_{1}$
7.: set the minimal required value $A_{m a x m i n}$ $(e . g ., t o \frac{2}{N - 1})$

Iteration:

8.: while $G r i d S i z e > D_{f i n a l} / 2$
9.: calculate $A (P)$ and $\bar{A} (P)$ for each tile in PromisingList
10.: $bestA : = \max (A (P))$
11.: $b e s t P : = \underset{P}{argmax} (A (P))$
12.: pruning 1: remove tiles with $\bar{A} (P) < A_{m a x m i n}$
13.: pruning 2: remove tiles with $\bar{A} (P) < b e s t A$
14.: $G r i d S i z e : = G r i d S i z e / 2$
15.: replace each tile with four smaller tiles of size GridSize
16.: end

Output: bestP

The inputs of F-HHT are the measurements, sensor positions, and the size of the grid on which the maximum is searched for (see steps 1–3). In step 4 of the initialization phase, the initial coarse grid is created. Then, in step 5, each tile is added to the list of promising tiles (PromisingList). The initial coarse grid size is set in step 6. The value

A_{m a x m i n}

in step 7 is the minimum required value of

\max A (P)

, i.e.,

\max (A (P)) > A_{m a x m i n}

. In a conservative solution

A_{m a x m i n} = \frac{2}{N - 1}

, since at least two hyperbolas intersect at the source position. When multiple sensors are available and there is a priori knowledge of the minimum number

R

of reliable measurements, then

A_{m a x m i n}

may be set close to, but below,

\frac{R}{N - 1}

The iteration is performed between steps 8 and 16. In step 9, the F-HHT algorithm calculates the exact HHT (25) and the upper bound (26) for each promising tile center. The current highest exact value is

b e s t A

(step 10) at positions

b e s t P

(step 11). A tile with center

P

is non-promising (the maximum cannot be inside the tile) if

\bar{A} (P) < A_{m a x m i n}

\bar{A} (P) < b e s t A

. In the pruning steps 12 and 13, the non-promising tiles are removed from the further search. In step 15, promising tiles are replaced by four new promising tiles with smaller size (see Figure 3). When the required grid size is reached, the center(s) of tile(s) with the highest HHT (

b e s t P

) form the output.

Note that the output

b e s t P

may contain multiple grid points (tile centers). In such cases, the location estimate can be calculated as the mean of bestP.

3.4. Global Convergence

In this subsection, Theorem 1 will be proven:

Theorem 1.

Let the global maximum of the HHT, calculated over search area

S

on fine grid

G_{f}

of grid size

D_{G_{f}}

, be

A_{G_{f}}

, such that the maximum is at grid point(s)

P_{0}

, i.e.,

A_{G_{f}} = A (P_{0}) = \max_{P \in G_{f}} A (P)

. Let us create a coarse grid

G_{c}

by omitting grid lines from

G_{f}

, such that the grid size of

G_{c}

D_{1} = 2^{k - 1} D_{G}

k \in N

. If the output of the F-HHT algorithm, executed over

S

with initial grid

G_{c}

, using

k

iterations, is

P_{0, F H H T}

, then

P_{0, F H H T} = P_{0}

The theorem ensures that the global maximum of the HHT, calculated on a fine grid with size

D_{G_{f}}

, can be found with the F-HHT in

k

iteration steps, starting from a much coarser grid of size

2^{k - 1} D_{G_{f}}

. Note that

P_{0}

may contain one or more grid points of

G_{f}

Lemma 1.

If there is a promising tile

T

with center

P

that contains one or more points of

P_{0}

, then this tile is preserved (not removed) during the pruning step.

Proof of Lemma 1.

For tile

T,

\bar{A} (P) \geq A (P_{0}) = A_{G_{f}}

. Since

A_{G_{f}} > A_{m a x m i n}

\bar{A} (P) > A_{m a x m i n}

, therefore

T

is not removed in step 12. Notice that each coarser grid during the iteration is part of the fine grid

G_{f}

. Thus, for the actual

b e s t A

calculated in step 10, it holds that

b e s t A \leq A_{G_{f}}

. This ensures that

T

is not removed in step 13, either. Thus

T

is preserved. □

Proof of Theorem 1.

At the start of the first iteration step, all tiles are marked promising, thus the tiles containing one or more points of

P_{0}

are also marked as promising. According to Lemma 1, these tiles are preserved in the pruning step and then they are replaced by four smaller promising tiles. These smaller tiles also contain the points of

P_{0}

. Thus, in the beginning of the second iteration step, promising tiles still contain all points of

P_{0}

. This is true for all iteration steps, thus at the beginning of the kth (last) iteration step, promising tiles contain all points of

P_{0}

Notice that in the last iteration step the set of search lists PromisingList is a (small) subset of the fine grid

G_{f}

(with grid size of

D_{1} / 2^{k - 1} = D_{G_{f}}

), where the HHT was evaluated. Thus, points of

P_{0}

are the center points of some promising tiles. Since all points in

P_{0}

are preserved, the maximum points

b e s t P

, computed in the last iteration step, are necessarily the same as the maximum points calculated by the HHT on the whole grid

G_{f}

. Therefore

P_{0, F H H T} = P_{0}

. □

4. Performance Evaluation

In this section, the performance of the F-HHT will be studied and comparisons with other methods and the theoretical limit will be provided. For comparison purposes, the standard and widely utilized least squares (LS) method was chosen, along with two earlier HHT methods: random HHT (R-HHT) and hybrid HHT (H-HHT). The performances of the localization methods will be investigated through simulations. First, the test setup will be introduced, followed by demonstrative examples to illustrate the operation of the F-HHT. The performance comparisons include accuracy, speed, computational cost, and also robustness. Finally, the performance of the proposed method will be illustrated with real measurements.

4.1. Test Simulation Setup

The simulation test setup, using 7 sensors, is shown in Figure 5, and the sensor locations are listed in Table 1.

In the tests, two source positions were used: a near-range position

P_{s} = (157.3 m, 113.9 m, 0 m)

and a long-range position

P_{l} = (250.2 m, 1280.4 m, 0 m) .

Notice that the source positions are known to be on the ground level, according to the 2.5D problem.

In the tests,

c = 340 m / s

was utilized (speed of sound), and the measurement errors were modeled by additive white noise. Two different noise levels were used:

σ_{t, i} = 210 μ s

and

σ_{t, i} = 2.1 m s,

resulting in

σ_{i, 1} = 0.1 m

and

σ_{i, 1} = 1 m

, respectively. The final resolution

D_{f i n a l}

of the F-HHT and H-HHT was

0.1 m

1 m

, and

5 m

. The parameters of the six scenarios are summarized in Table 2.

Scenarios 1–3 are near-range, while scenarios 4–6 are long-range. In both ranges, different noise levels and final grid sizes were applied, as shown in Table 2.

All of the tested methods were implemented in Matlab v. 9.11. Notes on the utilized algorithms:

The LS algorithm was Matlab’s built-in quasi-Newton solver in function fminunc. The gradient search was started from a random point within a 50m radius of the true source position.
The F-HHT algorithm used final grid size parameter $D_{f i n a l}$ . In the H-HHT, the same grid size was used in the second stage of the algorithm, where the search area was ${50 \cdot D}_{f i n a l} \times {50 \cdot D}_{f i n a l}$ around the result provided by the R-HHT in the first stage.
The search area for the HHT methods was set to [0, 0, 500, 500] in the near-range cases and [0, 0, 1500, 1500] in the long-range cases.
The R-HHT and H-HHT algorithms were implemented using the HHT as defined in section III-B.
The likelihood function (7) includes the measurement noise estimate $σ_{i, 1}$ . In practical cases, the exact value of the noise is often unknown. Therefore, it is recommended to overestimate the noise level and use the overestimated value in the HHT. To replicate this procedure in the simulations, we used $σ_{i, 1}$ to generate the measurements, and for the HHT algorithms, we provided $3 σ_{i, 1}$ as the noise level estimate.

4.2. Illustration of the Operation

In this subsection, the operation of the HHT and F-HHT will be illustrated. First, an example is shown in Figure 6, where the shape of the HHT can be observed. The HHT function values were calculated on a grid with resolution of 0.5 m. The insets show the functions around the global maximum, which is in close proximity to the true location, indicated by a blue cross. Figure 6a shows a near-range example using scenario #2, while Figure 6b shows a long-range example using scenario #5.

In the near-range example, one local minimum is visible (around location (160,20)), while in the long-range case several local minima can be observed. Thus, finding the global maximum is not a trivial problem. The impact of GDOP is also apparent in the figure: although the measurement noise is the same in both examples

(σ_{i, 1} = 1 m

), in the near-range case the peak is narrow and can be contained in a 20 m × 20 m box, while in the long-range case the peak is spread over a much larger area of approx. 100 m × 500 m.

Figure 7 illustrates the operation of the F-HHT in scenario #2. The figure shows steps 1, 2, 3, and 10 of the iteration. Red dots indicate the centers of the actual promising tiles, the figure titles also show their number

N_{t i l e}

. In the final step, the best grid point (bestP) is also shown, which is the position estimate.

In the experiment shown in Figure 7, the total number of processed tiles was approximately 1200, resulting in approximately 1200 evaluations of (25) and the same number of (26) (which has approximately the same computational complexity as (25)). This means that there were a total of around 2400 HHT evaluations. It is worth noting that the final grid size was 0.1 m, which would have required

25 \cdot 10^{6}

HHT calls of (25) in the brute force solution. Thus, in this example, the F-HHT resulted in an approximately 10,000-fold decrease in computational complexity.

Figure 8 shows the distribution of the F-HHT location estimates. Again, scenarios #2 and #5 were used; in both scenarios, 100 independent experiments were conducted; the estimates are shown by red dots. The estimates are scattered around the true source position. Despite both scenarios having the same noise level of

σ_{i, 1} = 1 m

, the long-range scenario has significantly higher variance due to the higher GDOP.

4.3. Error Analysis

The performance of the F-HHT will be compared to the theoretical limit of the Cramér–Rao lower bound (CRLB). The CRLB of our 2.5D TDOA problem, where the sensors are in 3D but the source position is searched for in a plane, can be derived using the method described in [27,47], as follows:

C R L B (x, y) = t r a c e (F^{- 1}),

(27)

where F is the Fisher information matrix

F = \frac{1}{c^{2}} G^{T} Q^{- 1} G,

(28)

with

Q = σ_{t, 1}^{2} 1_{N - 1} 1_{N - 1}^{T} + E_{N - 1} Σ_{N - 1}^{2},

(29)

Σ_{N - 1}^{2} = {[σ_{t, 2}^{2}, σ_{t, 3}^{2}, \dots, σ_{t, N}^{2}]}^{T},

(30)

G = [\begin{matrix} g_{2}^{T} - g_{1}^{T} \\ ⋮ \\ g_{N}^{T} - g_{1}^{T} \end{matrix}],

(31)

g_{i}^{T} = \frac{[x, y] - [x_{i}, y_{i}]}{{‖[x, y, 0] - [x_{i}, y_{i}, z_{i}]‖}_{2}}, 1 \leq i \leq N,

(32)

and

1_{k}

is the vector of length

k

containing 1 s and

E_{k}

is the identity matrix of size

k \times k

The results for all scenarios are shown in Table 3. The table shows the square root of the CRLB and the measured root mean square error (RMSE) for all algorithms, using 1000 independent experiments in each scenario.

In all cases, the measured errors of LS, H-HHT, and F-HHT were close to the theoretical limit. The LS algorithm was the most accurate, which is probably due to the fact that the search of the HHT was performed on a finite grid. The F-HHT provided results almost as accurate as the LS method, when the final grid size was set to a sufficiently small value.

The impact of grid size can be observed in experiments #2 and #3. The theoretical error limit was approximately 0.6m, so the grid size of 0.1m in experiment #2 was safely below it, while in experiment #3, the grid size of 1m exceeded the error limit. As expected, the estimation error in experiment #2 was closer to the theoretical limit, while in experiment #3 the error was significantly higher. The distinction between experiments #5 and #6 is less apparent, as both grid sizes (1 m and 5 m) were significantly below the theoretical error level (20 m).

Note that for the R-HHT and H-HHT algorithms, there is an additional parameter for the number of trials. This parameter defines the number of random test positions where (11) is evaluated. Increasing the number of trials decreases the error level, but the increase in accuracy is small near the theoretical level. Therefore, a small increase in accuracy requires a large increase in trial numbers. To ensure a fair performance comparison, the trial numbers of R-HHT and H-HHT were adjusted so that the error levels of R-HHT, H-HHT, and F-HHT were approximately equal (see Table 3). This allows for a comparison of the required computational power in the next section.

4.4. Evaluation Time

In this subsection, the evaluation time of F-HHT will be compared to that of other methods. The run times of the algorithms are listed in Table 4. As mentioned before, the trial numbers of the R-HHT and H-HHT were selected to achieve accuracy similar to that of F-HHT. The run times of the R-HHT and H-HHT were measured for settings corresponding to the accuracy levels, shown in Table 3. The run times were measured on a computer featuring an Intel i5 processor with a clock frequency of 1.6 GHz and 24 GB of RAM, using Matlab (without parallel processing).

The trial numbers are also shown in Table 5, along with the theoretical number of grid points

N_{g r i d}

for a full grid with a grid size of

D_{f i n a l}

In the tests, the LS solution converged in approximately the same amount of time in all scenarios. However, the evaluation times of HHT-based algorithms were spread over a wider range due to two effects: GDOP and measurement noise. The primary effect is GDOP: small GDOP results in sharp peaks in the HHT, while large GDOP widens the peaks (see Figure 6 for an illustration of this effect). The measurement noise also blurs the peaks, resulting in smaller peaks with gentler slopes.

The F-HHT algorithm performs very well when the peaks are sharp: in the near-range scenarios (#1, #2, #3) F-HHT was the fastest of all algorithms. In these cases, the elimination steps quickly excluded most of the search area, resulting in the processing of a small number of tiles and fast evaluation times.

In the long-range cases (with wider peaks), the elimination process of F-HHT was less efficient: a larger area was still active while the grid size decreased during the iteration, resulting in the processing of a large number of small tiles and thus longer evaluation times.

The evaluation time of the R-HHT is dependent on several factors, including the CRLB, the size of the search area, and the required accuracy, i.e., how close the required error is to the CRLB. Note that R-HHT takes random picks from the search area: the larger the search area, the more picks are needed to obtain a hit close to the optimum (the target area). The smaller the CRLB and the closer we try to get to the theoretical optimum, the smaller the target area, necessitating a higher number of trials. The results illustrate these effects. In scenarios #5 and #6, the CRLB and the required accuracy were the same, resulting in the same number of trials. In scenario #4, the CRLB was smaller, thus the trial number increased. In scenarios #1, #2, and #3, the effect of accuracy is also visible. The CRLB is smaller in #1 and larger in #2 and #3. As a result, the trial numbers in #2 and #3 are smaller than in the case of #1. In scenario #2, higher accuracy was required than in scenario #3 (while having the same CRLB), thus the trial number in #2 was higher. The impact of the search space can also be observed: e.g., in #4, the CRLB is smaller than in #1, but the increase in the search space compensates for this effect, resulting in both scenarios having approximately the same number of trials.

The evaluation time of H-HHT is closely related to that of R-HHT, but it greatly depends on the final grid size, too. The grid search, performed in the second stage, can be considered as an augmentation of the target area for the R-HHT: if the R-HHT in the first phase provides a position so that the grid search area around this position contains the optimal position, then the H-HHT will find a solution close to the optimum. Increasing the grid size results in a larger search area, which leads to faster evaluation. However, the accuracy decreases with larger grid sizes, thus the grid size cannot be set arbitrarily high. The effect is very apparent in scenarios #5 and #6: here, only the final grid size changed, and the larger grid size resulted in shorter run times. Both grid sizes were below the CRLB, thus both solutions provided good accuracy. In the experiments, the H-HHT demonstrated a 10–100-fold improvement over the R-HHT.

In most cases, the F-HHT significantly overperformed R-HHT and H-HHT in the experiments. Additionally, in the near-field scenarios, the F-HHT was faster than the LS, while in the far-field cases, the LS solution was more efficient.

4.5. Fault Tolerance

In this subsection, the fault tolerance of F-HHT is illustrated. For this purpose, a short-range scenario (#1) and a long-range scenario (#4) were utilized. In the experiments, various numbers of faulty sensors were used; the faulty sensors provided outlier measurements with a standard deviation 100 times larger than that of the correct sensors. The faulty sensors were selected randomly, but the reference sensor

S_{1}

was never faulty.

The measured RMSE values for the LS and F-HHT algorithms, as a function of the number

N_{o u l i e r}

of outliers, measured using 1000 independent experiments, are shown in Table 6. In the tests, the LS algorithm was found to be highly sensitive to outliers: the presence of even a single outlier significantly increased the error level. The F-HHT, however, proved to be much more robust: in the near-range scenario, the error did not significantly increase when up to two faulty sensors were present. In the far-range scenario, one faulty sensor was tolerated well (the error level increased from 2.1 m to 2.8 m), while two faulty sensors increased the error to 7.5 m. Notice that in the near-range case the LS solver behaved reasonably well in the presence of outliers, while in the far-range case its results were practically useless. The higher sensitivity to outliers in the far-range case was caused by the high GDOP.

In another experiment, a large number of sensors was utilized, as shown in Figure 9. The setup, containing

N = 57

sensors, was used in a real shooter localization system [33]. In the current simulation, the exact measurements were generated for source position (25, 20), shown by a green x in Figure 9, then white noise with

σ_{i, t} = 2 m s

(approximately equivalent to

σ_{i, 1} = 1 m

) was added to each measurement. Then, a random set of

N_{o u l i e r}

measurements was selected and a large noise with

σ_{i, t} = 200 m s

was added to them to create the outlier measurements. The measurement sets, including both correct measurements and outliers, were used as inputs for the LS and F-HHT algorithms. The number of outliers,

N_{o u l i e r},

was varied between 0 and 55, and for each value of

N_{o u l i e r}

1000 independent experiments were conducted. Figure 10 shows the theoretical error level

\sqrt{C R L B}

and the measured RMSE values for the LS and F-HHT algorithms on a logarithmic scale as a function of

N_{o u l i e r}

. Notice that the LS solver diverged in many cases, thus a constrained version (fmincon in MATLAB) was utilized during testing, which ensured that the solution remained within the search area depicted in Figure 9.

According to Figure 10, the LS solver did not tolerate any outliers, whereas the F-HHT method performed very well up to 30–40 outliers out of 57 measurements, which is remarkably robust behavior.

Note that in the experiments, the reference sensor was never faulty. If this cannot be guaranteed in practice, the role of the reference sensor must be rotated, and the solution with the highest HHT value must be selected. This process will be illustrated in the next subsection.

4.6. Real Measurements

In this subsection, real measurements will be used to validate the proposed method. The setup is identical to the one used in the simulations, as shown in Figure 9. In this real-world experiment, a gun was fired at positions listed in Table 7 and marked with red + symbols in Figure 9. Based on the sensor measurements, the location estimates were computed using the proposed F-HHT algorithm and the LS solver.

The measurements included several outliers and it was not possible to ensure the correctness of one specific reference sensor. Thus, in the tests, each sensor was selected as the reference sensor in a rotating fashion, and the F-HHT algorithm was run for each selection. The position estimate was determined based on the highest HHT value. The process is demonstrated by the measurement position (33.37, 48.24). In this experiment, 30 sensors provided measurements, and as a result, 30 independent F-HHTs were run, each using a different reference sensor. Figure 11 shows the highest values of the HHT as a function of the reference sensor index. For better visibility, the sensor indices were rearranged to display the HHT results in decreasing order: the best experiments are on the left-hand side and the worst are on the right-hand side. The corresponding positioning errors are also shown in Figure 11. Clearly, the highest HHT values correspond to small positioning errors, while small HHT indicates bad reference sensor measurement and large positioning error. In this experiment, the best HHT value obtained was 22.6, with an error of 0.43 m.

The position estimation errors are listed in Table 7 for the F-HHT and the LS solutions. The F-HHT performed well in all cases, producing errors around 0.5 m. The LS solution, due to the outlier measurements, had estimation errors in the range of 6–37 m.

Notice that in the real experiment, the exact measurement errors were not known but were estimated to be in the range of

σ_{i, 1} \approx 0.3 m

. The obtained error values correspond well with the error levels of the simulation tests.

5. Conclusions

A Hough-transform-based localization method was proposed to evaluate TDOA measurements for emitter localization purposes, where the target is located on a plane, with sensors deployed in the three-dimensional space. A new variant of the HHT was proposed along with a fast, branch, and bound-type calculation. The proposed F-HHT was proven to find the global maximum of the HHT on a predefined grid. The performance of the F-HHT was studied and compared to that of the LS solution, and earlier variants of fast HHT methods, using simulation examples.

The results showed that the error of the F-HHT is close to the theoretical limit (CRLB), similar to the LS solution and the other HHT methods. It was found that the F-HHT was faster in good GDOP situations, where the HT had a narrow peak around the global maximum. In such cases, the speed of the F-HHT was higher than that of the LS solution, and much higher than that of other HHT solutions. In cases when the GDOP was lower (i.e., the HT has wide peak around the global maximum), the F-HHT was slower than the LS solution but faster than the other HHT solutions. It is important to note that the F-HHT guarantees to find the global optimum, a property that neither the LS nor other HHT solutions can provide. The construction of further acceleration techniques, which maintain the global convergence property, is the topic of future research.

Another attractive property of F-HHT (and other HHT-based solutions) is that they are robust against outlier measurements. It was demonstrated that the F-HHT can tolerate faulty sensors (other than the reference sensor) well. In a large sensor system, F-HHT was shown to tolerate extreme situations, such as more than half of the sensors providing outlier measurements. The performance of the proposed method was also illustrated using real measurements.

Author Contributions

Conceptualization, G.S.; methodology, G.S.; software, G.S. and F.L.; validation, G.S. and F.L.; formal analysis, G.S.; writing, F.L. and G.S.; visualization, F.L.; supervision, G.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available in article.

Conflicts of Interest

The authors declare no conflict of interest.

References

Zhou, Z.; Feng, W.; Li, P.; Liu, Z.; Xu, X.; Yao, Y. A fusion method of pedestrian dead reckoning and pseudo indoor plan based on conditional random field. Measurement 2023, 207, 112417. [Google Scholar] [CrossRef]
Placed, J.A.; Strader, J.; Carrillo, H.; Atanasov, N.; Indelman, V.; Carlone, L.; Castellanos, J.A. A survey on active simultaneous localization and mapping: State of the art and new frontiers. IEEE Trans. Robot. 2023, 39, 1686–1705. [Google Scholar] [CrossRef]
Lin, H.; Zhan, J. GNSS-denied UAV indoor navigation with UWB incorporated visual inertial odometry. Measurement 2023, 206, 112256. [Google Scholar] [CrossRef]
Mercuri, M.; Sacco, G.; Hornung, R.; Zhang, P.; Visser, H.J.; Hijdra, M.; Liu, Y.H.; Pisa, S.; van Liempd, B.; Torfs, T. 2-D localization, angular separation and vital signs monitoring using a SISO FMCW radar for smart long-term health monitoring environments. IEEE Internet Things J. 2021, 8, 11065–11077. [Google Scholar] [CrossRef]
Wu, Y.; Yang, N. Study of rail damage diagnosis and localization method based on intelligent wireless acoustic sensor network. Measurement 2023, 217, 113032. [Google Scholar] [CrossRef]
Shafer, M.W.; Flikkema, P.G. Tracking small wildlife with minimal-complexity radio frequency transmitters: Near-optimal detection. IEEE Access 2023, 11, 40029–40037. [Google Scholar] [CrossRef]
Gao, S.; Liu, J.; Zong, Y.; Wang, M.; Jin, X.; Tian, G.; Dai, X. Blast source TDOA localization with time synchronization estimation based on spatial overpressure-monitoring network. Measurement 2022, 204, 112080. [Google Scholar] [CrossRef]
Li, J.; Lv, S.; Jin, Y.; Wang, C.; Liu, Y.; Liao, S. Geolocation and Tracking by TDOA Measurements Based on Space–Air–Ground Integrated Network. Remote Sens. 2023, 15, 44. [Google Scholar] [CrossRef]
Huang, S.; Pun, Y.-M.; So, A.M.-C.; Yang, K. A Provably Convergent Projected Gradient-Type Algorithm for TDOA-Based Geolocation Under the Quasi-Parabolic Ionosphere Model. IEEE Signal Process. Lett. 2020, 27, 1335–1339. [Google Scholar] [CrossRef]
Ho, K.C.; Chan, Y.T. Solution and performance analysis of geolocation by TDOA. IEEE Trans. Aerosp. Electron. Syst. 1993, 29, 1311–1322. [Google Scholar] [CrossRef]
Li, C.; Zhuang, W. Hybrid TDOA/AOA mobile user location for wideband CDMA cellular systems. IEEE Trans. Wirel. Commun. 2002, 1, 439–447. [Google Scholar] [CrossRef]
Ma, W.K.; Vo, B.N.; Singh, S.S.; Baddeley, A. Tracking an unknown time-varying number of speakers using TDOA measurements: A random finite set approach. IEEE Trans. Signal Process. 2006, 54, 3291–3304. [Google Scholar] [CrossRef]
Sun, J.; Wen, J. Target location method for pipeline pre-warning system based on HHT and time difference of arrival. Measurement 2013, 46, 2716–2725. [Google Scholar] [CrossRef]
Maroti, M.; Simon, G.; Ledeczi, A.; Sztipanovits, J. Shooter localization in urban terrain. Computer 2004, 37, 60–61. [Google Scholar] [CrossRef]
Still, L.; Oispuu, M.; Koch, W. Accuracy Study on Shooter Localization Using Incomplete Acoustic Measurements. In Proceedings of the 2020 IEEE 23rd International Conference on Information Fusion (FUSION), Rustenburg, South Africa, 6–9 July 2020; pp. 1–8. [Google Scholar] [CrossRef]
Damarla, T.; Kaplan, L.M.; Whipps, G.T. Sniper Localization Using Acoustic Asynchronous Sensors. IEEE Sens. J. 2010, 10, 1469–1478. [Google Scholar] [CrossRef]
Mahapatra, C.; Mohanty, A.R. Explosive sound source localization in indoor and outdoor environments using modified Levenberg Marquardt algorithm. Measurement 2022, 187, 110362. [Google Scholar] [CrossRef]
Weiss, A.J. Direct position determination of narrowband radio frequency transmitters. IEEE Signal Process. Lett. 2004, 11, 513–516. [Google Scholar] [CrossRef]
Weiss, A.J. Direct Geolocation of Wideband Emitters Based on Delay and Doppler. IEEE Trans. Signal Process. 2011, 59, 2513–2521. [Google Scholar] [CrossRef]
Hahn, W.; Tretter, S. Optimum processing for delay-vector estimation in passive signal arrays. IEEE Trans. Inf. Theory 1973, 19, 608–614. [Google Scholar] [CrossRef]
Knapp, C.; Carter, G. The generalized correlation method for estimation of time delay. IEEE Trans. Acoust. Speech Signal Process. 1976, 24, 320–327. [Google Scholar] [CrossRef]
Schmidt, R. Multiple emitter location and signal parameter estimation. IEEE Trans. Antennas Propag. 1986, 34, 276–280. [Google Scholar] [CrossRef]
Swindlehurst, A.L.; Stoica, P. Maximum likelihood methods in radar array signal processing. Proc. IEEE 1998, 86, 421–441. [Google Scholar] [CrossRef]
Foy, W.H. Position-Location Solutions by Taylor-Series Estimation. IEEE Trans. Aerosp. Electron. Syst. 1976, AES-12, 187–194. [Google Scholar] [CrossRef]
Zhou, Z.; Rui, Y.; Cai, X.; Lu, J. Constrained total least squares method using TDOA measurements for jointly estimating acoustic emission source and wave velocity. Measurement 2021, 182, 109758. [Google Scholar] [CrossRef]
Smith, J.; Abel, J. Closed-form least-squares source location estimation from range-difference measurements. IEEE Trans. Acoust. Speech Signal Process. 1987, 35, 1661–1669. [Google Scholar] [CrossRef]
Chan, Y.T.; Ho, K.C. A simple and efficient estimator for hyperbolic location. IEEE Trans. Signal Process. 1994, 42, 1905–1915. [Google Scholar] [CrossRef]
Mahajan, A.; Walworth, M. 3D position sensing using the differences in the time-of-flights from a wave source to various receivers. IEEE Trans. Robot. Autom. 2001, 17, 91–94. [Google Scholar] [CrossRef]
Gillette, M.D.; Silverman, H.F. A Linear Closed-Form Algorithm for Source Localization From Time-Differences of Arrival. IEEE Signal Process. Lett. 2008, 15, 1–4. [Google Scholar] [CrossRef]
Liu, N.; Xu, Z.; Sadler, B.M. Low-Complexity Hyperbolic Source Localization With A Linear Sensor Array. IEEE Signal Process. Lett. 2008, 15, 865–868. [Google Scholar] [CrossRef]
Mikhalev, A.; Ormondroyd, R.F. Comparison of Hough Transform and Particle Filter Methods of Emitter Geolocation using Fusion of TDOA Data. In Proceedings of the 2007 4th Workshop on Positioning, Navigation and Communication, Hannover, Germany, 22 March 2007; pp. 121–127. [Google Scholar] [CrossRef]
Khalaf-Allah, M. Particle Filtering for Three-Dimensional TDoA-Based Positioning Using Four Anchor Nodes. Sensors 2020, 20, 4516. [Google Scholar] [CrossRef]
Lédeczi, Á.; Nádas, A.; Völgyesi, P.; Balogh, G.; Kusy, B.; Sallai, J.; Pap, G.; Dóra, S.; Molnár, K.; Maróti, M.; et al. Countersniper system for urban warfare. ACM Trans. Sens. Netw. 2005, 1, 153–177. [Google Scholar] [CrossRef]
Mikhalev, A.; Hughes, E.J.; Ormondroyd, R.F. Comparison of Hough Transform and particle filter methods of passive emitter geolocation using fusion of TDOA and AOA data. In Proceedings of the 2010 13th International Conference on Information Fusion, Edinburgh, UK, 26–29 July 2010; pp. 1–8. [Google Scholar] [CrossRef]
Compagnoni, M.; Pini, A.; Canclini, A.; Bestagini, P.; Antonacci, F.; Tubaro, S.; Sarti, A. A Geometrical–Statistical Approach to Outlier Removal for TDOA Measurements. IEEE Trans. Signal Process. 2017, 65, 3960–3975. [Google Scholar] [CrossRef]
Apolinário, J.A.; Yazdanpanah, H.; Nascimento, A.S.; de Campos, M.L.R. A Data-selective LS Solution to TDOA-based Source Localization. In Proceedings of the ICASSP 2019—2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Brighton, UK, 12–17 May 2019; pp. 4400–4404. [Google Scholar] [CrossRef]
Hough, P.V.C. A Method and Means for Recognizing Complex Patterns. U.S. Patent 3,069,654, 18 December 1962. [Google Scholar]
Rosenfeld, A. Picture processing by computer. ACM Comput. Surv. 1969, 1, 147–176. [Google Scholar] [CrossRef]
Mukhopadhyay, P.; Chaudhuri, B.B. A survey of Hough Transform. Pattern Recognit. 2015, 48, 993–1010. [Google Scholar] [CrossRef]
Shyu, H.-J. Application of the Bearing Trace Hough Transform (BTHT) to Passive Shipping Lane Monitoring. In Naval Research Laboratory Report NRL/FR/5580--98-9877; Naval Research Laboratory: Washington, DC, USA, 1998. [Google Scholar]
Luo, Y.; Sallai, J.; Ledeczi, A. A Hough Transform based approach to acoustic shot localization. In Proceedings of the 2016 19th International Conference on Information Fusion (FUSION), Heidelberg, Germany, 5–8 July 2016; pp. 1420–1426. [Google Scholar]
Filos, J.; Canclini, A.; Thomas, M.R.P.; Antonacci, F.; Sarti, A.; Naylor, P.A. Robust inference of room geometry from acoustic measurements using the hough transform. In Proceedings of the 2011 19th European Signal Processing Conference, Barcelona, Spain, 9 August–2 September 2011; pp. 161–165. [Google Scholar]
Kiryati, N.; Eldar, Y.; Bruckstein, A.M. A probabilistic Hough Transform. Pattern Recognit. 1991, 24, 303–316. [Google Scholar] [CrossRef]
Xu, L.; Oja, E.; Kultanen, P. A new curve detection method: Randomized Hough transform (RHT). Pattern Recognit. Lett. 1990, 11, 331–338. [Google Scholar] [CrossRef]
Li, H.; Lavin, M.A.; Le Master, R.J. Fast Hough transform: A hierarchical approach. Comput. Vis. Graph. Image Process. 1986, 36, 139–161. [Google Scholar] [CrossRef]
Levanon, N. Lowest GDOP in 2-D scenarios. IEE Proc.-Radar Sonar Navig. 2000, 147, 149–155. [Google Scholar] [CrossRef]
Isaacs, J.T.; Klein, D.J.; Hespanha, J.P. Optimal sensor placement for time difference of arrival localization. In Proceedings of the 48th IEEE Conference on Decision and Control (CDC), Shanghai, China, 15–18 December 2009; pp. 7878–7884. [Google Scholar] [CrossRef]

Figure 1. (a) Measurement scenario with emitter

E

and sensors

S_{i}

; (b) estimation process.

Figure 1. (a) Measurement scenario with emitter

E

and sensors

S_{i}

; (b) estimation process.

Figure 2. HHT of a scenario with emitter

E

(blue cross) and sensors

S_{1}

(reference),

S_{2}, S_{3}

, and

S_{4}

(red crosses), placed on the same plane. Hyperbolas generated by

S_{i}

and

S_{1}

are denoted by

H_{i, 1}

. (a) Source at a position with good GDOP and (b) bad GDOP.

Figure 2. HHT of a scenario with emitter

E

(blue cross) and sensors

S_{1}

(reference),

S_{2}, S_{3}

, and

S_{4}

(red crosses), placed on the same plane. Hyperbolas generated by

S_{i}

and

S_{1}

are denoted by

H_{i, 1}

. (a) Source at a position with good GDOP and (b) bad GDOP.

Figure 3. Hierarchical calculation of the HHT. Red and grey circles indicate centers of promising and unpromising tiles, respectively. Grey squares are the tiles. (a) Start of iteration step

n

, (b) after the pruning in iteration step

n

, (c) start of iteration step

n

+ 1.

Figure 3. Hierarchical calculation of the HHT. Red and grey circles indicate centers of promising and unpromising tiles, respectively. Grey squares are the tiles. (a) Start of iteration step

n

, (b) after the pruning in iteration step

n

, (c) start of iteration step

n

+ 1.

Figure 4. The derivation of the upper bound in a tile of center

P

and size

D_{n}

Figure 4. The derivation of the upper bound in a tile of center

P

and size

D_{n}

Figure 5. Sensor placement in the simulation setup with 7 sensors.

Figure 6. The HHT. (a) Near-range example, (b) long-range example.

Figure 7. Operation of the F-HHT. Sensor positions are shown by blue circles. The source position and the estimated position are shown by a blue cross and a green x, respectively. The centers of promising tiles are shown by red dots.

n

: iteration,

D_{n}

: grid size,

N_{t i l e}

: number of promising tiles.

n

: iteration,

D_{n}

: grid size,

N_{t i l e}

: number of promising tiles.

Figure 8. The near-range and the long-range experiments, conducted with a distance noise standard deviation of

σ_{i, 1} = 1 m

Figure 8. The near-range and the long-range experiments, conducted with a distance noise standard deviation of

σ_{i, 1} = 1 m

Figure 9. Sensor setup with 57 sensors.

Figure 10. RMSE of LS and F-HHT as a function of number of outliers.

Figure 11. HHT and positioning error, as a function of reference sensor index.

Table 1. Sensor placement in the 7-sensor experiment.

Sensor ID	Position (m)
1	(10, 10, 100)
2	(430, 190, 50)
3	(450, 350, 200)
4	(70, 300, 150)
5	(400, 10, 150)
6	(200, 50, 50)
7	(250, 470, 20)

Table 2. Parameters of the simulations.

Scenario	Source Position	$σ_{i, 1}$	$D_{f i n a l}$
#1	(157, 114)	0.1 m	0.1 m
#2	(157, 114)	1 m	0.1 m
#3	(157, 114)	1 m	1 m
#4	$(250, 1280)$	0.1 m	0.1 m
#5	$(250, 1280)$	1 m	1 m
#6	$(250, 1280)$	1 m	5 m

Table 3. Comparison of localization errors.

Scenario	√CRLB	RMSE
Scenario	√CRLB	LS	R-HHT	H-HHT	F-HHT
#1	0.064 m	0.071 m	0.080 m	0.079 m	0.075 m
#2	0.64 m	0.69 m	0.72 m	0.71 m	0.70 m
#3	0.64 m	0.68 m	0.80 m	0.82 m	0.81 m
#4	1.98 m	2.03 m	2.06 m	2.32 m	2.04 m
#5	19.8 m	21.0 m	21.3 m	21.6 m	21.1 m
#6	19.8 m	21.0 m	21.3 m	21.8 m	21.2 m

Table 4. Comparison of evaluation times.

Scenario	Run Time (ms)
Scenario	LS	R-HHT	H-HHT	F-HHT
#1	5.4	$59 \times 10^{3}$	$4.3 \times 10^{3}$	2.6
#2	4.5	2.6 $\times 10^{3}$	66	3.4
#3	4.4	$539$	48	2.3
#4	4.4	$59 \times 10^{3}$	$839$	41
#5	3.7	513	30	16
#6	3.8	513	6.3	7.7

Table 5. Trial numbers of R-HHT, H-HHT, and F-HHT.

Scenario	$N_{g r i d}$		Trials
Scenario	$N_{g r i d}$	R-HHT	H-HHT	F-HHT
#1	$2.5 \times 10^{7}$	$7 \times 10^{7}$	$5 \times 10^{6}$	$2.5 \times 10^{3}$
#2	$2.5 \times 10^{7}$	$3 \times 10^{6}$	$4 \times 10^{4}$	$1.9 \times 10^{3}$
#3	$2.5 \times 10^{5}$	$6 \times 10^{5}$	$2 \times 10^{4}$	$1.3 \times 10^{3}$
#4	$2.25 \times 10^{8}$	$7 \times 10^{7}$	$1 \times 10^{6}$	$7.1 \times 10^{4}$
#5	$2.25 \times 10^{6}$	$6 \times 10^{5}$	$3 \times 10^{4}$	$2.8 \times 10^{4}$
#6	$9 \times 10^{4}$	$6 \times 10^{5}$	$6 \times 10^{3}$	$1.2 \times 10^{4}$

Table 6. Localization error of LS and F-HHT in the presence of outlier measurements.

Scenario	$N_{o u l i e r}$	RMSE (m)
Scenario	$N_{o u l i e r}$	LS	F-HHT
#1	0	0.07	0.07
	1	2.0	0.08
	2	2.9	0.10
	3	3.5	1.4
	4	4.1	11.8
#4	0	2.1	2.1
	1	80	2.8
	2	120	7.5
	3	143	41
	4	181	193

Table 7. Localization errors of LS and F-HHT, using real measurements.

True Position	Estimated Position		Estimation Error (m)
True Position	LS	F-HHT	LS	F-HHT
(33.37, 48.24)	(35.93, 41.94)	(33.44, 47.81)	6.80	0.43
(31.94, 57.34)	(39.92, 45.81)	(32.19, 57.81)	14.03	0.53
(34.61, 73.91)	(44.83, 52.51)	(34.69, 73.44)	23.71	0.48
(28.61, 75.62)	(34.18, 38.85)	(29.06, 75.31)	37.19	0.55

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Simon, G.; Leitold, F. Passive TDOA Emitter Localization Using Fast Hyperbolic Hough Transform. Appl. Sci. 2023, 13, 13301. https://doi.org/10.3390/app132413301

AMA Style

Simon G, Leitold F. Passive TDOA Emitter Localization Using Fast Hyperbolic Hough Transform. Applied Sciences. 2023; 13(24):13301. https://doi.org/10.3390/app132413301

Chicago/Turabian Style

Simon, Gyula, and Ferenc Leitold. 2023. "Passive TDOA Emitter Localization Using Fast Hyperbolic Hough Transform" Applied Sciences 13, no. 24: 13301. https://doi.org/10.3390/app132413301

APA Style

Simon, G., & Leitold, F. (2023). Passive TDOA Emitter Localization Using Fast Hyperbolic Hough Transform. Applied Sciences, 13(24), 13301. https://doi.org/10.3390/app132413301

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

Passive TDOA Emitter Localization Using Fast Hyperbolic Hough Transform

Abstract

1. Introduction

2. Related Work

3. Fast Hyperbolic Hough Transform

3.1. Problem Formulation

3.2. Hyperbolic Hough Transform

3.3. Fast HHT

3.4. Global Convergence

4. Performance Evaluation

4.1. Test Simulation Setup

4.2. Illustration of the Operation

4.3. Error Analysis

4.4. Evaluation Time

4.5. Fault Tolerance

4.6. Real Measurements

5. Conclusions

Author Contributions

Funding

Institutional Review Board Statement

Informed Consent Statement

Data Availability Statement

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI