Channel geometry from impulse response

Reference

This tutorial is adapted from Example C presented in:

Problem statement

Consider a scenario where a bottom-mounted sensor is deployed on a sub-surface mooring at an unknown altitude over the seabed. The exact location of the sensor is unknown, but we know the general area where it is deployed. The sensor is equipped with an acoustic transponder that we can query from a surface unit when we are within a 250 m range of it. We deploy the transducer of the surface unit from the boat on a 7 m rope with a weight attached. Due to currents, the transducer is not hanging perfectly vertically, and so it’s exact depth d1 is not known. The depth sounder of our boat tells us that the water depth h is 20 m, and the Captain of the boat assures us that the bathymetry is quite flat. We query the transponder and get a response. The broadband acoustic response from the transponder allows us to estimate the delays of 4 multipath arrivals. We wish to estimate the depth d2 of the sensor and the range r between the boat and the sensor using the multipath arrival delays.

Dataset

To illustrate this idea, let us generate a synthetic dataset with d1 = 7.2 m, d2 = 12.7 m, h = 20 m, and r = 97.3 m. Since we have an range-independent iso-velocity environment, we can use the PekerisRayTracer (otherwise we could use the RaySolver):

using UnderwaterAcoustics

function model(h, r, d1, d2)
  env = UnderwaterEnvironment(bathymetry=h)
  pm = PekerisRayTracer(env; max_bounces=4)
  tx = AcousticSource(r, -d2, 1000.0)
  rx = AcousticReceiver(0.0, -d1)
  arr = arrivals(pm, tx, rx)
  t = getfield.(arr, :t)        # get a.t for each arrival a in arr
  (t[2:end] .- t[1]) .* 1000    # relative time in milliseconds
end

data = model(20, 97.3, 7.2, 12.7)
8-element Vector{Float64}:
  1.2078734494638532
  1.2340398375722539
  6.4705086705396235
  3.755897382140966
 10.919625105605624
 10.987838758289167
 20.841074171559228
 16.304090813082887

Probabilistic model

We can model this as a Bayesian inference problem using Turing.jl and some reasonable priors:

using Turing
using Random

@model function pmodel(data)
  h  ~ Normal(20.0, 0.1)
  r  ~ Uniform(0.0, 250.0)
  d1 ~ Normal(7.0, 1.0)
  d2 ~ Uniform(0.0, 20.0)
  μ = model(h, r, d1, d2)
  data ~ MvNormal(μ, 0.1)
end

m = pmodel(data)
q,_ = vi(MersenneTwister(1), m, q_meanfield_gaussian(m), 10000)

# extract mean values of parameters h, r, d1 and d2
μ = mean(rand(q, 10000); dims=2)
4×1 Matrix{Float64}:
 20.009226482010614
 97.42325326194218
  7.200312567411043
 12.704686858499537

We see that the estimated values of h, r, d1 and d2 are fairly close to the ground truth values. We can also plot the uncertainty in the estimates:

using StatsPlots

plot(h -> pdf(q, [h, μ[2], μ[3], μ[4]]), 18, 22; xlabel="h", legend=false)
plot(r -> pdf(q, [μ[1], r, μ[3], μ[4]]), 50, 150; xlabel="r", legend=false)
plot(d1 -> pdf(q, [μ[1], μ[2], d1, μ[4]]), 5, 9; xlabel="d1", legend=false)
plot(d2 -> pdf(q, [μ[1], μ[2], μ[3], d2]), 0, 20; xlabel="d2", legend=false)