Inference of source signatures of merchant ships in shallow ocean environments^a)

The modeling and measurement of the frequency-dependent source level, SL(f), of merchant ships (Gray and Greeley, 1980; McKenna , 2012; Gassmann , 2017) continue to be of importance because merchant shipping represents the largest source of anthropogenic noise in the ocean. While the inference of SL(f) in deep ocean environments is relatively straightforward, to accomplish this in shallow water remains one the most challenging problems in ocean acoustics (Crocker , 2014; Jiang , 2020; Tollefsen , 2022; MacGillivray , 2023). The ability to characterize the sound propagation in bottom-limited regions is generally more complicated in part because the frequency-dependent propagation loss [PL(f)] becomes more sensitive to the geoacoustic properties of seabed. This in turn complicates the inference of SL(f), which requires knowledge of the PL(f). First, noise, model errors, and parameter ambiguities of a geoacoustic model limit the range accuracy of modeled PL(f). Second, because of the random temporal variability and aspect dependence of SL(f), any inference methodology needs to average SL(f) over some minimal time/range interval of observation. No simple approach exists to find an optimal solution to what we refer to as the shallow-water range dilemma.

The proposed solution to the range dilemma is a two-step process that utilizes a prior geophysical model, which assumes a two-layered seabed sediment parametrization of the viscous grain shearing (VGS) model (Buckingham, 2000, 2007) over a fixed deep seabed sediment. First, for each random sampling of the source depth, the closest position of approach (CPA) range, and the VGS hypothesis space $θ$⁠, the SL, SL $( θ , f , z j )$⁠, for $( j = 1 , 2 , … , N r − 1 , N r$⁠), is estimated implicitly (Dosso and Wilmut, 2006; Knobles, 2015), where z_j is the receiver depth of the jth hydrophone of the vertical line array (VLA), and N_r is the number of receivers. This SL is an average over the full time interval of a time-symmetric spectrogram. Second, the inversion method exploits the broadband coherent structure of the acoustic intensity at the CPA time. The main idea is the acoustic field at CPA over a large enough frequency band may possess sufficient information content about the properties of the source-receiver geometry and the acoustical properties of the seabed to permit a statistical inference of SL $( f , z j )$⁠. This concept may be viewed as a variant of a recently reported idea called feature-based inversion (Knobles , 2021; Knobles , 2022). The estimated time-averaged SL $( θ , z j , f k )$ is used to evaluate the error function at the CPA range (corresponding to $t i = t cpa$⁠). The statistics of the CPA-based error function are then utilized to estimate sediment parameters of the VGS model, the source depth, the CPA range, and the SL.

As a test of the efficacy of the proposed methodology, this study takes advantage of the placement of VLAs near a shipping lane in the New England Mudpatch during the Seabed Characterization Experiment (SBCEXP) in 2017 (Wilson , 2020). The measurements were made during benign water column conditions. Data samples were processed in the 200–750 Hz band, where the sound radiation of the ships is mostly characterized by continuous broadband noise generated by cavitation at the rotating propeller. Estimates of seabed and source parameters are obtained from optimal solutions of the parameter space and the properties of parameter marginal probability density functions (PDFs) derived from a maximum entropy approach. The PL(f) values for these solutions in the parameter space are used to estimate SL $( θ ̂ , f )$⁠, where $θ ̂$ is the parameter vector that minimizes an error function. The estimated SL $( θ ̂ , f )$ values for multiple ships are then compared to values for $SL W − H ( f )$ computed with the empirical Wales–Heitmeyer SL model (Wales and Heitmeyer, 2002).

The remainder of this paper is organized as follows. Section II discusses the acoustical measurements. Section III discusses prior assumptions about the ship motion and geophysical/geoacoustic characterization of the seabed and empirical constraints on the seabed model. Section IV defines the error function, inferred SLs, and PDFs. Sections V and VI provide the results and a conclusion, respectively.

Merchant ship noise data were recorded on two VLAs deployed on the New England continental shelf during a time when the sound speed of the water column was approximately constant with depth. The locations of the two VLAs, the shipping lane, and the track of the Viking Bravery that deviated from the southern shipping lane are shown in Fig. 1(a). The start and end times of the processing intervals of merchant ships considered for analysis, shown in Table I, place the CPA time at about the midpoint time of the processed data.

For this analysis, the spectrograms are processed in time intervals of about 20 min over the 200–750 Hz band, with a frequency spacing of 0.3815 Hz. An example of a processed spectrogram obtained from channel 05 on VLA 02 is shown in Fig. 1(b) for the merchant ship Kalamata. The peak and null structure at CPA time is a result of the coherent combination of direct, bottom, and surface reflected components.

An additional data sample from the ship Viking Bravery was collected at VLA 01 about 20 min after the recording was made on VLA 02. This data sample was held in reserve to test how well the implicit SL obtained from the data processing on VLA 02 generalized to data recorded on VLA 01. Only in the case of the Viking Bravery is a ship track reasonably aligned with the main sediment thickness map. The other tracks are skewed, and, thus, one would expect a greater lack of symmetry in the spectrograms for a single ship recorded at the two VLAs.

In this work, two upper sediment layers are parameterized by a hybrid VGS model, as shown in Fig. 3. Below these layers is a fixed deep-sediment layering that was inferred in Knobles et al. (Knobles , 2021; Knobles , 2022) to model the very low frequency characteristics of the waveguide invariant. The hybrid VGS model uses the VGS parameterization to establish the sound speed, the density, and the attenuation at the surface of each sediment layer and a non-linear depth dependence (NLDD) parameterization of the sediment sound speed as a function of depth in each layer. The NLDD parameterization was formulated to characterize the empirical nature of observed low-angle bottom loss for different classes of seabed environments (Spofford, 1980).

While K_g, ρ_g, n, and τ do not have analytical relationship to porosity, their empirical relationship to porosity can be estimated. The Wentworth scale (Wentworth, 1922) classifies sediments in terms of grain size or, equivalently, porosity. A Wentworth-inspired VGS sediment classification based on porosity is shown in Table II. The values for $K g ( N )$ and $ρ g ( N )$ in the different sediment classes are approximately established from the literature with the understanding that the values for K_g can have significant uncertainty. The values of n(N) and log[ $τ ( N )$] are deduced from collective attenuation measurements in the literature. For example, modeled values (not shown here) were compared to previously reported graphics of compressional sound speed ratio and attenuation of a medium sand (Zhou , 2009; Zhou , 2010; Richardson , 2001; Carey and Evans, 1988) and a clayey sandy silt (Wilson , 2020; Boyles, 1997). The modeled values were obtained from the VGS model, with the X(N) relationships shown in Table II for a medium sand and clayey sandy silt. In addition to these examples, reported sound speed and attenuation data for clay from measurements in the Gulf of Mexico (Rubano, 1980; Lynch , 1991; Collins , 1992; Knobles , 2003), a coarse silt sediment in the New England Bight (Potty , 2003), a thick fine silt sediment in the Gulf of Oman (Sagers and Knobles, 2014), and a mixed sediment region of the Yellow Sea (Rogers , 2000) were used to establish n(N) and $log 10 ( τ$⁠) for different sediment classifications in Table II.

In Table II, $log 10 ( τ )$ is characterized by the two asymptotic values for a very soft clay and sand granules of 3.05 and −3.92, respectively, with a sudden change in values that occurs between silty sand and clayey sand-silt. This change is a reflection of not having enough sound speed and attenuation dispersion data for sediment classes containing various mixtures of sands and silts to estimate empirically the width on a classification scale of the transition from the VGS to the grain shearing (GS) limit. Future work is needed to better understand the nature of viscous effects in mixed sediments. While most of the results presented use K_g, n, and τ values that are related to N, as specified in Table II, in Sec. V B, τ, K_g, and n in the top sediment layer are varied randomly, independent of N, as a form of error analysis.

The concept of Table II is not new or unique. Hamilton and Bachman (Hamilton and Bachman, 1982; Bachman, 1985; Hamilton, 1979, 1980), Richardson and Briggs (1993), and Vidmar (1983) have all used a similar approach in sediment classification. The most important aspect in the current analysis that is different from the previous attempts at sediment classification is the use of a causal model for the seabed that restricts the frequency dispersion of the sound speed and attenuation to a Kramers–Kronig relationship (Toll, 1956; Jackson, 1975). The causality constraint acts as an important means to reduce the dimensionality of the problem and constrain the multi-dimensional volume of the parameter space in a physical manner.

Advantages of using the VGS model with sediment classes based on porosity include (1) a reduced-order model, (2) physically realistic sediments, and (3) an ability to make comparison with measured geophysical ground truth. The VGS model that includes viscoelastic effects of a linear response is based on a simple mechanical system with a Zener model (Carcione, 2014). Since the VGS model is a causal model, the attenuation is not a free parameter; this fact allows for the general ambiguity of SLs and attenuation to be mitigated with data samples corresponding to longer ranges. While there are other physical models, such as poro-elastic models, that can also provide adequate fits to the dispersion characteristics of, for example, the sediment sound speed, such models also have a greater number of independent parameters. For example, permeability, creep, and tortuosity are important parameters for poro-elastic models (Chotiros, 2017, 2021), but how to relate such parameters in a scheme that orders classification based on porosity, such as in the current approach, remains unknown.

The inversion processing methodology, shown in Fig. 4, is initiated with a Monte Carlo sampling of the parameter space defined in Tables II and III. For each sampling of the parameter space $θ$⁠, the hybrid VGS model and NLDD model are used to compute the inputs for a normal mode propagation model (Westwood , 1996), which are used to calculate frequency-dependent modeled propagation loss, $PL M ( f , θ , z j )$(⁠ $j = 1 , 2 , … , N r$⁠). Then, using the data $D ℓ$⁠, a time-averaged SL, $SL ℓ ( θ , f , z j )$ is computed, which then permits a value of the modeled received level to be computed at the CPA time for each sampling of $θ$⁠.

The statistical optimization approach begins with defining the prior distribution. For this work, the priors include (1) identifying porosity as the independent variable and assuming that the other geophysical parameters with the VGS model are dependent parameters X(N), (2) using two sediment layers, and (3) establishing the upper and lower parameter bounds in Table III. The parameter vector $θ$ contains $[ N 1 , N 2 , T 1 , T 2 , z s , R cpa$], where N₁ and T₁ are the porosity and thickness of the first sediment layer, N₂ and T₂ are the porosity and thickness of the second sediment layer, z_s is an effective depth for the source, and R_cpa is the CPA range. The volume of $θ$ is defined by the bounds shown in Table II: the prior distributions are assumed to be uniform between these upper and lower bounds.

The results of the maximum entropy inversions are now presented for the five ships listed in Table I. Table IV shows the elements of the 5 × 5 matrix $E ( θ ̂ geo ( D ℓ ) , θ ̂ source ( D k ) ) , ⟨ E ℓ ⟩$⁠, and $b ℓ$ for all five ships. Then, using the PPD in Eq. (15), the PDFs for the six parameters in $θ$ are obtained for each data sample. The estimated PDFs are shown in Fig. 5 for $ℓ = 1 , 2 , 3 , 4 , 5$⁠.

Both Fig. 5 and the values in Tables V–IX provide insights into the parameters estimated by the maximum entropy method. Generally, the parameter values for N₁ and CPA range are well resolved. Parameter values for N₂ are not as well resolved as are those for N₁. The resolution for T₂ is generally poor. With the exception of the Matsuyama data sample, which has the largest CPA range, the $P ( θ peak )$ values for the porosity N₁ are near the ground truth value of 0.6012 from an analysis of piston cores taken by United States Geological Survey (USGS) (Chaytor , 2021). The marginal distributions for $R cpa$ show that $R cpa$ and the uncertainty σ are inversely related: The ships Hafnia Green and Matsuyama have the smallest and largest CPA ranges and range uncertainty, respectively. We may ascribe this to an increase in the uncertainty of the localization with increasing range. Finally, the average optimal values of z_s and σ for all the ships are 6.14 m and a 1.2 m, respectively. The small value of σ suggests that the Lloyd's mirror effect was being adequately modeled.

To further test the inferred parameter values, comparisons are shown in Figs. 6 and 7 and Figs. 8 and 9 of $PL mea$ and $PL M$ at the CPA times for all 16 hydrophone elements for the Kalamata and Tombarra data samples, respectively. These two data samples are separated in time by about 1 week. However, the CPA range and source speeds for the two ships are quite similar. The $PL M$ and $PL mea$ comparisons for the Kalamata and the Tombarra data samples are evaluated at the optimal values $θ ̂ ℓ = 1$ and $θ ̂ ℓ = 2$⁠, respectively. The model and data PL comparisons for $ℓ = 3 , 4 , 5$ (not shown here) are similar.

The $PL M$ captures the envelope of the peak and null structure of $PL mea$ but is unable to capture the observed highly irregular hash that is observed to be superimposed on the envelope. The optimal values of the $E ℓ$ occurs for values of N₂ of about 0.55 for both data samples. This value of N₂ is high enough to reduce the cause of the perturbations of the envelope structure by reducing the magnitude of the reflection coefficient at the first and second layer interface. Lower values of N₂ result in an irregular hash structure similar to the observed PL_mea; however, it is at the expense of an increase in the error function away from $E ( θ ̂ )$⁠.

The measured-model comparison of time-frequency propagation loss (PL) for channel 08 of the Kalamata and the Tombarra are shown in Figs. 10 and 11, respectively, over time windows of about 20 min. The $PL mea$ is computed using Eq. (12). Qualitatively the agreement for both data samples is good. In summary, the $PL M$ captures the envelope structure of $PL mea$⁠, which corresponds to the lower grazing angles. Such a model may not be adequate at shorter ranges or equivalently at high grazing angles.

In Sec. III B, K_g, $log 10 ( τ )$⁠, and n are considered fully dependent parameters of N and held fixed at discrete values of N (as listed in Table II) that essentially followed the Wentworth scale. Another series of optimizations was performed in which K_g, $log 10 ( τ )$⁠, and n in the top layer are varied independent of N with Monte Carlo sampling. The resulting nine-dimensional (9-D) $θ$ values, which include upper and lower parameter bounds for K_g, $log 10 ( τ )$⁠, and n, are shown in Table X. A statistical inversion was made for the 9-D $θ$ space, and the resulting expected and standard deviation of $θ ℓ$ for the Tombarra data sample are shown in Table XI. A comparison of Table XI to Table VI, which is for the six-dimensional (6-D) space, provides the ratios of σ from the PDFs for the 9-D to the 6-D $θ$⁠. For N₁, $T 1 , N 2 , T 2 , z s$⁠, and $R cpa$ the ratios are 0.286/0.011 = 26, 2.94/1.51 = 1.95, 0.619/0.073 = 8.5, 14.99/11.82 = 1.27, 1.07/1.34 = 0.80, and 69.31/49.70 = 1.39, respectively. With the exception of $z s$ the standard deviation of parameter values for the 9-D parameter space is larger than those for the 6-D space with constraints. Finally, the SL for the Tombarra optimal solutions for the 6-D and 9-D parameter space $θ$ are shown in Fig. 13. It is observed that, qualitatively, the variance of SL is about the same for Figs. 12(a) and 12(b).

A method was introduced where frequency-dependent SLs of merchant ships are inferred from acoustic measurements made in bottom-limited ocean environments. The method is based on the observation that the broadband acoustic features at CPA contain significant information content on both the source properties and the seabed characteristics. The inferred SLs are connected to the statistical inference of the PL, which is parameterized by the properties of an effective source depth, a CPA range, and a Wentworth classification scale based on a multi-layered VGS model with empirical constraints.

The method was tested using noise from merchant ships recorded on a vertical line array in the New England Mudpatch in 2017 that were analysed for the information content of both seabed geophysical and source parameters. The data collection was made on two VLAs in about 75 m of water with an approximate isospeed water column. A feature-based maximum entropy algorithm utilized acoustic data samples processed in the 200–750 Hz band from five merchant ships while they were located at the CPA to the VLAs. Features result from coherent effects from both sea surface and seabed. The physical constraints imposed by the hybrid VGS model allow for meaningful comparisons of measured and modeled PLs using estimated seabed parameters. The inferred implicit SLs for five merchant ships were then compared to the best fit of the Wales–Heitmeyer model. Above 200 Hz, the Wales–Heitmeyer model captures the inferred frequency dependence of the SLs.

A limited error analysis was performed on the use of the empirical parameter constraints enacted in Table I. Instead of using the empirical relationships for $K g ( N ) , τ ( N )$⁠, and n(N) in Table I for the first sediment layer, the parameter space was extended from six dimensions to nine dimensions. For the specific type of sediment where the data samples were collected, a clayey sand silt (mud), the optimal values in the 9-D space suggest that the error of determining SL incurred by using the empirical relationships X(N) appears small. Other sediment types may require additional study to establish suitable constraints for $K g ( N ) , τ ( N )$⁠, and n(N).

This work received funding from Office of Naval Research Contract No. N00014-22-P-2007.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from William Hodgkiss upon reasonable request.