Trans-dimensional geoacoustic inversion of wind-driven ambient noise

The acoustic study of seabed layering structure and composition has relied heavily on active-source techniques, although methods using naturally occurring noise^1,2 and man-made sources of opportunity^3–5 have been suggested. From these passive techniques, inversion of wind-driven ambient noise recorded at a vertical linear array (VLA) requires only simple hardware and deployment procedures, it has minimal environmental impact, and its generating mechanism is ubiquitous in the ocean,^1,6,7 making this technique suitable for exploring large geographic areas. A compelling demonstration of the seabed layering information carried by ambient noise is the passive fathometer,⁶ which has been shown to image seabed layering in terms of the time of arrival of acoustic reflections from sub-bottom layers. Significant research has been conducted to improve the passive fathometer's resolution of sub-bottom layers by adaptive array processing techniques,⁶ whereas analytical work has been carried out to understand the impact of discrete interferers.⁸ Postprocessing of passive fathometer seabed images has been proposed using a particle filter⁹ to extract the depth and the strength of acoustic reflectors, which generates sequential data from which other geoacoustic parameters could be estimated.

Unlike the passive fathometer, the aim of this work is full parameter and uncertainty estimation, including layer thicknesses (i.e., a true depth passive fathometer), sound speeds, densities, and attenuations by applying Bayesian inversion to the bottom-loss (BL) estimated from the ambient noise. The contribution of this letter is twofold: First, it compares the results of geoacoustic inversion using wind-driven ambient noise and active-source data for a realistic simulated environment with smooth variations in sound speed and density. Second, the impact on parameter uncertainties due to limited knowledge about appropriate parameterization of the environment is quantified.

Estimating uncertainty due to model selection is a challenging task, particularly for the highly nonlinear problems commonly found in geoacoustic inversion. The results obtained rely on the degree of understanding about the physical process that generates the observed data d. Knowledge of this process allows the formulation of mathematical models with enough complexity to faithfully capture the structure observed in measured data. Comparison between models I_k from a set of K candidates can be based on estimating the Bayesian evidence^4,10 given by the conditional probability $P(d|Ik)$⁠. Although evidence-based model selection has the potential to reduce parameter bias and improve the estimation of parameter uncertainties¹¹ it does not account for the uncertainty due to the choice of model parameterization, as inversions are carried out at fixed parameterizations I_k for k = 1 ,…, K. To address this, the most general method for data-driven model selection is to estimate the posterior probability density (PPD) from which uncertainties, parameter correlations, and other statistical quantities can be obtained. In this letter, trans-dimensional (trans-D) inversion^11,12 with parallel tempering¹³ is used for geoacoustic parameter estimation from BL data derived from ambient-noise measurements at a VLA in shallow water. The trans-D method is a general framework for data-driven inversion, with previous application to analysis of Earth's subsurface elastic properties,¹⁴ inversion of active-source spherical reflection coefficient data,¹¹ and matched-field geoacoustic inversion.¹³ In the trans-D formulation, data and prior information determine the geoacoustic parameters and uncertainties, and also provide a parsimonious parameterization (i.e., the number of parameters used in an inversion are consistent with data and prior information).

To gain insight into the algorithm's performance, the trans-D inversion is first applied to simulated data from a realistic seabed, and then to measured data from the MAPEX 2000 experiment.¹ The results are compared to previous work⁷ that utilized the Bayesian information criterion¹⁵ (BIC) for fixed-dimensional inversion.

Bayes' theorem gives the PPD as¹¹ 

where m_k (defined in the following) is a vector containing the geoacoustic parameters to be estimated, and P(I_k) and $P(mk|Ik)$ are the prior distributions of the parameterization and corresponding parameters, respectively. The likelihood function $P(d|mk,Ik)$ is defined here based on the assumption of Gaussian-distributed residuals d − d(m_k) as

where d(m_k) is a realization of the forward model and C_d is the covariance matrix of the residuals.

The PPD is sampled by a reversible-jump Markov chain Monte Carlo (rjMCMC) algorithm, which uses a generalized Metropolis–Hasting criterion to accept/reject candidate models m_k. Note that the length of m_k is determined here by the number of sediment layers, related to the corresponding parameterization I_k. In this formulation, the inversion parameter space spans multiple subspaces that can vary in dimension, which is referred to as a trans-dimensional space. Details of the sampling algorithm have been fully described elsewhere.^11,13

To provide a more complete search over the parameter space, parallel tempering¹³ is applied by running several interacting Markov chains at different sampling temperatures. Chains at high temperature explore large regions of the parameter space {m_k, I_k} and avoid becoming trapped in local maxima of the PPD. Chains at unit temperature explore local regions in detail, providing unbiased samples of the PPD. The advantage of parallel tempering is evident in sampling multimodal distributions where several families of models can fit the data, as observed in the simulated example presented in this letter. Parallel tempering also improves the acceptance rate (i.e., the number of accepted models relative to the total number of models generated while sampling the PPD) for trans-D jumps, resulting in a more efficient sampling of the PPD.¹³ 

In this work a partition model with k interfaces is used for the seabed. The vector

contains the sound speeds c_l, densities ρ_l, attenuations α_l, and layer thicknesses h_l for each of the k + 1 layers in the partition (the last layer is a half-space). The signal-to-noise ratio (SNR) parameters account for the strength of the wind-driven ambient-noise data (i.e., the useful signal) versus other unwanted sources of noise and are included as parameters to be estimated from the data at F frequencies. As in previous work,⁷ the marginal PPD of SNRs corresponding to simulated data are centered around the true (known) values used to generate the simulated data (not shown). In the case of experimental data, the SNRs depend on unknown factors,⁷ such as the frequency-dependent sensitivity of the array elements, sensor pre-amplifiers, and accuracy of the recording system. Therefore, the estimated values should not be used to infer the wind speed.

Sediments with features such as depth-dependent parameter gradients are represented by a series of layers. The forward model⁷ d(m_k) is based on a ray representation of the ambient-noise field developed by Harrison,¹⁶ and considers the distortion to the seabed reflection coefficient caused by beamforming when estimating the BL.

In this section, marginal probability profiles for the estimated sound speed, density, attenuation, and layering structure of the seabed are presented. The input to the inversion algorithm is the frequency- and angle-dependent BL, computed from the (simulated or experimental) ambient-noise field as the ratio of upward-to-downward energy fluxes.^1,7 In all cases the water column is 130 m deep with a sound-speed profile with a thermocline beginning at a depth of 45 m, as measured during the MAPEX-2000 experiment.⁷ The VLA consists of N = 32 elements spaced by 0.5 m, with the shallowest element at a depth of 88 m.

Inversion of experimental wind-driven ambient noise data can be affected by the presence of discrete interferers (e.g., ships), which violate the assumption of a large (infinite) layer of surface sources. It was shown⁸ that arrivals from discrete interferers with no more than one seabed reflection can obscure the seabed response given by the passive fathometer, although the impact of steeper arrivals (undergoing multiple bounces off the seabed) is reduced. This justifies omission of discrete interferers in the model used for inversion, as long as data at shallow grazing angles are not included. To minimize the potential harmful impact of distant shipping noise in the case of experimental data, in this letter the BL input to the trans-D algorithm is taken over 20 equally spaced angles from 14° to 90°, as inspection of the beamformer output [see Fig. 5(a) of Ref. 1] suggests that data arriving at shallower angles might be contaminated by distant shipping noise.

To take advantage of parallel computing, the PPD is sampled using seven independent Monte Carlo processes, where each process consists of eight interacting rjMCMC chains at different sampling temperatures ranging from T = 1 to 3.5. The results shown in this section are computed from the seven chains at T = 1, and convergence was judged by the stationarity of these chains (i.e., absence of trends in the likelihood and convergence of the chains to similar distributions). In all inversions, P(I_k) consists of a discrete uniform distribution for k = 1–25, whereas $P(mk|Ik)$ uses priors based on experimental laboratory measurements^7,17 of sound speed and intrinsic attenuation versus density given by

where [$cl−$⁠, $cl+$] and [$αl−$⁠, $αl+$] are the upper and lower bounds for sound speed and attenuation, respectively, for a given density $ρl∈[1250, 2050]$ kg/m³.

To gain insight into the performance of the trans-D approach in a complex environment, a simulated seabed with a total depth of 4 m was constructed using measurements of sound speed and density from cores extracted on the Malta Plateau.¹¹ Properties from these cores were partitioned into 120 layers of varying thickness to provide true sediment profiles with fine structure below the resolving power of BL data. Previous work¹¹ using this simulated environment has been carried out in the context of plane-wave reflection-coefficient inversion for active-source measurements. In this work the wavenumber-integration model OASES¹⁸ was used to compute the simulated ambient-noise field at the array. The wind speed (which determines the strength of the noise field) was taken to be 15 kn. Beamforming was performed to estimate the BL at eight frequencies from 300 to 1400 Hz, shown in Fig. 1(a). For this simulation, white Gaussian noise with standard deviation σ_f  = 0.5 dB was added to the BL, and σ_f at each frequency was treated as an unknown sampled implicitly in the inversion¹¹ (i.e., the matrix C_d is assumed diagonal with an unknown element $σf2$ at each frequency). Figure 1(b) shows the marginal PPD profiles for geoacoustic parameters along with the true sediment profiles. At most depths, the support of the marginal PPDs concentrates around the true parameters, indicating good agreement. As in previous work,¹¹ layer discontinuities represent smooth transitions in sound speed, density, and attenuation, resulting in models of lower complexity than the original simulated environment. One important observation is that more structure was resolved by using controlled-source data than with ambient-noise data: The first case yielded models with 6–11 layers (see Fig. 4 in Ref. 11), whereas the passive data suggests models with 4–7 layers as shown in Fig. 1(c). This decrease in resolution is likely related to the smearing effect of beamforming when preprocessing the ambient-noise field, reducing the data information content.

An interesting feature in Fig. 1(b) is the bimodal structure observed in the sound speed and the density over the top 1.5 m of sediment, where the dominant mode has a first interface at ≈0.65 m depth and ignores the pronounced gradient in density of the true model, whereas the secondary mode seems to be driven by this gradient and inserts an interface at ≈0.25 m depth. This secondary mode has been isolated in Fig. 2(a) by using only the samples for which the top acoustic interface is at depths below 0.3 m. In this case, the inversion seems strongly driven by the surface density gradient, while deviating from the true density for depths >0.3 m. The parameterization for the secondary mode is dominated by models with four to five layers [Fig. 2(b)].

Multimodal behavior can result from nonlinearity of the inverse problem (i.e., be real) or from ineffective sampling of the PPD (an artifact of the algorithm). In this work, the marginal PPDs for each of the seven independent chains used to generate Fig. 1(a) all exhibited the same bimodal structure, suggesting that the secondary mode is not the result of poor mixing when exploring the parameter space. This can also be observed by examining the likelihood chain in Fig. 2(c), where the dots and crosses correspond to samples of the main and the secondary mode, respectively. The right subpanel shows the histogram for each case, indicating large overlap between the likelihood of both modes of the PPD. The likelihood of the secondary mode has lower mean, which is expected since Fig. 2(a) is not as good representation of the true environment as Fig. 1(b). Nevertheless, both modes fit the data well, as evidenced by the 95% highest probability density credibility intervals in Fig. 1(a), which contain most data points in both cases.

Trans-D inversion was also applied to ambient-noise data collected at a moored VLA during the MAPEX 2000 experiment.¹ Previous work using the BIC approach to model selection⁷ yielded a three-layer sediment profile in good agreement with core samples of sound speed and density from the region, shown in Fig. 3(a). These results were also in agreement with inversions using an active-source technique based on a single hydrophone and towed source.¹⁹ Figure 3(b) shows the marginal PPDs obtained from the trans-D approach, which resulted in larger uncertainties indicated by wider support of the marginal PPDs. For this inversion, the covariance matrix C_d was obtained from Fig. 7 in Ref. 7. Underestimation of the parameter uncertainties using the BIC can be explained as a result of applying a point estimate for parameterization selection. With the trans-D approach, the parameterization is treated as a random variable with its own distribution determined by the information content of the data [Fig. 3(c)]. This distribution impacts geoacoustic parameter estimation by increasing the corresponding uncertainties and providing more meaningful estimates. For this particular data set, $P(Ik|d)$ indicates a high probability for models with three layers (as found with the BIC study), but it also assigns nonzero probabilities to models with four and five layers.

The trans-D method for model selection was applied to BL data estimated from the wind-driven ambient-noise field in a shallow-water waveguide. Algorithm performance was examined using passive (ambient-noise) data and compared to results from previous studies using active (controlled-source) data. For uniform comparison, the study was carried out with data at the same frequencies, except for 1600 Hz in Ref. 11, which was reduced to 1400 Hz to be below the design frequency of the array used in this work. The marginal PPDs obtained by the passive method have similar characteristics as those from active data (see Fig. 5 in Ref. 11), with good agreement of the sound speed and density with respect to the true profiles. In particular, the depths of interfaces of high acoustic contrast are consistent for both (passive and active) techniques.

Two main discrepancies between passive and active approaches were found: first, the passive study reveals a smaller number of layers (a decrease in depth resolution of geoacoustic features), likely caused by smearing of the seabed reflection coefficient through beamforming. If higher resolution was required, data could be collected with a larger aperture array to mitigate this smearing. Second, the PPD in the passive study has a bimodal structure, in which the main mode provides an overall fit to the true environment, whereas the secondary mode is driven by the near-surface density gradient. This was not the case in the active study, in which the only mode observed in the PPD includes this density gradient.

Application of the trans-D method to experimental data improved previous results in which model selection was addressed using a point estimate based on the BIC. The trans-D inversion gave marginal PPDs with wider support, which translates in more conservative (and realistic) estimates of geoacoustic parameter uncertainties. With this particular data set the trans-D results are not drastically different from the BIC-based results, as the data strongly support models consisting of three layers. Nevertheless, the trans-D method could play a major role with data that support a wider range of parameterizations.

The authors gratefully acknowledge the support of the Office of Naval Research postdoctoral fellowship and the Ocean Acoustics Program (ONR-OA Code 3211). We also acknowledge the NATO Undersea Research Centre for providing the MAPEX2000bis data and Dr. Charles W. Holland for providing the core measurements.