Statistical inference of seabed sound-speed structure in the Gulf of Oman Basin Jason D. Sagersa) and David P. Knobles Applied Research Laboratories, The University of Texas at Austin, P.O. Box 8029, Austin, Texas 78713-8029

(Received 7 February 2014; revised 3 April 2014; accepted 8 April 2014) Addressed is the statistical inference of the sound-speed depth profile of a thick soft seabed from broadband sound propagation data recorded in the Gulf of Oman Basin in 1977. The acoustic data are in the form of time series signals recorded on a sparse vertical line array and generated by explosive sources deployed along a 280 km track. The acoustic data offer a unique opportunity to study a deep-water bottom-limited thickly sedimented environment because of the large number of time series measurements, very low seabed attenuation, and auxiliary measurements. A maximum entropy method is employed to obtain a conditional posterior probability distribution (PPD) for the soundspeed ratio and the near-surface sound-speed gradient. The multiple data samples allow for a determination of the average error constraint value required to uniquely specify the PPD for each data sample. Two complicating features of the statistical inference study are addressed: (1) the need to develop an error function that can both utilize the measured multipath arrival structure and mitigate the effects of data errors and (2) the effect of small bathymetric slopes on the structure of the bottom C 2014 Acoustical Society of America. [http://dx.doi.org/10.1121/1.4873515] interacting arrivals. V PACS number(s): 43.30.Pc, 43.30.Cq [ZHM]

I. INTRODUCTION

Over the last 30 years matched-field based geoacoustic inversion methods have evolved from a maximum likelihood emphasis to a statistical inference or probability distribution emphasis.1 A cursory examination of the literature on geoacoustic inversion and statistical inference methods applied to ocean acoustics reveals that the vast majority of studies have been made with data recorded in shallow water or continental shelf environments. The foundations for seabed geoacoustic inversion were, however, established during the 1970s and early 1980s in the context of deep-water experiments designed to infer bottom reflection loss versus grazing angle. For instance, Mitchell et al. made a systematic study of analyzing broadband multipath arrivals in deep-water environments for the purpose of inferring bottom loss,2 and an example of a technique to infer geoacoustic parameters from bottom loss versus grazing angle measurements is provided by Spofford in Ref. 3. For thick soft sediments, commonly found in the deep basins of the Atlantic and Indian Oceans, upward refraction by a positive sound-speed gradient is the physical mechanism by which the majority of bottom interacting energy is returned to the water column.4–7 Hamilton developed geoacoustic models for such sediments that include a nearsurface sound-speed gradient and a curvature parameter such that the gradient decreases with increasing depth into the sediment.8 In bottom-limited waveguides, acoustic waves will interact with the seabed and upward refraction within the seabed can become an especially important propagation mechanism. Broadband acoustic data whose propagation is dominated by sediment penetrating refracting multipaths a)

Author to whom correspondence should be addressed. Electronic mail: [email protected]

J. Acoust. Soc. Am. 135 (6), June 2014

Pages: 3327–3337

over a sufficiently large spatial scale could be used in a statistical inference approach to estimate the probability distribution for parameters that describe the sound speed in the seabed. The focus of this study is a statistical inference of the sediment sound-speed structure using acoustic data, recorded on analog magnetic tapes in 1977, taken in the Gulf of Oman Basin, a bottom-limited, deep-water ocean environment with a thick, soft seabed. The original magnetic tapes, stored at The University of Texas Applied Research Laboratories (ARL:UT) for the past 36 years in sealed containers in climate controlled rooms, were recently recovered by a lengthy tape restoration process followed by an analog to digital conversion with a bandwidth of approximately 500 Hz. The subset of recovered data studied in this paper are in the form of calibrated acoustic signals generated by ship-deployed Signals Underwater Sound (SUS) explosive sources along a 280 km track that possesses small bathymetry variability. A maximum entropy9,10 (ME) methodology suited for the analysis of ocean acoustic experimental data11 is used to estimate the conditional posterior probability distribution (PPD) for the sound-speed ratio and the near-surface sound-speed gradient from the SUS data. The PPD quantifies the information content about the seabed that can be extracted from a given ensemble of acoustic measurements. The multiple SUS events allow for an estimation of an average error constraint value that uniquely specifies each PPD. An important aspect of any inference methodology is identifying and mitigating large model and data errors because the errors can affect the interpretation of the PPD. In this paper, mitigating a particular model error arising from the variable bathymetry improves the interpretation of the PPD for the sound-speed structure in the seabed. The remainder of this paper is organized as follows. Section II describes the exercise, the acoustic data, and the

0001-4966/2014/135(6)/3327/11/$30.00

C 2014 Acoustical Society of America V

3327

supporting information. Section III describes the statistical inference methodology. Section IV presents and discusses the results, and Sec. V provides a conclusion. II. EXERCISE DESCRIPTION A. Overview

Bearing Stake was a large-scale exercise employing U.S. Navy assets and occurred in the northwest Indian Ocean from January to May of 1977. One objective of the acoustic data collection was to understand the impact of sediment properties on propagation loss and ambient noise in a deep-water bottom-limited area. Although five experimental sites were studied throughout the region, this paper focuses on SUS data collected on a sparsely populated vertical array at Site 1B, shown in Fig. 1(a), located near the entrance to the Gulf of Oman in about 3350 m of water. SUS were deployed on a track that extended northwest from the array location. B. Acoustic source and receivers

The analysis in this paper utilizes a subset of the acoustic data recorded on the Acoustic Data Capsule (ACODAC) vertical line array (VLA) located at 23.5460 N, 61.0942 E (see Fig. 1). The geometry of the ACODAC consists of three clusters of four hydrophones spaced throughout the water column at depths of 300, 1685, and 3320 m plus a single hydrophone deployed on the seafloor. Channel 1 corresponds to the shallowest hydrophone and the channel number increased monotonically with depth. Four fixed gain settings, differing by 27 dB, were applied to the hydrophones within a cluster with

FIG. 1. (a) Map of the Bearing Stake Site 1B exercise region with bathymetry, SUS track, and ACODAC location superimposed, and (b) bathymetry along the SUS track. 3328

J. Acoust. Soc. Am., Vol. 135, No. 6, June 2014

the intent to provide sufficient dynamic range to capture both ambient noise and high-amplitude impulsive sound. The USNS Kingsport deployed both SUS and a towed continuous wave (CW) projector during the Site 1B portion of the Bearing Stake exercise. The SUS track originated near the ACODAC on 18 February 16:00 and proceeded radially outward on a course of 287 East of North until 19 February 05:30 (cf. Fig. 1). The speed of the Kingsport during this time was approximately constant at 6.45 m/s and the SUS track extended to a range of about 287 km. The range between each SUS event and the ACODAC was estimated from the arrival time in the acoustic record and the corresponding position of the Kingsport listed in the ship’s navigational log. The free fall time of the SUS before detonation and the acoustic propagation time between the SUS location and the ACODAC were carefully accounted for in the range estimation. C. Experimental data

The experimental data consist of calibrated pressure time series for each SUS event received on the ACODAC VLA. An initial screening of the data revealed several issues, which limited the number of usable channels. First, channels 4 and 6 were swapped, meaning that channel 4 was actually the second hydrophone in the 1685 m depth cluster instead of the last hydrophone in the shallow cluster. Second, channels 3, 6, 11, and 13 were discarded because of observed electronic crosstalk. Third, channels 1 and 5 were discarded due to low signal levels. Fourth, comparison of timeaveraged ambient noise levels to periodic electronic internal calibrations revealed that electrical noise dominated channels 8 and 12 throughout the experiment. Fifth, channels 7 and 9 had ambient noise levels that could neither be reconciled with the remaining channels nor with previously reported ambient noise levels.12 ACODAC channels 2, 4, and 10 were deemed usable after the initial data screening. Fortunately, due to the swapping of channel 4, there was one usable hydrophone at each cluster depth. Figure 2 shows the time series data over a range interval of about 170 km for the three cluster depths; a rich arrival structure is evident. In Fig. 2, the SUS events are time-aligned in a relative way, with direct path arrivals aligned at relative time equal to zero, with the first bottombounce arrival aligned relatively after that, and so forth. For the reasons that follow, only a subset of the available SUS events is processed in this paper. First, only the SUS events with a nominal source depth of 91 m were considered in the inference approach described in Sec. III. This choice was made on the basis of the number of SUS deployed for each of the four source depths and the desire to use a single source depth for the specialized error function that will be discussed in Sec. III C. Data from other source depths are reserved for validation purposes and will be presented in Sec. IV C. Second, Mitchel et al.13 noted that the distortion and dynamic range of the ACODAC recorder may have affected some of the waveforms. This is manifest as amplitude compression (due to saturation of the magnetic tape) of the leading arrivals in all the SUS events. However, it was observed that the number of leading arrivals affected by saturation decreased with increasing range, and beyond a J. D. Sagers and D. P. Knobles: Gulf of Oman inference

FIG. 2. Range evolution of the measured time series for the 91 m SUS source depth received at a depth of (a) 300 m, (b) 1685 m, and (c) 3320 m.

range of approximately 50 km the saturation artifacts were deemed manageable using the error function developed in Sec. III C. As a result, events with a range less than 50 km were not included. Third, a maximum range of 175 km was imposed. The first reason for this maximum range constraint is that for the time period corresponding to ranges greater than 175 km, the recovery of the data from the magnetic tape became more intermittent due to a deterioration of the time code channel. The second reason for the maximum range constraint is that the decrease in signal-to-noise ratio (SNR) started to eliminate the arrivals associated with the higher sediment grazing angles. In all, 28 SUS events were selected for analysis; the SUS event numbers and range estimates are listed in Table I. Figure 3 shows the received time series on ACODAC channel 4 for SUS event 175. Figure 3(a) shows multiple bottom-interacting arrivals occurring over 18 s where each successive arrival corresponds to an increasing number of seabed interactions. Earlier arrivals penetrate the sediment at lower grazing angles, have ray paths with shallower turning points in the sediment, and shorter total sediment path lengths than the later arrivals. Figure 3(b) is an expanded view of the first four arrivals of event 175. Within each arrival, both upward- and downward-going rays with the same number of bottom interactions can be present. Each arrival can also consist of rays that undergo specular reflection from the water/sediment interface and rays which penetrate the interface and then return to the water column via refraction. The combination of all of these ray multipaths results in a rather complicated structure within each arrival, whose details are sensitive to the sediment sound-speed profile. The inference methodology presented in Sec. III will describe how seabed soundspeed parameters are inferred from the arrival structure. J. Acoust. Soc. Am., Vol. 135, No. 6, June 2014

TABLE I. SUS event number and range. k

SUS No.

Rng. [km]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28

75 79 83 85 91 95 103 105 107 109 113 117 119 121 125 127 129 133 135 137 143 151 153 165 167 169 173 175

50.7 54.5 58.4 60.4 64.3 68.2 74.0 75.9 77.8 79.7 83.6 87.4 91.3 95.1 99.0 102.8 106.6 110.5 114.3 118.2 125.9 137.4 141.2 156.6 160.5 164.3 168.2 172.2

J. D. Sagers and D. P. Knobles: Gulf of Oman inference

3329

FIG. 3. Received time series on ACODAC channel 4 for SUS event 175; (a) 18 s time window showing 22 distinct arrivals, and (b) 1.8 s time window showing the first four arrivals.

D. Supporting information

We now discuss the available supporting information for the statistical inference analysis of Sec. IV. The prior or supporting information includes the bathymetry along the SUS track, measured sound-speed profiles (SSPs), and previous geoacoustic and geologic information. The bathymetry of the region was obtained from a National Oceanic and Atmospheric Administration (NOAA) bathymetry database.14 The interpolated bathymetry along the SUS track is presented in Fig. 1(b) and shows a nominal water depth around 3350 m with maximum excursions of about 125 m (3.7% of the nominal depth) over the entire track. This bathymetry information is used to quantify the sensitivity of and account for the effects of rangedependence on the joint marginal probability distribution for the sound-speed ratio and the near-surface gradient. Figure 4 shows an SSP measured near the ACODAC VLA.15 An additional SSP was also measured approximately 148 km from the VLA on the SUS track. A preliminary inversion analysis attempted to capture some of the effects of a range-variable SSP on propagation by including in the inversion model space a relative weighting factor of the two profiles. It was found that the inversion for seabed sound-speed parameters was not sensitive to such a weighted average of the two water-column profiles. As a result all the computations shown in Sec. IV employ the SSP measured near the array (and shown in Fig. 4) for an assumed range-independent

SSP. One should not conclude however that range variations in the SSP are not a potential source of model error. Previous physical measurements of the sound-speed depth structure of the seabed at Site 1B are limited, and there is not enough information to support a range-dependent geoacoustic profile. Sedimentation and the geologic structure of the seabed in the Gulf of Oman region have been previously described by Stewart et al.,16 and a seismic survey was reported by White and Klitgord17 at the entrance to the Gulf of Oman. The Site 1B location is contained within this seismic survey. Seismic reflection records suggest that the first sediment layer is around 1250 m thick and is underlaid by two 1000 m thick layers. Hamilton and Bachman refer to Ref. 17 in constructing geoacoustic models for the Bearing Stake region prior to the experimental measurements.18 They report that the Gulf of Oman Basin is filled with horizontally stratified turbidites with an underlying basement of volcanic rock. A shallow sediment core near the ACODAC array suggested that the near-surface sediment consisted of two thin layers. The first layer is composed of silty clay with a sound-speed ratio of about 0.98. The second layer is described as a silt sediment with a ratio of about 1.07. From this information Hamilton and Bachman concluded that the effective sound-speed ratio at the water/sediment interface was approximately 1.0. For a second core measurement the near-surface description was that of a silty clay with a ratio of about 0.98. Hamilton and Bachman also concluded that the depth dependence of the sound-speed profile (from twoway travel time measurements) show that the sound-speed gradient is approximately linear over the first 800 m of sediment and then slowly decreases. For the source/receiver ranges considered in the present study, the majority of ray paths have turning points at depths shallower than 800 m and thus the curvature of the sound-speed profile at the deeper depths is neglected and only a linear gradient is considered. Based on core data, Hamilton reported an in situ corrected value for the sediment density q of 1.58 g/cm3 (Ref. 18). Hamilton and Bachman hypothesized that the nearsurface attenuation a had a value around 0.1 dB/(m kHz). The issue of sediment attenuation for the Site 1B location was also studied by Mitchell et al.13 who estimated a from bottom-loss measurements to vary between 0.02 and 0.04 dB/(m kHz). The impact of q and a on the predicted acoustic field will be discussed further in Sec. III B. In summary, based on the available information and the data selected for analysis, it is assumed that the near-surface sound-speed structure of the seabed is range-independent and is characterized by two parameters: A sound-speed ratio and a linear sound-speed gradient. Further, it is assumed that the modeled SSP is range-independent and is the measured SSP made at the ACODAC array. The measured bathymetry will be used to compute an effective range-independent water depth, source depth, and receiver depth. III. INFERENCE METHODOLOGY

FIG. 4. A water column sound-speed profile measured near the ACODAC VLA. 3330

J. Acoust. Soc. Am., Vol. 135, No. 6, June 2014

The received time series of a SUS event is sensitive to the physical parameters of the waveguide. It is because of this sensitivity that information about the waveguide can, in J. D. Sagers and D. P. Knobles: Gulf of Oman inference

principle, be extracted from the measured waveforms. For example, the waveguide depth H, the sound-speed ratio at the water/sediment interface R, and the sediment sound-speed gradient g determine, in part, the relative timing between arrivals and the waveform oscillations within an arrival. In addition to the physical parameters representing the waveguide, the relative timing between multipaths is sensitive to the source/receiver geometry. The approach taken in this study is to utilize prior information on the source and receiver geometry, the bathymetry, and the sound-speed profile in the water column to mitigate uncertainty in the inference of a PPD for R and g. The analysis is performed in the context of an ME statistical inference approach. The ME approach does not presuppose a form for the distribution of data and model errors, but instead it requires an estimate of an average constraint value to uniquely specify a conditional PPD. The specifics of the ME approach, including a method to set the constraint value for ocean-acoustics problems, were previously described by Knobles et al.11 The basic steps of the ME approach include the selection of a data ensemble, the definition of the model space, the choice of an error function, a methodology to set the ME constraint, and the computation of conditional probability distributions for model space parameters. These steps are briefly discussed in Sec. III A through Sec. III F. A. Ensemble selection

The viewpoint that drives the selection of the ME approach to extract a PPD from the measured data is that, as a result of spatial and temporal inhomogeneities in typical ocean environments, a single acoustic measurement in the ocean should be viewed as a sample from many possible acoustic field microstates. From this view emerges the idea that geoacoustic parameters, such as those that represent the sound-speed structure of the seabed, should be probabilistic in nature. Consistent with this idea, the ME approach utilizes an ensemble of data samples CD consisting of individual measurements or samples Dk drawn from the same population. In this work, a single data sample k consists of 5 to 200 Hz bandpass filtered acoustic time series pmeas k; i ðtÞ on ACODAC channels i ¼ 2, 4, and 10 for a single SUS event, and the ensemble consists of the 28 SUS events listed in Table I. The common feature of the 28 data samples is that the propagation track is the same starting from a range of 50 km to the ACODAC array. The argument that all 28 data samples can be placed in the same population is that the unknown model errors induced by spatial and temporal fluctuations of the waveguide properties encountered by the propagation from all 28 SUS events are of a similar nature. B. Model space

An acoustic propagation model maps a set of environmental and source/receiver geometrical variables M into a predicted acoustic field PðM; f Þ. PðM; f Þ satisfies a linear response problem, i.e., PðM; f Þ ¼ GðM; f ÞSðM; f Þ, where G is the Green’s function for the waveguide, S is the source function, and f is the acoustic frequency. It is understood that G depends on the full set M whereas S, in our case a SUS explosive source, depends only on the depth of the J. Acoust. Soc. Am., Vol. 135, No. 6, June 2014

source. The specific propagation code utilized here to construct GðMÞ is a ray-tracing program named GAMA90.19,20 The inverse Fourier transform of the band-limited complex spectral output of GAMA90 multiplied by an analytic form of the SUS source spectrum21 produces acoustic time series at receiver depths zrec ¼ [300 1685 3320] m that are of a similar form as Dk . The choice of the GAMA90 code was made on the basis of both a need for accuracy and efficiency. In regard to accuracy, GAMA90 includes specular reflection at a layer interface, traces rays into the body of the sediment layers, and also includes multiple traversals in each of the layers. The need for efficiency arises because a numerical implementation of the ME approach required about 5  104 forward calls to GAMA90 to adequately sample GðMÞ in a Monte Carlo integration scheme. For a ray tracing model like GAMA90, the number of ray multipaths does increase with increasing range, but the computational time is still roughly 3 orders of magnitude faster than for a normal mode code because of the large number of modes required for a deep-water, frequency-dependent computation. The model vector M consists of parameters that are assumed known as well as a set of unknown hypothesis parameters H. It is the joint and marginal distributions for the parameters in H that are ultimately sought, but it must always be remembered that the final distributions reflect the a priori knowledge that defines the inference problem. The a priori knowledge is reflected in two ways: The choice of the parameter representation in H and the selection of lower and upper search bounds on the parameter values. In this analysis, the notable fixed parameters include the source depth zsrc , the receiver depths zrec , the water depth H, the range between source and receiver r, the sediment layer thickness Hsed , the sediment density at the water/sediment interface q, and the sediment layer attenuation a. The rationale for fixing these parameters is that they are adequately characterized by supporting information or that they have little effect on the acoustic inference for this environment. The specific zsrc for each SUS event was found to within 0.5 m by optimizing the fit between the Wakely model spectrum21 and the Fourier transform of pmeas k; i ðtÞ. Recovered experimental logs reported zrec . The water depth H is approximately independent of range but does have excursions around a mean value (cf. Fig. 1). While the bathymetry variability is small, it can affect the details of the multipath arrival structure with increasing range. Because GAMA90 is a range-independent model, the effective waveguide method of Harrison and Siderius22 was employed f to obtain effective values of the source depth zef src , the receiver f ef f , and the water depth H . Using the measured badepth zef rec thymetry, the effective values are obtained for each SUS event according to the range values in Table I. This was an important consideration to ensure that small shifts in arrival times of the multipaths were not caused by the misapplication of a rangeindependent forward model to an environment containing a slowly varying bathymetry. Performing the inference with a fixed average water depth can result in model error and the impact of this model error is explored in Sec. IV. The argument for the sound-speed representation of the seabed to be characterized only by R and g was previously made in Sec. II D. The parameter bounds on R were 0.95 to J. D. Sagers and D. P. Knobles: Gulf of Oman inference

3331

1.05 and the bounds on g were 0.5 to 5.0 s1 , consistent with the range of typical soft sediments found in the deep ocean. Additional geoacoustic parameters are fixed. The density of the thick first layer was fixed at the reported value of 1.58 g/cm3 . This is justified by the work of Rutherford and Hawker23 who established that density gradients can be ignored and only discontinuities in the density are important. The first layer is of sufficient thickness that all of the energetic rays for SUS ranges between 50 and 175 km have turning points within the layer, thus making the geoacoustic properties of the basement unimportant for this analysis. Preliminary modeling of the measured data demonstrated that an effective value of a ¼ 0:04 dB/(m kHz), independent of depth and with a linear frequency dependence, approximately gave the observed decay of the bottom interacting arrivals with time. While a does affect the decay of the arrival structure, it has a minimal impact on the timing of the arrival structure. It is important to emphasize that the focus of the present study is the extraction of information on R and g, and that the statistical inference of the probability distribution for aðzÞ and its frequency dependence lies outside the scope of the present paper. C. Error function

An error function EðH; Dk Þ quantifies the difference between the measured and modeled data. We note that the ME approach is valid for any selection of an error function. The PPD is of course sensitive to the choice of the error function, but ME places no constraints on its mathematical form. In this analysis, the error function was heuristically constructed to have sensitivity to two components of the time series: The relative arrival times and the waveform structure within each arrival. Because this is a multi-step error function, a detailed description is given. mod (1) Energy envelopes smeas k; i ðtÞ and sk; i ðtÞ are calculated by applying a sliding standard deviation function24 of period 60 ms. The usable data are truncated when smeas k; i ðtÞ falls below 2.5 times the standard deviation of the noise background. (2) A peak-finding algorithm finds the amplitude of each arrival envelope and normalized envelope functions s mod s^meas k; i ðtÞ and ^ k; i ðtÞ are computed. Equal amplitude weighting of all the arrivals is important so that a future application of the cross correlation function does not penalize the reduced SNR of the later arrivals, whose relative arrival times are very sensitive to R and g. (3) A coarse time alignment of measured and modeled data is accomplished by finding the peak cross correlation ^mod between s^meas k; i ðtÞ and s k; i ðtÞ. This coarse time alignment and pmod was attempted with pmeas k; i , but frequent misk; i alignment of arrival number occurred. The coarse alignment using the normalized envelope function proved extremely robust in this regard. (4) A sliding average of period 20 ms is applied to the pres^mod sure time series to produce p^meas k; i ðtÞ and p k; i ðtÞ. The averaging is necessary to allow the model to fit the major oscillations within an arrival and to deemphasize fitting every small oscillation. A fine time alignment is then accomplished via the maximum cross correlation between

3332

J. Acoust. Soc. Am., Vol. 135, No. 6, June 2014

these signals with a restricted time lag of 50 ms relative to the previous coarse alignment. This sets the final alignment s between the measured and modeled data. (5) The numerical value of the error function is computed according to pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi EðH; Dk Þ ¼ h Ek; i Fk; i ii ; (1) where h i mod ^ ðtÞ ? p ðtÞ ; Ek; i ¼ 1  p^meas k; i k; i

(2)

h i mod ^ Fk; i ¼ 1  s^meas ðtÞ ? s ðtÞ : k; i k; i

(3)

The symbol ? signifies the normalized cross correlation value at s and hii is an average over receiver. Equation (1) has sensitivity to both the timing of the energy of an arrival as well as the oscillation structure within an arrival. Equation (1) was carefully constructed to be insensitive to the amplitude of the arrival structure, which allows us to utilize a fixed value for a and helps to mitigate any saturated arrivals present in the measured data. The error function sensitivity will be discussed further in Sec. IV A and Fig. 6. Figure 5 shows the two components of the error function for a hypothesis vector of R ¼ 1:04 and g ¼ 4 for SUS event 175. This particular hypothesis vector leads to a modeled time series that compares poorly with the measured data. The SNR decrease is seen in Fig. 5(a) by the decrease of arrival amplitude with time and is seen in Fig. 5(b) by the increase of the noise floor in s^meas ðtÞ. The time series is truncated around 17 s because of the SNR criterion of 2.5 imposed on the measured data. It is observed from the figure that the first arrival near 0 s is well aligned in time and that the hypothesis seabed properties cause the later modeled arrivals to appear too early relative to the measured arrivals. For this data sample and H, the numerical values of the error function components are E28; 4 ¼ 0:864, pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi F28; 4 ¼ 0:381, and E28; 4 F28; 4 ¼ 0:574. When all three usable channels are averaged EðH; D28 Þ ¼ 0:551. D. Posterior probability distribution

The goal of a maximum entropy approach is to compute the most conservative posterior probability distribution

FIG. 5. (Color online) Comparison of (a) p^ðtÞ and (b) s^ðtÞ for SUS event 175 and channel 4 for measured and modeled data. J. D. Sagers and D. P. Knobles: Gulf of Oman inference

hEik ¼

h i K 1X ^ k Þ; Dj E HðD K j¼1

(6)

where the constraint equation for the average error is ð hEik ¼ dH PðHjDk Þ E ðH; Dk Þ :

(7)

XH

For each data sample Dk in CD , the value of H which minimizes EðH; Dk Þ is determined via a global optimization routine. The global solution that produces the minimum E is ^ Substituting the canonical distribution Eq. (4) denoted by H. and the partition function Eq. (5) into the constraint equation Eq. (7) allows one to solve for the value of bk given the value of hEik provided by Eq. (6). F. Computation of marginal distribution for R and g

FIG. 6. (Color online) Error function contours for the effective waveguide model and simulated data for four SUS events. Crosses indicate the location of the minimum E and circles indicate the true parameter value.

PðMjDÞ which is the probability that the model hypothesis M is correct given the data D. The ME PPD is P k ðHjDk Þ ¼ PðHÞ

exp ½bk E ðH; Dk Þ ; Zðbk Þ

(4)

where Zðbk Þ ¼

ð

dH PðHÞexp ½bk EðH; Dk Þ :

(5)

The PPD can be reduced to marginal distributions for specific components of H. For example, the marginal distribution for the parameter value h‘ for the kth element of CD is ð P k ðh‘ Þ ¼ dh1 dh2    dh‘1 dh‘þ1 dh‘þ2    dhN  Pðh1 h2    h‘1 h‘þ1 h‘þ2    hN jDk Þ:

(8)

The numerical implementation of the ME approach is one of quadrature for Eqs. (5), (7), and (8). For the specific case in this study, the PPD is the same as the joint probability distribution for R and g. Marginal distributions for R or g require only a single integration of the PPD.

XH

is the temperature, and P k ðHjDk Þ and By convention b1 k Zðbk Þ are known as the canonical distribution and the partition function, respectively. These functions are a cornerstone of modern statistical physics. The canonical distribution is the most conservative conditional probability distribution that satisfies the constraints. It is useful to note that Eq. (4) is Bayes’ rule for a conditional probability distribution and that P k ðHjDk Þ depends on choices of the model space H and the priors, PðHÞ. Prior distributions are simply a statement of what one believes to be true prior to the introduction of the data sample Dk . In fact, one may view a fixed parameter in M as a parameter with a delta function for a prior distribution. In the present study, we have only two parameters in H, specifically R and g. The existing geophysical information about the seabed for the region of interest is that of a soft, claylike sediment. Section III B discussed how the lower and upper bounds for these parameters were assigned based on prior information from the literature. In addition, it is assumed that PðHÞ for these parameters is uniform between the bounds, since there is no preexisting information to suggest otherwise. E. ME constraint

A methodology for estimating values for bk was the major focus of Ref. 11. The same methodology is employed here, namely, the constraint is estimated as J. Acoust. Soc. Am., Vol. 135, No. 6, June 2014

IV. RESULTS AND DISCUSSION A. Results with simulated data

Before discussing the inference of probability distributions with the measured acoustic data, we first study the inference of probability distributions with simulated acoustic data. The purpose of the simulated study is to investigate the impact of one type of model error on the inference process. For computational efficiency, a range-independent acoustic model is applied in the inference analysis. However, small bathymetric variations are known to exist along the SUS track. Because the error function is purposefully constructed to be sensitive to the timing of the arrival structure, the impact of neglecting the range-dependent bathymetry needs to be investigated. Simulated time series data for 28 SUS events were created using R ¼ 0.98, g ¼ 2.0 s1 , q ¼ 1.58 g/cm3 , a ¼ 0.04 dB/(m kHz), and the ranges listed in Table I. Values f ef f ef f were obtained from the nominal values for zef src , zrec , and H of zsrc ¼ 91 m, zrec ¼ [300 1685 3320] m, and H ¼ 3350 m using the Harrison effective waveguide theory22 based on the measured range-dependent bathymetry of Fig. 1(b). Because the simulated data were created using the effective waveguide theory, they would be similar to data produced by a range-dependent acoustic model with variable bathymetry. Random noise was then added to each simulated time series based on a time-averaged ambient noise power spectrum J. D. Sagers and D. P. Knobles: Gulf of Oman inference

3333

generated from the measured data. There exists a PPD for each of the 28 data samples that form the ensemble. Two studies were conducted with the simulated data: The first contained no model error and the second contained model errors arising from using fixed nominal values of zsrc , zrec , and H in the inference instead of the effective wavef ef f ef f guide parameters, zef src , zrec , and H . The purposes of the first study are (1) to evaluate the ability of the error function to locate the true geoacoustic solution in the absence of model error and (2) to understand the intrinsic parameter ambiguity between R and g. The purpose of the second study is to understand the biases introduced by one type of model error during the inference process. Henceforth, we will use the term nominal waveguide model to refer to the case where zsrc , zrec , and H are used to compute the modeled data and f effective waveguide model to refer to the case where zef src , f ef f , and H are used to compute the modeled data. zef rec Contours of E for the effective waveguide model are shown in Fig. 6 for SUS events 75, 119, 137, and 175, which are approximately equally spaced in range between 50 and 175 km. In all cases, including those SUS events not depicted in the figure, E had its minimum value at or very near the true value of R and g. In the absence of model error, the maximum bias in the solution for R and g were 0.0015 and 0.039 s1 , respectively. The conclusion is that the error function has adequate sensitivity to locate the true values of R and g in the absence of model errors. Figure 6 also illustrates the intrinsic parameter ambiguity between R and g and suggests that it evolves with range. This is seen by examining the steepening slope of the 0.15 contour ellipse with increasing range, indicating that at short ranges there is more uncertainty in determining R and less uncertainty in determining g. This is consistent with the physical explanation that rays from short range events will have larger grazing angles and thus penetrate deeper into the sediment before refracting into the water column, with the timing of ray arrivals being sensitive to g. However, the physical trade-off is that rays with larger grazing angles are not as sensitive to R. The figure also shows that the area encompassed by the 0.15 contour ellipse increases with range, suggesting that the uncertainty associated with extracting parameter values from the data could also depend on range. It is expected that this uncertainty would decrease with increasing SNR. The second study with simulated data examines potential biases introduced into the inference that result from using the nominal waveguide model when the simulated data were created using the effective waveguide model. The marginal probability distribution functions (PDFs) for both R and g were computed for each SUS event and are displayed in Fig. 7. The PDFs for R have maximum values that meander about the true value of 0.98 and the longer range events have slightly wider distributions than the shorter range events. Around 115 km the PDFs for g begin to deviate substantially from the true value. The mean values increase toward 3 s1 and a rapid increase in the width of the PDF is observed. Because the previous simulated study found that the error function is capable of locating the true geoacoustic solution in the absence of model errors, the only explanation 3334

J. Acoust. Soc. Am., Vol. 135, No. 6, June 2014

FIG. 7. (Color online) Marginal PPDs for (a) R and (b) g for the nominal waveguide model and simulated data for 28 SUS events. Dashed vertical lines indicate the true parameter value and circles indicate maximum likelihood values.

for the mean value shift is the presence of the model error during the inference process. This is an important observation because if the model error were unknown the inference results might erroneously be interpreted as a reduction in the resolving power of the data at longer ranges or the presence of range-dependent geoacoustic properties. A true interpretation of the results would be the introduction of a bias into the inference process due to unknown model error. However, the conundrum is obvious: The accuracy of the inference interpretation can, in some cases, depend on unknown model errors. This is a simple example of the care that must be taken when assigning physical meaning to statistical inference results. B. Results with measured data

The ME methodology was applied to the experimental data using both the nominal and the effective waveguide models. The marginal PDFs for R and g are displayed in Fig. 8 for both model assumptions. For the nominal waveguide model, the marginal distribution for R in Fig. 8(a) shows that the maximum likelihood values for R vary between 0.95 and 1.015. The nominal waveguide model shows maximum likelihood values for R in the vicinity of 0.98 for ranges less than 80 km, a strong tendency toward values of 1.01 between 80 km and 140 km, with a slightly reduced value at the longer ranges. For the nominal waveguide model, results of the marginal distribution for g shown in Fig. 8(b) indicate that the maximum likelihood values for g fluctuate around a value of 2 s1 for ranges less than 115 km. For longer ranges, J. D. Sagers and D. P. Knobles: Gulf of Oman inference

FIG. 8. (Color online) Marginal PDFs for nominal waveguide model (a)–(b) and the effective waveguide model (c)–(d). Circles indicate maximum likelihood values.

Fig. 8(b) displays the same trend as Fig. 7(b) in that the PDFs for g suddenly undergo both a mean shift and an increase in uncertainty. The similarity to the simulated data example using the nominal waveguide model suggests that employing effective waveguide parameters in the inference methodology is important when trying to infer g from these data, especially at longer ranges. For the effective waveguide model, results of the marginal distribution for R shown in Fig. 8(c) display a general trend of increasing maximum likelihood values of 0.96 to 1.01 for R from 50 to 100 km with a reduction in the value beyond 100 km. Results of the marginal distribution for g shown in Fig. 8(d) display maximum likelihood values which fluctuate around 2 s1 for all ranges and the large shift and increase in uncertainty are not present for ranges beyond 115 km. From the simulated and experimental results presented thus far, it appears that the model error introduced by not accounting for the range-dependent bathymetry can be mitigated through the effective waveguide approach, thus allowing longer-range data to be used in the inference. There are undoubtedly still other unknown model errors, such as the neglect of the possible range-dependence of the SSP and the parameterization of the geoacoustic profile, which could impact the inference results for R and g, but addressing other model errors is difficult in this case because of limited supporting information. A single joint PPD for the effective waveguide model, shown in Fig. 9(a), is constructed by averaging all 28 individual PPDs together and is interpreted as being the ensemble-average statistical inference solution. The ensemble-average maximum likelihood value for R and g is 0.998 and 1.744 s1 , respectively. Contour lines denoting J. Acoust. Soc. Am., Vol. 135, No. 6, June 2014

68%, 95%, and 99% cumulative probability areas are overlaid on the figure to quantify the uncertainty. The uncertainty is comprised of the intrinsic parameter ambiguity (cf. Fig. 6), uncertainty due to data noise, and uncertainty due to model error. The average marginal PPDs for R and g are shown in Figs. 9(b) and 9(c) and have mean values and standard deviations of 0.990 6 0.024 and 2.136 6 0.727 s1 , respectively. The widths of the marginal distributions

FIG. 9. (Color online) (a) Ensemble-average PPD for the effective waveguide model with 68, 95, an 99% probability contours overlaid. The cross marks the ensemble-average maximum likelihood solution. Marginal distributions for (b) R and (c) g with lines marking the mean value plus/minus one standard deviation. J. D. Sagers and D. P. Knobles: Gulf of Oman inference

3335

FIG. 10. (Color online) (a) Comparison between measured and modeled time series data generated from the ensemble-average maximum likelihood solution for four SUS depths and a nominal range of 85 km. (b) Same as (a) except the nominal range is 108 km. (c) Overlaid measured and modeled time series for seconds 7 to 9 of the 91 m source depth and 85 km range. (d) Same as (c) except the nominal range is 108 km.

reflect the uncertainty in the joint PPD but alone do not contain information about the intrinsic parameter ambiguity (i.e., the covariance between R and g). C. Solution verification

The inference for R and g was made with data from SUS events having a 91 m nominal source depth. One verification test that can be performed with the inference solution is to see how well it describes data with other source depths. The ensemble-average maximum likelihood solution for R and g was verified by comparing measured and modeled time series at two ranges (85 and 108 km) and four SUS source depths (18, 91, 244, and 457 m), as shown in Figs. 10(a) and 10(b). Many comparisons were investigated and the eight comparisons presented here are typical. In each comparison, the ensemble-average maximum likelihood solution produces a time series that recreates the relative timing of arrivals as well as many of the oscillations within an arrival. Expanded comparisons for the 91 m source depth are shown in Figs. 10(c) and 10(d), demonstrating the level of model agreement on short time scales. In all cases, the level of agreement is good. One additional validation is presented in Fig. 11 by comparing measured and modeled time series for the 244 m SUS source depth received on channel 4 of the ACODAC VLA. Figure 11(a) is similar to Fig. 2(b) except for the difference in SUS depth and the restriction of the data to the 50 to 178 km range aperture. The agreement for the timing of the arrival structure is very good for all ranges. The conclusion is that the ensemble-average maximum likelihood 3336

J. Acoust. Soc. Am., Vol. 135, No. 6, June 2014

FIG. 11. (Color online) Comparison between (a) measured and (b) modeled time series data for the 244 m SUS source depth received on ACODAC channel 4. J. D. Sagers and D. P. Knobles: Gulf of Oman inference

solution for the near-surface seabed properties when used in a propagation model will qualitatively reproduce many of the features of the impulse response of the Gulf of Oman Basin. V. SUMMARY

Calibrated time series data generated by SUS explosive sources in the Gulf of Oman Basin in 1977 were recently recovered from the original analog magnetic tapes. These acoustic data, along with a sufficient level of supporting measurements afforded an opportunity to explore the amount of information that could be obtained about the sound-speed structure of the seabed. A large number of data samples allowed for a maximum entropy statistical inference approach to be employed to extract probability distributions for the seabed sound-speed ratio at the water/sediment interface and the near-surface sound-speed gradient. The study is characterized by the identification of errors in both the data and the model as well as techniques used to mitigate these errors. Several data and model errors were identified that, left included in the analysis, would have or could have led to incorrect or erroneous interpretation of the range dependence of the marginal probability distributions for the ratio R and gradient g. Data errors were mitigated by excluding from the data ensemble short range data that were largely contaminated by saturation effects on the magnetic tapes and by excluding later arrivals that did not possess adequate SNR. Evidence was presented to show that small bathymetric variability along the source track could cause the marginal probability distributions for g to have a range dependence that could be misinterpreted as range dependence of the seabed parameters. The contribution of this paper can be summarized in three ways. First, the authors have applied a statistical inference methodology to infer distributions of two physical parameters that are important for acoustic propagation in bottom-limited, deep ocean waveguides bounded below by a thick, soft sediment. This represents an advancement over inversion-type results previously reported for similar environments. The work is generally applicable because of the similarity of the Gulf of Oman Basin to other deep-water basins in the ocean. Second, the error function introduced in this paper specifically exploited the relative timing between multipath arrivals to extract information about R and g. For reasons discussed in Sec. III C, utilizing the arrival structure in a robust manner was more complicated than simply crosscorrelating measured and modeled data. A natural consequence of investigating coherent features of the acoustic signal (as opposed to energy-only features) is that model errors might begin to have a larger impact on the interpretation of statistical inference results. A third contribution of this paper was a demonstration of how one specific error, arising through assumptions made about model bathymetry, could impact the marginal probability distributions. The ability to manage model error arising from the bathymetry allowed for the inclusion of longer-range SUS events in the ensemble-average PPD. For the inference of R and g, a larger range aperture for the ensemble average amounts to a wider sampling of ray grazing angles and J. Acoust. Soc. Am., Vol. 135, No. 6, June 2014

turning point depths in the sediment, which translates into a statistical inference solution that is representative of the larger experimental region. ACKNOWLEDGMENTS

This work was supported by ONR awards N00014-131-0109 and a grant from ARL:UT Internal Research and Development. The authors gratefully acknowledge the efforts of Jack Shooter and Thomas DeMary for their role in the tape digitization process. 1

N. R. Chapman and D. P. Knobles, “Perspectives on geoacoustic inversion,” AIP Conf. Proc. 1495, 202–221 (2012). 2 S. K. Mitchell, N. R. Bedford, and G. E. Ellis, “Multipath analysis of explosive source signals in the ocean,” J. Acoust. Soc. Am. 67, 1590–1597 (1980). 3 C. Spofford, “Inference of geo-acoustic parameters from bottom-loss data,” in Bottom-Interacting Ocean Acoustics, NATO Conference Series, edited by W. A. Kuperman and F. B. Jensen (Springer, New York, 1981), Vol. 5, pp. 159–171. 4 O. F. Hastrup, “Digital analysis of acoustic reflectivity in the Tyrrhenian Abyssal Plain,” J. Acoust. Soc. Am. 47, 181–190 (1970). 5 H. E. Morris, “Bottom-reflection-loss model with a velocity gradient,” J. Acoust. Soc. Am. 48, 1198–1202 (1970). 6 J. S. Hanna, “Short-range transmission loss and the evidence for bottomrefracted energy,” J. Acoust. Soc. Am. 53, 1686–1690 (1973). 7 R. E. Christensen, J. A. Frank, and W. H. Geddes, “Low-frequency propagation via shallow refracted paths through deep ocean unconsolidated sediments,” J. Acoust. Soc. Am. 57, 1421–1426 (1975). 8 E. L. Hamilton, “Geoacoustic modeling of the sea floor,” J. Acoust. Soc. Am. 68, 1313–1340 (1980). 9 E. T. Jaynes, “Information theory and statistical mechanics,” Phys. Rev. 106, 620–630 (1957). 10 E. T. Jaynes, “Information theory and statistical mechanics. II,” Phys. Rev. 108, 171–190 (1957). 11 D. P. Knobles, J. D. Sagers, and R. A. Koch, “Maximum entropy approach to statistical inference for an ocean acoustic waveguide,” J. Acoust. Soc. Am. 131, 1087–1101 (2012). 12 R. A. Wagstaff and J. W. Aitkenhead, “Ambient noise measurements in the northwest Indian Ocean,” IEEE J. Oceanic Eng. 30, 295–302 (2005). 13 S. K. Mitchell, K. C. Focke, J. J. Lemmon, and M. M. McSwain, “Analysis of acoustic bottom interaction in Bearing Stake,” Tech. Rep. ARL-TL-79-24 (Applied Research Laboratories, Austin, TX, 1979). 14 The bathymetry data was obtained from the National Oceanic and Atmospheric Administration, available at http://www.noaa.gov (Last accessed 23 May 2013). 15 D. F. Fenner and W. J. Cronin, Jr., “Bearing stake exercise: Sound speed and other environmental variability,” Technical report (Naval Ocean Research and Developmental Activity, NSTL Station, MS, 1978). 16 R. A. Stewart, O. H. Pilkey, and B. W. Nelson, “Sediments of the northern Arabian Sea,” Mar. Geol. 3, 411–427 (1965). 17 R. S. White and K. Klitgord, “Sediment deformation and plate tectonics in the Gulf of Oman,” Earth Planet. Sci. Lett. 32, 199–209 (1976). 18 E. L. Hamilton and R. T. Bachman, “Geoacoustic models of the sea floor: Gulf of Oman, Arabian Sea, and Somali Basin,” Technical Report 357 (Naval Ocean Systems Center, San Diego, CA, 1978). 19 E. K. Westwood and P. J. Vidmar, “Eigenray finding and time series simulation in a layered-bottom ocean,” J. Acoust. Soc. Am. 81, 912–924 (1987). 20 E. K. Westwood and C. T. Tindle, “Shallow water time-series simulation using ray theory,” J. Acoust. Soc. Am. 81, 1752–1761 (1987). 21 J. Wakeley, Jr., “Coherent ray tracing—Measured and predicted shallowwater frequency spectrum,” J. Acoust. Soc. Am. 63, 1820–1823 (1978). 22 C. Harrison and M. Siderius, “Effective parameters for matched field geoacoustic inversion in range-dependent environments,” IEEE J. Oceanic Eng. 28, 432–445 (2003). 23 S. R. Rutherford and K. E. Hawker, “Effects of density gradients on bottom reflection loss for a class of marine sediments,” J. Acoust. Soc. Am. 63, 750–757 (1978). 24 The specific algorithm employed in this work was written by John D’Errico, available at http://www.mathworks.com/matlabcentral/fileexchange/9428-moving-window-standard-deviation (Last accessed 1 April 2014). J. D. Sagers and D. P. Knobles: Gulf of Oman inference

3337

Copyright of Journal of the Acoustical Society of America is the property of American Institute of Physics and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.