Passive localization and tracking of a mobile emitter, and joint learning of its reverberant three-dimensional (3D) acoustic environment, where critical structural features are unknown, is a key open problem. Unaccounted-for occluders are potentially present, so that the emitter can lose line-of-sight to the receivers, and can only be observed through its reflected raypaths. The locations of reflective boundaries must therefore be jointly estimated with the emitter's position. A multistage global optimization and tracking architecture is developed to solve this problem with a relatively unconstrained model. Each stage of this architecture establishes domain knowledge such as synchronization and initial environment estimation, which are inputs for the following stages of more refined algorithms. This approach is generalizable to different physical scales and modalities and improves on methods that do not exploit the motion of the emitter. In one stage of this architecture, particle swarm optimization is used to simultaneously estimate the environment and the emitter location. In another stage, a Hough transform-inspired boundary localization algorithm is extended to 3D settings, to establish an initial estimate of the environment. The performance of this holistic approach is analyzed and its reliability is demonstrated in a reverberant watertank testbed, which models the shallow-water underwater acoustic setting.

Passive underwater acoustic localization in shallow-water settings is one of the most important and challenging localization problems (Fischell , 2019; Waterston , 2019). The emitter to be localized could be a surface vessel or underwater vehicle, and the receivers could be surreptitiously placed for passively monitoring a body of water. Receivers are typically omnidirectional tethered hydrophones, but could also be placed on a set of mobile underwater vehicles, as illustrated in Fig. 1. When the environment is accurately known, state-of-the-art matched-field processing (MFP) methods can be used for localization (Baggeroer , 1993). The standard MFP framework does not account for the presence of man-made features such as pier pilings or anchored vessels, or unsurveyed reefs (Chitre and Pelekanakis, 2014). Depending on their positions, these large and potentially unknown objects can act as additional reflective boundaries, or as occluders that block the line-of-sight (LOS) arrivals to receivers. MFP is sensitive to model mismatch (Le Gall , 2016; Xu , 2006) and its performance could therefore deteriorate due to such objects, motivating the formulation of a localization framework that incorporates their presence. A representative scenario is presented in Fig. 2, where the goal is to estimate the unknown source position at each transmission point along a trajectory (red circles), in the presence of a large occluder, in a scaled model of the shallow-water setting. This is an instance of an open problem that also arises for indoor acoustic or radio frequency localization (Dokmanic , 2015), where the non-line-of-sight (NLOS) arrivals due to multipath reflections must be judiciously exploited rather than mitigated, with a robust passive localization method.

FIG. 1.

(Color online) Typical shallow-water underwater acoustic localization problem setting [image is modified from the original in Kim (2019); licensed under a Creative Commons Attribution 4.0 License].

FIG. 1.

(Color online) Typical shallow-water underwater acoustic localization problem setting [image is modified from the original in Kim (2019); licensed under a Creative Commons Attribution 4.0 License].

Close modal
FIG. 2.

(Color online) Top-down view of the test environment with emitter's trajectory (a), and the actual experimental setup (b). This is a 3D localization and tracking problem with an occluding cylinder, and is more challenging than settings where the receivers surround the emitter (Setlur , 2012).

FIG. 2.

(Color online) Top-down view of the test environment with emitter's trajectory (a), and the actual experimental setup (b). This is a 3D localization and tracking problem with an occluding cylinder, and is more challenging than settings where the receivers surround the emitter (Setlur , 2012).

Close modal

Without restrictions on the localization setting, such problems are ill-posed and admit infinitely many solutions. The static nature of the environment, and the spatial continuity in the trajectory of a moving emitter, are key to unambiguous localization. Some of the important aspects of the setting are modeled as in Table I; the limited information here leaves some key parameters undetermined so that it is natural to adopt a multistage approach to first estimate them. Accordingly, we deploy a series of localization algorithms sequentially, as depicted in Fig. 3, using an architecture that we call passive end-to-end localization (PEEL). Each stage of PEEL produces approximate solutions to different aspects of the localization problem, which then become inputs for the following stages. We call this an end-to-end method because it starts from received waveforms that do not have a common clock with the emitter, and delivers a joint estimate of the full emitter trajectory and the environment. We ultimately demonstrate that PEEL achieves strong performance, even when there is modeling mismatch that violates the assumptions in Table I.

TABLE I.

Model of the localization setting.

Problem feature Modeling assumptions
Speed of sound, vs  Known, and constant within the environment 
Environment  Static; only the emitter position is changing 
Reflectors, η i  Planar, known number, unknown position and orientation 
Transmissions  Periodic transmissions of deterministic pulses, with known period 
Receivers, p r , j  At least 4 receivers deployed, with known positions 
Problem feature Modeling assumptions
Speed of sound, vs  Known, and constant within the environment 
Environment  Static; only the emitter position is changing 
Reflectors, η i  Planar, known number, unknown position and orientation 
Transmissions  Periodic transmissions of deterministic pulses, with known period 
Receivers, p r , j  At least 4 receivers deployed, with known positions 
FIG. 3.

(Color online) The sequence of algorithms in the PEEL method. The sections with green outlines incorporate novel developments or applications of those algorithms.

FIG. 3.

(Color online) The sequence of algorithms in the PEEL method. The sections with green outlines incorporate novel developments or applications of those algorithms.

Close modal

While there is currently no end-to-end passive localization method such as the one we propose, similar frameworks have been designed for different settings. Joint acoustic tracking and environment learning problems can be solved using simultaneous localization and mapping (SLAM) methods (Li and Krolik, 2013; Ribeiro , 2010) in active settings, or when the emitter and receivers form a network (Naseri and Koivunen, 2016). Similarly, while there are time-difference-of-arrival (TDOA) and time-of-arrival (TOA) methods for localization (Scheuing and Yang, 2008; Venkateswaran and Madhow, 2012; Alameda-Pineda and Horaud, 2014), they do not account for unknown occluders as PEEL does.

Our main contribution is an architecture for joint localization and environment learning. Among other tools, PEEL makes use of:

  • A global optimization method for joint boundary estimation and localization, which does not need to be tuned to a specific problem.

  • The mapping of unknown occluders after tracking, followed by iterations with a refined optimization metric, which has not been implemented for passive localization.

  • Closed-form derivations for some key steps of a boundary localization technique (Naseri and Koivunen, 2016), that improve computational efficiency and accuracy.

The paper is structured as follows. In Sec. II, we present the signal model and geometric considerations used for localization. We present the PEEL method in Sec. III. The results that are achieved in the scale watertank environment of Fig. 2 are given in Sec. IV. In Sec. V, we conclude with potential avenues for future research.

In our mathematical notation throughout the paper, lowercase bold variables such as p are vectors, and uppercase calligraphic characters such as P denote sets.

In this section, we introduce the notation and signal model, and present key geometrical considerations that are necessary to frame the localization problem. We specifically focus on how multipath arrivals can be modeled and used to improve localization performance.

Triangulation is a standard localization principle that models the emitter as the point of intersection of circles of equal time-of-flight (TOF) from the emitter to each receiver. Since these TOF arcs do not perfectly intersect in practice (e.g., due to noise), an estimation method has to be used; often, this method has to simultaneously solve a data association problem (Too , 2017). When the locations of reflectors are known, multipath can be leveraged for better accuracy. Reflections can be interpreted as producing ‘virtual receivers’ as illustrated in Fig. 4, which corresponds to the structure of our watertank setting.

FIG. 4.

(Color online) Illustration of virtual receivers produced by reflective planar boundaries in an occluder-free environment (a), and the corresponding noise-free simulated received signal for this scenario (b). A NLOS TOF to the real receiver is equivalent to the LOS TOF to the virtual receiver.

FIG. 4.

(Color online) Illustration of virtual receivers produced by reflective planar boundaries in an occluder-free environment (a), and the corresponding noise-free simulated received signal for this scenario (b). A NLOS TOF to the real receiver is equivalent to the LOS TOF to the virtual receiver.

Close modal
We make the simplification that reflective structures are (piecewise) approximated as planar boundaries (Deane, 1994), characterized by the range ρ, azimuth θ, and elevation ϕ of their normal vector relative to the origin. Thus, boundary i is parametrized as the vector η i = [ ρ i , θ i , ϕ i ] T. We denote the Cartesian coordinates of the emitter and receiver positions as p e 3 and p r 3, respectively. The LOS TOF, τ0, is calculated as
τ 0 = | | p r p e | | 2 v s ,
(1)
where vs is the constant speed of sound, which can be a valid approximation for well-mixed shallow-water setting at short ranges and high acoustic frequencies (Chitre, 2007; Too , 2017). For boundary i, we obtain the virtual receiver location p i by finding the corresponding reflection of p r. The TOF to the receiver from the NLOS arrival of boundary i is equal to the TOF from p e to the corresponding p i; thus
τ i = | | p i p e | | 2 v s .
(2)
We denote a single copy of the periodically emitted signal as s(t), and we merge the effects of attenuation and reflection as the equivalent attenuation coefficient αi for each path. The received signal r(t) is then
r ( t ) = i = 0 N α i s ( t τ i ( p e , η i ) ) + n ( t ) ,
(3)
where n(t) is an unknown random noise realization, and N is the (known) number of reflective surfaces. Note that αi can be zero, as when an occluder blocks a particular arrival, so Eq. (3) incorporates occlusion. We denote each receiver as p r , j, each of which has TOFs τ i , j from their own virtual receivers, and received signals r j ( t ) [defined as in Eq. (3)]. We only model the single-bounce reflections, which are typically the strongest multipath components. In practice, not only will there be reverberation, but the N that is assumed known may be inaccurate. A proposed localization method must be robust to the resulting discrepancies.

We model s(t) as pulsed and periodic, which limits our estimation scenarios to such man-made sources. This model allows us to record s(t) and to infer some key parameters: the approximate (truncated) channel impulse response length Tc (Sun and Wang, 2018), and the time period Tp between emissions. We require s(t) to meet the following constraints in order to have accurate estimates { τ ̂ i , j } of the { τ i , j }. First, we define Ts as the temporal duration of s(t), so that s(t) is strictly zero outside 0 t T s. Next, we define TA to be the full width at half maximum (FWHM) of the autocorrelation function of s(t), which is the time resolution of the matched-filtered emitted signal. Finally, we define Tm as the minimum time difference between the main multipath arrivals. The PEEL operation modes are then as follows:

  • If Ts < Tm, we have a pulse that is shorter than the multipath separation, as is the case with the Gaussian pulse in our experiments. Here, we do not need to know s(t) for PEEL; Ts is short enough that the multipath arrivals are roughly separable.

  • If Ts > Tm and TA < Tm, then after matched-filtering, the multipath peaks are roughly separable. Now, we do in fact need to know s(t) for this matched-filtering step, so that we are able to identify the distinct arrivals and thereby solve for the environment.

We also need Tp > Tc to avoid different periodic emissions being observed within the same r(t).1 With these introductory developments, we are ready to present an architecture for reliable localization in the presence of such errors.

We now present our proposed joint localization and structure learning method, PEEL. First, we use TDOA and TOA localization for initialization; then obtain an initial boundary localization estimate; and finally track a moving emitter as we simultaneously learn the environment, as summarized in Fig. 5.

FIG. 5.

(Color online) Conceptual flow chart for the proposed PEEL architecture.

FIG. 5.

(Color online) Conceptual flow chart for the proposed PEEL architecture.

Close modal

We first describe a standard TDOA localization algorithm (which does not require knowledge of { η i }) for synchronization with the emitter (Korhonen, 2008), then a TOA algorithm that leverages multipath from known boundaries { η i } after synchronization (Ribeiro , 2010; Brutti , 2010). Both algorithms use a grid-search for the emitter to obtain p ̂ e, where a set P of grid points p cand is defined over the target area.2 Each p cand P corresponds to a potential emitter location, with an associated set of { τ i , j , cand } to the receivers p r , j and the corresponding virtual receivers p i , j. Thus, each τ i , j , cand is a function of p r , j and of η i. The TDOA and TOA stages are used to initialize more refined techniques.

TDOA does not require a common clock with the emitter and is suitable to establish synchronization (Su , 2017; Qiao , 2013). Denoting the calculated cross correlation between the signals at receivers j1 and j2 as R j 1 , j 2 ( τ ), the TDOA metric that we use is
C TDOA ( τ cand ( p cand ) ) j 2 j 1 j 1 = 1 M | R j 1 , j 2 ( τ j 1 , cand ( p r , j 1 ) τ j 2 , cand ( p r , j 2 ) ) | .
(4)
Using the generalized cross correlation with phase transform (GCC-PHAT) method for better performance (Knapp and Carter, 1976), we estimate the emitter location p ̂ e , TDOA as
p ̂ e , TDOA = argmax p cand P C TDOA ( τ cand ( p cand ) ) .
(5)
The average difference between LOS peak delays and the expected TOFs from TDOA is then used to establish synchronization.3 One can then use more accurate TOA localization methodologies. We define the additive TOA metric C TOA ( τ cand ( p cand ) ) as
C TOA ( τ cand ( p cand ) ) j = 1 M i = 0 N | r j ( t τ i , j , cand ( p r , j , η i ) ) | ,
(6)
where the time indices ( t τ i , j , cand ( p r , j , η i ) ) are rounded to the nearest sample in a discrete-time implementation. The estimate p ̂ e , TOA is then obtained as
p ̂ e , TOA = argmax p cand P C TOA ( τ cand ( p cand ) ) .
(7)

The heuristic metric in Eq. (6) is maximized when each τ i , j corresponds to a peak location; such locations are the intersections of TOF arcs, as seen in Fig. 6(a). When the { η i } are known, and therefore the virtual receivers are labeled correctly, p ̂ e is obtained accurately.

FIG. 6.

(Color online) A TOA grid-search that localizes the emitter on real data (a), and a grid-search that fails at a different location due to the alternative arrival arc crossings (b), at ∼36 dB SNR.

FIG. 6.

(Color online) A TOA grid-search that localizes the emitter on real data (a), and a grid-search that fails at a different location due to the alternative arrival arc crossings (b), at ∼36 dB SNR.

Close modal

TDOA and TOA localization are more robust if the search space P can be limited to the vicinity of the correct location, p e. Consider the results in Fig. 6(a); although p e is estimated accurately, there are many local maxima in the search space. These can lead to large errors as in Fig. 6(b). Confining P to the vicinity of p e avoids such ambiguities. We summarize our TDOA and TOA localization algorithms in Algorithms 1 and 2.

Algorithm 1.

TDOA localization for initialization.

Require: Unoccluded LOS arrivals at each receiver. 
  Input: P , { r j ( t ) } j = 1 M , { p r , j } j = 1 M  
  for p cand P do  
   Calculate C TDOA [Eq. (4)
  end  
  Get p ̂ e , TDOA[Eq. (5)
  Result: p ̂ e , TDOA  
Require: Unoccluded LOS arrivals at each receiver. 
  Input: P , { r j ( t ) } j = 1 M , { p r , j } j = 1 M  
  for p cand P do  
   Calculate C TDOA [Eq. (4)
  end  
  Get p ̂ e , TDOA[Eq. (5)
  Result: p ̂ e , TDOA  
Algorithm 2.

TOA localization after initialization.

Require: Time synchronization after applying Algorithm 1; number of boundaries N is known.4 
  Input: P , { r j ( t ) } j = 1 M , { p r , j } j = 1 M , { η i } i = 1 N  
  for p cand P do  
   Calculate C TOA [Eq. (6)
  end  
  Get p ̂ e , TOA [Eq. (7)
  Result: p ̂ e , TOA  
Require: Time synchronization after applying Algorithm 1; number of boundaries N is known.4 
  Input: P , { r j ( t ) } j = 1 M , { p r , j } j = 1 M , { η i } i = 1 N  
  for p cand P do  
   Calculate C TOA [Eq. (6)
  end  
  Get p ̂ e , TOA [Eq. (7)
  Result: p ̂ e , TOA  

One advantage of Algorithm 2 is that it is a heuristic weighting of the arrivals. A localization method that is solely based on { τ i , j } will assign equal importance to each arrival so that there will be no automatic compensation for faulty TOF estimates. A TOF estimate is likelier to have a high error if the corresponding arrival is weak. In such cases, its contribution to C TOA in Algorithm 2 will be smaller, and its impact on p ̂ e , TOA in Eq. (7) will be reduced.

Algorithm 2 also allows us to compare different hypothesized environments. Given candidate boundaries { η i } cand, we can run Algorithm 2 for each { η i } cand, and pick the maximum metric solution as our joint estimate of p ̂ e , TOA and { η i }. C TOA is therefore a heuristic measure of the estimate's goodness. Recall that the { η i } cannot generally be assumed known; thus, we next propose a boundary localization algorithm to produce a suitable set of { η i } cand.

After TDOA and TOA localization, we have an initial estimate of the emitter's position and coarse synchronization. To jointly estimate the emitter and the environment in a computationally feasible manner, we need a good initial estimate of reflective boundary positions, which ultimately yields a tight search space for { η i }. We assume that a set of unlabeled { τ ̂ i , j } has been estimated (Demirli and Saniie, 2001), though some of these estimated arrival times might have large errors. We now develop a Hough transform-inspired method to estimate { η i } which is robust to errors in { τ ̂ i , j }, building on (Naseri and Koivunen, 2016).

As some background context, Euclidean distance matrices (EDM) (Dokmanic , 2015) and convex optimization methods (Naseri , 2014) have been used to solve for virtual emitter positions and thereby obtain boundary positions. These methods can suffer from large errors if an arrival τ ̂ i , j has been labeled incorrectly, and their performances are initialization-sensitive. An alternative approach relies on the fact that in two dimensions, these τ ̂ i , j each define an ellipse whose foci are p e and p r, as in Fig. 7(a); in three dimesions, they define spheroids. In Naseri and Koivunen (2016), boundaries are observed to be common tangents to such ellipses; by fitting tangent planes to these ellipses, we can estimate the boundaries while avoiding the error-prone initialization and labeling procedure in multipath environments, such as in Fig. 7(b). It is convenient to define a given η i by the spherical coordinates of a vector from the origin (arbitrarily selected) that is orthogonal to the boundary plane. Thus, each boundary is specified by a range ρ, azimuth θ, and elevation ϕ, where these local spherical coordinates correspond to a standard east-north-up Cartesian coordinate system.

FIG. 7.

(Color online) The NLOS arrivals define ellipses/spheroids of equidistance (a); we illustrate the ellipses defined by { τ i , j } for our experimental setting (depicted in 2D for convenience). Each boundary is a common tangent to a single ellipse due to each receiver, highlighted here by matching colors.

FIG. 7.

(Color online) The NLOS arrivals define ellipses/spheroids of equidistance (a); we illustrate the ellipses defined by { τ i , j } for our experimental setting (depicted in 2D for convenience). Each boundary is a common tangent to a single ellipse due to each receiver, highlighted here by matching colors.

Close modal

The transform for boundary localization in Naseri and Koivunen (2016), which for convenience we refer to as common tangents to spheroids (COTANS), is derived as follows. A boundary defined by ρ, θ, and ϕ can be conceptualized as a point ( ρ , θ , ϕ ) in a COTANS transform domain; working out the ( ρ , θ , ϕ ) expression of a plane is to take the plane's COTANS transform (Borrmann , 2011). When the space ρ × θ × ϕ is discretized as a 3-mode tensor, each hypothetical tangent has a finite set of ( ρ , θ , ϕ ) on which it may lie. Incrementing this “accumulator” tensor over every potential ( ρ , θ , ϕ ) entry yields a COTANS-domain image [as in Fig. 9(a)], with maxima at the true boundaries in the noiseless case.

Here, we extend the method of Naseri and Koivunen (2016) to three dimensions, and we further derive a direct solution for the tangents in the COTANS domain, which precludes the need for randomly sampling points ( ρ , θ , ϕ ) on the surface of an ellipse/spheroid. We define a grid over θ and ϕ, and obtain the COTANS transform for the tangents at each θ and ϕ. In two dimensions, for a given θ, we obtain the COTANS transform's ρ in Appendix A 1 as
ρ ( θ ) = a 2 cos 2 θ + b 2 sin 2 θ ,
(8)
and in three dimensions, the result for a given θ and ϕ is calculated in Appendix A 2 as
ρ ( θ , ϕ ) = c 2 cos 2 ϕ + a 2 sin 2 ϕ ,
(9)
where a, b, and c are standard axes of spheroids calculated from the { τ ̂ i , j }. This result is obtained for an origin-centered spheroid. To move a point cloud of ( ρ , θ , ϕ ) centered on the origin to the true emitter and receiver positions, we rotate the points to match the true spheroid's orientation (adjusting their θ and ϕ to be some θ rot and ϕ rot), and then translate the points to the true spheroid's position (yielding a final ρ COTANS , θ COTANS, and ϕ COTANS). These COTANS-domain operations are illustrated in Fig. 8.
FIG. 8.

(Color online) Illustration of the steps to obtain the COTANS transform of a particular tangent line (depicted in 2D). Description of one tangent to a standard ellipse (a); rotation of this origin-centered ellipse and its tangent (b); and the translation of this ellipse to its real position (c).

FIG. 8.

(Color online) Illustration of the steps to obtain the COTANS transform of a particular tangent line (depicted in 2D). Description of one tangent to a standard ellipse (a); rotation of this origin-centered ellipse and its tangent (b); and the translation of this ellipse to its real position (c).

Close modal

In practice, we first generate a set of points { ( ρ , θ , ϕ ) } by defining an azimuth-elevation grid, and obtaining the COTANS transforms for each of the corresponding tangents. The resulting points { ( ρ , θ , ϕ ) } COTANS are then rounded to a desired accuracy.5 We define a tensor over ρ, θ, and ϕ with this resolution, for each rounded { ( ρ , θ , ϕ ) } COTANS point, we increment the corresponding cell in the tensor by 1. We sum the tensors for each NLOS arrival at each receiver, to obtain a final tensor of superimposed COTANS curves such as in Fig. 9(a). Each ellipse contributes a curve; observe that the curves do not intersect at a single point because of errors in { τ ̂ i , j } and in the measured p e and p r, so we have multiple maxima. We therefore use a smoothing filter to replace the tensor values with local averages.6 Then, we extract the locations of as many maxima as there are boundaries as in Fig. 9(b), where we set the neighborhood of each maximum in the tensor to zero, to avoid picking the same boundary multiple times.7 This method is summarized in Algorithm 3.

FIG. 9.

(Color online) The COTANS accumulator for inconsistent { τ ̂ i , j } (a), and its boundary estimates (b). The image is periodic in azimuth, so the maximum near 360° is close to the correct location at 0°.

FIG. 9.

(Color online) The COTANS accumulator for inconsistent { τ ̂ i , j } (a), and its boundary estimates (b). The image is periodic in azimuth, so the maximum near 360° is close to the correct location at 0°.

Close modal
Algorithm 3.

COTANS transform-based boundary localization.

Require: Time synchronization from Algorithm 2; number of boundaries N is known. 
  Input: { τ ̂ i , j } , { p r , j } j = 1 M , p ̂ e , TDOA  
  for τ ̂ i , j { τ ̂ i , j } do  
    Calculate COTANS transform for spheroid, add to accumulator ( Appendix A
  end  
  Result: Apply smoothing to COTANS accumulator, find its N maxima as { η ̂ i }
Require: Time synchronization from Algorithm 2; number of boundaries N is known. 
  Input: { τ ̂ i , j } , { p r , j } j = 1 M , p ̂ e , TDOA  
  for τ ̂ i , j { τ ̂ i , j } do  
    Calculate COTANS transform for spheroid, add to accumulator ( Appendix A
  end  
  Result: Apply smoothing to COTANS accumulator, find its N maxima as { η ̂ i }

Modeling the environment as consisting of planar boundaries makes environment learning more computationally tractable. Note that the virtual receiver model is not limited to planar boundaries: the planar boundary estimates can be tangents to a non-planar boundary. In this case, however, the boundary estimates will be inaccurate for localization purposes when the emitter moves to a different position. Therefore, rather than using the initial COTANS boundary estimates for the entire emitter trajectory, we need evolving boundary estimates over time, utilizing the global optimization method which we now discuss.

At this stage, we have initial estimates of both the emitter position and the boundaries, allowing for simultaneous optimization over p e and { η ̂ i }. We assume candidate boundaries { η i } cand that are in the vicinity of the { η ̂ i }, and then carry out TOA localization (Algorithm 2) for each { η i } cand. These { η i } cand are drawn from a search space ρ [ ρ i ρ 0 , ρ i + ρ 0 ] , θ [ θ i θ 0 , θ i + θ 0 ] , ϕ [ ϕ i ϕ 0 , ϕ i + ϕ 0 ], where ρi, θi, and ϕ i are estimated from Algorithm 3, and ρ0, θ0, and ϕ 0 are error margins obtained as the typical maximum parameter estimation errors from COTANS, for a given setting. The { η ̂ i } associated with the highest maximum metric C TOA [Eq. (6)] yields the most likely p e and { η i }. Interlacing this algorithm with tracking over the emitter positions, we discover and map the occluders in the environment, and refine our estimates for joint localization and environment learning.

In this simultaneous optimization problem with measurement errors, it is advantageous to maintain multiple evolving estimates of the emitter and boundary locations, for which we use particle swarm optimization (PSO) (Clerc, 2010). In this global optimization algorithm, each { η i } cand is termed a “particle.” Each particle is initialized at a random location and given an initial “velocity” in the search space. Algorithm 2 is carried out for each particle, and the maximum C TOA and corresponding p ̂ e , cand is recorded. Particle velocities are modified at each PSO iteration towards the global maximum metric particle (with the maximum C TOA), and also to the maximum metric particle in a random “neighborhood” of particles. Thus, particles explore the region around the current best approximation.

PSO was used in PEEL because its metaparameters do not need to be adapted to a specific physical setting and because it can be run without re-initialization during tracking. There are alternative global optimization algorithms that can also be used in PEEL, however, and it would be advantageous in an application setting to fine-tune the metaparameters of PSO or another algorithm for that particular environment.

In order to be able to accommodate for occlusion, we rely on a tracking procedure over the periodic transmissions of the moving emitter. Kalman tracking with a constant-acceleration model, for example, can be applied. The p ̂ e obtained using PSO is used as the measurement vector to update the tracking. Critically, the grid-search is confined to a sufficiently small area P around the resulting prediction. This way, even if the emitter is temporarily occluded and the underlying C TOA metric is erroneous, p ̂ e will still be close to the actual p e. Thus, if the emitter enters an area that is unfavorable for localization, PEEL can potentially mitigate the resulting estimation errors. Our algorithm for PSO localization and boundary estimation is given in Fig. 10, and our overall approach is summarized in Algorithm 4.

FIG. 10.

(Color online) Flowchart for simultaneous boundary and emitter localization with PSO.

FIG. 10.

(Color online) Flowchart for simultaneous boundary and emitter localization with PSO.

Close modal
Algorithm 4.

Joint PSO localization and Kalman tracking.

Input: { r j ( t ) } j = 1 M , { p r , j } j = 1 M , { η i } cand  
Require: Moving emitter. 
 for each emitter position in trajectory do 
  for pre-specified number of PSO iterations do 
   Perform PSO, get p ̂ e , PSO , { η ̂ i } (Fig. 10,  Appendix B
  end for  
  Use p ̂ e , PSO to update Kalman tracking, get new P for next iteration. 
 end for  
   Result: p ̂ e , PSO , { η ̂ i } 
Input: { r j ( t ) } j = 1 M , { p r , j } j = 1 M , { η i } cand  
Require: Moving emitter. 
 for each emitter position in trajectory do 
  for pre-specified number of PSO iterations do 
   Perform PSO, get p ̂ e , PSO , { η ̂ i } (Fig. 10,  Appendix B
  end for  
  Use p ̂ e , PSO to update Kalman tracking, get new P for next iteration. 
 end for  
   Result: p ̂ e , PSO , { η ̂ i } 

From the results of Algorithm 4, we can determine whether each LOS path to receivers has been attenuated beyond expectation. Occlusion can be determined by assuming a 1 / R ̂ attenuation factor in the received signal magnitude, where R ̂ is a given estimated LOS range. Unoccluded arrivals allow estimation of the emitted power, by scaling the corresponding LOS magnitudes by R ̂. If the received LOS magnitude is substantially below the magnitude predicted by geometric spreading, we conclude that an occluder is present.8 Volumetric intersection methods can then be used to map the occluders (Srinivasan , 1990).

After mapping the occluders, we re-run Algorithm 4. We now know which LOS arrivals are blocked by the occluder, and we change the additive metric C TOA by leaving blocked paths out of the sum in Eq. (6). Our flexible cost function thereby allows us to compensate for model mismatch due to occlusion, and yields increasingly improved performance.

We now evaluate our PEEL method for its passive localization and environment learning capabilities, seeking to demonstrate and understand its performance and robustness. For this purpose, we use the convenient setting of the watertank depicted in Fig. 2. We first test PEEL on real data obtained from the watertank; then add synthetic noise to this data to analyze PEEL's performance compared to other localization methods and finally conduct trials in a simulation setting with a wider range of localization scenarios.

In our experiments, we emit a high-frequency Gaussian pulse
s ( t ) = e t 2 f 0 2 cos ( 2 π f 0 t ) ,
(10)
where the center frequency f0 was set to 280 kHz, corresponding to a wavelength of 5 mm. The boundaries used for localization are the water surface, whose position is known since the receiver depths are known; and the bottom and three plastic sides of the tank, whose ranges to the origin and azimuths are respectively known to within ±5 mm and ± 6 °, as obtained using COTANS. One side-boundary is deliberately left unaccounted for, which presents the difficulty of having an unmodeled set of arrivals in the data so that we have N = 5. We use M = 4 receivers. The speed of sound through water in the tank was measured as 1485 m/s.

For our PEEL implementation, we use 128 particles and 5000 PSO iterations for each emitter position. At each PSO iteration, we use a grid P of 0.5 mm resolution and sides of 6 mm centered on the expected emitter location. We perform ten experiments for each localization scenario and obtain the root mean square (rms) localization error. The rms error for the first pass of PEEL is 3.9 mm, as seen in Fig. 11(a). To put this into perspective, the average distance of the emitter from the receivers for this trajectory is 93.5 mm, and the rms error in this scenario is 4.2% of the mean distance to the emitter. The final boundary estimates in Fig. 11(a) are obtained using the full estimated trajectory, which leads to a more accurate solution compared to the results obtained from using just one emitter position. Running PEEL's second pass, we get the result in Fig. 11(b), where the rms error is reduced to 2.7 mm (2.9% of the mean distance). The occluder is reconstructed with sufficient accuracy to improve our localization result. To understand how much the unknown occluder's presence affects the result, we experimented with the same scenario but with no occluder, obtaining an rms error of 1.5 mm (1.6% of the mean distance), which is within our margin of measurement accuracy. In these trials, occlusion and interference due to reverberation are the key error factors; the SNR of unoccluded LOS arrivals is > 30 dB.

FIG. 11.

(Color online) First-pass localization performance, with an occluding cylinder present in the environment (a), and the improved second-pass result, using the modified metric obtained by taking the occluder into account (b). The boundary that is unaccounted-for to introduce model mismatch is the bottom boundary in these top-view figures.

FIG. 11.

(Color online) First-pass localization performance, with an occluding cylinder present in the environment (a), and the improved second-pass result, using the modified metric obtained by taking the occluder into account (b). The boundary that is unaccounted-for to introduce model mismatch is the bottom boundary in these top-view figures.

Close modal

If synchronization and an initial boundary estimate are given to us, single-receiver localization using Algorithm 4 becomes feasible. In an experiment with a single emitter at the origin and no occluder, we get 2.1 mm of rms error (2.2% of the mean distance). This strong performance is surprisingly only slightly worse than when we use four receivers.

In order to analyze the performances of the localization algorithms that we have examined (TOA localization with known boundaries, and PSO), we carried out localization trials with synthetic Gaussian noise added to experimental signals. For the first position in the trajectory, where there is LOS to every receiver, the same noise power in units of dBm was added to each recorded signal. In Fig. 12(a), we compare PSO's performance for different numbers of PSO iterations, vs that of TOA localization with known boundaries. Because PSO optimizes for the boundary locations as well, its high runtime led to us conduct 150 trials for each noise level, while 5000 trials were carried out for TOA. There is a performance gap between PSO and the “ground truth” TOA result, which is smaller for higher noise levels. Longer runtimes improve PSO's performance and help close this gap.

FIG. 12.

(Color online) Performance of PSO localization versus TOA with known boundaries, highlighting better accuracy with more PSO iterations (a); and rms boundary localization error performances for experimental results with synthetic added noise (b).

FIG. 12.

(Color online) Performance of PSO localization versus TOA with known boundaries, highlighting better accuracy with more PSO iterations (a); and rms boundary localization error performances for experimental results with synthetic added noise (b).

Close modal

We next evaluate the performance of environment learning by adding synthetic Gaussian noise to the signals measured at the last position in the trajectory (where two LOS arrivals are occluded) and obtaining the rms range error from the origin to all the boundaries, averaged over trials and the number of boundaries. We conducted 100 such trials for each noise level and for three scenarios: in the absence of the occluding cylinder as measured in a separate set of experiments, with the occluding cylinder present but not accounted for, and with the occluded arrivals removed in the revised metric (modeling the PEEL second pass). The results are given in Fig. 12(b), where we observe a gap in performance caused by the occlusion, beyond the performance deterioration caused by a lower effective SNR. Boundary estimation is improved by compensating for the occluder, though this improvement is less dramatic than the performance gains in emitter localization. Adding a noise power greater than 20 dBm causes much larger errors in the localization and boundary estimation tasks in this setting, as the noise peaks have much greater magnitude than the NLOS arrivals.

To test PEEL's performance on a wider range of experimental geometries, we generated synthetic data that modeled the watertank setting with randomized receiver locations. We implemented a simulation setting as per Deane (1994) that modeled a reverberant environment with multibounce arrivals, producing synthetic received signals with the mirror image method. In these trials, the emitter's trajectory, the boundary locations, the occluder, and the position of the receiver at the origin were kept constant, and PEEL was run for 10 different deployments of the other three receivers. These receivers were distributed uniformly over the region bounded by the positions of receivers in the real experimental setting. We used the same Gaussian pulse that was used in the watertank experiment, though the only distortion that these pulses underwent in simulation was attenuation. Ultimately, we obtained an rms error of 2.6 mm over ten trials, which was close to the watertank results.

Overall, we observed reliable performance across all of the trials for this particular demonstration. Once the unknown occluder is compensated-for, we only observe several millimeters of error, which is reasonable since there are already millimeter-scale errors in the assumed receiver positions. These results are summarized in Table II.

TABLE II.

Localization experiment results.

Scenario rms Error (mm)
Four receivers, no occluder  1.5 
Single receiver, no occluder  2.1 
Four receivers, occluder, first-pass  3.9 
Four receivers, occluder, final result  2.7 
Simulation, random receivers  2.6 
Scenario rms Error (mm)
Four receivers, no occluder  1.5 
Single receiver, no occluder  2.1 
Four receivers, occluder, first-pass  3.9 
Four receivers, occluder, final result  2.7 
Simulation, random receivers  2.6 

In this paper, we developed a method for localization and environment learning, which was tested in a three-dimensional (3D) reverberant underwater environment with occluding objects. Our PEEL method performed iteratively refined localization, with each solution yielding results that served as assumed knowledge for the following more sophisticated algorithm stages. This multi-stage approach handled a complex problem where successively establishing synchronization, boundary estimation, emitter tracking, and occluder estimation, made solving the overall task computationally feasible. Applying these stages separately is suboptimal as compared to a hypothetical joint optimization over all of the aspects of the problem. However, PEEL breaks down the localization task into subproblems with well-posed solutions. We ultimately achieved accurate simulation and experimental results in end-to-end localization.

The trajectory to which PEEL was applied resulted in a difficult localization problem which nevertheless proved to be solvable. There are many adversarial scenarios, however, that can defeat our method. The emitter may move behind an occluder and come to a halt, or it may take a sharp turn which throws off the tracking algorithm. One way to cope with such scenarios would be to run several instances of PEEL in parallel. A discrepancy between their results would indicate catastrophic errors. This parallelization strategy for PEEL is also important for choosing the number of boundaries, N. If one is not confident about the assumption on N, then running PEEL in parallel with different values of N becomes advantageous. Choosing the best N is an important open problem for PEEL.

Note that we do not require a complete deduction of the occluders for improved localization performance. In this passive setting, a complete occluder mapping requires the emitter to move in an arc around them, as is the case in tomography. It is not guaranteed that occluders can be mapped since the emitter's path is not under our control. As long as occlusion is detected, potential errors from summing over nonexistent arrivals are avoided.

When deploying PEEL in an open water setting, there are improvements that must be made to cope with this more difficult environment. To use as general a model as possible, we have used the simplest versions of each algorithm, such as GCC-PHAT for TDOA, basic PSO for optimization, and a constant acceleration model for Kalman tracking. Each of these algorithms can be enhanced if more domain knowledge is available. Applying an algorithm that tracks the multipath arrivals (Fleury , 1999; Salmi , 2008), for example, can help to further refine the boundary estimates. The underwater environment's bathymetric profile or the presence of occluding objects may already be known, and must then be taken into account. PEEL can also be applied to tasks such as indoor (acoustic or RF) localization, which features easier channels, but cluttered environments with different challenges.

PEEL currently incorporates global optimization over the emitter position and the locations of planar boundaries. A more sophisticated environmental model would allow it to work robustly with the complex non-planar boundaries in the real-life shallow-water environment. In the literature, global optimization with simulated annealing and genetic algorithms have been applied to MFP (Hermand, 1999) for geoacoustic inversion, using known emitters. Optimization over not just the locations of their boundaries but their physical properties, while more computationally expensive, could lead to improved results in field deployments.

There are broader questions in TOA localization that pertain to PEEL. Currently, PEEL exploits the emission of deterministic, pulsed signals. In practice, other types of emissions may arise, as with communications sources or vessel signatures. It would be challenging to establish time synchronization with such signals. Although we have experimented with performance in increasing noise, it would also be of interest to obtain a full bandwidth/noise characterization over different environments and signals. Finally, COTANS boundary estimation may be amenable to a machine learning method, to automate the smoothing and peak extraction tasks. COTANS transform results are two-dimensional (2D) or 3D images, so we can train a neural network on such data for boundary estimation in future experiments.

This work was supported, in part, by ONR under Grant Nos. N00014-19-1-2661, N00014-19-1-2662, and N00014-19-1-2665, and NSF under Grant No. CCF-1816209. The authors also thank Dr. Jim Preisig of JP Analytics LLC for valuable discussions and feedback on the work.

1. COTANS transform for tangent lines to ellipses in 2D

We present the derivation of the COTANS transform of tangent lines to an ellipse in 2D. First, we derive the tangent line to a standard ellipse centered around the origin, with a closed-form description in terms of the normal vector to the tangent. We then consider this ellipse to be rotated and translated to its actual position in 2D and modify the normal-vector description of the corresponding rotated and translated tangent. The final parametrization of the tangent in terms of its normal vector is ( ρ COTANS , θ COTANS ); this is a point in { ( ρ , θ ) }-space, where ρ is the length of the vector and θ is its direction as measured with respect to the x-axis. This description is by definition the COTANS transform of the tangent line.

First, we obtain a description of tangent lines to standard ellipses [as in Fig. 7(a)]. The equation for a standard ellipse centered on the origin is ( x 2 / a 2 ) + ( y 2 / b 2 ) = 1, where a and b are its major and minor axes, respectively. Consider a point ( x 0 , y 0 ) on the ellipse, and note that this point can be parametrized as ( a cos α , b sin α ), where the angle α is an intermediate parametrization. Differentiating the ellipse's equation, the slope of a tangent line is
2 x a 2 + 2 y b 2 d y d x = 0 d y d x = b 2 x a 2 y .
(A1)
From the fact that ( x 0 , y 0 ) satisfies both ( x 0 2 / a 2 ) + ( y 0 2 / b 2 ) = 1 and the equation of the tangent's slope, we obtain the equation of the tangent line at ( x 0 , y 0 ) as ( x x 0 / a 2 ) + ( y y 0 / b 2 ) = 1.

Now, we obtain the COTANS transform of this tangent line, by parametrizing it in terms of its Euclidean distance ρ from the origin and the azimuth θ of the normal line to the tangent going through the origin. Thus, a tangent to the ellipse at ( a cos α , b sin α ) will be expressed as a point ( ρ α , θ α ) in the COTANS domain. One derivation of this COTANS transform result is as follows. Our method will obtain ρ ( θ ), the distance to the origin of a given tangent plane as a function of azimuth.

The distance ρ α of the tangent line to the origin is the solution to
min | | x | | 2 s . t . w α T x = 1 ,
(A2)
where x = [ x y ] T is the point of intersection of the tangent line and its normal through the origin, and w α T = [ x 0 / a 2 y 0 / b 2 ] = [ cos α / a sin α / b ]. ρ α can be found by using Lagrange multipliers technique, where
L ( x , λ ) = 1 2 x T x λ ( w α T x 1 ) , x L ( x , λ ) = 0 x λ w α = 0 x = λ w α , λ L ( x , λ ) = 0 w α T x = 1 λ | | w α | | 2 2 = 1 λ = 1 | | w α | | 2 2 , x = λ w α = w α | | w α | | 2 2 = a 2 b 2 a 2 sin 2 α + b 2 cos 2 α [ cos α a sin α b ] T , ρ α = | | x | | 2 = 1 | | w α | | 2 = a b a 2 sin 2 α + b 2 cos 2 α .
(A3)
We now have ρ α in Eq. (A3), but our goal is to eliminate α and obtain ρ ( θ ). From Eq. (A1), we know that lines that are orthogonal to the tangent line have a slope of a 2 y / b 2 x, so
tan θ = a sin α b cos α = a b tan α , cos 2 α = 1 1 + tan 2 α = 1 1 + b 2 a 2 tan 2 θ ,
(A4)
sin 2 α = b 2 a 2 tan 2 θ 1 + b 2 a 2 tan 2 θ .
(A5)
Therefore, for a given θ, the corresponding ρ ( θ ) is obtained after substitution as
ρ ( θ ) = a b a 2 sin 2 α + b 2 cos 2 α = a 2 cos 2 θ + b 2 sin 2 θ .
(A6)

Equation (A6) allows us to obtain the set of points { ( ρ , θ ) } in the COTANS domain for a standard ellipse, but we wish to obtain such points for a general ellipse whose foci are p e and p r. We are given that the LOS distance between the emitter and receiver is d LOS and that the distance of the NLOS reflection from a given boundary is d NLOS. Our method is to first obtain ( ρ , θ ) for a standard ellipse with foci at p r = [ d LOS / 2 0 ] T and p e = [ d LOS / 2 0 ] T, with axes a = d NLOS / 2 , b = d NLOS 2 d LOS 2 / 2. One example of such an ellipse is given in Fig. 8(a). The ellipse corresponding to the true localization scenario is obtained by rotating this standard ellipse as in Fig. 8(b), and then translating it to its correct position as in Fig. 8(c). The ( ρ , θ ) of the tangent line must be correspondingly transformed to match the new position of the rotated and translated tangent that yields the desired ( ρ COTANS , θ COTANS ).

To modify ( ρ , θ ) as per Fig. 8, first consider the azimuth of the vector p e p r as θ rot. To align the standard ellipse with the target ellipse, we replace each ( ρ , θ ) with ( ρ , ( θ + θ rot ) mod 2 π ), which we term ( ρ , θ ). This rotation leaves the ρ-value of each tangent line unchanged, and only affects the azimuth, while the new foci are at p e and p r . Next, we calculate a translation vector p trans = p r p r , which would be added to any point on the rotated standard ellipse to obtain the target ellipse. To obtain the resulting ( ρ COTANS , θ COTANS ) pairs, we first calculate the dot product ρ proj = p trans · ρ ̂, where ρ ̂ = [ cos θ sin θ ] T is the unit vector pointing towards the tangent line. Thus, we project the translation vector p trans onto ρ ̂. If ρ proj 0, then we replace ( ρ , θ ) with ( ρ + ρ proj , θ ). If ρ proj < 0, then we replace ρ with | ρ | ρ proj | |. If | ρ proj | < ρ and ρ proj < 0, then we do not modify the azimuth θ ; else, we replace θ with ( θ + π ) mod 2 π. Carrying out these operations for each of the starting { ( ρ , θ ) } points, we obtain a final transformed set of points as { ( ρ COTANS , θ COTANS ) }.

2. COTANS transform for tangent planes to spheroids in 3D
The equation of an origin-centered spheroid is ( x 2 / a 2 ) + ( y 2 / b 2 ) + ( z 2 / c 2 ) = 1. A point ( x 0 , y 0 , z 0 ) on its surface can be parametrized as ( a cos α sin β , b sin α sin β , c cos β ), where α and β are intermediate parameters. Similarly to 2D, the tangent plane at ( x 0 , y 0 , z 0 ) is ( x x 0 / a 2 ) + ( y y 0 / b 2 ) + ( z z 0 / c 2 ) = 1, and we have an optimization similar to Eq. (A2), this time using w α , β T = [ x 0 / a 2 y 0 / b 2 z 0 / c 2 ] = [ cos α sin β / a sin α sin β / b cos β / c ] instead of w α T. A Lagrangian method as in Eq. (A3) yields
ρ α , β = 1 | | w α , β | | 2 = 1 cos 2 α sin 2 β a 2 + sin 2 α sin 2 β b 2 + cos 2 β c 2 ,
(A7)
x = w α , β | | w α , β | | 2 2 = ρ α , β 2 [ cos α sin β a sin α sin β b cos β c ] T .
(A8)
To eliminate α and β and express x in terms of the azimuth θ and elevation ϕ, we use
x = ρ α , β [ cos θ sin ϕ sin θ sin ϕ cos ϕ ] T .
(A9)
First, noting that tan ϕ = ( a / b ) tan α, we have the same results of Eq. (A4) for cos 2 α and Eq. (A5) for sin 2 α. Note that for the case of a spheroid where a = b, we have α = θ, since the x–y cross section of the ellipsoid is then a circle. To eliminate β, we observe that
tan 2 ϕ = c 2 a 2 cos 2 α cos 2 θ tan 2 β = c 2 a 2 tan 2 β , cos 2 β = 1 1 + a 2 c 2 tan 2 ϕ ,
(A10)
sin 2 β = a 2 c 2 tan 2 ϕ 1 + a 2 c 2 tan 2 ϕ .
(A11)
We substitute into Eq. (A7) our results in Eq. (A4), (A5), (A10), and (A11) to get
ρ ( θ , ϕ ) = c 2 cos 2 ϕ + a 2 sin 2 ϕ .
(A12)
Note that because of symmetry around the z-axis, there is no θ-dependence at this stage.

Using Eq. (A12), we can take the COTANS transform of tangent planes to a standard spheroid with a = b. We can rotate this spheroid and translate it, and correspondingly transform the points ( ρ , θ , ϕ ) that define the tangent planes. Our operations are similar to those in the 2D case. We first define our standard spheroid to have foci at p r = [ 0 0 d LOS / 2 ] T and p e = [ 0 0 d LOS / 2 ] T, with axes a = b = d NLOS 2 d LOS 2 / 2 and c = d NLOS / 2. We calculate the azimuth and elevation of the vector p e p r, and call these angles θ rot and ϕ rot, respectively. We rotate the spheroid by ϕ rot with the y-axis as the axis of rotation (thereby obtaining the correct elevation), then rotate the result by θ rot with the z-axis as the axis of rotation (thereby obtaining the correct azimuth), and finally, we translate it to its correct position.

Rotating the spheroid by ϕ rot will not change the ρ of a given ( ρ , θ , ϕ ) tangent plane, but may affect its azimuth and elevation. The unit vector orthogonal to ( ρ , θ , ϕ ) in Cartesian coordinates is given by ρ ̂ = [ cos θ sin ϕ sin θ sin ϕ cos ϕ ] T. To rotate this vector with the y-axis as the axis of rotation, we multiply it by the corresponding rotation matrix to obtain
[ cos ϕ rot 0 sin ϕ rot 0 1 0 sin ϕ rot 0 cos ϕ rot ] [ cos θ sin ϕ sin θ sin ϕ cos ϕ ] = [ cos ϕ rot cos θ sin ϕ + sin ϕ rot cos ϕ sin θ sin ϕ sin ϕ rot cos θ sin ϕ + cos ϕ rot cos ϕ ] .
(A13)
The azimuth and elevation of the vector obtained in Eq. (A13) are then substituted as the new θ and ϕ values of the plane, θ and ϕ . Next, we rotate the tangent plane by θ rot, which does not change ρ or ϕ but replaces θ with ( θ + θ rot ) mod 2 π = θ . At this stage, the point corresponding to the transformed value of p r is a point p r . Finally, we calculate the translation vector p trans = p r p r , and obtain ρ proj = p trans · ρ ̂. The transformation of the tangent plane in the COTANS domain is similar to the 2D case. If ρ proj 0, we replace ( ρ , θ , ϕ ) with ( ρ + ρ proj , θ , ϕ ). If ρ proj < 0, then we subtract the projection result from ρ, thus replacing ρ with | ρ | ρ proj | |. If | ρ proj | < ρ and ρ proj < 0, then we do not modify the azimuth θ or the elevation ϕ ; else, we replace θ with ( θ + π ) mod 2 π and ϕ with π ϕ . Carrying out this procedure for each of the starting { ( ρ , θ , ϕ ) } points, we obtain a final set of { ( ρ COTANS , θ COTANS , ϕ COTANS ) } points in the 3D COTANS domain.
In the standard PSO algorithm (Clerc, 2010), we define N particles x i ( k ), corresponding to points x i ( k ) = [ x i , 1 ( k ) x i , D ( k ) ] T in a D-dimensional search space at iteration k of the PSO algorithm. This search space is delimited by the boundaries [ x j , min , x j , max ] for each of the j { 1 D } parameters that we optimize over. Each particle has a velocity vector v i ( k ) = [ v i , 1 ( k ) v i , D ( k ) ] T, assigned randomly at initialization. Particles remember their best locations p i in the search space over iterations 1 through k, where their objective function f ( p i ) was the largest. Finally, each particle is in a set of M randomly chosen neighbors, over which we determine the best location p i g over past iterations. The PSO update step for the velocity and position of each particle is then given for each dimension d by
v i d ( k + 1 ) α v i d ( k ) + U ( 0 , β ) ( p i d x i d ( k ) ) + U ( 0 , β ) ( p i g x i d ( k ) ) ,
(B1)
x i d ( k + 1 ) x i d ( k ) + v i d ( k + 1 ) ,
(B2)
where α and β that define a particle's inertia and acceleration for updating its velocity (standard values are 0.7298 and 1.4986, respectively), and U ( 0 , β ) is a sample of the uniform random variable. Particles that exit the boundaries of the search space are reflected back.
1

In practice, if Tp > Tc, then Tp can be estimated, but this is an environment-specific problem outside the scope of our method, as is an estimation of Tm.

2

In our scale watertank setting, the regularly-spaced grid-squares have length 1 mm.

3

More sophisticated methods for extracting the LOS peak could also be used (Guvenc and Sahinoglu, 2005).

4

Algorithm 2 does not require leveraging all of the boundaries to yield an accurate result. The additional NLOS arrivals provide spatial diversity that helps to increase localization accuracy and mitigate occlusion.

5

We use a resolution of 0.1 mm for ρ and 0.1° for θ and ϕ in our experiments.

6

We heuristically use a uniform averaging filter that is 3 mm wide in ρ, and 3° wide in θ and ϕ.

7

Given a maximum of the accumulator result at some ( ρ max , θ max , ϕ max ), we zero out the θ and ϕ within ± 15 ° of the maximum, and all ρ values within this sector.

8

70% below the predicted magnitude, in our implementation.

1.
Alameda-Pineda
,
X.
, and
Horaud
,
R.
(
2014
). “
A geometric approach to sound source localization from time-delay estimates
,”
IEEE/ACM Trans. Audio. Speech. Lang. Process.
22
(
6
),
1082
1095
.
2.
Baggeroer
,
A. B.
,
Kuperman
,
W. A.
, and
Mikhalevsky
,
P. N.
(
1993
). “
An overview of matched field methods in ocean acoustics
,”
IEEE J. Oceanic Eng.
18
(
4
),
401
424
.
3.
Borrmann
,
D.
,
Elseberg
,
J.
,
Lingemann
,
K.
, and
Nüchter
,
A.
(
2011
). “
The 3D Hough transform for plane detection in point clouds: A review and a new accumulator design
,”
3D Res.
2
(
2
),
3
.
4.
Brutti
,
A.
,
Omologo
,
M.
, and
Svaizer
,
P.
(
2010
). “
Multiple source localization based on acoustic map de-emphasis
,”
EURASIP J. Audio Speech Music Process.
2010
,
1
17
.
5.
Chitre
,
M.
(
2007
). “
A high-frequency warm shallow water acoustic communications channel model and measurements
,”
J. Acoust. Soc. Am.
122
(
5
),
2580
2586
.
6.
Chitre
,
M.
, and
Pelekanakis
,
K.
(
2014
). “
Channel variability measurements in an underwater acoustic network
,” in
Proceedings of UComms 2014
,
September 3–5
,
Levante, Italy
, pp.
1
4
.
7.
Clerc
,
M.
(
2010
).
Particle Swarm Optimization
(
John Wiley & Sons
,
New York
).
8.
Deane
,
G. B.
(
1994
). “
A three-dimensional analysis of sound propagation in facetted geometries
,”
J. Acoust. Soc. Am.
96
(
5
),
2897
2907
.
9.
Demirli
,
R.
, and
Saniie
,
J.
(
2001
). “
Model-based estimation of ultrasonic echoes. Part I: Analysis and algorithms
,”
IEEE Trans. Ultrason, Ferroelect, Freq. Contr.
48
(
3
),
787
802
.
10.
Dokmanic
,
I.
,
Parhizkar
,
R.
,
Ranieri
,
J.
, and
Vetterli
,
M.
(
2015
). “
Euclidean distance matrices: Essential theory, algorithms, and applications
,”
IEEE Signal Process. Mag.
32
(
6
),
12
30
.
11.
Fischell
,
E. M.
,
Kroo
,
A. R.
, and
O'Neill
,
B. W.
(
2019
). “
Single-hydrophone low-cost underwater vehicle swarming
,”
IEEE Robot. Autom. Lett.
5
(
2
),
354
361
.
12.
Fleury
,
B. H.
,
Tschudin
,
M.
,
Heddergott
,
R.
,
Dahlhaus
,
D.
, and
I. Pedersen
,
K.
(
1999
). “
Channel parameter estimation in mobile radio environments using the SAGE algorithm
,”
IEEE J. Select. Areas Commun.
17
(
3
),
434
450
.
13.
Guvenc
,
I.
, and
Sahinoglu
,
Z.
(
2005
). “
Threshold selection for UWB TOA estimation based on kurtosis analysis
,”
IEEE Commun. Lett.
9
(
12
),
1025
1027
.
14.
Hermand
,
J.-P.
(
1999
). “
Broad-band geoacoustic inversion in shallow water from waveguide impulse response measurements on a single hydrophone: Theory and experimental results
,”
IEEE J. Oceanic Eng.
24
(
1
),
41
66
.
15.
Kim
,
J.
,
Kim
,
J.
,
Nguyen
,
L. T.
,
Shim
,
B.
, and
Hong
,
W.
(
2019
). “
Tonal signal detection in passive sonar systems using atomic norm minimization
,”
EURASIP J. Adv. Signal Process.
2019
,
43
.
16.
Knapp
,
C.
, and
Carter
,
G.
(
1976
). “
The generalized correlation method for estimation of time delay
,”
IEEE Trans. Acoust, Speech, Signal Process.
24
(
4
),
320
327
.
17.
Korhonen
,
T.
(
2008
). “
Acoustic localization using reverberation with virtual microphones
,” in
Proceedings of the IWAENC 2008
,
September 14–17
,
Seattle, WA
, pp.
211
223
.
18.
Le Gall
,
Y.
,
Socheleau
,
F.
, and
Bonnel
,
J.
(
2016
). “
Matched-field performance prediction with model mismatch
,”
IEEE Signal Process. Lett.
23
(
4
),
409
413
.
19.
Li
,
L.
, and
Krolik
,
J.
(
2013
). “
Simultaneous target and multipath positioning
,”
IEEE J. Sel. Top. Signal Process.
8
(
1
),
153
165
.
20.
Naseri
,
H.
,
Costa
,
M.
, and
Koivunen
,
V.
(
2014
). “
Multipath-aided cooperative network localization using convex optimization
,” in
Proceedings of the 48th Asilomar Conference on Signals, Systems and Computers
,
November 2–5
,
Pacific Grove, CA
, pp.
1515
1520
.
21.
Naseri
,
H.
, and
Koivunen
,
V.
(
2016
). “
Cooperative simultaneous localization and mapping by exploiting multipath propagation
,”
IEEE Signal Process. Mag
65
(
1
),
200
211
.
22.
Qiao
,
T.
,
Redfield
,
S.
,
Abbasi
,
A.
,
Su
,
Z.
, and
Liu
,
H.
(
2013
). “
Robust coarse position estimation for TDOA localization
,”
IEEE Wireless Commun. Lett.
2
(
6
),
623
626
.
23.
Ribeiro
,
F.
,
Zhang
,
C.
,
Florêncio
,
D. A.
, and
Ba
,
D. E.
(
2010
). “
Using reverberation to improve range and elevation discrimination for small array sound source localization
,”
IEEE Trans. Audio. Speech. Lang. Process.
18
(
7
),
1781
1792
.
24.
Salmi
,
J.
,
Richter
,
A.
, and
Koivunen
,
V.
(
2008
). “
Detection and tracking of MIMO propagation path parameters using state-space approach
,”
IEEE Trans. Signal Process.
57
(
4
),
1538
1550
.
25.
Scheuing
,
J.
, and
Yang
,
B.
(
2008
). “
Disambiguation of TDOA estimation for multiple sources in reverberant environments
,”
IEEE Trans. Audio. Speech. Lang. Process.
16
(
8
),
1479
1489
.
26.
Setlur
,
P.
,
Smith
,
G. E.
,
Ahmad
,
F.
, and
Amin
,
M. G.
(
2012
). “
Target localization with a single sensor via multipath exploitation
,”
IEEE Trans. Aerosp. Electron. Syst.
48
(
3
),
1996
2014
.
27.
Srinivasan
,
P.
,
Liang
,
P.
, and
Hackwood
,
S.
(
1990
). “
Computational geometric methods in volumetric intersection for 3D reconstruction
,”
Pattern Recogn.
23
(
8
),
843
857
.
28.
Su
,
Z.
,
Shao
,
G.
, and
Liu
,
H.
(
2017
). “
Semidefinite programming for NLOS error mitigation in TDOA localization
,”
IEEE Commun. Lett.
22
(
7
),
1430
1433
.
29.
Sun
,
W.
, and
Wang
,
Z.
(
2018
). “
Online modeling and prediction of the large-scale temporal variation in underwater acoustic communication channels
,”
IEEE Access
6
,
73984
74002
.
30.
Too
,
Y. M.
,
Chitre
,
M.
,
Barbastathis
,
G.
, and
Pallayil
,
V.
(
2017
). “
Localizing snapping shrimp noise using a small-aperture array
,”
IEEE J. Oceanic Eng.
44
(
1
),
207
219
.
31.
Venkateswaran
,
S.
, and
Madhow
,
U.
(
2012
). “
Localizing multiple events using times of arrival: A parallelized, hierarchical approach to the association problem
,”
IEEE Trans. Signal Process.
60
(
10
),
5464
5477
.
32.
Waterston
,
J.
,
Rhea
,
J.
,
Peterson
,
S.
,
Bolick
,
L.
,
Ayers
,
J.
, and
Ellen
,
J.
(
2019
). “
Ocean of Things: Affordable maritime sensors with scalable analysis
,” in
OCEANS 2019-Marseille
, pp.
1
6
.
33.
Xu
,
W.
,
Baggeroer
,
A. B.
, and
Schmidt
,
H.
(
2006
). “
Performance analysis for matched-field source localization: Simulations and experimental results
,”
IEEE J. Oceanic Eng.
31
(
2
),
325
344
.
Published open access through an agreement with Massachusetts Institute of Technology