Propagation & channel modeling

Quickstart guide

Propagation modeling

Let’s get started:

using UnderwaterAcoustics
using Plots

Creating an environmental description

We typically start with an environmental description:

env = UnderwaterEnvironment()

UnderwaterEnvironment(
  bathymetry = 100.0, 
  altimetry = 0.0, 
  temperature = 27.0, 
  salinity = 35.0, 
  soundspeed = 1538.9235842, 
  density = 1022.7198310217424, 
  seabed = RigidBoundary, 
  surface = PressureReleaseBoundary, 
)

If the defaults don’t suit our needs, we can customize the environment. For example, if we wanted an environment with 20 m water depth, sandy-clay seabed, and a smooth sound speed profile with 1500 m/s near the surface, 1490 m/s at 10 m depth, and 1520 m/s near the seabed, we could define:

env = UnderwaterEnvironment(
  bathymetry = 20.0,
  seabed = SandyClay,
  soundspeed = SampledField([1500, 1490, 1520]; z=0:-10:-20, interp=:cubic)
)

UnderwaterEnvironment(
  bathymetry = 20.0, 
  altimetry = 0.0, 
  temperature = 27.0, 
  salinity = 35.0, 
  soundspeed = SampledField(z-varying, 3 samples), 
  density = 1022.7198310217424, 
  seabed = FluidBoundary(ρ=1173.381, c=1504.9272, δ=0.00242), 
  surface = PressureReleaseBoundary, 
)

If we have Plots.jl installed, we can use plot recipes to plot the environment or the sound speed profile. For example:

plot(env.soundspeed)

Let us construct a range-independent iso-velocity environment with water depth of 20 m:

env = UnderwaterEnvironment(bathymetry=20, seabed=SandyClay)

UnderwaterEnvironment(
  bathymetry = 20, 
  altimetry = 0.0, 
  temperature = 27.0, 
  salinity = 35.0, 
  soundspeed = 1538.9235842, 
  density = 1022.7198310217424, 
  seabed = FluidBoundary(ρ=1173.381, c=1504.9272, δ=0.00242), 
  surface = PressureReleaseBoundary, 
)

Selecting a propagation model

Once we have an environment, we need to select a propagation model. Since we have a range-independent iso-velocity environment, we can use the PekerisRayTracer:

pm = PekerisRayTracer(env)

PekerisRayTracer(h=20, max_bounces=3)

Tip

We can get a list of all available propagation models by calling models(), and a shortlist of all models compatible with a given environment by calling models(env).

Setting up transmitters and receivers

Next, we need a source and a receiver:

tx = AcousticSource(0.0, -5.0, 1000.0)

TX[1000.0 Hz, 0 dB](0.0, 0.0, -5.0)

rx = AcousticReceiver(100.0, -10.0)

RX(100.0, 0.0, -10.0)

For improved readability, positions can be specified as tuples or named tuples, and the API fully supports Unitful.jl:

tx = AcousticSource((x=0u"m", z=-5u"m"), 1u"kHz")

TX[1000 Hz, 0 dB](0.0, 0.0, -5.0)

rx = AcousticReceiver(100u"m", -10u"m")

RX(100.0, 0.0, -10.0)

Tip

2-tuples are interpreted as (x, z) and 3-tuples as (x, y, z). The coordinate system has x and y axis in the horizontal plane, and z axis pointing upwards, with the nominal water surface being at 0 m. This means that all z coordinates in water are negative. If units are not specified, they are assumed to be S.I. units (meters for distances, Hz for frequency, etc).

We just defined an omnidirectional 1 kHz transmitter tx at a depth of 5 m at the origin, and an omnidirectional receiver rx at a range of 100 m and a depth of 10 m.

Running the model

Now that we have an environment, a propagation model, a transmitter and a receiver, we can modeling. First, we ask for all ray arrivals (eigenrays) between the transmitter and receiver:

rays = arrivals(pm, tx, rx)

7-element Vector{UnderwaterAcoustics.RayArrival{Float64, Float64, Float64, Float64, Vector{@NamedTuple{x::Float64, y::Float64, z::Float64}}}}:
 ∠ -2.9°  0⤒  0⤓ ∠  2.9° |  65.06 ms | -40.0 dB ϕ   0.0°↝
 ∠  8.5°  1⤒  0⤓ ∠  8.5° |  65.71 ms | -40.1 dB ϕ 180.0°↝
 ∠-14.0°  0⤒  1⤓ ∠-14.0° |  66.98 ms | -62.6 dB ϕ 170.9°↝
 ∠ 19.3°  1⤒  1⤓ ∠-19.3° |  68.85 ms | -74.3 dB ϕ -23.4°↝
 ∠-24.2°  1⤒  1⤓ ∠ 24.2° |  71.26 ms | -80.4 dB ϕ-145.6°↝
 ∠ 28.8°  2⤒  1⤓ ∠ 28.8° |  74.16 ms | -73.5 dB ϕ  10.8°↝
 ∠-33.0°  1⤒  2⤓ ∠-33.0° |  77.50 ms | -100.7 dB ϕ-167.2°↝

For each eigenray, this shows us the launch angle, number of surface bounces, number of bottom bounces, arrival angle, travel time, transmission loss along that ray, and phase change. The last “↝” symbol indicates that the complete ray path is also available. We can plot the ray paths:

plot(env; xlims=(-10,110))
plot!(tx)
plot!(rx)
plot!(rays)

The red star is the transmitter and the blue circle is the receiver. The stronger eigenrays are shown in blue, while the weaker ones are shown in red.

Often, we are interested in the arrival structure at a receiver. We generate an impulse response sampled at 48 kSa/s and plot it:

plot(impulse_response(pm, tx, rx, 48000))

The first arrival is at 0 ms, as the impulse response is generated with timing relative to the first arrival. If we wanted timings relative to the transmission time, we can set keyword argument abstime to true:

plot(impulse_response(pm, tx, rx, 48000; abstime=true))

We can also get the complex acoustic field or the transmission loss in dB:

acoustic_field(pm, tx, rx)

0.011324625487986332 + 0.013441467920062827im

transmission_loss(pm, tx, rx)

35.101536894459066

We can also pass in arrays of receivers, if we want to compute transmission loss at many locations simultaneously. Some models are able to compute transmission loss on a Cartesian grid very efficiently. This is useful to plot transmission loss as a function of space.

To define a 1000×200 Cartesian grid with 0.1 m spacing and compute the transmission loss over the grid:

rxs = AcousticReceiverGrid2D(1.0:0.1:100, -20:0.1:0)
x = transmission_loss(pm, tx, rxs)

991×201 Matrix{Float64}:
 19.1283  19.2328  19.5428  19.8492  …   8.82646  11.1099  16.4209  87.1227
 19.1329  19.2373  19.5471  19.8536      8.87496  11.1689  16.4855  87.123
 19.1381  19.2424  19.552   19.8586      8.92702  11.2323  16.5549  87.1233
 19.1442  19.2483  19.5577  19.8643      8.9824   11.2997  16.6288  87.1237
 19.151   19.255   19.5641  19.8709      9.04085  11.3708  16.7068  87.124
 19.1587  19.2626  19.5715  19.8784  …   9.10213  11.4454  16.7886  87.1244
 19.1674  19.2712  19.5799  19.8869      9.16604  11.5232  16.8739  87.1248
 19.1772  19.2809  19.5893  19.8965      9.23241  11.604   16.9625  87.1253
 19.1882  19.2918  19.6     19.9073      9.30109  11.6875  17.0541  87.1258
 19.2005  19.3039  19.6119  19.9194      9.37202  11.7737  17.1485  87.1263
  ⋮                                  ⋱                               ⋮
 37.3745  37.4387  37.4769  37.4921     57.8844   61.3578  67.2495  98.1624
 37.4079  37.4741  37.5144  37.5317     57.9056   61.3788  67.2688  98.1837
 37.4422  37.5104  37.5527  37.5722     57.9268   61.4001  67.2893  98.205
 37.4774  37.5476  37.592   37.6136  …  57.9481   61.4217  67.3109  98.2263
 37.5134  37.5857  37.6321  37.6558     57.9693   61.4435  67.3336  98.2477
 37.5504  37.6247  37.6732  37.6989     57.9903   61.4655  67.3573  98.269
 37.5882  37.6646  37.7151  37.7428     58.0112   61.4877  67.3819  98.2904
 37.6269  37.7053  37.7578  37.7876     58.0319   61.5099  67.4075  98.3118
 37.6665  37.7469  37.8014  37.8332  …  58.0524   61.5322  67.4338  98.3333

plot(env; xlims=(0,100))
plot!(rxs, x)

Channel modeling

Simulating a transmission

Once we have a propagation model, we can setup a channel simulation. Let’s say we had a 10 kHz source transmitting at a source level of 170 dB re 1 µPa @ 1m and we wanted to run our simulations at 192 kSa/s:

tx = AcousticSource((x=0u"m", z=-5u"m"), 10u"kHz"; spl=170)
ch = channel(pm, tx, rx, 192000)

SampledPassbandChannel(1×1, 192000.0 Hz)

This channel accepts one input acoustic channel and yields one output acoustic channel. To transmit a signal through the channel, we create a windowed pulse using SignalAnalysis.jl and transmit() it through the channel ch:

using SignalAnalysis

x = cw(10000, 0.001, 192000; window=(tukey, 0.5)) |> real
y = transmit(ch, x)

plot(
  plot(x; xlims=(0,5)),
  plot(y; xlims=(0,5));
  layout=(2,1)
)

The received signal is scaled to be in µPa, assuming the source level of the transmitter was specified in dB re 1 µPa @ 1 m. We can clearly see the multipath arrivals in the received signal. Had we specified keyword argument abstime = true, we would have also seen the signal delayed by 65 ms.

Adding noise

We can also simulate channels with noise. For example, if we wanted red Gaussian noise with standard deviation σ = 0.5 Pa, we can specify that when generating the channel:

ch = channel(pm, tx, rx, 192000; noise=RedGaussianNoise(0.5e6))
y = transmit(ch, x)

plot(
  plot(x; xlims=(0,5)),
  plot(y; xlims=(0,5));
  layout=(2,1)
)

We can see the noise in the timeseries. We can see the 1/f² variation (red Gaussian noise) in the power spectral density superimposed on the 10 kHz peak from the transmitted pulse:

psd(y)