Geoacoustic inversion with an acoustic profiler

References

This tutorial is adapted from Example B presented in:

A version of this example was also presented in the UComms 2020 webinar:

Problem statement

We consider a geoacoustic inversion problem where we have a static omnidirectional broadband acoustic source transmitting in a 5–7 kHz band, and a single omnidirectional receiver that records the signal at a fixed range. The receiver is able to profile the water column, and therefore makes transmission loss measurements at various depths. We wish to estimate seabed parameters from the transmission loss measurements. Do note that although we have acoustic measurements at various depths, they cannot be used for beamforming to separate out the bottom reflected arrival from other arrivals. We therefore only have transmission loss at each depth for our inversion.

Dataset

To illustrate this idea, let us generate a synthetic dataset for a known set of seabed parameters (density ρ = 1500 kg/m³, relative soundspeed c = 1850 m/s, and attenuation δ = 0.001). The environment is assumed to be an iso-velocity and with a constant depth of 20 m. The source is at a depth of 5 m. The receiver is at a range of 100 m from the source, and makes measurements at depths from 1 to 19 m in steps of 1 m.

Since we have an range-independent iso-velocity environment, we can use the PekerisRayTracer (otherwise we could use the RaySolver):

using UnderwaterAcoustics
using DataFrames

function 𝒴(θ)
  r, d, f, ρ, c, δ = θ
  env = UnderwaterEnvironment(
    bathymetry = 20.0,
    seabed = FluidBoundary(ρ, c, δ)
  )
  tx = AcousticSource(0.0, -5.0, f)
  rx = AcousticReceiver(r, -d)
  pm = PekerisRayTracer(env)
  transmission_loss(pm, tx, rx)
end

data = DataFrame([
  (depth=d, frequency=f, xloss=𝒴([100.0, d, f, 1500.0, 1850.0, 0.001]))
  for d ∈ 1.0:1.0:19.0 for f ∈ 5000.0:100.0:7000.0
])

Probabilistic model

We use some very loose uniform priors for ρ, c and δ, and estimate the transmission loss using the same model 𝒴, as used in the data generation, but without information on the actual seabed parameters. We assume that the measurements of transmission loss are normally distributed around the modeled transmission loss, with a covariance of 0.5 dB.

We define the probabilistic model as a Turing.jl model:

using Turing

# depths d, frequencies f, transmission loss measurements x
@model function geoacoustic(d, f, x)
  ρ ~ Uniform(1000.0, 2000.0)
  c ~ Uniform(750.0, 2000.0)
  δ ~ Uniform(0.0, 0.003)
  μ = [𝒴([100.0, d[i], f[i], ρ, c, δ]) for i ∈ 1:length(d)]
  x ~ MvNormal(μ, 0.5)
end

Variational inference

Once we have the model defined, we can run Bayesian inference on it. We could either use MCMC methods from Turing, or variational inference. Since our model is differentiable, we choose to use the automatic differentiation variational inference (ADVI):

# this may take a minute or two to run...
q = vi(
  geoacoustic(data.depth, data.frequency, data.xloss),
  ADVI(100, 1000)
)

The returned q is a 3-dimensional posterior probability distribution over the parameters ρ, c and δ. We can estimate the mean of the distribution by drawing random variates and taking the sample mean:

μ = mean(rand(q, 10000); dims=2)
3×1 Matrix{Float64}:
 1500.4884323728988
 1850.0874858495117
    0.001025261073074483

We see that the estimated parameter means for ρ, c and δ are quite close to the actual values used in generating the data.

We can also plot the conditional distributions of each parameter:

using StatsPlots

plot(ρ -> pdf(q, [ρ, μ[2], μ[3]]), 1300, 1700; xlabel="ρ", legend=false)
plot(c -> pdf(q, [μ[1], c, μ[3]]), 1650, 1950; xlabel="c", legend=false)
plot(δ -> pdf(q, [μ[1], μ[2], δ]), 0.0, 0.003; xlabel="δ", legend=false)