Fish size distribution and total fish biomass estimated by hydroacoustical methods: a statistical approach

Fisheries Research, 11 ( 1991 ) 41-73

41

Elsevier Science Publishers B.V., Amsterdam

Fish size distribution and total fish biomass estimated by hydroacoustical methods: a statistical approach B. Bjerkeng~, R. B o r g s t m m 2, ]~. Brabrand 3 and B. Faafeng l ~Norwegian Institute for Water Research, Box 69, Korsvoll, 0808 Oslo 8 (Norway) 2Department of Biology and Nature Conservation, 1432 As-NLH (Norway) 3Zoological Museum, University of Oslo, Sarsgt. 1, 0562 Oslo 5 (Norway) (Accepted for publication 20 August 1990 )

ABSTRACT Bjerkeng, B., Borgstrom, R., Brabrand, A. and Faafeng, B., 1991. Fish size distribution and total fish biomass estimated by hydroacoustical methods: a statistical approach. Fish. Res., 11: 41-73. Ten data sets of fish catches and echosounding recordings from two small Norwegian lakes were used to define the relationship between target strength (TS) and fish length for roach ( Rutilus rutilus ) and to estimate statistical errors. A maximum likelihood procedure was used, both separately for each data set, and jointly for all data, each with three alternatives regarding degrees of freedom. Results of the six different approaches have been compared with each other and to earlier estimates. Within estimated statistical errors our data confirmed an earlier estimate in the literature, T S = 20 Iog~0 L - 68, (where L = length) for a community of whitefish (Coregunus lavaretus) and smelt (Osmerus eperlanus). Two different parameter sets in the TS--L relationship were selected for use in biomass calculations. Results showed that fish biomass of pelagic areas in the two lakes varied considerably, depending on pelagic behaviour of roach. The two different parameter sets gave high estimates from 200 to 400 kg h a - ~ in Lake Gjersjoen, with standard deviation factors of 1.2 and 1.6 respectively. For Lake Arungen, high estimates were 300-500 and 900-2300 kg ha- ~. The main statistical errors seemed to be errors in the parameters of the T S - L relationship, errors due to the TS scale resolution of 2 dB, and errors in measured fish density in different TS classes.

INTRODUCTION

Recently, hydroacoustic equipment designed to operate from small boats has been developed for the assessment offish stocks in inland waters (Solli et al., 1982; Lindem, 1982; Lindem and Sandlund, 1984). Since the first measurement of fish target strength (Midtun and Hoff, 1962 ), considerable research has been carried out on marine fish stocks to establish a relationship between single fish size and target strength. However, relations between single fish size and target strength measurements are too incomplete to permit accurate sizing (Nakken and Olsen, 1977; Johannesson and Losse, 1977; Dickie 0165-7836/91/$03.50

© 1991 - - Elsevier Science Publishers B.V.

42

B. BJERKENG ET AL.

et al., 1983 ). Because target strength of fish also depends on fish behaviour and orientation (Foote, 1981 ), it is necessary to establish relations based on in situ techniques. Such information should be established empirically in natural and 'closed' systems, with only one fish species of distinct size classes present. In earlier echosounding studies (Lindem, 1982; Lindem and Sandlund, 1984; Rudstam et al., 1987 ), the relationship between target strength and fish size of smelt (Osmerus eperlanus) and ciscoes (Coregonus albula and Coregonus artedii) has been established. The last authors also showed a close relationship between depth distribution of ciscoe by echosounding and gill net catches. As portable single-transducer echosounders are practical in surveys of smaller lakes, it is important to carry out statistical analysis on target strength values and estimates on total fish biomass. We have carried out several parallel echosounding recordings and pelagic fishing in two small lakes (Fig. 1 ), in an attempt to estimate target strength-fish size relations for roach (Rutilus rutilus), which is by far the dominant fish species in many meso/ eutrophic Scandinavian lakes. The data are used to evaluate the accuracy of target strength measurements, and to calculate statistical variance of estimates of total fish biomass.

Arungen

Gjersjeen

) 0

500 m

'~~ . ,,

000 m

11"E

f Fig. I. Lake Gjersjoen and Lake ~rungen. Areas for echosounding analysis are indicated by shading.

HYDROACOUSTICAL

ESTIMATES

43

OF FISH SIZE DISTRIBUTION

MATERIAL AND METHODS

Sampling In the two investigated lakes, echosounding and captures were carried out during darkness, because at night the fish were mainly present in the epilimnetic part of the pelagic zone (Fig. 2 ). In order to identify the fish species and sizes, pelagic haul seine, floating gill nets and trawling were used.

Echosounder Echo signals from single fish were recorded along selected transects (Fig. l ) with a quantitative echosounder (SIMRAD EY-M). The working frequency of the system is 70 kHz, and the duration of the transmitted pulse is 0.6 ms, giving a depth resolution of about 0.5 m. The echo signals are shifted down from 70 kHz to 10 kHz at the calibrated signal output, to permit recording of the analog signal on a high fidelity tape recorder. All data have been recorded with a 40 log R, time-varied gain control, TVG, applied to the preamplifier. Echo signals were recorded on magnetic tape and analysed on a

Lake Arungen A E

'

i~'~,

~

"

,,,

'~/,7~

. .

,

,

.

,'! .... ~' ~',~'~..~'.

i

l

l

'

.

.

" ,

L

:,~,

i

~",~"~ , , '~' "."'~,.~.,~ ......

"

Ii . l*.i, i

,

Lake Gje~jeen 0-

41 81

121 161 E3

201 241

" "4

Fig. 2. Echogram showing single fish during darkness in Lakes Gjersj~en and /~rungen in November 1981.

44

B. BJERKENGET AL.

microcomputer by the software HADAS, according to a modified method described by Craig and Forbes ( 1969 ). The procedure tallies echo signals into 2 dB intervals, from - 38 dB down to - 56 dB. Weaker signals than this are considered to be noise, and are ignored. A more detailed description of the echosounder is given by Lindem ( 198 l, 1982).

Pelagicfish catches In both lakes, a surface haul seine of total length 50 m, depth 6 m and bag mesh size 4 m m was used in 1981, catching roach and perch (Perca fluviatills) with lengths down to 30 ram. A pelagic trawl with mesh size 5 m m in the cod end was used in Lake Gjersjoen in September 1980 and Lake/~rungen in November 1981. In September 1981 floating gill nets (mesh size: l0 mm, 19.5 mm, 22.5 m m ) were used in the pelagic zone of both lakes, covering depths from lake surface to 12 m. To catch under-yearling roach and perch in early summer, a plankton net with mesh size 0.5 m m was used in both lakes. In addition a beach seine with mesh size 5 m m was used in Lake/~rungen. The fish species present in the lakes are given in Table I. As roach was by far the most abundant or the only fish species in catches with plankton net, floating gill nets, surface haul seine and trawl in both lakes, we conclude that the echo signals reflect the pelagic part of the roach population. Echo signals from varying depth intervals between 1 and 15 m depth were analysed. In most cases an interval of 2-10 m below the transducer was analysed, as this covers the depth interval of fish sampling.

Description of the model The echosounding system is based on the fact that the mean echo intensity, or target strength (TS), from a single fish increases with fish size. The distribution of echoes over TS intervals of 2 dB will then be a measure of the distribution of fish in corresponding size classes. Target strength measurements therefore require sufficient resolution to separate single fish echoes. To correct for the dependence of echo amplitude on fish location in the TABLE 1 Some characteristics of Lakes/~rungen and Gjersjoen Lake

Area Ao (km 2)

Max. depth (m)

Altitude Fish species m a.s.i.

Arungen Gjersjoen

1.2 2.7

13.2 63.0

33 57

Roach, perch, pike, eel Roach, perch, pike, eel, bream, rudd; pike-perch, introduced in 1981

HYDROACOUSTICALESTIMATESOFFISHSIZEDISTRIBUTION

45

beam, the densities D~ of fish in n size classes are calculated from the measured number of echoes N~ by the following set of equations (Lindem, 1982 ):

Nt =Dr "At N2 =D2 "At +DI "A2

N~ =D~ "Al + D n - i "A2 + . . . + D i "An

( 1)

where At..-4n are areas given by the directivity pattern of the instrument (Lindem, 1982 ). The fish densities Di are calculated from the echosoundings as number of fish per hectare for selected depth intervals. The different size classes (in dB ) can be related to physical size (i.e. length) by using echosounding in parallel with gill netting or trawling. Available data indicate that there is a near linear relation between peak values of TS and lnL (Lindem and Sandlund, 1984; Rudstam et al., 1987). To estimate biomass from echosoundings, a relation between target strength and weight is needed. This is accomplished by combining the TS-L relation with an addition log-linear relation between length and weight, based on data available in literature and through the present study. The two equations used can be expressed as l n L = a + b - ( T S - r)

(2)

lnW=c+d- (lnL-2)

(3)

where TS = target strength of echo (dB), L = fish length (cm), W= fish weight (g), z = m e a n TS in the data set used in estimating Eqn. (2), 2 = m e a n lnL in the data set used in estimating Eqn. (3), a, b, c, d =coefficients to be estimated from the data. When the equations are written in this way, all coefficient estimates have independent statistical errors if ordinary linear regression assumptions are valid. Given a set of fish density values D (i) for target strength classes TS (i) ( i = 1.... ,n), the total biomass of fish per hectare can be expressed as

F= ~,=,

W(i,k)

(4)

i=t

where W(i,k) = biomass for a single fish in size class i as a stochastic variable. Expectation value and variance of F can be calculated for measured D (i) if the statistical properties of D (i) and W(i,k) are known.

46

B. 1KIERKENGET AL.

Statistical measures for biomass, with error estimates Error estimates for F in Eqn. (4) are based on the following error sources. ( l ) Statistical errors in estimated model parameters a, b, c, d. (2) Error due to tallying echosoundings into 2 dB intervals. Because length distributions can have quite narrow peaks, the mean length within a TS interval cannot be precisely determined: even for a large number of echoes within a TS interval, the mean length is not necessarily found near the centre of the corresponding length interval according to Eqn. (2). This error is considered to be independent of the number of echoes within the size class, and is referred to as 'TS scale resolution error' throughout this paper. ( 3 ) Stochastic variations for single fish around expectation values in Eqns. (2) and (3). (4) Statistical errors in the calculated fish densities. Statistical measures for F given by Eqn. (4) are derived as follows: first, an expression for the weight of a single fish is given, including all error sources as independent stochastic terms (/i..); second, the expression for F is rewritten, using the stochastic expression for single fish weight, and analysed as a combination of terms with different stochastic patterns (overall variation, between classes, between single fish); third, a theoretical expression for the expectation value of F is derived, expressed by the expectation values of the different terms in the stochastic expression; finally, an expression for the variance of F is found, in terms of expectation and variance values of the different terms in the stochastic expression.

A stochastic expression for weight ofsinglefish By combining Eqn. (2) and (3), and including stochastic variation terms, the weight W(i,k) for a single fish, numbered by index k, in size class i, is given by the model as

lnW( i,k )=c+d.(a+ b. [ T S ( i ) - z ] - 2 - / i 2 ( i ) + z f l n L ( i , k ) } + d l n W ( i , k ) (5) w h e r e / i t (i) = stochastic variation in mean value of In (L) for all fish in size class i, with identical distribution assumed for all i (TS scale resolution error)./ilnL (i,k),/iln W(i,k) = individual stochastic variations for each fish in lnL and In W respectively, with identical distribution assumed for all i, k. The model coefficients a, b, c, d are estimated from specific sets of data, and must be regarded as stochastic variables as well. From now on we let a, b, c, d denote expectation values of the coefficients, while the stochastic errors in the set of coefficient estimates are denoted by /ia,/ib,/ic,/id.

HYDROACOUSTICAL ESTIMATES OF FISH SIZE DISTRIBUTION

47

Throughout this paper E (x) denotes expected value (mean), of x, and V (x) denotes variance ofx. If the different stochastic terms vary independently, the expectation value oflnW(i,k) for a given i is given, identically for all k, by

E[lnW(i,k)]=E{ln[ W(i)] } = c + d. { a + b - I T S ( i ) - r 1 - 2 }

(6)

If all the different error terms are stochastically independent, only the first order terms will contribute to the variance of F. Expanding Eqn. (5) to the first order in all its error components results in the following equation

lnW(i,k) = E [ l n W ( i ) ] + d-/la+ d - S ( i ) - A b + z t c + [ a + b . S ( i ) - 2 ] .zld - d . A 2 ( i ) + d. AlnL ( i,k ) + .41nW ( i,k )

(7)

where

S(i)=TS(i)-r

(8)

A stochastic expression for biomass F The equation for F can be written

1 F=F,. i ~= l Fz(i)'F3(i)'g2(i)'k~=tF'(i'k)J

(9)

where Fl = e dn'+~c

(10)

F2(i) = e d'sti)'ab+ [a+b-S(i)--al'dd

(II )

F3(i ) = e d'aa(i)

(12)

F4 (i,k)

= e d''dlnL(i'k ) + /flnW(i,k)

( 13 )

are stochastic variables, and g2(i) = e Eu" w(i)I

(14)

is the geometric mean weight for fish in class no. i. Based on the previous assumptions about the stochastic variations, we can now state that (i) the factor F~ is common to all summed terms in F; (ii) the F2 (i) values vary from one size class to another, and have different distributions, but are not independent of each other since they are functions of common error sources; (iii) all F3(i) are independent of each other, and have identical distribution for all TS intervals; (iv) all F4(i,k) are independent of each other, and have identical distribution for all fish; (v) any F, is independent of any FK where t~ x.

48

B. BJERKENG ET AL.

Theoretical expressionfor expectation of F As Fis a sum of terms where each term consists of a product of independent stochastic variables, the expectation value of F is given by

E(F)=E(F~).E(F3)'E(F4). ~ E[F2(i) ].E[D(i) ].g2(i)

(15)

i=1

To estimate the expectation value of the different terms, we apply the general result that a function y = exp (x), where x is normally distributed, has an expectation value E (y) given by the expectation value E (x) and variance V (x) as

E(y) = E ( e x) =eE(x)+v(x)/2

(16)

Applying Eqn. ( 16 ) to Eqns. ( 10 ) - ( 13 ) results in the following expressions for expectation values of F~ to F4 E ( F I ) = e [d2-V(a) + V ( c ) ]/2 (17) E

(F2 (i)) =

E (F3)

e ~[d-S(i)]2.V(b) + [a+b.S(i) --.,I12-V (d)}/2

= e {d2"v ix(i) l }/2

E (F4) =

e {d2v [lnL(i,k) ] + v [lnW(i,k) ]}/2

(18) (19) (20)

where (i) is given by Eqn. (8) All d terms in Eqns. ( 10)-(13) have expectation equal to zero. When Eqns. (6) and ( 17 ) - ( 2 0 ) are inserted into Eqn. ( 15 ), E (F) is expressed only by expectations and variances which can be directly estimated from data.

Theoretical expressionfor variance of F The variance of F can be approximated by first order analysis. To take account of the exponential terms in the best possible way we first rewrite Eqn. (9) for F a s

F=F1 .G

(21)

where

G= ~ F2(i)'H(i)

(22)

i=1

D(i)

H(i) =F3(i)'Q(i) ~ F4(i,k) k=l

12(i) = e Eft"wt~)1= geometric mean weight for fish in class no. i.

(23)

HYDROACOUSTICAL ESTIMATES OF FISH SIZE DISTRIBUTION

49

The function FI can be regarded as a single stochastic variable independent from other terms in Eqn. (21). The functions H(i) can be regarded as a set of independent stochastic variables. The values F2 (i), on the other hand, are not independent of each other. For the expectation in Eqn. 15 this did not matter, as the expectation of a sum equals the sum of the expectations even if the terms are not independent, but for calculating the variance we have to consider the primary error factors Ab and 3d directly. The total variance V (F) of F thus can be expressed as a first order approximation using FI, Jb, 3d and H (i) as basic stochastic variables. To achieve this, F2 (i) given by Eqn. ( l 1 ) is linearized in Ab and ,Jd, thus using the approximation E (F2) = 1 instead of Eqn. ( 18 ). The result is then 2

where

The derivatives of G, all taken at expectation values over the stochastic variables Ab, dd and H(i), are given by ii

dG_ ~ {E[H(i) ].d.S(i)} ab i=l

(26)

~-~- ~1 {E [ H ( i ) ] - [ a + b - S ( i ) - 2 ] }

(27)

OG

- - OH(i)

1

(28)

where the expectation and variance of H ( i ) is given by E[H(i) ] =E(F3).E[D(i) ] .E(F4).I2(i)

(29)

V [ H ( i ) ] = ({E[D(i) ] -E(F4) }2.V(F3) + [E(F3) ]2-E[D(i) ] -V(F4) + [E(F3).E (F4) ]2-V[D(i) ] )-t2(i) 2

(30)

The expectation values of FI, F3 and F4 are already given by Eqns. (17)(20). The variances F~, F3 and F4 can be similarly expressed, by applying another general result for the distribution of a log-normally distributed function y=exp(x) V(e x) = [E(e x) ]2. (eVtx)_ 1 ) (31) with Eqn. (16) to V(F~ ), V(F3) and V(F4) in Eqns. (24) and (30).

50

B. BJERKENG ET AL.

By combining Eqns. (24)- ( 30 ), one finally gets for the variance of F 2

+

~ i=

o ( i ) - E [ D ( i ) ] - d . S ( i ,t

.V(b) 2

+ ~ {..Q(i)-E[D(i)]}z.{ [E(F3)]2_ 1} i=1

+ ~ {g2(i)2.E[D(i)]}.{[E(F4)]2 1} i=1

+ ,=,~{12(i)2.V[D(i)1})

(32)

When Eqns. (6) and (17)-(20) are inserted into Eqn. (32), V(F) is expressed only by expectations and variances which can be directly estimated from data.

Estimating the TS-length relationby peak identification By comparing density distributions over TS classes calculated from Eqn. ( 1 ) with distributions of length (L) from fish catches, one can identify peak values of TS and L in cases where the correspondence seems to be obvious, based on relative position and strength of peaks. The resulting set of value pairs can then be used to estimate L as a function ofTS. A log-linear relation between TS and length found in this way is presented by Lindem and Sandlund(1984). It might be argued that the method is subjective, and that the first peaks identified will inevitably influence subsequent identifications. However, if neither the catches nor the echosounding are significantly size selective, the method is expected to give reasonable results. If the correspondence between TS spectra and true length distributions are varying over time, between fish species or according to fish behaviour, this will contribute to the statistical error when length distributions are estimated from TS spectra. This variation cannot be expected to be included fully if the relation is based on the peak identification method, where only data with a clear correspondence between TS and length spectra will be included. The peak identification method requires clearly separated size classes and distinct

HYDROACOUSTICAL ESTIMATES OF FISH SIZE DISTRIBUTION

51

peaks. Furthermore, deviating data will probably tend to be discarded as erroneous. Thus the precision indicated by this method may be too high, because of the unavoidable subjective element.

Maximum likelihood estimation of the TS-length relation In contrast to the data of Lindem and Sandlund (1984), our data do not allow us to estimate the relation between TS and L by identifying corresponding peaks. Another method, based on maximum likelihood for the statistical distributions, is used as an alternative, in an attempt to use all available data to estimate the coefficients. The statistical errors calculated by this method will include all variability in the correspondence between TS spectra and length distribution of fish catches. This will include selectivity in the catch methods, and may exaggerate the errors between TS spectra and true length spectra. Thus, the two methods used together may indicate lower and upper bounds to the actual precision in estimating fish biomass from echosounding. Using this method, we estimate the parameters in the equation used by Lindem and Sandlund (1984), Rudstam et al. ( 1987 ) and other authors, that is: T S = A + B log(L)

(33)

The estimates will finally be inverted into the parameters of Eqn. (2): l n L = a + b ( T S - z)

(2')

At any given time there is assumed to be some unknown 'true' population of fish in the lake, with both the fish catches and the echosounding being statistical samples from this population. The fish catch data are given as semicontinuous distributions, for our data mostly in 1 mm body length intervals, and the echo distributions in 2 dB intervals. For any given set of K length intervals there will be a corresponding set of f values, with f = relative frequency of fish in interval number i, i = 1,...,K. It is assumed that a statistical sample of N fish is distributed among the K length groups following the multinomial distribution P ( n l ,n2 .... ,nK) m n l

N! f n If n 2 ['nK [n2[...nK[Ji J2 ...J K

(34)

where

N=nl +n2...+nx

(35)

If a set of values A, B are assumed for the coefficients of Eqn. (33), and the

52

a. BJERKENG ET AL.

length intervals are made to coincide with the 2 dB TS intervals, then the joint probability for the distribution of catch and echo data in these intervals is given by P=pc(nl,n2 ..... nx)'pe(ml,rn2, .... mr)

(36)

where the first factor (Pc) represents the catch sample, and the second factor (Pc) is the independent probability of the echosounding distribution, both given by formulas like Eqns. (34) and (35). The total number of echo signals is given by

M=ml +m2...+mr

(37)

The two probabilities are connected by using the same f values for both distributions. The maximum likelihood procedure ordinarily aims to find values for the parameters to be estimated so that a probability density function is maximized. These values are then considered the most likely values, giving the closest fit between the two distributions. At first it seems natural to apply this idea to the given problem by maximizing the joint probability P. The parameters to be varied are the unknown f values and the coefficients of the TS-L relation. It is easy to show analytically that for fixed values of n~and rn~,P is maximized by putting +m~

f -nN+ M

(38)

The task then would be reduced to finding the coefficients A, B which maximize P given by Eqns. ( 34 )- ( 38 ) as described above. Distribution among TS intervals of the echosoundings is fixed, and thus K, M and the rn~ values are constants. The n, values will vary as functions of the coefficients A and B, giving varying values of pe and Pc. This could work if the TS intervals covered the whole range of fish sizes in the distributions. In that ease we could restrict the variation in the coefficients to the eases where all fish caught were included in N. However, echo signals are only recorded in the range from T S = - 3 8 to T S = - 56 dB, and both big and small fish might give echoes outside this range. If this is allowed for, the probability Pc will tend to increase as a decreasing part of the catch is aligned with the echosounding data. To counteract this bias towards placing the catch outside the TS range, the probability P is multiplied by the density of sample points in the probability space, to give an approximate probability distribution function (pdf). This pdf is then maximized. We want to compare the density functions for different N' values, where N' is the total number of fish caught in the length intervals coinciding with the TS range. The observed relative frequencies nffN' in

HYDROACOUSTICALESTIMATESOF FISH SIZE DISTRIBUTION

53

each TS class, i = 1,...,k, are then used as stochastic variables going from 0 to 1 for any N'. In the resulting K-dimensional probability space, which has a constant 'volume' independent of N', the density of so-called sample points is then proportional to the number of such points, given by

V=[N'+(K-1)]/(K-1)

(39)

This can be easily verified by considering how many distinctly different ways there are to arrange N identical objects and K - l identical 'class separators' in a sequence. The procedure used maximizes the product pdf=V.P

(40)

(or rather log (pdf), which is equivalent). Because of the correction made by Eqn. ( 1 ), mi is calculated as the total number of echoes multiplied by relative fish density in each TS interval. It is also possible to extend the method to include variation in echo signals for fish of a given length. This was in fact done, as described below. The method described here has been programmed in FORTRAN, using a variant of the Nelder-Mead simplex optimization procedure (Nelder and Mead, 1965 ). The starting points for the optimization are found by selecting the best from a set of initial estimates of A and B based on mean and variance of the TS and length distribution, supplemented by a number of random parameter value sets, different for each run of the program. The initial selection is biased towards using parameter sets where most of the fish caught are included in the TS range from - 5 6 to - 3 8 dB, but then N' is allowed to be reduced unrestricted, if this leads to higher pdf values. The program is run repeatedly to get an indication of the stability of the results with respect to local maxima. The program was tested with artificial distribution pairs where a perfect fit existed for a given set of parameter values. The test showed that the method actually found these values. Three possible ways of defining the relation TS = A logL + B are applied. ( 1 ) The regression line is forced through the point TS = - 65.0 dB, L = 3.58 cm. This point is defined by the mean values of data for 7 small fish in Foote et al. ( 1987 ), where TS is calculated theoretically from the measured swim bladder size. Only A is then estimated, and B given as a function of A. (2) The regression line is not restricted, except that B is kept positive. Estimated parameters are A and B. Variation in TS for given length is not allowed for. ( 3 ) The same regression line as in (2) is used, but the echo strength of fish with a given length is assumed to be distributed normally, with an unknown standard deviation tr, estimated as a parameter in addition to A and B. The

54

B. BJERKENG ET AL.

optimization here proceeds in two steps, first keeping tr= O, and then allowing all three parameters to vary from this intermediate optimum. For each of these alternatives the maximum likelihood method is used in two ways. (i) Estimates are made for each of ten available data sets separately. The mean and standard deviation of these estimates are calculated with equal weight on each estimate, by robust 2-winzorized estimators, described for instance in Afifi and Azen ( 1979 ). The reason for using robust estimators is discussed later under 'Results'. A g-winzorized mean of a sample of n values is found by sorting the values, replacing g extreme values at each end of the sorted sequence with value no. g + 1 and n - g , respectively, and calculating the mean of this modified sample. The standard deviation estimate for the modified sample, multiplied with Q = ( n - 1 ) / ( n - 2 g - 1 ), is used as an approximate estimate for the standard deviation of the original distribution, in accordance with the confidence interval given by Afifi and Azen (1979). The standard deviation is divided by n to give the standard error of the mean estimate. The covariance matrix for the mean estimate of A, B and a is also calculated, using the same robust method, multiplying the covariance matrix of the modified sample by the factor Q squared, and dividing by n. (ii) A joint estimate is made for all data sets simultaneously, maximizing the product of the pdf's, with different fi values for each set, but with the same parameter values A, B. This means weighting by numbers, according to the theoretical multinomial distribution. Obviously there is no 'correct' way of weighting the data sets. The multinomial model is only an ideal model, and there can be systematic, large-scale differences between catch and echosounding sampling, more or less independently of sample sizes. On the other hand we will expect large samples to be better than small samples in eliminating stochastic errors. Alternatives (i) and (ii) above represent two extreme assumptions regarding weighting, and the differences in results for the two variants should indicate the possible error due to undefined weighting.

Statistical error in fish density calculations To get an estimate of the precision of the fish density measurements, we have used repeated echosoundings available from Lake Konnevesi and Lake Mjosa (Lindem, 1982 and personal communication, 1985 ). The values are given in Table 2. The total water volume 'seen' by the instrument varies with depth interval and the duration of the run; the ratio between number of echoes and estimated fish density varies accordingly. Measurements taken within an hour and during the same light conditions

HYDROACOUSTICALESTIMATESOF FISH SIZE DISTRIBUTION

55

TABLE 2 Repeated echosounding from Lakes Konnevesi and Mjosa. Data from Lindem ( 1982 ) Description

Depth interval (m )

Run A Echoes NA

Run B Fish h a - ~ DA

Echoes NB

Fish haDa

Konnevesi, Sept. 1980 Same hour, A: 91 pulses min -~ B: 182 pulses m i n - t Same hour, A: 91 pulses min -~ B: 182 pulses min -~ A: Sept. 7 B: Sept. 9 Different length of runs

1-10 10-20

19 701

275 1200

50 1530

262 1304

4-8 8-16 16-23 5-9 9-16 16-20

Il 97 1381 23 194 370

118 396 1472 369 986 597

28 313 2787 172 490 I 122

140 449 1484 642 905 756

Mjesa, Nov. 14 1980, within l hour

10-80

13189

4283

12185

4535

would be expected chiefly to reflect the statistical variation due to short term, more or less individual and random fish movements, while measurements separated by one or two days might include more of the variation due to correlated fish movements. Table 2 includes data sets of both types. The available data are not sufficient to describe such differences. In an attempt to use the available data to establish an empirical model of the statistical error as a function of echo number and fish density, we make the following assumptions. As a reasonable variance model, we have chosen the expression V ( D ) = [E(D)]2 ( ~ - + k o )

(41)

where D = estimated fish density (number per hectare), N = e c h o number, E ( D ) =expected value, or true mean of D, and ko, kl are unknown constants to be estimated. The variance consists of a sum of two terms, each proportional to the square of the fish density. The first term is the purely stochastic error inherent in the method, inversely proportional to number of echoes. With k~ = l this would be the variance of the Poisson distribution. The second term is independent of number of echoes, and can include both biological/ecological variance due to correlated behaviour, and other observational errors. For two parallel measurements, denoted run A and run B, of the same density D, the expected square of the difference between the results is given by E[ (DA --DB) 2 ] = V ( D ) A + V ( D ) B

(42)

56

B.BJERKENGETAL.

with V(O)~ = [E(D) ]2. [ (k~/N~) +ko];tz=A,B

(43)

If we assume that the differences between A and B in all parallel measurements are approximately normally distributed, then the variable (DA --DB) 2

(DA --DB) 2

U=E[ (DA --DB) 2] -- V(D)A + V(D)B

(44)

has an approximately chi-square distribution with 1 degree of freedom. The most precise estimate of D is found by weighting proportionally with the inverse of the variance /5_{ [NA/(k, + koNA) ]'DA} + { [NB/(k, +koNB) ]'DB} [NA/(kt +koNA) ] + [NB/(k, +koNB) ]

(45)

If this estimate is inserted in the expressions for V(D)A and V(D)B, we get quite precise estimates "~(D)A. and 9(D)B for these values if the variance model is valid, and the relative error of D is small. Thus we assume that even the estimates U'-

(DA --DB)2 V(D)A +V(D)B

(46)

have chi-square distributions. If we calculate 0 for each of the parallel measurements in Table 2, we get with o= 9 cases: E(i~l O~)=O= 9

(47a)

E(/~l/~/)= (0+2)0=99

(47b)

The coefficients kl and ko are found by iterating until both conditions are reasonably well fulfilled, minimizing the sum of the squared deviations from expected values of the sums, relative to theoretical variance according to the chi-square distribution. The data and/or the model is not sufficient to give reliable estimates of both coefficients simultaneously. By estimating one coefficient at a time, setting the other equal to zero, we get approximate upper limits for the coefficients ko =0.013 kt =3.6

HYDROACOUST1CALESTIMATESOF FISHSIZEDISTRIBUTION

57

The analysis above has been made for the accumulated density over all size classes. To apply it to each class, Eqn. ( 1 ) should be taken into account. For simplicity we use for size class i the relation V [ D ( i ) ] = {E[D(i)]}2.{

[k~/N(i) ] +ko}

(48)

where N ( i ) is calculated as

N(i) =Ntot .D(i)/~ [D(i) ]

(49)

I

This variance estimate should include variations due both to the method and to short-term biological movements. RESULTS A N D DISCUSSION

Estimating the TS-length relation by peak identification Table 3 gives data for separate peak values of TS and L resulting from parallel echosoundings and trawl samplings in Lake Mj~sa, where the dominant species are smelt, ciscoe and whitefish (Lindem and Sandlund, 1984). Each peak in L will in general represent one age class for small fish. As mentioned, the data indicate a linear relation between TS and log(L), expressed by Eqn. (33). The coefficients in this equation, with statistical errors, can be estimated by linear regression on the data in Table 3, and the results are: A = - 6 8 , B = 2 0 , var(A)=1.152, v a t ( B ) = 1.022, and cov(A,B) = - I. 136. The covariance of A and B expresses the fact that the data are centred on a mean value of log L different from zero; the relation is cov (A,B) = - var (B). (log L ) me,n If instead the data in Table 3 are used to estimate the coefficients in Eqn. (2), which is the inverse of Eqn. (33), the result is: a=2.51, V ( a ) =0.032, = - 46.5, b = 0.112, V ( b ) = 0 . 0 0 6 2 , and V ( residual lnL ) = 0.142. V ( residual lnL) can in this case be taken as a combination of TS scale resolution error and individual variation. The residuals vary quite systematically with predicted lnL. A plot of TS against log~oL (from Lindem and Sandlund, 1984 ) shows that most of the data points are found in two parallel nonlinear bands separated by 2 dB on the TS scale (Fig. 3 ). The two bands represent to a large extent two separate years. They are possibly an effect of the tallying of both TS spectra and length distributions into discrete intervals in the data used to extract the values in Table 3. As the residuals include systematic components, confidence intervals cannot be calculated from t distributions as usual for this regression. Echo target strength from fish of a given size will not only depend on fish size, but will also vary with the orientation of the fish relative to the direction

58

B. BJERKENG ET AL.

TABLE 3 Peaks from echosounding and net/trawl catches in Lake Mjosa, from diagram in Lindem and Sandlund (1984) TS n

L2 (cm)

TS

(dB)

L (cm)

TS (dB)

L (cm)

TS

(dB)

(dB)

L (cm)

-56 -54 -54 -52 -52

4 4 6 6 8

-50 -50 -48 --46 -48

8 10 10 ll 13

-46 -44 -44 -42 -42

15 17 17 20 21

-40 -40 -38 -38

21 28 31 32

Target strength. 2Fish length. _ TS = 19.72 Log L - 68.08 i~' -40- n = 19 ir~/'l r 2 = 0.96 //~-,4~

?/

-so .I

7f0~6 ' 018' 110 ' 112 ' 114 ' 1 . 6 log Length

,4 ? 8, lo 1,s 20 2,s, as Length (cm)

Fig. 3. Target strength and fish length relationship given by Lindem and Sandlund ( 1984 ) by linear regression (whole line) and in two non-linear bands (broken lines), largely representing two separate years.

of the beam. This means that a fish of a given size may give echoes over a range of TS values. This results in a statistical error when calculating In (L) from TS, an error which is difficult to estimate. Frequency distributions of length and TS (Lindem, 1982) indicate that the width of a TS peak is not much larger than given directly when the width of the length peak is converted to the TS scale by Eqn. (33). This would mean that the variation in signal strength for a fish of given size is not statistically important, and that V(lnL(i,k) ) can be neglected compared with the TS scale resolution error. On the other hand Nakken and Olsen ( 1977 ) found that considerable scatter can be expected due to variations in tilt angle and tail beats. This could mean that the variance could be of importance. The echo distribution is measured only for signals down to a strength of 56 dB. Weaker echo signals are considered to be noise and ignored. Thus, small fish are not included in the calculations. The SIMRAD EY-M echosounder operates with a frequency of 70 kHz, which makes it difficult to dis-

HYDROACOUSTICAL ESTIMATES OF FISH SIZE DISTRIBUTION

59

tinguish between different size classes for fish below 5 cm. This will reduce the advantage of the method when applied to populations dominated by small fish. Also size selectivity of the catches underestimates small fish. In both lakes, the O + fish in the length range 15-30 m m were not caught at all by haul seine (Fig. 4a, b).

Maximum likelihood estimation of TS-length relation In our data from Lakes Gjersjoen and Arungen there are in most cases no distinct TS peaks, but rather more or less continuous spectra over several dB intervals (Fig. 4a, b ). By comparison, the data in Table 3 from Lindem and Sandlund (1984) seem to be based on more well-defined peaks, separated by intervals of zero frequency, and a better visual correspondence between length distribution and TS spectra. One possible explanation could be that Lindem and Sandlund's data are based on observations of fish species with more distinct peaks in the length distribution. Compared with young stages of smelt, ciscoe and whitefish, the first year classes of roach have a small annual length increment, resulting in a small distance between the different peaks in the length distribution. Beyond l 0-12 cm, the length frequency distributions show no distinct peaks, due to the overlap in lengths between different year classes. Thus, for our data one cannot identify corresponding length peaks directly, and the maximum likelihood method as described above has been tried as an alternative method. Ten pairs of data sets (fish catches + echosounding recordings) were selected, for which it was believed that the catch gave a reasonably complete picture of the pelagic fish population. Table 4 lists information about the data sets used, and gives results from the maximum likelihood estimation of the parameters in Eqn. (33) made separately for each data set. The results were found by running the likelihood maximization procedure repeatedly, with randomly varying starting points, as described under Methods. When only A were allowed to vary freely, two repeated runs showed consistent results between runs. For the two and three parameter estimates, six repeated runs were made, with the two parameter set (A,B) in each case being the starting point of the three parameter set (A,B,a). The repetitions showed that the same data set could give convergence to different local maxima due to different starting conditions. The calculation procedures were for the most part successful. Typically more than 90% of the catch was included in the length intervals coinciding with the TS intervals. Occasionally likelihood maxima were found for parameter values with a very low fraction (1-5%) of the catch within the TS range; such values were accepted with the other results. The values in Table 4 represent for each data set the result with the highest maximum likelihood, chosen over the series of runs.

60

B. BJERKENGETAL.

(A) Lake Gjersj~n 10-

Lake Gjen~jeen: July 1981

60-

Lake Gjerlieen: July 1gel N = 1227 fith/ha

>

4( "='4u-2( 2-

D

Fish h m ~ h in mm

~

'48' '46' '44' '42' " ~

Target strength in +dB

Lake Gjenljeen: August 1981

60-

Lake Gjer~jeen: AuguSt 1981 N = 1171 fish/ha

~4( > u

=, u. 2(

,,

'1~'1do'~o Fi|h length in mm

L

'52' '50' '48' '46' Target strength in +00

'44'

'42'

'40'

Lake Gjeftj(Nm: November 1981

60-

Lake Gjersieen: November 1981 N = 20182 fish/ha

_ 4C

~2C

Fish length in mm

Target strength in +dB

Fig. 4. Length distribution of roach from pelagic areas (haul seine) and distribution of target strength of echo signals in Lake Gjersjoen (A) and Arungen (B) at three sampling times in 1981.

6]

HYDROACOUSTICAL ESTIMATES OF FISH SIZE DISTRIBUTION

(a) Lake Arungen 1t

Lake Arungen: July 1981

-[ 6"I~-"

Lake ~,rungen: July 1981 N = 5~o46 fish,ha

> 40

40

60

IC-

80 100 120 140 160 Fish length in mm

v,

190 200

'56'

'54'

'52' '50' '48' '46' '44' Target strength in +dB

'42'

'40'

Lake Arungen: A u ~ s t 1981

86O-

Lake Arungen: Augult 1981 N= 1712 fish/ha

40 ~4

g ~

40

10-

60

80 100 120 140 Fish length in mm

160 180 200

2(I

0

'56' Target strength in +d8

Lake Arungen: November 1981

860 4 I I

Lake Arungen: November 1981 N = 11141 fish/ha

i

i 4-

4°I

2-

I ' ,m

'eb' ~b' 16o' i~o'i~.o'i~o' 16o'26o Fish length in mm

! o~

'56'

'54'

'82' '5O' '48' '46' '44' Target strength in +dB

'42'

'40'

62

B. BJERKENG ET AL.

TABLE 4

Data sets used in maximum likelihood calculations of parameters in the relation TS = A + B*log~oL, with results of calculation ~for each data set separately, under three alternative assumptions described in the text Lake

Month Year Single fish echoes

Catch size (no. offish)

Number of estimated parameters 1

2

A

A

3 B

Gjersjoen

May Sep. Jul. Aug. Nov.

1980 884 1980 2217 1981 201 1981 231 1981 5418

397 220 687 498 76

-90.2 -97.3* -84.8* -82.9* -88.3

- 7 1 . 8 44.3* - 7 9 . 6 35.4w - 1 1 1 . 0 64.6* - 6 4 . 2 12.8" -55.9 6.2*

Arungen

May Jul. Aug. Nov. Sep.

1981 635 1981 4986 1981 221 1981 14342 1982 2606

809 1290 1290 8560 1600

-106.5" -93.9 -88.5 -87.1w -94.5w

- 7 9 . 9 24.5 - 8 4 . 0 23.7 - 7 2 . 3 23.5w - 8 2 . 3 34.9 - 7 2 . 8 26.8

A

B

a

- 7 1 . 8 44.3* - 7 4 . 0 30.5w - 6 3 . 2 ll.6w - 7 2 . 7 21.0 - 7 2 . 8 26.8

0.0 1.7 1.3 1.3 1.8

-56.2 -55.7 -72.3 -81.5 -72.8

2.5 4.5 0.0 0.1 0.0

3.2* 1.8" 23.5 34.1 26.8

~Calculation of robust winzorized means is indicated as follows: *, outliers in 2-winzorized means; w, replacement values (nos. 3 and 7 in sorted set) based on b values for alternatives 2 and 3. See methods for explanation of winzorization. 2Combination of two echosounding transects.

Inspection of Table 4 shows that values differ considerably between data sets. Some of the parameter sets in Table 4 deviate strongly from the rest, and there are some unlikely values. For instance, the three parameter estimates for Arungen in May and July 1981 give B values in the range 2-3, indicating the unreasonable result that echo strength barely depends on fish size et all. Instead, the TS spectra in these cases are adapted to the catch length distribution by a high value of tr. The repeated runs showed that results can be unstable, and that local maxima were sometimes found far from the best values, depending on initial guesses. Similar differences to those seen in Table 4 between data sets were encountered between runs for the same data set. The procedure would probably work better for data with somewhat more pronounced peaks. This is the reason for relying on robust 2-winzorized estimators as described under Methods. The extreme values are assumed not to represent real variations in T S - L relations, but rather to be a result of the mentioned instability. Both A and B were winzorized based on the values of B, that is, the two parameter sets with the most extreme B values in each end of the sorted sample were replaced with parameter sets nos. 3 and 7 in the sorted sequence, respectively.

HYDROACOUSTICAL ESTIMATES OF FISH SIZE DISTRIBUTION

63

Table 5 summarizes the statistics for the different estimation procedures, and includes joint estimates based on adapting all data sets simultaneously. The joint estimates are chosen from repeated runs in the same way as the single data set estimates in Table 4. When the TS-length relation is forced through the point given by Foote (1987), we find values A, B in the approximate range ( - 9 0 , 45) to ( - 8 5 , 40) (Alternative 1 in Table 5). This is not significantly different from the uncertain estimate A = - 98, B = 60, with standard error SE (B) = 17, which can be made from the data in Foote et al. ( 1987 ), when the statistical variation in measured swim bladder size is accounted for. The results of Alternative 1 are, however, significantly different from the values reported by Lindem and Sandlund ( 1984 ) for smelt, ciscoe and whitefish. When both A and B are allowed to vary, but without allowing for varying signal strength for given length (Alternative 2 in Table 5 ), the results indicate TABLE 5 Results o f m a x i m u m likelihood calculations o f the T S - L relation: I, separate data set estimates (Table 4 ); II, joint estimates based on adapting all data sets simultaneously Alternative 1: function TS = A + B log ( L ) forced through the point TS = - 65, L = 3.58 (from Foote et al., 1987) B -65-A

1. II.

Mean, 2-winzorized Std error o f mean ( w i n z o r i z e d ) Estimated values

A

= ~

- 90.5 1.9 - 86.4

46.1 3.4 38.6

Alternative 2: function TS = A + B log ( L ) , A unrestricted, B > 0 A Mean, 2-winzorized Std error o f mean ( w i n z o r i z e d ) Covariance estimates for mean values ( w i n z o r i z e d )

A

-77.5 2.6

28.7 3.3

6.7

-4.3 10.9

- 74.1

25.5

B

II.

Estimated values

Alternative 3: function TS = A + B log(L ), A unrestricted, B > 0, with o in TS for given length A B o I.

1I.

Mean, 2-winzorized Std error o f m e a n ( w i n z o r i z e d ) Covariance estimates for mean values ( w i n z o r i z e d )

Estimated values

A B a

- 70.2 2.8

22.4 4.6

7.7

-- 12.3 21.2

-71.0

22.8

1.2 0.4 0.008 0.173 0.145 2.4

64

B. BJERKENG ET AL.

a less steep line T S = f ( I n L ) , with A, B = ( - 7 6 , 27) as a reasonable overall mean. This is closer to the values of Lindem and Sandlund ( 1984 ), although still significantly different at the 5% level, judged from the estimated standard errors, which are approximately _+3 for both A and B. The robust mean of the separate data sets and the joint estimate for all data sets are consistent within estimated errors. When the variation in signal strength for a given length is finally taken into account through the third parameter a (Alternative 3 in Table 5 ), estimates are A, B = ( - 7 0 . 2 , 22.4) and ( - 7 l, 22.8 ), respectively, for robust mean of separate estimates and joint estimate over all data. This is not significantly different from the values of Lindem and Sandlund (1984), considering the estimated standard errors, which are larger than given by Lindem and Sandlund's data. Figure 5 shows the relation between TS and log L for Lindem and Sandlund's data, and for the m a x i m u m likelihood Alternatives 2 (I) and 3 (I). In this figure, the confidence intervals for the regression lines are given by the formula TS=A+B.IogL_+D where the half width of the confidence interval D is given by D = t ( p ) ~ . [var(A*) + ( l o g L - A ) 2 - v a r ( B ) Z ] i/2 with A= -cov(A,B)/var(B) A*=A+B-A For Fig. 5a, the line given by Lindem and Sandlund, this is the normal formula for the confidence interval of the regression of TS on logL, as given for instance by Afifi and Azen (1979). A is the sample mean of logL, and cov (A,B) has been derived as - A * v a r ( B ) . The adjusted intercept A* is the regression value of TS at the mean logL value .4, and by definition cov (A*,B) = 0. Figure 5b and c does not show ordinary regression lines. Here, estimates of var(A), var(B) and cov (A,B) are given by statistical analysis of the maxim u m likelihood coefficients for separate data sets, as given in Table 5. The values of.4 are derived by the formula above, ensuring that cov (A*,B) = 0. Thus, the confidence intervals are calculated by analogy with linear regression. Fig. 5. Target strength-fish length (L in cm) relationship for data based on (A) Lindem and Sandlund ( 1984): TS = 20. IogloL- 68, (B) Alternative2 (I ) this study: TS = 28.7.logjoL- 77.5, and (C) Alternative 3(I): TS=22.4-1og~oL- 70.2. Cases (B) and (C) are both estimated by maximum likelihood method. 95% confidencelimits are indicated. See text for explanation.

HYDROACOUST1CAL ESTIMATES OF FISH SIZE DISTRIBUTION

! A

..S~'./

-40-

/

/,

"0

O')

g -5c ~9

.,$

-60i

0.4 ' 016 ' 0.8

-4o1 -40

i

i

1.0

i

,

1.2

i

i

1'.4 1.'6

i

1.8

/.~

B

/

-5o

////,

-

.,o-////

0.4' oZ6 ' 0;8 ' 110 ' 1;2 ' 1'.4' 1;6

1.8

c

.

-40-

IDl

/

-50

-60-

j"

/

. . .0',8 0.4 . 016

4 I

1'.0 1 12 ' 1'.4 ' 1~6 log Length 6 8 10 15 2025 35 Length (cm) I

I

I

I

,

,

t

.8

65

66

B. BJERKENG ET AL.

The two-tailed t ( v ) . = 0 . 0 5 test statistic has been set equal to 2. l for Fig. 5a ( v = 1 7 ) and 2.25 for Fig. 5 b a n d c ( v = 9 ) .

Choosing between alternative estimates for the TS-L relation The choice of coefficients for Eqn. (2) is made as follows. Since the TS-L relation in Alternative 1 in Table 5 is forced through the mean point of the small data set given by Foote et al. ( 1987 ), and we have not considered the error in these data properly, we choose not to use this alternative in the calculations. Alternative 2 gives reasonable results, although with somewhat higher B values than Lindem and Sandlund's (1984) estimates. It is generally assumed that there is real variation in TS for fish of a given length. The effect of this should be to give broader TS spectra for a fish population with a given length distribution than given directly by a mapping through Eqn. (33). For Alternative 2, this should tend to bias B towards increased values when the maximum likelihood procedure is used, since this procedure will try to fit the two distributions to each other using such a mapping. Alternative 3, which allows explicitly for this variation, should on the other hand give unbiased values of B. The problem with Alternative 3 could be that it introduces too many degrees of freedom. For TS and L distributions with few specific features in the form of peaks the optimization problem may become too undetermined, as the results in Table 4 indicate. By using robust estimators, however, we have arrived at reasonable results, and the joint estimate (Alternative 3 (II) ) seems to be quite consistent between runs. In light of this, we have chosen the parameters estimated by Alternative 3, which are also closest to Lindem and Sandlund's estimates. As shown by Fig. 5, Alternative 3 also gives much more precise estimates in the most important TS ranges, compared with Alternative 2. To see the importance of the errors in estimating Eqn. (2) for the biomass calculations, we will use two alternative sets of values for the parameters a, b with variances V ( a ), V (b) and r: P 1. Parameter values estimated from Lindem and Sandlund (1984) by peak identification. P2. The coefficients for Alternative 3, inverted into coefficients of Eqn. (2). We have chosen to rely on the joint estimate for three parameter optimization (Alternative 3(1I )), but the separate estimates for each data set are used to calculate the r value, giving approximately independent errors in a and b as required by the derived formulae for E ( F ) and V ( F ) . The inversion is done as follows: all the parameter sets A, B from Table 4, Alternative 3, are inverted into values a*, b in the equation lnL = a* + b-TS, using the following formulas (note the transition from log~o to In )

HYDROACOUSTICAL ESTIMATES OF FISH SIZE DISTRIBUTION

67

a*= - l n ( 1 0 ) - A / B

(50)

b=ln(10)/B

(51)

The sample of inverted single data set parameters a*, b is then 2-winzorized according to the b values, and the value of z estimated by r = - c o v (a*,b)/var (b)

(52)

using the covariance matrix of the winzorized sample. Finally the inverted parameter set a*, b for the joint estimate over all data (Alternative 3(1I) ) is converted into the parameters a, b of Eqn. (2) for the calculated r, using the formula a=a*+b.r

(53)

As the robust mean and the joint estimates of Alternative 3 are so close together, the inverted parameters a and b can now be expected to be approximately stochastically independent for the winzorized sample of separate data set estimates. For both alternatives, P 1 and P2, an upper limit for the TS scale resolution error can be estimated from a rectangular distribution over a 2 dB TS interval, and is given by V[2(i) ] =b2Vl (TS)

(54)

where V1 ( T S ) = 0 . 6 7 = m e a n square TS scale resolution error over a 2 dB interval. The estimate of V [lnL(i,k)], the individual variance of lnL for given TS strength, can be calculated as

V[lnL(i,k) ] =b2V2 (TS)

(55)

where V2(TS) is the variance of individual TS signal strength for fish of a given length. The maximum likelihood method indicates that V2 (TS) is found in the range from 1.22 (Alternative 3(I) ) to 2.42 (Alternative 3(II) ). We will use the larger estimate for both Lindem and Sandlund's data and the maxim u m likelihood alternatives. The final estimates for the parameters in Eqn. (2) which are used in the calculations are given in Table 6, and in Fig. 5. Both alternatives give about the same value oflnL as a function of TS in the range - 56 to - 4 8 , that is for the length range 4 to 10 cm. For larger fish there is an increasing difference between the two estimates, with TS = - 38 corresponding to 32 and 38 cm for the P 1 and P2 estimates, respectively.

Estimating the length~weight relation To calculate total fish biomass, we used data on length and weight for 101 roach caught in Gjersjoen, covering the length range from 5 to 10 cm, which

68

a. BJERKENG ET AL.

TABLE 6 Alternative estimates for the coefficients in Eqn. (2), with statistical precision measures Parameter

Lindem and Sandlund (1984)

Max. likelihood, Alternative 3

a b z V(a) V(b) V(lnL(i,k) ) V(2(i) )

2.514 0.112 -46.5 0.032 0.0062 0.07 0.008

1.57 0.12 -55.6 0.0872 0.0312 0.08 0.0096

is often most abundant in the echosoundings, according to the TS-L relations. These data fit reasonably well with the data on mean length and weight for roach in different age classes given by Papageorgiou (1979), based on data for 120 larger fish, mainly from l 0 to 18 cm. A linear regression of In W on lnL for our own data and Papageorgiou's data combined indicates a nonlinearity, and the two data sets used separately give somewhat different equations W(g) = 0.0063 L (cm) 3.o9

( 56 )

For Pagageorgiou's data (weighted): W(g) = 0.0037 L (cm) 3.40

( 57 )

For our data ( single fish ):

The two equations give approximately the same values in the length range 4 to 10 cm, but for fish in the length range 20 to 30 cm, Eqn. (57) gives 5080% higher weight than Eqn. (56). This difference, and the nonlinearity in the combined regression, may be explained by the fact that Papageorgiou's data is given as means for each age class, presumably as arithmetic means. For a skewed distribution, like the log-normal distribution, this would give higher values than the geometric means represented by the single fish loglinear regression in Eqn. (56). If we had used Papageorgiou's data or the combined set in our calculations, the biomass would probably be overestimated, since we have made explicit corrections for the skewed distributions in our formulae. Therefore, we have chosen to use the relation between length and weight based only on our own data. The estimates for coefficients with variance estimates in Eqn. (3) corresponding to Eqn. (56) are: c = 1.15, V(c) =0.00752, 2=2.012, d = 3.09, V ( d ) =0.082, and V(lnW(i,k) ) =0.0055. The estimated variance of In W for single fish includes variation due to difference between sexes, and thus can be applied when the sex distribution is unknown. The residuals show random variation as a function of fish length (DurbinWatson statistic = 2.06), and are near normally distributed except for a few extreme points.

69

HYDROACOUSTICAL ESTIMATES OF FISH SIZE DISTRIBUTION

Estimating biomass Table 7 shows results of calculations based on measurements made along selected transects in the central pelagic part of Lake Gjersjoen and Lake Arungen. Fish biomass is calculated as outlined above, using the estimated coefficients, as fish h a - ~ over the depth interval indicated in the table, varying from 6-8 m to 1-15 m, with 2 - 1 0 m as the most frequently used interval. The variation in measurement depth must be taken into account when results are compared. Table 7 shows that total fish biomass in the pelagic zone varies considerably with time for each of the two alternative T S - L relations, clearly demonTABLE 7 Fish biomass (F) calculated from echosounding recordings from the pelagic parts of Lake Gjersjoen and Lake Arungen. Alternatives for inL=f(TS) as in Table 6. P1. Estimated by data from Lindem and Sandlund (1984). P2. Estimated by maximum likelihood analysis of data from Lakes Gjersjoen and Arungen (Alt. 3 (I) ), based on multinominal distribution Lake

Date

Measured depth interval ( m - m )

Biomass over given depth interval (kg ha- t ) PI

P2

F Gjersjoen

May Sep. Jul. Aug. Nov.

Arungen

May Jul. Aug. Nov. Sep.

SE(F)

F

SE(F)

1980 1980 1981 1981 1981

2-10 1.5-9.5 2-10 2-10 2-10

175 195 2.4 2.1 27

37 37 0.7 0.6 5

380 420 3.2 2.7 50

230 240 1.3 1.1 25

1981 1981 1981 1981 t 1982 Oct. 1984

8-10 2-10 2-10 6-8 1-12 1-13

6.0 394 12 26-48 330 500

1.4 67 3 10 60 130

8.6 970 23 100-110 920 2300

3.9 530 13 60 540 1200

Statistical error factor (70% conf. intervals)

Main sources of stat. error (% of total variance ) 3a ,Jb

AD(i) 3.~(i) t Results of two different transects (cf. Table 8 ).

1.2-1.3

5-30 0-20 10-70 13-50

1.4-1.6

5-40 20-90 1-15 2-20

70

B. BJERKENGET AL.

strating the variation in pelagic behaviour of roach. This must be taken into account when deciding the echosounding sampling strategy, and seems to be of major importance for semi-pelagic species like roach. Judging from the floating gill net catches, roach during summer often inhabit water strata close to lake surface, thus avoiding echosounding recording. However, during spring (before spawning) and autumn, roach more regularly stay in the depth strata 2 m below the transducer and deeper. In these periods, rather reasonable fish biomass seems to be calculated in Gjersjoen. In/~rungen there is also a low value in May 1981, but in this case the measurements are made in a very restricted depth interval between 6 and 8 m. In November 1981 the echosounding indicates low pelagic biomass in both lakes. The P 1 estimate shows high levels of 200 kg h a - ~in Gjersjoen for the depth interval 2 - l 0 m, and 400-500 kg ha-~ in/~rungen, for depth intervals varying from 2-10 m to 1-13 m. This difference was to be expected, since Gjersjoen is deeper and less eutrophic than ,/~rungen. Alternative P2 gives values roughly twice as large as P l, but the two alternatives are consistent with each other when the error estimates are considered. The P2 alternative gives an exceptional high value of 2300 kg h a - ~for Arungenin October 1984. As seen from Table 8, this is due to a high density of fish giving echo signals in the range - 4 4 to - 3 8 dB, in the length range where the two alternative length estimates differ most. In Lake Arungen total roach biomass (including littoral and pelagic areas ) was estimated by the Petersen mark and recapture method as 150 to 700 kg h a - ~in the period 1979-1983 (Borgstrom, 1984). We therefore suggest that the lower estimates given by Alternative P 1 are the more reasonable. In any TABLE 8 Fish density (fish ha-~ ) in different TS classes, as calculated from Eqn. ( 1 ) for echosounding transects (depth intervals, Table 7 ) TS (dB)

-38 -40 -42 -44 -46 -48 -50 -52 -54 -56

Lake Gjersjoen

Lake/~rungen

May 1980

Sep. 1980

Jul. 1981

Aug. 1981

Nov. 1981

May 1981

Jul. 1981

Aug. 1981

Nov. 1981

Sep. 1982

Oct. 1984

0 0 147 389 2723 3999 770 3643 0 0

0 35 140 458 i!34 5062 4518 8836 0 0

0 0 0 0 0 0 38 277 687 225

0 0 0 0 0 0 l0 265 588 307

0 0 6 42 190 488 874 1386 1265 1167

0 0 0 0 0 ll6 143 496 914 596

0 92 896 1297 1903 3385 4500 19284 11919 16370

0 0 0 49 75 201 380 566 440 0

0 0 32 144 529 913 678 2129 257 0

26 26 743 1714 3083 3451 1977 3860 0 7

181 1419 392 693 1038 866 1743 503 3542 2677

29 0 0 47 179 170 309 743 0 259

HYDROACOUSTICALESTIMATESOF FISH SIZE DISTRIBUTION

71

case, the n u m b e r of echoes recorded will represent a m i n i m u m of the fish present along the recorded transect at the time of measurement. The P l estimates of fish biomass are relatively precise, with an estimated mean square error of 20-30%. The P2 estimates have twice this relative error. Owing to the log-normal assumptions, a large standard error will also shift the arithmetic means somewhat towards larger values (Eqns. 16 and 31 ), and part of the difference between the two estimates is thus due to the higher standard errors in the m a x i m u m likelihood estimates.

Importance o f different statistical error sources The relative contributions from the main sources of statistical error in the results for total fish biomass are also shown in Table 7. For Alternative P2, that is if we do not rely on Lindem and Sandlund's data for calculating length as a function of echo strength, the statistical errors in coefficients a and b in Eqn. (2) are the d o m i n a n t error sources, with the error in b as the most important. For Alternative P l, where the length-echo target strength relation is based on the data from L i n d e m and Sandlund (1984), the most important sources of error are the uncertainty of the fish density D (i) ( 10-70% of variance), and the TS scale resolution error ( 13-50% of variance), but the statistical errors of parameters a and b in the T S - L relation can also be important, and each can account for up to 20-30% of the variance. The other error sources, i.e. the statistical errors in the length-weight relation and the individual, stochastic variation in single fish target strength, are not important error sources, contributing at most 5% of the total variance. CONCLUSION

Our results indicate that it is possible to use echosounding to estimate fish biomass in lakes with an estimated statistical error given by a factor in the range 1.25 to 1.5, depending on what conclusions can be made about the correspondence between TS values and fish length. If the data of Lindem and Sandlund (1984) for the relation between TS and length, based on distribution peak identification, are representative of the real variation in the relation between TS distribution and true length distribution, then the lower estimates of fish biomass given as Alternative P 1 in Table 7, with a precision of 20-30% are valid. This precision is quite satisfactory, and the P 1 estimates of fish biomass compare well with other independent measurements. Within the given precision, our alternative estimates, based on the maxim u m likelihood method, confirm the coefficients of Lindem and Sandlund.

72

B. BJERKENGEl"AL.

However, the maximum likelihood calculations indicate larger variation in the correspondence between TS distributions and length distributions. If some of this variation is due to real variation in the relation between TS spectra and true length distributions, the practical error factor may be around 1.5, and the higher estimates P2 in Table 7 may be more representative. One possible source of error is that the calculation scheme for transforming TS spectra into density of fish size classes (Eqn. l ) might give large errors in individual density values Di. This is not taken completely into account in our estimates of the statistical errors in the fish density D (i), but should show in the error estimates for the TS-L relation with the maximum likelihood method. In addition, there might be real variations through the year in the correspondence between TS and size, especially for small fish due to growth of O + fish into the lowest detectable TS group through the season, or due to correlated behaviour. As mentioned previously, the peak identification method of estimating the TS-L relation is based on selected peaks, while the maximum likelihood method should be able to include such variations in the statistical error estimate. Thus, it is not obvious that the error factor is close to the lower limit. On the other hand, the higher estimate of the error factor includes any error due to selectivity of fish catches, and this indicates that this estimate might be too high. As a tentative conclusion we suggest that the error factor could be closer to 1.25 than to 1.5. ACKNOWLEDGEMENTS

Thanks are due to Torfinn Lindem (University of Oslo) for valuable discussions and criticism of the manuscript. He has kindly made raw data from Lakes Mjosa and Konnevesi available for our analysis. This study forms part of the research programme on eutrophication of inland waters supported by the Norwegian Council for Scientific and Industrial Research (NTNF), the Norwegian Council for Agricultural Research (NLVF), the Norwegian Institute for Water Research and the University of Oslo.

REFERENCES Afifi, A.A. and Azen, S.P., 1979. Statistical Analysis. A Computer Oriented Approach. Academic Press, London, second edition, 442 pp. Borgstr~)m, R., 1984. Investigations of fish, zooplankton and benthos in Lake ,~rungen. The Agricuitural Research Council of Norway / Section Soil Pollution Research, Tech. Rep., 538: l - 16. (in Norwegian ). Craig, R.E. and Forbes, S.T., 1969. Design of a sonar for fish counting. Fiskeridiv. Skr. Ser. Havunders., 15: 210-219.

HYDROACOUSTICAL ESTIMATES OF FISH SIZE DISTRIBUTION

73

Dickie, L.M., Dowd, K.G. and Boudreau, P.R., 1983. An echo counting and logging system (ECOLOG) for demersal fish size distributions and densities. Can. J. Fish. Aquat. Sci., 40: 487-498. Foote, K.G., 1981. Analysis of some errors associated with the use of target strength-to-length regressions in estimating fish abundance. In: J.B. Soumala (Editor), Meeting Hydroacoust. Methods Estimat. Mar. Fish Populat., Cambridge, MA, June 1979. Draper Laboratory, Cambridge, MA, 964 pp.. Foote, K.G., Lindem, T. and Brabrand, ~., 1987. Target strength of tiny roach. International Council for the Exploration of the Sea (ICES), ( 13:6 ): 1-4. Johannesson, K.A. and Losse, G.F., 1977. Methodology of acoustic estimations of fish abundance in some UNDP/FAO resource survey projects. Rapp. P.-V. Reun. Cons. Int. Explor. Mer., 170: 296-318. Lindem, T., 1981. The application of hydroacoustical methods in monitoring the spawning migration of whitefish (Coregonus lavaretus) in Lake Rands0orden, Norway. In: J.B. Soumala (Editor), Meeting Hydroacoust. Methods Estimat. Mar. Fish Populat., Cambridge, MA, 2529 June 1979. Draper Laboratory, Cambridge, MA, 964 pp. Lindem, T., 1982. Successes with conventional in situ determinations of fish target strength. International Council for the Exploration of the Sea (ICES): Symp. Fish. Acoust., Bergen, Norway, 21-24 June 1982, Contr. 53:1-18 (Mimeo.). Lindem, T. and Sandlund, O.T., 1984. New methods in assessment of pelagic freshwater fish stocks - - coordinated use of echosounder, pelagic trawl and pelagic nets. ("Ekkoloddregistrering av pelagiske fiskebestander i innsjoer." ) Fauna (Oslo), 37:105-111. (in Norwegian, abstract in English). Midtun, L. and Hoff, I., 1962. Measurements of the reflection of sound by fish. Fiskevidir. Skr. Ser. Havunders., 13 ( 13 ): 1-18. Nakken, O. and Olsen, K., 1977. Target strength measurements offish. Rap. P.-V. Reun. Cons. Int. Explor. Mer, 170: 52-69. Fevrier 1977. Nelder, J.A. and Mead, R., 1965. A simplex method for function minimization. Comput. J., 7: 308-313. Papageorgiou, N.K., 1979. The length weight relationship, age, growth and reproduction of the roach Rutilus rutilus (L.) in Lake Volvi. J. Fish Biol., 14: 529-538. Rudstam, L.G., Clay, C.S. and Magnusson, J.J., 1987. Density and size estimates of ciscoe ( Coregonus artedii) using analysis of echo peak PDF from single-transducers Sonar. Can. J. Fish. Aquat. Sci., 44:811-821. Solli, H., Brede, R. and Lindem, T., 1982. A new instrument for in situ measurements offish target strength distributions. International Council for the Exploration of the Sea (ICES): Symp. Fish. Acoust., Bergen, Norway, 21-24 June 1982, Contr. 54:1-15 (Mimeo.).