Underwater acoustic propagation modeling with arlpy and BellhopĀ¶

The underwater acoustic propagation modeling toolbox (uwapm) in arlpy is integrated with the popular Bellhop ray tracer distributed as part of the acoustics toolbox. In this notebook, we see how to use arlpy.uwapm to simplify the use of Bellhop for modeling.

PrerequisitesĀ¶

You can install arlpy (in overall system) by typing this in your command prompt:

InĀ [Ā ]:
python3 -m pip install arlpy

Getting startedĀ¶

Start off with checking that everything is working correctly:

InĀ [1]:
import arlpy.uwapm as pm
import arlpy.plot as plt
import numpy as np
InĀ [2]:
pm.models()
Out[2]:
['bellhop']

The bellhop model should be listed in the list of models above, if everything is good. If it isn't listed, it means that bellhop.exe is not available on the PATH, or it cannot be correctly executed. Ensure that bellhop.exe from the acoustics toolbox installation is on your PATH (updated .profile or equivalent, if necessary, to add it in).

From here on we assume that the bellhop model is available, and proceed...

We next create an underwater 2D environment (with default settings) to model:

InĀ [3]:
env = pm.create_env2d()
pm.print_env(env)
                name : arlpy
   bottom_absorption : 0.1
      bottom_density : 1600
    bottom_roughness : 0
   bottom_soundspeed : 1600
               depth : 25
        depth_interp : linear
           frequency : 25000
           max_angle : 80
           min_angle : -80
              nbeams : 0
            rx_depth : 10
            rx_range : 1000
          soundspeed : 1500
   soundspeed_interp : spline
             surface : None
      surface_interp : linear
            tx_depth : 5
   tx_directionality : None
                type : 2D

We can see the default values above. Most numbers are in SI units. The only ones that aren't are:

  • bottom_absoption is in dB/wavelength
  • min_angle and max_angle are in degrees

The default environment is an isovelocity Pekeris waveguide with a water depth of 25 m, a omnidirectional transmitter at 5 m depth, and a receiver at 10 m depth 1 km away. We can visualize the environment (not to scale):

InĀ [5]:
pm.plot_env(env, width=900)

Let's simulate it and see what the eigenrays between the transmitter and receiver look like:

InĀ [6]:
rays = pm.compute_eigenrays(env)
pm.plot_rays(rays, env=env, width=900)

We can also compute the arrival structure at the receiver:

InĀ [7]:
arrivals = pm.compute_arrivals(env)
pm.plot_arrivals(arrivals, width=900)

The arrivals are returned in pandas data frame that we can query, if we like. For example, we can look up the time of arrival, angle of arrival, and the number of surface/bottom bounces for the first 10 arrivals:

InĀ [8]:
arrivals[arrivals.arrival_number < 10][['time_of_arrival', 'angle_of_arrival', 'surface_bounces', 'bottom_bounces']]
Out[8]:
time_of_arrival angle_of_arrival surface_bounces bottom_bounces
1 0.721796 22.538258 9 8
2 0.716791 -21.553932 8 8
3 0.709687 20.052084 8 7
4 0.705227 -19.034420 7 7
5 0.698960 17.484425 7 6
6 0.695070 -16.436064 6 6
7 0.689678 14.842228 6 5
8 0.686383 -13.766300 5 5
9 0.681901 12.133879 5 4
10 0.679223 -11.034208 4 4

We can also convert to a impulse response time series, if we want to use it for further signal processing:

InĀ [9]:
ir = pm.arrivals_to_impulse_response(arrivals, fs=96000)
plt.plot(np.abs(ir), fs=96000, width=900)

More complex environmentsĀ¶

So far, we have a simple isovelicty Pekeris waveguide. Let us next have something more interesting - an environment with some bathymetric structure and a more complicated soundspeed profile.

Bathymetry and soundspeed profileĀ¶

Let's first start off by defining our bathymetry, a steep up-slope for the first 300 m, and then a gentle downslope:

InĀ [10]:
bathy = [
    [0, 30],    # 30 m water depth at the transmitter
    [300, 20],  # 20 m water depth 300 m away
    [1000, 25]  # 25 m water depth at 1 km
]

and then our soundspeed profile:

InĀ [11]:
ssp = [
    [ 0, 1540],  # 1540 m/s at the surface
    [10, 1530],  # 1530 m/s at 10 m depth
    [20, 1532],  # 1532 m/s at 20 m depth
    [25, 1533],  # 1533 m/s at 25 m depth
    [30, 1535]   # 1535 m/s at the seabed
]

Now we can create our environment with a muddy bottom, and a transmitter that is at 20 m depth:

InĀ [12]:
env = pm.create_env2d(
    depth=bathy,
    soundspeed=ssp,
    bottom_soundspeed=1450,
    bottom_density=1200,
    bottom_absorption=1.0,
    tx_depth=15
)
InĀ [13]:
pm.print_env(env)
                name : arlpy
   bottom_absorption : 1.0
      bottom_density : 1200
    bottom_roughness : 0
   bottom_soundspeed : 1450
               depth : [[   0.   30.]
                        [ 300.   20.]
                        [1000.   25.]]
        depth_interp : linear
           frequency : 25000
           max_angle : 80
           min_angle : -80
              nbeams : 0
            rx_depth : 10
            rx_range : 1000
          soundspeed : [[   0. 1540.]
                        [  10. 1530.]
                        [  20. 1532.]
                        [  25. 1533.]
                        [  30. 1535.]]
   soundspeed_interp : spline
             surface : None
      surface_interp : linear
            tx_depth : 15
   tx_directionality : None
                type : 2D
InĀ [14]:
pm.plot_env(env, width=900)

We can also plot the soundspeed profile used:

InĀ [15]:
pm.plot_ssp(env)

Looks more interesting! Let's see what the eigenrays look like, and also the arrival structure:

InĀ [16]:
rays = pm.compute_eigenrays(env)
pm.plot_rays(rays, env=env, width=900)

We could also ignore the receiver, and plot rays launched at various angles:

InĀ [17]:
rays = pm.compute_rays(env)
pm.plot_rays(rays, env=env, width=900)

or place lots of receivers in a grid to visualize the acoustic pressure field (or equivalently transmission loss). We can modify the environment (env) without having to recreate it, as it is simply a Python dictionary object:

InĀ [18]:
env['rx_range'] = np.linspace(0, 1000, 1001)
env['rx_depth'] = np.linspace(0, 30, 301)
InĀ [19]:
tloss = pm.compute_transmission_loss(env)
pm.plot_transmission_loss(tloss, env=env, clim=[-60,-30], width=900)

We see a complicated interference pattern, but an interesting focusing at 800 m at a 15 m depth. The detailed interference pattern is of course sensitive to small changes in the environment. A less sensitive, but more averaged out, transmission loss estimate can be obtained using the incoherent mode:

InĀ [20]:
tloss = pm.compute_transmission_loss(env, mode='incoherent')
pm.plot_transmission_loss(tloss, env=env, clim=[-60,-30], width=900)

Source directionalityĀ¶

Now, let's use a directional transmitter instead of an omni-directional one:

InĀ [21]:
beampattern = np.array([
    [-180,  10], [-170, -10], [-160,   0], [-150, -20], [-140, -10], [-130, -30],
    [-120, -20], [-110, -40], [-100, -30], [-90 , -50], [-80 , -30], [-70 , -40],
    [-60 , -20], [-50 , -30], [-40 , -10], [-30 , -20], [-20 ,   0], [-10 , -10],
    [  0 ,  10], [ 10 , -10], [ 20 ,   0], [ 30 , -20], [ 40 , -10], [ 50 , -30],
    [ 60 , -20], [ 70 , -40], [ 80 , -30], [ 90 , -50], [100 , -30], [110 , -40],
    [120 , -20], [130 , -30], [140 , -10], [150 , -20], [160 ,   0], [170 , -10],
    [180 ,  10]
])
env['tx_directionality'] = beampattern
InĀ [22]:
tloss = pm.compute_transmission_loss(env)
pm.plot_transmission_loss(tloss, env=env, clim=[-60,-30], width=900)

Now you can see the directionality and the sidelobe structure of the transmitter.

Undulating water surfaceĀ¶

Finally, let's try adding a long wavelength swell on the water surface:

InĀ [23]:
surface = np.array([[r, 0.5+0.5*np.sin(2*np.pi*0.005*r)] for r in np.linspace(0,1000,1001)])
env['surface'] = surface
InĀ [24]:
tloss = pm.compute_transmission_loss(env)
pm.plot_transmission_loss(tloss, env=env, clim=[-60,-30], width=900)

Now, if I placed a receiver at 800 m, and 15 m depth, roughly where we see some focusing, what would the eigenrays and arrival structure look like?

InĀ [25]:
env['rx_range'] = 800
env['rx_depth'] = 15
InĀ [26]:
rays = pm.compute_eigenrays(env)
pm.plot_rays(rays, env=env, width=900)
InĀ [27]:
arrivals = pm.compute_arrivals(env)
pm.plot_arrivals(arrivals, dB=True, width=900)

We plotted the amplitudes in dB, as the later arrivals are much weaker than the first one, and better visualized in a logarithmic scale.