I like Sonars. I like computers. I like science. I'm a geek

Today was my last day at Reson, and I’m finally getting around to writing the Python version of my ray tracing script.  What I want to do is take sound speed profiles as input, and compute the difference in depth that would result from using one versus the other at different multibeam launch angles.  It’s just a little tool that I’d like to have to figure out how much error I could expect given two profiles.

So far, I’ve read in the profile (from a text file), and parsed the data:

filename = 'ssp.txt'

ssp = []
depth = []

file1 = open(filename)
for line in file1.readlines():
#    print line
    splitup = line.split()
    ssp.append(float(splitup[1]))
    depth.append(float(splitup[0]))

After this, I defined an array of launch angles (relative to horizontal, so +/- 75 degrees was actually 15 – 165), and an array of sound speed values interpolated from the sound speed profile to a certain depth increment (1m), and down to a given water depth (10m).

LaunchAngle = r_[15.:166.:5.] # array of launch angles in degrees
LaunchAngleRadians = LaunchAngle*pi/180
nadirdepth = 100 # nadir depth in meters
DepthIncrement = 0.5 # depth increment in m

if nadirdepth > max(depth):
    ssp.append(0.016*(nadirdepth - max(depth)))
    depth.append(nadirdepth)

if min(depth)>0:
    ssp = insert(ssp,0,ssp[0])
    depth = insert(depth,0,0)

depthvector = r_[0:nadirdepth:DepthIncrement]
interpfun = interp1d(depth,ssp)
sspvector = interpfun(depthvector)

And as if that much wasn’t suspicious enough, I decided that I wanted to do all of the calcs using matrix math. I figured I’d ditch all the loops that were in my original Matlab script, and replace them with some nice streamlined matrices. So efficient. The problem is that there are some exceptions that I need to take care of, that are difficult to deal with when doing everything in a series of matrices. For example, when the sound speed in a certain layer is the same as the sound speed in the first layer – then the gradient is zero, and the radius of curvature is infinite. Here’s what I’ve hacked together so far. It’s not pretty.

thetaLaunch = tile(LaunchAngleRadians,(size(depthvector),1))
depthmat = transpose(tile(depthvector,(size(LaunchAngleRadians),1)))
sspmat = transpose(tile(sspvector,(size(LaunchAngleRadians),1)))

g = (sspmat - sspvector[0]) / DepthIncrement

Rc = -sspvector[0]/(g*cos(thetaLaunch))

thetamat = arccos( ((Rc*cos(thetaLaunch)) - DepthIncrement ) / Rc )

dx = (Rc*sin(thetamat) - Rc*sin(thetaLaunch))
dxnandex = isnan(dx)
dx[dxnandex] = DepthIncrement/tan(thetaLaunch[dxnandex])
dx[0,] = 0
dxsum = cumsum(dx,axis=0)

Yuck. But at least I managed to get some rays that made sense.

The next step is to convert the ray paths to two-way travel times, and then compare the depth that would be calculated if the ray tracing were done with a second sound speed profile.

One of the most frustrating things that I run into when processing multibeam echosounder data is a smiling seafloor.  When the sound speed profile is incorrect, the computed seafloor becomes curved either up or down.  This is because the effect of an incorrect profile is more pronounced on the longer ranges (outer beams) than at nadir.  So the data is either smiling or frowning.  I should also note that the smiling/frowning could arise from incorrect beam steering angle, which can also be caused by having the wrong sound speed reading at the sonar head.

I spend far too much time messing around with sound speed profiles during post processing.  In the ocean, the speed of sound profile can change rapidy – both spatially and temporally.  If the profile changes while you’re surveying, your data will be wrong.  An interesting paper published by some of my favorite people can be found here:

IHO paper:  Estimation of Sounding Uncertaintly from Measurements of Water Mass Variability

So a lot of times I’ll end up having a profile that’s not quite right, and having to tweak it in my processing software, which takes a really long time.  I read an interesting paper recently by some folks at the Acoustic Remote Sensing Group in Delft University of Technology:

Underwater Acoustic Measurements Conference 2009: An efficient method for reducing the sound speed induced errors in multibeam echosounder bathymetric measurements

This paper discusses using the overlapping area between parallel lines to obtain an estimate of the sound speed throughout the water column.  The authors assumes that in shallow water the sound speed profile can be approximated by a single value.  The differences between overlapping lines are minimized using a Gauss-Newton method.

I’d like to look into this in more detail, but also maybe look at approaching the problem in a slightly different way – what if we looked at the overlap between crosslines instead of parallel lines?  The sound speed smiley error doesn’t show up in the along-track direction.  So if you did a beam-by-beam analysis of one line versus the perpendicular line, you would see the sound speed profile pop out.  If you average the difference in depth between each of the beams and an interpolated position on the crossline, you would end up with an average representation of the smile…  Okay, I might need to think about this a bit more.  But I’m imagining somehow running your entire survey, and then running a crossline across the whole thing, and computing profiles over the entire area.  And then maybe doing a whole checkerboard grid – type survey, and doing a full blown analysis on spatial and temporal sound speed variation.  I may be getting ahead of myself a bit.  But I think it’s worth pondering…