This is an author-deposited version published in: http://oatao.univ-toulouse.fr/
Eprints ID: 6859
To link to this article: DOI:
10.1121/1.4757740
URL:
http://dx.doi.org/10.1121/1.4757740
To cite this version:
Laplanche, Christophe Bayesian three-dimensional
reconstruction of toothed whale trajectories: Passive acoustics assisted
with visual and tagging measurements. (2012) The Journal of the
Acoustical Society of America, vol. 132 (n°5). pp. 3225-3233. ISSN 0001-
4966
Open Archive Toulouse Archive Ouverte (OATAO)
OATAO is an open access repository that collects the work of
Toulouse researchers and
makes it freely available over the web where possible.
Any correspondence concerning this service should be sent to the repository
administrator: staff-oatao@listes.diff.inp-toulouse.fr
Bayesian three-dimensional reconstruction of toothed whale
trajectories: Passive acoustics assisted with visual and tagging
measurements
Christophe Laplanche
a)
Universit
e de Toulouse, INP, UPS, CNRS, EcoLab (Laboratoire Ecologie Fonctionnelle et Environnement),
ENSAT, Avenue de l’Agrobiopole, 31326 Castanet Tolosan, France
The author
describes and evaluates a Bayesian method to reconstruct three-dimensional toothed
whale trajectories from a series of echolocation signals. Localization by using passive acoustic data
(time of arrival of source signals at receptors) is assisted by using visual data (coordinates of the
whale when diving and resurfacing) and tag information (movement statistics). The efficiency of
the Bayesian method is compared to the standard minimum mean squared error statistical approach
by comparing the reconstruction results of 48 simulated sperm whale (Physeter macrocephalus)
trajectories. The use of the advanced Bayesian method reduces bias (standard deviation) with
respect to the standard method up to a factor of 8.9 (13.6). The author provides open-source
software which is functional with acoustic data which would be collected in the field from any
three-dimensional receptor array design. This approach renews passive acoustics as a valuable tool
to study the underwater behavior of toothed whales.
I. INTRODUCTION
Researchers have brought three main approaches into
play
to explore the behavior of toothed whales in the field:
Visual, electronic tagging, and passive acoustics. Visual
methods use photo-identification to differentiate individuals,
map their surface movements, and catalogue their clustering
preferences ( Whitehead, 2003, pp. 206–285). Electronic tag-
gin
g consists of attaching embedded systems on whales and
record information on their subsequent behavior. Embedded
systems can contain diverse receptors (acoustic, accelerome-
ter, GPS, etc.) and provide as diverse information on whale
behavior (Johnson et
al., 2009). Passive acoustics consists of
recording whale sounds from dragged, hull-mounted, or
bottom-mounted receptors and real-time or post-process sig-
nals (
Houegnigan et
al., 2010; Miller and Dawson, 2009;
Nielsen and Mohl, 2006). Toothed whales profusely use
sound for communication and echolocation. All toothed
whale species probe their underwater environment by emit-
ting a series of transient, directive, high level clicks (
Madsen
and Wahlberg, 2007
). Passive acoustic outcomes go beyond
inferenc
e on whale acoustic behavior (
Teloni et al., 2008).
Passive acoustics also leads to: Source detection in ambient
noise (Sanchez–Garcia et al., 2010), separation of multiple
phonating individuals (Baggenstoss, 2011; Caudal and Glo-
tin, 2008), localization (Cranch et al., 2004; Hayes and Mel-
linger, 2000; Wahlberg et al., 2001), inference on whale
morphometry (
Growcott et al., 2011), and information on
swim orientation during predation (Laplanche et al., 2005,
2006; Nosal and Frazer, 2007). Each of the three latter
approaches has pros and cons, by providing distinct pieces
of information, with various equipment budget and time
cost, and with different degrees of contact with whales.
Whatsoever, all three approaches share a common feature:
The need to localize the whale as a fundamental step in
studying its behavior.
Passive acoustic localization is achieved by triangulat-
ing source signals on a synchronized array of receptors. Var-
ious designs of receptors have been operated to study the
behavior of toothed whales: One-dimensional (Thode et
al.,
2002), two-dimensional (Thode, 2004), or three-dimensional
arrays (Cranch et al., 2004; Hayes and Mellinger, 2000;
Wahlberg et al., 2001). Since times of emission of whale sig-
nals are unknown, triangulation is not achieved directly by
using times of arrival (TOA) at receptors, but by using times
of arrival differences (TOADs) at pairs of receptors. TOADs
are later processed with statistical software to compute loca-
tion estimates.
One popular option is to derive a whale tra-
jectory from minimum mean squared error (MMSE)
estimates. The drawback of the latter approach is a high sen-
sitivity to measurement errors resulting in broad inaccuracy
(bias) and uncertainty (variance) on the estimate (
Spies-
berger
, 2001; Wahlberg et al., 2001). Bias and variance can
be as large as to make localization results unhelpful.
The author presents an advanced statistical method of
processing TOAD data which aim is to compute localization
results of enhanced quality, that is to say of lower bias and var-
iance than the standard method. The essence of the advanced
method is to refine the processing of the acoustic data and to
use nonacoustic data. Refined processing of the acoustic data
will be achieved by reconstructing a whale trajectory while,
and not afterwards, processing acoustic data. The interest of
using nonacoustic data is to further improve the localization
a)
Author to whom correspondence should be addressed. Electronic mail:
laplanche@gmail.com
procedure (Davis and Pitre, 1995; Laplanche, 2007; Tiemann
et al., 2006
). The complexity of the statistical model (high
number of unknowns, non-linearities, heterogeneous sources
of data) prevents the use of standard statistical tools but Bayes-
ian modeling (
Congdon, 2003, pp. 1–457). Bayesian modeling
h
as already proven to be an efficient framework to address
advanced issues in passive acoustic localization (
Dosso and
Wilmut, 2011; Spiesberger, 2005; Tollefsen and Dosso, 2010).
The essence of Bayesian modeling is to (i) express (known)
measured variables as functions of (unknown) latent variables,
(ii) assign a prior distribution to the latent variables, (iii) calcu-
late a mathematical expression of the posterior distribution of
the latent variables, and (iv) use numerical methods to com-
pute posterior estimates of the latent variables. The mathemati-
cal expression of the posterior distribution as well as the
computation of posterior estimates are complex with complex
models, making Bayesian modeling difficult to apply for
researchers who are not familiar with computer programming
and Bayesian statistics. Recent user-friendly Bayesian model-
ing tools, however, such as BUGS (Bayesian inference Using
Gibbs Sampling), automatically calculate a mathematical
expression and simulate the posterior [steps (iii) to (iv)], leav-
ing only model formulation to users [steps (i) to (ii)], making
Bayesian modeling more accessible (Ntzoufras, 2009, pp. 83–
1
50).
The author first presents the standard localization proce-
dure. This method is reformulated into a Bayesian context,
before being extended, to get up to the full Bayesian local-
ization method. The efficiency of both methods is compared
by using simulated data, with the sperm whale (Physeter
marcocephalus) as an example. The full Bayesian localiza-
tion method could be of interest to study other toothed whale
species, which is discussed.
II. MATERIALS AND METHODS
A.
The standard statistical model
Sperm whales routinely undertake sever al hundred
meters deep, 30 to 60-min dives (Whitehead, 2003, pp. 78–
84, 156–168) interrupted by 10-min breathing breaks at the
sea surface (Watwood et al., 2006). Sperm whales are clearly
visible and identifiable when breathing. Let us consider the
full dive of some sperm whales, diving at time t
dive
and
resurfacing at time t
resurf
. Sperm whale underwater acoustic
activity is recorded on a synchronized array of receptors; let
R be the number of acoustic receptors and r [ {1, … , R} be
an index over receptors. Sperm whales emit several thou-
sands of clicks through their dive; let us consider only a sub-
sample of these clicks, where K denotes the number of
processed clicks and k [ {1, … , K} an index over clicks.
Let t
k
be the time of emission of click k, t
dive
< t
1
< < t
K
< t
resurf
. Let ðx
h
r;k
; y
h
r;k
; z
h
r;k
Þ be the Cartesian coordi-
nates of receptor r at time t
k
. The measured value of TOA of
click k on receptor rðr 2 f1; …; RgÞ is denoted TOA
r,k
. By
using the first receptor as a baseline, the measured value of
TOAD of click k on receptor r ðr 2 f2; …; RgÞ is denoted
TOAD
r,k
.
Let M
k
¼ ðx
k
; y
k
; z
k
Þ be the Cartesian coordinates of the
whale at time t
k
. By using the spherical propagation model, c
is the sound speed, the predicted value of TOA of click k on
receptor rðr 2 f1; …; RgÞ is
TOA
r;k
¼ t
k
þ
1
c
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
ðx
k
ÿ x
h
r;k
Þ
2
þ ðy
k
ÿ y
h
r;k
Þ
2
þ ðz
k
ÿ z
h
r;k
Þ
2
q
:
(1)
And by using the first receptor as baseline, the predicted
value
of TOAD of click k on receptor rðr 2 f2; …; RgÞ is
TOAD
r;k
¼ TOA
r;k
ÿ TOA
1;k
: (2)
The MMSE point estimate of (M
1
, … , M
K
) minimizes the
quadratic sum of the residuals S
2
¼
P
K
k¼1
S
2
k
where
S
2
k
¼
X
R
r¼2
ðTOAD
r;k
ÿ
TOAD
r;k
Þ
2
: (3)
The minimum of S
2
is actually reached by minimizing each
S
2
k
separately. The implications of the latter assertion are first
that the 3K-d imension optimization problem (minimizing S
2
)
can be handled without difficulty by partitioning it into K
three-dimension optimization problems (minimizing S
2
k
for
all k) and using a standard optimization method. Second,
acoustic data at time t
k
is only used to compute an estimate
of the location of the whale at the exact same time. In view
of whale inertia, measured values TOAD
r,k
actually contain
information on the location of the whale around time t
k
. The
advanced Bayesian model wh ich will be presented later will
process acoustic data in this perspective. The advanced
model is an extension of a Bayesian reformulation of the
standard model, which is presented below.
B. Bayesian formulation of the standard model
MMSE estimates are approximately equal to expecta-
tion a
posteriori (EAP) estimates by using a Bayesian statis-
tical model with independent, normally distributed residual
errors of equal variance and vague priors (Appendix
A).
Such
a statistical model is defined by Eqs.
(1) and (2) plus
following Eqs. (4) to (6). Let Normal(l , r
2
) denote some
Normal variate of expectation l and variance r
2
and
Gamma(a, b) some Gamma variate of shape a and scale b.
Measured TOAD values are modeled as independent, nor-
mally distributed variates of expectation predicted TOAD
values and of variance r
2
s
TOAD
r;k
Normalð
TOAD
r;k
; r
2
s
Þ; (4)
where r
2
s
is
the variance of the TOAD residual error. Vague
priors are as signed to the Cartesian coordinates of the whale
x
k
; y
k
; z
k
Normalð0; 10
8
Þ; (5)
and to the variance of the residual error
1=r
2
s
Gammað10
ÿ3
; 10
ÿ3
Þ: (6)
See Appendix A for a mathematical expression of the poste-
rior
of this model. Relationships between model variables
are illustrated with a Directed Acyclic Graph (DAG, Fig. 1).
The full model is an extension of the standard model by
including a trajectory model, visual data, speed statistics,
and acceleration statistics.
C. Whale trajectory
The underwater movement of the whale is modeled as a
contin
uous series of segments of uniform linear motion and
of equal duration (Fig.
2). Let I be the number of segments,
D
I
¼ (t
resurf
ÿ t
dive
)/I be the duration of the segments, and
t
i
¼ t
dive
þ iD
I
ði 2 f0; ::::; IgÞ. Let M
i
¼ ðx
i
; y
i
; z
i
Þ be the
location of the whale at time t
i
, and [M
i
M
iþ1
] be the seg-
ments which make up the modeled trajectory. With this
model, the location of the whale at time t
k
is M
k
¼
ðx
k
; y
k
; z
k
Þ with
x
k
¼ x
i
k
þ
x
i
k
þ1
ÿ x
i
k
t
i
k
þ1
ÿ t
i
k
ðt
k
ÿ t
i
k
Þ; (7)
and similar formulas for y
k
and z
k
, where i
k
2 f0; …; I ÿ 1g is
the index of the trajectory segment where the whale is located
at time t
k
. Predicted TOAD at time t
k
is still provided by Eqs.
(1) and (2) and the relationship between measured and pre-
d
icted TOAD values is still given by Eq. (4). Relationships
between model variables are illustrated in Fig. 1. Free parame-
t
ers of the model are the Cartesian coordinates ðx
i
; y
i
; z
i
Þ as
well as the variance of the residual error r
2
s
. A vague prior is
assigned to the Cartesian coordinates of the whale
x
i
; y
i
; z
i
Normalð0; 10
8
Þ; (8)
and a vague prior is assigned to the variance of the residual
error
[Eq. (6)].
FIG. 1. DAG of the full model.
Frames indicate levels: Receptor (r [
{1, … , R}), click ðk 2 f1; …; KgÞ,
and trajectory segment (i [ {0, … ,
I}). White rectangles: Latent varia-
bles; filled rectangles: observed
variables; circles: model components.
The standard model reduces to model
component 1. The standard model
connects whale coordinates ðx
k
; y
k
;
z
k
Þ,
receptor coordinates ðx
h
r;k
; y
h
r;k
;
z
h
r;k
Þ,
and acoustic data TOAD
r,k
to
each other. The full model is an
extension of the standard model by
adding a trajectory model (model
component 2), visual data (component
3), speed statistics (component 4), and
acceleration statistics (component 5).
FIG. 2. The whale movement is modeled as a series of segments of uniform
linear motion. The trajectory model is represented at two time steps, infini-
tesimal (full line, white and black dots) and at a larger time step (dashed
line, black dots only). At an infinitesimal time step, M
½j
i
¼ ðx
½j
i
; y
½j
i
; z
½j
i
Þ; V
½j
i
¼ ðdx
½j
i
=dt; dy
½j
i
=dt; dz
½j
i
=dtÞ and A
½j
i
¼ ðd
2
x
½j
i
=dt
2
; d
2
y
½j
i
=dt
2
; d
2
z
½j
i
=dt
2
Þ are
the
location, the speed, and the acceleration of the whale at time t
½j
i
. At a
larger time step, M
i
¼ ðx
i
; y
i
; z
i
Þ is the location of the whale at time t
i
,
V
i
¼ (dx
i
/dt, dy
i
/dt, dz
i
/dt) is the average speed of the whale on segment i,
and A
iþ1
¼ (d
2
x
iþ1
/dt
2
, d
2
y
iþ1
/dt
2
, d
2
z
iþ1
/dt
2
) is the acceleration of the whale
when passing from segment i to segment i þ 1. The modeled location of the
whale at click time t
k
by using the trajectory model with a non-infinitesimal
time step are also illustrated (M
k
, gray circles). Mathematical relationships
between M
i
; V
i
; A
i
; M
½j
i
; V
½j
i
; A
½j
i
are provided in Appendix C. The expression
of M
k
is given in the text.
D. Visual data
Sperm whales are clearly visible and identifiable
when breathing (
Whitehead, 2003, pp. 206–285). Let
ðx
dive
; y
dive
; 0Þ and ðx
resurf
; y
resurf
; 0Þ be the measured values
of the Cartesian coordinates of the whale at time t
dive
and t
re-
surf
. Deviations between predicted and measured values are
tolerated and are modeled as independent normally distrib -
uted errors of variance r
2
xy
x
dive
Normalðx
0
; r
2
xy
Þ;
(9)
and similar formulas for y
dive
, x
resurf
, and y
resurf
. The pre-
dicted depths when diving and resurfacing are forced to be
exactly equal to zero, z
0
¼ 0 and z
I
¼ 0.
E. Speed statistics
Sperm whales initiate and end dives by being silent and
by
swimming substantially vertically (Watwood et al.,
2006
). Let V
i
¼ (M
iþ1
ÿ M
i
)/D
I
be the average speed of the
whale on segment i (Fig. 2). The Cartesian coordinates of
V
i
are denoted (dx
i
/dt, dy
i
/dt, dz
i
/dt) with dx
i
=dt ¼
ðx
iþ1
ÿ x
i
Þ=D
I
(and similar formulas for dy
i
/dt and dz
i
/dt,
i [{0, … , I ÿ 1}). Let {0, … , i
start
ÿ 1} and {i
stop
ÿ 1, … ,
I ÿ 1} be the index of segments while the whale is silent at
the beginning and at the end of the dive, respectively. The
horizontal speed of the whale for i 2 f0; :::; i
start
ÿ 1g [
fi
stop
ÿ 1; :::; I ÿ 1g is modeled as independent, normally
distributed variates of expectation 0 and of variance r
2
t
=D
I
dxi
dt
;
dyi
dt
Normalð0; r
2
t
=D
I
Þ; (10)
where r
2
t
is the variance of the horizontal speed of the whale
which would be measured by using a speed recording device
at a 1 s time step (see Appendix C).
F. Acceleration statistic s
Sperm whale acceleration is limited due to hydrody-
namic
drag (Miller et al., 2004). Let A
i
¼ (V
i
ÿ V
iÿ1
)/D
I
be
the acceleration of the whale when passing from segment
i ÿ 1 to segment i (Fig. 2). The Cartesian coordinates of A
i
are noted ðd
2
x
i
=dt
2
; d
2
y
i
=dt
2
; d
2
z
i
=dt
2
Þ with d
2
x
i
=dt
2
¼
ðdx
i
=dt ÿ dx
iÿ1
=dtÞ=D
I
(and similar formulas for d
2
y
i
=dt
2
and d
2
z
i
=dt
2
; i 2 f1; …; I ÿ 1gÞ. The acceleration of the
whale is modeled as independent, normally distributed vari-
ates of expectation 0 and of variance r
2
a
=D
I
d
2
x
i
dt
2
;
d
2
y
i
dt
2
;
d
2
z
i
dt
2
Normalð0; r
2
a
=D
I
Þ;
(11)
where r
2
a
is
the variance of the acceleration of the whale
which would be measured by using an acceleration recording
device at a 1 s time step (see Appendix C).
G. Dataset
The author compares the efficiency of the standard and
the
full Bayesian methods with a simulated dataset. The
interest of using a simulated dataset is to have at one’s dis-
posal true values, and compare them to estimated values.
The author considers 48 simulated whale trajectories. An
example of a trajectory is illustrated in Figs. 3 and 4. The
wha
le dives at some arbitrary point (x
dive
¼ 500, y
dive
¼ 0 m)
at t
dive
¼ 0, starts clicking at t
start
¼ 2 min, stops clicking at
t
stop
¼ 30 min, and resurfaces at t
resurf
¼ 40 min. Trajectories
are randomly generated in accordance with the autoregres-
sive model of Appendix C with r
v
¼ 0.1 m/s, r
a
¼ 0.05 m/s
2
,
and D
t
¼ 1 s.
Sperm whale clicks are recorded on three hydrophones.
Hydrophones are initially located at (0, 0, ÿ30) m, (ÿ200, 0,
ÿ40) m, and (0, 200, ÿ50) m. All hydrophones drift North-
east at 0.1 m/s. Both direct and surface-reflected source sig-
nals are assumed to be detected on the receptors, leading to a
virtual array of six hydrophones (Skarsoulis and Kalogera-
kis,
2005
). The author investigates the consequences of the
variations of the quantity and the quality of the acoustic
dataset by considering two click rates as well as two noise
levels. Two series of predicted TOAD values are computed
for each trajectory, at a slow rate (1 click every D
K
¼ 30 s,
K ¼ 57) and at a high rate (1 click every D
K
¼ 5 s, K ¼ 337).
A white Gaussian noise of standard deviation r
s
¼ 0.1 ms or
r
s
¼ 1 ms is added to the predicted TOAD values, leading to
4 acoustical datasets for each trajectory. Furthermore, two
levels of trajectory resolution are compared by using the full
Bayesian model: Low-re solution trajectories, with D
I
¼ 60 s
segments (I ¼ 40), and smoother trajectories, with D
I
¼ 10 s
segments ( I ¼ 240). Acoustic data at a slow click rate are
processed with the standard model and with the full model at
a low-resolution and acoustic data at a high click rate are
processed with the standard model and with the full model at
a high-resolution. As a summary, a total of 384 simulations
are carried out: 48 trajectories, 2 noise levels, 2 click rates,
and 2 models.
H. Softwa re
Models are implemented in BUGS language by using
Ope
nBUGS, an open source version of WinBUGS (
Ntzouf-
ras,
2009
, pp. 1–492). The creation of input files for BUGS,
as well as the gathering of BUGS output files in order to
compute trajectory statistics and display, is achieved with R.
BUGS and R scripts are gathered within the open-source
software SBPLAsH version 2.0 (http://modtox.myftp.org/
software/sbplash). Users can provide input files and explore
simulation results through SBPLAsH graphical user inter-
face. SBPLAsH also creates Unix batch and portable batch
system scripts to perform parallel BUGS computations on a
UNIX desktop computer or a high performance computing
(HPC) resource. See Appendix
B for more computational
details.
I. Model comparisons
Models are compared in terms of goodness-of-fit and
comple
xity (Appendix B) as well as accuracy and uncer-
tainty. The average absolute bias
D
x
¼
P
K
k¼1
jx
k
ÿ
^
x
k
j=K
(with similar formulas for
D
y
and
D
z
) is used as a proxy of
model accuracy, where
^
x
k
denotes the estimate of x
k
. The
average standard deviation
R
x
¼
P
K
k¼1
^
r
x;k
=K(with similar
