2-D Tracking of Acoustic Sources from a Single Sensor

CAUTION! COGITATION IN PROCESS!
The following post will cover 2-D tracking of acoustic sources from a single acoustic sensor as suggested by Cato 1998. This concept is something that is often mentioned in theoretical instances but I have yet to see it applied. Below are my current thoughts on the matter and the issues I've run across while attempting to implement the method.

Passive acoustic monitoring is one way to  non-invasive study acoustically active animals (frogs, bats, snapping shrimp, codfish, toddlers etc.).

One of the most exciting things about passive acoustic monitoring  is the ability to track individual animals. What can be learned from acoustic tracking?

• Habitat preference
• Effects of anthropogenic (man-made) sounds
• Number of animals present
• Social affiliations (which animals spend more time together)

From a conservation standpoint, determining the total number of animals within a population is key in making management decisions and allocating limited resources (e.g. $$$). Sadly, there are not the resources to monitor the whole ocean, or even a small part of it at once. 

What is a scientist to do? As with anyone working on a limited budget, you make your data stretch by utilizing all the methods available. 

Once such way to get more out of your data is to use the first arrival and second arrivals (or the sound and it's echo) to gain meaningful information about the distance from the source to the receiver.

Using the time difference between the first arrival and second arrival as well as the difference in the relative amplitudes (volume) between the first and arrival and echo it is supposed to be possible to calculate the distance to and the depth of the caller, assuming you know how deep your recording device is. 

Adapted from Au and Hastings Principles of Marine Bioacoustics

For anyone interested here's how the math works.

Known variables:

τ- time difference of arrival between first arrival (R1) and echo (R2+R3)

D_hydrophone- How deep the hydrophone or acoustic detector is. In our case it's assumed to be on the ocean floor

I₁- Intensity of the first arrival (R1)

I₂- Intensity of the second arrival or the surface echo (R2+R3)

k-sqrt(I₁/I₂)

C- sounds speed in the water

Unknown variables: 

X-   Horizontal distance from source to hydrophone

R1- azimuthal distance from source to receiver

Dsource- Depth of the source

Assumptions:

1. The first recorded arrival is the direct path from the source to the receiver
2. The second recorded arrival is the echo from the surface
3. The surface is perfectly reflective (thought this can be modified)
4. K does not approach 1 
5. τ does not approach 0
6. The speed of sound (C) is uniform

Equations:

R1=(C*τ)/(k-1);

X=D_hydrophone-(R1.^2.*(k.^2-1))./(4*D_hydrophone);

Dsource=sqrt(R1^2-X^2);

Limitations:

Even under ideal circumstances, as below. This method has limitations. Very near the sensor and areas where k approaches 1 the equations break down. Below is an example of the errors that would be introduced if the sensor was in 8m of water. 

Range (R1) and source depth (Dsource) error at various distances from the receiver (red dot in lower left corner) 

Questions:

For the life of me I can make it work!  When real data are plugged into these equations, the output estimates are crazy (eg. sources 20m above the surface of the water, unlikely). 

I'm wondering if this may be, in part, due to the directional nature of the sounds I'm looking at. The origional author of this work, Dr. Cato (1998) used "whales" as an example when laying out this method. Believe it or not the term "whales" can be quite confusing. Was he referring to large whales only that tend to make low-frequency pseudo-omnidirectional sounds? Or was he referring to all animals in the cetacean family which includes everything from the blue whale to the tiny tiny (endangered) vaquita. Who's to say?

Interestingly, in Au and Hastings book they use dolphins, which produce highly directional sounds as an example, which suggests that it shouldn't be a problem. 

K-wave approximations:

To investigate whether the directionality of the source might be the issue I used K-wave software to model a dipole source with a reflective surface and measured the pressure at some distance away. 

The results from combining the K-wave pressure field with the above equations were less than ideal. The source position was way off.  Possibly this was due to the scale of the model, which is very tiny (violation of the assumptions) but it's impossible to parse out at the moment.

Pressure output from the sensor in the above k-wave model

Unfortunately, my current computer lacks the requisite *ahem* processing power to produce a scale model (meters rather than mm).

So, that's where things stand. It is presently unclear wither the errors I'm getting in my models are due to 

1. Violations of the assumptions
2. Directionality of the source
3. Something completely different 

If anyone has any brilliant ideas, I'm all ears. Thanks for reading!

---

Code

Matlab code for range error estimates below

function [tao, dis, r1, k, p1, p2]=simple_ray(c,hyd,dol, p, f)

%find the time difference of arrival between the first arrival and the surface

% bounce given speed of sound (c), hydrophone (x, y) and dolphin

% coordinates(x, y) and plot (binary)

if nargin <=4

    p = 0;

    alpha=0;

else

    T=8;

    D=hyd(2)

    [alpha_Db]=alpha_sea(f, T, D)

    alpha=10^(alpha_Db/20)

end

r1=sqrt((hyd(1)-dol(:,1)).^2+(hyd(2)-dol(:,2)).^2);

r2=sqrt((hyd(1)-dol(:,1)).^2+(hyd(2)+dol(:,2)).^2);

A=asin((hyd(2)-dol(:,2))./r1);

dis=r1.*cos(A);

tao=(r2-r1)/c;

p1=1./(r1);

p2=1./(r2);

k=r2./r1;

end

-------------------------------------------------------------------------------

function [r, x ,ds]=source_range(k, tao, dh)

% range from signal to receiver using single hydrophone and 1 surface

% reflection

% k -   Pressure ratio in pascals of the second and first arrival (P1/P2)

% tao-  Time difference between first and second arrival

% dh-   depth of hydrophone

c=1500;

r=(c*tao)./(k-1);

x=dh-(r.^2.*(k.^2-1))./(4*dh);

ds=sqrt(r.^2-x.^2);

end

--------------------------------------------------------------------------------------

%% Script used to  run above functions 

clear all; clc

c=1500;

hyd=[0 10]; % Hydrophone locations 

p=1; % Yes, make the figures

% make 1m by 1m grid

hab=zeros(hyd(2),50);

x=0:size(hab,2)-1;

y=0:size(hab,1)-1;

% Create a structure to store the k (p1/p2), tao, real horizontal distances, real horizontal ranges, range error (real-est range), and distance error

Trcking_mdl=struct('k_mask', hab, 'tao', hab, 'real_x', hab, 'rng', hab, 'rng_error', hab, 'dis_error', hab );

for ii=1:length(x)

     dol_x=ones(1,length(y))*(ii-1);

     dol_y=y;

     dol=[dol_x' dol_y'];

      % Call the simple ray function to get real horizontal distance, time difference of arrival (tao) and k (p1/p2) with input variables of speed of sound (c), hydrophone location (hyd), source/dolphin location (dol) and binary about whether to kick back figures (p)

     [tao, dis, r_true, k]=simple_ray(c,hyd,dol, p);

      % Fill in the the tracking model structure

     Trcking_mdl.rng(:,ii)=r_true;

     Trcking_mdl.k_mask(:,ii)=k;

     Trcking_mdl.real_x(:,ii)=dis;

     Trcking_mdl.tao(:,ii)=tao;

      % With known tao, k calculate the estimated range (r), horizontal distance (x_est) and source depth estimate (ds_est) using Cato 1998 method

     [r, x_est ,ds_est]=source_range(k, tao, hyd(2));

     % Calculate the error between the true range and the range obtained by Catos' method

     Trcking_mdl.rng_error(:,ii)=abs(r-r_true);

   
     % Calculate the horizontal distance error by taking the distance between the known source distances and the source distance obatined by Cato's method

     Trcking_mdl.dis_error(:,ii)=abs(dis-ds_est);

end

 %  Make some pictures!

figure (1)

subplot(2,1,1)

hold on

contourf(flipud(Trcking_mdl.rng_error), 20)

scatter(1,1, 'r', 'filled')

colormap (flipud(hot))

colorbar('EastOutside')

xlabel('Distance [m]')

ylabel('Depth [m]')

title('Range Error')

subplot(2,1,2)

hold on

contourf(flipud(abs(Trcking_mdl.dis_error)), 20)

scatter(1,1, 'r', 'filled')

colorbar('EastOutside')

xlabel('Distance [m]')

ylabel('Depth [m]')

title('Horizontal Distance Error')

figure(2)

subplot(2,1,1)

hold on

contourf(flipud(abs(Trcking_mdl.k_mask)), 20)

scatter(1,1, 'r', 'filled')

colorbar('EastOutside')

xlabel('Distance [m]')

ylabel('Depth [m]')

title('K (P1/P2)')

subplot(2,1,2)

hold on

contourf(flipud(abs(Trcking_mdl.tao)), 20)

scatter(1,1, 'r', 'filled')

colorbar('EastOutside')

xlabel('Distance [m]')

ylabel('Depth [m]')

title('Tao (sec)')