Synthetic array processing of ocean ambient noise for higher resolution seabed bottom loss estimation

The seabed bottom loss is an important quantity for predicting transmission loss in the ocean. There are a variety of methods used to estimate bottom loss, but a simple one was developed by Harrison and Simons and uses vertically beamformed measurements of ocean ambient noise.¹ This takes a ratio between averaged noise coming directly from the seabed with that coming from the surface. The ratio reveals the loss due to interaction with the seabed, which by definition, is the bottom loss. One of the advantages of the technique is that it produces bottom loss directly without requiring data inversion schemes. In theory, the bottom loss can be estimated exactly using this method, but this would require perfect beamforming and averaging, which implies an infinitely long hydrophone array (i.e., infinitely narrow beams) as well as infinite averaging time. With finite hydrophone arrays, the bottom loss estimate is somewhat smoothed due to beam widths, which is generally undesirable. The smoothing in the bottom loss estimate can shift the location of the critical angle or, if the seabed is layered, significantly reduce the level of interference fringes. When this estimated bottom loss is used directly in propagation models, this smearing can create errors in transmission loss estimates. It can also create errors if the estimated bottom loss is used in an inversion scheme to estimate geo-acoustic properties of the seabed.^2,3

In this paper, the Toeplitz (or approximately Toeplitz) property of the ambient noise cross-spectral-density matrix (CSDM) is used to reduce the degree of smearing caused by the finite beams. This property simply implies that the noise spatial coherence depends only on the distance between hydrophones and not their absolute position in the water column. Exploiting this property essentially provides higher resolution beamforming by making the array appear larger than the physical dimension. For surface generated ocean noise, the CSDM is theoretically expected to be Toeplitz as discussed by Buckingham⁴ as long as the array is not too near the boundaries.⁵ For practical measurements at frequencies of interest, the array can be expected to be several wavelengths from the boundaries so the CSDMs are expected to be Toeplitz, and the techniques developed here should provide higher resolution bottom loss estimates. Simulations are used to generate CSDMs and demonstrate the processing. Although the CSDMs are only approximately Toeplitz, this is sufficient to improve the bottom loss estimates.

The spatial noise coherence in the ocean (or cross-spectral density) is characterized in the frequency domain (ω) by the ensemble average of the product between the pressure field measured at one point with the complex conjugate of a measurement at a second point, and is denoted by

A formula for the surface generated ocean ambient-noise coherence was developed by Harrison⁵ and is appropriate for uniform noise source distributions such as that from wind and rain. Taking two hydrophones vertically separated by a distance z and using Eq. (8) from Harrison,⁵ the noise coherence reduces to

In the preceding, the same axisymmetric geometry with cylindrical coordinates as used in Harrison⁵ is adopted with θ being the angle measured at the receivers on the z axis (θ = π/2 toward the sea surface along the +z axis and θ = −π/2 toward the seabed on the −z axis). The bottom and surface angles (θ_b, θ_s) are defined between 0 and π/2, and for simplicity here, an iso-speed water column is assumed so that θ_b = θ_s = |θ|. Note that in general these angles are related through Snell's law, so this assumption is not necessary but simplifies the notation. A refracting sound speed is considered in the example in Sec. III. The quantity R(θ_b) is the bottom power reflection coefficient (note, the power reflection coefficient is the magnitude squared bottom reflection coefficient and the bottom loss is BL = −10 log R). R_s(θ_s) is the surface power reflection coefficient and c the water sound speed. The water volume absorption is given by a. To correctly add the appropriate amount of volume loss, the ray path cycle distances are used with s_c being a full ray cycle distance and s_p being a partial cycle from the receiver depth to the surface for one up going ray.

Neglecting the volume absorption terms (these are often very small) and combining the terms into the single function, G(θ) this formula can be written,

The neglected volume absorption terms can be added back in later, but this simplified equation helps illustrate the methods being presented here. Note that the simulation that will be done in Sec. III does not assume iso-speed water column or zero attenuation.

Substituting in Eq. (3) k = (ω/2πc) sin θ = (1/λ) sinθ where λ is the acoustic wavelength gives

where G(θ) and R(θ) become functions of k,$G̃(k)$ and $R̃(k)$ after the substitution. The following forward and inverse Fourier Transform pair can then be used,

to re-write Eq. (4) as

with the rectangular window function Π defined as,

and Π(k) = 0 for all other k.

Taking the Fourier transform of C(z) in Eq. (7),

The first term is windowed with the rectangular function such that it contains only terms in the positive part of the spatial frequency spectrum (i.e., 0 ≤ k ≤ 1/λ). Similarly, the second term is windowed such that it only has terms in the negative part of the spatial spectrum (i.e, −1/λ ≤ k ≤ 0). Therefore, the positive spectrum is written,

and the negative part of the spectrum,

Recognizing that by definition R(θ) = R(⁠$−θ$⁠) and R_s(θ) = R_s(⁠$−θ$⁠) so that G(θ) = G (⁠$−θ$⁠) and $R̃(k)=R̃(−k)$ and $G̃(k)=G̃(−k)$ so,

Taking the ratio,

The ratio gives the power reflection coefficient and is essentially equivalent to the result given by Harrison and Simons,¹ using energy flux and beamforming. Like the original derivation, effects such as surface losses R_s are contained in $G̃(k)$ and therefore cancel out when computing power reflection loss. The preceding is defined for 0 ≤ k ≤ 1/λ, which maps to angles between 0° and 90°. However, the derivation given here in wavenumber space is convenient to provide a description of the resolution of the bottom loss estimate and in the next section is used to illustrate the improvements possible for bottom loss estimation.

In practice, the Fourier transform of the coherence function C(z) is implemented using a discrete Fourier transform (DFT) because the measurements are on an array of M hydrophones that are typically in discrete locations (here, assumed equally spaced at Δz = λ/2 separation). The DFT across the hydrophone array is equivalent to conventional (i.e., delay-and-sum) plane-wave beamforming once the wavenumbers are mapped back into angle space, and therefore the resolution of the DFT operation will be the same as that for an angular plane-wave beamformer. The resolution of the DFT beamformer depends on the resolution in spatial wavenumber k. This also determines the resolution of estimated power reflection coefficient $R̃(k)$ from Eq. (13). This is given by the total aperture size because the separation between spatial wavenumbers is Δk = 2π/L where L is the total length of the aperture (array). For a physical aperture of length L, the coherence function C(z) can be estimated between the minimum separation C(0) and the maximum separation C(L).

However, note that in Eq. (4), the value for C(−z) is given by

The coherence function for negative values of z are just the complex conjugate of the values for positive z. Therefore the DFT of the coherence function C(z) can be taken from C(−L) to C(L), which is twice the length of the physical aperture L. This increases the resolution of the DFT beamformer and the estimate for reflection loss to Δk = π/L. The underlying assumption is the validity of the coherence function given in the preceding text. Implementing the DFT of C(z) is essentially equivalent to spectral estimation so in practice a shading window can be used to minimize spectral leakage. Section III will demonstrate this processing using a full wave treatment.

Typical data processing for the noise coherence starts by transforming measured time series data to the frequency domain followed by averaging to estimate the CSDM. The hydrophone data for each channel at angular frequency ω, are written as a column vector p = [p₀, p₁,…, p_M−1]^T for the M hydrophones (T indicates transpose operation). Each entry is determined through a discrete Fourier transform (DFT) of an ambient noise time series measured on each channel, p_m(ω) = $ℱ${p_m(t)}. The number of points in the DFT processing will be referred to as the snapshot size. A single snapshot of the CSDM $C̃n$ is formed as the outer product of the data vector,

where † indicates conjugate transpose operation. Multiple snapshots (N) can be averaged to better estimate the CSDM C,

The equivalent statement to C(−z) = C*(z) in Eq. (14) is that the CSDM C is Toeplitz. This implies the terms down each of the super- and subdiagonals (as well as the main diagonal) of the CSDM are all the same. For measured data, the coherence function can be estimated from the CSDM C by assuming Toeplitz and averaging the diagonals. This forms the symmetric coherence function and bottom loss can be determined using the DFT approach described in Sec. II A. For surface generated noise, Buckingham noted that the CSD matrices are Toeplitz as long as the frequency is high enough to support around 10 or more modes and the hydrophones are not too close to the boundaries.⁴ Harrison derived an expression for this distance from the boundary as approximately, z ≥ $34(λ/sin θc)$ where θ_c is the critical angle. For a somewhat typical critical angle of 20° and 3 kHz acoustic signal, this distance is about 1 m. If noise other than surface generated (such as from ships) dominates, then the CSDMs will not generally be Toeplitz. For practical application, these methods should be restricted to using surface generated noise rather than noise from other sources.

An alternative to estimating bottom loss from the coherence function using a DFT is with the original formulation from the paper by Harrison and Simons.¹ That is, using beamforming and dividing downward steered by upward steered beams. To beamform, each channel is multiplied by a complex weight to properly delay (phase shift) before summing all channels together. The weight for the mth hydrophone steered at angle θ is written, w_m = e^{−im(ω/c)Δzsinθ}, for plane waves arriving at grazing angle θ between the hydrophones separated by distance Δz. Therefore, a beam steered at angle θ is b(θ) = w^†(θ)p. The beam power is B(θ) = b(θ)b(θ)^* which is

According to the original derivation by Harrison and Simons, the bottom power reflection coefficient is estimated by dividing beams steered toward the seabed by beams steered toward the surface,¹ 

To envision the synthetic array processing using the original formulation from Harrison and Simons consider a CSDM with just three hydrophones (for simplicity). The coherences are denoted c_1,1 between hydrophone 1 and itself, c_1,2 between hydrophones 1 and 2, and so on to form the CSDM. Further, the coherence between hydrophones 2 and 1 is the conjugate of 1 and 2, $c1,2*$⁠. That is, by definition a CSDM matrix is always Hermitian, but it is not necessarily Toeplitz (e.g., when not from surface noise but from a signal). In the following text, the matrix on the left side is the most general form of a CSDM, C,

The matrix on the right is implied if the CSDM is also Toeplitz, which means only the difference in sensor spatial position matters and not the absolute position. This Toeplitz property shows the CSDM on the right above consists of just three complex numbers (and conjugates). This is compared with six complex numbers for the general CSDM (and conjugates).

Next, consider adding three imaginary hydrophones vertically below the original three (total of six hydrophones). Lacking any additional information, the CSDM is constructed from just the three real hydrophones, and most of the entries in the 6 × 6 CSDM would be unknown (left matrix in the following text). However, if it is known that the CSDM is Toeplitz, it implies the CSDM is the matrix on the right in the following text,

where here, zeros are entered for the unknown numbers in the CSDM. This implies that if the CSDM has no special properties, adding additional synthetic hydrophones does nothing. However, if Toeplitz, as expected with surface generated ambient noise, then synthesizing additional hydrophones allows a larger CSDM to be constructed with most entries non-zero. This new, larger CSDM can be beamformed and bottom loss estimated in the same way as a CSDM with only real hydrophones [i.e., using Eqs. (17) and (18)]. This same methodology can be used on arrays of sizes larger than three hydrophones and will be applied in the next section to simulated data on a 32 element array. The same shading window as used for the DFT processing can also applied here to the rows of the CSDM to minimize spectral leakage. In theory, this could be extended to construct even larger synthetic arrays (CSDMs) but the effect may be diminished as the number of unknown entries becomes too large.

A simulated noise coherence function is generated using the ocean noise model OASN.⁶ OASN is part of the oases acoustic propagation package that numerically implements a full wave solution producing a CSDM for surface noise in a horizontally stratified media using a spectral integration technique.⁷ Although the coherence function could be generated using Eq. (3), OASN is used because none of the simplifications made in the preceding derivations are assumed, and this methodology has been used in previous ocean noise and propagation studies.⁸ For these simulations, the acoustic frequency is 3500 Hz and the water depth is 200 m. The sound speed in the water column is 1500 m/s from the surface to 50 m depth and then is linearly downward refracting to 1490 m/s at the seabed. The seabed has a 0.75 m layer over a half-space. The layer has sound speed of 1550 m/s, density of 1.5 g/cm³, and attenuation of 0.2 dB/wavelength. The infinite half-space below has sound speed of 1600 m/s, density of 2.0 g/cm³, and attenuation of 0.15 dB/wavelength. The 0.75 m layer gives rise to an interference in the bottom reflection loss, which is also apparent in the beamformer output.

For the simulations, OASN produces a noise CSDM of size 32 × 32 from a 32 element array with the top hydrophone of the array located at a depth of 180 m with hydrophones spacing of Δz = 0.1875 m such that the total array length is L = 5.8125 m. Figure 1 shows the conventional beamforming on the 32 × 32 CSDM as a dashed line. Also shown as a gray line is the conventional beamforming on the synthetic array using the Toeplitz property to synthesize a size 64 × 64 CSDM. The solid gray line shows more depth in the nulls; this is an effect of the higher resolution. Because OASN is a full wave model, the output CSDM is only approximately Toeplitz (e.g., due to effects such as slight differences in absorption terms along the array). The new synthetic CSDM is formed by first averaging terms along diagonals as described previously to form a 64 element coherence function C(z) (from the original size 32 × 32 CSDM). This coherence function can then be expanded into a 64 × 64 CSDM by placing terms of C(z) along super- and subdiagonals of the CSDM and zeros where terms cannot be filled in. As with the DFT processing, C(z) is shaded prior to forming the CSDM to minimize spectral leakage.

The previously determined beam outputs are used to estimate the bottom loss BL(θ) = −10 log R(θ). Results are compared for BL(θ) using beamformer output from the original data contained in the 32 × 32 CSDM with results using beamformer output on the synthetic CSDM of size 64 × 64. These are both compared with the ground-truth bottom loss (exact solution can be determined in several ways, see for example, Jensen et al.).⁹ Figure 2 shows the results. The black solid line is the true bottom loss, the dashed black line is the conventional computation with the CSDM of size 32 × 32. The gray line is the synthetic array data with size 64 × 64 CSDM.

This paper describes a way to improve bottom loss estimates by exploiting the inherent property that the surface-generated-noise spatial coherence mainly depends on the distance between hydrophones and not their absolute position in the water column. This implies the noise CSDM is Toeplitz. With this property, additional entries in the CSDM can be added and this, effectively, creates synthetic hydrophones on an array. Simulations were used to demonstrate improvements in the resolution of bottom loss estimates that use vertical ambient noise directionality.

Support from the Office of Naval Research Ocean Acoustics Program is gratefully acknowledged. The first author would also like to thank the Centre for Marine Research and Experimentation for supporting this collaborative research during the summer of 2012.