Coherence extrapolation for underwater ambient noise Open Access

In shallow water, the ambient noise field resulting from wind-driven breaking waves and rain interacts with the bottom and hence can be inverted for remote exploration of seabed sediment properties.^1–8 The noise field is normally measured at a vertical line array (VLA), and the resolution of geoacoustic parameters inferred from such data is strongly affected by the array aperture. This paper describes a technique to extend the aperture of a VLA by extrapolating the noise coherence measured at an N-element VLA to approximate the coherence of an effective N_e-element VLA where N_e > N with the same element spacing.

For inference of seabed parameters from surface-generated ambient noise, two lines of processing focusing on different features of the seabed have been developed: The coherent approach⁵ (referred to as the passive fathometer) that uses the amplitude and phase information of surface-generated and seabed-reflected noise wavelets to produce an image of the seabed layering structure and the incoherent approach^3,4 that takes the ratio of the power of upward-to-downward noise wavelets to approximate the seabed power reflection coefficient. Both processing approaches apply beamforming to determine angular distribution of the noise.

Despite the advantage of larger arrays in terms of beamforming resolution, short arrays offer significant practical advantages in field measurements such as easy deployment, robustness to array deformation due to currents, as well as being amenable for integration into autonomous vehicles. For this reason, methods for synthetic extension of the array aperture are attractive. While a heuristic method based on zero-padding of the noise coherence has been suggested⁹ to enhance features of the estimated power reflection coefficient, the goal of this paper is to describe a formal mathematical approach for noise-covariance extrapolation.

It has been shown that under mild conditions regarding the separation of the VLA from the seabed or sea surface, the noise covariance matrix (required for beamforming) for an N-element array exhibits a Toeplitz-like structure³ with its first row given by the vector s_N = [s₁₁ s₁₂ ⋯ s_1N], where s_ab is the coherence of the ambient-noise field between the ath and bth array sensors (the frequency dependency of s_ab is omitted for brevity). The noise coherence can be written as⁹ 

where f is the frequency, Δ is the inter-element spacing, θ is the grazing angle measured at the VLA, $R̃(θb)$ is the seabed power reflection coefficient, θ_b = sin⁻¹(sin θ c_b/c_r) is the grazing angle at the seabed, c_r and c_b indicate the water sound speed at the VLA and the water-sediment interface, respectively, and

In Eq. (2), $R̃s(θs)$ is the power reflection coefficient at the air-water interface, θ_s = sin⁻¹ (sin θ c_s/c_r) is the grazing angle at the sea surface, and c_s is the water sound speed at the surface. With the change of variables k_z = (f/c_r) sinθ (where k_z is the vertical wavenumber), it can be shown⁹ that s_ab is the inverse Fourier transform of a band-limited signal with spectrum split into $G̃(kz)$ and $G̃(−kz)R̃(−kz)$ over positive (0 < k_z < W_f) and negative (−W_f < k_z < 0) spatial frequencies (wavenumbers), where W_f = f/c_r.

Extension of s_N to represent a longer, uniformly spaced N_e-element array (N_e > N) gives the coherence

where the upper index ^N emphasizes that the elements $s1bN$ (where b > N) are not observed (sampled) but rather are obtained from the s_1b (b ≤ N) according to a specific algorithm. In contrast, $sNe$ indicates a coherence measured with an actual N_e-element array. While the zero-padding method⁹ implemented by forcing $s1bN=0$ for N + 1 < b < N_e might be a reasonable approximation for fast-decaying coherences, a formal and more general extrapolation method for ambient noise coherence is described in the next section.

The extrapolation of a band-limited signal s(z) for $|z|>Z$ based on observations within the interval $|z|<Z$ can be achieved by projecting s(z) onto the prolate spheroid functions (PSFs).¹⁰ In this paper, the array top element corresponds to z = 0. The PSFs ψ_p(z) (p = 0, …, ∞) form a band-limited complete orthogonal basis with spectra entirely contained in −W_f < k_z < W_f. The PSFs are obtained as the eigenfunctions of the integral equation

where the eigenvalues μ_p can be thought as energy ratios that indicate the percentage of energy of the wavelet ψ_p(z) localized at $|z|<Z/2$ relative to its total energy. The PSFs are bi-orthogonal¹⁰ such that

In addition, ψ_p(z) are band-limited for all p, and their Fourier transform has the property¹¹ 

where ν_p is a complex constant that can be found by enforcing Eq. (5). Equation (6) implies that by knowing ψ_p(z) within the aperture, its Fourier transform can be constructed. Then taking the inverse Fourier transform gives ψ_p(z) for z outside the aperture.

Equation (5) allows expanding a function in terms of PSFs as

which, in theory, accomplishes perfect extrapolation of s(z) beyond $|z|<Z/2$ as P → ∞. In practice, this extrapolation can be affected by numerical artifacts for small μ_p, so a threshold for the number of terms P in Eq. (7) must be set. In this paper, computation of the PSFs is performed by discretizing Eq. (4); this yields the discrete prolate spheroidal sequences¹¹ for $|z|<Z/2$⁠. The discrete equation can be solved by commercial software packages such as the matlab script “dpss.m.” The property in Eq. (6) is then used with the inverse FFT to obtain the PSFs for $|z|>Z/2$⁠.

This section demonstrates the extrapolation approach in Eq. (7) applied to simulated and experimental noise coherence functions.

Table I shows the parameters of two environments with a seabed consisting of three and five sediment layers (including a bottom half-space). For both environments, the water depth is 130 m. These parameters are input to Eq. (1) to obtain the coherence function s(z) at f = 750 Hz with constant water sound speed c_r = 1500 m/s. Figures 1(a)–1(c) show the extrapolation of s(z) for the three-layer environment for increasing P. The simulated coherence (shown as solid lines) represents what would be “measured” with a 100-element array with inter-element spacing of 5 cm (4.95 m aperture). The real (and even) part of the coherence is shown for z > 0, while the imaginary (and odd) part is shown for z < 0. Notice that the coherence functions are normalized to unit maximum for display purpose. Because the extrapolated results are better for the imaginary part than for the real part of the coherence, zoomed-in views of the extrapolation are provided only for the latter. This convention will be followed throughout this paper. The extrapolation results are shown as dashed lines: Figs. 1(a)–1(c) correspond to P = 0, 9, and 11 in Eq. (7). For P = 0, the real part is coarsely approximated by ψ₀(z), and the representation of the true covariance is inaccurate even within the aperture. For P = 9 (which includes the PSFs with μ_p > 0.5), the extrapolated coherence matches well within the aperture, and it predicts the oscillatory behavior of the covariance for $|z|>Z/2$⁠. Adding two more PSFs (which corresponds to including all PSFs with μ_p > 0.06) to the summation gives even more accurate results shown in Fig. 1(c). In this paper, selection of the number of terms P was done based on minimizing the difference between the extrapolated coherence and the coherence measured at the aperture. Similar results can be obtained by plotting $∑p=1P|ap|2$ in a loglog scale against the L2 norm of the error within the aperture given by $∑t=1N|s1tN(P)−s1t|2$ [notice that $s1tN(P)$ emphasizes the dependence of the extrapolated coherence on P]. This plot results in an L-shaped trace that provides a graphical tool to select the value P_elbow corresponding to the L-shape elbow. The rationale of the approach inspired in the L-curve method is that even though the coefficients $|ap|$ tend to grow (due to normalization by decreasing eigenvalues μ_p), their contribution to minimizing the errors $|s1tN(P)−s1t|$ is not significant for P > P_elbow. In addition, as pointed out previously, it is a good practice to limit P because problems of numerical instability can occur by adding more terms to Eq. (7) due to inaccuracies in determining the projection coefficients a_p.

Figure 1(d) shows extrapolation results for a more challenging signal s(z), obtained from the coherence function at f = 750 Hz of the five-layer sediment profile described in Table I. In this case, the real part of the coherence exhibits more structure in its tail, as opposed to the decaying-envelope behavior of s(z) in Figs. 1(a)–1(c). Even though the extrapolation results are less accurate, the trend of oscillations is estimated correctly. Similar extrapolation results were obtained for simulated data corresponding to more sparse VLAs with Δ = 0.5 m inter-element spacing.

Throughout this paper, the reason for obtaining better extrapolation results for the imaginary part of the coherence can be seen from Eqs. (1) and (2): Consider that the Fourier transform of the coherence, $F(sab(z))$⁠, is real⁹ and can be expressed as a summation of even and odd functions. Because $G̃(kz)=G̃(−kz)$ and $R̃(kz)=R̃(−kz)$⁠, the even function [responsible for the real part of s_ab(z)] is strongly dependent on the seabed reflection coefficient, which introduces complicated features, while the effect of the seabed in the odd function [which gives the imaginary part of s_ab(z)] tends to cancel out due to symmetry, giving a smoother (easier to extrapolate) function.

The coherence extrapolation was applied to measured data from the BOUNDARY2003 experiment,¹² collected by the NATO-STO Centre for Maritime Research and Experimentation (CMRE formerly NATO Undersea Research Centre). The aperture used in the experiment consists of a 32-element drifting array of omni-directional pressure sensors with Δ = 0.18 m inter-element spacing. The data were collected in a region with water depth of approximately 131 m at a sampling rate of 12 kHz, and the frequency-dependent data covariance was computed from snapshots recorded over a 5-min interval. The 32 × 32 element covariance matrix was made Toeplitz by taking the mean along diagonals, and this is taken as a reference (i.e., ground truth) in this section on the understanding that it differs from the true covariance, which could only be measured accurately with a larger (and denser) aperture array.

Figure 2 shows the extrapolation for N = 12 at frequencies f = 703, 797, 891, 1008, 1102, and 1195 Hz. The extrapolated covariance with N_e = 32 is in agreement with the N_e-element reference covariance, and they exhibit similar behavior as the simulated results obtained in Sec. IV A. Similar to Figs. 1(c) and 1(d), the extrapolated imaginary part (z < 0) closely matches the reference, and the extrapolated real part (z > 0) follows the oscillatory pattern of the reference but is less accurate. In addition, the reconstructed coherence within the observation aperture (⁠$|z|<Z/2$⁠) closely matches the reference.

This paper developed a procedure for the extrapolation of spatially band-limited underwater ambient-noise measurements based on projections onto a prolate spheroid function basis. The algorithm offers a formal approach to synthetic array processing of surface-generated ambient noise as an alternative to more heuristic methods such as coherence zero-padding.⁹ The true coherence within the aperture can be accurately represented by P orthonormal prolate spheroid functions as shown here for experimental and simulated data. Outside the aperture, the algorithm successfully estimates the general oscillation pattern of the noise coherence.

A potential application of the extrapolated results that is subject of further research is on improving the resolution of seabed parameters estimated by geoacoustic inversion.^4,8 This is plausible as long as the behavior of the extrapolated coherence (or quantities derived from this coherence such as the seabed reflection coefficient input to inversion algorithms⁸) can be predicted and accounted for in the forward model used in the inversion method. This can be accomplished by applying the PSF basis not only to the measured data (i.e., for extrapolation) but also to the coherence model⁸ used to generate data replicas while searching the parameter space during execution of the inversion algorithm.

The authors gratefully acknowledge the support of the Office of Naval Research post-doctoral fellowship and the Ocean Acoustics Program (ONR-OA Code 3211) as well as the Natural Sciences and Engineering Research Council (NSERC) of Canada. We would also like to thank the NATO-STO Centre for Maritime Research and Experimentation (CMRE formerly NATO Undersea Research Centre), Chris Harrison, and Peter Nielsen, for providing the BOUNDARY2003 data.