The Intensity Vector Autonomous Recorder (IVAR) simultaneously measures acoustic particle velocity and pressure. IVAR was deployed during the 2017 Seabed Characterization Experiment (SBCEX) with the primary objective to study sound propagation in fine-grained, muddy sediments. In this study a Bayesian inversion framework is applied to ship underwater noise recorded by IVAR. The data are relative phase of pressure and vertical particle velocity, a quantity that is independent of the ship noise source spectrum. Inversion estimates for the sediment layer and underlying basement properties are in agreement with other reports from SBCEX.

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 methodologies1 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.7 These results are largely consistent with the extensive stratigraphic analysis of the area by Goff et al.8 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,9 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/m3 and mean water sound speed cw 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.

Fig. 1.

(Color online) (a) Potential and (b) kinetic energy density recorded during the transit of the M/V Alice Oldendorf. (c) Active intensity directional of arrival. (d) Phase cross-spectrum between pressure and vertical particle velocity revealing approximate cutoff frequencies for modes as suggested by horizontal dashed lines, e.g., mode 3 emerges at frequency 35 Hz. Time for (a)–(d) is relative to the CPA.

Fig. 1.

(Color online) (a) Potential and (b) kinetic energy density recorded during the transit of the M/V Alice Oldendorf. (c) Active intensity directional of arrival. (d) Phase cross-spectrum between pressure and vertical particle velocity revealing approximate cutoff frequencies for modes as suggested by horizontal dashed lines, e.g., mode 3 emerges at frequency 35 Hz. Time for (a)–(d) is relative to the CPA.

Close modal

The active intensity bearing10 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.11 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 approach9 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 zr is expressed as

P(f)S0(f)neiknRi(π/4)knRUn(zs,f)Un(zr,f),
(1)

where Un is the nth vertical mode function at frequency f, kn is the corresponding horizontal wavenumber (computed using the Kraken normal mode program12), zs is source depth, zr 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 zr and multiplied by (i/2πfρ).

From these quantities the potential Ep and kinetic Ek energy density are formed as

Ep(f)=14|P(f)|2/(ρcw2),
(2)
Ek(f)=14ρ(|Vr(f)|2+|Vz(f)|2),
(3)

and complex vertical intensity as

Izc=12P(f)Vz(f),
(4)

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.9 Another, used recently for inversion from an active source,13 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 Hsed, a constant sound speed within this layer csed, and sound speed of an underlying half-space cbase. 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/cm3, to represent a 1.5–1.7 g/cm3 range reported in gravity coring operations from SBCEX;14 sediment attenuation based on a fit to viscous grain shearing theory15 applicable to the inversion frequency range 30–120 Hz; a basement density that depends on cbase following Hamilton,16 translating to 1.7–2.1 g/cm3 for the range of cbase values used, and basement attenuation set to 0.6 dB/λ. Constant water depth and sound speed (cw) were taken as 75.5 m and 1468.3 m/s, respectively.

Table 1.

Inversion geoacoustic parameters list, uniform over given search bounds, and final inversion estimates. Final estimates represent the median MAP over the 21 different range inversions about the CPA and uncertainty σ is based on mixture PDFs as discussed in text. A basement density of 1.9 g/cm3 corresponds to the estimate cbase= 1800 m/s.

ParameterUnitSearch boundsMedian MAPσ
csed m/s [0.985,1.05]cw 1.007cw 0.009cw 
Hsed [7.0,11.5] 9.4 0.27 
cbase m/s [1700,2025] 1800 17 
ParameterUnitSearch boundsMedian MAPσ
csed m/s [0.985,1.05]cw 1.007cw 0.009cw 
Hsed [7.0,11.5] 9.4 0.27 
cbase m/s [1700,2025] 1800 17 

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 framework17 leading to the probability density function (PDF) P(m|d), representing the conditional PDF for model parameters m given the experimental data d.18 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.19 Key steps then follow closely Eqs. (15)–(17) from Bonnel et al.;20 here, put the model–data mismatch column vector r=dd(m), and find the data misfit function E(m)=(Nf/2)logerTr, with Nf = 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) estimates20 for Hsed (9.4 m), csed (1499 m/s), and cbase (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.

Fig. 2.

(a) Phase magnitude of the cross-spectrum (normalized by π) between pressure and vertical particle velocity at CPA (black circles) compared with modeled (gray line) MAP estimates for Hsed (9.4 m), csed (1499 m/s), and cbase (1825 m/s), along the MAP estimate for range 2380 m, and source depth 3 m. (b) Mixture 1D PDF for sediment thickness Hsed based on inversion of 21 ranges equally spaced over the 20-min observation period and weighted linear combination of the 21, 1D PDFs; listed are the average of the 21 MAP estimates and standard deviation computed from the mixture 1D PDF. (c) Same display for layer speed, csed. (d) Same display for cbase.

Fig. 2.

(a) Phase magnitude of the cross-spectrum (normalized by π) between pressure and vertical particle velocity at CPA (black circles) compared with modeled (gray line) MAP estimates for Hsed (9.4 m), csed (1499 m/s), and cbase (1825 m/s), along the MAP estimate for range 2380 m, and source depth 3 m. (b) Mixture 1D PDF for sediment thickness Hsed based on inversion of 21 ranges equally spaced over the 20-min observation period and weighted linear combination of the 21, 1D PDFs; listed are the average of the 21 MAP estimates and standard deviation computed from the mixture 1D PDF. (c) Same display for layer speed, csed. (d) Same display for cbase.

Close modal

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 mixture21 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,8 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.

Fig. 3.

(Color online) The ship track of the M/V Alice Oldendorff as it transits from east to west (dashed line is ground truth from AIS) and estimated from IVAR vector intensity data (green circles). The IVAR estimates include range via inversion, and bearing via inspection of active intensity, as per Fig. 1(c). Geoacoustic estimates are made at 21 locations (triangle symbols) to establish the mean environment (Table 1) for range estimation. Background color corresponds to the twtt measured during a chirp sonar survey of the area (Ref. 8), and inset shows detail of AIS and MAP range estimates over the inversion regime. Background data as plotted here provided courtesy John Goff.

Fig. 3.

(Color online) The ship track of the M/V Alice Oldendorff as it transits from east to west (dashed line is ground truth from AIS) and estimated from IVAR vector intensity data (green circles). The IVAR estimates include range via inversion, and bearing via inspection of active intensity, as per Fig. 1(c). Geoacoustic estimates are made at 21 locations (triangle symbols) to establish the mean environment (Table 1) for range estimation. Background color corresponds to the twtt measured during a chirp sonar survey of the area (Ref. 8), and inset shows detail of AIS and MAP range estimates over the inversion regime. Background data as plotted here provided courtesy John Goff.

Close modal

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 cw) 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.

1.
P. S.
Wilson
,
D. P.
Knobles
, and
T. B.
Neilsen
, “
Guest editorial: An overview of the seabed characterization experiment
,”
IEEE J. Oceanic Eng.
45
,
1
13
(
2020
).
2.
J.
Bonnel
,
S. E.
Dosso
,
D.
Eleftherakis
, and
N. R.
Chapman
, “
Trans-dimensional inversion of modal dispersion data on the New England mud patch
,”
IEEE J. Oceanic Eng.
45
,
116
130
(
2020
).
3.
L.
Wan
,
M.
Badiey
,
D. P.
Knobles
,
P. S.
Wilson
, and
J. A.
Goff
, “
Estimates of low-frequency sound speed and attenuation in a surface mud layer using low-order modes
,”
IEEE J. Oceanic Eng.
45
,
201
211
(
2020
).
4.
J.
Belcourt
,
C. W.
Holland
,
S. E.
Dosso
,
J.
Dettmer
, and
J. A.
Goff
, “
Depth-dependent geoacoustic inferences with dispersion at the New England mud patch via reflection coefficient inversion
,”
IEEE J. Oceanic Eng.
45
,
69
91
(
2019
).
5.
Y.-T.
Lin
,
J.
Bonnel
,
D. P.
Knobles
, and
P. S.
Wilson
, “
Broadband waveform geoacoustic inversions with absolute travel time
,”
IEEE J. Oceanic Eng.
45
,
189
200
(
2019
).
6.
D. P.
Knobles
,
P. S.
Wilson
,
J. A.
Goff
,
L.
Wan
,
M. J.
Buckingham
,
J. D.
Chaytor
, and
M.
Badiey
, “
Maximum entropy derived statistics of sound-speed structure in a fine-grained sediment inferred from sparse broadband acoustic measurements on the New England continental shelf
,”
IEEE J. Oceanic Eng.
45
,
161
173
(
2020
).
7.
D.
Tollefsen
,
S. E.
Dosso
, and
D.
Knobles
, “
Ship-of-opportunity noise inversions for geoacoustic profiles of a layered mud-sand seabed
,”
IEEE J. Oceanic Eng.
45
,
189
200
(
2019
).
8.
J. A.
Goff
,
A. H.
Reed
,
G.
Gawarkiewicz
,
P. W.
Wilson
, and
D. P.
Knobles
, “
Stratigraphic analysis of a sediment pond within the New England mud patch: New constraints from high resolution chirp acoustic reflection data
,”
Marine Geol.
412
,
81
94
(
2019
).
9.
P. H.
Dahl
and
D. R.
Dall'Osto
, “
Vector acoustic analysis of time-separated modal arrivals from explosive sound sources during the 2017 seabed characterization experiment
,”
IEEE J. Oceanic Eng.
45
,
131
143
(
2020
).
10.
D. R.
Dall'Osto
and
P. H.
Dahl
, “
Observations of water column and bathymetric effects on the incident acoustic field associated with shallow-water reverberation experiments
,”
IEEE J. Oceanic Eng.
42
,
1146
1161
(
2017
).
11.
V. A.
Shchurov
,
V. P.
Kuleshov
, and
A. V.
Cherkasov
, “
Vortex properties of the acoustic intensity vector in a shallow sea
,”
Acoust. Phys.
57
,
851
856
(
2011
).
12.
M. B.
Porter
and
E. L.
Reiss
, “
A numerical method for ocean acoustic normal modes
,”
J. Acoust. Soc. Am.
76
,
244
252
(
1984
).
13.
J.
Shi
,
S. E.
Dosso
,
D.
Sun
, and
Q.
Liu
, “
Geoacoustic inversion of the acoustic-pressure vertical phase gradient from a single vector sensor
,”
J. Acoust. Soc. Am.
146
,
3159
3173
(
2019
).
14.
M. S.
Ballard
,
K. M.
Lee
,
A. R.
McNeese
,
P. S.
Wilson
,
J. D.
Chaytor
,
J. A.
Goff
, and
A. H.
Reed
, “
In situ measurements of compressional wave speed during gravity coring operations in the New England mud patch
,”
IEEE J. Oceanic Eng.
45
,
26
38
(
2019
).
15.
M. J.
Buckingham
, “
On pore-fluid viscosity and the wave properties of saturated granular materials including marine sediments
,”
J. Acoust. Soc. Am.
122
,
1486
1501
(
2007
).
16.
E. L.
Hamilton
, “
Geoacoustic modeling of the sea floor
,”
J. Acoust. Soc. Am.
68
,
1313
1340
(
1980
).
17.
R.
Duan
,
N. R.
Chapman
,
K.
Yang
, and
Y.
Ma
, “
Sequential inversion of model data for sound attenuation in sediment at the new jersey shelf
,”
J. Acoust. Soc. Am.
139
,
70
84
(
2016
).
18.
N. R.
Chapman
, “
Perspectives on geoacoustic inversion of ocean bottom reflectivity data
,”
J. Mar. Sci. Eng.
4
,
61
(
2016
).
19.
C.-F.
Huang
,
P.
Gerstoft
, and
W. S.
Hodgkiss
, “
Uncertainty analysis in matched-field geoacoustic inversions
,”
J. Acoust. Soc. Am.
119
,
197
207
(
2006
).
20.
J.
Bonnel
,
S. E.
Dosso
, and
N. R.
Chapman
, “
Bayesian geoacoustic inversion of single hydrophone light bulb data using warping dispersion analysis
,”
J. Acoust. Soc. Am.
134
,
120
130
(
2013
).
21.
G. J.
McLachlan
,
S. X.
Lee
, and
S. I.
Rathnayake
, “
Finite mixture models
,”
Annu. Rev. Stat. Appl.
6
,
355
378
(
2019
).