Environmental information content of ocean ambient noise

Ocean ambient noise is generated in a variety of ways such as from waves breaking on the surface, ship sounds, and biologic activity. This acoustic noise propagates through the ocean and therefore contains information about the environment. For example, seabed reflectivity [or bottom loss (BL)] can be estimated from the vertical noise directionality.¹ The water depth and seabed sub-bottom layers have all been estimated using cross-correlations of ocean noise.^2–5 Ambient noise has also been used for determining the surface wind speed⁶ and rainfall amounts.⁷ Given these methods, a natural questions is: which environmental parameters can be estimated and how well? That is, which environmental parameters can usefully be determined and under what circumstances [e.g., required number of hydrophones and signal-to-noise ratio (SNR)]? In this paper, the information content contained in the ocean ambient noise field will be considered using the Cramér–Rao lower bounds (CRLBs).⁸ The CRLB provides a methodology to calculate the minimum variance that an unbiased estimator will have for a parameter of interest. The CRLB has been applied in underwater acoustics to isotropic noise,⁹ in the context of matched field processing^10,11 and tomography^12,13 using a known sound source but here, the focus will be on directional ambient noise as the data and the information about the seabed contained in those fields.

There are several methods that have been developed for determining seabed properties based on beamforming wind-driven¹ or ship^14,15 noise on a vertical hydrophone array. The seabed BL can be determined by taking the difference between the beamformed power coming from the surface to that coming from the seabed. BL can also be combined with inversion methods to estimate seabed properties such as sound speed and density.¹⁶ To use this beamforming methodology in a practical system to map seabed properties over a large area, measurements of ambient noise can be made from a drifting or slowly moving vertical array. For the moving array, data are averaged over a short period of time and an estimate of the seabed properties is made. The exact amount of time averaging is complicated and depends on factors such as seabed variability, frequency band used, SNR, and measurement geometry. For a given single measurement, the dependency of these factors on the accuracy of the parameter estimate can be difficult to determine. The CRLB can be used to compute the lower bound on the estimate variance for a set of unknown parameters. This can be used in the system design, the measurement strategy (e.g., time averaging), as well as in reporting the expected accuracy of results.

In this study, the framework for using the CRLB for ambient noise based parameter estimation is developed. The study site is from the Noise'09 experiment that used a moored vertical array with 16 hydrophones. The data set allowed for estimates of the seabed properties every 15 s (as might be done using a drifting array) and the variance for each these estimates computed over a total of about 1 h and 40 min. The variance of the 15 s estimates is compared with the CRLB. For an efficient, unbiased estimator, the variance should approach the CRLB. Although the estimator used is simple and robust, it will be shown that it does not attain the CRLB and suggests an alternative method for estimating these parameters (with lower variance) may exist. A Monte Carlo based simulation methodology is used to determine the performance of the beamformer based estimation and is compared to the experimental data.

This paper is organized as follows: In Sec. II, the algorithms for estimating BL and seabed parameters are described. In Sec. III the CRLBs are developed for seabed parameter estimation. This provides a measure of the information content that is contained in the ambient noise measurement. Section IV gives results from the Noise'09 controlled experiment and Sec. V uses Monte Carlo methods to simulate realizations of the Noise'09 ambient noise as a way to analyze the methodology with known parameters. Finally, these results and analyses are discussed in Sec. VI.

In this paper, the ambient noise field is considered to be dominated by breaking wave noise. This is not to say other sources could not be used for parameter estimation. Ship noise has been used for BL and seabed parameter estimation^14,15,17 as have biologic sources.¹⁸ Among the parameters that have been shown possible to estimate from ocean ambient noise:

• Seabed properties such as BL and sediment sound speed, density, and attenuation factor;^1,16

• Water column sound speed;¹⁹ 

• Water depth and seabed sub-bottom layering;^2–5 

• Wind speed;⁶ 

• Rainfall;⁷ and

• Hydrophone array parameters such as element locations or array tilt.²⁰ 

Here, only surface wave noise generated by wind is considered. The advantages of using wind driven noise is that it is ubiquitous and the estimation methodologies are relatively simple. The focus will mostly be on the seabed parameters and determining the information content in the noise field. However, the same approach can be applied to any parameter being estimated.

For the seabed properties, the estimation method is based on the research from Harrison and Simons.¹ The method first introduced by Harrison and Simons was not an inversion but a direct estimate of the seabed BL which is defined in terms of the seabed power reflection coefficient $|R(θ)|2$ as a function of grazing angle θ which varies from $0°–90°$⁠,

To beamform the data at frequency ω [adopting an $exp (−iωt)$ time convention], the complex pressure along a vertical hydrophone array is written as a vector, $p=[p1,…,pM]T$ for each of the M hydrophones (suppressing the frequency dependence and T indicates transpose operation). For conventional beamforming, the weight for the mth hydrophone at a depth of z_m with water column sound speed c, shading window h_m, and steered at angle ϕ,

where $k=ω/c$⁠. The vertical steering angle ϕ can vary from −90°–90° with 0° horizontal, and positive angles steered toward the surface and negative angles toward the seabed. Writing the steering weights as a vector, $w=[w1,…,wM]T$⁠, a beam at angle ϕ is $wHp$⁠, where H represents the Hermitian (conjugate transpose). The conventional beam power for a given direction ϕ is (suppressing the frequency dependence),

The expectation $⟨·⟩$ of the pressure field product $⟨ppH⟩$ is the data covariance matrix, K. For L independent snapshots of pressure field p_l, $K̃$ is the approximate value of K taken from averaging,

Forming the covariance allows for the complex data to be averaged before beamforming. The beam power can be written in terms of the approximated covariance matrix,

The seabed power reflection loss is estimated as the ratio of the positive and negative angles according to,

Sometimes these symmetric angles are referred to as conjugate angles since, for conventional beamforming, the steering weights are simply conjugates of each other. Typically, array shading is used to reduce undesired sidelobe levels in $B(ϕ)$⁠. This simply requires multiplying the steering vector by one of many options for shading vectors prior to the multiplication in Eq. (5). Choosing a shading vector depends on the application and here a Taylor window has been applied to all conventional beamforming (implemented using the “taylorwin.m” function in Matlab^21–23).

For parameter estimation, a model for the ambient noise field is used along with beamforming to produce a modeled power reflection loss (denoted $|Rm(θ)|2$⁠). Here, the wavenumber integration model OASES (Ref. 24) is used to generate the surface noise on the vertical array for a given set of environmental parameters. The OASES model directly produces the covariance matrix K(a) for a set of parameters in the vector a. This is an exact covariance matrix based on the theoretical model for surface noise developed by Kuperman and Ingineto.²⁵ 

To obtain an estimate for a, an exhaustive search is performed considering all possible parameter values (e.g., for sediment sound speed, density, and attenuation). For each possible a, a covariance K(a) is computed in OASES, beamformed, and then the power reflection coefficient $|Rm(θ)|2$ computed. The minimum mean squared error (MMSE) is then computed between the measured (estimated from data) and modeled for a set of discrete angles, θ_n,

where θ_n is a vector from 0° to 90° in 1° increments. The angular sampling increment was deemed sufficient based on the slowly varying nature of the reflection loss curves with respect to angle and that the beamforming is more than 10× oversampled for a 16 hydrophone array.

Adaptive beamforming is often used as an alternative to shading the steering vector to suppress sidelobes.^8,26 The adaptive steering vectors are computed using the data itself (⁠$K̃$⁠) and produces sidelobes that minimizes the beam power from directions other than the steering direction. This can be implemented in different ways, but here the minimum variance distortionless processor is used and adaptive weights $wA$ are computed from the conventional weights according to,

The adaptive weights are then used in place of the conventional weights in Eq. (5) and then to compute the power reflection loss.

Adaptive processing to compute BL from ambient noise has been proposed previously.^27,28 However, some caution is needed when directly using the power reflection coefficient or BL from an adaptive beamformer. The adaptive processor gives an undistorted response in the steering direction and minimizes the power in other directions with no constraint on symmetry between upward and downward steered beams. Since the BL algorithm depends on the symmetry of the beam powers (implicit in taking the ratio of upward and downward beams), distortions can occur due only to differences in the beams which can be misinterpreted as due to the seabed. These types of adaptive beamformer distortions of the BL estimate will be illustrated in Sec. IV and is more fully described in Ref. 29.

Even though the BL may contain artifacts, adaptive processing may still be advantageous for environmental parameter estimation. If the modeled BL is computed using the same adaptive processing then, presumably, the artifacts would be similar enough to allow a comparison with adaptively processed measured data. This implies that even if there are distortions introduced in the adaptive beamformer, the model will distort similarly. The minimum mean squared difference is again used between measured and modeled power reflection loss in Eq. (7) for the adaptively processed modeled and measured data. The adaptively processed power reflection loss curves can result in better estimates of the environmental parameters (in terms of the variance of the estimate) as will be demonstrated in Sec. IV.

The covariance K(a) is computed using OASES for a given set of environmental parameters and measurement geometry. The modeling output produces the covariance due only to surface waves so white, additive noise needs to be added. This would be due, for example, to uncorrelated sensor noise and is added as a scaled identity matrix I,

where the a dependency is suppressed. The scaling depends on the relative levels between the surface noise and the sensor noise and is a function of wind speed that changes the surface wave conditions. In the cases here, the SNR is the relevant quantity and the SNR (in decibels) is defined as

where M is the number of sensors (K has dimensions M × M) and tr indicates taking the trace of the matrix.

The information about the unknown parameters a is contained in the Fisher Information Matrix, $J$⁠, with elements defined by⁸ 

where L is the number of independent snapshots that are used to estimate K. The variance of any unbiased estimator of a is bounded by the inverse of the Fisher Information Matrix. For a given parameter estimate $ãi$ from the set of parameter estimates $ã$⁠, the variance is,

That is, the variance on the parameter estimate is bounded by the diagonal terms in the inverse of the Fisher Information Matrix.

Here, covariance matrices are computed using OASES and the derivatives in Eq. (11) are numerically determined using a central finite difference,

using calculations of the covariance matrix for small perturbations to the parameter a_i of $Δai$⁠.

A simple example problem is useful to not only show the methodology but to also test the numerical computation of K and its derivatives against a known analytic solution. To do this, the Cron–Sherman model can be used that calculates surface noise in a non-attenuating ocean without a seabed (i.e., an infinite half-space).^30,31 The covariance matrix for hydrophones separated by a vertical distance $dij=zi−zj$ using the Cron–Sherman model is

Since the covariance matrix is a function of the water column sound speed through $k=ω/c$⁠, the derivative with respect to sound speed, $∂K/∂c$ can be taken analytically. In this way, the CRLB for the water column sound speed can be computed using analytic expressions for the covariance and its derivatives and this compared to the numerical covariance computation (using OASES) and finite difference approximation. For this example, a vertical array is used with 16 hydrophones separated by 1 m (this configuration is chosen to match the experiments described in Sec. IV). The results are shown in Fig. 1 for various SNRs obtained by varying the level of added uncorrelated noise. The figure shows the CRLB on standard deviation rather than variance of the estimate. As expected, the standard deviation of the sound speed estimate approaches zero for high values of SNR. Further, the numerical computation is in good agreement with the analytic result.

The Noise'09 experiment took place off the coast of southern California (approximately 20 km southwest of Point Loma) between January 30 and February 10, 2009. Four 16-hydrophone vertical line arrays with 1 m element spacing were deployed at separation distances of about 500, 1000, and 2000 m. The arrays were moored on the Coronado Bank in about 150 m water depth with the deepest hydrophone 7 m from the seabed for each array. The arrays were self-recording with sampling frequency of 25 kHz and all analysis is done on a single array at a single frequency of 600 Hz. During the experiment, the wind speed varied from 0 to 14 m/s and therefore the surface wave noise also changed during the experiment. The data used here were chosen during a period with sustained wind speeds $≥10$ m/s on February 6, 2009 from 5:00 to 8:00 UTC (only 1 h 40 min of data were available during this 3 h period). The data processing consisted of dividing the time series data into non-overlapping snapshots of length 0.16384 s and transforming to the frequency domain using a 4096 sample length Fast Fourier Transform. The data snapshots were averaged using Eq. (4) with L = 90 (total time of 14.75 s).

Once the data covariance matrix is computed, these were beamformed and the power reflection loss and BL were estimated. One set of 90 snapshots will be referred to as a realization. Each realization was compared against the modeled data for each of the combinations of parameters. Three unknown parameters were considered, sediment sound speed, density, and attenuation and the search space and increments for each parameter are shown in Table I. The parameters that correspond to the MMSE for each realization is saved and results from all realizations were combined to produce a mean and standard deviation for the estimation. A total of 416 realizations of data were used and the histograms for the lowest MMSE for each estimate are shown in Fig. 2 using conventional beamforming. The histograms resulting from using the adaptive beamformer are shown in Fig. 3.

The CRLB for this experimental configuration was computed using Eq. (12). Table II shows the parameter estimates and the $CRLB$ (in terms of standard deviation rather than variance). The mean and standard deviation for the conventional beamforming process are denoted μ and σ and for the adaptive beamforming process as $μA$ and $σA$⁠.

It is worthwhile to consider the BL curves themselves since this is the basis for the seabed parameter estimates. Recall that for the estimates given in Table II, a short time average of 14.75 s of data was used to make a single estimate of sound speed density and attenuation. For the entire data set, this produces a total of 416 estimates, all at the same location (presented as histograms in Figs. 2 and 3). Since the array was moored, an alternative way to process is to average the entire data set of about 1 h and 40 min to form a single covariance matrix and that can then be used in the beamforming algorithm to produce what should be an excellent representation of BL. This long-average BL calculation is shown as the solid lines in Fig. 4. Both conventional beamforming (top panel) and adaptive beamforming (lower panel) of the 1 h 40 min of data are shown.

A comparison can be made between the long-average BL and the BL curves computed using the estimated seabed parameters given in Table II. To accomplish this, the estimated seabed parameters (i.e., the mean values) are input to the OASES noise model to produce a simulated covariance matrix that is beamformed using both conventional and adaptive beamforming and the up/down ratio taken to produce BL. One might think the mean seabed parameters could be put directly into a seabed reflection loss formula (e.g., as described in Ref. 26) to calculate BL directly for comparison with the data estimates, however, the additional steps of computing the covariance and beamforming provide a better comparison with the data that has been processed in a similar way. The main processing effect is due to the beamformer that produces beams of width that depend on the length of the array. When computing BL, the beamwidths cause a “smearing” of the calculated BL curve. A better comparison can therefore be made by smearing model results in the same way the data were processed. These modeled and beamformed-processed BL curves (using mean of the estimates) are also shown in Fig. 4 as the dashed-dotted lines. Note that the adaptive beamformer produces “wiggles” in the BL curve. While these wiggles can actually appear to be real (caused, for example, from layering in the seabed) they are often artifacts that are due to the adaptive beamformer. Note that the same wiggles appear even in the simulation where it is known that no layering exists. The appearance of similar wiggles in the model is why adaptive processing can still be used to estimate seabed parameters even if the BL curve has these artifacts (i.e., the same wiggles appear on data and model). More details about these artifacts are described in Muzi et al.²⁹ To illustrate the true (i.e., “unsmeared”) result, the actual BL is calculated using a theoretical seabed reflection loss model²⁶ along with Eq. (1) for BL. For the calculation of reflection loss, the same estimated mean parameters given in Table II are used. These curves are referred to as the estimation true BLs and are shown in Fig. 4 as the dashed lines. Note these curves are a bit sharper than the ones that have been “smeared” through beamforming.

One of the ways the CRLB can be used is for system design. There are typically many parameter trade-offs and some configurations may inherently be better than others in terms of the information contained in the measurement. To illustrate, the Noise'09 vertical array configuration can be compared to a hypothetical system that uses a horizontal array. In both cases, there are 16 hydrophones with 1 m spacing. The vertical array has hydrophones at depths from 128 to 143 m while for the horizontal array all hydrophones are at 128 m depth and vary in range from 0 to 15 m. Results showing the $CRLB$ are in Fig. 5. The figure shows there is a reasonable amount of information about the seabed even in the horizontal array data although the variance of the optimal estimator is higher than if the data were taken from a vertical array. Interestingly, the beamforming algorithm used with the vertical array would be useless when applied to the horizontal array data due to the orientation. Conjugate beams are perfectly symmetric on the horizontal array and therefore the calculation taking the ratio would not produce the power reflection coefficient (in fact, would result in 0 dB at all angles). Therefore, while the beamforming algorithm applied to a horizontal array would be pointless, the sensor level data does actually contain seabed information. A different estimator (from the ratio of beams) would be needed to extract this information.

The estimation method in Sec. II, and as applied to measurements in Sec. IV, showed that the variance of the estimate does not reach the CRLB. This is an indication that the estimator is not efficient as described by

which for the Noise'09 data the sound speed has an efficiency of approximately $γ=(7.70/32.12)2=5.8%$⁠. Although the estimation method is robust and simple to implement, it is not very efficient indicating there is information in the covariance matrix that does not translate into the estimate after beamforming and dividing conjugate beams. This raises the questions: Is there a more efficient estimator? Answering this is beyond the scope of this article as one would need to consider the complexity and robustness of a different estimator. That is, even if a more efficient estimator exists, it may not be preferred if it is not robust or not practical to implement. A second question that is considered in this section is: Are these estimation efficiencies inherent to the method, or is there something particular about this data set? For example, if there are biologic noise sources or ship noise in the data it may produce higher estimation variances leading to lower efficiencies even if the methodology itself is 100% efficient. Further, in cases such as this where it is not possible to attain the CRLB, how can the variance of the estimation be calculated? Monte Carlo methods can be used to help answer if the efficiency found here is inherent and it provides a framework for producing estimation variances (and biases as described in Sec. V A) in cases like this where the CRLB is not attained by the estimation method.

The Monte Carlo approach is to first create a synthetic data set similar to measurements. Here, the Noise'09 data set will be simulated and then parameters estimated to determine the variance for comparison against the actual measurements and the CRLB. The exact covariance matrix K is computed using a numerical model such as OASES (same as used to calculate the CRLB) and uncorrelated noise is added to this for whatever SNR is desired $KN=K+σN2I$⁠. Random snapshots of data can be computed using an eigenvalue decomposition procedure described in Ref. 26. Here, the covariance is decomposed into a vector of eigenvalues Λ and a matrix of eigenvectors, V,

To create random snapshots, a vector d containing a realization of complex, normally distributed random numbers are created such that $⟨ddH⟩=I$ where $⟨⟩$ indicates taking the expectation. The simulated covariance matrix for a single data snapshot is then,

This can be averaged over L snapshots as was done with the data to get an estimated covariance matrix,

which can then be processed in exactly the same way as the Noise'09 measured data.

A total of 416 realizations were computed using 90 snapshot averages for $K̃$⁠, followed by an estimate of the seabed parameters. The “true” parameters for the seabed properties are taken from the mean values found from conventional beamforming given in Table II. The added white noise was adjusted to produce an SNR of 15 dB. The histogram with results for conventional beamforming is shown in Fig. 6. The results are very similar to those found with the experimental data but it is worth noting the discrepancies. The sound speed of approximately 1653 m/s is slightly lower than the true value of 1667 m/s. However, in terms of a fraction of the standard deviation, the sound speed estimate is only off by 0.44σ. In comparison, the density estimated is too high and has an error which is a much larger factor of 2.12σ. This bias, particularly in the density estimate, will be explored further in Sec. V A.

The previous simulations showed a slight bias in the estimates. The sound speed and attenuation were slightly lower than the true value and the density somewhat higher. As previously mentioned, the CRLB gives the lower bound on the variance of an unbiased estimator. In addition to the beamforming BL algorithm not attaining the CRLB, it is also not unbiased. However, the bias depends on the SNR as well as many other factors including the specific seabed properties.

To determine the bias, Monte Carlo methods can be applied in a way similar to the previous examples. However, in this case, the simulated “data” is created for a variety of SNR values and for each an estimate of the seabed properties is made. One difference between these calculations and those simulating the Noise'09 scenario is the number of realizations is increased from 416 to 1000 (to better insure convergence). The results are shown in Fig. 7. Some features of the plot are worth noting. At low SNR values (below around 5 dB), many of the estimates hit the upper bound of the search space causing the bias curve to appear to flatten. This is an artifact and therefore very low SNR values are not shown. At low SNR values all three seabed properties are biased toward higher values. However, as SNR increases, the sound speed and attenuation then bias toward lower values. These slightly negative values for the sound speed and attenuation estimates are consistent with the simulations done in Sec. V with 416 realizations and 15 dB of SNR. All three seabed parameters rapidly approach zero bias as SNR increases above 15 dB.

Given the expected bias in the estimates, to first order, the bias can be taken away from the estimates in Table II (at 15 dB SNR). That results in a sound speed estimate for the Noise'09 data of 1678 m/s, a density estimate of 1.88 g/cm³, and attenuation of 0.58 dB/λ. As a check, when the bias values are compensated for on the simulated data, the sound speed is adjusted up from 1656 to 1667 m/s (the exact true value), the density is adjusted down from 2.48 to 2.2 g/cm³ (compared the true value of 2.16 g/cm³) and the attenuation adjusted from 0.57 to 0.53 dB/lambda (compared to the true value of 0.54 dB/lambda).

The ocean ambient noise field contains information about the environment and various estimation methods exist to extract quantities such as the seabed properties. The CRLBs) can be used to quantify the variance an unbiased estimator can reach for a given parameter. In this study, the framework for the CRLB from ambient noise measurements has been developed and the calculations were demonstrated. This methodology can be used for system design and for understanding various trade-offs.

For this study, the beamforming technique used to estimate BL was considered along with a parameter estimation for the seabed properties of sound speed, density, and attenuation. The estimation is based on taking short time averages of data and computing BL based on the beamforming algorithm. A model is then used to compute all possible BL curves over a selected search space. The minimum least squared error between the measurement and each of the modeled BL curves were computed and for each realization of data the minimum error gave the parameter estimates. This was repeated for many realizations to give a mean and standard deviation for the estimates for comparison with the CRLB. This is possible to do if the array is not moving as was the case for the Noise'09 data considered here. It was found that the beamforming BL method was not very efficient in terms of reaching the CRLB. Further, the estimates were shown to have a bias that can be significant at low values of SNR. The bias was computed and removed from the original estimates. For the Noise'09 data, the true values of the seabed properties are not known accurately enough to know if the adjustment helps. However, for the simulations, the seabed properties are known and the bias compensation brought the estimated values closer to the true value, making the final estimates extremely good. It was also shown that adaptive beamforming introduced artifacts into the BL curve. However, in the estimation process these artifacts are also introduced into the model. The end result was that the adaptive beamforming estimation of seabed properties had a lower variance than those estimated using conventional beamforming.

The Noise'09 experimental data was considered here and the array was moored; however, practical application of the algorithm would more likely have the array moving to survey the seabed over a desired region. At a given location, only a single realization of data would be used to produce an estimate of the seabed sound speed, density, and attenuation. However, since the estimator does not reach the CRLB, the variance is unknown and this is nearly as important as the estimate itself. To estimate the variance and bias of the beamforming based estimation method, Monte Carlo methods were developed and through simulations were shown to produce results very similar to the measured data. These methods can be used for this estimation method to determine, for example, the variance and bias a drifting array would have.

The authors would like to gratefully acknowledge support for this research by the Office of Naval Research Ocean Acoustics Program under Grant No. N00014-13-1-0632. We also would like to acknowledge the Marine Physical Laboratory of the Scripps Institution of Oceanography, University of California, San Diego for providing the Noise'09 experimental data.