Proceedings of the Institute of Acoustics

 

Synthesizing hydrophone data for applications in active sonar

 

S. T. H. Hartman, Norwegian Defence Research Establishment (FFI), PO Box 115, 3191 Horten, Norway

 

 

1 INTRODUCTION

 

The access to reliable high-quality data is often an important prerequisite for developing and testing novel processing techniques, and active sonar is no exception. However, carrying out at-sea experiments to obtain such datasets is an expensive and labor-intensive undertaking. Synthesizing data can therefore be a time and cost-effective alternative, with the added benefit that the data needs little pre-processing before it can be used. This is particularly useful in applications where large datasets need to be categorized, and has the added effect that all parameters of the dataset are known, which is not the case for real data. This can be utilized in many ways, for instance controlling the effect of certain parameters on the data, or constructing data of varying complexity. In this work we present a new tool under development for synthesizing realistic active sonar data at the hydrophone level called LASSIE (Long-range Active Sonar Simulations In the Environment). The theory and methods used, which include ray- and beam-tracing, eigenray analysis, modelling of reflections and scatterings, and the construction of the final synthetic hydrophone signals, are outlined in Section 2. As an example application of LASSIE, synthetic sonar data is generated and used as test data for a Markov Chain Monte Carlo (MCMC) sampling algorithm for extracting environmental information from the measured reverberation level, described in Section 3. In Section 4 we make concluding remarks about the current status of LASSIE and the results of the MCMC simulations.

 

2 SYNTHESIZING ACTIVE SONAR DATA

 

Generating realistic sonar data requires accurate modelling of the propagation of underwater sound from the source to the receiver, as well as the interactions with the environment that happen in between. Over the last several decades a wide range of methods have been developed for this purpose, such as normal mode expansion, wavenumber integration, parabolic equation methods, finite difference and element methods, and ray- and beam-tracing³. These methods vary in their implementation, complexity, and range of validity, but they all aim to solve the wave equation,

 

 

where P is the amplitude of pressure perturbations, c is the sound speed at ambient pressure and density, and f is a source term. The simplest and computationally most efficient of these methods is ray-tracing, which models the sound field by approximating the propagation of acoustic wavefronts with a finite number of rays that are integrated to give the position, amplitude, and phase as functions of time. The final sound field is then reconstructed at a given position and time by summing the contribution of nearby rays.

 

The ray tracing equations are obtained as the high-frequency limit of the wave equation, which in two dimensions, range r and elevation z, are³

 

 

where s is the path length. In the following, we assume the sound profile to be range independent, such that ρ is constant. The ray equations (2) and (3) are promoted to a beam tracer by including two additional equations for the curvature and width of the ray wavefronts, given by q and p, that are integrated along the ray. These are in general complex-valued, but in this work we use the simplifying initial conditions q(0) = 0 and p(0) = 1, such that they become real-valued ⁴,

 

 

From these we get the width W and geometric change in ray amplitude A = P/P _init as

 

 

where β is the take-off angle, and δβ is the initial angular spacing of the rays.

 

Only specular reflections are include d in LASSIE, and the amplitude of the reflected wave is given as a factor of the incoming wave,  where θ is both the angle of incident and reflection. The reflection factor R is also a function of position since the reflection amplitude can vary across the seabed and sea surface. For sea surface reflections, an extra factor − 1 is included to flip the sign of the wave.

 

To generate a synthetic reverberation signal, we need to find all eigenrays connecting the source and the receiver via scattering sites at the sea surface and seabed. In LASSIE, eigenrays are found by sorting rays into separate groups that can be considered as belonging to the same wave front². These groups are classified according to their ray history; the sign of the initial take-off angle, the number of reflections at the seabed, the number of reflections at the sea surface, the number of lower turning points, the number of upper turning points, and the type of scattering surface, i.e. sea surface or seabed. This approach is useful, as the number of cells on the numerical grid representing the sea surface and seabed (which function as scattering sites) is generally much larger than the number of rays, so interpolation is necessary for achieving sufficient coverage. In LASSIE, the properties of a given class of rays at grid points between scattering sites are interpolated if the beam widths overlap. At grid points with only one nearby ray, these properties are instead extrapolated, with the amplitude decaying to zero over a length scale on the order of the beam width. After determining the properties of each class of ray across the grid, the final eigenrays connecting the source and the receiver are found by pairing all incoming and outgoing rays at every point. The total geometric change in amplitude is then the product of the paired rays.

 

The amplitude of the scattered wave is given as a factor of the incoming wave,   where a ( x ) is the scattering area, while σ ( x , θ _i , θ _s , ϕ ) is the scattering amplitude and depends on the angle of incident θ _i, the angle of the scattered ray θ _s, and the angle ϕ of the scattered ray out of the plane of incident. As with reflections, the scattering area and the scattering amplitude are also functions of position.

 

Although the above beam tracer is only valid for a 2D slice, it can be used to approximate the reverberation in a 3D scene. This is done by splitting the area around the sonar into a fan of N wedges, each of which can be treated as a 2D slice for which the beam tracer is used to compute the acoustic field. This type of simplification of 3D geometry into separate 2D slices is often referred to as Nx2D, or 2.5D. In the monostatic case, where the set of incoming and outgoing rays at each scattering site is identical, a single 2D simulation per slice suffices. In the bistatic case, on the other hand, two separate Nx2D beam tracing simulations are needed, one centered on the source, and another centered on the receiver. The eigenrays connecting the source and receiver are then obtained by pairing the incoming rays at each scattering site from the source’s Nx2D fan with the outgoing rays at the corresponding (i.e. overlapping) scattering site in the receiver’s Nx2D fan.

 

2.1 Synthetic hydrophone signal

 

The final signal s_{( e,p,s )} near the receiver position from an eigenray e, scattering off point p, along slice s, can be expressed as

 

 

where I = (e, p, s) is a shorthand for the eigenray label, and S(t) is the acoustic signal emitted at the source at time t = 0. The eigenray has a path length s_I , a phase delay τ_I , a total geometric amplitude change A_I , and propagates in the direction of the unit vector dˆ_I towards the receiver. B_{T ,I} is the transmitted beam pattern, γ_I is the Doppler shift, α is the absorption coefficient in dB/m, and ∆x is the position relative to the receiver center. N_S,I and N_R,I are the number of reflections at the sea surface and the total number of reflections, respectively, while R_i,I are the reflection amplitudes of the eigenray. The scattering area and scattering amplitude are denoted by a_I and σ_I . Note that the arguments of the functions and quantities in Eq. (6), such as angles and positions, are suppressed for notational ease. Finally, η_I is a random factor for modelling stochastic scattering amplitudes, and ζ_I  is random phase delay to ensure that echoes from nearby scattering points combine incoherently. In this work we used a Rayleigh distribution with a mean of unity for η_I , and a uniform distribution over the interval [ −2T, 2T ] for ζ_I , where T is the period of the center frequency of the source signal. 

 

The total signal near the receiving position is the sum of echoes from all eigenrays, as well as ambient noise n_s (t, ∆x) = n_s (t + ∆x · d̂_s) from each slice direction d̂_s,

 

 

where the N_s is the number of slices, N_p = N_p (s) is the number of scattering points along slice s, and N_e = N_e (s, p) is the number of eigenrays connecting the source and receiver via scattering point p along slice s. The final synthetic hydrophone data is obtained by sampling s_R (t, ∆x) at the individual hydrophone positions. The ray-tracing simulation used to model the acoustic field is strictly valid only at the receiver center position, and is much more sensitive to changes in depth compared to changes in the horizontal directions. Because of this, the hydrophone positions relative to the receiver center should not be too large in depth.

 

2.2 Example of synthetic data

 

In this section, the theory and methods outlined above are used to generate a sample set of active sonar data that is used in the next section to infer reflection and scattering coefficients from the measured reverberation level using MCMC sampling.

 

 

Figure 1: The left panel shows the elevation map and the boundary between the two zones used for generating the synthetic sonar data. The right panel shows the matched filtered intensity of the beamformed synthetic data, centered on the receiver, which has a course along a bearing of 60^◦.

 

The elevation map used for the ray tracing simulations is shown in the left panel of Figure 1, which was segmented into 0.9^◦ wedges with 1.25m range resolution, centered on both the source and receiver. The source was omni-directional, and placed at the center of the map at a depth of 50m, emitting a 200dB and 0.5s frequency modulated (FM) pulse with center frequency 1000Hz and a bandwidth of 200Hz. The receiver center was positioned at the same depth as the source, and 100m behind it along a 60^◦ bearing course, although with zero speed. The synthetic hydrophones were configured in a 14.7m cardioid towed array with 50 triplets, and a triplet radius of 4cm. The rays were simulated with up to four reflections. 

 

The seabed consisted of two simple zones characterized by constant reflection, R_j (θ_i) = R_j , and scattering functions given by Lambert’s rule, σ_j (θ_i, θ_s) = µ_j sin (θ_i) sin (θ_s), where the index j ∈ {1, 2} denotes the two zones. Reflections at the sea surface were also taken to be constant, R_s (θ_i) = R_s , while surface scattering was neglected. The values used for these coefficients were µ₁ = −20dB, R₁ = −8dB, µ₂ = −15dB, R₁ = −3dB, and R_s = −1dB. A linear sound speed profile was used, decreasing from 1490m/s at the surface to 1480m/s at a depth of 1000m. Finally, white noise with 60dB/Hz was used for the ambient background, and the final synthetic data obtained by sampling Eq. (7) at the hydrophone positions at a rate of 4096Hz for 10s, corresponding to a sonar range of around 7.5km. The matched filtered intensity of the synthetic sonar data after cardioid beamforming is shown in the right panel of Figure 1. 

 

 

3 BAYESIAN INFERENCE WITH MARKOV CHAIN MONTE CARLO SAMPLING

 

MCMC sampling is a versatile tool for estimating probability distributions that are otherwise too difficult to study by analytic means. It is particularly useful in Bayesian inference, where the goal is to find the posterior distribution P(p|d̂) of a set of parameters p = (p₁ , p₂ , . . . , p_M), given a set of observations d̂  _{= (} d̂₁ , d̂₂ , . . . ,d̂_N ). Once the posterior distribution has been found, quantities such as expectation values and credible regions can be computed. MCMC methods estimate the posterior distribution via the likelihood function, L(d̂ |p), and Bayes’ rule¹_, P( p|d̂ ) ∝ L(d̂ |p) × π(p). The likelihood is, in a sense, the reverse of the posterior, i.e. the probability of measuring d̂ given that p are the actual parameters, although not exactly, since the posterior is also determined by the prior π(p), which quantifies the prior information or belief about the unknown parameters. A common choice for the likelihood, which is used in this work, is to assume that the observations d̂  are uncorrelated and drawn from separate Gaussian distributions with measurement error σ̂_i = (σ̂ ₁, σ̂ ₂, . . . , σ̂ _N ), such that the likelihood becomes

 

 

where d(p) = (d₁(p), d₂(p), . . . , d_N (p)) are the theoretically expected values of the measurements given the parameters p. The prior is assumed to be uniform in decibels over an interval much larger than the relevant parameter space. 

 

MCMC methods operate by generating samples from the target distribution, such that for a sufficiently large number of samples, the resulting histogram becomes a reasonable estimate of the distribution. This is done by combining the main properties of its namesakes; the next sample to be generated only depends on the previous sample, thus creating a Markov chain of samples, and it does so in a stochastic way, hence Monte Carlo. In this work we used the Metropolis algorithm¹.

 

3.1 Estimation of reflection and scattering coefficients

 

In this section MCMC is used to infer the scattering and reflection coefficients of the seabed and sea surface given the measured reverberation level of the synthetic data from Section 2.2. The data points used were the decibel levels of the average matched filtered intensity in 0.3s bins in 32 directions along the bearings 45^◦- 135^◦ and 225^◦- 315^◦, relative to the course. The first second of data was dominated by the direct pulse from the sonar, and therefore discarded. Bearings near end-fire were also discarded due to decreased angular resolution, and the expectation that for real towed arrays, these directions might have substantial noise from the platform the array is attached to. 

 

The theoretical values d_i (p) in the likelihood function (8) were computed using the same theory and methods as was used for generating the synthetic hydrophone data, but with two major differences in order to decrease computing time. First, the eigenrays were computed using a much coarser grid, with 30m range resolution instead of 1.25m, along 64 slices instead of 400. Second, the amplitude of the eigenrays in Eq. (6) were used to directly estimate the beamformed and matched filtered intensity in each bin. 

 

Marginal 1D posteriors, along with expectation values and uncertainties, are given in Figure 2, which all agree with the ground truth within a standard deviation. Credible regions for the pairwise 2D posteriors are also given in Figure 2, and shows, among other things, the correlations between the coefficients. Most of the coefficients are inside their respective 68% credible regions, expect for the pairs R₂−R₁ and R₂−R_s , which are just outside of their 95% credible regions. This discrepancy might indicate a modelling error from using a much coarser grid for computing the expected data d_i (p). For instance, around six slices were used per beamforming width to compute the synthetic data, while the fast computation of d_i (p) only used roughly one slice per beamforming width, and therefore might have produced less accurate estimates in directions where the seabed has relatively large variations from one slice to another. 

 

To get an indication of how modeling choices or incorrect environmental data might affect the estimates of the reflection and scattering coefficients, three more MCMC simulations were run; two with a reduced number of reflections during ray tracing, and one with an incorrect constant sound speed profile of 1500m/s. The number of reflections directly impacts the number of eigenrays that contribute to the final signal, and therefore also the computational cost of generating each MCMC sample. Hence, using the smallest number of reflections without significantly affecting the final result is desirable in order to have an efficient sampler. The sound speed profile, on the other hand, is a crucial piece of additional data needed for accurately simulating the sound paths, and must be measured during an experiment at sea, or reconstructed using either historical data or some other approach, both of which are unlikely to be exact.

 

 

Figure 2: The 95% (outer contour) and 68% (inner contour) credible regions of the 2D marginal posteriors, and the 1D posterior, of each seabed and sea surface coefficient. The expectation value and standard deviation of each parameter are given in the title of the 1D subplots, while the ground truth values are indicated with crosses and dotted lines.

 

 

Figure 3: Normalized 1D posteriors of each seabed and sea surface coefficient from different MCMC simulations. The MCMC simulations differ in the ray-tracing parameters used in the fast computation of the matched filter signal; using the same ray-tracing parameters as used for generating the synthetic data (maximum four reflections and a linear sound speed profile), using three reflections, two reflections, and a constant sound speed profile. The ground truth of the synthetic data is indicated with dotted lines.

 

The 1D posteriors of these MCMC simulations are shown in Figure 3. Reducing the number of reflections from four to three does not cause a significant change in the final estimates, whereas a reduction to two reflections causes the algorithm to fail, especially for the reflection coefficients. For eigenrays with an increasing number of reflections, the contribution to the final signal is expected to eventually become negligible compared to more direct paths, or even the ambient noise. At some point, there is simply no more accuracy to be gained from including more reflections, as seems to almost be the case for three reflections in Figure 3. Decreasing the number of reflections, on the other hand, eliminates increasingly important echoes, until the measured signal can no longer be faithfully reconstructed, leading to incorrect results. This type of error can, however, be mitigated through consistency checks, e.g. by finding the minimum number of reflections needed to yield unchanging results. Detecting an inaccurate sound speed profile can be much more challenging problem, if not impossible, although one way of mitigating the effect of this kind of error is to run several MCMC simulations, each with a different plausible sound speed profile, and using a weighted average of these for the final posterior. Doing so is likely to be much more computationally demanding, and result in weaker (albeit possibly more appropriate) constraints on the sampled coefficients.

 

4 CONCLUDING REMARKS

 

In this work we have outlined the theory and methods currently used in a new tool called LASSIE for synthesizing hydrophone data in active sonar, which can be used for, among other things, modelling sonar performance in realistic environments, creating test data for processing pipelines or novel processing techniques, or generating training data for machine learning. There are many avenues for further development, for instance, adding volume reverberation, multiple bistatic sources, a frequency dependent response in the attenuation and backscattering, and additional sources of noise, such as self noise, distant ship noise, and ambient noise from other directions than the horizontal. However, there are some fundamental limitations due to the approximations and simplifications used that compromise the fidelity of the final synthesized signal. Chief among these is the partitioning of the 3D scene into Nx2D wedges, which restricts the simulated rays to 2D slices, thereby losing out-of-plane reflections. Additionally, the model only accounts for specular reflections, since incorporating diffuse reflections would require creating multiple new rays in the diffuse directions at every reflection point, quickly leading to overwhelming computational demands.

 

Since ray and beam tracing is a high frequency approximation of the wave equation, every scattering and reflection event is treated as though there is zero penetration into the seabed. Because of this, subsurface scattering and the re-emission of subsurface waves back into the water volume are not included in the synthesized data, although these effects are not expected to be important in most applications of interest.

 

As an example application of LASSIE, synthetic data from an active triplet array was generated for a simple 3D scene consisting of two distinct seabed zones. The synthetic hydrophone data was then beamformed, matched filtered, and binned, to produce a set of measurements of the reverberation level. This data was input into a MCMC sampler, which reproduced the seabed and sea surface reflection and scattering coefficients within one standard deviation. The robustness of these results were checked by rerunning the MCMC simulations with both fewer reflections, and with an incorrect sound speed profile, illustrating that the estimates can be very sensitive to both poorly chosen modeling parameters and inaccurate input data. However, the former type of error can be mitigated by performing consistency checks on the simulated output, while the latter can be taken into account by marginalizing over the uncertain data. It should also be noted, that the MCMC sampler used the same theory and model for computing the expected data points d, as was used for generating the synthetic hydrophone data and measurements dˆ. Real data is much more complicated and messy, and contains physical effects that the model outlined in this work simply cannot reproduce.

 

Although the simple Metropolis algorithm used in this study worked well for the synthesized test data at hand, employing a more sophisticated algorithm may be beneficial or even necessary for cases involving a larger set of parameters. For example, an adaptive MCMC algorithm or Hamiltonian Monte Carlo could help address the curse of dimensionality and decrease the time required for the sampler to converge to the target distribution¹.

 

REFERENCES

 

1. S. Brooks, A. Gelman, G. Jones, and X.-L. Meng. Handbook of Markov Chain Monte Carlo. CRC Press, 2011.

2. J. H. Hovem. PlaneRay: An acoustic underwater propagation model based on ray tracing and plane-wave reflection coefficients. Technical Report 2008/00610, FFI, 2008.

3. F. B. Jensen, W. A. Kuperman, M. B. Porter, and H. Schmidt. Computational Ocean Acoustics. Springer Publishing Company, Incorporated, 2nd edition, 2011.

4. M. B. Porter. Beam tracing for two- and three-dimensional problems in ocean acoustics. The Journal of the Acoustical Society of America, 146(3):2016–2029, Oct. 2019.