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

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.


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][2][3][4][5][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.
Acoustic modems transmit series of acoustic (bandlimited) impulses, each containing some binary information, to relay messages between the source and the receiver [1][2][3][4][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, signalto-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 A 2 represents the acoustic energy of each impulse within a message, m Nakagami = E{A 2 } 2 / Var{A 2 } is the ratio of the squared expectation (statistical mean) of A 2 to the variance of A 2 .Here, the expectation is estimated using the average across all the impulses contained in a message.E{A 2 } is the average energy of the message.Var{A 2 } is the energy spread.The parameter m 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 m 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 m 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 m 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 m 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 threedimensional 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][13][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 computerefficient 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][19][20][21][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 ultrasound 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.

Acoustic Model
Consider the modeled duct shown in Figure 2, where the water mass density ρ and sound velocity c are constant.For now, we also assume that S 1 (x S , y S , z S ) is a point source producing a complex, harmonic spherical pressure wave p measured at the receiver location R(x R , y R , z R ).This complex, harmonic pressure wave is the solution to the Helmholtz equation [23]: Expressed in Cartesian coordinates and in the complex domain, the free-field harmonic solution to (1) is where k = ω c /c is the acoustic wave number, ω c = 2π f c is the angular frequency of the transmitted signal, and P 0 (in μPa•m) is the acoustic pressure times unity distance, measured at 1 m from S 1 .If we define R 1 as the distance between the source S 1 and the receiver as the complex pressure field produced by the point source becomes 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.

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. S 2 is located in plane Π 6 , which contains n 5 and s 1 r.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 S 1 and the pressure generated by the image S 2 , located at equal and opposite distance from the boundary: ( Advances in Acoustics and Vibration Figure 3: Three-dimensional geometry of image method with a rigid boundary and an open-end.R 1 is the distance between S 1 and the receiver.R 2 is the distance between S 2 and the receiver: ψ t l and ψ a l correspond to the angles of transmission and arrival unique to the source S l , where l = 1, 2 is the image index.Here l = 1 corresponds to the physical source S 1 : 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

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 S 1 , S 2 , and S 3 and the image of S 2 across Π 3 , noted S 4 , as shown in Figure 4.The index m corresponds to the group number containing the source S ma .The second index (a = 1, 2, 3, or 4) corresponds to the image number in each group m.The total pressure field is 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 The angles of transmission and arrival are defined as 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

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 Advances in Acoustics and Vibration boundary, where the pressure along the normal n 1 is null.The boundary Π 2 is a silt boundary with a mass density of ρ silt = 1500 kg/m 3 and a sound speed c silt = 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 .S nb corresponds to the image number b in each group n.Using these parameters, the reflection coefficient V Π2 (θ nb ) is obtained [23]: θ t nb is defined with respect to the positive z-axis and corresponds to the angle of transmission from a given image location S nb .The resulting pressure field is the sum of the source pressure field and the pressure fields corresponding to the images of the source S 01 .Since Π 2 is a partially reflecting boundary, the pressure field of the image S 02 is multiplied by the reflection coefficient of the boundary.Summing the pressure field of the image source S 02 to that of the pressure field of the source S 01 , we obtain Since Π 1 is a pressure release boundary, S 03 and S 04 are out of phase with the actual source.The combined acoustic pressure field due to these images is Advances in Acoustics and Vibration 7 The distance between the receiver and the corresponding image is Given a group number n, the total pressure field at the receiver is The distance between the receiver and each image is The angles of transmission and arrival are 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 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.
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 S lmanb : 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: The summations in (22) account for the reflected paths produced by every wall in the duct.The distance R lmanb between the image S lmanb and the receiver is given in (23): 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]: S is the salinity in parts per thousand (ppt), f c is the central operating frequency in Hertz, and f T is the temperaturedependent relaxation frequency given by f T = 21.9 × 10 9−(1520/(T+273)) Hz.
T is the temperature in degree centigrade.The absorption coefficient α in 1/m is 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 The corresponding distances R lmanb between the image S lmanb and the receiver are defined in (23).

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 f c and applied to a short, bandlimited impulse of carrier frequency f c .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 = 18 kHz, central frequency f c = 24 kHz, and duration T s = 13.54 ms.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 overestimation 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 b mdl (t): The transmitted signal is expressed in analytical form using a Hilbert transform operation H{}: Next, the normalized autocorrelation of this complex chirp q mdl (t) is computed to produce a band-limited impulse: The modeled source signal (complex acoustic pressure in μPa) is scaled using the source level SL: Using the approximation aforementioned, the impulse response h(t) of the duct is calculated for each source and receiver position: Advances in Acoustics and Vibration δ(t − τ lmanb ) represents a Dirac delayed by τ lmanb seconds.The in-band impulse response is defined as the convolution between (32) and (33): h mdl is computed at fixed time interval n τ Δτ, where n τ is defined as where Δτ = 1/F S and τ lmanb = R lmanb /c.The function int() corresponds to the integer part of the result, F S = 75600 Hz is the sampling frequency, τ lmanb is the time of arrival of the echo produced by image S lmanb , and c is the sound speed in the channel.

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.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 N r = 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 Δc/ΔH in the turning basin is negative: 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.

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: 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: Δτ = Δl/c is the time shift due to a receiver movement (Δl) in meters, and λ = c/ f c is the wavelength of the transmitted signal.In our case, the center frequency for the experimentation is f c = 24 kHz and the measured sound speed within the channel is c = 1540 m/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  delay of the modeled echo group n g from the time-of-arrival of the modeled direct path.Similarly, τ exp (n r , n g , n l , n reg ) represents the estimated time delay of the modeled echo group n g from the time of arrival of the same modeled direct path, including the measured direct path (n g = 1).
Each measured echo location is that of the largest echo |h exp (n r , n g , n l , n reg )| within a 9 ms search window centered on the corresponding modeled echo |h mdl (n r , n g , n l , n 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 (n r , n g , n l , n reg ) and in echo location ε loc (n r , n g , n l , n reg ) between the model and the measurement are The mean and standard deviation of these errors, computed across all the records at a specific location, region and group, are 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 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 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 n g = 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%, 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 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).

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

Figure 7 :
Figure 7: Three-Dimensional Plot of a Block of Thirty-Two Images.

Figure 8 :Figure 9 :
Figure 8: Aerial View of GPS Locations for Each Recorded Sample.

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

Figure 13 :Figure 14 :
Figure 13: Numbered Echo Groups for Measured and Modeled Impulse Response at Location 008.

Figure 15 :
Figure 15: Relative Error in Magnitude and Location for Region 2.

Figure 16 :
Figure 16: Relative Error in Magnitude and Location for Region 3.
of echo group location

Figure 17 :
Figure 17: Relative Error in Magnitude and Location for Region 4.
For each receiver location, one record contains N r = 5 impulse responses measured at fixed intervals T r = 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 n r (n r = 1, . . ., 5) and location n l within region n reg (n reg = 1, . . ., 4), τ mdl (n r , n g , n l , n reg ) represents the time