Sounding out life in the deep using acoustic data from ships of opportunity

The terminology used in this data descriptor follows Demer, et al.³, based mostly on Maclennan, et al.⁶⁴. All symbols signifying variables are italicized. Any symbol for a variable (x) that is not logarithmically transformed is in lower case. Any symbol for a logarithmically transformed variable, e.g. (X=1{0log }_{10}left(x/{x}_{ref}right)), with units of decibels referred to x_ref  (dB re x_ref) is capitalized.

Echosounder data

In a widely used Simrad echosounder (Table 1), the proprietary format raw data (.raw) from each transmission and reception cycle (here onwards ping) includes received echo power p_er (W), with the General Purpose Transceiver (GPT) settings: frequency f (kHz), transmit power p_et (W), pulse duration (tau ) (s), transducer on-axis gain G₀ (dB re 1), area backscattering coefficient s_a (m² m^–2) correction factor S_a corr (dB re 1), and equivalent two-way beam angle Ψ (dB re 1 sr) of the transducer. These data and associated settings were used to calculate and display volume backscattering strength ({S}_{v}) (dB re 1 m² m⁻³) for one or more frequency channels as³:

where ({P}_{er}) (dB re 1 W) is the received power, (r) (m) is the range to the target, ({alpha }_{a}) (dB m^–1) is the absorption coefficient, (lambda ) (m) is the wavelength, ({g}_{0}) (dimensionless) is the transducer on-axis gain, ({c}_{w}) (m s^–1) is the sound speed in water, (psi ) (sr) is the equivalent two-way beam angle, and the index i and j represent vertical sample number and horizontal ping number respectively.

Echosounder calibration

Echosounder calibration is a prerequisite for quantitative bioacoustic studies. The overall on-axis performance of echosounders installed on the participating platforms was routinely evaluated by established sphere calibration method^3,4. This method provides calibrated ({G}_{0}) and ({S}_{a{rm{c}}{rm{o}}{rm{r}}{rm{r}}}) required for standardizing ({S}_{v}) data (Eq. 1) collected by diverse platforms with a traceable calibration history. The sphere calibration also provides a check for transducer beam-pattern characteristics and related Ψ. The manufacturer-specified Ψ adjusting for the local sound speed variation at the calibration location was used due to the difficulty in obtaining an independent measurement of hull-mounted transducer beam pattern.

The raw data acquired using ES60 and ES70 echosounders were modulated with a triangle wave error sequence⁶⁵. The triangle wave error (with a 1 dB peak-to-peak amplitude and a 2720 ping period) was removed from calibration data before calculating ({G}_{0}) and ({S}_{a{rm{c}}{rm{o}}{rm{r}}{rm{r}}}). Open ocean transit (here onwards transect) data were not corrected for the triangle wave error due to data management and storage constraints at the start of the program. Generally, this error will average to zero over a full period of 2720 pings for normal operations and 1 km horizontal resolution of the processed data. To facilitate the processing of high-resolution data (e.g. 100 m horizontal resolution) and slow ping rate systems, transect data files were corrected for this error (if applicable) with associated metadata, since September 2020.

Data acquisition

Ensuring the operational need of participating platforms (e.g. fishing), the data acquisition settings in Table 2 were used to optimize quality and practical utility of collected data. The transmit power was selected based on the recommended²⁰ settings for commonly used Simrad echosounders. The pulse duration was chosen as a trade-off between sample resolution and acceptable signal-to-noise ratio (SNR, dB re 1) in the mesopelagic zone, and the logging range was set to provide robust estimates of echosounder background noise (dB re 1 W) levels⁶⁶.

Data registration and management

Depending on the primary purpose of participating platforms, raw data received from operators (Table 1) may cover transects and periods of fishing or scientific activities. A custom Java software suite was developed to assist data management and help identify transects for post processing (Fig. 5). These tools were used to create information (inf) files. The inf file is in plain text format that contains user-defined metadata (platform name, relevant platform call sign, voyage name, transect attributes, and relevant comments). It also includes key data acquisition settings extracted from the raw data files including frequency, transmit power, pulse duration, and echosounder details (GPT channel identifier and transducer model). The platform navigation details (total travel time, total distance covered, and average platform speed), temporal extent (start and end time of data volume), and geographic extent (limits of latitude and longitude) were also captured in the inf files. These inf files were checked for consistent data acquisition settings, transect selection, and excluding continental shelf water column backscatter data. Raw data files with inconsistent data acquisition or unknown calibration settings were not considered for further processing and archived locally.

Data processing routines

Data sets were initially processed using Echoview® software (Echoview Software Pty Ltd, Hobart, Tasmania, Australia) that includes a sequence of data processing filters⁵ designed to remove noise and improve data quality. Transect data files applying related time offset to Coordinated Universal Time (UTC) and calibration parameters were visualized (Eq. 1) as frequency-specific echograms in Echoview® for visual inspection, transducer motion correction, and filtering processes (Fig. 5). Subsequent processing and packaging were completed using MATLAB® software (MathWorks, Natick, Massachusetts, USA). All processing steps were semi-automated using a custom MATLAB® Graphical User Interface (GUI) integrated with Component Object Model (COM) objects controlling Echoview® software.

Visual inspection of data

Acoustic data quality from different platforms can vary significantly due to signal attenuation (i.e. attenuation of transmit and/or received signal to a level below the analysis threshold), and signal degradation due to combined transducer motion and noise. Data quality control involved visual inspection of echograms (Fig. 5), followed by marking the seabed (if present) and regions of bad data using echogram tools available in Echoview®. Pings with prolonged noise interference or signal attenuation were flagged as bad data. Data shallower than 10 m were removed to exclude echosounder transmit pulse and echoes in the transducer nearfield. Similarly, data deeper than the seabed (if present) were removed from the analyses. Additionally, regions of aliased seabed echoes (i.e. seabed reverberations from preceding pings coinciding with the current ping) were manually flagged as bad data. Valid high scattering from biological sources (e.g. pelagic fish schools that may occur between surface and 250 m depth) causing an apparent transition in backscatter intensities was manually preserved from the transient noise filter described below⁵.

Transducer motion correction

Echo-integration results will be biased if the change in orientation of transducer beam between the times of each ping is not accounted for. The effect of transducer motion on echo-integration was studied by Stanton⁶⁷ and later Dunford⁶⁸ developed a single correction function that can be applied for a wide range of circular transducers and related ({s}_{v}) data. To fully characterize platform movement, the Dunford⁶⁸ algorithm implemented in Echoview® requires motion data (i.e. pitch and roll of a platform) recorded at a rate above the Nyquist rate of platform’s angular motion⁶⁹ to avoid temporal aliasing due to an inadequate sampling rate. When platform motion data were available at a suitable sampling rate (see ‘Technical Validation’ section), transducer motion effects were corrected using Dunford⁶⁸ algorithm by ensuring time synchronization with recorded acoustic data (Fig. 5).

Data processing filters

Fishing vessels (FV) contributing to IMOS Bioacoustics sub-Facility were not purposely built for collecting high-quality bioacoustic data. Various factors including inclement weather and vessel design can affect data quality that could cause large biases in derived ({s}_{v}) values. To minimize these biases, data processing filters were applied to the raw data (Fig. 5). Transducer motion-corrected data were subject to a sequence of data processing filters⁵ designed to mitigate impulse noise, signal attenuation, transient noise, and background noise⁶⁶.

Data processing filters were applied to each ({S}_{v}) sample in an echogram, identified by a vertical sample number (i) and horizontal ping number (j). The ‘context window’ defined for filters include a current ping, and surrounding pings on either side of the current ping. Depending on the filter used, the context window either centres on the current ping or current sample, and slides over the entire echogram.

Impulse noise removal

Impulse noise affects discrete sections of the data with a duration of less than one ping, for example, transmit pulses originated from other unsynchronized acoustic systems. The impulse noise removal algorithm implemented in Echoview® (based on Ryan, et al.⁵) compares each ({S}_{v}) sample in a current ping to the adjacent ({S}_{v}) samples (at the same depth) in surrounding pings defined by a context window of specified width (W) (see details of context window in Table 3). A smoothed copy of original ({S}_{v}) values (i.e. unfiltered data) within the context window was used to identify impulse noise (see details of smoothing window in Table 3). The original ({S}_{v}) samples were identified as impulse noise if the corresponding smoothed ({S}_{v}) samples satisfy the condition:

where ({S}_{v}[i,j]) (dB re 1 m² m⁻³) represents smoothed copy of current ping with a vertical sample number (i) and horizontal ping number (j), (m) and (n) are the positive integer offsets from the current ping determined by the width ((W)) of context window, where (m,nin left{1,ldots ,frac{W-1}{2}right}) and (W) is an odd integer value in the range 3 to 9, and (delta ) (dB re 1 m² m⁻³) is an empirically determined impulse noise removal threshold value. Identified noise values were replaced as ‘no data’. The impulse noise removal parameters defined in Echoview® are given in Table 3.

Attenuated signal removal

Signal attenuation is generally caused by air bubbles beneath the transducer that may occur for one ping or can persist over multiple pings. The attenuated signal removal algorithm implemented in Echoview® (based on Ryan, et al.⁵) compares the percentile score of ({S}_{v}) samples in a current ping with the percentile score of ({S}_{v}) samples in surrounding pings defined by a context window (see details of context window in Table 4). The current ping was removed and replaced as ‘no data’ if:

where the symbol (p) denotes the desired percentile value, ({S}_{v}[i,j]) (dB re 1 m² m⁻³) is the current ping with a vertical sample number (i) and horizontal ping number (j), ({S}_{v}left[mtimes nright]) (dB re 1 m² m⁻³) represents ({S}_{v}) samples in the context window defined by (m) vertical samples and (n) horizontal pings, and (delta ) (dB re 1 m² m⁻³) is an empirically determined attenuated signal removal threshold value. The attenuated signal removal parameters defined in Echoview® are given in Table 4.

Transient noise removal

Transient noise is introduced to the received signal that can occur at irregular intervals and persists over multiple pings. The transient noise removal algorithm implemented in Echoview® (based on Ryan, et al.⁵) compares each ({S}_{v}) sample in a current ping with the percentile score of ({S}_{v}) samples in surrounding pings defined by a context window (see details of context window in Table 5). A smoothed copy of original ({S}_{v}) values (i.e. unfiltered data) within the context window was used to identify noise (see details of smoothing window in Table 5). The original ({S}_{v}) samples were identified as transient noise if the corresponding smoothed ({S}_{v}) samples satisfy the condition:

where the symbol (p) denotes the desired percentile value, ({S}_{v}[i,j]) (dB re 1 m² m⁻³) represents smoothed copy of current ping with a vertical sample number (i) and horizontal ping number (j), ({S}_{v}left[mtimes nright]) (dB re 1 m² m⁻³) represents smoothed copy of ({S}_{v}) samples in the context window defined by (m) vertical samples and (n) horizontal pings, and (delta ) (dB re 1 m² m⁻³) is an empirically determined transient noise removal threshold value. Identified noise values were replaced as ‘no data’. The transient noise removal parameters defined in Echoview® are given in Table 5.

Background noise removal

Background noise is introduced to the received signal that can vary in intensity and pattern (see section ‘Technical Validation’). According to De Robertis and Higginbottom⁶⁶, the calibrated ({S}_{v}) values (Eq. 1) can be expressed as the sum of contributions from the signal and noise as:

where ({S}_{{v}_{{rm{cal}}}}) (dB re 1 m² m⁻³) is the calibrated ({S}_{v}) samples derived from the raw data (i.e. Eq. 1), ({S}_{{v}_{{rm{signal}}}}) (dB re 1 m² m⁻³) is the calibrated ({S}_{v}) samples representing the contribution from signal, ({S}_{{v}_{{rm{noise}}}}) (dB re 1 m² m⁻³) is the calibrated ({S}_{v}) samples representing the contribution from noise, and the index i and j represent vertical sample number and horizontal ping number respectively.

To estimate background noise levels, calibrated received power ({P}_{e{r}_{{rm{cal}}}}[i,j]) (dB re 1 W) values were calculated from ({S}_{{v}_{{rm{cal}}}}[i,j]) values by subtracting the time-varied gain (TVG) function² (i.e. (2{0log }_{10}r+2{alpha }_{a}r)) from Eq. 1 as:

The calibrated ({P}_{e{r}_{{rm{cal}}}}[i,j]) values were averaged⁶⁶ (in linear domain) within an ‘averaging cell’ of (M) vertical samples (with an index (k)) and (N) horizontal pings (with an index (l)) to estimate noise as:

where (bar{{P}_{e{r}_{{rm{cal}}}}}[k,l]) (dB re 1 W) is the averaged ({P}_{e{r}_{{rm{cal}}}}[i,j]) values calculated for each averaging cell with a vertical sample interval (k) and horizontal ping interval (l), and (Noiseleft(lright)) (dB re 1 W) is the representative noise estimate for the ‘middle ping’ in each horizontal interval (l). Note that the averaging cell slides over the entire echogram (see details of averaging cell in Table 6).

An empirically determined maximum threshold (Nois{e}_{max}) (dB re 1 W) (see Table 6) was applied to (Noiseleft(lright)) values as an upper limit of background noise levels. Any (Noiseleft(lright)) values exceeding this threshold was replaced with the predefined (Nois{e}_{max}) value.

The (Noiseleft(lright)) value estimated for a given horizontal ping interval (l) was assigned to all individual pings constituting the interval to establish noise (Noiseleft(jright)) (dB re 1 W) estimate for each ping. The effect of TVG was added to the (Noiseleft(jright)) levels to compute ({S}_{{v}_{{rm{noise}}}}) for each vertical sample number (i) and horizontal ping number (j) as:

The background noise corrected volume backscattering strength ({S}_{{v}_{{rm{bnc}}}}[i,j]) (dB re 1 m² m⁻³) values for each vertical sample number (i) and horizontal ping number (j) were estimated as:

The SNR, a measure of the relative contribution of signal and noise was estimated as:

where (SNR[i,j]) (dB re 1) is the signal-to-noise ratio for each vertical sample number (i) and horizontal ping number (j).

An empirically determined threshold (Minimu{m}_{SNR}) (dB re 1) (see Table 6) was used as an acceptable SNR for background noise corrected ({S}_{{v}_{{rm{bnc}}}}[i,j]) data. The ({S}_{{v}_{{rm{bnc}}}}[i,j]) values with corresponding (SNR[i,j]) below this threshold were set to ‘−999’ dB re 1 m² m⁻³ (an approximation of zero in the linear domain). The background noise removal parameters defined in Echoview® are given in Table 6.

Residual noise removal

In the final stage, a 7 × 7 median filter was applied to remove residual noise retained in the core filtering stages (especially at far ranges). A median filter replaces the current ({S}_{v}) sample with the median value of ({S}_{v}) samples in a (M) × (M) neighbourhood. It is important to note that the output of 7 × 7 median filter was not directly used for echo-integration, rather it was used to flag residual noise retained from the core filtering process. A maximum data threshold of −50 dB re 1 m² m⁻³ and a time-varied threshold (TVT(r)) with the reference value of −160 dB re 1 m² m⁻³ (defined at 1 m range) was applied to the background noise corrected ({S}_{{v}_{{rm{bnc}}}}[i,j]) data before applying 7 × 7 median filter (see Table 3 for a description of time-varied threshold). ({S}_{{v}_{{rm{bnc}}}}[i,j]) values above the maximum threshold (i.e. −50 dB re 1 m² m^-3) were set to ‘−999’ dB re 1 m² m⁻³. Similarly, ({S}_{{v}_{{rm{bnc}}}}[i,j]) values below the calculated (TVT(r)) values were set to ‘−999’ dB re 1 m² m^-3 (note that median filter may replace ‘−999’ with the median of samples in the 7 × 7 neighbourhood). The output of the median filter was used to create a Boolean data range bitmap (between −998 to −20 dB re 1 m² m⁻³) with ‘true’ or ‘false’ values for each sample. This Boolean data range bitmap was applied to the background noise corrected ({S}_{{v}_{{rm{bnc}}}}[i,j]) data for removing any residual noise before echo-integration. ({S}_{{v}_{{rm{bnc}}}}[i,j]) values corresponding to ‘false’ values in the data range bitmap were set to ‘−999’ dB re 1 m² m⁻³.

Quality-controlled ({S}_{v}) data along with: (1) calibrated and motion corrected raw data, (2) transducer motion correction factor (i.e. difference between ‘motion corrected’ and ‘calibrated raw’ data), (3) background noise, and (4) SNR were exported from Echoview® as echo-integration cells (i.e. grid on an echogram) with a resolution of 1 km horizontal distance (i.e. ping-axis interval (p)) and 10 m vertical depth (i.e. range-axis interval (r)). Echo-integration values were stored as comma-separated values (CSV) files. Exported ({S}_{v}) data were converted to linear scale for further processing and packaging in MATLAB® (Fig. 5).

Secondary corrections for sound speed and absorption variation

Quality-controlled ({S}_{v}) data were echo-integrated and exported using a nominal sound speed ({c}_{w}) (m s⁻¹) and absorption coefficient ({alpha }_{a}) (dB m⁻¹) values estimated using the equations of Mackenzie⁷⁰ and Francois and Garrison⁷¹ respectively (see sound speed and absorption coefficient variables in Eq. 1 used for ({S}_{v}) calculation). However, open ocean transects pass through different hydrographical conditions, so a secondary range dependent correction was required to account for the changes in horizontal and vertical cumulative mean sound speed and absorption as:

or in linear terms:

where (bar{{s}_{{v}_{{rm{uncorr}}}}}[r,p]) (m² m⁻³) is the uncorrected (but filtered) volume backscattering coefficient values exported from Echoview® at the specified range-axis interval (r) (i.e. 10 m) and ping-axis interval (p) (i.e. 1 km), ({r}_{n}[r,p]) (m) is the regularly spaced depth values for each echo-integration cell, (bar{{c}_{w}}left[r,pright]=frac{{sum }_{r=1}^{n}{c}_{w}left[r,pright]}{n};,forall p) (m s⁻¹) is the cumulative mean sound speed values estimated at each echo-integration cell for the new range ({r}_{a}[r,p]={r}_{n}[r,p]frac{bar{{c}_{w}}[r,p]}{{c}_{w}}) (m) calculation, (bar{{alpha }_{a}}left[r,pright]=10,{{rm{log }}}_{10}left(frac{{sum }_{r=1}^{n}1{0}^{left(frac{{alpha }_{a}left[r,pright]}{10}right)}}{n}right);forall p) (dB m⁻¹) is the cumulative mean absorption coefficient values ‘interpolated’ at the new range ({r}_{a}[r,p]), and (bar{{s}_{{v}_{{rm{corr}}}}}[r,p]) (m² m⁻³) is the corrected volume backscattering coefficient values at the new range ({r}_{a}[r,p]).

Due to changes in cumulative mean sound speed, this correction step creates a grid with irregular ({r}_{a}[r,p]) values. Therefore, the (bar{{s}_{{v}_{{rm{corr}}}}}[r,p]) values at the new ranges ({r}_{a}[r,p]) were interpolated and reported at the regularly spaced ({r}_{n}[r,p]) values.

The sound speed and absorption coefficient values for secondary corrections were estimated using the equations of Mackenzie⁷⁰ and Francois and Garrison⁷¹ respectively. Francois and Garrison⁷¹ estimate their ‘total absorption equation’ to be accurate within 5% for ocean temperature values of −1.8–30 °C, frequencies of 0.4–1000 kHz, and salinity values of 30–35 PSU. The typical hydrographical conditions (temperature values of 0–27 °C and salinity values of 34–36 PSU) present along the open ocean transects are generally within the reliability limits of Francois and Garrison⁷¹ equation.

The temperature and salinity data for sound speed and absorption coefficient calculations were interpolated from either CSIRO Atlas of Regional Seas⁷² (CARS, http://www.marine.csiro.au/~dunn/cars2009/ version 2009) or Synthetic Temperature and Salinity (SynTS)⁷³ analyses (http://www.marine.csiro.au/eez_data/doc/synTS.html), but can also be derived from oceanographic reanalysis and ocean circulation models. CARS2009 is a digital climatology or atlas of seasonal ocean water properties. It is based on a comprehensive set of quality‐controlled vertical profiles of in situ ocean properties (i.e. temperature, salinity, oxygen, nitrate, silicate, and phosphate) collected between 1950 and 2008. CARS2009 NetCDF files contain a gridded mean of these ocean properties and average seasonal cycles generated from the collated observations. CARS2009 covers global oceans on a 0.5 × 0.5 degree grid spatial resolution, and are mapped onto 79 standard depth levels from the sea surface to 5500 m (from this vertical profiles of ocean properties along a bioacoustic transect can be extracted). SynTS is a daily three-dimensional (3D) temperature and salinity product generated by CSIRO, where the CARS temperature and salinity fields are adjusted with daily satellite sea surface temperature (SST) and gridded sea level anomaly (GSLA). SynTS has a 0.2 × 0.2 degree grid spatial resolution, and is mapped onto 66 standard depth levels from the sea surface to 2000 m. Due to limited spatial coverage (60°S–10°N and 90°E–180°E), the SynTS products may not always cover the transect region (e.g. Southern Indian Ocean), in that case CARS climatology values were used for the secondary corrections (Fig. 5).

Data review, packaging and submission routines

For each processed transect, secondary corrections applied (bar{{s}_{{v}_{{rm{corr}}}}}) data together with metrics of data quality and other auxiliary data variables were stored in Network Common Data Form (NetCDF, www.unidata.ucar.edu) file (NetCDF-4 format) with a resolution of 1 km horizontal distance (i.e. ping-axis interval) and 10 m vertical depth (i.e. range-axis interval) (see ‘Data Records’ section for data contents). This NetCDF file conforms standardized naming conventions and metadata content defined by the Climate and Forecast (CF)⁷⁴, IMOS⁷⁵, and International Council for the Exploration of the Sea (ICES)⁷⁶ published over the years (Fig. 6).

Processed NetCDF files were independently reviewed by both analyst and principal investigator to further investigate data quality. If suitable, the NetCDF file along with ancillary files: (1) acquired raw data (.raw files), (2) platform track in CSV format (containing date, time, latitude, longitude, and time offset to UTC), (3) platform motion data (if recorded) in CSV format (including date, time, pitch, and roll measurements), and (4) a snapshot of processed echogram as Portable Network Graphics (PNG) format were packaged and submitted to the publicly accessible AODN Portal (Fig. 5).

Source: Ecology - nature.com