Tracking the kinematics of caudal-oscillatory swimming: a comparison of two on-animal sensing methods

Many marine mammals and fish propel themselves with oscillations of the posterior body (termed caudal-oscillatory swimming which includes thunniform and sub-carangiform swimming styles). These undulations generate angular displacements of the body (body rotations, BRs) as well as specific accelerations (SAs), i.e. accelerations with respect to the body frame of the animal (Fish et al., 2003), both of which are measured by accelerometers in surface-attached tags (Johnson et al., 2009). The small size and low power consumption of micro-electro-mechanical system (MEMS) accelerometers have resulted in their widespread use to study animal locomotion, but their unavoidable sensitivity to both orientation and SA makes the resulting data difficult to interpret. In steady swimming, the mean body posture typically varies slowly relative to the stroking rate (Sato et al., 2007), allowing body posture to be removed from the accelerometer signal by high-pass filtering (Sato et al., 2003), resulting in the so-called dynamic acceleration (DA; Wilson et al., 2006). This filtered signal has been used both as a proxy for SA (Sato et al., 2011) and to calculate the overall dynamic body acceleration (ODBA), a widely used integrative indicator of activity level (Wilson et al., 2006; Gleiss et al., 2011). However, DA contains both SA and fast BR (Martín López et al., 2015), leading to biased estimates of animal-generated forces and complicating the relationship between ODBA and mechanical work done. In addition, fine-scale studies of caudal-oscillation swimming patterns and how they respond to different force loading require separate information about the BR (i.e. the instantaneous changes in orientation due to propulsor displacement) and the SA (i.e. the resulting net acceleration).

Inertial navigation systems (INS) in submarines and aeroplanes use triaxial gyroscopes to obtain an independent measure of the instantaneous orientation of the platform, allowing the SA component in the accelerometer signal to be estimated (Chatfield, 1997). Although MEMS rate gyroscopes are beginning to be used in animal tags (Noda et al., 2012) these have several drawbacks. MEMS gyroscopes measure the force needed to rotate a vibrating frame out of its plane of oscillation, producing a signal that is proportional to the turn rate. As sensitivity depends on vibration amplitude, there is an inherent trade-off between precision and power consumption. Gyroscopes require a stable oscillatory movement, meaning that they cannot be efficiently turned off between measurements, a technique often used with other sensors to save power. As orientation is derived by integrating the gyroscope signals, complex processing is needed to combat drift in the orientation estimates due to integration of sensor noise and temperature-induced changes in calibration (Luinge and Veltink, 2005). Drift can be combated both by periodic correction using another orientation sensor and by processing the sensor signals with a Kalman filter (Noda et al., 2012), non-linear estimator (Fourati et al., 2011b) or complementary filter (Fourati et al., 2011a) to track and remove the drift. The independent orientation estimate is typically obtained from accelerometers and magnetometers during intervals in which SA is presumed to be minimal, e.g. gliding (Fourati et al., 2009; Noda et al., 2012).

High sampling rate magnetometers can be used in place of gyroscopes to estimate BR in swimming, enabling an estimate of SA when combined with accelerometers (Martín López et al., 2015). Magnetometers have several advantages compared with gyroscopes. They can be turned off between measurements to save power, e.g. a magnetometer sampled at 100 Hz can be turned on for 1 ms to take a measurement and off for 9 ms, reducing power consumption to 10% at the expense of increased sensor noise. Moreover, these sensors are already implemented in many tags to estimate heading and are becoming available, integrated with accelerometers, as a low-power alternative to gyroscopes in consumer products (Delporte et al., 2012). A drawback of magnetometers is that, being a vector field sensor, they only measure two of the three rotational degrees of freedom, restricting the types of movements which they can track. However, magnetometers appear to capture well the largely planar swimming movements of cetaceans (Martín López et al., 2015) and may be applicable to other swimming and flying locomotion styles.

Here, we compare the performance of gyroscopes and magnetometers in concert with accelerometers to estimate kinematic parameters in Blainville's beaked whale (Mesoplodon densirostris) swimming. We first extend the magnetometer processing method developed in Martín López et al. (2015) to produce a comparable algorithm for gyroscopes. This simplified approach explicitly removes slowly varying forces by high-pass filtering the sensor data, and so estimates movements relative to the mean orientation and force balance of the animal. We then use data from a tag with synchronously sampled triaxial magnetometers, accelerometers and gyroscopes on a free-swimming beaked whale to compare the performance of gyroscopes and magnetometers in estimating SA and BR, and contrast the benefits and drawbacks of these sensors.

A DTAG logger was attached with suction cups to a male Blainville's beaked whale in May 2013 off the island of El Hierro, Canary Islands, Spain, following the methods described by Aguilar de Soto et al. (2012). The tag contained triaxial accelerometers (Kionix KXSC7-1050), magnetometers (Honeywell HMC1043) and gyroscopes (Invensense IDG-500 and ISZ500), and a pressure sensor, synchronously sampled at a rate of 200 Hz/sensor channel with 16-bit resolution. The DTAG also contained a temperature sensor located close to the gyroscopes, sampled at 33.3 Hz. The accelerometer, gyroscopes and pressure sensors were powered continuously and signals from these sensors were anti-alias filtered with a single pole filter at 50 Hz before sampling. The magnetometer was turned on for 1.67 ms of every 5 ms to reduce power and no anti-alias filter was used as the sampling rate was well above the rate of the rotational dynamics. Data from all sensor channels except temperature were decimated in post-processing to a 25 Hz sampling rate using identical symmetric finite impulse response (FIR) low-pass filters (10 Hz cut-off frequency). Sensors were calibrated in the laboratory using standard methods ( Appendix 1). The triaxial sensor signals were rotated to correct for the orientation of the tag on the whale, which was estimated at each surfacing from the stereotypical movements during respiration (Johnson and Tyack, 2003). Dives up to and including the second deep dive were analysed as the tag was partially detached thereafter.

The tagging method was approved by the University of St Andrews teaching and research ethics committee. Tagging was authorised by a cetacean research permit from the Spanish Ministry of Environment (MAGRAMA) to N.A.d.S. at Universidad de La Laguna.

In studying locomotion, we are primarily interested in the dynamics of orientation and SA associated with propulsor movement rather than the mean posture of the animal. Assuming that caudal-oscillation only produces signals at and above the stroking rate, these can be separated from slower postural changes by high-pass filtering each axis of A, M and Ω with identical filters to obtain ⁠, and ⁠, where A, M and Ω are the triaxial vectors measured by the accelerometer, magnetometer and gyroscope, respectively. A high-pass filter cut-off frequency of about 0.25–0.5 of the stroking rate (e.g. calculated as in Sato et al., 2007), is typically appropriate. Both and contain information about BR, while reflects both BR and SA. The procedure is to estimate the BR component in using either or and so obtain an estimate of SA from by subtraction (see  Appendix 2 for detailed algorithm descriptions). Using ⁠, the BR is obtained in only one axis, but this is sufficient to track movements in steady swimming (Martín López et al., 2015) if the BR can be assumed to occur primarily around a single axis (for a cetacean, this is the y or transverse axis). To monitor the validity of this assumption, we computed the R² statistic representing the quality of fit between the modelled and observed on a stroke-by-stroke basis (i.e. how well the assumed pitch-only BR fits the observed magnetometer measurements). Using ⁠, all three axes of the BR can be obtained by integration, allowing a more complete description of the motion. Because the integration is followed by high-pass filtering, there is no requirement for an absolute orientation reference or for drift estimation as normally needed with gyroscope processing.

Both the magnetometer and gyroscope methods were used to estimate BR and SA during steady swimming in dive descents and ascents. Glides as defined in Martín López et al. (2015) longer than 0.7 s were removed. A high-pass filter cut-off frequency of 0.4 times the average stroking rate of the analysed data for each dive phase was used. Applied to cetacean swimming, the magnetometer method estimates rotation about the y-axis. Rotations around the other two axes were set to zero for comparison with the gyroscope method. The strength of agreement between the two methods for BR and SA signals was estimated using Lin's concordance correlation coefficient (Lin, 1989, 2000) using epi.cc in the R package ‘epiR’ (R Development Core Team 2009). The similarity (R² statistic and scale factor) of the BR and SA estimates was also evaluated using least-squares linear regressions. Processing was performed in Matlab 7.0 (MathWorks).

Estimates of the BR and SA of a Blainville's beaked whale, calculated using the gyroscope and magnetometer methods, were compared over 47 min of data representing swimming in descent and ascent phases of two deep (>500 m) and three shallow dives (Tyack et al., 2006). Without an independent control, we cannot assess which of the methods was more accurate but we expect that the gyroscope will give a fair approximation of the true BRs, at least under conditions similar to those of bench calibration, whereas the magnetometer will only be able to track a subset of movements. However, the mean R² of the magnetometer method over all analysed strokes was 0.86, indicating that the pitch-only rotation model used in this method was a reasonable approximation to the BR dynamics. Visually, the two methods produced very similar BR and SA estimates, which tracked both the steady swimming and rapid gait changes characteristic of beaked whale ascent swimming (Fig. 1; Martín López et al., 2015). Over the 1085 strokes analysed, the BR estimates had a concordance of 0.97, a mean relative scale factor of 7% and an R² value between the two methods after correcting for the scale factor of 0.94. The scale factor between the magnetometer and gyroscope BR estimates varied with temperature (Fig. S1). In shallow dives and descents of deep dives where the temperature was >16°C and therefore close to the temperature during gyroscope calibration (22°C), the scale factor varied from 2% to 14% but tag temperatures down to 11°C in the ascents of deep dives gave scale factors of up to 22%.

The SA estimates from the two methods were also similar, with concordance of 0.95 and 0.97 for surge and heave SA, respectively (R² of 0.91 and 0.94 and scale factor of 0.96 and 0.97, respectively). The match for sway SA was relatively poor, with a concordance of 0.82 (R²=0.68, scale factor=0.90), probably due to small roll and yaw BRs during stroking that were unmodelled by the magnetometer method and so remain uncorrected in the magnetometer sway SA estimate (Martín López et al., 2015). However, the sway SA calculated with the gyroscope has 23% of the magnitude of heave and surge SA and so represents only 3% of the total SA magnitude. Thus, the absolute differences between the methods are small.

Although neither the magnetometer nor the gyroscope method allows direct measurement of slowly changing forces, the impact of these changes can be seen in caudal-oscillation rates and amplitudes. We used both methods to estimate BR and SA during the initial 30 m of a descent with steady swimming during which buoyancy forces from air in the lungs decrease as a result of increasing hydrostatic pressure causing progressive lung compression (Sato et al., 2002). The methods gave almost identical results in tracking the rapid changes in stroke rate, BR and SA as the whale dived (Fig. 2), demonstrating the utility of these methods in analysing gaits under changing loads. During this descent, speed estimated from the depth rate and pitch angle as in Miller et al. (2004) increased from 0.8 to 1.1 m s⁻¹ as the whale dived away from the surface (a net acceleration of about 0.007 m s⁻²) despite progressively lower stroking rate, BR (a proxy for propulsor amplitude), surge and heave SA, which decreased to 81%, 74%, 30% and 44% of their initial values, respectively (Fig. 2). Thus, despite a small increase in drag due to increased speed, the whale was able to reduce swimming effort as depth increased, due to rapidly decreasing buoyancy as the carried air volume was compressed. Assuming passive isothermal compression, gas volume should decrease as 1/(1+0.1×depth) for depth in metres. Thus, at 30 m depth, the gas volume should be 35% of its initial volume at 4 m, roughly comparable to the observed decrease in SA. The surge and heave SA decreased more linearly with relative air volume (R²=0.97 and 0.96, respectively) than with depth (R²=0.83 and 0.78), suggesting that this animal adjusted caudal-oscillation amplitude and rhythm so as to match its force production to the decreasing net buoyancy during the first 30 m of the descent.

Even though magnetometers and gyroscopes provided comparable estimates of BR and SA in the data examined here, there are several advantages and disadvantages with each of these sensors that help guide their use (Table S1). The primary advantage of triaxial gyroscopes over magnetometers is that rotations in all three axes can be estimated, making it possible to track rapid complex motions such as during escapes from predation (Noda et al., 2014) or foraging. For beaked whales, prey capture attempts are indicated by the presence of buzzes in the tag sound recording (Johnson et al., 2004). Blainville's beaked whales perform rapid manoeuvres while attempting to capture prey (Madsen et al., 2014), including rolls and heading changes, well tracked by the gyroscope but less well by the magnetometer method, leading potentially to large errors in the SA estimates (Fig. 3). The mean R² of the magnetometer method over this interval was 0.5, indicating that a pitch-only BR model is a poor fit and that the magnetometer-derived SA estimates will be inaccurate. Thus, the magnetometer method is most appropriate for caudal-oscillatory aquatic locomotion in which the BR largely occurs around a single axis, e.g. steady swimming of cetaceans, phocid seals and some fish. An additional advantage of gyroscopes is that they measure rotations equally well irrespective of their mean orientation whereas magnetometers measure orientation with respect to the Earth's magnetic field and so cannot detect rotations about an axis parallel to this field vector (Fig. 3; Martín López et al., 2015). For the descent and ascent data analysed here, only ∼0.2% of samples suffered from this unobservability problem. However, this percentage depends both on animal behaviour and on the inclination angle of the magnetic field in the tagging location, and could potentially be a significant limitation in some situations. It is in fact possible to extend the magnetometer method to track two BR angles (e.g. pitch and roll) simultaneously, more closely approximating the behaviour of a gyroscope, but this increases further the range of orientations over which BR is unobservable.

A drawback of gyroscopes for small biologging tags is their power consumption, which is 5–10 times greater than that of magnetometers in currently available components. Given that MEMS gyroscopes require a moving element to sense rotation while magnetometers do not, there may be little prospect for this gap to narrow. For MEMS gyroscopes, frame vibrations of 25–30 kHz are typical and coincide with the sensitive hearing range of many odontocetes (Southall et al., 2007), so gyroscope frame vibrations may also be audible to the animal to which the tag is attached. An additional disadvantage of gyroscopes is that their sensitivity can change considerably with temperature at least in the consumer grade devices used in tags, leading to errors in BR and SA if uncorrected. In comparison, the method used to derive BR from the magnetometer is ratiometric and so is insensitive to the gain of the magnetometer or the magnetic field strength (Fig. S1).

The processing methods presented here involve two simplifications compared with algorithms conventionally used with gyroscopes (Noda et al., 2012; Fourati et al., 2011a,b). The first is that low-frequency signals in the sensor data, such as from slow orientation changes and sensor drift, are removed by high-pass filtering. This bypasses the step of estimating orientation and sensor offset, thereby reducing the number of parameters and eliminating the need for a Kalman or complementary filter. The second simplification is the use of Euler angles instead of quaternions to represent BR. Four-parameter quaternions are used to linearise INS algorithms and avoid gimbal-lock at extreme orientations (Jahanchahi and Mandic, 2014). The orientation dynamics in thunniform and sub-carangiform swimming modes are of small amplitude (e.g. <10 deg for tags on cetaceans) and are centred on zero, allowing us to use the more familiar Euler angles, which represent directly the rotations of an animal around its three canonical body axes. More complex processing may be required for swimming modes such as anguilliform, which involve much larger BRs (Sfakiotakis et al., 1999). A drawback of the high-pass filter processing used here is that slowly varying forces such as mean drag, lift and buoyancy cannot be estimated. Instead, we estimate swimming dynamics relative to the mean orientation and force balance of the animal. Because of this, our SA estimates cannot be integrated to predict swimming velocity and position as in a conventional INS, i.e. they cannot be used to track an animal underwater. However, given the inaccuracy and drift (Fourati et al., 2011b; Shiau et al., 2012) of currently available low-power MEMS gyroscopes, it is questionable whether reliable velocity estimates, suitable for track reconstruction, could be obtained with a complete INS algorithm either.

The techniques developed here provide a consistent way to analyse movement dynamics using either magnetometers or gyroscopes, making it possible to compare results obtained with different tags and sensors. The above observations suggest that gyroscopes are best suited for studying complex motions in short-duration tags or when continuous monitoring of locomotion dynamics in all orientations is essential. The lower power consumption of magnetometers and their relative insensitivity to calibration errors suit them to longer-term tags on animals that have largely planar locomotion and where occasional unobservable intervals can be tolerated. In both cases, the ability to separate orientation dynamics in swimming from the resulting SA enables studies of gait adjustment, force production and swimming energetics. In particular, the estimated SA may be a more consistent and physically meaningful proxy for animal-generated forces than the DA used in integrative activity measures such as ODBA (Wilson et al., 2006). Thus, these new sensing methods provide a range of options for studying animal locomotion, opening the way to an improved understanding of how animals overcome and exploit forces to move efficiently.

Two Invensense sensors (IDG500 and ISZ500) were used to construct a three-axis gyroscope. Calibration was performed in each axis in turn using a turntable rotating at 0.5–12 rad s⁻¹. Turn rate was calculated by measuring the period of the resulting cyclical signals in the magnetometer. A linear fit was made between turn rate and the gyroscope output signals to deduce scale factor and offset.

The magnetometer in the DTAG was a Honeywell HMC1043. As the magnetometer method used here is ratiometric, the absolute sensitivity of the magnetometer is not required. However, each axis needs to have the same sensitivity, requiring calibration of the sensitivity of the y- and z-axes relative to the x-axis. This was achieved using the Blainville's beaked whale dataset and a least-squares procedure to adjust the y- and z-scale factors so that the measured magnetic field vectors lie close to a spherical manifold (Gebre-Egziabher et al., 2006). The magnetometer scale factor is relatively insensitive to temperature when driven with a current source as in the tag used.

The accelerometer in the DTAG was a Kionix KXSC7-1050. As the gravity field at the Earth's surface is a nearly constant 9.81 m s⁻², this sensor was calibrated using a similar least-squares fit of measured vectors to the expected spherical manifold. This was done initially in the laboratory and then with data from the field deployment. Acceleration vectors with high minimum SA (Simon et al., 2012) were rejected iteratively to select samples with low SA for this calibration procedure.

Sensor axes were nominally aligned by the construction of the tag with a maximum expected deviation of less than 10 deg between sensors. The alignment of the accelerometer (A) and magnetometer (M) was verified by computing the dot product of the measured vectors from these sensors after calibration and comparing with the local magnetic field inclination angle [A^TM should equal gbcos(π/2−i), where g and b are the gravity and magnetic field intensities, respectively, and i is the inclination angle in radians].

Following calibration, the full-scale ranges of the sensors were calculated to be ±20 m s⁻², ±600 μT and ±9 rad s⁻¹. The root mean square (RMS) noise of the decimated sensor signals in each axis over the 10 Hz analysis bandwidth, measured in a still tag, was 0.007 m s⁻², 0.12 μT and 0.0022 rad s⁻¹ for acceleration, magnetic field and angular velocity, respectively. The power consumption of the sensors was 0.85 mW for the accelerometer, 6.7 mW for the magnetometer and 37 mW for the gyroscope.

Estimation of BR using the gyroscope was achieved with the Matlab tool gyro_rot.m (available at http://dx.doi.org/10.17630/9e38ff35-341d-442f-9f04-c63b7a48bc8b).

where the components of correspond to the BRs around the x-, y- and z-axis, respectively, and H is a high-pass filter chosen to eliminate the low-frequency drift and posture. Here, we use a symmetric FIR high-pass filter with a cut-off frequency of 0.4 times the dominant stroke frequency f_r (Martín López et al., 2015).

Estimation of BR using the magnetometer was achieved with the Matlab tool magnet_rot.m (available at http://dx.doi.org/10.17630/9e38ff35-341d-442f-9f04-c63b7a48bc8b). For a full description, see Martín López et al. (2015).

where W_t is the pseudoinverse of ⁠, i.e. ⁠.

For SA estimation, we used the Matlab tool acc_sa.m (available at http://dx.doi.org/10.17630/9e38ff35-341d-442f-9f04-c63b7a48bc8b). For a full description, see Martín López et al. (2015).

The gyroscope method gives an estimate of all three BR angles and these can be substituted into Eqn A11 to give an estimate of the orientation component of DA. The magnetometer method only gives an estimate of and so the other two angles are set to zero in Eqn A11, i.e.:

where the components of in the longitudinal, lateral and dorso-ventral axes are termed surge, sway and heave accelerations, respectively.

We thank P. T. Madsen, A. Schiavi, J. Marrero, E. Morales, K. Gikopoulou, C. Yzoard, C. Auladell, E. Burgos, C. Canut, M. Gaona, B. García, N. García, A. Peña, S. Ravagni and many others for their help in the field in El Hierro. Thanks to R. Swift and F. Alier for building the calibration platform and to K. Shorter for performing the gyroscope calibration. Thanks to D. Russell, I. Paradinas and F. Jensen for statistical and programming advice.