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

Authors: Pierre-Philippe J. Beaujean; Matthew D. Staska
Journal: Advances in Acoustics and Vibration
 (2012)

---

## 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.


---

## Body

## 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, Figure1 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 Figure1, the authors use the Nakagami model [7, 9]. If A2 represents the acoustic energy of each impulse within a message, mNakagami=E{A2}2/Var{A2} is the ratio of the squared expectation (statistical mean) of A2 to the variance of A2. Here, the expectation is estimated using the average across all the impulses contained in a message. E{A2} is the average energy of the message. Var{A2} is the energy spread. The parameter mNakagami 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 parametermNakagami 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 mNakagami.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 parametermNakagami. 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 parametermNakagami 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 Figure2, where the water mass density ρ and sound velocity c are constant. For now, we also assume that S1(xS,yS,zS) is a point source producing a complex, harmonic spherical pressure wave p̃ measured at the receiver location R(xR,yR,zR). This complex, harmonic pressure wave is the solution to the Helmholtz equation [23]:
(1)Δp̃(xS,yS,zS,xR,yR,zR)+k2p̃(xS,yS,zS,xR,yR,zR)=0.
Expressed in Cartesian coordinates and in the complex domain, the free-field harmonic solution to (1) is
(2)p̃(xS,yS,zS,xR,yR,zR,t)=P0ei(k((xS-xR)2+(yS-yR)2+(zS-zR)2)-ωct)(xS-xR)2+(yS-yR)2+(zS-zR)2,
where k=ωc/c is the acoustic wave number, ωc=2πfc is the angular frequency of the transmitted signal, and P0 (in μPa·m) is the acoustic pressure times unity distance, measured at 1 m from S1. If we define R1 as the distance between the source S1 and the receiver as
(3)R1(xS,yS,zS,xR,yR,zR)=(xS-xR)2+(yS-yR)2+(zS-zR)2,
the complex pressure field produced by the point source becomes
(4)p̃(R1,t)=P̃(R1)e-iωct=P0eikR1R1e-iωct.
The geometry shown in Figure 2 contains a pressure release boundary Π1, a silt bottom Π2, and three rigid walls Π3,Π4, and Π5. The duct along the positive x-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Π5 and the Open-End

We first assume thatΠ5 is an infinite rigid boundary, as shown in Figure 3. S2 is located in plane Π6, which contains n̂5 and s1r⃗. Since Π5 is a rigid boundary, the pressure gradient along the normal n̂5 is null, so that (∂P̃/∂x)x=0=0. Based on the method of images [16], the pressure at the receiver is the sum of the pressure generated by the source S1 and the pressure generated by the image S2, located at equal and opposite distance from the boundary:
(5)P̃(xS,yS,zS,xR,yR,zR)=P̃(R1)+P̃(R2)=P0R1eikR1+P0R2eikR2.R1 is the distance between S1 and the receiver. R2 is the distance between S2 and the receiver:
(6)R1=(xS-xR)2+(yS-yR)2+(zS-zR)2R2=(xS+xR)2+(yS-yR)2+(zS-zR)2ψlt and ψla correspond to the angles of transmission and arrival unique to the source Sl, where l=1,2 is the image index. Here l=1 corresponds to the physical source S1:
(7)ψlt={-arctan((yS-yR)2+(zS-zR)2xR-xSl),ifl=1,arctan((yS-yR)2+(zS-zR)2xR-xSl),ifl=2,ψla=arctan((yS-yR)2+(zS-zR)2xR-xSl).
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)τl=Rlc.Figure 3
Three-dimensional geometry of image method with a rigid boundary and an open-end.


### 2.2. Method of Images Applied between BoundariesΠ3 and Π4

We now assume thatΠ3 and Π4 are infinite rigid boundaries, so that the pressure gradient along n̂3 and n̂4 is null, (∂P̃/∂y)y=0=(∂P̃/∂y)y=L=0. In this case, an infinite number of images are modeled. The images are grouped by four, where the very first group contains S1, S2, and S3 and the image of S2 across Π3, noted S4, as shown in Figure 4. The index m corresponds to the group number containing the source Sma. The second index (a=1,2,3,or4) corresponds to the image number in each group m. The total pressure field is
(9)P̃(xS,yS,zS,xR,yR,zR)=P0∑m=0∞[1Rm1eikRm1+1Rm2eikRm2+1Rm3eikRm3+1Rm4eikRm4]=P0∑m=0∞[∑a=14[1RmaeikRma]].
Figure 5 shows the first eight images (m=0,m=1) in plane Π7. Given a group number m, the distance between the receiver and the corresponding image is defined as
(10)Rm1=(xS-xR)2+(zS-zR)2+(2mL+yS-yR)2,Rm2=(xS-xR)2+(zS-zR)2+[2(m+1)L-yS-yR]2,Rm3=(xS-xR)2+(zS-zR)2+(2mL+yS+yR)2,Rm4=(xS-xR)2+(zS-zR)2+[2(m+1)L-yS+yR]2.
The angles of transmission and arrival are defined as
(11)ϕmat={-arctan((xS-xR)2+(zS-zR)2yR-ySma),ifa=1,4,arctan((xS-xR)2+(zS-zR)2yR-ySma),ifa=2,3,ϕmaa=arctan((xS-xR)2+(zS-zR)2yR-ySma).
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)τma=Rmac.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Π1 and Π2

We now assume thatΠ1 and Π2 are infinite nonrigid boundaries, as shown in Figure 6. Π1 is a pressure release boundary, where the pressure along the normal n̂1 is null. The boundary Π2 is a silt boundary with a mass density of ρsilt=1500kg/m3 and a sound speed csilt=0.985c relative to the sound speed in the channel c. The reflection coefficient VΠ2 is a function of the angle of incidence θnb with respect to n̂2. Snb corresponds to the image number b in each group n. Using these parameters, the reflection coefficient VΠ2(θnb) is obtained [23]:
(13)θnbt=sin-1(1-csilt2c2cos⁡2⁡(θnb)),VΠ2(θnb)=ρsiltcsilt/ρc-sin⁡(θnbt)/sin⁡(θnb)ρsiltcsilt/ρc+sin⁡(θnbt)/sin⁡(θnb).θnbt is defined with respect to the positive z-axis and corresponds to the angle of transmission from a given image location Snb. The resulting pressure field is the sum of the source pressure field and the pressure fields corresponding to the images of the source S01. Since Π2 is a partially reflecting boundary, the pressure field of the image S02 is multiplied by the reflection coefficient of the boundary. Summing the pressure field of the image source S02 to that of the pressure field of the source S01, we obtain
(14)P̃(xS,yS,zS,xR,yR,zR)=P0R01eikR01+VΠ2(θ02)P0R02eikR02,(15)R01=(xS-xR)2+(yS-yR)2+(zS-zR)2,R02=(xS-xR)2+(yS-yR)2+(2H-zS-zR)2.
Since Π1 is a pressure release boundary, S03 and S04 are out of phase with the actual source. The combined acoustic pressure field due to these images is(16)P̃(xS,yS,zS,xR,yR,zR)=P0R01eikR01+VΠ2(θ02)P0R02eikR02-P0R03eikR03-VΠ2(θ04)P0R04eikR04.
The distance between the receiver and the corresponding image is(17)R01=(xS-xR)2+(yS-yR)2+(zS-zR)2,R02=(xS-xR)2+(yS-yR)2+[2H-zS-zR]2,R03=(xS-xR)2+(yS-yR)2+(zS+zR)2,R04=(xS-xR)2+(yS-yR)2+[2H-zS+zR]2.
Given a group number n, the total pressure field at the receiver is
(18)P̃(xS,yS,zS,xR,yR,zR)=P0∑n=0∞(-1)n[(VΠ2(θn1))n1Rn1eikRn1+(VΠ2(θn2))(n+1)1Rn2eikRn2-(VΠ2(θn3))n1Rn3eikRn3-(VΠ2(θn4))(n+1)1Rn4eikRn4].
The distance between the receiver and each image is(19)Rn1=(xS-xR)2+(yS-yR)2+(2nH+zS-zR)2,Rn2=(xS-xR)2+(yS-yR)2+[2(n+1)H-zS-zR]2,Rn3=(xS-xR)2+(yS-yR)2+(2nH+zS+zR)2,Rn4=(xS-xR)2+(yS-yR)2+[2(n+1)H-zS+zR]2.
The angles of transmission and arrival are(20)θnbt={-arctan((xS-xR)2+(yS-yR)2zR-zSnb),ifb=1,4,arctan((xS-xR)2+(yS-yR)2zR-zSnb),ifb=2,3,θnba=arctan((xS-xR)2+(yS-yR)2zR-zSnb).
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)τnb=Rnbc.
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 sourceSlmanb: l corresponds to the image number in the x-direction (l=1,2), m and a correspond to the group number (m=0,1,2,…,∞) and the image number (a=1,2,3,4) in the y-direction, and n and b correspond to the group number (n=0,1,2,…,∞) and the image number (b=1,2,3,4) in the z-direction. Combining (5), (9), and (18), the acoustic pressure at the receiver is obtained:(22)P̃(xS,yS,zS,xR,yR,zR)=P0∑l=12[∑m=0∞[∑a=14[∑n=0∞(-1)n[(VΠ2(θlman1))n1Rlman1eikRlman1+(VΠ2(θlman2))(n+1)1Rlman2eikRlman2-(VΠ2(θlman3))n1Rlman3eikRlman3-(VΠ2(θlman4))(n+1)1Rlman4eikRlman4]]]].The summations in (22) account for the reflected paths produced by every wall in the duct. The distance Rlmanb between the image Slmanb and the receiver is given in (23):(23)Rlm1n1=[xS+(-1)lxR]2+(2mL+yS-yR)2+(2nH+zS-zR)2,Rlm1n2=[xS+(-1)lxR]2+(2mL+yS-yR)2+[2(n+1)H-zS-zR]2,Rlm1n3=[xS+(-1)lxR]2+(2mL+yS-yR)2+(2nH+zS+zR)2,Rlm1n4=[xS+(-1)lxR]2+(2mL+yS-yR)2+[2(n+1)H-zS+zR]2,Rlm2n1=[xS+(-1)lxR]2+[2(m+1)L-yS-yR]2+(2nH+zS-zR)2,Rlm2n2=[xS+(-1)lxR]2+[2(m+1)L-yS-yR]2+[2(n+1)H-zS-zR]2,Rlm2n3=[xS+(-1)lxR]2+[2(m+1)L-yS-yR]2+(2nH+zS+zR)2,Rlm2n4=[xS+(-1)lxR]2+[2(m+1)L-yS-yR]2+[2(n+1)H-zS+zR]2,Rlm3n1=[xS+(-1)lxR]2+(2mL+yS+yR)2+(2nH+zS-zR)2,Rlm3n2=[xS+(-1)lxR]2+(2mL+yS+yR)2+[2(n+1)H-zS-zR]2,Rlm3n3=[xS+(-1)lxR]2+(2mL+yS+yR)2+(2nH+zS+zR)2,Rlm3n4=[xS+(-1)lxR]2+(2mL+yS+yR)2+[2(n+1)H-zS+zR]2,Rlm4n1=[xS+(-1)lxR]2+[2(m+1)L-yS+yR]2+(2nH+zS-zR)2,Rlm4n2=[xS+(-1)lxR]2+[2(m+1)L-yS+yR]2+[2(n+1)H-zS-zR]2,Rlm4n3=[xS+(-1)lxR]2+[2(m+1)L-yS+yR]2+(2nH+zS+zR)2,Rlm4n4=[xS+(-1)lxR]2+[2(m+1)L-yS+yR]2+[2(n+1)H-zS+zR]2.For example, if l=1,2, m=n=0, and a,b=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]:
(24)αdB=0.0170SfTfc2fT2+fc2+0.0245fc2fTdB/km.S is the salinity in parts per thousand (ppt), fc is the central operating frequency in Hertz, and fT is the temperature-dependent relaxation frequency given by
(25)fT=21.9×109-(1520/(T+273))Hz.T is the temperature in degree centigrade. The absorption coefficient α in 1/m is
(26)α=αdB⋅ln⁡(10)20⋅1031/m.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,α≪k. 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(27)P̃(xS,yS,zS,xR,yR,zR)=P0∑l=12[∑m=0∞[∑a=14[∑n=0∞(-1)n[(VΠ2(θlman1))n1Rlman1e(ik-α)Rlman1+(VΠ2(θlman2))(n+1)1Rlman2e(ik-α)Rlman2-(VΠ2(θlman3))n1Rlman3e(ik-α)Rlman3-(VΠ2(θlman4))(n+1)1Rlman4e(ik-α)Rlman4]]]].The corresponding distances Rlmanb between the image Slmanb 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 frequencyfc and applied to a short, band-limited impulse of carrier frequency fc. Although this second approach is less accurate, it is retained as it is far less computer intensive. The simulated transmitted signal q(t) consists of a broadband chirp of bandwidth W=18kHz, central frequency fc=24kHz, and duration Ts=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 windowbmdl(t):(28)q(t)=bmdl(t)⋅sin⁡(2π(fc+W⋅t2Ts)t)μPa,(29)bmdl(t)=0.42-0.5cos⁡(2πtTs)+0.08cos⁡(4πtTs).
The transmitted signal is expressed in analytical form using a Hilbert transform operation H{}:(30)q̃mdl(t)=qmdl(t)+iH{qmdl(t)}=qmdl(t)+∫-∞∞iπt′qmdl(t-t′)dt′μPa.
Next, the normalized autocorrelation of this complex chirp q̃mdl(t) is computed to produce a band-limited impulse:
(31)r̃mdl(t)=[(1/2Ts)∫-TsTsq̃mdl(t)q̃mdl*(t+τ)dτ(1/2Ts)∫-TsTs|q̃mdl(t)|2dt].
The modeled source signal (complex acoustic pressure in μPa) is scaled using the source level SL:(32)s̃mdl(t)=10SL/20⋅r̃mdl(t).
Using the approximation aforementioned, the impulse response h̃(t) of the duct is calculated for each source and receiver position:(33)h̃(t)=∑l=12[∑m=0∞[∑a=14[∑n=0∞(-1)n[(VΠ2(θlman1))nδ(t-τlman1)Rlman1e(ik-α)Rlman1+(VΠ2(θlman2))(n+1)δ(t-τlman2)Rlman2e(ik-α)Rlman2-(VΠ2(θlman3))nδ(t-τlman3)Rlman3e(ik-α)Rlman3-(VΠ2(θlman4))(n+1)δ(t-τlman4)Rlman4e(ik-α)Rlman4]]]].δ(t-τlmanb) represents a Dirac delayed by τlmanb seconds. The in-band impulse response is defined as the convolution between (32) and (33):(34)h̃mdl(t)=s̃mdl(t)*h̃(t)=∑l=12[∑m=0∞[∑a=14[∑n=0∞(-1)n[(VΠ2(θlman1))ns̃mdl(t-τlman1)Rlman1e(ik-α)Rlman1+(VΠ2(θlman2))(n+1)s̃mdl(t-τlman2)Rlman2e(ik-α)Rlman2-(VΠ2(θlman3))ns̃mdl(t-τlman3)Rlman3e(ik-α)Rlman3-(VΠ2(θlman4))(n+1)s̃mdl(t-τlman4)Rlman4e(ik-α)Rlman4]]]].h̃mdl is computed at fixed time interval nτΔτ, where nτ is defined as(35)nτ=int⁡(τlmanbFS),
where Δτ=1/FS and τlmanb=Rlmanb/c. The function int() corresponds to the integer part of the result, FS=75600Hz is the sampling frequency, τlmanb is the time of arrival of the echo produced by image Slmanb, and c is the sound speed in the channel.


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

We first assume thatΠ5 is an infinite rigid boundary, as shown in Figure 3. S2 is located in plane Π6, which contains n̂5 and s1r⃗. Since Π5 is a rigid boundary, the pressure gradient along the normal n̂5 is null, so that (∂P̃/∂x)x=0=0. Based on the method of images [16], the pressure at the receiver is the sum of the pressure generated by the source S1 and the pressure generated by the image S2, located at equal and opposite distance from the boundary:
(5)P̃(xS,yS,zS,xR,yR,zR)=P̃(R1)+P̃(R2)=P0R1eikR1+P0R2eikR2.R1 is the distance between S1 and the receiver. R2 is the distance between S2 and the receiver:
(6)R1=(xS-xR)2+(yS-yR)2+(zS-zR)2R2=(xS+xR)2+(yS-yR)2+(zS-zR)2ψlt and ψla correspond to the angles of transmission and arrival unique to the source Sl, where l=1,2 is the image index. Here l=1 corresponds to the physical source S1:
(7)ψlt={-arctan((yS-yR)2+(zS-zR)2xR-xSl),ifl=1,arctan((yS-yR)2+(zS-zR)2xR-xSl),ifl=2,ψla=arctan((yS-yR)2+(zS-zR)2xR-xSl).
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)τl=Rlc.Figure 3
Three-dimensional geometry of image method with a rigid boundary and an open-end.


## 2.2. Method of Images Applied between BoundariesΠ3 and Π4

We now assume thatΠ3 and Π4 are infinite rigid boundaries, so that the pressure gradient along n̂3 and n̂4 is null, (∂P̃/∂y)y=0=(∂P̃/∂y)y=L=0. In this case, an infinite number of images are modeled. The images are grouped by four, where the very first group contains S1, S2, and S3 and the image of S2 across Π3, noted S4, as shown in Figure 4. The index m corresponds to the group number containing the source Sma. The second index (a=1,2,3,or4) corresponds to the image number in each group m. The total pressure field is
(9)P̃(xS,yS,zS,xR,yR,zR)=P0∑m=0∞[1Rm1eikRm1+1Rm2eikRm2+1Rm3eikRm3+1Rm4eikRm4]=P0∑m=0∞[∑a=14[1RmaeikRma]].
Figure 5 shows the first eight images (m=0,m=1) in plane Π7. Given a group number m, the distance between the receiver and the corresponding image is defined as
(10)Rm1=(xS-xR)2+(zS-zR)2+(2mL+yS-yR)2,Rm2=(xS-xR)2+(zS-zR)2+[2(m+1)L-yS-yR]2,Rm3=(xS-xR)2+(zS-zR)2+(2mL+yS+yR)2,Rm4=(xS-xR)2+(zS-zR)2+[2(m+1)L-yS+yR]2.
The angles of transmission and arrival are defined as
(11)ϕmat={-arctan((xS-xR)2+(zS-zR)2yR-ySma),ifa=1,4,arctan((xS-xR)2+(zS-zR)2yR-ySma),ifa=2,3,ϕmaa=arctan((xS-xR)2+(zS-zR)2yR-ySma).
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)τma=Rmac.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Π1 and Π2

We now assume thatΠ1 and Π2 are infinite nonrigid boundaries, as shown in Figure 6. Π1 is a pressure release boundary, where the pressure along the normal n̂1 is null. The boundary Π2 is a silt boundary with a mass density of ρsilt=1500kg/m3 and a sound speed csilt=0.985c relative to the sound speed in the channel c. The reflection coefficient VΠ2 is a function of the angle of incidence θnb with respect to n̂2. Snb corresponds to the image number b in each group n. Using these parameters, the reflection coefficient VΠ2(θnb) is obtained [23]:
(13)θnbt=sin-1(1-csilt2c2cos⁡2⁡(θnb)),VΠ2(θnb)=ρsiltcsilt/ρc-sin⁡(θnbt)/sin⁡(θnb)ρsiltcsilt/ρc+sin⁡(θnbt)/sin⁡(θnb).θnbt is defined with respect to the positive z-axis and corresponds to the angle of transmission from a given image location Snb. The resulting pressure field is the sum of the source pressure field and the pressure fields corresponding to the images of the source S01. Since Π2 is a partially reflecting boundary, the pressure field of the image S02 is multiplied by the reflection coefficient of the boundary. Summing the pressure field of the image source S02 to that of the pressure field of the source S01, we obtain
(14)P̃(xS,yS,zS,xR,yR,zR)=P0R01eikR01+VΠ2(θ02)P0R02eikR02,(15)R01=(xS-xR)2+(yS-yR)2+(zS-zR)2,R02=(xS-xR)2+(yS-yR)2+(2H-zS-zR)2.
Since Π1 is a pressure release boundary, S03 and S04 are out of phase with the actual source. The combined acoustic pressure field due to these images is(16)P̃(xS,yS,zS,xR,yR,zR)=P0R01eikR01+VΠ2(θ02)P0R02eikR02-P0R03eikR03-VΠ2(θ04)P0R04eikR04.
The distance between the receiver and the corresponding image is(17)R01=(xS-xR)2+(yS-yR)2+(zS-zR)2,R02=(xS-xR)2+(yS-yR)2+[2H-zS-zR]2,R03=(xS-xR)2+(yS-yR)2+(zS+zR)2,R04=(xS-xR)2+(yS-yR)2+[2H-zS+zR]2.
Given a group number n, the total pressure field at the receiver is
(18)P̃(xS,yS,zS,xR,yR,zR)=P0∑n=0∞(-1)n[(VΠ2(θn1))n1Rn1eikRn1+(VΠ2(θn2))(n+1)1Rn2eikRn2-(VΠ2(θn3))n1Rn3eikRn3-(VΠ2(θn4))(n+1)1Rn4eikRn4].
The distance between the receiver and each image is(19)Rn1=(xS-xR)2+(yS-yR)2+(2nH+zS-zR)2,Rn2=(xS-xR)2+(yS-yR)2+[2(n+1)H-zS-zR]2,Rn3=(xS-xR)2+(yS-yR)2+(2nH+zS+zR)2,Rn4=(xS-xR)2+(yS-yR)2+[2(n+1)H-zS+zR]2.
The angles of transmission and arrival are(20)θnbt={-arctan((xS-xR)2+(yS-yR)2zR-zSnb),ifb=1,4,arctan((xS-xR)2+(yS-yR)2zR-zSnb),ifb=2,3,θnba=arctan((xS-xR)2+(yS-yR)2zR-zSnb).
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)τnb=Rnbc.
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 sourceSlmanb: l corresponds to the image number in the x-direction (l=1,2), m and a correspond to the group number (m=0,1,2,…,∞) and the image number (a=1,2,3,4) in the y-direction, and n and b correspond to the group number (n=0,1,2,…,∞) and the image number (b=1,2,3,4) in the z-direction. Combining (5), (9), and (18), the acoustic pressure at the receiver is obtained:(22)P̃(xS,yS,zS,xR,yR,zR)=P0∑l=12[∑m=0∞[∑a=14[∑n=0∞(-1)n[(VΠ2(θlman1))n1Rlman1eikRlman1+(VΠ2(θlman2))(n+1)1Rlman2eikRlman2-(VΠ2(θlman3))n1Rlman3eikRlman3-(VΠ2(θlman4))(n+1)1Rlman4eikRlman4]]]].The summations in (22) account for the reflected paths produced by every wall in the duct. The distance Rlmanb between the image Slmanb and the receiver is given in (23):(23)Rlm1n1=[xS+(-1)lxR]2+(2mL+yS-yR)2+(2nH+zS-zR)2,Rlm1n2=[xS+(-1)lxR]2+(2mL+yS-yR)2+[2(n+1)H-zS-zR]2,Rlm1n3=[xS+(-1)lxR]2+(2mL+yS-yR)2+(2nH+zS+zR)2,Rlm1n4=[xS+(-1)lxR]2+(2mL+yS-yR)2+[2(n+1)H-zS+zR]2,Rlm2n1=[xS+(-1)lxR]2+[2(m+1)L-yS-yR]2+(2nH+zS-zR)2,Rlm2n2=[xS+(-1)lxR]2+[2(m+1)L-yS-yR]2+[2(n+1)H-zS-zR]2,Rlm2n3=[xS+(-1)lxR]2+[2(m+1)L-yS-yR]2+(2nH+zS+zR)2,Rlm2n4=[xS+(-1)lxR]2+[2(m+1)L-yS-yR]2+[2(n+1)H-zS+zR]2,Rlm3n1=[xS+(-1)lxR]2+(2mL+yS+yR)2+(2nH+zS-zR)2,Rlm3n2=[xS+(-1)lxR]2+(2mL+yS+yR)2+[2(n+1)H-zS-zR]2,Rlm3n3=[xS+(-1)lxR]2+(2mL+yS+yR)2+(2nH+zS+zR)2,Rlm3n4=[xS+(-1)lxR]2+(2mL+yS+yR)2+[2(n+1)H-zS+zR]2,Rlm4n1=[xS+(-1)lxR]2+[2(m+1)L-yS+yR]2+(2nH+zS-zR)2,Rlm4n2=[xS+(-1)lxR]2+[2(m+1)L-yS+yR]2+[2(n+1)H-zS-zR]2,Rlm4n3=[xS+(-1)lxR]2+[2(m+1)L-yS+yR]2+(2nH+zS+zR)2,Rlm4n4=[xS+(-1)lxR]2+[2(m+1)L-yS+yR]2+[2(n+1)H-zS+zR]2.For example, if l=1,2, m=n=0, and a,b=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]:
(24)αdB=0.0170SfTfc2fT2+fc2+0.0245fc2fTdB/km.S is the salinity in parts per thousand (ppt), fc is the central operating frequency in Hertz, and fT is the temperature-dependent relaxation frequency given by
(25)fT=21.9×109-(1520/(T+273))Hz.T is the temperature in degree centigrade. The absorption coefficient α in 1/m is
(26)α=αdB⋅ln⁡(10)20⋅1031/m.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,α≪k. 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(27)P̃(xS,yS,zS,xR,yR,zR)=P0∑l=12[∑m=0∞[∑a=14[∑n=0∞(-1)n[(VΠ2(θlman1))n1Rlman1e(ik-α)Rlman1+(VΠ2(θlman2))(n+1)1Rlman2e(ik-α)Rlman2-(VΠ2(θlman3))n1Rlman3e(ik-α)Rlman3-(VΠ2(θlman4))(n+1)1Rlman4e(ik-α)Rlman4]]]].The corresponding distances Rlmanb between the image Slmanb 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 frequencyfc and applied to a short, band-limited impulse of carrier frequency fc. Although this second approach is less accurate, it is retained as it is far less computer intensive. The simulated transmitted signal q(t) consists of a broadband chirp of bandwidth W=18kHz, central frequency fc=24kHz, and duration Ts=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 windowbmdl(t):(28)q(t)=bmdl(t)⋅sin⁡(2π(fc+W⋅t2Ts)t)μPa,(29)bmdl(t)=0.42-0.5cos⁡(2πtTs)+0.08cos⁡(4πtTs).
The transmitted signal is expressed in analytical form using a Hilbert transform operation H{}:(30)q̃mdl(t)=qmdl(t)+iH{qmdl(t)}=qmdl(t)+∫-∞∞iπt′qmdl(t-t′)dt′μPa.
Next, the normalized autocorrelation of this complex chirp q̃mdl(t) is computed to produce a band-limited impulse:
(31)r̃mdl(t)=[(1/2Ts)∫-TsTsq̃mdl(t)q̃mdl*(t+τ)dτ(1/2Ts)∫-TsTs|q̃mdl(t)|2dt].
The modeled source signal (complex acoustic pressure in μPa) is scaled using the source level SL:(32)s̃mdl(t)=10SL/20⋅r̃mdl(t).
Using the approximation aforementioned, the impulse response h̃(t) of the duct is calculated for each source and receiver position:(33)h̃(t)=∑l=12[∑m=0∞[∑a=14[∑n=0∞(-1)n[(VΠ2(θlman1))nδ(t-τlman1)Rlman1e(ik-α)Rlman1+(VΠ2(θlman2))(n+1)δ(t-τlman2)Rlman2e(ik-α)Rlman2-(VΠ2(θlman3))nδ(t-τlman3)Rlman3e(ik-α)Rlman3-(VΠ2(θlman4))(n+1)δ(t-τlman4)Rlman4e(ik-α)Rlman4]]]].δ(t-τlmanb) represents a Dirac delayed by τlmanb seconds. The in-band impulse response is defined as the convolution between (32) and (33):(34)h̃mdl(t)=s̃mdl(t)*h̃(t)=∑l=12[∑m=0∞[∑a=14[∑n=0∞(-1)n[(VΠ2(θlman1))ns̃mdl(t-τlman1)Rlman1e(ik-α)Rlman1+(VΠ2(θlman2))(n+1)s̃mdl(t-τlman2)Rlman2e(ik-α)Rlman2-(VΠ2(θlman3))ns̃mdl(t-τlman3)Rlman3e(ik-α)Rlman3-(VΠ2(θlman4))(n+1)s̃mdl(t-τlman4)Rlman4e(ik-α)Rlman4]]]].h̃mdl is computed at fixed time interval nτΔτ, where nτ is defined as(35)nτ=int⁡(τlmanbFS),
where Δτ=1/FS and τlmanb=Rlmanb/c. The function int() corresponds to the integer part of the result, FS=75600Hz is the sampling frequency, τlmanb is the time of arrival of the echo produced by image Slmanb, and c 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 Figure8, 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 ofNr = 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 Figure8. 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 Figure8. Figure 9 shows the average sound velocity versus depth for the six casts [26]. The velocity gradient Δc/ΔH in the turning basin is negative:(36)ΔcΔH=cmax⁡-cmin⁡ΔH=1541.3-154014=0.092857s-1.
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 σn=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. Figure11 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:(37)r̃1(t)=|h̃1|ejϕ1s̃mdl(t-τ1),r̃2(t)=|h̃2|ejϕ2s̃mdl(t-τ2).
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:
(38)ϕ=ϕboundary+ϕtravel,ϕtravel=2πfcΔτ.Δτ=Δl/c is the time shift due to a receiver movement (Δl) in meters, and λ=c/fc is the wavelength of the transmitted signal. In our case, the center frequency for the experimentation is fc=24kHz and the measured sound speed within the channel is c=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 |h̃1|=|h̃2|, |r̃1(t)+r̃2(t)| can vary between 0 and 2|r̃1(t)| 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 containsNr=5 impulse responses measured at fixed intervals Tr=4 seconds. The time origin t=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 t=0. The modeled impulse response is repeated every 4 seconds. For each record nr(nr=1,…,5) and location nl within region nreg(nreg=1,…,4), τmdl(nr,ng,nl,nreg) represents the time delay of the modeled echo group ng from the time-of-arrival of the modeled direct path. Similarly, τexp(nr,ng,nl,nreg) represents the estimated time delay of the modeled echo group ng from the time of arrival of the same modeled direct path, including the measured direct path (ng=1). Each measured echo location is that of the largest echo |hexp(nr,ng,nl,nreg)| within a 9 ms search window centered on the corresponding modeled echo |hmdl(nr,ng,nl,nreg)|. 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(nr,ng,nl,nreg) and in echo location εloc(nr,ng,nl,nreg) between the model and the measurement are
(39)εmag(nr,ng,nl,nreg)=(100×||h̃exp⁡(nr,ng,nl,nreg)|-|h̃mdl(nr,ng,nl,nreg)|||h̃mdl(nr,ng,nl,nreg)|)%,εloc⁡(nr,ng,nl,nreg)=(100×|τexp⁡(nr,ng,nl,nreg)-τmdl(nr,ng,nl,nreg)|τmdl(nr,ng,nl,nreg))%.
The mean and standard deviation of these errors, computed across all the records at a specific location, region and group, are
(40)εmag(ng,nl,nreg)¯=(1Nr∑i=1Nrεmag(nr,ng,nl,nreg))%,σεmag(ng,nl,nreg)=(1Nr-1∑i=1Nr(εmag(nr,ng,nl,nreg)-εmag(ng,nl,nreg)¯)2)%,εloc⁡(ng,nl,nreg)¯=(1Nr∑i=1Nrεloc⁡(nr,ng,nl,nreg))%,σεloc⁡(ng,nl,nreg)=(1Nr-1∑i=1Nr(εloc⁡(nr,ng,nl,nreg)-εloc⁡(ng,nl,nreg)¯)2)%.
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 x-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 y-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.Figure14 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).Figure15 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.Figure16 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 ng=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.Figure17 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).


## 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 Figure8, 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 ofNr = 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 Figure8. 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 Figure8. Figure 9 shows the average sound velocity versus depth for the six casts [26]. The velocity gradient Δc/ΔH in the turning basin is negative:(36)ΔcΔH=cmax⁡-cmin⁡ΔH=1541.3-154014=0.092857s-1.
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 σn=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. Figure11 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:(37)r̃1(t)=|h̃1|ejϕ1s̃mdl(t-τ1),r̃2(t)=|h̃2|ejϕ2s̃mdl(t-τ2).
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:
(38)ϕ=ϕboundary+ϕtravel,ϕtravel=2πfcΔτ.Δτ=Δl/c is the time shift due to a receiver movement (Δl) in meters, and λ=c/fc is the wavelength of the transmitted signal. In our case, the center frequency for the experimentation is fc=24kHz and the measured sound speed within the channel is c=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 |h̃1|=|h̃2|, |r̃1(t)+r̃2(t)| can vary between 0 and 2|r̃1(t)| 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 containsNr=5 impulse responses measured at fixed intervals Tr=4 seconds. The time origin t=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 t=0. The modeled impulse response is repeated every 4 seconds. For each record nr(nr=1,…,5) and location nl within region nreg(nreg=1,…,4), τmdl(nr,ng,nl,nreg) represents the time delay of the modeled echo group ng from the time-of-arrival of the modeled direct path. Similarly, τexp(nr,ng,nl,nreg) represents the estimated time delay of the modeled echo group ng from the time of arrival of the same modeled direct path, including the measured direct path (ng=1). Each measured echo location is that of the largest echo |hexp(nr,ng,nl,nreg)| within a 9 ms search window centered on the corresponding modeled echo |hmdl(nr,ng,nl,nreg)|. 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(nr,ng,nl,nreg) and in echo location εloc(nr,ng,nl,nreg) between the model and the measurement are
(39)εmag(nr,ng,nl,nreg)=(100×||h̃exp⁡(nr,ng,nl,nreg)|-|h̃mdl(nr,ng,nl,nreg)|||h̃mdl(nr,ng,nl,nreg)|)%,εloc⁡(nr,ng,nl,nreg)=(100×|τexp⁡(nr,ng,nl,nreg)-τmdl(nr,ng,nl,nreg)|τmdl(nr,ng,nl,nreg))%.
The mean and standard deviation of these errors, computed across all the records at a specific location, region and group, are
(40)εmag(ng,nl,nreg)¯=(1Nr∑i=1Nrεmag(nr,ng,nl,nreg))%,σεmag(ng,nl,nreg)=(1Nr-1∑i=1Nr(εmag(nr,ng,nl,nreg)-εmag(ng,nl,nreg)¯)2)%,εloc⁡(ng,nl,nreg)¯=(1Nr∑i=1Nrεloc⁡(nr,ng,nl,nreg))%,σεloc⁡(ng,nl,nreg)=(1Nr-1∑i=1Nr(εloc⁡(nr,ng,nl,nreg)-εloc⁡(ng,nl,nreg)¯)2)%.
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 x-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 y-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.Figure14 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).Figure15 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.Figure16 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 ng=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.Figure17 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.



---

*Source: 929174-2012-06-03.xml*