# Target motion analysis using angular measurements and underwater sound pressure levels

Enzo Iglésis, Luc de Montella

► To cite this version:

Enzo Iglésis, Luc de Montella. Target motion analysis using angular measurements and underwater sound pressure levels. 2025. hal-04796331v2

HAL Id: hal-04796331

https://inria.hal.science/hal-04796331v2

Preprint submitted on 25 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Target motion analysis using angular measurements and underwater sound pressure levels

Enzo Iglésis¹ and Luc de Montella*¹,²

¹Inria centre at the University of Bordeaux, Talence, 33405, France

²Naval Group, Bouguenais, 44340, France.

April 25, 2025

## Abstract

This paper addresses the problem of target motion analysis using data acquired by sonars operating in passive mode. Most theoretical and applied studies rely solely on bearing information, and various filtering algorithms, especially particle filter methods, have demonstrated strong performance in this context. To improve upon existing approaches, bearing information is combined with received underwater sound pressure levels. This method requires no additional sensors compared to bearings-only trackings, as hydrophones — already used for beamforming — can also measure sound pressure levels. Additionally, underwater transmission losses can be estimated through simulations based on commonly available environmental data. The effectiveness of the proposed method is demonstrated using simulated data, showing improved performance — especially in scenarios without observer maneuvers. To increase realism, the impact of replacing bearing measurements with conical angle measurements is also investigated. These measurements account for both elevation and bearing angles, thereby providing a fairly comprehensive characterization of bottom-bounce paths commonly encountered under real-world conditions. This enables exploration of scenarios with different immersion levels for both target and observer. The results further emphasize the benefits of incorporating sound pressure levels measurements into the tracking process.

## 1 Introduction

Target motion analysis (TMA) is a set of techniques used to estimate the position of a target based on data collected by passive sensors, such as passive sonars. These techniques analyze the direction of incoming sound over time and compare the movement of the target to that of the observer (Arulampalam et al., 2004; Farina, 1999). While passive sonars primarily provide directional information, they can also collect additional data such as frequency information (Maranda, 2008). Thus, in some cases, bearing measurements are fused with this data to improve tracking accuracy by taking advantage of the Doppler shift on radiated frequencies (Passerieux et al., 1989; Laneuville, 2007). Nevertheless, TMA remains challenging due to the inherent nonlinearities of the problem and the uncertainties associated with acoustic signal propagation in underwater environments.

Various filtering techniques have been developed to address the complexities of TMA, with particle filters (PFs) emerging as one of the most effective approaches (Arulampalam et al., 2004; Ristic et al.,

2019). PFs are well-suited to this problem due to their ability to cope with nonlinear estimation problems. However, despite their effectiveness, traditional bearings-only tracking (BOT) methods have limitations, particularly in scenarios where observer maneuvers are restricted due to the noise they generate — which compromises stealthiness (Ovalle and García-Pelaez, 2015) — and negatively affect the processing gain of passive sonar. In such conditions, bearing measurements alone often lack the necessary information to accurately resolve the target's motion (Nardone and Aidala, 1981).

To overcome these limitations, this paper introduces a novel approach that combines bearing information with received sound pressure level (SPL) measurements. Since sonar systems are essentially transducer arrays, they inherently capture these intensity levels (Urick, 1983, Section 3.4). By incorporating SPL data into the tracking process, it is possible to improve TMA accuracy without the need for additional sensors.

In order to fuse SPL information with bearing measurements, an estimation model for the underwater transmission losses (TLs) is required. This is achieved through computational simulations based on environ-

<footnote>*Luc de Montella would like to thank Naval Group for their financial support through the CIFRE PhD program.</footnote>
mental data such as the sound speed profile, the bathymetric profile and the frequency of the target's signal.

Since sound does not travel in a straight line underwater due to variations in sound speed with depth, different receiver positions can result in identical TL values, thereby introducing ambiguity. This multimodal behavior, combined with the nonlinear nature of the problem, leads to measurements that are both highly nonlinear and ambiguous. These challenges make PFs an ideal choice for addressing such issues. Moreover, this problem is analogous to terrain-aided navigation, where specialized PFs have been developed to handle similar complexities (Palmier et al., 2019; Merlinge et al., 2019), many of which could potentially be used to further improve the results presented in this paper.

However, the purpose of this study is not to determine the best filter for this type of problem — an area that has been extensively researched — but rather to demonstrate that incorporating previously unused information can improve tracking accuracy. To this end, the filter used in this paper is the regularized particle filter (RPF) (Musso et al., 2001), which has been applied to both BOT and TMA with bearing and SPL measurements, providing a fair basis for comparison.

The effectiveness of the approach is demonstrated through simulations, which show improved tracking performance compared to both traditional BOT methods and Doppler-bearing tracking. In particular, the numerical results highlight how the proposed method allows for a more accurate estimation of the target's position even without observer maneuvers, underscoring its potential for real-world applications.

Assuming the observer measures bearing angles is a widely used approximation for scenarios where the observer and target are at the same depth and the observer receives sound in a direct path. In addition, tracking based on bearing and SPL measurements highlights the ability of SPL to provide range information that bearing measurements alone cannot. However, as discussed in Oh et al. (2020), a more realistic approach is to replace the bearing information with a conical angle measurement and to account for the bottom bounce paths of the sound. Not only is this conical angle measurement more consistent with real-world scenarios, but it also provides richer data by deriving the measurement from a formula that includes both elevation and bearing information. However, depending on the sensor used, it may lead to left/right ambiguity (e.g., towed array). Additionally, uncertainties in the ray path (e.g., surface bounce, direct path, multiple bounces, ...) can introduce further difficulties in precisely estimating the target state. This paper is based on specific considerations — outlined in subsequent sections — that allow these difficulties to be overcome. A numerical study is then conducted, demonstrating how the inclusion of SPL, despite the use of conical angle measurements, still results in superior performance.

The remainder of this paper is organized as follows. Section 2 formally defines the TMA problem and provides a detailed description of BOT, conical angle measurements and frequency measurements. In Section 3, the acoustical emission model of the target and the computation of its transmission loss are detailed to provide a measurement function. Section 4 introduces the RPF that is used for both the BOT and TMA with combined bearing and SPL. The simulation setups and results are presented in Section 5, where the effectiveness of the approach is compared with traditional BOT methods and methods based on conical angle measurements. Finally, Section 6 concludes the paper by summarizing key findings and suggesting directions for future research.

# 2 Description of the Frameworks

## 2.1 Bearing-only tracking

Throughout this paper, it is assumed that the depths of both the observer and the target are known and constant. To consider an observer measuring bearing angles, it is assumed that both the ownship and target depths are equal, which allows for approximating the direct path of the acoustic signal using a simpler model. Therefore, the focus is on the two-dimensional version of the BOT problem, as presented in Arulampalam et al. (2004, Chapter 6).

A Cartesian coordinate system is considered, where the positions of the observer and the target at time $k$ are denoted, respectively, by $[x_k^t \ y_k^t]^\top$ and $[x_k^o \ y_k^o]^\top$. The velocity components are denoted, for the observer and the target at time $k$ respectively by $[\dot{x}_k^t \ \dot{y}_k^t]^\top$ and $[\dot{x}_k^o \ \dot{y}_k^o]^\top$. To simplify notation, the state vector $\mathbf{x}_k^t$ is defined as:

$$ \mathbf{x}_k^t = [x_k^t \ y_k^t \ \dot{x}_k^t \ \dot{y}_k^t]^\top \quad (1) $$

and the state vector $\mathbf{x}_k^o$ as:

$$ \mathbf{x}_k^o = [x_k^o \ y_k^o \ \dot{x}_k^o \ \dot{y}_k^o]^\top. \quad (2) $$

Let $\mathbf{x}_k$ be the vector representing the position of the target relative to the observer at time $k$, i.e.:

$$ \mathbf{x}_k = \mathbf{x}_k^t - \mathbf{x}_k^o := [x_k \ y_k \ \dot{x}_k \ \dot{y}_k]^\top. \quad (3) $$

The target is assumed to have a uniform linear motion. Denoting $\Delta T$ as the sampling interval, and $\mathbf{F}$ as the state transition matrix, the state equation can then be written as:

$$ \mathbf{x}_{k+1} = \mathbf{F}(\mathbf{x}_k + \mathbf{x}_k^o) - \mathbf{x}_{k+1}^o + \mathbf{\Gamma}\nu_k, \quad (4) $$

where $\mathbf{F}$ is transition matrix given by

$$ \mathbf{F} = \begin{bmatrix} 1 & 0 & \Delta T & 0 \\ 0 & 1 & 0 & \Delta T \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix}, \quad (5) $$
and $\Gamma$ is the process noise distribution matrix given by

$$ \Gamma = \begin{bmatrix} \frac{\Delta T^2}{2} & 0 \\ 0 & \frac{\Delta T^2}{2} \\ \Delta T & 0 \\ 0 & \Delta T \end{bmatrix}. \qquad (6) $$

The vector $\nu_k \sim \mathcal{N}(0, Q)$ is a white noise with covariance matrix $Q = \sigma_a I_2$, where $\sigma_a$ is the standard deviation of the acceleration process noise and $I_2$ stands for the $2 \times 2$ identity matrix.

At each time step, the observer obtains a new measurement of the angle between the relative vector $[x_k \ y_k]^\top$ and the y-axis. The measurement vector, denoted by $\mathbf{z}_k$, is given in the case of bearings-only tracking by

$$ \mathbf{z}_k = [\theta_k], \qquad (7) $$

where $\theta_k$ denoted the bearing measurement at time step $k$. It is given by:

$$ \theta_k = h_{\theta}(\mathbf{x}_k) + \omega_k, \qquad (8) $$

where $(\omega_k)_{k \in \mathbb{N}}$ is a family of zero-mean independent Gaussian noise with standard deviation $\sigma_\theta$. The function $h_\theta(\cdot)$ represents the true bearing and is defined as:

$$ h_{\theta}(\mathbf{x}_k) = \operatorname{atan2}(x_k, y_k), \qquad (9) $$

where atan2($x_k, y_k$) denotes the measure of the angle between the north axis (positive y-axis) and the observer-target vector.

The BOT problem is the estimation of the vectors $\mathbf{x}_k$ from the set $\mathbf{Z}_k = \{\mathbf{z}_i\}_{i \in [1, k]}$ of measurements up to time $k$.

## 2.2 Conical angle measurements

As mentioned previously, the measurement equation (9) represents a simplified model in which an observer positioned at the same depth as the target receives acoustic signals propagating in a straight line along a direct path. In reality, however, sound does not propagate in a straight line, and in some cases, nonlinearity prevents the existence of a direct path altogether. This makes it necessary to account for bottom bounce or surface bounce propagation, especially since bottom-bounce paths, which result from the reflection of acoustic waves off the ocean floor, can introduce bearing inaccuracies that affect BOT performance.

In this paper, the bathymetry is assumed to be flat with a depth of 2000 m, and the sound speed profile $C(z)$ follows a Munk's sound speed profile, representing the idealized sound channel described in Munk (1974):

$$ C(z) = C_0 \left[ 1 + \epsilon \left( \eta(z) + e^{-\eta(z)} - 1 \right) \right], \qquad (10) $$

where $z$ denote the immersion, $C_0$ the sound velocity at the channel axis $z_c$, $\epsilon$ the perturbation coefficient

and $\eta(z)$ a dimensionless distance between channel axis given by:

$$ \eta(z) = \frac{2(z - z_c)}{z_c}. \qquad (11) $$

From (10) and with $C_0 = 1500 \text{ m s}^{-1}$, $z_c = 400 \text{ m}$ and $\eta = 7.4 \times 10^{-3}$ the obtained sound speed profile is illustrated in FIG. 1

<figure><img src="image_2.png" /><figcaption>Figure 1: Sound speed profile obtained from (10) with $C_0 = 1500 \text{ m s}^{-1}$, $z_c = 400 \text{ m}$ and $\eta = 7.4 \times 10^{-3}$.</figcaption></figure>

This non-uniform sound speed profile causes sound to propagate in a nonlinear manner, curving its path. The sound paths generated using a ray tracing method based on Munk's sound speed profile from FIG. 1 are illustrated in FIG. 2.

<figure><img src="image_3.png" /><figcaption>Figure 2: Ray paths generated using bellhopcxx based on Munk's sound speed profile from FIG. 1 with a source at 10 m and a frequency of 60 Hz.</figcaption></figure>

In this paper, the software used to determine the sound path is bellhopcxx from Pisha et al. (2023), an open-source code based on BELLHOP (Porter and Bucker, 1987), a widely used beam-tracing model. By specifying environmental data, such as bathymetry and sound speed profiles, along with other parameters, the software computes the eigenrays at specified receiver positions. BELLHOP is capable of producing eigenrays in a 3D environment that accounts for variations in the sound speed profile. However, for the sake of brevity in Section 5 and the Appendix,
where detailed values are provided, and for computational efficiency, the acoustic environment in this paper is simplified to a flat-bottom 2D model. This 2D environment is then interpreted as a flat-bottom 3D environment by assuming that sound propagates uniformly in all directions. The input environmental file for BELLHOP, used to obtain FIG. 2 (except for the run-type set to 'R'), is provided in the Appendix.

FIG. 2 illustrates that when an observer operates at depths between 0 m and 500 m, and at ranges between 3 km and 10 km, no direct acoustic paths are received. This depth and range span represents the primary operational domain considered in this study, underscoring the importance of accounting for conical angles.

Note that acoustic signals can also bounce off the surface, and that FIG. 2 shows the ray paths as if the surface were absorbing 100% of the signal, which is not representative of real-world conditions. However, for clarity, the rays have been truncated at the surface because only direct and bottom-bounce rays are considered here, as these are expected to be the most energetic signals received by the observer in this scenario.

Having said that, a more realistic measurement model than (9) is needed. To this end, the model that best describes the measurements from the passive sonar is the conical angle measurement model, given by Oh et al. (2020):

$$\tilde{h}_{\theta}(\mathbf{x}_{k}^{t}, \mathbf{x}_{k}^{o}) = \arccos(\cos(\varphi_{b}(\mathbf{x}_{k}^{t}, \mathbf{x}_{k}^{o})) \cos(\mu(\mathbf{x}_{k}^{t}, \mathbf{x}_{k}^{o}))), \quad (12)$$

In this equation, $\mu(\cdot, \cdot)$ stands for the elevation angle in the vertical plane, and $\varphi_b(\cdot, \cdot)$ is obtained by subtracting the direction of the observer to the bearing angle:

$$\varphi_b(\mathbf{x}_k^t, \mathbf{x}_k^o) = h_\theta(\mathbf{x}_k) - \psi(\mathbf{x}_k^o), \quad (13)$$

with $\psi(\mathbf{x}_k^o)$ the heading of the ownship given by:

$$\psi(\mathbf{x}_k^o) = \operatorname{atan2}(\dot{x}_k^o, \dot{y}_k^o). \quad (14)$$

The $\mu(\cdot, \cdot)$ function is obtained using ray tracing software, such as the previously introduced bellhopccxx, with the run-type set to 'E' and retaining only the closest bottom-bounce ray.

The angular measurement $\theta_k$ in this situation is obtained by simply replacing $h_\theta$ with $\tilde{h}_\theta$ in (8).

All these angles are illustrated in FIG. 3.

Note that $\tilde{h}_\theta(\mathbf{x}_k^t, \mathbf{x}_k^o)$ can only take values from 0° to 180° due to (12), as the conical angle only indicates the magnitude of the measured angle relative to the heading direction of the line array. This uncertainty in determining the exact direction of the arriving signal is known as left/right ambiguity.

However, in real-life conditions, this ambiguity applies only to sonar arrays, such as towed arrays. For simplicity, and based on the scenario considered in this paper, it is assumed that the submarine is

<figure><img src="image_4.png" /><figcaption>Figure 3: Geometry between the observer and the target, with the angle used in the conical angle measurement equation (12). The scalars $z_k^o$ and $z_k^t$ denote the depths of the ownship and the target at time step $k$, respectively.</figcaption></figure>

equipped with two flank arrays, each having a 180° field of view. In this configuration, the starboard array cannot receive signals from the port side, and vice versa. As a result, left/right ambiguity does not exist in this scenario and is not considered in this paper. Therefore it is assumed that $\tilde{h}_\theta(\mathbf{x}_k^t, \mathbf{x}_k^o)$ can take any value from -180° to 180°.

## 2.3 Frequency measurements

Along with bearing measurements, frequency data is sometimes used to improve the accuracy of TMA. Taking these measurements into account extends the range of scenarios in which the target's position can be observed, such as when the observer remains stationary Passerieux et al. (1988); Jauffret and Pérez (2018). To provide a reference for comparison in the numerical simulations, the framework associated with these measurements is detailed hereafter.

As both ships are moving, the frequency emitted by the target undergoes a Doppler shift — except in the unlikely case that they are traveling at exactly the same speed and in exactly the same direction. At time $k$, assuming the target emits a pure single tone at an unknown frequency $\gamma_0$, the observer measures both the bearing $\theta_k$ and the frequency.

$$f_k = h_f(\mathbf{x}_k) + u_k, \quad (15)$$

where $(u_k)_{k \in \mathbb{N}}$ is a family of zero-mean independent Gaussian noise with standard deviation $\sigma_f$.

The function $h_f(\cdot)$ represents the true shifted frequency and is defined as:

$$h_f(\mathbf{x}_k) = \gamma_0 \left( 1 - \frac{\dot{x}_k \sin(h_\theta(\mathbf{x}_k)) + \dot{y}_k \cos(h_\theta(\mathbf{x}_k))}{C_0} \right). \quad (16)$$

The measurement vector at time $k$ becomes then

$$\mathbf{z}_k = [\theta_k \ f_k]. \quad (17)$$
The true emitted frequency is assumed to be constant over time. However, since it is unknown, filters using frequency data must estimate $\gamma_0$ along with the target's position and speed. Therefore, in these scenarios, an additional dimension in the state vector $\mathbf{x}_k$ is considered to account for the emitted frequency:

$$ \mathbf{x}_k = \mathbf{x}_k^t - \mathbf{x}_k^o := [x_k \quad y_k \quad \dot{x}_k \quad \dot{y}_k \quad \gamma_k]^\top. \quad (18) $$

The state equation remains unchanged for the position and speed, however it is extended with the following equation:

$$ \gamma_k = \gamma_{k-1} + s_k, \quad (19) $$

where $(s_k)_{k \in \mathbb{N}}$ is a family of zero-mean independent Gaussian noise with standard deviation $\sigma_s$.

## 3 Source level and transmission loss

As the name suggests, in bearings-only tracking, the observer relies solely on angular measurements to estimate the target's position, which presents a significant challenge due to the absence of direct range information. Despite this, traditional BOT disregards the received sound pressure level measurements, a potentially critical omission. While the exact distance between the observer and the target remains unknown, knowing the target's source level (SL) and the signal's transmission loss (TL) can significantly narrow down the target's potential positions. Indeed, the passive sonar equation in terms of intensity is given by (Maranda, 2008, Section 3):

$$ \mathrm{SE} = \mathrm{SL} - \mathrm{NL} + \mathrm{AG} - \mathrm{TL} - \mathrm{DT}, \quad (20) $$

where SE, NL, AG, and DT respectively denote the signal excess at the input to the detector, the noise level at the receiver, the array gain of the sensor system relative to an omnidirectional hydrophone, and the detection threshold.

In this paper, the sonar is assumed to be ideal, with no array gain (the signal is received with the same intensity regardless of directionality) and no detection threshold (the signal is always received and measured). Therefore, the signal excess corresponds to the received SPL, which is given by:

$$ \mathrm{SPL} = \mathrm{SL} + \mathrm{NL} - \mathrm{TL}. \quad (21) $$

To take full advantage of this insight, an accurate model representing both source level and transmission loss is essential. A simplified noise level model is used in this paper and is described later.

### 3.1 Target source level

Assuming that the target's noise propagates isotropically, the model used to estimate the source level is

that of Ross (1987, Section 8.6). This model, based on statistical observations of surface ship radiated levels, suggests a linear relationship between the source level and the logarithm of the ship's speed. Using a reference speed of $\approx 5.14 \, \mathrm{m s}^{-1}$ (10 kn), the resulting trend curve follows this established relationship:

$$ \mathrm{SL} (\mathbf{x}_k^t) = 170 + 53 \log \left( \frac{v_k^t}{5.14 \, \mathrm{m s}^{-1}} \right), \quad (22) $$

with $v_k^t \propto \sqrt{x_k^t + y_k^t}$.

Despite its apparent simplicity, this model has proven to be quite effective, and even recent, more sophisticated models retain this general form (see Wittekind, 2014; Jalkanen et al., 2018). The key difference is that these advanced models incorporate additional tuning parameters, allowing for better adaptation to the specific characteristics of the target.

It is also important to note that filtering methods are used here, rather than inverse methods, making minor inaccuracies in the model less critical. In fact, any errors in the estimation can be viewed as uncertainties within the model, which the filtering process can effectively handle Elfring et al. (2021).

### 3.2 Transmission loss

The transmission loss can be estimated numerically. As for the eigenrays described in Section 2.2, bellhopcxx is also used here to estimate the TL. BELLHOP provides three TL modes: coherent, semi-coherent, and incoherent. The incoherent mode is used in this paper, due to the instability of the detailed fluctuations in the coherent field. The assumption of a flat seabed also allows for the computation of 2D data, which can then be extrapolated to a 3D environment by assuming that sound propagates uniformly in all directions. Therefore, only the range and depth of the receiver are required to compute the transmission loss. The TL function is then denoted

$$ \mathrm{TL} (\mathbf{x}_k^t, \mathbf{x}_k^o). \quad (23) $$

Since bellhopcxx produce a TL array at finite position in both range and depth, in this paper, for sake of simplicity and computational efficiency a linear interpolation is performed to estimate values between each point.

Based on the sound-speed profile described in Section 2.2 and the provided environmental data file in the Appendix, the resulting transmission loss is illustrated in FIG. 4.

### 3.3 Measurement equation

At each time step, in addition to the angle measurement $\theta_k$ and the possible frequency measurement $f_k$, the observer receives another observation denoted $l_k$. The measurement vector $\mathbf{z}_k$ is then given by:

$$ \mathbf{z}_k = [\theta_k \quad l_k]^\top \quad \text{or} \quad \mathbf{z}_k = [\theta_k \quad f_k \quad l_k]^\top. \quad (24) $$
<figure><img src="image_5.png" /><figcaption>Figure 4: Side view of transmission loss with the environmental data of Appendix using incoherent Gaussian beams.</figcaption></figure>

The $l_k$ observation represents the sound pressure level received by the observer, which provides valuable information about the target’s emission. The measurement equation of the received SPL is given by:

$$l_k = h_l(\mathbf{x}_k^t, \mathbf{x}_k^o) + v_k, \quad (25)$$

where the function $h_l(\cdot)$ represents the true SPL received. Using (21) it is defined as:

$$h_l(\mathbf{x}_k^t, \mathbf{x}_k^o) = \mathrm{SL}(\mathbf{x}_k^t) - \mathrm{TL}(\mathbf{x}_k^t, \mathbf{x}_k^o). \quad (26)$$

The scalar $v_k$ represents the ambient noise level — previously denoted NL in (20) and (21) — which is modeled in this paper as Gaussian noise with zero mean and a standard deviation $\sigma_l$ (for more details on ambient ocean noise, see for example Curtis et al. (1999) or Urick (1983, Chapter 7)). Note that a more realistic noise could be handled here, since the PFs are not limited to Gaussian noise, unlike Kalman-based approaches. Unknown measurement noise could even be accounted for using application-specific PFs such as in Palmier et al. (2019).

This additional observation $l_k$ complements the bearing angle by providing indirect information about the distance between the target and the observer, which helps to reduce the uncertainty in the estimated position of the target.

# 4 Regularized particle filter

Due to its nonlinear nature, the TMA problem has no known exact filter solution. To address this challenge, several modifications of the Kalman filter have been proposed. Among these, the extended Kalman filter (EKF) (Anderson and Moore, 2012) and the unscented Kalman filter (UKF) (Julier et al., 2000) are commonly used. Both filters attempt to linearize the nonlinearities inherent to the BOT problem. However, they have limitations in terms of accuracy and robustness.

The most effective filters for this problem, tend to be variations of the particle filter (Sindhu et al., 2019; Arulampalam et al., 2004). Their ability to handle nonlinearities and non-Gaussian noise makes them particularly well-suited for the BOT problem. In addition, unlike most Kalman-based methods, PFs do not require the computation of the Jacobian of the observation function. This is particularly advantageous, as the Jacobians associated with the observation functions of both SPL and conical angle measurements—being based on numerical simulations—are unknown.

PFs are a class of sequential Monte Carlo methods that use a large number of discrete points, known as particles, to represent the posterior density $p(\mathbf{x}_k | \mathbf{z}_k)$. The total number of particles used is denoted $N$ in this paper.

The PF estimates the posterior density using two primary steps: the prediction step (a. k. a. mutation) and the update step (a. k. a. selection). The behavior of these steps, along with potential additional steps, depends on the version of the PF being used.

In this paper, the PF chosen is the RPF (Musso et al., 2001), which, in addition to the two primary steps mentioned above, includes a regularization-resampling step. This filter was selected because it is widespread and efficient without being application-specific.

It is important to note that some algorithms could be better suited to the problem at hand. For example, many algorithms developed for terrain-aided navigation have demonstrated improved performance in comparable scenarios. These algorithms, designed to match environmental data with real-time sensor inputs, could also prove highly effective for the application discussed in this paper (see Ma et al., 2023; Palmier et al., 2019; Merlinge et al., 2019).

However, evaluating and determining the optimal algorithm for this specific problem is beyond the scope of this paper. For this reason, the RPF was selected for the simulations as it offers a balanced approach between complexity and effectiveness. Its main steps are briefly presented.

## 4.1 Initialization

Using an a priori distribution, $N$ initial potential positions $(\mathbf{x}_0^{(i)})_{i \in [1, N]}$ of the target are drawn. To each of these positions, a weight $w_i$ is assigned, initially set to $\frac{1}{N}$.

## 4.2 Prediction step

The goal of this step is to predict the position of the target at the next measurement time. This is done by propagating the particles using the following probability transition density:
$$ \mathbf{x}_{k|k-1}^{(i)} \sim p(\mathbf{x}_{k|k-1} | \mathbf{x}_{k-1}^{(i)}), \quad (27) $$

This probability transition density is determined by the dynamics of the system. Thus, at this stage, each particle independently evolves according to the state equation (4).

## 4.3 Update step

The goal of this stage is to incorporate a new observation to improve the target's position estimation. This is achieved by adjusting the particles' weights according to the likelihood of the new observation. The weight of the i-th particle is given by

$$ \tilde{w}_{k}^{(i)} \propto w_{k-1}^{(i)} p(\mathbf{z}_{k} | \mathbf{x}_{k|k-1}^{(i)}). \quad (28) $$

## 4.4 Regularization-resampling step

If a few particles dominate the distribution, a resampling process is triggered to remove particles with small weights and duplicate those with large weights. This process occurs when the efficiency (Kong et al., 1994) falls below a user-defined threshold, denoted $N_{thr}$. During resampling, a new set of particles is drawn from the existing set, with probabilities proportional to their weights, and each new particle is assigned a weight of $1/N$.

To maintain diversity within the particle set, the RPF applies a regularization after each resampling. This step involves adding independent noise to each particle. The regularization equation is:

$$ \mathbf{x}_{k}^{(i)} = \hat{\mathbf{x}}_{k}^{(i)} + h\mathbf{D}_{k}\boldsymbol{\varepsilon}^{(i)}, \quad (29) $$

where $\mathbf{D}_k$ is the square root of the estimated covariance matrix, $\boldsymbol{\varepsilon}^{(i)}$ is distributed according to some kernel $K$ and $\hat{\mathbf{x}}$ is the output vector of the multino- mial draw. The Epanechnikov kernel $K$ is used, with the optimal bandwidth $h$ set according to the value specified by (Silverman, 2018).

## 4.5 Estimation output

The global estimate of the state vector $\mathbf{x}_k$, denoted $\hat{\mathbf{x}}_k$, and its associated estimated covariance matrix denoted, $\hat{\mathbf{P}}_k$, are respectively given by:

$$ \hat{\mathbf{x}}_k = \sum_{i=1}^{N} w_k^{(i)} \mathbf{x}_k^{(i)} \quad (30) $$

and

$$ \hat{\mathbf{P}}_k = \sum_{i=1}^{N} w_k^{(i)} (\mathbf{x}_k^{(i)} - \hat{\mathbf{x}}_k) (\mathbf{x}_k^{(i)} - \hat{\mathbf{x}}_k)^{\top}. \quad (31) $$

In this paper, the global estimate is computed after the regularization-resampling step, whether or not the step is performed.

## 4.6 Regularized particle filter for target motion analysis

The adaptation from the traditional BOT problem to TMA with combined bearing and SPL is straightforward, requiring only a modification of the likelihood function to incorporate the noise level data. Given the positions of both the target and the observer, the angle measurement and the noise level measurement can be considered independent. As a result, the likelihood can be computed by combining the two independent observations as follows:

$$ p(\mathbf{z}_k | \mathbf{x}_k, \mathbf{x}_k^o) = p(\theta_k | \mathbf{x}_k, \mathbf{x}_k^o) p(l_k | \mathbf{x}_k, \mathbf{x}_k^o), \quad (32) $$

with $p(\theta_k | \mathbf{x}_k, \mathbf{x}_k^o)$ and $p(l_k | \mathbf{x}_k, \mathbf{x}_k^o)$ respectively proportional to:

$$ p(\theta_k | \mathbf{x}_k, \mathbf{x}_k^o) \propto \exp \left( - \frac{(\theta_k - h_\theta(\mathbf{x}_k))^2}{2\sigma_\theta^2} \right) \quad (33) $$

and

$$ p(l_k | \mathbf{x}_k, \mathbf{x}_k^o) \propto \exp \left( - \frac{(l_k - h_l(\mathbf{x}_k + \mathbf{x}_k^o, \mathbf{x}_k^o))^2}{2\sigma_l^2} \right). \quad (34) $$

Since frequency measurements can be considered conditionally independent of the other types of measurements, their information can be fused in the same way.

# 5 Numerical results

## 5.1 Description of the scenarios

The scenarios used in the simulations are directly based on those described in Ristic et al. (2003, Chapter 6). However, in the scenarios described in subsections 5.1.2 and 5.1.3, the observer trajectory is changed and the time horizon is extended to provide a more complete visualization of the results.

In all scenarios, the target is initially positioned at [5000 800] in a Cartesian coordinate system where the coordinates are measured in meters, with a depth of 10 m. It maintains a constant speed of $\approx 2.57 \text{ m s}^{-1}$ (5 kn) and a steady course of $-140^\circ$.

The observer always starts at position [0 0]. However, his depth varies depending on the type of measurement considered. When the observer receives bearing measurements — as described in (9) —, it is assumed that his depth is the same as the target’s depth. In contrast, when considering conical angle measurements — as described in (12) —, it is assumed that the observer is at a depth of 100 m to present a more realistic scenario of a submarine tracking a surface ship.

The particles are initialized with a prior distribution for the range centered at 3 km beyond the true range, with a standard deviation of 3 km. This deliberately biased choice enables the assessment of the
robustness of the different methods. The initial di-
rection is determined by the first angle measurement,
with the angular standard deviation σ<sub>θ</sub> set to 1.5°.
The initial speed is drawn from a Gaussian distri-
bution centered at ≈1.54 m s<sup>-1</sup> (3 kn), with a stan-
dard deviation of ≈1.03 m s<sup>-1</sup> (2 kn). The initial
course is set to z<sub>0</sub>+180° with a standard deviation
of 180/&radic;12°. The complete initialization procedure
is detailed in Ristic et al. (2003, Section 6.4.1.1).

The filters use N = 5000 particles, with a resam-
pling threshold Nthr set to N/3. Measurements are
received every 60 s. The noise in the signal level is
generated with a standard deviation σl of 3 dB. The
standard deviation of the measured frequency is set
to σf=0.01 Hz. The process noise σs associated to
the emitted frequency is set to 0.001 Hz. The pro-
cess noise σa associated to the acceleration is set to
1.5 × 10<sup>-3</sup> m s<sup>-2</sup>.

The likelihoods within the filters are calculated un-
der the assumption that the standard deviations σθ,
σf and σl are 1.5 times larger than their actual val-
ues. On the one hand, this introduces more flexibility
into the algorithm, and on the other hand, it reflects a
more realistic scenario, since the exact real measure-
ment noises are often unknown.

The environmental data file used as input for bell-
hopcxx is provided in the Appendix.

For each simulation, the evolution of the root-mean-
square error (RMSE) for the target’s position over
1000 runs is presented, comparing angle-only mea-
surements with different combinations of frequency
and SPL measurements. A SPL-only tracking method
was also tested, but failed to track the target. Due to
its divergence, the RMSE for this method has been
excluded from the figures for clarity.

**5.1.1 Maneuvering observer**

In these simulations, the observation period lasts
30 min. The observer travels at a constant speed of
≈2.57 m s⁻¹ (5 kn). It starts with an initial course of
140° and executes a maneuver between 13 min and
17 min to attain a new course of 20°. The observer
maintains this new course for the remainder of the
observation period. This scenario is illustrated in
FIG. 5.

**5.1.2 Non-maneuvering observer**

In these simulations, the observation period lasts
60 min. The observer travels at a constant speed
of ≈2.57 m s⁻¹ (5 kn) on a course of −140°, travel-
ing parallel to the target. In this specific configura-
tion, the quality of the bearing data is significantly
degraded due to the near-zero bearing rate, making it
impossible for standard BOT methods to produce ac-
curate or reliable estimates. As the traditional BOT
problem is unobservable (see Nardone et al., 1984),
this scenario creates a significant challenge for TMA

<figure><img src="image_6.png" /><figcaption>Figure 5: Scenario with maneuvering ownship.</figcaption></figure>

with combined bearing and SPLs. The scenario is
illustrated in FIG. 6.

<figure><img src="image_7.png" /><figcaption>Figure 6: Scenario with non-maneuvering ownship.</figcaption></figure>

**5.1.3 Static observer**

In these simulations, the observation period extends
over 60 min, with the observer remaining stationary at
its initial position for the entire duration. Similar to
the non-maneuvering observer scenario, the bearing
information is of poor quality and the classical BOT
problem is unobservable. This scenario is illustrated
in FIG. 7.

All three of these scenarios are tested with both
bearing measurements in Section 5.2 and conical angle
measurements in Section 5.3.

**5.2 Performance analysis using bear-
ing measurements**

**5.2.1 Maneuvering observer**

The evolution of the RMSEs throughout this time
frame is illustrated in FIG. 8.

In this scenario, the integration of sound pressure
level and frequency data has only a limited influ-
ence on filter performance. The RMSEs follow sim-
<figure><img src="image_8.png" /><figcaption>Figure 7: Scenario with static ownership.</figcaption></figure>

ilar trends for every filter and maintain close values throughout the experiment.

Before the maneuver begins, the estimates from each filter show some divergence. While adding frequency data alone has little effect, combining frequency and SPL data yields the most accurate estimates. However, at this stage of the scenario, using SPL data alone slightly degrades the estimates. This performance could potentially be improved with specialized filters designed for map correlation, as discussed earlier. However, this degradation diminishes as the simulation progresses, demonstrating the resilience of the algorithm.

As soon as the maneuver begins, all filters show a sharp and significant improvement in performance, with comparable gains in all cases. After the maneuver, however, filters that incorporate SPL measurements achieve slightly better accuracy. Notably, the combination of SPL and frequency data consistently produces the best results throughout the scenario. However, the SPL contributes more significantly to the reduction of the RMSE than the frequency does.

Thus, this case shows that transmission losses become less critical when bearing measurements contain substantial information — in the information-theoretic sense Hammel *et al.* (1989) —, though they still offer some benefit. To further the analysis, the non-maneuvering scenario was examined, where the bearing measurements exhibited a lack of information.

### 5.2.2 Non-maneuvering observer

As shown in Fig. 9, traditional BOT methods struggle when data lacks sufficient information content — in the information-theoretic sense Hammel *et al.* (1989) —, resulting in consistent divergence and increasingly inaccurate estimates over time. This is consistent with theoretical results indicating that in this scenario, bearing information alone is insufficient for accurate target position estimation.

In this scenario, frequency alone does not enhance the estimate, underscoring the extremely low information content of the data available to the observer.

<figure><img src="image_9.png" /><figcaption>(a) Time Interval: 0 min to 30 min.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) Detailed View: 20 min to 30 min.</figcaption></figure>

<caption>Figure 8: Evolution of RMSE for a maneuvering observer and the simplified bearing model.</caption>

However, incorporating SPLs data helps mitigate these problems. By combining bearing measurements with SPLs information, the tracking system significantly slows the growth of the error. Although the error does not decrease over time, this scenario emphasizes the importance of considering SPLs data when bearing information is unreliable, at least to contain the error. In fact, despite identical initialization, neglecting SPLs data results in a 51.73% deterioration of the estimate after 60 min.

It is interesting to note that throughout this scenario, the distance between the observer and the target remains of the same order of magnitude. However, as illustrated in FIG. 4, a wide range of distances between the target and the observer results in similar transmission losses to those measured. This phenomenon is explained by the chosen seabed profile, as well as by the observer's depth and the order of magnitude of the target-observer distance. It is therefore likely that this uncertainty prevents the methods using SPL data from achieving better convergence. To test this hypothesis, the static scenario is examined, in which the lack of information in the bearing measurements persists, but the distance between the target and the observer varies more significantly.
<figure><img src="image_11.png" /><figcaption>Figure 9: Evolution of RMSE for a Non-Maneuvering Observer and the simplified bearing model.</figcaption></figure>

### 5.2.3 Static observer

<figure><img src="image_12.png" /><figcaption>Figure 10: Evolution of RMSE for a Static Observer and the simplified bearing model.</figcaption></figure>

As shown by the evolution of RMSEs in FIG. 10, all filters initially improve their estimation accuracy during the first 30 min. At this point, the filter using only bearing information reaches its lowest RMSE value of 2062 m, while the filter using sound pressure level (SPL) in addition to bearing data achieves a lower error of 264 m.

For the BOT method, the estimation error begins to increase steadily after this time, eventually leading to a complete loss of the target. In contrast, when SPLs data are included, the estimation accuracy continues to improve throughout the scenario, reaching its lowest RMSE at 60 min with a value of 110 m, while the BOT method has already deteriorated to 3960 m.

Consistent with the theoretical results mentioned earlier, it was observed that, in this scenario with a static observer, a filter combining bearing and frequency information provided an estimate of the target's position with progressively improving accuracy. However, while this combination is effective, the inclusion of SPL data further improves accuracy, resulting in a 76.72% improvement by the end of the scenario.

This significant improvement highlights the potential of the method and represents a promising direction for future applications. It also suggests that the incorporation of SPLs is most effective when noise levels vary significantly, as it helps to better distinguish the possible distances between the target and the observer.

## 5.3 Performance analysis using conical angle measurements

### 5.3.1 Maneuvering observer

The evolution of the RMSEs throughout this time frame is illustrated in FIG. 11.

<figure><img src="image_13.png" /><figcaption>Figure 11: Evolution of RMSE for BOT for a maneuvering observer and the bottom bounce model.</figcaption></figure>

As in the simplified bearing model, the incorporation of SPL data in this scenario has a limited effect on the performance of the filters, with RMSEs following similar trends for all filters. The overall behavior of the estimation error, however, is different. The information provided by the conical angle eliminates the initial increase in the error and replaces it with a sharp decrease, followed by a stable phase. This highlights how accounting for the bottom bounce path can reduce range ambiguity — assuming the bounce
is known — especially when the observer is not maneuvering.

Interestingly, in this case, the frequency information complements the angular measurements more effectively, while the SPL data does not seem to contribute to the reduction of the RMSE. Nevertheless, the filter that integrates all three sources of information maintains a consistently strong performance.

### 5.3.2 Non-maneuvering observer

<figure><img src="image_14.png" /><figcaption>Figure 12: Evolution of RMSE for Non-Maneuvering Observer and the bottom bounce model.</figcaption></figure>

In this scenario, the use of conical angle measurements has a significant effect on the results, influencing the overall outcome in a meaningful way as shown in FIG. 12. In contrast to the simplified bearing case, where the best achievable result was a significant reduction in error divergence, here the error rapidly decreases to very low values across all filters within a few minutes — an outcome made possible under the assumption that the received signal is a bottom bounce.

Initially, all filters maintain a stable error level. After about 30 minutes, however, the influence of the sound pressure level data becomes critical. Filters that ignore this information begin to see their errors increase, while those that incorporate SPL maintain stable accuracy. The effect is so significant that even the most favorable result without SPL leads to a 375.92% deterioration of the estimate after 60 min, compared to methods that incorporate SPL.

This scenario highlights both the effectiveness of conical angle measurements in resolving positional ambiguities when classical bearing data are unreliable, and the valuable complementary role of SPL in enhancing angular information.

### 5.3.3 Static observer

As shown by the evolution of RMSEs in FIG. 13, all filters initially improve in estimation quality during the first 8 min. Beyond this point, however, only those that incorporate SPL measurements continue to show a decreasing error, while the others begin to diverge.

<figure><img src="image_15.png" /><figcaption>Figure 13: Evolution of RMSE for a Static Observer and the bottom bounce model.</figcaption></figure>

Although the addition of frequency data alone does not improve the consistency of the target position's estimate, it proves to be valuable when combined with SPL data. In particular, only the filter using all three types of measurements — bearing, SPL, and frequency — exhibits a continuously decreasing error throughout the 60 min simulation, except for a slight increase towards the end, while consistently yielding the best results overall.

Once again, this scenario highlights the importance of SPL data and reinforces that the optimal approach for achieving both consistency and accuracy is to combine angle, sound pressure level, and frequency data.

## 6 Conclusion

This paper presented a comprehensive study on the impact of sound pressure level measurements in target motion analysis. By combining SPL data with two types of angle measurements—bearing and conical angle—along with Doppler-shifted frequency data, the filter performance was evaluated in three different scenarios: maneuvering observer, non-maneuvering observer, and static observer.

The results consistently demonstrate the critical role of SPL measurements in improving estimation accuracy. In most cases, especially when bearing measurements lacked information (e. g., due to a non-maneuvering observer), SPL data provided the only source of information ensuring reliable filter performances. However, despite these promising results, TMA cannot rely solely on SPL measurements.

As expected, frequency data were found to improve estimation accuracy in certain situations. However, its contribution was consistently less significant than that of SPL measurements. Nevertheless, frequency information remains valuable, as filters incorporating angle measurements, frequency data and SPL data proved to be the most effective and consistent.
Finally, conical angle measurements were tested and, under the given assumptions in this paper, inherently provided additional information compared to traditional bearing data, particularly with respect to range. However, even in this more realistic setting, where angle measurements alone provide good estimation accuracy, the integration of SPL data remained essential to ensure the lowest possible estimation error.

It is important to note that the study was conducted under several simplifying assumptions. For example, an isotropically radiating source was assumed, and the directivity index was neglected, which may affect the method's performance in more complex acoustic environments. Additionally, for the conical angle, the received signal path type (bottom bounce) was assumed to be known, and no left-right ambiguity existed with respect to the sensor used.

Another limitation of the study is the assumption that the depth of the target is known. While this assumption is not a major challenge when dealing with surface ships, it could introduce complexity in other scenarios. In addition, the analysis focused on targets in rectilinear motion at constant speed to avoid introducing filter-dependent parameters, thus making comparisons more straightforward. However, applying this approach to maneuvering targets using specialized filters could yield promising results and provide further insight into the robustness of the method.

## Author Declarations

The author Enzo Iglésis has no conflicts to disclose. Luc de Montella reports grants from Naval Group through the CIFRE PhD program during the conduct of the study.

## Data Availability

The data supporting the findings of this study are provided within the article.

*

# A BELLHOP Environmental data

The environmental data use as input of bellhopccx are the following:

<table><tr><td>Munk profile</td><td>! TITLE</td></tr><tr><td>60</td><td>! FREQ</td></tr><tr><td>1</td><td>! NMEDIA</td></tr><tr><td>CVM</td><td>! SSPOPT</td></tr><tr><td>51 0.0 2000</td><td>! DEPTH of bottom</td></tr><tr><td>0.0 1548.72</td><td>0.0 /</td></tr><tr><td>100.0 1522.00</td><td>0.0 /</td></tr><tr><td>200.0 1507.97</td><td>0.0 /</td></tr><tr><td>300.0 1501.65</td><td>0.0 /</td></tr><tr><td>400.0 1500.00</td><td>0.0 /</td></tr><tr><td>500.0 1501.18</td><td>0.0 /</td></tr></table>

<table><tr><td>600.0</td><td>1504.08</td><td>0.0 /</td></tr><tr><td>700.0</td><td>1508.03</td><td>0.0 /</td></tr><tr><td>800.0</td><td>1512.60</td><td>0.0 /</td></tr><tr><td>900.0</td><td>1517.56</td><td>0.0 /</td></tr><tr><td>1000.0</td><td>1522.75</td><td>0.0 /</td></tr><tr><td>1100.0</td><td>1528.09</td><td>0.0 /</td></tr><tr><td>1200.0</td><td>1533.50</td><td>0.0 /</td></tr><tr><td>1300.0</td><td>1538.97</td><td>0.0 /</td></tr><tr><td>1400.0</td><td>1544.47</td><td>0.0 /</td></tr><tr><td>1500.0</td><td>1550.00</td><td>0.0 /</td></tr><tr><td>1600.0</td><td>1555.53</td><td>0.0 /</td></tr><tr><td>1700.0</td><td>1561.07</td><td>0.0 /</td></tr><tr><td>1800.0</td><td>1566.61</td><td>0.0 /</td></tr><tr><td>1900.0</td><td>1572.16</td><td>0.0 /</td></tr><tr><td>2000.0</td><td>1577.70</td><td>0.0 /</td></tr><tr><td>A 0.0</td><td colspan="2">! Bottom model</td></tr><tr><td>2000.0</td><td>1600.0 0.0 1.8 0.8</td><td>! Halfspace</td></tr><tr><td>1</td><td colspan="2">! NSD</td></tr><tr><td>10 /</td><td colspan="2">! SD(1:NSD)</td></tr><tr><td>2000</td><td colspan="2">! NRD</td></tr><tr><td>0.0 1999.0 /</td><td colspan="2">! RD(1:NRD)</td></tr><tr><td>5000</td><td colspan="2">! NR</td></tr><tr><td>0.0 49.99 /</td><td colspan="2">! R(1:NR)</td></tr><tr><td>IG</td><td colspan="2">! Run type</td></tr><tr><td>0</td><td colspan="2">! NBEAMS</td></tr><tr><td>-90.0 90.0 /</td><td colspan="2">! ALPHA(1:NBEAMS)</td></tr><tr><td>0.0 2100.0 50.09</td><td colspan="2">! STEP, ZBOX, RBOX</td></tr></table>

## References

Anderson, B., and Moore, J. (2012). Dover Books on Electrical Engineering *Optimal Filtering* (Dover Publications).

Arulampalam, S., Ristic, B., Gordon, N., and Mansell, T. (2004). "Bearings-only tracking of manoeuvring targets using particle filters," EURASIP Journal on Advances in Signal Processing 2004, 1–15.

Curtis, K., Howe, B., and Mercer, J. (1999). "Low-frequency ambient sound in the North Pacific: Long time series observations," The Journal of the Acoustical Society of America 106(6), 3189–3200.

Elfring, J., Torta, E., and van de Molengraft, R. (2021). "Particle filters: A hands-on tutorial," Sensors 21(2).

Farina, A. (1999). "Target tracking with bearings-only measurements," Signal processing 78(1), 61–78.

Hammel, S. E., Liu, P. T., Hilliard, E. J., and Gong, K. F. (1989). "Optimal observer motion for localization with bearing measurements," Computers & Mathematics with Applications 18(1-3), 171–180.

Jalkanen, J.-P., Johansson, L., Liefvendahl, M., Bensow, R., Sigray, P., Östberg, M., Karasalo, I., Andersson, M., Peltonen, H., and Pajala, J. (2018). "Modelling of ships as a source of underwater noise," Ocean Science 14(6), 1373–1383.

Jauffret, C., and Pérez, A.-C. (2018). "Efficient pseudolinear estimation for bearings and frequencies target motion analysis," in 2018 21st International Conference on Information Fusion (FUSION), pp. 2556–2563.

Julier, S., Uhlmann, J., and Durrant-Whyte, H. (2000). "A new method for the nonlinear transformation of means and covariances in filters and estimators," IEEE Transactions on Automatic Control 45(3), 477–482.

Kong, A., Liu, J. S., and Wong, W. H. (1994). "Sequential imputations and bayesian missing data problems," Journal of the American statistical association 89(425), 278–288.
Laneuville, D. (2007). "Target motion analysis with passive data fusion,".

Ma, T., Ding, S., Li, Y., and Fan, J. (2023). "A review of terrain aided navigation for underwater vehicles," Ocean Engineering **281**, 114779.

Maranda, B. (2008). *Passive Sonar*, 1757–1781 (Springer New York, New York, NY).

Merlinge, N., Dahia, K., Piet-Lahanier, H., Brusey, J., and Horri, N. (2019). "A box regularized particle filter for state estimation with severely ambiguous and non-linear measurements," Automatica **104**, 102–110.

Munk, W. (1974). "Sound channel in an exponentially stratified ocean, with application to sofar," The Journal of the Acoustical Society of America **55**(2), 220–226.

Musso, C., Oudjane, N., and Le Gland, F. (2001). *Improving Regularised Particle Filters*, 247–271 (Springer New York, New York, NY).

Nardone, S., and Aidala, V. (1981). "Observability criteria for bearings-only target motion analysis," IEEE Transactions on Aerospace and Electronic Systems **AES-17**(2), 162–166.

Nardone, S., Lindgren, A., and G., K. (1984). "Fundamental properties and performance of conventional bearings-only target motion analysis," IEEE Transactions on Automatic Control **29**(9), 775–787.

Oh, R., Song, T., and Choi, J. (2020). "Batch processing through particle swarm optimization for target motion analysis with bottom bounce underwater acoustic signals," Sensors **20**(4), 1234.

Ovalle, D., and García-Pelaez, J. (2015). "Minimizing underwater noise generated by submarine maneuvering: an optimal control approach," IEEE Journal of Oceanic Engineering **41**(2), 362–372.

Palmier, C., Dahia, K., Merlinge, N., Del Moral, P., Laneuville, D., and Musso, C. (2019). "Adaptive approximate bayesian computational particle filters for underwater terrain aided navigation," in *2019 22th International Conference on Information Fusion (FUSION)*, IEEE, pp. 1–8.

Passerieux, J., Pillon, D., Blanc-Benon, P., and Jauffret, C. (1989). "Target motion analysis with bearings and frequencies measurements via instrumental variable estimator (passive sonar)," in *International Conference on Acoustics, Speech, and Signal Processing*, pp. 2645–2648 vol.4.

Passerieux, J., Pillon, D., Blanc-Benon, P., and Jauffret, C. (1988). "Target motion analysis with bearings and frequencies measurements," in *Twenty-Second Asilomar Conference on Signals, Systems and Computers*, Vol. 1, pp. 458–462.

Pisha, L., Snider, J., Jackson, K., and Jaffe, J. (2023). "Introducing bellhopcxx/bellhopcuda: Modern, parallel bellhop(3d)," The Journal of the Acoustical Society of America **153**(3 supplement), A218–A218.

Porter, M., and Bucker, H. (1987). "Gaussian beam tracing for computing ocean acoustic fields," The Journal of the Acoustical Society of America **82**(4), 1349–1359.

Ristic, B., Arulampalam, S., and Gordon, N. (2003). *Beyond the Kalman Filter: Particle Filters for Tracking Applications* (Artech House).

Ristic, B., Houssineau, J., and Arulampalam, S. (2019). "Robust target motion analysis using the possibility particle filter," IET Radar, Sonar & Navigation **13**(1), 18–22.

Ross, D. (1987). *Mechanics of Underwater Noise* (Peninsula Publishing).

Silverman, B. W. (2018). *Density estimation for statistics and data analysis* (Routledge).

Sindhu, B., Jayaraman, V., and Christopher, S. (2019). "Bear-ing only target tracking using single and multisensor: A re-view," Journal of Engineering Science and Technology Re-view **12**, 50–65.

Urick, R. (1983). *Principles of underwater sound — 2. ed.* (Peninsula Publishing).

Wittekind, D. (2014). "A simple model for the underwater noise source level of ships," Journal of Ship Production and Design **30**.