High-Frequency Underwater Acoustic Propagation in a Port Modeled as a Three-Dimensional Duct Closed at One End Using the Method of Images

Abstract

A computer-efficient model for underwater acoustic propagation in a shallow, three-dimensional rectangular duct closed at one end has been developed using the method of images. The duct simulates a turning basin located in a port, surrounded with concrete walls, and filled with sea water. The channel bottom is composed of silt. The modeled impulse response is compared with the impulse response measured between 15 kHz and 33 kHz. Despite small sensor-position inaccuracies and an approximated duct geometry, the impulse response can be modeled with a relative echo magnitude error of 1.62 dB at worst and a relative echo location error varying between 0% and 4% when averaged across multiple measurements and sensor locations. This is a sufficient level of accuracy for the simulation of an acoustic communication system operating in the same frequency band and in shallow waters, as time fluctuations in echo magnitude commonly reach 10 dB in this type of environment.

1. Introduction

Underwater vehicles and divers routinely operate in ports for security and maintenance operations. Communicating with underwater assets using acoustic modems is a critical feature whenever a tether cannot be used and remains very challenging due to large amounts of noise and fading [1–6]. The reader should note that underwater acoustic communication has become a mature field of research: this reference list is by no means exhaustive.

As in many marine operations, simulation tools play a critical role in optimizing a communication system performance and often in reducing the duration of field tests. Simulating underwater acoustic communications in a port is also a very challenging task, due in part to the complex geometry of the environment. A small displacement between acoustic sensors or a change in the channel geometry can result in a dramatic change in the measured impulse response when using high-frequency sound above 10 kHz. As a result, the channel response varies with time and cannot be modeled exactly. Consequently, stochastic models are used to estimate the performance of an underwater acoustic modem [1, 7]. In this case, the accuracy of the acoustic model is measured in terms of statistical moments rather than absolute accuracy in predicting the impulse response of the channel for a specific configuration. This in turn means that a computer-efficient model of the acoustic channel can provide acceptable results when averaged over a large number of simulations in the presence of small geometrical fluctuations.

To better understand the context of underwater acoustic simulation tools, Figure 1 shows an example of top-level architecture for an underwater acoustic network [7]. The purpose of this tool is to predict the behavior of one or more vehicles, each carrying an acoustic modem and completing a specific mission (helm). All the vehicles evolve in a world, which impacts the acoustic communication quality between any two acoustic modems. Here, each acoustic modem is represented as a protocol stack and a sensor within each vehicle. The passing of information through the acoustic channel is handled by a medium model within the world. In its simplest form, this model could be a vehicle and a boat, both equipped with an acoustic modem.

Figure 1 

Example of top-level architecture for an underwater acoustic network.

Acoustic modems transmit series of acoustic (band-limited) impulses, each containing some binary information, to relay messages between the source and the receiver [1–5]. State-of-the-art acoustic modems transmit hundreds or even thousands of impulses within a message, using either phase or frequency modulation (or a combination of both). The probability that this binary information contains errors is a function of the type of modulation, error coding, signal-to-noise ratio (SNR), and signal-to-multipath ratio (SMR) [1]. The SMR is the energy ratio of the direct echo (traveling directly from the source to the receiver) to the total energy in the reflected (or scattered) echoes measured at the receiver. This SMR is especially critical, as it indicates the amount of fading in the acoustic channel. Modeling the SMR is difficult, as it is a direct function of the acoustic channel response.

A critical issue is the amount of processing required to model the entire network operation, especially in terms of the acoustic channel model. As each acoustic source and receiver moves within the medium, the acoustic channel response changes. A very processor-intensive approach consists in using a powerful acoustic propagation model and recalculates the acoustic channel response given the source and receiver location. For example, Beaujean et al. [8] considered the Parabolic Equation (PE) model approach in a previous paper on a similar problem but realized that this approach was simply too processor-intensive for this application. The impulse response can be precomputed for each combination of source and receiver position, but the sheer number of combinations makes this approach impractical as well.

In contrast, stochastic models are an excellent trade-off between processor requirements and model accuracy, so long as a sufficient number of trials are performed to produce meaningful statistical averages. In the application shown in Figure 1, the authors use the Nakagami model [7, 9]. If 𝐴² represents the acoustic energy of each impulse within a message, 𝑚_Nakagami=𝐸{𝐴²}²/Var{𝐴²} is the ratio of the squared expectation (statistical mean) of 𝐴² to the variance of 𝐴². Here, the expectation is estimated using the average across all the impulses contained in a message. 𝐸{𝐴²} is the average energy of the message. Var{𝐴²} is the energy spread. The parameter 𝑚_Nakagami usually varies from 0.5 to 10.

The main difficulty with stochastic models is to reconcile the statistical parameter(s) with the actual environment in which the acoustic modems operate. In this case, the parameter 𝑚_Nakagami must be a realistic function of the source and receiver position within the medium. However, this parameter does not have to be extremely accurate either: two significant digits are sufficient to provide a reasonable binary error prediction. Therefore, a simple acoustic propagation model of limited accuracy may be sufficient to provide a realistic value for 𝑚_Nakagami.

A logical choice is the method of images applied to a specific 3D environment. Although the method is not novel, it may predict the amplitude and location of every echo with a sufficient level of accuracy to calculate the parameter 𝑚_Nakagami. A second option is to convolve the modeled acoustic response with the transmitted modem message to generate an artificial received signal. Once artificial ambient noise has been added, this artificial signal can be decoded. However, this second option is much more processor intensive.

In this context, the authors conduct a comparative analysis between the acoustic response predicted with the proposed 3D model and the acoustic response measured experimentally. A complete sensitivity analysis of the parameter 𝑚_Nakagami using the acoustic model and the field measurements is beyond the scope of this paper, as it requires a complete description of the actual acoustic message transmitted by the source. Instead, the band-limited impulse response of the acoustic channel is studied, using a pulse transmitted within the entire frequency band of the actual acoustic modem. The comparison is made in terms of the relative error in echo magnitude and time of arrival, across a large number of measurements.

The channel of interest is the south turning basin of Port Everglades, Florida, which is similar in shape to a three-dimensional duct open at one end. Unfortunately, most of the research conducted in underwater acoustic propagation in partial enclosures focuses on three-dimensional wedge geometry. Following the seminal work on horizontal refraction by Weston [10], Deane and Tindle [11] presented a model for the three-dimensional acoustic field in a wedge, leading to the calculation of a loss parameter and to the modeling of horizontal refraction. The results were also demonstrated experimentally in a wedge-shaped ocean [12–14]. To solve the three-dimensional Helmholtz equation, Buckingham [12] used normal mode theory to derive the mode shapes based on specific boundary conditions and found that the mode shapes varied with frequency and range within the mode coefficients. Borejko [15] created a representation of the image field in a perfect wedge using the ray integral method.

These powerful techniques become overly complex and computer intensive in the present case, due to the short acoustic wavelength and the geometry of the basin. If this duct is filled with seawater of uniform and constant properties and is partially enclosed between a still sea surface, a silt bottom, and three rigid vertical walls, the computer-efficient method of images [16, 17] can be used to model the channel response given specific acoustic sensor locations. The three-dimensional method of images is mostly used to model enclosed environments in airborne audio acoustics [18–22]. Allen and Berkley [18] developed a model to study the basics of room acoustics and were interested only in the point-to-point transmission between source and receiver. Viveiros and Gibbs [19] predicted the field performance of acoustic louvers using an image model compared with impulse measurements. Also using the method of images, Iu and Li [20] computed the acoustic channel in narrow street canyons, modeled as two parallel, infinitely long planes perpendicular to a horizontal ground. The geometry of the duct representing the turning basin is similar to that of the canyon modeled by Iu and Li [20]. However, significant differences exist, namely, in terms of the ultra-sound frequency band, sound absorption, and characteristic impedance of the medium and the boundaries. It should be noted that while the image method is not overly complicated for parallelepiped geometries, it can become very complex and processor-intensive for other geometries due to the screening for image sources.

In the following sections, the authors provide the detailed derivation of the method of images applied to the duct representing the basin, followed by a description of a set of experiments and a comparative study between the modeled and the measured channel response.

2. Acoustic Model

Consider the modeled duct shown in Figure 2, where the water mass density 𝜌 and sound velocity 𝑐 are constant. For now, we also assume that 𝑆₁(𝑥_𝑆,𝑦_𝑆,𝑧_𝑆) is a point source producing a complex, harmonic spherical pressure wave 𝑝 measured at the receiver location 𝑅(𝑥_𝑅,𝑦_𝑅,𝑧_𝑅). This complex, harmonic pressure wave is the solution to the Helmholtz equation [23]: Δ𝑝𝑥_𝑆,𝑦_𝑆,𝑧_𝑆,𝑥_𝑅,𝑦_𝑅,𝑧_𝑅+𝑘²𝑝𝑥_𝑆,𝑦_𝑆,𝑧_𝑆,𝑥_𝑅,𝑦_𝑅,𝑧_𝑅=0. (1) Expressed in Cartesian coordinates and in the complex domain, the free-field harmonic solution to (1) is 𝑝𝑥_𝑆,𝑦_𝑆,𝑧_𝑆,𝑥_𝑅,𝑦_𝑅,𝑧_𝑅=𝑃,𝑡₀𝑒^{𝑖(𝑘((𝑥𝑆−𝑥𝑅)2+(𝑦𝑆−𝑦𝑅)2+(𝑧𝑆−𝑧𝑅)2)−𝜔𝑐𝑡)}𝑥_𝑆−𝑥_𝑅²+𝑦_𝑆−𝑦_𝑅²+𝑧_𝑆−𝑧_𝑅², (2) where 𝑘=𝜔_𝑐/𝑐 is the acoustic wave number, 𝜔_𝑐=2𝜋𝑓_𝑐 is the angular frequency of the transmitted signal, and 𝑃₀ (in μPa·m) is the acoustic pressure times unity distance, measured at 1 m from 𝑆₁. If we define 𝑅₁ as the distance between the source 𝑆₁ and the receiver as 𝑅₁𝑥_𝑆,𝑦_𝑆,𝑧_𝑆,𝑥_𝑅,𝑦_𝑅,𝑧_𝑅=𝑥_𝑆−𝑥_𝑅²+𝑦_𝑆−𝑦_𝑅²+𝑧_𝑆−𝑧_𝑅², (3) the complex pressure field produced by the point source becomes 𝑝𝑅₁,𝑡=𝑃𝑅₁𝑒^{−𝑖𝜔𝑐𝑡}=𝑃₀𝑒^{𝑖𝑘𝑅1}𝑅₁𝑒^{−𝑖𝜔𝑐𝑡}. (4) The geometry shown in Figure 2 contains a pressure release boundary Π₁, a silt bottom Π₂, and three rigid walls Π₃,Π₄, and Π₅. The duct along the positive 𝑥-axis is open ended. The approach in developing the model is to break this three-dimensional problem down to a combination of two-dimensional acoustic models. For simplicity, all the calculations take place in Cartesian coordinates.

Figure 2 

Duct diagram.

2.1. Method of Images Applied between Boundary Π₅ and the Open-End

We first assume that Π₅ is an infinite rigid boundary, as shown in Figure 3. 𝑆₂ is located in plane Π₆, which contains 𝑛₅ and 𝑠₁𝑟⃗. Since Π₅ is a rigid boundary, the pressure gradient along the normal 𝑛₅ is null, so that (𝜕𝑃/𝜕𝑥)_𝑥=0=0. Based on the method of images [16], the pressure at the receiver is the sum of the pressure generated by the source 𝑆₁ and the pressure generated by the image 𝑆₂, located at equal and opposite distance from the boundary: 𝑃𝑥_𝑆,𝑦_𝑆,𝑧_𝑆,𝑥_𝑅,𝑦_𝑅,𝑧_𝑅=𝑃𝑅₁+𝑃𝑅₂=𝑃₀𝑅₁𝑒^{𝑖𝑘𝑅1}+𝑃₀𝑅₂𝑒^{𝑖𝑘𝑅2}. (5)𝑅₁ is the distance between 𝑆₁ and the receiver. 𝑅₂ is the distance between 𝑆₂ and the receiver: 𝑅₁=𝑥_𝑆−𝑥_𝑅²+𝑦_𝑆−𝑦_𝑅²+𝑧_𝑆−𝑧_𝑅²𝑅₂=𝑥_𝑆+𝑥_𝑅²+𝑦_𝑆−𝑦_𝑅²+𝑧_𝑆−𝑧_𝑅² (6)𝜓^𝑡_𝑙 and 𝜓^𝑎_𝑙 correspond to the angles of transmission and arrival unique to the source 𝑆_𝑙, where 𝑙=1,2 is the image index. Here 𝑙=1 corresponds to the physical source 𝑆₁: 𝜓^𝑡_𝑙=−arctan𝑦_𝑆−𝑦_𝑅²+𝑧_𝑆−𝑧_𝑅²𝑥_𝑅−𝑥_𝑆𝑙,if𝑙=1,arctan𝑦_𝑆−𝑦_𝑅²+𝑧_𝑆−𝑧_𝑅²𝑥_𝑅−𝑥_𝑆𝑙𝜓,if𝑙=2,_𝑎𝑙=arctan𝑦_𝑆−𝑦_𝑅²+𝑧_𝑆−𝑧_𝑅²𝑥_𝑅−𝑥_𝑆𝑙. (7) Note that the calculation of these angles is only useful if either the source or receiver is directional. The time of arrival for each image is given by 𝜏_𝑙=𝑅_𝑙𝑐. (8)

Figure 3 

Three-dimensional geometry of image method with a rigid boundary and an open-end.

2.2. Method of Images Applied between Boundaries Π₃ and Π₄

We now assume that Π₃ and Π₄ are infinite rigid boundaries, so that the pressure gradient along 𝑛₃ and 𝑛₄ is null, (𝜕𝑃/𝜕𝑦)_𝑦=0=(𝜕𝑃/𝜕𝑦)_𝑦=𝐿=0. In this case, an infinite number of images are modeled. The images are grouped by four, where the very first group contains 𝑆₁, 𝑆₂, and 𝑆₃ and the image of 𝑆₂ across Π₃, noted 𝑆₄, as shown in Figure 4. The index 𝑚 corresponds to the group number containing the source 𝑆_𝑚𝑎. The second index (𝑎=1,2,3,or4) corresponds to the image number in each group 𝑚. The total pressure field is 𝑃𝑥_𝑆,𝑦_𝑆,𝑧_𝑆,𝑥_𝑅,𝑦_𝑅,𝑧_𝑅=𝑃₀^∞_𝑚=01𝑅_𝑚1𝑒^{𝑖𝑘𝑅𝑚1}+1𝑅_𝑚2𝑒^{𝑖𝑘𝑅𝑚2}+1𝑅_𝑚3𝑒^{𝑖𝑘𝑅𝑚3}+1𝑅_𝑚4𝑒^{𝑖𝑘𝑅𝑚4}=𝑃₀^∞_𝑚=0⁴_𝑎=11𝑅_𝑚𝑎𝑒^{𝑖𝑘𝑅𝑚𝑎}. (9) Figure 5 shows the first eight images (𝑚=0,𝑚=1) in plane Π₇. Given a group number 𝑚, the distance between the receiver and the corresponding image is defined as 𝑅_𝑚1=𝑥_𝑆−𝑥_𝑅²+𝑧_𝑆−𝑧_𝑅²+2𝑚𝐿+𝑦_𝑆−𝑦_𝑅²,𝑅_𝑚2=𝑥_𝑆−𝑥_𝑅²+𝑧_𝑆−𝑧_𝑅²+2(𝑚+1)𝐿−𝑦_𝑆−𝑦_𝑅²,𝑅_𝑚3=𝑥_𝑆−𝑥_𝑅²+𝑧_𝑆−𝑧_𝑅²+2𝑚𝐿+𝑦_𝑆+𝑦_𝑅²,𝑅_𝑚4=𝑥_𝑆−𝑥_𝑅²+𝑧_𝑆−𝑧_𝑅²+2(𝑚+1)𝐿−𝑦_𝑆+𝑦_𝑅². (10) The angles of transmission and arrival are defined as 𝜙^𝑡_𝑚𝑎=−arctan𝑥_𝑆−𝑥_𝑅²+𝑧_𝑆−𝑧_𝑅²𝑦_𝑅−𝑦_𝑆𝑚𝑎,if𝑎=1,4,arctan𝑥_𝑆−𝑥_𝑅²+𝑧_𝑆−𝑧_𝑅²𝑦_𝑅−𝑦_𝑆𝑚𝑎𝜙,if𝑎=2,3,_𝑎𝑚𝑎=arctan𝑥_𝑆−𝑥_𝑅²+𝑧_𝑆−𝑧_𝑅²𝑦_𝑅−𝑦_𝑆𝑚𝑎. (11) Note that the calculation of these angles is only useful if either the source or receiver is directional. The times of arrival are given by 𝜏_𝑚𝑎=𝑅_𝑚𝑎𝑐. (12)

Figure 4 

Three-Dimensional Distribution of Images between Rigid Boundaries Π3 and Π4, Group m = 0.

Figure 5 

Two-Dimensional Distribution of Images between Rigid Boundaries Π3 and Π4, Groups m = 0 and m = 1.

2.3. Method of Images Applied between Boundaries Π₁ and Π₂

We now assume that Π₁ and Π₂ are infinite nonrigid boundaries, as shown in Figure 6. Π₁ is a pressure release boundary, where the pressure along the normal 𝑛₁ is null. The boundary Π₂ is a silt boundary with a mass density of 𝜌_silt=1500kg/m³ and a sound speed 𝑐_silt=0.985𝑐 relative to the sound speed in the channel 𝑐. The reflection coefficient 𝑉_Π2 is a function of the angle of incidence 𝜃_𝑛𝑏 with respect to 𝑛₂. 𝑆_𝑛𝑏 corresponds to the image number 𝑏 in each group 𝑛. Using these parameters, the reflection coefficient 𝑉_Π2(𝜃_𝑛𝑏) is obtained [23]: 𝜃^𝑡_𝑛𝑏=sin⁻¹𝑐1−^2silt𝑐²cos²𝜃_𝑛𝑏,𝑉_Π2𝜃_𝑛𝑏=𝜌_silt𝑐_silt/𝜌𝑐−sin𝜃^𝑡_𝑛𝑏/sin𝜃_𝑛𝑏𝜌_silt𝑐_silt/𝜌𝑐+sin𝜃^𝑡_𝑛𝑏/sin𝜃_𝑛𝑏. (13)𝜃^𝑡_𝑛𝑏 is defined with respect to the positive 𝑧-axis and corresponds to the angle of transmission from a given image location 𝑆_𝑛𝑏. The resulting pressure field is the sum of the source pressure field and the pressure fields corresponding to the images of the source 𝑆₀₁. Since Π₂ is a partially reflecting boundary, the pressure field of the image 𝑆₀₂ is multiplied by the reflection coefficient of the boundary. Summing the pressure field of the image source 𝑆₀₂ to that of the pressure field of the source 𝑆₀₁, we obtain 𝑃𝑥_𝑆,𝑦_𝑆,𝑧_𝑆,𝑥_𝑅,𝑦_𝑅,𝑧_𝑅=𝑃₀𝑅₀₁𝑒^{𝑖𝑘𝑅01}+𝑉_Π2𝜃₀₂𝑃₀𝑅₀₂𝑒^{𝑖𝑘𝑅02}, (14)𝑅₀₁=𝑥_𝑆−𝑥_𝑅²+𝑦_𝑆−𝑦_𝑅²+𝑧_𝑆−𝑧_𝑅²,𝑅₀₂=𝑥_𝑆−𝑥_𝑅²+𝑦_𝑆−𝑦_𝑅²+2𝐻−𝑧_𝑆−𝑧_𝑅². (15) Since Π₁ is a pressure release boundary, 𝑆₀₃ and 𝑆₀₄ are out of phase with the actual source. The combined acoustic pressure field due to these images is𝑃𝑥_𝑆,𝑦_𝑆,𝑧_𝑆,𝑥_𝑅,𝑦_𝑅,𝑧_𝑅=𝑃₀𝑅₀₁𝑒^{𝑖𝑘𝑅01}+𝑉_Π2𝜃₀₂𝑃₀𝑅₀₂𝑒^{𝑖𝑘𝑅02}−𝑃₀𝑅₀₃𝑒^{𝑖𝑘𝑅03}−𝑉_Π2𝜃₀₄𝑃₀𝑅₀₄𝑒^{𝑖𝑘𝑅04}. (16) The distance between the receiver and the corresponding image is𝑅₀₁=𝑥_𝑆−𝑥_𝑅²+𝑦_𝑆−𝑦_𝑅²+𝑧_𝑆−𝑧_𝑅²,𝑅₀₂=𝑥_𝑆−𝑥_𝑅²+𝑦_𝑆−𝑦_𝑅²+2𝐻−𝑧_𝑆−𝑧_𝑅²,𝑅₀₃=𝑥_𝑆−𝑥_𝑅²+𝑦_𝑆−𝑦_𝑅²+𝑧_𝑆+𝑧_𝑅²,𝑅₀₄=𝑥_𝑆−𝑥_𝑅²+𝑦_𝑆−𝑦_𝑅²+2𝐻−𝑧_𝑆+𝑧_𝑅². (17) Given a group number n, the total pressure field at the receiver is 𝑃𝑥_𝑆,𝑦_𝑆,𝑧_𝑆,𝑥_𝑅,𝑦_𝑅,𝑧_𝑅=𝑃₀^∞_𝑛=0(−1)^𝑛𝑉_Π2𝜃_𝑛1^𝑛1𝑅_𝑛1𝑒^{𝑖𝑘𝑅𝑛1}+𝑉_Π2𝜃_𝑛2^(𝑛+1)1𝑅_𝑛2𝑒^{𝑖𝑘𝑅𝑛2}−𝑉_Π2𝜃_𝑛3^𝑛1𝑅_𝑛3𝑒^{𝑖𝑘𝑅𝑛3}−𝑉_Π2𝜃_𝑛4^(𝑛+1)1𝑅_𝑛4𝑒^{𝑖𝑘𝑅𝑛4}. (18) The distance between the receiver and each image is𝑅_𝑛1=𝑥_𝑆−𝑥_𝑅²+𝑦_𝑆−𝑦_𝑅²+2𝑛𝐻+𝑧_𝑆−𝑧_𝑅²,𝑅_𝑛2=𝑥_𝑆−𝑥_𝑅²+𝑦_𝑆−𝑦_𝑅²+2(𝑛+1)𝐻−𝑧_𝑆−𝑧_𝑅²,𝑅_𝑛3=𝑥_𝑆−𝑥_𝑅²+𝑦_𝑆−𝑦_𝑅²+2𝑛𝐻+𝑧_𝑆+𝑧_𝑅²,𝑅_𝑛4=𝑥_𝑆−𝑥_𝑅²+𝑦_𝑆−𝑦_𝑅²+2(𝑛+1)𝐻−𝑧_𝑆+𝑧_𝑅². (19) The angles of transmission and arrival are𝜃^𝑡_𝑛𝑏=−arctan𝑥_𝑆−𝑥_𝑅²+𝑦_𝑆−𝑦_𝑅²𝑧_𝑅−𝑧_𝑆𝑛𝑏,if𝑏=1,4,arctan𝑥_𝑆−𝑥_𝑅²+𝑦_𝑆−𝑦_𝑅²𝑧_𝑅−𝑧_𝑆𝑛𝑏𝜃,if𝑏=2,3,_𝑎𝑛𝑏=arctan𝑥_𝑆−𝑥_𝑅²+𝑦_𝑆−𝑦_𝑅²𝑧_𝑅−𝑧_𝑆𝑛𝑏. (20) Note that the calculation of these angles is only useful if either the source or receiver is directional. The times of arrival are given by 𝜏_𝑛𝑏=𝑅_𝑛𝑏𝑐. (21) Although the objective of this paper is not an in-depth study of the method of images, the reader should be aware that some approximations are made in (13), (14), (16), and (18). While the method of images is perfectly accurate for impenetrable surfaces, either hard (Neumann boundary condition) or soft (pressure release or Dirichlet boundary condition), it loses accuracy in the presence of penetrable surfaces. These issues are covered in detail in [17]. Here we assume that one of the boundaries is made of a uniform, soft sediment, so that sound travels more slowly in this sediment than in water.

Figure 6 

Three-Dimensional Distribution of Images between Pressure Release Boundary Π1 and Soft Boundary Π2, Group m = 0.

We use the geometrical acoustics approximation [17] to the total reflected field to obtain (14), (16), and (18): the energy of the acoustic field reflected off the soft boundary is concentrated about the angle of specular reflection. In reality, a spherical wave incident to the soft bottom would produce weaker levels of sound reflected in every direction: we assume that this secondary sound field is negligible. To further reduce the complexity of the model, we assume that the soft bottom is perfectly flat, uniform, and infinitely deep, so that any possible scattering and diffraction effect is neglected.

The main consequence of these approximations is that the model will overestimate the strength of the reflected sound in the specular direction and omit secondary echoes. In other terms, the modeled channel response will contain strong echoes at times directly related to the specular angles, but will not contain any echo related to wave curvature, scattering, and diffraction effects.

2.4. Pressure Field in the Duct

The final step is to find the analytical expression for the pressure field in the duct. We use the following index notation for each source 𝑆_{𝑙𝑚𝑎𝑛𝑏}: 𝑙 corresponds to the image number in the 𝑥-direction (𝑙=1,2), 𝑚 and 𝑎 correspond to the group number (𝑚=0,1,2,…,∞) and the image number (𝑎=1,2,3,4) in the y-direction, and 𝑛 and 𝑏 correspond to the group number (𝑛=0,1,2,…,∞) and the image number (𝑏=1,2,3,4) in the z-direction. Combining (5), (9), and (18), the acoustic pressure at the receiver is obtained:𝑃𝑥_𝑆,𝑦_𝑆,𝑧_𝑆,𝑥_𝑅,𝑦_𝑅,𝑧_𝑅=𝑃₀²_𝑙=1^∞_𝑚=0⁴_𝑎=1^∞_𝑛=0(−1)^𝑛𝑉_Π2𝜃_{𝑙𝑚𝑎𝑛1}^𝑛1𝑅_{𝑙𝑚𝑎𝑛1}𝑒^{𝑖𝑘𝑅𝑙𝑚𝑎𝑛1}+𝑉_Π2𝜃_{𝑙𝑚𝑎𝑛2}^(𝑛+1)1𝑅_{𝑙𝑚𝑎𝑛2}𝑒^{𝑖𝑘𝑅𝑙𝑚𝑎𝑛2}−𝑉_Π2𝜃_{𝑙𝑚𝑎𝑛3}^𝑛1𝑅_{𝑙𝑚𝑎𝑛3}𝑒^{𝑖𝑘𝑅𝑙𝑚𝑎𝑛3}−𝑉_Π2𝜃_{𝑙𝑚𝑎𝑛4}^(𝑛+1)1𝑅_{𝑙𝑚𝑎𝑛4}𝑒^{𝑖𝑘𝑅𝑙𝑚𝑎𝑛4}. (22)The summations in (22) account for the reflected paths produced by every wall in the duct. The distance 𝑅_{𝑙𝑚𝑎𝑛𝑏} between the image 𝑆_{𝑙𝑚𝑎𝑛𝑏} and the receiver is given in (23):𝑅_{𝑙𝑚1𝑛1}=𝑥_𝑆+(−1)^𝑙𝑥_𝑅²+2𝑚𝐿+𝑦_𝑆−𝑦_𝑅²+2𝑛𝐻+𝑧_𝑆−𝑧_𝑅²,𝑅_{𝑙𝑚1𝑛2}=𝑥_𝑆+(−1)^𝑙𝑥_𝑅²+2𝑚𝐿+𝑦_𝑆−𝑦_𝑅²+2(𝑛+1)𝐻−𝑧_𝑆−𝑧_𝑅²,𝑅_{𝑙𝑚1𝑛3}=𝑥_𝑆+(−1)^𝑙𝑥_𝑅²+2𝑚𝐿+𝑦_𝑆−𝑦_𝑅²+2𝑛𝐻+𝑧_𝑆+𝑧_𝑅²,𝑅_{𝑙𝑚1𝑛4}=𝑥_𝑆+(−1)^𝑙𝑥_𝑅²+2𝑚𝐿+𝑦_𝑆−𝑦_𝑅²+2(𝑛+1)𝐻−𝑧_𝑆+𝑧_𝑅²,𝑅_{𝑙𝑚2𝑛1}=𝑥_𝑆+(−1)^𝑙𝑥_𝑅²+2(𝑚+1)𝐿−𝑦_𝑆−𝑦_𝑅²+2𝑛𝐻+𝑧_𝑆−𝑧_𝑅²,𝑅_{𝑙𝑚2𝑛2}=𝑥_𝑆+(−1)^𝑙𝑥_𝑅²+2(𝑚+1)𝐿−𝑦_𝑆−𝑦_𝑅²+2(𝑛+1)𝐻−𝑧_𝑆−𝑧_𝑅²,𝑅_{𝑙𝑚2𝑛3}=𝑥_𝑆+(−1)^𝑙𝑥_𝑅²+2(𝑚+1)𝐿−𝑦_𝑆−𝑦_𝑅²+2𝑛𝐻+𝑧_𝑆+𝑧_𝑅²,𝑅_{𝑙𝑚2𝑛4}=𝑥_𝑆+(−1)^𝑙𝑥_𝑅²+2(𝑚+1)𝐿−𝑦_𝑆−𝑦_𝑅²+2(𝑛+1)𝐻−𝑧_𝑆+𝑧_𝑅²,𝑅_{𝑙𝑚3𝑛1}=𝑥_𝑆+(−1)^𝑙𝑥_𝑅²+2𝑚𝐿+𝑦_𝑆+𝑦_𝑅²+2𝑛𝐻+𝑧_𝑆−𝑧_𝑅²,𝑅_{𝑙𝑚3𝑛2}=𝑥_𝑆+(−1)^𝑙𝑥_𝑅²+2𝑚𝐿+𝑦_𝑆+𝑦_𝑅²+2(𝑛+1)𝐻−𝑧_𝑆−𝑧_𝑅²,𝑅_{𝑙𝑚3𝑛3}=𝑥_𝑆+(−1)^𝑙𝑥_𝑅²+2𝑚𝐿+𝑦_𝑆+𝑦_𝑅²+2𝑛𝐻+𝑧_𝑆+𝑧_𝑅²,𝑅_{𝑙𝑚3𝑛4}=𝑥_𝑆+(−1)^𝑙𝑥_𝑅²+2𝑚𝐿+𝑦_𝑆+𝑦_𝑅²+2(𝑛+1)𝐻−𝑧_𝑆+𝑧_𝑅²,𝑅_{𝑙𝑚4𝑛1}=𝑥_𝑆+(−1)^𝑙𝑥_𝑅²+2(𝑚+1)𝐿−𝑦_𝑆+𝑦_𝑅²+2𝑛𝐻+𝑧_𝑆−𝑧_𝑅²,𝑅_{𝑙𝑚4𝑛2}=𝑥_𝑆+(−1)^𝑙𝑥_𝑅²+2(𝑚+1)𝐿−𝑦_𝑆+𝑦_𝑅²+2(𝑛+1)𝐻−𝑧_𝑆−𝑧_𝑅²,𝑅_{𝑙𝑚4𝑛3}=𝑥_𝑆+(−1)^𝑙𝑥_𝑅²+2(𝑚+1)𝐿−𝑦_𝑆+𝑦_𝑅²+2𝑛𝐻+𝑧_𝑆+𝑧_𝑅²,𝑅_{𝑙𝑚4𝑛4}=𝑥_𝑆+(−1)^𝑙𝑥_𝑅²+2(𝑚+1)𝐿−𝑦_𝑆+𝑦_𝑅²+2(𝑛+1)𝐻−𝑧_𝑆+𝑧_𝑅². (23)For example, if 𝑙=1,2, 𝑚=𝑛=0, and 𝑎,𝑏=1,2,3,4, thirty-two images are created, as shown in Figure 7. In seawater, a certain amount of the acoustic energy of the propagating signal is also lost in heat originating from viscosity and thermal relaxation. This absorption of sound represents a true loss of acoustic energy within the acoustic channel of propagation. The value of the sound absorption coefficient 𝛼_dB is given by Schulkin and Marsh [24]: 𝛼_dB=0.0170𝑆𝑓_𝑇𝑓²_𝑐𝑓²_𝑇+𝑓²_𝑐𝑓+0.0245^2𝑐𝑓_𝑇dB/km. (24)𝑆 is the salinity in parts per thousand (ppt), 𝑓_𝑐 is the central operating frequency in Hertz, and 𝑓_𝑇 is the temperature-dependent relaxation frequency given by 𝑓_𝑇=21.9×10^{9−(1520/(𝑇+273))}Hz. (25)𝑇 is the temperature in degree centigrade. The absorption coefficient α in 1/m is 𝛼𝛼=^dB⋅ln(10)20⋅10³1/m. (26)

Figure 7 

Three-Dimensional Plot of a Block of Thirty-Two Images.

To compute the final acoustic pressure field, we assume that the absorption coefficient (in salt water) per unit of wavelength is much smaller than the wave number, 𝛼≪𝑘. Using this assumption, we can apply the absorption coefficient to each image individually. Therefore, the final acoustic pressure field in the three-dimensional environment is𝑃𝑥_𝑆,𝑦_𝑆,𝑧_𝑆,𝑥_𝑅,𝑦_𝑅,𝑧_𝑅=𝑃₀²_𝑙=1^∞_𝑚=0⁴_𝑎=1^∞_𝑛=0(−1)^𝑛𝑉_Π2𝜃_{𝑙𝑚𝑎𝑛1}^𝑛1𝑅_{𝑙𝑚𝑎𝑛1}𝑒^{(𝑖𝑘−𝛼)𝑅𝑙𝑚𝑎𝑛1}+𝑉_Π2𝜃_{𝑙𝑚𝑎𝑛2}^(𝑛+1)1𝑅_{𝑙𝑚𝑎𝑛2}𝑒^{(𝑖𝑘−𝛼)𝑅𝑙𝑚𝑎𝑛2}−𝑉_Π2𝜃_{𝑙𝑚𝑎𝑛3}^𝑛1𝑅_{𝑙𝑚𝑎𝑛3}𝑒^{(𝑖𝑘−𝛼)𝑅𝑙𝑚𝑎𝑛3}−𝑉_Π2𝜃_{𝑙𝑚𝑎𝑛4}^(𝑛+1)1𝑅_{𝑙𝑚𝑎𝑛4}𝑒^{(𝑖𝑘−𝛼)𝑅𝑙𝑚𝑎𝑛4}. (27)The corresponding distances 𝑅_{𝑙𝑚𝑎𝑛𝑏} between the image 𝑆_{𝑙𝑚𝑎𝑛𝑏} and the receiver are defined in (23).

2.5. In-Band Channel Model

The acoustic model must accommodate broadband chirp signal transmission. One of two approaches is available for broadband analysis: (a) each frequency component of the source signal is analyzed individually to create the modeled transfer function of the acoustic channel, a computer intensive approach; (b) the echo location, magnitude, and phase are computed at the carrier frequency 𝑓_𝑐 and applied to a short, band-limited impulse of carrier frequency 𝑓_𝑐. Although this second approach is less accurate, it is retained as it is far less computer intensive. The simulated transmitted signal 𝑞(𝑡) consists of a broadband chirp of bandwidth 𝑊=18kHz, central frequency 𝑓_𝑐=24kHz, and duration 𝑇_𝑠=13.54ms.

Although a complete analysis of using (b) rather than (a) is beyond the scope of this paper, the impact of such an approximation should be briefly discussed. The major impact of using (b) is a distortion of the echo envelope (including the peak value), while the peak location for each echo remains fairly accurate. In addition, since sound absorption increases by 5 dB/km between 15 kHz and 33 kHz, the approximation in (b) leads to a small over-estimation of the pressure field above 24 kHz and inversely to an underestimation of the pressure field below 24 kHz. This distortion becomes more severe as the traveled distance associated with a specific echo increases.

The chirp is frequency modulated (linear sweep), and the envelope is a Blackman time window 𝑏_mdl(𝑡):𝑞(𝑡)=𝑏_mdl(𝑡)⋅sin2𝜋𝑓_𝑐+𝑊⋅𝑡2𝑇_𝑠𝑡𝜇𝑃𝑎, (28)𝑏_mdl(𝑡)𝑡=0.42−0.5cos2𝜋𝑇_𝑠𝑡+0.08cos4𝜋𝑇_𝑠. (29) The transmitted signal is expressed in analytical form using a Hilbert transform operation 𝐻{}:𝑞_mdl(𝑡)=𝑞_mdl(𝑡)+𝑖𝐻𝑞_mdl(𝑡)=𝑞_mdl(𝑡)+^∞_−∞𝑖𝑞𝜋𝑡′^mdl𝑡−𝑡′𝑑𝑡′𝜇Pa. (30) Next, the normalized autocorrelation of this complex chirp 𝑞_mdl(𝑡) is computed to produce a band-limited impulse: 𝑟̃_mdl(𝑡)=1/2𝑇_𝑠∫^𝑇𝑠_−𝑇𝑠𝑞_mdl(𝑡)𝑞^∗_mdl(𝑡+𝜏)𝑑𝜏1/2𝑇_𝑠∫^𝑇𝑠_−𝑇𝑠𝑞_mdl(𝑡)²𝑑𝑡. (31) The modeled source signal (complex acoustic pressure in μPa) is scaled using the source level SL:𝑠̃_mdl(𝑡)=10^SL/20⋅𝑟̃_mdl(𝑡). (32) Using the approximation aforementioned, the impulse response ℎ(𝑡) of the duct is calculated for each source and receiver position:ℎ(𝑡)=²_𝑙=1^∞_𝑚=0⁴_𝑎=1^∞_𝑛=0(−1)^𝑛𝑉_Π2𝜃_{𝑙𝑚𝑎𝑛1}^𝑛𝛿𝑡−𝜏_{𝑙𝑚𝑎𝑛1}𝑅_{𝑙𝑚𝑎𝑛1}𝑒^{(𝑖𝑘−𝛼)𝑅𝑙𝑚𝑎𝑛1}+𝑉_Π2𝜃_{𝑙𝑚𝑎𝑛2}^(𝑛+1)𝛿𝑡−𝜏_{𝑙𝑚𝑎𝑛2}𝑅_{𝑙𝑚𝑎𝑛2}𝑒^{(𝑖𝑘−𝛼)𝑅𝑙𝑚𝑎𝑛2}−𝑉_Π2𝜃_{𝑙𝑚𝑎𝑛3}^𝑛𝛿𝑡−𝜏_{𝑙𝑚𝑎𝑛3}𝑅_{𝑙𝑚𝑎𝑛3}𝑒^{(𝑖𝑘−𝛼)𝑅𝑙𝑚𝑎𝑛3}−𝑉_Π2𝜃_{𝑙𝑚𝑎𝑛4}^(𝑛+1)𝛿𝑡−𝜏_{𝑙𝑚𝑎𝑛4}𝑅_{𝑙𝑚𝑎𝑛4}𝑒^{(𝑖𝑘−𝛼)𝑅𝑙𝑚𝑎𝑛4}. (33)𝛿(𝑡−𝜏_{𝑙𝑚𝑎𝑛𝑏}) represents a Dirac delayed by 𝜏_{𝑙𝑚𝑎𝑛𝑏} seconds. The in-band impulse response is defined as the convolution between (32) and (33):ℎ_mdl(𝑡)=𝑠̃_mdl(𝑡)∗ℎ(𝑡)=²_𝑙=1^∞_𝑚=0⁴_𝑎=1^∞_𝑛=0(−1)^𝑛𝑉_Π2𝜃_{𝑙𝑚𝑎𝑛1}^𝑛𝑠̃_mdl𝑡−𝜏_{𝑙𝑚𝑎𝑛1}𝑅_{𝑙𝑚𝑎𝑛1}𝑒^{(𝑖𝑘−𝛼)𝑅𝑙𝑚𝑎𝑛1}+𝑉_Π2𝜃_{𝑙𝑚𝑎𝑛2}^(𝑛+1)𝑠̃_mdl𝑡−𝜏_{𝑙𝑚𝑎𝑛2}𝑅_{𝑙𝑚𝑎𝑛2}𝑒^{(𝑖𝑘−𝛼)𝑅𝑙𝑚𝑎𝑛2}−𝑉_Π2𝜃_{𝑙𝑚𝑎𝑛3}^𝑛𝑠̃_mdl𝑡−𝜏_{𝑙𝑚𝑎𝑛3}𝑅_{𝑙𝑚𝑎𝑛3}𝑒^{(𝑖𝑘−𝛼)𝑅𝑙𝑚𝑎𝑛3}−𝑉_Π2𝜃_{𝑙𝑚𝑎𝑛4}^(𝑛+1)𝑠̃_mdl𝑡−𝜏_{𝑙𝑚𝑎𝑛4}𝑅_{𝑙𝑚𝑎𝑛4}𝑒^{(𝑖𝑘−𝛼)𝑅𝑙𝑚𝑎𝑛4}. (34)ℎ_mdl is computed at fixed time interval 𝑛_𝜏Δ𝜏, where 𝑛_𝜏 is defined as𝑛_𝜏𝜏=int^{𝑙𝑚𝑎𝑛𝑏}𝐹_𝑆, (35) where Δ𝜏=1/𝐹_𝑆 and 𝜏_{𝑙𝑚𝑎𝑛𝑏}=𝑅_{𝑙𝑚𝑎𝑛𝑏}/𝑐. The function int() corresponds to the integer part of the result, 𝐹_𝑆=75600Hz is the sampling frequency, 𝜏_{𝑙𝑚𝑎𝑛𝑏} is the time of arrival of the echo produced by image 𝑆_{𝑙𝑚𝑎𝑛𝑏}, and 𝑐 is the sound speed in the channel.

3. Simulated and Experimental Results

3.1. Experimental Setup

A set of field experiments has been conducted on June 8, 2007, in the south turning basin of Port Everglades, Florida. The basin, shown in Figure 8, contains an unobstructed, flat bathymetry silt bottom surrounded with a vertical concrete wall to the south and vertically piled boulders on the west and north sides. The basin is 14 meters deep. The west wall is 255 meters long and the north wall is 290 meters long. Towards the east, the bottom slopes slowly upward towards the shore. The bottom on the east side is a relatively thick layer of mud and very fine silt. Because of the frequency of operations (15–33 kHz), the sound traveling east is for the most part absorbed by this thick layer of mud. Therefore, the acoustic channel can reasonably be approximated as a duct closed on the west end and open on the east end.

Figure 8 

Aerial View of GPS Locations for Each Recorded Sample.

The source is an FAU Gateway Buoy [25] configured as a pinger, equipped with an ITC-3460 transducer and Global Positioning System with Wide Area Augmentation System (GPS-WAAS) of ±1 meter accuracy. The source transducer is suspended in the water column at a fixed depth of 1.5 meters below the sea-surface and transmits the chirp signal given in (28). The ITC-6156 receiver transducer is mounted on an air-filled aluminum pressure vessel (6 inches in diameter, 30 inches in length) to simulate the acoustic shading produced by the underwater vehicle of interest. The source signal used in the acoustic model is the signal produced by the omnidirectional source transducer. Therefore, the frequency response of the transducer is taken into account in the model.

This pressure vessel is tethered to a small research vessel at a depth of 1.3 meters. The source is at a stationary location in the turning basin and the receiver produces a record of the impulse response within the basin at multiple locations. The receiver pressure vessel is parallel to the hull of the vessel, oriented broadside to the source at each record location. A total of 𝑁_𝑟 = 5 impulse response measurements are collected at each location.

Considering the large number of receiver locations tested, these locations have been grouped in separate regions, labeled 1 to 4 in Figure 8. Region 1 contains the data collected near the western boundary of the turning basin, composed of vertically piled boulders. Region 2 contains the data collected near the flat, concrete southern wall of the turning basin. Region 3 is the most distant from any walls and the deepest portion of the basin. The receiver is also located at fairly close range from the source in this region, so that the echoes originating from the walls are clearly separated from the surface and bottom bounces. Region four contains the data collected near the northern wall, composed of vertically piled boulders.

A series of vertical sound velocity measurements was performed at locations indicated in Figure 8. Figure 9 shows the average sound velocity versus depth for the six casts [26]. The velocity gradient Δ𝑐/Δ𝐻 in the turning basin is negative:Δ𝑐=𝑐Δ𝐻^max−𝑐_min=Δ𝐻1541.3−154014=0.092857s⁻¹. (36) Given the limited depth and range, we can reasonably approximate this channel as an iso-velocity channel. Figure 10 shows the actual noise power spectral density (PSD) recorded during the field experimentation [26]. The standard deviation of the in-band noise is 𝜎_𝑛=110148𝜇Pa between 15 kHz and 33 kHz.

Figure 9 

Average Sound Velocity Measured by the CTD.

Figure 10 

Ambient Noise PSD.

3.2. Echo Overlap

An issue in estimating the impulse response of an acoustic channel is the limited source bandwidth, which in turn limits the time resolution. Figure 11 shows the magnitude of the modeled impulse response for a receiver located at point 008, with a receiver depth of 1.3 meters and a source depth of 1.5 meters. Two echoes corresponding to the direct path and surface bounce are received between the observation time of 94 ms and 94.5 ms. In (31), the −6 dB pulse width for each echo is 172𝜇s, causing these two echo partially overlap in Figure 11. This in turn causes echo interference, so that the exact sample spacing and magnitude cannot be determined exactly. Consider two received (complex and band-limited) echoes in (34), noted as:𝑟₁̃(𝑡)=ℎ₁𝑒^𝑗𝜙1𝑠̃_mdl𝑡−𝜏₁,𝑟₂̃(𝑡)=ℎ₂𝑒^𝑗𝜙2𝑠̃_mdl𝑡−𝜏₂. (37) These echoes may partly overlap and interfere (as depicted in Figure 12), so that the actual peak location for each echo cannot be exactly known. This results in inaccuracies when estimating the time-of-arrival of a given echo. Each echo phase 𝜙 is the sum of the phase shift 𝜙_boundary due to boundary interactions and the phase shift 𝜙_travel due to the relative movement between the source and the receiver: 𝜙=𝜙_boundary+𝜙_travel,𝜙_travel=2𝜋𝑓_𝑐Δ𝜏. (38)Δ𝜏=Δ𝑙/𝑐 is the time shift due to a receiver movement (Δ𝑙) in meters, and 𝜆=𝑐/𝑓_𝑐 is the wavelength of the transmitted signal. In our case, the center frequency for the experimentation is 𝑓_𝑐=24kHz and the measured sound speed within the channel is 𝑐=1540m/s. A change in distance of a quarter of a wavelength (0.015 m) results in a change of π/2 radians in the relative phase between the two echoes. Overall, if |ℎ₁|=|ℎ₂|, |𝑟₁̃(𝑡)+𝑟₂̃(𝑡)| can vary between 0 and 2|𝑟₁̃(𝑡)| due to small fluctuations in the receiver location. Therefore, the limited time resolution of the experimental setup results in inaccurate measurements of the echo magnitude and location if interference occurs. This in turn means that clearly separated echoes are measured more precisely in terms of both magnitude and time of arrival. Because of the limited time resolution, two or more echoes arriving within a time window of two −6 dB pulse width (344 μs) are treated as a group of overlapping echoes in the field result analysis.

Figure 11 

Echo Group 1 of Example Modeled Impulse Response.

Figure 12 

Overlapping Echoes associated with the Direct Path and the first Sea Surface Bounce.

3.3. Model Sensitivity to Geometrical and Physical Inaccuracies

For each receiver location, one record contains 𝑁_𝑟=5 impulse responses measured at fixed intervals 𝑇_𝑟=4 seconds. The time origin 𝑡=0 of the record is the estimated time of arrival of the very first echo of the first impulse response. The modeled echo, corresponding to the direct path of the first record, is also located at 𝑡=0. The modeled impulse response is repeated every 4 seconds. For each record 𝑛_𝑟(𝑛_𝑟=1,…,5) and location 𝑛_𝑙 within region 𝑛_reg(𝑛_reg=1,…,4), 𝜏_mdl(𝑛_𝑟,𝑛_𝑔,𝑛_𝑙,𝑛_reg) represents the time delay of the modeled echo group 𝑛_𝑔 from the time-of-arrival of the modeled direct path. Similarly, 𝜏_exp(𝑛_𝑟,𝑛_𝑔,𝑛_𝑙,𝑛_reg) represents the estimated time delay of the modeled echo group 𝑛_𝑔 from the time of arrival of the same modeled direct path, including the measured direct path (𝑛_𝑔=1). Each measured echo location is that of the largest echo |ℎ_exp(𝑛_𝑟,𝑛_𝑔,𝑛_𝑙,𝑛_reg)| within a 9 ms search window centered on the corresponding modeled echo |ℎ_mdl(𝑛_𝑟,𝑛_𝑔,𝑛_𝑙,𝑛_reg)|. The specific 9 ms duration is based on the source and receiver position accuracy and the approximated geometry of the duct. As an example, Figure 13 shows the measured and modeled impulse response at location 008 for the first measured impulse response. Note that the location of each search window in Figure 13 is not exact. The relative errors in echo magnitude 𝜀_mag(𝑛_𝑟,𝑛_𝑔,𝑛_𝑙,𝑛_reg) and in echo location 𝜀_loc(𝑛_𝑟,𝑛_𝑔,𝑛_𝑙,𝑛_reg) between the model and the measurement are 𝜀_mag𝑛_𝑟,𝑛_𝑔,𝑛_𝑙,𝑛_reg=ℎ100×^exp𝑛_𝑟,𝑛_𝑔,𝑛_𝑙,𝑛_reg−ℎ_mdl𝑛_𝑟,𝑛_𝑔,𝑛_𝑙,𝑛_regℎ_mdl𝑛_𝑟,𝑛_𝑔,𝑛_𝑙,𝑛_reg𝜀%,_loc𝑛_𝑟,𝑛_𝑔,𝑛_𝑙,𝑛_reg=𝜏100×^exp𝑛_𝑟,𝑛_𝑔,𝑛_𝑙,𝑛_reg−𝜏_mdl𝑛_𝑟,𝑛_𝑔,𝑛_𝑙,𝑛_reg𝜏_mdl𝑛_𝑟,𝑛_𝑔,𝑛_𝑙,𝑛_reg%. (39) The mean and standard deviation of these errors, computed across all the records at a specific location, region and group, are 𝜀_mag𝑛_𝑔,𝑛_𝑙,𝑛_reg=1𝑁_𝑟^𝑁𝑟_𝑖=1𝜀_mag𝑛_𝑟,𝑛_𝑔,𝑛_𝑙,𝑛_reg𝜎%,_𝜀mag𝑛_𝑔,𝑛_𝑙,𝑛_reg=1𝑁_𝑟−1^𝑁𝑟_𝑖=1𝜀_mag𝑛_𝑟,𝑛_𝑔,𝑛_𝑙,𝑛_reg−𝜀_mag𝑛_𝑔,𝑛_𝑙,𝑛_reg²%,𝜀_loc𝑛_𝑔,𝑛_𝑙,𝑛_reg=1𝑁_𝑟^𝑁𝑟_𝑖=1𝜀_loc𝑛_𝑟,𝑛_𝑔,𝑛_𝑙,𝑛_reg𝜎%,_𝜀loc𝑛_𝑔,𝑛_𝑙,𝑛_reg=1𝑁_𝑟−1^𝑁𝑟_𝑖=1𝜀_loc𝑛_𝑟,𝑛_𝑔,𝑛_𝑙,𝑛_reg−𝜀_loc𝑛_𝑔,𝑛_𝑙,𝑛_reg²%. (40) Figures 14 to 17 show the mean and standard deviation of the relative error between the modeled impulse response and the measured impulse response. Each numbered point corresponds to the mean value of the relative error for each echo group. The error bars in each plot represent two standard deviations of the relative error for each echo group. The 𝑥-axis corresponds to the relative error for the first arriving echo of each group. The relative error in magnitude for the first arriving echo of each group is displayed on the 𝑦-axis. A complete convergence analysis of the model as a function of the number of images is beyond the scope of this paper. However, the authors performed a limited study in [27], which indicated that the model converged using a small number of images.

Figure 13 

Numbered Echo Groups for Measured and Modeled Impulse Response at Location 008.

Figure 14 

Relative Error in Magnitude and Location for Region 1.

Figure 14 shows the relative error for region 1. This region corresponds to the receiver locations close to the western wall of the turning basin. This wall is composed of boulders piled vertically in an irregular manner. A significant share of the error observed between the modeled and measured impulse response is attributed to the structure of this wall. Since the boulders create a boundary that is not perfectly flat, the transmitted signal will be scattered at the interaction with the boundary. This scattering differs for different locations even within this region. Across the groups, the mean error in echo location varies between 0.2% and 2.1%. The mean error in echo magnitude varies between 55% (−2.6 dB) and 145% (1.62 dB). The relative error in magnitude of the direct path (group 1) of 145% is due to a significant overlap between the direct path and the surface bounce, caused by the small difference in time of travel between the direct path and the first surface echo. This same observation is made in regions 2 and 4. The standard deviation of the magnitude error varies from 2% (group 7) to 120% (group 1). The standard deviation of the location error varies between 0.1% (group 1) and 1% (group 3).

Figure 15 shows the relative error for each of the eight groups of echoes in region 2. The best match between modeled and measured data is observed in this region, as the relative error in echo group location is less than 1% in most cases. This is due to the fact that the receiver is relatively close to the southern wall of the turning basin and further away from the source. This vertical wall is flat; therefore the flat rigid surface assumed for the model constitutes a fair approximation. This large vertical wall is also an excellent reflector and the scattering is minimal as compared to the rocky walls to the north and west. Echo group 3 produces the largest error in echo location (1.7%) as it includes reflections from the rocky western wall. Across the groups, the mean error in echo location varies between 0.05% and 1.7%. The mean error in echo magnitude varies between 65% (−1.87 dB) and 120% (0.79 dB). The standard deviation of the magnitude error varies from 20% (group 5) to 80% (group 2). The standard deviation of the location error varies between 0.1% (group 1) and 1% (group 3).

Figure 15 

Relative Error in Magnitude and Location for Region 2.

Figure 16 shows the relative error for each of the eight groups of echoes in region 3, where all of the receiver locations are relatively close to the source. The relative error in the direct path (group 1) magnitude is significantly lower for this region (55%) as compared with the other regions (115% to 145%). Since the source and receiver depths are constant, the reduced distance between the source and the receiver results in an increased difference in time of travel between the direct path and the first surface echo. As a result, the overlap between these two echoes is reduced and the model prediction for the magnitude of the direct path is more accurate. However, the mean error in echo location for the direct path is 1.7%, which is significantly larger than the corresponding error for regions 1, 2, and 4. This is easily explained by the fact that direct path and surface bounce occur within the same search window with very similar magnitudes, so that the peak location of group 𝑛_𝑔=1 is occasionally associated with the surface bounce. The overall accuracy of the model in region 3 is significantly better than that in the other regions. Across the groups, the mean error in echo location varies between 0.4% and 1.7%, while the mean error in echo magnitude varies between 30% (−5.23 dB) and 80% (−0.97 dB). The standard deviation of the magnitude error varies from 5% (group 6) to 30% (group 1). The standard deviation of the location error varies between 0.4% (group 4) and 1.8% (group 1).

Figure 16 

Relative Error in Magnitude and Location for Region 3.

Figure 17 

Relative Error in Magnitude and Location for Region 4.

Figure 17 shows the relative error for each of the eight groups of echoes in region 4. This region corresponds to the receiver locations that are close to the northern wall of the turning basin. This case is similar to region 1 since the nearest boundary is composed of vertically piled boulders. Across the groups, the mean error in echo location varies between 0.1% and 1.05%, while the mean error in echo magnitude varies between 70% (−1.55 dB) and 120% (0.79 dB). The standard deviation of the magnitude error does not exceed 20%, with the exception of the first group containing the direct path (120%). The standard deviation of the location error varies between 0.1% (group 1) and 1.1% (group 2).

4. Conclusion

At first glance, a mean relative error of up to 145% (1.62 dB) in estimating the echo magnitude appears significant. However, this is a very reasonable number from the standpoint of underwater acoustic communications in nonstationary fading channels. Indeed, fluctuations in echo magnitude of 10 dB are commonly observed in shallow waters [28]. The error in echo location is also very good overall, as it always remains within 4% of the measured echo location. These encouraging results do not mean that the approximations made in the present model (negligible effects of wave curvature, scattering, and diffraction) should be the norm: a more realistic model would be expected to produce more accurate results. Rather, the model may simply be an acceptable trade-off between accuracy and computational load for acoustic communication purposes. Since the performance estimation of underwater acoustic modems is measured in terms of probability of bit error averaged over a very large number of samples, the three-dimensional model presented provides a sufficient level of accuracy to be used in the simulation of an acoustic communication system operating between 15 kHz and 33 kHz, with the benefit of low computing requirements.

Acknowledgments

This work is sponsored by the Center for Coastline Security and Technology at Florida Atlantic University and the Office of Naval Research.