ZooScatR—An r package for modelling the scattering properties of weak scattering targets using the distorted wave Born approximation

The use of acoustic technologies as a standard tool for the monitoring and characterisation of marine organisms requires a detailed understanding of the sound scattering signatures of organisms within marine ecosystems. With the emergence of broadband acoustics, acoustic scattering models can help us to better understand the received acoustic signal at an improved resolution, when compared to narrowband acoustics. Among the scattering marine organisms, weakly scattering (fluid-like) zooplankton or crustacean species, such as euphausiids, salps, and copepods, make up a significant portion of the global marine biomass (Strömberg et al., 2009). To accurately estimate their biomass acoustically, it is essential to provide accurate acoustic scattering models that allow for fluctuations in acoustic scattering as a function of their diverse morphological, physiological, and behavioural aspects, and acoustic frequency (Amakasu et al., 2017; Bassett et al., 2016; Stanton and Chu, 2000). The availability of fast, standardised, and flexible sound scattering models, accessible to a broader audience hence becomes desirable.

Jech et al. (2015) discussed the usefulness of the distorted wave Born approximation (DWBA) and other scattering models in fisheries acoustics for a number of targeted organisms and compared model outcomes to analytical solutions. The DWBA has widely been used to model the scattering properties of weak scattering targets, such as zooplankton and small crustaceans, e.g., krill or decapod shrimp (Chu et al., 1993; Chu and Ye, 1999; Stanton et al., 1996, 1998; Stanton et al., 2000; Stanton and Chu, 2000). Based on previously developed matlab scripts (Chu, 1999), we developed an r (R Core Team, 2018) package called zooscatr, that allows for parameterisation through .dat or .csv plain text files, a web application (Fig. 1) or r scripting. zooscatr can deal with arbitrary cylindrically symmetric elongated shapes, with an arbitrarily varying radius along the cylindrical axis. Material properties can be variable along-axis and different lengths as well as target orientations can be defined. zooscatr is open-source and open-access. zooscatr (v 0.3.0) is provided in the supplementary material¹ with any future developments available through Github (Gastauer, 2019).

A detailed description of the applied phase-compensated DWBA (PC-DWBA) model and the application to scattering properties of zooplankton can be found in Stanton et al. (1998) and Chu and Ye (1999). This variation of the DWBA model was chosen because parameterization of the model is flexible (not restrictive in terms of scattering geometry, target shape, heterogeneous material properties, and acoustic frequency), computational speed and the resulting applicability to a wide range of weakly scattering organisms. Further, the model has been validated within the context of previous studies (e.g., Amakasu et al., 2017; Chu et al., 1993; Chu et al., 2000; Chu et al., 2003; Stanton et al., 2000; Stanton and Chu, 2000; Wiebe et al., 2009). Limitations of the DWBA approach have been discussed previously (Chu and Ye, 1999; Jech et al., 2015). The DWBA model is only valid for weak scatterers, with material properties (sound speed and density contrast) within 5% variation of the surrounding medium and shear waves are assumed to be absent (Stanton et al., 2000). Different variations of the DWBA have been proposed (i.e., Calise and Skaret, 2011; Demer and Conti, 2003; Jones et al., 2009), that provide comparable results with minor variations mainly caused by differences in how the model deals with phase and shape parametrisations. The DWBA for a body of an arbitrary shape and size for a given frequency can be written as (Chu et al., 1993)

where f_bs is the scattering amplitude; k is the wavenumber (defined as 2π/λ, where λ = f/c is the acoustic wavelength [m], with frequency f [Hz] and c the sound speed [m/s] in the surrounding fluid); the subscript s refers to the surrounding fluid, the subscript b refers to the scattering body; $rv→$ is the position vector along the body axis; v is the integration volume of the body; $γκ$ and γ_ρ are related to the body material properties [Chu et al., 1993; Stanton et al., 1998; also see Eq. (4)]. The relationship between the backscattering amplitude, the normalised backscattering amplitude (f_norm [ratio]), the scattering pressure (p_scat), the differential backscattering cross-section (σ_bs [m²]), target strength (TS [dB re m²]) and the reduced target strength (RTS [ratio]) are summarised in Table 1 in the supplementary material.¹

Within zooscatr, the DWBA model can be run through the bscat() function or the web application DWBAapp() (Fig. 1), which relies on bscat(). The model takes four categories of input parameters: shape (Table 2 in the supplementary material¹), orientation (Table 3 in the supplementary material¹), material properties (Table 4 in the supplementary material¹), and simulation parameters (Table 5 in the supplementary material¹). Templates for the input parameter files are included in the package and detailed instructions on how to run the model can be found in the package vignettes.

The modelling function runs through three phases: (1) the shape generation, where the coordinates are defined; (2) the modelling phase, where all other parameters are read into the model and the scattering model is computed; and (3) the conversion phase, where the backscattering amplitude is translated into the selected output format (Table 1 in the supplementary material¹).

All shape inputs of the scattering objects are required to have circular cross sections along the length axis [shapes are formed using the buildpos (build positions) function]. Shape information contains n_i points along the body axis (defined within the simulation parameters), normalized x and z coordinates, tapering information (i.e., y coordinates), the positional vector $rpos→$⁠, a distance vector dr₁, the angle γ_t [rad] which equals the angle of a hypothetical line passing through the coordinates and the origin as well as the angular position of the points θ_tilt (Fig. 2). zooscatr can be used to either model uniformly bent cylinders or arbitrary shapes made up of a series of adjacent, potentially tapered, cylinders.

A uniformly bent cylinder with regular tapering is the default shape computed, if a shape input is not given. Construction of a uniformly bent cylinder with regular tapering requires the ratio of the radius of curvature of the longitudinal axis (ρ_c) over L (ρ_c/L), the tapering order (n) [Eq. (3)] and the ratio of L [mm] over the cylindrical radius, i.e., half width of the cylindrical cross-section (a [mm]) (L/a) (Table 2 in the supplementary material¹). Arbitrary shapes are specified by x, z position (Fig. 2) and tapering information and optionally x_upper = x + tapering and x_lower = x − tapering.

When an arbitrary shape is used, the coordinates for the n_i points are normalised with the range of x and z coordinates minus the range average. The length L of the shape is computed as the cumulative sum of the two-dimensional distance of the coordinates x and z,

The final length normalised coordinates x and z are defined as $2x/L$ and $2z/L$⁠, respectively, to fit within a range of 0 and 1. The tapering at each integration point is computed through linear interpolation of the given tapering points along the z-axis. Tapering is controlled by the tapering order (n) where the radius for a given point is defined as

Note that n = 2 corresponds to a prolate spheroid. A smoother function can be applied to both, the body axis (shape complexity) and tapering. The smoother function is a simple symmetrical moving average, where the smoother factor n_s will be set to the number of input values (n_l) if it is larger than n_l.

Independent of the shape, the length of the positional vector $rpos→$ is computed as $x2+z2$⁠. The angle between the origin and the coordinates in radians is calculated as the two-argument arctangent of the coordinates [$γt=arctan2(z,x)$] (Fig. 2). The two-argument arctangent of the coordinates z, x represents the real and imaginary components of the complex number z + ix, allowing the user to define a unique phase in four quadrants (Fig. 2). The two components of the directional vector $dr→$⁠, dx and dz, are defined as $dx=xn+1−xn$ and $dz=zn+1−zn$⁠.

Besides orientation and length, the material properties, i.e., the sound speed contrast (h) and the density contrast (g) are important model input variables. These contrasts are explicitly contained within the general DWBA equation under the form of the term C_b,

$γκ$ is related to the compressibility $κ=1/ρc2$ and $γρ$ is defined through g. Both contrasts are the ratio of the material properties inside the body and the surrounding fluid, i.e., density (g = ρ_b/ρ_s) and sound speed (h = c_b/c_s), respectively. Ideally g and h should be estimated as accurately as possible for the given target at ambient conditions (Chu and Wiebe, 2005; Matsukura et al., 2009a). Further, for euphausiids, a life history dependency as well as seasonality of g and h has been reported (Chu and Wiebe, 2005; Køgeler et al., 1987; Matsukura et al., 2009b; Mikami et al., 2000). g can be measured by using volume displacement density bottles (Foote, 1990; Greenlaw, 1979; Mukai et al., 2004), a density increment method (Becker and Warren, 2014, 2015), or through dual density measures (Chu and Wiebe, 2005). The establishment of h can be determined either through a velocimeter probe (Chu and Wiebe, 2005; Foote, 1990; Matsukura et al., 2009a) or through computational tomography (CT) scans of the targets (Au and Fay, 2012; Gastauer et al., 2017). If g and h cannot be determined directly, literature values should be consulted and model estimates should be surrounded by an appropriate uncertainty estimate.

If a homogenous body along-axis is assumed, density (g) and sound speed (h) contrasts are directly read from the parameters. For the case of an inhomogeneous body along-axis, the input g and h values are treated as means where the number of segments and the standard deviation of g and h must be defined. If the number of segments is greater than 8, the correlation length (% of body length) must be provided as well. The correlation length defines a larger variability around the midsection, with the variability decreasing towards the ends. This correlation length is the standard deviation of a Gaussian function centred at the midsection of the body.

Solving the DWBA equation with the provided parameters for an elongated fluid object with a circular cross section leads to the following equation for each orientation angle j:

With ka the wavenumber (k)∗ cylindrical radius (a), t the tapering coefficient, and J₁ the Bessel function of the first kind of order 1. Prior to transforming the model outcome into the selected output variable format, if selected, the length and angular averages will be computed. If no average should be computed, only the mean value will be considered. For angular averages a uniform or a Gaussian PDF can be calculated

With m = the number of orientation points within a user defined angle range, and $σθ=sd(θtilt)$ the standard deviation of θ_tilt. The resulting scattering amplitude f₂ will then be computed as the square-root of the conjugate product of f multiplied by the PDF of the angular distribution as $f2=f f* PDF$ and the differential backscattering cross-section becomes $σbs=f22$⁠.

For length averages either a uniform or a Gaussian PDF can be computed, where only the mean length will be considered if no length average should be computed,

where $m=6[sd(L)/dL]$ with dL the length increasing step, and $σL=sd(L)$ the standard deviation of L After the length PDF has been defined, the updated σ_bs becomes the product of L², the length PDF and the linear interpolation of σ_bs over ka. The frequency-dependent differential backscattering cross section averaged over angle of orientation and length, $⟨σbs⟩$ can then be calculated as $Σσbsi/m$⁠. Following Table 1 in the supplementary material,¹ the target strength or reduced target strength can then be deduced.

Following methods described in Jech et al. (2015) the analytical solution for a weakly scattering fluid sphere based on Anderson (1950) and for a hypothetical fluid sphere with the material properties of a generic krill (Calise et al., 2011), ensonified by a plane wave, were used to benchmark the model. A summary of the sphere parameters can be found in Table 6 in the supplementary material.¹ A detailed description of the benchmarking process is included in zooscatr as a vignette.

The overall mean error was computed following Jech et al. (2015),

The overall mean error for the zooscatr implementation was 0.31 dB, which is comparable to results reported by Jech et al. (2015) for the DWBA model (0.38 dB) [Fig. 3(a)]. While the model performed very well in general, main variations to the analytical solution occurred at the location of the nulls, where the zooscatr model generally delivered lower TS values (Fig. 3). Similarly, when comparing the model output to the analytical solution for a fluid sphere with material properties of krill [Fig. 3(c)], the model generally performed well, detecting the troughs and peaks at the correct locations [Fig. 3(d)]. The overall average error was 1.93 dB, but if only TS values > −80 dB were considered, the error was reduced to 0.76 dB. Lower values, detected within the null locations, are considered to be of low importance, as these are likely to be linked to precision and the resolution at which the model is run rather than an error in the model itself.

Comparing modelling results to analytical solutions has shown that the DWBA model holds validity for weakly scattering targets. The next logical step would be the validation of the model with in situ or ex situ information. With the availability of more and more small and flexible acoustic sampling platforms with broadband capabilities, in combination with optical recordings, the availability of fast and flexible models could help to improve our understanding of the backscattering biomass and behaviour of weak-scattering layers found throughout the world's oceans. At the moment of writing, the inversion function is a simple lookup function, minimising a cost function defined as a sum of squares error function. This is a very basic function, which should be treated with care and seen as a first step towards more sophisticated inversion methods. A benefit of the provided implementation of the DWBA model within r is its capability of running the model at relatively high computation speeds, allowing the simulation of millions of targets within minutes (on a standard notebook with six cores). Inversion accuracy could be enhanced through the application of machine learning (MLA) and deep learning algorithms (DLA). The results of large numbers (millions) of DWBA simulations will provide information gain on the variability of the acoustic scattering properties of weakly scattering targets, in respect to their shape, orientation, internal structure and the selected frequencies and could form the basis for robust MLA and DLA training datasets.

Compared to previously developed tools [such as the online DWBA calculator developed by the NOAA Fisheries Southwest Fisheries Science Center (https://swfscdata.nmfs.noaa.gov/AST/SDWBA/sdwba.html)], major benefits of the present package are that it is open-source and open-access, providing the DWBA modelling capabilities to a wider audience. In order to make the package more accessible to the users, several vignettes have been included in the package, which can be read as tutorials, describing, among other things, how to use the integrated web application or how to use it from within the r command, how to define a parameter file, and how to run the model in parallel mode for the fast simulation of numerous targets.

Future developments of the package may include more sophisticated inversion methods such as an acoustic scattering model-based nonlinear inversion algorithm (Chu et al., 2016) or MLA/DLA algorithms, the inclusion of other scattering models, porting the model processing from the CPU to the GPU would increase the speed substantially.