Estimation of seabed properties and range from vector acoustic observations of underwater ship noise

The Seabed Characterization Experiment (SBCEX) was conducted over the so-called New England mud patch, located about 100 km south of Cape Cod, MA, during March–April 2017, with the primary objective to study sound propagation in fine-grained, muddy sediments. Emerging results from SBCEX based on differing acoustic measurement methodologies¹ are in agreement on key sediment properties: a layer of fine-grained sediment of order 10 m thick starting at the water–sediment interface with an upper sediment sound speed nearly equal to or slightly less than the sound speed at the water–sediment interface,^2–6 with layer thickness diminishing on the periphery.⁷ These results are largely consistent with the extensive stratigraphic analysis of the area by Goff et al.⁸ using chirp sonar.

In this study, vector acoustic properties of underwater noise from a cargo ship traversing the central region of the SBCEX area are analyzed and inverted for sediment characteristics using a Bayesian approach, for which the model space is limited to the aforementioned layer thickness and sound speed, and an underlying basement sound speed. Observations were made with a vector sensor; vector field properties relating to the phase of pressure and vertical particle velocity are used as input for the geoacoustic inversion.

Recordings from the Intensity Vector Autonomous Recorder (IVAR) made on March 8 (20:53–21:53 UTC) captured the transit of the M/V Alice Oldendorf (IMO number: 9183776, length: 190 m). This IVAR deployment from the R/V Endeavor during leg 1 of SBCEX was at 40.4716°N, 70.6061°W, where the mean water depth during the transit was 75.5 m as measured by IVAR.

Acoustic particle velocity along three axes and pressure were recorded coherently in four channels, at 1.25 m above the seabed. IVAR system details are given in Dahl and Dall'Osto,⁹ which discusses measurements from explosive sound sources made at the second SBCEX deployment site, about 3.2 km northwest from this location. Ship noise recordings were made on the four high-gain (30 dB) channels of IVAR, while those involving explosive sources used the four low-gain (0 dB) channels.

The data are expressed in the frequency domain, formulated by the Fourier transform of each channel over identical 1-s data windows. These define, for example, $Sp(f;t)$ and $Svz(f;t)$ as the complex spectra for pressure and the vertical component of particle velocity, respectively. The potential energy spectrum is thus $(1/2ρcw2)⟨|Sp(f;t)|2⟩$ based on water density ρ of 1027.6 kg/m³ and mean water sound speed c_w of 1468.3 m/s, with angle brackets denoting a 10-s moving time average. Similarly, the kinetic energy spectrum is $(ρ/2)∑i=13⟨|Svi(f;t)|2⟩$⁠, where i represents the x, y, z components of particle velocity.

Spectrograms of the potential [Fig. 1(a)] and kinetic [Fig. 1(b)] energy density corresponding to the transit of the M/V Alice Oldendorf each illustrate an interference field pattern. The time-vector is aligned for t = 0 s to be at the closest-point-of-approach (CPA) to IVAR, ∼2470 m, based on the ship's Automatic Identification System (AIS). The received noise level increases slightly as the vessel opens in range after CPA, presumably owing to greater stern-aspect noise radiation. Note that the 1-h time average of kinetic and potential energy density is within 1 dB over the approximate bandwidth ∼25–200 Hz.

The active intensity bearing¹⁰ is derived from the pressure and the horizontal (x,y) component of particle velocity, and the direction of arrival [Fig. 1(c)] is 180° from the arc-tangent of these two orthogonal active intensity measures, rectified by the IVAR internal compass (x-channel oriented along bearing 199°). A robust, single estimate of direction-of-arrival versus time is obtained by a frequency average over the band 30–120 Hz.

An additional detail is provided by the cross-spectrum $⟨Sp(f;t)Svz⋆(f;t)⟩$⁠, where ^⋆ indicates complex conjugation. The phase pattern of this quantity [Fig. 1(d)] reveals cutoff frequencies for modes, and has been shown to be informative in vector acoustic observations of ship noise.¹¹ These data are used in the inversion process, where phase is given by the inverse tangent of the ratio of imaginary to real parts of $⟨Sp(f;t)Svz⋆(f;t)⟩$⁠. For the inversion, however, the absolute value of phase is used because phase magnitude is a more robust quantity that avoids rapid $±π$ swings, which are difficult to model precisely. In the following, a data vector $d$ is defined as the absolute value of the phase (normalized by π) at a given time t for frequencies 30–120 Hz every 1 Hz. Frequencies 86–88 Hz are eliminated from $d$ owing to electronic contamination, setting the length $d$ equal to 88.

The forward model used in the Bayesian inversion is based on the sum of normal modes for acoustic pressure and particle velocity. An analogous approach⁹ is developed for single modes applicable to cases where geometric dispersion is sufficient to identify separated modes. The complex spectral component of pressure P(f) received at range R and receiver depth z_r is expressed as

where U_n is the nth vertical mode function at frequency f, k_n is the corresponding horizontal wavenumber (computed using the Kraken normal mode program¹²), z_s is source depth, z_r is receiver depth (fixed to the IVAR sensor depth of 1.25 m above the seabed), and the asymptotic representation for the Hankel function is used where R is range. An unknown source spectrum $S0(f)$ is shown for completeness, but is ultimately eliminated from the model upon taking a ratio.

Horizontal and vertical derivatives of Eq. (1) yield the corresponding horizontal $Vr(f)$ and vertical $Vz(f)$ particle velocities, respectively. For $Vr(f)$ the mode sum includes the factor $(kn/2πfρ)$⁠, where a second term in $Vr(f)$ associated with the horizontal derivative becomes negligible compared with the first term for ranges greater than a few wavelengths and is thus ignored. For $Vz(f)$ it is $Un(zr,f)$ in the mode sum that is replaced by a finite difference approximation for the vertical derivative evaluated at z_r and multiplied by $(i/2πfρ)$⁠.

From these quantities the potential E_p and kinetic E_k energy density are formed as

and complex vertical intensity as

from which the real part yields active vertical intensity $Iz(f)$ and the imaginary part reactive vertical intensity $Qz(f)$⁠.

There are multiple options to proceed toward a modeled quantity that is independent of the source spectrum. One is via the ratio of $Iz(f)$ to $Ep(f)+Ek(f)$⁠, which was developed as an inversion metric for observations that support resolved single mode arrivals.⁹ Another, used recently for inversion from an active source,¹³ is the vertical phase gradient of acoustic pressure, which is equivalent to the ratio $Iz(f)$ to $Ep(f)$ multiplied by $(2πf/4cw2)$⁠. The approach used here takes the absolute value of the phase of $Izc$ normalized by π, which corresponds to the data vector. This is computed for the same 88 frequencies as $d$⁠, based on a given geoacoustic parameter vector $m$⁠, and defines the model vector $d(m)$ representing a simulated replica of the data.

Geoacoustic inversion is undertaken for which the unknown parameters are fine sediment layer depth H_sed, a constant sound speed within this layer c_sed, and sound speed of an underlying half-space c_base. A combination of uniform prior bounds for these parameters (Table 1) and equidistant grid search space in each model dimension is used in the forward modeling. Additional required parameters were pre-set to streamline the inversion because experimentation suggests these parameters are poorly resolved. These are sediment density set to 1.6 g/cm³, to represent a 1.5–1.7 g/cm³ range reported in gravity coring operations from SBCEX;¹⁴ sediment attenuation based on a fit to viscous grain shearing theory¹⁵ applicable to the inversion frequency range 30–120 Hz; a basement density that depends on c_base following Hamilton,¹⁶ translating to 1.7–2.1 g/cm³ for the range of c_base values used, and basement attenuation set to 0.6 dB/λ. Constant water depth and sound speed (c_w) were taken as 75.5 m and 1468.3 m/s, respectively.

Additionally, range and source depth are unknown and this search space is assessed with little computational cost using the mode-based forward model. Source depth is uniformly sampled between 1 and 20 m in 1-m increments, and range is sampled between 1000 and 4000 m in 5-m increments. (This range extent can be reduced to ∼200 m upon invoking prior information on range provided by AIS data, but our interest is in keeping range as an unknown parameter.)

The inversion is placed within a Bayesian framework¹⁷ leading to the probability density function (PDF) $P(m|d)$⁠, representing the conditional PDF for model parameters $m$ given the experimental data $d$⁠.¹⁸ This equals the conditional PDF for the data given the model $P(d|m)$ multiplied by an a priori PDF for the model parameters $P(m)$⁠, assumed in this study to be a uniform distribution within the selected parameter bounds.¹⁹ Key steps then follow closely Eqs. (15)–(17) from Bonnel et al.;²⁰ here, put the model–data mismatch column vector $r=d−d(m)$⁠, and find the data misfit function $E(m)=(Nf/2) logerTr$⁠, with N_f = 88, which, given the established $P(m)$ leads directly to $P(m|d)$⁠. One-dimensional (1D) marginal PDFs for the geoacoustic model parameters (including the range and source depth), are found upon integrating $P(m|d)$ over the appropriate model dimensions.

Figure 2(a) presents estimated values (circles) of the phase magnitude of the cross-spectrum between pressure and vertical particle velocity at CPA compared with a model (Sec. 3.1) based on the maximum a posteriori (MAP) estimates²⁰ for H_sed (9.4 m), c_sed (1499 m/s), and c_base (1825 m/s). The MAP estimate for range is 2380 m and for source depth is 3 m. We note that the model–data errors in Fig. 2(a) are approximately uniform across the inversion frequency band, and can be approximated by a Gaussian distribution.

We next invert data from 21 ranges spaced equally over a 20-min period centered about the CPA, with each inversion providing a new $P(m|d)$ and set of MAP estimates. Each $P(m|d)$ yields 1D marginal PDFs for the inverted parameters. To evaluate uncertainty we derive mixture²¹ 1D PDFs [Figs. 2(b)–2(d)], based on a weighted linear combination of the 21 1D PDFs (representing component densities) with uniform mixing proportion equal to 1/21. The median values of the 21 MAP estimates, and the standard deviations based on the mixture 1D PDFs, constitute our final inversion estimates and uncertainty, respectfully (Table 1).

Finally, it is of interest to invert just for range (including source depth) using the final inversion estimates (Table 1) as known geoacoustic parameters. Here, a longer transit period is studied for which the range between the M/V Alice Oldendorf and IVAR varies from the CPA at ∼2.4 km to ∼7.5 km. Source depth is again uniformly sampled between 1 and 20 m, and range is now sampled between 1000 and 8000 m. The direction of arrival associated with each 10-s data vector [Fig. 1(c)] is combined with the inversions for range to construct the ship trajectory. The results (Fig. 3), plotted over a map of two-way travel time (twtt) of mud layer thickness,⁸ show correspondence between inverted range and computed bearing, and the AIS data. The 20-min period from which the geoacoustic parameters are estimated is indicated by the 21 triangle symbols, and here propagation paths between the ship and IVAR tend to traverse a relatively constant twtt regime of 12–13 ms. Interestingly, the propagation paths associated with the far west and east portions of the longer transit period appear to interact, respectively, with sediment characterized by a somewhat thinner and thicker mud layer, motivating a future study on seabed variation. The inversion necessarily yields an estimate of source depth with MAP estimates more variable but consistently confined to less than 10 m.

This paper applies a Bayesian framework for geoacoustic inversion of the relative phase of pressure and vertical particle velocity of ship noise, a quantity independent of noise source spectrum, observed during the SBCEX conducted about 100 km south of Cape Cod in waters of depth ∼75 m. A layer of fine-grained, muddy sediment is a distinctive feature of the SBCEX experimental area, and the focus of this study. The results reported here are in agreement with emerging results reported by SBCEX investigators, in particular the sediment layer thickness of ∼9.4 m and sound speed of ∼1478 m/s within the layer. A sound speed gradient could not be resolved in preliminary studies and thus was not included in the model space, but a corresponding sediment layer interval velocity of 1.007 (based on c_w) does not rule out a gradient. In addition, the inversion provides credible estimates of ship range and bearing (available through inspection of the angle of horizontal active intensity), both verified by AIS data.

This research was supported by the U.S. Office of Naval Research.