\newdateformat

monthyeardate\monthname[\THEMONTH] \THEDAY, \THEYEAR

Navigation in Shallow Water
Using Passive Acoustic Ranging

Junsu Jang and Florian Meyer

Email: {jujang, flmeyer}@ucsd.edu Scripps Institution of Oceanography and Department of Electrical and Computer Engineering University of California San Diego, La Jolla, CA

Abstract

Passive acoustics can provide a variety of capabilities with applications in oceanographic research and maritime situational awareness. In this paper, we develop a method for the navigation of autonomous underwater vehicles (AUVs) in shallow water. Our approach relies on passively recorded signals from acoustic sources of opportunities (SOOs). By making use of the waveguide invariant, a measurement of the range to the SOO is extracted from the spectrogram of a single hydrophone. Range extraction requires knowledge of the range rate, i.e., the radial velocity between the SOO and AUV, computed from the pressure fields at different time intervals. A particle-based navigation filter fuses the range measurements with the AUV’s internal velocity and heading measurements. As a result, the position error, which would otherwise increase over time, can be bounded. The ability to compute the range rate and range measurements from the pressure field measured using a single hydrophone is demonstrated on real data from the SWellEx-96 experiment. The capability of the developed navigation filter is shown based on synthetic data generated by the normal mode program KRAKEN.

Index Terms:

Navigation, Bayes Filtering, Passive Acoustics, Marine Robotics, Particle Filtering.

I Introduction

Autonomous underwater vehicles (AUVs) are widely used in scientific, commercial, and military applications due to their ability to perform complex tasks effectively. However, navigation of small and inexpensive AUVs remains a challenge. This is because underwater, the radio signals of the Global Positioning System (GPS) are unavailable. Navigation based on sensors that do not require wireless links, such as inertial sensors, relies on integrating acceleration or velocity measurements. This strategy, referred to as dead-reckoning, is subject to position errors that grow over time. As a result, the AUV resurfaces regularly to perform GPS localization [1]. There exist a variety of active acoustic underwater ranging techniques [2][3], but they have limited range and require deployments of multiple transponders.

In this paper, we investigate navigation by means of a passive acoustic ranging technique. Our approach only relies on a simple acoustic sensor hosted by the AUV. Since acoustic signals emitted by a source of opportunity (SOO) are used, no additional transponders are deployed. The proposed method is applicable in shallow water with frequent ship traffic, such as coastal waterways. The SOOs are typically commercial vessels whose propeller cavitation causes the ship noise to be used as acoustic signal [4]. In addition, this type of vessel is equipped with an automatic identification system (AIS) for broadcasting its current and predicted future position information[5].

Refer to caption — Figure 1: Bayesian network representing the proposed Bayes filter for AUV navigation. The state of the AUV at discrete time $n$ , $\bm{x}_{n}$ , consists of the AUV’s position, heading, and velocity in a global coordinate system. A measurement, $\bm{z}_{n}$ , consists of measurements of the range with respect to an SOO, the local velocity, and the heading. The range measurement is computed from the acoustic signal of an SOO in shallow water using the WI. The WI-based ranging relies on a range rate measurement computed in a preprocessing stage that follows the approach presented in [6].

Our method extracts range measurements from the acoustic signal recorded by a single hydrophone. In a shallow-water waveguide, the signal emitted from a broadband moving source leads to an interference pattern in the spectrogram of a single hydrophone. In particular, minima and maxima of signal intensity are slightly shifted with respect to frequency, $f$ , and range, $r$ [7]. Loci of constant intensities form a diagonal striation pattern characterized by a scalar parameter called the waveguide invariant (WI), $\beta$ [8]. The WI encapsulates the dispersive propagation characteristics of the waveguide[9].

Parameter estimation problems in underwater acoustics that rely on striation patterns described by the WI include ranging and source localization [10, 11, 12, 13], geoacoustic inversion[14], and internal wave characterization[15]. In [16], the WI is treated as a probability distribution, and the impacts of environmental parameters (depth of the source and receiver, sound speed, etc…) on the WI distribution are investigated.

The WI-based ranging with respect to SOOs that exhibit tonal signals has been explored in[17, 18] A single hydrophone is considered, and the application of the WI-based ranging to AUV navigation is discussed. Here, a statistical model for intensities along hypothetical striations is established, and the maximum likelihood (ML) estimation of the WI parameter and range is performed. A nonlinear least square-based approach for receiver localization is also proposed.

Range estimation based on the WI requires knowledge of the range rate, which can be time-varying. They are, however, assumed known in the aforementioned ranging approaches. A technique for estimating the range rate from the acoustic measurements of a single hydrophone has been presented in [6]. Assuming a narrowband source, this approach relies on normal mode theory[7] to extract range rate from the spectra of consecutive snapshots of data. A combination of the range rate estimation with the WI-based ranging is also discussed. The range rate estimation has been extended to scenarios with broadband sources in [19, 20].

The fundamental question addressed in this paper is the feasibility of underwater self-localization in shallow water based on the range measurements extracted from the data of a simple acoustic receiver. In the considered scenario, tonal signals are emitted by a moving SOO. We develop a self-localization method that relies on the WI-based ranging to bound the position error of an AUV in a shallow-water environment with frequent ship traffic. Our Bayesian filter fuses the range measurements with the velocity and heading measurements provided by the internal sensors of the AUV. It offers an opportunity for an AUV, which would otherwise solely perform dead-reckoning, to recalibrate its position information without resurfacing. The WI-based ranging relies on an acoustically range-independent underwater environment, passively recorded tonal signals emitted by the SOO, and the knowledge of the potentially time-varying range rate.

Within our approach, a range rate measurement is computed in a preparatory step [6] by using the same tonal signal as used for the WI-based ranging. Only a single hydrophone is required to obtain a range measurement if a single SOO is present. This makes the proposed method appealing for inexpensive AUV platforms with limited acoustic sensing capabilities. We illustrate the accuracy of the range rate and range measurements using real acoustic data recorded during the SWellEx-96 experiment[21]. We also demonstrate the capabilities of the proposed Bayesian navigation filter using acoustic data simulated in a realistic shallow water environment generated by the normal mode program KRAKEN. In particular, we show how the position error can be bounded and thus reduced compared to a reference method that relies on dead-reckoning.

The key contributions of this paper are summarized as follows:

•

We review the approaches for computing measurements of the range rate, the range, and the WI from the spectrogram recorded by a single hydrophone that is mobile with respect to an SOO.
•

We develop a Bayesian filter for navigation in shallow water that fuses information provided by the sensors for dead reckoning with the range measurements computed from the spectrogram.
•

We demonstrate passive acoustic ranging using real data from the SWellEx-96 experiment and present the capability of the developed navigation filter using synthetic data generated by the normal mode program KRAKEN.

Notation: Random variables are displayed in serif, italic fonts, and vectors and matrices are denoted by bold lowercase and uppercase letters, respectively. For example, a random variable and a random vector are denoted by $x$ and $\bm{x}$ , respectively. Furthermore, $\|\bm{x}\|$ and ${\bm{x}}^{\mathrm{T}}$ denote the Euclidean norm and the transpose of vector $\bm{x}$ , respectively; $\propto$ indicates equality up to a constant factor; $\mathpzc{p}(\bm{x})$ denotes the probability density function (PDF) of random vector $\bm{x}$ (this is a short notation for $\mathpzc{p}_{\bm{x}}(\bm{x})$ ); $\mathpzc{p}(\bm{x}|\bm{y})$ denotes the conditional PDF of random vector $\bm{x}$ conditioned on random vector $\bm{y}$ (this is a short notation for $\mathpzc{p}_{\bm{x}|\bm{y}}(\bm{x}|\bm{y})$ ). $\bm{x}_{0:k}$ is short for $[\bm{x}_{0}^{\mathrm{T}},\dots,\bm{x}_{k}^{\mathrm{T}}]^{\mathrm{T}}$ . $\bm{I}_{n}$ denotes the $n\hskip-0.85358pt\times\hskip-0.85358ptn$ identity matrix. $|\mathcal{S}|$ denotes the cardinality of set $\mathcal{S}$ . The operator ^∗ denotes the complex conjugate, and $\mathpzc{j}=\sqrt{-1}$ is the imaginary unit.

II WI-Based Ranging

In this section, we will present our approach to extracting the range measurements from acoustic data by exploiting the striation pattern described by the WI. First, we discuss the preprocessing stage that will provide the range rate measurements subsequently used for the WI-based ranging. We also establish an initialization stage for the computation of the WI. In what follows, it is assumed that consecutive snapshots of data samples have been extracted from the raw acoustic signal and that the discrete Fourier transform (DFT) of length $N_{\mathrm{DFT}}$ and overlap $\alpha$ has been performed for each snapshot. The time between snapshots is $t_{\Delta}=(1-\alpha)N_{\mathrm{DFT}}\hskip 0.85358ptT_{\Delta}$ , where $T_{\Delta}$ is the sampling interval of the raw acoustic signal. Furthermore, $n=0,1,2,\dots$ denotes the discrete-time indexes of snapshots.

II-A Computation of Range Rate Measurements

At each time step $n$ , we compute a measurement of the range rate, $\dot{r}_{n}$ , following the “difference” method proposed in [6]. Here, the complex DFT samples at source frequency $f$ and from consecutive snapshots of data result in the signal $p[n,f]$ . This signal represents the acoustic pressure field as a function of range, $r$ , and frequency, $f$ . To compute a range rate measurement, we consider two segments of snapshots before and after time $n$ , i.e., $p_{n,1}[k,f]=p[n\hskip-0.85358pt-\hskip-0.85358ptk,f]$ and $p_{n,2}[k,f]=p[n\hskip-0.85358pt+\hskip-0.85358ptk,f]$ with $k\hskip-0.85358pt\in\hskip-0.85358pt\{0,\dots,L\}$ . Both segments include the sample at time $n$ and are $L+1$ snapshots long, i.e., a total of $2L+1$ snapshots are considered.

A measurement of $\dot{r}_{n}$ can now be extracted from the real part of the point-wise product of $p_{n,1}[k,f]$ with $p^{*}_{n,2}[k,f]$ , i.e., $I_{\mathrm{c}}[k;n,f]=\mathcal{R}\big{\{}p_{n,1}[k,f]\hskip 0.85358ptp^{*}_{n,2}[k,f]\big{\}}$ . Following [6], the result of this product can be approximated as $I_{\mathrm{c}}[k;n,f]\hskip-0.85358pt\approx\hskip-0.85358pt\tilde{I}_{\mathrm{c}}[k;n,f]$ , where $\tilde{I}_{\mathrm{c}}[k;n,f]$ is given by

\tilde{I}_{\mathrm{c}}[k;n,f]\hskip 0.85358pt\propto\hskip 0.85358pt\cos{[wk]}.\vskip-1.42262pt

(1)

Here, we introduced $w\hskip-0.85358pt\triangleq\hskip-0.85358pt4\pi f\dot{r}_{n}t_{\Delta}/\overline{v_{\mathrm{p}}}$ where $\overline{v_{\mathrm{p}}}$ is the average horizontal phase speed [7]. Note that $\overline{v_{\mathrm{p}}}$ can be further approximated by the sound speed in the water, $v$ , which is assumed to be known [6]. Eq. (1) motivates the computation of a range rate measurement, $y_{\dot{r}_{n}}$ , from the location of largest magnitude in the frequency domain of $I_{c}[k;n,f]$ . In particular, let $F_{c}[w;n,f]$ be the DFT of $I_{c}[k;n,f]$ . A measurement of the range rate can then be obtained as

y_{\dot{r}_{n}}=\operatorname*{arg\,max}_{\dot{r}_{n}}|F_{\mathrm{c}}[4\pi f\dot{r}_{n}\hskip 0.85358ptt_{\Delta}/v;n,f]|^{2}.\vskip-4.2679pt

(2)

Note that the measurement $y_{\dot{r}_{n}}$ represents the average range rate over the $2L+1$ snapshots.

Unfortunately, the approach discussed above is not directly applicable to real-time navigation. This is because at time $n$ , the required signal $p_{n,1}[k,f]$ consists of samples that lie in the future. The most recent range rate measurement that can be computed at time $n$ is actually the outdated $y_{\dot{r}_{n-L}}$ . Thus, for real-time navigation, we use $y_{\dot{r}_{n-L}}$ to predict a real-time measurement at time $n$ , denoted as $z_{\dot{r}_{n}}$ . In particular, an average range acceleration measurement, $y_{\ddot{r}_{n-L}}$ , is extracted from the range rate measurements $y_{\dot{r}_{n-2L}},\dots,y_{\dot{r}_{n-L}}$ . The real-time range rate measurement $z_{\dot{r}_{n}}$ , is then obtained using a linear prediction, i.e.,

z_{\dot{r}_{n}}=y_{\dot{r}_{n-L}}+y_{\ddot{r}_{n-L}}t_{\Delta}\hskip 0.85358ptL.

(3)

Also note that the time interval $t_{\Delta}L$ has to be several minutes in length [6]. Only in this way a range rate measurement can be computed according to (2). Thus, It is assumed that the source signal is coherent within this time interval. Furthermore, note that a phase offset is introduced when the time interval between the snapshots, $t_{\Delta}$ , is not an integer multiple of $1/f$ [22]. This phase offset results in a shift of the peak location of $F_{\mathrm{c}}[w;n,f]$ . Hence, the phase compensation approach proposed in [22] is applied to obtain accurate range rate measurements.

II-B Computation of Range Measurements

In a shallow-water waveguide, the signal emitted from a broadband moving source leads to an interference pattern in the spectrogram of a single hydrophone. In particular, let us consider an observation interval $[t_{\mathrm{start}},t_{\mathrm{stop}}]$ and let $t_{\mathrm{start}}\leqslant t\leqslant t_{\mathrm{stop}}$ be an arbitrary time in the observation interval. We denote the ranges at times $t_{\mathrm{start}}$ , $t$ , and $t_{\mathrm{stop}}$ by $r_{\mathrm{start}}$ , $r$ , and $r_{\mathrm{stop}}$ , respectively. Assuming that the source’s range rate $\dot{r}$ is fixed and known during the observation interval, we can compute a spectrogram in the $(r,f)$ domain. Note that, since the range is typically unknown, the $r$ -axis of this spectrogram is a function of a reference range, e.g., with $r_{\mathrm{stop}}$ being the reference range, the range at time $t$ is equal to $r=r_{\mathrm{stop}}-(t_{\mathrm{stop}}-t)\dot{r}$ .

This $(r,f)$ spectrogram shows an inference pattern with striations of constant intensity. Each striation can be mathematically described as[7]

\frac{f_{0}}{f_{1}}=\Big{(}\frac{r_{0}}{r_{1}}\Big{)}^{\beta}\vskip 0.0pt

(4)

where $r_{0},r_{1}\in[r_{\mathrm{start}},r_{\mathrm{stop}}]$ and $f_{0},f_{1}$ are arbitrary source frequencies. According to (4), the same intensity as in $(r_{0},f_{0})$ can be found in any other point of the spectrogram, $(r_{1},f_{1})$ , where $r_{1}=r_{0}(f_{1}/f_{0})^{1/\beta}$ . The parameter $\beta$ is known as the WI. Since in typical shallow water waveguides, we have a $\beta$ close to $1$ , striations often appear as diagonal lines. An example with $\beta=1.07$ and a known $r_{\mathrm{stop}}$ is shown in Fig. 2(a). The convenience of exploiting the striation pattern described by the WI for ranging is that in an approximately range-independent waveguide, all parameters that characterize the environment are summarized by the single scalar WI, $\beta$ . As discussed in [17], striations are also visible when tonal signals are recorded, e.g., when a large container ship emits mechanical noise[23].

At time step $n$ , for the WI-based ranging, we aim to compute a measurement of the range between the source and the receiver, $r_{n}$ , by exploiting the WI[7][11][24]. Each column in the spectrogram corresponds to discrete time step $i\hskip-0.85358pt\in\hskip-0.85358pt\{0,\dots,n\}$ . By making use of the range rate measurements $z_{\dot{r}_{l}}$ , $l\in\{0,\dots,n-1\}$ , computed as discussed in the previous Sec. II-A, the following ranges can be assigned to each column of the spectrogram

r_{i}=r_{n}\hskip 0.85358pt\hskip 0.85358pt-\hskip-0.85358pt\hskip-0.85358pt\hskip-0.85358pt\hskip-0.85358pt\sum_{l=i}^{n-1}z_{\dot{r}_{l}}t_{\Delta},\hskip 14.22636pti\in\{0,\dots,n\}.

(5)

Here, we have used $r_{n}$ as the reference range.

Next, we develop an objective function for the computation of $r_{n}$ by pairwise comparing nonlinearly transformed intensity functions for different frequencies. In what follows, it is assumed that the WI, $\beta$ , is known. A method for the computation of a WI measurement $z_{\beta}$ will be discussed in Sec. II-C. Let $f_{0}$ and $f_{1}$ be the two source frequencies. For these two frequencies, the intensities as a function of the ranges, $r_{i}$ , $i\in\{0,\dots,n\}$ , are denoted as $I[r_{i};f_{0}]$ and $I[r_{i};f_{1}]$ , respectively. Each intensity function is represented by a single row in the spectrogram.

For the WI-based ranging, we compute cross-correlation coefficients for the most recent $N$ ranges, i.e., $r_{i},i\in\{n\hskip-0.85358pt-\hskip-0.85358ptN\hskip-0.85358pt+\hskip-0.85358pt1,\dots,n\}$ . Without loss of generality, let us assume that $f_{0}<f_{1}$ . Based on (4) and assuming a positive $\beta$ , the acoustic intensity at $f_{1}$ is a stretched function of that of $f_{0}$ . An example is shown in Fig. 2(b). First, we perform a similar stretching for $f_{0}$ by using (4) to convert the ranges of the last $M$ intensity samples as follows

r^{\prime}_{i}=r_{i}\Big{(}\frac{f_{1}}{f_{0}}\Big{)}^{1/\beta}\hskip-0.85358pt\hskip-0.85358pt\hskip-0.85358pt\hskip-0.85358pt,\hskip 14.22636pti\in\{n-M+1,\dots,n\}.\vskip 1.42262pt

(6)

This results in the stretched intensity function $I^{\prime}[r^{\prime}_{i};f_{0},\beta]$ . Note that $M\hskip-0.85358pt\geqslant\hskip-0.85358ptN$ has to be chosen such that the resulting stretched intensity function is defined on the entire interval $[r_{n-N+1},r_{n}]$ . Finally, we perform a linear interpolation to again obtain intensity for ranges $r_{i}$ , $i\in\{n-N+1,\dots,n\}$ . This interpolation implies a cutoff of ranges outside the interval $[r_{n-N+1},r_{n}]$ . As can easily be verified, the resulting intensity function depends on the value of the reference range $r_{n}$ and is thus denoted by $J[r_{i};f_{0},\beta,r_{n}]$ . An example of how an intensity function $J[r_{i};f_{0},\beta,r_{n}]$ matches $I[r_{i};f_{1}]$ after stretching and interpolation by the correct $\beta$ and $r_{n}$ , is shown in Fig. 2(c).

Given the correct $\beta$ and $r_{n}$ , the intensity pattern, e.g., the locations of minima and maxima of $J[r_{i};f_{0},\beta,r_{n}]$ , typically matches the intensity pattern of $I[r_{i};f_{1}]$ well. However, the absolute intensity values of the two intensity functions can be quite different, and there is a need for normalization. Therefore, the cross-correlation coefficient is used to measure similarity to develop an objective function. Let us denote the demeaned versions of $I[r_{i};f_{1}]$ and $J[r_{i};f_{0},\beta,r_{n}]$ as $\hat{I}[r_{i};f_{1}]$ and $\hat{J}[r_{i};f_{0},\beta,r_{n}]$ , respectively. As an objective function for the computation of a range measurement $z_{\dot{r}_{n}}$ , we now introduce the sum of the cross-correlation coefficients for all available pairs of demeaned intensity functions and the most recent $N$ samples, i.e.,

\displaystyle g(r_{n};\beta)=\sum_{(f_{0},f_{1})\in\mathcal{F}}\hskip 1.70717pt\frac{1}{Z}_{f_{0},f_{1}}\hskip 1.42262pt\sum_{i=n-N+1}^{n}\hskip-8.53581pt\hat{J}[r_{i};f_{0},\beta,r_{n}]\hskip 0.85358pt\hat{I}[r_{i};f_{1}],

where $\mathcal{F}$ is the set of all pairs of source frequencies with $f_{0}<f_{1}$ and $Z_{f_{0},f_{1}}$ is a normalization factor given by

Z_{f_{0},f_{1}}=\Biggl{[}\hskip 2.84526pt\sum_{i=n-N+1}^{n}\hat{J}^{2}[r_{i};f_{0},\beta,r_{n}]\sum_{i=n-N+1}^{n}\hskip-0.85358pt\hskip-0.85358pt\hat{I}^{2}[r_{i};f_{1}]\Biggr{]}^{\frac{1}{2}}\hskip-0.85358pt\hskip-0.85358pt\hskip-0.85358pt.

(8)

Finally, a range measurement is obtained by finding the location of the maximum of the objective function (LABEL:eq:ri), i.e.,

z_{r_{n}}=\operatorname*{arg\,max}_{r_{n}}g(r_{n};z_{\beta}).\vskip-1.42262pt

(9)

Here, $z_{\beta}$ is the WI measurement to be discussed in the next Sec. II-C. An example of the objective function maximizing is shown in Fig. 2(d). To reduce the computational complexity, unlike the likelihood function used in [18], the proposed objective function, (LABEL:eq:ri), does not consider noise statistics. However, as demonstrated in Sec. IV, it still leads to a competitive-ranging performance.

II-C Computation of the WI Measurement

A WI measurement, $z_{\beta}$ , can be computed in an initialization stage by following a similar approach as for the computation of $z_{r_{n}}$ . For initialization, it is assumed that the range with respect to the source, $r_{n}$ , is known. The initialization stage is performed while the AUV remains at the surface and can record GPS measurements. In particular, for a fixed $r_{n}$ , we obtain the WI measurement by finding the location of the maximum of $\beta$ in the objective function (LABEL:eq:ri), i.e.,

z_{\beta}=\operatorname*{arg\,max}_{\beta}g(r_{n};\beta).\vskip-2.84526pt

(10)

It is assumed that $\beta$ remains unchanged during the initialization stage and the next dive of the AUV. Depending on the sound speed profile, $\beta$ can change as a function of the receiver depth [16]. A method that can refine $\beta$ online is subject to future research.

III Statistical Model and Navigation Filter

In what follows, we introduce the statistical model for an AUV navigation and develop the proposed navigation filter. For notational simplicity, it is assumed that the navigation filter uses the same time scale as the computation of range rate and range measurements discussed in Sec. II, i.e., each discrete time step $n$ has a duration of $t_{\Delta}$ .

The state of the AUV at discrete time $n$ is denoted as $\bm{x}_{n}=\big{[}\bm{p}_{n}^{\mathrm{(W)}}\hskip 0.85358pt\hskip 0.85358pt\hskip 0.85358pt\bm{v}_{n}^{\mathrm{(W)}}\hskip 0.85358pt\hskip 0.85358pt\hskip 0.85358pt\stheta^{\mathrm{(W)}}_{n}\big{]}^{\mathrm{T}}\hskip-0.85358pt\hskip-0.85358pt\hskip-0.85358pt$ , where $\bm{p}_{n}^{\mathrm{(W)}}=\big{[}x_{n}^{\mathrm{(W)}}\hskip 0.85358pt\hskip 0.85358pt\hskip 0.85358pty_{n}^{\mathrm{(W)}}\big{]}^{\mathrm{T}}\hskip-0.85358pt\in\hskip-0.85358pt\mathbb{R}^{2}$ is position, $\bm{v}_{n}^{\mathrm{(W)}}=\big{[}\dot{x}_{n}^{\mathrm{(W)}}\hskip 0.85358pt\hskip 0.85358pt\hskip 0.85358pt\dot{y}_{n}^{\mathrm{(W)}}\big{]}^{\mathrm{T}}\hskip-0.85358pt\in\hskip-0.85358pt\mathbb{R}^{2}$ is the velocity, and $\theta^{\mathrm{(W)}}_{n}\hskip-0.85358pt\in\hskip-0.85358pt[0^{\circ},360^{\circ})$ is the heading (Fig. 3). All three quantities are defined with respect to a global (“world”) 2-D Cartesian coordinate system. In particular, the x- and y-axes are positive along the east and north directions, and the heading of 0° points east and increases counterclockwise. The superscripts (W) and (B) indicate the “world” and “body” reference frames, respectively. The heading is relevant for the measurement model where the local velocity measurement has to be rotated from the body reference frame to the world reference frame.

The state transition between time steps follows a nearly constant velocity motion model, i.e.,

\bm{x}_{n}=\begin{bmatrix}1&0&t_{\Delta}&0&0\\ 0&1&0&t_{\Delta}&0\\ 0&0&1&0&0\\ 0&0&0&1&0\\ 0&0&0&0&1\end{bmatrix}\bm{x}_{n-1}+\bm{u}_{n}\vskip 1.42262pt

(11)

where $t_{\Delta}$ is the length of each discrete time step and $\bm{u}_{n}\hskip-0.85358pt\in\hskip-0.85358pt\mathbb{R}^{5}$ is a zero-mean multivariate Gaussian random vector with covariance matrix given by

\bm{\varSigma}_{u}=\bm{U}\hskip-0.85358pt\hskip-0.85358pt\hskip-0.85358pt\hskip-0.85358pt\begin{bmatrix}\sigma_{\bm{u}_{1}}^{2}&0&0\\ 0&\sigma_{\bm{u}_{1}}^{2}&0\\ 0&0&\sigma_{\bm{u}_{2}}^{2}\end{bmatrix}\hskip-0.85358pt\hskip-0.85358pt\bm{U}^{\hskip-0.85358pt\mathrm{T}}\hskip-0.85358pt\hskip-0.85358pt\hskip-0.85358pt,\quad\bm{U}\hskip-0.85358pt=\hskip-0.85358pt\begin{bmatrix}\frac{t_{\Delta}^{2}}{2}&0&0\\ 0&\frac{t_{\Delta}^{2}}{2}&0\\ t_{\Delta}&0&0\\ 0&t_{\Delta}&0\\ 0&0&t_{\Delta}\end{bmatrix}\hskip-0.85358pt\hskip-0.85358pt\hskip-0.85358pt.\vskip 0.0pt

Here, $\sigma_{\bm{u}_{1}}^{2}$ and $\sigma_{\bm{u}_{2}}^{2}$ are the variance of a random linear acceleration and turn rate, respectively. The driving noise $\bm{u}_{n}$ is assumed statistically independent across time $n$ . The state-transition model in (11) induces the state-transition PDF $\mathpzc{p}(\bm{x}_{n}|\bm{x}_{n-1})$ that will be used in the prediction step of the proposed navigation filter. At the time $n\hskip-0.85358pt=\hskip-0.85358pt0$ , it is assumed that the prior PDF of the AUV state, i.e., $\mathpzc{p}(\bm{x}_{0})$ is available.

The measurements at time $n$ are denoted as $\bm{z}_{n}\hskip-0.85358pt=\hskip-0.85358pt\big{[}z_{r_{n}}\hskip 0.85358pt\hskip 0.85358pt\hskip 0.85358ptz_{\bm{v}_{n}^{\mathrm{(B)}}}\hskip 0.85358pt\hskip 0.85358ptz_{\varphi_{n}}\big{]}^{\mathrm{T}}$ , where $z_{r_{n}}$ is the range measurement that is computed as discussed in Sec. II, $z_{\bm{v}_{n}^{\mathrm{(B)}}}=\big{[}\dot{x}_{n}^{\mathrm{(B)}}\hskip 0.85358pt\hskip 0.85358pt\hskip 0.85358pt\dot{y}_{n}^{\mathrm{(B)}}\big{]}^{\mathrm{T}}$ is the velocity in the local “body” reference frame of the vehicle, and $z_{\varphi_{n}}$ is a compass reading, i.e., a noisy measurement of the AUV heading. As a statistical model for the range measurements, we use

\displaystyle z_{r_{n}}

\displaystyle=\|\bm{p}_{n}^{\mathrm{(W)}}-\bm{p}_{\mathrm{s},n}\|+\gamma_{n}.

(12)

Here, $\bm{p}_{\mathrm{s},n}$ is the position of the SOO in the global reference frame that is assumed known, and $\gamma_{n}$ is a zero-mean Gaussian measurement noise with variance $\sigma^{2}_{\gamma}$ . The measurement noise $\gamma_{n}$ is assumed statistically independent across $n$ . The initial position, $\bm{p}_{\mathrm{s},0}$ , and velocity, $\bm{p}_{v,0}$ , of the SOO are known, i.e., received over the automatic identification system (AIS) as the AUV starts to dive. Furthermore, it is assumed that the SOO stays on course, and thus, the position of the SOO at time $n$ can be predicted reliably, i.e., $\bm{p}_{\mathrm{s},n}=\bm{p}_{\mathrm{s},0}+n\hskip 0.85358pt\hskip 0.85358ptt_{\Delta}\hskip 0.85358pt\hskip 0.85358pt\bm{p}_{v,0}$ . Future work will include an improved model for passive acoustic ranging that consists of a statistical characterization of the processing performed in Sec. II and takes the uncertainty of knowing the position of the SOO into account.

The velocity measurements provided by the Doppler velocity log (DVL) of the AUV are modeled as

\displaystyle z_{\bm{v}_{n}^{\mathrm{(B)}}}

\displaystyle={}_{B\hskip-0.85358pt}R_{W}(\stheta^{\mathrm{(W)}}_{n})\hskip 0.85358pt\hskip 0.85358pt\bm{v}_{n}^{\mathrm{(W)}}+\bm{\xi}_{n},

(13)

where ${}_{B\hskip-0.85358pt}R_{W}\big{(}\stheta^{\mathrm{(W)}}_{n}\big{)}$ is the rotation matrix from the world to the body reference frame, i.e.,

{}_{B\hskip-0.85358pt}R_{W}\big{(}\stheta^{\mathrm{(W)}}_{n}\big{)}=\begin{bmatrix}\hskip 9.6739pt\cos{\big{(}\theta^{\mathrm{(W)}}_{n}\big{)}}&\sin{\big{(}\theta^{\mathrm{(W)}}_{n}\big{)}}\\ -\sin{\big{(}\theta^{\mathrm{(W)}}_{n}\big{)}}&\cos{\big{(}\theta^{\mathrm{(W)}}_{n}\big{)}}\end{bmatrix}\hskip-0.85358pt\hskip-0.85358pt\vskip 0.0pt

and $\bm{\xi}_{n}$ is a zero-mean multivariate Gaussian measurement noise with covariance matrix $\bm{I}_{2}\sigma_{\bm{\xi}}^{2}$ . The measurement noise is assumed to be statistically independent across $n$ .

Finally, the compass readings are modeled as

\displaystyle z_{\varphi_{n}}

\displaystyle=\theta_{n}+\eta_{n}

(14)

where $\eta_{n}$ is again a zero-mean Gaussian measurement noise with variance $\sigma^{2}_{\eta}$ . The measurement noise $\eta_{n}$ is also statistically independent across $n$ .

The measurement models in (12), (13), and (14) induce the following likelihood function at time $n$ , i.e.,

\displaystyle\mathpzc{p}\big{(}\bm{z}_{n}\big{|}\bm{x}_{n}\big{)}=\mathpzc{p}\big{(}z_{r_{n}}\big{|}\bm{p}_{n}^{\mathrm{(W)}}\big{)}\hskip 0.85358pt\mathpzc{p}\big{(}z_{\bm{v}_{n}^{\mathrm{(B)}}}\big{|}\bm{v}_{n}^{\mathrm{(W)}}\hskip-0.85358pt\hskip-0.85358pt,\theta_{n}^{\mathrm{(W)}}\big{)}\hskip 0.85358pt\mathpzc{p}\big{(}z_{\varphi_{n}}\big{|}\theta_{n}^{\mathrm{(W)}}\big{)}.

This likelihood function will be used in the update step of the proposed navigation filter.

III-A Estimation Problem and Navigation Filter

At time $n$ , we aim to estimate the state $\bm{x}_{n}$ of the AUV from all measurements $\bm{z}_{1:n}$ . Given the conditional PDF of the state, $p(\bm{x}_{n}|\bm{z}_{1:n})$ , the minimum mean-squared error (MMSE) estimate [25] of the state can be obtained as

\hat{\bm{x}}^{\mathrm{MMSE}}_{n}=\int\hskip-0.85358pt\bm{x}_{n}\hskip 0.85358pt\mathpzc{p}(\bm{x}_{n}|\bm{z}_{1:n})\hskip 0.85358pt\mathrm{d}\bm{x}_{n}.

(16)

A Bayes filter[26], which consists of prediction and update steps, is applied to compute the conditional PDF $p(\bm{x}_{n}|\bm{z}_{1:n})$ . The prediction step uses the Chapman-Kolomorogov equation, which involves the state-transition PDF. The update step is based on Baye’s rule and the likelihood function in (LABEL:eq:likelihood). Due to the nonlinearities in our measurement model, we follow the particle filtering approach [26] to compute a particle-based approximation of (16).

IV Results

The computed range rate and range results are presented based on the simulated and real data. Then, the navigation capability is demonstrated using the simulated data. Finally, the relevant hyperparameters to the simulation are summarized in Table. I.

Hyperparameters	Value
Prior Position STD (m)	$500$
Driving Noise STD, $\sigma_{\bm{u}_{1}}$ (m s^-1)	$0.02$
Turning Rate Noise STD, $\sigma_{\bm{u}_{2}}$ (° )	$0.50$
Measurement Noise (Heading) STD, $\sigma_{\eta}$ (° )	$0.50$
Measurement Noise (Range) STD, $\sigma_{\gamma}$ (m)	$150$
Measurement Noise (Velocity) STD, $\sigma_{\bm{\xi}}$ (m s^-1)	$0.02$
Average Acoustic SNR for navigation(dB)	$12$
Number of Particles	$500,\hskip-0.85358pt000$

TABLE I: Hyperparameters used to simulate the proposed AUV navigation method.

For the investigation using the real data, the acoustic recordings on the first element of the Horizontal Line Array South (HLAS) are processed from the S5 event of the SWellEx-96 experiment (Fig. 6(a)). The acoustic data from 23:52 to 00:29 in Greenwich Mean Time (GMT) are used. The sampling frequency is 3276.8 Hz, and the spectrogram was generated using a window length of 3276 samples (no overlap, Hanning window).

IV-A Simulation Environment and Data

As for the simulation, a scenario where a moving AUV, equipped with a single hydrophone, a DVL, and a compass, records a sound from a moving SOO in a shallow-water environment is considered. Initially, the SOO is at $(2000,0)$ moving north at 10 knots ( ${\sim}5.14$ m s^-1). The AUV starts at the origin, moving northeast at 3knots ( ${\sim}1.54$ m s^-1), i.e., $\bm{x}_{0}=[0,0,1.1,1.1,45]^{\mathrm{T}}$ , for 20 minutes. Its driving noise, $\sigma_{\bm{u}_{1}}$ , is 0.02 m s^-1, which accounts for random motion due to the control and the current in the ocean, and the turning rate noise, $\sigma_{\bm{u}_{2}}$ , is 0.5° .

The simulation environment follows the geoacoustic model from the S5 event of the SWellEx-96 experiment with the vertical line array (VLA)[21] (Fig. 4). The source depth was chosen as $z_{s}\hskip-0.85358pt=\hskip-0.85358pt9$ . This is the same depth as the towed source at the shallow depth in the SWellEx-96 experiment above the thermocline. Rouseff and Spindel [16] discuss the impact of different factors, including the source depth, on the “distribution” of the WI. The receiver depth was arbitrarily chosen close to the 34^th element of the VLA at a depth of $z_{r}\hskip-0.85358pt=\hskip-0.85358pt150$ m. The acoustic signals are simulated with an acoustic modal simulator, KRAKEN, which generates the pressure field in the range and frequency domain.

IV-B Range Rate and Range Computation

The range rate computation uses 4 min long segments, while the range computation uses a 10 min long signal. In the range rate computation, simulation and real data use the signals at $109$ Hz. For the range computation, the acoustic signals of four tones ( $109,127,145,163$ Hz) from a source at a shallow depth are used in both simulation and real data. When the first range rate measurement becomes available, another $10$ min long acoustic recording is necessary to compute the range. Recall that using the $4$ min long spectrogram, the first range rate can be computed after $2$ min. Hence, the first range measurement becomes available after $12$ min.

IV-B1 Simulation

Monte Carlo experiments with different signal-to-noise ratios (SNRs) of the acoustic recording under the scenario described in Sec. IV-A are performed to characterize the performance of the range rate computation and the range computation. The SNRs considered are $3$ , $5$ , $10$ , $15$ and $20$ dB, and $300$ Monte Carlo experiments were conducted for each SNR. Here, circularly symmetric Gaussian noise is added to the source signal.

The range rates computed from the simulated data are compared to the true range rate at every 30 s interval, and their root-mean-square errors (RMSEs) sampled for each SNR are calculated (Fig. 5(a) and (b)). The RMSE converges to a biased value of ${\sim}0.17$ m s^-1 as the SNR increases. This bias is introduced mainly by two factors. The first factor is the limited time interval sampling, which results in a range rate resolution of $0.11$ $\mathrm{m\,s^{-1}}$ . The second factor is a combination of the wavenumber variation [6] and the mismatch between the sound speed of the water, $v$ , and the horizontal phase speed, $\overline{v_{p}}$ .

Based on these computed range rates, the ranges between the source and receiver are calculated as described in Sec. II-B (Fig. 5(c)). The computed WI for the simulated environment is $z_{\beta,\mathrm{sim}}\hskip-0.85358pt=\hskip-0.85358pt1.11$ . The computed ranges are also compared to the true range at every 30 s interval to calculate the range RMSE over the Monte Carlo experiments for each SNR (Fig. 5(d)).

IV-B2 Real Data

The geoacoustic properties are equivalent to the simulation environment, except for the seafloor and receiver depths. The receiver is located on the seafloor, which is 198 m deep. The noise variance for the SNR calculation is obtained as the power at which the tones were not present.

The computed range rates are in good accordance with the true ranges rate but appear to be biased by $0.10$ m s^-1 (Fig. 6(b)). On top of the two factors for the bias discussed in Sec. IV-B1, there is also a model mismatch. The bathymetry along the event S5 slowly changes; therefore, it is a mildly range-dependent environment. The true range rate was computed by taking a moving average of the range difference between the receiver and the ship’s location, which was recorded with a GPS every minute. In addition, the computed range rate considers the phase offset discussed in [22]. In particular, given the sampling frequency of $3276.8$ Hz, the phase offset due to taking a segment of length $3276$ samples is $0.36$ m s^-1, which is removed from the computed range rates.

Using $z_{\beta,\mathrm{real}}\hskip-0.85358pt=\hskip-0.85358pt1.15$ , which was computed using the same data, the range RMSE compared to the actual range is $217$ m (Fig. 6(b)). The average SNR of the data used is calculated to be $8$ dB. Note also that the range computed using the average of the moving averaged true range rate ( $2.39$ m s^-1) yields an RMSE of $104$ m. In comparison, the range RMSE based on the simulation is approximately $109$ m at average SNR of $10$ dB in Fig. 5. Given the bias in the computed range rate, it is promising that the range computed using the SWellEx-96 event S5 data is close to the expected error.

Even though no real acoustic recordings from an AUV are available in this work, the following tracking simulation results show that our proposed method has a strong potential for bounding the error for self-localization in underwater navigation scenarios with just a single hydrophone in addition to the sensors for a dead-reckoning navigation.

IV-C Navigation

$300$ Monte Carlo experiments of the given the 20 min long scenario were performed to characterize the position error and covariance. The average SNR of the acoustic signal is $12$ dB. Assuming that the AUV has traveled based on dead-reckoning for some time, the initial position estimate is sampled from a normal distribution with a standard deviation of $500$ m. The corresponding WI parameter is computed as $z_{\beta,\mathrm{sim}}\hskip-0.85358pt=\hskip-0.85358pt1.11$ . To avoid particle degeneracy[27], $500,000$ particles were used, and a Gaussian kernel with a standard deviation (STD) of $10$ m was added to the position states of the particles at every resampling stage for roughening. Using a MacBook Pro with an Apple M1 Pro chip and $32$ GB memory, the particle-based estimation at each time step takes a median value of 0.11 s.

Each measurement noise STD reflects the error specifications of real sensors used in the AUV navigation. The STD of the heading and the speed over the ground are modeled after OceanServer OS5000[28] and Teledyne RD Instruments Explorer Doppler Velocity Log[29], respectively.

As discussed above, the first range measurement can be computed after $12$ min of acoustic data are available. As a result, sharp drops in the position error and the sample covariance of the particles after $12$ min are observed in Fig. 7. Furthermore, the navigation simulation results show that the position error is reduced once incorporated and keeps it bounded.

V Conclusion

We developed a Bayesian method for the navigation of autonomous underwater vehicles (AUVs) in shallow water. Our approach relies on passively recorded signals from acoustic sources of opportunities (SOOs). A measurement of the range with respect to the SOO is extracted from the spectrogram of a single hydrophone by exploiting striation patterns described by the waveguide invariant. A particle-based navigation filter processes the resulting range measurements and the AUV’s internal velocity and heading measurements. As a result, the position error, which would otherwise increase over time, can be bounded. The ability to compute the range rate and range measurements from the pressure field measured using a single hydrophone is demonstrated on real data from the SWellEx-96 experiment. The achieved ranging accuracy shows potential for real-world applications. In particular, by adding just a single hydrophone to a navigation system that only relies on dead-reckoning, the system can limit its position uncertainty if a suitable SOO is available. Future research directions include an extension to tracking of SOOs [30, 31, 32] and model adaptation using deep learning [33].

Acknowledgement

This work was supported by the Defense Advanced Research Projects Agency (DARPA) under Grant D22AP00151. We thank Dr. H. Akins, Dr. A. H. Young, Prof. S. T. Rakotonarivo, Prof. A. H. Harms and Prof. W. A. Kuperman for interesting discussions and comments.

References

[1] J. González-García, A. Gómez-Espinosa, E. Cuan-Urquizo, L. G. García-Valdovinos, T. Salgado-Jiménez, and J. A. E. Cabello, “Autonomous underwater vehicles: Localization, navigation, and communication for collaborative missions,” Appl. Sci., vol. 10, no. 4, 2020.
[2] R. P. Stokey and T. C. Austin, “Sequential long-baseline navigation for REMUS, an autonomous underwater vehicle,” in Proc. Information Systems for Navy Divers and Autonomous Underwater Vehicles Operating in Very Shallow Water and Surf Zone Regions, vol. 3711. SPIE, 1999, pp. 212–219.
[3] Y. Chen, D. Zheng, P. A. Miller, and J. A. Farrell, “Underwater inertial navigation with long baseline transceivers: A near-real-time approach,” IEEE Trans. Control Syst. Technol., vol. 24, no. 1, pp. 240–251, 2016.
[4] M. McKenna, D. Ross, S. Wiggins, and J. Hildebrand, “Underwater radiated noise from modern commercial ships,” J. Acoust. Soc. Am., vol. 131, pp. 92–103, 01 2012.
[5] H. M. Perez, R. Chang, R. Billings, and T. L. Kosub, “Automatic identification systems (AIS) data use in marine vessel emission estimation,” in Proc. 18th Annu. Int. Emission Inventory Conf., 2009.
[6] S. T. Rakotonarivo and W. A. Kuperman, “Model-independent range localization of a moving source in shallow water,” J. Acoust. Soc. Am., vol. 132, no. 4, pp. 2218–2223, 2012.
[7] F. B. Jensen, W. A. Kuperman, M. B. Porter, and H. Schmidt, Computational Ocean Acoustics, 2nd ed. Springer Publishing Company, Incorporated, 2011.
[8] S. Chuprov, “Interference structure of a sound field in a layered ocean,” Ocean Acoustics, Current State, pp. 71–99, 1982.
[9] G. L. D’Spain and W. A. Kuperman, “Application of waveguide invariants to analysis of spectrograms from shallow water environments that vary in range and azimuth,” J. Acoust. Soc. Am., vol. 106, no. 5, pp. 2454–2468, 1999.
[10] A. M. Thode, “Source ranging with minimal environmental information using a virtual receiver and waveguide invariant theory,” J. Acoust. Soc. Am., vol. 108, no. 4, pp. 1582–1594, 2000.
[11] K. L. Cockrell and H. Schmidt, “Robust passive range estimation using the waveguide invariant,” J. Acoust. Soc. Am., vol. 127, no. 5, pp. 2780–2789, 2010.
[12] A. Turgut, M. Orr, and D. Rouseff, “Broadband source localization using horizontal-beam acoustic intensity striations,” J. Acoust. Soc. Am., vol. 127, no. 1, pp. 73–83, 2010.
[13] C. Cho, H. C. Song, and W. S. Hodgkiss, “Robust source-range estimation using the array/waveguide invariant and a vertical array,” J. Acoust. Soc. Am., vol. 139, no. 1, pp. 63–69, 2016.
[14] C. Gervaise, B. G. Kinda, J. Bonnel, Y. Stéphan, and S. Vallez, “Passive geoacoustic inversion with a single hydrophone using broadband ship noise,” J. Acoust. Soc. Am., vol. 131, no. 3, pp. 1999–2010, 2012.
[15] B. Katsnelson, A. Lunkov, and I. Ostrovsky, “Interference pattern of the sound field in the presence of an internal kelvin wave in a stratified lake,” J. Acoust. Soc. Am., vol. 139, no. 2, pp. 881–889, 2016.
[16] D. Rouseff and R. Spindel, “Modeling the waveguide invariant as a distribution,” J. Acoust. Soc. Am., vol. 621, 06 2002.
[17] A. H. Young, J. Soli, and G. Hickman, “Self-localization technique for unmanned underwater vehicles using sources of opportunity and a single hydrophone,” in Proc. IEEE OCEANS, 2017, pp. 1–6.
[18] A. H. Young, H. A. Harms, G. W. Hickman, J. S. Rogers, and J. L. Krolik, “Waveguide-invariant-based ranging and receiver localization using tonal sources of opportunity,” IEEE J. Ocean. Eng., vol. 45, no. 2, pp. 631–644, 2020.
[19] X. Song, F. Liu, and Y. Shao, “Radial source velocity estimation using multiple line spectrum signals based on compressive sensing,” in Proc. IEEE WCSP, 2018.
[20] A. Zhao, P. Song, J. Hui, J. Li, and K. Tang, “Passive estimation of target velocity based on cross-spectrum histogram,” J. Acoust. Soc. Am., vol. 151, no. 5, pp. 2967–2974, 2022.
[21] J. Murray and D. Ensberg, “Swellex-96: S5 event,” http://swellex96.ucsd.edu/s5.htm, 1996, online; Accessed on 28-Feb-2023.
[22] C. Wang, J. Wang, and P. Du, “Estimation of moving target speed using weak line spectrum of single-hydrophon,” in Proc. IEEE Int. Conf. on Signal Process., Commun. and Comput., 2017, pp. 1–5.
[23] R. J. Urick, Principles of Underwater Sound, 3rd ed. Peninsula Publishing, 2013.
[24] A. Harms, J. L. Odom, and J. L. Krolik, “Ocean acoustic waveguide invariant parameter estimation using tonal noise sources,” in Proc. IEEE Int. Conf. Acoust., Speech, Signal Process., 2015, pp. 4001–4004.
[25] S. M. Kay, Fundamentals of Statistical Signal Processing: Estimation Theory. Upper Saddle River, NJ: Prentice-Hall, 1993.
[26] M. S. Arulampalam, S. Maskell, N. Gordon, and T. Clapp, “A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking,” IEEE Trans. Signal Process., vol. 50, no. 2, pp. 174–188, Feb. 2002.
[27] T. Li, M. Bolic, and P. M. Djuric, “Resampling methods for particle filtering: Classification, implementation, and strategies,” IEEE Signal Process. Mag., vol. 32, no. 3, pp. 70–86, 2015.
[28] OceanServer Technology, Inc., “Digital Compass Users Guide, OS5000 Series, Rev 4.0,” 2010.
[29] T. R. Instruments, “A Teledyne RD Instruments Navigation Datasheet,” Jul. 2013.
[30] F. Meyer, T. Kropfreiter, J. L. Williams, R. Lau, F. Hlawatsch, P. Braca, and M. Z. Win, “Message passing algorithms for scalable multitarget tracking,” Proc. IEEE, vol. 106, no. 2, pp. 221–259, 2018.
[31] J. Jang, F. Meyer, E. R. Snyder, S. M. Wiggins, S. Baumann-Pickering, and J. A. Hildebrand, “Bayesian detection and tracking of odontocetes in 3-D from their echolocation clicks,” J. Acoust. Soc. Am., vol. 153, no. 5, p. 2690, May 2023.
[32] F. Meyer and J. L. Williams, “Scalable detection and tracking of geometric extended objects,” IEEE Trans. Signal Process., vol. 69, pp. 6283–6298, 2021.
[33] M. Liang and F. Meyer, “Neural enhanced belief propagation for multiobject tracking,” IEEE Trans. Signal Process., 2023, to appear.

Navigation in Shallow Water Using Passive Acoustic Ranging