Proceedings of the Institute of Acoustics

 

Quantifying acoustic field horizontal variability and acoustic system performance in a canyon with internal tides

 

T.F. Duda, Woods Hole Oceanographic Institution, Woods Hole, MA, USA

B.D. Cornuelle, Scripps Institution of Oceanography, UC San Diego, La Jolla, CA, USA

Y.-T. Lin, Woods Hole Oceanographic Institution, Woods Hole, MA, USA

 

 

1 INTRODUCTION

 

The acoustic signals received by underwater systems have been recognized as complex since the inception of underwater acoustics. The modeling of underwater sound began in earnest in the mid twentieth century, with the purpose of predicting sound levels from sources, including shadow zones, seabed reverberation, and so on. Parallel with this was research into array processing¹ and signal processing² methodology to improve extraction of signal parameters such as direction of arrival (DOA) in the presence of white noise or unwanted interfering signals. Often the signal processing methods make simplifying assumptions about the properties of the signal of interest, for example as a cylindrical wave or plane wave (i.e. the signal is modeled), which can create a framework for parameter estimation but may introduce uncertainty or bias in any given circumstance because the actual signal is in fact more complex than the signal representation in the model.³

 

A theoretically elegant method to utilize the complexity of the signal is matched-field processing,⁴ where the signal model is extended to allow arbitrary complexity, and models that predict fields that match the observed fields are judged to be useful model for the sonic situation, which includes parameters describing the propagation environment as well as signal source parameters such as position. Application of this method has been restricted by the large parameter space that must be considered and the difficulty of making predictions for all parameter combinations. When fully successful, matched-field processing can give precise location information for sources in complex propagation environments.

 

Retreating from matched-field processing and adopting a simplified signal model, mitigating the effects of noise on parameter estimation has received much attention. Using conventional¹ or more recently refined⁵ beamforming approaches to obtain DOA for plane waves and evaluate errors is suitable for plane waves, but this is often too simple of a model for ocean acoustic propagation. In other words, the oceanic Green’s function G(x_s ,x_r) evaluated at many array positions x_r for source position x_s is inconsistent with a plane wave. Here, we present a simple formalism to evaluate source direction (DOA) error estimation for horizontally distributed vertical line arrays (i.e. array beamforming performance estimation) in a complex oceanic environment. After that, we show example DOA results and error estimates made using a time series of computationally simulated sound fields in a canyon, which differ in time through the influence of internal tides. The behavior of DOA estimates is examined and compared with expectations based on the complexity of the sound field. It is found that the DOA can be estimated to within a few degrees, in the simulation scenario, at high signal to noise ratio.

 

 

2 SOURCE DIRECTION ERROR ESTIMATION

 

2.1 Simple least squares estimation

 

The process of obtaining DOA estimates from synchronized acoustic data collected with a horizontal array with quantity M elements is usually an underdetermined inverse problem. The reason is that the range of possible source directions (angles) and the desired angular resolution give a space of size N > M for the model parameter vector, whose elements are plane wave energy at each angle. This can be written

                                                                           

 

where y is the Mx1 vector of complex demodulated data, A is the MxN steering matrix whose columns are the steering vectors at each of N angles, x is the Nx1 model parameter vector representing plane wave energy at each angle, and n is a noise vector the size of y. Recent efforts to efficiently solve this underdetermined system for cases where only a few plane wave signals are expected (i.e. the problem is of a sparse nature) in independent identically distributed complex Gaussian noise have exploited relationships between techniques that minimize the L₂ , L₁ and L₀ norms⁵. Those efforts are called compressive sensing. Here, we take a different approach, assuming that sources are smooth functions of angle, with no prior expectations of a small number of point sources.

 

As done by the compressive sensing authors and by virtually all predecessors, we model the signals as plane waves, so that A is conventional, and a constrained least squares approach is taken. The noise n is taken to be the sum of system noise, Gaussian noise from highly scattered sources that are not locatable, and the leftover sound from nearby sources that are potentially trackable after subtraction of the (hypothesized) major plane-wave constituent. In this case, the objective function is

                                                                       

 

where α2 is a trade-off parameter and μ is a Lagrange multiplier. Minimization of J gives a minimum-variance solution with minimum data/model misfit. The estimated solution is⁶

                                                                       

 

and the data-driven component of the uncertainty of the solution is

                                                                        

 

where I is the identity matrix and R_nn is the covariance of the collected noise processes. With α equal to zero this is a Moore-Penrose inverse solution. The uncertainty depends on the covariance of the noise, as is intuitive. For the case of uncorrelated noise, R_nn = σ_n² I, and a clear dependence on noise variance σ_n² is evident. A filled-out noise covariance matrix will alter the uncertainty for that case.

 

2.2 Gauss-Markov estimation

 

Gauss-Markov inversion formally allows the covariance properties of the noise to enter into the estimate. What is minimized is the dispersion of the inverse result about its true value. An important contributor to this is the covariance behavior of the field parameters, which for the DOA process for unknown targets is no correlation. The Gauss-Markov method explicitly recognizes that noise can exist and minimizes noise using the norm  

The relationship between parameters (vector x) and data (vector y) is again given by (1). The estimated solution is⁶

                                                                         

 

where P_xx is the source vector uncertainty covariance matrix (diagonal). The uncertainty of the solution including both data error and poor resolution is

                                                                     

 

With no prior information one may reasonably set  (no angle preference or extended angle sources, i.e. target correlation) giving

 

   

and 

   

 

2.3 The noise-free best-case scenarios

 

We take as a best-case scenario for a complex environment that there is no system noise and no scattered noise from distant sources that would mimic independent identically distributed complex Gaussian noise. This leaves the complexities in the signals themselves as the sole component of R_nn. This situation is exactly what is routinely modeled in ocean acoustics with propagation models. Thus, by computing the covariance R_ss with respect to the mean reception of signals arriving at synthetic arrays in computational domains, disregarding other contributors to R_nn by setting R_nn = R_ss into the C ~ or C^# estimates for each the estimation methods [(4) and (8)], we find the best-case DOA uncertainty for the complex environment. The matrix R_ss is a function of source position and array positions, and can be written in terms of G(x_s ,x_r) at all x_r in the array. In sections to follow we show how one would compute this for one source position and many array positions in a canyon.

 

The best case scenario can be computed directly from noise-free simulated 3D acoustic fields from single sources, which is the typical output of propagation models. Time series of simulated arrivals at arrays,  with n=1:N for an array of N elements and m=1:M indicating M time snapshots, can be analyzed to form the NxN covariance matrix R_ss. To study how independent noise on each sensor affects performance with respect to signal complexity described by R_ss, a scaled identity matrix can be added to R_ss before it is used in the expressions for the uncertainty C. In this methodology, the Gauss-Markov estimator gives results that depend on the spatially-dependent “noise” process (multipath) whose statistics are quantified by R_ss.

 

 

Figure 1: Output of the tidally forced dynamical model is shown. A snapshot in time of the displacement in of the isopycnal surface that on average resides at 44 dbar is shown in color. The contours show water depth in 100-m intervals. Coronado Canyon is indicated with the arrow. The ending of the word “Coronado” lies on the coastline at approximately the USA/Mexico border, with San Diego Bay shown to the north of this. The red color in Coronado Canyon and northwest of it is a humplike upward wave feature, which tapers to zero in the canyon at this time instant (10 days 16 hours).

 

 

3 SOUND FIELDS IN A CANYON WITH INTERNAL TIDES

 

We explore DOA estimation from simulated sound fields created for 300 Hz point sources at the head of Coronado Canyon. This canyon is southwest of San Diego, CA and lies in the domain of a computational study of internal tide dynamics.⁷ The dynamical study uses the MITgcm dynamical ocean model 8 to compute the internal tides within the model domain that result from boundary forcing with six constituents of barotropic tidal current and elevation. The horizontal grid spacing in the model is 1 km by 1 km, so that the bathymetry of the canyon and the detailed internal tidal response (the baroclinic waves that arise from the action of currents deflecting stratified water in areas of sloping seafloor) are not full resolved. Nonetheless, we use the time-dependent stratification from these computations to describe the water-column conditions in the canyon. The modeling of sound in the canyon uses bathymetry of much higher resolution, and the low-resolution time-dependent sound speed from the model, c(x,y,z,t), is fitted into the high-resolution canyon, which is formally inconsistent from a dynamical standpoint but provides us with representative c(x,y,z,t) at reasonable computational expense with which to examine acoustic behavior. Figure 1 shows an example of the tidally forced dynamical model output.

 

The acoustic modeling is done with our Cartesian three-dimensional parabolic equation (PE) model.⁹ Figure 2 shows time-snapshot model output signal level at two depths, with a point source of 0 dB at the (x,y) origin and 40 m depth in 200-m deep water. The model was set up with the source at the head of the Canyon and the x-coordinate increasing westward down the canyon. The model step in the x direction is one acoustic wavelength λ (5 m), the vertical grid interval is 0.1 λ, and the y grid interval is 0.2λ. The bathymetric grid is taken from the NOAA 3 arc-second Coastal Relief Model¹¹. The PE model was set up to ingest water-column sound speed data at 25-m intervals, interpolated from the 1 km by 1 km dynamical model grid. The bathymetry, also interpolated from the relief model, is input at each acoustic step. The dynamical model output was fitted to the acoustic domain by either truncating at the high-resolution seabed or extending to the seabed. The seabed geoacoustic properties used in the model are consistent with sand (c = 1650 m/s, attenuation 0.5 dB/wavelength, density 1.5 kg/m³).

 

                                                                                             

 

Figure 2: Output of the PE model for one snapshot in time is shown. Sound level from the 300 Hz point source at (x,y)=(0,0) is shown with cylindrical spreading removed (10 log₁₀ (r) is added to the dB level at each distance r from the source). The source at the origin at the left is near the head of the canyon, its eastern end (see Figure 1). Bathymetry in the canyon is contoured at 100-m intervals, with the 500 m contour shown with a wider line.

 

The complexity of the time-varying signals at synthetic arrays (99 m long, 34 elements, 3 m spacing) is depicted in Figure 3. Twenty-four hourly simulations are used. The data that are shown are for two array positions at x = 7700 m, 60 meters below the surface. Tidal periodicity is evident in the patterns of amplitude and phase variation. Figure 3 also shows R_ss for these data. The phase decorrelation rate as a function of inter-hydrophone distance, critical for DOA success, varies along each array. Mean amplitude and amplitude decorrelation scale also vary along each array. If one is to consider the departures of signal from the signal model of plane waves as noise, then the noise is a first-order effect, with the variance of the temporal variations, and of the variations from a fitted plane wave, exceeding the energy of the plane wave. The data are consistent with strong three-dimensional multipath in the canyon area. Before moving on to looking at DOA results and DOA estimation behavior, it is worthwhile noting that the information lies in the simulation results, and there is no substitute for looking at the fields to understand signal structure and multipath in terms of propagation physics and seabed interaction physics.

 

 

Figure 3: The left column shows time series of arrivals (24 at 1-hr interval) at a selected array at x=7700m, y=1165 m, magnitude above and phase below. 12-hour tidal periodicity can be seen. The signals have been phase adjusted (steered) toward the source, so that phase would be uniform for a plane wave coming from the source. The second column from the left shows R_ss for these data, magnitude above and phase below. The third column from the left shows signals, in the format of the first column, for another synthetic array with different center y, − 1435 m. The right-hand column shows R_ss for the y = − 1435 m data. The right-hand data set shows less phase trend in R_ss than the left-hand one, equivalent to a larger correlation scale, also called the coherence length scale.

 

 

4 DOA ESTIMATION AND ERROR

 

Simulated propagation data from a line of non-overlapping synthetic 99-m long 34-element horizontal arrays, also at 60 m depth and x = 7700 m were analyzed using three DOA estimators, the Moore-Penrose 

 

 

the “best-case” noiseless Gauss-Markov

  

and the conventional beamformer, which is the limit as α goes to infinity of

 

 

given by

 

 

All of the methods show results consistent with multiple arrivals. Because the steering matrix A for the line array geometry is consistent with look-direction cone at each steering angle, energy appears at relatively high angles, up to 20 degrees. The three estimators perform comparably. Note that there is bias of about −2 degrees in the average estimated directions. The A matrix that is used is consistent with 180 angles between -89 and 90 degrees.

 

The effect of R_ss on the Gauss-Markov output x^#B can be increased by multiplying R_ss by a factor greater than zero, which represents a reduction in signal to noise ratio for a plane wave arrival versus the true multipath arrival. A factor of 1000 does little to the x^#B. A factor 10 x 10¹⁰ removes most of the arrivals at y > 0, but the others are essentially unchanged.

 

                                                                                             

 

Figure 4 :  (a) Moore-Penrose DOA results x^% for 47 arrays at x=7700 m. (b) For the same data, DOA results x^#B for the Gauss-Markov estimator with array-dependent R_ss inserted as R_nn. (c) Conventional beamformer x*. (d) Time series of x* for one array, showing tidal periodicity.

 

The R_ss matrices for the data along the x = 7700 m, 60-m depth line are almost always singular and difficult to invert, making C^# difficult to compute, so the pseudoinverse is used to compute C^#. Values o C# for the collection of arrays are shown in Figure 5, C ~ are similar. The values are usually below the peak values for the arrival energy (Figure 4a), but do indicate significant correlated errors with angular spread that depends strongly on y.

 

Because Figures 4 and 5 show the broad angular nature of DOA output and the angular distribution of “errors” (the DOA beamformer outputs that differ from a plane wave at the source bearing), both fully predictable based on the R_ss characteristics (Figure 3), one may ask how a simple interpretation of the output would serve a user. A simple thing to do is to take the maximum value of the DOA output. Figure 6 shows this at the top for the conventional beamformer, for 24 hourly acoustic simulations. Clear 12-hour tidal modulation of the peak arrival direction is seen. The temporal mean value of the peaks is biased nearly everywhere, The variation indicated by the +/- one standard deviation markers is low for one region (where the canyon is between the source and receiver) and high for another ( − 500 m < y < 300 m, north of the northern escarpment of the canyon).

 

 

5 SUMMARY

 

Three-dimensional simulations of time-varying 300-Hz sound propagating over complex bathymetry with time-dependent water column properties show much complexity. Signals are consistent with multipath arrivals reflecting from the seabed that map into arrivals at angles up to about 25 degrees when three types of direction-of-arrival estimators (beamformers) are used to analyze time-snapshots along synthetic horizontal arrays. We believe that it is reasonable to use this approach for modeling variability in deep canyon or slope areas because high-frequency nonlinear internal waves are rare or absent in many of these areas, although there are notable exceptions.

 

                                                                                               

 

Figure 5: The Gauss-Markov mean-squared error matrix C^# is plotted for four positions. The best- case scenario is examined, where R_ss is used for R_nn in the formula, that is, the only noise that is allowed is from the multipath (Section 2.3). The center y position in meters is given above each plot. There are many areas of high uncertainty. The axes labels are in degrees relative to the source bearing.

 

                                                                                                                               

 

Figure 6: The peak of the x* beamformer output is shown at each time for each array, in color at top. Below is shown the mean +/ − standard dev. for each array. There is a trend at y < 0 (north) for arrival to be deflected (northerly). Tidal modulation can be seen in the y/time display at the top.

 

The covariance properties of the time-varying arrivals quantify the temporally-variable multipath and waveform distortions. The time-dependence of the data is used to compute data covariance matrices that quantify structures of amplitude and phase along the array that are not consistent with plane waves from the bearing of the source. Use of the covariance matrices in the Gauss-Markov estimator and the error calculations for the Gauss-Markov and simple least squares estimator indicate at this initial stage that these beamformers produce the features in the conventional beamformer, which produces no estimate of error and does not consider noise. Future work will examine other geometric scenarios in this and other sloped regions. Future work with the least- squares estimators to examine the relative effects of multipath, white noise, and unknown colored noise on the behavior of DOA estimates may allow better understanding of acoustic systems in regions like this. The results from the estimators shown here diverge as noise is added to the simulated fields (signal/noise ratio is reduced) prior to DOA estimation.

 

 

6 ACKNOWLEDGEMENTS

 

This work was partially supported by Office of Naval Research GrantsN00014-14-1-0223/N00014-16-1-2372 and N00014-11-1-0701.

 

 

7 REFERENCES

1. D.H. Johnson and D.E. Dudgeon, Array Signal Processing: Concepts and Techniques, PTR Prentice-Hall Inc. (1993).

2. H.L. Van Trees, K.L. Bell and Z. Tian, Detection Estimation and Modulation Theory, Part I: Detection, Estimation, and Filtering Theory, 2^nd Edition, Wiley (2013).

3. E.J. Sullivan, Model-Based Processing for Underwater Acoustic Arrays, Springer/ASA Press, 300–314 (2015).

4. A.B. Baggeroer, W.A. Kuperman, and P.N. Mikhalevsky, “An overview of matched field methods in ocean acoustics,” IEEE J. Oceanic Eng. 18, 401–424 (1993).

5. P. Gerstoft, A. Xenaki and C.F. Mecklenbraeuker, “Multiple and single snapshot compressive beamforming,” J. Acoust. Soc. Am. 138, 2003–2014 (2015).

6. C. Wunsch, The Ocean Circulation Inverse Problem, Cambridge (1996).

7. A.L. Ponte and B.D. Cornuelle, “Coastal numerical modelling of tides: Sensitivity to domain size and remotely generated internal tide,” Ocean Modelling 62, 17–26 (2013). doi: http://dx.doi.org/10.1016/j.ocemod.2012.11.007

8. J. Marshall, A. Adcroft, C. Hill, L. Perelman, and C. Heisey, “A finite-volume, incompressible Navier–Stokes model for studies of the ocean on parallel computers,” J. Geophys. Res. 102, 5753–5766 (1997).

9. Y.-T. Lin, T.F. Duda and A.E. Newhall, “Three-dimensional sound propagation models using the parabolic-equation approximation and the split-step Fourier method,” J. Comput. Acoust. 21, 1250018 (2013). doi: http://dx.doi.org/10.1142/S0218396X1250018X

10. NOAA National Centers for Environmental Information, U.S. Coastal Relief Model, 13 May 2011. http://www.ngdc.noaa.gov/mgg/coastal/crm.html