Estimating fin whale distribution from ambient noise spectra using Bayesian inversion

(1)

Estimating fin whale distribution from ambient noise spectra

using Bayesian inversion

1^ST JUNE, 2015

Author:

Sebastian Menze

University Supervisor:

Prof. Dr. Corinna Schrum First External Supervisor:

Dr. Hanne Sagen Second External Supervisor:

Prof. Dr. Halvor Hobæk

UNIVERSITY OF BERGEN, UNIVERSITY OFICELAND, UNIVERSITY OFAARHUS, UNIVERSITY OF THE FAROE ISLANDS

(2)

(3)

"Essentially, all models are wrong, but some are useful."

Box, G. E. P., and Draper, N. R. (1987) Empirical Model Building and Response Surfaces,John Wiley & Sons, New York, NY.

"Far better an approximate answer to the right question, which is often vague, than an exact answer to the wrong question, which can always be made precise."

John Tukey (1962) The future of data analysis. Annals of Mathematical Statistics 33 (1)

THIS THESIS WAS WRITTEN IN COOPERATION WITH THE NANSEN

ENVIRONMENTAL AND RÊMOTESÊNSING CÊNTER

(4)

(5)

Passive acoustic monitoring is increasingly used to study the distribution and migration of marine mammals. Marine mammal vocalizations are transient sounds, but the combined sound energy of a population continuously repeating a vocalization, adds up to a quasi-continuous chorus. Marine mammal choruses can be identified as peaks in ocean ambient noise spectra. In the North Atlantic, the fin whale chorus is commonly observed as peak at 20 Hz. This thesis proposes a method to estimate the distribution of vocalizing fin whales based on a set of fin whale chorus recordings. This is an extremely under-determined inverse problem. The method is based on Bayesian inverse theory and uses simulated annealing to estimate the most likely distribution of sound sources (vocalizing whales) on a geodesic grid. This includes calculating a transmission loss matrix connecting all grid nodes and recorders, using an arbitrary sound propagation model. Two models were successfully implemented: geometrical spreading and the ray trace model BELLHOP. The inversion method was tested under different scenarios. The results indicated that an imprecise transmission loss matrix is tolerated by the inversion method. The accuracy of the method depended mainly on the number and distribution of recorders. For the Norwegian sea, simulations showed that fin whale chorus inversion is possible using as few as 12 recorders between Iceland and Svalbard. An inversion based on data from published fin whale chorus observations indicated realistic winter distribution patterns. Existing methods to study marine mammal distribution are often confined to the summer months and a limited area. Future application of the proposed method admits automatic year-round monitoring of marine mammal distribution on a basin-wide scale.

(6)

1. Introduction

Recording and analyzing underwater sound is termed passive acoustic monitoring (PAM) in the scientific literature (Au and Hastings,2008). The term passive distinguishes the method from active acoustic monitoring, where sound signals are emitted actively (exam- ples are sonars or current measurement devices). Passive acoustic monitoring requires an array of sound receivers (hydrophones) that are either moored to the sea floor, towed behind a vessel or drifting freely. This thesis describes a new approach to analyze underwater sound recordings, based on inverse theory.

Underwater sound recordings contain signals from many different sources, ranging from breaking waves to marine mammals. The recorded sound signals can be divided into transient signals and a continuous acoustic background signal, the ambient noise.

In the ocean, most ambient noise stems from motion at the sea surface, such as wave breaking and oscillating bubble clouds (Carey and Evans,2011). Increasingly though, hu- man generated sound from marine traffic and seismic surveys dominates ambient noise between 25 and 200 Hz (Andrew et al.,2011,Frisk,2012). Vocalizing marine mammals are also a source of ambient noise. The combined sound of a population repeating the same vocalization many times, adds up to a constant marine mammal chorus, that can be identified as distinct peak in ambient noise spectra. These underwater choruses can be seen as oceanic analogue to the dawn chorus of birds and insects on land. The analysis of marine mammal choruses revealed extensive spatial as well as inter- and intra-annual variation (Burtenshaw et al.,2004,Nieukirk et al.,2012,Menze,2012).

The majority of PAM studies focus on the analysis of transient sounds from individual animals (Au and Hastings, 2008, Van Parijs et al., 2009, Mellinger et al., 2007). This thesis presents a different approach: Instead of determining the position of individual animals based on their transient vocalizations, the proposed method estimates the distribution of all vocalizing animals based on the chorus they produce. This is a highly under- determined inverse problem: the number of possible whale locations is much higher than the number of available recorders. In addition, the source levels and call rates vary between animals and seasons. To include these uncertainties into the distribution estimates, a numerical technique based on Bayesian inverse theory will be presented.

The method is developed for the chorus of North Atlantic fin whales (Balaenoptera physalus, Figure1), but can be applied to other marine mammal choruses as well. Fin whales are a highly migratory species, and little is known about their migration patterns and winter distribution (Víkingsson et al., 2013). Marine mammal distribution is tradi- tionally studied visually, using observers on ships and aircrafts (Marques,2009,Thomas et al.,2006). These surveys require considerable ship-time and personnel and are difficult to perform year-round (especially in darkness or sea ice). The inversion method described in this thesis can estimate the distribution of vocalizing fin whales year-round and on a basin scale. Thereby, the proposed method could close knowledge gaps in fin whale migration other methods cannot answer.

(9)

Figure 1: Fin whale surfacing near Svalbard in June 2012

The objectives of this thesis are summarized in the following research questions:

• Can we estimate the distribution of vocalizing animals from a set of chorus (their combined vocalizations) recordings?

• What are the capabilities and limitations of this method?

• How could it be implemented?

(10)

2. Marine mammals

The new method can be applied to all marine mammal species that produce a chorus.

Thus this section first introduces marine mammal distribution and acoustics in general and then focuses on fin whales. The term marine mammals includes cetaceans (toothed and baleen whales), seals, manatees and even polar bears (Perrin et al.,2009). Of these, only some baleen whale and seal species produce a marine mammal chorus. A chorus is only produced once a vocalization is repeated often enough, by enough animals and at a low enough frequency (due the reduced absorption of low frequency sound waves in seawater). A marine mammal chorus is a quasi-continuous signal, rendering it impossible to distinguish the vocalizations of individual animals. Table1summarizes the marine mammal species of which choruses have been observed and lists the respective references.

Table 1: Observed marine mammal choruses

Species Frequency

range Region References

Fin whales (Balaenoptera

physalus)

18 - 30 Hz North Atlantic

Simon et al.(2010), Nieukirk et al.(2012),

Klinck and Nieukirk (2012) Fin whales

(Balaenoptera physalus)

18 - 30 Hz and 85 -

105 Hz

Southern Ocean

Širovi´c et al.(2004), Gedamke et al.(2007),

Menze(2012),Dziak et al.(2015) Blue whales

(Balaenoptera musculus)

15 - 25 Hz North Pacific

Curtis et al.(1999), Burtenshaw et al.

(2004) Blue whales

(Balaenoptera musculus)

18 - 30 Hz and 26 - 27

Hz

Southern Ocean

Širovi´c et al.(2004), Gedamke et al.(2007),

Gavrilov et al.(2012), Menze(2012),Dziak

et al.(2015) Pygmy blue whales

(Balaenoptera musculus brevicauda)

65 - 75 Hz Southern Ocean Gedamke et al.(2007) Antarctic minke whales

(Balaenoptera bonaerensis)

100 - 300

Hz Southern Ocean Menze(2012)

Leopard seals (Hydrurga leptonyx)

320 - 350

Hz Southern Ocean Menze(2012)

(11)

In the ocean it is mainly baleen whales that produce a chorus. They are known for their extreme body sizes and belong to the top predators of most marine ecosystems (Perrin et al.,2009). Baleen whales have been exploited extensively in the last centuries, leaving several species close to extinction (Cressey,2015). To ensure sustainable management of marine mammal populations it is necessary to know each species population structure.

Because baleen whales are highly mobile, it has been difficult to distinguish between different stocks (Delarue et al.,2009,Castellote et al., 2004). Knowledge of migration patterns and the variability of baleen whale distribution can improve their management.

2.1. Distribution and migration

Baleen whales (Mysticeti) are known for their large scale migrations, that are commonly described as annual movement between high latitude feeding and low latitude breeding grounds (Perrin et al.,2009). The reasons for these large scale migrations are subject to speculation. Possible reasons range from the physiological effects of water temperature on calves to sea ice avoidance, inherited behavior or the avoidance of predation from killer whales (Corkeron and Connor,1999)

The exact distribution and migration patterns depend on each species’ ecological niche.

Baleen whales are most commonly observed in mid and high latitude waters, where they feed on aggregations of zooplankton or small fish. Baleen whale feeding is only efficient when the prey density is high enough (Goldbogen et al., 2011, 2007). In high latitude ecosystems, the spring or summer bloom of phytoplankton is followed by a respective increase in zooplankton biomass (Mann and Lazier,2013). Utilizing these seasonal peaks in prey availability is likely one of the major drivers of baleen whale migration and distribution. Their physiology allows them store large amount of energy as blubber, which also serves as insulation. This feature allows them to migrate over thousands of kilometers and renders the allocation of dense prey aggregations in their vast habitat a key element of their life history. For several baleen whales species it has been documented, that oceanic fronts or shelf slopes affect their distribution (Mauritzen et al.,2004,Ingram et al.,2007,Joiris et al.,2014). This is likely due to the increased productivity associated with frontal structures or upwelling along shelf or sea ice edges (Falk-Petersen et al., 2015).

Most studies on marine mammal distribution and behavior are based on visual surveys, either from the coast or using shipboard and aerial observers. Distance sampling using line transects provides statistically sound means to determine animal abundance in a pre- defined area (Marques,2009). However, these methods have a low temporal resolution since surveys are performed over a period of weeks or months. Due to the amount of ship-time and personnel required to observe the entire North Atlantic, these surveys are usually held every few years and provide little information on the migration patterns and intra-annual variation in marine mammal distribution.

With the rise of satellite technologies, tagging studies are increasingly used to inves- tigate migration routes. A recent study about right whales by Mate et al. (2015) documented the longest migration routes measured for any mammal so far (22511 km from Sakhalin Island to the Mexican coast). Tagging studies have proven useful to determine the behavior and migration patters of individual animals. They provide accurate informa-

(12)

2. Marine mammals

tion on an animals movement, but have the disadvantage of high cost, short sampling period (in the order of months) and very limited sample size. To study the migration of an entire population, it seem more feasible to combine visual and acoustic methods.

The next section will introduce how passive acoustic monitoring is used to study marine mammals.

2.2. Bioacoustics

Marine mammals utilize sound for foraging, orientation and communication. While toothed whales evolved high frequency vocalizations such as clicks and whistles, most baleen whales emit pulsed sounds between 10 to 1000 Hz (Au and Hastings,2008). Seals have been observed to emit a wide variety of sounds often described as grunts and moans between roughly 100 and 10000 Hz (Van Opzeeland,2010). Toothed whales utilize spe- cialized airs sacs and lips to produce clicks and whistles, whereas the sound production mechanisms of baleen whales remains largely unknown (Au and Hastings,2008).

Due to the low frequency of baleen whale vocalizations, their communication range is very large (hundreds of kilometers) compared to toothed whales. Whether baleen whales utilize these long transmission paths is seldom know. Due to their low frequency, it is un- likely that baleen whale use their vocalizations for echolocation. Humpback whales are known for their extensive acoustic repertoire, that are mainly thought of as male breeding display but may also serve other functions (Herman et al.,2013). Blue whales exhibit a far less diverse repertoire, that consists of breeding and feeding calls (Oleson and Calambokidis,2007). Blue and fin whales usually emit their calls in pulse trains, repeating a single call type over hours (Sirovi´c et al., 2007). Fin whale vocalizations will be discussed in the following section. Payne and Webb (1971) discuss the possibility that the monotonous calls of blue and fin whales act as long distance ranging signals, to lure possible mates to prey aggregations. Compared to the complex communication or echolocation found in other species, blue and fin whales likely use sound simply to find and attract mating partners in their vast habitat.

Baleen whale vocalizations are commonly recorded using moored or towed hydrophones.

Moored autonomous recorders are increasingly used to record the underwater sound- scape and marine mammal vocalizations over periods of years (Rettig et al.,2013). After recovering the moorings, the recordings are manually or automatically scanned for vocalizations (Towsey et al.,2013). From the detected vocalizations one can determine when calling animals are present in the detection range of the recorder. However, one cannot determine when the animals are absent, since they could simply not vocalize and still remain in the area.

When the recorder array consists of at least 3 hydrophones, that record the same vocalization simultaneously, the location of the animal can be determined from differences in the arrival time of the vocalization. This method is commonly termed triangulation (Clark, 1995,Au and Hastings,2008). Other methods of localizing vocalizing whales utilize the effects of multipath propagation (Valtierra et al.,2013, Newhall et al., 2012, McDonald and Fox,1999). Next to the use of hydrophones, whales have also been tracked using ocean bottom seismometers (Wilcock, 2012). However to determine marine mammal distribution it is necessary to know the approximate position of thousands of animals. To

(13)

track that many individuals on an ocean basin scale requires an unrealistically large number of recorders.

To connect the detection of vocalizations with animal density or distribution estimates, other approaches exist. These have to deal with the large uncertainties connected with animal behavior. For example, recording an increasing number of vocalizations does not automatically imply an increased abundance of whales, but could be caused by an increasing call rate or increasing detection probability (Helble et al.,2013b). Some of the most important variables to connect acoustic detection with animal abundance are the source pressure level (SL) of the animals, their call rate (CR) and the probability to detect the calls (Marques et al.,2013). Helble et al.(2013a) showed that ignoring the detection probability in passive acoustic surveys can introduce significant biases. The two most promising methods connecting acoustic detection and animal abundance, are distance sampling and mark-recapture methods (Marques et al.,2011, 2013,Küsel et al.,2011, Martin et al.,2013,Ward et al.,2012). However, they rely on recording individual whale calls. To determine the distribution of an entire population on an ocean basin-scale, this again requires an unrealistically large number of recorders. In this thesis, I present a different approach to estimate marine mammal distribution: Not by analyzing individual vocalizations but the chorus generated by all calling whales.

Depending on the local ambient noise levels from physical and anthropogenic sources, the marine mammal chorus can then be seen as a peak in recorded ambient noise spectra. Ambient noise spectra are calculated from spectral averages of sound recordings. To show how a spectrum changes over time, it is commonly visualized as a spectrogram. In a spectrogram, the time and frequency are represented by the x and y axis and the power spectral density by the color of each pixel. Figure2 shows a 3 year long ambient noise spectrogram from the Southern Ocean containing 3 whale choruses (Menze,2012). The horizontal line at 27 Hz is the blue whale chorus, the annually reoccurring line at 98 Hz the fin whale chorus and the annually reoccurring broad peak between 100 and 300 Hz is the Antarctic minke whale chorus. The spectrogram shows considerable intra- and inter- annual variation in the strength of the marine mammal choruses. This is connected to changes in the distribution and behavior of the vocalizing whales. How to use data from ambient noise spectra to estimate the distribution of calling whales will be presented in the methods section.

(14)

2. Marine mammals

Figure 2: Spectrogram of ambient noise recordings at 66^◦S, 0^◦E in the Southern Ocean.

X axis represents time, y axis frequency and color the power spectral density (PSD in dB re 1 µPa²Hz⁻¹). Sawtooth like pattern in the background represents noise variation caused by annual change in sea ice cover. Vertical spikes are caused by storm events and horizontal line represent the different marine mammal choruses as described in table1. FromMenze(2012).

2.3. Fin whales

Fin whales (Balaenoptera physalus) can be found from high to low latitudes, but are most common in cold and temperate waters. The global fin whale population can be divided into a Southern and Northern hemisphere stock, that rarely mix (Perrin et al., 2009).

During industrial whaling both populations were brought close to extinction. A recent study byRocha, Jr. et al.(2015) reports that about 874 068 fin whales were killed during industrialized whaling between 1909-1999. The Northern hemisphere stock can be separated into the North Atlantic and North Pacific stock. The management of whale stocks is an international political effort. Since the habitat of most whale species rarely falls into the jurisdiction of a single nation, international committees were formed. The international whaling commission (IWC) manages whales stocks on a global level. Several nations surrounding the North Atlantic, founded a separate commission to manage the marine mammal stocks in their jurisdiction: The North Atlantic Marine Mammal Commis- sion (NAMMCO). The IWC’s current estimate for the North Atlantic fin whale population is 26 500 animals (IWC,2015). The NAMMCO gives a similar number: Between 18,186 and 30,214 whales (Víkingsson et al.,2013). The IWC and NAMMCO divide the North Atlantic fin whale population into 5 stocks: West Greenland, East Greenland – Iceland, North Norway, West Norway and Faroe Islands and British Isles, Spain and Portugal (NAMMCO,2015,Donovan,1991). Figure3shows the boundaries of the putative stock

(15)

areas. The population structure adopted by the IWC and NAMMCO, is based on visual surveys and catch data from whaling operations (Donovan,1991). Given these sparse data, the stocks are only a rough approximation of the true population structure. A good example for this, is the fin whale population in the Mediterranean sea. They are thought to be separate from the North Atlantic stocks (Notarbartolo-di Sciara et al.,2003). However, an acoustic survey byCastellote et al. (2004), found that fin whales from the Icelandic stock seem to migrate through the Strait of Gibraltar into the Mediterranean. A genetic study on the other hand, roughly confirms the separate stocks as adopted by the IWC and NAMMCO (Bérubé et al.,1998).

Figure 3: Map showing the borders of the putative stock areas as used by the IWC and NAMMCO. FromDonovan(1991).

Since 1988, the NAMMCO conducts visual surveys of the fin whale distribution between Greenland and Europe (Víkingsson et al., 2013). These surveys found that fin whale sightings where clustered along the coast of Greenland, the Iceland-Faroe ridge as well as on the West coast of Svalbard (Figure4). The surveys were performed using several research vessels in June and July of each survey year, and indicate relatively sta- ble summer distribution patterns throughout the years.Víkingsson et al.(2013) estimated an increase in North Atlantic fin whale abundance from 17 482 animals in 1988 to 29 891 animals in 2001. In the Norwegian sea only, Øien(2004) estimated and abundance of 5,034 fin whales.

Other visual surveys on fin whale summer distribution report fin whales along the Span- ish coast and in the Bay of Biscay (Hammond et al.,2009), whereas comparatively few fin whales were sighted along the mid Atlantic ridge (Waring et al.,2008). All of the described surveys were performed in the summer months. Very little literature exists on

(16)

2. Marine mammals

sightings of fin whales during the three other seasons. A land based survey bySilva et al.

(2014) on the Azores, sighted fin whales mainly in spring and summer, but also in late autumn and winter. Fin whales that were tagged in spring off the Azores showed migratory behavior towards East Greenland and IcelandSilva et al.(2013).

Figure 4: Sightings of fin whales during the NASS ship surveys (1987 to 2001). Gray lines show the survey tracks and the circles fin whale sightings (symbol size is proportional to group size from 1 to 4+). FromVíkingsson et al.(2013)

The mating season of fin whales is between December and February (Perrin et al., 2009). Their winter distribution, presumably at their breeding grounds, remains largely unknown. Assuming a simple north-south migration between feeding and breeding grounds might oversimplify fin whale migration. It has been hypothesized that some populations do not migrate, or that one populations’ feeding ground is another populations’ breeding ground (Mizroch et al., 2009). Daylight and weather conditions render visual surveys difficult during winter, creating a lack of data on the distribution of fin whales during the

(17)

mating season. The calling activity of fin whales, which is mainly observed in winter (Watkins et al.,1987,Simon et al.,2010,Nieukirk et al.,2012), can help to fill this knowledge gap.

Fin whales emit pulsed vocalizations that have been recorded since the 1960s. Early underwater acoustic research was mainly led by the US navy, which was highly interested in this suspicious signal (Nichols,2005).Watkins et al.(1987) could eventually show that fin whales and not Russian submarines were the sources of these remarkably regular sounds. Pioneering work on underwater ambient noise led byWenz(1962) had already correctly identified whales as sources of ambient noise at certain frequencies.

Fin whale vocalizations are down-sweeps from 23 to 18 Hz, sometimes accompanied by a higher frequency down sweep (around 100 Hz) that varies between populations.

However, most energy is contained in the 20 Hz pulse. Figure5 shows an example fin whale pulse train and the fin whale chorus. The average fin whale source level has been estimated as 189 dB re 1 µPa (Weirathmueller et al.,2013,Watkins et al.,1987,Sirovi´c et al.,2007). The higher frequency component of the call can only be observed under low noise conditions, since it contains much less energy.Simon et al.(2010) have observed it around 130 Hz in West of Greenland. In the Southern Ocean it has been observed at 89 Hz west of the Antarctic peninsula (Sirovi´c et al.,2007) and at 98 Hz east of the Antarctic peninsula (Menze,2012).

Figure 5: Spectrogram of fin whale pulse train recored in Fram Strait: x axis represents time, y axis frequency and color the power spectral density. Fin whale calls visible as down-sweeps from 23 to 18 Hz. Fin whale chorus visible as red band at 20 Hz.

A study byCroll et al.(2002) in the Gulf of California found that only male fin whales vocalize. This supports the hypothesis that fin whale calls serve to attract mating partners, as proposed byPayne and Webb(1971). They furthermore observed a 1:1 sex ratio and that the whales foraged on dense patches of krill. Thus their conclusion was that male fin whales likely vocalize to attract females to favorable feeding locations (dense aggregation of prey).

The depth at which fin whales vocalize has been described as around 50 m byWatkins et al.(1987). The depth of the vocalizing animal influences the shape of the acoustic near field, due to interference with surface reflected sound waves.Weirathmueller et al.(2013)

(18)

2. Marine mammals

modeled this effect for fin whale calls and found a variability of±6 dB with depth in the near field. For the basin-scale distances (>100 km) between source and receiver considered in this thesis, near field interference due to source depth variation seems negligible.

The vocalizations are emitted in regular pulse trains, each lasting between hours and days (Watkins et al.,1987). The pulses are commonly emitted with a constant inter pulse interval (IPI). But it has also been observed that the pulse trains alternate between a shorter and a longer IPI (Watkins et al.,1987). In the North PacificWatkins et al.(1987) found IPIs ranging between 7-26 s with a mean of 12.48 s. Further studies showed that the IPI varies between populations, similar to the high frequency pulse component.

Delarue et al. (2009) found that the IPI recorded in the Gulf of Maine and Gulf of St.

Lawrence differs significantly and supports the notion of distinct populations. However, defining distinct populations based on their IPI remains difficult. Morano et al. (2012) found that the IPI of fin whales along the American coast changes seasonally: from 9.6 s in September–January to 15.1 s in March–May. A similar seasonal modulation of IPI exists in the North Pacific (Oleson et al.,2014).

Fin whale pulse trains have been frequently recorded in the North Atlantic, North Pacific and Southern Ocean (Nieukirk et al.,2012,Watkins et al.,2000,Širovi´c et al.,2004). The fin whale chorus (the combined pulse trains of all vocalizing animals) has been recorded and described in the North Atlantic and Southern Ocean. Figure6 shows a map of all published fin whale chorus observations.

Figure 6: Map showing fin whale chorus observations in the Northern (red) and Southern (blue) hemisphere.

The recorder array used byNieukirk et al.(2012) monitored the North Atlantic fin whale chorus with a high spatial resolution. The authors described the fin whales chorus with a measure termed fin index. The exact method will be discussed in section 5.1. Fig- ure7 displays the fin index time series recorded byNieukirk et al.(2012) using moored recorders along the mid Atlantic ridge. They found that the fin index is strongest north

(19)

of 32^◦N and increased during the late summer and fall months throughout the North At- lantic. Peak levels were reached in December and January and within the likely mating period.Simon et al.(2010) found that in the Davis strait (West of Greenland), the fin index peaked in November and December. The fin index in Davis strait was partly following a diel rhythm, being stronger during the dark hours and weaker during daylight hours. This indicates a link between calling activity and feeding. A similar diel marine mammal chorus rhythm can be found in Antarctic minke whales and is synchronous to the diel vertical migration of krill (Menze,2012). The fin index observed byKlinck and Nieukirk(2012) East of Greenland and in the Fram strait peaked in December and January. The northernmost fin whale chorus observation has been made byHope(2013) at 84^◦N, deep in the Arctic sea ice in September 2012.

Figure 7: Fin index observations in the mid North Atlantic byNieukirk et al.(2012)

(20)

3. Underwater sound propagation

3.1. Acoustic waves

Sound energy in a fluid propagates in the form of pressure waves. These can be described with the acoustic wave equation (Equation1,Lurton(2002)).

∇²p− 1 c²

∂²p

∂t² = 0 (1)

Here p represent the acoustic pressure,c the sound speed, x, y, z the spatial dimensions andt time. Solutions of varying complexity exist for equation 1. Two simple solutions are plane waves and spherical waves, which also form the basis of the transmission loss models described in section3.3. When assuming sinusoidal waves in one dimensional space with constant sound speed, plane waves (Eq. 2) are a solution to the acoustic wave equation.

p(x, t) =p₀e^(iω(t−^x^c⁾⁾ (2)

Hereωrepresents the circular frequency. This equation describes that surfaces of constant phase are planes normal to thexdirection.

If we on the other hand assume a point source in a homogeneous infinite three dimensional medium, spherical waves are a solution to equation1.

p(r, t) = p₀

r e^(iω(t−^r^c⁾⁾ (3)

These waves propagate simultaneously in all three dimensions. The surfaces of constant phase (wave front) take the form of spheres around the point source. The variable r represents the radius of this sphere. At ranges far away from the point source, the curvature of the wave front can be neglected. This allows us to approximately describe spherical waves using the plane wave equation.

The recording of underwater sound involves the use of hydrophones, that act as sound receivers. The sound pressure (emitted by a sound source) that arrives at a receiver is described by the sonar equation:

RL(x_r, y_r, z_r, t, f) = SL(x_s, y_s, z_s, t, f)−T L(x_s, y_s, z_s, x_r, y_r, z_r, t, f, τ)

−N L(x_r, yr, zr, t, f, τ) (4) In this equation,RLstands for received sound pressure level,SLfor the source sound pressure level,T L for the transmission loss and N Lfor the noise level at the receiver.

The variablesxs, ys, zs stand for the source position, those with index_r for receiver position,trepresents time,f frequency andτ varying environmental conditions (for example the sound speed profile). A wave propagating though the ocean is subject to attenuation, in the context of sound sources and receivers this is termed transmission loss. It is caused by absorption, scattering, refraction and reflection of the acoustic wave.

In sea water, absorption is an effect of the viscosity of water, chemical reactions and heat loss (Lurton,2002). To calculate the absorption loss depending on frequency, temperature, salinity and depth the empirical model ofFrancois and Garrison(1982) is widely

(21)

used. It describes how the viscosity of sea water, magnesium sulfate and boric acid ions absorb the energy of an acoustic wave. For the low frequencies considered in this thesis, the relaxation of boric acid ions and viscosity are the main causes for absorption. Figure8 shows the absorption coefficients for sea water over frequency. In general, the absorption coefficient increases with increasing frequency. In this thesis the Francois-Garrison model is part of all transmission loss estimates.

Figure 8: Absorption coefficients for sea water over frequency: for three temperatures values (0,10 and 20^◦C), a salinity of 35 PSU and a pH of 8. FromFrancois and Garrison(1982)

In addition to absorption, sound energy is attenuated by reflection and refraction. Re- flection occurs at the sea surface and the sea floor, but can also occur within the water column due to convex sound speed profiles.

The sound speedcin the ocean is not a constant, but a function of temperature, salinity and pressure. As with absorption, empirical models have been derived from measurements, leading to a standardized UNESCO reference model. In general sound speed increases mostly with temperature and pressure (Lurton,2002). Salinity has a minor but not negligible effect. In contrast to a medium with uniform sound speed, acoustic wave fronts in the ocean will thus not have the idealized spherical or plane shape. A sound wave is refracted when passing regions of variable sound speed. If the sound speed gradient is large enough this can even lead to the total internal reflection of a sound wave.

This process leads to the formation of so called sound channels, that are located at the

(22)

depth of sound speed minima (Lurton,2002). A wave entering such a channel will always be reflected back into it. This leads to reduced interaction with the sea floor and sea surface, and thus allows sound waves to propagate over large ranges and with little attenuation. In low and mid latitude oceans with warm surface layers a deep sound channel exists. It is often called the SOFAR (sound fixing and ranging) channel. In high latitude oceans, especially in ice covered regions, the sound speed minimum is located at or near the sea surface. This leads to a surface channel, where the wave fronts are repeatedly reflected towards the sea surface. Figure9illustrates the path of waves being trapped in a surface channel.

Figure 9: Propagation paths (rays) of plane waves given a sound speed profile with the sound speed minimum at the surface. The sound speed profile is depicted to the right. FromUrick(1983).

When the wave front hits the sea surface it is reflected according to Snells law (Lurton, 2002). However, a part of the sound energy is also scattered and absorbed by the rough sea surface or sea ice. How much the wave amplitude is reduced after each reflection varies with frequency and can be described with a reflection coefficient. Similarly to the absorption coefficient, higher frequencies are attenuated more than lower ones. Thus wave reflection at the sea surface or underside of sea ice, acts as a low pass filter (Dia- chok,1980). The magnitude of the transmission loss due to surface reflection changes with sea state and sea ice conditions. Reflection also occurs at the sea floor. However the reflection coefficient varies greatly with the material and structure of the sea floor.

Sediment that overlies the base rock can have similar sound speeds to water (Lurton, 2002). Depending on the structure of the sea floor, very complex wave propagation paths will occur. Implementing these correctly is one of the many challenges in underwater sound propagation modeling.

(23)

3.2. Geometrical spreading transmission loss

One of the simplest ways to calculate transmission loss is geometrical spreading. The attenuation due to sound speed variation as well as sea floor and surface reflection are neglected. It is assumed that initially, the sound wave propagates spherically according to equation 3. The intensity of the wave is defined as powerP passing through a unit area. Since energy is conserved, the intensityIof the wave will decrease with increasing radius of the sphere:

P = 4πr²Ir

P1m=Pr Ir

I1m = _4πr^4πr2² 1m

(5) If the absorption coefficientαis added and the intensity ratio is converted to the decibel scale, the equation for spherical transmission loss becomes:

T L(r) = 20log( r

r_1m) +αr (6)

This equation roughly holds until the sound wave hits the sea floor and sea surface.

From then on it can no longer propagate spherically but is trapped and forced to propagate in cylindrical form. The radius to where the wave propagates spherically is described by the radiusr_critical. To calculate the transmission loss for ranges larger than r_critical, one needs to assume cylindrical spreading. This results in the following equation where rcan be a range larger thanr_critical(Lurton,2002):

T L(r) = 20log(r_critical r1m

) + 10log( r rcritical

) +αr (7)

This simple description of the initial spherical and then cylindrical spreading of an acoustic wave in the ocean is a basic but still useful transmission loss estimate.

3.3. Ray trace sound propagation modeling

To include the effects of bathymetry and variable sound speeds on transmission loss, more sophisticated sound propagation models are required. A commonly used class of models for this purpose, are ray trace models (Lurton, 2002, Hovem, 2013). The key assumptions of raytracing is, that the wave front can be approximated as a fan of rays emerging from the source (Porter and Liu,1994). Figure9shows an example of multiple ray paths. Each of these rays represents the path of a plane wave (Eq.2). The path each ray takes is defined by Snells law of refraction and the sound speed gradient. If a ray hits the sea floor or sea surface, it is reflected according to Snells law and a specific reflection coefficient. The resulting sound pressure field can then be calculated by combining the sound pressure from each ray. Next to estimating the transmission loss, this method can also be used to calculate the impulse response for any given location (Porter and Liu, 1994).

Most available raytracing models are range dependent: They can calculate the sound pressure field depending on depth and radius from the sound source. Fully 3 dimensional raytracing models are theoretically possible, but not common due to their heavy computational requirements (Lurton,2002). Instead it is common to rotate the 2-D slices

(24)

created by a range depended model, over 360^◦ around the sources. This approach is often termed 2xN dimensional modeling. Unfortunately this methods does not implement ray reflections in the horizontal axis and is thus likely inaccurate around shelf slopes and sea mounts (Sturm et al.,2008). The rays interaction with the sea floor can be modeled with varying complexity, depending on the sea floor’s layering and properties. Often it is simply assumed as acousto-elastic half space with constant sound and shear wave speed.

Next to raytracing, other types of underwater sound propagation models exist. A large family of models is based on solving the wave equation using normal modes. These models can be used to simulate spreading of sound around islands and sea mounts, but rarely include range depend sound speed (Heaney et al.,2013,Ballard,2012). The most recent approaches include parabolic equation and fast field models, that have the advantage of being very accurate, but the disadvantage of heavy computational requirements (Lurton,2002). For the purpose of estimating the transmission loss between vocalizing fin whales and a set of recorders, ray trace models are considered sufficient enough (Mc- Donald and Fox,1999).

(25)

4. Inverse theory

4.1. Bayesian inversion

Inverse theory describes how to estimate a set of unknown parameters from indirect measurements (Tarantola,2006). It is almost always impossible, to obtain a perfect set of measurements describing all parameters of a system. We will, for example, never be able to record or sight all whales in the ocean. To estimate the whale distribution, an inverse theory which interprets the few recordings or sightings available is needed. The first step towards solving an inverse problem is, to know how the unknown parameters and observed data are connected. This part is termed the forward operator or forward model (Tarantola,2004). The forward operatorg describes which data vectordis produced by a parameter vectorm:

d = g(m) (8)

The forward operator g could for example stand for an acoustic propagation model, mfor a set of sound sources and d for sound pressure recordings. For many physical problems it is not possible to find the inverse ofg directly (Snieder and Trampert,1999).

This can be caused by the mathematical formulation of the forward problem or a lack of data d. In the case of marine mammal choruses we deal with sound pressure from many dispersed sources, so that many source combinations can produce similar sound pressure observations.

A special case occurs if the forward problem is linear, then g(m) is a function of the typed =Gm (Tarantola,2004). This simplifies the estimation of the correct solution to the inverse problem. However, even if the forward problem can be expressed as a set of linear equations, most inverse problems are ill-posed and ill-determined (Snieder and Trampert,1999). This means that there is not enough data pointsdas well as contradict- ing data, to determine a unique solutionm. Linear algebra determines that for a set of linear equations there are either infinitely many, one or no solutions at all. Similar dilem- mas can be found in nonlinear inverse theory.

Most progress in solving extremely under-determined inverse problems has been made in the field of geophysics (Snieder and Trampert, 1999). In this thesis I will follow the Bayesian approach described by Tarantola (2004). It embraces the fact, that the observed data and forward model are uncertain and can be described using probability distributions. Calculating with probability distributions has the advantage, that Bayesian reasoning can be applied (Tarantola,2004,Allmaras et al.,2013). This means that we can work with a priori (independent of observation, before reasoning) and a posteriori information (information acquired by observation, after reasoning).

The a priori information feeds into how the forward model is designed and describes the a priori probability distribution over the parameter space. The solution of the inverse problem is to find the posterior probability of the parametersmgiven the datad. To simplify the notation, I followTarantola(2004) and defineρ(m)as prior probability distribution over the parameter space andσ(m) as posterior probability distribution over the parameter space. Figure10shows a general scheme of an inverse problem including Bayesian reasoning. The steps to solve it include defining the a priori distributionρ(m), defining a

(26)

4. Inverse theory

suitable forward model to calculated_simulated fromm, appraising how well the simulated datad_simulated fit to the observed datad_observed and estimating the likely forward model parametersmgiven the observed data (the posterior distributionσ(m)).

Figure 10: General formulation of an inverse problem including Bayesian reasoning.

Adapted fromSnieder and Trampert(1999).

Every combination of forward model parametersm produces corresponding data as given by d = g(m). If we assume our forward model to be perfect, data and model will be connected by a line in the data and parameter space, as shown in figure 11.

Since knowing the perfect forward model is unrealistic, a better approach is to describe the output of the forward model with Gaussian uncertainties (with the standard deviation s). The forward model and data are then connected by the joint probability distribution θ(d|m), that describes the probability of the datadfor a given set of parametersm. It is shown in the right panel of figure11.

Figure 11: Comparison between a perfect forward model d = g(m) (left panel) and a forward model that produces the probability distributionθ(d|m)with the standard deviation s (right panel). X axis shows model parameter value, y axis data value. FromTarantola(2004).

(27)

The observations are given as probability distribution in data space: ρ(d). Normally the observed data has a very narrow probability distribution, so that it is often assumed perfect. Together the prior information in parameter and data space form the joint prior probability distributionρ(d,m). Figure12illustratesρ(d,m)as 2-D Gaussian distribution in the left panel. In this graphical example, model and data space have one dimension each. The marginal prior probability distribution in data and parameter space (ρ(d) and ρ(m)) are shown as one dimensional Gaussian to the left and below the panel.

Figure 12: General representation of a Bayesian solution of an inverse problem: x axis show model parameter value, y axis data value. The left panel displays the joint prior probability distribution ρ(d,m), the middle panel displays θ(d,m) (the probability the forward model assigns to the datad for a given parameter m), and the right panel displays the joint posterior probability distribution σ(d,m). FromTarantola(2004).

The joint posterior probability distributionσ(d,m)in data and parameter space can be calculated by combining the prior informationρ(d,m)with the information we have about the forward model θ(d|m). It is the conjunction of ρ(d,m) and θ(d|m), which can be seen in the right panel of figure12(Tarantola,2004). The solution of the inverse problem, the posterior probability distribution over the parameter spaceσ(m), is then a marginal probability distribution. It is shown below the right panel of figure12and can be calculated using the following equation (Tarantola,2004):

σ(m) =kρ(m) Z

D

ρ(d)θ(d|m)

µ(d) dd (9)

WhereDrepresents the data space,ka normalization constant andµ(d)the homogeneous probability density over the parameter space. If the parameter space is linear, the termR

D

ρ(d)θ(d|m)

µ(d) ddcan be replaced by the likelihood functionL(m) (Tarantola,2004).

It describes how well the parametersmfit to the datad, so that equation9becomes:

σ(m) =kρ(m)L(m) (10)

Herekis again a normalization constant. The likelihood functionL(m)is commonly expressed as exponential of the negative misfit functionS(m)(Mosegaard and Sambridge, 2002):

L(m) =ke^−S(m) (11)

(28)

4. Inverse theory

The uncertainty of the forward model (the standard deviations) needs to be considered when determining the misfit between observed and simulated data. The misfit is a measure of distance between the observed and simulated data points, also referred to as a norm in the literature (Tarantola,2004). The general equation for the misfitS(m)is:

S(m) = 1 p

ndata

X

i

|g_i(m)−di|^p

s^p_i (12)

Where p is the order of the norm ands the standard deviation of the forward model.

It is common to usel2 norm (p = 2) which is the sum of squared errors (SSE) between true and simulated data (Tarantola,2004). This has the advantage of being robust to the influence of outliers. However, it has the disadvantage of being sensitive to the location of the data points in data space, a property termed stability. Thel in the description of a norm stands for "least". Other measures of distance between data points change the order of the norm (betweenl₀ tol∞). In this thesis the misfit is determined using thel₂ norm:

S(m) = 1 2

ndata

X

i=1

(g_i(m)−d_i)²

s² (13)

To find the most likely parameters (for example the most likely distribution of vocalizing whales) that explain the observed datad, one needs to find the maximum of the posterior distributionσ(m), which equals finding the maximum of the likelihood function (Eq. 11):

max(σ(m)) =max(e^−S(m)) (14)

This approach is termed the maximum likelihood method (Allmaras et al.,2013). For well determined linear systems, analytical solutions to this problem exist. However for under-determined inverse problems (as most real world problems are), finding the maximum ofσ(m)is not a trivial task. Numerical techniques that sample the likelihood functions are needed. The next section presents an example problem with two forward model parameters to illustrate the misfit and likelihood functions.

(29)

4.2. Example problem with 2 parameters

This sections presents the misfit and likelihood functions for a simplified sound source and receiver problem. The example problem consists of 2 sound sources at fixed locations and a varying number of recorders. The emitted and received sound signals are assumed constant, so that time can be neglect. To simplify the problem, latitude, lon- gitude and depth of the sources and receiver are ignored and range is the only spatial dimension. The locations of the sound sources are at 50 and 80 km range and their source levels are 50 and 70 dB re 1 µPa respectively. Two recorder scenarios are used:

One where 10 recorders are distributed evenly over the range, and another where 100 recorders are also distributed evenly over the range. The received level the recorders measure is calculated using the sonar equation (Eq. 4) and the cylindrical spreading transmission loss equation (Eq.7). Figure13shows the locations as well as source and received levels for the two sources and recorder scenarios.

Figure 13: Example of parameter estimation problem including two fixed sources and cylindrical transmission loss over range. Blue dots represent the true source level for each source, lines with dots the received levels for a recorder at that location.

This example can now be reformulated as an inverse problem. The observed datad is given by the received levels of the recorders and the unknown parametersmare the source levels of the 2 sources. This means the likelihood and misfit functions will have 2 dimensions and are easy to visualize. The forward model is given by the sonar equation (Eq.4), that calculates the received level for each recorder, and the cylindrical spreading transmission loss equation (Eq. 7). This inverse problem is well-determined: There are more data pointsd(recorders) than parametersm(source level).

The misfit and likelihood values for each parameter combination between 0 and 100 dB re 1 µPa (source level combination) are calculated using equation 13 and 11. The forward model standard deviation is assumed as 5 dB. Figure14 shows the misfit and likelihood functions for both the 10 and 100 recorder scenario.

(30)

4. Inverse theory

Figure 14: Misfit and likelihood functions for the example problem: x and y axes show the SL of each source, Z axis and color the value of the misfit and likelihood function. Plots a and c show the likelihood and misfit for the 10 recorder scenario, plots b and d for the 100 recorder scenario. The white and black dots represent the minimum and maximum of the misfit and likelihood function. It is located at the true values of the parameters (true source levels).

The misfit functions show little optically discernible difference between the 10 and 100 recorder scenario. The likelihood functions however, show a more pronounced peak at the true source level values when 100 instead of 10 recorders are used. For both the 10 and 100 recorder scenario, the maximum of the likelihood function can be correctly identified. Using 10 instead of 100 recorders, adds regions of quasi-constant likelihood to the likelihood function. For this example problem, having a smaller number of observations dreduces the gradients of the likelihood functions. This will render it difficult to find the maximum of the likelihood functions with gradient based methods. The effect of forward model uncertainty on the likelihood function is shown in figure15. It can be seen that the peak of the likelihood function becomes broader with increasing forward model uncertainty. However, the maximum of the likelihood function stays identifiable up to 20 db forward model standard deviation.

Since the example problem had only two parameters, we could easily calculate the misfit and likelihood values for each parameter combination (100² combinations). This will become impossible for inverse problems with many parameters. Unfortunately most real world inverse problems include hundreds of parameters. Instead of calculating the

(31)

likelihood of all parameter combinations one can resort to finding the maximum of the likelihood function using numerical techniques. The next section describes these techniques.

Figure 15: Comparison of likelihood function using 100 recorders, thel2 norm and different forward model uncertainties: x and y axes show SL of each recorder, z axis and color the value of the likelihood function. Black points mark the true parameters and maximum of likelihood function.

4.3. Monte Carlo Markov chains

A useful method of finding the maximum of the posterior distributionσ(m), when the parameter space has many dimensions, are Monte Carlo Markov chains (Metropolis et al., 1953). The term Monte Carlo stands for the technique of describing a probability distribution with a set of discrete samples (Mosegaard and Sambridge,2002), whereas Markov chains describe a process where a transition from one state to another can be defined with a probability (Tarantola,2004). In the context of Bayesian inversion, this means that a set of forward model parameter samples are drawn from the prior probability distribution ρ(m). Each of these samples is a parameter vectorm(a solution of the inverse problem).

The samples are moved iteratively across parameter space to find the parameter values with maximum likelihood. For each iteration, the parameters are randomly altered. The change is accepted if it increases the likelihood of the solution, otherwise it is neglected and a new random change is produced. The consequent movement of the Monte Carlo Markov chain (MCMC) through parameter space, is commonly termed a random walk (Mosegaard and Sambridge,2002). When the likelihood function has only one maximum (similar to the 2 parameter example), the MCMC will find it after a sufficient number of

(32)

4. Inverse theory

iterations. Figure 16 displays an example random walk across the parameter space of the two parameter example problem.

Figure 16: Monte Carlo Markov chain performing a random walk through the parameter space of the example problem from section 4.2: x and y axis represent the two parameter space dimensions, z axis and color the value of the likelihood function. Black dots show the location of the Monte Carlo Markov chain at each iteration and arrows indicate the direction of its movement.

With the Monte Carlo Markov chain method it is possible to find the maximum of a function with hundreds of dimensions. However, many functions have more than one maximum. One can then distinguish between the global maximum and other local maxima (Mosegaard and Sambridge,2002). To find the maximum likelihood point, we need to identify the global maximum of the likelihood function. The key to finding the global maximum is to determine which moves of the MCMC are accepted and which are rejected.

This is not a trivial task: If one only accepts increases in likelihood (better solutions) and neglects all moves that lead to a likelihood decrease (worse solutions), the MCMC will get "stuck" at a local maximum.

A commonly used acceptance rule is given by the Metropolis-Hastings algorithm: it always allows increases in likelihood, and randomly accepts moves that decrease the likelihood (Metropolis et al., 1953). This results in a constant movement of the MCMC through parameter space and prohibits the MCMC from getting stuck at local maxima.

After a sufficient number of iterations the MCMC will converge towards the global maximum. However, when an inverse problem is extremely under-determined, the likelihood function can have large regions of almost equal likelihood as well as many local maxima.

This renders it difficult to find global maximum using the Metropolis-Hastings algorithm (Mosegaard and Vestergaard,1991). For these inverse problems, an algorithm termed

"simulated annealing" describes a useful acceptance rule.

The term annealing stems from metallurgy and describes the gradual cooling of a metal, which results in metals with small grain size and high ductility (Kirkpatrick et al., 1982). In the context of Monte Carlo Markov chains, simulated annealing is a modifi- cation of the Metropolis-Hastings algorithm (Kirkpatrick et al., 1982). In this case, the temperature is used to determine the acceptance of moves that decrease the likelihood.

(33)

Similar to metallurgic annealing, the temperature is reduced exponentially over the iterations (time). In the beginning of the algorithm’s run-time, the MCMC will move randomly through parameter space. However, over the iterations the acceptance rule is altered, and will eventually only accept moves that increase the likelihood. Analogous to the cooling of metal, the variable controlling the acceptance rule is termed the temperature. It has been shown byGranville et al.(1994) that this algorithm will converge to the global maximum within a finite amount of time. For the extremely under-determined inverse problems found in geophysics, parameter estimation with simulated annealing is widely used (Mosegaard and Vestergaard,1991).

(34)

5. Methods

This section will initially describe how to measure the fin whale chorus received level, and then explain the inversion method used to estimate fin whale distribution from these measurements.

5.1. Measuring chorus received level 5.1.1. Processing of spectral averages

The major part of a recorded underwater sound signal is noise (Lurton,2002,Carey and Evans,2011). This means, that the ambient noise spectrum can be calculated as power spectral density (PSD) of a recording (Duennebier et al.,2012). The power spectral density describes how the power of a signal is distributed over frequency. Basically, it can be expressed as the squared Fourier transformation of a time series, divided by the time series length (Smith,2003). Multiple methods exists to calculate the PSD, but Welch’s method (Welch,1967) is the most used (Smith,2003). It divides a time series into a sub- set of overlapping sections (Hamming windows) and then averages their squared Fourier transformations. All power spectra presented in this thesis have been calculated with this method.

In the literature, the strength of the fin whale chorus is often described with a measure termed the fin index (Nieukirk et al.,2012,Klinck and Nieukirk,2012,Simon et al.,2010).

Nieukirk et al.(2012) defines the fin index as the sum of the normalized PSD between f_min = 19.0Hz andf_max = 22.0Hz:

F in index(t) =

fmax

X

fmin

P SD(t, f)−P SD_median(t) (15) Here P SDmedian describes the median PSD between 0 and 55 Hz (Nieukirk et al., 2012). The sound pressure described by the fin index contains energy both from the fin whale chorus and other noise sources. This renders it unsuitable for an inversion.

Instead, a different method is proposed to calculate the received level of the fin whale chorus. It is based on Parsevals relation (Eq. 16,Smith(2003)). This relation states that the energy of a signal in the time domain (x(t)) equals its energy in the frequency domain (X(t)).

Z ∞

−∞

|x(t)|²dt= Z ∞

−∞

|X(f)|²df (16) Every underwater sound recording will contain energy from biotic, abiotic and anthropogenic sources. In the time domain we cannot differentiate between these noise sources. However, in the frequency domain we can utilize the fact that marine mammal choruses are usually confined to a narrow frequency band. From Parcevals relation, we can derive the recorded sound pressure for a given frequency band as:

precorded= s

Z fmax

fmin

F F Trecorded(f)df (17)

(35)

Here FFT stands for the Fourier transformation of the sound pressure signal. For the discrete values of the PSD, the integral can be written as sum over the desired frequency band:

precorded= v u u t

fmax−fmin

X

i=fmin

P SDrecorded(f) (18)

To now calculate the chorus sound pressure, we need to subtract the abiotic sound pressure from the recorded sound pressure. It can be interpolated from the PSD of frequency bands slightly below and above the chorus frequency band. In this thesis it was interpolated with an 8th degree polynomial from the recorded spectrum (excluding the fin whale chorus band between 15 and 25 Hz). Figure17shows an example ambient noise spectrum with distinct fin whale chorus peak. The interpolated abiotic PSD is displayed in red.

Figure 17: Spectrogram and ambient noise spectrum of a recording taken on the 3.10.2008 in Fram strait. Fin whale chorus is visible as a red band at 20 Hz in the spectrogram (upper panel). In the ambient noise spectrum, the fin whale chorus is visible as a peak at 20 Hz (lower panel). Here the black line represents the recorded spectrum and the red line the interpolated abiotic spectrum.

(36)

5. Methods

After determining the abiotic sound pressure, the chorus sound pressure and respective chorus received level can be calculated as follows:

p_chorus =p_recorded−p_abiotic (19)

RLchorus = 20log10(pchorus) (20)

How much the chorus differs from the abiotic noise is defined by the chorus signal-to- noise ratio (SNR):

SN R_chorus = 20log10(p_recorded)−20log10(p_abiotic) (21) The chorus signal-to-noise ratio gives similar values to the fin index introduced by Nieukirk et al.(2012). Both describe the intensity of the fin whale chorus relative to the abiotic noise. As an example, the next section presents the fin whale chorus recorded by a hydrophone array in Fram strait.

5.1.2. Example recordings from the DAMOCLES array

The DAMOCLES acoustic array was moored in Fram strait (Figure 18) between 2008 and 2009 (Sagen et al.,2010b,Mikhalevsky et al.,2015). DAMOCLES stands for: Devel- oping Arctic Modeling and Observing Capabilities for Long-term Environmental Studies.

The array consisted of a sound source that emitted 90 s sweeps between 200 and 300 Hz and a set of 8 hydrophones in a vertical string. The recording unit, the Simple Tomog- raphy Acoustic Receiver system (STAR), sampled with 1000 Hz. The recordings were kindly provided by the Nansen Environmental and Remote Sensing Center. The original purpose of the acoustic array was to monitor the mean ocean temperature of Fram strait (Sagen et al.,2008). In this thesis however, the recordings are used to calculate the fin whale chorus received levels.

Figure 18: Map showing the location of the DAMCLOES hydrophone array as a red dot.

The STAR unit recorded 100 s snippets every 3 h. However, due to technical problems many recordings were either missing, too short or contained electronic noise. These

(37)

recordings were sorted out on the basis of their low pass filtered amplitude (electronic noise had unrealistically high amplitude). In the remaining recordings, visual analysis revealed sparse occurrence of fin whale calls, but reoccurring presence of the fin whale chorus. The ambient noise PSD was calculated from the first 10 s of each recording using Welch’s method (Welch,1967) with a 1 s window length. This results in a frequency resolution of 1 Hz. The consequently obtained ambient noise spectrogram can be seen in figure19. The fin whale chorus can be seen as a horizontal band at 20 Hz between August and April. In June and July increased noise levels from seismic surveys dominate the spectrum between 10 and 100 Hz. This agrees with the analysis of recordings from Fram strait byMoore et al.(2011) andKlinck and Nieukirk(2012).

Figure 19: Longterm spectrogram of ambient noise recorded by the DAMOCLES tomog- raphy array, PSD of 10 s snippets using Welch’s method with 1 s window length, 50 % overlap. Fin whale chorus visible at 20 Hz.

The fin whale index, chorus received level and signal to noise ratio were calculated from the ambient noise spectra as described in the previous section. The results are displayed in figure20and indicate that the fin whale chorus is detectable from September to February in the Fram strait. Peak chorus received levels are reached in November. The few fin whale chorus received level observations larger than 100 dB re 1 µPa after De- cember 2008 are likely caused by individual vocalizing whales near the recorder. The fin index shows different patterns than the chorus received level. Upon closer examination, it becomes evident that some patterns in the fin index correspond to an increase in signal to noise ratio. Whilst other patterns correspond to an increase in chorus received level.

This demonstrates that the method used to describe a marine mammal chorus, significantly influences the results and interpretation of a study. Variation in signal to noise ratio might explain why the fin index varies strongly between 2 close recorders in the study of Simon et al.(2010).

For the inversion of marine mammal chorus observations, it is necessary to calculate the chorus received level and not the fin index. The fin index obtained from the DAMO- CLES array differs from observations byKlinck and Nieukirk(2012). They observed fin index values between 0 and 3. This could either be caused by inter-annual variations, or differences in the way the fin index was calculated.

(38)

5. Methods

Figure 20: Time series of fin whale chorus recorded by the DAMOCLES array.

(39)

5.2. Overview of inversion method

The aim of the inversion method is to find the most likely distribution of vocalizing fin whales, for a given set of chorus received level observations. Therefore we need to describe the unknown whale distribution with a set of parametersm. The datadis a vec- tor describing the received levels of the fin whale chorus recorded at different locations:

RL_observed. It could for example be (in dB re 1 µPa):

d=RL_observed = [93.45, 95.34, 91.45, ...94.28] (22) The inversion method can be divided into two processes: The forward model g(m) and the parameter estimation algorithm. Combined they form the inversion algorithm.

The flowchart in figure21shows the basic structure of the inversion algorithm. All input variables (received level data, a priori information and auxiliary variables) are displayed in yellow. The modules describing the forward model are displayed in blue, whereas the parameter estimation modules are shown in gray and white. The inversion algorithm uses Monte Carlo Markov chains and simulated annealing (introduced in section 4.3) to find the maximum of the posterior distribution (max(σ(m))). The output of the inversion algorithm is a set of parameters, that are solutions to the inverse problem. The most likely parameters (corresponding to the maximum of the posterior distribution) can be averaged from these solutions. A detailed discussion of the different modules will be given in the following sections.

All code developed for this thesis was written in Matlab 2014b. The source code can be found in the Appendix. Several functions are also available as binaries for Windows and Linux. Due to the high computational effort of the acoustic propagation model and inver- sions, most computations were run on a Unix cluster owned by the University of Bergen.

However all software is also able to run on the 2 cores of a normal consumer computer.

Future implementation of the inversion algorithm in C++ or Fortran will likely reduce the computational costs.

(40)

5. Methods

Figure 21: Flow chart showing the basic functionality of the inversion algorithm. Yellow

Estimating fin whale distribution from ambient noise spectra using Bayesian inversion