Comparison of two different methods to include a beam pattern in parabolic equation models

 

Alexandra Schäfke and Jan Ehrlich

 

Citation: Proc. Mtgs. Acoust. 47, 070006 (2022); doi: 10.1121/2.0001593

 

View online: https://doi.org/10.1121/2.0001593

 

View Table of Contents: https://asa.scitation.org/toc/pma/47/1

 

Published by the Acoustical Society of America

 

ARTICLES YOU MAY BE INTERESTED IN

 

Accuracy of numerically predicted underwater sound of a ship-like structure

 

Proceedings of Meetings on Acoustics 47, 070001 (2022); https://doi.org/10.1121/2.0001565

 

Correlations between sound speed and density in seabed sediment cores collected in Norwegian waters

 

Proceedings of Meetings on Acoustics 47, 070004 (2022); https://doi.org/10.1121/2.0001588

 

Analysis of hydroacoustic time series by parametric predictive modelling

 

Proceedings of Meetings on Acoustics 47, 055001 (2022); https://doi.org/10.1121/2.0001596

 

Exploring the use of AI in marine acoustic sensor management

 

Proceedings of Meetings on Acoustics 47, 070008 (2022); https://doi.org/10.1121/2.0001601

 

Experimental investigation of a virtual planar array for MIMO sonar systems

 

Proceedings of Meetings on Acoustics 47, 070005 (2022); https://doi.org/10.1121/2.0001591

 

Estimation of uncertainties in underwater sound measurements of ships

 

Proceedings of Meetings on Acoustics 47, 070002 (2022); https://doi.org/10.1121/2.0001571

 

 

 

Comparison of two different methods to include a beam pattern in parabolic equation models

 

Alexandra Schafke and Jan Ehrlich

 

Bundeswehr Technical Centre for Ships and Naval Weapons, Maritime Technology and Research (WTD 71), Eckernforde, Schleswig-Holstein, 24340, GERMANY; alexandraschaefke@bundeswehr.org; alexandra.schaefke@gmail.com; janehrlich@bundeswehr.org

 

While it is straightforward to include the beam pattern of a sonar in ray models and normal mode models, things become more challenging for parabolic equation models. Two methods are compared in this paper. The first method uses an array starter for the parabolic equation model, derived by MacGillivray and Chapman as a generalization of Collins’ self-starter for a point source. Range propagation of the depth-dependent pressure field generated by this array starter then contains the effects of the beam pattern of the respective array. The second method inverts the computation direction so that the beam pattern can be applied to the simulation result for the complex sound pressure. This is done by performing a windowed Fourier transform into vertical wavenumber space at the sonar depth, thereby decomposing the sound pressure field into propagation angles allowing for multiplication with the beam pattern. This method is not limited to parabolic equation models. Both methods have been realized for WTD 71’s parabolic equation code PESSim (Parabolic Equation Sound Simulation). Using a homogeneous half-space with an analytic reference solution, the two methods are validated and compared. Both methods perform well while each has its advantages and disadvantages.

 

1. INTRODUCTION

 

The inclusion of a beam pattern of a source or receiver into a propagation model becomes important when contributions to or from steep propagation angles cannot be neglected. While it is relatively straightforward to implement a beam pattern into a ray or normal mode model, it is more challenging when dealing with a parabolic equation (PE) model as the contributions for different propagation angles are not directly accessible. In this paper, two methods to include a beam pattern in a PE model are examined. The first method is the array starter by MacGillivray and Chapman, 1 the second method consists of a decomposition into propagation angles inspired by similar approaches in the context of modeling scattering and reverberation from surfaces within PE models. 2,3 The performance of the two methods is compared using a test case with an analytic solution.

 

2. METHOD 1: ARRAY STARTER

 

A PE model propagates an initial pressure field over depth p(r₀,z) at a range r 0 in discrete steps to larger ranges. There are different methods to generate the initial pressure field of which Collins’ PE self-starter^4,5 for a point source has become somewhat of a standard. In principle, for an array of point sources, a separate PE simulation can be run for each of them. Combining the individual pressure fields of the point sources then yields the wanted result for an array. MacGillivray and Chapman introduced instead an array starter,¹ which necessitates only a single simulation run. Section 2.1 covers briefly the theory behind the self-starter for a point source,^4,5 before Section 2.2 describes its generalization to an array source.¹

 

A. SELF-STARTER FOR A POINT SOURCE

 

The self-starter for a point source was developed by Collins.^4,5 For a point source at depth z_s in a cylinder-symmetrical, range-independent scenario, the pressure field can be approximated by the normal mode solution⁶

 

 

The normal modes Ψ_m (z) fulfill the equation

 

 

The above equation has already been written with the operator X used in PE models. The reference wave number k₀ originates from the PE ansatz which splits off a cylindrical wave with wave number k₀ from the pressure field,

 

 

For a general analytical function f , the normal mode equation (2) can be generalized to

 

 

Together with the completeness condition Σ_m Ψ^∗_m (z_s)Ψ_m (z) = δ (z − z_s) , this yields for the normal mode solution (1) a formal solution for the field ψ from the PE ansatz (3) using pseudo-differential operators,

 

 

This expression can be solved numerically with the same techniques that are employed for the PE propagator to provide a starting field ψ ( r₀,z) at a range r₀ close to the point source.

 

B. GENERALIZATION FOR AN ARRAY SOURCE

 

In the context of simulating an airgun array, MacGillivray and Chapman generalized the self-starter to an array starter. 1 The pressure field of an array of N point sources with locations ( ⃗r_n , z_n) , n = 1 . . . N (using cylindrical coordinates ( ⃗r, z) = ( r, ø, z ) with a two-dimensional range vector ⃗r instead of ( x, y, z )) can be written as a superposition of the pressure fields of the point sources,

 

 

Due to the symmetry of the individual point sources, p depends only on the norm of the horizontal radial vector. The complex amplitude A_n takes into account different source strengths and phases of the individual point sources. The simulation runs within a computational (r, z) plane defined by its azimuth angle ø . Far away from the array, r ≫| ⃗r_n| , components of  ⃗r_n perpendicular to the computational plane can be neglected,

 

 

where r_n is the projection of  ⃗r_n into the computational plane,

 

 

Using the normal mode solution (1) for the point sources and  the normal mode solution for the pressure field of the array in the computational domain becomes

 

 

A s for the self-starter for a point source, the use of (4) for the reduced field ψ a (r, z) defined by p a (r, z) ≈  yields the starting field

 

 

for the array with the point source starting fields ψ from (5).

 

3. METHOD 2: ANGLE DECOMPOSITION

 

The basic idea of the second method is the decomposition of the complex pressure field into contributions from different propagation angles. This method was inspired by similar approaches within the context of scattering and reverberation in PE models.^2,3 It is not restricted to PE models but can be applied as a post-processing step to any coherent propagation model.

 

A. TRANSFORMATION TO VERTICAL WAVE NUMBER SPACE

 

This method requires the reciprocal computation of the sound pressure field starting from a potential target towards the transmitting or receiving sonar whose beam pattern is to be included. Let p(z) be the sound pressure field over depth at the sonar range r_s .

 

 

Figure 1: Window function w (z) and its transform ŵ(κ) for a rectangular window with width d (blue) and for the smoothed version (12) (red) using f = 500 Hz, z s = 50 m for d = 5λ (left) and d = 10λ (right). The dashed lines mark the region |κ| ≤ k where the vertical wave number κ corresponds to a propagation angle θ. For the smaller window width, the window transform has a wider main lobe. The smooth window function suppresses the side lobes.

 

Within a depth window around the sonar depth, located inside the water column and described by the window function w(z), this pressure field can be transformed from position space with depth coordinate z into vertical wave number space with the vertical wave number κ,

 

 

The window function w(z) yields a convolution of the windowless Fourier transform (κ) with the Fourier transform ŵ(z) of the window function. Examples for window functions and their transforms are shown in Fig. 1. A smooth window function like e.g.

 

 

instead of a rectangular one suppresses side lobes. The parameter d determines the window width. For |κ| ≤ k , where k = ω/c(z_s) is the wave number at sonar depth z_s , a propagation angle θ , relative to horizontal, can be attributed to the vertical wave number via the relation

 

 

The narrower the window function w(z) is around the sonar depth, the wider spread is ŵ(κ) in κ-space. The wider the window transform is spread in κ-space, the more values from the full transform  around a certain κ contribute to the windowed transform  _w (κ) . The result is a certain fuzziness in κ-space.

 

B. BEAM PATTERN INCLUSION FOR PRESSURE FIELD

 

The back transform of (11) is given at sonar depth z_s , where w(z_s) = 1 , by

 

 

To include a beam pattern given by a function B(θ) = |b(θ)|², this can be modified to

 

 

The minus in the beam pattern function results from the reciprocal computation direction. Eq. 15 contains a restriction to |κ| ≤ k which is in general admissible as |κ| > k correspond to non-propagating evanescent modes. But if the window w(z) has been chosen too narrow, the resulting fuzziness in κ-space leads to parts of propagating modes being cut-off also. This leads to an unwanted loss of sound intensity. Therefore, if using the back transform, it is best to define the window function w(z) as wide as possible within the water column.

 

C. BEAM PATTERN INCLUSION FOR INTENSITY

 

If the phase-exact beam pattern function b(θ) is not known but only the intensity beam pattern B(θ) = |b(θ)|², a different approach can be used. For the pressure squares, Parseval’s theorem states

 

 

The right hand side is, up to a normalization, an average of the squared pressure over the window,

 

 

By restricting the integration to |κ| ≤ k and changing the integration variable (dκ = k cosθ dθ), the left side can be rewritten as

 

 

This allows for an inclusion of the intensity beam pattern function B(θ). Inserting (17) and (18) with B(θ) into (16) yields

 

 

Again, the minus in the beam pattern function results from the reciprocal computation direction. Eq. 19 is an approximation of the resulting intensity with beam pattern, containing an inherent averaging over the window width.

 

In this case, the window width should be chosen as narrow as possible while still making sure that the width of the main lobe of the window transform ŵ(κ) is well smaller than 2k. In fact, it should even be smaller than the width of the main lobe of the beam pattern B(θ), as otherwise sound intensity that actually belongs within the angular range of the main lobe would be cut off leading to an underestimation of the sound pressure level. This effect is illustrated in Fig. 2.

 

D. NUMERICAL REALIZATION

 

When implementing the transform via an FFT, the pressure field has to be known with a sufficient depth resolution ∆z. For the resolution in κ-space follows ∆κ = 2π/ (∆zN_FFT ) and for the maximal vertical wave number κ_max ≈ π/∆z. The finer the depth resulution, the more vertical wave numbers and therefore propagation angles are included.

 

 

Figure 2: Effect of choosing the window width too small, illustrated by an example scenario (shallow water, summer sound speed profile, 1 kHz): (a) Comparison of receiver beam pattern B(θ) and normalized squared window transform for several window widths. (b) Relative error of the resulting signal level compared to the array starter solution as a reference. The resulting signal level is being underestimated when the window width is of the same order or larger than the width of the main lobe of the beam pattern.

 

To include all propagation angles, it has to be κ_max ≥ k or ∆z ≤ π/k = λ/2, a reasonable choice is ∆z ≈ λ/8. A sufficiently large N_FFT then ensures a good resolution for the vertical wave number and therefore the propagation angle. As a minimum requirement, N_FFT∆z has to cover the full width where the window function w(z) is non-zero (i.e. 2d for the window function in (12)). Although a smaller window width leads to a lower computational demand, the risks of choosing it too small have been stated above. Window sizes of several wavelengths are recommended.

 

4. COMPARISON OF METHODS

 

To validate and compare the two methods, a scenario with an analytic solution is considered. For an unsteered array of N identical point sources in an iso-halfspace, the normalized sound pressure is given by

 

 

where ( x_n, y_n, z_n ) , n = 1 . . . N are the coordinates of the point sources and the azimuth angle ø determines the computational (r, z) plane. For the following example, a linear array in the (y, z) plane with coordinates (angle α between array and y-axis)

 

 

and parameters f = 500 Hz, c = 1500 m/s, z_s = 50 m, N = 10 , δ = λ/2 = 1 . 5 m and α = 60^◦ has been chosen. The resulting beam pattern in the far field is

 

 

Figure 3a shows this beam pattern for ø = 20^◦. The main lobe is pointing slightly upwards and has its maximum for θ ≈−11.17^◦ .

 

The array starter method from Section 2 was implemented in WTD 71’s PE model PESSim and configured with the coordinates given in (21). The result for propagation loss is shown in Fig. 4 (upper right) in comparison to the exact solution (upper left). The mean relative error over range and depth ⟨|(PL−PL_ref)/PL_ref|⟩_r,z is about 0.0044.

 

 

Figure 3: (a) Far-field beam pattern (22) of the linear array used in the test case. The main lobe is centered around θ ≈− 11.17^◦ , pointing slightly upwards. (b) Propagation loss at z = 80 m over range for the test case: analytic solution (blue), PE solution (PESSim) with array starter (red, dashed) and reciprocal PE solution (PESSim) for a point source combined with angle decomposition (15) using phase exact beam pattern (yellow, dotted) show excellent agreement, angle decomposition (19) using intensity beam pattern (purple) contains the expected inherent depth averaging.

 

 

Figure 4: Propagation loss over range and depth for the test case: analytic solution (upper left), PE calculation (PESSim) with array starter (upper right), reciprocal analytic solution for a point source combined with angle decomposition (15) using phase exact beam pattern (lower left) and angle decomposition (19) using intensity beam pattern (lower right). The agreement is excellent, for the method of angle decomposition using intensity beam pattern the expected inherent depth averaging leads to a loss in detail.

 

For the angle decomposition method, the result of the reciprocal calculation from a target at depth z to the sonar situated at depth z_s at a range r from the target is needed. This can be obtained either by separate PE runs for each target depth (with a point source self-starter) or via the analytic solution

 

 

With the window function (12) and the beam pattern b(θ) from (22) the pressure field p_b(r, z_s) results from (15). Alternatively, using B(θ) the squared pressure results from (19). Using ∆z = λ/8 , d = 10λ and N_FFT = 1024, the results for the propagation loss are shown in Fig. 4 (lower left) and (lower right). The mean relative errors are 0.0027 and 0.0377. In the squared pressure version of the angle decomposition, the inherent averaging over the depth window is visible. Figure 3b compares the results for propagation loss at a depth of 80 m.

 

5. CONCLUSION

 

Two methods to include a beam pattern into a PE model have been compared, an array starter to be used as initial starting field for the PE propagation (Sec. 2) and a decomposition into propagation angles used as a post-processing step on a pressure field obtained by a reciprocal simulation run for a point source (Sec. 3).

 

While both methods show an excellent performance for the test case (Sec. 4), they have their individual advantages and disadvantages depending on the scenario. For the angle decomposition, the required reciprocal computation towards the sonar necessitates individual simulation runs for different target depths and, in the case of range-dependent scenarios, even for different target ranges. The array starter works with an array of point sources, for which positions, amplitudes and phases have to be known while the angle decomposition requires a beam pattern (either phase exact or for intensity) which can easily be calculated from array data but not necessarily vice versa. If only an intensity beam pattern is provided, the angle decomposition method allows still for a solution albeit one containing some inherent averaging.

 

REFERENCES

 

1. A. O. MacGillivray, N. R. Chapman, “Modeling underwater sound propagation from an airgun array using the parabolic equation method”, Canadian Acoustics 40 (1), 19-25 (2012).
2. M. E. Moore-Head, W. Jobst, E. Holmes, “Parabolic-equation modeling with angle-dependent surface loss”, J. Acoust. Soc. Am. 86 (1), 247-251 (1989).
3. H. G. Schneider, “Surface loss, scattering and reverberation with the split-step parabolic equation model”, J. Acoust. Soc. Am. 93 (2), 770-775 (1993).
4. M. D. Collins, “A self-starter for the parabolic equation method”, J. Acoust. Soc. Am. 92 (4), 2069-2074 (1992).
5. M. D. Collins, “The stabilized self-starter for the parabolic equation method”, J. Acoust. Soc. Am. 106 (4), 1724-1726 (1999).
6. F. B. Jensen, W. A. Kuperman, M. B. Porter, H. Schmidt, Computational Ocean Acoustics 2nd edition, Springer, 2011