A thorough understanding of the scattering characteristics of marine organisms is a prerequisite for robust quantitative fisheries acoustic data processing or interpretation. Target strength models, such as the distorted wave Born approximation (DWBA) can be used to improve the understanding of field recordings of weakly scattering targets. With acoustic methods now being used by a wide audience, allowing access to such models becomes a necessity. To ease access to the DWBA model, an r package (zooscatr) which includes a web application and the ability to parameterise the model either through the web application, text files, or pure scripting has been developed and is now freely available on Github.
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 material1 with any future developments available through Github (Gastauer, 2019).
2. The DWBA model
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 fbs 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; 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 (fnorm [ratio]), the scattering pressure (pscat), the differential backscattering cross-section (σbs [m2]), target strength (TS [dB re m2]) and the reduced target strength (RTS [ratio]) are summarised in Table 1 in the supplementary material.1
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 material1), orientation (Table 3 in the supplementary material1), material properties (Table 4 in the supplementary material1), and simulation parameters (Table 5 in the supplementary material1). 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 material1).
3. Shape generation
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 ni points along the body axis (defined within the simulation parameters), normalized x and z coordinates, tapering information (i.e., y coordinates), the positional vector , a distance vector dr1, 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 material1). Arbitrary shapes are specified by x, z position (Fig. 2) and tapering information and optionally xupper = x + tapering and xlower = x − tapering.
When an arbitrary shape is used, the coordinates for the ni 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 and , 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 ns will be set to the number of input values (nl) if it is larger than nl.
Independent of the shape, the length of the positional vector is computed as . The angle between the origin and the coordinates in radians is calculated as the two-argument arctangent of the coordinates  (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 , dx and dz, are defined as and .
4. Modelling phase
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 Cb,
is related to the compressibility 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 = cb/cs), 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 J1 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 the standard deviation of θtilt. The resulting scattering amplitude f2 will then be computed as the square-root of the conjugate product of f multiplied by the PDF of the angular distribution as and the differential backscattering cross-section becomes .
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 with dL the length increasing step, and the standard deviation of L After the length PDF has been defined, the updated σbs becomes the product of L2, 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, can then be calculated as . Following Table 1 in the supplementary material,1 the target strength or reduced target strength can then be deduced.
5. Model validation
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.1 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.
6. Conclusions and outlook
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.
See supplementary material at https://doi.org/10.1121/1.5085655 E-JASMAN-145-504901 for the zooscatr package binary (v.0.3.0). Support and package updates will be provided through Github (Gastauer, 2019). The package is distributed under the GNU Public License as published by the Free Software Foundation, either version 3 of the license, or (at the reader's option) any later version. A definition of the scattering metrics used within ZooScatR and their relationship to the scattering amplitude is given in Table 1. A description of the model shape input (Table 2), orientation (Table 3), material properties (Table 4), and simulation (Table 5) parameters are provided. The used model input parameters for a weakly scattering fluid sphere can be found in supplementary materials (Table 6).